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