Heat Diffusion in an Iron Sphere — A FEniCSx Simulation

The Problem

An iron sphere with a hot core sitting in air. The core starts at 800°C and the surrounding material begins at ambient temperature. How does heat spread outward? How fast does the core cool? How does the surface temperature respond?

This is a classical heat diffusion problem — but solving it correctly in 3D, on a curved domain, with physically realistic boundary conditions, requires more care than it first appears.

The Governing PDE

The temperature field $T(\mathbf{x}, t)$ inside the sphere satisfies the **heat equation**:

where $\Omega$ is the spherical domain, $\nabla^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2}$ is the Laplacian, and $\alpha = k / (\rho c_p)$ is the **thermal diffusivity**. For iron:

• Thermal conductivity: $k = 80$ W/(m·K)
• Density: $\rho = 7874 \, \text{kg/m}^3$
• Specific heat: $c_p = 449$ J/(kg·K)
• Thermal diffusivity: $\alpha \approx 2.26 \times 10^{-5} \, \text{m}^2/\text{s}$

Iron is an excellent thermal conductor. This high diffusivity means heat propagates rapidly — a 1 m sphere equilibrates internally within a few hours.

Boundary and Initial Conditions

Two physically distinct boundaries are present.

**Initial condition** — the core (radius $R_\text{core} = 0.3$ m) starts at 800°C; the rest of the sphere starts at 20°C:

**Outer surface BC — Robin condition** — the sphere loses heat to air by natural convection. This is a **Robin (mixed) boundary condition**, not a Dirichlet condition. Fixing the surface at 20°C would mean the sphere is actively cooled by an infinite heat sink — physically wrong. The correct condition is:

where $h_\text{conv} = 10$ W/(m²·K) is the convective heat transfer coefficient for natural convection in air, $T_\infty = 20°\text{C}$ is ambient temperature, and $\mathbf{n}$ is the outward normal. This condition allows the surface temperature to rise freely as heat arrives from the interior.

Finite Element Discretisation

The domain $\Omega$ (a sphere of radius $R = 1$ m) is meshed using **tetrahedral elements** generated with Gmsh. The mesh contains approximately 20,000 tetrahedra. The right panel of the animation below shows a cross-sectional slice of the element mesh in green ($\texttt{\#00ff00}$), giving a direct view of how the domain is partitioned.

Each tetrahedron is a **Lagrange $p = 1$ element** — piecewise linear — meaning the temperature is approximated as a linear function inside each cell. The solution is continuous across element boundaries but not differentiable. Using $p = 1$ keeps the system sparse and fast to solve; higher-order elements (e.g. $p = 2$) would improve accuracy at the cost of more degrees of freedom.

The **weak (variational) form** of equation (1) with the Robin condition (3) reads:

for all test functions $v$ in the finite element space. The Robin term appears naturally as a surface integral — no special treatment is needed once the mesh surface is tagged.

Time is discretised with **backward Euler** (implicit), unconditionally stable, at time steps of $\Delta t = 30 \, \text{s}$ over a 3-hour simulation window. Each time step requires solving one sparse linear system, handled by a direct LU factorisation.

The FEniCSx Implementation

The simulation uses **DOLFINx**, the Python interface to FEniCSx — a modern finite element library designed for scalable PDE solving on unstructured meshes.

from dolfinx import fem, mesh
from dolfinx.io import gmshio
from dolfinx.fem.petsc import assemble_matrix, assemble_vector
from petsc4py import PETSc
import ufl, gmsh
import numpy as np
from mpi4py import MPI

# Physical parameters (iron)
k, rho, cp = 80.0, 7874.0, 449.0
alpha = k / (rho * cp)       # thermal diffusivity ~2.26e-5 m^2/s
h_conv, T_amb = 10.0, 20.0

# Build sphere mesh with Gmsh
gmsh.initialize()
gmsh.model.add("sphere")
gmsh.option.setNumber("General.Verbosity", 0)
ov = gmsh.model.occ.addSphere(0, 0, 0, 1.0)
gmsh.model.occ.synchronize()
surfs = gmsh.model.getBoundary([(3, ov)], oriented=False)
gmsh.model.addPhysicalGroup(2, [s[1] for s in surfs], tag=2)
gmsh.model.addPhysicalGroup(3, [ov], tag=3)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.10)
gmsh.model.mesh.generate(3)
msh, _, facet_tags = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)
gmsh.finalize()

# Function space (Lagrange p=1)
V = fem.functionspace(msh, ("Lagrange", 1))
coords = V.tabulate_dof_coordinates()
r_dofs = np.sqrt((coords**2).sum(axis=1))

# Initial condition: 800 degrees inside core (r < 0.3), 20 outside
T_n = fem.Function(V)
T_n.x.array[:] = np.where(r_dofs < 0.3, 800.0, T_amb)
T_n.x.scatter_forward()

# Variational form with Robin BC on outer surface
ds_out = ufl.Measure("ds", domain=msh, subdomain_data=facet_tags, subdomain_id=2)
T_, v_ = ufl.TrialFunction(V), ufl.TestFunction(V)
dt_c = fem.Constant(msh, 30.0)          # 30-second time step
al_c = fem.Constant(msh, alpha)
hc   = fem.Constant(msh, alpha * h_conv / k)
Ta   = fem.Constant(msh, T_amb)

a = fem.form(
    T_*v_*ufl.dx
    + dt_c*al_c*ufl.inner(ufl.grad(T_), ufl.grad(v_))*ufl.dx
    + dt_c*hc*T_*v_*ds_out
)
L = fem.form(T_n*v_*ufl.dx + dt_c*hc*Ta*v_*ds_out)

# Assemble and factorise once (matrix is constant)
A = assemble_matrix(a, bcs=[]); A.assemble()
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A)
ksp.setType(PETSc.KSP.Type.PREONLY)
ksp.getPC().setType(PETSc.PC.Type.LU)
ksp.setUp()

T_h = fem.Function(V)
T_h.x.array[:] = T_n.x.array

# Time loop: 360 steps x 30s = 3 hours
for step in range(360):
    b = assemble_vector(L)
    b.ghostUpdate(addv=PETSc.InsertMode.ADD,
                  mode=PETSc.ScatterMode.REVERSE)
    ksp.solve(b, T_h.x.petsc_vec)
    T_h.x.scatter_forward()
    T_n.x.array[:] = T_h.x.array

A key design decision is that the system matrix $A$ is assembled **once** and factorised once — because the geometry, time step, and material parameters do not change. Each time step then reduces to a single triangular solve, which is extremely fast. The Robin boundary condition enters the bilinear form as a surface integral over the tagged outer boundary, appearing naturally in the weak formulation without any special treatment.

What the Simulation Shows

The animation captures the full 3-hour window at 30-second intervals.

In the first 20 minutes, the heat front is entirely inside the sphere. The outer surface remains at 20°C while the core drops from 800°C to below 300°C — iron conducts heat rapidly. By 30 minutes, the front reaches the surface and the convective loss begins. The surface temperature climbs through 22°C, 25°C, 30°C. By 90 minutes the gradient across the sphere has largely flattened, with the core near 60°C and the surface near 36°C. At 3 hours, the sphere sits close to thermal equilibrium at roughly 40°C, slowly surrendering energy to the surrounding air.

The diffusion time scale $\tau = R^2 / \alpha \approx 12 \, \text{h}$ tells us that full equilibration with the environment takes much longer than the simulation window — and the data confirms it.

Takeaway

Heat diffusion in a bounded domain combines two physics: fast internal equilibration driven by conduction, and slow external cooling driven by convection. Iron’s high diffusivity means the first process dominates early and finishes quickly. The Robin condition — often substituted by the simpler Dirichlet condition in textbook examples — is essential for physical realism. Fixing the surface temperature removes the convective resistance and produces a sphere that cools far too quickly and completely.

Need a reliable numerical implementation for your problem? Get in touch.