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