Heat Conduction in a Slab, a Cylinder and a Sphere

Prepared by: Prof. Alexander Dowling (2026)

This notebook is the PDE companion to Lecture 7. We take one transient conduction problem, write it for three geometries at once, discretize space with the method of lines, and solve the resulting differential-algebraic system with Pyomo.dae.

The physical problem and the idea of indexing the three geometries with a single shape parameter are adapted from J. C. Kantor, ND Pyomo Cookbook, notebook 5.03 “Heat Conduction in Various Geometries”. The derivation, the code and the numerical experiments below are our own.

Learning Objectives¶

• Nondimensionalize a transient conduction PDE and recognize the Fourier number as its only time scale

• Write one model that covers a slab, a cylinder and a sphere through a shape index $s$

• Apply the method of lines: discretize space first, then hand a DAE in time to Pyomo.dae

• Compare finite difference, finite volume and finite element spatial discretizations, and say where each belongs

• Handle the $r = 0$ coordinate singularity correctly

• Run the IDAES Diagnostics Toolbox on a Pyomo.dae model and read its report -- including the difference between a warning and a caution

• Recognize a dangling variable: a derivative that appears in no discretization equation, which IPOPT will happily declare optimal

The Physical Problem¶

A solid body, initially at uniform temperature $T_0$, is dropped at $t=0$ into a bath held at $T_\infty$. The surface reaches the bath temperature immediately, and we want the temperature everywhere inside as a function of time.

With constant thermal conductivity $k$, density $\rho$ and heat capacity $C_p$, and no source term, the energy balance is

Divide by $\rho C_p$ and define the thermal diffusivity $\alpha = k / (\rho C_p)$, which has units of length$^2$/time:

One Equation, Three Geometries¶

The three classical bodies -- an infinite slab of half-thickness $R$, an infinite cylinder of radius $R$, and a sphere of radius $R$ -- are all one-dimensional: temperature depends on a single space coordinate $r$ measured from the center. In each case the Laplacian of a radially symmetric field is

where the shape index $s$ is

For the slab, $s=0$ and the extra term vanishes -- Cartesian coordinates. For $s = 1$ and $s = 2$ the geometry contributes the $\tfrac{s}{r}\partial T/\partial r$ term, which is singular at the center. That singularity is the interesting part of this notebook, so hold onto it.

Nondimensionalization¶

Three scales are available: the temperature difference the body must cross, the body’s own size, and the time the heat needs to cross it. Define

Substituting, with $T = T_0 + (T_\infty - T_0) T'$, $r = R r'$ and $t = R^2 t' / \alpha$:

Every term carries the same factor $\alpha (T_\infty - T_0)/R^2$, which therefore cancels:

with

Two things are worth noticing.

1. No parameters survive. The dimensionless problem depends only on $s$. Solve it once per geometry and you have solved it for every material, every size and every pair of temperatures.

2. $t'$ is the Fourier number, $\mathrm{Fo} = \alpha t / R^2$: the ratio of the time elapsed to the time heat needs to diffuse across the body. $\mathrm{Fo} \approx 1$ is “thoroughly heated.”

The last condition is the symmetry (or no-flux) condition. It is not an extra physical assumption -- it follows from the body being symmetric about $r'=0$, so no heat can cross the center. It is, however, the boundary condition that will cause all the trouble.

From here on we drop the primes and write $T$, $r$, $t$ for the dimensionless quantities.

The Method of Lines¶

We have one equation in two independent variables. Pyomo.dae is comfortable with either, but the standard strategy for a parabolic PDE is to break the symmetry between them:

Method of lines. Discretize the spatial derivatives on a fixed mesh $0 = r_0 < r_1 < \dots < r_N = 1$. What is left is a system of ODEs (here, a DAE, because the boundary rows are algebraic) in time alone, one equation per mesh point -- the “lines” of the name.

Writing $T_i(t) \approx T(t, r_i)$, the semi-discrete system is

with $\mathcal{L}_i$ the discrete approximation to $\partial^2 T/\partial r^2 + (s/r)\,\partial T/\partial r$ at node $i$, plus algebraic rows at $i = 0$ (symmetry) and $i = N$ (the bath temperature).

Why this order?

• The time direction is an initial-value problem, and Lecture 7’s machinery -- orthogonal collocation on finite elements -- is built exactly for that. We can hand the time direction to dae.collocation and get a high-order, stiffly stable integrator for free.

• The space direction is a boundary-value problem in disguise. All of its structure -- the boundary conditions, the coordinate singularity, the conservation properties -- is decided by how we build $\mathcal{L}_i$. Doing it deliberately, rather than letting a transformation do it for us, is the whole game.

In Pyomo.dae both directions are ContinuousSets and both get a TransformationFactory, so the method of lines is a way of thinking about the model rather than a different API. But it tells you which transformation to worry about.

Three Ways to Discretize Space¶

All three families replace $\mathcal{L}$ by a sparse linear operator on nodal values. They differ in what they ask the approximation to reproduce exactly, and that choice is what makes one appropriate and another awkward for a given problem.

Finite difference (FD): match the Taylor series¶

Pick a stencil and choose coefficients so the truncation error is as high order as possible. On a uniform mesh of spacing $h$,

• Cheap and transparent. Every row is three numbers. This is what dae.finite_difference builds.

• Appropriate when the geometry is simple, the mesh is structured, and the coefficients ($k$, $\rho C_p$) are smooth.

• Awkward when you need a guarantee. FD approximates the differential equation pointwise; it makes no promise that the discrete system conserves energy, and it has nothing to say at a point where the differential equation itself is singular -- such as $r = 0$ here.

Finite volume (FV): match the integral balance¶

Partition the domain into control volumes (here, shells) and integrate the conservation law over each one. Using the self-adjoint form $\partial_t T = r^{-s}\,\partial_r ( r^{s} \partial_r T )$ and integrating $r^{s}\,dr$ over the cell $[r_{i-1/2},\,r_{i+1/2}]$,

Approximating each face gradient by a two-point difference gives a scheme in which the flux leaving one cell is identically the flux entering the next, so the sum over cells telescopes: the scheme conserves energy to machine precision on any mesh, coarse or fine.

• Appropriate when conservation matters (reactors, flow, shocks), the mesh is unstructured, or the coefficients jump at an interface -- FV puts the interface on a face, where the flux is single valued.

• Note what it does at $r = 0$: the inner face of the first cell has area $r_{1/2}^{\,s} \to 0^{s} = 0$ for $s > 0$. The symmetry condition is structural: no special case, no boundary condition to write. We will use this below.

Finite element (FE): match a weighted residual¶

Expand $T(t,r) \approx \sum_j T_j(t)\,\phi_j(r)$ in basis functions and require the residual to be orthogonal to a set of test functions. Galerkin FE integrates the weak form by parts, which lowers the smoothness the solution must have and produces symmetric, provably stable operators; orthogonal collocation, the variant used throughout Lecture 7, instead forces the residual to vanish exactly at a set of collocation points inside each element.

• Appropriate when the geometry is complicated, the solution is smooth enough to reward high order, or you want error estimates and adaptivity on a rigorous footing.

• In Pyomo.dae this is dae.collocation. We use it in time below. It can also be applied in space -- and, as we will see, it inherits exactly the same $r=0$ trap as the finite difference scheme.

The Coordinate Singularity at $r = 0$¶

For $s > 0$ the term $\dfrac{s}{r}\dfrac{\partial T}{\partial r}$ is undefined at $r = 0$. The solution is perfectly well behaved there -- by symmetry $\partial T/\partial r \to 0$ at least as fast as $r$, so L’Hopital’s rule gives $\lim_{r \to 0} \frac{s}{r}\frac{\partial T}{\partial r} = s\,\frac{\partial^2 T}{\partial r^2}$ -- but the equation as written cannot be evaluated at $r=0$. Two consequences for the discrete model:

1. Never write the PDE row at $r = 0$. In Pyomo, return pyo.Constraint.Skip. The row that takes its place is the symmetry condition $\partial T/\partial r|_{r=0} = 0$, which is what makes the center temperature well posed.

2. Skip the corner $(t, r) = (0, 1)$. The initial condition asks for $T = 0$ there and the boundary condition asks for $T = 1$. Both cannot hold. We keep the boundary condition (the bath is imposed from $t=0^+$ onward) and skip the initial condition at $r=1$. Physically this is the discontinuity that makes the surface flux infinite at $t = 0^+$; numerically it is simply one row we must not write twice.

The second point is routine bookkeeping. The first one hides a trap.

A Pyomo.dae Model¶

Here is the model, written once for all three geometries. Read the three constraint rules first -- they are a direct transcription of the boxed PDE and its three side conditions.

A Pyomo.dae Model¶

Here is the model, written once for all three geometries. Read the three constraint rules first -- they are a direct transcription of the boxed PDE and its three side conditions. Then read the block after the transformations, which supplies the one equation the spatial discretization does not write for us. The next section explains why it is needed and what happens without it.

The Trap: A Derivative That Appears in No Equation¶

Everything in build_fd up to the transformations is a transcription of the mathematics. The block after them is not, and it is the part that is easy to leave out. Here is the natural way to write the symmetry condition:

@m.Constraint(m.t)
def bc_center(m, t):
    return m.dTdr[t, 0] == 0        # DO NOT DO THIS AND NOTHING ELSE

It reads exactly like $\partial T/\partial r|_{r=0} = 0$. It is also wrong, for a reason that has nothing to do with the mathematics and everything to do with what scheme="CENTRAL" actually generates.

A central difference needs a neighbour on both sides. At an interior node it has one:

At $i = 0$ there is no $T_{-1}$, so Pyomo.dae writes no equation at all for $\texttt{dTdr}[t, 0]$. It is not set to a one-sided difference; it is not set to anything. It is a variable that appears in no discretization equation -- a dangling variable. The constraint above would then not impose symmetry on the temperature field; it would assign a value to a number that touches nothing else in the model, and the center temperature would be determined by nothing at all.

We can see the hole directly, without building a model at all. Discretize a bare DerivativeVar and ask which mesh points the generated defining equations cover.

Read the CENTRAL row: neither endpoint gets a defining equation for the first derivative, and neither gets one for the second. At $r=1$ that is harmless -- the PDE row is skipped there and the bath temperature is imposed directly. At $r=0$ it is fatal, because the PDE row is skipped there too.

The other two rows are the reason “just use a one-sided scheme everywhere” is not a fix: FORWARD and BACKWARD each close one hole and open another somewhere else. They move the problem rather than solve it. See the first table below.

So build_fd writes both halves of the symmetry condition itself: it fixes $\texttt{dTdr}[t,0]$ to zero, which is the physics, and it adds bc_center_def, which is the equation the scheme never wrote. Together with $\texttt{dTdr}[t,0] = 0$ the three-point stencil collapses to $3T_0 = 4T_1 - T_2$, the second-order statement that the profile is flat at the center.

The Reference Solution¶

Each of the three geometries has a classical separation-of-variables solution, so we can measure error rather than eyeball it. For the dimensionless problem above:

Checking the Model Before You Trust It¶

This is a simulation: the temperature field is completely determined by the data, so a correctly posed model must have zero degrees of freedom. That is the cheapest check available, and it happens before the solver is ever called.

The IDAES Diagnostics Toolbox¶

Counting degrees of freedom tells you how many variables are unaccounted for. It does not tell you which. On a model with tens of thousands of rows that difference is the whole game, so run the structural analysis rather than the arithmetic.

The IDAES Diagnostics Toolbox does exactly this. report_structural_issues() needs no solution and no initialization -- it works on the incidence pattern of the model, so you can run it the moment the model is built. Make it the first thing you do to any Pyomo.dae model, before you ever look at a number.

Zero warnings. No structural singularity, no under-constrained set, no over-constrained set, degrees of freedom consistent with the constraint count. That is what a well-posed Pyomo.dae model looks like.

The two cautions are worth understanding, because they are not defects and you will see them on every model of this shape.

• 13 variables fixed to 0 -- these are $\texttt{dTdr}[t,0]$ at the 13 time points, fixed by build_fd as statement (1) of the symmetry condition. The toolbox flags variables fixed to zero because a stray fix(0) is a common way to accidentally delete physics. Here it is the physics.

• 49 unused variables -- variables that appear in no activated constraint. All 49 are inert: 9 are $\texttt{dTdt}$ at $t = 0$, which Radau collocation never references; 26 are $\texttt{d2Tdr2}$ on the two boundary rows $r = 0$ and $r = 1$, where the PDE is skipped; and 14 are $\texttt{dTdr}$ at the surface and at $(t,r) = (0,0)$, likewise unreferenced. Pyomo.dae creates a full grid of derivative variables and only some of them are used. They influence no result.

A caution is an invitation to look, not a finding. Knowing which cautions your model should produce is what lets you notice the day it produces a different one.

Alternative Discretizations at $r = 0$, and What They Do¶

The table below documents the discretizations we did not ship, so that the failure modes are on record without a broken model being a copy-paste away.

All numbers were measured on the slab ($s = 0$) unless stated, at $\mathrm{Fo} = t_f = 0.5$, with nfe_t=30, ncp=3 Radau collocation in time, and IPOPT 3.13.2 with ma27. “Center error” is $|T(t_f, 0) - T_\text{exact}(t_f, 0)|$ with $T_\text{exact}(0.5, 0) = 0.6292$; the five values run over $h = 1/10,\, 1/20,\, 1/40,\, 1/80,\, 1/160$ and “order” is $\log_2$ of the ratio of the last two.

The Answer Does Not Depend on the Initial Guess¶

A model with unconstrained degrees of freedom and no objective has no unique answer, and the quickest way to see that is to change the initial guess and re-solve. Run the same check on the correct model: if the answer moves, something is not pinned.

Identical to every digit printed across a range of guesses that includes a physically impossible one, in every geometry and on every mesh. The remaining discrepancy against exact is truncation error, and unlike a structural defect it goes away under refinement -- which is what the next section measures.

Mesh Refinement Tells the Truth¶

A structural defect does not announce itself in a single solve. It announces itself under refinement: a consistent scheme’s error must go to zero with $h$, and an underdetermined one’s does not have to. This is the check that turns “the model ran” into evidence.

Read the order column: it is $\log_2$ of the ratio of the last two errors, so 2.0 means second order in $h$, 1.0 means first order, and $\approx 0$ would mean the error is not going anywhere. All three geometries land at 1.9--2.0, which is what the interior stencil and the center row were each built to deliver -- and it is only true because they deliver it at the same order. The first-order center row in the table above produces the same picture with slope 1, and the dangling model produces a flat line.

Three straight lines of slope 2, over four mesh halvings and more than two decades of error. That is what a consistent, correctly posed discretization looks like on a log-log plot, and it is the picture the dangling model cannot produce on the slab no matter how fine the mesh gets.

Reproducibility and archived figure data¶

The two figures printed in Lecture 7 are generated from committed numerical data, not from a second implementation of this PDE. figures/results/heat-convergence.json stores the raw mesh- refinement errors for all three geometries, the initial-guess experiment, and the expected structural-diagnostics summary. figures/results/heat-slab-profile.json stores every profile and analytical-series comparison for the slab, cylinder, and sphere; the handout uses the slab subset.

A maintainer re-renders either handout figure without Pyomo or IPOPT with:

python figures/render_from_notebook.py heat-convergence heat-slab-profile

The archived runs used Pyomo 6.10.1, IPOPT 3.13.2 (IDAES build with ma27), nfe_t=30, and three Radau collocation points for the refinement experiment; the profile data use nfe_r=80, nfe_t=50, and three Radau points. Exact iteration counts and timings are intentionally not used as scientific evidence: they depend on the operating system, solver build, linear solver, and CPU. The structural claims do not: the degree-of-freedom count, expected diagnostics cautions, initial- guess invariance, and second-order convergence are the checks to reproduce.

The Correct Model, Geometry by Geometry¶

The builder is the same in all three cases. What changes is why each piece is correct.

Slab, $s = 0$¶

The PDE is $\partial_t T = \partial_r^2 T$ -- no singular term at all. We could legitimately write the PDE row at $r=0$ and drop the symmetry condition; the slab is only “half” of a symmetric plate of thickness $2R$, and $r=0$ is its midplane. We do not, for two reasons. scheme="CENTRAL" supplies no second-derivative equation at $r=0$ either, so that row would need a one-sided stencil of its own; and keeping the same skip-plus-symmetry structure for all three geometries means one code path, tested three ways. So: skip the PDE at $r=0$, impose $\partial_r T|_0 = 0$ with its defining equation.

Cylinder, $s = 1$¶

Now the skip at $r = 0$ is mandatory: $\frac{1}{r}\partial_r T$ cannot be evaluated there. The symmetry condition is doing real work -- it is the only equation that determines the axis temperature, and without a defining equation for $\texttt{dTdr}[t,0]$ the axis is unconstrained. Every interior row $i = 1, \dots, N-1$ has $r_i = i h > 0$, so the $1/r$ coefficient is finite everywhere it is used; it is large ($1/h$) at the first interior node, which is why the discrete operator is stiff near the axis and why an $\mathcal{O}(h^2)$ symmetry stencil is worth having.

Sphere, $s = 2$¶

Identical structure, with $2/r$. The $1/r$ weighting is stronger, so the near-center rows dominate the Jacobian even more -- which is also why the dangling model of the first table is least visibly wrong here. Nothing about the correct model changes.

In all three cases the count works out the same way, per time point: $N$ unknowns $T_0, \dots, T_{N-1}$ (the surface value $T_N$ is imposed), $N-1$ PDE rows at the interior nodes, and one center row. $N$ equations, $N$ unknowns, zero degrees of freedom -- which is what the DOF table above confirms.

The lecture handout carries the slab profile and the convergence plot above. The slab is the geometry where the dangling-variable defect is loud, so it is the one worth printing. The cylinder and sphere are on this page only — see the three-panel figure above.

At a fixed Fourier number the sphere is furthest along and the slab is furthest behind. That is the surface-to-volume ratio: per unit volume, a sphere presents $3/R$ of surface, a cylinder $2/R$, a slab $1/R$. The shape index $s$ is doing exactly that bookkeeping, and it is the reason $\mathrm{Fo}$ alone is not enough to say “how hot is it inside” without also saying which shape you meant.

The Finite Volume Alternative¶

Everything above was finite difference. It is worth building the finite volume version of the same problem, because it makes a structural point that no amount of care with stencils can:

In the finite volume formulation the symmetry condition at $r = 0$ is not a boundary condition at all. It is a consequence of the geometry, and there is nothing to forget.

Place $N$ shells with faces at $r_{i+1/2} = (i+1)h$, $h = 1/N$, and unknowns at the cell centers. The balance on shell $i$ is the discrete form of the integral law from earlier,

For the innermost shell, $r_{-1/2} = 0$ and therefore $A_{-1/2} = 0^{s} = 0$ for $s > 0$: the inner face has no area, so no heat can cross it. For the slab, $s = 0$ makes $A_{-1/2} = 1$, and we set the inner flux to zero by symmetry -- one line, and the only geometry-specific line in the model.

At the outer face of the last shell the temperature is known, and the half-cell distance from the center of that shell to the surface is $h/2$, so the flux is $A_{N-1/2}\,(1 - T_{N-1}) / (h/2)$.

Note also that there is no DerivativeVar in space at all. The method of lines is explicit here: the only ContinuousSet is time.

Zero degrees of freedom in every geometry, second-order convergence, and no symmetry constraint anywhere in the model. The Diagnostics Toolbox says the same thing structurally.

Zero warnings again, and this time a single caution: the 8 unused variables are $\texttt{dTdt}$ at $t = 0$, one per cell, which Radau collocation never references. There is no “fixed to 0” caution, because there is nothing to fix -- and that is the argument for finite volume in one line. Not that it is more accurate (on this uniform mesh it is not dramatically so), but that the condition we spent the first half of this notebook getting right cannot be gotten wrong: the inner face of the first shell has zero area, so no heat can cross it, structurally and on every mesh.

The price is that a finite volume model is yours to write. Pyomo.dae will not build it for you, you do not get a DerivativeVar to inspect, and the cell-center unknowns do not sit on the boundary, so reporting a surface value takes an extra half-cell extrapolation. Which is the right trade depends on the problem, and that is the point of knowing all three families.

The mistakes move with the method. Finite volume cannot lose the symmetry condition, but it has its own short list -- and, symmetrically with the first table, its worst failure is the one the slab cannot see.

Alternative Finite Volume Discretizations, and What They Do¶

Same convention as the first table: these were measured, not reasoned about, and none of them ships as a runnable option above. Conditions are $\mathrm{Fo} = t_f = 0.5$, nfe_t=30, ncp=3 Radau, IPOPT 3.13.2 with ma27, and $N = 10,\, 20,\, 40,\, 80$ cells. The error is taken at the first cell center, $r = h/2$, against the series solution evaluated there; “order” is $\log_2$ of the ratio of the last two.

Take-Home Messages¶

1. Nondimensionalize first. Three geometries, any material, any size and any temperature pair reduce to one equation in one parameter $s$, with the Fourier number as the clock.

2. The method of lines is a decision, not a syntax. Discretizing space first is what lets you hand the time direction to Lecture 7’s collocation machinery and concentrate on the spatial operator, where the physics lives.

3. Finite difference, finite volume and finite element differ in what they reproduce exactly -- the Taylor series, the integral balance, or a weighted residual. Reach for FV when conservation or an interface matters, FE when the geometry or the accuracy demands it, FD when neither does.

4. 🔴 A Neumann boundary condition written on a DerivativeVar may constrain nothing. CENTRAL finite differences and Radau collocation both leave the derivative at the left endpoint of the domain with no defining equation. Writing m.dTdr[t, 0] == 0 and stopping there assigns a value to a dangling variable, and the center temperature is unconstrained. Write the defining equation too.

5. optimal is not correct. Every failure in the two tables above terminated cleanly, several of them with plausible-looking numbers. Three things catch this class of defect: a degree-of-freedom count -- a simulation must have zero -- report_structural_issues(), which says which variables are loose rather than how many, and mesh refinement, which turns a structural defect into a visible plateau. Do all three.

6. Run the diagnostics before you run the solver. report_structural_issues() needs no solution and no initialization; it costs seconds and it reads the incidence pattern you cannot read by eye. Learn what cautions your model should emit, so that a new one means something.

7. Test on the geometry where the bug is loud. The finite difference defect was obvious in the slab and nearly invisible in the sphere; the finite volume defect is the other way round. A model verified only on the forgiving case is not verified.

Try It Yourself¶

1. Predict, before running it, which endpoint loses its defining equation under scheme="BACKWARD" and what that breaks in build_fd. Check your answer with stencil_coverage("BACKWARD") -- and then say what one-sided row you would owe the model to repair it.

2. Replace dae.finite_difference in space with dae.collocation in a copy of build_fd and confirm, with the toolbox rather than by eye, that the derivative at $r=0$ still needs a defining equation of your own.

3. Add a Robin (convective) boundary condition at the surface, $-\partial_r T|_{r=1} = \mathrm{Bi}\,(T(t,1) - 1)$, where $\mathrm{Bi} = hR/k$ is the Biot number. Which endpoint does CENTRAL leave undefined now, and what does that imply about which one-sided stencil you owe the model? Recover the results above as $\mathrm{Bi} \to \infty$.

4. Grade the mesh toward $r=1$ so the short-time boundary layer is resolved, and redo the $\mathrm{Fo}=0.01$ comparison. ContinuousSet(initialize=[...]) accepts an explicit list of points.

5. Reproduce one row of either table for yourself in a scratch copy of the builder -- then delete it. Confirming a documented failure mode is worth an hour; leaving a runnable broken model in a notebook someone else will copy is not.