Transient Conduction: Method

Transient Conduction

When the surroundings of a solid change suddenly, its temperature does not follow at once: heat has to diffuse in from the surface. This tool answers the usual questions of that process (how hot is the centre after a given time, how long until a point reaches a temperature, how much heat has been exchanged) by five analyses that share one set of inputs. All of them assume constant properties and no change of phase. The thermal quantities are entered in SI units in both unit systems (k in W/(m⋅K), ρ in kg/m³, c_p in J/(kg⋅K), h in W/(m²⋅K), heat flux in W/m², temperatures in °C, time in s); lengths follow the unit system.

Symbols

α = k/(ρ c_p) is the thermal diffusivity. For a body of half-thickness or radius L the Biot number Bi = h L/k compares the resistance inside the body with the resistance of the surface film, and the Fourier number Fo = α t/L² is the dimensionless time. θ* = (T − T_∞)/(T_i − T_∞) is the fraction of the initial temperature difference still present: 1 at the start, 0 when the body has reached the fluid temperature. Cooling and heating are the same problem.

Lumped capacitance

If the body conducts so well that it is nearly at one temperature, an energy balance on the whole body gives ρ c_p V dT/dt = −h A (T − T_∞), whose solution is the exponential θ* = exp(−t/τ) with the time constant τ = ρ c_p V/(h A). The model holds when the Biot number formed with the length L_c = V/A is below about 0.1 (L_c is L for the wall, L/2 for the cylinder, L/3 for the sphere); the tool prints the check and, beside the lumped value, the exact centre, surface and mean temperatures of the same body so that the error is visible.

With uniform internal generation q‴ the balance stays linear, θ = θ_i exp(−t/τ) + (q‴V/(hA))(1 − exp(−t/τ)), and the body tends to T_∞ + q‴V/(hA). With radiation to large surroundings the term ε σ A (T⁴ − T_sur⁴) makes the equation non-linear; the tool integrates it by the fourth-order Runge–Kutta method with a step of one fiftieth of the shortest time constant. For radiation alone the equation separates, and the time to cool from T_i to T is the closed form (absolute temperatures)

t = [ρ c_p V/(4 ε σ A T_sur³)] { ln|(T_sur + T)/(T_sur − T)| − ln|(T_sur + T_i)/(T_sur − T_i)| + 2 [tan⁻¹(T/T_sur) − tan⁻¹(T_i/T_sur)] },

which is printed as a comparison when a target temperature is entered. The optional second node is an inner mass of capacity C₂ joined to the body by a conductance G (a core, a fitting, a filling): the two linear equations have two exponential modes, a fast one in which the two masses equalise and a slow one in which they cool together, and the tool gives both time constants and both temperatures. Radiation and generation are not carried into the two-node model.

Plane wall, long cylinder and sphere: the exact series

For a wall cooled on both faces, a long cylinder or a sphere, initially uniform and exposed to a fluid with a constant coefficient, separation of variables gives

θ* = Σ C_n exp(−ζ_n² Fo) f_n(x/L),

Body	Roots ζn of	Cn	fn	Volume mean of fn
Plane wall	ζ tan ζ = Bi	4 sin ζ/(2ζ + sin 2ζ)	cos(ζ x/L)	sin ζ/ζ
Long cylinder	ζ J1(ζ)/J0(ζ) = Bi	(2/ζ) J1(ζ)/(J0²(ζ) + J1²(ζ))	J0(ζ r/L)	2 J1(ζ)/ζ
Sphere	1 − ζ cot ζ = Bi	4 (sin ζ − ζ cos ζ)/(2ζ − sin 2ζ)	sin(ζ r/L)/(ζ r/L)	3 (sin ζ − ζ cos ζ)/ζ³

The roots are bracketed between the poles and zeros of the characteristic function (for the cylinder between successive zeros of J₁ and J₀) and found by bisection to the last bit; the Bessel functions are evaluated from Bessel's integral below an argument of 25 and from Hankel's asymptotic series above it. At least fifty terms are summed, and more as Fo falls (the number is chosen so that the last term is below 10⁻¹³, up to 3,000), so the series is used at any time for which the body is not better treated as semi-infinite. The heat transferred up to time t, as a fraction of Q_o = ρ c_p V (T_i − T_∞), is 1 minus the volume mean of θ*. The time at which a point reaches a given temperature is solved by bisection on the series, which is safe because the temperature at every point moves monotonically towards the fluid temperature.

One-term approximation. All terms after the first die out quickly, and for Fo > 0.2 the first term alone, θ* ≈ C₁ exp(−ζ₁² Fo) f₁, is within 2 % at the centre (the largest errors at Fo = 0.2 are 1.8 % for the wall, 1.6 % for the cylinder and 1.0 % for the sphere, over all Biot numbers). The tool computes ζ₁ and C₁ for the entered Biot number instead of interpolating a table, prints the one-term values beside the series and flags them when Fo is below 0.2. At Bi = 1 the values are ζ₁ = 0.8603, C₁ = 1.1191 (wall), 1.2558 and 1.2071 (cylinder), 1.5708 and 1.2732 (sphere).

Charts. Heisler's charts (centre temperature against Fo for lines of 1/Bi, and a position correction) and Gröber's chart (Q/Q_o against Bi² Fo) were the working form of this solution before computers. The tool redraws all three from the full series for the chosen body, adds the line of the entered Biot number and marks the case. The position-correction chart is the first shape function and therefore shares the Fo > 0.2 limit; the other two are exact at every time drawn.

Check. Each run of the series is accompanied by an independent finite-difference solution of the same problem (200 cells, Crank–Nicolson), printed in the temperature table. A wall insulated on one face is the same problem with L its full thickness.

Semi-infinite solid

Early in any transient, and for all time in a very thick body, only a layer near the surface has changed and the far side plays no part. With η = x/(2√(α t)):

• Surface temperature stepped to T_s: (T − T_s)/(T_i − T_s) = erf η; surface flux k (T_s − T_i)/√(π α t); heat absorbed per unit area 2 k (T_s − T_i) √(t/(π α)).
• Constant surface flux q_s: T − T_i = (2 q_s/k) √(α t/π) exp(−η²) − (q_s x/k) erfc η.
• Surface convection: (T − T_i)/(T_∞ − T_i) = erfc η − exp(h x/k + h² α t/k²) erfc(η + h √(α t)/k). The second term is computed as exp(−η²) times the scaled function exp(z²) erfc z, which stays finite however large the time.
• Two solids brought into contact: the interface jumps to T_s = (e_A T_A + e_B T_B)/(e_A + e_B) with e = √(k ρ c_p) the effusivity, and stays there; each solid follows the first solution. This is why a metal at room temperature feels colder than wood.
• Periodic surface temperature T_m + A sin ω t (sustained state): T = T_m + A exp(−x/d) sin(ω t − x/d) with the damping depth d = √(2α/ω). The wave is reduced by e and delayed by 1/(2π) of a period for every d of depth; 4√(α/ω) is the usual penetration depth. The depth below which the temperature never drops under a threshold T_f is d ln[A/(T_m − T_f)]. With the yearly wave and 0 °C this is a first estimate of frost depth; it ignores the latent heat of freezing soil water, snow cover and the difference between air and ground-surface temperature, so it is not a design frost depth.

The depth at which the change is 1 % of a step at the surface is 3.643 √(α t) (10 % at 2.326 √(α t)); a wall thicker than that can be treated as semi-infinite up to the time t.

Product solutions

A short cylinder is the intersection of a long cylinder and a wall; a bar is the intersection of two walls and a brick of three. When the initial temperature is uniform and every face sees the same fluid, the solution in two or three dimensions is the product of the one-dimensional solutions, each with its own Bi and Fo: θ* = θ*₁ θ*₂ (θ*₃). The tool lists the factors and the temperature at the centre, at the middle of each face and edge and at the corner. The fraction of heat transferred is 1 − Π(mean of θ*_i), which is Langston's rule (Q/Q_o)₁ + (Q/Q_o)₂[1 − (Q/Q_o)₁] + … written compactly. The product rule does not hold for different coefficients on opposite faces of one direction, for a non-uniform start or for a prescribed flux.

Multilayer wall by finite differences

Each layer is divided into equal cells with a node on every cell boundary, so that a node sits on each face and each interface. The energy balance of the control volume around node i is C_i dT_i/dt = K_i−1(T_i−1 − T_i) + K_i(T_i+1 − T_i) (+ the boundary term), where K = k/Δx of the cell between two nodes and C_i is the heat capacity of the two half cells around the node; an interface node therefore carries half a cell of each material and no special interface equation is needed. A face may convect, be held at a temperature, receive a heat flux or be insulated.

The explicit scheme evaluates the right-hand side at the old time. It is stable only if the time step does not exceed C_i/ΣK at every node, which for an interior node is the familiar Δx²/(2α) (cell Fourier number ½) and for a convecting face Δx²/(2α(1 + h Δx/k)); the tool prints the limit, takes 0.8 of it by default and reduces a larger entered step. Thin, conductive layers make the limit very small, and the run is refused above 150,000 steps. The Crank–Nicolson scheme averages the right-hand side between the old and the new time, is stable for any step and second-order accurate; the tool starts it with two fully implicit steps so that a step change at a face does not ring. Both schemes conserve energy exactly, which the test script checks.

When the wall is a single layer with one insulated and one convecting face the tool prints the exact series beside the finite-difference result. The default (25 mm of steel, the half-plate of the default quench) agrees to 0.005 °C with 40 cells.

Material list

The list fills k, ρ and c_p with typical room-temperature values from two public documents: Table 2-3 of NUREG-1805 (U.S. Nuclear Regulatory Commission, 2004; the table credits Klote and Milke) for building materials, pure aluminum and carbon steel, and the room-temperature physical-property lines of the MIL-HDBK-5J design tables for the alloys whose density, specific heat and conductivity are printed as numbers (converted from lb/in³, Btu/(lb⋅°F) and Btu/(h⋅ft⋅°F)). Properties change with temperature, strongly for steels above a few hundred degrees; for a wide temperature range enter values at the mean temperature of the process. The convection coefficient is also taken constant, which a real quench (film boiling, then nucleate boiling) does not satisfy; the default case is illustrative.

Limits

Constant properties and coefficient; no phase change, no contact resistance between layers; uniform initial temperature (except the periodic case); the same coefficient and fluid temperature on every face in the series and product analyses; one-dimensional conduction in the multilayer wall. The method was checked against Tables 4.19, 4.20 and 9.7 of Abramowitz and Stegun (roots of the three characteristic equations), against scipy, against a method-of-lines solution and against two-dimensional numerical solutions of a bar and a short cylinder.