A cooling curve, integrated step by step

A part comes out of a furnace and cools in still air. If it is small enough and conductive enough that its inside and its surface stay at the same temperature, the whole thing is one lumped heat capacity losing heat to the room, and its temperature obeys a first-order differential equation with a known solution. That makes it the right problem to check a numerical integrator against: the answer is already known, so the question is whether the method reproduces it.

use steel

use constants

digits 4

The part, and whether the method applies

L_{side}

50 mm

V

{L_{side}}^{3}

{(50 mm)}^{3}

125000 {mm}^{3}

A_{surf}

6 \cdot {L_{side}}^{2}

6 \cdot {(50 mm)}^{2}

15000 {mm}^{2}

m_{part}

ρ_{steel} \cdot V

7850 m⁻³·kg \cdot 125000 {mm}^{3}

0.9813 kg

c_{steel}

\frac{480 J}{kg \cdot K}

h_{air}

\frac{25 W}{m^{2} \cdot K}

k_{steel}

\frac{45 W}{m \cdot K}

The lumped assumption holds when conduction inside the part is fast compared with convection off its surface — the Biot number, on the characteristic length that is volume over area. Below about a tenth, the inside and the surface differ by less than the other approximations here.

L_{char}

\frac{V}{A_{surf}}

\frac{125000 {mm}^{3}}{15000 {mm}^{2}}

8.333 mm

Bi

\frac{h_{air} \cdot L_{char}}{k_{steel}}

\frac{25 kg·s⁻³·K⁻¹ \cdot 8.333 mm}{45 m·kg·s⁻³·K⁻¹}

0.00463

check

Bi \leq 0.1

0.00463 \leq 0.1

pass

The equation

The time constant is the heat the part stores divided by the rate it loses it at one kelvin of excess.

τ

\frac{m_{part} \cdot c_{steel}}{h_{air} \cdot A_{surf}}

\frac{0.9813 kg \cdot 480 m²·s⁻²·K⁻¹}{25 kg·s⁻³·K⁻¹ \cdot 15000 {mm}^{2}}

1256 s

T_{room}

300 K

T_{start}

500 K

Absolute temperatures throughout. An offset scale has no place in this arithmetic — the equation subtracts one temperature from another and scales the difference, and °C is refused for it — so the worksheet works in kelvin and converts at the end if a reader wants degrees.

fn cooling defined

Integrating it

rk4 names its method and takes its step count, because both change the answer. Fifty steps over five time constants is a step of a tenth of a time constant, which for a smooth exponential is far finer than it needs to be.

t_{end}

5 \cdot τ

5 \cdot 1256 s

6280 s

steps

50

history

rk4 (cooling, T_{start}, 0 s, t_{end}, steps)

rk4 (cooling, 500 K, 0 s, 6280 s, 50)

=[51x2 table]

Against the answer that is already known

fn exact defined

T_{final_rk4}

{history}_{steps + 1, 2}

{([51x2 table])}_{50 + 1, 2}

301.3 K

T_{final_exact}

exact (\frac{t_{end}}{s})

exact (\frac{6280 s}{s})

301.3 K

{error}_{K}

| T_{final_rk4} - T_{final_exact} |

| 301.3 K - 301.3 K |

6.104 \times 10^{−6} K

check

{error}_{K} \leq 0.001 K

6.104 \times 10^{−6} K \leq 0.001 K

pass

Five hundred kelvin falling to three hundred and one, and the two methods agree to a thousandth of a kelvin. That is the fourth-order method earning its name: halving the step divides the error by sixteen.

The step count is part of the answer

The same integration over five steps rather than fifty — a step of a whole time constant, which is the sort of thing that looks reasonable and is not.

coarse

rk4 (cooling, T_{start}, 0 s, t_{end}, 5)

rk4 (cooling, 500 K, 0 s, 6280 s, 5)

=[[0 s, 500 K], [1256 s, 375 K], [2512 s, 328.1 K], [3768 s, 310.5 K], [5024 s, 304 K], [6280 s, 301.5 K]]

T_{final_coarse}

{coarse}_{6, 2}

{([[0 s, 500 K], [1256 s, 375 K], [2512 s, 328.1 K], [3768 s, 310.5 K], [5024 s, 304 K], [6280 s, 301.5 K]])}_{6, 2}

301.5 K

{error}_{coarse}

| T_{final_coarse} - T_{final_exact} |

| 301.5 K - 301.3 K |

0.1356 K

Still close, because an exponential decay is the easiest case there is. A worksheet integrating something stiffer would find the same step size wrong by a great deal more, and nothing here would say so — which is why the step count is written where a reader can see it rather than chosen by the engine.

The curve

axis y limits

plot (history)