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