DIRECT NUMERICAL SIMULATION STUDY OF WATER DROPLETS FREEZING ON A HORIZONTAL PLATE

There are many experiments on solidification (i.e. freezing) of water droplets on the cold plate but the numerical investigation of its freezing process is very little. So, to provide an understanding of that, we present a numerical method and its results to describe the solidification process of water droplets on a cold plate at different wetting angles (o) in the range of 30 155. The Navier-Stokes and energy equations are used and solved by an axisymmetric fronttracking/finite difference technique. The distinct phases are separated by the interfaces represented by connected elements. So, there is a three-junction point among three phases: solid, liquid, and gas. The water droplets are assumed as a spherical cap and placed on the cold plate which is kept at a subfreezing temperature. At the end of solidification, we obtain a small protrusion shape at the top of the frozen water droplet and its height is also higher than that of the initial water droplet. That is explained by the increase in volume upon freezing because of the difference in densities of water and ice. The frozen water droplets are also compared with the corresponding experimental ones reported in the literature.


INTRODUCTION
The solidification of water droplets can find in nature and industrial processes such as the freezing water droplets on wind turbines, aircraft, air conditioning, fridge and any cooling types of equipment [1 -4]. It may cause serious hazards as reducing the performance of machines or decreasing their lifetime. Providing energy to remove the ice on the surface is one of the good ideas but it seems to be expensive and influences the performance of machines, too. Using chemical compounds is another choice but it may cause bad effects on the environment. Therefore, deep understandings of the freezing process would be a significant task for researchers to find the solution to avoid the damage as well as enhance the efficiency of machines. The spreading and solidification process of water drop was investigated at different contact angles (i.e. wetting angels) (o) by Huang et al. [5] as well as by Pan et al. [6] in the range of 76 o -155 o and 77 o -145 o , respectively. In the same work [6], Pan et al. also studied the solidification of water droplet on an inclined surface (φ = 30 o ). The freezing process of a water droplet was monitored on different cold concave surfaces with radii (R) in the range of 10 mm -25 mm and R = ꚙ [7] as well as on spherical surfaces with radii (R) ranging from 15 mm to 30 mm. In this study, we only consider the solidification of a water droplet on a cold flat surface.
Numerically, Shetabivash et al. [9] presented a multi-level-set approach to simulate the solidification process of water droplet on a flat surface. In another research, Vahab et al. [10] used the level-set approach and the moment of fluid method to consider the solidification of a water droplet on a substrate in order to apply to aircraft icing. Ajaev and Davis [11] applied the boundary integral method [12] to simulate the solidification process of droplet. However, these researches did not consider the solidification of a water droplet at different wetting angles (o). To fill this gap, we represent a front-tracking method [13] to numerically simulate the freezing process of water droplets on a horizontal plate.

NUMERICAL MODEL AND METHOD
The problem? is described in . The shape of the water droplet is assumed to have a spherical cap. Because of the symmetry of the droplet, we only simulated a half of the droplet which is placed on a cold plate (To) in a computational domain W×H. The interfaces separate three phases: solid, liquid, and gas. Thus, there are three interfaces, namely solid-liquid, solid-gas and liquid-gas. Accordingly, there is the presence of a three-junction point. We use the Navier-Stokes and energy equations to solve the problem, and the front-tracking method [14] is used to represent the phase boundaries. The liquid and gas in this study are supposed to be incompressible, immiscible and Newtonian. In each phase, density ρ, viscosity µ, heat capacity Cp and thermal conductivity k are constant, the Navier-Stokes and continuity equations are: .
here, u is the velocity vector, p is the pressure, σ and κ are the interfacial tension coefficient and twice mean curvature, respectively. The Dirac delta function is denoted by δ, x is the position vector, subscript f expresses the interface, n is the unit normal vector, force f is used to apply the non-slip condition on the solid interface [15 -18], g is the gravitational acceleration, Lh is the latent heat of fusion, q is thermal flux at the solid-liquid interface, subscripts s and l represent solid and liquid, respectively.
In this study, the front-tracking technique with the finite difference scheme, previously presented in some works [13,19], is applied, where the interfacesinclude the points xf connected to each other, and each point moves with the velocity Vn which is interpolated from the nearest grid points: with the time evolution, t and t + t correspond to steps n and n + 1. Thanks to the points xf, we can compute the interfacial tension forces acting on the interfaces and hence build the indicator functions to identify the properties of the phase. To determine the distinct phases, we use two indicator functions I1 and I2 reconstructed from the positions of the interfaces. The value of each indicator function is 1 in a fluid and 0 in another one. The subscripts s, l, and g are for solid, liquid, and gas, respectively. Thus, the properties of the phases at every position in the computational domain are given as [20,21]: More details about our method can be found in some previous works [13,19]. The density difference between ice and water can be seen in Figure 2b, where the density in the region under the freezing interface is smaller than that above the freezing interface, which reflects a property of water (the density of ice is smaller than that of water). Besides, a velocity field is presented in Figure 2b. Over time, at  = 0.72 (Figure 2c), the region of the solid state of the water droplet is increased, this is in contrast to the decrease of the region of the fluid state of the water droplet. Similarly, the distinct regions of the water droplet, density field, velocity field can be seen in Figure 2c. Finally, at  = 1.243 (Figure 2d), the solidification process of the water droplet is almost finished. The protrusion of the frozen water droplet's top appeared; this feature can be explained by the volume expansion. We compare our study with the experiment of Pan et al. [6] at a wetting angle o = 78 o . In Figure 2d, our study is in the right frame and the experiment of Pan et al. [6] is in the left frame. As reported by Pan et al. [6], the water droplet was initially a part of a sphere (see the inset in the top-left corner of the left frame in Figure 2d) because of surface tension acting on the water-gas interface, and it became a frozen water droplet with an apex at the droplet top as shown in the left frame ( Figure 2d). As can be seen in Figure 2d, our numerical model agrees well with the experiment of Pan et al. [6].  Finally, at  = 2.74 (Figure 3d), the solidification process of the water droplet is almost complete. Similar to the previous case, we can see a shape like a horn on top of the water droplet. This feature is induced by the volume expansion. This case is compared with the experiment of Huang et al. [5] at a wetting angle o =124 o . In Figure 3d, our study is in the right frame and the experiment of Huang et al. [5] is in the left one. As reported by Huang and coworkers [5], at the beginning of the freezing process, the water droplet has a spherical cap (see the inset in the top-left corner of the left frame in Figure 3d) due to the surface tension force. The droplet became a frozen one with an apex at the droplet top at the end of freezing, as shown in the left frame (Figure 3d). Figure 3d confirms that the computed frozen droplet and that reported in Huang et al. [5] match very well.  Comparison between the numerical simulation (right) and the experiment Huang et al [5] (left). In (d), the inset in the top-left corner shows the initial water droplet on a cold plate in the experiment.

Water droplet at a wetting angle o = 155 o
interface of the droplet moves from the bottom to the top. Along with the movement of the freezing interface, from the bottom of the drop to its top, the normalized temperature has got low to high temperatures, the solidification happens at the location where it's temperature below the freezing temperature of water. Then, at  = 4.0 (Figure 4c), the normalized temperature expands from the bottom to the top of the domain, more and more area of the droplet is frozen. Hence, the liquid region of the droplet is decreased, this is in contrast to the increase in the solid region when time proceeds. Finally, at  = 5.79, as we can see from Figure 4d, the solidification of the water droplet is almost finished. A shape of the frozen droplet with a horn is formed. Similar to the previous case, we compare this case (at a wetting angle o = 155 o ) to the experiment of Huang et al. [5]. In Figure 4d, our study is shown in the right frame and the left one is the experiment of Huang et al. [5]. Like the previous cases (Figure 2d and Figure 3d), the water droplet with a spherical cap at the beginning of the freezing process (see the inset in the top-left corner of the left frame of Figure 4d) [5] became a frozen one with a conical top surface at the end of freezing (the left frame of Figure 4d). There is not much difference between the simulated and experimental droplets as shown in Figure 4d.

Variations of freezing water droplets on a cold plate at different wetting angles
Experimentally, Satunkin [22] studied the solidification of molten silicon and germanium droplets at a wetting angle of about 33 o . For water droplets, Huang et al. [5] conducted experiments on their solidification at wetting angles in the range of 76 o -154.9 o . In the work of Pan et al. [6], the wetting angle of water droplets was varied in the range of 77 o -145 o . It is found that the wetting angle of a liquid droplet on a cold plate studied in the aforementioned works is in the range of 30 o -155 o . To our knowledge, no further investigation on a water droplet at a wetting angle out of this range has been conducted so far. Accordingly, basing on our literature survey, we vary the wetting angle in the range of 30 o -155 o to investigate the solidification process of water droplets.  Figure 5a shows the shapes of the droplets after the complete solidification process, where the height of the ice droplets increases with increasing the wetting angle. We can see the increasing heights of the droplets in Figure 5b. After finishing the freezing process, the evolution height, h -ho, of the droplet increases with increasing the wetting angle in the range of 30 o -90 o . In contrast, further increasing the wetting angle from 90 o to 155 o causes the evolution height of the droplet to be decreased. Interestingly, the evolution heights of the droplets at wetting angles in the range of 90 o -155 o all have got the negative values during the initial stage of the freezing process. This can be explained by the effect of the gravitational acceleration. As mentioned, the initial droplets were assumed to be spherical and have an initial volume Vo. The droplets having a great height and a small area of contact with the cold plate are more affected by gravity. Therefore, the height of the water droplet is pulled down during the initial stage of the freezing process, resulting in the negative value of hho as shown in Figure  5b. Figure 5c presents the temporal variation of the ratio, in terms of %, of the volume V of the droplet to its initial volume Vo. It is clear that the volume of the droplet increases as the freezing process proceeds. After complete solidification process, the volume rises about 10 % as compared to the initial volume because of the volume expansion of the water drop upon freezing. In addition, the solidification time increases with an increase in the wetting angle in the range of 30 o -155 o . This can be explained by the contact area of the droplet with the cold plate. The larger the contact area of the droplet, the faster the solidification. Figure 5d shows the variation with respect to time of the mean value of the freezing interface height of the droplet. The solidification process finished with the maximum average interface height at a wetting angle o = 155 o . It indicates that decreasing the wetting angle will enhance the freezing process.

CONCLUSIONS
The numerical model was presented to simulate the freezing of water droplets on a cold plate under the influence of the wetting angles o. We consider the temporal evolution of the freezing process and compare with the experiments of Pan et al. [6] (o = 78 o ) and Huang et al.
[5] (o = 124 o and o = 155 o ). The results of the numerical simulations agree well with the experiments. The shape of the droplet after complete solidification process, the evolution height of the droplet, the percentage of the droplet volume compared to its initial volume, and the average freezing interface height were presented. After ending the solidification process, we can see the top of the droplets like a horn and that can be explained by the volume expansion of the water droplet upon freezing.
Declaration of competing interest. The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.