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