Two-Dimensional Simulation of Frost Formation on the NACA0012 Airfoil under Strong Convection by Using the P-VOF Method

: A newly developed frosting simulation method, pseudo-VOF (p-VOF) method, has applied to simulate the dynamic frost formation on the NACA0012 airfoil under strong convection. The p-VOF method is a pseudo volume of fraction simulation method of the multiphase flow with phase change. By solving a simplified mass conservation equation explicitly instead of the original volume fraction equations in CFD software, the efficiency and robustness of calculation are greatly improved. This progress makes it possible to predict a long-time frost formation. The p-VOF method was successfully applied to the simulation of dynamic frosting on the two-dimensional NACA0012 airfoil under strong convection conditions with constant frost physical properties. The simulation result shows that the average thickness of the frost layer increases, and the frost bulges and flow separation appear earlier, when the airfoil surface temperature decreases or the air humidity increases. The frost bulges and flow separation appear earlier, when the air velocity is faster, the growth rate of the frost layer at the early stage is greater, but the final frost layer is thinner.


Introduction
When The aircraft cruise at high altitude where the ambient temperature is well below freezing, the fuel in the wing tanks could be cooled to below freezing point [1]. If the aircraft flies into a warm, humid airspace immediately, the water vapor in the air may desublimate on the surface of the wing. This frost phenomenon usually occurs when the aircraft descends from the cruise altitude to land. The frosting on the wing is different from the icing on the wing. The icing on the wing is the process of liquid water turning into solid ice when the supercooled water drops hit the wing [2][3][4], while the frosting on the wing is the process of water vapor desublimate into solid frost on the cold wing. The wing frosting has not received enough attention and there is a few public researches of the desublimation frosting on the wing [1,5], although there have been many studies on aircraft icing.
The frost on the wing surface will change the roughness of the skin and the shape of the wing, resulting in change of the lift-drag ratio even the airflow separation behavior [6], thereby affecting the aerodynamic performance of the aircraft. For energy saving and flight safety requirements, the importance of wing frosting studies has become increasingly prominent. It is necessary to predict the frost formation on the wing in different flight environments and evaluate its possible adverse effects firstly. The wing frosting is the problem of the desublimation frost on the cold surface under strong convection conditions.
In order to predict the frosting behavior on cold surface, several numerical models have been developed, which can be divided into three categories. The first type is onedimensional model, which uses empirical correlations obtained from experiments to calculate the heat transfer coefficient and mass transfer coefficient on the frost surface [7][8][9]. Nevertheless, the one-dimensional models cannot reflect the characteristics of complex surface shape, and are usually applied to the researches of frosting on cold plate [10]. Moreover, this type models do not solve the coupled flow and heat transfer problems in air and frost, so it cannot be used to predict the frost formation under changing flows with frosting. The second type of models solves the flow, heat and mass transfer within the air, the heat and mass transfer within the frost layer, which is regarded as a solid region [11][12][13]. This type models solves the governing equations in air domain and frost domain respectively. The third type of models uses the multi-phase flow method to solve the flow, heat and mass transfer simultaneously in humid air and frost layer with the same set of governing equations. Current researches often adopt the VOF model [14,15] and the Eulerian model [16,17]. The multiphase flow model requires the least empirical correlations and assumptions, so it is more suitable for the simulation of frosting in actual working conditions. However, those researches only studied the simulation of frosting under natural convection or low velocity flow conditions.
According to the frosting experiment on the cold plate under strong convection [18], it was found that when the flow velocity is faster than 10m/s that was much faster than the flow velocity in the refrigeration systems, the frost layer is compact. Under strong convection, the new desublimated water vapor mainly contributes to the increase of the frost thickness, rather than increases the density of frost layer. However, the desublimated water vapor increases both the thickness and the density of the frost layer under low airflow velocity [10]. As the frosting progresses, the frost surface temperature gradually increases, and when the water vapor content in the humid air is equal to the saturated water vapor content at the frost surface temperature, the water vapor in the humid air will stop to change into frost. If the water vapor content in the humid air is lower than the saturated water vapor content at 273.15K, the frost surface temperature cannot reaches to 273.15K, and no liquid water will appear on the frost surface. This kind of desublimation frosting processes under these humidity conditions was called as dry mode frosting [18]. This paper developed a new multi-phases flow method, p-VOF method, to simulate the dry mode frosting on the NACA0012 airfoil under strong convection conditions.

P-VOF method for frosting simulation
The p-VOF method is a pseudo volume of fraction method of multiphase flow simulation, which solves a simplified mass conservation equation to obtain the volume fraction change of frost phase explicitly instead of the original volume of fraction (VOF) equations in CFD software to improve the calculation efficiency and robustness. This progress makes it possible to predict a long-time frost formation. This method does not solve the volume fraction equation in original VOF method, so it is called pseudo-VOF (p-VOF) method.
According to the observation of experiments [18], the assumptions for simulation are as follows:  The humid air is an incompressible Newtonian fluid with constant properties.  The phase change between water vapor and frost occurs only on the frost surface (or the airfoil surface at the initial moment).  When the frost surface temperature is lower than the freezing point, if the water vapor content is greater than the saturated water vapor content on the frost surface, water vapor desublimates into frost, and vice versa.  The mass transfer inside the frost is neglected because the frost layer formed under strong convection is dense. Humid air is set to be the primary phase and frost is the secondary phase. Humid air is multi-component which composes of dry air and water vapor. The VOF model calculate the mass, momentum, energy transport and water vapor diffusion by equations as shown in Eqs. (1), (2), (3) and (4) respectively, where  is the density, p is the pressure, u v is the velocity vector, t is the frosting time, α is the volume fraction, Sm is the mass source/sink term due to phase change. subscript i is a or f , representing the humid air phase or the frost phase. The sum of volume fractions of two phases is 1, as shown in Eq. (5), and F v is the body force or momentum source, μ is the viscosity and keff is the effective thermal conductivity of the mixture, , μ and keff can be calculated by sum of each phase weighted by volume fraction. T is the temperature and Sh is the heat released or absorbed during phase change. E is the internal energy of the two-phase mixture, which is mass-averaged by phase fraction. w is the mass fraction of the water vapor in humid air, 2 HO D is the diffusion coefficient of water vapor in humid air.
Because the frost phase does not move after it generated and the frost surface changes only when desublimation or sublimation occurs, the convective term (the second term on the left) of volume fraction in Eq. (1) can be ignored. Eq. (1) can be simplified to At the phase interface, the source term can be calculated by the mass transfer rate from water vapor in humid air to frost, which is determined by the kinetic theory [19], where σ is the coefficient of correction and 0.03 is suggested for water [19]; v is the volume of the mesh cell. Rw is the water vapor gas constant, which equals to 461.5 J/(Kg·K); ρv is the density of water vapor, ρs is the saturated density of water vapor at the phase interface corresponding to T at the phase interface. Source term Sh in energy Eq. (3) is the latent heat released or absorbed due to phase change, which is not zero only at the phase interface where γ is the latent heat of phase change.
The transient two-dimensional frosting process is simulated by using ANSYS FLUENT with VOF multiphase flow model. The flow chart of simulation is shown in Fig. 1.

Fig. 1 Flow chart of simulation
The PISO algorithm is chosen to couple the pressure and velocity. The pressure is discretized by the PRESTO! method. At each time step, the change of frost volume fraction at the phase interface is determined only by explicitly solving which is transformed from Eq. is set to 0.
The calculation of source terms, searching for phase interface cells, and the evolution of volume fraction are implemented by employing UDFs (user-defined functions), which are hooked to ANSYS FLUENT. It should be noted that since the original VOF equation of FLUENT is not solved, the software has data transmission problems in the assignment of the thermal conductivity and viscosity of the mixture. Therefore, it is necessary to use UDFs to assign the thermal conductivity and viscosity of the mixture to the humid air phase in order to simulate frosting correctly.
In order to verify the applicability of the p-VOF method for dynamic desublimation frosting on airfoil under strong convection conditions, simulations on NACA0012 are carried out under different air velocity, air humidity and airfoil surface temperature conditions with constant frost physical properties. The frost density uses 560 kg/m 3 and frost thermal conductivity uses 0.2 W/m·K. The chord length of airfoil is 10mm. The simulation conditions are in Table 1. The maximum Reynolds number of these conditions is 4×10 4 , so the laminar model is used.  Fig. 2. The half model is used to reduce the burden of computation. The grids is refined near the airfoil surface. The length (X direction) of the mesh cell on the surface of airfoil is 1 mm. A boundary layer is attached to the airfoil surface with first layer of 0.01 mm and increasing ratio of 1.007 (Y direction). The airfoil surface is set to isothermal wall condition.

Dynamic desublimation frosting on airfoil
The baseline condition is set for the simulation: the airfoil surface temperature 250 K; the air velocity 10 m/s; the air temperature 310K and the relative humidity 7%. The frosting duration is 30 minutes. The average thickness curve of the frost layer is shown in Fig. 3, and the average thickness of the frost layer at the end of 30 minutes is 0.43mm. Under the baseline condition, the frost morphology at different instant is shown in Fig. 4. The flow direction is from left to right. The frost layer thickness is relatively uniform on the airfoil in the early stage of frosting. Only the frost layer close to the leading edge of the wing is slightly thicker than the frost on the other parts of the airfoil. As the frosting progresses, the frost layer gradually thickens. At the 15th minute, a slight bulge appears in the front of the frost layer. Then the bulge becomes more and more obvious, and some new smaller bulges continue to form on the frost surface afterwards.

The influence of air velocity on frosting
The dynamic frosting behaviors on the two-dimensional airfoil under different air velocities are simulated and compared. The variation of the average thickness of the frost layer under different air velocities is shown in Fig. 9. The higher the air velocity, the faster the growth rate of the frost layer at early stage, but the earlier the frosting rate begins to decrease significantly, and the thinner the final average thickness of the frost layer will be.

Fig. 9 Variation of the average thickness of the frost layer under different air velocities
The change of the frost layer morphology on the airfoil surface under different flow velocity conditions are shown in Fig. 10. The frost layer morphology difference under different air velocities are quite obvious. The first bulge position of the frost layer under 50m/s is more rearward than that under 10m/s, as can be seen from the comparison of the frost layer morphology at the 5th minute in the figure. The morphology of the frost layer at 50 m/s did not change significantly after 10 minutes, while the morphology of the frost layer kept changing at a lower speed of 10 m/s.

The influence of air humidity on frosting
The frosting behaviors on the airfoil under different air humidity were simulated successfully. The change in average frost thickness with time under different humidity conditions is shown in Fig. 13. It can be seen from the curves that as the air humidity increases, the growth rate of the frost layer increases, and the average thickness of the frost layer increases in the same moment. The morphology change of frost layer with time under different air humidity conditions is shown in Fig. 14. Higher air humidity leads to greater frosting rate, and the thicker frost layer, which is consistent with the change of the average frost thickness in Fig. 13 . Moreover, the higher the huimidity is, the earlier the obvious frost bulges appears, and the frost morphology change is more complicated. The velocity vector field varying with time under different air humidity are shown in Fig. 15. From the vortex sequence in the velocity vector field figure, it can be seen more clearly that the higher the incoming flow humidity, the earlier the flow separation occurs, which caused by the faster frost formation rate on the airfoil surface. The change of temperature field with time under different humidity conditions are shown in Fig. 16. The higher the air humidity, the earlier the occurrence of the lowtemperature air masses caused by the vortex, which is a result of the morphology change of the frost layer. The frosting simulations using the p-VOF method were carried out successfully under different airfoil temperatures. Fig. 17 shows the average thickness curves of the frost layer under different airfoil surface temperatures. The lower the airfoil surface temperature, the faster the frost layer grows, and the thicker the frost layer is in the same moment. The development of the frost layer morphology under different airfoil surface temperatures is shown in Fig. 18. When the airfoil surface temperature is lower, the frost formation rate is greater, and the frost layer is thicker on the airfoil. Moreover, when the temperature of the airfoil surface is lower, the frost bulges appears earlier and more obvious, and the bulges locate more close to the leading edge. After the bulges formed under the airfoil surface temperature of 230K, there are some holes wrapped in the frost layer at the position behind the first bulge location. This may be due to the greater frost formation rate and the relative lower water vapor replenishment by diffusion under such frosting condition, leading to the advantageous concentration effect of desublimation frost, which may be similar to porous structure of frost formed under natural convection and low speed flow conditions. The influence of the frost layer morphology caused by different airfoil surface temperatures on flow separation can be seen clearly from the velocity vector field (shown in Fig. 19). It can be seen that the flow separation and the generation of vortices appear earlier due to the faster frost formation when the temperature of the airfoil surface is lower. Fig. 20 shows the temperature field under different airfoil surface temperatures. When the airfoil temperature is lower, the low temperature air zones caused by the airfoil flow vortex appears earlier. It can be found by comparing the temperature fields under different airfoil surface temperatures that the location of the frost layer bulges move forward, as well as the position where the vortices begin to separate when the surface is lower.

Conclusion
Through the above simulation, the following conclusions can be obtained: 1) The p-VOF method can be applied to conduct a stable and efficient numerical solution of dynamic desublimation frosting on the airfoil under different conditions.
2) The dynamic frosting process on the airfoil can be obtained clearly and vividly by using the p-VOF method. This method can be used to predict the frosting behavior on the airfoil.
3) The simulation result of frosting on the NACA0012 shows that the average thickness of the frost layer increases, and the frost bulges and flow separation appear earlier, with decrease of airfoil surface temperature or the increases of air humidity. However, the frost bulges and flow separation appear earlier, when the air velocity is faster, the growth rate of the frost layer at the early stage is greater, but the final frost layer is thinner. Figure 1 Flow chart of simulation Mesh for frosting simulation on airfoil The morphology of the frost layer on the airfoil varying with time    Water vapor content change in the air with the frosting on airfoil