pistonEngine.py
You can view and download this file on Github: pistonEngine.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Create piston engine with variable number of pistons, crank and piston angles;
5# Showing unbalance and harmonics of unbalance
6#
7# Model: Generic piston engine
8#
9# Author: Johannes Gerstmayr
10# Date: 2020-12-20
11#
12# 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.
13#
14# *clean example*
15#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
16
17## import basic libaries
18import exudyn as exu
19from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
20import numpy as np
21import exudyn.graphics as graphics
22from math import sin, cos, asin, acos, pi, exp, log, tan, atan, radians
23
24
25## some simulation parameters for
26createFigures = False
27useLogY = False
28showSolutionViewer = True
29omegaDrive = 2*pi #angular velocity
30tEnd = 2.5+40 #simulation time
31nodeType = exu.NodeType.RotationEulerParameters
32fixedSpeed = True #if false, the speed is given only for first 1 second
33
34
35## a class to store engine parameters and with geometric functions for piston engine
36class EngineParameters:
37 def __init__(self, crankAnglesDegrees=[], pistonAnglesDegrees=[]):
38 #parameters in m, s, kg, rad, ...
39 self.crankAnglesDegrees = crankAnglesDegrees
40 if pistonAnglesDegrees == []:
41 self.pistonAnglesDegrees = list(0*np.array(crankAnglesDegrees))
42 else:
43 self.pistonAnglesDegrees = pistonAnglesDegrees
44
45 crankAngles = pi/180*np.array(crankAnglesDegrees)
46 self.crankAngles = list(crankAngles)
47
48 pistonAngles = pi/180*np.array(self.pistonAnglesDegrees)
49 self.pistonAngles = list(pistonAngles)
50
51 densitySteel = 7850
52 #kinematics & inertia & drawing
53 fZ = 1#0.2
54 self.pistonDistance = 0.08
55 self.pistonMass = 0.5
56 self.pistonLength = 0.05
57 self.pistonRadius = 0.02
58
59 self.conrodLength = 0.1 #X
60 self.conrodHeight = 0.02*fZ#Y
61 self.conrodWidth = 0.02*fZ #Z
62 self.conrodRadius = 0.012*fZ #Z
63
64 self.crankArmLength = 0.04 #X
65 self.crankArmHeight = 0.016 #Y
66 self.crankArmWidth = 0.01*fZ #Z width of arm
67 self.crankBearingWidth = 0.012*fZ #Z
68 self.crankBearingRadius = 0.01
69
70 self.conrodCrankCylLength = 0.024*fZ #Z; length of cylinder (bearing conrod-crank)
71 self.conrodCrankCylRadius = 0.008 #radius of cylinder (bearing conrod-crank)
72
73 self.pistonDistance = self.crankBearingWidth + 2*self.crankArmWidth + self.conrodCrankCylLength #Z distance
74
75 self.inertiaConrod = InertiaCuboid(densitySteel, sideLengths=[self.conrodLength, self.conrodHeight, self.conrodWidth])
76
77 eL = self.Length()
78 #last bearing:
79 densitySteel2 = densitySteel
80 self.inertiaCrank = InertiaCylinder(densitySteel2, self.crankBearingWidth, self.crankBearingRadius, axis=2).Translated([0,0,0.5*eL-0.5*self.crankBearingWidth])
81
82
83
84 for cnt, angle in enumerate(self.crankAngles):
85 A = RotationMatrixZ(angle)
86 zOff = -0.5*eL + cnt*self.pistonDistance
87 arm = InertiaCuboid(densitySteel2, sideLengths=[self.crankArmLength, self.crankArmHeight, self.crankArmWidth])
88 cylCrank = InertiaCylinder(densitySteel2, self.crankBearingWidth, self.crankBearingRadius, axis=2)
89 cylConrod = InertiaCylinder(densitySteel2, self.conrodCrankCylLength, self.conrodCrankCylRadius, axis=2)
90 #add inertias:
91 self.inertiaCrank += cylCrank.Translated([0,0,zOff+self.crankBearingWidth*0.5])
92 self.inertiaCrank += arm.Rotated(A).Translated(A@[self.crankArmLength*0.5,0,zOff+self.crankBearingWidth+self.crankArmWidth*0.5])
93 self.inertiaCrank += cylConrod.Translated(A@[self.crankArmLength,0,zOff+self.crankBearingWidth+self.crankArmWidth+self.conrodCrankCylLength*0.5])
94 self.inertiaCrank += arm.Rotated(A).Translated(A@[self.crankArmLength*0.5,0,zOff+self.crankBearingWidth+self.crankArmWidth*1.5+self.conrodCrankCylLength])
95
96 # self.inertiaCrank = InertiaCylinder(1e-8*densitySteel, length=self.pistonLength,
97 # outerRadius=self.pistonRadius, innerRadius=0.5*self.pistonRadius, axis=2)
98
99 self.inertiaPiston = InertiaCylinder(densitySteel, length=self.pistonLength,
100 outerRadius=self.pistonRadius, innerRadius=0.5*self.pistonRadius, axis=0)
101
102 def Length(self):
103 return self.pistonDistance*len(self.crankAngles) + self.crankBearingWidth
104
105 def MaxDimX(self):
106 return self.crankArmLength + self.conrodLength + self.pistonLength
107
108## compute essential geometrical parameters for slider-crank with crank angle, piston angle, crank length l1 and conrod length l2
109def ComputeSliderCrank(angleCrank, anglePiston, l1, l2):
110 phi1 = angleCrank-anglePiston
111 h = l1*sin(phi1) #height of crank-conrod bearing
112 phi2 = asin(h/l2) #angle of conrod in 2D slider-crank, corotated with piston rotation
113 angleConrod = anglePiston-phi2
114 Acr = RotationMatrixZ(angleConrod)
115 dp = l1*cos(phi1) + l2*cos(phi2) #distance of piston from crank rotation axis
116 return [phi1,phi2, angleConrod, Acr, dp]
117
118
119## function to create multibody system for certain crank and piston configuration
120def CreateEngine(P):
121
122 colorCrank = graphics.color.grey
123 colorConrod = graphics.color.dodgerblue
124 colorPiston = graphics.color.brown[0:3]+[0.5]
125 showJoints = True
126
127 ## set up ground object
128 gravity = [0,-9.81*0,0]
129 eL = P.Length()
130 oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [])))
131 nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
132
133 gEngine = [graphics.Brick(centerPoint=[0,0,0], size=[P.MaxDimX()*2, P.MaxDimX(), eL*1.2],
134 color=[0.6,0.6,0.6,0.1], addEdges=True,
135 edgeColor = [0.8,0.8,0.8,0.3], addFaces=False)]
136
137 ## create rigid body for housing; this body allows to measure support forces and torques
138 oEngine = mbs.CreateRigidBody(referencePosition=[0,0,0],
139 inertia=InertiaCuboid(1000, sideLengths=[1,1,1]), #dummy engine inertia
140 nodeType = nodeType,
141 graphicsDataList = gEngine
142 )
143 nEngine = mbs.GetObjectParameter(oEngine, 'nodeNumber')
144
145 ## create joint between engine and ground to measure forces
146 oEngineJoint = mbs.CreateGenericJoint(itemNumbers=[oEngine, oGround],
147 position=[0,0,0],
148 constrainedAxes=[1,1,1, 1,1,1],
149 show=False)
150
151 ## add sensors for
152 sEngineForce = mbs.AddSensor(SensorObject(objectNumber=oEngineJoint, storeInternal=True,
153 outputVariableType=exu.OutputVariableType.ForceLocal))
154 sEngineTorque = mbs.AddSensor(SensorObject(objectNumber=oEngineJoint, storeInternal=True,
155 outputVariableType=exu.OutputVariableType.TorqueLocal))
156
157 ## loop over all slider-cranks in n-piston engine
158 bConrodList = []
159 bPistonList = []
160 gCrank = []
161 for cnt, angleCrank in enumerate(P.crankAngles):
162 anglePiston = P.pistonAngles[cnt]
163 Ac = RotationMatrixZ(angleCrank)
164 Ap = RotationMatrixZ(anglePiston)
165 [phi1,phi2, angleConrod, Acr, dp] = ComputeSliderCrank(angleCrank, anglePiston, P.crankArmLength, P.conrodLength)
166
167 zOff = -0.5*eL + cnt*P.pistonDistance
168 zAdd = 0
169 if cnt>0: zAdd = P.crankArmWidth
170
171 ### create graphics for crank part
172 gCrank += [graphics.Cylinder(pAxis=[0,0,zOff-zAdd], vAxis=[0,0,P.crankBearingWidth+P.crankArmWidth+zAdd],
173 radius=P.crankBearingRadius, color=graphics.color.red)]
174 ### create graphics for crank arm1
175 arm1 = graphics.Brick([P.crankArmLength*0.5,0,zOff+P.crankArmWidth*0.5+P.crankBearingWidth],
176 size=[P.crankArmLength,P.crankArmHeight,P.crankArmWidth], color=colorCrank)
177 gCrank += [graphics.Move(arm1, [0,0,0], Ac)]
178
179 ### create graphics for conrod bearing
180 gCrank += [graphics.Cylinder(pAxis=Ac@[P.crankArmLength,0,zOff+P.crankBearingWidth+P.crankArmWidth*0],
181 vAxis=[0,0,P.conrodCrankCylLength+2*P.crankArmWidth], radius=P.conrodCrankCylRadius, color=colorCrank)]
182
183 ### create graphics for crank arm2
184 arm2 = graphics.Brick([P.crankArmLength*0.5,0,zOff+P.crankArmWidth*1.5+P.crankBearingWidth+P.conrodCrankCylLength],
185 size=[P.crankArmLength,P.crankArmHeight,P.crankArmWidth],
186 color=colorCrank)
187 gCrank += [graphics.Move(arm2, [0,0,0], Ac)]
188
189 if cnt == len(P.crankAngles)-1:
190 gCrank += [graphics.Cylinder(pAxis=[0,0,zOff+P.crankArmWidth+P.crankBearingWidth+P.conrodCrankCylLength], vAxis=[0,0,P.crankBearingWidth+P.crankArmWidth],
191 radius=P.crankBearingRadius, color=graphics.color.red)]
192
193 #++++++++++++++++++++++++++++++++++++++
194 ### create graphics for conrod
195 gConrod = [ graphics.RigidLink(p0=[-0.5*P.conrodLength, 0, 0], p1=[0.5*P.conrodLength,0,0], axis0= [0,0,1], axis1= [0,0,1],
196 radius= [P.conrodRadius]*2,
197 thickness= P.conrodHeight, width=[P.conrodWidth]*2, color= colorConrod, nTiles= 16)]
198
199 ### create rigid body for conrod
200 bConrod = mbs.CreateRigidBody(inertia = P.inertiaConrod,
201 nodeType = nodeType,
202 referencePosition=Ac@[P.crankArmLength,0,0] + Acr@[0.5*P.conrodLength,0,
203 zOff+P.crankArmWidth+P.crankBearingWidth+0.5*P.conrodCrankCylLength],
204 referenceRotationMatrix=Acr,
205 gravity = gravity,
206 graphicsDataList = gConrod
207 )
208 bConrodList += [bConrod]
209 #++++++++++++++++++++++++++++++++++++++
210 ### create graphics for piston
211 gPiston = [graphics.Cylinder(pAxis=[-P.conrodRadius*0.5,0,0],
212 vAxis=[P.pistonLength,0,0], radius=P.pistonRadius, color=colorPiston)]
213 ### create rigid body for piston
214 bPiston = mbs.CreateRigidBody(inertia = P.inertiaPiston,
215 nodeType = nodeType,
216 referencePosition=Ap@[dp,0,
217 zOff+P.crankArmWidth+P.crankBearingWidth+0.5*P.conrodCrankCylLength],
218 referenceRotationMatrix=Ap,
219 gravity = gravity,
220 graphicsDataList = gPiston
221 )
222 bPistonList += [bPiston]
223
224 ## create rigid body for crankshaft
225 bCrank = mbs.CreateRigidBody(inertia = P.inertiaCrank,
226 nodeType = nodeType,
227 referencePosition=[0,0,0],
228 gravity = gravity,
229 graphicsDataList = gCrank
230 )
231 nCrank = mbs.GetObjectParameter(bCrank, 'nodeNumber')
232
233 ## add sensor for crank angular velocity
234 sCrankAngVel = mbs.AddSensor(SensorNode(nodeNumber=nCrank, storeInternal=True,
235 outputVariableType=exu.OutputVariableType.AngularVelocity))
236
237 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
238 ## create revulute joint between engine and crankshaft
239 oJointCrank = mbs.CreateRevoluteJoint(itemNumbers=[oEngine, bCrank],
240 position=[0,0,-0.5*eL],
241 axis=[0,0,1],
242 show=showJoints,
243 axisRadius=P.crankBearingRadius*1.2,
244 axisLength=P.crankBearingWidth*0.8)
245
246 ## loop over all slider cranks to create joints
247 for cnt, angleCrank in enumerate(P.crankAngles):
248 anglePiston = P.pistonAngles[cnt]
249 Ac = RotationMatrixZ(angleCrank)
250 Ap = RotationMatrixZ(anglePiston)
251 [phi1,phi2, angleConrod, Acr, dp] = ComputeSliderCrank(angleCrank, anglePiston, P.crankArmLength, P.conrodLength)
252
253 zOff = -0.5*eL + cnt*P.pistonDistance
254 #zOff = 0
255
256 ### create revolute joint between crankshaft and conrod
257 oJointCC = mbs.CreateRevoluteJoint(itemNumbers=[bCrank, bConrodList[cnt]],
258 position=Ac@[P.crankArmLength,0,zOff + P.crankBearingWidth+P.crankArmWidth+0.5*P.conrodCrankCylLength],
259 axis=[0,0,1],
260 show = showJoints,
261 axisRadius=P.crankBearingRadius*1.3,
262 axisLength=P.crankBearingWidth*0.8)
263
264 ### create revolute joint between conrod and piston
265 pPiston = Ap@[dp,0,zOff + P.crankBearingWidth+P.crankArmWidth+0.5*P.conrodCrankCylLength]
266 oJointCP = mbs.CreateRevoluteJoint(itemNumbers=[bConrodList[cnt], bPistonList[cnt]],
267 position=pPiston,
268 axis=[0,0,1],
269 show=showJoints,
270 axisRadius=P.crankBearingRadius*1.3,
271 axisLength=P.crankBearingWidth*0.8)
272
273 ### create prismatic joint between piston and engine, using a generic joint
274 mbs.CreateGenericJoint(itemNumbers=[bPistonList[cnt], oEngine],
275 position=[0,0,0],
276 constrainedAxes=[0,1,0, 0,0,1],
277 useGlobalFrame=False,
278 show=True,
279 axesRadius=P.conrodRadius*1.4,
280 axesLength=0.05)
281
282 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
283 ## define user function for crankshaft angle (not used, because velocity level is used):
284 def UFoffset(mbs, t, itemNumber, lOffset):
285 return 0
286
287 ## define user function for crankshaft angular velocity:
288 def UFoffset_t(mbs, t, itemNumber, lOffset): #time derivative of UFoffset
289 return SmoothStep(t, 0, 0.5, 0, omegaDrive)
290
291 ## create coordinate constraint for crankshaft velocity
292 mCrankRotation = mbs.AddMarker(MarkerNodeRotationCoordinate(nodeNumber=nCrank, rotationCoordinate=2))
293 mNodeEngine = mbs.AddMarker(MarkerNodeRotationCoordinate(nodeNumber=nEngine, rotationCoordinate=2))
294 oRotationConstraint = mbs.AddObject(CoordinateConstraint(markerNumbers=[mNodeEngine, mCrankRotation],
295 velocityLevel=True,
296 offsetUserFunction=UFoffset,
297 offsetUserFunction_t=UFoffset_t,
298 visualization=VCoordinateConstraint(show=False)))
299
300 return [oEngine, oEngineJoint, sEngineForce, sEngineTorque, sCrankAngVel, oRotationConstraint, nCrank, bCrank]
301
302## define engine parameters for certain case
303# engine = EngineParameters([0]) #R1
304# engine = EngineParameters([0,180]) #R2
305# engine = EngineParameters([0,180,180,0]) #R4 straight-four engine, Reihen-4-Zylinder
306# engine = EngineParameters([0,90,270,180]) #R4 in different configuration
307engine = EngineParameters([0,180,180,0],[0,180,180,0]) #Boxer 4-piston perfect mass balancing
308
309# engine = EngineParameters([0,120,240]) #R3
310# engine = EngineParameters(list(np.arange(0,5)*144))] #R5
311# engine = EngineParameters([0,120,240,240,120,0]) #R6
312# engine = EngineParameters([0,0,120,120,240,240],[-30,30,-30,30,-30,30]) #V6
313# 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
314
315# engine = EngineParameters([0,90,180,270,270,180,90,360]) #R8
316# engine = EngineParameters([0,0,90,90,270,270,180,180], [-45,45,-45,45, 45,-45,45,-45]) #V8
317
318SC = exu.SystemContainer()
319mbs = SC.AddSystem()
320
321[oEngine, oEngineJoint, sEngineForce, sEngineTorque, sCrankAngVel, oRotationConstraint,
322 nCrank, bCrank] = CreateEngine(engine)
323
324## add prestep user function to turn off drive in case fixedSpeed=False
325def PreStepUF(mbs, t):
326 u = mbs.systemData.GetODE2Coordinates()
327
328 if not fixedSpeed and t >= 1: #at this point, the mechanism runs freely
329 mbs.SetObjectParameter(oRotationConstraint, 'activeConnector', False)
330
331 return True
332
333## add prestep user function
334mbs.SetPreStepUserFunction(PreStepUF)
335
336## assemble system
337mbs.Assemble()
338
339## setup simulation parameters
340stepSize = 0.002
341simulationSettings = exu.SimulationSettings() #takes currently set values or default values
342
343simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
344simulationSettings.timeIntegration.endTime = tEnd
345simulationSettings.timeIntegration.verboseMode = 1
346
347simulationSettings.timeIntegration.realtime.active = True
348
349simulationSettings.solution.file.writePeriod=0.01
350simulationSettings.solution.file.write = True
351simulationSettings.solution.sensors.writePeriod = stepSize
352simulationSettings.solution.file.writeInitialValues = False #otherwise values are duplicated
353simulationSettings.solution.file.name = 'solution/coordinatesSolution.txt'
354
355simulationSettings.timeIntegration.newton.useModifiedNewton = True
356simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
357simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = False
358
359simulationSettings.timeIntegration.generalizedAlpha.lieGroupAddTangentOperator = False
360simulationSettings.linearSolver.solverType=exu.LinearSolverType.EigenSparse
361
362simulationSettings.solution.file.information = "Piston engine"
363
364SC.visualizationSettings.general.graphicsUpdateInterval = 0.01
365SC.visualizationSettings.view0.scene.drawWorldBasis = True
366SC.visualizationSettings.view0.scene.worldBasisSize = 0.1
367
368SC.visualizationSettings.loads.show = False
369SC.visualizationSettings.nodes.show = False
370SC.visualizationSettings.connectors.show = False
371
372SC.visualizationSettings.openGL.multiSampling = 4
373SC.visualizationSettings.openGL.lineWidth = 3
374SC.visualizationSettings.view0.camera.perspective = 0.5
375SC.visualizationSettings.openGL.light0.position = [0.25,1,3,0]
376SC.visualizationSettings.view0.window.renderWindowSize = [1600,1200]
377
378#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
379
380## start visualization and solve
381SC.visualizationSettings.general.autoFitScene = False #use loaded render state
382SC.renderer.Start()
383SC.renderer.RestoreSavedState()
384
385SC.renderer.DoIdleTasks()
386
387mbs.SolveDynamic(simulationSettings)
388
389SC.renderer.Stop() #safely close rendering window!
390
391
392## import plot tools and plot some sensors
393from exudyn.plot import PlotSensor,PlotSensorDefaults
394
395PlotSensor(mbs, closeAll=True)
396PlotSensor(mbs, [sCrankAngVel], components=[2], title='crank speed', sizeInches=[2*6.4,2*4.8])
397PlotSensor(mbs, sEngineForce, components=[0,1,2], title='joint forces')
398PlotSensor(mbs, sEngineTorque, components=[0,1,2], title='joint torques')