Diffusion Phenomena: Cases and Studies: Second Edition

Index

1

The Diffusion Equation

We begin with a model for diffusion: the isotropic one-dimensional random walk.¹ – ⁴ It is so simple that the basic physical processes cannot elude us. It also has a continuum limit, the diffusion equation, whose solutions are our main concern here and some of whose properties we then examine. Conversely, this model forms the basis for numerical methods of solution. We then discuss the diffusion equation’s form in higher dimensions and other physical instances where that equation offers a realistic description. This chapter ends with a brief account of the origin of conservation principles and constitutive relations that pervade transport phenomena.

1.1. The Isotropic One-Dimensional Random Walk

Consider points on a line, as shown in Fig. 1.1, choose an arbitrary origin, and label these points through integers i = 0, ±1, + 2, Attribute particles to each point, henceforth called sites, namely, Ni is the number of particles at site i. Let us assume that each of these particles can jump to adjacent sites with a frequency Γ that does not depend on the site i. ΓNi particles jump, per unit time, from site i to each of its nearest-neighbors.†

Fig. 1.1. Particle numbers Ni at sites i, and Ji+1/2, which is a minimum at the midpoint.

Particles that jump from site to site cause the numbers Ni to change in the course of time. To obtain a rate equation, assume that, sitting on the ith site, an observer evaluates all possible transitions to and away from that site

Collecting terms, we get a rate equation for the particle distribution Ni(t)

This equation involves three sites: the observer’s station i and his two nearest neighbors. Just as a second-order differential equation can be decomposed into two first-order equations, Eq. (1.1) is decomposed in the following way:

where

and

The quantities Ji±1/2 represent the rate of exchange, the net flux, between adjacent sites. The subscript notation i , admittedly somewhat obscure (the N’s are defined only at integer sites), will later become clear. For the moment it simply indicates that the particle fluxes (1.3) depend on two adjacent sites, in addition to time. These fluxes vanish when Ni = Ni + 1 for all i; in other words, the flow of matter vanishes when mass is evenly distributed among sites. This is nothing but a necessary condition for equilibrium.

Equation (1.1) gives the rate of change in terms of the nearest-neighbor distribution around site i, whereas Eq. (1.2) expresses that change as a difference of fluxes. In a sense, this last expression is the more fundamental, for it is a statement of a conservation law. Imagine a cell centered at site i, whose boundaries are somewhere in the intervals [i – 1, i] and [i, i + 1]. Then, Eq. (1.2) merely expresses the change in the cell’s content caused by the difference of the stuff in Ji – 1/2 and the stuff out Ji + 1/2. In other words, this equation is nothing but a statement of (local) mass conservation, and it would still be true even if the fluxes had forms different from (1.3), for example polynomial functions of the differences Ni – Ni ± 1, or functions involving next nearest-neighbor jumps. Let us pause here for a few remarks as follows:

Remarks

a. Equation (1.1) is linear in the unknown functions Ni(t). Structurally, it is a linear difference-differential equation, second-order in the discrete variable i and first-order in the continuous variable t. Given an initial configuration [the initial condition Ni(0)] and boundary conditions (on i), this initial-value problem can be solved analytically. We will not develop this here. However, Eq. (1.1) does provide a numerical scheme for its continuum analog, as Sect. 1.4 will show.

b. The particular simplicity of our random walk rests on the constancy of the transition frequency Γ. Recall that Γ is assumed independent of the site index i and of time. In addition, Γ does not depend on the actual particle distribution in the neighborhood of an ith site. Succinctly, this random walk model is homogeneous, an assumption that must often be relaxed, for example in vacancy diffusion: a jump can occur only if adjacent sites are vacant. As pointed out in the first footnote, the model is also isotropic Γ represents either frequency if Γ is the total jump freqency. Later, in Chap. 3, we shall see that external forces introduce a measure of anisotropy.

c. Equations (1.1)–(1.3) are topological in the sense that neither the distance between points nor the shape of the lattice plays any role whatsoever. Only the order of points matters: One must always be able to tell which of two adjacent points is the first. For example, these equations are valid for points distributed along a space curve — a helix, say — but points of self-intersection, as in a figure eight, require special attention. In brief, these equations are valid for sites embedded in spaces that are topologically equivalent to the real line.

d. The reader will note the similarity between Eqs. (1.1)–(1.3), applied here to mass transport, and rate equations that describe nuclear decay, chemical kinetics, and stochastic processes. In fact, similar equations describe the time evolution of any property over an ordered list of states i. If Pi(t) is that property, and Γi,k is the transition frequency from state i to k, then

is a generalization of Eq. (1.1), often called a master equation. Its probabilistic interpretation is interesting: If Pi is the α priori occupation probability of state i and Γi,k is the conditional probability of the i → k transition (given that i is occupied), then sums such as ΣίΓί,kPi represent the total probability of transition to a given final state k, regardless of the initial state i. Both sums of Eq. (1.4) can be interpreted accordingly.

EXERCISE 1.1. Steady-state solutions of the rate equations: As will be emphasized in i

1.2. Continuum Equations

In many instances we may be unconcerned with the discrete nature of matter. It is then possible to find a continuum analog of Eqs. (1.1)–(1.3). The essential step requires that the sites i be points of the real line, thereby introducing a distance, or metric, between points. Assume, then, that the sites are equidistant, a being the lattice (or jump) distance. Therefore, the ith site has the coordinate xi = ia.

Now introduce any (x, t) that interpolates the previous functions Ni(t), namely, that satisfies the condition

at sites x = xi, but that is quite arbitrary elsewhere. This is shown schematically in Fig. 1.1. Omitting the time variable for the moment, we expand the interpolating function in a Taylor series around the central site xi. For adjacent sites we get†

which, with the condition (1.5), becomes

Introducing this equation into Eq. (1.1), we get

If we agree to neglect terms of fourth order in the jump distance a, and if we assume the truncated equation’s validity at all points x (and not merely at sites), then Eq. (1.7) reduces to a single . It replaces the denumerable (possibly infinite) coupled system of difference-differential equations (1.1).N1

It is somewhat more delicate to analyze the fluxes, for it is a priori not at all evident how a continuous analog J(x, t) should be defined.

Again introducing Eqs. (1.6) into Eqs. (1.3), we find that

In other words, both discrete fluxes converge to the same value

to first order in a, and the second order terms are equal and opposite. It follows that = 0 + O(a²): We would have no continuous analog of the discrete conservation law (1.2). Therefore, let us agree to define a continuous flux function

. In other words, a flux function can be calculated for all x (x, twill ever pass exactly through the discrete fluxes (1.3), no matter where these are placed. The gap shown in Fig. 1.1 can be minimized, however, if we measure the difference between Ji + 1/2, say, and expression (1.9) evaluated at arbitrary points in the interval [i, i + 1]. Specifically, we define an arbitrary intermediate point xθ = (i + θ)a, θ 1, and expand Ni and Ni+1 around xθ to form the difference Ji (xθ, t). The result is given in the exercise that follows.

EXERCISE 1.2. To justify the notation Ji±1/2: Show that the above difference in fluxes is optimal when θ . What else can you say about the series in θ? Then develop the same arguments for the flux Ji–1/2 in the interval [i – 1, i], and show that it has the same formal expansion around the appropriate intermediate point as does Ji + 1/2. Finally, form the difference Ji + 1/2 – Ji–1/2, expand it around the current lattice site xi, and show that it is equal to

plus fourth-order terms in a.

It follows from this exercise that the discrete conservation law (1.2) has a continuum analog

We have now reached our goal, for, if C /a is the average concentration per cell, and if we introduce the definition

for the diffusion coefficient (or diffusivity for short), then Eqs. (1.9), (1.10), and (1.7) become

To summarize briefly, we obtain the continuum equations (1.12)–(1.14) from the discrete equations (1.1)–(1.3) by: (i) embedding the sites i in the real line; (ii) interpolating the discrete distribution Ni(t(x, t; (iv) truncating the resulting equations and requiring their validity at all points; (v) defining the concentration C ; (vi) defining a diffusivity D Γa² which remains constant as Γ → ∞ and a → 0.

EXERCISE 1.3. Next-nearest-neighbor jumps: Let Γ1 and Γ2 be the total jump frequencies to nearest and next-nearest neighbors, respectively. Show that

Remarks:

a. Equations (1.12) and (1.14) are sometimes called Fick’s first and second laws.³ This is but nomenclature. It is important to remember, however, that Eq. (1.13), the so-called continuity equation, is a precise statement of mass conservation, as we shall soon see. On the other hand, the linear relation (1.12) between the flux J and the thermodynamic force ∂C/∂x, sometimes called a constitutive relation, does not have the same generality. Its form depends on the details of the diffusion mechanism. Note, nonetheless, that its minus sign means that material must flow from regions of high concentration to regions of lower concentration, a property that we ordinarily associate with diffusion. We will refer, henceforth, to Eq. (1.14) as the diffusion equation; it is a consequence of the two previous equations if the diffusivity is a constant. Indeed, any one of Eqs. (1.12)–(1.14) is a consequence of the two other equations.

b. The diffusion equation is sometimes also called the heat equation,⁵,⁶ because the evolution of energy also obeys a conservation law. From a mathematical point of view, it is the simplest so-called parabolic" equation (more nomenclature). From the theory of partial differential equations⁷ we learn that it must be solved as an initial- value problem: Given initial data on C and boundary conditions on the spatial variable x, we can then find the solution by marching forward in time. This will soon become clear in Sects. 1.4 and 1.5 on numerical methods.

c. From Eq. (1.11) we learn that the constancy of D is related to the constancy of the jump distance and frequency. If, for any reason, these are not constants (see the next exercise), then Eq. (1.14) must be replaced by the more general

d. Physical dimensions are important, especially if we wish to cast a given problem in as few variables as possible.N2 We will always denote dimensions by square brackets. Thus, from Eq. (1.11), we learn that [D] = cm²/s.† The dimensions of C and J depend on that choice for Νi . For example, if the particles are evenly distributed along planes perpendicular to the direction of diffusion, the usual case of volume diffusion, then we have [Ni] = cm–2, [C] = cm–3, and [J] = cm–2s–1. Can you find these dimensions if the particles are evenly distributed along lines (surface diffusion) or restricted only to sites (edge diffusion)? Note also that mass or molar units, rather than numbers of particles, can be used to measure Ni, with concomitant changes for C and J.

e. Order-of-magnitude estimates of D are instructive. If, for Γ, we take the Debye frequency in a solid, about 10¹³s–1, and if we take roughly 2Å for the jump distance a, then Eq. (1.11) tells us that D 2 × 10–3cm²/s. This is many orders of magnitude larger than what is actually measured in solids.N3 The moral is simply that not all atomic oscillations lead to successful jumps: Only big ones succeed, as is true in other random walks of life.

f. The careful reader will have noticed that the hat has come off the symbol J for the flux somewhere between Eqs. (1.9) and (1.12). This is not an oversight. There are two points to bear in mind regarding fluxes: First, they are directed quantities; the function J must then be the component of a flux vector by dividing by appropriate lengths or areas. This remark is related to remark (d), above; a fuller account is provided in Sect. 1.6 and in Appendix A.

EXERCISE 1.4. A model for variable diffusivity: Write first the master equation (1.4) restricted to nearest-neighbor jumps. Then assume that the transition frequencies depend on the occupation of the final sites: Γi,i Γfi ± 1 for all i, where the dimensionless and arbitrary function fi±1 depends only on the occupation numbers Ni±1 of adjacent sites. Show that Eq. (1.15) follows from Eq. (1.4), with a concentration-dependent diffusivity D(CΓa²f(1 – d ln f /d ln C

1.3. Elementary Properties of the Diffusion Equation

The diffusion equation enjoys several interesting properties⁴–⁷ that are useful, either to characterize the behavior of the solutions or to derive new solutions from known ones. These properties can be subdivided into three categories: conservation, smoothing, and invariance.

A.1 Conservation in Continuous Regions

The conservative nature of Eqs. (1.13) and (1.14)–(1.15), or of their discrete analogs (1.1)–(1.2), has already been emphasized. We are now ready for a disarmingly simple proof. We consider a diffusion process in a stationary slab, and denote by a and b the coordinates of its left- and right-hand boundaries, respectively. The total mass† contained in this region is then

Taking its time derivative and using Eq. (1.13), we get

This string of equalities can be viewed in two ways: Either the container (the slab’s boundary) is impermeable, in which case the two boundary fluxes in the above expression are zero, and the total mass M is constant, or the constancy of M implies, at the very least, that whatever enters the container must also leave. This characterizes the steady state, to which we shall return in the next chapter. The very same operations can be performed on the discrete form (1.2). Note that they are independent of the constitutive relations (1.3) and (1.12). At the end of this chapter we shall see a converse statement, namely, global mass conservation in continuous regions implies the continuity equation (1.13).

EXERCISE 1.5. Fluxes are useful: Consider a semi-infinite region {x > 0}, initially at the uniform concentration C∞, into which material (e.g., a dopant) can diffuse from its free surface. Call M(t) the additional mass that has accrued at time t. Show that this mass can be calculated as the integral

where J0 is the surface flux. Note that this result is independent of the boundary conditions at x

A.2 Conservation at Fixed Surfaces of Discontinuity

The previous calculations assumed that there were no points of discontinuity within the slab. What if that slab consists of two distinct materials or of two distinct phases of the same substance, each characterized by its own diffusivity? If x = ξ is such a point separating the two materials or phases, then it is only necessary that we break up the integral (1.16) into two parts:

For simplicity, we assume that the fluxes at the external boundaries, J(1)(a, t) and J(2)(b, t) either vanish or are equal and that the total mass is constant.† Then, manipulations similar to those leading to Eq. (1.17) yield‡

In other words, the flux is always continuous at stationary internal boundaries as long as these boundaries are not the seat of other processes that change the number of particles. The concentration distribution is not smooth, however. With Eq. (1.12) we have

which shows that the concentration profile suffers a change in slope if the diffusivities are unequal.

A.3 Conservation at a Moving Boundary

Not all boundaries are fixed. In fact, some of the more interesting ones execute motions, as we will see in the examples of Chaps. 2, 4, and 5. It is then required, in the last calculations, to relax one assumption: The internal boundary point’s position is now an arbitrary function of time, x = ξ(ή. The integrals (1.18) then have variable limits, and, using Leibniz’s ruleN4 for such cases, we easily get a generalization of Eqs. (1.19):

C – J . Concentration fields often do suffer discontinuities at phase boundaries, although, in general, the temperature field for analogous melting and freezing problems does not.L, where L is the latent heat per unit volume. Relations such as Eq. (1.20) are called Stefan conditions; they are, in fact, equations of motion for moving boundaries. Any failure to satisfy Stefan conditions is equivalent to a violation of mass or energy conservation.N6

EXERCISE 1.6. Geometric interpretation of flux continuity: Show that the flux continuity condition

B. Smoothing Property

The diffusion equation (1.14) has an interesting geometric interpretation if we think of the concentration distribution as a surface C = C(x, t) over the plane of the independent variables (x, t). The time derivative ∂C/∂t is proportional to the local curvature ∂²C/∂x² along isochronal lines (we usually call these profiles) of the surface.

Fig. 1.2. Profile of a bumpy concentration surface. It tends to smooth out with time.

Therefore, the distribution’s behavior can be inferred from C(x, t + ∆t) ≈ C(x, t) + ∆t∂C/∂t. Figure 1.2 shows a snapshot of that surface. Now, because the diffusivity is positive, bumps with negative curvature will tend to move down in time, and positive curvature bumps will move up. In other words, a bumpy profile becomes smoother as time progresses.

Figure 1.3 shows a particular concentration surface that we will discuss at length in Chap. 4. The vectors (0, 1, ∂C/∂t), tangent to lines of constant x and directed toward increasing values of t, indeed point downward where the curvature is negative. Can you make similar statements for Eq. (1.15)?

C.1 Linearity of the Diffusion Equation

Any constant is a trivial solution of the diffusion equation (1.14), as is any first-degree polynomial in x. More generally, the diffusion equation satisfies a principle of superposition: If C1 and C2 are two solutions of this equation and if α1 and α2 are two arbitrary constants, then α1C1 + α2C2 is also a solution, provided that the initial and boundary conditions of the full problem can be decomposed into those of the partial problems for C1 and Cis a linear operator acting on a solution C, such as differentiation or integration with respect to the independent variables (x, t(C) is also a solution. For example, the flux J, derived from any solution C, is also a solution of the diffusion equation. These facts are easy to verify directly from Eq. (1.14) because the operations in question are linear and they commute. Can you make similar statements for Eq. (1.15)?

Fig. 1.3. Concentration distribution represented as a surface over the (x, t)-plane. This is, in fact, the Gaussian (4.3) of Chapter 4. Note the essential singularity at t = 0. [See Note 2 of Chap. 4.]

EXERCISE 1.7. Integrating solutions: Contrary to what was stated above, the integral of a solution is not always a solution. Let u(x, t) be any solution of Eq. (1.14), and construct the new function

where a is some fixed coordinate. Calculate its derivatives, and find what additional condition one must impose so that U will also be a solution of the diffusion

C.2 Symmetries: Invariance Under Reflections and Translations

The diffusion equation does not change its form under reflection operations in space. More precisely, if C(x, t) is a solution and if x → x = – x, then C(x, t) is also a solution, provided that the boundary conditions are unchanged under that transformation. This is easily verified (using the chain rule) because Eq. (1.14) is second order in x. Similarly, this equation is invariant under translations of the independent variables, i.e., the addition of an arbitrary constant to x and t does not affect its form. A reversal of the time arrow, however, does not preserve the form of the diffusion equation because it is only first order in t. This is the mathematical content of irreversibility, and diffusion is a prime example of an irreversible process. These results are amplified in Appendix Β and in the exercise that follows.

EXERCISE 1.8. Affine transformations of the diffusion equation: Consider transformations of the independent variables (x, t) of the form

where all the coefficients α, β, γ, δ, x0, and t

The diffusion equation has many other interesting properties among which are the maximum principle and the proof of the existence and uniqueness of solutions. This last, in particular, is not merely of esthetic appeal. In practice, if we have found a solution to a particular diffusion problem — by whatever means — then we are assured of having found the solution of the problem in question. The interested reader will find a full account in Refs. 6 and 7, at least for linear problems, but the properties listed here are sufficient for our purposes. We shall