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