driveTrainTest.py
You can view and download this file on Github: driveTrainTest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test of Node1D, Mass1D and Rotor1D for drivetrains;4-piston compressor
5# Uses different constraints to realize the drivetrain;
6# includes joints to connect 3D bodies and 1D drivetrain;
7# the crank is elastically supported (except z-direction);
8# this makes a direct coupling of the rotation angle to the drivetrain more difficult
9#
10# Author: Johannes Gerstmayr
11# Date: 2020-04-22
12#
13# 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.
14#
15#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
16
17import exudyn as exu
18from exudyn.utilities import *
19from exudyn.FEM import *
20
21import numpy as np
22
23useGraphics = True #without test
24#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
25#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
26try: #only if called from test suite
27 from modelUnitTests import exudynTestGlobals #for globally storing test results
28 useGraphics = exudynTestGlobals.useGraphics
29except:
30 class ExudynTestGlobals:
31 pass
32 exudynTestGlobals = ExudynTestGlobals()
33#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
34# useGraphics = False
35
36SC = exu.SystemContainer()
37mbs = SC.AddSystem()
38
39color = [0.1,0.1,0.8,1]
40s = 0.1 #width of cube
41sx = 3*s #length of cube/body
42
43background0 = GraphicsDataRectangle(-1,-1,1,1,color4white)
44oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
45 visualization=VObjectGround(graphicsData= [background0])))
46
47L0 = 0.2 #crank length
48L1 = 0.5 #connecting rod length
49L2 = 0.1 #cylinder length
50a = 0.025 #general width dimension
51c = 0.05 #piston radius
52
53discBig = 0.15
54discSmall = 0.05
55
56rho=7850 #steel density
57
58kJoint = 1e4 #for elastic support of crankshaft
59dJoint = 2e2 #for elastic support of crankshaft
60
61#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
62#add crank:
63inertiaCrank = InertiaCuboid(density=rho, sideLengths=[L0,a,a])
64
65
66omega0 = [0,0,2*pi*100*0] #initial angular velocity of bodies
67#ep0 = eulerParameters0
68p0 = [0,0,0] #reference position / COM of crank
69v0 = [0,0,0] #initial translational velocity
70
71color0= color4grey
72gGraphics0 = GraphicsDataOrthoCube(-0.5*L0,-a*0.9,-a*0.9,0.5*L0,a*0.9,a*0.9, color0)
73gGraphics0b = GraphicsDataOrthoCube(-0.5*L0,-a*0.9,-a*0.9+4*a,0.5*L0,a*0.9,a*0.9+4*a, color0)
74gGraphics0c = GraphicsDataCylinder([0.5*L0,0,-a],[0,0,6*a],a, color4darkgrey)
75gGraphics0d = GraphicsDataCylinder([-0.5*L0,0,3*a],[0,0,4*a],a, color4darkgrey)
76
77oRB0 = mbs.CreateRigidBody(inertia=inertiaCrank,
78 referencePosition=p0,
79 referenceRotationMatrix=np.eye(3),
80 initialVelocity=v0,
81 initialAngularVelocity=omega0,
82 gravity=[0.,-9.81*0,0.],
83 graphicsDataList=[gGraphics0,gGraphics0b,gGraphics0c,gGraphics0d])
84
85nRB0 = mbs.GetObject(oRB0)['nodeNumber']
86
87mGround0 = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oGround, localPosition = [0,0,-a]))
88mRigid0 = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oRB0, localPosition = [0,0,-a]))
89mRigid0RotationCoordinate = mbs.AddMarker(MarkerNodeRotationCoordinate(nodeNumber=nRB0, rotationCoordinate=2))#z-axis
90
91useElasticSupport=True
92if not useElasticSupport:
93 mbs.AddObject(GenericJoint(markerNumbers = [mRigid0,mGround0],
94 constrainedAxes=[1,1,1, 1,1,0],
95 visualization= VObjectJointGeneric(axesRadius=a,axesLength=4*a)))
96else:
97 mGround0b = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oGround, localPosition = [0,0,6*a]))
98 mRigid0b = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oRB0, localPosition = [0,0,6*a]))
99 mbs.AddObject(GenericJoint(markerNumbers = [mRigid0,mGround0],
100 constrainedAxes=[0,0,1, 1,1,0],
101 visualization= VObjectJointGeneric(axesRadius=a,axesLength=4*a)))
102 mbs.AddObject(CartesianSpringDamper(markerNumbers = [mRigid0,mGround0],
103 stiffness=[kJoint,kJoint,0],
104 damping = [dJoint,dJoint,0]))
105 mbs.AddObject(CartesianSpringDamper(markerNumbers = [mRigid0b,mGround0b],
106 stiffness=[kJoint,kJoint,0],
107 damping = [dJoint,dJoint,0]))
108
109sCrankPos = mbs.AddSensor(SensorNode(nodeNumber=nRB0, #fileName="solution/sensorCrankPos.txt",
110 storeInternal=True,
111 outputVariableType=exu.OutputVariableType.Position))
112sCrankAngVel = mbs.AddSensor(SensorNode(nodeNumber=nRB0, #fileName="solution/sensorCrankAngVel.txt",
113 storeInternal=True,
114 outputVariableType=exu.OutputVariableType.AngularVelocity))
115sCrankAngle = mbs.AddSensor(SensorNode(nodeNumber=nRB0, #fileName="solution/sensorCrankAngle.txt",
116 storeInternal=True,
117 outputVariableType=exu.OutputVariableType.Rotation))
118
119#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
120#add connecting rods and pistons:
121for i in range(4):
122 inertiaConrod = InertiaCuboid(density=rho, sideLengths=[L1,a,a])
123 phi = i/4*(2*pi) #rotation of subsystem
124 A = RotationMatrixZ(phi) #transformation of subsystem
125
126
127 A2 = RotationMatrixZ(0)
128 v2 = np.array([0,0,0])
129 offsetPiston=0
130 if i==1 or i==3:
131 phi2=np.arctan2(0.5*L0,L1) #additional conrod angle
132 A2 = RotationMatrixZ(phi2) #additional transformation of conrod
133 v2 = A@np.array([-0.5*L0,-0.25*L0,0])
134 offsetPiston = L1*(1.-np.cos(phi2)) + 0.5*L0
135
136 offZ = 0
137 if i < 2:
138 Acrank = RotationMatrixZ(0)
139 else:
140 Acrank = RotationMatrixZ(pi)
141 offZ = 4*a
142
143 color1= color4grey
144 gGraphics1 = GraphicsDataOrthoCube(-0.5*L1,-a*0.9,-a*0.9,0.5*L1,a*0.9,a*0.9, color1)
145 omega1 = [0,0,0] #initial angular velocity of bodies
146
147 p1 = [0.5*L0+0.5*L1,0,2*a+offZ] #reference position / COM of crank
148 p1 = list(v2 + A @ np.array(p1))
149
150 v1 = [0,0,0] #initial translational velocity
151
152 oRB1 = mbs.CreateRigidBody(inertia=inertiaConrod,
153 referencePosition=p1,
154 referenceRotationMatrix=A@A2,
155 initialAngularVelocity=omega1,
156 initialVelocity=v1,
157 gravity=[0.,-9.81*0,0.],
158 graphicsDataList=[gGraphics1])
159
160
161 locPos0 = list(Acrank @ np.array([0.5*L0,0,a+offZ]))
162 mbs.CreateGenericJoint(bodyNumbers=[oRB0,oRB1],
163 position=locPos0, constrainedAxes=[1,1,1, 1,1,0],
164 useGlobalFrame=False,
165 axesRadius=0.5*a,axesLength=2*a)
166 #alternative, using markers and objects:
167 # mRigid01 = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oRB0, localPosition = locPos0)) #connection crank->conrod
168 # mRigid10 = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oRB1, localPosition = [-0.5*L1,0,-a]))#connection conrod->crank
169 # mbs.AddObject(GenericJoint(markerNumbers = [mRigid01,mRigid10],
170 # constrainedAxes=[1,1,1, 1,1,0],
171 # visualization= VObjectJointGeneric(axesRadius=0.5*a,axesLength=2*a)))
172
173 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
174 #add piston as Mass1D:
175
176 axPiston = list(A @ np.array([L2,0,0]))
177 refPosPiston = list(A @ np.array([0.5*L0+L1-offsetPiston,0,2*a+offZ]))
178 gGraphicsPiston = GraphicsDataCylinder(pAxis=[0,0,0],vAxis=axPiston, radius=2*a, color=color4steelblue)
179 pistonMass = 0.2
180 if i==3:
181 pistonMass *=1.5 #add disturbance into system ...
182 gGraphicsPiston = GraphicsDataCylinder(pAxis=[0,0,0],vAxis=axPiston, radius=2*a, color=color4red)
183
184 n1D1 = mbs.AddNode(Node1D(referenceCoordinates=[0]))
185 oPiston1 = mbs.AddObject(Mass1D(physicsMass = pistonMass,
186 nodeNumber = n1D1,
187 referencePosition=refPosPiston,
188 referenceRotation=A,
189 visualization=VObjectMass1D(graphicsData=[gGraphicsPiston])))
190
191 mbs.CreateSphericalJoint(bodyNumbers=[oPiston1,oRB1], position=[0,0,0],
192 constrainedAxes=[1,1,0], useGlobalFrame=False,
193 jointRadius=1.5*a)
194 #alternatively, using markers and objects:
195 # mRigid11 = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oRB1, localPosition = [ 0.5*L1,0,0])) #connection to piston
196 # mPiston = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oPiston1, localPosition = [ 0,0,0]))
197 # mbs.AddObject(SphericalJoint(markerNumbers=[mRigid11,mPiston],
198 # constrainedAxes=[1,1,0],
199 # visualization=VObjectJointSpherical(jointRadius=1.5*a)))
200
201
202#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
203#drive train
204inertiaDiscBig = InertiaCylinder(density=rho, length=2*a, outerRadius=discBig, axis=2)
205inertiaDiscSmall = InertiaCylinder(density=rho, length=2*a, outerRadius=discSmall, axis=2)
206#print("Jzz=", inertiaDiscBig.GetInertia6D()[2], 0.5*rho*pi*discBig**4*(2*a))
207
208gGraphicsDiscBig0a = GraphicsDataCylinder([0,0,-a],[0,0,2*a], discBig, color4lightred, 64)
209gGraphicsDiscBig0b = GraphicsDataOrthoCube(0,-0.25*a,-a*1.01, discBig, 0.25*a, a*1.01, color4lightgrey) #add something to the cylinder to see rotation
210gGraphicsDiscSmall0a = GraphicsDataCylinder([0,0,-a],[0,0,2*a], discSmall, color4lightred, 32)
211gGraphicsDiscSmall0b = GraphicsDataOrthoCube(0,-0.25*a,-a*1.01, discSmall, 0.25*a, a*1.01, color4lightgrey) #add something to the cylinder to see rotation
212
213#Gear0:
214nDT0 = mbs.AddNode(Node1D(referenceCoordinates = [0]))
215oDT0 = mbs.AddObject(Rotor1D(nodeNumber = nDT0,
216 physicsInertia=inertiaDiscBig.GetInertia6D()[2],
217 referencePosition = [0,0,-2*a],
218 visualization=VObjectRotationalMass1D(graphicsData=[gGraphicsDiscBig0a,gGraphicsDiscBig0b])))
219
220mDT0Rigid = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oDT0, localPosition=[0,0,a]))
221mDT0Coordinate = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nDT0, coordinate=0)) #coordinate for rotation
222
223mbs.AddObject(ObjectJointGeneric(markerNumbers=[mRigid0,mDT0Rigid],
224 constrainedAxes=[0,0,0, 0,0,1],
225 visualization= VObjectJointGeneric(axesRadius=0.6*a,axesLength=1.95*a)))
226
227#Gear1:
228nDT1 = mbs.AddNode(Node1D(referenceCoordinates = [0]))
229oDT1 = mbs.AddObject(Rotor1D(nodeNumber = nDT1,
230 physicsInertia=inertiaDiscSmall.GetInertia6D()[2],
231 referencePosition = [discBig+discSmall,0,-2*a],
232 visualization=VObjectRotationalMass1D(graphicsData=[gGraphicsDiscSmall0a,gGraphicsDiscSmall0b])))
233
234mDT1Coordinate = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nDT1, coordinate=0)) #coordinate for rotation
235mbs.AddObject(CoordinateConstraint(markerNumbers=[mDT0Coordinate,mDT1Coordinate],
236 factorValue1=-discSmall/discBig,
237 visualization=VObjectConnectorCoordinate(show=False)))
238
239#Gear2:
240gGraphicsDiscAxis2 = GraphicsDataCylinder([0,0,-2*a],[0,0,7*a], a, color4grey)
241nDT2 = mbs.AddNode(Node1D(referenceCoordinates = [0]))
242oDT2 = mbs.AddObject(Rotor1D(nodeNumber = nDT2,
243 physicsInertia=inertiaDiscBig.GetInertia6D()[2],
244 referencePosition = [discBig+discSmall,0,-5*a],
245 visualization=VObjectRotationalMass1D(graphicsData=[gGraphicsDiscAxis2,gGraphicsDiscBig0a,gGraphicsDiscBig0b])))
246
247mDT2Coordinate = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nDT2, coordinate=0)) #coordinate for rotation
248mbs.AddObject(CoordinateConstraint(markerNumbers=[mDT1Coordinate,mDT2Coordinate],factorValue1=1,
249 visualization=VObjectConnectorCoordinate(show=False)))
250
251#Gear3:
252gGraphicsDiscAxis3 = GraphicsDataCylinder([0,0,-2*a],[0,0,4*a], a, color4grey)
253nDT3 = mbs.AddNode(Node1D(referenceCoordinates = [0]))
254oDT3 = mbs.AddObject(Rotor1D(nodeNumber = nDT3,
255 physicsInertia=inertiaDiscSmall.GetInertia6D()[2],
256 referencePosition = [(discBig+discSmall)*2,0,-5*a],
257 visualization=VObjectRotationalMass1D(graphicsData=[gGraphicsDiscAxis3,gGraphicsDiscSmall0a,gGraphicsDiscSmall0b])))
258
259mDT3Coordinate = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nDT3, coordinate=0)) #coordinate for rotation
260mbs.AddObject(CoordinateConstraint(markerNumbers=[mDT2Coordinate,mDT3Coordinate],
261 factorValue1=-discSmall/discBig,
262 visualization=VObjectConnectorCoordinate(show=False)))
263
264#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
265#flywheel, connected with MarkerNodeRotationCoordinate:
266gGraphicsDiscFlyWheel0a = GraphicsDataCylinder([0,0,-2*a],[0,0,2*a], discBig, [0.4,0.9,0.4,0.5], 64)
267gGraphicsDiscFlyWheel0b = GraphicsDataOrthoCube(0,-0.25*a,-a*2.01, discBig, 0.25*a, a*0.01, [0.7,0.7,0.7,0.5]) #add something to the cylinder to see rotation
268nDT4 = mbs.AddNode(Node1D(referenceCoordinates = [0]))
269oDT4 = mbs.AddObject(Rotor1D(nodeNumber = nDT4,
270 physicsInertia=5*inertiaDiscBig.GetInertia6D()[2],
271 referencePosition = [0,0,9*a],
272 visualization=VObjectRotationalMass1D(graphicsData=[gGraphicsDiscAxis3,gGraphicsDiscFlyWheel0a,gGraphicsDiscFlyWheel0b])))
273
274mDT4Coordinate = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nDT4, coordinate=0)) #coordinate for rotation
275
276mbs.AddObject(CoordinateConstraint(markerNumbers=[mDT4Coordinate,mRigid0RotationCoordinate],
277 velocityLevel = True, #needed to securly compute multiple rotations
278 visualization=VObjectConnectorCoordinate(show=False)))
279
280sFlyWheelAngVel = mbs.AddSensor(SensorBody(bodyNumber=oDT4, #fileName="solution/sensorFlyWheelAngVel.txt",
281 storeInternal=True,
282 outputVariableType=exu.OutputVariableType.AngularVelocity))
283sFlyWheelAngle = mbs.AddSensor(SensorNode(nodeNumber=nDT4, #fileName="solution/sensorFlyWheelRotation.txt",
284 storeInternal=True,
285 outputVariableType=exu.OutputVariableType.Coordinates))
286
287#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
288#add torque (could also use LoadTorqueVector() on mDT0Rigid)
289def UFLoad(mbs, t, load):
290 if t < 0.25:
291 return load
292 else:
293 return 0
294mbs.AddLoad(LoadCoordinate(markerNumber=mDT3Coordinate, load=100, loadUserFunction=UFLoad))
295
296mbs.Assemble()
297#exu.Print(mbs)
298
299simulationSettings = exu.SimulationSettings() #takes currently set values or default values
300
301tEnd = 0.1
302h=1e-4
303if useGraphics:
304 tEnd = 2
305
306simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
307simulationSettings.timeIntegration.endTime = tEnd
308simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/1000
309simulationSettings.solutionSettings.sensorsWritePeriod = h
310simulationSettings.timeIntegration.verboseMode = 1
311
312simulationSettings.timeIntegration.newton.useModifiedNewton = True
313simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
314simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
315#simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.6 #0.6 works well
316
317simulationSettings.solutionSettings.solutionInformation = "rigid body tests"
318SC.visualizationSettings.nodes.defaultSize = 0.025
319SC.visualizationSettings.nodes.drawNodesAsPoint = False
320SC.visualizationSettings.nodes.showBasis = True
321
322#simulationSettings.displayComputationTime = True
323#simulationSettings.displayStatistics = True
324
325#simulationSettings.solutionSettings.recordImagesInterval = 0.005
326#SC.visualizationSettings.exportImages.saveImageFileName = "images/frame"
327SC.visualizationSettings.window.renderWindowSize = [1920,1080]
328SC.visualizationSettings.openGL.multiSampling = 4
329
330if useGraphics:
331 exu.StartRenderer()
332 if 'lastRenderState' in vars():
333 SC.SetRenderState(lastRenderState) #load last model view
334 mbs.WaitForUserToContinue()
335
336mbs.SolveDynamic(simulationSettings)
337
338#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++
339phiCrank = mbs.GetSensorValues(sCrankAngle)[2]
340phiFlyWheel = mbs.GetSensorValues(sFlyWheelAngle) #scalar coordinate!
341
342phiCrankData = mbs.GetSensorStoredData(sCrankAngle)[:,2:4] #y,z values
343
344exu.Print("phiCrank",phiCrank)
345exu.Print("phiFlyWheel",phiFlyWheel)
346u = phiCrank-phiFlyWheel
347exu.Print("solution of driveTrainTest=", u)
348exudynTestGlobals.testError = u - (0.8813172426357362 - 0.8813173353288565) #2020-05-28: 0.8813172426357362 - 0.8813173353288565
349exudynTestGlobals.testResult = u
350
351if useGraphics:
352 SC.WaitForRenderEngineStopFlag()
353 exu.StopRenderer() #safely close rendering window!
354
355 lastRenderState = SC.GetRenderState() #store model view for next simulation
356
357mbs.PlotSensor(sensorNumbers=[sFlyWheelAngle],
358 components=[0], closeAll=True, offsets=-phiCrankData,
359 labels=['crank angle - flywheel angle'])
360
361if useGraphics:
362 mbs.PlotSensor(sensorNumbers=[sCrankPos, sCrankAngVel, sCrankAngle, sFlyWheelAngVel, sFlyWheelAngle],
363 components=[0,2,2,2,0], markerStyles=['^ ','o ','H ','x','v '],closeAll=True,markerSizes=12,
364 labels=['crank position','crank angular velocity','crank angle','flywheel angular velocity', 'flywheel angle'])