Abstract
A method to simulate ice accretion on an aircraft wing using a three-dimensional compressible Navier-Stokes solver, a Eulerian droplet flow field model, a mesh morphing model, and a thermodynamic model, is presented in this paper. The above models are combined together into one solver and implemented in OpenFOAM. Two-way coupling is achieved between airflow field calculation and ice simulation. The density-based solver rhoEnergyFoam is used to calculate the airflow field. The roughness wall function is proposed to simulate the roughness effect caused by ice accretion. For droplet flow field calculation, the Eulerian model is applied and the permeable wall boundary condition is used on the wing to simulate the droplet impingement. The icing thermodynamic model is built based on the Messinger model. The mesh morphing model adjusts the wing’s shape every time step based on the amount of accreted ice so that the airflow field is updated during the simulation. The effect of the ice accretion on the airflow is studied by comparing the aerodynamic performance—with and without ice. The ice accretion on the ONERA M6 wing model under a specific condition has been simulated to validate the solver’s performance and investigate the effect of the accreted ice on the aerodynamic performance.
1. Introduction
Aircraft can experience icing when encountering a cloud that contains supercooled water droplets. It is one of the main hazards to flying due to the degradation of aerodynamic performance because the ice accretion on the wing alters the originally designed aerodynamic configuration. Investigation on aircraft icing mechanism and effects is the basis of the anti/deicing technique and the establishment of the flight and operation rules in icing condition.
Considerable research has been done, and several approaches have been used to investigate the ice accretion, including the flight test, experimental study, and numerical simulation. Since it is expensive to use the flight test and experimental simulations, numerical simulation is adopted widely. Due to the significant advances of high performance computing, numerical methods have obtained considerable development. Several codes have been developed to simulate the ice accretion process, and the representative ones are as follows: LEWICE code of NASA [1], TRAJICE code of DRA [2], FENSAP-ICE code [3, 4], the code of ONERA [5], and the code of CIRA [6, 7].
It has been demonstrated that LES modeling produces a more accurate prediction for aerodynamic performance, especially in the study of the near field of the aircraft wake compared to RANS modeling [8, 9]. However, when applied to turbulent boundary layers, the use of LES is much more computational expensive because of the wall treatment. It is well known that the number of grid points required by LES in the boundary layer exponentially increases with the Reynolds number, which makes LES not feasible on realistic wall-bounded flows. The alternative is to use wall-modeled LES or hybrid RANS-LES. The former imposes a mean velocity profile in the boundary layer whereas the latter couples the RANS equations of calculating the mean flow close to the wall while the LES resolves the turbulence and vortex dynamics far from the wall. In this work, we apply the second method. The Menter’s K-Omega SST model [10] is adopted in the current work, and this paper would be an effort to the development of the advanced icing solver based on hybrid RANS-LES modeling. In this paper, the numerical simulation consists of four subsolvers: airflow solver, droplet motion solver, ice thermodynamic solver, and mesh motion solver.
The airflow field is obtained by using a rhoEnergyFoam solver [11] to solve the compressible Navier-Stokes equations. The rhoEnergyFoam solver is a newly released solver in OpenFOAM open-source community [12], which has shock-capturing capabilities, and it does not require the burden associated to the calculation of the eigenstates of Navier-Stokes equations as in characteristics-based schemes, although it maintains the same level of accuracy of these schemes. The roughness wall function is proposed to simulate the roughness effect caused by ice accretion.
The Eulerian two-phase model is adopted in the droplet motion solver coupling with the airflow solver. The droplet-phase governing equations are established by considering the droplets in the airflow as a pseudo-fluid and expressed as a continuity equation and a momentum equation. The droplet collection efficiency is obtained by solving the governing equations. In addition, the wall boundary conditions have to be modified to represent the droplet impingement properly. In the current work, the permeable wall boundary condition is proposed to mimic the impingement of droplets on the wing.
The thermodynamic model in the ice thermodynamic solver is based on the classical Messinger model [13] employed by most of the ice accretion solvers. By solving the mass balance and energy balance equations, the ice amount generated in every time step is obtained. Then, based on the assumption that the ice only grows in the normal direction to the wing’s surface, ice shape is calculated. Therefore, the wing’s shape is changing during the simulation because of the ice accretion and airflow field is updating together with it. The droplet collection efficiency is also recalculated, and the new ice shape can be built by repeating the process.
This paper first presents the four subsolvers used for the simulation and describes the basic equations. Second, the validation of the solver is performed on a two-dimensional airfoil icing case. Third, the ice accretion on the three-dimensional ONERA M6 wing is simulated to study the solver’s performance and the results are discussed.
2. Governing Equations of the Two-Phase Flow
2.1. Airflow Solver
For the simulation of unsteady compressible flows, the rhoEnergyFoam [11] solver is used. rhoEnergyFoam is a newly released solver in the OpenFOAM open-source community. It is based on AUSM flux splitting, which has shock-capturing properties, and it does not require the burden associated to the calculation of the eigenstates of Navier-Stokes equations as in characteristics-based schemes, although it maintains the same level of accuracy of these schemes. In addition, it discretely preserves kinetic energy; indeed, it has been shown to be less dissipative in canonical compressible flow simulations (such as isotropic turbulence and forward-facing step) and also in aerodynamics flow simulations in complex geometries as those targeted in this project.
The governing equations of the rhoEnergyFoam solver are the Navier-Stokes equations for a compressible ideal gas, integrated over an arbitrary control volume . where is the outward normal and , , and represent the conservative variable vector, the associated Eulerian, and viscous fluxes, respectively. Here, is the density of the air, is the velocity component in the -th coordinate direction, is the total energy per unit mass, is the internal energy per unit mass, is the total enthalpy, is the viscous stress tensor, and is the heat flux vector.
Then, the Eulerian fluxes are split into the convective and pressure contribution:
The convective and pressure fluxes both consist of the central part and the diffusive part. The central part of the convective flux is evaluated through the midpoint interpolation scheme which discretely preserves kinetic energy of the flow from convection, and the pressure flux is evaluated through the standard central interpolation. In order to maintain the stability of the simulation, some amount of numerical diffusion is necessary. A shock sensor is used in the rhoEnergyFoam solver to judge the smoothness of the numerical solution. When capturing shock waves, the artificial diffusion terms provided by the AUSM [14] scheme are applied to pressure and convective fluxes.
The RANS (Reynolds-Averaged Navier-Stokes) models usually underpredict the convective heat flux because they consider the wall as smooth while the icing surfaces are quite rough. Therefore, the thermal wall function derived by Lima da Silva et al. [15] has been implemented to study the roughness effect, in which the roughness effect is considered by using the momentum and heat transfer analogy factor. Here, is the turbulent thermal diffusivity. where is the turbulent viscosity, is the turbulent Prandtl number, is the rough skin friction, is the shear velocity, and is the equivalent sand-grain roughness. , , and are defined by Stefanini et al. [16] as three constants as -0.45, -0.8, and 1.42, respectively.
2.2. Droplet Motion Solver
The governing equations for the droplet phase are established based on the following assumptions: (1)The distribution of the water droplets is uniform, and they are simplified as sphere with a median volumetric diameter(2)The physical parameters of the droplets does not change by assuming that there is no heat or mass transfer between the droplets and air(3)The droplet collision, splashing, and bouncing effects are neglected(4)The airflow viscosity has no effect on the droplets where represents the droplet volume fraction and and are the droplet density and velocity, respectively. In the momentum equation, is the gravity factor and is the drag force caused by airflow, which can be calculated as follows: where is the Reynolds number, , and is the droplet diameter. can be calculated by
The Gauss’s theorem is firstly applied to convert the spatial integral to surface integral. Then, the first-order upwind scheme is used for divergence terms, and the time integration of the resulting ordinary differential equations system is carried out by a third-order, four-stage Runge-Kutta algorithm. where represents the conservative variables and , is the time step size, and is the spatial derivative.
Based on the previous assumption that the droplet splashing and bouncing effects are neglected, the droplets would adhere to the wall after impingement, which means that the high velocity of droplets comes to zero instantaneously upon impingement. However, since the viscosity is not considered in the droplet phase, the no-slip wall boundary condition applied in the airflow field cannot be used in the droplet phase. On the other hand, the droplet volume fraction cannot be set as zero on the wall because it will prevent the occurrence of any impingement. Therefore, a specific wall boundary condition is applied on the wing to satisfy the mass conservation while at the same time it allows for droplet impingement. For simplicity, it is illustrated in Figure 1 in which the droplet velocity component normal to the wall is represented by . Upon the impingement, the zero-gradient boundary condition is applied for both droplet velocity and droplet volume fraction. Otherwise, if the droplet splashing is detected, the fixed value boundary condition is applied so the droplet volume fraction is assigned with a minimal value and the droplet velocity is set to be zero to prevent it from splashing on the wall. In order to differentiate between the impingement area and the non-impingement area, in the control volumes adjacent to the wall, we check the direction of the droplet velocity component in the normal direction of the wall. This is illustrated in Figure 2: when , no impingement occurs while when , droplet impingement on the wall is expected.
(a)
(b)
The droplet volume fraction and droplet velocity are obtained by solving the governing equations. Then, the nondimensional local collection efficiency is calculated by where is , represents the normalized droplet volume fraction, represents the normal direction, and is the droplet velocity in the free stream.
3. Thermodynamic Modeling
Different icing conditions cause different types of ice accretion. Rime ice accretion occurs when supercooled water droplets freeze instantaneously upon impact on the aircraft wings. If the environmental temperature is not low enough, the droplets may runback as liquid water to form glaze ice. The current thermodynamic model is based on the Messinger model [13] which was employed by a few commercial icing codes such as LEWICE [1], TRAJICE [2], and ONERA [5]. However, these icing codes assume that the freezing fraction stays constant during the icing simulation, which is not accurate. More recently, Cao et al. [17, 18] include the freezing fraction as a changing variable in the icing simulation. In this paper, the freezing fraction is updating throughout the simulation and the ice height is obtained by solving the mass balance equation and energy balance equation.
3.1. Mass Balance
As shown in Figure 3, in the control volume on the wing’s surface, the mass balance equation can be written as follows: where represents the mass flux of water impacting the wing’s surface, is the mass flux of generated ice, represents the mass flux of evaporation or sublimation, and and denote the mass flux of water flow into the control volume and flow out of the control volume, respectively. where is the liquid water content and is the freezing fraction.
3.2. Energy Balance
Figure 4 shows the energy balance equation which is necessary to obtain the balance temperature and the freezing fraction. Its establishment and solution are based on the following assumptions: (1)There is no runback water in the control volume at the stagnation point, and any runback water flowing out of the control volume flows along the direction away from the stagnation point(2)The heat and mass transfer only happens in the direction normal to the wing’s surface(3)In the mixture of water and ice, a balance temperature is reached
Then, the energy balance equation can be written as follows: where the four energy terms represent the convective heat, impingement heat, latent heat, and sensible heat, respectively. The descriptions of these energy terms are given below.
The convective heat can be calculated as where and represent the balance temperature and recovery temperature in the control volume and is the convective heat transfer coefficient.
The impingement heat can be calculated as where is the water impingement mass flux and is the water droplet impinging velocity. Both of them are obtained from the droplet solver solutions.
The latent heat can be calculated as the sum of latent heat of freezing and latent heat of sublimation or evaporation. where and represent the latent heat of freezing and the latent heat of sublimation or evaporation, respectively.
The sensible heat consists of all the heat exchanged by water without changing state. Firstly, the temperature of the impinging water is increased to freezing temperature at which the state change happens. Then, in the control volume, the newly formed ice and the liquid water reach the balance temperature . where and represent the specific heat of water and ice, respectively. and represent the temperature of the impinging water droplets and the temperature of runback water from the previous control volume, respectively.
3.3. Solution of the Mass and Energy Balance Equations
By combining equation (15) with equation (13), an equation with unknown variables and can be obtained. The value of must be between 0 and 1 due to the icing physical process, which provides a constraint for this equation. The solving procedure is a predictor-corrector methods which starts from the stagnation point and along the upper and lower surfaces of the wing. Firstly, it is assumed that both the water and ice exist in the control volume and the balance temperature is reached. Then, the freezing fraction can be calculated, which have three possible scenarios: (1): the assumption is correct. Both water and ice exist in the control volume and the calculation will be continued to the next control volume based on this (2): there is no water in the control volume and the value of will be assigned with 1 to continue the calculation(3): there is no ice accretion
3.4. Ice Shape Reconstruction
Based on the assumption that the ice grows only in the normal direction to the wing’s surface, the ice height during time step can be written as where is the density of the ice.
The irregular ice shape represents a major challenge in numerical simulation of long-time icing. Because the ice accretion changes the wing’s shape, the previously solved airflow field and droplet flow field are also influenced. Therefore, it is necessary to recalculate the airflow field and droplet flow field in order to obtain the new ice shape accurately based on the updated mesh. However, manual remeshing is a time-consuming procedure. In the current work, a mesh movement scheme is presented.
At each time step, the mesh nodes next to the wing move in the normal direction to the wing’s surface with the distance calculated from equation (20). Therefore, the ice shape is constructed by the first layer nodes. The movement of the rest of the mesh nodes is decided by solving the Laplace equation with variable diffusivity. where is the diffusivity and is the node motion displacement. The diffusivity field is calculated based quadratically on the inverse of the cell center distance to the nearest boundary which is a wing in this case. Therefore, the nodes which are far from the wing have smaller motion displacement than the ones which are close to the wing.
4. Results and Discussions
In this section, the developed solver is tested against various two- and three-dimensional cases. Firstly, icing simulation is performed on two-dimensional NACA0012 airfoil for validation. Then, the ice accretion on three-dimensional ONERA M6 wing and the influence of the ice accretion on the airflow field are investigated.
4.1. Ice Accretion on NACA0012 Airfoil
Figure 5 shows the computational mesh for NACA0012 airfoil, which contains 6200 cells. The chord of the airfoil is , and the computational far field domain extends to . The angle of attack is 4°, the free stream velocity is , and the pressure is . The static temperature . For the droplet phase, and the median volumetric diameter (MVD) is 20 μm. The icing time is 360 seconds.
Figure 6 shows the ice shape comparison between the predicted results in this paper, experimental data [19], and predicted results from Cao et al. [17]. The typical rime ice accretion happens in this case due to the relatively low environmental temperature. From the comparison plot, a good agreement is obtained.
4.2. Ice Accretion on ONERA M6 Wing
The mesh of ONERA swept M6 is shown in Figure 7 and contains 671160 cells. The mean aerodynamic chord is , the semispan is , and the computational domain extends to . The angle of attack is 6°, the free stream Mach number is , the static temperature , and the Reynolds number is [6]. For the droplet phase, and . The icing time is 200 seconds.
The calculated collection efficiency is presented in Figure 8 which shows that its value is higher in the lower region of the wing. Since the wing is at 6° angles of attack, the main impingement region is the lower surface of the wing. The same phenomenon is also observed in the ice height distribution presented in Figure 9. The predicted ice shape at three sections: section A 90% span, section B 60% span, and section C 20% span, is shown in Figures 10(a)–10(c). It can be seen that the thickness of the ice accretion increases from the root to the tip of the wing which agrees with the distribution of the local collection efficiency.
(a) Section A
(b) Section B
(c) Section C
In order to study the impact of ice accretion on the aerodynamic performance of the wing, we compare the computed distributions of pressure coefficient with and without ice, which are shown in Figures 11(a)–11(c). and denote the coordinates of the leading edge and trailing edge of each wing section, respectively. In the three sections, the ice accretion’s influence on pressure coefficient at section A is the most significant one because the ice thickness increases from the root to the tip. And the irregular ice shape severely changed the pressure coefficient distribution.
(a) Section A
(b) Section B
(c) Section C
Also, in order to study and understand the impact of ice accretion on the tip vortex, Figure 12 shows the comparison of peak vorticity magnitude along the free stream direction with and without ice. Ice accretion on the wing reduced the strength of the vortex close to the wing tip. It also caused some disturbance to the peak vorticity along the free stream direction due to the irregular shape of the ice accreted on the wing. The peak vorticity keep decreasing as it moves away from the wing tip when ice formation is not considered. However, icing produces local overshoots as the nonuniform growth on the wing surface. These oscillations eventually smooth out at sufficiently large distance from the trailing edge. The same phenomenon is also observed when comparing vorticity magnitude distribution with and without ice at the three sections after the wing tip: , 0.22, and 0.5 as shown in Figure 13. It can be seen that the vorticity distribution for the iced wing is much more irregular and the vortex with medium strength is more widely spread. This could have potential impact on the safety space between aircraft at takeoff and landing [20] because it has been shown in this study that the ice accretion on the wing will affect the aircraft wakes.
(a)
(b)
5. Conclusions
A three-dimensional model is established to predict the ice accretion on the wing as well as the ice’s impact on the aerodynamic performance of the wing. The ice accretion has been simulated on a 2-D NACA0012 airfoil and a 3-D ONERA M6 wing. All the results show that this solver is feasible.
The model solves the icing problem by a two-way coupled procedure. The airflow field is solved by using density-based solver rhoEnergyFoam with the AUSM flux splitting scheme. The droplet motion field is solved explicitly in an Eulerian field with Runge-Kutta time-advancing scheme. The specific impinging boundary condition and the roughness wall function effectively help with simulating the icing phenomenon. The remeshing model based on a Laplace equation functions well to capture the ice shape and prepare for the study of the updated airflow field.
The study on the impact of ice accretion on the aerodynamic performance of the wing shows that the pressure coefficient near the tip of the wing is greatly affected by the ice. Also, it is noted that the ice accretion has important effects on both the tip vorticity distribution and the peak vorticity along the free stream direction. This work represents a first step of a larger project that is aimed at developing a hybrid RANS-LES solver to simulate icing over aerodynamic surfaces. Next steps will be, first, to test detached-eddy simulations based on (1) or (2) equation turbulence models and then to couple the ice formation/growth model into those solvers.
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 they have no conflicts of interest.
Acknowledgments
Computer time was provided by ACCC at UIC on the university cluster “Extreme.” This work was supported by Argonne National Laboratory through grant #ANL 4J-30361-0030A “Multiscale Modeling of Complex Flows.”