RattlebackKane.html   (MotionGenesis input/output).
   (1) % MotionGenesis file:  RattlebackKane.txt
   (2) % Copyright (c) 2009 Motion Genesis LLC.  All rights reserved.
   (3) %--------------------------------------------------------------------
   (4) NewtonianFrame  N
   (5) RigidBody       B            % Rattleback.
   (6) RigidFrame      S            % Ellipsoidal surface of B (So is centroid of surface).
   (7) Point           Bn(B)        % Point of B in contact with N.
   (8) %--------------------------------------------------------------------
   (9) Variable  wx', wy', wz'           % Measures of B's angular velocity in N.
   (10) Variable  q1', q2', q3'           % Orientation angles of B in N.
   (11) Variable  sx', sy', sz', Lambda'  % Variables that describe the surface of B.
   (12) Constant  h = 1.0 cm              % Distance between So and Bcm.
   (13) Constant  g = 9.81 m/sec^2        % Gravity.
   (14) B.SetMassInertia( m = 1.0 kg,  Ixx = 17 kg*cm^2, Iyy = 2 kg*cm^2, Izz = 16 kg*cm^2, 0, Iyz = 0.2 kg*cm^2, 0 )
   (15) SetGeneralizedSpeed( wx, wy, wz )
   (16) %--------------------------------------------------------------------
   (17) %   Rotational kinematics.
   (18) B.SetAngularVelocityAcceleration( N,  wx*Bx> + wy*By> + wz*Bz> )
-> (19) w_B_N> = wx*Bx> + wy*By> + wz*Bz>
-> (20) alf_B_N> = wx'*Bx> + wy'*By> + wz'*Bz>

   (21) B.SetRotationMatrixODE( N,  BodyXYZ,  q1, q2, q3 )
-> (22) B_N[1,1] = cos(q2)*cos(q3)
-> (23) B_N[1,2] = sin(q3)*cos(q1) + sin(q1)*sin(q2)*cos(q3)
-> (24) B_N[1,3] = sin(q1)*sin(q3) - sin(q2)*cos(q1)*cos(q3)
-> (25) B_N[2,1] = -sin(q3)*cos(q2)
-> (26) B_N[2,2] = cos(q1)*cos(q3) - sin(q1)*sin(q2)*sin(q3)
-> (27) B_N[2,3] = sin(q1)*cos(q3) + sin(q2)*sin(q3)*cos(q1)
-> (28) B_N[3,1] = sin(q2)
-> (29) B_N[3,2] = -sin(q1)*cos(q2)
-> (30) B_N[3,3] = cos(q1)*cos(q2)
-> (31) q1' = -(sin(q3)*wy-cos(q3)*wx)/cos(q2)
-> (32) q2' = sin(q3)*wx + cos(q3)*wy
-> (33) q3' = wz + tan(q2)*(sin(q3)*wy-cos(q3)*wx)

   (34) SetDt( C11 = N_B[1,1] );  SetDt( C12 = N_B[1,2] );  SetDt( C13 = N_B[1,3] )
-> (35) C11 = cos(q2)*cos(q3)
-> (36) C12 = -sin(q3)*cos(q2)
-> (37) C13 = sin(q2)

   (38) SetDt( C21 = N_B[2,1] );  SetDt( C22 = N_B[2,2] );  SetDt( C23 = N_B[2,3] )
-> (39) C21 = sin(q3)*cos(q1) + sin(q1)*sin(q2)*cos(q3)
-> (40) C22 = cos(q1)*cos(q3) - sin(q1)*sin(q2)*sin(q3)
-> (41) C23 = -sin(q1)*cos(q2)

   (42) SetDt( C31 = N_B[3,1] );  SetDt( C32 = N_B[3,2] );  SetDt( C33 = N_B[3,3] )
-> (43) C31 = sin(q1)*sin(q3) - sin(q2)*cos(q1)*cos(q3)
-> (44) C32 = sin(q1)*cos(q3) + sin(q2)*sin(q3)*cos(q1)
-> (45) C33 = cos(q1)*cos(q2)

   (46) B_N := Transpose( [C11, C12, C13;  C21, C22, C23;  C31, C32, C33] )
-> (47) B_N = [C11, C21, C31;  C12, C22, C32;  C13, C23, C33]

   (48) %--------------------------------------------------------------------
   (49) %   Position vectors.
   (50) So.SetPosition( Bcm,  h*Bx> )
-> (51) p_Bcm_So> = h*Bx>

   (52) Bn.SetPosition( So,   sx*Bx> + sy*By> + sz*Bz> )
-> (53) p_So_Bn> = sx*Bx> + sy*By> + sz*Bz>

   (54) %--------------------------------------------------------------------
   (55) %   Surface Geometry: The outward normal to the surface S is parallel
   (56) %   to the gradient of the function that characterizes S's surface.
   (57) %   And the outward normal to S at Bn is in the -Nx> direction.
   (58) %   Lambda is a convenient scalar quantity (intermediate variable).
   (59) Constant  a+ = 2.0 cm,  b+ = 20 cm,  c+ = 3.0 cm
   (60) Surface = (sx/a)^2 + (sy/b)^2 + (sz/c)^2 - 1
-> (61) Surface = -1 + sx^2/a^2 + sy^2/b^2 + sz^2/c^2

   (62) Gradient> = D(Surface, sx)*Bx> + D(Surface, sy)*By> + D(Surface, sz)*Bz>
-> (63) Gradient> = 2*sx/a^2*Bx> + 2*sy/b^2*By> + 2*sz/c^2*Bz>

   (64) GradientIsNegativeNxDirection> = Gradient> + 2*Lambda*Nx>
-> (65) GradientIsNegativeNxDirection> = 2*sx/a^2*Bx> + 2*sy/b^2*By> + 2*sz/c^2*Bz> + 2*Lambda*Nx>

   (66) GradientEquations = Dot( GradientIsNegativeNxDirection>,  [Bx>; By>; Bz>] )
-> (67) GradientEquations = [2*sx/a^2 + 2*Lambda*C11;  2*sy/b^2 + 2*Lambda*C12;  2*sz/c^2 + 2*Lambda*C13]

   (68) SolveDt( GradientEquations = 0,  sx, sy, sz )
-> (69) sx = -a^2*Lambda*C11
-> (70) sy = -b^2*Lambda*C12
-> (71) sz = -c^2*Lambda*C13
-> (72) sx' = -a^2*(C11*Lambda'-Lambda*(sin(q2)*cos(q3)*q2'+sin(q3)*cos(q2)*q3'))
-> (73) sy' = -b^2*(C12*Lambda'+Lambda*(sin(q2)*sin(q3)*q2'-cos(q2)*cos(q3)*q3'))
-> (74) sz' = -c^2*(C13*Lambda'+Lambda*cos(q2)*q2')

   (75) %--------------------------------------------------------------------
   (76) %   By inspecting the rattleback schematic, it can be seen that when q1 = q2 = q3 = 0,
   (77) %   point Bn is below So so sx must be negative. Since sx = a^2*C11*Lambda and
   (78) %   since C11 = 1 when q1 = q2 = q3 = 0, Lambda must be the positive root.
   (79) FactorQuadratic( Surface, Lambda )
-> (80) Surface = -1 + Lambda^2*(a^2*C11^2+b^2*C12^2+c^2*C13^2)

   (81) SolveQuadraticPositiveRootDt( Surface,  Lambda )
-> (82) Lambda = 1/sqrt(a^2*C11^2+b^2*C12^2+c^2*C13^2)
-> (83) Lambda' = Lambda*(a^2*C11*(sin(q2)*cos(q3)*q2'+sin(q3)*cos(q2)*q3')-c^2
        *cos(q2)*C13*q2'-b^2*C12*(sin(q2)*sin(q3)*q2'-cos(q2)*cos(q3)*q3'))/(a^2
        *C11^2+b^2*C12^2+c^2*C13^2)

   (84) %--------------------------------------------------------------------
   (85) %   Translational kinematics.
   (86) Bn.SetVelocity( N, 0> )           % Automatically imposes motion constraint.
-> (87) v_Bn_N> = 0>

   (88) Bcm.SetVelocity( N, Bn )
-> (89) v_Bcm_N> = (sy*wz-sz*wy)*Bx> + (sz*wx-(h+sx)*wz)*By> + ((h+sx)*wy-sy*
        wx)*Bz>

   (90) %--------------------------------------------------------------------
   (91) %   Relevant forces on B for generalized force calculation.
   (92) System.AddForceGravity( -g*Nx> )
-> (93) Force_Bcm> = -m*g*Nx>

   (94) %--------------------------------------------------------------------
   (95) %   Form Kane's equations of motion.
   (96) Dynamics = System.GetDynamicsKane()
-> (97) Dynamics[1] = wy*(Iyz*wy+Izz*wz) + m*(sy*((h+sx)*wx*wz+wy*(sy*wz-sz*wy)
        +wx*sy'-wy*sx')+sz*(wz*(sy*wz-sz*wy)+wx*sz'-(h+sx)*wx*wy-wz*sx'))
        + (Ixx+m*(sy^2+sz^2))*wx' - m*g*(sy*C13-sz*C12) - wz*(Iyy*wy+Iyz*wz)
        - m*sy*(h+sx)*wy' - m*sz*(h+sx)*wz'

-> (98) Dynamics[2] = wx*(Ixx*wz-Iyz*wy-Izz*wz) + m*(sz*(sy*wx*wy+wz*(sz*wx-(h+
        sx)*wz)+wy*sz'-wz*sy')-(h+sx)*(sy*wy*wz+wx*sy'-wx*(sz*wx-(h+sx)*wz)-wy*
        sx')) + (Iyz-m*sy*sz)*wz' + (Iyy+m*(sz^2+(h+sx)^2))*wy' - m*g*(sz*C11-(
        h+sx)*C13) - m*sy*(h+sx)*wx'

-> (99) Dynamics[3] = m*g*(sy*C11-(h+sx)*C12) + (Iyz-m*sy*sz)*wy' + (Izz+m*(sy^2
        +(h+sx)^2))*wz' - wx*(Ixx*wy-Iyy*wy-Iyz*wz) - m*(sy*(sz*wx*wz+wy*(sy*
        wx-(h+sx)*wy)+wy*sz'-wz*sy')+(h+sx)*(wx*(sy*wx-(h+sx)*wy)+wx*sz'-sz*wy*
        wz-wz*sx')) - m*sz*(h+sx)*wx'

   (100) %--------------------------------------------------------------------
   (101) %   Expressions for quantities to be output by ODE command.
   (102) delta = AngleBetweenUnitVectors( Nx>, Bx> )   % Wobble/inclination angle
-> (103) delta = acos(C11)

   (104) %--------------------------------------------------------------------
   (105) %   Input numerical integration parameters and initial values.
   (106) Input  tFinal = 5 sec, tStep = 0.01 sec,  absError = 1.0E-07,  relError = 1.0E-07
   (107) Input  q1 = 0.0 deg,    q2 = 0.5 deg,      q3 = -0.5 deg
   (108) Input  wx = 5 rad/sec,  wy = 0.0 rad/sec,  wz = 0.0 rad/sec
   (109) %--------------------------------------------------------------------
   (110) %   List output quantities and solve ODEs (or write MATLAB, C, Fortran, ... code).
   (111) OutputPlot  t seconds,  q1 degrees,  delta degrees
   (112) ODE( Dynamics = 0,  wx', wy', wz' )  RattlebackKane.m

   (113) %--------------------------------------------------------------------
   (114) %   Note: With the given initial values and parameters, the simulation
   (115) %   predicts spin reversal and the following values at t = 5.0 seconds.
   (116) %   t (sec)        q1 (degees)    delta (degrees)
   (117) %   5.000000       157.5726       2.813488
   (118) %--------------------------------------------------------------------
Saved by Motion Genesis LLC.   MotionGenesis 6.3 Professional.
Command names and syntax: Copyright (c) 2009-2023 Motion Genesis LLC. All rights reserved.