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