# Computation of Trajectories and Displacement Fields in a Three-Dimensional Ternary Diffusion Couple: Parabolic Transform Method

Authors: Marek Danielewski; Henryk Leszczyński
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The problem of Kirkendall’s trajectories in finite, three- and one-dimensional ternary diffusion couples is studied. By means of the parabolic transformation method, we calculate the solute field, the Kirkendall marker velocity, and displacement fields. The velocity field is generally continuous and can be integrated to obtain a displacement field that is continuous everywhere. Special features observed experimentally and reported in the literature are also studied: (i) multiple Kirkendall’s planes where markers placed on an initial compositional discontinuity of the diffusion couple evolve into two locations as a result of the initial distribution, (ii) multiple Kirkendall’s planes where markers placed on an initial compositional discontinuity of the diffusion couple move into two locations due to composition dependent mobilities, and (iii) a Kirkendall plane that coincides with the interphase interface. The details of the deformation (material trajectories) in these special situations are given using both methods and are discussed in terms of the stress-free strain rate associated with the Kirkendall effect. Our nonlinear transform generalizes the diagonalization method by Krishtal, Mokrov, Akimov, and Zakharov, whose transform of diffusivities was linear.


---

## Body

## 1. Introduction

The Kirkendall effect in solids involves diffusion and distortion [1]. The distortion of a solid material is the change in position with time of a set of inert markers (real or imagined) fixed in the material. The Kirkendall effect means the changes in the observed marker positions as they move in reaction to the locally nonbalanced diffusion fluxes. During one-dimensional diffusion the material deforms unidirectionally. The markers that are initially fixed in a plane perpendicular to the diffusion direction are likely to move all together, defining marker planes. As shown by Stephenson [1], the deformation that accompanies diffusion is described by a stress-free strain rate.Essential for the assessment of the Kirkendall effect is the definition of the reference frame (RF) for the components diffusion fluxes (intrinsic diffusion). Darken postulated defining it as a lattice RF [2]. Inert (or imagined) ideal markers do not take part in the diffusion and their velocity equals the local lattice velocity. It implies that intrinsic diffusion fluxes can be defined with respect to this set of markers. Although the Kirkendall effect has been known for 60 years [2, 3], a number of interesting theoretical and experimental effects continue to be revealed which warrant further considerations, including (i) multiple Kirkendall’s planes, (ii) stability of individual Kirkendall’s planes and (iii) trajectory of Kirkendall’s plane in three-dimensional diffusion couples.Diffusion couples are usually used for experimental and theoretical studies. For a ternary alloy maintained at a constant temperature, the initial state consists of two end members joined with an initial planar discontinuity at a position we take to bez(0). The end members are usually uniform and of different compositions that may also be of different phases. When the couple is held at a fixed temperature, diffusion occurs over time in the z-direction. Experimentally, markers are typically placed only on the initial discontinuity of the diffusion couple; however, they can be placed throughout the sample to measure the entire field of deformation field [4, 5]. The rate of change of the position of the entire set of markers in the sample determines the Kirkendall velocity as a field defined throughout the sample.In the analysis of diffusion couples, there are several quantities of interest, mainly associated with the initial planar discontinuity. The Matano plane is usually defined by proper equal-area constructions established on a concentration profile [6] and remains at z(0) for all t≥0. Kirkendall’s plane(s) problem from the origin of coordinates at t=0 represents the motion of inert markers tracking the deformation; that is, a Kirkendall plane is a material surface. The Matano plane is not a material surface. The Kirkendall plane(s), the material trajectories, and the Matano interface are generally associated with different locations in the sample once diffusion takes place.This paper is stimulated by a variety of theoretical and experimental efforts. The Kirkendall markers, that is, those markers placed initially in the planez(0), are occasionally found on two different planes after heat treatment. The possibility of such multiple Kirkendall’s planes was first described theoretically for single-phase diffusion couples by Cornet and Calais [7], using similarity solutions for the composition and Kirkendall’s velocity fields in a binary diffusion couple. A necessary factor for multiple Kirkendall’s planes is a strong variation of intrinsic diffusion coefficients with composition; for example, multiple Kirkendall’s planes can occur if the difference in diffusivities changes sign with composition across the diffusion couple. Stability of the Kirkendall planes was studied by van Dal et al. [5, 8, 9]. They relate stability or instability to the sign, negative or positive, respectively, of the spatial divergence of the Kirkendall velocity at the position of the Kirkendall plane. The criterion established by Cornet and Calais [7] was verified experimentally for solid solutions of NiPd and FePd with single Kirkendall’s planes [5]. Double Kirkendall’s planes were found during interdiffusion in β-NiAl [10]. The term “stability” refers to whether markers placed on the plane of the initial discontinuity of a diffusion couple remain sharply concentrated on a plane or whether they become dispersed in the diffusion direction. Computation of the trajectories in solid solutions has been performed by Höglund and Ågren [11].The parabolic transform method will be presented and used in this paper for a finite-length diffusion couple. The stress generation due to diffusion in constrained geometries will not be discussed in the paper for many reasons. It was shown that in the case, analyzed here (all diffusion fluxes are linearly correlated), it is possible to obtain a symmetrical matrix of transport coefficients for the interdiffusion fluxes (even for a general case of a multicomponent system) and the expression for the entropy production. Thus, the drift generation is a stress-free deformation of the solid in one-dimensional diffusion couples. In complex geometries stresses play an important role and definitely are nonnegligible. It is possible to develop this approach, through more adequate formulations of the free energy functional [12, 13].The purpose of this paper is to describe the peculiarities of the deformation in three- and one-dimensional binary diffusion couples under various conditions. We compute the trajectories (position dependent on time) of a full set of markers, focusing on cases with multiple Kirkendall’s planes and sole Kirkendall’s planes that coincide with the moving initial contact interface. Numerical experiments are compared with Belova et al. [14].


## 2. Model

### 2.1. Diffusion and Material Velocity

To permit applications to nonideal solid solutions, we present expressions for the standard Darken approach in terms of gradients of potentials of the componentsμi in the three-dimensional space R3. The mass conservation laws have the form(1)∂ci∂t=-div⁡J~ifori=1,2,3,where ci=ci(t,x) and J~i=J~i(t,x) are the component concentrations and their fluxes measured in the laboratory RF. The most convenient RF in the case of a finite (bounded) diffusion couple considered in this work is, for example, the boundary of the diffusion couple.The central Darken postulate is defining the component diffusion fluxesJ~i(t,x) via a drift velocity field vD(t,x) that is measured in the laboratory reference frame. It is also calledconvection velocity [1],Darken velocity [15],lattice velocity [2], ormaterial velocity [16, 17]. It characterizes the solid solution motion and equals the velocity of a set of markers fixed in the material. If the material is crystalline, the markers move with the lattice, and the process occurs due to the existence of vacancy sources and sinks, for example, one- and two-dimensional defects. The diffusion fluxes J~i are called intrinsic (or lattice) fluxes and are related to those measured in the laboratory frame by the formula(2)J~i=Ji+civDfori=1,2,3.The overall molar volume Ω and mole fractions Ni(i=1,2,3) are defined by(3)Ni=ciΩ,Ω=1c=∑i=13ΩiNi,where c=∑i=13‍ci and Ωi is the ith component partial molar volume. From (3), we get(4)∑i=13Ωici=1.In the multicomponent solid solutions the following relations hold:(5)∑i=13Ωidci=0,∑i=13cidΩi=0,SdT+Ωdσ+∑i=13Nidμi=0.For isothermal processes, when all partial molar volumes are constant and stress effects are negligible, the relations (5) reduce to(6)∑i=13Ωidci=0,∑i=13Nidμi=0.Variable partial molar volumes in binary systems were treated by Sauer and Freise [18] and are not considered here.


### 2.2. Kirkendall Velocity and Displacement

Multiplying the conservation law for each species, (1), by matching partial molar volume, results in the volume continuity formula(7)∑i=13div⁡ΩiJ~i=0.The intrinsic fluxes, governed by the Nernst-Planck formula, are assumed to be proportional to the gradients of their chemical potentials (diffusion potentials in general [19]):(8)Ji=-ciBi∇μifori=1,2,3,where the intrinsic (or lattice) mobilities Bi may depend on the concentration. From relations (2), (3), (7), and (8), the Darken velocity field (vD Darken’s drift) is given by(9)div⁡vD=-∑i=13div⁡ΩiJi=∑i=13div⁡ΩiciBi∇μi.If the velocity field is not one-dimensional, additional equations might be required to determine the flow [16]. The mass conservation law can be written as(10)∂ci∂t=div⁡ciBi∇μi-civDfori=1,2,3,where the drift velocity fulfills (9). In the following, we will consider a ternary system (i=1,2,3) and constant partial molar volumes such that Ωi=Ω for i=1,2,3. Consequently, using the relations from (3), (9) and (10) become(11)div⁡vD=∑i=13div⁡NiBi∇μi,∂Ni∂t=div⁡NiBi∇μi-NivDfori=1,2,3.The Darken velocity field characterizes a deformation of the material during diffusion. If one considers diffusion in a crystalline material, the markers are positions in an imperfect lattice in which production and/or annihilation of vacancies occur at rates that are sufficient to maintain an equilibrium chemical potential of vacancies everywhere; that is, the local thermodynamical equilibrium is sustained in the material [20].The same equations apply to diffusion in fluids and amorphous solids, where the drift velocity describes the convective motion of the substance as indicated by the motion of small inert particles in a fluid as a result of different diffusion mobilities. The effective diffusion coefficient in gases, in the manner of Darken [2], describing the net convective motion, was developed already in 1899 in the kinetic theory of gases [20].


### 2.3. Diffusion Potentials and Diffusion Coefficients

The termdiffusion potential was introduced by Larché and Cahn [21]. They derived a potential as a sum of the elastic energy and chemical potential. In general, the diffusion potential is a sum all potentials that account for the storage of energy due to the addition of a diffusant [19]. The regular solutions can be treated by analogous formalism. In this work we focus on an ideal solid solution:(12)μi=μi0+RTlnNifori=1,2,3,where μi0 is the reference potential; R is the gas constant; T is the temperature. From relations (3), (8), and (12), the intrinsic diffusion fluxes Di are defined by(13)Ji=DiΩ∇Nifori=1,2,3.Equations (11) become(14)div⁡vD=∑i=13div⁡Di∇Ni,(15)∂Ni∂t=div⁡Di∇Ni-NivDfori=1,2,3.From (14) in the closed and bounded system the drift velocity is given by [22](16)vD=∑i=13Di∇Ni.Upon denoting u=N1, v=N2, and w=N3, from (15) and (16) we have(17)∂u∂t=D1∇2u-∇·uvD,∂v∂t=D2∇2v-∇·vvD,∂w∂t=D3∇2w-∇·wvDwith u+v+w=1. Values of u,v,w belong to the interval [0,1]. The initial conditions u(0,x)=u0(x), v(0,x)=v0(x), and w(0,x)=w0(x) are imposed for x=(x1,x2,x3)∈Q≔[-L1,L1]×[-L2,L2]×[-L3,L3]. Boundary Neumann conditions are given by(18)∂u∂n=∂v∂n=∂w∂n=0on the lateral boundary x=(x1,x2,x3)∈∂Q, t>0, where n=n(x) is the outward normal vector and ∂/∂n is the operator of normal derivative; see [23].


## 2.1. Diffusion and Material Velocity

To permit applications to nonideal solid solutions, we present expressions for the standard Darken approach in terms of gradients of potentials of the componentsμi in the three-dimensional space R3. The mass conservation laws have the form(1)∂ci∂t=-div⁡J~ifori=1,2,3,where ci=ci(t,x) and J~i=J~i(t,x) are the component concentrations and their fluxes measured in the laboratory RF. The most convenient RF in the case of a finite (bounded) diffusion couple considered in this work is, for example, the boundary of the diffusion couple.The central Darken postulate is defining the component diffusion fluxesJ~i(t,x) via a drift velocity field vD(t,x) that is measured in the laboratory reference frame. It is also calledconvection velocity [1],Darken velocity [15],lattice velocity [2], ormaterial velocity [16, 17]. It characterizes the solid solution motion and equals the velocity of a set of markers fixed in the material. If the material is crystalline, the markers move with the lattice, and the process occurs due to the existence of vacancy sources and sinks, for example, one- and two-dimensional defects. The diffusion fluxes J~i are called intrinsic (or lattice) fluxes and are related to those measured in the laboratory frame by the formula(2)J~i=Ji+civDfori=1,2,3.The overall molar volume Ω and mole fractions Ni(i=1,2,3) are defined by(3)Ni=ciΩ,Ω=1c=∑i=13ΩiNi,where c=∑i=13‍ci and Ωi is the ith component partial molar volume. From (3), we get(4)∑i=13Ωici=1.In the multicomponent solid solutions the following relations hold:(5)∑i=13Ωidci=0,∑i=13cidΩi=0,SdT+Ωdσ+∑i=13Nidμi=0.For isothermal processes, when all partial molar volumes are constant and stress effects are negligible, the relations (5) reduce to(6)∑i=13Ωidci=0,∑i=13Nidμi=0.Variable partial molar volumes in binary systems were treated by Sauer and Freise [18] and are not considered here.


## 2.2. Kirkendall Velocity and Displacement

Multiplying the conservation law for each species, (1), by matching partial molar volume, results in the volume continuity formula(7)∑i=13div⁡ΩiJ~i=0.The intrinsic fluxes, governed by the Nernst-Planck formula, are assumed to be proportional to the gradients of their chemical potentials (diffusion potentials in general [19]):(8)Ji=-ciBi∇μifori=1,2,3,where the intrinsic (or lattice) mobilities Bi may depend on the concentration. From relations (2), (3), (7), and (8), the Darken velocity field (vD Darken’s drift) is given by(9)div⁡vD=-∑i=13div⁡ΩiJi=∑i=13div⁡ΩiciBi∇μi.If the velocity field is not one-dimensional, additional equations might be required to determine the flow [16]. The mass conservation law can be written as(10)∂ci∂t=div⁡ciBi∇μi-civDfori=1,2,3,where the drift velocity fulfills (9). In the following, we will consider a ternary system (i=1,2,3) and constant partial molar volumes such that Ωi=Ω for i=1,2,3. Consequently, using the relations from (3), (9) and (10) become(11)div⁡vD=∑i=13div⁡NiBi∇μi,∂Ni∂t=div⁡NiBi∇μi-NivDfori=1,2,3.The Darken velocity field characterizes a deformation of the material during diffusion. If one considers diffusion in a crystalline material, the markers are positions in an imperfect lattice in which production and/or annihilation of vacancies occur at rates that are sufficient to maintain an equilibrium chemical potential of vacancies everywhere; that is, the local thermodynamical equilibrium is sustained in the material [20].The same equations apply to diffusion in fluids and amorphous solids, where the drift velocity describes the convective motion of the substance as indicated by the motion of small inert particles in a fluid as a result of different diffusion mobilities. The effective diffusion coefficient in gases, in the manner of Darken [2], describing the net convective motion, was developed already in 1899 in the kinetic theory of gases [20].


## 2.3. Diffusion Potentials and Diffusion Coefficients

The termdiffusion potential was introduced by Larché and Cahn [21]. They derived a potential as a sum of the elastic energy and chemical potential. In general, the diffusion potential is a sum all potentials that account for the storage of energy due to the addition of a diffusant [19]. The regular solutions can be treated by analogous formalism. In this work we focus on an ideal solid solution:(12)μi=μi0+RTlnNifori=1,2,3,where μi0 is the reference potential; R is the gas constant; T is the temperature. From relations (3), (8), and (12), the intrinsic diffusion fluxes Di are defined by(13)Ji=DiΩ∇Nifori=1,2,3.Equations (11) become(14)div⁡vD=∑i=13div⁡Di∇Ni,(15)∂Ni∂t=div⁡Di∇Ni-NivDfori=1,2,3.From (14) in the closed and bounded system the drift velocity is given by [22](16)vD=∑i=13Di∇Ni.Upon denoting u=N1, v=N2, and w=N3, from (15) and (16) we have(17)∂u∂t=D1∇2u-∇·uvD,∂v∂t=D2∇2v-∇·vvD,∂w∂t=D3∇2w-∇·wvDwith u+v+w=1. Values of u,v,w belong to the interval [0,1]. The initial conditions u(0,x)=u0(x), v(0,x)=v0(x), and w(0,x)=w0(x) are imposed for x=(x1,x2,x3)∈Q≔[-L1,L1]×[-L2,L2]×[-L3,L3]. Boundary Neumann conditions are given by(18)∂u∂n=∂v∂n=∂w∂n=0on the lateral boundary x=(x1,x2,x3)∈∂Q, t>0, where n=n(x) is the outward normal vector and ∂/∂n is the operator of normal derivative; see [23].


## 3. Reduction to Parabolic Equations inR3

In three-component models with three distinct diffusion coefficientsD1>0,D2>0,D3>0 the third variable w can be expressed in terms of the other variables u,v, w=1-u-v. Denote(19)D1′≔D1-D3,D2′≔D2-D3.Assume that D1′,D2′>0, which means that D3 is the smallest diffusion coefficient (it is convenient to assume that D1>D2>D3). Thus, the ternary system (17) reduces to the following:(20)∂u∂t=D1∇2u-uD1′∇2u+D2′∇2v-∇u·D1′∇u+D2′∇v,∂v∂t=D2∇2u-vD1′∇2u+D2′∇2v-∇v·D1′∇u+D2′∇vwith ∂u/∂n=0, ∂v/∂n=0 on the lateral boundary x∈∂Q, t>0 and u(0,x)=u0(x) and v(0,x)=v0(x) for x∈Q. We reduce this to a parabolic system by virtue of a change of dependent variables. Take(21)G≔Γ∘u,v,Gt,x=Γut,x,vt,xand denote(22)Γ1=∂Γ∂u,Γ2=∂Γ∂v;then(23)∂G∂t=Γ1∂u∂t+Γ2∂v∂t=Γ1∇2uD1-uD1′-uD2′∇2v+Γ2-vD1′∇2u+D2-vD2′∇2v-Γ1∇u+Γ2∇v·D1′∇u+D2′∇v.Similarly, one calculates(24)∂G∂xj=Γ1∂u∂xj+Γ2∂v∂xjfor j=1,2,3 and(25)∂2G∂xjl2=Γ1∂2u∂xjl2+Γ2∂2v∂xjl2+Γ11∂u∂xj∂u∂xl+Γ12∂u∂xj∂v∂xl+Γ21∂u∂xl∂v∂xj+Γ22∂v∂xj∂v∂xlfor j,l=1,2,3 where Γ11,Γ12,Γ21, and Γ22 are the second-order derivatives of Γ in u,v. Consider the parabolic operator(26)∂G∂t-λ∇2G=Γ1D1-uD1′-λ-Γ2vD1′∇2u+-Γ1uD2′+Γ2D2-vD2′-λ∇2v+“lower  order  terms”,where “lower order terms” will be specified later. The second-order terms vanish, provided that λ is an eigenvalue(27)D1-uD1′-λ-vD1′-uD2′D2-vD2′-λ=0,and (Γ1,Γ2) is an eigenvector. The characteristic polynomial(28)λ2-λD1+D2-uD1′-vD2′+D1D2-D1vD2′-D2uD1′has positive and distinct roots, because u+v+w=1; thus(29)1>u1-D3D1+v1-D3D2,and Δ0(u,v)>0, where(30)Δ0u,v=D1-D2-uD1′+vD2′2+4uD1′vD2′.The eigenvalues are equal to(31)λ1,2u,v=12D1+D2-uD1′-vD2′±12Δ0u,v.Their eigenvectors fulfill the algebraic equation(32)Γ1D1-uD1′-λ1,2u,v-Γ2vD1′=0.Hence, we have the hyperbolic PDE(33)∂Γ/∂ξηD1′-∂Γ/∂ηD1-ξD1′-λ1,2ξ,η=0.Its characteristics(34)ηD1′dξ+D1-ξD1′-λ1,2ξ,ηdη=0have solutions Γ(ξ,η)=const. Thus, we get two families of curves Γ(1)=const. for λ1 and Γ(2)=const. for λ2. Introduce new dependent variables(35)G1=Γ1u,v,G2=Γ2u,vand the inverse functions(36)u=UG1,G2,v=VG1,G2,wherever they exist. We find their gradients from the system(37)∇G1=∂Γ1∂u∇u+∂Γ1∂v∇v,∇G2=∂Γ2∂u∇u+∂Γ2∂v∇v,which can be inverted as(38)∇u=U1∇G1+U2∇G2,∇v=V1∇G1+V2∇G2,where Ui=UiG1,G2 and Vi=ViG1,G2 for i=1,2. The ultimate parabolic PDE system looks as follows:(39)∂Gi∂t-λiΔGi=-λiΓ11iU1∇G1+U2∇G22+2Γ12iU1∇G1+U2∇G2·V1∇G1+V2∇G2+Γ22iV1∇G1+V2∇G22for i=1,2, where λi=λi(U,V). The right-hand side of (39) consists of the “lower order terms” from (26).


## 4. Trajectories inR3

An arbitrary smooth trajectoryz=z(t) is governed by the ordinary differential problem(40)ddtz=gt,z,z0=x∈Q,where g depends on u,v,w and ∇u,∇v,∇w; confer [23]. Recall that x=(x1,x2,x3). It is obvious that trajectories of material points are governed by ordinary differential equations. The only question is as follows: how to derive these ODEs? We draw the reader’s attention to the microscopic setting based on stochastic diffusion equations in [24].


### 4.1. Material Trajectories

In Darken’s method the functiong is equal to the Darken (material) velocity vD:(41)vD=D1∇u+D2∇v+D3∇w.Suppose that the marker at t=0 is placed at a surface F0(x)=0 for x=(x1,x2,x3). We want to trace its trajectory for t≥0 which constitutes planes F(t,x)=0 for x=(x1,x2,x3). Since the characteristics are governed by (40), F fulfills the hyperbolic initial-value problem(42)∂F∂t+∇F·g=0,F0,x=F0x,because(43)ddtFt,zt=∂F∂t+∇F·dzdt=0.In particular, Darken’s trajectory fulfills the equation(44)dzdt=vDt,z,and Kirkendall’s plane satisfies the PDE(45)∂F∂t+∇F·vD=0,F0,x=F0x,where x=(x1,x2,x3); see [23].The Kirkendall plane given by (45) with characteristics coinciding with Darken’s trajectories (44) fulfills the following elegant property.Proposition 1.
For any regular regionΩt with the boundary ∂Ωt defined by F(t,x)=0 where F is the solution of the PDE (45), or equivalently being constituted by the characteristics (t,z(t)) obeying (44) starting from the initial plane F0(x)=0, the total mass of each component u,v,w plus its overall inflow through the boundary is constant:(46)∫Ωtudx-∫0t∫∂ΩτD1∂u∂ndSdτ=∫Ω0u0dx.Analogous identities are valid for v and w.Proof.
Denote byLH=LH(t) the left-hand side of (46). Then, we have(47)LH=∫Ωtudx-∫0t∫ΩτD1∇2udxdτ.
We differentiateLH as follows:(48)LH′=∫∂Ωtn·uz′dS+∫Ωtutdx-∫ΩtD1∇2udx=∫∂Ωtn·uz′dS-∫Ωt∇·uvDdx=∫∂Ωtn·uz′dS-∫∂Ωtn·uvDdS=0.Since LH′=0, we have LH≡const.=∫Ω0‍u0dx. This proves the assertion.This property proves that in ternary and higher solid solutions the Kirkendall plane allows finding the Matano plane that is identical for all components.


### 4.2. Displacement Field

Since the velocity fieldv is described in terms of laboratory coordinates, the displacement field u is a function of the laboratory coordinates and time as well. The relation between the reference and laboratory coordinates is defined by a mapping Z=Z(t,z) which gives a reference (initial) position of the marker that is currently at position z; the inverse mapping is z=z(t,Z) providing the current position z of the marker initially placed at position Z. The displacement field u satisfies the relation(49)ut,z=z-Zt,z.That is, u(t,z) is the difference between the current and initial positions of the marker which is now at position z in the laboratory frame. The displacement and velocity fields satisfy the equation(50)vt,z=∂ut,z∂t+vt,z∂ut,z∂z.


### 4.3. Component Trajectories

Consider any regular subdomainΩ0⊂Q in R3 with a smooth boundary F0(x)=0. If we are interested in the evolution of u starting from Ω0 then we obtain a family of subdomains Ωt described by the equations F(t,x)=0. The mass balance leads to the integral identity(51)ddt∫Ωtudx=0.The integration by parts leads to the formula(52)∫∂Ωtu∂F∂t+n·D1∇u-uvDdS=0,where n=n(x) is the unit normal vector, proportional to ∇F, and dS is the surface element. A similar relation holds for v and w. We can get the evolution of a region which contains some amount of u,v:(53)ddt∫Ωtu+vdx=0,∫∂Ωtu+v∂F∂t+n·D1∇u+D2∇v-u+vvDdS=0.(a) Balance of Mass in 1D. Consider the one-dimensional case with the central position separating substances at x=0. Thus, F0(x)=x∈R. Suppose that the separating position z=z(t) moves according to the mass-balance law(54)∫-Lzu+vdx=const.which means that the total mass of u+v is preserved on the left of the separating position z=z(t). If we differentiate the identity (54) with respect to t we obtain(55)z′u+v+∫-Lzut+vtdx=0.Taking into account the governing equations, we get(56)z′u+v+∫-LzD1uxx+D2vxx-u+vvDxdx=0;hence,(57)z′u+v+D1ux+D2vx-u+vvD=0.In this case, we arrive at the equation of the separating position trajectory(58)z′=vD-D1ux+D2vxu+vor equivalently(59)z′1-w+wvD=D3wx.The respective plane is a solution of the PDE(60)∂F∂t+∂F∂xvD-D1ux+D2vxu+v=0,F0,x=x.(b) Tracer Trajectory. Let us return to the trajectory of w which starts from a spike-shaped droplet. Suppose that we measure the middle position z3=z3(t) of this marker; hence, we have the equation(61)∫-Lz3wdx=12∫-LLwdx=const.Alike (54), we arrive at the differential equation(62)z3′=vD-D3wxw,z30=0.The ideal marker means a fictitious observer (indicator) within solid moving with the local drift velocity [2] and not altering the diffusive mass flow. In practical realizations nonideal markers affect diffusion for the reason that they are often metal wires or submicron inclusions of a stable oxide embedded in the diffusion zone. The method presented by Belova et al. [14] and analyzed in the present paper has an advantage of a negligible impact on the diffusive mass transport (when concentrations of marker atoms are low). An obvious drawback of both methods (Belova’s and ours) is diffusion of the marker itself. This process can be to some extent delayed by an introduction of marker atoms whose diffusivity in the material is low; see Figure 1. Observe that the trajectory of z3 is very close to the ideal marker, since the maximum of w(t,·) is attained at a point x such that wx(t,x)=0. We draw a useful conclusion: z3′≈vD. We can also derive similar equations to (62) for u and v. Having in mind a severe singularity at x=0, we can analyse the front movement of u and v as follows. Suppose that u is initially concentrated on the left x∈[-L,0]. We study an evolution z1=z1(t) defined by(63)∫-Lz1udx=∫-L-εu0dx=const.with some small ε>0 (to avoid the singularity). Hence, we get the differential equation(64)z1′=vD-D1uxu,z10=-ε.Since u(t,·) is decreasing, we have ux≤0 and z1 moves to the right. In an analogous way z2=z2(t) initially placed at x=+ε moves to the left according to the following:(65)z2′=vD-D2vxv,z20=+ε.Figure 1
Distribution att=0.01 for D1=1, D2=0.7, and D3=0.01. The initial distribution of the component w is a very thin nanolayer around x=0 where w=0.2.


## 4.1. Material Trajectories

In Darken’s method the functiong is equal to the Darken (material) velocity vD:(41)vD=D1∇u+D2∇v+D3∇w.Suppose that the marker at t=0 is placed at a surface F0(x)=0 for x=(x1,x2,x3). We want to trace its trajectory for t≥0 which constitutes planes F(t,x)=0 for x=(x1,x2,x3). Since the characteristics are governed by (40), F fulfills the hyperbolic initial-value problem(42)∂F∂t+∇F·g=0,F0,x=F0x,because(43)ddtFt,zt=∂F∂t+∇F·dzdt=0.In particular, Darken’s trajectory fulfills the equation(44)dzdt=vDt,z,and Kirkendall’s plane satisfies the PDE(45)∂F∂t+∇F·vD=0,F0,x=F0x,where x=(x1,x2,x3); see [23].The Kirkendall plane given by (45) with characteristics coinciding with Darken’s trajectories (44) fulfills the following elegant property.Proposition 1.
For any regular regionΩt with the boundary ∂Ωt defined by F(t,x)=0 where F is the solution of the PDE (45), or equivalently being constituted by the characteristics (t,z(t)) obeying (44) starting from the initial plane F0(x)=0, the total mass of each component u,v,w plus its overall inflow through the boundary is constant:(46)∫Ωtudx-∫0t∫∂ΩτD1∂u∂ndSdτ=∫Ω0u0dx.Analogous identities are valid for v and w.Proof.
Denote byLH=LH(t) the left-hand side of (46). Then, we have(47)LH=∫Ωtudx-∫0t∫ΩτD1∇2udxdτ.
We differentiateLH as follows:(48)LH′=∫∂Ωtn·uz′dS+∫Ωtutdx-∫ΩtD1∇2udx=∫∂Ωtn·uz′dS-∫Ωt∇·uvDdx=∫∂Ωtn·uz′dS-∫∂Ωtn·uvDdS=0.Since LH′=0, we have LH≡const.=∫Ω0‍u0dx. This proves the assertion.This property proves that in ternary and higher solid solutions the Kirkendall plane allows finding the Matano plane that is identical for all components.


## 4.2. Displacement Field

Since the velocity fieldv is described in terms of laboratory coordinates, the displacement field u is a function of the laboratory coordinates and time as well. The relation between the reference and laboratory coordinates is defined by a mapping Z=Z(t,z) which gives a reference (initial) position of the marker that is currently at position z; the inverse mapping is z=z(t,Z) providing the current position z of the marker initially placed at position Z. The displacement field u satisfies the relation(49)ut,z=z-Zt,z.That is, u(t,z) is the difference between the current and initial positions of the marker which is now at position z in the laboratory frame. The displacement and velocity fields satisfy the equation(50)vt,z=∂ut,z∂t+vt,z∂ut,z∂z.


## 4.3. Component Trajectories

Consider any regular subdomainΩ0⊂Q in R3 with a smooth boundary F0(x)=0. If we are interested in the evolution of u starting from Ω0 then we obtain a family of subdomains Ωt described by the equations F(t,x)=0. The mass balance leads to the integral identity(51)ddt∫Ωtudx=0.The integration by parts leads to the formula(52)∫∂Ωtu∂F∂t+n·D1∇u-uvDdS=0,where n=n(x) is the unit normal vector, proportional to ∇F, and dS is the surface element. A similar relation holds for v and w. We can get the evolution of a region which contains some amount of u,v:(53)ddt∫Ωtu+vdx=0,∫∂Ωtu+v∂F∂t+n·D1∇u+D2∇v-u+vvDdS=0.(a) Balance of Mass in 1D. Consider the one-dimensional case with the central position separating substances at x=0. Thus, F0(x)=x∈R. Suppose that the separating position z=z(t) moves according to the mass-balance law(54)∫-Lzu+vdx=const.which means that the total mass of u+v is preserved on the left of the separating position z=z(t). If we differentiate the identity (54) with respect to t we obtain(55)z′u+v+∫-Lzut+vtdx=0.Taking into account the governing equations, we get(56)z′u+v+∫-LzD1uxx+D2vxx-u+vvDxdx=0;hence,(57)z′u+v+D1ux+D2vx-u+vvD=0.In this case, we arrive at the equation of the separating position trajectory(58)z′=vD-D1ux+D2vxu+vor equivalently(59)z′1-w+wvD=D3wx.The respective plane is a solution of the PDE(60)∂F∂t+∂F∂xvD-D1ux+D2vxu+v=0,F0,x=x.(b) Tracer Trajectory. Let us return to the trajectory of w which starts from a spike-shaped droplet. Suppose that we measure the middle position z3=z3(t) of this marker; hence, we have the equation(61)∫-Lz3wdx=12∫-LLwdx=const.Alike (54), we arrive at the differential equation(62)z3′=vD-D3wxw,z30=0.The ideal marker means a fictitious observer (indicator) within solid moving with the local drift velocity [2] and not altering the diffusive mass flow. In practical realizations nonideal markers affect diffusion for the reason that they are often metal wires or submicron inclusions of a stable oxide embedded in the diffusion zone. The method presented by Belova et al. [14] and analyzed in the present paper has an advantage of a negligible impact on the diffusive mass transport (when concentrations of marker atoms are low). An obvious drawback of both methods (Belova’s and ours) is diffusion of the marker itself. This process can be to some extent delayed by an introduction of marker atoms whose diffusivity in the material is low; see Figure 1. Observe that the trajectory of z3 is very close to the ideal marker, since the maximum of w(t,·) is attained at a point x such that wx(t,x)=0. We draw a useful conclusion: z3′≈vD. We can also derive similar equations to (62) for u and v. Having in mind a severe singularity at x=0, we can analyse the front movement of u and v as follows. Suppose that u is initially concentrated on the left x∈[-L,0]. We study an evolution z1=z1(t) defined by(63)∫-Lz1udx=∫-L-εu0dx=const.with some small ε>0 (to avoid the singularity). Hence, we get the differential equation(64)z1′=vD-D1uxu,z10=-ε.Since u(t,·) is decreasing, we have ux≤0 and z1 moves to the right. In an analogous way z2=z2(t) initially placed at x=+ε moves to the left according to the following:(65)z2′=vD-D2vxv,z20=+ε.Figure 1
Distribution att=0.01 for D1=1, D2=0.7, and D3=0.01. The initial distribution of the component w is a very thin nanolayer around x=0 where w=0.2.


## 5. Results

In the case of a binary, one-dimensional system by using the Gibbs-Duhem equation, (5), one can define the interdiffusion mobility B~:(66)B~=N1N2N1B2+N2B1,and diffusion equation has the form [17](67)∂c2∂t=div⁡B~Ω1Ω2∇μ2-μ1.Consider the one-dimensional case with Ωt=[-L,z(t)], where z=z(t) satisfies (44). We formulate a simple consequence of Proposition 1 for the Kirkendall plane.Proposition 2.
Ifz=z(t) satisfies (44) with z(0)=0, the total mass of each component u,v,w plus its total diffusive flux through the boundary of the interval [-L,z(t)] is constant:(68)∫-Lzudx-∫0tD1∫-Lzτuxxdxdτ=∫-L0u0dx.Analogous identities are valid for v and w.Proof.
Denote byLH=LH(t) the left-hand side of (68). Then, we have(69)LH≔∫-Lzudx-∫0tD1∫-Lzτuxxdxdτ.If we differentiate this expression with respect to t, we get(70)z′u+∫-Lzut-D1∫-Lzuxxdx=z′u-∫-LzuvDxdx=0.This means that LH′=0; hence, LH≡const. and the assertion follows.Remark 3.
We draw the reader’s attention to the conclusion of Proposition2: the same trajectory, given by (44) with z(0)=0, is feasible for each component u,v,w. This fact allows for effective computations of Kirkendall’s planes.Corollary 4 (gradients on the Kirkendall plane).
Under the assumptions of Proposition2 we have the following gradient representation on the Kirkendall plane z=z(t):(71)D1ux=vDu+∫-Lzutdx.Analogous representations have D2vx and D3wx.Proof.
If we differentiate (68), we get(72)uz′+∫-Lzutdx-D1∫-Lzuxxdx=0.Because z′=vD and the last term is equal to -D1ux(z), we obtain the assertion.


### 5.1. Up-Hill Diffusion

It was shown in [25] that uphill diffusion of w is indicated by the condition(73)∂w∂t>0.This inequality can be locally achieved by means of suitable distributions of u and v, provided that D1≠D2.We return to this discussion because it is possible to simplify the reasoning. Taking into account (17) and u+v=1-w we get(74)-∂w∂t=D1∇2u+D2∇2v-∇·1-wvD.Hence the uphill requirement (73) becomes(75)D1-1-wD1-D3∇2u+D2-1-wD1-D3∇2v+∇w·vD<0.Consider the case w≡const.∈(0,1); see Figure 4. Then, one gets ∇w≡0 and u+v≡const.; that is, ∇u+∇v=0. The uphill condition (75) is as follows:(76)D1-1-wD1-D3∇2u-D2-1-wD1-D3∇2u<0,or equivalently(77)D1-D2∇2u<0.If D1>D2 then any region of concavity of u yields ∇2u<0; hence, the concentration w increases as it is illustrated in Figure 4.


## 5.1. Up-Hill Diffusion

It was shown in [25] that uphill diffusion of w is indicated by the condition(73)∂w∂t>0.This inequality can be locally achieved by means of suitable distributions of u and v, provided that D1≠D2.We return to this discussion because it is possible to simplify the reasoning. Taking into account (17) and u+v=1-w we get(74)-∂w∂t=D1∇2u+D2∇2v-∇·1-wvD.Hence the uphill requirement (73) becomes(75)D1-1-wD1-D3∇2u+D2-1-wD1-D3∇2v+∇w·vD<0.Consider the case w≡const.∈(0,1); see Figure 4. Then, one gets ∇w≡0 and u+v≡const.; that is, ∇u+∇v=0. The uphill condition (75) is as follows:(76)D1-1-wD1-D3∇2u-D2-1-wD1-D3∇2u<0,or equivalently(77)D1-D2∇2u<0.If D1>D2 then any region of concavity of u yields ∇2u<0; hence, the concentration w increases as it is illustrated in Figure 4.


## 6. Numerical Results and Conclusions

### 6.1. Numerical Experiments

We apply standard finite difference approximations to system (17) based on classical theories [26]. Their stability and convergence can be justified by the same parabolic transform which leads to the parabolic system. In our numerical experiments we take the interval Q=[-1,1]. We performed computations by means of explicit FDMs with many discretization parameters. The results presented here are done for time step ht=1/500000 and space step hx=1/500. The initial condition for Figures 1 and 4 are the same as in [27]; the initial condition for Figure 5 are taken from [14] in a graphical mode. The nonlinear parts are approximated in the same way as in [27].(i)
Figure1 presents the cases D1=1, D2=0.7, D3=0.01, and the marker at t=0 (droplet). The initial conditions are as follows: the droplet w(0,x) is equal to 0 outside the interval (-0.01,0.01), and it is a polynomial on (-0.01,0.01), namely, 20000000·x4-4000·x2+0.2. The functions u(0,x), v(0,x) satisfy the following properties: v(0,x)=u(0,-x) and u(0,x)+v(0,x)=1-w(0,x); u(0,x) is nonincreasing and v(0,x) is nondecreasing.(ii)
Figures2 and 3 show characteristics for λ1, λ2 and the regions mapped by the transformation (u,v)↦(Γ1,Γ2) for D1=1, D2=0.7, and D3=0.0001. In order to show the range of characteristics we take a coarse mesh with N=10, h=0.1 (Figure 2), and a dense mesh N=50, h=0.02 (Figure 3). Computations are made on the dense mesh, whereas the figure shows only characteristics which start from the coarse mesh.(iii)
Figure4 shows the cases D1=1, D2=0.1, and D3=0.01 and a uniformly distributed substance w at t=0. One can observe uphill effects. The medium w(0,x) is uniformly distributed, w(0,x)≡0. The functions u(0,x),v(0,x) satisfy the following properties: v(0,x)=u(0,-x) and u(0,x)+v(0,x)=1-w(0,x); u(0,x) is nonincreasing, v(0,x) is nondecreasing, and u(0,x), v(0,x) are polynomials ±200000·x3∓60·x+0.4 on (-0.01,0.01).(iv)
Figure5 demonstrates the agreement with Belova’s results from [14]. Our computations are made with the same data as [14].Figure 2
Characteristics forλ1 and λ2 with N=10 and D1=1, D2=0.7, and D3=0.0001.
(a)
(b)Figure 3
Regions mapped by the parabolic transformation.N=50, D1=1, D2=0.7, and D3=0.0001.
(a)
(b)Figure 4
Uphill diffusion, distributions ofu,v,w at t=0 and t=0.01. Component w is uniformly distributed for t=0, w≡0.2.Figure 5
Uphill diffusion, an experiment with data from Belova [14].


### 6.2. Conclusions

(1)
In the literature one-dimensional Kirkendall’s planes are treated. We show a three-dimensional evolution of the Kirkendall plane; confer (45). Our considerations of Darken’s trajectories passing through all points of the domain lead to the most convenient definition of this object; namely, these trajectories are characteristic curves of a transport equation that uniquely defines a family of planes.
(2)
The proposed parabolic transform is justified for smooth initial data until they remain in[0,1]. In fact these densities become more regular as time goes by. Moreover, they obey a kind of maximum principle which implies the mass conservation and [0,1]-invariance.
(3)
Observe that the mass balance matching extremum detection leads to a feasible criterion of experimental finding the Kirkendall plane by introducing an isotope thin layer on an initial contact interface; compare (62).
(4)
We propose a new method of overcoming the problems caused by an initial singularity at the contact surface of the diffusing substances; see (64) and (65).
(5)
Our discussion concerning uphill diffusion is concluded by an intuitive observation: ifw is uniformly distributed and D1>D2, then uphill diffusion of w is driven by any concavity of u; see (75). However, if the marker is a thin layer then an uphill effect is peculiar, because a fast dissipation of this substance is expected rather than any increase of this spike. Is it possible to provoke an uphill effect for a thin layer of the marker w? Contrary to the above expected behavior, one can find such initial distributions of u,v that an uphill phenomenon occurs. Suppose that w is an isotope and its maximum can be measured; then we have ∇w=0 and ∇2w≤0 at this maximal point. Hence, uphill condition (75) is fulfilled provided that we have at this maximal point the inequality(78)D1-D2w∇2u<D2-1-wD2-D3∇2w.Since ∇2w≤0 at this maximum, then uphill diffusion can be generated by concavity of u at this point under the mild requirement D3+w(D2-D3)≥0, which holds for the case D1>D2>D3.
(6)
There is an alternative view of the gradient∇Γ in our parabolic transform; namely, one can look for a nondegenerate function Γ=Γ(u,v) whose gradient ∇Γ=(Γ1,Γ2) is proportional to an eigenvector E=E(u,v)=(E1,E2) corresponding to λ=λ(u,v). This reads(79)∇Γ=ρE,ρ=ρu,v.Then ρ>0 satisfies the equation(80)ρuE2+ρE2,u=ρvE1+ρE1,v.Characteristics of this equation fulfill the systems u′=E2, v′=-E1, and Γ can be calculated explicitly as(81)Γu,v=∫0,0u,vρE1du+ρE2dv.This generalizes the diagonalization method by Krishtal et al. [28]. Their transform of diffusivities was linear, whereas our transform is nonlinear.
(7)
Observe that the relationD1>D2>D3 leads to finding the eigenvalues λ1,λ2 in the intervals (D2,D1), (D3,D2), respectively. One can prove an analogous statement for multicomponent systems with n>4; however, it is not obvious whether suitable parabolic transforms can be derived, because the resulting PDE system might well be overdetermined.


## 6.1. Numerical Experiments

We apply standard finite difference approximations to system (17) based on classical theories [26]. Their stability and convergence can be justified by the same parabolic transform which leads to the parabolic system. In our numerical experiments we take the interval Q=[-1,1]. We performed computations by means of explicit FDMs with many discretization parameters. The results presented here are done for time step ht=1/500000 and space step hx=1/500. The initial condition for Figures 1 and 4 are the same as in [27]; the initial condition for Figure 5 are taken from [14] in a graphical mode. The nonlinear parts are approximated in the same way as in [27].(i)
Figure1 presents the cases D1=1, D2=0.7, D3=0.01, and the marker at t=0 (droplet). The initial conditions are as follows: the droplet w(0,x) is equal to 0 outside the interval (-0.01,0.01), and it is a polynomial on (-0.01,0.01), namely, 20000000·x4-4000·x2+0.2. The functions u(0,x), v(0,x) satisfy the following properties: v(0,x)=u(0,-x) and u(0,x)+v(0,x)=1-w(0,x); u(0,x) is nonincreasing and v(0,x) is nondecreasing.(ii)
Figures2 and 3 show characteristics for λ1, λ2 and the regions mapped by the transformation (u,v)↦(Γ1,Γ2) for D1=1, D2=0.7, and D3=0.0001. In order to show the range of characteristics we take a coarse mesh with N=10, h=0.1 (Figure 2), and a dense mesh N=50, h=0.02 (Figure 3). Computations are made on the dense mesh, whereas the figure shows only characteristics which start from the coarse mesh.(iii)
Figure4 shows the cases D1=1, D2=0.1, and D3=0.01 and a uniformly distributed substance w at t=0. One can observe uphill effects. The medium w(0,x) is uniformly distributed, w(0,x)≡0. The functions u(0,x),v(0,x) satisfy the following properties: v(0,x)=u(0,-x) and u(0,x)+v(0,x)=1-w(0,x); u(0,x) is nonincreasing, v(0,x) is nondecreasing, and u(0,x), v(0,x) are polynomials ±200000·x3∓60·x+0.4 on (-0.01,0.01).(iv)
Figure5 demonstrates the agreement with Belova’s results from [14]. Our computations are made with the same data as [14].Figure 2
Characteristics forλ1 and λ2 with N=10 and D1=1, D2=0.7, and D3=0.0001.
(a)
(b)Figure 3
Regions mapped by the parabolic transformation.N=50, D1=1, D2=0.7, and D3=0.0001.
(a)
(b)Figure 4
Uphill diffusion, distributions ofu,v,w at t=0 and t=0.01. Component w is uniformly distributed for t=0, w≡0.2.Figure 5
Uphill diffusion, an experiment with data from Belova [14].


## 6.2. Conclusions

(1)
In the literature one-dimensional Kirkendall’s planes are treated. We show a three-dimensional evolution of the Kirkendall plane; confer (45). Our considerations of Darken’s trajectories passing through all points of the domain lead to the most convenient definition of this object; namely, these trajectories are characteristic curves of a transport equation that uniquely defines a family of planes.
(2)
The proposed parabolic transform is justified for smooth initial data until they remain in[0,1]. In fact these densities become more regular as time goes by. Moreover, they obey a kind of maximum principle which implies the mass conservation and [0,1]-invariance.
(3)
Observe that the mass balance matching extremum detection leads to a feasible criterion of experimental finding the Kirkendall plane by introducing an isotope thin layer on an initial contact interface; compare (62).
(4)
We propose a new method of overcoming the problems caused by an initial singularity at the contact surface of the diffusing substances; see (64) and (65).
(5)
Our discussion concerning uphill diffusion is concluded by an intuitive observation: ifw is uniformly distributed and D1>D2, then uphill diffusion of w is driven by any concavity of u; see (75). However, if the marker is a thin layer then an uphill effect is peculiar, because a fast dissipation of this substance is expected rather than any increase of this spike. Is it possible to provoke an uphill effect for a thin layer of the marker w? Contrary to the above expected behavior, one can find such initial distributions of u,v that an uphill phenomenon occurs. Suppose that w is an isotope and its maximum can be measured; then we have ∇w=0 and ∇2w≤0 at this maximal point. Hence, uphill condition (75) is fulfilled provided that we have at this maximal point the inequality(78)D1-D2w∇2u<D2-1-wD2-D3∇2w.Since ∇2w≤0 at this maximum, then uphill diffusion can be generated by concavity of u at this point under the mild requirement D3+w(D2-D3)≥0, which holds for the case D1>D2>D3.
(6)
There is an alternative view of the gradient∇Γ in our parabolic transform; namely, one can look for a nondegenerate function Γ=Γ(u,v) whose gradient ∇Γ=(Γ1,Γ2) is proportional to an eigenvector E=E(u,v)=(E1,E2) corresponding to λ=λ(u,v). This reads(79)∇Γ=ρE,ρ=ρu,v.Then ρ>0 satisfies the equation(80)ρuE2+ρE2,u=ρvE1+ρE1,v.Characteristics of this equation fulfill the systems u′=E2, v′=-E1, and Γ can be calculated explicitly as(81)Γu,v=∫0,0u,vρE1du+ρE2dv.This generalizes the diagonalization method by Krishtal et al. [28]. Their transform of diffusivities was linear, whereas our transform is nonlinear.
(7)
Observe that the relationD1>D2>D3 leads to finding the eigenvalues λ1,λ2 in the intervals (D2,D1), (D3,D2), respectively. One can prove an analogous statement for multicomponent systems with n>4; however, it is not obvious whether suitable parabolic transforms can be derived, because the resulting PDE system might well be overdetermined.



---

*Source: 650452-2015-06-02.xml*