Combined probability distributions of random-walks: A new method to simulate diffusion processes

Computational Materials Science 34 (2005) 254–263 www.elsevier.com/locate/commatsci

Combined probability distributions of random-walks: A new method to simulate diffusion processes ˚ gren Henrik Larsson *, John A Division of Physical Metallurgy, Department of Materials Science and Engineering, Royal Institute of Technology (KTH), Brinellv. 23, SE-100 44 Stockholm, Sweden Received 25 August 2004; received in revised form 24 January 2005; accepted 7 February 2005

Abstract Two related methods to simulate diffusion processes are presented. Both are based on conceiving diffusion as a random-walk process. Several example simulations are presented: single phase diffusion couples, diffusion controlled growth and prediction of Kirkendall porosity. Comparisons with experimental results and simulation software based on established technique (DICTRA) show good agreement with results obtained from the presented models. Ó 2005 Elsevier B.V. All rights reserved. PACS: 02.60.Cb; 05.60.-k; 66.10.Cb; 81 Keywords: Diffusion; Simulation; Multi-component; Interface movement; Phase transformation; Kirkendall effect

1. Introduction The purpose of the present work is to show that the random-walk technique can be used to solve many different types of diffusion problems. In addition, the technique opens up a somewhat different perspective on diffusion and diffusion controlled processes that may yield additional insight. Most results presented in this paper have *

Corresponding author. Tel.: +46 8 790 8308; fax: +46 8 207 681. E-mail address: [email protected] (H. Larsson).

been reported earlier, see Refs. [1–3], but in a rather condensed and divided manner. Thus it is hoped that this report will show that the technique is both simple and quite versatile, although many aspects remain to be investigated. In the present work only one-dimensional, isothermal cases will be considered. Furthermore, with regards to phase transformations, we will only discuss systems without interstitial species and with the molar volume Vm approximated to be constant. In order to clarify the simulation models discussed in the present paper, we start with a brief

0927-0256/$ - see front matter Ó 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2005.02.004


gren / Computational Materials Science 34 (2005) 254–263 H. Larsson, J. A

overview of diffusion in general and random-walk in particular.

The mean square displacement of an atom in an ideal solution is given by EinsteinÕs relation Z 2 ¼ 2D t

2. Diffusion and random-walk The variation of energy with position of an atom has been schematically drawn in Fig. 1. Let us first consider an ideal solution, cf. the top left graph in Fig. 1. The corresponding flux of k-atoms is given simply by the amount by which the diffusing species differ over a certain length times the so-called tracer diffusivity, viz. J ideal ¼
k

1  oxk D V m k oz

ð1Þ

The mobility Mk is approximately related to the tracer diffusivity Dk by Dk ¼ M k RT

ð2Þ

Thus J ideal k

1 oðRT ln xk Þ ¼
M k xk Vm oz

ð3Þ

255

ð4Þ

In a non-ideal solution there can, beside a gradient in ideal configurational entropy o(RT ln xk)/ oz, exist gradients in (lk
RT ln xk), cf. the top right graph in Fig. 1, and in mobility, Mk, cf. the bottom left graph in Fig. 1. Regarding the gradient in (lk
RT ln xk) it is easily shown, see e.g. Bardeen [4], to yield a total flux Jk ¼


1 ol M k xk k Vm oz

ð5Þ

if the gradient is small and any vacancy chemical potential gradient is disregarded. Eq. (5) is commonly used in the so-called lattice-fixed frame of reference. A gradient in mobility does not affect the flux; it is taken into account when the equation of continuity 1 oxk o ¼ ð
J k Þ V m ot oz

ð6Þ

is invoked.

Fig. 1. Schematic view of how the energy of an atom may vary with its position in the lattice. x-axis is distance and y-axis is energy. E* is an activation energy for diffusion. Top left graph shows an ideal solution. Top right indicates the presence of a gradient o(lk
RT ln xk)/oz; the energy difference between neighbouring planes is DE. Bottom left indicates the presence of a gradient in mobility oMk/oz; the difference in activation energy in different directions from a certain plane is DE*. In the bottom right graph both types of gradients are present. After Le Claire [5].


gren / Computational Materials Science 34 (2005) 254–263 H. Larsson, J. A

256

The gradients depicted in Fig. 1 give rise to an atomic mean drift velocity, see e.g. Refs. [5,6]. Regarding the gradient in (lk
RT ln xk) it yields a contribution to the drift proportional to ðJ k
J ideal Þ, as the mean displacement in ideal difk fusion is zero. The drift caused by the mobility gradient is easily shown to be RT oMk/oz, if the fluxes from a given plane are considered. Thus the total mean drift velocity of atom k obeys: vk ¼

Vm oM k ðJ k
J ideal Þ þ RT k xk oz

ð7Þ

The overall movement of an atom may be emulated by an ideal random-walk superimposed on which is a drift velocity. This behaviour corresponds to an Itoˆ stochastic differential equation (SDE) pffiffiffiffiffi dZ k ¼ Ak dt þ Bk dW ð8Þ where Zk is the position of a certain atom of species k along the z-axis, Ak and Bk are the drift and diffusion parameters, respectively, and W is a stochastic variable. The Itoˆ SDE yields a probability distribution pk which can be obtained by solving the corresponding Fokker–Planck equation [7] opk o 1 o2 ¼
ðAk pk Þ þ ðBk pk Þ oz 2 oz2 ot

1 o2 ð2M k RTxk Þ 2 oz2 ð10Þ

By comparing Eqs. (9) and (10) the parameters Ak and Bk may be identified as Ak  vk ¼ RT

oM k o
M k ðlk
RT ln xk Þ oz oz

Bk ¼ 2M k RT ¼ 2Dk

where n is a normal distributed, stochastic variable with mean zero and variance one. Z is the position of a walker. The subscript k was dropped for clarity. It must be noted that the technique outlined above is not atomistic. The walkers represent some arbitrary amount of mass of a certain species and they are not confined to any specific locations such as lattice sites.

ð9Þ

However, in the limit of a large system this probability distribution is identical to the mole fraction xk. By combining Eqs. (5) and (6) and rewriting in Fokker–Planck form we obtain    oxk o oM k o RT ¼

M k ðlk
RT ln xk Þ xk oz oz ot oz þ

the combined Eqs. (5) and (6) in Fokker–Planck form we obtain the drift velocity vk of Eq. (7) as Ak of Eq. (11) and the diffusion parameter Bk corresponding to Eq. (4). Simulations may be performed by introducing point carriers of mass, which we refer to as walkers, and solving the Itoˆ SDE individually for a sufficiently large number of walkers. It is clear that, in the limit of infinitely many walkers and as the time-step tends to zero, an exact solution of the diffusion equation is obtained. However, in order to evaluate the parameters it is also necessary to estimate the local composition, which the parameters in general depend upon. During simulations, the Itoˆ SDE is solved by the Euler forward scheme, viz. pffiffiffiffiffiffiffiffi Z iþ1 ¼ Z i þ ADt þ BDtni ð13Þ

ð11Þ ð12Þ

It is clear that the parameter Ak represents a deviation from ideality. Note that by simply rewriting

3. Example simulation using walkers: The darken experiment In 1948 Darken performed this famous experiment [8]. Two pieces of steel with similar carbon content but dissimilar silicon contents were joined together in a diffusion couple. After heat-treatment for 13 days at 1050 °C it was found that carbon, diffusing much quicker than silicon, had redistributed extensively, exhibiting so-called uphill diffusion. A corresponding simulation was performed using only walkers representing carbon; no significant redistribution of silicon occurs. Local mean compositions were estimated during each timestep, by dividing up the system into a number of ‘‘cells’’ and simply counting the number of walkers present in each; the silicon content in each cell was constant during the simulation, 3.8 mass% in the


gren / Computational Materials Science 34 (2005) 254–263 H. Larsson, J. A

left half and 0.05 in the right. The initial carbon content was 0.49 mass% in the left half and 0.45 in the right. Gradients were evaluated by finite differences. The code was written in Matlab [9]. Independently assessed thermodynamic and kinetic parameters were obtained from the Thermo-calc program [10] by means of a specifically developed interface connecting the two softwares [11]. The carbon concentration profile resulting from the simulation is presented in Fig. 2 together with the corresponding experimental data. As can be seen, the agreement is satisfactory. The ‘‘wiggles’’ in the curve is due to the stochastic nature of the technique.

4. Random-walk and diffusion controlled growth In this section we will consider a reaction b ! a controlled by diffusion in the bulk, i.e. the situation in the interface will be neglected. In ‘‘sharp interface methods’’ local equilibrium is often assumed to hold at the phase interface. The interfacial velocity is then obtained by solving the flux balance equation, i.e. the Stefan condition

 v  a xk
xbk ¼ J ak
J bk Vm

257

ð14Þ

which must hold for all species in order to conserve mass. v is the interfacial velocity and xak and xbk are the mole fractions of k adjacent to the interface on the a and b side, respectively. Eq. (14) must obviously yield the same velocity for all species. In a binary system the composition at the interface is given directly by the phase diagram, if local equilibrium is assumed, but in ternary and higher systems there are an infinite number of possible tie-lines. In those cases the operating tie-line must be found by some interative procedure. In a system containing only substitutional elements of equal partial molar volume, the interfacial velocity is given directly by, see e.g. Ref. [12], v ¼
V m J m

ð15Þ

m

where J is the net flux across the phase interface. Eq. (15) is usually not directly applicable because it is impossible to directly evaluate J m using traditional approaches. However, by using point carriers of mass—walkers—as described above, the net flux J m is given simply by counting the net passage of walkers during a time-step.

Fig. 2. Result from a simulation of a diffusion couple consisting of two pieces of steel with similar carbon contents but dissimilar silicon contents. y-axis shows carbon content in mass%. x-axis shows distance [m]. Experimental results (circles) from Darken [8].


gren / Computational Materials Science 34 (2005) 254–263 H. Larsson, J. A

258

5. Example simulation using walkers: diffusion controlled growth Diffusion controlled growth was simulated for a hypothetical ternary A–B–C regular solution system. The chemical potential of component A is given by 

lA ¼ lA þ RT ln xA þ xB ðxB þ xC ÞLAB þ xC ðxB þ xC ÞLAC
xB xC LBC þ ð1
2xA ÞLABC ð16Þ and similarly for the other components. For simplicity the mobility M was assumed to be constant and equal for all species, in both phases. The system size was set to 100 lm and it was assumed that a layer of a 1 lm thick had nucleated on one side of the system of which the rest was b. One kind of walker was defined for each species. During each time-step the system was divided into a number of ‘‘a’’- and ‘‘b’’-cells; the Gibbs energy of each phase are given by separate functions. Note that no thermodynamic condition is imposed on the phase interface; the situation at the interface comes out as a result of the simulation and it tends towards local equilibrium. As mentioned above, the displacement of the phase interface during a time-step was given directly by the net passage of walkers across it. The phase boundary position as a function of time is presented in Fig. 3. For comparison, the result from a corresponding simulation using the DICTRA software [13] has been included.

6. Combined probability distributions of random walks—‘‘the distribution function method’’ The method based on walkers—point carriers of mass—discussed above has many attractive features, it is for example very easy to program. However, for some types of simulations it is of interest to obtain noise-free solutions. In order to retain the basic features of the ‘‘discrete model’’, but still make it possible to obtain smooth and continuous curves, the following method is suggested. As indicated by the title of this section, it is based on considering the system as made up of distributions of

the participating species and evaluating the temporal evolution of these. It is well known that the probability distribution of an ideal random-walk is Gaussian. Specifically, if a species k is initially concentrated at z = 0 in an amount of s moles per unit area, the concentration profile is given by  2 1 s
z exp xk ¼ pffiffiffiffiffiffiffiffiffiffiffiffi ð17Þ Vm 4Dk t p4Dk t Comparing the Gaussian, Eq. (17), with the Itoˆ SDE, see Eqs. (8), (11) and (12), it is seen that the parameter Bk equals, as before, 2Dk . In the non-ideal case, the drift yields the expectancy, i.e. the drift velocity Ak multiplied with t yields the displacement of the distribution and should be included in the exponential. Thus, instead of the discrete walkers used in the example simulations presented above, continuous probability distributions may be used. In the ideal case, the probability distributions are Gaussians with a variance equal to 2Bkt. Deviations from ideality gives rise to a non-zero expectancy equal to Akt. For an ideal case when a species k is initially evenly distributed between z = 0 and infinity a solution may be obtained by integrating over Gaussians, a so-called error function solution ! Z 1 2 xk0
ðz
~zÞ pffiffiffiffiffiffiffiffiffiffiffiffi xk ¼ d~z ð18Þ exp 4Dk t p4Dk t 0 If, during each time-step, Bk and the product Akt are considered as constants in each in a series of Gaussians, then, as indicated above, a solution may be written as Z zmax xk0 ðzÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi xk ðzÞ ¼ p2Bk ðzÞt zmin ! 2
ðz
Ak ðzÞt
~zÞ  exp d~z ð19Þ 2Bk ðzÞt in analogy with Eq. (18). The species is initially distributed betweeen zmin and zmax as given by the distribution function xk0(z). In Fig. 4 the division of species into these functions have been indicated; in the interior of phases overlapping triangular distribution functions are used. At


gren / Computational Materials Science 34 (2005) 254–263 H. Larsson, J. A

259

Fig. 3. Results from simulations of a b ! a transformation in a hypothetical ternary system using a model described in this work (—) and a corresponding simulation using the DICTRA software [13] (- - -). x-axis shows time [s]. y-axis shows the position of the phase interface [m].

Fig. 4. Schematic view of how species are divided up into distribution functions, the roughly triangular areas indicated in the figure.

phase interfaces sharp transitions between distribution functions are used (not shown in figure). The chemical potentials and mobilities at the center of each distribution function is estimated from local mean compositions xk;i as indicated in Fig. 5. Estimates of gradients are made by finite differences using the discretisation zi. The temporal

evolution of a distribution function during a timestep is schematically drawn in Fig. 6. When evaluating the integral Eq. (19) the parameters Ak and Bk are assumed to vary in a piecewise linear fashion between the centers of the distribution functions. Mass-balance is preserved by checking, and, if necessary, making a small adjustment to the mass content of each distribution function. The size of the time-step is chosen such that the ‘‘mean interaction distance’’ is sufficiently smaller than the distance between the centers of neighbouring distribution functions. For moving phase boundary problems the principle for obtaining the interface displacement during a time-step is similar to the discrete approach. With reference to Fig. 7 the net mass flow J m across the interface is given by the difference in size of the shaded areas. A simulation time-step proceed as follows. First the integral Eq. (19) is solved numerically for all species and distribution functions, yielding a ‘‘preliminary’’ (this is discussed below) concentration profile. The concentration profile for species k at z = zm is thus given by summing over all

260


gren / Computational Materials Science 34 (2005) 254–263 H. Larsson, J. A

Fig. 5. Schematic view of evaluation of local mean compositions xk;i and discretisation zi for calculation of parameters Ak and Bk.

Fig. 6. Schematic view of the temporal evolution of a distribution function during a time-step.

distribution functions i, the contribution from each distribution function being given approximately by the sum over j, viz.

xprel k ðz

The profile obtained by summing together all mole-fractions may differ from unity due to the Kirkendall effect and non-zero interface velocities, thus the notion of the ‘‘preliminary’’ profile. The profile is therefore normalised by expanding it where the sum is larger than one and contracting it where it is less. If the diffusion problem being solved involves phase interfaces, the interface displacement is obtained ‘‘automatically’’ in the process. This normalisation procedure can also be said to correspond to a transformation from the lattice fixed frame of reference to the volume fixed. An estimate of the Kirkendall shift can be obtained by tracking the movement of one or more reference points as the grid is expanded and contracted. In the bulk of a phase, contraction of the grid corresponds to annihilation of vacancies. Thus a measure of the local probability for formation of Kirkendall porosity is obtained by keeping a record of the fraction of vacancies annihilated as a function of position. The sum over all time-steps

" ! 2 XX1 xk0;ij
ðzm
Ak;j Dt
zj Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp ¼ zm ; t þ DtÞ ’ 2 2Bk;j Dt p2Bk;j Dt i j 2 xk0;iðjþ1Þ
ðzm
Ak;jþ1 Dt
zjþ1 Þ þpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp 2Bk;jþ1 Dt p2Bk;jþ1 Dt

!#  ðzjþ1
zj Þ

ð20Þ


gren / Computational Materials Science 34 (2005) 254–263 H. Larsson, J. A

261

much smaller spatial step-size than the discretisation (indexed with i in Fig. 5) used for evaluating the parameters at the center of distribution functions.

7. Example simulation using distribution functions: Prediction of Kirkendall porosity

Fig. 7. Schematic view of how the net mass flow across an interface can be evaluated. The position of the interface is given by the vertical line. The net mass flow is given by the difference in size of the shaded areas.

yields an upper limit on the fraction of pores formed due to the Kirkendall effect. Finally, it must be noted that the discretisation used when evaluating the integral (indexed with j in the trapezoid approximation of Eq. (20)) has a

Nesbitt and Heckel [14] studied interdiffusion in the Ni–Cr–Al system. A large number of diffusion couples were investigated. In the majority of the specimens Kirkendall porosity formed. The temporal evolution of two Ni–Cr–Al diffusion couples were simulated using the distribution function method as described above. Independently assessed thermodynamic and mobility data were used as described in Section 3. The simulations correspond to anneals at 1200 °C for 100 h. A concentration profile obtained from one of the simulations along with experimental data due to Nesbitt and Heckel is shown in Fig. 8. The obtained integrated (summed over all time-steps) vacancy profiles are shown in Fig. 9. Approximate measures, as inferred from micrographs in Nesbitt

Fig. 8. Result from a simulation of a ternary Ni–Cr–Al diffusion couple showing concentration profiles of chromium and aluminum (solid lines). x-axis shows distance [m]. y-axis shows mole fractions. Experimental data (symbols) from Nesbitt and Heckel [14].

262


gren / Computational Materials Science 34 (2005) 254–263 H. Larsson, J. A

Fig. 9. Results from simulations of ternary diffusion couples showing integrated vacancy fraction (y-axis) as function of distance (x-axis [m]), which correspond to upper bounds of the Kirkendall porosity. Thick lines show approximate extent of experimentally observed Kirkendall porosity as inferred from micrographs in Nesbitt and HeckelÕs article [14].

Fig. 10. Results from simulations of a b ! a transformation using a model described in this work (—) and a corresponding simulation using the DICTRA software (- - -). x-axis shows distance [m]. y-axis shows mole fraction of specie A.

and HeckelÕs work, of the width of the regions in which Kirkendall porosity formed are indicated.

The simulation profiles yield upper bounds for the amount of Kirkendall porosity. If it is assumed


gren / Computational Materials Science 34 (2005) 254–263 H. Larsson, J. A

that a vacancy fraction of about 0.007–0.008 may annihilate without causing porosity there is good agreement between simulations and experiments.

8. Example simulation using distribution functions: Diffusion controlled growth Diffusion controlled growth of a into b was simulated for a hypothetical binary system using the distribution function method as described above. Both phases were modelled as regular solutions, i.e. the chemical potentials were given by expressions similar to Eq. (16). The mobility was assumed constant and equal for both species in both phases. The concentration profile after a certain time is presented in Fig. 10. For comparison, result from a corresponding simulation using the DICTRA software [13] is included. As can be seen, the agreement is satisfactory although it is not perfect.

263

˚ grenÕs work on Ostwald ripening Schwind and A [15]. Both methods show interesting properties. It is interesting to note that free boundary problems can be solved without imposing any type of thermodynamic condition on the interface. Especially the discrete method is intuitive in this respect. Perhaps the most interesting feature with the distribution function method is that the diffusion problem is solved directly in the lattice fixed frame of reference. For simulations of phase transformations both methods essentially require that local equilibrium is a good description of the state at the interface, i.e. that the reaction is diffusion controlled. As diffusion of species is modelled as always occurring down its chemical potential gradient it will not be possible to model e.g. solute trapping. For those cases some additional mechanism for mass transfer across the interface must be included.

Acknowledgment 9. Conclusions Two related techniques which can be used to simulate diffusion controlled processes have been presented: random-walk and ‘‘the distribution function method’’. Random-walk can hardly be called a new method, but it has, hopefully, been shown that it has a broader range of applicability than it has historically been used for. The distribution function method is essentially a straightforward development of discrete randomwalks, when matter is considered as a continuum, but with different properties than more conventional ‘‘continuum methods’’. Both methods have in common that they are explicit and very easy to program. On the other hand an explicit technique means that stiff problems would be awkward to simulate; the techniques are multi- but not all-purpose. The method based on discrete point carriers of mass is very well suited for some types of 3D problems, especially when certain simplifying assumptions can be applied. A recommended article is

Funding from the Swedish Research Council is gratefully acknowledged. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]

˚ gren, Scripta Mater. 49 (2003) 521. H. Larsson, J. A ˚ gren, Scripta Mater. 51 (2004) 137. H. Larsson, J. A H. Strandlund, H. Larsson, Acta Mater. 52 (2004) 4695. J. Bardeen, Phys. Rev. 76 (1949) 1403. A.D. Le Claire, Phil. Mag. 3 (1958) 921. J.R. Manning, Phys. Rev. 124 (1961) 470. C.W. Gardiner, Handbook of Stochastic Methods, second ed., Springer, Berlin, Germany, 1985. L.S. Darken, Trans. AIME 180 (1948) 430. Matlab is a trademark of the Mathworks, Inc. J.-O. Andersson, T. Helander, L. Ho¨glund, P. Shi, B. Sundman, Calphad (26) (2002) 273. H. Strandlund, Master thesis, Trita-Na-E0050, KTH, Stockholm, Sweden, 2000. M. Hillert, Acta Mater. 47 (1999) 4481. ˚ gren, J. A. Borgenstam, A. Engstro¨m, L. Ho¨glund, J. A Phase Equil. 21 (2000) 269. J.A. Nesbitt, R.W. Heckel, Metall. Trans. A 18 (1987) 2061. ˚ gren, Acta Mater. 49 (2001) 3821. M. Schwind, J. A