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), and cylinder_transient() and sphere_transient() do the same for his cylinder and sphere charts. This example redraws the centre chart, measures the one-term error, times the quench of a steel plate, and compares the three shapes.

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

from engineeringkit.heat_transfer import biot_number, cylinder_transient, fourier_number, plane_wall_transient, sphere_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")

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

Slab, cylinder and sphere

A sphere has the most surface per unit volume, a slab the least, so at the same Biot and Fourier numbers (on the half-thickness or radius) the sphere’s centre cools first. For small Biot numbers all three become lumped bodies, \(\theta = \exp[-(V/A)\,hAt/(\rho c V)]\), which on the radius is \(e^{-Bi\,Fo}\), \(e^{-2Bi\,Fo}\) and \(e^{-3Bi\,Fo}\). The example steel ball is a 50 mm bearing ball in the same oil bath.

Bi_small, Fo_long = 0.01, 50.0
for transient, factor, name in ((plane_wall_transient, 1, "slab"), (cylinder_transient, 2, "cylinder"), (sphere_transient, 3, "sphere")):
    exact = transient(Fo_long, Bi_small)
    print(f"{name:8s}: centre {exact:.4f}, lumped {np.exp(-factor * Bi_small * Fo_long):.4f}")
    assert abs(exact / np.exp(-factor * Bi_small * Fo_long) - 1) < 0.01

Bi_ball = biot_number(h, L, k)
Fo_ball = brentq(lambda F: sphere_transient(F, Bi_ball) - theta_target, 0.01, 50)
print(f"a 50 mm ball reaches 600 K at its centre after {Fo_ball * L**2 / alpha:.0f} s, against {t_target:.0f} s for the plate")
assert Fo_ball < Fo_target

fig, ax = plt.subplots(figsize=(6, 4))
for transient, name in ((plane_wall_transient, "slab"), (cylinder_transient, "long cylinder"), (sphere_transient, "sphere")):
    ax.semilogy(Fo, transient(Fo, 1.0, 0.0, terms=200), label=name)
ax.set(xlabel="Fourier number", ylabel="centre excess temperature", ylim=(1e-3, 1.1), title="Bi = 1: the three Heisler shapes")
ax.legend()
fig.tight_layout()

plt.show()

slab    : centre 0.6085, lumped 0.6065
cylinder: centre 0.3697, lumped 0.3679
sphere  : centre 0.2245, lumped 0.2231
a 50 mm ball reaches 600 K at its centre after 42 s, against 125 s for the plate