Mathematical Modeling of the Solid–Liquid Interface Propagation by the Boundary Integral Method with Nonlinear Liquidus Equation and Atomic Kinetics

: In this paper, we derive the boundary integral equation (BIE), a single integrodifferential equation governing the evolutionary behavior of the interface function, paying special attention to the nonlinear liquidus equation and atomic kinetics. As a result, the BIE is found for a thermodiffusion problem of binary melt crystallization with convection. Analyzing this equation coupled with the selection criterion for a stationary dendritic growth in the form of a parabolic cylinder, we show that nonlinear effects stemming from the liquidus equation and atomic kinetics play a decisive role. Namely, the dendrite tip velocity and diameter, respectively, become greater and lower with the increasing deviation of the liquidus equation from a linear form. In addition, the dendrite tip velocity can substantially change with variations in the power exponent of the atomic kinetics. In general, the theory under consideration describes the evolution of a curvilinear crystallization front, as well as the growth of solid phase perturbations and patterns in undercooled binary melts at local equilibrium conditions (for low and moderate Péclet numbers). In addition, our theory, combined with the unsteady selection criterion, determines the non-stationary growth rate of dendritic crystals and the diameter of their vertices.


Introduction
The dynamics of curved phase transition fronts moving in undercooled binary mixtures completely characterizes the birth of material structures that appear as a result of the phase transformation process. This dynamic behavior of such fronts under the influence of external forces also defines the evolution of morphological perturbations, the appearance and growth of interfacial ridges, and the formation of a two phase layer that divides the material phases [1][2][3][4][5][6]. An important point is that the curved phase interface is usually considered to be mesoscopically or macroscopically acute.
Let us introduce the interface function z interface = ζ(x, t) characterizing the position and velocity of a curvilinear phase transition front (see also Figure 1). This function satisfies the boundary-integral equation deduced in [7,8] for purely thermal problem and then extended in [9,10] for some solid-solid and fluid-fluid systems. Here, x = (x, y) and t stand for the vector of spatial coordinates and the time variable. Note that this method enables obtaining the boundary-integral equation for the interface function when considering the phase interface propagation in purely thermal or concentration problems in the absence of convection. Considering the thermal and concentration fields synchronously (an undercooled binary mixture), the authors of [11][12][13] extended the aforementioned approach for slow, moderate, and fast crystallization phenomena. On the other hand, the purely thermal theory [7,8] has been extended to the case of interface propagation in one-component and binary undercooled melts in the presence of a forced convective flow [14,15]. In this paper, the boundary integral theory is developed for the case of the nonlinear liquidus equation containing the kinetic contribution, which plays an important role with the increasing phase interface velocity. Strictly speaking, a new boundary-integral equation (BIE) for the interface function that describes the propagation of a curvilinear front into a binary undercooled mixture is derived. This BIE also describes the evolution of morphological perturbations of the solid-liquid phase interface arising as a result of temperature fluctuations or some external reasons [16]. What is more, the BIE under question characterizes how solid phase ridges or patterns can evolve into an undercooled melt and form complex material structures.

The Heat and Mass Transfer Model
The heat transfer equation for the liquid-solid (bulk) phases and the solute diffusion equation in the liquid phase can be written as where T d is the melt temperature, C d is the solute concentration, u is the melt velocity, D T is the thermal conductivity, D C is the diffusion coefficient, and V is a constant velocity of the Cartesian coordinate system x, y, and z (steady-state growth rate). The temperature, solute concentration, and fluid velocity are fixed far from the phase transition interface, T d = T d∞ , C d = C d∞ , u = −U ∞ , and the boundary conditions at the solid-liquid interface take the form of where d c = d 0 [1 − α d cos(ιθ)] is the anisotropy capillary length, d 0 is the capillary constant, α d 1 is the stiffness parameter, θ is the angle expressing the orientation between the normal and growth direction, ι is the order of the crystalline symmetry, K d is the average curvature of the interface, Q is the latent heat of crystallization, m 1 is the liquidus line slope, m 2 is the deviation of the liquidus equation from a linear form [17], c p is the heat capacity, β k is the kinetic coefficient, n is the power exponent of the atomic kinetics, T d0 is the phase transition temperature for the flat front, ds is the surface area vector element directed toward the liquid phase, k 0 is the solute partition coefficient, and ζ is the interface function. The average curvature of the interface K d in the two-dimensional case can be determined as In this section, we derive an integrodifferential equation for the interface function ζ(x, t) in the case of the thermo-diffusion problem with convection (1)-(4). For convenience, let us introduce the characteristic length 2D T /V and time 4D T /V 2 scales, as well as dimensionless temperature T = T d c p /Q, solute concentration C = m 1 c p C d /Q, and fluid velocity v = u/U ∞ .
Thus, the dimensionless model of convective heat and mass transfer following from (1) and (2) has the form: where v is the dimensionless velocity. The boundary conditions read as ∇T solid − ∇T liquid · ds = 2 + ∂ζ ∂t where β = β k c p V n /(2 n Q) and subscript i corresponds to the interface.

The Green's Function Technique
Let us use Green's functions G T (p|p 1 ) and G C (p|p 1 ) for the heat transfer and diffusion equations without flow to solve the convective problem (5)-(10): (11) where p = (x, z, t) and p 1 = (x 1 , z 1 , t 1 ).
Green's functions for heat G T (p|p 1 ) and chemical G C (p|p 1 ) problems can be obtained using the Fourier transform (see, among others, [9][10][11]18]) and read as By multiplying Equations (5) and (6) at p 1 by G T (p|p 1 ) and G C (p, p 1 ) and Equation (11) by T(p 1 ) and C(p 1 ), respectively, subtracting one from the other and integrating the result over p 1 , we obtain the temperature and solute concentration in liquid: where Λ 1 contains the point (x, z) and ε represents a small positive parameter. The integrals of the time derivatives in Equations (14) and (15) turn to zero after integration over time t 1 . In this case, Green's functions also tend to zero, and the temperature and concentration are bounded at the lower limit of integration. Moreover, the Green function is equal to zero at the upper limit because of the causality condition [18]. In addition, the integrals of the z 1 -derivative over z 1 in Equations (14) and (15) give a constant temperature T ∞ and concentration C ∞ . Next, applying Green's theorem, substituting the heat and mass boundary conditions (8) and (9) into (14) and (15), extending the point p to the phase transition surface, and taking ε → 0, we arrive at Here, the superscript ζ defines the phase interface. Let us write down the Gibbs-Thomson condition for the isothermal growth with the concentration gradient in the form of The solution to this equation relative to C i yields two real roots. One of them determining the interfacial concentration is positive and is given by Thus, Equations (17) and (19) together determine the growth of an isothermal dendrite from solution. For the combined thermal and concentration problem, let us now substitute T i and C i expressed in Equations (16) and (17) into the Gibbs-Thomson Equation (7). The result reads as where C i is defined by Equation (17). Thus, Equation (20) represents an integrodifferential equation for the interface function ζ(x, t) satisfying the convective problem (5)-(10) and describing the phase interface dynamics for purely thermal and thermo-concentration problems with convection with allowance for nonlinear liquidus equation and atomic kinetics. Let us especially note that the second term on the left-hand side of Equation (20) takes the flow into account.

Stationary Growth
In the case of stationary growth, the BIE transforms to Here, ∆ = T 0 − T ∞ is the dimensionless undercooling and the solute concentration in solid is constant. In the case of stationary growth, the only time-dependent function is Green's function, and we can integrate it over time for the thermal problem: where K 0 is the modified Bessel function. For the concentration problem, the integral over time can be written as Now, we obtain a general solution for the growth of any stationary form by substituting the thermal and concentration integrals of Green's function: where

The Parabolic Cylinder Reference Frame
Let us introduce the two-dimensional parabolic coordinates [14,19]: where p T = ρV/(2D T ) is the Péclet number and ρ is the crystal tip diameter. The integrals of Green's function were previously calculated in [12] as Next, we use the Jacobian J and replace the variables z 1 and x 1 in the convective integral of Equation (21) with the variables w 1 and s 1 defined by Equation (25): Therefore, replacing the variables in (21), we obtain Note that the convective integral contribution stemming from the solute concentration has the same form. The next step in our calculations is to obtain the temperature gradient. To do this, let us rewrite the heat conduction equation in parabolic coordinates as Integrating Equation (28) and taking the heat balance boundary condition into account, we express the convective term as The hydrodynamic field was found in [20,21] and reads as where Re = ρU ∞ /ν is the Reynolds number and ν is the kinematic viscosity. Now, substituting (29) into (27), we obtain The last step is to integrate the Green's function derived in [15]. The result reads as for the thermal and for the concentration contributions. Finally, we obtain the convective boundary integral equation for a dendrite having the form of a parabolic cylinder and growing with a constant velocity in an undercooled binary melt: As a special note, Equation (34) contains the fluid flow in terms of ν 2 , which is dependent on the Reynolds number Re(U ∞ ) (see Equation (30)).

Numerical Examples
The selection criterion σ * = 2d 0 D T /(ρ 2 V) for the 2D thermo-solutal dendritic growth with convection and surface tension anisotropy is defined by [21] where σ 0 is the selection constant, b is the selection parameter responsible for fluid flow, a 1 = (8σ 0 /7) 1/2 (3/56) 3/8 , a 2 = √ 2a 1 , α d = 15ε c , and ε c is the surface energy anisotropy. The selection criterion (35) describes the stable growth mode of a dendritic crystal with a constant velocity in an undercooled binary melt with a forced flow. As this takes place, the first term in square brackets represents the thermal contribution; the second term is the contribution from the solute concentration; the term containing coefficient b is responsible for the incoming flow; the terms containing the Péclet numbers p T describe crystallization with an increased velocity. In general, Formula (35) works for low and moderate Péclet numbers in local-equilibrium conditions. When dealing with the local-nonequilibrium scenario that occurs for high crystal growth velocities (Péclet numbers), the corresponding high-speed selection criterion should be used [13,22]. Figure 2 shows the influence of the nonlinear liquidus slope (coefficient m 2 ) on the crystal velocity-melt undercooling and tip diameter melt undercooling curves. As would be expected, the velocity V builds up, and the diameter ρ reduces with increasing the melt undercooling. As this takes place, V and ρ, respectively, become greater and lower with increasing the deviation m 2 of the liquidus equation from a linear form. Therefore, for example, choosing the melt undercooling ∆ = 75 K (right panel in Figure 2), we see that the dendrite tip diameter changes within the range of 20% when the tip radius varies between ∼4 × 10 −6 m (m 2 = 8 K (at.%) −2 ) and ∼5 × 10 −6 m (m 2 = 0). This significantly changes the alloy microstructure (e.g., the primary interdendritic spacing [23][24][25]). The influence of the power exponent of the atomic kinetics n on V(∆) and ρ(∆)/2 is illustrated in Figure 3. As is easily seen, the crystal growth velocity V can change up to 60% (compare the red and green curves at ∆ = 200 K). This means a great impact of the atomic kinetic law on the crystallization rate.
In general, the BIE developed in this paper describes a broader class of undercooled crystallizing binary systems with a moving interphase boundary than previously known theories. This is explained by taking into account the nonlinear dependence of the phase transition temperature on the dissolved impurity concentration and the power law for the kinetics of particle attachment to the solidification front.

Conclusions
In summary, in this study, we developed the boundary integral theory with allowance for the nonlinear liquidus equation with a kinetic contribution playing a significant role with the increasing phase interface velocity. The boundary integral equation, a single integrodifferential equation for the interface function, was derived for the thermo-concentration problem, which describes a solidification process of an undercooled binary melt in the presence of convection. This equation characterizes the motion of the curved solid-liquid interface (solidification front) propagating into an undercooled liquid phase in cases of small and moderate crystallization velocities (Péclet numbers) [26]. Our calculations of the boundary integral based on a parabolic tip shape of the dendritic crystal showed that a nonlinear contribution to the liquidus equation (the term proportional to the square of the solute concentration) and the power exponent of the atomic kinetics can essentially change the interface velocity-melt undercooling curve for moderate values of ∆T ∼ 50-250 K. With the further increase of the undercooling, the developed theory has to be generalized to the case of the hyperbolic impurity diffusion equation [27,28].
As a special note, the boundary integral equation under consideration has great potential applications for the further development of crystallization theory. For example, it can be used for simulating the pattern formation and crystal growth in undercooled liquids [12,29], studying the development of morphological instability [2,9,10,30], and deriving a selection criterion for a stable mode of dendritic growth [1,21,30]. An important task of crystallization theory is to derive the boundary integral while taking into account the nucleation and growth of particles in an undercooled liquid ahead of the phase transformation boundary. This problem can be solved by combining the present analysis and the theory of crystal growth at the intermediate stage of phase transformation [31][32][33][34][35][36][37][38]. Another important task is to determine the boundary integral equation for the directional crystallization process with a two-phase region. Such a problem can be attempted by combining the theory under consideration and the crystallization theory with a quasi-equilibrium two phase region [39][40][41][42][43][44]. The method of the boundary integral equation developed in this paper can also be applied to describe moving boundaries in biologically active media [45,46], as well as in complex heterogeneous systems of different physical natures [47,48].  Data Availability Statement: All data generated or analyzed during this study are included in this published article.

Conflicts of Interest:
The authors declare no conflict of interest.