Compendium About Multicomponent Interdiffusion in Two Dimensions

Interdiffusion between dissimilar solids can change the properties of joined materials. Although much work has been done to study experimentally and theoretically interdiffusion in one-dimensional (1-D) diffusion couples, studying interdiffusion in two-dimensional (2-D) or three-dimensional (3-D) solids remains a challenge. In this article, we report an experiment and develop a model to study interdiffusion in a multicomponent system of 2-D geometry. The results (concentration maps and profiles) are compared with data obtained by modeling and numerical simulations. It is assumed that the system satisfies Vegard’s rule and diffusion coefficients are composition dependent. To model the multidimensional diffusion with a drift, we take benefit of the concept of the drift potential. A nonlinear parabolic-elliptic system of strongly coupled differential equations is formulated and the implicit difference method, preserving Vegard’s rule, is applied in the simulations.

Interdiffusion between dissimilar solids can change the properties of joined materials. Although much work has been done to study experimentally and theoretically interdiffusion in one-dimensional (1-D) diffusion couples, studying interdiffusion in two-dimensional (2-D) or three-dimensional (3-D) solids remains a challenge. In this article, we report an experiment and develop a model to study interdiffusion in a multicomponent system of 2-D geometry. The results (concentration maps and profiles) are compared with data obtained by modeling and numerical simulations. It is assumed that the system satisfies Vegard's rule and diffusion coefficients are composition dependent. To model the multidimensional diffusion with a drift, we take benefit of the concept of the drift potential. A nonlinear parabolic-elliptic system of strongly coupled differential equations is formulated and the implicit difference method, preserving Vegard's rule, is applied in the simulations.

I. INTRODUCTION
INTERDIFFUSION between two dissimilar solids has long been a subject of study because of both academic interest and industrial applications. Interdiffusion plays an important role in many processes in materials and it is a basis for numerous technologies such as oxidation, brazing, or coatings. Therefore, its understanding is of considerable importance in a variety of applications and extensive efforts have been made to solve the diffusion problems. Key models of the interdiffusion base on linear phenomenological laws together with the Onsager reciprocity relations. The Onsager approach is a generalization of the Nernst-Planck flux formula to imply so-called cross effects. [1][2][3][4][5] A choice of a symmetric matrix of the phenomenological coefficients, however, can be somewhat ambiguous. Therefore, in computations, one often uses the Darken bi-velocity method, which is commonly accepted in materials science. [6][7][8] In References 9 and 10, we have shown consistency between Onsager and Darken formalisms. The Darken approach allows quantifying the Kirkendall effect and it refers to the macroscopic material movement in response to the diffusion. Kirkendall's discovery, further quantified by Darken, plays an important role in formulating the basis of the theory of interdiffusion in multicomponent material, i.e., the vacancy-dependent diffusion mechanism and convection. [11][12][13] In the Darken method, an isothermal isobaric process of diffusional mixing of the components in a closed system is described. It is assumed that the diffusion takes place in one direction, diffusion in one dimension. In a simple case, an experiment is considered in which two flat plates made of pure metals (A and B) are in contact at a common surface. They form a diffusion couple. Due to the difference of terminal compositions, the diffusion of the components (interdiffusion), driven by the concentration gradient, begins in the direction of decreasing concentrations. The species of various components can move with different average velocities. An imbalance of intrinsic fluxes, vacancy generation, and annihilation results in experimentally observed marker movement and porosity in the reaction zone of the diffusion couple. Changes in mass distribution in the diffusion couple imply such a movement of the system that a position of the center of mass would be in accordance with the principles of dynamics. In the case of the diffusion without a contribution of external forces, the position of the system can change so that the position of the overall mass center remains unchanged for an external observer (in the reference system in which the sum of the forces acting on the system is zero). This consideration provides an essence of the Darken model in which existence of a ''free'' boundary in the system such that it can move even at the absence of the external forces is postulated. A material shift, so-called drift, has been proved by numerous authors. [6,7,14,15] The diffusion fluxes arising in the system, J d i ; i ¼ 1; . . . ; s, are coupled and a certain common velocity, called drift velocity (Darken velocity (t D )), is generated. As a result, two velocities and two fluxes are assigned to each component, one due to diffusion, superscript d, and another describing Darken drift, superscript D: where s is a number of the diffusing species. If a process is governed by diffusion, the system of equations (mass continuity) has to be solved for all species: and initial and boundary conditions have to be established. A detailed list of all symbols is available in the Appendix. In many engineering applications, there are problems that involve the diffusion in solids of arbitrary three-dimensional (3-D) geometry. Numerical modeling of the transport phenomena in three dimensions, however, presents one of the most challenging applied mathematical problems. There are only a few studies dedicated to the diffusion in solids with the arbitrary shape and they are predominantly limited to Fickian diffusion. [16][17][18][19][20][21] In several studies, the geometry of the body is approximated as spheres, cylinders, or spheroids. A mathematical model to account for leaching of a solute during the heating of a parallelepiped solid was analyzed and the numerical solutions were given by Hough and Irigoyen in Reference 16. The authors considered the temperature-dependent diffusion coefficients. Carmo and Lima [17] presented a two-dimensional (2-D) mathematical model to describe moisture transport inside oblate spheroidal solids. In a series of works, da Silva and co-workers focused on modeling of the 3-D diffusion in solids of arbitrary geometry in terms of application to drying. [18][19][20] Interdiffusion and stress formation in a spherical core-shell system (lithium-ion batteries) were modeled by Verbrugge et al. [21] The authors simulated Fickian fluxes but they unacceptably neglected the interdiffusion effects. Although they calculated stress evolution in the system, they did not compare the results with the experiment. Yu et al. studied vacancy diffusion in two dimensions and plastic deformation induced by the Kirkendall effect. [12] The authors reformulated the diffusion equation to the smoothed boundary form and solved it for a deforming body. Danielewski and Wierzba modeled the 2-D interdiffusion between silicon and germanium. [14] The authors used the bi-velocity method in 1-D approximation and derived a set of partial differential equations describing the mass transport. They did not compare the results with the experiment.
It is necessary to notice here that the modeling of the interdiffusion with the drift in two or three dimensions presents a qualitatively different task compared to the 1-D case. In the 2-D and 3-D problems, the drift velocity is a vector that contains two or three components, respectively. Having one or two more unknowns, the system of equations is not closed and it is necessary to generalize the method. It has been proposed by Wierzba and Danielewski in Reference 22 to express the drift velocity by a scalar drift potential. The authors limited calculations to the interdiffusion in a 2-D arbitrary binary system. The necessity to introduce the drift potential was also mentioned by Wierzba in Reference 23, but a discussion of the model was not conclusive. Taking benefit of the concept of the drift potential, we formulated the nonlinear parabolic-elliptic system of strongly coupled differential equations and constructed the implicit difference methods to simulate the interdiffusion in a multicomponent system. [24][25][26] Exemplary numerical calculations were made for ternary systems of 2-D geometry. In Reference 25, the component diffusivities were assumed composition independent. In Reference 26, the composition-dependent diffusivities were used, and for the first time, the numerical results were compared with the experiment.
In Section III, we describe an experiment and develop a mathematical model together with a numerical scheme to simulate interdiffusion in a multicomponent system of 2-D geometry. Our article advances the field of interdiffusion of dissimilar metals, an important area having a lot of practical implications. Simulated results show acceptable agreement with experimental results for the ternary Co-Fe-Ni system annealed at 1373 K for 200 hours. We also developed a mathematical model in 3-D (Section III). Moreover, we showed that in one dimension, a well-known mathematical model is easier to construct because the drift, t D , can be calculated as a function of the first space derivatives of concentrations (Section III). Let us stress that, in the literature, there is a lack of comprehensive research of the multicomponent interdiffusion in more than one dimension.

II. DIFFUSION MULTIPLE AND EXPERIMENT
The diffusion multiple investigated in this work is an assembly of three metal blocks, Co, Fe, and Ni, subjected to annealing at a temperature of 1373 K under the vacuum of 10 À4 Pa to allow interdiffusion. The metals have been chosen not to allow the two-phase zone to grow. As it results from the phase diagram, at this temperature, the three metals show unlimited solubility.
The metals provided by Alfa Aesar (Haverhill, MA) with purity of at least 99.95 pct were used to fabricate a multiple. To achieve a thorough contact among all these pieces, the materials were ground and polished. The polishing was performed with the use of Struers MDclothes with diamond suspension gradation decreasing gradually from 9 to 1 lm. The multiple was made in two steps. First, a Fe|Ni couple was prepared by annealing for 30 minutes at 1373 K and quenching in an Al tube in a vacuum while the metal blocks were placed in the molybdenum holder to achieve intimate interfacial contact among the pieces. The prepared couple was then cut and joined in the holder with Co plate. Next, it was subjected to annealing for 30 minutes and quenching. The assembly of the three metals was cut into slices 4 mm 9 4 mm, 2-mm thick. One of the slices, arbitrarily chosen, was subjected to 199.5 hours annealing in a vacuum furnace, making the total diffusion exposure 200 hours. The diffusion time was chosen to develop diffusion profiles measurable over lengths of~200 lm. The annealed multiple was thinned and polished for scanning electron microscope (SEM) characterization and energy dispersive X-ray analysis (EDXA). We used an FEI Versa 3-D SEM equipped with an EDAX Apollo XP Silicon Drift detector. Concentration maps were measured with a resolution of 1024 9 926 points. In Figure 1, we show an SEM image of the assembled, not annealed multiple together with energy-dispersive spectroscopy (EDS) mapping; here, its fragment subjected to examination. The geometry of the studied fragment was chosen so that the three metals would have equal surfaces.
Since the mapping procedure provides a low signal-to-noise ratio, it is fraught with an error larger than line-scan analysis. To evaluate the errors, the 1-D concentration profiles from the line-scan analysis measured in the binary Fe-Ni couple, i.e., far from Co, were compared with the profile retrieved from the map. The comparison is shown in Figure 2. It is seen that the error is clearly below 10 at. pct. It appears because the signals from Co are very close to those of Fe and Ni, all three being neighbors in the periodic table. Therefore, the signals from Fe, Ni, and Co can interfere and produce errors such as an additional signal from Co present in the results from the mapping.

III. MODELING
A. Bi-Velocity in One, Two, and Three Dimensions In References 24 through 26, we have presented a basic formulation of the interdiffusion models in one, two, and three dimensions that is a generalization of the Darken bi-velocity method. Let an open and bounded set X Ì R n , n = 1, 2, 3 with the smooth boundary ¶X be given. Moreover, let t p means the processing time. In two and three dimensions, a key assumption of the model is that there exists a scalar potential function, F = F(t, x), such that the vector drift velocity is expressed by the minus of its gradient: The existence of such potential is appropriate, even above the Tammann temperature (two-thirds melting temperature), because the viscosity of solids is very high and rot t D ¼ 0. In Reference 25, the model has been limited to the interdiffusion of the species having concentration-independent diffusivities. In Reference 26 the next step has been made and the system of the species having concentration-dependent diffusivities has been studied. The numerical results have been compared with the experiment. Here, we present a compendium about 1-D and multidimensional multicomponent interdiffusion based on our previous articles. The system under consideration has s ‡ 2 components and it is represented by molar concentrations, c 1 ; . . . ; c s ; c i ¼ c i t; x ð Þ and molar ratios We assume incompressible transport. In such case, the Euler theorem allows expressing the molar volume of a mixture of the components, X m , having partial molar volumes, X i = const, i = 1, ..., s, as: which is accepted as Vegard's rule and equivalently reads: The flux of the ith species has two components, diffusion and drift, and we assume: where the diffusivity is composition dependent, D i = D i (c 1 , ..., c s ). Note that Eq. [7] is a generalization of the Fick flux formula. For simplicity, we will apply in the following the diffusivities represented in terms of mole fractions, D i = D i (N i ). The flux J i enters the continuity equation (Eq. [3]). Multiplying Eq. [3] by X i , adding by its sides, and using the Vegard rule (Eq. [6]), we obtain the volume continuity equation: Define the initial condition on the concentrations: where n is the outside normal to the boundary ¶ X.
Assume the Vegard rule on the initial concentrations: We first formulate the 1-D model. [27][28][29] Let X = ( À K, K) Ì R. It follows from Eq. [8] that: where K is an arbitrary function. By Eqs. [10] and [12], we get the unique: KðtÞ ¼ X s i¼1 X i j i ðt; KÞ ¼ À X s i¼1 X i j i ðt; ÀKÞ: ½13 The second equality in Eq. [13] can be treated also as an assumption on the boundary evolutions j i . On the other hand, Eqs. [6] and [7] imply: In consequence, from Eqs. [9], [13], and [14], we obtain the formula on the drift: Finally, we get the strongly coupled nonlinear parabolic differential system: and the nonlinear coupled boundary conditions: Now we formulate 2-D and 3-D models. [24][25][26] Let X Ì R n , n = 2, 3. In these cases, the volume continuity equation (Eq. [8]) cannot be applied for finding the drift because it does not imply independency with respect to x of the field under divergence. Assuming that the drift velocity can be written as minus of the gradient of the scalar field F(t, x) (Eq. [4]), we get the strongly coupled nonlinear parabolic-elliptic differential system: .., s. The first of the preceding equations is the continuity equation (mass conservation) (Eq. [3]). The second one is immediately generated by Eq. [8] and it is an elliptic equation for potential F. The third equation represents the uniqueness property for the second Poisson equation for F together with the boundary condition defined later. For the system Eq. [18], we want to find the boundary conditions. Multiplying Eq. [7] by X i and n, adding by its sides, and using the Vegard rule (Eq. [6]) and Eq. [10], we have: Hence, the nonlinear coupled boundary conditions are implied: À D i @c i @n À c i @U @n ¼ j i t; x ð Þ on ½0; t p Â @X; @U @n ¼ À X s k¼1 X k D k @c k @n þ j k ðt; xÞ on ½0; t p Â @X;

> > < > > :
½20 for i = 1, ..., s. If an incompressible transport is assumed (constant overall volume), the compatibility condition on the boundary evolutions follows from the Gauss theorem: Let us stress that in 2-D and 3-D cases, the drift t D has two or three coordinates, respectively. Hence, we cannot describe the diffusion model by the differential algebraic system (Eqs. [3] and [6]), because a number of unknowns is larger than a number of equations. In both cases, we have only (s + 1) equations: s continuity equations and one Poisson equation. By elimination of the vectorial drift-velocity component (Eq. [4]), we have the same number of unknowns: the concentrations of s components and the scalar potential. Obviously, in one dimension, this number is also the same. However, from a numerical and analytical point of view, it is better to study the equivalent differential system (Eq. [16] in one dimension and Eq. [18] in two and three dimensions). It is interesting that the parabolic-elliptic model Eqs. [18], [9], and [20] can also be used in one dimension. Moreover, it leads, after integration of the elliptic equation for F, to the parabolic model (Eqs. [16], [9], and [17]). [25] B. Difference Scheme in Two Dimensions To find numerical solutions for the model, we use the implicit finite difference method (FDM). Let us focus on the 2-D case, X = (À K, K) 9 (À K, K) Ì R 2 , x = (x 1 , x 2 ). It is assumed that the system is closed, e.g.: The initial concentrations are given (Eq. [9]). The differential system (Eq. [18]) and the boundary conditions (Eq. [20]) can be written in the following form, respectively: and for j = 1, 2, i = 1, ..., s-1. The concentration of the sth species, c s , is calculated as: The modification with the use of Eq. [25] allowed construction of the Vegard rule preserving the implicit finite difference method, generated by some linearization and splitting ideas used, to present the numerical solutions. Numerical methods that preserve some futures of a continuous model are very important and interesting from a physical and mathematical point of view. Because of the nonlinearity of the diffusion coefficients, we applied a middle points discretization that admits a special subtle approximation. In fact, two implicit linear difference schemes were defined for the system (Eqs. [23] and [24]), with the initial condition (Eq. [9]): 1. for the elliptic subsystem on the potential F; 2. for the parabolic subsystem on the concentrations c i , i = 1, ..., s-1.
The Gauss elimination method has been applied to solve them. For details concerning the existence and uniqueness of the discrete solutions and other properties of the numerical method, we refer the readers to our previous articles. [25,26] C. Data Given the model and numerical method, the computations were made to simulate the interdiffusion provided in the experiment. The necessary data are given in Table I. The initial concentrations given in the table include the noises interpreted by the EDS software as additional signals from the remaining elements. The partial molar volumes of the components were estimated from the density of the pure metals. [30] The concentration-dependent diffusion coefficients of the species in the ternary Co-Fe-Ni mixture were calculated using the modified Boltzmann-Matano method proposed by Wierzba and Skibinski [31] with further modification by Zajusz et al. [32] (Table I, Figure 3).

IV. RESULTS AND DISCUSSION
We begin presentation of the results with a comparison of experimental and computed concentration maps for the Co-Fe-Ni multiple subjected to annealing (interdiffusion) at 1373 K for 200 hours (Figures 4(a) and  (b)). To reduce noise, the experimental data were smoothed with the Savitzky-Golay filter. [33] Good agreement of both sets of data is seen. Small exceptions can be due to the Frenkel porosity in iron that disrupts the experimental results. It can also be seen that the bottom of Figure 4(b), i.e., far from Co, is actually a binary diffusion couple of Fe|Ni. Similarly, the left and right sides, far from Ni and Fe, respectively, are from the Fe|Co and Ni|Co couples. The positions of the Matano planes drawn in the graphs were determined from the concentration profiles at the binary cross sections, Fe|Ni, Co|Ni, and Co|Fe, using the procedure by Leszczyn´ski et al. [34] Below the maps, we show vector flux fields J i , i = Co, Fe, Ni visualized as a collection of arrows with a given magnitude and direction. In Figure 5, the binary cross sections, SEM images, and concentration profiles at the end of the processing are shown. We show experimental concentration profiles and the results of simulations made using the 1-D and 2-D models. Note that the deviations from zero are due to the error of the EDS-2-D mapping procedure ( Figure 2). For all the couples, the 1-D and 2-D models give overlapping profiles, which confirms the correctness of our 2-D model and 2-D simulations. Moreover, for the Co|Ni couple, the experimental and simulated concentration profiles agree perfectly. A good agreement is also seen for the Co|Fe couple. In the case of the Fe|Ni couple, the simulated data do not thoroughly overlap with the experimental data. Close to the Matano plane (red dashed line), the data differ, which is probably due to the void formation in iron.
The detailed SEM photography for the binary cross sections is shown in Figure 6 together with the elemental mapping. It confirms lack of cracks, additional phase precipitations, and other than void elements, which can have a significant influence on diffusion and further interpretation. It also confirms the claim that interdiffusion is planar.
The void formation (the Frenkel effect) and the resulting porosity are common for the interdiffusion at intermediate temperature and are usually related to the difference of diffusivities of the species, which is the case here. To study the effect, we calculated the drift velocity and its divergence ( Figure 5(b)). According to expectations, the positive drift velocity decreases and its extremes move away from the Matano plane over time. The porosity in Fe in Co|Fe and Fe|Ni couples overlaps with the positive drift velocity divergence at the 25 hours, divt D , range. Not all results, however, are unambiguous, such as the lack of porosity in Ni, in the Fe|Ni couple, despite positive divt D , or porosity on both sides of the Matano plane in the Co|Ni couple, although divt D is negative in Ni only. Therefore, it is difficult to draw reliable conclusions concerning the formation of voids. Many experiments confirm that the porosity varies from one specimen to another and the most reasonable explanation is that there may be a difference in the size of the nuclei for pores in various specimens. [35] The results of our experiment confirm this reasoning.  [30][31][32] Mixture Co-Fe-Ni X ( À K, K) 9 ( À K, K), K = 0.042 cm In the next set of results (Figures 7 and 8), we show the concentration profiles at the four other cross sections, as indicated in Figure 9. Together with the concentration profiles, the diffusion paths are shown for various processing times. Every diffusion path connects the terminal (initial) compositions (here, pure metals) and traces the changes in average composition of the diffusion zone.
In Figures 7(a) and (b), the results along two diagonals of the cross sections of the studied rectangle are presented. The results for the time t = t p = 200 hours are a comparison of numerical results from the 1-D and 2-D models with experimental data. Clearly, better agreement in comparison to the experiment is seen in the case of the 2-D model, which shows its advantage over the well-known 1-D model. The time evolution presented as a series of the diffusion paths shows progressive mixing of the species. At the beginning, t = 0, a system consists of two diffusion couples, respectively, Co|Fe, Fe|Ni in Figure 7(a) and Fe|Ni, Ni|Co in Figure 7(b). With time progression, the diffusion zone grows and the species mix in the center area of the couple while the terminal compositions remain the same.
The advantage of the 2-D model can be seen even better when the cross sections pass through a point of contact of the three metals. The 1-D model in this case cannot be applied as it ignores the third component. Focusing on the benefits of modeling in two dimensions, the results confirm entrance of the third species into the couple, i.e., Fe into Co|Ni and Ni into Fe|Co in Figures 8(a) and (b), respectively. The gradually growing amount of the third species in the initial two-component system affects the diffusion paths. A worse agreement in the case of the subcouple Fe|Co can be due to the voids.
The presented models of interdiffusion are limited to the high-temperature range, above the Tammann temperature. At intermediate temperature voids, formations, stresses, and grain boundaries are decisive and our models are not adequate.

V. SUMMARY
Many processes that occur in materials are diffusion controlled; therefore, diffusion is a subject of numerous theoretical (modeling) and experimental studies. Considerable progress in the understanding of diffusion in materials started after Kirkendall's discovery of an effect referred to by his name and later given a mathematical description by Darken. Modeling of the interdiffusion in multicomponent systems presents a difficult task and is usually simplified to the 1-D diffusion in diffusion couples.
In the present work, we go beyond this scheme and present a method to study experimentally and theoretically (model and numerical simulations) the interdiffusion in multidimensional systems. The model is based on the assumption that the drift velocity can be replaced by minus of the gradient of the scalar potential. In this way, a number of unknowns is reduced to a number of equations, i.e., (s+1), and we get the closed, mathematically well-posed parabolic-elliptic system of differential equations with initial-boundary conditions. If in the 2-D and 3-D cases X is a connected region and rott D = 0, then there exists the potential F of the drift t D . Nearly all of the diffusion processes in solids show negligible turbulence; thus, it can be neglected and the postulate about rott D = 0 is appropriate. A series of numerical simulations for the ternary Co-Fe-Ni multiple and comparison with the experimental results confirm the reliability of the method, as well as the validity and the admissibility of the assumption about the drift potential. The numerical simulations provide time evolution, which can be useful, e.g., in designing new materials. A very important achievement in our research is that we solve the nontrivial mathematical problem. While in one dimension the interdiffusion is given by the parabolic system of the differential equations, in more dimensions, the system is parabolic-elliptic. Both theory and numerical methods to solve such systems are difficult and rather poorly researched. Therefore, making numerical computations for such a system is nontrivial.

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source,  provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat ivecommons.org/licenses/by/4.0/.