ObjectFFRFconvergenceTestBeam.py

You can view and download this file on Github: ObjectFFRFconvergenceTestBeam.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test for meshing with NETGEN and import of FEM model;
  5#           Model is a simple flexible pendulum meshed with tet elements;
  6#           Note that the model is overly flexible (linearized strain assumption not valid),
  7#           but it should serve as a demonstration of the FFRFreducedOrder modeling
  8#
  9# Author:   Johannes Gerstmayr
 10# Date:     2021-02-05
 11#
 12# 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.
 13#
 14#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 15
 16
 17import exudyn as exu
 18from exudyn.itemInterface import *
 19from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 20import exudyn.graphics as graphics #only import if it does not conflict
 21from exudyn.FEM import *
 22from exudyn.graphicsDataUtilities import *
 23
 24SC = exu.SystemContainer()
 25mbs = SC.AddSystem()
 26
 27import numpy as np
 28
 29import time
 30import timeit
 31
 32import exudyn.basicUtilities as eb
 33import exudyn.rigidBodyUtilities as rb
 34import exudyn.utilities as eu
 35
 36import numpy as np
 37
 38useGraphics = True
 39fileName = 'testData/netgenLshape' #for load/save of FEM data
 40
 41
 42#+++++++++++++++++++++++++++++++++++++++++++++++++++++
 43#netgen/meshing part:
 44fem = FEMinterface()
 45#standard:
 46a = 0.1 #height/width of beam
 47b = a
 48h = 0.2*a
 49L = 1     #Length of beam
 50nModes = 10
 51
 52#plate:
 53# a = 0.025 #height/width of beam
 54# b = 0.4
 55# L = 1     #Length of beam
 56# h = 0.6*a
 57# nModes = 40
 58
 59rho = 1000
 60Emodulus=1e8
 61nu=0.3
 62#analytical solution due to self-weight:
 63g=9.81
 64EI = Emodulus*b*a**3/12
 65rhoA = rho*b*a
 66uTip = rhoA*g * L**4/(8*EI) #Gieck, 29th edition, 1989, page 166 (P13)
 67
 68doStatic = True
 69
 70meshCreated = False
 71
 72
 73
 74#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 75if True: #needs netgen/ngsolve to be installed to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
 76
 77    import ngsolve as ngs
 78    from netgen.geom2d import unit_square
 79    import netgen.libngpy as libng
 80    from netgen.csg import *
 81
 82    geo = CSGeometry()
 83
 84    block = OrthoBrick(Pnt(0,-a*0.5,-b*0.5),Pnt(L,a*0.5,b*0.5))
 85    geo.Add(block)
 86
 87    #Draw (geo)
 88
 89    mesh = ngs.Mesh( geo.GenerateMesh(maxh=h))
 90    mesh.Curve(1)
 91
 92    if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
 93        import netgen.gui
 94        Draw (mesh)
 95        netgen.Redraw()
 96
 97    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 98    #Use fem to import FEM model and create FFRFreducedOrder object
 99    fem.ImportMeshFromNGsolve(mesh, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
100    meshCreated  = True
101    if (h==a): #save only if it has smaller size
102        fem.SaveToFile(fileName)
103
104    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
105if True: #now import mesh as mechanical model to EXUDYN
106    if not meshCreated: fem.LoadFromFile(fileName)
107
108    #fix left plane
109    supportMidPoint = [0,0,0]
110
111    nodeListSupport = fem.GetNodesInPlane([0,0,0], [1,0,0])
112    lenNodeListSupport = len(nodeListSupport)
113    weightsNodeListSupport = np.array((1./lenNodeListSupport)*np.ones(lenNodeListSupport))
114
115    nodeListTip= fem.GetNodesInPlane([L,0,0], [1,0,0])
116    lenNodeListTip= len(nodeListTip)
117    weightsNodeListTip= np.array((1./lenNodeListTip)*np.ones(lenNodeListTip))
118
119    print("nNodes=",fem.NumberOfNodes())
120
121    strMode = ''
122    if False: #pure eigenmodes
123        print("compute eigen modes... ")
124        start_time = time.time()
125        fem.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
126        print("eigen modes computation needed %.3f seconds" % (time.time() - start_time))
127        print("eigen freq.=", fem.GetEigenFrequenciesHz())
128
129    else:
130        strMode = 'HCB'
131        boundaryList = [nodeListSupport]
132        #boundaryList = [nodeListTip,nodeListSupport]
133        #boundaryList = [nodeListSupport,nodeListTip] #gives approx. same result as before
134
135        print("compute HCB modes... ")
136        start_time = time.time()
137        fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList,
138                                      nEigenModes=nModes,
139                                      useSparseSolver=True,
140                                      computationMode = HCBstaticModeSelection.RBE2)
141
142        print("eigen freq.=", fem.GetEigenFrequenciesHz())
143        print("HCB modes needed %.3f seconds" % (time.time() - start_time))
144
145
146
147
148    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
149    #compute stress modes:
150    varType = exu.OutputVariableType.Displacement
151    if False:
152        mat = KirchhoffMaterial(Emodulus, nu, rho)
153        varType = exu.OutputVariableType.StressLocal
154        #varType = exu.OutputVariableType.StrainLocal
155        print("ComputePostProcessingModes ... (may take a while)")
156        start_time = time.time()
157        fem.ComputePostProcessingModes(material=mat,
158                                       outputVariableType=varType)
159        print("--- %s seconds ---" % (time.time() - start_time))
160
161    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
162    print("create CMS element ...")
163    cms = ObjectFFRFreducedOrderInterface(fem)
164
165    objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
166                                                  initialVelocity=[0,0,0], initialAngularVelocity=[0,0,0],
167                                                  color=[0.1,0.9,0.1,1.])
168
169
170    #add gravity (not necessary if user functions used)
171    oFFRF = objFFRF['oFFRFreducedOrder']
172    mBody = mbs.AddMarker(MarkerBodyMass(bodyNumber=oFFRF))
173    mbs.AddLoad(LoadMassProportional(markerNumber=mBody, loadVector= [0,-g,0]))
174
175    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
176    #add markers and joints
177    nodeDrawSize = 0.0025 #for joint drawing
178
179
180    #++++++++++++++++++++++++++++++++++++++++++
181    nTip = fem.GetNodeAtPoint([L,-a*0.5,-b*0.5]) #tip node
182
183    if True:
184        oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
185
186        #altApproach = True
187        lockedAxes=[1,1,1,1,1*1,1]
188
189        mSupport = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
190                                                        meshNodeNumbers=np.array(nodeListSupport), #these are the meshNodeNumbers
191                                                        weightingFactors=weightsNodeListSupport))
192        mGroundSupport = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround,
193                                                    localPosition=supportMidPoint,
194                                                    visualization=VMarkerBodyRigid(show=True)))
195        mbs.AddObject(GenericJoint(markerNumbers=[mGroundSupport, mSupport],
196                                    constrainedAxes = lockedAxes,
197                                    visualization=VGenericJoint(show=False, axesRadius=0.1*b, axesLength=0.1*b)))
198
199
200
201    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
202    fileDir = 'solution/'
203    sTip = mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'],
204                                     meshNodeNumber=nTip, #meshnode number!
205                             fileName=fileDir+'nMidDisplacementCMS'+str(nModes)+'Test.txt',
206                             outputVariableType = exu.OutputVariableType.Displacement))
207
208    mbs.Assemble()
209
210    simulationSettings = exu.SimulationSettings()
211
212    SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
213    SC.visualizationSettings.nodes.drawNodesAsPoint = False
214    SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
215
216    SC.visualizationSettings.nodes.show = False
217    SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
218    SC.visualizationSettings.nodes.basisSize = 0.12
219    SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
220
221    SC.visualizationSettings.openGL.showFaceEdges = True
222    SC.visualizationSettings.openGL.showFaces = True
223
224    SC.visualizationSettings.sensors.show = True
225    SC.visualizationSettings.sensors.drawSimplified = False
226    SC.visualizationSettings.sensors.defaultSize = 0.01
227    SC.visualizationSettings.markers.drawSimplified = False
228    SC.visualizationSettings.markers.show = False
229    SC.visualizationSettings.markers.defaultSize = 0.01
230
231    SC.visualizationSettings.loads.drawSimplified = False
232
233    # SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
234    # SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
235    SC.visualizationSettings.contour.reduceRange=False
236    SC.visualizationSettings.contour.outputVariable = varType
237    SC.visualizationSettings.contour.outputVariableComponent = 1 #y-component
238
239    simulationSettings.solutionSettings.solutionInformation = "ObjectFFRFreducedOrder test"
240
241    h=0.25e-3
242    tEnd = 0.12
243
244    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
245    simulationSettings.timeIntegration.endTime = tEnd
246    simulationSettings.solutionSettings.writeSolutionToFile = False
247    simulationSettings.timeIntegration.verboseMode = 1
248    #simulationSettings.timeIntegration.verboseModeFile = 3
249    simulationSettings.timeIntegration.newton.useModifiedNewton = True
250
251    simulationSettings.solutionSettings.sensorsWritePeriod = h
252
253    simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8 #SHOULD work with 0.9 as well
254    #simulationSettings.displayStatistics = True
255    #simulationSettings.displayComputationTime = True
256
257    #create animation:
258    # simulationSettings.solutionSettings.recordImagesInterval = 0.005
259    # SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
260    SC.visualizationSettings.window.renderWindowSize=[1920,1080]
261    SC.visualizationSettings.openGL.multiSampling = 4
262
263    useGraphics = True
264    if useGraphics:
265        SC.visualizationSettings.general.autoFitScene=False
266
267        exu.StartRenderer()
268        if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view
269
270        mbs.WaitForUserToContinue() #press space to continue
271
272
273    if doStatic:
274        mbs.SolveStatic(simulationSettings=simulationSettings, showHints=True)
275        uTipNum = -mbs.GetSensorValues(sTip)[1]
276        print("uTipNumerical=", uTipNum, ", uTipAnalytical=",uTip)
277        #HCB:
278        #h=0.2*a:
279        #uTipNumerical= 0.013870128561063066 , uTipAnalytical= 0.014714999999999999
280        #h=0.1*a:
281        #uTipNumerical= 0.014492581916470945 , uTipAnalytical= 0.014714999999999999
282
283        #10 modes HCB (two interfaces:support/tip):
284        #uTipNumerical= 0.013862260226352854
285        #10 modes HCB (two interfaces:tip/support):
286        #uTipNumerical= 0.013867428098277693 (nearly identical with other case)
287    else:
288        mbs.SolveDynamic(#solverType=exu.DynamicSolverType.TrapezoidalIndex2,
289                          simulationSettings=simulationSettings)
290        uTipNum = -mbs.GetSensorValues(sTip)[1]
291        print("uTipNumerical=", uTipNum)
292        #10 eigenmodes:
293        #uTipNumerical= 0.005782728750346744
294        #100 eigenmodes:
295        #uTipNumerical= 0.020578363592264157
296        #2 modes HCB:
297        #uTipNumerical= 0.022851728744898644
298        #10 modes HCB:
299        #uTipNumerical= 0.022998972747996865
300
301
302    if useGraphics:
303        SC.WaitForRenderEngineStopFlag()
304        exu.StopRenderer() #safely close rendering window!