MGDiskRollingDownInclinedPlaneEnergy.html   (MotionGenesis input/output).
   (1) % MotionGenesis file:  MGDiskRollingDownInclinedPlaneEnergy.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) %--------------------------------------------------------------------
   (24) %   Translational kinematics.
   (25) Bcm.SetPositionVelocity( No,  x*Nx> + r*Ny> )
-> (26) p_No_Bcm> = x*Nx> + r*Ny>
-> (27) v_Bcm_N> = x'*Nx>

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

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

   (35) Solve( RollingConstraint = 0,  x' )
-> (36) x' = r*theta'

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

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

   (43) %--------------------------------------------------------------------
   (44) %   Form B's kinetic energy in N.
   (45) KE = B.GetKineticEnergy()        % Valid when B slides or rolls on N.
-> (46) KE = 0.5*Izz*theta'^2 + 0.5*m*x'^2

   (47) Explicit( KE,  m, r, theta' )    % Valid when B rolls on N.
-> (48) KE = 0.5*(Izz+m*r^2)*theta'^2

   (49) %--------------------------------------------------------------------
   (50) %   Calculate a gravitational potential energy PEgravity.
   (51) PEgravity = System.GetForceGravityPotentialEnergy( -g*Ey>, No )
-> (52) PEgravity = m*g*(r*cos(psi)-sin(psi)*x)

   (53) %--------------------------------------------------------------------
   (54) %   For this system, mechanical energy (kinetic energy plus potential energy) is conserved.
   (55) MechanicalEnergy = Explicit( KE + PEgravity,   m, r, theta' )
-> (56) MechanicalEnergy = m*g*(r*cos(psi)-sin(psi)*x) + 0.5*(Izz+m*r^2)*theta'^2

   (57) %--------------------------------------------------------------------
   (58) %   Determine the initial value of MechanicalEnergy.
   (59) InitialMechanicalEnergy = Evaluate( MechanicalEnergy,  x = 0,  theta' = 0 )
-> (60) InitialMechanicalEnergy = m*g*r*cos(psi)

   (61) %--------------------------------------------------------------------
   (62) %   Use the fact that MechanicalEnergy is constant to solve for theta'.
   (63) eqnToSolve = Explicit( MechanicalEnergy - InitialMechanicalEnergy,  m, r, theta' )
-> (64) eqnToSolve = 0.5*(Izz+m*r^2)*theta'^2 - m*g*sin(psi)*x

   (65) SolveQuadraticPositiveRoot( eqnToSolve = 0,  theta' )
-> (66) theta' = 1.414214*sqrt(m*g*sin(psi)*x/(Izz+m*r^2))

   (67) %--------------------------------------------------------------------
   (68) %   Determine w = theta' for various shapes of B (solid/hollow disk or sphere).
   (69) wSolidDisk = Evaluate( theta',  Izz = 1/2*m*r^2 )
-> (70) wSolidDisk = 1.154701*sqrt(g*sin(psi)*x/r^2)

   (71) wHollowDisk = Evaluate( theta',  Izz = m*r^2 )
-> (72) wHollowDisk = 1*sqrt(g*sin(psi)*x/r^2)

   (73) wSolidSphere = Evaluate( theta',  Izz = 2/5*m*r^2 )
-> (74) wSolidSphere = 1.195229*sqrt(g*sin(psi)*x/r^2)

   (75) wHollowSphere = Evaluate( theta',  Izz = 2/3*m*r^2 )
-> (76) wHollowSphere = 1.095445*sqrt(g*sin(psi)*x/r^2)

   (77) %--------------------------------------------------------------------
   (78) %   Form power on B (consider all forces -- e.g., gravity and normal/friction force).
   (79) %   Note: Due to contact, Fy does not show up in power (normal force is a "workless" force).
   (80) %   Note: When B rolls,   Fx does not show up in power (friction force is "workless" force when rolling).
   (81) PowerGravity = Dot( Bcm.GetResultantForce(), Bcm.GetVelocity(N) )
-> (82) PowerGravity = m*g*sin(psi)*x'

   (83) PowerContact = Dot(  BN.GetResultantForce(),  BN.GetVelocity(N) )
-> (84) PowerContact = Fx*(x'-r*theta')

   (85) Explicit( PowerContact )
-> (86) PowerContact = 0

   (87) Power = Explicit( PowerGravity + PowerContact,  x' )
-> (88) Power = m*g*sin(psi)*x'

   (89) %--------------------------------------------------------------------
   (90) %   Optional: Verify that Power = -Dt( PEgravity ).
   (91) isEqualPower = IsSimplifyEqual( Power,  -Dt(PEgravity) )
-> (92) isEqualPower = true

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