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

The problem of Kirkendall’s trajectories in finite, threeand one-dimensional ternary diffusion couples is studied. By means of the parabolic transformationmethod, 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 bothmethods 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.


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 onedimensional 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 be (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 -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 (0) for all  ≥ 0. Kirkendall's plane(s) problem from the origin of coordinates at  = 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 plane (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].

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   in the three-dimensional space R 3 .The mass conservation laws have the form where   =   (, ) and J = J (, ) 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 fluxes J (, ) via a drift velocity field V  (, ) that is measured in the laboratory reference frame.It is also called convection velocity [1], Darken velocity [15], lattice velocity [2], or material 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 twodimensional defects.The diffusion fluxes J are called intrinsic (or lattice) fluxes and are related to those measured in the laboratory frame by the formula J =   +   V  for  = 1, 2, 3. ( The overall molar volume Ω and mole fractions   ( = 1, 2, 3) are defined by where  = ∑ 3 =1   and Ω  is the th component partial molar volume.From (3), we get In the multicomponent solid solutions the following relations hold: Mathematical Problems in Engineering 3 For isothermal processes, when all partial molar volumes are constant and stress effects are negligible, the relations (5) reduce to Variable partial molar volumes in binary systems were treated by Sauer and Freise [18] and are not considered here.

Kirkendall Velocity and Displacement.
Multiplying the conservation law for each species, (1), by matching partial molar volume, results in the volume continuity formula 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]): where the intrinsic (or lattice) mobilities   may depend on the concentration.From relations (2), ( 3), (7), and ( 8), the Darken velocity field (V  Darken's drift) is given by 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 where the drift velocity fulfills (9).In the following, we will consider a ternary system ( = 1, 2, 3) and constant partial molar volumes such that Ω  = Ω for  = 1, 2, 3. Consequently, using the relations from (3), ( 9) and ( 10) become 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].

Diffusion Potentials and Diffusion
Coefficients.The term diffusion 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: where  0  is the reference potential;  is the gas constant;  is the temperature.From relations (3), (8), and ( 12), the intrinsic diffusion fluxes   are defined by Equations ( 11) become From ( 14) in the closed and bounded system the drift velocity is given by [22] Upon denoting  =  1 , V =  2 , and  =  3 , from ( 15) and ( 16) we have with on the lateral boundary  = ( 1 ,  2 ,  3 ) ∈ ,  > 0, where  = () is the outward normal vector and / is the operator of normal derivative; see [23].
We reduce this to a parabolic system by virtue of a change of dependent variables.Take and denote then Similarly, one calculates for  = 1, 2, 3 and for ,  = 1, 2, 3 where Γ 11 , Γ 12 , Γ 21 , and Γ 22 are the secondorder derivatives of Γ in , V. Consider the parabolic operator where "lower order terms" will be specified later.The secondorder terms vanish, provided that  is an eigenvalue and (Γ 1 , Γ 2 ) is an eigenvector.The characteristic polynomial has positive and distinct roots, because  + V +  = 1; thus and Δ 0 (, V) > 0, where The eigenvalues are equal to Their eigenvectors fulfill the algebraic equation Hence, we have the hyperbolic PDE Its characteristics have solutions Γ(, ) = const.Thus, we get two families of curves Γ (1) = const.for  1 and Γ (2) = const.for  2 .Introduce new dependent variables and the inverse functions wherever they exist.We find their gradients from the system which can be inverted as where The ultimate parabolic PDE system looks as follows: for  = 1, 2, where   =   (, ).The right-hand side of (39) consists of the "lower order terms" from (26).

Trajectories in R 3
An arbitrary smooth trajectory  = () is governed by the ordinary differential problem where  depends on , V,  and ∇, ∇V, ∇; confer [23].
Recall that  = ( 1 ,  2 ,  3 ).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].
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 Ω  with the boundary Ω  defined by (, ) = 0 where  is the solution of the PDE (45), or equivalently being constituted by the characteristics (, ()) obeying (44) starting from the initial plane  0 () = 0, the total mass of each component , V,  plus its overall inflow through the boundary is constant: Analogous identities are valid for V and .
Proof.Denote by LH = LH() the left-hand side of (46).Then, we have We differentiate LH as follows: Since LH  = 0, we have LH ≡ const.= ∫ Ω 0  0 .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.
Mathematical Problems in Engineering 4.2.Displacement Field.Since the velocity field k 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  = (, ) which gives a reference (initial) position of the marker that is currently at position ; the inverse mapping is  = (, ) providing the current position  of the marker initially placed at position .The displacement field u satisfies the relation u (, ) =  −  (, ) . (49) That is, u(, ) is the difference between the current and initial positions of the marker which is now at position  in the laboratory frame.
The integration by parts leads to the formula where  = () is the unit normal vector, proportional to ∇, and  is the surface element.A similar relation holds for V and .We can get the evolution of a region which contains some amount of , V: (a) Balance of Mass in 1D.Consider the one-dimensional case with the central position separating substances at  = 0. Thus,  0 () =  ∈ R. Suppose that the separating position  = () moves according to the mass-balance law which means that the total mass of  + V is preserved on the left of the separating position  = ().If we differentiate the identity (54) with respect to  we obtain Taking into account the governing equations, we get In this case, we arrive at the equation of the separating position trajectory or equivalently The respective plane is a solution of the PDE (b) Tracer Trajectory.Let us return to the trajectory of  which starts from a spike-shaped droplet.Suppose that we measure the middle position  3 =  3 () of this marker; hence, we have the equation Alike (54), we arrive at the differential equation 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  3 is very close to the ideal marker, since the maximum of (, ⋅) is attained at a point  such that   (, ) = 0. We draw a useful conclusion:   3 ≈ V  .We can also derive similar equations to (62) for  and V. Having in mind a severe singularity at  = 0, we can analyse the front movement of  and V as follows.Suppose that  is initially concentrated on the left  ∈ [−, 0].We study an evolution  1 =  1 () defined by with some small  > 0 (to avoid the singularity).Hence, we get the differential equation Since (, ⋅) is decreasing, we have   ≤ 0 and  1 moves to the right.In an analogous way  2 =  2 () initially placed at  = + moves to the left according to the following: (65)

Results
In the case of a binary, one-dimensional system by using the Gibbs-Duhem equation, (5), one can define the interdiffusion mobility B: and diffusion equation has the form [17] Consider the one-dimensional case with Ω  = [−, ()], where  = () satisfies (44).We formulate a simple consequence of Proposition 1 for the Kirkendall plane.
Proposition 2. If  = () satisfies (44) with (0) = 0, the total mass of each component , V,  plus its total diffusive flux through the boundary of the interval [−, ()] is constant: Analogous identities are valid for V and .
Proof.Denote by LH = LH() the left-hand side of (68).Then, we have If we differentiate this expression with respect to , we get This means that LH  = 0; hence, LH ≡ const.and the assertion follows.
Remark 3. We draw the reader's attention to the conclusion of Proposition 2: the same trajectory, given by (44) with (0) = 0, is feasible for each component , V, .This fact allows for effective computations of Kirkendall's planes.
Corollary 4 (gradients on the Kirkendall plane).Under the assumptions of Proposition 2 we have the following gradient representation on the Kirkendall plane  = (): Analogous representations have  2 V  and  3   .
Proof.If we differentiate (68), we get Because   = V  and the last term is equal to − 1   (), we obtain the assertion.

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

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  = [−1, 1].We performed computations by means of explicit FDMs with many discretization parameters.The results presented here are done for time step ℎ  = 1/500000 and space step ℎ  = 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].
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.
25 −1 −0.75 −0.5 −0.25 0 0.25 0.5 0.75 1 1.25 x (m) Consider any regular subdomain Ω 0 ⊂  in R 3 with a smooth boundary  0 () = 0.If we are interested in the evolution of  starting from Ω 0 then we obtain a family of subdomains Ω  described by the equations (, ) = 0.The mass balance leads to the integral identity