NGsolveModalAnalysis.py
You can view and download this file on Github: NGsolveModalAnalysis.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Example for linear modal analysis with Hurty-Craig-Bampton modes
5#
6# Author: Johannes Gerstmayr
7# Date: 2024-10-30
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.utilities import * #includes itemInterface and rigidBodyUtilities
16import exudyn.graphics as graphics
17from exudyn.FEM import *
18import sys
19
20SC = exu.SystemContainer()
21mbs = SC.AddSystem()
22
23import numpy as np
24import time
25
26import numpy as np
27
28useGraphics = True
29fileName = OutputFilePath('solution/modalAnalysisFEM', 'NGsolveModalAnalysis') #generated output, not input: goes to the ignored solution/ (#2491)
30#(the path no longer needs '../Examples' to survive a run out of TestModels)
31
32
33
34#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
35#netgen/meshing part:
36fem = FEMinterface()
37
38
39#Create mesh for crane-like structure on 4 piles
40hFoot = 1 #z-dir
41wFoot = 0.4
42
43wxBody = 2.5
44wyBody = 6
45wzBody = 0.3
46
47wArm = 0.5
48wHole = 0.4
49lArm = 16
50angleArm = 60/180*np.pi
51
52nModes = 8
53meshSize = wArm*0.5
54
55rho = 7800
56Emodulus=2.1e11
57nu=0.3
58meshOrder = 1 #use order 2 for higher accuracy, but more unknowns
59
60meshCreated = False
61
62pFoot0 = [-0.5*wxBody-wFoot,-0.5*wyBody,0]
63pFoot1 = [0.5*wxBody,-0.5*wyBody,0]
64pFoot2 = [-0.5*wxBody-wFoot,0.5*wyBody-wFoot,0]
65pFoot3 = [0.5*wxBody,0.5*wyBody-wFoot,0]
66pTip = [0,lArm*np.cos(angleArm),lArm*np.sin(angleArm)+hFoot]
67dirTip = [0,np.cos(angleArm),np.sin(angleArm)]
68
69#list of points and directions for finding boundary nodes:
70dirZ = [0,0,1]
71pList = [pFoot0,pFoot1,pFoot2,pFoot3,pTip]
72dirList = [dirZ,dirZ,dirZ,dirZ,dirTip]
73
74#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
75if True: #needs netgen/ngsolve to be installed with pip install; to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
76
77 import ngsolve as ngs
78 from netgen.meshing import *
79 from netgen.csg import *
80
81 geo = CSGeometry()
82
83 #for i in range(2):
84
85 block1 = OrthoBrick(Pnt(*pFoot0),
86 Pnt(-0.5*wxBody,-0.5*wyBody+wFoot,hFoot))
87 block2 = OrthoBrick(Pnt(*pFoot1),
88 Pnt(0.5*wxBody+wFoot,-0.5*wyBody+wFoot,hFoot))
89 block3 = OrthoBrick(Pnt(*pFoot2),
90 Pnt(-0.5*wxBody,0.5*wyBody,hFoot))
91 block4 = OrthoBrick(Pnt(*pFoot3),
92 Pnt(0.5*wxBody+wFoot,0.5*wyBody,hFoot))
93
94 body = OrthoBrick( Pnt(-0.5*wxBody,-0.5*wyBody,hFoot-wzBody),Pnt(0.5*wxBody,0.5*wyBody,hFoot))
95
96
97 arm1 = Cylinder(Pnt(0,0,hFoot),
98 Pnt(*pTip),
99 wArm*0.5)
100 arm1in = Cylinder(Pnt(0,0,hFoot),
101 Pnt(0,lArm*np.cos(angleArm),lArm*np.sin(angleArm)+hFoot),
102 wHole*0.5)
103 arm2 = Plane(Pnt(0,0,hFoot*0.99),Vec(0,0,-1)) #Half-space plane points to exterior of included space
104 arm3 = Plane(Pnt(0,lArm*np.cos(angleArm),lArm*np.sin(angleArm)+hFoot),
105 Vec(*dirTip) )
106
107 geo.Add(block1+block2+
108 block3+block4+
109 body+(arm1-arm1in)*arm2*arm3)
110
111 mesh = ngs.Mesh( geo.GenerateMesh(maxh=meshSize))
112 mesh.Curve(1)
113
114 #%%
115 if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
116 # import netgen
117 import netgen.gui
118 ngs.Draw (mesh)
119 for i in range(200):
120 netgen.Redraw() #this makes the window interactive
121 time.sleep(0.05)
122 sys.exit()
123 meshCreated = True
124 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
125 #Use fem to import FEM model and create FFRFreducedOrder object
126 [bfM, bfK, fes] = fem.ImportMeshFromNGsolve(mesh, density=rho,
127 youngsModulus=Emodulus, poissonsRatio=nu,
128 meshOrder=meshOrder)
129 fem.SaveToFile(fileName)
130
131#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
132#compute Hurty-Craig-Bampton modes
133if not meshCreated: #now import mesh as mechanical model to EXUDYN
134 fem.LoadFromFile(fileName)
135
136if True:
137 boundaryNodesList = []
138 boundaryWeightsList = []
139 for i in range(len(pList)):
140 if i < len(pList)-1:
141 nodesBoundary = fem.GetNodesInCube(np.array(pList[i])-[wFoot,wFoot,1e-3],
142 np.array(pList[i])+[wFoot,wFoot,1e-3])
143 else:
144 #take care: GetNodesInPlane takes the whole set of nodes in the infinite plane!
145 nodesBoundary = fem.GetNodesInPlane(pList[i], dirList[i])
146 print('nodes B'+str(i),':',nodesBoundary)
147 weightsBoundary = fem.GetNodeWeightsFromSurfaceAreas(nodesBoundary)
148 boundaryNodesList.append(nodesBoundary)
149 boundaryWeightsList.append(weightsBoundary)
150
151
152 print("total number of nodes=",fem.NumberOfNodes())
153
154 print("compute HCB modes... ")
155 start_time = time.time()
156 fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryNodesList, #use boundaryNodesList[:-1] to exclude boundary of arm tip surface!
157 nEigenModes=nModes,
158 useSparseSolver=True,
159 excludeRigidBodyMotion=False, #for modal analysis, we must set this to False
160 computationMode = HCBstaticModeSelection.RBE2)
161 print("HCB modes needed %.3f seconds" % (time.time() - start_time))
162
163 print('eigen frequencies:',np.round(fem.GetEigenFrequenciesHz(),4))
164 #==> if tip is not fixed, gives: 1.8202 1.8223 11.3098 11.3469 31.2923 31.3375 46.6563 53.1684 Hz
165
166 #alternatives:
167 #fem.ComputeEigenModesWithBoundaryNodes(boundaryNodes=nodesLeftPlane, nEigenModes=nModes, useSparseSolver=False)
168 #fem.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
169 #print("eigen freq.=", fem.GetEigenFrequenciesHz())
170
171
172 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
173 #compute stress modes for postprocessing (inaccurate for coarse meshes, just for visualization):
174 if True and meshCreated: #if True, use meshOrder=2 for second order elements!
175 mat = KirchhoffMaterial(Emodulus, nu, rho)
176 varType = exu.OutputVariableType.StressLocal
177 #varType = exu.OutputVariableType.StrainLocal
178 print("ComputePostProcessingModes ... (may take a while)")
179 start_time = time.time()
180 if True: #faster with ngsolve; requires fes
181 fem.ComputePostProcessingModesNGsolve(fes, material=mat,
182 outputVariableType=varType)
183 else:
184 fem.ComputePostProcessingModes(material=mat,
185 outputVariableType=varType)
186 print(" ... needed %.3f seconds" % (time.time() - start_time))
187 SC.visualizationSettings.contour.reduceRange=False
188 SC.visualizationSettings.contour.outputVariable = varType
189 SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
190 else:
191 SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
192 SC.visualizationSettings.contour.outputVariableComponent = 1
193
194 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
195 #Now create a simulation model with the modal body - Component Mode Synthesis (CMS)-element
196 #for this case, we use the AddObjectFFRFreducedOrder, but lock all rigid body DOF of this object
197 print("create CMS element ...")
198 cms = ObjectFFRFreducedOrderInterface(fem,
199 rigidBodyNodeType=exu.NodeType.RotationRxyz, #use this node to allow fixing all rigid body coordinates
200 )
201
202 alphaDamping = 0.01 #set to 0 to eliminate damping (for higher modes, this should be smaller!)
203 objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
204 initialVelocity=[0,0,0],
205 initialAngularVelocity=[0,0,0],
206 gravity=[0,0,-9.81], #applied to all nodes of the body
207 color=[0.1,0.9,0.1,1.],
208 stiffnessProportionalDamping=alphaDamping, #alpha*K
209 )
210
211 #fix rigid body motion of object (otherwise, it could also move!)
212 nodeFFRF = objFFRF['nRigidBody']
213 nCoords = len(mbs.GetNode(nodeFFRF)['referenceCoordinates'])
214 nGround = mbs.AddNode(NodePointGround(visualization=VNodePointGround(show=False)))
215 mnGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
216 for i in range(nCoords):
217 mCoord = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nodeFFRF, coordinate=i))
218 mbs.AddObject(ObjectConnectorCoordinate(markerNumbers=[mnGround, mCoord]))
219
220
221 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
222 #animate modes (only for visualization! => set False for computation!)
223 if False:
224 from exudyn.interactive import AnimateModes
225 mbs.Assemble()
226 SC.visualizationSettings.nodes.show = False
227 SC.visualizationSettings.view0.scene.showFaceEdges = True
228 SC.visualizationSettings.openGL.multiSampling=4
229 SC.visualizationSettings.view0.window.renderWindowSize = [1920,1200]
230 # SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
231 # SC.visualizationSettings.contour.outputVariableComponent = 0 #component
232
233
234 #%%+++++++++++++++++++++++++++++++++++++++
235 #animate modes of ObjectFFRFreducedOrder (only needs generic node containing modal coordinates)
236 SC.visualizationSettings.general.autoFitScene = False #otherwise, model may be difficult to be moved
237
238 nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
239 AnimateModes(SC, mbs, nodeNumber, scaleAmplitude=20)
240 import sys
241 sys.exit()
242
243
244 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
245 #fix 4 feet:
246 #add ground object:
247 oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
248
249 for i in range(len(pList)-1):
250 p = np.array(pList[i]) + [0.5*wFoot,0.5*wFoot,0] #mid point of foot!
251 mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=p))
252
253 mFoot = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
254 meshNodeNumbers=np.array(boundaryNodesList[i]), #these are the meshNodeNumbers
255 weightingFactors=boundaryWeightsList[i]))
256 #add joint on foot i; rotation is left free
257 mbs.AddObject(GenericJoint(markerNumbers=[mGround, mFoot],
258 constrainedAxes = [1,1,1, 0,0,0], #fix displacements x,y,z; leave rotations x,y,z free!
259 visualization=VGenericJoint(axesRadius=0.1*wFoot, axesLength=0.1*wFoot)))
260
261
262 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
263 #add arm tip force:
264 mTip = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
265 meshNodeNumbers=np.array(boundaryNodesList[-1]),
266 weightingFactors=boundaryWeightsList[-1]))
267 #user function to add arbitrary loads:
268 def UFload(mbs,t, loadVector):
269 if t <= 10: #activate load after 10 seconds (wait for decay of initial disturbances)
270 return [0,0,0]
271 else: #dynamic load:
272 #for linear mesh, meshSize=wArm*0.5, resonances are approx:
273 # 1.82, 11.3, 31.3, 46.6, 53.2 Hz
274 # => for larger f, you should increase deformationScaleFactor below to 10 or larger to visualize vibrations
275 f = 1.8
276 return np.sin(t*2*pi*f) * np.array(loadVector)
277 #alternatively: add load as step:
278 # return np.array(loadVector)
279
280 mbs.AddLoad(LoadForceVector(markerNumber=mTip,
281 loadVector=[10000,0,-20000], #some large load to see vibrations
282 loadVectorUserFunction=UFload) )
283
284 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
285 #add sensor for tip motion:
286 sensTipDispl = mbs.AddSensor(SensorMarker(markerNumber=mTip,
287 fileName='solution/armTipDisplacement',
288 storeInternal=True,
289 outputVariableType = exu.OutputVariableType.Displacement))
290
291 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
292 mbs.Assemble()
293
294
295 #+++++++++++++++++++++++++++++++++++++
296 #some simulation and visualization settings
297 simulationSettings = exu.SimulationSettings()
298
299 simulationSettings.solution.file.information = "Modal analysis example"
300
301 stepSize=5e-3 #step size can be large, because system is linear!
302 tEnd = 20
303
304 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
305 simulationSettings.timeIntegration.endTime = tEnd
306 simulationSettings.timeIntegration.newton.useModifiedNewton = True #faster simulation
307 simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
308
309 # simulationSettings.solution.file.write = False
310 simulationSettings.timeIntegration.verboseMode = 1
311 simulationSettings.timeIntegration.newton.useModifiedNewton = True
312
313 simulationSettings.solution.sensors.writePeriod = stepSize
314 simulationSettings.solution.file.writePeriod = stepSize*2
315 #simulationSettings.show.computationTime = True
316
317 #++++++++++++++++
318 SC.visualizationSettings.nodes.defaultSize = meshSize*0.05
319 SC.visualizationSettings.nodes.drawNodesAsPoint = False
320 SC.visualizationSettings.connectors.defaultSize = meshSize*0.05
321
322 SC.visualizationSettings.nodes.show = False
323 SC.visualizationSettings.nodes.showBasis = False #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.view0.scene.showFaceEdges = True
328 SC.visualizationSettings.view0.scene.showFaces = True
329
330 SC.visualizationSettings.sensors.show = True
331 SC.visualizationSettings.sensors.drawSimplified = False
332 SC.visualizationSettings.sensors.defaultSize = 0.01
333 SC.visualizationSettings.markers.drawSimplified = False
334 SC.visualizationSettings.markers.show = False
335 SC.visualizationSettings.markers.defaultSize = 0.01
336
337 SC.visualizationSettings.loads.drawSimplified = False
338
339 SC.visualizationSettings.view0.window.renderWindowSize=[1920,1080]
340 SC.visualizationSettings.openGL.multiSampling = 4 #set to 1 for less powerful graphics cards!
341
342 useGraphics=True
343 if True:
344 if useGraphics:
345 SC.visualizationSettings.general.autoFitScene=False
346
347 SC.renderer.Start()
348 SC.renderer.RestoreSavedState() #load last model view
349
350 # SC.renderer.DoIdleTasks() #press space to continue
351
352 #we could also perform a static analysis at the beginning,
353 # then starting dynamic problem from static equilibrium
354 #mbs.SolveStatic(simulationSettings=simulationSettings,updateInitialValues=True)
355
356 mbs.SolveDynamic(#solverType=exu.DynamicSolverType.TrapezoidalIndex2,
357 simulationSettings=simulationSettings)
358
359 #uTip = mbs.GetSensorValues(sensTipDispl)[1]
360 #print("nModes=", nModes, ", tip displacement=", uTip)
361
362 if useGraphics:
363 #SC.renderer.DoIdleTasks()
364 SC.renderer.Stop() #safely close rendering window!
365
366 mbs.PlotSensor(sensTipDispl,components=[0,1,2],title='arm tip displacements',closeAll=True)
367
368 #this loads the solution after simulation and allows visualization and storing animation frames!
369 mbs.SolutionViewer()