objectFFRFreducedOrderTest.py

You can view and download this file on Github: objectFFRFreducedOrderTest.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test for ObjectFFRFreducedOrder with python user function for reduced order equations of motion
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2020-05-13
  8#
  9# 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.
 10#
 11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 12
 13import exudyn as exu
 14from exudyn.utilities import *
 15from exudyn.FEM import *
 16
 17useGraphics = True #without test
 18#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 19#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 20try: #only if called from test suite
 21    from modelUnitTests import exudynTestGlobals #for globally storing test results
 22    useGraphics = exudynTestGlobals.useGraphics
 23except:
 24    class ExudynTestGlobals:
 25        pass
 26    exudynTestGlobals = ExudynTestGlobals()
 27#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 28
 29SC = exu.SystemContainer()
 30mbs = SC.AddSystem()
 31
 32import numpy as np
 33
 34
 35#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 36#Use FEMinterface to import FEM model and create FFRFreducedOrder object
 37fem = FEMinterface()
 38inputFileName = 'testData/rotorDiscTest' #runTestSuite.py is at another directory
 39#if useGraphics:
 40#    inputFileName = 'testData/rotorDiscTest'        #if executed in current directory
 41
 42nodes=fem.ImportFromAbaqusInputFile(inputFileName+'.inp', typeName='Instance', name='rotor-1')
 43
 44fem.ReadMassMatrixFromAbaqus(inputFileName+'MASS1.mtx')
 45fem.ReadStiffnessMatrixFromAbaqus(inputFileName+'STIF1.mtx')
 46fem.ScaleStiffnessMatrix(1e-2) #for larger deformations, stiffness is reduced to 1%
 47
 48#nodeNumberUnbalance = 9  #on disc, max y-value
 49nodeNumberUnbalance = fem.GetNodeAtPoint(point=[0. , 0.19598444, 0.15])
 50#exu.Print("nodeNumberUnbalance =",nodeNumberUnbalance)
 51unbalance = 0.1
 52fem.AddNodeMass(nodeNumberUnbalance, unbalance)
 53#print(fem.GetMassMatrix()[8*3:11*3,:])
 54
 55nModes = 8
 56fem.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
 57#print("eigen freq.=", fem.GetEigenFrequenciesHz())
 58
 59cms = ObjectFFRFreducedOrderInterface(fem)
 60
 61objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
 62                                        initialVelocity=[0,0,0], initialAngularVelocity=[0,0,50*2*pi],
 63                                        color=[0.1,0.9,0.1,1.])
 64
 65#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 66#add markers and joints
 67nodeDrawSize = 0.0025 #for joint drawing
 68
 69pLeft = [0,0,0]
 70pRight = [0,0,0.5]
 71nMid = fem.GetNodeAtPoint([0,0,0.25])
 72#print("nMid=",nMid)
 73
 74mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
 75oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
 76
 77mGroundPosLeft = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=pLeft))
 78mGroundPosRight = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=pRight))
 79
 80#torque on reference frame:
 81#mbs.AddLoad(Torque(markerNumber=mRB, loadVector=[0,0,100*2*pi]))
 82
 83if False: #OPTIONAL: lock rigid body motion of reference frame (for tests):
 84    mbs.AddObject(GenericJoint(markerNumbers=[mGround, mRB], constrainedAxes=[1,1,1, 1,1,0]))
 85
 86#++++++++++++++++++++++++++++++++++++++++++
 87#find nodes at left and right surface:
 88nodeListLeft = fem.GetNodesInPlane(pLeft, [0,0,1])
 89nodeListRight = fem.GetNodesInPlane(pRight, [0,0,1])
 90#nLeft = fem.GetNodeAtPoint(pLeft)
 91#nRight = fem.GetNodeAtPoint(pRight)
 92
 93
 94lenLeft = len(nodeListLeft)
 95lenRight = len(nodeListRight)
 96weightsLeft = np.array((1./lenLeft)*np.ones(lenLeft))
 97weightsRight = np.array((1./lenRight)*np.ones(lenRight))
 98
 99addSupports = True
100if addSupports:
101    k = 2e8     #joint stiffness
102    d = k*0.01  #joint damping
103
104    useSpringDamper = True
105
106    mLeft = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=objFFRF['oFFRFreducedOrder'],
107                                                    meshNodeNumbers=np.array(nodeListLeft), #these are the meshNodeNumbers
108                                                    weightingFactors=weightsLeft))
109    mRight = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=objFFRF['oFFRFreducedOrder'],
110                                                    meshNodeNumbers=np.array(nodeListRight), #these are the meshNodeNumbers
111                                                    weightingFactors=weightsRight))
112    if useSpringDamper:
113        oSJleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mLeft, mGroundPosLeft],
114                                            stiffness=[k,k,k], damping=[d,d,d]))
115        oSJright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mRight,mGroundPosRight],
116                                            stiffness=[k,k,0], damping=[d,d,d]))
117    else:
118        oSJleft = mbs.AddObject(SphericalJoint(markerNumbers=[mGroundPosLeft,mLeft], visualization=VObjectJointSpherical(jointRadius=nodeDrawSize)))
119        oSJright= mbs.AddObject(SphericalJoint(markerNumbers=[mGroundPosRight,mRight], visualization=VObjectJointSpherical(jointRadius=nodeDrawSize)))
120
121
122#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
123fileDir = 'solution/'
124sDisp=mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'], meshNodeNumber=nMid, #meshnode number!
125                         storeInternal=True,#fileName=fileDir+'nMidDisplacementCMS'+str(nModes)+'Test.txt',
126                         outputVariableType = exu.OutputVariableType.Displacement))
127
128sAngVel=mbs.AddSensor(SensorNode(nodeNumber=objFFRF['nRigidBody'],
129                         storeInternal=True,#fileName=fileDir+'nRigidBodyAngVelCMS'+str(nModes)+'Test.txt',
130                         outputVariableType = exu.OutputVariableType.AngularVelocity))
131
132mbs.Assemble()
133
134simulationSettings = exu.SimulationSettings()
135
136SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
137SC.visualizationSettings.nodes.drawNodesAsPoint = False
138SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
139
140SC.visualizationSettings.nodes.show = True
141SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
142SC.visualizationSettings.nodes.basisSize = 0.12
143SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
144
145SC.visualizationSettings.openGL.showFaceEdges = True
146SC.visualizationSettings.openGL.showFaces = True
147
148SC.visualizationSettings.sensors.show = True
149SC.visualizationSettings.sensors.drawSimplified = False
150SC.visualizationSettings.sensors.defaultSize = 0.01
151SC.visualizationSettings.markers.drawSimplified = False
152SC.visualizationSettings.markers.show = True
153SC.visualizationSettings.markers.defaultSize = 0.01
154
155SC.visualizationSettings.loads.drawSimplified = False
156
157SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
158SC.visualizationSettings.contour.outputVariableComponent = 1 #y-component
159
160simulationSettings.solutionSettings.solutionInformation = "ObjectFFRFreducedOrder test"
161
162h=1e-4
163tEnd = 0.01
164#if useGraphics:
165#    tEnd = 0.1
166
167simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
168simulationSettings.timeIntegration.endTime = tEnd
169simulationSettings.solutionSettings.solutionWritePeriod = h
170simulationSettings.timeIntegration.verboseMode = 1
171#simulationSettings.timeIntegration.verboseModeFile = 3
172simulationSettings.timeIntegration.newton.useModifiedNewton = True
173
174simulationSettings.solutionSettings.sensorsWritePeriod = h
175simulationSettings.solutionSettings.coordinatesSolutionFileName = "solution/coordinatesSolutionCMStest.txt"
176simulationSettings.solutionSettings.writeSolutionToFile=False
177
178simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5 #SHOULD work with 0.9 as well
179#simulationSettings.displayStatistics = True
180#simulationSettings.displayComputationTime = True
181
182#create animation:
183#simulationSettings.solutionSettings.recordImagesInterval = 0.0002
184#SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
185
186if useGraphics:
187    exu.StartRenderer()
188    if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view
189
190    mbs.WaitForUserToContinue() #press space to continue
191
192mbs.SolveDynamic(simulationSettings)
193
194
195# data = np.loadtxt(fileDir+'nMidDisplacementCMS'+str(nModes)+'Test.txt', comments='#', delimiter=',')
196data = mbs.GetSensorStoredData(sDisp)
197result = abs(data).sum()
198#pos = mbs.GetObjectOutputBody(objFFRF['oFFRFreducedOrder'],exu.OutputVariableType.Position, localPosition=[0,0,0])
199exu.Print('solution of ObjectFFRFreducedOrder=',result)
200
201#factor 0.05: make error smaller, as there are small changes for different runs (because of scipy sparse eigenvalue solver!)
202exudynTestGlobals.testError = 0.01*(result - (0.5354530110580623)) #2020-05-26(added EP-constraint): 0.5354530110580623; 2020-05-17 (tEnd=0.01, h=1e-4): 0.535452257303538
203exudynTestGlobals.testResult = 0.01*result
204
205if useGraphics:
206    SC.WaitForRenderEngineStopFlag()
207    exu.StopRenderer() #safely close rendering window!
208    lastRenderState = SC.GetRenderState() #store model view for next simulation
209
210#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
211#plot results
212if useGraphics:
213
214
215    mbs.PlotSensor([fileDir+'nMidDisplacementCMS8.txt',sDisp,fileDir+'nMidDisplacementFFRF.txt'],
216               components=1, closeAll=True)