MGRollingDiskDynamics.html   (MotionGenesis input/output).
```   (1) % MotionGenesis file:  RollingDiskDynamics.txt
(3) % Problem: Simulation of a rolling disk.
(4) %--------------------------------------------------------------------
(5) NewtonianFrame  A                 % Horizontal plane.
(6) RigidFrame      B, C              % Heading and tilt reference frames.
(7) RigidBody       D                 % Rolling disk.
(8) Point           DA( D )           % Point of D in contact with A.
(9) %--------------------------------------------------------------------
(10) Variable  wx',  wy',  wz'         % Disk D's angular velocity in A.
(11) Variable  x'',  y''               % Locates contact point DA from No.
(12) Variable  qH',  qL',  qS'         % Heading angle, lean angle, spin angle.
(13) Constant  r = 0.343 meters        % Radius of disk (13.5 inches).
(14) Constant  g = 9.8   m/s^2         % Earth's gravitational acceleration.
(15) Constant  b = 0 N*m*s/rad         % Aerodynamic damping on disk.
(16) D.SetMass( m = 2 kg )
(17) D.SetInertia( DCm,  C,  J = 0.25*m*r^2,  I = 2*J,  J )
-> (18) J = 0.25*m*r^2
-> (19) I = 2*J

(20) %--------------------------------------------------------------------
(21) %   Rotational kinematics.
(22) B.RotatePositiveZ( A, qH )
-> (23) B_A = [cos(qH), sin(qH), 0;  -sin(qH), cos(qH), 0;  0, 0, 1]
-> (24) w_B_A> = qH'*Bz>

(25) C.RotateNegativeX( B, qL )
-> (26) C_B = [1, 0, 0;  0, cos(qL), -sin(qL);  0, sin(qL), cos(qL)]
-> (27) w_C_B> = -qL'*Cx>
-> (28) w_C_A> = qH'*Bz> - qL'*Cx>

(29) D.RotatePositiveY( C, qS )
-> (30) D_C = [cos(qS), 0, -sin(qS);  0, 1, 0;  sin(qS), 0, cos(qS)]
-> (31) w_D_C> = qS'*Dy>
-> (32) w_D_A> = -qL'*Cx> + (qS'-sin(qL)*qH')*Cy> + cos(qL)*qH'*Cz>

(33) %---------------------------------------------------------------
(34) %   For efficient dynamics, make a change of variable.
(35) wDA> = D.GetAngularVelocity( A )
-> (36) wDA> = -qL'*Cx> + (qS'-sin(qL)*qH')*Cy> + cos(qL)*qH'*Cz>

(37) ChangeVariables[1] = wx - Dot( wDA>, Cx> )
-> (38) ChangeVariables[1] = wx + qL'

(39) ChangeVariables[2] = wy - Dot( wDA>, Cy> )
-> (40) ChangeVariables[2] = wy + sin(qL)*qH' - qS'

(41) ChangeVariables[3] = wz - Dot( wDA>, Cz> )
-> (42) ChangeVariables[3] = wz - cos(qL)*qH'

(43) Solve( ChangeVariables = 0,  qH', qL', qS' )
-> (44) qH' = wz/cos(qL)
-> (45) qL' = -wx
-> (46) qS' = wy + wz*tan(qL)

(47) %---------------------------------------------------------------
(48) %   Use the simpler version of angular velocity from here on.
(49) D.SetAngularVelocity( A,  wx*Cx> + wy*Cy> + wz*Cz> )
-> (50) w_D_A> = wx*Cx> + wy*Cy> + wz*Cz>

(51) %--------------------------------------------------------------------
(52) %   Form angular acceleration (differentiating in D is more efficient).
(53) alf_D_A> = Dt( D.GetAngularVelocity(A), D )
-> (54) alf_D_A> = (wx'-wz*qS')*Cx> + wy'*Cy> + (wz'+wx*qS')*Cz>

(55) %---------------------------------------------------------------
(56) %   Translational kinematics.
(57) Dcm.Translate( Ao,  x*Ax> + y*Ay> + r*Cz> )
-> (58) p_Ao_Dcm> = x*Ax> + y*Ay> + r*Cz>
-> (59) v_Dcm_A> = x'*Ax> + y'*Ay> - r*sin(qL)*qH'*Cx> + r*qL'*Cy>
-> (60) a_Dcm_A> = x''*Ax> + y''*Ay> - r*sin(qL)*qH'^2*By> - r*(2*cos(qL)^3*qH'
*qL'+sin(qL)*(cos(qL)*wz'+wz*sin(qL)*qL'))/cos(qL)^2*Cx> - r*wx'*Cy> - r*qL'^2*Cz>

(61) DA.SetPositionVelocity( Dcm,  -r*Cz> )
-> (62) p_Dcm_DA> = -r*Cz>
-> (63) v_DA_A> = x'*Ax> + y'*Ay> - r*(wy+sin(qL)*qH')*Cx> + r*(wx+qL')*Cy>

(64) %--------------------------------------------------------------------
(65) %   v_DA_A> is simpler if explicit in wx, wy, wz.
(66) Explicit( v_DA_A>,   wx, wy, wz )
-> (67) v_DA_A> = x'*Ax> + y'*Ay> - r*(wy+wz*tan(qL))*Cx>

(68) %---------------------------------------------------------------
(69) %   Motion constraints (D rolls on A at point DA).
(70) RollingConstraint[1] = Dot( DA.GetVelocity(A), Bx> )
-> (71) RollingConstraint[1] = sin(qH)*y' + cos(qH)*x' - r*(wy+wz*tan(qL))

(72) RollingConstraint[2] = Dot( DA.GetVelocity(A), By> )
-> (73) RollingConstraint[2] = cos(qH)*y' - sin(qH)*x'

(74) SolveDt( RollingConstraint = 0,  x', y' )
-> (75) x' = r*cos(qH)*(wy+wz*tan(qL))
-> (76) y' = r*sin(qH)*(wy+wz*tan(qL))
-> (77) x'' = r*cos(qH)*(wy'+tan(qL)*wz'+wz*qL'/cos(qL)^2) - qH'*y'
-> (78) y'' = qH'*x' + r*sin(qH)*(wy'+tan(qL)*wz'+wz*qL'/cos(qL)^2)

(79) %--------------------------------------------------------------------
(80) %   Motion constraints make v_Dcm_A> and a_Dcm_A> simpler if expressed in C.
(81) v_Dcm_A> := Express( Explicit(v_Dcm_A>, wx, wy, wz),  C )
-> (82) v_Dcm_A> = r*wy*Cx> - r*wx*Cy>

(83) a_Dcm_A> := Express( Explicit(a_Dcm_A>, wx, wy, wz),  C )
-> (84) a_Dcm_A> = r*(wx*wz+wy')*Cx> + r*(wy*wz-wx')*Cy> - r*(wx^2-wy*wz*tan(
qL))*Cz>

(85) %---------------------------------------------------------------
(86) %   Add relevant contact and distance force/torques.
(87) %   Note: Aerodynamic damping affects wx most, wy least.
-> (89) Force_Dcm> = -m*g*Az>

(90) D.AddTorque(  -b*( 2*wx*Cx> + 0.2*wy*Cy> + wz*Cz> ) )
-> (91) Torque_D> = -2*b*wx*Cx> - 0.2*b*wy*Cy> - b*wz*Cz>

(92) %--------------------------------------------------------------------
(93) %   Equations of motion with free-body-diagrams (moments about point DA).
(94) Dynamics[1] = Dot( Cx>,  D.GetDynamics(DA) )
-> (95) Dynamics[1] = 2*b*wx + m*g*r*sin(qL) + J*(wx'-wz*qS') - (I-J)*wy*wz
- m*r^2*(wy*wz-wx')

(96) Dynamics[2] = Dot( Cy>,  D.GetDynamics(DA) )
-> (97) Dynamics[2] = 0.2*b*wy + I*wy' + m*r^2*(wx*wz+wy')

(98) Dynamics[3] = Dot( Cz>,  D.GetDynamics(DA) )
-> (99) Dynamics[3] = b*wz + (I-J)*wx*wy + J*(wz'+wx*qS')

(100) %------------------------------------------------------------
(101) %   Optional: Form and simplify Kane's equations of motion.
(102) SetGeneralizedSpeed( wx, wy, wz )
(103) KaneDynamics = System.GetDynamicsKane()
-> (104) KaneDynamics[1] = m*g*r*sin(qL) + 2*b*wx + (J+m*r^2)*wx' - wz*(m*r^2*
wy+(I-J)*wy+J*qS')
-> (105) KaneDynamics[2] = 0.2*b*wy + m*r^2*wx*wz + (I+m*r^2)*wy'
-> (106) KaneDynamics[3] = b*wz + wx*((I-J)*wy+J*qS') + J*wz'

(107) FactorQuadratic( KaneDynamics,  wx, wy, wz )
-> (108) KaneDynamics[1] = m*g*r*sin(qL) + 2*b*wx + (J+m*r^2)*wx' - J*tan(qL)*wz^2
- (I+m*r^2)*wy*wz
-> (109) KaneDynamics[2] = 0.2*b*wy + m*r^2*wx*wz + (I+m*r^2)*wy'
-> (110) KaneDynamics[3] = b*wz + I*wx*wy + J*tan(qL)*wx*wz + J*wz'

(111) %--------------------------------------------------------------------
(112) %   Solve equations of motion for wx', wy', wz'.
(113) Solve( Dynamics = 0,  wx', wy', wz')
-> (114) wx' = -(m*g*r*sin(qL)+2*b*wx-m*r^2*wy*wz-(I-J)*wy*wz-J*wz*qS')/(J+m*r^2)
-> (115) wy' = -0.2*(b*wy+5*m*r^2*wx*wz)/(I+m*r^2)
-> (116) wz' = -(b*wz+(I-J)*wx*wy)/J - wx*qS'

(117) %--------------------------------------------------------------------
(118) %   Integration parameters and initial values.
(119) Input  tFinal = 12 sec,  tStep = 0.05 sec,  absError = 1.0E-08
(120) Input  x = 0 m,  y = 0 m,  qH = 0 deg,  qL = 10 deg,  qS = 0 deg