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:])