gyroStability.py

You can view and download this file on Github: gyroStability.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Example showing stable and unstable rotation axes of gyros
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2022-01-19
  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.itemInterface import *
 15from exudyn.utilities import * #includes graphics and rigid body utilities
 16import numpy as np
 17
 18SC = exu.SystemContainer()
 19mbs = SC.AddSystem()
 20
 21
 22#%%++++++++++++++++++++++++++++++++++++++++++++++++++++
 23#physical parameters
 24g = [0,-9.81*0,0]           #gravity
 25L = 0.1                     #length
 26w = 0.03                    #width
 27r = 0.5*w
 28
 29#origin of bodies:
 30pList =    [[0,0,0], [3*L,0,0], [6*L,0,0], ]
 31#initial angular velocities of bodies:
 32angVel = 5
 33eps = 2e-4
 34omegaList =    [[angVel,eps*angVel,eps*angVel],
 35                [eps*angVel,angVel,eps*angVel],
 36                [eps*angVel,eps*angVel,angVel]]
 37sAnglVelList= []
 38sAngleList  = []
 39nodeList = []
 40#ground body
 41oGround = mbs.AddObject(ObjectGround())
 42
 43doImplicit = True
 44#create several bodies
 45for i, p0 in enumerate(pList):
 46    omega0 = omegaList[i]
 47    #%%++++++++++++++++++++++++++++++++++++++++++++++++++++
 48    iCyl0 = InertiaCylinder(density=1000, length=2*L, outerRadius=r, axis=0, innerRadius=r*0.5)
 49    iCyl1 = InertiaCylinder(density=1000, length=L, outerRadius=r, axis=1, innerRadius=r*0.5)
 50    iCylSum = iCyl1.Translated([0,0.5*L,0]) + iCyl0
 51    com = iCylSum.COM()
 52    iCylSum = iCylSum.Translated(-com)
 53
 54    # print(iCyl0)
 55    # print(iCyl1)
 56    # print(iCylSum)
 57
 58    #graphics for body
 59    gCyl0 = GraphicsDataSolidOfRevolution (pAxis=[0,0,0]-com, vAxis=[2*L,0,0],
 60                                           contour = [[-L,-r],[L,-r],[L,-0.5*r],[-L,-0.5*r],],
 61                                           color= color4steelblue, nTiles= 32, smoothContour=False)
 62    gCyl1 = GraphicsDataCylinder(pAxis=[-L,0,0]-com, vAxis=[2*L,0,0],radius=0.5*w, color=color4steelblue, nTiles=32)
 63    gCyl2 = GraphicsDataCylinder(pAxis=[0,0,0]-com, vAxis=[0,L,0],radius=0.5*w, color=color4steelblue, nTiles=32)
 64    gCOM = GraphicsDataBasis(origin=[0,0,0],length = 1.25*L)
 65
 66    nodeType = exu.NodeType.RotationRotationVector
 67    if doImplicit:
 68        nodeType = exu.NodeType.RotationEulerParameters
 69
 70    [n0,b0]=AddRigidBody(mainSys = mbs,
 71                         inertia = iCylSum, #includes COM
 72                         nodeType = nodeType,
 73                         position = p0,
 74                         rotationMatrix = np.diag([1,1,1]),
 75                         angularVelocity = omega0,
 76                         gravity = g,
 77                         graphicsDataList = [gCyl1, gCyl2, gCOM])
 78
 79    sAngVel = mbs.AddSensor(SensorNode(nodeNumber=n0, fileName='solution/gyroAngVel'+str(i)+'.txt',
 80                                       outputVariableType=exu.OutputVariableType.AngularVelocityLocal))
 81    sAngle = mbs.AddSensor(SensorNode(nodeNumber=n0, fileName='solution/gyroAngle'+str(i)+'.txt',
 82                                       outputVariableType=exu.OutputVariableType.Rotation))
 83
 84    nodeList += [n0]
 85    sAnglVelList += [sAngVel]
 86    sAngleList += [sAngle]
 87
 88#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++
 89#assemble system before solving
 90mbs.Assemble()
 91
 92simulationSettings = exu.SimulationSettings() #takes currently set values or default values
 93
 94tEnd = 20 #simulation time
 95h = 0.5e-3 #step size
 96simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
 97simulationSettings.timeIntegration.endTime = tEnd
 98simulationSettings.timeIntegration.verboseMode = 1
 99simulationSettings.solutionSettings.solutionWritePeriod = 0.04 #store every 5 ms
100# simulationSettings.displayComputationTime = True
101
102SC.visualizationSettings.window.renderWindowSize=[1600,1080]
103SC.visualizationSettings.openGL.multiSampling = 4
104SC.visualizationSettings.general.autoFitScene = False
105SC.visualizationSettings.general.graphicsUpdateInterval = 0.02
106SC.visualizationSettings.general.showSolutionInformation = 0
107SC.visualizationSettings.general.showSolverInformation = 0
108# SC.visualizationSettings.general.showSolverTime = 0
109SC.visualizationSettings.general.renderWindowString = 'gyro stability for rotation about axis with \nsmallest (red), middle (green), largest (blue) moment of inertia\ninitial angular velocity = '+str(angVel)+' rad/s, disturbed by '+str(eps*angVel)+' rad/s'
110SC.visualizationSettings.general.textSize = 16
111
112SC.visualizationSettings.nodes.drawNodesAsPoint=False
113SC.visualizationSettings.nodes.showBasis=False
114
115#h=5e-4:
116# omega 0 =  [ 4.99821534  0.17224996 -0.02919272]
117# omega 1 =  [1.14695373 4.853952   0.82419239]
118# omega 2 =  [-0.08742566  0.11224089  4.99987917]
119useGraphics = False
120
121if useGraphics:
122    simulationSettings.timeIntegration.simulateInRealtime = True
123    exu.StartRenderer()
124    if 'renderState' in exu.sys: #reload old view
125        SC.SetRenderState(exu.sys['renderState'])
126
127    mbs.WaitForUserToContinue() #stop before simulating
128
129if doImplicit:
130    mbs.SolveDynamic(simulationSettings = simulationSettings,
131                     solverType=exu.DynamicSolverType.TrapezoidalIndex2)
132else:
133    mbs.SolveDynamic(simulationSettings = simulationSettings,
134                     solverType=exu.DynamicSolverType.RK44)
135
136if useGraphics:
137    SC.WaitForRenderEngineStopFlag() #stop before closing
138    exu.StopRenderer() #safely close rendering window!
139
140
141for i, n in enumerate(nodeList):
142    om = mbs.GetNodeOutput(n,exu.OutputVariableType.AngularVelocityLocal)
143    print('omega '+str(i)+' = ',om)
144
145if not useGraphics and True:
146    sol = LoadSolutionFile('coordinatesSolution.txt')
147
148    mbs.SolutionViewer(sol)
149
150if False:
151
152
153    for i, omega in enumerate(omegaList):
154        if i==1:
155            s='unstable'
156        else:
157            s='stable'
158        mbs.PlotSensor([sAnglVelList[i]]*3,[0,1,2],
159                   yLabel='omega init = '+str(omega),colorCodeOffset=0*i*3,
160                   #fileName='plots/gyroRotAxis'+str(i)+'-'+s+'-Eps'+str(eps*100)+'percent.png',
161                   fontSize=12,
162                   newFigure=True, closeAll=(i==0))