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