Simulating Thermal Systems

Thermal simulation is one of the most accessible and practically valuable applications of numerical ODE solving. Nearly every engineering system — an electronic component cooling down after shutdown, a metal billet being quenched in a bath, a room losing heat through its walls overnight — can be described by a differential equation relating stored thermal energy to the rate at which heat flows in or out. When the geometry and material are simple enough that we can treat the entire body as having a single, uniform temperature at any instant, we arrive at what is called a lumped parameter thermal model. This topic develops that model thoroughly, derives and simulates Newton's Law of Cooling, extends the approach to networks of thermally connected nodes, and closes the loop by comparing simulation output against real measured data.

Lumped Parameter Thermal Models

The central assumption of the lumped model is spatial uniformity: temperature is treated as a function of time only, not of position within the body. This is a significant simplification. In reality, if you drop a cold steel ball into hot oil, the surface heats up before the core does. But if the ball is small enough and conducts heat internally much faster than heat flows across its surface, then the temperature gradient inside the ball is negligible compared with the temperature difference between the ball and the fluid. In that regime, a single number — the ball's temperature T(t) — describes its thermal state entirely.

The criterion that justifies this assumption is the Biot number:

Bi = h · L_c / k

where h is the convective heat transfer coefficient (W/m²·K), L_c is a characteristic length (volume divided by surface area for a convex body), and k is the thermal conductivity of the solid (W/m·K). When Bi < 0.1, the internal conduction resistance is at least ten times smaller than the external convection resistance, so the temperature inside the body is nearly uniform at every instant and the lumped model is valid. A copper sphere in slow air flow easily satisfies this (copper has very high k); a large block of insulating foam in the same flow would not.

The thermal analogy to an electrical RC circuit is exact and deeply useful. Define:

• Thermal capacitance: C = m · cp (J/K), where m is mass (kg) and cp is specific heat at constant pressure (J/kg·K). This quantity represents the body's ability to store thermal energy — just as electrical capacitance stores charge.
• Thermal resistance: R = 1 / (h · A) (K/W), where A is the surface area (m²). This quantity opposes the flow of heat just as electrical resistance opposes the flow of current. A large surface area or a high heat transfer coefficient both lower the resistance.

The product R·C = m·cp / (h·A) has units of seconds and is the system's time constant τ. This single parameter controls the speed of the thermal response, exactly as τ = RC controls the speed of an electrical transient.

Newton's Law of Cooling

With the lumped assumption in place, an energy balance on the body is straightforward. The rate of change of stored thermal energy equals the net rate of heat transfer to the body. Convection to an ambient fluid at temperature T∞ removes heat at a rate h·A·(T − T∞), so:

m · cp · dT/dt = -h · A · (T - T∞)

This is Newton's Law of Cooling expressed as a first-order linear ODE. The negative sign indicates that if the body is hotter than the ambient (T > T∞), heat flows out and the body cools; if the body is colder, heat flows in and the body warms. Dividing both sides by m·cp isolates the derivative:

dT/dt = -(h · A) / (m · cp) · (T - T∞)
      = -(1/τ) · (T - T∞)

where τ = m·cp / (h·A). This equation is separable. Introducing the excess temperature θ = T − T∞ (which is constant since T∞ is constant), we get dθ/dt = −θ/τ, which integrates immediately to:

T(t) = T∞ + (T₀ - T∞) · exp(-t / τ)

Here T₀ is the body's initial temperature. This is an exponential decay (or rise) toward T∞. After one time constant the body has covered 63.2% of the gap to T∞; after five time constants it is within 0.7% of T∞ and is effectively at equilibrium. A small τ (small mass, small specific heat, large surface area, or vigorous convection) means fast equilibration; a large τ means sluggish response.

As a concrete example, consider a 0.5 kg aluminum block (cp = 900 J/kg·K) with a surface area of 0.05 m² cooled by natural convection with h = 10 W/m²·K, starting at T₀ = 150 °C in a room at T∞ = 25 °C:

τ = (0.5 × 900) / (10 × 0.05) = 450 / 0.5 = 900 s  (15 minutes)

After 15 minutes the block is at 25 + (150−25)·e⁻¹ ≈ 71 °C. After an hour (≈ 4τ) it is within about 2 °C of room temperature.

Simulating Cooling and Heating Transients

Even though this ODE has an exact analytical solution, implementing it numerically is the foundation for everything that follows — because more complex systems (multiple nodes, nonlinear radiation, time-varying boundary conditions) do not have closed-form solutions. The workflow has four steps.

First, define all physical parameters as named Python constants. This makes the code readable and makes parameter sweeps trivial:

import numpy as np
import matplotlib.pyplot as plt

# Physical parameters
m   = 0.5       # mass, kg
cp  = 900.0     # specific heat, J/(kg·K)
h   = 10.0      # heat transfer coefficient, W/(m²·K)
A   = 0.05      # surface area, m²
T_inf = 25.0    # ambient temperature, °C
T0    = 150.0   # initial temperature, °C

tau = m * cp / (h * A)
print(f"Time constant τ = {tau:.1f} s")

Second, implement the derivative function. Its signature must accept the current time and current state (temperature) and return dT/dt:

def dTdt(t, T):
    return -(h * A) / (m * cp) * (T - T_inf)

Third, run the numerical integration loop. Here we use the classical fourth-order Runge-Kutta (RK4) method, which is far more accurate than Euler for the same step size:

dt   = 10.0          # time step, s
t_end = 5 * tau      # integrate for 5 time constants
t_arr = np.arange(0, t_end + dt, dt)
T_arr = np.zeros(len(t_arr))
T_arr[0] = T0

for i in range(len(t_arr) - 1):
    t = t_arr[i]
    T = T_arr[i]
    k1 = dTdt(t,        T)
    k2 = dTdt(t + dt/2, T + dt/2 * k1)
    k3 = dTdt(t + dt/2, T + dt/2 * k2)
    k4 = dTdt(t + dt,   T + dt   * k3)
    T_arr[i+1] = T + (dt/6) * (k1 + 2*k2 + 2*k3 + k4)

Fourth, compute the analytical solution at every time point and overlay it on the plot to validate accuracy:

T_analytical = T_inf + (T0 - T_inf) * np.exp(-t_arr / tau)

plt.figure()
plt.plot(t_arr / 60, T_arr,        label='RK4 simulation')
plt.plot(t_arr / 60, T_analytical, '--', label='Analytical')
plt.xlabel('Time (min)')
plt.ylabel('Temperature (°C)')
plt.legend()
plt.title('Lumped-Parameter Cooling Transient')
plt.tight_layout()
plt.show()

With a step size of 10 s and τ = 900 s, RK4 produces results visually indistinguishable from the analytical curve. The maximum absolute error across 5τ will typically be well under 0.01 °C — a reassuring validation before moving to more complex models.

Multi-Node Thermal Networks

Many real systems cannot be treated as a single lump. A printed circuit board has a processor die, a heat spreader, a thermal interface material, a heat sink base, and fins, each with its own thermal mass and connected to neighbors through thermal resistances. The appropriate model is a network of N nodes, each with its own capacitance Cᵢ and exchanging heat with adjacent nodes through conductances Gᵢⱼ = 1/Rᵢⱼ.

The energy balance on node i is:

Cᵢ · dTᵢ/dt = Σⱼ Gᵢⱼ · (Tⱼ - Tᵢ) + Qᵢ

where the sum runs over all nodes j connected to i, and Qᵢ is any external heat source at node i (e.g., power dissipation). Collecting all nodes into a state vector T = [T₁, T₂, …, Tₙ]ᵀ, this system of ODEs takes the matrix form:

C · dT/dt = A · T + b

where C is the diagonal capacitance matrix (Cᵢᵢ = Cᵢ), A is the conductance matrix (Aᵢⱼ = Gᵢⱼ for i ≠ j; Aᵢᵢ = −Σⱼ Gᵢⱼ), and b contains boundary conditions and heat sources. To put it in standard form for numerical integration, multiply both sides by C⁻¹:

dT/dt = C⁻¹ · (A · T + b)

Consider a simple three-node model: a heat source (node 1) attached to a spreader (node 2) attached to ambient through a convective resistance (node 3 = ambient, fixed). With NumPy this becomes:

import numpy as np

# Capacitances (J/K)
C1, C2 = 0.5, 2.0
C_mat = np.diag([C1, C2])

# Conductances (W/K)
G12 = 2.0   # between node 1 and node 2
G2a = 1.0   # between node 2 and ambient (node fixed at T_inf)
T_inf = 25.0
Q1   = 5.0  # heat source at node 1, W

# Conductance matrix A (2×2, node 3 = ambient absorbed into b)
A_mat = np.array([[-G12,  G12],
                  [ G12, -(G12 + G2a)]])

# Forcing vector b
b_vec = np.array([Q1, G2a * T_inf])

C_inv = np.linalg.inv(C_mat)

def dTdt_multi(t, T_vec):
    return C_inv @ (A_mat @ T_vec + b_vec)

# Initial conditions and RK4 loop
T0_vec = np.array([25.0, 25.0])
dt = 1.0
t_end = 300.0
t_arr = np.arange(0, t_end + dt, dt)
T_history = np.zeros((len(t_arr), 2))
T_history[0] = T0_vec

for i in range(len(t_arr) - 1):
    t  = t_arr[i]
    Tv = T_history[i]
    k1 = dTdt_multi(t,        Tv)
    k2 = dTdt_multi(t + dt/2, Tv + dt/2 * k1)
    k3 = dTdt_multi(t + dt/2, Tv + dt/2 * k2)
    k4 = dTdt_multi(t + dt,   Tv + dt   * k3)
    T_history[i+1] = Tv + (dt/6) * (k1 + 2*k2 + 2*k3 + k4)

The derivative function now returns a NumPy array of length 2 — one derivative per node — and the matrix-vector multiply @ operator handles the coupling cleanly. Steady state can be verified analytically: set dT/dt = 0, so A·T* + b = 0, giving T* = −A⁻¹·b. For the parameters above this yields T₁* ≈ 32.5 °C and T₂* ≈ 30.0 °C, values that the simulation should approach asymptotically.

As the network grows to dozens or hundreds of nodes (a full finite-element thermal model of a PCB, for instance), building A by hand becomes impractical. The standard approach is to assemble the conductance matrix programmatically from a list of node connections, then use scipy.integrate.solve_ivp with a stiff solver such as method='Radau' or method='BDF', since multi-node thermal networks with very different time constants can be stiff systems where explicit methods like RK4 require impractically small time steps.

Comparison with Analytical and Experimental Data

A simulation that has never been compared to ground truth is of limited value. Validation — checking that the model reproduces known behavior — is an essential engineering step.

For a single-node system the analytical solution is known exactly. Plotting it alongside the numerical result is the first check. To go further, compute the root-mean-square error (RMSE) between the two:

rmse = np.sqrt(np.mean((T_arr - T_analytical)**2))
print(f"RMSE between RK4 and analytical: {rmse:.4f} °C")

A small RMSE (typically much less than 0.1 °C for reasonable step sizes) confirms that the numerical scheme is implemented correctly and the step size is adequate.

When experimental data exist — say, a thermocouple reading from a cooling test — the comparison is more informative. Load the data, align the time axes, interpolate if the sampling rate differs from the simulation's, then compute RMSE between simulation and measurement. A large discrepancy points to an incorrect parameter (wrong h, wrong cp, unexpected radiation losses) or a violated assumption (Bi ≥ 0.1, multi-dimensional conduction).

Parameter sensitivity analysis answers the question: which parameter matters most? Perturb each parameter by ±10% in isolation and observe how much the temperature response changes. For the single-node model, τ = m·cp/(h·A), so temperature at any fixed time depends on all four parameters, but because h and A always appear together as the product h·A, they are indistinguishable from cooling data alone — an important identifiability insight.

Curve fitting using scipy.optimize.curve_fit can extract thermal properties directly from measured data. If you trust the analytical form, fit T∞, T₀, and τ to the measured time series:

from scipy.optimize import curve_fit

def cooling_model(t, T_inf, T0, tau):
    return T_inf + (T0 - T_inf) * np.exp(-t / tau)

# t_meas and T_meas are 1-D arrays from the experiment
popt, pcov = curve_fit(cooling_model, t_meas, T_meas,
                       p0=[25.0, 150.0, 900.0])
T_inf_fit, T0_fit, tau_fit = popt
print(f"Fitted τ = {tau_fit:.1f} s")

Once τ is known, and if you independently measure mass and surface area, you can back-calculate h — a non-trivial experimental heat transfer result obtained from a simple cooling test.

Connecting Data Analysis to Thermal Simulation

Real experiments produce CSV files, not perfect analytical signals. The full workflow of loading, cleaning, and using measured data as a simulation input ties data analysis skills directly to physical modeling.

Load and inspect the data with pandas:

import pandas as pd

df = pd.read_csv('cooling_test.csv')
print(df.head())
print(df.dtypes)
df.dropna(inplace=True)
df['time_s'] = pd.to_numeric(df['time_s'], errors='coerce')
df['T_sample'] = pd.to_numeric(df['T_sample'], errors='coerce')
df['T_ambient'] = pd.to_numeric(df['T_ambient'], errors='coerce')
df.dropna(inplace=True)

A critical subtlety: the ambient temperature T∞ is often not perfectly constant during an experiment — the room heats up, or the bath temperature drifts. Rather than assuming T∞ is fixed, use the measured ambient as a time-varying boundary condition. The ODE becomes:

dT/dt = -(h*A)/(m*cp) * (T - T_inf(t))

where T∞(t) is now a function of time. Implement this with NumPy interpolation so the derivative function can look up the ambient temperature at any instant during integration:

from scipy.interpolate import interp1d

T_inf_func = interp1d(df['time_s'].values, df['T_ambient'].values,
                      kind='linear', fill_value='extrapolate')

def dTdt_varying(t, T):
    return -(h * A) / (m * cp) * (T - T_inf_func(t))

Pass this function to scipy.integrate.solve_ivp for clean, robust integration:

from scipy.integrate import solve_ivp

t_span = (df['time_s'].iloc[0], df['time_s'].iloc[-1])
t_eval = df['time_s'].values
sol = solve_ivp(dTdt_varying, t_span, [T0], t_eval=t_eval,
                method='RK45', rtol=1e-6, atol=1e-8)
T_sim = sol.y[0]

Now plot simulation and measurement on the same axes:

plt.figure()
plt.plot(df['time_s'] / 60, df['T_sample'], 'o', ms=3, label='Measured')
plt.plot(df['time_s'] / 60, T_sim,          '-',       label='Simulation')
plt.xlabel('Time (min)')
plt.ylabel('Temperature (°C)')
plt.legend()
plt.tight_layout()
plt.show()

If the curves diverge significantly, the next step is parameter iteration — adjusting h (or equivalently τ) to minimize RMSE. This can be automated with curve_fit or scipy.optimize.minimize, using the simulation output (not the analytical formula) as the model. This closes the full engineering loop: physical model → numerical simulation → comparison with data → parameter refinement → improved model → re-simulate. Each iteration produces a model that better represents reality, and the RMSE provides an objective, quantitative measure of progress.