Abstract
Considering the transient heat and mass transfer process of the impinged water droplets during aircraft icing, an unsteady thermodynamic model was established to simulate the dynamic developments of the water film and the ice layer on aircraft surfaces. The unsteady model was discretized in an implicit scheme with a corresponding solution method. Icing simulations were performed for a NACA0012 airfoil, and the results show acceptable agreement with the data in the literature. Water film first appears near the stagnation point, and then, the film thickness increases, and the runback water region expands with time, affecting the icing rate, the surface temperature, and the ice type. The development of the water film is rapid, and the thickness and range of the film, along with the icing rate, reach a steady state in a short time. The stable characteristics obtained by the unsteady model are consistent with those of the Messinger steady model. Despite that the unsteady and steady models can obtain similar ice shapes in icing simulations, the dynamic developments of the water film and the ice layer should be considered at the initial stage of ice accretion or in the short-time icing simulations.
1. Introduction
Supercooled water droplets in clouds, of which the temperature is below 0°C but still unfrozen, may impinge on aircraft components and cause ice accretion on their windward surfaces [1]. The ice layer would change the aerodynamic configuration and affect the lift and drag, leading to a deterioration of the flight performance, maneuverability, and stability. Furthermore, it can even cause flight accidents [2]. Research on aircraft icing to accurately predict the ice shape could help to analyze the influence of the ice layer on aerodynamic performances [3] and guide the design of aircraft ice protection systems [4], which is very important for flight safety.
When supercooled water droplets hit the aircraft surface, the thickness of the formed water film and the height of the ice layer both increase dynamically with time, whereas the surface temperature changes accordingly [5, 6]. Aircraft icing is an unsteady heat and mass transfer process. Simulation model was firstly applied to predict the icing process in 1950s and has been widely used since then [7]. A series of mature icing analysis codes have been developed so far, including LEWICE [8] in the United States, FENSAP-ICE [9] in Canada, ONERA [10] in France, TRAJICE [11] in the United Kingdom, and CIRAAMIL [12] in Italy. During aircraft icing simulations, it is essential to analyze the mass and energy transfer process of the runback water flow and the phase transition, and establish a simplified thermodynamic model to obtain the icing rate and finally the ice shape. Various icing thermodynamic models with different simplifications have been developed and applied in the mentioned codes, of which the most popular ones are the Messinger model, the shallow-water model, and the Myers model.
The classic Messinger model [7] is the base of almost all the thermodynamic models for aircraft ice accretion and has been widely applied in the icing codes such as LEWICE. This model establishes the mass and energy balance equations of a control volume (CV) on the icing surface according to the conservation law. It makes the assumption that all the unfrozen water on the surface would flow to adjacent CVs, meaning that there is no resident liquid water (or water film) in any CV [13]. The Messinger model is a steady one, in which the icing rate, the surface temperature, and the runback water flux are constant with a stable growth rate of the ice layer. Therefore, it could not depict the water film development and the temperature change of the ice layer.
Based on the Messinger theory, a shallow-water model was established by FENSAP-ICE to consider the unsteady icing effects [14]. It is assumed that the unfrozen water can form a liquid film upon the aircraft skin. The thickness of the water film changes with time when it runs back along the surface and freezes with the heat loss to the external airflow. The film thickness and icing rate are obtained in each time step by the simultaneous solution of the momentum and energy conservation equations of the water film. In FENSAP-ICE, the air shear stress is the only factor that drives the water runback, and thus, the velocity of the water film is linearly distributed along the normal direction. Furthermore, the model focuses on the water film development, but the effects of the accumulated ice are ignored without considering the heat conduction and the temperature change of the ice layer.
Both the Messinger and shallow-water models assume that the aircraft surface is adiabatic. Myers [15] applied a one-dimensional (1D) heat conduction method to consider the freezing latent heat transferring to the external airflow through water film and to the inner skin through the ice layer. An unsteady icing model was then established based on the Stefan phase transition theory, classifying the icing process into glaze ice and rime ice. The original Myers model only considered the phase transition process, ignoring the water film flow [15]. Researchers then applied a flow model [16, 17] similar to the shallow-water one for this heat transfer theory and even extended the Myers model to the three-dimensional (3D) cases [18]. Since the heat transfer process in aircraft icing is comprehensively considered in the Stefan phase transition theory, the development and extension of the Myers model attracts more and more attention [19–21], although the solution cost is expensive. In those models, the latent heat transferred to the internal skin is determined by the ice thickness and the temperature of the substrate. Since the thermal conductivity of ice is very small, the heat transfer through the ice layer would decrease to a very low level when the ice thickness increases to a certain extent. Then, it can be considered that the freezing latent heat would not conduct to the skin, and the ice layer can be simplified as an adiabatic boundary, which is consistent with the Messinger and shallow-water models. As for the thin ice layer, the aircraft skin temperature plays a very important role [21], and the icing process is determined by the coupling of the Myers model and the heat conduction in the substrate [22]. However, the coupling solution would take lots of time, and few of those models for aircraft icing take this into account [15–21].
Generally, ice shape is determined by the mass flow, heat transfer, and phase transition processes of the supercooled droplets impinged on the aircraft surface. The dynamic developments of the water film and the ice layer are simplified to varying degrees in the mentioned icing thermodynamic models, with the research trend mainly from steady-state model to dynamic calculation. The latent heat conducted to the substrate has also been included. When icing model is becoming more and more complex, the cost and computing resources of the solution process are increasing. However, the prediction accuracy of ice shape is equivalent or only slightly improved, compared with the classical models. Furthermore, the icing codes that are widely used in practice are still LEWICE with Messinger model and FENSAP-ICE with shallow-water model. Therefore, it is necessary to analyze the factors affecting the icing process and select a suitable model with reasonable simplification at the least cost.
The present paper focuses on the dynamic developments of the water film and the ice layer, as well as the temperature change during aircraft icing process. The freezing latent heat transferred to the skin described in the Myers type models is not considered, since the substrate temperature changes with time and is difficult to determine. Based on the basic physical process of aircraft icing, the Messinger steady model and an unsteady thermodynamic model similar to the shallow-water one are derived separately according to different assumptions, as illustrated in the next section. The icing simulation procedure is presented in Section 3. The validation of the unsteady model and the mechanism analysis of the water film and ice layer developments are carried out in Section 4, along with the comparisons of the Messinger steady model and the unsteady model.
2. Mathematical Model
When there are supercooled water droplets impact or runback water flows into a certain control volume (CV) on the icing surface, the water would be in the form of solid ice, liquid water, or both, if the evaporation rate is less than the water inflow rate. Without considering the heat conducted to the inner substrate, the heat and mass transfer of the total liquid water and solid ice in the CV is analyzed, as shown in Figure 1. In addition to the convection terms caused by the flow velocity of the water , the water quality and energy changes are affected by various mass and heat sources [23]. The mass sources in the CV consist of the impinged water flow rate and the water evaporation rate , while the energy sources include the heat flow rate of the impinged droplets , the evaporative heat flow rate , and the convective heat flow rate .
Assuming that the thicknesses of the water film and the ice layer are both small on the icing surface, the temperature of the whole water is uniform in the CV and considered equal to the value of the outer skin surface. In addition, the aircraft skin is regarded as adiabatic, and the lateral heat conduction between the CVs along the icing surface is ignored, meaning that there is no diffusion term in the water energy Equation (23). As for the convection terms, the mass flow rate of the runback water entering the current CV is with the heat flow rate of , and the flowing out water rate is with the heat flow rate of . Therefore, the mass and energy equations of the whole liquid water and solid ice in the CV can be obtained by the conservation law: where is the total amount of the water in the CV, is the time, and is the total energy of the water in the CV.
The first term of the conservation equations is the unsteady term, the second one in the brackets is the convective term, and the ones on the right side are the mass and energy sources. Firstly, the convective term is a function of water film thickness, flow velocity, and surface temperature. Secondly, the liquid water and the solid ice might coexist in the CV, leading to a complex heat and mass transfer process of the total water. Therefore, the momentum equation of the water cannot be established as the parameters in the conservation equations are difficult to determine. The mass and heat transfer model of the whole water cannot be solved directly [23], and additional assumptions are needed to simplify and solve the conservation equations.
2.1. Simplification with Messinger Theory
As described in the introduction, the Messinger model assumes that all the unfrozen water would flow downstream, and the icing process is steady. The dynamic changes of the surface temperature and icing rate are not considered, and the unsteady terms can be written as: where is the mass of ice layer in the CV, is the icing rate, is the energy of the ice layer, and is the latent heat flow rate of water freezing.
In addition, since there is no resident liquid water in the CV, the convective term can be determined by the inflow and outflow relationship between adjacent CVs. The velocity and momentum equation of water is not needed. In 2D cases, there are only one inflow and one outflow boundaries for a CV, and the conservation equations of Equation (1) become:
It can be seen that Equation (3) is consistent with the equations in the icing codes such as LEWICE [8]. The model could be solved with the help of the definition of the freezing fraction by supplementing the relationship between the CV temperature and the water/ice phase state. The freezing fraction is defined as the ratio of the ice quality to the total water in the CV [13]:
Therefore, when , liquid water and solid ice coexist in the CV, and the CV temperature is at the phase transition value. When , there is no water freeze, and . And when all the inflow water freezes, , and . Additionally, since there is no runback water entering the stagnation point, the solutions of the above conservation equations can be initiated at the stagnation CV and continue in turn in the surface direction [8].
2.2. Unsteady Model with Developments of Water Film and Ice Layer
Considering the dynamic developments of the water film and the ice layer in Equation (1), an unsteady thermodynamic icing model similar to the shallow-water one is established in this paper. Assuming that the unfrozen water can always form a continuous water film with the thickness of , part of the liquid water would flow under the effect of shear stress, pressure, and gravity, while the rest would stay as the film form in the CV. Since the solid ice layer would not move, the movement of the total water in the CV is determined by the water film, and its momentum equation is established based on the law of motion under force, as expressed as where is the density of liquid water, is the velocity vector of water film, is the shear stress tensor, is gravitational acceleration, and is the pressure consisting of the air pressure and the surface tension: where is the surface tension coefficient of the water film and is the curvature. The water film is quite thin (thickness of 10-4 m or less), and its thickness gradient and the curvature are also very small [9]. Therefore, the effect of surface tension is much less than the shear stress and air pressure and is negligible. In addition, the velocity of water film is always small (10-1 m/s or less) in ice accretion conditions, and the water flow can be regarded as laminar flow [24]. At the normal direction of the icing surface, there is: where is the viscosity of liquid water.
For a thin water film, the velocity and its gradient along tangential direction are very small, and the momentum diffusion is negligible [24]. In addition, the thin water film would rapidly develop into a stable flow. Therefore, the momentum equation of the water film can be simplified as:
The bottom of the water film is in contact with the skin surface or ice layer, which is a nonslip boundary, and the top is under the effect of shear stress . There is
In a 2D case with the tangential direction of , the momentum equation becomes: where is the velocity at the direction and is the -component of the gravitational acceleration. By integrating the above equation, the velocity distribution of water film in the normal direction is obtained:
It can be seen from Equation (11) that when the pressure and gravity are ignored, the velocity distribution would be linear correlated with , which is consistent with the shallow-water modelin FENSAP-ICE [24]. If considered, the velocity distribution would be squarely correlated with , and the average velocity of the water film is derived by integrating Equation (11):
In a 3D case, the velocity along the other tangential direction on the skin surface can be derived in a similar way. Then, the average velocity would be applied in the conservation equations in Equation (1), and the mass equation of the water film is expressed as: where is the bottom surface area of the CV.
The mass change of the ice layer depends on the icing rate in the CV:
Since this unsteady model considers the development of the ice layer, the icing rate varies with time, which is different from the Messinger model in Equation (2) where the icing rate is constant. Meanwhile, the thickness of the ice layer can be obtained by: where is the density of solid ice.
In addition, since the liquid water is incompressible, Equation (13) becomes:
It can be found that the above mass conservation equation of the water film is consistent with that of the shallow-water model in EENSAP-ICE [24]. However, the present unsteady model adds Equation (14) to depict the dynamic growth of the ice layer on the icing surface.
As for the energy equation of the whole water (Equation (1)), by assuming that the enthalpy of the liquid water at 0°C is zero, and applying the average velocity to the energy conservation equation, there is: where and are the specific heats of the liquid water and the solid ice, is the surface temperature (°C), and is the latent heat of water freezing.
Considering that the mass and temperature of the ice layer vary with time when ice grows, the energy equation become:
Different from the energy equation of the shallow-water model in EENSAP-ICE [24], the present unsteady model includes the heat flow rate caused by the temperature variation of the ice layer, as shown by the last term of the above equation. Since this energy conservation equation takes the solid ice and the liquid water into account, it is applicable to all the water states on the icing surface, i.e., rime ice, glaze ice, and liquid water film state.
To sum up, when the physical property parameters of water and ice are constant in 2D cases, besides the momentum equation of Equation (12) and the dynamic mass equation of the ice layer (Equation (14)), the unsteady thermodynamic model also includes the transient mass and energy conservation equations of the water film:
The source terms in the above equations are calculated as follows.
The impinging water flow rate can be calculated by: where is the droplet collection efficiency, is the air velocity in the free stream, and is the liquid water content.
The heat flow rate of impinged droplets consists of the kinetic energy and the water enthalpy: where is the ambient temperature (°C).
The water evaporation rate is calculated by the Chilton-Colburn analogy method of mass and heat transfer [23]: where is the convective heat transfer coefficient, is the specific heat of air, is the Prandtl number, is the Schmidt number, and are the molecular weights of the water vapor and the air, is the saturated vapor pressure at the temperature , is the vapor pressure at the edge of boundary layer, and is the static pressure there.
With water evaporation, the energy flow rate during evaporation can be obtained by: where is the latent heat of sublimation, is the specific heat of ice, and is the latent heat of evaporation.
In addition, the convective heat flux is caused by the temperature difference and can be calculated by [23]: where is the reference surface temperature calculated with an adiabatic wall (°C).
2.3. Discretization and Solution of Unsteady Model
An implicit finite difference scheme is applied in the discretization of the unsteady terms in the mass and energy conservation equations: where superscript denotes the time step .
The convection terms are discretized in a forward difference scheme: the previous provides the data of the velocity, thickness, and temperature of the water film entering the current , and the current provides those of the water film runback to the downstream .Then, there is:
In addition, the CV parameters in the source terms of the conservation and movement equations and the icing rate should use the data of the current time step and CV. Then, the discretization of the mass conservation equation is obtained:
Meanwhile, the discrete form of the energy conservation equation becomes:
Note that the parameters without superscript are the data in the current time step, and those with no subscript are of the current CV. The parameters in Equation (12) are also the data of the current time step and CV.
The solution of the unsteady equations starts from the stagnation point at when the initial condition is known and proceeds downstream along the icing surface. When calculating the , the parameters of the inflow water are known from those of the , according to the inflow-outflow relationship determined by the velocity direction of the water film. The solution of each CV on the surface can be therefore obtained at and then comes to the next time step, of which the calculation also starts from the stagnation point.
It can be seen that there are 5 unknown parameters: , , , , and in Equation (12), Equation (14), Equation (27), and Equation (28). Another constraint equation is needed here. The freezing fraction described in the Messinger model is then extended to the unsteady thermodynamic model, which is redefined as the ratio of the ice mass to the total mass of the liquid water and ice layer in the CV at the current time step:
In a specific solution iteration for a CV, assume the surface temperature , solve the above equations to obtain the icing rate , and then, check:
If , the CV is under an ice-water mixed phase state. The assumption of is rational, and the calculated and are the right icing rate and water film thickness for the CV
If , all the water is liquid phase. Set , and resolve the conservation equations to obtain the actual temperature and the water film thickness
If , all the water freezes and no liquid water film runback to adjacent control volume. Then, set , and resolve the conservation equations to obtain the actual temperature and the icing rate
3. Simulation Procedure
Numerical simulation of aircraft icing always consists of the solutions of the airflow field, droplet impingement characteristics and icing thermodynamic model, and the geometric reconstruction for the ice shape. The present paper focuses on the heat and mass transfer process when the water film and the ice layer develop on the icing surface. Therefore, a single-step method is applied where the air and droplet flow fields are regarded as steady. Based on those results, the solution of the thermodynamic model is implemented by the commercial CFD software ANSYS FLUENT-19.1 with its user-defined functions (UDFs) [25].
The flow and temperature fields of the air around the aircraft are obtained by solving the Reynolds-averaged Navier-Stokes equations (RANS), in which a Spalart-Allmaras turbulent model is applied. The roughness on the icing surface is simulated by the equivalent roughness height method using Shin’s empirical formula [26], considering the effects of air velocity, , droplet diameter, and ambient temperature. The of the grid is less than 1 to ensure the grid independence and meet the boundary layer solution requirement of the turbulence model. The droplet flow field is one-way coupled with the airflow fields, and an Eulerian method [13] is applied to simulate the droplet impingement characteristics, obtaining the local droplet collection efficiency. The solutions of the air and droplet flow fields are conducted using the SIMPLE algorithm and second-order upwind scheme.
Based on the air and droplet flow field simulation, the results of the Messinger model can be obtained with UDFs. The solution of the unsteady icing thermodynamic model started from the stagnation point at when the icing surface is clean with no water film or ice layer. The unsteady solution method in Section 2.3 is then applied, and the ice shape is constructed.
4. Results and Discussions
4.1. Geometry and Conditions
A 2D NACA0012 airfoil with a chord of 0.5334 m is simulated to validate the unsteady icing thermodynamic model and analyze the development characteristics of the water film and the ice layer. The results are compared with the icing tunnel test data from NASA Glenn Research Center and the LEWICE simulations [27].
Aircraft icing is mainly affected by droplet diameter, , ambient temperature, flight velocity, angle of attack, and so on. Under different conditions, the ice accretion may form into glaze ice, rime ice, or mixed rime-glaze ice [28]. Since there is no water film formed for rime ice, the rime ice shape predicted by the unsteady model must be consistent with that of the traditional steady model. Therefore, only some glaze and mixed ice cases are selected from Ref. 27 for the icing simulations, as listed in Table 1. In the cases of R401 and R402, the ambient temperature and LWC are relatively low, leading to a mixed ice case; while the temperature and LWC in R421 and R422 are relatively high, resulting in a glaze ice state.
4.2. Effect of Time Step Selection
The selection of time step may affect the simulation results of the dynamic icing process. Three values of the time step (0.1 s, 0.01 s, and 0.001 s) are selected to conduct the unsteady calculation for the icing cases of R401 and R421, analyzing the effect of the time step selection on the model solution and validating the independence of time step. Figure 2 shows the variations of the maximum water film thickness and average ice thickness with time on the icing surface.
(a) Maximum water film thicknesses
(b) Average ice thickness
As shown in Figure 2(a), the maximum water film thickness increases with time, but its growth rate decreases gradually. After about 1.5 s of ice accretion, the film thickness becomes stable, and the result of R421 is larger than that of R401 due to the higher temperature and larger LWC. It can be found from Figure 2(b) that the average ice thickness increases roughly linearly with time, and the growth rate of R421 is slower than that of R401.The developments of the water film and the ice layer will be studied in detail in Section 4.4. Generally, the effect of the time step selection is appreciable on the water film development but not significant on the growth of the ice layer. The solution of the unsteady thermodynamic model can achieve the parameter independence of time step at about .
4.3. Unsteady Model Validation
Figure 3 shows the comparisons of the ice shapes with the literature data for R401 and R402. With a positive angle of attack, the stagnation point locates on the lower airfoil surface. When supercooled water droplets impinge on the icing surface, part of them freeze around the stagnation point due to high collection efficiency [5], and the rest would form a water film. The water film then flows backward and freezes on the downstream CV, ending at the runback water limit where all the liquid water turns into solid ice. In these two cases, the thickness of the ice layer increases along the airfoil surface from the stagnation point and then drops sharply around the runback water limits. After that, it decreases gently until reaches zero at the droplet impingement limits. The runback water limits locate within the droplet impingement range. Therefore, glaze ice is formed around the stagnation point with runback water, and rime ice is formed behind the runback limits without water film covered, indicating that the two cases are under a mixed rime-glaze ice state. Additionally, when the ambient temperature is lower in R402, the heat dissipation is larger, leading to a thicker ice layer around the stagnation point and a smaller runback region. Because the icing area in the mixed ice state is determined by the droplet impingement range and not affected by the ambient temperature, R401 and R402 obtain almost the same icing extent.
(a) R401
(b) R402
In addition, the ice thicknesses of the present model match well with the LEWICE results around the stagnation point, and both simulations are in acceptable agreement with the experimental data there. Difference is mainly in the ice shape near the runback water limits. At the upper limit, the ice thickness obtained by LEWICE is larger, and an ice horn is formed, which is in good agreement with the experimental data. The position of the upper runback limit is also well predicted by the present model, but the ice thickness is smaller with the smoother ice shape there. The reason is that the effect of the ice layer on the air and droplet flow fields is not considered in the present one-step icing simulations, which is critical for the prediction of ice horn. As for the lower runback water limit, the limit location and the ice shape of the present model show better agreement with the experimental data than those of LEWICE. In general, the present model could effectively predict the icing range, the runback water limits, and the ice shape for mixed rime-glaze ices.
The ice shapes of R421 and R422 are shown in Figure 4. Since the temperature and LWC are relatively high, the unfrozen droplets would also form a liquid water film at the stagnation point and then flow backwards. The water film covers the whole icing range, and the ice shape is smooth without any ice horn, which means these cases are under a glaze ice state. When the temperature is lower in R422, the runback water region and the icing range is smaller than those of R421, but the ice thickness is larger due to higher icing rate. The ice shapes of the present model show better agreement with the experimental data than those of LEWICE, especially at the lower airfoil surface.
(a) R421
(b) R422
To sum up, the present model could predict the ice shape acceptably, and the comparisons of the simulated ice shapes with the experimental data validate the unsteady model for the icing simulations under mixed and glaze ice conditions.
4.4. Developments of Water Film and Ice Layer
The dynamic developments of the water film and the ice layer are studied in detail with the present unsteady model. Figure 5(a) shows the water film thickness variation with time on the icing surface for R401, where is the curve length with at the stagnation point and at the lower surface of the airfoil. The thicknesses of the water film and the ice layer are both zero at the initial state. When water droplets impinge on the surface, the water film is formed firstly around the stagnation point since the impinging droplets cannot freeze totally due to the large collection efficiency there. Then, the film thickness increases and the runback water region expands over time. At about , the water film develops to a stable state, and the runback water limits are located within the droplet impingement range, also indicating that the ice accretion of R401 is under a mixed ice state.
(a) Water film thickness
(b) Ice accretion rate
(c) Surface temperature
The variation of the icing rate for R401 is shown in Figure 5(b), which is determined by both the impinged water flow rate and the heat dissipation capacity. Around the stagnation point, the heat dissipation capacity is less than the total freezing latent heat of the impinged droplets due to large collection efficiency, and there is water film formed with a constant icing rate. On the other hand, the heat dissipation capacity is much larger than the latent heat near the impingement limits due to small collection efficiency there, and all the impinged droplets freeze immediately without any water runback. Therefore, the icing rate is also constant near the impingement limits, and the value is the same with the impinged water flow rate. At the transition part between those two zones, as the runback water region expands, the icing rate and the water film thickness increase with time until the heat dissipation and latent heat balance at about , although the heat dissipation capacity changes slightly during this process.
Figure 5(c) shows the temperature variation during the water film and ice layer development. The glaze ice is formed with the water film around the stagnation point, and the temperature is right at freezing point, while the rime ice is formed downstream near the impingement limits, and the temperature there is lower than the freezing point. As the runback water flows downstream, the liquid film and glaze ice range expends to the position that was covered by the rime ice, and the surface temperature increases at the same time. This indicates that the present unsteady model could capture the icing type change in a CV and the effect of the temperature variation of the ice layer on the heat transfer characteristics. Additionally, the surface temperature also gets stable at about .
Figure 6 shows the changes of the water film thickness, the icing rate, and the temperature for R421. It can be seen that the water film thickness also increases, and the runback water region extends with time, leading to a change of the icing rate and surface temperature. Meanwhile, the glaze ice range expends, and the rime ice area shrinks as the time goes by. Compared with the results of R401, the water film thickness and glaze ice range are larger with smaller icing rate at the same time, since the ambient temperature and of R421 are larger than those of R401.The runback water region and the icing range on the upper surface are smaller than those of the lower surface due to the positive angle of attack. At , the runback water limit on the upper surface exceeds the droplet impingement limit, and the same situation occurs on the lower surface at . When the runback water and icing extent is larger than the impingement range, all the CVs on the icing surface change to the glaze ice state. In addition, the water film development also reaches stable at , and the water film thickness, runback region, surface temperature, and icing rate are constant after that.
(a) Water film thickness
(b) Ice accretion rate
(c) Surface temperature
In both cases, a relatively large value of the water film thickness is observed at the stagnation point. The reason is that the shear stress and pressure gradient are both very small there, leading to an accumulated liquid water and a thick film with low velocity. The air shear stress becomes larger along the airfoil surface, and then, the thickness of the water film decreases until it reaches zero at the runback water limits.
In general, water film would be formed in glaze ice and mixed ice cases. The thickness and runback region of the water film increase with time, and the icing rate and surface temperature vary accordingly. Rime ice occurs in the downstream area, and it may translate to glaze ice when the water film runs back over it. The water film development is rapid and would reach steady in a short period (3 s in present cases). Then, the surface temperature and icing rate are constant. The development of the water film and the variation of the surface temperature are also observed in the icing experiments in Ref. 5. The time for the surface temperature to reach a steady state is less than 10 s in the mixed and glaze icing cases, which are longer than those predicted in the present model. This might be due to the heat conduction effect of the substrate, and the time would change a lot with different substrate materials [6].
4.5. Comparison with Messinger Steady Model
The icing process can reach steady after the water film development. The predicted stable parameters such as the icing rate, the runback water flow rate of the film , and the surface temperature are compared with the results of the Messinger steady model mentioned in Section 2.1, as shown in Figure 7.
(a) Ice accretion rate
(b) Runback water flow rate
(c) Surface temperature
For the stable state in the unsteady model (), the results of the icing and temperature characteristics show good agreement with those of the Messinger model for both R401 and R421. Due to the distributions of the impinged water flow rate and the heat dissipation capacity, the icing rate increases from the stagnation point along the airfoil surface and reaches its maximum at the runback water limits. Meanwhile, the runback water flow rate increases at first due to the droplet impingement. As the collection efficiency decreases and the water film freezes, the runback water flow rate then decreases and gets to 0 at the runback limits. In the runback water region, liquid water and solid ice coexist, and the surface temperature keeps at the freezing point, while the temperature is lower in the rime ice range. Comparing with that of R421, the air velocity of R401 is larger and the ambient temperature is lower. Therefore, the runback water region and the mass flow rate of R401 are smaller with larger icing rate.
The growth processes of the ice layer predicted by the present unsteady model and the Messinger model are compared for R401 and R421 in Figure 8. The difference mainly locates at the runback water limits. At the early stage of the water film development in the unsteady model, the film is thin with a small runback region, leading to a small extent of the glaze ice. As the time goes by, the runback water region, as well as the glaze ice extent, gradually increases with the movement of the runback limits. On the other hand, the water film development is not considered in the Messinger steady model, and the icing rate and the runback water limits are constant from the very beginning. Therefore, the thickness of the ice layer predicted by the unsteady model is smaller than that of the Messinger model, especially at the runback water limits of the Messinger model since there is little water inflow until the end of the expansion of the water film in the unsteady model. Since the icing rate of the unsteady model increases over time to approach the value of the Messinger model, the deviation of the ice shapes would diminish gradually. At in these two cases, the relative error of the ice layer thickness is less than 5%.
(a) R401
(b) R421
From the above analyses and discussions, it is understandable that different icing thermodynamic models can predict similar ice shapes with the equivalent accuracy when comparing with experimental data, and both LEWICE with Messinger model and FENSAP-ICE with shallow-water model are accepted by airworthiness certification departments. However, at the initial stage of ice accretion or in the short-time icing simulations, it is necessary to consider the dynamic developments of the water film and the ice layer and use an unsteady thermodynamic model to predict the aircraft icing process. In addition, the dynamic developments of the water film and the ice layer might be longer considering the heat exchange between the water and the substrate, since the heat conduction in the aircraft skin might take a lot of time to stabilize according to its thermal properties.
As for the computational efficiency, the ice accretion characteristics could be obtained by solving the steady Messinger conservation equations only once from the stagnation CV to the end CVs on the icing surface, while the calculation cycles of the unsteady model are determined by the total icing time and the time step used in the simulation. It can be seen that with a small time step and a long icing time, the unsteady model costs much more computing resource than the steady Messinger one.
5. Conclusions
An unsteady icing thermodynamic model was established to consider the dynamic developments of the water film and the ice layer in the aircraft icing process. The corresponding implicit discretization and solution method were developed. Icing cases of a NACA0012 airfoil were numerically analyzed, and the results show acceptable agreement with the experimental data in the literature, which validates the present unsteady model. In addition, the developments of the water film and the ice layer were studied, and the main conclusions include the following: (1)The unsteady simulation would be affected by the selection of the time step, and the parameter independence of time step is achieved at (2)The unsteady model could effectively predict the development of the water film, which affects the surface temperature, the icing rate, and further the ice shape. The rime ice could change to the glaze ice due to the extension of the runback region(3)The water film develops rapidly, and the film thickness, the icing rate, and the surface temperature reach steady quickly. The stable values of the icing parameters in the unsteady model are consistent with those of the Messinger steady model(4)The dynamic developments of the water film and the ice layer play an important role in the aircraft icing simulations at the initial stage or in the short-time condition, and the effects are insignificant and even can be ignored in a long-time icing simulation(5)The temperature variation of the ice layer is considered in the unsteady model, but the heat transfer coupled with the skin substrate has not yet been taken into account, which will be studied in future research referring to the Myers model
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (no. 51806008), and the Open Fund of Rotor Aerodynamics Key Laboratory (no. RAL202104-2).