Sprung and unsprung mass

%% Sprung and unsprung mass: the quarter car
%
% Companion to the documentation page "Sprung and unsprung mass". Runs unchanged in GNU
% Octave (7 or later) and MATLAB. Edit the car in example_car.m and run:
%
%   >> quarter_car
%
% One corner of the car as two masses: the sprung mass m_s (its share of the body) on the
% suspension spring k_s and damper c_s, and the unsprung mass m_u (wheel, upright, brake)
% on the tyre k_t. The road z_r moves the tyre; z_u and z_s are the wheel and the body.
%
%   m_s z_s'' = -k_s (z_s - z_u) - c_s (z_s' - z_u')
%   m_u z_u'' =  k_s (z_s - z_u) + c_s (z_s' - z_u') - k_t (z_u - z_r) - c_t (z_u' - z_r')
%
% Two modes: the body bouncing on the spring and tyre in series (ride), and the wheel
% bouncing between the two (wheel hop). The ride frequency is the textbook one; the
% damping ratio is read off the ride mode itself, because the damper compresses by the
% spring's share of the body's motion, not all of it.
%
% Reference: T. D. Gillespie, Fundamentals of Vehicle Dynamics (SAE, 1992), ch. 5.
%
% Part of the PaceDynamics documentation, https://pacedyn.dev/docs
% SPDX-License-Identifier: MIT

[car, app] = example_car();
app = app.mass;

g = 9.81;                                           % m/s^2
series = @(a, b) a .* b ./ (a + b);

%% Masses
m_u  = 4 * car.unsprungPerCorner;                   % kg, the four corners
m_s  = car.massTotal - m_u;                         % kg, sprung
m_sf = m_s * car.frontMassFrac;                     % kg, sprung, on the front axle
m_sr = m_s * (1 - car.frontMassFrac);
W    = car.massTotal * g;                           % N
static = W * [car.frontMassFrac, 1 - car.frontMassFrac] / 2;   % N, one wheel, front and rear

pd_compare('Sprung mass [kg]',              m_s,       app.sprungMass);
pd_compare('Sprung mass, front axle [kg]',  m_sf,      app.sprungFront);
pd_compare('Sprung mass, rear axle [kg]',   m_sr,      app.sprungRear);
pd_compare('One front wheel, standing [N]', static(1), app.staticFront);
pd_compare('One rear wheel, standing [N]',  static(2), app.staticRear);

%% A quarter car per axle
% The suspension at the wheel in heave: corner spring and heave element, corner damper
% and heave damper (a roll element does not move when both wheels rise together).
axle = {'front', 'rear'};
ms = [m_sf, m_sr] / 2;                                                   % kg, one corner
ks = [car.wheelRateFront + car.heaveSpringFront, car.wheelRateRear + car.heaveSpringRear] * 1e3;   % N/m
cs = [car.damperFront + car.heaveDamperFront, car.damperRear + car.heaveDamperRear] * 1e3;         % N*s/m
kt = [car.tireRateFront, car.tireRateRear] * 1e3;                        % N/m
ct = car.tireDamping * 1e3;                                              % N*s/m
mu = max(0.1, car.unsprungCornerMass);                                   % kg

for i = 1:2
  a = app.(axle{i});
  % Ride: the body on the spring and the tyre in series. Wheel hop: the wheel on both at once.
  f_ride = sqrt(series(ks(i), kt(i)) / ms(i)) / (2 * pi);
  f_hop  = sqrt((ks(i) + kt(i)) / mu) / (2 * pi);

  % The two modes of the state matrix, state [z_s, z_s', z_u, z_u'].
  A = [0, 1, 0, 0;
       -ks(i) / ms(i), -cs(i) / ms(i),  ks(i) / ms(i),            cs(i) / ms(i);
       0, 0, 0, 1;
        ks(i) / mu,     cs(i) / mu,    -(ks(i) + kt(i)) / mu,   -(cs(i) + ct) / mu];
  [V, D] = eig(A);
  lambda = diag(D);
  % The ride mode: the oscillating one whose motion is mostly the body's (by kinetic energy).
  best = 0;  zeta = NaN;
  for j = 1:4
    if imag(lambda(j)) <= 0, continue; end
    es = ms(i) * abs(V(1, j))^2;  eu = mu * abs(V(3, j))^2;
    share = es / (es + eu);
    if share >= 0.5 && share > best
      best = share;  zeta = -real(lambda(j)) / abs(lambda(j));
    end
  end

  % Transmissibility: how far the body moves for each millimetre of road, at frequency f.
  w = 2 * pi * a.gainFreq;
  spring = ks(i) + 1i * cs(i) * w;                  % the suspension's complex stiffness
  tyre   = kt(i) + 1i * ct * w;
  H = spring .* tyre ./ ((spring - ms(i) * w.^2) .* (spring + tyre - mu * w.^2) - spring.^2);

  pd_compare(sprintf('Ride frequency, %s [Hz]', axle{i}),      f_ride, a.rideFreq);
  pd_compare(sprintf('Wheel-hop frequency, %s [Hz]', axle{i}), f_hop,  a.wheelHopFreq);
  pd_compare(sprintf('Ride damping ratio, %s [-]', axle{i}),   zeta,   a.dampingRatio);
  for j = 1:numel(w)
    pd_compare(sprintf('  |z_s/z_r| at %4.1f Hz, %s', a.gainFreq(j), axle{i}), abs(H(j)), a.gain(j));
  end
end
pd_compare();