sliderCrankFloatingTest.py
You can view and download this file on Github: sliderCrankFloatingTest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Slider crank model with verification in MATLAB for machine dynamics course
5# optionally, the slider crank is mounted on a floating frame, leading to vibrations
6# if the system is unbalanced
7# This example features 3D graphics of the links
8#
9# Author: Johannes Gerstmayr
10# Date: 2019-12-07
11#
12# 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.
13#
14#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
15
16import exudyn as exu
17from exudyn.utilities import *
18
19import numpy as np
20
21useGraphics = True #without test
22#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
23#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
24try: #only if called from test suite
25 from modelUnitTests import exudynTestGlobals #for globally storing test results
26 useGraphics = exudynTestGlobals.useGraphics
27except:
28 class ExudynTestGlobals:
29 pass
30 exudynTestGlobals = ExudynTestGlobals()
31#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32
33if useGraphics:
34 import matplotlib.pyplot as plt
35 import matplotlib.ticker as ticker
36
37SC = exu.SystemContainer()
38mbs = SC.AddSystem()
39
40#++++++++++++++++++++++++++++++++
41#ground object/node:
42
43#background = GraphicsDataRectangle(-0.5, -0.5, 1, 0.5, color=[1,1,1,1.]) #invisible background
44##background2 = GraphicsDataOrthoCube(-1, -1, -1, 2, -0.8, -0.8, color=[0.3,0.5,0.5,1.])
45#background2 = GraphicsDataCylinder(pAxis=[0,0.5,0],vAxis=[0,0,1],radius=0.3, color=[0.3,0.5,0.5,1.],
46# nTiles=16, angleRange=[0,pi*1.2], lastFace=True, cutPlain=True)
47#
48#background2 = GraphicsDataSphere(point=[0,0.5,0],radius=0.3,color=[0.3,0.5,0.5,1.],nTiles=8)
49#
50#background2 = GraphicsDataRigidLink(p0=[0,0.5,0],p1=[1,0.5,0], axis0=[0,0,1], axis1=[0,0,1],
51# radius=[0.1,0.1],thickness=0.2, width=[0.2,0.2],color=[0.3,0.5,0.5,1.],nTiles=16)
52
53solutionSliderCrankIndex2 = 0
54
55rangeTests = range(1,2) #(0,1): fixed frame, (1,2):floating frame
56if exudynTestGlobals.performTests: #consider shorter integration time
57 rangeTests = range(0,2)
58
59for testCases in rangeTests:
60
61 nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
62 mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
63
64
65 #++++++++++++++++++++++++++++++++
66 #floating body to mount slider-crank mechanism
67 constrainGroundBody = (testCases == 0) #use this flag to fix ground body
68
69 #graphics for floating frame:
70 #gFloating = GraphicsDataRectangle(-0.25, -0.25, 0.8, 0.25, color=[0.7,0.4,0.4,1.])
71 gFloating = GraphicsDataOrthoCube(-0.25, -0.25, -0.1, 0.8, 0.25, -0.05, color=[0.3,0.3,0.3,1.])
72
73 if constrainGroundBody:
74 floatingRB = mbs.AddObject(ObjectGround(referencePosition=[0,0,0], visualization=VObjectGround(graphicsData=[gFloating])))
75 mFloatingN = mbs.AddMarker(MarkerBodyPosition(bodyNumber = floatingRB, localPosition=[0,0,0]))
76 else:
77 nFloating = mbs.AddNode(Rigid2D(referenceCoordinates=[0,0,0], initialVelocities=[0,0,0]));
78 mFloatingN = mbs.AddMarker(MarkerNodePosition(nodeNumber=nFloating))
79 floatingRB = mbs.AddObject(RigidBody2D(physicsMass=2, physicsInertia=1, nodeNumber=nFloating, visualization=VObjectRigidBody2D(graphicsData=[gFloating])))
80 mRB0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nFloating, coordinate=0))
81 mRB1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nFloating, coordinate=1))
82 mRB2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nFloating, coordinate=2))
83
84 #add spring dampers for reference frame:
85 k=5000 #stiffness of floating body
86 d=k*0.01
87 mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mGround,mRB0], stiffness=k, damping=d))
88 mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mGround,mRB1], stiffness=k, damping=d))
89 mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mGround,mRB2], stiffness=k, damping=d))
90
91
92
93 #++++++++++++++++++++++++++++++++
94 #nodes and bodies
95 omega=2*pi/60*300 #3000 rpm
96 L1=0.1
97 L2=0.3
98 s1=L1*0.5
99 s2=L2*0.5
100 m1=0.2
101 m2=0.2
102 m3=0.4
103 M=0.1 #torque (default: 0.1)
104 #lambda=L1/L2
105 J1=(m1/12.)*L1**2 #inertia w.r.t. center of mass
106 J2=(m2/12.)*L2**2 #inertia w.r.t. center of mass
107
108 ty = 0.05 #thickness
109 tz = 0.05 #thickness
110 #graphics1 = GraphicsDataRectangle(-0.5*L1,-0.5*ty,0.5*L1,0.5*ty,color4steelblue)
111 #graphics1 = GraphicsDataOrthoCube(-0.5*L1,-0.5*ty,-tz,0.5*L1,0.5*ty,0,color4steelblue)
112 graphics1 = GraphicsDataRigidLink(p0=[-0.5*L1,0,-0.5*tz],p1=[0.5*L1,0,-0.5*tz],
113 axis0=[0,0,1], axis1=[0,0,1],radius=[0.5*ty,0.5*ty],
114 thickness=0.8*ty, width=[tz,tz], color=color4steelblue,nTiles=16)
115
116 #graphics2 = GraphicsDataRectangle(-0.5*L2,-0.5*ty,0.5*L2,0.5*ty,color4lightred)
117 #graphics2 = GraphicsDataOrthoCube(-0.5*L2,-0.5*ty,0,0.5*L2,0.5*ty,tz,color4lightred)
118 graphics2 = GraphicsDataRigidLink(p0=[-0.5*L2,0,0.5*tz],p1=[0.5*L2,0,0.5*tz],
119 axis0=[0,0,1], axis1=[0,0,1],radius=[0.5*ty,0.5*ty],
120 thickness=0.8*ty, width=[tz,tz], color=color4lightred,nTiles=16)
121
122 #crank:
123 nRigid1 = mbs.AddNode(Rigid2D(referenceCoordinates=[s1,0,0],
124 initialVelocities=[0,0,0]));
125 oRigid1 = mbs.AddObject(RigidBody2D(physicsMass=m1,
126 physicsInertia=J1,
127 nodeNumber=nRigid1,
128 visualization=VObjectRigidBody2D(graphicsData= [graphics1])))
129
130 #connecting rod:
131 nRigid2 = mbs.AddNode(Rigid2D(referenceCoordinates=[L1+s2,0,0],
132 initialVelocities=[0,0,0]));
133 oRigid2 = mbs.AddObject(RigidBody2D(physicsMass=m2,
134 physicsInertia=J2,
135 nodeNumber=nRigid2,
136 visualization=VObjectRigidBody2D(graphicsData= [graphics2])))
137
138
139 #++++++++++++++++++++++++++++++++
140 #slider:
141 c=0.025 #dimension of mass
142 graphics3 = GraphicsDataOrthoCube(-c,-c,-c*2,c,c,0,color4grey)
143
144 #nMass = mbs.AddNode(Point2D(referenceCoordinates=[L1+L2,0]))
145 #oMass = mbs.AddObject(MassPoint2D(physicsMass=m3, nodeNumber=nMass,visualization=VObjectMassPoint2D(graphicsData= [graphics3])))
146 nMass = mbs.AddNode(Rigid2D(referenceCoordinates=[L1+L2,0,0]))
147 oMass = mbs.AddObject(RigidBody2D(physicsMass=m3, physicsInertia=0.001*m3, nodeNumber=nMass,visualization=VObjectRigidBody2D(graphicsData= [graphics3])))
148
149 #++++++++++++++++++++++++++++++++
150 #markers for joints:
151 mR1Left = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRigid1, localPosition=[-s1,0.,0.])) #support point # MUST be a rigidBodyMarker, because a torque is applied
152 mR1Right = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid1, localPosition=[ s1,0.,0.])) #end point; connection to connecting rod
153
154 mR2Left = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid2, localPosition=[-s2,0.,0.])) #connection to crank
155 mR2Right = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid2, localPosition=[ s2,0.,0.])) #end point; connection to slider
156
157 mMass = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oMass, localPosition=[ 0.,0.,0.]))
158 mG0 = mFloatingN
159
160 #++++++++++++++++++++++++++++++++
161 #joints:
162 mbs.AddObject(RevoluteJoint2D(markerNumbers=[mG0,mR1Left]))
163 mbs.AddObject(RevoluteJoint2D(markerNumbers=[mR1Right,mR2Left]))
164 mbs.AddObject(RevoluteJoint2D(markerNumbers=[mR2Right,mMass]))
165
166 #++++++++++++++++++++++++++++++++
167 #markers for node constraints:
168 #mNodeSlider = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nMass, coordinate=1)) #y-coordinate is constrained
169 #coordinate constraints for slider (free motion in x-direction)
170 #mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mNodeSlider]))
171
172
173 #prismatic joint:
174 mRigidGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber = floatingRB, localPosition = [L1+L2,0,0]))
175 mRigidSlider = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oMass, localPosition = [0,0,0]))
176
177 mbs.AddObject(PrismaticJoint2D(markerNumbers=[mRigidGround,mRigidSlider], constrainRotation=True))
178
179
180 #user function for load; switch off load after 1 second
181 def userLoad(mbs, t, load):
182 if t <= 2: return load
183 return 0
184
185 #loads and driving forces:
186 mRigid1CoordinateTheta = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nRigid1, coordinate=2)) #angle coordinate is constrained
187 mbs.AddLoad(LoadCoordinate(markerNumber=mRigid1CoordinateTheta, load = M, loadUserFunction=userLoad)) #torque at crank
188 #mbs.AddLoad(Torque(markerNumber = mR1Left, loadVector = [0, 0, M])) #apply torque at crank
189
190 #++++++++++++++++++++++++++++++++
191 #assemble, adjust settings and start time integration
192 mbs.Assemble()
193
194 simulationSettings = exu.SimulationSettings() #takes currently set values or default values
195
196 simulationSettings.timeIntegration.numberOfSteps = 50000 #1000 steps for test suite/error
197 simulationSettings.timeIntegration.endTime = 3 #1s for test suite / error
198
199 if exudynTestGlobals.performTests: #consider shorter integration time
200 simulationSettings.timeIntegration.numberOfSteps = 5000 #1000 steps for test suite/error
201 simulationSettings.timeIntegration.endTime = 0.3 #1s for test suite / error
202
203 #simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8 #10000
204 simulationSettings.timeIntegration.verboseMode = 1 #10000
205
206 simulationSettings.solutionSettings.solutionWritePeriod = 2e-4
207 simulationSettings.timeIntegration.newton.useModifiedNewton = True
208 simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8
209 simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-8
210 simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
211
212 #++++++++++++++++++++++++++++++++++++++++++
213 #solve index 2 / trapezoidal rule:
214 simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
215 simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
216
217 dSize = 0.02
218 SC.visualizationSettings.nodes.defaultSize = dSize
219 SC.visualizationSettings.markers.defaultSize = dSize
220 SC.visualizationSettings.bodies.defaultSize = [dSize, dSize, dSize]
221 SC.visualizationSettings.connectors.defaultSize = dSize
222
223 #data obtained from SC.GetRenderState(); use np.round(d['modelRotation'],4)
224 SC.visualizationSettings.openGL.initialModelRotation = [[ 0.87758, 0.04786, -0.47703],
225 [ 0. , 0.995 , 0.09983],
226 [ 0.47943, -0.08761, 0.8732]]
227 SC.visualizationSettings.openGL.initialZoom = 0.47
228 SC.visualizationSettings.openGL.initialCenterPoint = [0.192, -0.0039,-0.075]
229 SC.visualizationSettings.openGL.initialMaxSceneSize = 0.4
230 SC.visualizationSettings.general.autoFitScene = False
231 #mbs.WaitForUserToContinue()
232
233 if useGraphics:
234 exu.StartRenderer()
235
236 mbs.SolveDynamic(simulationSettings)
237
238 if useGraphics:
239 #+++++++++++++++++++++++++++++++++++++
240 #animate solution
241# mbs.WaitForUserToContinue
242# fileName = 'coordinatesSolution.txt'
243# solution = LoadSolutionFile('coordinatesSolution.txt')
244# AnimateSolution(mbs, solution, 10, 0.025, True)
245 #+++++++++++++++++++++++++++++++++++++
246
247 SC.WaitForRenderEngineStopFlag()
248 exu.StopRenderer() #safely close rendering window!
249
250 u = mbs.GetNodeOutput(nMass, exu.OutputVariableType.Position) #tip node
251 exu.Print('sol =', abs(u[0]))
252 solutionSliderCrankIndex2 += abs(u[0]) #x-position of slider
253
254
255exu.Print('solutionSliderCrankIndex2=',solutionSliderCrankIndex2)
256exudynTestGlobals.testError = solutionSliderCrankIndex2 - 0.5916491633788333 #2020-01-15: 0.5916491633788333(corrected PrismaticJoint); 2019-12-26: 0.5916499441339551; 2019-12-15: 0.591689710999802 (absTol: 1e-8 now; 1e-2 before); before 2019-12-15: 0.5896009710727431
257exudynTestGlobals.testResult = solutionSliderCrankIndex2
258
259
260#plotResults = True#constrainGroundBody #comparison only works in case of fixed ground
261plotResults = useGraphics#constrainGroundBody #comparison only works in case of fixed ground
262if plotResults:
263 dataIndex2 = np.loadtxt('coordinatesSolution.txt', comments='#', delimiter=',')
264 #dataMatlab = np.loadtxt('slidercrankRefSolM0.1_tol1e-4.txt', comments='#', delimiter=',') #this is quite inaccurate
265 dataMatlab2 = np.loadtxt('slidercrankRefSolM0.1_tol1e-6.txt', comments='#', delimiter=',')
266
267 vODE2=mbs.systemData.GetODE2Coordinates()
268 nODE2=len(vODE2) #number of ODE2 coordinates
269
270 nAngle = mbs.systemData.GetObjectLTGODE2(oRigid1)[2] #get coordinate index of angle
271 plt.plot(dataIndex2[:,0], dataIndex2[:,1+nAngle], 'b-') #plot angle of crank;
272 plt.plot(dataIndex2[:,0], dataIndex2[:,1+nODE2+nAngle], 'r-') #plot angular velocity of crank
273 #plt.plot(dataMatlab[:,0], dataMatlab[:,2], 'g-') #plot angular velocity of crank from MATLAB
274 plt.plot(dataMatlab2[:,0], dataMatlab2[:,2], 'k-') #plot angular velocity of crank from MATLAB
275
276 #plt.plot(dataIndex3[:,0], dataIndex3[:,1+globalIndex], 'b-') #plot x-coordinate of slider
277
278 ax=plt.gca() # get current axes
279 ax.grid(True, 'major', 'both')
280 ax.xaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
281 ax.yaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
282 plt.tight_layout()
283 plt.show()