In the thermal domain, the natural effort variable is temperature \( T(t) \), and the natural flow variable is heat rate (power) \( \dot{Q}(t) \). These choices align with the electrical/hydraulic analogies already introduced earlier in this chapter: temperature differences drive heat flow through thermal resistances, and thermal capacitances store internal energy.

The main heat-transfer mechanisms used in lumped network models are:

Conduction (Fourier’s law, 1D):

\[ \dot{Q} = -kA \frac{dT}{dx} \]

For a uniform slab of thickness \( L \), constant conductivity \( k \), and cross-sectional area \( A \), steady conduction yields the equivalent thermal resistance \( R_{\mathrm{cond}} \):

\[ \dot{Q} = \frac{T_1 - T_2}{R_{\mathrm{cond}}},\qquad R_{\mathrm{cond}} = \frac{L}{kA} \]

Convection (Newton’s law of cooling):

\[ \dot{Q} = hA\,(T_s - T_{\infty}) \;=\; \frac{T_s - T_{\infty}}{R_{\mathrm{conv}}},\qquad R_{\mathrm{conv}} = \frac{1}{hA} \]

Here \( h \) is the convection coefficient, \( A \) is surface area, \( T_s \) is surface temperature, and \( T_{\infty} \) is the ambient fluid temperature.

Thermal capacitance (energy storage):

If a body can be treated as isothermal (lumped), its stored thermal energy relative to a reference is approximately linear in temperature:

\[ E(T) = \int m c_p \, dT \;=\; m c_p T + \text{constant} \]

Therefore the thermal capacitance is \( C_{\mathrm{th}} = m c_p \), and conservation of energy gives

\[ \dot{E}(t) = C_{\mathrm{th}} \frac{dT}{dt} = \dot{Q}_{\mathrm{in}}(t) - \dot{Q}_{\mathrm{out}}(t) \]

Lumped modeling assumes the body has negligible internal temperature gradients, i.e. \( T(\mathbf{x},t) \approx T(t) \). A standard dimensionless test compares internal conduction resistance to boundary convection resistance.

Define the characteristic length \( L_c \) by \( L_c = \frac{V}{A_s} \), where \( V \) is volume and \( A_s \) is the convecting surface area. The Biot number is:

\[ \mathrm{Bi} = \frac{h L_c}{k} \]

Interpretation: \( \mathrm{Bi} \) is approximately the ratio \( R_{\mathrm{cond,internal}} / R_{\mathrm{conv,boundary}} \). A commonly used condition for lumped validity is \( \mathrm{Bi} \lt 0.1 \).

Derivation sketch (dimensionless resistance ratio):

Model internal conduction by a resistance scale \( R_{\mathrm{cond,int}} \sim \frac{L_c}{kA_s} \), and convection by \( R_{\mathrm{conv}} = \frac{1}{hA_s} \). Then

\[ \frac{R_{\mathrm{cond,int}}}{R_{\mathrm{conv}}} \sim \frac{\frac{L_c}{kA_s}}{\frac{1}{hA_s}} = \frac{h L_c}{k} = \mathrm{Bi}. \]

Thus \( \mathrm{Bi} \ll 1 \) implies convection dominates, and internal gradients are small, supporting the lumped approximation.

Consider a lumped body with temperature \( T(t) \), thermal capacitance \( C_{\mathrm{th}} \), exchanging heat with an ambient node at temperature \( T_{\infty}(t) \) through an equivalent thermal resistance \( R_{\mathrm{th}} \). Let an external heat input (heater, dissipation, etc.) be \( \dot{Q}_u(t) \).

The heat leaving the body through the resistance is \( \dot{Q}_{\mathrm{loss}}(t) = \frac{T(t) - T_{\infty}(t)}{R_{\mathrm{th}}} \). Energy balance yields:

\[ C_{\mathrm{th}}\frac{dT}{dt} = \dot{Q}_u(t) - \frac{T(t) - T_{\infty}(t)}{R_{\mathrm{th}}}. \]

Rearranging into a standard first-order linear ODE:

\[ \frac{dT}{dt} + \frac{1}{R_{\mathrm{th}}C_{\mathrm{th}}}T = \frac{1}{R_{\mathrm{th}}C_{\mathrm{th}}}T_{\infty}(t) + \frac{1}{C_{\mathrm{th}}}\dot{Q}_u(t). \]

The natural time constant is \( \tau = R_{\mathrm{th}}C_{\mathrm{th}} \). For the important special case of constant ambient \( T_{\infty}(t)=T_{\infty} \) and constant input \( \dot{Q}_u(t)=\dot{Q}_0 \), the solution is:

\[ T(t) = T_{\infty} + R_{\mathrm{th}}\dot{Q}_0 + \Big(T(0)-\big(T_{\infty}+R_{\mathrm{th}}\dot{Q}_0\big)\Big)\, e^{-t/(R_{\mathrm{th}}C_{\mathrm{th}})}. \]

Proof (integrating factor):

With constant \( T_{\infty} \) and \( \dot{Q}_0 \), write the ODE as \( \frac{dT}{dt} + \frac{1}{\tau}T = \frac{1}{\tau}(T_{\infty}+R_{\mathrm{th}}\dot{Q}_0) \). The integrating factor is \( \mu(t)=e^{t/\tau} \). Multiply both sides:

\[ e^{t/\tau}\frac{dT}{dt} + \frac{1}{\tau}e^{t/\tau}T = \frac{1}{\tau}e^{t/\tau}(T_{\infty}+R_{\mathrm{th}}\dot{Q}_0). \]

The left-hand side is \( \frac{d}{dt}\big(e^{t/\tau}T(t)\big) \), hence

\[ \frac{d}{dt}\big(e^{t/\tau}T(t)\big) = \frac{1}{\tau}e^{t/\tau}(T_{\infty}+R_{\mathrm{th}}\dot{Q}_0). \]

Integrate from \( 0 \) to \( t \), solve for \( T(t) \), and obtain the stated exponential form.

Thermal resistances combine exactly like electrical resistances under the temperature–heatflow analogy. The equivalences are not mere analogies; they follow from the same algebraic structure of conservation laws.

Series combination: Two resistances \( R_1 \) and \( R_2 \) in series carry the same heat flow \( \dot{Q} \), and temperature drops add:

\[ T_1 - T_3 = (T_1 - T_2) + (T_2 - T_3) = \dot{Q}R_1 + \dot{Q}R_2 = \dot{Q}(R_1+R_2). \]

Hence the equivalent resistance is:

\[ R_{\mathrm{eq}} = R_1 + R_2. \]

Parallel combination: Two resistances connect the same two nodes, hence share the same temperature difference \( \Delta T \), but their heat flows add:

\[ \dot{Q} = \dot{Q}_1 + \dot{Q}_2 = \frac{\Delta T}{R_1} + \frac{\Delta T}{R_2} = \Delta T\left(\frac{1}{R_1}+\frac{1}{R_2}\right). \]

Therefore:

\[ \frac{1}{R_{\mathrm{eq}}} = \frac{1}{R_1} + \frac{1}{R_2}. \]

It is often convenient to define \( G = 1/R \) as a thermal conductance, making parallel addition linear: \( G_{\mathrm{eq}} = G_1 + G_2 \).

Many engineering systems require more than one lumped temperature state. For nodes \( i=1,\dots,n \) with temperatures \( T_i(t) \) and thermal capacitances \( C_i \), resistive couplings between nodes \( i \) and \( j \) are captured by \( R_{ij} \) (or conductance \( G_{ij}=1/R_{ij} \)). Let \( \dot{Q}_{u,i}(t) \) be net external heat injected into node \( i \).

For each node, conservation of energy gives:

\[ C_i \frac{dT_i}{dt} = \dot{Q}_{u,i}(t) + \sum_{j \in \mathcal{N}(i)} \frac{T_j(t) - T_i(t)}{R_{ij}}. \]

In matrix form, define the diagonal capacitance matrix \( \mathbf{C}=\mathrm{diag}(C_1,\dots,C_n) \), the temperature vector \( \mathbf{T}(t) \), and an input vector \( \mathbf{q}(t) \) collecting \( \dot{Q}_{u,i}(t) \) plus any affine terms from ambient nodes (treated as known signals). Then:

\[ \mathbf{C}\frac{d\mathbf{T}}{dt} = -\mathbf{G}\mathbf{T} + \mathbf{q}(t), \]

where \( \mathbf{G} \) is a symmetric “conductance Laplacian”:

\[ G_{ij} = \begin{cases} -\frac{1}{R_{ij}} & \text{if } i \neq j \text{ and nodes } i,j \text{ are connected} \\ \sum_{k \in \mathcal{N}(i)} \frac{1}{R_{ik}} & \text{if } i=j \\ 0 & \text{otherwise.} \end{cases} \]

flowchart TD
  Tinf["Ambient node: T_inf(t)"] --- Rth["Thermal resistance R_th"] --- T1["Lumped node: T(t)"]
  Qdot["Heat input Qdot_u(t)"] --> T1
  T1 --- Cth["Thermal capacitance C_th"]
  note1["Dynamics: C_th * dT/dt = Qdot_u - (T - T_inf)/R_th"]
        

Stability/passivity insight (no new control theory required):

With zero external heat input and constant ambient references, these networks dissipate stored energy through resistances. This can be formalized by an energy function.

Consider an autonomous thermal network with fixed ambient references absorbed into a shifted temperature vector. For clarity, assume no external heat injection: \( \mathbf{q}(t)=\mathbf{0} \). Then

\[ \mathbf{C}\frac{d\mathbf{T}}{dt} = -\mathbf{G}\mathbf{T}. \]

Define an energy-like storage function

\[ V(\mathbf{T}) = \frac{1}{2}\mathbf{T}^\top \mathbf{C}\mathbf{T}. \]

Since \( \mathbf{C} \) is diagonal with \( C_i \gt 0 \), we have \( V(\mathbf{T}) \ge 0 \) and \( V(\mathbf{T})=0 \) only at \( \mathbf{T}=\mathbf{0} \). Differentiate along trajectories:

\[ \dot{V}(\mathbf{T}) = \frac{1}{2}\left(\dot{\mathbf{T}}^\top \mathbf{C}\mathbf{T} + \mathbf{T}^\top \mathbf{C}\dot{\mathbf{T}}\right) = \mathbf{T}^\top \mathbf{C}\dot{\mathbf{T}} = -\mathbf{T}^\top \mathbf{G}\mathbf{T}. \]

For any resistive network, \( \mathbf{G} \) is symmetric positive semidefinite, implying \( \mathbf{T}^\top \mathbf{G}\mathbf{T} \ge 0 \) and thus \( \dot{V}(\mathbf{T}) \le 0 \).

Why is \( \mathbf{G} \) positive semidefinite?

Using the Laplacian structure, one can show

\[ \mathbf{T}^\top \mathbf{G}\mathbf{T} = \frac{1}{2}\sum_{(i,j)\in \mathcal{E}} \frac{(T_i - T_j)^2}{R_{ij}} \;\;\ge\; 0, \]

where \( \mathcal{E} \) is the set of edges (resistive connections). This identity is obtained by expanding \( \mathbf{T}^\top \mathbf{G}\mathbf{T} \) and grouping terms by edges. Therefore, stored energy never increases without injection, which matches the physical expectation that resistances dissipate heat flow.

We implement the single-node model \( C_{\mathrm{th}}\frac{dT}{dt} = \dot{Q}_u(t) - \frac{T - T_{\infty}}{R_{\mathrm{th}}} \) with a step heat input \( \dot{Q}_u(t)=\dot{Q}_0 \) and constant ambient \( T_{\infty} \).

8.1 Python (SciPy + optional symbolic check)

import numpy as np
from scipy.integrate import solve_ivp

# Parameters
C_th = 120.0      # J/K
R_th = 2.0        # K/W
T_inf = 25.0      # degC
Q0 = 10.0         # W

tau = R_th * C_th

def dTdt(t, T):
    Q_u = Q0  # step input
    return (Q_u - (T - T_inf)/R_th) / C_th

t_span = (0.0, 2000.0)
t_eval = np.linspace(t_span[0], t_span[1], 400)
sol = solve_ivp(dTdt, t_span, y0=[25.0], t_eval=t_eval)

T = sol.y[0]
T_ss = T_inf + R_th * Q0
print("tau =", tau, "s")
print("Steady-state T_ss =", T_ss, "degC")
print("T(t_end) =", T[-1], "degC")

Library notes (system-dynamics oriented): scipy.integrate provides robust time-domain simulation; sympy can be used to derive and verify closed-form solutions for linear ODEs; and later (Chapter 6), transfer-function representations can be built from this ODE if desired.

8.2 C++ (Boost.ODEInt)

#include <iostream>
#include <vector>
#include <cmath>
#include <boost/numeric/odeint.hpp>

using state_type = double;
using namespace boost::numeric::odeint;

struct ThermalODE {
    double C_th, R_th, T_inf, Q0;
    ThermalODE(double C, double R, double Tinf, double Q) : C_th(C), R_th(R), T_inf(Tinf), Q0(Q) {}

    void operator()(const state_type &T, state_type &dTdt, const double /*t*/) const {
        double Q_u = Q0; // step input
        dTdt = (Q_u - (T - T_inf)/R_th) / C_th;
    }
};

int main() {
    double C_th = 120.0, R_th = 2.0, T_inf = 25.0, Q0 = 10.0;
    ThermalODE sys(C_th, R_th, T_inf, Q0);

    state_type T = 25.0;
    double t0 = 0.0, t1 = 2000.0, dt = 1.0;

    runge_kutta4<state_type> stepper;
    for(double t = t0; t <= t1; t += dt) {
        if ((static_cast<int>(t) % 200) == 0) {
            std::cout << t << " " << T << std::endl;
        }
        stepper.do_step(sys, T, t, dt);
    }
    return 0;
}

Library notes: Boost.ODEInt provides production-quality ODE integration commonly used for dynamic system simulation in C++. For multi-node models, \( \mathbf{T} \) becomes a vector state and the same integrators apply.

8.3 Java (Apache Commons Math ODE Integrators)

import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
import org.apache.commons.math3.ode.sampling.StepHandler;
import org.apache.commons.math3.ode.sampling.StepInterpolator;

public class ThermalSingleNode {
    static class ThermalODE implements FirstOrderDifferentialEquations {
        final double Cth, Rth, Tinf, Q0;
        ThermalODE(double Cth, double Rth, double Tinf, double Q0) {
            this.Cth = Cth; this.Rth = Rth; this.Tinf = Tinf; this.Q0 = Q0;
        }
        public int getDimension() { return 1; }
        public void computeDerivatives(double t, double[] y, double[] yDot) {
            double T = y[0];
            double Qu = Q0; // step input
            yDot[0] = (Qu - (T - Tinf)/Rth) / Cth;
        }
    }

    public static void main(String[] args) {
        double Cth = 120.0, Rth = 2.0, Tinf = 25.0, Q0 = 10.0;
        ThermalODE ode = new ThermalODE(Cth, Rth, Tinf, Q0);

        double t0 = 0.0, t1 = 2000.0;
        double[] y0 = new double[] {25.0};

        DormandPrince853Integrator integrator =
            new DormandPrince853Integrator(1.0e-6, 10.0, 1.0e-9, 1.0e-9);

        integrator.addStepHandler(new StepHandler() {
            public void init(double t0, double[] y0, double t) {}
            public void handleStep(StepInterpolator interpolator, boolean isLast) {
                double t = interpolator.getCurrentTime();
                if (((int)t) % 200 == 0) {
                    double[] y = interpolator.getInterpolatedState();
                    System.out.println(t + " " + y[0]);
                }
            }
        });

        integrator.integrate(ode, t0, y0, t1, y0);
    }
}

Library notes: Apache Commons Math offers standard ODE solvers suitable for system dynamics simulations in Java. For networks, you implement a higher-dimensional ODE with a state vector \( \mathbf{T}(t) \).

8.4 MATLAB + Simulink (ODE simulation and block implementation)

% Parameters
C_th = 120;   % J/K
R_th = 2;     % K/W
T_inf = 25;   % degC
Q0 = 10;      % W
tau = R_th*C_th;

% ODE: dT/dt = (Q0 - (T - T_inf)/R_th)/C_th
f = @(t,T) (Q0 - (T - T_inf)/R_th)/C_th;

[t,T] = ode45(f, [0 2000], 25);

T_ss = T_inf + R_th*Q0;
disp(['tau = ', num2str(tau), ' s'])
disp(['T_ss = ', num2str(T_ss), ' degC'])

% --- Simulink mapping (conceptual) ---
% Implement: dT/dt = (1/C_th)*Q0 - (1/(R_th*C_th))*T + (1/(R_th*C_th))*T_inf
% Blocks:
% 1) Constant Q0, Gain 1/C_th
% 2) Integrator -> T
% 3) Feedback path: T -> Gain (1/(R_th*C_th)) -> subtract from input
% 4) Constant T_inf -> Gain (1/(R_th*C_th)) -> add to input

Toolbox notes: For higher-fidelity thermal component modeling, MATLAB/Simulink ecosystems provide Simscape Thermal blocks; however, the lumped ODE model above is the foundational representation that remains essential for analysis and control-oriented modeling.

8.5 Wolfram Mathematica (symbolic + numeric solution)

Cth = 120.0;
Rth = 2.0;
Tinf = 25.0;
Q0 = 10.0;

(* ODE: Cth T'(t) == Q0 - (T(t)-Tinf)/Rth *)
ode = Cth*T'[t] == Q0 - (T[t] - Tinf)/Rth;
ic = T[0] == 25.0;

solSym = DSolve[{ode, ic}, T, t][[1]];
Tsol[t_] = T[t] /. solSym;

tau = Rth*Cth;
Tss = Tinf + Rth*Q0;

{tau, Tss, Simplify[Tsol[t]]}

(* Numeric evaluation *)
Plot[Evaluate[Tsol[t]], {t, 0, 2000}, PlotRange -> All]

Problem 1 (Lumped validity and time constant): A solid sphere of radius \( r \) has conductivity \( k \) and is cooled by convection with coefficient \( h \) over its entire surface. (a) Compute \( L_c \) and \( \mathrm{Bi} \). (b) State the lumped validity condition. (c) If lumped is valid, write the time constant in terms of \( m,c_p,h,A_s \).

Solution:

(a) For a sphere, \( V=\frac{4}{3}\pi r^3 \) and \( A_s=4\pi r^2 \), so \( L_c=\frac{V}{A_s}=\frac{r}{3} \). Hence

\[ \mathrm{Bi} = \frac{h L_c}{k} = \frac{h r}{3k}. \]

(b) Lumped validity (common criterion): \( \mathrm{Bi} \lt 0.1 \). (c) With convection only, \( R_{\mathrm{th}}=\frac{1}{hA_s} \) and \( C_{\mathrm{th}}=mc_p \), thus \( \tau = R_{\mathrm{th}}C_{\mathrm{th}} = \frac{mc_p}{hA_s} \).

Problem 2 (Series conduction + convection): A wall conducts heat through a slab of thickness \( L \), area \( A \), and conductivity \( k \), and then convects to air with coefficient \( h \) over area \( A \). Find the equivalent thermal resistance between the left wall temperature \( T_1 \) and ambient \( T_{\infty} \).

Solution:

The conduction resistance is \( R_{\mathrm{cond}}=\frac{L}{kA} \) and the convection resistance is \( R_{\mathrm{conv}}=\frac{1}{hA} \). They are in series, so

\[ R_{\mathrm{eq}} = \frac{L}{kA} + \frac{1}{hA}. \]

Problem 3 (Step heating and 63.2% rule): A lumped body with \( R_{\mathrm{th}} \), \( C_{\mathrm{th}} \) is initially at ambient \( T_{\infty} \). At \( t=0 \), a constant heater input \( \dot{Q}_0 \) is applied. (a) Derive \( T(t) \). (b) Show that at \( t=\tau \) the temperature has reached \( 1-e^{-1} \) of its final rise.

Solution:

(a) From Section 3 with \( T(0)=T_{\infty} \), steady-state is \( T_{\mathrm{ss}}=T_{\infty}+R_{\mathrm{th}}\dot{Q}_0 \). Thus

\[ T(t) = T_{\infty}+R_{\mathrm{th}}\dot{Q}_0\left(1-e^{-t/(R_{\mathrm{th}}C_{\mathrm{th}})}\right). \]

(b) Evaluate at \( t=\tau=R_{\mathrm{th}}C_{\mathrm{th}} \):

\[ T(\tau) = T_{\infty}+R_{\mathrm{th}}\dot{Q}_0\left(1-e^{-1}\right). \]

The fraction of final rise is \( 1-e^{-1}\approx 0.632 \).

Problem 4 (Two-node lumped model): Two lumped masses with temperatures \( T_1(t) \), \( T_2(t) \) have capacitances \( C_1 \), \( C_2 \) and are connected by thermal resistance \( R_{12} \). Node 2 convects to ambient \( T_{\infty} \) through \( R_{2\infty} \). Node 1 receives heat input \( \dot{Q}_u(t) \). Write the coupled ODEs.

Solution:

Heat flow from 1 to 2 is \( \frac{T_1-T_2}{R_{12}} \). Heat flow from 2 to ambient is \( \frac{T_2-T_{\infty}}{R_{2\infty}} \). Energy balances:

\[ C_1\frac{dT_1}{dt} = \dot{Q}_u(t) - \frac{T_1-T_2}{R_{12}}, \qquad C_2\frac{dT_2}{dt} = \frac{T_1-T_2}{R_{12}} - \frac{T_2-T_{\infty}}{R_{2\infty}}. \]

Problem 5 (Parameter identification from a step test): A lumped thermal system is modeled as \( C_{\mathrm{th}}\frac{dT}{dt}=\dot{Q}_0-\frac{T-T_{\infty}}{R_{\mathrm{th}}} \) with step input \( \dot{Q}_0 \). A test yields steady-state rise \( \Delta T_{\mathrm{ss}} = T_{\mathrm{ss}}-T_{\infty} \) and the time \( t_{63} \) at which \( T(t) \) reaches 63.2% of its total rise. Express \( R_{\mathrm{th}} \) and \( C_{\mathrm{th}} \) in terms of measured quantities.

Solution:

From steady-state, \( \Delta T_{\mathrm{ss}} = R_{\mathrm{th}}\dot{Q}_0 \), hence

\[ R_{\mathrm{th}} = \frac{\Delta T_{\mathrm{ss}}}{\dot{Q}_0}. \]

The 63.2% time equals the time constant: \( t_{63}=\tau=R_{\mathrm{th}}C_{\mathrm{th}} \). Therefore

\[ C_{\mathrm{th}} = \frac{t_{63}}{R_{\mathrm{th}}} = \frac{t_{63}\dot{Q}_0}{\Delta T_{\mathrm{ss}}}. \]

We introduced thermal effort/flow variables, derived lumped thermal-capacitance ODE models from first principles, established lumped validity via the Biot number, and proved that passive thermal RC networks dissipate stored energy using an explicit energy function argument. We also implemented canonical simulations in Python, C++, Java, MATLAB/Simulink, and Mathematica. These models directly parallel the RLC and hydraulic network modeling patterns from earlier lessons in this chapter, preparing you for multi-domain couplings in Lesson 5.