## Abstract

The flight of many birds and bats, and their robotic counterparts, occurs over a range of chord-based Reynolds numbers from 1 × 10^{4} to 1.5 × 10^{5}. It is precisely over this range where the aerodynamics of simple, rigid, fixed wings becomes extraordinarily sensitive to small changes in geometry and the environment, with two sets of consequences. The first is that practical lifting devices at this scale will likely not be simple, rigid, fixed wings. The second is that it becomes non-trivial to make baseline comparisons for experiment and computation, when either one can be wrong. Here we examine one ostensibly simple case of the NACA 0012 aerofoil and make careful comparison between the technical literature, and new experiments and computations. The agreement (or lack thereof) will establish one or more baseline results and some sensitivities around them. The idea is that the diagnostic procedures will help to guide comparisons and predictions in subsequent more complex cases.

## 1. Introduction

### 1.1. The aerodynamics of small wings

The interface between natural and artificial flyers is becoming blurred as robotic devices either exist or have been proposed for almost the entire range of speeds and length scales of birds, bats and many insects [1–5]. Many of these examples have quite complex flapping and/or flexible wing structures and kinematics, though the aerodynamic principles upon which the designs rest are typically those inherited from standard aeronautics where the sizes and speeds are much larger. In the limit of simple fixed wing flight at steady speed *U*, with characteristic length scale in the flightwise direction, *c*, the Reynolds number is
1.1where ν is the kinematic viscosity. *Re* shows the relative importance of inertial versus viscous terms in the Navier–Stokes governing equations, and in classical aerodynamics is a large number, perhaps ranging from 10^{6} to 10^{8}. In such cases, the effects of viscosity are either easily ignored entirely or limited to a thin region close to the lifting surface, the boundary layer, where the outer potential flow adjusts to the presence of a solid body. When *Re* = 10^{5}, it is still a large number, but if based on length and velocity scales inside the boundary layer, is no longer overwhelmingly large. As remarked in [6], the performance of wings at moderate *Re* is then substantially dictated by the dynamics of the laminar boundary layer, and most particularly, whether it separates, and perhaps then reattaches. When and if the separated shear layer does reattach, at least in some mean sense, then the resulting laminar separation bubble (LSB) dynamics can be quite influential in determining global force coefficients, both instantaneous and time-averaged. There is a range of *Re* where LSB formation is likely because *Re* is low enough to make separation, at some chordwise location, inevitable, and *Re* is also high enough so that transition to turbulence can occur quickly in the separated shear layer, which increases the likelihood of reattachment. In this sense, we use the term ‘moderate *Re*’ to describe a regime that is dominated neither by viscosity nor by inertial terms, but rather by the balance between the two, which can vary greatly with boundary conditions and locally in the flow.

The global force coefficients are expressed in the standard way, so the resultant net aerodynamic force on an aerofoil section, or wing, is described through the lift and drag components, normal and parallel to the mean flow:
1.2where *L′*, *D′* are the lift and drag per unit span for a two-dimensional (2D) or infinite wing, *c* is the chord length and *q* = ½ *ρU*^{2} is the dynamic pressure. The equivalent expressions for a finite wing, or spanwise section of a wing, are
1.3where *S* is a planform area. For a given section geometry, the aerodynamic force coefficients on a wing section are then functions of only two parameters,
1.4the Reynolds number and geometric angle of attack, *α*.

The sensitivity to boundary layer separation and reattachment can make accurate and repeatable measurement, even of time-averaged *C*_{l} and *C*_{d} difficult, as shown in figure 1 [7], where the similarity between measured *C*_{l}(*C*_{d}) polars breaks down as *Re* drops below 10^{5}. The flow over an Eppler 387 aerofoil is sensitive to small variations in geometry and environmental conditions. The Eppler 387 is not designed to operate at such small *Re*, but the case does illustrate a rather common phenomenon for a number of different section geometries [8,9]. The extreme sensitivity can not only be problematic, but also can present an opportunity if it can be reliably exploited for control, and a number of publications by Yang & Spedding [10–12] show how the sensitivity to either passive or active acoustic forcing can be used to produce local sectional changes in *L*/*D* of up to 80% (figure 2).

### 1.2. The NACA 0012 as a standard test case

There is great uncertainty over rather basic quantities such as time-averaged lift and drag for wings and wing sections, simply because the Reynolds number lies in a particular range. This unfortunate state is a consequence of the high sensitivity to small variations in conditions, and also to the fact that most testing facilities are not well suited to measuring such small forces with small tolerances. Furthermore, until the recent advent of small robotic flying devices, there had been little practical incentive to spend resources in this niche problem. If a sound technical database is required, then there is much work to do. Not only that, but it turns out that many sensitive problems of laminar boundary layer stability, separation, transition and reattachment combine in complex ways that are intrinsically interesting. The *Re* range from *Re* = 1 × 10^{4} to 1.5 × 10^{5} has a number of subdomains, and now from the lower *Re* upwards, can be accessible to direct computational approaches. We, therefore, make a start by selecting a common and simple aerofoil section, the NACA 0012, and attempt a comprehensive summary.

A thorough review of experimental studies in more than 40 wind tunnel test facilities by McCroskey [13] found that there were quite substantial disagreements between experimental results even at *Re* > 10^{6}, and even though criteria were suggested for reasonable agreement at higher *Re*, for lower *Re*, the data were almost absent and agreement worse. In a similar vein, Ladson [14] found variations with *Re* and Ma (Mach number) in the Langley test facility for *Re* ≥ 2 × 10^{6}, but noted that when roughness fixed the transition point, differences were far less evident. The particular difficulties of aerodynamic force measurement at low *Re* have been pointed out before [15–18]. Added to the sensitivities to freestream conditions and surface geometry, there are compounding factors of surface finish/polish, acoustic environment [19] and non-negligible influence of viscous corner flows on endplate configurations [20]. A number of experimental studies [21–23] on the NACA 0012 near *Re* = 5 × 10^{4} have shown nonlinearities in the curves of *C*_{l}(*α*), even at small *α* when thin aerofoil theory would predict a linear relation of *C*_{l} = 2*πα*. The nonlinear shapes have not been the same from study to study, and no consistent explanation has been given for the phenomenon based on flow physics. An exception is [23], but their results will be seen to be quite different from those given here, perhaps because of the relatively higher turbulence levels. The nonlinearities may sometimes be evident only in small *α* intervals, and many studies lack such resolution [24,25]. This latter category includes all those operating at the standard 1 degree *α* interval.

Aeronautical flows at transitional *Re* are becoming accessible to computations of varying degrees of fidelity, though full direct numerical simulations (DNS) remain expensive/lengthy, so approximations are sought. In the context of this paper, we confine our attention to time-averaged force coefficients and pressure distributions. Even this modest goal is not easy, for the same reasons that physical experiments are not. The most rapid calculations are based on inviscid solutions with varying degrees of viscous correction. Of these, the XFOIL panel code [26] is the best known and most accessible, and also has a quite elaborate viscid–inviscid coupling that allows both transonic and low Reynolds number flows to be modelled [27] with some success. (One of the test cases was the E387 at *Re* = 6 × 10^{4}, Langley data plotted in figure 1.) Computing the forces with a surface-based method is efficient for a simple shape in a uniform flow. Vorticity originates only at the fluid-surface boundary, and the majority of the far field is well approximated by potential flow.

However, computational fluid dynamics methods are well packaged and readily adaptable to these cases too. The Navier–Stokes equations (or a convenient, discrete form of them) can be solved directly on structured or unstructured grids, with varying degrees of approximation. Lee *et al*. [28] compared various numerical methods for the NACA 0012 at moderate *Re*, including results from 2D laminar codes (with no turbulence model), 2D Reynolds-averaged Navier–Stokes (RANS) with turbulence modelling (Baldwin–Lomax) and an implicit three-dimensional (3D) large eddy simulation (iLES). The computations were described as satisfactory, but one may only be able to say so after the fact, when it is clear which model works best. For example, *C*_{l} at *α* = 4.5° varied from 0.22 (2D RANS) to 0.52 (3D iLES)—for reference, 2*πα* = 0.49 for *α* = 4.5°. The dangers of running a turbulence model at lower *Re* than its design point have been pointed out in [29], who stressed that successful matching of the observed physics will then also require some kind of transition model. In cases of relaminarization in transitional flows, the predictions are often inconveniently sensitive to model parameter selection. More recently, the detailed effect of numerical dissipation in immersed boundary/LES methods has been examined in [30,31]. The numerical dissipation may exceed viscous or sub-grid scale model dissipation in transitional flows with laminar separation bubbles, so turbulence quantities and dissipation rates may agree with other computations or experiment only after suitable tuning of model coefficients.

The most extensive calculations to date (and part of the reason for this test case choice) are the DNS of Jones *et al*. [32], who used a fourth-order central difference scheme (with no modelling) to compute the flow over an NACA 0012 at *Re* = 5 × 10^{4} and *α* = 5°. The full 3D calculations were initialized by expanding initial 2D simulations, and the naturally observed process of transition was stimulated by a 3D forcing, which was then relaxed after some time. The 2D *C*_{L} values fluctuated around 0.5, and then rose to 0.6 when forcing was applied in the 3D case. After a transient, *C*_{L} varied between 0.6 and 0.64. The unforced 3D flow was not steady, but transition around the periphery of the separation bubble was self-sustaining with persistent pressure fluctuation amplitudes. The separated shear layer from the bubble was identified as a possible source of absolute instability that could account for the self-sustained turbulence. These computations were performed at a single *α*, with no appropriate experimental data for comparison. The same case (NACA 0012 at *α* = 5°) was re-examined by Almutairi *et al*. [33] who compared the original DNS results with filtered DNS and LES models and with a viscous–inviscid interaction model that allowed coupling of an outer potential flow solution with boundary layer models, quite similar to the XFOIL formulation except including unsteady terms. In the LES, the bubble length, growth and bursting and fluctuating and mean *C*_{l} were all sensitive to the selected spanwise domain length of 0.2*c* (as in the original [32]) or 0.5*c*. It was concluded that this sensitivity would be expected whenever the LSB occupies a significant fraction of the chord. The comparisons of unsteady mixed viscous–inviscid models with LES were promising, suggesting that considerably cheaper computations could be made at this *Re*, hence over a range of *α*. Again there were no experiments and model results were mostly restricted to *α* ≥ 9°.

### 1.3. Objectives

Experimental measurements and computations of aeronautical flows in the transitional Reynolds number regime are extremely challenging, characterized by enormous sensitivity to small disturbances or variations in boundary conditions. Few have paid attention to the scarcity of reliable information until recently as studies of natural and artificial flyers overlap this *Re* domain. The purpose of this paper is to present reliable data for one specific case that we propose as a canonical test case for codes and experiments that aim to make predictions of aerodynamics in such a range of *Re*. Although ostensibly simple, the boundary layer and separation dynamics are quite subtle, leading to shapes of time-averaged force coefficient curves that would not readily have been predicted. The second purpose is to investigate how simple computations perform and why, with a view to suggesting the most productive future routes for design codes that must run fast.

## 2. Material and methods

### 2.1. Wind tunnel experiments

#### 2.1.1. Model and wind tunnel

An NACA 0012 wing was milled from solid aluminium, with chord, *c* = 7.5 cm and span, *b* = 48 cm, for an aspect ratio AR = *b*/*c* = 6.4. (Most tests reported here are for the 2D configuration where *AR* is not explicitly a parameter.) Tests were carried out in the closed-loop Dryden Wind Tunnel at the University of Southern California. The tunnel has a contraction ratio of 8 : 1 and an octagonal test section measuring 1.37 m wall to wall. Possible blockage effects of the finite volume model in the tunnel cross-section were estimated and found to be small. A correction suggested by Barlow *et al.* [34] for unusually shaped objects where the model volume includes wing, supports, endplates and shroud at maximum *α* leads to a correction in the effective freestream, *U* of 0.4%. An empirical procedure described in [7] yields a blockage correction factor, *ɛ*_{sb} = 0.0044. An equivalent estimate of the wake blockage, *ɛ*_{wb} = 0.0024, and that due to streamline curvature is approximately 0.0007. The combined sum of the maximum likely blockage effects yields a maximum effect on force coefficients (*C*_{l} and *C*_{d} are treated separately) of 1.8%. Twelve screens reduced the turbulence levels (*T* = *q*/*U*, where ) in the test section to less than 0.03% for spectral frequencies between 2 and 200 Hz over the speed range of 5–26 m s^{−1} [35,36].

For 2D tests, an infinite aspect ratio was approximated by placing endplates at either tip of the wing. The endplates were 9.48*c* (in *x*) × 4.75*c* (in *z*) × 0.17*c* thick, with sharp leading edges and were aligned carefully, parallel with the flow. The ends of the model were kept within about 1 mm of the endplates. This is within the separation distance of 0.005*b* (2.4 mm) recommended by [34], and less than the estimated laminar boundary layer thickness on the plates themselves at the leading edge of the model, *δ* = 5.2*x*/*Re*^{½} = 2.7 mm. The model support rod passed through a hole in the bottom endplate 1.13*c* behind the leading edge and equidistant from the edges in *z*. The bottom end of the model itself is 0.85*b* above the wind tunnel floor.

#### 2.1.2. Force balance set-up and calibration and test procedures

The business of making lift and drag estimates in transitional flow regimes is sensitive to numerous influences, and great care was taken in model set-up, data acquisition and analysis to assure that data were reliable and that experimental uncertainties were well characterized.

The force balance measurements were performed with a custom, three-component, cruciform-shaped force balance with a parallel plate sandwich design [36]. A new static calibration was performed before each test, generating a 3 × 4 calibration matrix. The three most recent matrices were averaged to generate the final calibration matrix used during the test. Precautions were taken to assure that sensitivities of the estimated drag to off-diagonal terms in the calibration matrix were correctly controlled, as noted below. The uncertainty of lift and drag measurements is estimated to be less than 8 mN, which is 0.13 of the minimum expected drag force on this sized model.

The force balance is located below the tunnel and is connected to the model via a sting that extends through the tunnel floor and is shielded by an aerodynamic shroud. The model support rod was inserted into the sting and secured with a screw. Lift and drag measurements were zeroed with forces corresponding to the empty sting, the weight of the model and free stream flow interaction with the sting. Because the sting was shielded by a shroud, the free stream flow interaction component was much smaller than the others, generally less than 5 mN. For 3D tests, in order to reduce the interference of the shroud with the lower wing tip vortex, the length of the shroud was reduced so that the lower end of the model was raised 1.5*c* above the top of the shroud, exposing the support rod to the air flow. To correct for this, the drag force on a matching rod was also measured and subtracted from the test results.

Each force balance test consisted of 10 sweeps, five forward and five backward, through an angle of attack (*α*) range of −5° ≤ *α* ≤ 9° in increments of 0.5°. After each *α* step, the flow was allowed to settle for 10 s before 10 s of data were collected at 1 kHz and averaged to produce one time-averaged measurement. The 10 sweeps produced 10 time-averaged measurements for each *α*, which were averaged to yield a single value. The uncertainty at each *α* was taken as the standard deviation of the 10 measurements at that *α*.

Although the curve shapes were consistent from test to test, drag results were sometimes found to be non-symmetric. It was discovered through repeated testing that the drag results are extremely sensitive to the off-diagonal terms of the calibration matrix, which were caused by small alignment errors during calibration. Because a symmetric model was being used, non-symmetric drag results were identified as incorrect and discarded. The lift curve was relatively insensitive to these small calibration errors, and there was generally a negligible difference between lift results from test to test. Tests were generally performed on separate days, but even when two tests were performed on the same day, the entire procedure, including calibration and collection of zeroing forces, was repeated for each test. In this case of a symmetric model, curves were shifted by small *α* (less than 0.3°) to ensure zero net lift force at zero *α*. Close to stall, the model begins to oscillate considerably, so measurements were not taken for post-stall *α*.

#### 2.1.3. Particle image velocimetry measurements

Particle image velocimetry (PIV) tests were carried out for 0° ≤ *α* ≤ 8°. The tunnel was filled with glycerin-based smoke with a typical particle diameter of 0.2–0.3 µm and a laser sheet parallel to the flow direction (in {*x*, *z*}) was generated by a Quantel EverGreen double-pulsed Nd:YAG laser. An Image Pro X 2M CCD camera (1600 × 1200 pixel, 14 bit) imaged particle fields on a cross-section of 2.5 cm (0.05b) above mid-span in {*x*, *z*} with a Nikon 70–210 mm f/4–5.6 NIKKOR AF lens. Because the model was symmetric, the suction and pressure sides were illuminated by rotating the model in positive and negative *α*. To increase spatial resolution, the flow field on each side of the aerofoil was split into two, slightly overlapping sub-regions that were imaged in separate experiments. Sequences of 200 image pairs were captured for each sub-region at a sample rate of 9.6 Hz, and the time delay between images in an image pair (*δt* = 8–30 µs) was carefully tuned to maximize the dynamic range of observable displacements. Minimizing the effect of peak locking errors tends to increase *δt*, while minimizing the frequency of occurrence of untrackable shear deformations of correlation boxes inside the separation bubble tends to decrease *δt*. The optimum *δt* is different for each *α*. The images were processed with LaVision's DaVis software to produce velocity field estimates {*u*, *w*} in the streamwise and vertical directions {*x*, *z*}, on a uniform grid using a multi-pass algorithm, with initial 64 × 64 pixel interrogation windows reducing to 16 × 16 pixels by the final pass. A 50% overlap gave a final spatial resolution of 8 pixels, which is 0.27 mm, or 0.0036*c*.

All 200 instantaneous velocity fields were averaged to produce one time-averaged velocity field for each sub-region. A built in Matlab thin-plate smoothing spline with a single smoothing parameter, similar to the spline used in [37], was applied to the averaged results in order to reduce random noise. The average change in either velocity component in a sub-region due to smoothing was always less than 0.4% of the maximum value of that component in the sub-region. The spanwise component of vorticity, *ω _{y}* = ∂

*w*/∂

*x*− ∂

*u*/∂

*z*, was calculated at each grid location from the derivatives of the smoothing spline coefficients. All four averaged and smoothed sub-region velocity and vorticity fields were finally stitched together to form one composite velocity/vorticity field, on both sides of the aerofoil, for each

*α*.

### 2.2. Numerical experiments

The commercial package Star-CCM+ was used to run RANS simulations for a 2D NACA 0012 geometry, as a worked example of how such codes behave given reasonable user choices. For numerical convenience, the aerofoil section was truncated at 0.99*c* at the trailing edge. The aerofoil was modelled in two ways; a circular 2D plane as shown in figure 3*a* and a thin slice 3D, hereafter referred to as pseudo-2D or P2D, domain with symmetry planes on either side of the domain and spanwise thickness of 0.55*c*, as shown in figure 3*b*. Polyhedral mesh shapes were used with 22 cells across the boundary layer. The cell size on the surface of the aerofoil was 0.0005*c* which expanded into the outer domain to 0.1*c*, except in the wake where it was 0.025*c* (figure 4). The sufficiency of the mesh resolution was checked through a standard grid convergence index on a sequence of three increasingly fine resolutions, ending with a final mesh count of 10^{6} and 2 × 10^{6} cells, respectively.

The flow was assumed to be steady and incompressible, with a constant velocity inlet boundary condition. The shear-stress transport *k*–ω turbulence model [38] was used with an additional *γ–Re _{θ}* transition model [39] specifically formulated for unstructured CFD codes to predict laminar to turbulent transition, known to be a requirement in problems with laminar separation and possible reattachment. A sensitivity analysis for the 2D model on the various physics parameters required (a turbulence intensity, a turbulence viscosity ratio) showed less than 0.5% variation in total aerodynamic forces for the NACA 0012 aerofoil at

*Re*= 10

^{5}and

*α*= 5°.

## 3. Results

### 3.1. Wind tunnel measurements

Figure 5*a,b* shows *C*_{l}(*α*) and *C*_{d}(*α*) for the wind tunnel experiments. There are a number of features of *C*_{l}(*α*) alone that are notable. First, about *α* = 0°, the lift slope, d*C*_{l}/d*α*, or *C*_{l}* _{,α}*, is negative.

*C*

_{l}reaches a local minimum value at

*α*= 0.5° and then increases with a slope significantly above the theoretical thin aerofoil result (

*C*

_{l,}

*= 2*

_{α}*π*) up to

*α*= 3°. At this point

*C*

_{l}exceeds the 2D theoretical value. With further increase in

*α*up to 9°,

*C*

_{l,}

*< 2*

_{α}*π*, and at higher

*α*the aerofoil begins to stall, but not abruptly. The resolution in

*α*is only just sufficient to show the negative

*C*

_{l,}

*about*

_{α}*α*= 0°, but the result is robust and repeatable. The inset of figure 5

*a*shows a separate set of experiments from −1.6° to +1.6° in steps of 0.2°. Table 1 shows slopes for linear least-squares fits through the data over the three characteristic regions. The slopes are different from each other and always different from 2

*π*, the inviscid, thin aerofoil value. The data agree with the

*C*

_{l,}

*= 2*

_{α}*π*line only coincidentally, at the three points where the curves intersect.

*C*

_{d}(

*α*) (figure 5

*b*) does not have discontinuous regions similar to

*C*

_{l}. It is also symmetric about

*α*= 0°, and could be reasonably fit within uncertainties with a smooth function.

Since *C*_{l} is negative at small positive *α*, so is *L*/*D* (figure 6*a*). There is a broad maximum in *L*/*D* from about *α* = 3–7°. Its value, from 13 to 15, is respectable for a moderate *Re* wing (Lissaman [6] shows data where expected (*L*/*D*)_{max} might be approximately 10 at this *Re*, but where the variance can be high), mainly because the lift is higher than expected, and at the same time, *D* is not significantly over-estimated. The lift-drag polar (figure 6*b*) has a characteristic loop with *C*_{l} = 0 at three different *C*_{d} points. d*C*_{l}/d*C*_{d} is also quite steep, up to *C*_{l} = 0.4. All observations exceed experimental uncertainty, and the curves are symmetric within those bands about *α* = 0°.

The observations from figures 5 and 6 are not peculiar to the 2D case (as simulated with endplates), but are just as evident for the finite wing, AR = 6.4 geometry shown in the equivalent figures 7 and 8. *C*_{L} is measurably and repeatably negative for small positive *α*; the three, almost-linear slopes all differ from the classical, inviscid value,
3.1where the 2D value, *C*_{l,α} = 2*π* is decreased as AR decreases (we ignore the correction for span efficiency, which is close to 1), though it remains a good average value over all *α*. From *α* = 3° to *α* = 6°, *C*_{L} significantly exceeds the theoretical value. *C*_{D} is not, in general, appreciably higher than *C*_{d} so *L*/*D* (figure 8*a*) again has a broad peak in *α* = 3–6°. The loop in the lift-drag polar (figure 8*b*) remains as distinct as for the 2D case.

PIV-derived, time-averaged fields of |*u*|(*x*, *z*), *ω*_{y}(*x*, *z*) (figure 9) and *u*(*x*, *z*), *w*(*x*, *z*) (figure 10) explain the force balance observations. Figures 9 and 10 are near mid-span sections through the AR = 6.4 wing. Based on the similar shapes of the force balance data, for example, *C*_{l}(*C*_{d}) (figure 6*b*) and *C*_{L}(*C*_{D}) (figure 8*b*), the 2D and finite wing flow fields are not expected to be significantly different far from the tips, and the bubble dynamics are not appreciably different, consistent with observations in [35] for an E387 at similar *Re*.

At *α* = 0°, the flow about the NACA 0012 is symmetric and separates (on both sides) before the trailing edge. At *α* = 0.5°, the separation point has moved forward on the upper (suction) surface, but has moved aft in the lower (pressure) surface (figure 9*a*, row 2). The contour of zero spanwise vorticity no longer leaves the trailing edge straight, but is deflected upwards (figure 9*b*, row 2). Regions of *u* < 0 are more prominent on the upper surface (figure 10*a*, row 2), and the distribution of *w*(*x*, *z*) aft of the trailing edge is asymmetric, with stronger, positive *w* on the lower side (in figure 10*b*, row 2, *w* = 0.25 m s^{−1} beneath the aerofoil, and −0.2 m s^{−1} above it). At this *α*, the lowest part of the aerofoil is below the trailing edge and a slightly favourable pressure gradient allows the streamlines to follow the curvature into an upward direction. The net acceleration is upwards, as the lower streamlines have higher curvature, and the net lift is negative. At the same time, the laminar separation before the upper surface trailing edge assures that streamlines here are also deflected slightly away from the upper surface. The effect here is similar to an aerofoil with trailing edge reflex, but it also includes an important contribution from the upward sweep around the aerofoil thickness on the pressure side.

At *α* = 2° (row 3 in figures 9 and 10), the streamlines are not deflected strongly upwards on the lower surface, and though the separation line has moved further forward on the upper surface, the net flow has returned downwards. The trailing wake has its smallest streamwise extent, like an attached recirculation bubble. (The flow is highly unsteady here, and explanations on time-averaged fields need careful interpretation.) At *α* = 4°, (row 4 in figures 9 and 10), all fields show signatures of a separation bubble that reattaches on the upper surface, close to where *w*(*x*, *z*) has its highest negative value.

The flow acts as if the aerofoil had a higher convex curvature on the upper surface, formed by the combination of aerofoil surface and separation bubble. The changed effective camber accounts for the better than inviscid thin aerofoil *C*_{l} around this *α* (figure 5*a*). The drag cost is comparatively small and this *α*-range marks the beginning of the broad maximum in *L*/*D* (figure 6*a*). If the separation bubble can be termed an LSB in the classical sense, then the LSB is associated with improved *L*/*D* over this range of *α*, which is contrary to most literature interpretations based on observations at higher *Re*. As the LSB moves forward on the upper surface with further increases in *α* (6° and 8°, rows 5, 6, in figures 9 and 10) the LSB shrinks in streamwise extent and the spanwise vorticity becomes more strongly negative at the outer shear layer. This phase is associated with reduced *C*_{l,α}, but continued high *L*/*D*. There is a broad downwash region in *w*(*x*, *z*) from the mid-point of the LSB to the aerofoil trailing edge.

### 3.2. Numerical simulations

Neither 2D nor P2D simulation agrees with each other, or with the theoretical curve for *C*_{l}(*α*) (figure 11*a*). Normally, *C*_{l}(*α*) is considered a simple calculation, determined by the pressure field, which at higher *Re* is well approximated by potential flow. Here, the lack of agreement shows that the RANS performance is largely determined by boundary layer and model coefficients, and not by the outer potential flow. Figure 11*b* shows the friction and pressure components of the drag for the 2D and P2D calculations. The friction drag decreases slowly as *α* increases in both cases, with similar magnitudes. The decrease in *C*_{f} with *α* is presumably associated with the thickening boundary layers and reduced shear stress at the wall. *C*_{p}, on the other hand, increases strongly with *α* and the two computations do not agree by how much.

There is another problem with a RANS-type calculation in that the true flow field is not steady, even at small *α*. A RANS estimate will converge to some solution, but that state is not necessarily a good measure of a time-average. Figure 12 shows the variation in *C*_{l} with iteration number, first for the 2D case and then for P2D. In P2D, the 3D flow field can stabilize the solution, and the fluctuation amplitudes are reduced, but the flow remains unsteady. The complicated three-dimensional structure of the RANS solution is seen in figure 13, where the surface distribution of *C*_{f} is mostly 2D (uniform in span) up until the separation point at about mid-chord. After, there are streamwise streaks with spanwise wavelength small compared with the chord, and then towards the trailing edge there is a much more irregular footprint.

The mean streamwise velocity fields for the *α* range {0°, 0.5°, 2°, 4°, 6°, 8°} (figure 14) are different for the 2D and P2D case. Except at *α* = 8°, the 2D flow field has numerous vortex-like structures that cannot be steady and cannot be time-averages. The P2D fields look more likely to be converged on an average. Qualitatively, the *u*(*x*, *z*) fields do not look similar to the equivalent experiments (cf. Figure 9) as the LSB signature is much reduced in size, or absent. In this respect, though much care has been taken in experiment to not over-smooth the data during spline interpolation, recall that all PIV data are spatial averages over the smallest correlation box, which is 16 pixels, or two grid points (0.54 mm) in size. This compares with an estimated LSB height (normal to the surface) of between 7 and 11 grid points for *α* = 4–11°, so is small but may smear out the apparent height by a factor of 2/7 at worst.

In numerical results, the first separation and reattachment points can be defined by successive zero crossings for the wall shear stress, and these locations are superimposed on the maps of spanwise vorticity in figure 15. The two computations differ, as in figure 14, but the marked separation and reattachment points differ from each other and from experiment (cf. Figure 9*a,b*). The 2D simulation is not steady, at any *α*, and the P2D will be unsteady for *α* ≥ 4°, as indeed shown in figure 12. Regardless, the fields in figures 14 and 15, if not true averages, are still solutions to the RANS equations even if they are best regarded as snapshots of an underlying unsteady flow field.

## 4. Discussion

### 4.1. Negative lift slope at *α* = 0°

The lift curve slope, *C*_{l,α} < 0 around *α* = 0° for the NACA 0012 at *Re* = 5 × 10^{4}. This is a robust and quite surprising result, the opposite of any existing theoretical treatment. The correct explanation involves *Re*, the aerofoil thickness and curvature, and the natural laminar separation point at small *α*. At small positive *α*, the flow around the lower aerofoil surface is accelerated upwards because the flow turns around the lowest part of the wing, and the separation point moves back while the corresponding separation point moves forward on the upper surface, decreasing the downward-induced velocities there at the same time. In a check of existing literature, it was found that the negative lift had been seen before, but in isolated instances and with no general explanation. Most particularly [40] showed a similar shaped *C*_{l}(*α*) curve (their figure 5), with initial negative *C*_{l} for 0 ≤ *α* ≤ 3°, and high *C*_{l,α} thereafter. This pattern was observed for the symmetric, 18% thick NACA 66_{3}-018 at *Re* = 1.3 × 10^{5}, and smoke visualization experiments suggested the same fore and aft separation point variation described here. The wind tunnel turbulence level was about 0.1%. The phenomenon was completely removed by addition of surface roughness, and also by acoustic excitation. The effect of thickness in imparting a positive upward acceleration of the air at small *α* was not noted, perhaps because the visualizations were qualitative in nature. This peculiar result has stood alone in the literature—some product of the aerofoil shape and low to moderate *Re*. Here, we propose that the robust appearance of this same result for the common NACA 0012 shape implies that in careful experiments, most smooth, symmetric aerofoils with sufficient thickness (*t*/*c* ≥ 10%) will show negative *C*_{l,α} at some *Re*, depending also on the surface finish and the ambient turbulence levels. The other relevant literature studies do not contradict this contention, though examples are isolated. Tsuchiya *et al.* [22] measured *C*_{l}(*α*) for an NACA 0012 at *Re* = 4.7 × 10^{4}, and reported nonlinear *C*_{l}(*α*) curves but no negative lift. The resolution in *α* was only 1° and *T* = 0.5%. Negative *C*_{l} at *α* = 1° for *Re* = 2.5 × 10^{4} (*T* = 1%) can be seen in their data but is not discussed. In 2D Navier–Stokes calculations, Yonemoto *et al.* [41] investigated the negative lift observed in experiments by Ohtake *et al.* [42], and proposed that a reversed flow at the trailing edge upper surface could accelerate the corresponding boundary layer on the lower surface. Unique challenges with 2D computations have been noted above.

### 4.2. Laminar separation bubble improves aerodynamic performance

The behaviour of the LSB observed here—an LSB that extends over a large portion of the chord but moves forward and shortens with increasing *α*—agrees, in general, with the literature reports on wing and aerofoil performance at low *Re* [43–46]. All support the claim by Lissaman [6] that bubbles covering a large percentage of the chord can have a significant effect on aerofoil performance by altering the outer potential streamlines. It is not always acknowledged, however, that the separation and reattachment of the flow to form an LSB, while associated with lower *C*_{l,}* _{α}*, is also associated with the highest

*L*/

*D*, or aerodynamic efficiency. The increase in

*L*/

*D*can be deliberately exploited, for example, by acoustic forcing [10,12]. In both 2D and 3D cases reported here,

*L*/

*D*has a broad maximum over 3° ≤

*α*≤ 7°, when the LSB reattaches and moves forward on the upper surface (figures 9 and 10, rows 4–6). Huang

*et al.*[45] and Huang & Lee [23] made thorough associations between surface oil traces and aerodynamic performance as functions of

*Re*and

*T*, and the documentation of separation locations is generally consistent with the findings here. However, the lowest

*T*was 0.11% in [45] and 0.2% in [23], and no negative

*C*

_{l,}

*was reported at*

_{α}*Re*= 7.6 or 5 × 10

^{4}, respectively.

### 4.3. Poor agreement with previous experiments

Figure 16 compares the more complete available experimental data at moderate *Re* with the current results. Our data have uncertainty bands, but none of the other data sources do. A strict interpretation of the line plots is that none of the literature data are the same as ours. The only data points that lie inside the uncertainty range come from Tsuchiya *et al.* [22], which differ greatly at low *α*, and do not have the negative *C*_{l}. Both Tsuchiya *et al.* [22] and Huang & Lee [23] have steep initial C_{l,}* _{α}* (but no negative values). These two references agree on the higher than thin-aerofoil lift at moderate

*α*, though [25] does not. At

*α*= 2°, the variations in reported

*C*

_{l}are about 100% of the nominal theoretical value.

*C*

_{d}(

*α*) of one study [23] is significantly different, probably because of the comparatively high turbulence levels, which increase

*C*

_{d}at these

*Re*, as documented by the same authors. In general, the specific reasons for the disagreements are not clear. Table 2 summarizes the various experimental conditions. High turbulence levels in some facilities are likely to have an influence.

### 4.4. Poor agreement with computations

The issues in performing numerical simulations have been discussed in that section. Here we briefly compare the RANS results with wind tunnel equivalents, and with an XFOIL simulation. The point of the simulations is not to make the best possible (even a cursory reading of [32] and [33] will give an idea of the serious technical challenges to be overcome), but to find the outcomes when readily available commercial codes are applied to this problem. The RANS code is run on a fine mesh, and spends much time computing flows in cells that are dominated by pressure terms in basically Euler codes. The critical part is in the viscous boundary layer, and here the results are actually strongly affected by modelling of the turbulence transition, of the boundary layer itself and of the turbulence. XFOIL is in some respects a much simpler code. The outer potential flow is inviscid and a viscous boundary layer model with explicit accommodation for separation and transition is coupled with it. Both simulations are, therefore, quite dependent on the boundary layer modelling. The comparison of *C*_{l}(*α*) and *C*_{d}(*α*) is given in figure 17. The wind tunnel experiments only agree well with the RANS in *C*_{l} for 4° ≤ *α* ≤ 8°. This range of *α* is where the flow is dominated by the LSB itself, while small-*α* results are mainly influenced by details of trailing edge separation. For that part, XFOIL has much better qualitative agreement, and the three almost-linear *C*_{l,}* _{α}* segments seen in experiment are reproduced in XFOIL. However, the extra

*C*

_{l}above the theory line is over predicted, and over most

*α*, XFOIL does not give the same

*C*

_{l}as experiment.

The advanced computations of Jones *et al.* [32] were one of the main motivations for making the NACA 0012 comparison. In figure 17, the 2D computation gives *C*_{l} below the wind tunnel result, while the 3D forced and unforced cases lie above it. A strict interpretation of this figure is that no simulation has results that are the same as experiment. In *C*_{d}, there is no agreement between RANS, XFOIL and experiment. The simulations persistently predict *C*_{d} lower than experiment, and the current experiments have mostly lower *C*_{d} than others in the literature (cf. Figure 16*b*). By contrast, *C*_{d} in the 3D unforced case of Jones *et al.* [32] is equal to the wind tunnel result. This is the single clear point of agreement, but absent any others, it is not possible to assert that any of these results is a closer approximation to some baseline truth.

The overall agreement between experiment and DNS with varying degrees of sophistication is not much better than with a simpler and faster viscous boundary model coupled with outer potential flow, such as XFOIL, and also discussed in [33]. Progress in complex, high amplitude flapping kinematics typical of many natural flyers might, therefore, be effectively achieved on a modified inviscid basis, such as the models in [47,48] and then developed in [49–51].

### 4.5. How general are results at one *Re*?

We conclude that the experiments reported here do not agree in a meaningful way with any others reported in the literature, and that none agree with existing computations either. The result has been comprehensively established at one single Reynolds number for one single aerofoil shape. Of what significance then is it?

For more than 80 years, the NACA 0012 has been used as a test case, a canonical example of a smooth aerofoil. It has no special design features that make it suitable for use at moderate *Re*, though as described in the Introduction, the agreement between facilities at *Re* = 10^{6} is not what one might hope for either. Since it continues to be used as a test case for the emerging computations of moderate *Re* aerodynamics, establishing a reliable set of baseline measurements is a matter of some urgency. The difficulty in doing so is symptomatic of the challenges at moderate *Re*, for which it remains a sensitive diagnostic and test case. The mechanism proposed for the curious and counterintuitive result of the negative lift does not depend on any special features of the NACA 0012. The effective camber and reflex caused by the separation bubble, combined with the upwards flow induced by the thickness at small *α* could occur in many symmetric, smooth-surfaced, thick aerofoils, at some *Re* and in a low-turbulence environment.

The objection could be raised that no practical aircraft flies in perfectly smooth, turbulence-free conditions. That is true, of course, but there is no possibility, even in principle, of establishing a meaningful baseline and comparison between experiments themselves and between experiments and computations if these factors are not removed, controlled or decreased in influence. What remains are the sensitive dynamics of separation, transition and reattachment that lie at the heart of low- or moderate *Re* aerodynamics. Properly understanding these phenomena is what will lead to generalizable results, for birds or bats or aerial robots.

## 5. Conclusion

Serendipitously, the venerable NACA 0012 provides, at a moderate *Re* far below its original design point, a delicate and sensitive laboratory for the study of the viscous–inviscid balance that can have far-reaching effects in global aerodynamics. At present, we cannot claim to have agreement on even the integrated aerodynamics of the time-averaged, steady, rigid case. A reasonable criterion for agreement would be when any two studies show overlap of their aerodynamic coefficients within experimental or numerical uncertainty. There is one such data point in this paper, and a technically sound baseline would require more.

## Authors' contributions

The manuscript was jointly written by G.R.S. and J.T. L.S. ran, and wrote about, the RANS computations. J.T. performed all the wind tunnel experiments, and extracted and analysed the data.

## Competing interests

We declare we have no competing interests.

## Funding

We are grateful to the Air Force Office of Scientific Research for equipment funding and partial funding of J.T. during this work under grants FA9550-15-1-0255 and FA9550-16-1-0392, both under the management of Doug Smith. L.S.'s work is supported by South African National Aerospace Center, and by a Research Completion Grant from the University of Pretoria, South Africa.

## Acknowledgement

We thank Tyler Davis for invaluable assistance in wind tunnel experiments and numerous conversations on this work.

## Footnotes

One contribution of 19 to a theme issue ‘Coevolving advances in animal flight and aerial robotics’.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.