Multiphase field modelling of alloy solidification

Computational Materials Science 171 (2020) 109085

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Multiphase field modelling of alloy solidification

T

P.C. Bollada , P.K. Jimack, A.M. Mullis ⁎

University of Leeds, United Kingdom

ARTICLE INFO

ABSTRACT

Keywords: Multiphase field Alloy solidification Crystal growth

We present an approach to alloy solidification modelling that incorporates binary interface energies in a manner that correctly reproduces the associated theoretical angles at triple junctions in eutectic solidification. We find that simply applying the principle that the correct binary junction behaviour is recovered when only two phases are present is insufficient. Previous research (Toth, 2016) recommends a modification of the surface energy by adding an energy barrier at the triple junction, and we explore alternative models that would benefit from this approach. The main approach we recommend here, though, is to extend the minimal model of Folch and Plapp (2003, 2005), which, without modification, is limited to 120° junction angles. This is achieved by a linear transformation of this formulation, and facilitated by an analytical multiphase solution presented here for the first time.

1. Notation

preserves the mass of each species and introduces the solute diffusivity, ) D liq + Dsol]/(RT ) D, which is usually of the form D = c (1 c )[(1 D liq , where R is the molar gas constant. with Dsol The most complete temperature formulation was established in [5], with the simplest approximation to formulation found there being

We begin with a short section to establish the notation. 1.1. Single phase formulation and notation A phase field model for binary alloy solidification uses three field variables: the phase, ; the concentration, c; and temperature, T. The [0, 1], where phase field is designed to vary smoothly in the region = 0 is liquid and = 1 is solid, and similarly the alloy component of one species varies in the range c [0, 1]. Values of away from the bulk extremes indicates a transition region lying in a smoothed out representation of the interface. Once the variables are given, the free energy, F in terms of a free energy density, f, is constructed. This typically does not contain gradients of c or T, but it must contain gradients of to keep the interface smooth:

F=

f(

, , c, T )dxdy

(1)

where we choose without loss of generality a two dimensional domain, . Evolution of is then prescribed by

=

M

F

,

(2)

where the mobility, M, is prescribed, for example as in [4], where it may be a function of c. The concentration equation

c=

⁎

·D

F c

(3)

Cp T =

·

(4)

T+L

with latent heat of fusion, L, heat capacity, Cp , and thermal conductivity, , are constant. In this form the heat generated by solidification ( > 0 ) is clearly seen. The free energy density may be split into two contributions: surface free energy density, fS ( , ) (and perhaps also a function of c and T), and bulk free energy density, fB ( , c, T ) : (5)

f = fS + fB .

Crucially, the surface free energy density, fS , is a function of gradients of and is formed such that this term vanishes in the bulk ( = 0, 1). The formulation of fS given in [4] is of the form

fS =

1 2

2

·

+W

2 (1

)2

(6)

where = 6 , W = 12 / , is the liquid-solid interfacial energy, and 1 thus 2 / W = 2 2 . Here has the dimensions of length and can be taken as the interface width. The bulk energy, fB , can be formed by combining free energy densities for each phase, fliq (c, T ) and fsol (c, T ) , with the phase, , serving as a weight to interpolate between the two as will be discussed in more detail in this paper. fliq and fsol may be obtained from a database or 2

Corresponding author. E-mail address: [email protected] (P.C. Bollada).

https://doi.org/10.1016/j.commatsci.2019.109085 Received 11 March 2019; Received in revised form 13 June 2019; Accepted 14 June 2019 0927-0256/ © 2019 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/BY/4.0/).

Computational Materials Science 171 (2020) 109085

P.C. Bollada, et al.

· . As we will see, there are a number (companion) gradient term, of attempts to generalise this term, but, unlike the double well potentials on offer, there would appear to be no way of comparing them. We introduce a way of visualising and comparing rival gradient terms over a simplex in this paper.

Table 1 The notational differences between single and multiphase quantities. The single index is associated with the bulk and the double indices are associated with the interface. Physical quantity

Single phase

Phase

Mobility

Multi phase i,

M

Energy per length

2

Energy height

2 ij ,

W

Interface width

Interfacial energy Bulk Energies

fliq , fsol

Surface energy density

fS (

i

Mij , i, j

[1, n]

i, j

[1, n]

Wij, i, j

[1, n]

ij,

i, j

[1, n]

ij ,

i, j

[1, n]

fi , i

There are two main ingredients to modelling alloy solidification under the principle of energy minimisation: the diffusion parameters; and the free energy specification. Under single phase growth the mobility is the diffusion parameter associated with phase change, and, for alloys, there is an additional solute diffusion parameter associated with solute diffusion. For multiphase and/or multi-component systems these diffusion parameters become matrices. But here is the central problem in generalising from single phase to multiphase or multi-component: in single phase one variable represents a combination of two states. That [0, 1] has at its extreme values, two is, a single phase variable, states of matter (the same issue arises with the generalisation of binary alloy concentration, c). So the question arises: do we adopt variables for each bulk phase, i , or instead choose variables representing transitions between two phases of matter, ij . The latter would appear to be more closely associated with the single phase modelling, and yet the former, i , is almost universally adopted. This causes a mismatch between physically derived quantities such as surface energies and mobilities, which are intrinsically binary in origin, and the chosen single index variables of the model, e.g. i . The most reliable guide to forming a consistent multiphase model is that when we set 1 + 2 = 1 the system reproduces the single phase model. Inevitably, there are a number of possible multiphase models that have this feature. Let us assume that given Mij and Wij a multiphase model reduces correctly to the single phase model, then this implies that the equations for i and j cannot depend on any other indices except i , j . That is, for (i, j) = (1, 2) there can only be dependency on 1, 2 and also W12 (it can, however, depend on M11, M22, M12 and M21). An ideal multiphase model has a double index dependence (or more) in the surface energy and single index dependence in bulk behaviour. Thus given database energy density functions, fi (c, T ) for each phase, i , we might form the bulk energy as

[1, n]

fS ( ¯ , ¯ ), where ¯ = [ 1,

, )

3. Background

[1, n]

2,

…,

n]

CALPHAD calculation. 1.2. To multiphase notation To generalise the single phase formulation to an n-phase multiphase system involves defining the slightly richer notation in Table 1: All double index quantities are symmetric and there is a constraint on the n phases, i , such that i = 1 i = 1. Other multiphase models which do not impose this constraint will be considered in Section 4.2. 2. Motivation We now present the current state of multiphase modelling of alloy solidification in an attempt to extract a definitive model. We build upon the prior work of [6] which lays out a series of formal constraints on the form such a model must have. For example, clearly one measure of success in any postulated multiphase model is that it must reduce to an accepted single phase model when only two phases (the solid and the melt) are present in some spatial region. Similarly [1] investigates models with respect to the reproduction of the theoretically correct angles associated with interface energy terms, Wij (as well as aspects of flow, which are outside the scope of this paper). The model advocated by [6,1] incorporates interface energy into the model by creating scalar (non indexed) functions, W, of Wij (and also with 2 and ij2 ) and the phases, i in the following manner

fB =

W ( ¯) =

2 2 i j

(7) n j=2

j 1 i = 1 ),

(here i < j means the double sum, where it can be verified that when, for example, 1 + 2 = 1 we have W = W12 . In addition, for n = 3 phases, the energy is supplemented by a function that vanishes in the bulk and on a binary interface, h¯ = | 1 || 2 || 3 |. It is this combination that appears to allow energy preference to select the particular angle at triple junctions associated with the three binary values, Wij . We observe, in addition to the expression Eq. (7) not being well defined in the bulk (e.g. 1 = 1), that there may be more computation than is perhaps necessary in their proposed formulation, following the variational procedure:

F i

=

·

f

+ i

gi ( ¯) fi

fB = i

(10)

where gi , such that gi = 1, must be a function of all phases, i.e. they must be properly constructed weights. However, there are other approaches, for example, [7,8]. The latent heat of fusion is also naturally a double index quantity, because it is associated with a transition between two phases (even potentially solid-solid). In single phase models, heat generation is associated with phase change and a term (see, for example, [9] for pure metal solidification)

f i

(9)

where g is some monotonic function that interpolates between 0 and 1. Yet, even in this relatively simple term, Eq. (9) is incorrect: g ( i ) cannot be set equal to i since this results in a driving force independent of phase. Also it cannot be any other function because g ( i) 1 in general. A correct interpolation could read

i
i
g ( i ) fi i

Wij i2 j2

(8)

L

when f contains this term (this formulation will be stated later in this paper). How much weight to give to h¯ in the formulation and what the generalisation is to n phases is also not rigorously established. There is another strong motive to our paper, and this is that although much effort has been made in seeking a proper multiphase ) 2 , used in single generalisation of the quartic potential well, 2 (1 phase, there has been a less cohesive effort in generalising the

(11)

but, again, it is by no means obvious how a quantity, Lij can be incorporated into a multiphase model. Following a thermodynamically consistent approach developed by [10], and consistent with [11], this problem was resolved in [5], and, in fact, a single index quantity is advocated there for the heat source term. As stated in Eq. (2), the mobility parameter, M is incorporated into a 2

Computational Materials Science 171 (2020) 109085

P.C. Bollada, et al.

single phase model by

=

M

being slightly complicated also suffer from the awkward feature that they are not defined in the bulk (see also [37]). The scalar functions, W, apart from being slightly complicated also suffer from the awkward feature that they are not defined in the bulk. These two facts suggest that there may yet be a more elegant solution to this problem. Designing a potential term for more than two phases is at least readily analysable (and easily illustrated for three phases), so that one can see at a glance the changes in phase implied by the gradients. This is not the situation with the gradient energy, being a function of, at least, i and possibly also i . Perhaps because of this, rival formulations for this term have not been the subject of comparison or significant discussion. The approach of [6] to the potential term does lend itself to the postulation of the gradient term, and it certainly appears that this combination of potential and gradient energy terms is the most consistent in use today. That said, we believe there are other avenues of investigation that have not been tried. One of these is discussed here in Section 5, where we construct a multiphase solution analogous to the tanh profile 1D solution adopted in single phase modelling. Using this, we are able to inspect the gradient term as a function on the simplex (i.e. not only the potential term). The results are revealing. For example, we know that the potential term in [12] exhibits a lower triple junction energy, but, surprisingly, inspection of the gradient term used in [12] also has this feature. Also, using this visualisation, the model advocated by [6] does not appear to have the same energy shape as exhibited in their own potential. All the above seem to suggest that the optimum multiphase model has yet to be found. One way of discussing the candidate gradient energy models is to note that they are all quadratic in gradients and so fit the general form

F (12)

where F is the free energy functional. Initial attempts to generalise include i

=

M

F (13)

i

but it was observed that this will not guarantee the constraint i = 1 (see [12]), which lead to the introduction of a Lagrange multiplier approach. Looked at more generally the Lagrange multiplier approach is part of a family of formulations of the form i

=

Mij

F j

(14)

where, in the case of a Lagrange multiplier, the matrix Mij has constant entries and guarantees that i = 0 . This formulation correctly reduces to single phase, but [13] formally recognised and resolved, a fact well known in the modelling community, that this model exhibits n phase dependence. This is the disadvantageous modelling trait that on any njunction, other phases play a part in the dynamics (other that the n phases present). Perhaps unique in multiphase modelling is the model advocated by [14]. This model is n-phase independent and exhibits a model construction that incorporates double index quantities naturally. The formal objection to the model is that it does not construct a total free energy. The model was recognised at the time by its author as having shortcomings, and more recently [6] revealed that free energy minimisation was not guaranteed. It also does not lend itself to incorporating a temperature field by the methods of [5]. That said, there seems to be no formal argument that this model is incorrect, largely due to the non-physical nature of the interface region. Furthermore, this model has significant advantage computationally since it has fewer complex terms, and does not generate spurious phases. One critical feature of phase field modelling is the balance between gradient terms and potential terms that form the surface energy density. The premise is that a finite gradient in the transition region (phasephase interface) evolves due to the balance of these two terms. The double well potential term is readily visualised, as it is a function of phase, and has minima at the two extremes and a maximum at the junction. Thus this term attempts to create a steep gradient between the two phases. Typically a quartic function of phase is used as being the lowest order polynomial that has zero derivatives at the extremes. In passing, it may be observed that this function has a finite second derivative which also plays a part in the dynamics of solidification. Extending this function to three phases and more is not obvious. In fact, it is our view that this term is still not optimal in models to date despite recent improvements in understanding, [6]. The problem is to construct a potential that both increases away from an n-junction towards an (n + 1) -junction, and also incorporates the double index surface energies in a natural way. That is, the potential needs to provide a barrier to the formation of another phase from bulk; a further barrier to the formation of three phases from two, etc. The solution in [6] is to first assume equal surface energies and successfully find a polynomial that has all the right characteristics. They then combine the double index quantities into the single scalar, W, given in Eq. (7), which multiplies this function of i . The function W has the property that it equals the average of the binary heights at the triple junction which, in general, must be lower than one of the binary heights, Wij . This discrepancy can be compensated for, as shown in [1], by introducing a special triple junction term 1 2 3. Although it is no doubt correct, for generality, to incorporate a triple term, such a term may also be added to other models to achieve the same effect at less computational cost, e.g. the potential in [12], i < j Wij i2 j2 . The scalar functions, W, apart from

fG =

i·

ij

j

(15)

ij

where ij can be functions of, ¯ (i.e. i , i [1, n]). Applying the principle that the multiphase model reduces to single phase when n = 2 appears to narrow the models found in the literature to a handful of candidates. We can also use ij to postulate new models that necessarily reduce to single phase when n = 2 . This is compelling in itself, but we also find later in the paper, using a constructed three phase solution, further models which have extra desirable features - not least the ability to incorporate binary surface energies, ij . 4. Multiphase field modelling We discuss present approaches to multiphase modelling, for example, binary alloys with three distinct phases of matter represented by n = 1 where n = 3 in 1 (liquid), 2 (solid 1) and 3 (solid 2). Clearly i i this case. As in the single phase formulation, Eq. (1), the total free energy on the domain, is defined in terms of its density:

F=

f dxdy

(16)

where, without loss of generality we assume 2 spatial dimensions (2D) in two Cartesian coordinates x , y . Given Eq. (14), a form for Mij stems from the work described in [13]:

Mij =

M (c )

i j

(1

i )(1

j)

,i

j,

(17)

Mij (no implied summation with diagonal entries given by Mii = j i on the left). Eq. (17) is easily generalised to accommodate binary interface mobilities, mij by Mij =

M (c )

mij (1

i j

i )(1

j)

,i

j.

(18)

This formulation allows, for example, kinetic anisotropy to be included into a multiphase formulation by associating an i , j anisotropy, Aij with mij , but an exploration of this topic is beyond the scope of this paper. Associated with minimisation of free energy is the optimal increase 3

Computational Materials Science 171 (2020) 109085

P.C. Bollada, et al.

of entropy, which can in turn be related to temperature changes. In [10] the authors advocate the application of a dissipative bracket for modelling relaxation phenomena. By applying this approach to multiphase alloy solidification [5] gives the temperature field as

CT =

·j

1

T

F T

i

1

i

T

T

F c c

In [1] more work is done on the double well potential than in Eq. (25) by including a triple junction term which is designed to give a higher potential barrier for three phase coexistence. There is a perhaps more natural way to incorporate the double index quantity, Wij , into the formulation used in [28,12], given by

(19)

where heat capacity is related to the free energy density by

C=

2f T2

T

T+D

f c

f . c

·

T

1

T

T

i

i.

i

i

fB i

where

i=1

2 i (1

i)

2

.

2 j| ,

(28)

j.

fS = ·

fS

=

fS

k xa

xa

(29)

fS k

k

fS

a k

xa

(30)

k

a = 1. .d are the cartesian coordinates and a k , b x

ab k

,

(31)

a i

=

ab i

+ i

b i

,

(32)

to write

Typically, the bulk contribution, fB is constructed using a thermodynamic database [15] and using the phase variables, i , to interpolate [4,16]. Other methods of constructing a phase field model from database functions are found in [17] and put in context by [7,8]. For single phase alloy models the correction for interface width effects advocated in [18–20] has been generalised to multiphase, for example in [21] for dilute alloys, but is not readily generalisable to a full thermodynamic database approach. Thus one of the motivations for adopting the type of approach advocated in [17] is that there is a degree of interface width independence built in, so the generalisation to multiphase is more readily accomplished [22–26]. An alternative to [18] that builds on [4], allows computationally convenient interface widths, and the direct use of database is found in [27] though this is not yet generalised to multiphase. The surface energy, on the other hand, serves the double role of allowing the incorporation of surface quantities and controlling the finite interface between phases. It is safe to say that, at present, there is no concensus in modelling the surface free energy in a multiphase setting. As suggested in [6], the surface energy, fS , can be formed by inserting the surface terms, W, outside the interface controlling terms. Amongst the simplest ways of achieving this is to combine this idea with a form originating in [2,3], namely, as a single sum over phases:

+

(27)

We then use the chain rule

(24)

2 i|

2 i| .

j

i·

x a,

a k

The remaining modelling is done within the specification of the free energy as a function of phases, i , concentration, c, and temperature, T. The free energy is decomposed into separate surface, fS , and bulk, fB , contributions

1 | 2

j

i

k

4.1. Free energy construction

3

(26)

The gradient term Eq. (29) is more computationally convenient than Eq. (27) and appears to have no obvious drawbacks. To illustrate the potential complications in this term, we evaluate the variational derivative:

where the free energy density, f only contains the bulk contribution, fB . Of course there may well be scenarios where this approximation is not valid.

fS = W ( )

i

1 2

rij =

(23)

f = fS + fB .

2 2 i j)

for some combination of constants C. Another form for this term is introduced in [30]:

(22)

i

1 | 2

rij = C|

For computational simplicity all models we are aware of adopt an approximation in the last term of Eq. (22) which ignore terms in the free energy relating to the surface and use

F

j=1

More “natural” in the sense that the binary barrier heights are more directly associated with two phases. Other forms for rij are also available including [29],

Though formally present this term can be neglected in many modelling scenarios with the simplification to Eq. (19) being:

CT =

i=2

Wij ( 2rij +

rij =

(21)

F

i 1

where

(20)

and where, in general, the heat flux, j, has in addition to gradient of temperature, a term of entropic nature associated with the solute field:

j=

n

fS =

fS

fS

a k

xa

fS

a i

=

k

a k

i

+

ab i

fS

fS

a k

b i

.

(33)

k

Applying this expression to the formulation in [28] tends to give quite large expressions when compared with [30], i.e. Eqs. (26), (29) and also Eq. (25). To get a more general view of the possibilities for the surface energy term consider the more general quadratic form n

n

i

j

fG =

ij (

¯)

i·

j,

(34)

where ij is, in general, a function of 1, …, n . Using summation notation, and comma notation for derivatives of ij we find the useful expression for implementing a variety of models:

FG

=2

j·

ik, j

k

+2

2

ik

ij, k

i

i·

j

+2

ik

2

i

= (2

ik , j

ij, k )

i·

j

(35)

i.

Returning to the expression Eq. (34), a constraint on the coefficient ij is that the formulation should reduce to single phase with only two phases present, i.e. when n + m = 1. From Eq. (34) set (with summation on the left side only, i.e. not over n or m which are specific values)

fG

(25) 4

ij

i·

j

=

2

2

Wnm

n·

n.

(36)

Computational Materials Science 171 (2020) 109085

P.C. Bollada, et al.

Using

+

n

2

nn

+

nm

= 0 we have a constraint on the coefficients

m

2

=

mm

Wnm.

2

4. Introduced as a candidate here in Eq. (45) (37)

using Eq. (37), let us check the form used in [29] (adapted from Eq. (51) on page 173 in [29]) n

n

i

j

C 2Wij |

fG =

2 j|

i

C 2Wni,

C 2Wmi,2

mm=

i n

nm

2C 2Wnm

=

2C 2Wmn

2

+

nm

= 6CWnm + otherterms,

mm

nn=

2

i n

Wni i2

2

=

2

Wnm m2 , mm=

2 i m

2

Wmi i2

=

2

2

(49)

2 2 i j

2 (1

)2,

(40)

n=3 i 1 i=2

(50)

n=3 i 1

2 2 i j

Wij

<

j=1

Wij i=2

1 = 2 = 3 = 1/3

n=3 i 1

2

nn

+

nm

mm=

2

2

2

=

2

2 m

Wnm (

+2

Wnm

m n

+

2 n )=

2

2

m

+

i=2

2 n)

mm=0, nm

2

=

4

n

nn= i n

4

Wni i ,

nm

= 0, n

m,

1 + 12

n

2

j

2

i j

Wij

j

i·

n

i=2

1 Wij 2| 2

(45)

j=1

i

j

j

i

i 1

i=2

j=1

1 Wij 2

2

i·

j.

2 2 i j

i
i
2 2 i j

3 i=1

fDW = (46)

3

2

(52)

i)

2

i

.

(53)

1 = 2 = 1/2

1 2

2 2 i j

i
(54)

2 2 i j

1 + 12

2 2 i j

n i=1

4 i

3 i

4

3

+

1 2

2 2 i j i< j

. (55)

2 2 i j

i< j 2 2 i j

1 2

2 i (1

i)

2.

i

(56)

In Fig. 1 we present a generalised double well potential, fDW , that retains zero derivatives normal to the boundary as well as allowing varying barrier heights and quartic potential barriers on each binary junction. This is constructed in the following manner. First we define a new coordinate system on the triangular simplex defined by i , i = 1, 2, 3

(47)

i·

+

i< j

i< j

(u , v )

1 2

. 1= 2 = 1/2

For n = 3 phases this is identical to the simpler (which can be verified by setting 3 = 1 1 2 in both cases):

3. Advocated by [6]

Wij

4

i< j

2. Used in [30] n

2 i (1

>

3 i

Wij |2 .

2

4 i

i=1

1. Used in [28,12] i 1

i)

1 = 2 = 3 = 1/3

fDW =

In summary, by discounting Eq. (28), there are four gradient energy candidates that satisfy the 2 phase reduction correctly. We summarise these here:

n

n

2 i (1

Wij

i.

j =1

which has the practical feature that higher order junctions have a higher barrier height, but combining this directly with the binary junction height, Wij is not possible (due to the single index sum). So [6] postulate, in common with their treatment of the gradient term,

(44)

giving a form which seems like a generalisation of the form given in [2], but which incorporates double index terms naturally:

fG =

Wij | i || j | i=2

1 = 2 = 3 = 1/3

Toth et al. [6], present a more general (n phase) double well potential,

(43)

We might try to create new formulations that satisfy the constraint Eq. (37) as follows: 2

>

j=1

i

Wnm.

(51)

and with the single index sum

(42)

Similarly, for [30] we have:

=

. 1= 2 = 1/2

n=3 i 1

Wij | i || j |

Wnm (

2 2 i j

j =1

which is easily verified when, say, Wij = 1. This is in contrast to the double obstacle found in [30]

(41)

since

nn

i.

which implies a barrier between the two remaining phases, it does not provide a higher barrier still (even with equal Wij ), at the coexistence of a higher number of phases. For example, for three phase only

Wnm n2 , nm

2W nm n m ,

=

i·

j

single phase

which shows that this form does not reduce correctly. On the other hand, for [28,30] constraint Eq. (37) is satisfied. In [28] 2

j i

(39)

So in this case nn

i=1

1 Wij 2

In [12], there is a problem with the double well potential fDW = Wij i2 j2 . Although it reduces to single phase when i, j k = 0, k

i m

4C 2Wnm.

=

n

4.2. Double well potential

(38)

Then nn=

n

[

1

= uv,

2

= u (1

v ),

3

=1

u].

(57)

This consists of lines parallel to constant 3 (u = constant, 0 < v < 1) and radial lines from 3 = 1 (v = constant, 0 < u < 1) to the opposite side. We then construct a function

i.

(48)

fDW (u, v ) = 16u4 (1 5

u/2) 4v 2 (1

v)2,

(58)

Computational Materials Science 171 (2020) 109085

P.C. Bollada, et al.

Fig. 1. A generalisation of the Folch potential (for 3 phases) over a triangular simplex with bulk phases located at the three vertices. This construction retains zero derivative normal to the boundary, and has a high point at the triple junction, like the Folch potential. It also has variable barrier heights on the binary junctions (in a ratio, here, of 1:2:3. Such a potential provides gradient driving forces that never leave the simplex.

Fig. 2. Shows the potential for multiphase field models without a constraint. The potential is instead modified so that bulk phases are energetically preferred.

which has the property of being zero on two sides of the triangular simplex and quartic on the remaining side, and also zero derivative in the u direction by design on the non-zero edge. We chose the interpolating function (a monotonic function that maps [0, 1] to itself), 16u4 (1 u/2) 4 , over other candidates because it leads to a final polynomial expression (Eq. (60)), which another interpolation may not, e.g. 3u2 2u3 (see later for the effect of using this function). When transforming back to i coordinates using

u=

1

+

2 ,v

1

= 1

+

quantities naturally. However, we have not addressed the choice of gradient terms. It is not so easy to assess the gradient terms because of the difficulty of visualising functions of both i and i. The relations between 2 , W and given in Wheeler [4] is 2

2 2 1 2 (2

2)

1

4,

Wij fij .

a1 2

i

2 i

+

a2 4

4 i

+

a3 2

i
2 2 i j

=

M

x /2

=

1 1+e

x

.

(66)

(67)

The interface, when 2

(x ) =

=

1 , 2

1

exp(x / ) =

is thus located at x = 0 . Then

(1

)

(68)

and so

(62)

1 1 ( )2 = 2 2 2

2 (1

) 2.

(69)

Further differentiation of Eq. (68) gives

1

=

F i

(65)

x

1 . 1 + exp(x / )

=

where a1, a2, a3 are chosen so that there is a set of minima in each bulk phase, e.g. i = 1, e.g. a1 = a2 = 1, a3 = 2. An example of this is shown in Fig. 2 for just two phase fields. There are minima at 1 = 1, 2 = 0 and 1 = 0, 2 = 1. The multiphase field equations take the form i

x

So let us define a flat surface, normal to the x direction. by defining

Other multiphase models that do not impose the constrain i i = 1, include [31–35]. Central to this model design is a double well potential typically of the form

fDW =

ex e ex + e

1 e x /2 (1 + tanhx /2) = x /2 2 e +e

(61)

i
(64)

and consequently

(60)

and cyclically for the other indices. Note that when j + j = 1, Eq. (60) reduces to fij = i2 j2 . We finally set the double well potential of all three contributions using:

fDW =

12

W=

tanh(x )

we obtain the polynomial function

f12 =

,

which we wish to reproduce by having a fresh look at the 1D solution. The tanh function has the identity

(59)

2

=6

(1

2 )

=

1 2

(1

2 )(1

)

(70)

or

(63)

=

where by implication the mobility matrix is diagonal: M = M I . Hence, the constraint is not imposed and indeed is not necessarily met in these models, instead the model construction seeks to make it energetically favourable to be in any one bulk phase.

E=

4.3. Gradient Energy

is

1 2 2

[

2 (1

) 2]

(71)

Noting that the equilibrium equation for the energy functional

E

We have narrowed the choice of double well potential by seeking a construction that is simple, yet still incorporates double index 6

=

1 2

2|

2

2

|2 +

+

1 2

1 2

2 (1

[

2 (1

)2,

) 2] = 0

(72)

(73)

Computational Materials Science 171 (2020) 109085

P.C. Bollada, et al.

Fig. 3. Mapping of a series of circles on the x , y plane (left) to the triangular simplex (right), using Eq. (83). The origin of the x , y plane maps to the centre of the [ 1, 1]: x tanhx . simplex. Points on the simplex boundary are mapped from infinity in the x , y plane. Such a mapping is a 2D generalisation of the mapping R

we see that satisfying Eq. (70) also satisfies Eq. (73) and thus the latter has a solution Eq. (67). Introducing dimensional energy parameters, and W, we can relate physical surface energy, , to these parameters using our 1D equilibrium solution:

1 2

= =

2| 2

1

2

0

+

2 x|

+

W 2

1 W 2

2 (1

(1

where we use the notation this is to set 2

=6

,

W=

2

)2dx=

)d = x

2

1 2 W + 6 2 2

+

2

W 2

2 (1

,

5. A multiphase solution Having considered the manipulation of a 1D solution to produce a single phase model. We now construct an analogous 2D solution for three multiphase fields, i . Our motivation is to use these phase fields to illustrate and analyse the gradient function term in the free energy. For this mapping, 1 (x , y ), 2 (x , y ) and 3 (x , y ) , we can, for example, evaluate gradient fields such as q (x , y ) = i i· i . If, further, the mapping is one-one and invertible, we can associate q for every value of i , i = 1. .3 in the simplex. This allows us to see the relation q ( ¯) q (x ( ¯ )) , i.e. q as a function of 1, 2 and 3 . We now postulate a 2D multiphase solution, which is a generalisation of the tanh profile. Our only guide being: to use something analogous to the 1D case for each field; the solution has a triple junction; behaves like a two phase field away from the triple junction; and that the sum of these fields must be unity. Define three equal length normals to three interfaces such that their sum, n12 + n23 + n31 = 0 (an equilateral triangle, but we relax the equal length constraint later). Using x ij x ·nij , which implies x12 + x23 + x31 = 0 and x ij = xji ,we postulate a set of phase fields:

)2dx (74)

(x ) . The most natural way to satisfy

6

(75)

where we find that we do not reproduce Wheeler’s result, [4], Eq. (64), and so use Eq. (75) from hereon. Using this, multiphase generalises as 2 ij

=6

Wij =

6

ij

=6

Wij =

6

1

(77)

ij

A constant interface width model has 2 ij

(76)

ij ij,

ij

=

(78)

1

ij

.

(79)

= W ij /6,

(80)

where W is a constant with units of energy density. This gives 2 ij

=W

2 ij

Wij = W ,

+

2

+

3

=1

(83) 2,

3.

These fields fortunately (84)

without further manipulation. A visualisation of the mapping Eq. (83) (and cyclically) is shown in Fig. 3. On the left the axes are the cartesian x , y and on the right the 3 simplex so that any value of ¯ is located within. Any point x 0 = (x 0 , y0 ) maps to a unique point ¯0 in the 3 simplex on the right, and the circles on the left map to closed curves on the right. The diagram reveals that the mapping is one-one and invertible, with the origin mapping to the centre of the triangle and infinity mapping to the boundary of the simplex. To get a better feel for the functions, Eq. (83), consider the limiting case where x ·n12 is near zero, and x ·n13 , then

On the other hand, a model that we advocate in Section 6 is where surface energy is proportional to interface width, ij

1 1 + e x12 + e x13

and cyclically for the other phase fields, have the property

so that

ij

=

(81) (82)

all constant and equal. We cannot, and do not, state that physical interface widths, ij are proportional to their respective binary surface energies. We do find, though, in Section 6, that making this assumption allows us to control triple angles in a relatively natural manner. This is by establishing a formulation for equal triple junction angles and then transforming the formulation to arbitrary angles. This not only changes the angles at the junction, but also changes the interface widths, ij .

1

1 1 + e x12

(85)

as in Eq. (66). Thus binary junctions are found, in (x , y ) space, well away from the origin in a direction normal to any of the sides of the triangle defined by the three normals. 7

Computational Materials Science 171 (2020) 109085

P.C. Bollada, et al.

Fig. 6. A proposed gradient energy construction, Eq. (49), has a maximum at the triple junction, but minima near the binary junctions.

Fig. 4. Tiaden-Nestler-Wheeler gradient energy, Eq, (46), has a severe low value at the triple junction. The black triangle at the base represents the 3 simplex 1, 2 , 3 such that i i = 1.

5.1. Visualising and constructing the gradient term Using the multiphase field, Eq. (83), we visualise the gradient term for different models, and also construct a gradient term that matches a potential. Using the constructed solution we can plot each candidate gradient energy as a function of i to inspect their characteristics. We do this in Figs. 4–7, being the gradient terms of the models detailed in Eqs. (46)–(49) and a model to be developed subsequently, respectively. Here we see the energy profile for each model. In particular the model Eq. (46) displays a local minimum at the triple point precisely analogous to the energy profile of the double well potential, i < j i2 j2 . One observation is that 47 and 48 have identical energy profiles (for constant n 2 Wij ) as they should since one can add multiples of | i i | = 0 to any gradient term. We note that these do not have global maximum at the triple junction. The model postulated here for no other reason than simplicity, Eq. (49), has a more attractive energy profile, and finally, a gradient energy model that matches the Folch exactly is shown in Fig. 7. We should emphasise that the equilibrium solution for each choice of gradient term will not match our analytical solution and therefore the actual shape of the gradient energy function as a function of phases, i will not be as shown in Figs. 4–7, which we can only take as indicative. We naturally wish to constrain our free energy constructions so that they reduce correctly on binary junctions, but agreement on the simplex edge is insufficient to guarantee a correct multiphase construction, which we feel must have similar properties to the potential term, e.g. a maximum at the triple junction and zero normal gradients on the

Fig. 7. New proposed gradient energy construction, Eq. (95), has a maximum at the triple junction, correct saddle points at the binary junctions. In fact this is 2 designed to match the Folch-Plapp potential, i i2 (1 i) .

boundary of the simplex. To this end we examine again the Folch-Plapp double well potential 3

fDW = i=1

2 i (1

i)

2 /2.

(86)

This has all the correct properties and has the benefit of being a low order polynomial. Certainly, for say, 3 = 0, 1 = = (1 2) ,

fDW =

2 (1

(87)

) 2.

It is instructive to use our multiphase solution Eq. (83) to plot this potential over x , y space as shown in Fig. 8, where it is worth observing the higher energy over the triple junction. But it is also useful to inspect its properties along the line defined by )/2 . This is given by 3 = , 1 = 2 = (1

1 ( 1 + )2 (9 16

2

+ 2 + 1)

(88)

and shown in Fig. 9. Using this function, and inspecting the relations between phases of Eq. (83) we were able to construct the following gradient function

fG = i
1 ( 2

i

+

j )[4

3(

i

+

j )]

i·

j

(89)

which equals the Folch-Plapp double well potential, Eq. (86), throughout the simplex (for three phases only). At this point we can postulate an alternative surface energy model to [1] that avoids the functions of ¯ outside the summation:

Fig. 5. Folch-Plapp-Toth-Steinbach-Piazola model, Eqs. (47) and (48) has a small local maximum at the triple point), but a pronounced maximum at the binary junction. 8

Computational Materials Science 171 (2020) 109085

P.C. Bollada, et al. n

fDW =

Wij

3

2 (

i
i

2

i

+

j)

j

2 2 i j

2

(92)

and so the more general formulation is n

1 2

f= i
2 ij ( i

+

n

+

Wij

j )[4

3(

3

2 j

2 i

+

j )]

i·

j

2 2 i j

( i + j )2

i
i

(93)

where, for constant interface width, we have n

1

f= i
2 0( i

Wij { 2

+

+

3

2 i

j )[4 2 j

3(

i

+

j )]

i·

j

2 2 i j}

( i + j )2

(94)

Despite this model reducing correctly across each binary junction, we find, computationally, that it does not reproduce the theoretical angles at the triple junction. That said we postulate that such a model may well benefit from the addition of a hump potential such as advocated in [1] (see Appendix B for an indication of the properties of such a function). We now focus the paper on a more rigorously derived formulation that necessarily produces the correct triple junction angles. We seek to generalise the Folch-Plapp formulation, [2], which we write in the form

Fig. 8. This is a plot of the Folch-Plapp potential over the multiphase solution Eq. (83).

3 2

FS =

i· i

i

+

2 i (1

i)

2.

(95)

which has a 1D equilibrium solution with only two phases present: 1

=1

2=

1 , 1 + exp(x / )

3

= 0.

(96)

6. Energy formulation to control triple junction angle We seek a generalisation with a stronger mathematical underpinning: one that mathematically guarantees that the theoretically correct triple junction angles are recovered. Before doing this we wish to illustrate a key transformation in a simpler 1D setting. Consider two 1D solutions, 0 and , that differ only in their interface width: 0

=

1 1 + exp(x )

(97)

and Fig. 9. View of the Folch-Plapp double well potential from the bulk of one phase to the middle of the binary junction of the other two phases. This function has zero gradients at the extremes and a maximum at = 1/3. 3

f= i
1 2

2 ij ( i

=

j )[4

3(

i

+

j )]

i·

j

+ i
i

j

)4

d

(90)

2u3) v 2 (1

v )2

(99)

0

=

0

x

x X d . X

(100)

At a point x = X = 0

At this stage, we only state this formulation as being a plausible alternative to that found in [1]. On inspecting the double well potential for higher values of n = 3 we find that n = 4 gives a higher potential for i = 1/4, i = 1, 4 , but this feature does not extend to higher values of n > 4 . This issue is rectified by choosing a lower order interpolating function for f (u, v ) in Eq. (58), i.e.

fDW (u, v ) = (3u2

1 , 1 + exp(X )

we find by the chain rule that

3

Wij i2 j2 (2

(98)

By writing X = x/ so that

=

+

1 . 1 + exp(x / )

0

x

= 1/ 0

X 0

(101)

and so at x = X = 0

d

(91)

0

=

x d = d . X

(102)

Following the above as a template for multiphase, we use the solution Eq. (83) as a reference solution, and introduce the notation i0 to distinguish it from the general solution, i . Thus, in a slightly different

This gives the double well potential 9

Computational Materials Science 171 (2020) 109085

P.C. Bollada, et al.

from Eq. (83):

1 0 1

Now the term in the bracket can be written

1

= 1 + exp(x12) + exp(x13),

= 1 + exp(x23) + exp(x21),

0 2

1

0 1

= 1 + exp(x + exp(y

y ) + exp(x x ),

1 0 3

1

z ),

0 2

= 1 + exp(z

= 1 + exp(y x ) + exp(z

0 1

= 1 + exp(x

y ) + exp(x

+ 1/exp(x

y )exp(y

z ),

1 0 2

= 1 + exp(y

(104)

exp(x

y) =

z)

d d d

xij = ln( j0 )

ln(

1/3 2/3 1/3

= 1 + exp(x /

0 i ).

(107)

dXi =

0 1

ln(

1/3 1/3 2/3

ln( ln(

0 1) 0 2) 0 3)

0 2

1

2/3 1/3 1/3

X Y = Z

Y ) + exp(X

Y

Xij

=

i (x ,

1

z / 3)

=

i (X ,

W=

Z)

ln( 2) . ln( 3)

(111)

2 ij

(112)

2 i

j

xa

=

XA

(121)

i

(122)

0 j

0T

0

T DDT

=

(123)

W 4

[

2 i

i·

i

+

i

2 i (1

i)

2]

(124)

6

12

(125)

=

6

ij

(126)

W

=

2 i

+

2 j

i

to the double indexed,

ij .

(127)

XA j

j

xa

.

=

1 ( 2

2 ij

+

2 ik

2 jk ),

i

j

k.

(128)

We need the above relation because the binary surface energies, ij relate directly to the binary interface widths by Eq. (80). Eq. (128) appears to present us with a problem in that it is possible that i2 0 . The critical angle, when this occurs is when one of the angles at the triple junction is /2. However, computational results demonstrate that this is not a problem, where we find, for example, that even 2 1/3, 22 = 32 = 1, can be simulated successfully. 1 =

(114)

Y , Z)

j

3d

which gives on inversion

ln( i )

(113)

xb

(120)

Finally, we need to relate the single index, This is readily done since by definition

ln( 1)

1/3 1/3 2/3

y , z ),

0 i b x

(119)

and

(110)

ij

x ·Nij = ln( j )

0 i

dx dy dz

12

where Xi = x i / i . That is, they can be represented by the same function but on a different space. Using this examine 0 i = a x

1 1 2

where W can be taken from any of the barrier heights, e.g.

(109)

then the general multiphase field is given, using the same function, by i

1 2 1

1 dxi 3

=

fS =

where Nij represents a triple junction where angles are not necessarily equal. If we define the reference function 0 i

(118)

Thus we postulate the surface energy as being

Using this we find that

X

0 i

0 i ·

etc, we have

1/3 2/3 1/3

dz)

This means that (108)

1

where X = x/

z)(dx

= 1/3, x = y = z = 0 we have

i

2 1 1

1 9

dy) + exp(x

Consequently from Eq. (116) we have, at the origin:

and cyclically. By writing Eq. (109) in the equivalent form

= 1 + exp(X

y )(dx

R = D.

.

1

1

(117)

and so the inverse relation

y / 2) + exp(x /

1

= exp(x

=

d

Eq. (104) and its inverse, Eq. (108) relate to an equal angle triple junction. We now consider, the more arbitrary solution where all three angles can vary by tilting the equal angle solution as follows:

1

0 1 0 2 0 3

(106)

=1 By using the substitution and fixing a plane through (x + y + z = 0) , we find the inverse transformation Eq. (104): 2/3 1/3 1/3

0

If we further use the knowledge that dx + dy + dz = 0 we find that

0 3

x y = z

0 1 2 1

0 2 0 1

and in general

x ·nij

(116)

0 0 . 3.

2

At the midpoint,

and so we find that

exp(x12 )

0 0

d

(105)

y)

Jib DbA (J 1)Aj.

From Eq. (104) we have

which can be written in terms of just two coupled variables on the righthand side, (x y ) and (y z ) ,

1

0

1

D

z) y)

xb XA XA j

Since the left hand side and right hand side of the final term are inverse matrices of each other and the middle matrix D is a diagonal matrix

and consider that this is defined in a plane that passes through the 3D e1 [1, 0, 0], e2 [0, 1, 0], e3 [0, 0, 1], points with normals, nij = ej ei . In this coordinate system we can write this as

1

0 i b x

=

j

(103)

= 1 + exp(x31) + exp(x32)

0 i

Rij

0 3

6.1. Some failed models It may be useful to the reader to note that we found other variants of the above model to have only limited applicability in simulation. We

(115) 10

Computational Materials Science 171 (2020) 109085

P.C. Bollada, et al.

find that a formulation of the form

We finish off with a restatement of the proposed model for the surface energy density

3 2 i Ui

fS =

i·

2 i (1

+

j

i=1

i

i)

2,

(129)

fS =

W 4

where Ui = i + 3 j k , or Ui = i develops an instability as the triple junction angle approaches /2 (from larger angles). Since there are reports on smaller angles than /2 such as the groove angle at solidsolid-liquid triple junction [36,37], there remains a limitation with this model. We also add, in passing, that the model 2 ij

fS =

i·

j

+

i
i

2 i (1

i)

2,

i=1

i)

=

2 i

1 ( 2 1 = ( 2 =

(131)

+

2 31

+

2 2 12), 2

2 12

+

2 23

+

2 31).

=

1 ( 2

2 31

+

2 12

+

=

2 1

2 2 23), 3

(133) (134)

2 4,

+

we determine the remaining single index coefficient 2 2 4 = 14

2 1

=

2 14

1 ( 2

2 23

+

2 31

+

2 12 ).

2 jk |i j k

(136)

1 , 1 + expx12 + expx13 1 = , 1 + expx31 + expx32

2

=

1 , 1 + expx23 + expx21

3

(137)

D

fS dxdy

(138)

to give us the total energy associated with these fields for any chosen set of normals, nij (see Fig. 10). By choosing the domain, D to be a large enough circular disc about the origin, and varying one or more of the normals the total energy over D depends on the normals, and thus the angles at the triple junction. What we seek in doing this is the set of angles that gives an energy minimum. This then can be compared with the theoretical set of angles associated with the choice of i2 . We limit this test to that in which two angles are equal so that the energies will be a function of the remaining angle, , only, i.e. FS ( ) . In Fig. 11 the sub Fig. 11a–f each give a range of FS ( ) for a different choice 1 = 2 = 1 and 3 = 0.7, 0.8, 0.9, 1.0, 2.0, 4.0 , which correspond to the theoretical angles: 109°, 113°, 117°, 120°, 143° and 160° respectively. For example, 1 = 1, 2 = 1, 3 = 2 corresponds to a triangle with vertices, in 3D, [1, 0, 0], [0, 1, 0], [0, 0, 2], with length sides 2 , 5 , 5 , with smallest angle 2sin 1 (1/ 10 ) = 37 degrees and therefore the largest angle in the triple junction being 180 37 = 143 degrees (used in Fig. 11e). For a wide range of cases the theoretical angle is well approximated by the minimum of FS ( ) .

and then again using Eq. (128), 2 14

2 ij

=

=

FS =

(132)

2 23

12

6 ij W

where x ij = x ·nij and n12 + n23 + n31 = 0 . By setting W = 1, and prescribing i2 for each phase we can evaluate fS in Eq. (136) at any point. Furthermore, we can integrate fS over any domain, D,

For n = 4 , for example, the distances from the origin in 4D, 1, 2, 3, 4 , define the vertices of a tetrahedron and so we determine one of the sides by (using the definition Eq. (128)): 2 1

=

6 12

We test directly the model Eq. (136) at a constructed triple junction. We prescribe a triple junction field at the origin in 2D by

2

2 j.

+

2]

7.1. Angle dependent energy minimum

where the binary interface energies are (still) related by 2 ij

i)

7. Tests and results

(130)

1

2 i (1

+

2 i (1

+

j i

n i

i

i

2 i

A more general n > 3 has not been explored in this paper, but clearly a good starting point is the model i·

i·

ij

6.2. Generalisation to n > 3

2 i

2 i

W=

was stable in simulation but fails to deliver the correct angle at the triple junction despite affecting the interface width in a similar manner to that of Eq. (124) (although this model might be rescued by the addition of a triple junction barrier as in [1]).

fS =

[

(135)

We speculate that only n binary interface energies are independent for n > 2. 6.3. Summary An overview of what we have done: 1. By recognising that a more general solution with arbitrary angles at the triple junction can be produced by an invertible mapping, we rewrite the equal angle formulation, which is specified in terms of the equal angle solution i0 , in terms of the general angle solution i . This results in a general multiphase model which controls angles at triple junctions via the binary surface energies, ij . 2. The mapping makes convenient use of a two dimensional surface in 3D, which gives a more direct relation between the three phases, i0 and the three Cartesian coordinates, x , y , z (than into the 2D coordinates x , y ). Moreover, the generalisation to a plane which cuts the coordinate axes at arbitrary points, governed by the interface widths, follows comfortably. 3. The formulation appear computationally stable for all angles at the triple junction apart from the extremes, 0° and 180°.

Fig. 10. Circular disc containing three phases, 1, 2 and 3 , at the triple junction. The angle, is made to vary by keeping |n31| = |n23| fixed and varying n12 . By integrating over the disc we obtain an energy, FS ( ) as a function of . 11

Computational Materials Science 171 (2020) 109085

P.C. Bollada, et al.

(a) Theoretical angle 109◦ ; computed angle 106 ± 1◦ ; δ3 = δ2 = 1, δ1 = 0.7

(b) Theoretical angle 113◦ ; computed angle 113 ± 1◦ ; δ3 = δ2 = 1, δ1 = 0.8

(c) Theoretical angle 117◦ ; computed angle 117 ± 1◦ ; δ3 = δ2 = 1, δ1 = 0.9

(d) Theoretical angle 120◦ ; computed angle 120 ± 1◦ ; δ3 = δ2 = 1, δ1 = 1.0

(e) Theoretical angle 143◦ ; computed angle 143 ± 1◦ ; δ3 = δ2 = 1, δ1 = 2.0

(f) Theoretical angle 160◦ ; computed angle 160 ± 1◦ ; δ3 = δ2 = 1, δ1 = 4.0

Fig. 11. (a)–(f) Show the energy as the angle sweeps around the theoretical angle.

Clearly, almost any formulation that treats all phases equally will return an optimum angle as = 2 /3 = 120°. However, one formulation [28] actually had a maximum energy at this angle. We have only given our model’s results for these sweeps as none of the other models came close when the angle varied from 120°.

7.2. Phase field implementation We now examine how our postulated model behaves in phase field under dynamic conditions close to equilibrium, where we expect there to be some flattening or sharpening of the eutectic interface curves due 12

Computational Materials Science 171 (2020) 109085

P.C. Bollada, et al.

Parameter

Value

M =

0.01 0.4 Multiples of 3 0.05 32 × 16 x = 0.125 1 nm 0.67 ns 3 (x-direction) 2 × 16 (y-direction) c1 = 1/2, c2 = 1/4, c3 = 3/4

3

2

1

T Domain size Mesh size Length unit Time unit Initial eutectic length Initial lamellae width Common tangents:

(139)

fB = g1 f1 + g2 f2 + g3 f3

Table 2 Table of values for the phase field simulation.

where gi = g ( i )/ functions, and

g ( i ) with g ( ) = 3

2

2 3 , are interpolation

f1 = [c

c1 + (c1

c2) g2 + (c1

c3) g3]2 + g1 T ,f2

= [c

c2 + (c2

c3) g3 + (c2

c1) g1]2 + g1 T ,f3

= [c

c3 + (c3

c1) g1 + (c3

c2) g2]2 + g1 T .

(140)

The concentrations ci are from common tangent construction and are given in Table 2 along with other parameters for the simulation. Note that when g1 = 1, we have f1 = (c c1) 2 + T ; when g2 = 1 we have f2 = (c c2) 2 ; and when g3 = 1 we have f3 = (c c3) 2 . All three of which are the bulk free energies of the pure three phases. Thus the combination in Eq. (139) is a phase dependent interpolation. We adopt this interpolation technique over the WBM method to minimise spurious phase generation between any two phases (see [7,8] for further discussion). Fig. 12(a)–(d) show solute profile contours from the multiphase field simulation. Superimposed on the plot is the theoretically correct equilibrium angle corresponding the input, i . These results show the effect of differing binary surface energy terms on the manner of growth in this dynamic setting with a clear flattening off of the curved surface

to our method of implementation of the binary surface energies, ij . The phase field implementation for these tests uses Model Eq. (136) with mobility model given by Eq. (17), and a bulk force to be described. The model used is a simplest example model where each of the three free bulk energies are represented by quadratics and the diffusion throughout the domain is constant. We also adopt the following bulk free energy:

Fig. 12. (a)–(d) Show simulation results using the new model Eq. (124). The contours of the solute profile and the theoretical angles (60°, 90°, 120° and 60°, respectively), are shown superimposed for direct comparison. The values for 12 are 1/3, 0, 1, 16 respectively. A distinctive feature of the model is the relation between angle and interface width, so that, for example,the liquid solid interface increases with angle. 13

Computational Materials Science 171 (2020) 109085

P.C. Bollada, et al.

Fig. 13. (a) and (b) show simulation results for 143°, and also with Ui =

i

+3

j k

(rather than Ui = 1 as in Eq. (124)) and with increasingly finer interface widths.

the correct angle on the eutectic simulation results in Fig. 12 clearly show the model works well even in this dynamic setting. Further narrowing of the interface width has the effect shown in Fig. 13, which suggests that in the limit of near equilibrium and zero interface width the theoretical angle is approached at the triple junction as the interfaces become sharp. 8. Conclusion Commencing with a review we have narrowed the field of possible multiphase models to a handful. We have developed some new tools that not only narrow the field of candidates, but create new ones. In seeking a consistent model that allows binary interface energies to be incorporated into the model in a way that correctly reduces to a single phase formulation, we postulated the model Eq. (90). This incorporated a potential well that allows barrier heights proportional to the binary surface energies, which in turn facilitated the future development of a model that keeps interface widths equal. However, unlike the model found in [1], computation with this model does not return correct theoretical angles. This illustrates that correct reduction across any binary junction is not sufficient to reproduce triple junction angles correctly, and so modification by the addition of triple junction terms becomes mandatory. In Section 6 we postulate a more rigorous theoretically justified multiphase field model for alloy solidification, Eq. (136). This incorporates binary surface energies, ij , in a manner that mathematically guarantees the correct reproduction of triple junction angles between phases. In deriving this model we have made use of an analytical multiphase solution, Eq. (104). The analytical solution is modified by changing the underlying arbitrary triangle, which suggests that an affine transformation of an equal angle formulation will deliver the sought after general formulation. The only potential down side of this is that the relative interface widths, ij vary with the binary interface energies, ij . Implementing the model Eq. (136) together with a simple bulk energy model we find that the effect of the binary surface on the eutectic growth is largely in line with the angles expected at equilibrium. This observation is further backed up by integrating the energy field as a function of underlying fields with varying angles, and finding that the correct angle implied by the binary surface energies is at a minimum at that angle - see Fig. 11.

Fig. 14. The evolution of the angle at the triple junction in time in units of 0.67 ns. The initial condition sets the correct angle but the symmetric boundary condition and growth changes the profile.

approaching the tangent indicated by the equilibrium correct angle. These eutectics grow from the left, and though they have evolved towards a steady velocity state they are not in an equilibrium state. We choose a small value of T to give near equilibrium behaviour, but the boundary conditions being symmetric guarantee that the liquid-solid interfaces are curved. Because of the curvature at the interface it is difficult to prescribe a quantitative measure of angle. One such measure is given by 1

(141)

= cos 1 (n12· n13)

(and cyclically) where

nij =

|

(

j

i

i

j) H (

¯ )dxdy

(

j

i

i

j) H (

¯)dxdy|

(142)

and H ( ¯ ) is some function that limits the integration area to the vicinity of the triple junction. If H is too local then the angle will be too small and, conversely, if too large (e.g. H = 1) the angle will be too large. Using H = 1 2 3 we found we could reproduce the angles between 100 and 143 within a couple of degrees, see for example Fig. 14, but were less successful outside this region. That said the superposition of

14

Computational Materials Science 171 (2020) 109085

P.C. Bollada, et al.

CRediT authorship contribution statement

Declaration of Competing Interest

P.C. Bollada: Formal analysis, Conceptualization, Validation, Software, Writing - original draft, Investigation, Methodology, Visualization. P.K. Jimack: Writing - review & editing. A.M. Mullis: Supervision, Funding acquisition, Writing - review & editing, Project administration.

The authors declare that they have no conflict of interest. Acknowledgements This research was funded by EPSRC Innovative Manufacturing Research Hub in Liquid Metal Engineering (LiME), Grant No. EP/ N007638/1.

Appendix A. Plotting gradient functions over the simplex We present this section to give a clear exposition of how to convert scalar gradient functions of phase to functions of phase only using Eq. (104). Beginning with 104, but dropping the superscript 0, we find that

d

1 =exp(x 2 1

y )(dx

dy) + exp(x

z)(dx

dz)=

2

(dx

dy) +

1

3

(dx

dz)

(A.1)

1

so that

d

1

=

2 1 (dx

dy) +

3 1 (dx

(A.2)

dz).

We then make the cartesian identification dx = i , dy = j, dz = k and d 1

=

2 1 (i

j) +

3 1 (i

1

=

1

to write (A.3)

k).

and similarly 2= 3 2 ( j

k) +

1 2 (j

i)

3

=

1 3 (k

i) +

2 3 (k

(A.4)

j).

Hence, for example, 1·

1

=2

2 2 1 2

+2

2 1 2 3

+2

2 2 1 3

(A.5)

and consequently

1 2

3 i·

i =2

i=1

i
2 2 i j

+

2 1 2 3

+

2 2 3 1

+

2 3 1 2 , =2

i
2 2 i j

+

1 2 3( 1

+

2

+

3)

=2 i
2 2 i j

+

1 2 3

(A.6)

Appendix B. Gradient energy equivalent of the hump function Somewhat tangential to the main argument of our paper, but of interest, is an observation regarding a “hump” function, which is similar to that used in [1] to control angles ([1] use its square root). This is given by

fH = (

1 2 3)

(B.1)

2,

which has the property of being zero and has zero gradient on the simplex boundary. A plot of how this energy varies as a function of one normal varying is given in Fig. B.15 together with an equivalent cross product gradient term (that agrees at 120 degrees)

Fig. B.15. Energy as a function of one angle around a triple junction for both the hump function (green) and cross term (cyan). By construction the two energies agree at the 120 deg angle. 15

Computational Materials Science 171 (2020) 109085

P.C. Bollada, et al.

Fig. B.16. Energy as a function of one angle around a triple junction for both the double well potential (red) and the equivalent gradient term (blue). By construction the gradient energy agrees with the double well term at the 120 deg angle.

fC =

4 9

|

i

×

2 j| .

(B.2)

i< j

For comparison we show the Folch-Plapp potential with the equivalent gradient term in Fig. B.16. The suggestion here is that just as the potential and gradient term balance, so should the hump function balance with the cross term, Eq. (B.2). However, inclusion of the cross term means having an inconvenient forth order surface energy term.

[19] K. Glasner, Solute trapping and the non-equilibrium phase diagram for solidification of binary alloys, Physica D 151 (2001) 253–270. [20] B. Echebarria, Quantitative phase-field model of alloy solidification, Phys. Rev. E 70 (2004) 061604. [21] N. Opoku, A quantitative multi-phase field model of polycrystalline alloy solidification, Acta Mater. 58 (6) (2010) 2155–2164. [22] Seong Gyoon Kim, Won Tae Kim, Toshio Suzuki, Machiko Ode, J. Crystal Growth 261 (2004) 135–158. [23] J. Eiken, Multiphase-field approach for multicomponent alloys with extrapolation scheme for numerical application, Phys. Rev. E 73 (2006) 066122. [24] A. Choudhury, B. Nestler, Grand-potential formulation for multicomponent phase transformations combined with thin-interface asymptotics of the double-obstacle potential, Phys. Rev. E 85 (2012) 21602. [25] I. Steinbach, L. Zhang, M. Plapp, Phase-field model with finite interface dissipation, Acta Mater. (2012) 2689–2701. [26] A. Choudhury, M. Kellner, B. Nestler, A method for coupling the phase-field model on a grand-potential formalism to thermodynamic databases, Curr. Opin. Solid State Mater. Sci. 19 (2015) 287–300. [27] P. Bollada, P.K. Jimack, A.M. Mullis, A numerical approach to compensate for phase field interface effects in alloy solidification, Comput. Mater. Sci. 151 (2018) 338–350. [28] J. Tiaden, B. Nestler, H.J. Diepers, I. Steinbach, The muiltiphase-field model with an integrated concept for modelling diffusion, Physica D 115 (1998) 73–86. [29] A. Basak, V.I. Levitas, Nanoscale multiphase phase field approach for stress - and temperature – induced martensitic phase transformations with interfacial stresses at finite points, J. Mech. Phys. Solids 113 (2018) 162–196. [30] I. Steinbach, F. Pezzolla, A generalized field method for multiphase transformations using interface fields, Physica D 134 (1999) 385–393. [31] L.Q. Chen, W. Yang, Phys. Rev. B 50 (21) (1994) 15752. [32] D. Fan, L.-Q. Chen, Acta Mater. 45 (1996) 611–622. [33] A. Kazaryan, Y. Wang, S. Dregia, B.R. Patton, Grain growth in systems with anisotropic boundary mobility: analytical model and computer simulation, Physica Rev. B 63 (2000) 184102. [34] N. Moelans, B. Blanpain, P. Wollants, Phys. Rev. B 78 (2008) 024113. [35] Ofori-Opoku, N. Provatas, A quantitative multi-phase field model of polycrystalline alloy solidification, Acta Mater. 58 (2010) 2155–2164. [36] M. Gunduz, J.D. Hunt, Acta Metall. 33 (9) (1985) 1651–1672. [37] M.E. Glicksman, K. Ankit, J. Mater. Sci. 53 (2018) 10955–10978.

References [1] G. Toth, Phase-field modeling of isothermal quasi-incompressible multicomponent liquids, Physica Rev. E 94 (2016) 033114. [2] R. Folch, M. Plapp, Towards a quantitative phase-field model of two-phase solidification, Phys. Rev. E 68 (2003) 010602. [3] R. Folch, M. Plapp, Quantitative phase-field modeling of two-phase growth, Phys. Rev. E 72 (2005) 011602. [4] A. Wheeler, W. Boettinger, G. McFadden, Phase-field model for isothermal phase transitions in binary alloys, Phys. Rev. A 45 (10) (1992) 7424–7440. [5] P.C. Bollada, P.K. Jimack, A.M. Mullis, Bracket formalism applied to phase field models of alloy solidification, Comput. Mater. Sci. 126 (2017) 426–437. [6] G.I. Toth, T. Pusztai, L. Granasy, Consistent multiphase-field theory for interface driven multidomain dynamics, Phys. Rev. B 92 (2015) 184105. [7] M. Plapp, Unified derivation of phase-field models for alloy solidification from a grand-potential functional, Phys. Rev. E 84 (2011) 31601. [8] P.C. Bollada, P.K. Jimack, A.M. Mullis, Free energy vs. grand potential energy formulations in phase field modelling of alloy solidification, Proceedings of the 7th International Conference on Solidification & Gravity, 2018, pp. 47–51. [9] A. Karma, W. Rappel, Quantitative phase-field modeling of dendritic growth in two and three dimensions, Phys. Rev. E 57 (4) (1998) 4323–4349. [10] A.N. Beris, B.J. Edwards, Thermodynamics of Flowing Systems, Oxford University Press, New York, 1994. [11] O. Penrose, Thermodynamically consistent models of phase-field type for the kinetics of phase transitions, Physica D 43 (1990) 44–62. [12] B. Nestler, A.A. Wheeler, A multi-phase-field model of eutectic and peritectic alloys: numerical simulation of growth structures, Physica D 138 (2000) 114–133. [13] P.C. Bollada, P.K. Jimack, A.M. Mullis, A new approach to multi-phase formulation for the solidification of alloys, Physica D 241 (2012) 816–829. [14] I. Steinbach, A phase field concept for multiphase systems, Physica D 94 (1996) 135–147. [15] J. Groebner, H.L. Lukas, F. Aldinger, CALPHAD 20 (1996) 2247–2254. [16] A. Wheeler, W. Boettinger, G. McFadden, Phase-field model of solute trapping during solidification, Phys. Rev. E 47 (3) (1993) 1893–1909. [17] S.G. Kim, W.T. Kim, T. Suzuki, Phase-field model for binary alloys, Phys. Rev. E 60 (6) (1999) 7186–7197. [18] A. Karma, Phase-field formulation for quantitative modeling of alloy solidification, Phys. Rev. Lett. 87–11 (2001) 115701.

16