- Email: [email protected]
Microstructure dependence of diffusional transport
Computational Materials Science 20 (2001) 37±47
www.elsevier.com/locate/commatsci
Microstructure dependence of diusional transport Jingzhi Zhu a,*, Long-Qing Chen a, Jie Shen b, Veena Tikare c a
Department of Materials Science and Engineering, 231 Steidle Building, Penn State University, University Park, PA 16802, USA b Department of Mathematics, Penn State University, University Park, PA 16802, USA c Materials Modeling and Simulation, Sandia National Laboratory, Albuquerque, NM 87185-1411, USA Received 2 February 2000; received in revised form 10 May 2000; accepted 30 May 2000
Abstract A simple and eective numerical method is proposed for simulating the temporal diusive mass transport process through a microstructure with arbitrary complexity described by a phase-®eld approach. The mass diusion through a given microstructure is modeled by a diusion equation with a variable diusion coecient, which is solved by an ecient and accurate semi-implicit spectral method. It is shown that it is possible to extract the eective diusion coecient for any given microstructure from the temporal concentration pro®les. The method is used to simulate the grain boundary diusion in a single-phase polycrystalline grain structure and the heterogeneous diusion in a twophase microstructure with dierent diusion coecient in each phase. Results are compared with existing analytical theories and computer simulations. Ó 2001 Elsevier Science B.V. All rights reserved. PACS: 66.30.-h; 02.60.Cb; 66.10.Cb Keywords: Diusion; Microstructures; Phase-®eld; Spectral method; Grain boundary diusion; Eective diusion coecient
1. Introduction Diusion in a solid with a multi-phase microstructure or with defects such as grain boundaries and dislocations can be highly heterogeneous. For example, in a polycrystalline grain structure, the diusion coecients along a grain boundary and through a bulk lattice can be orders of magnitude dierent. As a result, although the volume fraction of grain boundaries is small, grain boundary diffusion may dominate a number of processes such as grain growth, micro-electronics device failure,
* Corresponding author. Tel.: +1-814-863-9957; fax: +1-814867-0476. E-mail address: [email protected] (J. Zhu).
and sintering [1]. Macroscopically, the eective diusion coecient of a grain structure depends on the relative magnitudes of the grain boundary and bulk diusion coecients, the average grain size and grain boundary width. Although it is possible to derive analytical solutions for very special microstructures, it is impossible to derive the eective properties of a microstructure exactly. Therefore, most of the statistical theories are concerned with the property bounds [2]. Recently, there has been increasing interest in computing the eective properties from a microstructure either from a digitization of an experimental microstructure or from a microstructure simulation model [3]. However, most of these approaches are based on steady-state equations and thus cannot describe the temporal mass diusive transport behavior.
0927-0256/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 7 - 0 2 5 6 ( 0 0 ) 0 0 1 2 3 - 3
38
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
The purpose of this paper is to develop a simple and eective numerical method to simulate the temporal diusive transport through microstructures with arbitrary complexity. A microstructure is described by ®eld variables and generated by the phase-®eld method. A semi-implicit Fourier± Chebyshev spectral method was proposed to eciently solve the diusion equation with a variable diusion coecient. It will be shown that, using this method, both the temporal diusion pro®le through a microstructure and the macroscopic eective diusivity can be obtained. It is used to simulate the diusion through a single-phase grain structure and two-phase microstructures with different diusive properties in each phase. 2. Computer modeling 2.1. Microstructure description To model the diusion through an arbitrary microstructure, we used a set of space- and timedependent ®eld variables to describe the microstructure, similar to the phase-®eld approach [4]. For example, the concentration dierence throughout the microstructure is described by a concentration ®eld C r. A single-phase polycrystalline microstructure can be described by a set of order-parameter ®eld variables, g1 r, g2 r,. . ., gp r where p is the number of grain orientations in the system [5]. To simplify the problem, we studied diusion in a static microstructure which remained unchanged during the diusion process. The solid state reactions between the diusant atoms and the microstructure such as segregation, which may occur in some practical applications, are not considered in this work. Therefore, only temporal and spatial evolution of one variable, the concentration of diusing species, is needed to investigate the diffusive transport in a microstructure described by a set of static ®eld variables. The spatial dependence of diusivity is introduced through its dependence on the ®eld variables C r or g r, i.e., we describe the diusion coecient D as a function of C or g written as D f C; g1 ; g2 ; . . . ; gp . For example, in a single-phase polycrystalline material, the diu-
sion coecient throughout the system can P be expressed in a scaled form as D D 1 ÿ a pi1 g2i , and a are positive constants. The values where P D 2 of gi in the grain bulk are higher than those in the grain boundaries, and therefore the dierence in grain boundary and bulk diusion coecients can be easily described by this simple equation. Although the choice of this equation was quite arbitrary, we believed it would not aect the general results of the diusive transport process. One important advantage is the fact that it avoids the speci®cation of boundary conditions at the interfaces between dierent regions with dierent diffusivites. Since there was no microstructural P 2 evolution in our study, the variable gi was only a function of spatial coordinates r, and thus the diusion coecient was mapped to each lattice point as D r. It should also be pointed out that the ®eld variables are not required to be diuse as they are in the phase-®eld model. 2.2. Numerical method Neglecting the segregation of diusing species, the temporal diusive transport process through the microstructure can be described by Fick's second law written as oC r; t r D rrC r; t; ot
1
where C is the concentration of diusant, t the time, r the spatial coordinate, and D r is the diffusion coecient distribution in a given microstructure. Our numerical simulation was performed in two dimensions with a periodic boundary condition in one-dimension and a ®xed boundary condition in the other which has a ®xed diusant source and sink. The initial concentration of diusant throughout the system is zero. Analytical solutions to the above equation are only possible with highly idealized assumptions about the diusion coecient distribution D r. For studying a distribution with arbitrary complexity, numerical simulations have to be employed. Most of the existing numerical simulations for this problem were performed using explicit forward Euler method in time and ®nite-dierence in space. To maintain the stability of the scheme
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
and to achieve high accuracy, the time step and spatial grid size have to be very small, which seriously limit the system size and time duration of a simulation. Spectral methods, which were widely used in the ®eld of ¯uid dynamics [6], oered us an accurate and ecient alternative to solve partial dierential equations. Recently, Chen and Shen [7] used a semi-implicit Fourier spectral method for solving the time-dependent Ginzburg±Landau and Cahn±Hilliard equations. They showed that the semi-implicit Fourier spectral method is signi®cantly more ecient and accurate than a conventional explicit ®nite dierence method. Consider the forward Euler method for Eq. (1) C n1 ÿ C n r DrC n : Dt
2
The above scheme with a ®nite dierence or Fourier approximation in space has a severe time step constraint of the form Dt 6 k=N 2 [6], where k is a constant and N is the number of lattice points to be used in each space dimension. The constraint becomes Dt 6 k=N 4 [6] when the problem has a ®xed boundary and a spectral approximation is used in space. To avoid such a severe restriction, we consider a semi-implicit treatment. The idea is to split the variable diusion coecient D r into A and D r ÿ A, where A is a suitable constant, and treat them separately to obtain C n1 ÿ C n ÿ Ar2 C n1 ÿr A ÿ DrC n : Dt
3
write u C n1 , a 1= ADt, and f C n = ADtÿ r 1 ÿ D=ArC n , and then Eq. (3) can be simpli®ed as au ÿ r2 u f ;
4
with a periodic boundary condition in one direction and a ®xed boundary condition in the other. At each time step, this equation can be solved accurately and eciently by a Fourier±Chebyshev Galerkin algorithm which we brie¯y describe below. Expanding u and f in Fourier series in the x direction in which a periodic boundary P condition m is applied, and substituting u x; y 1 mÿ1 u y eimx (same for f) in Eq. (4), we obtain
39
a m2 um ÿ umyy f m ; 5
m 0; 1; 2; . . . ; um 1 0:
Then, we solve the above equation for each m by using the Chebyshev±Galerkin algorithm developed by Shen [8]. More precisely, let Tk y be the kth degree Chebyshev polynomial, we look for an approximation of um in the form umN P N ÿ2 m ^k /k y where /k y Tk y ÿ Tk2 y are k0 u the Galerkin basis functionsR satisfying /k 1 0. 1 Let us denote u; vx ÿ1 uvx dy where x 2 ÿ 1=2 1 ÿ y is the Chebyshev weight function. N ÿ2 Then, the coecients f^ umk gk0 (hence the approxm imation uN ) are determined by the weighted variational formulation for Eq. (5) a m2 umN ; /k x ÿ umN yy ; /k x IN f m ; /k x ;
k 0; 1; 2; . . . ; N ÿ 2;
6
where IN is a polynomial interpolation operator based on the Chebyshev±Gauss±Lobatto (CGL) points: yj cos jp=N ; j 0; 1; . . . ; N . More precisely, for any continuous function g on ÿ1; 1, IN g is the unique polynomial of degree less or equal than N such that IN g yj g yj ;
j 0; 1; . . . ; N :
Furthermore, the transform between the function values ff yj gNj0 and the Chebyshev coecients of IN f can be evaluated in O N log N operations by using the fast Fourier transform (FFT). Now, we describe how to determine the coeN ÿ2 cients f^ umk gk0 from Eq. (6). Let us denote Z 1 T IN f /k x dy; f f0 ; f1 ; . . . ; fN ÿ2 ; fk ÿ1 Z 1 /00j /k x dy; S skj k;j0;1;...;N ÿ2 ; skj ÿ ÿ1 Z 1 /j /k x dy; B bkj k;j0;1;...;N ÿ2 : bkj ÿ1
Then, it is easy to see that Eq. (6) is equivalent to the following linear system u f; a m2 B S
7
40
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47 T
where u ^ um0 ; u^m1 ; . . . ; u^mN ÿ2 . It is shown in [8] that the entries of the matrices S and B are given by:
3. Results and discussion
skk 2p k 1 k 2;
By numerically solving the diusion equation with a variable coecient, we can obtain the temporal evolution of the concentration pro®le through any arbitrary microstructure. The ®rst example is the grain boundary diusion, which has been studied for many years since the early 50s. For an idealized grain boundary geometry, analytical solutions are available for the concentration pro®les [9±11]. These solutions can be employed to analyze the grain boundary diusion coecients from experimental measurements [12]. In recent years, some numerical models were developed to overcome the simple geometry and topology restrictions. Gui et al. [13] analyzed the grain boundary diusion process in thin ®lms by transmission-line-matrix modeling . More recently, Swiler et al. [14,15] investigated some important heterogeneous diusion eects in polycrystalline microstructures obtained from Potts model simulations. We ®rst tested our approach by studying the diusion through an ideal grain boundary described by the Fisher's model where the analytical results were available [9]. The grain boundary is treated as a high-diusivity, semi-in®nite isotropic slab of uniform width, embedded in a low-diusivity, semi-in®nite perfect crystal which is normal to the surface that carries the diusant shown in Fig. 1(a). However, in numerical simulations we can approximate a grain boundary as a ®nite slab embedded in a ®nite crystal. This approximation can represent a semi-in®nite system very well if the grain boundary thickness d is far less than the system size L (d L). A simple way to describe such a grain boundary in Fisher's model is to choose two ®eld variables g1 , g2 . In grain 1, g1 1, g2 0. In grain 2 g1 0, g2 1. The phase ®eld model applied to describe the microstructure usually involves a diusive interface, where the ®eld variables g1 , g2 change continuously from 0 to 1 across the grain boundary, shown in Fig. 1(b). The solid line in the Fig. 1(b) plotted the diusion coecient D across the grain boundary described by D D 1 ÿ a g21 g22 where D 1
k 0; 1; . . . ; N ÿ 2;
skj 4p k 1; j k 2; k 4; k 6; . . . ; skj 0; j < k or j k odd; and bkk p;
k 6 0; and p bkj ÿ ; j k 2; 2 bkj 0; otherwise:
b00 2p;
Although S is not sparse but it is a special upper-triangular matrix whose non-zero elements (except the diagonal element) at each row are all the same. Therefore, one can design a special Gaussian elimination process which takes into account this special structure to solve Eq. (7) in O N operations (more precisely, at a cost comparable to that of solving a penda-diagonal matrix). We refer to [8] for more details on this matter. Thus, the overall operation counts for solving Eq. (5) is O N log N for each m. Hence, assuming M Fourier coecients are retained in the Fourier expansion, the total cost of the Fourier± Chebyshev Galerkin method for solving Eq. (4) is O NM log NM which is quasi optimal with respect to the number of points (lattices) used. Note that the total cost is only slightly more than O NM operations which are required by the conventional ®nite dierence/®nite element methods, but the Fourier±Chebyshev Galerkin method is capable of producing much more accurate results using signi®cantly less number of points (lattices). On the other hand, the semiimplicit treatment will allow us to use much larger time steps than that is allowed by an explicit method (about two orders of magnitude larger from our numerical experiments). Hence, the resulting semi-implicit Fourier±Chebyshev Galerkin method is an extremely ecient method for solving diusion equation with a variable diusion coecient.
3.1. Grain boundary diusion
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
41
The kinetics of grain boundary diusion are usually classi®ed to type A, B, and C according to the parameters of the model, such as the ratio of grain boundary diusion coecient to the bulk lattice diusion coecient, the diusion time, and the ratio of grain boundary width to the grain size [16]. We have carried out simulations in all three regimes and compared the concentration pro®le at dierent times obtained from our simulations with the analytical solutions in each regime. The A regime is for long diusion time t and small grain size d with the condition: Dv t1=2 d. The concentration pro®les in type A kinetics appear to obey Fick's law for a homogeneous system with an eective coecient Deff which can be written as y hCi erfc p ; 8 2 Deff t
Fig. 1. Diusion through an idealized grain boundary: (a) Fisher's model for grain boundary diusion; (b) ®eld variables g1 , g2 and diusion coecient D across a grain boundary in a phase-®eld description; (c) an isoconcentration contour plot with eight levels through a grain boundary after a certain diffusion time, D D 1 ÿ a g21 g22 , D 1 and a 0:999.
and a 0:999. Fig. 1(c) shows an isoconcentration contour plot through a grain boundary after a certain diusion time. The diusion coecient within the grain boundary being higher than in the crystal, the diusion penetrated deeper along the grain boundary than in the grain bulk. Consequently, it started leaking through the two walls of the grain boundary into the grain bulk.
where Deff fDgb 1 ÿ f Dv and f is the volume fraction of grain boundaries. The comparison of the average concentration hCi at dierent times in A regime was shown in Fig. 2(a) where our numerical results (shown in symbols) agree very well with analytical solutions from Eq. (8) (shown in lines). Type C kinetics is the opposite of type A kinetics where volume diusion in the grain bulk is 1=2 negligible, i.e., Dv t d. Therefore the diusion is only one-dimensional in Fig. 1(a) (along y direction). With the de®ned boundary condition, the concentration along the grain boundary in Fig. 1(a) is ! y 9 C y; t erfc p : 2 Dgb t Fig. 2(b) showed a good agreement of the numerically calculated and analytical results of the concentration along the grain boundary in regime C. Type B kinetics is of the intermediate case between B regime and C regime with the condition d Dv t1=2 < d=2. Most experimental conditions fall into this regime. The analytical solution in the B regime was obtained by Whipple in a complex form [10]. The distribution of concentration through the microstructure is described by
42
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
Fig. 2. Comparison of the concentration pro®le at dierent diusion distance y for a grain boundary diusion at dierent diusion times: (a) type A kinetics, hCi vs y; (b) type C kinetics, C y; t vs y; (c) type B kinetics, hCi vs y 6=5 .
Z D g 1 g2 p exp ÿ 2 4r 2 p 1 r3=2 " r # 1 Dÿ1 rÿ1 n dr: erfc 2 Dÿr b
C erfc
g
10
Several dimensionless parameters pare de®ned in the above equation as: g y= Dv t , D p Dgb =Dv , p b D ÿ 1d= 2 Dv t, and n x ÿ d=2= Dv t. Eq. (10) was evaluated with all parameters the same as those chosen in our simulations. For example, the thickness of the grain boundary d should be determined from the diuse-interface ®eld variable pro®le [7]. Fig. 2(c) shows one example of the average concentration hCi as a function of y 6=5 at 3 dierent time steps. We observed a good agreement between the simulation results and the analytical Whipple's solution, especially at long times in a type-B kinetics where the Whipple's solution is valid [17]. However, at short times t 2000 there is a little dierence between the simulated hCi and the calculated hCi because 1=2 the condition d Dv t is not well satis®ed at short times. The plot of the logarithm of the average concentration versus y 6=5 was the basis for experimental measurements of grain boundary diusion coecients [12]. Based on the simple examples discussed above, the numerical model works quite well in all three regimes of grain boundary diusion. Although for the idealized Fisher's model, analytical solutions are possible, it is usually not the case for polycrystals, particularly if one is interested in the spatial distribution of diusing species in the microstructure. The numerical model proposed in this work can be applied to diusion through a polycrystal without any further complication. Fig. 3 shows an example of the temporal evolution of concentration during diusion through a polycrystalline grain structure. The grain microstructure, which was generated from a phase ®eld simulation of grain growth [5], was mapped to the x-y plane. Heterogeneous diusion eect is clear that the concentration values in the grain boundaries are much higher than those in the grain bulk. As can be seen from the concentration distribution in Fig. 3(a) (t 50,000) and Fig. 3(b) (t 500,000), more species diused into the microstructure at longer diusion times. The diusion coecient in the P grain microstructure is written as D D 1 ÿ a i g2i where D and a are positive constants. Dierent ratio of grain boundary diusion coecient Dgb and bulk diusion coecient Dv can be easily obtained by
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
43
Fig. 3. Temporal evolution of concentration through a polycrystalline microstructure at two dierent times, D P D 1 ÿ a i g2i , D 1 and a 0:9999: (a) t 50,000; (b) t 500,000.
Fig. 4. Concentration pro®le in a polycrystal for two ratios ofPDgb =Dv at the same diusion time t 10,000, D D 1ÿ a i g2i , D 1: (a) a 0:8; (b) a 0:98.
adjusting the value of a. Fig. 4 shows the concentration pro®le through a polycrystalline material for two ratios of Dgb =Dv at the same diusion time. The ratio in Fig. 4(b) (Dgb =Dv 20) is higher than in Fig. 4(a) (Dgb =Dv 2:5), and thus the concentration of the diusing atoms in the grain boundary diers more considerably from that in the grain bulk.
conductivity and diusion coecient, were often approximated using the so-called mixture rules [18,19] X n Vi Kin ; ÿ1 6 n 6 1; 11 Keff
3.2. Eective diusion coecients of microstructures The eective properties of a microstructure, such as electrical conductivity, dielectric constant, magnetic permeability, elastic moduli, thermal
i
where Keff is the eective property of the microstructure, Ki and Vi are the corresponding property and volume fraction of the ith phase. More rigorous theories can predict the upper and lower bounds for the properties [2,18]. Recently, there has been increasing interest in directly calculating eective properties of digitized experimental microstructures or microstructures generated from
44
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
computer models [3]. Most of these numerical approaches aimed at approximate solutions of the steady-state equations using ®nite element or ®nite dierent methods [3]. It is easy to show that one can derive the effective properties of a microstructure from the temporal concentration pro®les by solving the non-steady-state diusion equation. We employed two approaches to extract the eective diusion coecient Deff for any given microstructure. One is studying the average concentration pro®le hCi as a function of diusion distance y at short diusion times. Assuming Eq. (8) describes the relationship between hCi and y, we can approximately calculate Deff by ®tting diusion data to Eq. (8). Another approach is to derive Deff from the steady-state solution of the non-steady-state diusion equation at longer intervals of time. The ¯ux J for long diusion time is related to the concentration gradient by Fick's ®rst law, J ÿDeff rC:
12
The total ¯ux J can be obtained easily from the diusion coecient distribution D r and the concentration distribution C r through the microstructure at the steady state. rC in Eq. (12) is determined by the applied concentration gradient. Our numerical work shows Deff computed from Eqs. (8) and (12) has a dierence less than 2%, which can satisfy most of the practical applications. One of the main advantages of this approach is that it produces both the time-dependent
concentration pro®les in a microstructure and the eective properties of the microstructure. In addition, it avoids any complicated boundary conditions and thus allows one to study the diusion transport through any arbitrary microstructures. As test examples, we ®rst computed Deff for two of the simplest possible microstructures. The system is composed of parallel slabs of A and B with dierent diusion coecients DA and DB , shown in Fig. 5. In Fig. 5(a), the ¯ux is parallel to the plane of the slabs, and the arrangement is equivalent to a parallel electrical circuit. Each slab has the same concentration gradient, and most of the ¯ux is through the component with higher diusivity. Deff is given by Deff DA VA DB VB , where VA and VB are the volume fraction of each component. This corresponds to n 1 in Eq. (11). In contrast, for diusion ¯ux perpendicular to the plane of the slabs in Fig. 5(b), it is equivalent to a series electrical circuits. The ¯ux through each slab is equal, but the concentration gradients are dierent. In this case, n ÿ1 in Eq. (11) and Deff DA DB = DA VB DB VA . Fig. 6 shows the comparison of calculated and theoretical Deff as a function of volume fraction for these two parallel slab microstructures. Slab arrangements in Fig. 5(a) and (b) are indicated by jj and respectively. As can be seen from the graph, the eective diusion coecients computed agree well with those of theoretical solutions for both cases, thus validating our methods. Generally Deff for the two parallel slab microstructures will cover the complete range of
Fig. 5. Idealized two-phase microstructures ± parallel slabs: (a) ¯ux parallel to the plane of slabs jj; (b) ¯ux perpendicular to the plane of slabs ( ).
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
Fig. 6. Computed and theoretical eective diusion coecient Deff as a function of volume fraction of phase A for parallel slab microstructures, DA 1.0, DB 0.2.
eective properties between the upper and lower bounds as shown in the solid line and dashed line in Fig. 6. It would be expected that for arbitrary microstructures, Deff was within this range. For a polycrystalline material, diusion was carried out as described in the previous section. In particular, we examined the eect of the average grain size d of a grain structure on the eective diusion coecient. The results are plotted as the ®lled circles shown in Fig. 7, where Deff;simulation is calculated from steady state concentration pro®le. Deff;simulation varies inversely with the average grain
Fig. 7. Eective diusion coecient Deff as a function of average grain size d for a single-phase polycrystalline microP structure, D 1 ÿ 0:999 i g2 . Deff;upper and Deff;low are the eective diusion coecients predicted by mixture rules.
45
size d. This result is in good agreement with that of sharp-interface simulations by Swiler et al. [14]. The driving force for grain growth is the reduction of grain boundary energy. As the grain size becomes larger, the volume fraction of grain boundaries decreases, resulting in a decrease in the eective diusion coecient. In Fig. 7, we also plotted the theoretical eective diusion coecients Deff;upper and Deff;low as the upper and low bounds of the eective diusion coecients which can be predicted by rule of mixing as Dgb fgb Dv 1 ÿ fgb and Dgb Dv =Dgb 1 ÿ fgb Dv fgb , respectively, where fgb is the volume fraction of grain boundaries. However, for a diuse-interface description of the microstructure, there is some ambiguity in determining the exact values of Dgb , Dv and fgb . Recognizing that Deff predicted by rule of mixing is essentially a P volume average, we calculated Deff;upper as 1=N r D r where N is the total number of latticePpoints. Similarly Deff;low was computed as N = 1=D r. As expected, for all the grain sizes studied, Deff;simulation is between Deff;upper and Deff;low , which are the eective diusion coecients for idealized grain structures where grain boundaries can be treated as parallel slabs embedded in the grains. The dierence between Deff;simulation and Deff;upper was explained by Swiler et al. [15] as a result of triple junctions between grains as bottlenecks to diusion. Their sharp-interface simulations showed that Deff;simulation 0:27Dgb d=d [14], where d is the grain boundary width and d=d is the approximated volume fraction of grain boundaries. As a ®nal example, diusion through two-phase microstructures was studied. The two-phase microstructures in Fig. 8 were obtained from our phase ®eld simulations of phase separation of an initially homogeneous system. The two-phase system forms an inter-connected or a dispersed microstructure at dierent volume fractions. Fig. 9 shows the change of Deff (obtained from the steady state concentration pro®le) as a function of volume fraction for these microstructures (plotted with circle symbols). As the volume fraction of phase A with high diusion coecient DA DB increases, the shape of A droplets changes, and Deff increases gradually. After VA reaching about 50%, there is a sharp increase in Deff when the
46
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
Fig. 8. Two-phase microstructures for dierent volume fraction of phase A (shown in white color) showing the percolation transition, (a)±(e) for VA 16%, 32%, 50%, 62%, and 87%, respectively.
Fig. 9. Eective diusion coecient Deff as a function of volume fraction. Circles were calculated from the microstructures in Fig. 8. Dashed line and solid line were plotted from Eq. (11) representing the upper and lower bounds of Deff . DA =DB 1000.
microstructure changes from isolated phase A in B to inter-connected A and B phases. The ratio of DA =DB is three orders of magnitude in Fig. 9. The dashed line and solid line represent the upper and low bounds of Deff determined from Eq. (11) as n 1 and ÿ1, respectively. It is usually observed that the larger the distinction between the properties of dierent phases in heterogeneous systems, the more obvious is the percolation transition [18]. Our numerical experiments have demonstrated that the Fourier±Chebyshev spectral method is very ecient for solving the time-dependent diffusion equation in which the diusivity distribution is a smooth function throughout the
microstructure, especially for microstructures described by the diuse-interface model. The excellent agreement between our diuse-interface simulation and sharp-interface analytical or numerical methods showed that diuse-interface model is very effective for describing diusional transport in mesoscale microstructures. In recent years, the diuse-interface phase-®eld approach has been recognized as a powerful mathematical model for simulating the mesoscale morphological pattern formation such as grain growth and Ostwald ripening. By coupling the phase-®eld modeling of microstructure evolution, the approach developed in this paper allows us to model the eective diusivity evolution in an evolving microstructure. More details will be reported in the near future.
4. Conclusion We have simulated the temporal mass diusion transport through an arbitrary microstructure based on a phase-®eld approach. Diusion coecients were calculated as a function of the order parameters describing a single-phase grain structure where grain boundary diusion was much higher than bulk diusion and a two-phase microstructure where the diusivities varied with the phase. The semi-implicit Fourier±Chebyshev spectral method was shown to be accurate and ecient in solving the variable diusion equations. The accuracy of the model was demonstrated by reproducing the classic diusion results of grain boundary diusion in all three kinetic regimes
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
(type A, B and C) as well as the parallel slabs geometries. We demonstrated that one can extract eective diusion coecient of a microstructure from the temporal diusion pro®les. Furthermore, the sharp interface results of Swiler et al. were reproduced to show that the eective diusivity varies inversely with grain size. This model was used to predict the eective diusivity of twophase percolating and non-percolating structures where each phase has very dierent diusivities. This approach can be also used to study similar transport problems such as heat transport by recognizing the similarity between mass diusion and thermal diusion equations. Acknowledgements The authors are grateful for the ®nancial support from the Sandia National Laboratory and the National Science Foundation under Grant No. DMR 96-33719 and DMS 9721413. References [1] I. Kaur, Y. Mishin, W. Gust, Fundamentals of Grain and Interphase Boundary Diusion, third ed., Wiley, Chichester, 1995.
47
[2] S. Torquato, Int. J. Solids Struct. 35 (1998) 2385. [3] E.J. Garboczi, D.P. Bentz, Constr. Building Mater. 10 (1996) 293. [4] L.Q. Chen, Y.Z. Wang, JOM 48 (1996) 13. [5] D.N. Fan, L.Q. Chen, S.P. Chen, Mater. Sci. Eng. A 238 (1997) 78. [6] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid Dynamics, Springer, Berlin, 1987. [7] L.Q. Chen, J. Shen, Comput. Phys. Commun. 108 (1998) 147. [8] J. Shen, SIAM J. Sci. Comput. 16 (1995) 74. [9] J.C. Fisher, J. Appl. Phys. 22 (1951) 74. [10] R.T.P. Whipple, Philos. Mag. 45 (1954) 1225. [11] T. Suzuoka, J. Phys. Soc. Jpn. 19 (1964) 839. [12] H.S. Levine, C.J. MacCallum, J. Appl. Phys. 31 (1960) 595. [13] X. Gui, S.K. Dew, M.J. Brett, J. Appl. Phys. 74 (1993) 7173. [14] T.P. Swiler, V. Tikare, E.A. Holm, Mater. Sci. Eng. A 238 (1997) 85. [15] T.P. Swiler, E.A. Holm, Ceram. Trans. 69 (1997) 203. [16] L.G. Harrison, Trans. Faraday Soc. 57 (1961) 1191. [17] J. Philebert, Atom Movements ± Diusion and Mass Transport in Solids, Les Editions de Physique, Les Ulis, 1991. [18] C.W. Nan, Prog. Mater. Sci. 37 (1993) 1. [19] L.E. Nielsen, Predicting the Properties of Mixtures: Mixture Rules in Science and Engineering, Marcel Dekker, New York, 1978.
www.elsevier.com/locate/commatsci
Microstructure dependence of diusional transport Jingzhi Zhu a,*, Long-Qing Chen a, Jie Shen b, Veena Tikare c a
Department of Materials Science and Engineering, 231 Steidle Building, Penn State University, University Park, PA 16802, USA b Department of Mathematics, Penn State University, University Park, PA 16802, USA c Materials Modeling and Simulation, Sandia National Laboratory, Albuquerque, NM 87185-1411, USA Received 2 February 2000; received in revised form 10 May 2000; accepted 30 May 2000
Abstract A simple and eective numerical method is proposed for simulating the temporal diusive mass transport process through a microstructure with arbitrary complexity described by a phase-®eld approach. The mass diusion through a given microstructure is modeled by a diusion equation with a variable diusion coecient, which is solved by an ecient and accurate semi-implicit spectral method. It is shown that it is possible to extract the eective diusion coecient for any given microstructure from the temporal concentration pro®les. The method is used to simulate the grain boundary diusion in a single-phase polycrystalline grain structure and the heterogeneous diusion in a twophase microstructure with dierent diusion coecient in each phase. Results are compared with existing analytical theories and computer simulations. Ó 2001 Elsevier Science B.V. All rights reserved. PACS: 66.30.-h; 02.60.Cb; 66.10.Cb Keywords: Diusion; Microstructures; Phase-®eld; Spectral method; Grain boundary diusion; Eective diusion coecient
1. Introduction Diusion in a solid with a multi-phase microstructure or with defects such as grain boundaries and dislocations can be highly heterogeneous. For example, in a polycrystalline grain structure, the diusion coecients along a grain boundary and through a bulk lattice can be orders of magnitude dierent. As a result, although the volume fraction of grain boundaries is small, grain boundary diffusion may dominate a number of processes such as grain growth, micro-electronics device failure,
* Corresponding author. Tel.: +1-814-863-9957; fax: +1-814867-0476. E-mail address: [email protected] (J. Zhu).
and sintering [1]. Macroscopically, the eective diusion coecient of a grain structure depends on the relative magnitudes of the grain boundary and bulk diusion coecients, the average grain size and grain boundary width. Although it is possible to derive analytical solutions for very special microstructures, it is impossible to derive the eective properties of a microstructure exactly. Therefore, most of the statistical theories are concerned with the property bounds [2]. Recently, there has been increasing interest in computing the eective properties from a microstructure either from a digitization of an experimental microstructure or from a microstructure simulation model [3]. However, most of these approaches are based on steady-state equations and thus cannot describe the temporal mass diusive transport behavior.
0927-0256/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 7 - 0 2 5 6 ( 0 0 ) 0 0 1 2 3 - 3
38
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
The purpose of this paper is to develop a simple and eective numerical method to simulate the temporal diusive transport through microstructures with arbitrary complexity. A microstructure is described by ®eld variables and generated by the phase-®eld method. A semi-implicit Fourier± Chebyshev spectral method was proposed to eciently solve the diusion equation with a variable diusion coecient. It will be shown that, using this method, both the temporal diusion pro®le through a microstructure and the macroscopic eective diusivity can be obtained. It is used to simulate the diusion through a single-phase grain structure and two-phase microstructures with different diusive properties in each phase. 2. Computer modeling 2.1. Microstructure description To model the diusion through an arbitrary microstructure, we used a set of space- and timedependent ®eld variables to describe the microstructure, similar to the phase-®eld approach [4]. For example, the concentration dierence throughout the microstructure is described by a concentration ®eld C r. A single-phase polycrystalline microstructure can be described by a set of order-parameter ®eld variables, g1 r, g2 r,. . ., gp r where p is the number of grain orientations in the system [5]. To simplify the problem, we studied diusion in a static microstructure which remained unchanged during the diusion process. The solid state reactions between the diusant atoms and the microstructure such as segregation, which may occur in some practical applications, are not considered in this work. Therefore, only temporal and spatial evolution of one variable, the concentration of diusing species, is needed to investigate the diffusive transport in a microstructure described by a set of static ®eld variables. The spatial dependence of diusivity is introduced through its dependence on the ®eld variables C r or g r, i.e., we describe the diusion coecient D as a function of C or g written as D f C; g1 ; g2 ; . . . ; gp . For example, in a single-phase polycrystalline material, the diu-
sion coecient throughout the system can P be expressed in a scaled form as D D 1 ÿ a pi1 g2i , and a are positive constants. The values where P D 2 of gi in the grain bulk are higher than those in the grain boundaries, and therefore the dierence in grain boundary and bulk diusion coecients can be easily described by this simple equation. Although the choice of this equation was quite arbitrary, we believed it would not aect the general results of the diusive transport process. One important advantage is the fact that it avoids the speci®cation of boundary conditions at the interfaces between dierent regions with dierent diffusivites. Since there was no microstructural P 2 evolution in our study, the variable gi was only a function of spatial coordinates r, and thus the diusion coecient was mapped to each lattice point as D r. It should also be pointed out that the ®eld variables are not required to be diuse as they are in the phase-®eld model. 2.2. Numerical method Neglecting the segregation of diusing species, the temporal diusive transport process through the microstructure can be described by Fick's second law written as oC r; t r D rrC r; t; ot
1
where C is the concentration of diusant, t the time, r the spatial coordinate, and D r is the diffusion coecient distribution in a given microstructure. Our numerical simulation was performed in two dimensions with a periodic boundary condition in one-dimension and a ®xed boundary condition in the other which has a ®xed diusant source and sink. The initial concentration of diusant throughout the system is zero. Analytical solutions to the above equation are only possible with highly idealized assumptions about the diusion coecient distribution D r. For studying a distribution with arbitrary complexity, numerical simulations have to be employed. Most of the existing numerical simulations for this problem were performed using explicit forward Euler method in time and ®nite-dierence in space. To maintain the stability of the scheme
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
and to achieve high accuracy, the time step and spatial grid size have to be very small, which seriously limit the system size and time duration of a simulation. Spectral methods, which were widely used in the ®eld of ¯uid dynamics [6], oered us an accurate and ecient alternative to solve partial dierential equations. Recently, Chen and Shen [7] used a semi-implicit Fourier spectral method for solving the time-dependent Ginzburg±Landau and Cahn±Hilliard equations. They showed that the semi-implicit Fourier spectral method is signi®cantly more ecient and accurate than a conventional explicit ®nite dierence method. Consider the forward Euler method for Eq. (1) C n1 ÿ C n r DrC n : Dt
2
The above scheme with a ®nite dierence or Fourier approximation in space has a severe time step constraint of the form Dt 6 k=N 2 [6], where k is a constant and N is the number of lattice points to be used in each space dimension. The constraint becomes Dt 6 k=N 4 [6] when the problem has a ®xed boundary and a spectral approximation is used in space. To avoid such a severe restriction, we consider a semi-implicit treatment. The idea is to split the variable diusion coecient D r into A and D r ÿ A, where A is a suitable constant, and treat them separately to obtain C n1 ÿ C n ÿ Ar2 C n1 ÿr A ÿ DrC n : Dt
3
write u C n1 , a 1= ADt, and f C n = ADtÿ r 1 ÿ D=ArC n , and then Eq. (3) can be simpli®ed as au ÿ r2 u f ;
4
with a periodic boundary condition in one direction and a ®xed boundary condition in the other. At each time step, this equation can be solved accurately and eciently by a Fourier±Chebyshev Galerkin algorithm which we brie¯y describe below. Expanding u and f in Fourier series in the x direction in which a periodic boundary P condition m is applied, and substituting u x; y 1 mÿ1 u y eimx (same for f) in Eq. (4), we obtain
39
a m2 um ÿ umyy f m ; 5
m 0; 1; 2; . . . ; um 1 0:
Then, we solve the above equation for each m by using the Chebyshev±Galerkin algorithm developed by Shen [8]. More precisely, let Tk y be the kth degree Chebyshev polynomial, we look for an approximation of um in the form umN P N ÿ2 m ^k /k y where /k y Tk y ÿ Tk2 y are k0 u the Galerkin basis functionsR satisfying /k 1 0. 1 Let us denote u; vx ÿ1 uvx dy where x 2 ÿ 1=2 1 ÿ y is the Chebyshev weight function. N ÿ2 Then, the coecients f^ umk gk0 (hence the approxm imation uN ) are determined by the weighted variational formulation for Eq. (5) a m2 umN ; /k x ÿ umN yy ; /k x IN f m ; /k x ;
k 0; 1; 2; . . . ; N ÿ 2;
6
where IN is a polynomial interpolation operator based on the Chebyshev±Gauss±Lobatto (CGL) points: yj cos jp=N ; j 0; 1; . . . ; N . More precisely, for any continuous function g on ÿ1; 1, IN g is the unique polynomial of degree less or equal than N such that IN g yj g yj ;
j 0; 1; . . . ; N :
Furthermore, the transform between the function values ff yj gNj0 and the Chebyshev coecients of IN f can be evaluated in O N log N operations by using the fast Fourier transform (FFT). Now, we describe how to determine the coeN ÿ2 cients f^ umk gk0 from Eq. (6). Let us denote Z 1 T IN f /k x dy; f f0 ; f1 ; . . . ; fN ÿ2 ; fk ÿ1 Z 1 /00j /k x dy; S skj k;j0;1;...;N ÿ2 ; skj ÿ ÿ1 Z 1 /j /k x dy; B bkj k;j0;1;...;N ÿ2 : bkj ÿ1
Then, it is easy to see that Eq. (6) is equivalent to the following linear system u f; a m2 B S
7
40
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47 T
where u ^ um0 ; u^m1 ; . . . ; u^mN ÿ2 . It is shown in [8] that the entries of the matrices S and B are given by:
3. Results and discussion
skk 2p k 1 k 2;
By numerically solving the diusion equation with a variable coecient, we can obtain the temporal evolution of the concentration pro®le through any arbitrary microstructure. The ®rst example is the grain boundary diusion, which has been studied for many years since the early 50s. For an idealized grain boundary geometry, analytical solutions are available for the concentration pro®les [9±11]. These solutions can be employed to analyze the grain boundary diusion coecients from experimental measurements [12]. In recent years, some numerical models were developed to overcome the simple geometry and topology restrictions. Gui et al. [13] analyzed the grain boundary diusion process in thin ®lms by transmission-line-matrix modeling . More recently, Swiler et al. [14,15] investigated some important heterogeneous diusion eects in polycrystalline microstructures obtained from Potts model simulations. We ®rst tested our approach by studying the diusion through an ideal grain boundary described by the Fisher's model where the analytical results were available [9]. The grain boundary is treated as a high-diusivity, semi-in®nite isotropic slab of uniform width, embedded in a low-diusivity, semi-in®nite perfect crystal which is normal to the surface that carries the diusant shown in Fig. 1(a). However, in numerical simulations we can approximate a grain boundary as a ®nite slab embedded in a ®nite crystal. This approximation can represent a semi-in®nite system very well if the grain boundary thickness d is far less than the system size L (d L). A simple way to describe such a grain boundary in Fisher's model is to choose two ®eld variables g1 , g2 . In grain 1, g1 1, g2 0. In grain 2 g1 0, g2 1. The phase ®eld model applied to describe the microstructure usually involves a diusive interface, where the ®eld variables g1 , g2 change continuously from 0 to 1 across the grain boundary, shown in Fig. 1(b). The solid line in the Fig. 1(b) plotted the diusion coecient D across the grain boundary described by D D 1 ÿ a g21 g22 where D 1
k 0; 1; . . . ; N ÿ 2;
skj 4p k 1; j k 2; k 4; k 6; . . . ; skj 0; j < k or j k odd; and bkk p;
k 6 0; and p bkj ÿ ; j k 2; 2 bkj 0; otherwise:
b00 2p;
Although S is not sparse but it is a special upper-triangular matrix whose non-zero elements (except the diagonal element) at each row are all the same. Therefore, one can design a special Gaussian elimination process which takes into account this special structure to solve Eq. (7) in O N operations (more precisely, at a cost comparable to that of solving a penda-diagonal matrix). We refer to [8] for more details on this matter. Thus, the overall operation counts for solving Eq. (5) is O N log N for each m. Hence, assuming M Fourier coecients are retained in the Fourier expansion, the total cost of the Fourier± Chebyshev Galerkin method for solving Eq. (4) is O NM log NM which is quasi optimal with respect to the number of points (lattices) used. Note that the total cost is only slightly more than O NM operations which are required by the conventional ®nite dierence/®nite element methods, but the Fourier±Chebyshev Galerkin method is capable of producing much more accurate results using signi®cantly less number of points (lattices). On the other hand, the semiimplicit treatment will allow us to use much larger time steps than that is allowed by an explicit method (about two orders of magnitude larger from our numerical experiments). Hence, the resulting semi-implicit Fourier±Chebyshev Galerkin method is an extremely ecient method for solving diusion equation with a variable diusion coecient.
3.1. Grain boundary diusion
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
41
The kinetics of grain boundary diusion are usually classi®ed to type A, B, and C according to the parameters of the model, such as the ratio of grain boundary diusion coecient to the bulk lattice diusion coecient, the diusion time, and the ratio of grain boundary width to the grain size [16]. We have carried out simulations in all three regimes and compared the concentration pro®le at dierent times obtained from our simulations with the analytical solutions in each regime. The A regime is for long diusion time t and small grain size d with the condition: Dv t1=2 d. The concentration pro®les in type A kinetics appear to obey Fick's law for a homogeneous system with an eective coecient Deff which can be written as y hCi erfc p ; 8 2 Deff t
Fig. 1. Diusion through an idealized grain boundary: (a) Fisher's model for grain boundary diusion; (b) ®eld variables g1 , g2 and diusion coecient D across a grain boundary in a phase-®eld description; (c) an isoconcentration contour plot with eight levels through a grain boundary after a certain diffusion time, D D 1 ÿ a g21 g22 , D 1 and a 0:999.
and a 0:999. Fig. 1(c) shows an isoconcentration contour plot through a grain boundary after a certain diusion time. The diusion coecient within the grain boundary being higher than in the crystal, the diusion penetrated deeper along the grain boundary than in the grain bulk. Consequently, it started leaking through the two walls of the grain boundary into the grain bulk.
where Deff fDgb 1 ÿ f Dv and f is the volume fraction of grain boundaries. The comparison of the average concentration hCi at dierent times in A regime was shown in Fig. 2(a) where our numerical results (shown in symbols) agree very well with analytical solutions from Eq. (8) (shown in lines). Type C kinetics is the opposite of type A kinetics where volume diusion in the grain bulk is 1=2 negligible, i.e., Dv t d. Therefore the diusion is only one-dimensional in Fig. 1(a) (along y direction). With the de®ned boundary condition, the concentration along the grain boundary in Fig. 1(a) is ! y 9 C y; t erfc p : 2 Dgb t Fig. 2(b) showed a good agreement of the numerically calculated and analytical results of the concentration along the grain boundary in regime C. Type B kinetics is of the intermediate case between B regime and C regime with the condition d Dv t1=2 < d=2. Most experimental conditions fall into this regime. The analytical solution in the B regime was obtained by Whipple in a complex form [10]. The distribution of concentration through the microstructure is described by
42
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
Fig. 2. Comparison of the concentration pro®le at dierent diusion distance y for a grain boundary diusion at dierent diusion times: (a) type A kinetics, hCi vs y; (b) type C kinetics, C y; t vs y; (c) type B kinetics, hCi vs y 6=5 .
Z D g 1 g2 p exp ÿ 2 4r 2 p 1 r3=2 " r # 1 Dÿ1 rÿ1 n dr: erfc 2 Dÿr b
C erfc
g
10
Several dimensionless parameters pare de®ned in the above equation as: g y= Dv t , D p Dgb =Dv , p b D ÿ 1d= 2 Dv t, and n x ÿ d=2= Dv t. Eq. (10) was evaluated with all parameters the same as those chosen in our simulations. For example, the thickness of the grain boundary d should be determined from the diuse-interface ®eld variable pro®le [7]. Fig. 2(c) shows one example of the average concentration hCi as a function of y 6=5 at 3 dierent time steps. We observed a good agreement between the simulation results and the analytical Whipple's solution, especially at long times in a type-B kinetics where the Whipple's solution is valid [17]. However, at short times t 2000 there is a little dierence between the simulated hCi and the calculated hCi because 1=2 the condition d Dv t is not well satis®ed at short times. The plot of the logarithm of the average concentration versus y 6=5 was the basis for experimental measurements of grain boundary diusion coecients [12]. Based on the simple examples discussed above, the numerical model works quite well in all three regimes of grain boundary diusion. Although for the idealized Fisher's model, analytical solutions are possible, it is usually not the case for polycrystals, particularly if one is interested in the spatial distribution of diusing species in the microstructure. The numerical model proposed in this work can be applied to diusion through a polycrystal without any further complication. Fig. 3 shows an example of the temporal evolution of concentration during diusion through a polycrystalline grain structure. The grain microstructure, which was generated from a phase ®eld simulation of grain growth [5], was mapped to the x-y plane. Heterogeneous diusion eect is clear that the concentration values in the grain boundaries are much higher than those in the grain bulk. As can be seen from the concentration distribution in Fig. 3(a) (t 50,000) and Fig. 3(b) (t 500,000), more species diused into the microstructure at longer diusion times. The diusion coecient in the P grain microstructure is written as D D 1 ÿ a i g2i where D and a are positive constants. Dierent ratio of grain boundary diusion coecient Dgb and bulk diusion coecient Dv can be easily obtained by
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
43
Fig. 3. Temporal evolution of concentration through a polycrystalline microstructure at two dierent times, D P D 1 ÿ a i g2i , D 1 and a 0:9999: (a) t 50,000; (b) t 500,000.
Fig. 4. Concentration pro®le in a polycrystal for two ratios ofPDgb =Dv at the same diusion time t 10,000, D D 1ÿ a i g2i , D 1: (a) a 0:8; (b) a 0:98.
adjusting the value of a. Fig. 4 shows the concentration pro®le through a polycrystalline material for two ratios of Dgb =Dv at the same diusion time. The ratio in Fig. 4(b) (Dgb =Dv 20) is higher than in Fig. 4(a) (Dgb =Dv 2:5), and thus the concentration of the diusing atoms in the grain boundary diers more considerably from that in the grain bulk.
conductivity and diusion coecient, were often approximated using the so-called mixture rules [18,19] X n Vi Kin ; ÿ1 6 n 6 1; 11 Keff
3.2. Eective diusion coecients of microstructures The eective properties of a microstructure, such as electrical conductivity, dielectric constant, magnetic permeability, elastic moduli, thermal
i
where Keff is the eective property of the microstructure, Ki and Vi are the corresponding property and volume fraction of the ith phase. More rigorous theories can predict the upper and lower bounds for the properties [2,18]. Recently, there has been increasing interest in directly calculating eective properties of digitized experimental microstructures or microstructures generated from
44
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
computer models [3]. Most of these numerical approaches aimed at approximate solutions of the steady-state equations using ®nite element or ®nite dierent methods [3]. It is easy to show that one can derive the effective properties of a microstructure from the temporal concentration pro®les by solving the non-steady-state diusion equation. We employed two approaches to extract the eective diusion coecient Deff for any given microstructure. One is studying the average concentration pro®le hCi as a function of diusion distance y at short diusion times. Assuming Eq. (8) describes the relationship between hCi and y, we can approximately calculate Deff by ®tting diusion data to Eq. (8). Another approach is to derive Deff from the steady-state solution of the non-steady-state diusion equation at longer intervals of time. The ¯ux J for long diusion time is related to the concentration gradient by Fick's ®rst law, J ÿDeff rC:
12
The total ¯ux J can be obtained easily from the diusion coecient distribution D r and the concentration distribution C r through the microstructure at the steady state. rC in Eq. (12) is determined by the applied concentration gradient. Our numerical work shows Deff computed from Eqs. (8) and (12) has a dierence less than 2%, which can satisfy most of the practical applications. One of the main advantages of this approach is that it produces both the time-dependent
concentration pro®les in a microstructure and the eective properties of the microstructure. In addition, it avoids any complicated boundary conditions and thus allows one to study the diusion transport through any arbitrary microstructures. As test examples, we ®rst computed Deff for two of the simplest possible microstructures. The system is composed of parallel slabs of A and B with dierent diusion coecients DA and DB , shown in Fig. 5. In Fig. 5(a), the ¯ux is parallel to the plane of the slabs, and the arrangement is equivalent to a parallel electrical circuit. Each slab has the same concentration gradient, and most of the ¯ux is through the component with higher diusivity. Deff is given by Deff DA VA DB VB , where VA and VB are the volume fraction of each component. This corresponds to n 1 in Eq. (11). In contrast, for diusion ¯ux perpendicular to the plane of the slabs in Fig. 5(b), it is equivalent to a series electrical circuits. The ¯ux through each slab is equal, but the concentration gradients are dierent. In this case, n ÿ1 in Eq. (11) and Deff DA DB = DA VB DB VA . Fig. 6 shows the comparison of calculated and theoretical Deff as a function of volume fraction for these two parallel slab microstructures. Slab arrangements in Fig. 5(a) and (b) are indicated by jj and respectively. As can be seen from the graph, the eective diusion coecients computed agree well with those of theoretical solutions for both cases, thus validating our methods. Generally Deff for the two parallel slab microstructures will cover the complete range of
Fig. 5. Idealized two-phase microstructures ± parallel slabs: (a) ¯ux parallel to the plane of slabs jj; (b) ¯ux perpendicular to the plane of slabs ( ).
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
Fig. 6. Computed and theoretical eective diusion coecient Deff as a function of volume fraction of phase A for parallel slab microstructures, DA 1.0, DB 0.2.
eective properties between the upper and lower bounds as shown in the solid line and dashed line in Fig. 6. It would be expected that for arbitrary microstructures, Deff was within this range. For a polycrystalline material, diusion was carried out as described in the previous section. In particular, we examined the eect of the average grain size d of a grain structure on the eective diusion coecient. The results are plotted as the ®lled circles shown in Fig. 7, where Deff;simulation is calculated from steady state concentration pro®le. Deff;simulation varies inversely with the average grain
Fig. 7. Eective diusion coecient Deff as a function of average grain size d for a single-phase polycrystalline microP structure, D 1 ÿ 0:999 i g2 . Deff;upper and Deff;low are the eective diusion coecients predicted by mixture rules.
45
size d. This result is in good agreement with that of sharp-interface simulations by Swiler et al. [14]. The driving force for grain growth is the reduction of grain boundary energy. As the grain size becomes larger, the volume fraction of grain boundaries decreases, resulting in a decrease in the eective diusion coecient. In Fig. 7, we also plotted the theoretical eective diusion coecients Deff;upper and Deff;low as the upper and low bounds of the eective diusion coecients which can be predicted by rule of mixing as Dgb fgb Dv 1 ÿ fgb and Dgb Dv =Dgb 1 ÿ fgb Dv fgb , respectively, where fgb is the volume fraction of grain boundaries. However, for a diuse-interface description of the microstructure, there is some ambiguity in determining the exact values of Dgb , Dv and fgb . Recognizing that Deff predicted by rule of mixing is essentially a P volume average, we calculated Deff;upper as 1=N r D r where N is the total number of latticePpoints. Similarly Deff;low was computed as N = 1=D r. As expected, for all the grain sizes studied, Deff;simulation is between Deff;upper and Deff;low , which are the eective diusion coecients for idealized grain structures where grain boundaries can be treated as parallel slabs embedded in the grains. The dierence between Deff;simulation and Deff;upper was explained by Swiler et al. [15] as a result of triple junctions between grains as bottlenecks to diusion. Their sharp-interface simulations showed that Deff;simulation 0:27Dgb d=d [14], where d is the grain boundary width and d=d is the approximated volume fraction of grain boundaries. As a ®nal example, diusion through two-phase microstructures was studied. The two-phase microstructures in Fig. 8 were obtained from our phase ®eld simulations of phase separation of an initially homogeneous system. The two-phase system forms an inter-connected or a dispersed microstructure at dierent volume fractions. Fig. 9 shows the change of Deff (obtained from the steady state concentration pro®le) as a function of volume fraction for these microstructures (plotted with circle symbols). As the volume fraction of phase A with high diusion coecient DA DB increases, the shape of A droplets changes, and Deff increases gradually. After VA reaching about 50%, there is a sharp increase in Deff when the
46
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
Fig. 8. Two-phase microstructures for dierent volume fraction of phase A (shown in white color) showing the percolation transition, (a)±(e) for VA 16%, 32%, 50%, 62%, and 87%, respectively.
Fig. 9. Eective diusion coecient Deff as a function of volume fraction. Circles were calculated from the microstructures in Fig. 8. Dashed line and solid line were plotted from Eq. (11) representing the upper and lower bounds of Deff . DA =DB 1000.
microstructure changes from isolated phase A in B to inter-connected A and B phases. The ratio of DA =DB is three orders of magnitude in Fig. 9. The dashed line and solid line represent the upper and low bounds of Deff determined from Eq. (11) as n 1 and ÿ1, respectively. It is usually observed that the larger the distinction between the properties of dierent phases in heterogeneous systems, the more obvious is the percolation transition [18]. Our numerical experiments have demonstrated that the Fourier±Chebyshev spectral method is very ecient for solving the time-dependent diffusion equation in which the diusivity distribution is a smooth function throughout the
microstructure, especially for microstructures described by the diuse-interface model. The excellent agreement between our diuse-interface simulation and sharp-interface analytical or numerical methods showed that diuse-interface model is very effective for describing diusional transport in mesoscale microstructures. In recent years, the diuse-interface phase-®eld approach has been recognized as a powerful mathematical model for simulating the mesoscale morphological pattern formation such as grain growth and Ostwald ripening. By coupling the phase-®eld modeling of microstructure evolution, the approach developed in this paper allows us to model the eective diusivity evolution in an evolving microstructure. More details will be reported in the near future.
4. Conclusion We have simulated the temporal mass diusion transport through an arbitrary microstructure based on a phase-®eld approach. Diusion coecients were calculated as a function of the order parameters describing a single-phase grain structure where grain boundary diusion was much higher than bulk diusion and a two-phase microstructure where the diusivities varied with the phase. The semi-implicit Fourier±Chebyshev spectral method was shown to be accurate and ecient in solving the variable diusion equations. The accuracy of the model was demonstrated by reproducing the classic diusion results of grain boundary diusion in all three kinetic regimes
J. Zhu et al. / Computational Materials Science 20 (2001) 37±47
(type A, B and C) as well as the parallel slabs geometries. We demonstrated that one can extract eective diusion coecient of a microstructure from the temporal diusion pro®les. Furthermore, the sharp interface results of Swiler et al. were reproduced to show that the eective diusivity varies inversely with grain size. This model was used to predict the eective diusivity of twophase percolating and non-percolating structures where each phase has very dierent diusivities. This approach can be also used to study similar transport problems such as heat transport by recognizing the similarity between mass diusion and thermal diusion equations. Acknowledgements The authors are grateful for the ®nancial support from the Sandia National Laboratory and the National Science Foundation under Grant No. DMR 96-33719 and DMS 9721413. References [1] I. Kaur, Y. Mishin, W. Gust, Fundamentals of Grain and Interphase Boundary Diusion, third ed., Wiley, Chichester, 1995.
47
[2] S. Torquato, Int. J. Solids Struct. 35 (1998) 2385. [3] E.J. Garboczi, D.P. Bentz, Constr. Building Mater. 10 (1996) 293. [4] L.Q. Chen, Y.Z. Wang, JOM 48 (1996) 13. [5] D.N. Fan, L.Q. Chen, S.P. Chen, Mater. Sci. Eng. A 238 (1997) 78. [6] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid Dynamics, Springer, Berlin, 1987. [7] L.Q. Chen, J. Shen, Comput. Phys. Commun. 108 (1998) 147. [8] J. Shen, SIAM J. Sci. Comput. 16 (1995) 74. [9] J.C. Fisher, J. Appl. Phys. 22 (1951) 74. [10] R.T.P. Whipple, Philos. Mag. 45 (1954) 1225. [11] T. Suzuoka, J. Phys. Soc. Jpn. 19 (1964) 839. [12] H.S. Levine, C.J. MacCallum, J. Appl. Phys. 31 (1960) 595. [13] X. Gui, S.K. Dew, M.J. Brett, J. Appl. Phys. 74 (1993) 7173. [14] T.P. Swiler, V. Tikare, E.A. Holm, Mater. Sci. Eng. A 238 (1997) 85. [15] T.P. Swiler, E.A. Holm, Ceram. Trans. 69 (1997) 203. [16] L.G. Harrison, Trans. Faraday Soc. 57 (1961) 1191. [17] J. Philebert, Atom Movements ± Diusion and Mass Transport in Solids, Les Editions de Physique, Les Ulis, 1991. [18] C.W. Nan, Prog. Mater. Sci. 37 (1993) 1. [19] L.E. Nielsen, Predicting the Properties of Mixtures: Mixture Rules in Science and Engineering, Marcel Dekker, New York, 1978.