Abstract

We report a comprehensive scaling law and novel lift generation mechanisms relevant to the aerodynamic functions of structural flexibility in insect flight. Using a Navier–Stokes equation solver, fully coupled to a structural dynamics solver, we consider the hovering motion of a wing of insect size, in which the dynamics of fluid–structure interaction leads to passive wing rotation. Lift generated on the flexible wing scales with the relative shape deformation parameter, whereas the optimal lift is obtained when the wing deformation synchronizes with the imposed translation, consistent with previously reported observations for fruit flies and honeybees. Systematic comparisons with rigid wings illustrate that the nonlinear response in wing motion results in a greater peak angle compared with a simple harmonic motion, yielding higher lift. Moreover, the compliant wing streamlines its shape via camber deformation to mitigate the nonlinear lift-degrading wing–wake interaction to further enhance lift. These bioinspired aeroelastic mechanisms can be used in the development of flapping wing micro-robots.

1. Introduction

Insects are among the most commonplace animals in daily life, but aerodynamically their flight characteristics and performance are stunning yet inadequately understood. To date, existing theories cover unsteady aerodynamics beyond the steady, stationary wing theories. In particular, clap and fling [1], delayed stall via prolonged leading-edge vortices (LEVs) [2], wake-capture of a wing during the return stroke [3,4], and rotational forces due to combined pitching and plunging [3,5] offer intriguing insight into the lift generation required for an insect to stay aloft. These unsteady physical processes are particularly important for small flyers, such as fruit flies, which operate at a Reynolds number Re = Uc/ν around O(102), based on the midstroke velocity U, wing chord length c and kinematic viscosity of air ν.

However, the aforementioned theories are based on rigid wing studies, neglecting the effects of wing flexibility. Flexible structures interacting with surrounding fluid are not rare in nature [6]. For organisms, such as tree leaves [7] or benthic organisms [8], the compliance nature can lead to drag reduction, which may be favourable for their survival [6,9,10]. The Reynolds number for species such as hummingbirds is Re = O(104), while for smaller insects, such as fruit flies or honeybees, Re = O(102 − 103). These small insect flyers have flexible wings that experience large deformations just like birds and bats [11,12]. It has been, for example, shown that wing deformations [13,14] can enhance the force generation and efficiency of a locust operating at Re ≈ 4 × 103 [13], or a hawkmoth at Re ≈ 6.3 × 103 [15]. Moreover, several studies have suggested that insect wing rotations may be passive [16,17], meaning that the resulting rotation is due to a dynamic balance between the wing inertial force, elastic restoring force and fluid dynamic force.

It is established that shape adaptation associated with flexibility can affect the effective angle of attack and hence the aerodynamic outcome [18]. However, our understanding of fluid physics and the resulting structural dynamics has been insufficient to explain all the salient features of this coupled fluid–structure system. Furthermore, the relationship between structural flexibility and flapping kinematics has not been fully explored. The rigid wing studies [3,19] indicate that the timing between translation and rotation of the wing can be critical to lift generation. Specifically, when the wing is in advanced rotational mode, i.e. the wing rotates before it reaches the end-of-the-stroke, the highest lift coefficients are obtained [3,19]. Moreover, tethered flying Drosophila actively control left/right wing rotation timing to initiate yaw turns [20]. However, hovering fruit flies [21] and honeybees [22], beetles [11] and hymenopterans in general exhibit symmetric rotation [11], which means that they rotate their wings around the same time as they reverse the stroke.

The effects of wing flexibility can be explored from the static and dynamic perspectives. For the static approach, we can study the equilibrium shape compliance in accordance with the stress distributions around the body. The resulting modification of the angle of attack and the wing shape can potentially better adapt to the surrounding flow, reducing tendencies towards massive flow separation and stall. Alternatively, structural dynamics during transient processes can offer further time-dependent behaviour other than static shape compliance, which is the dynamic approach. The interplay of the fluid flow with the shape compliance and transient response is nonlinear and intriguing even at Re = 100, which is addressed in this study. In this Reynolds number regime, the fluid flow can be considered as laminar and the computational accuracy of the Navier–Stokes equation solver employed in this study is satisfactory [19]. The instantaneous flow field is synchronized with the dynamic displacement field of the wing structure in a carefully validated fully coupled aeroelastic framework [23].

We show that by adjusting the frequency ratio f/f1, which is the ratio between the motion frequency f and the first natural frequency of the wing f1, passive rotation can also yield advanced, symmetric or delayed rotation modes, which are the key elements for the aerodynamic performance of rigidly rotating wings [3,5]. The natural frequency f1 is measured in the chordwise direction between the leading-edge (LE) and trailing-edge (TE). Moreover, we assess the effects of wing flexibility by dissecting the aeroelastic response into the time evolution of the wing motion and the camber deformation of the wing shape.

For our purpose, it is sufficient to address only the chordwise flexibility and represent the flexible wing as a homogeneous elastic flat plate [24]. While three-dimensional effects, such as spanwise flow that seem to stabilize the LEVs [25] or LEV–tip-vortex interaction [26] are notable in general, the chordwise flexibility is essential and warrants an independent study [18]. Furthermore, at Re = O(102), the effects of spanwise flow on the overall aerodynamics are less important than at higher Reynolds numbers [25,27]. Also, the characteristics of the LEVs in two dimensions for plunging motions are representative of three-dimensional flapping wings as long as the stroke-to-chord ratio is within the range of typical insects, i.e. around four to five [28–30], which we consider in this study.

2. Material and methods

2.1. Non-dimensional parameters and case set-up

For a complex physical system, dimensional analysis via the Buckingham Π-theorem and scaling analysis reduces the number of non-dimensional parameters by recognizing specific properties of the system in consideration. This can be helpful even if a complete mathematical solution is missing [31]. We approximate the insect wing with an elastic flat plate [24] with a chord of c made of a material with density ρs, Young's modulus E, and thickness hs. Furthermore, we consider a flapping motion with an amplitude of ha and frequency of f in a fluid of density ρf and viscosity ν, resulting in lift per unit length L. The dimensional analysis leads to a set of non-dimensional parameters, which we choose as the density ratio ρ* = ρs/ρf, thickness ratio hs* = hs/c, reduced frequency k = 2πfc/(2U), frequency ratio f/f1 = 2πf/[k12{Ehs*2/(12ρfc2)}1/2], Reynolds number Re, and the coefficient of lift CL = L/(0.5ρfU2c).

Based on high-fidelity numerical computations of a wide range of wing motion, wing structure and geometry, and surrounding fluid, we found that a relative shape deformation parameter, γ, plays a pivotal role in characterizing the force and propulsion characteristics of a flapping, flexible wing [23]. Specifically, γ is defined as2.1which includes the effects of f/f1, k, hs* and ρ*, such that the individual effects of these materials and flapping parameters are integrated into a single parameter.

Here, we consider a non-dimensional flow field with unit density initiated by a hovering two-dimensional flat plate with unit chord, thickness ratio of 0.02 and flat edges [19]. The flow field is governed by the unsteady Navier–Stokes equations with constant density and viscosity as2.2for the velocity u, pressure p and time t. The variables are non-dimensionalized with the reference velocity U, inverse of the motion frequency 1/f and chord c as the velocity, time and length scales, respectively. The superscript asterisk (*) indicates non-dimensional variables. In absence of a freestream in hovering wings, we take the maximum translation velocity U of the flat plate at the LE as the reference velocity, such that U = 2πfha = 1 [19,22]. Note that reduced frequency in hover becomes a geometric quantity: k = c/(2ha) [19]. The Reynolds number Re = 100 corresponds to conditions at level of a fruit fly.

A sinusoidal plunge motion is imposed on the LE of the wing as2.3To isolate the sources of aerodynamic force modification, we consider a flexible plunging wing operated in translational motion only by prescribing a sinusoidal displacement on the LE of the wing (figure 1). The flexible wing structure is modelled locally by2.4where w is the displacement due to bending motion and ff is the distributed transverse fluid force on the wing per unit span. Note that in the case of a spanwise flexible wing, the effects of the aspect ratio of the wing need to be taken into account [23]. Consequently, the resulting wing camber deformations and rotation α(t*), the angle between the TE and LE, are purely passive due to dynamic balance between the aerodynamic loading, elastic restoring force and inertia of the wing. For dynamic similarity of equation (2.4), both the stiffness and the aerodynamic force term, normalized by the inertia term, need to be matched. The former results in f/f1 and the latter can be characterized by a combination of k, ρ* and . The parameter γ includes all three, i.e. inertia, stiffness and aerodynamic force terms.

Schematic of the abstraction of the wing motion. (a) Top view of a Drosophila melanogaster. (b) Front view of a Drosophila melanogaster wing and schematics of two-dimensional cross sections of the wing, indicated by red solid lines. The LE of the wing is indicated by a red dot and wing shapes are shown at seven equally placed time instants during a flapping stroke to the left, illustrated with dashed arrows. (c) Two-dimensional model of the flapping stroke shown in (b). The wing is represented as a flat plate with constant thickness. A sinusoidal plunge motion with amplitude of ha is imposed at the LE. Wing deformations lead to a passive pitch angle of α. The angles at the mid and end of the strokes are αm and αe, respectively. The directions of the lift and drag are normal and parallel to the imposed plunge motion, respectively The pictures of the Drosophila melanogastor are used with permission of Nicolas Gompel and Benjamin Prud'homme.

The density ratio is ρ* = 7.8, similar to steel in water or a light material in air. The current case set-up is motivated by the rigid wing experiments in mineral oil [3–5,32]. For insects, wing bending is mainly due to wing inertia as the density ratio ρ* is typically of the order O(103) [23,24]. The mineral oil has a density that is 7×102 times greater than air, such that the acceleration–reaction force, which is often linearized by the added mass, becomes the dominant force that deforms the wing [23]. Still, the study of the aeroelastic response of insect wings in such a low density ratio system is justified as long as the resulting aeroelasticity matches the actual motion [32] as illustrated by a rigid wing with ρ* = 1.4 [3], crane fly [17] or shown by the γ-scaling analysis [23].

The frequency ratio f/f1 and the reduced frequency k are varied by changing E and ha, respectively, to probe their influence on the resulting aerodynamics and the structural deformations. It is reported that the natural flyers operate at a frequency ratio less than 0.8 [33] and in this study, we consider 0.06 < f/f1 < 0.82. The range selection of ha is motivated by biological flyers [28], such that 0.25 < k < 3.75 (see the electronic supplementary material, table S1). In the main presentation, we replace f/f1 by γ given by equation (2.1) for better clarity.

The time-averaging operator is defined as2.5for example, for CL. To avoid initial transient effects, we take the time-average between the third and the fifth cycle, i.e. m = 3. The time histories of forces are plotted by phase averaging the six forward and backward strokes within cycles 3–5. Moreover, to account for the extension of the wing chord length, measured from the LE to TE, owing to the use of the linear elastic model for the wing at large deformations, we use the instantaneous chord length c(t*) in the non-dimensionalization of the lift coefficient. As a result, the resulting lift coefficient is always lower than when c is used in the non-dimensionalization.

In order to elicit the effects of the nonlinear transient dynamics, we also consider an approximation of α with a first-order harmonic (FH) as2.6

The two unknowns, the angular amplitude αa and the phase lag φ between the plunge and angular motions, are calculated by measuring αm = α(0) and αe = α(1/4) (figure 1), the angles at the mid and end-of-the-strokes, respectively. Matching αm = αFH(0) and αe = αFH(1/4) yields a system of two equations for the two unknowns αa and φ.

We systematically assess the resulting interplay between the fluid physics and structural response by constructing two rigid flapping wing models based on the resulting angular motion of the flexible wing. Specifically, we impose α and αFH on a rigid flat plate without camber. A rigid wing rotating with α has an angle-of-attack identical to that of a flexible wing, but without the wing shape variation in camber. The wing with αFH represents the transient dynamics without either shape deformation or the nonlinear response of the wing motion. We decompose the effects of flexibility into (i) the nonlinear transient dynamics by comparing the lift generation on the two rigid wings without camber deformation and (ii) the shape compliance by comparing the flexible wing and the rigid wing with the same time history of pitching angle α.

2.2. Numerical models

Full details of the fluid–structure interaction (FSI) framework and careful validation analysis against well-documented experimental results can be found in our previous work [23]. The governing equations for the fluid given by equation (2.2) are solved using an in-house three-dimensional, unstructured, pressure-based finite volume solver [34,35] written in a rule-based framework [36]. The geometric conservation law [37] is satisfied [38]. We use the radial basis function interpolation method to deform the computational mesh at each time step [39]. The spatial and temporal sensitivity studies as well as computational set-up are shown in appendix A.

The dynamics of the wing structure given by equation (2.4) is solved using a finite-element representation [23]. The structural damping is not considered in this study. Two degrees of freedom, namely the displacement and bending, are allowed at each node. Computations performed for an aerofoil composed of a rigid teardrop and elastic flat plate at higher Reynolds numbers and for various motion frequencies [23] showed that a linear Euler–Bernoulli beam is sufficient for qualitative analysis of the FSI coupling. The flat plate is modelled with 51 nodes equally distributed over it.

The FSI coupling is a time-domain partitioned solution process. The governing partial differential equations for the fluid and the structure are solved independently and coupled spatially through the interface between the fluid and the structure. At each time step, we iterate the flow field and the structural displacement field such that sufficient convergence on the deformed wing displacement is achieved, before advancing to the next time step. In order to accelerate and ensure the convergence of the FSI, the Aitken relaxation method has been incorporated [23].

3. Results and discussion

3.1. Lift scales with the relative shape deformation parameter γ

There exists a scaling relationship [23] between the time-averaged lift coefficient and γ as,3.1which is derived by considering the non-dimensional energy balance on the deforming wing for . When a sinusoidal plunge motion is imposed on the LE of the wing, the amplitude of the TE displacement relative to the LE is proportional to γ. Based on this scaling relationship, the power input and the propulsive efficiency could be scaled using γ [23].

One of the conditions made in the previous derivation is that . As a consequence, the scaling for the TE velocity became 2πf/f1γ [23]. When we apply equation (3.1) to the current cases, only the cases with f/f1 < 0.6 follow equation (3.1), while 0.6 < f/f1 < 0.8 deviate. To account for 0.6 < f/f1 < 0.8, we revise the scaling for the TE velocity to 2πγ, which can be interpreted as relative TE velocity scale being proportional to the ratio between the relative tip deflection γ and the non-dimensional motion period. Then, without making further approximations, the non-dimensional energy balance leads to3.2where a decreases from 2.0 for γ < 1 to 0.34 for γ > 2. These scaling exponents were determined using a linear fit to the current data in a least-squares sense. The coefficients of determinant were 0.94 and 0.71, respectively. When the wing deformations become more significant at γ > 0, the correlation worsens. This indicates that other factors may start to play a role as equation (3.2) only considers FH as insects usually flap with f/f1 < 1 [33] and assumes aerodynamics forces from a high-frequency small-deformation motion [23]. Compared with equation (3.1), equation (3.2) has a corrector for the frequency ratio effects. Note also that when f/f1 → 0, the normalization factor of lift coefficient in equation (3.2) approaches to that in equation (3.1). When f/f1 is small, the elastic, restoring effects dominate over the inertia effects implying the work done on the wing is mostly in balance with the restoring elastic force. As f/f1 increases, the influence from the wing inertia cannot be neglected anymore and a correction needs to be made.

The outcome is that the time-averaged lift coefficients versus γ covers a wide variety of wing flexibility and flapping characteristics relevant to insects as well as artificially devised flexible wings (figure 2a). Moreover, γ can be interpreted as a non-dimensional TE displacement relative to the LE. The chordwise wing deformations act as passive pitch with an angle of attack α, which is a key measure of the aerodynamics and directly affects lift (figure 2b,c). Furthermore, figure 2d,e shows that the normalized power required [23] and the propulsive efficiency defined as the ratio between the lift and the power required for the current cases scale with γ.

γ-scaling relationships. (a) Scaling given by equation (3.2) of the normalized lift by γ. The current dataset (plus symbols) for the chordwise flexible aerofoils in hover at Re = 1.0 × 102 is plotted together with wide range of cases reported by Kang et al. [23]: plunging two-dimensional chordwise aerofoils at Re = 9 × 103 in forward flight in water (circles), plunging three-dimensional spanwise wings at Re = 3 × 104 in forward flight in water (diamonds), flapping three-dimensional isotropic wing in hover in air at Re = 1.5 × 103 (squares) and insects (cross symbols), where the lift for the insects is estimated by their weight [23]. The letters c and s indicate chordwise and spanwise flexibilities, respectively. (b,c) Correlations between , maximum angle normalized in the same manner as for , and γ. (d,e) Normalized aerodynamic performance against γ from the current computations for (d) time-averaged power required and (e) propulsive efficiency , which is the ratio between the time-averaged lift and the time-averaged power required. In (b,c), the colours indicate advanced (φ > 95°, blue), symmetric (85 < φ < 95°, red) and delayed (φ < 85°, black) rotation modes.

From equation (3.2) it follows that optimal lift is obtained when f/f1 = 0.5 with other parameters fixed. For a self-propelled insect model [33], the best performance was observed at f/f1 = 0.7 and the performance decreased by a factor of more than 4 at f/f1 = 1 consistent with the current prediction. The current discussion also suggests that the maximum lift may be obtained at a frequency ratio well below the resonance condition because increasing the TE deformation relative to the LE results in delayed rotation mode, deteriorating the lift generation. Finally, we see that the parameters used in this study, i.e. ρ* = 7.8, and k around 0.6 result in a factor introduced in the definition γ given by equation (2.1), which is in a comparable range for the typical values for fruit flies with ρ* = 1100 [23], [23], and k = 0.212 [18]. This means that the effects of f/f1 on γ given by equation (2.1) in this study will be similar to that of a fruit fly.

3.2. Symmetric rotation is optimal

Additional forces arise at the stroke reversals for rigid wings. When actively rotating wings exert force on the surrounding fluid, the fluid adds to or subtracts from the lift due to translational mechanisms [3,40]. The benefits of the rotational mechanism depend on the phase difference between translation and rotation [3]. Previous studies have shown that the advanced rotation is optimal for rigid wings [3,19].

Depending on the wing structure and imposed motion, the passive wing rotations also result in the three phase modes (figure 3a). Moreover, φ strongly correlates to f/f1, yielding the advanced, symmetric and delayed rotation modes with increasing f/f1 (figure 3b). However, the large pressure differentials that exist on actively rotating rigid wings near the TE due to lift-enhancing rotational effects [40] are relaxed by the compliant nature of flexible wings. Instead of generating rotational forces, the wing streamlines its shape such that the wing shape and motion is in equilibrium with the fluid forces, similar to the drag-reducing reconfiguration of flexible bodies [6]. Passive rotational angle α is purely due to deformation where the amplitude and the phase scale with γ and f/f1. Moreover, αm is measured at the instant of maximum translational velocity and is therefore indicative of the forces related to the translational mechanisms [3]. For advanced rotations (f/f1 < 0.25), αm is greater than 70°. The wing orientation is almost vertical, producing little lift. On the other hand, deformations can become substantial for delayed rotations modes with f/f1 > 0.4, but since the translation is out-of-phase with the passive rotation, the resulting αm and hence the lift remain smaller than for symmetric rotations. Figure 3c also shows that the resulting lift for the delayed rotations can be higher than for the advanced rotation modes, unlike the actively rotating rigid wings [3]. The translational forces peak when αm is around 45°–50° [3], which corresponds to 0.25 < f/f1 < 0.4 (figure 3b). The resulting lift is of similar magnitude to that of rigid wings with active rotations [3]. However, the lift is now optimal for symmetric rotation modes (figure 3c). This is also observed in the wing kinematics of fruit flies [21] and honeybees [22].

3.3. Flexible wings outperform their rigid counterparts

To further assess the effects of flexibility, we compare the lift generation on the flexible wing to the lift on the rigid wing rotating with αFH. Both wings generate lift that correlates with γ (figure 4a). For the cases considered, most flexible wings yield higher lift with a difference that depends on γ/k2 (figure 4b). As γ increases, wing deformations become larger and will eventually cause the observed differences. The dependence on k2 can be explained by considering the relative contribution of the added mass force and the forces induced by the vorticity in the flow field [22,23], which scale as k and 1/k, respectively, for hover [23]. The ratio between these two forces is then k/(1/k) = k2 [22,23] and the added mass gains its relative contribution to the total lift with increasing k.

Difference in lift between flexible and rigid wings, rotating with αFH. (a) on the flexible and rigid wings as a function of γ, (b) r.m.s. difference in the lift as a function of log10γ/k2. Time histories of lift for the cases and schematics of wing shapes with the smallest and largest differences in lift are shown in (c,d) and (e,f), respectively. Both cases exhibit symmetric rotations with (c) a high k = 3.05 (γ = 1.23) and with (e) a low k = 0.6 (γ = 2.59) for the flexible (red), rigid (blue) wings and their instantaneous wing shapes (d,f). The integrated aerodynamic forces are indicated with the blue arrow.

As such, at higher k the added mass force, which is linearly proportional to the wing acceleration, has a larger contribution to lift [22]. The impact of the added mass on the total aerodynamic force for both wings remains similar, because the essence of wing motion is consistent between flexible and rigid wings. Consequently, the resulting lift is close () as shown for k = 3.05, γ = 1.23 in figure 4c,d. Hence, a high k corresponds to small stroke amplitude, high-frequency motion as illustrated in figure 4d.

When k is lowered, the aerodynamic force induced by vorticity in the flow field, such as translational forces (equation (B 2)) and nonlinear wing–wake interactions [3,4], starts to dominate the lift [22,23]. These nonlinearities in the unsteady aerodynamics lead to intriguing consequences caused by small differences in the resulting nonlinear wing motion and the camber deformation, which will be discussed further in §3.5. For example, figure 4e illustrates the case (k = 0.6, γ = 2.59) corresponding to the largest . Lift on the flexible wing is considerably higher during the midstroke compared with that of the rigid wing (figure 4e,f), such that the lift generated by the flexible wing is superior: 1.2 versus 0.88. The effects of wing flexibility are further assessed in the next sections by decomposing the resulting aeroelastic response into the nonlinear temporal evolution and camber deformation of the wing.

3.4. Lift enhancement due to nonlinear wing motion

The simple representation αFH in the study of flapping wings is a popular approximation for rigid wings [19,28], despite the observations of hovering fruit flies and honeybees [21,22] using wing strokes that are clearly far from sinusoidal. Similar to the stroke used by these insects, the resulting flexible wing motion significantly deviates from being a sinusoid, depending on the wing properties (figure 5a). To assess the effects of the resulting nonlinear structural motion without the influences from camber deformations, we compare the lift generation on two rigid wings that follow the full motion α and the first-order approximation αFH, respectively. The former means that the modified angle-of-attack is identical to that of the flexible wing, but without the wing shape variation in camber. The latter represents the structural dynamics without either shape deformation or the nonlinear response of the wing motion.

In general, the lift is higher for rigid wings with α. The dynamic balance between the fluid and the wing structure results in a rotation with larger deformation indicated by αmax = max(90 − α) near the midstroke (figure 5a). Correspondingly, the instantaneous lift on the rigid wing with α results in higher magnitudes during the midstroke (figure 5b). The overall effect is that the rigid wing with α generates higher lift that correlates to the difference between αmax and αmax,FH, as illustrated in figure 5c for most cases. Inversely, this means that better performance can be obtained by departing from a simple sinusoidal rotational motion.

3.5. Streamlining of wing deformation enhances lift

The wing–wake interaction provides an important mechanism responsible for lift enhancement over a flexible wing compared with its rigid counterpart, as highlighted in figure 6. Here, we consider a flexible wing and its rigid counterpart, rotating with α, to focus solely on the effects of camber deformations. The LEV and the trailing-edge vortex (TEV) shed in the previous motion stroke, denoted as LEV0 and TEV0, respectively, form a vortex pair and induce a downward wake around the midstroke [4,19], which in turn interacts with the wing during its return stroke. The outcome of the rigid wing–wake interactions varies. Under favourable conditions, added momentum causes the lift to increase during the first portion of the stroke [3]. The wing then experiences two peaks in lift: the wake-capture and the delayed stall [3,19] peaks as shown in figure 4e. However, the wake can also form a downward flow, causing the lift to drop significantly during the midstroke [4,19].

Effects of streamlining for the symmetric rotation mode (k = 0.6; γ = 2.59) with the largest difference in lift compared with the flow around the rigid wing rotating with α. (a) Representative vorticity and vertical velocity fields at selected time instants during a backward stroke. The wing thickness is exaggerated for clarity. (b) Schematic of the wing–wake interaction. The subscript numbers indicate the stroke at which the vortices are shed: current (1), previous (0) and two strokes ago (−1). The arrows between the LEV0 and TEV0 illustrate the downwash. The solid lines illustrate higher downwash magnitude than the dashed lines. (c) Time history of lift normalized by c and c(t*). (d) as a function of γ. (e) Time histories of the vorticity of TEV0, downwash and LEV1.

A flexible wing, on the other hand, can reconfigure its shape, adjusting its camber to better streamline itself to the surrounding flow. Consequently, the formation of the TEVs reduces for the flexible wing and the negative impact of the induced downwash is mitigated (figure 6a,b), leading to higher lift during the midstroke. For a locust, positive camber was measured in the hind wing due to the compression of the radial veins of the vannus via TE tension [14]. The compression causes the radial veins to undergo Euler buckling [41]. The resulting positive camber is called the umbrella effect [14,41], enhancing lift in the downstroke [13,14]. Compared with the locust, whose Reynolds number is Re ≈ 4 × 103 [13], figure 6 focuses on Re = 100 and different detailed wing characteristics. These lead to different observations: the camber deformation is purely passive as a result of the dynamic balance between the fluid, elastic restoring and inertial forces. While the resulting negative camber is associated with less favourable aerodynamic performance as shown from the inviscid steady state thin aerofoil theories [42], this may not be the case for unsteady flapping flexible wings at low Reynolds numbers. A flexible wing is able to generate translational lift with higher efficiency via the streamlining mechanism, while at stroke reversals, it lessens rotation-related force enhancement. As a result, the time history lift changes from the two peak shape to one peak in the middle of the stroke (figure 6c). The lift enhancement increases with increasing deformations, characterized by γ (figure 6d). Note that the lift coefficient is even greater when is normalized by the original chord c instead of the instantaneous chord c(t). To quantify this streamlining process, we estimated the magnitudes of the non-dimensional downwash vd and vorticity ∇* × u* of the vortices TEV0 and LEV1 based on the computational results (see the electronic supplementary material, figure S1). The vorticity magnitude of TEV0 for the rigid wing is higher, which indeed correlates to a stronger downwash (figure 6e). Figure 6e also shows that the stronger downwash delays the increase of the LEV for the rigid wing.

The angle of attack is a key measure of the aerodynamics and can directly affect features such as LEV and lift. Since the flow field surrounding a hovering wing is complicated, one needs to define an effective angle of attack, which we model as a combination of the translational velocity of the wing and the downwash of the surrounding fluid as in equation (B 4) (figure 7a) [4], in order to better characterize the aerodynamics. Furthermore, to estimate the effects of the downwash on the lift, we use the effective angle of attack for the translational force term of the quasi-steady model (equation (B 3)). The correlation of the resulting quasi-steady lift using equation (B 1) to the Navier–Stokes solutions substantially improves compared with the estimation without the correction for the downwash (figure 7b,c). The improved agreement suggests that the drop in lift during the midstroke is indeed caused by the wing–wake interaction for the rigid wing. Quasi-steady analysis based on an effective angle of attack, composed of the downwash and the instantaneous translational velocity supports the lift enhancement via streamlining (figure 7). The flexible wing alleviates the downwash by streamlining its wing shape and outperforms its rigid counterpart with a difference in lift coefficient.

Effects of downwash illustrated with quasi-steady model predictions for the case shown in figure 6. (a) Angles of attack and (b,c) quasi-steady model predictions for the uncorrected angle of attack based on the wing's speed and orientation (black) and for the angles corrected with the downwash shown in figure 6e for the flexible (b) and the rigid (c) wings (blue). In (b,c), the lift calculations from the full Navier–Stokes (N–S) computations are indicated by the red lines.

4. Concluding remarks

The analysis presented in this paper represents a comprehensive assessment of the effects of structural flexibility on the aerodynamic performance of flapping wings. The Reynolds number of Re = 100 considered in this study is relevant to small insect flyers, such as fruit flies. The density ratio of ρ* = 7.8 is rather more relevant to aquatic swimmers, while for insect flyers ρ* = O(103). However, as long as the resulting wing motion, which is better characterized by the integrated parameter γ than by ρ* alone, remains the same, the resulting aerodynamics should remain the same [32]. Hence, the current aeroelastic analysis is relevant to insect flight.

However, there are caveats. This study only includes the role of chordwise flexibility and passive pitch in two-dimensional plunge motion. A real insect wing is intriguingly complex with complicated vein patterns, non-uniform anisotropic wing structures, three-dimensional wing kinematics with wing rotation at the wing root, etc. For example, the difference between the positive camber in a locust [14] and the negative camber in this study demonstrates that more complexity may need to be incorporated in the model to simulate insect aerodynamics.

Nonetheless such studies often result in sufficient information. This bioinspired, low Reynolds number, flexible flapping wing study highlights various aspects of the resulting FSI, applicable for the design of micro-air vehicles. By focusing on a canonical problem, the results are not limited to a specific flyer and are relevant to a wide range of low Reynolds number flyers.

The broadened scaling shows that the lift generation can be scaled with a priori known parameter γ that is a measure of the angle of attack due to passive pitch. Moreover, we have decomposed the resulting aerodynamics of the flexible wing into the nonlinear transient wing dynamics and the wing camber deformation by comparing with its rigid counterparts. Both components show enhanced performance for the flexible wing. In particular, the streamlining of the flexible wing can lead to lift enhancement by mitigating the lift-degrading wing–wake interaction. Finally, this study shows that optimal lift is obtained when the passive pitch is symmetric with the imposed translation. This observation is consistent with previously reported measurements of small insect flyers [21,22], which suggest that the current findings may be highly relevant to the understanding of how insects fly.

Acknowledgements

The research leading to these results has received funding from the Air Force Office of Scientific Research's Multidisciplinary University Research Initiative grant no. FA9550-07-1-0547.

Appendix A. Spatial and temporal sensitivity study

The fluid flow governed by equation (2.2) is solved on a body-fitted mesh around the flat plate as shown in figure 8. A sinusoidal plunge motion given by equation (2.3) is imposed at the LE of the flat plate. A free boundary condition is applied at the TE of the flat plate. The boundary conditions are the incompressible inlet with zero velocity at the outer boundary of the computational domain and noslip on the flat plate surface. The outer boundary is located at 63 chords away from the flat plate. Zero velocity and pressure are assigned as the initial condition.

Computational set-up. Computational domain around the flat plate (a) and imposed boundary conditions (b) for the fluid flow. The wing is placed at the centre of the computational mesh in (b). Grid (c) and time step (d) sensitivity studies for a plunging rigid flat plate at Re = 100 with reduced frequency of 0.33 for lift and drag. (c) Three resolutions are used for the grid sensitivities: 31 × 5 (red), 61 × 9 (black), 121 × 17 (blue) with 480 time steps per cycle. (d) On the 61 × 9 grids, the temporal resolution is varied by taking 240 (blue), 480 (black) and 960 (red) time steps per cycle.

To assess the grid sensitivity, in total four grids are employed with increasing numbers of cells on the chord and the edges being 31×5, 61×9, 121×17 and 241×33. The resulting time histories of lift coefficient using 480 time steps per motion period for flat plates at the reduced frequency of 0.33 are compared in figure 8. We consider the flow around a rigid flat plate, because the aerodynamic response for the parameters considered is more sensitive for a purely plunging rigid wing than for a flexible wing: the flow structures remain symmetric at first without generating lift, which changes abruptly after three cycles breaking the symmetry (figure 8) [43]. On the other hand, the symmetry is perturbed for flexible wings due to passive wing rotation. The mean and the r.m.s. values of the lift on the 61×9 and the 121×17 show converging trends (see also table 1). Based on these results, all computations presented in this study are performed on the 61×9 grid. To investigate temporal sensitivity, 240, 480 and 960 time steps are used on the 61×9 grid. Figure 8 also shows that the computations using 480 time steps per cycle are sufficient to obtain grid and time-step independent solution.

Appendix B. Quasi-steady model for hovering wing

The quasi-steady model [5] used in this study has three components asB1where the translational force component, which estimates the delayed stall, is computed asB2where Ut is the instantaneous translational velocity. The coefficient of lift [3] isB3where the angle of attack αi can either be the geometric angle of attack α or the effective angle of attack αe defined asB4where vd is the downwash upstream of the wing.

The rotational force term, which is only used for actively rotated rigid wings and represents the coupling term between the translational velocity and the angular velocity, isB5where Cr is an empirical coefficient that depends on the angular velocity. Here, we use Cr = 0.75 based on the suggestions by Sane & Dickinson [5].

Finally, the added mass force consists of two components, i.e.B6where the non-circulatory term can be written asB7where ah = −1 when the pivot point is located at LE [3]. The circulatory term isB8

. 2009Deformable wing kinematics in the desert locust: how and why do camber, twist and topography vary through the stroke?J. R. Soc. Interface6, 735–747.doi:10.1098/rsif.2008.0435 (doi:10.1098/rsif.2008.0435)