Modelling the dynamics of the flow within freezing water droplets

The flow within freezing water droplets is here numerically modelled assuming fixed shape throughout freezing. Three droplets are studied with equal volume but different contact angles and two cases are considered, one including internal natural convection and one where it is excluded, i.e. a case where the effects of density differences is not considered. The shape of the freezing front is similar to experimental observations in the literature and the freezing time is well predicted for colder substrate temperatures. The latter is found to be clearly dependent on the plate temperature and contact angle. Including density differences has only a minor influence on the freezing time, but it has a considerable effect on the dynamics of the internal flow. To exemplify, in the vicinity of the density maximum for water (4 ∘C) the velocities are about 100 times higher when internal natural convection is considered for as compared to when it is not.


Introduction
Problems associated with the build-up of ice, e.g. on wind turbines, airplane wings, and roads, calls for a better understanding of the freezing process for water droplets on cold surfaces. Previous research has identified a number of factors important to the freezing process such as the temperature of the cooling surface, T p [1], the size of the droplet [2] the impact of free and forced convection [3,4], the roughness and wettability of a surface [5], the freezing on superhydrophobic surfaces, e.g. [6], the effect of an inclined surface [7][8][9], internal heat transfer [10,11] and internal flow [12]. A number of studies have also experimentally shown that the contact angle, θ, has a strong influence on the water droplet freezing time, t f [13][14][15].
Hao et al. [5] found that t f was dependent on the size of the contact area, A and the thermal conductivity of the surface. A smaller A, i.e. a larger θ , resulted in a longer t f . The overall knowledge of the freezing process is still fragmented and there is a need for more research in this area, especially regarding numerical models of the freezing process. Whilst a wide range of studies has used experimental methods to study the freezing process, e.g. [16][17][18] only a few have been using numerical methods. Numerical models has been used to study the heat transfer and phase change in the droplet. To exemplify, Kavanami et al. [12] used a model considering both surface tension and the density maximum at 4 • C and derived that these mechanisms are important for the internal flow whilst Chaudhary and Li [19] modelled the freezing process of water droplets on surfaces with different wettability. Numerical models has also been proposed for the geometrical phenomena that occur in the last stages of the freezing process, when a pointy shape is appearing. Anderson [21][22][23][24] proposed numerical models to predict the angle of the conical tip and capture the volume expansion of the droplet. In this work a droplet with fixed shape is used instead of a moving boundary (e.g. using a Volume of Fluid approach (VOF) where the free surface is tracked and located). The benefit of this approach is a simpler model where the focus is only on the transport of heat inside the droplet. The main concern is to see whether it is possible to capture the main features of the freezing process even though the model is not covering all aspects of the process, such as variation in θ. To fulfil the aim, a numerical model of the freezing process is created and the internal flow is studied for θ = 77, 84 and 90 • , whilst the volume is kept constant. Two cases are considered, one where internal natural convection is included in the model (case 2), and one where it is not (case 1). The experimental work done by Jin et al. [4] is used to validate the numerical model.

Theory
The continuity equation, energy and momentum equations are given by ∂ρ ∂t The mixture viscosity is given as μ = x l μ l + x s μ s . Note that the specific heat, specific volume and the conductivity also follow the aforementioned correlation. The viscosity for the solid is introduced to handle the zero velocities in the solid zone. The mixture enthalpy is calculated as ΔH mix = c p(mix) ΔT . The liquid and solid mass fractions, x l and x s follows the relationship x l + x s = 1. At the water/ice temperature T > 273 K, x l = 1 and at T < 273 K, x l = 0.
To determine the mass fraction of water at T = 273 K, the Lever rule is used The density is given by ρ(T ) = f (T ), where f (T ) is provided by Andersland et al. [25]. H is the static enthalpy given as [26] The latent heat released in the freezing process is modelled as the difference in enthalpy between water and ice. Here, the reference enthalpy for ice is 0 and the total reference enthalpy is then given by  [27,28]. The source term in Eq. 3 is given by, where the first term is used to regulate the velocity in the droplet as it freezes, meaning that K p is small in the water and large in the ice, and the second term is used to model the buoyancy effect with ρ ref = 990 kg/m 3 and β = 67.34 μK −1 . For case 1, i.e. excluding internal natural convection, the density of the water is kept constant and for case 2, i.e. including internal natural convection, the density of the water is allowed to vary with temperature. Please note that for both cases, since ice has another density than water there will always be a density difference between the two materials (about 8%).

Method
The simulations are set-up and carried out with the hybrid Finite Volume/Finite Element CFD software ANSYS CFX 15 in a similar manner as done in Ljung et al. [29][30][31] for evaporating droplets.

Geometry and grid generation
The geometries investigated in the simulations are based on the experiments by Jin et al. [4], where the droplets have an initial volume of V d = 9.32 μL. In the simulations, the droplets are kept at this volume throughout the freezing process to ensure that the correct amount of water is being cooled. For sufficiently small droplets, the effect of gravity can be neglected and therefore the droplets can be assumed to be half spherical. Then, a relation between droplet height (the height of the apex) a, surface radius r and θ exist such that, Experimentally it has been shown that the surface contact area, A, does not change with time during the freezing process [32]. The total V d will in its turn increase during freezing due the specific volume of the ice, leading to a increased θ and a. To investigate the influence of a, r and θ on the freezing process, three representative θ are approximated from the experiments by Jin et al. [4]. For the three times, t = 0, 6 and 12 s, a ≈ 0.80r, 0.86r and r are retrieved. The corresponding θ are 77, 84 and 90 • , see Table 1 where a, r and θ are calculated based on the attained ratio following the expression, Please note that the value of r attained in the experiments is displayed at t = 0 in Table 1. Both the geometries and the meshes are created using the software ICEM. Due to symmetry reasons, the droplets are modelled as a slice of a hemisphere, which is one element thick and extruded one degree. The meshes constructed are unstructured hexahedral meshes using the O-grid method. The geometry and mesh for the θ = 90 • droplet can be found in Fig. 1.

Simulation settings
The droplets are cooled from below at T p = −11.2, −8.2 and −5.3 • C, in accordance with the work by Jin et al. [4]. At t = 0 s there is only water in the model at a temperature of T w = 21 • C, which has no initial velocity and zero pressure. An adiabatic boundary condition in applied at the water-air interface. The model does not account for a subcooled liquid. Due to low velocities and the small scale, the freezing process is considered laminar. A transient approach is considered here and small time steps are required to reach a converged solution. A second order advection scheme is used, by setting the specified blend factor to 1.0. However, due to boundedness problems a high resolution scheme is used for calculations of mass fraction and energy [33]. When using this scheme, the blend factor will vary throughout the domain. This factor will be close to 1.0 (second order) in regions with low variable gradients and close to 0.0 (first order) in regions where the gradients change rapidly to prevent overshoots and undershoots and maintain robustness.

Numerical accuracy
Discretisation, iterative and modelling errors are considered in this paper. To investigate the discretisation errors a grid independence study is performed based on four subsequent grid sizes. The chosen key variable is t f . Due to the similar geometries, and since at this point only the numerical accuracy of the model is investigated, it is considered sufficient to study any of the three θ listed in Table 1 and θ = 90 • is selected for this study. From Table 2 it can be seen that the difference in t f between the grids is small, therefore the 5163 nodes grid is used for further studies. This is based on a balance between computer power and numerical accuracy around the freezing front. A time step analysis is also performed, revealing only a difference of 0.3% if a time step of 1e-05 s is used instead of 1e-04 s, suggesting that the larger time step can be used. The residual target is set to RMS = 1e-04, since the use of a smaller target (RMS = 1e-05) only gives a 0.3% difference in t f . The iterative errors are therefore investigated and not considered an issue here.

Results and discussion
The numerical results are first presented for the different angles according to Table 1, then they are compared with experimental results found in the literature and finally the impact of internal flow is investigated.

Influence of geometry and plate temperature on freezing time
To enable a later comparison to experiments in Jin et al. [4], simulations are performed for the droplets in Table 1 with V d = 9.32 μL that are placed on a surface having three different temperatures. For θ = 90 • and T p = −8.2 • C the simulations yield that the volume fraction ice, V ice in the droplet as a function of time is practically independent on the internal transport model applied, i.e. with or without internal natural convection, see Fig. 2.
Here, the experimental results from Jin et al. [4] is also included in the figure. The solid (green) line is simulations excluding internal natural convection, the dashed (blue) line is simulations including internal natural convection and the solid (black) line with dots is experimental data from Jin et al. [4]. Note that the dashed (blue) line is partly covered by the solid (green) line. The volume expansion of the droplet is included in the experimental data and therefore the volume fraction of ice will be larger than in the simulations. The interesting part, and also where the focus should be at this point, is that the behaviour of the ice formation is the same for both the experiments and the simulations. The ice formation is faster in the beginning of the freezing process due to the fact that the surface cooling the water is larger, but as more ice is formed a smaller area is cooling the ice, resulting in a slower ice formation. This can also be explained by the Stefan problem  [4]. Note that the dashed (blue) line is partly covered by the solid (green) line which suggests that the release of latent heat at the freezing front, the heat conduction through the ice and the small temperature gradients leads to smaller growth velocities (less ice is formed) closer to the end of the process. When, in Section 4.3, scrutinising the actual flow and temperature of the water in the droplet, distinct differences will be disclosed between the two cases of modelling but this is not reflected in the motion of the solid-liquid interface within the droplet. Additional simulations of the droplets presented in Table 1 are compared yielding that t f show a clear dependence on both θ and T p , see Table 3. An increase in t f of around 100% between minimum and maximal T p is observed independent of θ, indicating that there is no apparent interaction between T p and θ in the simulations. The influence of θ and A on t f is further investigated in Table 4, where it can be seen that changes in t f can be directly related to changes in A. There is however a decreased correspondence between the change in t f and A when θ = 77 and 90 • are compared. A possible reason for this is the increase of a with θ, i.e. if a is increased there is a larger distance for the heat to travel from top to bottom. Also, an increase in a corresponds to a smaller A. The results thus indicate that an increase in t f is expected when a is increased (and θ is increased). Note that the volume of the droplet is constant.

Model validation-comparison with experiments
The numerical simulations including internal natural convection (case 2) show that when the solid-liquid interface moves up in the droplet it has its highest position at the edge of the droplet and its lowest position in the middle of it, see Fig. 3 where water is indicated with blue colour (grey in non colour print) and ice with white. Here, water is assumed when x l ≥ 0.5. This is in qualitative agreement with the shape of the solid-liquid interfaces obtained in experiments in Jin et al. [4], see the drawn (white) lines in Fig. 3 being estimates from images in Jin et al. [4]. In Fig. 4 the height of the freezing front from both experiments and simulations can be seen. The solid (green) line represents a rough estimate of the ice water interface as derived from experimental images presented in Jin et al. [4] (taken in the middle of the droplet) and the dashed (blue) line is derived from simulations for θ = 90 • and T p = −8.2 • C when x l = 0.5 (taken furthest to the left in the droplet). From Figs. 3 and 4 it can be seen that the largest differences in the height of the ice front appear closer to the end of the freezing process. This is due to the volume expansion that accelerates rapidly in the beginning up to about half way of the freezing process, resulting in larger differences in the solid-liquid interfaces closer to the end. However, the similarity in shape and ratio of water and ice in the simulations and the experiments indicates that the droplets with fixed θ and shape in the model behave similar to real droplets.
When doing a quantitative comparison between simulations of droplets with θ = 77, 84 and 90 • , and experiments in Jin et al. [4] it is found that the difference in freezing time t f vary with temperature according to Table 3. The differences are both positive and negative and note that comparisons are done for all angles since it is difficult to determine θ in the images in Jin et al. [4]. Results for a droplet with θ = 77 • , i.e. corresponding to the initial values of the experiments presented by Jin et al. [4] and thus the most realistic angle, show a difference of less than 7% for T p = −11.2 • C, around 13% for T p = −8.2 • C and a maximum difference of around 25% for T p = −5.3 • C, see Table 3. The smallest difference is thus for the lowest temperature. A possible explanation is that this case is less influenced by the surrounding conditions because of the shorter freezing time. This explanation is supported by experimental results in Jin et al. [4], where the smallest dif-  Water is assumed when x l ≥ 0.5. The blue colour (grey in a non colour print) is water and white is ice as derived in the simulations. The drawn (white) lines are rough estimates of the ice-water interface as derived from experimental images presented in Jin et al. [4] ference between natural and forced convection is observed at the lowest temperature investigated T p = −11.2 • C. The results also indicate that conduction through the contact surface between substrate and droplet is dominating the heat transfer at the lower temperature, when compared to for To further investigate the influence of external heat transfer, a simulation when the adiabatic boundary condition at the droplet water-air surface is replaced with a conduction condition for T a = 21 • C. This is performed for the case of T p = −5.3 • C. The method for finding the heat transfer coefficient and how it was implemented in the model is outlined in Appendix A. Results from simulations show an increase in t f of less than 5% when applying heat transfer over the water-air surface, 22 s compared to 21 s. Although the increase in surface area (around 20%) is not accounted for in the simulation, the results point towards a minor effect of heat transfer from conduction at the droplet surface. Other surrounding conditions such as free and forced convection might however effect the results, but is not investigated here. The influence of conductivity in the substrate surface is furthermore considered negligible due to the high conductivity of red copper. All in all the agreement between simulation and experiments is good and it is therefore of interest to further investigate the influence of internal transport on the freezing droplet.

Internal flow investigation
Next, a comparison is made between simulations excluding internal natural convection (case 1) and including internal natural convection (case 2) with respect to average velocity, u and average temperature,T f . For all θ at all T p the difference inū in the water is the same during most part of the freezing process, except for a certain period of time when simulations with case 2 results in a much higherū. This is exemplified for a droplet with θ = 90 • cooling at T p = −8.2 • C in Fig. 5. The maximumū, which is about 100 times larger in case 2 as compared to case 1, occurs at t = 3 s for θ = 77 and 84 • and at t = 3.5 s for θ = 90 • , for all T p . The reason to the increased velocity can be traced to the shift in density gradient of water at T w = 4 • C, see Fig. 6. Hence, for case 2, gravity has a large effect on the internal flow. This conclusion is in accordance with the work by Kavanami et al. [12], and this behaviour has also been shown experimentally by Jin et al.,Enriquez et al. and Snoeijer and Brunet [4,17,21].
For case 2, the mixing is larger due to higher velocities, which results in a more uneven flow. This is also reflected in Fig. 7 where T w is very similar for the two cases, but in the temperature range T w = 2 − 4 • C the temperature is more unsteady for case 2 due the uneven flow mentioned. To further illustrate this, in Fig. 8 the velocity contours and direction of flow for case 1 and case 2 respectively are shown (for θ = 90 • , cooling at T p = −8.2 • C at t = 3.5 s). In case 1, the flow is only driven by the volume change that occurs as water turns to ice due to the density difference between the two materials. This transformation occur at times close to t = 0 near the bottom of the droplet   and later on close to the freezing front. This give rise to a pressure gradient in the direction from the ice to the top of the droplet resulting in a flow of water in this direction, as seen in the top of Fig. 8. In case 2, the flow is moving in a circular motion due to the large ΔT w inside the droplet. Since warmer water tends to flow to areas where there is colder water, in this case from the top of the droplet down to the freezing front, this give rise to the flow seen in the bottom of Fig. 8. Note that this circular flow pattern is only visible for T w = 2 − 4 • C, during other times the flow resembles more the flow in case 1.

Conclusions
The impact of internal flow and a fixed boundary on a water droplet freezing on a cold surface has been investigated numerically in this work. Droplets with three different contact angles, θ = 77, 84 and 90 • , but of equal volume, have been chosen based on experimental data. Two cases has been studied, one where internal natural convection has been included in the model and one where it has not. The results show that the droplets behave similar to droplets in experiments, and the freezing time, t f is predicted well for colder substrate temperature, T p (here, T p = −11.2 and −8.2 • C), but for higher substrate temperatures (here, T p = −5.3 • C) the disagreement is larger. This may be explained by influence from surrounding conditions in the experiments. By using these simplified models the complexity reduces if more realistic droplets want to be studied, like a droplet on a rotor blade exposed to external winds. The disadvantage is the lack of detailed information about the volume increase i.e. tracking and locating the free surface during freezing. A clear dependence of T p and θ on t f is furthermore observed in the numerical result. For the studied conditions, simulations indicate that the most important geometrical parameter to account for is the contact area between droplet and substrate. When investigating the internal flow it can be concluded that for all droplets the effects of gravity plays an important role for temperatures close to T w = 4 • C where the average velocity in the water is as highest. The results suggest that the gravity effects do not have to be accounted for with respect to t f , but for the internal flow, internal natural convection should be included. Further studies of the internal flow mechanisms including for example Marangoni convection are therefore recommended as well as more detailed experimental work.

Compliance with Ethical Standards
Conflict of interest On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.