The Fourier Heat Equation
The heat (or diffusion) equation is one of the three fundamental second-order linear PDEs of mathematical physics � alongside the wave equation and Laplace equation. Derived by Joseph Fourier in 1822, it describes how temperature diffuses through a material over time:
Where T(x, y, t) is temperature, t is time, and a is thermal diffusivity (m�/s). The Laplacian ?�T captures the local curvature of the temperature field � regions where temperature curves upward (surrounded by hotter neighbours) heat up, and regions where it curves downward cool down.
Origin: Fourier's Law of Conduction
The heat equation follows directly from two physical laws combined:
- Fourier's Law: q = -k ?T � heat flux is proportional to the negative temperature gradient (flows from hot to cold)
- Energy conservation: ?(?cT)/?t = -?�q � rate of heat change equals net flux in
Together: ?c ?T/?t = k ?�T ? ?T/?t = (k/?c) ?�T = a ?�T.
Thermal Diffusivity of Real Materials
| Material | k (W/m�K) | ? (kg/m�) | c (J/kg�K) | a (�10?6 m�/s) |
|---|---|---|---|---|
| Copper | 401 | 8960 | 385 | 116 |
| Aluminium | 237 | 2700 | 900 | 97.5 |
| Steel (carbon) | 50 | 7850 | 470 | 13.6 |
| Glass | 1.0 | 2500 | 840 | 0.48 |
| Water (20 �C) | 0.60 | 998 | 4182 | 0.143 |
| Wood (oak) | 0.17 | 700 | 1700 | 0.143 |
| Air (20 �C) | 0.026 | 1.2 | 1005 | 21.5 |
Finite-Difference Method (FDTD)
The simulator solves the heat equation numerically using the explicit Euler finite-difference scheme, discretising space into a 200 � 200 grid with spacing ?x = ?y = 1 (normalised), and advancing time in steps ?t.
Discretisation Stencil
The 5-point star (Laplacian stencil): each cell's new temperature depends on itself plus its four orthogonal neighbours. This is the spatial second derivative by central differences: ?�T/?x� � (T[i+1,j] - 2T[i,j] + T[i-1,j]) / ?x�.
Stability � CFL Condition
The explicit scheme is only conditionally stable. If the stability parameter r = a�?t/?x� exceeds 0.25 (in 2-D), small errors grow exponentially � the grid "explodes". The simulator automatically clamps r = 0.24 at all alpha-slider values, so you can safely explore the full diffusivity range.
Boundary Conditions Implemented
| Type | Mathematical form | Physical meaning | Used in preset |
|---|---|---|---|
| Dirichlet | T|_boundary = T_bc | Fixed temperature at boundary | All presets (edges / fins) |
| Neumann (zero flux) | ?T/?n|_boundary = 0 | Insulated surface, no heat leaves | Default edge treatment |
| Interior source | T[i,j] = T_source | Constant-temperature heat source/sink within domain | Heat Sink (CPU die) |
| Zone diffusivity | a_local � a_global | Low-conductivity insulation layer | Insulated Wall |
Preset Scenarios Explained
🌡️ Rod Heating
The classic 1-D Fourier experiment: a long metal rod held at 600 �C at the left end and 20 �C at the right end. At steady state (?T/?t = 0), the temperature profile becomes linear � T(x) = T_hot - (T_hot - T_cold) � x/L. The simulator shows the transient approach to this linear profile, with the heat front advancing from left to right.
🔥 Plate Cooling
A square plate initially at 800 �C cools through its four edges (held at 20 �C). The solution is a 2-D Fourier series. Symmetry means the centre stays hottest longest. The characteristic "bull's-eye" cooling pattern is directly visible in the fire colour map.
💨 Heat Sink
Models a CPU die (constant 200 �C central source) with aluminium cooling fins along the top edge. Fins are represented as thin Dirichlet columns held at 25 �C. Real heat-sink fins work identically: the fin efficiency ? = tanh(mL)/(mL), where m = v(hP/kA) (h = convective coefficient, P = fin perimeter).
🧱 Insulated Wall
Two zones: hot interior (600 �C, left) and cold exterior (20 �C, right), separated by an insulation layer with a reduced to 5% of the bulk value. The temperature jump across the thin insulation layer demonstrates thermal resistance R = L/(kA): thicker or lower-k insulation gives a steeper T-drop.
🔄 Phase Change
A heat front propagates from the left boundary (900 �C) into a cold domain (20 �C), analogous to the advance of a solidification or melting front. In a real material undergoing phase change, a Stefan condition ?L dS/dt = k ?T/?x|_front must be added (where L is latent heat and S is the front position). This preset shows the pure diffusion component.
⚠️ Thermal Shock
A ceramic disc at 900 �C is suddenly quenched in cold water (20 �C, pinned at all edges). The large, rapid thermal gradient creates biaxial thermal stress s_thermal = E�a_thermal�?T/(1-?), which can exceed material fracture toughness. This is why ceramics, glass, and graphite crack when thermally shocked. The simulation shows the steep radial gradient propagating inward.
Analytical Solutions
Steady-State Rod (1-D, Dirichlet BCs)
Cooling Plate (2-D, Fourier Series)
Semi-Infinite Solid (Step-Change BC at x=0)
This is the exact solution for the Phase Change preset's early-time behaviour: the thermal penetration depth d � 2v(at), growing as the square root of time.
Curriculum Connections
| Topic | Qualification | Concepts Covered |
|---|---|---|
| Heat conduction and convection | GCSE / A-Level Physics | Fourier's law, thermal conductivity, temperature gradient |
| Thermal physics | A-Level Physics / IB HL | Specific heat, thermal diffusivity, boundary conditions |
| PDEs and finite differences | A-Level Further Maths / Uni Maths | 2-D heat equation, CFL stability, FDTD discretisation |
| Numerical methods | Engineering & Computing degrees | Explicit vs implicit schemes, stability analysis, von Neumann |
| Materials science | A-Level / Undergraduate | Thermal conductivity, diffusivity, thermal shock, phase change |
| Building physics | Architecture / Engineering | U-values, thermal bridges, insulation R-values |