V08 — Steady and transient heat conduction

← Back to the validation dossier

Tier A (analytic) - priority P0 Reference frozen on 2026-08-13. Nabla 0.1.0, solver licence mode enforcing, run in 48.0 s.

1. Problem

Six heat-conduction problems, each solved twice - once with linear (P1) and once with quadratic (P2) elements. Twelve models, no coils, no magnets, no current anywhere: this is the chapter where the thermal module is the only thing under test.

Sub-model	Formulation	What it is	What only it can fail on
slab_p1 / slab_p2	planar	slab, uniform generation, both faces held (BC 20)	the Cartesian kernel, the volumetric source, Dirichlet
conv_p1 / conv_p2	planar	the same slab cooled by a film (BC 22)	the Robin term, which enters the matrix and the load
annulus_p1 / annulus_p2	axisymmetric	hollow cylinder wall, both surfaces held	the r weight and the 2*pi depth override, read in watts
solid_p1 / solid_p2	axisymmetric	solid cylinder, generation + film at r = R	the axis, which carries no boundary condition
trans_p1 / trans_p2	planar, transient	slab quenched at both faces from T0	the theta scheme against a Fourier series
lumped_p1 / lumped_p2	planar, transient	body with a film and a ramped q(t) table	the time-dependent heat-source path

The parameters are the single source of truth in build.py; the whole set:

	slab / conv	annulus	solid	trans	lumped
geometry	L = 20 mm, H = 20 mm	r1 = 10, r2 = 40, H = 20 mm	R = 25, H = 20 mm	L = 25, H = 20 mm	20 x 20 mm
k [W/(m K)]	20	15	30	40	200
q [W/m^3]	2.0e6	-	1.0e6	-	2000 t (ramp)
boundary	T_s = 40 degC / h = 500, T_inf = 25 degC	T1 = 100, T2 = 20 degC	h = 400, T_inf = 20 degC	T_s = 20 degC	h = 50, T_inf = 20 degC
transient	-	-	-	T0 = 100 degC, dt = 1 s, 80 steps, CN	T0 = 20 degC, dt = 10 s, 60 steps, BE

Two model-building points that are the physics and not the plumbing:

The unmarked boundaries are the problem statement. Top and bottom of every slab, and the axis of the solid cylinder, carry no boundary condition at all. An unassigned thermal boundary is adiabatic - natural Neumann - and that is what makes these 1-D problems solved on a 2-D mesh. It is only safe to rely on because marker 1 is reserved as "no BC": Triangle runs without -B and rewrites every unmarked boundary segment to marker 1, so before that reservation an unmarked edge silently inherited the first BC in the model (the thermal solver specification section 2).

The axis of solid is deliberately bare. Zero radial flux at r = 0 is natural in the thermal weak form, and the thermal driver must not reuse the magnetic formulation's axis handling, which clamps the unknown to zero there (the thermal solver specification section 1.3). If it did, the on-axis temperature would be 0 degC instead of 56.4583 degC, and this row is what would say so.

2. Reference

Closed form throughout - Fourier conduction with constant properties, derived in the case's reference notes and evaluated by harness/analytic.py at run time. A disagreement between the closed form and the frozen numbers is reported as DRIFT and fails the case, so the reference cannot move under the chapter.

slab      T(x)  = T_s + q (L^2 - x^2)/(2k)            q''(face) = q L
conv      T(x)  = T_inf + q L/h + q (L^2 - x^2)/(2k)
annulus   T(r)  = T1 + (T2 - T1) ln(r/r1)/ln(r2/r1)   Q = 2 pi k H (T1-T2)/ln(r2/r1)
solid     T(r)  = T_inf + q R/(2h) + q (R^2 - r^2)/(4k)
trans     (T - T_s)/(T0 - T_s) = SUM (4/((2n+1)pi)) (-1)^n
                                 cos((2n+1) pi x/2L) exp(-(2n+1)^2 pi^2 Fo/4)
lumped    theta' + theta/tau = q(t)/(rho cp),  tau = rho cp V/(h A)

Sources: Carslaw & Jaeger (2nd ed., ch. 3.3 and ch. 7) for the series and the cylinder; Incropera & DeWitt (6th ed., ch. 3.5 and ch. 5.1-5.2) for generation with a convective surface and for the lumped-capacitance criterion. Neither is needed to check the numbers: every formula above is two lines of algebra and is written out in harness/analytic.py.

What is judged, and why it is a rise

Every judged profile is a temperature difference from the level the boundary condition sets - T - T_s, T - T_inf, T - T_outer. 0.5 % of the 60 degC centre reading of the held slab is mostly 0.5 % of the 40 degC the Dirichlet faces were handed as an input, and would be met by a solver that got the physics half right. 0.5 % of the 20 K rise is a real tolerance. The Celsius readings are published as their own rows - they are simply not where the accuracy claim is made, and the scalar temperature rows are judged in kelvin against the same fraction of their rise.

The energy balance

Each steady problem also carries a conservation check, and it is read across an interior surface at 80 % of the way out rather than at the boundary. There both sides are exactly known - q_x = q x in the slab, q pi r^2 H inside the rod, 2 pi k H dT/ln(r2/r1) through any cylinder of the annulus - and the recovered heat flux, which is element-constant before it is nodally averaged, is not being read one element from a boundary where that averaging is one-sided. It is the same check as "what leaves equals what was generated", and a harder one to pass by accident.

3. Nabla model

Built entirely from nabla_api by build.py: geometry as four lines per rectangle, one region, thermal properties through setThermalRegionProperty, thermal boundaries through assignBC(20) and assignBC(22), and runThermalSolver() in place of the harness's default electromagnetic pipeline. The planar models declare an axial length of 1 m so that every reported watt is a watt per metre; the axisymmetric ones declare none, because the solver replaces the depth with 2*pi in axisymmetric mode and the model's own z extent is then the height.

Sub-model	Element order	Nodes	Elements
slab_p1	P1	1555	2967
slab_p2	P2	6076	2967
conv_p1	P1	1555	2967
conv_p2	P2	6076	2967
annulus_p1	P1	1770	3400
annulus_p2	P2	6939	3400
solid_p1	P1	1513	2893
solid_p2	P2	5918	2893
trans_p1	P1	1235	2342
trans_p2	P2	4811	2342
lumped_p1	P1	783	1461
lumped_p2	P2	3026	1461

The two orders share a mesh topology - the same element count, the P2 run carrying the midside nodes - because the only thing changed between them is polyOrder, and the mesh is regenerated so Triangle is called with -o2.

4. Results

Quantity	Metric	Nabla	Reference	Error	Tolerance	Verdict
slab_rise_profile_p1_K	profile	min=1.943, max=20, rms=14.8 K	min=1.95, max=20, rms=14.8 K	L2 +0.043 %, max +0.068 %	L2 0.500 %, max 0.750 %	PASS
slab_centre_p1_degC	abs	59.9972 degC	60 degC	-0.002832	0.1	PASS
slab_face_flux_p1_W_per_m2	rel	40018.3 W/m^2	40000 W/m^2	+0.046 %	0.500 %	PASS
slab_balance_p1_W	rel	1280.25 W	1280 W	+0.020 %	0.500 %	PASS
conv_rise_profile_p1_K	profile	min=81.94, max=100, rms=93.85 K	min=81.95, max=100, rms=93.85 K	L2 +0.007 %, max +0.014 %	L2 0.500 %, max 0.750 %	PASS
conv_centre_p1_degC	abs	124.997 degC	125 degC	-0.002832	0.5	PASS
conv_surface_p1_degC	abs	105.491 degC	105.497 degC	-0.006347	0.4025	PASS
conv_balance_p1_W	rel	1280.25 W	1280 W	+0.020 %	0.500 %	PASS
annulus_rise_profile_p1_K	profile	min=2.207, max=71.95, rms=36.38 K	min=2.206, max=71.93, rms=36.37 K	L2 +0.024 %, max +0.043 %	L2 0.500 %, max 0.750 %	PASS
annulus_heat_flow_p1_W	rel	108.791 W	108.777 W	+0.013 %	0.500 %	PASS
solid_rise_profile_p1_K	profile	min=31.76, max=36.46, rms=34.87 K	min=31.76, max=36.46, rms=34.87 K	L2 +0.003 %, max +0.010 %	L2 0.500 %, max 0.750 %	PASS
solid_centre_p1_degC	abs	56.4658 degC	56.4583 degC	0.007514	0.1823	PASS
solid_surface_p1_degC	abs	51.3533 degC	51.3536 degC	-0.000327	0.1568	PASS
solid_balance_p1_W	rel	25.1309 W	25.1327 W	-0.007 %	0.500 %	PASS
trans_centre_p1_K	profile	min=3.012, max=74.89, rms=30.81 K	min=3.012, max=75.09, rms=30.8 K	L2 +0.116 %, max +0.267 %	L2 1.000 %, max 1.500 %	PASS
trans_half_p1_K	profile	min=2.13, max=57.81, rms=22.24 K	min=2.13, max=57.52, rms=22.2 K	L2 +0.233 %, max +0.505 %	L2 1.000 %, max 1.500 %	PASS
trans_startup_p1_K	report	[79.91124252435219, 79.57696223261244, 79.07817649161628, 78.17774796292977, 76.77731711153336] K	- K	-	-	REPORT
lumped_centre_p1_K	profile	min=1.556, max=75.82, rms=40.09 K	min=1.367, max=75.51, rms=39.81 K	L2 +0.861 %, max +0.509 %	L2 1.000 %, max 1.500 %	PASS
lumped_startup_p1_K	report	[0.07911904817776971, 0.23423331699050465, 0.4623414543636635, 0.7605606326405798, 1.1261218787358622] K	- K	-	-	REPORT
lumped_final_p1_K	rel	75.8209 K	75.5144 K	+0.406 %	1.000 %	PASS
slab_balance_signed_p1_W	report	[-639.9625037361353, 640.2921506019713] W	- W	-	-	REPORT
conv_balance_signed_p1_W	report	[-639.9620790993392, 640.2919500944433] W	- W	-	-	REPORT
annulus_heat_flow_signed_p1_W	report	108.791 W	- W	-	-	REPORT
solid_balance_signed_p1_W	report	25.1309 W	- W	-	-	REPORT
slab_rise_profile_p2_K	profile	min=1.948, max=20, rms=14.8 K	min=1.95, max=20, rms=14.8 K	L2 +0.008 %, max +0.012 %	L2 0.100 %, max 0.200 %	PASS
slab_centre_p2_degC	abs	59.9986 degC	60 degC	-0.001376	0.02	PASS
slab_face_flux_p2_W_per_m2	rel	40000 W/m^2	40000 W/m^2	+0.000 %	0.500 %	PASS
slab_balance_p2_W	rel	1280 W	1280 W	+0.000 %	0.500 %	PASS
conv_rise_profile_p2_K	profile	min=81.95, max=100, rms=93.85 K	min=81.95, max=100, rms=93.85 K	L2 +0.001 %, max +0.002 %	L2 0.100 %, max 0.200 %	PASS
conv_centre_p2_degC	abs	124.999 degC	125 degC	-0.001376	0.1	PASS
conv_surface_p2_degC	abs	105.495 degC	105.497 degC	-0.002009	0.0805	PASS
conv_balance_p2_W	rel	1280 W	1280 W	+0.000 %	0.500 %	PASS
annulus_rise_profile_p2_K	profile	min=2.206, max=71.94, rms=36.37 K	min=2.206, max=71.93, rms=36.37 K	L2 +0.004 %, max +0.009 %	L2 0.100 %, max 0.200 %	PASS
annulus_heat_flow_p2_W	rel	108.778 W	108.777 W	+0.001 %	0.500 %	PASS
solid_rise_profile_p2_K	profile	min=31.76, max=36.46, rms=34.87 K	min=31.76, max=36.46, rms=34.87 K	L2 +0.000 %, max +0.001 %	L2 0.100 %, max 0.200 %	PASS
solid_centre_p2_degC	abs	56.4583 degC	56.4583 degC	-6.889e-05	0.03646	PASS
solid_surface_p2_degC	abs	51.3536 degC	51.3536 degC	-6.026e-05	0.03135	PASS
solid_balance_p2_W	rel	25.1327 W	25.1327 W	-0.000 %	0.500 %	PASS
trans_centre_p2_K	profile	min=3.016, max=74.9, rms=30.83 K	min=3.012, max=75.09, rms=30.8 K	L2 +0.153 %, max +0.263 %	L2 1.000 %, max 1.500 %	PASS
trans_half_p2_K	profile	min=2.133, max=57.86, rms=22.26 K	min=2.13, max=57.52, rms=22.2 K	L2 +0.301 %, max +0.593 %	L2 1.000 %, max 1.500 %	PASS
trans_startup_p2_K	report	[79.91025509892577, 79.57441363989808, 79.07413303924466, 78.17331729319503, 76.77529168231351] K	- K	-	-	REPORT
lumped_centre_p2_K	profile	min=1.556, max=75.82, rms=40.09 K	min=1.367, max=75.51, rms=39.81 K	L2 +0.861 %, max +0.509 %	L2 1.000 %, max 1.500 %	PASS
lumped_startup_p2_K	report	[0.07911886125125989, 0.23423276618251876, 0.4623403722958699, 0.7605588609824423, 1.1261192677498144] K	- K	-	-	REPORT
lumped_final_p2_K	rel	75.8207 K	75.5144 K	+0.406 %	1.000 %	PASS
slab_balance_signed_p2_W	report	[-640.0000000000131, 640.0000000000549] W	- W	-	-	REPORT
conv_balance_signed_p2_W	report	[-640.0000000001744, 640.0000000002199] W	- W	-	-	REPORT
annulus_heat_flow_signed_p2_W	report	108.778 W	- W	-	-	REPORT
solid_balance_signed_p2_W	report	25.1327 W	- W	-	-	REPORT
transient_alpha_m2_per_s	report	1.11483e-05 m^2/s	- m^2/s	-	-	REPORT
lumped_tau_s	report	243 s	- s	-	-	REPORT
lumped_biot	report	0.00125 -	- -	-	-	REPORT

• trans_startup_p1_K: published without a verdict, by declaration
• lumped_startup_p1_K: published without a verdict, by declaration
• slab_balance_signed_p1_W: published without a verdict, by declaration
• conv_balance_signed_p1_W: published without a verdict, by declaration
• annulus_heat_flow_signed_p1_W: published without a verdict, by declaration
• solid_balance_signed_p1_W: published without a verdict, by declaration
• trans_startup_p2_K: published without a verdict, by declaration
• lumped_startup_p2_K: published without a verdict, by declaration
• slab_balance_signed_p2_W: published without a verdict, by declaration
• conv_balance_signed_p2_W: published without a verdict, by declaration
• annulus_heat_flow_signed_p2_W: published without a verdict, by declaration
• solid_balance_signed_p2_W: published without a verdict, by declaration
• transient_alpha_m2_per_s: published without a verdict, by declaration
• lumped_tau_s: published without a verdict, by declaration
• lumped_biot: published without a verdict, by declaration

Three things in that table are worth reading twice.

P2 is not "more accurate" here; on three of the four steady problems it is exact. The slab and the cooled slab have parabolic solutions and the solid cylinder's is a parabola in r; a quadratic element reproduces them to the linear solver's tolerance and nothing else. The energy balance says it plainly: 1280 W against a generated 1280 W, an error of +0.000 %. That is not a tolerance being met, it is the discretization error being zero. The annulus, whose solution is logarithmic and therefore not in the P2 space, is the one that still shows a finite error - +0.004 % L2 against +0.024 % at P1.

The axisymmetric heat flow is the load-bearing row. A profile can be right while the r weight in the kernel and the 2*pi depth override are both wrong - they do not appear in the temperature at all. Only an integrated power sees them, and annulus_heat_flow is that number: 108.778 W against the closed form's 108.777 W (+0.001 %).

The transient is held to a tolerance the scheme can meet, and says so. See the discussion below.

5. Discussion

The transient is run with Crank-Nicolson, and that is a statement about Backward Euler

The validation plan V08 asks for 1 % after the first five steps. The slowest mode of the quenched slab has lambda_1 = alpha (pi/2L)^2 = 0.044 1/s, and a first-order scheme accumulates roughly lambda^2 dt t/2 on it - about 16 % by t = 80 s at dt = 2 s. That is Backward Euler being first order, not Nabla being wrong, and holding the code to 1 % under it would have meant either 620 time steps or a tolerance chosen to be met rather than to mean something.

So the quenched slab runs Crank-Nicolson with a two-step Rannacher startup, and lands at +0.153 % L2 (P2) / +0.116 % (P1) over records 6 to 80. Backward Euler is still exercised - by the lumped sub-model, where it happens to be exact (below). The order study itself is V19's subject, and this case's numbers are one of its inputs rather than a substitute for it.

The first five records are published without a verdict, on purpose. The initial condition is discontinuous - the interior at 100 degC, the faces at 20 degC from t = 0+ - so every mode the mesh can carry is excited at once, the first two steps are Backward Euler by design to damp them, and this is where a time-integration bug would show. Showing it is more useful than judging it against a series that is itself slowest to converge exactly there.

Backward Euler is exact on the lumped ramp, which is why the ramp was chosen

The lumped body is driven by a two-point heat-source table, q(t) = 2000 t. Two things follow. The solver interpolates the table linearly in time, so a ramp is represented exactly and the reference stays a closed form rather than becoming a quadrature. And Backward Euler applied to theta' + theta/tau = a t reproduces the continuous particular solution theta_p = a tau (t - tau) term for term - substitute theta^n = A t_n + B and the discrete equations give A = a tau, B = -a tau^2, the same coefficients. The scheme's whole error therefore sits in the decaying homogeneous part, whose amplitude is 48.6 K at t = 0 and 4.1 K at the end.

What is left over is the reference's error, not the solver's. The measured final rise is 75.8209 K against the lumped ODE's 75.5144 K, +0.406 % high - and it should be high, because the probe is at the centre of a body that is not actually isothermal. With Bi = 0.00125 the centre sits a few tenths of a kelvin above the lumped mean, which is what the +0.3 K is. That is why lumped_biot is published as a row of its own: a document that quotes a lumped answer without quoting its Biot number is quoting an unbounded error. P1 and P2 agree to the last digit here, which is the other half of the same statement - the residual is not discretization.

One harness trap, recorded because it cost a run

The energy balance is integrated as mean(q_n) * length, not with a trapezoid over the arclengths the contour probe returns. The probe subdivides a path into nSamples - 1 segments and returns one value per segment tagged with that segment's starting arclength, so the reported abscissae stop one segment short of the end. A trapezoid over them integrates over 1 - 1/(n-1) of the contour: a clean -0.5 % at 200 segments, which is exactly the tolerance this case's conservation check is held to. On the first run all four balance rows came in at -0.48 to -0.51 % and one of them failed - for a reason with nothing to do with the solver. The segment mean times the geometric length is the Riemann sum those samples actually are, and after the fix the P2 balances are exact to thirteen digits.

6. Limitations

What this case does not prove:

• Nothing is coupled. The thermal module is uncoupled from the electromagnetic solver: no loss field is passed to it here, no temperature is fed back to a material. The loss-to-thermal handoff is V23's subject.
• Every material is linear and constant. No temperature-dependent k, no contact resistance, no radiation (BC type 23 is reserved and unimplemented).
• No convection is modelled, only its boundary condition. h is an input; nothing here validates a heat-transfer correlation.
• The geometries are 1-D solutions on 2-D meshes - straight or circular boundaries, smooth solutions, no re-entrant corner. That is deliberate (it is what makes exactness the expected answer, and P2's exactness meaningful), and it means this case says nothing about a singular corner.
• The transient tolerance is a statement about a scheme at a step size, not a general accuracy claim. Read it with the discussion above and with V19.
• BC type 21 (imposed heat flux) is not exercised here. Types 20 and 22 are. A flux boundary appears in the thermal solver specification's own V3 check; adding it to this register is open work.
• Periodic thermal boundaries (types 3 and 7) are not exercised.

7. Reproduce

cd validation
python -m harness.run_case V08

Needs a built solver (make solver), triangle.exe (make triangle; 32-bit on Windows), and a licence carrying the Thermal feature - the thermal module and P2 elements are both gated, and without them the case cannot run at all rather than running degraded. The twelve models are written under artifacts/models/ and can be opened in the GUI.

To regenerate the frozen expectations from the closed forms alone:

python cases/V08_heat_conduction/reference/freeze_expected.py

Frozen-reference changelog

Date	Change	Reason
2026-08-13	frozen at case creation	Every reference value is a closed form evaluated from the parameters in build.py, and every tolerance is transcribed from the validation plan V08. No number here was read off a Nabla run; reference/freeze_expected.py regenerates the file from the same two sources.