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

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.

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 be . 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 for all . Kirkendall’s plane(s) problem from the origin of coordinates at 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 , 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 in the three-dimensional space . The mass conservation laws have the form where and 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 via a drift velocity field 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 two-dimensional defects. The diffusion fluxes are called intrinsic (or lattice) fluxes and are related to those measured in the laboratory frame by the formula The overall molar volume and mole fractions are defined by where and is the th component partial molar volume. From (3), we get In the multicomponent solid solutions the following relations hold: 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.

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 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 ( 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 () and constant partial molar volumes such that for . 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].

2.3. 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 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 , , and , from (15) and (16) we have with . Values of belong to the interval . The initial conditions , , and are imposed for . Boundary Neumann conditions are given by on the lateral boundary , , where is the outward normal vector and is the operator of normal derivative; see [23].

3. Reduction to Parabolic Equations in

In three-component models with three distinct diffusion coefficients the third variable can be expressed in terms of the other variables , . Denote Assume that , which means that is the smallest diffusion coefficient (it is convenient to assume that ). Thus, the ternary system (17) reduces to the following: with , on the lateral boundary , and and for . We reduce this to a parabolic system by virtue of a change of dependent variables. Take and denote then Similarly, one calculates for and for where , and are the second-order derivatives of in . Consider the parabolic operator where “lower order terms” will be specified later. The second-order terms vanish, provided that is an eigenvalue and is an eigenvector. The characteristic polynomial has positive and distinct roots, because ; thus and , where The eigenvalues are equal to Their eigenvectors fulfill the algebraic equation Hence, we have the hyperbolic PDE Its characteristics have solutions . Thus, we get two families of curves . for and . for . Introduce new dependent variables and the inverse functions wherever they exist. We find their gradients from the system which can be inverted as where and for . The ultimate parabolic PDE system looks as follows: for , where . The right-hand side of (39) consists of the “lower order terms” from (26).

4. Trajectories in

An arbitrary smooth trajectory is governed by the ordinary differential problem where depends on and ; confer [23]. Recall that . 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 function is equal to the Darken (material) velocity : Suppose that the marker at is placed at a surface for . We want to trace its trajectory for which constitutes planes for . Since the characteristics are governed by (40), fulfills the hyperbolic initial-value problem because In particular, Darken’s trajectory fulfills the equation and Kirkendall’s plane satisfies the PDE where ; 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 with the boundary defined by where is the solution of the PDE (45), or equivalently being constituted by the characteristics obeying (44) starting from the initial plane , the total mass of each component plus its overall inflow through the boundary is constant: Analogous identities are valid for and .

Proof. Denote by the left-hand side of (46). Then, we have We differentiate as follows: Since , we have . 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 field is described in terms of laboratory coordinates, the displacement field 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 satisfies the relation That is, is the difference between the current and initial positions of the marker which is now at position in the laboratory frame. The displacement and velocity fields satisfy the equation

4.3. Component Trajectories

Consider any regular subdomain in with a smooth boundary . If we are interested in the evolution of starting from then we obtain a family of subdomains described by the equations . The mass balance leads to the integral identity 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 and . We can get the evolution of a region which contains some amount of :

(a) Balance of Mass in 1D. Consider the one-dimensional case with the central position separating substances at . Thus, . Suppose that the separating position moves according to the mass-balance law which means that the total mass of 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 hence, 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 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 is very close to the ideal marker, since the maximum of is attained at a point such that . We draw a useful conclusion: . We can also derive similar equations to (62) for and . Having in mind a severe singularity at , we can analyse the front movement of and as follows. Suppose that is initially concentrated on the left . We study an evolution defined by with some small (to avoid the singularity). Hence, we get the differential equation Since is decreasing, we have and moves to the right. In an analogous way initially placed at moves to the left according to the following:

Figure 1 

Distribution at for , , and . The initial distribution of the component is a very thin nanolayer around where .

5. Results

In the case of a binary, one-dimensional system by using the Gibbs-Duhem equation, (5), one can define the interdiffusion mobility : 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 , the total mass of each component plus its total diffusive flux through the boundary of the interval is constant: Analogous identities are valid for and .

Proof. Denote by the left-hand side of (68). Then, we have If we differentiate this expression with respect to , we get This means that ; hence, . 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 , is feasible for each component . 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 and .

Proof. If we differentiate (68), we get Because and the last term is equal to , we obtain the assertion.

5.1. 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 , provided that .

We return to this discussion because it is possible to simplify the reasoning. Taking into account (17) and we get Hence the uphill requirement (73) becomes Consider the case ; see Figure 4. Then, one gets and .; that is, . The uphill condition (75) is as follows: or equivalently If then any region of concavity of yields ; hence, the concentration 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 . We performed computations by means of explicit FDMs with many discretization parameters. The results presented here are done for time step and space step . 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)Figure 1 presents the cases , , , and the marker at (droplet). The initial conditions are as follows: the droplet is equal to outside the interval , and it is a polynomial on , namely, . The functions , satisfy the following properties: and ; is nonincreasing and is nondecreasing.(ii)Figures 2 and 3 show characteristics for , and the regions mapped by the transformation for , , and . In order to show the range of characteristics we take a coarse mesh with , (Figure 2), and a dense mesh , (Figure 3). Computations are made on the dense mesh, whereas the figure shows only characteristics which start from the coarse mesh.(iii)Figure 4 shows the cases , , and and a uniformly distributed substance at . One can observe uphill effects. The medium is uniformly distributed, . The functions satisfy the following properties: and ; is nonincreasing, is nondecreasing, and , are polynomials on .(iv)Figure 5 demonstrates the agreement with Belova’s results from [14]. Our computations are made with the same data as [14].

(a)

(b)

Figure 2 

Characteristics for and with and , , and .

(a)

(b)

Figure 3 

Regions mapped by the parabolic transformation. , , , and .

Figure 4 

Uphill diffusion, distributions of at and . Component is uniformly distributed for , .

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 . 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 -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: if is uniformly distributed and , then uphill diffusion of is driven by any concavity of ; 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 ? Contrary to the above expected behavior, one can find such initial distributions of that an uphill phenomenon occurs. Suppose that is an isotope and its maximum can be measured; then we have and at this maximal point. Hence, uphill condition (75) is fulfilled provided that we have at this maximal point the inequality Since at this maximum, then uphill diffusion can be generated by concavity of at this point under the mild requirement , which holds for the case .(6)There is an alternative view of the gradient in our parabolic transform; namely, one can look for a nondegenerate function whose gradient is proportional to an eigenvector corresponding to . This reads Then satisfies the equation Characteristics of this equation fulfill the systems , , and can be calculated explicitly as 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 relation leads to finding the eigenvalues in the intervals , , respectively. One can prove an analogous statement for multicomponent systems with ; however, it is not obvious whether suitable parabolic transforms can be derived, because the resulting PDE system might well be overdetermined.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported by the National Science Center (Poland) Decision no. DEC-2011/02/A/ST8/00280. Numerical experiments have been done by A. A. Majewski and M. Matusik. The authors express their thanks to J. Dąbrowa for his technical support.