You can view and download this file on Github: driveTrainTest.py
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details: Test of Node1D, Mass1D and Rotor1D for drivetrains;4-piston compressor
# Uses different constraints to realize the drivetrain;
# includes joints to connect 3D bodies and 1D drivetrain;
# the crank is elastically supported (except z-direction);
# this makes a direct coupling of the rotation angle to the drivetrain more difficult
#
# Author: Johannes Gerstmayr
# Date: 2020-04-22
#
# 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.
#
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
import exudyn as exu
from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
import exudyn.graphics as graphics #only import if it does not conflict
from exudyn.FEM import *
import numpy as np
useGraphics = True #without test
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
try: #only if called from test suite
from modelUnitTests import exudynTestGlobals #for globally storing test results
useGraphics = exudynTestGlobals.useGraphics
except:
class ExudynTestGlobals:
pass
exudynTestGlobals = ExudynTestGlobals()
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# useGraphics = False
SC = exu.SystemContainer()
mbs = SC.AddSystem()
color = [0.1,0.1,0.8,1]
s = 0.1 #width of cube
sx = 3*s #length of cube/body
background0 = GraphicsDataRectangle(-1,-1,1,1,graphics.color.white)
oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
visualization=VObjectGround(graphicsData= [background0])))
L0 = 0.2 #crank length
L1 = 0.5 #connecting rod length
L2 = 0.1 #cylinder length
a = 0.025 #general width dimension
c = 0.05 #piston radius
discBig = 0.15
discSmall = 0.05
rho=7850 #steel density
kJoint = 1e4 #for elastic support of crankshaft
dJoint = 2e2 #for elastic support of crankshaft
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#add crank:
inertiaCrank = InertiaCuboid(density=rho, sideLengths=[L0,a,a])
omega0 = [0,0,2*pi*100*0] #initial angular velocity of bodies
#ep0 = eulerParameters0
p0 = [0,0,0] #reference position / COM of crank
v0 = [0,0,0] #initial translational velocity
color0= graphics.color.grey
gGraphics0 = graphics.BrickXYZ(-0.5*L0,-a*0.9,-a*0.9,0.5*L0,a*0.9,a*0.9, color0)
gGraphics0b = 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)
gGraphics0c = graphics.Cylinder([0.5*L0,0,-a],[0,0,6*a],a, graphics.color.darkgrey)
gGraphics0d = graphics.Cylinder([-0.5*L0,0,3*a],[0,0,4*a],a, graphics.color.darkgrey)
oRB0 = mbs.CreateRigidBody(inertia=inertiaCrank,
referencePosition=p0,
referenceRotationMatrix=np.eye(3),
initialVelocity=v0,
initialAngularVelocity=omega0,
gravity=[0.,-9.81*0,0.],
graphicsDataList=[gGraphics0,gGraphics0b,gGraphics0c,gGraphics0d])
nRB0 = mbs.GetObject(oRB0)['nodeNumber']
mGround0 = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oGround, localPosition = [0,0,-a]))
mRigid0 = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oRB0, localPosition = [0,0,-a]))
mRigid0RotationCoordinate = mbs.AddMarker(MarkerNodeRotationCoordinate(nodeNumber=nRB0, rotationCoordinate=2))#z-axis
useElasticSupport=True
if not useElasticSupport:
mbs.AddObject(GenericJoint(markerNumbers = [mRigid0,mGround0],
constrainedAxes=[1,1,1, 1,1,0],
visualization= VObjectJointGeneric(axesRadius=a,axesLength=4*a)))
else:
mGround0b = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oGround, localPosition = [0,0,6*a]))
mRigid0b = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oRB0, localPosition = [0,0,6*a]))
mbs.AddObject(GenericJoint(markerNumbers = [mRigid0,mGround0],
constrainedAxes=[0,0,1, 1,1,0],
visualization= VObjectJointGeneric(axesRadius=a,axesLength=4*a)))
mbs.AddObject(CartesianSpringDamper(markerNumbers = [mRigid0,mGround0],
stiffness=[kJoint,kJoint,0],
damping = [dJoint,dJoint,0]))
mbs.AddObject(CartesianSpringDamper(markerNumbers = [mRigid0b,mGround0b],
stiffness=[kJoint,kJoint,0],
damping = [dJoint,dJoint,0]))
sCrankPos = mbs.AddSensor(SensorNode(nodeNumber=nRB0, #fileName="solution/sensorCrankPos.txt",
storeInternal=True,
outputVariableType=exu.OutputVariableType.Position))
sCrankAngVel = mbs.AddSensor(SensorNode(nodeNumber=nRB0, #fileName="solution/sensorCrankAngVel.txt",
storeInternal=True,
outputVariableType=exu.OutputVariableType.AngularVelocity))
sCrankAngle = mbs.AddSensor(SensorNode(nodeNumber=nRB0, #fileName="solution/sensorCrankAngle.txt",
storeInternal=True,
outputVariableType=exu.OutputVariableType.Rotation))
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#add connecting rods and pistons:
for i in range(4):
inertiaConrod = InertiaCuboid(density=rho, sideLengths=[L1,a,a])
phi = i/4*(2*pi) #rotation of subsystem
A = RotationMatrixZ(phi) #transformation of subsystem
A2 = RotationMatrixZ(0)
v2 = np.array([0,0,0])
offsetPiston=0
if i==1 or i==3:
phi2=np.arctan2(0.5*L0,L1) #additional conrod angle
A2 = RotationMatrixZ(phi2) #additional transformation of conrod
v2 = A@np.array([-0.5*L0,-0.25*L0,0])
offsetPiston = L1*(1.-np.cos(phi2)) + 0.5*L0
offZ = 0
if i < 2:
Acrank = RotationMatrixZ(0)
else:
Acrank = RotationMatrixZ(pi)
offZ = 4*a
color1= graphics.color.grey
gGraphics1 = graphics.BrickXYZ(-0.5*L1,-a*0.9,-a*0.9,0.5*L1,a*0.9,a*0.9, color1)
omega1 = [0,0,0] #initial angular velocity of bodies
p1 = [0.5*L0+0.5*L1,0,2*a+offZ] #reference position / COM of crank
p1 = list(v2 + A @ np.array(p1))
v1 = [0,0,0] #initial translational velocity
oRB1 = mbs.CreateRigidBody(inertia=inertiaConrod,
referencePosition=p1,
referenceRotationMatrix=A@A2,
initialAngularVelocity=omega1,
initialVelocity=v1,
gravity=[0.,-9.81*0,0.],
graphicsDataList=[gGraphics1])
locPos0 = list(Acrank @ np.array([0.5*L0,0,a+offZ]))
mbs.CreateGenericJoint(bodyNumbers=[oRB0,oRB1],
position=locPos0, constrainedAxes=[1,1,1, 1,1,0],
useGlobalFrame=False,
axesRadius=0.5*a,axesLength=2*a)
#alternative, using markers and objects:
# mRigid01 = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oRB0, localPosition = locPos0)) #connection crank->conrod
# mRigid10 = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oRB1, localPosition = [-0.5*L1,0,-a]))#connection conrod->crank
# mbs.AddObject(GenericJoint(markerNumbers = [mRigid01,mRigid10],
# constrainedAxes=[1,1,1, 1,1,0],
# visualization= VObjectJointGeneric(axesRadius=0.5*a,axesLength=2*a)))
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#add piston as Mass1D:
axPiston = list(A @ np.array([L2,0,0]))
refPosPiston = list(A @ np.array([0.5*L0+L1-offsetPiston,0,2*a+offZ]))
gGraphicsPiston = graphics.Cylinder(pAxis=[0,0,0],vAxis=axPiston, radius=2*a, color=graphics.color.steelblue)
pistonMass = 0.2
if i==3:
pistonMass *=1.5 #add disturbance into system ...
gGraphicsPiston = graphics.Cylinder(pAxis=[0,0,0],vAxis=axPiston, radius=2*a, color=graphics.color.red)
n1D1 = mbs.AddNode(Node1D(referenceCoordinates=[0]))
oPiston1 = mbs.AddObject(Mass1D(physicsMass = pistonMass,
nodeNumber = n1D1,
referencePosition=refPosPiston,
referenceRotation=A,
visualization=VObjectMass1D(graphicsData=[gGraphicsPiston])))
mbs.CreateSphericalJoint(bodyNumbers=[oPiston1,oRB1], position=[0,0,0],
constrainedAxes=[1,1,0], useGlobalFrame=False,
jointRadius=1.5*a)
#alternatively, using markers and objects:
# mRigid11 = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oRB1, localPosition = [ 0.5*L1,0,0])) #connection to piston
# mPiston = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oPiston1, localPosition = [ 0,0,0]))
# mbs.AddObject(SphericalJoint(markerNumbers=[mRigid11,mPiston],
# constrainedAxes=[1,1,0],
# visualization=VObjectJointSpherical(jointRadius=1.5*a)))
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#drive train
inertiaDiscBig = InertiaCylinder(density=rho, length=2*a, outerRadius=discBig, axis=2)
inertiaDiscSmall = InertiaCylinder(density=rho, length=2*a, outerRadius=discSmall, axis=2)
#print("Jzz=", inertiaDiscBig.GetInertia6D()[2], 0.5*rho*pi*discBig**4*(2*a))
gGraphicsDiscBig0a = graphics.Cylinder([0,0,-a],[0,0,2*a], discBig, graphics.color.lightred, 64)
gGraphicsDiscBig0b = 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
gGraphicsDiscSmall0a = graphics.Cylinder([0,0,-a],[0,0,2*a], discSmall, graphics.color.lightred, 32)
gGraphicsDiscSmall0b = 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
#Gear0:
nDT0 = mbs.AddNode(Node1D(referenceCoordinates = [0]))
oDT0 = mbs.AddObject(Rotor1D(nodeNumber = nDT0,
physicsInertia=inertiaDiscBig.GetInertia6D()[2],
referencePosition = [0,0,-2*a],
visualization=VObjectRotationalMass1D(graphicsData=[gGraphicsDiscBig0a,gGraphicsDiscBig0b])))
mDT0Rigid = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oDT0, localPosition=[0,0,a]))
mDT0Coordinate = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nDT0, coordinate=0)) #coordinate for rotation
mbs.AddObject(ObjectJointGeneric(markerNumbers=[mRigid0,mDT0Rigid],
constrainedAxes=[0,0,0, 0,0,1],
visualization= VObjectJointGeneric(axesRadius=0.6*a,axesLength=1.95*a)))
#Gear1:
nDT1 = mbs.AddNode(Node1D(referenceCoordinates = [0]))
oDT1 = mbs.AddObject(Rotor1D(nodeNumber = nDT1,
physicsInertia=inertiaDiscSmall.GetInertia6D()[2],
referencePosition = [discBig+discSmall,0,-2*a],
visualization=VObjectRotationalMass1D(graphicsData=[gGraphicsDiscSmall0a,gGraphicsDiscSmall0b])))
mDT1Coordinate = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nDT1, coordinate=0)) #coordinate for rotation
mbs.AddObject(CoordinateConstraint(markerNumbers=[mDT0Coordinate,mDT1Coordinate],
factorValue1=-discSmall/discBig,
visualization=VObjectConnectorCoordinate(show=False)))
#Gear2:
gGraphicsDiscAxis2 = graphics.Cylinder([0,0,-2*a],[0,0,7*a], a, graphics.color.grey)
nDT2 = mbs.AddNode(Node1D(referenceCoordinates = [0]))
oDT2 = mbs.AddObject(Rotor1D(nodeNumber = nDT2,
physicsInertia=inertiaDiscBig.GetInertia6D()[2],
referencePosition = [discBig+discSmall,0,-5*a],
visualization=VObjectRotationalMass1D(graphicsData=[gGraphicsDiscAxis2,gGraphicsDiscBig0a,gGraphicsDiscBig0b])))
mDT2Coordinate = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nDT2, coordinate=0)) #coordinate for rotation
mbs.AddObject(CoordinateConstraint(markerNumbers=[mDT1Coordinate,mDT2Coordinate],factorValue1=1,
visualization=VObjectConnectorCoordinate(show=False)))
#Gear3:
gGraphicsDiscAxis3 = graphics.Cylinder([0,0,-2*a],[0,0,4*a], a, graphics.color.grey)
nDT3 = mbs.AddNode(Node1D(referenceCoordinates = [0]))
oDT3 = mbs.AddObject(Rotor1D(nodeNumber = nDT3,
physicsInertia=inertiaDiscSmall.GetInertia6D()[2],
referencePosition = [(discBig+discSmall)*2,0,-5*a],
visualization=VObjectRotationalMass1D(graphicsData=[gGraphicsDiscAxis3,gGraphicsDiscSmall0a,gGraphicsDiscSmall0b])))
mDT3Coordinate = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nDT3, coordinate=0)) #coordinate for rotation
mbs.AddObject(CoordinateConstraint(markerNumbers=[mDT2Coordinate,mDT3Coordinate],
factorValue1=-discSmall/discBig,
visualization=VObjectConnectorCoordinate(show=False)))
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#flywheel, connected with MarkerNodeRotationCoordinate:
gGraphicsDiscFlyWheel0a = graphics.Cylinder([0,0,-2*a],[0,0,2*a], discBig, [0.4,0.9,0.4,0.5], 64)
gGraphicsDiscFlyWheel0b = 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
nDT4 = mbs.AddNode(Node1D(referenceCoordinates = [0]))
oDT4 = mbs.AddObject(Rotor1D(nodeNumber = nDT4,
physicsInertia=5*inertiaDiscBig.GetInertia6D()[2],
referencePosition = [0,0,9*a],
visualization=VObjectRotationalMass1D(graphicsData=[gGraphicsDiscAxis3,gGraphicsDiscFlyWheel0a,gGraphicsDiscFlyWheel0b])))
mDT4Coordinate = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nDT4, coordinate=0)) #coordinate for rotation
mbs.AddObject(CoordinateConstraint(markerNumbers=[mDT4Coordinate,mRigid0RotationCoordinate],
velocityLevel = True, #needed to securly compute multiple rotations
visualization=VObjectConnectorCoordinate(show=False)))
sFlyWheelAngVel = mbs.AddSensor(SensorBody(bodyNumber=oDT4, #fileName="solution/sensorFlyWheelAngVel.txt",
storeInternal=True,
outputVariableType=exu.OutputVariableType.AngularVelocity))
sFlyWheelAngle = mbs.AddSensor(SensorNode(nodeNumber=nDT4, #fileName="solution/sensorFlyWheelRotation.txt",
storeInternal=True,
outputVariableType=exu.OutputVariableType.Coordinates))
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#add torque (could also use LoadTorqueVector() on mDT0Rigid)
def UFLoad(mbs, t, load):
if t < 0.25:
return load
else:
return 0
mbs.AddLoad(LoadCoordinate(markerNumber=mDT3Coordinate, load=100, loadUserFunction=UFLoad))
mbs.Assemble()
#exu.Print(mbs)
simulationSettings = exu.SimulationSettings() #takes currently set values or default values
tEnd = 0.1
h=1e-4
if useGraphics:
tEnd = 2
simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
simulationSettings.timeIntegration.endTime = tEnd
simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/1000
simulationSettings.solutionSettings.sensorsWritePeriod = h
simulationSettings.timeIntegration.verboseMode = 1
simulationSettings.timeIntegration.newton.useModifiedNewton = True
simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
#simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.6 #0.6 works well
simulationSettings.solutionSettings.solutionInformation = "rigid body tests"
SC.visualizationSettings.nodes.defaultSize = 0.025
SC.visualizationSettings.nodes.drawNodesAsPoint = False
SC.visualizationSettings.nodes.showBasis = True
#simulationSettings.displayComputationTime = True
#simulationSettings.displayStatistics = True
#simulationSettings.solutionSettings.recordImagesInterval = 0.005
#SC.visualizationSettings.exportImages.saveImageFileName = "images/frame"
SC.visualizationSettings.window.renderWindowSize = [1920,1080]
SC.visualizationSettings.openGL.multiSampling = 4
if useGraphics:
exu.StartRenderer()
if 'lastRenderState' in vars():
SC.SetRenderState(lastRenderState) #load last model view
mbs.WaitForUserToContinue()
mbs.SolveDynamic(simulationSettings)
#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++
phiCrank = mbs.GetSensorValues(sCrankAngle)[2]
phiFlyWheel = mbs.GetSensorValues(sFlyWheelAngle) #scalar coordinate!
phiCrankData = mbs.GetSensorStoredData(sCrankAngle)[:,2:4] #y,z values
exu.Print("phiCrank",phiCrank)
exu.Print("phiFlyWheel",phiFlyWheel)
u = phiCrank-phiFlyWheel
exu.Print("solution of driveTrainTest=", u)
exudynTestGlobals.testError = u - (0.8813172426357362 - 0.8813173353288565) #2020-05-28: 0.8813172426357362 - 0.8813173353288565
exudynTestGlobals.testResult = u
if useGraphics:
SC.WaitForRenderEngineStopFlag()
exu.StopRenderer() #safely close rendering window!
lastRenderState = SC.GetRenderState() #store model view for next simulation
mbs.PlotSensor(sensorNumbers=[sFlyWheelAngle],
components=[0], closeAll=True, offsets=-phiCrankData,
labels=['crank angle - flywheel angle'])
if useGraphics:
mbs.PlotSensor(sensorNumbers=[sCrankPos, sCrankAngVel, sCrankAngle, sFlyWheelAngVel, sFlyWheelAngle],
components=[0,2,2,2,0], markerStyles=['^ ','o ','H ','x','v '],closeAll=True,markerSizes=12,
labels=['crank position','crank angular velocity','crank angle','flywheel angular velocity', 'flywheel angle'])