next up previous
Next: About this document ...

Last Time
Periodic or Finite Domains: Separation of Variables

Short Wavelengths Disappear Rapidly

Laplace Transform Method

Anisotropic Diffusion $ \rightarrow$ Isotropic Case

3.21 Spring 2001: Lecture 12

Solution Strategies

A number of different methods of finding a closed-form solution of the diffusion equation have been briefly described. Each one has a particular type of domain, boundary condition, or initial condition to which it is particularly suited.

Green's function solutions
Generally applicable to infinite domains, zero-flux conditions at infinity.

Separation of Variables
Finite or infinitely periodic domains with factorable boundary and initial conditions.

Laplace Transforms
Time-dependent boundary conditions on a semi-infinite (with zero flux at $ \infty$) or finite domains.

Method of Images
Extension or extraction of a solution from a known solution on different domain with a judicious choice of image source and sink initial conditions for a particular set of boundary conditions.

However, many problems will not fall into these categories--and even if they do the solution may be either too difficult or impossible to obtain. It may be said with some confidence that most problems that you as professional materials scientists will face may be described thus.

The following five step program is designed to aid you in your quest.1

  1. Draw a picture
  2. Integrate by parts
  3. Attempt to solve the parts using an appropriate method
  4. Go to the Library and copy an existing solution
  5. Use brute force

Brute Force: Numerical Methods

A closed-form solution is always preferable. From a closed-form solution, one can distinguish important physical parameters from less important ones; one can analyze limiting cases and deconstruct a solution to visualize its physical implications. However, closed-form solutions are sometimes time-consuming or impossible.

Numerical methods can be useful as a last (brute force) resort. Furthermore, numerical methods can also be a useful part of finding a closed-form solution--the visualization of a numerical solution will often provide clues of useful approximations or solution strategies.

At any rate, the tools and techniques of numerical methods have become part of the requisite apparatus of a materials scientist. A few methods of numerically solving the diffusion equation are introduced below. These are not necessarily the best methods--or the fastest (which can be an entirely different thing altogether). However, they are sufficiently general as to serve as a starting point for many numerical methods.

A typical first approach is to discretize space so that the solution at any one time can be represented by values at a set of points $ c(\vec{x}_i , t)$ and the values that may be interpolated between them.

In the finite difference approach, the values on the grid can be used and a discrete version of a spatial derivative is obtained in terms of differences at the grid points.

For instance, in two dimensions, a discrete version of the Laplacian on a square grid of size $ h$, $ \nabla^2 = \frac{\partial^2}{\partial x^2} +
\frac{\partial^2 }{\partial y^2}$, may derived directly from the definition of a derivative:

\begin{displaymath}\begin{array}{lcccc} & & & +c(x,y+h) & \\  \nabla^2 c(x,y) & ...
...h , y) & -4 c(x,y) & +c(x+h,y)\\  & & & +c(x,y-h) & \end{array}\end{displaymath} (12-1)

In the finite element approach, a fitting function such as a polynomial is obtained for the regions between the points and the coefficients of the fitting function are determined by the values at specified points (nodes).

In each case, the spatial derivative can be used to iterate the solution in time.

In the forward differencing method, the definition of the time derivative is used to infer the solution at some time $ \Delta \tau$ later, e.g. for the diffusion equation:

$\displaystyle c(x,y,t + \Delta \tau) = c(x,y,t) + D (\nabla^2_{\mbox{discrete}} c(x,y,t)) \Delta \tau$ (12-2)

However, this method becomes numerically unstable for time-steps larger than $ \Delta \tau > h^2/D$.

Another method, called the Implicit or Crank-Nicholson method, uses part of the previous solution and part of an implicit future solution to extend the stability condition. This increases the numerical stability at the expense of having to invert a set of linear equations to obtain the `future' solution.

Figure 12-1: Illustration of the forward difference and implicit time-iteration methods. The implicit method uses values from the inferred solution to stabilize the numerical time-step.

More details may be found in the textbook and its cited references.

Diffusion with Moving Interfaces

The methods for solving the diffusion equation above were presented for cases of fixed boundary conditions. However, there many examples of kinetic processes in materials where boundaries (e.g. interfaces, phase boundaries) move in response or because of diffusion. Below, methods to treat such problems will be shown to be straightforward extensions of the diffusion equation-the additional physics is a conservation principle relating the velocity of the moving interface the rate at which a conserved quantity is consumed per unit area of the interface. While exact solutions are difficult to obtain, a few general results and approximations can be obtained and applied to materials processes.

The analysis of the moving interface problem originates with Stefan who was developing a model for the rate of melting of the polar ice-caps and icebergs. This problem remains as one of the biggest alloy solidification problems. Heat must be conducted from the oceans to the melting interface to to provide the latent heat of melting and salt must be supplied as well since the equilibrium concentrations salt in the liquid and solid differ.

Interface Motion due to Heat Absorption at the Interface

To simplify the analysis of the problem, consider the heat-flux problem independently; specifically, consider freezing a liquid-solid mixture by extraction of heat:

Figure 12-2: Schematic illustration of a temperature distribution resulting in the freezing and motion of the liquid/solid interface. The velocity of the interface will depend on the difference in enthalpy density in each phase.

Assuming density, $ \rho$, is same in each phase, and equating the volume swept out with heat required for the phase change:

\begin{displaymath}\begin{split}v \Delta t A \rho \Delta h_{L \rightarrow S} = \...
...h_{L \rightarrow S} = (J_{in} - J_{out}) A \Delta t \end{split}\end{displaymath} (12-3)

where $ h$ is the enthalpy per unit volume, therefore

$\displaystyle \ensuremath{\frac{d {X(t)}}{d {t}}}\rho \Delta h_{L \rightarrow S...
...rtial{x}}}\right)\ensuremath{\left.\mbox{\rule{0pt}{16pt}}\right\vert}_{x=X(t)}$ (12-4)

Equation 12-4 is known as the ``Stefan Condition,'' $ X(t)$ is the position of the (assumed planar) interface.

It is probably wise to check for wayward minus signs. Consider the usual case, $ \Delta h_{L \rightarrow S}/T_m = \Delta s_{L
\rightarrow S} < 0$, and suppose the thermal diffusivity in the solid phase is zero (i.e. all heat is absorbed by the interface and supplied by the liquid reservoir); does the velocity of the interface have the expected sign?

Therefore the thermal diffusion equations become:

\begin{displaymath}\begin{split}\ensuremath{\frac{\partial{T_S}}{\partial{t}}}= ...
...rtial{x}^2}}\mbox{\hspace{0.5in}} X(t) < x < \infty \end{split}\end{displaymath} (12-5)

$\displaystyle J_s(x=0,t) = J_o \;\; T_S(x=X(t),t) = T_m \;\; T_L(x=X(t),t) = T_m \;\; J_s(x=\infty,t) = 0$ (12-6)

with the new unknown function, the interface position $ X(t)$, to be determined by the subsidiary Stefan condition:

$\displaystyle \ensuremath{\frac{d {X(t)}}{d {t}}}\rho \Delta h_{L \rightarrow S...
...rtial{x}}}\right)\ensuremath{\left.\mbox{\rule{0pt}{16pt}}\right\vert}_{x=X(t)}$ (12-7)

Mass Diffusion in an Alloy

The Stefan condition relates the velocity of the interface to the ``jump'' in the density of an extensive quantity. For the case of heat above, that quantity was the enthalpy density. Next, the diffusion of chemical species will be coupled to the jump in alloy composition (amount/volume) at a moving interface--an analogous Stefan condition results.

Consider a diffusion couple between two alloys at different compositions for a system with multiple phases in equilibrium at a given temperature.

Figure 12-3: Schematic illustration of diffusion in an alloy with more than one equilibrium phase at a given temperature

The mass balance at the moving interface is related to the phase diagram:

Figure 12-4: Illustration of the composition difference at an interface in local equilibrium.

$\displaystyle dN_B = (c_{\beta}^{eq} - c_{\alpha}^{eq}) A v \Delta t$ (12-8)

This must be balance by the amount going in:

$\displaystyle dN^{in}_B = -D_\alpha \ensuremath{\frac{\partial{C_\alpha}}{\partial{x}}}A \Delta t$ (12-9)

minus the amount going out

$\displaystyle dN^{out}_B = -D_\beta \ensuremath{\frac{\partial{C_\beta}}{\partial{x}}}A \Delta t$ (12-10)

Therefore, the Stefan condition is:

$\displaystyle (c_{\beta}^{eq} - c_{\alpha}^{eq}) \ensuremath{\frac{d {X(t)}}{d ...
...rtial{x}}}\right)\ensuremath{\left.\mbox{\rule{0pt}{16pt}}\right\vert}_{x=X(t)}$ (12-11)

Simple Stefan Example

A limiting case for the mass diffusion case is developed below; the result that the interface grows as $ \sqrt{t}$ is derived. This result, as shown in the textbook, is a general one for the Stefan problem with uniform diffusivity in each phase. Therefore, this result is applicable to materials processes where material must diffuse through a growing phase towards an interface where is can react and form new material--such as oxidation of a surface.

The coupled diffusion equations are:

...),t) = c_{\beta}^{eq}\;\; c_{\beta}(x=\infty,t) = 1 \end{split}\end{displaymath} (12-12)

With the simplifying assumptions that $ \tilde{D}_\beta \gg \tilde{D}_\alpha$ and a steady-state profile applies in $ \alpha$-phase, the concentration profiles become:

\begin{displaymath}\begin{split}c_{\alpha}(x,t) = c^{eq}_{\alpha}\frac{x}{X(t)} ...
...eq}_{\beta} \mbox{\hspace{0.5in}} X(t) < x < \infty \end{split}\end{displaymath} (12-13)

Incorporating this limiting case into the Stefan condition and integrating,

$\displaystyle X^2(t) - X^2(t=0) = -\frac{2 \tilde{D}_\alpha c_{\alpha}^{eq}}{(c_{\beta}^{eq} - c_{\alpha}^{eq})}$ (12-14)

next up previous
Next: About this document ...
W. Craig Carter 2001-03-07