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