Derivation of the Heat Equation
We will now derive the form of the heat equation, (3) of the last tutorial. Consider a one-dimensional rod in the \(x\)-direction. The flow of heat in a medium is down-gradient and follows Fourier’s Law.
\begin{equation}
q = -kA \frac{\partial T}{\partial x} \tag{1}
\end{equation}
where \( k \) is the thermal conductivity (SI unit: \(\text{W}\cdot\text{m}^{-1}\cdot\text{K}^{-1}\)), \( A \) is the cross-sectional area of the rod, and \( \partial T/\partial x \) is the gradient of the temperature \(T\) in the \(x\)-direction. We have taken the convention that a positive flux is directed to the right. Now consider an infinitesimally small slice of the rod with a length of \( \Delta x\). The inflow of heat at the left boundary is
\begin{equation}
q(x) = -kA \frac{\partial T}{\partial x}|_{x} \tag{2}
\end{equation}
Similarly, the outflow of heat at the right boundary is
\begin{equation}
q(x+\Delta x) = -kA \frac{\partial T}{\partial x}|_{x+\Delta x} \tag{3}
\end{equation}
The difference between (2) and (3) is then the net rate of change in heat within this slice of rod:
\begin{align}
\Delta q &= q(x) -q(x+\Delta x) \\
&= -kA \frac{\partial T}{\partial x}|_{x} -(-kA \frac{\partial T}{\partial x}|_{x+\Delta x}) = kA(\frac{\partial T}{\partial x}|_{x+\Delta x} -\frac{\partial T}{\partial x}|_{x}) \tag{4}
\end{align}
The total amount of heat stored/lost during a short time period \( \Delta t \) is just \( Q = \Delta q \Delta t \). By considering the specific heat relationship:
\begin{align}
Q = mc\Delta T \tag{5}
\end{align}
where \( m = \rho A \Delta x \) is the mass, with \( \rho \) being the density, and \( c \) is the specific heat capacity (SI unit: \(\text{J}\cdot\text{kg}^{-1}\cdot\text{K}^{-1}\)), and plugging in (4), we obtain
\begin{align}
kA(\frac{\partial T}{\partial x}|_{x+\Delta x} -\frac{\partial T}{\partial x}|_{x}) \Delta t &= mc\Delta T \\
kA(\frac{\partial T}{\partial x}|_{x+\Delta x} -\frac{\partial T}{\partial x}|_{x}) &= (\rho A\Delta x) c\frac{\Delta T}{\Delta t} \\
\frac{k}{\rho c}\frac{\frac{\partial T}{\partial x}|_{x+\Delta x} -\frac{\partial T}{\partial x}|_{x}}{\Delta x}&= \frac{\Delta T}{\Delta t} \tag{6}
\end{align}
Taking the limit \( \Delta x, \Delta t \rightarrow 0 \), we arrive at the heat equation:
\begin{align}
\frac{\partial T}{\partial t} = \frac{k}{\rho c}\frac{\partial^2 T}{\partial x^2} = \kappa \frac{\partial^2 T}{\partial x^2} \tag{7}
\end{align}
where \( \kappa = k/(\rho c) \) is defined to be the thermal diffusivity (SI unit: \(\text{m}^2 \cdot\text{s}^{-1}\)).
Separation of Variables for PDE (including Heat Equation)
To solve the heat equation in a bounded domain, we can usually approach it with the method of Separation of Variables, inspired by its ODE counterpart. We assume the solution may be decomposed as the product
\begin{align}
u(x,t) = X(x)T(t) \tag{8}
\end{align}
where here \(X(x)\) is a function of the space \(x\) only, \(T(t)\) is a function of the time \(t\) only, with \( u \) in place of temperature. Substituting this into the heat equation gives
\begin{align}
\frac{\partial}{\partial t}(X(x)T(t)) &= \kappa \frac{\partial^2}{\partial x^2}(X(x)T(t)) \\
X(x)T'(t) &= \kappa X^{\prime\prime}(x)T(t) \tag{9}
\end{align}
We now separate the parts of \( x \) and \( t \) as
\begin{equation}
\frac{X^{\prime\prime}(x)}{X(x)} = \frac{T'(t)}{\kappa T(t)} \tag{10}
\end{equation}
The key observation is that since the L.H.S. depends on \(x\) only while the R.H.S. depends on \(t\) only, the remaining possibility is that they are both equal to a constant that cannot depend on either \(x\) or \(t\), let’s say \(-\lambda\) (the negative sign is customary to facilitate the subsequent steps):
\begin{equation}
\frac{X^{\prime\prime}(x)}{X(x)} = \frac{T'(t)}{\kappa T(t)} = -\lambda \tag{11}
\end{equation}
Then, we can instead solve two ODEs:
\begin{align}
\frac{X^{\prime\prime}(x)}{X(x)} &= -\lambda \\
X^{\prime\prime} + \lambda X &= 0 \tag{12} \\
\frac{T'(t)}{\kappa T(t)} &= -\lambda \\
T’ + \lambda \kappa T &= 0 \tag{13}
\end{align}
(12) is a Sturm-Liouville equation (see this tutorial) that may be solved given appropriate boundary conditions, and let’s say its (countably infinite) eigenfunctions are \( X_{\lambda_1}, X_{\lambda_2}, \ldots \), indexed by the allowed eigenvalues \(\lambda_j\). Plugging in these eigenvalues, (13) can then be easily solved as a first-order linear/separable ODE:
\begin{align}
\int \frac{T’}{T}dt &= -\int \lambda_j \kappa dt \\
\ln |T| &= -\lambda_j\kappa t + C \\
T &= Ae^{-\lambda_j\kappa t} \tag{14}
\end{align}
where \( A = \pm e^C \) is some integration constant. Hence, for any \(j\) and \(\lambda_j\), \( X_{\lambda_j}e^{-\lambda_j\kappa t} \) is a possible solution to the heat equation (7). Notice that the form of (14) confirms that \( \lambda_j > 0 \), else \(T\) will grow exponentially and become unphysical. And by the linearity of the heat equation itself, the general solution will be the superposition of this over all possible \(j\), that is
\begin{align}
u(x,t) &= A_1X_{\lambda_1}(x)e^{-\lambda_1\kappa t} + A_2X_{\lambda_2}(x)e^{-\lambda_2\kappa t} + \cdots \\
&= \sum_{j=1}^{\infty} A_jX_{\lambda_j}(x)e^{-\lambda_j\kappa t} \tag{15}
\end{align}
To determine the values of the coefficients \( A_j \), we need an initial condition that specifies the distribution of temperature at the beginning, \(u_0(x) = u(x,t=0)\), then
\begin{align}
u_0(x) &= \sum_{j=1}^{\infty} A_jX_{\lambda_j}(x)e^{-\lambda_j\kappa (0)} \\
&= \sum_{j=1}^{\infty} A_jX_{\lambda_j}(x) = A_1X_{\lambda_1}(x) + A_2X_{\lambda_2}(x) + \cdots \tag{16}
\end{align}
By the Spectral Theorem (see this tutorial), \( A_j \) is then computed as the orthogonal projections of \( u_0 \) onto the eigenfunctions \( X_{\lambda_j} \), and the full solution (15) is done.
Example
A metal rod of length \( L \) initially has a temperature distribution of \( u(x,t=0) = x \), \( 0 < x < L \). The temperature at the two ends is fixed at zero, \( u(x=0,t) = u(x=L,t) = 0 \). Derive the evolution of the temperature profile.
Assume the solution is in the form of (8), \( X(x)T(t)\). To satisfy the Dirichlet boundary conditions above, we also need \( X(0) = X(L) = 0 \). Hence the solution and the eigenvalues/eigenfunctions for (12) are
\begin{align}
X = \sin(\sqrt{\lambda} x) = \sin(\frac{n\pi}{L}x), n = 1,2,\ldots \tag{17}
\end{align}
where \( \lambda = n^2\pi^2/L^2 \), and (14) becomes
\begin{align}
T = Ae^{-\frac{n^2\pi^2\kappa}{L^2} t} \tag{18}
\end{align}
The general solution will then be, according to (15):
\begin{align}
u &= A_1\sin(\frac{\pi}{L}x)e^{-\frac{\pi^2\kappa}{L^2} t} + A_2\sin(\frac{2\pi}{L}x)e^{-\frac{4\pi^2\kappa}{L^2} t} + \cdots \\
&= \sum_{n=1}^{\infty} A_n\sin(\frac{n\pi}{L}x)e^{-\frac{n^2\pi^2\kappa}{L^2} t} \tag{19}
\end{align}
And (16) means that
\begin{align}
u_0(x) = x &= \sum_{n=1}^{\infty} A_n\sin(\frac{n\pi}{L}x) \\
&= A_1\sin(\frac{\pi}{L}x) + A_2\sin(\frac{2\pi}{L}x) + \cdots \tag{20}
\end{align}
The coefficients \( A_n \) are then evaluated as the Fourier sine series of \( u_0(x) = x \) via orthogonal projections:
\begin{align}
A_n &= \frac{\int_0^{L} x \sin(\frac{n\pi}{L}x) dx}{\int_0^{L} \sin^2(\frac{n\pi}{L}x) dx} \\
&= \frac{-\frac{L}{n\pi}\int_0^{L} x d(\cos(\frac{n\pi}{L}x))}{L/2} \\
&= -\frac{2}{n\pi} \left( [x\cos(\frac{n\pi}{L}x)]_0^{L} -\int_0^{L} \cos(\frac{n\pi}{L}x) dx \right) \\
&= -\frac{2}{n\pi} ( L\cos(n\pi) -0) \\
&= -\frac{2L(-1)^n}{n\pi} \tag{21}
\end{align}
So the full solution will be
\begin{align}
u &= -\sum_{n=1}^{\infty} \frac{2L(-1)^n}{n\pi}\sin(\frac{n\pi}{L}x)e^{-\frac{n^2\pi^2\kappa}{L^2} t} \tag{22}
\end{align}
For a higher value of \(n\), the sinusoidal oscillation is more rapid, and the decay rate of temperature is faster (known as the smoothing property).
Exercise
Solve the heat equation, now with non-homogeneous Neumann boundary conditions over \( 0 \leq x \leq L \), where \( u_x(0,t) = u_x(L,t) = \gamma \), \( \gamma \) is a fixed constant, and the initial temperature profile is \( u(x,0) = x -L/2 \). What is the physical interpretation of the Neumann boundary conditions here? Also, what happens when \( \gamma = 1 \)? Note: we need a change of variable so as to transform the B.C. into a homogeneous one.
Answer
We can let \( u = v + \gamma x \) so that \( u_x(0,t) = u_x(L,t) = \gamma \) becomes homogeneous, \( v_x(0,t) = v_x(L,t) = 0 \). The heat equation remains the same after replacing \( u \) by \( v \):
\begin{align}
v_{t} = \kappa v_{xx}
\end{align}
and the initial condition becomes \( v(x,0) = (1-\gamma)x -L/2 \). For the Neumann B.C., (12) implies that
\begin{align}
X = \cos(\frac{n\pi}{L}x), n = 0,1,2,\ldots
\end{align}
particularly, \(n=0\) stands for the constant term. \( T \) has the same form as in (18), with an extra constant term \( T=1 \) for \( n=0 \). Therefore,
\begin{align}
v &= A_0 + A_1\cos(\frac{\pi}{L}x)e^{-\frac{\pi^2\kappa}{L^2} t} + A_2\cos(\frac{2\pi}{L}x)e^{-\frac{4\pi^2\kappa}{L^2} t} + \cdots \\
&= \sum_{n=0}^{\infty} A_n\cos(\frac{n\pi}{L}x)e^{-\frac{n^2\pi^2\kappa}{L^2} t}
\end{align}
Subsequently, we compute the coefficients \( A_n \) as a Fourier cosine series just like (21):
\begin{align}
A_n &= \frac{\int_0^{L} ((1-\gamma)x -\frac{L}{2}) \cos(\frac{n\pi}{L}x) dx}{\int_0^{L} \cos^2(\frac{n\pi}{L}x) dx} \\
&= \frac{\frac{L}{n\pi}\int_0^{L} ((1-\gamma)x -\frac{L}{2}) d(\sin(\frac{n\pi}{L}x))}{L/2} \\
&= \frac{2}{n\pi} \left( [((1-\gamma)x -\frac{L}{2}) \sin(\frac{n\pi}{L}x)]_0^L -(1-\gamma)\int_0^L \sin(\frac{n\pi}{L}x) dx \right) \\
&= \frac{2}{n\pi} ((0) +(1-\gamma)[\frac{L}{n\pi}\cos(\frac{n\pi}{L}x)]_0^L) \\
&= \frac{2L(1-\gamma)((-1)^n-1)}{n^2\pi^2}
\end{align}
and
\begin{align}
A_0 &= \frac{\int_0^{L} ((1-\gamma)x -\frac{L}{2}) dx}{\int_0^{L} dx} \\
&= \frac{[(1-\gamma)\frac{x^2}{2} -\frac{L}{2}x]_0^{L}}{L} \\
&= \frac{-\gamma\frac{L^2}{2}}{L} = -\frac{\gamma L}{2}
\end{align}
Therefore, the required answer is
\begin{align}
v &= -\frac{\gamma L}{2} + \sum_{n=1}^{\infty} \frac{2L(1-\gamma)((-1)^n-1)}{n^2\pi^2} \cos(\frac{n\pi}{L}x)e^{-\frac{n^2\pi^2\kappa}{L^2} t} \\
\Rightarrow u &= \gamma(x-\frac{L}{2}) + \sum_{n=1}^{\infty} \frac{2L(1-\gamma)((-1)^n-1)}{n^2\pi^2} \cos(\frac{n\pi}{L}x)e^{-\frac{n^2\pi^2\kappa}{L^2} t}
\end{align}
The Neumann B.C. represents the flow of heat across the boundary, and here the inflow is equal to the outflow. When \( \gamma = 1 \), we simply have \( u = x-\frac{L}{2} = u_0 \), i.e. the temperature profile is in equilibrium to begin with under this case.
Leave a Reply