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