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.
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
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.2
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:
(09-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:
(09-2) |
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.
Another method is called the spectral method and is very powerful for solving the diffusion equation. They work somewhat in the same way that the Laplace transform method works--a derivative is transformed into an algebraic operation.
In particular, the Fourier spectral method is commonly used because of the existence of very fast techniques to perform the transformation operations (i.e., Fast Fourier Transforms--FFTs).
In a nutshell, the Fourier spectral method is implemented by representing the spatial part of the field as a series of trigonometric functions:
(09-3) |
The diffusion equation can be written as independent ordinary differential equations (one for each coefficient of each wavelength):
(09-4) |
Each of the ODEs is stable (increasingly more stable with large or short wavelengths). A finite difference approach can be applied to each of the coefficients:
(09-5) |
So, the prediction of the future value of each coefficient can be solved simply and implicitly:
(09-6) |
The real space time-dependent solution can be obtained by application of an FFT.
More details may be found in the textbook and its cited references.
Diffusion as a Molecular Process
The derivation of the diffusion equation followed from consideration of a continuum distribution of conserved particles. A treatment of individual molecular or atomic motions as part of a statistical ensemble of diffusing particles will prove to be not only consistent with the continuum approach, but will provide new physical concepts that enhance understanding of diffusive processes.
Below, the fundamental diffusive process--an atomic or molecular displacement--is treated as a random jump occurring with an average frequency . It is instructive to consider how this simple model of random atomic jumps corresponds to the various diffusion systems, or diffusivities discussed above.
First of all, if a jump is truly random, then the probability of a given jump giving a net displacement must be exactly the same as jump with displacement . Therefore, the expected displacement, , of a particle after one or more jumps must be identically zero. Furthermore, the expected velocity of the particle, , must also be zero. Considering the Einstein drift equation , it must be concluded that the random walk treatment corresponds to the physical experiment that , which is precisely the case that is intended for the measurement of the self-diffusivity ; i.e. a diffusion couple of two alloys of the same composition, but differing in isotopic concentration.
Secondly, if the jumps are assumed to take place on a crystalline lattice, then each random displacement will be selected from a small subset of all the possible displacements that would be available in the more general case of diffusion in a gas. Treating the general case first, as is done below, leaves the diffusion on the lattice as particular case.
Third, for a process to be completely random, there can be no restriction on the availability of a nearby empty site--the probability of an empty site cannot appear in the simple model of random walk. In the simple random walk process, two atoms may occupy the same space. Therefore, the simple random walk treatment should not correspond to substitutional diffusion on a crystal lattice. However, substitutional crystalline diffusion can be treated as a sequence of random hops on a lattice if the ``particle'' that hops is one of a infinitely dilute concentration of vacancies-a vacancy will always have an occupied neighboring site to hop into in the dilute limit.
In the treatment of the random walk, the average rate at which a walker (i.e. an atom or molecule) takes steps is an empirical frequency
(09-7) |
Diffusion as a Random Walk Process
Let be the average rate that a random walker (an atom or molecule) `take steps.' is an average frequency. be the position of a particular random walker:
(09-8) |
What is meant by a ``random jump'' is not specified, but could be represented by a number of different processes:
In either case, the expected value ; this is an indication that no particular direction is favored and that the driving force for a random walk.
However, it should be intuitively clear that the distance between two particles should increase with time. This distance is measured by the spread, or second moment, of the distribution :
This can be rewritten as a sum along the diagonal and off-diagonal terms:
Using the relationship , where is the angle between the and jump-vectors. Therefore,
(09-12) |
If the lengths of the jumps do not correlate to the jump-angles:
(09-13) |
Define the Correlation Factor, :
Therefore,
This is an important equation, it relates the spread of the expected square displacements (or the second moment of the probability distribution) per unit time with three factors that are dictated by the mechanism of diffusion: the rate of successful hops , the average square hop distance , and a factor related to the correlation between hops, . If every hop is uncorrelated and the probability of positive displacement is exactly balanced by an equal probability then .
Diffusion as a Time-Dependent Probability Distribution
In the last section, the expected position and the expected square displacement of a particle was calculated by considering the statistics of a sequence of random jumps. The concept of statistics can be coupled with that of a probability distribution. The probability of finding a particle after a time will be determined and will be related to the form of a point source solution of the diffusion equation.
Consider a sequence of random jumps on a one-dimensional lattice. If the random walker is at position after jumps (after starting at ), then the number of jumps to the right and to the left must satisfy:
(09-16) |
The number of different ways that a random jumper could land at site from the origin is given by the binomial coefficient:
Therefore the probability of getting to site after jumps is:
(09-18) |
If the probability of jumping right the probability of jumping left, then
(09-20) |
(09-21) |
Using and :
(09-22) |
(09-23) |
(09-24) |
Comparing this with the fundamental point solution in 1D:
(09-25) |
The probability distribution can be used to identify the diffusivity with a random process:
(09-26) |
Relation of the Self-Diffusivity to a Random Walk
The notion of the continuum limit for the concentration, of particles is consistent with an interpretation of the concentration related to a probability of finding a particle of a given type within a small distance of the point at a time . Treating the concentration as a probability distribution, the second moment of the distribution is related to the mean-square displacement from the average displacement:
Suppose all of the particles are located at the origin at time , then if diffusion takes place in three dimensions:
(09-28) |
(09-29) |
(09-30) |
(09-31) |
Using the point source in two dimensions, the analogous result for diffusion in one dimension is:
(09-32) |
and in general for a point source in dimensions:
(09-33) |