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