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

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 [1][2][3][4][5][6] . As this takes place, the growth mechanisms of dendrites, their shape and interaction determine the characteristics of the internal microstructure of the crystallized substance [7][8][9] . 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 q depending on the melt undercooling DT. 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 q as functions of DT 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 [10][11][12][13] 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 q. Section 3 is concerned with the background of the enthalpy method and our main results following from simulations of dendritic growth. The main outcomes of our analysis and future directions are discussed in Sect. 4.

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 14 (Fig. 1).
The heat balance condition at the dendritic interface can be written as where T Q is the hypercooling, v Á 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 where w is the fluid velocity, t is the time variable, and r 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 [15][16][17] Here, U is the fluid velocity far from the dendritic surface (see Fig. 1), p is the pressure, q l is the density of liquid, and m is the kinematic viscosity. The Gibbs-Thomson equation at the solid-liquid boundary holds where T int is the phase transition temperature at the dendrite interface, T m is the melting temperature for the pure system, and K is the interface curvature, dðh; /Þ, andbðh; /Þ are the capillary length and the function of anisotropic kinetics with the spherical angles h and /, which define the orientation of the normal to the dendrite interface and its growth direction.
In the case of cubic symmetry, dðh; /Þ andbðh; /Þ are described by where d 0 is the capillary constant, a d ( 1 stands for the stiffness, which depends on a small anisotropy parameter e c of surface energy, b 0 is the kinetic constant, and a b ( 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 / 18 and written out for the case of n-fold symmetry: where h d and h b designate the angles between the directions of growth and minimal functions dðhÞ and bðhÞ.
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 n and g connected with the Cartesian ones, x and z, by means of the familiar expressions where q=2 represents the dendritic tip radius. The equation g ¼ 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 g and u n in a parabolic reference frame take the form 17,19  where f ðgÞ ¼ 2ðU þ VÞ ffiffi ffi g p À 2UgðgÞ; with the Reynolds number R ¼ qU=m. 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 where T 1 is the temperature in the liquid phase far from the dendrite surface, and P g ¼ qV=ð2D T Þ and P f ¼ qU=ð2D T Þ stand for the growth and flow Péclet numbers.

Solvability Condition
Pelcé, Bouissou and Bensimon 3,19,20 showed that the solvability condition, which gives a unique combination between q and V, describes a stable dendritic growth mode as 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 wavenumbers k m entering in the solvability integral (13), the linear stability analysis should be carried out (see, among others, 14 ). It leads to the marginal mode of the wavenumber k m (see, for details, 19,21 ), which is determined by the following cubic equation An analytical solution of the cubic equation (14) has been found using Cardano's formula (d 0 6 ¼ 0) for the fourfold symmetry of dendritic growth as 22 where

Solvability Criterion for the Thermally Controlled Growth
Let us consider the purely thermal mode of dendritic solidification when a d ( 1, a b ( 1, h d ¼ 0 and h b 6 ¼ 0. Substituting the analytical solution for the fourfold symmetry of the dendritic growth wavenumber (15) into the solvability integral (13), we arrive at where The wave-number in the kinetically limited regime can be obtained under the assumption that h b ¼ 0 and h d 6 ¼ 0: 14 We calculate the solvability integral (16) in two stages as detailed in 22 . At the first stage, we neglect the kinetic contribution (proportional to b 0 in (16)) and pay our attention to the case known as ''thermally controlled'' crystal growth. Setting b 0 ¼ 0, we come to the selection criterion (see also 19,21,23 ) where and r 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 b 0 ). By doing that, we neglect the summand proportional to the growth Péclet number in (16) and arrive at (by analogy with 22 ) where a 0 1 represents a new constant. Sewing together the obtained limiting criteria (18) and (19), we come to the generalized criterion of the form where a 0 1 ¼ a 1 d 0 and d 0 represents a fitting constant. Note that the selection constants r 0 and b for each n might be found experimentally or from the phase-field modeling.

Undercooling Balance
The second relation for V and q 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 q, at a unique undercooling DT. The total undercooling balance connects the melting temperature T m of a one-component liquid and the far-field temperature T 1 as DT ¼ T m À T 1 and introduces the first model equation, which consists of several contributions: where DT T represents the thermal contribution, DT R ¼ 2d 0 T Q =R is the undercooling due to the Gibbs-Thomson effect, and DT K ¼ V=l k is the kinetic undercooling (where l k stands for the kinetic coefficient). The thermal contribution DT 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: where the Ivantsov function depends on the growth P g ¼ qV=ð2D T Þ and flow P g ¼ qU=ð2D T Þ Péclet numbers. Expression (21) represents a function of two variables, velocity V and diameter q, through the functions DT T , DT R and DT K for a given full undercooling DT. To obtain V and q simultaneously, we use the second equation in the form of stability criterion (20). These two equations allow us to obtain V and q for a given DT.

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 where t is time, K is the thermal conductivity, which is assumed to be the same for both liquid and solid phases, 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 where f liq is the liquid fraction phase variable and is 0 for fully solid and 1 for fully liquid. At the solidliquid interface, the Gibbs-Thompson condition gives the temperature for a pure material neglecting any kinetic effects as where T f is the fusion temperature for a surface energy anisotropy, c, of the form c ¼ d 0 ð1 þ 4 cos 4hÞ, the surface stiffness CðhÞ given by and j is the interface curvature that can be captured from the liquid fraction gradients as Here, f x and f xx represent the first and the second derivatives of f liq with respect to the coordinate x.
The interface orientation h, which is the angle between the interface normal and the x-axis, is given by Lattice Boltzmann Method The momentum transport for an incompressible flow is given by the Navier-Stokes equations (NSEs) as where q l is density, p is pressure, and m 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 where the left-hand side represents the streaming process and the right-hand side describes collisions.
Here, x is the lattice node coordinate, c i are the discretized lattice velocities, Dt is the LBM time step, s is a non-dimensional relaxation parameter, and f i eq is the equilibrium PDF given by where the asterisks (Ã) mark the non-dimensional lattice Boltzmann variables, w i is the lattice weight coefficient, and c ¼ Dx=Dt is the lattice speed. The lattice spacing, Dx, time stepping, Dt and equilibrium density, q Ã , 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 Other physical quantities like the fluid density and viscosity can be calculated from the following expressions: For stability purposes, the two-relaxation-time collision scheme 28 is used in the model. The momentbased boundary method 29,30 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 microscales. 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 facetedlike 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 DT ¼ 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 10 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. 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 q as functions of the melt undercooling DT. 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 r 0 and b listed in Table I. Analyzing a series of curves shown in Fig. 4, we conclude that the dendrite tip radius q=2 decreases and the tip velocity V increases with increasing the fluid velocity U at a fixed value of undercooling DT. 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  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 terrestial 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.  Table I.