# A New Ice Accretion Model for Aircraft Icing Based on Phase-Field Method
## 2 Materials and Methods
### Two-Dimensional Form of Governing Equation
\(\dot{m}_{exp}\) is the evaporation rate, which is obtained by analogy with the convection heat transfer term [17],  
\[V=s\cdot h,V_{w}=s\cdot h_{w},V_{i}=s\cdot h_{i}, \tag{2}\]6Among them, \(Cp_{a}\) is the specific heat capacity of air at constant pressure, which can be taken as \(1003.5\) J/\((\text{kg}\cdot\text{K})\); \(\rho_{a}\) is the air density; \(R_{v}\) is the gas constant of water vapor, which can be taken as \(461\) J/\((\text{kg}\cdot\text{K})\); _Lew_ is the analogy criterion number of heat and mass transfer, which is related to the content of liquid water in the air, and can be taken as \(1\) in general icing conditions; \(e_{s}(T_{w})\) and \(e_{s}(T_{\infty})\) are the saturated steam pressure corresponding to the water temperature \(T_{w}\) of the control volume and the ambient temperature \(T_{\infty}\), respectively [20]. When \(T_{w}\geq T_{m}\),  
\[ \begin{split} e_{s}(T_{w})&=611.011+44.481(T_{w}-T_ {m})+1.419(T_{w}-T_{m})^{2}\\ &=+0.0239(T_{w}-T_{m})^{3}+1.744\times 10^{-4}(T_{w}-T_{m})^{4}, \end{split}  \tag{17}\]  
and when \(T_{w}<T_{m}\),  
\[ \begin{split} e_{s}(T_{w})&=609.603+49.495(T_{w}-T _{m})+1.739(T_{w}-T_{m})^{2}\\ &=+0.031(T_{w}-T_{m})^{3}+2.292\times 10^{-4}(T_{w}-T_{m})^{4}, \end{split}  \tag{18}\]  
\(T_{m}\) is the the phase transition equilibrium temperature of water, \(T_{m}=273.15\) K.  
The energy equation is given by:  
\[\frac{\partial(\rho hH)}{\partial t}+\
abla\cdot(\rho hH\textbf{U})=\dot{Q}_{ rad}+\dot{Q}_{imp}+\dot{Q}_{cep}+\dot{Q}_{htc}+\dot{Q}_{cond}. \tag{19}\]  
Among them, \(H\) is the internal energy of the mixture, which can be obtained as:  
\[H=CpT+(1-C)L=(1-C)Cp_{w}T_{w}+CCp_{i}T_{m}+(1-C)Lf, \tag{20}\]  
where \(Cp_{i}\) is the specific heat capacity of ice and water, respectively. \(\dot{Q}_{rad}\) is the radiation heat flux of the water film. \(\dot{Q}_{imp}\) is the collected water energy and includes the internal energy and kinetic energy of the droplet. \(\dot{Q}_{cep}\) is the energy density taken away by evaporation or sublimation on the surface of the water film. \(\dot{Q}_{htc}\) is the convective heat transfer density between the water film and the air flow. \(\dot{Q}_{cond}\) is the energy transmitted by the aircraft skin. The unknown terms are specified in the form of [22]:  
\[\dot{Q}_{rad}=\sigma_{r}(T_{w}^{4}-T_{\infty}^{4}), \tag{21}\]  
\[\dot{Q}_{imp}=\dot{n}_{imp}[Cp_{w}(T_{w}-T_{\infty})-\frac{u_{d}^{2}}{2}], \tag{22}\]  
\[\dot{Q}_{cep}=\dot{n}_{cep}Le, \tag{23}\]  
\[\dot{Q}_{htc}=htc(T_{w}-T_{rec}). \tag{24}\]  
Figure 2: Diagram of mass and energy transfer in the control volume. As shown in Figure 3, the water film flow can be simplified as a Couette flow under the action of air shear stress, and the velocity of the water film can be approximately described in the form of a linear distribution [21]:  
\[\mathbf{U}_{w}(s,y)=\frac{y}{\mu_{w}}\tau_{a}(s), \tag{25}\]  
where \(s\) is the arc direction of the surface and \(y\) is the normal direction of the surface. Then, the average flow rate of the water film can be defined as:  
\[\overline{\mathbf{U}}_{w}(s,y)=\frac{h_{w}}{\mu_{w}}\tau_{a}(s), \tag{26}\]  
in which \(\tau_{a}(s)\) is the shear stress of airflow on the water film. In icing simulations, the convective heat transfer process \(htc\) has a great influence on the ice shape, and is the key to solving the convective heat transfer process. At present, most icing simulations use the boundary layer integral functions proposed by LEWICE, which consider surface roughness and velocity variation to solve the convection heat transfer coefficient [11]. \[htc=\left\{ \begin{array}{ll}0.296\cdot\frac{\lambda_{s}}{\mu_{s}^{3.5}}(V_{e }^{-2.87}\int_{0}^{s}V_{e}^{1.87}ds)^{-0.5},&\text{Re}_{k}<\text{Re}_{c}\\ \frac{C_{f}/Z}{Prt+\sqrt{C_{f}}/St_{k}}\cdot\rho_{a}\cdot V_{e}\cdot Cp_{a},& \text{Re}_{k}\geq\text{Re}_{c}\end{array} \right. \tag{27}\]  
where \(\lambda_{s}\) is the thermal conductivity of air. \(v_{a}\) is the kinematic viscosity of air. \(V_{e}\) is the velocity of air at the outer boundary of the boundary layer.