- Email: [email protected]
Phase-field modeling of void formation and growth under irradiation
Available online at www.sciencedirect.com
Acta Materialia 60 (2012) 6112–6119 www.elsevier.com/locate/actamat
Phase-field modeling of void formation and growth under irradiation A.A. Semenov b, C.H. Woo a,⇑ a
Department of Physics and Materials, City University of Hong Kong, Kowloon, Hong Kong, China b Institute for Nuclear Research, Russian Academy of Sciences, Moscow, Russia Received 23 February 2012; received in revised form 21 May 2012; accepted 19 July 2012 Available online 23 August 2012
Abstract The application of phase-field modeling to void formation and growth under irradiation is analyzed. Here our focus is on the energy of the diffuse interfaces of the void system, which is customarily modeled by a Ginzburg-type gradient energy term with a parameterized coefficient. To correctly emulate the kinetics of void formation, the corresponding free-energy change is calibrated against the classical nucleation model, where the void–matrix boundary is well defined. We find that the gradient-energy coefficient may be conveniently treated as a constant independent of void size, even for voids down almost to the atomic size. We also find that the diffuse void interface in phase-field simulation has to be treated on an ultrafine spatial scale to reproduce the void energetic characteristics and growth behavior derived from the sharp boundary approach. Ó 2012 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Phase-field theory; Interfacial energy; Gradient energy coefficient; Void critical radius
1. Introduction Recent renewed interest in studying the effects of neutron-damage accumulation in aging reactors has brought to attention modeling techniques such as the phase-field model (PFM). The PFM has the advantage of allowing the kinetics of void nucleation and growth to be studied without moving boundaries as a spatially continuous process involving the entire vacancy concentration field. To properly model the evolution of a void ensemble in the framework of a phase-field approach, one must recognize that the formation of stable voids due to irradiation damage is physically a first-order transition process in the supersaturated solution of single vacancies. It is well known within the classical nucleation theory that void nucleation occurs when the free-energy barrier separating the metastable and stable states is overcome due to the action of stochastic fluctuations [1]. The associated statistics dictates that the nucleation rate of thermally stable voids typically depends exponentially on temperature and ⇑ Corresponding author.
E-mail address: [email protected] (C.H. Woo).
critical size of void embryo [2–5]. Void nucleation occurring this way is stochastic and local, and is fundamentally very different from that due to spinodal decomposition, for which the phase-field methodology is initially developed [6–11]. Indeed, spinodal decomposition is a deterministic and global process that occurs when the homogeneous solid solution becomes unstable, not metastable. There is no energy barrier to overcome. The absence of an activation barrier corresponds to a zero critical size for the new phase embryo. This means that, unlike the case of metastable solution, in which the detailed behavior of the stochastic fluctuations plays a fundamental and quantitative role, in the spinodal decomposition, stochastic fluctuations only serve as a noise in the vacancy concentration field to trigger the initiation of phase separation process. As a result, the phase separation takes place globally throughout the initially homogeneous solution [10]. In this regard, the adoption of the phase-field approach in modeling irradiation damage processes must be conducted with due care. In a first-order process, the probability of nucleation depends exponentially on the height of the activation barrier, which contains contributions from the energy of the
1359-6454/$36.00 Ó 2012 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.actamat.2012.07.049
A.A. Semenov, C.H. Woo / Acta Materialia 60 (2012) 6112–6119
interface between the critical nucleus and the metallic matrix [2–5]. Even after nucleation, the survival probability of a supercritical void still depends exponentially on the interfacial energy through the rate of vacancy emission [12]. It is thus obvious that an accurate account of this energy contribution is of paramount importance in the usefulness of the phase-field method to model irradiation damage effects. Within the phase-field methodology, the interfacial energy in the free-energy functional is customarily described using a Ginzburg-type gradient energy term with a constant coefficient [6–11]. However, it was shown [1] that, because of the diffusive nature of the void–matrix boundary, the “interfacial” energy must also contain a contribution from the bulk terms in the free energy functional in order to correctly emulate the void thermodynamics. In practical computations, the bulk contribution does depend on the void-size [1], which itself is not a well defined quantity in the phase-field model. This treatment then runs into a self-consistency problem, particularly for small voids. The aim of the present work is to find a practical solution for this issue. 2. Interfacial energy In the phase-field approach, vacancies are assumed to be continuously distributed everywhere in space, including voids, which are treated as specific regions of the vacancy concentration field cv(x), where the local concentration of single vacancies approaches unity, i.e., cv(x) ffi 1. The behavior of the concentration field is assumed to be governed by the thermodynamic potential, conventionally taken to be the phenomenological Cahn–Hilliard-type functional /(cv(x)) [6–11]: Z kBT 2 ½/ðcv ðxÞÞ þ j2 ðrcv Þ dx ð1Þ Uðcv ðxÞÞ ¼ X V where X is the atomic volume, T the absolute temperature, kB the Boltzmann constant. The integral is conducted over the total volume V. The Ginzburg-type gradient energy term with the empirical j-coefficient accounts for the free energy associated with the diffuse void–matrix interface, and is a key ingredient in the phase field model. The bulk free energy density /(cv) in Eq. (1) can be expressed without loss of generality as the dilute ideal solution limit cv lnðcv =ecv1 Þ [13] minus a deviation denoted by u(cv): /ðcv Þ ¼ cv lnðcv =ecv1 Þ uðcv Þ
ð2Þ
Here cv1 is the equilibrium concentration of single vacancies. In the dilute solution limit, i.e., for cv ! 0, u(cv) must satisfy juðcv Þ=cv lnðcv =ecv1 Þj ! 0. In practice, u is usually represented as a phenomenological function with parameters fitted to experimental data. To ensure the existence of a stable “void” phase, /(cv) must also have a global minimum at cv ffi 1.
6113
We have previously shown that in an initially uniform dilute vacancy solution of concentration c0v , i.e., uðc0v Þ ffi 0, the change in the thermodynamic potential DU associated with the formation of a “void” occupying a volume V0 is given by [1] Z kBT DU ¼ ½cv ðxÞ lnðc0v =cv1 Þ þ ð/ðcv ðxÞÞ /ð1ÞÞdx X V0 Z 2 þj2 ðrcv Þ dx ð3Þ V
where /(1) is set equal to (c0v ) to ensure a global minimum of /(cv) at cv = 1. Eq. (3) has been derived assuming that the deviation dcv(x) = cv(x)c0v is small, i.e., |dcv(x)| c0v , everywhere except inside the void volume V0 V, and the law of mass conservation is obeyed, i.e., Z Z dcv ðxÞdx þ dcv ðxÞdx ¼ 0 ð4Þ V V 0
V0
In a supersaturated solution (c0v > cv1) the density of small subcritical voids, serving as embryos for the supercritical voids, is proportional to exp(DU/kBT) [2–5]. In addition, the supercritical voids, which are growing directly from the vacancy solution without the assistance of stochastic fluctuations, nucleate with a rate that is also proportional to exp(DU/kBT), DU being the free energy of formation of a critical void [5]. Thus, to ensure the correct emulation of the kinetics of void precipitation in the phasefield framework, the j-coefficient has to be defined such that DU in Eq. (3) is equal to the corresponding free energy change in the sharp boundary approach. In the latter case, the free energy change associated with the formation of a void with volume V0 is given by [13]: DU ¼ ðV 0 =XÞk B T lnðc0v =cv1 Þ þ rR
ð5Þ
where r is the surface tension (i.e., surface energy density), and R is the void surface area. The volume of a void within the phase-field framework, however, cannot be defined uniquely because of a diffuse character of the interface. This means that DU in Eqs. (3) and (5) do not necessarily correspond to the equivalent voids. Indeed, noting that kBT ln(c0v =cv1 ) is the chemical potential of a vacancy in the dilute ideal solution, the first term of Eq. (5) gives the bulk free energy reduction when a void is formed with m = V0/X vacancies. Similarly, the first term of Eq. (3) may also be interpreted as the free energy change caused R by the formation of a void with volume 0 0 V 0 ¼ m X ¼ V 0 cv ðxÞdx. The spatial layer between the volumes V0 and V 00 may then be treated as part of the diffusive interface between the “void” formed by m0 vacancies and the matrix. Then, using V 00 and R0 (where R0 is the surface area of the volume V 00 ) in Eq. (5) instead of V0 and R, the interfacial energy may be written as Z Z kBT 2 0 2 rR ¼ ð/ðcv ðxÞÞ /ð1ÞÞdx þ j ðrcv Þ dx X V0 V ð6Þ
6114
A.A. Semenov, C.H. Woo / Acta Materialia 60 (2012) 6112–6119
Alternatively, equating Eqs. (3) and (5), we arrive at Z kBT rR ¼ ½ð1 cv ðxÞÞ lnðc0v =cv1 Þ þ ð/ðcv ðxÞÞ /ð1ÞÞ X V0 Z 2 2 dx þ j ðrcv Þ dx ð7Þ V
The diffuse character of the interface also allows Eq. (3) to be considered as the free energy change associated with the formation of a void of size intermediate between m0 and m. Similar to Eqs. (6) and (7), the corresponding expression for the interfacial energy can be easily written down as well. Note that, similar to the gradient term, contributions to the integrals over V0 in Eqs. (6) and (7) mainly come from the interfacial region, where c0v < cv(x) 6 1. These contributions may therefore be considered as part of the interfacial energy. The critical size for large voids, which corresponds to a maximum of DU in Eq. (5), can only be recovered within the phase-field framework when both the gradient term and the bulk term in Eqs. (6) and (7) are taken into account [1]. In such cases, the difference between Eqs. (6) and (7) has been shown to be insignificant [1]. The case of small voids is considered in the next section. As mentioned in the foregoing, equalities like Eqs. (6) and (7) can be used for evaluating the gradient energy coefficient j, which is a model parameter in the phase-field approach. Obviously, values of this coefficient are generally void-size-dependent and, moreover, are not necessarily equal to each other. Nevertheless, there is an alternative treatment, which allows the corresponding dependence to be determined in a unique way. This is the subject of the following section. 3. Gradient energy coefficient In this section, we establish the condition of equivalence between voids in the two approaches and find out the phase-field vacancy distribution cv(x) that corresponds to a void of radius R in the sharp boundary case. Let us consider a vacancy solution in equilibrium with the said void. In this regard, the void is a critical one with a radius that maximizes DU in Eq. (5) [13]. Thus, the solution concentration c0v equals cv1exp(2rX/kBTR), and DU in Eq. (5) can be rewritten in a form that does not explicitly depend on the void radius R, i.e., DU ¼
4pR2 r 16pr3 X2 ¼ 2 3 3ðk B T Þ ln2 ðc0v =cv1 Þ
ð8Þ
Within the framework of the phase-field model, we focus our attention on the vacancy concentration field cv(x), which is in equilibrium with vacancy solution of the same concentration c0v . In equilibrium, the chemical potential l(cv(x)) of the vacancy solution is constant throughout the entire volume, i.e., lðcv ðxÞÞ ¼ k B T lnðc0v =cv1 Þ. Expressing the chemical potential as the functional derivative of the thermodynamic potential in Eq. (1), we may write down the following equation for cv(x):
d/ðcv Þ 2j2 Dcv ðxÞ ¼ lnðc0v =cv1 Þ dcv subjected to the following boundary conditions: ðrcv ðxÞÞx¼0 ¼ 0;
cv ðxÞjjxj!1 ! c0v
ð9Þ
ð10Þ
We note that the void–matrix boundary is not a mathematical boundary in the phase-field model. Assuming spherical symmetry, Eq. (9) reduces to ( ) 2 2 d 4j2 dcv 2 dcv j þ U ðcv ðrÞÞ ¼ ð11Þ dr dr r dr where U ðcv Þ ¼ cv lnðc0v =cv1 Þ /ðcv Þ ¼ cv lnðcv =ec0v Þ þ uðcv Þ 2
ð12Þ
2
Expressing j (dcv/dr) from Eq. (11) and taking into account the boundary conditions (10), DU from Eq. (3) becomes Z 1 8pk B T DU ¼ ðU ðcv ðrÞÞ c0v Þr2 dr ð13Þ X 0 Using the boundary condition (10) at infinity and the definition of the volume V0, outside of which |dcv(r)| is considered to be much smaller than c0v , we may calculate the integrals in Eq. (3) in the limit of r ! 1. The explicit void-volume dependence of Eq. (13) disappears in this case. In a similar way, it can be shown that Eq. (13) can also be derived by substituting the solution cv(r) of the boundary value problem (10)–(12) directly into Eq. (1). Equating DU in Eqs. (8) and (13), we obtain an explicit expression for the j-coefficient: !1=3 !1=3 2p ln c0v =cv1 2rX 2p j¼ ¼R k B T 3Iðc0v ;cv1 Þln2 ðc0v =cv1 Þ 3I c0v ; cv1 ð14Þ where Iðc0v ; cv1 Þ
¼2
Z
1
ðU ðcv ðqÞÞ c0v Þ4pq2 dq
ð15Þ
0
The last equality in Eq. (14) is derived from the expression for the critical void radius: R = 2rX/[kBTln(c0v =cv1 )]. Making the transformation r = jq in Eqs. (10)–(12), it can be seen that the solution cv(q) of the boundary value problem has no explicit dependence on j, and consequently, neither does the integral I(c0v , cv1). The gradient energy coefficient j can then be calculated directly using Eq. (14). Indeed, in terms of the variable q, Eq. (11) can be transformed into the following equation by integrating once: " # " # 2 2 dcv dcv ¼ ½U ðcv ðqiþ1 ÞÞ dq dq q¼qiþ1
q¼qi
U ðcv ðqi ÞÞ
Z
qiþ1 qi
2 4 dcv dq0 q0 dq0
ð16Þ
A.A. Semenov, C.H. Woo / Acta Materialia 60 (2012) 6112–6119 10 0
c v∞ = 10
0
5 0
0
0 -10
-φ
ð17Þ
ð19Þ
To proceed further, the explicit form of the free energy density /(cv) has to be specified. Following Ref. [1], we write down this function as /ðcv Þ ¼ cv lnðcv =ecv1 Þ þ ½n lnðcv1 Þ þ ðn þ 1Þð1 c0v Þcnv ½ðn 1Þ lnðcv1 Þ þ nð1 c0v Þcnþ1 v
-9
Cv = 10
0
For the asymptotic cv(q) to satisfy Eq. (16), it must fulfill the condition 2
-8
Cv = 10
Cv = 10
Here q0 is a point such that 0 < dcv ðq0 Þ=c0v 1. Differentiating Eq. (17), we obtain qffiffiffiffiffiffiffi i dcv ðq Þ h ðdcv =dqÞq¼q0 ¼ pffiffiffiffiffiffi0ffi 1 þ 2c0v =q0 2c0v rffiffiffiffi rffiffiffiffi c0v dcv ðq0 Þ c0v ffi ð18Þ << 0 2 cv 2
ðdcv =dqÞq¼q0 >> juðc0v Þj
-7
Cv = 10
dcv ðqÞ ¼ cv ðqÞ c0v
qffiffiffiffiffiffiffi dcv ðq0 Þq0 expððq q0 Þ= 2c0v Þ > 0 ffi q
-6
Cv = 10
−10
U
The latter provides a recursive relation via which dcv/dq and hence cv(q) can be readily obtained as a function of cv(0) using a finite difference scheme, with dcv/dq = 0 at q = 0. The value of cv(0) can then be varied to obtain the solution that satisfies the asymptotic behavior of cv(q) at q ! 1. In the immediate vicinity of c0v (i.e., for jcv c0v j c0v ), where U(cv) maximizes, U(cv) can be expanded in the form U ðcv Þ ffi c0v ðdcv Þ2 =ð2c0v Þ. The solution of Eqs. (10) and (11) then gives the following asymptotic description of cv(q):
6115
ð20Þ
where n P 2, which provides the correct dilute ideal solution limit. For inequality (19) to be valid even for the vacancy concentrations near the melting temperatures, i.e., for c0v 105–104, the exponent n in Eq. (20) should satisfy n P 2.5. 4. Results and discussions Profiles of the “potential” U(cv) calculated for various values of c0v are shown in Fig. 1. For larger values of n, the curvature of U(cv) near the maximum at cv = 1 increases, its minimum deepens [1], and the diffuse interface becomes thinner, similar to the case of the planar surface considered earlier [1]. Thus, to minimize the number of grid points for the numerical calculation, we take the smallest possible value for n, which is 2.5. For n > 2, the behavior of the free energy density /(cv) starts to deviate significantly from the ideal solution (/ cvln(cv/ecv1)) when cv > 102 [1]. For weaker solutions, i.e., when c0v cv(q) < 102, we have shown in the Appendix that the following relation holds for the spatial behavior of cv(q):
-5 0.0
0.1
0.2 0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
cv Fig. 1. “Potential” U(cv) at cv1 = 1010 as a function of vacancy concentration for the exponent n = 2.5 and various values of vacancy supersaturation c0v . Arrow indicates the drop in the “potential energy” when the “motion” takes place from cv ffi 1 at r = 0 to cv = c0v at r ! 1.
Dq ¼ q0 q sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffi cv ðqÞ 1 0 ð1 þ ln 2 ðc ðqÞ=c ÞÞ 2c0v v v lnðcv ðqÞ=c0v Þ lnðdcv ðq0 Þ=c0v Þ
ð21Þ
which describes a very sharp, albeit diffuse, interface. Indeed, within a spatial range of Dr = jDq < 101j, where j has a value of at most a few lattice constants [1], the vacancy concentration in the interface region could drop by eight orders of magnitude from a value of 102. Accumulation of vacancies in the void is driven by the difference between the actual concentration of vacancies and the concentration of the equilibrium vacancy solution c0v [4]. Under reactor operating conditions, c0v is often smaller than 1010 [1]. If we assume, for the sake of argument, that a feasible upper bound for the actual vacancy concentration could be as large as the melting-temperature equilibrium concentration of 105–104, then, from Eq. (21), the vacancy concentration in the interface has to drop five to six orders of magnitude within a very narrow spatial range Dq 103–102. Thus, to correctly represent in the phase-field approach the concentration gradient encountered in irradiation-induced void growth, numerical calculations would have to be performed on an extremely fine sub-nanometer spatial scale, a formidable task for the computation. Eq. (16) is numerically solved using a finite difference scheme, with the left-hand-side boundary value cv(0) set near 1, and the derivative (dcv/dq)q=0 = 0. The solution is repeated, varying values of cv(0), until the asymptotic boundary condition (18) is fulfilled. To ensure that the calculated vacancy concentration field cv(q) is the solution of
6116
A.A. Semenov, C.H. Woo / Acta Materialia 60 (2012) 6112–6119 0
0
10
10
-1
10
-1
10
C v∞ = 10
-2
C v∞ = 10
-2
−12
10
−10
10
-3
10
cv
-6
Cv= 10
0
-7
Cv= 10
0
-8
Cv= 10
0
-9
Cv = 10
cv
0
-3
10
0
0
-9
-8
Cv= 10
1.0
1.5
Cv= 10
-4
10
0
0
-10
-11
Cv = 10
Cv= 10
-5
10
-4
10
-6
10 -5
10
-7
10
-8
-6
10
0.0
0.5
1.0
1.5
-2
2.0
2.5
10
3.0
0.0
0.5
-3
10
10
10
-4
10
-4
10
3.0
-7
10
-3
2.5
-6
-3
10
10
2.0
-5
10 -4
0
10
-6
-6
10
Cv = 10
-8
-5
0
10
-9
Cv = 10
0
-8
Cv = 10
0
-11
Cv = 10
-9
-6
10
10
-7
10
10
-5
10
-10
-7
10
10
-8
10
-8
-6
1.00
1.02
1.04
1.06
1.08
10 1.10
1.05
2.908 2.912 2.916 2.920 2.924
Eq. (16) with boundary conditions (10), the consistency of the free energy changes D/ calculated using both Eqs. (3) and (13) is verified. In Figs. 2 and 3, we present the calculated concentration profiles for two different equilibrium concentrations cv1 and various levels of vacancy supersaturation c0v =cv1 . For copper and molybdenum, parameters are chosen to give values of critical radii between 0.3 and 1.5 nm, which are typical for these materials under irradiation for a wide range of temperatures [12]. The spatial positions of the sharp void–matrix boundary can be calculated from Eq. (14). They are found to be within the diffuse interface, which is consistent with our general treatment given in Section 2. Accordingly, the sharp boundary position R00 can also be estimated directly from the phase-field calculations as 1=3 3 0 n cv ; cv1 R00 =j ¼ ð22Þ 4p
nðc0v ; cv1 Þ ¼ 4p
Z
1
cv ðqÞq2 dq
ð23Þ
0
is the number of vacancies in the “void” in the phase-field approach. The ratio of R00 =R is presented in Table 1, from
1.06
1.07
1.08
10 3.1400 3.1402 3.1404 3.1406 3.1408 3.1410
Fig. 3. Same as Fig. 2, but at cv1 = 1012.
Fig. 2. Concentration profiles cv(r/j) at cv1 = 1010 for various values of vacancy supersaturation c0v . Arrows show spatial positions of the sharp boundary calculated with Eq. (14). Two inserts demonstrate the behavior of concentration profiles on finer spatial scales.
where
-11
10
-9
10
Table 1 Ratio R00 =R. cv1 10
10 1012
c0v =cv1 101
102
103
104
0.992 0.993
0.999 0.995
1.022 1.007
1.066 1.031
which one may conclude that Eq. (22) in the phase-field model does provide a good approximation for the sharp boundary position. From the inserts in both figures the validity of the approximation (21) can be also verified. In Table 2, values of the j-coefficient calculated with Eq. (14) are listed. The value of j is relatively constant, almost independent of the void size R. Depending on temperature, it is maintained at about one or two lattice constants a. This is consistent with our estimation in the foregoing section, from which we have concluded that the diffuse void interface is very sharp and has to be realized on a very fine spatial scale. For such thin interfaces the integrals in Eq. (3) can be separated with sufficient accuracy into bulk and surface contributions. This explains why values of R and R00 are close to each other. It also explains why the j-value only weakly depends on the void size, and why the j-value obtained in the plane boundary case (c0v =cv1 ¼ 1) [1] can be used in studies even for voids with
A.A. Semenov, C.H. Woo / Acta Materialia 60 (2012) 6112–6119 Table 2 Gradient term coefficient j/(2rX/kBT). cv1 10
10 1012
c0v =cv1 1.0
101
102
103
104
0.173 0.156
0.169 0.154
0.168 0.152
0.169 0.151
0.173 0.153
radii as small as one lattice constant. Indeed, when c0v =cv1 ¼ 106 , which corresponds to voids formed with just a few vacancies in the case of molybdenum and copper, the width of the interface becomes comparable with the void radius. Computer calculations show that only then the jcoefficient will start to deviate significantly from the plane boundary value. It is worth noting that in the phase-field theory of solidification, the Cahn–Hilliard analysis also often relies on the constancy of the j-coefficient [14–20]. Physically, the gradient term in the total free energy expression represents spatial correlation effects, and the j-coefficient in this term is related to the second derivatives of the direct correlation function of the liquid evaluated at the reciprocal lattice vectors of the corresponding crystal [14,21,22]. Although the j-coefficient is constant when the interpolated bulk free energy density / is assumed to be in a polynomial form (piecewise parabolic [15–17], triple-parabolic [18], quartic [19,20]), yet the corresponding effective solid–liquid interfacial tension r depends on the size of the critical nucleus [17–20]. This dependence is especially strong at high undercoolings, when the critical nucleus has a size comparable to the thickness of the interface. The sharp-interface assumption of the classical nucleation theory is invalid in such cases. As a result, in agreement with the experimental observations, the nucleation rates of crystalline structures in highly undercooled liquids calculated within the diffuse interface approach is orders of magnitude higher than those from the classical nucleation theory [17,19,20]. It also reflects the fact that the solid– liquid interface has a finite width and is diffuse on the molecular scale, extending over several molecular layers [11]. On the other hand, the interface considered in the present work remains thin compared to the size of the nucleus, even for vacancy supersaturations as high as c0v =cv1 ¼ 104 . As a result, the constancy of the surface tension in the sharp interface approach translates into the constancy of the gradient energy coefficient, and vice versa. Similar to Refs. [14,15], the diffuse nature of the interface in the present case does not significantly affect the interfacial tension because the spinodal point where vacancy solution becomes unstable (as opposed to being metastable) remains sufficiently far away. Indeed, the second derivative d2// (dcv)2 becomes negative above the spinodal point [1]. At the same time, even under in-reactor conditions with continuous production of vacancies, the vacancy solution is still dilute and ideal. Thus, the concentrations c0v of vacancy solutions considered here are far below those (c0v > 102 [1]) for which polynomial-type deviation u(cv) from the free
6117
energy density cv lnðcv =ecv1 Þ of the dilute solution is important. We will show in a forthcoming paper that the dynamical behavior of an individual void, well-known from the ratetheory treatment of void evolution [23], can be quite accurately reproduced in the framework of the phase-field model considered here. However, the characteristic spatial scales of the void–metal diffuse interface found in the present analysis presents a challenge to modeling the evolution of a void ensemble under irradiation. As a possible solution to this problem, one may consider smearing out the interface and then using the method of matched asymptotic expansion [11,24] to relate the phase-field model to the sharp interface approach. Since the vacancy concentration in the solution is very low, the void growth is very slow compared to the relaxation time of the vacancy diffusion field, so that the diffusion field in front of the interface can be considered at steady state. Thus, in the case of voids the interface Pe´clet number [11,24] is vanishingly small, which will significantly simplify the corresponding matched asymptotic analysis. 5. Summary and conclusions In this paper, we consider phase-field modeling of void development as a first-order phase transition in a metastable vacancy solution. The thermodynamic barrier against the nucleation of thermally stable growing voids is caused by the formation of the void–metal interface. Under the conventional framework, in which the Ginzburg-type gradient energy term is used to represent the interfacial energy, the associated coefficient is a model parameter to ensure the equivalence of the physical and the emulated void. Establishing equivalence between voids in the phase-field model with the diffuse interface and physical voids with the sharp boundary allows the diffuse model to elucidate the concentration of vacancies and structure of voids in a thermodynamic context, something that molecular dynamics cannot do. This forms the focus of the present study. In general, it is expected that the solution of the voidgrowth kinetic equation would be complicated by the void-size dependence of the gradient-energy coefficient. However, we found that for a fixed temperature this parameter may be treated, to a good approximation, as a constant, not only for large voids, but even for voids with sizes comparable with a lattice constant. This result can be traced to the sharpness of the diffuse interface, which has a characteristic width of less than one lattice constant, despite its diffuse character. Thus, the actual spatial position of the sharp boundary can be interpolated very well by its diffuse analogue, even in the case of very small voids. The sharpness of the diffuse interface is a direct result of the behavior of the thermodynamic potential in the dilute ideal solution limit, and the requirement that the Cahn–Hilliard-type interpolative functional must satisfy this limit. This requirement is absolutely necessary because under real experimental conditions, even under
6118
A.A. Semenov, C.H. Woo / Acta Materialia 60 (2012) 6112–6119
continuous production of vacancies by irradiation, the supersaturated vacancy solution in which void nucleation and growth take place is still dilute and ideal. Acknowledgement This project is initiated and funded by Grant 534409 from the Hong Kong Research Grant Commission, to which the authors are thankful.
and using the relation M(a, a; z) = exp(z), Eq. (A5) can be expressed through the tabulated function M(1/2, 1/2; z): sffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi" cv ðqÞ c0v q0 q ffi 2 1 M 0 cv ðqÞ lnðcv ðqÞ=cv Þ qffiffiffiffiffiffiffi 1 cv ðqÞ dcv ðq0 Þ 0 ln 2c 1=2; 1=2; ln v 2 c0v c0v
ðA7Þ
Using the asymptotical behavior CðcÞ z ac 1 e z ½1 þ Oðjzj Þ CðaÞ
Appendix A
Mða; c; zÞ ¼
In the one-dimensional case Eq. (9) can be written as 2 dcv þ U ðcv ðqÞÞ ¼ const ðA1Þ dq
where C(a) is the gamma function, we arrive at Eq. (17). In the case of spherical symmetry Eq. (9) has the form
The boundary condition at q ! 1 gives that const ¼ c0v . Thus, the solution can be presented in the following integral form: Z cv1 dcv pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q0 q1 ¼ 0 cv U ðcv Þ cv0 Z cv1 1=2 ffi ðc0v þ cv lnðcv =ec0v ÞÞ dcv ðA2Þ cv0
Here cv0 = cv(q0), and cv1 = cv(q1). The second equality in Eq. (A2) is valid when cv 6 102. If both d1 = dcv1/c0v and d0 = dcv0/c0v 1, Eq. (A2) reduces to qffiffiffiffiffiffiffi q0 q1 ¼ 2c0v lnðd1 =d0 Þ ðA3Þ On the other hand, for dcv1/c0v 1 and c0v cv(q) 6 102, a major contribution to the integral (A2) comes from the interval of integration, where the difference (c0v U(cv)) can be approximated by cvln(cv/c0v ). Then, using the substitution y = 0.5ln(cv/c0v ), we have that qffiffiffiffiffiffiffi Z 0:5 lnðcv ðqÞ=c0v Þ y e dy q1 q ffi 2c0v pffiffiffi y 0:5 lnð1þd1 Þ qffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi lnðcv ðqÞ=c0v ÞMð1=2; 3=2; 0:5 lnðcv ðqÞ=c0v ÞÞ ¼ 2 c0v p lnð1 þ d1 ÞMð1=2; 3=2; 0:5 lnð1 þ d1 ÞÞ ðA4Þ where M(a, c; z) is the confluent (degenerate) hypergeometric function [25,26]. Noting that the second term in Eq. (A4) is much less than the first one, and that the correct expression for the sum of Eqs. (A3) and (A4) has to be independent of d1, we may write down the interpolation for both equations as follows: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q0 q ¼ 2 c0v lnðcv ðqÞ=c0v ÞMð1=2; 3=2; 0:5 qffiffiffiffiffiffiffi lnðcv ðqÞ=c0v ÞÞ 2c0v lnðdcv ðq0 Þ=c0v Þ ðA5Þ From the recurrence relation [15], cMða; c; zÞ cMða 1; c; zÞ ¼ zMða; c þ 1; zÞ
ðA6Þ
2
d 2 cv 4 dcv dU ¼ þ dcv dq2 q dq
ðA8Þ
ðA9Þ
It can be approximated by the foregoing one-dimensional equation when dU 4 dcv ðA10Þ q dq dcv According to Eq. (A4), the derivative dcv/dq is approximately given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 dcv 4 cv ðqÞ lnðcv ðqÞ=c0v Þ ffi q dq q 2Dq lnðcv ðqÞ=c0v Þ ffi ðA11Þ q Since dU/dcv = ln(cv/c0v ) for c0v < cv 6 102, the inequality (A10) is fulfilled when the void radius is much larger than the characteristic width of the interface. References [1] Semenov AA, Woo CH. J Nucl Mater 2011;411:144. [2] Zeldovich YaB. Sov J Exp Theor Phys 1942;12:525. [3] Frenkel YaI. Kinetic theory of liquids. Moscow: Academy of Science of the USSR; 1959. [4] Lifshitz EM, Pitaevskii LP. Physical kinetics. Oxford: Pergamon Press; 1981 [Chapter XII]. [5] Slezov VV, Schmelzer J, Tkatch YaY. J Chem Phys 1996;105:8340. [6] Cahn JW, Hilliard JE. J Chem Phys 1958;28:258. [7] Hillert M. Acta Metall 1961;9:525. [8] Cahn JW. Acta Metall 1961;9:795. [9] Cahn JW. Acta Metall 1962;10:179. [10] Skripov VP, Skripov AV. Uspekhi Fiz Nauk 1979;128:193 [in Russian]. [11] Emmerich H. Adv Phys 2008;57:1. [12] Semenov AA, Woo CH. Phys Rev B 2002;66:024118. [13] Landau LD, Lifshitz EM. Statistical physics. Oxford: Pergamon Press; 1980. [14] Harrowell P, Oxtoby DW. J Chem Phys 1984;80:1639. [15] Bagdassarian CK, Oxtoby DW. J Chem Phys 1994;100:2139. [16] Wild R, Harrowell P. Phys Rev E 1997;56:3265. [17] Gra´na´sy L, Iglo´i F. J Chem Phys 1997;107:3634. [18] Gra´na´sy L, Oxtoby DW. J Chem Phys 2000;112:2399. [19] Gra´na´sy L, Bo¨rzo¨nyi T, Pusztai T. Phys Rev Lett 2002;88:206105. [20] Gra´na´sy L, Pusztai T, James PF. J Chem Phys 2002;117:6157. [21] Haymet ADJ, Oxtoby DW. J Chem Phys 1981;74:2559.
A.A. Semenov, C.H. Woo / Acta Materialia 60 (2012) 6112–6119 [22] Haymet ADJ, Oxtoby DW. J Chem Phys 1982;76:6262. [23] Brailsford AD, Bullough R. J Nucl Mater 1972;44:121. [24] Edler KR, Grant M, Provatas N, Kosterlitz JM. Phys Rev E 2001;64:021604.
6119
[25] Gradshteyn IS, Ryzhik IM. Table of Integrals, series and products. New York: Academic Press; 1965. p. 1058. [26] Abramowitz M, Stegan IA, editors. Handbook of mathematical functions. New York: Dover; 1965. p. 504.
Acta Materialia 60 (2012) 6112–6119 www.elsevier.com/locate/actamat
Phase-field modeling of void formation and growth under irradiation A.A. Semenov b, C.H. Woo a,⇑ a
Department of Physics and Materials, City University of Hong Kong, Kowloon, Hong Kong, China b Institute for Nuclear Research, Russian Academy of Sciences, Moscow, Russia Received 23 February 2012; received in revised form 21 May 2012; accepted 19 July 2012 Available online 23 August 2012
Abstract The application of phase-field modeling to void formation and growth under irradiation is analyzed. Here our focus is on the energy of the diffuse interfaces of the void system, which is customarily modeled by a Ginzburg-type gradient energy term with a parameterized coefficient. To correctly emulate the kinetics of void formation, the corresponding free-energy change is calibrated against the classical nucleation model, where the void–matrix boundary is well defined. We find that the gradient-energy coefficient may be conveniently treated as a constant independent of void size, even for voids down almost to the atomic size. We also find that the diffuse void interface in phase-field simulation has to be treated on an ultrafine spatial scale to reproduce the void energetic characteristics and growth behavior derived from the sharp boundary approach. Ó 2012 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Phase-field theory; Interfacial energy; Gradient energy coefficient; Void critical radius
1. Introduction Recent renewed interest in studying the effects of neutron-damage accumulation in aging reactors has brought to attention modeling techniques such as the phase-field model (PFM). The PFM has the advantage of allowing the kinetics of void nucleation and growth to be studied without moving boundaries as a spatially continuous process involving the entire vacancy concentration field. To properly model the evolution of a void ensemble in the framework of a phase-field approach, one must recognize that the formation of stable voids due to irradiation damage is physically a first-order transition process in the supersaturated solution of single vacancies. It is well known within the classical nucleation theory that void nucleation occurs when the free-energy barrier separating the metastable and stable states is overcome due to the action of stochastic fluctuations [1]. The associated statistics dictates that the nucleation rate of thermally stable voids typically depends exponentially on temperature and ⇑ Corresponding author.
E-mail address: [email protected] (C.H. Woo).
critical size of void embryo [2–5]. Void nucleation occurring this way is stochastic and local, and is fundamentally very different from that due to spinodal decomposition, for which the phase-field methodology is initially developed [6–11]. Indeed, spinodal decomposition is a deterministic and global process that occurs when the homogeneous solid solution becomes unstable, not metastable. There is no energy barrier to overcome. The absence of an activation barrier corresponds to a zero critical size for the new phase embryo. This means that, unlike the case of metastable solution, in which the detailed behavior of the stochastic fluctuations plays a fundamental and quantitative role, in the spinodal decomposition, stochastic fluctuations only serve as a noise in the vacancy concentration field to trigger the initiation of phase separation process. As a result, the phase separation takes place globally throughout the initially homogeneous solution [10]. In this regard, the adoption of the phase-field approach in modeling irradiation damage processes must be conducted with due care. In a first-order process, the probability of nucleation depends exponentially on the height of the activation barrier, which contains contributions from the energy of the
1359-6454/$36.00 Ó 2012 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.actamat.2012.07.049
A.A. Semenov, C.H. Woo / Acta Materialia 60 (2012) 6112–6119
interface between the critical nucleus and the metallic matrix [2–5]. Even after nucleation, the survival probability of a supercritical void still depends exponentially on the interfacial energy through the rate of vacancy emission [12]. It is thus obvious that an accurate account of this energy contribution is of paramount importance in the usefulness of the phase-field method to model irradiation damage effects. Within the phase-field methodology, the interfacial energy in the free-energy functional is customarily described using a Ginzburg-type gradient energy term with a constant coefficient [6–11]. However, it was shown [1] that, because of the diffusive nature of the void–matrix boundary, the “interfacial” energy must also contain a contribution from the bulk terms in the free energy functional in order to correctly emulate the void thermodynamics. In practical computations, the bulk contribution does depend on the void-size [1], which itself is not a well defined quantity in the phase-field model. This treatment then runs into a self-consistency problem, particularly for small voids. The aim of the present work is to find a practical solution for this issue. 2. Interfacial energy In the phase-field approach, vacancies are assumed to be continuously distributed everywhere in space, including voids, which are treated as specific regions of the vacancy concentration field cv(x), where the local concentration of single vacancies approaches unity, i.e., cv(x) ffi 1. The behavior of the concentration field is assumed to be governed by the thermodynamic potential, conventionally taken to be the phenomenological Cahn–Hilliard-type functional /(cv(x)) [6–11]: Z kBT 2 ½/ðcv ðxÞÞ þ j2 ðrcv Þ dx ð1Þ Uðcv ðxÞÞ ¼ X V where X is the atomic volume, T the absolute temperature, kB the Boltzmann constant. The integral is conducted over the total volume V. The Ginzburg-type gradient energy term with the empirical j-coefficient accounts for the free energy associated with the diffuse void–matrix interface, and is a key ingredient in the phase field model. The bulk free energy density /(cv) in Eq. (1) can be expressed without loss of generality as the dilute ideal solution limit cv lnðcv =ecv1 Þ [13] minus a deviation denoted by u(cv): /ðcv Þ ¼ cv lnðcv =ecv1 Þ uðcv Þ
ð2Þ
Here cv1 is the equilibrium concentration of single vacancies. In the dilute solution limit, i.e., for cv ! 0, u(cv) must satisfy juðcv Þ=cv lnðcv =ecv1 Þj ! 0. In practice, u is usually represented as a phenomenological function with parameters fitted to experimental data. To ensure the existence of a stable “void” phase, /(cv) must also have a global minimum at cv ffi 1.
6113
We have previously shown that in an initially uniform dilute vacancy solution of concentration c0v , i.e., uðc0v Þ ffi 0, the change in the thermodynamic potential DU associated with the formation of a “void” occupying a volume V0 is given by [1] Z kBT DU ¼ ½cv ðxÞ lnðc0v =cv1 Þ þ ð/ðcv ðxÞÞ /ð1ÞÞdx X V0 Z 2 þj2 ðrcv Þ dx ð3Þ V
where /(1) is set equal to (c0v ) to ensure a global minimum of /(cv) at cv = 1. Eq. (3) has been derived assuming that the deviation dcv(x) = cv(x)c0v is small, i.e., |dcv(x)| c0v , everywhere except inside the void volume V0 V, and the law of mass conservation is obeyed, i.e., Z Z dcv ðxÞdx þ dcv ðxÞdx ¼ 0 ð4Þ V V 0
V0
In a supersaturated solution (c0v > cv1) the density of small subcritical voids, serving as embryos for the supercritical voids, is proportional to exp(DU/kBT) [2–5]. In addition, the supercritical voids, which are growing directly from the vacancy solution without the assistance of stochastic fluctuations, nucleate with a rate that is also proportional to exp(DU/kBT), DU being the free energy of formation of a critical void [5]. Thus, to ensure the correct emulation of the kinetics of void precipitation in the phasefield framework, the j-coefficient has to be defined such that DU in Eq. (3) is equal to the corresponding free energy change in the sharp boundary approach. In the latter case, the free energy change associated with the formation of a void with volume V0 is given by [13]: DU ¼ ðV 0 =XÞk B T lnðc0v =cv1 Þ þ rR
ð5Þ
where r is the surface tension (i.e., surface energy density), and R is the void surface area. The volume of a void within the phase-field framework, however, cannot be defined uniquely because of a diffuse character of the interface. This means that DU in Eqs. (3) and (5) do not necessarily correspond to the equivalent voids. Indeed, noting that kBT ln(c0v =cv1 ) is the chemical potential of a vacancy in the dilute ideal solution, the first term of Eq. (5) gives the bulk free energy reduction when a void is formed with m = V0/X vacancies. Similarly, the first term of Eq. (3) may also be interpreted as the free energy change caused R by the formation of a void with volume 0 0 V 0 ¼ m X ¼ V 0 cv ðxÞdx. The spatial layer between the volumes V0 and V 00 may then be treated as part of the diffusive interface between the “void” formed by m0 vacancies and the matrix. Then, using V 00 and R0 (where R0 is the surface area of the volume V 00 ) in Eq. (5) instead of V0 and R, the interfacial energy may be written as Z Z kBT 2 0 2 rR ¼ ð/ðcv ðxÞÞ /ð1ÞÞdx þ j ðrcv Þ dx X V0 V ð6Þ
6114
A.A. Semenov, C.H. Woo / Acta Materialia 60 (2012) 6112–6119
Alternatively, equating Eqs. (3) and (5), we arrive at Z kBT rR ¼ ½ð1 cv ðxÞÞ lnðc0v =cv1 Þ þ ð/ðcv ðxÞÞ /ð1ÞÞ X V0 Z 2 2 dx þ j ðrcv Þ dx ð7Þ V
The diffuse character of the interface also allows Eq. (3) to be considered as the free energy change associated with the formation of a void of size intermediate between m0 and m. Similar to Eqs. (6) and (7), the corresponding expression for the interfacial energy can be easily written down as well. Note that, similar to the gradient term, contributions to the integrals over V0 in Eqs. (6) and (7) mainly come from the interfacial region, where c0v < cv(x) 6 1. These contributions may therefore be considered as part of the interfacial energy. The critical size for large voids, which corresponds to a maximum of DU in Eq. (5), can only be recovered within the phase-field framework when both the gradient term and the bulk term in Eqs. (6) and (7) are taken into account [1]. In such cases, the difference between Eqs. (6) and (7) has been shown to be insignificant [1]. The case of small voids is considered in the next section. As mentioned in the foregoing, equalities like Eqs. (6) and (7) can be used for evaluating the gradient energy coefficient j, which is a model parameter in the phase-field approach. Obviously, values of this coefficient are generally void-size-dependent and, moreover, are not necessarily equal to each other. Nevertheless, there is an alternative treatment, which allows the corresponding dependence to be determined in a unique way. This is the subject of the following section. 3. Gradient energy coefficient In this section, we establish the condition of equivalence between voids in the two approaches and find out the phase-field vacancy distribution cv(x) that corresponds to a void of radius R in the sharp boundary case. Let us consider a vacancy solution in equilibrium with the said void. In this regard, the void is a critical one with a radius that maximizes DU in Eq. (5) [13]. Thus, the solution concentration c0v equals cv1exp(2rX/kBTR), and DU in Eq. (5) can be rewritten in a form that does not explicitly depend on the void radius R, i.e., DU ¼
4pR2 r 16pr3 X2 ¼ 2 3 3ðk B T Þ ln2 ðc0v =cv1 Þ
ð8Þ
Within the framework of the phase-field model, we focus our attention on the vacancy concentration field cv(x), which is in equilibrium with vacancy solution of the same concentration c0v . In equilibrium, the chemical potential l(cv(x)) of the vacancy solution is constant throughout the entire volume, i.e., lðcv ðxÞÞ ¼ k B T lnðc0v =cv1 Þ. Expressing the chemical potential as the functional derivative of the thermodynamic potential in Eq. (1), we may write down the following equation for cv(x):
d/ðcv Þ 2j2 Dcv ðxÞ ¼ lnðc0v =cv1 Þ dcv subjected to the following boundary conditions: ðrcv ðxÞÞx¼0 ¼ 0;
cv ðxÞjjxj!1 ! c0v
ð9Þ
ð10Þ
We note that the void–matrix boundary is not a mathematical boundary in the phase-field model. Assuming spherical symmetry, Eq. (9) reduces to ( ) 2 2 d 4j2 dcv 2 dcv j þ U ðcv ðrÞÞ ¼ ð11Þ dr dr r dr where U ðcv Þ ¼ cv lnðc0v =cv1 Þ /ðcv Þ ¼ cv lnðcv =ec0v Þ þ uðcv Þ 2
ð12Þ
2
Expressing j (dcv/dr) from Eq. (11) and taking into account the boundary conditions (10), DU from Eq. (3) becomes Z 1 8pk B T DU ¼ ðU ðcv ðrÞÞ c0v Þr2 dr ð13Þ X 0 Using the boundary condition (10) at infinity and the definition of the volume V0, outside of which |dcv(r)| is considered to be much smaller than c0v , we may calculate the integrals in Eq. (3) in the limit of r ! 1. The explicit void-volume dependence of Eq. (13) disappears in this case. In a similar way, it can be shown that Eq. (13) can also be derived by substituting the solution cv(r) of the boundary value problem (10)–(12) directly into Eq. (1). Equating DU in Eqs. (8) and (13), we obtain an explicit expression for the j-coefficient: !1=3 !1=3 2p ln c0v =cv1 2rX 2p j¼ ¼R k B T 3Iðc0v ;cv1 Þln2 ðc0v =cv1 Þ 3I c0v ; cv1 ð14Þ where Iðc0v ; cv1 Þ
¼2
Z
1
ðU ðcv ðqÞÞ c0v Þ4pq2 dq
ð15Þ
0
The last equality in Eq. (14) is derived from the expression for the critical void radius: R = 2rX/[kBTln(c0v =cv1 )]. Making the transformation r = jq in Eqs. (10)–(12), it can be seen that the solution cv(q) of the boundary value problem has no explicit dependence on j, and consequently, neither does the integral I(c0v , cv1). The gradient energy coefficient j can then be calculated directly using Eq. (14). Indeed, in terms of the variable q, Eq. (11) can be transformed into the following equation by integrating once: " # " # 2 2 dcv dcv ¼ ½U ðcv ðqiþ1 ÞÞ dq dq q¼qiþ1
q¼qi
U ðcv ðqi ÞÞ
Z
qiþ1 qi
2 4 dcv dq0 q0 dq0
ð16Þ
A.A. Semenov, C.H. Woo / Acta Materialia 60 (2012) 6112–6119 10 0
c v∞ = 10
0
5 0
0
0 -10
-φ
ð17Þ
ð19Þ
To proceed further, the explicit form of the free energy density /(cv) has to be specified. Following Ref. [1], we write down this function as /ðcv Þ ¼ cv lnðcv =ecv1 Þ þ ½n lnðcv1 Þ þ ðn þ 1Þð1 c0v Þcnv ½ðn 1Þ lnðcv1 Þ þ nð1 c0v Þcnþ1 v
-9
Cv = 10
0
For the asymptotic cv(q) to satisfy Eq. (16), it must fulfill the condition 2
-8
Cv = 10
Cv = 10
Here q0 is a point such that 0 < dcv ðq0 Þ=c0v 1. Differentiating Eq. (17), we obtain qffiffiffiffiffiffiffi i dcv ðq Þ h ðdcv =dqÞq¼q0 ¼ pffiffiffiffiffiffi0ffi 1 þ 2c0v =q0 2c0v rffiffiffiffi rffiffiffiffi c0v dcv ðq0 Þ c0v ffi ð18Þ << 0 2 cv 2
ðdcv =dqÞq¼q0 >> juðc0v Þj
-7
Cv = 10
dcv ðqÞ ¼ cv ðqÞ c0v
qffiffiffiffiffiffiffi dcv ðq0 Þq0 expððq q0 Þ= 2c0v Þ > 0 ffi q
-6
Cv = 10
−10
U
The latter provides a recursive relation via which dcv/dq and hence cv(q) can be readily obtained as a function of cv(0) using a finite difference scheme, with dcv/dq = 0 at q = 0. The value of cv(0) can then be varied to obtain the solution that satisfies the asymptotic behavior of cv(q) at q ! 1. In the immediate vicinity of c0v (i.e., for jcv c0v j c0v ), where U(cv) maximizes, U(cv) can be expanded in the form U ðcv Þ ffi c0v ðdcv Þ2 =ð2c0v Þ. The solution of Eqs. (10) and (11) then gives the following asymptotic description of cv(q):
6115
ð20Þ
where n P 2, which provides the correct dilute ideal solution limit. For inequality (19) to be valid even for the vacancy concentrations near the melting temperatures, i.e., for c0v 105–104, the exponent n in Eq. (20) should satisfy n P 2.5. 4. Results and discussions Profiles of the “potential” U(cv) calculated for various values of c0v are shown in Fig. 1. For larger values of n, the curvature of U(cv) near the maximum at cv = 1 increases, its minimum deepens [1], and the diffuse interface becomes thinner, similar to the case of the planar surface considered earlier [1]. Thus, to minimize the number of grid points for the numerical calculation, we take the smallest possible value for n, which is 2.5. For n > 2, the behavior of the free energy density /(cv) starts to deviate significantly from the ideal solution (/ cvln(cv/ecv1)) when cv > 102 [1]. For weaker solutions, i.e., when c0v cv(q) < 102, we have shown in the Appendix that the following relation holds for the spatial behavior of cv(q):
-5 0.0
0.1
0.2 0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
cv Fig. 1. “Potential” U(cv) at cv1 = 1010 as a function of vacancy concentration for the exponent n = 2.5 and various values of vacancy supersaturation c0v . Arrow indicates the drop in the “potential energy” when the “motion” takes place from cv ffi 1 at r = 0 to cv = c0v at r ! 1.
Dq ¼ q0 q sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffi cv ðqÞ 1 0 ð1 þ ln 2 ðc ðqÞ=c ÞÞ 2c0v v v lnðcv ðqÞ=c0v Þ lnðdcv ðq0 Þ=c0v Þ
ð21Þ
which describes a very sharp, albeit diffuse, interface. Indeed, within a spatial range of Dr = jDq < 101j, where j has a value of at most a few lattice constants [1], the vacancy concentration in the interface region could drop by eight orders of magnitude from a value of 102. Accumulation of vacancies in the void is driven by the difference between the actual concentration of vacancies and the concentration of the equilibrium vacancy solution c0v [4]. Under reactor operating conditions, c0v is often smaller than 1010 [1]. If we assume, for the sake of argument, that a feasible upper bound for the actual vacancy concentration could be as large as the melting-temperature equilibrium concentration of 105–104, then, from Eq. (21), the vacancy concentration in the interface has to drop five to six orders of magnitude within a very narrow spatial range Dq 103–102. Thus, to correctly represent in the phase-field approach the concentration gradient encountered in irradiation-induced void growth, numerical calculations would have to be performed on an extremely fine sub-nanometer spatial scale, a formidable task for the computation. Eq. (16) is numerically solved using a finite difference scheme, with the left-hand-side boundary value cv(0) set near 1, and the derivative (dcv/dq)q=0 = 0. The solution is repeated, varying values of cv(0), until the asymptotic boundary condition (18) is fulfilled. To ensure that the calculated vacancy concentration field cv(q) is the solution of
6116
A.A. Semenov, C.H. Woo / Acta Materialia 60 (2012) 6112–6119 0
0
10
10
-1
10
-1
10
C v∞ = 10
-2
C v∞ = 10
-2
−12
10
−10
10
-3
10
cv
-6
Cv= 10
0
-7
Cv= 10
0
-8
Cv= 10
0
-9
Cv = 10
cv
0
-3
10
0
0
-9
-8
Cv= 10
1.0
1.5
Cv= 10
-4
10
0
0
-10
-11
Cv = 10
Cv= 10
-5
10
-4
10
-6
10 -5
10
-7
10
-8
-6
10
0.0
0.5
1.0
1.5
-2
2.0
2.5
10
3.0
0.0
0.5
-3
10
10
10
-4
10
-4
10
3.0
-7
10
-3
2.5
-6
-3
10
10
2.0
-5
10 -4
0
10
-6
-6
10
Cv = 10
-8
-5
0
10
-9
Cv = 10
0
-8
Cv = 10
0
-11
Cv = 10
-9
-6
10
10
-7
10
10
-5
10
-10
-7
10
10
-8
10
-8
-6
1.00
1.02
1.04
1.06
1.08
10 1.10
1.05
2.908 2.912 2.916 2.920 2.924
Eq. (16) with boundary conditions (10), the consistency of the free energy changes D/ calculated using both Eqs. (3) and (13) is verified. In Figs. 2 and 3, we present the calculated concentration profiles for two different equilibrium concentrations cv1 and various levels of vacancy supersaturation c0v =cv1 . For copper and molybdenum, parameters are chosen to give values of critical radii between 0.3 and 1.5 nm, which are typical for these materials under irradiation for a wide range of temperatures [12]. The spatial positions of the sharp void–matrix boundary can be calculated from Eq. (14). They are found to be within the diffuse interface, which is consistent with our general treatment given in Section 2. Accordingly, the sharp boundary position R00 can also be estimated directly from the phase-field calculations as 1=3 3 0 n cv ; cv1 R00 =j ¼ ð22Þ 4p
nðc0v ; cv1 Þ ¼ 4p
Z
1
cv ðqÞq2 dq
ð23Þ
0
is the number of vacancies in the “void” in the phase-field approach. The ratio of R00 =R is presented in Table 1, from
1.06
1.07
1.08
10 3.1400 3.1402 3.1404 3.1406 3.1408 3.1410
Fig. 3. Same as Fig. 2, but at cv1 = 1012.
Fig. 2. Concentration profiles cv(r/j) at cv1 = 1010 for various values of vacancy supersaturation c0v . Arrows show spatial positions of the sharp boundary calculated with Eq. (14). Two inserts demonstrate the behavior of concentration profiles on finer spatial scales.
where
-11
10
-9
10
Table 1 Ratio R00 =R. cv1 10
10 1012
c0v =cv1 101
102
103
104
0.992 0.993
0.999 0.995
1.022 1.007
1.066 1.031
which one may conclude that Eq. (22) in the phase-field model does provide a good approximation for the sharp boundary position. From the inserts in both figures the validity of the approximation (21) can be also verified. In Table 2, values of the j-coefficient calculated with Eq. (14) are listed. The value of j is relatively constant, almost independent of the void size R. Depending on temperature, it is maintained at about one or two lattice constants a. This is consistent with our estimation in the foregoing section, from which we have concluded that the diffuse void interface is very sharp and has to be realized on a very fine spatial scale. For such thin interfaces the integrals in Eq. (3) can be separated with sufficient accuracy into bulk and surface contributions. This explains why values of R and R00 are close to each other. It also explains why the j-value only weakly depends on the void size, and why the j-value obtained in the plane boundary case (c0v =cv1 ¼ 1) [1] can be used in studies even for voids with
A.A. Semenov, C.H. Woo / Acta Materialia 60 (2012) 6112–6119 Table 2 Gradient term coefficient j/(2rX/kBT). cv1 10
10 1012
c0v =cv1 1.0
101
102
103
104
0.173 0.156
0.169 0.154
0.168 0.152
0.169 0.151
0.173 0.153
radii as small as one lattice constant. Indeed, when c0v =cv1 ¼ 106 , which corresponds to voids formed with just a few vacancies in the case of molybdenum and copper, the width of the interface becomes comparable with the void radius. Computer calculations show that only then the jcoefficient will start to deviate significantly from the plane boundary value. It is worth noting that in the phase-field theory of solidification, the Cahn–Hilliard analysis also often relies on the constancy of the j-coefficient [14–20]. Physically, the gradient term in the total free energy expression represents spatial correlation effects, and the j-coefficient in this term is related to the second derivatives of the direct correlation function of the liquid evaluated at the reciprocal lattice vectors of the corresponding crystal [14,21,22]. Although the j-coefficient is constant when the interpolated bulk free energy density / is assumed to be in a polynomial form (piecewise parabolic [15–17], triple-parabolic [18], quartic [19,20]), yet the corresponding effective solid–liquid interfacial tension r depends on the size of the critical nucleus [17–20]. This dependence is especially strong at high undercoolings, when the critical nucleus has a size comparable to the thickness of the interface. The sharp-interface assumption of the classical nucleation theory is invalid in such cases. As a result, in agreement with the experimental observations, the nucleation rates of crystalline structures in highly undercooled liquids calculated within the diffuse interface approach is orders of magnitude higher than those from the classical nucleation theory [17,19,20]. It also reflects the fact that the solid– liquid interface has a finite width and is diffuse on the molecular scale, extending over several molecular layers [11]. On the other hand, the interface considered in the present work remains thin compared to the size of the nucleus, even for vacancy supersaturations as high as c0v =cv1 ¼ 104 . As a result, the constancy of the surface tension in the sharp interface approach translates into the constancy of the gradient energy coefficient, and vice versa. Similar to Refs. [14,15], the diffuse nature of the interface in the present case does not significantly affect the interfacial tension because the spinodal point where vacancy solution becomes unstable (as opposed to being metastable) remains sufficiently far away. Indeed, the second derivative d2// (dcv)2 becomes negative above the spinodal point [1]. At the same time, even under in-reactor conditions with continuous production of vacancies, the vacancy solution is still dilute and ideal. Thus, the concentrations c0v of vacancy solutions considered here are far below those (c0v > 102 [1]) for which polynomial-type deviation u(cv) from the free
6117
energy density cv lnðcv =ecv1 Þ of the dilute solution is important. We will show in a forthcoming paper that the dynamical behavior of an individual void, well-known from the ratetheory treatment of void evolution [23], can be quite accurately reproduced in the framework of the phase-field model considered here. However, the characteristic spatial scales of the void–metal diffuse interface found in the present analysis presents a challenge to modeling the evolution of a void ensemble under irradiation. As a possible solution to this problem, one may consider smearing out the interface and then using the method of matched asymptotic expansion [11,24] to relate the phase-field model to the sharp interface approach. Since the vacancy concentration in the solution is very low, the void growth is very slow compared to the relaxation time of the vacancy diffusion field, so that the diffusion field in front of the interface can be considered at steady state. Thus, in the case of voids the interface Pe´clet number [11,24] is vanishingly small, which will significantly simplify the corresponding matched asymptotic analysis. 5. Summary and conclusions In this paper, we consider phase-field modeling of void development as a first-order phase transition in a metastable vacancy solution. The thermodynamic barrier against the nucleation of thermally stable growing voids is caused by the formation of the void–metal interface. Under the conventional framework, in which the Ginzburg-type gradient energy term is used to represent the interfacial energy, the associated coefficient is a model parameter to ensure the equivalence of the physical and the emulated void. Establishing equivalence between voids in the phase-field model with the diffuse interface and physical voids with the sharp boundary allows the diffuse model to elucidate the concentration of vacancies and structure of voids in a thermodynamic context, something that molecular dynamics cannot do. This forms the focus of the present study. In general, it is expected that the solution of the voidgrowth kinetic equation would be complicated by the void-size dependence of the gradient-energy coefficient. However, we found that for a fixed temperature this parameter may be treated, to a good approximation, as a constant, not only for large voids, but even for voids with sizes comparable with a lattice constant. This result can be traced to the sharpness of the diffuse interface, which has a characteristic width of less than one lattice constant, despite its diffuse character. Thus, the actual spatial position of the sharp boundary can be interpolated very well by its diffuse analogue, even in the case of very small voids. The sharpness of the diffuse interface is a direct result of the behavior of the thermodynamic potential in the dilute ideal solution limit, and the requirement that the Cahn–Hilliard-type interpolative functional must satisfy this limit. This requirement is absolutely necessary because under real experimental conditions, even under
6118
A.A. Semenov, C.H. Woo / Acta Materialia 60 (2012) 6112–6119
continuous production of vacancies by irradiation, the supersaturated vacancy solution in which void nucleation and growth take place is still dilute and ideal. Acknowledgement This project is initiated and funded by Grant 534409 from the Hong Kong Research Grant Commission, to which the authors are thankful.
and using the relation M(a, a; z) = exp(z), Eq. (A5) can be expressed through the tabulated function M(1/2, 1/2; z): sffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi" cv ðqÞ c0v q0 q ffi 2 1 M 0 cv ðqÞ lnðcv ðqÞ=cv Þ qffiffiffiffiffiffiffi 1 cv ðqÞ dcv ðq0 Þ 0 ln 2c 1=2; 1=2; ln v 2 c0v c0v
ðA7Þ
Using the asymptotical behavior CðcÞ z ac 1 e z ½1 þ Oðjzj Þ CðaÞ
Appendix A
Mða; c; zÞ ¼
In the one-dimensional case Eq. (9) can be written as 2 dcv þ U ðcv ðqÞÞ ¼ const ðA1Þ dq
where C(a) is the gamma function, we arrive at Eq. (17). In the case of spherical symmetry Eq. (9) has the form
The boundary condition at q ! 1 gives that const ¼ c0v . Thus, the solution can be presented in the following integral form: Z cv1 dcv pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q0 q1 ¼ 0 cv U ðcv Þ cv0 Z cv1 1=2 ffi ðc0v þ cv lnðcv =ec0v ÞÞ dcv ðA2Þ cv0
Here cv0 = cv(q0), and cv1 = cv(q1). The second equality in Eq. (A2) is valid when cv 6 102. If both d1 = dcv1/c0v and d0 = dcv0/c0v 1, Eq. (A2) reduces to qffiffiffiffiffiffiffi q0 q1 ¼ 2c0v lnðd1 =d0 Þ ðA3Þ On the other hand, for dcv1/c0v 1 and c0v cv(q) 6 102, a major contribution to the integral (A2) comes from the interval of integration, where the difference (c0v U(cv)) can be approximated by cvln(cv/c0v ). Then, using the substitution y = 0.5ln(cv/c0v ), we have that qffiffiffiffiffiffiffi Z 0:5 lnðcv ðqÞ=c0v Þ y e dy q1 q ffi 2c0v pffiffiffi y 0:5 lnð1þd1 Þ qffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi lnðcv ðqÞ=c0v ÞMð1=2; 3=2; 0:5 lnðcv ðqÞ=c0v ÞÞ ¼ 2 c0v p lnð1 þ d1 ÞMð1=2; 3=2; 0:5 lnð1 þ d1 ÞÞ ðA4Þ where M(a, c; z) is the confluent (degenerate) hypergeometric function [25,26]. Noting that the second term in Eq. (A4) is much less than the first one, and that the correct expression for the sum of Eqs. (A3) and (A4) has to be independent of d1, we may write down the interpolation for both equations as follows: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q0 q ¼ 2 c0v lnðcv ðqÞ=c0v ÞMð1=2; 3=2; 0:5 qffiffiffiffiffiffiffi lnðcv ðqÞ=c0v ÞÞ 2c0v lnðdcv ðq0 Þ=c0v Þ ðA5Þ From the recurrence relation [15], cMða; c; zÞ cMða 1; c; zÞ ¼ zMða; c þ 1; zÞ
ðA6Þ
2
d 2 cv 4 dcv dU ¼ þ dcv dq2 q dq
ðA8Þ
ðA9Þ
It can be approximated by the foregoing one-dimensional equation when dU 4 dcv ðA10Þ q dq dcv According to Eq. (A4), the derivative dcv/dq is approximately given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 dcv 4 cv ðqÞ lnðcv ðqÞ=c0v Þ ffi q dq q 2Dq lnðcv ðqÞ=c0v Þ ffi ðA11Þ q Since dU/dcv = ln(cv/c0v ) for c0v < cv 6 102, the inequality (A10) is fulfilled when the void radius is much larger than the characteristic width of the interface. References [1] Semenov AA, Woo CH. J Nucl Mater 2011;411:144. [2] Zeldovich YaB. Sov J Exp Theor Phys 1942;12:525. [3] Frenkel YaI. Kinetic theory of liquids. Moscow: Academy of Science of the USSR; 1959. [4] Lifshitz EM, Pitaevskii LP. Physical kinetics. Oxford: Pergamon Press; 1981 [Chapter XII]. [5] Slezov VV, Schmelzer J, Tkatch YaY. J Chem Phys 1996;105:8340. [6] Cahn JW, Hilliard JE. J Chem Phys 1958;28:258. [7] Hillert M. Acta Metall 1961;9:525. [8] Cahn JW. Acta Metall 1961;9:795. [9] Cahn JW. Acta Metall 1962;10:179. [10] Skripov VP, Skripov AV. Uspekhi Fiz Nauk 1979;128:193 [in Russian]. [11] Emmerich H. Adv Phys 2008;57:1. [12] Semenov AA, Woo CH. Phys Rev B 2002;66:024118. [13] Landau LD, Lifshitz EM. Statistical physics. Oxford: Pergamon Press; 1980. [14] Harrowell P, Oxtoby DW. J Chem Phys 1984;80:1639. [15] Bagdassarian CK, Oxtoby DW. J Chem Phys 1994;100:2139. [16] Wild R, Harrowell P. Phys Rev E 1997;56:3265. [17] Gra´na´sy L, Iglo´i F. J Chem Phys 1997;107:3634. [18] Gra´na´sy L, Oxtoby DW. J Chem Phys 2000;112:2399. [19] Gra´na´sy L, Bo¨rzo¨nyi T, Pusztai T. Phys Rev Lett 2002;88:206105. [20] Gra´na´sy L, Pusztai T, James PF. J Chem Phys 2002;117:6157. [21] Haymet ADJ, Oxtoby DW. J Chem Phys 1981;74:2559.
A.A. Semenov, C.H. Woo / Acta Materialia 60 (2012) 6112–6119 [22] Haymet ADJ, Oxtoby DW. J Chem Phys 1982;76:6262. [23] Brailsford AD, Bullough R. J Nucl Mater 1972;44:121. [24] Edler KR, Grant M, Provatas N, Kosterlitz JM. Phys Rev E 2001;64:021604.
6119
[25] Gradshteyn IS, Ryzhik IM. Table of Integrals, series and products. New York: Academic Press; 1965. p. 1058. [26] Abramowitz M, Stegan IA, editors. Handbook of mathematical functions. New York: Dover; 1965. p. 504.