SliderCrank.py
You can view and download this file on Github: SliderCrank.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Slider crank test model
5#
6# Author: Johannes Gerstmayr
7# Date: 2019-11-01
8#
9# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
10#
11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12
13import exudyn as exu
14from exudyn.itemInterface import *
15from exudyn.utilities import *
16
17import numpy as np
18import matplotlib.pyplot as plt
19import matplotlib.ticker as ticker
20
21SC = exu.SystemContainer()
22mbs = SC.AddSystem()
23
24#++++++++++++++++++++++++++++++++
25#slider-crank
26# test nonlinear model; index2 and index3-formulation for ConnectorCoordinate and RevoluteJoint2D
27# crank is mounted at (0,0,0); crank length = 2*a0, connecting rod length = 2*a1
28
29#++++++++++++++++++++++++++++++++
30#ground object/node:
31
32background = GraphicsDataRectangle(-1, -2, 3, 2, color=[0.9,0.9,0.9,1.])
33
34oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background])))
35nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
36
37#++++++++++++++++++++++++++++++++
38#nodes and bodies
39g = 9.81 # gravity
40
41a0 = 0.25 #half x-dim of body
42b0 = 0.05 #half y-dim of body
43massRigid0 = 2
44inertiaRigid0 = massRigid0/12*(2*a0)**2
45graphics0 = GraphicsDataRectangle(-a0,-b0,a0,b0)
46
47a1 = 0.5 #half x-dim of body
48b1 = 0.05 #half y-dim of body
49massRigid1 = 4
50inertiaRigid1 = massRigid1/12*(2*a1)**2
51graphics1 = GraphicsDataRectangle(-a1,-b1,a1,b1)
52
53nRigid0 = mbs.AddNode(Rigid2D(referenceCoordinates=[a0,0,0],
54 initialVelocities=[0,0,0]));
55oRigid0 = mbs.AddObject(RigidBody2D(physicsMass=massRigid0,
56 physicsInertia=inertiaRigid0,
57 nodeNumber=nRigid0,
58 visualization=VObjectRigidBody2D(graphicsData= [graphics0])))
59
60nRigid1 = mbs.AddNode(Rigid2D(referenceCoordinates=[2*a0+a1,0,0], initialVelocities=[0,0,0]));
61oRigid1 = mbs.AddObject(RigidBody2D(physicsMass=massRigid1, physicsInertia=inertiaRigid1,nodeNumber=nRigid1,visualization=VObjectRigidBody2D(graphicsData= [graphics1])))
62
63c=0.05 #dimension of mass
64sliderMass = 1
65graphics2 = GraphicsDataRectangle(-c,-c,c,c)
66
67nMass = mbs.AddNode(Point2D(referenceCoordinates=[2*a0+2*a1,0]))
68oMass = mbs.AddObject(MassPoint2D(physicsMass=sliderMass, nodeNumber=nMass,visualization=VObjectRigidBody2D(graphicsData= [graphics2])))
69
70#++++++++++++++++++++++++++++++++
71#markers for joints:
72mR0Left = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRigid0, localPosition=[-a0,0.,0.])) #support point # MUST be a rigidBodyMarker, because a torque is applied
73mR0Right = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid0, localPosition=[ a0,0.,0.])) #end point; connection to connecting rod
74
75mR1Left = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid1, localPosition=[-a1,0.,0.])) #connection to crank
76mR1Right = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid1, localPosition=[ a1,0.,0.])) #end point; connection to slider
77
78mMass = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oMass, localPosition=[ 0.,0.,0.]))
79mG0 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=[0,0,0.]))
80
81#++++++++++++++++++++++++++++++++
82#joints:
83mbs.AddObject(RevoluteJoint2D(markerNumbers=[mG0,mR0Left]))
84mbs.AddObject(RevoluteJoint2D(markerNumbers=[mR0Right,mR1Left]))
85mbs.AddObject(RevoluteJoint2D(markerNumbers=[mR1Right,mMass]))
86
87#++++++++++++++++++++++++++++++++
88#markers for node constraints:
89mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
90mNodeSlider = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nMass, coordinate=1)) #y-coordinate is constrained
91
92#++++++++++++++++++++++++++++++++
93#coordinate constraints
94mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mNodeSlider]))
95
96#loads and driving forces:
97mbs.AddLoad(Torque(markerNumber = mR0Left, loadVector = [0, 0, 10])) #apply torque at crank
98
99#++++++++++++++++++++++++++++++++
100#assemble, adjust settings and start time integration
101mbs.Assemble()
102
103#now as system is assembled, nodes know their global coordinate index (for reading the coordinate out of the solution file):
104#deprecated: globalIndex = mbs.CallNodeFunction(nMass, 'GetGlobalODE2CoordinateIndex')
105globalIndex = mbs.GetNodeODE2Index(nMass)
106print('global ODE2 coordinate index of mass:', globalIndex)
107#alternatively: use mbs.systemData.GetObjectLTGODE2(oMass)[0] to obtain e.g. first coordinate index of sliding mass object
108
109simulationSettings = exu.SimulationSettings() #takes currently set values or default values
110
111simulationSettings.timeIntegration.numberOfSteps = 2*100000 #1000 steps for test suite/error
112simulationSettings.timeIntegration.endTime = 2 #1s for test suite / error
113#simulationSettings.timeIntegration.newton.relativeTolerance = 1e-10 #10000
114simulationSettings.timeIntegration.verboseMode = 1 #10000
115
116simulationSettings.solutionSettings.solutionWritePeriod = 1e-3
117
118simulationSettings.timeIntegration.newton.useModifiedNewton = True
119
120# simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
121# simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
122simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
123
124
125exu.StartRenderer()
126mbs.WaitForUserToContinue()
127
128#++++++++++++++++++++++++++++++++++++++++++
129#solve generalized alpha / index3:
130mbs.SolveDynamic(simulationSettings)
131
132SC.WaitForRenderEngineStopFlag()
133exu.StopRenderer() #safely close rendering window!
134
135
136u = mbs.GetNodeOutput(nMass, exu.OutputVariableType.Position) #tip node
137errorSliderCrankIndex3 = u[0] - 1.3513750614331235 #x-position of slider
138print('error errorSliderCrankIndex3=',errorSliderCrankIndex3)
139
140dataIndex3 = np.loadtxt('coordinatesSolution.txt', comments='#', delimiter=',')
141
142#++++++++++++++++++++++++++++++++++++++++++
143##solve index 2 / trapezoidal rule:
144#simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
145#simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
146#
147#mbs.SolveDynamic(simulationSettings)
148#
149#u = mbs.GetNodeOutput(nMass, exu.OutputVariableType.Position) #tip node
150#errorSliderCrankIndex2 = u[0] - 1.3528786319585837 #x-position of slider
151#print('error errorSliderCrankIndex2=',errorSliderCrankIndex2)
152#
153#dataIndex2 = np.loadtxt('coordinatesSolution.txt', comments='#', delimiter=',')
154#plt.plot(dataIndex2[:,0], dataIndex2[:,1+globalIndex], 'r-') #plot x-coordinate of slider
155
156plt.plot(dataIndex3[:,0], dataIndex3[:,1+globalIndex], 'b-') #plot x-coordinate of slider
157
158ax=plt.gca() # get current axes
159ax.grid(True, 'major', 'both')
160ax.xaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
161ax.yaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
162plt.tight_layout()
163plt.show()
164
165##animate solution
166#fileName = 'coordinatesSolution.txt'
167#solution = LoadSolutionFile('coordinatesSolution.txt')
168#AnimateSolution(mbs, solution, 10, 0.05)