Diffuse interface models of solidification with convection: The choice of a finite interface thickness

The thin interface limit aims at minimizing the effects arising from a numerical interface thickness, inherent in diffuse interface models of solidification and microstructure evolution such as the phase field model. While the original formulation of this problem is restricted to transport by diffusion, we consider here the case of melt convection. Using an analysis of the coupled phase field-fluid dynamic equations, we show here that such a thin interface limit does also exist if transport contains both diffusion and convection. This prediction is tested by comparing simulation studies, which make use of the thin-interface condition, with an analytic sharp-interface theory for dendritic tip growth under convection.


Introduction
Letting aside critical phenomena, physical interfaces often have a width in the nanometer range. For problems on the mesoscale (i.e., dealing with micrometers or larger scales), this thickness is negligible and the physical interface can safely be approximated to as a mathematically sharp boundary separating the phases of interest. The major aim of modeling at the mesoscale is thus to solve problems involving a sharp interface (SI). On the other hand, in the past thirty years, the so-called diffuse interface models such as the phase field (PF) approach [1,2] have proved quite powerful in studying solidification and microstructure evolution. These models involve a finite interface thickness, W , which, in view of the above mentioned fact, is of a numerical nature. Ideally, one would like to minimize the effect of this numerical parameter. The equivalence of a PF model of solidification to the SI formulation was established by Caginalp [3] as the diffuse interface becomes progressively narrow a e-mail: fathollah.varnik@rub.de (W → 0). This ideal limit, however, is numerically quite expensive and often impractical. A major advancement was achieved in the mid 1990s by Karma and Rappel [4] with the so-called thin-interface limit for problems involving diffusive transport. They found that instead of vanishing interface thickness, it is sufficient to have W small compared to the diffusion length of the solidification problem. This diffusion length is defined as L d = D V , where D is the thermal diffusivity and V is the normal velocity of the interface.

Thin interface analysis in the presence of flow
To account for transport due to convection in the solidification phenomena, a couple of melt flow and PF couplings have been proposed and analyzed. Anderson et al. performed [5] a sharp interface asymptotics of a PF model where the viscosity of liquid melt-solid interface diverges while approaching solid end of the interface (known as the variable viscosity model [6][7][8]). Beckermann et al. [9] proposed a dissipative drag force ansatz that acts a momentum sink within the liquid melt-solid interface.
The strength of such a dissipative force, that is suitable for wide ranges of interface width to characteristics flow length ratio to ensure no-slip boundary condition, is then termed as an optimum coupling parameter h * . Due to numerical simplicity of this approach, many researchers have employed it for the simulation studies [10][11][12], notwithstanding the fact that, a formal thin interface limit, for both of these coupling, has not been established. The goal of present work is to summarize our findings on the existence of a thin interface limit in such a case.
To keep the analysis tractable, anisotropies of the surface energy and kinetic coefficient are neglected. Diffusion coefficients and densities of the liquid and solid phases are assumed to be identical. Due to this equal density assumption, the melt velocity in a direction normal to the interface vanishes, thus simplifying the analysis. The growing solid is assumed to be stationary and does not move under the forces exerted by melt flow. Special attention is paid to ensure the no-slip boundary condition. We introduce the following notation: u is the reduced temperature field u = T −Tm L/Cp , T is temperature, T m is melting temperature, L/C p is the so-called hypercooling limit, δ is capillary length, β is kinetic coefficient, w is the melt velocity, ρ is density, µ l is dynamic viscosity of the melt, p is pressure, g is acceleration due to gravity. The phase field (ϕ) and reduced temperature field equations are, where f (ϕ) = −ϕ + ϕ 3 is the well-known double well potential corresponding to phase field with values ϕ = −1, +1 in the liquid and solid phases, respectively. g (ϕ) = (1 − ϕ 2 ) 2 is an interpolating function that is non-zero only inside the interface, A 1 is a numerical constant and τ is the relaxation time. For the melt flow, we first proceed with an improved version of the drag force model that ensures Galilean invariance of the melt flow equations [13]. With this choice, the Navier-Stokes equations read, where γ is a coefficient related to thermal expansion, and H(ϕ) is an interpolating polynomial with H(±1) = 0. The description of melt flow dynamics is completed with the continuity equation, The small parameter for the asymptotic expansion is identified as a ratio of interface width to diffusion length, ε = W V D = W L d . A curvilinear orthogonal system of coordinates that is attached to the moving interface, with unit vectorsr (normal to the interface) andŝ (tangential) is chosen to analyze the coupled set of equations. The scaled length in a direction normal to the interface, r ε , is denoted by η. The limit η → ±∞ corresponds to the liquid and solid side of the interface, respectively.
Melt flow is expanded for inner w (microscopic) and outerw (macroscopic) variables, up to second order in ε as w ≈ w 0 + εw 1 + ε 2 w 2 andw ≈w 0 + εw 1 + ε 2w 2 . A similar expansion is used for ϕ and u, where u n , ϕ n denotes order of approximation in ε for integer n. The macroscopic melt velocity can be Taylor expanded in the direction normal to the interface around the position of a hypothetical sharp interface at r = 0 as follows [14], In the present case, the no-slip boundary condition at the liquid-solid interface can be written asw − (0) = 0. The superscript − denotes the quantity evaluated at the interface when approached from the solid side of the phase.w(0) reminds us that these macroscopic quantities are evaluated at the sharp interface position, r = 0, which coincides with the center of the diffuse interface. From the no-slip boundary condition we conclude that ∂ kw− 0 (0) ∂r k = 0, for positive natural integer k. We denote the normal and tangential components of the melt velocity by w s and w r . We write the continuity and momentum balance equations for the melt flow dynamics as, where κ is interface curvature and g r , g s are normal and tangential components of the gravity. Noting that the densities of the melt and solid are the same, the variation of the normal component of the melt velocity w r across the interface is neglected. With this assumption, we proceed with analysis of equations (6) and (8) at successive orders of ε. At the second order, in combination with equation (5), we obtain ∂ 2ws− 2 . This means that the matching condition on the solid side of the interface is not satisfied in the second order of ε-expansion for the drag force model. In view of this result, we also examined the variable viscosity model [7]. For this choice,

450
The European Physical Journal Special Topics through a similar analysis, we obtain ∂w − 1 (0) ∂r = 0, ∂ 2w− 0 (0) ∂r 2 = 0. An extensive account of this approach and additional issues faced by the original drag force model [9] due to violation of the Galilean invariance can be found in reference [14]. We just here conclude that, the variable viscosity model can satisfy matching condition for inner and outer velocity fields, even in the presence of body forces like gravity. For both of these couplings the relation between phase field parameters τ = W 2 β δ + M D W δ , that was originally devised for the diffusive transport [4], remains valid. This relation is necessary to comply with macroscopic energy balance and Gibbs-Thomson relation at the interface. The constant M depends on the chosen forms of f (ϕ) and g (ϕ) in equation (1).

Numerical simulations
To test the above analysis, we perform numerical simulations of a 2D dendrite, growing in the direction opposite to an externally imposed melt flow. We include a fourfold surface energy anisotropy in equation (1) as described in reference [4]. The phase field equation in this cases is, where τ = τ 0 a(n) 2 , W = W 0 a(n), a(n) = 1+ cos(4θ) and θ is the angle between normal to the interface and some fixed direction. τ 0 and W 0 are the reference relaxation time and interface width. denotes the strength of surface energy anisotropy and a positive value of is a necessary condition to achieve a steady state. This construction effectively ensures the capillary length δ to be of the form δ = δ 0 (1 − 15 cos(4θ)). In addition to equation (9), the heat transport and melt flow dynamics are solved with equation (2), equations (3) and (4).
To compare the simulated growth velocity with the corresponding sharp-interface solution, we refer to a recently developed analytic Alexandrov-Galenko (AG) theory [15], which predicts, Here: b, a 1 , σ 0 are numerical constants, and R is the tip radius of resulting steady state parabola (|w ∞ | is the far field melt velocity). V g is the steady state growth velocity, Re = ρR|w∞| µ l is the Reynolds number and P g = V g R/(2D) is the growth Péclet number.
To predict the behaviour of dendritic growth and compare it with the AG theory, the constants of computations have been taken from reference [14]. Figure 1a shows a typical growth of a dendrite along with iso-temperature curves and velocity vector arrows surrounding the dendrite. Figure 1b compares the scaled velocityṼ g versus P g for the PF simulations and the AG theory, showing excellent agreement between the two. Simulations with higher flow velocities confirm this agreement further (results not shown). This agreement suggests that neglecting anisotropic terms in thin interface asymptotics does not alter the main conclusion regarding the independence of the simulation results on the interface thickness.

Conclusion
A thin-interface analysis of the phase field equations in the presence of melt convection is provided. It is shown that, as in the case of diffusive transport, the thickness of the diffuse interface can be chosen such that its effects on the obtained results are minimized. This prediction is verified by a comparison of the numerical simulation results for dendritic tip velocity and an analytic theory, which accounts for flow effects. As an outlook for further work, it shall be noted that, unlike the temperature or the solute fields, the melt velocity identically vanishes in the solid domain.
The non-vanishing normal gradients of the tangential velocity contribute to the shear stress tensor that generates an equal and opposite resulting force on the growing solid. The solid structure, when allowed to move, in turn, influences the melt flow field and thereby transport of heat and solute. In such a case, the unbalanced shear stresses on solid might play an important role.