In Lessons 1–4, each physical domain was modeled with lumped storage and dissipation elements (e.g., inductance/capacitance/resistance; inertance/capacitance/resistance; thermal capacitance/resistance). Multi-domain systems arise when a variable in one domain appears in the constitutive equations of another domain, producing bidirectional interactions and shared power/energy pathways.

A multi-domain lumped model typically combines:

• \( \textbf{Storage laws} \) (dynamic): e.g., \( v = L \, \frac{di}{dt} \), \( q = C_h \, \frac{dp}{dt} \), \( \dot{Q} = C_{th} \, \frac{dT}{dt} \).
• \( \textbf{Dissipation laws} \) (algebraic): e.g., \( v = R i \), \( \Delta p = R_f q \), \( \dot{Q} = \frac{T - T_{amb}}{R_{th}} \).
• \( \textbf{Coupling laws} \) (cross-domain): e.g., back-EMF/torque constants, convective heat transport by flow, Joule heating.

flowchart LR
  E["Electrical subsystem\n(states: i, ... )"] --> C1["Electromechanical coupling\n'eb = ke*omega' and 'tau = kt*i'"]
  C1 --> M["Mechanical subsystem\n(states: omega, ... )"]

  F["Fluid subsystem\n(states: p, q, ... )"] --> C2["Thermo-fluid coupling\n'q moves enthalpy'\n'Tdot depends on q'"]
  C2 --> TH["Thermal subsystem\n(states: T, ... )"]

  E2["Electrical subsystem\n(states: i or v, ... )"] --> C3["Electro-thermal coupling\n'P = i^2*R(T)'\n'R depends on T'"]
  C3 --> TH2["Thermal subsystem\n(states: T, ... )"]
        

Throughout, we will validate models using power/energy consistency. For example, in the motor, the electrical power converted into mechanical form appears as a coupling term that cancels in the total energy derivative (it transfers energy but does not create or destroy it).

Consider a standard armature-controlled DC motor. The electrical side contains an inductance \( L_a \) and resistance \( R_a \), driven by an applied voltage \( v_a(t) \). The mechanical side has inertia \( J \) and viscous friction \( b \), with angular speed \( \omega(t) \) and load torque \( \tau_L(t) \).

Coupling laws (motor constants):

• Back electromotive force (back-EMF): \( e_b(t) = k_e \, \omega(t) \)
• Electromagnetic torque: \( \tau_m(t) = k_t \, i_a(t) \)

Electrical ODE (KVL):

\[ v_a(t) = L_a \frac{di_a}{dt} + R_a i_a(t) + k_e \omega(t). \]

Mechanical ODE (rotational Newton’s law):

\[ J \frac{d\omega}{dt} = k_t i_a(t) - b\,\omega(t) - \tau_L(t). \]

These two equations already form a closed electromechanical model with two dynamic states (current and angular speed). Writing them as a first-order system (without formal “state-space” terminology, which will be developed later), define \( x_1(t)= i_a(t) \) and \( x_2(t)=\omega(t) \).

\[ \begin{aligned} \dot{x}_1(t) &= -\frac{R_a}{L_a}x_1(t) - \frac{k_e}{L_a}x_2(t) + \frac{1}{L_a}v_a(t),\\ \dot{x}_2(t) &= \frac{k_t}{J}x_1(t) - \frac{b}{J}x_2(t) - \frac{1}{J}\tau_L(t). \end{aligned} \]

Energy consistency proof (storage + dissipation + supplied power):

Define stored energy as the sum of magnetic energy in the armature inductance and kinetic energy in the rotor: \( W(t) = \frac{1}{2}L_a i_a(t)^2 + \frac{1}{2}J \omega(t)^2 \). Differentiate:

\[ \dot{W}(t) = L_a i_a \frac{di_a}{dt} + J \omega \frac{d\omega}{dt}. \]

Substitute the ODEs (multiply electrical equation by \( i_a \), mechanical equation by \( \omega \)):

\[ \begin{aligned} i_a v_a &= L_a i_a \frac{di_a}{dt} + R_a i_a^2 + k_e \omega i_a,\\ J \omega \frac{d\omega}{dt} &= k_t i_a \omega - b\omega^2 - \tau_L \omega. \end{aligned} \]

Add the two equalities and rearrange using \( \dot{W} = L_a i_a \frac{di_a}{dt} + J \omega \frac{d\omega}{dt} \):

\[ \dot{W}(t) = i_a v_a - R_a i_a^2 - b \omega^2 - \tau_L \omega + (k_t - k_e)\, i_a \omega. \]

In consistent SI motor modeling (with compatible definitions of \( k_t \) and \( k_e \)), we have \( k_t = k_e \), and the coupling term cancels:

\[ \boxed{ \dot{W}(t) = i_a(t)\,v_a(t) - R_a i_a(t)^2 - b\,\omega(t)^2 - \tau_L(t)\,\omega(t). } \]

Interpretation: electrical input power \( i_a v_a \) is partitioned into stored energy rate \( \dot{W} \), resistive loss \( R_a i_a^2 \), mechanical friction loss \( b\omega^2 \), and mechanical output power to the load \( \tau_L \omega \). This is a rigorous internal consistency check for the coupled model.

Special case (unforced, no load): if \( v_a(t)=0 \) and \( \tau_L(t)=0 \), then \( \dot{W}(t) = -R_a i_a^2 - b\omega^2 \le 0 \), so stored energy is non-increasing and the system dissipates toward rest (a stability insight derived purely from physics).

We combine a (lumped) hydraulic subsystem and a (lumped) thermal subsystem by recognizing that fluid flow transports energy. A standard modeling assumption (appropriate here) is a well-mixed tank temperature \( T(t) \).

Hydraulic subsystem (Lessons 3 concepts):

• Flow inertance \( I_f \) and resistance \( R_f \) in the feed line
• Tank hydraulic capacitance \( C_h \) relating tank pressure \( p(t) \) to net inflow
• Outlet hydraulic resistance \( R_o \) so \( q_{out}(t)=\frac{p(t)}{R_o} \)

Let \( q(t) \) be the inflow through the line and \( p_s(t) \) a supply pressure. The line pressure drop satisfies:

\[ p_s(t) - p(t) = I_f \frac{dq}{dt} + R_f q(t). \]

The tank capacitance equation is the flow balance into the capacitance element:

\[ C_h \frac{dp}{dt} = q(t) - \frac{p(t)}{R_o}. \]

Thermal subsystem (Lesson 4 concepts):

Let \( C_{th} \) be the lumped thermal capacitance of the tank contents (units J/K) and \( R_{th} \) the thermal resistance to ambient temperature \( T_{amb} \). An external heater supplies heat rate \( \dot{Q}_h(t) \). Inflow has temperature \( T_{in} \) (assumed constant for this lesson).

The key coupling: inflow \( q(t) \) transports enthalpy. Under constant density \( \rho \) and specific heat \( c_p \), the net convective heat flow into a well-mixed volume is \( \rho c_p q(t)\,(T_{in}-T(t)) \). Thus:

\[ C_{th}\frac{dT}{dt} = \rho c_p q(t)\bigl(T_{in}-T(t)\bigr) - \frac{T(t)-T_{amb}}{R_{th}} + \dot{Q}_h(t). \]

This model is coupled because \( q(t) \) is not a prescribed input but is produced by the hydraulic dynamics through \( (p,q) \). Collecting, we obtain the closed ODE system in \( (q,p,T) \):

\[ \begin{aligned} \dot{q} &= \frac{1}{I_f}\Bigl(p_s(t) - p - R_f q\Bigr),\\ \dot{p} &= \frac{1}{C_h}\Bigl(q - \frac{p}{R_o}\Bigr),\\ \dot{T} &= \frac{\rho c_p}{C_{th}}\,q\,(T_{in}-T) - \frac{1}{R_{th}C_{th}}\,(T-T_{amb}) + \frac{1}{C_{th}}\dot{Q}_h(t). \end{aligned} \]

Proof of uniqueness and boundedness for fixed flow (a useful sub-result):

If \( q(t)=\bar{q} \ge 0 \) is constant and \( \dot{Q}_h(t)=\bar{Q} \) is constant, then the temperature equation becomes linear time-invariant:

\[ \dot{T} + a\,T = b, \quad a = \frac{\rho c_p}{C_{th}}\bar{q} + \frac{1}{R_{th}C_{th}}, \quad b = \frac{\rho c_p}{C_{th}}\bar{q}\,T_{in} + \frac{1}{R_{th}C_{th}}T_{amb} + \frac{\bar{Q}}{C_{th}}. \]

Since \( a > 0 \) when \( C_{th} > 0 \), \( R_{th} > 0 \), and \( \bar{q}\ge 0 \), the unique solution is:

\[ T(t) = T(0)\,e^{-a t} + \frac{b}{a}\bigl(1-e^{-a t}\bigr), \quad \text{so } T(t) \; → \; \frac{b}{a} \text{ as } t \to \infty. \]

This establishes (in a fully elementary way) existence, uniqueness, and convergence for the thermal sub-dynamics under constant flow. In the fully coupled model, \( q(t) \) becomes dynamic and the thermal equation becomes bilinear; nevertheless, the structure above remains valuable for sanity-checks and parameter intuition.

Electro-thermal coupling appears whenever electrical losses heat a component and its temperature feeds back into its electrical parameters. A minimal lumped model uses:

• Electrical algebraic law (Ohm): \( v(t)=i(t)\,R(T(t)) \)
• Thermal dynamics: \( C_{th}\frac{dT}{dt} = P_{J}(t) - \frac{T-T_{amb}}{R_{th}} \)
• Joule power: \( P_{J}(t)=i(t)^2 R(T(t)) \)

A common first-order temperature coefficient model is: \( R(T)=R_0\bigl(1+\alpha(T-T_0)\bigr) \), where \( \alpha \) can be positive (PTC, typical metals) or negative (NTC-like behavior).

If we drive the resistor by a prescribed voltage \( v(t)=\bar{v} \) (constant), then \( i(t)=\frac{\bar{v}}{R(T(t))} \), and the thermal ODE becomes a single nonlinear ODE in \( T(t) \):

\[ C_{th}\frac{dT}{dt} = \frac{\bar{v}^2}{R(T)} - \frac{T-T_{amb}}{R_{th}} = \frac{\bar{v}^2}{R_0(1+\alpha(T-T_0))} - \frac{T-T_{amb}}{R_{th}}. \]

Equilibrium temperature: solve \( \frac{dT}{dt}=0 \), i.e.

\[ \frac{T^\star - T_{amb}}{R_{th}} = \frac{\bar{v}^2}{R_0(1+\alpha(T^\star-T_0))}. \]

This can be rearranged into a quadratic in \( T^\star \). Define \( \Delta = T^\star - T_0 \) and \( \Delta_{amb} = T_{amb}-T_0 \). Then \( T^\star - T_{amb} = \Delta - \Delta_{amb} \) and:

\[ (\Delta-\Delta_{amb})(1+\alpha\Delta) = \frac{R_{th}\bar{v}^2}{R_0}. \]

Expanding:

\[ \alpha \Delta^2 + (1-\alpha \Delta_{amb})\Delta - \Delta_{amb} - \frac{R_{th}\bar{v}^2}{R_0} = 0. \]

Local stability test (scalar ODE):

Write the scalar dynamics as \( \dot{T} = f(T) \). For a differentiable scalar ODE, an equilibrium \( T^\star \) is locally asymptotically stable if \( f'(T^\star) < 0 \). Here:

\[ f(T) = \frac{1}{C_{th}}\left(\frac{\bar{v}^2}{R_0(1+\alpha(T-T_0))} - \frac{T-T_{amb}}{R_{th}}\right). \]

Differentiate:

\[ f'(T) = \frac{1}{C_{th}}\left( -\frac{\bar{v}^2}{R_0}\cdot\frac{\alpha}{(1+\alpha(T-T_0))^2} -\frac{1}{R_{th}} \right). \]

If \( \alpha \ge 0 \), then the first term is non-positive, so \( f'(T) \le -\frac{1}{C_{th}R_{th}} < 0 \). Therefore, any equilibrium that satisfies \( 1+\alpha(T^\star-T_0) > 0 \) is locally asymptotically stable. If \( \alpha < 0 \), the electrical heating term contributes positive feedback and instability/thermal runaway can occur depending on parameters.

Across the three examples, a repeatable modeling pattern emerges:

1. Choose domain state variables from storage elements: currents (inductors), voltages (capacitors), pressures (hydraulic capacitances), flows (hydraulic inertances), temperatures (thermal capacitances), velocities (mechanical inertias).
2. Write conservation laws (KVL/KCL analogs, flow continuity, Newton’s law, thermal balance).
3. Insert coupling laws that relate variables across domains.
4. Eliminate purely algebraic internal variables (e.g., back-EMF, Joule power) to obtain a closed ODE.
5. Validate using sign/units checks and (when possible) an energy/dissipation identity.

Implementations in any language follow the same computational workflow: define parameters, implement the right-hand side of the ODE, integrate numerically, and validate against physics-based invariants (e.g., energy balance) and limiting cases.

flowchart TD
  S["Start: define subsystem elements\n(R, L, C, inertance, thermal C/R, J, b)"] --> X["Choose dynamic variables\n(i, omega, p, q, T)"]
  X --> EQ["Write equations\n(conservation + constitutive + coupling)"]
  EQ --> CL["Close the model\n(eliminate algebraic vars, define inputs)"]
  CL --> CK["Consistency checks\n(units, signs, limiting cases)"]
  CK --> NUM["Numerical integration\n(Euler/RK4/ODE solvers)"]
  NUM --> VAL["Validation\n(energy/dissipation, steady-state, sanity tests)"]
  VAL --> REF["Refine parameters/assumptions\n(if mismatch)"]
        

Recommended libraries for system dynamics workflows in Python include: numpy (arrays), scipy.integrate (ODE solvers), and sympy (symbolic checks). (Control-specific libraries will be used more heavily once transfer functions and state-space are introduced later.)

import numpy as np
from scipy.integrate import solve_ivp

# ----------------------------
# 1) Electromechanical DC motor
# ----------------------------
def dc_motor_rhs(t, x, p):
    i, w = x
    La, Ra, ke, kt, J, b = p["La"], p["Ra"], p["ke"], p["kt"], p["J"], p["b"]
    va = p["va"](t)
    tauL = p["tauL"](t)

    di = (va - Ra*i - ke*w) / La
    dw = (kt*i - b*w - tauL) / J
    return [di, dw]

# ----------------------------
# 2) Thermo-fluid (q, p, T)
# ----------------------------
def thermofluid_rhs(t, x, p):
    q, pres, T = x
    If, Rf, Ch, Ro = p["If"], p["Rf"], p["Ch"], p["Ro"]
    rho, cp, Cth, Rth = p["rho"], p["cp"], p["Cth"], p["Rth"]
    Tin, Tamb = p["Tin"], p["Tamb"]

    ps = p["ps"](t)
    Qh = p["Qh"](t)

    dq = (ps - pres - Rf*q) / If
    dpres = (q - pres/Ro) / Ch
    dT = (rho*cp*q*(Tin - T) - (T - Tamb)/Rth + Qh) / Cth
    return [dq, dpres, dT]

# ----------------------------
# 3) Electro-thermal (T only), voltage driven
# ----------------------------
def electrothermal_rhs(t, x, p):
    T = x[0]
    R0, alpha, T0 = p["R0"], p["alpha"], p["T0"]
    Cth, Rth, Tamb = p["Cth"], p["Rth"], p["Tamb"]
    v = p["v"](t)

    R = R0*(1.0 + alpha*(T - T0))
    # guard against invalid R (e.g., negative when alpha<0 and T too high)
    if R <= 1e-9:
        R = 1e-9
    dT = (v*v/R - (T - Tamb)/Rth) / Cth
    return [dT]

# ----------------------------
# Example runs
# ----------------------------
if __name__ == "__main__":
    # DC motor parameters
    pm = {
        "La": 5e-3, "Ra": 1.2, "ke": 0.08, "kt": 0.08,
        "J": 2e-4, "b": 1e-4,
        "va": lambda t: 12.0*(t >= 0.05),   # step at 0.05 s
        "tauL": lambda t: 0.0
    }
    sol_m = solve_ivp(lambda t, x: dc_motor_rhs(t, x, pm), [0, 0.5], [0.0, 0.0], max_step=1e-3)

    # Thermo-fluid parameters
    ptf = {
        "If": 2.0e3, "Rf": 8.0e7, "Ch": 2.5e-10, "Ro": 1.2e8,  # illustrative
        "rho": 1000.0, "cp": 4180.0, "Cth": 3.0e4, "Rth": 0.08,
        "Tin": 330.0, "Tamb": 295.0,
        "ps": lambda t: 2.0e5,                # constant supply pressure
        "Qh": lambda t: 500.0*(t >= 10.0)     # heater step at 10 s
    }
    sol_tf = solve_ivp(lambda t, x: thermofluid_rhs(t, x, ptf), [0, 60], [0.0, 0.0, 295.0], max_step=0.05)

    # Electro-thermal parameters
    pet = {
        "R0": 10.0, "alpha": 0.004, "T0": 293.15,
        "Cth": 40.0, "Rth": 2.0, "Tamb": 293.15,
        "v": lambda t: 12.0
    }
    sol_et = solve_ivp(lambda t, x: electrothermal_rhs(t, x, pet), [0, 200], [293.15], max_step=0.2)

    # Basic numerical outputs (plotting intentionally omitted in HTML lesson content)
    print("Motor final current, speed:", sol_m.y[0, -1], sol_m.y[1, -1])
    print("Thermo-fluid final q, p, T:", sol_tf.y[0, -1], sol_tf.y[1, -1], sol_tf.y[2, -1])
    print("Electro-thermal final T:", sol_et.y[0, -1])
      

Suggested verification exercises: (i) check the motor energy balance numerically by computing \( \frac{d}{dt}\left(\frac{1}{2}L_ai^2+\frac{1}{2}J\omega^2\right) \) and comparing it to \( i v_a - R_ai^2 - b\omega^2 - \tau_L \omega \); (ii) set \( R_{th}\to \infty \) in the electro-thermal model to confirm temperature grows as dictated purely by Joule heating.

Common C++ libraries relevant to system dynamics work include: Eigen (linear algebra), boost::numeric::odeint (ODE solvers), and plotting via external tools. Below is a self-contained RK4 integrator for the DC motor ODEs.

#include <iostream>
#include <vector>
#include <functional>

struct Params {
  double La, Ra, ke, kt, J, b;
};

static std::vector<double> rhs_motor(double t,
                                     const std::vector<double>& x,
                                     const Params& p,
                                     const std::function<double(double)>& va,
                                     const std::function<double(double)>& tauL) {
  double i = x[0];
  double w = x[1];
  double di = (va(t) - p.Ra*i - p.ke*w) / p.La;
  double dw = (p.kt*i - p.b*w - tauL(t)) / p.J;
  return {di, dw};
}

static std::vector<double> rk4_step(double t, double h,
                                   const std::vector<double>& x,
                                   const std::function<std::vector<double>(double,const std::vector<double>&)>& f) {
  auto k1 = f(t, x);

  std::vector<double> x2(x.size());
  for (size_t i=0;i<x.size();++i) x2[i] = x[i] + 0.5*h*k1[i];
  auto k2 = f(t + 0.5*h, x2);

  std::vector<double> x3(x.size());
  for (size_t i=0;i<x.size();++i) x3[i] = x[i] + 0.5*h*k2[i];
  auto k3 = f(t + 0.5*h, x3);

  std::vector<double> x4(x.size());
  for (size_t i=0;i<x.size();++i) x4[i] = x[i] + h*k3[i];
  auto k4 = f(t + h, x4);

  std::vector<double> xn(x.size());
  for (size_t i=0;i<x.size();++i) {
    xn[i] = x[i] + (h/6.0)*(k1[i] + 2.0*k2[i] + 2.0*k3[i] + k4[i]);
  }
  return xn;
}

int main() {
  Params p{5e-3, 1.2, 0.08, 0.08, 2e-4, 1e-4};

  auto va = [](double t){ return (t >= 0.05) ? 12.0 : 0.0; };
  auto tauL = [](double /*t*/){ return 0.0; };

  double t0 = 0.0, tf = 0.5, h = 1e-4;
  std::vector<double> x{0.0, 0.0}; // i, omega

  auto f = [&](double t, const std::vector<double>& xcur){
    return rhs_motor(t, xcur, p, va, tauL);
  };

  int steps = static_cast<int>((tf - t0)/h);
  double t = t0;
  for (int k=0;k<steps;++k) {
    x = rk4_step(t, h, x, f);
    t += h;
    if (k % 1000 == 0) {
      std::cout << "t=" << t << "  i=" << x[0] << "  omega=" << x[1] << "\n";
    }
  }
  return 0;
}
      

For larger coupled systems, replace the hand-coded integrator with boost::numeric::odeint and manage vectors/matrices using Eigen. Energy-balance checks (as derived in Section 2) remain strongly recommended regardless of implementation choice.

Java libraries relevant to system dynamics include Apache Commons Math (ODE integration, linear algebra), and EJML (matrix computation). Below is a minimal from-scratch RK4 for the voltage-driven electro-thermal ODE.

import java.util.function.DoubleUnaryOperator;

public class ElectroThermalRK4 {

  static class Params {
    double R0, alpha, T0;
    double Cth, Rth, Tamb;
    DoubleUnaryOperator v;
  }

  static double rhs(double t, double T, Params p) {
    double R = p.R0 * (1.0 + p.alpha * (T - p.T0));
    if (R <= 1e-9) R = 1e-9; // guard
    double v = p.v.applyAsDouble(t);
    return (v*v / R - (T - p.Tamb)/p.Rth) / p.Cth;
  }

  static double rk4Step(double t, double h, double T, Params p) {
    double k1 = rhs(t, T, p);
    double k2 = rhs(t + 0.5*h, T + 0.5*h*k1, p);
    double k3 = rhs(t + 0.5*h, T + 0.5*h*k2, p);
    double k4 = rhs(t + h, T + h*k3, p);
    return T + (h/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4);
  }

  public static void main(String[] args) {
    Params p = new Params();
    p.R0 = 10.0;
    p.alpha = 0.004;     // PTC example
    p.T0 = 293.15;
    p.Cth = 40.0;
    p.Rth = 2.0;
    p.Tamb = 293.15;
    p.v = (t) -> 12.0;

    double t = 0.0, tf = 200.0, h = 0.01;
    double T = 293.15;

    int steps = (int)((tf - t)/h);
    for (int k=0;k<steps;k++) {
      T = rk4Step(t, h, T, p);
      t += h;
      if (k % 5000 == 0) {
        System.out.printf("t=%.2f  T=%.4f K%n", t, T);
      }
    }
  }
}
      

Extension exercise: replace constant v with a step or sinusoid and compare transient temperature responses. In later chapters, you can connect this thermal state to other subsystems (e.g., temperature-dependent motor winding resistance).

MATLAB toolchain options include: ode45 for ODE simulation, Simulink for block-diagram simulation, and (optionally) Simscape for multi-domain physical networks (introduced as a tooling option, not as a new theory).

MATLAB script (ODE45):

function thermoFluidDemo()
  % Parameters
  p.If  = 2.0e3;
  p.Rf  = 8.0e7;
  p.Ch  = 2.5e-10;
  p.Ro  = 1.2e8;

  p.rho = 1000.0;
  p.cp  = 4180.0;
  p.Cth = 3.0e4;
  p.Rth = 0.08;

  p.Tin  = 330.0;
  p.Tamb = 295.0;

  % Inputs
  ps = @(t) 2.0e5;
  Qh = @(t) (t >= 10.0)*500.0;

  % ODE
  f = @(t,x) rhs(t,x,p,ps,Qh);

  % Initial conditions: q, p, T
  x0 = [0.0; 0.0; 295.0];

  % Solve
  tspan = [0 60];
  opts = odeset('MaxStep',0.05);
  [t,x] = ode45(f,tspan,x0,opts);

  % Display final values
  disp(['Final q = ', num2str(x(end,1))]);
  disp(['Final p = ', num2str(x(end,2))]);
  disp(['Final T = ', num2str(x(end,3))]);

  % Optional plotting (can be enabled in your environment)
  % figure; plot(t,x(:,3)); xlabel('t'); ylabel('T');
end

function dx = rhs(t,x,p,ps,Qh)
  q = x(1); pres = x(2); T = x(3);

  dq    = (ps(t) - pres - p.Rf*q)/p.If;
  dpres = (q - pres/p.Ro)/p.Ch;
  dT    = (p.rho*p.cp*q*(p.Tin - T) - (T - p.Tamb)/p.Rth + Qh(t))/p.Cth;

  dx = [dq; dpres; dT];
end
      

Simulink block outline (no new theory required):

• Use three Integrator blocks to represent \( q \), \( p \), \( T \).
• Implement \( \dot{q}=(p_s-p-R_f q)/I_f \) with Sum and Gain blocks.
• Implement \( \dot{p}=(q-p/R_o)/C_h \) similarly.
• Implement \( \dot{T}=(\rho c_p q (T_{in}-T) - (T-T_{amb})/R_{th} + \dot{Q}_h)/C_{th} \) using Product blocks (for the coupling term \( q(T_{in}-T) \)).

Mathematica is particularly effective for symbolic checks (e.g., differentiating the energy function) and direct numerical integration with NDSolve.

(* DC motor parameters *)
La = 5*10^-3; Ra = 1.2; ke = 0.08; kt = 0.08;
J  = 2*10^-4; b  = 1*10^-4;

va[t_] := If[t >= 0.05, 12.0, 0.0];
tauL[t_] := 0.0;

eqs = {
  i'[t] == (va[t] - Ra*i[t] - ke*ω[t])/La,
  ω'[t] == (kt*i[t] - b*ω[t] - tauL[t])/J,
  i[0] == 0.0,
  ω[0] == 0.0
};

sol = NDSolve[eqs, {i, ω}, {t, 0, 0.5}];

(* Energy function *)
W[t_] := 1/2 La (i[t] /. sol[[1]])^2 + 1/2 J (ω[t] /. sol[[1]])^2;

(* Optional: evaluate final states and energy *)
finalI = i[0.5] /. sol[[1]];
finalW = W[0.5];
{finalI, finalW}
      

Extension exercise: define \( \tau_L(t) \) as a step load torque and observe the trade-off between speed regulation and armature current.

Problem 1 (DC Motor Coupled ODE Derivation): Starting from KVL and rotational dynamics, derive the coupled first-order ODEs for \( i_a(t) \) and \( \omega(t) \) given \( e_b = k_e \omega \) and \( \tau_m = k_t i_a \).

Solution:

KVL gives \( v_a = L_a \frac{di_a}{dt} + R_a i_a + k_e \omega \), hence \( \frac{di_a}{dt} = \frac{1}{L_a}(v_a - R_a i_a - k_e \omega) \). Rotational dynamics gives \( J\frac{d\omega}{dt} = k_t i_a - b\omega - \tau_L \), hence \( \frac{d\omega}{dt} = \frac{1}{J}(k_t i_a - b\omega - \tau_L) \).

Problem 2 (Energy Balance and Dissipation Inequality): For the DC motor with \( k_t=k_e \), show that with \( v_a(t)=0 \) and \( \tau_L(t)=0 \), the energy \( W(t)=\frac{1}{2}L_a i_a^2 + \frac{1}{2}J\omega^2 \) satisfies \( \dot{W}(t) \le 0 \).

Solution:

From Section 2, when \( k_t=k_e \) we have:

\[ \dot{W}(t) = i_a v_a - R_a i_a^2 - b\omega^2 - \tau_L \omega. \]

Setting \( v_a=0 \) and \( \tau_L=0 \) yields \( \dot{W}(t) = -R_a i_a^2 - b\omega^2 \le 0 \) since \( R_a > 0 \) and \( b > 0 \). Therefore energy decreases monotonically (or stays constant only at \( i_a=\omega=0 \)).

Problem 3 (Thermo-Fluid Steady-State Temperature Under Constant Flow): Assume constant flow \( q(t)=\bar{q}\ge 0 \) and constant heater power \( \dot{Q}_h(t)=\bar{Q} \). Derive the steady-state temperature \( T^\star \) for: \( C_{th}\dot{T}=\rho c_p \bar{q}(T_{in}-T) - \frac{T-T_{amb}}{R_{th}} + \bar{Q} \).

Solution:

At steady-state, \( \dot{T}=0 \) so:

\[ 0 = \rho c_p \bar{q}(T_{in}-T^\star) - \frac{T^\star-T_{amb}}{R_{th}} + \bar{Q}. \]

Collect terms in \( T^\star \):

\[ \left(\rho c_p \bar{q} + \frac{1}{R_{th}}\right)T^\star = \rho c_p \bar{q} \, T_{in} + \frac{1}{R_{th}}T_{amb} + \bar{Q}. \]

Therefore:

\[ \boxed{ T^\star = \frac{\rho c_p \bar{q} \, T_{in} + \frac{1}{R_{th}}T_{amb} + \bar{Q}} {\rho c_p \bar{q} + \frac{1}{R_{th}}}. } \]

Problem 4 (Electro-Thermal Equilibrium and Stability for PTC): Consider the voltage-driven electro-thermal model \( C_{th}\dot{T}=\frac{\bar{v}^2}{R_0(1+\alpha(T-T_0))}-\frac{T-T_{amb}}{R_{th}} \) with \( \alpha \ge 0 \). Show that any equilibrium satisfying \( 1+\alpha(T^\star-T_0) > 0 \) is locally asymptotically stable.

Solution:

Define \( \dot{T}=f(T) \). From Section 4:

\[ f'(T) = \frac{1}{C_{th}}\left( -\frac{\bar{v}^2}{R_0}\cdot\frac{\alpha}{(1+\alpha(T-T_0))^2} -\frac{1}{R_{th}} \right). \]

If \( \alpha \ge 0 \) and \( 1+\alpha(T^\star-T_0) > 0 \), then the fraction is nonnegative and the first term is non-positive. Hence: \( f'(T^\star) \le -\frac{1}{C_{th}R_{th}} < 0 \), proving local asymptotic stability.

Problem 5 (Dimension/Units Consistency Check): Verify that each term in the motor energy identity \( \dot{W} = i_a v_a - R_a i_a^2 - b\omega^2 - \tau_L\omega \) has units of power.

Solution:

• \( i_a v_a \): A·V = W (watts).
• \( R_a i_a^2 \): (V/A)·A² = V·A = W.
• \( b\omega^2 \): since \( b\omega \) is torque (N·m), multiply by \( \omega \) (rad/s) gives N·m/s = W.
• \( \tau_L\omega \): torque (N·m) times angular speed (rad/s) gives W.

We constructed three canonical coupled multi-domain lumped models and showed how coupling laws (back-EMF/torque, convective transport, and Joule heating with parameter feedback) convert single-domain element equations into closed ODE systems. We emphasized an essential validation tool: energy/dissipation identities, which provide rigorous internal consistency checks for coupled models. Implementations in Python, C++, Java, MATLAB/Simulink, and Mathematica demonstrated that multi-domain dynamics are handled computationally by the same ODE integration workflow used for single-domain systems.