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