- Email: [email protected]
A phase-field model with antitrapping current for multicomponent alloys with arbitrary thermodynamic properties
Acta Materialia 55 (2007) 4391–4399 www.elsevier.com/locate/actamat
A phase-field model with antitrapping current for multicomponent alloys with arbitrary thermodynamic properties Seong Gyoon Kim
*
Department of Materials Science and Engineering, Kunsan National University, Kunsan 573-701, Republic of Korea Received 16 March 2007; received in revised form 3 April 2007; accepted 3 April 2007 Available online 25 May 2007
Abstract The development of accurate and efficient computational tools for phase transformation is a prerequisite in order to predict the microstructure evolution during phase transformation in engineering alloys. Even though the antitrapping phase-field model (PFM) has been developed and is increasingly being used for this purpose, it has been limited to dilute binary alloys. In this study, the antitrapping PFM was extended to multicomponent systems with arbitrary thermodynamic properties. We showed that at the condition of vanishing chemical potential jump, the antitrapping term in the multicomponent diffusion equation remains unchanged with the case of binary dilute alloys. We then derived the relationship between the phase-field mobility and the real interface mobility at the thin interface limit, and showed that the multicomponent effect in the relationship can be expressed in a concise matrix form. Ó 2007 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Phase-field modeling; Multicomponent solidification; Multicomponent diffusion; Phase transformation kinetics
1. Introduction Recently, the phase-field model (PFM) has been increasingly used as a powerful computational tool for the prediction of microstructure formation. This model adopts a diffuse interface approach, where the phase state changes gradually from one phase into another across an interfacial region with a finite width. Any point within the interfacial region is assumed to be a mixture of the two phases, whose fractions vary gradually from zero to one across the interface. All the thermodynamic and kinetic variables then are assumed to follow the relevant mixture rules. Though the PFM was originally devised for studying dendritic pattern formation during the solidification of pure undercooled melts [1–6], it was quickly extended to modeling alloy solidification [7–12], solid-state phase transformation [13–17], multiphase transformations [10,18–26] and multicompo-
*
Tel./fax: +82 63 469 4732. E-mail address: [email protected]
nent transformations in alloy systems with arbitrary phase diagrams [27–36]. PFMs for transformations between two phases consist of a phase-field equation and a diffusion equation. Three parameters exist in the phase-field equation, which have definite relationships with the interface’s energy, width and mobility. For quantitative computations, the relationships should be precisely determined so that the interface dynamics of the PFM corresponds to that of the real interface. A simple way to determine the relationships is to set the interface width in the PFM as the real interface width. In this case, however, the computational grid size needs to be smaller than the real interface width of about 1 nm. Mesoscale computations then become almost impossible because of the small grid size. This stringent restriction of the interface width was overcome by Karma and Rappel’s remarkable findings [37–39]. For the solidification of pure substances with the same heat capacities and thermal diffusivities in solid and liquid forms, they showed that the dynamics of an interface with a vanishingly small width (that is, classical sharp interface dynamics) can be correctly described by a PFM with a thin, but finite, interface width
1359-6454/$30.00 Ó 2007 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.actamat.2007.04.004
4392
S.G. Kim / Acta Materialia 55 (2007) 4391–4399
if a new relationship between the phase-field mobility and the real interface mobility is adopted. In their thin interface PFM, the interface width needs to be much smaller than the characteristic length scales of the diffusion field as well as the interface curvature. But this restriction appears to be significantly less severe than that in previous PFMs, and the thin interface PFM enables computation of the microstructure evolution on practical scales. The idea of the thin interface model for pure substance solidification has been applied to the PFM of alloy solidification in order to improve the computational efficiency and accuracy. The first improved model was obtained by adopting the condition of equal chemical potentials and finding a new relationship between the phase-field mobility and the real interface mobility with the thin interface condition [12]. However, this model suffered from anomalous interface diffusion and/or anomalous chemical potential jump at the interface which induces an exaggerated solute-trapping effect [40,41]. Such anomalous interface effects become significant with increasing interface width in computation. Karma again found the solution [41,42]: all the anomalous interface effects could be suppressed by introducing an antitrapping current term into the diffusion equation which acts against the solute-trapping current driven by the chemical potential gradient. Following Karma’s finding, two similar methods [43,44] were proposed independently. Both methods are based on the fact that all the anomalous interface effects originate from the finite interface width in the diffusion equation, not in the phase-field equation. The anomalous interface effects can then be suppressed by decoupling the interface width 2nd in the diffusion equation from the width 2np in the phase-field equation and taking the limit 2nd ! 0. In numerical computations, however, the minimum 2nd value which one can take appears to be slightly above the grid size Dx [43,44]. Therefore the solute-trapping phenomenon persists to a certain extent, as long as Dx in the computation is much thicker than the real interface width. This is in contrast with the antitrapping PFM, where the solute-trapping phenomenon can be controlled independently of the grid size. Antitrapping PFMs are being increasingly used for quantitative computation of the solidification of dilute binary alloys [25,45–49] because of their high accuracy and numerical simplicity. Predicting the microstructure evolution during phase transformation in alloys has long been a dream of materials designers. Realizing this dream requires the development of accurate and efficient tools for numerical computations. The antitrapping PFM can be a good choice of such a tool, but it is still limited to binary dilute alloys. Most engineering alloys are multicomponent systems and their thermodynamic properties are usually described by a complex database. Therefore, the extension of the antitrapping PFM to multicomponent systems with arbitrary thermodynamic properties is an important step towards the computational design of materials.
In this study, we first derive the phase-field and diffusion equations for multicomponent systems with arbitrary thermodynamic properties under the equal chemical potential condition. We then introduce the antitrapping current in the multicomponent diffusion equation and show that the antitrapping term under the condition of vanishing chemical potential jump remains unchanged in the case of binary dilute alloys. We also derive a relationship between the phase-field mobility and the real interface mobility at the thin interface limit. Lastly, we briefly discuss a computational issue in obtaining phase-field mobility. 2. Governing equations 2.1. Variables and free energy functional The multicomponent system in this study consists of (n + 1) components in solid and liquid phases. The mole fractions of the ith solute ði ¼ 1; 2; . . . ; nÞ in the solid and liquid phases is denoted by ciS and ciL, respectively. The (n + 1)th component is the solvent atom. The free energy density of a phase is a function of the solute concentrations: f p ¼ ðc1p ; c1p ; . . . ; cip ; . . . ; cnp Þ
ð1Þ
where p = S for the solid phase or p = L for the liquid phase. The free energy density f of a point in the interfacial region which is defined as a mixture of the solid and liquid phases is assumed to follow a mixture rule f ¼ hð/Þf S þ ½1 hð/Þf L
ð2Þ
where / is the phase field, which is defined as / = 0 in bulk liquid, / = 1 in bulk solid and 0 < / < 1 in the interfacial region between them. The interpolation function hð/Þ is a monotonous function satisfying hð0Þ ¼ 0 and hð1Þ ¼ 1. The total free energy of the system is given by a functional Z F ¼ ððe2 =2Þjr/j2 þ -gð/Þ þ hð/Þf S þ ½1 hð/Þf L Þ dV V
ð3Þ where e is the gradient energy coefficient and gð/Þ ¼ /2 ð1 /Þ2 is the double-well potential and - is the potential height. The concentration ci of the ith solute at a given point is ci ¼ ciS in bulk solid and ci ¼ ciL in bulk liquid, and in the interfacial region it follows the mixture rule ci ¼ hð/ÞciS þ ½1 hð/ÞciL
ð4Þ
The compositions ciS and ciL in a given point are not independent of each other, but restricted by the equal chemical potential condition [12,24,34] of S of L ~i ¼ l ociS ociL
ð5Þ
By taking the gradient, this can be written into an equivalent form n X j¼1
fijS rcjS ¼
n X j¼1
fijL rcjL ¼ r~ li
ð6Þ
S.G. Kim / Acta Materialia 55 (2007) 4391–4399
where we used the notation fijp ¼ o2 f p =ocip ocjp . Imposing the equal chemical potential condition upon the solid and liquid phases at a point of the system has two advantages over the traditional equal composition condition [7,8]. The first is the relaxation of the restriction on the interface width in computation [12,24]. The second is that the profile of the equilibrium phase-field gradient becomes symmetric, which suppresses the anomalous nonlinear part in the Gibbs–Thomson effect [40] in the thin interface limit. 2.2. Phase-field and diffusion equations The phase-field and diffusion equations for multicomponent systems have recently been derived in a rigorous way [34]. For a clearer and more consistent presentation, we briefly derive the equations again, but in a slightly different way. The concentration ci and phase field / are obtained at every time step of the computation as the solutions of the diffusion and phase-field equations. Then the variables ciS and ciL can be obtained from the nonlinear simultaneous Eqs. (4) and (5). Thus ciS and ciL are functions of ci and /. Differentiating Eq. (4) with respect to cj and / yields hð/Þ
ociS ociL þ ½1 hð/Þ ¼ dij ocj ocj
ð7Þ
ociS ociL þ ½1 hð/Þ ¼ h0 ð/ÞðciL ciS Þ o/ o/
ð8Þ
and hð/Þ
where dij is the Kronecker delta. Then the thermodynamic driving force of solidification becomes
S
S
ðciS ciL Þ~ li Þ
ð12Þ where we used Eq. (5). Thus Eq. (11) becomes n X oci ¼r M ik r~ lk ot k¼1
ð13Þ
We assume that the mobility Mij follows a mixture rule M ik ¼ hð/ÞM Sik þ ½1 hð/ÞM Lik
ð14Þ
where M Sik and M Lik are the diffusion mobilities in the solid and liquid phases, respectively. Inserting Eq. (14) into Eq. (13) and using Eq. (6), we rewrite Eq. (11) as n X oci lk ¼ r hð/ÞM Sik þ½1hð/ÞM Lik r~ ot k¼1 ! n n n X X X S S L L ¼ r hð/Þ M ik fkj rcjS þ½1hð/Þ M ik fkj rcjL k¼1
j¼1
j¼1
ð16Þ
where Dpij is the diffusivity in solid (p = S) or liquid (p = L) phase, we get the diffusion equation ! n n X X oci S L ¼ r hð/Þ Dij rcjS þ ½1 hð/Þ Dij rcjL ot j¼1 j¼1
¼ h ð/Þðf f Þ n X ociS ociL ~i þ þ ½1 hð/Þ hð/Þ l o/ o/ i¼1 n X
dF oF of S of L ¼ ¼ hð/Þ þ ½1 hð/Þ dck ock ock ock n S X of ocjS of L ocjL ¼ hð/Þ þ ½1 hð/Þ ocjS ock ocjL ock j¼1 n n X X ocjS ocjL ~k ~j ¼ l ~j ¼ hð/Þ þ ½1 hð/Þ djk l ¼ l ock ock j¼1 j¼1
k¼1
L
¼ h0 ð/Þðf S f L þ
where Mik is the diffusion mobility. The derivative dF =dck can be modified as following:
If we put n X Dpij ¼ M pik fkjp
L
¼ h ð/Þðf f Þ n X of S ociS of L ociL þ þ ½1 hð/Þ hð/Þ ociS o/ ociL o/ i¼1 0
The time evolution of the diffusion field under the mass conservation condition is assumed to follow n X oci dF ¼r M ik r ð11Þ dc ot k k¼1
ð15Þ
of of S of L ¼ h0 ð/Þðf S f L Þ þ hð/Þ þ ½1 hð/Þ o/ o/ o/ 0
4393
ð17Þ
ð9Þ
i¼1
where Eqs. (2), (5) and (8) are used in the first, third and fourth lines, respectively. The time evolution of the phase field is assumed to follow o/=ot ¼ M / dF =d/, which is explicitly written as n X 1 o/ dgð/Þ dhð/Þ S ðciS ciL Þ~ li Þ ¼ e2 r2 / ðf f L þ M / ot d/ d/ i¼1
ð10Þ
where M/ is the phase-field mobility.
3. Chemical potential jump As shown by Almgren [40], several anomalous interface effects appear when one takes a finite interface width in an alloy PFM. They become significant as one increases the width for enhanced computational efficiency. Even though each of the anomalous interface effects can be effectively suppressed by adopting the relevant interpolation function with a specific symmetry, it appears that not all of them can be suppressed simultaneously. For dilute binary alloys with
4394
S.G. Kim / Acta Materialia 55 (2007) 4391–4399
DS DL , Karma solved this problem by introducing an antitrapping term in the diffusion equation [41,42]. We here extend the antitrapping method to the cases of the arbitrary multicomponent alloys. The assumption DS DL is, however, maintained. Apart from its reasonable validity in most alloys, this allows the concentration (or chemical potential) profile at steady state during solidification to be determined unambiguously, which is indispensable for the present work. Under this assumption, we rewrite the phase-field equation 1 o/ dgð/Þ ¼ e2 r 2 / M / ot d/ n X dhp ð/Þ fS fL þ ðciS ciL Þ~ li d/ i¼1
! ð18Þ
and introduce an antitrapping term into the diffusion equation: n X oci o/ r/ ð19Þ ¼ r ½1 hd ð/Þ DLij rcjL þ r ai ot jr/j ot j¼1 where ai is a function of ciS and ciL, and ci is given by ci ¼ hr ð/ÞciS þ ½1 hr ð/ÞciL
ð20Þ
The original interpolation functions hð/Þ in Eqs. (4), (10) and (15) were discriminated by the subscripts (p, d and r) used in Eqs. (18)–(20). Even though a single function hð/Þ must be adopted in the rigorous thermodynamic derivation of the previous section, it is not the case in mapping the diffuse interface model onto the classical sharp interface model [41,42]; one rather can choose their detailed forms within 0 < / < 1. However, a specific symmetry in their functional forms must be imposed in order to suppress the anomalous interface effects [40–42], such as interface diffusion and interface stretching. This symmetry condition is equivalent to the requirement that the positions of the effective sharp interfaces for the driving force action (hp), diffusivity change (hd) and solute partitioning (hr) must be in accordance with that of the effective Gibbs–Thomson interface which is the symmetry axis position of the poten2 tial gð/Þ. For the potential gð/Þ ¼ /2 ð1 /Þ having a symmetry axis at / ¼ 1=2, the simple interpolation functions such as /, /2 ð3 2/Þ and /3 ð10 15/ þ 6/2 Þ satisfy the symmetry condition. Even after the choice of functions, there remains an anomalous interface effect: the chemical potential jump at the effective sharp interface. The chemical potential jump can be understood by looking into the composition or chemical potential profile around the interfacial region. Let us consider a one-dimensional solidifying system at an instantaneous steady state with an interface velocity V. We assume that the interface width is sufficiently smaller than the diffusion boundary layer width in liquid, that is, the thin interface condition. The composition profile across the interface ðn < x < nÞ is shown schematically in Fig. 1, where the composition profile ciS ðxÞ of the ith solute in the solid, ciL ðxÞ in the liquid and the mixture composition ci ðxÞ defined by Eq. (20) are
Fig. 1. Composition profiles across the interface ðn < x < nÞ. The composition profile ciS ðxÞ of ith solute in the solid, ciL ðxÞ in the liquid and the mixture composition ci ðxÞ are denoted by the lower thick curve, the upper thick curve and the dotted curve, respectively. The origin (x = 0) was defined as the position with / ¼ 0:5. The dashed lines are the extrapolations of the linear parts in ciS ðxÞ and ciL ðxÞ into the interfacial region.
denoted by the lower thick curve, the upper thick curve and the dotted curve, respectively. The origin (x = 0) is defined as the position with / ¼ 0:5. The chemical potential ~i ðxÞ, which is not shown in the figure, appears to be profile l of a similar pattern to ciS ðxÞ or ciL ðxÞ. One should note that ~i can be defined on all the points in the system ciS, ciL and l by virtue of the condition (5). If one of them is known at a given point, the other two at that point are fixed by the condition (5). We take ciL as the independent variable in the whole space of the system. There exist two straight parts in the profile of ciL ðxÞ. One is at the bulk solid side and the other at the bulk liquid side near the interface. We can extrapolate these two straight parts into the interfacial regions, as shown by the dashed lines in Fig. 1. Then the dynamics of the diffuse interface with the composition profiles of the thick curves represents effectively that of the classical sharp interface with the composition profiles of the dotted lines. These two extrapolated lines from the solid and liquid sides intersect the vertical axis at ciL ¼ c iL and ciL ¼ cþ . For an interface with a finite width, there exist a iL þ finite difference between c and c . This makes a correiL iL sponding difference in chemical potential, which has been called the chemical potential jump [41,42]. For multicomponent systems with arbitrary thermodynamic properties, as for dilute binary alloys [41,42], the chemical potential jump can be suppressed by a suitable choice of the interpolation functions and aðciS ; ciL Þ, that is, by balancing the anomalous solute trapping arising from the diffusion through the thick interface with the antitrapping current. The procedure for balancing the solute trapping with the antitrapping current is straightforward. First, we find the composition profile ciL ðxÞ by solving the steady-state diffusion equation. Next, we extract the straight part from the þ profile ciL ðxÞ and then get c iL and ciL . Lastly, by putting
S.G. Kim / Acta Materialia 55 (2007) 4391–4399 þ c iL ¼ ciL , we determine the interpolation functions and aðciS ; ciL Þ for the condition of vanishing chemical potential jump. We start by writing the steady-state form of the diffusion Eq. (19): ! n X dci d d d/ L dcjL V Dij ½1 hd ð/Þ ai ¼ þV dx dx dx dx dx j¼1
ð21Þ for the one-dimensional system shown in Fig. 1. Integrating this equation yields n X dcjL d/ ð22Þ þ V ai V ðci ci Þ ¼ ½1 hd ð/Þ DLij dx dx j¼1 where ci is an integration constant. It then follows that ci ¼ ci ðnÞ because hd ð/Þ ! 1 at / ! 1 and d/=dx ! 0 at x ! n. Following the matrix notations of Ref. [50], Eq. (22) can be written as V L dcL ½D ½BÞ ð23Þ ¼ 1 hd ð/Þ dx where the elements of the square matrix ½DL , the column vectors ½dcL =dxÞ and ½BÞ are given by DLij , dcjL =dx and d/ Bi ¼ ci ci þ ai dx
ð24Þ
respectively. We then get the composition gradient dcL V 1 ½DL ½BÞ ¼ 1 hd ð/Þ dx
ð25Þ
~ 1, where D ~ is the At the thin interface limit P nV =D average diffusivity at the interfacial region and P is the interface Pecklet number, the gradient can be approximated as dcL V 1 ½DL ½B0 Þ ð26Þ 1 hd ð/0 Þ dx 0
within the first order of P, where /0 and ½B Þ denote / and ½BÞ in the equilibrium state with V = 0. The equilibrium phase-field profile /0 is the solution of e2 d2 /0 =dx2 ¼ -gð/0 Þ or equivalently pffiffiffiffiffiffi pffiffiffiffiffiffi d/0 2- pffiffiffiffiffiffiffiffiffiffiffiffi 2gð/0 Þ ¼ /0 ð1 /0 Þ ¼ ð27Þ e e dx
4395
where d Lij and B0j are the element of ½DL 1 and ½B0 Þ, respectively. Integrating this equation yields the composition profile Z x ciL ðxÞ ¼ c V Qi ð/0 Þ dx ð30Þ iL n
c iL
is denoted in Fig. 1 and we define n X 1 d L B0 Qi ð/0 Þ 1 hd ð/0 Þ j¼1 ij j where
ð31Þ
Let us look at the function Qi ð/0 Þ. Its denominator vanishes at /0 ! 1. Nevertheless, we assume that Qi ð/0 Þ ! 0 at /0 ! 1. It can be true by suitable choices of hd ð/0 Þ, hr ð/0 Þ and aei , as shown later. Also we see that n X d Lij ðcejl cejS Þ ð32Þ Qi ð0Þ ¼ j¼1
because all of hd ð/0 Þ, hr ð/0 Þ and gð/0 Þ vanish at /0 ! 0. The profile of Qi ð/0 Þ across the interface is shown schematically in Fig. 2. At n < x D=V , the integration in Eq. (30) can be approximated as Z x Qi ð/0 Þ dx ¼ S 1 þ xQi ð0Þ S 2 ð33Þ n
where S1 and S2 are the areas of the shaded regions at x < 0 and x > 0 in Fig. 2, respectively. Then the linear part clin iL ðxÞ at n < x D=V in the profile ciL ðxÞ, which is the upper dashed line in Fig. 1, can be found as clin iL ðxÞ ¼ ciL V ½S 1 þ xQi ð0Þ S 2
ð34Þ
clin iL ðxÞ
When the straight line is extrapolated into the interfacial region, it intersects the x = 0 axis at the position denoted as cþ iL in Fig. 1. It then follows that lin cþ iL ciL ð0Þ ¼ ciL V ðS 1 S 2 Þ
ð35Þ
The condition of vanishing chemical potential jump, þ c iL ¼ ciL , is satisfied when S 1 ¼ S 2 , that is
Also, using Eqs. (20) and (27), ½B0 Þ is given by B0i ¼ cei ceiS þ aei ¼
ðceiL
ceiS Þ½1
d/0 dx hr ð/0 Þ
aei
pffiffiffiffiffiffi 2/0 ð1 /0 Þ e
ð28Þ
where ceiS and ceiL are the equilibrium compositions of the solid and liquid, respectively, and aei is the antitrapping coefficient in the equilibrium state. Eq. (26) may be rewritten into an explicit form n X dciL V ¼ d L B0 1 hd ð/0 Þ j¼1 ij j dx
ð29Þ
Fig. 2. Profile of the integrand Qi ð/0 Þ in Eq. (30) across the interface. S1 and S2 are the areas of the shaded regions at x < 0 and x > 0, respectively.
4396
Z
S.G. Kim / Acta Materialia 55 (2007) 4391–4399
0
Qi ð/0 Þ dx ¼
Z
1
1
½Qi ð0Þ Qi ð/0 Þ dx
ð36Þ
0
Considering the symmetry of /0 in Eq. (27) with respect to x, a simple choice to satisfy the condition (36) is Qi ð/0 Þ ¼ ð1 /0 ÞQi ð0Þ
ð37Þ
From B0i , Qi ð/0 Þ and Qi ð0Þ, given by Eqs. (28), (31) and (32), respectively, Eq. (37) can be written as ! pffiffiffiffiffiffi n X 1 /0 2e e 1 hr ð/0 Þ e /0 aj d ij ðcjL cjS Þ 1 hd ð/0 Þ e 1 hd ð/0 Þ j¼1 ¼ ð1 /0 Þ
n X
d ij ðcejL cejS Þ
ð38Þ
j¼1
There exist combinations of hr ð/0 Þ, hd ð/0 Þ and aei which satisfy the above condition. One set of the simplest p choices ffiffiffiffiffiffi is as follows; hr ð/0 Þ ¼ hd ð/0 Þ ¼ /0 and aei ¼ ð 2-=eÞ ðceiL ceiS Þ. (With this combination, one can see that Qi ð/0 Þ ! 0 at /0 ! 1, which was assumed before.) Therefore the interpolation functions and antitrapping function in the diffusion Eq. (19) are given by hr ð/Þ ¼ hd ð/Þ ¼ /
ð39Þ
and pffiffiffiffiffiffi 2ðciL ciS Þ ai ¼ ð40Þ e for vanishing chemical potential jump. Remarkably, Eq. (40) for ai which was derived for multicomponent systems with arbitrary thermodynamic properties is basically identical to that for a dilute binary system [41,42]. The apparent difference is due to the different definition of the phase field in Refs. [41,42], where the phase-field values of the bulk solid and liquid phases were defined as / ¼ 1 and / = 1, respectively. 4. Phase-field mobility We determined the interpolation functions hd ð/Þ and hr ð/Þ and the antitrapping function ai for suppressing the anomalous interface effects in the thin interface PFM of multicomponent systems with arbitrary thermodynamic properties. However, the parameters e, - and M/ in the phase-field Eq. (18) remain undetermined. The parameters e and - can be found from their relationships with the interface width 2n and the interface energy r in the equilibrium state: Z /b e d/0 e / ð1 /a Þ ¼ pffiffiffiffiffiffi ln b ð41Þ 2n ¼ pffiffiffiffiffiffi / ð1 / Þ 2- /a 2- /a ð1 /b Þ 0 0 and 2 pffiffiffiffi Z 1 d/0 e 2 r¼e dx ¼ pffiffiffi ð42Þ dx 3 2 1 respectively, where we use Eq. (27) and define 2n as the width over which /0 changes from /a to /b.
The phase-field mobility M/ has a relationship with the interface mobility m defined by the ratio between the driving force and the interface velocity. In this section, the relationship will be found at the thin interface limit. Although the procedure to find M/ appears to be somewhat lengthy, the basic strategy is simple: first we acquire the composiþ tion profile ciL ðxÞ under the condition c iL ¼ ciL for vanish~i ðxÞ follow ing chemical potential jump. Then ciS ðxÞ and l from the equal chemical potential condition (5) if necessary. We next insert the composition and chemical potential profiles into the driving force term of the phase-field equation and extract the new driving force for the effective sharp interface with the straight composition profiles in the interfacial region as the dashed lines in Fig. 1. We then get the relationship between the new driving force and the interface velocity, which yields a relationship between the physical interface mobility m and the phase-field mobility M/ at the thin interface limit. With the interpolation functions (39) and the antitrapping function (40) for vanishing chemical potential jump, the function B0i in Eq. (28) becomes B0i ¼ ðceiL ceiS Þð1 /0 Þ2 and so the composition gradient (29) is written as n X dciL ¼ V ð1 /0 Þ d Lij ðcejL cejS Þ dx j¼1
ð43Þ
ð44Þ
In the phase-field Eq. (18), the thermodynamic driving force Df is given by n X Df ¼ f S f L ðciS ciL Þ~ li ð45Þ i¼1 S
L
~i in this form all vary from position to f ; f ; ciS ; ciL and l position. In order to find the interface velocity by solving the phase-field equation, all the variables may be needed as functions of position. However, an assumption of the thin interface greatly simplifies the situation so that one needs only the gradient (44) according to the following process. ~ 1, the thermodyAt the thin interface limit P ¼ nV =D namic state of the interface should be close to the local equilibrium state. Within the first order of P, we expand the variables around their equilibrium values: ciS ðxÞ ¼ ceiS þ dciS ðxÞ ciL ðxÞ ¼ ceiL þ dciL ðxÞ ~i ðxÞ ¼ l ~ei þ d~ l li ðxÞ n n X of L e X L L;e dciL ¼ f L;e þ ~ei dciL f ¼f þ l oc iL i¼1 i¼1 n n S e X X of dciS ¼ f S;e þ ~ei dciS f S ¼ f S;e þ l ociS i¼1 i¼1
ð46Þ ð47Þ ð48Þ ð49Þ ð50Þ
where we use Eq. (5) and the superscript e hereafter denotes the equilibrium value. Then the thermodynamic driving force within the first order of P is given by
S.G. Kim / Acta Materialia 55 (2007) 4391–4399
Df ¼ f L;e f S;e þ
n X
~ei ðdciL dciS Þ l
f L;e f S;e
i¼1
n X
ð51Þ
i¼1
¼ f L;e f S;e ¼ f L;e f S;e
i¼1 n X
ðceiL ceiS Þð~ lei þ d~ li Þ
ð52Þ
ðceiL ceiS Þ~ li ðxÞ
ð53Þ
f
where f is defined as n n n X X X f¼ ðceiL ceiS Þ fijL;e d Ljk ðcekL cekS Þ
f
j¼1
S;e
! n X dhp ð/Þ e e ðciL ciS Þ~ li ðxÞ d/ i¼1 ð54Þ
at instantaneous steady state. Integrating this equation after multiplying d/=dx on both sides results in Z 1 2 n X V d/ dx ¼ f L;e þ f S;e ðceiL ceiS Þ M / 1 dx i¼1 Z 1 dhp ð/Þ ~i ðxÞ dx ð55Þ l dx 1 because hp ð/Þ ! hp ð1Þ ¼ 1 at x ! 1 and hp ð/Þ ! hp ð0Þ ¼ 0 at x ! 1. The last integration on the right-hand side of Eq. (55) can be performed by parts: Z 1 Z 1 x¼þ1 dhp ð/Þ d~ li ~i ðxÞhp ð/Þx¼1 ~i ðxÞ dx ¼ l hp ð/Þ dx l dx 1 1 dx Z 0 n X dcjL dx d/ ð56Þ ¼ ~ l0i hp ð/Þ fijL;e d/ dx 1 j¼1 where we used a modified form of Eq. (6), n n d~ li X dcjL X dcjL ¼ fijL fijL;e dx dx dx j¼1 j¼1
ð57Þ
This can be written into a compact form 1
ðDc ¼ ceiL L 1
ð62Þ T ceiS , ½DcÞ ¼ ðDc , L 1 L 1
with the matrix notations L;e L L L ½G ¼ fij , ½M ¼ M ij and ½D ¼ ½G ½M from Eq. (16). It should be noted that the left-hand side of Eq. (60) is the thermodynamic driving force Df sharp for the ~0i . effective sharp interface with a chemical potential l Therefore Eq. (60) shows that the interface velocity V is proportional to the driving force for the effective sharp interface; V ¼ mDf sharp , where the interface mobility m is given by a relationship pffiffiffiffi 1 1 e pffiffiffi a2 pffiffiffiffiffiffi f ¼ ð63Þ m M / 3 2e 2The second term on the right-hand side of Eq. (63) is the first-order correction on the phase-field mobility at the thin interface limit. With the relationship (63), the phase-field mobility M/ can be determined from the experimentally measured m value so that the interface velocity of the PFM mimics that of the sharp interface in the solidification of real materials. In particular, for the infinite interface mobility, we obtain the phase-field mobility M/ ¼ 2 ð64Þ 3e a2 f at which the local equilibrium condition is effectively maintained at the interface. It is quite remarkable that the multicomponent effect in the relationship between M/ and m can be put into such a concise form as Eq. (62). It is worthwhile considering some simple cases. For binary alloys, the parameter f becomes 2
f¼ ~0i ¼ l ~0i ð1Þ denotes the chemical potential correand l sponding to the constant composition c iL in Fig. 1. Inserting Eq. (44) for dcjL =dx and Eq. (27) for d/=dx d/0 =dx into Eq. (56), we get Z 1 n n X dhp ð/Þ e X ~i ðxÞ l fijL;e d Ljk ðcekL cekS Þ dx ¼ ~ l0i þ a2 V pffiffiffiffiffiffi dx 2- j¼1 1 k¼1 ð58Þ
where a2 is a constant given by Z 1 Z 1 ð1 /0 Þhp ð/0 Þ hp ð/0 Þ pffiffiffiffiffiffiffiffiffiffiffiffi a2 ¼ d/0 d/0 ¼ /0 gð/ Þ 0 0 0
ð61Þ
k¼1
1
V d/ d2 / dgð/Þ ¼ e2 2 M / dx dx d/ þ
pffiffiffiffi 1 e pffiffiffi a2 pffiffiffiffiffiffi f M / 3 2e 2-
f ¼ ðDc½GL ½DL ½DcÞ ¼ ðDc½ML ½DcÞ
With this driving force, the phase-field equation for the one-dimensional solidification system is given by
L;e
ð60Þ
i¼1
i¼1
ðceiL ceiS Þ~ l0i ¼ V
i¼1
ðceiL þ dciL ceiS dciS Þð~ lei þ d~ li Þ n X
n X
4397
ðce1L ce1S Þ f11L;e DL11
The mobility M/ from Eq. (65) is similar, but slightly different from that in the previous study [12], which has also been derived for the binary alloy with arbitrary thermodynamic properties. This originates from the different forms of hd ð/Þ in the present study and Ref. [12]. For dilute binary alloys, the thermodynamic driving force Df in the phase-field equation is given by Df ¼
ð59Þ
If we insert Eq. (59) into Eq. (55) and use Eq. (27) for the integral on the left-hand side of Eq. (55), it follows that
ð65Þ
ð1 k SL ÞRT e ðc1L c1L Þ vm
ð66Þ
and the parameter f becomes f¼
ce1L ð1 k SL Þ2 RT vm DL11
ð67Þ
4398
S.G. Kim / Acta Materialia 55 (2007) 4391–4399
because f11L;e ¼ ðRT =vm Þð1=c1L Þ and k SL c1S =c1L ¼ ce1S =ce1L from the equal chemical potential condition under the dilute alloy approximation. The mobility M/ from Eqs. (67) and (63) is identical to that in Refs. [41,42]. Although they apparently look different from each other, it is due to the different definitions of concentration, phase-field and interface mobility in Refs. [41,42]. Another interesting case is when the off-diagonal terms in the diffusion mobility matrix vanish: M Lij ¼ 0 for i 6¼ j. In this case, Eq. (62) yields 2 n X ðceiL ceiS Þ f¼ ð68Þ M Lii i¼1 which is the sum of the independent contributions from all components. 5. Discussion The present PFM of solidification consists of the phasefield equation (18) with the parameter relationships (41), (42) and (63) for e, - and M/, respectively, the diffusion equation (19) with the mixture composition (20), and the equal chemical potential condition (5). As the PFM for dilute binary alloys [41,42], the present PFM for multicomponent solidification into a solid phase is mapped onto the classical sharp interface model under the conditions of ~ DS DL , n D=V and n Ri , where Ri is the radius of the interface curvature. A problem in getting the phase-field mobility M/ during numerical computation must be mentioned: in order to get the f value in Eq. (62) and thereby the M/ value from Eq. (63) or (64), the local equilibrium concentrations ceiS and ceiL of all components at the interface are needed. For binary alloys under a constant pressure these values are uniquely determined at a given temperature, corresponding to two end-concentrations of a tie line in the phase diagram. The situation in multicomponent alloys is very different: an infinite number of tie lines exist at a given temperature. The local equilibrium concentrations at the interface, that is, the tie line which the interface takes, cannot be uniquely determined, but depends on the diffusion history of all the components and varies with time towards the global equilibrium concentrations of the system. Thus, finding ceiS and ceiL at each grid point at every time step during numerical computation is not possible. However, this problem can be circumvented by simply replacing ceiS and ceiL with ciS and ciL in Eq. (62). In this case, f during the numerical computation is calculated on each grid point by using the instantaneous concentrations ciS and ciL on that grid point, without calling the equilibrium concentrations. The term including f in Eq. (63) is the first-order correction to M/ for the small, but finite P (interface Pecklet number) effect. Replacing the equilibrium concentrations with the instantaneous ones changes the relationship between m and M/ as the second order, which has been neglected not only in Eq. (63), but also in all previous treatments of the present study.
All the equations in this study have been derived for an isothermal and isotropic system. If the release of the latent heat can be ignored as in usual directional solidification experiments, however, the nonisothermal system can also be treated by taking the thermodynamic data as it is given by a function of temperature. The anisotropy in the interface energy and/or interface mobility can also be introduced by following the standard method [6,8]. Neither the temperature gradient nor anisotropy changes the parameter relationships (41), (42) and (63), as in the model for dilute binary alloys [41,42]. The PFM developed in this study is restricted to solidification into a single solid phase. Solidification of useful multicomponent alloys often accompanies two or more solid phases. Recently Folch and Plapp [25] developed an antitrapping PFM for eutectic solidification where two solid phases form simultaneously from a liquid phase. However, their model is restricted to dilute binary alloys. Developing an antitrapping PFM for multicomponent and multiphase solidification with arbitrary thermodynamic properties remains an important task for the future. 6. Conclusion In order to predict microstructure evolution during phase transformation in an alloy with a given composition, the development of accurate and efficient computational tools for phase transformation is prerequisite. Even though the antitrapping PFM has been developed and widely used as such a tool, it has been limited to dilute binary alloys. Most engineering alloys are multicomponent systems and their thermodynamic properties are usually given by a complex database. In this study, the antitrapping PFM for dilute binary alloys was extended to multicomponent systems with arbitrary thermodynamic properties. We showed that at the condition of vanishing chemical potential jump, the antitrapping term in the multicomponent diffusion equation remains unchanged from the case of binary dilute alloys. We then derived the relationship between the phase-field mobility and the real interface mobility at the thin interface limit, and showed that the multicomponent effect in the relationship can be put into a concise matrix form. Even though the present study is an important step towards the computational design of engineering materials, further development of the model to include multiphase systems and to relax the one-sided diffusion condition ðDS DL Þ remains an important task for the future. References [1] Langer JS. Directions in condensed matter. Singapore: World Scientific; 1986. [2] Collins JB, Levine H. Phys Rev B 1985;31:1669. [3] Caginalp G. Phys Rev A 1989;39:5887. [4] Kobayashi R. Physica D 1993;63:410. [5] Penrose O, Pife PC. Physica D 1990;43:44. [6] McFadden GB, Wheeler AA, Braun RJ, Coriell SR, Sekerka RF. Phys Rev E 1993;48:2016.
S.G. Kim / Acta Materialia 55 (2007) 4391–4399 [7] Wheeler AA, Boettinger WJ, McFadden GB. Phys Rev A 1992;45:7424. [8] Warren JA, Boettinger WJ. Acta Metall Mater 1995;43:689. [9] Steinbach I, Pezzolla F, Nestler B, Seeßelberg M, Prieler R, Schmitz GJ, et al. Physica D 1996;94:135. [10] Tiaden I, Nestler B, Dipers H-J, Steinbach I. Physica D 1998;115:73. [11] Kim SG, Kim WT, Suzuki T. Phys Rev E 1998;58:3316. [12] Kim SG, Kim WT, Suzuki T. Phys Rev E 1999;60:7186. [13] Dy LI, Chen LQ. Acta Mater 1998;46:2573. [14] Rubin G, Khchaturyan AG. Acta Mater 1995;47:1995. [15] Wen YH, Wang Y, Bendersky LA, Chen LQ. Acta Mater 2000;48:4125. [16] Cha P-R, Yeon D-H, Chung S-H. Scripta Mater 2005;52:1241. [17] Zhang W, Jin YM, Khachaturyan AG. Acta Mater 2007;55:565. [18] Chen LQ, Yang W. Phys Rev B 1994;50:15752. [19] Karma A. Phys Rev E 1994;49:2245. [20] Elder KR, Drotlet F, Kosterlitz JM, Grant M. Phys Rev Lett 1994;72:677. [21] Wheeler AA, McFadden GB, Boettinger JB. Proc R Soc Lond A 1996;452:495. [22] Fan D, Chen LQ. Acta Mater 1997;45:611. [23] Steinbach I, Pezzolla F. Physica D 1999;134:385. [24] Kim SG, Kim WT, Suzuki T, Ode M. J Cryst Growth 2004;261: 135. [25] Folch R, Plapp M. Phys Rev E 2005;72:011602. [26] Kim SG, Kim DI, Kim WT, Park YB. Phys Rev E 2006;74:061605. [27] Cha P-R, Yeon D-H, Yoon J-K. Acta Mater 2001;49:3295. [28] Zhu JZ, Liu ZK, Vaithyanathan V, Chen LQ. Scripta Mater 2002;46:401.
4399
[29] Kobayashi H, Ode M, Kim SG, Kim WT, Suzuki T. Scripta Mater 2003;48:689. [30] Qin RS, Wallach ER. Acta Mater 2003;51:6199. [31] Zhu JZ, Wang T, Zhou SH, Liu ZK, Chen LQ. Acta Mater 2004;52:833. [32] Cha P-R, Yeon D-H, Yoon J-K. J Cryst Growth 2005;274:281. [33] Qin RS, Wallach ER, Thomson RC. J Cryst Growth 2005;279: 163. [34] Eiken J, Bo¨ttger B, Steinbach I. Phys Rev E 2006;73:066122. [35] Bo¨ttger B, Eiken J, Steinbach I. Acta Mater 2006;54:2697. [36] Zhang R, Jing T, Jie W, Liu B. Acta Mater 2006;54:2235. [37] Karma A, Rappel WJ. Phys Rev Lett 1996;77:4050. [38] Karma A, Rappel WJ. Phys Rev E 1996;53:R3017. [39] Karma A, Rappel WJ. Phys Rev E 1998;57:4323. [40] Almgren RF. SIAM J Appl Math 1999;59:2086. [41] Karma A. Phys Rev Lett 2001;87:115701. [42] Echebarria B, Folch R, Karma A, Plapp M. Phys Rev E 2004;70:061604. [43] Amberg G. Phys Rev Lett 2003;91:265505. [44] Kim SG, Kim WT, Suzuki T. J Cryst Growth 2004;263:620. [45] Greenwood M, Haataja M, Provatas N. Phys Rev Lett 2004;93:246101. [46] Ramirez JC, Beckermann C, Karma A, Diepers H-J. Phys Rev E 2004;69:051607. [47] Ramirez JC, Beckermann C. Acta Mater 2005;53:1721. [48] Badillo A, Beckermann C. Acta Mater 2006;54:2015. [49] Pons AJ, Karma A, Akamatsu S, Newey M, Pomerance A, Singer H, et al. Phys Rev E 2007;75:021602. [50] Chen H, Morral JE. Acta Mater 1999;47:1175.
A phase-field model with antitrapping current for multicomponent alloys with arbitrary thermodynamic properties Seong Gyoon Kim
*
Department of Materials Science and Engineering, Kunsan National University, Kunsan 573-701, Republic of Korea Received 16 March 2007; received in revised form 3 April 2007; accepted 3 April 2007 Available online 25 May 2007
Abstract The development of accurate and efficient computational tools for phase transformation is a prerequisite in order to predict the microstructure evolution during phase transformation in engineering alloys. Even though the antitrapping phase-field model (PFM) has been developed and is increasingly being used for this purpose, it has been limited to dilute binary alloys. In this study, the antitrapping PFM was extended to multicomponent systems with arbitrary thermodynamic properties. We showed that at the condition of vanishing chemical potential jump, the antitrapping term in the multicomponent diffusion equation remains unchanged with the case of binary dilute alloys. We then derived the relationship between the phase-field mobility and the real interface mobility at the thin interface limit, and showed that the multicomponent effect in the relationship can be expressed in a concise matrix form. Ó 2007 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Phase-field modeling; Multicomponent solidification; Multicomponent diffusion; Phase transformation kinetics
1. Introduction Recently, the phase-field model (PFM) has been increasingly used as a powerful computational tool for the prediction of microstructure formation. This model adopts a diffuse interface approach, where the phase state changes gradually from one phase into another across an interfacial region with a finite width. Any point within the interfacial region is assumed to be a mixture of the two phases, whose fractions vary gradually from zero to one across the interface. All the thermodynamic and kinetic variables then are assumed to follow the relevant mixture rules. Though the PFM was originally devised for studying dendritic pattern formation during the solidification of pure undercooled melts [1–6], it was quickly extended to modeling alloy solidification [7–12], solid-state phase transformation [13–17], multiphase transformations [10,18–26] and multicompo-
*
Tel./fax: +82 63 469 4732. E-mail address: [email protected]
nent transformations in alloy systems with arbitrary phase diagrams [27–36]. PFMs for transformations between two phases consist of a phase-field equation and a diffusion equation. Three parameters exist in the phase-field equation, which have definite relationships with the interface’s energy, width and mobility. For quantitative computations, the relationships should be precisely determined so that the interface dynamics of the PFM corresponds to that of the real interface. A simple way to determine the relationships is to set the interface width in the PFM as the real interface width. In this case, however, the computational grid size needs to be smaller than the real interface width of about 1 nm. Mesoscale computations then become almost impossible because of the small grid size. This stringent restriction of the interface width was overcome by Karma and Rappel’s remarkable findings [37–39]. For the solidification of pure substances with the same heat capacities and thermal diffusivities in solid and liquid forms, they showed that the dynamics of an interface with a vanishingly small width (that is, classical sharp interface dynamics) can be correctly described by a PFM with a thin, but finite, interface width
1359-6454/$30.00 Ó 2007 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.actamat.2007.04.004
4392
S.G. Kim / Acta Materialia 55 (2007) 4391–4399
if a new relationship between the phase-field mobility and the real interface mobility is adopted. In their thin interface PFM, the interface width needs to be much smaller than the characteristic length scales of the diffusion field as well as the interface curvature. But this restriction appears to be significantly less severe than that in previous PFMs, and the thin interface PFM enables computation of the microstructure evolution on practical scales. The idea of the thin interface model for pure substance solidification has been applied to the PFM of alloy solidification in order to improve the computational efficiency and accuracy. The first improved model was obtained by adopting the condition of equal chemical potentials and finding a new relationship between the phase-field mobility and the real interface mobility with the thin interface condition [12]. However, this model suffered from anomalous interface diffusion and/or anomalous chemical potential jump at the interface which induces an exaggerated solute-trapping effect [40,41]. Such anomalous interface effects become significant with increasing interface width in computation. Karma again found the solution [41,42]: all the anomalous interface effects could be suppressed by introducing an antitrapping current term into the diffusion equation which acts against the solute-trapping current driven by the chemical potential gradient. Following Karma’s finding, two similar methods [43,44] were proposed independently. Both methods are based on the fact that all the anomalous interface effects originate from the finite interface width in the diffusion equation, not in the phase-field equation. The anomalous interface effects can then be suppressed by decoupling the interface width 2nd in the diffusion equation from the width 2np in the phase-field equation and taking the limit 2nd ! 0. In numerical computations, however, the minimum 2nd value which one can take appears to be slightly above the grid size Dx [43,44]. Therefore the solute-trapping phenomenon persists to a certain extent, as long as Dx in the computation is much thicker than the real interface width. This is in contrast with the antitrapping PFM, where the solute-trapping phenomenon can be controlled independently of the grid size. Antitrapping PFMs are being increasingly used for quantitative computation of the solidification of dilute binary alloys [25,45–49] because of their high accuracy and numerical simplicity. Predicting the microstructure evolution during phase transformation in alloys has long been a dream of materials designers. Realizing this dream requires the development of accurate and efficient tools for numerical computations. The antitrapping PFM can be a good choice of such a tool, but it is still limited to binary dilute alloys. Most engineering alloys are multicomponent systems and their thermodynamic properties are usually described by a complex database. Therefore, the extension of the antitrapping PFM to multicomponent systems with arbitrary thermodynamic properties is an important step towards the computational design of materials.
In this study, we first derive the phase-field and diffusion equations for multicomponent systems with arbitrary thermodynamic properties under the equal chemical potential condition. We then introduce the antitrapping current in the multicomponent diffusion equation and show that the antitrapping term under the condition of vanishing chemical potential jump remains unchanged in the case of binary dilute alloys. We also derive a relationship between the phase-field mobility and the real interface mobility at the thin interface limit. Lastly, we briefly discuss a computational issue in obtaining phase-field mobility. 2. Governing equations 2.1. Variables and free energy functional The multicomponent system in this study consists of (n + 1) components in solid and liquid phases. The mole fractions of the ith solute ði ¼ 1; 2; . . . ; nÞ in the solid and liquid phases is denoted by ciS and ciL, respectively. The (n + 1)th component is the solvent atom. The free energy density of a phase is a function of the solute concentrations: f p ¼ ðc1p ; c1p ; . . . ; cip ; . . . ; cnp Þ
ð1Þ
where p = S for the solid phase or p = L for the liquid phase. The free energy density f of a point in the interfacial region which is defined as a mixture of the solid and liquid phases is assumed to follow a mixture rule f ¼ hð/Þf S þ ½1 hð/Þf L
ð2Þ
where / is the phase field, which is defined as / = 0 in bulk liquid, / = 1 in bulk solid and 0 < / < 1 in the interfacial region between them. The interpolation function hð/Þ is a monotonous function satisfying hð0Þ ¼ 0 and hð1Þ ¼ 1. The total free energy of the system is given by a functional Z F ¼ ððe2 =2Þjr/j2 þ -gð/Þ þ hð/Þf S þ ½1 hð/Þf L Þ dV V
ð3Þ where e is the gradient energy coefficient and gð/Þ ¼ /2 ð1 /Þ2 is the double-well potential and - is the potential height. The concentration ci of the ith solute at a given point is ci ¼ ciS in bulk solid and ci ¼ ciL in bulk liquid, and in the interfacial region it follows the mixture rule ci ¼ hð/ÞciS þ ½1 hð/ÞciL
ð4Þ
The compositions ciS and ciL in a given point are not independent of each other, but restricted by the equal chemical potential condition [12,24,34] of S of L ~i ¼ l ociS ociL
ð5Þ
By taking the gradient, this can be written into an equivalent form n X j¼1
fijS rcjS ¼
n X j¼1
fijL rcjL ¼ r~ li
ð6Þ
S.G. Kim / Acta Materialia 55 (2007) 4391–4399
where we used the notation fijp ¼ o2 f p =ocip ocjp . Imposing the equal chemical potential condition upon the solid and liquid phases at a point of the system has two advantages over the traditional equal composition condition [7,8]. The first is the relaxation of the restriction on the interface width in computation [12,24]. The second is that the profile of the equilibrium phase-field gradient becomes symmetric, which suppresses the anomalous nonlinear part in the Gibbs–Thomson effect [40] in the thin interface limit. 2.2. Phase-field and diffusion equations The phase-field and diffusion equations for multicomponent systems have recently been derived in a rigorous way [34]. For a clearer and more consistent presentation, we briefly derive the equations again, but in a slightly different way. The concentration ci and phase field / are obtained at every time step of the computation as the solutions of the diffusion and phase-field equations. Then the variables ciS and ciL can be obtained from the nonlinear simultaneous Eqs. (4) and (5). Thus ciS and ciL are functions of ci and /. Differentiating Eq. (4) with respect to cj and / yields hð/Þ
ociS ociL þ ½1 hð/Þ ¼ dij ocj ocj
ð7Þ
ociS ociL þ ½1 hð/Þ ¼ h0 ð/ÞðciL ciS Þ o/ o/
ð8Þ
and hð/Þ
where dij is the Kronecker delta. Then the thermodynamic driving force of solidification becomes
S
S
ðciS ciL Þ~ li Þ
ð12Þ where we used Eq. (5). Thus Eq. (11) becomes n X oci ¼r M ik r~ lk ot k¼1
ð13Þ
We assume that the mobility Mij follows a mixture rule M ik ¼ hð/ÞM Sik þ ½1 hð/ÞM Lik
ð14Þ
where M Sik and M Lik are the diffusion mobilities in the solid and liquid phases, respectively. Inserting Eq. (14) into Eq. (13) and using Eq. (6), we rewrite Eq. (11) as n X oci lk ¼ r hð/ÞM Sik þ½1hð/ÞM Lik r~ ot k¼1 ! n n n X X X S S L L ¼ r hð/Þ M ik fkj rcjS þ½1hð/Þ M ik fkj rcjL k¼1
j¼1
j¼1
ð16Þ
where Dpij is the diffusivity in solid (p = S) or liquid (p = L) phase, we get the diffusion equation ! n n X X oci S L ¼ r hð/Þ Dij rcjS þ ½1 hð/Þ Dij rcjL ot j¼1 j¼1
¼ h ð/Þðf f Þ n X ociS ociL ~i þ þ ½1 hð/Þ hð/Þ l o/ o/ i¼1 n X
dF oF of S of L ¼ ¼ hð/Þ þ ½1 hð/Þ dck ock ock ock n S X of ocjS of L ocjL ¼ hð/Þ þ ½1 hð/Þ ocjS ock ocjL ock j¼1 n n X X ocjS ocjL ~k ~j ¼ l ~j ¼ hð/Þ þ ½1 hð/Þ djk l ¼ l ock ock j¼1 j¼1
k¼1
L
¼ h0 ð/Þðf S f L þ
where Mik is the diffusion mobility. The derivative dF =dck can be modified as following:
If we put n X Dpij ¼ M pik fkjp
L
¼ h ð/Þðf f Þ n X of S ociS of L ociL þ þ ½1 hð/Þ hð/Þ ociS o/ ociL o/ i¼1 0
The time evolution of the diffusion field under the mass conservation condition is assumed to follow n X oci dF ¼r M ik r ð11Þ dc ot k k¼1
ð15Þ
of of S of L ¼ h0 ð/Þðf S f L Þ þ hð/Þ þ ½1 hð/Þ o/ o/ o/ 0
4393
ð17Þ
ð9Þ
i¼1
where Eqs. (2), (5) and (8) are used in the first, third and fourth lines, respectively. The time evolution of the phase field is assumed to follow o/=ot ¼ M / dF =d/, which is explicitly written as n X 1 o/ dgð/Þ dhð/Þ S ðciS ciL Þ~ li Þ ¼ e2 r2 / ðf f L þ M / ot d/ d/ i¼1
ð10Þ
where M/ is the phase-field mobility.
3. Chemical potential jump As shown by Almgren [40], several anomalous interface effects appear when one takes a finite interface width in an alloy PFM. They become significant as one increases the width for enhanced computational efficiency. Even though each of the anomalous interface effects can be effectively suppressed by adopting the relevant interpolation function with a specific symmetry, it appears that not all of them can be suppressed simultaneously. For dilute binary alloys with
4394
S.G. Kim / Acta Materialia 55 (2007) 4391–4399
DS DL , Karma solved this problem by introducing an antitrapping term in the diffusion equation [41,42]. We here extend the antitrapping method to the cases of the arbitrary multicomponent alloys. The assumption DS DL is, however, maintained. Apart from its reasonable validity in most alloys, this allows the concentration (or chemical potential) profile at steady state during solidification to be determined unambiguously, which is indispensable for the present work. Under this assumption, we rewrite the phase-field equation 1 o/ dgð/Þ ¼ e2 r 2 / M / ot d/ n X dhp ð/Þ fS fL þ ðciS ciL Þ~ li d/ i¼1
! ð18Þ
and introduce an antitrapping term into the diffusion equation: n X oci o/ r/ ð19Þ ¼ r ½1 hd ð/Þ DLij rcjL þ r ai ot jr/j ot j¼1 where ai is a function of ciS and ciL, and ci is given by ci ¼ hr ð/ÞciS þ ½1 hr ð/ÞciL
ð20Þ
The original interpolation functions hð/Þ in Eqs. (4), (10) and (15) were discriminated by the subscripts (p, d and r) used in Eqs. (18)–(20). Even though a single function hð/Þ must be adopted in the rigorous thermodynamic derivation of the previous section, it is not the case in mapping the diffuse interface model onto the classical sharp interface model [41,42]; one rather can choose their detailed forms within 0 < / < 1. However, a specific symmetry in their functional forms must be imposed in order to suppress the anomalous interface effects [40–42], such as interface diffusion and interface stretching. This symmetry condition is equivalent to the requirement that the positions of the effective sharp interfaces for the driving force action (hp), diffusivity change (hd) and solute partitioning (hr) must be in accordance with that of the effective Gibbs–Thomson interface which is the symmetry axis position of the poten2 tial gð/Þ. For the potential gð/Þ ¼ /2 ð1 /Þ having a symmetry axis at / ¼ 1=2, the simple interpolation functions such as /, /2 ð3 2/Þ and /3 ð10 15/ þ 6/2 Þ satisfy the symmetry condition. Even after the choice of functions, there remains an anomalous interface effect: the chemical potential jump at the effective sharp interface. The chemical potential jump can be understood by looking into the composition or chemical potential profile around the interfacial region. Let us consider a one-dimensional solidifying system at an instantaneous steady state with an interface velocity V. We assume that the interface width is sufficiently smaller than the diffusion boundary layer width in liquid, that is, the thin interface condition. The composition profile across the interface ðn < x < nÞ is shown schematically in Fig. 1, where the composition profile ciS ðxÞ of the ith solute in the solid, ciL ðxÞ in the liquid and the mixture composition ci ðxÞ defined by Eq. (20) are
Fig. 1. Composition profiles across the interface ðn < x < nÞ. The composition profile ciS ðxÞ of ith solute in the solid, ciL ðxÞ in the liquid and the mixture composition ci ðxÞ are denoted by the lower thick curve, the upper thick curve and the dotted curve, respectively. The origin (x = 0) was defined as the position with / ¼ 0:5. The dashed lines are the extrapolations of the linear parts in ciS ðxÞ and ciL ðxÞ into the interfacial region.
denoted by the lower thick curve, the upper thick curve and the dotted curve, respectively. The origin (x = 0) is defined as the position with / ¼ 0:5. The chemical potential ~i ðxÞ, which is not shown in the figure, appears to be profile l of a similar pattern to ciS ðxÞ or ciL ðxÞ. One should note that ~i can be defined on all the points in the system ciS, ciL and l by virtue of the condition (5). If one of them is known at a given point, the other two at that point are fixed by the condition (5). We take ciL as the independent variable in the whole space of the system. There exist two straight parts in the profile of ciL ðxÞ. One is at the bulk solid side and the other at the bulk liquid side near the interface. We can extrapolate these two straight parts into the interfacial regions, as shown by the dashed lines in Fig. 1. Then the dynamics of the diffuse interface with the composition profiles of the thick curves represents effectively that of the classical sharp interface with the composition profiles of the dotted lines. These two extrapolated lines from the solid and liquid sides intersect the vertical axis at ciL ¼ c iL and ciL ¼ cþ . For an interface with a finite width, there exist a iL þ finite difference between c and c . This makes a correiL iL sponding difference in chemical potential, which has been called the chemical potential jump [41,42]. For multicomponent systems with arbitrary thermodynamic properties, as for dilute binary alloys [41,42], the chemical potential jump can be suppressed by a suitable choice of the interpolation functions and aðciS ; ciL Þ, that is, by balancing the anomalous solute trapping arising from the diffusion through the thick interface with the antitrapping current. The procedure for balancing the solute trapping with the antitrapping current is straightforward. First, we find the composition profile ciL ðxÞ by solving the steady-state diffusion equation. Next, we extract the straight part from the þ profile ciL ðxÞ and then get c iL and ciL . Lastly, by putting
S.G. Kim / Acta Materialia 55 (2007) 4391–4399 þ c iL ¼ ciL , we determine the interpolation functions and aðciS ; ciL Þ for the condition of vanishing chemical potential jump. We start by writing the steady-state form of the diffusion Eq. (19): ! n X dci d d d/ L dcjL V Dij ½1 hd ð/Þ ai ¼ þV dx dx dx dx dx j¼1
ð21Þ for the one-dimensional system shown in Fig. 1. Integrating this equation yields n X dcjL d/ ð22Þ þ V ai V ðci ci Þ ¼ ½1 hd ð/Þ DLij dx dx j¼1 where ci is an integration constant. It then follows that ci ¼ ci ðnÞ because hd ð/Þ ! 1 at / ! 1 and d/=dx ! 0 at x ! n. Following the matrix notations of Ref. [50], Eq. (22) can be written as V L dcL ½D ½BÞ ð23Þ ¼ 1 hd ð/Þ dx where the elements of the square matrix ½DL , the column vectors ½dcL =dxÞ and ½BÞ are given by DLij , dcjL =dx and d/ Bi ¼ ci ci þ ai dx
ð24Þ
respectively. We then get the composition gradient dcL V 1 ½DL ½BÞ ¼ 1 hd ð/Þ dx
ð25Þ
~ 1, where D ~ is the At the thin interface limit P nV =D average diffusivity at the interfacial region and P is the interface Pecklet number, the gradient can be approximated as dcL V 1 ½DL ½B0 Þ ð26Þ 1 hd ð/0 Þ dx 0
within the first order of P, where /0 and ½B Þ denote / and ½BÞ in the equilibrium state with V = 0. The equilibrium phase-field profile /0 is the solution of e2 d2 /0 =dx2 ¼ -gð/0 Þ or equivalently pffiffiffiffiffiffi pffiffiffiffiffiffi d/0 2- pffiffiffiffiffiffiffiffiffiffiffiffi 2gð/0 Þ ¼ /0 ð1 /0 Þ ¼ ð27Þ e e dx
4395
where d Lij and B0j are the element of ½DL 1 and ½B0 Þ, respectively. Integrating this equation yields the composition profile Z x ciL ðxÞ ¼ c V Qi ð/0 Þ dx ð30Þ iL n
c iL
is denoted in Fig. 1 and we define n X 1 d L B0 Qi ð/0 Þ 1 hd ð/0 Þ j¼1 ij j where
ð31Þ
Let us look at the function Qi ð/0 Þ. Its denominator vanishes at /0 ! 1. Nevertheless, we assume that Qi ð/0 Þ ! 0 at /0 ! 1. It can be true by suitable choices of hd ð/0 Þ, hr ð/0 Þ and aei , as shown later. Also we see that n X d Lij ðcejl cejS Þ ð32Þ Qi ð0Þ ¼ j¼1
because all of hd ð/0 Þ, hr ð/0 Þ and gð/0 Þ vanish at /0 ! 0. The profile of Qi ð/0 Þ across the interface is shown schematically in Fig. 2. At n < x D=V , the integration in Eq. (30) can be approximated as Z x Qi ð/0 Þ dx ¼ S 1 þ xQi ð0Þ S 2 ð33Þ n
where S1 and S2 are the areas of the shaded regions at x < 0 and x > 0 in Fig. 2, respectively. Then the linear part clin iL ðxÞ at n < x D=V in the profile ciL ðxÞ, which is the upper dashed line in Fig. 1, can be found as clin iL ðxÞ ¼ ciL V ½S 1 þ xQi ð0Þ S 2
ð34Þ
clin iL ðxÞ
When the straight line is extrapolated into the interfacial region, it intersects the x = 0 axis at the position denoted as cþ iL in Fig. 1. It then follows that lin cþ iL ciL ð0Þ ¼ ciL V ðS 1 S 2 Þ
ð35Þ
The condition of vanishing chemical potential jump, þ c iL ¼ ciL , is satisfied when S 1 ¼ S 2 , that is
Also, using Eqs. (20) and (27), ½B0 Þ is given by B0i ¼ cei ceiS þ aei ¼
ðceiL
ceiS Þ½1
d/0 dx hr ð/0 Þ
aei
pffiffiffiffiffiffi 2/0 ð1 /0 Þ e
ð28Þ
where ceiS and ceiL are the equilibrium compositions of the solid and liquid, respectively, and aei is the antitrapping coefficient in the equilibrium state. Eq. (26) may be rewritten into an explicit form n X dciL V ¼ d L B0 1 hd ð/0 Þ j¼1 ij j dx
ð29Þ
Fig. 2. Profile of the integrand Qi ð/0 Þ in Eq. (30) across the interface. S1 and S2 are the areas of the shaded regions at x < 0 and x > 0, respectively.
4396
Z
S.G. Kim / Acta Materialia 55 (2007) 4391–4399
0
Qi ð/0 Þ dx ¼
Z
1
1
½Qi ð0Þ Qi ð/0 Þ dx
ð36Þ
0
Considering the symmetry of /0 in Eq. (27) with respect to x, a simple choice to satisfy the condition (36) is Qi ð/0 Þ ¼ ð1 /0 ÞQi ð0Þ
ð37Þ
From B0i , Qi ð/0 Þ and Qi ð0Þ, given by Eqs. (28), (31) and (32), respectively, Eq. (37) can be written as ! pffiffiffiffiffiffi n X 1 /0 2e e 1 hr ð/0 Þ e /0 aj d ij ðcjL cjS Þ 1 hd ð/0 Þ e 1 hd ð/0 Þ j¼1 ¼ ð1 /0 Þ
n X
d ij ðcejL cejS Þ
ð38Þ
j¼1
There exist combinations of hr ð/0 Þ, hd ð/0 Þ and aei which satisfy the above condition. One set of the simplest p choices ffiffiffiffiffiffi is as follows; hr ð/0 Þ ¼ hd ð/0 Þ ¼ /0 and aei ¼ ð 2-=eÞ ðceiL ceiS Þ. (With this combination, one can see that Qi ð/0 Þ ! 0 at /0 ! 1, which was assumed before.) Therefore the interpolation functions and antitrapping function in the diffusion Eq. (19) are given by hr ð/Þ ¼ hd ð/Þ ¼ /
ð39Þ
and pffiffiffiffiffiffi 2ðciL ciS Þ ai ¼ ð40Þ e for vanishing chemical potential jump. Remarkably, Eq. (40) for ai which was derived for multicomponent systems with arbitrary thermodynamic properties is basically identical to that for a dilute binary system [41,42]. The apparent difference is due to the different definition of the phase field in Refs. [41,42], where the phase-field values of the bulk solid and liquid phases were defined as / ¼ 1 and / = 1, respectively. 4. Phase-field mobility We determined the interpolation functions hd ð/Þ and hr ð/Þ and the antitrapping function ai for suppressing the anomalous interface effects in the thin interface PFM of multicomponent systems with arbitrary thermodynamic properties. However, the parameters e, - and M/ in the phase-field Eq. (18) remain undetermined. The parameters e and - can be found from their relationships with the interface width 2n and the interface energy r in the equilibrium state: Z /b e d/0 e / ð1 /a Þ ¼ pffiffiffiffiffiffi ln b ð41Þ 2n ¼ pffiffiffiffiffiffi / ð1 / Þ 2- /a 2- /a ð1 /b Þ 0 0 and 2 pffiffiffiffi Z 1 d/0 e 2 r¼e dx ¼ pffiffiffi ð42Þ dx 3 2 1 respectively, where we use Eq. (27) and define 2n as the width over which /0 changes from /a to /b.
The phase-field mobility M/ has a relationship with the interface mobility m defined by the ratio between the driving force and the interface velocity. In this section, the relationship will be found at the thin interface limit. Although the procedure to find M/ appears to be somewhat lengthy, the basic strategy is simple: first we acquire the composiþ tion profile ciL ðxÞ under the condition c iL ¼ ciL for vanish~i ðxÞ follow ing chemical potential jump. Then ciS ðxÞ and l from the equal chemical potential condition (5) if necessary. We next insert the composition and chemical potential profiles into the driving force term of the phase-field equation and extract the new driving force for the effective sharp interface with the straight composition profiles in the interfacial region as the dashed lines in Fig. 1. We then get the relationship between the new driving force and the interface velocity, which yields a relationship between the physical interface mobility m and the phase-field mobility M/ at the thin interface limit. With the interpolation functions (39) and the antitrapping function (40) for vanishing chemical potential jump, the function B0i in Eq. (28) becomes B0i ¼ ðceiL ceiS Þð1 /0 Þ2 and so the composition gradient (29) is written as n X dciL ¼ V ð1 /0 Þ d Lij ðcejL cejS Þ dx j¼1
ð43Þ
ð44Þ
In the phase-field Eq. (18), the thermodynamic driving force Df is given by n X Df ¼ f S f L ðciS ciL Þ~ li ð45Þ i¼1 S
L
~i in this form all vary from position to f ; f ; ciS ; ciL and l position. In order to find the interface velocity by solving the phase-field equation, all the variables may be needed as functions of position. However, an assumption of the thin interface greatly simplifies the situation so that one needs only the gradient (44) according to the following process. ~ 1, the thermodyAt the thin interface limit P ¼ nV =D namic state of the interface should be close to the local equilibrium state. Within the first order of P, we expand the variables around their equilibrium values: ciS ðxÞ ¼ ceiS þ dciS ðxÞ ciL ðxÞ ¼ ceiL þ dciL ðxÞ ~i ðxÞ ¼ l ~ei þ d~ l li ðxÞ n n X of L e X L L;e dciL ¼ f L;e þ ~ei dciL f ¼f þ l oc iL i¼1 i¼1 n n S e X X of dciS ¼ f S;e þ ~ei dciS f S ¼ f S;e þ l ociS i¼1 i¼1
ð46Þ ð47Þ ð48Þ ð49Þ ð50Þ
where we use Eq. (5) and the superscript e hereafter denotes the equilibrium value. Then the thermodynamic driving force within the first order of P is given by
S.G. Kim / Acta Materialia 55 (2007) 4391–4399
Df ¼ f L;e f S;e þ
n X
~ei ðdciL dciS Þ l
f L;e f S;e
i¼1
n X
ð51Þ
i¼1
¼ f L;e f S;e ¼ f L;e f S;e
i¼1 n X
ðceiL ceiS Þð~ lei þ d~ li Þ
ð52Þ
ðceiL ceiS Þ~ li ðxÞ
ð53Þ
f
where f is defined as n n n X X X f¼ ðceiL ceiS Þ fijL;e d Ljk ðcekL cekS Þ
f
j¼1
S;e
! n X dhp ð/Þ e e ðciL ciS Þ~ li ðxÞ d/ i¼1 ð54Þ
at instantaneous steady state. Integrating this equation after multiplying d/=dx on both sides results in Z 1 2 n X V d/ dx ¼ f L;e þ f S;e ðceiL ceiS Þ M / 1 dx i¼1 Z 1 dhp ð/Þ ~i ðxÞ dx ð55Þ l dx 1 because hp ð/Þ ! hp ð1Þ ¼ 1 at x ! 1 and hp ð/Þ ! hp ð0Þ ¼ 0 at x ! 1. The last integration on the right-hand side of Eq. (55) can be performed by parts: Z 1 Z 1 x¼þ1 dhp ð/Þ d~ li ~i ðxÞhp ð/Þx¼1 ~i ðxÞ dx ¼ l hp ð/Þ dx l dx 1 1 dx Z 0 n X dcjL dx d/ ð56Þ ¼ ~ l0i hp ð/Þ fijL;e d/ dx 1 j¼1 where we used a modified form of Eq. (6), n n d~ li X dcjL X dcjL ¼ fijL fijL;e dx dx dx j¼1 j¼1
ð57Þ
This can be written into a compact form 1
ðDc ¼ ceiL L 1
ð62Þ T ceiS , ½DcÞ ¼ ðDc , L 1 L 1
with the matrix notations L;e L L L ½G ¼ fij , ½M ¼ M ij and ½D ¼ ½G ½M from Eq. (16). It should be noted that the left-hand side of Eq. (60) is the thermodynamic driving force Df sharp for the ~0i . effective sharp interface with a chemical potential l Therefore Eq. (60) shows that the interface velocity V is proportional to the driving force for the effective sharp interface; V ¼ mDf sharp , where the interface mobility m is given by a relationship pffiffiffiffi 1 1 e pffiffiffi a2 pffiffiffiffiffiffi f ¼ ð63Þ m M / 3 2e 2The second term on the right-hand side of Eq. (63) is the first-order correction on the phase-field mobility at the thin interface limit. With the relationship (63), the phase-field mobility M/ can be determined from the experimentally measured m value so that the interface velocity of the PFM mimics that of the sharp interface in the solidification of real materials. In particular, for the infinite interface mobility, we obtain the phase-field mobility M/ ¼ 2 ð64Þ 3e a2 f at which the local equilibrium condition is effectively maintained at the interface. It is quite remarkable that the multicomponent effect in the relationship between M/ and m can be put into such a concise form as Eq. (62). It is worthwhile considering some simple cases. For binary alloys, the parameter f becomes 2
f¼ ~0i ¼ l ~0i ð1Þ denotes the chemical potential correand l sponding to the constant composition c iL in Fig. 1. Inserting Eq. (44) for dcjL =dx and Eq. (27) for d/=dx d/0 =dx into Eq. (56), we get Z 1 n n X dhp ð/Þ e X ~i ðxÞ l fijL;e d Ljk ðcekL cekS Þ dx ¼ ~ l0i þ a2 V pffiffiffiffiffiffi dx 2- j¼1 1 k¼1 ð58Þ
where a2 is a constant given by Z 1 Z 1 ð1 /0 Þhp ð/0 Þ hp ð/0 Þ pffiffiffiffiffiffiffiffiffiffiffiffi a2 ¼ d/0 d/0 ¼ /0 gð/ Þ 0 0 0
ð61Þ
k¼1
1
V d/ d2 / dgð/Þ ¼ e2 2 M / dx dx d/ þ
pffiffiffiffi 1 e pffiffiffi a2 pffiffiffiffiffiffi f M / 3 2e 2-
f ¼ ðDc½GL ½DL ½DcÞ ¼ ðDc½ML ½DcÞ
With this driving force, the phase-field equation for the one-dimensional solidification system is given by
L;e
ð60Þ
i¼1
i¼1
ðceiL ceiS Þ~ l0i ¼ V
i¼1
ðceiL þ dciL ceiS dciS Þð~ lei þ d~ li Þ n X
n X
4397
ðce1L ce1S Þ f11L;e DL11
The mobility M/ from Eq. (65) is similar, but slightly different from that in the previous study [12], which has also been derived for the binary alloy with arbitrary thermodynamic properties. This originates from the different forms of hd ð/Þ in the present study and Ref. [12]. For dilute binary alloys, the thermodynamic driving force Df in the phase-field equation is given by Df ¼
ð59Þ
If we insert Eq. (59) into Eq. (55) and use Eq. (27) for the integral on the left-hand side of Eq. (55), it follows that
ð65Þ
ð1 k SL ÞRT e ðc1L c1L Þ vm
ð66Þ
and the parameter f becomes f¼
ce1L ð1 k SL Þ2 RT vm DL11
ð67Þ
4398
S.G. Kim / Acta Materialia 55 (2007) 4391–4399
because f11L;e ¼ ðRT =vm Þð1=c1L Þ and k SL c1S =c1L ¼ ce1S =ce1L from the equal chemical potential condition under the dilute alloy approximation. The mobility M/ from Eqs. (67) and (63) is identical to that in Refs. [41,42]. Although they apparently look different from each other, it is due to the different definitions of concentration, phase-field and interface mobility in Refs. [41,42]. Another interesting case is when the off-diagonal terms in the diffusion mobility matrix vanish: M Lij ¼ 0 for i 6¼ j. In this case, Eq. (62) yields 2 n X ðceiL ceiS Þ f¼ ð68Þ M Lii i¼1 which is the sum of the independent contributions from all components. 5. Discussion The present PFM of solidification consists of the phasefield equation (18) with the parameter relationships (41), (42) and (63) for e, - and M/, respectively, the diffusion equation (19) with the mixture composition (20), and the equal chemical potential condition (5). As the PFM for dilute binary alloys [41,42], the present PFM for multicomponent solidification into a solid phase is mapped onto the classical sharp interface model under the conditions of ~ DS DL , n D=V and n Ri , where Ri is the radius of the interface curvature. A problem in getting the phase-field mobility M/ during numerical computation must be mentioned: in order to get the f value in Eq. (62) and thereby the M/ value from Eq. (63) or (64), the local equilibrium concentrations ceiS and ceiL of all components at the interface are needed. For binary alloys under a constant pressure these values are uniquely determined at a given temperature, corresponding to two end-concentrations of a tie line in the phase diagram. The situation in multicomponent alloys is very different: an infinite number of tie lines exist at a given temperature. The local equilibrium concentrations at the interface, that is, the tie line which the interface takes, cannot be uniquely determined, but depends on the diffusion history of all the components and varies with time towards the global equilibrium concentrations of the system. Thus, finding ceiS and ceiL at each grid point at every time step during numerical computation is not possible. However, this problem can be circumvented by simply replacing ceiS and ceiL with ciS and ciL in Eq. (62). In this case, f during the numerical computation is calculated on each grid point by using the instantaneous concentrations ciS and ciL on that grid point, without calling the equilibrium concentrations. The term including f in Eq. (63) is the first-order correction to M/ for the small, but finite P (interface Pecklet number) effect. Replacing the equilibrium concentrations with the instantaneous ones changes the relationship between m and M/ as the second order, which has been neglected not only in Eq. (63), but also in all previous treatments of the present study.
All the equations in this study have been derived for an isothermal and isotropic system. If the release of the latent heat can be ignored as in usual directional solidification experiments, however, the nonisothermal system can also be treated by taking the thermodynamic data as it is given by a function of temperature. The anisotropy in the interface energy and/or interface mobility can also be introduced by following the standard method [6,8]. Neither the temperature gradient nor anisotropy changes the parameter relationships (41), (42) and (63), as in the model for dilute binary alloys [41,42]. The PFM developed in this study is restricted to solidification into a single solid phase. Solidification of useful multicomponent alloys often accompanies two or more solid phases. Recently Folch and Plapp [25] developed an antitrapping PFM for eutectic solidification where two solid phases form simultaneously from a liquid phase. However, their model is restricted to dilute binary alloys. Developing an antitrapping PFM for multicomponent and multiphase solidification with arbitrary thermodynamic properties remains an important task for the future. 6. Conclusion In order to predict microstructure evolution during phase transformation in an alloy with a given composition, the development of accurate and efficient computational tools for phase transformation is prerequisite. Even though the antitrapping PFM has been developed and widely used as such a tool, it has been limited to dilute binary alloys. Most engineering alloys are multicomponent systems and their thermodynamic properties are usually given by a complex database. In this study, the antitrapping PFM for dilute binary alloys was extended to multicomponent systems with arbitrary thermodynamic properties. We showed that at the condition of vanishing chemical potential jump, the antitrapping term in the multicomponent diffusion equation remains unchanged from the case of binary dilute alloys. We then derived the relationship between the phase-field mobility and the real interface mobility at the thin interface limit, and showed that the multicomponent effect in the relationship can be put into a concise matrix form. Even though the present study is an important step towards the computational design of engineering materials, further development of the model to include multiphase systems and to relax the one-sided diffusion condition ðDS DL Þ remains an important task for the future. References [1] Langer JS. Directions in condensed matter. Singapore: World Scientific; 1986. [2] Collins JB, Levine H. Phys Rev B 1985;31:1669. [3] Caginalp G. Phys Rev A 1989;39:5887. [4] Kobayashi R. Physica D 1993;63:410. [5] Penrose O, Pife PC. Physica D 1990;43:44. [6] McFadden GB, Wheeler AA, Braun RJ, Coriell SR, Sekerka RF. Phys Rev E 1993;48:2016.
S.G. Kim / Acta Materialia 55 (2007) 4391–4399 [7] Wheeler AA, Boettinger WJ, McFadden GB. Phys Rev A 1992;45:7424. [8] Warren JA, Boettinger WJ. Acta Metall Mater 1995;43:689. [9] Steinbach I, Pezzolla F, Nestler B, Seeßelberg M, Prieler R, Schmitz GJ, et al. Physica D 1996;94:135. [10] Tiaden I, Nestler B, Dipers H-J, Steinbach I. Physica D 1998;115:73. [11] Kim SG, Kim WT, Suzuki T. Phys Rev E 1998;58:3316. [12] Kim SG, Kim WT, Suzuki T. Phys Rev E 1999;60:7186. [13] Dy LI, Chen LQ. Acta Mater 1998;46:2573. [14] Rubin G, Khchaturyan AG. Acta Mater 1995;47:1995. [15] Wen YH, Wang Y, Bendersky LA, Chen LQ. Acta Mater 2000;48:4125. [16] Cha P-R, Yeon D-H, Chung S-H. Scripta Mater 2005;52:1241. [17] Zhang W, Jin YM, Khachaturyan AG. Acta Mater 2007;55:565. [18] Chen LQ, Yang W. Phys Rev B 1994;50:15752. [19] Karma A. Phys Rev E 1994;49:2245. [20] Elder KR, Drotlet F, Kosterlitz JM, Grant M. Phys Rev Lett 1994;72:677. [21] Wheeler AA, McFadden GB, Boettinger JB. Proc R Soc Lond A 1996;452:495. [22] Fan D, Chen LQ. Acta Mater 1997;45:611. [23] Steinbach I, Pezzolla F. Physica D 1999;134:385. [24] Kim SG, Kim WT, Suzuki T, Ode M. J Cryst Growth 2004;261: 135. [25] Folch R, Plapp M. Phys Rev E 2005;72:011602. [26] Kim SG, Kim DI, Kim WT, Park YB. Phys Rev E 2006;74:061605. [27] Cha P-R, Yeon D-H, Yoon J-K. Acta Mater 2001;49:3295. [28] Zhu JZ, Liu ZK, Vaithyanathan V, Chen LQ. Scripta Mater 2002;46:401.
4399
[29] Kobayashi H, Ode M, Kim SG, Kim WT, Suzuki T. Scripta Mater 2003;48:689. [30] Qin RS, Wallach ER. Acta Mater 2003;51:6199. [31] Zhu JZ, Wang T, Zhou SH, Liu ZK, Chen LQ. Acta Mater 2004;52:833. [32] Cha P-R, Yeon D-H, Yoon J-K. J Cryst Growth 2005;274:281. [33] Qin RS, Wallach ER, Thomson RC. J Cryst Growth 2005;279: 163. [34] Eiken J, Bo¨ttger B, Steinbach I. Phys Rev E 2006;73:066122. [35] Bo¨ttger B, Eiken J, Steinbach I. Acta Mater 2006;54:2697. [36] Zhang R, Jing T, Jie W, Liu B. Acta Mater 2006;54:2235. [37] Karma A, Rappel WJ. Phys Rev Lett 1996;77:4050. [38] Karma A, Rappel WJ. Phys Rev E 1996;53:R3017. [39] Karma A, Rappel WJ. Phys Rev E 1998;57:4323. [40] Almgren RF. SIAM J Appl Math 1999;59:2086. [41] Karma A. Phys Rev Lett 2001;87:115701. [42] Echebarria B, Folch R, Karma A, Plapp M. Phys Rev E 2004;70:061604. [43] Amberg G. Phys Rev Lett 2003;91:265505. [44] Kim SG, Kim WT, Suzuki T. J Cryst Growth 2004;263:620. [45] Greenwood M, Haataja M, Provatas N. Phys Rev Lett 2004;93:246101. [46] Ramirez JC, Beckermann C, Karma A, Diepers H-J. Phys Rev E 2004;69:051607. [47] Ramirez JC, Beckermann C. Acta Mater 2005;53:1721. [48] Badillo A, Beckermann C. Acta Mater 2006;54:2015. [49] Pons AJ, Karma A, Akamatsu S, Newey M, Pomerance A, Singer H, et al. Phys Rev E 2007;75:021602. [50] Chen H, Morral JE. Acta Mater 1999;47:1175.