objectFFRFreducedOrderAccelerations.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test accelerations output for ObjectFFRFreducedOrder
  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
 14#a local run must not open a plot window when windows are suppressed (#2477)
 15if exu.special.userInterface.suppressPlots: import matplotlib; matplotlib.use('Agg')
 16from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 17import exudyn.graphics as graphics
 18from exudyn.FEM import *
 19
 20import numpy as np
 21
 22testIsActive = exu.sys.get('testIsActive', False)
 23
 24SC = exu.SystemContainer()
 25mbs = SC.AddSystem()
 26testIsActive = True #without test
 27
 28
 29#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 30#Use FEMinterface to import FEM model and create FFRFreducedOrder object
 31fem = FEMinterface()
 32#inputFileName = 'C:/DATA/cpp/EXUDYN_git/main/pythonDev/TestModels/testData/rotorDiscTest' #runTestSuite.py is at another directory
 33inputFileName = 'testData/rotorDiscTest' #runTestSuite.py is at another directory
 34#if not testIsActive:
 35#    inputFileName = 'testData/rotorDiscTest'        #if executed in current directory
 36
 37nodes=fem.ImportFromAbaqusInputFile(inputFileName+'.inp', typeName='Instance', name='rotor-1')
 38
 39fem.ReadMassMatrixFromAbaqus(inputFileName+'MASS1.mtx')
 40fem.ReadStiffnessMatrixFromAbaqus(inputFileName+'STIF1.mtx')
 41fem.ScaleStiffnessMatrix(6e-4) #for larger deformations, stiffness is reduced to 1%
 42
 43#nodeNumberUnbalance = 9  #on disc, max y-value
 44nodeNumberUnbalance = fem.GetNodeAtPoint(point=[0. , 0.19598444, 0.15])
 45#exu.Print("nodeNumberUnbalance =",nodeNumberUnbalance)
 46unbalance = 0.1
 47fem.AddNodeMass(nodeNumberUnbalance, unbalance)
 48
 49nModes = 8
 50fem.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = False)
 51#print("eigen freq.=", fem.GetEigenFrequenciesHz())
 52
 53cms = ObjectFFRFreducedOrderInterface(fem)
 54
 55#user functions should be defined outside of class:
 56def UFmassFFRFreducedOrder(mbs, t, itemIndex, qReduced, qReduced_t):
 57    return cms.UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
 58
 59def UFforceFFRFreducedOrder(mbs, t, itemIndex, qReduced, qReduced_t):
 60    return cms.UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
 61
 62objFFRF = cms.AddObjectFFRFreducedOrderWithUserFunctions(exu, mbs, positionRef=[0,0,0], rotationMatrixRef=np.eye(3), 
 63                                              initialVelocity=[0,0,0], initialAngularVelocity=[0,0,50*2*pi],
 64                                              gravity = [0,-0*9.81,0],
 65                                              UFforce=UFforceFFRFreducedOrder, UFmassMatrix=UFmassFFRFreducedOrder,
 66                                              color=[0.1,0.9,0.1,1.])
 67
 68#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 69#add markers and joints
 70nodeDrawSize = 0.0025 #for joint drawing
 71
 72pLeft = [0,0,0]
 73pRight = [0,0,0.5]
 74nMid = fem.GetNodeAtPoint([0,0.05,0.25])
 75#exu.Print("nMid=",nMid)
 76
 77mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
 78oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
 79
 80mGroundPosLeft = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=pLeft))
 81mGroundPosRight = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=pRight))
 82
 83#torque on reference frame:
 84#mbs.AddLoad(Torque(markerNumber=mRB, loadVector=[0,0,100*2*pi])) 
 85
 86if False: #OPTIONAL: lock rigid body motion of reference frame (for tests):
 87    mbs.AddObject(GenericJoint(markerNumbers=[mGround, mRB], constrainedAxes=[1,1,1, 1,1,0]))
 88
 89#++++++++++++++++++++++++++++++++++++++++++
 90#find nodes at left and right surface:
 91nodeListLeft = fem.GetNodesInPlane(pLeft, [0,0,1])
 92nodeListRight = fem.GetNodesInPlane(pRight, [0,0,1])
 93#nLeft = fem.GetNodeAtPoint(pLeft)
 94#nRight = fem.GetNodeAtPoint(pRight)
 95
 96
 97lenLeft = len(nodeListLeft)
 98lenRight = len(nodeListRight)
 99weightsLeft = np.array((1./lenLeft)*np.ones(lenLeft))
100weightsRight = np.array((1./lenRight)*np.ones(lenRight))
101
102addSupports = True
103if addSupports:
104    k = 2e8     #joint stiffness
105    d = k*0.01  #joint damping
106
107    useSpringDamper = True
108
109    mLeft = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=objFFRF['oFFRFreducedOrder'], 
110                                                    meshNodeNumbers=np.array(nodeListLeft), #these are the meshNodeNumbers
111                                                    weightingFactors=weightsLeft))
112    mRight = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=objFFRF['oFFRFreducedOrder'], 
113                                                    meshNodeNumbers=np.array(nodeListRight), #these are the meshNodeNumbers 
114                                                    weightingFactors=weightsRight))
115    if useSpringDamper:
116        oSJleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mLeft, mGroundPosLeft],
117                                            stiffness=[k,k,k], damping=[d,d,d]))
118        oSJright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mRight,mGroundPosRight],
119                                            stiffness=[k,k,0], damping=[d,d,d]))
120    else:
121        oSJleft = mbs.AddObject(SphericalJoint(markerNumbers=[mGroundPosLeft,mLeft], visualization=VObjectJointSpherical(jointRadius=nodeDrawSize)))
122        oSJright= mbs.AddObject(SphericalJoint(markerNumbers=[mGroundPosRight,mRight], visualization=VObjectJointSpherical(jointRadius=nodeDrawSize)))
123                                                    
124
125#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
126fileDir = 'solution/'
127#rigid body node:
128mbs.AddSensor(SensorNode(nodeNumber=objFFRF['nRigidBody'], storeInternal=True,#fileName=fileDir+'rbDisplacement.txt', 
129                         outputVariableType = exu.OutputVariableType.Displacement))
130
131mbs.AddSensor(SensorNode(nodeNumber=objFFRF['nRigidBody'], storeInternal=True,#fileName=fileDir+'rbVelocity.txt', 
132                         outputVariableType = exu.OutputVariableType.Velocity))
133
134mbs.AddSensor(SensorNode(nodeNumber=objFFRF['nRigidBody'], storeInternal=True,#fileName=fileDir+'rbAcceleration.txt', 
135                         outputVariableType = exu.OutputVariableType.Acceleration))
136
137mbs.AddSensor(SensorNode(nodeNumber=objFFRF['nRigidBody'], storeInternal=True,#fileName=fileDir+'nRigidBodyAngVelCMS.txt', 
138                         outputVariableType = exu.OutputVariableType.AngularVelocity))
139
140
141#FFRF object, selected node:
142sCMSdisp=mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'], meshNodeNumber=nMid, #meshnode number!
143                         storeInternal=True,#fileName=fileDir+'nMidDisplacementCMS.txt', 
144                         outputVariableType = exu.OutputVariableType.Displacement))
145
146sCMSvel=mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'], meshNodeNumber=nMid, #meshnode number!
147                         storeInternal=True,#fileName=fileDir+'nMidVelocityCMS.txt', 
148                         outputVariableType = exu.OutputVariableType.Velocity))
149
150sCMSacc=mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'], meshNodeNumber=nMid, #meshnode number!
151                         storeInternal=True,#fileName=fileDir+'nMidAccelerationCMS.txt', 
152                         outputVariableType = exu.OutputVariableType.Acceleration))
153
154mbs.Assemble()
155
156#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
157simulationSettings = exu.SimulationSettings()
158
159SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
160SC.visualizationSettings.nodes.drawNodesAsPoint = False
161SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
162
163SC.visualizationSettings.nodes.show = True
164SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
165SC.visualizationSettings.nodes.basisSize = 0.12
166SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
167
168SC.visualizationSettings.view0.scene.showFaceEdges = True
169SC.visualizationSettings.view0.scene.showFaces = True
170
171SC.visualizationSettings.sensors.show = True
172SC.visualizationSettings.sensors.drawSimplified = False
173SC.visualizationSettings.sensors.defaultSize = 0.01
174SC.visualizationSettings.markers.drawSimplified = False
175SC.visualizationSettings.markers.show = True
176SC.visualizationSettings.markers.defaultSize = 0.01
177
178SC.visualizationSettings.loads.drawSimplified = False
179
180SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
181SC.visualizationSettings.contour.outputVariableComponent = 1 #y-component
182
183simulationSettings.solution.file.information = "ObjectFFRFreducedOrder test"
184
185h=1e-4
186tEnd = 0.001
187#testIsActive = True
188if not testIsActive:
189    tEnd = 0.1
190    #if not testIsActive:
191#    tEnd = 0.1
192
193simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
194simulationSettings.timeIntegration.endTime = tEnd
195simulationSettings.solution.file.writePeriod = h
196simulationSettings.timeIntegration.verboseMode = 1
197#simulationSettings.timeIntegration.verboseModeFile = 3
198simulationSettings.timeIntegration.newton.useModifiedNewton = True
199
200simulationSettings.solution.sensors.writePeriod = h
201simulationSettings.solution.file.name = "solution/coordinatesSolutionCMStest.txt"
202
203simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5 #SHOULD work with 0.9 as well
204simulationSettings.solution.file.write = False
205#simulationSettings.show.statistics = True
206#simulationSettings.show.computationTime = True
207
208#create animation:
209#simulationSettings.solution.recordImagesInterval = 0.0002
210#SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
211
212if not testIsActive:
213    SC.renderer.Start()
214    SC.renderer.RestoreSavedState() #load last model view
215
216    SC.renderer.DoIdleTasks() #press space to continue
217
218mbs.SolveDynamic(simulationSettings)
219    
220
221data=mbs.GetSensorStoredData(sCMSacc)
222#data = np.loadtxt(fileDir+'nMidAccelerationCMS.txt', comments='#', delimiter=',')
223result = abs(data).sum()
224exu.Print('solution of ObjectFFRFreducedOrderAccelerations=',result)
225
226exu.sys['testResult'] = result/(10*61576.266114362006)
227exu.Print('ObjectFFRFreducedOrderAccelerations test result=',exu.sys['testResult'])
228
229if not testIsActive:
230    #SC.renderer.DoIdleTasks()
231    SC.renderer.Stop() #safely close rendering window!
232
233
234#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
235#plot results
236if not testIsActive:
237        
238    import matplotlib.pyplot as plt
239    import matplotlib.ticker as ticker
240    
241    from exudyn.signalProcessing import FilterSensorOutput, FilterSignal
242    # 
243
244    cList=['r-','g-','b-','k-','c-','r:','g:','b:','k:','c:']
245 
246    # data = np.loadtxt(fileDir+'nMidDisplacementCMS.txt', comments='#', delimiter=',') #new result from this file
247    # # plt.plot(data[:,0], data[:,2], cList[1],label='displ mid') #numerical solution, 1 == x-direction
248
249    # dataV = np.loadtxt(fileDir+'nMidVelocityCMS.txt', comments='#', delimiter=',') #new result from this file
250    # #plt.plot(dataV[:,0], dataV[:,2], cList[0],label='vel mid') 
251
252    # dataA = np.loadtxt(fileDir+'nMidAccelerationCMS.txt', comments='#', delimiter=',') #new result from this file
253    
254    data = mbs.GetSensorStoredData(sCMSdisp)
255    dataV = mbs.GetSensorStoredData(sCMSvel)
256    dataA = mbs.GetSensorStoredData(sCMSacc)
257    plt.plot(dataA[:,0], dataA[:,2], cList[0],label='acc mid') 
258
259    # der = FilterSensorOutput(data, 0, 3, 1)
260    # plt.plot(der[:,0], der[:,2], cList[1],label='vMid,num diff cent') 
261
262    # der = FilterSensorOutput(data, 0, 3, 1, False)
263    # plt.plot(der[:,0], der[:,2], cList[2],label='vMid,num diff left') 
264
265    # der = FilterSensorOutput(data, 5, 3, 1)
266    # plt.plot(der[:,0], der[:,2], cList[3],label='vMid,savgol, window=5, p=3') 
267
268    # der = FilterSensorOutput(data, 7, 5, 1)
269    # plt.plot(der[:,0], der[:,2], cList[4],label='vMid,savgol, window=7, p=5') 
270
271    der = FilterSensorOutput(data, filterWindow=0, polyOrder=3, derivative=2)
272    plt.plot(der[:,0], der[:,2], cList[1],label='diffdiff(displ) mid, direct') 
273
274    #der = FilterSensorOutput(data, filterWindow=5, polyOrder=3, derivative=2)
275    #plt.plot(der[:,0], der[:,2], cList[3],label='diffdiff(displ) mid, savgol, w=5, p=3') 
276    der0 = data[:,0]
277    der2 = FilterSignal(data[:,2], samplingRate=data[1,0]-data[0,0], filterWindow=5, polyOrder=3, derivative=2)
278    plt.plot(der0, der2, cList[3],label='diffdiff(displ) mid, savgol, w=5, p=3') 
279
280    ax=plt.gca() # get current axes
281    ax.grid(True, 'major', 'both')
282    ax.xaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
283    ax.yaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
284    plt.tight_layout()
285    plt.legend()
286
287    plt.figure()
288    # data = np.loadtxt(fileDir+'nMidDisplacementCMS.txt', comments='#', delimiter=',') #new result from this file
289    # # plt.plot(data[:,0], data[:,2], cList[1],label='uMid,Test') #numerical solution, 1 == x-direction
290
291    # dataV = np.loadtxt(fileDir+'nMidVelocityCMS.txt', comments='#', delimiter=',') #new result from this file
292    # #plt.plot(dataV[:,0], dataV[:,2], cList[0],label='vMid,Test') 
293
294    # dataA = np.loadtxt(fileDir+'nMidAccelerationCMS.txt', comments='#', delimiter=',') #new result from this file
295    plt.plot(dataA[:,0], dataA[:,2], cList[0],label='rigid node, acc') 
296
297    der = FilterSensorOutput(dataV, filterWindow=5, polyOrder=3, derivative=1)
298    plt.plot(der[:,0], der[:,2], cList[1],label='rigid node, diff(vel)') 
299
300    der = FilterSensorOutput(data, filterWindow=5, polyOrder=3, derivative=2)
301    plt.plot(der[:,0], der[:,2], cList[3],label='rigid node, diffdiff(displ)') 
302
303    ax=plt.gca() # get current axes
304    ax.grid(True, 'major', 'both')
305    ax.xaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
306    ax.yaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
307    plt.tight_layout()
308    plt.legend()
309
310    plt.show()