(1) % MotionGenesis file:  MGStickSlipMassSpringDamperPulledBySpecifiedForceA.txt
   (2) % Copyright (c) 2017-2021 Motion Genesis LLC.  All rights reserved.
   (3) % Purpose: Simulate mass-spring-damper on rough table.
   (4) %-------------------------------------------------------------------------------
   (5) NewtonianFrame  N       % Earth (ground).
   (6) Particle        B       % Block that is connected by a spring/damper to N.
   (7) B.SetMass( m = 1 kg )
   (8) %-------------------------------------------------------------------------------
   (9) Variable   x''                    % B's horizontal displacement from No.
   (10) Variable   Fn                     % Resultant normal   force on B.
   (11) Variable   Ff                     % Resultant friction force on B.
   (12) Constant   g = 9.8 m/s^2          % Earth's gravitational acceleration.
   (13) Constant   k = 400 N/m            % Linear spring constant.
   (14) Constant   Ln = 4 m               % Natural length of spring.
   (15) Constant   b = 4 N*s/m            % Linear damper constant (zeta = 0.1).
   (16) Constant   muK = 0.4 noUnits      % Coefficient of kinetic friction.
   (17) Constant   epsilonV = 1.0E-5 m/s  % Represents a small sliding speed.
   (18) Specified  Fx = 20 * cos(t)       % Specified horizontal force on B.
-> (19) Fx = 20*cos(t)

   (20) %-------------------------------------------------------------------------------
   (21) %   Position, velocity, acceleration.
   (22) B.Translate(  No,  x*Nx>  )
-> (23) p_No_B> = x*Nx>
-> (24) v_B_N> = x'*Nx>
-> (25) a_B_N> = x''*Nx>

   (26) %-------------------------------------------------------------------------------
   (27) %   Forces on B.
   (28) B.AddForce( -m*g*Ny>  )                          % Gravity force.
-> (29) Force_B> = -m*g*Ny>

   (30) stretch = x - Ln                                 % Spring stretch.
-> (31) stretch = x - Ln

   (32) B.AddForce( -(k*stretch + b*Dt(stretch))*Nx> )   % Spring/damper force.
-> (33) Force_B> = (-k*stretch-b*x')*Nx> - m*g*Ny>

   (34) B.AddForce( Fn*Ny> + Ff*Nx> )                    % Normal and friction forces.
-> (35) Force_B> = (Ff-k*stretch-b*x')*Nx> + (Fn-m*g)*Ny>

   (36) B.AddForce( Fx*Nx> )                             % Specified force.
-> (37) Force_B> = (Fx+Ff-k*stretch-b*x')*Nx> + (Fn-m*g)*Ny>

   (38) %-------------------------------------------------------------------------------
   (39) %   Equations of motion via F = m*a.
   (40) Dynamics[1] = Dot( Nx>,  B.GetDynamics() )
-> (41) Dynamics[1] = k*stretch + b*x' + m*x'' - Fx - Ff

   (42) Dynamics[2] = Dot( Ny>,  B.GetDynamics() )
-> (43) Dynamics[2] = m*g - Fn

   (44) %-------------------------------------------------------------------------------
   (45) %   Equation governing Ff when B is sliding on N.
   (46) %   Herein we use the Continuous Friction Law to simulate both sticking and sliding.
   (47) %   This uses a small positive number for epsilonV to avoid divide-by-zero problems.
   (48) %   The advantage of the Continuous Friction Law is that it is easy to use and accurately
   (49) %   simulates both sticking and sliding.  The disadvantage is that the Continuous Friction
   (50) %   law is slow (stiff) when simulating sticking.  For better simulation speed, use the
   (51) %   Unified Friction Formulation in MGStickSlipMassSpringDamperPulledBySpecifiedForceB.txt.
   (52) magVelocity = B.GetSpeed( N )
-> (53) magVelocity = abs(x')

   (54) magVelocityPlusEpsilon = magVelocity + epsilonV
-> (55) magVelocityPlusEpsilon = epsilonV + magVelocity

   (56) Ff = Dot( -muK*Fn*B.GetVelocity(N) / magVelocityPlusEpsilon, Nx> )
-> (57) Ff = -muK*Fn*x'/magVelocityPlusEpsilon

   (58) %-------------------------------------------------------------------------------
   (59) %   Solve sliding equation of motion for x''.
   (60) Solve( Dynamics = 0,  x'', Fn )
-> (61) x'' = (Fx-k*stretch-b*x')/m - g*muK*x'/magVelocityPlusEpsilon
-> (62) Fn = m*g

   (63) %-------------------------------------------------------------------------------
   (64) %   Set initial values for variables (for subsequent ODE command).
   (65) Input  x = Input(Ln) m,  x' = 0 m/s
   (66) %-------------------------------------------------------------------------------
   (67) %   List output quantities from the ODE command.
   (68) Output  t sec,   x m,   x' m/s,   x'' m/s^2,  Ff Newton
   (69) %-------------------------------------------------------------------------------
   (70) %   Set numerical integration parameters and solve ODEs.
   (71) Input  tFinal = 12 sec,  tStep = 0.01 sec,  absError = 1.0E-07
   (72) ODE()  MGStickSlipMassSpringDamperPulledBySpecifiedForceA

   (73) %-------------------------------------------------------------------------------