Towards a quantitative phase-field model of two-phase solidification

0.49 MB·English

by  R. Folch

Towards a quantitative phase-field model of two-phase solidification
∗
R. Folch and M. Plapp
Laboratoire de Physique de la Mati`ere Condens´ee,
CNRS/E´cole 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ˆi·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ˆi = ~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ˆi·∇~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´andez-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).