## 3 Numerical Method
### Gas-Liquid Multiphase Model  
Given incompressible and laminar flow, the governing equations of continuity and momentum are as below:  
\[\
abla\times\overset{\rightarrow}{u}=0 \tag{10}\]  
\[\frac{\partial}{\partial t}[\overset{\rightarrow}{u}]+\
abla[\overset{ \rightarrow}{u}\overset{\rightarrow}{u}]=-\frac{1}{\rho(\varphi)}\
abla p +\frac{1}{\rho(\varphi)}\
abla\times\overset{\rightarrow}{\tau}+\overset{ \rightarrow}{g}-\overset{\rightarrow}{F}_{sf}+S_{y} \tag{11}\]  
where \(\overset{\rightarrow}{u}\) is the velocity vector, \(\varphi\) is the level-set function, \(\rho\) is the density, \(\overset{\rightarrow}{\tau}\) is the stress tensor, \(\overset{\rightarrow}{F}_{sf}\) is the force arising from surface tension effects, and \(S_{y}\) is the additional sink term induced by the solidification model. The coupled level-set and VOF (CLSVOF) model was specifically designed for two-phase flows, a popular interface-tracking method ([PERSON], 2000). It can precisely capture the gas-liquid interface and be characterized by good conservation properties. The interface was captured and tracked by the level-set function \(\varphi\), defined as a signed distance from the interface. Then, the continuity equation becomes:  
\[\frac{\partial\varphi}{\partial t}+\
abla\cdot(\overset{\rightarrow}{u} \varphi)=0 \tag{12}\]  
The \(\overset{\rightarrow}{F}_{sf}\) in Equation (5) specifically represents capillary forces and is written as:  
\[\overset{\rightarrow}{F}_  
and its value was between 0 and 1. The effect due to the reduced porosity in the mushy area was incorporated in the momentum equation sink term Sy, expressed as:  
\[S_{y}=\frac{\left(1-\beta\right)^{2}}{\left(\beta^{3}+\varepsilon\right)}A_{ mushy}\overset{\rightarrow}{u} \tag{16}\]  
where \(\varepsilon\) is a small number (0.001) introduced to prevent division to zero. Amush is the mushy zone constant, which reveals the degree of difficulty of the liquid solidification process. It was set as 108, based on [PERSON] et al. ([PERSON] et al., 2014 [37])  
For the enthalpy-porosity technique, the energy equation is expressed as follows:  
\[\frac{\partial}{\partial t}[\rho H]+\
abla\cdot[\rho\overset{\rightarrow}{u} H]=\
abla\cdot\left(\lambda\cdot\
abla T\right)+S_{h} \tag{17}\]  
where \(\lambda\) is the thermal conductivity and \(H\) is the total enthalpy per control volume. \(H\) is computed as the sum of the sensible enthalpy \(h\) and the latent heat \(\triangle H\), expressed as:  
\[H=h+\Delta H \tag{18}\]  
where \(h\) is defined as \(h=h_{ref}+\int_{T_{ref}}^{T}c_{p}dT\), \(h_{ref}\) is the sensible enthalpy at the reference temperature \(T_{ref}\), and \(c_{p}\) is the specific heat at constant pressure. The latent heat \(\triangle H\) is defined as \(\beta\times L\), where \(L\) is the latent heat of the material. The volume of all the droplets in this study was 10 \(\mu\)L. When using the solidification/melting model, there are two important parameters. One is the mushy zone parameter, which affects the phase transition rate and the accuracy of the calculation during the freezing process. Another important parameter is \(\Delta T=T_{liquid}-T_{solid}\), which is the difference between the liquid temperature and the solid temperature, which affects the stability and the thickness of the phase interface during icing calculations. Therefore, it is very important to choose appropriate parameters for a simulation of the droplet freezing process. This study used the test value as a reference for verification; the final \(A_{mesh}\) value was \(B\), and \(\Delta T\) was 0.5 \({}^{\circ}\)C.  
In the method of verifying the numerical simulation calculation, according to the experimental measurements, the bottom diameter of the droplet at \(-15\) \({}^{\circ}\)C was 4.4 mm, and the bottom diameter of the droplet at \(-5\) \({}^{\circ}\)C was also 4.4 mm. The droplet contour model for the numerical simulation was established by pushing back the contour parameters of the droplet by the Young-Laplace equation.