Heisler’s transient conduction charts (1947)

Separation of variables gives the cooling of a plane wall as a Fourier series whose eigenvalues solve \(\zeta\tan\zeta = Bi\). Before computers, nobody summed it by hand. In 1947 M. P. Heisler, working on induction heating, noticed that after a short time only the first term matters, and plotted that term as charts of centre temperature against the Fourier number for each Biot number. Generations of students read them with a ruler. plane_wall_transient() sums the series (or its first term, terms=1); this example redraws the centre chart, measures the one-term error, and times the quench of a steel plate.

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import brentq

from engineeringkit.heat_transfer import biot_number, fourier_number, plane_wall_transient

The centre-temperature chart

Heisler’s chart is labelled by \(1/Bi\); the curves are straight on a log scale once the first term dominates.

Fo = np.linspace(0, 4, 201)
fig, ax = plt.subplots(figsize=(6, 4))
for inv_bi in [0.0001, 0.5, 1, 2, 5, 10]:
    theta = plane_wall_transient(Fo, 1 / inv_bi, 0.0, terms=200)
    ax.semilogy(Fo, theta, label=f"1/Bi = {inv_bi:g}")
ax.set(xlabel="Fourier number", ylabel="centre excess temperature", ylim=(1e-3, 1.1), title="Heisler's plane-wall centre chart")
ax.legend(ncol=2)
fig.tight_layout()

How good is one term?

At the centre the one-term approximation stays within 2% of the full series for \(Fo > 0.2\), whatever the Biot number. Very early, many terms are needed: at \(Fo = 10^{-3}\) the centre has not yet felt the surface.

worst = 0.0
Fo_late = np.linspace(0.2, 3, 60)
for Bi in [0.01, 0.1, 1, 10, 100]:
    full = plane_wall_transient(Fo_late, Bi, 0.0)
    worst = max(worst, np.max(np.abs(plane_wall_transient(Fo_late, Bi, 0.0, terms=1) / full - 1)))
print(f"largest one-term error at the centre for Fo > 0.2: {worst:.2%}")
assert worst < 0.02
early = plane_wall_transient(1e-3, 10.0, 0.0, terms=200)
print(f"centre at Fo = 0.001 (200 terms): {early:.6f}")
assert np.isclose(early, 1.0, atol=1e-6)

largest one-term error at the centre for Fo > 0.2: 1.53%
centre at Fo = 0.001 (200 terms): 1.000000

Quenching a steel plate

A 50 mm steel plate (\(k = 40\) W/(m K), \(\alpha = 1.1\times10^{-5}\) m²/s) at 1100 K is quenched in oil at 330 K with \(h = 1000\) W/(m² K). When does its centre reach 600 K?

k, alpha, L, h, Ti, Tinf, target = 40.0, 1.1e-5, 0.025, 1000.0, 1100.0, 330.0, 600.0
Bi = biot_number(h, L, k)
theta_target = (target - Tinf) / (Ti - Tinf)
Fo_target = brentq(lambda F: plane_wall_transient(F, Bi, 0.0) - theta_target, 0.01, 50)
t_target = Fo_target * L**2 / alpha
print(f"Bi = {Bi:.3f}, Fo = {Fo_target:.3f}: the centre reaches 600 K after {t_target:.0f} s")
assert np.isclose(fourier_number(alpha, t_target, L), Fo_target)
surface = Tinf + (Ti - Tinf) * plane_wall_transient(Fo_target, Bi, 1.0)
print(f"the surface is then at {surface:.0f} K")

plt.show()

Bi = 0.625, Fo = 2.196: the centre reaches 600 K after 125 s
the surface is then at 534 K