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
18from exudyn.FEM 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 = 2*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
92 import ngsolve as ngs
93 import netgen
94 from netgen.meshing import *
95
96 from netgen.geom2d import unit_square
97 #import netgen.libngpy as libng
98 from netgen.csg import *
99
100 geo = CSGeometry()
101
102 #plate
103 block = OrthoBrick(Pnt(0,-0.5*h, -0.5*w),Pnt(L, 0.5*h, 0.5*w))
104
105 geo.Add(block)
106
107 mesh = ngs.Mesh( geo.GenerateMesh(maxh=meshH))
108 mesh.Curve(1)
109
110 if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
111 # import netgen
112 import netgen.gui
113 ngs.Draw(mesh)
114 netgen.Redraw()
115
116 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
117 #Use femInterface to import femInterface model and create FFRFreducedOrder object
118 eigenModesComputed = False
119 femInterface.ImportMeshFromNGsolve(mesh, density=rho, youngsModulus=Emodulus, poissonsRatio=nu,
120 computeEigenmodes=eigenModesComputed, verbose=False, excludeRigidBodyModes = 6,
121 numberOfModes = nModes, maxEigensolveIterations=20)
122
123 P1=[0,0,0]
124 nodesLeft = femInterface.GetNodesInPlane(P1, [1,0,0])
125 #print("boundary nodes bolt=", nodesLeft)
126 nodesLeftLen = len(nodesLeft)
127 nodesLeftWeights = np.array((1./nodesLeftLen)*np.ones(nodesLeftLen))
128
129 P2=[L,0,0]
130 nodesRight = femInterface.GetNodesInPlane(P2, [1,0,0])
131 #print("boundary nodes bolt=", nodesRight)
132 nodesRightLen = len(nodesRight)
133 nodesRightWeights = np.array((1./nodesRightLen)*np.ones(nodesRightLen))
134
135 print("nNodes=",femInterface.NumberOfNodes())
136
137 strMode = ''
138 boundaryList = [nodesLeft, nodesRight]
139
140 print("compute HCB modes... ")
141 start_time = time.time()
142 femInterface.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList,
143 nEigenModes=nModes,
144 useSparseSolver=True,
145 computationMode = HCBstaticModeSelection.RBE2)
146
147 print("eigen freq.=", femInterface.GetEigenFrequenciesHz())
148 print("HCB modes needed %.3f seconds" % (time.time() - start_time))
149
150
151
152 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
153 #compute stress modes for postprocessing (inaccurate for coarse meshes, just for visualization):
154 if False:
155 mat = KirchhoffMaterial(Emodulus, nu, rho)
156 varType = exu.OutputVariableType.StressLocal
157 #varType = exu.OutputVariableType.StrainLocal
158 print("ComputePostProcessingModes ... (may take a while)")
159 start_time = time.time()
160 femInterface.ComputePostProcessingModes(material=mat,
161 outputVariableType=varType)
162 print(" ... needed %.3f seconds" % (time.time() - start_time))
163 SC.visualizationSettings.contour.reduceRange=False
164 SC.visualizationSettings.contour.outputVariable = varType
165 SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
166 else:
167 SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
168 SC.visualizationSettings.contour.outputVariableComponent = 1
169
170 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
171 print("create CMS element ...")
172 cms = ObjectFFRFreducedOrderInterface(femInterface)
173
174 objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
175 initialVelocity=[0,0,0],
176 initialAngularVelocity=[0,0,0],
177 color=[0.9,0.9,0.9,1.],
178 )
179
180 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
181 #add markers and joints
182 nodeDrawSize = 1 #for joint drawing
183
184
185 #mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
186
187 if True:
188
189 oGround = mbs.AddObject(ObjectGround(referencePosition= P1))
190
191 #compute offset of nodes at plane (because the average does not give [0,0,0]):
192 pOff = np.zeros(3)
193 nodes = femInterface.nodes['Position']
194 for i in nodesLeft:
195 pOff += nodes[i]
196 pOff *= 1/len(nodesLeft)
197
198 #create marker:
199 altApproach = True
200 mLeft = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
201 meshNodeNumbers=np.array(nodesLeft), #these are the meshNodeNumbers
202 offset=-pOff,
203 useAlternativeApproach=altApproach,
204 weightingFactors=nodesLeftWeights))
205
206 lockedAxes=[1,1,1,1,1,1]
207 if True:
208
209 mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround,
210 localPosition=P1,
211 visualization=VMarkerBodyRigid(show=True)))
212 mbs.AddObject(GenericJoint(markerNumbers=[mGround, mLeft],
213 constrainedAxes = lockedAxes,
214 visualization=VGenericJoint(show=False, axesRadius=0.1*w, axesLength=0.1*h)))
215
216
217 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
218 #animate modes
219 SC.visualizationSettings.markers.show = True
220 SC.visualizationSettings.markers.defaultSize=1
221 SC.visualizationSettings.markers.drawSimplified = False
222
223 SC.visualizationSettings.loads.show = False
224 SC.visualizationSettings.loads.drawSimplified = False
225 SC.visualizationSettings.loads.defaultSize=10
226 SC.visualizationSettings.loads.defaultRadius = 0.1
227 SC.visualizationSettings.openGL.multiSampling=4
228 SC.visualizationSettings.openGL.lineWidth=2
229
230 if False: #activate to animate modes
231 from exudyn.interactive import AnimateModes
232 mbs.Assemble()
233 SC.visualizationSettings.nodes.show = False
234 SC.visualizationSettings.view0.scene.showFaceEdges = True
235 SC.visualizationSettings.openGL.multiSampling=4
236 SC.visualizationSettings.openGL.lineWidth=2
237 SC.visualizationSettings.view0.window.renderWindowSize = [1600,1080]
238 SC.visualizationSettings.contour.advanced.showColorBar = False
239 SC.visualizationSettings.view0.window.globalFontSize = 16
240
241 #%%+++++++++++++++++++++++++++++++++++++++
242 #animate modes of ObjectFFRFreducedOrder (only needs generic node containing modal coordinates)
243 SC.visualizationSettings.general.autoFitScene = False #otherwise, model may be difficult to be moved
244
245 nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
246 AnimateModes(SC, mbs, nodeNumber, period=0.1, showTime=False, renderWindowText='Hurty-Craig-Bampton: 2 x 6 static modes and 8 eigenmodes\n')
247 import sys
248 sys.exit()
249
250 oFFRF = objFFRF['oFFRFreducedOrder']
251 if useGravity:
252 #add gravity (not necessary if user functions used)
253 mBody = mbs.AddMarker(MarkerBodyMass(bodyNumber=oFFRF))
254 mbs.AddLoad(LoadMassProportional(markerNumber=mBody, loadVector= gForce))
255 else:
256 mRight = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
257 meshNodeNumbers=np.array(nodesRight), #these are the meshNodeNumbers
258 useAlternativeApproach=altApproach,
259 weightingFactors=nodesRightWeights))
260 mbs.AddLoad(LoadForceVector(markerNumber=mRight, loadVector=[0,-F,0]))
261
262 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
263 fileDir = 'solution/'
264 # sensBolt = mbs.AddSensor(SensorMarker(markerNumber=mBolt,
265 # fileName=fileDir+'hingePartBoltPos'+str(nModes)+strMode+'.txt',
266 # outputVariableType = exu.OutputVariableType.Position))
267 # sensBushing= mbs.AddSensor(SensorMarker(markerNumber=mBushing,
268 # fileName=fileDir+'hingePartBushingPos'+str(nModes)+strMode+'.txt',
269 # outputVariableType = exu.OutputVariableType.Position))
270 nTip = femInterface.GetNodeAtPoint([L,0.5*h,0.5*w])
271 sensTip = mbs.AddSensor(SensorSuperElement(bodyNumber=oFFRF,
272 meshNodeNumber=nTip,
273 fileName=fileDir+'displacementTip.txt',
274 outputVariableType = exu.OutputVariableType.DisplacementLocal))
275
276 # print("tip0=",mbs.GetSensorValues(sensTip))
277 mbs.Assemble()
278
279 simulationSettings = exu.SimulationSettings()
280
281 SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
282 SC.visualizationSettings.nodes.drawNodesAsPoint = False
283 SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
284
285 SC.visualizationSettings.nodes.show = False
286 SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
287 SC.visualizationSettings.nodes.basisSize = 0.12
288 SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
289
290 SC.visualizationSettings.view0.scene.showFaceEdges = True
291 SC.visualizationSettings.view0.scene.showFaces = True
292
293 SC.visualizationSettings.sensors.show = True
294 SC.visualizationSettings.sensors.drawSimplified = False
295 SC.visualizationSettings.sensors.defaultSize = 2
296 SC.visualizationSettings.loads.defaultSize = 10
297
298
299 simulationSettings.solution.file.information = "pendulum verification"
300
301 h=0.25e-3*4
302 tEnd = 0.25*8
303
304 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
305 simulationSettings.timeIntegration.endTime = tEnd
306 simulationSettings.solution.file.write = False
307 simulationSettings.timeIntegration.verboseMode = 1
308 #simulationSettings.timeIntegration.verboseModeFile = 3
309 simulationSettings.timeIntegration.newton.useModifiedNewton = True
310
311 simulationSettings.solution.sensors.writePeriod = h
312
313 simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
314 #simulationSettings.show.statistics = True
315 #simulationSettings.show.computationTime = True
316
317 #create animation:
318 # simulationSettings.solution.recordImagesInterval = 0.005
319 # SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
320 SC.visualizationSettings.view0.window.renderWindowSize=[1920,1080]
321 SC.visualizationSettings.openGL.multiSampling = 4
322
323 useGraphics=False
324 if useGraphics:
325 if useGraphics:
326 SC.visualizationSettings.general.autoFitScene=False
327
328 SC.renderer.Start()
329 SC.renderer.RestoreSavedState() #load last model view
330
331 SC.renderer.DoIdleTasks() #press space to continue
332
333 if False:
334 mbs.SolveDynamic(simulationSettings=simulationSettings)
335 else:
336 mbs.SolveStatic(simulationSettings=simulationSettings)
337
338 # print("tip1=",mbs.GetSensorValues(sensTip))
339
340 if useGraphics:
341 SC.renderer.DoIdleTasks()
342 SC.renderer.Stop() #safely close rendering window!
343
344data = np.loadtxt(OutputFilePath('solution/displacementTip.txt'), comments='#', delimiter=',')
345print("tip disp=", data[-1,1:])