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.
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) |
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.
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:
Assuming density, , is same in each phase, and equating the volume swept out with heat required for the phase change:
(12-3) |
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) |
(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.
The mass balance at the moving interface is related to the phase diagram:
(12-8) |
This must be balance by the amount going in:
(12-9) |
(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,
(12-14) |