Last Time

Steady-State Solutions to the Diffusion Equation

Scaling Solutions for the Infinite Domain

Characteristic Diffusion Times and Lengths

Solution for Step Function IC's

Fundamental Point Source Solution in Infinite 1D Domain

3.21 Spring 2002: Lecture 08

Fundamental Solutions with Different Dimensionalities

A solution for a point source initial condition with diffusion in one dimension was derived in the last lecture. The same solution would naturally apply to a line source diffusing into a plane or for a plane source diffusion into three-dimensional space.

The same method can be used to find fundamental solutions for cases where the finite source is a object of a given dimensionality embedded in an infinite space of larger dimension. This is called the ``co-dimension'' and it is the dimensionality of the infinite region minus the dimensionality of the embedded object.

All the pertinent examples are contained in the following table:

Co-dimension	Example	Symmetric Part of $\nabla^2$	Fundamental Solution
	Point on Line
1	Line on Plane	$\ensuremath{\frac{\partial^2{}}{\partial{x}^2}}$	$\frac{c(x=0,t=0) dx}{\sqrt{4 \pi D t}} e^{\left(-\frac{x^2}{4 D t}\right)}$
	Plane in 3D
2	Point in Plane	$\frac{1}{r}\ensuremath{\frac{\partial{}}{\partial{r}}}\left( r \ensuremath{\frac{\partial{}}{\partial{r}}} \right)$	$\frac{c(x=0,t=0) 2 \pi r dr}{{4 \pi D t}} e^{\left(-\frac{r^2}{4 D t}\right)}$
	Line in 3D
3	Point in 3D	$\frac{1}{r^2}\ensuremath{\frac{\partial{}}{\partial{r}}}\left( r^2 \ensuremath{\frac{\partial{}}{\partial{r}}} \right)$	$\frac{c(x=0,t=0) 4 \pi r^2 dr}{(4 \pi D t)^{3/2}} e^{\left(-\frac{r^2}{4 D t}\right)}$

For example, the fundamental solution for co-dimension 2 can be used for arbitrary initial conditions for diffusion in an infinite thin film or the solution of co-dimension 3 can be used for initial conditions in three-dimensions.

However, each case above is strictly applicable to an infinite domain because that was an underlying assumption in the derivation of the point-source solution.

Method of Images or Symmetric Extension

While the method of superposition of point-source solutions is generally always useful for infinite domains, it is possible to finesse the use of the fundamental solutions for other domains. The technique usually involves extending the computational domain by employing symmetries suggested by the boundary conditions.

A sufficiently complicated example serves to explain the use and generalization of the technique.

Figure 8-1: Example of initial and boundary conditions for a problem on the positive quadrant of the infinite plane. There is a specified source in the corner with Dirichlet conditions on the positive $y$-axis and Neumann conditions on the positive $x$-axis.

Consider the positive $x$-axis. The existence of an additional source term that is the image from a mirror-reflection across the $x$-axis would satisfy the no-flux condition along the $x$-axis. This is due to the fact that each bit of the two source would contribute in equal and opposite directions to the net flux at the boundary.

Next consider the positive $y$-axis: if an sink term were present on the negative $x$ half plane so that every bit of diffusant would flow towards the axis in the infinite problem would be removed by the sink, then the concentration would not change. It follows that if a ``negative source'' was present at the image of the mirror reflected across the $y$-axis, then the fixed concentration condition would be satisfied along the positive $y$-axis.

The method is illustrated in the figure below:

Figure 8-2: Example of using image sources and sinks to replace the boundary conditions in the preceeding figure with an infinite system. The solution for the quadrant is the sum from each of the images.

This method could be extended to a finite domain, but then one would need to sum over an infinite number of (periodic) images. Fortunately, there is a simpler method available for finite domains that has the infinite sums built in.

Periodic or Finite Domains: Separation of Variables

The Separation of Variables method is applicable for finite domains where the boundary conditions and initial conditions can be factored, i.e. separated by factors for each of the variables:

$$c(x,t) = X(x) T(t)$$ (08-1)

or

$$c(x,y,z,t) = X(x) Y(y) Z(z) T(t)$$ (08-2)

or

$$c(r,\theta,t) = R(r) \Theta(\theta) T(t)$$ (08-3)

For the separation of variables, it is assumed that a factored solution exists, when (or if) it is found, uniqueness guarantees that it is the correct one for the specified initial and boundary conditions.

For example on the finite one-dimensional domain, $0 < x < L$, suppose for Dirichlet BCs:

$$c(x=0,t) = 0 = c(x=L,t)$$ (08-4)

and ICs

$$c(x,t=0) = c_0$$ (08-5)

Assume that $c(x,t) = X(x) T(t)$ and insert into the diffusion equation:

$$\frac{1}{D T} \frac{d T}{d t} = \frac{1}{X} \frac{d^2 X}{d x^2}$$ (08-6)

Note that the left hand side depends only on $t$ and the right hand side depends only on $x$.

Figure 8-3: Suppose a function of two variables is plotted as a surface. The surface on the left is a function only of one variable and the surface on the right is a function of only the other variable. If the two functions (surfaces) are equal, the surfaces can only be constant-i.e. a function of neither variable.

A solution exists only if both sides equal the same constant. Let that constant be $-\lambda$.

Therefore, there are two ordinary differential equations to solve:

$$\frac{d^2 X}{d x^2} + \lambda X = 0$$ (08-7)

with BCs $X(0) = X(L) = 0$ and

$$\frac{d T}{d t} = - \lambda D T$$ (08-8)

which must get its integration constant from the initial conditions.

The first ODE has well-known solutions:

$$X(x) = \left\{ \begin{array}{lr} A \sin(\sqrt{\lambda} x) + B \co... ...-\lambda} x} & (\lambda < 0)\\ A'' x + B'' & (\lambda = 0) \end{array} \right.$$ (08-9)

and the second,

$$T(t) = \tau e^{-\lambda D t}$$ (08-10)

If the concentration is to remain finite as $t \rightarrow \infty$, then he only reasonable values for the separation constant for the given BCs and ICs are $\lambda > 0$.

Applying the BCs to Eq. 8-7, $B=0$, but there are an infinite number of other values that will satisfy the boundary conditions:

$$\lambda_n = \frac{n^2 \pi^2}{L^2}$$ (08-11)

where $n$ is any integer.

The general solutions (for the BCs and ICs given above) are:

$$\begin{split}X_n(x) = A_n \sin \left( \frac{n \pi x}{L} \righ... ...au_n e^{-\left( \frac{ D n^2 \pi^2 t}{L^2} \right)} \end{split}$$ (08-12)

However, to satisfy the initial conditions the general solution must be taken as a linear superposition of each possible solution--and with the weighting factor at $t=0$ chosen to satisfy the initial conditions. Therefore, seek $a_n$ such that:

$$c_o = \sum_{n=1}^{\infty} a_n \sin(n \pi \frac{x}{L} )$$ (08-13)

which is a Fourier sine series representation of $c_o$.¹

Remarks on the Fourier Method of Solution

(A)

From the solution above, a periodic solution in $-\infty < x < \infty$ can be constructed for $c(x=NL,t) = 0$ and $c(x,t=0) = \pm c_0$ for all integer $N$ (use the $-c_0$ in $-L < x < 0$, etc).

Figure 8-4: Illustration of equivalent initial and boundary conditions for a particular solution.

(B)
Note how each term in the Fourier series decays much faster than the previous term--the rate at which amplitudes get smaller is exponential with smaller wavelengths.

This dependence on wavelengths is a useful intuitive device for understanding the diffusion equation.

Figure 8-5: Illustration of how small wavelengths will decay much faster than longer ones for initial conditions given in Eq. 8-14.

For example, consider two different initial conditions on the infinite domain $-\infty < x < \infty$:

$$\begin{split}c_1 (x, t=0) = a \left[ 1 - \sin \left( \frac{2 ... ...\sin \left( \frac{4 \pi x}{\lambda} \right) \right] \end{split}$$ (08-14)

The solutions $c_1(x,t)$ and $c_2(x,t)$ are obviously,

$$\begin{split}c_1 (x, t=0) = a \left[ 1 - \sin \left( \frac{2 ... ...\right) e^{-\frac{(4 \pi)^2 t}{D \lambda^2}}\right] \end{split}$$ (08-15)

Therefore shorter wavelengths, $\lambda_{short}$, decay $e^{\left(\frac{\lambda_{long}}{\lambda_{short}}\right)^2}$ times faster than longer wavelengths $\lambda_{long}$.

Time Dependent Boundary Conditions, Semi-Infinite Domains

The final analytic method for obtaining an analytic solution to to the diffusion equation is the method of Laplace transforms. This method is applicable to semi-infinite domains and is especially effective for problems with time-dependent boundary conditions.

The method of Laplace transforms is an example of operator calculus. In this case, the operator is a transformation of time derivatives to algebraic expressions in a transformation variable. The transformed diffusion equation becomes an inhomogeneous ordinary differential equation in the spatial variable. The ordinary differential equation is solved for the transformed boundary conditions and then the transformation is reversed--usually through a table of Laplace transform pairs.

The Laplace Transform is defined as the linear operator:

$${\ensuremath{\cal L}}\{f(x,t)\} \equiv \hat{f}(x,p) = \int_0^\infty e^{-pt} f(x,t) dt$$ (08-16)

that maps $f(x)$ to another function $\hat{f}(p)$.

The essential property of the Laplace transform is its operation on time derivatives:

$${\ensuremath{\cal L}}\{\frac{\partial f}{\partial t}\} = \int_0^\infty e^{-pt} \frac{\partial f(x,t)}{\partial t} dt$$ (08-17)

which can be integrated by parts,

$${\ensuremath{\cal L}}\{\frac{\partial f}{\partial t}\} = p {\ensuremath{\cal L}}\{f\} - f(x,t=0) = p \hat{f(p)} - f(x,t=0)$$ (08-18)

so the initial condition gets integrated into the the transformed derivative.

The transform has no effect on spatial derivatives, therefore the transformed diffusion equation becomes:

$$p \hat{c}(x,p) - c(x,t=0) = D \ensuremath{\frac{\partial^2{\hat{c}(x,p)}}{\partial{x}^2}}$$ (08-19)

which is an ordinary differential equations since there are only derivatives with respect to $x$.

Two Examples of the Laplace Transform Method

Figure 8-6: Semi-infinite, one-dimensional domain with BCs $c(x=0,t) = C_0$ and $J(x=\infty ,t) = 0$ and uniform initial conditions, $c(x>0,t=0) = C_i$.

For this initial condition, the transformed diffusion equation, Eq. 8-19, becomes

$$D \ensuremath{\frac{\partial^2{\hat{c}(x,p)}}{\partial{x}^2}}- p \hat{c(x,p)} = - c(x,t=0) = -C_i$$ (08-20)

The general solution, $\hat{c}$, to this differential equation is the sum of two different kinds of solutions. The first kind of solution is the one that satisfies Eq. 8-20 for a finite value of the right-hand-side ( $C_i \neq 0$, called the particular solution). Secondly, a sum of the particular solution and any solution for which the right-hand-side is zero (called the homogeneous solution) is also a solution to Eq. 8-20.

For this problem, the homogeneous solution is:

$$\hat{c}(x,p) = a_1 e^{qx} + a_2 e^{-qx}$$ (08-21)

where

$$q \equiv \sqrt{\frac{p}{D}}$$ (08-22)

The particular solution is:

$$\hat{c}(x,p) = \frac{C_i}{p}$$ (08-23)

so the general solution is given by:

$$\hat{c}(x,p) = a_1 e^{qx} + a_2 e^{-qx} + \frac{C_i}{p}$$ (08-24)

The transformed boundary conditions become

$$\hat{c}(x=0,p) = \int_0^\infty C_0 e^{-pt} dt = \frac{C_0}{p}$$ (08-25)

and

$$\ensuremath{\frac{\partial{\hat{c}}}{\partial{x}}}(x=\infty,p) = 0$$ (08-26)

This equation immediately above implies that $a_1 = 0$ and the other determines $a_2$,

$$\hat{c}(x,p) = \frac{C_0 - C_i}{p} e^{-qx} + \frac{C_i}{p}$$ (08-27)

This equation requires transformation back to the variable $t$ through the inverse of a Laplace transform. The inverses are usually much more difficult to find than the Laplace transforms. Fortunately, Laplace transforms and their inverses are usually tabulated in math handbooks. For example,

Selected Laplace Transform Pairs
$\hat{c}(x,p) = \int_0^\infty e^{-pt} c(x,t) dt$
$q \equiv \sqrt{(p/D)}$
$\frac{1}{p}$	$1$
$\frac{1}{p^{\nu+1}} \; \; \; \nu>-1$	$\frac{t^\nu}{\Gamma(\nu+1)}$
$\frac{1}{p + \alpha}$	$e^{-\alpha t}$
$\frac{\omega}{p^2 + \omega^2}$	$\sin \omega t$
$\frac{p}{p^2 + \omega^2}$	$\cos \omega t$
$e^{-qx}$	$\frac{x e^{-\frac{x^2}{4Dt}}}{\sqrt{4 \pi D t^3}}$
$\frac{e^{-qx}}{q}$	$\sqrt{\frac{D}{\pi t}} e^{-\frac{x^2}{4Dt}}$
$\frac{e^{-qx}}{p}$	erfc$\left( \frac{x}{\sqrt{4 D t}} \right)$
$\frac{e^{-qx}}{pq}$	$\sqrt{\frac{4 D t}{\pi}} e^{-\frac{x^2}{4Dt}} - x$   erfc$\left( \frac{x}{\sqrt{4 D t}} \right)$

So, the solution can be obtained through the use of the above table:

$$c(x,t) = (C_0 - C_i) \left[ 1 - \ensuremath{\mbox{erf}}\left( \ensuremath{\frac{x}{\sqrt{4 D t}}}\right) \right] + C_i$$ (08-28)

which could have been guessed from the symmetry of the infinite domain solution.

For a second example with time-dependent solutions,

Figure 8-7: Semi-infinite, one-dimensional domain with BCs $J(x=0,t) = J_0$ and $J(x=\infty ,t) = 0$ and uniform initial conditions, $c(x>0,t=0) = C_i$. Because, material is constant flowing into the domain, the concentration at the boundary $x=0$ will be a function of time.

For this case, the transformed diffusion equation is the same as that calculated in the previous example:

$$\hat{c}(x,p) = a_1 e^{qx} + a_2 e^{-qx} + \frac{C_i}{p}$$ (08-29)

and the zero-flux boundary condition at $x=\infty$ again implies that $a_1 = 0$. The transformed constant flux boundary condition is:

$${\ensuremath{\cal L}}\{\ensuremath{\frac{\partial{c}}{\partial{x}... ...artial{x}}}\ensuremath{\left.\mbox{\rule{0pt}{16pt}}\right\vert}_{x=0} = -a_2 q$$ (08-30)

Solving for the constant $a_2$,

$$\hat{c}(x,p) = \frac{J_0}{D p q} e^{-qx} + \frac{C_i}{p}$$ (08-31)

Using the table of Laplace transform pairs,

$$c(x,t) = C_i + \frac{J_0}{D} \left[ \sqrt{\frac{4Dt}{\pi}} e^{\frac{-x^2}{4Dt}} - x \; \mbox{\em {erfc}}(\frac{x}{\sqrt{4Dt}})\right]$$ (08-32)

Therefore, the surface concentration changes as $\sqrt{t}$:

$$c(0,t) = C_i + 2 J_O \sqrt{\frac{t}{\pi D}}$$ (08-33)

where $J_O$ is the constant surface flux.

Anisotropic Diffusion Coefficients

The methods of solution that have been treated do not account for the general case of an anisotropic diffusion coefficient. However, a simple strategem can be used to reconstruct the diffusion equation from its tensor form into the scalar form that is treated above.

Recall form of the the flux relationship in a particular frame of reference:

$$\left( \begin{array}{c} J_x \\ J_y \\ J_z \end{array} \right) =... ...ac{\partial }{\partial y}\\ \frac{\partial }{\partial z} \end{array} \right) c$$ (08-34)

The diffusion equation becomes:

$$\ensuremath{\frac{\partial{c}}{\partial{t}}}= \left( \frac{\parti... ...ac{\partial }{\partial y}\\ \frac{\partial }{\partial z} \end{array} \right) c$$ (08-35)

This would be difficult to solve in the general case. However, there is a trick for the case where the anisotropic diffusivity does not depend on $c$.

First, find the rotation $a_{ij}$ which diagonalizes $\ensuremath{\underline{D}}$:

$${\ensuremath{\underline{\hat{D}}}} = \left( \begin{array}{ccc} \h... ...11} & 0 & 0\\ 0 & \hat{D}_{22} & 0\\ 0 & 0 & \hat{D}_{33} \end{array} \right)$$ (08-36)

The $\hat{D}_{ii}$ are the eigenvalues of $\ensuremath{\underline{D}}$ and the coordinate frame of reference of ${\ensuremath{\underline{\hat{D}}}}$ have axes parallel to the eigenvalues of $\ensuremath{\underline{D}}$.

The diffusion equation is much simpler in the eigenframe:

$$\frac{\partial c}{\partial t} = \hat{D}_{11} \frac{\partial^2 c}{... ...{\partial \hat{x_2^2}} + \hat{D}_{33} \frac{\partial^2 c}{\partial \hat{x_3^2}}$$ (08-37)

Define an effective diffusivity $D_{eff}$:

$$D_{eff} = (\ensuremath{\underline{D}})^{1/3} = (\ensuremath{\underline{\hat{D}}})^{1/3} = (\hat{D}_{11} \hat{D}_{22} \hat{D}_{33})^{1/3}$$ (08-38)

and rescale the lengths of the coordinate axes with:

$$\hat{x_i} = \frac{\hat{D}_{11}^{1/2}} {D_{eff}^{1/2}} \xi_i$$ (08-39)

Using the $D_{eff}$ and the rescaled lengths,

$$\frac{\partial c}{\partial t} = {D_{eff}}\nabla^2_{\vec{\xi}} \cd... ... \frac{\partial^2 c}{\partial \xi_2^2} + \frac{\partial^2 c}{\partial \xi_3^2})$$ (08-40)

which is the simple scalar isotropic diffusion equation to which the solution methods (scaling, method of images, Fourier sum, and Laplace transform) were applied.