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))