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