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