# A Stable Dendritic Growth with Forced Convection: A Test of Theory Using Enthalpy-Based Modeling Methods

A. Kao, L. V Toropova, I. Krastins, Gilles Demange, D. V Alexandrov, P. K Galenko

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

## To cite this version:

A. Kao, L. V Toropova, I. Krastins, Gilles Demange, D. V Alexandrov, et al., A Stable Dendritic Growth with Forced Convection: A Test of Theory Using Enthalpy-Based Modeling Methods. JOM Journal of the Minerals, Metals and Materials Society, 2020, 72 (9), pp.3123-3131. 10.1007/s11837-020-04292-4 . hal-02943168

HAL Id: hal-02943168

https://normandie-univ.hal.science/hal-02943168

Submitted on 20 Jul 2021

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Stable Dendritic Growth with Forced Convection: A Test of Theory Using Enthalpy-Based Modeling Methods

A. KAO<sup>1,7</sup> L.V. TOROPOVA,<sup>2,8</sup> I. KRASTINS,<sup>3</sup> G. DEMANGE,<sup>4</sup>
D.V. ALEXANDROV,<sup>5</sup> and P.K. GALENKO<sup>6</sup>

1.—Centre for Numerical Modelling and Process Analysis, University of Greenwich, Old Royal Naval College, Park Row, London SE10 9LS, UK. 2.—Laboratory of Mathematical Modeling of Physical and Chemical Processes in Multiphase Media, Department of Theoretical and Mathematical Physics, Ural Federal University, Lenin ave., 51, Ekaterinburg, Russian Federation 620000. 3.—Institute of Physics, University of Latvia, Miera iela 32, Salaspils LV-2169, Latvia. 4.—GPM, CNRS-UMR 6634, University of Rouen Normandy, 76801 Saint Étienne Du Rouvray, France. 5.—Laboratory of Multi-Scale Mathematical Modeling, Department of Theoretical and Mathematical Physics, Ural Federal University, Lenin ave., 51, Ekaterinburg, Russian Federation 620000. 6.—Physikalisch-Astronomische Fakultät, Friedrich-Schiller-Universität-Jena, 07743 Jena, Germany. 7.—e-mail: a.kao@gre.ac.uk. 8.—e-mail: l.v.toropova@urfu.ru

The theory of stable dendritic growth within a forced convective flow field is tested against the enthalpy method for a single-component nickel melt. The growth rate of dendritic tips and their tip diameter are plotted as functions of the melt undercooling using the theoretical model (stability criterion and undercooling balance condition) and computer simulations. The theory and computations are in good agreement for a broad range of fluid velocities. In addition, the dendrite tip diameter decreases, and its tip velocity increases with increasing fluid velocity.

## INTRODUCTION

It is well known that the growth of dendritic crystals takes place in many areas of modern science ranging from materials physics, geophysics and atmosphere physics to the chemical industry, biophysics and life science<sup>1–6</sup>. As this takes place, the growth mechanisms of dendrites, their shape and interaction determine the characteristics of the internal microstructure of the crystallized substance<sup>7–9</sup>. These mechanisms, in turn, depend on heat and mass transfer processes complicated by hydrodynamic and convective fluid flows, the presence of dissolved impurities and various crystal growth symmetries. To determine the stable dendritic growth mode, as well as to establish the boundaries of morphologic transitions of the internal structure in solidified materials, it is necessary to independently determine the growth rate V of the dendrite tip and its diameter ρ depending on the

melt undercooling $\Delta T$. This problem can be solved using the microscopic solvability theory together with the sharp interface model, which lead to two transcendental equations for V and ρ as functions of $\Delta T$ and other physical parameters of dendritic growth. Such a theoretical approach has been recently tested against experimental data and computations in a series of works<sup>10–13</sup> in the absence of a forced convection. The present study compares the theory of stable dendritic growth in the presence of a forced convection with computer simulations by the enthalpy method.

Our article is organized as follows. Section 2 summarizes the main outcomes following from the solvability theory and the sharp interface model keeping in mind the fourfold crystalline symmetry and a forced convective flow. Here, we present the thermally controlled model of anisotropic dendritic growth as well as its final analytical solution, which consists of two transcendental equations for V and ρ. Section 3 is concerned with the background of the enthalpy method and our main results following

<footnote>(Received April 24, 2020; accepted July 13, 2020;</footnote>
from simulations of dendritic growth. The main outcomes of our analysis and future directions are discussed in Sect. 4.

## SHARP INTERFACE MODEL AND THE SOLVABILITY THEORY

### Governing Equations

Consider the steady-state growth of a two-dimensional thermal dendritic crystal along the spatial axis $z$ in the presence of a forced convective flow coming from the opposite direction¹⁴ (Fig. 1).

The heat balance condition at the dendritic interface can be written as

$$T_Q(\mathbf{v} \cdot \mathbf{n}) = D_T(\nabla T_s - \nabla T_l) \cdot \mathbf{n}, \quad (1)$$

where $T_Q$ is the hypercooling, $\mathbf{v} \cdot \mathbf{n}$ is the growth velocity, $D_T$ is the thermal diffusivity, $T_s$ and $T_l$ are the temperatures in solid and liquid phases, respectively.

The heat transport equations in the liquid and solid phases take the form

$$\frac{\partial T_l}{\partial t} + (\mathbf{w} \cdot \nabla) T_l = D_T \nabla^2 T_l, \quad \frac{\partial T_s}{\partial t} = D_T \nabla^2 T_s, \quad (2)$$

where $\mathbf{w}$ is the fluid velocity, $t$ is the time variable, and $\nabla$ is the differential nabla operator. Note that we consider the case of equal thermal diffusivities in both the phases, because this hypothesis does not change the selection criterion (it can change the selection constant only).

To describe the hydrodynamic flows, we use the linearized Oseen model for a viscous flow¹⁵⁻¹⁷

$$U \frac{\partial \mathbf{w}}{\partial z} = - \frac{1}{\rho_l} \nabla p + \nu \nabla^2 \mathbf{w}, \quad \nabla \cdot \mathbf{w} = 0. \quad (3)$$

Here, $U$ is the fluid velocity far from the dendritic surface (see Fig. 1), $p$ is the pressure, $\rho_l$ is the density of liquid, and $\nu$ is the kinematic viscosity.

The Gibbs-Thomson equation at the solid-liquid boundary holds

$$T_{\text{int}} = T_l = T_s = T_m - T_Q d(\theta, \phi) \mathcal{K} - \tilde{\beta}(\theta, \phi) v_n, \quad (4)$$

<figure><img src="image_3.png" /><figcaption>Fig. 1. A scheme of growing dendritic crystal in a convective flow.</figcaption></figure>

where $T_{int}$ is the phase transition temperature at the dendrite interface, $T_m$ is the melting temperature for the pure system, and $\mathcal{K}$ is the interface curvature, $d(\theta, \phi)$, and $\tilde{\beta}(\theta, \phi)$ are the capillary length and the function of anisotropic kinetics with the spherical angles $\theta$ and $\phi$, which define the orientation of the normal to the dendrite interface and its growth direction.

In the case of cubic symmetry, $d(\theta, \phi)$ and $\tilde{\beta}(\theta, \phi)$ are described by

$$d(\theta, \phi) = d_0 \{1 - \alpha_d [\cos^4 \theta + \sin^4 \theta (1 - 2 \sin^2 \phi \cos^2 \phi)]\}, \quad (5)$$

$$\tilde{\beta}(\theta, \phi) = \beta_0 \{1 - \alpha_{\beta} [\cos^4 \theta + \sin^4 \theta (1 - 2 \sin^2 \phi \cos^2 \phi)]\}, \quad (6)$$

where $d_0$ is the capillary constant, $\alpha_d \ll 1$ stands for the stiffness, which depends on a small anisotropy parameter $\epsilon_c$ of surface energy, $\beta_0$ is the kinetic constant, and $\alpha_\beta \ll 1$ is the kinetic anisotropy parameter.

Considering the case of needle-like crystal in the form of a paraboloid of revolution, Eq. 5 can be reduced by averaging over $\phi$¹⁸ and written out for the case of $n$-fold symmetry:

$$d(\theta) = d_0 \{1 - \alpha_d \cos[n(\theta - \theta_d)]\}, \quad (7)$$

$$\tilde{\beta}(\theta) = \beta_0 T_Q \{1 - \alpha_\beta \cos[n(\theta - \theta_\beta)]\}, \quad (8)$$

where $\theta_d$ and $\theta_\beta$ designate the angles between the directions of growth and minimal functions $d(\theta)$ and $\tilde{\beta}(\theta)$.

The convective heat transfer problem (1)–(8) written out for the case of $n$-fold symmetry of crystal growth enables us to find the generalized solvability criterion including convection and the effects of kinetics.

### Analytical Solution

In the first instance, we define the corresponding parabolic coordinates $\xi$ and $\eta$ connected with the Cartesian ones, $x$ and $z$, by means of the familiar expressions

$$x = \rho \sqrt{\xi \eta}, \quad z = \frac{\rho}{2} (\eta - \xi), \quad (9)$$

where $\rho/2$ represents the dendritic tip radius. The equation $\eta = 1$ corresponds to the solid/liquid surface of the dendrite.

The Oseen hydrodynamic Eq. 3 supplemented by the corresponding no-slip boundary conditions for the fluid velocity components $u_\eta$ and $u_\xi$ in a parabolic reference frame take the form¹⁷,¹⁹

$$u_{\eta} = - \frac{f(\eta)}{2\sqrt{\xi} + \eta}, \quad u_{\xi} = \frac{\sqrt{\xi} \eta}{\sqrt{\xi} + \eta} \frac{df(\eta)}{d\eta}, \quad (10)$$
where

$$
\begin{align}
f(\eta) &= 2(U + V)\sqrt{\eta} - 2Ug(\eta), \nonumber \\
g(\eta) &= \sqrt{\eta} \frac{erfc\sqrt{\eta\Re/2}}{erfc\sqrt{\Re/2}} + \frac{\sqrt{2/(\pi\Re)}}{erfc\sqrt{\Re/2}} \tag{11} \\
&\quad \times \left[ \exp\left(-\frac{\Re}{2}\right) - \exp\left(-\frac{\eta\Re}{2}\right) \right] \nonumber
\end{align}
$$

with the Reynolds number $\mathcal{R} = \rho U/v$.

Now rewriting the heat transfer Eq. 2 as well as the boundary condition (1) in parabolic coordinates, and then integrating them, we find the temperature distribution in the two-dimensional geometry as

$$
T_l(\eta) = T_i + (T_\infty - T_i) \frac{I_T(\eta)}{I_T(\infty)}, \quad (12)
$$

where

$$
\begin{align*}
I_T(\eta) &= \int_1^\eta \exp \left[ P_f \int_1^{\eta'} \frac{g(\eta'')}{\sqrt{\eta''}} d\eta'' - P_0 \eta' \right] \frac{d\eta'}{\sqrt{\eta'}}, \\
T_i &= T_\infty + T_Q P_g \exp(P_0) I_T(\infty), \quad P_0 = P_g + P_f.
\end{align*}
$$

$T_\infty$ is the temperature in the liquid phase far from the dendrite surface, and $P_g = \rho V/(2D_T)$ and $P_f = \rho U/(2D_T)$ stand for the growth and flow Péclet numbers.

**Solvability Condition**

Pelcé, Bouissou and Bensimon<sup>3,19,20</sup> showed that the solvability condition, which gives a unique combination between $\rho$ and $V$, describes a stable dendritic growth mode as

$$
\int_{-\infty}^{\infty} G[X_0(l)] Y_m(l) dl = 0, \quad Y_m(l) = \exp \left[ i \int_{0}^{l} k_m(l_1) dl_1 \right], \quad (13)
$$

where $G$ is the curvature operator, $k_m(l)$ designates the marginal wavenumber mode, $i$ represents the imaginary unit, and $X_0(l)$ is a continuum of solutions from which the dependence $k_m(l)$ can be derived.

Furthermore, to obtain the marginal wavenum-
bers $k_m$ entering in the solvability integral (13), the
linear stability analysis should be carried out (see,
among others,¹⁴). It leads to the marginal mode of
the wavenumber $k_m$ (see, for details, ¹⁹,²¹), which is
determined by the following cubic equation

$$
\begin{equation}
\begin{split}
k_m^3 ={}& \frac{V \exp(i\theta)}{2d(\theta)D_T} k_m + \frac{i a U \sin\theta \cos\theta}{8\rho D_T} k_m \\
& - \frac{iV \sin\theta}{2D_T} k_m^2 + \frac{V^2 \cos\theta \exp(i\theta)}{4d(\theta)D_T^2} + \frac{iV \tilde{\beta}(\theta) \sin\theta}{d(\theta)T_Q} k_m^2.
\end{split}
\tag{14}
\end{equation}
$$

An analytical solution of the cubic equation (14) has
been found using Cardano’s formula ($d_0 \neq 0$) for the
fourfold symmetry of dendritic growth as<sup>22</sup>

$$
k_m = k_{mBP} + \frac{V \exp(-i\theta)}{4D_T} + \frac{iV\tilde{\beta}(\theta)\sin\theta}{2d(\theta)T_Q}, \quad (15)
$$

where

$$
k_{mBP} = - \sqrt{\frac{V \exp(i\theta)}{2d(\theta)D_T}} \sqrt{1 + \frac{i a U d(\theta)}{8 \rho V}} \sin(2\theta) \exp(-i\theta).
$$

**Solvability Criterion for the Thermally**
**Controlled Growth**

Let us consider the purely thermal mode of
dendritic solidification when $\alpha_d \ll 1$, $\alpha_\beta \ll 1$, $\theta_d = 0$
and $\theta_\beta \neq 0$. Substituting the analytical solution for
the fourfold symmetry of the dendritic growth
wavenumber (15) into the solvability integral (13),
we arrive at

$$
\int_{-\infty}^{\infty} d\eta G[X_0(\eta)] \exp \left\{ -i \int_0^{\eta} \left[ \frac{P_g(1 - i\eta_1)}{2} + \frac{i\rho V \beta_0 \eta_1}{2d_0} - \sqrt{\frac{(1 + i\eta_1)(1 + \eta_1^2)^{\lambda_n} + i\alpha\eta_1 B_n(\eta_1)}{\sigma^* B_n(\eta_1)}} \right] d\eta_1 \right\} = 0, \quad (16)
$$

where

$$
\begin{align*}
l_1 &= -\frac{\rho}{2} \left[ \frac{\tan \theta}{\cos \theta} + \ln \left( \frac{1}{\cos \theta} + \tan \theta \right) \right], \quad \eta = \tan \theta, \\
\alpha &= \frac{aUd_0}{4\rho V}, \quad \lambda_n = \frac{n+1}{2}, \quad \sigma^* = \frac{2d_0D_T}{\rho^2V}, \\
B_n(\eta) &= (1 + \eta^2)^{n/2} - \alpha_d \sum_{k=0}^{n} \binom{n}{k} \eta^{n-k} \cos \frac{(n-k)\pi}{2}.
\end{align*}
$$

The wave-number in the kinetically limited regime
can be obtained under the assumption that $\theta_\beta = 0$
and $\theta_d \neq 0$:¹⁴
$$ \int_{-\infty}^{\infty} d\eta G[X_0(\eta)] \exp \left\{ -i \int_0^{\eta} \left[ P_g + \frac{i\rho(1+i\eta_1)(1+\eta_1^2)^{\lambda_n}}{2D_T\beta_0\eta_1 B'_n(\eta_1)} \right] d\eta_1 \right\} = 0, \\ B'_n(\eta) = (1+\eta^2)^{n/2} - \alpha_\beta \sum_{k=0}^{n} \binom{n}{k} \eta^{n-k} \cos \frac{(n-k)\pi}{2}. \quad (17) $$

We calculate the solvability integral (16) in two stages as detailed in$^{22}$. At the first stage, we neglect the kinetic contribution (proportional to $\beta_0$ in (16)) and pay our attention to the case known as “thermally controlled” crystal growth. Setting $\beta_0 = 0$, we come to the selection criterion (see also$^{19,21,23}$)

$$ \sigma^* = \frac{2d_0 D_T}{\rho^2 V} = \frac{\sigma_0 \alpha_d^{7/n} A_n^{7/n}}{\left(1 + a_1 \alpha_d^{2/n} A_n^{2/n} P_g\right)^2 (1 + b \tau_n^{\nu_n})}, \quad (18) $$

where

$$ \tau_n = \alpha \alpha_d^{-3/n} A_n^{-3/n}, \quad \nu_n = \frac{n+7}{2(n+3)}, $$

$$ A_n = 2^{-3n/4} \sum_{k=0}^{n} \binom{n}{k} i^{n-k} \cos \frac{(n-k)\pi}{2}, $$

and $\sigma_0$ and $b$ are the constants.

The second stage to evaluate the integral expression (16) is to analyze the dendrite growth mode controlled by the kinetic contribution (proportional to $\beta_0$). By doing that, we neglect the summand proportional to the growth Péclet number in (16) and arrive at (by analogy with$^{22}$)

$$ \sigma^* = \frac{\sigma_0 \alpha_d^{7/n} A_n^{7/n}}{\left[1 + a'_1 \alpha_d^{2/n} P_g D_T \beta_0 A_n^{2/n} / d_0\right]^2 (1 + b \tau_n^{\nu_n})}, \quad (19) $$

where $a'_1$ represents a new constant.

Sewing together the obtained limiting criteria (18) and (19), we come to the generalized criterion of the form

$$ \sigma^* = \frac{2d_0 D_T}{\rho^2 V} = \frac{\sigma_0 \alpha_d^{7/n} A_n^{7/n}}{\left[1 + a_1 \alpha_d^{2/n} A_n^{2/n} P_g \left(1 + \delta_0 D_T \beta_0 / d_0\right)\right]^2 \left(1 + b \tau_n^{\nu_n}\right)}, \quad (20) $$

where $a'_1 = a_1 \delta_0$ and $\delta_0$ represents a fitting constant. Note that the selection constants $\sigma_0$ and $b$ for each $n$ might be found experimentally or from the phase-field modeling.

## Undercooling Balance

The second relation for $V$ and $\rho$ is called the undercooling balance. This second relation and selection criterion (20) allow us to obtain a pair of most important parameters of primary solidification, $V$ and $\rho$, at a unique undercooling $\Delta T$. The total undercooling balance connects the melting

temperature $T_m$ of a one-component liquid and the far-field temperature $T_\infty$ as $\Delta T = T_m - T_\infty$ and introduces the first model equation, which consists of several contributions:

$$ \Delta T = \Delta T_T + \Delta T_R + \Delta T_K, \quad (21) $$

where $\Delta T_T$ represents the thermal contribution, $\Delta T_R = 2d_0 T_Q/R$ is the undercooling due to the Gibbs-Thomson effect, and $\Delta T_K = V/\mu_k$ is the kinetic undercooling (where $\mu_k$ stands for the kinetic coefficient).

The thermal contribution $\Delta T_T$ can be written out using the Ivantsov function $Iv_T$, which describes the temperature field around the growing steady-state dendrite of a parabolic form:

$$ \Delta T_T = T_Q Iv_T(P_g, P_f), \quad (22) $$

where the Ivantsov function

$$ Iv_T(P_g, P_f) = P_g \exp(P_0) I_T(\infty), \quad P_0 = P_g + P_f \quad (23) $$

depends on the growth $P_g = \rho V/(2D_T)$ and flow $P_g = \rho U/(2D_T)$ Péclet numbers.

Expression (21) represents a function of two variables, velocity $V$ and diameter $\rho$, through the functions $\Delta T_T$, $\Delta T_R$ and $\Delta T_K$ for a given full undercooling $\Delta T$. To obtain $V$ and $\rho$ simultaneously, we use the second equation in the form of stability criterion (20). These two equations allow us to obtain $V$ and $\rho$ for a given $\Delta T$.

## METHODOLOGY

In this work, an enthalpy-based method that captures solidification is coupled to a lattice Boltzmann method (LBM) that resolves fluid flow. The methods are linked by the convective thermal transport equation and solid fraction, although in other work they can also be linked by body forces.$^{24}$ The enthalpy-based method is based on the work of Voller,$^{25}$ and before that Tacke,$^{26,27}$ and it is written using a finite difference scheme. The LBM uses a discretized form of the Boltzmann equation. This section outlines the two methods and then describes the problem setup and typical solution procedure for generating steady-state solutions for dendrite tip velocity and tip radius for a given bulk undercooling and mean flow.

## Enthalpy Method

Convective transport of enthalpy is governed by

$$ \frac{\partial H}{\partial t} = \nabla^2(KT) - \mathbf{u}(\nabla H), \quad (24) $$

where $t$ is time, $K$ is the thermal conductivity, which is assumed to be the same for both liquid and solid phases, $\mathbf{u}$ is velocity, and $H$ is enthalpy, which is a
function of the heat capacity $c_p$, temperature $T$ and the latent $L$. The volumetric enthalpy is defined as a sum of sensible and latent heats

$$H = c_p T + f_{\text{liq}} L, \qquad (25)$$

where $f_{\text{liq}}$ is the liquid fraction phase variable and is 0 for fully solid and 1 for fully liquid. At the solid-liquid interface, the Gibbs-Thompson condition gives the temperature for a pure material neglecting any kinetic effects as

$$T^i = T_f - \frac{\Gamma(\theta)T_f}{L}, \qquad (26)$$

where $T_f$ is the fusion temperature for a surface energy anisotropy, $\gamma$, of the form $\gamma = d_0(1 + \epsilon_4 \cos 4\theta)$, the surface stiffness $\Gamma(\theta)$ given by

$$\Gamma(\theta) = \gamma + \frac{\partial^2 \gamma}{\partial \theta^2} = d_0(1 - \epsilon_4 \cos 4\theta). \qquad (27)$$

and $\kappa$ is the interface curvature that can be captured from the liquid fraction gradients as

$$\kappa = \nabla \cdot \frac{\nabla f_{\text{liq}}}{|\nabla f_{\text{liq}}|} = (f_y^2 f_{xx} - 2 f_x f_y f_{xy} + f_x^2 f_{yy}) \cdot (f_x^2 + f_y^2)^{-\frac{3}{2}}. \qquad (28)$$

Here, $f_x$ and $f_{xx}$ represent the first and the second derivatives of $f_{\text{liq}}$ with respect to the coordinate $x$. The interface orientation $\theta$, which is the angle between the interface normal and the $x$-axis, is given by

$$\theta = \arctan \frac{f_x}{f_y}. \qquad (29)$$

## Lattice Boltzmann Method

The momentum transport for an incompressible flow is given by the Navier-Stokes equations (NSEs) as

$$\nabla \cdot \mathbf{u} = 0, \qquad (30)$$

$$\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = -\frac{1}{\rho_l} \nabla p + \nu \nabla^2 \mathbf{u}, \qquad (31)$$

where $\rho_l$ is density, $p$ is pressure, and $\nu$ is the kinematic viscosity. To avoid non-linearity in the convective term and to improve convergence of the problem, the lattice Boltzmann method is employed. It describes the evolution of a particle distribution function (PDF), $f_i$, and it has been shown to recover the NSEs in the low Mach number limit via the Chapman and Enskog multi-scale analysis. The implementation is 3D, using a D3Q19 lattice, which

in this two-dimensional analysis reduces to the equivalent D2Q9 lattice. The lattice Boltzmann equation generally can be written as

$$f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) - f_i(\mathbf{x}, t) = -\frac{1}{\tau} (f_i(\mathbf{x}, t) - f_i^{\text{eq}}(\mathbf{x}, t)), \qquad (32)$$

where the left-hand side represents the streaming process and the right-hand side describes collisions. Here, $\mathbf{x}$ is the lattice node coordinate, $\mathbf{c}_i$ are the discretized lattice velocities, $\Delta t$ is the LBM time step, $\tau$ is a non-dimensional relaxation parameter, and $f_i^{\text{eq}}$ is the equilibrium PDF given by

$$f_i^{\text{eq}} = \rho^* w_i \left( 1 + \frac{3\mathbf{c}_i \cdot \mathbf{u}^*}{c^2} + \frac{(3\mathbf{c}_i \cdot \mathbf{u}^*)^2}{2c^4} - \frac{3\mathbf{u}^{*2}}{2c^2} \right), \qquad (33)$$

where the asterisks (*) mark the non-dimensional lattice Boltzmann variables, $w_i$ is the lattice weight coefficient, and $c = \Delta x/\Delta t$ is the lattice speed. The lattice spacing, $\Delta x$, time stepping, $\Delta t$ and equilibrium density, $\rho^*$, are all defined as unity. Fluid properties such as density and fluid velocities can be calculated from the PDF by taking the velocity moments as

$$\rho^* = \sum_i f_i, \quad \rho^* \mathbf{u}^* = \sum_i f_i \mathbf{c}_i. \qquad (34)$$

Other physical quantities like the fluid density and viscosity can be calculated from the following expressions:

$$p^* = \frac{\rho^*}{3}, \quad v^* = \frac{1}{3} \left( \tau - \frac{1}{2} \right). \qquad (35)$$

For stability purposes, the two-relaxation-time collision scheme<sup>28</sup> is used in the model. The moment-based boundary method<sup>29,30</sup> describes the flat domain boundaries while the bounce-back scheme describes the advancing solid front of the growing crystal.

## Problem Setup and Solution Procedure

The main advantage of the enthalpy method is that it is derived without explicitly using the atomic scale interface thickness, in contrast to traditional phase-field methods. A single-cell-thick interface is still required to track phase change; however, this method can be used at both the macro- and micro-scales. In the context of this work, the method can encompass a wide range of undercoolings. However, the drawback is that the method suffers from two errors that are a consequence of using a single cell to represent the interface and related to the cell size and the problem considered. If the cell size is too large, then the error analogous to using too coarse of a mesh becomes significant; however, if the cell size is too small, the so-called narrow band error feature appears. The latter has a significant effect on
calculating interfacial features, namely curvature, where with refinement the local features become either flat or a corner representing either a zero or a very large curvature, respectively. Consequently, this causes the formation of an unrealistic faceted-like structure. To address this problem, the authors developed an adaptive cell size method for the enthalpy method, where the tip radius was defined by a fixed number of cells. However, as the tip radius is unknown, so too is the cell size, and it becomes a solved variable. For brevity, the full method is not repeated here, but in summary, through an iterative approach, a steady-state solution for a given undercooling was found when the cell size, tip radius (hence curvature) and tip velocity all became constant in time. The same approach is used here, but with the addition of forced convection. Aside from resolving the fluid flow, no significant changes are necessary for this adaptive cell size method. There are three key solvers that are weakly coupled with each solved sequentially. A single time step consists of first resolving the thermal transport (24), which determines the local free energy for solidification, which is resolved by the enthalpy equation with input from the Gibbs-Thompson condition (26), and finally fluid flow is calculated. The fluid velocity then feeds into the transport equation. Therefore, while fluid flow has a significant effect on solidification, it does not directly feed into the calculation of tip velocity, curvature or the update of the cell size.

The numerical model represents a square domain of 1200 × 1200 cells with a dendrite tip growing from the origin along the positive x-axis. Due to the symmetry of the problem along the x-axis, only half of the problem is solved with the negative y boundary acting as a symmetry plane. The far field boundary conditions (positive x and positive y boundaries) are set to the bulk undercooling and mean fluid flow values. The negative x boundary acts as a zero Neumann condition. To prevent the

dendrite tip from growing too close and being influenced by the far field boundaries, a moving mesh technique is used that keeps the dendrite in the same relative position within the domain. Both the numerical model and analytic solution use the same physical parameters given in Table I.

## RESULTS AND CONCLUSIONS

An example of a steady-state solution for 1.08 m/s incident flow is given in Fig. 2, showing the dendrite tip morphology and the x and y components of the velocity profile (u and v in Fig. 2a and b respectively). As the dendrite is solid, fluid flow must go around, and this leads to a stagnation point at the very tip with higher velocities inside the thermal boundary layer. However, as the relative velocity at the interface is zero, a viscous boundary layer forms between the interface and the higher velocity region.

Figure 3 shows a comparison of the thermal field for $\Delta T = 106$ K with an incident flow of 0 m/s, 1.08 m/s and 8.66 m/s. In the presence of an incident fluid flow, advective thermal transport causes the isotherms to bunch up, leading to an increase in thermal gradient in the liquid. This in turn allows for a higher diffusion rate of temperature back into the liquid and therefore an increase in free energy. In all cases with increasing flow velocity, there is an increase in tip velocity and consequently a reduction in tip radius; the latter leads to an increase in tip undercooling, highlighted by Fig. 3d, e and f. In Figs. 2 and 3, the axes are in terms of cell numbers, where the adaptive cell size method¹⁰ has the tip radius defined as eight cells, and therefore with increasing velocity the steady-state cell size decreases. Figures 2 and 3 are also focused on the lower left corner of the computational domain encompassing the first 360 cells in x and y. The full domain is much larger so that the far field boundaries do not have any direct influence on the tip dynamics.

<table><caption>**Table I.** Physical and numerical parameters used in the present analytical calculations and numerical simulations for the pure nickel</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Calc.</th><th>EnthM</th><th>Units</th></tr></thead><tbody><tr><td>Melting temperature</td><td>T<sub>m</sub></td><td>1728</td><td>1728</td><td>K</td></tr><tr><td>Adiabatic temperature</td><td>T<sub>Q</sub></td><td>435</td><td>435</td><td>K</td></tr><tr><td>Thermal diffusivity</td><td>D<sub>T</sub></td><td>1</td><td>1</td><td>·10<sup>-5</sup> m<sup>2</sup>/s</td></tr><tr><td>Density</td><td>ρ<sub>l</sub></td><td>7.9</td><td>7.9</td><td>·10<sup>3</sup> kg/m<sup>3</sup></td></tr><tr><td>Capillary length</td><td>d<sub>0</sub></td><td>4</td><td>4</td><td>·10<sup>-10</sup> m</td></tr><tr><td>Anisotropy strength</td><td>β</td><td>0.018</td><td>0.018</td><td>—</td></tr><tr><td>Selection parameter</td><td>a<sub>1</sub></td><td>0.5</td><td>—</td><td>—</td></tr><tr><td>Selection constant<sup>a</sup></td><td>b</td><td>0.1</td><td>—</td><td>—</td></tr><tr><td>Selection constant<sup>a</sup></td><td>σ<sub>0</sub></td><td>0.03</td><td>—</td><td>—</td></tr></tbody></table>

Here, Calc. means the calculated sharp interface model, and EnthM means the enthalpy methodª Note that the selection constant b is chosen on the basis of our previous works²¹,³¹ while another constant $\sigma_0$ represents a fitting parameter
<figure><img src="image_4.png" /><figcaption>Fig. 2. Incident fluid flow velocity of 1.08 m/s onto the dendrite tip for $\Delta T = 106$ K. **a** Absolute of x-component of velocity. **b** y-component of velocity.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Thermal field for $\Delta T = 106$ K with incident flow velocity of **a**, **d** 0 m/s, **b**, **e** 1.08 m/s and **c**, **f** 8.66 m/s. **d–f** Reduced temperature scaled to highlight tip undercooling. The circle in **d** corresponds to the defined tip radius for all cases.</figcaption></figure>

The stability criterion (20) and undercooling balance condition (21) represent a set of two nonlinear equations. Solving these two simulatneous equations gives the tip velocity $V$ and tip diameter $\rho$ as functions of the melt undercooling $\Delta T$. The solution of these is compared with the enthalpy method in Fig. 4 for pure nickel at different flow velocities $U$. The results show a good agreement between the theoretical and numerical models over a broad range of fluid velocities for fixed values of selection constants $\sigma_0$ and $b$ listed in Table I.

Analyzing a series of curves shown in Fig. 4, we conclude that the dendrite tip radius $\rho/2$ decreases and the tip velocity $V$ increases with increasing the fluid velocity $U$ at a fixed value of undercooling $\Delta T$. This means that the dendritic crystals become narrower and grow faster as a forced flow of melt rises.

This article tests the convective theory of stable dendritic growth with computer simulation data for pure (impurity-free) melts. The results highlight that the enthalpy-based method, along
<figure><img src="image_6.png" /><figcaption>Fig. 4. **a** Dendrite tip velocity and **b** tip radius for pure nickel versus the melt undercooling for different flow velocities *U*. Physical parameters are listed in Table I.</figcaption></figure>

with an adaptive cell size approach, is capable of accurately predicting the behavior of dendritic growth in the presence of forced convection. In practice, fluid flow is generally present in terrestrial experiments. Consequently, further understanding its effect is necessary to bridge knowledge between solidification under terrestrial and microgravity conditions.

The next step is to test the theory for binary systems with the modified stability criterion and undercooling balance written out with allowance for impurity effects (see, for details, the theory developed in Refs. 14,22,23,32 for a binary melt). Furthermore, the effects of non-equilibrium crystallization occurring at high growth rates should be tested against experimental data according to the theory developed in Ref. 33. Finally, a full three-dimensional approach will allow for direct comparison of the numerical model, analytic theory and experimental evidence.

## ACKNOWLEDGEMENT

D.V. Alexandrov and P.K. Galenko acknowledge the support from the Russian Science Foundation (Grant No. 16-11-10095) and support from the German Space Center Space Managment under contract no. 50WM1941.

## CONFLICT OF INTEREST

The authors declare that they have no conflict of interest.

## OPEN ACCESS

This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## REFERENCES

1. A.A. Chernov, *Mod. Crystallogr. III.* (Springer, Berlin, 1984).

2. D.A. Kessler, J. Koplik, and H. Levine, *Adv. Phys.*, **37**, 255 (1988).

3. P. Pelcé, *Dynamics of Curved Fronts*. (Academic Press, Boston, 1988).

4. D.V. Alexandrov and P.K. Galenko, *J. Phys. A Math. Theor.* **46**, 195101 (2013).

5. K. Libbrecht, *Snowflakes*. (Voyageur Press, Minneapolis, 2004).

6. V.Y. Shur and A.R. Akhmatkhanov, *Philos. Trans. R. Soc. A* **376**, 20170204 (2018).

7. O.P. Fedorov, *Crystal Growth Processes: Kinetics, Formation, Inhomogeneities*. (Naukova Dumka, Kiev, 2010).

8. A.J. Melendez and C. Beckermann, *J. Cryst. Growth* **340**, 175 (2012).

9. D.V. Alexandrov and P.K. Galenko, *Philos. Trans. R. Soc. A* **378**, 20190243 (2020).

10. A. Kao, L.V. Toropova, D.V. Alexandrov, G. Demange, and P.K. Galenko, *J. Phys. Condens. Matter* **32**, 194002 (2020).

11. L.V. Toropova, D.V. Alexandrov, M. Rettenmayr, and P.K. Galenko, *J. Cryst. Growth* **535**, 194002 (2020).

12. L.V. Toropova, P.K. Galenko, D.V. Alexandrov, G. Demange, A. Kao, and M. Rettenmayr, *Eur. Phys. J. Spec. Top.* **229**, 275 (2020).

13. D.V. Alexandrov and P.K. Galenko, *J. Phys. Chem. Solids* **108**, 98 (2017).

14. D.V. Alexandrov, P.K. Galenko, and L.V. Toropova, *Philos. Trans. R. Soc. A* **376**, 20170215 (2018).

15. H. (Sir) Lamb, *Hydrodynamics*. (Dover Publication, New York, 1945).

16. N.E. Kochin, I.A. Kibel, and N.V. Roze, *Theoretical Hydromechanics*. (Interscience, New York, 1964).

17. S.K. Dash and W.N. Gill, *Int. J. Heat Mass Transf.* **27**, 1345 (1984).
18. A. Barbieri and J.S. Langer, *Phys. Rev.* A 39, 5314 (1989).

19. Ph. Bouissou and P. Pelcé, *Phys. Rev.* A 40, 6673 (1989).

20. P. Pelcé and D. Bensimon, *Nucl. Phys. B* 2, 259 (1987).

21. D.V. Alexandrov and P.K. Galenko, *Phys. Usp.* 57, 771 (2014).

22. D.V. Alexandrov and P.K. Galenko, *Phys. Chem. Chem. Phys.* 17, 19149 (2015).

23. D.V. Alexandrov and P.K. Galenko, *Phys. Rev. E* 87, 062403 (2013).

24. A. Kao, J. Gao, and K. Pericleous, *Philos. Trans. R. Soc. A* 376, 20170206 (2018).

25. V.R. Voller, *Int. J. Heat Mass Transf.* 51, 823 (2008).

26. K.-H. Tacke, K.-H. Hofmann, J. Sprekels (Eds.), *Free Boundary Problems: Theory and Applications*. (Longman, Harlow, 1990), 636.

27. K.-H. Tacke, A. Harnisch, *Proceedings of the International Conference on Computational Modeling of Free and Moving Boundary Problems*. (WIT press, Southampton, 1991).

28. I. Ginzburg, F. Verhaeghe, and D. d'Humieres, *Commun. Comput. Phys.* 3, 427 (2008).

29. S. Bennett, P. Asinari, and P. Dellar, *Int. J. Numer. Methods Fluids*, 69, 171 (2012).

30. S. Mohammed and T. Reis, *Arch. Mech. Eng.* 64, 57 (2017).

31. J. Gao, M. Han, A. Kao, K. Pericleous, D.V. Alexandrov, and P.K. Galenko, *Acta Mater.*, 103, 184 (2016).

32. D.V. Alexandrov and P.K. Galenko, *J. Phys. Condens. Matter* 30, 105702 (2018).

33. D.V. Alexandrov and P.K. Galenko, *Proceedings of IUTAM Symposium on Recent Advances in Moving Boundary Problems in Mechanics*, ed. S. Gutschmidt, J.N. Hewett, M. Sellier (Springer, Berlin, 2019), Vol. 34, pp. 203–215.

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.