Heat equation: A slab cooling

Your work doesn't need to be as detailed as this.

Source code for the images: code.wl
Return to the home page.

Problem

\gdef \TI {T_\mathrm{I}} \gdef \TE {T_\mathrm{E}}

Consider a uniform slab of thickness 2L, thermal conductivity k, specific heat capacity c, and density \rho.

Suppose the slab is initially at temperature \TI throughout, and thereafter loses heat to the environment at (lower) temperature \TE via Newtonian cooling through its two faces, with a heat transfer coefficient of h.

What happens?

Diagram for slab cooling example

Quantities

Take the slab to be the region -L \le x \le L. By symmetry, the temperature profile will be an even function of x, so we need only consider the right half 0 \le x \le L by using a zero-slope (or insulating) boundary condition at x = 0.

We have the following:

Independent variables
x	position (0 \le x \le L)
t	time (t \ge 0)
Dependent variable
T	temperature
Constants
L	semi-thickness
k	thermal conductivity
c	specific heat capacity
\rho	density
\TI	initial temperature
\TE	environment temperature
h	heat transfer coefficient

Defining equations

Partial differential equation (PDE)

Heat equation in T = T (x, t), throughout the slab, for all time:

\rho c \frac{\pd T}{\pd t} = k \frac{\pd^2 T}{{\pd x}^2}

Boundary conditions (BCs)

Symmetry along the central plane x = 0:

\eval{\frac{\pd T}{\pd x}}_{x = 0} = 0

Newtonian cooling along the right-hand face x = L (flux through the surface is proportional to the difference between the surface temperature and the environment temperature):

\eval{-k \frac{\pd T}{\pd x}}_{x = L} = \eval{h (T - \TE)}_{x = L}

Both sides of the equation have dimensions of power per area. The heat transfer coefficient h has dimensions of power per area per temperature.

Initial condition (IC)

Temperature is \TI throughout the slab, initially:

\eval{T}_{t = 0} = \TI

Scaling

\gdef \unscaled #1 {\colg{#1}} \gdef \scaled #1 {\colv{#1'}}

Temperature

As in the rod problem we are given two temperatures, in this case the initial temperature \TI and the environment temperature \TE. We put

\unscaled{T} = \TE + (\TI - \TE) \scaled{T}

so that the scaled initial temperature is \scaled{T} = 1 and the scaled environment temperature is \scaled{T} = 0.

Position

\unscaled{x} = L \scaled{x}

Time

As in the rod problem the time scale isn't immediately obvious, so we put

\unscaled{t} = \tau \scaled{t},

with the time scale \tau yet to be determined, i.e. free.

Dimensionless groups

\gdef \group #1 {\colr{\squarebr{#1}}}

After moving from unscaled to scaled variables, the PDE, boundary conditions, and initial condition become the following:

\begin{gathered} \frac{\pd\scaled{T}}{\pd\scaled{t}} = \group{\frac{k \tau}{\rho c L^2}} \frac{\pd^2 \scaled{T}}{{\pd\scaled{x}}^2} \\[\tallspace] \begin{aligned} \eval{\frac{\pd\scaled{T}}{\pd\scaled{x}}}_{\scaled{x} = 0} &= 0 \\[\tallspace] \eval{-\frac{\pd\scaled{T}}{\pd\scaled{x}}}_{\scaled{x} = 1} &= \eval{\group{\frac{h L}{k}} \scaled{T}}_{\scaled{x} = 1} \\[\tallspace] \eval{\scaled{T}}_{\scaled{t} = 0} &= 1 \end{aligned} \end{gathered}

Since there are two dimensionless groups but only one free scale \tau, only one of the groups can be eliminated. In particular, \tau only appears in the PDE group, so only it can be eliminated:

\group{\frac{k \tau}{\rho c L^2}} = 1

\tau = \frac{\rho c L^2}{k}

(Quick mental check: Thicker slab implies longer time scale? Better conductor implies shorter time scale? Yes to both.)

The second dimensionless group \group{h L / k} (called the Biot number) cannot be eliminated because h, L, and k are all given, and NOT free. To save writing, let us call it \gamma for short:

\gamma = \group{\frac{h L}{k}}

Finally we drop the primes:

\begin{gathered} \frac{\pd T}{\pd t} = \frac{\pd^2 T}{{\pd x}^2} \\[\tallspace] \begin{aligned} \eval{\frac{\pd T}{\pd x}}_{x = 0} &= 0 \\[\tallspace] \eval{-\frac{\pd T}{\pd x}}_{x = 1} &= \eval{\gamma T}_{x = 1} \\[\tallspace] \eval{T}_{t = 0} &= 1 \end{aligned} \end{gathered}

Equilibrium solution

Unlike the rod example, here the PDE and the two boundary conditions are already homogeneous. There is no need to subtract out the equilibrium solution because the equilibrium solution is identically zero (corresponding to the slab completely cooling to the environment temperature). So we may proceed directly to separation of variables:

Separation of variables

\gdef \pos #1 {\colg{#1}} \gdef \time #1 {\colv{#1}} \gdef \con #1 {\colb{#1}}

After putting

T (x, t) = \pos{X (x)} \time{Y (t)}

and following the usual procedure, we obtain the form

T (x, t) = \ee ^ {\con{-\lambda^2} t} \squarebr{A \cos (\con{\lambda} x) + B \sin (\con{\lambda} x)}.

Boundary conditions

Central plane

Along the plane of symmetry x = 0 we have symmetry (zero slope),

\begin{aligned} \eval{\frac{\pd T}{\pd x}}_{x = 0} &= \ee ^ {\con{-\lambda^2} t} \eval{ \squarebr{ -\con{\lambda} A \sin (\con{\lambda} x) + \con{\lambda} B \cos (\con{\lambda} x) } }_{x = 0} \\ &= \con{\lambda} B \ee ^ {\con{-\lambda^2} t} \\ &= 0, \end{aligned}

which implies B = 0. Therefore

T = A \ee ^ {\con{-\lambda^2} t} \cos (\con{\lambda} x).

Right face

Along the right face we have the Newtonian cooling condition. The left hand side is

\begin{aligned} \eval{-\frac{\pd T}{\pd x}}_{x = 1} &= A \ee ^ {\con{-\lambda^2} t} \cdot \eval{\con{\lambda} \sin (\con{\lambda} x)}_{x = 1} \\ &= A \ee ^ {\con{-\lambda^2} t} \cdot \con{\lambda} \sin\con{\lambda} \end{aligned}

and the right hand side is

\eval{\gamma T}_{x = 1} = A \ee ^ {\con{-\lambda^2} t} \cdot \gamma \cos\con{\lambda}.

Therefore \con{\lambda} \sin\con{\lambda} = \gamma \cos\con{\lambda}, or

\tan\con{\lambda} = \frac{\gamma}{\con{\lambda}}.

Eigenvalues

\gdef \lam #1 {\lambda_{\con{#1}}}

The equation for \con{\lambda} is transcendental, so the roots can only be determined numerically.

Eigenvalues for gamma equals 2

From Sturm–Liouville theory we already know that the roots will be \lam{1} < \lam{2} < \dots to infinity, and by looking at the equation (or a plot of \tan\con{\lambda} and \gamma / \con{\lambda}) we see that each \lam{n} lies between the \con{n}th root and the \con{n}th pole of tan, i.e.

(\con{n} - 1) \pi < \lam{n} < (\con{n} - 1/2) \pi.

In Mathematica you will want to use FindRoot, and the bound above suggests (\con{n} - 3/4) \pi as a reasonable initial guess for \lam{n} (although for large \con{n} the lower bound would be better). E.g. for \gamma = 2,

Module[{gamma, nMax, initialGuess, root},
  gamma = 2;
  nMax = 5;
  Table[
    initialGuess = (n - 3/4) Pi;
    root = FindRoot[Tan[#] - gamma / # &, {initialGuess}];
    {n, root}
  , {n, nMax}
  ] // TableForm
]

gives the following:

\con{n}	\lam{n}
1	1.07687
2	3.6436
3	6.57833
4	9.62956
5	12.7223

Since FindRoot is more or less a black box, there is no guarantee that each root found corresponds to the specified \con{n}. If the initial guess is poor, it could have jumped to a root corresponding to a different value of \con{n}. Therefore, always check the result of numerical root finding with a plot to see if it makes sense.

Fourier series

Thus given any \gamma we can find the \lam{n} numerically. By linearity, we have

T (x, t) = \sum_{\con{n} = 1}^{\infty} A_{\con{n}} \ee ^ {-{\lam{n}}^2 t} \cos (\lam{n} x).

The coefficients are determined from the initial condition,

\eval{T}_{t = 0} = \sum_{\con{n} = 1}^{\infty} A_{\con{n}} \cos (\lam{n} x) = 1.

Therefore

A_{\con{n}} = \frac{ \int_0^1 1 \cdot \cos (\lam{n} x) \td x }{ \int_0^1 \cos^2 (\lam{n} x) \td x },

and we are done.

How long to reach equilibrium?

We restore the dropped primes lest we confuse ourselves. At large times the solution will be dominated by the \con{n} = 1 term (corresponding to the smallest eigenvalue \lam{1}). Setting the exponent to 4, the dimensionless equilibrium time is given by

\begin{aligned} {\lam{1}}^2 \scaled{t} &\sim 4 \\ \scaled{t} &\sim \frac{4}{{\lam{1}}^2}. \end{aligned}

The unscaled equilibrium time is therefore

\unscaled{t} \sim \frac{4\tau}{{\lam{1}}^2} = \frac{4 \rho c L^2}{{\lam{1}}^2 k}.

Note that \lam{1} is a function of \gamma = \group{h L / k}.

END

Return to the home page.