Influence of flexibility on the aerodynamic performance of a hovering wing

Over the past decade, insect flight has attracted a lot of interest in a variety of disciplines in science and engineering. As a result, many experimental investigations (e.g. Ellington et al., 1996; Dickinson et al., 1999) as well as computational investigations (e.g. Liu and Kawachi, 1998; Sun and Tang, 2002; Ramamurti and Sandberg, 2002; Ramamurti and Sandberg, 2006; Wang,2000a; Wang,2000b; Wang et al.,2004) have been reported in the literature. The aim of most of these studies has been to understand the complex, unsteady mechanisms that enable the generation of aerodynamic forces for hovering and maneuvering. Insect wings are complex structures that during flapping undergo deformations due to aerodynamic forces and elastic forces, as well as inertial forces due to the accelerations experienced by the system mass. The wing structural behavior depends, to a large extent, on the internal distribution of compliant components and mechanisms (Wootton,1999). It is important to note that insect wings lack internal muscles and, hence, there are no actuators to realize internal control forces(Wootton et al., 2003).

In a variety of species, the roles of inertial, elastic and aerodynamic forces during flapping flight have been the focus of many investigations (see,for example, Ellington, 1984b; Ennos, 1989; Lehman and Dickinson, 1997; Sun and Tang,2002; Daniel and Combes,2002; Combes and Daniel,2003; Song et al.,2001). It is difficult to make direct comparisons between the different studies, not only because the studies usually involve different species but also because different approaches have been used to compute the forces. For example, Combes and Daniel assessed the relative contributions of aerodynamic, inertial and elastic forces to the wing deformation of the Manduca sexta hawkmoth (Combes and Daniel, 2003). They concluded that the wing motion of this particular insect is mostly determined by the wing inertial and elastic forces with the aerodynamic loads providing damping. During hovering, the typical ratio of wing inertial force to aerodynamic force was found to be about seven. This result was obtained by using scaling arguments and assuming a weight balance to get a fluid–force estimate. In other species, this ratio has been found to be much lower. Ennos, for example, showed that for several species of Diptera, the magnitudes of inertial bending moments are about twice the aerodynamic moments during harmonic flapping(Ennos, 1989). Also, in this case, the analysis was based on the weight-balance assumption and harmonic kinematics. However, unlike Combes and Daniel(Combes and Daniel, 2003),Ennos considered the effect of the virtual or added mass of the surrounding fluid. It should be noted that in the studies of Ennos(Ennos, 1989) and Combes and Daniel (Combes and Daniel,2003), the aerodynamic forces were underestimated, since the drag component of the fluid force was neglected.

With increases in computational power, computational models of insect flight have become more sophisticated. Two-dimensional computations (e.g. Wang, 2000a; Wang, 2000b; Miller and Peskin, 2005; Miao and Ho, 2006) and three-dimensional computations (e.g. Liu and Kawachi, 1998; Ramamurti and Sandberg, 2002; Ramamurti and Sandberg, 2006; Sun and Tang, 2002) with various degrees of complexity have been reported. In general, for low ratios of inertial to aerodynamic forces, one can expect complex aeroelastic interactions to occur. An interesting open question within this context is the following: how does structural flexibility affect the aerodynamic performance of a given flapping wing and what is the effect of the Reynolds number? The present study attempted to address this question by using computational investigations. To the best of our knowledge, no prior computational studies addressing this question have been carried out.

In the present study, in order to explore a fairly wide parametric regime in a cost-efficient manner, we limited ourselves to studies in two dimensions. A representative section of the wing (two-dimensional foil) was used, and spanwise bending and torsion flexibility were discarded. A two-link structure connected with a torsion spring was used to account for deformation in the chordwise direction. This system has four degrees of freedom, which are effectively reduced to one by prescribing harmonic hovering motions of one of the links. The links were considered to be rigid in the present work, and they are currently being extended to flexible beams in ongoing efforts. The large angular deformations of the links gave rise to cubic and higher order odd non-linearities in the governing equations like those seen in equations governing a pendulum as well as flexible beams (e.g. Anderson et al., 1994). In a sense, one could consider the two links as a double pendulum with a torsion spring. Fluid non-linearities were also considered here. Different values of the torsion spring stiffness were considered at the Reynolds numbers of 75,250 and 1000, and the results obtained are reported in the form of mean lift force, mean drag force, ratio of lift to drag, and ratio of mean lift coefficient to total power input. The performance of the hovering wing when excited at a non-linear resonance of the structural system was also examined.

In the following section, a description of the system is provided along with the computational formulation. Then, the parametric space and system kinematics are detailed. Next, Results and Discussion sections follow, with a Conclusions section at the end.

As aforementioned, we considered a section of a three-dimensional wing and accounted for chordwise deformations, but not for spanwise bending and torsion flexibility. As shown in Fig. 1, the structural system under consideration consists of two rigid links A and B, which are joined at the center b by a pin that contains a linear torsion spring. In this model, flexibility is concentrated at one discrete location of the structural system, and inclusion of elastic links will allow chordwise variations of stiffness and mass to be accounted for. For computational purposes, the two links are covered by a set of aerodynamic surfaces that define the boundary between the airfoil and the fluid, and deform as the angle between the two links changes. The aerodynamic surfaces consist of two rigid segments, R_Sa and R_Sb (see Fig. 1B), and two segments that dynamically deform according to the angle between the two links. The deformation is prescribed by fitting the Hermite interpolation polynomials c₁–c₂ and c₃–c₄. We have found that this modeling is robust and helps maintain the smoothness of the surface even for large values of the angle between the plates (large deformation configurations).

Based on the adopted normalization, the problem is completely defined by the density ratio ρ_b/ρ_f, the frequency ratio ω_f/ω_n and the Reynolds number Re. Here, ρ_b is the density of the wing's material, and ω_n=√(k/I_A)is the linear natural frequency of the oscillator(Eqn 4). The frequency ratioω _f/ω_n is used to characterize the flexibility of the wing section.

Three Reynolds numbers were considered (Re=75, 250 and 1000) to investigate the effect of the reduction in viscous dissipation on the system dynamics. The mass ratio was set toρ _A/ρ_f=25, as this value provided a ratio close to 2 for the maximum translational inertial force over maximum drag force at Re=75 for the chosen geometry and kinematics. The above ratio was determined through numerical experiments with the rigid wing. To compute the maximum horizontal translational inertial force, the total wing mass was multiplied by the maximum acceleration determined from the second derivative of x(t) in Eqn 5. The value of peak drag force, on the other hand, was obtained from the rigid wing simulation at Re=75. The wing has a thickness of 10% of the undeformed chord length and circularly formed edges. For the simulations conducted at Re=75, the frequency ratioω _f/ω_n was set to 1/2, 1/2.5, 1/3, 1/3.5, 1/4 and 1/6. For Re=250 and 1000, this ratio was set to 1/2, 1/3, 1/4 and 1/6. The resulting range of maximum deflection angles varied from 10 to 70 deg. Also, the rigid wing problem (no angular deformation between the links)was run for all of these Reynolds numbers. It should be noted that for frequency ratios below 1/2, the computations would fail since the two plates collide during rotation. This limitation arises from the fact that the flexibility in the present model is concentrated at the hinge point and the distributed chordwise variations of stiffness and mass are not accounted for. A flexible beam model and/or inclusion of structural damping may help to address this issue and enable computations with frequency ratios of about one.

Eqns 1, 2 and 4 governing the dynamics of the fluid–structure system are numerically solved by using a strongly coupled, embedded-boundary formulation. The overall approach is a mixed Lagrangian–Eulerian formulation, where Eqns 1 and 2 governing the fluid flow are solved on a fixed Cartesian grid, which is not aligned with the wing surface,and the non-slip conditions are enforced via local reconstructions of the solution near the solid interface (see, for example, Balaras, 2004; Uhlmann, 2005; Yang and Balaras, 2006). The fluid and the structure are treated as elements of a single dynamical system,and all governing equations are integrated simultaneously and interactively in the time domain by using a predictor–corrector scheme. Further details on the coupling scheme and the overall fluid–structure interaction algorithm can be found in the work of Yang and colleagues(Yang et al., 2008).

In this section, computations of aerodynamic forces are presented for three different Reynolds numbers and different wing flexibility values. Comparisons are also made between the results obtained for the rigid wing and flexible wing cases and with the results obtained by Wang and colleagues(Wang et al., 2004).

The computational grid was carefully selected to resolve the thin boundary layers and detached shear layers on the moving links and the wake vortical structures for the different Reynolds numbers considered in this study. The rigid wing was set to move in the center of a box with dimensions of 30L_c×30L_c, in order to minimize interference effects from the far-field boundaries. A near-uniform grid zone was generated near the center, where the motions of the two-link system took place, and this zone was stretched towards the boundaries. For the Re=75 simulations, the uniform grid zone had a local cell size ofΔ x=Δy=0.0038L_c, and the total number of points was 1229×551 along the x and ydirections, respectively. Through grid refinement studies, we found that the above resolution was sufficient to capture all flow features. In Fig. 2, computationally obtained aerodynamic forces are shown for approximately half the resolution throughout the computational domain (total number of points was 664×400)and the same forcing conditions (i.e. τ=0 in Eqn 5) and Reynolds number as that for the baseline rigid wing computation. The corresponding lift and drag coefficients determined in the computations of Wang and colleagues(Wang et al., 2004), where a hovering ellipse with the same kinematics is considered instead of a plate,are also included in the figure. The agreement between the results obtained with the two different grids is good, with a maximum difference of around 3%. Despite the differences in the wing-section shapes, after the initial transients (t/T>2), good agreement is also seen with the results obtained by Wang and colleagues(Wang et al., 2004).

For the results of Fig. 2,within the boundary layers on the link surface, we estimated the number of grid points to be approximately 8 and 16 for the coarse and fine grids,respectively. As the Reynolds number increases, the boundary layer thickness is expected to decrease as 1/√Re, and in order to keep the resolution within the above range, a grid with 1320×1038 points was found to be sufficient for both Re=250 and Re=1000 cases. All the results presented in this article were obtained with a 1229×551 grid for the Re=75 simulations and a 1320×1038 grid for the Re=250 and Re=1000 simulations. The governing equations were integrated for a time length of 14 periods, 21 periods and 15 periods, for Re=75, 250 and 1000, respectively. The time-averaged quantities were computed over the last 7 periods, 13 periods and 10 periods, for Re=75, 250 and 1000, respectively.

In Fig. 3, the variations of the mean values of C_L and C_D and the aerodynamic performance ratios C_L/C_Dand C_L/C_PW with respect to the frequency parameter ω_f/ω_n are shown. The rigid wing has zero torsion stiffness or equivalently ω_n=0. For all cases, the lift and drag coefficients exhibit a peak at a frequency ratio of ω_f/ω_n=1/3. The performance ratio C_L/C_D also exhibits a prominent peak at this frequency ratio. For Re=75, 250 and 1000, increases of about 28%, 23% and 21% over those obtained for the rigid wing were observed,respectively. The variations of the aerodynamic quantities with respect to the frequency ratio show similar characteristics for all three Reynolds numbers. However, it is interesting to note the striking difference between the graph of C_L/C_D obtained for the Re=75 case and those obtained for the higher Reynolds numbers. For the lowest Reynolds number and ω_f/ω_n=1/4,the above ratio is over 13% higher than that obtained for the rigid wing,while for Re=250 it is only increased by 0.5%. In Fig. 3C, it can be seen that for ω_f/ω_n=1/3, the performance ratio C_L/C_PW is 39% and 28% higher than that obtained for the rigid wing for the Re=75 and Re=250 cases,respectively. Interestingly, this measure is only about 13% higher than that obtained for the rigid case at Re=1000.

In Fig. 4, the time histories of the lift and drag coefficients are shown for all Reynolds numbers at three selected frequency ratios. The effects of flexibility are noticeable on the lift-force peaks at the initiation of the stroke (indicated with a black arrow in Fig. 4A). For Re=75 and ω_f/ω_n=1/2, corresponding to the most flexible foil, this peak is negative, while for the rigid case the coefficient of lift peaks at 0.5. For all cases in between, the enhancement of the mean lift force seen in Fig. 3 comes from the gradual increase of this peak, which is at 0.83 and 1.28 for ω_f/ω_n=1/4 andω _f/ω_n=1/3, respectively. For the latter frequency ratio, where a structural non-linear resonance occurs, the lift peak is also delayed and nearly coincides with the translational lift peak. This is translated into a larger area under the lift curve per period and a larger time-averaged C_L value. The temporal variations of the lift and drag coefficients for Re=250 and 1000 are more complex than the variations seen at Re=75, and the periodicity is practically lost. Still, in an average sense, negative lift peaks after stroke reversal and larger translational lift peaks are seen whenω _f/ω_n=1/2. Also, a widened two-peak lift curve is observed when ω_f/ω_n=1/3.

In order to relate the temporal variations of the lift and drag forces to specific flow structures, we carefully examined several realizations of the instantaneous flow fields. In Fig. 5, vorticity isolines are shown for the rigid andω _f/ω_n=1/3 cases at Re=75. For clarity, the lift coefficient variation has been added(Fig. 5K), together with the temporal variation of the phase-averaged circulation of the most important vortical structures generated during a flapping cycle. These are the leading edge vortex (LEV) shown in Fig. 5A, the end of stroke vortex (ESV) shown in Fig. 5C and trailing edge vortex (TEV) shown in Fig. 5E. The circulation of each of these vortices was computed as a function of time by direct integration of the vorticity within a given threshold contour around each vortex. The selection of the threshold contour, although arbitrary, was consistently taken to be the lowest closed vorticity isoline in the vicinity of the given vortex.

As the flexible wing approaches the end of the stroke in Fig. 5A, it exhibits different rotation velocities on the two components A and B. The driven link B rotates with the prescribed angular velocity

Another effect of the higher rotation speeds at the trailing edge for the flexible wings is the formation of a stronger shear layer, which rolls up into a stronger ESV (see Fig. 5C). On examining the ESV circulation plots(Fig. 5K), one finds that the strength and life span are significantly enhanced when compared with those of a rigid wing. The ESV pinches off later, forming a pair of counter-rotating vortices together with the LEV. This vortex pair generates flow directed towards the wing, enhancing the wake-capturing effects. This is more clearly reflected in the lift coefficient evolution shown in Fig. 5K. In contrast to the rigid wing, where the lift curve reaches a maximum (point H) and starts to decrease, for the flexible wing the production of increased lift continues for longer (point C).

Once a flexible wing's deflection has reached its maximum, the elastic energy stored in the torsion spring is released to generate a restoring motion, whose timing again depends on the degree of flexibility of the structure (see Fig. 5D and 5E). This restoring motion produces a dynamical change in the wing's camber with a resulting increase in the fluid forces. This also affects the formation and growth of the TEV. It is well established that this flow structure generates a low-pressure zone, which translates into increasing forces up to the pinch-off time (see, for example, Wang,2000b). The time at which the TEV pinches off is correlated with the translational force peak; this peak is advanced in theω _f/ω_n=1/3 case when compared with that for the rigid wing.

In Fig. 6, a quantitative comparison of the LEV, ESV and TEV dynamics is shown for different flexibilities at Re=75 by determining their average circulations as a function of time. The maximum averaged circulation for each vortex and the time at which it occurs with respect to the stroke reversal are provided in Table 1. As expected, from what was observed in Fig. 5, the LEV dynamics is similar for all frequency ratios in terms of both strength and timing. The TEV on the other hand, attains a higher maximum circulation as the wing becomes more flexible. However, the time it takes to reach this maximum circulation is shortest for ω_f/ω_n=1/3, where the best aerodynamic performance is seen. For the ESV vortex, the maximum circulation increases for the frequency ratiosω _f/ω_n=1/3 and 1/4. The peak circulation forω _f/ω_n=1/3 is 20% lower than that obtained for ω_f/ω_n=1/4, but the deflection at stroke reversal in terms of the maximum deformation angle is 90% larger whenω _f/ω_n=1/3 (56 deg. forω _f/ω_n=1/3 and 29 deg. forω _f/ω_n=1/4). This is translated into a larger projected area contributing to the lift force. Also, as seen from Fig. 6, the time delay of the peak circulation of the ESV vortex is increased as the wing becomes more flexible.

A more direct illustration of the above-mentioned vortex evolutions is given in Fig. 7, where instantaneous vorticity isolines are shown for eight characteristic instances during a flapping cycle. For ω_f/ω_n=1/3 and 1/4, it is clear that the enhanced ESV vortices produce an oblique-shaped TEV vortex. For ω_f/ω_n=1/2, an excessive negative camber is produced at stroke reversal, which then generates a high suction zone on the lower side of the wing leading to the negative peak in the C_L curve seen in Fig. 4. The Reynolds number effects on the temporal evolution of the lift and drag forces seen in Fig. 4 can also be observed in Fig. 8, where the instantaneous vorticity isolines are shown for Re=250 for all frequency ratios. Clearly, as the viscous damping is decreased, the system dynamics system ceases to be periodic and the important vortices are stronger and are not dissipated as quickly as seen in the Re=75 case (see Fig. 7). For the case of a rigid wing, for example, the LEV from a given stroke interacts with the shear layer being generated in the next stroke, and this induces a premature formation of the new LEV. This process is not periodic, which is also reflected in the evolution of lift and drag forces. A similar interaction is observed atω _f/ω_n=1/2 (see last three frames in Fig. 8B).

Insect wings are flexible structures that undergo large displacements and deformations during flapping, as the wing structures interact with the surrounding flow. There has been speculation that many insects flap their wings at frequencies close to the natural frequency of the structure. For example, analysis of Manduca sexta wings (see Combes and Daniel, 2003; Wootton et al., 2003) has shown that the wing's first natural frequency is close to the driving frequency in normal flapping motion. This suggests that insects may be taking advantage of a structural resonance to reduce energy consumption and enhance aerodynamic performance. Despite the significance of such a hypothesis, only a limited number of studies have addressed the problem due to the exceedingly complex fluid–structure interactions that are encountered in experimental or numerical work. Although three-dimensional computations of the Navier–Stokes system coupled with a wing structural system are within the reach of today's computers, one still needs to develop appropriate mathematical models and tools to capture all important phenomena in this complex system. In this regard, the present study extends computational work that has been conducted before with simplified two-dimensional rigid wings to include the effects of flexibility. The wing is represented by two rigid links, which are joined at the center by a pin that contains a torsion spring. The kinematics of one of the links is prescribed, while the motion of the other link is determined by the fluid–structure interactions. Although a fairly wide range of Reynolds numbers and frequency ratios was examined, we found that the computations would fail at forcing frequencies close to the linear resonance frequency due to an excessive degree of flexibility. This limitation is due to the concentration of flexibility at a discrete point, and replacement of rigid links with elastic links modeled as elastic beams is expected to help in overcoming this limitation.

As mentioned earlier, the two-link structural system can be perceived as a double pendulum with a common hinge. In particular, when one of the link motions is prescribed, the other link behaves as a pendulum subjected to a constraint arising from the prescribed motion and complex fluid–structure interactions. Eqn 4 does resemble the equation of a pendulum driven in a fluid. A straightforward perturbation analysis (e.g. Nayfeh and Balachandran, 1995)shows that the structural system can exhibit non-linear resonances atω _f=1/3ω_n andω _f=3ω_n. These resonances are expected in systems with cubic non-linearities; for example, in the equations governing local oscillations of a pendulum about an equilibrium position and elastic systems such as beams (e.g. Anderson et al., 1994). Operation of the flexible wing at the non-linear superharmonic resonance ω_f=1/3ω_n was found to be beneficial for aerodynamic performance. Inclusion of fluid effects will give rise to quadratic non-linearities and additional non-linear resonances.

For the specific set of kinematics that we considered, most of the benefits of having a flexible wing are associated with the stroke reversal phase of the cycle. Especially for the optimal flexibility cases(ω_f=1/3ω_n), the strength and timing of the ESV, as well as the dynamical changes of the wing's camber due to structural deformations, are responsible for the performance enhancement. The overall enhancement mechanism is analogous to that produced by orientation advancement in rigid wings (see, for example, Wang et al., 2004). It is noted that the present computations cover a wide range of frequency ratios and, consequently, wing deflections range from a few degrees to very large values. For example, in the case of highly stiff wings(e.g. ω_f/ω_n=1/6), the maximum deflection angle between the links was about 11, 13 and 16 deg. for Re=75, 250 and 1000, respectively, while for highly flexible wings (e.g.ω _f/ω_n=1/2) the corresponding numbers were 67, 68 and 91 deg., respectively. In insects, the wing deformation magnitudes increase as the body size and mass increase, and it is conceivable that deformations seen in this study at the aerodynamically preferred frequency ratio of ω_f/ω_n=1/3 could be possible in some species. On the other hand, for small insects such as Drosophila only small magnitude wing deformations have been observed. The computations of this study show that as the wing is made stiffer, the performance enhancements are marginal when compared with a rigid foil. For example, at Re=75 and the highly stiff case ω_f/ω_n=1/6, C_L/C_D is approximately 6% higher than that of a rigid foil. For the higher Reynolds numbers Re=250 and 1000 there is actually no enhancement, and the performance is worse than that obtained with a rigid foil. The above results indicate that low Reynolds number regimes might benefit performance even at small chordwise distortions.

The force histories, in particular for the low Reynolds numbers, appear to reach a periodic steady state after the initial transients for all of the frequency ratios that were considered, suggesting that quasi-steady models might be able to reproduce this behavior. Such models have been reported in the literature, and have been adapted for flapping flight based on models developed for high Reynolds number fixed wing aeroelasticity studies by including wing rotation along with translation(Ellington, 1984a; Ellington, 1999), and forces due to added mass (Sane and Dickinson,2002). In a more recent study by Wang and colleagues(Wang et al., 2004) the unsteady forces from experiments and with two-dimensional computations were compared with the quasi-steady model predictions. They pointed out that the force predictions, which were made by using models based on potential flow theory (Munk, 1925) for a constant pitching amplitude and constant translating speed wing, deviated substantially from the experimentally determined unsteady forces. The force predictions from a semi-empirical model based on numerical results from steady translating wings at a fixed angle of attack were in broad qualitative agreement with the unsteady forces. However, detailed comparisons revealed that, depending on the kinematics, the unsteady effects can reduce the total lift by a factor of 2 to 3. In the present case, due to the wing's flexibility, the identification of the quasi-steady contributions is more complex as additional new states have been included.

In the present work, the influence of flexibility on the aerodynamic performance of a two-dimensional hovering wing section was numerically studied. The wing model consists of two rigid links that are joined at the center with a linear torsion spring. By prescribing the kinematics of the top link, the structural system is effectively reduced to a single degree-of-freedom non-linear oscillator. The viscous flow around this structure is described by the incompressible form of Navier–Stokes equations. The combined set of equations describing the fluid and structural dynamics are integrated in time by using a scheme that can capture strongly coupled fluid–structure interactions.

The results obtained in this study demonstrate that flexibility can be beneficial in terms of enhancing aerodynamic performance. Furthermore, we found that in the frequency range below the first natural frequency, the best performance is achieved when the wing is driven at a frequency close to one of the non-linear resonances (a superharmonic resonance of order three) of the system. This behavior is common to all of the Reynolds numbers studied. In terms of the flow physics, the wake capture mechanism is enhanced partially due to a stronger flow around the wing at stroke reversal. However, it needs to be noted that the cases where the wing is driven at or close to the first natural frequency of the system were not considered in this study, and it is possible that a better aerodynamic performance may be achieved with a linear resonance and this remains to be explored.

The study also leads to the following open questions. (i) Why is there a performance enhancement when the system is excited at a flapping wing's non-linear resonance and would one achieve a better performance with a non-linear resonance compared with a linear resonance? (ii) Which kinematics is preferable from an aerodynamic efficiency standpoint? The interplay between wing flexibility and kinematics together with qualitative changes (and bifurcations) in the system dynamics as a function of the Re number require further investigation.