Phase-field model for multicomponent alloy solidification

Phase-field model for multicomponent alloy solidification

ARTICLE IN PRESS Journal of Crystal Growth 274 (2005) 281–293 www.elsevier.com/locate/jcrysgro Phase-field model for multicomponent alloy solidificati...

396KB Sizes 2 Downloads 78 Views

ARTICLE IN PRESS

Journal of Crystal Growth 274 (2005) 281–293 www.elsevier.com/locate/jcrysgro

Phase-field model for multicomponent alloy solidification Pil-Ryung Chaa,, Dong-Hee Yeonb, Jong-Kyu Yoonb a

School of Advanced Materials Engineering, Kookmin University, 861-1, Chongnung-dong, Songbuk-gu, Seoul 136-702, Republic of Korea b School of Materials Science and Engineering, Seoul National University, Seoul 151-742, Republic of Korea Received 21 May 2004; accepted 1 October 2004 Communicated by G.B. McFadden Available online 18 November 2004

Abstract A phase-field model is presented for studying the solidification of multicomponent alloys, which can consider the different diffusion mechanisms between the substitutional and the interstitial solutes appropriately. By defining the interfacial region as a mixture of solid and liquid with the same inter-diffusion potential and determining and matching the phase-field parameters to the alloy properties at a thin-interface limit condition, the limit of interface thickness due to the contribution of the chemical free energy density in the interfacial region is remarkably relieved and a large-scale calculation domain is permitted. One-dimensional numerical calculations for the solidification of an dilute Fe–Mn–C alloy system demonstrate that the model could provide a good description for the equilibrium and kinetic properties for isothermal solidification of multicomponent alloys. Two-dimensional computations for the evolution of large-scale dendritic microstructures is performed and solute distributions in both the solid and liquid phase and the dendritic shapes are compared at different temperatures. r 2004 Elsevier B.V. All rights reserved. PACS: 64.70.Dv; 81.30.Fb; 81.10.Aj; 05.70.Ln Keywords: A1. Growth models; A1. Morphological instability; A1. Phase-field model; A1. Solidification; B1. Multicomponent alloys

1. Introduction The phase-field model has attracted many scientists due to its ability in describing the complex pattern formations in phase transitions Corresponding author. Tel.: +82 29104656;

fax: +82 28885284. E-mail address: [email protected] (P.-R. Cha).

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,

0022-0248/$ - see front matter r 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jcrysgro.2004.10.002

ARTICLE IN PRESS 282

P.-R. Cha et al. / Journal of Crystal Growth 274 (2005) 281–293

fð~ r; tÞ; characterizes the physical state of the system at each position and time: f ¼ 1 for the solid, f ¼ 0 for the liquid, and 0ofo1 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 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 1980s, 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 the singular perturbation method [6], Caginalp [7] proved that the phase-field model could be reduced to the Stefan problem in the limit that the thickness of the interface approaches zero and Kobayashi [8] studied the dendritic growth of a pure melt by using the phase-field model for pure materials. The models which were proposed for simulating dendritic growth in pure undercooled melt [8–13] have been extended to modeling of binary alloy solidification [14–23] and recently multicomponent alloy solidification [24]. The phase-field models for alloy solidification may be classified into four groups depending on the definition of free energy density for interfacial region and how they were derived [23]. The first is a model by Wheeler, Boettinger, and McFadden (hereafter abbreviated by WBM) model [14–16,20], which has been used most widely. In WBM model, the interfacial region is assumed to be a mixture of solid and liquid with the same composition, but with different chemical potentials, which leads a certain limit in the interface thickness due to chemical energy contribution to the interfacial energy and this limitation makes the large-scale simulation difficult [20]. Caginalp and Xie [18] have also proposed a similar model. The second is a model by Steinbach and coworkers

[19], which uses a different definition that the interfacial region is assumed to be a mixture of solid and liquid with different compositions, but constant in their ratio. It has less limit in the interface thickness but the model is thermodynamically correct for a dilute alloy only. The third model has been developed by Losert and coworkers [21], and Lo¨wen and coworkers [22], which is based on two unrealistic assumptions: One is that the liquidus and solidus lines in the phase diagram are assumed to be parallel and the other that the solute diffusivity is constant in the whole space of the system. The fourth was proposed by Kim and coworkers [23] which is free from the limit for the interface thickness in the WBM model [14–16] by assuming the interfacial region as a mixture of solid and liquid with the same chemical potential, but with different compositions. In studying alloy solidification, however, there still remains a major problem with the models discussed above [14–16,19–21,23]. 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 [24,25]. It can be said that the application of the phase-field model for solidification to practical technology is closely linked to our ability to model microstructural development in multicomponent alloys (three or more components) because the most alloys for practical use are multicomponent. The existing models [26,27], which are the sharp-interface model and can describe the multicomponent alloy solidification, are applied to only one-dimensional problems (for example, DICTRA code) and hence these models can not describe the pattern formation, solute trapping, solute drag and paraequilibrium phenomena in multicomponent alloy solidification. In addition, there have been a few researches for the solidification of multicomponent alloys (three or more components) [28,29], but most of them have the limit of the interface thickness or can be applied only to substitutional alloy systems. Therefore, it is very important to develop a phase-field model for the solidification of multicomponent alloy systems with both substitutional and interstitial elements.

ARTICLE IN PRESS P.-R. Cha et al. / Journal of Crystal Growth 274 (2005) 281–293

283

Recently, Cha and coworkers [24] presented a phase-field model for the isothermal solidification of multicomponent alloys. However, the model has a severe limit of the interface thickness due to the variation of the chemical free energy in the interfacial region. For example, the interface thickness should be less than 85 nm for an Fe— 1.0 at% Mn—0.5 at% C alloy system [24]. This thickness limit constrains the size of calculation domain and makes it very difficult numerical simulation at a low undercooling, where the curvature of the dendrite tip would become small and thus a large computation domain would be required. In this study, we present a phase-field model for multicomponent alloy solidification which relieves the limit of the interface thickness. For this purpose, 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 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 [30]. Therefore, the phase-field mobility is determined using the chemical rate theory for solidification [30–33,35]. In Section 3, the evolution equations of the phase-field and concentration fields are rewritten in the dilute solution limit and the steady-state solutions of the concentration profiles across the interface are deduced. In Section 4, numerical calculations for one-dimensional solidification of a ternary alloy and for the large-scale two-dimensional dendritic growth of dilute Fe–Mn–C alloy system are presented.

diffusion mechanism between an interstitial and a substitutional species (see [24,25] for detailed discussion). By employing the number of moles in a unit volume as the concentration variable and defining the interfacial region as a mixture of solid and liquid with the same chemical potential but with different compositions, a new phase-field model is here presented for studying the solidification of a multicomponent alloy system containing both interstitial and substitutional alloying elements, which does not suffer from the limit of interface thickness due to the contribution of chemical free energy in the interfacial region [23,24]. Consider a solution with n components, which may exist as either a solid or a liquid in a domain, O: 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) both the solvent and all substitutional alloying elements have an identical partial molar volume, D; (ii) there is no contribution of interstitial element to the alloy volume, (iii) the difference in the molar volume between the solid and liquid phase is neglected, and (iv) the amount of substitutional vacancies are also negligible. As the number of moles of the substitutional and interstitial lattice sites per unit volume is conserved, respectively, the new concentration variables have the P P following relationships: c ¼ 1=D and i2S i i2I ci þ cIE ¼ const:; 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 total free energy of the system is defined as

2. Phase-field model to multicomponent alloys

Fðf; c1 ; . . . ; cn1 Þ   Z 2 ¼ dO f ðf; c1 ; . . . ; cn1 Þ þ jrfj2 ; 2 O

Most phase field approaches for alloy solidification using a mole fraction as the concentration variable [14–23] are not useful for an alloy system with interstitial solutes, due to the difference in

ð1Þ

where f is the free energy density of the system,  represents gradient corrections to the free energy

ARTICLE IN PRESS P.-R. Cha et al. / Journal of Crystal Growth 274 (2005) 281–293

284

density, and f is a phase field ranging from zero in the liquid to one in the solid. In this work, the definition of the interfacial region proposed by Kim et al. [23] is adopted to remove the limitation of the interface thickness due to the contribution of the chemical free energy in the interfacial region. Any point within the interface region is considered as a mixture of solid and liquid with a same inter-diffusion potential. Then, the free energy density and concentration of ith component are expressed as f ðf; c1 ; . . . ; cn1 Þ ¼ wgðfÞ þ hðfÞf S ðcS1 ; . . . ; cSn1 Þ þ ½1  hðfÞ f L ðcL1 ; . . . ; cLn1 Þ

ð2Þ

and ci ¼ hðfÞcSi þ ½1  hðfÞ cLi ;

(3)

with f ScS ðcS1 ; . . . ; cSn1 Þ ¼ f LcL ðcL1 ; . . . ; cLn1 Þ; i

(4)

i

where f S and f L are the free energy densities of the homogeneous solid and liquid phases, respectively, and cSi and cLi are the mole numbers of the ith component in the solid and liquid phases, respectively. f ScS and f LcL are qf S =qcSi and qf L =qcLi ; i i respectively. gðfÞ is a double-well potential, for example, one of the simplest forms is f2 ð1  fÞ2 and w is its height. The interpolation function hðfÞ represents a normalized area under the potential, gðfÞ: In order to derive a kinetic model, it is assumed that the system evolves in time so that its total free energy decreases monotonically. The evolution equation for phase field becomes: 1 qf dF ¼  ¼ r 2 rf  f f M qt df ¼ r 2 rf  wg0 ðfÞ þ h0 ðfÞ " # n1 X L S L S L ðci  ci Þf cL ;

f f  i¼1

i

ð5Þ

while the concentration of the kth component should evolve with: n1 n1 X X qck dF ¼r M ki r ¼r M ki rf ci dci qt i¼1 i¼1

¼r

n1 X n1 X

M ki f ci cj rcj

i¼1 j¼1

þr

n1 X

M ki f ci f rf;

ð6Þ

i¼1

where M is the phase-field mobility and M ki represents a component of the mobility matrix of ~ f f ; f c ; f c f ; and f c c are qf =qf; solute atoms, M: i i i j 2 qf =qci ; q f =qci qf; and q2 f =qci qcj ; respectively. M ki isPdetermined from the diffusivity matrix, 0 0 Dkj ¼ n1 i¼1 M ki f ci cj : h ðfÞ and g ðfÞ represent dhðfÞ=df and dgðfÞ=df; respectively. In order to find the composition (c0k ðxÞ) and phase field (f0 ðxÞ) profiles at equilibrium and the relationships between material properties and phase-field parameters ( and w) of the present phase-field model, a one-dimensional solidification problem is considered at equilibrium (see [24] for detailed discussion). Under boundary conditions of f0 ¼ 1 (solid) at x ! 1 and f0 ¼ 0 (liquid) at x ! þ1; Eqs. (5) and (6) for a one-dimensional equilibrium state yield using the same approach introduced in Ref. [24]:  pffiffiffiffi  w 1 f0 ðxÞ ¼ 1  tanh pffiffiffi x (7) 2 2 and L;e 0 c0i ¼ hðf0 ðxÞÞcS;e i þ ½1  hðf ðxÞÞ ci ;

cL;e i

(8)

cS;e i

where and are the ith concentrations in the liquid and solid phases at equilibrium, respectively. Using the composition and phase-field profiles at equilibrium, direct evaluation of the interface energy, s; and interface thickness, 2l; gives pffiffiffiffi  w s ¼ pffiffiffi ; (9) 3 2 pffiffiffi  2l ¼ a 2 pffiffiffiffi ; w

(10)

ARTICLE IN PRESS P.-R. Cha et al. / Journal of Crystal Growth 274 (2005) 281–293

where a is a constant which is dependent on the definition of the interface thickness, e.g., a ’ 2:2 when f0 changes from 0.1 to 0.9 at loxol; and a ’ 2:94 when f0 changes from 0.05 to 0.95 [24,25]. It should be noted that the parameter relationships (9) and (10) are identical to those obtained for pure materials in the sharp interface limit [16,24]. The phase-field mobility, M in Eq. (5), can be determined using the chemical rate theory for solidification [30–33,35]. Through the same approach introduced in appendix of Ref. [24], the phase-field mobility at the thin-interface limit is determined from the following equation: RTctot s  ¼  pffiffiffiffiffiffi M2 V0 2w L;e S;e S;e

xðcL;e 1 ; . . . ; cn1 ; c1 ; . . . ; cn1 Þ;

ð11Þ

P where ctot ¼ 1=D þ n1 i¼1;i2I ci ; V 0 is the limiting velocity at an infinite driving force and S;e xðcL;e i ; . . . ; ci ; . . .Þ is a temperature-dependent function defined by L;e S;e S;e xðcL;e 1 ; . . . ; cn1 ; c1 ; . . . ; cn1 Þ



n1 X

S;e ðcL;e j  cj Þ

j¼1

Z

0

Z

S;e ðcL;e k  ck Þ

k¼1 0

1

n1 X

f0

Bejk

½1  hðfÞ

dfh0 ðf0 Þ df0 fð1  fÞ

ð12Þ

3. Kinetics of the multicomponent alloy solidification The dilute solution limit is often useful in both engineering applications and finding out fundamentals. Here, it is shown that a steady-state solution of the concentration profile across the diffuse interface can be obtained as a function of the interface velocity and that this solution reduces to the solution obtained in the classical sharpinterface model when the thickness of the moving interface is much smaller than that for the diffusion boundary layer in front of the interface. Throughout this section, the following assumptions are made: (i) the off-diagonal components of diffusivity matrix are neglected, (ii) each diagonal diffusivity term, Dii ; in both the interfacial region and liquid is constant, and (iii) diffusivities in the solid phase are negligible. In the dilute solution limit, the inter-diffusion potential of a phase U is expressed as follows: U fU cU ¼ A i þ i

and is a component of the inverse matrix ~ at the equilibrium state. In the case of a of M dilute binary alloy system, the relationship between V 0 and the interface kinetic coefficient [34] reduces the relationship (11) to the existing relationship between the interface kinetic coefficient and the phase-field mobility [14,15,23]. In the sharp interface limit that the thickness of the moving pffiffiffiffiffiffi solid/liquid interface vanishes, i.e. = 2w ! 0; the phase-field mobility reduces RTctot =V 0 ¼ s=ðM2 Þ; which is only valid for the case that the inter-diffusion potential of each solute atom is constant across the interfacial region.

RT lnðXU i Þ; D

(13)

where AU i is a constant for the ith component U element in phase U: XU i is defined as X i for U U interstitial element and X i =X n for substitutional U element in U phase, where X U i and X n represent the mole fractions of the solute atom and solvent atom, respectively, and X U is defined by i Pn1 U U ci =ð1=D þ i¼1;i2I ci Þ: Eq. (13) produces that f LcL cL with iak is negligible in comparison with i

Bejk

285

k

f LcL cL when both ith and kth components are i

i

substitutional solute elements and that f LcL cL with i

k

iak becomes zero when the ith or the kth component is interstitial. Eq. (4) yields the following relationships for both substitutional and interstitial solute at the limit that all n  1 compositions go to zero, i.e. XU n ’ 1: cSi cS;e i ¼ ¼ kei ; cLi cL;e i

(14)

where kei represents the equilibrium partition coefficient of the ith alloy component.

ARTICLE IN PRESS P.-R. Cha et al. / Journal of Crystal Growth 274 (2005) 281–293

286

In a one-dimensional steady state with the interfacial velocity, V ; the governing equation for concentration is transformed to d L f L dx ci

fcc ¼  i i V ðci  cS;i i Þ; Dii

(15)

where cS;i k indicates the composition of the kth element at the solid side (x ¼ l) of the interface. At the limit where all n  1 compositions vanish, we can get the following approximations: d L dcL RT 1 dcLi f cL ¼ f LcL cL i  ; i i dx dx i D cLi dx

(16)

ci  ½1  ð1  kei ÞhðfÞ cLi

(17)

and f ci ci 

RT : D½1  ð1  kei ÞhðfÞ cLi

(18)

Therefore, Eq. (15) becomes dcLi V L V cS;i i þ : ci ¼ Dii 1  ð1  kei ÞhðfÞ dx Dii

(19)

Under the boundary condition cLi ¼ cSi =kei ¼ e cS;i i =ki at x ¼ l; the solution of this first-order inhomogeneous equation is cLi ¼

cS;i V S;i Vx=Dii i eV ðxþlÞ=Dii þ c e Dii i kei Z x 0 eVx =Dii dx0 :

e l 1  ð1  k i ÞhðfÞ

ð20Þ

For the case that the thickness of the moving interface is much smaller than that for the diffusion boundary layer in front of the interface, it can be easily shown that cLi ðxÞ at xbl converges to the following solution obtained in the classical sharp-interface model with the same diffusivity in liquid as Dii :

1 Vx=Dii S;i S;i L : (21) ci ¼ ci  ci 1  e e ki

From Eqs. (17) and (20), the concentration profile ci ðxÞ across the interface is given by cLi ¼ ½1  ð1  kei ÞhðfÞ

" cS;i V S;i Vx=Dii c e

i e eV ðxþlÞ=Dii þ Dii i ki 3 Zx 0 eVx =Dii dx05:

1  ð1  kei ÞhðfÞ

ð22Þ

l

Eq. (22) shows very interesting phenomenon, which occurs only in multicomponent alloy system. There can arise a regime of solidification, that can be characterized by the absence of partitioning of the slow diffusing element. This type of transformation is called para-equilibrium, which was first introduced by Hultgren in 1947 [38]. Under para-equilibrium, fast diffusing solute atoms diffuse at an appreciable rate, but slow diffusing ones diffuse so slow that the growing phase inherits the same content of slow diffusing solutes from the parent phase. To examine this phenomenon, an ternary alloy system is selected, where the diagonal diffusivity terms and the equilibrium partition coefficient of both solute elements are: DL11 ¼ 2 109 m2 =s; DL22 ¼ 2 108 ; ke1 ¼ 0:5; and ke2 ¼ 0:5: Concentration profiles are calculated from Eq. (22). Because the diffusion rate of the second solute is faster than that of the first one, the first element is expected to be trapped earlier than the second one as the interface velocity increases. Concentration profiles of the solute elements for various interface velocities are displayed in Fig. 1. For increasing interface velocity, both solute elements shows decrease of partitioning and at 1 m/s of the interface velocity, the slow diffusing element is trapped while the fast diffusing one has appreciable partitioning. Thus solidification with more than 1 m/s of the migration rate shows para-equilibrium and if the interface velocity further increases, the fast diffusing solute is expected to be also trapped.

4. Results and discussion In this section, several numerical simulations for solidification of a ternary alloy are presented. Fe—

ARTICLE IN PRESS P.-R. Cha et al. / Journal of Crystal Growth 274 (2005) 281–293

287

Fig. 1. Concentration profiles of the first and second solutes for various interface velocity. The equilibrium partition coefficients of both solutes are the same and set at 0.5 and the diffusivity of the second solute atom is larger by order of one than that of the first one.

0.1 at% Mn—1.0 at% C alloy system, which has not only the substitutional solute atoms but also interstitial solute atoms, is selected for this numerical study. This alloy system solidifies to dferrite under the liquidus temperature 1793.5 K. The equilibrium volume fraction of solid, f V ; and the equilibrium concentrations of both solid and liquid phases are: f V ¼ 0:515; X L;e Mn ¼ 1:172 S;e 2 103 ; X L;e ¼ 1:738

10 ; X ¼ 8:389 104 ; Mn C S;e 3 and X C ¼ 3:091 10 : With the numerical study of one-dimensional isothermal solidification, the phase-field approach is shown to predict reasonably well the equilibrium properties and the kinetics for a multicomponent alloy system. Some preliminary results for the morphological evolution of a two-dimensional dendrite are presented. 4.1. Equilibrium and kinetics in one-dimensional isothermal solidification To examine equilibrium properties and kinetics at the solidification of the selected ternary system, the computer simulation of the one-dimensional isothermal solidification is performed. In this study, the selected ternary alloy system is regarded as a dilute solution and under this assumption the

Eqs. (5) and (6) reduce into: 1 qf q2 f RT X L X S;e n ¼ 2 2 þ h0 ðfÞ ln nS L;e  wg0 ðfÞ; M qt qx D X nXn (23) n1 qf ci qck q X ¼ ; M ki qx i¼1 qt qx

f ci ¼

RT ðhðfÞ ln XSi þ ð1  hðfÞÞ ln XLi Þ: D

(24)

(25)

Eqs. (23) and (24) 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, rn f ¼ rn ci ¼ 0; is imposed at both ends of the computation domain. In order to test the equilibrium properties at solidification, solidification process from a small solid seed is simulated. The solid seed has the same composition as the liquid and is located at one end of the computation domain. The material properties used in simulations are: s ¼ 0:4 J=m2 ; 2l ¼ 600 nm and a ¼ 2:94: For the test of the equilibrium properties, the selection of the input

ARTICLE IN PRESS P.-R. Cha et al. / Journal of Crystal Growth 274 (2005) 281–293

288

parameters related with kinetics are not important. Here, V 0 is set to 0.5 m/s and only diagonal terms of the solute diffusivity matrix are considered and are set to 2 108 m2 =s for both phases. Note that the diffusivity matrix of each component for each phase could be obtained in a general way [36,37]. In computer simulation, the total grid number in the system is fixed at 200 and grid spacing is set to be 100 nm such that an interfacial region could be covered with six grid spacings. The time advancing step, Dt; is set at ðDxÞ2 =5DLCC ; well below the upper numerical instability limit, ðDxÞ2 =2DLCC : Time-dependent concentration profiles for Mn and C are displayed in Fig. 2(a) and (b), respectively. The sharp concentration profile at the initial stage can be explained with transient effect due to the initial uniform concentration distribution over whole computation domain. After solidification proceed, the solute concentrations in both phases reaches equilibrium at about 0.24 s. At equilibrium, the predicted concentrations are 0.0837 at% Mn and 0.3091 at% C in the solid phase with a volume fraction of 0.517, and 0.1170 at% Mn and 1.7385 at% C in the liquid phase: they are in excellent agreement with the given thermodynamic data. In order to check the ability of the model to consider the kinetics of the solidification, the numerical solution by the phase-field model for a steady-state solidification in one-dimension is compared with exact one obtained with the classical sharp interface model. In the case that diffusion occurs in liquid only and there is a solute sink in the liquid phase away from the solid–liquid interface such that the interface moves with constant velocity V ; the concentration profile obtained with the classical sharp-interface model is (24): L

1 cLi ðxÞ ¼ c1 i þ ci

ð1  kei ÞðeVx=Dii  eV x 1  ð1 

kei Þð1





=DLii

 L eV x =Dii Þ

Þ

; (26)

where x is the prescribed distance between the interface and the solute sink in liquid and cL;i c1 i i are the concentration of the ith component at the liquid side of the interface and the concentration

Fig. 2. Concentration profiles for (a) Mn and (b) C in onedimensional Fe—0.1 at% Mn—1.0 at% C alloy. The times are 0; 1 102 ; 2 102 ; 3 102 ; 5 102 ; 0.24 s.

far from the interface, respectively. Eq. (26) has the unique value at any position if V and x are given. The selected alloy system at 1780 K is tested for comparison between the exact solution and the phase-field model. The following additional material properties are used: DSMnMn ¼ 1:0

ARTICLE IN PRESS P.-R. Cha et al. / Journal of Crystal Growth 274 (2005) 281–293

289

1014 m2 =s; DSCC ¼ 1:0 109 ; DLMnMn ¼ e 9 8 L 5:0 10 ; DCC ¼ 2:0 10 ; kMn ¼ 0:716; and keC ¼ 0:178: The interface thickness, 2l is set at 6 nm and grid size of 1 nm is used. The phase-field mobility is assumed to be 0.5. 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 Eq. (26). In Fig. 3, the concentrations for Mn and C at the liquid side of the interface are plotted as a function of x : An excellent agreement is shown for the two results with less than 1% error. 4.2. 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 [8] is used:  ¼ ¯Z ¼ ¯½1 þ ge cosð4yÞ ;

(27)

where ¯ and ge are constants and y ¼ tan1 ðfy =fx Þ: Consequently, the system has a four-fold symmetry and its governing equations for a two-dimensional dendritic growth become 1 qf ¼ ¯ 2 ZZ0 ½sinð2yÞðfyy  fxx Þ þ 2 cosð2yÞfxy

M qt  ð1=2Þ¯2 ½Z02 þ ZZ00

½2 sinð2yÞfxy  r2 f  cosð2yÞ

ðfyy  fxx Þ þ ¯2 Z2 r2 f  f f

ð28Þ

and n1 X qck ¼r M ki rf ci ; qt i¼1 0

(29) 00

2

2

where Z ¼ dZ=dy and Z ¼ d Z=dy : Using the Eqs. (28) and (29) two-dimensional computer simulations were performed for the isothermal solidification of the selected ternary alloy system (Fe—0.1 at% Mn—1.0 at% C alloy system) at 1770 and 1780 K. In computer simulation, this alloy system is regarded as a dilute solution. The interface thickness, 2l; is set at 1200 nm with a grid size equal to 200 nm, the interfacial energy, s; is equal to 0.4 J/m2. Other

Fig. 3. 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 circles represent the phase-field results and the solid curves are from the analytic solution of Eq. (26). The unit of each concentration is mole fraction.

ancillary input data are: ge ¼ 0:06; DSMnMn ¼ 1:0 1014 m2 =s; DSCC ¼ 1:0 109 ; DLMnMn ¼ 9 8 L 5:0 10 ; and DCC ¼ 2:0 10 : As before, the off-diagonal diffusivity terms are assumed to be

ARTICLE IN PRESS 290

P.-R. Cha et al. / Journal of Crystal Growth 274 (2005) 281–293

zero. When the phase-field mobility is determined, it is assumed that V 0 is independent of the temperature and is set to 5.5 m/s. The phase-field mobilities, which are determined by the Eq. (11), are 0.16 at 1770 K and 0.0582 at 1780 K. During the numerical calculation for Eq. (28), 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 [39]. Equiaxed dendrite shape at 1770 K is pictured at the evolution time of 1:25 104 s in Fig. 4. 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 corner of the bottom-left, that is at the origin. Therefore, the dendrite picture in Fig. 4 is that obtained by reflecting the calculated results about the y-axis. No-flux conditions are used at the right and top boundaries of the computational cell. The CPU hours are about 48 on an Intel Pentium IV—2 GHz computer. The growth pattern in Fig. 4 shows many features consistent with those observed in the work of Warren and Boettinger [39]. 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. These features are commonly observed in real dendrites. A close examination reveals that spines of some secondary branches show a concentration as low as that along the spine of the primary stalk. 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 dendrite arm. In the ternary Fe–Mn–C 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. 4(a), the Mn concentration along the primary spine, i.e. along the vertical line at the center,

Fig. 4. (a) Mn concentration and (b) C concentration profiles in the dendritic growth of the Fe—0.1 at% Mn—1.0 at% C alloy at 1770 K and at the evolution time of 1:25 104 s:

remains low all the way to the original seed place, while the carbon concentration along the primary spine in Fig. 4(b) is already homogenized. Fig. 5 shows the equiaxed dendrite shape at 1780 K with that at 1770 K. The spatial distributions of Mn concentration at the same evolution time for 1780 and 1770 K are presented in Fig. 5(a) and (b), respectively. In comparison with the

ARTICLE IN PRESS P.-R. Cha et al. / Journal of Crystal Growth 274 (2005) 281–293

291

Fig. 5. The spatial distributions of Mn concentration at the same evolution time of 1:25 104 s for (a) 1780 K and (b) 1770 K. Later stage of the growth of the dendrite at 1780 K is shown in (c) and (d), which represents the concentration profile of Mn and C at 1:0 103 s; respectively.

dendrite at 1770 K, the velocity of the solid–liquid interface at 1780 K is slow and the tip radius of the dendrite is larger due to the low undercooling. Due to the low velocity of the solid–liquid interface, the partitioning of Mn atoms becomes clearly and this is consistent with the results in Section 3. Later stage of the growth of the dendrite is shown in Fig. 5(c) and (d), which represents the concentration profile of Mn and C at 1:0 103 s; respectively. From Fig. 5(c), the homogenization of Mn atoms is progressed conspicuously in solid phase, which can be explained by partitioning due to low interface velocity and enough time to diffuse in solid phase. The carbon concentration is already homogenized as a matter of course. Until now our numerical simulations were presented only for a dilute ternary alloy. In numerical simulations of non-dilute multicomponent alloys, there are some complexities related to

solve the condition (4) of the equal chemical potentials in solving a set of Eqs. (5), (6), (3), and (4). Contrary to the dilute solution approximation, there is no analytic expression to solve the constraint equations (4) in general non-dilute multicomponent alloys nor even in non-dilute binary alloys. In binary alloy systems, Kim et al. have proposed a useful and efficient numerical method to solve it [23] but it is very difficult to extend their method to general multicomponent alloy systems (for ternary alloys, it is possible to use their method by making the table of compositions cSi and cLi of the solid and liquid phases as functions of the chemical potentials f ScS ¼ f LcL ¼ i i f ci ; but for non-dilute multicomponent alloys with more than three components, this method is not available.). We have used multi-dimensional Newton–Raphson method for nonlinear systems of equations (MDNR) [40] to solve the constraint

ARTICLE IN PRESS 292

P.-R. Cha et al. / Journal of Crystal Growth 274 (2005) 281–293

equations (4). In general, as the free energy function is well defined (positive definite) one of the solute compositions, the convergence rate of MDNR is not too bad. We have simulated the dendritic growth of the same system as shown in Fig. 4 with MDNR. The calculation time was increased by  350%; when compared with the simulation using the dilute solution approximation (with MDNR, it takes one week to calculate the system shown in Fig. 4 on Intel Pentium IV— 2 GHz. Recently, as the computing performance is more advanced, it is acceptable). Although MDNR is time-consuming method, we can consider it as one of the numerical methods to solve the constraint equations (4) of the equal chemical potentials for non-dilute alloys.

5. Conclusions A phase-field model is presented for solidification of a multicomponent alloy system. The model could be applicable to any alloy systems with both substitutional and interstitial elements by using the number of moles as concentration variable and it relieves the limit of interface thickness due to the contribution of the chemical free energy density in the interfacial region which seriously restricts the calculation domain size [24]. One-dimensional numerical calculations for solidification of an Fe–Mn–C alloy system demonstrate that the model could provide a good description for the equilibrium and kinetic properties for isothermal ternary solidification. It is also demonstrated that the model could be extended to simulate the evolution of large-scale dendritic microstructures and to predict the pattern of the solute distribution both in the solid and in the liquid phase.

Acknowledgements This work was supported by the internal research fund 2004 of Kookmin University in Korea.

References [1] Van der Waals, Originally published in Verhandel. Konink. Akad. Weten. Amsterdam (Section 1) 1 (1893) 56 (in Dutch), transl. J.S. Rowlinson, J. Stat. Phys. 20 (1979) 197 (in English). [2] V.L. Ginzburg, L.D. Landau, J. Exp. Theor. Phys. (USSR) 20 (1950) 1064. [3] J.W. Cahn, J.E. Hilliard, J. Chem. Phys. 28 (1958) 258. [4] S.M. Allen, J.W. Cahn, Acta Metall. 27 (1979) 1085. [5] J.S. Langer, in: G. Grinstein, G. Mazenko (Eds.), Directions in Condensed Matter Physics, World Scientific, Singapore, 1986, p. 164. [6] R.E. O’Malley Jr., Singular Perturbation Methods for Ordinary Differential Equations, Springer, New York, 1991. [7] G. Caginalp, Phys. Rev. A 39 (1989) 5887. [8] R. Kobayashi, Physica D 63 (1993) 410. [9] S.-L. Wang, R.F. Sekerka, Phys. Rev. E 53 (1996) 3760. [10] S.-L. Wang, R.F. Sekerka, A.A. Wheeler, B.T. Murray, S.R. Coriell, R.J. Braun, G.B. McFadden, Physica D 69 (1993) 189. [11] O. Penrose, P.C. Fife, Physica D 43 (1990) 44. [12] A. Karma, W.-J. Rappel, Phys. Rev. E 53 (1996) R3017. [13] A. Karma, W.-J. Rappel, Phys. Rev. E 57 (1998) 4323. [14] A.A. Wheeler, W.J. Boettinger, G.B. McFadden, Phys. Rev. A 45 (1992) 7424. [15] A.A. Wheeler, W.J. Boettinger, G.B. McFadden, Phys. Rev. E 47 (1993) 1893. [16] G.B. McFadden, A.A. Wheeler, R.J. Braun, S.R. Coriell, R.F. Sekerka, Phys. Rev. E 48 (1993) 2016. [17] A.A. Wheeler, G.B. McFadden, W.J. Boettinger, Proc. R. Soc. London Ser. A 452 (1996) 495. [18] G. Caginalp, W. Xie, Phys. Rev. E 48 (1993) 1897. [19] J. Tiaden, B. Nestler, H.J. Diepers, I. Steinbach, Physica D 115 (1998) 73. [20] S.G. Kim, W.T. Kim, T. Suzuki, Phys. Rev. E 58 (1998) 3316. [21] W. Losert, D.A. Stillman, H.Z. Cummins, P. Kopczyn´ski, W.J. Rappel, A. Karma, Phys. Rev. E 58 (1998) 7492. [22] H. Lo¨wen, J. Bechhoefer, L.S. Tuckerman, Phys. Rev. A 45 (1992) 2399. [23] S.G. Kim, W.T. Kim, T. Suzuki, Phys. Rev. E 60 (1999) 7186. [24] P.-R. Cha, D.-H. Yeon, J.-K. Yoon, Acta Mater. 49 (2001) 3295. [25] D.-H. Yeon, P.-R. Cha, J.-K. Yoon, Scripta Mater. 45 (2001) 661. [26] S. Crusius, G. Inden, U. Knoop, L. Hoglund, J. A˚gren, Z. Metallkd. 83 (1992) 673. [27] B.-J. Lee, K.H. Oh, Z. Metallkd. 87 (1996) 195. [28] L.-Q. Chen, Ann. Rev. Mater. Res. 32 (2002) 113. [29] W.J. Boettinger, J.A. Warren, C. Beckermann, A. Karma, Ann. Rev. Mater. Res. 32 (2002) 163. [30] A. Ludwig, Physica D 124 (1998) 271. [31] M.J. Aziz, T. Kaplan, Acta Metall. 36 (1988) 2335. [32] D. Turnbull, J. Phys. Chem. 66 (1962) 609.

ARTICLE IN PRESS P.-R. Cha et al. / Journal of Crystal Growth 274 (2005) 281–293 [33] M.J. Aziz, W.J. Boettinger, Acta Metall. 42 (1994) 527. [34] W. Kurz, D. Fisher, Fundamentals of Solidification, Trans. Tech. Publications, Aedermannsdorf, Schweiz, 1992. [35] R. Willnecker, D.M. Herlach, B. Feuerbacher, Phys. Rev. Lett. 62 (1989) 2707. [36] J.-O. Andersson, J. A˚gren, J. Appl. Phys. 72 (1992) 1350.

293

[37] J.S. Kirkaldy, D.J. Young, Diffusion in the Condensed State, Institute of Metals, London, 1987. [38] A. Hultgren, Trans. Am. Soc. Met. 39 (1947) 915. [39] J.A. Warren, W.J. Boettinger, Acta Metall. Mater. 43 (1995) 689. [40] W.H. Press, S.A. Teokoisky, W.T. Vettering, B.P. Flannery, The Numerical Recipies in Fortran, Cambridge University Press, Cambridge, 1992, p. 372.