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