Table Of Content

Towards a quantitative phase-field model of two-phase solidification ∗ R. Folch and M. Plapp Laboratoire de Physique de la Matière Condensée, CNRS/Ećole Polytechnique, 91128 Palaiseau, France We construct a diffuse-interface model of two-phase solidification that quantitatively reproduces 4 theclassic free boundaryproblem on solid-liquid interfaces inthethin-interfacelimit. Convergence 0 tests and comparisons with boundary integral simulations of eutectic growth show good accuracy 0 for steady-state lamellae, but theresults for limit cycles depend on theinterface thickness through 2 thetrijunction behavior. This raises thefundamental issue of diffusemultiple-junction dynamics. n a PACSnumbers: 64.70.Dv,81.30.Fb,05.70.Ln J 4 1 Complex microstructures that ariseduring alloysolid- tweencapillarityanddiffusivebulktransportbetweenad- ification are a classical example of pattern formation [1] jacentsolidphasescangiverisetomorecomplexpatterns ] i and influence the mechanical properties of the finished and nonlinear phenomena such as bifurcations, limit cy- c material [2]. A long-standing challenge is to understand cles, solitary waves, and spatiotemporal chaos [11]. s - the pattern selection starting from the basic ingredients: A two-phase solidification front consists of (i) solid- l r bulk transport, solute and heat rejection on the solidifi- liquidinterfacesand(ii)trijunctionpointswhereallthree t m cation front, and the front’s local response. Simple as it phases meet. Our strategy is to construct a phase-field may seem,this free boundary problem(FBP) accurately modelthatallowsustoanalyzethethin-interfacebehav- . at describes many experimental features, but has few ana- ior of (i) separately from (ii). We quantitatively repro- m lyticsolutions,sothatnumericalmodelingismandatory. duce the correct FBP on (i); (ii) satisfy Young’s law at The phase-field method [3] has become the method of equilibrium. WetestconvergenceinW/ℓforlamellareu- - d choice for simulating solidification fronts [4], and more tecticgrowthatexperimentallyrelevantparameters,and n generally for tackling FBPs and interfacial pattern for- compareour results to boundary integral(BI) [12] simu- o mation phenomena, e.g. in materials science [5] and lations and other phase-field models. For steady states, c [ fluid flow [6]. Its main advantage (essential in three di- weachievegoodagreementwiththe BI anda drastically mensions) is that it circumvents front tracking by using improved, fast convergence compared to previous mod- 2 phase fields to locate the fronts. These fields interpo- els. In contrast, convergence is slow for limit cycles, due v 9 late betweendifferentconstantvaluesineachbulk phase to a trijunction behavior affecting the overall dynamics. 4 through interfacial regions of thickness W. The model We use one phase fieldp to indicate presence(p =1) i i 4 is then required to reproduce the FBP in the sharp- or absence(p =0)of eachphase i=α,β,L in the spirit i 3 interface limit, in which the extra length scale W van- of volume fractions [13], which requires 0 ishes. 3 p +p +p =1. (1) 0 Inpractice,simulationshavetoresolvethevariationof α β L / the phase fields through the interfaces, so that W must t Thephase fieldsevolveintime to minimize afreeenergy a stay finite. Their results generally depend on the ratio functional ofp~ (p ,p ,p ),thesoluteconcentration, m α β L W/ℓ, where ℓ is a relevant length scale of the FBP. Ex- F ≡ and temperature, - plicitcorrectionstotheoriginalFBPtofirstorderinW/ℓ d have been calculated by a so-called thin-interface analy- ∂p 1 δ n i = F i, (2) co sciosminpleatefewcanccaeslelsa,tiaonnd, asocmhieevceadncfoerledsinogulte-[p6h,a7s,e8s,o9li]d.ifiA- ∂t −τ(p~) δpi(cid:12)(cid:12)(cid:12)pα+pβ+pL=1 ∀ v: cation [7, 9], means that results become independent of where τ(p~) is a phase-depend(cid:12)ent relaxation time. This i W/ℓforsomefinite valueofW. ThecorrectFBPisthen classicalproblem ofminimizing a functional subject to a X reproduced already at that value, much larger than the constraint is treated by the method of Lagrange multi- r thickness ofrealinterfaces,enabling quantitative contact pliers;(δ /δp ) =δ /δp (1/3) δ /δp a F i |pα+pβ+pL=1 F i− j F j inthreedimensions betweensimulations,theory,andex- for three phases, where the functional derivatPives on the periments in reasonable simulation times [10]. r.h.s. are now taken as if all p were independent. i Here,weextendtheseadvancestotwo-phasesolidifica- Todistinguishbetweenphases,earlierphase-fieldmod- tion, which already includes the most widespread solid- elsoftwo-phasesolidificationusedeithertheusualsolid– ification microstructures after dendrites: eutectic com- liquid phase field and the local concentration [14] or in- posites. They consist of alternate lamellae of two solids troduced a second, α–β phase field [15]. Across a solid– (α and β) or of rods of one solid embedded in the other, liquid interface, both fields must vary, so that their dy- growing from a melt L near a eutectic point, where all namics are coupled, which complicates a thin-interface three phases coexist at equilibrium. The interplay be- analysis. The same is true for a generic choice of in F 2 Eq. (2). However, if on an i–j interface we can assure concentration in phase i coexisting with phase j. A eu- that the third phase field p is exactly zero, p or p can tectic phase diagram with constant concentration gaps k i j be eliminated using Eq. (1), so that the interface can be and straight liquidus and solidus lines is generated by describedintermsofasingleindependentvariable. This A =c c(C ) and B =c (T T )/(m ∆C), with m i i i i i E i i ≡ − was recently achieved using a free energy with cusp-like the (signed) liquidus slopes, i=α,β. Non-constant con- minima[16],butnothin-interfaceanalysisisavailablefor centration gaps and peritectic phase diagrams can also that model. We also achieve absence of the third phase, be treated. Without loss of generality, A =B =0. L L butusingasmoothfreeenergy,byrequiringp =0tobe In order for µ=µij to keep the balance all across the k eq a stable solutionfor p ofEqs. (2) for eachi–j interface: i–j interface as p goes from 0 to 1, we require k i g (p ,p ,0)=1 g (p ,p ,0) i. (8) i i j i j i − ∀ δ F =0 k, (3a) Otherwise, several thin-interface corrections arise [8, 9]. δp (cid:12) ∀ k(cid:12)(cid:12)pα+pβ+pL=1,pk=0 The simplest choice satisfying also Eq. (3a) is gi = δ2 (cid:12) p2 15(1 p )[1+p (p p )2]+p (9p2 5) /4. F >0 k. (3b) i{ − i i− k− j i i − } δp2 (cid:12) ∀ The evolution of µ is obtained from its definition and k (cid:12)(cid:12)pα+pβ+pL=1,pk=0 massconservation,∂ c+~ J~=0, J~= Dp ~µ+J~ : (cid:12) t ∇· − L∇ AT The advantage is that the simplest choice for yields a F model that turns out to coincide with the quantitative ∂µ ∂h model of Ref. [9] on those i–j interfaces. =D~ pL~µ Ai i ~ J~AT, (9) To construct our free energy, we split it into parts, ∂t ∇·(cid:16) ∇ (cid:17)−Xi ∂t −∇· where Dp ~µistheusualdiffusioncurrent,withadif- = f +f +λ˜f . (4) − L∇ grad TW c fusivity that varies from D in the liquid to 0 in the solid F Z V (one-sided model), and J~ is an extension of the anti- AT The first is a free energy penalty trapping current introduced in [9] that counterbalances spurious solute trapping, W2 2 fgrad = 2 Xi (cid:12)(cid:12)∇~pi(cid:12)(cid:12) (5) J~AT ≡−nˆL2W√2 Ai∂∂pti(nî·nˆL), (10) (cid:12) (cid:12) i=Xα,β for the gradients of the phase fields that provides the interfacethicknessW. The nextis atriple-wellpotential where nî = ~pi/ ~pi are unit vectors normal to i–L −∇ |∇ | interfaces, and nˆ nˆ prevents solute exchange between i L · the two solids. The model is not variational, because of fTW =Xi p2i (1−pi)2 (6) tuhseether=mpJ~A,Twahnicdhbaelcloawussefoµr6=a ∂cofac/rs∂ecr,dbiusctreentiazbalteisonus[7t]o. i i Ourmodel[Eqs. (2)and(9)]hasstable interfacesolu- that generates the basic “landscape”: one well per pure tions connecting two coexisting phases i and j: µ=µij, phase and “valleys” with double-well profiles along each eq pk = 0 cut, separated by a potential barrier on tri- pi = 1−pj = {1±tanh[r/(W√2)]}/2 (with r the dis- junctions pα = pβ = pL = 1/3. The last part has a tance to the interface), pk =0. Since these solutions are strength λ˜ (a constant that controls convergence) and identical for all i-j pairs, so are the i-j surface tensions. Unequalsurfacetensionscanbe obtainedby addingnew couples the phase fields p to the temperature T and the i terms in Eq. (4) that shift the i–j free energy barriers. solute concentration C through c(C) (C C )/∆C, E ≡ − Remarkably,onsolid–liquid(i–L)interfaces,assuming where ∆C C C , C and C are the limits of the β α α β ≡ − aweakdependenceoftheA ,B onT,andτ(p~)=τ ,the eutectic plateau, and (C ,T ) is the eutectic point, i i i E E change of variables φ =p p , u=(µiL µ)/A maps i i− L eq − i Eqs. (2)and(9)tothequantitativemodelwithconstant f = g (p~)[B (T) µA (T)], (7) c i i i Xi − concentrationgapin[9],uptonumericalprefactors. The thin-interface limit can hence be deduced by inspection where we have introduced the chemical-potential-like and yields the classic FBP on i–L interfaces, variable µ c A (T)h , and g (p~) and h (p~) (given below) inte≡rpol−atPe bietwi een 0i for pii=0 and1ifor pi =1. ∂tc = D∇2c, (11a) The term fc drives the system out of equilibrium by −Dnî·∇~c = vn(ciiL−ciLL), (11b) unbalancing the pure phase free energies: Each well i is T T E µshi=fteµdiejqby=an(Bajm−ouBnit)/B(Ai −j −µAAii.) Tgihveeseqeuquilaiblrsiuhmiftsvaalnude c = ∓(cid:18)|m−i|∆C +diκ+βivn(cid:19), (11c) hence restores the balance between phases i and j; from where Eq. (11a) holds in the liquid and the others are the definition of µ, we obtain cij = A + µij for the boundary conditions on the interface that has normal i i eq 3 velocityvnandcurvatureκ;theminus(plus)referstoi= (a) α(β),andthecapillarylengthsd andkineticcoefficients i -0.003 β read in terms of our model parameters i W d = a , (12) i 1|Ai|λ˜ )/lT τi Ai W Vt βi = a1(cid:20) Ai λ˜W −a2| D| (cid:21), (13) (z- | | witha =√2/3anda =1.175. Theconstantλ˜ W/d -0.0035 1 2 i ∝ inEqs.(4),(12)and(13)controlstheconvergencetothe original FBP. Any set of β can be treated with suitable i boundary integral τi. Weconsiderhereβα =ββ =0,whichisachievedwith qualitative model τ = a A2λ˜W2/D. The different τ for A = A (e.g. -0.0026 (b) present without antitrapping i 2 i i α 6 β present with antitrapping different concentration gaps) are interpolated by τ(p~) = τ¯+(1/2)(τ τ )(p p )/(p +p ),τ(p +p =0)=τ¯, α β α β α β α β − − with τ¯=(τ +τ )/2. T α β /l We test our model in directional solidification with Vt) >/lT0.003 Tan=d VTE>+0G,(tzh−eVput)ll,inwghesrpeeeGd,>b0otihstdhieretchteerdmaalolnggratdhieenzt (z- z-Vt) <( axis. Halfaeutecticlamellaepairoftotalwidthλissim- -0.0028 ulatedintwodimensions(xandz)withno-fluxboundary -0.0031 0 32 λ6/4W 96 128 conditions in the midline of each lamella, using a finite- 0 x/λ 1 difference Euler scheme with a grid spacing ∆x = 0.8W (coarser far into the liquid to improve efficiency). We FIG. 1: Steady-state lamellae pair profiles (dimensionless adopt l /d¯= 51200 and ¯l /l = 4, where l D/V D T D D undercooling vs. x/λ) for different models. Four curves at ≡ is the diffusion length, lTi ≡ |mi|∆c/G are the thermal λ/W =32, 64, 96 and 128 shown per model; curves closer to lengths,andd¯ (d +d )/2,¯l (lα+lβ)/2. Thesecor- theboundary integral: larger λ/W. [λ/W = 64–128 collapse respond to typ≡icalαexpeβrimenTta≡l valTuesTG 100K/cm, forthepresentmodelwithantitrappingcurrentin(a)]. Phase V 1µm/s for CBr -C Cl , an organic≈eutectic for diagram used: (a) symmetric; (b) close to CBr4–C2Cl6. See 4 2 6 ≈ parameters in the text. Inset: Averaged undercooling in (b) which accurate experimental data exist [11]. We use vs. λ/W,compared to that without antitrappingcurrent. m = m , c = c (a symmetric phase diagram) α β α β − − or m /m = 2, c /c = d /d = 2.5 (one close to β α β α α β − − CBr4-C2Cl6). In both cases µ(z + ) = 0 (eutectic [18]; in this situation, several thin-interface corrections → ∞ composition). We test convergence to the thin-interface to the FBP occur simultaneously [8, 9]. limit with decreasing W by conversely increasing λ/W Results are similar for the phase diagram close to whilekeepingalltheratiosaboveandλ/λminfixed,where CBr4-C2Cl6 [Fig. 1(b)]. The convergence is somewhat λmin d¯lD is the minimal undercooling spacing [17]. slower, since one of the lamellae is thinner and needs to This i∝s apchieved by varying the constant λ˜ in Eq. (12). be properly resolved. Some small deviation from the BI Figure 1 shows the solid–liquid interfaces of a steady- persists, probably due to the trijunction behavior (see state lamellae pair calculated by different phase-field below). In the inset, we plot the average undercooling models and the boundary integral method (BI) [12] for vs. λ/W. This is a less stringent test, as shown by the λ λ . For the symmetric phase diagram [Fig. 1(a)], fact that results for our model are convergedalready for min ou≈rmodel(thinsolidlines)agreeswellwiththeBI(thick λ/W = 32. However, those for the model with J~ =~0 AT solid line). Moreover, the curves at λ/W = 64, 92 and still depend on λ/W at λ/W = 128, which illustrates 128 are indistinguishable. This means that the results how all corrections need to be canceled before quantita- are independent of λ/W for λ/W 64, the signature of tive results can be achieved. ≥ a quantitative model. In contrast,if we remove the anti- Next, we increase λ to 2.2λ , close above the min trapping current in our model, J~ =~0, which leads to threshold λ 2λ [12] for≈the bifurcation from steady AT min ≈ solute trapping and finite interface kinetics, the results lamellae to oscillatory limit cycles, a situation in which dependonλ/W foralltherangefrom32(bottomdashed the oscillation amplitude is very sensitive to all param- line) to 128 (top one). The convergence of models not eters. Indeed, for the symmetric phase diagram and backed by a thin-interface analysis can even be slower, λ/W = 64, the qualitative model of Ref. [18] still yields as shown by the dotted curves for a qualitative version lamellae, whereas the present model correctly produces of our model with h =g violating Eq. (8) and J~ =~0 cycles, which are shown in Fig. 2(a). However, the i i AT 4 (a) 0.14 (b) yield a substantial efficiency gain for small curvatures of trijunction trajectories, which makes it a promising λ tool for three-dimensional simulations. Second, the free A/ boundary problem to converge to should also be recon- sidered. It was shown elsewhere that Young’s condition λz/ 1 0.1 64 96 128 160 192 ontheanglesbetweeninterfacesisviolatedoutofequilib- λ/W rium for kinetically limited growth [19]; here, the global trijunction rotation was found to be fairly independent oftheinterfacethickness,sothatitmightpersistforreal nanometric interfaces. These effects should be further investigated, possibly by atomistic simulations. 0 (c) 0 x/λ 1 We thank S. Akamatsu and G. Faivre for discussions, FIG. 2: Limit cycles. (a) Superimposed snapshots of the A.KarmafortheBIcode,andCentreNationald’E´tudes interfaces at constant time intervals for λ/W = 64. Thicker Spatiales (France) for support. R. F. also acknowledges lines: α–β interfaces. (b) Amplitude of the trijunction oscil- a European Community Marie Curie Fellowship. lation in units of λ vs. λ/W. The line is a fit that yields A(λ/W →∞)/λ=0.142. (c)Blowup of6.4W×6.4W. Solid lines: trijunction passage; dashed line: later α–β interface. ∗ amplitude of the trijunction oscillation A/λ, defined as Presentaddress: UniversiteitLeiden,Postbus9506,2300 its maximal displacement in x/λ, strongly depends on RA Leiden, The Netherlands. [1] E.BenJacob andH.Levine,Adv.Phys.49,395(2000); λ/W, as shown in Fig. 2(b). An extrapolation yields W. J. Boettinger et al. ActaMater. 48, 42 (2000). A(λ/W ) = 0.142λ, not far from the BI result → ∞ [2] W.KurzandD.J.Fisher,FundamentalsofSolidification A = 0.139λ, but the results are still not converged for (Trans Tech, Aedermannsdorf, Switzerland, 1992). λ/W =192,instrongcontrasttothesteady-statebehav- [3] J. S. Langer, in Directions in Condensed Matter, edited ior. This suggeststhat some correction(s)to the FBP in by G. Grinstein and G. Mazenko (World Scientific, Sin- W/λ remain in our model. Since solid–liquid interfaces gapore, 1986), p.164; J.B.Collins andH.Levine,Phys. are controlled, we turn to the the trijunctions. Rev. B 31, 6119 (1985); G. Caginalp and P. Fife, Phys. Rev. B 33, 7792 (1986). Thesolid(dashed)linesinFig.2(c)showafirst(later) [4] W. J. Boettinger, J. A. Warren, C. Beckermann and A. snapshot of the interfaces close to a turning point of the Karma, Annu.Rev.Mater. Res. 32, 163 (2002). trijunctiontrajectory. Inthelateronethetrijunctionhas [5] L.-Q. Chen, Annu.Rev.Mater. Res. 32, 113 (2002). moved away and only the α–β interface remains, which [6] R.Folch,J.Casademunt,A.Hernández-MachadoandL. has slightly moved sideways. In the one-sided FBP, (i) Ram´ırez-Piscina,Phys.Rev.E60,1724(1999);60,1734 the α–β interface cannot move, so it is the trace left by (1999); T. Biben and C. Misbah, Eur. Phys. J. 29, 311 the trijunction, and (ii) its direction close to the trijunc- (2002); Phys.Rev.E 67, 031908 (2003). [7] A. Karma and W.-J. Rappel, Phys. Rev. E 53, R3017 tion approaches that of the trijunction velocity. In a (1996); 57, 4323 (1998). diffuse-interfacemodel,thediffusivitybehindthetrijunc- [8] R. F. Almgren, SIAMJ. Appl. Math. 59, 2086 (1999). tion point pα = pβ = pL = 1/3 falls to zero on the scale [9] A.Karma,Phys.Rev.Lett.87,115701(2001);R.Folch, of W, so that (i) and (ii) do not hold. We consistently A. Karma, M. Plapp and B. Echebarria (unpublished). observe the displacement to be a fraction of W fairly [10] J.Bragard,A.Karma,Y.H.Lee,andM.Plapp,Interface independent of λ/W, and the whole trijunction to be Science 10, 121 (2002). slightly rotated with respect to its velocity, features also [11] For instance M. Ginibre, S. Akamatsu, and G. Faivre, Phys. Rev.E 56, 780 (1997). observed for the steady state in Fig. 1b. This effect ex- [12] A. Karma and A. Sarkissian, Metall. Trans. A 27, 635 plains the remaining mismatch between phase-field and (1996); BI in Fig. 1b and the slow convergence of A/λ here. [13] I. Steinbach et al.,Physica D 94, 135 (1996). We have presented a phase-field model of two-phase [14] K. R. Elder, F. Drolet, J. M. Kosterlitz, and M. Grant, solidificationthatcoincides withthe best models to date Phys. Rev. Lett. 72, 677 (1994); A. Karma, Phys. Rev. [7,9]onsolid–liquidinterfaces,whosedynamicsarecom- E 49, 2245 (1994). pletelycontrolled. Thishasallowedustoidentifytherole [15] A. A. Wheeler, G. B. McFadden, and W. J. Boettinger, Proc. R. Soc. London, Ser. A 452, 495 (1996). of diffuse trijunctions in the convergence of the results. [16] H. Garcke, B. Nestler, and B. Stoth, SIAM Understanding their dynamics is both a fundamental is- J. Appl.Math. 60, 295 (1999). sueandaprerequisiteforafullyquantitativemodelingof [17] K.A.JacksonandJ.D.Hunt,Trans.Metall.Soc.AIME multiphase solidification: First, a thin-interface analysis 236, 1129 (1966). of the trijunction region in the phase-fieldmodel is lack- [18] R. Folch and M. Plapp, cond-mat/0206237 (2002). ing. Even so, our model is expected to be precise and [19] C. Caroli and C. Misbah, J. Phys. I 7, 1259 (1997).