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