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