For what types of unsteady heat conduction problems can the separation of variables method be used?

It is applicable to one-dimensional unsteady heat conduction problems in uniform cross-section plates, infinite cylinders, and spheres. Under conditions of isotropic, temperature-independent properties and linear boundary conditions, an exact infinite series solution can be obtained.

What happens when the Fourier number (Fo) exceeds 0.2?

For Fo > 0.2, the temperature distribution can be approximated with an error of less than 0.1% using only the first term (n=1) of the eigenfunction expansion. This is called the one-term approximation and is widely used in practical engineering calculations.

How is the separation of variables solution utilized for FEM verification?

The analytical solution from the separation of variables method is adopted in benchmarks like NAFEMS T2. It serves as a reference solution for verifying FEM time integration schemes and mesh convergence. The error between the FEM result and the analytical solution is confirmed to be below a predetermined threshold.

How is the eigenvalue equation determined?

The eigenvalue equation varies depending on the type of boundary condition. For a flat plate with a convective boundary condition (Robin condition), the equation is ζn·tan(ζn)=Bi, where ζn is the n-th positive root of this transcendental equation, found numerically.

1D Transient Heat Conduction — Separation of Variables

Category: 熱解析 > 非定常熱伝導 | Integrated 2026-04-12

変数分離法による1次元非定常温度分布の時間発展 — 初期の矩形分布が固有関数の重ね合わせで指数的に減衰する様子

Theory and Physics

Overview — What is the Separation of Variables Method?

🧑‍🎓

For what kind of heat conduction problems is the separation of variables method used?

🎓

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.

🧑‍🎓

When you say exact solution, does that mean it gives the precise answer, not an approximation like FEM?

🎓

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.

🧑‍🎓

Specifically, in what situations is it used?

🎓

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.

🧑‍🎓

Huh, so it's like an "answer sheet" to check if FEM is correct, right?

🎓

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.

Governing Equation and Non-Dimensionalization

🧑‍🎓

First, please teach me the basic equation.

🎓

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

$$ \frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2} $$

Here, $\alpha = k / (\rho c_p)$ is the thermal diffusivity [m²/s], $k$ is the thermal conductivity, $\rho$ is the density, and $c_p$ is the specific heat.

🧑‍🎓

How much does the value of $\alpha$ differ between materials?

🎓

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?

To handle the analysis uniformly, we introduce the dimensionless temperature $\theta$ and dimensionless numbers:

$$ \theta = \frac{T - T_\infty}{T_i - T_\infty}, \quad \text{Fo} = \frac{\alpha t}{L^2} \;(\text{Fourier number}), \quad \text{Bi} = \frac{hL}{k} \;(\text{Biot number}) $$

Here, $T_i$ is the initial temperature, $T_\infty$ is the ambient fluid temperature, $L$ is the characteristic length (half-thickness for a plate), and $h$ is the heat transfer coefficient. $\theta=1$ corresponds to the initial state, and $\theta=0$ corresponds to the steady state (matching the ambient temperature).

Procedure for Separation of Variables

🧑‍🎓

Now we get to the core of "separation of variables." How is it actually done?

🎓

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

$$ \theta(x, t) = X(x) \cdot \Gamma(t) $$

Substituting this into the heat conduction equation gives:

$$ \frac{1}{\alpha \Gamma} \frac{d\Gamma}{dt} = \frac{1}{X} \frac{d^2 X}{dx^2} = -\lambda^2 $$

🧑‍🎓

The left side is a function of $t$ only, and the right side is a function of $x$ only, yet they are equal... Does that mean both must be constant for it to make sense?

🎓

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.

This yields two ordinary differential equations:

$$ \frac{d\Gamma}{dt} + \alpha \lambda^2 \Gamma = 0 \quad \Rightarrow \quad \Gamma(t) = e^{-\alpha \lambda^2 t} $$

$$ \frac{d^2 X}{dx^2} + \lambda^2 X = 0 \quad \Rightarrow \quad X(x) = A\cos(\lambda x) + B\sin(\lambda x) $$

🧑‍🎓

The time part is a simple exponential decay. The space part is a trigonometric function... it kind of looks like the wave equation, doesn't it?

🎓

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.

Boundary Conditions and Eigenvalue Equation

🧑‍🎓

How is the value of $\lambda$ determined?

🎓

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.

For a plate of thickness $2L$, with a symmetry condition at the center $x=0$ and various boundary conditions applied at both surfaces $x=\pm L$:

Boundary Condition	Eigenvalue Equation	Eigenfunction $X_n$
Type 1 (Constant surface temperature)	$\cos(\zeta_n) = 0$ i.e., $\zeta_n = (2n-1)\pi/2$	$\cos(\zeta_n x/L)$
Type 2 (Constant surface heat flux)	$\sin(\zeta_n) = 0$ i.e., $\zeta_n = n\pi$	$\cos(\zeta_n x/L)$
Type 3 (Convection: $-k \partial T/\partial x = h(T-T_\infty)$)	$\zeta_n \tan(\zeta_n) = \text{Bi}$	$\cos(\zeta_n x/L)$

🧑‍🎓

$\zeta_n \tan(\zeta_n) = \text{Bi}$ for Type 3... can that be solved? With tan involved, hand calculation seems impossible...

🎓

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$.

Using the orthogonality of the eigenfunctions, the expansion coefficients $C_n$ are uniquely determined from the initial conditions. For a plate with a uniform initial temperature $T_i$:

$$ C_n = \frac{4 \sin(\zeta_n)}{2\zeta_n + \sin(2\zeta_n)} $$

The final solution is expressed as the following infinite series:

$$ \boxed{\theta(x, t) = \sum_{n=1}^{\infty} C_n \cos\!\left(\zeta_n \frac{x}{L}\right) \exp(-\zeta_n^2 \, \text{Fo})} $$

Fourier Number and One-Term Approximation

🧑‍🎓

An infinite series... how many terms do we need to sum?

🎓

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%.

🧑‍🎓

What, just one term? Do the higher-order terms vanish that quickly?

🎓

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.

The equation for the center temperature $\theta_0$ ($x=0$) with the one-term approximation applied:

$$ \theta_0 = C_1 \exp(-\zeta_1^2 \, \text{Fo}) \quad (\text{Fo} > 0.2) $$

🧑‍🎓

Conversely, when Fo is small, meaning the initial stage of heating/cooling, many terms are needed, right?

🎓

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.

Extension to Cylinders and Spheres

🧑‍🎓

Can it be used for shapes other than plates?

🎓

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

Shape	Coordinate	Eigenfunction	Eigenvalue Equation (Type 3 BC)
Plate	Cartesian coordinate $x$	$\cos(\zeta_n x/L)$	$\zeta_n \tan(\zeta_n) = \text{Bi}$
Infinite Cylinder	Cylindrical coordinate $r$	$J_0(\zeta_n r/r_0)$	$\zeta_n J_1(\zeta_n)/J_0(\zeta_n) = \text{Bi}$
Sphere	Spherical coordinate $r$	$\sin(\zeta_n r/r_0)/(\zeta_n r/r_0)$	$1 - \zeta_n \cot(\zeta_n) = \text{Bi}$

🧑‍🎓

What is $J_0$? I've never seen that symbol before.

🎓

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.

Physical Meaning of Each Term

Thermal storage term $\rho c_p \partial T/\partial t$: Rate of thermal energy storage per unit volume. An iron frying pan is hard to heat and cool ($\rho c_p$ is large), but an aluminum pot heats and cools easily. The larger this product, the slower the temperature change.
Heat conduction term $k \partial^2 T/\partial x^2$: Heat transport based on Fourier's law. Heat moves proportionally to the curvature (concavity/convexity) of the temperature distribution. The handle of a metal spoon gets hot when placed in a hot pot because $k$ is large.
Thermal diffusivity $\alpha = k/(\rho c_p)$: Represents "how quickly temperature equalizes." The larger it is, the faster temperature diffuses. There's a difference of 3 orders of magnitude between copper ($1.1 \times 10^{-4}$) and wood ($1.3 \times 10^{-7}$).
Fourier number $\text{Fo} = \alpha t / L^2$: Dimensionless time. The square of the ratio of the distance heat has diffused to the object size. Fo=1 roughly corresponds to the state where "heat has spread throughout the object."
Biot number $\text{Bi} = hL/k$: Ratio of surface convection resistance to internal conduction resistance. If Bi < 0.1, the interior is nearly uniform temperature (lumped capacitance method can be used). If Bi > 100, the surface temperature immediately becomes $T_\infty$ (Type 1 BC approximation).

Assumptions and Applicability Limits

Material properties are independent of temperature: $k$, $\rho$, $c_p$ are constant. For large temperature differences (several hundred degrees or more), temperature dependence cannot be ignored, and the separation of variables method is not applicable. FEM is needed.
Isotropic material: Thermal conductivity does not depend on direction. For composite materials like CFRP, other methods considering anisotropy are necessary.
No internal heat generation: If there is Joule heating or chemical reaction heat, it may be possible to handle by superposition with a steady-state solution.
Linear boundary conditions: Cannot be directly applied to boundary conditions involving radiation (proportional to $T^4$). Linearization approximation is needed.
1D problem: For 2D or higher problems, it can be extended as the product of solutions in each direction (product solution) for Cartesian coordinates.

Dimensional Analysis and Unit Systems

Variable	SI Unit	Typical Values / Notes
Temperature $T$	K or °C	Temperature difference calculations are the same for K or °C. Absolute temperature is mandatory for radiation calculations.
Thermal conductivity $k$	W/(m·K)	Copper: 400, Steel: 50, Concrete: 1.4, Air: 0.026
Thermal diffusivity $\alpha$	m²/s	Copper: $1.1\times10^{-4}$, Steel: $1.2\times10^{-5}$, Water: $1.4\times10^{-7}$
Heat transfer coefficient $h$	W/(m²·K)	Natural convection: 5–25, Forced convection: 25–250, Water quenching: 1,000–10,000
Fourier number Fo	Dimensionless	One-term approximation applicable for Fo > 0.2. The most frequently used criterion in practice.
Biot number Bi	Dimensionless	Lumped capacitance method if Bi < 0.1. The separation of variables method comes into play for 0.1 < Bi < 100.

Coffee Break Trivia Corner

Fourier's Passion and Napoleon

Joseph Fourier (1768-1830), who established the separation of variables method, accompanied Napoleon on his Egyptian expedition. It is said that experiencing the intense desert heat and extreme nighttime cold obsessed him with the question, "Why does temperature change like this?" In his 1822 masterpiece "Théorie analytique de la chaleur" (The Analytical Theory of Heat), he presented the separation of variables method and Fourier series, but contemporaries like Lagrange and Laplace did not believe the claim that "any function can be expressed as a sum of trigonometric functions." Fourier was ultimately proven correct, and this debate became the starting point for modern functional analysis.

Numerical Solution and Implementation

Procedure for Numerical Calculation of the Analytical Solution

🧑‍🎓

To actually calculate the separation of variables solution on a computer, what procedure should I follow?

🎓

The procedure has 4 steps:

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.
Calculate expansion coefficients: Compute $C_n = 4\sin(\zeta_n) / (2\zeta_n + \sin(2\zeta_n))$ for each $n$.
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})$.
Truncation criterion: Terminate the series when $|C_n \exp(-\zeta_n^2 \text{Fo})| < \epsilon$ (e.g., $10^{-10}$).

🧑‍🎓