Heat Transfer Numerics: Finite Differences and Finite Elements

1. The Heat Equation and Why We Discretize

Fourier's law states that heat flux is proportional to the negative temperature gradient, q = −k∇T. Combined with conservation of energy, it gives the heat diffusion equation, a parabolic PDE:

Analytic solutions exist for infinite rods, semi-infinite slabs, or simple Fourier-series expansions on rectangles with trivial boundary conditions. A real object — a heat sink, a CPU die, a window pane — has irregular geometry and mixed boundary conditions, so we replace the continuous field T(x,y,t) with values on a discrete grid and replace derivatives with finite differences (or, for irregular shapes, with finite element basis functions).

2. Finite Difference Method: Explicit Scheme

Discretize space into a grid with spacing Δx and time into steps Δt. The second spatial derivative becomes a three-point stencil, and the time derivative becomes a forward difference — the explicit (FTCS) scheme:

Every new value depends only on old values — no linear system to solve, trivial to vectorize, and exactly what the site's own Heat Equation 2D simulation does on an 80×80 grid every animation frame. The catch is stability.

3. Stability: The CFL / Von Neumann Condition

Explicit FTCS is only conditionally stable. Von Neumann stability analysis (expanding the error as a Fourier mode and checking its amplification factor stays ≤ 1) gives a hard limit on the time step:

4. Implicit and Crank-Nicolson Schemes

The implicit (BTCS) scheme evaluates the spatial stencil at the new time level, which requires solving a linear system every step but removes the stability limit entirely:

BTCS is unconditionally stable for any Δt but only first-order accurate in time; Crank-Nicolson is second-order accurate and still unconditionally stable, at the cost of solving a tridiagonal (1D) or banded (2D, via ADI — alternating direction implicit) system each step. For 2D grids, ADI splits one full step into two half-steps, each solving a cheap tridiagonal system along one axis, which keeps the cost close to explicit while removing the CFL limit.

5. Finite Element Method: The Basic Idea

FDM needs a regular grid; the finite element method handles arbitrary geometry — an L-shaped bracket, a gear tooth, a car body — by tiling the domain with triangles (2D) or tetrahedra (3D) and representing temperature as a piecewise-linear function over the mesh, weighted by nodal values.

Once assembled, M and K are sparse matrices and the time derivative is still handled with an explicit, implicit, or Crank-Nicolson scheme exactly as in §2–4 — FEM changes how space is discretized, not how time is stepped. The payoff is that mesh density can vary locally (fine near sharp corners, coarse in open regions), something a uniform FDM grid cannot do without heavy modification.

6. Boundary Conditions in Practice

7. JavaScript Implementation

// Explicit FTCS solver for the 2D heat equation with Dirichlet edges
function stepHeat2D(T, next, nx, ny, alpha, dt, dx) {
  const r = alpha * dt / (dx * dx);
  if (r > 0.24) throw new Error(`Unstable: r=${r.toFixed(3)} exceeds 0.25 (2D CFL limit)`);

  for (let j = 1; j < ny - 1; j++) {
    for (let i = 1; i < nx - 1; i++) {
      const idx = j * nx + i;
      const lap = T[idx + 1] + T[idx - 1] + T[idx + nx] + T[idx - nx] - 4 * T[idx];
      next[idx] = T[idx] + r * lap;
    }
  }
  // Dirichlet edges: hold boundary rows/cols at their existing value
  for (let i = 0; i < nx; i++) { next[i] = T[i]; next[(ny-1)*nx+i] = T[(ny-1)*nx+i]; }
  for (let j = 0; j < ny; j++) { next[j*nx] = T[j*nx]; next[j*nx+nx-1] = T[j*nx+nx-1]; }
  return next;
}

// Thomas algorithm — solves a tridiagonal system in O(N), used for 1D implicit/Crank-Nicolson
function thomasSolve(a, b, c, d) {
  const n = d.length;
  const cp = new Float64Array(n), dp = new Float64Array(n);
  cp[0] = c[0] / b[0]; dp[0] = d[0] / b[0];
  for (let i = 1; i < n; i++) {
    const m = b[i] - a[i] * cp[i-1];
    cp[i] = c[i] / m;
    dp[i] = (d[i] - a[i] * dp[i-1]) / m;
  }
  const x = new Float64Array(n);
  x[n-1] = dp[n-1];
  for (let i = n-2; i >= 0; i--) x[i] = dp[i] - cp[i] * x[i+1];
  return x;
}

8. Applications and Limits

• Electronics cooling: heat-sink fin design is dominated by 3D FEM solves with Robin boundary conditions modelling convection to moving air.
• Building physics: whole-building thermal models use coarse implicit FDM grids stepped in hour-long increments — stability, not accuracy, is usually the binding constraint at that timescale.
• Casting and welding: phase-change (latent heat) adds a nonlinear source term, usually handled with the "enthalpy method" so solidification fronts don't need explicit tracking.
• GPU acceleration: explicit FDM is embarrassingly parallel and maps directly onto compute shaders — exactly why the browser demo above can run interactively on a laptop GPU.
• When implicit isn't worth it: for real-time visual simulations where a small time step is already needed for visual smoothness, explicit FTCS with a modest safety margin is usually simpler and fast enough.