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