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