Numerical Methods in Heat Transfer

Numerical heat conduction

Discretization turns the heat equation into algebraic equations. Finite differences approximate derivatives on a grid, finite volumes balance energy over cells, and finite elements use a weak formulation and interpolation over elements. Boundary conditions, material properties and discretization errors matter as much as the solver.

Explicit one-dimensional scheme

For a uniform grid, constant thermal diffusivity alpha = k/(rho c), no heat generation and one-dimensional conduction:

T_i^(n+1) = T_i^n + Fo_step [T_(i-1)^n - 2T_i^n + T_(i+1)^n],

where Fo_step = alpha Δt/(Δx)². Subscript i is spatial position and superscript n is time level. Update all new values using only the old time level; do not overwrite old temperatures while stepping through the nodes.

For this simple explicit interior scheme, stability requires Fo_step <= 1/2. Boundary treatments and multidimensional grids can impose different limits. A stable result is not necessarily accurate: check mesh and time-step convergence.

Worked example

A 1 m rod has k = 200 W/(m K), rho = 7800 kg/m³ and c = 500 J/(kg K). Its lateral surface is insulated. At t = 0 the interior is uniformly 100°C and both ends are suddenly set to 0°C. Use Δx = 0.1 m and Δt = 1 s to estimate the midpoint temperature after 10 s.

There are 11 nodes i = 0 through 10. Boundary nodes 0 and 10 stay at 0°C; all interior nodes start at 100°C. Nodes 4 and 6 are not boundaries.

alpha = 200/(7800×500) = 5.128205e-5 m²/s. Fo_step = 0.005128205, comfortably below 0.5.

At the first step, T_1 = 100 + Fo_step(0-200+100) = 99.487179°C and T_9 is equal by symmetry. Nodes 2 through 8 remain 100°C at this discrete step.

Repeating ten simultaneous updates gives approximately:

x (m)	T after 10 s (°C)
0	0
0.1	95.100589
0.2	99.887931
0.3	99.998466
0.4	99.999986
0.5	99.99999983

The remaining half is symmetric. The computed midpoint is essentially 100°C on this coarse grid at this short time. Do not treat the displayed digits as physical accuracy; the thermal penetration length sqrt(alpha t) ≈ 0.0226 m is much smaller than the 0.1 m spacing. Refinement is necessary to resolve the near-end transient accurately.

Boundary conditions and validation

Specify temperature, heat flux or convection at every boundary. Verify units and energy balance, compare a simple limiting case with an analytic result, and refine the grid and time step. For this fixed-end problem the temperatures should remain between 0 and 100°C. Implicit methods can remove this explicit stability restriction but still need time-step accuracy checks.

Quick check

If Δx is halved while Δt stays fixed, Fo_step increases fourfold. To preserve Fo_step, divide Δt by four. A node index is not a time index.