superElementRigidJointTest.py
You can view and download this file on Github: superElementRigidJointTest.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
18import numpy as np
19
20testIsActive = exu.sys.get('testIsActive', False)
21
22SC = exu.SystemContainer()
23mbs = SC.AddSystem()
24
25#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
26#Use FEMinterface to import FEM model and create FFRFreducedOrder object
27fem = FEMinterface()
28inputFileName = 'testData/rotorDiscTest' #runTestSuite.py is at another directory
29
30nodes=fem.ImportFromAbaqusInputFile(inputFileName+'.inp', typeName='Instance', name='rotor-1')
31
32fem.ReadMassMatrixFromAbaqus(inputFileName+'MASS1.mtx')
33fem.ReadStiffnessMatrixFromAbaqus(inputFileName+'STIF1.mtx')
34fem.ScaleStiffnessMatrix(1e-3) #for larger deformations, stiffness is reduced to 1%
35
36nodeNumberUnbalance = 9 #on disc, max y-value
37unbalance = 0.1
38fem.AddNodeMass(nodeNumberUnbalance, unbalance)
39
40nModes = 20
41fem.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
42#print("eigen freq.=", fem.GetEigenFrequenciesHz())
43
44cms = ObjectFFRFreducedOrderInterface(fem)
45
46#user functions should be defined outside of class:
47def UFmassFFRFreducedOrder(mbs, t, itemIndex, qReduced, qReduced_t):
48 return cms.UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
49
50def UFforceFFRFreducedOrder(mbs, t, itemIndex, qReduced, qReduced_t):
51 return cms.UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
52
53objFFRF = cms.AddObjectFFRFreducedOrderWithUserFunctions(exu, mbs, positionRef=[0,0,0], rotationMatrixRef=np.eye(3),
54 initialVelocity=[0,0,0], initialAngularVelocity=[0,0,0*50*2*pi],
55 gravity = [0,-9.81,0],
56 UFforce=UFforceFFRFreducedOrder, UFmassMatrix=UFmassFFRFreducedOrder,
57 color=[0.1,0.9,0.1,1.])
58cms2 = ObjectFFRFreducedOrderInterface(fem)
59#user functions should be defined outside of class:
60def UFmassFFRFreducedOrder2(mbs, t, itemIndex, qReduced, qReduced_t):
61 return cms2.UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
62
63def UFforceFFRFreducedOrder2(mbs, t, itemIndex, qReduced, qReduced_t):
64 return cms2.UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
65
66objFFRF2 = cms2.AddObjectFFRFreducedOrderWithUserFunctions(exu, mbs, positionRef=[0,0,0.5], rotationMatrixRef=np.eye(3),
67 initialVelocity=[0,0,0], initialAngularVelocity=[0,0,0*50*2*pi],
68 gravity = [0,-9.81,0],
69 UFforce=UFforceFFRFreducedOrder2, UFmassMatrix=UFmassFFRFreducedOrder2,
70 color=[0.1,0.9,0.1,1.])
71
72#draw one mode:
73#mbs.SetNodeParameter(nodeNumber = objFFRF['nGenericODE2'], parameterName='initialCoordinates', value=[0.1]+[0]*(nModes-1))
74# mbs.SetNodeParameter(nodeNumber = objFFRF2['nRigidBody'], parameterName='initialCoordinates', value=[0,0,0.5]+[0]*4)
75
76#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
77#add markers and joints
78nodeDrawSize = 0.0025 #for joint drawing
79
80pLeft = [0,0,0]
81pRight = [0,0,0.5]
82pLeft2 = [0,0,0.5]
83nMid = fem.GetNodeAtPoint([0,0,0.25])
84#print("nMid=",nMid)
85
86mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
87oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
88
89mGroundPosLeft = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=pLeft))
90mGroundPosRight = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=pRight))
91
92#torque on reference frame:
93#mbs.AddLoad(Torque(markerNumber=mRB, loadVector=[0,0,100*2*pi]))
94
95if False: #OPTIONAL: lock rigid body motion of reference frame (for tests):
96 mbs.AddObject(GenericJoint(markerNumbers=[mGround, mRB], constrainedAxes=[1,1,1, 1,1,0]))
97
98#++++++++++++++++++++++++++++++++++++++++++
99#find nodes at left and right surface:
100nodeListLeft = fem.GetNodesInPlane(pLeft, [0,0,1])
101nodeListRight = fem.GetNodesInPlane(pRight, [0,0,1])
102#print("nodeListLeft=",nodeListLeft)
103#nLeft = fem.GetNodeAtPoint(pLeft)
104#nRight = fem.GetNodeAtPoint(pRight)
105
106
107lenLeft = len(nodeListLeft)
108lenRight = len(nodeListRight)
109weightsLeft = np.array((1./lenLeft)*np.ones(lenLeft))
110weightsRight = np.array((1./lenRight)*np.ones(lenRight))
111
112addSupports = True
113if addSupports:
114 k = 2e8*10*0.01 #joint stiffness
115 d = k*0.01 #joint damping
116
117 useGenericJoint = True
118
119# mLeft = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=objFFRF['oFFRFreducedOrder'],
120# meshNodeNumbers=np.array(nodeListLeft), #these are the meshNodeNumbers
121# weightingFactors=weightsLeft))
122 mRight = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=objFFRF['oFFRFreducedOrder'],
123 meshNodeNumbers=np.array(nodeListRight), #these are the meshNodeNumbers
124 weightingFactors=weightsRight))
125 mLeft2 = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=objFFRF2['oFFRFreducedOrder'],
126 meshNodeNumbers=np.array(nodeListLeft), #these are the meshNodeNumbers
127 weightingFactors=weightsLeft))
128
129 mLeftRigid = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
130 #referencePosition=pLeft, #deprecated
131 useAlternativeApproach=True,
132 meshNodeNumbers=np.array(nodeListLeft), #these are the meshNodeNumbers
133 weightingFactors=weightsLeft))
134 mRightRigid = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'],
135 #referencePosition=pRight, #deprecated
136 useAlternativeApproach=True,
137 meshNodeNumbers=np.array(nodeListRight), #these are the meshNodeNumbers
138 weightingFactors=weightsRight))
139 mLeftRigid2 = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF2['oFFRFreducedOrder'],
140 #referencePosition=pLeft, #deprecated
141 useAlternativeApproach=True,
142 meshNodeNumbers=np.array(nodeListLeft), #these are the meshNodeNumbers
143 weightingFactors=weightsLeft))
144 if useGenericJoint:
145# oSJleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mLeftRigid, mGroundPosLeft],
146# stiffness=[k,k,k], damping=[d,d,d]))
147# oSJright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mRight,mGroundPosRight],
148# stiffness=[k,k,0], damping=[d,d,d]))
149# oSJleft = mbs.AddObject(RigidBodySpringDamper(markerNumbers=[mLeftRigid, mGroundPosLeft],
150# stiffness=0.1*k*np.eye(6), damping=0.01*d*np.eye(6)))
151 oSJleft = mbs.AddObject(GenericJoint(markerNumbers=[mLeftRigid, mGroundPosLeft], constrainedAxes=[1,1,1,1,1,1],
152 visualization=VGenericJoint(axesRadius=0.02)))
153 oSJleft2 = mbs.AddObject(GenericJoint(markerNumbers=[mRightRigid, mLeftRigid2], constrainedAxes=[1,1,1*1,1,1,1],
154 visualization=VGenericJoint(axesRadius=0.02)))
155
156 # oSJright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mRightRigid, mLeftRigid2],
157 # stiffness=[k,k,k], damping=[d,d,d]))
158 # oSJright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mRight, mLeft2],
159 # stiffness=[k,k,k], damping=[d,d,d]))
160
161 else:
162 oSJleft = mbs.AddObject(SphericalJoint(markerNumbers=[mGroundPosLeft,mLeft], visualization=VObjectJointSpherical(jointRadius=nodeDrawSize)))
163 oSJright= mbs.AddObject(SphericalJoint(markerNumbers=[mGroundPosRight,mRight], visualization=VObjectJointSpherical(jointRadius=nodeDrawSize)))
164
165
166fileDir = 'solution/'
167sDisp=mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'], meshNodeNumber=nMid, #meshnode number!
168 storeInternal=True,#fileName=fileDir+'nMidDisplacementCMS'+str(nModes)+'Test.txt',
169 outputVariableType = exu.OutputVariableType.Displacement))
170
171sAngVel=mbs.AddSensor(SensorNode(nodeNumber=objFFRF['nRigidBody'],
172 storeInternal=True,#fileName=fileDir+'nRigidBodyAngVelCMS'+str(nModes)+'Test.txt',
173 outputVariableType = exu.OutputVariableType.AngularVelocity))
174
175mbs.Assemble()
176
177simulationSettings = exu.SimulationSettings()
178
179SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
180SC.visualizationSettings.nodes.drawNodesAsPoint = False
181SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
182
183SC.visualizationSettings.nodes.show = True
184SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
185SC.visualizationSettings.nodes.basisSize = 0.12
186SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
187
188SC.visualizationSettings.view0.scene.showFaceEdges = True
189SC.visualizationSettings.view0.scene.showFaces = True
190
191SC.visualizationSettings.sensors.show = True
192SC.visualizationSettings.sensors.drawSimplified = False
193SC.visualizationSettings.sensors.defaultSize = 0.01
194SC.visualizationSettings.markers.drawSimplified = False
195SC.visualizationSettings.markers.show = True
196SC.visualizationSettings.markers.defaultSize = 0.01
197
198SC.visualizationSettings.loads.drawSimplified = False
199
200SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
201SC.visualizationSettings.contour.outputVariableComponent = 1 #y-component
202#SC.visualizationSettings.contour.automaticRange = False
203SC.visualizationSettings.contour.reduceRange = False
204#SC.visualizationSettings.contour.maxValue = 0
205#SC.visualizationSettings.contour.minValue = -1
206
207simulationSettings.solution.file.information = "ObjectFFRFreducedOrder test"
208
209h=1e-3
210tEnd = 0.005 #standard:0.005
211if testIsActive:
212 #test suite:
213 simulationSettings.solution.file.write = False
214 tEnd = 0.005
215 h=1e-3
216
217simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
218simulationSettings.timeIntegration.endTime = tEnd
219simulationSettings.solution.file.writePeriod = h
220simulationSettings.timeIntegration.verboseMode = 1
221#simulationSettings.timeIntegration.verboseModeFile = 0
222simulationSettings.timeIntegration.newton.useModifiedNewton = True
223
224simulationSettings.solution.sensors.writePeriod = h
225simulationSettings.solution.file.name = "solution/coordinatesSolutionCMStest.txt"
226
227useIndex2 = False
228simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = useIndex2
229simulationSettings.timeIntegration.generalizedAlpha.useNewmark = useIndex2
230
231simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5 #SHOULD work with 0.9 as well
232#simulationSettings.show.statistics = True
233#simulationSettings.show.computationTime = True
234
235#create animation:
236#simulationSettings.solution.recordImagesInterval = 0.0002
237#SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
238
239if not testIsActive:
240 SC.renderer.Start()
241 if 'lastRenderState' in vars():
242 SC.renderer.SetState(lastRenderState) #load last model view
243
244 SC.renderer.DoIdleTasks() #press space to continue
245
246mbs.SolveDynamic(simulationSettings)
247
248if not testIsActive:
249 SC.renderer.DoIdleTasks()
250 SC.renderer.Stop() #safely close rendering window!
251 lastRenderState = SC.renderer.GetState() #store model view for next simulation
252
253
254#data = np.loadtxt(fileDir+'nMidDisplacementCMS'+str(nModes)+'Test.txt', comments='#', delimiter=',')
255data=mbs.GetSensorStoredData(sDisp)
256result = abs(data).sum()
257#pos = mbs.GetObjectOutputBody(objFFRF['oFFRFreducedOrder'],exu.OutputVariableType.Position, localPosition=[0,0,0])
258exu.Print('solution of superElementRigidJointTest=',result)
259
260exu.sys['testResult'] = result
261
262##++++++++++++++++++++++++++++++++++++++++++++++q+++++++
263#plot results
264if not testIsActive:
265
266
267 mbs.PlotSensor(sDisp, components=1, closeAll=True, labels=['uMid,linear'])