As briefly mentioned in Section 4-1, it is speculated that the droplets of pure water observed in the present study were cooled down to ~234–238 K or less before freezing. In order to clarify the relationship between temperature and freezing time, a cooling curve (a change in temperature) of a droplet of pure water was numerically simulated.
A cooling curve of a droplet of pure water was simulated as follows. First, I assumed that the surface of the droplet was rapidly cooled by collision-free evaporation.
This evaporation process should occur when a mean free path of a water molecule is longer than the radius of the evaporating droplet [14]. The radius of the droplet employed in this experiment was 25–35 µm, which corresponds to a mean free path of a water molecule in the saturated water vapor at 265–270 K. The temperature of 265–270 K was lower than the nozzle temperature of 277–278 K, i.e., the temperature before droplet evaporation started. However, the temperature of the droplet was expected to become lower during the first several hundreds of microseconds, because the evaporative cooling rate of micrometer-sized water droplet was known to be faster than 104 K/s [14,15]. Therefore, I adopted collision-free evaporation without diffusion limited process. Second, I took into account the temperature gradient inside the droplet because evaporative cooling at the surface of the droplet was considered to be faster than the thermal conduction inside the droplet. To deal with the effect of the thermal conduction, the droplet was virtually divided into 100 shells with a uniform thickness, as shown in Figure 4-2. The inner-most and surface shells of the droplet were defined as the 1st and 100th shells, respectively.
Under these assumptions, the rate of the temperature change of the surface due to evaporative cooling is expressed by the following differential equation:
dTsurf
dt =(dm dt⁄ ) ∆Hvap(Tsurf)+&dQ99⁄ 'dt
Cp(Tsurf)Vsurf ρ(Tsurf) (4.1) with dm dt⁄ =−AsurfγP(Tsurf)%M⁄(2πRTsurf), where Tsurf is the surface temperature of the droplet, t the time elapsed after generation, m the mass of the droplet, Asurf the surface area of the droplet, g the evaporation coefficient, P(T) the vapor pressure, M the molecular weight, R the gas constant, ΔHvap(T) the enthalpy of vaporization, Cp(T) the heat capacity at constant pressure, Vsurf the volume of the surface shell of the droplet, r(T) the density of a droplet, and dQ99/dt the heat flow from the 99th shell to the surface shell. The change in mass of the droplet dm/dt was derived from the Knudsen equation, which is commonly applied to the collision-free evaporation without diffusion limited processes.
The evaporation coefficient g, which was used to correct for the difference between the experimental results and theoretical ones, was assumed to be 1 for water by referring to a recent study that measured the temperature of micrometer-sized droplets in a vacuum [15].
The heat flow dQ /dt was ignored at the initial state, i.e., t = 0 s, but added after the Figure 4-2 The droplet was virtually divided into 100 shells with a uniform thickness. All the shells of the droplet were numbered in the order of the shells from inside to outside.
temperature of the surface became lower than that of inside. It was calculated by using Fick’s first law of the diffusion equation:
dQn
dt =−Anκ(Tn)Tn+1 −Tn
r⁄100 (4.2)
where Qn is the heat flow from the n-th shell to the (n+1)-th shell, An the surface area of the n-th shell, k(T) the thermal conductivity, Tn the temperature of the n-th shell, and r the radius of the droplet. The heat flows other than dQ99/dt are used to calculate changes in temperature of the shells inside the droplet by using the following equation:
dTn+1
dt = dQn+1−dQn
Cp(Tn)Vnρ(Tn) (4.3)
where Vn is the volume of the n-th shell. I integrated a set of these differential equation with a 10-ns time step. As well as temperature, the radius r of the droplet was calculated from the mass m of the droplet at each time step by using the following equation:
m =6Vnρ(Tn)
99
n = 0
= 64πr3 3 78n+1
1009
3
− : n
100;3<ρ(Tn)
99
n = 0
. (4.4)
In the simulation, P,ΔHvap, Cp, r, and k were treated as functions of temperature. P and ΔHvap were calculated by using the equations proposed by Murphy and Koop [23], which were derived from the exponential fit to the data of Cp measured by Archer and Carter [24]. The exponential fit was not shown in their report. Thus, I used the following equation: Cp(T) / Jmol−1K−1 = 75.61 + 21.6 exp[− (T / K − 236) / 8.6]. r was calculated by using the sixth-order polynomial reported by Hare and Sorensen [25]. The thermal conductivity k was calculated by using the equation proposed by the International Association for the Properties of Water and Steam (IAPWS) [26].
Figure 4-3 shows temperatures of the surface and 10 selected shells out of 100 calculated for a droplet with an initial diameter and temperature of 49.2 µm and 277 K,
respectively. Even at 7 ms, the temperature difference between the inner-most and surface shells was about 2 K. Therefore, I defined the mass-averaged temperature, Tave
= Sn=099 [Tn Vnr(Tn)]/m, as a representative temperature of the droplet. An error in Tave
arising from the 2-µm resolution of size measurement was estimated to be ±0.4 K. On the other hand, an error caused by the initial temperature of the droplet was much smaller;
when I adjust the initial temperature of the droplet higher or lower by 1 K, Tave changed by less than ±0.1 K. Therefore, I ignored the uncertainty in the initial temperature. I further estimated errors arising from four thermodynamic parameters: ΔHvap, Cp, r, and k. I first recalculated Tave of the 49-µm droplet for the following cases: values of four parameters were fixed at T = 273.15 K and compared the result of numerical simulation with the result shown in Figure 4-3. Tave at t = 7 ms changed by +0.3 K when ΔHvap was fixed to 45051 Jmol−1; changed by +1.5 K when Cp was fixed to 75.90 Jmol−1K−1; changed by +0.1 K when r was fixed to 0.9999 kgm−3; and changed by −0.3 K when k was fixed to 0.5556 Wm−1K−1. The variance of Tave calculated for fixed ΔHvap, r, and k was smaller than the error in Tave arising from the uncertainty in droplet’s size. Since the variance of Tave should become much smaller if these three parameters were treated as functions of temperature, I concluded that the errors in Tave arising from ΔHvap, r, and k were negligible. Therefore, only the accuracy of Cp was taken into account to estimate the error in Tave. The data of Cp are reported not only by Archer and Carter but also by Angell et al. [27] and Tombari et al. [28]. All three data showed an exponential increase as the water temperature decreases, but the increases in the data of Angell et al. and Tombari et al. were slightly faster than that observed in the data of Archer and Carter.
Tave calculated with the data of Angell et al. and Tombari et al. was higher by +0.4 K and +0.9 K, respectively, than Tave calculated with the data of Archer and Carter. I concluded
that the error in Tave was mainly caused by the limited resolution of size measurement and by the accuracy of Cp. Tave in the present numerical simulation might be higher by +1.3 K or lower by −0.4 K.