1. Introduction

Recent research in unmanned aerial vehicles (UAVs) deals with experimental measurements in which the wings present different aerodynamic characteristics as the chord Reynolds number or the aspect ratio are varied. It is also of great interest the aerodynamic performance regarding the maximum lift, lift curve slope, and polar curves [1]. Typically, UAVs operate in the range between 50,000 and 150,000 chord Reynolds number. There is also a great interest to study 3D shape-change and optimization frameworks where the aspect ratio plays a significant role [2-3]. This is the motivation of our experimental study: to analyse the aerodynamic characteristics at low-to-moderate Reynolds numbers in order to compare our results with other works with different aspect ratios.

A wing profile is a surface that might be designed to provide lift force and the minimum drag. The relationship between both forces is determined by the wing cross section aerodynamic features [4]. In finite wings at low Reynolds numbers, drag and lift coefficient variations are mainly due to three mechanisms: Wingtip vortex [5-9]; laminar boundary layer separation leading to the formation of a laminar separation bubble (LSB) and the subsequent turbulent separated shear layer [10-12] and finally the vortex shedding in the wake behind the wing [13-15]. Most of these numerical or experimental investigations were performed analysing specifically the flow behaviour. Hence, C_D and C_L coefficients involve a whole fluid-structure scenario, giving us an overview of the interaction between the wing and the flow that passes over it.

The NACA0012 airfoil has been extensively studied and its aerodynamic features are well known as well as the comparison between numerical and experiments by using Large Eddy Simulation [16]. Abbot and von Doenhoff [4] presented a large experimental data summary for different 2D airfoils. The maximum C_L observed by these authors was 1.1 - 1.6 for chord based Reynolds numbers ranging from 3·10⁶ to 9·10⁶, and with stall angles between 12º and 16º. For lower angles than stall, the lift coefficient increases linearly with a constant slope ΔC_L/Δα. The minimum C_D depends on the roughness of the wing surface. Two decades later, Sheldahl and Klimas [17] tested a 2D NACA0012 model for chord based Reynolds numbers between 3.6·10⁵ ≤ Re_c ≤ 7·10⁵ and AoA in the range 0º < α< 180º. Their results showed a reduction of the maximum lift coefficient (C_Lmax = 0.9 - 1.0) and the stall angle (α_stall = 10º - 12º) in comparison to the study of Abbot and von Doenhoff [4], though the slope of the lift curves as function of α for small AoA were in agreement. It was also observed by Sheldahl and Klimas [17] that the slope decreases for the AoA greater than 5º - 6º. Nevertheless, if the aspect ratio (AR) between the chord and the span of the wing decreases up to a certain limit value, three-dimensional effects on the wingtip dynamics play a key role. This modifies the expected aerodynamics performance of the wing and there is a meandering movement of the vortex core that is complex to identify [18, 19], though low frequencies have been measured in this random movement for the same NACA0012 airfoil studied here by means of flow visualizations and Proper Orthogonal Decomposition [20]. Besides, the vortex core location has been adjusted while doing PIV measurements due to the presence of the vortex wandering [21]. Actually, wingtip vortex on finite wings generates a lift coefficient reduction respect to the 2D infinite airfoil (called “downwash”), as well as the increment of the drag coefficient. This behaviour is even more significant if the aspect ratio reduces its value. Laitone [22] performed tests with a rectangular finite wing with a NACA0012 cross section profile, with AR=6, chord based Reynolds numbers below 7·10⁴, and a free stream turbulence intensity 0.02% - 0.1%. These results showed a discontinuity at small AoA in the slope on the linear region of the lift curves for Re_c = 2.07·10⁴ and Re_c = 4.21·10⁴. This discontinuity decreases as the free stream velocity increases together with the smooth increment of the slope ΔC_L\Δα. Likewise, the maximum lift coefficient for a constant Reynolds number decreases as the free stream turbulence intensity increases. Laitone also reported a slight reduction of the minimum drag coefficient C__Dmin (null AoA) as Re_c increases. The C_Dmin evolution with Reynolds number has a curve fitting C__Dmin=0.35·Re_c^-0.25. Moreover, Nngo and Barlow [23] conducted measurements in a low turbulence wind tunnel of a rectangular NACA0012 (AR≈4) in order to study the mechanism that reduces the drag force induced by the wingtip vortex for Re_c = 4.8·10⁵. They obtained a value of C_Dmin = 0.048 for α = 0º, while C_Lmax is close 0.6 for the critical (stall) angle α_stall = 11º. These values are reviewed to highlight how the aerodynamic characteristics of the NACA0012 airfoil strongly depends on the geometry (e.g. aspect ratio) and the presence of edge effects.

Finally, Mueller [5] measured experimentally drag and lift coefficients for different wing models in a wind tunnel. These models had semi-span aspect ratios (sAR) between 0.5 and 3 and they were machined using thin flat plates with rounded edges. The Reynolds numbers tested were in the range between 6·10⁴ and 2·10⁵, and AoA between -15º and +25º. Mueller found a direct relationship between the wing aspect ratio and the drag and lift polar profiles. Specifically, for a rectangular flat plate model and decreasing the sAR from 3 to 0.5, the lift coefficient is reduced by 50% in comparison to the rounded edge model. However, as the semi-aspect ratio is decreased, the linear region of C_L is extended to greater values of α, as well as an increment of the stall angle. In other work [6], the same author reported a value of C_Lmax = 0.8 with a stall angle of 18º for the (flat plate) rectangular wing model with aspect ratio AR = 2, together with other experimental observation: an increase of the drag force due to the occurrence of edge effects, not moving from squared to rounded edges but due to wingtip vortex formation.

The outline of this manuscript is as follows. Details of the experimental setup are given in Section 2 including a brief description on aerodynamic coefficient computations. Results and discussions are shown in order to achieve a better understanding of the general curve fittings for several aspect ratios in Section 3, including an estimation of sources of experimental errors in different flow regimes. Finally, we draw the main findings in Section 4.

2. Experimental Arrangement

Experimental tests are performed in the Vehicle Aero-Hydrodynamics Laboratory (VAHL) at Málaga University, with a closed (return-flow) low-speed wind tunnel, which has a 4 m long closed test section, with one squared meter cross-section. The free-stream velocity ranging from U_∞ = 4 to 30 m/s. The turbulence level (Turbulence Intensity I [%]) is shown in Table 1. The calibration of the free-stream velocity was done by means of Laser Doppler Anemometry (LDA). The percentage of power %P supplied to the four wind tunnel fans (relative to their nominal power) is the control parameter. This parameter produces a wind velocity U_∞ with a given deviation δU_{∞ LDA}. The values %P, U_{∞ LDA} and δU_{∞ LDA} are also shown in Table 1. We use four different free-stream velocities. The chord based Reynolds number Re_c is defined as Re_c = U_∞·c/ν; ν being the temperature-dependent kinematic viscosity of the air. The blockage ratio is lower than 1\%, so this effect is negligible. The model is a solid rectangular wing NACA0012 airfoil. The wing model has a chord and a length of c=100 mm and l=200 mm, respectively, so the aspect ratio is AR = l/c = 2. The NACA0012 airfoil is made from an aluminium alloy with rounded wingtip. The maximum thickness is 12 mm and it is located at 30 mm from the profile leading edge [20, 21].

The wing is firmly fixed to a cylindrical base, also made of aluminium, allowing the coupling of the wing to a precision force sensor placed under the test section floor (see the sketch in Figure 1 and photographs in Figure 2). This digital device measures 3D forces and it is levelled taking into account the location of the Z-axis in the vertical position respect to the measurement section floor, in other words, Z-axis is parallel to the gravitational force. The X-Y force sensor plane is perpendicular to the longitudinal axis of the wing, the X-axis being parallel to U_∞. The force sensor is coupled to an automatic rotation system which allows to vary the wing-base system orientation in the range of α = [-180º, 180º]. The rotation mechanism was previously calibrated to ensure its accuracy [24].

Experiments are carried out during T = 200 s, and for the AoA ranging from α = 0º to α = +35º constraining the rotation angle to 1 or 2 degrees increments, depending on the case.

Figure 1. Schematics of experimental set-up with the wing model mounted

Table 1. Percentage of power, wind velocities and their variations, chord based Reynolds numbers and levels of turbulence intensity

We record forces in each orthogonal direction in the range 0 V ≤ V_output ≤ 5 V, with a sampling frequency of f_s = 250 points/s. The conversion factor between the force and the electrical output is k = 32N/5V. The rotation system is also calibrated in order to find the null angle where the lift force is equal to zero and the drag force presents minimum values.

This is done to find the angular position in a relative coordinate system (respect to the X-axis of the force sensor) at which the airfoil is aligned to the free-stream direction. Hence, a set of tests are conducted for the AoA between -10º ≤ α ≤ +10º (not shown) in order to detect the symmetry axis of the as function of α, that corresponds to the angle where the force F presents its minimum value. Besides, and following an accurate experimental procedure [24], a set of experiments for different AoA are carried out without velocity in the tunnel in order to determine the force offsets (F_X0,F_Y0) in the plane X-Y, and possible deviations in the Z axis measuring the weight of the wing. Therefore, the expressions for the different angles and forces are:

(1)

(2)

Where α is the AoA, α_dev is the angle between the X force sensor axis and the null AoA, and α_axis is the angle between the X force sensor axis and the drag direction. F_X, F_Y correspond to raw data from the force sensor in the axis X and Y, respectively, and F_Xnet, F_Ynet are net forces in X and Y directions. Net forces are post-processed using Matlab^® to obtain the lift and drag forces as follows:

(3)

(4)

D, L being the drag and lift forces on the wing model, which are used to compute the drag and lift non-dimensional coefficients. We calculate these coefficients taking into account the air density as function of temperature ρ(T). Tests are conducted continuously, acquiring the air temperature inside the tunnel for each experiment. Fluctuations in the air temperature are within 0.2 K for a single test. The equations to compute drag and lift coefficients are

(5)

(6)

ρ being the air density inside the wind tunnel, and A = l·c the aerodynamic area of the rectangular model (200 · 100 mm²). We repeat three times each experiment to obtain average values. The relaxation time between measurements at each AoA was large enough to avoid any kind of perturbation. Besides, the unsteady force was monitored at each AoA to ensure that the force reached a settled value.

A photograph of the experimental set up real is shown in Figure 2. The photograph above has been obtained from inside the tunnel, and shows the NACA 0012. The photograph below shows the digital force sensor coupled to an automatic rotation system which allows to vary the wing-base system orientation.

Figure 2. Up: NACA 0012 profile inside the wind tunnel; Down: digital force sensor coupled to an automatic rotation system which allows to vary the wing-base system orientation in the range of α = [-180º, 180º]

3. Aerodynamic Characteristics: Results and Discussion

Drag coefficients as function of α are depicted in Figure 3. One can observe in the inset that our results also show a slight deviation in comparison to those obtained by Ngo [23]. Besides, our values of C_Dmin are in agreement with those reported by Mueller [6] for a flat wing with the same aspect ratio and Reynolds numbers. However, our results differ from those given by Sheldahl [17] for a bidimensional NACA0012 airfoil. The main cause of this discrepancy is the appearance of effects in the trailing edge, e.g. wingtip vortex formation, which dominate the wing aerodynamic performance due to the low-aspect-ratio of the model, AR = 2. In addition, the tip vortex constrains the laminar separation bubble at the leading edge. This effect was reported in Yen [25]. For this reason, the drag coefficient induced by three-dimensional effects at null AoA is approximately five times greater than the value of C_Dmin for a 2D profile measured by Sheldahl [17].

Regarding the ratio between drag coefficient and AoA, ΔC_D/Δα, we observe an increment of its value up to the stall angle, so that the wing loses its aerodynamic effectiveness. Different slopes are found at each Reynolds number, increasing the values of the slope with the AoA up to the stall angle. These ratios are of great interest in the case of the lift coefficient, as it will be discussed below.

Figure 3. CD vs α for all Reynolds numbers tested. The inset represents a zoom of CD values for low angles of attack and high Reynolds numbers, together with those reported by Ngo and Barlow [23], Mueller and Torres [6] and Sheldahl and Klimas [17]

The lift coefficient versus AoA is shown in Figure 4, and it presents a reduction of 40% in the maximum lift coefficient C_Lmax in comparison to the values measured by Sheldahl [17] for a 2D model. The maximum lift coefficient C_Lmax has values between 0.52 and 0.61 for any value of the Reynolds numbers, and the stall angles stall appear between 12 and 14 degrees. These values show a slight deviation compared to those published by Ngo [23] with AR = 4. As expected, our stall angles are slightly greater than Ngo's results due to the aspect ratio reduction (AR = 2 in our case). It is worth mentioning that the values of C_L are lower than those obtained by Mueller [6] for a flat plate with the same aspect ratio and Reynolds numbers. Furthermore, β_L defines the ratio between the lift coefficient and the AoA, ΔC_L/Δα, see the inset in Figure 4.

In this work, we suggest a classification of the wing aerodynamics defining β_L as function of AoA. We observe two different slopes, so we split β_L into two, depending on the value of α: β_L1 (for null AoA up to α₁ ≈ 3º) and β_L2 (from α₁ till the stall angle). β_L1 and β_L2 for the two highest Reynolds numbers are plotted in the inset of Figure 4 together with the results from Prandtl's lifting line theory for symmetric, rectangular wing with no twist and finite span, C_L ≈ mα, with m = 2π/[1 + 1/(sAR)] [26]. The slopes β_L1 for highest Reynolds numbers are 0.059 deg^-1 and 0.063 deg^-1 for Re_c = 1·10⁵ and Re_c = 1.33·10⁵, respectively. These slopes at the origin are far from those reported by Sheldahl [17] for a 2D model and Re_c of the same order of magnitude (β_L1 = 0.11 deg^-1). However, these data are practically in agreement with those reported by Mueller [6] for small AoA and for a flat plate with AR = 2 (β_L1 = 0.05 deg^-1) and with those obtained by Ngo [23] for a NACA0012 airfoil with AR = 4 (β_L1 = 0.075 deg^-1).

Figure 4. CL vs α for all Reynolds numbers together with those reported by Ngo and Barlow [23], Mueller and Torres [6] and Sheldahl and Klimas [17]. The inset shows detail of βL1 and βL2 for the two highest Reynolds numbers in comparison to the result from Prandtl’s lifting line theory for symmetric, rectangular wing with no twist and finite span [26]

For the two lowest Reynolds numbers, the linear region of the lift curve (before stall) has a discontinuity at small AoA, see Figure 5. The values of β_L1 correspond to the slope before the discontinuity appears, so the values are 0.032 deg^-1 and 0.048 deg^-1 for Re_c = 3.33·10⁴ and 6.67·10⁴, respectively. On the one hand, these slopes show an excellent agreement with those reported by Laitone [22] for similar Re_c, in a NACA0012 wing with AR = 6 and AoA close to zero. On the other hand, the lift coefficients curves show a similar behaviour to those achieved by Laitone: as the free stream velocity increases, the discontinuity of slope β_L1 is smoothly reduced, so the final tendency is a constant slope at small AoA (0º ≤ α ≤ 3º). Once parameters β_L1, β_L2 are defined and characterized, we will discuss later our curve fitting as function of Re_c and any value of the aspect ratio AR. We first analyse the minimum value of C_D.

The minimum drag coefficient decreases as Re_c increases. The experimental points depicted in Figure 6 (a) are adjusted using the following equation C_Dmin = 0.55·Re_c^-0.21. Laitone [22] obtained the equation C_Dmin = 0.35·Re_c^-0.25 to correlate the minimum drag coefficient for 2·10⁴ ≤ Re_c ≤ 7·10⁴. The main cause of this difference, especially in the range of C_Dmin, are the different aspect ratios used in both experiments and the different range of Reynolds numbers. Nevertheless, we suggest a general correlation using the following expression for any value of AR: C_Dmin·AR = 2·Re_c^-0.25. This curve fitting is plotted in Figure 6 (b). There is a reasonable good agreement for the data reported by other authors [22, 14, 17 and 9] for different AR and Re_c. We present only a disagreement with those results reported by Ngo [23] and Yen [25]. The discrepancy may be explained in terms of edge effects, surface roughness or rounded shape, the minimum offset used in the force sensor or the turbulence intensity in the wind tunnel. Nevertheless, it is worth mentioning that four different experimental arrangements give data that collapse in one curve fitting.

Figure 5. CL vs αfor Rec = 3.33·104 and Rec = 6.67·104 with AoA values lower than the stall angle (AR = 2), together with Laitone’s data [22] (AR = 6, I = 0.02%)

Figure 6. (a) CDmin vs Rec for 3.33·104≤Rec≤1.33·105 and AR = 2 compared with the results of Ngo and Barlow [23], Yen and Huang [25], Laitone [22], Lee and Pereira [9] and Sheldahl and Klimas [17]. (b) CDmin · AR vs Rec for the universal NACA 0012 wings fitting law CDmin = 2/AR·Rec−0.25

Other general correlations for β_L1 and β_L2 as function of Re_c for NACA0012 finite wings and any AR are shown in Figure 7. For Re_c lower than 1·10⁵ the β_L1 slope is equal to β_L1·AR^-0.5=0.00012·Re_c^0.51, and β_L1·AR^-0.5 ≈ 0.042 for Re_c greater than 1·10⁵ whilst β_L2·AR^-0.5 = 0.056·Re_c^-0.07 for any value of Re_c. In the case of our curve fitting β_L1·AR^-0.5, there is an excellent agreement for all the cases reported. The function β_L1·AR^-0.5 presents a linear dependence on the Reynolds number up to Re_c ≤ 1·10⁵, but for greater values the function β_L1·AR^-0.5 remains almost constant. It is clear that there is a critical value of the chord based Reynolds number at Re_c = 1·10⁵, where β_L1 finds a saturation value which strongly depends on the aspect ratio. Regarding β_L2·AR^-0.5, there is a reasonable good agreement with all the data reported for several AR.

Finally, the polar curve C_L/C_D against AoA is plotted in Figure 8 together with C_D as function of C_L. The maximum value C_L/C_D increases with Re_c. However, once Re_c is greater than 1·10⁵, there is a variation in the upward trend, so that the maximum value decreases with Re_c instead of increasing its value. Obviously, this critical Reynolds number tends to coincide with the value at which the ratio β_L1 saturates, giving validity and reliability to the curve fitting presented in Figure 9 (a). The maximum efficiency is found in the polar curve between 10º and 12º, depending on Re_c. In addition, C_Lmax/C_Dmin is plotted against Re_c in Figure 9 (a) together with those data provided by other authors [23, 25, 22, 9 and 17]. Again, we suggest a general correlation for NACA0012 airfoils and any value of the aspect ratio, AR, the final expression being C_Lmax/C_Dmin = 0.27·AR·Re_c^0.25. There is a reasonable good agreement, so most data collapse in one curve fitting, see Figure 9 (b).

Figure 7. Curve fitting for βL1 ·AR-0.5 vs Rec (a) and curve fitting for βL2·AR-0.5 vs Rec (a) for 3.33·104 ≤Rec ≤1.33·105 and AR = 2 (b) together with those results reported by Ngo and Barlow [23], Yen and Huang [25], Laitone [22] and Lee and Pereira [9]

Figure 8. Lift-Drag ratio CL/CD vs angle of attack α. The inset shows CD vs CL for all Reynolds numbers tested

Figure 9. (a) CLmax/CDmin vs Rec for 3.33·104 ≤ Rec ≤ 1.33·105 and AR = 2 compared with the results of Ngo and Barlow [23], Yen and Huang [25], Laitone [17], Lee and Pereira [9] and Sheldahl and Klimas [17]. (b) CLmax/CDmin·AR−1 versus Rec for the fitting law CLmax/CDmin=0.27·AR·Rec0.25

4. Conclusions

The precise dynamic response of low-aspect-ratio NACA0012 airfoil (aspect ratio, AR = 2) has been characterized by means of a digital force sensor at low-to-moderate Reynolds numbers. We compute drag and lift coefficients through the temporal evolution of force measurements. We show in this experimental study that wingtip vortex generation and low-aspect-ratio are the main causes of the lift force reduction in comparison to the NACA0012 (bidimensional) infinite profile, being this reduction almost 40% for any value of the Reynolds number. In addition, the edge effect generates an induced drag force which provides a value six times greater than the (bidimensional) infinite profile. Finally, and regarding the aerodynamic characteristics, the critical (stall) angle, α_stall, takes place at 12 degrees approximately, and this angle increases with the Reynolds number slightly.

On the other hand, the relation between C_Dmin and Re_c has been experimentally adjusted and verified with previous studies for any AR by the expression C_Dmin = 2/AR·Re_c^-0.25. Besides, we pay our attention on the lift coefficient, so we define the ratio between the increment of lift coefficient and the angle of attack, ΔC_L/Δα, namely β_L1 for 0º ≤ α ≤ 3º and β_L2 for 0º ≤ α ≤ α_stall. Our main finding of this experimental study is the general curve fittings for β_L1 and β_L2, and any value of AR. We find that β_L1· AR^-0.5 strongly depends on Re_c up to Re_c ≈ 1·10⁵ (β_L1·AR^-0.5=0.00012·Re_c^0.51), being β_L1·AR^-0.5 almost constant for Re_c greater than 1·10⁵ (β_L1·AR^-0.5 ≈ 0.042). On the other hand, β_L2·AR^-0.5 is also constant for all Re_c considered in this work (β_L2·AR^-0.5 = 0.056·Re_c^-0.07). In addition, we propose a general correlation for the relation C_Lmax/C_Dmin and any AR: C_Lmax/C_Dmin = 0.27·AR·Re_c^0.25.

Finally, we do observe that the critical chord based Reynolds number Re_c ≈ 1·10⁵ at which β_L1 saturates its value is connected to the change of the upward trend in the maximum value of the polar curve. More effort must be done regarding forced vibration response. This is a work in progress.

ACKNOWLEDGEMENTS

This research has been supported by Proyecto Junta de Excelencia de la Junta de Andalucía, Grant number P11-TEP-7776.