Computational investigation of insect-based flapping wings for micro air vehicle applications

Abstract

In this study, a computational fluid dynamics analysis was performed on bio-inspired micro aerial vehicle scale rigid flapping wings. The computational fluid dynamics analysis used a compressible unsteady Reynolds Averaged Navier–Stokes solver with low-Mach number preconditioning to study the complex, highly vortical, three-dimensional flow of low aspect ratio flapping wings at micro aerial vehicle-scale Reynolds numbers. The wing was flapped at a constant 5 Hz flap frequency at a mean chord reference Reynolds number of 25,000. The flapping and pitching kinematics were set to match those of a previous experimental study resulting in a constant flap stroke of 107° at translational pitching angles of 40°, 50°, and 60°. The force and flowfield measurements of the previous flapping-wing experiment were used for the validation of the 3D computational fluid dynamics model. The objectives of this effort were to understand the unsteady aerodynamic mechanisms and their relation to force production and aerodynamic efficiency on a rigid wing undergoing an insect-type flapping motion with passive pitching kinematics. Overall, the computational fluid dynamics results showed good agreement with the measured experimental force data. Additionally, the computational fluid dynamics simulation was able to adequately predict the process of leading edge vortex formation and shedding observed during experimentation. A vorticity summation approach used to calculate the strength of the leading edge vortex from the experimental measurements and from the computational fluid dynamics predicted flowfields showed comparable results. The computational fluid dynamics results were utilized to further analyze the differences in the flowfield and leading edge vortex formation for the three pitch angles tested as well as the instantaneous aerodynamic loads and aerodynamic power.

Introduction

As interest and technological potential for the widespread use of micro air vehicles (MAVs) continue to grow, it has become increasingly imperative to assess the performance and efficiency of a given MAV design. Current MAV systems based on conventional aircraft designs, such as fixed wing and rotary wing designs, suffer from decreased lifting force generation and degraded aerodynamic efficiency.^1,2 The major factor for the lowered aerodynamic performance can be attributed to low-operating Reynolds numbers effects (O∼10⁴), which many studies have found to cause degraded lift-to-drag ratios, increased drag coefficient magnitudes, and increased susceptibility to flow transition and separation. ³ Attention is focused on flapping wing-based MAV designs because of the potential for increased lift capability, aerodynamic performance, gust tolerance, and maneuverability exhibited by biological flapping-wing flyers, which operate in the same aerodynamic regime. The ability of flapping-wing designs to generate large forces is achieved by utilizing unsteady aerodynamic mechanisms such as the leading edge vortex (LEV). ⁴ Other unsteady flight mechanisms, such as rotational lift and wake capture, have been postulated to aid in the enhanced aerodynamic capabilities of flapping-wing flight; ⁵ however, these phenomena only happen for a short duration of the wing stroke and depend heavily on the timing of wing rotation. Therefore, these mechanisms do not contribute much to the time-averaged wing forces. The LEV is thought to be the most prominent feature responsible for aerodynamic force augmentation,^6,7 but there is still limited understanding of the complex, unsteady vortical flow inherent to flapping wings and how the LEV influences aerodynamic force production over time.

Some initial studies sought to analyze the vortex formation on wings undergoing pitching and plunging kinematics using 2D analysis techniques.^8–10 While these studies did play an important role in building the current understanding of how unsteady mechanisms influence the production of aerodynamic force on pitching and plunging airfoils, their analysis focused on forward flight as opposed to the hover flight condition. Yuan et al. ¹¹ conducted experimental investigations on a 2D pitching and plunging airfoil in hover in accord with 2D numerical analysis. Results showed that shed vortical structures during flapping may remain near the airfoil surface and can potentially interfere with new vortex formation affecting overall aerodynamic force production. While 2D analysis can provide some significant insight into the unsteady fluid dynamics of the flapping-wing problem, the wings of many biological flyers and man-made flapping MAVs have a low aspect ratio (AR < 5). Thus, it is expected that 3D effects could have a significant influence on the aerodynamics about the wing and should be appropriately accounted for within analysis.

Many previous 3D experiments were performed on rigid, dynamically scaled mechanical models in oil or water tanks.^12–14 Conducting the experiments in oil or water allows for better matching of the Reynolds numbers associated with natural flyers and MAVs while operating at much lower frequencies. However, using a lowered operating frequency decreases the generated inertial force and its influence on the total force measured during experimentation. The prior 3D experimental studies have aided tremendously in gaining a better understanding of the vortical structures produced on flapping wings and how various unsteady mechanisms and flap kinematics affect aerodynamic force production. However, the aforementioned 3D experiments were limited to using rigid wings at very low Reynolds numbers (O∼10³) (which is below that for many currently feasible and practical MAV designs) and lacked accompanying computational studies.

There have been several notable studies that include 3D simulations with experiments. Liu and Kawachi ¹⁵ published the first 3D Navier–Stokes simulation of insect flapping wings. Their simulation modeled the vortex structure over a scaled mechanical flapping hawkmoth ¹⁶ and showed good qualitative agreement with smoke flow visualization. They also examined the aerodynamic characteristics of the hawkmoth and fruit fly wing shapes, focusing on the spanwise flow inside the core of the LEV. ¹⁷ Other notable examples include studies by Sun and Tang ¹⁸ and Ramamurti and Sandberg. ¹⁹ Particularly, the Robofly mechanical flapping apparatus has provided a wealth of force data to aid in the validation of 3D computational fluid dynamics (CFD) studies. However, there has not been much comparison of the flowfield. More recently, studies by Badrya et al. ²⁰ investigated the influence of Reynolds number, wing kinematics, and AR on the aerodynamic force and power of a modeled blow fly wing (similar in shape to the Robofly wing). Their simulation results showed that Reynolds number and wing pitch angle have a strong effect on peak aerodynamic force production and aerodynamic power. Despite all the insight gained from these studies, it is still not clear whether the same unsteady force augmentation mechanisms and conclusions on flow physics recognized at these low Reynolds numbers apply at higher, MAV-scale Reynolds numbers.

Many notable flapping-wing studies at the University of Maryland have been conducted in air on a biomimetic flapping-wing device.^21–23 The flapping mechanism was operated at a mean chord Reynolds number of approximately 15,500. While these studies aided to further the understanding of flapping wings in hover, the number of tests was limited and an accompanying 3D numerical analysis was not conducted. Recently, extensive flowfield studies, in tandem with a high-fidelity aeroelastic analysis, were conducted on a flapping-wing test setup undergoing avian-like flapping kinematics.^24,25 While the computational results correlated well with the measured experimental data, the tests focused on a flapping wing at a fixed root pitch angle in forward flight as opposed to hover.

Recently, Benedict et al. ²⁶ performed instantaneous force measurements and flowfield studies on a rigid wing undergoing active insect-like flapping kinematics with passive pitch rotations. Given that one of the main objectives of the study was to validate the applicability of an inertial force subtraction method in air, only a limited set of flowfield experiments were performed. An important area for expansion of this study was the inclusion of computational results to further analyze the flow. In addition, a high-fidelity computational study would provide added confidence in the feasibility of using an inertial force subtraction technique to aid in analyzing the aerodynamic forces in air, especially when the inertial forces account for a majority of the total force measured.

The present work evaluates the aerodynamic properties and performance of a hover-capable flapping-wing MAV concept with a high-fidelity unsteady CFD model using previous force and flowfield experiments for validation. The predictive capability of the CFD model is assessed against instantaneous aerodynamic force time-history and particle image velocimetry (PIV) flowfield measurements. Focus is on resolving the temporal and spatial development of the LEV and the effects of high angles of attack on LEV evolution over time. By further understanding the complex and unsteady 3D flow environment and overall flapping MAV aerodynamic efficiency, a better understanding of the influence of the LEV can be obtained.

Computational methodology and validation experiment description

Description of validation experiment

To validate the results from the numerical simulation, the force and flowfield data from the experiment conducted on a rigid, bio-inspired low AR wing by Benedict et al. ²⁶ were utilized for comparison. The wing used throughout the experiment, depicted in Figure 1, consists of a stiff carbon fiber frame covered by a thin Mylar membrane film resulting in a total wing weight of 3 g. Wing length (b) and thickness (h) were 15.24 cm and 1.6 mm, respectively. The wing had a straight leading edge from root to tip with a tapered trailing edge defined by the equation

$c(x)=a_1{r}^3+a_2{r}^2+a_{3}r+a_4$

(1)

where a₁ = −0.5764, a₂ = 0.1503, a₃ = −0.9443, a₄ = 9.1091, c is the wing chord and r is the spanwise location from the wing root in centimeters. Wing planform area (S) is approximately 97.85 cm² giving the wing an AR of 2.38 with a mean aerodynamic chord ( $\overline{c}$ ) of 6.42 cm. Note that there is an 8.89 cm offset between the flapping axis of rotation and the wing root and a 1.016 cm offset from the leading edge to the pitching axis of rotation. All of the experimental tests were carried out at an operating flap frequency of 5 Hz and a flapping stroke amplitude of 107°. The mean chord Reynolds number (Re) for all tests was approximately 25,000 based on an average wing maximum tip speed (V_tip) of 6.3 m/s. Tests were conducted at three separate translational pitch angles of 40°, 50°, and 60°. Figure 2(a) and (b) shows the variation of the flap and pitch kinematics. Flowfield measurements were taken using time-resolved, two-component PIV with a double-pulsed Nd-YLF laser (Litron LDY304, 30 mJ/pulse, 10 kHz max). The laser sheet was set to illuminate the region of interest at the 25%, 50%, and 75% span locations, measured from the wing root, at various azimuthal positions within the flap cycle. The estimated uncertainty of the measured velocity within the experimental flowfield measurements is 0.667%. A more detailed description of the experimental set-up and procedure can be found in Benedict et al. ²⁶

Figure 1.

Wing pitching mechanism. (a) Schematic of flapping wing with coordinate system and (b) passive wing pitching mechanism.

Figure 2.

Variation of pitch angle for the 40°, 50°, and 60° translational pitch cases. (a) Pitch angle variation in vacuum and (b) pitch angle variation in air.

CFDs solver

The 3D simulations of the flapping wing were computed using an in-house overset transonic unsteady Reynolds Averaged Navier–Stokes (URANS) solver called OVERTURNS. ²⁷ To solve the URANS equations, the solver utilizes the diagonal form of the implicit approximate factorization method developed by Pulliam and Chaussee with second-order accuracy in time. A time accurate, dual-time stepping scheme, in conjunction with low-Mach pre-conditioning, is used to improve simulation convergence and solution accuracy. A third-order accurate MUSCL scheme utilizing a Korens limiter is used for spatial reconstruction. The inviscid and viscous terms are respectively computed using a Roe flux difference splitting scheme and a second-order central differencing scheme. The Spalart-Allmaras turbulence model with rotational correction is employed for closure of the URANS equations and to more accurately account for the formation of vortical flow.^28,29

Computational model

An O-O type mesh, displayed in Figure 3, was utilized with 231 × 161 × 128 nodes in the wrap-around, spanwise, and normal directions, respectively. The grid size was chosen based on the conclusions of a grid sensitivity study using the nominal grid, a coarser (155 × 107 × 82) grid, and a finer (315 × 187 × 143) grid. Simulation results showed that there was less than 3.0% difference between the predicted mean lift and drag values using the nominal and finer grids. The wing planform of the CFD mesh was created to match the experimental test wing described earlier. In the simulation, the wing has a thickness-to-chord ratio (h/ $\overline{c}$ ) of 2.5% and was offset from the flapping axis and pitching axis by 1.38 $\overline{c}$ and 0.16 $\overline{c}$ , respectively, to match that of the experimental set-up. A higher concentration of grid points were placed near the leading edge, trailing edge, and wing tip to better capture the highly vortical flow expected near those regions.

Figure 3.

O-O wing mesh.

To capture the highly vortical flow in the wake, a Cartesian background mesh was implemented with 252 × 252 × 92 nodes in the x, y, and z directions, respectively. The background mesh was made up of approximately six million nodes with farfield boundary conditions applied 10 chord lengths away from the wing surface. All simulations were run for at least four consecutive flap cycles using 2880 iterations per flap cycle with eight sub-iterations to minimize factorization errors and improve the unsteady computational solution accuracy. The flap kinematics were set to vary harmonically, as in the experiment, with the flap amplitude set to ± 53.5°. Given the nature of the experiment, the measured pitch kinematics in air, shown in Figure 2(b), were read directly into the simulation in order to correctly account for the pitch angle variation over the flap cycle for the nominal 40°, 50°, and 60° pitch cases studied.

Comparison of experimental and computational results

Force-time history validation

This section focuses on the comparison of the aerodynamic force-time histories and flowfields acquired from the experimental tests and CFD simulations. Primarily, the experimental force measurements and PIV flowfield data are used to validate the predictive capability of the unsteady CFD simulations. Figure 4(a) and (b) is a comparison of the measured and predicted aerodynamic (lift and drag) force-time histories acquired from the experiment and CFD analysis for the 50° translational pitch case. The experimental force trends shown consist of an average of 10 trials where each trial is an average of 90 flap cycles. The shaded region encompassing the experimental force trend represents the standard deviation from the average of the 10 trials. The average is shown as a black line. Figure 4(a) shows the variation of the lift over the course of one flap cycle. The lifting force remains positive throughout a majority of the forward and backward strokes, but some negative lift was produced at stroke reversal (t/T = 0.0, 0.5, and 1.0). This behavior in the lift force trend, described in literature,^30,31 is expected because of the delayed rotation kinematics inherent to the passive pitching mechanism used. Nevertheless, the overall variation in lift predicted by the CFD model agrees well with the measured force trend. There is a slight phase difference between the experimental and computational force-time histories, but the predicted instantaneous force magnitudes are within the bounds of uncertainty for the experimental results.

Figure 4.

Instantaneous variation of experimental and computational aerodynamic forces. (a) Lift force variation (50° case) and (b) drag force variation (50° case).

Figure 4(b) displays the variation in drag over the course of one flap cycle. Note that the sign value of the instantaneous drag designates the direction in which it acts. Positive and negative values of drag indicate that the drag force acted in the positive and negative x-directions, respectively. Throughout the flap cycle, drag opposed the wing flapping motion. During the flap cycle, there was a slight difference between the measured and predicted drag force values in the instance at which maximum instantaneous drag occurred. The peak magnitudes of the CFD drag force-time history occur at mid-forward stroke and mid-backward stroke (t/T = 0.25 and 0.75). Peak magnitude in the experimentally measured drag values occurred slightly after and slightly before mid-forward stroke and mid-backward stroke, respectively. However, it is important to note that the standard deviation in the experimental data becomes quite large near the peak values of drag in comparison to that near the end of a stroke.

The large deviation in the drag force peak magnitude is due to the challenge of separating the large inertial loads in the flap plane from the aerodynamic loads generated by the wing. Nevertheless, the predicted drag force-time history from the CFD lies within the experimental bounds of uncertainty. Overall, the variation in drag predicted by the CFD model correlates well with the measured force trend. This gives confidence in the accuracy and feasibility of using vacuum force subtraction techniques to experimentally measure instantaneous aerodynamic forces in air, especially when the inertial forces dominate the total force measurements.

Flowfield comparison: Spatial variation validation

In order to compare the PIV and CFD flowfields, the vorticity field was used, where vorticity ( $\omega$ ) is equal to $\nabla \times V$ . The velocity gradients from the velocity field were calculated using a second-order accurate “Least Squares” differencing method. Note that overlaid velocity vectors are plotted within the vorticity contour plots with every other velocity vector displayed for clarity. In addition, the magnitude of vorticity is normalized by the maximum tip velocity experienced by the wing (V_tip) and the mean aerodynamic chord ( $\overline{c}$ ). Figure 5 provides a comparison of the measured PIV and computed CFD vorticity field near the leading edge of the wing for the 50° translational pitch case. Figure 5(a) to (c) and (d) to (f) shows the results at the 25%, 50%, and 75% span locations, for the PIV and CFD, respectively, at the mid-forward stroke (t/T = 0.25). This comparison displays the capability of the CFD model to predict the spatial variation of the vorticity within the LEV. The mid-forward stroke position (t/T = 0.25) was selected because at this instance in the flap cycle, the wing experienced the highest rotational velocity and vortex strength was expected to be near its greatest magnitude.

Figure 5.

Spatial variation of PIV and CFD vorticity contours at t/T = 0.25. (a) PIV, r/R = 25%, (b) PIV, r/R = 50%, (c) PIV, r/R = 75%, (d) CFD, r/R = 25%, (e) CFD, r/R = 50%, and (f) CFD, r/R = 75%.

The LEV near the wing root in Figure 5(a) and (d) is tightly formed and attached to the wing. The CFD solver was able to accurately model the LEV size and location in comparison to the PIV results. In Figure 5(b), the LEV has burst, becoming more diffused and begins to separate from the wing surface. Again, the CFD solver is able to model this behavior, shown in Figure 5(e), displaying a larger, more diffused LEV. For the PIV (Figure 5(c)) at the 75% span location, the LEV has decreased in size and is clearly not attached to the wing. The separated LEV remains in relatively close proximity to the wing surface and a smaller secondary LEV has begun to form at the leading edge. While the CFD results at the 75% span location in Figure 5(f) also predict a decrease in LEV size and separation from the wing surface, there are noticeable differences between the PIV and CFD. Mainly, the CFD solver over-predicts the degree by which the LEV has separated from the wing surface. Also, the general formation of the predicted LEV differs from that measured by the PIV. Qualitatively, the LEV from CFD is less diffused and has a more oblong shape in comparison to the LEV from PIV. Largely, there is good correlation between the PIV and CFD results, and many of the vortical flow variations in the LEV were predicted by the CFD model. There is some deviation between the PIV and CFD vorticity fields at the more outboard sections, but the salient flow features are still predicted.

Flowfield comparison: Temporal variation validation

Figure 6(a) to (c) and (d) to (f) shows the results at the t = 0.15T, 0.25T, and 0.35T time instances of the flap cycle, for the PIV and CFD, respectively, at the 50% span location of the wing. The contours in Figure 6 validate the temporal evolution of the LEV predicted by the CFD against that measured during the experiment. In Figure 6(a), the LEV has formed and rotation of the flow is visible via the overlaid velocity vectors. In Figure 6(b), the LEV begins to form, but clear rotation of the flow has not yet begun. Figure 6(b) and (e) show that the LEV has burst, become more diffused and then begins to separate from the wing surface. From the PIV data in Figure 6(c), the flow has become highly separated with a small amount of vorticity generated at the leading edge. Also, there is a small region of counter-rotating vorticity located below the region of separation. The CFD results shown in Figure 6(f) exhibit largely separated flow, but the concentration of vorticity is not as diffused as that seen in Figure 6(c). The CFD model is also predicting a small region of counter-rotating vortical flow similar to the PIV. However, the location of the counter-rotating vorticity is much closer to the leading edge in comparison to that seen in Figure 6(c). Overall, the CFD solver has demonstrated an acceptable capability to model the temporal evolution of the LEV and predict the salient variations of the LEV over time seen in the PIV.

Figure 6.

Temporal variation of the PIV and CFD contours at the 50% span location. (a) PIV, t/T = 0.15, (b) PIV, t/T = 0.25, (c) PIV, t/T = 0.35, (d) CFD, t/T = 0.15, (e) CFD, t/T = 0.25, and (f) CFD, t/T = 0.35.

LEV circulation strength validation

The PIV and CFD vorticity contours in the previous sections provided a qualitative comparison of the LEV from the measured and predicted results. In addition, the velocity field can be used to quantitatively compare the estimated LEV circulation strength in both the CFD and PIV. The circulation is obtained by choosing a suitable integration contour (i.e., around the LEV in this case) and numerically evaluating the closed-loop velocity integral. This loop must completely enclose the LEV, but obviously without intersecting any other extraneous circulation. The equation for circulation is given below as

$\Gamma =-\oint V·{\mathrm{ds}}^{\rightharpoonup }=-\int {\int }_A^{(}\nabla \overrightarrow{V})\text{d}S$

(2)

where $V$ is the velocity vector, ${\mathrm{ds}}^{\rightharpoonup }$ is the directed line segment at a point on a predefined contour, $\nabla \times V$ is the vorticity, and ${\text{d}S}^{\rightharpoonup }$ is the area enclosed by a curve. For this analysis, a vorticity threshold was set (−400 < ω < 400) so that low magnitude extraneous vorticity would not be included in the integration contour that was used to enclose the LEV. Next, the vorticity field is scanned in order to locate pockets of vorticity of a certain size (3 × 3 grid square) within that threshold. These regions are then stored as grid locations. Line contours are drawn around these concentrations of vorticity so that an area integral of vorticity can be assessed to obtain the circulation for the vortices in the PIV experiment flowfield and those generated by the CFD solver.

Figure 7 provides an illustration of how the circulation is computed using a representative phase-averaged vorticity field for the 50° translational pitch case. The white dots represent the calculated center of vorticity for a given region of vorticity while the black boxes represent the discretized vorticity field. By summing the area integral of vorticity within the vorticity field, the circulation strength of the LEV can be determined. The circulation strength of the LEV is found using the following equation

$\Gamma =\sum _{i=1}^n{\Gamma }_i=\sum _{i=1}^n\int {\int }_A^{(}\nabla \overrightarrow{V})\text{d}{S}_i$

(3)

where Γ _i is the incremental circulation within the incremental surface area $\text{d}{S}_i$ and n is the total number of incremental surface areas within the discretized vorticity field.

Figure 7.

Representative illustration of circulation summation method.

Figure 8 shows the circulation values for the shedding process at the 50% span location for the 50° translational pitch case. From the PIV analysis, it can be seen that the circulation of the LEV at the 50% span location increases up to the midstroke location (t/T = 0.25) and then drops. However, for the CFD, the circulation increases past the midstroke location to about t/T = 0.30 and then begins to decrease. Also, the decrease in circulation after midstroke can be correlated to the LEV burst within the vorticity contours (see Figures 5 and 6).

Figure 8.

Comparison of PIV, with error bars, and CFD LEV circulation strength.

Comparison of computational results

Computational aerodynamic force-time histories

In this section, the simulation results from the CFD analysis will be used to further investigate the difference in aerodynamic force production, flow physics, and aerodynamic efficiency for the three cases tested. Note that the flap and pitch kinematics used during simulation for the three translational pitch cases (40°, 50°, and 60°) can be seen in Figure 2(b).

Figure 9(a) shows the variation in instantaneous lift force over the course of one flap cycle for the three translational pitch cases at a constant flap frequency of 5 Hz. While the 50° and 60° translational pitch cases have similar force trends, the 40° case has a noticeably lower maximum peak in lift, especially during the forward stroke. Additionally, the force trend for the 40° case has a distinct difference in phase in comparison to the 50° and 60° cases. For all three cases, the aerodynamic lift produced is positive for the majority of the flap cycle. Negative lift is produced at the moments of stroke reversal (t/T = 0.0, 0.5 and 1.0) due to the delayed rotation of the prescribed pitch kinematics recorded from the experiment.

Figure 9.

Instantaneous force variation over a flap cycle for 40°, 50°, and 60° pitch cases. (a) Aerodynamic lift force variation and (b) aerodynamic drag force variation.

Figure 9(b) displays the variation in drag over a flap cycle for the three translational pitch angle cases. As anticipated, the peak magnitude of the drag force increases with increased translational pitch angle. The peak in drag force for all cases occurs at the mid-forward stroke and mid-backward stroke. This is expected due to the fact that the rate of change of the flap angle ( $\stackrel{·}{\phi }$ ) is greatest at t/T = 0.25 and 0.75, and the drag force is proportional to the square of rotational velocity.

Temporal variation of computational sectional vorticity contours

During experimentation, only a limited amount of PIV flowfield measurements were conducted. However, the predicted flowfield from the CFD simulation can be used to further understand the flow about the flapping wing. Focus will be on analyzing the vorticity near the leading edge for all three translational pitch cases to provide more insight into how the evolution of the LEV varies with pitch angle and how that influences aerodynamic force production.

Figure 10 contains the normalized vorticity contours, with overlaid velocity vectors, at the 50% span location for all of the translational pitch angle cases at various instances of the flap cycle. Figure 10(a), (d), and (g) shows representative normalized vorticity contours for the 40° case at t/T = 0.15, 0.25, and 0.35, respectively. Similarly, Figure 10(b), (e), and (h) and 10(c), (f), and (i) shows the vorticity contours for the 50° and 60° cases, respectively, at the aforementioned time instances. In Figure 10(a) to (c), the LEV has just begun to form and the formation of the vortex at t/T = 0.15 is more coherent with increasing translational pitch angle. This is because at this point in the flap cycle, for all three cases, the wing still is undergoing stroke reversal pitch rotation and is operating at an extremely high instantaneous angle of attack. Note that in Figure 10(a) to (c), the instantaneous pitch angle for the 40° case is greater than the 50° case, which, in turn, is greater than the 60° case. Thus, as the pitch angle decreases toward the set translational pitch value, the flow reattaches to the airfoil and a new LEV is formed.

Figure 10.

Temporal variation of CFD vorticity contours at the 50% span location 40°, 50°, and 60° translation pitch cases. (a) t/T = 0.15, 40° case, (b) t/T = 0.15, 50° case, (c) t/T = 0.15, 60° case, (d) t/T = 0.25, 40° case, (e) t/T = 0.25, 50° case, (f) t/T = 0.25, 60° case, (g) t/T = 0.35, 40° case, (h) t/T = 0.35, 50° case, and (i) t/T = 0.35, 60° case.

In Figure 10(d) to (f), a LEV is clearly seen near the wing surface for all three cases, but there are significant differences between the vortex characteristics. For the 40° pitch case (Figure 10(d)), the LEV is comparatively small, tightly formed, and relatively far from the wing’s surface. The LEV is completely formed in the 50° case (Figure 10(e)) and is quite diffused. The LEV is also relatively close to the wing’s surface. For the 60° case (Figure 10(f)), the LEV has clearly burst and a region of counter-rotating vorticity has formed near the wing surface. Note that the LEV is significantly larger than those occurring in the 40° and 50° cases, encompassing a large area of the wing surface and inducing much stronger vortical flow about the wing.

At t/T = 0.35, the flows for the 50° and 60° cases, Figure 10(h) and (i), respectively, show signs of separation. The flow has become more incoherent and exhibits small formations of shed vorticity. On the contrary, the LEV for the 40° case (Figure 10(g)) has grown in size and become more diffused, but remains attached to the wing’s surface. The flow features shown in Figure 10 help to clarify the difference in the force trends in Figure 9. The delayed formation of the LEV in the 40° case, as compared to the 50° and 60° cases, is responsible for the noticeable phase difference in the lift force trend in comparison to the other two pitch angle cases. Additionally, this explains why the drag force trend for the 40° case has a shallow peak while the 50° and 60° pitch cases exhibit steeper, more distinct peaks.

Computational surface pressure contours

In the previous section, we compared the evolution in time of the LEV for all three translational pitch cases. The temporal, as well as spatial, variation of the LEV strongly affects the net surface pressure distribution along the wing. This, in turn, determines the aerodynamic forces generated by the wing. Analyzing and comparing the pressure distribution across the wing surface can give a more global view on the influence of the LEV on the wing.

Figure 11 shows the spanwise distribution of pressure coefficient (C_p) contours along the top and bottom surface of the wing at t/T = 0.25 for all three pitch cases. Note that in the absence of a freestream velocity, pressure coefficient is define as

$C_p=p-p_{\infty }{p}_{\infty }$

(4)

where p is the static pressure of some point in the flow and p_∞ is the freestream static pressure. A C_p of zero indicates that the pressure is equal to the freestream pressure. Positive C_p values indicate that the local pressure is greater than p_∞ while negative C_p values indicate that the local pressure is less than p_∞. At t/T = 0.25, the wing is at its highest rotational velocity and the local velocities along the wing span are at a maximum. For all three cases, the pressure contours along the bottom surface of the wing are quite similar. The 40° case (Figure 11(a)) has slightly lower pressure on the bottom surface, in comparison to the 50° and 60° cases (Figure 11(b) and (c), respectively), due to its smaller nominal translational pitch angle. Most of the noticeable difference in pressure distribution along the wing occurs on the top surface which is subject to the LEV.

Figure 11.

Pressure coefficient contours on the wing surface for the 40°, 50°, and 60° pitch cases. (a) 40° case, t/T = 0.25, (b) 50° case, t/T = 0.25, and (c) 60° case, t/T = 0.25.

In the 40° case (Figure 11(a)), the global influence of the LEV is comparatively weak, with the lowest pressure regions occurring along the leading edge near the wing root. At midstroke, most of the wing’s top surface is at, or slightly lower, than the freestream static pressure. For the 50° case (Figure 11(b)), the entire top surface of the wing is at a lower pressure than the freestream static pressure. The LEV has a much stronger global effect on the wing with the significant low pressure region (denoted in green) growing until about the 70% span location. After approximately the 70% span location, the LEV has burst and the wing is subject to separated flow resulting in slightly higher pressure near the wing tip. The 60° case, shown in Figure 11(c), has a similar pressure distribution to the 50° case with the region of low pressure caused by the LEV also growing until about the 70% span location. However, one visible difference is the size of the lowest pressure region (denoted in blue) at the leading edge near the wing root. This small, concentrated region of low pressure, due to an attached LEV, is much larger for the 60° case, in comparison to the 50° case, extending to approximately the 10% span location.

From Figure 11, we see that the LEV at midstroke has a much stronger influence on the net surface pressure along the wing for the 50° and 60° translational pitch cases in comparison to the 40° pitch case. How the LEV and the low pressure region it creates grow throughout the flap cycle not only impacts the aerodynamic force production but also the aerodynamic power as well.

Aerodynamic power and efficiency

As the design and development of future MAV concepts progresses, it is important to assess not only the aerodynamic performance but also the aerodynamic efficiency of a particular design. In conventional aircraft design, typical measures of aerodynamic efficiency include power (P), the lift-to-drag ratio (L/D), and power loading (L/P). Note that throughout this section, references to power infer aerodynamic power. The inertial power of the flapping mechanism was not investigated in this study.

Aerodynamic power is calculated by summing the dot product of the local aerodynamic force vector with the local velocity vector along the span of the wing. Figure 12 shows the variation of instantaneous power over a flap cycle for the 40°, 50°, and 60° translational pitch cases. As expected, the peak magnitude of aerodynamic power increases with increased translational pitch angle and approaches a magnitude of zero during stroke reversal. Slight differences in the peak magnitude and variation of the instantaneous power trend line occur between the forward and backward stroke for all of the translational pitch cases. This is because of fluctuations in the measured wing pitch angles, which were discussed in previous sections. If perfectly symmetric pitch kinematics were used, as opposed to the experimentally measured values, it is expected that the results of the forward and backward strokes would be symmetric.

Figure 12.

Variation of instantaneous power over a flap cycle (40°, 50°, and 60° cases).

Two metrics will be utilized in determining the aerodynamic efficiency of the wing kinematics used in this study: lift-to-drag ratio (L/D) and lift per unit power (L/P). For this study, the instantaneous variation of L/D and L/P, only the translational phases of the forward stroke and backward stroke (t/T = 0.1—0.4 and 0.6—0.9) will be analyzed. The instances of the flap cycle where stroke reversal is occurring are ignored because the wing produces negative lift at those times and was not operating at its highest possible efficiency. Figure 13 shows the variation of L/D over the translational portion of the forward and backward strokes. In calculating L/D, the absolute value of instantaneous drag is used because the sign of the drag value only denotes the direction in which drag acts for the given coordinate system. The variation in L/D is symmetric between the forward stroke and backward stroke for all the cases. For the 50° and 60° cases, L/D is fairly consistent for a majority of the translational phase. However, L/D in the 40° case shows a large amount of variation as the wing is still undergoing stroke reversal pitch rotation.

Figure 13.

Lift-to-drag ratio during the translational phase of forward and backward stroke for the 40°, 50°, and 60° cases. (a) Translational phase: forward stroke and (b) translational phase: backward stroke.

The power loading variation over the translational phase of the flap cycle for the forward stroke and backward stroke is shown in Figure 14. Figure 14(a) also contains experimental data for instantaneous lift-to-power ratio presented in Benedict et al. ²⁶ Note that the available experimental data are only for the forward stroke from approximately t/T = 0.2−0.4, which omits regions of negative lift in the lift-to-power ratio trends. Similar to Figure 13, there is little difference between the computational trends predicted in the forward stroke in comparison to the backward stroke. The CFD 50° and 60° cases have a relatively constant L/P during the translational phase similar to the experimental data for those two cases. However, the computational trend lines initially under-predict the magnitude of L/P (especially for the 50° case) and then over-predict L/P toward the end of the translational phase. For the 40° case, the CFD L/P curve is monotonically increasing while it is relatively constant for the respective experimental L/P curve. Initially, the instantaneous L/P for the computational 40° case is negative because of the negative lift produced at the beginning of the specified segment of the flap cycle. However, toward the end of the translational phase, the magnitude of the lift-to-power ratio for the 40° case far exceeds that of the other two cases. In general, the CFD model predicts a similar trend in lift-to-power ratio in comparison to the 50° and 60° cases. However, it tends to under-predict the magnitude of L/P compared to the experiment. There is a noticeable discrepancy between the experimental and CFD-predicted L/P trends for the 40° case and needs further investigation. It is believed that due to the larger variation in pitch angle for the 40° case, the CFD model may be prolonging the diffusion of the LEV and consequently over-predicting the lift-to-power ratio toward the end of the stroke.

Figure 14.

Lift per unit power during the translational phase of the forward and backward stroke for the 40°, 50°, and 60° cases. (a) Translational phase: forward stroke (with experimental data) and (b) translational phase: backward stroke.

Table 1 contains the time-averaged values of CFD-predicted lift, drag, power, lift-to-drag ratio, and power loading for the three translational pitch cases tested. The values in Table 1 are the averaged instantaneous values during the translational phase of the forward stroke. The 50° case had the best lift-to-power ratio at 0.19 Newton/Watt followed by the 40° case, and the 60° case being the worst. While the 50° case also had the best lift-to-drag ratio at 0.75, the 60° had a slightly better L/D in comparison to the 40° case. The low L/D for the 40° case is due to the negative lift generated at the beginning of the translational phase.

Table 1.

Time-averaged CFD force, power, and efficiency values for the three translational pitch cases.

Pitch angle	40°	50°	60°
Lavg (g)	7.03	13.02	13.65
Davg (g)	12.44	17.38	23.50
Pavg (W)	0.44	0.68	0.93
Lavg/Davg	0.56	0.75	0.58
Lavg/Pavg (N/W)	0.16	0.19	0.14

CFD: computational fluid dynamics.

When taking advantage of the lift enhancement capabilities of the LEV, the lift-to-drag ratio using the given wing kinematics was less than 1.0. While rigid flapping wings have demonstrated some initial promise, further improvements in their aerodynamic efficiency is necessary to make the flapping-wing platform more practical for future MAV applications. One way to improve aerodynamic performance would require the variation of the wing pitch angle along the span throughout the flap cycle. Proper spatial and temporal variation of the pitch angle (dynamic twist) would enable the wing to sustain a stable LEV over a greater region of the wing surface. Dynamic twisting can be achieved using a flexible wing with some torsional compliance and an appropriate degree of spanwise and chordwise flexibility.

Conclusions

The goal of the present study was to validate the predictive capability of a CFD solver to calculate the unsteady aerodynamic mechanisms produced on a rigid insect-based flapping wing at MAV-scale Reynolds numbers. Data from previous experiments conducted in air at a constant flap frequency of 5 Hz at translational pitch angles of 40°, 50°, and 60° were used for validation purposes. The experimental aerodynamic forces were measured using subtraction of the inertial force from the total force, and PIV studies were conducted to investigate the flowfield in the immediate vicinity of the wing’s leading edge. Variations in the wing flapping and pitching angles measured during experimentation were used in the CFD model to replicate the wing kinematics. Analysis of the experimental results allowed for the validation of the CFD by comparing the unsteady instantaneous aerodynamic forces, evolution of the vorticity field over time, and LEV strength. The validated CFD model was further used to understand the unsteady aerodynamic force production on a flapping wing and the influence of the passive pitch kinematics on the generation of the LEV. Additionally, the CFD results were used to investigate the aerodynamic efficiency of the flapping wing. The following specific conclusions are drawn from the study:

The CFD model was able to resolve the salient flow features present for the flapping wing, particularly toward the inboard section of the wing. While the degree of flow separation and vortex diffusion were not predicted as well toward the wing tip, the CFD model was able to both spatially and temporally predict many of the general flowfield characteristics. Satisfactory correlation of the CFD and PIV flowfields manifested in good agreement between the measured and calculated unsteady aerodynamic force trends. This provides confidence in the use of an inertial force subtraction technique to determine the instantaneous aerodynamic forces generated by a flapping wing in air. Overall, the CFD results showed good agreement with the experimental data.

The flowfield was highly susceptible to changes in the pitch angle resulting in corresponding changes in the predicted aerodynamic force trends. Specifically, the delayed pitch rotation kinematics, inherent to the passive pitching mechanism utilized throughout the experiment, led to the production of negative lift during stroke reversal. In addition, the delayed rotation kinematics resulted in a noticeable phase shift in the force-time history for the 40° translational pitch case in comparison to the 50° and 60° cases, lowering the magnitude of lift produced during the translational phases of the flap cycle.

The 40° translational pitch kinematics led to more stable formation of the LEV during the translational phase of the flap stroke. At 50° and 60° translational pitch angles, LEV formation occurred earlier in the flap cycle resulting in increased LEV size and concentration. However, LEV convection, diffusion, and bursting also occurred earlier in the translational phase of the flap stroke especially for the 60° case. For all translational pitch cases, the LEV temporal behavior, including creation, residency time, and destruction, resembled that of dynamic stall.

The aerodynamic performance and efficiency of the various wing kinematics tested was quantified using lift-to-drag ratio and power loading. CFD-predicted instantaneous power loading curves show similar trends to the experimentally acquired data for the 50° and 60° translational pitch cases. However, there was a noticeable discrepancy between the experimental and computational power loading curves for the 40° pitch case. Maximum L/D during the translational phase of the flap stroke was 0.75 for the CFD 50° translational pitch case. The CFD 50° pitch case also had the highest power loading of 0.19 N/W followed by 0.16 N/W and 0.14 N/W for the 40° and 60° translational pitch cases, respectively.