- Original article
- Open access
- Published:
An analysis of two classes of phase field models for void growth and coarsening in irradiated crystalline solids
Materials Theory volume 2, Article number: 1 (2018)
Abstract
A formal asymptotic analysis of two classes of phase field models for void growth and coarsening in irradiated solids has been performed to assess their sharp-interface kinetics. It was found that the sharp interface limit of type B models, which include only point defect concentrations as order parameters governed by Cahn-Hilliard equations, captures diffusion-controlled kinetics. It was also found that a type B model reduces to a generalized one-sided classical Stefan problem in the case of a high driving thermodynamic force associated with the void growth stage, while it reduces to a generalized one-sided Mullins-Sekerka problem when the driving force is low in the case of void coarsening. The latter case corresponds to the famous rate theory description of void growth. Type C models, which include point defect concentrations and a non-conserved order parameter to distinguish between the void and solid phases and employ coupled Cahn-Hilliard and Allen-Cahn equations, are shown to represent mixed diffusion and interfacial kinetics. In particular, the Allen-Cahn equation of model C reduces to an interfacial constitutive law representing the attachment and emission kinetics of point defects at the void surface. In the limit of a high driving force associated with the void growth stage, a type C model reduces to a generalized one-sided Stefan problem with kinetic drag. In the limit of low driving forces characterizing the void coarsening stage, however, the model reduces to a generalized one-sided Mullins-Sekerka problem with kinetic drag. The analysis presented here paves the way for constructing quantitative phase field models for the irradiation-driven nucleation and growth of voids in crystalline solids by matching these models to a recently developed sharp interface theory.
Introduction
Irradiation drives complex microstructure evolution in solids by producing large non-equilibrium densities of point defects (Olander, 1976; Was, 2017; Brailsford & Bullough, 1972). Features such as dislocation loops and voids are often observed in irradiated solids. The presence of voids, the microstructure type of interest here, strongly affects the thermal and mechanical properties of irradiated materials. As such, many investigations were conducted to understand the void formation and growth processes (Olander, 1976; Was, 2017; Brailsford & Bullough, 1972; Krishan, 1982; Dubinko et al., 1989; El-Azab et al., 2014; Hochrainer & El-Azab, 2015; Millet & Tonks, 2011a). In analogy to the theory of sintering (Lifshitz & Slyozov, 1961; Rahaman, 2003), theoretical models of void formation and growth usually treat voids as a second phase of vacancies that nucleate in a metastable supersaturated solid matrix with the metastability sustained by the production and accumulation of point defects. These models thus consider the process as an example of non-equilibrium first-order phase transition in driven systems. The nucleation models predict a nucleation barrier for the void formation that decreases as the supersaturation of vacancies (interstitials) increases (decreases) (Olander, 1976; Was, 2017; Katz & Wiedersich, 1971; Russell, 1971; Mayer & Brown, 1980). The void growth models can be classified into two types; the rate theory models (Olander, 1976; Was, 2017; Brailsford & Bullough, 1972; Krishan, 1982; Dubinko et al., 1989) and the field-theoretic, spatiotemporal models (El-Azab et al., 2014; Hochrainer & El-Azab, 2015; Millet & Tonks, 2011a). In the former, a test void in an effective homogeneous medium is considered to represent the average kinetics of the system. Moreover, the process of void evolution is treated as diffusion-controlled, and the void growth equation is usually derived assuming quasistatic diffusion in the solid matrix (Olander, 1976; Was, 2017; Brailsford & Bullough, 1972; Krishan, 1982; Dubinko et al., 1989). Some rate theory models known as cluster dynamics models take into consideration the size distribution of the evolving void system (Was, 2017). These models will not be mentioned further in the current discussion. In field theoretic models, local balance equations and void growth are treated within the principles of irreversible thermodynamics and the system of voids is resolved in space and time (De Groot & Mazur, 1962).
The field-theoretic models can be further classified into two categories, sharp- and diffuse-interface models. In a sharp-interface model developed recently by Hochrainer and El-Azab (2014; 2015), the void surface is treated as a singular surface across which jump conditions are applied. The diffuse-interface models, also called the phase field models, consider the void surface to be diffuse, i.e., has a finite width, in order to circumvent the numerical difficulties associated with solving the moving boundary problem of evolving voids (El-Azab et al., 2014; Hochrainer & El-Azab, 2015). However, the void surface is atomically-sharp and hence one must ensure that diffuse-interface models recover the sharp-interface description as the interface thickness vanishes. As of now, two different classes of phase field models for void growth have appeared in literature (Yu & Lu, 2005; Hu et al., 2009; Hu & Henager, 2009; Hu & Henager, 2010; Li et al., 2010; Semenov & Woo, 2012; Rokkam et al., 2009; Millett et al., 2009; Millett et al., 2011c; Millett et al., 2011b; Li et al., 2013; Xiao et al., 2013). Following the terminology of the field-theoretic approach set by physicists for modeling heterogeneous materials (Hohenberg & Halperin, 1977; Binder, 1987; Provatas & Elder, 2010), these models are conveniently classified as models of type B and type C. Such models, while intended to represent one and the same phenomenon, void nucleation and growth, do appear physically and mathematically different. Therefore, an assessment of these phase field models showing their relations to the sharp-interface physics and the classical rate theory models is desirable for making further progress.
The method of matched asymptotic expansions (Provatas & Elder, 2010; Pego, 1989; Dai & Du, 2012; Elder et al., 2001; Fife, 1992; Emmerich, 2008; Garcke et al., 2004; Ahmed et al., 2016; Caginalp, 1989) is used here to deduce the sharp-interface limits of type B and type C phase field models for void evolution in irradiated solids. The asymptotic analysis serves two purposes: it demonstrates the consistency of the phase field models in terms of their physics content, and, in the process, it helps to relate their parameters to the regular thermodynamic and kinetic parameters that appear in the sharp-interface models. This is critical for making the phase field models quantitative. The asymptotic analyses performed here can be considered as generalizations of those by Pego (1989), Dai and Du (2012), and Elder et al. (2001) to the case of driven multi-component systems.
Our analysis concludes that phase field models of type B, which utilize the point defect concentrations as the only order parameters, can describe diffusion-controlled kinetics. Moreover, in the low driving forces limit (coarsening stage) they are equivalent to the rate theory models. On the other hand, phase field models of type C, which couple Cahn-Hilliard type equations (Cahn, 1961) for the local balance of point defects with an Allen-Cahn equation (Allen & Cahn, 1979) for the motion of the void surface, are able to take into account the attachment kinetics of point defects to the void surface. The attachment kinetics of the point defects to the void surface affects the overall void growth rate, as was shown in the numerical simulations of the sharp-interface model (El-Azab et al., 2014; Hochrainer & El-Azab, 2015). It is shown that the additional time-dependent Allen-Cahn (Ginzburg-Landau) equation in model C acts as the interfacial constitutive law that ensures positive interfacial entropy production associated with the surface motion due to its reactions with point defects (El-Azab et al., 2014; Hochrainer & El-Azab, 2015). The connections between the coarsening limits of phase field models B and C and the classical mean-field Lifshitz-Slyozov-Wagner theories of Ostwald ripening are further addressed here (Lifshitz & Slyozov, 1961; Rahaman, 2003; Mullins & Sekerka, 1963; Niethammer, 2000; Dai et al., 2010). These limits are all important in fixing the phase field model parameters and understanding their results.
The theoretical models of void growth in irradiated solids are reviewed first. The asymptotic analysis of the different phase field models for void growth is then presented, followed by a summary and concluding remarks.
Models for void growth in irradiated solids
Rate theory description
Rate theory assumes uniform background concentrations of point defects and sinks in the solid matrix in which an average void is immersed. The void growth process is viewed as a diffusion-controlled process. A typical rate theory model consists of three coupled ordinary differential equations (Eqs. (1.1) and (1.2)), the first two of which represent the balance of point defect concentrations and the third is the diffusion-controlled void growth equation (Olander, 1976; Was, 2017; Brailsford & Bullough, 1972; Krishan, 1982; Dubinko et al., 1989),
In the above, cv and ci are the average (fractional) vacancy and interstitial concentrations in the solid, with the superposed dot referring to their time rates of change, Pv and Pi are the respective production terms, kvi is a rate constant for vacancy-interstitial recombination, kvs and kis are the rate constants for defects reaction with sinks of average concentration cs, R is the void radius, Ω is the atomic volume, γ is the surface energy, \( {c}_{\mathrm{v}}^{\mathrm{R}} \) and \( {c}_{\mathrm{i}}^{\mathrm{R}} \)are vacancy and interstitial concentrations in equilibrium with a void of radius R, \( {c}_{\mathrm{v}}^{\mathrm{eq}} \) and \( {c}_{\mathrm{i}}^{\mathrm{eq}} \)are the equilibrium vacancy and interstitial concentrations with a flat surface, and Dv and Di are the diffusion coefficients of vacancies and interstitials, respectively and k and T are Boltzman constant and temperature. The rate theory has two major shortcomings. The first is the treatment of the solid as a homogeneous medium. Under irradiation, high gradients of point defect concentrations exist in the solid matrix, particularly near the surfaces of sinks such as free surfaces (or voids), grain boundaries, and dislocations. Describing such an effective homogeneous medium is complicated and never exact. The second shortcoming pertains to the assumption of diffusion-controlled growth (Eq. (1.2)). The process of void growth is not necessarily diffusion-controlled, and the point defect concentrations at the void surface usually deviate from their equilibrium values (Eq. (1.3)). In fact, as was recently shown by Hochrainer and El-Azab (2014; 2015), the overall growth kinetics depends on the interaction between the point defects and the void surface.
Sharp-interface description
Spatiotemporal models of void evolution (El-Azab et al., 2014; Hochrainer & El-Azab, 2015; Millet & Tonks, 2011a) treat voids as a second phase that nucleates in a supersaturated matrix (Lifshitz & Slyozov, 1961; Rahaman, 2003). The sharp-interface description of void growth may thus be considered a generalization of the classical sharp-interface models of particle growth from a supersaturated matrix, as in solidification and precipitation. However, the classical precipitation/solidification models are examples of phase transitions in non-driven systems, i.e., systems relaxing toward equilibrium from close thermodynamic states, while the void nucleation and growth is an example of non-equilibrium phase transition under an external driver, irradiation. It is the ongoing production of point defects that renders the solid matrix metastable all the time and causes the nucleation and growth of voids. The lack of sharp-interface models for void growth under irradiation motivated Hochrainer and El-Azab (2014; 2015) to introduce an elaborate one. Nevertheless, it is useful to first review the classical models of precipitation since they represent limiting cases of the generalized sharp- and diffuse-interface models for void growth, before introducing the sharp interface model of voids reported in (El-Azab et al., 2014; Hochrainer & El-Azab, 2015).
The classical sharp-interface models of precipitation consider the interface between the precipitate and the parent phase to be a singular surface (Lifshitz & Slyozov, 1961; Rahaman, 2003). Therefore, in addition to the regular balance and constitutive laws for phases around the interfaces, interfacial balance and constitutive laws are required (Provatas & Elder, 2010; Pego, 1989; Dai & Du, 2012; Elder et al., 2001; Fife, 1992; Emmerich, 2008; Garcke et al., 2004; Ahmed et al., 2016; Caginalp, 1989). For a particle of phase α growing from a supersaturated matrix of phase β by mass diffusion, the isothermal growth kinetics is captured by
In the above, c k is the fractional concentration of species k, satisfying the condition \( \sum \limits_k{c}_k=1 \), J k is the flux of the corresponding species, μ k is its chemical potential, M k is its mobility. f is the Helmholtz free energy density, κ is the interface curvature, and λ is the interface relaxation constant (λ−1 is the interface kinetic coefficient); it is also sometimes called the coefficient of kinetic drag (Niethammer, 2000; Dai et al., 2010). The notation \( {\left[\bullet \right]}_{\alpha}^{\beta } \) refers to the jump across the interface, Γ, of the quantity in brackets, and n is the unit normal to the interface, pointing from Ω− to Ω+, which, respectively refer to α and β phases. Eqs. (2.1) and (2.2) represent the mass balance and constitutive laws for each species in both phases. Eq. (2.3) is a statement of the continuity of the chemical potential at the interface. Eq. (2.4) and Eq. ((2.5) are the interfacial balance and constitutive laws, respectively. Eq. (2.4) is also known as the Stefan jump condition while Eq. ((2.5) is called the dynamical Gibbs-Thompson equation (Provatas & Elder, 2010; Pego, 1989; Dai & Du, 2012; Elder et al., 2001; Fife, 1992; Emmerich, 2008; Garcke et al., 2004; Ahmed et al., 2016; Caginalp, 1989). The system of eqs. (2.1)-(2.5) is of course to be supplemented with appropriate initial conditions and boundary conditions on the outer boundary ∂Ω. For a growing particle with fixed concentration, the dynamical system becomes one-sided, i.e., diffusion takes place only in the metastable parent phase.
It is worth noting that the dynamical Gibbs-Thompson relation (Eq. ((2.5)) reduces to the static (or equilibrium) Gibbs-Thompson relation for a stationary interface. Moreover, for the case of a flat interface in a binary system, it recovers the common tangent (Maxwell) construction rule. While the dynamical Gibbs-Thompson condition was first introduced in an ad hoc manner based on experimental data (Langer, 1980), it is now known that it arises as a consequence of the non-negative entropy production expression of the second law of thermodynamics. Specifically, the dynamical Gibbs-Thompson relation can be obtained by requiring the bulk and interfacial entropy production to be non-negative independently and assuming a linear relationship between the normal interface velocity (the normal fluxes at the interface) and the corresponding driving thermodynamic forces (El-Azab et al., 2014; Hochrainer & El-Azab, 2015). For the case of infinitely fast interface kinetics (λ−1 → ∞, λ → 0), the velocity term vanishes and we get the classical Gibbs-Thompson condition as the interfacial constitutive law indicating a diffusion-controlled growth. In other words, the diffusion-controlled models ignore the interface kinetics and assume zero interfacial entropy production. When the velocity term is considered, and hence the interface kinetics is taken into account, the dynamical Gibbs-Thompson relation states that the chemical potentials and concentrations of each species deviate from their equilibrium thermodynamic values at the interface.
The classical sharp-interface models of particle growth based on the dynamical system (2.1)-((2.5) are often called Stefan or Mullins-Sekerka models (Provatas & Elder, 2010; Caginalp, 1989; Mullins & Sekerka, 1963; Niethammer, 2000; Dai et al., 2010). Stefan models are valid for the high supersaturation case and hence describe the growth stage. Mullins-Sekerka models are valid for the low supersaturation case since they assume quasistatic diffusion in the bulk phases and are hence more suitable for describing the coarsening stage. Furthermore, each model is sometimes called classical or modified (with kinetic drag) depending on whether the classical equilibrium Gibbs-Thompson condition or the dynamical Gibbs-Thompson condition is considered as the interfacial constitutive law. A schematic illustration of these different limits is shown in Table 1. As will be shown below, the limits of the sharp- and diffuse interface models of void growth under irradiation can be considered as generalizations of these classical models.
As mentioned previously, the lack of a general sharp-interface model of particle growth in driven systems has motivated the development of such model for single component solids (El-Azab et al., 2014; Hochrainer & El-Azab, 2015). This model serves as the sharp interface limit of the diffuse interface model for the case of voids. Ignoring the contribution of surface diffusion to the motion of the void surface, the model is summarized here. The local balance equations for the point defect concentrations are
In the above, c α (α = i, v) denotes the site fraction of interstitials (i) and vacancies (v). Interstitials are assumed to appear in dumbbell (or split) configuration where two atoms share a regular lattice site. Hence, a specific lattice site is occupied with a regular atom, a vacancy, or a dumbbell interstitial such that ca + cv + ci = 1, where ca is the regular atom site fraction. Therefore, in the fixed lattice frame, only the point defect fluxes are independent, i.e., the atom diffusion is mediated by the diffusion of point defects. P α is production rate of defects due to irradiation. It is proportional to the irradiation dose rate and is stochastic in space and time. Moreover, it usually has a characteristic spatial structure (cascade structure), commonly assumed to be a core-shell structure with vacancies concentrated in the cascade center and interstitial are distributed at the periphery (Olander, 1976; Was, 2017). Riv is the recombination rate of interstitials and vacancies, which is proportional to the product of the concentrations of the two reacting species as in Eq. (2.1). However, more general forms can be utilized as was discussed in (El-Azab et al., 2014; Hochrainer & El-Azab, 2015). J α refers to the diffusive fluxes of defects,
with μ α being the coresponding chemical potentials and the mobility M α being generally dependent on concetrations. These chemical potentials are obtained as the derivatives of the free energy with respect to the corresponding species concentration. The Helmholtz free energy density (per lattice site) of the defective lattice is given by (Hochrainer & El-Azab, 2015)
f a is the free energy of an atom in the perfect (defect-free) crystal, fi and fv are the free energy of formation of interstitials and vacancies, respectively, k B is the Boltzmann constant, and T is the temperature.
The balance of mass across the void surface combines the Stefan (jump) conditions for vacancies and interstitials. This gives the surface velocity as follows
The kinetics of the point defects attachment to the void surface is included as interfacial constitutive laws (or boundary conditions). These boundary conditions require the normal fluxes of point defects at the void surface to be matched to the reaction rates of point defects with the void surface. Using the transition state theory expressions for the defect-surface reactions, the normal fluxes at the void surface are thus expressed as
Here, δ is the surface layer thickness over which the reaction takes place, which is on the order of a lattice parameter, c α is the limiting value of the point defect concentration at the surface, ν α is the attempt frequency, Δg α is the surface kinetic barrier, f is the limiting value of the free energy at the surface, and κ is the local surface curvature.
It is clear that the sharp-interface formulation, Eqs. (3) - (7), is a generalization of the classical particle growth models, Eqs. (2.1)-(2.5), that is suitable for irradiation-driven void evolution. The canonical form of the surface boundary condition, Eq. (7), derived from transition state theory can be considered as a generalized interfacial constitutive law. By expanding the exponential term in Eq. (7) up to first order by assuming (μ α ± ΔE)/k B T < < 1 and using Eq. (6), one recovers the linear dynamical Gibbs-Thompson relation.
Hochrainer and El-Azab conducted numerical simulations of their model for different scenarios of void growth and shrinkage (El-Azab et al., 2014; Hochrainer & El-Azab, 2015). They demonstrated that the point defect reaction with the void surface has a strong effect on the overall kinetics of the problem. In particular, they showed that, as the surface barrier increases, the rate of void growth diminishes and the deviation of the point defect concentrations from their equilibrium values increases. This is of course consistent with a kinetic type interfacial constitutive law (such as Eq. ((2.5) or Eq. (7)) that accounts for kinetic drag. However, the diffusion-controlled models that assume local equilibrium at the surface cannot predict such effect. Therefore, the sharp-interface formulation of the problem of void growth overcomes the limitations of the rate theory, and hence provides a reference for constructing the corresponding diffuse-interface models of the problem.
Diffuse-interface description
Phase field (diffuse-interface) models of void growth were introduced (El-Azab et al., 2014; Hochrainer & El-Azab, 2015; Millet & Tonks, 2011a) with the goal to avoid the numerical difficulties encountered with sharp interface models. In accordance with the famous classification of Hohenberg and Halperin for the models of dynamical critical phenomena in condensed-matter physics (Hohenberg & Halperin, 1977), existing phase field models for voids are of type B or type C. In the phase field models of type B, the point defect concentrations are the only order parameters, and hence the kinetic evolution equations are given by generalized Cahn-Hilliard equations governing the local balance of point defect concentrations (Yu & Lu, 2005; Hu et al., 2009; Hu & Henager, 2009; Hu & Henager, 2010; Li et al., 2010; Semenov & Woo, 2012; Li et al., 2013; Xiao et al., 2013). In the phase field models of type C, an additional non-conserved order parameter is used to distinguish between voids and matrix (Rokkam et al., 2009; Millett et al., 2009; Millett et al., 2011c; Millett et al., 2011b). Earlier type B models for void growth considered vacancies only (Yu & Lu, 2005; Hu & Henager, 2009; Hu & Henager, 2010; Semenov & Woo, 2012; Xiao et al., 2013) but were later generalized to account for vacancies and interstitials (Hu et al., 2009; Li et al., 2010). In model B, the free energy of the system is expressed in the form
We do not specify the form of f(cv, ci) but require it to have two global minimizers that correspond to the two homogeneous phases, namely the void phase in which cv = 1 and ci = 0 and the solid phase with \( {c}_{\mathrm{v}}={c}_{\mathrm{v}}^{\mathrm{eq}} \) and \( {c}_{\mathrm{i}}={c}_{\mathrm{i}}^{\mathrm{eq}} \). ε and q are constants related to the surface energy and diffuse-interface width. The chemical potentials are obtained from (8) as
and the defect fluxes are expressed in terms of the chemical potentials by
The point defect balance equations are written as follows
where α = i, v. As mentioned earlier, for consistency with the sharp-interface description, we require P α = 0, M α = 0, and Riv = 0 when cv = 1 and ci = 0. The balance equations are to be solved in conjunction with the following initial and boundary conditions,
Here we take the void phase to be in Ω−, and consider a zero flux boundary condition on the outer boundary ∂Ω, which has a unit outward normal m. The results of the asymptotic analysis presented here also hold for Dirichlet or periodic boundary conditions. We further require Cv(x) and Ci(x) to be outside of the spinodal (unstable) regime. This means that Cv(x) and Ci(x) locally satisfy the condition that the Hessian matrix of the Helmholtz free energy density is positive-definite at any point in the solid phase. Note that Cv(x) and Ci(x) have also to be compatible with Eq. (12.3).
Model C version presented below is a generalization of phase field model C presented by Bi and Sekerka (Bi & Sekerka, 1998) that ensures thermodynamic consistency, i.e., non-negative entropy production. In the model, the total free energy of the heterogeneous system is given by a Ginzburg-Landau (Cahn-Hilliard) functional (Ginzburg & Landau, 1950; Cahn & Hilliard, 1958) as
In the above, η is the non-conserved order parameter such that η = ηM in the solid matrix and η = ηV in the void. g(η) is a regular double-well potential in the non-conserved order parameter. f(η, cv, ci) is the generalized Helmholtz free energy density (per unit volume) with two global minimizers, namely, (ηV, cv = 1, ci = 0) and \( \left({\eta}^{\mathrm{M}},{c}_{\mathrm{v}}={c}_{\mathrm{v}}^{\mathrm{eq}},{c}_{\mathrm{i}}={c}_{\mathrm{i}}^{\mathrm{eq}}\right) \) representing the void phase and the matrix phase, respectively. ω, ε, q, and h are constants that determine the surface energy and the diffuse-interface width. The non-classical (non-local) chemical potentials are the functional derivatives of the free energy functional with respect to defect concentrations:
We also define a generalized driving force that evolves the non-conserved order parameter as
The diffusive fluxes of the point defects are then expressed as
The point defect balance and the phase field evolution equations can then be written as
Here, α = i, v and L is the Ginzburg-Landau (Allen-Cahn) mobility (Allen & Cahn, 1979), which is considered to be independent of all order parameters for simplification. Absorption of defects at sinks can be easily added to this description by assuming smeared sinks or discrete sinks. Since the kinetic eqs. (17.1) and (17.2) guarantee a reduction of the free energy along the evolution path, the free energy functional (Eq. (13)) is a Lyapunov functional and the two global minimizers (homogeneous phases) are Lyapunov stable.
In current phase field models of void growth (El-Azab et al., 2014; Rokkam et al., 2009; Millett et al., 2009; Millett et al., 2011c; Millett et al., 2011b), and for consistency with the corresponding sharp-interface description, the term P α and the parameters M α and Riv are required to vanish in the void phase. This is usually accomplished by expressing them as functions of the order parameters such that the values of these functions are identically zero when the order parameters take on their equilibrium values that define the void phase. We will show here using the method of matched asymptotic expansions that this is indeed necessary for the phase field models to recover the desired sharp-interface limit.
Asymptotic limits of phase field models for void growth
Formal asymptotic analysis of phase field models for voids is crucial for proving thermodynamic and mathematical consistency of the models and for determining their parameters from sharp-interface counterpart. Several analyses of phase field models of type B and C for other phenomena exist in literature (Provatas & Elder, 2010; Pego, 1989; Dai & Du, 2012; Elder et al., 2001; Fife, 1992; Emmerich, 2008; Garcke et al., 2004; Ahmed et al., 2016; Caginalp, 1989). The ones relevant to our situation are those by Pego (1989), Dai and Du (2012) for the case of model B and Elder et al. (2001) for the case of model C. Pego (1989) has shown that the Cahn-Hilliard equation with a constant mobility in a binary system recovers the classical Stefan problem at fast time scale (growth stage) and the classical Mullins-Sekerka problem at slow time scale (coarsening stage). Dai and Du (2012) have generalized the analysis to the case of highly dissimilar mobility in a two-phase system, which is relevant to voids since the diffusive mobility is zero in the void phase. Our asymptotic analysis for phase field models B for void growth generalizes their work to driven ternary systems. For the case of model C, Elder et al. (2001) presented an asymptotic analysis for the low supersaturation case in a non-driven binary system.
We present here an elaborate formal asymptotic analysis of phase field models C for void growth for both cases of high and low driving forces in driven ternary systems. For the case of low driving force limits of models B and C, we discuss the connection between the coarsening limits and the classical diffusion- or interface-controlled Lifshitz-Slyozov-Wagner theories of Ostwald ripening (Lifshitz & Slyozov, 1961; Rahaman, 2003; Mullins & Sekerka, 1963; Niethammer, 2000; Dai et al., 2010).
Asymptotic analysis of model B
Model B is summarized in Eqs. (8) through (12). First, all fields are expanded in terms of the small parameter, ε, appearing in the free energy expression, which is proportional to the interface width, both in the bulk domains (outer expansion) and in the close vicinity of the interface (inner expansion). The two expansions are then matched at the interface in the limit ε → 0. From the matching, one recovers both the boundary conditions on the front required for solving the field equations in both phases, and the equation of motion of the interface (void surface). The geometry relevant to this analysis is shown in Fig. 1. The left part gives a schematic of an irradiated solid with a void ensemble. Part of a void surface is amplified to the right as a diffuse interface of a finite thickness that partitions the domain into an inner (diffuse interface) and outer regions. Values of the order parameters across the void surface are shown in Ω− (void) and matrix (Ω+). ∂Ω is the outer boundary and Γ is the sharp void surface.
A schematic of the diffuse interface and the local coordinate system used in the asymptotic analysis of model B
In the outer region, and for α = i, v, we expand the fields in the form
Note that superscripts on ε denote powers while superscripts on the field quantities such as c α , μ α , etc., denote the order in the perturbation expansion. Explicit expressions of \( {\mu}_{\alpha}^0 \), \( {\mathbf{J}}_{\alpha}^0 \), etc., can be obtained by expanding Eqs. (9) and (10) in ε and equating terms of the same order. Moreover, any function of c α must be expanded in Taylor series around \( {c}_{\alpha}^0 \). By doing so, we arrive at
We only need terms up to first order in our analysis. The different orders of the outer problem can now be obtained by plugging in the outer expansions in the balance equations. The leading order (ε0) outer equation is now given by
Similar equations hold for the higher order terms. Therefore, one recovers the regular balance equations in the relevant phases. It is to be noted here that requiring P α = 0, M α = 0, and Riv = 0 when cv = 1 and ci = 0, i.e., in the void phase, makes Eq. (20) trivial. Boundary conditions on the front necessary for solving these equations will be derived from the inner expansion.
At the physical interface, we define a local orthogonal coordinate system (r, s), where r is the normal distance from the point x in Ω to the interface Γ(t), such that r > 0 in Ω+ and r < 0 in Ω−, and s being the (d-1) coordinates perpendicular to r and tangent to Γ, where d is the dimensionality of the problem; see Fig. 1. One immediately obtains
Here, n is the unit normal to Γ pointing toward Ω+, κ is the mean curvature (sum of principal curvatures) of Γ, positive when the center of curvature lies within Ω−, and v is the normal velocity of the interface, positive when the interface moves toward Ω+(for a growing void). The sign convention here is the opposite of the sharp-interface model of Hochrainer and El-Azab (2014; 2015). Since the gradients of the fields are much higher near the interface, we further introduce a stretched variable, z = r/ε. Hence, in the moving coordinate system (z, s) the spatial and time derivatives transform as follow
In the above, ∇ s , ∇ s ⋅ and \( {\nabla}_s^2 \) are the (d-1) surface gradient, divergence and Laplacian operators, respectively.
In the inner region, we expand the fields as (tilde is used to distinguish inner expansion of all fields),
with similar expressions for other field quantities as in Eq. (18). The expressions for the different orders of \( {\tilde{\mu}}_{\alpha } \)are then given by
where σv = 1 and σi = q. Now matching conditions to connect the outer and inner expansions at the interface must be found. The derivation of such conditions was presented before (Provatas & Elder, 2010; Pego, 1989; Dai & Du, 2012; Elder et al., 2001; Fife, 1992; Emmerich, 2008; Garcke et al., 2004; Ahmed et al., 2016; Caginalp, 1989), and hence we only summarize the final relations here. Given any relevant field quantity ϕ, the following matching conditions hold:
By the substitution of the relations of Eq. (22) in Eq. (11), the inner equations can now be written as
These inner equations must be solved simultaneously. For the leading order (ε0) we have
Requiring bounded solutions and taking into account phase change across the interface, the only possible solution is
For α = i, v, Eqs. (28) are the Euler-Lagrange equations that give the equilibrium planar profiles for the point defect concentrations. Therefore, the interfacial region is in equilibrium and we have the classical equilibrium thermodynamic conditions on the front, namely (by using the matching condition Eq. (25.1)),
For the next-to-the leading order (ε), we obtain (recall that \( {\tilde{\mu}}_{\alpha}^0=0 \))
which, by a direct integration, yields
By employing the matching conditions Eq. (25.1) and Eq. (25.2) and considering the flux definition (Eq. (19.5)), we arrive at the Stefan (jump) conditions,
Now if one requires M α to be zero when cv = 1 and ci = 0, that is the defect mobility vanishes in the void phase, Eq. (32) reduces to
where \( {\left[{\mathbf{J}}_{\alpha}^0\cdot \widehat{\mathbf{n}}\right]}^{+} \)represents the normal flux of the corresponding specie calculated at the interface from the solid side.
According to Eqs. (19.1), (19.5), (20), (29), (33), phase field model B for void growth under irradiation reduces to a generalized one-sided classical Stefan problem. Another limiting case of interest for model B is the one for low driving forces. This limiting case represents the coarsening behavior of the system when the driving forces have been exhausted. As shown below, it is in this limit that one recovers a diffusion-controlled growth equation for the void radius as the one used in the rate theory models (Eq. (1.2)). To derive this limit, we consider all the driving forces to be of order ε. Specifically, we require
The last condition states that the initial supersaturation in the solid matrix is of order ε. The other initial and boundary conditions are the same as in Eq. (12).
In this limit, one should study the evolution of the system on the slow time scale, t1 = εt. The outer expansions are the same as in Eq. (18), with t replaced by t1. The relations of Eq. (19) still hold. The outer equations are then given by,
For the leading order (ε0), we have,
Eq. (36) indicates that the leading order terms in the bulk phases are in steady state. Strictly speaking, we need boundary conditions on the front to be able to solve Eq. (36). However, since by construction (recall the initial conditions) we require the bulk phases to be close to equilibrium, then
We will see from the inner expansion that the leading order terms on the front take on their equilibrium values, and hence the outer solution is in fact given by Eq. (37).
For the next-to-the-leading order (ε), we have
where \( {\mathbf{J}}_{\alpha}^1 \)is now given by \( {\mathbf{J}}_{\alpha}^1={M}_{\alpha}^0\nabla {\mu}_{\alpha}^1 \) since \( {\mu}_{\alpha}^0=0 \) (see Eq. (19)). Therefore, the first order terms in the bulk phases are in steady states. As before, boundary conditions on the front required for solving Eq. (38) will be deduced from the inner expansion via the matching conditions. Note that all the terms in Eq. (38) vanish in the void phase by construction, and hence Eq. (38) is trivially satisfied there. Therefore, the next order in the outer expansions (ε2) must be considered to determine \( {c}_{\alpha}^1 \) and \( {\mu}_{\alpha}^1 \) in the void phase. However, this is not necessary since the front velocity will be determined solely from the data in the solid side, similar to the previous limit (see Eq. (33), as will be shown below.
In the inner region, we expand the fields as in Eq. (23), with t replaced by t1 . The relations of Eq. (24) still hold. The relation in Eq. (22.3) becomes
where \( {v}_1={\partial}_{t_1}r \) is the front velocity measured with respect to the slow time scale. Hence the inner equations can now be written as,
For the leading order (ε0) we get
Again, this gives \( {\tilde{\mu}}_{\alpha}^0=0 \) for α = i, v as in Eq. (28), and hence we have the same classical equilibrium conditions on the front as in Eq. (29). Moreover, in the present case this indicates that equilibrium holds everywhere in the domain to the leading order, i.e.,
For the next-to-leading order (ε), we obtain (recall that\( {\tilde{\mu}}_{\alpha}^0=0 \))
For Eq. (43) to hold everywhere in the interfacial zone, we must have
where A α (s, t) is a function that does not depend on z. Note that this indicates that the chemical potentials are constants in the interfacial region (continuous in the sharp-interface sense). In order to determine this function, we substitute for \( {\tilde{\mu}}_{\alpha}^1 \) from Eq. (24.2) in Eq. (44), and then multiply Eq. (44) by \( {\partial}_z{\tilde{c}}_{\alpha}^0 \) and integrate in z from −∞ to +∞. In doing so, we perform integration by parts on some terms taking into account the fact that \( {\tilde{\mu}}_{\alpha}^0=0 \) for α = i, v and \( {\partial}_z\left[{\partial}_{{\tilde{c}}_{\alpha }}f\left({\tilde{c}}_{\mathrm{v}}^0,{\tilde{c}}_{\mathrm{i}}^0\right)\right]={\partial}_{{\tilde{c}}_{\alpha}}^2f\left({\tilde{c}}_{\mathrm{v}}^0,{\tilde{c}}_{\mathrm{i}}^0\right)\;{\partial}_z{\tilde{c}}_{\alpha}^0+{\partial}_{{\tilde{c}}_{\beta }{\tilde{c}}_{\alpha}}^2f\left({\tilde{\eta}}^0,{\tilde{c}}_{\mathrm{v}}^0,{\tilde{c}}_{\mathrm{i}}^0\right){\partial}_z{\tilde{c}}_{\beta}^0 \) for (α, β) = (v, i), (i, v), we arrive at
In the above, \( B\left(s,t\right)=\underset{-\infty }{\overset{+\infty }{\int }}\left({\tilde{c}}_{\mathrm{v}}^1{\partial}_z{\tilde{c}}_{\mathrm{i}}^0-{\tilde{c}}_{\mathrm{i}}^1{\partial}_z{\tilde{c}}_{\mathrm{v}}^0\right)\kern0.24em {\partial}_{{\tilde{c}}_{\mathrm{v}}{\tilde{c}}_{\mathrm{i}}}^2f\left({\tilde{c}}_{\mathrm{v}}^0,{\tilde{c}}_{\mathrm{i}}^0\right)\kern0.24em dz \), \( {\gamma}_{c_{\mathrm{i}}}=\underset{-\infty }{\overset{+\infty }{\int }} q\varepsilon {\left({\partial}_z{\tilde{c}}_{\mathrm{i}}^0\right)}^2 dz \), and\( {\gamma}_{c_{\mathrm{v}}}=\underset{-\infty }{\overset{+\infty }{\int }}\varepsilon {\left({\partial}_z{\tilde{c}}_{\mathrm{v}}^0\right)}^2 dz \) such that the total surface energy is given as \( \gamma ={\gamma}_{c_{\mathrm{v}}}+{\gamma}_{c_{\mathrm{i}}} \). Using the matching conditions (25.1) and (25.2) and eliminating B(s, t), Eq. (45) reduces to the more familiar form
which is the static (equilibrium) Gibbs-Thompson condition. Therefore, at order (ε) in the inner expansion, one deduces the interfacial constitutive law (boundary condition) required for solving the outer problem Eq. (38).
In order to derive a relation for the front velocity, one must consider the next order (ε2) in the inner expansion. At that order, the inner equations are given by (where we utilized the fact that \( {\tilde{\mu}}_{\alpha}^0={\partial}_z{\tilde{\mu}}_{\alpha}^1=0 \))
By direct integration, we deduce
which after using the matching conditions (Eq. (25.1) and Eq. (25.2)) becomes
Since by construction M α = 0 when cv = 1 and ci = 0, Eq. (49) reduces to
The front velocity is thus determined solely from the normal fluxes coming to the interface from the solid matrix. The analysis so far shows that, for the low driving force limit that is suitable for describing the coarsening stage, the phase field model B for void growth under irradiation reduces to a generalized one-sided classical Mullins-Sekerka problem (Eqs. (38), (46), and (50)).
In model B the interface is always in local equilibrium since the chemical potentials and concentrations at the interface take on their equilibrium thermodynamic values. This means that the attachment kinetics of point defects to the void surface is ignored and the void growth process is diffusion-controlled. In addition, the curvature correction term (Gibbs-Thompson condition) appears only as a first order correction term (Eq. (46)) and hence the curvature effect is absent from the leading order terms. While this approximation could be tolerated for large voids (R > 100 nm), it certainly breaks down for small voids(R < 10 nm). As we will show later, model C can obviate such shortcomings of model B. Moreover, if one assumes that Mi = ε−1Mv, the contribution of the interstitial to the normal velocity of the interface will be of order ε, and hence can be ignored as in the phase field models that consider vacancies only (Yu & Lu, 2005; Hu & Henager, 2009; Hu & Henager, 2010; Semenov & Woo, 2012; Xiao et al., 2013). This could be a reasonable approximation at low temperature where the interstitial are highly more mobile than vacancies (interstitial migration energy is usually much smaller than its vacancy counterpart in most solids), but it does not hold at high temperature.
Now we show how the coarsening limit of model B just discussed reduces to the famous rate theory approximation. For simplicity, we consider the case where interstitials can be ignored. Moreover, we assume constant vacancy mobility in the solid matrix which is assumed to be initially supersaturated with a uniform supersaturation of O(ε) (see Fig. 2). Furthermore, we assume the system is weakly driven such that the production and reaction terms are of order ε2 or lower. From Eqs. (38), (46), and (50), we then have
A schematic illustration of a spherical void growing in a solid matrix with a uniform low supersaturation (low driving force limit of model B)
In the above, we took \( 1-{c}_{\mathrm{v}}^{\mathrm{eq}}\approx 1 \). Hence, for weakly driven systems, we recover the regular Mullins-Sekerka problem as for the case of non-driven systems. For a spherical void of radius R embedded in an infinite solid matrix, see Fig. 2, Eqs. (51) reduce to:
The solution of Eq. (52.1) is simply given by
Eq. (52.2) thus yields
which is the well-known growth equation for diffusion-controlled processes (Olander, 1976; Was, 2017; Brailsford & Bullough, 1972; Krishan, 1982; Dubinko et al., 1989; Lifshitz & Slyozov, 1961; Rahaman, 2003). One can further define a critical radius at which the void does not grow or shrink (v = 0) by
This confirms that in the low driving force limit, model B recovers the rate theory description or, more precisely, the classical diffusion-controlled particle growth (Olander, 1976; Was, 2017; Brailsford & Bullough, 1972; Krishan, 1982; Dubinko et al., 1989; Lifshitz & Slyozov, 1961; Rahaman, 2003; Mullins & Sekerka, 1963; Niethammer, 2000; Dai et al., 2010). Xiao et al. (2013) have recently shown using numerical simulations that phase field model B predictions agree well with the rate theory predictions.
It is well-known that the classical Mullins-Sekerka problem exhibits kinetic scaling such that R(t) ∝ t1/3 (Lifshitz & Slyozov, 1961; Rahaman, 2003; Mullins & Sekerka, 1963; Niethammer, 2000; Dai et al., 2010). Lifshitz and Slyozov were the first to show that a collection of second phase particles follows such cubic growth kinetics (Lifshitz & Slyozov, 1961). The connection between the phase field model B (Cahn-Hilliard equation) and the diffusion-controlled LSW theory was discussed before by Bray (1994). However, only weakly driven systems can be considered to behave asymptotically as non-driven systems. In general, such kinetic scaling laws will not hold for driven systems since the supersaturation is not decaying with time as in non-driven systems.
Asymptotic analysis of model C
Eqs. (13) through (17) summarize the main equations of type C phase field models for void growth in irradiated solids (Rokkam et al., 2009; Millett et al., 2009; Millett et al., 2011c; Millett et al., 2011b), where a non-conserved order parameter η is considered to distinguish the void and matrix phases. We analyze here a generic form in which the Ginzburg-Landau (Cahn-Hilliard) free energy of the system is given by Eq. (13). Recall that we only require g(η) to have two global minimizers at ηVand ηM representing the void and matrix phases, respectively. Similarly, f(η, cv, ci) has two global minimizers that correspond to the void phase (ηV, cv = 1, and ci = 0) and the solid phase (ηM, \( {c}_{\mathrm{v}}={c}_{\mathrm{v}}^{\mathrm{eq}} \), and \( {c}_{\mathrm{i}}={c}_{\mathrm{i}}^{\mathrm{eq}} \)). The chemical potentials, the generalized driving force for the non-conserved order parameter, and the point defect diffusive fluxes are given by Eq. (14), Eq. (15), and Eq. (16), respectively. The point defect balance equations and the phase field evolution equation are given by Eqs. (17). We will use the following scaling for some of the parameters that appear in the above mentioned equations:
The reasoning behind this scaling is as follows. The first scaling ensures that the leading order terms of the non-conserved order parameters take on their equilibrium values in the bulk phases regardless of the supersaturation. The second guarantees that the curvature effect appear at leading order, in contrast to model B where it was a first order correction. The last relation can be explained if one defines a non-dimensionalized Peclet number as, \( pe=\frac{a_v}{D_{\alpha }} \), where D α is the defect diffusivity and a is the lattice parameter. Hence, Peclet number represents the ratio between the velocity of points on the surface and the speed of diffusion in the bulk. Clearly, the reaction of the point defect with the surface is relevant when Peclet number is small; otherwise the process is diffusion-controlled. In phase field model C, one can define Peclet number as \( pe=\frac{L{\varepsilon}^2}{D_{\alpha }} \). Now, if we assume pe = O(ε), we get L = O(ε−1) as in Eq. (56.3). It is worth noting that if one assumes that pe = O(1) or higher and hence L = O(ε−2) or higher, the asymptotic limits of model C reduce to their counterparts of model B.
The dynamical system of Eq. (17) becomes
where α = i, v, with the following initial and boundary conditions (see Fig. 3),
A schematic illustration of the local coordinate system used in the asymptotic analysis of model C showing values of the order parameters across the void (free) surface
Again, similar to the case of model B, we require the initial conditions in the solid matrix to be consistent with the fact that the solid matrix is metastable, i.e., to avoid the unstable (spinodal) region.
In the outer region we expand the fields as in Eq. (18). Note that u has now to start at O(ε−1) , i.e., u(x, t) = ε−1u−1(x, t) + u0(x, t) + εu1(x, t)(x, t) + ⋯. The analogue of Eq. (19) is now given by
The different orders of the outer problem can be obtained by substituting Eq. (59) in Eq. (57). For the leading order (ε−1) we have
For non-zero L−1 , the only possible solution of Eq. (60) that is compatible with the initial and boundary conditions (Eq. (58) is
For the next-to-the leading order (ε0) we have,
From Eq. (62.2), one has u0 = 0. Since away from the interface the gradients of the non-conserved order parameter should vanish (regardless of the supersaturation in the point defect concentrations), one must require \( {\partial}_{\eta }f\left({\eta}^0={\eta}^{\mathrm{M}},{c}_{\mathrm{v}}^0,{c}_{\mathrm{i}}^0\right)={\partial}_{\eta }f\left({\eta}^0={\eta}^{\mathrm{V}},{c}_{\mathrm{v}}^0,{c}_{\mathrm{i}}^0\right)=0 \), so that η1(x) = 0, for x ∈ Ω±. As usual, the boundary conditions on the front required to solve Eq. (62.1) will be imposed from the matching with the inner expansion.
In the inner region, we expand the fields as in Eq. (23). We can then find explicit expressions for the different orders of \( {\tilde{\mu}}_{\alpha } \) and \( \tilde{u} \) as follows,
Here σv = 1 and σi = q. Hence the inner equations can now be written as,
These inner equations must be solved simultaneously. For the leading order (ε−1), we have
This is the Euler-Lagrange equation that gives an equilibrium planar profile (curvature effect not appearing) for the non-conserved order parameter. Therefore, we have \( {\tilde{\eta}}^0\left(-\infty \right)={\eta}^{\mathrm{V}} \) and \( {\tilde{\eta}}^0\left(+\infty \right)={\eta}^{\mathrm{M}} \). For the next-to-the leading order (ε0), one has
In contrast to model B, the above coupled equations indicate that the leading order chemical potential (and hence point defect concentration) profiles in the interfacial zone are generally not the equilibrium planar profiles but rather are steady-state profiles consistent with a specific interface velocity. In order to demonstrate that, we first note that from Eq. (66.1) we deduce
Hence the leading order chemical potentials are constant in the interfacial region (continuous in the sharp-interface sense). We then multiply Eq. (66.2) by \( {\partial}_z{\tilde{\eta}}^0 \) and integrate inz from −∞ to +∞, namely,
Using integration by parts and observing Eq. (65), one can show that\( \underset{-\infty }{\overset{+\infty }{\int }}{\partial}_z{\tilde{\eta}}^0\left({\omega}^{-1}\;{d}_{\tilde{\eta}}^2g\left({\tilde{\eta}}^0\right)\;{\tilde{\eta}}^1\hbox{-} {h}^{-1}{\partial}_z^2{\tilde{\eta}}^1\right)\kern0.5em dz=0 \). Furthermore, to proceed with the integration on the right hand side, we recall that from the chain rule that we have
By plugging Eq. (69) into Eq. (68) and taking into account Eq. (67), Eq. (63.1), and the matching condition Eq. (25.1), we finally arrive at the dynamical Gibbs-Thompson equation:
In the above, \( {\gamma}_{\eta }=\underset{-\infty }{\overset{+\infty }{\int }} h\varepsilon {\left({\partial}_z{\tilde{\eta}}^0\right)}^2 dz \) is the surface energy and δ = h ε2/γ η is the diffuse interface width that is on the order of a surface layer thickness.
The interfacial balance (Stefan/jump) condition is obtained from the next order. At order (ε) we have (recall that \( {\partial}_t{\tilde{\eta}}^0={\partial}_z{\tilde{\mu}}^0=0 \)),
At this level, it is sufficient to consider only the first equation, i.e., the second equation gives a higher order correction to the dynamical Gibbs-Thompson condition (Eq. (70)). Direct integration of Eq. (71.1) gives
which after using the matching condition Eq. (25.1) becomes
Since according to the dynamical Gibbs-Thompson condition the point defect concentrations on the front may deviate from their equilibrium values, the Stefan condition can be restricted to be one-sided only if we require that M α = 0 at η = ηV, cv = 1, and ci = 0. In this case, Eq. (73) immediately reduces to
According to Eqs. (62), (70), (74), phase field model C for void growth under irradiation reduces to a generalized one-sided Stefan problem with kinetic drag.
Now we embark on the task of determining the low driving forces limit of model C as was done before for model B. Again, we assume all the driving forces are of O(ε) as in Eq. (35). Since the leading orders of the order parameters in the bulk phases are in equilibrium, using different scaling for the non-conserved order parameter as before is unnecessary and we only require
The chemical potentials are then still given by Eq. (14). However, the generalized force for evolving the non-conserved order parameter (Eq. (15) reduces to
As before, for such limit one should study the evolution of the system on the slow time scale, t1 = εt. The dynamical system becomes
The initial and boundary conditions are the same as in Eq. (58) with the exception that the supersaturation in the solid matrix is of O(ε), specifically, \( {c}_{\alpha}\left(x,0\right)={c}_{\alpha}^{\mathrm{eq}}+\varepsilon {C}_{\alpha}^1(x) \) for x ∈ Ω+. In the outer region we expand the fields as before. Note that u now starts normally at O(ε0), i.e., u(x, t) = u0(x, t) + εu1(x, t)(x, t) + ε2u2(x, t) + ⋯, with similar expression for other field quantities. The relations for the different orders of μ α (Eq. (59.1) and Eq. (59.2)) still hold but their counterparts for u are now given by
For the leading order (ε0) of the outer problem, we have
For both of the above equations to be satisfied simultaneously and given the initial and boundary conditions, we conclude that the leading order terms of the order parameters take on their equilibrium values in the bulk phases, namely,
For the next-to-the-leading order (ε), one has
From Eq. (81.2), one has u1 = 0. Again in order to ensure that η1(x) = 0, for x ∈ Ω±, we require (see Eq. (78.2))
For the point defect concentrations, Eq. (81.1) indicates that the first order terms in the bulk phases are in steady state. We will deduce the boundary conditions on the front required for solving Eq. (81.1) from the inner expansion via the matching conditions as before.
We expand the field in the inner region as before. The relations of Eq. (63.1) and Eq. (63.2) are still applicable. The corresponding relations for \( \tilde{u} \) are now given by
The inner equations can now be written as,
Hence, for the leading order (ε0), we have
Again, these equation give \( {\tilde{u}}^0=0 \) and \( {\tilde{\mu}}_{\alpha}^0=0 \). Therefore, equilibrium holds everywhere to leading order, specifically,
For the next-to-the leading order (ε), we obtain (recall that \( {\tilde{\mu}}_{\alpha}^0=0 \))
Again, the above coupled equations indicate that the first order chemical potentials (and hence point defect concentrations) are in steady-state. Their profiles in the interfacial zone are consistent with specific interface velocity as we will show below. From Eq. (87.1), we deduce
Here A α (s, t) is a constant of integration, which is to be determined. According to Eq. (88), the first order chemical potentials are continuous across the interface. In order to determine A α (s, t), we follow the same procedure as with model B, i.e., we multiply Eq. (88) by \( {\partial}_z{\tilde{c}}_{\alpha}^0 \) and integrate in z from −∞ to +∞, we get
In the above,
and
We apply the same procedure for Eq. (87.2) (in this case, we multiply by \( {\partial}_z{\tilde{\eta}}^0 \)) and we arrive at
From Eq. (89) and Eq. (90) and using the matching conditions (25.1) and (25.2), one can deduce
such that the total surface energy is given as\( \gamma ={\gamma}_{\eta }+{\gamma}_{c_{\mathrm{v}}}+{\gamma}_{c_{\mathrm{i}}} \). (Recall that \( {\gamma}_{\eta }=\underset{-\infty }{\overset{+\infty }{\int }} h\varepsilon {\left({\partial}_z{\tilde{\eta}}^0\right)}^2 dz \) in Eq. (70)). This is a milder version of the dynamical Gibbs-Thompson condition suitable for low driving forces since it indicates that the deviation from equilibrium at the interface is of O(ε) .
The Stefan (jump) conditions will be recovered from the next order in the expansion. For the next order (ε2), we have (recall that \( {\tilde{\mu}}^0={\partial}_z{\tilde{\mu}}^1=0 \))
It is sufficient here to consider only the first equation. By direct integration, we deduce
Which, after using the matching conditions (Eq. (25.1) and Eq. (25.3)), becomes
Again, since by construction M α = 0 at η = ηV, cv = 1, and ci = 0, Eq. (94) reduces to the one-sided Stefan (jump) condition
Hence in the coarsening stage phase field model C reduces to a generalized version of the one-sided Mullins-Sekerka problem with kinetic drag (Eqs. (81), (91), and (95)).
This is the first time where a scaling limit of model C is shown to recover the Mullins-Sekerka model with a kinetic drag. Note that this limit is different from the one obtained by Elder et al. (2001). In their asymptotic analysis, they obtained a time-dependent diffusion equation for the first order term of the concentration field in the bulk phases, i.e., a modified Stefan problem for the first order term. This is however inconsistent with the assumption of low driving force where one expects the diffusion to be quasistatic. The reason that they did not arrive at the same limit obtained here is the fact that they did not study the interface motion at the slow time scale suitable for the low driving force limit. Our derivations of the asymptotic limits of model C for high and low driving forces are consistent with and can be viewed as generalizations of their model B counterparts.
Based on the limits of model C derived above, it is clear that model C is able to incorporate the surface attachment kinetics into the overall growth kinetics, as opposed to model B which can only describe diffusion-controlled growth. Model C accomplishes this through the extra Allen-Cahn (time-dependent Ginzburg-Landau) equation that acts as an interfacial constitutive law ensuring positive interfacial entropy production associated with the interface motion. Hence, it is obvious that the Allen-Cahn mobility in the diffuse-interface formalism is the counterpart of the interface kinetic coefficient in the sharp-interface formalism. In particular, from Eq. (70) and (2.5), we have
However, a general nonlinear (canonical) form of the interfacial constitutive law as in Hochrainer and El-Azab model (Hochrainer & El-Azab, 2015) can only be obtained if one assumes a similar nonlinear form for the Allen-Cahn equation.
Similar to model B, for the case of low driving forces, one can deduce a compact expression for the growth of a void of radius R embedded in an infinite solid matrix. We use the same assumptions made for model B (see Fig. 4). The problem then reduces to a regular Mullins-Sekerka model with kinetic drag (Mullins & Sekerka, 1963; Niethammer, 2000; Dai et al., 2010); specifically,
A schematic of a spherical void growing in a solid matrix with a uniform supersaturation for the case of model C. The order parameter η takes on its equilibrium values ηM and ηV in the matrix and void, respectively
Note that for the solidification case, the solution of the Mullins-Sekerka model with kinetic drag was presented before by Niethammer (2000). It is straightforward to show that the solution of Eq. (97) is
In the above, Leff is an effective mobility. Note that, from Eq. (98.1) or from Eq. (97.2) and Eq. (98.2), the vacancy chemical potential at the interface is given by
The difference between the coarsening limits of model B and C, e.g., the regular Mullins-Sekerka problem and the Mullins-Sekerka with kinetic drag problem is obvious. By comparing Eq. (98.2) and Eq. (54), we deduce that model C predicts lower coarsening rates for growing/shrinking voids than model B (since from Eq. (98.3), we always have \( {L}^{\mathrm{eff}}<{M}_{\mathrm{v}}^0 \)) . Also from Eq. (99) and Eq. (52.2), the vacancy chemical potential (and hence concentration) at the void surface obtained from model C is higher than the value predicted from model B (from the equilibrium Gibbs-Thompson condition) for a growing void. The situation is reversed for a shrinking void; both limits give of course the same result for a stationary void. Moreover, according to Eq. (99), for the case of a growing void the vacancy chemical potential (and hence concentration) at the void surface increases as the Allen-Cahn mobility decreases (or equivalently, as the interface kinetic coefficient decreases or the surface barrier increases). All these predictions of model C are in good agreement with the recent simulation results obtained by Hochrainer and El-Azab from their sharp-interface model (El-Azab et al., 2014; Hochrainer & El-Azab, 2015).
One can also conclude from Eq. (98) and Eq. (99) that the void growth equation obtained from model C is a general growth equation, from which two limiting cases can immediately be identified. When \( R\;\delta L>>{M}_{\mathrm{v}}^0 \), we recover the diffusion-controlled growth given by
Therefore, for the diffusion-controlled growth, the interface is in local equilibrium, as anticipated. On the other hand, when \( {M}_{\mathrm{v}}^0>>R\;\delta L \), we recover the attachment-controlled growth given by
Hence, as one should expect, for infinitely fast bulk diffusion kinetics, the vacancy chemical potential/concentration at the void surface is the same as in the bulk.
It is clear from Eq. (98.2) that a transition from interface-controlled kinetics to diffusion-controlled kinetics may always take place. For a particular material, the thermodynamic and kinetic parameters are given and the only factor that determines the prevailing kinetics is the particle radius. For small particles \( \left(R<<{M}_{\mathrm{v}}^0/L\;\delta \right) \), interface-controlled kinetics dominates. For large particles \( \left(R>>{M}_{\mathrm{v}}^0/L\;\delta \right) \) diffusion-controlled kinetics dominates. Hence, while a small particle grows, its kinetics changes from interface- to diffusion-controlled. Therefore, in the mean-field limit a collection of second phase particles in weakly- or non-driven systems coarsen such that the average particle size follows parabolic growth (R(t) ∝ t1/2) at short times and cubic growth (R(t) ∝ t1/3) at long times. Such behavior was recently captured by the numerical simulations of Dai et al. (2010). Therefore, the mean-field limit of the Mullins-Sekerka model with kinetic drag gives rise to a generalized Lifshitz-Slyozov-Wagner (LSW) theory from which the diffusion- and interface-controlled theories emerge as limiting cases. Moreover, Kockelkoren and Chaté reported similar trend from their numerical simulations of the late stage of coarsening in phase field model C (Kockelkoren & Chaté, 2002). Their results are then in agreement with our asymptotic analysis presented here. Nevertheless, as we mentioned before, such kinetic scaling laws are not expected to hold in driven systems.
Summary and conclusions
The diffuse-interface void growth models of type B and C were analyzed by deriving their sharp-interface limits. The limits obtained from the asymptotic analyses are summarized below followed by concluding remarks. A sample simulation is included in the Appendix.
Sharp-interface limits of phase field model B
High driving force limit (growth stage)
In the high driving force limit, suitable for describing the growth stage, all the driving forces are assumed to be of O(1) and the system evolves on the regular time scale, t = O(1). In this limit, the phase field model B reduces to a generalized one-sided classical Stefan problem consisting of:
-
Bulk balance laws and constitutive laws:
-
Interfacial balance laws and constitutive laws:
The quantity \( {\left[{\mathbf{J}}_{\alpha}^0\cdot \mathbf{n}\right]}^{+} \)is the normal flux of the defect species at the solid side of the void surface.
Low driving force limit (coarsening stage)
In the low driving force limit, suitable for describing the coarsening stage, all the driving forces are assumed to be of O(ε) and the system evolve on the slow time scale, t = O(ε−1). In this limit, phase field model B reduces to a generalized one-sided classical Mullins-Sekerka problem. In this limit, the leading order terms of the point defect concentrations and chemical potentials take on their equilibrium values and we solve for a first order correction. The governing equations are:
-
Bulk balance laws and constitutive laws:
-
Interfacial balance laws and constitutive laws:
Again \( {\left[{\mathbf{J}}_{\alpha}^0\cdot \mathbf{n}\right]}^{+} \)is the normal flux of the defect species at the solid side of the void surface.
Sharp-interface limits of phase field model C
High driving force limit (growth stage)
In the high driving force limit, suitable for describing the growth stage, all the driving forces are assumed to be of O(1) and the system evolve on the regular time scale t = O(1). In this limit, phase field model C reduces to a generalized one-sided Stefan problem with kinetic drag. The governing equations are:
-
Bulk balance laws and constitutive laws:
where the leading order terms of the non-conserved order parameters take on their equilibrium values in the bulk phases, namely,
-
Interfacial balance laws and constitutive laws:
Low driving force limit (coarsening stage)
In the low driving force limit, suitable for describing the coarsening stage, we assume that all the driving forces are of O(ε) and study the evolution of the system on the slow time scale, t = O(ε−1). In this limit, phase field model C reduces to a generalized one-sided Mullins-Sekerka problem with kinetic drag. Hence, the leading order terms of the point defect concentrations and chemical potentials take on their equilibrium values and we solve for the first order correction. The governing equations are:
-
Bulk balance laws and constitutive laws:
-
Interfacial balance laws and constitutive laws:
According to the asymptotic analyses conducted here, a phase field model of type B recovers the diffusion-controlled growth and coarsening. This model thus ignores the surface attachment kinetics. On the other hand, a phase field model of type C is able to account for the surface attachment kinetics via the extra Allen-Cahn equation, which ensures non-negative interfacial entropy production associated with the interface motion. Therefore, for constructing models of voids/bubbles growth and coarsening, phase field models of type C must be used.
It is obvious that sharp- and diffuse-interface models are able to relax all the assumptions of the rate theory. In particular, they account for the heterogeneity of the void microstructure and consider the effect of surface attachment of point defects to the void surface on the overall kinetics. One can, however, to some extent account for the latter in the rate theory models by replacing the diffusion-controlled growth equation in the rate theory (Eq. (1.2)) by a general growth equation as Eq. (98.2). One also should replace the equilibrium Gibbs-Thompson condition (Eq. 1.3), by the dynamical Gibbs-Thompson condition which leads to Eq. (99). Nevertheless, this is only valid for the low driving force case where the assumption of quasistatic diffusion in the solid matrix is valid.
In passing, we recall here the recommendation by Hochrainer and El-Azab that a proper treatment of defect attachement to void surfaces, which is rigorously treated in their development of the sharp interface model for void evolution (Hochrainer & El-Azab, 2015), must be incorporated into phase field model construction to gurantee an accurate account of defect removal by reactions at the void surface. In earlier models of type C (El-Azab et al., 2014), such reactions were treated indirectly by amplifying the vacancy-interstitial recombination term, Riv, in the diffuse interface region. The form of this term, however, does not impact the results presented here.
While the analysis presented here was tailored for void growth and coarsening, it can be easily generalized to particle growth and coarsening by mass or/and heat diffusion in multi-component and multi-phase driven systems. Moreover, thanks to the varitional formulation of phase field models, one can account for long-range interactions such as elastic or electrostatic in a systematic fashion, as was done several times in literature (Provatas & Elder, 2010). The kinetic trends based on the analysis presented here will not change. One will only have to work with a generalized chemical potential that account for all the driving forces of short-range or long-range instead of the usual chemical potential used here. Therefore, our analysis represents a clear way of constructing diffuse-interface models and deducing their sharp-interface limits for investigating the growth and coarsening stages of non-equilibrium phase transitions in a general driven system.
References
K Ahmed, T Allen, A El-Azab, J. Mater. Sci. 51, 1261 (2016)
K Ahmed, M Tonks, Y Zhang, B Biner, A El-Azab, Comput. Mater. Sci. 134, 25 (2017)
SM Allen, JW Cahn, Acta Metall. 27, 1085 (1979)
L Amirouche, M Plapp, Acta Mater. 57, 237 (2009)
Z Bi, R Sekerka, Physica A 261, 95 (1998)
K Binder, Rep. Prog. Phys. 50, 783 (1987)
A Brailsford, R Bullough, J. Nucl. Mater. 44, 121 (1972)
AJ Bray, Adv. Phys. 43, 357 (1994)
G Caginalp, Phys. Rev. A 39, 5887 (1989)
JW Cahn, Acta Metall. 9, 795 (1961)
JW Cahn, JE Hilliard, J. Chem. Phys. 28, 258 (1958)
S Dai, Q Du, SIAM Appl. Math. 72, 1818 (2012)
S Dai, B Niethammer, RL Pego, Proc. Roy. Soc. Edinburgh A140, 553 (2010)
S De Groot, P Mazur, Non-Equilibrium Thermodynamics (North-Holland Publishing Company, Amsterdam, 1962)
VI Dubinko, PN Ostapchuk, VV Slezov, J. Nucl. Mater. 161, 239 (1989)
A El-Azab, K Ahmed, S Rokkam, T Hochrainer, Curr. Opin. Solid State Mater. Sci. 18, 90 (2014)
K Elder, M Grant, N Provatas, JM Kosterlitz, Phys. Rev. E 64, 021604 (2001)
H Emmerich, Adv. Phys. 57, 1 (2008)
PC Fife, Dynamics of Internal Layers and Diffuse Interfaces (SIAM, Philadelphia, 1992)
H Garcke, B Nestler, B Stinner, SIAM Appl. Math. 64, 775 (2004)
V Ginzburg, L Landau, Sov. Phys. JETP 20, 1064 (1950)
T Hochrainer, A El-Azab, Phil. Mag. 95, 948 (2015)
PC Hohenberg, BI Halperin, Rev. Mod. Phys. 49, 435 (1977)
S Hu, C Henager, H Heinisch, M Stan, M Baskes, S Valone, J. Nucl. Mater. 392, 292 (2009)
S Hu, H Henager, J. Nucl. Mater. 394, 155 (2009)
S Hu, H Henager, Acta Mater. 58, 3230 (2010)
J Katz, H Wiedersich, J. Chem. Phys. 55, 1414 (1971)
J Kockelkoren, H Chaté, Physica D 168-169, 80 (2002)
K Krishan, Radiat. Eff. 66, 121 (1982)
JS Langer, Rev. Mod. Phys. 52, 1 (1980)
Y Li, S Hu, R Montgomery, F Gao, X Sun, Nucl. Instr. Meth. in Phys. Res B 303, 62 (2013)
Y Li, S Hu, X Sun, F Gao, H Henager, M Khaleel, J. Nucl. Mater. 407, 119 (2010)
I Lifshitz, V Slyozov, J. Phys. Chem. Solids 19, 35 (1961)
R Mayer, L Brown, J. Nucl. Mater. 95, 58 (1980)
P Millet, M Tonks, Curr. Opin. Solid State Mater. Sci. 15, 125 (2011a)
P Millett, A El-Azab, S Rokkam, M Tonks, D Wolf, Comput. Mater. Sci. 50, 949 (2011b)
P Millett, A El-Azab, D Wolf, Comput. Mater. Sci. 50, 960 (2011c)
P Millett, S Rokkam, A El-Azab, M Tonks, D Wolf, Mod. Sim. Mater. Sci. Eng. 17, 064003 (2009)
W Mullins, R Sekerka, J. Appl. Phys. 34, 323 (1963)
B Niethammer, Proc. R. Soc. Edinb. A130, 1337 (2000)
D Olander, Fundamental Aspects of Nuclear Reactor Fuel Element (Technical Information Service, Springfield, VA, 1976)
RL Pego, Proc. Roy. Soc. London. Ser. A. 422, 261 (1989)
N Provatas, K Elder, Phase-Field Methods in Materials Science and Engineering (Wily-VCH, Weinheim, Germany, 2010)
M Rahaman, Ceramic Processing and Sintering (Marcel Dekker, New York, 2003)
S Rokkam, A El-Azab, P Millett, D Wolf, Mod. Sim. Mater. Sci. Eng. 17, 064002 (2009)
K Russell, Acta Metall. 19, 753 (1971)
A Semenov, C Woo, Acta Mater. 60, 6112 (2012)
G Was, Fundamentals of Radiation Materials Science-Metals and Alloys (Springer, Berlin, 2017)
Z Xiao, A Semenov, C Woo, SQ Shi, J. Nucl. Mater. 439, 25 (2013)
H Yu, W Lu, Acta Mater. 53, 1799 (2005)
Acknowledgements
The authors wish to thank Srujan Rokkam and Thomas Hochrainer for useful discussions.
Funding
This material is based upon work supported as part of the Center for Materials Science of Nuclear Fuel, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Sciences, Office of Basic Energy Sciences under award number FWP 1356, through subcontract number 00122223 at Purdue University.
Availability of data and materials
Not applicable.
Author’s contributions
Both authors of this paper have contributed equally to the research presented in this paper, as well as to the writing of the paper. Both authors read and approved the final manuscript.
Author information
Authors and Affiliations
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendix
Appendix
Sample void growth simulation
For simplicity and direct comparison with analytical solutions, we ignore interstitials and consider only vacancies. The main goal here is to demonstrate that model C captures both the diffusion-controlled (DC) and interface-controlled (IC) kinetics. This will be accomplished here by investigating the kinetics of void growth from a supersaturated matrix.
The simplest free energy of Model C that satisfies the requirements obtained from the asymptotic analysis has the form,
Here, c is a normalized vacancy concentration such that it equals 1 in the void and 0 in the matrix, e.g., \( c=\left({c}_{\mathrm{v}}-{c}_{\mathrm{v}}^{\mathrm{eq}}\right)/\left(1-{c}_{\mathrm{v}}^{\mathrm{eq}}\right) \). Therefore, in the void phase, the order parameters take on the values c = 1, η = 1, and in the matrix phase, the order parameters take on the values c = 0, η = 0. Note that the second term in the free energy functional represents a parabolic approximation of the thermodynamic free energy of the phases. This specific form of the free energy functional has the advantage of decoupling the bulk thermodynamic free energy from the interfacial energy as was shown by Amirouche and Plapp (2009). Specifically, the interfacial energy and interface width are given by
The parabolic free energy, that approximates the bulk free energy, is parameterized such that the resulting chemical potential in the matrix equals its counterpart from the ideal solution free energy (Eq. 5), e.g.,
The phase field kinetic coefficients are calculated as follows. The Cahn-Hilliard (chemical) mobility of vacancies is given by
where Dv is the vacancy diffusion coefficient. This guarantees that the Cahn-Hilliard equation recovers Fick’s law in the matrix. The Allen-Cahn mobility can now, thanks to the asymptotic analysis presented here, be determined as function of the surface kinetic barrier that appears in the corresponding sharp-interface model (Eq. 7), namely,
This result was obtained from comparing Eq. (91) with a linearized version of Eq. (7).
Now that all the phase field model parameters have been identified in terms of their sharp-interface parameters, quantitative simulations of void growth in irradiated solids can be conducted. Here, we chose copper as our material system as in the sharp-interface analysis presented in (Hochrainer & El-Azab, 2015). The phase field kinetic equations are solved here using a fully-coupled, fully-implicit scheme implemented in the finite-element framework MOOSE. The discretization and integration schemes and the nonlinear solver options are the same as in the work presented by Ahmed et al. (2017).
In this investigation, we conduct 2D simulations of void growth in supersaturated copper at 800 K. The formation and migration energies of vacancies in copper are 1 eV and 0.8 eV, respectively (Hochrainer & El-Azab, 2015). The attempt frequency and atomic volume are taken after (Hochrainer & El-Azab, 2015) as υ = 72300ns−1 and Ω = 0.012 nm3. At 800 K, the equilibrium vacancy concentration is 5 × 10−7and the vacancy diffusion coefficient is 22.5 × 10−3 nm2/ns. In order to compare the numerical results with the analytical predictions for infinite systems presented above, a domain size much larger than the initial void radius is used in conjunction with Dirichlet boundary conditions at the outer boundary. Specifically, we use a square domain of 512 nm × 512 nm and an initial void radius of 10 nm. A diffuse-interface width of 2 nm was assumed. The copper matrix is slightly supersaturated with vacancies, e.g., c = 10−5.
A few simulations were conducted to study the different kinetics of void growth in copper. The surface kinetic barrier for the vacancies was assumed to be multiple of their bulk migration energy as in the sharp-interface study (Hochrainer & El-Azab, 2015). The effect of surface kinetic barrier on the overall kinetics of void growth is captured in Fig. A1, Fig. A2, Fig. A3. Figure A1 shows snapshots of diffusion-controlled void growth where the surface kinetic barrier (Δg) was the same as the bulk migration energy (Em). Figure A2 and Fig. A3 demonstrate quantitatively the effect of the surface kinetic barrier on the growth process. As is clear from Fig. A2, for the case of diffusion-controlled, high gradients of the vacancy chemical potential are established across the domain with steepest gradients at the void-matrix interface. The value of the vacancy chemical potential inside the void and at the interface is set by the void radius and surface energy. On the other hand, for the case of interface-controlled kinetics (with Δg = 2.1Em), the chemical potential inside the void and at the surface is almost equal to the chemical potential in the bulk and hence no gradients develop in this case. All these results are in good agreement with the analytical predictions of the asymptotic analysis (Eqs. 98-101). Moreover, as evident from Fig. A3, the overall growth rate is strongly dependent on the value of the surface kinetic barrier. Higher values of the surface kinetic barrier lead to lower growth rates in agreement with the analytical predictions derived here and the sharp-interface results presented in (Hochrainer & El-Azab, 2015).
ᅟ
ᅟ
ᅟ
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Ahmed, K., El-Azab, A. An analysis of two classes of phase field models for void growth and coarsening in irradiated crystalline solids. Mater Theory 2, 1 (2018). https://doi.org/10.1186/s41313-017-0008-y
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s41313-017-0008-y