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