ObjectFFRFconvergenceTestBeam.py
You can view and download this file on Github: ObjectFFRFconvergenceTestBeam.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test for meshing with NETGEN and import of FEM model;
5# Model is a simple flexible pendulum meshed with tet elements;
6# Note that the model is overly flexible (linearized strain assumption not valid),
7# but it should serve as a demonstration of the FFRFreducedOrder modeling
8#
9# Author: Johannes Gerstmayr
10# Date: 2021-02-05
11#
12# 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.
13#
14#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
15
16
17import exudyn as exu
18from exudyn.itemInterface import *
19from exudyn.utilities import *
20from exudyn.FEM import *
21from exudyn.graphicsDataUtilities import *
22
23SC = exu.SystemContainer()
24mbs = SC.AddSystem()
25
26import numpy as np
27
28import time
29import timeit
30
31import exudyn.basicUtilities as eb
32import exudyn.rigidBodyUtilities as rb
33import exudyn.utilities as eu
34
35import numpy as np
36
37useGraphics = True
38fileName = 'testData/netgenLshape' #for load/save of FEM data
39
40
41#+++++++++++++++++++++++++++++++++++++++++++++++++++++
42#netgen/meshing part:
43fem = FEMinterface()
44#standard:
45a = 0.1 #height/width of beam
46b = a
47h = 0.2*a
48L = 1 #Length of beam
49nModes = 10
50
51#plate:
52# a = 0.025 #height/width of beam
53# b = 0.4
54# L = 1 #Length of beam
55# h = 0.6*a
56# nModes = 40
57
58rho = 1000
59Emodulus=1e8
60nu=0.3
61#analytical solution due to self-weight:
62g=9.81
63EI = Emodulus*b*a**3/12
64rhoA = rho*b*a
65uTip = rhoA*g * L**4/(8*EI) #Gieck, 29th edition, 1989, page 166 (P13)
66
67doStatic = True
68
69meshCreated = False
70
71
72
73#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
74if True: #needs netgen/ngsolve to be installed to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
75
76 import ngsolve as ngs
77 from netgen.geom2d import unit_square
78 import netgen.libngpy as libng
79 from netgen.csg import *
80
81 geo = CSGeometry()
82
83 block = OrthoBrick(Pnt(0,-a*0.5,-b*0.5),Pnt(L,a*0.5,b*0.5))
84 geo.Add(block)
85
86 #Draw (geo)
87
88 mesh = ngs.Mesh( geo.GenerateMesh(maxh=h))
89 mesh.Curve(1)
90
91 if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
92 import netgen.gui
93 Draw (mesh)
94 netgen.Redraw()
95
96 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
97 #Use fem to import FEM model and create FFRFreducedOrder object
98 fem.ImportMeshFromNGsolve(mesh, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
99 meshCreated = True
100 if (h==a): #save only if it has smaller size
101 fem.SaveToFile(fileName)
102
103 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
104if True: #now import mesh as mechanical model to EXUDYN
105 if not meshCreated: fem.LoadFromFile(fileName)
106
107 #fix left plane
108 supportMidPoint = [0,0,0]
109
110 nodeListSupport = fem.GetNodesInPlane([0,0,0], [1,0,0])
111 lenNodeListSupport = len(nodeListSupport)
112 weightsNodeListSupport = np.array((1./lenNodeListSupport)*np.ones(lenNodeListSupport))
113
114 nodeListTip= fem.GetNodesInPlane([L,0,0], [1,0,0])
115 lenNodeListTip= len(nodeListTip)
116 weightsNodeListTip= np.array((1./lenNodeListTip)*np.ones(lenNodeListTip))
117
118 print("nNodes=",fem.NumberOfNodes())
119
120 strMode = ''
121 if False: #pure eigenmodes
122 print("compute eigen modes... ")
123 start_time = time.time()
124 fem.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
125 print("eigen modes computation needed %.3f seconds" % (time.time() - start_time))
126 print("eigen freq.=", fem.GetEigenFrequenciesHz())
127
128 else:
129 strMode = 'HCB'
130 boundaryList = [nodeListSupport]
131 #boundaryList = [nodeListTip,nodeListSupport]
132 #boundaryList = [nodeListSupport,nodeListTip] #gives approx. same result as before
133
134 print("compute HCB modes... ")
135 start_time = time.time()
136 fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList,
137 nEigenModes=nModes,
138 useSparseSolver=True,
139 computationMode = HCBstaticModeSelection.RBE2)
140
141 print("eigen freq.=", fem.GetEigenFrequenciesHz())
142 print("HCB modes needed %.3f seconds" % (time.time() - start_time))
143
144
145
146
147 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
148 #compute stress modes:
149 varType = exu.OutputVariableType.Displacement
150 if False:
151 mat = KirchhoffMaterial(Emodulus, nu, rho)
152 varType = exu.OutputVariableType.StressLocal
153 #varType = exu.OutputVariableType.StrainLocal
154 print("ComputePostProcessingModes ... (may take a while)")
155 start_time = time.time()
156 fem.ComputePostProcessingModes(material=mat,
157 outputVariableType=varType)
158 print("--- %s seconds ---" % (time.time() - start_time))
159
160 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
161 print("create CMS element ...")
162 cms = ObjectFFRFreducedOrderInterface(fem)
163
164 objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
165 initialVelocity=[0,0,0], initialAngularVelocity=[0,0,0],
166 color=[0.1,0.9,0.1,1.])
167
168
169 #add gravity (not necessary if user functions used)
170 oFFRF = objFFRF['oFFRFreducedOrder']
171 mBody = mbs.AddMarker(MarkerBodyMass(bodyNumber=oFFRF))
172 mbs.AddLoad(LoadMassProportional(markerNumber=mBody, loadVector= [0,-g,0]))
173
174 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
175 #add markers and joints
176 nodeDrawSize = 0.0025 #for joint drawing
177
178
179 #++++++++++++++++++++++++++++++++++++++++++
180 nTip = fem.GetNodeAtPoint([L,-a*0.5,-b*0.5]) #tip node
181
182 if True:
183 oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
184
185 #altApproach = True
186 lockedAxes=[1,1,1,1,1*1,1]
187
188 mSupport = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
189 meshNodeNumbers=np.array(nodeListSupport), #these are the meshNodeNumbers
190 weightingFactors=weightsNodeListSupport))
191 mGroundSupport = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround,
192 localPosition=supportMidPoint,
193 visualization=VMarkerBodyRigid(show=True)))
194 mbs.AddObject(GenericJoint(markerNumbers=[mGroundSupport, mSupport],
195 constrainedAxes = lockedAxes,
196 visualization=VGenericJoint(show=False, axesRadius=0.1*b, axesLength=0.1*b)))
197
198
199
200 #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
201 fileDir = 'solution/'
202 sTip = mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'],
203 meshNodeNumber=nTip, #meshnode number!
204 fileName=fileDir+'nMidDisplacementCMS'+str(nModes)+'Test.txt',
205 outputVariableType = exu.OutputVariableType.Displacement))
206
207 mbs.Assemble()
208
209 simulationSettings = exu.SimulationSettings()
210
211 SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
212 SC.visualizationSettings.nodes.drawNodesAsPoint = False
213 SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
214
215 SC.visualizationSettings.nodes.show = False
216 SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
217 SC.visualizationSettings.nodes.basisSize = 0.12
218 SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
219
220 SC.visualizationSettings.openGL.showFaceEdges = True
221 SC.visualizationSettings.openGL.showFaces = True
222
223 SC.visualizationSettings.sensors.show = True
224 SC.visualizationSettings.sensors.drawSimplified = False
225 SC.visualizationSettings.sensors.defaultSize = 0.01
226 SC.visualizationSettings.markers.drawSimplified = False
227 SC.visualizationSettings.markers.show = False
228 SC.visualizationSettings.markers.defaultSize = 0.01
229
230 SC.visualizationSettings.loads.drawSimplified = False
231
232 # SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
233 # SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
234 SC.visualizationSettings.contour.reduceRange=False
235 SC.visualizationSettings.contour.outputVariable = varType
236 SC.visualizationSettings.contour.outputVariableComponent = 1 #y-component
237
238 simulationSettings.solutionSettings.solutionInformation = "ObjectFFRFreducedOrder test"
239
240 h=0.25e-3
241 tEnd = 0.12
242
243 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
244 simulationSettings.timeIntegration.endTime = tEnd
245 simulationSettings.solutionSettings.writeSolutionToFile = False
246 simulationSettings.timeIntegration.verboseMode = 1
247 #simulationSettings.timeIntegration.verboseModeFile = 3
248 simulationSettings.timeIntegration.newton.useModifiedNewton = True
249
250 simulationSettings.solutionSettings.sensorsWritePeriod = h
251
252 simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8 #SHOULD work with 0.9 as well
253 #simulationSettings.displayStatistics = True
254 #simulationSettings.displayComputationTime = True
255
256 #create animation:
257 # simulationSettings.solutionSettings.recordImagesInterval = 0.005
258 # SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
259 SC.visualizationSettings.window.renderWindowSize=[1920,1080]
260 SC.visualizationSettings.openGL.multiSampling = 4
261
262 useGraphics = True
263 if useGraphics:
264 SC.visualizationSettings.general.autoFitScene=False
265
266 exu.StartRenderer()
267 if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view
268
269 mbs.WaitForUserToContinue() #press space to continue
270
271
272 if doStatic:
273 mbs.SolveStatic(simulationSettings=simulationSettings, showHints=True)
274 uTipNum = -mbs.GetSensorValues(sTip)[1]
275 print("uTipNumerical=", uTipNum, ", uTipAnalytical=",uTip)
276 #HCB:
277 #h=0.2*a:
278 #uTipNumerical= 0.013870128561063066 , uTipAnalytical= 0.014714999999999999
279 #h=0.1*a:
280 #uTipNumerical= 0.014492581916470945 , uTipAnalytical= 0.014714999999999999
281
282 #10 modes HCB (two interfaces:support/tip):
283 #uTipNumerical= 0.013862260226352854
284 #10 modes HCB (two interfaces:tip/support):
285 #uTipNumerical= 0.013867428098277693 (nearly identical with other case)
286 else:
287 mbs.SolveDynamic(#solverType=exu.DynamicSolverType.TrapezoidalIndex2,
288 simulationSettings=simulationSettings)
289 uTipNum = -mbs.GetSensorValues(sTip)[1]
290 print("uTipNumerical=", uTipNum)
291 #10 eigenmodes:
292 #uTipNumerical= 0.005782728750346744
293 #100 eigenmodes:
294 #uTipNumerical= 0.020578363592264157
295 #2 modes HCB:
296 #uTipNumerical= 0.022851728744898644
297 #10 modes HCB:
298 #uTipNumerical= 0.022998972747996865
299
300
301 if useGraphics:
302 SC.WaitForRenderEngineStopFlag()
303 exu.StopRenderer() #safely close rendering window!