# Parameter Sensitivity of the Microdroplet Vacuum Freezing Process

Authors: Zhijun Zhang; Yuekai Zhang; Lili Zhao; Wenhui Zhang; Shuangshuang Zhao
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The vacuum freezing process of microdroplets (<100μm in diameter) is studied by dynamic mesh method. The mass transfer coefficient was studied using the results of related papers that considered droplet diameters exceeding 1 mm. The diameter, initial temperature, and vacuum chamber pressure effects are also discussed. To estimate parameter sensitivity, the effects of material density, specific heat, and thermal conductivity in 20% scope, as well as latent evaporation/sublimation in 5%, were simulated. The results show that the mass transfer coefficient K is essentially different between microdroplets (<100 μm) and macrodroplet (>1 mm). Pressure and droplet diameter have an effect on cooling and freezing stages, but initial temperature only affects the cooling stage. The thermal conductivity coefficient k
l affected the cooling stage, whereas k
i affected the freezing stage. Heat capacity C
l affected the cooling stage, but C
i has virtually no effect on all stages. The actual latent heat of freezing Δ
H was also affected. Higher density corresponds to lower cooling rate in the cooling stage.


---

## Body

## 1. Introduction

The use of evaporative supercooling method for ice production presents an immensely popular topic for research [1, 2]. Droplet vacuum freezing involves a rapid decrease in temperature in which the phase changes from liquid to solid, from liquid to vapor, and from solid to vapor [3]. The basic principle is as follows. Water evaporates under vacuum conditions. If the vapor pressure of the atmosphere is below 611 Pa (vapor pressure at 0°C), then the water continues to evaporate. Heat must be provided during evaporation to decrease droplet temperature. Thus, the temperature reaches below 0°C and then ice is produced.The heat and mass transfer of droplet vacuum freezing are complicated processes and have been studied by numerical and experiment methods. Kim et al. [4, 5] conducted theoretical and experimental studies on the use of water spray evaporation method to produce ice particles. In their study, the conditions for the formation of ice particles were theoretically investigated using a diffusion-controlled evaporation model. The experiment included pure water and 7% ethylene glycol. Li et al. [6] experimentally studied the cooling/freezing phenomena of a water droplet due to evaporation in an evacuated chamber, as well as the heat transfer dominating the evaporation/freezing phenomena, to estimate the pressure in the evaporator. Asaoka et al. [7, 8] also proposed an alternative vacuum freezing method to produce ice slurry by evaporating ethanol solution instead of water. This system improves ice forming stability but requires a large amount of energy for vacuum production. These models and simulations have been used as analytical model, and finite element method has been used to examine the process [9]. This method is based on the diffusion model through a porous medium, and radius reduction attributed to water evaporation is disregarded. We have studied the droplet vacuum freezing process using the dynamic mesh method [10], in which the droplet diameter was reduced during the freezing process. The diameter of droplet was set as 7.5, 10.5, and 12.5 mm for comparison with experiment results [9]. To continuously produce ice from water droplets, the diameter of the droplet should be less than 100 μm because of the height limitation of the spray vacuum chamber [3].In the present study, a microdroplet (<100μm diameter, compared with the conventional >1 mm) vacuum freezing process with dynamic mesh method is studied. The mass transfer coefficient was studied along with the results of related studies [3, 10]. The diameter, initial temperature, and vacuum chamber pressure effect were also studied. To estimate parameter sensitivity, the effect of material density, specific heat, thermal conductivity, and latent heat of evaporation was simulated in 20% or 5% scope.


## 2. Material and Method

### 2.1. Material

The microdroplet is spherically shaped and its diameter is reduced during vacuum freezing process because of surface evaporation. To simplify computations, the two-dimensional axis symmetry model (Figure1) was adopted. The solid and liquid properties of the microdroplet are provided in Table 1. To estimate sensitivity, parameters including ρ, C
l, C
i, k
i, and k
l were changed by ±20%. By contrast, the parameters λ
i and λ
l were just changed by ±5%. These adjustments are adopted because λ
i must be larger than λ
l and these two parameters are highly sensitive. These reasons will be discussed later.Table 1
Parameters used in the simulation process.
Parameter
Initial value
+20% or 5% value
−20% or 5% value
Unit
ρ
998
1197.6
789.4
Kg m−3
C
l
3420
4104
2736
J kg−1 K−1
C
i
1720
2064
1376
J kg−1 K−1
k
l
0.6
0.72
0.48
W m−1 K−1
k
i
2
2.4
1.6
W m−1 K−1
λ
l
2.5× 106
2.625× 106
2.375× 106
J kg−1
λ
i
2.8× 106
2.94× 106
2.66× 106
J kg−1Figure 1
2D axis symmetry model of the microdroplet.


### 2.2. Heat Transfer Model

By considering a hypothesis of local thermal equilibrium, energy conservation is reduced to a unique equation:(1)
ρ
C
∂
T
∂
t
-
∇
k
·
∇
T
=
0
.The boundary conditions for1 on the axis symmetric surface are as follows:
(2)
k
∇
T
=
0
.The boundary conditions for1 on the outer surface are defined as
(3)
k
∂
T
∂
n
=
-
λ
·
I
˙
.The initial condition for1 is
(4)
T
-
i
=
T
0
.


### 2.3. Mass Transfer Model

Mass transfer occurs through a phase change of water from liquid or solid to gas on the droplet surface. However, mass transfer is not considered inside the droplet. A general expression of nonequilibrium evaporation rate used for phase change modeling in porous medium is given by [11, 12]
(5)
I
˙
=
K
ρ
v
,
eq
-
ρ
v
.This equation could also be expressed as [13]
(6)
I
˙
=
K
m
v
P
sat
⁡
-
P
v
R
T
.In vacuum cooling, the equation could be expressed as [14–18]
(7)
I
˙
=
4
φ
d
h
m
P
sat
⁡
-
P
.Here, phase change occurs based on the movement of the pressure difference between the saturation vapor pressure in the microdroplet surface and vacuum chamber pressure [11], which is also considered the mass transfer boundary condition:
(8)
I
˙
=
K
P
sat
⁡
-
P
v
.The mass transfer coefficient,K, is a rate constant parameter with the dimension of reciprocal time in which phase change occurs. A large value of K signifies that phase change occurs within a short amount of time. The K value is conventionally obtained from experiments. In the present paper, we decided to use experimental data from [3]. Detailed information will be provided later.Water partial pressure,P
v, is assumed to be equal to the total vacuum chamber pressure. The saturation pressure, P
sat
⁡, is related to the surface temperature of the microdroplet. This value is estimated using the equation when the surface temperature is higher than freezing temperature [9] but uses another equation for estimation when the surface temperature is lower than freezing temperature [9]. The details could be found in related papers [9, 19].The initial condition for5 is set using the initial diameter of the droplet:
(9)
d
=
d
0
.
Thus, as in the dynamic mesh method, the boundary conditions of 8 could be rewritten as
(10)
v
˙
d
=
K
P
sat
⁡
-
P
v
ρ
.


### 2.4. Resolution and Validation

COMSOL Metaphysics 3.5a was used to solve the set of equations. COMSOL is an advanced software used for modeling and simulating any physical process described by partial derivative equations. The set of equations introduced above was solved using the corresponding relative initial and boundary conditions. COMSOL offers three possibilities for writing the equations:(
1
) using a template (Fick’s law or Fourier’s law), (
2
) using the coefficient form (for mildly nonlinear problems), and (
3
) using the general form (for most nonlinear problems). Differential equations in the coefficient form were written using a nonsymmetrical pattern multifrontal method. We used a direct solver for sparse matrices (UMFPACK), which involves significantly more complicated algorithms than those solvers used for dense matrices. The need for efficient handling of the fill-in factors L and U is considered the method’s major complication.A two-dimensional (2D) axis symmetry grid was used to solve the equations using COMSOL Metaphysics 3.5a. The mesh consists of 1068 elements (2D), and time stepping is freely taken by the solver. A backward differentiation formula was used to solve time-dependent variables. Relative tolerance was set at 1 × 10−4, whereas absolute tolerance was set at 1 × 10−6. The simulations were performed using a Tongfang PC running on Windows 7, with an Intel Core 2 Duo processor clocked at 3.0 GHz processing speed and 4096 MB of RAM.Several grid sensitivity tests were conducted to determine the sufficiency of the mesh scheme and to ensure that the results are grid independent. Figures2 and 3 present the grid test results. The vacuum chamber pressure would be higher or lower than 600 Pa, depending on the different physical process. In the transition involved in Figure 2, liquid water is changed to vapor. In Figure 3, liquid water changes into solid ice, and then the conversion from ice to vapor is attributable to the physical property difference. The test results in Figures 2 and 3 show that all of the three elements yielded accurate results, although certain numerical results have step-changed. Simulation times of 0.5, 1, and 3 h presented 267, 1068, and 4602 elements in Figure 3. In the following simulation, 1068 elements were used.Figure 2
Grid sensitivity tests with 267, 1068, and 4272 elements at 600 Pa vacuum chamber pressure.Figure 3
Grid sensitivity tests with 267, 1068, and 4272 elements at 400 Pa vacuum chamber pressure.As presented above, the mass transfer coefficient or rate constantK should be obtained in the experiment. In the present simulation, the experimental data from [3] was used for decision-making. In Figure 4, different K values were used, but the results do not show the difference when K exceeds 0.1 because, at K
≥
0.1, mass and heat transfer is limited by inner microdroplet process. In addition, the obtained simulation results were used to approximate the experimental data. The calculation results in the reference are also shown in Figure 4. The bias between the experimental present simulation results and reference results [3] is attributed to the fact that the microdroplet is not spherical during vacuum freezing. The droplet surface is in the boiling state, and thus the increased surface area results in larger mass and heat transfer than in the simulation. Certain information on the inner controlled results could be obtained later.Figure 4
Simulation results with differentK values (10, 1, 0.1, 0.01, and 0.001 kg Pa−1 s−1) and the experiment and calculation results from [3].In our previous study [10], the rate constant K is 4 × 10−6, which is considerably lower than 0.1 in the present study. The former diameter of droplet is 7.5 mm or larger. However, the later diameter is less than 100 μm. The droplet dimension was reduced, thereby altering the mass and heat transfer process.


## 2.1. Material

The microdroplet is spherically shaped and its diameter is reduced during vacuum freezing process because of surface evaporation. To simplify computations, the two-dimensional axis symmetry model (Figure1) was adopted. The solid and liquid properties of the microdroplet are provided in Table 1. To estimate sensitivity, parameters including ρ, C
l, C
i, k
i, and k
l were changed by ±20%. By contrast, the parameters λ
i and λ
l were just changed by ±5%. These adjustments are adopted because λ
i must be larger than λ
l and these two parameters are highly sensitive. These reasons will be discussed later.Table 1
Parameters used in the simulation process.
Parameter
Initial value
+20% or 5% value
−20% or 5% value
Unit
ρ
998
1197.6
789.4
Kg m−3
C
l
3420
4104
2736
J kg−1 K−1
C
i
1720
2064
1376
J kg−1 K−1
k
l
0.6
0.72
0.48
W m−1 K−1
k
i
2
2.4
1.6
W m−1 K−1
λ
l
2.5× 106
2.625× 106
2.375× 106
J kg−1
λ
i
2.8× 106
2.94× 106
2.66× 106
J kg−1Figure 1
2D axis symmetry model of the microdroplet.


## 2.2. Heat Transfer Model

By considering a hypothesis of local thermal equilibrium, energy conservation is reduced to a unique equation:(1)
ρ
C
∂
T
∂
t
-
∇
k
·
∇
T
=
0
.The boundary conditions for1 on the axis symmetric surface are as follows:
(2)
k
∇
T
=
0
.The boundary conditions for1 on the outer surface are defined as
(3)
k
∂
T
∂
n
=
-
λ
·
I
˙
.The initial condition for1 is
(4)
T
-
i
=
T
0
.


## 2.3. Mass Transfer Model

Mass transfer occurs through a phase change of water from liquid or solid to gas on the droplet surface. However, mass transfer is not considered inside the droplet. A general expression of nonequilibrium evaporation rate used for phase change modeling in porous medium is given by [11, 12]
(5)
I
˙
=
K
ρ
v
,
eq
-
ρ
v
.This equation could also be expressed as [13]
(6)
I
˙
=
K
m
v
P
sat
⁡
-
P
v
R
T
.In vacuum cooling, the equation could be expressed as [14–18]
(7)
I
˙
=
4
φ
d
h
m
P
sat
⁡
-
P
.Here, phase change occurs based on the movement of the pressure difference between the saturation vapor pressure in the microdroplet surface and vacuum chamber pressure [11], which is also considered the mass transfer boundary condition:
(8)
I
˙
=
K
P
sat
⁡
-
P
v
.The mass transfer coefficient,K, is a rate constant parameter with the dimension of reciprocal time in which phase change occurs. A large value of K signifies that phase change occurs within a short amount of time. The K value is conventionally obtained from experiments. In the present paper, we decided to use experimental data from [3]. Detailed information will be provided later.Water partial pressure,P
v, is assumed to be equal to the total vacuum chamber pressure. The saturation pressure, P
sat
⁡, is related to the surface temperature of the microdroplet. This value is estimated using the equation when the surface temperature is higher than freezing temperature [9] but uses another equation for estimation when the surface temperature is lower than freezing temperature [9]. The details could be found in related papers [9, 19].The initial condition for5 is set using the initial diameter of the droplet:
(9)
d
=
d
0
.
Thus, as in the dynamic mesh method, the boundary conditions of 8 could be rewritten as
(10)
v
˙
d
=
K
P
sat
⁡
-
P
v
ρ
.


## 2.4. Resolution and Validation

COMSOL Metaphysics 3.5a was used to solve the set of equations. COMSOL is an advanced software used for modeling and simulating any physical process described by partial derivative equations. The set of equations introduced above was solved using the corresponding relative initial and boundary conditions. COMSOL offers three possibilities for writing the equations:(
1
) using a template (Fick’s law or Fourier’s law), (
2
) using the coefficient form (for mildly nonlinear problems), and (
3
) using the general form (for most nonlinear problems). Differential equations in the coefficient form were written using a nonsymmetrical pattern multifrontal method. We used a direct solver for sparse matrices (UMFPACK), which involves significantly more complicated algorithms than those solvers used for dense matrices. The need for efficient handling of the fill-in factors L and U is considered the method’s major complication.A two-dimensional (2D) axis symmetry grid was used to solve the equations using COMSOL Metaphysics 3.5a. The mesh consists of 1068 elements (2D), and time stepping is freely taken by the solver. A backward differentiation formula was used to solve time-dependent variables. Relative tolerance was set at 1 × 10−4, whereas absolute tolerance was set at 1 × 10−6. The simulations were performed using a Tongfang PC running on Windows 7, with an Intel Core 2 Duo processor clocked at 3.0 GHz processing speed and 4096 MB of RAM.Several grid sensitivity tests were conducted to determine the sufficiency of the mesh scheme and to ensure that the results are grid independent. Figures2 and 3 present the grid test results. The vacuum chamber pressure would be higher or lower than 600 Pa, depending on the different physical process. In the transition involved in Figure 2, liquid water is changed to vapor. In Figure 3, liquid water changes into solid ice, and then the conversion from ice to vapor is attributable to the physical property difference. The test results in Figures 2 and 3 show that all of the three elements yielded accurate results, although certain numerical results have step-changed. Simulation times of 0.5, 1, and 3 h presented 267, 1068, and 4602 elements in Figure 3. In the following simulation, 1068 elements were used.Figure 2
Grid sensitivity tests with 267, 1068, and 4272 elements at 600 Pa vacuum chamber pressure.Figure 3
Grid sensitivity tests with 267, 1068, and 4272 elements at 400 Pa vacuum chamber pressure.As presented above, the mass transfer coefficient or rate constantK should be obtained in the experiment. In the present simulation, the experimental data from [3] was used for decision-making. In Figure 4, different K values were used, but the results do not show the difference when K exceeds 0.1 because, at K
≥
0.1, mass and heat transfer is limited by inner microdroplet process. In addition, the obtained simulation results were used to approximate the experimental data. The calculation results in the reference are also shown in Figure 4. The bias between the experimental present simulation results and reference results [3] is attributed to the fact that the microdroplet is not spherical during vacuum freezing. The droplet surface is in the boiling state, and thus the increased surface area results in larger mass and heat transfer than in the simulation. Certain information on the inner controlled results could be obtained later.Figure 4
Simulation results with differentK values (10, 1, 0.1, 0.01, and 0.001 kg Pa−1 s−1) and the experiment and calculation results from [3].In our previous study [10], the rate constant K is 4 × 10−6, which is considerably lower than 0.1 in the present study. The former diameter of droplet is 7.5 mm or larger. However, the later diameter is less than 100 μm. The droplet dimension was reduced, thereby altering the mass and heat transfer process.


## 3. Results and Discussion

The process of microdroplet vacuum freezing includes three stages: cooling, freezing with phase change, and frozen stages. In the cooling stage, the temperature of the microdroplet is lowered to the phase change temperature from the initial temperature. In the freezing stage, the liquid water is changed to solid ice, and the latent heat of the phase change is released. The temperature was maintained at a stable value. When all liquid water is turned into solid ice, the temperature of the microdroplet is rapidly decreased to the saturation temperature corresponding to the vacuum chamber pressure. This state is the frozen stage.When the external conditions, such as pressure, initial temperature, and microdroplet diameter, or internal conditions, such as the property of material, were changed, the three stages of microdroplet vacuum freezing are also changed. The pressure is the only condition that could be controlled in the vacuum freezing process. The saturation temperature corresponds to the pressure, with lower pressure corresponding to lower temperature and vice versa. However, the pressure not only affects the end temperature but also affects the process. When the pressure is lowered, the freezing time is evidently shortened (Figure5). Except at 600 and 500 Pa, the three stages presented above were all altered. The phase change is a process. The change from liquid water to solid ice occurs at a range of temperature. Pressures at 600 and 500 Pa are located near the triple-phase temperature, and thus the process occurs differently for other pressures. In the present simulation, the temperature scope of the phase change is 2 K from triple-phase temperature point. Results show that the freezing time is evidently different when the pressure is lower than 500 Pa. The three stages are affected by the pressure in the process.Figure 5
Effect of vacuum chamber pressure (600, 500, 400, 300, 200, and 100 Pa) with initial temperature 286 K and microdroplet diameter of 50μm.However, the initial temperature effect is different from the pressure effect in Figure6. The initial temperature mainly affects the cooling stage, almost exerting no effect on the freezing and frozen stages. Based on the results, if pressure is not changed, then the end temperature of the microdroplet remains the same regardless of the initial temperature of the material or microdroplet.Figure 6
Effect of initial temperature (286, 306, 326, 346, and 366 K) at a vacuum chamber pressure 300 Pa and microdroplet initial diameter 50μm.Similar to pressure and initial temperature differences, the initial diameter of microdroplet has a considerable effect on the freezing process. The diameter of microdroplet ranges from 10μm to 100 μm and is a one-order dimension. However, the freezing time exceeds two orders in Figure 7. In reality, a lower initial diameter results in flash evaporation [5] and a decreasing diameter increases this probability.Figure 7
Effect of microdroplet diameter (10, 20, 30, 40, 50, 60, 70, 80, 90, and 100μm) at a vacuum chamber pressure of 300 Pa and microdroplet initial temperature of 286 K.Unlike rate constantK, the process is externally controlled. When K
>
0.1, the process is internally controlled. These changes indicate process complexity.Unlike external conditons, material property is related to the temperature. However, the accuracy of the property parameters is the basis for the validity of simulation and models. Thermal conductivity is related to temperature and varies with the process. The manner by which this parameter affects the process should be presented. Thermal conductivity is related to both temperature and physical state, varying in the liquid and the solid states. To estimate parameter sensitivity, the value of parameter is changed in the 20% scope, except for latent heat of phase change. The results show the evident effect of thermal conductivityk
l and k
i, which are different between Figures 8 and 9. Larger k
l results in faster cooling rate during the cooling stage. By contrast, the frozen stage is less affected. The value of k
i is less affected in the cooling stage and is evidently affected in the freezing stage. Larger k
i results in decreased freezing time.Figure 8
Effect of thermal conductivityk
l at an initial temperature of 286 K, vacuum chamber pressure 300 Pa, and microdroplet initial diameter of 50 μm.Figure 9
Effect of thermal conductivityk
i at initial temperature of 286 K, vacuum chamber pressure of 300 Pa, and microdroplet initial diameter of 50 μm.As seen in Figures10 and 11, the heat capacities C
l and C
i have different effects on the freezing process. C
i almost has no effect on the freezing process. By contrast, C
l evidently affects the first stage (or the cooling stage). Larger C
l results in a higher cooling rate. However, the end time is nearly similar.Figure 10
Effect of heat capacityC
l at initial temperature of 286 K, vacuum chamber pressure of 300 Pa, and microdroplet initial diameter of 50 μm.Figure 11
Effect of heat capacityC
i at an initial temperature of 286 K, vacuum chamber pressure of 300 Pa, and microdroplet initial diameter of 50 μm.As the main parameter, the latent heatsλ
l and λ
i present an evident effect on all processes, as seen in Figures 12 and 13. However, their effects are different on the freezing process. In fact, in the second stage, namely, the freezing stage, liquid water is changed into solid microdroplets. The released latent heat of freezing, Δ
H, is calculated as
(11)
Δ
H
=
λ
i
-
λ
l
.Figure 12
Effect of latent heatλ
l at an initial temperature of 286 K, vacuum chamber pressure of 300 Pa, and microdroplet initial diameter of 50 μm.Figure 13
Effect of latent heatλ
i at an initial temperature of 286 K, vacuum chamber pressure 300 Pa, and microdroplet initial diameter of 50 μm.By increasingλ
l, the Δ
H value decreases. Given that the latent heat of freezing Δ
H is decreased, cooling and freezing time decreased (Figure 12) and vice versa. The latent heat of sublimation λ
i also varies. Increasing λ
i results in the increase in latent heat of freezing Δ
H, and then the cooling rate in cooling stage is decreased. Thus, freezing time is also increased because of the increased Δ
H and vice versa.In the heat loss in3, given the increase in latent heat values λ
l and λ
i, the cooling rates in the cooling stage shown in Figures 12 and 13 are both increased. As discussed above, the rate constant K is fixed, and the process is internally controlled. Here, 3 is unaffected by the heat transfer process by the changed λ
l and λ
i, obtained as the highest value as in K
>
0.1. These heat values represent the latent heat of freezing Δ
H. The latent heat values λ
l and λ
i changed in the 5% scope because of higher sensitivity.Finally, density is discussed. In the phase change from liquid water to solid ice, density changes. As shown in Figure14, higher density is correlated with lower cooling rate in the cooling stage, and freezing time is consequently shortened. This trend is true vice versa.Figure 14
Effect of densityρ at an initial temperature of 286 K, vacuum chamber pressure of 300 Pa, and microdroplet initial diameter of 50 μm.


## 4. Conclusion

The mass transfer on the droplet surface varies between microdroplets and macrodroplets, as represented by the mass transfer coefficientK. The pressure of the vacuum chamber and microdroplet diameter affect the cooling and freezing stages. However, the initial temperature of the microdroplet only affects the cooling stage. The thermal conductivity coefficient k
l affected the cooling stage, whereas k
i affected the freezing stage. Heat capacity C
l affected the cooling stage, but C
i almost has no effect in all stages. Latent heat of λ
l and λ
i affected cooling and freezing stages, but their effects are different. Both latent heat values affected the latent heat of freezing Δ
H. Finally, higher density results in decreased cooling rate during the cooling stage.



---

*Source: 370159-2015-08-11.xml*