An enhanced phase field model for micron-scale droplet impact with solidification microstructure formation

In this work, a novel diffuse interface model combining droplet impact with solidification microstructure formation was developed. A number of numerical models simulating droplet impact with solidification have been invented, but few are capable of unveiling the nucleation and growth of polycrystalline crystals on a micro scale. This paper thus aims to propose a diffuse interface model to simulate droplet impact, and moreover, solidification microstructure formation. To implicitly track the evolving liquid-gas interface, the Cahn-Hilliard equation is coupled with the Navier-Stokes equation. A phase field model involving polycrystalline growth is responsible for the capturing of solid-liquid interface and grain-grain boundaries. The current model is discretized explicitly such that it lends itself to shared-memory parallelism like OpenMP. A parallel SOR scheme based on Red/Black ordering is used to solve the pressure velocity coupling. The model was applied to various impact conditions in plasma spraying, including both single and multiple droplet impact, to show its robustness. Besides, with the model being modified slightly, non-equilibrium effect induced by rapid solidification was investigated as well.


Introduction
Droplet impact with solidification is common in industrial applications, such as three dimensional (3D) printing and plasma spraying. In 3D printing, metallic droplets, often millimeters in diameters, are ejected through a nozzle onto a cold surface, and then undergo solidification. The impacting speed is low and a part is manufactured layer by layer. Besides, to obtain better surface quality or uniform metal trace, it is necessary to optimize printing step and ejection/deposition temperature. In plasma spraying, on the contrary, the impacting speed is rather large, sometimes even supersonic and thus bringing rich fluid behaviors such as splashing. Droplet size appears much smaller than in 3D printing, usually in microns or submicrons. Another marked difference between 3D printing and plasma spraying lies in the solidification behavior. A droplet of polycrystalline materials under plasma spraying conditions is greatly supercooled before solidification, with solidification microstructure often taking on a columnar shape.
Numerical studies regarding droplet impact with solidification are extensive. There are two moving interfaces in this process: liquid-gas and solid-liquid interfaces. The moving liquid-gas interface is captured by volume of fluid method [1][2][3], lattice Boltzmann method [4,5], level set method [6,7], phase field method [8], and so on. While the method used in tracking solid-liquid interface in droplet impact simulations is confined to a couple of ones. Pasandideh-Fard et al [9], via the enthalpy transformation model by Cao et al [10], altered the heat balance equation into one depending only on the enthalpy, and rested on the value of enthalpy to distinguish solid and liquid phases. A variant of this method, which introduces the notion of liquid fraction, simply adds a source term to the original heat equation with temperature as the sole dependent variable, accommodating readily to well-developed commercial/private codes and thus embracing wider application. Such a method relies on the value of liquid fraction to differentiate solid and liquid phases and usually needs iterative solutions. With this method, Tembely et al [11] studied supercooled micro-droplet impact and freezing on superhydrophobic surfaces. Also, Xiong and Cheng [12] examined air entrapment in a molten and solidifying droplet colliding with a cold substrate. Shen et al [8] utilized such a source based method to simulate tin droplet impact with solidification under 3D printing conditions. Other variants such as apparent heat capacity method, effective heat capacity method, and heat integration method, can be found in [13]. They are all fixed grid methods and applied to both solid and liquid phases.
While the above mentioned methods are popular, they can deal only with macroscopic solidification behaviors. Microscopic solidification behaviors, such as crystal nucleation, grain impingement and grain coalescence, can by no means be restored by such methods. Nevertheless, phase field method can also serve the purpose [14]. By solving an extra governing equation for the order parameter, a phase indicator that plays the role of liquid fraction, the solid-liquid interface is implicitly tracked. The modified heat balance equation is, however, amenable to explicit finite difference methods. Moreover, the mechanic properties of 3D printed products or plasma sprayed coatings are directly related to their solidification microstructures. Therefore, introducing phase field method for solidification into droplet impact simulations is of practical importance, since it could be employed to predict how solidification microstructure is formed in such processes and what impacting parameters influence its formation and growth. Nevertheless, such research efforts are relatively few [15,16].
The paper is thus aimed to present an enhanced phase field model for droplet impact with solidification microstructure formation. The two interfaces are both implicitly tracked by the phase field method. It is able not only to predict droplet spreading, but also to monitor solidification microstructure formation. The paper is organized as follows. First, the mathematical statement is given, where detailed discretization and solution algorithm are provided. Second, the model is applied to several benchmark tests to show its robustness.

Mathematical statement
2.1. Governing equations for the two-phase flow An alternative to the traditional sharp interface model is a diffuse interface one, which assumes that the interface is of finite thickness and that physical quantities vary smoothly and continuously across the interface. The interface is advected by the Cahn-Hilliard equation [17].
where M is the mobility and d d = G f c / the chemical potential, withd dc (·) / being the variational derivative with respect to c.  is the gradient operator and  2 the Lapacian operator.
is the free energy density, where x measures interfacial thickness and g surface tension.
( ) ( ) / takes on a double well shape and represents specific free energy at equilibrium. The chemical potential in two dimensions can be obtained via the Euler equation [17], where Y¢ means seeking derivative with respect to c.
When the interface is at equilibrium in the absence of fluid flow, the diffusive flux vanishes and the interfacial profile can be found making = G 0. In one dimension, it takes the shape where y denotes the coordinate perpendicular to the interface. Besides, at equilibrium, the surface tension force can be interpreted as the integral of the excess free energy per unit surface area across the interface. (3) into equation (4) yields a = 6 2.In this paper = c 1 is to denote a gas and = c 0 a liquid, with the value in between indicating the diffuse interface.

Substitution of equation
The coupling of the Cahn-Hilliard equation with the Navier-Stokes equation for two immiscible and incompressible fluids with large density and viscosity ratios is a nontrivial thing. A model H, proposed by Hohenberg and Halperin [18], has been widely used for matched densities, and its thermodynamic consistency has been proved by Gurtin et al [19]. As for two fluids with large density ratio, a so-called modified H model was proposed by Jacqmin [20], who directly replaced the constant density r with r c , ( ) a function of the order parameter, while keeping the velocity field divergence free. This use of the phase field modeling renders thermophysical properties as functions of c (e.g., r r r r = -+ c c g l l ( ) ( ) where g denotes gas and l liquid). The resulted equation therefore applies to the whole two-phase flow region, taking the form of.
Therefore, mc is the volumetric counterpart of surface tension. It can be derived from different approaches, such as the least action principle [22] and visual work method [23]. A recent review can be found in [24].
To the authors' best knowledge, the thermodynamic consistency of equation (5) has not been proved. A benchmark test, however, conducted by Aland and Viogt [25] showed little difference between the modified H model and another thermodynamically consistent model by Abels et al [26], which is much more complicated. The difference dwindles as interface thickness decreases. Moreover, both models are divergence free.
Also, as for the coupling between the Cahn-Hilliard equation with heat energy balance equation, a number of researchers [27,28], have modified the density r and thermal conductivity k as functions of the order parameter c directly without further proof. This paper is not intended to give a full derivation of the heat balance equation coupled with the Cahn-Hilliard equation so that we simply adopt their approach, with the heat balance equation given by S T is a source term and c P heat capacity at constant pressure. To summarize, the two-phase flow with heat transfer, coupled with the Cahn-Hilliard equation, is governed by equations (1), (5), (11), forming a system of equations previously used by [29], and equation (12). These four equations, along with the phase field equation for solidification microstructure formation, are the focus of the paper.
To model solidification, a liquid fraction f l is introduced, which is related to the other order parameter f as f =f 1 2 .
l ( )/ As discussed below, a diffuse interface is assumed between solid and liquid, with f = -1 indicating liquid and f = 1 solid. The value in between stands for the diffuse interface, or the mushy region. Such a linear transformation will make f l approach zero when the content in a control volume is completely solidified and f l be unity when no solidification occurs at all. The latent heat of solidification will dissipate through the newly solidified portion and S T in equation (12) accounts for this latent heat, its form resembling that of the enthalpy of mixing [30] Note that the subscript l denotes liquid properties, and that L l flags latent heat per unit mass, and that f ¶ ¶t, / as introduced below, will evolve only where  c 0.5. f ¶ ¶t / is forced to vanish where > c 0.5. The current method treats the mushy zone as a porous medium. If the liquid fraction of a computational cell is equal to zero, the velocities within will also become zero. This is effected through the source term S m [31].

= --
where d is the mushy region constant to dampen the velocity in the solid, while b is a quite small number to avoid division by zero.

Governing equations for solidification microstructure formation
Kobayashi and coworkers [32][33][34] pioneered the orientation-based phase field model that can track solid-liquid interface and grain boundaries. The starting point is a free energy functional in terms of q, f, u, with q representing local orientation of crystals with respect to a fixed coordinate, f the solid-liquid order parameter, where T m is the melting point. Note that f is another order parameter, with f = -1 denoting a fluid and f = 1 a solid. The free energy functional is given by )is the bulk free energy density, and the second term accounts for the energy penalty within the interface, with e j f ( )being the anisotropic gradient energy coefficient. f h( ) is so chosen that it activates only in solid and vanishes in liquid. s and e q are coupling constants.
is the inclination angle of the interface. The kinetic equations for f and q are given by [35,36] where t f and t q are interface attachment kinetic times, l is another coupling constant, and e q m e  = + - If anisotropy matters, then e f and t f have fourth-fold symmetry. e 4 and e f are the strength of anisotropy and characteristic interface thickness, respectively. e j e e j q = + -

Boundary conditions
The boundary conditions are given in figure 1. Note that only haft computational domain enters calculation. When computing the thermal field, all boundaries are adiabatic; y contains all field variables, except for the normal or all velocity components at the axis of symmetry and at the wall, respectively, and for the order parameter c at substrate surface that satisfies where n is the unit normal pointing to the wall. Equation (17) may be derived as follows. The total free energy now consists of two parts, The prime means derivative with respect to the order parameter c. The first term on the RHS can be replaced using Green's first identity the Lapacian operator, applying the divergence theorem yields.
where again n is the unit normal pointing outwards to the wall. Note that (20) into equation (18) and rearrangement gives c is the bulk chemical potential and the surface chemical potential. To minimize the total free energy we make which also prescribes the wetting condition for the order parameter c. In this work on other boundaries, the latter of which leads to a homogenous Neumann condition by equation (17), or more precisely, a contact angle of  90 . A proper form of the wall free energy density f c w ( ) must comply with the following requirements. First, = f c 0 w ( )gives the liquid-solid interfacial tension g w1 and = f c 1 w ( )produces the gas-solid interfacial tension g , w2 the two determining the static contact angle q S through Young's equation g g g q -= cos .
Second, = f c 0.  series of ghost nodes are added to enforce boundary conditions. In the explicit method, the value of a field variable at the current time level, like f + , n 1 are evaluated with available values at the previous time level, here denoted with superscript n. Take the phase field equation, equation (15) and equation (16), for an instance. Their explicit discretization scheme takes the form. t f where S is a source term so that f + n 1 can be evaluated directly from the available values at the previous time level. The same strategy is applied to equations (12) and (1).
Note that f in equation (24), or q in equation (23), may use the value at the current or previous time level according to the programmer, and that such minor difference has little impact on final outcome since the time step satisfying numerical stability is rather small, as demonstrated in section 2.4.3. The discretization of convective and diffusion terms is described below, which is a little bit complicated. But before that, attention is paid to the incorporation of thermal contact resistance into the energy balance equation and to the discretization of the wetting boundary condition. Notice that the no slip condition is applied at row = j 1. The wetting boundary condition equation (17) can thus be discretized as / withn the unit normal pointing downward into the wall. Thermal contact resistance, accounted for in row = j 1 and = j 0, comes into the heat equation as the modified thermal conductivity k e in equation (28).
In addition, in discretizing the Poisson equation / where u* is the interim velocity, harmonic interpolation could be used to obtain a better and more consistent result. Take x direction for an example. For any node i j , ( ) The convective terms like y  u · is discretized using the first order upwind scheme, which is acceptable when the flow direction is parallel to the grid. Take x direction for an example.  (23) and (24). (23)    where the superscript + n 1 indicates the current time level. It is easy to verify that summing equations (37) and (38) restores equation (5). The pressure is sought through which follows directly from equation (38) after applying the divergence operator and remembering  = Due to its elliptic nature, equation (39) is solved iteratively using a parallel version of SOR as mentioned above. To sum up, the velocity is first updated using equation (37) without considering the pressure effect, and having solved equation (39) for the pressure at the current time level, correct the velocity via equation (38).

Numerical stability of the explicit scheme
For explicit methods to be stable, care must be exercised when choosing the time marching step. There are four major fields of concern in the current study: the phase field for solidification, the orientation field, the heat convection field, and the flow field. For the former two fields, the time step is set to satisfy where a is the thermal diffusivity of the larger one between YSZ and air, and u and v are defined as the maximum velocity that may be attained. Finally, the time step in the flow field is restricted by the Courant-Friedrichs-Lewy u max   is the norm of the largest possible velocity vector. Accordingly, the minimum of the three will be selected as the time marching step. Throughout the paper, the time step D~t 10 s. 13 The flow chart of the solution algorithm is displayed in figure 3.

Parallel programing using OpenMP
To accelerate computation, parallel programing based on OpenMP is applied. Remarks are made about the Poisson equation that is tackled using a parallel version of SOR as in figure 4(a). While the domain partition strategy in figure 4(b) still works and converges, the convergence rate is relatively low due to data race: when updating line 3, the value in line 2 may not be the latest.

Results and discussion
This section begins with a Nickel droplet impact onto a preheated surface, and then discusses the effect of the phase field mobility M on YSZ droplet spreading, followed by the sequential YSZ droplet impact and the investigation of rapid solidification. The thermophysical properties are listed in and that d 0 is the thermal capillary length and D thermal diffusivity. L c is a characteristic length, for instance, droplet diameter. m m m = e g l refers to the effective viscosity, with g denoting gas and l liquid. The model parameters are so chosen that the diffuse model would converge to the traditional sharp interface model. Throughout the paper, unless otherwise stated, D =´x 4 10 m, 8 which is fine enough for numerical outcome to be mesh independent [15]. Grid independence test and experimental validation can also be found in [15]. The model parameters in table 2, unless otherwise stated, is used throughout the paper.

Eliminating lattice effect
Before moving to case studies, attention is directed to eliminate lattice effect in the selection of crystalline orientations. When a control volume centered on node P in figure 5 turns solid, its crystalline orientation is randomly selected from those of its neighbors that are already solidified with definite crystalline orientations. This can be done in two ways, one staring from the node WS either counterclockwise or clockwise and the other from the nearest black nodes to the yellow ones, each stage being without order. More specific, the newly solidified node with its sum of coordinates being odd (even) proceeds counterclockwise (clockwise).
A simple test was conducted. A computational box with adiabatic boundaries is filled with supercooled YSZ at 2100 K. Initially a number of nuclei are seeded on the bottom. The grid is301 101, with a uniform spatial step of D = D =´x y 4 10 m.

10
, 2˜l = 14.53. The outcome is shown in figure 6. As demonstrated in figure 6(a), the approach in figure 5(a) can effectively reproduce grain impingement and competitive growth, while the method in figure 5(b) leads to unphysical grain boundaries that align themselves with grid axis in a nearly fixed angle.

A subsonic Ni droplet impact
Among the materials employed in plasma spraying, Ni and YSZ stands out due to their excellent thermophysical properties and improved thermodynamic compatibility with the oxide. The two, however, differ greatly in thermal conductivities by almost an order of magnitude. According to [15], solidification may be safely assumed to commence after droplet spreading, for YSZ droplet impact. In other words, convection effect is limited in determining solidification microstructure formation. But, this may not be the case for Ni droplet impact, for Ni has a higher thermal conductivity, which may enhance solidification by extracting more heat during early stages of spreading, thereby arresting fluid flow in advance. This section is thus aimed at investigating Ni droplet impact with solidification microstructure formation. The thermophysical properties and model parameters are listed in tables 1 and 2, respectively. But attention is paid to the choice of t q and e q , both of which are varied in the coming simulations to check their effect on orientation field evolution. A Ni droplet of diameter m 4.8 m impinges at -200 m s 1 onto a substrate of stainless steel that is preheated to 423 K. The droplet and surrounding is at 1800 K initially. The contact angle is set to  120 and interfacial heat transfer coefficient W 10 m K.
and consuming part of kinetic energy. It is not long before the spreading front ceases to move, as shown in the spreading factor in figure 8. The discontinuity in figure 8 is caused by a tiny daughter-droplet ejected at m 0.03 s. Subsequently, the floating liquid above the solidified, pushed by inertia force and rounded by surface tension, gradually settles down, taking a final shape as shown at m 0.14 s when the whole droplet is turned solid. Another matter of concern is the solidification microstructure. Careful observation in column (b) and (c) finds out that with the increase of t q and e q , grain boundary motion and grain rotation is greatly invigorated, with curvature smoothened and grain coalesced as shown in the last row of column (b) and (c).
To summarize, marked differences exist between Ni and YSZ droplet impact regarding the interaction between droplet spreading and solidification. Because of a higher thermal conductivity, it may no longer be reasonable to stick to the principle of spreading first and solidification second. Besides, though the precise value of t q is not available, it appears its choice doesn't affect droplet spreading significantly, for the grain boundary motion, as well as grain rotation, and droplet spreading are definitely not on the same time scale.
3.3. Effect of the phase field mobility on supersonic YSZ droplet impact In this section, let us return to YSZ and check one of the subtleties in the Cahn-Hilliard equation, the phase field mobility M, which realizes contact line motion in the presence of the no slip condition. Empirically, it should be large enough for diffusion to maintain a more or less constant interfacial thickness, but small enough to avoid excessive damping [23]. Consequently, this parameter appears to be material independent from this point of view, since it could be adjusted artificially. But it is actually a material dependent parameter [37], although not amendable to experimental measurement. To make the Cahn-Hilliard phase field method, which originates as a  phenomenological method, of quantitative importance, many research efforts have been directed to the study of its sharp interface limit. There exists two dominating standing points on how to realize it. Yue et al [38] thought the sharp interface limit could be achieved for a fixed mobility M in the limit of vanished interface thickness x. This, however, contradicts some other views [37,39] that the mobility M, although material dependent, should be adjusted according to a power scaling law due to an arbitrary interface thickness x. The debate, however, has not ended. But it is still of practical significance to examine the effect of the phase field mobility M on supersonic impact of a molten YSZ droplet under plasma spraying.
The numerical configurations are as follows. A YSZ droplet, m 4.8 m in diameter, impinges at -518 m s 1 onto a solid surface that is preheated to 423 K. The initial temperature is 3000 K and contact angle  120 . Interfacial heat transfer coefficient is set to 10 W m K. = Cn 0.0052, l = 72.68. Figure 9 shows the snapshot sequences of the numerical results. Note that the mobility is chosen following an empirical formula proposed by Yue et al [38]  / to approach the sharp interface limit. Its effect on splat morphology is evident, but not so on solidification microstructure formation as seen in figure 10. The whole solidification time is nearly the same, but there are clear satellite droplets ejected in the case with smaller phase field mobility. Besides, the spread factor appears to be larger. This may be induced by contact line pining, that is, immobilized contact lines along the substrate surface, due to a smaller diffusion rate. Although the contact line is inactive, the fluid nearby has sufficient kinetic energy, sharpening the spreading front while pushing it forward. Finally, the satellite

Sequential YSZ droplet impact
In this section, sequential YSZ droplet impact is modelled to demonstrate the robustness of the current model. The second droplet is introduced when the first is solidified. The mesh size is751 351. = Cn 0.0104, l = 145.36, and =´--M 2.2 10 m s kg . 13 3 1 · Note that the impacting velocities are -518 m s 1 for both droplets, but that the anisotropic strength e 4 has been reduced to 0.01 to check its effect. Besides, the nuclei density is increased compared with that in [16], where other numerical configurations can be found. The numerical outcome is shown below in figure 11.
The second droplet sets in at m 0.4 s, when the first has been frozen completely. The second droplet, driven by a radial pressure gradient, spreads on the first. However, owing to the uneven surface of the first splat as at m 0.4 s, a small droplet splashes at m 0.6 s. Surface tension curls up the droplet as shown at m 0.7 s, and then splat shape undergoes few changes. It takes around m 1 s for the second drop to solidify. Figure 12 gives solidification microstructure evolution. Of interest is the grain structure in the overlapping region where the grain boundary has been heavily skewed, taking a zig-zag shape as shown at 1.4 ms. The grain boundaries on the bottom bends rightwards, which is most likely caused by drop spreading, while those on the top bends leftwards, which may result from the curling up of the interface induced by surface tension.

·
The red represents solid, the blue liquid, and the green gas.  where + n 1 denotes the current time level while n andn 1 signifies one and two time steps back, respectively. At the beginning, f n is given and f + n 1 is the one to be sought. The numerical configurations are the same as in section 3.2, except =´-M m s 4 10 kg 13 3 · / and the initial temperature 3000 K. t is varied from´-0 4 10 s. 15 The results are shown in figure 13.