Available online at www.sciencedirect.com
Scripta Materialia 62 (2010) 85–88 www.elsevier.com/locate/scriptamat
Precipitate growth with composition-dependent diffusivity: Comparison between theory and phase field simulations R. Mukherjee,a,* T.A. Abinandanana and M.P. Gururajanb a
b
Department of Materials Engineering, Indian Institute of Science, Bangalore 560012, India Department of Metallurgical Engineering and Materials Science, Indian Institute of Technology-Bombay, Powai, Mumbai 400076, India Received 30 August 2009; revised 21 September 2009; accepted 22 September 2009 Available online 1 October 2009
We consider the growth of an isolated precipitate when the matrix diffusivity depends on the composition. We have simulated precipitate growth using the Cahn–Hilliard model, and find good agreement between our results and those from a sharp interface theory for systems with and without a dilatational misfit. With misfit, we report (and rationalize) an interesting difference between systems with a constant diffusivity and those with a variable diffusivity in the matrix. Ó 2009 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Phase field modeling; Phase transformation kinetics; Precipitation
In a recent paper [1] (hereafter referred to as Paper I), we used a phase field model to study growth of an isolated precipitate in one and two dimensions, and critically compared the simulation results with those from the classical sharp interface theory due to Zener [2] and Frank [3]. We also studied the effect of a dilatational misfit between the precipitate and matrix phases, and used our simulations as “computer experiments” to validate the sharp interface model of Laraia et al. [4]. In the phase field model in Paper I, the precipitate phase differs from the matrix in both composition (a conserved parameter) and an order parameter (a nonconserved parameter). We used this model so that its thermodynamic and kinetic features approximate closely the assumptions made in the theory of Zener and Frank, and that of Laraia et al.; specifically, it ensures a constant diffusivity in the matrix phase, and enabled us to focus our study on the effects due to a dilatational misfit and interface curvature. In this paper, we extend the work in Paper I to the classical Cahn–Hilliard model. Indeed, this model is simpler than the one used in Paper I in that it uses only one field variable, composition, to distinguish the precipitate and matrix phases. However, this model introduces a * Corresponding author. Tel.: +91 80 2293 2676; fax: +91 80 2360 0472; e-mail:
[email protected]
complication in that the diffusivity in the matrix phase is not a constant; therefore, even the results of the sharp interface theory need to be obtained numerically. Further, the extensive use of the Cahn–Hilliard model for studying precipitate coarsening (in addition to its deployment in studies of spinodal decomposition) demands that it be validated quantitatively against sharp interface results for both unstressed and stressed systems. In the present work, we consider a two-phase isothermal binary alloy system, where an isolated precipitate (p) is growing into a supersaturated matrix (m). The microstructure is described in terms of composition c(x, t) in a periodic domain X. In the absence of capillarity and elastic stress effects, the equilibrium compositions of precipitate and matrix phases are cpe ¼ 1 and cme ¼ 0, respectively. All the variables in this paper are in their scaled or non-dimensional form; we refer the reader to Ref. [5] for details of the non-dimensionalization procedure; see also the text below Eq. (6). The evolution of the composition field c(x, t) is governed by the Cahn–Hilliard equation [6]: @c ¼ r Mrl; @t
ð1Þ
where M is the atomic mobility and l is the chemical potential, defined as the variational derivative of the total free energy per atom, F , with respect to the local composition c,
1359-6462/$ - see front matter Ó 2009 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.scriptamat.2009.09.030
86
R. Mukherjee et al. / Scripta Materialia 62 (2010) 85–88
dðF =N V Þ : ð2Þ dc N V is the number of atoms per unit volume. For a stressed system, F is given by: Z Z 1 2 rij eel ð3Þ F ¼ N V ½f ðcÞ þ jðrcÞ dX þ ij dX; 2 X X l¼
where j is the gradient energy coefficient and f ðcÞ is the bulk free energy per atom, assumed to have the form: f ðcÞ ¼ c2 ð1 cÞ2 . eel ij is the elastic strain field and rij is the stress field that satisfies the equation of mechanical equilibrium: rij;j ¼ 0:
ð4Þ
Elastic stress and strain are related to each other through the generalized Hooke’s law: rij ¼ C ijkl eel kl ;
ð5Þ
where C ijkl is the elastic modulus tensor, assumed to be independent of composition and hence, the same everywhere. The elastic strain eel ij in Eq. (3) is given by: T eel ij ¼ eij bðcÞe dij ;
ð6Þ
where eij ¼ 0:5ðui;j þ uj;i Þ, in which ui is the displacement vector field that satisfies the equation of mechanical equilibrium, Eq. (4), and bðcÞeT dij is the (dilatational) composition-dependent eigenstrain, measured with respect to the strain-free matrix of composition c ¼ 0 as the reference state. bðcÞ ¼ c3 ð10 15c þ 6c2 Þ is a sigmoidal function that interpolates the matrix eigenstrain (which is assumed to be zero) to that in the precipitate phase. For the non-dimensionalization, characteristic length, energy and time were chosen as L0 0:4 nm, E0 1020 J and T 0 0:2 s, using interfacial energy = 0:08 J=m2 , free energy barrier height = 400 J=mol and interdiffusivity = 1:6 1018 m2 =s. With these choices, a non-dimensional shear modulus of 400 corresponds to about 250 GPa. The simulations of precipitate growth start with a particle of initial radius R ¼ 5 at the centre of a 1024 1024 grid with Dx ¼ Dy ¼ 1:0 (one-dimensional (1-D) simulations used 2 104 grid points with Dx ¼ 1:0); the composition inside the precipitate is set to cp ¼ 1:0 and that in the matrix is set to c1 ; we have studied precipitate growth in two alloys, with c1 ¼ 0:1 and 0:2. Microstructural evolution is studied by numerical integration of the Cahn–Hilliard equation, Eq. (1) (with M ¼ 1:0 and j ¼ 1:0). The integration involves the following two steps. (i) At a given time t, we first solve the equation of mechanical equilibrium, Eq. (4) (with shear modulus G=N V ¼ 400, Poisson’s ratio m ¼ 0:3 and precipitate–matrix misfit parameter eT ¼ 0:01), using a Fourier spectral technique [7]. (ii) Using the stress and strain fields, we evaluate the chemical potential field (see Eqs. (2) and (3)), and use it to integrate Eq. (1), using a semi-implicit Fourier spectral technique [8], to arrive at cðt þ DtÞ, on which steps (i) and (ii) are repeated. The Fourier transforms are carried out using the FFTW package developed by Frigo and Johnson [9].
For the Cahn–Hilliard model (with M ¼ 1) used in this study, the composition-dependent diffusivity in the matrix is given by: d 2f ¼ 2Mð6c2 6c þ 1Þ: ð7Þ dc2 Thus, DðcÞ decreases from Dðc ¼ 0Þ ¼ 2:0 in the saturated matrix phase to Dðc ¼ cs Þ ¼ 0 when the composition is on the chemical spinodal cs ¼ 0:5–120:5 0:2113; indeed, for c > cs , the diffusivity is negative. In all our simulations, therefore, the far-field composition is less than cs . During the simulations, we track the particle radius R (defined as the location at which c ¼ 0:5) with time, and compute the instantaneous growth coefficient a defined by R2t R2ts ¼ a2 s (see Eq. (9) below), where s ¼ 1000Dt. The instantaneous values of a, plotted against particle size R for a 2-D non-misfitting system with c1 ¼ 0:1 and 0:2 in Figure 1, show clearly the role of particle size on growth rate: the interface curvature raises the matrix interfacial composition through the well-known Gibbs–Thomson effect by DcmI and decreases the effective supersaturation for growth. The extent of this effect can be obtained from the expression DcmI ¼ vc=ðcpe cme ÞWm , where v is the mean interface curvature, c is the interfacial energy and Wm ¼ N V ½@ 2 f =@c2 cme . At R ¼ 10, for example, DcmI ¼ 0:0165; at R ¼ 50, we get DcmI ¼ 0:0033, which is a larger fraction (3:3%) of the supersaturation when c1 ¼ 0:1 than that (1:65%) when c1 ¼ 0:2. Thus, at R ¼ 50, the a values from the simulations are closer to that from sharp interface theory (horizontal lines) for the alloy with a greater c1 . We turn now to an outline of this sharp interface theory. The composition profile in the matrix around a growing precipitate obeys the standard diffusion equation: @c 1 @ k1 @c ¼ DðcÞr ; ð8Þ @t rk1 @r @r DðcÞ ¼ M
where r is the distance from the precipitate centre, k indicates the dimensionality, t is the time and DðcÞ is the composition-dependent diffusivity.
Figure 1. Growth coefficient a as a function of particle radius R for 2-D systems without misfit.
R. Mukherjee et al. / Scripta Materialia 62 (2010) 85–88
Zener [2] used dimensional arguments to show that the particle size R follows the parabolic growth law: R ¼ ak tð1=2Þ ;
ð9Þ
where a is the growth coefficient. Following Frank [3], we use the substitution w ¼ ðr=t1=2 Þ to transform Eq. (8) to the following form: k d dc w dc wk1 DðcÞ ¼ : ð10Þ dw dw 2 dw The first boundary condition is given by the interfacial equilibrium (neglecting capillarity effects): cðw ¼ aÞ ¼ cmi ;
ð11Þ
where cmi is the matrix interfacial composition; ignoring capillarity effects, it takes a value of cme ¼ 0 or cmE (see Eq. (13) below) for, respectively, the non-misfitting or the misfitting case. The second boundary condition is given by the growth rate: dc ðcp cme Þa a ¼ e ¼ ; ð12Þ dw w¼a 2D0 2D0 where D0 ¼ Dðcmi Þ is the matrix diffusivity at the interface. The third boundary condition is the far-field composition: cðw ¼ 1Þ ¼ c1 . While Eq. (10) can be solved analytically [2,3] for constant D, it has to be solved numerically when D is composition-dependent (Eq. (7)). Starting with the values of c (Eq. (11)) and dc=dw (Eq. (12)) at the interface (w ¼ a), Eq. (10) is integrated numerically to yield both c and ðdc=dwÞ at ðw þ nDwÞ, where Dw is the grid size used for discretizing w, and n is a positive integer. The integration procedure is repeated until dc=dw reaches zero, and c reaches a constant value of c1 . This procedure yields a (numerical) relation between the growth rate a and far-field composition c1 . In Figure 2 we have plotted the dependence of a on supersaturation c1 for 1-D, 2-D and 3-D systems without misfit. We have also presented in this figure a values obtained from our Cahn–Hilliard simulations in 1-D and 2-D systems. In the 1-D systems, there is no
Figure 2. Growth coefficient a as a function of supersaturation c1 in non-misfitting systems.
87
curvature effect, and the agreement between simulations and theory is excellent, being within 0:7% and 0:4% for alloy compositions c1 ¼ 0:1 and 0:2, respectively. In the 2-D systems, the data points represent the instantaneous value of a at a particle size of R ¼ 50; again, the agreement between simulations and theory is very good, within 1:8% and 0:5% for c1 ¼ 0:1 and 0:2, respectively (see Table 1). In other words, at the same size, curvature effects are smaller for more concentrated alloys. A more critical comparison between the phase field simulations and sharp interface theory is presented in Figure 3, which shows the scaled composition profiles for a non-misfitting system. It is clear that the match between the matrix composition profiles is very good, except for the region near the (compositionally diffuse) interface. We now turn to the effect of a dilatational misfit between the precipitate and matrix phases, presented in Figure 4. First, we note that the agreement between simulations and theory is good; it is within 6% and 0:6% for c1 ¼ 0:1 and 0:2, respectively (see Table 1). Moreover, the a vs. c1 curve for the misfitting case in Figure 4 is shifted to the right of the curve for the non-misfitting case. In other words, at the same supersaturation, the growth coefficient is lower in systems with misfit than in those without misfit. This may be explained by considering the role of misfit in elevating the matrix interfacial composition cmE : cmE ¼ cme þ DcmE ;
ð13Þ
Table 1. Comparison of 2-D growth coefficients from theory (at ) and simulations (as ). System
n
nE
at
as
No misfit
0.1 0.2
0.1 0.2
0.453 0.652
0.445 0.649
Misfit
0.1 0.2
0.072 0.172
0.337 0.532
0.318 0.529
Figure 3. A comparison of composition profiles from theory and simulations in two dimensions for a non-misfitting system with c1 ¼ 0:2. Distance is scaled by the particle radius R.
88
R. Mukherjee et al. / Scripta Materialia 62 (2010) 85–88
Figure 4. Growth coefficient a as a function of matrix supersaturation c1 in a misfitting system. The solid curve, for a non-misfitting system, is also shown for comparison.
where DcmE is the elevation in matrix interfacial composition due to misfit strains, given by [10] rmij ðemij epij Þ þ ½rpij ðepij eT dij Þ rmij emij =2 ; ðcpe cme ÞWm
ð14Þ
where superscripts m and p refer to the matrix and the precipitate phase, r and e are as depicted in Eq. (5) and (6). Thus, the effect of misfit is similar to the capillarity effect, in that it reduces the effective supersaturation from n ¼ c1 cme to nE ¼ c1 cmE (see Table 1), and therefore, for the same alloy composition, the growth rate is smaller for misfitting systems than that for non-misfitting systems. However, in contrast to the capillarity effect (which decreases with particle size), DcE is a constant, independent of particle size; it does depend on particle shape, but shape changes are not considered in our model. More importantly, the effect of misfit strain introduces an interesting difference between systems with constant diffusivity and those with a variable diffusivity. In the former (the case studied in Paper I), it can be shown that the rightward offset between the a vs c1 curves for the misfitting case from that for the non-misfitting case is constant at every value of a; it is not constant (see Fig. 4) when D depends on composition. The reason is as follows: in the Cahn–Hilliard model, D decreases with matrix composition. Since misfit strains elevate the composition at the interface, D0 ¼ DðcmE Þ is lower. Thus, for the same growth rate a,
the supersaturation has to be higher (with a lower farfield D) in the misfitting system than in the non-misfitting one. Thus, in Figure 4, the rightward offset between the two curves keeps increasing with increasing a. In our study, we have considered only one source of stress: particle–matrix misfit. In their sharp interface study [4], Laraia et al. also considered compositional stress in the matrix, arising from the composition dependence of matrix lattice parameter. While the misfit leads to a decrease in the effective supersaturation (Eq. (14)), the compositional strain may also change the diffusivity in the matrix. We may estimate the level of compositional strain in our study as follows: both misfit and compositional strains enter into the Cahn–Hilliard model through the eigenstrain (Eq. (6), whose composition dependence is given by the interpolation function, bðcÞ ¼ c3 ð10 15c þ 6c2 Þ. The compositional strain is only 0:00058 (or less than 6% of the precipitate–matrix misfit) even for c1 ¼ 0:2. As a result, the magnitude of compositional strain is too small to affect the growth kinetics in our simulations. In conclusion, we have studied precipitate growth into a matrix in which diffusivity is not constant using simulations based on the Cahn–Hilliard model and a sharp interface theory. Our results show that the precipitate growth behaviour in the simulations is in very good agreement with those from sharp interface theory in systems with and without misfit. The authors thank the Volkswagen Foundation for financial support for this work. They also thank Prof. V. Jayaram and Dr. H. Ramanarayan for discussions. [1] R. Mukherjee, T.A. Abinandanan, M.P. Gururajan, Acta Materialia 57 (2009) 3947–3954. [2] C. Zener, Journal of Applied Physics 20 (1949) 950–953. [3] F.C. Frank, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 201 (1950) 586–599. [4] V.J. Laraia, W.C. Johnson, P.W. Voorhees, Journal of Materials Research 3 (1988) 257–266. [5] M.P. Gururajan, T.A. Abinandanan, Acta Materialia 55 (2007) 5015–5026. [6] J.W. Cahn, Acta Metallurgica 9 (1961) 795–801. [7] A.G. Khachaturyan, Theory of Structural Transformations in Solids, John Wiley & Sons, New York, 1983. [8] L.Q. Chen, J. Shen, Computer Physics Communications 108 (1998) 147–158. [9] M. Frigo, S.G. Johnson, in: Proceedings of the 1998 IEEE International Conference on Acoustics, Speech and Signal Processing, vol. 3, 1998, pp. 1381–1384. [10] W.C. Johnson, Metallurgical Transactions A 18 (1987) 233–248.