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'])