General Principles and Applications of PF Modeling

The first step in PF modeling lies in choosing and defining the fields of interest. These continuous variables are functions of space and time, and they are in most cases scalar fields, such as temperature, or the concentration of some chemical species of interest. In systems with solid-liquid interfaces, a phenomenological field variable is introduced in such a way that it varies continuously from 0 to 1 as one goes from a fully solid to a fully liquid phase. Multidimensional fields can be used as well, for instance, to describe the local composition of a multicomponent alloy, the local degree of chemical order, or the local crystallographic orientation of grains. These multidimensional fields may transform like vectors under symmetry operations, thus leading to a vectorial representation of the system and ten — sorial expressions for mobilities (as will be discussed later), but there are cases for which the multidimen­sional fields cannot be reduced to vectors.27 In all cases, an averaging procedure is necessary to define continuous field variables for systems that are intrin­sically discrete at the atomic scale. Various averages can be used, including (1) a spatial average over representative volume elements, which will corre­spond to the cells used for evolving the PF variables; (2) a spatial and temporal average; or (3) a spatial and ensemble average. The spatial averaging method is used most often, although in many cases the exact conditions of the averaging procedure are not defined. Section 1.15.3 will cover a model where this coarse-graining is performed explicitly and rigor­ously. The last two averaging procedures are rarely explicitly invoked, although one oftheir advantages is that a smaller volume can be used for the spatial average, thanks to the additional averaging per­formed either in time or in the configuration space of a system ensemble.

Turning now to the kinetic equations used to describe the evolution of these field variables, an important distinction is whether the field variable is conserved or nonconserved. For the sake of simplic­ity, the following discussion focuses on alloy systems. Let us consider two simple examples, one where the field variable is the local composition, C(r, t), in a binary A-B alloy system, and a second example, also for a binary alloy system, but this time with a fixed composition and where chemical ordering takes place. The degree of chemical order is described by the field S(r, t). For the sake of convenience, one

Подпись: d C (r,t) dt Подпись:Подпись: MVПодпись: [2]Подпись: [3]Подпись: F{C(r, t)}Подпись: [4]Подпись: dVПодпись:Подпись: amay normalize that field such that S(r, t) = 0 corre­sponds to a fully disordered state and S(r, t) = ±1 to a fully ordered state. The first field variable C(r, t) is globally conserved — assuming here that the system of interest is not exchanging matter with its environ­ment. This imposes the constraint that the time evolu­tion of the field variable at r is balanced by the divergence of the flux of species exchanged between the representative volume centered on r and the remainder of the system:

Чг =-VJ (V) [1]

One then makes use of linear response theory in the context of thermodynamics of irreversible processes28 to linearly relate the flux J(r, t) to the driving force responsible for this flux. Here, this driving force is the gradient of the chemical potential m(r, t) = 8F/8C(r, t), where F is the free energy of the system for a compositional field given by C(r, t). The resulting evolution equation is thus

‘ 8 F 8C(r, t)

where M is a mobility coefficient. In contrast, for the nonconserved order parameter S(r, t), its evolution is directly related to the free energy change as S(r, t) varies, so that by making use of linear response theory again

d S (r, t) r 8 F

dt = ~L 8S(r, t)

where L is the mobility coefficient for the nonconserved field S(r, t). Two important consequences of eqns [2] and [3] are worth noticing. First, although all extrema of the system free energy (i. e., minima, maxima, saddle points) are stationary states, often in practice only the minima can be obtained at steady state due to numerical errors. Second, the stationary state reached from some initial state may not correspond to the absolute minimum of the free energy. In order to overcome this problem, noise can be added to transform these deterministic equations into stochastic (Langevin) equations, as will be discussed in Section 1.15.3.

Following the work of Cahn and coworkers2-4 and Landau and Ginzburg,6 the free energy F is decomposed into a homogeneous contribution and an heterogeneous contribution. Treating the inhomo­geneity contribution as a perturbation of a homoge­neous state, one finds that, in the limit of small amplitude and long wavelength for this perturbation, the lowest order correction to the homogeneous free
energy is proportional to the square of the gradient of the field variable. For instance, returning to the sim­ple example of an alloy described by the concentra­tion field C(r, t), the total free energy can be written as dV [f (C) + k(VC)2] where f (C) is the free-energy density of a homoge­neous alloy for the composition C, and к the gradient energy coefficient, which is positive for an alloy system with a positive heat of mixing. A similar expression can be used in the case of a nonconserved order parameter, for example, S(r, t), or more generally, in the case of an alloy described by nC conserved order parameters and nS nonconserved order parameters f (Cb C2 ••• CnC ; Sb S2 • •• SnS )

nS

+ Kp (VCp) + Vij ViSqVjSq

p=1 q=1

An implicit summation over the indices i and j is assumed in the last term of eqn [5]. The number of nonzero and independent gradient energy coefficients V j for the nonconserved order parameters is dictated by the symmetry of the ordered phase. Specific exam­ples, for instance for the L12 ordered structure, can be found in Braun eta/.27 and Wang eta/.29 The free energy can also be augmented to include other contributions, in particular those coming from elastic fields using the elasticity theory of multiphase coherent solids pio­neered by Khachaturyan,30 in the homogeneous mod­ulus case approximation. This makes it possible to take into account the effect of coherent strains imposed by phase transformations or by a second phase, for exam­ple, a substrate onto which a thin film is deposited.31

Two important interfacial quantities, the excess interface free energy and the interface width, can be derived from eqn [5] for a system at equilibrium. We follow here the derivation given by Cahn and Hilliard.2 Considering the case of a binary alloy where two phases may form, referred to as a and p, and with respective B atom concentrations Ca and Cp, the existence of an interface between these two phases results in an excess free energy a dv[f(c) + k(vc)2 — CmB — (1 — C)mA] [6] where mA and mB are the chemical potential of A and B species when the two phases a and p coexist at

Подпись: dV K(V С)2 V Подпись:Подпись: sПодпись: Cp 2 dC KVC ca Подпись:Подпись: [12]Подпись: 2С - 1 Cap Подпись:Подпись:equilibrium. At equilibrium, this excess free energy is minimum. A homogeneous free energy Dfreferenced to the equilibrium mixture of a — and p-phases is introduced as

Df(c)=f(c)-[CmB + (1 — c)mA]

= сMC) — mB] + (1 — c)[ma(c) — mA] [7]

(Note that the ‘D’ symbol in Df in eqn [7] does not refer to a Laplacian.) The variational derivative of this excess energy with respect to the concentration field is given by

8s dDf дк 2 r і

8c = 1>C — 2kDC — дс(гс) [8]

At equilibrium, the excess free energy s is minimum, and the concentration field must be such that 8s/8C = 0. Thus,

= 2kDC + (vc)2 = (k(VC)2) [9]

дС 3l ; дС

Equation [9] must hold locally for any value of the concentration field along the equilibrium profile join­ing the a — and p-phases, and this can only be satisfied if

k(VC)2 = Df [10]

for all values of C(r). It is interesting to note that eqn [10] means that the equilibrium concentration profile is such that, at any point on this profile, the homoge­neous and inhomogeneous contributions to the total free energy are equal. The interfacial excess free energy is thus given by 2 dV Df = 2

V

This last integral over the spatial coordinates can be rewritten as an integral over the concentration field. Assuming a one-dimensional (1D) system for simplicity,

Cp

2 dCy/KAf

Ca

In order to proceed further, it is necessary to assume a functional shape for the concentration profile or for the homogeneous free energy Df. Expanding the free energy near the critical point Tc yields a symmetric double-well potential for the homogeneous free energy,2 which we write here as 2С — 1

Cap
with Cap = Cp — Ca = 1 — 2Ca = 2Cp — 1. Using eqn [10], the equilibrium concentration profile can now be obtained:

ОД="Т +1 [14|

Integration along this equilibrium profile from eqn [11] yields the interfacial energy

s 3Cap у8 KA/max [15]

Furthermore, the width of the equilibrium profile we, which is defined as the length scale entering the argu­ment of the hyperbolic tangent function in eqn [14], is given by

w=crvлЕ [16]

In this conventional approach to PFMs, Dfmax and к are phenomenological coefficients. Equations [15] and [16] play an important role in assigning values to these coefficients for a specific alloy system. The excess interfacial energy may be known experimen­tally or it may be calculated separately, for instance by ab initio calculations.32 If the interfacial width we is also known, one can obtain Dfmax and к from an inverse solution of eqns [15] and [16] (note that Cap is given by the equilibrium phase diagram). Even if we is not known, values for Dfmax and к can be chosen to yield a prescribed value for s. In all cases, it is important to recognize that any microstructural fea­ture that develops during the simulations is expressed in units of we. At elevated temperatures, as T! Tc, s vanishes2 while we goes to infinity, and therefore, at high enough temperatures, interfaces are diffuse, thus meeting this essential requirement underlying the PF method.

The PF eqns [2] and [3] are usually solved numer­ically on a uniform mesh with an explicit time inte­gration, using periodic boundary conditions, when surface effects are not of interest. When the free energy contains an elastic energy contribution, it is quite advantageous to use semi-implicit Fourier — spectral algorithms (see Chen9 and Feng et a/.33 for details). Variable meshing can also be employed, in particular to better resolve interfaces when they tend to be sharp, for instance at low temperatures.

A few examples selected from the literature serve to illustrate the capacity of PFMs to successfully reproduce a wide range of phenomena. In particular, Khachaturyan30 and his collaborators34-38 proposed a microelasticity theory of multiphase coherent solids,
which has been widely used to include a strain energy in the overall free energy. A method for systems with strong elastic heterogeneity has been proposed by Hu and Chen,39 which includes higher order terms that are usually neglected in Khachaturyan’s approach. Figure 1 illustrates the anisotropic morphology of Al2Cu precipitates growing in an Al-rich matrix.3 Bulk-free energies were calculated using a mixed- space cluster expansion technique, with input from first-principle calculations for about 40 different ordered structures with full atomic relaxations. Interfacial energies were calculated at T = 0 K from first-principle calculations as well, using configura­tions where the Al-rich solid solution and the tetrag­onal 0′-Al2Cu coexist. For the elastic strain energy calculations, the elastic constants of 0′-Al2Cu were calculated ab initio. An important feature of this sys­tem is that both elastic and interfacial energies are strongly anisotropic, and the PF approach makes it possible to include these anisotropies. Furthermore, when the high-aspect-ratio 0′-phase forms, its growth kinetics will be anisotropic as well, which can be included in a phenomenological way by introducing a dependence of the mobility on the orientation of the precipitate-matrix interface. Figure 1 illustrates that these three anisotropies, interfacial, elastic, and kinetic, are required to reproduce the morphology of 0′ precipitates.

Figure 2 illustrates another effect of coherency stress on microstructural evolution, this time for an A1-L10 order-disorder transition in a Co-Pt alloy.40 The tetragonal distortion accompanying the ordering reaction leads to the formation of self-organized tweed patterns of coexisting (cubic) A1- and (tetrago­nal) L10-phases. As seen from Figure 2, the agreement between experimental and simulated microstructures is remarkable. Ni and Khachaturyan proposed recently that, in order to minimize elastic energy during trans­formations involving symmetry changes and lattice strain, a pseudospinodal decomposition is likely to take place, leading to 3D chessboard patterns.41

PF modeling has also been used extensively to study martensitic transformations,34-38,42-44 phase transformations in ferroelectrics45-57 (see also the recent review by Chen58 on that topic), transforma­tions in thin films,47,59-65 grain growth and recrystal­lization,66-81 and microstructural evolution in the presence of cracks or voids.82-84 A recent extension of PFMs has been the inclusion of dislocations in the models,85-87 by taking advantage of the equivalence between dislocation loops and coherent misfitting platelet inclusions.88 This approach has been applied, for instance, to study the interaction between moving dislocations and solute atoms,89 or to study the influ­ence of dislocation arrays on spinodal decomposition in thin films.61 Rodney et a/.87 have pointed out,

image206

Подпись: IsotropicПодпись: Interface onlyimage929Strain only

(c)

image930
Int +strain + kinetics Experiment

(b) (c)

Подпись:Time

(d) (e) (f)

Figure 2 Comparison between transmission electron microscopy experimental observations (a-c) and phase field modeling (d-f) of formation of chessboard pattern in Co395Pt605 cooled from 1023 K to (a) 963 K, (b) 923 K, and (c) 873 K. The scale bar corresponds to 30 nm. Reproduced from Le Bouar, Y.; Loiseau, A.; Khachaturyan, A. G. Acta Mater. 1998,46(8), 2777-2788.

image932however, that the artificially wide dislocation cores required by the above approach lead to weak short — range interactions. These authors have introduced a different PFM for dislocations, which allows for nar­row dislocation cores. As an illustration of that model, Figure 3 shows the development of dislocation loops and their interaction with hard precipitates in a 3D y/y’ single crystal. It is interesting to note that dislo­cation loop initially expands by gliding in the soft y channels, until the local stresses are large enough for the dislocation to shear the hard y’-phase.

The above presentation of the PF equations leaves certain questions open. First, the maximum homoge­neous free energy difference А/max involves the free energy of the unstable state separating the two minima at Ca and Cp. It is thus been questioned90 whether this quantity can be rigorously defined from thermody­namic principles. If one employs mean field techni­ques such as the cluster variation method (CVM)91-95 to derive the homogeneous free energy of an alloy, А/max is in fact very sensitive to the approximation used, and generally decreases as the size of the largest cluster used in the CVM increases.96 Kikuchi97-99 has argued that, in order to resolve this paradox, А/ should not be considered as the free energy of any homogeneous state, but that it should be understood as the local contribution to the free energy of the system along the equilibrium composition profile.

The second set of questions relates to the gradient energy coefficients к and q in eqn [5]. In many applications of PF modeling, these coefficients are taken as phenomenological constants that can be adjusted at will, as long as the microstructures are scaled in units of к1/2 or q1/2. Such an approach, however, is problematic for many reasons. First, when one scalar field variable is employed, for instance C(r/), a regular solution model,2, 0 or equiv­alently a Bragg-Williams approximation,101 estab­lishes that the gradient energy coefficient is not arbitrary but that it is directly proportional to the interaction energy between atoms, that is, to the heat of mixing of the alloy. Furthermore, in the most general case, к should in fact be composition and temperature dependent. Starting from an atomistic model, rigorous calculations of к are possible by monitoring the intensity of composition fluctuations as a function of their wave vector, and using the fluctuation-dissipation theorem.1 0 In the case of a simple Ising-like binary alloy, it is observed that к varies as C(1 — C), where C is the local composition of the alloy. 0 Furthermore, when more than one field variable is employed, care should be taken to consider all possible contributions of field heterogeneities to the free energy of the system, as the different fields may be coupled. Symmetry considerations are important to identify the nonvanishing terms, but it

may remain challenging to assign values to these nonvanishing terms that are consistent with the ther­modynamics of the alloy considered.

Another important point is that interfacial ener­gies are in general anisotropic. In order to obtain realistic morphological evolution, it is often impor­tant, and sometimes even absolutely necessary, to include this anisotropy, for example, in the modeling of dendritic solidification. The symmetry of the mesh chosen for numerically solving the PF equations introduces interfacial anisotropy but in an unphysical and uncontrolled way. One possible approach to introduce interfacial anisotropy is to let к vary with the local orientation of the interface with respect to crystallographic directions.11, Another approach is to rely on symmetry constraints27,30 to determine the number of independent coefficients in a general expression of the inhomogeneity term, see eqn [5]. In both approaches, the different coefficients entering
the interfacial anisotropy can be fitted to experiments or to atomistic simulations.

Let us now return to the mobility coefficients M and L introduced in eqns [2] and [3]. For the sake of simplicity, many PF calculations are performed while assigning an arbitrary constant value to these coeffi­cients. An improvement can be made by relating the mobility to a diffusion coefficient. In the case of M, for instance, in order to make eqn [1] consistent with Fick’s second law for an ideal binary alloy sys­tem, one should choose

C)

-Б where C is the average solute concentration and D the interdiffusion coefficient. In both cases, the simulated times are expressed in arbitrary units of M-1 or L — , thus precluding a direct connection with experimental kinetics. This problem is also directly related to the lack of absolute physical length scales in these simulations. Moreover, using a 1D Bragg-Williams model composed of atomic planes, Martin101 showed that M is not a constant but is in fact a function of the local composition along the equilibrium profile. A complete connection between atomistic dynamics and M will be made in Section 1.15.3. Similar to the discussion on coupling between various fields for the gradient energy terms, kinetic coupling is also expected in general. The kinetic couplings between composi­tion (a conserved order parameter) and chemical ordering (a nonconserved order parameter) are revealed by including sublattices into Martin’s 1D model and deriving the macroscopic evolution of the fields from the microscopic dynamics. In that case, atoms jump between adjacent planes.102,103 As a result, instead of the mere superposition of eqns [2] and [3], the kinetic evolution of coupled con­centration and chemical order in a binary alloy is given by

dC(r, t) dt

+ V( M-’V( Щ5))]

+ V(L-’V

+ V L3 V where L1 is a mobility coefficient, and L2, L3, M1, and M2 are second-rank mobility tensors, since they

image422

relate diffusional fluxes (vectors) to chemical potential gradients (vectors). In the case of cubic crystalline phases, second-rank tensors reduce to scalars, but in many ordering reactions, noncubic phases form, thus leading to anisotropic mobility. Vaks and coworkers104 have also derived PFMs for simultaneous ordering and decomposition starting from microscopic models. These works, however, illustrate the fact that it would be quite difficult, especially for multidimensional field variables, to assign correct values to the kinetic coefficients for a given alloy system by relying solely on a phenomeno­logical approach.