MGRollingDiskDynamics.html   (MotionGenesis input/output).
   (1) % MotionGenesis file:  RollingDiskDynamics.txt
   (2) % Copyright (c) 2018 Motion Genesis LLC.  All rights reserved.
   (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.
   (88) Dcm.AddForce( -m*g*Az> )
-> (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
   (121) Input  wx = 0 rad/sec,  wy = 5 rad/sec,  wz = 0 rad/sec
   (122) %---------------------------------------------------------------
   (123) %   List output quantities and solve ODEs.
   (124) Output  t sec,  x meters,  y meters,  qL degrees,  qH degrees
   (125) ODE() RollingDiskDynamics

   (126) %--------------------------------------------------------------------
Saved by Motion Genesis LLC.   Command names and syntax: Copyright (c) 2009-2021 Motion Genesis LLC. All rights reserved.