pendulumVerify.py

You can view and download this file on Github: pendulumVerify.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#
  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.itemInterface import *
 16from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 17import exudyn.graphics as graphics #only import if it does not conflict
 18from exudyn.FEM import *
 19from exudyn.graphicsDataUtilities import *
 20
 21SC = exu.SystemContainer()
 22mbs = SC.AddSystem()
 23
 24import numpy as np
 25
 26import time
 27
 28import exudyn.basicUtilities as eb
 29import exudyn.rigidBodyUtilities as rb
 30import exudyn.utilities as eu
 31
 32import numpy as np
 33
 34useGraphics = True
 35fileName = 'testData/pendulumSimple' #for load/save of FEM data
 36
 37#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 38#netgen/meshing part:
 39femInterface = FEMinterface()
 40
 41#geometrical parameters:
 42L = 400  #Length of plate (X)
 43ff=1 #factor for higher thickness
 44h = 8*ff #height of plate (Y)
 45w = 4*ff #width of plate  (Z)
 46nModes = 10
 47meshH = 2*ff
 48
 49
 50print("mesh h=", meshH)
 51
 52#steel:
 53rho = 2.7e-9 #tons/mm^3
 54Emodulus=70e3 #N/mm^2
 55g = 9810                #[mm/s^2]
 56nu=0.35
 57gForce=[0,-g,0]
 58
 59V = L*h*w
 60mass = V*rho
 61print("V=", V, " (mm^3), mass=", mass*1000, "(kg)")
 62
 63useGravity = True
 64A=w*h
 65q=rho*A*g
 66EI = Emodulus*w*h**3/12.
 67F=1 #N
 68
 69if useGravity:
 70    tipDisp = q*L**4/(8*EI)
 71    print("tip displ (weight)=", tipDisp)
 72else:
 73    tipDisp = F*L**3/(3*EI)
 74    print("tip displ (tip F )=", tipDisp)
 75
 76
 77#h*w=8*4, nu=0.35, E=70e3:
 78#F=1
 79#meshH=2:   w = 1.5989535
 80#meshH=1:   w = 1.73735541
 81#meshH=0.5: w = 1.77126479
 82#analytical:w = 1.78571428
 83#weight:
 84#meshH=2:   w = 0.20311787
 85#meshH=1:   w = 0.220752717
 86#analytical:w = 0.2270314285714286
 87
 88meshCreated = False
 89
 90#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 91if True: #needs netgen/ngsolve to be installed to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
 92
 93    import ngsolve as ngs
 94    import netgen
 95    from netgen.meshing import *
 96
 97    from netgen.geom2d import unit_square
 98    #import netgen.libngpy as libng
 99    from netgen.csg import *
100
101    geo = CSGeometry()
102
103    #plate
104    block = OrthoBrick(Pnt(0,-0.5*h, -0.5*w),Pnt(L, 0.5*h, 0.5*w))
105
106    geo.Add(block)
107
108    mesh = ngs.Mesh( geo.GenerateMesh(maxh=meshH))
109    mesh.Curve(1)
110
111    if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
112        # import netgen
113        import netgen.gui
114        ngs.Draw(mesh)
115        netgen.Redraw()
116
117    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
118    #Use femInterface to import femInterface model and create FFRFreducedOrder object
119    eigenModesComputed = False
120    femInterface.ImportMeshFromNGsolve(mesh, density=rho, youngsModulus=Emodulus, poissonsRatio=nu,
121                              computeEigenmodes=eigenModesComputed, verbose=False, excludeRigidBodyModes = 6,
122                              numberOfModes = nModes, maxEigensolveIterations=20)
123
124    P1=[0,0,0]
125    nodesLeft = femInterface.GetNodesInPlane(P1, [1,0,0])
126    #print("boundary nodes bolt=", nodesLeft)
127    nodesLeftLen = len(nodesLeft)
128    nodesLeftWeights = np.array((1./nodesLeftLen)*np.ones(nodesLeftLen))
129
130    P2=[L,0,0]
131    nodesRight = femInterface.GetNodesInPlane(P2, [1,0,0])
132    #print("boundary nodes bolt=", nodesRight)
133    nodesRightLen = len(nodesRight)
134    nodesRightWeights = np.array((1./nodesRightLen)*np.ones(nodesRightLen))
135
136    print("nNodes=",femInterface.NumberOfNodes())
137
138    strMode = ''
139    boundaryList = [nodesLeft, nodesRight]
140
141    print("compute HCB modes... ")
142    start_time = time.time()
143    femInterface.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList,
144                                  nEigenModes=nModes,
145                                  useSparseSolver=True,
146                                  computationMode = HCBstaticModeSelection.RBE2)
147
148    print("eigen freq.=", femInterface.GetEigenFrequenciesHz())
149    print("HCB modes needed %.3f seconds" % (time.time() - start_time))
150
151
152
153    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
154    #compute stress modes for postprocessing (inaccurate for coarse meshes, just for visualization):
155    if False:
156        mat = KirchhoffMaterial(Emodulus, nu, rho)
157        varType = exu.OutputVariableType.StressLocal
158        #varType = exu.OutputVariableType.StrainLocal
159        print("ComputePostProcessingModes ... (may take a while)")
160        start_time = time.time()
161        femInterface.ComputePostProcessingModes(material=mat,
162                                       outputVariableType=varType)
163        print("   ... needed %.3f seconds" % (time.time() - start_time))
164        SC.visualizationSettings.contour.reduceRange=False
165        SC.visualizationSettings.contour.outputVariable = varType
166        SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
167    else:
168        SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
169        SC.visualizationSettings.contour.outputVariableComponent = 1
170
171    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
172    print("create CMS element ...")
173    cms = ObjectFFRFreducedOrderInterface(femInterface)
174
175    objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
176                                                  initialVelocity=[0,0,0],
177                                                  initialAngularVelocity=[0,0,0],
178                                                  color=[0.9,0.9,0.9,1.],
179                                                  )
180
181    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
182    #add markers and joints
183    nodeDrawSize = 1 #for joint drawing
184
185
186    #mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
187
188    if True:
189
190        oGround = mbs.AddObject(ObjectGround(referencePosition= P1))
191
192        #compute offset of nodes at plane (because the average does not give [0,0,0]):
193        pOff = np.zeros(3)
194        nodes = femInterface.nodes['Position']
195        for i in nodesLeft:
196            pOff += nodes[i]
197        pOff *= 1/len(nodesLeft)
198
199        #create marker:
200        altApproach = True
201        mLeft = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
202                                                      meshNodeNumbers=np.array(nodesLeft), #these are the meshNodeNumbers
203                                                      offset=-pOff,
204                                                      useAlternativeApproach=altApproach,
205                                                      weightingFactors=nodesLeftWeights))
206
207        lockedAxes=[1,1,1,1,1,1]
208        if True:
209
210            mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround,
211                                                        localPosition=P1,
212                                                        visualization=VMarkerBodyRigid(show=True)))
213            mbs.AddObject(GenericJoint(markerNumbers=[mGround, mLeft],
214                                        constrainedAxes = lockedAxes,
215                                        visualization=VGenericJoint(show=False, axesRadius=0.1*w, axesLength=0.1*h)))
216
217
218    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
219    #animate modes
220    SC.visualizationSettings.markers.show = True
221    SC.visualizationSettings.markers.defaultSize=1
222    SC.visualizationSettings.markers.drawSimplified = False
223
224    SC.visualizationSettings.loads.show = False
225    SC.visualizationSettings.loads.drawSimplified = False
226    SC.visualizationSettings.loads.defaultSize=10
227    SC.visualizationSettings.loads.defaultRadius = 0.1
228    SC.visualizationSettings.openGL.multiSampling=4
229    SC.visualizationSettings.openGL.lineWidth=2
230
231    if False: #activate to animate modes
232        from exudyn.interactive import AnimateModes
233        mbs.Assemble()
234        SC.visualizationSettings.nodes.show = False
235        SC.visualizationSettings.openGL.showFaceEdges = True
236        SC.visualizationSettings.openGL.multiSampling=4
237        SC.visualizationSettings.openGL.lineWidth=2
238        SC.visualizationSettings.window.renderWindowSize = [1600,1080]
239        SC.visualizationSettings.contour.showColorBar = False
240        SC.visualizationSettings.general.textSize = 16
241
242        #%%+++++++++++++++++++++++++++++++++++++++
243        #animate modes of ObjectFFRFreducedOrder (only needs generic node containing modal coordinates)
244        SC.visualizationSettings.general.autoFitScene = False #otherwise, model may be difficult to be moved
245
246        nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
247        AnimateModes(SC, mbs, nodeNumber, period=0.1, showTime=False, renderWindowText='Hurty-Craig-Bampton: 2 x 6 static modes and 8 eigenmodes\n')
248        import sys
249        sys.exit()
250
251    oFFRF = objFFRF['oFFRFreducedOrder']
252    if useGravity:
253        #add gravity (not necessary if user functions used)
254        mBody = mbs.AddMarker(MarkerBodyMass(bodyNumber=oFFRF))
255        mbs.AddLoad(LoadMassProportional(markerNumber=mBody, loadVector= gForce))
256    else:
257        mRight = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
258                                                      meshNodeNumbers=np.array(nodesRight), #these are the meshNodeNumbers
259                                                      useAlternativeApproach=altApproach,
260                                                      weightingFactors=nodesRightWeights))
261        mbs.AddLoad(LoadForceVector(markerNumber=mRight, loadVector=[0,-F,0]))
262
263    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
264    fileDir = 'solution/'
265    # sensBolt = mbs.AddSensor(SensorMarker(markerNumber=mBolt,
266    #                                       fileName=fileDir+'hingePartBoltPos'+str(nModes)+strMode+'.txt',
267    #                                       outputVariableType = exu.OutputVariableType.Position))
268    # sensBushing= mbs.AddSensor(SensorMarker(markerNumber=mBushing,
269    #                                       fileName=fileDir+'hingePartBushingPos'+str(nModes)+strMode+'.txt',
270    #                                       outputVariableType = exu.OutputVariableType.Position))
271    nTip = femInterface.GetNodeAtPoint([L,0.5*h,0.5*w])
272    sensTip = mbs.AddSensor(SensorSuperElement(bodyNumber=oFFRF,
273                                                     meshNodeNumber=nTip,
274                                          fileName=fileDir+'displacementTip.txt',
275                                          outputVariableType = exu.OutputVariableType.DisplacementLocal))
276
277    # print("tip0=",mbs.GetSensorValues(sensTip))
278    mbs.Assemble()
279
280    simulationSettings = exu.SimulationSettings()
281
282    SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
283    SC.visualizationSettings.nodes.drawNodesAsPoint = False
284    SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
285
286    SC.visualizationSettings.nodes.show = False
287    SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
288    SC.visualizationSettings.nodes.basisSize = 0.12
289    SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
290
291    SC.visualizationSettings.openGL.showFaceEdges = True
292    SC.visualizationSettings.openGL.showFaces = True
293
294    SC.visualizationSettings.sensors.show = True
295    SC.visualizationSettings.sensors.drawSimplified = False
296    SC.visualizationSettings.sensors.defaultSize = 2
297    SC.visualizationSettings.loads.defaultSize = 10
298
299
300    simulationSettings.solutionSettings.solutionInformation = "pendulum verification"
301
302    h=0.25e-3*4
303    tEnd = 0.25*8
304
305    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
306    simulationSettings.timeIntegration.endTime = tEnd
307    simulationSettings.solutionSettings.writeSolutionToFile = False
308    simulationSettings.timeIntegration.verboseMode = 1
309    #simulationSettings.timeIntegration.verboseModeFile = 3
310    simulationSettings.timeIntegration.newton.useModifiedNewton = True
311
312    simulationSettings.solutionSettings.sensorsWritePeriod = h
313
314    simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
315    #simulationSettings.displayStatistics = True
316    #simulationSettings.displayComputationTime = True
317
318    #create animation:
319    # simulationSettings.solutionSettings.recordImagesInterval = 0.005
320    # SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
321    SC.visualizationSettings.window.renderWindowSize=[1920,1080]
322    SC.visualizationSettings.openGL.multiSampling = 4
323
324    useGraphics=False
325    if useGraphics:
326        if useGraphics:
327            SC.visualizationSettings.general.autoFitScene=False
328
329            exu.StartRenderer()
330            if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view
331
332            mbs.WaitForUserToContinue() #press space to continue
333
334    #SC.RedrawAndSaveImage()
335    if False:
336        mbs.SolveDynamic(simulationSettings=simulationSettings)
337    else:
338        mbs.SolveStatic(simulationSettings=simulationSettings)
339
340    # print("tip1=",mbs.GetSensorValues(sensTip))
341
342    if useGraphics:
343        SC.WaitForRenderEngineStopFlag()
344        exu.StopRenderer() #safely close rendering window!
345
346data = np.loadtxt('solution/displacementTip.txt', comments='#', delimiter=',')
347print("tip disp=", data[-1,1:])