# Comparing mesoscopic models for dendritic growth

Damien Tourret, Laszlo Sturz, Alexandre Viardin, Miha Založnik

► To cite this version:

Damien Tourret, Laszlo Sturz, Alexandre Viardin, Miha Založnik. Comparing mesoscopic models for dendritic growth. IOP Conference Series: Materials Science and Engineering, 2020, 861, pp.012002. 10.1088/1757-899X/861/1/012002. [https://doi.org/10.1088/1757-899X/861/1/012002](https://doi.org/10.1088/1757-899X/861/1/012002). [hal-02978745](https://hal-univ-lorraine.fr/hal-02978745)

HAL Id: hal-02978745

https://hal.univ-lorraine.fr/hal-02978745

Submitted on 23 Dec 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comparing mesoscopic models for dendritic growth

D Tourret¹, L Sturz², A Viardin² and M Založnik³

¹ IMDEA Materials Institute, Getafe, Madrid, Spain

² Access e.V., Aachen, Germany

³ Université de Lorraine, CNRS, IJL, F-54000 Nancy, France

E-mail: damien.turret@imdea.org

**Abstract.** We present a quantitative benchmark of multiscale models for dendritic growth simulations. We focus on approaches based on phase-field, dendritic needle network, and grain envelope dynamics. As a first step, we focus on isothermal growth of an equiaxed grain in a supersaturated liquid in three dimensions. A quantitative phase-field formulation for solidification of a dilute binary alloy is used as the reference benchmark. We study the effect of numerical and modeling parameters in both needle-based and envelope-based approaches, in terms of their capacity to quantitatively reproduce phase-field reference results. In light of this benchmark, we discuss the capabilities and limitations of each approach in quantitatively and efficiently predicting transient and steady states of dendritic growth. We identify parameters that yield a good compromise between accuracy and computational efficiency in both needle-based and envelope-based models. We expect that these results will guide further developments and utilization of these models, and ultimately pave the way to a quantitative bridging of the dendrite tip scale with that of entire experiments and solidification processes.

## 1. Introduction

In metallic alloys obtained by solidification processing, dendritic microstructures are common [1]. From a fundamental standpoint, dendritic growth theory and modeling stands as a challenge to combine phenomena across a wide range of scales, from microscopic capillarity at dendritic tips to macroscopic transport of heat and solute in the melt [2,3]. Due to this multiscale aspect, numerous approaches have emerged over the years that aim at bridging length scales in solidification modeling [4-6]. Each of these approaches operates within a different range of scales, which makes them valuable tools from a technological innovation perspective, since a key hurdle on the way to effective ICME (Integrated Computational Materials Engineering) implementations relies upon our ability to couple models at different length scales [7,8]. Yet, most modeling approaches to dendritic growth have thus far been developed separately, with little effort to compare them and to discuss their capabilities and limitations on a quantitative basis.

In this article, we focus on models of solidification that operate from the microscopic to an intermediate, or mesoscopic, scale between that of the microscopic solid-liquid interface pattern and that of the macroscopic solidification process. At this intermediate scale, interactions between grains occur that determine microstructural features such as grain size, shape, morphology, and crystal texture. Mesoscopic models are more recent than their microscopic and macroscopic counterparts, such that a rigorous quantitative assessment of their relative advantages and limitations remains to be performed. We also focus on moderate solute supersaturation, $\Omega \le 0.25$, which is a common regime for most processes that do not involve rapid solidification.

<figure><img src="image_2.png" /></figure>

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Here, we focus on three prominent models for dendritic crystal growth — namely phase-field, dendritic needle network, and mesoscopic grain envelope models. Using a set of benchmark simulations for isothermal equiaxed growth in a binary alloy at given undercoolings (i.e. at different solute supersaturations), we compare the predicted growth velocities of the dendrite tips in both steady and transient growth regimes. The objective of this comparison is to highlight the advantages and limitations of the different methods, and to suggest how to select model and numerical parameters in order to achieve quantitative predictions while retaining profitable computational efficiency. While this study is still ongoing, the preliminary results presented here already provide useful insight into the capabilities of the models and the choice of parameters.

## 2. Models

### 2.1. Phase-field

The phase-field (PF) method is a powerful tool for the resolution of free-boundary problems. The method can handle moving interfaces of arbitrary morphological complexity. Interfaces and boundaries are not tracked explicitly, but instead different phases or grains are represented by a “phase field” that varies continuously through an interface of finite width W (figure 1(a)).

The free energy of the system is formulated, typically using a Ginzburg-Landau functional [9,10], which describes the energetic contribution of interfaces, transformations, and local state (e.g. temperature, concentration, stresses, etc.) [11-16]. The time evolution of the fields derives directly from the functional derivatives of the constructed free energy functional, using the so-called Cahn-Hilliard equation [17] for a conserved field, such as the solute concentration, or a time-dependent Ginzburg-Landau (also referred to as Allen-Cahn [18]) equation for non-conserved fields. A number of PF models have been proposed, including some that consider a grand potential functional in place of a free energy description [19,20], but the basic recipe typically follows these same common traits for every phase-field model. Hence, what makes the distinction between different PF models is essentially the choice of the number of phase fields (or order parameters) necessary to describe the system and the transformation(s), the mathematical/physical formulation of the functional, the interpolation functions used for different properties between different phases, and the matching of the model parameters to physical parameters relevant to the phenomena one aims to describe.

In the context of solidification, the PF method has been particularly successful, in part due to the presence of well-defined sharp interface problems to be matched by the PF formulations [12,14]. In particular, “thin-interface” asymptotic analyses have permitted using a diffuse interface width W much higher than the actual width of the solid-liquid interface [21,22]. The introduction of a so-called “anti-trapping current” has allowed quantitative simulations of alloys by eliminating spurious numerical solute trapping for wide interfaces [22,23]. In the model considered here, the temperature is assumed homogeneous and constant in time, solid-state diffusion and kinetic undercooling are neglected, a binary alloy is considered within its dilute limit where its phase diagram can be linearized with a constant interface solute partition coefficient k, and solute transport is limited to diffusion in the liquid with a constant diffusivity D. These assumptions result in the following evolution equations [21-24]

$$a_s(\mathbf{n})^2 \frac{\partial \psi}{\partial t} = \nabla \cdot [a_s(\mathbf{n})^2 \nabla \psi] + \varphi \sqrt{2} - \lambda \sqrt{2} (1 - \varphi^2) U + \sum_{m=x,y,z} \left[ \frac{\partial_m (\nabla \psi)^2 a_s(\mathbf{n})}{\partial (\partial_m \psi)} \right] \quad (1)$$

$$[1 + k - (1 - k) \varphi] \frac{\partial U}{\partial t} = D' \nabla \cdot [(1 - \varphi) \nabla U] + [1 + (1 - k)U] \frac{(1 - \varphi^2)}{\sqrt{2}} \frac{\partial \psi}{\partial t} \\ + \nabla \cdot \left\{ [1 + (1 - k)U] \frac{(1 - \varphi^2)}{2} \frac{\partial \psi}{\partial t} \frac{\nabla \psi}{|\nabla \psi|} \right\}, \quad (2)$$

where the dimensionless solute field is

$$U = \frac{1}{1 - k} \left[ \frac{c/c_l^0}{(1 - \varphi)/2 + k(1 + \varphi)/2} - 1 \right], \quad (3)$$
with *c* the solute concentration field, $c_l^0$ the equilibrium solute concentration of a flat interface at the considered temperature *T*₀. In these equations, space is scaled with the diffuse interface width *W* and time with the relaxation time $\tau_0$ at *T*₀. Non-dimensional values of the liquid diffusion coefficient $D' = D\tau_0/W^2 = a_1a_2W/d_0$ and the coupling factor $\lambda = a_1W/d_0$, with $a_1 = 5\sqrt{2}/8$, and $a_2 = 47/75$ from thin-interface asymptotic analysis [21]. We use the standard cubic symmetry of the surface tension with $a_s(n) = (1 - 3\varepsilon_4) + 4\varepsilon_4(n_x^4 + n_y^4 + n_z^4)$, where $n_x, n_y$, and $n_z$ are the components of the unit vector **n** normal to the interface, and $\varepsilon_4$ is the strength of the surface tension anisotropy. The classical phase field $\varphi$, with $\varphi = +1$ in the solid and $\varphi = -1$ in the liquid phase, was replaced by a preconditioned phase field $\psi$, defined by $\varphi = \tanh(\psi/\sqrt{2})$. This change of variable stabilizes the numerical resolution of the equations, thus allowing the use of a coarser grid spacing with negligible loss in accuracy [25]. The model is solved using finite differences on a homogeneous grid of cubic elements and an explicit Euler time stepping scheme. The code is written in cuda (Compute Unified Device Architecture) programming language, which allows its acceleration through multi-threading on Nvidia graphics processing units (GPUs).

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic illustration of the phase field (a), dendritic needle network (b), and envelope (c) models, illustrating characteristic model parameters, namely the PF diffuse interface width *W*, the DNN integration radius $r_i$ and bounding radius $r_{max}$, and the envelope stagnant-film thickness $\delta$.</figcaption></figure>

## 2.2. Dendritic Needle Network

The DNN approach aims at bridging the scale of the dendrite tip radius $\rho$ and the larger scale of transport phenomena, e.g. the solutal diffusion length $l_D = D/V$, with $V$ the crystal growth velocity. The model, rigorously derived in the limit of small Péclet numbers, $P = \rho/(2l_D) \ll 1$, represents a dendritic crystal as a network of thin needles that approximate the primary stems and higher order sidebranches of the grain. In earlier descriptions of the model, in two dimensions, the needle-like branches were treated as infinitely sharp [26]. In later versions of the model, branches were given a parabolic shape, which lead to the same governing equations, while also allowing a three-dimensional extension of the model [27].

The DNN model was shown to quantitatively reproduce primary dendritic spacings measured in directional solidification experiments in different Al-based alloys [27-29]. Implementations have been developed that include fluid flow in the liquid [30,31] and the interaction between nucleation and growth [32,33].

The evolution of the network during solidification is described by the growth dynamics of each needle tip and the generation of new sidebranches. For a binary alloy, typical assumptions of the model include a linear phase diagram with an interface partition coefficient *k*, negligible fluid flow in the liquid phase (at least in the vicinity of the dendrite tips [30]), no solute transport in the solid phase, and negligible release of latent heat during solidification. Under these conditions, for an isothermal domain at a given solute supersaturation $\Omega$, in three dimensions, the needle tip radius, $\rho$, and its velocity, $V$, are calculated by combining the following set of equations [27]:
a solvability condition at the needle tip scale $\approx \rho$

$$ \rho^2 V = 2D d_0 / \sigma, \qquad (4) $$

a solute conservation at intermediate scale much larger than $\rho$ and much smaller than $l_D$

$$ \rho V = D F, \qquad (5) $$

and solute diffusion in the bulk liquid at a scale $\approx l_D$

$$ \partial_t u = D \nabla^2 u, \qquad (6) $$

with

$$ u = \frac{c_l^0 - c}{(1 - k) c_l^0} \qquad (7) $$

the normalized solute concentration, $c_l^0$ the liquid equilibrium concentration and $d_0$ the chemical capillary length, both at the fixed temperature $T_0$, and $\sigma$ the tip selection parameter. The flux intensity factor $F$ is a measure of the solute flux towards the needle in the vicinity of its tip (figure 1(b)), and captures intra- and intergranular solutal interaction. It is defined as

$$ F = \frac{1}{2\pi a} \int_{\Gamma'} \frac{\partial u}{\partial n} d\Gamma' \approx \frac{1}{2\pi a} \left( \int_{\Gamma} \frac{\partial u}{\partial n} d\Gamma + \frac{V}{D} \int_{\Sigma} \frac{\partial u}{\partial x'} d\Sigma \right), \qquad (8) $$

where $a$ is a chosen integration distance behind the tip, typically a few tip radii, and $\Gamma'$ is the resulting surface of the paraboloid up to a distance $a$ from the tip (figure 1(b)). In equation (8), the last equality assumes that the field $u$ is Laplacian in a moving frame of velocity $V$, i.e. $D\nabla^2 u = -V \partial u / \partial x'$, at the scale of the integration domain, with $x'$ the principal growth direction of the branch. This assumption allows calculating the normal solute flux along the surface $\Gamma'$ using the normal solute flux through any surface $\Gamma$ that intersects the parabolic tip at a length $a$ behind the tip, with an enclosed volume $\Sigma$ (figure 1(b)) [27]. In earlier model implementations [26-30], the outer surface $\Gamma$ was set as a square or a cuboid for their ease of implementation onto a square numerical grid. In later works, a circular or spherical shape of the integration domain is preferred [31,32], which facilitates the implementation of grains of arbitrary orientations regardless of the numerical grid.

With $c_\infty$ the alloy nominal solute concentration, the supersaturation $\Omega = (c_l^0 - c_\infty) / [(1-k)c_l^0]$ serves as boundary condition far from the interface, where $u \to \Omega$, while thermodynamic equilibrium neglecting capillary correction is applied along the dendritic branches, with $u = 0$.

We solve the equations numerically on a 3D finite difference mesh with a homogeneous spacing $\Delta x$, using an explicit time stepping scheme with a constant time step $\Delta t$. At each time step, we integrate the value of $F(t)$ using equation (8) with a spherical integration domain of radius $r_i$ centered on the tip of the needle [30], which provides the value of the tip velocity $V(t)$ and radius $\rho(t)$ from equations (4) and (5).

Far behind its tip, we bound the maximum thickness of a branch to a cylinder of radius $r_{\max}$ (figure 1(b)). While it does not affect the paraboloid shape in the tip vicinity, this allows a tip to slow down and stop with $V \to 0$ without numerical instabilities that would otherwise be induced by the constancy of $\rho^2 V$. Sidebranches normal to the growth direction (here assuming $\langle 100 \rangle$ growth directions) are generated at positions $\delta = (N_{br}+\delta_{br}) \times \rho_s$ behind the tip in all four normal directions, when the tip has advanced by a distance $N_{br} \times \rho_s$ since the last branching event. A random fluctuation $\delta_{br}$ is generated for every branching event within a range of amplitude $\Delta_{br}$. The average distance between sidebranches $N_{br}$ and its fluctuation $\Delta_{br}$, both expressed in units of the theoretical steady-state tip radius $\rho_s$, are user input parameters. Both truncating and branching methods are similar to those presented and discussed in previous works [26-32].
## 2.3. Mesoscopic Grain Envelope Model

The Mesoscopic Grain Envelope Model (GEM) [34,35] bridges the scales from the dendrite tip, over the internal grain structure, to the mesoscopic scale of an ensemble of grains that interact by diffusion and convection of solute. The core idea of the GEM is that solutal interactions at the mesoscopic scale can be accurately simulated without a detailed representation of the ramified structure of the solid-liquid interface. Instead, the description of a dendritic grain can be simplified and represented by its envelope, an artificial smooth surface that links the tips of the actively growing dendrite branches (figure 1(c)). Individual branches are not represented by the model, although, by definition, the growth of the envelope depends on the growth speed of the branch tips. The growth of the tips is controlled by the solute flux that they reject into the liquid and is therefore determined by the local solute concentration of the surrounding liquid. In the envelope model a constitutive model of tip kinetics is used that gives the tip growth speed as a function of the solute concentration field adjacent to the envelope.

The model was shown to be capable of quantitatively predicting grain envelope shapes, growth velocities, and internal solid fraction for a single grain growing into an infinite melt [34,36,37] and for transient interactions among several grains [35,38]. The envelope model has been successfully used for simulations of in situ synchrotron X-ray imaging experiments [38-40], for modeling spacing adjustments and growth competition in columnar growth [37,41], and for characterizing grain shapes and growth kinetics under the influence of solutal interactions used in upscaling to macroscopic models [42]. The model was also extended to include fluid flow and solute advection [40,43].

Two key assumptions are made: (i) that the phenomena controlling the growth of a dendrite tip are universal and are therefore valid for tips of any order (primary, secondary, tertiary,...) and (ii) that the characteristic time of the interactions at the small (tip) scale is much smaller than that of the interactions at the large (grain) scale. In this work, we further assume an isothermal domain at a given solute supersaturation, $\Omega$, in three dimensions, and we formulate the model as follows. The normalized concentration at the tip interface is that of thermodynamic equilibrium, neglecting capillary corrections: $u = 0$. A stagnant-film formulation of the Ivantsov solution is written that relates the normalized concentration difference, $u_\delta$, between the liquid at the tip and that at a finite distance $\delta$ from the tip [44] to the growth Péclet number of the tip, $P = \rho V/(2D)$, as

$$u_{\delta} = P \exp(P) \left\{ E_{1}(P) - E_{1}\left(P\left[1 + \frac{2\delta}{\rho}\right]\right) \right\}. \qquad (9)$$

Equation (9) is used to solve for $P$. The concentration $u_\delta$ is obtained from the concentration field in the liquid around the envelope, which is resolved numerically as explained below. The tip selection criterion $\rho^2V = 2Dd_0/\sigma$, is a supplementary relation that gives the tip radius, $\rho = d_0/(\sigma P)$, needed in Equation (9), and the tip speed, $V = 2\sigma P^2 D/d_0$. Dendrite tips are assumed to grow in predefined growth directions. In the present study, the cubic crystal dendrites are given six possible growth directions perpendicular to one another. The normal envelope growth velocity, $V_n$, is calculated from the local tip speed, $V$, by the relation $V_n = V n \cos(\theta)$, where $\theta$ is the angle between the outward normal to the envelope, $n$, and the preferential growth direction that forms the smallest angle with $n$ (figure 1(c)). To track the envelope on a numerical grid we use a phase-field-like interface tracking method [45], combined with a surface reconstruction method for better accuracy [36]. Combining these methods, one can avoid explicit tracking of the envelope, which is defined throughout the whole domain with a level set (iso-surface) of a continuous indicator field. Details are given in references [36,45].

Solute transport at the mesoscopic scale, here by diffusion, is described by volume-averaged equations that are valid in the whole domain, both inside and outside the envelopes. Solidification inside the envelope uses the Gulliver-Scheil assumptions of thermodynamic equilibrium at the solid-liquid interface, negligible diffusion in the solid, and instantaneous diffusion in the liquid [46]. The resulting liquid fraction field, $g_l$, in the interior of the envelopes, was shown to realistically represent the distribution of solid and liquid within an actual dendritic grain envelope [34,36,37]. The liquid
inside the envelope is at the normalized equilibrium concentration, $u_l = 0$. These assumptions lead to the conservation equation for the solute in the liquid phase

$$g_l \partial_t u_l = D \nabla \cdot (g_l \nabla u_l) + \partial_t g_l. \qquad (10)$$

The solution of equation (10) gives the liquid concentration $u_l$ outside the envelope and $g_l$ inside the envelope. Outside the envelope, $g_l = 1$ and equation (10) reduces to a single-phase diffusion equation. Inside the envelope, $u_l = 0$ is the known equilibrium concentration and equation (10) gives the evolution of the liquid fraction inside the envelope. The concentration of the solid growing inside the envelope at constant $\Omega$ is $u_s = 1$.

The supersaturation across the stagnant film, $u_\delta$, needed to calculate the envelope growth velocity (equation (9) and additional relations), is obtained from the concentration field in the liquid around the envelope, which is fully resolved numerically (equation (10)). This means that the locally valid Ivantsov solution for diffusion around a dendrite tip is matched to the solution of the mesoscopically valid solute transport equation at a distance $\delta$ from the envelope (figure 1(c)).

The solute conservation equation (10) and the equation of the evolution of the indicator field used for envelope tracking are solved by a finite volume solver with the OpenFOAM® toolbox. Second and higher order discretization schemes are used for the spatial operators. Implicit Euler time-stepping is used for the solute conservation equations and explicit time-stepping for the envelope indicator field equation. Adaptive time-stepping is employed and is useful during fast transients, such as the initial transient after nucleation. During other growth stages the time step is essentially limited by the characteristic diffusion time, $\Delta x^2/D$, at the scale of the mesh cell, $\Delta x$. In the present work uniform meshes with cubic cells are used. The mesh size needs to be adapted for an accurate description of the diffusion layer adjacent to the envelope, i.e. $\Delta x \lesssim 0.2D/V$, where $V$ is a typical primary tip speed. A special search and interpolation algorithm is designed to provide the concentration at the distance $\delta$ from the envelope, needed for the calculation of the envelope velocity [36].

### 3. Benchmarks

Here, we compare the predictions of the three methods described above with a similar benchmark test case for undercooled isothermal equiaxed growth of a binary alloy. The phase field simulation uses well-converged numerical parameters, i.e. small enough grid spacing and diffuse interface width $W$, such that the PF results will be used as a reference for the discussion of the effect and choice of parameters in the other two methods.

We set a solid seed (about as small as the method allows) in the corner of a cubic domain initialized at a given solute supersaturation $\Omega$. Thus, we simulate the growth of 1/8 of an equiaxed grain, using no-flux (mirror symmetry) Neumann boundary conditions on every face of the cube. We choose a value of the domain size $L$ (edge length of the cube) that is large enough to reach (or closely approach) a steady growth velocity, but also small enough to make simulations computationally achievable within a few days (less than two weeks for the longest) using quantitative phase-field using the current single-GPU-parallelized implementation. Hence, we selected $L/l_D \approx 2, 4, 6, 8$, and 10 for $\Omega = 0.05, 0.10, 0.15, 0.20$, and 0.25, respectively.

In the PF simulations, we use a small interface energy anisotropy $\epsilon_4 = 0.007$ and a solute partition coefficient $k = 0.1$. In the DNN and GEM simulation we consider a tip selection constant $\sigma = 0.0264$. This value was estimated from preliminary PF simulations of isothermal growth for $0.05 \le \Omega \le 0.25$ in a domain with $512^3$ grid points with tip velocities and radii measured as soon as the diffusion field gets affected at the side of the domain opposite to the dendrite. This value agrees with results from joint directional solidification simulation with similar $\epsilon_4$ and $k$. We noted a typical uncertainty of the order of 10 to 15% depending upon the procedure chosen to extract the dendrite tip radius [47]. Here, we follow the method described in reference [48] (figure 5 therein), i.e. by 2D cross-section least-square fitting to a second order parabola over a fitting range that spans a distance of one tip radius behind the tip. A deeper analysis in terms of selection parameter is underway, such that the value used here only provides a reasonable estimation of $\sigma$ to use consistently in DNN and GEM simulations.
The objective is to compare predicted growth velocities, not only in the steady-state, which should closely match the Ivantsov solution, but also and most importantly the initial transient, as the velocity tends towards the free dendrite steady-state, and the final transient, as the velocity decreases to zero when the dendrite approaches the no-flux boundary wall.

For a given solute supersaturation $\Omega$, the theoretical Péclet number $P_s = \rho_s V_s / (2D)$ at steady-state (hence marked with a subscript “s”) is given by the Ivantsov solution [49]. It directly gives the ratio between the steady-state diffusion length $D/V_s$ and dendrite tip radius $\rho_s$. In this paper, most results are presented scaled with respect to the steady-state tip radius $\rho_s$ and velocity $V_s$. For each value of $\Omega$ considered here, table 1 summarizes the correspondence between these theoretical steady-state values and material parameters in physical units, such as the liquid solute diffusivity $D$ and the interface capillarity length $d_0$. Within the selected supersaturation range, there is an order of magnitude of variation in the theoretical steady-state values of the Péclet and tip radius, and two orders of magnitude of variation of the tip velocity.

In the phase-field simulations, the finite difference grid element size is set equal to the diffuse interface width with $\Delta x = W \approx \rho_s / 10$, which was found to yield well-resolved converged results. The initial spherical seed has a radius $r_0 = 2.5 \Delta x$, with a preconditioned phase field $\psi = r - r_0$ with $r = (x^2 + y^2 + z^2)^{1/2}$ and a solute field $U = -\Omega$ everywhere, which is equivalent to setting $c = c_\infty$ (i.e. $u = \Omega$) in the liquid and $c = k c_\infty$ in the solid.

In the DNN simulations, we initialize a solid seed with $u = 0$ in the corner of the domain and set $u = \Omega$ in the liquid. The seed consists of three quarters of paraboloids along the three normal grid directions with initial tip velocity $V = 0$, radius $\rho = \rho_s$, and length $L_0$ larger than $r_i$ and $r_{\max}$. In particular, we study the effect of mesh spacing, $\Delta x$, and truncation radius, $r_{\max}$. During the final transient, the integration sphere may be truncated if it exceeds the limits of the domain.

In the GEM simulations, we initialize a spherical envelope seed with $u = 0$ and $g_l = 1 - \Omega$ in the corner of the domain and set $u = \Omega$ in the liquid. The size of the envelope seed is given by the smallest radius of curvature that can be accurately represented by the envelope tracking method on a mesh with spacing $\Delta x$ [36,45]: $r_0 = 4.46 \Delta x$. The mesh spacing and time step needed for well-resolved converged results were determined by systematic convergence studies. Remaining mesh dependence, due to the initial seed size, which can be reduced on finer grids, is perceptible only during the initial growth transient and is discussed in section 3.2.

<table><caption>**Table 1.** Theoretical steady-state Péclet number $P_s$, tip radius $\rho_s$, and velocity $V_s$, for different $\Omega$.</caption><thead><tr><th>Ω</th><th>P<sub>s</sub> = ρ<sub>s</sub>V<sub>s</sub> / (2D)</th><th>ρ<sub>s</sub> / d<sub>0</sub> = 1 / (σP<sub>s</sub>)</th><th>V<sub>s</sub>d<sub>0</sub> / D = 2σP<sub>s</sub><sup>2</sup></th></tr></thead><tbody><tr><td>0.05</td><td>0.0131</td><td>2895</td><td>9.03·10<sup>-6</sup></td></tr><tr><td>0.10</td><td>0.0341</td><td>1111</td><td>6.14·10<sup>-5</sup></td></tr><tr><td>0.15</td><td>0.0624</td><td>607</td><td>2.06·10<sup>-4</sup></td></tr><tr><td>0.20</td><td>0.0988</td><td>383</td><td>5.15·10<sup>-4</sup></td></tr><tr><td>0.25</td><td>0.1446</td><td>261</td><td>1.10·10<sup>-3</sup></td></tr></tbody></table>

## 4. Results and discussions

The results presented here were obtained using model and numerical parameters that gave a good fit to the time evolution of the growth velocities obtained from PF simulations, while conserving a substantially lower computational cost. Since we used a wide range of computers (e.g. CPU/GPU) and methods (e.g. explicit/implicit time stepping), we did not perform a systematic performance analysis and focused on the potential of each method to yield quantitative predictions. Yet, in the following subsection, we provide values used for grid spacing and time steps, which provide a consistent picture of the relative computational workload to each method.
## 4.1. Steady-state growth

First, we compare steady-state growth predictions by all models for the range of solute supersaturation $0.05 \le \Omega \le 0.25$. Steady-state growth Péclet numbers predicted by PF, DNN, and GEM, are compared to the theoretical Ivantsov solution in figure 2. Values reported in figure 2 are averaged in time during steady-state growth, in order to integrate numerical oscillations. Such oscillation may occur as a dendrite tip crosses successive grid points and are due to the coarse grids used in DNN and GEM. The fluctuations in tip velocity and radius (see e.g. figure 4 discussed later) lead to a negligible fluctuation in Péclet number (within the symbol size used in figure 2). In summary, all models are capable of solving numerically the Ivantsov problem, and they thus appropriately predict the steady growth state.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Predicted steady-state Péclet number for different imposed supersaturation $\Omega$ compared to the analytical Ivantsov solution [49].</figcaption></figure>

Model parameters used for these simulations are summarized in table 2. In order to provide a consistent illustration of the computational cost for all three methods, we also summarize the numerical parameters, namely the grid element size and time step, for each $\Omega$ in table 3.

<table><caption>**Table 2.** Model parameters used with the three methods.</caption><thead><tr><th rowspan="2">Ω</th><th>PF</th><th colspan="2">DNN</th><th>GEM</th></tr><tr><th>W/d<sub>0</sub></th><th>r<sub>i</sub>/ρ<sub>s</sub></th><th>r<sub>max</sub>/ρ<sub>s</sub></th><th>δ/(D/V<sub>s</sub>)</th></tr></thead><tbody><tr><td>0.05</td><td>289</td><td>2.6</td><td>3, 10</td><td>0.06</td></tr><tr><td>0.10</td><td>110</td><td>2.4</td><td>3</td><td>0.20</td></tr><tr><td>0.15</td><td>60</td><td>1.6</td><td>3</td><td>0.40</td></tr><tr><td>0.20</td><td>38</td><td>1.9</td><td>3</td><td>0.60</td></tr><tr><td>0.25</td><td>26</td><td>0.6</td><td>2</td><td>0.80</td></tr></tbody></table>

<table><caption>**Table 3.** Normalized mesh element size and time step used with the three methods.</caption><thead><tr><th rowspan="2">Ω</th><th colspan="3">Δx/(D/V<sub>s</sub>)</th><th colspan="3">Δt/(D/V<sub>s</sub><sup>2</sup>)</th></tr><tr><th>PF</th><th>DNN</th><th>GEM</th><th>PF</th><th>DNN</th><th>GEM</th></tr></thead><tbody><tr><td>0.05</td><td>0.0026</td><td>0.0065</td><td>0.02</td><td>10<sup>-6</sup></td><td>7·10<sup>-6</sup></td><td>10<sup>-4</sup></td></tr><tr><td>0.10</td><td>0.0068</td><td>0.017</td><td>0.04</td><td>7·10<sup>-6</sup></td><td>4·10<sup>-5</sup></td><td>4·10<sup>-4</sup></td></tr><tr><td>0.15</td><td>0.0125</td><td>0.031</td><td>0.06</td><td>2·10<sup>-5</sup></td><td>10<sup>-4</sup></td><td>9·10<sup>-4</sup></td></tr><tr><td>0.20</td><td>0.0198</td><td>0.049</td><td>0.08</td><td>6·10<sup>-5</sup></td><td>4·10<sup>-4</sup></td><td>10<sup>-3</sup></td></tr><tr><td>0.25</td><td>0.0289</td><td>0.072</td><td>0.10</td><td>10<sup>-4</sup></td><td>8·10<sup>-4</sup></td><td>2·10<sup>-3</sup></td></tr></tbody></table>
The shape of the predicted grain when the length of the primary trunk reaches one half of the domain size is shown in figure 3 for all three models for $\Omega = 0.05$ (top) and $\Omega = 0.25$ (bottom). For the sake of illustration, the DNN microstructure is shown for a simulation performed without sidebranches and no thickness bounding (i.e. with $r_{max} = \infty$) for $\Omega = 0.05$ and with sidebranches and $r_{max} = 2\rho_s$ for $\Omega = 0.25$. In the latter case, the effect of a finite $r_{max}$ becomes apparent from the cylindrical shape of the branches. While at this stage the velocity has already started decreasing due to the effect of boundaries, the shape of the microstructures is barely distinguishable from that at steady-state, for all three models.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Grain shape (gray) and iso-concentration lines, from $u = 0$ (red) to $u = 0.9\Omega$ (blue) with steps of $0.1\Omega$, along the ($x = 0$), ($y = 0$), and ($z = 0$) planes when the length of the primary dendritic trunk reaches one half of the simulation domain for $\Omega = 0.05$ (a) and 0.25 (b) simulated using phase-field (left), DNN (middle), and Envelope (right) methods. The total size of the domain, illustrated by the visible bottom left corner of each simulation, is $2 D/V_s \approx 221400 d_0$ for $\Omega = 0.05$ and $10 D/V_s \approx 9060 d_0$ for $\Omega = 0.25$.</figcaption></figure>

## 4.2. Transient growth

Next, we compare the predicted transient evolution of grain growth velocities. We focus on $\Omega = 0.05$ and $\Omega = 0.25$, which are illustrative of results obtained at intermediate $0.05 < \Omega < 0.25$. The evolution of the velocity from initial to final transient appears in figure 4 for $\Omega = 0.05$ (a) and $\Omega = 0.25$ (b). The velocity is shown as a function of the tip location normalized by the domain length. This way, we filter out the possible different onset times for the final decrease of the velocity — such differences arise when the steady-state is maintained at a different velocity. The right-hand-side panels, showing the evolution of the primary trunk length versus time, illustrate that the oscillations appearing in the left-hand-side velocity plot have a negligible effect on the length of the dendrite, once integrated over time. These results show that, for appropriately chosen sets of parameters, both DNN and GEM models can reproduce PF results with reasonable accuracy.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Evolution of the primary dendrite tip velocity from the initial to the final transient state (left) and primary trunk length versus time (right) for Ω = 0.05 (a) and Ω = 0.25 (b). In order to allow comparing models all together at different supersaturations, time and space are scaled with respect to the theoretical steady-state velocity *V*<sub>s</sub> and tip radius ρ<sub>s</sub>. The velocity evolution is plotted as a function of the tip location normalized by the domain length. Note that left-hand-side velocity plots have been smoothed to reduce numerical oscillations and improve readability (right-hand-side length plots are raw).</figcaption></figure>

DNN simulations were performed with Δx = 0.25ρs, rᵢ = 2.6ρs, rmax = 3ρs for Ω = 0.05, and with
Δx = 0.25ρs, rᵢ = 0.6ρs, rmax = 2ρs for Ω = 0.25 (see table 2). We chose a resolution providing a
compromise between accuracy and simulation time. The chosen grid size is smaller than the typical
value used in DNN simulations, i.e. Δx ≈ ρs or higher. Other parameters are close to typical parameters
used in the DNN model, i.e. fulfilling rᵢ ≪ lD and rmax chosen as to not bound the needle thickness
within the flux integration domain. Interestingly, we notice very little influence of sidebranching on
the velocity of the primary tip, both in steady and transient states. It is also worth noting that the
integration radius rᵢ can be taken lower than the tip radius, as long as the integration domain contains a
sufficient number of integration points.
GEM simulations were performed with a stagnant film thickness of $\delta = 0.06 D/V_s$ and $0.80 D/V_s$, and with a mesh spacing of $\Delta x = 0.02 D/V_s$ and $0.10 D/V_s$, for $\Omega = 0.05$ and $0.25$, respectively. The $\delta$ and the mesh used at high $\Omega$ are close to those already recommended by Souhar et al. [36] for steady-state equiaxed growth. At low $\Omega$, significantly smaller $\delta$ were needed to accurately match both the initial and the final transient of the tip evolution. The same conclusion was recently made by Olmedilla et al. [38] when comparing the GEM to experiments on transient equiaxed growth at comparably low undercoolings. The mesh spacing is dictated by the diffusion length at high $\Omega$ ($\Delta x \approx 0.1 D/V_s$ here), however at low $\Omega$ the small stagnant-film thickness, $\delta$, becomes the limiting factor. For an accurate calculation and interpolation of $u_\delta$ the mesh spacing should be $\Delta x < 0.5\delta$.

### 4.3. Effect of DNN parameters

Figure 5 shows the tip velocity evolution predicted by DNN simulations for different combinations of grid spacing $\Delta x$ and needle truncation radius $r_{\text{max}}$, for a flux integration radius $r_i = 2.6 \sim 2.7 \rho_s$ and $\Omega = 0.05$.

Decreasing the grid spacing $\Delta x$ for identical values of $r_{\text{max}}$ and $r_i$ results mainly in (i) less fluctuation of the velocity and (ii) higher velocities, in better agreement with PF results. On the other hand, the effect of increasing the radius of the stem truncation $r_{\text{max}}$ for identical values of $\Delta x$ and $r_i$ is to decrease the steady-state tip velocity.

The effect of the initial length of the seed branches needs to be discussed. The initial length $L_0$ of the branches differs among the simulations because the simulations are initialized with paraboloids of length $L_0 = r_{\text{max}} + \Delta x$. The difference in initial length most notably affects the growth kinetics during the initial and final transient stages. All four simulations nonetheless stabilize close to the theoretical steady-state velocity given by the Ivantsov solution. This plateau region, however, lasts for a shorter time for the simulations with larger $r_{\text{max}}$ and hence longer $L_0$. At this low $\Omega$, in order for PF simulations to be achievable within a few days, we limited the domain size to $2l_D$. For this reason, the amount of time spent in steady-state is limited, and for high $L_0$ values, the final transient sets in before steady-state is fully maintained. In the last stage of the final transient, a velocity plateau may be observed (not represented here but visible by interrupted plots) due to a numerical restriction set to the tip velocity when the distance to the wall becomes lower than $r_i$. This limitation might be accommodated by either (i) truncating the integration domain, (ii) integrating boundary conditions within the integration of $F$, or (iii) using an adaptive $r_i$ that decreases when approaching a boundary.

<figure><img src="image_7.png" /></figure>

**Figure 5.** Evolution of the primary dendrite tip velocity for $\Omega = 0.05$ from DNN simulations, showing the effect of grid spacing $\Delta x$ and stem truncation radius $r_{\text{max}}$. In all runs, the radius of the integration sphere $r_i/\rho_s$ was kept to 2.6 or 2.7. For readability, symbols show raw data whereas lines show smoothed data (filtering numerical oscillations as the tip progresses between successive grid points).

### 4.4. Effect of envelope parameters

The main model parameter of the GEM is the stagnant film thickness, $\delta$. Table 2 summarizes the $\delta$ that gave the best match to the tip evolution predicted by phase field. The best $\delta$ clearly increases with $\Omega$. The dependence of the transient growth velocity upon the stagnant film thickness $\delta$ at $\Omega = 0.05$ is shown in figure 6(a) for a fixed grid spacing $\Delta x/(D/V_s) = 0.02$. The most striking observation is the
effect on the steady-state velocity (already pointed out by Souhar *et al.* [36]), which affects the whole evolution. The de- and accelerations during the initial and final transient do not significantly depend on $\delta$.

Figure 6(b) shows the influence of the mesh spacing on the transient tip growth at $\Omega = 0.05$. The influence of the mesh spacing is twofold. In addition to its direct impact on the solution of the PDEs, $\Delta x$ also limits the smallest radius of the initial nucleus, $r_0$, and thus the initial conditions (see section 2.4). Typically, $r_0 \ge 4.46 \Delta x$. The direct impact of $\Delta x$, excluding that of $r_0$, is shown by four simulations in which $r_0$ was kept fixed and only the mesh spacing was varied from 0.013 $D/V_s$ to 0.04 $D/V_s$. The initial radius was $r_0 = 0.18 D/V_s$ which is the smallest $r_0$ that can be resolved on the coarsest of these four grids. The differences between the four solutions are strictly limited to the initial transient, and perfectly overlap at larger arm lengths. The reason for this is that the thin diffusion layer around the grain during the initial stages of growth requires a fine mesh resolution to be accurately resolved. Once the tip decelerates and the diffusion layer grows wider, coarse meshes are sufficient.

However, the most important influence related to the mesh resolution is actually that of the initial nucleus size. This is shown through a comparison of the evolution of tip velocity for $r_0 = 0.18 D/V_s$ and $0.09 D/V_s$ with the same $\Delta x = 0.02 D/V_s$ (open-red versus full-black square symbols in figure 6(b)). We note a much slower deceleration during the initial transient at smaller $r_0$. Furthermore, the steady growth stage is replaced by a smooth transition to the final transient, which is not significantly affected.

In summary, accurate simulation of transient equiaxed growth requires a stagnant film thickness that roughly scales as $\delta \sim \Omega D/V_s$. This does not contradict, but goes beyond prior results [36] that established a recommendation for steady-state growth: $\delta \sim D/V_s$. If compared, the new recommendation is slightly less accurate in steady-state (with a difference of < 10%) but is more accurate during transients.

The mesh spacing and the initial nucleus radius have a significant influence on the results only during the initial growth transient, but not during steady-state and the final transient. Note that because of the relatively small domain size ($L = 2 D/V_s$) in the case presented here ($\Omega = 0.05$), the impact of the initial transient on the whole evolution is prominent. This is not the case at higher $\Omega$, where larger domains are used and the initial transient is a smaller portion of the growth. Generally, $\Delta x \sim 0.1 D/V_s$ is sufficiently fine for steady-state and for the final growth transient, but should be refined if good accuracy during initial transients is required. Furthermore, $\Delta x < 0.5\delta$ is an additional constraint that is important at low $\Omega$, where smaller $\delta/(D/V_s)$ must be used.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Evolution of the primary dendrite tip velocity for $\Omega = 0.05$ from GEM simulations, showing the effect of the model and numerical parameters. (a) The effect of the stagnant film thickness, $\delta$; all simulations with a mesh $\Delta x/(D/V_s) = 0.02$. (b) The effect of the mesh spacing $\Delta x$, and of the size of the initial nucleus, $r_0$; all simulations with $\delta/(D/V_s) = 0.08$.</figcaption></figure>
**5. Summary and perspectives**

We have compared the results of phase-field, dendritic needle network, and grain envelope models for a simulation of the equiaxed growth of an isothermally undercooled binary alloy. Using well-converged PF results as the benchmark reference, we discussed the effect and choice of model and numerical parameters in the DNN and GEM methods upon their accurate prediction of steady and transient dendrite growth velocities.

As expected, all models can appropriately predict steady-state growth conditions prescribed by microscopic solvability and Ivantsov solutions. Moreover, in the transient regimes, both GEM and DNN can accurately describe growth velocities, if appropriate model and numerical parameters are selected. The final transients (solutal interactions between grains in practical applications) are typically accurate with standard parameters. The early-stage growth transient, on the other hand, requires a fine tuning of model parameters in order to achieve quantitative predictions, such that the initial seed size and mesh spacing need to be finer than typical values recommended in previous works, in both DNN and GEM approaches.

While the present study provides already some new insight and guidance for the use of these relatively new modeling techniques, these results remain preliminary, and the study is ongoing. Work in progress on the current benchmarks include: (i) the extraction of scaling laws for initial and final transient growth stages, (ii) further systematic exploration of model parameters (e.g. $r_i$ for DNN), (iii) comparison of envelope shapes in equiaxed growth, which is important in the scope of grain interactions, and (iv) exploring the advantages and limitations of models at higher supersaturations. Perspectives for future benchmarks include: (i) benchmarks for directional solidification, in particular in the scope of spacing selection, and (ii) configurations with convection in the liquid phase.

From the results of these studies, we expect that new numerical techniques will improve the predictive capabilities of the mesoscopic models. Such techniques could include, for instance, the introduction of adaptive parameters to ensure better treatment of transient regimes, e.g. using adaptive $r_i$ in DNN or $\delta$ in GEM that could depend upon the growth acceleration/deceleration. We also expect that the present work, and those to follow, will provide useful guidance to potential users and developers of these promising tools for bridging length and time scales in the modeling of dendritic solidification.

**Acknowledgements**

D.T. acknowledges support from the European Union's Horizon 2020 research and innovation programme under a Marie Skłodowska-Curie Individual Fellowship (Grant Agreement 842795). L.S. and A.V. acknowledge support by the German Ministry of Economy via the Deutsche Zentrum für Luft- und Raumfahrt, DLR (Contract 50WM1743). M.Z. acknowledges support by the French State through the program "Investment in the future" operated by the National Research Agency (ANR) and referenced by ANR-11 LABX-0008-01 (LabEx DAMAS).

**References**

[1] Flemings MC 1974 *Solidification Processing* (New York: McGraw Hill)

[2] Langer JS 1980 *Rev. Mod. Phys.* **52** 1

[3] Trivedi R and Kurz W 1994 *Acta Metall. Mater.* **42** 15

[4] Boettinger WJ, Coriell SR, Greer A, Karma A, Kurz W, Rappaz M and Trivedi R 2000 *Acta Mater.* **48** 43

[5] Asta M, Beckermann C, Karma A, Kurz W, Napolitano R, Plapp M, Purdy G, Rappaz M and Trivedi R 2009 *Acta Mater.* **57** 941

[6] Karma A and Tourret D 2016 *Curr. Opin. Solid St. Mater. Sci.* **20** 25

[7] Allison J 2011 *JOM* **63** 15

[8] Panchal JH, Kalidindi SR and McDowell DL 2013 *Comput. Aided Design* **45** 4

[9] Ginzburg V and Landau L 1950 *Soviet Phys. JETP* **20** 1064

[10] Halperin B, Hohenberg P and Ma S 1974 *Phys. Rev. B* **10** 139
[11] Chen LQ 2002 *Annu. Rev. Mater. Res.* **32** 113

[12] Boettinger WJ, Warren JA, Beckermann C and Karma A 2002 *Annu. Rev. Mater. Res.* **32** 163

[13] Moelans N, Blanpain B and Wollants P 2008 *Calphad* **32** 268

[14] Steinbach I 2009 *Model. Simul. Mater. Sci. Eng.* **17** 073001

[15] Plapp M 2016 *J. Indian Inst. Sci.* **96** 179

[16] Tonks MR and Aagesen LK 2019 *Annu. Rev. Mater. Res.* **49** 79

[17] Cahn JW and Hilliard JE 1958 *J. Chem. Phys.* **28** 258

[18] Allen SM and Cahn JW 1979 *Acta Metall.* **27** 1085

[19] Plapp M 2011 *Phys. Rev. E* **84** 031601

[20] Choudhury A and Nestler B 2012 *Phys. Rev. E* **85** 021602

[21] Karma A and Rappel WJ 1998 *Phys. Rev. E* **57** 4323

[22] Echebarria B, Folch R, Karma A and Plapp M 2004 *Phys. Rev. E* **70** 061604

[23] Karma A 2001 *Phys. Rev. Lett.* **87** 115701

[24] Tourret D and Karma A 2015 *Acta Mater.* **82** 64

[25] Glasner K 2001 *J. Comput. Phys.* **174** 695

[26] Tourret D and Karma A 2013 *Acta Mater.* **61** 6474

[27] Tourret D and Karma A 2016 *Acta Mater.* **120** 240

[28] Tourret D, Clarke AJ, Imhoff SD, Gibbs PJ, Gibbs JW and Karma A 2015 *JOM* **67** 1776

[29] Tourret D, Karma A, Clarke AJ, Gibbs PJ and Imhoff SD 2015 *IOP Conf. Ser.: Mater. Sci. Eng.* **84** 012082

[30] Tourret D, Francois MM and Clarke AJ 2019 *Comput. Mater. Sci.* **162** 206

[31] Isensee T and Tourret D 2020 *IOP Conf. Ser.: Mater. Sci. Eng. current issue*

[32] Sturz L and Theofilatos A 2016 *Acta Mater.* **117** 356

[33] Geslin P and Karma A 2016 *TMS 2016 Annual Meeting, Supplemental Proceedings, Symposium Frontiers in Solidification* 29

[34] Steinbach I, Beckermann C, Kauerauf B, Li Q and Guo J 1999 *Acta Mater.* **47** 971

[35] Steinbach I, Diepers HJ and Beckermann C 2005 *J. Cryst. Growth* **275** 624

[36] Souhar Y, De Felice VF, Beckermann C, Combeau H and Založnik M 2016 *Comput. Mater. Sci.* **112** 304

[37] Viardin A, Založnik M, Souhar Y, Apel M and Combeau H 2017 *Acta Mater.* **122** 386

[38] Olmedilla A, Založnik M and Combeau H 2019 *Acta Mater.* **173** 249

[39] Delaleau P, Beckermann C, Mathiesen RH and Arnberg L 2010 *ISIJ Int.* **50** 1886

[40] Olmedilla A, Založnik M, Cisternas Fernández M, Viardin A and Combeau H 2019 *IOP Conf. Ser.: Mater. Sci. Eng.* **529** 012040

[41] Li Y, Olmedilla A, Založnik M, Zollinger J, Dembinski L and Mathieu A 2019 *Proceedings 5th International Conference on Advances in Solidification Processes (ICASP-5)*

[42] Torabi Rad M, Založnik M, Combeau H and Beckermann C 2019 *Materialia* **5** 100231

[43] Delaleau P 2011 *Mesoscale modeling of dendritic growth during directional solidification of aluminium alloy* Ph.D. thesis NTNU, Trondheim, Norway

[44] Cantor B and Vogel A 1977 *J. Cryst. Growth* **41** 109

[45] Sun Y and Beckermann C 2007 *J. Comput. Phys.* **220** 626

[46] Scheil E 1942 *Z. Metallkd.* **34** 70

[47] Karma A, Lee YH and Plapp M 2000 *Phys. Rev. E* **61** 3996

[48] Clarke A, Tourret D, Song Y, Imhoff S, Gibbs P, Gibbs J, Fezzaa K and Karma A 2017 *Acta Mater.* **129** 203

[49] Ivantsov GP 1947 *Dokl. Akad. Nauk. SSSR* **58** 567