# Numerical Simulation of Airfoil Aerodynamic Penalties and Mechanisms in Heavy Rain

Authors: Zhenlong Wu; Yihua Cao; M. Ismail
Journal: International Journal of Aerospace Engineering
 (2013)

---

## Abstract

Numerical simulations that are conducted on a transport-type airfoil, NACA 64-210, at a Reynolds number of2.6×106 and LWC of 25 g/m3 explore the aerodynamic penalties and mechanisms that affect airfoil performance in heavy rain conditions. Our simulation results agree well with the experimental data and show significant aerodynamic penalties for the airfoil in heavy rain. The maximum percentage decrease in CL is reached by 13.2% and the maximum percentage increase in CD by 47.6%. Performance degradation in heavy rain at low angles of attack is emulated by an originally creative boundary-layer-tripped technique near the leading edge. Numerical flow visualization technique is used to show premature boundary-layer separation at high angles of attack and the particulate trajectories at various angles of attack. A mathematic model is established to qualitatively study the water film effect on the airfoil geometric changes. All above efforts indicate that two primary mechanisms are accountable for the airfoil aerodynamic penalties. One is to cause premature boundary-layer transition at low AOA and separation at high AOA. The other occurs at times scales consistent with the water film layer, which is thought to alter the airfoil geometry and increase the mass effectively.


---

## Body

## 1. Introduction

The aerodynamic penalty of aircraft flight through heavy rain has been deemed to be a critical cause in many severe aviation accidents. The Eastern Flight 066 accident at Kennedy Airport (NTSB, 1976) is a telling example, though the factor of heavy rain was not taken into consideration at that time [1]. Three years later, another Flight named 693, a Boeing 727-25 civil airplane, suffered from an intense rainfall associated with wind shears on its eventual routine to the Atlanta International Airport [2]. Several severe aviation accidents in 1981 aroused people’s consciousness of the seriousness of rain influence on aircraft flight [3].Investigation of rain effect on aircraft flight was begun with the wind tunnel test, and the earliest was conducted by Rhode in 1941 [4]. It dealt with the situation of an aircraft encountering heavy rain at moderate cruising altitude of about 5000 ft and concluded that, the heavy rain exposure time is not sufficient to force the aircraft to the ground. In 1982, Haines and Luers did a research concerning the frequency and intensity of very heavy rains and their effects on a landing aircraft [1]. In 1987, Hansman and Criag compared the aerodynamic performance degradation of NACA 64-210, NACA 0012, and Wortman FX 67-K170 airfoils under the low Reynolds numbers in heavy rain conditions and explored the various mechanisms underlying by forcing boundary-layer to transition [5]. In other similar wind tunnel experiments, laminar flow airfoils were also found to experience performance degradation approximately equivalent to that caused by tripping the boundary layer to turbulence [6–8]. In 1992, Bezos et al. determined the severity of rain effect, the aerodynamic penalty over a range of rain intensities, and the importance of surface tension interactions of water as a scaling parameter [9]. Thompson and his team examined another NACA 4412 airfoil in moderate rain in wind tunnel. They primarily placed emphasis on the correlation of surface-film behavior including rivulet formation [10]. Subsequently, they went on a further examination on the aerodynamic efficiency of the same airfoil in moderate rain [11]. Comparisons with different flow patterns showed the aerodynamic degradation depended on the location of rivulet formation and the diameter of these rivulets. The latter factor was found to be more important to aerodynamic performance.Numerical simulation approach was introduced and developed with the development of computer technology. In 1995, Valentine and Decker studied the NACA 64-210 airfoil aerodynamic performance [12] and the track of raindrops in flow over the airfoil [13] by numerical simulation. In 1999, Thompson and Marrochello calculated the location of the oneset of rivulet formation in the surface-water flow over a wing with a NACA 4412 airfoil and compared the results with wind-tunnel experiments [14]. In 2003, Wan and Wu also conducted the numerical simulation of heavy rain effect on airfoil [15]. The water film layer and vertical rain mass flow rate on the airfoil upper surface were added, thus increasing the airfoil roughening effects. In 2010, Wan and Pan studied the cruise and high-lift NACA 64-210 airfoil aerodynamic efficiency in heavy rain via a two-phase flow approach [16]. Later, he reinvestigated the high-lift NACA 64-210 with the consideration of proper modeling of discrete water droplets, shear flow between airfoil elements [17], and studied aerodynamic performance of a 3D blended-wing-body aircraft under severe rain through a two-phase flow approach [18]. Zhang and Cao and Ismail et al. studied aerodynamic characteristics of the NACA64-210 and NACA 0012 airfoils in rain and preliminarily explored the mechanism [19, 20].As is presented above, present research methods of rain effect on aircraft aerodynamic performance are mainly wind-tunnel test and numerical simulation. However, for the former method in one hand, the result of model test in wind tunnel may not be directly used for large-scale aircraft due to the complexity of multiphase environment. For example, the Weber number of water cannot be maintained as uniform meanwhile. On the other hand, if the scale of model is decreased, then the testing velocity must be increased proportionally to a uniform Reynolds number; thus the wind tunnel test will be much more expensive and complex to manipulate. So numerical simulation method may as well be accepted as a good approach to deal with the involved problem.In numerical simulations, two approaches have been used to model multiphase flows, that is, the Eulerian approach and the Lagrangian approach. They have been reviewed by Valentine and Decker [13]. The Eulerian approach treats the continuous fluid phase and the dispersed particle phase as continuum, while the Lagrangian approach solves time-averaged Navier-Stokes equations for the continuous fluid phase first and then integrates the Lagrangian motion equation for the dispersed phase; this model is called discrete phase model (DPM). There are two models including a one-way coupled model and a two-way coupled model in the Lagrangian approach. The former model assumes that the particles motion is affected by the continuous phase, but the continuous phase is not affected by the presence of the dispersed phase. The latter takes the two-way exchange of mass, momentum, and energy between the two phases into consideration.Though rain effect on airfoil aerodynamic performance and some mechanisms have been exposed by wind-tunnel tests, little work of numerically studying the mechanisms has been done. For example, Wan and others mainly discussed airfoils aerodynamic performance in rain conditions at low angles of attack. The rain effect of airfoil at high angles of attack is little concerned, neither is the mechanism researched.The present study uses the discrete phase model (DPM) in Fluent to study the typical commercial transport airfoil NACA64-210 aerodynamic performance in heavy rain, but mainly it places emphasis on the mechanism exploration. The raindrops in our study are assumed to be nonevaporating, nondeforming spheres. The objective of our study is fourfold: first, to determine the aerodynamic penalties of the airfoil over a wider range of attack angle; second, to investigate the premature boundary-layer transition and separation mechanism leading to the dynamic penalty at low and high angles of attack, respectively, using a CFD method; third, to track the particulate trajectory at various angles of attack; fourth, to acquire the water-film-layer influence on the geometric changes of airfoil in different rain rates and at different angles of attack.


## 2. Numerical Approach

### 2.1. Fluid Phase

Previous wind-tunnel experiments performed by Hansmann and Barsotti [8] suggested that premature boundary-layer transition causes the aerodynamic penalty of natural laminar airfoils at the low Reynolds numbers in heavy rain conditions. To investigate this hypothesis, Hansman Jr. and Craig [5] placed boundary-layer transition elements on the suction and pressure surfaces of airfoil to force airfoil boundary layer to transition. Trip strips 0.25 in. wide and madeup of sand grains ranging in diameter from 0.025 to 0.040 inch were placed at the 5, 25, 50, or 75% chordwise station on the top and bottom surface of airfoil. Trip strips on the lower surface of the airfoils resulted in minimal performance changes, whereas forcing boundary-layer transition on the upper surface resulted in fairly significant performance changes. Therefore, the location of the lower trip strip was generally fixed at 5% chord, and the upper trip strip location was varied. To investigate the hypothesis from a perspective of computational fluid dynamics, we develop a CFD (Computational Fluid Dynamics) method to emulate the forced boundary-layer transition, that is, to give a normal and outward increment as  modeling the trip strip, like a step, to the specified-position section of the airfoil, as shown in Figure 1(d). This is a new method, so we name it numerical boundary-layer-tripped technique in this paper. We also put the trip at 5, 25, 50, and 75% chordwise stations on the top surface and at 5% on the bottom surface of a NACA 64-210 airfoil for our study. The width and height of the increment have the same ratio to the chord length as that of the trip in Hansman’s experiment.(a) O-type computational domain of meshes for both the original and numerically tripped NACA 64-210 airfoils; (b) near mesh of the original airfoil; (c) near mesh of the airfoil tripped at 5% chordwise station; (d) close view of the tripped position on the 5% chordwise tripped airfoil.
(a)
(b)
(c)
(d)The computational domain consists of an extruded O-type grid around the NACA 64-210 airfoil. The first step in computational fluid dynamics is to examine the grid dependency on numerical results. Generally, the more nodes are distributed, the more accurate solution will be acquired, however the more expensive computational memory and time will be required. Table1 shows the effect of mesh size on lift and drag coefficients for the original airfoil at angles of attack (AOA) of 0° and 10° in the dry and rain conditions. It indicates that the aerodynamic penalties due to rain are not greatly affected by the number of grid cells. Considering our computational limits, we choose the mesh with 12,000 grid cells for a further study, as shown in Figures 1(a) and 1(b). The height of the first cell adjacent to the original airfoil surface is set to 1 × 10−5 m, and the stretching ratio of the mesh is 1.1. The mesh of 5% chordwise position tripped NACA 64-210 airfoil is shown in Figures 1(c) and 1(d), and the number of cells is 34, 632. The height of the first cell adjacent to the tripped airfoil surface is also set to 1 × 10−5 m as well as the stretching ratio of 1.1. The origin is set at the airfoil leading edge with the X-axis pointing to the right and Y-axis upward. The circular outer boundary of the two O-type computational domains is set as the velocity-inlet condition and the airfoil surface as no-slip wall condition. The chord length is 0.762 m and the free-stream Reynolds number of air is set to 2.6 × 106 equal to a dynamic pressure of 30 psf in order to be consistent with parameters of theory of wing section [9]. It is stated here that the results of which airfoil is not stated explicitly in all below figures are about the original NACA 64-210 airfoil, including the experimental and the numerical ones.Table 1
Lift and drag coefficients for different sizes of the original airfoil mesh.
Mesh size
C
L AOA 0°
C
L AOA 10°
LWC 0 g/m3
25 g/m3
ΔCL
LWC 0 g/m3
25 g/m3
ΔCL
8
0
×
100
0.161
0.155
0.006
1.081
1.029
0.052
100
×
100
0.157
0.152
0.005
1.077
1.028
0.049
100
×
120
0.156
0.151
0.005
1.076
1.027
0.049
120
×
120
0.155
0.151
0.004
1.075
1.027
0.048
C
D AOA 0°
C
D AOA 10°
LWC 0 g/m3
25 g/m3
ΔCD
LWC 0 g/m3
25 g/m3
ΔCD
80
×
100
0.020
0.022
0.002
0.078
0.082
0.006
100
×
100
0.019
0.020
0.001
0.075
0.080
0.005
100
×
120
0.019
0.020
0.001
0.074
0.080
0.006
120
×
120
0.019
0.020
0.001
0.074
0.079
0.005
LWC: liquid water content; AOA: angle of attack;CL: lift coefficient; CD: drag coefficient.The incompressible air flow field is solved by FLUENT, a common commercial flow field solver, the details of which can be referred to in the help literature [21] and will not be repeated here. For a Reynolds number of 2.6 × 106, the flow characteristic is considered as turbulent, so turbulence model is added to solve the Navier-Stokes equations. During the calculation of the original and tripped airfoils in dry condition, the steady pressure-based solver is chosen, of which the segregated SIMPLE algorithm is adopted to discretize the pressure-velocity coupling term. The pressure term uses second-order scheme, and the QUICK scheme is used in the momentum term discretization.Airfoil aerodynamic performance is measured by lift and drag coefficients in this research, which are defined, respectively as follows:(1)CL=L(1/2)ρav∞2c,CD=D(1/2)ρav∞2c,
where CL is the lift coefficient and CD is the drag coefficient, L is the lift, D is the drag, ρa is the density of air, v∞ is the air free-stream velocity, and c is the chord length of the airfoil of interest.


### 2.2. Particulate Phase

#### 2.2.1. Scaling of Rain Model

To study the heavy rain effect, first of all, it is necessary to measure the intensity and frequency of heavy rain. To study the heavy rain effect, first of all, it’s necessary to measure the intensity and frequency of heavy rain. Usually the rainfall rate,R in millimeter per hour or the Liquid Water Content, LWC in gram per cubic meter is chosen to categorize different intensities of rain. A rainfall of 100 mm/h or greater is often deemed as heavy. A rainfall of 100 mm/h or greater is often deemed as heavy.The LWC can be written as a function ofN(D) as follows:
(2)LWC=∫0+∞ρwπ6D3N(D)dD,
where ρw is the density of water. Integrating the above formula, we may attain the correlation of LWC and R by
(3)LWC=0.054R0.84.Subsequently, it is necessary to establish the size distribution of water droplets under different rain rates. Many researchers like Best [22], Ulbrich [23], and so on have established various raindrop size distribution formulas for various situations. Marshall and Palmer developed the classic formula of drop size distribution in 1948 based on massive experimental data [24]. It is shown as follows:
(4)N(D)=N0EXP(-ID)(0≤D≤Dmax),
where N(D) (m−3mm−1) is the number density of spherical raindrops of diameter D (mm) per cubic meter of air, Dmax is the maximum drop diameter, and N0 and I (mm−1) are parameters of N(D) and have different values for different types of rain. For storm-type heavy rainfall, I varies with rainfall rate R as I=3R-0.21, and N0 has the constant value: N0=1400 m−3mm−1 [25] correspondingly.Here, it is assumed that raindrops have been with uniform velocity before hitting the aircraft surface, that is, without acceleration. So it is important to determine the terminal velocity of raindrops. It has been developed by Markowitz [26] as follows:
(5)V(D)=9.58{1-EXP[-(D1.77)1.147]},
where V(D) is the terminal velocity. A correction for it aloft is given by Markowitz as
(6)V(D)=V0(D)(ρ0ρa)0.4,
where V0(D) is the terminal velocity consistent with the density of air aloft ρ0.


#### 2.2.2. Wall-Film Model

In our study, the wall-film model in Fluent [21] is mainly adopted to model the interaction of particle and wall surface. It allows a single-component liquid drop to impinge upon a boundary surface of arbitrary configuration and form a thin liquid film. The major physical processes that affect the liquid film include mass and momentum contributions to the film thanks to drop impingement, droplet splashing effects, evaporation, shear forces on the film, dynamic pressure effects, gravity driven flow, convective heat and mass transfer, flow separation, and sheet breakup, as shown in Figure 2(a). In present study, we ignore the film evaporation to simplify our solution, so it is unnecessary to consider the effects of the thin liquid film on the air flow.(a) Mechanisms of mass, momentum, and energy transfers for the wall-film model; (b) the four impingement regimes included in the model; (c) simplified decision for wall interaction criterion.
(a)
(b)
(c)The main assumptions for the film model are as follows.(i)
The layer is thin, less than 500 microns in thickness due to the assumption of a linear velocity profile in the film.(ii)
The temperature in the film particles changes relatively slow due to the use of an analytical integration scheme.(iii)
The film temperature is always below the boiling temperature for the liquid.(iv)
Film particles are assumed to be in direct contact with the wall surface and the heat transfer from the wall to the film occurs through conduction.The wall interaction regimes are calculated for a drop-wall interaction based on local information. The four regimes including stick, rebound, spread, and splash are based on the impact energy and wall temperature, as shown in Figures2(b) and 2(c) (Tb is the liquid boiling temperature and Tw is the wall face temperature). Below the liquid boiling temperature, the impinging droplet can either stick, spread, or splash, while above the boiling temperature, the particle can either rebound or splash. As to our case of which the temperature is below the boiling point, particles stick, spread, and splash, resulting in aerodynamic efficiency degradation of an aircraft.The criteria by which the regimes are partitioned are based on the impact energy and the boiling temperature of the liquid. The impact energyE is defined by
(7)E2=ρwVr2Dσw(1min(h0/D,1)+δbl/D),
where ρw is the water density, Vr is the particle relative velocity in the frame of the wall (i.e., Vr=Vp-Vwall), σw is the water surface tension, and h0 is the total height of the wall film. δbl denotes the thickness of the boundary layer and is defined by
(8)δbl=μDρwVr.The sticking regime is applied when the dimensionless energyE is less than 16 and the particle velocity is set equal to the wall velocity. In the spreading regime, the probability of the drop having a particular direction along the surface is given by an analogy of an inviscid liquid jet with an empirically defined radial dependence for the momentum flux. If the wall temperature is above the boiling temperature of the liquid, impingement events below a critical impact energy (E) result in the particles rebounding from the wall. Splashing occurs when the impingement energy is above a critical energy threshold, defined as Ecr = 57.7. Besides, in our study, we sample a cumulative probability distribution function (CPDF), which is acquired from the Weibull distribution function and fitted to the data from Mundo et al. [27], to determine the different diameter of each splashed parcel. The equation of the cumulative probability distribution function can be expressed as
(9)pdf(diDp)=2diDp2exp[-(diDp)2]
and it represents the probability of finding drops of diameter di in a sample of splashed raindrops.Bilanin [28] has investigated the evaporation of the particles near the surface and found that evaporation does not affect the airfoil aerodynamic efficiency, so in our study we ignore the vaporization of water film. Since the film particle vaporization is ignored, only the momentum and energy conservation equations remain in the conservation equations for wall-film particles.Conservation equation of momentum is as follows:(10)ρahdV⃑pdt+h(∇spf)c=τgt⃑g+τwt⃑w+P⃑˙imp,c-M˙imp,cV⃑p+F⃑˙n,c+ρah(g⃑-a⃑w),
where c denotes the current face where the particle locates, h denotes the current film height where the particle locates, ∇s is the gradient operator as to the surface, pf is the film pressure, τg and τw are, respectively, the magnitude of the shear stress of the air flow on the film surface and the magnitude of the stress of the film exerted by the wall. t⃑g and t⃑w are the unit vectors in the direction of the relative motion of the film and the wall, P⃑˙imp,c denotes the impingement pressure on the film surface, M˙imp,c is the impingement momentum source, F⃑˙n,c is the force to keep the film on the surface, and ρah(g⃑-a⃑w) is the body force term in which a⃑w is the acceleration of the relative wall motion.Conservation equation of energy is as follows:(11)mpCpdTpdt=κAph(Tw-Tp)+hfAp(Ta-Tp),
where mp is the particle mass, Cp is the equilibrium concentration of the rain droplet, κ is the liquid thermal conductivity, h is the film height where the particle locates, hf is the film heat transfer coefficient, Ap is the area of the film particle, Tp, Tw, and Ta, respectively, denote the temperature of the particle, wall, and air. The first term on the right-hand side denotes the thermal conduction from the wall to the film particle, and the second term denotes the thermal convection from the top surface to the film particle.


#### 2.2.3. Mathematic Modeling for Water Film

In this part, a numerical model for the water film based on [1] is adopted to calculate water film thickness by airfoil location and rainfall rate.  In each grid cell, the film thickness is resulted from water remaining after the gains and losses due to the incoming rain and film flow. After a time, an equilibrium thickness is completed in each cell. Here we will mainly examine the effect of water film on the airfoil geometric changes qualitatively, so for simplicity the pressure gradient and thermodynamic influences are neglected. It is assumed that the water film develops normally to the airfoil surface, as shown in Figure 3. The main program is acquired and shown as follows.Figure 3
Illustration of the water film development.The raindrop local collection efficiency within every grid cell, symbolized asβ(i), can be expressed as
(12)β(i)=n(i)n0Ymax-Ymins(i),
where n(i) is the number of raindrops within cell i, n0 is the total number of raindrops which can have an impact on the airfoil surface. Ymax and Ymin represent the maximum and minimum vertical position, respectively. s(i) is the arc length of cell i and  is shown as follows:
(13)s(i)=(x(i+1)-x(i))2+(y(i+1)-y(i))2,
where x(i) and y(i) represent position of cell i in X and Y direction, respectively. Consider
(14)β¯=∑1Nβ(i)N,
where β¯ is the average collection efficiency at the airfoil surface, N is the total number of cells adjacent to the airfoil surface.The air shear stress within every cell,τ(i) can be expressed as
(15)τ(i)=12Cf(i)×ρa×v¯2,
where Cf(i) is the air coefficient of friction and v¯ is the average velocity of cells at the airfoil surface as follows:
(16)v¯=∑1Nv(i)N,
where v(i) is the velocity of air in cell i.The resulting film is a balance between water runback and rain water reception. For the film in each cell, the equilibrium condition of the mass balance can be expressed as(17)dmidt=ρwA∫z=0x=iΔxz=hudz-ρwA∫z=0x=(i+1)Δxz=hudz+∂mi∂t,
where u=τz/μw is the film velocity at height z, μw is the molecular viscosity of water. The first and second terms at the right side represent, respectively, mass flux into the ith cell (x=iΔx) of Area A and out of the cell (x=(i+1)Δx) due to film flow. The third term represents mass flux into the ith cell due to rain.The change rate of raindrop mass in celli, (dm/dt)(i), is calculated from (12):
(18)dmdt(i)=12ρw×τ(i)[((N-(i+1))h0(i)N)2-((N-i)h0(i)N)2]·1μws(i)+LWC×β(i)×v∞,
where h0(i) is the initial thickness of water film in cell i and set as
(19)h0(i)={0,β(i)=0,0.0005,β(i)≠0.After all, a new water film thickness in celli, h(i), results from [1]:
(20)h(i)=h0(i)+1ρwdmdt(i)Δt,
where Δt is the time step.Experimental research shows that, once the water droplets impact the thin water film layer, scallops will form, and the water remaining in the scallops will form convex film layers, thus increasing the roughness of the original water film. Based on Macklin’s research [29], the present program assumes the convex film layers to be cylindrical so as to take the water film roughness into consideration.From Driling [30], the correlation equations for equivalent sandgrain roughness developed from experimental data can be expressed as the function of the average height of the convex film; that is
(21)ksk¯={0.0164Λ3.7,Λ≤4.93,139Λ-1.9,Λ≥4.93,
where ks is equivalent sandgrain roughness, k¯ is the average height of convex film, Λ is the correlating parameter, and both can be referred to [1].The original NACA 64-210 airfoil in wet condition of LWC = 25 g/m3 is simulated and the results are compared with the wind-tunnel experimental ones in [9]. The wind tunnel experiment was conducted by Bezos et al. in the NASA Langley 14- by 22-Foot Subsonic Tunnel in 1992 to determine the aerodynamic penalty associated with a simulated heavy rain encounter. The model was comprised of a NACA 64-210 airfoil section with a chord of 2.5 ft, a span of 8 ft and was mounted on the tunnel centerline between two large endplates. The rain simulation system manifold, which was located 10 chord lengths upstream of the model, produced liquid water contents ranging from 16 to 46 g/m3. Aerodynamic measurements in and out of the simulated rain environment were obtained to measure the airfoil aerodynamic coefficients. In our rain calculations, the pressure-based unsteady solver with the first-order implicit scheme is adopted; the calculation condition of rain is shown in Table 2.Table 2
Calculated rain condition.
Variable
Value
LWC/(g/m3)
25
Raindrop arithmetic mean diameter/(mm)
1
Air freestream velocity/(m/s)
48.4632
Raindrop vertical velocity/(m/s)
−4
Air density/(kg/m3)
1.225
Airfoil chord length/(m)
0.762
Air freestream Reynolds number
2.6
×
10
6
Reference pressure/(Pa)
1.01325
×
10
5
Reference temperature/(K)
273


## 2.1. Fluid Phase

Previous wind-tunnel experiments performed by Hansmann and Barsotti [8] suggested that premature boundary-layer transition causes the aerodynamic penalty of natural laminar airfoils at the low Reynolds numbers in heavy rain conditions. To investigate this hypothesis, Hansman Jr. and Craig [5] placed boundary-layer transition elements on the suction and pressure surfaces of airfoil to force airfoil boundary layer to transition. Trip strips 0.25 in. wide and madeup of sand grains ranging in diameter from 0.025 to 0.040 inch were placed at the 5, 25, 50, or 75% chordwise station on the top and bottom surface of airfoil. Trip strips on the lower surface of the airfoils resulted in minimal performance changes, whereas forcing boundary-layer transition on the upper surface resulted in fairly significant performance changes. Therefore, the location of the lower trip strip was generally fixed at 5% chord, and the upper trip strip location was varied. To investigate the hypothesis from a perspective of computational fluid dynamics, we develop a CFD (Computational Fluid Dynamics) method to emulate the forced boundary-layer transition, that is, to give a normal and outward increment as  modeling the trip strip, like a step, to the specified-position section of the airfoil, as shown in Figure 1(d). This is a new method, so we name it numerical boundary-layer-tripped technique in this paper. We also put the trip at 5, 25, 50, and 75% chordwise stations on the top surface and at 5% on the bottom surface of a NACA 64-210 airfoil for our study. The width and height of the increment have the same ratio to the chord length as that of the trip in Hansman’s experiment.(a) O-type computational domain of meshes for both the original and numerically tripped NACA 64-210 airfoils; (b) near mesh of the original airfoil; (c) near mesh of the airfoil tripped at 5% chordwise station; (d) close view of the tripped position on the 5% chordwise tripped airfoil.
(a)
(b)
(c)
(d)The computational domain consists of an extruded O-type grid around the NACA 64-210 airfoil. The first step in computational fluid dynamics is to examine the grid dependency on numerical results. Generally, the more nodes are distributed, the more accurate solution will be acquired, however the more expensive computational memory and time will be required. Table1 shows the effect of mesh size on lift and drag coefficients for the original airfoil at angles of attack (AOA) of 0° and 10° in the dry and rain conditions. It indicates that the aerodynamic penalties due to rain are not greatly affected by the number of grid cells. Considering our computational limits, we choose the mesh with 12,000 grid cells for a further study, as shown in Figures 1(a) and 1(b). The height of the first cell adjacent to the original airfoil surface is set to 1 × 10−5 m, and the stretching ratio of the mesh is 1.1. The mesh of 5% chordwise position tripped NACA 64-210 airfoil is shown in Figures 1(c) and 1(d), and the number of cells is 34, 632. The height of the first cell adjacent to the tripped airfoil surface is also set to 1 × 10−5 m as well as the stretching ratio of 1.1. The origin is set at the airfoil leading edge with the X-axis pointing to the right and Y-axis upward. The circular outer boundary of the two O-type computational domains is set as the velocity-inlet condition and the airfoil surface as no-slip wall condition. The chord length is 0.762 m and the free-stream Reynolds number of air is set to 2.6 × 106 equal to a dynamic pressure of 30 psf in order to be consistent with parameters of theory of wing section [9]. It is stated here that the results of which airfoil is not stated explicitly in all below figures are about the original NACA 64-210 airfoil, including the experimental and the numerical ones.Table 1
Lift and drag coefficients for different sizes of the original airfoil mesh.
Mesh size
C
L AOA 0°
C
L AOA 10°
LWC 0 g/m3
25 g/m3
ΔCL
LWC 0 g/m3
25 g/m3
ΔCL
8
0
×
100
0.161
0.155
0.006
1.081
1.029
0.052
100
×
100
0.157
0.152
0.005
1.077
1.028
0.049
100
×
120
0.156
0.151
0.005
1.076
1.027
0.049
120
×
120
0.155
0.151
0.004
1.075
1.027
0.048
C
D AOA 0°
C
D AOA 10°
LWC 0 g/m3
25 g/m3
ΔCD
LWC 0 g/m3
25 g/m3
ΔCD
80
×
100
0.020
0.022
0.002
0.078
0.082
0.006
100
×
100
0.019
0.020
0.001
0.075
0.080
0.005
100
×
120
0.019
0.020
0.001
0.074
0.080
0.006
120
×
120
0.019
0.020
0.001
0.074
0.079
0.005
LWC: liquid water content; AOA: angle of attack;CL: lift coefficient; CD: drag coefficient.The incompressible air flow field is solved by FLUENT, a common commercial flow field solver, the details of which can be referred to in the help literature [21] and will not be repeated here. For a Reynolds number of 2.6 × 106, the flow characteristic is considered as turbulent, so turbulence model is added to solve the Navier-Stokes equations. During the calculation of the original and tripped airfoils in dry condition, the steady pressure-based solver is chosen, of which the segregated SIMPLE algorithm is adopted to discretize the pressure-velocity coupling term. The pressure term uses second-order scheme, and the QUICK scheme is used in the momentum term discretization.Airfoil aerodynamic performance is measured by lift and drag coefficients in this research, which are defined, respectively as follows:(1)CL=L(1/2)ρav∞2c,CD=D(1/2)ρav∞2c,
where CL is the lift coefficient and CD is the drag coefficient, L is the lift, D is the drag, ρa is the density of air, v∞ is the air free-stream velocity, and c is the chord length of the airfoil of interest.


## 2.2. Particulate Phase

### 2.2.1. Scaling of Rain Model

To study the heavy rain effect, first of all, it is necessary to measure the intensity and frequency of heavy rain. To study the heavy rain effect, first of all, it’s necessary to measure the intensity and frequency of heavy rain. Usually the rainfall rate,R in millimeter per hour or the Liquid Water Content, LWC in gram per cubic meter is chosen to categorize different intensities of rain. A rainfall of 100 mm/h or greater is often deemed as heavy. A rainfall of 100 mm/h or greater is often deemed as heavy.The LWC can be written as a function ofN(D) as follows:
(2)LWC=∫0+∞ρwπ6D3N(D)dD,
where ρw is the density of water. Integrating the above formula, we may attain the correlation of LWC and R by
(3)LWC=0.054R0.84.Subsequently, it is necessary to establish the size distribution of water droplets under different rain rates. Many researchers like Best [22], Ulbrich [23], and so on have established various raindrop size distribution formulas for various situations. Marshall and Palmer developed the classic formula of drop size distribution in 1948 based on massive experimental data [24]. It is shown as follows:
(4)N(D)=N0EXP(-ID)(0≤D≤Dmax),
where N(D) (m−3mm−1) is the number density of spherical raindrops of diameter D (mm) per cubic meter of air, Dmax is the maximum drop diameter, and N0 and I (mm−1) are parameters of N(D) and have different values for different types of rain. For storm-type heavy rainfall, I varies with rainfall rate R as I=3R-0.21, and N0 has the constant value: N0=1400 m−3mm−1 [25] correspondingly.Here, it is assumed that raindrops have been with uniform velocity before hitting the aircraft surface, that is, without acceleration. So it is important to determine the terminal velocity of raindrops. It has been developed by Markowitz [26] as follows:
(5)V(D)=9.58{1-EXP[-(D1.77)1.147]},
where V(D) is the terminal velocity. A correction for it aloft is given by Markowitz as
(6)V(D)=V0(D)(ρ0ρa)0.4,
where V0(D) is the terminal velocity consistent with the density of air aloft ρ0.


### 2.2.2. Wall-Film Model

In our study, the wall-film model in Fluent [21] is mainly adopted to model the interaction of particle and wall surface. It allows a single-component liquid drop to impinge upon a boundary surface of arbitrary configuration and form a thin liquid film. The major physical processes that affect the liquid film include mass and momentum contributions to the film thanks to drop impingement, droplet splashing effects, evaporation, shear forces on the film, dynamic pressure effects, gravity driven flow, convective heat and mass transfer, flow separation, and sheet breakup, as shown in Figure 2(a). In present study, we ignore the film evaporation to simplify our solution, so it is unnecessary to consider the effects of the thin liquid film on the air flow.(a) Mechanisms of mass, momentum, and energy transfers for the wall-film model; (b) the four impingement regimes included in the model; (c) simplified decision for wall interaction criterion.
(a)
(b)
(c)The main assumptions for the film model are as follows.(i)
The layer is thin, less than 500 microns in thickness due to the assumption of a linear velocity profile in the film.(ii)
The temperature in the film particles changes relatively slow due to the use of an analytical integration scheme.(iii)
The film temperature is always below the boiling temperature for the liquid.(iv)
Film particles are assumed to be in direct contact with the wall surface and the heat transfer from the wall to the film occurs through conduction.The wall interaction regimes are calculated for a drop-wall interaction based on local information. The four regimes including stick, rebound, spread, and splash are based on the impact energy and wall temperature, as shown in Figures2(b) and 2(c) (Tb is the liquid boiling temperature and Tw is the wall face temperature). Below the liquid boiling temperature, the impinging droplet can either stick, spread, or splash, while above the boiling temperature, the particle can either rebound or splash. As to our case of which the temperature is below the boiling point, particles stick, spread, and splash, resulting in aerodynamic efficiency degradation of an aircraft.The criteria by which the regimes are partitioned are based on the impact energy and the boiling temperature of the liquid. The impact energyE is defined by
(7)E2=ρwVr2Dσw(1min(h0/D,1)+δbl/D),
where ρw is the water density, Vr is the particle relative velocity in the frame of the wall (i.e., Vr=Vp-Vwall), σw is the water surface tension, and h0 is the total height of the wall film. δbl denotes the thickness of the boundary layer and is defined by
(8)δbl=μDρwVr.The sticking regime is applied when the dimensionless energyE is less than 16 and the particle velocity is set equal to the wall velocity. In the spreading regime, the probability of the drop having a particular direction along the surface is given by an analogy of an inviscid liquid jet with an empirically defined radial dependence for the momentum flux. If the wall temperature is above the boiling temperature of the liquid, impingement events below a critical impact energy (E) result in the particles rebounding from the wall. Splashing occurs when the impingement energy is above a critical energy threshold, defined as Ecr = 57.7. Besides, in our study, we sample a cumulative probability distribution function (CPDF), which is acquired from the Weibull distribution function and fitted to the data from Mundo et al. [27], to determine the different diameter of each splashed parcel. The equation of the cumulative probability distribution function can be expressed as
(9)pdf(diDp)=2diDp2exp[-(diDp)2]
and it represents the probability of finding drops of diameter di in a sample of splashed raindrops.Bilanin [28] has investigated the evaporation of the particles near the surface and found that evaporation does not affect the airfoil aerodynamic efficiency, so in our study we ignore the vaporization of water film. Since the film particle vaporization is ignored, only the momentum and energy conservation equations remain in the conservation equations for wall-film particles.Conservation equation of momentum is as follows:(10)ρahdV⃑pdt+h(∇spf)c=τgt⃑g+τwt⃑w+P⃑˙imp,c-M˙imp,cV⃑p+F⃑˙n,c+ρah(g⃑-a⃑w),
where c denotes the current face where the particle locates, h denotes the current film height where the particle locates, ∇s is the gradient operator as to the surface, pf is the film pressure, τg and τw are, respectively, the magnitude of the shear stress of the air flow on the film surface and the magnitude of the stress of the film exerted by the wall. t⃑g and t⃑w are the unit vectors in the direction of the relative motion of the film and the wall, P⃑˙imp,c denotes the impingement pressure on the film surface, M˙imp,c is the impingement momentum source, F⃑˙n,c is the force to keep the film on the surface, and ρah(g⃑-a⃑w) is the body force term in which a⃑w is the acceleration of the relative wall motion.Conservation equation of energy is as follows:(11)mpCpdTpdt=κAph(Tw-Tp)+hfAp(Ta-Tp),
where mp is the particle mass, Cp is the equilibrium concentration of the rain droplet, κ is the liquid thermal conductivity, h is the film height where the particle locates, hf is the film heat transfer coefficient, Ap is the area of the film particle, Tp, Tw, and Ta, respectively, denote the temperature of the particle, wall, and air. The first term on the right-hand side denotes the thermal conduction from the wall to the film particle, and the second term denotes the thermal convection from the top surface to the film particle.


### 2.2.3. Mathematic Modeling for Water Film

In this part, a numerical model for the water film based on [1] is adopted to calculate water film thickness by airfoil location and rainfall rate.  In each grid cell, the film thickness is resulted from water remaining after the gains and losses due to the incoming rain and film flow. After a time, an equilibrium thickness is completed in each cell. Here we will mainly examine the effect of water film on the airfoil geometric changes qualitatively, so for simplicity the pressure gradient and thermodynamic influences are neglected. It is assumed that the water film develops normally to the airfoil surface, as shown in Figure 3. The main program is acquired and shown as follows.Figure 3
Illustration of the water film development.The raindrop local collection efficiency within every grid cell, symbolized asβ(i), can be expressed as
(12)β(i)=n(i)n0Ymax-Ymins(i),
where n(i) is the number of raindrops within cell i, n0 is the total number of raindrops which can have an impact on the airfoil surface. Ymax and Ymin represent the maximum and minimum vertical position, respectively. s(i) is the arc length of cell i and  is shown as follows:
(13)s(i)=(x(i+1)-x(i))2+(y(i+1)-y(i))2,
where x(i) and y(i) represent position of cell i in X and Y direction, respectively. Consider
(14)β¯=∑1Nβ(i)N,
where β¯ is the average collection efficiency at the airfoil surface, N is the total number of cells adjacent to the airfoil surface.The air shear stress within every cell,τ(i) can be expressed as
(15)τ(i)=12Cf(i)×ρa×v¯2,
where Cf(i) is the air coefficient of friction and v¯ is the average velocity of cells at the airfoil surface as follows:
(16)v¯=∑1Nv(i)N,
where v(i) is the velocity of air in cell i.The resulting film is a balance between water runback and rain water reception. For the film in each cell, the equilibrium condition of the mass balance can be expressed as(17)dmidt=ρwA∫z=0x=iΔxz=hudz-ρwA∫z=0x=(i+1)Δxz=hudz+∂mi∂t,
where u=τz/μw is the film velocity at height z, μw is the molecular viscosity of water. The first and second terms at the right side represent, respectively, mass flux into the ith cell (x=iΔx) of Area A and out of the cell (x=(i+1)Δx) due to film flow. The third term represents mass flux into the ith cell due to rain.The change rate of raindrop mass in celli, (dm/dt)(i), is calculated from (12):
(18)dmdt(i)=12ρw×τ(i)[((N-(i+1))h0(i)N)2-((N-i)h0(i)N)2]·1μws(i)+LWC×β(i)×v∞,
where h0(i) is the initial thickness of water film in cell i and set as
(19)h0(i)={0,β(i)=0,0.0005,β(i)≠0.After all, a new water film thickness in celli, h(i), results from [1]:
(20)h(i)=h0(i)+1ρwdmdt(i)Δt,
where Δt is the time step.Experimental research shows that, once the water droplets impact the thin water film layer, scallops will form, and the water remaining in the scallops will form convex film layers, thus increasing the roughness of the original water film. Based on Macklin’s research [29], the present program assumes the convex film layers to be cylindrical so as to take the water film roughness into consideration.From Driling [30], the correlation equations for equivalent sandgrain roughness developed from experimental data can be expressed as the function of the average height of the convex film; that is
(21)ksk¯={0.0164Λ3.7,Λ≤4.93,139Λ-1.9,Λ≥4.93,
where ks is equivalent sandgrain roughness, k¯ is the average height of convex film, Λ is the correlating parameter, and both can be referred to [1].The original NACA 64-210 airfoil in wet condition of LWC = 25 g/m3 is simulated and the results are compared with the wind-tunnel experimental ones in [9]. The wind tunnel experiment was conducted by Bezos et al. in the NASA Langley 14- by 22-Foot Subsonic Tunnel in 1992 to determine the aerodynamic penalty associated with a simulated heavy rain encounter. The model was comprised of a NACA 64-210 airfoil section with a chord of 2.5 ft, a span of 8 ft and was mounted on the tunnel centerline between two large endplates. The rain simulation system manifold, which was located 10 chord lengths upstream of the model, produced liquid water contents ranging from 16 to 46 g/m3. Aerodynamic measurements in and out of the simulated rain environment were obtained to measure the airfoil aerodynamic coefficients. In our rain calculations, the pressure-based unsteady solver with the first-order implicit scheme is adopted; the calculation condition of rain is shown in Table 2.Table 2
Calculated rain condition.
Variable
Value
LWC/(g/m3)
25
Raindrop arithmetic mean diameter/(mm)
1
Air freestream velocity/(m/s)
48.4632
Raindrop vertical velocity/(m/s)
−4
Air density/(kg/m3)
1.225
Airfoil chord length/(m)
0.762
Air freestream Reynolds number
2.6
×
10
6
Reference pressure/(Pa)
1.01325
×
10
5
Reference temperature/(K)
273


## 2.2.1. Scaling of Rain Model

To study the heavy rain effect, first of all, it is necessary to measure the intensity and frequency of heavy rain. To study the heavy rain effect, first of all, it’s necessary to measure the intensity and frequency of heavy rain. Usually the rainfall rate,R in millimeter per hour or the Liquid Water Content, LWC in gram per cubic meter is chosen to categorize different intensities of rain. A rainfall of 100 mm/h or greater is often deemed as heavy. A rainfall of 100 mm/h or greater is often deemed as heavy.The LWC can be written as a function ofN(D) as follows:
(2)LWC=∫0+∞ρwπ6D3N(D)dD,
where ρw is the density of water. Integrating the above formula, we may attain the correlation of LWC and R by
(3)LWC=0.054R0.84.Subsequently, it is necessary to establish the size distribution of water droplets under different rain rates. Many researchers like Best [22], Ulbrich [23], and so on have established various raindrop size distribution formulas for various situations. Marshall and Palmer developed the classic formula of drop size distribution in 1948 based on massive experimental data [24]. It is shown as follows:
(4)N(D)=N0EXP(-ID)(0≤D≤Dmax),
where N(D) (m−3mm−1) is the number density of spherical raindrops of diameter D (mm) per cubic meter of air, Dmax is the maximum drop diameter, and N0 and I (mm−1) are parameters of N(D) and have different values for different types of rain. For storm-type heavy rainfall, I varies with rainfall rate R as I=3R-0.21, and N0 has the constant value: N0=1400 m−3mm−1 [25] correspondingly.Here, it is assumed that raindrops have been with uniform velocity before hitting the aircraft surface, that is, without acceleration. So it is important to determine the terminal velocity of raindrops. It has been developed by Markowitz [26] as follows:
(5)V(D)=9.58{1-EXP[-(D1.77)1.147]},
where V(D) is the terminal velocity. A correction for it aloft is given by Markowitz as
(6)V(D)=V0(D)(ρ0ρa)0.4,
where V0(D) is the terminal velocity consistent with the density of air aloft ρ0.


## 2.2.2. Wall-Film Model

In our study, the wall-film model in Fluent [21] is mainly adopted to model the interaction of particle and wall surface. It allows a single-component liquid drop to impinge upon a boundary surface of arbitrary configuration and form a thin liquid film. The major physical processes that affect the liquid film include mass and momentum contributions to the film thanks to drop impingement, droplet splashing effects, evaporation, shear forces on the film, dynamic pressure effects, gravity driven flow, convective heat and mass transfer, flow separation, and sheet breakup, as shown in Figure 2(a). In present study, we ignore the film evaporation to simplify our solution, so it is unnecessary to consider the effects of the thin liquid film on the air flow.(a) Mechanisms of mass, momentum, and energy transfers for the wall-film model; (b) the four impingement regimes included in the model; (c) simplified decision for wall interaction criterion.
(a)
(b)
(c)The main assumptions for the film model are as follows.(i)
The layer is thin, less than 500 microns in thickness due to the assumption of a linear velocity profile in the film.(ii)
The temperature in the film particles changes relatively slow due to the use of an analytical integration scheme.(iii)
The film temperature is always below the boiling temperature for the liquid.(iv)
Film particles are assumed to be in direct contact with the wall surface and the heat transfer from the wall to the film occurs through conduction.The wall interaction regimes are calculated for a drop-wall interaction based on local information. The four regimes including stick, rebound, spread, and splash are based on the impact energy and wall temperature, as shown in Figures2(b) and 2(c) (Tb is the liquid boiling temperature and Tw is the wall face temperature). Below the liquid boiling temperature, the impinging droplet can either stick, spread, or splash, while above the boiling temperature, the particle can either rebound or splash. As to our case of which the temperature is below the boiling point, particles stick, spread, and splash, resulting in aerodynamic efficiency degradation of an aircraft.The criteria by which the regimes are partitioned are based on the impact energy and the boiling temperature of the liquid. The impact energyE is defined by
(7)E2=ρwVr2Dσw(1min(h0/D,1)+δbl/D),
where ρw is the water density, Vr is the particle relative velocity in the frame of the wall (i.e., Vr=Vp-Vwall), σw is the water surface tension, and h0 is the total height of the wall film. δbl denotes the thickness of the boundary layer and is defined by
(8)δbl=μDρwVr.The sticking regime is applied when the dimensionless energyE is less than 16 and the particle velocity is set equal to the wall velocity. In the spreading regime, the probability of the drop having a particular direction along the surface is given by an analogy of an inviscid liquid jet with an empirically defined radial dependence for the momentum flux. If the wall temperature is above the boiling temperature of the liquid, impingement events below a critical impact energy (E) result in the particles rebounding from the wall. Splashing occurs when the impingement energy is above a critical energy threshold, defined as Ecr = 57.7. Besides, in our study, we sample a cumulative probability distribution function (CPDF), which is acquired from the Weibull distribution function and fitted to the data from Mundo et al. [27], to determine the different diameter of each splashed parcel. The equation of the cumulative probability distribution function can be expressed as
(9)pdf(diDp)=2diDp2exp[-(diDp)2]
and it represents the probability of finding drops of diameter di in a sample of splashed raindrops.Bilanin [28] has investigated the evaporation of the particles near the surface and found that evaporation does not affect the airfoil aerodynamic efficiency, so in our study we ignore the vaporization of water film. Since the film particle vaporization is ignored, only the momentum and energy conservation equations remain in the conservation equations for wall-film particles.Conservation equation of momentum is as follows:(10)ρahdV⃑pdt+h(∇spf)c=τgt⃑g+τwt⃑w+P⃑˙imp,c-M˙imp,cV⃑p+F⃑˙n,c+ρah(g⃑-a⃑w),
where c denotes the current face where the particle locates, h denotes the current film height where the particle locates, ∇s is the gradient operator as to the surface, pf is the film pressure, τg and τw are, respectively, the magnitude of the shear stress of the air flow on the film surface and the magnitude of the stress of the film exerted by the wall. t⃑g and t⃑w are the unit vectors in the direction of the relative motion of the film and the wall, P⃑˙imp,c denotes the impingement pressure on the film surface, M˙imp,c is the impingement momentum source, F⃑˙n,c is the force to keep the film on the surface, and ρah(g⃑-a⃑w) is the body force term in which a⃑w is the acceleration of the relative wall motion.Conservation equation of energy is as follows:(11)mpCpdTpdt=κAph(Tw-Tp)+hfAp(Ta-Tp),
where mp is the particle mass, Cp is the equilibrium concentration of the rain droplet, κ is the liquid thermal conductivity, h is the film height where the particle locates, hf is the film heat transfer coefficient, Ap is the area of the film particle, Tp, Tw, and Ta, respectively, denote the temperature of the particle, wall, and air. The first term on the right-hand side denotes the thermal conduction from the wall to the film particle, and the second term denotes the thermal convection from the top surface to the film particle.


## 2.2.3. Mathematic Modeling for Water Film

In this part, a numerical model for the water film based on [1] is adopted to calculate water film thickness by airfoil location and rainfall rate.  In each grid cell, the film thickness is resulted from water remaining after the gains and losses due to the incoming rain and film flow. After a time, an equilibrium thickness is completed in each cell. Here we will mainly examine the effect of water film on the airfoil geometric changes qualitatively, so for simplicity the pressure gradient and thermodynamic influences are neglected. It is assumed that the water film develops normally to the airfoil surface, as shown in Figure 3. The main program is acquired and shown as follows.Figure 3
Illustration of the water film development.The raindrop local collection efficiency within every grid cell, symbolized asβ(i), can be expressed as
(12)β(i)=n(i)n0Ymax-Ymins(i),
where n(i) is the number of raindrops within cell i, n0 is the total number of raindrops which can have an impact on the airfoil surface. Ymax and Ymin represent the maximum and minimum vertical position, respectively. s(i) is the arc length of cell i and  is shown as follows:
(13)s(i)=(x(i+1)-x(i))2+(y(i+1)-y(i))2,
where x(i) and y(i) represent position of cell i in X and Y direction, respectively. Consider
(14)β¯=∑1Nβ(i)N,
where β¯ is the average collection efficiency at the airfoil surface, N is the total number of cells adjacent to the airfoil surface.The air shear stress within every cell,τ(i) can be expressed as
(15)τ(i)=12Cf(i)×ρa×v¯2,
where Cf(i) is the air coefficient of friction and v¯ is the average velocity of cells at the airfoil surface as follows:
(16)v¯=∑1Nv(i)N,
where v(i) is the velocity of air in cell i.The resulting film is a balance between water runback and rain water reception. For the film in each cell, the equilibrium condition of the mass balance can be expressed as(17)dmidt=ρwA∫z=0x=iΔxz=hudz-ρwA∫z=0x=(i+1)Δxz=hudz+∂mi∂t,
where u=τz/μw is the film velocity at height z, μw is the molecular viscosity of water. The first and second terms at the right side represent, respectively, mass flux into the ith cell (x=iΔx) of Area A and out of the cell (x=(i+1)Δx) due to film flow. The third term represents mass flux into the ith cell due to rain.The change rate of raindrop mass in celli, (dm/dt)(i), is calculated from (12):
(18)dmdt(i)=12ρw×τ(i)[((N-(i+1))h0(i)N)2-((N-i)h0(i)N)2]·1μws(i)+LWC×β(i)×v∞,
where h0(i) is the initial thickness of water film in cell i and set as
(19)h0(i)={0,β(i)=0,0.0005,β(i)≠0.After all, a new water film thickness in celli, h(i), results from [1]:
(20)h(i)=h0(i)+1ρwdmdt(i)Δt,
where Δt is the time step.Experimental research shows that, once the water droplets impact the thin water film layer, scallops will form, and the water remaining in the scallops will form convex film layers, thus increasing the roughness of the original water film. Based on Macklin’s research [29], the present program assumes the convex film layers to be cylindrical so as to take the water film roughness into consideration.From Driling [30], the correlation equations for equivalent sandgrain roughness developed from experimental data can be expressed as the function of the average height of the convex film; that is
(21)ksk¯={0.0164Λ3.7,Λ≤4.93,139Λ-1.9,Λ≥4.93,
where ks is equivalent sandgrain roughness, k¯ is the average height of convex film, Λ is the correlating parameter, and both can be referred to [1].The original NACA 64-210 airfoil in wet condition of LWC = 25 g/m3 is simulated and the results are compared with the wind-tunnel experimental ones in [9]. The wind tunnel experiment was conducted by Bezos et al. in the NASA Langley 14- by 22-Foot Subsonic Tunnel in 1992 to determine the aerodynamic penalty associated with a simulated heavy rain encounter. The model was comprised of a NACA 64-210 airfoil section with a chord of 2.5 ft, a span of 8 ft and was mounted on the tunnel centerline between two large endplates. The rain simulation system manifold, which was located 10 chord lengths upstream of the model, produced liquid water contents ranging from 16 to 46 g/m3. Aerodynamic measurements in and out of the simulated rain environment were obtained to measure the airfoil aerodynamic coefficients. In our rain calculations, the pressure-based unsteady solver with the first-order implicit scheme is adopted; the calculation condition of rain is shown in Table 2.Table 2
Calculated rain condition.
Variable
Value
LWC/(g/m3)
25
Raindrop arithmetic mean diameter/(mm)
1
Air freestream velocity/(m/s)
48.4632
Raindrop vertical velocity/(m/s)
−4
Air density/(kg/m3)
1.225
Airfoil chord length/(m)
0.762
Air freestream Reynolds number
2.6
×
10
6
Reference pressure/(Pa)
1.01325
×
10
5
Reference temperature/(K)
273


## 3. Results and Discussions

### 3.1. Validation

To get a better numerical accuracy, two turbulence models, namely the Spalart-Allmaras model (abbreviated as SA model) and K-E model, are chosen to get the lift and drag coefficients of the airfoil over a wider range of attack angles in the dry condition, and the results of both are compared with the experimental data from wing section theory [9], as is shown in Figures 4 and 5. It can be seen that both the turbulence models have the same accuracy of lift coefficient at low angles of attack, while the SA model has a better tendency at high angles of attack. What is more, the SA model has fairly satisfactory results of drag coefficient at low angles of attack. Hence, the SA model is chosen as the major turbulence model, and all the results below are achieved by it.Figure 4
Lift coefficient comparison between numerical and experimental results.Figure 5
Drag coefficient comparison between numerical and experimental results.


### 3.2. Lift and Drag Force Data

Figures6 and 7 present the lift and drag coefficient polar data in dry and wet conditions for numerical and wind tunnel experimental approach, respectively. It shows a decrease in lift and an increase in drag due to the heavy rain effect for both approaches. For our simulation, the maximum percentage decrease in CL is reached by 13.2% and the maximum percentage increase in CD by 47.6% for the airfoil in the rain condition. And for the experiment, the maximum percentage decrease in CL is reached by 13.7% and the maximum percentage increase in CD by 56.3%. Our simulation shows good agreement with the experiment. Although the numerical lift data in both conditions are larger than the experimental data, the trend matches well. The numerical drag coefficients correspond well with the wind tunnel experimental data before AOA = 13 degree but smaller than the experimental data at higher angles of attack.Figure 6
Lift coefficient versus AOA for numerical and experimental results.Figure 7
Drag coefficient versus AOA for numerical and experimental results.Figures8 and 9 show the pressure contour and the streamlines around the NACA 64-210 airfoil at AOA = 14 degree in the dry and wet conditions, respectively. The operating pressure is 101325 Pa. In the dry condition, the boundary layer has no rain-induced changes at the moment. While in the wet condition, the boundary layer has separated due to rain. The two contours reveal that premature boundary-layer separation may happen for airfoil in rain at high angles of attack; this conclusion is consistent with that in [5]. Besides, the difference of the pressure on the upper and lower surfaces of the airfoil leading edge decreases with rain; thus the lift decreases in the rain.Figure 8
Pressure distribution and streamlines around the airfoil in the dry condition.Figure 9
Pressure distribution and streamlines around the airfoil in the rain condition.


### 3.3. Forced Boundary-Layer Transition Results

To prove our hypothesis of premature boundary layer transition at low angles of attack, the NACA 64-210 airfoil with tripped boundary layer is emulated in the dry condition. Figure10 shows the numerical lift coefficient comparison for the original airfoil in dry and wet (LWC 25 g/m3) conditions and the tripped airfoil with different tripping positions in the dry condition. Evidently, trip strips at 5% chord on the top surface best modeled the wet conditions. However, the high-AOA behavior of the NACA 64-210 airfoil is not emulated by the tripping technique, possibly because the tripped position is aft of the leading-edge separation point. The overall ability to model aerodynamic performance degradation in heavy rain condition with our numerically tripped boundary-layer technique suggests that the aerodynamic degradation results from premature boundary-layer transition at low angles of attack, which is also consistent with that in [5]. Figures 11 and 12 show the component of air velocity in the X direction, namely X-velocity distribution near the leading edge of the original airfoil and the 5% chordwise tripped airfoil at the same AOA of 8 degree in the dry condition. It can be clearly seen that the negative X-velocity appears at the tripped location compared with the smooth original one, which manifests that the premature boundary-layer transition is gained.Figure 10
Numerical lift coefficient comparison for different tripped positions.Figure 11
X-velocity distribution of air around the original airfoil.Figure 12
X-velocity distribution of air around the 5% chordwise tripped airfoil.


### 3.4. Tracking the Particulate Trajectories

From Figures13, 14, 15, 16, 17, and 18 the raindrop traces colored by particulate X-velocity at AOA = 0, 4, 8, 12, 16, and 20 degrees at time of t=0.125 s are shown. At AOA = 0 degree, raindrops concentrate on both the top and the bottom surfaces of the airfoil. At AOA = 4 degree, there have been few raindrops on the top surface but still many flowing from the stagnation point to downstream on the bottom. At AOA = 8 degree, the discrete phase has developed a full separation on the bottom surface, so there are almost no raindrops. Meanwhile, the discrete phase begins to separate on the top surface. As the angle of attack increases further, the separation extent becomes larger coupling with the presence of the continuous phase separation. At AOA = 20 degree, a maximum degree of coupled separation is reached and large-scale vortexes appear resulting in the runback of rivulets on the upper surface. The patterns should be viewed qualitatively; however, our grid is nonuniform and resolution decreases away from the airfoil surface. Turbulent dispersion of the small splashed-back droplets is not modeled, and thus the contours show mean trajectories.Figure 13
Raindrop traces at AOA = 0 degree.Figure 14
Raindrop traces at AOA = 4 degree.Figure 15
Raindrop traces at AOA = 8 degree.Figure 16
Raindrop traces at AOA = 12 degree.Figure 17
Raindrop traces at AOA = 16 degree.Figure 18
Raindrop traces at AOA = 20 degree.Figure19 shows the droplet concentration within the second cell layer from the airfoil surface at AOA of 8 degrees. Since the freestream raindrop number density is of the order of 103 raindrops/m3 and is well below the minimum value of 1.21 × 105 of the droplet concentration, it can be assumed that the vast majority of droplets in the second cell layer is owing to splashback. For an angle of attack of 8 degrees, the stagnation point is slightly below the leading edge, so more splashed-back raindrops from impacts at the leading edge are carried over the upper surface than the lower surface. What is more, the downward component of raindrop velocity leads to more impacts on the upper surface, especially for larger raindrops which have higher vertical velocities.Figure 19
Droplet concentration within the second cell layer from the airfoil surface at AOA = 8 degree.


### 3.5. Geometric Changes Based on Water Film Layer

Figure20 shows the leading-edge water film thickness comparison between LWC = 25 g/m3, 39 g/m3, and the dry condition at AOA of 14 degrees and time of 0.125 s. It can be concluded that as the LWC increases, the water film thickness increases simultaneously at the same angle of attack, which conforms to the objective fact. For the condition of LWC = 25 g/m3, the resulting average film thickness for the top and bottom surfaces of a 1-meter-chord-length airfoil is, respectively, 0.3502 and 0.3167 mm. And for the condition of LWC = 39 g/m3, it increases up to 0.3794 and 0.3469 mm, respectively. As analyzed in Subsection 4 of this part, for an angle of attack of 14 degrees more raindrops have an impact on the top surface and certainly form a thicker film. Besides, the film thickness increases with the increasing rain rate (or liquid water content). The water film increases the weight, changes the geometry of surface, and may cause vital aerodynamic performance degradation of aircraft.Figure 20
Water film thickness comparison between different LWC conditions at AOA = 14 degree.Figure21 shows the leading-edge water film thickness comparison for AOA = 0 and 14 degrees and LWC = 25 g/m3 at the same time. From the figure we can see that, in the same LWC conditions, water films have approximately the same thickness at the stagnation point at different angles of attack. But at higher angles of attack, water films on the bottom surface get thicker, this is because that at higher angles of attack, water collection efficiency gets larger on the bottom surface and water on both surfaces runs faster to downstream.Figure 21
Water film thickness comparison at different angles of attack and LWC = 25 g/m3 condition.


## 3.1. Validation

To get a better numerical accuracy, two turbulence models, namely the Spalart-Allmaras model (abbreviated as SA model) and K-E model, are chosen to get the lift and drag coefficients of the airfoil over a wider range of attack angles in the dry condition, and the results of both are compared with the experimental data from wing section theory [9], as is shown in Figures 4 and 5. It can be seen that both the turbulence models have the same accuracy of lift coefficient at low angles of attack, while the SA model has a better tendency at high angles of attack. What is more, the SA model has fairly satisfactory results of drag coefficient at low angles of attack. Hence, the SA model is chosen as the major turbulence model, and all the results below are achieved by it.Figure 4
Lift coefficient comparison between numerical and experimental results.Figure 5
Drag coefficient comparison between numerical and experimental results.


## 3.2. Lift and Drag Force Data

Figures6 and 7 present the lift and drag coefficient polar data in dry and wet conditions for numerical and wind tunnel experimental approach, respectively. It shows a decrease in lift and an increase in drag due to the heavy rain effect for both approaches. For our simulation, the maximum percentage decrease in CL is reached by 13.2% and the maximum percentage increase in CD by 47.6% for the airfoil in the rain condition. And for the experiment, the maximum percentage decrease in CL is reached by 13.7% and the maximum percentage increase in CD by 56.3%. Our simulation shows good agreement with the experiment. Although the numerical lift data in both conditions are larger than the experimental data, the trend matches well. The numerical drag coefficients correspond well with the wind tunnel experimental data before AOA = 13 degree but smaller than the experimental data at higher angles of attack.Figure 6
Lift coefficient versus AOA for numerical and experimental results.Figure 7
Drag coefficient versus AOA for numerical and experimental results.Figures8 and 9 show the pressure contour and the streamlines around the NACA 64-210 airfoil at AOA = 14 degree in the dry and wet conditions, respectively. The operating pressure is 101325 Pa. In the dry condition, the boundary layer has no rain-induced changes at the moment. While in the wet condition, the boundary layer has separated due to rain. The two contours reveal that premature boundary-layer separation may happen for airfoil in rain at high angles of attack; this conclusion is consistent with that in [5]. Besides, the difference of the pressure on the upper and lower surfaces of the airfoil leading edge decreases with rain; thus the lift decreases in the rain.Figure 8
Pressure distribution and streamlines around the airfoil in the dry condition.Figure 9
Pressure distribution and streamlines around the airfoil in the rain condition.


## 3.3. Forced Boundary-Layer Transition Results

To prove our hypothesis of premature boundary layer transition at low angles of attack, the NACA 64-210 airfoil with tripped boundary layer is emulated in the dry condition. Figure10 shows the numerical lift coefficient comparison for the original airfoil in dry and wet (LWC 25 g/m3) conditions and the tripped airfoil with different tripping positions in the dry condition. Evidently, trip strips at 5% chord on the top surface best modeled the wet conditions. However, the high-AOA behavior of the NACA 64-210 airfoil is not emulated by the tripping technique, possibly because the tripped position is aft of the leading-edge separation point. The overall ability to model aerodynamic performance degradation in heavy rain condition with our numerically tripped boundary-layer technique suggests that the aerodynamic degradation results from premature boundary-layer transition at low angles of attack, which is also consistent with that in [5]. Figures 11 and 12 show the component of air velocity in the X direction, namely X-velocity distribution near the leading edge of the original airfoil and the 5% chordwise tripped airfoil at the same AOA of 8 degree in the dry condition. It can be clearly seen that the negative X-velocity appears at the tripped location compared with the smooth original one, which manifests that the premature boundary-layer transition is gained.Figure 10
Numerical lift coefficient comparison for different tripped positions.Figure 11
X-velocity distribution of air around the original airfoil.Figure 12
X-velocity distribution of air around the 5% chordwise tripped airfoil.


## 3.4. Tracking the Particulate Trajectories

From Figures13, 14, 15, 16, 17, and 18 the raindrop traces colored by particulate X-velocity at AOA = 0, 4, 8, 12, 16, and 20 degrees at time of t=0.125 s are shown. At AOA = 0 degree, raindrops concentrate on both the top and the bottom surfaces of the airfoil. At AOA = 4 degree, there have been few raindrops on the top surface but still many flowing from the stagnation point to downstream on the bottom. At AOA = 8 degree, the discrete phase has developed a full separation on the bottom surface, so there are almost no raindrops. Meanwhile, the discrete phase begins to separate on the top surface. As the angle of attack increases further, the separation extent becomes larger coupling with the presence of the continuous phase separation. At AOA = 20 degree, a maximum degree of coupled separation is reached and large-scale vortexes appear resulting in the runback of rivulets on the upper surface. The patterns should be viewed qualitatively; however, our grid is nonuniform and resolution decreases away from the airfoil surface. Turbulent dispersion of the small splashed-back droplets is not modeled, and thus the contours show mean trajectories.Figure 13
Raindrop traces at AOA = 0 degree.Figure 14
Raindrop traces at AOA = 4 degree.Figure 15
Raindrop traces at AOA = 8 degree.Figure 16
Raindrop traces at AOA = 12 degree.Figure 17
Raindrop traces at AOA = 16 degree.Figure 18
Raindrop traces at AOA = 20 degree.Figure19 shows the droplet concentration within the second cell layer from the airfoil surface at AOA of 8 degrees. Since the freestream raindrop number density is of the order of 103 raindrops/m3 and is well below the minimum value of 1.21 × 105 of the droplet concentration, it can be assumed that the vast majority of droplets in the second cell layer is owing to splashback. For an angle of attack of 8 degrees, the stagnation point is slightly below the leading edge, so more splashed-back raindrops from impacts at the leading edge are carried over the upper surface than the lower surface. What is more, the downward component of raindrop velocity leads to more impacts on the upper surface, especially for larger raindrops which have higher vertical velocities.Figure 19
Droplet concentration within the second cell layer from the airfoil surface at AOA = 8 degree.


## 3.5. Geometric Changes Based on Water Film Layer

Figure20 shows the leading-edge water film thickness comparison between LWC = 25 g/m3, 39 g/m3, and the dry condition at AOA of 14 degrees and time of 0.125 s. It can be concluded that as the LWC increases, the water film thickness increases simultaneously at the same angle of attack, which conforms to the objective fact. For the condition of LWC = 25 g/m3, the resulting average film thickness for the top and bottom surfaces of a 1-meter-chord-length airfoil is, respectively, 0.3502 and 0.3167 mm. And for the condition of LWC = 39 g/m3, it increases up to 0.3794 and 0.3469 mm, respectively. As analyzed in Subsection 4 of this part, for an angle of attack of 14 degrees more raindrops have an impact on the top surface and certainly form a thicker film. Besides, the film thickness increases with the increasing rain rate (or liquid water content). The water film increases the weight, changes the geometry of surface, and may cause vital aerodynamic performance degradation of aircraft.Figure 20
Water film thickness comparison between different LWC conditions at AOA = 14 degree.Figure21 shows the leading-edge water film thickness comparison for AOA = 0 and 14 degrees and LWC = 25 g/m3 at the same time. From the figure we can see that, in the same LWC conditions, water films have approximately the same thickness at the stagnation point at different angles of attack. But at higher angles of attack, water films on the bottom surface get thicker, this is because that at higher angles of attack, water collection efficiency gets larger on the bottom surface and water on both surfaces runs faster to downstream.Figure 21
Water film thickness comparison at different angles of attack and LWC = 25 g/m3 condition.


## 4. Conclusion

In summary, numerical experiments at a Reynolds Number of 2.6 × 106 and LWC of 25 g/m3 have been conducted in the dry and wet conditions to compare the quantitative and qualitative aerodynamic penalties of NACA 64-210 airfoil with wind-tunnel experiments in heavy rain conditions. A number of mechanisms underlying are found by numerical simulation approach. The innovative points of present study are fivefold as follows.(i)
Aerodynamic penalties of NACA 64-210 airfoil over a wider range of AOA up to 20 degrees are determined. Airfoil aerodynamic performance in rain conditions at high AOA more than 15 degree has little been studied by numerical simulation before. The numerical simulation shows a shift of stall angle from 15 degrees in the dry condition to 13 degrees in the wet condition. For our simulation, the maximum percentage decrease inCL is reached by 13.2% and the maximum percentage increase in CD by 47.6%, for the airfoil in the rain condition.  The numerical drag coefficients correspond well with the wind tunnel experimental data before AOA = 13 degrees.(ii)
A newly creative boundary-layer-tripped technique is developed to prove that the airfoil aerodynamic degradation results from premature boundary-layer transition in heavy rain conditions at low angles of attack, which is consistent with the corresponding reference.(iii)
The pressure distribution contour and the streamlines around the airfoil are attained by using flow-visualization technique so as to illustrate the hypothesis that premature boundary-layer separation may happen for airfoil in heavy rain conditions at high angles of attack, which is also consistent with the corresponding reference.(iv)
The trajectories of raindrops at various angles of attack are tracked by flow-visualization technique, and the corresponding flow characteristics are discussed in detail. Impact and splash of raindrops, flow and runback of rivulets, and continuous phase and discrete phase coupled separation are displayed. Besides, a detailed analysis of raindrop splashback is given.(v)
A mathematic water-film model is established and the corresponding program is developed to get the water film on the airfoil surface, and the effects of LWC and AOA on water film thickness are discussed. The water film alters the airfoil geometry effectively, but this effect is most likely exaggerated in our study due to the small scale.It is recognized that because of the small-scale test case and specified flight conditions, the results of our tests may not be extrapolated to larger scales cases and various real flight conditions. However, the physical mechanisms of rain effects and their importance may be similar. So it is believed that all efforts in this study can be a preliminary evaluation of potential heavy rain mechanisms induced at larger scales and an important reference for civil aviation community.



---

*Source: 590924-2013-11-07.xml*