A phase field model for isothermal solidification of multicomponent alloys

Acta mater. 49 (2001) 3295–3307 www.elsevier.com/locate/actamat

A PHASE FIELD MODEL FOR ISOTHERMAL SOLIDIFICATION OF MULTICOMPONENT ALLOYS P. -R. CHA†, D. -H. YEON and J. -K. YOON Seoul National University, School of Materials Science and Engineering, Seoul 151-742, South Korea ( Received 16 February 2001; received in revised form 2 May 2001; accepted 4 May 2001 )

Abstract—With the ability to model the kinetics and the pattern formation for solidification, a phase field model has been studied by many scientists. Currently available models, however, are restricted not only to binary alloys but also to those with substitutional solute elements. In this work, a new phase field model is developed to study solidification of a multicomponent alloy containing substitutional as well as interstitial solute elements. By employing the number of moles per unit volume as the concentration variable, the evolution equations of both the phase field and the concentration fields are derived from the free energy functional in the thermodynamically consistent way. In the model, the interfacial region is assumed to be a mixture of solid and liquid with the same composition, but with different chemical potentials. Based on this assumption, the phase field parameters are matched to the alloy properties and an interface thickness limitation is also deduced. Using the chemical rate theory, the phase field mobility is determined at a thin-interface limit condition under the assumption of negligible diffusivity in the solid phase. Another advantage of the model is that any thermodynamic database available in the literature can be directly ported to the model such that quantitative results for solidification of the real alloy systems could be made. As an example, a dendritic growth in an Fe-Mn-C ternary alloy is examined with the thermodynamic data from the commercial software Thermo-Calc code.  2001 Published by Elsevier Science Ltd on behalf of Acta Materialia Inc. Keywords: Solidification; Microstructure; Theory and modeling (kinetics, diffusion); Multicomponent alloys; Phase field model

1. INTRODUCTION

The phase field model has attracted many scientists due to its ability in describing the complex pattern formation in phase transitions such as dendrites. It is a useful method for realistically simulating microstructural evolution involving diffusion, coarsening of dendrites and the curvature and kinetic effects on the moving solid–liquid interface. It is efficient, especially in numerical treatment, because all the governing equations are written in a unified form without distinguishing the interface from the solid or from the liquid phase. In the model, the phase field, f(x,t), characterizes the physical state of the system at each position and time: f = 1 for the solid, f = 0 for the liquid, and 0⬍f⬍1 at the interface. The phase field variable, f, changes steeply but smoothly at the solid–liquid interface region, which avoids direct tracking of the interface position. Therefore, the model can be regarded as a type of a diffuse interface model, which assumes that the interface has a finite

† To whom all correspondence should be addressed. Fax: +82-2-888-5284. E-mail address: [email protected] (P. -R. Cha)

thickness and that the physical properties of the system vary smoothly through the interface. A diffuse interface model was first developed by Van der Waals [1], who considered fluid density as an order parameter. Thereafter, by the mid of 1980’s, the diffuse interface model was applied to the equilibrium properties of the interface [3], antiphase boundary migration by curvature [4], and to the second order transition [2] but not to the first order transition. Langer [5] proposed that the diffuse interface model could be applied to solidification phenomena. By using a singular perturbation method, Caginalp [6] proved that the phase field model could be reduced to the Stefan problem in the limit that the thickness of the interface approaches zero. Kobayashi [7] studied the dendritic growth of a pure melt and Wheeler, McFadden and Boettinger (hereafter abbreviated by WBM) [8] proposed the phase field model for isothermal solidification of a binary alloy. Warren and Boettinger [9] investigated the dendritic growth of a binary alloy with WBM model [8]. The phase field models mentioned above suffer from two limitations which severely restrict their range of application. The first limitation is that the models cannot simulate the case where the kinetic

1359-6454/01/$20.00  2001 Published by Elsevier Science Ltd on behalf of Acta Materialia Inc. PII: S 1 3 5 9 - 6 4 5 4 ( 0 1 ) 0 0 1 8 4 - 7

3296

CHA et al.: MODEL FOR ISOTHERMAL SOLIDIFICATION OF MULTICOMPONENT ALLOYS

undercooling is negligibly small compared to the curvature undercooling (local equilibrium condition) [10]. The second is that the temperature variation in the finite interface region is negligibly small compared to the interface kinetic undercooling [10]. This severely restricts the size of the calculation domain and makes the reliable simulation of dendritic growth possible only at a large undercooling [11]. Karma and Rappel [10, 12] relieved these limitations and showed that it is feasible to determine the parameters in the phase field model at the thin interface limit (finite interface thickness condition) for the solidification of a pure melt using the concept that the temperature gradient in the thin interface region could be linearly approximated. Following the work of Karma and Rappel, Kim et al. [13, 14] demonstrated that the parameters in the phase field model for solidification of binary alloys can be determined at the thin interface limit by linearly approximating chemical potentials in the thin interface region. In studying alloy solidification, however, there still remains a major problem with WBM or other modified models discussed above [8, 9, 13–16]. The problem is that all the existing models are restricted to binary alloys with substitutional solutes. For an alloy with interstitial solutes, these models are not applicable because a mole fraction as a concentration variable can cause a serious error. Furthermore, a model should be able to treat multicomponent alloys, as most commercial alloys are multicomponent. There were several models [17, 18] which can describe the multicomponent alloy solidification. But they were limited to one-dimensional problems (for example, DICTRA code) and cannot describe pattern formation, solute trapping, solute drag, and paraequilibrium phenomena in multicomponent alloy solidification. Therefore, a new phase field model is needed for studying the solidification of multicomponent alloy systems with both substitutional and interstitial elements. This work is organized as follows. In Section 2, a phase field model is developed for multicomponent alloys with both substitutional and interstitial alloy elements. The parameters in the phase field equation are then matched with alloy physical properties, and the restriction on the solid–liquid interface thickness is deduced in Section 3. In Section 4, the phase field mobility is determined at the thin interface limit where the inability to deal with local equilibrium and the restriction on the calculation domain are relieved. In multicomponent alloy solidification, no general relationships exist among temperature, interface velocity and concentration [19]. Therefore, the phase field mobility is determined using the chemical rate theory. In Section 5, numerical calculations for one-dimensional solidification of both Fe-Mn-C and Fe-Cr-Ni alloy systems and for the two-dimensional dendritic growth of Fe-Mn-C alloy and Fe-C alloy systems are presented.

2. PHASE FIELD MODEL TO MULTICOMPONENT ALLOYS

Currently, most phase field approaches for alloy solidification [8, 9, 13–16] use a mole fraction as the concentration variable and are limited to the solidification of binary alloys. Thus, they are not useful for an alloy system with interstitial solutes, due to the difference in diffusion mechanism between an interstitial and a substitutional species (see Section 5.1 for detailed discussion). By employing a new concentration variable, a new phase field model is here presented for studying the solidification of a multicomponent alloy system containing both interstitial and substitutional alloying elements. Consider an isothermal solution with n components, which may exist as either a solid or a liquid in a domain, ⍀. Let mi be the chemical potential of the ith component, which is equal to the partial molar free energy, ¯fi. We now choose the number of moles of the ith component in a unit volume, ci, as the concentration variable. For mathematical simplicity, however, the following assumptions are made: (i) the difference in the molar volume between the solid and liquid phase is neglected, (ii) both the solvent and all substitutional alloying elements have an identical partial molar volume, ⌬, (iii) no interstitial element contributes to the alloy volume, that is, the partial molar volume of any interstitial alloying element is zero, and (iv) the amount of substitutional vacancies is also negligible. As the number of moles of the substitutional and interstitial lattice sites per unit volume is conserved respectively, the new concentration variable provides the following relationship:

冘

1 ci ⫽ , ⌬ i苸S

(1)

and

冘

ci ⫹ cIE ⫽ const.

(2)

i苸I

where the symbols, S and I, indicate a substitutional and an interstitial element, respectively, and cIE is the mole number of empty interstitial sites per unit volume. If the solvent is designated as the nth element, the concentrations of both the solvent and the empty interstitial sites are conveniently removed, leaving the number of the independent concentration variables at n⫺1. The analysis begins with the free energy functional given by Wheeler, Boettinger and McFadden [8]. The total free energy of the system is defined as F[f,c1,…,cn⫺1] ⫽ ⫹

2

冕

f(f,c1,…,cn⫺1)

⍀

e 兩ⵜf兩2d⍀ 2

(3)

CHA et al.: MODEL FOR ISOTHERMAL SOLIDIFICATION OF MULTICOMPONENT ALLOYS

where f is the free energy density of the system, e represents gradient corrections to the free energy density, and f is a phase field ranging from zero in the liquid to one in the solid. It is required that f(1,c1,%,cn⫺1) and f(0,c1,%,cn⫺1) are the thermodynamic free energy densities for a homogeneous solid and a homogeneous liquid phase, respectively. In order to derive a kinetic model, it is assumed that the system evolves in time so that its total free energy decreases monotonically. Since concentration is the conserved field variable and phase field variable is not, the composition field is determined by solving the Cahn–Hilliard diffusion equation [20], whereas the phase field variable is determined by solving the Ginzburg–Landau equation [4]. This case corresponds to Model C of Hohenberg and Halperin [21]. Thus, the evolution equation for phase field becomes: dF ∂f ⫽ ⫺M ∂t df

where f S and f L are the free energy densities of the homogeneous solid and liquid phases, respectively. g(f) is a double-well potential in the form of f2(1⫺ f)2 and w is its height. The function h(f) is equal to f3(6f2⫺15f + 10) and represents a normalized area under the potential, g(f). With equation (8), the derivative term, ff, of equation (6) becomes ff ⫽ h⬘(f)[f s(c1,…,cn⫺1)⫺f L(c1,…,cn⫺1)] ⫹ wg⬘(f)

(9)

where h⬘(f) and g⬘(f) represent dh(f)/df and dg(f)/df, respectively. 3. PARAMETERS FOR THE PHASE FIELD MODEL

冘

∂ck dF Mki(f,c1,…,cn⫺1)ⵜ ⫽ ⵜ· ∂t dci i⫽1 n⫺1

(5)

where M is the phase-field mobility and Mki represents a component of the mobility matrix of solute ˜ . Equation (5) can be derived directly from atoms, M the classical linear irreversible thermodynamics [22]. From equations (4) and (5), the governing equations of the system become: ∂f ⫽ M(e2ⵜ2f⫺ff) ∂t

(6)

and

冘

(8)

(4)

while the concentration of the kth component should evolve with:

冘冘

∂ck ⫽ ⵜ· Mkiⵜfci ⫽ ⵜ· Mkifcicjⵜcj ∂t i⫽1 j ⫽ 1i ⫽ 1 n⫺1

f(f,c1,…,cn⫺1) ⫽ h(f)f S(c1,…,cn⫺1) ⫹ [1 ⫺h(f)]f L(c1,…,cn⫺1) ⫹ wg(f)

3297

In this section, we will find the composition and phase field profiles at equilibrium and the relationships between material properties and phase field parameters (e and w) of the present phase field model. At equilibrium, both the phase field and the concentrations become independent of time. Therefore, for a one-dimensional solidification problem, equations (6) and (7) yield: e2

d dx

d2f0(x) ⫽ ff dx2

冉冘

冊

n⫺1

dfci Mki ⫽0 dx i⫽1

(10)

(11)

Solutions of equations (10) and (11) should provide the concentration c0k(x) and the phase field f0(x) at the thermodynamic equilibrium state. The solution of the diffusion equation (11) in a closed system then becomes

n⫺1 n⫺1

冘

(7)

fci(f0(x),c01(x),…,c0n⫺1(x)) ⫽ const. ⫽ f eci

(12)

n⫺1

⫹ ⵜ·

Mkifcifⵜf

i⫽1

where ff, fci, fcif, and fcicj are ∂f/∂f, ∂f/∂ci, ∂2f/∂ci∂f, and ∂2f/∂ci∂cj, respectively. Mki is determined from

冘

n⫺1

the diffusivity matrix, Dkj =

Mkifcicj. Both equa-

i=1

tions (6) and (7) can be also derived in the thermodynamically consistent way from the entropy functional of Wang et al. [23] and Penrose and Fife [24]. In this work, the free energy density is defined as follows:

where f eci is the chemical potential of the ith component at equilibrium. Equation (12) can be rewritten as

冘

n⫺1

fcif0df0 ⫹

fcicjdc0j ⫽ 0

(13)

j⫽1

If the free energy density, f(f,c1,…,cn⫺1) of equation (8) is given for an alloy system, the equilibrium concentration profile can be expressed as a function of f0(x): c0k(x) ⫽ c˜k(f0)

(14)

3298

CHA et al.: MODEL FOR ISOTHERMAL SOLIDIFICATION OF MULTICOMPONENT ALLOYS

For the determination of c0k(x), the appropriate at f0 = 1 and boundary conditions are: c0k = cS,e k L,e S,e at f = 0, where c and c are the equilibc0k = cL,e 0 k k k rium compositions of the liquid and solid phase, respectively. Let ⌶ be the free energy density function relative to the tie line for the free energy curves at equilibrium,

冘

冕

(21)

0

The interface thickness can also be examined with equation (17). For example, if the interface is defined as the region over from f0 = 0.1 to f0 = 0.9, its thickness is equal to e

2l ⫽

n⫺1

⌶(f,c1,…,cn⫺1) ⫽ f(f,c1,…,cn⫺1)⫺

1

s ⫽ √2e √⌶(f0)⫺⌶(0)df0

e ci i

f c (15)

√2

冕

0.9

0.1

df0 √⌶(f0)⫺⌶(0)

(22)

i⫽1

Then, ff(f0,c01,…,c0n⫺1) = d⌶(f0,c01,…,c0n⫺1)/df because ⌶ci(f0,c01,…,c0n⫺1) = 0 by equation (12). In terms of ⌶, equation (10) becomes e2

d2f0 d⌶(f0,c˜ 1,…,c˜ n⫺1) ⫽ dx2 df0

(16)

For a given s and l, the two phase field parameters, e and w (via ⌶), can be ascertained backward through a numerical integration of equations (21) and (22) and Newton–Raphson method. Combination of equations (21) and (22) and use of wg(f)ⱖ0 provide the following limit for the interface thickness:

Integrating equation (16) after multiplying both sides by df0/dx results in √2 df0 ⫽ ⫺ √⌶(f0)⫺⌶(0) dx e

(17)

where the negative sign is to take account for the decrease in f0(x) with an increase in x in the onedimensional solidification bar (with the liquid phase on the right end). Integrating equation (17) yields f0(x) when the relevant thermodynamic data for an alloy system and the values of e and w are given. The function, ⌶(f0)⫺⌶(0) may be expressed as the sum of two potentials: ⌶(f0)⫺⌶(0) ⫽ wg(f0) ⫹ [Z(f0)⫺Z(0)]

(18)

where Z(f0) is equal to Z(f0) ⫽ h(f0)f S(c˜ 1(f0),…,c˜ n⫺1(f0)) ⫹ [1 ⫺h(f0)]f L(c˜ 1(f0),…,c˜ n⫺1(f0))

冘

(19)

n⫺1

⫺

f ecic˜ i(f0)

i⫽1

The function, Z(f0)⫺Z(0), also has a double-well potential form with two minima at f0 = 0 and f0 = 1 and imposes a restriction on the interface thickness, as will be shown later. Since the solid/liquid interface energy at equilibrium is given by [13, 16]

冕冉 冊

s ⫽ e2

⬁

⫺⬁

df0 2 dx dx

use of equation (17) provides

2lⱕ

s 2

冕

0.9

[1/√Z(f0)⫺Z(0)]df0

0.1

冕

1

(23) √Z(f0)⫺Z(0)df0

0

This limiting condition on the interface thickness can cause a serious problem in the choice of a grid spacing in a numerical application, especially for a two- or a three-dimensional solidification problem. For example, consider an Fe-1.0 at% Mn-0.5 at% C system. Use of equation (23) with the free energy functions for Fe-Mn-C alloy [37] and the thermodynamic data from Thermo-Calc which uses the same free energy functions provides 2lⱕ0.085 µm at 1780 K. For a reasonable accuracy, four or more grids are typically required across an interface, and thus with a grid spacing of about 0.02 µm, a domain of only 20×20 µm2 can be covered even with 106 grids, a practical upper limit for a phase field computation. The grid size problem can get worse at a low undercooling where the curvature of the dendrite tip would become small and thus a large computation cell would be required. If a sufficiently thin interface with a small l is selected while keeping a fixed value for s, w should become very large, while e gets very small. In such a case, the second term, Z(f0)⫺Z(0), in equation (18) may be neglected compared with the first potential well term, and the solution of equation (16) finds

冉 冊

1 √w (24) f0(x) ⫽ [1⫺tanh x ] 2 e√2 The interface thickness and energy then become

(20)

2l ⫽ 2.2√2

and

e √w

(25)

CHA et al.: MODEL FOR ISOTHERMAL SOLIDIFICATION OF MULTICOMPONENT ALLOYS

s⫽

e√w

(26)

3√2

which are identical to those obtained for pure materials in the sharp interface limit [25]. 4. THIN INTERFACE LIMIT

In this section, the phase field mobility, M in equation (4), is examined at the thin interface limit for the solidification of a multicomponent alloy. Here, the “thin interface limit” indicates a situation where the thickness of the moving solid/liquid interface is much smaller than that for the diffusion boundary layer in front of the interface, but still has a finite value. The assumption allows the phase field model to treat a large computation. With the determination of the phase field mobility at the thin interface limit, the two major limitations described in Introduction (the stringent restriction on the interface thickness and the inability to deal with a local equilibrium condition) can be alleviated. Karma and Rappel [10, 12] showed that the phase field mobility can be determined at the thin interface thickness condition for the solidification of pure materials with the idea that the temperature gradient in the thin interface region could be approximated linearly. Inspired by the work of Karma and Rappel, Kim et al. [13, 14] demonstrated that the phase field mobility could also be determined at a thin interface thickness condition for solidification of a binary alloy: they used the idea that the chemical potential can also be linearly approximated in the thin interface region. An approach similar to the work of Kim et al. is followed in this work. The chemical rate theory for solidification [19, 26– 28] describes interface movement by assuming thermally activated hopping process for individual atoms across the interface, and thus the growth velocity is expressed as

冉 冊

V ⫽ V0(T)·[1⫺exp

⌬G ] RT

冘

ci) for the solidifi-

i苸I

cation. In the low velocity limit, i.e. V
V0, equation (27) becomes, by a Taylor series expansion, V⬵⫺V0

⌬G ⌬Fs/ctot ⫽ ⫺V0 RT RT

relationship (dilute binary solution) or temperature– interface velocity relationship (pure materials) [29] and these second response functions have been proved experimentally [30]. If a growth process is “diffusion limited”, V0 is often approximated as DL/di, where DL is the diffusivity of the solute in liquid and di is the thickness of the interface. Many alloy solidifications are known as a diffusion-limited growth, but in general there are many uncertainties for the determination of V0 [28]. The driving force for the solidification without the solute-drag effect in a multicomponent system is given by s L ⌬FS ⫽ f S(cS1,…,cn⫺1 )⫺f L(cL1 ,…,cn⫺1 ) n⫺1 L ∂f (cjL⫺cSj ) L ⫹ ∂cj j⫽1

冘

(29)

With a steady-state interfacial velocity, V, the governing equations for phase field and concentration may be written as ⫺

V df ⫽ e2ⵜ2f⫺ff M dx

(30)

and d dck ⫽ ⫺V dx dx

冉冘

n⫺1

冊

d Mki fci dx i⫽1

(31)

Integrating equation (31), in the solid phase side,

冘

n⫺1

d MSki fci ⫹ A dx i⫽1

⫺VcSk ⫽

(32)

where A is integral constant and in the liquid phase side,

冘

n⫺1

(27)

where V0(T) is the limiting velocity at an infinite driving force and ⌬G is the driving force per mole ( = ⌬FS/ctot, where ctot is the total number of moles per unit volume, equal 1/⌬ +

3299

(28)

In the case of dilute binary solution or pure material, equation (28) can be reduced to the existing temperature–interface velocity–concentration

d MLki fci ⫹ A dx i⫽1

⫺VckL ⫽

(33)

Subtracting equation (33) from equation (32) and using the definition for the mobility matrix:

冘冉

n⫺1

V(ckL⫺cSk) ⫽

j⫽1

冊

dcSj dcLj DSki ⫺DLkj ⫽ JLk ⫺JSk dx dx

(34) Equation (34) simply displays mass balance, i.e. conservation of a solute element at the moving interface. When the diffusion in solid phase is negligible compared with that in liquid phase, the phase field mobility at the thin interface limit can be found from s RTctot ⫽ V0 Me2

(35)

3300

CHA et al.: MODEL FOR ISOTHERMAL SOLIDIFICATION OF MULTICOMPONENT ALLOYS

⫹

e √2

冕 冘冘

carbon atoms per unit volume. The diffusion equation of the presented model is then reduced to

0.1 n⫺1 n⫺1

(⌳)jifcif0bj(f0)df0

0.9 j ⫽ 1i ⫽ 1

where ⌳ is the inverse matrix of fcicj, and βj(f0) is a complicated integral over from f = f0 to f = 0.1. The detail of the derivation for equation (35) is given in Appendix A. In the sharp interface limit of l→0, the second term of the right hand side of equation (35) vanishes and equation (35) becomes s RTctot ⫽ V0 Me2

∂ ∂cC ∂cC ⫽ D ∂t ∂x ∂x

where D is the diffusivity of carbon in Fe. The wellknown error function solution is

冉 冊

cC ⫽ cs⫺cserf (36)

In the case of a dilute binary alloy system, it can be easily shown by using the relationship between V0 and the interface kinetic coefficient [29] that the above relationships (35) and (36) reduce to the existing relationships between the interface kinetic coefficient and the phase field mobility in a thin interface limit condition [13] and in a sharp interface limit condition [8], respectively. 5. RESULTS AND DISCUSSION

As noted before, there can arise a difficulty when a mole fraction is used as a concentration variable for an alloy with interstitial elements. Through an elementary diffusion problem, the difficulty is examined and its implication on a phase field model is discussed. Within the limit of one-dimensional isothermal solidification, the phase field approach is then shown to predict reasonably well the equilibrium properties and the kinetics for a multicomponent alloy system. Finally, some preliminary results for the morphological evolution of a two-dimensional dendrite are presented for both a binary and a ternary alloy with carbon as interstitial element. 5.1. Difficulty with a mole fraction as a concentration variable for an interstitial alloy A mole fraction is a concentration variable widely used in many fields. It can, however, cause a significant error in a phase field model for solidification of alloys with interstitial solutes. If substitutional, the diffusion of a solute requires the back diffusion of the solvent and the number of moles per unit volume is conserved at 1/⌬. If interstitial, however, the back diffusion of the solvent is neither required nor the number of moles is conserved. Therefore, care should be exercised for a phase field model involving interstitial diffusion with a mole fraction as the concentration variable. As an example, let us consider a carburization process of a steel plate. The surface concentration at x = 0 is maintained at a constant value, while carbon atoms diffuse from the surface into the steel. The boundary conditions are cC(at x = 0) = cs and cC(at x = ⬁) = 0, where cC is the number of moles for

(37)

x

2√Dt

(38)

Now, in terms of the mole fraction for carbon atoms, XC, the diffusion equation of some phase field models has the same form as equation (37), and is equal to ∂XC ∂ ∂XC ⫽ D ∂t ∂x ∂x

(39)

The boundary conditions appropriate for equation (39) are XC(at x = 0) = Xs = cs/(cFe + cs) and XC(at x = ⬁) = 0, and its solution is XC ⫽ Xs⫺Xserf

冉 冊 x

2√Dt

(40)

If carbon atoms were substitutional, cFe + cC would be a constant and equations (38) and (40) would become identical to each other. However, cFe + csolute depends on cC which is a function of position, and equations (38) and (40) must be different. A relative error associated with equation (40) may be expressed as (vC⫺XC)/vC = Xserf(x/2√Dt) where XC of equation (40) and vC is the correct mole fraction, cC/(cFe + cC), evaluated with cC from equation (38). Unfortunately, the relative error increases with x. At a typical diffusion distance of x⬇√Dt, the error is equal to one half of the surface mole fraction, 苲Xs/2, for example, yielding about 2% when Xs = 0.04. The difference, (vC⫺XC)/Xs, is displayed for a few values of Xs in Fig. 1, where the peak error is shown to occur at x⬇√Dt. Imagine that a phase field model is applied for a solidification problem with a mole fraction for an interstitial component. The diffusion field around the solid–liquid interface would then involve such an error, which would in turn alter the concentration profile of the solute element affecting incorrectly the kinetics of solidification. So far, a binary system is considered. Certainly, if a care is not taken, a serious error is easily conceivable when a ternary or a higherorder alloy system is studied [31]. 5.2. Equilibrium concentrations in one-dimensional isothermal solidification To examine equilibrium properties at solidification, two alloy systems are considered: one is a substi-

CHA et al.: MODEL FOR ISOTHERMAL SOLIDIFICATION OF MULTICOMPONENT ALLOYS

3301

system and all the ancillary thermodynamic data are taken from Thermo-Calc, which uses the same free energy functions as the before-mentioned free energy functions. The small thickness limit of the Fe-Mn-C calls for a grid size smaller than that for the Fe-CrNi system. Thus 2l is set at 60 nm for the Fe-Cr-Ni system and at 6 nm for the Fe-Mn-C system, respectively. Grid sizes of both 10 nm and 1 nm are used such that an interfacial region could be covered with six grid spacings. The total grid number in a system is fixed at 200. The parameters of the phase field equation, e and w were obtained from equations (21) and (22). With the assumption of a diffusion-limited growth, V0 = VD = 1 m/s. The diagonal terms of the solute diffusivity matrix are assumed to be the same for each phase and are taken to be 3×10⫺9 m2/s. The off-diagonal terms are set at zero for convenience, although the diffusivity matrix of each component for each phase could be obtained in a general way [33, 34]. Equations (6) and (7) are discretized by a secondorder differencing scheme for the spatial derivatives and a simple explicit Euler scheme for the time derivative. The no-flux boundary condition, Fig. 1. Absolute error, (vC⫺XC)/Xs versus diffusion depth in the analysis of a carburization process.

tutional ternary Fe-9.0 at% Cr-10.0 at% Ni system and the other is an interstitial ternary Fe-1.0 at% Mn0.5 at% C alloy system. The Fe-Cr- Ni system solidifies to austenite at 1757 K with an undercooling from the liquidus temperature of 3.52 K, whereas the FeMn-C system solidifies to d-ferrite at 1780 K with an undercooling from the liquidus temperature of 17.65 K. The equilibrium volume fractions and concentrations of both solid and liquid phases for each of the two alloys are listed in Table 1. When the interface energy, s, is set at 0.4 J/m2, equation (23) yields 2lⱕ0.973 µm for the Fe-Cr-Ni system at 1757 K and 2lⱕ0.085 µm for the Fe-MnC system at 1780 K: the thermodynamic free energy functions of each phase are obtained from [35, 36] for the Fe-Cr-Ni system and [37] for the Fe-Mn-C Table 1. Equilibrium volume fractions and concentrations of both solid and liquid phases for Fe-Cr-Ni and Fe-Mn-C alloys [32] System Phase Equilibrium concentration Equilibrium volume fraction System Phase Equilibrium concentration Equilibrium volume fraction

1757 K, Fe-9.0 at% Cr-10.0 at %Ni Liquid 10.03 at% Cr 10.38 at% Ni 0.143

Austenite (FCC) 8.83 at% Cr 9.94 at% Ni 0.857

1780 K, Fe-1.0 at% Mn-0.5 at% C Liquid 1.27 at% Mn 1.43 at% C 0.218

d-ferrite (BCC) 0.923 at% Mn 0.242 at% C 0.782

∂f ∂ci ⫽ ⫽0 ∂n ∂n

(41)

is imposed at both ends of the computation cell, where n indicates the outward normal. Solidification is initiated by placing a small solid seed of the same composition as the liquid at one end of the computation cell. The time advancing step, ⌬t, is set at (⌬x)2/5DLCC, well below the upper numerical instability limit, (⌬x)2/2DLCC. Time-dependent concentration profiles for Cr and Ni in the Fe-Cr-Ni system are displayed in Fig. 2(a) and (b), respectively, where the evolution time from 0 to 10⫺2 set is examined. In the initial stage up to 3.75×10⫺4 s, the solid concentrations are shown to sharply decrease, while the liquid concentrations adjacent to the interface rapidly increase: this transient effect should be, at least in part, attributed to the initial uniform concentration in both the solid seed and the liquid phase. As the solid–liquid interface advances, the solute concentrations in both phases increase, reaching equilibrium at about 10⫺2 s. At equilibrium, the predicted concentrations are 8.829 at% Cr and 9.936 at% Ni in the solid phase with a volume fraction of 0.8564, and 10.031 at% Cr and 10.38 at% Ni in the liquid phase with a volume fraction of 0.1436: they are in excellent agreement with the data from Thermo-Calc (see Table 1). Figure 3(a) and (b) show solute concentration profiles for the Fe-Mn-C alloy system at 1780 K. At 5.88×10⫺5 s, an equilibrium state is reached in the system with 0.933 at% Mn and 0.263 at% C in solid and 1.263 at% Mn and 1.444 at% C in liquid. The equilibrium volume fractions of d-ferrite and liquid are 0.797 and 0.203, respectively. Again, the pre-

3302

CHA et al.: MODEL FOR ISOTHERMAL SOLIDIFICATION OF MULTICOMPONENT ALLOYS

Fig. 2. Concentration profiles for (a) Cr and (b) Ni in one-dimensional Fe-9.0 at% Cr-10.0 at% Ni alloy. The times are 0, 3.75×10⫺4, 7.50×10⫺4, 1.25×10⫺3, 1.625×10⫺3, 2.375×10⫺4, 1.0×10⫺2 s.

Fig. 3. Concentration profiles for (a) Mn and (b) C in one-dimensional Fe-1.0 at% Mn-0.5 at% C alloy. The times are 0, 2.96×10⫺6, 5.92×10⫺6, 1.18×10⫺5, 1.77×10⫺5, 2.35×10⫺5, 5.88×10⫺5 s.

dicted equilibrium properties are in good agreement with the data from Thermo-Calc as shown in Table 1. Both systems show rapid solute rejections out of the solid seeds in their initial stage of evolution. It is interesting to note that the evolving solid phases display solute concentrations well below the equilibrium ones. 5.3. Steady-state, one-dimensional isothermal solidification For a steady state solidification in one-dimension, exact solute concentrations can be obtained with the classical sharp interface model. Thus, at least under a steady state condition, it is possible to test a phase field model by comparing its numerical solutions with the exact values. Consider a system having initially a uniform composition, c⬁ at a temperature below the melting point. Imagine that a solid phase starts to

grow from one end of the system and then develops a steady state because of a solute sink in the liquid phase away from the solid–liquid interface [13, 14]. A dilute solution is also assumed such that all the offdiagonal diffusivities can be neglected. In the classical sharp interface model with diffusion in liquid only, the situation can be described by [13, 14] d2cLi dcLi ⫽ DiiL 2 ⫺V dx dx

(42)

冉 冊

L V(1⫺kei )cL,i i ⫽ ⫺Dii

ciL(x∗) ⫽ c⬁i

dcLi dx

(43)

i

(44)

CHA et al.: MODEL FOR ISOTHERMAL SOLIDIFICATION OF MULTICOMPONENT ALLOYS

where kei is the equilibrium partition coefficient of the ith component, x* is the prescribed distance between is the interface and the solute sink in liquid, and cL,i i the concentration of the ith component at the interface. The exact solution is then given by cLi (x) ⫽ c⬁i ⫹ c⬁i

e i

(1⫺k )(e

⫺Vx/DL ii

⫺e

e i

⫺Vx∗/DL ii

)

⫺Vx∗/DL ii

1⫺(1⫺k )(1⫺e

)

(45)

where the interface is located at x = 0. Equation (45) has the unique value at any position if V and x* are given. The same interstitial alloy with Fe-1.0 at% Mn-0.5 at% C at 1780 K is tested for comparison between the exact solution and the phase field model. The following additional material properties are used: m2/s, DCCS = 6.0×10⫺9, DSMnMn = 3.0×10⫺13 L = 2.0×10⫺8, keMn = 0.72677, DLMnMn = 1.0×10⫺9, Dcc and keC = 0.16923. For the phase field calculation, a long time behavior of its full dynamics is closely examined to determine a steady state velocity, V, which is in turn used for the analytical solution of equation (45). In Fig. 4, the interface concentrations for Mn and C are plotted as a function of x*. The

3303

solid curves stand for the numerical results from the phase field model, while the solid circles indicate the analytical solution. An excellent agreement is shown for the two results with less than 1% error. 5.4. Two dimensional dendritic growth of a ternary alloy The phase field model is here extended to a twodimensional solidification problem. In order to see the effect of anisotropic interfacial energy, an energy form adopted by Kobayashi [7] is used: e ⫽ e¯ h ⫽ e¯ [1 ⫹ gecos(4q)]

(46)

where e¯ and ge are constants and q = tan⫺1(fy/fx). Consequently, the system has a 4-fold symmetry and its governing equations for a two-dimensional dendritic growth become: 1 ∂f ⫽ e¯ 2hh⬘[sin(2q)(fyy⫺fxx) M ∂t ⫹ 2cos(2q)fxy]⫺(1/2)e¯ 2[h⬘2 ⫺hh⬘⬘][2sin(2q)fxy⫺ⵜ2f⫺cos(2q)(fyy ⫺fxx)] ⫹ e¯ 2h2ⵜ2f⫺fφ

(47)

冘

(48)

and ∂ck Mkiⵜfci ⫽ ⵜ· ∂t i⫽1 n⫺1

冘冘

n⫺1 n⫺1

⫽ ⵜ·

j ⫽ 1i ⫽ 1

Fig. 4. Interface concentration versus solute sink distance x*. The top figure stands for the variation of Mn concentration while the bottom represents C concentration in the liquid phase at the interface. The solid curves represent the phase field results and the solid circles are from the analytic solution of equation (45).

冘

n⫺1

Mkifcicjⵜcj ⫹ ⵜ·

Mkifcifⵜf

i⫽1

where h⬘ = dh/dq and h⬘⬘ = d2h/dq2. The interface thickness, 2l, is set at 50 nm with a grid size equal to 10 nm, the interfacial energy, s, is equal to 0.4 J/m2, and e¯ is determined from equations (21) and (22). Other ancillary input data are: ge = 0.04, VD = 10 m/s, DsMnMn = 3.0×10⫺13 m2/s, s L = 6.0×10⫺9, DLMnMn = 1.0×10⫺9, DCC = 2.0×10⫺8, DCC and T = 1775 K. As before, the off-diagonal diffusivity terms are assumed to be zero. During the numerical calculation for equation (47), a stochastic noise is added at the interface region in order to include fluctuations at the interface, which gave rise to the initiation of the morphological instability [9]. A Fortran code is developed and run on Intel Pentium chip incorporated computer. Equiaxed dendrite shapes for an Fe-O.5 at% C binary and an Fe-0.5 at% C-0.005 at% Mn ternary alloy are pictured at the evolution time of 2.125×10⫺5 s in Fig. 5. Taking advantage of the four-fold morphological symmetry, only a quarter of the entire dendrite is chosen as the computational cell, and a solid seed is placed at the comer of the bottom-left, that is at the origin. Therefore, the dendrite pictures in Fig. 5

3304

CHA et al.: MODEL FOR ISOTHERMAL SOLIDIFICATION OF MULTICOMPONENT ALLOYS

ary arms are shown to better develop, that is, more in an orderly form than those of the binary alloy, and slow-diffusing Mn atoms appears responsible as they can get amassed in front of the interface, enhancing the interface instability to a certain degree. The growth patterns show many features consistent with those observed in the work of Warren and Boettinger [9]. The spine of the primary stalk has a low concentration, but the regions between the secondary arms have the highest. There is a tendency for the sidebranches to grow initially not perpendicular to the primary stalk, but then later to grow into the perpendicular direction, Behind the dendrite tip, the roots of secondary arms are usually narrower than their midriffs, and this is demonstrated more clearly in Fig. 6, where solute concentration profiles for the Fe-0.5 at% C-O.005 at% Mn alloy are shown at the evolution time of 2.5×10⫺5 s. These features are commonly observed in real dendrites. A close examination reveals that spines of some

Fig. 5. Carbon concentration profiles in the dendritic growth of (a) an Fe-0.5 at% C and (b) an Fe-0.005 at% Mn-0.5 at% C alloy at 1775 K and at the evolution time of 2.125×10⫺5 s. The calculation was done on a 500×1000 grid. The colors range from blue (lowest concentration) to green (central concentration) to red (highest concentration).

are those obtained by reflecting the calculated results about the y-axis. The periodic boundary conditions along the x- and the y-axis are imposed during the computation. No-flux conditions of equation (41) are used at the right and top boundaries of the computational cell. The undercoolings from the liquidus temperature are 27.65 K for the Fe-C binary and 27.63 K for the Fe-C-Mn ternary alloy. The CPU hours are about 35 for the Fe-C and 78 for the Fe-MnC system on an Intel Pentium III-450 MHz computer. The dendritic growth rate for the binary alloy is found to be greater than that for the ternary alloy: it may be reasoned that addition of a third element lowers the equilibrium liquidus and solidus lines and reduces the degree of the supersaturation and the driving force of the solidification (see Appendix B for detailed discussion). In the ternary alloy, the second-

Fig. 6. (a) Mn concentration and (b) C concentration profiles in the dendritic growth of the Fe-0.005 at% Mn-0.5 at% C alloy at 1775K and at the evolution time of 2.5×10⫺5 s.

CHA et al.: MODEL FOR ISOTHERMAL SOLIDIFICATION OF MULTICOMPONENT ALLOYS

secondary branches show a concentration as low as that along the spine of the primary stalk (see Fig. 6(b)). Such a low concentration along the spine of a primary stalk appears due to capillary effect at the growing tip and it is natural to see such a spine even in a secondary dendritic arm. In the ternary Fe-MnC system, the carbon diffusivity in the solid phase is almost as fast as the diffusivity of Mn in the liquid phase. Consequently, carbon homogenization process is expected to occur even in the solid phase. Thus it is not a surprise to observe that in Fig. 6(a), the Mn concentration along the primary spine, i.e. along the vertical line at the center, remains low all the way to the original seed place, while the carbon concentration along the primary spine in Fig. 6(b) increases as the seed place is approached. In addition, an inset figure at the left-upper corner of Fig. 6(b) shows the close-up of two small liquid pockets near the seed, having nearly a 4-fold symmetry, consistent with the effects of the anisotropic interfacial energy given in equation (46). 6. CONCLUSIONS

A phase field model is presented for solidification of a multicomponent alloy system with a low interface velocity condition and a thin interface limit condition. The model could be applicable to any alloy systems with both substitutional and interstitial elements, and is different from other phase field models, which are restricted for substitutional binary alloys. It is found that treatment of the diffusion field for an interstitial alloy system requires caution in choosing correct concentration variables. One-dimensional numerical calculations for solidification of an Fe-Mn-C and an Fe-Cr-Ni alloy system demonstrate that the model could provide a good description for the equilibrium and kinetic properties for isothermal solidification. It is also demonstrated that the model could be extended to simulate isothermal growth of large-scale dendritic microstructures in a highly supersaturated liquid, and to predict the pattern of the solute distribution both in the solid and in the liquid phase. Acknowledgements—The authors extend their thanks to Prof. Kyu-Hwan Oh at Seoul National University for the help in thermodynamic data evaluation and to Prof. Seong-Gyoon Kim at Kunsan National University and Prof. Jong K. Lee at Michigan Technological University for valuable discussions. They also appreciate greatly the financial support of Brain Korea 21 program supported by Ministry of Education, Korea.

3305

3. Cahn, J. W. and Hilliard, J. E., Acta metall. J. Chem. Phys., 1958, 28, 258. 4. Allen, S. M. and Cahn, J. W., Acta metall., 1979, 27, 1085. 5. Langer, J. S., Directions in Condensed Matter. World Scientific, Singapore, 1986. 6. Caginalp, G., Phys. Rev. A, 1989, 39, 5887. 7. Kobayashi, R., Physica D, 1993, 63, 410. 8. Wheeler, A. A., Boettinger, W. J. and McFadden, G. B., Phys. Rev. A, 1992, 45, 7424. 9. Warren, J. A. and Boettinger, W. J., Acta metall. mater., 1995, 43, 689. 10. Karma, A. and Rappel, W. -J., Phys. Rev. E, 1996, 53, R3017. 11. Wang, S. -L. and Sekerka, R. F., Phys. Rev. E, 1996, 53, 376O. 12. Karma, A. and Rappel, W. -J., Phys. Rev. E, 1998, 57, 4323. 13. Kim, S. G., Kim, W. T. and Suzuki, T., Phys. Rev. E, 1998, 58, 3316. 14. Kim, S. G., Kim, W. T. and Suzuki, T., Phys. Rev. E, 1999, 60, 7186. 15. Tiaden, J., Nestler, B., Diepers, H. J. and Steinbach, I., Physica D, 1998, 115, 73. 16. Wheeler, A. A., Boettinger, W. J. and McFadden, G. B., Phys. Rev. E, 1993, 47, 1893. 17. Crusius, S., Inden, G., Knoop, U., Hoglund, L. and Agren, J., Z. Metallkd., 1992, 83, 673. 18. Lee, B. -J. and Oh, K. H., Z. Metallkd., 1996, 87, 195. 19. Ludwig, A., Physica D, 1998, 124, 271. 20. Cahn, J. W., Acta metall., 1961, 9, 759. 21. Hohenberg, P. C. and Halperin, B. I., Rev. Mod. Phys., 1977, 49, 435. 22. de Groot, S. R. and Mazur, P., Non-equilibrium Thermodynamics. North-Holland Publishing Company, Amsterdam, 1963. 23. Wang, S. -L., Sekerka, R. F., Wheeler, A. A., Murray, B. T., Coriell, S. R., Braun, R. J. and McFadden, G. B., Physica D, 1993, 69, 189. 24. Penrose, O. and Fife, P. C., Physica D, 1990, 43, 44. 25. McFadden, G. B., Wheeler, A. A., Braun, R. J., Coriell, S. R. and Sekerka, R. F., Phys. Rev. E, 1993, 48, 2016. 26. Aziz, M. J. and Kaplan, T., Acta metall., 1988, 36, 2335. 27. Turnbull, D., J. Phys. Chem., 1962, 66, 609. 28. Aziz, M. J. and Boettinger, W. J., Acta metall., 1994, 42, 527. 29. Kurz, W. and Fisher, D., Fundamentals of Solidification. Trans Tech Publications, Aedermannsdorf, Schweiz, 1992. 30. Willnecker, R., Herlach, D. M. and Feuerbacher, B., Phys. Rev. Lett., 1989, 62, 2707. 31. Yeon, D.-H., Cha, P.-R. and Yoon, J.-K., Scripta Materialia, (in press). 32. Thermodynamic equilibrium data for Fe-Cr-Ni and Fe-MnC alloy systems were calculated from thermodynamic software Thermo-Calc. 33. Andersson, J. -O. and Agren, J., J. Appl. Phys., 1992, 72, 1350. 34. Kirkaldy, J. S. and Young, D. J., Diffusion in the Condensed State. Institute of Metals, London, 1987. 35. Hillert, M. and Qiu, C., Metall. Trans. A, 1990, 21A, 1673. 36. Lee, B. -J., J. Kor. Inst. Met. Mater., 1993, 31, 480. 37. Huang, W., Metall. Trans. A, 1990, 21A, 2115. 38. Rappaz, M. and Boetinger, W. J., Acta mater., 1999, 47, 3205. APPENDIX A

REFERENCES 1. Van der Waals, originally published (in Dutch) in Verhandel. Konink. Akad. Weten. Amsterdam (Sect. 1), 1893, 1, 56: Translated by J.S. Rowlinson (in English) in J. Stat. Phys., 1979, 20, 197. 2. Ginzburg, V. L. and Landau, L. D., J. Exptl. Theoret. Phys. (USSR), 1950, 20, 1064.

When the diffusivity in solid is negligible, from equation (32), A = ⫺VcSk. Integrating equation (31) yields,

冘

n⫺1

⫺V(ck⫺cSk) ⫽

d Mki fci dx i⫽1

(A1)

3306

CHA et al.: MODEL FOR ISOTHERMAL SOLIDIFICATION OF MULTICOMPONENT ALLOYS

This can be transformed to d f ⫽ dx ci

冘

a⫽

n⫺1

⫺BikV(ck⫺cSk)

(A2)

k⫽1

where Bik is a component of the inverse matrix of ˜ . Integrating equation (A1) at the interface provides M the chemical potential profile across the interface:

冕 冘

n⫺1

x

fci ⫽ f Sci⫺V

Bik(ck⫺cSk)dx⬘

(A3)

V M

冕冉冊 冕 l

冕

l df 2 d2f df dx ⫽ e2 dx 2 dx ⫺l dx dx l ∂f df dx ⫺ ∂f dx ⫺l

⫺l

(A4)

l

冘

n⫺1

⫺f S(c1,…,cn⫺1)⫺

冘冕 冕冘 cjL,i

n⫺1

V

S,i

cj

j⫽1

(cLi ⫺cSi )f Lci

a⫽

(A5)

⫹

j⫽1

cj

S,i

冕 冘冘

0.1 n⫺1 n⫺1

(⌳)jifcif0bj(f0)df0

0.9 j ⫽ 1i ⫽ 1

where the matrix, ⌳, is the inverse matrix of fcicj and

冕 冘

f0 k ⫽ 1

Bejk(c0k⫺cS,e k ) √⌶(f)⫺⌶(0)

df

(A10)

It should be noted that the left hand side of equation (A7) is the same as the driving force for the solidification per unit volume. Thus, with equations (22) and (27), the phase field mobility in the thin interface limit can be determined from the following equation s RTctot ⫽ V0 Me2

Bjk(ck⫺cSk)dx⬘dcj

xk ⫽ 1

冋冕冉冊

1 ⫽V M

冘冕 冕 冘 cjL,i

Bejk(c0k⫺cS,e k )dx⬘dcj

xk ⫽ 1

e s ⫹ Me2 √2

l n⫺1

冘

n⫺1

n⫺1

S,e

i⫽1

f L(c1,…,cn⫺1)⫺f S(c1,…,cn⫺1)⫺ ⫺c )f

cj

(A8)

(A9)

Therefore, equation (A4) can be written as

L ci

⫺l l n⫺1

df0 2 dx dx

c0k and f0 are the kth composition and phase field at equilibrium, respectively, and are the solutions of equations (11) and (17). a may be modified into more tractable form by using equations (13), (17) and (20) as follows:

⫹

S i

l

0.1 n⫺1

∂f df dx ⫽ f L(c1,…,cn⫺1) ∂f dx ⫺l

⫺

j⫽1

bj(f0) ⫽

By equation (A3), the last term on the right-hand side of equation (A4) becomes

冕

冕冉 冊

冘冕 冕 冘

n⫺1

cjL,e

⫺lk ⫽ 1

Multiplying df/dx on the both side of equation (30) and integrating from x = ⫺l to x = l yields ⫺

⫹

1 M

l

⫺l

i⫽1

(A6)

l n⫺1

Bjk(ck⫺cSk)dx⬘dcj

xk ⫽ 1

册

冘

n⫺1

⫺cSi )f cLi ⫽ aV

(cLi

0.1 n⫺1 n⫺1

(⌳)jifcif0bj(f0)df0

0.9 j ⫽ 1i ⫽ 1

ctot ⫽ 1/⌬ ⫹

冘

cS,e i .

i苸I

APPENDIX B

Now consider a situation when the interface velocity is low. To the first order in the Pe´ clet number, ˜ where D ˜ is average interface diffusivity. P = 2lV/D Equation (A6) then becomes f L(c1,…,cn⫺1)⫺f S(c1,…,cn⫺1)⫺

√2

冕 冘冘

where

(cLi

df 2 dx dx

e

(A11)

(A7)

i⫽1

where a is a constant for a given temperature and is equal to

The supersaturation associated with solute element i L,e S,e i can be defined as (cL,e i ⫺ci)/(ci ⫺ci ) [38], where i ci is an initially given concentration of solute i. If this definition is used, the supersaturation associated with carbon of the Fe-0.5 at% C system is 0.9184, that in the Fe-0.5at% C-0.005at% Mn system is 0.9180 and the supersaturation associated with Mn is also 0.9180: the equilibrium concentrations of each phase in the Fe-C and Fe-C-Mn systems are obtained from the software Thermo-Calc. Hence, by adding small amount of Mn, the supersaturation decreases. Although the difference of supersaturations between the Fe-C and Fe-C-Mn systems is small, the differ-

CHA et al.: MODEL FOR ISOTHERMAL SOLIDIFICATION OF MULTICOMPONENT ALLOYS

ence will increase around the interface because slowdiffusing Mn atoms can enrich in front of the interface. The difference of the growth rates between the binary and ternary alloys can be also explained in terms of the driving force of the solidification. The

3307

driving force for the solidification can be directly calculated from equation (29). The driving force for the Fe-0.5at% C system is 80.68 J/mole and that of the Fe-0.5at% C-0.005at% Mn system is 80.56 J/mole at 1775 K. The driving force for the solidification also decreases by addition of the third element.