By the end of this lesson, you should be able to:

• Apply torque balance to derive rotational equations of motion for lumped inertias.
• Compute the mass moment of inertia from discrete sums and continuous integrals.
• Prove the parallel-axis (Huygens–Steiner) theorem and use it to shift inertia about axes.
• Model torsional springs and dampers and derive their stored/dissipated energy relations.
• Simulate a torsional oscillator in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

Assumption set (linear lumped rotational model):

• Rigid bodies with a single dominant rotation about a fixed axis.
• Small angular deflections so linear torsional constitutive laws are valid.
• Viscous torsional damping (torque proportional to angular velocity).
• Parameters \( J \), \( k_\theta \), \( b_\theta \) are constant (time-invariant).

Let \( \theta(t) \): angular displacement (rad), \( \omega(t) \): angular velocity (rad/s), \( \alpha(t) \): angular acceleration (rad/s²).

\[ \omega(t) = \frac{d\theta(t)}{dt}, \qquad \alpha(t) = \frac{d\omega(t)}{dt} = \frac{d^2\theta(t)}{dt^2}. \]

For a point at distance \( r \) from the axis, the tangential speed is:

\[ v_t(t) = r\,\omega(t), \qquad a_t(t) = r\,\alpha(t). \]

These relations allow us to connect rotational motion to the translational kinetic energy of particles at radius \( r \).

The rotational analog of Newton’s second law (used in Chapter 3) is that the net applied torque about the axis equals the time rate of change of angular momentum. For planar rotation about a fixed axis, this reduces to: \( \sum T(t) = J\,\alpha(t) \). (Torque is often denoted \( \tau \) in textbooks; here we use \( T \).)

Derivation from particle kinetics (rigorous, lumped axis):

Consider a rigid body composed of particles \( i=1,\dots,N \) with masses \( m_i \) at radii \( r_i \) from the rotation axis. Under rigid rotation with angular speed \( \omega \), particle \( i \) has speed \( v_i = r_i\omega \), so its kinetic energy is \( \tfrac{1}{2}m_i v_i^2 \). Summing over particles gives:

\[ K = \sum_{i=1}^N \frac{1}{2} m_i (r_i\omega)^2 = \frac{1}{2}\left(\sum_{i=1}^N m_i r_i^2\right)\omega^2. \]

Define the mass moment of inertia about the axis: \( J \):

\[ J \triangleq \sum_{i=1}^N m_i r_i^2. \]

The rotational work-energy identity is \( dW = T\,d\theta \), hence power is:

\[ P = \frac{dW}{dt} = T\,\omega. \]

But also \( P = \tfrac{dK}{dt} \). Using \( K=\tfrac{1}{2}J\omega^2 \):

\[ \frac{dK}{dt} = \frac{d}{dt}\left(\frac{1}{2}J\omega^2\right)=J\omega\frac{d\omega}{dt} = J\omega\alpha. \]

Equate power expressions (assuming \( \omega \neq 0 \) and by continuity also for \( \omega=0 \)):

\[ T\,\omega = J\omega\alpha \;\;\Rightarrow\;\; T = J\alpha. \]

flowchart TD
  Tin["Applied torque: \nT_in(t)"] --> SUM["Sum torques about axis"]
  K["Torsional spring torque: \nT_k = k*theta"] --> SUM
  B["Torsional damper torque: \nT_b = b*omega"] --> SUM
  SUM --> LAW["Dynamics: J*theta_ddot = T_in - T_k - T_b"]
  LAW --> ODE["ODE: J*theta_ddot + b*theta_dot + k*theta = T_in(t)"]
        

For a continuous rigid body with mass density \( \rho(\mathbf{r}) \), the discrete sum becomes an integral. If \( r \) denotes perpendicular distance to the rotation axis:

\[ J = \int r^2 \, dm = \int_V r^2 \rho(\mathbf{r})\, dV. \]

Example 1: Thin ring (radius \( R \), mass \( m \))

All mass lies at \( r=R \), so:

\[ J = \int r^2\,dm = \int R^2\,dm = R^2\int dm = mR^2. \]

Example 2: Solid disk (radius \( R \), mass \( m \))

Use polar area element \( dA = r\,dr\,d\phi \) and surface density \( \sigma = \frac{m}{\pi R^2} \). Then \( dm = \sigma\,dA \) and:

\[ J = \int_0^{2\pi}\int_0^R r^2(\sigma r\,dr\,d\phi) = \sigma\int_0^{2\pi}\int_0^R r^3\,dr\,d\phi = \sigma(2\pi)\left(\frac{R^4}{4}\right) = \frac{1}{2}mR^2. \]

Example 3: Slender rod (length \( L \), mass \( m \)) about center, axis perpendicular to rod

Let coordinate \( x \in [-L/2,\,L/2] \), uniform line density \( \lambda = m/L \), and distance to axis is \( r=|x| \):

\[ J = \int x^2\,dm = \int_{-L/2}^{L/2} x^2 \lambda\,dx = \lambda\left[\frac{x^3}{3}\right]_{-L/2}^{L/2} = \frac{m}{L}\cdot \frac{2}{3}\left(\frac{L}{2}\right)^3 = \frac{1}{12}mL^2. \]

Many rotational models need inertia about an axis that does not pass through the center of mass. The parallel-axis theorem provides a closed form.

Let \( J_C \): inertia about an axis through the center of mass, and consider a parallel axis displaced by distance \( d \). Then the inertia about the displaced axis is:

\[ J = J_C + m d^2. \]

Proof (continuous mass):

Let \( r \) be distance from the displaced axis, and \( r_C \) be distance from the center-of-mass axis. In planar geometry, the squared distance expands as:

\[ r^2 = r_C^2 + d^2 + 2d\,x, \]

where \( x \) is the signed coordinate of the mass element measured along the displacement direction (so that the cross term is linear in \( x \)). Integrate:

\[ J = \int r^2\,dm = \int r_C^2\,dm + d^2\int dm + 2d\int x\,dm. \]

The first term is \( J_C \), the second term is \( md^2 \). The third term vanishes because the center of mass is the origin for the \( x \)-coordinate: \( \int x\,dm = 0 \). Therefore:

\[ J = J_C + md^2. \]

Immediate application: A rod about one end (distance \( d=L/2 \) from its center axis) has:

\[ J_{\text{end}} = \frac{1}{12}mL^2 + m\left(\frac{L}{2}\right)^2 = \frac{1}{3}mL^2. \]

A linear torsional spring produces a restoring torque proportional to angular deflection: \( T_k(t) \): spring torque.

\[ T_k(t) = k_\theta\,\theta(t), \qquad [k_\theta] = \mathrm{N\cdot m/rad}. \]

Stored energy (proof by work integral):

The incremental work to twist by \( d\theta \) is \( dW = T_k\,d\theta \), hence:

\[ W(\theta) = \int_0^\theta k_\theta \vartheta\, d\vartheta = \frac{1}{2}k_\theta \theta^2. \]

Equivalent torsional stiffness: (direct rotational analog of Chapter 3, Lesson 3)

Parallel: same angle \( \theta \), torques add \( T = \sum T_i \):

\[ T = \sum_i k_{\theta i}\theta = \left(\sum_i k_{\theta i}\right)\theta \;\;\Rightarrow\;\; k_{\theta,\mathrm{eq}} = \sum_i k_{\theta i}. \]

Series: same torque \( T \), angles add \( \theta = \sum \theta_i \) with \( \theta_i = T/k_{\theta i} \):

\[ \theta = \sum_i \frac{T}{k_{\theta i}} = T\sum_i \frac{1}{k_{\theta i}} \;\;\Rightarrow\;\; k_{\theta,\mathrm{eq}} = \left(\sum_i \frac{1}{k_{\theta i}}\right)^{-1}. \]

Engineering formula: torsional stiffness of a uniform circular shaft

For a circular shaft of length \( L \) and shear modulus \( G \), linear torsion gives:

\[ T = \frac{G J_p}{L}\,\theta \;\;\Rightarrow\;\; k_\theta = \frac{G J_p}{L}, \]

where \( J_p \) is the polar second moment of area (a geometric property of the shaft cross-section). For a solid circular cross-section of radius \( R \):

\[ J_p = \int_A r^2\,dA = \int_0^{2\pi}\int_0^R r^2(r\,dr\,d\phi) = \frac{\pi R^4}{2}. \]

A viscous torsional damper produces a torque proportional to angular velocity: \( T_b(t) \): damping torque.

\[ T_b(t) = b_\theta\,\omega(t) = b_\theta\,\dot{\theta}(t), \qquad [b_\theta] = \mathrm{N\cdot m\cdot s/rad}. \]

Power dissipation (always nonnegative):

Damper power is \( P_b = T_b\omega \), so:

\[ P_b(t) = (b_\theta\omega)\omega = b_\theta \omega^2(t) \ge 0 \quad \text{for } b_\theta \ge 0. \]

Thus viscous damping removes mechanical energy monotonically (it cannot create energy).

Equivalent damping: The same series/parallel formulas as stiffness apply under the same structural assumptions:

• Parallel dampers: \( b_{\theta,\mathrm{eq}} = \sum_i b_{\theta i} \).
• Series dampers: \( b_{\theta,\mathrm{eq}} = \left(\sum_i \frac{1}{b_{\theta i}}\right)^{-1} \).

Consider a single inertia \( J \) connected to ground by a torsional spring \( k_\theta \) and torsional damper \( b_\theta \), driven by an external torque \( T_\mathrm{in}(t) \). Using torque balance from Section 3:

\[ J\ddot{\theta}(t) = T_\mathrm{in}(t) - b_\theta\dot{\theta}(t) - k_\theta\theta(t), \]

or equivalently the standard second-order ODE:

\[ J\ddot{\theta}(t) + b_\theta\dot{\theta}(t) + k_\theta\theta(t) = T_\mathrm{in}(t). \]

For numerical simulation, rewrite as a first-order system (without invoking state-space theory yet). Define \( x_1(t)=\theta(t) \) and \( x_2(t)=\dot{\theta}(t)=\omega(t) \). Then:

\[ \dot{x}_1(t) = x_2(t), \qquad \dot{x}_2(t) = \frac{1}{J}\left(T_\mathrm{in}(t) - b_\theta x_2(t) - k_\theta x_1(t)\right). \]

flowchart TD
  P1["Choose parameters: J, k, b"] --> P2["Specify input torque T_in(t)"]
  P2 --> P3["Write ODE: J*theta_ddot + b*theta_dot + k*theta = T_in"]
  P3 --> P4["Convert to first-order: x1=theta, x2=theta_dot"]
  P4 --> P5["Numerically integrate (RK4/ODE solver)"]
  P5 --> P6["Check units + sanity: energy, limiting cases"]
        

Sanity checks (highly recommended): (i) If \( k_\theta=0 \), the model becomes a damped rotational integrator. (ii) If \( b_\theta=0 \), energy is conserved in free motion. (iii) If \( J \to 0 \), the model tends to an algebraic constraint (nonphysical for rigid inertia, but useful for detecting modeling errors).

We simulate the torsional oscillator with a step input torque: \( T_\mathrm{in}(t)=T_0 \) for \( t \ge 0 \), with initial conditions \( \theta(0)=0 \), \( \omega(0)=0 \).

9.1 Python (SciPy + from-scratch RK4)

Core libraries used in system-dynamics simulation workflows: \( \texttt{NumPy} \) (arrays), \( \texttt{SciPy} \) (ODE solvers), and later (in future chapters) \( \texttt{python-control} \) for LTI analysis.

import numpy as np
from math import isfinite
from scipy.integrate import solve_ivp

# Parameters
J = 0.02          # kg*m^2
k = 1.5           # N*m/rad  (k_theta)
b = 0.08          # N*m*s/rad (b_theta)
T0 = 0.3          # N*m

def T_in(t):
    return T0  # step torque for t >= 0

def f(t, x):
    theta, omega = x
    dtheta = omega
    domega = (T_in(t) - b*omega - k*theta) / J
    return [dtheta, domega]

# Simulate using a robust adaptive solver
t_span = (0.0, 5.0)
x0 = [0.0, 0.0]
sol = solve_ivp(f, t_span, x0, max_step=1e-2, rtol=1e-8, atol=1e-10)

theta = sol.y[0]
omega = sol.y[1]
t = sol.t

# Quick sanity checks
assert np.all(np.isfinite(theta)) and np.all(np.isfinite(omega))
assert isfinite(theta[-1])

print("Final theta (rad):", theta[-1])
print("Final omega (rad/s):", omega[-1])

# From-scratch fixed-step RK4 (educational)
def rk4(f, t0, x0, tf, h):
    t = t0
    x = np.array(x0, dtype=float)
    ts = [t]
    xs = [x.copy()]
    while t < tf - 1e-12:
        k1 = np.array(f(t, x))
        k2 = np.array(f(t + 0.5*h, x + 0.5*h*k1))
        k3 = np.array(f(t + 0.5*h, x + 0.5*h*k2))
        k4 = np.array(f(t + h, x + h*k3))
        x = x + (h/6.0)*(k1 + 2*k2 + 2*k3 + k4)
        t = t + h
        ts.append(t)
        xs.append(x.copy())
    return np.array(ts), np.array(xs)

ts, xs = rk4(lambda tt, xx: f(tt, xx), 0.0, x0, 5.0, 1e-3)
print("RK4 final theta (rad):", xs[-1, 0])
      

9.2 C++ (from-scratch RK4)

Common C++ ecosystem choices for system dynamics: \( \texttt{Eigen} \) (linear algebra), \( \texttt{Boost::odeint} \) (ODE integration), plus custom solvers for embedded constraints.

#include <iostream>
#include <array>
#include <cmath>

struct Params {
  double J;
  double k;
  double b;
  double T0;
};

static inline double T_in(double /*t*/, const Params& p) {
  return p.T0; // step torque
}

static inline std::array<double, 2> f(double t, const std::array<double,2>& x, const Params& p) {
  const double theta = x[0];
  const double omega = x[1];
  const double dtheta = omega;
  const double domega = (T_in(t, p) - p.b*omega - p.k*theta) / p.J;
  return {dtheta, domega};
}

static inline std::array<double,2> add(const std::array<double,2>& a, const std::array<double,2>& b) {
  return {a[0] + b[0], a[1] + b[1]};
}

static inline std::array<double,2> mul(double s, const std::array<double,2>& a) {
  return {s*a[0], s*a[1]};
}

int main() {
  Params p{0.02, 1.5, 0.08, 0.3};
  double t = 0.0;
  const double tf = 5.0;
  const double h = 1e-3;

  std::array<double,2> x{0.0, 0.0}; // [theta, omega]

  while (t < tf - 1e-12) {
    auto k1 = f(t, x, p);
    auto k2 = f(t + 0.5*h, add(x, mul(0.5*h, k1)), p);
    auto k3 = f(t + 0.5*h, add(x, mul(0.5*h, k2)), p);
    auto k4 = f(t + h, add(x, mul(h, k3)), p);

    x = add(x, mul(h/6.0, add(add(k1, mul(2.0, k2)), add(mul(2.0, k3), k4))));
    t += h;
  }

  std::cout << "Final theta (rad): " << x[0] << "\n";
  std::cout << "Final omega (rad/s): " << x[1] << "\n";
  return 0;
}
      

9.3 Java (from-scratch RK4)

Common Java ecosystem choices: \( \texttt{Apache Commons Math} \) (ODE solvers, linear algebra) and \( \texttt{EJML} \) (matrices). Below is a minimal RK4 implementation to keep dependencies explicit.

public class TorsionalOscillatorRK4 {

  static class Params {
    double J, k, b, T0;
    Params(double J, double k, double b, double T0) {
      this.J = J; this.k = k; this.b = b; this.T0 = T0;
    }
  }

  static double torqueIn(double t, Params p) {
    return p.T0; // step torque
  }

  // x[0]=theta, x[1]=omega
  static double[] f(double t, double[] x, Params p) {
    double theta = x[0];
    double omega = x[1];
    double dtheta = omega;
    double domega = (torqueIn(t, p) - p.b*omega - p.k*theta) / p.J;
    return new double[]{dtheta, domega};
  }

  static double[] add(double[] a, double[] b) {
    return new double[]{a[0] + b[0], a[1] + b[1]};
  }

  static double[] mul(double s, double[] a) {
    return new double[]{s*a[0], s*a[1]};
  }

  public static void main(String[] args) {
    Params p = new Params(0.02, 1.5, 0.08, 0.3);

    double t = 0.0;
    double tf = 5.0;
    double h = 1e-3;

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

    while (t < tf - 1e-12) {
      double[] k1 = f(t, x, p);
      double[] k2 = f(t + 0.5*h, add(x, mul(0.5*h, k1)), p);
      double[] k3 = f(t + 0.5*h, add(x, mul(0.5*h, k2)), p);
      double[] k4 = f(t + h, add(x, mul(h, k3)), p);

      double[] incr = add(add(k1, mul(2.0, k2)), add(mul(2.0, k3), k4));
      x = add(x, mul(h/6.0, incr));
      t += h;
    }

    System.out.println("Final theta (rad): " + x[0]);
    System.out.println("Final omega (rad/s): " + x[1]);
  }
}
      

9.4 MATLAB + Simulink

MATLAB’s \( \texttt{ode45} \) is the standard entry point for continuous-time ODE simulation. Simulink provides block-based integration using Integrator blocks and algebraic sums/gains (useful for later block-diagram chapters).

% Parameters
J = 0.02;      % kg*m^2
k = 1.5;       % N*m/rad
b = 0.08;      % N*m*s/rad
T0 = 0.3;      % N*m

Tin = @(t) T0; % step torque

% First-order form: x1=theta, x2=omega
f = @(t,x) [ x(2);
            (Tin(t) - b*x(2) - k*x(1))/J ];

tspan = [0 5];
x0 = [0; 0];

opts = odeset('RelTol',1e-8,'AbsTol',1e-10);
[t,x] = ode45(f, tspan, x0, opts);

theta = x(:,1);
omega = x(:,2);

fprintf('Final theta (rad): %.6f\n', theta(end));
fprintf('Final omega (rad/s): %.6f\n', omega(end));
      

Simulink build notes (block-level):

• Use two Integrator blocks in series to integrate \( \ddot{\theta} \) to \( \dot{\theta} \) and then to \( \theta \).
• Compute \( \ddot{\theta} = (T_\mathrm{in} - b\dot{\theta} - k\theta)/J \) using Sum and Gain blocks.
• Input block: Step (amplitude \( T_0 \)) represents \( T_\mathrm{in} \).
• Scopes for \( \theta \) and \( \dot{\theta} \).

Optional: programmatic Simulink construction (MATLAB script)

model = 'torsional_oscillator';
new_system(model); open_system(model);

% Blocks
add_block('simulink/Sources/Step', [model '/Tin']);
add_block('simulink/Math Operations/Sum', [model '/SumTorque']);
add_block('simulink/Math Operations/Gain', [model '/GainInvJ']);
add_block('simulink/Math Operations/Gain', [model '/GainB']);
add_block('simulink/Math Operations/Gain', [model '/GainK']);
add_block('simulink/Continuous/Integrator', [model '/IntOmega']);
add_block('simulink/Continuous/Integrator', [model '/IntTheta']);
add_block('simulink/Sinks/Scope', [model '/ScopeTheta']);
add_block('simulink/Sinks/Scope', [model '/ScopeOmega']);

% Parameters
set_param([model '/GainInvJ'], 'Gain', num2str(1/0.02));
set_param([model '/GainB'], 'Gain', num2str(0.08));
set_param([model '/GainK'], 'Gain', num2str(1.5));
set_param([model '/Tin'], 'Time', '0', 'Before', '0', 'After', '0.3');

% Sum block signs: +Tin -b*omega -k*theta
set_param([model '/SumTorque'], 'Inputs', '+--');

% Wiring (positions omitted for brevity; can be arranged after adding)
add_line(model,'Tin/1','SumTorque/1');
add_line(model,'IntOmega/1','GainB/1');
add_line(model,'GainB/1','SumTorque/2');
add_line(model,'IntTheta/1','GainK/1');
add_line(model,'GainK/1','SumTorque/3');
add_line(model,'SumTorque/1','GainInvJ/1');
add_line(model,'GainInvJ/1','IntOmega/1');
add_line(model,'IntOmega/1','IntTheta/1');
add_line(model,'IntTheta/1','ScopeTheta/1');
add_line(model,'IntOmega/1','ScopeOmega/1');

set_param(model, 'StopTime', '5');
save_system(model);
      

9.5 Wolfram Mathematica (NDSolve)

Mathematica’s \( \texttt{NDSolve} \) handles stiff and non-stiff ODEs. For symbolic manipulation, \( \texttt{DSolve} \) can sometimes produce closed-form results for linear ODEs.

(* Parameters *)
J = 0.02;
k = 1.5;
b = 0.08;
T0 = 0.3;

Tin[t_] := T0;

(* ODE: J theta'' + b theta' + k theta == Tin(t) *)
sol = NDSolve[
  {
    J*theta''[t] + b*theta'[t] + k*theta[t] == Tin[t],
    theta[0] == 0,
    theta'[0] == 0
  },
  theta,
  {t, 0, 5}
];

thetaFinal = theta[5] /. sol[[1]];
omegaFinal = theta'[5] /. sol[[1]];

Print["Final theta (rad): ", thetaFinal];
Print["Final omega (rad/s): ", omegaFinal];
      

Problem 1 (Disk inertia): A uniform solid disk of mass \( m \) and radius \( R \) rotates about its symmetry axis. Derive \( J=\frac{1}{2}mR^2 \).

Solution: Use the surface density \( \sigma = m/(\pi R^2) \), with \( dm=\sigma r\,dr\,d\phi \):

\[ J=\int r^2\,dm=\int_0^{2\pi}\int_0^R r^2(\sigma r\,dr\,d\phi) =\sigma(2\pi)\frac{R^4}{4}=\frac{1}{2}mR^2. \]

Problem 2 (Parallel-axis theorem application): A slender rod (mass \( m \), length \( L \)) rotates about an axis through one end, perpendicular to the rod. Compute \( J \).

Solution: From Section 4, about the center: \( J_C=\frac{1}{12}mL^2 \). Distance to end axis is \( d=L/2 \). By Section 5:

\[ J = J_C + md^2 = \frac{1}{12}mL^2 + m\left(\frac{L}{2}\right)^2 = \frac{1}{3}mL^2. \]

Problem 3 (Equivalent torsional stiffness): Two torsional springs \( k_{\theta1} \) and \( k_{\theta2} \) are connected (a) in parallel and (b) in series. Derive \( k_{\theta,\mathrm{eq}} \).

Solution:

• Parallel: same \( \theta \), torques add: \( T=(k_{\theta1}+k_{\theta2})\theta \Rightarrow k_{\theta,\mathrm{eq}}=k_{\theta1}+k_{\theta2} \).
• Series: same \( T \), angles add: \( \theta=\frac{T}{k_{\theta1}}+\frac{T}{k_{\theta2}} \Rightarrow k_{\theta,\mathrm{eq}}=\left(\frac{1}{k_{\theta1}}+\frac{1}{k_{\theta2}}\right)^{-1} \).

Problem 4 (Energy identity and damping): For free motion \( T_\mathrm{in}(t)=0 \), define \( E(t)=\frac{1}{2}J\omega^2(t)+\frac{1}{2}k_\theta\theta^2(t) \). Show \( \dot{E}(t)=-b_\theta\omega^2(t)\le 0 \).

Solution: Differentiate and substitute the ODE.

\[ \dot{E}=J\omega\dot{\omega}+k_\theta\theta\dot{\theta} =\omega(J\dot{\omega}+k_\theta\theta). \]

From the free-motion ODE \( J\dot{\omega}+b_\theta\omega+k_\theta\theta=0 \), we get \( J\dot{\omega}+k_\theta\theta=-b_\theta\omega \). Hence:

\[ \dot{E}=\omega(-b_\theta\omega)=-b_\theta\omega^2 \le 0 \quad \text{for } b_\theta \ge 0. \]

Problem 5 (Steady deflection under constant torque): If \( T_\mathrm{in}(t)=T_0 \) is constant and the system reaches equilibrium, find the equilibrium angle \( \theta_\infty \).

Solution: At equilibrium, \( \dot{\theta}=0 \) and \( \ddot{\theta}=0 \), so the ODE reduces to:

\[ k_\theta \theta_\infty = T_0 \;\;\Rightarrow\;\; \theta_\infty = \frac{T_0}{k_\theta}. \]

We derived the fundamental rotational dynamics law \( \sum T = J\alpha \) from kinetic energy and power, defined the mass moment of inertia for discrete and continuous bodies, proved the parallel-axis theorem, and introduced linear torsional spring and viscous damper elements with energy and dissipation interpretations. These components combine into the canonical torsional oscillator ODE, which we implemented across five computational environments.