complexEigenvaluesTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test for computation of complex eigenvalues with ODE2 equations + algebraic joint constraints
  5#
  6# Author:   Michael Pieber, Johannes Gerstmayr
  7# Date:     2024-05-03
  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 #only import if it does not conflict
 16import numpy as np
 17import sys
 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
 31u = 0
 32#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 33SC = exu.SystemContainer()
 34mbs = SC.AddSystem()
 35
 36L=0.5
 37mass = 1.6          #mass in kg
 38spring = 4000       #stiffness of spring-damper in N/m
 39damper = 8          #damping constant in N/(m/s)
 40
 41omega0 = np.sqrt(spring/mass)     #eigen frequency of undamped system
 42dRel = damper/(2*np.sqrt(spring*mass)) #dimensionless damping
 43omega = omega0*np.sqrt(1-dRel**2) #eigen frequency of damped system
 44#note: (without offset!!!)
 45#decay = q_k+1 / q_k = np.exp(-2*pi/np.sqrt(1/dRel**2-1))
 46#D = 1/np.sqrt(4*pi**2/(np.log(decay)**2)+1)
 47
 48exu.Print('dRel =',dRel, ', decay=',np.exp(-2*pi/np.sqrt(1/dRel**2-1)))
 49
 50u0=-0.08            #initial displacement
 51v0=1                #initial velocity
 52f =80               #force on mass
 53x0=f/spring         #static displacement
 54
 55exu.Print('resonance frequency = '+str(np.sqrt(spring/mass)))
 56exu.Print('static displacement = '+str(x0))
 57
 58#node for 3D mass point:
 59n1=mbs.AddNode(Point(referenceCoordinates = [L,0,0],
 60                     initialCoordinates = [u0,0,0],
 61                     initialVelocities= [v0,0,0]))
 62
 63#ground node
 64nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
 65
 66#add mass point (this is a 3D object with 3 coordinates):
 67massPoint = mbs.AddObject(MassPoint(physicsMass = mass, nodeNumber = n1))
 68
 69#marker for ground (=fixed):
 70groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
 71
 72#marker for node coordinates to be used with spring-damper and constraint:
 73nodeMarker0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 0))
 74nodeMarker1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 1))
 75nodeMarker2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 2))
 76
 77#spring-damper between two marker coordinates
 78oCSD = mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundMarker, nodeMarker0],
 79                                          stiffness = spring, damping = damper))
 80
 81if True: #False: 3 eigen values, True: only one eigenvalue
 82    oCC1 = mbs.AddObject(CoordinateConstraint(markerNumbers=[groundMarker, nodeMarker1]))
 83    oCC2 = mbs.AddObject(CoordinateConstraint(markerNumbers=[groundMarker, nodeMarker2]))
 84
 85#add load:
 86mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker0,
 87                                         load = f))
 88
 89#add sensor:
 90sForce = mbs.AddSensor(SensorObject(objectNumber=oCSD, storeInternal=True,
 91                           outputVariableType=exu.OutputVariableType.Force))
 92sPos = mbs.AddSensor(SensorNode(nodeNumber=n1, storeInternal=True,
 93                           outputVariableType=exu.OutputVariableType.Position))
 94
 95mbs.Assemble()
 96
 97tEnd = 1     #end time of simulation
 98h = 0.001    #step size; leads to 1000 steps
 99
100simulationSettings = exu.SimulationSettings()
101simulationSettings.solutionSettings.solutionWritePeriod = 5e-3 #output interval general
102simulationSettings.solutionSettings.sensorsWritePeriod = 5e-3  #output interval of sensors
103simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h) #must be integer
104simulationSettings.timeIntegration.endTime = tEnd
105
106simulationSettings.timeIntegration.verboseMode = 1             #show some solver output
107# simulationSettings.displayComputationTime = True               #show how fast
108
109simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
110
111#add some drawing parameters for this example
112SC.visualizationSettings.nodes.drawNodesAsPoint=False
113SC.visualizationSettings.nodes.defaultSize=0.1
114
115
116#start solver:
117mbs.SolveDynamic(simulationSettings)
118
119if useGraphics:
120    mbs.PlotSensor(sPos, closeAll=True, title='linear mass-spring-damper')
121
122[eigenValues, eVectors] = mbs.ComputeODE2Eigenvalues(
123                                                  computeComplexEigenvalues=True,
124                                                  useAbsoluteValues=False)
125eigenFreqHz = abs(eigenValues[0].imag) / (2*np.pi)
126eigenDRel = abs(eigenValues[0].real)
127freqAnalytic = omega0*(-dRel+sqrt(1-dRel**2)*1j)
128
129exu.Print('MSD analytical eigenvalues:',freqAnalytic)
130exu.Print('MSD complex eigenvalues:',eigenValues)
131exu.Print('MSD numerical eigenfrequency in Hz:',eigenFreqHz)
132exu.Print('numerical damping dRel           :',eigenDRel/abs(eigenValues[0]))
133exu.Print('*********************\n')
134
135u += eigenFreqHz-omega0/(2*pi) + eigenDRel/abs(eigenValues[0])
136# sys.exit()
137
138#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
139#rotating rigid body:
140SC = exudyn.SystemContainer()
141mbs = SC.AddSystem()
142
143beamL=0.1 #in m
144beamW=0.01
145beamH=0.001
146rho=5000 #kg/m**3
147
148dRel = 0.0001
149springL = 0.02 #in m
150springK = 10   #in N/m
151
152oGround = mbs.AddObject(ObjectGround())
153
154inertiaCuboid=InertiaCuboid(density=rho,
155                        sideLengths=[beamL,beamH,beamW])
156thetaZZ=inertiaCuboid.Translated([-beamL/2,0,0]).Inertia()[2,2]
157
158#damping:
159springD = dRel * (2*np.sqrt(springK*beamL**2/thetaZZ)) #undamped!
160
161bBeam = mbs.CreateRigidBody(inertia = inertiaCuboid,
162                        referencePosition = [beamL*0.5,0,0],
163                        #gravity = [0,-9.81*0,0],
164                        graphicsDataList = [graphics.Brick(size=[beamL,beamH,beamW],
165                        color=graphics.color.orange)])
166mBeamRight = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bBeam, localPosition=[beamL*0.5,0,0]))
167
168mbs.CreateGenericJoint(bodyNumbers= [oGround,bBeam], position= [0.,0.,0.],
169                       rotationMatrixAxes= np.eye(3), constrainedAxes= [1,1,1,1,1,0],
170                       axesRadius=0.001, axesLength= 0.01, color= graphics.color.default)
171
172markerToConnect = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[beamL,-springL,0]))
173
174mbs.AddObject(CartesianSpringDamper(markerNumbers=[markerToConnect,mBeamRight],
175                                    stiffness=[0,springK,0],
176                                    damping=[0,springD,0],
177                                    offset=[0,springL*0.9,0], #this induces damped oscillations
178                                    visualization=VObjectConnectorCartesianSpringDamper(show=True,drawSize=0.01)
179                                    ))
180sPos = mbs.AddSensor(SensorBody(bodyNumber=bBeam, storeInternal=True, localPosition=[beamL*0.5,0,0],
181                           outputVariableType=exu.OutputVariableType.Position))
182
183mbs.Assemble()
184[eigenValues, eVectors] = mbs.ComputeODE2Eigenvalues(computeComplexEigenvalues=True,
185                                                  useAbsoluteValues=False)
186
187evNumerical = np.abs(eigenValues[0]) / (2*np.pi)  #abs(omega0*(-dRel+sqrt(1-dRel**2)*1j) ) gives omega0!!!
188evNumericalDampedImag = np.abs(eigenValues[0].imag) / (2*np.pi) #this gives the damped case!
189evNumericalDampedReal = np.abs(eigenValues[0].real) #this gives the damped case!
190
191#use generalized eigenvalue solver to compare with undamped case
192[eigenValuesGE, eVectorsGE] = mbs.ComputeODE2Eigenvalues(computeComplexEigenvalues=False,
193                                                  useAbsoluteValues=True)
194
195
196evAnalytical = np.sqrt( springK*beamL**2/thetaZZ ) / (2*np.pi)
197
198u += (evAnalytical-abs(evNumerical))/evAnalytical
199exu.Print('numerical eigenfrequency (in Hz) :',evNumerical)
200exu.Print('numerical eigenfrequency damped  :',evNumericalDampedImag)
201exu.Print('numerical damping dRel           :',evNumericalDampedReal/(evNumerical*2*pi))
202exu.Print('numerical eigenfrequency GE      :',np.sqrt(eigenValuesGE[0])/(2*pi))
203exu.Print('analytical eigenfrequency (in Hz):',evAnalytical)
204exu.Print('error eigenvalues:', u, '\n*********************\n')
205
206#decay measured: 1.130mm vs 0.777mm => decay= 1.454
207#      ==> dRel measured  = 0.0595
208#      ==> dRel numerical = 0.05999850003749906
209
210tEnd = 1.5     #end time of simulation
211h = 0.002    #step size; leads to 1000 steps
212
213simulationSettings = exu.SimulationSettings()
214simulationSettings.solutionSettings.sensorsWritePeriod = h  #output interval of sensors
215simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h) #must be integer
216simulationSettings.timeIntegration.endTime = tEnd
217
218mbs.SolveDynamic(simulationSettings)
219mbs.PlotSensor(sPos,components=[1],title='bar with spring at tip')
220
221#sys.exit()
222
223#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
224#mechanism
225SC = exudyn.SystemContainer()
226mbs = SC.AddSystem()
227
228beamL=0.1 #in m
229beamW=0.01
230beamH=0.001
231rho=5000 #kg/m**3
232springK=1e3 #in N/m
233
234oGround = mbs.AddObject(ObjectGround())
235
236inertiaCuboid=InertiaCuboid(density=rho,
237                        sideLengths=[beamL,beamH,beamW])
238
239p0 = np.array([beamL*0.5,0,0])
240b0 = mbs.CreateRigidBody(inertia = inertiaCuboid,
241                         referencePosition = p0,
242                         gravity = [0,-9.81,0],
243                         graphicsDataList = [graphics.Brick(size=[beamL,beamH,beamW],
244                         color=graphics.color.orange)])
245
246R1 = RotationMatrixZ(-0.25*pi)@RotationMatrixY(0.25*pi)
247p1 = 2*p0 + R1@p0
248b1 = mbs.CreateRigidBody(inertia = inertiaCuboid,
249                         referencePosition = p1,
250                         referenceRotationMatrix = R1,
251                         gravity = [0,-9.81,0],
252                         graphicsDataList = [graphics.Brick(size=[beamL,beamH,beamW],
253                         color=graphics.color.dodgerblue)])
254
255mbs.CreateGenericJoint(bodyNumbers= [oGround,b0], position= [0.,0.,0.],
256                       constrainedAxes= [1,1,1,1,1,0],
257                       axesRadius=beamH*2, axesLength=beamW*1.05)
258
259mB0 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=b0, localPosition=p0))
260mB1 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=b1, localPosition=-p0))
261
262mbs.AddObject(GenericJoint(markerNumbers=[mB1,mB0], constrainedAxes=[1,1,1, 1,0,0],
263                           rotationMarker0=np.eye(3),
264                           rotationMarker1=np.eye(3),
265                           # rotationMarker1=R1.T,
266                           visualization=VGenericJoint(axesRadius=beamH*2, axesLength=beamW*1.05)))
267
268mbs.CreateCartesianSpringDamper(bodyOrNodeList=[b1, oGround],
269                                localPosition0=p0,
270                                localPosition1=2*p0 + R1@(2*p0),
271                                stiffness=[springK]*3,
272                                damping=[springK*1e-4]*3,
273                                drawSize = beamW
274                                )
275# mbs.CreateGenericJoint(bodyNumbers= [b0, b1], position= 2*p0,
276#                        constrainedAxes= [1,1,1,1,0,0],
277#                        axesRadius=beamH, axesLength=beamW)
278
279sPos = mbs.AddSensor(SensorBody(bodyNumber=b1, localPosition=p0,
280                                storeInternal=True,
281                                outputVariableType=exu.OutputVariableType.Displacement
282                                ) )
283
284mbs.Assemble()
285SC.visualizationSettings.loads.show=False
286SC.visualizationSettings.openGL.multiSampling=4
287simulationSettings = exu.SimulationSettings()
288simulationSettings.solutionSettings.sensorsWritePeriod = 1e-3
289simulationSettings.timeIntegration.numberOfSteps=1000
290
291#useAbsoluteValues is used to compare with GE, thus having double eigen values
292[eigenValues, eVectors] = mbs.ComputeODE2Eigenvalues(computeComplexEigenvalues=True,
293                                                  useAbsoluteValues=True)
294
295evNumerical = eigenValues / (2*np.pi)
296exu.Print('numerical eigenvalues in Hz:',evNumerical)
297#compute with generalized eigenvalue solver without damping:
298[eigenValues, eVectors] = mbs.ComputeODE2Eigenvalues(computeComplexEigenvalues=False,
299                                                  useAbsoluteValues=True)
300
301evNumerical = np.sqrt(eigenValues) / (2*np.pi)
302exu.Print('numerical eigenvalues GE:',evNumerical)
303
304if useGraphics:
305    mbs.SolveDynamic(simulationSettings=simulationSettings)
306    mbs.PlotSensor(sPos)
307
308u += evNumerical[0]/100
309exu.Print('result of computeODE2AEeigenvaluesTest2:', u)
310
311exudynTestGlobals.testError = u - 0.38811732950413347 #should be zero
312exudynTestGlobals.testResult = u