MGVehicleSkidKaneEmbedded.html   (MotionGenesis input/output).
   (1) % MotionGenesis file:  MGVehicleSkidKaneEmbedded.txt
   (2) % Copyright (c) 2009 Motion Genesis LLC.  All rights reserved.
   (3) % Problem: Analysis of Skidding Vehicle.
   (4) %--------------------------------------------------------------------
   (5) %       Physical declarations.
   (6) NewtonianFrame N        % Newtonian reference frame (ground).
   (7) RigidBody      A        % Vehicle chassis (Ao is center of axle).
   (8) RigidBody      B, C     % Vehicle's rolling wheels.
   (9) Point          BN( B )  % Point of wheel B in contact with N.
   (10) Point          CN( C )  % Point of wheel C in contact with N.
   (11) %--------------------------------------------------------------------
   (12) %       Mathematical declarations (constants and variables).
   (13) Constant  b = 0.75 m    % Distance between Ao and Bcm (or Ao and Ccm).
   (14) Constant  R = 0.35 m    % Wheel radii.
   (15) Constant  a = 1.64 m    % Distance from Ao to Acm (A's mass center).
   (16) Constant  g = 9.8 m/s^2 % Earth's gravitational acceleration.
   (17) Variable  v'            % Ax> measure of Ao's velocity in N.
   (18) Variable  wB'           % Ay> measure of B's angular velocity in A.
   (19) Variable  wC'           % Ay> measure of C's angular velocity in A.
   (20) Variable  theta''       % Vehicle heading angle (car yaw).
   (21) Variable  FBx, FBy      % Rolling contact forces on wheel B (Ax>, Ay> directions).
   (22) Variable  FCx, FCy      % Rolling contact forces on wheel C (Ax>, Ay> directions).
   (23) SetGeneralizedSpeed( v, theta' )
   (24) %--------------------------------------------------------------------
   (25) %       Mass and inertia.
   (26) A.SetMass( mA = 640 kg )
   (27) B.SetMass( m = 30 kg )
   (28) C.SetMass( m )
   (29) A.SetInertia( Acm,    IAxx, IAyy, IAzz = 166.6 kg*m^2)
   (30) B.SetInertia( Bcm, A, K, J = 2.0 kg*m^2, K = 1.0 kg*m^2 )
   (31) C.SetInertia( Ccm, A, K, J, K  )
   (32) %--------------------------------------------------------------------
   (33) %       Rotational kinematics.
   (34) A.RotateZ( N, theta )
-> (35) A_N = [cos(theta), sin(theta), 0;  -sin(theta), cos(theta), 0;  0, 0, 1]
-> (36) w_A_N> = theta'*Az>
-> (37) alf_A_N> = theta''*Az>

   (38) B.SetAngularVelocityAcceleration( A, wB*Ay> )
-> (39) w_B_A> = wB*Ay>
-> (40) w_B_N> = wB*Ay> + theta'*Az>
-> (41) alf_B_A> = wB'*Ay>
-> (42) alf_B_N> = -wB*theta'*Ax> + wB'*Ay> + theta''*Az>

   (43) C.SetAngularVelocityAcceleration( A, wC*Ay> )
-> (44) w_C_A> = wC*Ay>
-> (45) w_C_N> = wC*Ay> + theta'*Az>
-> (46) alf_C_A> = wC'*Ay>
-> (47) alf_C_N> = -wC*theta'*Ax> + wC'*Ay> + theta''*Az>

   (48) %--------------------------------------------------------------------
   (49) %       Translational kinematics.
   (50) Ao.SetVelocityAcceleration( N, v*Ax> )
-> (51) v_Ao_N> = v*Ax>
-> (52) a_Ao_N> = v'*Ax> + theta'*v*Ay>

   (53) Acm.Translate( Ao,  a*Ax> )
-> (54) p_Ao_Acm> = a*Ax>
-> (55) v_Acm_N> = v*Ax> + a*theta'*Ay>
-> (56) a_Acm_N> = (v'-a*theta'^2)*Ax> + (theta'*v+a*theta'')*Ay>

   (57) Bcm.Translate( Ao, -b*Ay> )
-> (58) p_Ao_Bcm> = -b*Ay>
-> (59) v_Bcm_N> = (v+b*theta')*Ax>
-> (60) a_Bcm_N> = (v'+b*theta'')*Ax> + theta'*(v+b*theta')*Ay>

   (61) Ccm.Translate( Ao,  b*Ay> )
-> (62) p_Ao_Ccm> = b*Ay>
-> (63) v_Ccm_N> = (v-b*theta')*Ax>
-> (64) a_Ccm_N> = (v'-b*theta'')*Ax> + theta'*(v-b*theta')*Ay>

   (65) BN.SetPositionVelocity(  Bcm,  -R*Az> )
-> (66) p_Bcm_BN> = -R*Az>
-> (67) v_BN_N> = (v+b*theta'-R*wB)*Ax>

   (68) CN.SetPositionVelocity(  Ccm,  -R*Az> )
-> (69) p_Ccm_CN> = -R*Az>
-> (70) v_CN_N> = (v-R*wC-b*theta')*Ax>

   (71) %--------------------------------------------------------------------
   (72) %       Rolling constraints (relates wB and wC to v and theta').
   (73) RollingConstraint[1] = Dot( BN.GetVelocity(N), Ax> )
-> (74) RollingConstraint[1] = v + b*theta' - R*wB

   (75) RollingConstraint[2] = Dot( CN.GetVelocity(N), Ax> )
-> (76) RollingConstraint[2] = v - R*wC - b*theta'

   (77) SolveDt( RollingConstraint = 0,  wB, wC )
-> (78) wB = (v+b*theta')/R
-> (79) wC = (v-b*theta')/R
-> (80) wB' = (v'+b*theta'')/R
-> (81) wC' = (v'-b*theta'')/R

   (82) %--------------------------------------------------------------------
   (83) %       Add forces (although none of these appear in DynamicEquations).
   (84) %       Note: Gravity and normal contact forces also do not appear in dynamics.
   (85) System.AddForceGravity( -g*Az> )
-> (86) Force_Acm> = -mA*g*Az>
-> (87) Force_Bcm> = -m*g*Az>
-> (88) Force_Ccm> = -m*g*Az>

   (89) BN.AddForce( FBx*Ax> + FBy*Ay> )
-> (90) Force_BN> = FBx*Ax> + FBy*Ay>

   (91) CN.AddForce( FCx*Ax> + FCy*Ay> )
-> (92) Force_CN> = FCx*Ax> + FCy*Ay>

   (93) %--------------------------------------------------------------------
   (94) %       Form Kane's equations of motion.
   (95) DynamicEquations = System.GetDynamicsKane()
-> (96) DynamicEquations[1] = (mA+2*m+2*J/R^2)*v' - mA*a*theta'^2
-> (97) DynamicEquations[2] = mA*a*theta'*v + (IAzz+2*K+mA*a^2+2*m*b^2+2*J*b^2/R^2)
        *theta''

   (98) %--------------------------------------------------------------------
   (99) %       Differential equations governing x and y.
   (100) Variable x' = Dot( Ao.GetVelocity(N),  Nx> )
-> (101) x' = cos(theta)*v

   (102) Variable y' = Dot( Ao.GetVelocity(N),  Ny> )
-> (103) y' = sin(theta)*v

   (104) %--------------------------------------------------------------------
   (105) %       Integration parameters and initial values.
   (106) Input  tFinal = 1 sec,  tStep = 0.01 sec,  absError = 1.0E-07
   (107) Input  theta = 0 deg,  x = 0 m,  y = 0 m
   (108) Input  theta' = 0.01 rad/sec           % Perturbation in angular rate.
   (109) Input  v = -25 m/s                     % Approximately 55 miles/hour.
   (110) %--------------------------------------------------------------------
   (111) %       List output quantities and solve ODEs.
   (112) Output  t sec,  x m,  y m,  v m/sec,  theta deg,  theta' rad/sec
   (113) ODE( DynamicEquations,  v',  theta'' )  MGVehicleSkidKaneEmbedded

   (114) % Plot MGVehicleSkidKaneEmbedded.1
   (115) %--------------------------------------------------------------------
   (116) %       At times it is helpful, particularly for control system design
   (117) %       and stability analysis, to produce linearized equations of motion.
   (118) %       For this system, useful information about the stability of certain
   (119) %       motions can be obtained be linearizing the equations of motion
   (120) %       in dv, dtheta', perturbations of v and theta', about a nominal (possible)
   (121) %       solution. One nominal solution is v=Vo (a constant) and theta'=0.
   (122) Constant    Vo               % Forward Speed of Cart v=Vo
   (123) Variable    dv', dw'         % Perturbations of v, theta'
   (124) %--------------------------------------------------------------------
   (125) %       Linearize equations of motion about the nominal solution
   (126) LinearizedDynamics = Taylor( DynamicEquations, 1, v=Vo:dv, v'=0:dv', theta'=0:dw, theta''=0:dw')
-> (127) LinearizedDynamics[1] = (mA+2*m+2*J/R^2)*dv'
-> (128) LinearizedDynamics[2] = mA*a*Vo*dw + (IAzz+2*K+mA*a^2+2*m*b^2+2*J*b^2/R^2)*dw'

   (129) %--------------------------------------------------------------------
   (130) %       Solve for dv' and dw' and  factor the constants in the resulting expression.
   (131) Solve( LinearizedDynamics = 0,  dv',  dw' )
-> (132) dv' = 0
-> (133) dw' = -mA*a*Vo*dw/(IAzz+2*K+mA*a^2+2*m*b^2+2*J*b^2/R^2)

   (134) Zee( dw',  dw, Vo )
-> (135) z1 = mA*a/(IAzz+2*K+mA*a^2+2*m*b^2+2*J*b^2/R^2)
-> (136) dw' = -Vo*z1*dw

   (137) %--------------------------------------------------------------------
   (138) %       Closed-form solutions for linearized differential equations.
   (139) Constant    dv0, dw0         % Initial values of dv and dw
   (140) dv = dv0
-> (141) dv = dv0

   (142) dw = dw0*exp(-z1*Vo*t)
-> (143) dw = dw0*exp(-Vo*z1*t)

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