You can view and download this file on Github: pistonEngine.py
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details: Create piston engine with variable number of pistons, crank and piston angles;
# Showing unbalance and harmonics of unbalance
#
# Model: Generic piston engine
#
# Author: Johannes Gerstmayr
# Date: 2020-12-20
#
# 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.
#
# *clean example*
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## import basic libaries
import exudyn as exu
from exudyn.utilities import *
from math import sin, cos, asin, acos, pi, exp, log, tan, atan, radians
## some simulation parameters for
createFigures = False
useLogY = False
showSolutionViewer = True
omegaDrive = 2*pi #angular velocity
tEnd = 2.5+40 #simulation time
nodeType = exu.NodeType.RotationEulerParameters
fixedSpeed = True #if false, the speed is given only for first 1 second
## a class to store engine parameters and with geometric functions for piston engine
class EngineParameters:
def __init__(self, crankAnglesDegrees=[], pistonAnglesDegrees=[]):
#parameters in m, s, kg, rad, ...
self.crankAnglesDegrees = crankAnglesDegrees
if pistonAnglesDegrees == []:
self.pistonAnglesDegrees = list(0*np.array(crankAnglesDegrees))
else:
self.pistonAnglesDegrees = pistonAnglesDegrees
crankAngles = pi/180*np.array(crankAnglesDegrees)
self.crankAngles = list(crankAngles)
pistonAngles = pi/180*np.array(self.pistonAnglesDegrees)
self.pistonAngles = list(pistonAngles)
densitySteel = 7850
#kinematics & inertia & drawing
fZ = 1#0.2
self.pistonDistance = 0.08
self.pistonMass = 0.5
self.pistonLength = 0.05
self.pistonRadius = 0.02
self.conrodLength = 0.1 #X
self.conrodHeight = 0.02*fZ#Y
self.conrodWidth = 0.02*fZ #Z
self.conrodRadius = 0.012*fZ #Z
self.crankArmLength = 0.04 #X
self.crankArmHeight = 0.016 #Y
self.crankArmWidth = 0.01*fZ #Z width of arm
self.crankBearingWidth = 0.012*fZ #Z
self.crankBearingRadius = 0.01
self.conrodCrankCylLength = 0.024*fZ #Z; length of cylinder (bearing conrod-crank)
self.conrodCrankCylRadius = 0.008 #radius of cylinder (bearing conrod-crank)
self.pistonDistance = self.crankBearingWidth + 2*self.crankArmWidth + self.conrodCrankCylLength #Z distance
self.inertiaConrod = InertiaCuboid(densitySteel, sideLengths=[self.conrodLength, self.conrodHeight, self.conrodWidth])
eL = self.Length()
#last bearing:
densitySteel2 = densitySteel
self.inertiaCrank = InertiaCylinder(densitySteel2, self.crankBearingWidth, self.crankBearingRadius, axis=2).Translated([0,0,0.5*eL-0.5*self.crankBearingWidth])
for cnt, angle in enumerate(self.crankAngles):
A = RotationMatrixZ(angle)
zOff = -0.5*eL + cnt*self.pistonDistance
arm = InertiaCuboid(densitySteel2, sideLengths=[self.crankArmLength, self.crankArmHeight, self.crankArmWidth])
cylCrank = InertiaCylinder(densitySteel2, self.crankBearingWidth, self.crankBearingRadius, axis=2)
cylConrod = InertiaCylinder(densitySteel2, self.conrodCrankCylLength, self.conrodCrankCylRadius, axis=2)
#add inertias:
self.inertiaCrank += cylCrank.Translated([0,0,zOff+self.crankBearingWidth*0.5])
self.inertiaCrank += arm.Rotated(A).Translated(A@[self.crankArmLength*0.5,0,zOff+self.crankBearingWidth+self.crankArmWidth*0.5])
self.inertiaCrank += cylConrod.Translated(A@[self.crankArmLength,0,zOff+self.crankBearingWidth+self.crankArmWidth+self.conrodCrankCylLength*0.5])
self.inertiaCrank += arm.Rotated(A).Translated(A@[self.crankArmLength*0.5,0,zOff+self.crankBearingWidth+self.crankArmWidth*1.5+self.conrodCrankCylLength])
# self.inertiaCrank = InertiaCylinder(1e-8*densitySteel, length=self.pistonLength,
# outerRadius=self.pistonRadius, innerRadius=0.5*self.pistonRadius, axis=2)
self.inertiaPiston = InertiaCylinder(densitySteel, length=self.pistonLength,
outerRadius=self.pistonRadius, innerRadius=0.5*self.pistonRadius, axis=0)
def Length(self):
return self.pistonDistance*len(self.crankAngles) + self.crankBearingWidth
def MaxDimX(self):
return self.crankArmLength + self.conrodLength + self.pistonLength
## compute essential geometrical parameters for slider-crank with crank angle, piston angle, crank length l1 and conrod length l2
def ComputeSliderCrank(angleCrank, anglePiston, l1, l2):
phi1 = angleCrank-anglePiston
h = l1*sin(phi1) #height of crank-conrod bearing
phi2 = asin(h/l2) #angle of conrod in 2D slider-crank, corotated with piston rotation
angleConrod = anglePiston-phi2
Acr = RotationMatrixZ(angleConrod)
dp = l1*cos(phi1) + l2*cos(phi2) #distance of piston from crank rotation axis
return [phi1,phi2, angleConrod, Acr, dp]
## function to create multibody system for certain crank and piston configuration
def CreateEngine(P):
colorCrank = color4grey
colorConrod = color4dodgerblue
colorPiston = color4brown[0:3]+[0.5]
showJoints = True
## set up ground object
gravity = [0,-9.81*0,0]
eL = P.Length()
oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [])))
nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
gEngine = [GraphicsDataOrthoCubePoint(centerPoint=[0,0,0], size=[P.MaxDimX()*2, P.MaxDimX(), eL*1.2],
color=[0.6,0.6,0.6,0.1], addEdges=True,
edgeColor = [0.8,0.8,0.8,0.3], addFaces=False)]
## create rigid body for housing; this body allows to measure support forces and torques
oEngine = mbs.CreateRigidBody(referencePosition=[0,0,0],
inertia=InertiaCuboid(1000, sideLengths=[1,1,1]), #dummy engine inertia
nodeType = nodeType,
graphicsDataList = gEngine
)
nEngine = mbs.GetObjectParameter(oEngine, 'nodeNumber')
## create joint between engine and ground to measure forces
oEngineJoint = mbs.CreateGenericJoint(bodyNumbers=[oEngine, oGround],
position=[0,0,0],
constrainedAxes=[1,1,1, 1,1,1],
show=False)[0]
## add sensors for
sEngineForce = mbs.AddSensor(SensorObject(objectNumber=oEngineJoint, storeInternal=True,
outputVariableType=exu.OutputVariableType.ForceLocal))
sEngineTorque = mbs.AddSensor(SensorObject(objectNumber=oEngineJoint, storeInternal=True,
outputVariableType=exu.OutputVariableType.TorqueLocal))
## loop over all slider-cranks in n-piston engine
bConrodList = []
bPistonList = []
gCrank = []
for cnt, angleCrank in enumerate(P.crankAngles):
anglePiston = P.pistonAngles[cnt]
Ac = RotationMatrixZ(angleCrank)
Ap = RotationMatrixZ(anglePiston)
[phi1,phi2, angleConrod, Acr, dp] = ComputeSliderCrank(angleCrank, anglePiston, P.crankArmLength, P.conrodLength)
zOff = -0.5*eL + cnt*P.pistonDistance
zAdd = 0
if cnt>0: zAdd = P.crankArmWidth
### create graphics for crank part
gCrank += [GraphicsDataCylinder(pAxis=[0,0,zOff-zAdd], vAxis=[0,0,P.crankBearingWidth+P.crankArmWidth+zAdd],
radius=P.crankBearingRadius, color=color4red)]
### create graphics for crank arm1
arm1 = GraphicsDataOrthoCubePoint([P.crankArmLength*0.5,0,zOff+P.crankArmWidth*0.5+P.crankBearingWidth],
size=[P.crankArmLength,P.crankArmHeight,P.crankArmWidth], color=colorCrank)
gCrank += [MoveGraphicsData(arm1, [0,0,0], Ac)]
### create graphics for conrod bearing
gCrank += [GraphicsDataCylinder(pAxis=Ac@[P.crankArmLength,0,zOff+P.crankBearingWidth+P.crankArmWidth*0],
vAxis=[0,0,P.conrodCrankCylLength+2*P.crankArmWidth], radius=P.conrodCrankCylRadius, color=colorCrank)]
### create graphics for crank arm2
arm2 = GraphicsDataOrthoCubePoint([P.crankArmLength*0.5,0,zOff+P.crankArmWidth*1.5+P.crankBearingWidth+P.conrodCrankCylLength],
size=[P.crankArmLength,P.crankArmHeight,P.crankArmWidth],
color=colorCrank)
gCrank += [MoveGraphicsData(arm2, [0,0,0], Ac)]
if cnt == len(P.crankAngles)-1:
gCrank += [GraphicsDataCylinder(pAxis=[0,0,zOff+P.crankArmWidth+P.crankBearingWidth+P.conrodCrankCylLength], vAxis=[0,0,P.crankBearingWidth+P.crankArmWidth],
radius=P.crankBearingRadius, color=color4red)]
#++++++++++++++++++++++++++++++++++++++
### create graphics for conrod
gConrod = [ GraphicsDataRigidLink (p0=[-0.5*P.conrodLength, 0, 0], p1=[0.5*P.conrodLength,0,0], axis0= [0,0,1], axis1= [0,0,1],
radius= [P.conrodRadius]*2,
thickness= P.conrodHeight, width=[P.conrodWidth]*2, color= colorConrod, nTiles= 16)]
### create rigid body for conrod
bConrod = mbs.CreateRigidBody(inertia = P.inertiaConrod,
nodeType = nodeType,
referencePosition=Ac@[P.crankArmLength,0,0] + Acr@[0.5*P.conrodLength,0,
zOff+P.crankArmWidth+P.crankBearingWidth+0.5*P.conrodCrankCylLength],
referenceRotationMatrix=Acr,
gravity = gravity,
graphicsDataList = gConrod
)
bConrodList += [bConrod]
#++++++++++++++++++++++++++++++++++++++
### create graphics for piston
gPiston = [GraphicsDataCylinder(pAxis=[-P.conrodRadius*0.5,0,0],
vAxis=[P.pistonLength,0,0], radius=P.pistonRadius, color=colorPiston)]
### create rigid body for piston
bPiston = mbs.CreateRigidBody(inertia = P.inertiaPiston,
nodeType = nodeType,
referencePosition=Ap@[dp,0,
zOff+P.crankArmWidth+P.crankBearingWidth+0.5*P.conrodCrankCylLength],
referenceRotationMatrix=Ap,
gravity = gravity,
graphicsDataList = gPiston
)
bPistonList += [bPiston]
## create rigid body for crankshaft
bCrank = mbs.CreateRigidBody(inertia = P.inertiaCrank,
nodeType = nodeType,
referencePosition=[0,0,0],
gravity = gravity,
graphicsDataList = gCrank
)
nCrank = mbs.GetObjectParameter(bCrank, 'nodeNumber')
## add sensor for crank angular velocity
sCrankAngVel = mbs.AddSensor(SensorNode(nodeNumber=nCrank, storeInternal=True,
outputVariableType=exu.OutputVariableType.AngularVelocity))
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## create revulute joint between engine and crankshaft
[oJointCrank, mBody0Crank, mBody1Crank] = mbs.CreateRevoluteJoint(bodyNumbers=[oEngine, bCrank],
position=[0,0,-0.5*eL],
axis=[0,0,1],
show=showJoints,
axisRadius=P.crankBearingRadius*1.2,
axisLength=P.crankBearingWidth*0.8)
## loop over all slider cranks to create joints
for cnt, angleCrank in enumerate(P.crankAngles):
anglePiston = P.pistonAngles[cnt]
Ac = RotationMatrixZ(angleCrank)
Ap = RotationMatrixZ(anglePiston)
[phi1,phi2, angleConrod, Acr, dp] = ComputeSliderCrank(angleCrank, anglePiston, P.crankArmLength, P.conrodLength)
zOff = -0.5*eL + cnt*P.pistonDistance
#zOff = 0
### create revolute joint between crankshaft and conrod
[oJointCC, mBody0CC, mBody1CC] = mbs.CreateRevoluteJoint(bodyNumbers=[bCrank, bConrodList[cnt]],
position=Ac@[P.crankArmLength,0,zOff + P.crankBearingWidth+P.crankArmWidth+0.5*P.conrodCrankCylLength],
axis=[0,0,1],
show = showJoints,
axisRadius=P.crankBearingRadius*1.3,
axisLength=P.crankBearingWidth*0.8)
### create revolute joint between conrod and piston
pPiston = Ap@[dp,0,zOff + P.crankBearingWidth+P.crankArmWidth+0.5*P.conrodCrankCylLength]
[oJointCP, mBody0CP, mBody1CP] = mbs.CreateRevoluteJoint(bodyNumbers=[bConrodList[cnt], bPistonList[cnt]],
position=pPiston,
axis=[0,0,1],
show=showJoints,
axisRadius=P.crankBearingRadius*1.3,
axisLength=P.crankBearingWidth*0.8)
### create prismatic joint between piston and engine, using a generic joint
mbs.CreateGenericJoint(bodyNumbers=[bPistonList[cnt], oEngine],
position=[0,0,0],
constrainedAxes=[0,1,0, 0,0,1],
useGlobalFrame=False,
show=True,
axesRadius=P.conrodRadius*1.4,
axesLength=0.05)
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## define user function for crankshaft angle (not used, because velocity level is used):
def UFoffset(mbs, t, itemNumber, lOffset):
return 0
## define user function for crankshaft angular velocity:
def UFoffset_t(mbs, t, itemNumber, lOffset): #time derivative of UFoffset
return SmoothStep(t, 0, 0.5, 0, omegaDrive)
## create coordinate constraint for crankshaft velocity
mCrankRotation = mbs.AddMarker(MarkerNodeRotationCoordinate(nodeNumber=nCrank, rotationCoordinate=2))
mNodeEngine = mbs.AddMarker(MarkerNodeRotationCoordinate(nodeNumber=nEngine, rotationCoordinate=2))
oRotationConstraint = mbs.AddObject(CoordinateConstraint(markerNumbers=[mNodeEngine, mCrankRotation],
velocityLevel=True,
offsetUserFunction=UFoffset,
offsetUserFunction_t=UFoffset_t,
visualization=VCoordinateConstraint(show=False)))
return [oEngine, oEngineJoint, sEngineForce, sEngineTorque, sCrankAngVel, oRotationConstraint, nCrank, bCrank]
## define engine parameters for certain case
# engine = EngineParameters([0]) #R1
# engine = EngineParameters([0,180]) #R2
# engine = EngineParameters([0,180,180,0]) #R4 straight-four engine, Reihen-4-Zylinder
# engine = EngineParameters([0,90,270,180]) #R4 in different configuration
engine = EngineParameters([0,180,180,0],[0,180,180,0]) #Boxer 4-piston perfect mass balancing
# engine = EngineParameters([0,120,240]) #R3
# engine = EngineParameters(list(np.arange(0,5)*144))] #R5
# engine = EngineParameters([0,120,240,240,120,0]) #R6
# engine = EngineParameters([0,0,120,120,240,240],[-30,30,-30,30,-30,30]) #V6
# engine = EngineParameters([0,0,120,120,240,240,240,240,120,120,0,0],[-30,30,-30,30,-30,30,30,-30,30,-30,30,-30]) #V12
# engine = EngineParameters([0,90,180,270,270,180,90,360]) #R8
# engine = EngineParameters([0,0,90,90,270,270,180,180], [-45,45,-45,45, 45,-45,45,-45]) #V8
SC = exu.SystemContainer()
mbs = SC.AddSystem()
[oEngine, oEngineJoint, sEngineForce, sEngineTorque, sCrankAngVel, oRotationConstraint,
nCrank, bCrank] = CreateEngine(engine)
## add prestep user function to turn off drive in case fixedSpeed=False
def PreStepUF(mbs, t):
u = mbs.systemData.GetODE2Coordinates()
if not fixedSpeed and t >= 1: #at this point, the mechanism runs freely
mbs.SetObjectParameter(oRotationConstraint, 'activeConnector', False)
return True
## add prestep user function
mbs.SetPreStepUserFunction(PreStepUF)
## assemble system
mbs.Assemble()
## setup simulation parameters
stepSize = 0.002
simulationSettings = exu.SimulationSettings() #takes currently set values or default values
simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
simulationSettings.timeIntegration.endTime = tEnd
simulationSettings.timeIntegration.verboseMode = 1
simulationSettings.timeIntegration.simulateInRealtime = True
simulationSettings.solutionSettings.solutionWritePeriod=0.01
simulationSettings.solutionSettings.writeSolutionToFile = True
simulationSettings.solutionSettings.sensorsWritePeriod = stepSize
simulationSettings.solutionSettings.writeInitialValues = False #otherwise values are duplicated
simulationSettings.solutionSettings.coordinatesSolutionFileName = 'solution/coordinatesSolution.txt'
simulationSettings.timeIntegration.newton.useModifiedNewton = True
simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = False
simulationSettings.timeIntegration.generalizedAlpha.lieGroupAddTangentOperator = False
simulationSettings.linearSolverType=exu.LinearSolverType.EigenSparse
simulationSettings.solutionSettings.solutionInformation = "Piston engine"
SC.visualizationSettings.general.graphicsUpdateInterval = 0.01
SC.visualizationSettings.general.drawWorldBasis = True
SC.visualizationSettings.general.worldBasisSize = 0.1
SC.visualizationSettings.loads.show = False
SC.visualizationSettings.nodes.show = False
SC.visualizationSettings.connectors.show = False
SC.visualizationSettings.openGL.multiSampling = 4
SC.visualizationSettings.openGL.lineWidth = 3
SC.visualizationSettings.openGL.perspective = 0.5
SC.visualizationSettings.openGL.light0position = [0.25,1,3,0]
SC.visualizationSettings.window.renderWindowSize = [1600,1200]
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## start visualization and solve
SC.visualizationSettings.general.autoFitScene = False #use loaded render state
exu.StartRenderer()
if 'renderState' in exu.sys:
SC.SetRenderState(exu.sys[ 'renderState' ])
mbs.WaitForUserToContinue()
exu.SolveDynamic(mbs, simulationSettings)
exu.StopRenderer() #safely close rendering window!
## import plot tools and plot some sensors
from exudyn.plot import PlotSensor,PlotSensorDefaults
PlotSensor(mbs, closeAll=True)
PlotSensor(mbs, [sCrankAngVel], components=[2], title='crank speed', sizeInches=[2*6.4,2*4.8])
PlotSensor(mbs, sEngineForce, components=[0,1,2], title='joint forces')
PlotSensor(mbs, sEngineTorque, components=[0,1,2], title='joint torques')