Modelling investigation of water droplet evaporative freezing in artificial snowmaking

In this paper, a modelling investigation of water droplet evaporative freezing was conducted in order to better understand the snowmaking process and hence to optimise the design of the artificial snowmaking device. To this end, mass and heat transfer theoretical models of a single water droplet cooling in an air space were established and implemented in a numerical model developed using the software COMSOL Multiphysics. The effects of the air temperature, relative humidity and velocity and the water droplet initial diameter and temperature on this process were identified and analysed, and their appropriate ranges for the snowmaking were determined.


Introduction
Winter snowfall progressively heralds all manner of fun and games, such as snowball fights, skiing, tobogganing and snowboarding. The main challenge however is that in some countries the snowfall is not enough and significantly changes with nature. It is therefore judicious to use the artificial snowmaking. There are three common modes of artificial snowmaking, which are the snowmaking by compressed air, high pressure water and ice-cutting. For the snowmaking by compressed air, the water is mixed with high pressure air and blown using a nozzle as a form of tiny droplets, which are then frozen in cold air to form snowflakes. For the snowmaking by high pressure water, the water is sprayed as tiny droplets thanks to an atomising nozzle in a cold air, forming thus snowflakes. For the snowmaking by ice-cutting, ice blocks are scraped with blades to form extremely thin ice sheets as snowflakes.
In this context, many researchers carried out extensive works on snowflakes in three research directions. The first direction is to study the impact of the snow distribution on the ecological environment and the natural disasters [1][2][3][4]. The second direction is to study the nature of snow, mainly its physical properties such as density, viscosity, porosity, structure [5,6], as well as its optical properties such as reflectivity [7][8][9][10]. The third direction is to study the growth characteristics of ice crystals, including their nucleation and polymerisation processes, under different operating conditions [11][12][13][14]. In 1930s, Nakaya [15] conducted an in-depth study on the growth morphology of ice crystals for the first time. Carignano et al. [16] investigated the growth of ice crystals on the prism and base * e-mail: georges.elachkar@tjcu.edu.cn surfaces by molecular dynamics simulation. The kinetics and time evolution of pure water and brine freezing were studied. The experiments showed that the total freezing time of the brine is longer than that of the pure water. The authors attributed this phenomenon to the influence of the ion clusters on the cluster formation ability. Nelson [17] provided a model that was improved by estimating the diffusion field adjacent to the growing ice crystal and using various ice crystal growth data. He experimentally obtained a large amount of data to establish an ice crystal growth prediction model to analyse the ice below -22 • C. Evans [18] studied the growth and fragmentation of ice crystals under electric field, especially the growth rate of ice crystals in the direction of electric field under the supersaturation temperature. He found that, when a live wire was placed in a cold cloud, it would produce thousands of ice crystals. These crystals are formed by the growth of ice crystals on the wires, which are then broken by the stress of the strong electric field of the wires.
There are many forms of ice crystal, but there is no clear model for their growth process due to its great complexity. The growth state of ice crystal is greatly affected by the environment and slight changes in environmental conditions result in different ice crystal morphologies. This sensitivity is one of the signs of ice crystal growth. Even at a qualitative level, it is difficult to accurately explain the growth mechanism of ice crystal under different conditions.
Despite the numerous aforementioned studies, less attention was paid to the evaporative freezing process of water droplets from the design point of view of the artificial snowmaking device. Therefore, a modelling investigation of this process was carried out and presented in this work.
The effects of the air temperature, relative humidity and velocity and the water droplet initial diameter and temperature on this process were identified in order to get their appropriate ranges for the snowmaking, which can provide some elements for the design of the artificial snowmaking device.

Modelling
The evaporative freezing process of water droplets is essentially the cooling process of these droplets from their initial temperature until their transformation into snowflakes (0 • C) under the action of the surrounding air. In order to characterise this process, mass and heat transfer theoretical models of a single water droplet freezing in an air space were established and implemented in a numerical model developed using the software COMSOL Multiphysics.

Mass transfer
According to Fick first law, the water at the droplet interface evaporates due to the difference of concentration between the droplet and the air. The droplet evaporation rate (ṁ d ), equal to the droplet mass (m d ) loss over time (t), is expressed as follows: Where A i is the droplet interface area, k w is the water mass transfer coefficient, M w is the water molar mass, C w,sv and C w,∞ are respectively the water saturated vapour concentration and vapour concentration in air, p w,s and p w,∞ are respectively the water saturation pressure and vapour pressure in air, T d and T a are respectively the droplet and air temperatures and R=8.3145 J.mol −1 .K −1 is the universal gas constant. The p w,s can be obtained by the following empirical formula [19]: and c 6 =6.5459673 are polynomial coefficients.
The k w can be obtained by the droplet Sherwood number (S h d ) formula as follows: Where D d is the droplet diameter, γ w is the diffusion coefficient of water vapour in air, Re d and S c d are respectively the droplet Reynolds and Schmidt numbers, ρ w and µ w are respectively the water density and dynamic viscosity and u d is the droplet velocity. The evaporation rate of the droplet can also be expressed as a function of its diameter as follows: Where V d is the droplet volume. According to Eqs. 1&8, the variation over time of the droplet diameter due to the evaporation can be written as follows:

Heat transfer
By applying the energy balance on the droplet during its movement in the air, the energy equation is expressed as follows: The sensible heat Q sens resulting from the droplet temperature variation is defined as follows: The latent heat Q lat resulting from the droplet evaporation is defined as follows: Where v,w is the water latent heat of vaporisation. The convective heat transfer Q conv at the droplet interface is defined as follows: Where h is the convective heat transfer coefficient, Nu d is the droplet Nusselt number and Pr w , c p,w and λ w are respectively the water Prandtl number, specific heat and thermal conductivity. The radiative heat transfer Q rad between the droplet and its surrounding air is defined as follows: Where w is the water emissivity and σ is the constant of The Q rad is neglected because the temperature difference between the droplet and the air is relatively small (4−8 • C). Using the Eqs. 1,5,11−14, the energy equation (Eq. 10) becomes: The integration of Eq. 9 gives: By substituting Eq. 19 into Eq. 18, the relationship between the droplet temperature and other parameters can be obtained as follows:

Numerical model
In order to develop the numerical model, the laminar flow, transport of diluted species and heat transfer modules were selected according to the physical phenomena involved in the droplet evaporative freezing process.

Physical model
The considered physical model consists of a water droplet placed in the center of an air space, as depicted in Fig. 1.

Mathematical model
The used governing equations are written as follows: The inlet boundary conditions are as follows: The outlet boundary condition is as follows: −n.q = 0 (28)

Assumptions
In order to simplify the model, the following assumptions were considered: • The air space is infinite and the effect of the droplet evaporation on the air concentration is neglected; • The droplet is perfectly spherical with constant thermophysical properties; • The radiative heat transfer between the droplet and the air is neglected.

Numerical simulations
Different numerical simulations were carried out by varying the values of the air temperature (-4−4 • C), relative humidity (20−100 %) and speed (2−10 m.s −1 ), and the droplet initial diameter (100−500 µm) and initial temperature (4−12 • C). The data processing and analysis are presented in the next section. According to Eq. 20, the mass and heat transfers at the droplet interface depend on many parameters, such as the air temperature, relative humidity and velocity, the droplet initial diameter and temperature, the snowmaking duration, etc. If some of these parameters are fixed as contraint conditions for the snowmaking device, the other parameters can thus be calculated. For instance, if the snowmaking duration is 2 s, the droplet initial temperature is 4 • C and the air speed is 2 m.s −1 , the evolutions of the water droplet critical diameter, below which the droplet temperature can reach 0 • C, as a function of the air temperature and relative humidity are shown in Figs. 2−4.

Water droplet critical diameter (µm)
Air temperature (°C) Water droplet critical diameter (µm) Air relative humidity (%) The water droplet critical diameter decreases with the increase of the air temperature regardless of the air relative humidity, and the higher the air relative humidity is the more obvious the air temperature influence is (Fig. 3).
In addition, the water droplet critical diameter decreases with the increase of the air relative humidity regardless of the air temperature, and the higher the air temperature is the more obvious the air relative humidity influence is (Fig. 4). When the air temperature or relative humidity is high, the water droplet critical diameter is relatively small, which requires a higher atomisation ability for the snowmaking device. Therefore, the water spraying parameters should be automatically adjusted according to the air temperature and relative humidity to generate droplets with appropriate diameters. It should be noted that, owing to the evaporation effect, the water droplet temperature can reach 0 • C even if the air temperature is higher than 0 • C, but strictly for a low range of air relative humidity and water droplet initial diameter. Moreover, when the air temperature is below 0 • C, the water droplet initial diameter also needs to be carefully chosen in order to reduce the snowmaking duration. For instance, for a snowmaking duration of 2 s with air temperature of -4 • C and relative humidity of 40 %, the water droplet critical diameter is 760 µm.

Droplet freezing time
The droplet freezing time is defined as the time required for cooling the droplet from its initial temperature to its freezing temperature (0 • C).    5 shows the evolutions of the water droplet temperature as a function of time for droplet initial temperature of 4 • C and diameter of 100 µm and air speed of 2 m.s −1 , relative humidity of 20 % and different temperatures. The freezing time decreases with the air temperature decrease due to the increase of the temperature difference between the droplet and the air, leading to an increase of the convective heat transfer. However, as the freezing time decreases, the droplet evaporated mass decreases, leading to the droplet cooling slowdown, which explains the nonlinear variation of the freezing time as a function of the air temperature. If the air relative humidity is increased, the droplet evaporation rate decreases regardless of the air temperature, leading to an increase of the freezing time.   6 shows the evolutions of the water droplet temperature as a function of time for droplet initial temperature of 4 • C and diameter of 100 µm and air temperature of -2 • C, relative humidity of 20 % and different speeds. The freezing time decreases with the increase of the air speed due to the increase of the heat transfer coefficient at the droplet interface, leading to an enhancement of the convective heat transfer. When the relative humidity is equal to 20 %, 40 %, 60 %, 80 % and 100 %, the difference of the freezing times corresponding to the minimum and maximum air speeds (i.e. 2 and 10 m.s −1 ) is 0.007 s, 008 s, 0.01 s, 0.012 and 0.016 s, respectively. Hence, the convective heat transfer has a relatively weak influence on the evaporative freezing process of water droplet. Therefore, it is not necessary to apply an excessive air speed in the artificial snowmaking.    7 shows the evolutions of the water droplet temperature as a function of time for air speed of 2 m.s −1 , temperature of -2 • C and relative humidity of 20 % and droplet initial temperature of 4 • C and different initial diameters. The freezing time increases with the increase of the droplet initial diameter. Indeed, when the droplet initial diameter increases, the droplet interfacial area increases, leading to an increase of the convective heat transfer and hence to a decrease of the freezing time. However, as the droplet mass increases, a greater heat quantity should be evacuated and hence a longer freezing time is required. The competition between these two factors results in an increase of the freezing time. Water droplet temperature (°C) Time (s) Figure 8. Evolutions of the water droplet temperature as a function of time for air speed of 2 m.s −1 , temperature of -2 • C and relative humidity of 20 % and droplet initial diameter of 100 µm and different initial temperatures. Fig. 8 shows the evolutions of the water droplet temperature as a function of time for air speed of 2 m.s −1 , temperature of -2 • C and relative humidity of 20 % and droplet initial diameter of 100 µm and different initial temperatures. The freezing time increases with the droplet initial temperature due to the increase of the evacuated heat quantity. However, as the convective heat transfer is greater for higher droplet initial temperature, the increasing rate of the freezing time decreases.
From Figs. 5−8, one can easy notice that the water droplet temperature variations can be roughly divided into three stages: they firstly decrease steadily, then decrease rapidly and finally decrease linearly. In the first stage, the initial temperature of the droplet progressively decreases from the interface towards the centre. Therefore, the whole droplet temperature variation is delayed, resulting in a low cooling rate in this stage. In the second stage, a greater cooling rate is involved. However, this cooling rate decreases and gradually becomes stable in the last stage. The reason is that the saturated vapour pressure at the droplet interface decreases due to the lower droplet temperature, leading to a reduction of the heat absorbed by the evaporation. In addition, the decrease of the convective heat transfer also reduces the droplet cooling rate.