SUMMARY
In the present study, a computational investigation was carried out to understand the influence of flexibility on the aerodynamic performance of a hovering wing. A flexible, twodimensional, twolink model moving within a viscous fluid was considered. The Navier–Stokes equations governing the fluid dynamics were solved together with the equations governing the structural dynamics by using a strongly coupled fluid–structure interaction scheme. Harmonic kinematics was used to prescribe the motions of one of the links, thus effectively reducing the wing to a single degreeoffreedom oscillator. The wing's flexibility was characterized by the ratio of the flapping frequency to the natural frequency of the structure. Apart from the rigid case, different values of this frequency ratio (only in the range of 1/2 to 1/6) were considered at the Reynolds numbers of 75, 250 and 1000. It was found that flexibility can enhance aerodynamic performance and that the best performance is realized when the wing is excited by a nonlinear resonance at 1/3 of the natural frequency. Specifically, at Reynolds numbers of 75, 250 and 1000, the aerodynamic performance that is characterized by the ratio of lift coefficient to drag coefficient is respectively increased by 28%, 23% and 21% when compared with the corresponding ratios of a rigid wing driven with the same kinematics. For all Reynolds numbers, the lift generated per unit driving power is also enhanced in a similar manner. The wake capture mechanism is enhanced, due to a stronger flow around the wing at stroke reversal, resulting from a stronger end of stroke vortex at the trailing edge. The present study provides some clues about how flexibility affects the aerodynamic performance in low Reynolds number flapping flight. In addition, it points to the importance of considering nonlinear resonances for enhancing aerodynamic performance.
 flapping wing
 fluid–structure interactions
 finitedifference method
 wing flexibility
 nonlinear resonance
 low Reynolds numbers
INTRODUCTION
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 weightbalance 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. Twodimensional computations (e.g. Wang, 2000a; Wang, 2000b; Miller and Peskin, 2005; Miao and Ho, 2006) and threedimensional 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 costefficient manner, we limited ourselves to studies in two dimensions. A representative section of the wing (twodimensional foil) was used, and spanwise bending and torsion flexibility were discarded. A twolink 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 nonlinearities 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 nonlinearities 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 nonlinear 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.
MATERIALS AND METHODS
System description and computational formulation
As aforementioned, we considered a section of a threedimensional 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_{1}–c_{2} and c_{3}–c_{4}. 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).
System equations
The flapping motions of the chosen twodimensional configuration in a fluid are governed by a coupled system of equations describing the respective fluid and structural mechanics. The fluid dynamics is governed by the Navier–Stokes equations for an incompressible flow; that is: (1) (2) where t is time, x_{i} (i=1,2) is the spatial coordinate in the ith Cartesian direction, u_{i} is the corresponding velocity, p is the pressure and f_{i} is an external body force. The above equations have been made dimensionless by using the chord length of the undeformed plate, L_{c}, as the reference length scale, and the maximum translational velocity at the junction of the links, U_{c}, as the reference velocity scale. The Reynolds number is defined as Re=ρ_{f}L_{c}U_{c}/μ, where ρ_{f} and μ are the fluid density and viscosity, respectively.
The dimensionless form of the equations governing the motion of the structural system shown in Fig. 1 can be derived as: (3) where x(t), y(t) and θ(t) are, respectively, the joint horizontal motion, joint vertical motion and orientation angle of link B measured from an inertial reference frame as shown in Fig. 1A, andα (t) is the deflection angle between links A and B. Here, m_{i} is the total mass of the ith link (i=A,B), η_{i} is the distance from the junction to the center of mass of bar i as shown in Fig. 1A, and I_{i} is the moment of inertia of link i with respect to the hinge point b. Also, Q_{x} and Q_{y} are the fluid forces along the x and y directions, respectively, and Q_{θ} and Q_{α} are the fluid moments associated with the generalized coordinates θ(t) and α(t), respectively. The quantities g_{x}, g_{y}, g_{θ} and g_{α} are the corresponding contributions of centrifugal, elastic and gravitational forces. The reference length scale and velocity defined above, together with the fluid density are used to make Eqn 3 dimensionless. The fluid forces and moments are determined from Eqns 1 and 2.
In the numerical experiments conducted in this study, the translational motions of the junction as well as the orientation of link B are prescribed. With these prescribed motions, the four degrees of freedom of the system can be effectively reduced to one; that is, the deflection angleα (t) between plates A and B. Thus, the overall deformation of the wing section is determined by the deflection angleα (t), which is governed by the following reduced form of Eqn 3: (4) Eqn 4 resembles the equation governing a harmonic oscillator with forcing due to the prescribed kinematics and the fluid forces (e.g. Nayfeh and Balachandran, 1995). The nonlinearities arise from the sin(θ+α) term due to the kinematics and the fluid forcing. For this particular study, we only take into account the fluid damping which arises through the fluid moment Q_{α}. It should be noted that selecting a proper structural damping model is far from trivial, and this is an active research topic in structural biomechanics. Damping models for insect wings are relatively few [e.g. classical viscous damping model used by Herbert (Herbert, 2002) and the viscoelastic model used by Bao and colleagues (Bao et al., 2006)] and the existing models require a fair amount of empirical information.
Prescribed kinematics, parameter values and computational formulation
To prescribe the translational motions of the junction and the orientation of link B, we define the states x(t), y(t) and θ(t) as: (5) where A_{0} is the stroke length of the pin point,θ _{0} is the mean orientation angle for link B, γ is the rotation amplitude, ω_{f} is the frequency of the prescribing or forcing oscillation and ϕ is the phase angle between x(t) and θ(t). The exponential terms were used in order to reduce transient effects (Combes and Daniel, 2003). The time constant was chosen as τ=1.6×π/ω_{f} because 99.8% of the prescribed amplitude was reached after a time length of 5 periods. The following parameters corresponding to symmetric hovering were selected (Wang et al., 2004): (6)
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, embeddedboundary 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 nonslip 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).
RESULTS
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).
Computational setup
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 farfield boundaries. A nearuniform grid zone was generated near the center, where the motions of the twolink 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 y directions, 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 wingsection 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 timeaveraged quantities were computed over the last 7 periods, 13 periods and 10 periods, for Re=75, 250 and 1000, respectively.
Aerodynamic quantities
For the flexible wings considered in this study, the lift and drag coefficients were defined as: (7) where the quantities and are dimensional quantities and Q_{x}(t) and Q_{y}(t) are nondimensional quantities. Once the equation for the deformationα (t) is solved at each iteration, the driving forces in the prescribed generalized coordinates x(t), y(t) and θ(t) are computed from Eqn 3. The total power input, which is the sum of horizontal translational power and rotational power, can be computed from: (8) where P_{tr} and P_{rot} are the translational and rotational power inputs at the hinge b, R_{x}(t) is the driving force in the x direction and R_{θ}(t) is the driving moment in the θ(t) angular direction. In an ideal case were the driving mechanism is perfectly elastic and the negative power provided to the mechanism can be stored as potential energy for later use, this power will enhance the wing's aerodynamic efficiency. Here, following Berman and Wang (Berman and Wang, 2007), a conservative approach was used and it was assumed that the negative power was not available for reuse; that is: (9) The power coefficient is defined as: (10) where and are dimensional quantities and P_{tr}(t) and P_{rot}(t) are nondimensional quanitites.
In Fig. 3, the variations of the mean values of C_{L} and C_{D} and the aerodynamic performance ratios C_{L}/C_{D} and 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 liftforce 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 nonlinear 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 timeaveraged 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 twopeak lift curve is observed when ω_{f}/ω_{n}=1/3.
Vortex structures
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 phaseaveraged 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 , and the lower link A rotates with an angular speed . The added angular speed affects the overall dynamics at stroke reversal. First, the camber generated by the angular deformation α(t) at the end of the stroke (see Fig. 5B) reorients the zero lift direction on the wing and enhances wake capture effects. This enhancement mechanism is analogous to the one produced by orientation advancement in rigid wings (e.g. Wang et al., 2004). It is important to note that an excessive degree of flexibility (low frequency ratios) produces a large camber at stroke reversal, which has a negative effect on the lift production (e.g. atω _{f}/ω_{n}=1/2 in Fig. 4). The evolution and strength of the LEV on the other hand (see Fig. 5A) is only a weak function of the wing's flexibility. The formation time as well as the maximum circulation shown in Fig. 5L are approximately the same for the rigid and the flexible wings.
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 counterrotating vortices together with the LEV. This vortex pair generates flow directed towards the wing, enhancing the wakecapturing 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 lowpressure zone, which translates into increasing forces up to the pinchoff 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 abovementioned 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 obliqueshaped 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).
DISCUSSION
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 threedimensional 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 twodimensional 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 twolink 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 nonlinear resonances atω _{f}=1/3ω_{n} andω _{f}=3ω_{n}. These resonances are expected in systems with cubic nonlinearities; 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 nonlinear superharmonic resonance ω_{f}=1/3ω_{n} was found to be beneficial for aerodynamic performance. Inclusion of fluid effects will give rise to quadratic nonlinearities and additional nonlinear 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 quasisteady 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 twodimensional computations were compared with the quasisteady 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 semiempirical 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 quasisteady contributions is more complex as additional new states have been included.
Conclusions
In the present work, the influence of flexibility on the aerodynamic performance of a twodimensional 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 degreeoffreedom nonlinear 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 nonlinear 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 nonlinear resonance and would one achieve a better performance with a nonlinear 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.
FOOTNOTES

The authors gratefully acknowledge the support received through AFOSR Grant No. FA95500610093 and ARO Grant No. W911NF0610369. Support received from the Minta Martin Foundation is also acknowledged. We are thankful to Mr N. Beratlis for providing the vortex circulation computing algorithm.