DIGM: MOTION RESULTING FROM DIFFUSION WITHIN AN INTERFACE

J. W. CAHN, P. C. FIFE, and O. PENROSE

Abstract. A well-known but little understood experimental phenomenon called diffusion induced grain boundary motion is given a theory based on a phase-field model with order parameter and concentration of a diffusing solute as field variables. Suggested physical mechanisms for producing the motion are investigated and assessed in the framework of the model.

1 Introduction

Interfaces separating crystalline grains in polycrystalline solids, such as metal and ceramics, are sometimes observed to migrate, causing some grains to grow at the expense of others. When this is motion by curvature, it is called grain growth. But in the presence of chemical composition gradients, there is another cause for the motion that has been found to be strong enough to overwhelm motion by curvature, when curvatures are less than approximately 10^{6}m^{−1}. This phenomenon, called diffusion-induced grain boundary motion (DIGM), was first discovered by Hillert and Purdy in 1977 [4]. A particularly simple geometry for observing this phenomenon is a thin metallic film containing transverse planar grain boundaries, immersed in a vapor containing a substance that can be incorporated in the crystals as a solute in a solid solution. Although DIGM has been studied extensively and various explanations have been proposed, the true mechanism underlying this type of boundary migration remains something of a mystery.

It is the purpose of this note to report on progress made on a theoretical basis for the observed coupling of the motion of such a grain boundary with the solute diffusion field within the boundary itself. We consider a phase-field model for an abstract diffuse interface that has many of the attributes of a grain boundary. We take the boundary to be a thin mobile zone of disorder (grain boundary) in which the solute diffuses rapidly, between two highly ordered regions (crystalline grains) with low diffusivities. The ordered regions support stress, but this stress is not transmitted through the zone where the ordering is low. Since grain boundaries are surfaces between grains of the same phase, differing only in the orientation of the crystal axes, the properties of both crystalline grains are represented by only the solute concentration and the magnitude of an order parameter.

Two coupling mechanisms are considered within the context of our phase field model: the dependence of the stress in the ordered grains on the concentration of the diffusing element (Sec. 6) and the dependence of the interfacial free energy on this same concentration (Sec. 5). Both effects are found to produce boundary motion under certain conditions, and the velocity is estimated in terms of the other parameters and data entering the model, such as the diffusivity, the initial concentration in the crystal, and that in the surrounding vapor.

2 A phase field approach

We use a diffuse-interface model with an order parameter ϕ and a concentration u as the only field variables. Moreover, the interface we consider is straight and represents a boundary between two grains of the same phase. In particular, we exclude at the outset any influence that curvature of the boundary or temperature variation may have on the motion. In reality, these boundaries are usually curved; that additional effect has just not as yet been included in our model.

Typically, the concentration u could represent the fraction of available lattice sites which are filled by solute atoms, and so we shall take its range to be the interval [0, 1]. Similarly −1 ≤ ϕ ≤ 1, the two extremes ±1 representing the two perfect grains on either side.

As is usual with diffuse-interface models (e.g. [2]), the postulated evolution system can be understood as a gradient flow with respect to a functional. In this case, the functional is appropriately named the free energy F [ϕ, u], since temperature will be constant. It has contributions coming individually from ϕ and u, plus a small contribution resulting from the interaction between the two fields:

F [ϕ, u] = F_{1}[ϕ] + F_{2}[u] + ϵF_{12}[ϕ, u].  (1)

The first term will be of standard van der Waals-Cahn-Hilliard form

F_{1}[ϕ] = A ∫_{Ω}[w(ϕ) + \frac{1}{2}δ^{2}|∇ϕ|^{2}]dv.  (2)

Here Ω is the domain in space occupied by the solid material. Rather, assuming everything is independent of one space coordinate, we take it to be the strip {−∞ < x < ∞} × {−h < y < h}, representing a 2D cross section of the metal film. The grain boundary will be a subregion of Ω consisting of a slim rectangle Ω_{0} bounded on top and bottom by the two faces {y = ±h}. The law of motion of this rectangle is our primary object of interest. The parameter δ^{2} is some measure of the relative weight attached to gradients of ϕ in the expression for the free energy (we shall see that δ is also a measure of the thickness of the interface Ω_{0}), and w is approximately the free energy density if ϕ were constant in space. We take it to be the following double-well potential with barrier at |ϕ| = 1 (see, e.g. [1]):

w(ϕ) = \frac{1}{2}(1 − ϕ^{2}) if |ϕ| < 1
and = +∞ if |ϕ| > 1  (3)

For the values |ϕ| = 1, w(ϕ) is left undefined.

We take

F_{2}[u] = ∫_{Ω} W (u)dv,  (4)

where W is a convex function representing the free energy density of the solid due to the presence of the solute. A typical form for W is

W (u) = const · [u ln u + (1 − u) ln(1 − u)].  (5)

Other effects of the type we are seeking can be obtained by including terms involving ∇u in (4); we shall omit them for simplicity.

Our object will be to discover interaction terms F_{12} which can serve as impetus for the motion of the interface and which have physical interpretations. They will be of the assumed form

F_{12} = ϵA ∫_{Ω} p(ϕ, u)dv.  (6)

The small dimensionless normalizing parameter ϵ is chosen so that the derivatives of p are of order unity. Setting ϵ̄ = \frac{ϵA}{W′′(1/2)}, we assume the smallness relations

|ϵ| ≪ 1; |ϵ̄| ≪ 1.  (7)

Our chosen kinetic equations take the form

τ \frac{∂ϕ}{∂t} = −\frac{δF}{δϕ} , \frac{∂u}{∂t} = ∇ · D(ϕ)∇u(1 + O(ϵ̄))  (8)

(the error term depends on p and will not be given explicitly), where

\frac{δF}{δϕ} = −δ^{2}Δϕ + w′(ϕ) + ϵA\frac{∂p}{∂ϕ},

D(ϕ) represents the diffusivity of the solute as a function of the degree (measured by ϕ) to which the material is ordered, and τ is a relaxation time.

One interpretation and/or derivation of the system (8) involves its characterization as a gradient (steepest descent) flow for the functional F [ϕ, u] with respect to a certain D-dependent scalar product, at least in the case that D > 0 and W is given by (5).

We are mostly interested, however, in the case when D(ϕ(x, t)) vanishes in much of the domain of interest; the gradient nature of the flow then degenerates. It is well known, in fact, that the diffusivity is much larger in the grain boundary than in a perfect crystal. We shall make the simplifying assumption that D(ϕ) = 0 in the crystals, that is when |ϕ| = 1. Our specific choice for D is

D(ϕ) = D(0)(1 − ϕ^{2}),  (9)

where D(0) is a constant.

By the choice (3) of w, the system (8) now becomes

τ \frac{∂ϕ}{∂t} = ϕ + δ^{2}∇^{2}ϕ − ϵ∂p/∂ϕ if |ϕ| ≠ 1,  (10)

\frac{∂u}{∂t} = ∇ · D(ϕ)∇u(1 + O(ϵ̄)).  (11)

To account for the possibility that |ϕ| = 1, we extend (10) as a variational inequality (see [1] for an analogous situation). The result is that in places where |ϕ| = 1, (10) need not hold; all that is necessary there is for ϕ and its gradient to be continuous.

As boundary conditions on ∂Ω, physical considerations suggest that

\frac{∂ϕ}{∂n} = 0 on ∂Ω,  (12)

where ∂/∂n denotes differentiation in the direction normal to ∂Ω, and (neglecting ϵ̄)

u = u^{∗} on ∂Ω,  (13)

where u^{∗} is a given value depending on the properties of the external vapor. Conditions at x = ±∞ are also needed. We require

ϕ = ±1 and \frac{∂u}{∂t} = 0  (14)

in regions to the left and right in our strip domain Ω where the pure crystals reside, i.e. where D = 0. These regions are not known at the outset.

3 Travelling wave solutions

To model a moving grain boundary we look for a leftward travelling wave solution of the equations, in which ϕ and u depend on x and t only through the combination x + ct, where the positive velocity c is to be determined. We shall use the notation z = x + ct and look for solutions of (10) - (14) in the form ϕ = ϕ(z, y), u = u(z, y). We also prescribe the resident value u_{0} of u in the eroding grain on the left. In all, we get, in addition to (12) and (13),

cτϕ_{z} = ϕ + δ^{2}(ϕ_{zz} + ϕ_{yy}) − ϵ(∂p/∂ϕ),  (15)
cu_{z} = ∇ · D(ϕ)∇u(1 + O(ϵ̄));  (16)

u(z, y) = u_{0}, ϕ(z, y) = −1, z large negative ,
ϕ(z, y) = 1, z large positive  (17)

Assuming that a travelling wave exists, we see that (15) can be used to obtain an expression for c in terms of p(ϕ, u), hence in terms of u, which will be useful in the following sections. We set

g(z, y) = −\frac{∂p}{∂ϕ}[ϕ(z, y), u(z, y)],  (18)

multiply (15) by ϕ_{z}, integrate with respect to z and y and use the boundary conditions (12), (17) to obtain

cτ ∫ ∫ ϕ_{z}^{2}dzdy = ϵ ∫ ∫ ϕ_{z}gdzdy.  (19)

Physically, this is surprising: when the coupling ϵg is absent, no traveling wave motion is possible, even though there would be a decrease in the total free energy if the boundary were to move.

4 An approximate solution

We formally simplify the problem (15), (16), (12), (17) by utilizing the smallness of ϵ and ϵ̄. First, (15) with ϵ = 0, c = 0 gives the following lowest order approximation for ϕ :

ϕ^{(0)} = S(x/δ),  (20)

where the function S(x) is defined to equal sin x for |x| < π/2, and is continued to equal ±1 for x outside that interval.

One can now obtain the lowest order approximation to u in terms of the unknown velocity c (which will be determined later). We use (20) in (16) with the error term neglected. By (20) and (9),

D(ϕ(z)) ≈ D(ϕ^{(0)}(z)) = D(S(z/δ))
= \{ \begin{matrix}D(0) cos^{2}(z/δ) if − πδ/2 < z < πδ/2 \\ 0 otherwise . \end{matrix} (21)

The equation to be solved is

cu_{z} = (D(z)u_{z})_{z} + D(z)u_{yy},  (22)

where D(z) stands for D(ϕ^{(0)}(z)).

The domain Ω is split into three parts as shown:

Ω_{0} Ω_{−}  Ω_{+} y = h y = −h z > πδ/2 z < −πδ/2

Fig. 1: The regions Ω_{−}, Ω_{0}, Ω_{+}

In the left grain Ω_{−}, D = 0, so that u_{z} = 0, and from (17),

u(z, y) = u_{0} in Ω_{−}  (23)

In the region Ω_{+} we again have u_{z} = 0, and hence

u(z, y) = u(πδ/2, y) in Ω_{+}  (24)

where u(πδ/2, y) is not yet known (this is the concentration left in the wake of the moving boundary). In the remaining region, Ω_{0} = {−πδ/2 < z < πδ/2, −h < y < h}, the full equation (22) is to be solved. The boundary conditions, determined on {y = ±h} by (13) and on {z = ±πδ/2} by the condition that u and its gradient should be continuous inside the metal, are

u = u∗  (y = ±h)
u = u_{0} and u_{z} = 0 (z = −πδ/2)
u continuous;  u_{z} = 0  (z = πδ/2)  (25)

The number of boundary conditions appears to make the problem overdetermined, but in fact some of them follow automatically from others. This circumstance is due to the very singular nature of (22): on the two lateral sides of Ω_{0}, D and its derivative both vanish.

It turns out that a useful approximation to the solution of (22), (25) can be obtained under the assumption that the following parameter η is small:

η ≡ \sqrt{\frac{δc}{D(0)}} ≪ 1.  (26)

The approximation is

u(z, y) ≈ u^{(0)}(z, y) = u_{0} + (u^{∗} − u_{0})f(z)\frac{cosh yη/δ}{cosh hη/δ} ,  (27)

where f(z) is a positive function with the properties

f(z) ≈ e^{−η^{2} tan \frac{|z|}{δ}} (−πδ/2 ≤ z ≤ 0)
and ≈ 1 (0 ≤ z ≤ πδ/2).  (28)

We shall use this approximate expression for u in all of the following.

Note that if hη/δ ≫ 1, then for most values of z ∈ (−πδ/2, πδ/2) the function u^{(0)}(z, y) as given by (27) changes abruptly from u^{∗} at the top and bottom (y = ±h) to near u_{0} within a distance of the order of δ/η. On the other hand if hη/δ ≪ 1, then the profile is approximately parabolic :

u(z, y) ≈ u^{∗} + \frac{η^{2}}{δ^{2}} (u_{0} − u^{∗})(h^{2} − y^{2}).  (29)

Finally, with knowledge of an approximation u = u^{(0)} obtained this way, a leading order approximation for c may also be obtained. In (19) we set ϕ = ϕ^{(0)} as given by (20) and u = u^{(0)} as given by (27). Using the fact that ∫_{−∞}^{∞}(ϕ_{z}^{(0)})^{2} dz = π/2δ by (20), we obtain

c ≈ \frac{2ϵ}{πτ} ∫_{−πδ/2}^{πδ/2} ḡ(z) cos \frac{z}{δ} dz,  (30)

where ḡ(z) is the y-average of g(z, y). This formula is the main tool in our search for mechanisms to make c > 0. In the following sections it will be applied to different interaction functions g.

5 Coupling by a concentration-dependent interfacial free energy

It was suggested by Meyrick [6] that interface motion in some contexts could be caused by an interaction in which the interfacial free energy depends on the concentration of solute atoms. This interaction is analogous to the Marangoni effect [5, 7], in which an interface moves because the surface tension of a fluid depends on concentration or temperature.

In our context, the interfacial free energy should be interpreted as the term F_{1} in (1). A simple way to model this mechanism is to make the coefficient of w(ϕ) in the formula (2) depend linearly on u. This is equivalent to adding a term of the form (6) in which

p = uw(ϕ).  (31)

so that, by (3), g(u, ϕ) = −∂p/∂ϕ = uϕ for |ϕ| < 1.

The coefficient ϵ in (6) can be of either sign. If ϵ > 0, the interaction free energy density ϵAp increases with u at fixed ϕ, and so it is expected that the interface will tend to move towards smaller values of u. If ϵ < 0, it would tend to move the other way.

With ϕ^{(0)} given by (20) and u^{(0)} by (27), formula (30) now yields

c ≈ \frac{ϵδη^{2}}{2τ} (u^{∗} − u_{0})\frac{tanh \sqrt{2/π}hη/δ}{\sqrt{2/π}hη/δ} .  (32)

Our setup requires c ≥ 0. For the Meyrick mechanism, the boundary can therefore move only if u^{∗} −u_{0} has the same sign as ϵ. Even if the signs are right, there is no guarantee that (32) can be solved for c. The right side also depends on c through η (26), resulting in an implicit equation for c. If a solution with c > 0 exits, one factor “c” may be eliminated from both sides. The equation for c so obtained has a solution only if (u^{∗} − u_{0})ϵδ^{2}/τD(0) is large enough.

However, experiments [3] have shown that for many alloy systems, it is possible for the grain boundary to move for either sign of u^{∗} −u_{0} (The direction of motion of the boundary depends on the sign of u^{∗} − u_{0}). So a term of the form (31) cannot provide the sole driving force for DIGM.

6 Coupling by stress energy.

In this model, we depend on the assumption that the travelling wave is fully developed, so that the dissolving grain, occupying the region {z < −πδ/2}, has a given limiting concentration u_{0} as z→ − ∞, and the grain being formed in the wake of the moving interface has a limiting concentration profile u_{+}(y) as z→ + ∞, which is unknown at the moment.

At any location, the equilibrium lattice spacing of a perfect crystal is generally determined by the solute concentration u, and therefore a spatial variation of u in either grain will induce stress, whose energy should be incorporated into the free energy functional F. In the dissolving grain, the stress energy density in our approximation will take the form ϵAP_{−}(u − u_{0})^{2}, with ϵAP_{−} positive. In the advancing grain, which we suppose to be disconnected as explained below from the other grain, it is ϵAP_{+}(u−u_{+})^{2}. In the interface where any crystalline structure is imperfect, the transmission of stress due to concentration differences is weakened. This fact is the source of a coupling between the two fields u and ϕ, and hence of a stress-induced term F_{12}.

We model the ability to transmit stress in this way as an all-or-none phenomenon. No stress can be transmitted to neighboring portions of the metal across a sufficiently disordered layer, say when |ϕ| < ϕ_{1} for some ϕ_{1} ∈ (0, 1), but stresses due to nonconstant u in the dissolving (advancing) grain occur as though it were a perfect crystal occupying the entire region {z : ϕ(z) < −ϕ_{1} (> ϕ_{1})}.

Since there will be no stress energy in the region {z : |ϕ(z)| < ϕ_{1}}, the scaled coupling free energy density p (see (6)) can be written

p(z, y) = P_{−}(u − u_{0})^{2}H(−ϕ − ϕ_{1}) + P_{+}(u − u_{+})^{2}H(ϕ − ϕ_{1}),  (33)

where H is the Heaviside function (unit step function).

Expressing the derivative of H as a delta-function, we obtain

ḡ(z) = −\overline{\frac{δp}{δϕ}} = δ(ϕ + ϕ_{1})P_{−}\overline{(u − u_{0})^{2}} − δ(ϕ − ϕ_{1})P_{+}\overline{(u − u_{+})^{2}}  (34)

where the bars indicate averages over −h < y < h at fixed z.

After some tedious calculations using (20), (27), and (30), we find under certain assumptions on ϕ_{1} and the other parameters that

c ≈ \frac{2ϵδ}{πτ} P_{−}(u_{0} − u^{∗})^{2}.  (35)

Other expressions are obtained from different assumptions, but in all cases, c is an even function of (u^{∗} − u_{0}), so that a reversal of the sign of (u^{∗} − u_{0}) has no effect on the velocity. This appears to conform to some experiments, but more investigation is needed. We conclude that this effect is a feasible impetus for the motion.

A possibly more realistic case is that in which both effects in this and the last section occur, with different values of ϵ. Then an expression for c depending on these two values can easily be obtained.

7 Acknowledgements

We thank Gary Purdy for many stimulating discussions which caught the attention of PF and OP during an ICMS workshop at Heriot-Watt University, and for his continued interest and critical comments during this research. The workshop was supported by the SERC,

the European Community and the London Mathematical Society. JC is grateful for the invitation to come to Heriot-Watt and for support from ARPA. PF is grateful to Heriot Watt. His research was also supported by NSF Grant 9201714.

References

[1] J. F. Blowey and C. M. Elliott, “The Cahn-Hilliard gradient theory for phase separation with non-smooth free energy. Part I: Mathematical analysis,” Euro. J. Appl. Math. 2, 233-279 (1991).

[2] P. C. Fife, “Barrett Lecture Notes, Univ. of Tennessee,” 1991.

[3] C. A. Handwerker, “Diffusion-induced grain boundary motion in thin films,” in “Diffusion phenomena in thin films and microelectronic materials,” edited by D. Gupta and P. S. Ho, Noyes Publishing Co., Park Ridge, NJ, 1988, 245-322.

[4] M. Hillert & G. R. Purdy, “Chemically induced grain boundary motion,” Acta Metall. 26 333-240 (1978).

[5] Carlo G. M. Marangoni, Ann. Phys. Poggendorf 143, 337 (1871).

[6] G. Meyrick, “On the initiation of discontinuous precipitation,” Scripta Metall. 10 649-653 (1976).

[7] L. E. Scriven and C. V. Sternling, “The Marangoni Effects,” Nature 187, 186-188 (1960).

J. W. Cahn NIST, Bldg 223, Rm A153, Gaithersburg, MD 20899, USA

P. C. Fife Mathematics Department, University of Utah Salt Lake City, UT 84112, USA

O. Penrose Mathematics Department, Heriot-Watt University Riccarton, Edinburgh EH14 4AS, U.K.