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.

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

Draw a picture

Integrate by parts

Attempt to solve the parts using an appropriate method

Go to the Library and copy an existing solution

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
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 ,
, may derived directly
from the definition of a derivative:

(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
later, e.g. for the diffusion
equation:

(12-2)

However, this method becomes numerically unstable for time-steps
larger than
.

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, , is same in each phase, and equating the
volume swept out with heat required for the phase change:

(12-3)

where is the enthalpy per unit volume, therefore

(12-4)

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

It is probably wise to check for wayward minus signs.
Consider the usual case,
,
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:

(12-5)

(12-6)

with the new unknown function, the interface position ,
to be determined by the subsidiary Stefan condition:

(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.

(12-8)

This must be balance by the amount going in:

(12-9)

minus the amount going out

(12-10)

Therefore, the Stefan condition is:

(12-11)

Simple Stefan Example

A limiting case for the mass diffusion case is developed below; the
result that the interface grows as 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:

(12-12)

With the simplifying assumptions that
and
a steady-state profile applies in -phase, the concentration profiles
become:

(12-13)

Incorporating this limiting case into the Stefan condition and integrating,