NGsolvePistonEngine.py
You can view and download this file on Github: NGsolvePistonEngine.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: generate a piston engine with finite element mesh
5# created with NGsolve and with variable number of pistons
6#
7# Author: Johannes Gerstmayr
8# Date: 2020-06-12
9#
10# 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.
11#
12#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
13
14
15import sys
16import exudyn as exu
17
18from exudyn.itemInterface import *
19from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
20import exudyn.graphics as graphics #only import if it does not conflict
21from exudyn.rigidBodyUtilities import *
22from exudyn.FEM import *
23
24import time
25
26try: #if installed, may help to improve performance (depending on solver)
27 import mkl
28 mklThreads = 8
29 mkl.set_num_threads(mklThreads)
30 exu.Print('using',mklThreads,'mkl threads ...')
31except: pass
32
33from ngsolve import *
34from netgen.geom2d import unit_square
35
36#import netgen.libngpy as libng
37
38netgenDrawing = False #set true, to show geometry and mesh in NETGEN
39#if netgenDrawing, uncomment the following line and execute in external terminal, not in spyder (see preferences "Run"):
40#import netgen.gui
41
42from netgen.csg import *
43
44import numpy as np
45import timeit
46
47verbose = True
48meshSize = 0.005*4 #fast: 0.005*4 with order 1; standard:0.005; h=0.004 with meshOrder 2 gives 258k unknowns; fine: 0.0011: memory limit (96GB) for NGsolve; < 0.0015 makes problems with scipy eigensolver
49meshOrder = 1 #2 for stresses!
50showStresses = True #may take very long for large number of modes/nodes
51saveMesh = False
52
53#++++++++++++++++++++++++++++++++++++
54#helper functions (copied from EXUDYN):
55def RotationMatrixZ(angleRad):
56 return np.array([ [np.cos(angleRad),-np.sin(angleRad), 0],
57 [np.sin(angleRad), np.cos(angleRad), 0],
58 [0, 0, 1] ]);
59
60def VAdd(v0, v1):
61 if len(v0) != len(v1): print("ERROR in VAdd: incompatible vectors!")
62 n = len(v0)
63 v = [0]*n
64 for i in range(n):
65 v[i] = v0[i]+v1[i]
66 return v
67
68def VSub(v0, v1):
69 if len(v0) != len(v1): print("ERROR in VSub: incompatible vectors!")
70 n = len(v0)
71 v = [0]*n
72 for i in range(n):
73 v[i] = v0[i]-v1[i]
74 return v
75
76def NormL2(vector):
77 value = 0
78 for x in vector:
79 value += x**2
80 return value**0.5
81
82def Normalize(v):
83 v2=[0]*len(v)
84
85 fact = NormL2(v)
86 fact = 1./fact
87 for i in range(len(v2)):
88 v2[i]=fact*v[i]
89 return v2
90#++++++++++++++++++++++++++++++++++++
91startTotal = timeit.default_timer()
92#parameters
93
94#crank:
95b1 = 0.012 #width of journal bearing
96r1 = 0.012 #radius of journal bearing
97dk = 0.015 #crank arm width (z)
98bk = 0.032 #crank arm size (y)
99
100l3 = 0.030
101l4 = 0.040
102#l4x= 0.005 #offset of counterweight
103lk = 0.030 #l4*0.5+l3 #crank arm length (x)
104bm = 0.065
105dBevel = dk*0.5
106#shaft:
107r0 = 0.012 #0.012
108d0 = 0.020 #shaft length at left/right support
109d1 = 0.012 #shaft length at intermediate support
110
111#distance rings:
112db = 0.002 #width of distance ring
113rdb0 = r0+db #total radius of distance ring, shaft
114rdb1 = r1+db #total radius of distance ring, crank
115
116#conrod:
117bc = 0.024 #height of conrod
118dc = 0.012 #width of conrod
119lc = 0.080 #length of conrod (axis-axis)
120r1o= r1+0.006 #outer radius of conrod at crank joint
121r2 = 0.008 #radius of piston journal bearing
122r2o= r2+0.006 #outer radius of conrod at piston joint
123
124cylOffZ=0.010 #z-offset of cylinder cut out of conrod
125cylR = 0.008 #radius of cylinder cut out of conrod
126
127angC = 4*np.pi/180
128
129#piston:
130dpb = r2o-0.000 #axis inside piston
131r2p = r2o+0.004 #0.018
132lp = 0.034
133bp = 0.050
134lpAxis = dc+2*db
135lOffCut = 0.011 #offset for cutout of big cylinder
136
137#total length of one segment:
138lTotal = db+dk+db+b1+db+dk+db+d1
139
140#eps
141eps = 5e-4 #added to faces, to avoid CSG-problems
142
143#++++++++++++++++++++++++++++++++++++
144#points
145pLB = [0 ,0,-d0]
146p0B = [0 ,0,0]
147p1B = [0 ,0,db]
148#p2B = [0, 0,db+dk]
149p21B =[lk,0,db+dk]
150p31B = [lk,0,db+dk+db]
151p41B = [lk,0,db+dk+db+b1]
152p51B =[lk,0,db+dk+db+b1+db]
153p6B = [0 ,0,db+dk+db+b1+db+dk]
154p7B = [0 ,0,db+dk+db+b1+db+dk+db]
155p8B = [0 ,0,lTotal]
156
157def CSGcylinder(p0,p1,r):
158 v = VSub(p1,p0)
159 v = Normalize(v)
160 cyl = Cylinder(Pnt(p0[0],p0[1],p0[2]), Pnt(p1[0],p1[1],p1[2]),
161 r) * Plane(Pnt(p0[0],p0[1],p0[2]), Vec(-v[0],-v[1],-v[2])) * Plane(Pnt(p1[0],p1[1],p1[2]), Vec(v[0],v[1],v[2]))
162 return cyl
163
164def CSGcube(pCenter,size):
165 s2 = [0.5*size[0],0.5*size[1],0.5*size[2]]
166 p0 = VSub(pCenter,s2)
167 p1 = VAdd(pCenter,s2)
168 brick = OrthoBrick(Pnt(p0[0],p0[1],p0[2]),Pnt(p1[0],p1[1],p1[2]))
169 return brick
170
171
172#transform points
173def TransformCrank(p, zOff, zRot):
174 p2 = RotationMatrixZ(zRot) @ p
175 pOff=[0,0,zOff]
176 return VAdd(p2,pOff)
177
178#cube only in XY-plane, z infinite
179def CSGcubeXY(pCenter,sizeX,sizeY,ex,ey):
180 #print("pCenter=",pCenter)
181 pl1 = Plane(Pnt(pCenter[0]-0.5*sizeX*ex[0],pCenter[1]-0.5*sizeX*ex[1],0),Vec(-ex[0],-ex[1],-ex[2]))
182 pl2 = Plane(Pnt(pCenter[0]+0.5*sizeX*ex[0],pCenter[1]+0.5*sizeX*ex[1],0),Vec( ex[0], ex[1], ex[2]))
183
184 pl3 = Plane(Pnt(pCenter[0]-0.5*sizeY*ey[0],pCenter[1]-0.5*sizeY*ey[1],0),Vec(-ey[0],-ey[1],-ey[2]))
185 pl4 = Plane(Pnt(pCenter[0]+0.5*sizeY*ey[0],pCenter[1]+0.5*sizeY*ey[1],0),Vec( ey[0], ey[1], ey[2]))
186
187 return pl1*pl2*pl3*pl4
188
189
190#create one crank face at certain z-offset and rotation; side=1: left, side=-1: right
191def GetCrankFace(zOff, zRot, side=1):
192 ex = RotationMatrixZ(zRot) @ [1,0,0]
193 ey = RotationMatrixZ(zRot) @ [0,1,0]
194 #print("zOff=",zOff, "zRot=", zRot, "side=", side,"ex=", ex)
195 pLeft = [0,0,zOff]
196 pRight = [0,0,zOff+dk]
197 pMid = [0,0,zOff+0.5*dk]
198
199 pcLeft=VAdd(pLeft,lk*ex)
200 pcRight=VAdd(pRight,lk*ex)
201 f=0.5**0.5
202 cyl1pl = Plane(Pnt(pcLeft[0],pcLeft[1],pcLeft[2]+0.5*dk-side*dk),Vec(f*ex[0],f*ex[1],f*ex[2]-side*f))
203 cyl1 = Cylinder(Pnt(pcLeft[0],pcLeft[1],pcLeft[2]-1), Pnt(pcRight[0],pcRight[1],pcRight[2]+1), 0.5*bk)*cyl1pl
204
205 #cone2 = Cylinder(Pnt(pcLeft[0],pcLeft[1],pcLeft[2]-1), Pnt(pcRight[0],pcRight[1],pcRight[2]+1), lk+l4)
206 cone2 = Cone(Pnt(pcLeft[0],pcLeft[1],pcLeft[2]-side*dBevel+0.5*dk), Pnt(pcLeft[0],pcLeft[1],pcLeft[2]+side*dBevel+0.5*dk), lk+l4-1.5*dBevel, lk+l4-0.5*dBevel)
207 cube1 = CSGcubeXY(VAdd(pMid,0.49*l3*ex),1.02*l3,bk,ex,ey) #make l3 a little longer, to avoid bad edges
208 cube2 = CSGcubeXY(VAdd(pMid,-0.5*l4*ex),1.0*l4,bm,ex,ey)*cone2
209
210 pc3a = VAdd(pLeft,0.*l3*ex+(0.5*bk+0.4*l3)*ey)
211 cyl3a = Cylinder(Pnt(pc3a[0],pc3a[1],pc3a[2]-1), Pnt(pc3a[0],pc3a[1],pc3a[2]+1), 0.42*l3)
212 pc3b = VAdd(pLeft,0.*l3*ex+(-0.5*bk-0.4*l3)*ey)
213 cyl3b = Cylinder(Pnt(pc3b[0],pc3b[1],pc3b[2]-1), Pnt(pc3b[0],pc3b[1],pc3b[2]+1), 0.42*l3)
214 #cube3a = (CSGcubeXY(VAdd(pMid,0.26*l3*ex+(0.5*bk+0.26*l3)*ey),0.5*l3,0.5*l3,ex,ey)-cyl3a)
215
216 return ((cube1+cube2+cyl1)-(cyl3a+cyl3b))*Plane(Pnt(0,0,pLeft[2]),Vec(0,0,-1))*Plane(Pnt(0,0,pRight[2]),Vec(0,0,1))
217 #return (cube1+cube2+cyl1)*Plane(Pnt(0,0,pLeft[2]),Vec(0,0,-1))*Plane(Pnt(0,0,pRight[2]),Vec(0,0,1))
218
219#generate one crank, rotated around z-axis in radiant
220def GenerateCrank(zOff, zRot):
221 pL = TransformCrank(pLB,zOff, zRot)
222 p0 = TransformCrank(p0B,zOff, zRot)
223 p1 = TransformCrank(p1B,zOff, zRot)
224
225 p21 = TransformCrank(p21B,zOff, zRot)
226 p31 = TransformCrank(p31B,zOff, zRot)
227 p41 = TransformCrank(p41B,zOff, zRot)
228 p51 = TransformCrank(p51B,zOff, zRot)
229
230 p6 = TransformCrank(p6B,zOff, zRot)
231 p7 = TransformCrank(p7B,zOff, zRot)
232 p8 = TransformCrank(p8B,zOff, zRot)
233
234 crank0 = CSGcylinder(pL,[p0[0],p0[1],p0[2]+eps],r0)
235 crank1 = CSGcylinder(p0,[p1[0],p1[1],p1[2]+eps],rdb0)
236
237 #conrod bearing:
238 crank3 = CSGcylinder([p21[0],p21[1],p21[2]-eps],p31,rdb1)
239 crank7 = CSGcylinder(p31,p41,r1)
240 crank8 = CSGcylinder(p41,[p51[0],p51[1],p51[2]+eps],rdb1)
241
242 crank9 = CSGcylinder([p6[0],p6[1],p6[2]-eps],p7,rdb0)
243 crank10 = CSGcylinder([p7[0],p7[1],p7[2]-eps],p8,r0)
244
245 #return crank0+crank1+crank3+crank4+crank5+crank6+crank7+crank8+crank4b+crank5b+crank6b+crank9+crank10
246 if zOff==0:#add first shaft
247 crank1 = crank1+crank0
248 return crank1+GetCrankFace(db+zOff,zRot,1)+crank3+crank7+crank8+GetCrankFace(db+2*db+dk+b1+zOff,zRot,-1)+crank10+crank9
249
250
251geoCrank = CSGeometry()
252
253#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
254#choose configuration for crankshaft:
255#crankConfig = [0] #1-piston
256#crankConfig = [np.pi/2] #1-piston
257#crankConfig = [0,np.pi] #2-piston
258#crankConfig = [0,np.pi*2./3.,2.*np.pi*2./3.] #3-piston
259#crankConfig = [0,np.pi,np.pi,0] #4-piston
260crankConfig = [0,np.pi*2./3.,2.*np.pi*2./3.,2.*np.pi*2./3.,np.pi*2./3.,0] #6-piston
261#crankConfig = crankConfig*2 #12-piston
262
263nPistons = len(crankConfig)
264
265crank = GenerateCrank(0, crankConfig[0])
266zPos = lTotal
267for i in range(len(crankConfig)-1):
268 angle = crankConfig[i+1]
269 crank += GenerateCrank(zPos, angle)
270 zPos += lTotal
271
272# crank = (GenerateCrank(0, 0) + GenerateCrank(lTotal, np.pi*2./3.) + GenerateCrank(2*lTotal, np.pi*2.*2./3.)+
273# GenerateCrank(3*lTotal, np.pi*2.*2./3.) + GenerateCrank(4*lTotal, np.pi*2./3.))
274
275geoCrank.Add(crank)
276
277#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
278#conrod model:
279def GenerateConrod(zOff):
280 ey0 = [0,1,0] #top/bottom face vector of conrod
281 ey1 = [0,-1,0]
282
283 ex0 = [1,0,0] #top/bottom face vector of conrod
284 ex1 = [1,0,0]
285
286 ey0 = RotationMatrixZ(-angC)@ey0
287 ey1 = RotationMatrixZ(angC)@ey1
288 ex0 = RotationMatrixZ(-angC)@ex0
289 ex1 = RotationMatrixZ(angC)@ex1
290
291
292 pl1 = Plane(Pnt(0, 0.5*bc,0),Vec(ey0[0],ey0[1],ey0[2]))
293 pl2 = Plane(Pnt(0,-0.5*bc,0),Vec(ey1[0],ey1[1],ey1[2]))
294
295 pl3 = Plane(Pnt(-0.5*lc,0,0),Vec(-1,0,0))
296 pl4 = Plane(Pnt( 0.5*lc,0,0),Vec( 1,0,0))
297
298 pl5 = Plane(Pnt( 0,0,-0.5*dc+zOff),Vec( 0,0,-1))
299 pl6 = Plane(Pnt( 0,0, 0.5*dc+zOff),Vec( 0,0, 1))
300
301
302 cylC1 = Cylinder(Pnt(-0.5*lc,0,-1), Pnt(-0.5*lc,0,1), r1)
303 #cylC1o = Cylinder(Pnt(-0.5*lc,0,-1), Pnt(-0.5*lc,0,1), r1o)
304 cylC1o = Sphere(Pnt(-0.5*lc,0,zOff), r1o) #in fact is a sphere
305
306 cylC2 = Cylinder(Pnt( 0.5*lc,0,-1), Pnt( 0.5*lc,0,1), r2)
307 #cylC2o = Cylinder(Pnt(0.5*lc,0,-1), Pnt( 0.5*lc,0,1), r2o)
308 cylC2o = Sphere(Pnt(0.5*lc,0,zOff), r2o) #in fact is a sphere
309
310 cylSideA = (Cylinder(Pnt(-0.5*lc+r1o,0,cylOffZ+zOff), Pnt(0.5*lc-r2o,0,cylOffZ+zOff), cylR)*
311 Plane(Pnt(-0.5*lc+r1o-0.002,0,0),Vec(-1,0,0))*
312 Plane(Pnt( 0.5*lc-r2o+0.002,0,0),Vec( 1,0,0)))
313
314 cylSideB = (Cylinder(Pnt(-0.5*lc+r1o,0,-cylOffZ+zOff), Pnt(0.5*lc-r2o,0,-cylOffZ+zOff), cylR)*
315 Plane(Pnt(-0.5*lc+r1o-0.002,0,0),Vec(-1,0,0))*
316 Plane(Pnt( 0.5*lc-r2o+0.002,0,0),Vec( 1,0,0)))
317
318
319 return ((pl1*pl2*pl3*pl4+cylC1o+cylC2o)-cylC1-cylC2)*pl5*pl6-cylSideA-cylSideB
320 #return pl1*pl2*pl3*pl4*pl5*pl6
321
322geoConrod = CSGeometry()
323conrod = GenerateConrod(0)#db+dk+db+0.5*b1
324geoConrod.Add(conrod)
325
326# if netgenDrawing:
327# Draw(geoCrank)
328
329#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
330#conrod model:
331def GeneratePiston(zOff):
332 p0 = [-dpb,0,zOff]
333 p1 = [-dpb+lp,0,zOff]
334 cylPo = CSGcylinder(p0, p1, 0.5*bp) #piston outside
335 cylPaxis= CSGcylinder([0,0,-0.5*lpAxis-eps+zOff], [0,0, 0.5*lpAxis+eps+zOff], r2) #piston axis
336 cylPaxis0= CSGcylinder([0,0,-0.5*lpAxis-eps+zOff], [0,0,-0.5*lpAxis+db+zOff], r2+db) #piston axis
337 cylPaxis1= CSGcylinder([0,0, 0.5*lpAxis-db+zOff], [0,0, 0.5*lpAxis+eps+zOff], r2+db) #piston axis
338 cylPin = CSGcylinder([0,0,-0.5*lpAxis+zOff], [0,0, 0.5*lpAxis+zOff], r2p) #piston inner cutout
339
340 #box = CSGcube([0,0,zOff], [dpb+r2p,2*(r2p),lpAxis])
341 box = CSGcube([-0.5*dpb,0,zOff], [dpb,2*(r2p)-0.002,lpAxis-0.000])
342
343 cylCut = CSGcylinder([-(l4+l3+lOffCut),0,-bp+zOff], [-(l4+l3+lOffCut),0, bp+zOff], l4+l3) #piston inner cutout
344
345 return (cylPo-box-cylCut-cylPin)+cylPaxis+cylPaxis0+cylPaxis1
346
347geoPiston = CSGeometry()
348piston = GeneratePiston(0)#db+dk+db+0.5*b1
349geoPiston.Add(piston)
350
351if verbose: print("Generate meshes ...")
352#do meshing, if geometry is successful
353if True:
354 ngMeshCrank = geoCrank.GenerateMesh(maxh=meshSize)
355 meshCrank = Mesh(ngMeshCrank)
356 meshCrank.Curve(1)
357 if netgenDrawing:
358 Draw(meshCrank)
359
360 if False:
361 #save mesh to file:
362 meshCrank.ngmesh.Export('testData/crankshaft.mesh','Neutral Format')
363
364if True:
365 ngMeshConrod = geoConrod.GenerateMesh(maxh=meshSize) #in videos 0.003
366 meshConrod = Mesh(ngMeshConrod)
367 meshConrod.Curve(1)
368 if netgenDrawing:
369 Draw(meshConrod)
370 if False:
371 meshConrod.ngmesh.Export('testData/conrod.mesh','Neutral Format')
372 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++
373
374if True:
375 ngMeshPiston = geoPiston.GenerateMesh(maxh=meshSize+0.001*0)
376 meshPiston = Mesh(ngMeshPiston)
377 meshPiston.Curve(1)
378 if netgenDrawing:
379 Draw(meshPiston)
380 if False:
381 meshPiston.ngmesh.Export('testData/piston.mesh','Neutral Format')
382 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++
383
384#here starts the EXUDYN part
385if True:
386 SC = exu.SystemContainer()
387 mbs = SC.AddSystem()
388
389 #crankshaft and piston mechanical parameters:
390 density = 7850
391 youngsModulus = 2.1e11 *1e-1
392 poissonsRatio = 0.3
393 fRotorStart = 20 #initial revolutions per second, only crankshaft
394
395 totalFEcoordinates = 0 #accumulated FE-mesh coordinates
396 #%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
397 #import crankshaft mesh into EXUDYN FEMinterface
398 femCrank = FEMinterface()
399 eigenModesNGsolve=True
400 nModes=8
401
402 [bfM, bfK, fes] = femCrank.ImportMeshFromNGsolve(meshCrank, density, youngsModulus, poissonsRatio,
403 verbose = True, meshOrder = meshOrder)
404 # computeEigenmodes=eigenModesNGsolve, excludeRigidBodyModes = 6,
405 # numberOfModes = nModes, maxEigensolveIterations=20)
406
407 nModes = 20
408 excludeRigidBodyModes = 6
409 if verbose: print("number of coordinates crank =", femCrank.NumberOfCoordinates())
410 if verbose: print("Compute eigenmodes crank ....")
411
412 if not eigenModesNGsolve:
413 startCrank = timeit.default_timer()
414 femCrank.ComputeEigenmodes(nModes, excludeRigidBodyModes = excludeRigidBodyModes, useSparseSolver = True)
415 stopCrank = timeit.default_timer()
416 print("\ncrank eigen analysis time=", stopCrank-startCrank)
417 else:
418 start_time = time.time()
419 femCrank.ComputeEigenmodesNGsolve(bfM, bfK, nModes=nModes,
420 excludeRigidBodyModes=excludeRigidBodyModes, maxEigensolveIterations=20)
421 print("NGsolve mode computation needed %.3f seconds" % (time.time() - start_time))
422
423 totalFEcoordinates+=femCrank.NumberOfCoordinates()
424 print("eigen freq. crank=", femCrank.GetEigenFrequenciesHz()[0:nModes])
425
426 #+++++++++++++++++++++++++++++++++++++++++++++++++++++
427 #compute stress modes:
428 SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
429 mat = KirchhoffMaterial(youngsModulus, poissonsRatio, density)
430 varType = exu.OutputVariableType.DisplacementLocal
431 #varType = exu.OutputVariableType.StrainLocal
432 if showStresses:
433 print("ComputePostProcessingModes femCrank ... ")
434 start_time = time.time()
435 varType = exu.OutputVariableType.StressLocal
436 femCrank.ComputePostProcessingModesNGsolve(fes, material=mat,
437 outputVariableType=varType)
438 print("--- %s seconds ---" % (time.time() - start_time))
439
440 SC.visualizationSettings.contour.outputVariable = varType
441
442 #print("Create CMS object and matrices ....")
443 cmsCrank = ObjectFFRFreducedOrderInterface(femCrank)
444
445 objFFRFcrank = cmsCrank.AddObjectFFRFreducedOrderWithUserFunctions(exu, mbs,
446 positionRef=[0,0,0],
447 eulerParametersRef=eulerParameters0,
448 initialVelocity=[0,0,0], initialAngularVelocity=[0,0,1*fRotorStart*2*pi],
449 gravity = [0,-0*9.81,0],
450 color=[0.1,0.9,0.1,1.])
451 mbs.SetObjectParameter(objFFRFcrank['oFFRFreducedOrder'],'VshowNodes',False)
452
453
454 if False:#animate eigenmodes of crankshaft
455 from exudyn.interactive import AnimateModes
456 mbs.Assemble()
457
458 SC.visualizationSettings.general.textSize = 16 #30 for cover figure
459 SC.visualizationSettings.general.useGradientBackground = True
460 SC.visualizationSettings.openGL.lineWidth = 2
461 SC.visualizationSettings.openGL.showFaceEdges = True
462 SC.visualizationSettings.openGL.showFaces = True
463 SC.visualizationSettings.openGL.multiSampling = 4
464 SC.visualizationSettings.nodes.show = False
465 SC.visualizationSettings.window.renderWindowSize = [1600,1080]
466
467 SC.visualizationSettings.contour.outputVariableComponent = 0
468
469 SC.visualizationSettings.general.autoFitScene=False
470
471 AnimateModes(SC, mbs, 1, period=0.2)
472 exit()
473
474 #%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
475 #import conrod and piston mesh into EXUDYN FEMinterface and compute eigenmodes
476 nModes = 8
477 excludeRigidBodyModes = 6
478 femConrod = FEMinterface()
479 # femConrod.ImportMeshFromNGsolve(meshConrod, density, youngsModulus, poissonsRatio, verbose = False)
480 [bfM, bfK, fes] = femConrod.ImportMeshFromNGsolve(meshConrod, density, youngsModulus, poissonsRatio,
481 verbose = False, meshOrder = meshOrder)
482 # computeEigenmodes=eigenModesNGsolve, excludeRigidBodyModes = 6,
483 # numberOfModes = nModes, maxEigensolveIterations=20)
484 if verbose: print("number of coordinates conrod =", femConrod.NumberOfCoordinates())
485 if verbose: print("Compute eigenmodes conrod ....")
486
487 if not eigenModesNGsolve:
488 femConrod.ComputeEigenmodes(nModes, excludeRigidBodyModes = excludeRigidBodyModes, useSparseSolver = True)
489 else:
490 femConrod.ComputeEigenmodesNGsolve(bfM, bfK, nModes=nModes, excludeRigidBodyModes=excludeRigidBodyModes)
491
492 totalFEcoordinates+=femConrod.NumberOfCoordinates()
493 if verbose: print("eigen freq. conrod=", femConrod.GetEigenFrequenciesHz()[0:nModes])
494
495 if showStresses:
496 print("ComputePostProcessingModes femConrod ... ")
497 start_time = time.time()
498 femConrod.ComputePostProcessingModesNGsolve(fes, material=mat,
499 outputVariableType=varType)
500 print("--- %s seconds ---" % (time.time() - start_time))
501
502 #%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
503 #import piston mesh into EXUDYN FEMinterface
504 femPiston = FEMinterface()
505 #femPiston.ImportMeshFromNGsolve(meshPiston, density, youngsModulus, poissonsRatio, verbose = False)
506 [bfM, bfK, fes] = femPiston.ImportMeshFromNGsolve(meshPiston, density, youngsModulus, poissonsRatio, verbose = False, meshOrder = meshOrder)
507
508 if verbose: print("number of coordinates piston =", femPiston.NumberOfCoordinates())
509 if verbose: print("Compute eigenmodes piston ....")
510
511 if not eigenModesNGsolve:
512 femPiston.ComputeEigenmodes(nModes, excludeRigidBodyModes = excludeRigidBodyModes, useSparseSolver = True)
513 else:
514 femPiston.ComputeEigenmodesNGsolve(bfM, bfK, nModes=nModes, excludeRigidBodyModes=excludeRigidBodyModes)
515
516 totalFEcoordinates+=femPiston.NumberOfCoordinates()
517 if verbose: print("eigen freq. Piston=", femPiston.GetEigenFrequenciesHz()[0:nModes])
518
519 if showStresses:
520 print("ComputePostProcessingModes femPiston ... ")
521 start_time = time.time()
522 femPiston.ComputePostProcessingModesNGsolve(fes, material=mat,
523 outputVariableType=varType)
524 print("--- %s seconds ---" % (time.time() - start_time))
525
526 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
527 #import multiple conrods and pistons
528
529 conrodsRotList = []
530 conrodsPosList = []
531 pistonsRotList = []
532 pistonsPosList = []
533
534 objFFRFconrodList=[]
535 cmsConrodList=[]
536 objFFRFpistonList=[]
537 cmsPistonList=[]
538 pkList = []
539 pcList = []
540 ppList = []
541 zOffsetList = []
542 for iCrank in range(len(crankConfig)):
543 zOffset = db+dk+db + lTotal*iCrank #left end of conrod, for multiple conrods in a loop
544 zOffsetList.append(zOffset)
545 #compute crank (pK), conrod (pC) and piston position (pP) for any crank angle:
546 phi = crankConfig[iCrank]
547 pK = np.array([lk*np.cos(phi),lk*np.sin(phi),0])
548 alpha=np.arcsin(pK[1]/lc)
549 pC = pK + np.array([0.5*lc*np.cos(alpha),-0.5*lc*np.sin(alpha),0])
550 pP = pK + np.array([lc*np.cos(alpha),-lc*np.sin(alpha),0])
551 pkList.append(pK)
552 pcList.append(pC)
553 ppList.append(pP)
554 #print("pK=",pK)
555 #print("pC=",pC)
556 #print("pP=",pP)
557
558 eulerParametersInit = RotationMatrix2EulerParameters(RotationMatrixZ(-alpha))
559 #pRef = [lk+0.5*lc,0,zOffset+0.5*b1] #0-degree
560 pRef = pC + [0,0,zOffset+0.5*b1]
561 conrodsRotList.append(RotationMatrixZ(-alpha).tolist())
562 conrodsPosList.append(pRef.tolist())
563
564 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
565 #import conrod CMS
566 cmsConrod = ObjectFFRFreducedOrderInterface(femConrod)
567 cmsConrodList.append(cmsConrod)
568 objFFRFconrod = cmsConrod.AddObjectFFRFreducedOrderWithUserFunctions(exu, mbs,
569 positionRef=pRef,
570 eulerParametersRef=eulerParametersInit,
571 initialVelocity=[0,0,0],
572 initialAngularVelocity=[0,0,0*fRotorStart*2*pi],
573 gravity = [0,-0*9.81,0],
574 color=[0.1,0.9,0.1,1.])
575 mbs.SetObjectParameter(objFFRFconrod['oFFRFreducedOrder'],'VshowNodes',False)
576 objFFRFconrodList.append(objFFRFconrod)
577
578 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
579 #import piston CMS
580 cmsPiston = ObjectFFRFreducedOrderInterface(femPiston)
581 cmsPistonList.append(cmsPiston)
582
583 pRefPiston = pP+[0,0,zOffset+0.5*b1]
584
585 pistonsRotList.append(RotationMatrixZ(0).tolist())
586 pistonsPosList.append(pRefPiston.tolist())
587
588 objFFRFpiston = cmsPiston.AddObjectFFRFreducedOrderWithUserFunctions(exu, mbs,
589 positionRef=pRefPiston,
590 eulerParametersRef=eulerParameters0,
591 initialVelocity=[0,0,0], initialAngularVelocity=[0,0,0*fRotorStart*2*pi],
592 gravity = [0,-0*9.81,0],
593 color=[0.1,0.9,0.1,1.])
594 mbs.SetObjectParameter(objFFRFpiston['oFFRFreducedOrder'],'VshowNodes',False)
595 objFFRFpistonList.append(objFFRFpiston)
596
597 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
598 if True: #connect bodies:
599 k = 1e6 #joint stiffness
600 d = k*0.002 #joint damping
601 nMarkerPerPiston = 10 #number of markers per crank/conrod/piston part
602
603 genMarkerPos = [[0,0,-d0],[0,0,lTotal*nPistons]]
604 genMarkerR = [r0,r0]
605 genMarkerFEM = [femCrank,femCrank]
606 genMarkerObject = [objFFRFcrank,objFFRFcrank]
607
608 for iCrank in range(len(crankConfig)):
609 genMarkerPos += [pkList[iCrank]+[0,0,zOffsetList[iCrank]],pkList[iCrank]+[0,0,zOffsetList[iCrank]+b1],
610 [-0.5*lc,0,-0.5*dc],[-0.5*lc,0, 0.5*dc],[0.5*lc,0,-0.5*dc],[0.5*lc,0, 0.5*dc],
611 [0,0,-0.5*dc],[0,0,0.5*dc], [-dpb,0,0],[lp-dpb,0,0]]
612 genMarkerR += [r1,r1,
613 r1,r1,r2,r2,
614 r2,r2,0.5*bp,0.5*bp]
615 genMarkerFEM += [femCrank,femCrank,
616 femConrod,femConrod,femConrod,femConrod,
617 femPiston,femPiston,femPiston,femPiston]
618 genMarkerObject += [objFFRFcrank,objFFRFcrank,
619 objFFRFconrodList[iCrank],objFFRFconrodList[iCrank],objFFRFconrodList[iCrank],objFFRFconrodList[iCrank],
620 objFFRFpistonList[iCrank],objFFRFpistonList[iCrank],objFFRFpistonList[iCrank],objFFRFpistonList[iCrank]]
621
622 markerList = []
623 #generate markers for joints:
624 for i in range(len(genMarkerPos)):
625 p = genMarkerPos[i]
626 nodeList=[]
627 if p[2] != 0:
628 nodeList= genMarkerFEM[i].GetNodesOnCircle(p, [0,0,1], genMarkerR[i])
629 else:
630 nodeList= genMarkerFEM[i].GetNodesOnCircle(p, [1,0,0], genMarkerR[i])
631 #print("nodeList"+str(i)+":", nodeList)
632 lenNodeList = len(nodeList)
633 weights = np.array((1./lenNodeList)*np.ones(lenNodeList))
634
635 markerList += [mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=genMarkerObject[i]['oFFRFreducedOrder'],
636 meshNodeNumbers=np.array(nodeList), #these are the meshNodeNumbers
637 weightingFactors=weights))]
638
639 oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
640
641 mGroundPosLeft = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=genMarkerPos[0]))
642 mGroundPosRight = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=genMarkerPos[1]))
643
644
645 #joints for crankshaft/ground
646 oSJleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mGroundPosLeft, markerList[0]],
647 stiffness=[k,k,k], damping=[d,d,d]))
648 oSJright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mGroundPosRight, markerList[1]],
649 stiffness=[k,k,k], damping=[d,d,d]))
650
651 for iCrank in range(len(crankConfig)):
652 mOff = nMarkerPerPiston*iCrank
653 #joints for crankshaft/conrod:
654 oJointCCleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[markerList[mOff+2], markerList[mOff+4]],
655 stiffness=[k,k,k], damping=[d,d,d]))
656 oJointCCright= mbs.AddObject(CartesianSpringDamper(markerNumbers=[markerList[mOff+3], markerList[mOff+5]],
657 stiffness=[k,k,k], damping=[d,d,d]))
658
659 #joints for conrod/piston:
660 oJointCPleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[markerList[mOff+6], markerList[mOff+8]],
661 stiffness=[k,k,k], damping=[d,d,d]))
662 oJointCPright= mbs.AddObject(CartesianSpringDamper(markerNumbers=[markerList[mOff+7], markerList[mOff+9]],
663 stiffness=[k,k,k], damping=[d,d,d]))
664
665 mGroundPosPiston = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround,
666 localPosition=[ppList[iCrank][0],0,zOffsetList[iCrank]+0.5*b1]))
667 oJointPGleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mGroundPosPiston, markerList[mOff+10]],
668 stiffness=[0,k,k], damping=[0,d,d]))
669 oJointPGright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mGroundPosPiston, markerList[mOff+11]],
670 stiffness=[0,k,k], damping=[0,d,d]))
671
672
673 stopTotal = timeit.default_timer()
674 print("\ntotal elapsed time=", stopTotal-startTotal)
675 mbs.Assemble()
676
677 if saveMesh:
678 #%%%
679 dictMesh = {}
680 [points, triangles, normals] = graphics.NGsolveMesh2PointsAndTrigs( ngMesh=ngMeshCrank)
681 dictMesh['ngMeshCrank'] = {'points':points, 'triangles':triangles, 'normals':normals}
682 [points, triangles, normals] = graphics.NGsolveMesh2PointsAndTrigs( ngMesh=ngMeshConrod)
683 dictMesh['ngMeshConrod'] = {'points':points, 'triangles':triangles, 'normals':normals}
684 [points, triangles, normals] = graphics.NGsolveMesh2PointsAndTrigs( ngMesh=ngMeshPiston)
685 dictMesh['ngMeshPiston'] = {'points':points, 'triangles':triangles, 'normals':normals}
686 dictMesh['conrodsRotList'] = conrodsRotList
687 dictMesh['conrodsPosList'] = conrodsPosList
688 dictMesh['pistonsRotList'] = pistonsRotList
689 dictMesh['pistonsPosList'] = pistonsPosList
690
691 fileName = 'testData/pistonEngineNGmesh.hdf5'
692 SaveDictToHDF5(fileName, dictMesh)
693
694
695 #%%now simulate model in exudyn:
696 #%%+++++++++++++++++++++
697 if True:
698 print("totalFEcoordinates=",totalFEcoordinates)
699
700 simulationSettings = exu.SimulationSettings()
701
702 nodeDrawSize = 0.0005
703 SC.visualizationSettings.general.textSize = 14 #30 for cover figure
704 SC.visualizationSettings.general.useGradientBackground = True
705 SC.visualizationSettings.openGL.lineWidth = 2
706
707 SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
708 SC.visualizationSettings.nodes.drawNodesAsPoint = False
709 SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
710 SC.visualizationSettings.connectors.show = False
711
712 SC.visualizationSettings.nodes.show = False
713 SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
714 SC.visualizationSettings.nodes.basisSize = 0.12
715 SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
716
717 SC.visualizationSettings.openGL.showFaceEdges = True
718 SC.visualizationSettings.openGL.showFaces = True
719 SC.visualizationSettings.openGL.multiSampling = 4
720
721 SC.visualizationSettings.sensors.show = True
722 SC.visualizationSettings.sensors.drawSimplified = False
723 SC.visualizationSettings.sensors.defaultSize = 0.01
724 SC.visualizationSettings.markers.drawSimplified = False
725 SC.visualizationSettings.markers.show = False
726 SC.visualizationSettings.markers.defaultSize = 0.01
727
728 SC.visualizationSettings.loads.drawSimplified = False
729
730 #SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
731 SC.visualizationSettings.contour.outputVariableComponent = -1
732 SC.visualizationSettings.contour.reduceRange = True
733 #SC.visualizationSettings.contour.automaticRange = False
734 #SC.visualizationSettings.contour.maxValue = 3e7
735 # SC.visualizationSettings.contour.minValue = -0.0003
736 # SC.visualizationSettings.contour.maxValue = 0.0003
737
738 simulationSettings.solutionSettings.solutionInformation = "NGsolve/NETGEN engine test"
739
740 h=0.05e-3
741 tEnd = 2
742
743 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
744 simulationSettings.timeIntegration.endTime = tEnd
745 simulationSettings.solutionSettings.solutionWritePeriod = h*10 #writing already costs much time
746 simulationSettings.timeIntegration.verboseMode = 1
747 #simulationSettings.timeIntegration.verboseModeFile = 3
748 simulationSettings.timeIntegration.newton.useModifiedNewton = True
749
750 simulationSettings.solutionSettings.sensorsWritePeriod = h
751 #simulationSettings.solutionSettings.coordinatesSolutionFileName = "solution/coordinatesSolution.txt"
752 simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse #faster, because system size already quite large
753
754 simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5 #SHOULD work with 0.9 as well
755 simulationSettings.displayStatistics = True
756 #simulationSettings.displayComputationTime = True
757 SC.visualizationSettings.general.autoFitScene = False #for reloading of renderState to work
758
759 #create animation:
760 if False:
761 simulationSettings.solutionSettings.recordImagesInterval = 0.001
762 SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
763 SC.visualizationSettings.window.renderWindowSize=[1920,1080]
764
765 SC.renderer.Start()
766 if 'renderState' in exu.sys: SC.renderer.SetState(exu.sys['renderState']) #load last model view
767
768 SC.renderer.DoIdleTasks() #press space to continue
769
770 simulate = True #set false to show last stored solution
771 if simulate:
772 mbs.SolveDynamic(simulationSettings)
773 else:
774 SC.visualizationSettings.general.autoFitScene = False
775 sol = LoadSolutionFile('coordinatesSolution.txt')
776 if False: #directly show animation
777 AnimateSolution(mbs, solution=sol, rowIncrement = 1, timeout=0.01,
778 createImages = False, runLoop = True)
779 else: #interact with animation
780
781 mbs.SolutionViewer(sol, rowIncrement=1, timeout=0.02)
782
783
784 if False: #draw with matplotlib, export as pdf
785 SC.visualizationSettings.exportImages.saveImageFormat = "TXT"
786 SC.visualizationSettings.exportImages.saveImageAsTextTriangles=True
787 SC.renderer.RedrawAndSaveImage() #uses default filename
788
789 from exudyn.plot import LoadImage, PlotImage
790
791 # plot 2D
792 # data = LoadImage('images/frame00000.txt', trianglesAsLines=True)
793 # PlotImage(data, HT=HomogeneousTransformation(RotationMatrixZ(0.5*pi)@RotationMatrixX(0.5*pi), [0,0,0]),
794 # lineWidths=0.5, lineStyles='-', title='', closeAll=True, plot3D=False,
795 # fileName='images/test.pdf')
796
797 data = LoadImage('images/frame00000.txt', trianglesAsLines=False)
798 PlotImage(data, HT=HomogeneousTransformation(2.5*RotationMatrixZ(0.5*pi)@RotationMatrixY(-0.5*pi), [0,1,0.25]),
799 lineWidths=0.5, lineStyles='-', triangleEdgeColors='black', triangleEdgeWidths=0.25, title='', closeAll=True, plot3D=True,
800 fileName='images/test3D.pdf')
801
802 SC.renderer.DoIdleTasks()
803 SC.renderer.Stop() #safely close rendering window!
804 lastRenderState = SC.renderer.GetState() #store model view for next simulation