- Email: [email protected]
The lattice shift generated by two dimensional diffusion process
Computational Materials Science 95 (2014) 192–197
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
The lattice shift generated by two dimensional diffusion process Bartek Wierzba a,⇑, Marek Danielewski b a b
Rzeszow University of Technology, Faculty of Mechanical Engineering and Aeronautics, al. Powstancow Warszawy 12, 35-959 Rzeszow, Poland AGH University of Science and Technology, Faculty Materials Science and Ceramics, Al. Mickiewicza 30, 30-059 Krakow, Poland
a r t i c l e
i n f o
Article history: Received 13 February 2014 Received in revised form 16 June 2014 Accepted 9 July 2014
Keywords: Interdiffusion Bi-velocity Darken Two-dimensions Drift Potential field
a b s t r a c t The Poisson equation is used to calculate the drift velocity in the two-dimensional diffusion couple. This approach is based on the bi-velocity (Darken) method which combines the Darken and Brenner concepts that the volume velocity is essential in defining the local material velocity in multicomponent mixture at non-equilibrium. As an example the arbitrary binary system is considered. It is shown that (1) the two dimensional calculations should be applied with the stochastization method and (2) the drift term in mass conservation law does not affect the calculations. Ó 2014 Elsevier B.V. All rights reserved.
@NA : @x
1. Introduction
tdrift ¼ ðDA DB Þ
Modelling of interdiffusion phenomena is far from being unified. The specificity of different approaches lies in the arbitrary choice of the drift (convection) velocity. In computational solid mechanic it is defined as being equal to the mass average velocity [1]. Such a method allows for simplified formulation of an interdiffusion strain rate tensor and allows the Lagrangian approach for modeling interdiffusion phenomena in solid metals. In gases and fluids the drift is defined basing on the volume average velocity [2,3]. In material science the drift definition base on the Darken method [4]. This method is not widely accepted in physics and consequently it is common to postulate drift as an average mass velocity [1] or neglect drift (convection) entirely and assume simple Fick’ian diffusion only [4]. Such simplified approach to the mass transport was used in three-dimensional model of interdiffused quantum dots where an automatic solution to the Fick’s diffusion equation was implemented. Obviously it does not allow considering the interdiffusion effects [4]. The entirely new concept concerning the mass transfer started with Darken method. His basic assumption was the postulate of drift velocity in binary solid solution [5]. The drift velocity, in fact is the vacancy velocity generated during interdiffusion process caused by difference in the intrinsic diffusion coefficients. Darken determine the drift velocity from mass conservation laws of A and B component as:
Thus, the overall velocity in binary mixture was divided in two parts (1) the diffusion velocity, tdi determined by Nernst–Planck flux equation and (2) drift velocity tdrift. The drift velocity concept nowadays is still generalized. It allows to estimate: (1) the interdiffusion process in multicomponent and multiphase systems [6], (2) the Kirkendall phenomena (the position of the kirkendall planes) [7], (3) the interactions between the self stress and interdiffusion [8–10], (4) the influence of the external pressure field and entropy during diffusion process [11], (5) a mathematical model was developed to simulate the morphological change of the sintering neck by its growth and the Kirkendall void formation in the Cu–Ni alloy [12] and many others. In 2012 we have shown how to calculate the entropy and entropy production during diffusion process [13,11]. This concept was further used to calculate the interdiffusion process in multiphase binary and ternary couples. The use of the entropy production rate in multiphase calculations is an additional factor which allow to calculate the evolution of the phases in multicomponent systems was recently demonstrated [14]. The practical applicability causes that in spite of mathematical and numerical difficulties the three dimensional interdiffusion problems are still of interest. Recently Gusak et al. suggested 3D generalization of Martin–Erdelyi–Beke model of simultaneous interdiffusion and ordering in binary diffusion couples [15]. The objective of this work is to analyze the diffusion is solid solutions (binary alloy) and show the differences between the one- and two-dimensional interdiffusion calculations. The presented
⇑ Corresponding author. E-mail address: [email protected] (B. Wierzba). http://dx.doi.org/10.1016/j.commatsci.2014.07.015 0927-0256/Ó 2014 Elsevier B.V. All rights reserved.
193
B. Wierzba, M. Danielewski / Computational Materials Science 95 (2014) 192–197
approach base on the bi-velocity method where the volume velocity is essential in defining the local material velocity at nonequilibrium [8,16]. The model is formulated for an arbitrary, three-dimensional system. Discussion here is limited to the two-dimensional binary solution. 2. The two-dimensional bi-velocity method In this paragraph the two dimensional generalized Darken method will be formulated. The main law is the mass conservation law for each reacting component in the diffusion couple in the sub-set of the volume occupied by the mixture |X| X:
Z jXj
@ci þ divJ i dx ¼ 0 i ¼ 1; 2; . . . ; r: @t
ð1Þ
Above equation was derived from the Liouville theorem. The flux, Ji is defined after Darken as a sum of the diffusion ji and drift fluxes, jdrift. The ci denote the concentration of the ith element in the mixture. All of the above variables are a functions of the position in three dimensional space and time. The diffusion flux can be defined from Nernst–Planck flux equation [17,18]:
ji ¼
Di grad kT
li ;
ð2Þ
where the Di and li denote the diffusion coefficient and diffusion potential of the ith component, respectively. In this work for simplification we will apply the ideality sweeping statement (i.e. ai = ci). Thus the diffusion flux can be rewritten as follow:
ji ¼ Di grad ci ;
ð3Þ
jXj
" r X @ci i¼1
@t
# r X þ div J i dx ¼ 0
ð4Þ
i¼1
assuming that the overall concentration, c is constant and that the drift overall flux is a sum of the diffusion and drift fluxes, Ji ¼ ji þ j the following relation holds:
Z
jXj
i¼1
¼ f ðx; tÞ
ð7Þ
drift
where u denote the unknown drift potential. In order to solve the problem numerically we need to replace the second order partial derivatives with second-order finite difference approximations: drift drift udrift k1;l 2uk;l þ ukþ1;l
Dx2 k ¼ 2; . . . ; M;
þ
drift drift udrift k;l1 2uk;l þ uk;lþ1
Dy 2
¼ fk;l ;
l ¼ 2; . . . ; N
ð8Þ
where M and N denote number of nodes in x and y directions, respectively. The unknowns are located strictly in the interior of grid since udrift is known at the boundaries from the boundary conditions, and hence there are (M 1)(N 1) unknowns. The iterative update of the Jacobi iteration can now be written as: ðnÞ drift 2 ðnÞ drift 2 2 2 uk1;l þ ðnÞ udrift uk;l1 þ ðnÞ udrift kþ1;l Dy þ k;lþ1 Dx Dx Dy fk;l ðnþ1Þ drift uk;l ¼ ; 2ðDx2 þ Dy2 Þ ð9Þ Moreover, to calculate the drift velocity the Poisson equation must be solved: r X divðtdrift Þ ¼ divðgrad udrift þ rot udrift Þ ¼ divðXi ci tdi Þ i¼1
¼ f ðx; tÞ
ð10Þ
drift
After Darken the drift velocity will be calculated from the mass conservation laws of each component. The sum of Eq. (1) give the following expression:
Z
r X divðtdrift Þ ¼ divðgrad udrift þ rot udrift Þ ¼ divðXi ci tdi Þ
r X drift div ji þ ji dx ¼ i¼1
Z
jXj
r X div ji þ ci tdrift dx ¼ 0
ð5Þ
i¼1
The above set of the equations allows to calculate the interdiffusion process in two dimensional domain.
where u denote the unknown drift potential. In order to solve the problem numerically we need to replace the second order partial derivatives with second-order finite difference approximations: drift drift udrift k1;l 2uk;l þ ukþ1;l
Dx2 k ¼ 2; . . . ; M;
þ
drift drift udrift k;l1 2uk;l þ uk;lþ1
l ¼ 2; . . . ; N
Dy 2
¼ fk;l ; ð11Þ
where M and N denote number of nodes in x and y directions, respectively. The unknowns are located strictly in the interior of grid since udrift is known at the boundaries from the boundary conditions, and hence there are (M 1)(N 1) unknowns. The iterative update of the Jacobi iteration can now be written as:
3. Solution To solve above model two independent solvers should be applied. Mainly, e.g. the lines method [19] to solve the problem in space and Runge–Kutta–Fehlberg method to solve the ordinary differential equations in time. First, the uniform grid, contained M and N mesh points along the x and y direction respectively, was generated and the concentrations and drift velocity were defined at points xk,l. The equations was discretize in space to obtain the ordinary differential equation. The discretized overall flux of the ith component can be written as follow: drift
Ji ðtÞ :¼ ji ðxk;l ; yk;l ; tÞ þ j ðxk;l ; yk;l ; tÞ 2 3 kþ1;l udrift k1;l udrift kþ1;l k1;l k;l Di x ci x ci þ k;l ci x i x i ; kþ1;l k1;l kþ1;l k1;l 6 7 4 5; k;lþ1 udrift k;l1 udrift k;lþ1 c k;l1 c k;l k;l i i i i þ ci y y Di y y k;lþ1
k;l1
k;lþ1
i ¼ 1; 2; . . . ; r
k;l1
ð6Þ the udrift denote the drift potential. To estimate the drift velocity the i Poisson equation must be solved:
Fig. 1. The flow-chart depicting the various steps involved in the calculations.
194
B. Wierzba, M. Danielewski / Computational Materials Science 95 (2014) 192–197
ðnþ1Þ drift uk;l
¼
ðnÞ drift uk1;l
2 þ ðnÞ udrift kþ1;l Dy þ
ðnÞ drift uk;l1
2ðDx2
þ
2 2 2 þ ðnÞ udrift k;lþ1 Dx Dx Dy fk;l
Dy2 Þ
;
ð12Þ
To solve this model in time the Runge–Kutta–Fehlberg method was used. This method base on the Runge–Kutta 4th degree method. The main difference is the adaptative time step which allow to minimize time complexity of the algorithm [20]. Concluding the following calculations should be performed in each time step (see Fig. 1): (1) Solve Poisson equation, Eq. (7) to calculate drift potential. The drift velocity, tdrift is than a gradient of the drift potential. The Poisson equation allows determining the drift potential of the spatial distribution and hence the drift velocity can be calculated in the two-dimensional calculations. (2) Solve the conservation laws, Eq. (4). (3) Draw the results. 4. Results In this section, we present the two-dimensional evolution of the arbitrary A–B diffusion couple. We consider an ideal solid solution
(ai = ci) in which the diffusivities of components differ, D1 ¼ 1:0 109 and D2 ¼ 1:0 1010 ðcm2 s1 Þ. Other parameters used in the calculations are: temperature 1273 K; simulation time 100 h; sample length/width 0.1/0.1 cm. The initial concentration of the first (A) component was:
c0A ðt; xÞ ¼
20%; x 2 ½0; 0:05; y 2 ½0; 0:1 100%; x 2 ½0:05; 0:1; y 2 ½0; 0:1
:
In Fig. 2, the impact of the presented method is shown. The different colors present the range of the composition. It is worth to notice the bugle in the middle (along y-direction) of the presented calculations. Fig. 3 shows the computed drift potential profile in binary diffusion couple with the Poisson method presented in Section 3. It is shown, that the drift potential plot changing the sign from the minus (minimum) to the plus (maximum). The minimum and maximum are due to the difference of intrinsic diffusivities, DA > DB . Figs. 4 and 5 present the drift velocity profiles calculated from the known drift potential. It is shown that even in y-direction the drift profile is not zero and should be not neglected during the calculations.
Fig. 2. The concentration profile of the first component (A) in the binary diffusion couple after 100 h (2D simulation).
Fig. 3. The drift potential profile in the binary diffusion couple, t = 100 h.
B. Wierzba, M. Danielewski / Computational Materials Science 95 (2014) 192–197
195
Fig. 4. The drift velocity along the x-direction in the binary diffusion couple, t = 100 h.
Fig. 5. The drift velocity along the y-direction in the binary diffusion couple, t = 100 h.
Fig. 6 shows the lattice shift in the A–B diffusion couple after 100 h by 2D simulation. The lattice shift correspond to the drift velocity. Finally, in Fig. 7 we show the comparison of the one and two dimensional calculations for the binary diffusion couple after
Fig. 6. The lattice shift in the binary diffusion couple, t = 100 h.
100 h. The concentration profile is mapped from 2D calculations at y = 0.05 cm. One-dimensional profile was calculated by onedimensional Darken approximation (no Poisson equation was used). It can be noted that the boundary conditions in the 2D
Fig. 7. The comparison of the one and two dimensional calculations of the binary diffusion couple, t = 100 h. Concentration profile is mapped from 2D calculations at y = 0.05 cm.
196
B. Wierzba, M. Danielewski / Computational Materials Science 95 (2014) 192–197
Fig. 8. The stochastization method used in calculations. The left picture present the initial concentration of A component, the right one the final concentration after 100 h of diffusion.
approximation decrease the rate of the interdiffusion. More realistic analysis must allow the evolution of the couple shape (external boundary of the volume occupied by the system), e.g., X ¼ XðtÞ and jXj ¼ const. 5. Discussion The two and three-dimensional calculations may also include some concentration stochastization (stochastic distribution) as an initial condition at the Matano interface (initial contact interface). Such ‘‘stochastization’’ allows to consider the interface geometry. The convoluted interface can be a results of experimental preparation of the samples (polishing, preparation of the diffusion couple, etc.). The stochastization can be easily introduced to the calculations by ex. Gauss distribution (one can introduce the fractals). At the beginning in each node at the Matano interface (in the middle of the sample) the probability of being component A or B is calculated. The results of such treatment are presented in Fig. 8. Fig. 8 shows that the stochastisation method after large time do not influence the final results. However the shorter times can be affected by the initial interface. It is of interest that the concentration evolution during the process practically does not depend on the drift velocity but on its divergence only. The mass conservation law can be rewritten in the equivalent form:
@ci þ ci div @t
tdi þ tdi grad ci þ ci div tdrift þ tdrift grad ci ¼ 0:
In the analyzed diffusion couple the tdrift grad ci term is almost zero and can be neglected in presented calculations, only the divergence of the drift velocity influence the diffusion. Fig. 9 shows the comparison of the calculated diffusion profiles when (1) the complete mass conservation law was analyzed and (2) when the tdrift grad ci term was ignored. It is shown that both results converge. This observation might be related to some unsolved problems that depend on the drift velocity, e.g., the bifurcation of the Kirkendall plane and requires further analysis. The Kirkendall plane shift depends on the difference of the intrinsic diffusion coefficients. However, the Darken analysis of the binary diffusion couple introduces the interdiffusion coefficient that defines the diffusion flux(es) relatively to the system boundary (laboratory reference frame). In laboratory frame of reference the drift velocity ‘‘equals zero’’ and in such a case Kirkendall effect may not be estimated.
Fig. 9. The comparison of the calculated diffusion profiles when (1) the complete mass conservation law is considered and (2) when the tdrift grad ci term is ignored in the continuity equation.
The experimental procedures that allow to estimate the intrinsic diffusivities instead of interdiffusion coefficient were recently presented [21]. It may allow to redefine the drift velocity concept. Such redefinition is one of the most important scientific task to improve diffusion calculations. Acknowledgement This work is supported by the National Science Centre (NCN, Poland), decision number DEC-2011/02/A/ST8/00280. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]
E. Feulvarch, C. R. Mecanique 340 (2012) 695. H. Brenner, Int. J. Eng. Sci. 54 (2012) 67. S.D. Kokou, Phys. Lett. A 376 (2012) 3223. O. Gunawan, H.S. Djie, B.S. Ooi, Phys. Rev. B 71 (2005) 205319. L.S. Darken, Trans. AIME 174 (1948) 184. B. Wierzba, Physica A 392 (2013) 2860. O. Kirkendall, Trans. AIME 147 (1942) 104. M. Danielewski, B. Wierzba, Acta Mater. 58 (2010) 6717. G.B. Stephenson, Acta Metall. Mater 36 (1988) 2663. D.L. Beke, Defect Diffus. Forum 129–130 (1996) 9. B. Wierzba, Physica A 391 (2012) 56.
B. Wierzba, M. Danielewski / Computational Materials Science 95 (2014) 192–197 [12] [13] [14] [15] [16]
T. Yamashita, T. Uehara, R. Watanabe, Mater. Trans. 46 (2005) 88. B. Wierzba, M. Danielewski, Philos. Mag. 91 (2011) 3228. B. Wierzba, K. Tkacz-Smiech, Physica A 392 (2013) 1100. N.V. Storozhuk, K.V. Sopiga, A.M. Gusak, Philos. Mag. 93 (2013) 1999. B. Wierzba, M. Danielewski, K. Tkacz-Smiech, J. Sieniawski, Chem. Vap. Deposition 18 (2012) 267. [17] W. Nernst, Z. Phys. Chem. (N F) 4 (1889) 129.
197
[18] M. Planck, Annu. Rev. Phys. Chem. 40 (1890) 561. [19] W.E. Schiesser, The Numerical Method of Lines, Academic Press, San Diego, 1991. [20] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press, 1992. [21] I.V. Belova, N.S. Kulkarni, Y.H. Sohnc, G.E. Murch, Philos. Mag. (2013), http:// dx.doi.org/10.1080/14786435.2013.813982.
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
The lattice shift generated by two dimensional diffusion process Bartek Wierzba a,⇑, Marek Danielewski b a b
Rzeszow University of Technology, Faculty of Mechanical Engineering and Aeronautics, al. Powstancow Warszawy 12, 35-959 Rzeszow, Poland AGH University of Science and Technology, Faculty Materials Science and Ceramics, Al. Mickiewicza 30, 30-059 Krakow, Poland
a r t i c l e
i n f o
Article history: Received 13 February 2014 Received in revised form 16 June 2014 Accepted 9 July 2014
Keywords: Interdiffusion Bi-velocity Darken Two-dimensions Drift Potential field
a b s t r a c t The Poisson equation is used to calculate the drift velocity in the two-dimensional diffusion couple. This approach is based on the bi-velocity (Darken) method which combines the Darken and Brenner concepts that the volume velocity is essential in defining the local material velocity in multicomponent mixture at non-equilibrium. As an example the arbitrary binary system is considered. It is shown that (1) the two dimensional calculations should be applied with the stochastization method and (2) the drift term in mass conservation law does not affect the calculations. Ó 2014 Elsevier B.V. All rights reserved.
@NA : @x
1. Introduction
tdrift ¼ ðDA DB Þ
Modelling of interdiffusion phenomena is far from being unified. The specificity of different approaches lies in the arbitrary choice of the drift (convection) velocity. In computational solid mechanic it is defined as being equal to the mass average velocity [1]. Such a method allows for simplified formulation of an interdiffusion strain rate tensor and allows the Lagrangian approach for modeling interdiffusion phenomena in solid metals. In gases and fluids the drift is defined basing on the volume average velocity [2,3]. In material science the drift definition base on the Darken method [4]. This method is not widely accepted in physics and consequently it is common to postulate drift as an average mass velocity [1] or neglect drift (convection) entirely and assume simple Fick’ian diffusion only [4]. Such simplified approach to the mass transport was used in three-dimensional model of interdiffused quantum dots where an automatic solution to the Fick’s diffusion equation was implemented. Obviously it does not allow considering the interdiffusion effects [4]. The entirely new concept concerning the mass transfer started with Darken method. His basic assumption was the postulate of drift velocity in binary solid solution [5]. The drift velocity, in fact is the vacancy velocity generated during interdiffusion process caused by difference in the intrinsic diffusion coefficients. Darken determine the drift velocity from mass conservation laws of A and B component as:
Thus, the overall velocity in binary mixture was divided in two parts (1) the diffusion velocity, tdi determined by Nernst–Planck flux equation and (2) drift velocity tdrift. The drift velocity concept nowadays is still generalized. It allows to estimate: (1) the interdiffusion process in multicomponent and multiphase systems [6], (2) the Kirkendall phenomena (the position of the kirkendall planes) [7], (3) the interactions between the self stress and interdiffusion [8–10], (4) the influence of the external pressure field and entropy during diffusion process [11], (5) a mathematical model was developed to simulate the morphological change of the sintering neck by its growth and the Kirkendall void formation in the Cu–Ni alloy [12] and many others. In 2012 we have shown how to calculate the entropy and entropy production during diffusion process [13,11]. This concept was further used to calculate the interdiffusion process in multiphase binary and ternary couples. The use of the entropy production rate in multiphase calculations is an additional factor which allow to calculate the evolution of the phases in multicomponent systems was recently demonstrated [14]. The practical applicability causes that in spite of mathematical and numerical difficulties the three dimensional interdiffusion problems are still of interest. Recently Gusak et al. suggested 3D generalization of Martin–Erdelyi–Beke model of simultaneous interdiffusion and ordering in binary diffusion couples [15]. The objective of this work is to analyze the diffusion is solid solutions (binary alloy) and show the differences between the one- and two-dimensional interdiffusion calculations. The presented
⇑ Corresponding author. E-mail address: [email protected] (B. Wierzba). http://dx.doi.org/10.1016/j.commatsci.2014.07.015 0927-0256/Ó 2014 Elsevier B.V. All rights reserved.
193
B. Wierzba, M. Danielewski / Computational Materials Science 95 (2014) 192–197
approach base on the bi-velocity method where the volume velocity is essential in defining the local material velocity at nonequilibrium [8,16]. The model is formulated for an arbitrary, three-dimensional system. Discussion here is limited to the two-dimensional binary solution. 2. The two-dimensional bi-velocity method In this paragraph the two dimensional generalized Darken method will be formulated. The main law is the mass conservation law for each reacting component in the diffusion couple in the sub-set of the volume occupied by the mixture |X| X:
Z jXj
@ci þ divJ i dx ¼ 0 i ¼ 1; 2; . . . ; r: @t
ð1Þ
Above equation was derived from the Liouville theorem. The flux, Ji is defined after Darken as a sum of the diffusion ji and drift fluxes, jdrift. The ci denote the concentration of the ith element in the mixture. All of the above variables are a functions of the position in three dimensional space and time. The diffusion flux can be defined from Nernst–Planck flux equation [17,18]:
ji ¼
Di grad kT
li ;
ð2Þ
where the Di and li denote the diffusion coefficient and diffusion potential of the ith component, respectively. In this work for simplification we will apply the ideality sweeping statement (i.e. ai = ci). Thus the diffusion flux can be rewritten as follow:
ji ¼ Di grad ci ;
ð3Þ
jXj
" r X @ci i¼1
@t
# r X þ div J i dx ¼ 0
ð4Þ
i¼1
assuming that the overall concentration, c is constant and that the drift overall flux is a sum of the diffusion and drift fluxes, Ji ¼ ji þ j the following relation holds:
Z
jXj
i¼1
¼ f ðx; tÞ
ð7Þ
drift
where u denote the unknown drift potential. In order to solve the problem numerically we need to replace the second order partial derivatives with second-order finite difference approximations: drift drift udrift k1;l 2uk;l þ ukþ1;l
Dx2 k ¼ 2; . . . ; M;
þ
drift drift udrift k;l1 2uk;l þ uk;lþ1
Dy 2
¼ fk;l ;
l ¼ 2; . . . ; N
ð8Þ
where M and N denote number of nodes in x and y directions, respectively. The unknowns are located strictly in the interior of grid since udrift is known at the boundaries from the boundary conditions, and hence there are (M 1)(N 1) unknowns. The iterative update of the Jacobi iteration can now be written as: ðnÞ drift 2 ðnÞ drift 2 2 2 uk1;l þ ðnÞ udrift uk;l1 þ ðnÞ udrift kþ1;l Dy þ k;lþ1 Dx Dx Dy fk;l ðnþ1Þ drift uk;l ¼ ; 2ðDx2 þ Dy2 Þ ð9Þ Moreover, to calculate the drift velocity the Poisson equation must be solved: r X divðtdrift Þ ¼ divðgrad udrift þ rot udrift Þ ¼ divðXi ci tdi Þ i¼1
¼ f ðx; tÞ
ð10Þ
drift
After Darken the drift velocity will be calculated from the mass conservation laws of each component. The sum of Eq. (1) give the following expression:
Z
r X divðtdrift Þ ¼ divðgrad udrift þ rot udrift Þ ¼ divðXi ci tdi Þ
r X drift div ji þ ji dx ¼ i¼1
Z
jXj
r X div ji þ ci tdrift dx ¼ 0
ð5Þ
i¼1
The above set of the equations allows to calculate the interdiffusion process in two dimensional domain.
where u denote the unknown drift potential. In order to solve the problem numerically we need to replace the second order partial derivatives with second-order finite difference approximations: drift drift udrift k1;l 2uk;l þ ukþ1;l
Dx2 k ¼ 2; . . . ; M;
þ
drift drift udrift k;l1 2uk;l þ uk;lþ1
l ¼ 2; . . . ; N
Dy 2
¼ fk;l ; ð11Þ
where M and N denote number of nodes in x and y directions, respectively. The unknowns are located strictly in the interior of grid since udrift is known at the boundaries from the boundary conditions, and hence there are (M 1)(N 1) unknowns. The iterative update of the Jacobi iteration can now be written as:
3. Solution To solve above model two independent solvers should be applied. Mainly, e.g. the lines method [19] to solve the problem in space and Runge–Kutta–Fehlberg method to solve the ordinary differential equations in time. First, the uniform grid, contained M and N mesh points along the x and y direction respectively, was generated and the concentrations and drift velocity were defined at points xk,l. The equations was discretize in space to obtain the ordinary differential equation. The discretized overall flux of the ith component can be written as follow: drift
Ji ðtÞ :¼ ji ðxk;l ; yk;l ; tÞ þ j ðxk;l ; yk;l ; tÞ 2 3 kþ1;l udrift k1;l udrift kþ1;l k1;l k;l Di x ci x ci þ k;l ci x i x i ; kþ1;l k1;l kþ1;l k1;l 6 7 4 5; k;lþ1 udrift k;l1 udrift k;lþ1 c k;l1 c k;l k;l i i i i þ ci y y Di y y k;lþ1
k;l1
k;lþ1
i ¼ 1; 2; . . . ; r
k;l1
ð6Þ the udrift denote the drift potential. To estimate the drift velocity the i Poisson equation must be solved:
Fig. 1. The flow-chart depicting the various steps involved in the calculations.
194
B. Wierzba, M. Danielewski / Computational Materials Science 95 (2014) 192–197
ðnþ1Þ drift uk;l
¼
ðnÞ drift uk1;l
2 þ ðnÞ udrift kþ1;l Dy þ
ðnÞ drift uk;l1
2ðDx2
þ
2 2 2 þ ðnÞ udrift k;lþ1 Dx Dx Dy fk;l
Dy2 Þ
;
ð12Þ
To solve this model in time the Runge–Kutta–Fehlberg method was used. This method base on the Runge–Kutta 4th degree method. The main difference is the adaptative time step which allow to minimize time complexity of the algorithm [20]. Concluding the following calculations should be performed in each time step (see Fig. 1): (1) Solve Poisson equation, Eq. (7) to calculate drift potential. The drift velocity, tdrift is than a gradient of the drift potential. The Poisson equation allows determining the drift potential of the spatial distribution and hence the drift velocity can be calculated in the two-dimensional calculations. (2) Solve the conservation laws, Eq. (4). (3) Draw the results. 4. Results In this section, we present the two-dimensional evolution of the arbitrary A–B diffusion couple. We consider an ideal solid solution
(ai = ci) in which the diffusivities of components differ, D1 ¼ 1:0 109 and D2 ¼ 1:0 1010 ðcm2 s1 Þ. Other parameters used in the calculations are: temperature 1273 K; simulation time 100 h; sample length/width 0.1/0.1 cm. The initial concentration of the first (A) component was:
c0A ðt; xÞ ¼
20%; x 2 ½0; 0:05; y 2 ½0; 0:1 100%; x 2 ½0:05; 0:1; y 2 ½0; 0:1
:
In Fig. 2, the impact of the presented method is shown. The different colors present the range of the composition. It is worth to notice the bugle in the middle (along y-direction) of the presented calculations. Fig. 3 shows the computed drift potential profile in binary diffusion couple with the Poisson method presented in Section 3. It is shown, that the drift potential plot changing the sign from the minus (minimum) to the plus (maximum). The minimum and maximum are due to the difference of intrinsic diffusivities, DA > DB . Figs. 4 and 5 present the drift velocity profiles calculated from the known drift potential. It is shown that even in y-direction the drift profile is not zero and should be not neglected during the calculations.
Fig. 2. The concentration profile of the first component (A) in the binary diffusion couple after 100 h (2D simulation).
Fig. 3. The drift potential profile in the binary diffusion couple, t = 100 h.
B. Wierzba, M. Danielewski / Computational Materials Science 95 (2014) 192–197
195
Fig. 4. The drift velocity along the x-direction in the binary diffusion couple, t = 100 h.
Fig. 5. The drift velocity along the y-direction in the binary diffusion couple, t = 100 h.
Fig. 6 shows the lattice shift in the A–B diffusion couple after 100 h by 2D simulation. The lattice shift correspond to the drift velocity. Finally, in Fig. 7 we show the comparison of the one and two dimensional calculations for the binary diffusion couple after
Fig. 6. The lattice shift in the binary diffusion couple, t = 100 h.
100 h. The concentration profile is mapped from 2D calculations at y = 0.05 cm. One-dimensional profile was calculated by onedimensional Darken approximation (no Poisson equation was used). It can be noted that the boundary conditions in the 2D
Fig. 7. The comparison of the one and two dimensional calculations of the binary diffusion couple, t = 100 h. Concentration profile is mapped from 2D calculations at y = 0.05 cm.
196
B. Wierzba, M. Danielewski / Computational Materials Science 95 (2014) 192–197
Fig. 8. The stochastization method used in calculations. The left picture present the initial concentration of A component, the right one the final concentration after 100 h of diffusion.
approximation decrease the rate of the interdiffusion. More realistic analysis must allow the evolution of the couple shape (external boundary of the volume occupied by the system), e.g., X ¼ XðtÞ and jXj ¼ const. 5. Discussion The two and three-dimensional calculations may also include some concentration stochastization (stochastic distribution) as an initial condition at the Matano interface (initial contact interface). Such ‘‘stochastization’’ allows to consider the interface geometry. The convoluted interface can be a results of experimental preparation of the samples (polishing, preparation of the diffusion couple, etc.). The stochastization can be easily introduced to the calculations by ex. Gauss distribution (one can introduce the fractals). At the beginning in each node at the Matano interface (in the middle of the sample) the probability of being component A or B is calculated. The results of such treatment are presented in Fig. 8. Fig. 8 shows that the stochastisation method after large time do not influence the final results. However the shorter times can be affected by the initial interface. It is of interest that the concentration evolution during the process practically does not depend on the drift velocity but on its divergence only. The mass conservation law can be rewritten in the equivalent form:
@ci þ ci div @t
tdi þ tdi grad ci þ ci div tdrift þ tdrift grad ci ¼ 0:
In the analyzed diffusion couple the tdrift grad ci term is almost zero and can be neglected in presented calculations, only the divergence of the drift velocity influence the diffusion. Fig. 9 shows the comparison of the calculated diffusion profiles when (1) the complete mass conservation law was analyzed and (2) when the tdrift grad ci term was ignored. It is shown that both results converge. This observation might be related to some unsolved problems that depend on the drift velocity, e.g., the bifurcation of the Kirkendall plane and requires further analysis. The Kirkendall plane shift depends on the difference of the intrinsic diffusion coefficients. However, the Darken analysis of the binary diffusion couple introduces the interdiffusion coefficient that defines the diffusion flux(es) relatively to the system boundary (laboratory reference frame). In laboratory frame of reference the drift velocity ‘‘equals zero’’ and in such a case Kirkendall effect may not be estimated.
Fig. 9. The comparison of the calculated diffusion profiles when (1) the complete mass conservation law is considered and (2) when the tdrift grad ci term is ignored in the continuity equation.
The experimental procedures that allow to estimate the intrinsic diffusivities instead of interdiffusion coefficient were recently presented [21]. It may allow to redefine the drift velocity concept. Such redefinition is one of the most important scientific task to improve diffusion calculations. Acknowledgement This work is supported by the National Science Centre (NCN, Poland), decision number DEC-2011/02/A/ST8/00280. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]
E. Feulvarch, C. R. Mecanique 340 (2012) 695. H. Brenner, Int. J. Eng. Sci. 54 (2012) 67. S.D. Kokou, Phys. Lett. A 376 (2012) 3223. O. Gunawan, H.S. Djie, B.S. Ooi, Phys. Rev. B 71 (2005) 205319. L.S. Darken, Trans. AIME 174 (1948) 184. B. Wierzba, Physica A 392 (2013) 2860. O. Kirkendall, Trans. AIME 147 (1942) 104. M. Danielewski, B. Wierzba, Acta Mater. 58 (2010) 6717. G.B. Stephenson, Acta Metall. Mater 36 (1988) 2663. D.L. Beke, Defect Diffus. Forum 129–130 (1996) 9. B. Wierzba, Physica A 391 (2012) 56.
B. Wierzba, M. Danielewski / Computational Materials Science 95 (2014) 192–197 [12] [13] [14] [15] [16]
T. Yamashita, T. Uehara, R. Watanabe, Mater. Trans. 46 (2005) 88. B. Wierzba, M. Danielewski, Philos. Mag. 91 (2011) 3228. B. Wierzba, K. Tkacz-Smiech, Physica A 392 (2013) 1100. N.V. Storozhuk, K.V. Sopiga, A.M. Gusak, Philos. Mag. 93 (2013) 1999. B. Wierzba, M. Danielewski, K. Tkacz-Smiech, J. Sieniawski, Chem. Vap. Deposition 18 (2012) 267. [17] W. Nernst, Z. Phys. Chem. (N F) 4 (1889) 129.
197
[18] M. Planck, Annu. Rev. Phys. Chem. 40 (1890) 561. [19] W.E. Schiesser, The Numerical Method of Lines, Academic Press, San Diego, 1991. [20] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press, 1992. [21] I.V. Belova, N.S. Kulkarni, Y.H. Sohnc, G.E. Murch, Philos. Mag. (2013), http:// dx.doi.org/10.1080/14786435.2013.813982.