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, p-VOF method, is applied to simulate the dynamic frost formation on the NACA0012 airfoil under strong convection. The p-VOF method is a pseudo VOF method of the multiphase flow simulation 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 cruises 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 desublimating 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 first. 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 a cold surface, several numerical models have been developed, which can be divided into three categories. The first type is the one-dimensional 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 research of frosting on cold plates [10]. Moreover, this type of models does 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][14]. This type of models solves the governing equations in the air domain and frost domain respectively, and couples them through the heat and mass transfer relation at the phase interface, and uses dynamic meshes technique to move the phase interface. The third type of models uses the multi-phase flow method to solve the flow, heat and mass transfer simultaneously in the humid air and frost layer with the same set of governing equations. Current researches often adopt the VOF model [15,16] and the Eulerian model [17][18][19][20]. This multiphase flow model can couple the heat and mass transfer between the phases and move the phase interface without modifying the mesh. 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 [21], it was found that when the flow velocity is faster than 10 m/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 the 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 changing into frost. If the water vapor content in the humid air is lower than the saturated water vapor content at 273.15 K, the frost surface temperature cannot reach to 273.15 K, 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 [21]. 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 VOF 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 fraction 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 equations in original VOF method, so it is called pseudo-VOF (p-VOF) method.
According to the observation of experiments [21], 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 is composed of dry air and water vapor. The VOF model calculates 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 * is the velocity vector, t is the frosting time, α is the volume fraction, and S m 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 * is the body force or momentum source, μ is the viscosity and k eff is the effective thermal conductivity of the mixture. ρ, μ and k eff can be calculated by sum of each phase weighted by volume fraction. T is the temperature and S h 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, and D H 2 O is the diffusion coefficient of water vapor in humid air.
Because the frost phase does not move after it is 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 [22], where σ is the coefficient of correction and 0.03 is suggested for water [22]; v is the volume of the mesh cell. R w is the water vapor gas constant, which equals to 461.5 J/ (Kg·K); ρ v is the density of water vapor, and ρ s is the saturated density of water vapor at the phase interface corresponding to T at the phase interface. Source term S h 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. 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. (6), where Δt is the time step, and Δα f is the variation of the frost volume fraction during Δt. By explicitly solving Eq. (9) instead of the original volume fraction equations of FLUENT (Eq. (1)), the efficiency and robustness of calculation can be greatly improved. If the phase interface temperature T > 273.15 K, there is no growth of frost, and Δα f is set to 0. The time step is 0.5 s for the simulation and the maximum iteration in every time step is 20. The residual is set to 10 − 3 for convergence of velocity, continuity and vapor diffusion, and 10 − 6 for convergence of energy.
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 volume fraction 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.

Frosting simulation on the two-dimensional airfoil
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 is set to 560 kg/m 3 and the frost thermal conductivity is set to 0.2 W/m·K. The chord length of the airfoil is 10 mm. The simulation conditions are in Table 1. The maximum Reynolds number of these conditions is 4 × 10 4 , so the laminar model is used.
The mesh of the computation region is shown in Fig. 2. The half model is used to reduce the burden of computation. The grid is refined near the airfoil surface. The length (X direction) of the mesh cell on the airfoil surface 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. The flow direction is from left to right. The right boundary is set to the velocity inlet condition and the right boundary is set to the pressure outlet. The top and bottom boundaries are set to the symmetry 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 310 K and the relative humidity 7%. The frosting duration is 30 min. The independence of mesh is verified first. Two meshes are used for comparison. One mesh is refined by 2 times, and the other mesh is coarsened to half. The average thickness curve of the frost layer is shown in Fig. 3. The average thickness of the frost layer is calculated by dividing the cross-sectional area of the frost layer by the perimeter of the airfoil. The relative error of the average thickness of the frost layer of different meshes is within 5%, indicating that the simulation results are insensitive to the mesh used. The mesh is suitable for the simulation of frosting by using the p-VOF. The average thickness of the frost layer at 15 min is 0.32 mm, and the average thickness at 30 min is 0.43 mm. The average thickness only increases 0.11 mm in the last 15 min.
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 velocity contour with streamlines (upper) and temperature contour (lower) varying with time are shown in Fig. 5. At the 5th minute, the flow on the airfoil surface still attaches to the surface. At the 10th minute, flow separation has occurred near the rear of the airfoil, and vortices appear. At the 20th and 30th minute, the flow separation and vortices on the airfoil surface are more obvious and regular. A series of small to large vortices can be clearly distinguished. At the 5th minute, the temperature boundary layer on the airfoil still maintains an attached flow. At the 10th minute, the temperature boundary layer at the rear of the airfoil thickens and develops some lowtemperature air zones in the shape of eddy current at the tail of the airfoil. At the 20th and 30th minute, from the position of the frost bulge, some clear vortices of lowtemperature air appear.
The moisture content change of the air near the frost layer is shown in Fig. 6. Before 10 min, there are obvious areas with low water vapor content near the frost layer, especially near the rear of the airfoil. After 20 min, there are only some low humidity areas on the leeward side of the frost layer bulges, indicating that the frost layer areas in other locations have almost reached the mass transfer balance.

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. 7. 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. The change of the frost layer morphology on the airfoil surface under different flow velocity conditions is shown in Fig. 8. The frost layer morphology difference under different air velocities is quite obvious. The first bulge position of the frost layer under 50 m/s is more rearward than that under 10 m/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 min, while the morphology of the frost layer kept changing at a lower speed of 10 m/s. There is a relationship of heat transfer balance in frosting, which is that the sum of the convective heat transfer from the air and the latent heat release from the phase change is equal to the heat transferred to the cold surface through the frost layer. As the frosting progresses, the frost layer becomes thicker, the frost surface temperature will rise, the convective heat transfer becomes smaller, and the phase change rate  becomes smaller. According to Assumption 3) (that is, Eq. (7)), when the saturated water vapor content corresponding to the frost surface temperature is equal to the saturated water vapor content of the incoming air, the frosting will stop. After this moment, the phase change latent heat release is zero, the convective heat transfer is equal to the heat transfer through the frost layer, and the frost layer thickness reaches the maximum. The faster the speed, the greater the heat transfer by convection and the greater the heat transfer from the frost layer when the heat transfer is balanced. The greater the heat transfer, the smaller the thickness of the frost layer, in the case of constant thermal conductivity of the frost layer. In addition, the faster the speed, the greater the frosting rate in the early stage of frosting. Therefore, when the air velocity is faster, the frosting will reach the balance earlier, that is, the frost layer will reach the maximum earlier.
The velocity contour with streamlines varying with time under different air velocities is shown in Fig. 9. It can be clearly seen from the figure that the flow on the frost surface has an obvious vortex sequence at the 5th minute under the incoming airflow velocity of 50 m/s.
The change in temperature field under different air velocities is shown in Fig. 10. Compared with the case of 10 m/s that maintains the attached temperature boundary layer at the initial stage of frosting, the flow on the surface of the frost layer has shown an obvious low-temperature vortices sequence at the 5th minute under 50 m/s condition, which indicates the separation of flow.

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. 11. 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 at the same moment. The morphology change of the frost layer with time under different air humidity conditions is shown in Fig. 12. Higher air humidity leads to greater frosting rates and thicker frost layers, which is consistent with the change of the average frost thickness in Fig. 11. Moreover, the higher the huimidity is, the earlier the obvious frost bulges appear, and the frost morphology change is more complicated.
The velocity contour with streamlines varying with time under different air humidity is shown in Fig. 13. It can be seen clearly that the higher the incoming flow humidity, the earlier the flow separation occurs, which is caused by the faster frost formation rate on the airfoil surface.
The change of temperature field with time under different humidity conditions is shown in Fig. 14. 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 influence of airfoil surface temperature on frosting
The frosting simulations using the p-VOF method were carried out successfully under different airfoil temperatures. Figure 15 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 at the same moment.
The development of the frost layer morphology under different airfoil surface temperatures is shown in Fig. 16. 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 appear earlier and more obvious, and the bulges are located more close to the leading edge. After the bulges are formed under the airfoil surface temperature of 230 K, 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 relatively lower water vapor replenishment by diffusion under such frosting conditions, 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 figures of velocity contour with streamlines (shown in Fig. 17). 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. Figure 18 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 appear earlier. It can be found by comparing the temperature fields under different airfoil surface temperatures that the location of the frost layer bulges moves forward, as well as the position where the vortices begin to separate when the surface temperature is lower.

Conclusion
Through the above simulations, 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 the decrease of airfoil surface temperatures or the increase 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.

Declaration
Competing interests