In Lessons 1–3, we used ideal rotational elements and rigid kinematic constraints. Real drive trains violate rigidity: shafts twist, teeth deform, bearings dissipate energy, and gears have clearance. These effects introduce additional dynamic states (torsional modes) and piecewise behavior (backlash), which can dominate performance and stability margins in servo drives.

We will build models that satisfy: (i) \( \text{power consistency} \) across gear ratios, and (ii) \( \text{energy consistency} \) (stored energy is nonnegative; damping dissipates energy).

flowchart TD
  M["Motor inertia Jm"] --> S["Flexible shaft: ks, cs"]
  S --> G["Gear mesh: kg, cg, losses"]
  G --> B["Backlash clearance: width 2b"]
  B --> L["Load inertia Jl + load torque"]
  M --> U["Input torque tau_in"]
  L --> D["Disturbance/load torque tau_L"]
        

Notation used throughout: \( \theta_m \): motor angle, \( \theta_l \): load angle, \( \omega_m = \dot{\theta}_m \): motor speed, \( \omega_l = \dot{\theta}_l \): load speed, \( \tau_{in}(t) \): applied motor torque, \( \tau_L(t) \): external load torque.

Consider a uniform circular shaft of length \( L \), shear modulus \( G \), and polar second moment of area \( J_p \). Let the end angles be \( \theta_1 \) and \( \theta_2 \), so the twist is \( \phi = \theta_1 - \theta_2 \).

Under Saint-Venant torsion (small twist), the shear strain at radius \( r \) is \( \gamma(r) \):

\[ \gamma(r) = \frac{r}{L}\,(\theta_1 - \theta_2). \]

With linear elastic shear stress \( \sigma_{sh}(r) = G\gamma(r) \), the transmitted torque \( \tau_s \) is obtained by integrating the stress moment over the cross-section:

\[ \tau_s = \int_A \sigma_{sh}(r)\,r\,dA = \int_A \left(G\frac{r}{L}(\theta_1-\theta_2)\right) r\,dA = \frac{G}{L}(\theta_1-\theta_2)\int_A r^2\,dA. \]

By definition, \( J_p = \int_A r^2\,dA \), hence:

\[ \tau_s = \underbrace{\frac{GJ_p}{L}}_{k_s}\,(\theta_1-\theta_2), \qquad k_s = \frac{GJ_p}{L}. \]

For a solid circular shaft of diameter \( d \), \( J_p = \frac{\pi d^4}{32} \). For a hollow shaft with outer diameter \( d_o \) and inner diameter \( d_i \), \( J_p = \frac{\pi}{32}(d_o^4 - d_i^4) \).

To include material/internal damping (Kelvin–Voigt torsion), we add a torsional damper proportional to relative twist rate:

\[ \tau_s = k_s(\theta_1-\theta_2) + c_s(\dot{\theta}_1-\dot{\theta}_2), \qquad c_s > 0. \]

This model is lumped: it is accurate when the torsional wavelength is much larger than \( L \) in the operating frequency range (students saw this modeling philosophy in Chapter 1).

The canonical flexible drive train model is two inertias connected by a torsional spring–damper. Let \( J_m \) and \( J_l \) be motor and load inertias; the coupling torque is \( \tau_s \) as in Section 2.

\[ \begin{aligned} J_m \ddot{\theta}_m &= \tau_{in}(t) - \tau_s, \\ J_l \ddot{\theta}_l &= \tau_s - \tau_L(t), \\ \tau_s &= k_s(\theta_m-\theta_l) + c_s(\dot{\theta}_m-\dot{\theta}_l). \end{aligned} \]

Decoupling into rigid-body and torsional modes (derivation). Define the relative twist \( q \): and the inertia-weighted average angle \( \psi \):

\[ q = \theta_m - \theta_l, \qquad \psi = \frac{J_m\theta_m + J_l\theta_l}{J_m+J_l}. \]

Differentiate twice. Using the two equations of motion and noting that internal torque cancels in the sum:

\[ (J_m+J_l)\ddot{\psi} = J_m\ddot{\theta}_m + J_l\ddot{\theta}_l = \tau_{in}(t) - \tau_L(t). \]

Thus \( \psi \) is a rigid-body coordinate driven only by external torques. For the relative coordinate, subtract the load equation from the motor equation after scaling to eliminate \( \tau_s \). A clean way is to write:

\[ \ddot{q}=\ddot{\theta}_m-\ddot{\theta}_l = \frac{\tau_{in}-\tau_s}{J_m} - \frac{\tau_s-\tau_L}{J_l}. \]

Rearranging gives:

\[ \underbrace{\frac{J_mJ_l}{J_m+J_l}}_{J_{red}}\ddot{q} + c_s\dot{q} + k_s q = \frac{J_l}{J_m+J_l}\tau_{in}(t) + \frac{J_m}{J_m+J_l}\tau_L(t), \]

where \( J_{red} = \frac{J_mJ_l}{J_m+J_l} \) is the reduced inertia. In free vibration (\( \tau_{in}=\tau_L=0 \)), the torsional mode obeys:

\[ J_{red}\ddot{q} + c_s\dot{q} + k_s q = 0. \]

Therefore the undamped torsional natural frequency and damping ratio are:

\[ \omega_n = \sqrt{\frac{k_s}{J_{red}}}, \qquad \zeta = \frac{c_s}{2\sqrt{k_sJ_{red}}}. \]

Energy dissipation proof (unforced stability). For \( \tau_{in}=\tau_L=0 \), define total stored energy:

\[ V = \frac{1}{2}J_m\omega_m^2 + \frac{1}{2}J_l\omega_l^2 + \frac{1}{2}k_s(\theta_m-\theta_l)^2. \]

Differentiate and substitute dynamics:

\[ \begin{aligned} \dot{V} &= J_m\omega_m\dot{\omega}_m + J_l\omega_l\dot{\omega}_l + k_s(\theta_m-\theta_l)(\omega_m-\omega_l) \\ &= \omega_m(-\tau_s) + \omega_l(\tau_s) + k_s(\theta_m-\theta_l)(\omega_m-\omega_l) \\ &= -(\omega_m-\omega_l)\tau_s + k_s(\theta_m-\theta_l)(\omega_m-\omega_l). \end{aligned} \]

Using \( \tau_s = k_s(\theta_m-\theta_l)+c_s(\omega_m-\omega_l) \) yields:

\[ \dot{V} = -c_s(\omega_m-\omega_l)^2 \le 0. \]

Hence the coupling damper dissipates energy monotonically; this is the key physical constraint your models should preserve.

From Lesson 2, an ideal gear pair acts like a power-preserving transformer. Define gear ratio \( g \): motor speed divided by load speed: \( g = \frac{\omega_m}{\omega_l} \) (typically \( g > 1 \) for reduction). Neglecting sign conventions, ideal kinematics imply:

\[ \omega_l = \frac{1}{g}\omega_m, \qquad \theta_l = \frac{1}{g}\theta_m \; (+\;\text{constant}). \]

Power preservation (no loss) implies:

\[ \tau_m\omega_m = \tau_l\omega_l \quad \Rightarrow \quad \tau_l = g\, \tau_m. \]

Reflection of inertia (proof). Consider a load inertia \( J_l \) on the slow side. Load kinetic energy is: \( T_l = \frac{1}{2}J_l\omega_l^2 \). Express it in terms of motor-side speed \( \omega_m \) using \( \omega_l=\omega_m/g \):

\[ T_l = \frac{1}{2}J_l\left(\frac{\omega_m}{g}\right)^2 = \frac{1}{2}\underbrace{\left(\frac{J_l}{g^2}\right)}_{J_{l,ref}}\omega_m^2. \]

Therefore the reflected inertia seen at the motor shaft is: \( J_{l,ref} = \frac{J_l}{g^2} \). This result is a direct consequence of energy consistency.

Gear mesh compliance and damping. Tooth contact deforms elastically and dissipates energy. A common lumped model uses a spring–damper in the mesh deflection:

\[ \delta = \theta_m - g\theta_l, \qquad \tau_{mesh} = k_g\delta + c_g\dot{\delta}, \qquad k_g>0,\;c_g>0. \]

Here \( \delta \) is the relative angular “closure” of the mesh after accounting for ratio. If you also have a flexible shaft between motor and gear, then you will have two elastic elements (shaft and mesh). In many drive trains, mesh stiffness is significantly higher than shaft stiffness, but it is rarely infinite.

Efficiency / losses (simple dynamical inclusion). If the gear train is not lossless, a minimal consistent model is to include viscous loss torques (or an equivalent mesh damping increase). A simple static efficiency \( \eta \in (0,1] \) can be approximated by scaling transmitted torque:

\[ \tau_l \approx \eta\,g\, \tau_m, \qquad 0 < \eta \le 1. \]

In dynamics, it is generally preferable (and more physically transparent) to represent losses via damping terms rather than a discontinuous efficiency map.

Backlash is a finite clearance between mating teeth. When the relative mesh deflection \( \delta=\theta_m-g\theta_l \) lies within a dead band, teeth do not contact and no elastic torque is transmitted. Let half-width clearance be \( b > 0 \), so the free-play region is \( \lvert \delta \rvert \le b \).

A widely used backlash with mesh stiffness model is:

\[ \tau_{mesh}(\delta,\dot{\delta}) = \begin{cases} k_g(\delta-b) + c_g\dot{\delta}, & \delta > b, \\ 0, & \lvert \delta \rvert \le b, \\ k_g(\delta+b) + c_g\dot{\delta}, & \delta < -b. \end{cases} \]

Key properties you should internalize:

• Non-smooth stiffness: the slope switches between \( 0 \) and \( k_g \).
• Hybrid contact: the dynamics alternate between “free” (no torque) and “engaged” (spring–damper).
• Energy consistency: in engagement, the elastic energy \( \frac{1}{2}k_g(\delta \mp b)^2 \) is nonnegative; damping dissipates energy through \( c_g\dot{\delta}^2 \).

Local linearization (engaged). If the motion remains on one engaged side (e.g., \( \delta > b \)), then defining \( \tilde{\delta} = \delta - b \) yields a purely linear relation \( \tau_{mesh} = k_g\tilde{\delta} + c_g\dot{\tilde{\delta}} \). The offset \( b \) shifts the equilibrium but does not change stiffness in engagement.

A compact and practical model for many servo systems places elasticity and backlash at the gear mesh. Let \( J_m \) be motor-side inertia and \( J_l \) load-side inertia. The torque transmitted through mesh is \( \tau_{mesh} \); the load sees amplified torque by \( g \).

\[ \begin{aligned} J_m \ddot{\theta}_m &= \tau_{in}(t) - \tau_{mesh}(\delta,\dot{\delta}) - b_m\dot{\theta}_m, \\ J_l \ddot{\theta}_l &= g\, \tau_{mesh}(\delta,\dot{\delta}) - \tau_L(t) - b_l\dot{\theta}_l, \\ \delta &= \theta_m - g\theta_l, \qquad \dot{\delta} = \dot{\theta}_m - g\dot{\theta}_l, \end{aligned} \]

where \( b_m \ge 0 \) and \( b_l \ge 0 \) represent viscous losses at motor and load bearings. If you also have a flexible shaft between motor and gear, introduce an intermediate angle and another spring–damper element (this produces a 3-inertia torsional chain).

Energy dissipation in engaged mode (proof sketch). Assume the mesh is engaged on one side so that: \( \tau_{mesh} = k_g\tilde{\delta} + c_g\dot{\tilde{\delta}} \) with \( \tilde{\delta}=\delta-b \) or \( \tilde{\delta}=\delta+b \). Consider the storage function:

\[ V = \frac{1}{2}J_m\omega_m^2 + \frac{1}{2}J_l\omega_l^2 + \frac{1}{2}k_g\tilde{\delta}^2. \]

For the unforced case (\( \tau_{in}=\tau_L=0 \)), differentiating and substituting dynamics yields:

\[ \dot{V} = -b_m\omega_m^2 - b_l\omega_l^2 - c_g\dot{\tilde{\delta}}^2 \le 0. \]

This is the same physical conclusion as in Section 3: damping mechanisms guarantee non-increasing energy. When backlash enters the free-play region (\( \tau_{mesh}=0 \)), the elastic term disappears and the remaining viscous losses still dissipate kinetic energy.

Backlash introduces piecewise dynamics. Robust simulation requires a clear switching rule: compute \( \delta \), decide whether contact exists, then compute \( \tau_{mesh} \). For many educational and engineering analyses, the static dead-band law in Section 5 is sufficient.

flowchart TD
  X["Compute delta = theta_m - g*theta_l"] --> A["Check |delta| <= b ?"]
  A -->|yes| Z["tau_mesh = 0"]
  A -->|no| B1["If delta > b"]
  A -->|no| B2["If delta < -b"]
  B1 --> T1["tau_mesh = kg*(delta-b) + cg*delta_dot"]
  B2 --> T2["tau_mesh = kg*(delta+b) + cg*delta_dot"]
  Z --> O["Integrate ODEs one step"]
  T1 --> O
  T2 --> O
        

Practical note: switching creates non-smooth acceleration. Use sufficiently small steps (fixed-step solvers) or a variable-step integrator with tight tolerances (adaptive solvers).

We simulate the unified model of Section 6 using scipy.integrate.solve_ivp. Recommended libraries for System Dynamics in Python at this stage: numpy (arrays), scipy (ODE integration).

import numpy as np
from scipy.integrate import solve_ivp

def backlash_tau(delta, delta_dot, kg, cg, b):
    # Piecewise backlash torque
    if delta > b:
        return kg*(delta - b) + cg*delta_dot
    elif delta < -b:
        return kg*(delta + b) + cg*delta_dot
    else:
        return 0.0

def tau_in(t):
    # Example: step torque + small sinusoid
    return 1.0 + 0.2*np.sin(10.0*t)

def tau_L(t):
    # Example load disturbance
    return 0.3

def odefun(t, x, params):
    theta_m, omega_m, theta_l, omega_l = x
    Jm, Jl, g, kg, cg, b, bm, bl = params

    delta = theta_m - g*theta_l
    delta_dot = omega_m - g*omega_l

    tau_mesh = backlash_tau(delta, delta_dot, kg, cg, b)

    domega_m = (tau_in(t) - tau_mesh - bm*omega_m) / Jm
    domega_l = (g*tau_mesh - tau_L(t) - bl*omega_l) / Jl

    return [omega_m, domega_m, omega_l, domega_l]

# Parameters
Jm = 2.0e-4
Jl = 8.0e-3
g  = 15.0
kg = 50.0
cg = 0.02
b  = 2.0e-3
bm = 1.0e-4
bl = 2.0e-3

params = (Jm, Jl, g, kg, cg, b, bm, bl)

# Initial conditions: small mismatch inside backlash
x0 = [0.0, 0.0, 0.0, 0.0]
t_span = (0.0, 0.8)

sol = solve_ivp(lambda t, x: odefun(t, x, params), t_span, x0,
                method="RK45", rtol=1e-8, atol=1e-10, max_step=1e-4)

t = sol.t
theta_m, omega_m, theta_l, omega_l = sol.y

# Post-processing: compute mesh torque history
tau_mesh_hist = np.zeros_like(t)
for i in range(len(t)):
    delta = theta_m[i] - g*theta_l[i]
    delta_dot = omega_m[i] - g*omega_l[i]
    tau_mesh_hist[i] = backlash_tau(delta, delta_dot, kg, cg, b)

print("Sim done. Final speeds (omega_m, omega_l) =", omega_m[-1], omega_l[-1])
print("Fraction of time in contact (|delta| > b):",
      np.mean(np.abs(theta_m - g*theta_l) > b))
      

Suggested exercises: (i) vary \( b \) to see how clearance increases “free motion” time; (ii) vary \( k_g \) to push torsional resonance higher; (iii) set \( c_g=0 \) and observe that oscillations persist longer (energy is not dissipated at the mesh).

In C++, common System Dynamics-friendly choices include Boost.Odeint for ODE integration and Eigen for linear algebra. Below is a minimal from-scratch fixed-step RK4 integrator, which is sufficient for backlash switching when you use a small time step.

#include <iostream>
#include <cmath>

struct Params {
  double Jm, Jl, g, kg, cg, b, bm, bl;
};

static double tau_in(double t) {
  return 1.0 + 0.2*std::sin(10.0*t);
}

static double tau_L(double /*t*/) {
  return 0.3;
}

static double backlash_tau(double delta, double delta_dot, const Params& p) {
  if (delta > p.b) {
    return p.kg*(delta - p.b) + p.cg*delta_dot;
  } else if (delta < -p.b) {
    return p.kg*(delta + p.b) + p.cg*delta_dot;
  } else {
    return 0.0;
  }
}

// State x = [theta_m, omega_m, theta_l, omega_l]
static void f(double t, const double x[4], double dxdt[4], const Params& p) {
  const double theta_m = x[0], omega_m = x[1];
  const double theta_l = x[2], omega_l = x[3];

  const double delta     = theta_m - p.g*theta_l;
  const double delta_dot = omega_m - p.g*omega_l;

  const double tau_mesh = backlash_tau(delta, delta_dot, p);

  dxdt[0] = omega_m;
  dxdt[1] = (tau_in(t) - tau_mesh - p.bm*omega_m) / p.Jm;
  dxdt[2] = omega_l;
  dxdt[3] = (p.g*tau_mesh - tau_L(t) - p.bl*omega_l) / p.Jl;
}

static void rk4_step(double t, double h, double x[4], const Params& p) {
  double k1[4], k2[4], k3[4], k4[4], xtmp[4];

  f(t, x, k1, p);
  for (int i=0; i<4; ++i) xtmp[i] = x[i] + 0.5*h*k1[i];

  f(t + 0.5*h, xtmp, k2, p);
  for (int i=0; i<4; ++i) xtmp[i] = x[i] + 0.5*h*k2[i];

  f(t + 0.5*h, xtmp, k3, p);
  for (int i=0; i<4; ++i) xtmp[i] = x[i] + h*k3[i];

  f(t + h, xtmp, k4, p);

  for (int i=0; i<4; ++i) {
    x[i] = x[i] + (h/6.0)*(k1[i] + 2.0*k2[i] + 2.0*k3[i] + k4[i]);
  }
}

int main() {
  Params p;
  p.Jm = 2.0e-4; p.Jl = 8.0e-3; p.g = 15.0;
  p.kg = 50.0; p.cg = 0.02; p.b = 2.0e-3;
  p.bm = 1.0e-4; p.bl = 2.0e-3;

  double x[4] = {0.0, 0.0, 0.0, 0.0};

  const double t0 = 0.0, tf = 0.8;
  const double h = 1.0e-5; // small step to handle switching
  const int N = static_cast<int>((tf - t0)/h);

  double t = t0;
  for (int k=0; k<N; ++k) {
    rk4_step(t, h, x, p);
    t += h;
  }

  std::cout << "Final omega_m = " << x[1] << "\n";
  std::cout << "Final omega_l = " << x[3] << "\n";
  return 0;
}
      

Engineering note: backlash creates high-frequency content at switching; fixed-step RK4 is acceptable if \( h \) is small. For higher efficiency, use Boost.Odeint with adaptive stepping.

In Java, Apache Commons Math provides ODE solvers, and libraries like EJML support matrix operations. Below is a lightweight RK4 implementation consistent with the model of Section 6.

public class BacklashDriveTrainRK4 {

  static class Params {
    double Jm, Jl, g, kg, cg, b, bm, bl;
  }

  static double tauIn(double t) {
    return 1.0 + 0.2*Math.sin(10.0*t);
  }

  static double tauL(double t) {
    return 0.3;
  }

  static double backlashTau(double delta, double deltaDot, Params p) {
    if (delta > p.b) {
      return p.kg*(delta - p.b) + p.cg*deltaDot;
    } else if (delta < -p.b) {
      return p.kg*(delta + p.b) + p.cg*deltaDot;
    } else {
      return 0.0;
    }
  }

  // x = [theta_m, omega_m, theta_l, omega_l]
  static void f(double t, double[] x, double[] dxdt, Params p) {
    double thetaM = x[0], omegaM = x[1];
    double thetaL = x[2], omegaL = x[3];

    double delta = thetaM - p.g*thetaL;
    double deltaDot = omegaM - p.g*omegaL;

    double tauMesh = backlashTau(delta, deltaDot, p);

    dxdt[0] = omegaM;
    dxdt[1] = (tauIn(t) - tauMesh - p.bm*omegaM) / p.Jm;
    dxdt[2] = omegaL;
    dxdt[3] = (p.g*tauMesh - tauL(t) - p.bl*omegaL) / p.Jl;
  }

  static void rk4Step(double t, double h, double[] x, Params p) {
    double[] k1 = new double[4];
    double[] k2 = new double[4];
    double[] k3 = new double[4];
    double[] k4 = new double[4];
    double[] xt = new double[4];

    f(t, x, k1, p);
    for (int i=0; i<4; i++) xt[i] = x[i] + 0.5*h*k1[i];

    f(t + 0.5*h, xt, k2, p);
    for (int i=0; i<4; i++) xt[i] = x[i] + 0.5*h*k2[i];

    f(t + 0.5*h, xt, k3, p);
    for (int i=0; i<4; i++) xt[i] = x[i] + h*k3[i];

    f(t + h, xt, k4, p);

    for (int i=0; i<4; i++) {
      x[i] = x[i] + (h/6.0)*(k1[i] + 2.0*k2[i] + 2.0*k3[i] + k4[i]);
    }
  }

  public static void main(String[] args) {
    Params p = new Params();
    p.Jm = 2.0e-4; p.Jl = 8.0e-3; p.g = 15.0;
    p.kg = 50.0; p.cg = 0.02; p.b = 2.0e-3;
    p.bm = 1.0e-4; p.bl = 2.0e-3;

    double[] x = new double[] {0.0, 0.0, 0.0, 0.0};

    double t0 = 0.0, tf = 0.8;
    double h = 1e-5;
    int N = (int)((tf - t0)/h);

    double t = t0;
    for (int k=0; k<N; k++) {
      rk4Step(t, h, x, p);
      t += h;
    }

    System.out.println("Final omega_m = " + x[1]);
    System.out.println("Final omega_l = " + x[3]);
  }
}
      

If you migrate to Apache Commons Math, the same right-hand side f(t,x) can be used with adaptive integrators (Dormand–Prince), which is valuable for piecewise systems.

MATLAB offers ode45 for variable-step simulation. In Simulink, you can implement backlash using the standard Backlash block (or a custom dead-zone + memory logic), and mesh compliance via a spring–damper.

function backlash_drivetrain_demo()
  % Parameters
  Jm = 2.0e-4;  Jl = 8.0e-3;  g = 15.0;
  kg = 50.0;    cg = 0.02;    b = 2.0e-3;
  bm = 1.0e-4;  bl = 2.0e-3;

  x0 = [0; 0; 0; 0]; % [theta_m; omega_m; theta_l; omega_l]
  tspan = [0 0.8];

  opts = odeset('RelTol',1e-8,'AbsTol',1e-10,'MaxStep',1e-4);
  sol = ode45(@(t,x) rhs(t,x,Jm,Jl,g,kg,cg,b,bm,bl), tspan, x0, opts);

  t = linspace(tspan(1), tspan(2), 2000);
  X = deval(sol, t);

  theta_m = X(1,:); omega_m = X(2,:);
  theta_l = X(3,:); omega_l = X(4,:);

  delta = theta_m - g*theta_l;
  delta_dot = omega_m - g*omega_l;

  tau_mesh = arrayfun(@(d,dd) backlash_tau(d,dd,kg,cg,b), delta, delta_dot);

  disp(['Final omega_m = ', num2str(omega_m(end))]);
  disp(['Final omega_l = ', num2str(omega_l(end))]);
  disp(['Contact fraction = ', num2str(mean(abs(delta) > b))]);
end

function dx = rhs(t,x,Jm,Jl,g,kg,cg,b,bm,bl)
  theta_m = x(1); omega_m = x(2);
  theta_l = x(3); omega_l = x(4);

  delta = theta_m - g*theta_l;
  delta_dot = omega_m - g*omega_l;

  tau_mesh = backlash_tau(delta, delta_dot, kg, cg, b);

  tau_in = 1.0 + 0.2*sin(10.0*t);
  tau_L  = 0.3;

  dx = zeros(4,1);
  dx(1) = omega_m;
  dx(2) = (tau_in - tau_mesh - bm*omega_m)/Jm;
  dx(3) = omega_l;
  dx(4) = (g*tau_mesh - tau_L - bl*omega_l)/Jl;
end

function tm = backlash_tau(delta, delta_dot, kg, cg, b)
  if delta > b
    tm = kg*(delta - b) + cg*delta_dot;
  elseif delta < -b
    tm = kg*(delta + b) + cg*delta_dot;
  else
    tm = 0.0;
  end
end
      

Simulink guidance (block-level): Use integrators for \( \omega_m,\theta_m,\omega_l,\theta_l \). Compute \( \delta=\theta_m-g\theta_l \) and \( \dot{\delta}=\omega_m-g\omega_l \), then implement the piecewise map for \( \tau_{mesh} \) using Switch blocks (or a Backlash block feeding a spring–damper). Feed \( -\tau_{mesh} \) into the motor acceleration and \( g\tau_{mesh} \) into the load acceleration.

Mathematica’s NDSolve handles piecewise right-hand sides effectively when tolerances are set carefully.

ClearAll["Global`*"];

(* Parameters *)
Jm = 2.0*10^-4; Jl = 8.0*10^-3; g = 15.0;
kg = 50.0; cg = 0.02; b = 2.0*10^-3;
bm = 1.0*10^-4; bl = 2.0*10^-3;

tauIn[t_] := 1.0 + 0.2*Sin[10.0*t];
tauL[t_] := 0.3;

tauMesh[delta_, deltaDot_] := Piecewise[{
  {kg*(delta - b) + cg*deltaDot, delta > b},
  {kg*(delta + b) + cg*deltaDot, delta < -b}
}, 0.0];

eqns = {
  θm'[t] == ωm[t],
  θl'[t] == ωl[t],

  Jm*ωm'[t] == tauIn[t] - tauMesh[θm[t] - g*θl[t], ωm[t] - g*ωl[t]] - bm*ωm[t],
  Jl*ωl'[t] == g*tauMesh[θm[t] - g*θl[t], ωm[t] - g*ωl[t]] - tauL[t] - bl*ωl[t],

  θm[0] == 0, ωm[0] == 0,
  θl[0] == 0, ωl[0] == 0
};

sol = NDSolve[eqns, {θm, ωm, θl, ωl}, {t, 0, 0.8},
  Method -> {"TimeIntegration" -> "Adams"},
  AccuracyGoal -> 12, PrecisionGoal -> 12, MaxStepSize -> 10^-4
];

ωmFinal = ωm[0.8] /. sol[[1]];
ωlFinal = ωl[0.8] /. sol[[1]];

Print["Final omega_m = ", ωmFinal];
Print["Final omega_l = ", ωlFinal];
      

Extension: compute the contact indicator \( \mathbf{1}_{\{|\delta|>b\}} \) over time to estimate how often the gear train is transmitting torque versus “free-wheeling” in clearance.

Problem 1 (Shaft stiffness from torsion): A hollow shaft has length \( L \), shear modulus \( G \), outer diameter \( d_o \), inner diameter \( d_i \). Derive \( k_s \) such that \( \tau_s = k_s(\theta_1-\theta_2) \).

Solution: From Section 2, \( k_s = \frac{GJ_p}{L} \). For a hollow circular cross-section:

\[ J_p = \int_A r^2\,dA = \frac{\pi}{32}(d_o^4-d_i^4) \quad \Rightarrow \quad k_s = \frac{G}{L}\frac{\pi}{32}(d_o^4-d_i^4). \]

Problem 2 (Torsional natural frequency of two inertias): Two inertias \( J_m \) and \( J_l \) are connected by torsional spring \( k_s \) and damper \( c_s \). With no external torques, derive the torsional mode equation and give \( \omega_n \), \( \zeta \).

Solution: Using Section 3, define \( q=\theta_m-\theta_l \). Then:

\[ J_{red}\ddot{q} + c_s\dot{q} + k_s q = 0, \qquad J_{red}=\frac{J_mJ_l}{J_m+J_l}. \]

Hence:

\[ \omega_n=\sqrt{\frac{k_s}{J_{red}}},\qquad \zeta=\frac{c_s}{2\sqrt{k_sJ_{red}}}. \]

Problem 3 (Inertia reflection across gear ratio): An ideal gear train has ratio \( g=\omega_m/\omega_l \). Show that a load inertia \( J_l \) appears as \( J_l/g^2 \) when viewed from the motor side.

Solution: Energy consistency gives:

\[ T_l=\frac{1}{2}J_l\omega_l^2=\frac{1}{2}J_l\left(\frac{\omega_m}{g}\right)^2 =\frac{1}{2}\left(\frac{J_l}{g^2}\right)\omega_m^2. \]

Therefore the motor-side equivalent inertia that yields the same kinetic energy for the same \( \omega_m \) is \( J_{l,ref}=\frac{J_l}{g^2} \).

Problem 4 (Backlash engagement condition and torque): Let \( \delta(t)=A\sin(\Omega t) \) and let backlash half-width be \( b \). (i) For what condition on \( A \) does engagement occur? (ii) Assuming \( c_g=0 \), write \( \tau_{mesh}(t) \) explicitly using the piecewise law.

Solution: (i) Engagement requires \( \lvert \delta(t)\rvert > b \) for some \( t \), so \( A > b \). (ii) With \( c_g=0 \):

\[ \tau_{mesh}(t)= \begin{cases} k_g(A\sin(\Omega t)-b), & A\sin(\Omega t) > b, \\ 0, & \lvert A\sin(\Omega t)\rvert \le b, \\ k_g(A\sin(\Omega t)+b), & A\sin(\Omega t) < -b. \end{cases} \]

Problem 5 (Minimum shaft diameter to push torsional mode above a target): A solid shaft (diameter \( d \), length \( L \), shear modulus \( G \)) connects two inertias \( J_m \), \( J_l \). Choose \( d \) so that the undamped torsional natural frequency satisfies \( \omega_n \ge \omega_{min} \).

Solution: For a solid shaft, \( k_s=\frac{GJ_p}{L}=\frac{G}{L}\frac{\pi d^4}{32} \). Using \( \omega_n=\sqrt{k_s/J_{red}} \) with \( J_{red}=\frac{J_mJ_l}{J_m+J_l} \), the constraint \( \omega_n^2 \ge \omega_{min}^2 \) becomes:

\[ \frac{1}{J_{red}}\frac{G}{L}\frac{\pi d^4}{32} \ge \omega_{min}^2 \quad \Rightarrow \quad d \ge \left(\frac{32LJ_{red}\omega_{min}^2}{G\pi}\right)^{1/4}. \]

This design rule is frequently used as a first pass before refining the model with additional elastic elements (gear mesh, couplings) and damping identification.

We derived flexible-shaft torque from torsion mechanics, built the two-inertia torsional resonance model, and proved energy dissipation under positive damping. We then extended ideal gear transformers to include mesh compliance, losses, and backlash via a piecewise torque law, yielding a physically consistent drive-train model. Finally, we implemented simulations across Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica, emphasizing robust handling of backlash switching.