slidercrankWithMassSpring.py

You can view and download this file on Github: slidercrankWithMassSpring.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Slider crank with additional mass and spring (flexible slidercrank);
  5#           Example of paper Arnold, Brüls, 2007, Convergence of the generalized-[alpha] scheme for constrained mechanical systems, Multibody System Dynamics
  6#
  7# Author:   Johannes Gerstmayr
  8# Date:     2019-12-28
  9#
 10# 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.
 11#
 12#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 13
 14import exudyn as exu
 15from exudyn.itemInterface import *
 16from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 17import exudyn.graphics as graphics #only import if it does not conflict
 18
 19import numpy as np
 20import matplotlib.pyplot as plt
 21import matplotlib.ticker as ticker
 22
 23SC = exu.SystemContainer()
 24mbs = SC.AddSystem()
 25
 26
 27useGraphics = True
 28
 29nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 30mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 31
 32
 33#++++++++++++++++++++++++++++++++
 34#floating body to mount slider-crank mechanism
 35#constrainGroundBody = True #use this flag to fix ground body
 36
 37#graphics for floating frame:
 38gFloating = GraphicsDataRectangle(-0.25, -0.25, 0.8, 0.25, color=[0.95,0.95,0.95,1.])
 39#gFloating = graphics.BrickXYZ(-0.25, -0.25, -0.1, 0.8, 0.25, -0.05, color=[0.3,0.3,0.3,1.])
 40
 41oGround = mbs.AddObject(ObjectGround(referencePosition=[0,0,0], visualization=VObjectGround(graphicsData=[gFloating])))
 42mG0 = mbs.AddMarker(MarkerBodyPosition(bodyNumber = oGround, localPosition=[0,0,0]))
 43#mFloatingN = mbs.AddMarker(MarkerBodyPosition(bodyNumber = floatingRB, localPosition=[0,0,0]))
 44
 45
 46#++++++++++++++++++++++++++++++++
 47#nodes and bodies
 48omega=2*pi/60*300 #3000 rpm
 49L1=0.3
 50L2=0.6
 51L3=0.2
 52s1=L1*0.5
 53s2=L2*0.5
 54m1=0.36
 55m2=0.15
 56m3=0.1
 57m4=0.7
 58M=1 #torque (default: 0.1)
 59#lambda=L1/L2
 60J1=(m1/12.)*L1**2*1e-10 #inertia w.r.t. center of mass
 61J2=(m2/12.)*L2**2*1e-10 #inertia w.r.t. center of mass
 62
 63ty = 0.05    #thickness
 64tz = 0.05    #thickness
 65#graphics1 = GraphicsDataRectangle(-0.5*L1,-0.5*ty,0.5*L1,0.5*ty,graphics.color.steelblue)
 66#graphics1 = graphics.BrickXYZ(-0.5*L1,-0.5*ty,-tz,0.5*L1,0.5*ty,0,graphics.color.steelblue)
 67graphics1 = graphics.RigidLink(p0=[-0.5*L1,0,-0.5*tz],p1=[0.5*L1,0,-0.5*tz],
 68                                  axis0=[0,0,1], axis1=[0,0,1],radius=[0.5*ty,0.5*ty],
 69                                  thickness=0.8*ty, width=[tz,tz], color=graphics.color.steelblue,nTiles=16)
 70
 71#graphics2 = GraphicsDataRectangle(-0.5*L2,-0.5*ty,0.5*L2,0.5*ty,graphics.color.lightred)
 72#graphics2 = graphics.BrickXYZ(-0.5*L2,-0.5*ty,0,0.5*L2,0.5*ty,tz,graphics.color.lightred)
 73graphics2 = graphics.RigidLink(p0=[-0.5*L2,0,0.5*tz],p1=[0.5*L2,0,0.5*tz],
 74                                  axis0=[0,0,1], axis1=[0,0,1],radius=[0.5*ty,0.5*ty],
 75                                  thickness=0.8*ty, width=[tz,tz], color=graphics.color.lightred,nTiles=16)
 76
 77#crank:
 78nRigid1 = mbs.AddNode(Rigid2D(referenceCoordinates=[s1,0,0],
 79                              initialVelocities=[0,0,0]));
 80oRigid1 = mbs.AddObject(RigidBody2D(physicsMass=m1,
 81                                    physicsInertia=J1,
 82                                    nodeNumber=nRigid1,
 83                                    visualization=VObjectRigidBody2D(graphicsData= [graphics1])))
 84
 85#connecting rod:
 86nRigid2 = mbs.AddNode(Rigid2D(referenceCoordinates=[L1+s2,0,0],
 87                              initialVelocities=[0,0,0]));
 88oRigid2 = mbs.AddObject(RigidBody2D(physicsMass=m2,
 89                                    physicsInertia=J2,
 90                                    nodeNumber=nRigid2,
 91                                    visualization=VObjectRigidBody2D(graphicsData= [graphics2])))
 92
 93
 94#++++++++++++++++++++++++++++++++
 95#slider:
 96c=0.025 #dimension of mass
 97graphics3 = graphics.BrickXYZ(-c,-c,-c*2,c,c,0,graphics.color.grey)
 98
 99nMass3 = mbs.AddNode(Point2D(referenceCoordinates=[L1+L2,0]))
100oMass3 = mbs.AddObject(MassPoint2D(physicsMass=m3, nodeNumber=nMass3,visualization=VObjectRigidBody2D(graphicsData= [graphics3])))
101
102nMass4 = mbs.AddNode(Point2D(referenceCoordinates=[L1+L2+L3,0]))
103oMass4 = mbs.AddObject(MassPoint2D(physicsMass=m4, nodeNumber=nMass4,visualization=VObjectRigidBody2D(graphicsData= [graphics3])))
104
105#++++++++++++++++++++++++++++++++
106#markers for joints:
107mR1Left = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid1, localPosition=[-s1,0.,0.])) #support point # MUST be a rigidBodyMarker, because a torque is applied
108mR1Right = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid1, localPosition=[ s1,0.,0.])) #end point; connection to connecting rod
109
110mR2Left = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid2, localPosition=[-s2,0.,0.])) #connection to crank
111mR2Right = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid2, localPosition=[ s2,0.,0.])) #end point; connection to slider
112
113mMass3 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oMass3, localPosition=[ 0.,0.,0.]))
114
115mMass4 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oMass4, localPosition=[ 0.,0.,0.]))
116
117
118#++++++++++++++++++++++++++++++++
119#joints:
120mbs.AddObject(RevoluteJoint2D(markerNumbers=[mG0,mR1Left]))
121mbs.AddObject(RevoluteJoint2D(markerNumbers=[mR1Right,mR2Left]))
122mbs.AddObject(RevoluteJoint2D(markerNumbers=[mR2Right,mMass3]))
123
124#++++++++++++++++++++++++++++++++
125#markers for node constraints:
126mNodeSliderX = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nMass3, coordinate=0)) #y-coordinate is constrained
127mNodeSliderY = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nMass3, coordinate=1)) #y-coordinate is constrained
128mNodeSliderX2= mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nMass4, coordinate=0)) #y-coordinate is constrained
129mNodeSliderY2= mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nMass4, coordinate=1)) #y-coordinate is constrained
130#coordinate constraints for slider (free motion in x-direction)
131mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mNodeSliderY]))
132mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mNodeSliderY2]))
133
134#add spring between mass 3 and 4
135mbs.AddObject(ObjectConnectorCoordinateSpringDamper(markerNumbers = [mNodeSliderX, mNodeSliderX2],
136                                                    stiffness = 1000))
137
138#+++++++++++++++++++++++++++++++++++++++++
139#loads and driving forces:
140
141mRigid1CoordinateTheta = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nRigid1, coordinate=2)) #angle coordinate is constrained
142constraintCrankAngle = mbs.AddObject(CoordinateConstraint(markerNumbers=[mRigid1CoordinateTheta, mGround], offset = -1.*np.pi/2.))
143
144mbs.AddLoad(LoadCoordinate(markerNumber=mRigid1CoordinateTheta, load = M)) #torque at crank
145
146
147
148#++++++++++++++++++++++++++++++++
149#assemble, adjust settings and start time integration
150mbs.Assemble()
151if useGraphics:
152    exu.StartRenderer()
153    #mbs.WaitForUserToContinue()
154
155simulationSettings = exu.SimulationSettings() #takes currently set values or default values
156
157initCrank = True
158if initCrank:
159    #turn crank to 90° as enforced by constraintCrankAngle
160    mbs.SolveStatic(simulationSettings)
161
162    #use static solution as initial conditions for dynamic solution
163    currentState = mbs.systemData.GetSystemState()
164    mbs.systemData.SetSystemState(currentState, configuration=exu.ConfigurationType.Initial)
165
166    mbs.SetObjectParameter(constraintCrankAngle, 'activeConnector', False)
167    #mbs.WaitForUserToContinue()
168
169h = 5e-3   #5e-3 in paper of Arnold and Bruls
170T = 1
171simulationSettings.timeIntegration.endTime = T               #1s for test suite / error
172simulationSettings.timeIntegration.numberOfSteps = int(T/h)  #1000 steps for test suite/error
173
174#simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8 #10000
175simulationSettings.timeIntegration.verboseMode = 1 #10000
176
177simulationSettings.solutionSettings.solutionWritePeriod = 1e-3
178#simulationSettings.timeIntegration.newton.useModifiedNewton = False
179simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.6 #0.7 in paper of Arnold and Bruls
180
181#++++++++++++++++++++++++++++++++++++++++++
182#solve index 2 / trapezoidal rule:
183simulationSettings.timeIntegration.generalizedAlpha.useNewmark = False
184simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = False
185
186dSize = 0.02
187SC.visualizationSettings.nodes.defaultSize = dSize
188SC.visualizationSettings.markers.defaultSize = dSize
189SC.visualizationSettings.bodies.defaultSize = [dSize]*3
190SC.visualizationSettings.connectors.defaultSize = dSize
191
192mbs.SolveDynamic(simulationSettings)
193
194if useGraphics:
195    #+++++++++++++++++++++++++++++++++++++
196    #animate solution
197#        mbs.WaitForUserToContinue
198#        fileName = 'coordinatesSolution.txt'
199#        solution = LoadSolutionFile('coordinatesSolution.txt')
200#        AnimateSolution(mbs, solution, 10, 0.025, True)
201    #+++++++++++++++++++++++++++++++++++++
202
203    #SC.WaitForRenderEngineStopFlag()
204    exu.StopRenderer() #safely close rendering window!
205
206u = mbs.GetNodeOutput(nMass4, exu.OutputVariableType.Position) #tip node
207print('sol =', abs(u[0]))
208solutionSliderCrank = abs(u[0]) #x-position of slider
209
210
211print('solutionSliderCrankIndex2=',solutionSliderCrank)
212
213
214plotResults = useGraphics#constrainGroundBody #comparison only works in case of fixed ground
215if plotResults:
216    data = np.loadtxt('coordinatesSolution.txt', comments='#', delimiter=',')
217
218    vODE2=mbs.systemData.GetODE2Coordinates()
219    nODE2=len(vODE2) #number of ODE2 coordinates
220
221    nAngle = mbs.systemData.GetObjectLTGODE2(oRigid1)[2] #get coordinate index of angle
222    nM3 = mbs.systemData.GetObjectLTGODE2(oMass3)[0] #get X-coordinate of mass 4
223    nM4 = mbs.systemData.GetObjectLTGODE2(oMass4)[0] #get X-coordinate of mass 4
224    print("nAngle=", nAngle)
225    print("nM3=", nM3)
226    print("nM4=", nM4)
227
228    plt.plot(data[:,0], data[:,1+nAngle], 'b-') #plot angle of crank;
229    #plt.plot(data[:,0], data[:,1+nM3], 'g-')    #Y position of mass 3
230    plt.plot(data[:,0], data[:,1+nM4], 'r-')    #Y position of mass 4
231
232    ax=plt.gca() # get current axes
233    ax.grid(True, 'major', 'both')
234    ax.xaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
235    ax.yaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
236    plt.tight_layout()
237    plt.show()