NGsolveCMStutorial.py
You can view and download this file on Github: NGsolveCMStutorial.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test for Hurty-Craig-Bampton modes using a simple flexible pendulum meshed with Netgen
5#
6# Author: Johannes Gerstmayr
7# Date: 2021-04-20
8#
9# 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.
10#
11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12
13
14import exudyn as exu
15from exudyn.itemInterface import *
16from exudyn.utilities import *
17from exudyn.FEM import *
18from exudyn.graphicsDataUtilities import *
19
20SC = exu.SystemContainer()
21mbs = SC.AddSystem()
22
23import numpy as np
24import time
25
26#import timeit
27
28import exudyn.basicUtilities as eb
29import exudyn.rigidBodyUtilities as rb
30import exudyn.utilities as eu
31
32
33useGraphics = True
34fileName = 'testData/netgenHinge' #for load/save of FEM data
35
36#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
37#netgen/meshing part:
38fem = FEMinterface()
39
40#geometrical parameters:
41L = 0.4 #Length of plate (X)
42w = 0.04 #width of plate (Y)
43h = 0.02 #height of plate (Z)
44d = 0.03 #diameter of bolt
45D = d*2 #diameter of bushing
46b = 0.05 #length of bolt
47nModes = 8
48meshH = 0.01 #0.01 is default, 0.002 gives 100000 nodes and is fairly converged;
49#meshH = 0.0014 #203443 nodes, takes 1540 seconds for eigenmode computation (free-free) and 753 seconds for postprocessing on i9
50
51#steel:
52rho = 7850
53Emodulus=2.1e11
54nu=0.3
55
56#test high flexibility
57Emodulus=2e8
58# nModes = 32
59
60
61#helper function for cylinder with netgen
62def CSGcylinder(p0,p1,r):
63 v = VSub(p1,p0)
64 v = Normalize(v)
65 cyl = Cylinder(Pnt(p0[0],p0[1],p0[2]), Pnt(p1[0],p1[1],p1[2]),
66 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]))
67 return cyl
68
69meshCreated = False
70
71#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
72if True: #needs netgen/ngsolve to be installed to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
73 import ngsolve as ngs
74 import netgen
75 from netgen.meshing import *
76
77 from netgen.geom2d import unit_square
78 #import netgen.libngpy as libng
79 from netgen.csg import *
80
81 geo = CSGeometry()
82
83 #plate
84 block = OrthoBrick(Pnt(0, 0, -0.5*h),Pnt(L, w, 0.5*h))
85
86 #bolt
87 bolt0 = CSGcylinder(p0=[0,w,0], p1=[0,0,0], r=1.6*h)
88 bolt = CSGcylinder(p0=[0,0.5*w,0], p1=[0,-b,0], r=0.5*d)
89
90 #bushing
91 bushing = (CSGcylinder(p0=[L,w,0], p1=[L,-b,0], r=0.5*D) -
92 CSGcylinder(p0=[L,0,0], p1=[L,-b*1.1,0], r=0.5*d))
93
94 geo.Add(block+bolt0+bolt+bushing)
95
96 curvaturesafety = 5
97 if meshH==0.04:
98 curvaturesafety = 1.2#this case is for creating very small files ...
99
100 mesh = ngs.Mesh( geo.GenerateMesh(maxh=meshH, curvaturesafety=curvaturesafety))
101 mesh.Curve(1)
102
103 if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
104 # import netgen
105 import netgen.gui
106 ngs.Draw(mesh)
107 for i in range(10000000):
108 netgen.Redraw() #this makes the netgen window interactive
109 time.sleep(0.05)
110
111 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
112 #Use fem to import FEM model and create FFRFreducedOrder object
113 [bfM, bfK, fes] = fem.ImportMeshFromNGsolve(mesh, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
114 meshCreated = True
115 if (meshH==0.04):
116 print('save file')
117 fem.SaveToFile(fileName)
118
119
120#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
121#compute Hurty-Craig-Bampton modes
122if True: #now import mesh as mechanical model to EXUDYN
123 if not meshCreated: fem.LoadFromFile(fileName)
124
125 boltP1=[0,0,0]
126 boltP2=[0,-b,0]
127 nodesOnBolt = fem.GetNodesOnCylinder(boltP1, boltP2, radius=0.5*d)
128 #print("boundary nodes bolt=", nodesOnBolt)
129 nodesOnBoltLen = len(nodesOnBolt)
130 nodesOnBoltWeights = np.array((1./nodesOnBoltLen)*np.ones(nodesOnBoltLen))
131
132 bushingP1=[L,0,0]
133 bushingP2=[L,-b,0]
134 nodesOnBushing = fem.GetNodesOnCylinder(bushingP1, bushingP2, radius=0.5*d)
135 #print("boundary nodes bushing=", nodesOnBushing)
136 nodesOnBushingLen = len(nodesOnBushing)
137 nodesOnBushingWeights = np.array((1./nodesOnBushingLen)*np.ones(nodesOnBushingLen))
138
139 print("nNodes=",fem.NumberOfNodes())
140
141 strMode = ''
142 if True: #pure eigenmodes
143 print("compute eigen modes... ")
144 start_time = time.time()
145
146 if False: #faster but not so accurate
147 fem.ComputeEigenmodesNGsolve(bfM, bfK, nModes, excludeRigidBodyModes = 6)
148 else:
149 fem.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
150 print("eigen modes computation needed %.3f seconds" % (time.time() - start_time))
151 print("eigen freq.=", fem.GetEigenFrequenciesHz())
152
153 else:
154 strMode = 'HCB'
155 #boundaryList = [nodesOnBolt, nodesOnBolt, nodesOnBushing] #for visualization, use first interface twice
156 boundaryList = [nodesOnBolt, nodesOnBushing]
157
158 print("compute HCB modes... ")
159 start_time = time.time()
160 fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList,
161 nEigenModes=nModes,
162 useSparseSolver=True,
163 computationMode = HCBstaticModeSelection.RBE2)
164
165 print("eigen freq.=", fem.GetEigenFrequenciesHz())
166 print("HCB modes needed %.3f seconds" % (time.time() - start_time))
167
168
169
170 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
171 #compute stress modes for postprocessing (inaccurate for coarse meshes, just for visualization):
172 if True:
173 mat = KirchhoffMaterial(Emodulus, nu, rho)
174 varType = exu.OutputVariableType.StressLocal
175 #varType = exu.OutputVariableType.StrainLocal
176 print("ComputePostProcessingModes ... (may take a while)")
177 start_time = time.time()
178 #without NGsolve:
179 if True: #faster with ngsolve
180 fem.ComputePostProcessingModesNGsolve(fes, material=mat,
181 outputVariableType=varType)
182 else:
183 fem.ComputePostProcessingModes(material=mat,
184 outputVariableType=varType)
185 print(" ... needed %.3f seconds" % (time.time() - start_time))
186 SC.visualizationSettings.contour.reduceRange=True
187 SC.visualizationSettings.contour.outputVariable = varType
188 SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
189 else:
190 varType = exu.OutputVariableType.DisplacementLocal
191 SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
192 SC.visualizationSettings.contour.outputVariableComponent = 0
193
194 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
195 print("create CMS element ...")
196 cms = ObjectFFRFreducedOrderInterface(fem)
197
198 objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
199 initialVelocity=[0,0,0],
200 initialAngularVelocity=[0,0,0],
201 color=[0.9,0.9,0.9,1.],
202 )
203
204 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
205 #add markers and joints
206 nodeDrawSize = 0.0025 #for joint drawing
207
208
209 #mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
210
211 if True:
212 boltMidPoint = 0.5*(np.array(boltP1)+boltP2)
213
214 oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
215
216 altApproach = True
217 mBolt = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
218 meshNodeNumbers=np.array(nodesOnBolt), #these are the meshNodeNumbers
219 #referencePosition=boltMidPoint,
220 useAlternativeApproach=altApproach,
221 weightingFactors=nodesOnBoltWeights))
222 bushingMidPoint = 0.5*(np.array(bushingP1)+bushingP2)
223
224 #add marker for visualization of boundary nodes
225 mBushing = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
226 meshNodeNumbers=np.array(nodesOnBushing), #these are the meshNodeNumbers
227 #referencePosition=bushingMidPoint,
228 useAlternativeApproach=altApproach,
229 weightingFactors=nodesOnBushingWeights))
230
231 lockedAxes=[1,1,1,1,1*0,1]
232 if True:
233
234 mGroundBolt = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround,
235 localPosition=boltMidPoint,
236 visualization=VMarkerBodyRigid(show=True)))
237 mbs.AddObject(GenericJoint(markerNumbers=[mGroundBolt, mBolt],
238 constrainedAxes = lockedAxes,
239 visualization=VGenericJoint(show=False, axesRadius=0.1*b, axesLength=0.1*b)))
240
241 else:
242
243 mGroundBushing = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=bushingMidPoint))
244 mbs.AddObject(GenericJoint(markerNumbers=[mGroundBushing, mBushing],
245 constrainedAxes = lockedAxes,
246 visualization=VGenericJoint(axesRadius=0.1*b, axesLength=0.1*b)))
247
248
249 if False:
250 cms = ObjectFFRFreducedOrderInterface(fem)
251
252 objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
253 initialVelocity=[0,0,0],
254 initialAngularVelocity=[0,0,0],
255 color=[0.9,0.9,0.9,1.],
256 )
257
258 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
259 #animate modes
260 SC.visualizationSettings.markers.show = True
261 SC.visualizationSettings.markers.defaultSize=0.0075
262 SC.visualizationSettings.markers.drawSimplified = False
263
264 SC.visualizationSettings.loads.show = False
265 SC.visualizationSettings.loads.drawSimplified = False
266 SC.visualizationSettings.loads.defaultSize=0.1
267 SC.visualizationSettings.loads.defaultRadius = 0.002
268
269 SC.visualizationSettings.openGL.multiSampling=4
270 SC.visualizationSettings.openGL.lineWidth=2
271
272 if False: #activate to animate modes
273 from exudyn.interactive import AnimateModes
274 mbs.Assemble()
275 SC.visualizationSettings.nodes.show = False
276 SC.visualizationSettings.openGL.showFaceEdges = True
277 SC.visualizationSettings.openGL.multiSampling=4
278 SC.visualizationSettings.openGL.lineWidth=2
279 SC.visualizationSettings.window.renderWindowSize = [1600,1080]
280 SC.visualizationSettings.contour.showColorBar = False
281 SC.visualizationSettings.general.textSize = 16
282
283 #%%+++++++++++++++++++++++++++++++++++++++
284 #animate modes of ObjectFFRFreducedOrder (only needs generic node containing modal coordinates)
285 SC.visualizationSettings.general.autoFitScene = False #otherwise, model may be difficult to be moved
286
287 nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
288 AnimateModes(SC, mbs, nodeNumber, period=0.1, showTime=False, renderWindowText='Hurty-Craig-Bampton: 2 x 6 static modes and 8 eigenmodes\n',
289 runOnStart=True)
290 # import sys
291 # sys.exit()
292
293 #add gravity (not necessary if user functions used)
294 oFFRF = objFFRF['oFFRFreducedOrder']
295 mBody = mbs.AddMarker(MarkerBodyMass(bodyNumber=oFFRF))
296 mbs.AddLoad(LoadMassProportional(markerNumber=mBody, loadVector= [0,0,-9.81]))
297
298
299 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
300 fileDir = 'solution/'
301 # sensBolt = mbs.AddSensor(SensorMarker(markerNumber=mBolt,
302 # fileName=fileDir+'hingePartBoltPos'+str(nModes)+strMode+'.txt',
303 # outputVariableType = exu.OutputVariableType.Position))
304 # sensBushing= mbs.AddSensor(SensorMarker(markerNumber=mBushing,
305 # fileName=fileDir+'hingePartBushingPos'+str(nModes)+strMode+'.txt',
306 # outputVariableType = exu.OutputVariableType.Position))
307 sensBushingVel= mbs.AddSensor(SensorMarker(markerNumber=mBushing,
308 fileName=fileDir+'hingePartBushingVel'+str(nModes)+strMode+'.txt',
309 outputVariableType = exu.OutputVariableType.Velocity))
310 sensBushing= mbs.AddSensor(SensorMarker(markerNumber=mBushing,
311 fileName=fileDir+'hingePartBushing'+str(nModes)+strMode+'.txt',
312 outputVariableType = exu.OutputVariableType.Position))
313
314 mbs.Assemble()
315
316 simulationSettings = exu.SimulationSettings()
317
318 SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
319 SC.visualizationSettings.nodes.drawNodesAsPoint = False
320 SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
321
322 SC.visualizationSettings.nodes.show = False
323 SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
324 SC.visualizationSettings.nodes.basisSize = 0.12
325 SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
326
327 SC.visualizationSettings.openGL.showFaceEdges = True
328 SC.visualizationSettings.openGL.showFaces = True
329
330 SC.visualizationSettings.sensors.show = True
331 SC.visualizationSettings.sensors.drawSimplified = False
332 SC.visualizationSettings.sensors.defaultSize = 0.01
333
334
335 simulationSettings.solutionSettings.solutionInformation = "CMStutorial "+str(nModes)+" "+strMode+"modes"
336
337 h=0.25e-3*4
338 tEnd = 0.25*8
339
340 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
341 simulationSettings.timeIntegration.endTime = tEnd
342 simulationSettings.solutionSettings.writeSolutionToFile = True
343 simulationSettings.timeIntegration.verboseMode = 1
344 #simulationSettings.timeIntegration.verboseModeFile = 3
345 simulationSettings.timeIntegration.newton.useModifiedNewton = True
346
347 simulationSettings.solutionSettings.sensorsWritePeriod = h
348
349 simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
350 #simulationSettings.displayStatistics = True
351 simulationSettings.displayComputationTime = True
352
353 #create animation:
354 # simulationSettings.solutionSettings.recordImagesInterval = 0.005
355 # SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
356 SC.visualizationSettings.window.renderWindowSize=[1920,1080]
357 SC.visualizationSettings.openGL.multiSampling = 4
358
359 useGraphics=True
360 if True:
361 if useGraphics:
362 SC.visualizationSettings.general.autoFitScene=False
363
364 exu.StartRenderer()
365 if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view
366
367 mbs.WaitForUserToContinue() #press space to continue
368
369 #SC.RedrawAndSaveImage()
370 if True:
371 # mbs.SolveDynamic(solverType=exu.DynamicSolverType.TrapezoidalIndex2,
372 # simulationSettings=simulationSettings)
373 mbs.SolveDynamic(simulationSettings=simulationSettings)
374 else:
375 mbs.SolveStatic(simulationSettings=simulationSettings)
376
377 # uTip = mbs.GetSensorValues(sensTipDispl)[1]
378 # print("nModes=", nModes, ", tip displacement=", uTip)
379
380 if varType == exu.OutputVariableType.StressLocal:
381 mises = CMSObjectComputeNorm(mbs, 0, exu.OutputVariableType.StressLocal, 'Mises')
382 print('max von-Mises stress=',mises)
383
384 if useGraphics:
385 SC.WaitForRenderEngineStopFlag()
386 exu.StopRenderer() #safely close rendering window!
387
388 if False:
389
390 mbs.PlotSensor(sensorNumbers=[sensBushingVel], components=[1])
391
392#%%
393if False:
394 import matplotlib.pyplot as plt
395 import matplotlib.ticker as ticker
396 import exudyn as exu
397 from exudyn.utilities import *
398 CC = PlotLineCode
399 comp = 1 #1=x, 2=y, ...
400 var = ''
401 # data = np.loadtxt('solution/hingePartBushing'+var+'2.txt', comments='#', delimiter=',')
402 # plt.plot(data[:,0], data[:,comp], CC(7), label='2 eigenmodes')
403 # data = np.loadtxt('solution/hingePartBushing'+var+'4.txt', comments='#', delimiter=',')
404 # plt.plot(data[:,0], data[:,comp], CC(8), label='4 eigenmodes')
405 data = np.loadtxt('solution/hingePartBushing'+var+'8.txt', comments='#', delimiter=',')
406 plt.plot(data[:,0], data[:,comp], CC(9), label='8 eigenmodes')
407 data = np.loadtxt('solution/hingePartBushing'+var+'16.txt', comments='#', delimiter=',')
408 plt.plot(data[:,0], data[:,comp], CC(10), label='16 eigenmodes')
409 data = np.loadtxt('solution/hingePartBushing'+var+'32.txt', comments='#', delimiter=',')
410 plt.plot(data[:,0], data[:,comp], CC(11), label='32 eigenmodes')
411
412 data = np.loadtxt('solution/hingePartBushing'+var+'2HCB.txt', comments='#', delimiter=',')
413 plt.plot(data[:,0], data[:,comp], CC(1), label='HCB + 2 eigenmodes')
414 data = np.loadtxt('solution/hingePartBushing'+var+'4HCB.txt', comments='#', delimiter=',')
415 plt.plot(data[:,0], data[:,comp], CC(2), label='HCB + 4 eigenmodes')
416 data = np.loadtxt('solution/hingePartBushing'+var+'8HCB.txt', comments='#', delimiter=',')
417 plt.plot(data[:,0], data[:,comp], CC(3), label='HCB + 8 eigenmodes')
418 data = np.loadtxt('solution/hingePartBushing'+var+'16HCB.txt', comments='#', delimiter=',')
419 plt.plot(data[:,0], data[:,comp], CC(4), label='HCB + 16 eigenmodes')
420 data = np.loadtxt('solution/hingePartBushing'+var+'32HCB.txt', comments='#', delimiter=',')
421 plt.plot(data[:,0], data[:,comp], CC(5), label='HCB + 32 eigenmodes')
422 data = np.loadtxt('solution/hingePartBushing'+var+'64HCB.txt', comments='#', delimiter=',')
423 plt.plot(data[:,0], data[:,comp], CC(6), label='HCB + 64 eigenmodes')
424 data = np.loadtxt('solution/hingePartBushing'+var+'128HCB.txt', comments='#', delimiter=',')
425 plt.plot(data[:,0], data[:,comp], CC(7), label='HCB + 128 eigenmodes')
426
427
428 ax=plt.gca() # get current axes
429 ax.grid(True, 'major', 'both')
430 ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
431 ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
432 #
433 plt.xlabel("time (s)")
434 plt.ylabel("y-component of tip velocity of hinge (m)")
435 plt.legend() #show labels as legend
436 plt.tight_layout()
437 plt.show()