computeODE2AEeigenvaluesTest.py
You can view and download this file on Github: computeODE2AEeigenvaluesTest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test for computation of eigenvalues with ODE2 equations + algebraic joint constraints
5#
6# Author: Michael Pieber, Johannes Gerstmayr
7# Date: 2023-06-08
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
16
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
30#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31#rotating rigid body:
32SC = exudyn.SystemContainer()
33mbs = SC.AddSystem()
34
35beamL=0.1 #in m
36beamW=0.01
37beamH=0.001
38rho=5000 #kg/m**3
39springL=0.02 #in m
40springK=1e1 #in N/m
41
42oGround = mbs.AddObject(ObjectGround())
43
44inertiaCuboid=InertiaCuboid(density=rho,
45 sideLengths=[beamL,beamH,beamW])
46
47bBeam = mbs.CreateRigidBody(inertia = inertiaCuboid,
48 referencePosition = [beamL*0.5,0,0],
49 gravity = [0,-9.81*0,0],
50 graphicsDataList = [GraphicsDataOrthoCubePoint(size=[beamL,beamH,beamW],
51 color=color4orange)])
52mBeamRight = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bBeam, localPosition=[beamL*0.5,0,0]))
53
54mbs.CreateGenericJoint(bodyNumbers= [oGround,bBeam], position= [0.,0.,0.],
55 rotationMatrixAxes= np.eye(3), constrainedAxes= [1,1,1,1,1,0],
56 axesRadius=0.001, axesLength= 0.01, color= color4default)
57
58markerToConnect = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[beamL,-springL,0]))
59
60mbs.AddObject(CartesianSpringDamper(markerNumbers=[markerToConnect,mBeamRight],
61 stiffness=[0,springK,0],
62 damping=[0,0,0],
63 offset=[0,springL,0],
64 visualization=VObjectConnectorCartesianSpringDamper(show=True,drawSize=0.01)
65 ))
66mbs.Assemble()
67[eigenValues, eVectors] = mbs.ComputeODE2Eigenvalues()
68
69evNumerical = np.sqrt(eigenValues[0]) / (2*np.pi)
70
71thetaZZ=inertiaCuboid.Translated([-beamL/2,0,0]).Inertia()[2,2]
72evAnalytical = np.sqrt( springK*beamL**2/thetaZZ ) / (2*np.pi)
73
74u = (evAnalytical-evNumerical)/evAnalytical
75exu.Print('numerical eigenvalues in Hz:',evNumerical)
76exu.Print('analytical eigenvalues in Hz:',evAnalytical)
77exu.Print('error eigenvalues:', u)
78
79#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
80#mechanism
81SC = exudyn.SystemContainer()
82mbs = SC.AddSystem()
83
84beamL=0.1 #in m
85beamW=0.01
86beamH=0.001
87rho=5000 #kg/m**3
88springK=1e3 #in N/m
89
90oGround = mbs.AddObject(ObjectGround())
91
92inertiaCuboid=InertiaCuboid(density=rho,
93 sideLengths=[beamL,beamH,beamW])
94
95p0 = np.array([beamL*0.5,0,0])
96b0 = mbs.CreateRigidBody(inertia = inertiaCuboid,
97 referencePosition = p0,
98 gravity = [0,-9.81,0],
99 graphicsDataList = [GraphicsDataOrthoCubePoint(size=[beamL,beamH,beamW],
100 color=color4orange)])
101
102R1 = RotationMatrixZ(-0.25*pi)@RotationMatrixY(0.25*pi)
103p1 = 2*p0 + R1@p0
104b1 = mbs.CreateRigidBody(inertia = inertiaCuboid,
105 referencePosition = p1,
106 referenceRotationMatrix = R1,
107 gravity = [0,-9.81,0],
108 graphicsDataList = [GraphicsDataOrthoCubePoint(size=[beamL,beamH,beamW],
109 color=color4dodgerblue)])
110
111mbs.CreateGenericJoint(bodyNumbers= [oGround,b0], position= [0.,0.,0.],
112 constrainedAxes= [1,1,1,1,1,0],
113 axesRadius=beamH*2, axesLength=beamW*1.05)
114
115mB0 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=b0, localPosition=p0))
116mB1 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=b1, localPosition=-p0))
117
118mbs.AddObject(GenericJoint(markerNumbers=[mB1,mB0], constrainedAxes=[1,1,1, 1,0,0],
119 rotationMarker0=np.eye(3),
120 rotationMarker1=np.eye(3),
121 # rotationMarker1=R1.T,
122 visualization=VGenericJoint(axesRadius=beamH*2, axesLength=beamW*1.05)))
123
124mbs.CreateCartesianSpringDamper(bodyOrNodeList=[b1, oGround],
125 localPosition0=p0,
126 localPosition1=2*p0 + R1@(2*p0),
127 stiffness=[springK]*3,
128 damping=[springK*1e-5]*3,
129 drawSize = beamW
130 )
131# mbs.CreateGenericJoint(bodyNumbers= [b0, b1], position= 2*p0,
132# constrainedAxes= [1,1,1,1,0,0],
133# axesRadius=beamH, axesLength=beamW)
134
135sPos = mbs.AddSensor(SensorBody(bodyNumber=b1, localPosition=p0,
136 storeInternal=True,
137 outputVariableType=exu.OutputVariableType.Displacement
138 ) )
139
140mbs.Assemble()
141SC.visualizationSettings.loads.show=False
142SC.visualizationSettings.openGL.multiSampling=4
143simulationSettings = exu.SimulationSettings()
144simulationSettings.solutionSettings.sensorsWritePeriod = 1e-3
145simulationSettings.timeIntegration.numberOfSteps=1000
146
147[eigenValues, eVectors] = mbs.ComputeODE2Eigenvalues()
148evNumerical = np.sqrt(eigenValues) / (2*np.pi)
149exu.Print('numerical eigenvalues in Hz:',evNumerical)
150
151if useGraphics:
152 mbs.SolveDynamic(simulationSettings=simulationSettings)
153 mbs.PlotSensor(sPos)
154 period=0.521/20 #measured 20 peaks of oscillation in plot sensor
155 f = 1./period
156 exu.Print('frequency simulated=',f)
157
158 # mbs.SolutionViewer()
159
160u += evNumerical[0]/100
161exu.Print('result of computeODE2AEeigenvaluesTest:', u)
162
163exudynTestGlobals.testError = u - 0.38811732950413347 #should be zero
164exudynTestGlobals.testResult = u