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