1D Transient Heat Conduction — Separation of Variables

It's a method to exactly determine the transient temperature distribution in plates, cylinders, and spheres of uniform cross-section. It separates the partial differential equation into the product of a "function of space only" and a "function of time only," obtaining an infinite series solution via eigenfunction expansion.

Exactly. However, it's conditional: "material properties are independent of temperature," "geometry is simple (plate/cylinder/sphere)," and "boundary conditions are linear." As long as these conditions are met, the obtained solution is mathematically exact. That's precisely why it's extremely important as a benchmark for verifying FEM.

Let me list three typical applications.

• Internal temperature prediction for quenching (rapid cooling) processes — After immersing steel in water, predicting to the second when the center temperature passes through the martensite transformation point. This can be predicted using the separation of variables method.
• Verification of FEM solutions (V&V) — The NAFEMS T2 benchmark problem uses the exact solution from the separation of variables method as the reference value.
• Rough estimation in early design stages — Estimating the solar response of a concrete wall or the temperature rise of an electronic board upon sudden power application by hand calculation before running FEM.

Exactly that. Relying solely on FEM without knowing the analytical solution from the separation of variables method is like trusting a mental arithmetic answer without a calculator.

The one-dimensional heat conduction equation is as follows, for the case of no internal heat generation and constant material properties:

Roughly speaking, copper is $\alpha \approx 1.1 \times 10^{-4}$ m²/s, steel is $\alpha \approx 1.2 \times 10^{-5}$ m²/s, and concrete is $\alpha \approx 7 \times 10^{-7}$ m²/s. There's a difference of more than two orders of magnitude. A block of copper heats up in an instant, but a concrete wall takes many hours, right?

We separate the temperature into the product of a "function of space only $X(x)$" and a "function of time only $\Gamma(t)$":

Exactly! That's the key to the separation of variables method. This constant $-\lambda^2$ is called the separation constant. The minus sign is added to reflect the physics that temperature decays (does not diverge) over time.

Good observation. Mathematically, it's the same form of eigenvalue problem. The difference is that waves continue to oscillate, but heat conduction decays over time due to $e^{-\alpha\lambda^2 t}$ and settles into a steady state.

Applying boundary conditions restricts the possible values of $\lambda$ to discrete eigenvalues $\lambda_n$. If we non-dimensionalize as $\zeta_n = \lambda_n L$, the eigenvalue equation is determined for each type of boundary condition.

It's a transcendental equation, so it can't be solved analytically. We find the roots numerically using Newton's method or Brent's method. Have you seen tables of $\zeta_1$ values for different Biot numbers in textbook appendices? Those are pre-calculated tables. For example, if Bi=1, then $\zeta_1 \approx 0.8603$.

This is the beautiful part of this method. Look at the term $\exp(-\zeta_n^2 \text{Fo})$. As $n$ increases, $\zeta_n$ also increases, so it decays exponentially faster. For Fo > 0.2, just the first term ($n=1$) gives accuracy within 0.1%.

Let's calculate concretely. For Type 1 BC, $\zeta_1 = \pi/2 \approx 1.571$, $\zeta_2 = 3\pi/2 \approx 4.712$. For Fo=0.2:

• First term: $\exp(-\zeta_1^2 \times 0.2) = \exp(-0.493) \approx 0.611$ — still significant
• Second term: $\exp(-\zeta_2^2 \times 0.2) = \exp(-4.44) \approx 0.012$ — already about 1%
• Third term: $\exp(-\zeta_3^2 \times 0.2) = \exp(-12.3) \approx 4.5 \times 10^{-6}$ — effectively zero

That's why the one-term approximation is sufficient. Heisler charts are precisely graphs of this one-term approximation.

Yes. For the very initial stage, around Fo < 0.05, sometimes over 100 terms are needed. In that case, using the semi-infinite body solution (error function solution) is more efficient. We'll cover that in detail in another topic.

The same method can be used for infinite cylinders and spheres. Only the coordinate system and eigenfunctions change.

It's the Bessel function of the first kind of order zero. When writing the Laplacian in cylindrical coordinates, Bessel functions appear instead of $\cos$ or $\sin$ like in plates. In Python, you can calculate it instantly with scipy.special.j0(x), so in practice, it's nothing to fear.

The procedure has 4 steps:

1. Calculate eigenvalues: Solve $\zeta_n \tan(\zeta_n) = \text{Bi}$ using Newton's method. The $n$-th root always exists in the interval $((n-1)\pi, \; (n-0.5)\pi)$, so bracketing with Brent's method is reliable.
2. Calculate expansion coefficients: Compute $C_n = 4\sin(\zeta_n) / (2\zeta_n + \sin(2\zeta_n))$ for each $n$.
3. Sum the series: For the desired $x$ and $t$ (Fo), compute $\theta = \sum C_n \cos(\zeta_n x/L) \exp(-\zeta_n^2 \text{Fo})$.
4. Truncation criterion: Terminate the series when $|C_n \exp(-\zeta_n^2 \text{Fo})| < \epsilon$ (e.g., $10^{-10}$).

In that case, pre-compute and cache the eigenvalues and coefficients to avoid recalculating them. Create a lookup table for $\zeta_n$ and $C_n$ indexed by Biot number. For moderate Bi (0.1 to 10), 10-20 eigenvalues are usually enough. Only recalculate when you change material or geometry.

There are specialized libraries like DPLeak and commercial codes built on separation of variables solutions. For academic work, scipy.optimize.brentq (root finding) and scipy.special (special functions like Bessel) are your best friends. Combining these two, you can implement the full algorithm in under 200 lines of Python.