Abstract
The aerodynamic characteristics around a micro air vehicle wing with an inverse-Zimmerman configuration are numerically investigated by an in-house programmed solver particularly dedicated for aircrafts operating in low Reynolds number regime. The complex three-dimensional aerodynamic performance was investigated in terms of force generation and flow structures visualization. Results show that the flow around the low aspect ratio MAV wing is characterized by complex three-dimensional separation-dominated flow. The flow fields exhibit separation, reattachment, secondary separation, secondary reattachment, and strong interaction between the separated boundary layer and wingtip vortices. In addition, the effect of tip-attached vertical stabilizers on flow structure and aerodynamic forces is addressed in this paper. The stabilizers significantly influence both the flow structure and aerodynamic forces via reducing the strength of wingtip vortices and shedding and interacting of wingtip vortices. Eventually, the unsteadiness of the aerodynamics revealed that higher angle of attack will result in stronger unsteady phenomena as demonstrated by the oscillating forces.
Keywords
Introduction
Micro air vehicles (MAVs), defined as a class of unmanned aircraft with a maximal dimension of 15 cm and a flight speed less than 15 m/s, attract an increasing attention due to their potential military and civilian applications. In terms of configuration, existing MAVs can be categorized into three concepts: i.e. fixed-, rotary-, and flapping-wing MAVs. The fixed-wing concept has received substantial research interest in the past years, as they can be directly benefited from well-established conventional aircraft design theories and technologies. However, the inherent small size and low flight speed make MAVs operate in a low Reynolds number regime (typically 104–105) which is much lower than the conventional aircrafts, 1 whereas the low Reynolds number aerodynamics are still not fully explored.
One unfavorable flow behavior in low Reynolds number airfoil aerodynamics is the presence of a laminar separation bubble 2 on the suction surface of the airfoil, which can bring more drag and hence decrease the lift-to-drag ratio. In recent decades, extensive two-dimensional (2D) experimental and numerical studies3–5 have revealed many details of the low Reynolds number airfoil aerodynamics, such as: lift/drag characteristics; static surface pressure distribution; transient behavior of separated laminar boundary layers; and formation, evolution and burst processes of laminar separation bubbles. However, in the three-dimensional (3D) case, as low aspect ratio (LAR) wings are generally employed for MAVs to meet the operational and weight requirement, the wingtip vortex occupies a large proportion on the wing suction side. The wingtip vortex and its interaction with boundary layer separation may induce strong 3D flow phenomena which do not exist on the 2D airfoils. It therefore stimulates further studies with regard to the aforementioned flow phenomena of LAR.
In recent years, some attention has been paid to examine the low Reynolds number flow characteristics of LAR wings. Experimental measurements on the lift and drag characteristics or pressure distributions of LAR wings at low Reynolds number are well documented,6–9 which can provide promising data regarding the various aspects of the effect of wing geometry, camber, aspect ratio, and Reynolds number on the performances of LAR wings. Pelletier and Mueller 6 experimentally investigated the lift, drag, and pitching moment characteristics on a series of thin and cambered LAR wings at low Reynolds number. They concluded that the cambered plates can offer better aerodynamic performance. Later on, Torres and Mueller 7 further studied the effects of wing planform and aspect ratio on the aerodynamics of LAR wings at low Reynolds numbers (range from 7 × 103 to 2 × 105). Their results indicated that LAR can enhance the nonlinearity of lift-curve. The wing-tip vortex can significantly influence the flow over the wings with AR below 1.5 and lift is provided by vortex flow and attached flow on such LAR wings. In views of experimental studies, flow visualization using particle image velocimetry (PIV) of LAR wings have been conducted. Gursul and Taylor 10 performed a PIV measurement in the water tunnel to visualize the flow for several LAR wings with different planforms and observed strong unsteady interactions of vortical structure and shear layers. Tang and Zhu 11 visualized the flow around an accelerating elliptical wing. Kaplan and Altman 12 measured the trailing vortex of LAR wings at Reynolds numbers of 8 × 103 and 2.4 × 104. Very recent PIV studies on the flow structures of LAR wing at low Reynolds numbers can be found in the work of Khambatta and Ukeiley 13 and Gamble and Reeder. 14 Numerically, 3D simulations for MAVs with LAR wings have been performed tremendously.8,15–19 The research group of Florida University also numerically studied the low Reynolds number flow over a fixed membrane wing in a series of work20–23 and unveiled a few aspects of flow characteristics of such flexible membrane wings.
To extend the previous studies, low Reynolds number flows over a LAR MAV wing are studied by numerical simulations in the present paper. The lift/drag performance, flow structure and the interaction between boundary layer and wingtip vortex of this MAV wing are investigated to better understand the aerodynamics of LAR wings at low Reynolds number. Additionally, the effect of vertical stabilizer on the whole aerodynamic performance and the unsteady phenomena at a high angle of attack are addressed in this paper.
Computational model and methodology
MAV model
The fixed-wing MAV used in the present study is shown in Figure 1, which was designed and built at Nanjing University of Aeronautics and Astronautics (NUAA). The MAV uses a similar inverse-Zimmerman LAR wing with two elevators at the trailing edge for the sake of controlling the pitch and yaw. To reduce the rolling maneuverability, two vertical stabilizers were attached to the wingtip below the wing as seen in Figure 1. The total wingspan and the root chord length are both set at 0.15 m. The cross-section of the wing is Eppler-636 airfoil, which was chosen because of its good performance at low Reynolds regime and the thickness allows equipping the onboard electrical and propulsive devices.
The fixed-wing MAV and its schematic views.
Numerical method
The computations in the present study were implemented using an in-house developed finite volume based flow solver. The flow solver solves the compressible Reynolds-averaged Navier–Stoke (RANS) equations, and a low Mach number preconditioning technique is applied to extent its application in low-speed regime by avoiding stiff problem caused the large disparity between the acoustic and convective wave speeds.
24
The compressible Navier–Stokes equations can be written in an integral formulation as
The basic idea of low Mach number preconditioning is to modify the governing equations without altering the steady solution, so that the eigenvalues are all the same order. With preconditioned time derivative introduced, the governing equations of equation (1) become
Note that, when Ma∞→1, the preconditioned equations will be degraded into a non-preconditioned one. The preconditioned Navier–Stokes equations are discretized by a second-order upwind finite-volume scheme in spatial domain and solved by an implicit matrix-free LU-SGS iteration on unstructured hybrid meshes. For unsteady computation, a dual-time stepping scheme is employed to obtain meaningful solution for unsteady flow. A detail description of the flow solver can be found in Xiao et al. 24
Concerning the turbulent flow, the one-equation Spalart-Allmaras (S-A) turbulence model is selected to close the equation system. The S-A turbulence model including the damping functions is a low Reynolds number type turbulence model which resolves the viscous-affected region using the damping functions to predict viscous sublayer behavior. Also, the selection of the S-A turbulence model is supported by the later work of Cosyn and Vierendeels, 25 which evidenced that the S-A turbulence model has a good performance in predicting lift and drag level for low Reynolds number flow simulations.
However, as summarized by Langtry and Menter, 26 the low Reynolds number turbulence models cannot be trusted to predict transition accurately. So we make use of the e n -method which is a pragmatic prediction tool for predicting laminar-turbulent transition. The e n model modified by Lian and Shyy 27 is employed and coupled with the Navier–Stokes code. Having computed the velocity field using the Navier–Stokes solver, the boundary layer integral parameters, such as momentum thickness, displacement thickness, and boundary-layer edge velocity, can be extracted to determine the point of transition. We conduct the implement procedure of e n -method in conjunction with RANS code following the work of Nebel et al. 28 In this paper, the value of n = 9 is chosen for all the computations.
Verification and validation
Grid convergence study
The computational domain of the MAV wing is set at 20 L × 20 L × 20 L cubic where L is the mean chord length. Note that only half of the MAV is employed since it is a symmetric MAV along with the midline of the fuselage. Three unstructured meshes around the MAV wing with different spatial resolutions were generated for the grid convergence study (see Figure 2). For each grid, a triangular surface mesh is constructed first around the CAD model of the wing. Then, starting from the surface mesh, 20 layers of stretched triangular prism elements with a target y+ value of 1.0 are constructed to ensure an appropriate resolution of the viscous sublayer, and the remainder of computational domain is filled with tetrahedrons. The refinement of grids on the leading and trail edges and near the wingtip can be observed in Figure 2. The detail information about the three grids is summarized in Table 1.
Unstructured viscous grids around the MAV wing: (a) coarse grid; (b) medium grid; and (c) fine grid. Information about the unstructured grids around the MAV wing.
The primary objective of the grid convergence study is to estimate the ordered discretization error. The Richardson’s extrapolation
29
is used here to evaluate the computed data by plotting them as a function of N−2/3, where N is the total number of grid cells. The 2/3 power is based on the order of accuracy (second order) of the numerical method used to compute the results. The simulation was computed on the coarse, medium and fine grids with an incident angle at 8°, Re = 1.15 × 105. The solution convergence histories of drag coefficient are shown in Figure 3. The three grids have rendered good convergence rate by which the drag has settled at a constant value after about 4000, 5000, and 5500 iterations, respectively. Figure 4 shows the computed drag coefficient against N−2/3. The linear variation in drag with N−2/3 indicates that the results computed by the current code vary monotonically with grid refinement and asymptotic grid convergence is achieved.
Solution convergence history of drag coefficient (α = 8°, Re = 1.15 × 105). Grid convergence of drag versus number of grid cells.

Validation
The flow solver is validated by comparing the measurement data from Sathaye andYuan, 30 where they experimentally studied the aerodynamic performance of LAR NACA0012 wings at low Reynolds number. The computational grid around the NACA0012 wing is generated according to the medium grid density used in the grid convergence study in the earlier section. The case of Re = 8.4 × 104 and α = 6° was calculated.
Figure 5 shows the computed surface pressure distribution at the two span locations of y/(b/2) = 0.6875, 0.8125, where y is spanwise position and b is full wing span, in comparison with the measured data of Sathaye.
30
As can be seen, a good agreement is observed between the computed results and the experimental data, indicating that the CFD code used in this study can accurately predict the low Reynolds number flow in the MAV operating region. Also, this test case and the foregoing grid convergence study depicted that the medium grid resolution is sufficient for accurate simulations of these type flows.
Comparison of pressure distribution on the NACA0012 wing.
Results and discussion
Flow separation and its interaction with wingtip vortices
For a low aspect ratio and highly swept back wing as used in the present study, flow separation behavior is more complex than on a conventional airplane. This section will focus on the flow separation and its interaction with wingtip vortices. To obtain a better understanding of the flow separation phenomena, the topological approach is employed. The near-surface streamline patterns at several angles of attack from 0° to 18°, with an increment at 2°, are presented in Figure 6 which clearly demonstrates the complex 3D flows on the upper surface such as separated flow, reattached flow, and span-wise flow.
Near-surface streamline patterns varied with angles of attack from 0° to 18° (without stabilizers, Re = 1.15 × 105).
Although the streamline patterns at different angles of attack are not absolutely identical, the basic flow structures are similar as can be seen from Figure 6. Here, a pattern at one of the angles of attack (α = 8°) is selected for characterizing the fluid dynamics of the current LAR wing as illustrated in Figure 7. It can be observed that four confluence lines appear on the upper surface (see Figure 7(a)) with a pattern of twice separation and twice reattachment. Figure 7(b) indicates the pressure coefficient distributions on different span-wise sections.
Illustrations of flow structure and pressure distribution due to boundary layer separation: (a) skin friction streamline and (b) pressure coefficient distribution (α = 8°).
Starting from the leading edge, the flow accelerates along the chordwise on the upper surface and induces a pressure drop as expected. The attached boundary layer separates when undergoing a relative large adverse pressure gradient and the separation points can be revealed from the first streamline confluence line (separation line). When the adverse pressure gradient drops, the separated flow will reattach in a short distance of chord length, which causes the second streamline confluence line (reattachment line) formed on the upper surface. Effected by the sweepback of the LAR wing and absorbed by the wingtip vortex, the separation region between the first and the second confluence lines moves downstream along a span-wise direction from the root to the wingtip and finally merges into the wingtip vortex. As an indication of the separated flow, a plateau in the surface pressure distribution as seen in Figure 7(b) exists on the upper surface for all five span-wise locations.
Then, the flow on the upper surface experiences a second separation and reattachment indicated by the third and the fourth streamline confluence lines (see Figure 7(a)). Two types of flow are observed in the region between the third and the fourth confluence lines on the upper surface, i.e. separated flow in the root area and attached span-wise flow in the outer portion. The secondary separation region exists near the root characterized by the secondary plateau in pressure distribution of y/b = 0.01, 0.1, 0.2, 0.3 as presented in Figure 7(b). The plateau was found to have a similar width near the root (y/b = 0.01, 0.1, 0.2) but then reduces in size outward (y/b = 0.3). At the span-wise position of y/b = 0.4, this plateau completely vanishes indicating the disappearance of the secondary separation region. The reason is that, in the outer portion of the upper surface, the separated flow is absorbed and energized by the tip vortices, and gradually becomes attached span-wise flow.
The flow structures with different angles of attack are quite similar except some slight variations as observed in Figure 6. Increasing the angle of attack, the first streamline confluence line moves forward to the leading-edge due to a steeper adverse pressure gradient, the secondary separation region reduces in size and the attached span-wise flow is enhanced because of the stronger tip vortices.
Considering the influence of Reynolds number on the flow structure of the current LAR wing, the computations were conducted under different Reynolds number at angle of attack of 8°. Figure 8 plots the computed near-surface streamline patterns at Reynolds number between 5.75 × 104 and 1.70 × 105. As can be seen, higher Reynolds number will delay the first separation on the wing surface and the region between second and third confluence lines are narrowed with an increasing Reynolds number. Overall, the flow structures at different Reynolds numbers appear to have a similar flow topology.
Near-surface streamline patterns at different Reynolds numbers (without stabilizers, α = 8°): (a) Re = 5.75 × 104; (b) Re = 8.63 × 104; (c) Re = 1.15 × 105; and (d) Re = 1.70 × 105.
Effect of wing platform
As mentioned above, the flow topology for a LAR swept back wing experiences a complex flow behavior, i.e. twice separation and twice reattachment. A question was raised up that if such flow behavior is unique for the specific MAV in our study. To clarify the aforementioned question, another fixed wing MAV with a trapezoidal wing platform as shown in Figure 9 (hereafter referred to as “Trap-MAV”) is analyzed and compared with the inverse-Zimmerman platform MAV. Note that the two MAVs have the same chord and span length.
Figures 10 and 11 show the computed near-surface streamline patterns on the upper surface for the Trap-MAV at different angles of attack and different Reynolds numbers, respectively. It is clear that the flows over the upper surface of the Trap-MAV wing also experience separation, reattachment, secondary separation, and secondary reattachment, similarly to the above investigated phenomena, which is not different as the layout or platform changes.
Fixed-Wing MAV with a trapezoidal wing platform. Near-surface streamline patterns at different angles of attack (Trap-MAV without stabilizers, Re = 1.15 × 105). Near-surface streamline patterns at different Reynolds numbers (Trap-MAV without stabilizers, α = 8°): (a) Re = 5.75 × 104; (b) Re = 8.63 × 104; (c) Re = 1.15 × 105; and (d) Re = 1.70 × 105.


Effect of vertical stabilizer
The vertical stabilizers are placed at both wingtips to guarantee a stable flight by reducing the rolling effect. Moreover, the presence of such stabilizers will also significantly change the flow phenomena on the upper surface. The current section would emphasize the effect of the stabilizer on the aerodynamic performance of the fixed wing MAV. Besides the numerical simulation, a wind tunnel experiment of the 1:1 scale MAV was also conducted at the NUAA Aerodynamic Lab, where the aerodynamic forces are used to validate computational results. The experimental setup is shown in Figure 12. Note that the propeller was removed to avoid its rotating effect. The close-loop low-speed wind tunnel has a 3 m × 2.5 m × 6 m test section and the turbulent intensity is around 0.1%–0.14% within the wind speed range of 5–90 m/s. The MAV model was directly mounted on a miniature six-component force transducer to instantaneously record the forces and moments. The measurement uncertainty for the force sensor is less than 1% with a 95% confidence level. All the data were acquired with a recording frequency at 20 kHz for averaging.
The MAV in the wind tunnel.
The computed coefficients of lift, drag, and pitching moment (moment center is located at the quarter-chord of the airfoil at the symmetrical plane) are presented in Figure 13 as well as the experimental data. A good agreement was revealed from the experimental and numerical results before stall (α = 28°) which indicate that the numerical method developed in the present study can accurately predict pre-stall lift, drag, and pitching moment. The numerical model slightly over-predicts the stalling angle and maximum lift coefficient. Mild nonlinearities can be seen in the lift curve, presumably due to a growth in the low pressure cells at the wing tips of LAR wing,7,13 Negative slopes CMα at all angles of attack are observed in the pitching moment curve. This would imply that the present wing is statically stable. The lift-curve slopes CL
α
of the wing model with stabilizers obtained by wind-tunnel test and by computed data are 0.0332/° and 0.0329/°, respectively. Here we compare the CL
α
of the present wing with theoretical values estimated by a theoretical formula for thin wings of different aspect ratios,
6
Comparison of aerodynamic forces between the two models with/without stabilizers (Re = 1.15 × 105).

The computational results with and without stabilizers almost collapse on top of each other at a very small angle of attack (say α < 6°) where the tip vortex is also small. A visible disparity was found when the incident angle was between 6° and 25°, where the lift of the MAV with stabilizer is higher than that of without stabilizers. To reveal the effect of the stabilizer when the tip vortex dominates the flow structures on the upper wing, the velocity vectors on three chord-wise oriented slices (x/L = 0.8, 0.9, 1.0) when α = 14° are compared in Figure 14. As can be seen, for the model without stabilizers, the tip vortex is relatively strong and the circulatory motion of the flow over the wingtip causes obvious span-wise flow below the wing near the wingtip. However, for the model with stabilizers, by comparison, the stabilizer basically deters the flow from the high-pressure region below the wing surface to reach the low-pressure region above the wing and no clearly span-wise flow occurs below the wing, which results in a weaker tip vortex compared with the model with stabilizer.
Comparison of computed velocity vector field at different chord-wise oriented slices between the two models without (top row)/with (bottom row) stabilizers (Re = 1.15 × 105, α = 14°).
Respectively, Figures 15 and 16 show the comparisons of near-surface streamline pattern and surface pressure coefficients contours between the two numerical configurations, i.e. without/with stabilizers. The presence of the vertical stabilizers does not influence the upstream where the first separation lines stay at the same location on the upper surface as depicted in Figure 15(a). A small pressure drop is found at the rear portion of wingtip (Figure 16(a)). On the lower surface, for the wing without stabilizers, obvious span-wise flow can be seen near the wingtip (Figure 15(b)). However, for the wing with stabilizers, they are cleared away. Due to the deterring of stabilizers, more fluid momentum below the wing is transferred into pressure instead of being shed as vortex at the wingtip, resulting in higher pressure in the hind portion of lower surface as indicated by Figure 16(b); hence, a larger lift and pitching down moment was created as shown in Figure 13.
Comparisons of near-surface streamline patterns between the two models with/without stabilizers (Re = 1.15 × 105, α = 14°): (a) upper surface; and (b) lower surface. Comparisons of surface pressure contours between the two models with/without stabilizers (Re = 1.15 × 105, α = 14°): (a) upper surface; and (b) lower surface.

The reducing of pressure drop on the upper surface and the strengthening of high pressure below the wing, caused by the stabilizers, can also be quantificationally demonstrated by Figure 17 which gives the comparison of pressure coefficient between the two models at the chord-wise location of x/L = 0.8 at angles of attack of 6°, 8°, 10°, 12°.
Comparison of pressure coefficient distributions between the two models with/without stabilizers at angle of attack of 6°, 8°, 10°, 12° (Re = 1.15 × 105, α = 14°, x/L = 0.8): (a) upper surface; and (b) lower surface.
The unsteadiness aerodynamics of the MAV
As discussed above, the low Reynolds number flow around the LAR wing is very complex in terms of the 3D flow structures and the interaction between tip vortices and separation flow. The shedding and interacting of these vortices will be expected to cause unsteadiness in aerodynamic performance especially at high angles of attack. A straightforward behavior of the unsteady phenomena is the oscillation of aerodynamic forces.
Figure 18(a) shows the steady-state solution convergence history of lift coefficient versus iteration number at several different angles of attack. As can be seen, after around 2000 iterations, the lift coefficients for smaller incident angles experience smaller oscillating amplitude when compared with the higher angle case, which evidences a higher unsteadiness for higher angle of attack. For the purpose of comparison, unsteady simulations were performed at the angles of attack of 16°, 22°, and 28°. The time histories of the lift coefficients plot in Figure 18(b) indicate the unsteadiness of flows at high angles of attack.
Comparisons of lift coefficient convergence histories between steady and unsteady computations: (a) Steady computations; and (b) unsteady computations.
Mean square roots of lift coefficients by steady and unsteady computations.
The previous study of Cummings and Morton 31 in which the Reynolds number is higher than that of the MAV regime reported that the unsteady computations predicated noticeably lower lift coefficients than did the steady computations at large angles of attack. However, cases change under MAV flight conditions. In the present study, the difference of average lift between the steady and unsteady computations is quite small, as can be seen in Table 2.
Conclusions
The low Reynolds number aerodynamics of a LAR MAV wing has been numerically investigated. Compared with experimental data, the computations can accurately predict pre-stall aerodynamic forces. The complex 3D flow-field around the wing exhibits a variety of separation, reattachment, secondary separation, secondary reattachment, and strong interaction between the separated flow and wingtip vortices. The stabilizers attached to the wingtip influence mildly both the flow structure and aerodynamic forces by reducing the strength of wingtip vortices. At high angles of attack, the shedding and interacting of wingtip vortices and separated flows cause aerodynamic unsteadiness of which the straightforward behavior is the oscillating forces computed by both steady and unsteady simulations.
Footnotes
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work has been supported by the National Science Foundation of China (grant no. 11002072) and the China Scholarship Council (grant no. 201303070173 and 2011683005).
