CMSexampleCourse.py
You can view and download this file on Github: CMSexampleCourse.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Create flexible multibody system
5#
6# Author: Johannes Gerstmayr
7# Date: 2021-09-10
8# Python version: Python 3.7, 64bits, Anaconda3 + ngsolve + webgui_jupyter_widgets
9# Jupyter: requires upgrade of scipy and uninstall and install tk (tkinter)
10# 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.
11#
12#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
13
14
15import exudyn as exu
16from exudyn.utilities import *
17from exudyn.FEM import *
18import numpy as np
19from math import sqrt, sin, cos, pi
20
21doMeshing = True #set false if mesh shall be loaded
22useHCBmodes = True #Hurty-Craig-Bampton modes
23computeStresses = False #takes some time
24loadStresses = False #set True only if already computed previously; may lead to severe problems if wrong modes are loaded!!!!
25
26fileName = 'testData/FMBStest1' #for load/save of FEM data
27
28
29# # Parameter Definition
30# ![example geometry](exampleBody.png "Geometry")
31
32#flexible body dimensions:
33femInterface = FEMinterface()
34#geometrical parameters
35ri = 0.010 #radius of hole/bolt
36ra = 0.016 #outer radius
37t = 0.020 #part thickness, bolt length
38f = 0.004 #radius of part rounding
39f2 = 0.002 #radius of bolt rounding
40L = 0.250 #distance hole/bolt
41
42hp = ra*sqrt(0.5)-f*(1-sqrt(0.5)) #height of bar
43d1 = (ra+1*f)*sqrt(0.5) #length of rounded part
44L1 = L-2*d1 #length of rectangular section
45
46#number of eigenmodes
47nModes = 8
48
49rho = 1000
50Emodulus=1e8
51nu=0.3
52
53
54# # Create FEM mesh in Netgen
55
56if doMeshing: #needs netgen/ngsolve to be installed to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
57
58 from netgen.occ import *
59 #from ngsolve.webgui import Draw #in Jupyter
60 from ngsolve import Mesh, Draw
61
62 #create part geometry
63 wp = WorkPlane(gp_Ax3(p=(0,0,-1.5*t), n=Z, h=X))
64
65 wp.Rotate(90).MoveTo(-ra,0).Arc(ra,-135).Arc(f,45).Line(0.5*L1).Line(0.5*L1).Arc(f,45).Arc(ra,-135)
66 wp.Arc(ra,-135).Arc(f,45).Line(L1).Arc(f,45).Arc(ra,-135)
67 wp.Close().Reverse()
68
69 p1 = wp.Face().Extrude(t)
70
71 #bolt:
72 wp2 = WorkPlane(Axes(p=(0,0,-0.5*t), n=X, h=Y))
73 f22 = np.sqrt(2)*f2
74 #wp2.MoveTo(0,0).LineTo(ri,0).LineTo(ri,t).LineTo(0,t).LineTo(0,0) #without rounding
75 wp2.MoveTo(0,0).Line(ri+f2).Rotate(180).Arc(f2,-90).Line(t-2*f2).Rotate(45).Line(f22).Rotate(45).Line(ri-f2).Rotate(90).Line(t).Close() #with rounding+chamfer
76 axis2 = Axis((0,0,0),Z)
77 p2 = wp2.Face().Revolve(axis2,360)
78
79 #hole:
80 wp3 = WorkPlane(Axes(p=(L,0,-1.5*t), n=X, h=Y))
81 #wp2.MoveTo(0,0).LineTo(ri,0).LineTo(ri,t).LineTo(0,t).LineTo(0,0) #without rounding
82 wp3.MoveTo(0,0).Line(ri+f2).Rotate(135).Line(f22).Rotate(-45).Line(t-2*f2).Rotate(-45).Line(f22).Rotate(135).Line(ri+f2).Rotate(90).Line(t).Close() #with rounding
83 axis3 = Axis((L,0,0),Z)
84 p3 = wp3.Face().Revolve(axis3,360)
85
86 p1 = p1 + p2
87 p1 = p1 - p3
88
89 #for geometry check:
90 # box = Box((0,0,0), (ri,ri,ri)) #show (0,0,0)
91 # p1 = p1+ box
92 geo = OCCGeometry( p1 )
93
94 #Jupyter, webgui, draw geometry
95 #NEEDS: pip install webgui_jupyter_widgets
96 from netgen.webgui import Draw as DrawGeo
97 #DrawGeo(geo.shape)
98
99 #generate mesh:
100 from ngsolve.webgui import Draw
101 mesh = Mesh(geo.GenerateMesh(maxh=1.5*f2))
102 #Jupyter, webgui, draw mesh
103 Draw(mesh)
104
105
106
107# # Import mesh into Exudyn
108SC = exu.SystemContainer()
109mbs = SC.AddSystem()
110
111if doMeshing: #needs netgen/ngsolve to be installed to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
112 #save FEM mesh
113 femInterface.ImportMeshFromNGsolve(mesh, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
114 femInterface.SaveToFile(fileName)
115else:
116 femInterface.LoadFromFile(fileName)
117print("number of nodes = ", femInterface.NumberOfNodes())
118
119
120# In[12]:
121
122
123femInterface.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
124if False: #activate to animate modes
125 from exudyn.interactive import AnimateModes
126 mbs.Reset()
127 cms = ObjectFFRFreducedOrderInterface(femInterface)
128
129 objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
130 initialVelocity=[0,0,0],
131 initialAngularVelocity=[0,0,0],
132 color=[0.1,0.9,0.1,1.],
133 )
134 mbs.Assemble()
135 SC.visualizationSettings.nodes.show = False
136 SC.visualizationSettings.openGL.showFaceEdges = True
137 SC.visualizationSettings.openGL.multiSampling=4
138 SC.visualizationSettings.openGL.lineWidth=2
139 SC.visualizationSettings.window.renderWindowSize = [1600,1080]
140
141 #%%+++++++++++++++++++++++++++++++++++++++
142 SC.visualizationSettings.general.autoFitScene = False #otherwise, model may be difficult to be moved
143
144 nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
145 AnimateModes(SC, mbs, nodeNumber, period=0.1,
146 scaleAmplitude = 0.02,
147 showTime=False, renderWindowText='Show modes\n',
148 runOnStart=True)
149
150
151# # Define interfaces
152addSensors = True
153pLeft = [0,0,0] #midpoint of bolt
154pRight = [L,0,-t] #midpoint of hole
155pMid = [0.5*L,hp,-0.5*t] #midpoint of bar
156pTip = [L+ra,0,-0.5*t] #midpoint of bar
157
158#%%
159if addSensors:
160 nMid = femInterface.GetNodeAtPoint(pMid, tolerance=1e-2) #tip node (do not use midpoint, as this may not be a mesh node ...)
161 print("pMid=",pMid,", nMid=",nMid)
162 nTip = femInterface.GetNodeAtPoint(pTip, tolerance=1e-2) #tip node (do not use midpoint, as this may not be a mesh node ...)
163 print("pTip=",pTip,", nTip=",nTip)
164
165tV = np.array([0,0,0.5*t])
166nodesLeft = femInterface.GetNodesOnCylinder(pLeft-tV, pLeft+tV, ri)
167# print('nodesLeft=',nodesLeft)
168nodesRight = femInterface.GetNodesOnCylinder(pRight-tV, pRight+tV, ri)
169# print('nodesRight=',nodesRight)
170
171lenNodesLeft = len(nodesLeft)
172weightsNodesLeft = np.array((1./lenNodesLeft)*np.ones(lenNodesLeft))
173
174lenNodesRight = len(nodesRight)
175weightsNodesRight = np.array((1./lenNodesRight)*np.ones(lenNodesRight))
176
177boundaryList = [nodesLeft, nodesRight] #second boudary (right plane) not needed ...
178
179
180# # Compute eigenmodes
181#remark: ComputeEigenmodes requires upgrade of scipy (python -m pip install --upgrade scipy) as compared to Anaconda installation...
182import time
183
184print("compute modes... ")
185start_time = time.time()
186
187if useHCBmodes:
188 femInterface.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList,
189 nEigenModes=nModes,
190 useSparseSolver=True,
191 computationMode = HCBstaticModeSelection.RBE2)
192else:
193 femInterface.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
194
195print("computation of modes needed %.3f seconds" % (time.time() - start_time))
196print("eigen freq.=", femInterface.GetEigenFrequenciesHz())
197
198
199
200
201# # Compute stresses
202
203femModesName = fileName+'modes'
204if useHCBmodes:
205 femModesName+='HCB'
206
207varType = exu.OutputVariableType.StressLocal
208
209if computeStresses:
210 mat = KirchhoffMaterial(Emodulus, nu, rho)
211 #varType = exu.OutputVariableType.StrainLocal
212 print("ComputePostProcessingModes ... (may take a while)")
213 start_time = time.time()
214 femInterface.ComputePostProcessingModes(material=mat,
215 outputVariableType=varType)
216 print(" ... needed %.3f seconds" % (time.time() - start_time))
217 SC.visualizationSettings.contour.reduceRange=False
218 SC.visualizationSettings.contour.outputVariable = varType
219 SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
220 #save modes + stresses
221 femInterface.SaveToFile(femModesName)
222else:
223 if loadStresses:
224 femInterface.LoadFromFile(femModesName)
225 SC.visualizationSettings.contour.outputVariable = varType
226 SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
227
228
229# # Setup flexible body in exudyn
230cms = ObjectFFRFreducedOrderInterface(femInterface)
231
232objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
233 initialVelocity=[0,0,0],
234 initialAngularVelocity=[0,0,0],
235 color=[0.1,0.9,0.1,1.],
236 )
237
238
239# # Visualize modes
240if False:
241 from exudyn.interactive import AnimateModes
242 mbs.Assemble()
243 SC.visualizationSettings.nodes.show = False
244 SC.visualizationSettings.openGL.showFaceEdges = True
245 SC.visualizationSettings.openGL.multiSampling=4
246 #SC.visualizationSettings.window.renderWindowSize = [1600,1080]
247 SC.visualizationSettings.general.autoFitScene = False #otherwise, model may be difficult to be moved
248
249 nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
250 AnimateModes(SC, mbs, nodeNumber, scaleAmplitude = 0.1, runOnStart = True)
251
252
253# # Add gravity
254
255# In[11]:
256
257
258#add gravity (not necessary if user functions used)
259oFFRF = objFFRF['oFFRFreducedOrder']
260mBody = mbs.AddMarker(MarkerBodyMass(bodyNumber=oFFRF))
261mbs.AddLoad(LoadMassProportional(markerNumber=mBody, loadVector= [0,-9.81,0]))
262
263
264# # Add joint constraint
265
266# In[12]:
267
268
269#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
270#add markers and joints
271
272#mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
273oGround = mbs.AddObject(ObjectGround(referencePosition = [0,0,0]))
274mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oGround,
275 localPosition = pLeft))
276
277#add marker
278#NOTE: offset is added in order to compensate for small errors in average node position
279# because mesh is not fully symmetric and the average node position does not match
280# the desired (body-fixed) joint position [0,0,0]; if not used, a small initial jump
281# happens during simulation when the body moves to the constrained position
282mLeft = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
283 meshNodeNumbers=np.array(nodesLeft), #these are the meshNodeNumbers
284 weightingFactors=weightsNodesLeft,
285 offset=-femInterface.GetNodePositionsMean(nodesLeft),
286 ))
287oJoint = mbs.AddObject(GenericJoint(markerNumbers=[mGround, mLeft],
288 constrainedAxes = [1,1,1,1,1,0],
289 visualization=VGenericJoint(axesRadius=0.05*ri, axesLength=1.5*t)))
290# oJoint = mbs.AddObject(RevoluteJointZ(markerNumbers=[mGround, mLeft],
291# visualization=VRevoluteJointZ(axisRadius=0.05*ri, axisLength=1.5*t)))
292
293if False: #if this is used, remove offset in MarkerSuperElementRigid above
294 #alternative to offset above: compensate joint offset by computation of current displacement in joint: (if not done in MarkerSuperElementRigid)
295 mbs.Assemble() #initialize system to compute joint offset
296 jointOffset = mbs.GetObjectOutput(oJoint,exu.OutputVariableType.DisplacementLocal)
297 print('jointOffset=',jointOffset)
298
299 mbs.SetMarkerParameter(mLeft.GetIndex(), 'offset', list(-jointOffset)) #compensate offset; mLeft.GetIndex() because of BUG751
300
301 #now check new offset:
302 mbs.Assemble() #initialize system to compute joint offset
303 jointOffset = mbs.GetObjectOutput(oJoint,exu.OutputVariableType.DisplacementLocal)
304 print('jointOffset=',jointOffset)
305
306
307# # Add sensors
308fileDir = 'solution/'
309if addSensors:
310 sMidDispl = mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'],
311 meshNodeNumber=nMid, #meshnode number!
312 fileName=fileDir+'uMid'+str(nModes)+'modes.txt',
313 outputVariableType = exu.OutputVariableType.Displacement))
314 sTipDispl = mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'],
315 meshNodeNumber=nTip, #meshnode number!
316 fileName=fileDir+'uTip'+str(nModes)+'modes.txt',
317 outputVariableType = exu.OutputVariableType.Displacement))
318
319
320# # Set up visualization
321# (not needed)
322nodeDrawSize = 0.0025 #for joint drawing
323SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
324SC.visualizationSettings.nodes.drawNodesAsPoint = False
325SC.visualizationSettings.connectors.defaultSize = nodeDrawSize
326
327SC.visualizationSettings.nodes.show = False
328SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
329SC.visualizationSettings.nodes.basisSize = t*4
330SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
331
332SC.visualizationSettings.openGL.showFaceEdges = True
333SC.visualizationSettings.openGL.showFaces = True
334
335SC.visualizationSettings.sensors.show = True
336SC.visualizationSettings.sensors.drawSimplified = False
337SC.visualizationSettings.sensors.defaultSize = nodeDrawSize*2
338SC.visualizationSettings.markers.drawSimplified = False
339SC.visualizationSettings.markers.show = False
340SC.visualizationSettings.markers.defaultSize = nodeDrawSize*2
341
342SC.visualizationSettings.loads.drawSimplified = False
343SC.visualizationSettings.loads.defaultSize = t*3
344SC.visualizationSettings.loads.defaultRadius = 0.05*t
345
346SC.visualizationSettings.window.renderWindowSize=[1280,720]
347SC.visualizationSettings.openGL.multiSampling = 4
348
349#create animation:
350# simulationSettings.solutionSettings.recordImagesInterval = 0.005
351# SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
352
353
354# # Set up simulation
355mbs.Assemble() #initialize bodies, assemble system; necessary to simulate
356
357simulationSettings = exu.SimulationSettings()
358simulationSettings.solutionSettings.solutionInformation = "ObjectFFRFreducedOrder test"
359
360h=1e-3
361tEnd = 2
362
363simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
364simulationSettings.timeIntegration.endTime = tEnd
365simulationSettings.solutionSettings.writeSolutionToFile = True
366simulationSettings.solutionSettings.solutionWritePeriod = h
367simulationSettings.timeIntegration.verboseMode = 1
368#simulationSettings.timeIntegration.verboseModeFile = 3
369simulationSettings.timeIntegration.newton.useModifiedNewton = True
370
371simulationSettings.solutionSettings.sensorsWritePeriod = h
372
373simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
374#simulationSettings.displayStatistics = True
375simulationSettings.displayComputationTime = True
376
377
378# # Start renderer and Simulate
379
380lifeVisualization = True
381
382if lifeVisualization:
383 SC.visualizationSettings.general.autoFitScene=False #if reloaded view settings
384 exu.StartRenderer()
385 if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view
386 mbs.WaitForUserToContinue() #press space to continue
387
388mbs.SolveDynamic(#solverType=exu.DynamicSolverType.TrapezoidalIndex2,
389 simulationSettings=simulationSettings)
390
391if addSensors:
392 uTip = mbs.GetSensorValues(sMidDispl)
393 print("nModes=", nModes, ", mid displacement=", uTip)
394
395if lifeVisualization:
396 SC.WaitForRenderEngineStopFlag()
397 exu.StopRenderer() #safely close rendering window!
398
399
400# # 3D rendering of FMBS
401if False: #use this to reload the solution and use SolutionViewer
402 SC.visualizationSettings.general.autoFitScene=False #if reloaded view settings
403
404
405 mbs.SolutionViewer() #can also be entered in IPython ...
406
407
408# # Plot sensor
409
410if addSensors:
411
412 mbs.PlotSensor(sensorNumbers=[sMidDispl,sMidDispl,sMidDispl], components=[0,1,2])