objectGenericODE2Test.py
You can view and download this file on Github: objectGenericODE2Test.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test for ObjectGenericODE2 with python user function for linear and rotor dynamics
5# This test represents a FEM model of a rotor, which has elastic supports and rotation is locked
6# We compute eigenmodes, compute the linear response as well as the response of the rotor with gyroscopic terms
7#
8# Author: Johannes Gerstmayr
9# Date: 2020-05-13
10#
11# 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.
12#
13#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
14
15import exudyn as exu
16from exudyn.utilities import *
17from exudyn.FEM import *
18
19useGraphics = True #without test
20#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
21#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
22try: #only if called from test suite
23 from modelUnitTests import exudynTestGlobals #for globally storing test results
24 useGraphics = exudynTestGlobals.useGraphics
25except:
26 class ExudynTestGlobals:
27 pass
28 exudynTestGlobals = ExudynTestGlobals()
29#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
30
31SC = exu.SystemContainer()
32mbs = SC.AddSystem()
33
34import numpy as np
35accumulatedError = 0
36#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
37#Use FEMinterface to import FEM model and create FFRFreducedOrder object
38fem = FEMinterface()
39inputFileName = 'testData/rotorDiscTest' #runTestSuite.py is at another directory
40
41#useGraphics = False
42
43nodes=fem.ImportFromAbaqusInputFile(inputFileName+'.inp', typeName='Instance', name='rotor-1')
44nNodes = len(nodes)
45nODE2 = nNodes*3 #total number of ODE2 coordinates in FEM system; size of M and K
46
47fem.ReadMassMatrixFromAbaqus(inputFileName+'MASS1.mtx')
48fem.ReadStiffnessMatrixFromAbaqus(inputFileName+'STIF1.mtx')
49fem.ScaleStiffnessMatrix(1e-2) #for larger deformations, stiffness is reduced to 1%
50
51#+++++++++++ add elastic supports to fem ==> compute correct eigen frequencies
52pLeft = [0,0,0]
53pRight = [0,0,0.5]
54nLeft = fem.GetNodeAtPoint(pLeft)
55nRight = fem.GetNodeAtPoint(pRight)
56kJoint = 2e8 #joint stiffness
57dJoint = kJoint*0.01 #joint damping
58
59fem.AddElasticSupportAtNode(nLeft, springStiffness=[kJoint,kJoint,kJoint])
60fem.AddElasticSupportAtNode(nRight, springStiffness=[kJoint,kJoint,kJoint])
61
62#+++++++++++ compute eigenmodes for comparison
63nModes = 8
64fem.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = False)
65exu.Print("eigen freq.=", fem.GetEigenFrequenciesHz()[0:6+nModes].round(4)) #mode 0 is rigid body mode (free rotation)!
66exu.Print("eigen freq. first mode =", fem.GetEigenFrequenciesHz()[1]) #mode1 with supports: 57.6317863976366Hz; free-free mode6: sparse: 104.63701326020315, dense: 104.63701326063597
67accumulatedError += 1e-2*(fem.GetEigenFrequenciesHz()[1]/57.63178639764625 - 1.) #check mode (with supports); this is subject to small variations between 32 and 64bit! ==>*1e-2
68
69#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
70#create generic object for rotor:
71forwardWhirl = True #test: True; switch this flag to turn on rotordynamics effects
72backwardWhirl = False #test: False; switch this flag to turn on rotordynamics effects
73excitationSign = 1
74if backwardWhirl: excitationSign = -1
75
76forceVector = [0,-1.,0]*nNodes #static force vector, optional; add some force in y-direction; could also use node mass to compute force due to weight
77fUnbalance = 2000 #fictive force, not depending on frequency
78
79nForce = fem.GetNodeAtPoint([0,0,0.15])#position where unbalance (excitation) force acts
80exu.Print("excitation node=", nForce)
81
82fRotX = np.array([0]*nODE2)
83fRotX[nForce*3+0] = fUnbalance
84fRotY = np.array([0]*nODE2)
85fRotY[nForce*3+1] = fUnbalance
86
87#sweep parameters:
88t1 = 5 #end time of sweep (t0=0)
89f0 = 50 #start frequency
90f1 = 250 #terminal frequency
91omega1=np.pi*2.*f1
92
93G = fem.GetGyroscopicMatrix(rotationAxis=2, sparse=False) #gyroscopic matrix for rotordynamics effects
94
95useSparseG=True #speed up computation, using sparse matrix in user function
96if useSparseG:
97 from scipy.sparse import csr_matrix
98 G=csr_matrix(G) #convert to sparse matrix
99
100def UFforce(mbs, t, itemIndex, q, q_t):
101 #print("UFforce")
102 omega = 2.*np.pi*FrequencySweep(t, t1, f0,f1)
103 fact = omega/omega1
104 force = (fact*SweepCos(t, t1, f0, f1))*fRotY + (excitationSign*fact*SweepSin(t, t1, f0, f1))*fRotX #excitationSign: '+' = forward whirl, '-' = backward whirl
105
106 if forwardWhirl or backwardWhirl:
107 #omegaQ_t = omega * np.array(q_t)
108 force -= G @ (omega * np.array(q_t)) #negative sign: because term belongs to left-hand-side!!!
109 #force -= omega*(G @q_t) #negative sign: because term belongs to left-hand-side!!!
110 return force
111
112#add nodes:
113nodeList = [] #create node list
114for node in fem.GetNodePositionsAsArray():
115 nodeList += [mbs.AddNode(NodePoint(referenceCoordinates=node))]
116
117#now add generic body built from FEM model with mass and stiffness matrix (optional damping could be added):
118oGenericODE2 = mbs.AddObject(ObjectGenericODE2(nodeNumbers = nodeList,
119 massMatrix=fem.GetMassMatrix(sparse=False),
120 stiffnessMatrix=fem.GetStiffnessMatrix(sparse=False),
121 forceVector=forceVector, forceUserFunction=UFforce,
122 visualization=VObjectGenericODE2(triangleMesh = fem.GetSurfaceTriangles(),
123 color=color4lightred)))
124
125#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
126#add markers and joints
127nodeDrawSize = 0.0025 #for joint drawing
128
129nMid = fem.GetNodeAtPoint([0,0,0.25])
130nTopMid = fem.GetNodeAtPoint([0., 0.05, 0.25]) #lock rotation of body
131exu.Print("nMid=",nMid)
132exu.Print("nTopMid=",nTopMid)
133
134nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0], visualization = VNodePointGround(show=False))) #ground node for coordinate constraint
135mGroundCoordinate = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
136
137#add constraint to avoid rotation of body
138mTopRight = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nTopMid, coordinate=0)) #x-coordinate of node at y-max
139mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundCoordinate, mTopRight]))
140
141oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
142
143mGroundPosLeft = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=pLeft))
144mGroundPosRight = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=pRight))
145
146#++++++++++++++++++++++++++++++++++++++++++
147#find nodes at left and right surface:
148nodeListLeft = fem.GetNodesInPlane(pLeft, [0,0,1])
149nodeListRight = fem.GetNodesInPlane(pRight, [0,0,1])
150
151
152lenLeft = len(nodeListLeft)
153lenRight = len(nodeListRight)
154weightsLeft = np.array((1./lenLeft)*np.ones(lenLeft))
155weightsRight = np.array((1./lenRight)*np.ones(lenRight))
156
157addDampers = True
158if addDampers:
159 for i in range(3):
160 mLeft = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nLeft, coordinate=i))
161 mRight = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nRight, coordinate=i))
162
163 mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mGroundCoordinate,mLeft],
164 stiffness=0, damping=dJoint))
165 if i != 2: #exclude double constraint in z-direction (axis)
166 mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mGroundCoordinate,mRight],
167 stiffness=0, damping=dJoint))
168
169
170addJoint = False #set this flag, if not adding supports to stiffness matrix in fem
171if addJoint:
172
173 useSpringDamper = True
174
175 mLeft = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=oGenericODE2,
176 meshNodeNumbers=np.array(nodeListLeft), #these are the meshNodeNumbers
177 weightingFactors=weightsLeft))
178 mRight = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=oGenericODE2,
179 meshNodeNumbers=np.array(nodeListRight), #these are the meshNodeNumbers
180 weightingFactors=weightsRight))
181
182 oSJleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mLeft, mGroundPosLeft],
183 stiffness=[kJoint,kJoint,kJoint], damping=[dJoint,dJoint,dJoint]))
184 oSJright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mRight,mGroundPosRight],
185 stiffness=[kJoint,kJoint,0], damping=[dJoint,dJoint,0]))
186
187
188
189fileDir = 'solution/'
190sDisp=mbs.AddSensor(SensorSuperElement(bodyNumber=oGenericODE2, meshNodeNumber=nMid, #meshnode number!
191 storeInternal=True,#fileName=fileDir+'nMidDisplacementLinearTest.txt',
192 outputVariableType = exu.OutputVariableType.Displacement))
193
194mbs.Assemble()
195
196simulationSettings = exu.SimulationSettings()
197
198SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
199SC.visualizationSettings.nodes.drawNodesAsPoint = False
200SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
201
202SC.visualizationSettings.bodies.show = True
203#SC.visualizationSettings.connectors.show = False
204
205SC.visualizationSettings.bodies.deformationScaleFactor = 10 #use this factor to scale the deformation of modes
206if SC.visualizationSettings.bodies.deformationScaleFactor !=1:
207 SC.visualizationSettings.nodes.show = False #nodes are not scaled
208
209SC.visualizationSettings.openGL.showFaceEdges = True
210SC.visualizationSettings.openGL.showFaces = True
211
212#SC.visualizationSettings.sensors.show = True
213#SC.visualizationSettings.sensors.drawSimplified = False
214SC.visualizationSettings.sensors.defaultSize = 0.01
215#SC.visualizationSettings.markers.drawSimplified = False
216#SC.visualizationSettings.markers.show = True
217SC.visualizationSettings.markers.defaultSize = 0.01
218
219SC.visualizationSettings.loads.drawSimplified = False
220
221SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
222SC.visualizationSettings.contour.outputVariableComponent = 1 #y-component
223
224simulationSettings.solutionSettings.solutionInformation = "ObjectGenericODE2 test"
225simulationSettings.solutionSettings.writeSolutionToFile=False
226
227h=1e-3
228tEnd = 0.05
229#if useGraphics:
230# tEnd = 0.1
231
232simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
233simulationSettings.timeIntegration.endTime = tEnd
234simulationSettings.solutionSettings.solutionWritePeriod = h
235simulationSettings.timeIntegration.verboseMode = 1
236#simulationSettings.timeIntegration.verboseModeFile = 3
237simulationSettings.timeIntegration.newton.useModifiedNewton = True
238
239simulationSettings.solutionSettings.sensorsWritePeriod = h
240
241simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5 #SHOULD work with 0.9 as well
242simulationSettings.displayStatistics = False
243simulationSettings.displayComputationTime = False
244
245#create animation:
246#simulationSettings.solutionSettings.recordImagesInterval = 0.0002
247#SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
248
249#useGraphics = True
250if useGraphics:
251 exu.StartRenderer()
252 if 'lastRenderState' in vars():
253 SC.SetRenderState(lastRenderState) #load last model view
254
255 #mbs.WaitForUserToContinue() #press space to continue
256
257
258mbs.SolveDynamic(simulationSettings)
259
260if useGraphics:
261 SC.WaitForRenderEngineStopFlag()
262 exu.StopRenderer() #safely close rendering window!
263 lastRenderState = SC.GetRenderState() #store model view for next simulation
264
265accumulatedError += mbs.GetNodeOutput(nMid,exu.OutputVariableType.Position)[0] #take x-coordinate of position
266
267exu.Print('solution of ObjectGenericODE2=',accumulatedError)
268
269exudynTestGlobals.testError = accumulatedError - (-2.2737401292182432e-05) #2020-05-18: -2.2737401292182432e-05
270exudynTestGlobals.testResult = accumulatedError
271
272##++++++++++++++++++++++++++++++++++++++++++++++q+++++++
273#plot results
274if useGraphics:
275
276
277 mbs.PlotSensor(sDisp, components=1, closeAll=True, labels=['uMid,linear'])