Towards a physically consistent phase-ﬁeld model for alloy solidiﬁcation

: We summarise contributions made to the computational phase-ﬁeld modelling of alloy 1 solidiﬁcation from the University of Leeds spoke of the LiME project. We begin with a general 2 introduction to phase-ﬁeld, and then reference the numerical issues that arise from solution of the 3 model, before detailing each contribution to the modelling itself. These latter contributions range 4 from controlling and developing interface-width independent modelling; controlling morphology 5 in both single and multiphase settings; generalising from single to multi-phase models; and 6 creating a thermodynamic consistent framework for modelling entropy ﬂow and thereby postulate 7 a temperature ﬁeld consistent with the concepts of, and applicable in, multiphase and density- 8 dependent settings. 9

. A path down a hill side following the path of steepest descent. The descending particle moves more quickly when the slope is steep (towards the top in this case) but slows down as it approaches a local minimum (near the bottom in this case).
to give insight into the role of free energy minimisation in the resulting morphology of  (1) or, in component form, Here x 1 , x 2 are the two Cartesian directions in space, (more conventionally denoted x, y 46 and (z for 3D), which we will use interchangeably with x i where convenient);ẋ is the 47 velocity; h is the height in coordinates x 1 , x 2 ; and the parameter, M > 0, which may be a 48 function of x 1 and x 2 , controls the response to the driving gradient ∇h. Also notice the 49 negative sign that enforces reduction rather than an increase in height. Such a descent 50 is illustrated in Figure 1 following the dashed line which is always perpendicular to 51 contours of equal height. The phase equation for the evolution of φ is very similar in form: where M controls the speed of descent -the mobility. F is the free energy of the domain 54 as a whole. In 3D this is 55 F = f (φ, c) dx 1 dx 2 dx 3 (4) where the free energy density, f , is given in general by gradients and functions of φ and 56 c (and also T if necessary). Variational techniques (using the notation discussed in [5]) 57 convert the variational notation into functions of phase, φ, and its gradients, ∇φ: In Eq. 1, the path has two degrees of freedom. On the other hand, since the function, is used for solute evolution (typically, gradients of solute, c, are not used in the formation 73 of free energy). Here D is the solute diffusion which is always (for substitutional 74 diffusion) several orders smaller in the solid than in the liquid, and written (in the two

77
Once the free energy density is specified, the evolution equations are known given 78 initial conditions. For binary alloy models the phase-field initial condition consists of an 79 initial nucleus, without which the system would remain in a permanent metastable state 80 of undercooled liquid. This nucleus must be the solid phase (not just a surface such as 81 a boundary or another phase). The initial condition for solute consists of a boundary 82 value and often a value which differs in the solid region.

83
The evolution equations can rarely be solved analytically, even in 1D. Considerable can be prohibitive, see [6] for a detailed discussion of numerical techniques brought to 92 bear on model of [3].

93
There are a number of numerical techniques in operation in the solution of such 94 equations. One in particular is the treatment of non-linearity, see [7]. This typically, 95 involves linearising the system locally, using the derivative of terms in the evolution 96 equations. This can lead to complicated analytical expression which are model depen-97 dent. That is, these terms will often change when a new material is introduced, say, from 98 AlCu to AlFe -or even changing from a physically realistic data based model to a more possibly thousands of nodes, becomes an issue, see [8]. Note that the contributions of 104 [7,8], whilst using phase-field simulation as the motivation, are primarily around our 105 development and understanding of numerical methods that can be applied to a much 106 wider class of time-dependent problems -see [9].

108
For simplicity of exposition we focus now on systems with only two phases present.

109
Multiphase will be revisited briefly in Sec. 6 and subsequently Sec 9.

110
Formation of the free energy density, f , is where much of the modelling takes place.

111
The free energy density is a sum of two parts: the surface, f S , and the bulk, f B : For the simplest binary alloy model, f S is given by where δ has the units of length and is associated with the thickness of the transition  Considering just the surface contribution in 1D so that then setting the variational derivative of this to zero, Eq. 12 is solved by the phase profile which is near zero for x 0 and near unity for x 0 with a transition region in 129 proportion to δ about x = 0. The profile is illustrated in Figure 2 with δ = 1, and means, 130 in the absence of any driving force the shape of the phase-field satisfying just the surface 131 contributions to the phase-field, is given by Eq. 13. This may prompt the question: How 132 is the phase profile controlled in general, say to match a specific profile? Eq. 13 may be By considering a generalisation of Eq. 14 in the form of it becomes possible to use the coefficients a 1 , a 2 and a 3 to fit a desired profile to data 136 from, say, a molecular dynamics simulation or ab initio calculation. As shown in [10], the 137 form of the potential well, Ω, necessary for a skewed profile is given by and the form for the interpolation polynomial, g, is given by Derivation of these Equations is given in [10], where an application of the method to real 140 data is shown in Figure 3.

142
Anisotropy is typically introduced into phase-field models by the addition of an 143 anisotropy term, A(∇φ) in the following manner The anisotropy, A, is dimensionless, but because of the presence of gradients in its form of A is that of the 2D four-fold anisotropy (not withstanding its limited use in 149 nature, since 2D systems are rare): where the new variable, θ, is related to the phase variable where φ x , φ y is shorthand for the x and y derivative of φ respectively. The effect of    The role of the parameter, , in the anisotropy is interesting. When = 0 we 160 return isotropy. However, there is a value, = 1/15 in Eq. 19, beyond which the 161 resulting predicted equilibrium morphology in unstable. At the critical value, the 162 function develops points of inflection (see Figure 5). Now, it is widely known (and our 163 paper, [12], discusses this ) that a point of inflection in the anisotropy is associated with 164 a discontinunuity in the associated equilibrium shape (also called the Wulff shape). The 165 equilibrium shape, W, is given parametrically by and is illustrated with three values of in Figure 6 167 This has not prevented some researchers from using values of 1/15, with some 168 stabilisation techniques in order to simulate faceted growth, [13][14][15]. After first following 169 standard literature to explore the growth of a regular hexagon in [16] to investigate the 170 transition from pure faceted to dendritic morphology (see Figure 7), we argue, in our 171 paper [12] that there is a cleaner, more flexible approach to facet modelling. To introduce 172 faceting into the model we first specify the vertices and give the following formula for 173 anisotropy: where p i are the n vertices of the resulting required equilibrium morphology. The (see [17]) where ∆ i , a i are given constants. A traditional method for combining these 189 separate energies into a bulk free energy, f B , is given by where, assuming a suitable interpolation function (i.e. monotonically increasing between 191 g(0) = 0 and g(1) = 1, and with zero gradients at the extremes, g (0) = g (1) = 0), prescribes a set of normalized functions (∑ i g i = 1), for multiphase interpolation.

193
For two phase, the formulation of the bulk free energy reads   Figure 10.

209
This problem was largely solved by [18] and later put on a sound theoretical basis 210 by [20]. The method of [18] as illustrated in Figure 11 is most readily adopted when 211 quadratic approximation is used, and is the method we adopt in [12] when simulating 212 faceted growth. The most distinguishing feature of the method of [18], illustrated in 213 Figure 11 is that all the intervening curves have a common tangent value, i.e. there 214 is a line which is tangent to all curves. Noticing this latter feature of the [21] model, 215 we suggested in [19], an alternative method for forming the intervening curves, which 216 involves a straightforward linear translation of the curves illustrated in Figure 12. There 217 are two advantages of the latter method, as compared with [18]: one, it readily applies 218 to non-quadratic approximations; two, the intervening curves do not extend outside of 219 the unit range, c = 0 to 1 if the two bulk curves do not extend outside this range. For 220 method [19], the bulk free energy is given by: Figure 11. As in Figure 10 but with a different interpolation scheme. The transition between these two pure phases (in brown, as we move across the boundary between the two) is shown with the intervening curves. This is the method advocated by [18] for interpolating the two curves The transition between these two pure phases (in brown, as we move across the boundary between the two) is shown with the intervening curves. This is an alternative method suggested by [19] for interpolating the two curves That is, we use the equilibrium common tangent values c L , c S to specify the translations.

223
Although this method does transition between solid and liquid smoothly between the 224 original curves, it implies a very rapid change in the 2nd derivative near the liquid 225 value, φ = 1. This feature can be compensated for, but is (naturally) not present in the 226 method of [18] illustrated in Figure 11 and, together with its theoretical underpinning 227 [20], makes this the method of choice where possible. The previous section concludes that [18] is the method of choice for constructing 230 the bulk free energy for intermetallic alloys. For less extreme free energy functions of 231 fully soluble solids (i.e. not intermetallic), either the methods of [18] or [17] can be used 232 successfully. That said, both methods still present an unresolved problem, for which 233 this section puts forward a solution, applicable for both (and any) interpolation method.

234
A key parameter in phase-field modelling is δ, the interface width. The value of δ may 235 be given a physically real value, possibly arising from a molecular dynamic simulation.
where Ω is the double well potential, e.g. Ω = 1 2 φ 2 (1 − φ) 2 , and derive equations of and where f c ≡ ∂ f ∂c . We now assume that we have a phase profile of the form that solves the 1D equation, Eq. 30 for some V and δ. Note that the slope of the

259
For two convex curves, there is a unique point on each curve such that So, by choosing c S as one boundary value, we find, as V decreases that the extreme  Using the value of V and δ at the tip of a growing dendrite in 2D or 3D, together with 272 β, then allows us define a compensating current for use in two and three dimensions: to be used with a modified solute equation: where the square matrix, M is given as a multiple of This guarantees that and so with consistent initial conditions, Eq. 42 is satisfied for all time, t. Unfortunately, 295 such a mobility matrix, Eq. 44, and its generalisation to n > 3 fails to reduce properly when only two of the n > 2 phases are present. A spurious phase will likely be generated.

297
This is because, for example, with φ 1 = 0: the contribution from the second and third terms can be non zero. An alternative 299 formulation, that we propose in [24], is to specify the mobility matrix as phase dependent and To understand this form of matrix it helps to note that when only two phases are present ). An illustration of the the mobility matrix is given in Figure 16 (model B -a more 305 complicated, but similar, model -is also discussed in [24]).

306
For single phase the formation of the surface energy density for isotropic growth is 307 given by Eq. 18. So, it is natural to postulate a multiphase surface energy density by the 308 equation postulated in [23]: where, indeed, this form reduces correctly when two phases only are present. As energies, σ ij ∝ δ ij , to be included, and also consistent with two phases, was proposed in 316 [25] and [26] involving a double sum: It is not clear that the formulation Eq. 50 does actually control the angle at the junction.

318
Some formulations introduce angle control by introducing an extra potential in the form 319 of (for three phases), φ 1 φ 2 φ 2 or its square, [27]. In [28] we propose a new formulation 320 which gives direct control of the angle at the triple junction: where as a consequence of a definition the single indexed quantities may be written as a function of the double indexed quanti-323 ties by Here, δ ij , i = j is the surface energy between phase φ i and φ j . The relation of the surface 325 energy to the angle at the triple junction is shown in Figure 17, which shows three 326 lines containing three desired angles and then a triangle such that each side, n ij crosses 327 the lines at right angles. The relative lengths of the triangle sides is proportional to 328 the quantities, δ ij . To get a feel for the formulation set the two solid interface widths 329 δ 2 = δ 3 = 1, then when δ 1 (the liquid) is equal to unity the formulation Eq. 49 is returned, 330 but when δ 1 > 1 the angle in the liquid area begins to increase past 120 • , with the other 331 two equal regions decreasing accordingly. Using a surface energy defined by Eq. 51 in 332 phase-field gives the desired angle control as shown in Figure 18. other postulated models as shown in Figure 19, where the surface energy at the centre 337 of the simplex, φ i = 1/3, i = 1, 2, 3, is lower than at the binary junctions on the left but, 338 conversely, higher for the model advocated in our paper, [28]. The other feature present  The centre of the simplex is formally the triple junction point where, ideally, the surface energy barrier should be highest. The left plot shows the models advocated by [25] and [26] containing a hollow at the junction, whereas the model proposed by [28] (and also [23] in the special case of equal angles) has the formally correct barrier at the triple point as well as the binary junctions.
in our model (the right of Figure 19) is that the gradient of surface energy is zero at the 340 boundary of the simplex, but this is not the case in the models advocated in [25] and 341 [26]. A low triple junction barrier generates a spurious phase at a binary junction. The 342 model advocated in [28] avoids both spurious phase generation and also, due to the 343 zero gradient, does not give preference to a value outside the simplex, i.e. non-physical 344 negative values of φ i . role in alloy solidification such as growth restriction [29], but also when the partition 352 coefficient, k E → 1 (e.g. for almost complete solute trapping at very high growth rate), 353 whence thermal forces exceed diffusion effects. One way of processing this is to impose 354 a temperature field on the other fields, but this does not take into account the heat 355 generation due to the volume latent heat of fusion.

356
In the early 2000s fully coupled single phase formulations were in existence, [3,6], where C v is heat capacity at constant volume, k conductivity and L the latent heat.

364
The term Lφ generates heat as the material changes phase. Now consider a plausible 365 multiphase generalisation given by This fails because ∑ n iφ i = 0. A compromise might be to single out the liquid phase as 367 special and write but this undermines an implicit assumption in phase-field modelling that no phase is see [30]. These models and the concern that the resulting models may not be thermody-377 namically consistent was addressed in [5], which proposed a formalism for modelling 378 complex fluids that closely resembled phase-field models of solidification. That is, to and multiphase phase models. The formalism also shows how to recover the entropy 384 from which it is a short step to generating a temperature equation. Then it becomes 385 clear that the only additional parameter required to specify a temperature equation is 386 the conductivity, k. Not only is the heat capacity and latent heat prescribed (from the 387 free energy specification), but the addition of a term involvingċ, similar to latent heat of 388 fusion, also falls out of the formalism, and treatment of density change together with 389 necessary fluid flow also is prescribed.

390
The n-phase formulation for the temperature field of a binary alloy, postulated in

391
[31] is given as: where D is the solute diffusivity. If we neglect the contributions due to solute change 397 and gradients we have This, with hindsight, seems a reasonable guess, but Eq. 59 is required to complete the 399 definition, not only making the latent heat terms, L i , non-constant but even differing in 400 sign.

401
Finally, the general bracket formalism also poses the question as to why there 402 is no further coupling between phase and solute concentration as would naturally 403 occur if solute, c and the phase-fields were treated equally as part of a general vector 404 φ = [c, φ 1 , φ 2 , φ 3 ], something that has been noted and used by other authors such as 405 [32]. Figure 20 shows the general scheme for constructing evolution equations, that 406 we described in [31] . The key novel part of which is that one constructs the bracket 407 alongside the free energy. The bracket includes the diffusion parameters, which fields 408 are conserved or not, and draws a distinction between energy conserving fields such 409 as mass, and the dissipative dynamics as is the staple of relaxation phenomena such as

413
Here we summarise our main contributions to phase-field computational modelling 414 1.
Non-linear adaptive mesh multi-grid solver for multi-scale problems [ phase-field formulation with correct thermal field and entropy generation [31] 427