MGDiskRollingDownInclinedPlaneDynamics.html   (MotionGenesis input/output).
   (1) % MotionGenesis file:  MGDiskRollingDownInclinedPlaneDynamics.txt
   (2) % Copyright (c) 2020 Motion Genesis LLC.
   (3) %--------------------------------------------------------------------
   (4) NewtonianFrame  N       % Inclined plane with Nx> downhill and Ny> perpendicular to the plane.
   (5) RigidFrame      E       % Earth with Ex> horizontally-right, Ey> vertically-upward.
   (6) RigidBody       B       % Cylindrical disk whose symmetry axis parallel to Nz>.
   (7) Point           BN( B ) % Point of B in contact with inclined plane N.
   (8) %--------------------------------------------------------------------
   (9) Constant  r             % Cylinder's radius.
   (10) Constant  g             % Earth's gravitational acceleration.
   (11) Constant  psi           % Angle inclined plane makes with Earth's horizontal.
   (12) Variable  theta''       % B's rotation angle from Ey> with sense -Nz>.
   (13) Variable  x''           % Nx> measure of Bcm's position from No.
   (14) Variable  Fx, Fy        % Friction and normal force on B at BN.
   (15) B.SetMassInertia( m,  Ixx, Iyy, Izz )
   (16) %--------------------------------------------------------------------
   (17) %   Rotational kinematics.
   (18) E.SetRotationMatrixZ( N,  psi )
-> (19) E_N = [cos(psi), sin(psi), 0;  -sin(psi), cos(psi), 0;  0, 0, 1]

   (20) B.RotateNegativeZ( E,  theta )
-> (21) B_E = [cos(theta), -sin(theta), 0;  sin(theta), cos(theta), 0;  0, 0, 1]
-> (22) w_B_E> = -theta'*Bz>
-> (23) alf_B_E> = -theta''*Bz>

   (24) %--------------------------------------------------------------------
   (25) %   Translational kinematics.
   (26) Bcm.Translate( No,  x*Nx> + r*Ny> )
-> (27) p_No_Bcm> = x*Nx> + r*Ny>
-> (28) v_Bcm_N> = x'*Nx>
-> (29) a_Bcm_N> = x''*Nx>

   (30) BN.SetPositionVelocity( Bcm,  -r*Ny> )
-> (31) p_Bcm_BN> = -r*Ny>
-> (32) v_BN_N> = (x'-r*theta')*Nx>

   (33) %--------------------------------------------------------------------
   (34) %   Motion constraints (rolling at BN).
   (35) RollingConstraint = Dot( BN.GetVelocity(N),  Nx> )
-> (36) RollingConstraint = x' - r*theta'

   (37) SolveDt( RollingConstraint = 0,  x' )
-> (38) x' = r*theta'
-> (39) x'' = r*theta''

   (40) %--------------------------------------------------------------------
   (41) %   Add relevant contact/distance forces on body B.
   (42) B.AddForceGravity( -g*Ey> )      % Gravity force on B.
-> (43) Force_Bcm> = -m*g*Ey>

   (44) BN.AddForce( Fx*Nx> + Fy*Ny> )   % Friction and normal force on B at BN.
-> (45) Force_BN> = Fx*Nx> + Fy*Ny>

   (46) %--------------------------------------------------------------------
   (47) %   Method 1: Form equation of motion by taking moments about BN.
   (48) Dynamics = Dot( B.GetDynamics(BN),  -Nz> )
-> (49) Dynamics = Izz*theta'' + m*r*x'' - m*g*r*sin(psi)

   (50) %--------------------------------------------------------------------
   (51) %   Method 2: Form equation of motion with Kane's method.
   (52) SetGeneralizedSpeed( theta' )
   (53) KaneDynamics = System.GetDynamicsKane()
-> (54) KaneDynamics = [(Izz+m*r^2)*theta'' - m*g*r*sin(psi)]

   (55) %--------------------------------------------------------------------
   (56) %   Method 3: Form equations of motion with Newton/Euler method.
   (57) NewtonEuler[1] = Dot(  B.GetDynamics(),     Nx> )
-> (58) NewtonEuler[1] = m*r*theta'' - m*g*sin(psi) - Fx

   (59) NewtonEuler[2] = Dot(  B.GetDynamics(),     Ny> )
-> (60) NewtonEuler[2] = m*g*cos(psi) - Fy

   (61) NewtonEuler[3] = Dot(  B.GetDynamics(Bcm),  Nz> )
-> (62) NewtonEuler[3] = -r*Fx - Izz*theta''

   (63) %--------------------------------------------------------------------
   (64) %   Method 4 (optional): Form equation of motion via power/energy-rate principle.
   (65) Power = System.GetPower()
-> (66) Power = m*g*r*sin(psi)*theta'

   (67) KE = System.GetKineticEnergy()
-> (68) KE = 0.5*(Izz+m*r^2)*theta'^2

   (69) PowerEnergyRatePrinciple = Dt(KE) - Power
-> (70) PowerEnergyRatePrinciple = (Izz+m*r^2)*theta'*theta'' - Power

   (71) DynamicsViaPowerEnergyRate = Explicit( PowerEnergyRatePrinciple ) / theta'
-> (72) DynamicsViaPowerEnergyRate = (Izz+m*r^2)*theta'' - m*g*r*sin(psi)

   (73) %--------------------------------------------------------------------
   (74) %   Record input together with responses.
Saved by Motion Genesis LLC.   Command names and syntax: Copyright (c) 2009-2021 Motion Genesis LLC. All rights reserved.