， ，

1. College of Energy and Power Engineering, Nanjing University of Aeronautics and Astronautics,Nanjing 210016, P.R. China；2. Key Laboratory of UAV’s Advanced Technology(Nanjing University of Aeronatics and Astronautics),Ministry of Industry and Information Technlogy, Nanjing 210016, P.R.China

Abstract: Under the condition of thermal anti-icing, the liquid water on the leading edge of the airfoil that would flow to the downstream non-protective zone will produce ridge ice, thus endangering flight safety. Based on the existing three-dimensional (3D) icing model which considers the water film flow on the ice layer, an icing model with thermal boundary condition is introduced. With the boundary conditions of none anti-icing and thermal anti-icing, glaze ice accretion and ridge ice accretion are simulated on a simplified airfoil of unmanned aerial vehicle(UAV), and then the lift coefficient and drag coefficient are calculated and compared with the smooth airfoil under the same conditions. The results show that the lift-drag ratio obviously decreases after glaze ice occurred on the leading edge under the condition of none anti-icing; and that after setting the condition of anti-icing heat flux in the impingement area, the glaze ice on the leading edge becomes thinner and the ridge ice occurs in the non-protective zone, so the airfoil with this icing characteristic gets a lower lift-drag ratio.

Key words: glaze ice; ridge ice; water film flow；thermal anti-icing; lift coefficient

0 Introduction

When airplanes pass through clouds which contain supercooled water droplets, ice will accrete on the surfaces of the upwind components, such as the airfoils and the inlet blades of the engine[1]. Ice reduces the lift as much as 35%[2], and the drag can be increased by three times compared with a smooth airfoil. Ice layer also destroys the aerodynamic characteristics of components and causes the stalling of aircraft and the surge of engine. Obviously, ice accretion affects the flight performance and flight safety of the aircraft, and even causes the plane to crash[3].

Ice accretion caused by supercooled water droplets first occurs on the leading edge of the component[4]. Due to the different environmental conditions of ice accretion, it mainly consists of rime ice and glaze ice[5]. Rime ice is loose and the ice shape is closer to the aerodynamic shape. Glaze ice is firm and behaves double horn shape, so the impact on the flight performance is more serious. In order to reduce the danger of ice accretion, heat flux is generally provided to the front ice accretion zone. On this occasion, it is difficult for all the supercooled water to freeze in this zone. The water film would flow downstream and freeze in the zone without heat flux[6]. This is the ridge ice and it is very difficult to be found. It has the same characteristic as the glaze ice, which has firm structure and poor aerodynamic shapes. Research shows that ridge ice has serious influence on the aerodynamic performance of airfoil[7]. So the study on ridge ice is very important to flight performance and safety.

Numerical simulation is an important way to study the ice accretion of aircraft. Several icing software has been developed abroad, such as LEWICE[8](USA), CIRA[9](Italy), TRJICE[10](UK) and ONERA[11](France). All these software is developed on the Messinger icing model[12], which assumes that all the unfrozen water flows into the next control volume. The assumption is impractical on the airfoil surface, and this leads to some simulation deviation, especially for the glaze ice. Moreover, the “next” control volume is not unique on the 3D surface. Thus the application of the Messinger icing model is limited in the 3D icing simulation. In another icing software FENSAP-ICE[13](Canada), the water film flow is considered and the velocity is proportional to the shear force. In domestic research, Chen[14]simplified the N-S equations by dimensional analysis. So the one-dimensional governing equations for water film flow were obtained and the two-dimensional icing model was developed. Good results on the airfoil were obtained with this icing model. However, the model only considers the flow of the water film in the chord direction，while the flow in the span direction is ignored. Cao et al.[15]developed a 3D icing model in which the two-dimensional governing equations of water film flow were obtained. Some typical 3D ice accretions were simulated preliminarily and certificated.

After the crash of ATR-72, the research on ridge ice was mainly focused on the formation of ridge ice and its influence on the aerodynamic characteristics[16]. Tan et al.[17]obtained the location and height of ridge ice under certain supercooled large droplets. Bragg[18]found that ridge ice decreased the lift and increased the drag of the airfoil by experiments. In the domestic research, Yi[19]improved the thermodynamic model and proposed a simple method to simulate the 3D ridge ice. Xiao et al.[20]conducted experiment on the formation process of ridge ice under electric heat flux. Wang et al.[21]simulated the ridge ice with the Messinger icing model. Until now, the research on the ridge ice was mainly carried out in icing tunnel experiments or simple two-dimensional simulation.

This paper will develop the 3D icing model established in Ref.[15]. The heat flux for anti-icing is introduced into the model. The glaze ice and ridge ice on the 3D low-speed airfoil are simulated with the new icing model. Then, the lift and drag are computed and compared to the smooth airfoil. The study in this paper will benefit the ridge ice research in the future.

1 Mathematical Model and Calculation Method for Ridge Ice Accretion

1.1 Assumptions

According to the characteristics and complexity of 3D ice accretion, the assumptions are proposed as follows:

(1) The incoming flow condition is assumed to be constant.

(2) The thermophysical properties of air and supercooled water are assumed to be constant. Once the water film is freezing, the thermophysical properties change to be the same as ice instantly.

(3) The surface tension of the water film and the influence of the water droplets impact are ignored. The water film is considered to be laminar flow.

(4) The evaporation of the water film and the sublimation of the ice layer are ignored. The convection term in the film heat transfer equation is also ignored.

(5) The roughness of the airfoil surface is assumed to be constant.

(6) The influence of water droplets on the lift coefficient is ignored.

1.2 3D ice accretion mathematical model

The control volume of ice accretion is shown in Fig.1. There are three layers from top to bottom: air-supercooled water droplets flow, thin water film flow and ice layer.

Fig.1 Control volumn

The flow of the thin water film directly affects the ice accretion. In Fig.1:mimpis the mass flux of supercooled water droplets impingement;Qis the given heat flux on the airfoil substrate for anti-icing;HwandTware the thickness and temperature of the water film, respectively;uandvare the flow speeds of the water film inxandydirections, respectively;HiandTiare the thickness and temperature of ice layer. Based on the icing model in the Ref.[15], which coupled ice accretion and water film flow, the mathematical model that simulates ridge ice accretion will be developed. The model contains seven equations as follows: mass equation, momentum equations, energy equations, and the Stefan equation.

(1)

(2)

(3)

(4)

(5)

(6)

(7)

Eqs.(1)—(4) are the flow equations of unfrozen thin water film. Eqs.(5) and (6) are the energy equations of water film layer and ice layer. Eq.(7) is the growth rate equation for ice layer. In these equations,ρiandρware the density of ice and water, respectively;λiandλware the thermal conductivity of ice and water;υ,gandLfare the dynamic viscosity, acceleration of gravity and latent heat of solidification, respectively.

The boundary conditions for Eqs.(2)—(4) are the same with that in Ref.[15]. And the boundary conditions for Eq.(5) are[15]

Tw|z=0=Tf

(8)

(9)

whereTfis the water freezing temperature，andCpwis the constant-pressure specific heat capacity in these equations.

Without heat flux in the leading edge area, glaze ice accretion occurs. The boundary conditions for Eq.(6) are

Ti|z=0=Tf

(10)

(11)

When heat flux is introduced in the leading edge area, the glaze ice layer will be thinner. And the unfrozen water film will flow to downstream and freeze to the ridge ice, Eq.(11) in the heating region is modified as follow

(12)

According to the results in Ref.[22], the anti-icing thermal load is the lowest when the heating region is the same with the water droplets impingement area. Therefore, this paper sets the impingement area as the heating region. A given heat flux is introduced in this region to decrease the thickness of glaze ice layer.

From the mathematical model of ice accretion,Hw，Hi，TwandTiare coupled with each other. However, whenQis given, it is possible to obtainHifirst from Eq.(7). And then all the equations are solved orderly[15].

1.3 Calculation method of lift and drag in 3D airfoil

In this paper, the lift coefficientCland drag coefficientCdin ANSYS-FLUENT are used to describe the lift and drag on the 3D airfoil.

(13)

(14)

whereLandDare the lift and drag, respectively, which are calculated by ANSYS-FLUENT;ρandV∞are the flow density and speed, respectively;Sis the reference area, which means the projected area of the airfoil on the horizontal plane.

1.4 Wall roughness setting

In the ice accretion simulation, the roughness of the airfoil or the ice layer surface has a great influence on the convective heat transfer coefficient, which will affect the ice accretion very much. In the lift and drag calculation, the roughness also has a great influence on the properties of the boundary layer, which will affect the lift and drag. Therefore, the roughness setting of the surface is one of the most important concern in this paper.

In fact, when the incoming icing conditions change, the roughness of the ice layer surface varies. Research is still underway on the mechanisms and characteristics of roughness generation[23]. Currently, most of the studies on ice accretion use the equivalent roughness height to describe the roughness of ice layer. The empirical equations[8]introduced by the NASA Lewice Center are the most widely used. In these empirical equations, the equivalent roughness height is related to the size of the airfoil and some incoming icing conditions.

(15)

0.571 4+0.245 7LWC+1.257 1LWC2(16)

(17)

(19)

whereksis the equivalent roughness height,cis the chord length of airfoil, parameters with subscript “base” represent the reference values, and parameters with subscripts LWC，T∞and MVD represent the effects of liquid water content, temperature, and maximum droplet diameter, respectively.

1.5 Calculation software and its settings

In this paper, ANSYS-CFX is used to calculate the air-water droplets flow around the airfoil. The continuous-discrete fluid model based on Euler-Euler method is used. The air is set as continuous fluid phase, while the water droplets is set as discrete fluid phase. The sum of the volume fractions of air and water droplets is 1. Flow turbulence is treated byk-εturbulence model.

The lift and drag of airfoil are calculated by ANSYS-FLUENT. Since the velocity of incoming flow is low, pressure-based solver is used and the flow turbulence is also treated byk-εturbulence model.

2 Simulation and Analysis

2.1 Physical model

The 3D airfoil of an unmanned aerial vehicle (UAV) is simplified as shown in Fig.2. Letcrepresent the average chord length,Lthe span length, AOA the angle of attack, andχthe sweepback angle, respectively.

Fig.2 Swept airfoil

The calculated domain is shown in Fig.3. The fuselage is set as a large flat wall and the winglet is ignored. The upstream area is set as a semicircular region with radius of 6c, while the downstream area is set as 8c×12crectangular region. The surface “Wall” (airfoil and fuselage) is set as non-sliding. The surface “In” (surfaces upstream and beside the airfoil, the surface on the winglet) is set as velocity-inlet. The surface “Out” (surface downstream the airfoil) is set as pressure outlet. The parameters for the airfoil are shown in Table 1.

Fig.3 Computational domain

c / mL / mAOA / (°)χ / (°)0.217 51.523

The simulation for 3D ice accretion can be mainly divided into five parts: (1) physical modeling and meshing, (2) air/supercooled flow simulation, (3) water droplet impingement calculation, (4) ice accretion and water film flow simulation, (5) ice layer shape and mesh reconstruction. The total ice accretion time is 600 s. During the process of simulation, the meshing and the calculation for two-phase flow are updated every 60 s. The lift and drag can be obtained after the ice layer shape is updated.

2.2 Validation example

The glaze ice on the NACA0012 airfoil is simulated to validate the rationality of the 3D icing model developed in this paper. The physical model and computational parameters (as shown in Table 2) are the same as the example in Lewice[8].

Table 2 Computational parameters

The results are shown in Fig.4. It shows that both the area and the volume of the ice accretion simulated in this paper are close to the results in Ref.[8]. This means that the 3D icing model is reliable.

Fig.4 Comparison results

2.3 Calculation model and meshing

Hexahedral structured mesh is generated by ANSYS-ICEM. The air/supercooled flow and the lift/drag are both simulated on this mesh. The water film flow and ice accretion are both calculated in the first layer mesh. Thus the mesh near the wall should be dense and well-orthogonal.

The mesh used in this paper is shown in Fig.5. The height of the first layer mesh is about 0.1 mm. The growth factor is set as 1.2. The first layer mesh shows good orthogonality to the airfoil surface.

Fig.5 Mesh on the airfoil surface

2.4 3D simulation of ice accretion

The incoming icing conditions are shown in Table 3.

Table 3 Icing conditions

Without heat flux (Q= 0) in the leading edge area, the 3D state of glaze ice accretion on the airfoil after 600 s is shown in Fig.6.

Fig.6 3D glaze ice after 600 s

When heat flux (Q=2 500 W/m2) is introduced in the leading edge area, the 3D state of ridge ice accretion on the airfoil after 600 s is shown in Fig.7.

Fig.7 3D ridge ice after 600 s

Since the boundary effect near the wall is significant, the results in the cross section of the intermediate airfoil are selected for comparison as shown in Fig.8 and Fig.9. Fig.8 shows that the glaze ice occurs on the leading edge without anti-icing. Fig.9 shows that when heat flux is introduced to the impingement area, the glaze ice becomes thinner and the ridge ice occurs in the downstream area without anti-icing. In Fig.9,sis the arc length. The reason is that: Without the heat flux in the leading edge area, the icing ability is so strong that all water film are frozen; when the heat flux is introduced to the leading edge area, the icing ability in this area is weakened, so that the unfrozen water film flows to the downstream and ridge ice occurs.

Fig.8 2D glaze/ridge ice after 600 s

Fig.9 Water film thickness distribution

2.5 Calculation of lift and drag

The lift coefficient and drag coefficient are calculated after 600 s and the results are compared to the smooth airfoil (0 s) as shown in Fig.10.

Figs.10(a) and (b) show that the aerodynamic performance of icing airfoil is worse over time. From 2 min to 10 min, the lift coefficient of the airfoil with glaze ice is reduced from 97% to 57% compared to the smooth airfoil while the drag coefficient increases from 113% to 319%. The lift coefficient of the airfoil with ridge ice is reduced from 95% to 41% compared with the smooth airfoil while drag coefficient increases from 120% to 449%. So the airfoil with ridge ice gets a less lift but a larger drag.

Fig.10 Aerodynamic performance for smooth and iced airfoil

Fig.10(c) shows that the lift-drag ratio of the airfoil with ice is lower than the smooth airfoil. And it is continuously reduced over time. From 2 min to 10 min, the lift-drag ratio of the airfoil with glaze ice is reduced from 86% to 18% compared to the smooth airfoil, while the lift-drag ratio of the airfoil with ridge ice is reduced from 79% to 9%. It means that if the anti-icing heat flux is arranged improperly, the glaze ice and ridge ice would both occur and the aerodynamic performance will get worse. Therefore, the area and the heat flux should be optimized in the anti-icing design of aircraft, which will be studied in the future.

2.6 Flow field analysis

The flow field near the airfoil with ice is shown in Fig.11. It shows that when the air passes through the airfoil with glaze ice, a recirculation zone would appear behind the ice horn, which reduces the lift and increases the drag. Moreover, when the air passes through the airfoil with ridge ice, the recirculation zone is larger. This is why the aerodynamic performance get worse.

Fig.11 Effect of ice shape on flow field

3 Conclusions

Based on the 3D icing model established in Ref.[15], the ridge icing model and its calculation methodology are developed under the thermal anti-icing condition. The glaze ice accretion and ridge ice accretion are simulated. The lift coefficient and drag coefficient are compared to the smooth airfoil. Some results are concluded as follows:

(1) Based on the 3D icing model and its calculation methodology developed in Ref.[15], ridge ice accretion can be simulated by changing the thermal boundary conditions from adiabatic to heat flux.

(2) When the heat flux for anti-icing is introduced, some water film would flow from the leading edge to the downstream area, so glaze ice and ridge ice are both occurred.

(3) Compared with the smooth airfoil, both glaze ice and ridge ice would decrease the lift and increase the drag.And this change would be worse over icing time. After 10 min of ice accretion, the lift-drag ratio of the airfoil with glaze ice drops to 5.23, while the lift-drag ratio of the airfoil with ridge ice drops to 2.66. This means that the protective area and the heat flux should be reasonably arranged in the anti-icing design.

Acknowledgements

This work was financially supported by the Natural Science Foundation of Jiangsu Province (No.BK20150740)，and the National Natural Science Foundation of China (No.51506084).