Sage Journals: Discover world-class research

Abstract

It has been widely investigated before that both heavy rain and icing conditions can cause severe aerodynamic performance penalties to aircraft in flight; however, aircraft aerodynamics in the presence of both heavy rain and ice accretion has been little explored so far. In this article, a numerical simulation is performed to study the aerodynamic performance of an iced NACA 0012 airfoil in simulated heavy rain condition. A two-way momentum-coupled Eulerian–Lagrangian two-phase flow approach is developed to simulate raindrop trajectories in the airflow field. The analysis indicates that both rain and ice can separately cause severe aerodynamic performance degradation to the airfoil at most angles of attack of interest. In the presence of both ice and rain, a further severe aerodynamic penalty can be produced by rain based on the ice-induced aerodynamic performance degradation. Premature boundary-layer separation and boundary-layer momentum losses by splashed-back raindrops are also observed.

Keywords

Numerical simulation heavy rain ice accretion premature boundary-layer separation boundary-layer momentum loss

Introduction

Aircraft performance degradation due to flying in adverse weather conditions remains a constant concern within the aviation community. Herein, two factors, that is, heavy rain and icing, have attracted particular attention due to their undetermined detrimental effects on aircraft aerodynamic performance.^1–3 Although there have been many research achievements on the effects of icing^4–6 and heavy rain^7–9 on aircraft aerodynamic characteristics, none is found to consider the coupling effects of the two. However, a situation that an aircraft with an initial leading-edge ice accretion is flying through a cumulonimbus cloud or rain condition actually exists, especially in the Northeastern America where the weather is always cold and sometimes with rainfall in winter. Additionally, aircraft icing also occurs in freezing drizzle or rainfall where the coupling effects of rain and icing exist. Thus, it is necessary and significant to study the coupling effects of rainfall and ice accretion on aircraft aerodynamics.

In this article, a two-way-momentum-coupled Eulerian–Lagrangian two-phase flow approach is adopted to study the aerodynamic characteristics of a NACA 0012 airfoil with post-leading-edge ice accretion in a simulated heavy rain encounter on the basis of our earlier work.⁹ Regarding the improvement of the previous method, the dispersion of particles due to turbulence in the continuous phase is taken into consideration in this study by incorporating a random walk dispersion model. Both the individual and coupling effects of heavy rain and ice accretion on the aerodynamic performance of the airfoil are investigated. Considering the fact that there is no open report on the current topic, our study can provide a new insight into the effects of the two common adverse weather factors on aircraft aerodynamics.

Methodology description

Continuous phase

The steady-state incompressible Reynolds-averaged Navier–Stokes equations of mass, momentum, and energy are solved for the airflow field. The pressure-based segregated SIMPLE algorithm is adopted for pressure–velocity coupling calculation. The convection terms were discretized using the second-order central scheme, while the viscous terms using the QUICK scheme.¹⁰ The viscous coefficients were solved by associating with the S–A turbulence model. The governing equations for the airflow field are written as follows

\frac{\partial {\bar{u}}_{i}}{\partial x_{i}} = 0

(1)

ρ_{a} \frac{\partial {\bar{u}}_{j}}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [- \bar{p} + μ (\frac{\partial {\bar{u}}_{i}}{\partial x_{j}} + \frac{\partial {\bar{u}}_{j}}{\partial x_{i}}) - \bar{{u^{'}}_{i} {u^{'}}_{j}}] - M_{e x_{i}}

(2)

where $u_{i}$ is the Eulerian air velocity vector, $x_{i}$ is the position vector, $ρ_{a}$ stands for the air density, p for the air pressure, and $μ$ for the air molecular turbulent viscosity. Here, overbars and primes denote ensemble averages and fluctuations, respectively. $- \bar{u'_{i} u'_{j}}$ is related to the Reynolds stresses. $M_{e x_{i}}$ is the component of the interphase momentum exchange $M_{ex}$ vector in equation (14) in the i direction.

The conversation equations for the momentum and turbulence after discretization at a certain cell “C” can be generally written as

a_{c} ϕ_{c} = \sum_{n b} a_{n b} ϕ_{n b} + b

(3)

where $ϕ = [\begin{matrix} u_{i} \\ u_{j} \\ u_{k} \\ \tilde{ν} \end{matrix}]$ is a general variable representing the componential velocity in each axis direction and the eddy viscosity $\tilde{ν}$ in the S–A turbulence equation.

The total residual $R_{ϕ}^{N}$ for all cells at the Nth time step can be written as

R_{ϕ}^{N} = \sum_{cells} | \sum_{nb} a_{nb} ϕ_{nb} + b - a_{c} ϕ_{c} |

(4)

The convergence criteria for our calculation are the residuals of the continuity, momentum, and turbulence equations which are reduced by five orders of magnitude relative to the initial residuals. For the momentum and turbulence equations, that is

\frac{R_{ϕ}^{N}}{R_{ϕ}^{0}} < 10^{- 5}

(5)

where $R_{ϕ}^{0}$ and $R_{ϕ}^{0}$ are the residuals at the Nth time step and the first time step, respectively. And for the continuity equation, it is

\frac{\sum_{cells}^{N} ρ u \cdot n}{\sum_{cells}^{0} ρ u \cdot n} < 10^{- 5}

(6)

where $ρ u \cdot n$ is the net mass flux rate at the neighboring faces.

Discrete phase

Raindrop tracking algorithm

The raindrops in our study are assumed to be non-interacting, non-deforming, and non-evaporating spheres only subject to drag and gravity forces. For raindrops whose density is much greater than air density, the dimensional equations of motion for a particle can be written as

\frac{d u_{p}}{dt} = \frac{3 μ C_{D} R e_{p}}{4 ρ_{p} D_{p}^{2}} (u - u_{p}) + g

(7)

\frac{d x_{p}}{dt} = u_{p}

(8)

where $x_{p} = (x_{p}, y_{p}, z_{p})$ is the raindrop position in the Cartesian coordinate frame where x-axis is horizontal and positive backward, y-axis is vertical and positive upward, and z-axis is determined by the right-handed rule. $u_{p} = (u_{p}, v_{p}, w_{p})$ and $u = (u, v, w)$ denote raindrop particle and air velocity vectors, respectively. $g = (0, - g, 0)$ stands for gravity acceleration vectors, as well as $ρ_{p}$ and D_p for particle density and diameter, respectively. $R e_{p}$ is the particle Reynolds number defined by

R e_{p} = \frac{ρ_{a} D_{p} | u - u_{p} |}{μ}

(9)

The drag coefficient $C_{D}$ is given by¹¹

C_{D} = \frac{K_{1}}{R e_{p}} + \frac{K_{2}}{{Re}_{p}^{2}} + K_{3}

(10)

In this work, the second-order accurate implicit trapezoidal discretization is adopted to solve equations (7) and (8). The final iterative formulas are given as follows

u_{p}^{n + 1} = \frac{[u_{p}^{n} (1 - \frac{3 μ C_{D} R e_{p}}{8 ρ_{p} D_{p}^{2}} Δ t) + {\frac{3 μ C_{D} R e_{p}}{4 ρ_{p} D_{p}^{2}} Δ t} (u^{n} + \frac{Δ t}{2} u_{p}^{n} \cdot \nabla u^{n})]}{[(1 + \frac{3 μ C_{D} R e_{p}}{8 ρ_{p} D_{p}^{2}} Δ t)]}

(11)

x_{p}^{n + 1} = x_{p}^{n} + \frac{1}{2} Δ t (u_{p}^{n} + u_{p}^{n + 1})

(12)

where the superscript “n” denotes the current time step and “n + 1” the next time step.

Random walk dispersion model

In this study, the dispersion of raindrop particles due to turbulence in the air phase is taken into consideration. In equation (7), the air velocity can be represented as $u = \bar{u} + u' (t)$ , where $\bar{u}$ denotes the mean velocity and $u' (t)$ the instantaneous fluctuating air velocity. Then, particle turbulent dispersion can be predicted using the instantaneous air velocity along the particle trajectory. The instantaneous air velocity is determined by the Markov-chain random walk dispersion model^12,13 as

u'_{n + 1} = r_{u} u'_{n} + \sqrt{(1 - r_{u}^{2})} σ'_{u} ζ_{u} + (1 - r_{u}) T_{Lu} \frac{\partial τ'_{uv}^{2}}{\partial z}

(13)

v'_{n + 1} = r_{u} v'_{n} + \sqrt{(1 - r_{v}^{2})} σ'_{v} ζ_{v} + (1 - r_{v}) T_{Lv} \frac{\partial σ'_{v}^{2}}{\partial z}

(14)

w'_{n + 1} = r_{w} u'_{w} + \sqrt{(1 - r_{w}^{2})} σ'_{w} ζ_{w} + (1 - r_{w}) T_{Lw} \frac{\partial τ'_{wv}^{2}}{\partial z}

(15)

where $r_{i} = \exp (- Δ t / T_{Li}), i = u, v, w$ is the typical correlation for air turbulence. $σ'_{u}^{2} = \bar{{u'}^{2} (t)}$ , $σ'_{v}^{2} = \bar{{v'}^{2} (t)}$ , $σ'_{w}^{2} = \bar{{w'}^{2} (t)}$ , $τ'_{uv}^{2} = \bar{{u'}^{2} (t) v' (t)}$ , and $τ'_{wv}^{2} = \bar{{w'}^{2} (t) v' (t)}$ are time-averaged turbulence velocity variances for air. $ζ_{u}$ , $ζ_{v}$ , and $ζ_{w}$ are random number sequences with Gaussian distributions (0 for mean value, and 1 for variance). $T_{Lu}, T_{Lv}, and T_{Lw}$ are Lagrangian integral time scales for air velocity, which can be approximated as

T_{Li} = \int_{0}^{\infty} \frac{u'_{i} (t) u'_{i} (t + s)}{u'_{i}^{2}} ds

(16)

where s denotes the particle path.

Water film model

To simulate the interaction between raindrop and airfoil wall surface, the water film model in Stanton and Rutland¹⁴ is adopted in our study, which allows a single-component liquid drop to impinge upon a surface of arbitrary configurations and form a thin water film. The main physical processes that affect the liquid film are shown in Figure 1(a). The governing equations for water film can be referred to in our previous work.⁹

Figure 1.

(a) Mechanisms of mass, momentum, and energy transfers for the thin water film and (b) the four impingement regimes included in the film model.¹⁴

In our work, we do not consider the energy issue, since the current flight is a low-speed movement (130 mile/h, about 58 m/s). With regard to this, convective heat transfer from water film to air, heat conduction from wall to water film and raindrop evaporation are not considered. There are four regimes for determining the different interactions between raindrop and wall, including stick, rebound, spread, and splash, as shown in Figure 1(b). The criteria by which the regimes are partitioned are based on the raindrops’ Weber number defined by

W_{e} = \frac{ρ_{p} u_{r}^{2} D_{p}}{σ_{p}}

(17)

where $u_{r}$ is the raindrop relative velocity in the frame of the wall (i.e. $u_{r} = u_{p} - u$ ), and $σ_{p}$ is the raindrop surface tension.

The stick regime occurs when an impinging drop with extremely low impact energy adheres to the film surface in, approximately, a spherical form and is triggered at $W_{e} < 5$ . Next, the rebound regime occurs when the impinging drop with low impact energy bounces off the film using the lost energy of the air layer trapped between raindrop and water film and is triggered at $5 < W_{e} < 10$ . The third regime, spread, occurs at higher $W_{e}$ , where the drop merges with the liquid film upon impact at the transition condition of $10 < W_{e} < 18 . 0^{2} D_{p} (ρ_{p} / σ_{p})^{1 / 2} u_{r}^{1 / 4} f^{3 / 4}$ , where f is the frequency of the impinging drops. The final regime, splashing, occurs at high impact energy. In this circumstance, a raindrop with effectively high energy impacts airfoil surface and splashes and breaks up into a cloud of secondary droplets with a crown formed at the periphery. The transition criterion for the splashing regime is $18 . 0^{2} D_{p} (ρ_{p} / σ_{p})^{1 / 2} u_{r}^{1 / 4} f^{3 / 4} < W_{e}$ . In this study, u_r = 1.87 m/s and f = 1000 Hz. Therefore, the range of the raindrops’ Weber number characterizing the four regimes for this study is shown in Table 1.

Table 1.

Range of the raindrops’ Weber number characterizing the four regimes for this study.

Regime	Weber number (W_e)
Stick	W_e < 5
Rebound	5 < W_e < 10
Spread	10 < W_e < 3672
Splash	W_e > 3672

Interphase coupling

In our study, a two-way momentum coupling process is used to account for the interphase interaction between the air and raindrop phases. The momentum transfer from air to raindrop acts as a momentum sink in the air-phase momentum balance in the subsequent solution for the airflow field. It originates from the interphase drag and is explicitly accounted for in the Lagrangian particle trajectory (equation (7)), which is expressed as

M_{ex} = \sum \frac{3 μ C_{D} R e_{p}}{4 ρ_{p} D_{p}^{2}} (u - u_{p}) {\overset{\cdot}{m}}_{p} Δ t

(18)

where $M_{ex}$ is the momentum change of a particle as it passed through each control volume and is subtracted from the right-handed source terms of the air momentum equation to account for the effect of raindrop on air. The effect of airflow on the raindrop trajectory is naturally incorporated through the air velocity u in equation (7). ${\overset{\cdot}{m}}_{p}$ is the particle mass flow rate per unit volume.

Calculation process

In the process of rain simulation, first, the dry airflow field is calculated. Then, the raindrop trajectory is calculated based on the obtained air velocity field. Afterward, the momentum coupling “body force” term calculated by equation (18) is added to the source term of the air momentum equation to account for the raindrop effect on the air phase. So far, a loop of calculation is completed. This loop is repeated until both the air and raindrop phases reach a steady state, as shown in Figure 2. The convergence criteria for our calculation are as follows: first, the residuals of the continuity, momentum, and turbulence equations of the continuous phase are reduced by five orders of magnitude relative to the initial residuals, which denotes the convergence of the air phase; and second, the final results of lift and drag coefficients were independent of the raindrop number in the model; that is, further increase in particle number in the simulation did not remarkably affect the final results, which is deemed as the equilibrium of the raindrop flowfield, since raindrops are ejected from the upstream side and exit from the downstream outlet.

Figure 2.

Flowchart of the calculation process.

Results and discussions

The first step in our study is the modeling of the ice shapes obtained in the Ohio State University subsonic wind tunnel.¹⁵ In this study, the No. 2 NACA 0012 airfoil is adopted because of its available geometry and experimental aerodynamic data from the experiment. The airfoil chord length is 0.5334 m. The icing conditions that produced this ice shape are angle of attack of 4°, free-stream velocity of 58 m/s (corresponding to Reynolds number of 1.5 × 10⁶), median volume diameter of 20 µm, liquid water content (LWC) of 2.1 g/m³, total temperature of 265 K, and 5 min spray time.

The two-dimensional (2D) grids are used in our study for modeling the computational region around the clean and iced airfoils consisting of two blocks, where the inner and outer blocks have an overlapping interface at 0.6 c, as shown in Figure 3. Both blocks have a C-type structured grid, while the inner one has a wake cut downstream of the trailing edge and a much denser distribution of grid points than the outer block. Both blocks have a boundary at 10-c upstream and 15-c downstream from the leading edge. This mesh strategy of constructing two-block grids not only allows for easier control of the grid generation process and better quality of the resulting grid, especially near the complex iced surfaces, but also reduces the mesh size of the region far from the airfoil surface which is not necessarily large for a high equality of acquired flowfield simulation, thus saves much computational resource.

Figure 3.

Two-block grid strategy.

Grid independence of the solutions is assessed by constructing a new grid with additional refinement for the clean airfoil configuration. As shown in Figure 4, the simulated lift coefficient C_L obtained from the clean airfoil grid resolution with a total of 340 nodes on the surface of the airfoil and 80 nodes in the normal direction of the inner block, with a first layer thickness of 2.0 × 10⁻⁵ c, guarantees a grid-independent solution for, as an 18% and 25% refinement in the streamwise and normal directions, respectively, changing the results by less than 1%. Besides, the simulated results from this grid solution shows a maximum error of only about 1.72% compared to the wind tunnel experimental results.¹⁵

Figure 4.

Grid sensitivity test in the normal direction: lift.

For rain simulations, a heavy rain condition of LWC of 30 g/m³, free-stream velocity of 58 m/s (corresponding to Reynolds number of 1.5 × 10⁶), and a range of raindrop diameter from 0.1 to 0.9 mm which is equivalent to an arithmetic mean diameter of 0.4672 mm are selected in this study.

Figure 5(a) and (b) presents the numerical results of lift and drag coefficients for the airfoil in the four cases (i.e. clean airfoil in dry air, clean airfoil in rain, iced airfoil in dry air, and iced airfoil in rain). On one hand, both the heavy rain and the ice bring significant aerodynamic performance degradations to the clean airfoil, although it can be seen that ice induces a heavier aerodynamic penalty than heavy rain due to a very severe glazed ice accretion selected. Thus, an adverse effect of heavy rain equivalent to that of icing can be expected, if a lighter icing condition is encountered. On the other hand, rain brings a further consistent decrease in lift as well as increase in drag to the iced airfoil at angles of attack below 9°, although no significant changes occurred at higher angles of attack. Regarding the effect of rain, a maximum percentage decrease in C_L by 24% and 7% and increase in C_D up to 108% and 6% are caused by rain for the clean and iced airfoils, respectively. Considering aircrafts’ cruise (usually at 0°–2°), take off, and landing (usually at 4°–8°) at low angles of attack, these findings can be of much importance to guide both engineering design and pilot flight training.

Figure 5.

(a) Lift comparison for iced airfoil in rain and (b) drag comparison for iced airfoil in rain.

Previous studies have found that rain can cause premature separations to airfoil boundary layer.^8,9 The streamlines around both the clean and iced NACA 0012 configurations in the dry and heavy rain conditions at α = 12° are presented through Figure 6(a)–(d). From Figure 6(a), it can be seen that there is little airflow separation near the trailing edge of the clean configuration. However, clear separation can be seen for the clean configuration in the heavy rain condition at the same angle of attack in Figure 6(b), suggesting a premature boundary-layer separation induced by heavy rain to this airfoil. Besides, the leading-edge ice shape causes more severe boundary-layer separation in a large region over the trailing edge, as shown in Figure 6(c). It can be seen from Figure 6(d) that heavy rain causes a larger separation zone above the trailing edge of the iced airfoil, making it more severe. To quantify the effect of rain on the separation characteristics, a plot of the shear stress along the horizontal direction, x-shear stress, as a function of the normalized chord position for the upper wall of the airfoil is presented in Figure 7. For the clean configuration in dry condition, no significant separation can be quantified. While in rain, a clean separation emerges with a separation point starting at about 80% chord. For the iced configuration, rain does not bring any advance in the separation point, since ice has caused a full separation on the whole chord of the upper airfoil wall. Despite that it can be seen that there exists a small decrease in the magnitude of the shear stress at the leading edge and the front of the airfoil, indicating a decreased friction drag for the iced airfoil in rain. It further suggests that the boundary-layer velocity of the iced airfoil is decelerated by the splashback of raindrop at this time, since friction drag is reduced with decreasing boundary-layer air velocity, which will be presented more in detail in the next discussion.

Figure 6.

Streamlines around the clean airfoil (a) in dry air at α = 12° and (b) in the rain condition at α = 12°. Streamlines around the iced airfoil (c) in dry air at α = 12° and (d) in the rain condition at α = 12°.

Figure 7.

x-shear stress on the upper wall of the airfoil.

As raindrops strike on the airfoil surface, some of them are splashed back and break up into secondary drops. The acceleration of these splashed-back droplets by the airflow field acts as a momentum sink for the boundary layer, resulting in a decrease in the air velocity field. De-energization of the boundary layer can induce a loss of lift and L/D, increase in drag, premature stall, transition and separation.⁹ Figure 8 shows the boundary-layer velocity profiles at four chordwise positions on the upper surface of the iced airfoil in the heavy rain condition at α = 6°. Near the leading edge, more severe deceleration of the boundary layer by splashed-back raindrops relative to the boundary without rain can be seen. Moving downstream, the boundary layer appears to recover.

Figure 8.

Boundary-layer velocity profiles at four chordwise positions on the upper surface of the iced airfoil at α = 6°. Position is measured from the leading edge. Dimensionless velocity $\bar{u} = u / u_{\infty}$ ( u _∞ denotes free-stream velocity) is plotted versus dimensionless distance $\bar{y} = y / c$ ( c denotes airfoil chord length) from the airfoil surface.

Conclusion

In this article, a post-ice-accretion computational fluid dynamics (CFD) numerical simulation is performed on the aerodynamic performance of an iced NACA 0012 airfoil under a simulated heavy rain condition. This study is conducted to obtain some quantitative data of the coupling effects of heavy rain and ice accretion on airfoil aerodynamics and to provide a new insight into the aerodynamic penalties caused by such two adverse weather conditions as icing and heavy rain that may lead to a control upset. Aerodynamic performance degradations were observed for the clean airfoil in the conditions of heavy rain and ice accretion. Heavy rain can cause premature boundary-layer separation to both the clean and iced airfoils. Boundary-layer momentum losses by splashed-back raindrops were also observed on the iced NACA 0012 airfoil, which is a typical consequence of the adverse effects of rain imparting airfoils/wings. As the requirement for commuting airplanes is rising, our study can be significant for future design of commuting airplanes and pilot training to safely fly in such adverse weather conditions of icing and heavy rain.

Footnotes

Appendix 1

Academic Editor: Roslinda Nazar

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Lynch

Khodadoust

Effects of ice accretions on aircraft aerodynamics. Prog Aerosp Sci 2001; 37: 669–767.

Bragg

Broeren

Blumenthal

LA.

Iced-airfoil aerodynamics. Prog Aerosp Sci 2005; 41: 323–362.

Cao

Effects of rainfall on aircraft aerodynamics. Prog Aerosp Sci 2014; 71: 85–127.

Koomullil

Thompson

Soni

BK.

Iced airfoil simulation using generalized grids. Appl Numer Math 2003; 46: 319–330.

Broeren

Bragg

MB.

Flowfield measurements about an airfoil with leading edge ice shapes. J Aircraft 2006; 43: 1226–1234.

Cao

Chen

Sheridan

Flowfield simulation and aerodynamic performance analysis of complex iced aerofoils with hybrid multi-block grid. Proc IMechE, Part G: J Aerospace Engineering 2008; 222: 416–422.

Yeom

Chang

Baek

SW.

Numerical prediction of airfoil characteristics in a transonic droplet-laden air flow. Int J Heat Mass Tran 2012; 55: 453–469.

Valentine

Decker

AR.

A Lagrangian-Eulerian scheme for flow around an airfoil in rain. Int J Multiphas Flow 1995; 21: 639–648.

Cao

Aerodynamic study of aerofoil and wing in simulated rain environment via a two-way coupled Eulerian-Lagrangian approach. Aeronaut J 2014; 118: 643–668.

10.

Leonard

Mokhtari

ULTRA-SHARP nonoscillatory convection schemes for high-speed steady multidimensional flow. Technical report. NASA TM 102568, April 1990. Washington, DC: NASA.

11.

Morsi

Alexander

AJ.

An investigation of particle trajectories in two-phase flow systems. J Fluid Mech 1972; 55: 193–208.

12.

Legg

Raupach

MR.

Markov-chain simulation of particle dispersion in inhomogeneous flows: the mean drift velocity induced by a gradient in Eulerian velocity variance. Bound: Lay Meteorol 1982; 24: 3–13.

13.

Legg

BJ.

Movement of plant pathogens in the crop canopy [and discussion]. Philos T Roy Soc B 1983; 302: 559–574.

14.

Stanton

Rutland

CJ.

Multi-dimensional modeling of thin liquid films and spray-wall interactions resulting from impinging sprays. Int J Heat Mass Tran 1998; 41: 3037–3054.

15.

Bragg

MB.

An experimental study of the aerodynamics of a NACA 0012 airfoil with a simulated glaze ice accretion. Technical report. NASA-CP-179897, November 1986. Washington, DC: NASA.

Numerical simulation of airfoil aerodynamic performance under the coupling effects of heavy rain and ice accretion

Abstract

Keywords

Introduction

Methodology description

Continuous phase

Discrete phase

Raindrop tracking algorithm

Random walk dispersion model

Water film model

Interphase coupling

Calculation process

Results and discussions

Conclusion

Footnotes

Appendix 1

Declaration of conflicting interests

Funding

References