Transient heat conduction

Governing equation

Heat diffusion, one dimension, no generation

∂T/∂t = α ∂²T/∂x²

Thermal diffusivity

α = k / (ρ c_p) [m² s⁻¹]

• k — thermal conductivity, W m⁻¹ K⁻¹: how much heat flows for a given gradient.
• ρ c_p — volumetric heat capacity, J m⁻³ K⁻¹: how much heat it takes to move a cubic metre one kelvin.
• α — their ratio: how fast a temperature shape spreads. It says nothing about how much heat flows.

Boundary conditions offered here

T(0, t) = T_L

T(L, t) = T_R | ∂T/∂x|_L = 0 | −k ∂T/∂x|_L = h (T_L,end − T_∞)

Dirichlet, Neumann (adiabatic) and Robin (convective) respectively.

Dimensionless groups

Fo = α t / L² · Bi = h L / k · Fo_mesh = α Δt / Δx²

• Fo — dimensionless time. Fo ≈ 1 means the disturbance has crossed the rod; every material behaves identically when plotted against Fo. That is the whole content of the Heisler charts.
• Bi — surface resistance versus internal resistance. Bi < 0.1 means the internal gradient is negligible and a single lumped thermal capacitance is a valid model; above that you have to keep the profile.
• L²/α — the diffusion time scale, the natural clock of the problem. It is 13 minutes for a 0.3 m copper rod and 8 days for the same rod in oak.

Steady state

d²T/dx² = 0 ⇒ T(x) is a straight line

T_end = (T_L + Bi · T_∞) / (1 + Bi)

Derived by equating the conducted flux k(T_L − T_end)/L with the convected flux h(T_end − T_∞) and solving for T_end. As Bi → 0 the end sits at T_L; as Bi → ∞ it is pinned to T_∞.

Exact solution, two fixed ends, uniform start

T(x,t) = T_ss(x) + Σ_n≥1 B_n sin(nπx/L) · exp(−α (nπ/L)² t)

B_n = (2/nπ) [ (T_i − T_L)(1 − (−1)ⁿ) + (T_R − T_L)(−1)ⁿ ]

Separation of variables on u = T − T_ss, which has homogeneous ends. Mode n decays with time constant L²/(α n²π²), so the wrinkles vanish first and the last thing left is the half-sine — which is why almost every transient conduction chart is a one-term approximation for Fo > 0.2.

How it is solved

Space is divided into 200 control volumes with 201 nodes, so Δx = L/200. Interior nodes get the usual central second difference. Both end nodes sit on the boundary and own a half-width control volume — the detail that matters.

Energy balance on the boundary node's half control volume

ρ c_p (Δx/2) dT_N/dt = k (T_N−1 − T_N)/Δx − h (T_N − T_∞)

The ρ c_p (Δx/2) term is the boundary node's own thermal capacitance. Imposing the Robin condition through a mirrored ghost node instead drops that term, which is a first-order error concentrated exactly where the interesting physics is.

Crank–Nicolson (θ = ½) in time

(T_iⁿ⁺¹ − T_iⁿ)/Δt = (α/2)(δ²Tⁿ⁺¹ + δ²Tⁿ)/Δx²

Each step is a tridiagonal system, solved by the Thomas algorithm in O(N). The matrix is strictly diagonally dominant, so no pivoting is needed.

The first law, and the two flux estimates

The conservation check

E(t) − E(0) = Q_in − Q_out, E = ρ c_p A ∫₀^L (T − T_ref) dx

q″ = −k ∂T/∂x ≈ −k (−3T₀ + 4T₁ − T₂) / (2Δx)

The stored energy is a trapezoid integral over the profile; the heats are time integrals of the end fluxes. Nothing links them in the code, so if the solver were wrong they would drift apart. The residual readout is |E − E₀ − (Q_in − Q_out)| / max(|Q_in|, 1).

Two different flux estimates live on this page, and the difference between them is worth understanding:

• The displayed q″ at each end is the second-order one-sided difference above — the better estimate of the true surface gradient, and the number an engineer would quote.
• The ledger integrates the control-volume face flux k(T₀ − T₁)/Δx. That is the flux the discrete scheme actually conserves, so the balance closes to rounding error rather than to discretisation error.

Their difference is exactly (k/2Δx)(T₀ − 2T₁ + T₂) ≈ ρc_p(Δx/2)·dT₁/dt — half a cell's worth of stored energy. The "gap between the two q″ estimates" readout shows it live: large in the first instants, when the thermal layer is thinner than one cell, and a few hundredths of a percent thereafter. Both ledgers are exposed in the page's test hook (residual and residualSurface); neither is hidden.

Material properties

Incropera & DeWitt, Fundamentals of Heat and Mass Transfer, Tables A.1, A.3 and A.6 at 300 K. Note the fibreglass blanket: 37× less conductive than plate glass, yet almost twice the diffusivity, because it has almost no mass to heat. Conductivity governs the steady heat flow; diffusivity governs the transient.

What to try

• The diffusivity race. Material Copper, ghost Oak, L = 0.3 m, left end 200 °C, right end fixed at 20 °C. Copper's α is 891× oak's, so copper reaches Fo = 1 in about 13 minutes and oak needs about 8 days. On the shared clock the oak line barely leaves 20 °C while copper has already become a straight line. Read the two "diffusion time L²/α" figures and check that their ratio is the α ratio.
• Low conductivity is not the same as slow. Set material Fibreglass blanket and ghost Glass, plate. The blanket carries 37× less heat at steady state — check q″ at x = 0 once the profile is a straight line — yet its shape settles about twice as fast, because ρc_p is tiny. This is why insulation works and why it warms up so quickly to the touch.
• Insulate the far end. Right end → Insulated. The slope at x = L collapses to zero, Q_out stays exactly 0, and the whole rod climbs to T_L. On the ledger plot the stored-energy line and the Q_in line lie on top of each other — with nothing leaving, everything that enters is stored.
• Walk the Biot number. Right end → Convection, T_∞ = 20 °C, copper, L = 0.3 m. At h = 5 W/m²·K, Bi ≈ 0.0037: the surface resistance dominates completely, the rod is nearly isothermal at 200 °C, and the "lumped model valid" flag is on. Slide h up to 10 000: Bi ≈ 7.5, the end is dragged down near 40 °C and the profile becomes a steep straight line. Check the far-end reading against the "Steady far-end T" figure at each end of the sweep.
• Beat the explicit-scheme limit. Watch the two "mesh Fourier number" readouts. The requested one sits at about 67; the integrated one at about 3.9 — still nearly eight times the 0.5 ceiling an explicit FTCS scheme needs to stay stable, and the solution is perfectly smooth. That gap is the accuracy clamp doing its work: the solver sub-divides seventeen times per displayed step, not because it would go unstable — Crank–Nicolson cannot — but because its amplification factor approaches −1 for the shortest mesh waves and a very large step would ring on the step initial condition.
  Now change the material or the rod length and watch both numbers hold. That is not a coincidence and it is not the solver: the clock itself runs in Fourier number of the primary rod, so Δt = τ/600 and Fo_mesh = αΔt/Δx² = 40000/600 for any α and any L. The discrete problem is scale-invariant in Fo, which is also why the series check below reads bit-identically for copper and for oak.
  To see the requested number actually move, put a different material in the ghost rod: its Fourier number scales with its diffusivity ratio to the primary, and past about 1.44 its integration scheme changes — the "Its integration scheme" readout says when.
• Check the numerics against the exact answer. With two fixed ends, tick Show the Fourier-series solution. The dotted green curve is the closed form; it should sit on the numerical line at every time. The "vs Fourier-series solution" readout gives the worst deviation as a fraction of the temperature span — expect a few hundredths of a percent from Fo ≈ 0.01 onward. Change the material and it does not move: the two solutions agree because they are the same problem, not because one was tuned to the other.

Assumptions and limits

• Strictly one-dimensional. No lateral heat loss, no radial gradient. Real bar stock in still air loses heat along its length; if that matters, the fin equation replaces this one.
• Constant properties. k, ρ and c_p are frozen at their 300 K values. Copper's conductivity falls about 10 % between 300 K and 700 K; AISI 304's rises about 40 % over the same range. Once the swing exceeds a couple of hundred kelvin the equation becomes non-linear and this model is only indicative.
• No radiation. Above roughly 200 °C an exposed end radiates as much as it convects into still air. Here h is a constant, temperature-independent number, whereas a real natural-convection coefficient scales roughly as ΔT^1/4.
• No internal generation, no phase change. Add ohmic heating (a busbar under fault current), latent heat (a freezing pipe) or a moving front and the equation gains terms this solver does not have.
• Perfect contact at the ends. The fixed-temperature end is a zero-resistance joint to an infinite reservoir. A real bolted or clamped joint has a thermal contact resistance around 10⁻⁴ to 10⁻³ m²·K/W, which behaves like an extra Robin condition in series.
• Water is here as a conduction-only medium. Its diffusivity is right, but any still-water layer thicker than a few millimetres heated from the side will convect long before conduction matters. Treat it as a reference number, not a model of a water column.
• A step initial condition. The rod is uniform and the end temperature is imposed instantly at t = 0. Real boundaries ramp. During the first few instants — while the thermal layer is thinner than one cell, Δx — the surface gradient is the least trustworthy number on the page, which is exactly what the "gap between the two q″ estimates" readout is showing you.
• Discretisation. 201 nodes, second order in space and (at θ = ½) in time. Reported to four significant figures throughout, which is more precision than the material data deserves — the property values themselves are good to two or three figures at best.

Why an engineer cares

This one equation carries a surprising amount of routine practice. Short-circuit heating of a busbar or a cable conductor is this problem with a generation term, and the adiabatic limit of it is where the familiar I²t conductor-sizing rules come from. Concrete pour temperature control, quench schedules in heat treatment, and the fire rating of an insulated steel column are all 1-D transient conduction with a prescribed surface temperature or a convective face. In electrical work specifically, switchgear enclosure warm-up and transformer winding hot-spot models are almost always built first as lumped thermal capacitances and resistances — and the Biot number readout on this page is precisely the test of whether that shortcut is allowed. Bi < 0.1 and the lumped model is defensible; above it, you owe the reviewer a profile.