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