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