Incorporating the CALPHAD sublattice approach of ordering into the phase-field model with finite interface dissipation

Available online at www.sciencedirect.com

ScienceDirect Acta Materialia 88 (2015) 156–169 www.elsevier.com/locate/actamat

Incorporating the CALPHAD sublattice approach of ordering into the phase-field model with finite interface dissipation ⇑

Lijun Zhang,a,b Matthias Stratmann,a Yong Du,b Bo Sundmanc and Ingo Steinbacha, a

Interdisciplinary Centre for Advanced Materials Simulation (ICAMS), Ruhr-Universita¨t Bochum, D-44780 Bochum, Germany b State Key Laboratory of Powder Metallurgy, Central South University, 410083 Changsha, PR China c Institut National des Sciences et Techniques Nucle´aires (INSTN), CEA Saclay, France Received 1 May 2014; revised 17 November 2014; accepted 23 November 2014

Abstract—A new approach to incorporate the sublattice models in the CALPHAD (CALculation of PHAse Diagram) formalism directly into the phase-field formalism is developed. In binary alloys, the sublattice models can be classified into two types (i.e., “Type I” and “Type II”), depending on whether a direct one-to-one relation between the element site fraction in the CALPHAD database and the phase concentration in the phase-field model exists (Type I), or not (Type II). For “Type II” sublattice models, the specific site fractions, corresponding to a given mole fraction, have to be established via internal relaxation between different sublattices. Internal minimization of sublattice occupancy and solute evolution during microstructure transformation leads, in general, to a solution superior to the separate solution of the individual problems. The present coupling technique is validated for Fe–C and Ni–Al alloys. Finally, the model is extended into multicomponent alloys and applied to simulate the nucleation process of VC monocarbide from austenite matrix in a steel containing vanadium. Ó 2014 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Phase-field model; Computational thermodynamics; Kinetics of phase transformations; CALPHAD modelling

1. Introduction The phase-field method, being based on a diffuse representation of the phase boundaries, has emerged as a powerful tool to simulate the microstructural evolution in various materials processes during their lifetime [1–6]. Besides a consistent theoretical basis and good numerical solutions, the input thermophysical parameters are decisive for quantitative simulations related to real materials. Among various thermophysical input parameters needed for quantitative simulations, the bulk free energy and its derivatives, diffusion potential and thermodynamic factor, are key parameters of all phase-field simulations. These entities determine the driving forces for phase-field evolution and diffusion. Nowadays, the coupling to the CALPHAD (CALculation of PHAse Diagrams) technique has become a good standard [7–9] for providing reasonable thermodynamic information needed in the phase-field model. Up to now, a variety of coupling approaches have been reported for phase-field simulations of solidification and solid-state transformation. In general, those approaches roughly split into two classes. One class can be referred to as indirect methods, including i) the link to a data file containing the necessary thermodynamic information such

⇑ Corresponding author; e-mail: [email protected]

as Gibbs energy, chemical potential and phase equilibria calculated with the help of an external CALPHAD software package [10,11]; ii) the approximation of Gibbs energy curve/surface from the CALPHAD calculations via some simple functions [12,13], e.g., parabolic equations and Landau polynomials, and iii) the use of locally linearized phase diagrams. The local linearization, or partitioning, is provided by running a CALPHAD software package in parallel to the phase-field simulation, e.g. in the commercial package MICRESS (MICRostructure Evolution Simulation Software) [14,15]. The other class can be addressed as “direct methods” in which the Gibbs energy functions with CALPHAD thermodynamic parameters stored in thermodynamic databases are incorporated directly into the phase-field model. This method has been applied in substitutional solution phases where the mole fractions as used in CALPHAD-type Gibbs energy functions are directly related to concentrations in phase-field models via a given molar density [16]. This method, however, is not applicable for intermetallic compounds and interstitial solution phases, which are usually described by the sublattice models [17,18]. The sublattice models are a subset of the Compound Energy Formalism (CEF) and use site fractions to describe the variation of constituents on different sublattices. The site fractions are denoted y is where i is the constituent

http://dx.doi.org/10.1016/j.actamat.2014.11.037 1359-6462/Ó 2014 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.

L. Zhang et al. / Acta Materialia 88 (2015) 156–169

and s is the sublattice. For each sublattice the sum of the site fractions are unity X y is ¼ 1: ð1Þ i

The mole fractions can be calculated from the site fractions using P P s as j bij y js xi ¼ P P P ; ð2Þ a s k s j bkj y js where xi is a mole fraction of a component i; as is the number of sites on sublattice s; bij is the stoichiometric factor of component i in the constituent j. The Gibbs energy function for each phase is modelled as a function of T ; P and the site fractions, y is . For a given set of mole fractions one must in general minimize the Gibbs energy in order to find the site fractions. If this minimization is made by an external procedure, it will influence the numerical performance and block the direct connection of the phase-field models and the CALPHAD modelling. The concentration, ci , of component i, used to define the phase-field equations can be transformed to a mole fraction using the following relation: ci ¼

N i NNi xi ¼V ¼ ; V V m N

ð3Þ

where N i is the number of moles of component i; V the volume, N the total amount of moles and V m the molar volume. The main advantage of using mole fractions instead of concentrations is that they are independent of T and P whereas concentrations are not. Even if V m is frequently assumed constant in most phase-field simulations, the introduction of mole fractions in the phase-field equations will avoid problems when extending the phase-field technique to cases when one has variable volumes or allows the volume to depend on T. In this case, and if a phase transformation is associated to large Bain strain, the mechanical equilibrium has to be considered as well. An approach to cope with large strain has been recently developed [19]. Another advantage of using mole fractions instead of concentrations is that all CALPHAD modelling of the Gibbs energy use mole fractions or site fractions. Zhu et al. [20] were the first who incorporated the site fraction of a four-sublattice model into the phase-field model by converting the element site fraction to an order parameter for the c=c0 solid state transformation in the Ni–Al system. This approach has been further extended by Kitashima et al. [8,21] to multi-component Ni-based superalloys. This kind of coupling approach also avoids the equilibrium calculations and energy minimization during the simulation. The approach, however, had been limited up to now to the c=c0 order/disorder transition. The purpose of the present treatment is to provide a general framework of incorporating CALPHAD sublattice modeling into a phase-field model. It must be stated clearly that only a common minimization of all thermodynamic variables under the constraint of the local state of the microstructure (e.g. curvature of interfaces) leads to a comprehensive solution. Separate minimization, e.g. sublattice occupancy and depletion of solute against equilibrium in a nano-scale precipitate due to high curvature, leads to an erroneous result, in general.

157

Very recently, a phase-field model with finite interface dissipation was developed with major contributions from two of the present authors [22,23] in the framework of the multi-phase-field (MPF) formalism. This approach provides the description of various kinetic processes at the mesoscopic scale without restriction to the type of transformation ranging from the chemical equilibrium to strongly non-equilibrium phase transformations. The novel feature of the model is that each phase concentration is assigned by a kinetic equation to account for finite interface dissipation instead of applying an extra condition for solute partitioning between the phases as in traditional models: the condition of a given partitioning [24,25] or the condition of equal diffusion potentials [15,26]. With such a novel feature, the external equilibrium calculation for the partitioning at the interface can be avoided in phase-field simulations. The thermodynamic potentials can be directly incorporated from a CALPHAD thermodynamic database. In our previous work, the coupling of the phase-field model with finite interface dissipation and thermodynamic parameters in the CALPHAD formalism has been demonstrated for substitutional solution phases [22,23,27]. In the present work, a generalization to materials which are described by the sublattice formalism is developed and applied to several different alloys and types of transformations. In the following section, we will give a brief introduction of the phase-field model with finite interface dissipation as well as the sublattice model in CALPHAD formalism. In Section 3, the approach for coupling of the phase-field model with finite interface dissipation and different sublattice models in binary alloys is presented. After that, the validation of the present coupling approach is made in both Fe–C and Ni–Al alloys in Section 4. In Section 5, the coupling technique is extended into multicomponent compounds, and applied to simulate the nucleation process of VC monocarbide from austenite matrix in Fe–V–C alloys.

2. Model description 2.1. Phase-field model with finite interface dissipation Assuming a multicomponent alloy with a ¼ 1 . . . q phases and i ¼ 1 . . . n components, the total free energy a functional F is as a function  a of the phase fields f/ g and the phase concentrations ci . The phase-fields fulfil the sum convention and the phase concentrations sum up to the mixture concentration: q X a¼1 q X

/a ¼ 1;

ð4Þ

/a cai ¼ ci :

ð5Þ

a¼1

Here, in order to be consistent with the standard notation in thermodynamics, the phase will always be given as a superscript. F integrates the interfacial energy density s and the chemical energy density f over a volume X: Z F ¼ ð6Þ fs þ f gd 3 v; X

where the integrand within the domain X at time t, with q phases coexisting, is a sum of the gradients of the phase fields and the potential function:

158

L. Zhang et al. / Acta Materialia 88 (2015) 156–169

sðf/a gÞ ¼

q1 X q  X a¼1 b¼aþ1

 2  4rab g a b a b  2 r/  r/ þ j/ / j ; g p

ð7Þ

where rab is the interfacial energy between the phases/ grains of a and b in multiphase junction with q phases/ grains. g is the interface width, and a chemical free energy density given by: f ¼

q X /a f a ðxai Þ;

ð8Þ

a¼1

where we expressed to compositions using mole fractions, xai ¼ V m cai . The f a ðxai Þ is the bulk volume free energy density of phase a evaluated for the phase composition xai . As already noted in the introduction, it is preferably to use composition variables that are independent of T and, in order to use the CALPHAD modelling, the free energy as a molar Gibbs energy: f a ðxai Þ ¼

Gam ðT ; P ; xa Þ V cm

ð9Þ

~ai and l ~bi for solute and we define the diffusion potentials l component i in a and b phases, relative to the solvent, component n [28]:  a  a @Gm @Gm ~ai ¼  ; ð10Þ l @xi T ;P ;xj–i @xn T ;P ;xj–n where Gam is the molar Gibbs energy for a. The standard mathematical notation is used to indicate the variables kept constant when calculating the derivative. The phase-field evolution is derived from the principle of energy minimization [23,29]:   q q X 4g mab d d 1 X _ ab F ¼ /_ a ¼  w ;  ð11Þ a d/ q b¼1 qp2 d/b b¼1 ab

where m is the mobility of the interface between phase a and b. This generalized Onsager relation is based on the decomposition of a multi-phase interaction in a junction between three or more phases into pairwise contributions w_ ab [29]. It can be easily proven that this formulation reduces to the standard dual phase-field equation for dual interfaces. In order to evaluate the pair interactions consistently with the mole fractions of solute in the restricted phase space of the pair, we introduce the ‘pair-composition’ a a b b xab i ¼ / xi þ / xi , where we have in multiple junctions a b / þ / < 1 in general. The demand of conservation of the pair-composition follows from the demand that all pairs shall be varied independently while keeping all other phases and phase compositions constant. We therefore introduce the Lagrange function for the pair Lab Z  ab ab  Lab ¼ F þ k ðxi  /a xai  /b xbi Þ d 3 v: ð12Þ X

The Lagrange multiplier kab is evaluated from the demand of conservation of the pair compositions as [23] kab i ¼

~ai þ /b l ~bi w_ ab ðxa  xb Þ /a l  ab i a i b : a b / þ/ qP i ð/ þ / Þ

ð13Þ

P ab is the so-called interface permeability for the dual i interaction a–b, i.e. the inverse resistivity of the interface

against redistribution of solute. If no special resistivity against redistribution of solute at the interface is given as, e.g. in electrochemical transformations [30], it can be estimated from the bulk chemical mobilities P ab i ¼

~ a þ /b V b M ~ bi Þ 8ð/a V am M m i ; ag

ð14Þ

~ a=b where M is the chemical mobility for component i in the i a / b phase respectively. a is a microscopic length scale, i.e., width of the interface on atomistic scale. One might include an excess volume of the phase boundary into the definition of the permeability. Then the phase or grain boundary motion is also related to mechanical distortion, as treated in [31,32]. The special model for the permeability will in general depend on many details of the system as well as the phase boundary orientation. An alternative treatment based on a maximum entropy production principle is presented in [33]. The phase-field evolution Eq. (11) can now be evaluated (" # ) q ab X p2 ab ab a b bc ac c _wab ¼ K r ðs s Þ þ ðr  r Þs þ Dg q 4g c¼1;c–a;c–b ð15Þ with nthe independent variations of the surface terms o 2 a p2 a a s ¼ r / þ g2 / . K ab ¼

4qgð/a þ /b Þmab ; Pn1 ðxai xbi Þ2 a b ab 2 4qgð/ þ / Þ þ m p i¼1 P ab

ð16Þ

i

and Dgab ¼ Gbm  Gam 

n1 a a X ~ þ /b l ~bi /l i

i¼1

/a þ /b

ðxbi  xai Þ;

ð17Þ

where mab is the interfacial mobility between a and b; I a is the capillary contribution for phase a; Dgab is the driving force for the phase transformation, and K ab is the renormalized interfacial mobility considering the influence of partitioning. The evolution equations for phase compositions expressed as mole fractions xai (independent of T and P) are in diagonal approximation neglecting cross terms in the chemical mobilities [34]: ~ a r~ /a x_ ai ¼ V am rð/a M lai Þ þ i þ

q X

q X /a /b ~ai Þ P ab ð~ lbi  l i a b / þ / b¼1

/a w_ ab ðxbi  xai Þ: a b qð/ þ / Þ b¼1

ð18Þ

Three mechanisms contribute to the concentration evolution for each phase: spatial concentration (or chemical potential) gradient in each phase driving diffusion fluxes, differences in diffusion potential driving partitioning, and phase transformation affecting the local reference configuration. Actually, the summation of Eq. (18) over all q phases will lead to the standard Fick’s law for the mixture concentration ci . The nontrivial point is that this suspends the need to employ an extra condition to fix the concentrations for each phase. Instead, one can use the separated concentration evolution equations of each phase for the iteration, which is applicable for arbitrary initial conditions.

L. Zhang et al. / Acta Materialia 88 (2015) 156–169

2.2. Sublattice model in CALPHAD formalism In the framework of CALPHAD, the Gibbs energy of each phase is described by an appropriate thermodynamic model dependent on its physical and chemical properties, for instance, crystallography, type of bonding, order–disorder transitions, as well as magnetic properties. For intermetallic compounds and interstitial solution phases, the sublattice model is usually the choice.1 For a fictitious phase h described by the sublattice model ðA; BÞk ðD; E; F Þl ðG; H Þm , it means that A and B mix on the first sublattice, D; E and F mix on the second one, and G and H mix on the third one. The coefficients k; l and m are P the number of sites on the sublattices we assume that k þ l þ m ¼ 1, i.e. there are totally 1 mol of sites. The general form for the total molar Gibbs energy of phase h is given as [17]. Ghm ¼ srf Ghm  T cnf S hm þ E Ghm þ mag Ghm ; srf

ð19Þ

Ghm

where represents the molar Gibbs energy of an unreacted mixture of the end-members, I of the phase: XY srf h y is2I GI ; Gm ¼ ð20Þ I

where the summation over I is for all end-members. An end-member I specifies one constituent i in each sublattice and the notation y is2I means the site fraction of constituent i in sublattice s as specified by I. As already given in Eq. (1), the sum of the site fractions in each sublattice is unity.  GI is the molar Gibbs energy of the end-member I, and its value can be either assessed from the experimental phase equilibria and/or thermodynamic properties or directly calculated from the first-principles calculations [35,36]. In the second term of Eq. (19) T is the absolute temperature, and the term cnf S hm denotes the molar configurational entropy of the phase, and is based on the number of possible arrangements of the constituents in the phase h. In CEF one assumes random mixing on each sublattice which it is expressed as X X cnf h S m ¼ R as y is lnðy is Þ; ð21Þ s

i

where R is the gas constant. The term E Ghm in Eq. (19) is the molar excess Gibbs energy and it can have many different forms. These as well as the magnetic contribution to the molar Gibbs energy, mag Gm are described in detail in the book by Lukas et al. [36]. In the general case a constituent i can represent any type of species, i.e., atoms, molecule, ion or vacancy, and thus the total number of atoms in phase h, should be P s as ð1  y Va;s Þ with y Va;s as the site fraction of vacancies on sublattice s. Here, as is the general notation for the number of sites of the sublattice s, viz. k; l or m. It is important to realize that the mole fractions, xai , are for the components of the system whereas the constituent fractions, y ais , are for the constituents and one can have many more constituents than one has components. In the case each constituent consists of a single component, but can be present in several sublattices, the relation, Eq. (2), between the mole fraction and the site fraction becomes 1

159

P

One exception is the so-called stoichiometric or “line” compound with no or extremely negligible solubility. This kind of description is not included here.

Nh ah y h xhi ¼ P i h ¼ P h s s is h : s as ð1  y Va;s Þ jN j

ð22Þ

2.3. Chemical potentials from molar Gibbs energy models In a typical CALPHAD database, all the parameters needed to calculate Ghm according to Eq. (19) as a function of T ; P and the site fractions are stored and from these the equilibrium state can be calculated as well as various thermodynamic properties. With some care this can also be made for a single phase in a metastable state. A property that is of great interest for phase-field simulations is the chemical potential which can be defined from the Gibbs energy G and the total amount of the components N i as:  @G li ¼ ; ð23Þ @N i T ;P ;N j–i where the T ; P and the amount of all other components must be kept constant. For a model of the molar Gibbs energy for a phase h with several sublattices the definition of the chemical potential for a component can be transformed [17,18] to the chemical potential of an end-member I using the molar Gibbs energy Ghm : ! ! X @Gh XX @Gh h h m m lI ¼ Gm þ  ; ð24Þ @y i;s2I @y j;s s s j T ;P ;y j–i;s

T ;P ;y k–j

where the first summation is for one constituent in each sublattice as defined by the end-member I and the second double summation is for all constituents in all sublattices. Note that the derivatives are calculated keeping all other constituent fractions constant. Even using a model with mole fractions, xi , the chemical potential must be calculated summing several partial derivatives:  h X@Gh @Gm m lhi ¼ Ghm þ  : ð25Þ @xi T ;P ;xj–i @x j T ;P ;xk–j j The reason why one must sum several derivatives is that mole fractions as well as constituent fractions are not independent variables whereas the total amounts of components are. Thus one cannot calculate a derivative of the molar Gibbs energy with respect to a single fraction variable to find the chemical potential, but must take into account the effect of the other fractions variables. A diffusion potential for two components, like in Eq. (10), can be written as a difference between two partial Gibbs energies according to Eq. (25) as many terms cancel. But the fact that  h  h @Gm @Gm li  lj ¼  ð26Þ @xi T ;P ;xk–i @xj T ;P ;xk–i

h m does not mean that li ¼ @G . @xi T ;P ;xk–i

An end-member I represents a compound with fixed stoichiometry as given by the sublattices. Provide that the constituents on the sublattices are the same as the components one can at equilibrium write the Gibbs energy for this compound as the following sum of the chemical potentials of the components:

160

GI ¼

L. Zhang et al. / Acta Materialia 88 (2015) 156–169

X

as lis2I :

ð27Þ

s

Even when there is a small deviation from equilibrium this relation holds. Applied to a model like ðA; BÞk ðD; E; F Þl ðG; H Þm this means GA:D:G ¼ klA þ llE þ mlG

ð28Þ

3. Coupling technique in different binary sublattice models As stated in Section 1, the direct coupling of the phasefield model with finite interface dissipation and the CALPHAD databases relies on substituting the volume free energy densities, like f a and f b , as well as the diffusion potentials, ~ai and l ~bi , in the evolution Eqs. (11)–(18) by the molar like l Gibbs energies and the chemical potentials shown in Sections 2.2 and 2.3 with the same quantities based on sublattices and molar fractions. This requirement is easy to achieve in substitutional solution phases because the one-to-one relation between the phase concentration cai in the phase-field model and the mole fraction xai in the CALPHAD database can be simply established by assuming the constant molar volume for each species, as defined in Eq. 3. However, for the phases described by the sublattice models, the situation becomes complex because the site fractions are the variables in CALPHAD database. Though the relation between the phase concentration and the site fractions can be established via both Eqs. (3) and (22), in some cases the site fractions corresponding to the given phase concentrations are unknown before the equilibrium calculations by means of the energy minimization procedure. An external minimization of the sublattice occupancy, e.g. via the minimizer integrated in a thermodynamic software, is only a substitute, since it neglects the influence of the morphology of the phases on sublattice occupancy. In fact, we will demonstrate how the curvature of interfaces in small precipitates influences the sublattice occupancy significantly, the effect of which should not be neglected. Therefore a direct and efficient coupling technique needs to be developed. We will start out from simple binary systems. In a fictitious A–B binary alloys, the frequently used sublattice models are ðA; BÞk ðAÞl , ðAÞk ðB; VaÞl , and ðA; BÞk ðA; BÞl . Actually, the three types of sublattice models can be classified into two groups. In the first two sublattice models, a one-to-one relation between the phase concentration in the phase-field model and the site fraction in the CALPHAD database can be established via Eqs. (3) and (22). These models are named as “Type I” sublattice models. As for the third sublattice model, the establishment of such one-to-one relations needs a energy minimization procedure in addition to Eqs. (3) and (22), and this kind of model is named “Type II” sublattice model. In the following, A is generally assumed to be the solvent and B is the solute. 3.1. “Type I” sublattice models in binary alloys 3.1.1. Model ðA; BÞk ðAÞl This is the prototype of a series of intermetallic compounds with the homogeneity extending into A-rich side. The one-to-one relation between the phase concentration chB in the phase-field model and the site fractions, y uA and y uB , in CALPHAD database according to Eqs. (3) and (22) is given as

8 V m chB > > < y uB ¼ k ð29Þ h > > u : y ¼ 1  V m cB A k by considering that ðk þ lÞ ¼ 1. With Eq. (29), the site fraction for each species in model ðA; BÞk ðAÞl is known for a given phase concentration chB . Then, the Gibbs energy shown in Eqs. (19) can be directly used in the phase-field models via Eq. (9). Now, the remaining are the diffusion potentials that are needed in Eqs. (17) and (18). Though the diffusion potential can be obtained by the first derivative of the Gibbs energy over the phase concentration, this method is still complex, and even much more complex in multicomponent alloys. The alternative way is to employ the approach used in the CALPHAD software packages [37]. According to Eq. (28), we have lhBk Al  lhAk Al ¼ kluB þ llvA  kluA  llvA ¼ kðluB  luA Þ  h  h @Gm @Gm ¼  : u @y B yu @y uA y u lhB

ð30Þ

B

A

Furthermore, we have the constraint that  lhA ¼ luB  luA

ð31Þ

when the system is in internal equilibrium [18]. Based on Eqs. (30) and (31), we can obtain the diffusion potential for phase h described by model ðA; BÞk ðAÞl in phase-field notation as lhB  lhA luB  luA ¼ Vm Vm (  h ) h 1 @Gm @Gm : ¼  u V mk @y B yu @y uA yu

~hB ¼ l

ð32Þ

B

A

Based on Eq. (29), the diffusion potentials needed in Eqs. (17) and (18) can be directly taken from Eq. (32). 3.1.2. Model ðAÞk ðB; VaÞl This model is widely used to describe the interstitial solution phase, like austenite and ferrite in steel. In this case, B is usually the interstitial atom, i.e., C or N. The one-to-one relation between the phase concentration chB in the phase-field model and the site fractions, y vB and y vVa , according to Eqs. (3) and (22) is given as 8 kV m chB > > > y vB ¼ > < lð1  V m chB Þ ð33Þ > > l  V m chB > v > : y Va ¼ lð1  V m chB Þ by considering that ðk þ lÞ ¼ 1. With Eq. (33), the site fraction for each species in model ðAÞk ðB; VaÞl will be known for a given phase concentration chB , and thus the Gibbs energy showing in Eqs. (19) can be directly used in the phase-field models via Eq. (9). According to Eq. (28), one can obtain  h  h ! X @Gm h h u v v @Gm lAk Val ¼ klA þ llVa ¼ Gm þ  yi ; ð34Þ @y vVa y v @y vi y v i j

B

lhAk Bl

¼ kluA þ llvB

¼ Ghm þ



@Ghm @y vB

  y vVa

X i

y vi



@Ghm @y vi

 ! : y vj

ð35Þ

L. Zhang et al. / Acta Materialia 88 (2015) 156–169

Based on the condition that lvVa ¼ 0, which exists at thermal equilibrium [18], we can get luA as (  h  h  !) X 1 @Gm h u v @Gm Gm þ : ð36Þ lA ¼  yi k @y vVa y v @y vi yv i j

B

Furthermore, the difference between Eq. (34) and Eq. (35) leads to  h  h @Gm @Gm lhAk Bl  lhAk Val ¼ lðlvB  lvVa Þ ¼  ; ð37Þ v @y B y v @y vVa y v Va

B

from which one can easily obtain lvB as (   h ) 1 @Ghm @Gm v : lB ¼  l @y vB y v @y vVa y v Va

ð38Þ

B

On the basis of Eqs. (37) and (38), the diffusion potential for phase h described by model ðAÞk ðB; VaÞl in phase-field notation can be obtained as (   h ) lhB  lhA lvB  luA 1 @Ghm @Gm h ~B ¼ ¼ ¼  l V ml Vm Vm @y vB yv @y vVa y v Va B ( !)  h   h X 1 @G @G m m Ghm þ : ð39Þ  y vi  V mk @y vVa yv @y vi y v i j

B

With Eq. (33), the diffusion potential shown in Eq. (39) can be directly used in the phase-field evolution Eqs. (17) and (18). 3.2. “Type II” sublattice models in binary alloys Phases described by the sublattice model ðA; BÞk ðA; BÞl represent a large number of intermetallic compounds with certain homogeneity range in CALPHAD formalism. Its Gibbs energy can be described by the general energy Eqs. (19). According to Eq. (28), we can obtain the following equations  h @Gm h h u v h lAk Al ¼ klA þ llA ¼ lA ¼ Gm þ @y uA y s i  h   ! XX @Gm @G m ; ð40Þ þ  y si @y vA y s @y si ys s i i

j

lhBk Bl ¼ kluB þ llvB ¼ lhB ¼ Ghm þ  ! XX s @Gm :  yi @y si y s s i 



h

@Gm @y uB

 þ y si

@Ghm @y vB

 y si

ð41Þ

j

From the difference between Eq. (40) and Eq. (41), we get the diffusion potential for phase h described by model ðA; BÞk ðA; BÞl in phase-field notation, lhB  lhA Vm (   h  h  h ) 1 @Ghm @Gm @Gm @Gm : ¼  þ  Vm @y uB y s @y uA y s @y vB y s @y vA y s

l ~hB ¼

i

i

i

ð42Þ

i

Based on the constraint when the system is in internal equilibrium [18], that is, lhB  lhA ¼ luB  luA ¼ lvB  lvA , we alternatively write

161

 h  h !  h  h ! 1 @Gm @Gm 1 @Gm @Gm h l ~B ¼   ¼ : V mk vm l @y uB y s @y uA y s @y vB y s @y vA y s i

i

i

i

ð43Þ

Though the free energy densities and diffusion potentials in the phase-field models can be related to the Gibbs energy and chemical potentials in CALPHAD database, the oneto-one relation between the phase concentration chB and the site fractions cannot be directly derived for the sublattice model ðA; BÞk ðA; BÞl , as did in “Type I” sublattice models. That is because in the phase described by ðA; BÞk ðA; BÞl we have chB ¼

ky uB þ ly vB : Vm

ð44Þ

In order to determine the pair of y uB and y vB for a given chB , the minimization process for the Gibbs energy in the phase is needed. In the present approach, however, we allow relaxation of the sublattice occupancies (site fractions) with finite rate. In our previous work [22,23], the space was formally divided into distinct reference volumes RVs. A further assumption is that the internal processes inside a RV and the exchange of atoms between adjacent RVs are independent and can be superimposed. By considering the equilibration between two phases inside a RV, the Lagrange multiplier is then determined, from which the total energy in arbitrary thermodynamic states can be constructed. Following the same idea, we will demonstrate the new way to determine the pair of y uB and y vB for a given chB . For a two-phase binary system consisting of a substitutional solution phase a, and an intermetallic compound b described by the sublattice model ðA; BÞk ðA; BÞl , its RV is schematically demonstrated in Fig. 1. In each RV, the internal partition processes between a and b phases and the exchange of atoms between adjacent RVs are also independent and can be superimposed. Furthermore, the intermetallic compound b is divided into two sub-RVs, which correspond to the two sublattices in the model, respectively. In the first sublattice ðA; BÞk ; A is the majority species, while in the second sublattice ðA; BÞl ; B is the majority species, as is the most common case for such two-sublattice model in CALPHAD databases. In sub-RVs of the phase b, we assume that the atoms can exchange between the first sublattice and the second sublattice to achieve the internal equilibrium state, but the atoms can only interact with those in a phase as an integral phase. It means that the entire partitioning process can be divided into two separate steps. In the first step, the atoms in a phase and integral b phase redistribute via the evolution Eqs. (11)–(18) in the phase-field models. Then, new phase concentrations can be obtained for both a and b phases. While in the second step, the atoms in both sublattices of b phase exchange with each other to achieve the internal equilibrium for the known phase concentration from the first step, that is, the one-to-one relation between the phase concentration and the site fractions is established. At this step, the sub-RVs are treated as insulated from outside, and the solute contents (i.e., site fractions) in each sublattice may adjust with finite rates for equilibrium values. It should be noted that the redistribution of the components on the different sublattices in a single phase within a RV with known overall composition temperature and pressure can be solved by the standard equilibrium calculation technique as explained by Hillert [38]. However, in order to take into account that

162

L. Zhang et al. / Acta Materialia 88 (2015) 156–169

RV

Aatom Batom

l_y vB ¼ 

αphase AB

ð49Þ

Summing up the above two equations, we can obtain  k y_ u þ l_y vB M @Gbm @Gbm b c_ bB ¼ B ¼ 2 þ  k : ð50Þ Vm @y vB V m @y uB In the isolated sub-RVs, the phase concentration cbB at the second step is independent of time, c_ bB ¼ 0. Then we can solve Eq. (50) for kb

1stsublattice

2ndsublattice

βphase

kb ¼

@Gbm @Gbm þ v : @y uB @y B

ð51Þ

Considering the fact that only the site fraction of solute in each sublattice is independent in the phase-field notation, and carefully comparing Eq. (51) and Eq. (42), we find

ABAB

Fig. 1. Schematic demonstration of reference volume RV in a twophase binary system consisting of a substitutional solution phase a, and an intermetallic compound b described by the sublattice model ðA; BÞk ðA; BÞl . The intermetallic compound b is divided into two subRVs, corresponding to the two sublattices in the model, respectively.

this process is not instantaneous one may introduce a timedependent relation as explained below. In this picture, the total volume chemical free energy in b phase can be expressed as    Z  k u l v dX: ð45Þ Fb ¼ f b þ kb cbB  yB þ yB Vm Vm X Here, A is assumed to be the dominant solvent in b phase in accordance with the phase-field notation, while f b is the chemical free energy density of b phase, and can be related b

to the molar Gibbs energy Gbm via f b ¼ GV mm . The molar Gibbs energy Gbm is described by Eqs. (19) in Section 2.2. The Lagrange multiplier kb , as a generalized diffusion potential, is introduced to ensure solute conservation. cbB is the phase concentration of solute B in b phase, while y uB and y vB are the site fractions of B in the first and the second sublattice, respectively. At the second step during the entire partitioning process, the rate of change of the site fractions, which are treated as independent functional variables in the sense of the Lagrange formalism, are derived from a variational principle for non-conserved order parameters  b  d @f k b M u @Gbm b ¼  ;  k  kk y_ uB ¼ M u u F b ¼ M u dy B @y uB V m V m @y uB ð46Þ y_ vB ¼ M v

 M @Gbm b :  lk V m @y vB





b

v

d b @f l M F ¼ M v  kb ¼  dy vB @y vB V m Vm



@Gbm b  lk ; @y vB ð47Þ

u

v

with the relaxation coefficients M and M . For simplicity we assume that the redistribution of solutes between the two sub-RVs only depends on the diffusion potentials in each sublattice with a common relaxation coefficient M, which is defined as kM u ¼ lM v ¼ M. M has the dimensions 3 . This yields the rate equations as of cm Js  M @Gbm b u ð48Þ k y_ B ¼   kk ; V m @y uB

~bB : kb ¼ V m l

ð52Þ b

It indicates that k is exactly the diffusion potential in b phase. Substituting Eq. (52) into Eqs. (48) and (49), we finally obtain the evolution equations for site fractions  1 @Gbm b ð53Þ y_ uB ¼ M  l ~ B ; V m k @y uB y_ vB ¼ M



1 @Gbm b ~  l B : V m l @y vB

ð54Þ

With a given phase concentration cbB from the phase-field simulation, the corresponding site fractions y uB and y vB can be obtained via Eqs. (53) and (54). Therefore, by combining Eqs. (53) and (54) with those evolution Eqs. (11)–(18), phase-field simulation of microstructure evolution in systems with compounds described as “Type II” sublattice models can be performed.

4. Test case: binary alloys To demonstrate the applicability of the model in connection with established CALPHAD databases, we choose typical test cases in different binary model alloys. The technically important alloys Fe–C and Ni–Al are chosen. In Fe–C alloys the ferrite and austenite phases are described by the “Type I” sublattice model, while in Ni–Al alloys ordering of the c0 phase in comparison to the disordered c phases is described by the “Type II” sublattice model. 4.1. Fe–C alloys Steels with low carbon concentration solidify in a peritectic transformation, and the solidification of steel Fe–C in binary approximation has been deserving numerous experimental and theoretical investigation [24,39]. The peritectic reaction in the Fe–C alloys involves the motion of a triple junction among a liquid (denoted as “Liquid” or “L”) and two solid phases, i.e., high-temperature ferrite (d) and austenite (c). The calculated phase diagram for the Fe–C binary system in Fe-rich region according to the thermodynamic database developed by Gustafson [40] is shown in Fig. 2(a), from which it is known that the peritectic reaction among L, d and c occurs at 1767.8 K. According to Gustafson [40], liquid phase is described as the

L. Zhang et al. / Acta Materialia 88 (2015) 156–169

substitutional solution phase, in which the one-to-one relation between the variable in CALPHAD database (i.e., the mole fraction of C, xLC ) and the variable in the phase-field model (i.e., cLC ) exists via Eq. (3) by assuming constant molar volume. d and c phases are described by the “Type I” sublattice models ðFeÞ0:25 ðC; VaÞ0:75 and ðFeÞ0:5 ðC; VaÞ0:5 , respectively. The one-to-one relation between the site fractions in CALPHAD database (i.e., y vC and y vVa ) and the phase concentrations (i.e., cdC or ccC ) can be established on the basis of Eq. (33). With the assessed thermodynamic parameters by Gustafson [40], the molar Gibbs energy for both d and c phases can be constructed via Eqs. (19) and then their chemical free energy density. Based on Eq. (39), the diffusion potential of C in d or c phase can be computed. With the known phase concentrations from the phase-filed simulation at each iteration, the corresponding mole fraction, site fractions, free energy densities, as well as the diffusion potentials can be calculated, and these values are utilized for a new iteration in the phase-field simulation. Before two-dimensional (2-D) simulation of the peritectic reaction, a benchmark test in a one-dimensional (1-D) system needs to be performed first. This 1-D system initially consists of a d with concentration of Fe–0.00325C (in mole fraction) and a liquid phase with concentration of Fe– 0.033C. Between the d and liquid phases, a c phase exists. The entire simulation domain is 20 lm with grid spacing as 0.1 lm. The initial thicknesses of the d; c and liquid phases are set to be 8, 5 and 7 lm, respectively. The left and right boundaries for the phase fields, the overall concentration, and the phase concentrations are all set as insulation conditions. The simulation is performed under isothermal condition at 1760 K, which is lower than the peritectic temperature by 7.8 K. The concentration-dependent diffusivities calculated from the mobility database by Liu et al. [41] are employed in the present simulation, while the diffusivity in liquid phase is assumed to be a constant, 2 1  105 cms [42]. The other thermophysical parameters used in the present phase-field simulation are: rdc ¼ rcL ¼ 2:04  104 cmJ 2 for interface energy [24], 4 mdc ¼ mcL ¼ 2:0  102 cm for physical interface mobility Js 3 [4], and P dc ¼ P cL ¼ 1:0  104 cm for interface permeability Js [22]. The initial phase concentrations as well as the average concentration are superimposed on the calculated Fe–C phase diagram (see Fig. 2(a)), from which one may anticipate that the phase concentration of the three phases will approach to the equilibrium ones, and the final state of this 1-D system will consist of equilibrium c and liquid phases. The simulations are started with a fixed distribution of carbon concentration as indicated in Fig. 2(b). The system starts to equilibrate at constant temperature of 1760 K by simultaneously adjusting the c/liquid and c/d interface. The c phase starts to grow gradually, resulting in the outward movement of both d=c and c=L interfaces. As displayed in Fig. 2(b), the d=c interface always moves towards bulk d phase till the d phase disappears fully. In contrast, the c/liquid interface moves first towards the bulk c phase, and then moves back to bulk liquid phase. Moreover, the d=c interface moves much faster than the c=L interface, indicating that the c phase grows at the major expense of the d phase before the disappearance of d phase. After the disappearance of d phase, the fractions of c and liquid phases will adjust to the equilibrium ones for the average concentration in consistency with the phase dia-

163

gram. At 0.2 s, the final equilibrium state is achieved, i.e., that the c and liquid phases with respective equilibrium concentrations coexist according to the equilibrium phase fractions. Next, we perform 2-D phase-filed simulation of the peritectic reaction in Fe–C alloys by coupling to the same CALPHAD database [40]. Considering that the free energy constructed from the CALPHAD database usually covers the entire concentration and temperature range, it is convenient to vary the concentration and temperature by keeping the same CALPHAD database. The initial setting for 2-D system is schematically shown in Fig. 3. All the boundaries for the phase fields as well as the left and right boundaries for the concentrations are set as insulation conditions, while the upper and lower boundaries for the concentrations are set as symmetric conditions. The d and liquid phases initially coexist with a planar d=L interface. A semicircular c phase is inserted between the d and liquid phases. The simulation temperature is chosen as 1758 K, which is close to, but slightly different from 1760 K in 1-D simulation. Again, the concentration-dependent diffusivities calculated from the mobility database by Liu et al. [41] are employed in the 2-D simulation, while the diffusivity in 2 liquid phase is assumed to be a constant, 1  105 cms [42]. All the other thermophysical parameters are kept as the previous values in 1-D simulation because the simulation temperature is just 2 K lower than 1760 K. In addition, we also assume that the new interface, d=L, in the 2-D settings has the same thermophysical parameters as the c/liquid and c/d interface for simplification. Then, we have rdL ¼ rdc ¼ rcL ¼ 2:04  104 cmJ 2 for interface energy, 4 mdL ¼ mdc ¼ mcL ¼ 2:0  102 cm for physical interface Js 3 mobility, and P dL ¼ P dc ¼ P cL ¼ 1:0  104 cm for interface Js permeability. The evolution of the concentration field of C in this Fe–C alloy according to the present 2-D phasefield simulation is presented in Fig. 4. As shown in the figure, the c phase grows gradually along the d=L interface. Due to the solute partitioning from the difference of diffusion potentials over the interface, the concentration in each phase adjust rapidly. Furthermore, it is also observed that the boundaries close to the triple junction become curved due to the curvature contribution in the phase-field equation. The phase-field contours ranging from 0.9 to 0.5 for each phase at time ¼ 9:50 s are shown in Fig. 5(a). The contour with / ¼ 0:5 is represented by the solid lines, by which a hatched area is surrounded. This hatched area is actually the triple junction region, in which the phase fields in each phase are smaller than 0.5. The evolution of the contours with / ¼ 0:5 as well as the trip junctions due to the present phase-field simulation are shown in Fig. 5(b). 4.2. Ni–Al alloys Ni-based super alloys, mainly consisting of disordered c matrix phase and ordered c0 precipitate phase, usually possess excellent mechanical properties, and are widely used as high-temperature structural materials in aerospace and other fields. As the basis for multicomponent Ni-based superalloys, Ni–Al alloys are chosen as the example for validating the present coupling technique in phases described by “Type II” sublattice model. In Ni–Al alloys, c is usually modeled as the substitutional phase, while c0 is described by a “Type II” two-sublattice model (Al,Ni)0.25(Al,Ni)0.75.

164

L. Zhang et al. / Acta Materialia 88 (2015) 156–169

1850

δ

1750

1700

1767.8K

T=1760K

γ

Average concentration Initial concentration of γ phase Initial concentration of liquid phase Initial concentration of δ phase (a)

1650 0

Mole fraction C

Temperature, K

1800

Fe-0.00325C Fe-0.033C

γ

δ

0.03

Liquid

Liquid

Fe-0.0125C

T = 1760 K Time=0.2s 0.02 Time=0.05s Time=0.025s Time=0.003s Time=0s

γ/L interface

0.01 δ/γ interface (b)

0 0

0.01 0.02 0.03 0.04 0.05 Mole fraction C

4

8 12 16 Distance (μm)

20

Fig. 2. (a) Calculated phase diagram of the Fe–C binary system in Fe-rich region according to the thermodynamic database by Gustafson [40]. The initial phase concentrations as well as the average concentration for a three-phase d=c=L Fe–C diffusion couple at 1760 K are superimposed; (b) Evolution of the mixture concentration of C in the three-phase d=c=L diffusion couple annealed at 1760 K according to the present phase-field simulation.

Liquid

D=3μm Wy=8μm Δy=0.1μm

Fe00078C

γLinterface

Fe00256C

γ

δγinterface δγinterface

δ

Fe00045C

Wx=20μm Δx=0.1μm

y x

Fig. 3. Schematic illustration of the initial settings for the 2-D simulation of the peritectic reaction in the Fe–C system.

Here, the thermodynamic parameters are taken from Du and Clavaguera [43], from which the molar Gibbs energies for both c and c0 phases can be then constructed. Based on Eq. (42) or (43), the diffusion potential of Al in c0 phase can be computed. As described above, there is no direct one-toone relation between the variable in CALPHAD database and the variable in the phase-field model. Thus, there is a

Time 008ms

need to employ the two evolution equations ((53) and (54)) for internal equilibrium, from which the one-to-one relation can be established. After that, the phase-field simulation can be started. A simple 1-D phase-field simulation is made in binary Ni–Al alloys. The 1-D system initially consists of a c0 with concentration of Ni–0.25Al (in mole fraction) and a c phase with concentration of Ni–0.10Al. The entire simulation domain is 8 lm with grid spacing as 0.1 lm. The initial thicknesses of the c0 and c phases are set to be 5 and 3 lm, respectively. The left and right boundaries for the phase fields, the concentration and the phase concentrations are all set as insulation conditions, i.e., no flux conditions. The simulation is performed under isothermal condition at 1473 K. At this temperature, the equilibrium site fractions for c0 calculated from the thermodynamic database by Du and Clavaguera [43] are: y uAl ¼ 0:8686 and y vAl ¼ 0:0044. For c phase, the concentration-dependent diffusivities calculated from the mobility database from Zhang et al. [44] are employed in the present simulation. While the diffusivity in c0 phase is taken as the average value over the entire homogeneity range of c0 phase, 2 2:4  109 cms , as done in [44] for phase-field simulation of diffusion couples in the Ni–Al system. The other

Time 150ms

Liquid

γ

Liquid

γ δ

δ

Time 400ms

Time 950ms

Liquid

γ

Liquid

γ δ

0005

001

0015

δ 002

0025 003 0005

001

0015

002

0025

003

Fig. 4. Evolution of concentration field of C in the peritectic reaction of the Fe–C system according to the present 2-D phase-field simulation.

φγ=φL=0.5 φL=0.6 φL=0.7 φL=0.8 φL=0.9

L. Zhang et al. / Acta Materialia 88 (2015) 156–169

y (μm)

Liquid

4.6 4.2

7

γ

3.8

δ

3.4

φδ= φγ =0.5 φδ=0.6 φδ=0.7 φδ=0.8 φδ=0.9

3.6 4.0 4.4 4.8 Spatial coordinate, x (μm)

Time=1.50ms

2

Liquid

Time=0.08ms

γ δ

3 Time=4.00ms Time=9.50ms

(a)

3.0 3.2

φ L =φ δ =φ γ =0.5

6

φδ=φL=0.5 5 φδ=0.6 φδ=0.7 φδ=0.8 φδ=0.9 4

Spatial coordinate,

Spatial coordinate, y (μm)

5.0

165

1 0

1 2 3 4 5 6 Spatial coordinate, x (μm)

(b)

7

Fig. 5. (a) Contour lines (dashed lines) of phase fields for each phase near the triple junction in the peritectic reaction of the Fe–C alloys at 9.50 ms according to the present phase-field simulation. The interval of the contour lines is D/ ¼ 0:1. The hatched region indicates the region where / 6 0:5 for each phase. The solid lines represent the level 0.5 contour lines of the phase fields; (b) Evolution of level 0.5 contours of phase fields in the peritectic reaction of the Fe–C alloys according to the present phase-field simulation.

thermophysical parameters used in the present phase-field 0 simulation are: rcc ¼ 2:0  104 cmJ 2 for interface energy, 8 cm4 cc0 m 0 ¼ 8:7  10 J s for physical interface mobility, and 3 P cc ¼ 10 cm for interface permeability. Js During the simulation, the iteration of phase field and phase concentrations is made first by solving Eqs. (11) and (18), resulting in a new set of phase field and phase concentrations. With the new phase concentration for c0 , the iteration of site fractions are evaluated by means of Eqs. (53) and (54). Then, the new set of site fractions corresponding to the new phase concentration can be obtained. Here, the relaxation coefficient M is taken to be equal to the 3 interface permeability P, i.e., M ¼ 10 cm which ensures Js smooth convergence of the set of equations. The initial phase concentrations as well as the average concentration are superimposed on the calculated Ni–Al phase diagram (see Fig. 6(a)), from which one may anticipate that the phase concentrations of the two phases will approach to the equilibrium ones, and the final state of this 1-D system will consist of equilibrium c and c0 phases. The evolution of the concentration in this 1-D system is presented in Fig. 6(b) according to the phase-field simulation. As can be seen in the figure, the concentrations of the two phases gradually adjust to the respective equilibrium values. Meanwhile, the interface between c and c0 phases moves towards the c0 phase. At 400 s, the final equilibrium state is achieved, i.e., that the c and c0 phases with respective equilibrium concentrations coexist according to the equilibrium phase fractions, which are in consistency with the phase diagram. One much more important feature is presented in Fig. 7, in which the evolution of the overall concentration, and the corresponding phase concentrations of c0 as well as their site fractions due to the phase-field simulation over the interface are given. In addition, the calculated site fractions corresponding to the phase-field simulated phase concentrations of c0 from the CALPHAD database are also appended for direct comparison. As can be seen in the figure, the phase-field simulated site fractions are exactly the same as the those calculated ones from the CALPHAD database. This validates the present coupling

technique for “Type II” sublattice models. Of note is that the relaxation coefficient M set here is large enough, and can guarantee that the exchange between the sublattices achieves the internal equilibrium states rapidly. If a smaller value is given to M, the internal equilibrium may not be achieved so rapidly, and the obtained site fractions will deviate from the equilibrium ones. This kind of nonequilibrium case will subject to our future study.

5. Extension to multicomponent alloys 5.1. Coupling technique in multicomponent alloys Most of the technical alloys are multicomponent, such as steel, Ni-based superalloys, light Al and Mg alloys, and so on. In multicomponent alloys, the situation is much more complex, but the coupling techniques demonstrated in Section 3 are also applicable. In fact, the phases described by the sublattice models in multi-component alloys can be roughly classified into three types, i.e., “Type I”, “Type II” and the mixed type. For “Type I” sublattice models in multicomponent alloys, such as ðA; B; C; . . .Þk ðAÞl , ðA; B; C; . . .Þk ðD; VaÞl , etc., the direct one-to-one relation between the phase concentrations and the site fractions can be directly obtained by Eqs. (3) and (22). Among different kinds of “Type I” sublattice models, ðA; B; C; . . .Þk ðD; VaÞl with D as the interstitial atom, i.e., C or N, is the most representative one. It is the prototype of austenite or ferrite in steels, for instance, ðFe; MnÞk ðC; VaÞl . As for “Type II” sublattice models in multicomponent alloys, like ðA; B; C; . . .Þk ðA; B; C; . . .Þl , one needs the evolution equations similar to (53) and (54) for each solute, from which the relation between the phase concentrations and the site fractions can be obtained. This kind of sublattice model is commonly used to describe the ordered phase in Ni-based superalloys. A more complex case lies in that more than two sublattices are introduced in the phases. For example, the four-sublattice (4SL) model is frequently

166

L. Zhang et al. / Acta Materialia 88 (2015) 156–169 2000 Liquid

1800

Time=10s

0.22 Mole fraction Al

Temperature, K

Initial concentration of γ phase

β

1600 T=1473K

1400

Average concentration

1200

γ

Initial concentration of γ' phase

0.14

γ'

(a)

1000 0 Ni

0.18

Time=400s Time=100s Time=50s

0.1 0.2 0.3 Mole fraction Al

0.4

0.10 0

Time=1s Time=0s

Ni-0.25Al Ni-0.10Al

γ

γ'

T = 1473 K (b)

2

4 6 Distance (μm)

8

Fig. 6. (a) Calculated phase diagram of the Ni–Al system in Ni-rich region according to Du and Clavaguera [43]. The configuration of the diffusion couple Ni–0.25Al/Ni–0.10Al at 1473 K is also superimposed; (b) evolution of the mixture concentration of Al in the diffusion couple Ni–0.25Al/Ni– 0.10Al annealed at 1473 K according to the present phase-field simulation.

Mole fraction / Site fraction Al

1.0 0.8

Time = 1s

Time = 10s

y Alu

Time = 50s

y Alu

Time = 400s

y Alu

Site fraction Al

0.6

Phase-field simulation CALPHAD calculation

0.4 0.2

Time = 100s

y Alu

y Alu

c Alγ' v m

c Alγ' v m

c Alγ' v m

c Alγ' v m

c Alγ' v m

yv

yv

yv

yv

yv

Al Al Al Al Al 0 4.5 5.0 4.0 4.5 3.5 3.0 3.5 3.0 3.5 Distance (μm) Distance (μm) Distance (μm) Distance (μm) Distance (μm)

Fig. 7. Evolution of the phase concentration of c0 as well as its site fractions due to the phase-field simulation. The calculated site fractions corresponding to the phase-field simulated phase concentrations of c0 from the CALPHAD database are also appended for direct comparison.

utilized to model the ordered c0 phase nowadays instead of the two-sublattice (2SL) model due to its superiority in physical meanings. For the 4SL model, e.g., ðA; B; C; . . .Þk ðA; B; C; . . .Þl ðA; B; C; . . .Þm ðA; B; C; . . .Þn , four sub-RVs corresponding to the sublattices should be divided in the phases, as done in Fig. 1. Following the same derivation process from Eqs. (45)–(54) described in Section 3.2, the evolution equations similar to (53) and (54) can be then obtained for accounting for the internal relaxation between different sublattices. In addition to “Type I” and “Type II” sublattice models as in binary alloys, there exists another kind of phases described by the so-called mixed type sublattice models in multicomponent alloys. In fact, most of the intermetallic compounds may fall into this type because the model for compound depends on its crystal structure. For instance, ðA; BÞk ðCÞl ðA; DÞm ðA; BÞn is one of such mixed type sublattice models. For this mixed type sublattice model, the direct one-to-one relation between the phase concentration and the site fraction for species D can be established using equations similar to (3) and (22), as in “Type I” sublattice models. However, the direct one-to-one relation between the the phase concentration and the site fraction for species B does not exist. One needs the evolution equations similar to (53) and (54) to evaluate the relation between the the phase concentration and the site fraction, as done in “Type II” sublattice models.

5.2. Case study on nucleation growth of VC carbide in Fe–V– C alloys Carbides, as the important constituents of the microstructure in steels, play a positive effect on the grain size in ferritic or dual phase steels. This kind of effect comes from two factors. One is the pinning of carbides on grain boundaries, while the other is that carbides act as nucleation sites for ferrite during cooling, and thereby promote fine grains. Vanadium is a common microalloying element in steel. Fine VC monocarbides tend to precipitate from austenite matrix in steels containing vanadium after prolonged tempering, and have a significant effect on the grain boundary mobility and yield strength. We investigate here the growth process of a single vanadium carbide from the austenite matrix after an initial nucleus is set with zero fraction. In fact, during the nucleation process, VC monocarbides share the same lattice structure (B1) and orientation (cube-on-cube) with the surrounding austenite matrix. Thus, the simulation of this special type of coherent nucleation relies significantly on the thermodynamic properties. Relevant differences for a simulation are the alloying composition and resulting lattice parameters. The thermodynamic descriptions of austenite (c) and VC carbide phases are directly taken from the work by Huang [45]. According to Huang, austenite and VC carbide are described as one fcc phase by using a 2SL model

L. Zhang et al. / Acta Materialia 88 (2015) 156–169 2000

167

-5.0 Liquid Liquid+VC

1600 γ

T = 1173 K

G

1200 α+γ+VC

1000

-6.0 -6.5 -7.0

γ,VC

4

γ+VC initial concentration

1400

-1

-5.5 in 10 J mol

Temperature in Kelvin

1800

-7.5 -8.0

α+VC+C

(a)

800 0

-8.5

0.05 0.10 0.15 0.20 Mole fraction of vanadium, carbon

(b) 0

0.1 0.2 0.3 0.4 0.5 Mole fraction of vanadium, carbon

Fig. 8. Calculated (a) quasi binary vertical section along xC ¼ xV ¼ 0:5, and (b) Gibbs energy for fcc phase (both austenite and VC phases) at 1173 K according to the thermodynamic descriptions by Huang [45]. Here, L denotes liquid, c austenite, and a ferrite.

7RWDOUDGLXVRISUHFLSLWDWHLQQP

 yV - Vanadium

\&&DUERQ

γ 9&

QP

(a) 





    (b)  











7LPHLQV

Fig. 9. (a) A snapshot of the simulated concentration fields (here, site fraction is used) of C and V at 6 s; (b) The radius of VC carbide as function of time according the phase-field simulation.

0.95

−6.2

r = 12

interface interface

9

0.94

0.93

t = 6s

−6.4

0.85

carbon

0.84

vanadium

Sitefraction carbon

Sitefraction vanadium

r=5

Diffusion potential of C 10 J m

-3

0.86

−6.6 −6.8

VC

γ

γ

t = 28s

−7

−7.2 −7.4

(a)

0.92 0.83 4 9 14 19 24 Total radius of precipitate in nm

−7.6 0

(b)

32

64 96 X coordinate in nm

128

Fig. 10. (a) Calculated size-dependent site fractions of vanadium and carbon in the center of the VC carbide; (b) The evolution of diffusion potentials of vanadium.

ðFe; V Þ1 ðC; VaÞ1 , which belongs to the “Type I” sublattice model. The coupling technique has been demonstrated in Section 5.1, and can be utilized here. A quasi binary vertical section along xC ¼ xV ¼ 0:5 is calculated on the basis of thermodynamic parameters from Huang [45], and presented in Fig. 8(a). The initial alloy composition for austenite matrix is chosen as Fe–0.05V–0.05C (in mole fraction), and the simulation temperature is fixed at 1173 K, as superimposed in Fig. 8(a). The calculated Gibbs energy for fcc phase (both austenite and VC) is shown in Fig. 8(b). As

can be seen, a metastable miscibility gap of fcc phase should exist, which can be nicely treated in our developed phase-field model with finite interface dissipation [23]. The other thermophysical parameters are: rc;VC ¼ 0:4 mJ2 4 [46] for isotropic interface energy, mc;VC ¼ 1015 mJ s for c 27 VC 2 1 interface mobility, M V ¼ M V ¼ 6:0  10 m s for atomic mobility of vanadium, M cC ¼ M VC C ¼ 1:5 1026 m2 s1 for atomic mobility of carbon, and 5 m3 P cVC ¼ P cVC for interface permeability. Here, only V C ¼ 10 Js the main diffusivities are considered, and the cross ones are

168

L. Zhang et al. / Acta Materialia 88 (2015) 156–169

omitted. In addition, the diffusivities of C in both austenite and VC phases have been lowered by several orders to be not much larger than those of V in the two phases. With this treatment, the present phase-field simulation will be accelerated without significantly changing the simulation outcome. That is because the growth of the VC carbide is limited by the diffusion of vanadium. The simulation domain is set to be 128  128  128 with grid spacing as 1 nm. The interface width is fixed as 5 nm. For the beginning, the matrix is fully occupied by austenite, and a VC carbide nucleus is planted in the center of the cubic domain. A snapshot of the simulated concentration fields for C and V in both matrix austenite and precipitate VC at 6 s is visualised in Fig. 9(a). Here, site fractions of V and C are used to represent concentration. The site fraction of C is displayed on the left side of VC precipitate, while that of V is on the right side. From such a side-by-side plot, one can see that a noticeable depletion exists on the site fraction profile of V ahead of the precipitate due to the lower diffusivity of V. The radius of VC carbide as a function of time according to the phase-field simulation is shown in Fig. 9(b). As can be seen, the VC carbide experiences growth until it begins to establish an equilibrium, with no significant growth after 30 s. As the spherical phase grows to approach the chemical equilibrium according to equations (11) and (15), the only size-dependent parameter in Eq. (11), the capillary contribution I a , is responsible for a deviation between the simulated equilibrium and the pure chemical equilibrium. As the capillary contribution is dependent on the inverse of the total radius, a growing particle approaches the chemical equilibrium depending on the increasing radius. The results are size-dependent site fractions, if assuming that the phase mobility is low enough to neglect kinetic effects on the composition. The result is shown in Fig. 10(a), where the relationship between size and composition is clearly visible, as the site fractions of vanadium and carbon are plotted in the midside node over the value of the radius in the growth phase. In the growth stage below a total radius of 5 nm, the interface with the austenite matrix is not yet fully established, which limits the effective curvature value. The carbide phase enriches with V and C. After a total radius of 5 nm, the size dependence governs the alloy composition. Notably, the site fraction of vanadium increases after a radius greater than 13 nm in contrast to the site fraction of carbon. A kinetic effect must be excluded as explanation for the sinking site fraction of vanadium between 5 nm and 13 nm. The vanadium site fraction in the matrix around the precipitate is not decreasing as an effect of insufficient diffusion transport of vanadium in the matrix. The carbon site fraction is decreasing too, even after a radius larger than 13 nm. The atomic mobility of carbon, which decreases constantly with an increasing radius, is 10 times greater than the atomic mobility of vanadium and the initial site fraction in austenite as well as the permeability is the same for vanadium and carbon. In addition, simulations with a lower interface mobility lab show the same behaviour of the site fractions. Another consequence of the formation of an equilibrium is the convergence of the diffusion potentials, provided that no strain is present [47]. This effect is induced by the chemical diffusion governed by Eq. (18). The change of the diffusion potentials over time is shown in Fig. 10(b) for the example of vanadium. The result of the depletion of the

austenite phase and the change of the diffusion potential is clearly visible too.

6. Conclusions The phase-field model with finite interface dissipation, the sublattice model in CALPHAD formalism, as well as different coupling techniques of the phase-field model and CALPHAD databases available in the literature are briefly summarized. Based on the mesoscopic phase-field model with finite interface dissipation a novel approach to direct coupling of the phase-field model and the sublattice model in CALPHAD formalism is developed. In binary alloys, the sublattice models can be classified into “Type I” and “Type II”. For “Type I” sublattice model, a direct one-to-one relation between the element site fraction in the CALPHAD database and the phase concentration in the phase-field model is found to exist. For “Type II” sublattice models the relation between site and mole fractions is established via an internal relaxation between different sublattices. With the internal relaxation, the energy and potential information from the sublattice model in CALPHAD formalism can be directly used in the phase-field simulation. The present coupling technique in binary alloys is validated in the Fe–C alloys and Ni–Al alloys. Furthermore, the extension of the present coupling techniques into multicomponent alloys is also made, and successfully applied to simulate the nucleation of VC from the austenite matrix in a target Fe–V–C alloy. Acknowledgements Lijun Zhang acknowledges the financial support from the Alexander von Humboldt Foundation, Bonn, Germany, the National Natural Science Foundation of China (Grant Nos. 51301208, 51474239), and the Shenghua Scholar Program of Central South University, Changsha, P.R. China. Yong Du acknowledges the financial support from the National Basic Research Program of China (Grant Nos. 2011CB610401, 2014CB644002). Bo Sundman is grateful to the Alexander von Humboldt Foundation for a senior research grant and the distinguished professor program released by ministry of education of China and the State Administration of Foreign Experts Affairs of China. This work is also supported by Sino-German cooperation group “Microstructure in Al alloys” (Grant No. GZ755), Sino-German Center for Promotion of Science, National Nature Science Foundation of China and German Science Foundation (DFG) and the collaborative research centre SFB-TR 103 “Superalloy Single Crystals” projects C5 and C6. The authors would like to thank Dr. Qing Chen from Thermo-Calc Software AB in Stockholm, Sweden for the discussion on the sublattice model, and Dr. Oleg Shchyglo from ICAMS, Ruhr-Universita¨t Bochum for helpful discussion on this paper.

References [1] W.J. Boettinger, J.A. Warren, C. Beckermann, A. Karma, Annu. Rev. Mater. Res. 32 (2002) 163. [2] L.Q. Chen, Annu. Rev. Mater. Res. 32 (2002) 113. [3] N. Moelans, B. Blanpain, P. Wollants, CALPHAD 32 (2008) 268. [4] I. Steinbach, Modell. Simul. Mater. Sci. Eng. 17 (2009) 073001. [5] Y. Wang, J. Li, Acta Mater. 58 (2010) 1212. [6] I. Steinbach, Annu. Rev. Mater. Res. 43 (2013) 89.

L. Zhang et al. / Acta Materialia 88 (2015) 156–169

[7] I. Steinbach, B. Bo¨ttger, J. Eiken, N. Warnken, S.G. Fries, J. Phase Equilib. Diffus. 28 (2007) 101. [8] T. Kitashima, Philos. Mag. 88 (2008) 1615. [9] S.G. Fries, B. Boettger, J. Eiken, I. Steinbach, Int. J. Mater. Res. (2009) 128. [10] M. Ode, T. Abe, H. Murakami, Y. Yamabe-Mitarai, H. Onodera, J. Metals 59 (2007) 209. [11] T. Wang, G. Sheng, C. Shen, Z.K. Liu, L.Q. Chen, Acta Mater. 56 (2008) 5544. [12] K. Wu, Y.A. Chang, Y. Wang, Scr. Mater. 50 (2004) 1145. [13] K. Wu, N. Zhou, X. Pan, J.E. Morral, Y. Wang, Acta Mater. 56 (2008) 3854. [14] U. Grafe, B. Bo¨ttger, J. Tiaden, S.G. Fries, Scr. Mater. 42 (2000) 1179. [15] J. Eiken, B. Bo¨ttger, I. Steinbach, Phys. Rev. E 73 (2006) 066122. [16] Q. Chen, N. Ma, K. Wu, Y. Wang, Scr. Mater. 50 (2004) 471. [17] M. Hillert, J. Alloys Compd. 320 (2001) 161. [18] M. Hillert, Phase Equilibria, Phase Diagrams and Phase Transformations: Their Thermodynamic Basis, second ed., Cambridge University Press, New York, 2008. [19] E. Borukhovich, P.S. Engels, T. Bo¨hlke, O. Shchyglo, I. Steinbach, Modell. Simul. Mater. Sci. Eng. 22 (2014) 034008. [20] J.Z. Zhu, Z.K. Liu, V. Vaithyanathan, L.Q. Chen, Scr. Mater. 46 (2002) 401. [21] T. Kitashima, H. Harada, Acta Mater. 57 (2009) 2020. [22] I. Steinbach, L. Zhang, M. Plapp, Acta Mater. 60 (2012) 2689. [23] L. Zhang, I. Steinbach, Acta Mater. 60 (2012) 2702. [24] J. Tiaden, B. Nestler, H.J. Diepers, I. Steinbach, Physica D 115 (1998) 73. [25] A. Karma, Phys. Rev. Lett. 87 (2001) 115701. [26] S.G. Kim, W.T. Kim, T. Suzuki, Phys. Rev. E 60 (1999) 7186. [27] L. Zhang, E.V. Danilova, I. Steinbach, D. Medvedev, P.K. Galenko, Acta Mater. 61 (2013) 4155.

169

[28] J.O. Andersson, T. Helander, L. Ho¨glund, P. Shi, B. Sundman, CALPHAD 26 (2002) 273. [29] I. Steinbach, F. Pezzolla, Physica D 134 (1999) 385. [30] U. Preiss, E. Borukhovich, N. Alemayehu, I. Steinbach, F. LaMantina, Model. Simul. Mater. Sci. Eng. 21 (2013) 074006. [31] I. Steinbach, X. Song, A. Hartmaier, Philos. Mag. 90 (2010) 485. [32] R.D. Kamachali, J. Hua, I. Steinbach, A. Hartmaier, Int. J. Mater. Res. 22 (2010) 1332. [33] H. Wang, F. Liu, G.J. Ehlen, D.M. Herlach, Acta Mater. 61 (2013) 2617. ˚ gren, J. Appl. Phys. 72 (1992) 1350. [34] J.O. Andersson, J. A [35] L. Zhang, J. Wang, Y. Du, R. Hu, P. Nash, X.G. Lu, C. Jiang, Acta Mater. 57 (2009) 5324. [36] H. Lukas, S.G. Fries, B. Sundman, Computational Thermodynamics – The Calphad Method, Cambridge University Press, 2007. [37] http://www.thermocalc.se (assessed 16.03.2014). [38] M. Hillert, Physica B 103 (1981) 31. [39] M. Ohno, K. Matsuura, Acta Mater. 58 (2010) 6134. [40] P. Gustafson, J. Scand. Metall. 14 (1985) 289. [41] Y. Liu, L. Zhang, Y. Du, D. Liang, CALPHAD 33 (2009) 614. [42] L. Zhang, Y. Du, I. Steinbach, Q. Chen, B. Huang, Acta Mater. 58 (2010) 3664. [43] Y. Du, N. Clavaguera, J. Alloys Compd. 237 (1996) 20. [44] L. Zhang, I. Steinbach, Y. Du, Int. J. Mater. Res. (formerly Z. Metallkd.) 102 (2011) 371. [45] W. Huang, Z. Metallkd. 82 (1991) 391. [46] Z.G. Yang, M. Enomoto, Metall. Mater. Trans. A 32 (2001) 267. [47] R.D. Kamachali, E. Borukhovich, O. Shchyglo, I. Steinbach, Philos. Mag. Lett. 93 (2013) 680.