How to Solve the Heat Conduction Equation

If a plate of thickness a is prepared at temperature T₀ and left alone, how does its temperature distribution evolve? The first half follows the lecture notes: we derive the heat conduction equation from Fourier’s law and solve it by separation of variables and a Fourier series. The second half goes beyond the notes and solves for the surroundings of the plate (air or water) as well. What appears there is a quantity different from the thermal diffusivity κ² — the thermal effusivity b = √(Kρc). The answer to “why does metal feel cold to the touch?” lies in this b. Finally, we replace the boundary condition with Newton’s law of cooling to include convection and radiation. The eigenvalue condition changes from ωa = nπ to the transcendental equation μ tan μ = Bi, and you can see that when the boundary condition changes, the whole eigenvalue problem changes with it — what happens around an electric furnace for firing ceramics comes out directly.

1. From Fourier’s Law

Place a rod of cross-sectional area S along the x axis. Let the temperatures at its two ends be T(0) and T(a), with T(0) > T(a). The amount of heat Q that flows through a cross-section during a time t should be proportional to t, to S, and to the temperature gradient. If we call the constant of proportionality K, then

Q = K t S \frac{T (0) - T (a)}{a - 0}

(1)

This K is the thermal conductivity. In differential form, the heat flux (the heat per unit time per unit area) is q = −K ∂T/∂x. The minus sign expresses that “heat flows from the higher temperature to the lower.”

2. From the Balance in a Small Interval to the Equation

Let A and B be the two cross-sections at x and x + dx. The heat entering through A and the heat leaving through B are

d Q_{A} = K S {(- \frac{\partial T}{\partial x})}_{A} d t, d Q_{B} = K S {(- \frac{\partial T}{\partial x})}_{B} d t

Expanding the gradient at B in a Taylor series around A gives (∂T/∂x)_B ≈ (∂T/∂x)_A + (∂²T/∂x²)_A dx. Therefore the heat that accumulates in this interval is

d Q = d Q_{A} - d Q_{B} = K S \frac{\partial^{2} T}{\partial x^{2}} d x d t

(2)

On the other hand, if that heat raised the temperature by dT, then in terms of the density ρ and the specific heat c we have

d Q = c ρ S d x d T

(3)

Equating Eqs. (2) and (3) and dividing by S, dx, and dt, we obtain the heat conduction equation.

\frac{\partial T}{\partial t} = κ^{2} \frac{\partial^{2} T}{\partial x^{2}}, κ^{2} \equiv \frac{K}{ρ c}

(4)

κ² is a quantity of “speed” The thermal diffusivity κ² = K/(ρc) is how readily heat is conducted, K, divided by how readily it is stored, ρc. The smaller the denominator, the larger the value — the less a material stores heat and the more it lets heat pass straight through, the faster it is. Its unit is m²/s, that is, length squared divided by time. From this we can read the quantity that appears later, √(κ²t), as “the depth to which heat has soaked in.”

Material ValuesK is from the table in the lecture notes; ρ and c are standard literature values. The last column is a guide to the time a 1 cm plate takes to even out, a²/κ²
Material	K [W/(m·K)]	ρ [kg/m³]	c [J/(kg·K)]	κ² [m²/s]	b [J/(m²·K·√s)]	1 cm
Silver	428	10490	235	1.74×10⁻⁴	32482	0.58 s
Copper	403	8960	385	1.17×10⁻⁴	37285	0.86 s
Aluminum	236	2700	897	9.74×10⁻⁵	23907	1.0 s
Iron	83.5	7874	449	2.36×10⁻⁵	17182	4.2 s
Air	0.0241	1.293	1005	1.86×10⁻⁵	5.6	5.4 s
Ice	2.2	917	2050	1.17×10⁻⁶	2034	85 s
Glass	0.75	2500	750	4.00×10⁻⁷	1186	250 s
Water	0.58	999.7	4192	1.38×10⁻⁷	1559	723 s

What to notice in this table Note that the κ² of air is about the same as that of iron. In terms of how fast heat “spreads” alone, air is on a par with metals. Why then does it not feel that way? The answer lies not in κ² but in b, the second column from the right (§7).

3. The Steady State First — Lecture Notes (2)

Suppose the temperatures of the two faces are held at T₀ and T_a, and that enough time has passed for things to settle down. Since ∂T/∂t = 0, Eq. (4) gives ∂²T/∂x² = 0; that is, T is a linear function of x.

T (x) = \frac{T_{a} - T_{0}}{a} x + T_{0}

(5)

A straight line. This is enough when we consider the insulation of a wall. Things get interesting once time enters the picture.

4. Separation of Variables — Lecture Notes (3)

We solve the case where both faces of the plate are held at T = 0 and the initial distribution is given. Putting T(x,t) = X(x)T(t), substituting into Eq. (4), and dividing by XT, we get

\frac{\dot{T}}{T} \frac{1}{κ^{2}} = \frac{X^{″}}{X} = α (constant)

(6)

The left side is a function of t only and the right side a function of x only. Since the two are equal, both must be a constant. What remains is to consider cases according to the sign of α.

[i] α > 0: X = Ae^{√α x} + Be^{−√α x}. Imposing X(0) = X(a) = 0, the determinant is e^−√αa − e^√αa ≠ 0, so A = B = 0. Only the trivial solution remains.
[ii] α = 0: X = Ax + B. With both ends at 0, this also gives X ≡ 0.
[iii] α < 0: putting α = −ω², X = A cos ωx + B sin ωx. From X(0) = 0 we get A = 0, and from X(a) = 0 we get sin ωa = 0, that is, ωa = nπ.

The key point The boundary condition has discretized the allowed values of ω. Only the “discrete values” ω_n = nπ/a are allowed. — The structure is exactly the same as that of a particle confined in a potential well, whose energy takes only discrete values. The heat conduction equation and the Schrödinger equation differ, as far as the form of the equation goes, only in whether an i is present. — Still, that single i decides whether the solution smooths out and dies away, or whether the phase rotates and produces interference. The similarity extends only as far as the Laplacian eigenvalue problem.

For the time part, the left-hand equality of Eq. (6) gives Ṫ/T = −n²π²κ²/a², that is,

T (t) = e^{- \frac{n^{2} π^{2} κ^{2} t}{a^{2}}}

(7)

5. Determining the Fourier Coefficients

We superpose the allowed solutions X_n = sin(nπx/a) and fit them to the initial condition.

T (x, 0) = \sum_{n = 1}^{\infty} B_{n} \sin \frac{n π x}{a}, B_{n} = \frac{2}{a} \int_{0}^{a} T (x, 0) \sin \frac{n π x}{a} d x

(8)

We use the initial condition of the notes, T(x,0) = (4T₀/a²)·x(a−x) (a parabola that equals T₀ at the center), and compute the coefficients. Only two integrals are needed.

\int_{0}^{a} x \sin \frac{n π x}{a} d x = \frac{a^{2}}{n π} {(- 1)}^{n + 1}

\int_{0}^{a} x^{2} \sin \frac{n π x}{a} d x = \frac{a^{3}}{n π} {(- 1)}^{n + 1} - \frac{2 a^{3}}{n^{3} π^{3}} [1 - {(- 1)}^{n}]

Forming a·(first) − (second), the (−1)ⁿ⁺¹ terms cancel completely, and only one term remains.

B_{n} = \frac{8 T_{0}}{a^{3}} \cdot \frac{2 a^{3}}{n^{3} π^{3}} [1 - {(- 1)}^{n}] = {\begin{cases} \frac{32 T_{0}}{n^{3} π^{3}} & (n : odd) \\ 0 & (n : even) \end{cases}

(9)

Three checks ① It is natural that the even modes vanish. The initial distribution is symmetric about x = a/2, but sin(2πx/a) is antisymmetric. Expanding a symmetric function in an antisymmetric one always gives a coefficient of 0. If an even term survives, the calculation is wrong — a handy marker.
② It is also natural that the coefficients fall off as n³. The smoother the initial distribution (the fewer corners it has), the faster the high components shrink. For the “uniform T₀,” which has corners, the coefficients fall off only as n¹, as we will see later.
③ Add up the terms at the center. At x = a/2, B₁ − B₃ + B₅ − … = (32T₀/π³)·β(3), and since β(3) = π³/32 this is exactly T₀.

The case where the whole plate is uniformly at T₀ can be computed with the same formula, and it is even simpler.

B_{n} = \frac{2 T_{0}}{n π} [1 - {(- 1)}^{n}] = \frac{4 T_{0}}{n π} (n : odd)

(10)

To summarize, the solution we are looking for is

T (x, t) = \sum_{n : odd} B_{n} e^{- n^{2} π^{2} κ^{2} t / a^{2}} \sin \frac{n π x}{a}

(11)

Fig. 1: How many terms of the series to add. This shows the distribution at t = 0 built up by adding terms in order from n = 1. The parabola is matched quickly, but the uniform T₀ has corners, so the ends ripple — the Gibbs phenomenon. However many terms you add, the overshoot remains; only its position moves closer to the edge. Its size is about 9 % of the whole step (2T₀ when the uniform T₀ is extended as an odd function), or about 18 % in terms of T₀ — if you switch the initial distribution to “Uniform T₀,” the “Maximum overshoot” below reads about 18 K. As time passes, the higher-order terms decay and the ripples disappear.

Number of terms: 5 (odd only) Initial distribution Time κ²t/a² = 0

6. Reading the Solution — n² and a²

Look at the exponent in Eq. (11). It contains n². That is, the n = 3 mode dies out 9 times faster than n = 1, and the n = 5 mode 25 times faster.

That is why heat conduction acts to “smooth things out” The finer the bumps (the higher the n), the faster they disappear, so however jagged the initial distribution, it soon becomes just the single hump of n = 1. After that it keeps that shape and shrinks exponentially. In the middle figure of the simulator (the bar graph of the modes), you can see the bars disappear in order from the right.

The slowest mode, n = 1, determines the time scale of the whole process. Its time constant is

τ = \frac{a^{2}}{π^{2} κ^{2}}

(12)

It is proportional to the square of a. Doubling the thickness makes it take four times as long. This is the basic law of diffusion: “traveling a distance x takes a time t ∼ x²/κ².” Conversely, the depth to which heat has soaked in by time t is √(κ²t). We will use this later.

τ for a 2 cm Thick PlateChanging the material changes only “how fast the clock runs”
Material	τ = a²/(π²κ²)	Ratio to copper
Copper	347 ms	1
Aluminum	416 ms	1.2
Iron	1.72 s	5
Air	2.19 s	6
Ice	34.6 s	100
Glass	101 s	292
Water	293 s	844

7. Beyond the Notes — Solving the Surroundings Too

So far the boundary condition has been “hold the surface of the plate at T_s.” In reality, however, the temperature of the surrounding air changes too. So we set up the same Eq. (4) for x < 0 and x > a as well. Only K, ρ, and c differ.

The boundary condition is no longer “specify the temperature”; it becomes a condition that joins the two regions.

\begin{array}{ll} Temperature is continuous & T_{1} (0^{-}) = T_{2} (0^{+}) \\ Heat flux is continuous & K_{1} {(\frac{\partial T}{\partial x})}_{0^{-}} = K_{2} {(\frac{\partial T}{\partial x})}_{0^{+}} \end{array}

(13)

The first says “no temperature jump is left at the interface,” and the second says “no heat accumulates at the interface.” Because K differs, it is the slope of the temperature that bends. This is the kink in the line at the interface (dashed line) in the simulator’s figure.

Two semi-infinite bodies can be solved by hand

As long as the plate is thick enough that the influence of the opposite face has not yet arrived (t ≪ a²/κ²), the problem becomes that of two semi-infinite bodies brought into contact at x = 0. This can be solved exactly.

For a semi-infinite medium whose surface is held at T_c, the solution can be written with the error function.

T (x, t) = T_{c} + (T_{i} - T_{c}) \erf \frac{| x |}{2 \sqrt{κ^{2} t}}

Differentiating this, the magnitude of the heat flux through the surface is

| q | = K | {(\frac{\partial T}{\partial x})}_{0} | = \frac{\sqrt{K ρ c} | T_{i} - T_{c} |}{\sqrt{π t}}

— What appeared is neither K nor κ² but the combination √(Kρc). If we choose this T_c so that the heat fluxes on the two sides balance (the second condition of Eq. (13)), we get

T_{c} = \frac{b_{1} T_{1} + b_{2} T_{2}}{b_{1} + b_{2}}, b \equiv \sqrt{K ρ c}

(14)

b is a quantity of “strength” The thermal effusivity b = √(Kρc) is a completely different quantity from κ² = K/(ρc). κ² is “how fast heat spreads,” while b is “how strongly a material can exchange heat.”
The meaning of Eq. (14) is striking — the interface temperature is fixed the instant the bodies touch, and stays unchanged for some time afterward, because t does not appear anywhere in it. Moreover, its value is pulled toward the body with the larger b.

Fig. 2: Bringing two bodies into contact. The kink in the line at the interface (x = 0) is Eq. (13). The green dashed line is T_c of Eq. (14), and it does not move as time passes. If you choose copper and air, the interface sticks almost to the temperature of the copper, and the air is dragged along by the copper. With copper and water, the b of water is much larger than that of air, so the interface drops more noticeably than it does with air.

Left body (100 °C) Right body (0 °C) Time 1 s

Air is “fast but weak”

Look at the table once more. The κ² of air is about the same as that of iron (1.9×10⁻⁵ versus 2.4×10⁻⁵ m²/s). Yet b is 5.6 versus 17182 — 1/3000. Because its density is only 1/6000 as large, it can hardly store or carry any heat.

So if a copper plate at 100 °C is placed in air, Eq. (14) gives

T_{c} = \frac{37285 \times 373 + 5.6 \times 293}{37285 + 5.6} = 372.99 K

The copper hardly cools at all. It is the air that cools. (In reality convection comes into play, so it is not this extreme. See §9.)

Why metal feels cold The b of skin is about 1000. When you touch a metal at 20 °C (b = 37000), Eq. (14) shows that the interface sticks almost to the metal side, and the surface of the skin is pulled down to nearly 20 °C in an instant. With wood at the same 20 °C (b ≈ 400), the interface lies closer to the skin side, so it hardly feels cold. What you feel by touch is not the other object’s temperature but its b. — This is also why we need thermometers.

8. How to Solve It Numerically

When the plate and its surroundings are combined, K jumps at the interface, so a neat series solution cannot be written. We solve it numerically. We used the finite volume method — instead of discretizing the differential equation, we write down the energy balance of each cell as it is, and this way the conservation law is automatically respected.

{(ρ c)}_{i} h \frac{d T_{i}}{d t} = K_{i + 1 / 2} \frac{T_{i + 1} - T_{i}}{h} - K_{i - 1 / 2} \frac{T_{i} - T_{i - 1}}{h}

(15)

The K at a cell face is taken as the harmonic mean.

K_{i + 1 / 2} = \frac{2 K_{i} K_{i + 1}}{K_{i} + K_{i + 1}}

Why not the arithmetic mean? Two neighboring cells have thermal resistances h/(2K) connected in series. For a series connection it is the resistances that add. So K is averaged through its reciprocal — the harmonic mean. With the arithmetic mean, when K differs by a factor of about 17000, as for copper and air, the resistance on the air side is effectively ignored and too much heat flows. With the harmonic mean, the smaller K properly becomes the rate-limiting one.

Time is advanced by the backward Euler method. Writing Eq. (15) with Tⁿ⁺¹ as the unknown gives a tridiagonal system of simultaneous equations, which the Thomas algorithm solves in O(N) operations. With forward Euler, the stability condition κ²Δt/h² ≤ 1/2 would bind us, and for copper Δt would have to be absurdly small. Backward Euler is unconditionally stable, so there is no such worry.

Checking the Numerical SolutionResults of extracting the shipped code as it is and running it
What was checked	Result
Agreement with the series solution (11) (copper, glass, water; both initial distributions)	Maximum difference 0.008 K (0.08 % of the peak)
Agreement of the interface temperature with Eq. (14) (5 combinations)	Difference 0.002–0.13 K
Agreement with the erf solution	Maximum difference 0.05–0.29 K
Energy conservation in an insulated box	Relative error 10⁻⁸ or less
Final temperature reached	Within 10⁻⁴ K of the average weighted by the volumetric heat capacity

9. Putting Convection into the Boundary Condition — Newton’s Law of Cooling

The boundary condition so far has been “hold the surface at T_s.” When a plate is actually placed in air, no one is clamping down the surface temperature. What happens is this: the heat carried by conduction from inside the plate to the surface is carried straight off into the outside air. This handoff is described by Newton’s law of cooling, which assumes that the amount carried away is proportional to the temperature difference, and since no heat accumulates at the surface, the following condition results.

- K \frac{\partial T}{\partial x} = h (T - T_{\infty})

(16)

The left side is the heat flux arriving by conduction, and the right side is the heat flux leaving by convection. What differs from before is that it imposes not the temperature itself but a relation between the temperature and the slope.

9.1 The eigenvalue problem changes

This condition changes the eigenvalue problem that follows separation of variables completely. Since both the initial distribution and the boundary condition are symmetric about x = a/2, let us take the center as the origin, ξ = x − a/2, and put the half thickness L = a/2. Only symmetric solutions remain, so we need to consider only X(ξ) = cos(λξ), and substituting this into Eq. (16) at ξ = L gives

K λ \sin λ L = h \cos λ L ⟹ μ \tan μ = Bi, μ \equiv λ L, Bi \equiv \frac{h L}{K}

(17)

This is where it differs most from before ωa = nπ could be solved immediately just by comparing the two sides. Eq. (17), however, is a transcendental equation, and μ_n cannot be written in closed form. The only option is to find them numerically. Still, the structure is clear: on each interval ((n−1)π, (n−1)π + π/2), μ tan μ increases monotonically from 0 to +∞, so there is exactly one root in each such interval. The simulator finds them by the bisection method and shows μ₁, μ₂, and μ₃ in its readout.

Fig. 3: Roots of the transcendental equation μ tan μ = Bi. The blue curve is y = μ tan μ and the green horizontal line is y = Bi. The intersections are the eigenvalues μ_n. As Bi is increased, the roots approach the asymptotes on their right, (2n−1)π/2 (gray dashed lines), and return to the answer in the lecture notes. As Bi is decreased, μ₁ falls toward 0, and only the slowest mode lives extremely long — this is the lumped-capacitance system. In the readout below, see which of the two formulas the time constant approaches in each case.

Biot number Bi = 1 Plate material Thickness a = 20 mm

The expansion coefficients, as usual, follow from orthogonality. Noting that the denominator ∫cos²(λξ)dξ equals L(2μ + sin 2μ)/(2μ), and putting θ = T − T_∞ and θ₀ = T₀ − T_∞, we get

\begin{array}{l} Uniform θ_{0} : C_{n} = \frac{4 \sin μ_{n}}{2 μ_{n} + \sin 2 μ_{n}} \\ Parabola θ_{0} (1 - ξ^{2} / L^{2}) : C_{n} = \frac{8 (\sin μ_{n} / μ_{n}^{2} - \cos μ_{n} / μ_{n})}{2 μ_{n} + \sin 2 μ_{n}} \end{array}

(18)

T (x, t) = T_{\infty} + θ_{0} \sum_{n = 1}^{\infty} C_{n} e^{- μ_{n}^{2} κ^{2} t / L^{2}} \cos \frac{μ_{n} ξ}{L}

9.2 In two limits we return to known answers

That this solution is correct can be checked by taking the limits at both ends.

The Two Limits of the Biot NumberBi = hL/K = (thermal resistance inside the plate) ÷ (thermal resistance beyond the surface)
	Bi → ∞ (the surface can take away heat without limit)	Bi → 0 (the surface is the bottleneck)
Eigenvalues	μ_j → (2j−1)π/2	μ₁² → Bi
Coefficients (uniform)	4(−1)^j−1/((2j−1)π)	C₁ → 1, the others 0
Coefficients (parabola)	32(−1)^j−1/((2j−1)³π³)	C₁ → 2/3
Time constant	τ → a²/(π²κ²)	τ → ρcL/h
What it really is	Exactly the problem solved in §5	A lumped-capacitance system (the whole plate cools exponentially while staying uniform)

The Bi → ∞ coefficients carry the factor ${(- 1)}^{j - 1}$ but the answer of §5, $4 T_{0} / n π$ (n odd), carries none. This is merely a difference in the choice of basis. §5 uses a sin basis with its origin at the edge, and here we use a cos basis with its origin at the center. With the correspondence $n = 2 j - 1$ , the relation $\sin (n π x / a) = {(- 1)}^{j - 1} \cos (μ_{j} ξ / L)$ holds, so the sign has merely moved, and the summed temperature distribution is the same.

What “hold the surface at T_s” really is Look at the Bi → ∞ column. The problem solved in the lecture notes was the limit of “infinitely strong convection.” The surface temperature can be fixed only when the outside accepts as much heat as needed. If h is finite, the surface does not drop all the way, and the temperature difference inside the plate also becomes smaller.
At the opposite end (Bi → 0), conduction inside the plate is so fast that the temperature difference disappears, and the problem becomes an ordinary differential equation rather than a partial differential equation — ρca dT/dt = −2h(T − T_∞). Note that κ² does not appear in τ. What limits the rate is not the inside of the plate.

Bi is the ratio of “the thermal resistance inside the plate, L/K,” to “the thermal resistance beyond the surface, 1/h.” Which of the two is larger determines the shape of the temperature distribution — it is the same story as in an electric circuit: the temperature difference falls across the larger resistance.

The initial condition need not satisfy the boundary condition (but there is a price) In §5, the parabolic initial distribution happened to equal T_s exactly at the edges. Eq. (16), however, is a “relation between the slope and the temperature difference,” and neither the uniform distribution nor the parabola satisfies it. The expansion is still possible — because the eigenfunctions form a complete set. However, the coefficients decay only as 1/μ², and the convergence is slow only at the edges at t = 0. After even a little time, exponential decay takes over and the slow convergence disappears at once.

10. Radiation, and What Happens in an Electric Furnace

At high temperatures radiation matters more than convection. The radiative heat flux is εσ(T⁴ − T_∞⁴) (σ = 5.670×10⁻⁸ W/(m²·K⁴), and ε is the emissivity). As the fourth power it is not in the form of Newton’s law of cooling, but factoring it gives

ε σ (T^{4} - T_{\infty}^{4}) = \underset{h_{rad}}{\underset{⏟}{ε σ (T^{2} + T_{\infty}^{2}) (T + T_{\infty})}} (T - T_{\infty})

(19)

This is not an approximation but an identity (we simply used A⁴−B⁴ = (A²+B²)(A+B)(A−B) as it stands). However, since h_rad depends on T, the boundary condition becomes nonlinear. Superposition then fails, and neither separation of variables nor the Fourier series can be used — the answer here is that we have no choice but to solve it numerically. At each time step, the simulator rebuilds h_rad from the current surface temperature.

Radiation Grows Roughly as T³ For a surface with ε = 0.9 (ceramics, painted surfaces) facing room temperature, 293 K
Surface temperature	h_rad [W/(m²·K)]	Natural convection h	Which dominates
50 °C	6.0	5–10	About the same
100 °C	7.7	5–10	About the same
500 °C	37	5–10	Radiation
1000 °C	136	5–10	Almost entirely radiation

At 1000 °C, h_rad = 136 against 5–10 for natural convection — even taking h = 10, radiation alone carries 14 times as much heat as convection, and 27 times as much if h = 5. At firing temperatures, what determines how the sample and the furnace cool is radiation, not convection.

A red-hot object can be “seen” because part of its radiation reaches the visible range. Most of the energy carried away is infrared (at 1000 °C the peak is around 2.3 μm), but if the object is hot enough to be seen by visible light, radiation as a whole has become the main player in its cooling — that is the telltale sign. Polished metal has only ε ≈ 0.03, so even at the same temperature it hardly cools by radiation (this is why the inside of a thermos flask is a mirror).

10.1 Reading the three boundary conditions through an electric furnace

Every boundary condition that has appeared so far corresponds to a situation that actually occurs around an electric furnace for firing ceramics.

Boundary Conditions and the Electric FurnaceThey differ only in how heat is handled at the outer edge
Boundary condition	Equation	In terms of an electric furnace
Insulated	∂T/∂x = 0	Inside the furnace. Since heat does not escape, the sample and the furnace wall settle to the average temperature weighted by the heat capacity
Fixed temperature	T = T_∞	An infinite heat bath. Quenching in a water tank, or placing the sample on a large metal surface plate. In the language of §7, the case where the other body’s b is larger by orders of magnitude
Newton’s law of cooling	−K∂T/∂x = h(T−T_∞)	Taking the sample out of the furnace / heat escaping through the furnace wall to room temperature. Real cooling is almost always this

Correspondence with the simulator’s symbols The T_∞ on this page is the outside temperature seen by the boundary condition. The simulator’s ② “Solve the surroundings too” also solves the material around the plate (air, water, firebrick) itself, so the initial temperature of the surroundings T_∞ and, beyond the outer edge of that material, the outside temperature T_out are separate inputs. When a furnace cools, the furnace wall stays hot (T_∞ is high) and only its outside is at room temperature (293.15 K). The T_∞ in the table above and in Eq. (16) corresponds to T_out in ②.

Why furnace walls are thick The K of firebrick is about 0.3 W/(m·K). For a thickness d = 50 mm, the thermal resistance is d/K = 0.17 m²·K/W, which is larger than that of the outer surface (1/h ≈ 0.05–0.1). Most of the resistance is inside the wall, so making the wall thicker pays off directly. Conversely, as the wall is made thinner, the thermal resistance falls to the surface’s 1/h, and levels off there (for a temperature difference ΔT, the heat flux is at most hΔT) — the usual story of resistances connected in series.

Bi is also the reason slow cooling (annealing) is needed Fired ceramics and glass crack if they are cooled suddenly. The cause is not the temperature itself but the temperature gradient — if only the surface shrinks while the inside does not, tensile stress arises. What determines the size of the gradient is Bi = hL/K.
The K = 0.75 of glass is only 1/40 of the value 30 for alumina, so even with the same h its Bi is 40 times larger. The temperature difference between the inside and the surface is correspondingly larger. This is why glass needs long annealing, whereas a metal (Bi ≪ 1) cools uniformly to the inside and rarely cracks from thermal shock alone — though in the quenching of steel, because of the volume change accompanying the martensitic transformation and residual stresses, quench cracks can occur.
In other words, slow cooling is the operation of keeping h small to lower Bi, so that the temperature falls while the inside of the plate stays uniform.

11. Summary

Do not confuse the three quantities κ² = K/(ρc) (thermal diffusivity) — determines the time. τ = a²/(π²κ²), and the depth to which heat soaks in, √(κ²t).
b = √(Kρc) (thermal effusivity) — determines the interface temperature just after contact. T_c = (b₁T₁ + b₂T₂)/(b₁ + b₂). What determines where the system settles after a sufficiently long time is not b but the volumetric heat capacity ρc.
Bi = hL/K (Biot number) — determines the shape of the temperature distribution. Bi ≪ 1 gives a uniform plate; for Bi ≫ 1 the surface cools first.

Air has a κ² as large as iron’s, yet its b is 1/3000 — fast but weak. Once you can think of these three separately, “why metal feels cold,” “why insulation traps air,” “why thick walls work,” and “why glass needs slow cooling” all come out of the same single equation.

When the boundary condition changes, the whole eigenvalue problem changes This is probably the one thing I would most like you to take away from this material. With “hold the surface at T_s,” we could solve immediately to get ωa = nπ. Merely changing that to Newton’s law of cooling turns it into μ tan μ = Bi, a transcendental equation that cannot be written in closed form. Adding radiation makes it nonlinear, and the method of separation of variables itself can no longer be used.
Even when the equation is the same, it is the treatment of the boundary with the outside that determines how it is solved.

This is the explanation page for the thermodynamics teaching material of the Mochizuki Research Group (heat conduction equation). §1–§6 follow the lecture notes. §7 onward — solving the surroundings too, thermal effusivity, Newton’s law of cooling and radiation, and electric furnaces — goes beyond them.