Efficient and quantitative phase field simulations of polycrystalline solidification using a vector order parameter
Tatu Pinomaa, Nana Ofori-Opoku, Anssi Laukkanen, Nikolas Provatas
EEfficient and quantitative phase field simulations of polycrystalline solidification using a vectororder parameter
Tatu Pinomaa , Nana Ofori-Opoku , Anssi Laukkanen , and Nikolas Provatas ICME group , VTT Technical Research Centre of Finland Ltd , Finland Computational Techniques Branch, Canadian Nuclear Laboratories, Chalk River, Canada Department of Physics and Centre for the Physics of Materials , McGill University , Montreal , Canada (Dated: February 23, 2021)A vector order parameter phase field model derived from a grand potential functional is presented as a newapproach for modeling polycrystalline solidification of alloys. In this approach, the grand potential density isdesigned to contain a discrete set of finite wells, a feature that naturally allows for the growth and controlledinteraction of multiple grains using a single vector field. We verify that dendritic solidification in binary alloysfollows the well-established quantitative behavior in the thin interface limit. In addition, it is shown that grainboundary energy and solute back-diffusion are quantitatively consistent with earlier theoretical work, with grainboundary energy being controlled through a simple solid-solid interaction parameter. Moreover, when consid-ering polycrystalline aggregates and their coarsening, we show that the kinetics follow the expected parabolicgrowth law. Finally, we demonstrate how this new vector order parameter model can be used to describe nucle-ation in polycrystalline systems via thermal fluctuations of the vector order parameter, a feature which has beenthus far lacking from multi-phase or multi-order parameter based phase field models. The presented vector orderparameter model serves as a practical and efficient computational tool for simulating polycrystalline materials.We also discuss the extension of the order parameter to higher dimensions as a simple method for modelingmultiple solid phases.
I. INTRODUCTION
It is well established that the properties of materials andtheir performance are ultimately determined by their mi-crostructure. One of the main hallmarks of “microstructure”in materials such as metal alloys is their polycrystalline na-ture. One of the quintessential goals of materials processingis therefore to understand, predict, and control the polycrys-talline structure during the solidification stage of their materi-als processing. In the last few decades, this task has been moreand more undertaken using computational modeling based onso-called multi-phase field or multi-order parameter models[1–3]. The approaches use a discrete set of fields to modeldifferent orientations of solids in the liquid, each driven byits own equation of motion that is coupled to diffusion of so-lutes and/or heat, and interacting with other grains in somephenomenological manner designed into the equations. Whilesuch approaches have been very beneficial in pushing the fron-tiers of computational modeling, they suffer from some limi-tations in theoretical consistency such as the incorporation ofnucleation, thermal fluctuations, orientational invariance andrequiring the use of many number of fields to simulate poly-crystalline materials.For accurate assessment and prediction of polycrystallinesolidification, including nucleation, growth and grain merg-ing, it is necessary for a simulation approach to have a self-consistent description of the magnitude of the growth kinet-ics and crystal anisotropy (capillary and kinetic), as well as aproper description of solute segregation and diffusion. Uponmerger of grains, when a cohesive solid network has formed,coarsening ensues, which is dominated by grain boundary mi-gration determined by grain boundary energetics and mobil-ity. Various modeling methods have been developed to studysolidification, predict polycrystalline structure formation, andthe subsequent coarsening. Cellular automaton and kinetic Monte Carlo methods aretwo popular techniques that enable large-scale simulationsof polycrystalline solidification. However, user interventionis required to properly define the growth kinetics, and thesemethods are prone to numerical issues and they have consid-erable limitations in addressing the simulation of alloys anddiffusion controlled kinetics of solidification.The phase field (PF) methodology is typically the mostwidely used state-of-the-art numerical scheme to describe mi-crostructural solidification. Various PF approaches exist forpolycrystalline simulations. One family of models, pioneeredby Kobayashi [4], Granasy [5], Warren [6], and co-workers, ischaracterized by the use both orientation ( θ ) and order param-eter ( φ ) fields (coined θ − φ models) to model solidificationof any number of crystal orientations. This approach allowsfor an efficient representation of a large number of grain ori-entations. It successfully describes basic features of polycrys-talline solidification: controllable grain boundary energy ver-sus misorientation, polycrystalline growth, grain merger, andcoarsening. However, these models are not without their chal-lenges to implement numerically. Moreover, matched thin in-terface asymptotics has not been performed for this family ofmodels, and it is likely to be challenging due to the inclusionof singular terms in the model description, e.g. |∇ θ | gradientterm in the free energy. Additionally, these models have beenshown to exhibit unphysical topological defects, that have re-cently been, at least partially, addressed in 2D simulations [7].Finally, the extension to 3D of this class of models is not im-mediately straightforward, although some schemes have beenexplored [8].Another family of PF models is called multi-phase field(MPF) models [9]. These deviate from traditional PF mod-els in that the phase field φ represents the volume fraction ofa phase, which implies that even the liquid is assigned a phasefield. As with standard PF models, these models compute the a r X i v : . [ c ond - m a t . m t r l - s c i ] F e b energy of an interface, and anisotropy, by tuning appropri-ate gradient energy terms [10]. In addition, interaction en-ergy between solid grains includes various polynomial inter-action terms between phase fractions. However, MPF modelscan suffer from unphysical non-local adsorption effects dueto a complicated free energy landscape wherein the sum ofthe phase fraction fields must be maintained by a Lagrangemultiplier. These need to be eliminated in a somewhat ad-hocmanner by adding higher order polynomial penalties betweenphase fraction fields. Also, it is not clear how to integratethermal fluctuations into the distinct phase fraction fields self-consistently. Another challenge with using MPF models isthat, in principle, they nominally require as many fields asthere are crystal grains, thus requiring significant computingresources to simulate realistic materials domains. To circum-vent this problem requires dynamic grain re-assignment algo-rithms to be implemented [11].Yet another family of PF models is called multi-order pa-rameter PF models. As the name suggests, in this class ofmodels, a phase field φ maintains the traditional interpreta-tion as a physical order parameter that distinguishes betweenordered and disordered states of a crystal. As with MPF mod-els, one order parameter is required for each crystal grain. Anearlier approach developed by Ofori-Opoku et al. [3] useda quadratic repulsion energy term between order parameters,which is all that is required since this class of models doesnot suffer from non-local adsorption effects. This approach ofRef. [3] allowed for a quantitative accounting of solid-liquidinterface kinetics and qualitative description of solid-solid in-terfaces in metals. A weakness with this implementation isthat the grain boundary energy cannot be easily tuned, re-sulting in grain boundary energies corresponding to high an-gle/energy grain boundaries. As with MPF, for multi-order pa-rameter models the Langevin noise based nucleation is prob-lematic to implement when multiple order parameters are fluc-tuating in the same material point due to the use of polynomial(here quadratic) repulsive interaction terms.In this work, we present an alternative phase field modelingscheme to describe solidification and coarsening of multiplecrystallographic orientations in the context of a physical or-der parameter description. We accomplish this by introduc-ing a two dimensional vector order parameter as in Ref. [12],which is then integrated with a previous quantitative modelingscheme for solidification based on a grand potential energyfunctional [13]. This main innovations of this approach arethe trivial generation of a Landau free energy landscape cor-responding to any number of discrete orientations in the con-text of a single physical order parameter field, simple controlof grain boundary energy, and self-consistent incorporation ofthermal fluctuations into dynamical simulations.The paper is organized as follows. In Section II, we in-troduce the vector order parameter model and review its pa-rameterization. Section III examines the model’s capacity toreproduce some important benchmarks quantitatively. Theseinclude, free solidification of a single dendrite, grain bound-ary energy, grain merger and back-diffusion, grain growthand coarsening, grain boundary segregation in polycrystallinegrowth, and noise-induced nucleation. Section IV further dis- cusses some of the qualities of our model and potential furtherextensions. We conclude in Section V. II. METHODOLOGY
This section develops the formalism of the vector phasefield model. In this work, the focus is to demonstrate theproperties of this phase field methodology, and therefore welimit ourselves to binary alloys. We will formulate the modelformally in the grand potential ensemble, meaning that thedependent fields are the order parameter and the chemical po-tential. We begin by describing the grand potential energyfunctional, highlighting its vectorized Landau landscape. Fol-lowing this, we discuss the dynamics of the vector order pa-rameter and the chemical potential. We then highlight the pa-rameter relations necessary to achieve quantitative correspon-dence with the sharp interface model of solidification.
A. Grand potential functional
We start by defining a vector order parameter defined as (cid:126)φ ( (cid:126)x, t ) = ( φ X ( (cid:126)x, t ) , φ Y ( (cid:126)x, t )) . We also project (cid:126)φ into polarcoordinates R and θ according to φ X = R cos( θ ) (1) φ Y = R sin( θ ) , (2)where R can be interpreted as the “traditional” solid-liquidorder parameter, and θ is an additional degree of freedom thatis used to design a Landau functional whose bulk part containsa discrete number of solid wells. This is given by f Landau ( R, θ ) = R − b cos ( N wells θ )1 + b R + R , (3)where N wells sets the number of discrete solid wells corre-sponding to unique crystalline orientations, and b controls thefree energy barrier between neighboring solid wells. Eq. (3)is usually referred to as the “double well function” in con-ventional PF models, where there is one solid well and oneliquid well. An example of the Landau landscape is shownin Fig. 1 with N wells = 8 and b = 0 . . Our approach isanalogous to Morin et al. [12], who used a different orderpolynomial in R to describe the Landau free energy function(polynomial orders R n , n = 2 , , ). Here, we modify theirapproach to a polynomial which is a natural extension of thesingle order parameter PF models typically used for solidifi-cation, thus adopting polynomial orders of R n , n = 2 , , .It is noted that for 2D cubic lattices, the physical crystallineorientation is one fourth of the Landau angle θ . This is illus-trated in Fig. 1. Thus, for example, θ = π (180 o ) correspondsto the maximum crystal misorientation, which for a cubic 2Dlattice is π/ (45 o ).To the multi-well function in Eq. (3) we add a chemical Well o Liquid θ Well o Well o Well
Crystal 22.5 o Well
Crystal 78.75 o Well o φ X φ Y FIG. 1: Landau free energy landscape in the space of a vectororder parameter with N wells = 8 and a solid-solid interactionstrength parameter b = 0 . . The energy minimum in thecenter corresponds to the liquid well, while each of the N wells correspond to a distinct solid orientation. In this work,the dynamics of the order parameter are written in term of itsCartesian projections φ X = R cos( θ ) and φ Y = R sin( θ ) .energy contribution of the form ω ch ( µ, R ) = P ( R ) ω s ( µ, T ) + (cid:16) − P ( R ) (cid:17) ω (cid:96) ( µ, T ) , (4)where µ = µ ( (cid:126)x, t ) is the chemical potential field, ω s and ω (cid:96) are the grand potential densities of the solid and liquid phases,respectively, and P ( R ) interpolates the thermodynamic po-tentials across bulk states. It satisfies P (0) = 0 and P (1) = 1 .Eq. (4) provides the chemical driving force for solidificationwritten in terms of the chemical potential, which is the naturalthermodynamic variable of the grand potential functional.The last contribution to the grand potential energy requiresgradient energy terms to account for non-local interactions. Interms of the Cartesian representation of the order parameter,we also define gradient interaction term as according to f int = 12 W X |∇ φ X | + 12 W Y |∇ φ Y | , (5)where { W X , W Y } can be shown to be proportional to gra-dient energy and inversely proportional to the square root ofthe nucleation barier [14]. The coefficients must be madeanisotropic in order to self-consistently reproduce the sym-metries of crystal growth in solidification. The anisotropy ofthese coefficients will be discussed further below.Combining Eqs. (3), (4) and (5) yields a grand potential functional of the form H Ω[ (cid:126)φ, µ, T ] = (cid:90) V (cid:32) W X |∇ φ X | + 12 W Y |∇ φ Y | + R − b cos ( N wells θ )1 + b R + R + 1 H P ( R ) ω s ( µ, T )+ 1 H (cid:16) − P ( R ) (cid:17) ω (cid:96) ( µ, T ) (cid:33) dV, (6)where H is the nucleation barrier, T is the temperature and P ( R ) is a function that interpolates the grand potential be-tween its bulk solid to liquid values. An explicit form of P ( R ) introduced by Wheeler-Boettinger-McFadden [15] is P ( R ) = 10 R − R + 6 R , although several other formsare possible. B. Phase field dynamics
The dynamical evolution for the vector order parameteris governed by the usual dissipative relaxation dynamics ineach component field, while solute evolution follows the usualmass conservation dynamics. The explicit dynamical equa-tions for these fields are given by ∂φ K ∂t = − ∂ Ω ∂φ K + η K , K = { X, Y } (7) ∂c∂t = ∇ · (cid:16) M ( (cid:126)φ, c ) ∇ µ (cid:17) , (8)where K = { X, Y } denote the two components of (cid:126)φ , η K arenoise sources in each component (discussed further below),while c ( (cid:126)x, t ) is the solute concentration field of a binary alloy,which in the grand potential ensemble is determined via thechemical potential field µ ( (cid:126)x, t ) according to the relation c = − δ Ω δµ = P ( R ) c s ( µ ) + (cid:0) − P ( R ) (cid:1) c (cid:96) ( µ ) , (9)where c s = ∂ω s /∂µ and c (cid:96) = ∂ω (cid:96) /∂µ . Eq. (9) can be seen asa constraint between c and µ . Combining Eq. (9) with Eq. (8)yields the following evolution equation for µ ( (cid:126)x, t ) directly, χ ( µ ) ∂µ∂t = ∇· (cid:16) M ( (cid:126)φ, c ) ∇ µ (cid:17) − P (cid:48) ( R ) (cid:2) c s ( µ ) − c (cid:96) ( µ ) (cid:3) ∂R∂t , (10)where the susceptibility field χ ( µ ( (cid:126)x, t )) is given by χ = P ( R ) ∂c s ( µ ) ∂µ + [1 − P ( R )] ∂c (cid:96) ( µ ) ∂µ (11)Either Eq. (10) or Eq. (8) can be coupled to the order parame-ter component dynamcis of Eq. (7).In what follows, we will be focusing on examining the ef-ficiencies and self-consistencies afforded when using a vectororder parameter in favor of multiple scalar order parametersor phase fractions in phase field modeling. As a result, wewill limit the system studied to a dilute binary alloy. We willsimulate phase field dynamics using Eqs. (8) and (7). Fur-thermore, in order to achieve quantitative solidification re-sults in the limit of a diffuse interface, we will incorporatean anti-trapping flux in Eq. (8), whose role is to suppress spu-rious kinetic behavior cause by the use of a diffuse interface[16, 17]. It is noted that the dynamical evolution equationfor mass transport no longer follows a variational formulationwhen anti-trapping is used.Taking into account the above considerations, we proceedas follows to arrive at an evolution equation for the vector or-der parameter components of our model: first we substitutethe grand potential densities for the solid and liquid phasesof an ideal binary alloy into Eq. (6); secondly, we carry outthe variational derivatives as specified in the first of Eqs. (7);thirdly, we simplify the chemical driving force analogously toRef. [13]. This yields, τ K ∂φ K ∂t = ∇ · (cid:18) ∂f int ∂ ( ∇ φ K ) (cid:19) − ∂f int ∂φ K + (cid:0) R − R (cid:1) φ K + 6 b ( cos ( N wells θ ) − b Rφ K + 2 bN wells sin ( N wells θ )1 + b R ∂θ∂φ K − λ − k e (cid:18) e u − − T (cid:96) − T | m e(cid:96) | c (cid:19) (1 − R ) φ K + η K for K ∈ { X, Y } , (12)where u is a dimensionless chemical potential differencerelative to the equilibrium chemical potential for a bi-nary alloy, expressed in terms of concentration as u =ln [ c/ ( c · (1 − (1 − k e ) R ))] [14]. The first two terms onthe right-hand side of Eq. (12) represent the interface energypenalty, whose expressions are rather lengthy, and are thusshown in the Appendix. τ K is the anisotropic time constant, λ is the coupling constant, which is used to control the model’sconvergence onto the key results of the sharp interface modeland k e is the equilibrium partition coefficient. We have alsotacitly appended into the chemical driving force of Eq. (12) atemperature component, where T (cid:96) is the liquidus temperature, T is the local temperature, m e(cid:96) is the liquidus slope, and c isa reference alloy concentration. The last term in Eq. (12), η K ,is thermal noise, discussed further below.The concentration evolution follows from Eq. (8) and mapsidentically onto the dilute binary alloy model of Refs. [16, 17],with R replacing the solid-liquid order parameter, but where R is scaled here to the limits [0 , . Explicitly, the concentra- tion equation becomes ∂c∂t = ∇· (cid:34)(cid:32) h ( R ) D s k e + D (cid:96) (cid:16) − h ( R ) (cid:17)(cid:33) ∇ e u + W a t c (1 − k e ) e u ∂R∂t ∇ R |∇ R | (cid:33)(cid:35) , (13)where h ( R ) = R is an interpolation function, D S ( D L ) is thesolute diffusion coefficient in solid (liquid), and a t is the anti-trapping current prefactor. Its value is listed in Table I andis determined by conducting a matched asymptotic analysisusing the methods in Refs. [14, 17].A key feature of quantitative phase field modeling is theuse of matched interface asymptotic analysis to determine theparameter relations that allow the model to emulate the in-terfacial kinetic of the classical sharp interface model of so-lidification. There have been several important works on thistopic for models using a scalar order parameter [14, 17–19].The outcomes are similar; they lead to relationships betweenmicroscopic phase field parameters ( τ , W and λ ) and themacroscopic material parameters of the sharp interface model,namely, the capillary length ( d ) and kinetic coefficient ( β ).It turns out that we can directly apply the aforementionedmatched asymptotic analysis methods, developed for scalarorder parameter PF models, in order to inter-relate the param-eters τ K , W K and λ in our model. The rationale in doing sois as follows: if we ignore anisotropy (as is done in the theabove-referenced approaches) and only consider single indi-vidual order parameter profiles going from the liquid well intoany of the solid wells of our model (see Fig. 1), we can map τ K → τ and W K → W . Further, in this situation, the Lan-dau angle variable becomes a constant everywhere, equal to θ = m π/N wells where m is an integer. This guarantees theelimination of the cubic terms in the Landau contribution tothe energy. This reduces Eq. (12)) identically to the traditionalform of the scalar order parameter evolution equation in bothorder parameter components, φ K , K = X, Y . As a result, thematched interface asymptotics in Refs. [14, 17, 19] is exactlyapplicable. It should be noted that this assumes that the orderparameter profile can “ride the rails” along a line between theliquid well to any of the solid wells. This is easily satisfied,and any error is negligible in our single crystal benchmarktests presented in the Section III B. These considerations leadto the following relations, d = a W λ (14) β = a τ λ W (cid:20) − a λW τ D (cid:96) (cid:21) , (15)where a and a are asymptotic analysis constants listed inTable I, where it is notable that solid-liquid ordering is in theinterval [0 , , while perhaps a more common practice is toscale the order parameter for the interval [ − , .To make the model anisotropic, the phase field parameters W K and τ K are appended with an anisotropy correction ac-cording to W K = W a K ( φ X , φ Y ) , and assuming that β = 0 , τ K = τ a K ( φ X , φ Y ) (assuming β = 0 ), where a k is ananisotopy function. For 2D cubic crystals, a k is defined by a K ( φ X , φ Y ) = 1 + (cid:15) c cos (4 θ ∇ φ K − θ ) , K ∈ { X, Y } , (16)where (cid:15) c is the anisotropy strength. It should be noted thatin Eq. (16), the Landau angle θ is not multiplied by four toachieve the correct rotation for a cubic 2D lattice, as visu-alized in Fig. 1 . Here, θ ∇ φ K is the local solid-liquid inter-face normal angle for vector order parameter component φ K ( K ∈ { X, Y } , given by θ ∇ φ K = arctan (cid:18) ∂ y φ K ∂ x φ K (cid:19) , K ∈ { X, Y } (17)where θ sets the reference Landau angle, i.e., θ = arctan (cid:18) φ Y φ X (cid:19) . (18)The last term, η K for K ∈ { X, Y } , represents thermal fluc-tuations, which must satisfy the fluctuation-dissipation theo-rem in order to assure convergence of the phase field equa-tions to thermodynamics equilibrium. We follow the approachgiven in Ref. [20] to scale these noise sources quantitatively,which yields (cid:104) η K , η K (cid:105) = 2 k B T Hτ δ ( x − x (cid:48) ) δ ( t − t (cid:48) )= 2 k B T τ a W Jd c g (cid:96)c,c x t , (19)where k B is Boltzmann’s constant. On the second line we ap-proximate the delta functions by ∆ x and ∆ t following Ref.[20]. Furthermore, in the second line of Eq. (19), we substi-tuted H with an expression from the matched interface asymp-totics for dilute binary alloys [14, 17], where J = 1 / is ascaling factor (specific to the interpolation function we use forthe driving force), and g (cid:96)c,c is the liquid free energy curvature.In addition, the noise term in Eq. (19) is filtered in space andtime so that the fluctuations occur approximately in the lengthscale of W and time scale τ , as these can be considered asthe smallest physical scales in solidification that the currentphase field model can resolve. More details about this filter-ing scheme can be found in the Supplemental material [21]. C. Alternate gradient formulation
It is noted that as in Morin et al. [12], we chose a gradientenergy penalty in Eq. (5) that is a sum of the contributing com-ponents, and our coefficients are made anisotropic to satisfy acrucial criterion required for solidification of metals. Anotherperhaps more consistent choice would be to represent the en-ergy penalty in terms of polar coordinates R and θ , in the form W |∇ R | + 12 W R |∇ θ | . (20) While this form reproduces some features of solidification andcoarsening qualitatively, this R − θ representation suffers fromlattice pinning due to topological defects, analogous to thoseseen in the orientation PF models [7]. Furthermore, it is un-clear how to implement thermodynamically consistent noisefor the θ field in our vector order parameter model. As suchthis alternate formulation is not considered further in this workand will be pursued in future work. III. MODEL BENCHMARKS
In what follows, we demonstrate that our vector order pa-rameter model reproduces certain established benchmark re-sults of solidification and solid-state coarsening phenomena.In particular, we first analyze single crystal dendrite growthkinetics. We then conduct an analysis of grain boundary pro-files, their energies, and their coalescence. Following this,we verify that the polycrystalline coarsening for Model A dy-namics (i.e. no concentration diffusion) is consistent with the-oretical expectations. Finally, we demonstrate nucleation andgrowth of a polycrystalline alloy system using thermal noisefluctuations of the vector order parameter.
A. Numerical implementation
All phase field simulations we present are performed usingthe computational platform presented in Ref. [22], which con-tains adaptive mesh refinement (AMR) and distributed mem-ory parallelization based on MPI.Concentration evolution in Eq. (13) was implemented us-ing a finite volume method, as well as the divergence termsin the order parameter equations of Eq. (12). Other gradi-ent terms were expressed via simple five stencil finite differ-encing. Neumann, i.e. zero flux, boundary conditions wereapplied, unless otherwise stated. We use forward Euler timeintegration, where it noted that the forward Euler time step ofthe stochastic noise term η K in Eqs. (12) and (19) is propor-tional to √ ∆ t . B. Convergence analysis for dendritic solidification
Since our vector order parameter model relies on the samediffuse interface analysis as Refs. [16, 17], we verified that itsdynamics reproduces the dendritic solidification results pre-dicted by the classical sharp interface model of solidification.To do so, we initialized a seed with radius d in a uniformsupersaturation Ω = c/ ( c (1 − (1 − k e ) R ) = 0 . . We thenperformed a reference simulation using the scalar order pa-rameter model presented in Ref. [16] and compared it to asimulation conducted with our model with N wells = 8 anda solid-solid interaction barrier parameter b = 4 . . We set W = d / . and ∆ x = 0 . W in our simulations.The solid seed is initialized into the second well of ourLandau energy landscape (Well o relative to x-axis.TABLE I: Phase field model parameters for Al-4.5at%Cu.Partition coefficient k e m e(cid:96) -5.3 K/wt%Alloy concentration c Γ × − K mSolutal capillary length d g (cid:96)c,c × J/m a Liquid diffusion coefficient D (cid:96) × − m /sSolid diffusion coefficient D s × − m /s b Kinetic coefficient β (cid:15) c W d / . Asymptotics constant a a a t √ Mesh spacing ∆ x W / 0.2 W / 0.6 W c Time step size ∆ t ∆ x / (6 D (cid:96) ) a : Value from Thermo-Calc TCAL database. b : D s is non-zero only for back-diffusion analysis in Fig. 5. c : We use different ∆ x depending on the sections: ∆ x = 0 . W for the single crystal solidification benchmark (Section III B), ∆ x = 0 . W for 1D grain boundary analyses(Section III C), and ∆ x = 0 . W for polycrystallinesimulations (Sections III D and III E).We applied the same lattice rotation to the standard phasefield model of Karma [16]. The resulting dendrite morphol-ogy and center-line concentration through to the dendrite tipare shown in Fig. 2 at a time t = 480 τ . The results showexcellent agreement between the standard scalar order param-eter alloy model and our vector order parameter model. It isnoted that these results do not change with a decrease of thesolid-solid barrier parameter to b = 0 . . We also verified thatthe tip speeds are indistinguishable (not shown here). It isnoted that an excessively large interaction parameter b leadsto unphysical oscillations inside the solid. These can be sup-pressed by decreasing ∆ t . Thus, for N wells = 8 , simulationswith b = 0 . are stable for ∆ t = 0 . x / (6 D (cid:96) ) , while for b = 4 . t = 0 . x / (6 D (cid:96) ) is required to maintain stability. C. 1D grain boundary analysis
This section analyses grain boundaries and their energy inthe context of the vector order parameter model. We per-form simulations in one spatial dimensional to model a 1Dgrain boundary analysis using so-called Model A dynamics,i.e. only the vector order parameter evolution of Eq. (12) isconsidered, while the concentration dynamics in Eq. (13) isswitched off.We applied a fixed value boundary condition so that the farends of the 1D grain boundary geometry satisfy φ X + φ Y = 1 ,and such that θ = arctan( φ Y /φ X ) corresponds to the ap-propriate crystal orientations. The grid spacing was set to ∆ x = 0 . W , even though ∆ x = 0 . W also yields conver-gent results. It noted that the use of larger grid spacing, such Distance from center [dx] C o n c e n t r a t i o n [ c ] N=8, b=4.0Karma01.0 ]
50 100 150 200 250
FIG. 2: Binary alloy dendrite growth benchmark. Centre-lineconcentration computed with the scalar PF model of Ref. [16](dashed line) compared to corresponding result from the newvector order parameter model with N wells = 8 and b = 4 . (solid line). In both cases, the grain is rotated 11.25 o , W = d / . , the initial seed radius is d and the initialsupersaturation Ω = 0 . . (Bottom right inset) dendritesolute contours. The profiles are recorded at time t = 480 τ .as ∆ x = 0 . W , may lead to lattice pinning in the multidi-mensional polycrystal coarsening simulations, which are dis-cussed in Section III D. Grain boundary equilibration was con-ducted by initializing two 1D “solid blocks” corresponding totwo different solid wells, with hyperbolic tangent tails over-lapping at R = 0 . . Approximately time steps were re-quired to equilibrate the profiles into a static grain boundary.The final equilibrium profiles were found to be independentof the details of the profile initialization.The resulting profiles are shown in Fig. 3 for a solid-solidbarrier parameter b = 4 . and b = 0 . , along with the cor-responding vector order parameter trajectories defining eachgrain boundary in the Landau free energy landscape (super-imposed as a black solid line). For each choice of b , the driv-ing force was set to either zero, i.e. ∆ := ( T l − T ) / ((1 − k e ) | m el | c ) = 0 , or ∆ = 0.55. In all cases, solid-liquid or-dering, represented by R , develops a typical cusp-like mor-phology in the grain boundary region consistent with earlierpolycrystalline PF models, such as in orientation field models[7]. Moreover, the solid-solid grain boundary ordering de-creases as the misorientation increases, approaching a “wet”grain boundary at high misorientation, as expected. The useof a smaller solid-solid barrier parameter b = 0 . increasesthe ordering at the grain boundary, which is consistent witha smaller energy barrier between the solid wells. Finally, ap-plying a driving force ( ∆ = 0.55) increases the ordering at thegrain boundary as expected.In Fig. 3d, one can see slight oscillations in the solid-liquidordering across the grain boundaries GB1,3 and GB1,6. Thisis expected for the small solid-solid barrier coefficient ( b =0.1) and large quench ( ∆ = 0.55), as it is energetically fa-vorable for the order parameter profile to approach the inter-mediate solid well, as shown in the grain boundary trajecto-ries (black lines) in Fig. 3d. Curiously, this type of solid-liquid ordering oscillation is also found in phase field crys-tal simulations of Mellenthin et al [23]. However, we expectthat these solid-liquid ordering oscillations are uncommon formetals and should be generally avoided by setting the solid-solid barrier to a sufficiently large value, especially when thecomputational interface width W is increased to thin inter-face conditions.We next evaluated grain boundary energies, the results ofwhich are shown in Fig. 4. The grain boundary energy is de-fined as the excess energy at the interface, which for 1D sys-tems can be shown to have the form γ GB = (cid:90) ∞∞ (cid:32) (cid:15) (cid:16)(cid:16) ∂∂x φ X (cid:17) + (cid:16) ∂∂x φ Y (cid:17) (cid:17) + H (cid:16) R − b cos( N wells θ φ )1 + b R + R (cid:17) + P ( R )∆ (cid:33) dx, (21)where H = a W d c g (cid:96)c,c ,(cid:15) = W √ H, ∆ = T l − T (1 − k e ) | m el | c . (22)For an ideal amorphous grain boundary, an estimate of theenergy can be derived analytically; the details are providedin the Supplemental material [21]. An ideal amorphous grainboundary assumes that the Landau angle θ is constant exceptat the center of the grain boundary, where the ordering goesto zero. Without a driving force ( ∆ = 0 ), this energy has aclosed-form expression that is twice the solid-liquid interfaceenergy given by γ amorpGB = 2 · γ SL = W H √ ≈ . J/m , (23)where parameters are taken from Table I. When an undercool-ing (here ∆ = 0 . ) is applied, the grain boundary energyneeds to be numerically computed, yielding γ amorpGB (∆ =0 . ≈ . J/m (see the Supplemental material [21]for more details). These estimated grain boundary energies( γ amorpGB ) provide ideal limits for comparison and are shownas horizontal lines in Fig. 4. These energies are comparedfor the Model A-equilibrated grain boundary profiles shownin Fig. 3, where the excess grain boundary energies are cal-culated by numerically evaluating Eq. (21), and the resultingenergies are shown by the scatter data in Fig. 4.The grain boundary energies obtained from our model arein good agreement with previous phase field crystal simula-tions of Mellenthin et al. [23], which are characterized by three main features. First, without a driving force (filled mark-ers in Fig. 4), the smallest misorientations (up to misorienta-tion ∼ . o ), produce slightly attractive grain boundaries.Second, large misorientations up to the symmetry point forcubic crystals leads to grain boundary energies identical tothe ideal grain boundary energy, corresponding to zero ap-plied driving force ( γ amorpGB (∆ = 0 )). Third, adding a drivingforce of ∆ = 0 . (hollow markers in Fig. 4) increases thegrain boundary energy in a gradual manner as a function ofmisorientation.As expected, an increase in solid-solid barrier parameter b leads to a corresponding increase in grain boundary energy forall cases, except for those cases where the energy has alreadysaturated to it’s limiting value, i.e. when the grain boundaryhas saturated to a large misorientation (for ∆ = 0 misorien-tations 22.50-45 o and its corresponding crystalline symmetricvalues, while for ∆ = 0.55 similar saturation is not found).When we consider a thermodynamic driving force appliedto the grain boundary, the simulated Model A energies at themaximum misorientation are smaller than our analytical esti-mate (dashed horizontal line in Fig. 4). This discrepancy isdue to the fact that in the analytical estimate, we assume thatsolid-liquid order is zero in the grain boundary region, and that θ varies only at zero ordering, neither of which are well satis-fied for ∆ = 0.55, particularly for small b ( b = 0.1), which, asdiscussed above, leads to small interfacial oscillations withinthe interface.
1. Grain boundary back-diffusion
This section demonstrates that the solute evolution, lead-ing to the development of cohesive solid networks and grainboundaries, follows behavior consistent with the seminal ofRappaz and co-workers [24]. For this part of the study, weactivate the concentration dynamics of Eq. (13), and apply aconstant cooling rate × K/s. In addition, we allow forenhanced back-diffusion by assuming that the solid diffusioncoefficient is D s = D (cid:96) / . Here, the two solid regions wereinitialized with a thick liquid layer in between them to emu-late, in 1D, the conditions of grain-grain merger during solid-ification.Back diffusion was simulated in this grain boundary sys-tem having the smallest misorientation allowed by our model(labeled GB 1,2 in Fig. 4), corresponding to a misorientationof 11.25 o . The data was collected at two solid-solid barrierparameters, b = 0.1 and b = 4.0. The snapshots of the phasefield profiles and the corresponding concentration profiles areshown in Fig. 5a, where time increases from top to bottom.Fig. 5b plots the peak grain boundary concentration as a func-tion of time superimposed on the dilute Al-Cu phase diagram.At early times, when the tails of the two grains are not in con-tact, the concentration follows the liquidus line. As the twograins begin to “feel” one another through their order param-eter tails, and start to interact, the concentration starts to droptowards the solidus concentration. Due to back-diffusion, theyasymptotically start to approach the average alloy concentra-tion (vertical dashed line). This behavior is consistent with a) b) G B , G B , G B , b=4.0, Δ =0 b = 0.1, Δ =0 G B , G B , G B , c) d) G B , G B , G B , b=4.0, Δ =0.55 b = 0.1, Δ =0.55 G B , G B , G B , Position [ μ m] Position [ μ m] --- - - - -- - -- - ------- - - -- - - FIG. 3: Grain boundary profiles and vector order parameter trajectories in Landau landscapes with zero driving force ∆ = 0 (a-b) and a relatively large driving force ∆ = 0 . (c-d). Two different solid-solid barrier heights are considered ( b = 0 . and b = 0 . ). Concentration dynamics are not considered here.the back-diffusion results reported by Rappaz et al. [24] andOfori-Opoku and Provatas [3]. Using a large solid-solid bar-rier parameter, b = 4.0 (red circles in Fig. 5b), decreases thesolid-solid ordering around the grain boundary for a given un-dercooling, and thereby increases the concentration peak, asexpected.Overall, the above 1D grain boundary back-diffusion studyshows the expected behavior in terms of solid-solid order- ing, grain boundary energy behavior, grain boundary mergerand solute distribution at the grain boundary. It is also note-worthy that the properties reproduced by the proposed vec-tor order parameter model have also all been exhibited by themore fundamental phase field crystal model, as documentedin Ref. [23], as well as in more traditional phase field mod-els [3, 24]. GB Labels 𝛾 𝐺𝐵𝑎𝑚𝑜𝑟𝑝 (Δ = 0) 𝛾
𝐺𝐵𝑎𝑚𝑜𝑟𝑝 (Δ = 0.55)
FIG. 4: Grain boundary energy versus grain boundarymisorientations for the case of N wells = 8 . Grain boundarieswere grown with Model A dynamics (no concentrationdynamics) for two values of undercooling, ∆ = b = 0 . or b = 4 . . Analytical grain boundary energies, γ amorpGB (∆) , are shown as horizontal lines, where the solidline corresponds to ∆ = 0 and the dashed line to ∆ = 0 . . D. Polycrystalline grain coarsening
In this section, we verify that the Model A dynamics of theproposed vector order parameter model satisfies the theoret-ically established ∼ t / coarsening kinetics. For this study,we set the grid size ∆ x = 0 . W as we found that larger valueof mesh size (e.g., ∆ x = 0 . W ) can lead to lattice pinning,especially when a large driving force is applied. Simulationswere performed using a driving force ∆ = 0 . , an interac-tion barrier strength of b = 4 . , and we applied zero gradientboundary conditions. For this study, we initialized the sim-ulation domain with a random Voronoi tessellation with fullordering R = φ X + φ Y = 1 within the bulk of each grain,as shown in the first column of Fig. 6a. The other columnsof Fig. 6a show the time evolution and grain coarsening inthe polycrystalline system, where each row highlights differ-ent measures of the vector order parameter. It is noteworthythat the regions corresponding to high-angle grain boundarieslead to the emergence of solid-liquid ordering.To verify that the grain coarsening follows the well-knownresults of linear (parabolic) growth with area (radius) versustime, we evaluated the average grain size using the so-calledintercept method, which is typically used in experimental set-tings to determine grain sizes. In this method, a line is drawnrandomly across the microstructure and the size of the grainsintercepting this line is measured. A sufficiently large num-ber of these line-based measurements, at various orientationsthrough the microstructure data, leads to a good estimate ofthe grain size d grain . We converted the so-calculated aver-age grain size to an equivalent circle via π ( d grain / . Thegrain area versus time is shown in Fig. 6b. Mean-field argu- ments, stating that grains grow by mean curvature, dictate thatthe grains should coarsen such that the average grain area in-creases linearly with time. This linear coarsening regime canbe seen in Fig. 6b, where a dashed line is drawn as a guide forthe eye for the linear regime. As the grain size approaches thesystem size, the mean field assumption breaks down and thescaling relation no longer holds. E. Nucleation and polycrystalline solidification in an alloy
This section demonstrates noise-induced nucleation andpolycrystalline solidification in a binary alloy using the pro-posed vector order parameter model. The aim here is to high-light the robustness of the model and the self-consistencyof applying thermal noise acting on the vector order param-eter defined everywhere. To prevent interference with theevolution of actual ordered crystals, we threshold the noiseterms such that it is only active in the vicinity where R = (cid:112) φ X + φ Y < . . We additionally filter the noise in spaceand time, so that the fluctuations occur roughly over a lengthscale of W (not ∆ x ), and over the time scale τ (not ∆ t ). Thedetails of these filtering procedures are given in the Supple-mental material [21]. These filtering procedures were neces-sary to ensure the numerical stability of the system, and fluc-tuations below the length and time scales ( W and τ ) have noclear physical meaning in coarse grained field theories. Forthe simulations shown in this section we used periodic bound-ary conditions and a cooling rate of . × K/s, and initial-ized the system as a uniform liquid.Fig. 7 shows the nucleation and subsequent time evolutionof a polycrystalline network of grains in a binary alloy, start-ing from a liquid state and with thermal fluctuations appliedover time. Time goes from left to right as the correspondingundercooling ( ∆ ) grows. The rows display the concentrationfield and different measures of the vector order parameter. Thedata shows that the thermal fluctuations grow over time as thedriving force increases, leading to homogeneous nucleation ofglobular structures, i.e. the nuclei, as shown in the first col-umn of Fig. 7. The nuclei then start to merge and coalesceinto a polycrystalline network, which starts to coarsen overtime. The grain boundaries in this set of simulations appearjagged and noisy due to thermal fluctuations, but eventuallydo smooth out (coarsen) due to surface energy minimizationat later times. IV. DISCUSSION
In the proposed vector order parameter phase field model,the number of solid wells, i.e. the number of unique grains,is practically limited to a value of roughly N wells = 10, how-ever this can be increased easily if the grid spacing ∆ x is de-crease. A more practical and simple technique for increasingthe number of grains that avoids the excessive computationalcost associated with decreasing ∆ x would be to increase thevector order parameter dimensionality from d = φ X and φ Y ) to e.g. d = 3 or 4 ( φ K , K = 1, ..., d ).0 a) t = 0 μ st = 20 μ st = 250 μ s Distance [W ] b = 0.1 b = 4.0 b = 4.0 Average b) FIG. 5: Grain boundary coalescence for a cooling rate K/s, showing the evolution of the corresponding order parameter andconcentration profiles. a) Solid-liquid ordering and concentration at three time steps for b =4.0, where concentration (C) is inunits of the average composition c = 4.5 at%Cu; b) Grain boundary peak concentration at various time steps for b = b = N wells )needs to be rather small to describe phases with notably dif-ferent solubilities, such as an α FCC aluminium phase and anintermetallic Al Cu phase.In order to reduce computational costs associated with theuse of multiple order parameters to represent distinct phases,multi-phase field model or multi-order parameter phase fieldmodels can use so-called re-indexing algorithms as a wayto re-assign grains to avoid artificial melding of neighboringgrains [11]. This re-indexing might not be possible for thecurrent vector order parameter model, as the order parametersof neighboring grains are directly interacting with each otherand relax to profiles such as those shown in Fig. 3. However,the use of higher dimensional vector order parameter fields,coupled with adaptive mesh refinement should be able to pro-vide ample degrees of freedom to simulate polycrystalline net-works relatively efficiently in the context of a single physicalorder parameter.The anisotropy of grain boundaries can be an importantfeature in polycrystalline coarsening, characterized in 3D by three misorientation angles and two inclination angles [26]—in a 2D system by two misorientation angles and a single incli-nation angle. In its current form, the present model dependsmostly on the misorientation, i.e. on which solid wells arein contact. There is a minor inclination dependence throughthe anisotropic solid-liquid interface energy. We expect thatthe inclination dependence can be increased by, for example,adding an anisotropic energy term proportional to | θ | for theenergy expression in Eq. (6), analogously to how it is done forthe solid-liquid interface energy.Physically, grain boundaries move significantly slower thansolid-liquid interface. As grain boundaries are necessarily as-sociated with a change in orientation, the mobility decreasecould be implemented by adding an orientation gradient |∇ θ | dependence to the phase field time scale, for example in theform τ K → (1 − exp( − A |∇ θ | )) τ K , where A controls howstrongly the mobility decays around grain boundaries.Numerically, the present vector order parameter model re-quires a relatively dense grid, with roughly ∆ x ∼ . W , inorder to resolve the grain boundary profiles and to avoid meshpinning. This becomes more arduous when the driving forceor undercooling is increased since this leads to steeper pro-files, a feature also present to some extent in all phase fieldmodels. This requires further investigation to find possiblenumerical remedies.The vector order parameter model can be more elegantly1 μs μs μs μs R a) b) 512 dx (4.4 μ m) FIG. 6: Polycrystalline coarsening using the proposed vector order parameter model, constrained to Model A (no solutedynamics), with an undercooling of ∆ = 0 . and solid-solid barrier b = 4 . . a) Time evolution of polycrystalline grainnetwork, where time increases from left to right. The rows show different measures of the order parameter; b) Solid fractionversus time, with the linear scaling regime is shown by the dashed line as a guide for the eye.formulated in polar coordinates R, θ (see Fig. 1). However, inour tests, the system is prone to mesh pinning due to topolog-ical defects analogical to those seen by Korbuly et al. [7] forthe orientation phase field model. For those readers interestedin purusing this formulation, Korbuly et al. present techniquesand remedies to address these issues that are expected to alsobe applicable to the present model.
V. CONCLUSION
The vector order parameter phase field model presented inthis work is introduced as a practical tool for the study ofmicrostructure evolution in polycrystalline solidification, andwhich reproduces key benchmarks quantitatively. The specificbenchmarks that were demonstrated in this work are: dendritesolidification [16, 17], grain boundary energy versus misori-entation and undercooling, coalescence and segregation asso-ciated with back-diffusion and nucleation-induced solidifica-tion and coarsening in binary alloys.Future work with this new phase field model will explore2 μ s 61 μ s 205 μ s 409 μ s 819 μ s Δ = 0.6 Δ = 1.0 Δ = 2.1 Δ = 3.6 Δ = 6.61-11-110 R
512 dx (4.4 μ m) C FIG. 7: Thermally induced homogeneous nucleation, polycrystalline solidification and coarsening in a binary alloy simulatedwith the proposed vector order parameter phase field model. Model parameters are taken from Table I ( ∆ x = 0 . W ). Thecooling rate is set to . × K/s. Time flows from left to right as the corresponding undercooling increases. The rowsindicate the fields of the model.the addition of higher dimensional vector order parameterfield in order to increase the efficiency increasing the num-ber of unique grains, and to introduce new thermodynamicallydistinct phases, thus making possible the efficient simulationof multi-phase solidification, such as for example the simula-tion of intermetallic phases in the Al-Cu system.
ACKNOWLEDGEMENT
NP acknowledges the National Science and EngineeringResearch Council of Canada (NSERC) and the Canada Re- search Chairs (CRD) Program. TP and AL acknowledgethe support of Academy of Finland through the HEADFOREproject, Grant No. 333226.
Appendix A: Variational derivative of interface energy term
This section shows the explicit expression for the interfaceenergy term for the vector order parameter model, Eq. (5).Let’s first separate the anisotropy function from interfacewidth as W K = W a K . Then the gradient penalty free en-3ergy contribution for order parameter model is (cid:90) V f int dV = (cid:90) V (cid:88) K = X,Y W | a K ( θ ∇ φ K , θ ) ∇ φ K | dV, (A1)where a K ( φ X , φ Y , ∇ φ K ) = 1 + (cid:15) c cos (4 θ ∇ φ K − θ ) , (A2)where the local interface normal angle is θ ∇ φ K = arctan (cid:18) ∂ y φ K ∂ x φ K (cid:19) , (A3)and the reference angle of order parameter vector (cid:126)φ =( φ X , φ Y ) is θ = arctan (cid:18) φ Y φ X (cid:19) . (A4)For order parameter vector components K ∈ { X, Y } , the to-tal variational derivative of this term is then δδφ K (cid:88) K (cid:48) = X,Y | a K (cid:48) ∇ φ K (cid:48) | = (cid:18) ∂∂φ K − ∇ · ∂∂ ∇ φ K (cid:19) (cid:88) K (cid:48) = X,Y | a K (cid:48) ∇ φ K (cid:48) | = a X ∂a X ∂φ K |∇ φ X | + a Y ∂a Y ∂φ K |∇ φ Y | − ∇ · (cid:18) a K ∇ φ K + a K |∇ φ K | ∂a K ∂ ∇ φ K (cid:19) , (A5) where ∂∂ ∇ φ K := (cid:80) x,y,zi ˆ e i ∂∂ i ( ∂φ K ) , and the following iden-tities are used: ∂θ∂φ X = − φ Y | φ | , ∂θ∂φ Y = φ X | φ | , (A6) ∂θ ∇ φ K ∂ ( ∂ x φ K ) = − ∂ y φ K |∇ φ K | , , ∂θ ∇ φ K ∂ ( ∂ x φ K ) = ∂ x φ K |∇ φ K | , (A7) a (cid:48) K ( θ ∇ φ K , θ ) = − (cid:15) c sin(4 θ ∇ φ K − θ ) , (A8) ∂a K ∂ ( ∂ x φ K ) = a (cid:48) K ( θ ∇ φ K , θ ) ∂θ ∇ φ K ∂ ( ∂ x φ K )= − a (cid:48) K ( θ ∇ φ K , θ ) ∂ y φ K |∇ φ K | , (A9) ∂a K ∂ ( ∂ y φ K ) = a (cid:48) K ( θ ∇ φ K , θ ) ∂θ ∇ φ K ∂ ( ∂ y φ K )= a (cid:48) K ( θ ∇ φ K , θ ) ∂ x φ K |∇ φ K | , (A10) ∂a K ∂φ K = (cid:15) c sin (4 ( θ ∇ φ K − θ/ ∂θ∂φ K = − a (cid:48) K ( θ ∇ φ K , θ ) ∂θ∂φ K . (A11) [1] I. Steinbach, F. Pezzolla, Physica D 134 (1996) 385.[2] B. Bottger, J. Eiken, I. Steinbach, Acta Materialia 54 (2006)2697.[3] N. Ofori-Opoku, N. Provatas, A quantitative multi-phase fieldmodel of polycrystalline alloy solidification, Acta Materialia58 (6) (2010) 2155 – 2164.[4] R. Kobayashi, J. A. Warren, W. C. Carter, Vector-valued phasefield model for crystallization and grain boundary formation,Physica D: Nonlinear Phenomena 119 (3-4) (1998) 415–423.[5] L. Gr´an´asy, T. B¨orzs¨onyi, T. Pusztai, Crystal nucleation andgrowth in binary phase-field theory, Journal of Crystal Growth237 (2002) 1813–1817.[6] J. A. Warren, R. Kobayashi, A. E. Lobkovsky, W. C. Carter,Acta Materialia 51 (2003) 6035.[7] B. Korbuly, M. Plapp, H. Henry, J. A. Warren, L. Gr´an´asy,T. Pusztai, Topological defects in two-dimensional orientation-field models for grain growth, Physical Review E 96 (5) (2017)052802.[8] T. Pusztai, G. Bortel, L. Gr´an´asy, Phase field theory of poly-crystalline solidification in three dimensions, EPL (EurophysicsLetters) 71 (1) (2005) 131.[9] I. Steinbach, F. Pezzolla, A generalized field method for mul- tiphase transformations using interface fields, Physica D: Non-linear Phenomena 134 (4) (1999) 385–393.[10] H. Salama, J. Kundin, O. Shchyglo, V. Mohles, K. Marquardt,I. Steinbach, Role of inclination dependence of grain bound-ary energy on the microstructure evolution during grain growth,Acta Materialia (2020).[11] C. E. K. III, L.-Q. Chen, Acta Materialia 50 (2002) 3057.[12] B. Morin, K. R. Elder, M. Sutton, M. Grant, Model of the ki-netics of polymorphous crystallization, Physical Review Letters75 (11) (1995) 2156.[13] M. Plapp, Unified derivation of phase-field models for alloy so-lidification from a grand-potential functional, Physical ReviewE 84 (3) (2011) 031601.[14] N. Provatas, K. Elder, Phase-Field Methods in Materials Sci-ence and Engineering, Wiley-VCH Verlag GmbH & Co. KGaA,2010. doi:10.1002/9783527631520.ch1 .[15] W. J. B. A. A. Wheeler, G. B. McFadden, Physical Review A45 (1992) 7424–7439.[16] A. Karma, Phase-field formulation for quantitative modeling ofalloy solidification, Phys. Rev. Lett 87 (2001) 115701.[17] B. Echebarria, R. Folch, A. Karma, M. Plapp, Quantitativephase-field model of alloy solidification, Physical Review E 704