Solving the Heat Equation with Oscillatory B.C.
In this tutorial, we will address the heat equation problem under special boundary conditions. The first case is a semi-infinite domain with an oscillatory B.C., \( x \geq 0 \), and
\begin{align}
u(x=0,t) = u_0 \cos(\omega t) \tag{1}
\end{align}
where \( \omega \) is the frequency of temperature variation at the edge. We can also express (1) using a complex exponential:
\begin{align}
u(x=0,t) = u_0 e^{i\omega t} \tag{2}
\end{align}
This will be more convenient and afterwards we just need to take the real part of the solution. The B.C. at large distance \( x \to \infty \) is physically decaying, \( u \to 0 \). This problem can be approached by the method of separation of variables as well: let \( u(x,t) = X(x)T(t) \) as usual, then similar to (9) of the last tutorial, we have
\begin{align}
X(x)T'(t) &= \kappa X^{\prime\prime}(x)T(t) \tag{3}
\end{align}
but now we do the separation as follows:
\begin{align}
\frac{T'(t)}{T(t)} &= \frac{\kappa X^{\prime\prime}(x)}{X(x)} = c_1 \tag{4}
\end{align}
For the time part, we have
\begin{align}
\int \frac{T’}{T} dt &= \int c_1 dt \\
\ln |T| &= c_1 t + c_2 \\
T &= Ae^{c_1t} \tag{5}
\end{align}
where \( A = \pm e^{c_2} \). By comparing this to (2), we identify with \(A = u_0 \) and \( c_1 = i\omega \). Then the spatial part of (4) becomes
\begin{align}
\frac{\kappa X^{\prime\prime}}{X} &= i\omega \\
X^{\prime\prime} -\frac{i\omega}{\kappa}X &= 0 \tag{6}
\end{align}
Recalling the solution for second-order constant coefficient ODEs, we know that
\begin{align}
X = c_3e^{qx} + c_4e^{-qx} \tag{7}
\end{align}
where
\begin{align}
q = \sqrt{\frac{i\omega}{\kappa}} = (1+i) \sqrt{\frac{\omega}{2\kappa}} \tag{8}
\end{align}
Since the B.C. requires \( u \to 0 \) as \( x \to \infty \), it must be \(c_3 = 0\). To match (2), we simply take \(c_4 = 1\), so that the full solution will be
\begin{align}
u = XT &= u_0e^{-qx}e^{i\omega t} \\
&= u_0 e^{-\sqrt{\frac{\omega}{2\kappa}}x} e^{i(\omega t-\sqrt{\frac{\omega}{2\kappa}}x)} \tag{9}
\end{align}
or in the real-valued form:
\begin{align}
u = u_0 e^{-\sqrt{\frac{\omega}{2\kappa}}x} \cos(\omega t-\sqrt{\frac{\omega}{2\kappa}}x) \tag{10}
\end{align}
The skin depth is the depth at which the amplitude of the periodic perturbation falls to \( 1/e \) compared to the surface, and inferring from (10), it is
\begin{align}
L = \sqrt{\frac{2\kappa}{\omega}} \tag{11}
\end{align}
The higher the frequency of the temperature oscillation, the less penetrating the skin depth.
Solving the Heat Equation over an Infinite Domain
The heat equation over a one-dimensional infinite domain \( -\infty < x < \infty \) can be solved via the heat kernel:
\begin{align}
G(x,t) = \frac{1}{\sqrt{4\pi\kappa t}}\exp(-\frac{x^2}{4\kappa t}) \tag{12}
\end{align}
To show this, we simply substitute (12) into both sides of the heat equation and see if they are equal:
\begin{align}
\frac{\partial G}{\partial t} &= -\frac{1}{2\sqrt{4\pi\kappa}t^{3/2}}\exp(-\frac{x^2}{4\kappa t}) + \frac{x^2}{4\kappa t^2}\frac{1}{\sqrt{4\pi\kappa t}}\exp(-\frac{x^2}{4\kappa t}) \tag{13}
\end{align}
\begin{align}
\frac{\partial G}{\partial x} &= -\frac{x}{2\kappa t}\frac{1}{\sqrt{4\pi\kappa t}}\exp(-\frac{x^2}{4\kappa t}) \\
\kappa\frac{\partial^2 G}{\partial x^2} &= \kappa(-\frac{1}{2\kappa t}\frac{1}{\sqrt{4\pi\kappa t}}\exp(-\frac{x^2}{4\kappa t}) + \frac{x^2}{4\kappa^2 t^2} \frac{1}{\sqrt{4\pi\kappa t}}\exp(-\frac{x^2}{4\kappa t})) \\
&= -\frac{1}{2\sqrt{4\pi\kappa}t^{3/2}}\exp(-\frac{x^2}{4\kappa t}) + \frac{x^2}{4\kappa t^2}\frac{1}{\sqrt{4\pi\kappa t}}\exp(-\frac{x^2}{4\kappa t}) = \frac{\partial G}{\partial t} \tag{14}
\end{align}
Moreover, notice that if \( t = 0 \), then (12) produces a Dirac delta function. First, the integral of (12) over \(x\) at any \(t\) is \(1\):
\begin{align}
\int_{-\infty}^{\infty} G(x,t) dx &= \int_{-\infty}^{\infty} \frac{1}{\sqrt{4\pi\kappa t}}\exp(-\frac{x^2}{4\kappa t}) dx \\
&= \frac{1}{\sqrt{4\pi\kappa t}} \int_{-\infty}^{\infty} \sqrt{4\kappa t} \exp(-(\frac{x}{\sqrt{4\kappa t}})^2) d(\frac{x}{\sqrt{4\kappa t}}) \\
&= \frac{1}{\sqrt{\pi}} (\sqrt{\pi}) = 1 \tag{15}
\end{align}
where we have used the value of the Gaussian integral \(\int_{-\infty}^{\infty} e^{-u^2} du = \sqrt{\pi} \). It means that \( G(x,t) \) represents the spread of one unit of heat over space and time. At the initial time \( t = 0 \) and any location \(x \neq 0\) away from the origin, it is clear that \( G(x,t) = 0 \), and so this unit of heat must concentrate at \(x = 0\), resulting in a Dirac delta.
Now we can proceed to solve the heat equation if the initial temperature distribution over the infinite domain is specified, let’s say \( u(x,t=0) = u_0(x) \). The solution is subsequently constructed by a convolution:
\begin{align}
u(x,t) = \int_{-\infty}^{\infty} G(x-z,t) u_0(z) dz \tag{16}
\end{align}
We now show that (16) solves the heat equation and it matches the initial condition. For the first part, we substitute (16) into the heat equation accordingly:
\begin{align}
\frac{\partial u}{\partial t} &= \frac{\partial}{\partial t}\int_{-\infty}^{\infty} G(x-z,t) u_0(z) dz \\
&= \int_{-\infty}^{\infty} (\frac{\partial}{\partial t}G(x-z,t)) u_0(z) dz \\
&= \int_{-\infty}^{\infty} (\kappa\frac{\partial^2}{\partial x^2}G(x-z,t)) u_0(z) dz \\
&= \kappa\frac{\partial^2}{\partial x^2} \int_{-\infty}^{\infty} G(x-z,t) u_0(z) dz = \kappa \frac{\partial^2 u}{\partial x^2} \tag{17}
\end{align}
where we have used (14) and assume that the partial derivatives can be passed across the integral. For the second part, we simply note that (16) when \( t = 0 \) becomes
\begin{align}
u(x,t=0) &= \int_{-\infty}^{\infty} G(x-z,0) u_0(z) dz \\
&= \int_{-\infty}^{\infty} \delta(x-z) u_0(z) dz = u_0(x) \tag{18}
\end{align}
since the heat kernel reduces to a Dirac delta as shown above, and we recall the defining property of the Dirac delta function. Physically, the temperature profile at a later time \(t\) is found as the sum of evolving heat kernels (response) from the initial heat distribution (signal) over all the places.
A Special Case of Semi-infinite Domain
Although the above solution is derived for an infinite domain, it can be applied to a special case with a semi-infinite domain through a smart modification. The setting of the problem is
\begin{align}
\left\{\begin{aligned}
u(x=0,t) &= A \\
u(x,t=0) &= 0
\end{aligned}\right. \tag{19}
\end{align}
where \(A\) is a fixed constant and \( x \geq 0 \). First, for convenience, we make a change of variable, \( v = A-u \), so that the B.C. at \(x = 0\) is now homogeneous:
\begin{align}
\left\{\begin{aligned}
v(x=0,t) &= 0 \\
v(x,t=0) &= A
\end{aligned}\right. \tag{20}
\end{align}
The trick is to mirror the problem to an infinite one by using odd symmetry: \( v(x=0,t) \) is still always \(0\), but now \( -\infty < x < \infty \) and
\begin{align}
v(x,t=0) &= \begin{cases}
A & x > 0 \\
-A & x < 0
\end{cases} \tag{21}
\end{align}
Then we can apply (16) in our case:
\begin{align}
v(x,t) &= \int_{-\infty}^{\infty} G(x-z,t) v_0(z) dz \\
&= A(\int_{0}^{\infty} G(x-z,t) dz -\int_{-\infty}^{0} G(x-z,t) dz) \tag{22}
\end{align}
For the second integral, we make another change of variable \( z’ = 2x -z \) so that it becomes
\begin{align}
-\int_{-\infty}^{0} G(x-z,t) dz &= \int_{\infty}^{2x} G(z’-x,t) dz’ \\
&= -\int_{2x}^{\infty} G(x-z’,t) dz’ \tag{23}
\end{align}
by noting that \( G(x-z) = G(z-x) \). The integrals in (22) then cancel out and what remains is
\begin{align}
v(x,t) &= A\int_{0}^{2x} G(x-z,t) dz \\
&= A\int_{0}^{2x} \frac{1}{\sqrt{4\pi\kappa t}}\exp(-\frac{(x-z)^2}{4\kappa t}) dz \\
&= 2A\int_{0}^{x} \frac{1}{\sqrt{4\pi\kappa t}}\exp(-\frac{z^2}{4\kappa t}) dz \\
&= 2A\int_{0}^{z=x} \frac{1}{\sqrt{\pi}}\exp(-\frac{z^2}{4\kappa t}) d(\frac{z}{\sqrt{4\kappa t}}) \\
&= A \text{erf}(\frac{x}{\sqrt{4\kappa t}}) \tag{24}
\end{align}
where we have again exploited symmetry in the third line and \( \text{erf} \) is the error function. Then, the final answer will be \( u = A(1-\text{erf}(\frac{x}{\sqrt{4\kappa t}})) \).
Exercise
Revisit the heat equation problem with the oscillatory B.C., and consider a soil layer with \( k = 2.5 \text{ W} \text{ m}^{-1} \text{ K}^{-1} \), \( c = 1000 \text{ J} \text{ kg}^{-1} \text{ K}^{-1} \), \( \rho = 2500 \text{ kg} \text{ m}^{-3} \). Compute the thermal diffusivity \( \kappa \) and evaluate the skin depth if the period of the temperature variation is \( 1 \) day (diurnal), according to (11). How about if the oscillation period is \( 1 \) year now (annual)?
Answer
The thermal diffusivity is
\begin{align}
\kappa = \frac{k}{\rho c} = \frac{2.5}{(1000)(2500)} = 1 \times 10^{-6} \text{ m}^2 \text{ s}^{-1}
\end{align}
The angular frequency of the diurnal oscillation is
\begin{align}
\omega = \frac{2\pi}{1 \text{ day}} = 7.272 \times 10^{-5} \text{ s}^{-1}
\end{align}
By (11), the skin depth is therefore
\begin{align}
L = \sqrt{\frac{2\kappa}{\omega}} = \sqrt{\frac{2(1 \times 10^{-6})}{(7.272 \times 10^{-5})}} = 0.166 \text{ m}
\end{align}
For the annual case, we may recalculate \( \omega \), or use the proportionality relationship indicated by (11). Either way should give \( 3.17 \text{ m} \) (\(\sqrt{365} \approx 19.1 \text{ times}\)).
Leave a Reply