%Levico xxFeb 3 \documentstyle{article} % To remove the boxes showing the equation codes, move the percent % sign from the lower to the upper of the following pair of lines. %\newcommand{\lbl}[1]{\qquad\qquad\framebox{\scriptsize #1} \label{#1}} \newcommand{\lbl}[1]{\label{#1}} \def\be{\begin{equation}} \def\bea{\begin{eqnarray}} \def\ee{\end{equation}} \def\eea{\end{eqnarray}} \def\half{{1 \over 2}} \def\bF{{\bar{F}}} \def\bD{{\bar{D}}} \def\by{{\bar{y}}} \def\bx{{\bar{x}}} \def\bu{{\bar{u}}} \def\ds{{\displaystyle}} \def\div{{\rm div}} \def\grad{{\bf grad}} \def\1{{\pi}} \def\q{{\gamma}} \def\Q{{\Gamma}} \def\w{{\delta}} \def\W{{\Delta}} \def\e{{\epsilon}} \def\r{{\theta}} \def\R{{\Theta}} \def\t{{\tau}} \def\y{{\ypsilon}} \def\Y{{\Ypsilon}} \def\u{{\xi}} \def\U{{\Xi}} \def\i{{\iota}} \def\p{{\rho}} \def\a{{\alpha}} \def\A{{\nabla}} \def\s{{\sigma}} \def\S{{\Sigma}} \def\d{{\phi}} \def\D{{\Phi}} \def\g{{\lambda}} %\def\G{{\Lambda}} \def\h{{\eta}} \def\l{{\omega}} \def\L{{\Omega}} \def\z{{\zeta}} \def\Z{{\approx}} \def\X{{\equiv}} \def\c{{\psi}} \def\C{{\Psi}} \def\b{{\beta}} \def\B{{\inf}} \def\l{{\lambda}} \def\n{{\nu}} \def\N{{\nonumber \\}} \def\m{{\mu}} \def\M{{\partial}} \def\ss{{\cal }} \def\hh{{\hat{ }}} \def\-{{\bar}} \def\.{{\dot}} \def\tt{{\tilde}} \def\inf{{\infty}} \def\bbr{{{\Bbb{R}}}} \def\ra{{\rightarrow}} \def\da{{\downarrow}} \def\ua{{\uparrow}} \def\'{{\prime}} %\includeonly{cfp1} %\includeonly{cfp2} %\include=input \begin{document} \vspace*{2cm} \begin{center} {\large \bf DIGM: MOTION RESULTING FROM DIFFUSION WITHIN AN INTERFACE} \vspace{1cm} {\sc J. W. CAHN, P. C. FIFE, and O. PENROSE} \end{center} \vspace{3cm} {\bf 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. \vspace{1.5cm}\\ \section{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}$. %may or may not be %accompanied by changes in the chemical composition; in any case, it %can have important practical consequences. This phenomenon, called diffusion-induced grain boundary motion (DIGM), was first discovered by Hillert and Purdy in 1977 \cite{HP}. A particularly %interesting type %of boundary motion is sometimes induced, or greatly accelerated, simple geometry for observing this phenomenon is a thin metallic film containing transverse planar grain boundaries, immersed in a vapor containing a substance that %which 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 %will consider a phase-field model for an abstract diffuse interface that %will have has many of the attributes of a grain boundary. We take the boundary to %It will 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. \ref{cbse}) and the dependence of the interfacial free energy on this same concentration (Sec. \ref{cdife}). 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. \section{A phase field approach} \lbl{kebc} We use a diffuse-interface model with an order parameter $\d$ 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 \le \d \le 1$, the two extremes $\pm 1$ representing the two perfect grains on either side. As is usual with diffuse-interface models (e.g. \cite{fife}), 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[\d,u]$, since temperature will be constant. It has contributions coming individually from $\d$ and $u$, plus a small contribution resulting from the interaction between the two fields: \be F [ \phi, u ] = F_1 [ \phi ] + F_2[ u ] + \epsilon F_{12} [ \phi, u ]. \lbl{F} \ee The first term will be of standard van der Waals-Cahn-Hilliard form \be F_1 [ \phi] = A\int_{\Omega} [ w( \phi) + \frac{1}{2} \w^2 | \nabla \phi |^2 ] dv. \lbl{F1} \ee Here $\L$ 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 $\{-\B < x < \B\}\times \{-h < y < h\}$, representing a 2D cross section of the metal film. The grain boundary will be a subregion of $\L$ consisting of a slim rectangle $\L_0$ bounded on top and bottom by the two faces $\{y = \pm h\}$. The law of motion of this rectangle is our primary object of interest. The parameter $\w^2$ is some measure of the relative weight attached to gradients of $\d$ in the expression for the free energy (we shall see that $\w$ is also a measure of the thickness of the interface $\L_0$), and $w$ is approximately the free energy density if $\d$ were constant in space. We take it to be the following double-well potential with barrier at $|\d| = 1$ (see, e.g. \cite{elliott}): \bea w(\phi) & = & {1 \over 2} (1 - \phi^2) \hbox{~~if~~} | \phi| < 1 \nonumber \\ \hbox{and} & = & + \infty \hbox{~~if~~} |\phi| > 1 \lbl{w} \eea For the values $|\phi| = 1, ~w(\phi)$ is left undefined. We take \be F_2[u] = \int_{\Omega}W(u) dv, \lbl{F2} \ee 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 \be W(u) = \mbox{const}\cdot [u \ln u + (1-u) \ln (1-u)]. \lbl{W2} \ee Other effects of the type we are seeking can be obtained by including terms involving $\A u$ in (\ref{F2}); 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 \be F_{12} = \epsilon A \int_{\Omega} p(\phi, u) dv. \lbl{F12} \ee The small dimensionless normalizing parameter $\epsilon$ is chosen so that the derivatives of $p$ are of order unity. Setting $\bar{\e} = \frac{\e A}{W''(1/2)}$, we assume the smallness relations \be |\epsilon| \ll 1;~~|\bar{\e}| \ll 1. \lbl{eps} \ee Our chosen kinetic equations take the form \be \t\frac{\M\d}{\M t} = - \frac{\w F}{\w \d},~~\frac{\M u}{\M t} = \A\cdot D(\d)\A u(1 + O(\bar{\e})) \lbl{ev}\ee (the error term depends on $p$ and will not be given explicitly), where $$\frac{\w F}{\w \d} = -\w^2\W\d + w'(\d) + \e A \frac{\M p}{\M\d}, $$ $D(\d)$ represents the diffusivity of the solute as a function of the degree (measured by $\d$) to which the material is ordered, and $\t$ is a relaxation time. One interpretation and/or derivation of the system (\ref{ev}) involves its characterization as a gradient (steepest descent) flow for the functional $F[\d,u]$ with respect to a certain $D$-dependent scalar product, at least in the case that $D > 0$ and $W$ is given by (\ref{W2}). We are mostly interested, however, in the case when $D(\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(\phi) = 0$ in the crystals, that is when $|\phi| = 1$. Our specific choice for $D$ is \be D(\phi) = D(0)(1-\phi ^2), \lbl{Dphi} \ee where $D(0)$ is a constant. By the choice (\ref{w}) of $w$, the system (\ref{ev}) now becomes \be \tau \frac{\partial \phi}{\partial t} = \phi + \w^2 \nabla^2 \phi - \epsilon \partial p/ \partial \phi ~~\hbox{ if }~~| \phi| \ne 1, \lbl{dphi/dt} \ee \be \frac{\partial u}{\partial t} = \A\cdot D(\d) \A u (1 + O(\bar{\e})). \lbl{du/dt2} \ee To account for the possibility that $|\d| = 1$, we extend (\ref{dphi/dt}) as a variational inequality (see \cite{elliott} for an analogous situation). The result is that in places where $| \phi| = 1$, (\ref{dphi/dt}) need not hold; all that is necessary there is for $\phi$ and its gradient to be continuous. As boundary conditions on $\M\L$, physical considerations suggest that \be \frac{\partial \phi}{\partial n} = 0 {\rm ~on~} \partial \Omega, \lbl{phibc} \ee where $\partial / \partial n$ denotes differentiation in the direction normal to $\partial \Omega$, and (neglecting $\bar{\e}$) \be u = u^* \mbox{ on }\M\L, \lbl{vapor} \ee where $u^*$ is a given value depending on the properties of the external vapor. Conditions at $x = \pm \B$ are also needed. We require \be \d = \pm 1\mbox{ and } \frac{\M u}{\M t} = 0 \lbl{inf} \ee in regions to the left and right in our strip domain $\L$ where the pure crystals reside, i.e. where $D = 0$. These regions are not known at the outset. \section{Travelling wave solutions} \lbl{tws} To model a moving grain boundary we look for a leftward travelling wave solution of the equations, in which $\phi$ 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 (\ref{dphi/dt}) - (\ref{inf}) in the form $\phi=\phi(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 (\ref{phibc}) and (\ref{vapor}), \bea c\tau \phi_z & = & \phi + \w^2(\phi_{zz}+\phi_{yy}) - \epsilon (\partial p / \partial \phi ), \lbl{twphi}\\ cu_z & = & \A\cdot D(\d)\A u(1 + O(\bar{\e})); \lbl{twu} \eea \bea u(z,y) = u_0,~~\d(z,y) = -1,~~~~z\mbox{ large negative }, \nonumber \\ \d(z,y) = 1,~~~~z\mbox{ large positive } \lbl{infbc2} \eea Assuming that a travelling wave exists, we see that (\ref{twphi}) can be used to obtain an expression for $c$ in terms of $p(\d,u)$, hence in terms of $u$, which will be useful in the following sections. We set \begin{equation} g(z,y) = -\frac{\partial p}{\partial \phi}[\d (z,y), u(z,y)], \lbl{ar} \end{equation} multiply (\ref{twphi}) by $\d_z$, integrate with respect to $z$ and $y$ and use the boundary conditions (\ref{phibc}), (\ref{infbc2}) to obtain \begin{equation} c\t\int\int \d_z^2 dzdy = \epsilon \int\int\d_zgdzdy.\lbl{ao} \end{equation} Physically, this is surprising: when the coupling $\e 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. \section{An approximate solution} $\lbl{as}$ We formally simplify the problem (\ref{twphi}), (\ref{twu}), (\ref{phibc}), (\ref{infbc2}) by utilizing the smallness of $\e$ and $\bar{\e}$. First, (\ref{twphi}) with $\epsilon = 0, ~c=0$ gives the following lowest order approximation for $\phi$ : \be \phi^{(0)} = S(x/\delta), \lbl{phi0} \ee where the function $S(x)$ is defined to equal $\sin x$ for $|x| < \pi/2$, and is continued to equal $\pm 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 (\ref{phi0}) in (\ref{twu}) with the error term neglected. By (\ref{phi0}) and (\ref{Dphi}), \bea D(\phi(z)) \approx D(\phi^{(0)}(z)) & = & D(S(z/\delta)) \nonumber \\ & = & \left\{ \begin{array}{l} D(0) \cos^2(z/\delta) \hbox{ if } -\pi\w/2 < z < \pi\w/2 \\ 0 \hbox{ otherwise }. \\ \end{array} \right. \lbl{D1} \eea The equation to be solved is \be cu_z = (D(z)u_z)_z + D(z)u_{yy}, \lbl{mueqn} \ee where $D(z)$ stands for $D(\phi^{(0)}(z))$. The domain $\L$ is split into three parts as shown: \vskip 40 pt \begin{picture}(200,100)(-50,-20) \thicklines \put(0,0){\line(1,0){200}} \put(0,80){\line(1,0){200}} \put(90,0){\line(0,1){80}} \put(110,0){\line(0,1){80}} \put(95,40){$\Omega_0$} \put(0,40){$\Omega_-$} \put(170,40){$\Omega_+$} \put(85,85){$y = h$} \put(80,-10){$y = -h$} \put(115,40){$z > \pi\w/2$} \put(30,40){$z < -\pi\w/2$} %the coordinate axes %\thinlines %\put(-20,40){\line(1,0){240}} %\put(100,-20){\line(0,1){120}} %\put(220,40){\line(-1,1){8}} %\put(220,40){\line(-1,-1){8}} %\put(230,40){$z$-axis} %\put(100,100){\line(1,-1){5}} %\put(100,100){\line(-1,-1){5}} %\put(900,105){$y$-axis} \end{picture} Fig. 1: The regions $\Omega_-, \Omega_0, \Omega_+$ \vskip 20 pt In the left grain $\Omega_-,~D = 0$, so that $u_z=0$, and from (\ref{infbc2}), \be u(z,y) = u_0 \hbox{ in } \Omega_- \lbl{uinOm-} \ee In the region $\Omega_+ $ we again have $u_z=0$, and hence \be u(z,y)=u(\pi\w/2,y) \hbox{ in } \Omega_+ \lbl{uinOm+} \ee where $u(\pi\w/2,y)$ is not yet known (this is the concentration left in the wake of the moving boundary). In the remaining region, $\Omega_0 = \{-\pi\w/2 0$. In the following sections it will be applied to different interaction functions $g$. \section{Coupling by a concentration-dependent interfacial free energy} \label{cdife} It was suggested by Meyrick \cite{Meyrick} 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 \cite{marangoni,scriven}, 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 (\ref{F}). A simple way to model this mechanism is to make the coefficient of $w(\phi)$ in the formula (\ref{F1}) depend linearly on $u$. This is equivalent to adding a term of the form (\ref{F12}) in which \be p = uw(\d). \lbl{Mey} \ee so that, by (\ref{w}), $g(u,\phi) = -\partial p/ \partial \phi = u\phi$ for $|\d| < 1$. The coefficient $\epsilon$ in (\ref{F12}) can be of either sign. If $ \e > 0$, the interaction free energy density $\epsilon A p$ increases with $u$ at fixed $\phi$, and so it is expected that the interface will tend to move towards smaller values of $u$. If $\epsilon < 0$, it would tend to move the other way. With $\phi^{(0)}$ given by (\ref{phi0}) and $u^{(0)}$ by (\ref{v}), formula (\ref{c}) now yields \be c \approx {\epsilon \w \h^2 \over 2 \tau}(u^*-u_0) {\tanh { \sqrt{2/\pi} h \h/\w} \over { \sqrt{2/\pi} h \h/\w}}. \lbl{c3M} \ee Our setup requires $c \ge 0$. For the Meyrick mechanism, the boundary can therefore move only if $ u^*-u_0 $ has the same sign as $\epsilon $. Even if the signs are right, there is no guarantee that (\ref{c3M}) can be solved for $c$. The right side also depends on $c$ through $\h$ (\ref{eta}), 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)\e\w^2/\t D(0)$ is large enough. However, experiments \cite{handwerker} 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 (\ref{Mey}) cannot provide the sole driving force for DIGM. \section{Coupling by stress energy.} \label{cbse} In this model, we depend on the assumption that the travelling wave is fully developed, so that the dissolving grain, occupying the region $\{z<-\pi\w/2\}$, has a given limiting concentration $u_0$ as $z \ra -\B$, and the grain being formed in the wake of the moving interface has a limiting concentration profile $u_+(y)$ as $z \ra +\B$, 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 $\epsilon A P_-(u - u_0)^2$, with $\e A P_-$ positive. In the advancing grain, which we suppose to be disconnected as explained below from the other grain, it is $\epsilon A P_+(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 $\d$, 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 $|\d| < \d_1$ for some $\d_1 \in (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 : \d(z) < -\d_1~(> \d_1)\}$. Since there will be no stress energy in the region $\{z : | \d(z)| < \d_1\}$, the scaled coupling free energy density $p$ (see (\ref{F12})) can be written \begin{equation} p(z,y) = P_-(u -u_0)^2 H(-\d -\d_1) + P_+(u - u_+)^2H(\d - \d_1), \lbl{fa} \end{equation} where $H$ is the Heaviside function (unit step function). Expressing the derivative of $H$ as a delta-function, we obtain \begin{equation} \bar{g}(z) = -\overline{\frac{\w p}{\w \d}} = \mbox{ \boldmath $ \delta$}(\d + \d_1)P_-\overline{(u -u_0)^2} - \mbox{ \boldmath $ \delta$} (\d - \d_1)P_+\overline{(u -u_+)^2} \lbl{fb} \end{equation} where the bars indicate averages over $-h