MGGyrostatSpinStability.html  (MotionGenesis input/output).
   (1) % MotionGenesis file:  MGGyrostatSpinStability.txt 
   (2) % Copyright (c) 2009 Motion Genesis LLC.  All rights reserved.
   (3) %--------------------------------------------------------------------
   (4) NewtonianFrame  N          % Newtonian reference frame
   (5) RigidBody       A          % Carrier
   (6) RigidFrame      B          % Rotor
   (7) %--------------------------------------------------------------------
   (8) Constant  Omega = -20 rad/sec  % B's' angular speed  in A
   (9) Constant  J = 0.07634 kg*m^2   % B's moment of inertia about spin axis
   (10) Variable  wx',  wy',  wz' 
   (11) %--------------------------------------------------------------------
   (12) %       Mass and inertia   (attribute G's inertia to A)
   (13) A.SetMassInertia( m,  Ix = 1.25 kg*m^2,  Iy = 4.25 kg*m^2,  Iz = 5 kg*m^2 ) 
   (14) %--------------------------------------------------------------------
   (15) %       Rotational and translational kinematics.
   (16) %       Note: Velocity/acceleration of G's mass center is 0>.
   (17) A.SetAngularVelocityAcceleration( N, wx*Ax> + wy*Ay> + wz*Az> )
-> (18) w_A_N> = wx*Ax> + wy*Ay> + wz*Az>
-> (19) alf_A_N> = wx'*Ax> + wy'*Ay> + wz'*Az>

   (20) B.SetAngularVelocityAcceleration( A, Omega*Az> )
-> (21) w_B_A> = Omega*Az>
-> (22) w_B_N> = wx*Ax> + wy*Ay> + (Omega+wz)*Az>
-> (23) alf_B_A> = 0>
-> (24) alf_B_N> = (Omega*wy+wx')*Ax> + (wy'-Omega*wx)*Ay> + wz'*Az>

   (25) Acm.SetVelocityAcceleration(  N,  0>  )
-> (26) v_Acm_N> = 0>
-> (27) a_Acm_N> = 0>

   (28) %--------------------------------------------------------------------
   (29) %       Form Kane's equations of motion.
   (30) SetGeneralizedSpeed( wx, wy, wz )
   (31) Dynamics = System.GetGeneralizedForce() + Frstar() + Gyrostat(FrStar,CYLINDER,A,B,J)
-> (32) Dynamics[1] = (Iy-Iz)*wy*wz - J*Omega*wy - Ix*wx'
-> (33) Dynamics[2] = J*Omega*wx - (Ix-Iz)*wx*wz - Iy*wy'
-> (34) Dynamics[3] = (Ix-Iy)*wx*wy - Iz*wz'

   (35) %--------------------------------------------------------------------
   (36) %       Calculate gyrostat's angular momentum .
   (37) H> = System.GetAngularMomentum(Acm) + Gyrostat(Angmom,CYLINDER,A,B,J)
-> (38) H> = Ix*wx*Ax> + Iy*wy*Ay> + (J*Omega+Iz*wz)*Az>

   (39) H = GetMagnitude( H> )
-> (40) H = sqrt(Ix^2*wx^2+Iy^2*wy^2+(J*Omega+Iz*wz)^2)

   (41) phi = AngleBetweenUnitVectors(  GetUnitVector( H> ),  Az> )
-> (42) phi = acos((J*Omega+Iz*wz)/sqrt(Ix^2*wx^2+Iy^2*wy^2+(J*Omega+Iz*wz)^2))

   (43) %--------------------------------------------------------------------
   (44) %       Integration parameters and initial values.
   (45) Input  tFinal = 20 sec,  tStep = 0.1 sec,  absError = 1.0E-07
   (46) Input  wx = 0 rad/sec,  wy = 0.02 rad/sec,  wz = 1 rad/sec
   (47) %--------------------------------------------------------------------
   (48) %       List output quantities and solve ODEs.
   (49) Output  t sec,  phi degs,  H kg*m^2/s
   (50) ODE( Dynamics,  wx', wy', wz' )  MGGyrostatSpinStability

   (51) %********************************************************************
   (52) %       STABILITY ANALYSIS
   (53) %********************************************************************
   (54) %       Linearization: Perturbation vars + nominal solution parameters.
   (55) Variable  dwx',  dwy',  dwz'      % Perturbations of wx, wy, wz.
   (56) Constant  nwz = 1 rad/sec         % Nominal solution for wz.
   (57) %--------------------------------------------------------------------
   (58) %       Check nominal solution satisifies the equations of motion.
   (59) Check = Evaluate( Dynamics,  wx=0, wx'=0,  wy=0, wy'=0,  wz=nwz, wz'=0 )
-> (60) Check = [0;  0;  0]

   (61) %--------------------------------------------------------------------
   (62) %       Linearize equations of motion about nominal solution.
   (63) Perturb = Linearize1( Dynamics, wx = 0 : dwx,  wx' = 0 : dwx',  wy = 0 : dwy, &
                            wy' =0 : dwy', wz = nwz: dwz,  wz' = 0 : dwz' )
-> (64) Perturb[1] = -(J*Omega-(Iy-Iz)*nwz)*dwy - Ix*dwx'
-> (65) Perturb[2] = (J*Omega-(Ix-Iz)*nwz)*dwx - Iy*dwy'
-> (66) Perturb[3] = -Iz*dwz'

   (67) Solve( Perturb,  dwx', dwy', dwz' )
-> (68) dwx' = -(J*Omega-(Iy-Iz)*nwz)*dwy/Ix
-> (69) dwy' = (J*Omega-(Ix-Iz)*nwz)*dwx/Iy
-> (70) dwz' = 0

   (71) %--------------------------------------------------------------------
   (72) %       Form, X, X', and A matrices in the matrix equation X' = A * x
   (73) Xm = [ dwx;  dwy;  dwz ]
-> (74) Xm = [dwx;  dwy;  dwz]

   (75) Xp = Dt( Xm )
-> (76) Xp = [dwx';  dwy';  dwz']

   (77) Am = D(  Xp,  Transpose(Xm)  )
-> (78) Am[1,1] = 0
-> (79) Am[1,2] = -(J*Omega-(Iy-Iz)*nwz)/Ix
-> (80) Am[1,3] = 0
-> (81) Am[2,1] = (J*Omega-(Ix-Iz)*nwz)/Iy
-> (82) Am[2,2] = 0
-> (83) Am[2,3] = 0
-> (84) Am[3,1] = 0
-> (85) Am[3,2] = 0
-> (86) Am[3,3] = 0

   (87) %--------------------------------------------------------------------
   (88) %       To find eigenvalues of Am symbolically, find the roots of the
   (89) %       equation found by setting determinant( Lambda * I - A ) = 0.
   (90) Variable Lambda
   (91) det = Determinant( Lambda * GetIdentityMatrix(3) - Am )
-> (92) det = Lambda*((J*Omega-(Ix-Iz)*nwz)*(J*Omega-(Iy-Iz)*nwz)/(Ix*Iy)+Lamb
        da^2)

   (93) det /= Lambda                     % Inspection of det shows Lambda = 0 is a root
-> (94) det = (J*Omega-(Ix-Iz)*nwz)*(J*Omega-(Iy-Iz)*nwz)/(Ix*Iy) + Lambda^2

   (95) %--------------------------------------------------------------------
   (96) %       Find values of Omega which result in Lambda > 0.
   (97) det := EvaluateAtInput( det,  Omega = Omega )
-> (98) det = 0.001096997*(9.824469+Omega)*(49.12235+Omega) + Lambda^2

   (99) %--------------------------------------------------------------------
Saved by Motion Genesis LLC.   Portions copyright (c) 2009-2017 Motion Genesis LLC. Rights reserved. Only for use with MotionGenesis.