## SUMMARY

Of the insects that have been filmed in flight, those that are 1 mm in length or less often clap their wings together at the end of each upstroke and fling them apart at the beginning of each downstroke. This `clap and fling'motion is thought to augment the lift forces generated during flight. What has not been highlighted in previous work is that very large forces are required to clap the wings together and to fling the wings apart at the low Reynolds numbers relevant to these tiny insects. In this paper, we use the immersed boundary method to simulate clap and fling in rigid and flexible wings. We find that the drag forces generated during fling with rigid wings can be up to 10 times larger than what would be produced without the effects of wing–wing interaction. As the horizontal components of the forces generated during the end of the upstroke and beginning of the downstroke cancel as a result of the motion of the two wings, these forces cannot be used to generate thrust. As a result, clap and fling appears to be rather inefficient for the smallest flying insects. We also add flexibility to the wings and find that the maximum drag force generated during the fling can be reduced by about 50%. In some instances, the net lift forces generated are also improved relative to the rigid wing case.

## INTRODUCTION

A vast body of research has described the complexity of flight in insects ranging from the fruit fly, *Drosophila melanogaster*, to the hawk moth, *Manduca sexta* (Sane,2003; Wang, 2005). Over this range of scales, flight aerodynamics as well as the relative lift and drag forces generated are surprisingly similar(Birch et al., 2004; Hedrick et al., 2009). The smallest flying insects have received far less attention,although previous work has shown that flight kinematics and aerodynamics may be significantly different (Miller and Peskin, 2004; Miller and Peskin, 2005; Sunada et al., 2002; Weis-Fogh, 1973). These insects are on the order of 1 mm in length or smaller, fly at Reynolds numbers(*Re*) near 10 or below and include ecologically and agriculturally important species such as parasitoid wasps, thrips and haplothrips. At these low *Re*, lift forces relative to drag forces decrease substantially(Miller and Peskin, 2004; Wang,2000).

While quantitative data on wing beat kinematics is not available, it is thought that many tiny insects clap their wings together at the end of the upstroke and fling them apart at the beginning of the downstroke(Ellington, 1999) (see supplementary material Fig. S1). In fact, all tiny insects that have been filmed to date appear to use `clap and fling'. Weis-Fogh first described how the tiny wasp *Encarsia formosa* uses this motion and speculated that it augments the lift forces generated by strengthening the bound vortex on each wing (Weis-Fogh, 1973). The clap and fling motion has also been reported in the greenhouse whitefly *Trialeurodes vaporariorum*(Weis-Fogh, 1975) and *Thrips physapus* (Ellington,1984a). Miller and Hedrick have observed clap and fling kinematics in high-speed videos of the tiny parasitoid wasps *Muscidifurax raptor*and the jewel wasp *Nasonia vitripennis* (see supplementary materialFig. S2). Fling has also been observed in a few medium and larger insects such as butterflies and moths (Marden,1987) and the tethered flight of *Drosophila melanogaster*(Vogel, 1967).

A combination of experimental, theoretical and numerical work supports that insects augment the lift forces generated during flight by the clap and fling mechanism (Ellington, 1984a; Ellington, 1984b; Lehmann and Pick, 2007; Lehmann et al., 2005; Lighthill, 1973; Miller and Peskin, 2005; Weis-Fogh, 1973). Lift is augmented during wing rotation (fling) and subsequent translation by the formation of two large leading edge vortices(Miller and Peskin, 2005; Sun and Xin, 2003). Much less attention has been given to the total force required to actually clap the wings together and fling the wings apart. Miller and Peskin found that relatively large forces are required to fling rigid wings apart at *Re*below 20 (Miller and Peskin,2005). Furthermore, the lift to drag ratios produced during clap and fling are lower than the ratios for the corresponding one wing case,although the absolute lift forces generated are larger. For the smallest insects, the forces required to fling the wings are so large that it raises the question of why tiny insects clap and fling in the first place.

One major assumption in previous clap and fling studies using physical models (Lehmann and Pick,2007; Lehmann et al.,2005; Maxworthy,1979; Spedding and Maxworthy,1986; Sunada et al.,1993), mathematical models(Lighthill, 1973) and numerical simulations (Chang and Sohn,2006; Miller and Peskin,2005; Sun and Xin,2003) is that the wings are rigid. It seems likely that wing flexibility could allow the wings to reconfigure to lower drag profiles during clap and fling. In this case, the fling might appear more like a `peel', and the clap might be thought of as a reverse peel(Ellington, 1984a; Ellington, 1984b). A number of studies support the idea that flexibility allows for reconfiguration of biological structures, which results in reduced drag forces experienced by the organisms (e.g. Alben et al.,2002; Alben et al.,2004; Denny, 1994; Etnier and Vogel, 2000; Koehl, 1984; Vogel, 1989). The basic idea in these cases is that the force on the body produced by the moving fluid causes the flexible body to bend or reconfigure, which, in turn, reduces the force felt by the body.

A number of recent computational and experimental studies have explored the role of wing flexibility in augmenting aerodynamic performance in insects but most of this work has focused on single wings at `higher' *Re*(>75). For example, Vanella et al. explored the influence of flexibility on the lift to drag ratio for *Re* ranging from 75 to 1000 using two-dimensional, two-link models of a single wing(Vanella et al., 2009). They found that aerodynamic performance was enhanced when the wing was flapped at a 1/3 of the natural frequency. Ishihara et al. studied passive pitching by modeling a rigid wing that was free to pitch and were able to generate sinusoidal motions that produced enough lift to support some Diptera(Ishihara et al., 2009). A number of other studies have explored the role of wing flexibility in avian flight (Heathcote et al.,2008; Kim et al.,2008) and thrust generation(Alben, 2008; Lauder et al., 2006; Mittal, 2006).

*Drosophila virilis*(Vogel, 1967) and Lepidoptera (Norberg, 1972). The wings of the tiny wasps

*Muscidifurax raptor*and

*Nasonia vitripennis*definitely flex during flight but it is difficult to determine the extent of peeling from available videos. This motion can be approximated as a `flat peel'. Ellington(Ellington, 1984b) claims that if the wings `unzip' with a velocity

*u*(

_{z}*t*), then the circulation Γ around the separated section of the wing may be approximated as:

*x*is the length of the exposed section of the wing,and β is the constant half-angle between the wings and

_{e}*f*(β)≈(β –π/2)

^{2}+2,0≤β≤π. According to this model, Ellington(Ellington, 1984b) calculated that the circulation generated by the flat peel is 2.6 times greater than the corresponding case of inviscid fling calculated by Lighthill(Lighthill, 1973).

*c*

^{2}. This air mass gains velocity

*u*in the duration

_{c}*c*/

*u*during the clap. The mean lift force (

*F̄*) per unit span will be approximately equal to the change in momentum of the air mass:

In this paper, we have used computational fluid dynamics to study the effects of wing flexibility on the forces produced during clap and fling. The immersed boundary method was used to model pairs of rigid and flexible wings performing a two-dimensional clap and fling stroke at *Re*=10. The cases of flexible clap and peel were compared with the cases of rigid clap and fling. Lift and drag coefficients were calculated as functions of time and related to the bending of the wing and the relative strengths of the leading and trailing edge vortices.

## MATERIALS AND METHODS

As both viscosity and the interactions of the flexible wing with the air are important at this scale, a direct numerical simulation of the fully coupled fluid–structure interaction problem is appropriate. In this paper, the immersed boundary method(Peskin, 2002; Zhu and Peskin, 2002) is used to simulate two flexible wings immersed in a viscous, incompressible fluid. The immersed boundary method has been used successfully to model a variety of problems in biological fluid dynamics. Such problems usually involve the interactions between incompressible viscous fluids and deformable elastic boundaries. Some examples of biological problems that have been studied with the immersed boundary method include aquatic animal locomotion(Fauci, 1990; Fauci and Fogelson, 1993; Fauci and Peskin, 1988),cardiac blood flow (Kovacs et al., 2001; McQueen and Peskin, 1997; McQueen and Peskin, 2000; McQueen and Peskin, 2001) and ciliary driven flows (Grunbaum et al.,1998).

**u**is the fluid velocity,

*p*is the pressure,

**F**is the force per unit area applied to the fluid by the immersed wing, ρ is the density of the fluid, μ is the dynamic viscosity of the fluid,

*t*is dimensional time and

**x**is the position. The dimensionless variables are

**u**′,

**x**′,

*p*′ and

*t*′, which represent the non-dimensional velocity, position,pressure and time, respectively.

**F**′ is the non-dimensional force per unit volume,

*L*is a characteristic length (such as the chord length of a wing) and

*U*is a characteristic velocity (such as the average wing tip velocity).

*Re*is then given by Eqn 4 and may be thought of as being roughly proportional to the ratio of inertial to viscous forces in the fluid.

*m*′ is the dimensionless mass,

*T*is the tension,

*k*′ is the dimensionless bending stiffness,

**f**′ is the dimensionless force acting on the wing,

*n̂*is the unit vector pointing normal to the wing, and

*r̂*is the unit vector tangent to the wing. As above, time is scaled by the chord of the wing divided by the characteristic velocity. The dimensionless wing mass and bending stiffness may then be written as:

_{s}is the mass per unit length of the wing and

*k*

_{beam}is the dimensional flexural stiffness.

In this paper, the *Re* is set to 10, corresponding to the case of some of the smallest flying insects such as *Thripidae frankliniella*(Sunada et al., 2002). The non-dimensional bending stiffness is varied from about 0.25 to 4. The mass of a thrips is about 6.0 × 10^{–8} kg(Tanaka, 1995) and its wing length is about 0.75 mm. Assuming that a wing weighs about 3.0 ×10^{–9} kg (this is 1/20 of the total mass, and the feathery wings probably weigh much less) and the wing length is 0.75 mm,ρ _{s} would be roughly equal to 4.0 × 10^{–6}kg m^{–1}. Given a wing chord of about 0.25 mm and an air density of 1.2 kg m^{–3}, *m*′ would equal 0.016. Because this value is much less than one, the effects of inertia on the wing are ignored. The wing is modeled as a massless boundary that resists bending and stretching.

The basic idea behind the immersed boundary numerical method is as follows.

At each time step, calculate the forces that the boundaries impose on the fluid. These forces are determined by the deformation of the elastic boundaries. Additional external forces used to drive the motion of the boundary may also be applied to the fluid.

Spread the force from the Lagrangian grid describing the position of the boundaries to the Cartesian grid used to solve the Navier–Stokes equations.

Solve the Navier–Stokes equations for one time step.

Use the new velocity field to update the position of the boundary. The boundary is moved at the local fluid velocity, enforcing the no-slip condition.

Further details of the immersed boundary method and its discretization may be found in the Appendix.

### Numerical simulations

The two-dimensional numerical simulations of flight in this paper were constructed to be similar to the physical experiments of Dickinson and Götz (Dickinson and Götz,1993) and previous two-dimensional numerical simulations of clap and fling (Miller and Peskin,2005). Dickinson and Götz used an aluminium wing with a chord of 5 cm immersed in a sucrose solution with a kinematic viscosity of 0.0000235 m^{2} s^{–1} (about 20 times that of water) moving with a characteristic velocity in the range of 0.04–0.12 m s^{–1}. The dimensions of the sucrose tank used in the physical experiment were 1 m in length by 0.4 m in width. The same parameters as listed for this physical experiment were used in all of the following numerical experiments with two exceptions: (1) the size of the computational tank was increased to reduce wall effects at lower *Re* and (2) the translational velocity was changed to simulate *Re*=10. In the following simulations, we use a computational tank that is 1 m × 1 m in size. For numerical convenience, we place this tank within a slightly larger periodic domain, of size (1 m+30 *h*) × (1 m+30 *h*),where *h*=Δ*x*=Δ*y* is the mesh width of our fluid grid. The edges of the computational tank are made of immersed boundary points that are linked by stiff springs to stationary target points. The region within the four walls is called the `computational tank'. The Navier–Stokes equations were solved on a 1230×1230 Cartesian grid,and each wing was discretized on a Lagrangian array of 120 points. Miller and Peskin (Miller and Peskin,2005) presented a convergence analysis that showed that this mesh size is within the range of convergence for the two-wing problem at *Re* below 100.

Unless otherwise stated, the motion of the flexible wings was prescribed by attaching target points to the top 1/5 of the boundary along the leading edge of the wing with springs (Fig. 2). The target points moved with the prescribed motion and applied a force to the boundary proportional to the distance between the target and corresponding boundary points. The bottom 4/5 of the wing (trailing edge) was free to bend. This has the effect of modeling a flexible wing with a rigid leading edge. In the `nearly-rigid' case, springs were attached to target points along the entire length of the wing, which prevented any significant deformation. The dimensional stiffness coefficient of the springs that attach the boundary points to the target points was *k*_{targ}=1.44×10^{5} kg s^{–2},and the stiffness coefficient for the tension or compression of the wing was also set to *k*_{str}=1.44×10^{5} kg s^{–2}. These values were chosen to prevent any significant stretching or deformation of the wings in the rigid case.

The dimensional flexural or bending stiffness of the wing was varied from *k*_{beam}=0.125 κ to 2 κ Nm^{2}, whereκ=5.5459×10^{–6} Nm^{2}. As stated above,this range of bending stiffnesses corresponds to dimensionless values ranging from 0.25 to 4. We chose this value of κ to represent the case where deformations during translation are small but bending does occur when the wings are close. This range of values was found to produce deformations that are qualitatively similar to those observed in flight videos.

*C*

_{L}is the lift coefficient,

*C*

_{D}is the drag coefficient,

*S*is the surface area per unit length of the model wing,

*F*

_{D}is the drag force per unit length and

*F*

_{L}is the lift force per unit length. As

*C*

_{L}and

*C*

_{D}vary with time, they may also be thought of as normalized lift and drag forces. The normalization does not change with time, so the actual force plots would maintain the same shape. It should be noted that these definitions are designed for high

*Re*,and in this intermediate range lift and drag coefficients become functions of

*Re*.

For the smallest insects that probably transport themselves on gusts of wind, we use the ratio of the average lift force produced during the stroke to the average drag force as a simple measure of flight efficiency. Lift (the vertical component of the force) is generated to help keep the insects afloat while drag (the horizontal component of the force) is of less importance since any thrust produced would likely be swamped by wind gusts. This metric is particularly appropriate when the wings are close to each other because the horizontal component of the force acting on each wing cancels. It is worthwhile to note that a number of other measures of efficiency are used in the literature, many of which are likely to be more appropriate for medium to large insects (e.g. Wang,2004).

### Kinematics of the clap and fling strokes

Quantitative descriptions of the wing beat kinematics of the smallest flying insects are currently unavailable. This is partially due to the fact that these insects are extremely difficult to film. Many of these insects may be as small as 0.25 mm in length and flap their wings at frequencies of 200 Hz or greater (Dudley, 2000). The authors know of several unpublished high-speed videos of *Thysanoptera,Muscidifurax raptor* and *Nasonia vitripennis* shot from a single camera, but quantitative reconstruction of the wingbeat kinematics is not possible from one camera angle. The simplified flight kinematics of clap and fling used in this paper are similar to those used by Lighthill(Lighthill, 1973), Bennett(Bennett, 1977), Maxworthy(Maxworthy, 1979), Spedding and Maxworthy (Spedding and Maxworthy,1986), Sun and Xin (Sun and Xin, 2003), Miller and Peskin(Miller and Peskin, 2005), and Chang and Sohn (Chang and Sohn,2006) to investigate flight aerodynamics using mathematical and physical models.

Either a single clap upstroke or a single fling downstroke was simulated. This simplification was made because the influence of the wake produced by the previous stroke is small for *Re* of 10 and below. Chang and Sohn(Chang and Sohn, 2006) found changes in the strength of the leading edge vortices at *Re* on the order of 100 when isolated clap and fling motions were compared with cyclical clap and fling motions; however, the overall aerodynamics and forces acting on the wings at *Re*=10 were quite similar.

*V*is the maximum translational velocity during the stroke,

*v*(τ) is the translational velocity at dimensionless time τdefined by Eqn 12,

*t*is the actual time,

*c*is the chord length of the wing,τ

_{accel}is the dimensionless time when translational acceleration begins and Δτ

_{accel}is the dimensionless duration of translational acceleration. Δτ

_{accel}was set to 1.3. For the clap strokes, τ

_{accel}was set to 0. For the fling strokes, τ

_{accel}was set to 0 (translation starts at the beginning of wing rotation), 0.435, 0.87, 1.305 or 1.74 (translation starts at the end of rotation). After acceleration, the translational velocity of the wing was fixed as

*V.*

_{decel}is the dimensionless time when translational deceleration begins and Δτ

_{decel}is the dimensionless duration of translational deceleration. The translational velocity during the fling stroke is symmetric to the downstroke and may be constructed similarly. Unless otherwise noted, τ

_{final}was taken to be 6.87 (this gives a translation of about 3.8–5.5 chords depending upon the kinematics),where τ

_{final}is the duration of each half stroke.Δτ

_{decel}was taken to be 1.3, and

*V*was set to 4.7×10

^{–3}ms

^{–1}. For the clap strokes,τ

_{decel}was set to 3.83 (rotation starts with deceleration),4.265, 4.7, 5.135 or 5.57 (rotation starts at the end of translation).

*x*-axis (the origin is defined as the intersection of the wing with the

*x*-axis at the initial time). The angular velocity of the left wing during the rotational phase at the end of the upstroke is given by:

_{rot}is a constant determined by the total angle of rotation and by the duration of the rotational phase in Eqn 15, ω(τ) is the angular velocity as a function of dimensionless time, τ

_{turn}is the dimensionless time wing rotation begins, Δτ

_{rot}is the dimensionless duration of the rotational phase and Δθ is the total angle through which rotation occurs. Δθ was set to 45 deg., andΔτ

_{rot}was set to 1.74 in all simulations. Rotation at the beginning of the downstroke was constructed similarly.

## RESULTS

### Translational/rotational overlap

For this set of simulations, we vary the translational-rotational overlap during wing rotation for clap and fling strokes and consider nearly rigid and flexible wings (*k*_{beam}=κ). For the clap simulations,the wings come to rest at the beginning, first quarter, middle quarter, third quarter or end of rotation. This corresponds to 0, 0.25, 0.5, 0.75 and 1.0 overlap between rotation and translation (τ_{decel}=3.83, 4.265,4.7, 5.135 and 5.57). For the fling simulations, translation begins at the beginning, first quarter, middle quarter, third quarter or end of rotation. Similarly, this corresponds to 1.0, 0.75, 0.5, 0.25 and 0 overlap between rotation and translation (τ_{accel}=0, 0.435, 0.87, 1.305 and 1.74).

Lift coefficients and drag coefficients as functions of time (expressed as the fraction of the stroke) for the four cases (rigid clap, rigid fling,flexible clap and flexible fling) are shown in Figs 4 and 5. For the clap strokes, the first peak in the force coefficients corresponds to the acceleration of the wing from rest. The force coefficients quickly reach steady values until the wings begin to decelerate and rotate. The next peak in the force coefficients corresponds to the rotation of the wings together (clap). Large lift and drag forces are produced as the air between the wings is squeezed downward. Note that lift coefficients produced during the rigid and flexible clap strokes are comparable, but the maximum drag coefficients produced during the clap are significantly lower in the flexible case.

For the fling strokes, the first peak in the force coefficients corresponds to the rotation of the wings apart (fling). The peak lift forces in the flexible cases are higher than in the corresponding rigid cases. There are several possible explanations for this phenomenon. Peel delays the formation of the trailing edge vortices, thereby maintaining vortical asymmetry and augmenting lift for longer periods (Miller and Peskin, 2005). Similarly, the suppression of the formation of the trailing edge vortex would also reduce the Wagner effect at the beginning of translation (Weis-Fogh,1973). Lift augmentation by `peel' *versus* `fling' was also predicted using simple analytical models of inviscid flows around rigid wings (Ellington, 1984b). As the wings are peeled apart, the wings are deformed and store elastic energy. As the wings translate apart, the wings `straighten' and push down on the fluid, causing an upwards lift force. Not only is the lift greater but the drag is also lower in the flexible case; thus, further enhancing the advantage of a flexible wing.

For a comparison of the average forces generated in each case, lift and drag coefficients were averaged over wing rotation and 2.5 chord lengths of travel as shown in Fig. 6. Please note that this is equivalent to take the ratio of the average lift force to the average drag force produced. For the fling cases, average lift coefficients increase as the amount of rotational/translational overlap increases. In the case of 100% overlap, average lift forces are larger in the flexible case than in the rigid case. For other degrees of overlap, average lift coefficients for rigid and flexible wings are comparable. For all overlapping fling cases, the average drag coefficients produced by flexible wings are lower than those produced by rigid wings. Moreover, the average ratio of lift to drag produced during fling is higher for flexible wings than for rigid wings. This effect increases as the degree of rotational/translational overlap is increased. For the clap cases, average lift coefficients are comparable for rigid and flexible wings, and these forces increase as the degree of rotational/translational overlap increases. The average drag forces produced in the rigid wing cases, however, are higher than the corresponding flexible wing cases. As a result, average lift over drag is higher for flexible wings. In summary, these results suggest that wing flexibility could improve the efficiency of hovering flight at low *Re*.

In order to answer the age-old question, what is the sound of one wing clapping, we now consider a single wing moving according to the same kinematics as in the two-wing case. The force coefficients as functions of time are shown in Fig. 7. During the fling motion, the first peak in the force coefficients corresponds to the forces generated during rotation, and the second peak corresponds to the forces generated during translational acceleration. Wing deformation during one-winged clap and fling is minimal, and these force traces are similar to the rigid wing case (data not shown). By comparing the scale bar between Figs 5 and 7, one can easily see the that maximum drag forces generated during two-winged fling are 10 times higher than the maximum generated by one wing (either rigid or flexible). During the clap motion, the first peak in the force coefficients corresponds to the forces generated during translational acceleration. The first dip in the force coefficients towards the end of the stroke corresponds to the deceleration of the wing. As the wing rotates towards the end of the stroke, another peak in the lift and drag forces force is observed. These forces quickly drop as the wing decelerates to rest. Force coefficients for one-winged clap and fling averaged over rotation and 2.5 chord lengths of translation are shown in Fig. 8.

Streamline plots of the fluid flow around wings performing a fling with 100% overlap between the translational and rotational phases are shown in Fig. 9. The streamlines are curves that have the same direction as the instantaneous fluid velocity **u**(**x**,*t*) at each point. The curves were drawn by making a contour map of the stream function because the stream function is constant along streamlines. The density of the streamlines is proportional to the speed of the flow. Color has been added to the streamline plots to help the reader distinguish individual streamlines and vortices. Regions of negative vorticity appear as warm colors and positive vorticity appear as cool colors. Close up views of the deformation of the left wing for each case at four points in time are shown in Fig. 10.

Two-winged fling with rigid wings is shown in Fig. 9A. During wing rotation,two large leading edge vortices begin to form as the wings fling apart(i–ii). As translation begins, a pair of trailing edge vortices forms and begins to grow in strength (ii–iv). Two-winged fling with flexible wings (*k*_{beam}=κ) is shown in Fig. 9B. The wings move with the same motion as in Fig. 9A(100% rotational/translational overlap). As the wings move apart, the point of separation moves from the leading edge to the trailing edge of the wing(i–ii). The formation of the trailing edge vortices occurs later in the stroke, and the trailing edge vortices are relatively weaker than the rigid wing case (ii–iv). One flexible wing (*k*_{beam}=κ)moving with the same fling motion as A and B is shown in Fig. 9C. Because of the smaller aerodynamic forces acting on the wing, its deformation is negligible and is very close to the rigid wing case (Fig. 10).

Streamline plots of the flow around wings performing a clap with 100%overlap between the translational and rotational phases are shown in Fig. 11. Two-winged clap with rigid wings is shown in Fig. 11A. As the wings clap together, the fluid is pushed out between the trailing edges causing an upwards lift force (i–iv). Two-winged clap with flexible wings (*k*_{beam}=κ) using the same motion is shown in Fig. 11B. Towards the end of the stroke, the wings bend as they are clapped together, reducing the peak drag forces generated (ii–iv). In addition, the point of`attachment' moves from the leading edge to the trailing edge of the wing. One flexible wing (*k*_{beam}=κ) moving with the same clap motion is shown in Fig. 11C. Because of the smaller aerodynamic forces acting on the wing, its deformation is negligible.

### Varying wing flexibilities

In this set of simulations, the flexural stiffness of the wings was varied from 0.25 κ to 2 κ, and the translational/rotational during clap and fling overlap was set to 100%. The lift and drag coefficients as functions of time for clap and fling cases are shown in supplementary material Fig. S3. The largest lift forces are produced for the case where *k*_{beam}=1.25 κ. Lift coefficients decrease as the bending stiffness of the wing increases or decreases from this value. For flexible two-winged clap, the lift coefficients are comparable for all five values of the bending stiffness. During flexible fling, the drag coefficients increase with increasing bending stiffness. Drag coefficients generated during two-winged flexible clap also increase with increasing bending stiffness.

Average lift and drag during wing rotation and a 2.5 chord translation are shown in Fig. 12. Average lift coefficients for flexible clap and fling are shown in Fig. 12A. The average lift coefficient generated during fling was greatest when the bending stiffness was set to 1.25 κ. The average lift generated during clap was relatively constant for this range of values. The average drag coefficients for flexible clap and fling are shown in Fig. 12B. For the both clap and fling, the average drag coefficients increase with increasing bending stiffness. Average lift over drag ratios for flexible clap and fling are shown in Fig. 12C. Lift over drag increases with decreasing bending stiffness.

### Varying the rigid section of the wing

To investigate the effect of wing stiffness asymmetries on the forces produced during flight, the rigid section of the flexible wing (1/5 of the chord length) was moved from the leading to the trailing edge of the wing in five steps. In all cases, the flexural stiffness of the wings was set to 1.0κ. The rotational/translational overlap was set to 100%. Combes and Daniels measured the flexural stiffness of *M. sexta* wings as a function of distance along the chord and found that the bending stiffness decreases from the leading to the trailing edge of the wing(Combes and Daniels, 2003). A quick look at the wing morphology of most insect wings suggests that this is true for many species, making the assumption that flexural stiffness is proportional to wing thickness. There might be a few exceptions to this rule,however. For example, if the bristles of thrips' wings are more flexible than the solid portion, then there may be some variation in the location of the stiffest portion of the wing as a function of distance from the leading to the trailing edge.

Lift coefficients and drag coefficients as functions of time (expressed as the fraction of the stroke) for clap and fling with the various wing designs are shown in supplementary material Fig. S4. In all cases, the rotational/translational overlap was set to 100% and the flexural stiffness of the wing was set to *k*_{beam}=κ. The rigid portion of the wing was also the location where the external force used to move the wing was applied. During flexible fling, the largest lift forces were produced in the case where the middle of the wing was made rigid. The smallest lift forces where produced in the case where the trailing edge of the wing was rigid. During flexible clap, the lift forces produced for all wing designs are comparable. The largest drag forces during flexible fling were produced when the middle portion of the wing was rigid. These forces are significantly larger than all other cases. During clap, the smallest peak drag forces were produced when the leading edge of the wing was rigid.

Average lift and drag coefficients during rotation and a 2.5 chord translation are shown in Fig. 13. Average lift coefficients as a function of the wing design(location of the rigid portion of the wing) for flexible clap and fling are shown in Fig. 13A. The average lift coefficient generated during fling was greatest when the middle part of the wing was rigid. The lowest lift forces were generated when the trailing edge of the wing was rigid. The average lift generated during clap was relatively constant for each wing design. Average drag coefficients are shown in Fig. 13B. In the case of fling, the largest average drag coefficients were generated when the middle part of the wing was rigid. These forces dropped as the rigid part was moved to either the leading or trailing edge of the wing. In the case of clap, the forces generated were relatively constant. Average lift over drag for clap and fling as a function of the location of the rigid part of the wing are shown in Fig. 13C. Average lift over drag increases as the rigid part is moved towards the leading edge of the wing.

Streamline plots of the flow around wings during fling for three wing designs are shown in supplementary material Fig. S5, and plots of the wing configurations at four points in time are shown in Fig. 14. When the leading edge is rigid, two large leading edge vortices begin to form as the wings fling apart. As translation begins, a pair of trailing edge vortices forms and begins to grow in strength. When the rigid wing is in the middle of the wing,large deformations of the wings occur along the leading and trailing edges as the wings are pulled apart. When the leading and trailing edges finally separate, large leading edge vortices are formed. When the rigid section is located at the trailing edge, the point of separation moves from the trailing to the leading edge of the wing as the wings are pulled apart. Large trailing edge vortices are formed at the beginning of the stroke, and the leading edge vortices are small.

Streamline plots of the flow around wings performing clap with three different wing designs are shown in supplementary material Fig. S6. For the case of a rigid leading edge towards the end of the upstroke, the wings bend as they are clapped together, reducing the peak drag generated. When the rigid section is in the middle or trailing edge of the wing, the fluid dynamics and forces generated are similar to the rigid wing case since wing deformations are minimal.

## DISCUSSION

The results of this study suggest that wing flexibility may be important for reducing drag forces generated during clap and fling at low *Re*,and some lift augmenting effects may also be produced. This is significant since the drag forces generated during two-wing fling may be as much as 10 times higher than those generated for one wing moving with the same motion. For flexible wings, the fling part of the stroke is more like a peel, and the clap part of the stroke is more like a reverse peel. When wing translation begins with wing rotation (100% overlap), lift forces are higher for flexible wings. As the wings peel apart, the point of separation travels from the leading to the trailing edge of the wing. This delays the formation of the trailing edge vortices, which reduces the Wagner effect and sustains vortical asymmetry (large leading edge vortex and small trailing edge vortex) for a larger portion of the stroke (please see the discussion below). Increased lift production by wing peel *versus* fling was also predicted by Ellington(Ellington, 1984b).

Asymmetries in flexural stiffness along the wing chord also influence aerodynamic performance. Wings that are more rigid along the trailing edge of the wing maximize the average lift/drag forces produced during fling. Wings that are more rigid in the middle maximize the average lift force produced during fling. This result suggests that differences in wing design could reflect different performance parameters that an insect has `maximized'.

Average lift forces were greatest during fling when the flexural stiffness of the wing was set to *k*_{beam}=1.25 κ. Using this value for the flexural stiffness of the wings, deformations during translation or single wing flapping are minimal, and the force coefficients produced are comparable with the rigid wing case. These results agree with the qualitative observations of Ellington (Ellington,1980) who oscillated single thrips' wings at the frequency and amplitude characteristic of flight. As the forces generated during the clap and fling portion of the strokes are so much larger than the single wing case because of wing–wing interactions, wing deformations in the numerical simulations presented in this paper are significant and greatly influence the aerodynamics. Measurements of the actual flexural stiffness of tiny insect wings and experimental work that considers wing–wing interactions at low *Re* are needed to verify this effect.

### Relating the wake to the forces generated

*R*_{n}*1*and

*R*_{p}*1*) of equal strength and opposite sign were formed and remain attached to the wing. During translation, two small trailing edge vortices of equal strength and opposite sign begin to form and grow in strength (

*R*_{n}*2*and

*R*_{p}*2*). Let the rest of the fluid domain

**be of negligible vorticity. Note that the subscript**

*R*_{f}*n*denotes regions of negative (clockwise) vorticity and

*p*denotes regions of positive (counterclockwise) vorticity. The total lift acting on both wings can then be defined as follows(Miller and Peskin, 2005):

*Re*of the wing in (Fig. 15A) is larger than that of(Fig. 15B) but both wings translate at

*Re*<30 (this prevents the separation of the trailing edge vortex). During translation, a leading edge vortex(

**) and a trailing edge vortex(**

*R*_{n}**) of equal strength and opposite sign were formed and remain attached to the wing. Let the rest of the fluid domain**

*R*_{p}**be of negligible vorticity. The vortices formed are more diffuse for the lower**

*R*_{f}*Re*. The total drag force acting on both wings can then be defined as follows:

*y*-direction increases the relative drag forces generated. In other words, larger vortices that move away from the wing in the vertical direction generate larger drag forces than compact vortices.

### Limitations of the model

Although the lift forces generated for flexible clap and fling using the particular wing kinematics described in this paper are in some cases larger than the corresponding lift forces generated for the rigid wing, one cannot extend these results to flexible wings in general. In fact, a flexible single wing using the same motion generates lower lift than a rigid wing. Because the wing kinematics studied in this paper do not represent the optimal case, it is not necessarily true that the optimal flexible wing stroke would outperform the optimal rigid wing stroke. What can be concluded from this work is that(1) the drag forces generated from wing–wing interactions can be an order of magnitude larger than a single wing, and (2) the addition of flexibility can reduce the drag but the maximum and average drag forces are still substantially larger than the single wing case.

Ideally, one would like to be able to explore a wide parameter space of wing beat kinematics and find the optimal rigid and flexible wing strokes,similar to the studies presented by Berman and Wang(Berman and Wang, 2007) for rigid wings and Alben (Alben,2008) for flexible appendages. Unfortunately, one cannot make the quasi-steady or inviscid assumptions at *Re*=10 to make this sort of analysis feasible. It may be possible to determine optimal rigid and flexible wing strokes with the use of physical models, but new custom experimental systems will need to be designed to measure the small forces generated by flapping appendages at this scale.

### Implications for bristled wings

Although wing flexibility reduces the amount of force needed to clap the wings together and fling them apart, these forces are still significantly larger than the forces generated during single-wing translation. In addition,lift over drag ratios are lower for two-winged clap and fling than for one-winged translation. It could be the case that tiny insects sacrifice aerodynamic efficiency for increased lift. Another possibility for some tiny insects is that wing fringing further reduces the force required to clap the wings together and fling the wings apart. During wing rotation in the clap and fling, there could be some flow between the wings' bristles, which would reduce the aerodynamic forces generated. If the spacing of the bristles and the *Re* is near the transition where the bristled appendages act either as leaky rakes or solid paddles(Cheer and Koehl, 1987), then it could be possible for the wings to act as solid plates during translation. This might allow for lift to be preserved during flight.

## APPENDIX

**u**(

**x**,

*t*) is the fluid velocity,

*p*(

**x**,

*t*) is the pressure,

**F**(

**x**,

*t*) is the force per unit area applied to the fluid by the immersed wing, ρ is the density of the fluid and μ is the dynamic viscosity of the fluid. The independent variables are the time

*t*and the position

**x**. Note that bold letters represent vector quantities. Eqns A1 and A2 are the Navier–Stokes equations for viscous flow in Eulerian form. Eqn A2 is the condition that the fluid is incompressible.

**f**(

*r, t*) is the force per unit length applied by the wing to the fluid as a function of Lagrangian position and time, δ(

**x**)is a two-dimensional delta function,

**X**(

*r, t*) gives the Cartesian coordinates at time

*t*of the material point labeled by the Lagrangian parameter

*r*. Eqn A3 applies force from the boundary to the fluid grid, and Eqn A4 evaluates the local fluid velocity at the boundary. The boundary is then moved at the local fluid velocity and this enforces the no-slip condition. Each of these equations involves a two-dimensional Dirac delta function δ, which acts in each case as the kernel of an integral transformation. These equations convert Lagrangian variables to Eulerian variables and

*vice versa*.

These equations describe the forces applied to the fluid by the boundary in Lagrangian coordinates. Eqn A5describes the force applied to the fluid as a result of the difference between the actual position of the wing and the position of a target boundary; the motion of which serves to drive the motion of the wing in this work. The function **f**_{targ}(*r, t*) gives the external force per unit length applied to the wing, *k*_{targ} is a stiffness coefficient and **Y**(*r, t*) gives the desired position of the wing as a function of time. Eqn A6describes the force applied to the fluid as a result of the deformation of the actual boundary, which here is modeled as a beam. The function **f**_{beam}(*r, t*) gives the force per unit area of wing that results from the bending of the beam and *k*_{beam} is the corresponding stiffness coefficient. Eqn A7 describes the force applied to the fluid per unit area of the wing as a result of the resistance to stretching by boundary given as **f**_{str}(*r, t*), where *k*_{str} is the corresponding stiffness coefficient. Finally, Eqn A8 describes the total force applied to the fluid per unit length, **f**(*r, t*), as a result of both the target boundary and the deformation of the boundary.

*t*be the duration of the time step, let Δ

*r*be the length of a boundary spatial step, let

*n*be the time step index and let

*r*define the location of a boundary point in the Lagrangian framework(

*r*=

*m*Δ

*r, m*is an integer). Also, let

**X**

*(*

^{n}*r*)=

**X**(

*r,n*Δ

*t*),

**Y**

*(*

^{n}*r*)=

**Y**(

*r,n*Δ

*t*),

**u**

*(*

^{n}**x**)=

**u**(

**x**,

*n*Δ

*t*) and

*p*(

^{n}**x**)=

*p*(

**x**,

*n*Δ

*t*). For any function ϕ(

*r*), let(D

_{r}ϕ)(

*r*) be defined by the following equation:

**f**

_{targ}, defined by Eqn A5, is discretized as follows:

**f**

_{beam}, defined by Eqn A8, is discretized as follows:

*n*is the total number of grid points in the boundary andδ

_{f}*is the Kronecker symbol, which is defined as follows:*

_{ml}*to collapse the summations to individual terms but it is actually advantageous not to do this. The reason is that Eqn A11 automatically handles the otherwise special cases*

_{ml}*l*=1, 2,

*n*2 and

_{f}–*n*1. An efficient implementation of Eqn A11 which handles these cases seamlessly is to loop over

_{f}–*m*rather than

*l*and, for each

*m*, to compute the contribution to all of the relevant values of(

*f*

_{beam})

*. For more details on the discretization of the bending force, see Zhu and Peskin(Zhu and Peskin, 2002). The stretch force,*

_{l}**f**

_{str}, defined by Eqn A7, is discretized as follows:

*, used in these calculations, is given by the following equations:*

_{h}*h*is the size of the spatial step in the Cartesian grid and

**x**=(

*x, y*).

*, Eqns A3, A4 can be discretized as follows:*

_{h}*r*=

*m*Δ

*r*, and

**x**=(

*ih,jh*), where

*i, j, m*are integers and

*h*is the mesh width.

**D**

^{0}is the central difference approximation to ▿. It is defined as follows:

_{1}, e

_{2}} is the standard basis of

**R**

^{2}. The operators

*x*

_{α}. They are defined as follows:

*S*(

_{h}**u**) is a skew-symmetrical difference operator,which serves as a difference approximation of the non-linear term

**u•▿u**. This skew-symmetrical difference operator is defined as follows:

For further discussion, see Peskin(Peskin, 2002). Since the equations are linear in the unknowns **u**^{n}^{+1} and p^{n}^{+1}, the Fast Fourier Transform (FFT)algorithm was used to solve for **u**^{n}^{+1} and p^{n}^{+1} from **u*** ^{n}*,p

*and*

^{n}**F**

*.*

^{n}**LIST OF ABBREVIATIONS**

- c
chord length of the wing

*C*_{D}drag coefficient

*C*_{L}lift coefficient

**D**^{0}central difference approximation to ▿

- \(D_{{\alpha}}^{{\pm}}\)
forward and backward difference approximations of∂/∂

*x*_{α} - f′
dimensionless force acting on the wing

- F̄
mean lift force

- F
force per unit area applied to the fluid by the immersed wing

- F′
non-dimensional force per unit volume

*F*_{D}drag force per unit length

*F*_{L}lift force per unit length

**F**(**x**,*t*)force per unit area applied to the fluid by the immersed wing

**f**_{beam}(*r*,*t*)force per unit area of wing that results from the bending of the beam and

*k*_{beam}is the corresponding stiffness coefficient**f**(*r*,*t*)force per unit length applied by the wing to the fluid as a function of Lagrangian position and time

*t***f**_{str}(*r*,*t*)force applied to the fluid per unit area of the wing as a result of the resistance to stretching by boundary

**f**_{targ}(*r*,*t*)external force per unit length applied to the wing

- h=Δx=Δy
mesh width of fluid grid

- k′
dimensionless bending stiffness

*k*_{beam}dimensional flexural stiffness

*k*_{str}stiffness coefficient for the tension or compression of the wing

*k*_{targ}dimensional stiffness coefficient of the target springs

- L
characteristic length

- m′
dimensionless mass

- n̂
unit vector pointing normal to the wing

- n
as subscript is negative (clockwise) vorticity

- n
time step index

- n
_{f}total number of grid points in the boundary

- p
pressure

- p
as subscript is positive (counterclockwise) vorticity

- p′
non-dimensional pressure

*p*(**x**,*t*)pressure

- R
_{n}leading edge vortex

- R
_{p}trailing edge vortex

- R
_{f}fluid domain

- r
location of a boundary point in the Lagrangian framework

- r̂
unit vector tangent to the wing

- S
surface area per unit length of the model wing

*S*(_{h}**u**)skew-symmetrical difference operator

- t
dimensional time

- t′
non-dimensional time

- T
tension

- U
velocity vector in Lagrangian framework

- U
characteristic velocity

- u
fluid velocity

- u
represents the velocity of the separation point

- u′
non-dimensional velocity

- u
_{c}average velocity of air mass

- unzipping
velocity of air mass between wings

- u
_{z}(t)unzipping velocity of air mass between wings

**u**(**x**,*t*)fluid velocity in Eulerian framework

- V
maximum translational velocity during the stroke

**X**(*r*,*t*)Cartesian coordinates at time

*t*of the material point labeled by the Lagrangian parameter*r*- x
position in Cartesian coordinates

- x
exposed portion of the wing

- x′
non-dimensional position

- x
_{e}length of the exposed section of the wing

**Y**(*r*,*t*)desired Eulerian position of the wing

- β
constant half-angle between the wings

- Γ
circulation

- Δr
length of a boundary spatial step

- Δt
duration of the time step

- Δθ
total angle through which rotation occurs.

- Δτaccel
dimensionless duration of translational acceleration

- Δτdecel
dimensionless duration of translational deceleration

- Δτrot
dimensionless duration of the rotational phase

- δ(x)
two-dimensional delta function

- δ
(_{h}**x**)discretization of the two-dimensional delta function

- δml
Kronecker symbol

- \({\hat{{\kappa}}}\)
curvature

- v(τ)
translational velocity at dimensionless time τ

- μ
dynamic viscosity of the fluid

- ρ
density of the fluid

- ρs
mass per unit length of the wing

- \(\begin{array}{c}\\{\Sigma}\\r\end{array}\)
sum over all of the discrete collection of points of the form

*r*=*m*Δ*r* - \(\begin{array}{c}\\{\Sigma}\\x\end{array}\)
sum over all of the discrete points of the form

*x*=(*ih*,*jh*), where*i*,*j*,*m*are integers - τaccel
dimensionless time when translational acceleration begins

- τdecel
dimensionless time when translational deceleration begins

- τfinal
dimensionless duration of each half stroke

- τturn
dimensionless time wing rotation begins

- ωrot
constant determined by the total angle of rotation and by the duration of the rotational phase in Eqn 15

- ω(τ)
angular velocity as a function of dimensionless time

- |ω|
absolute value of the vorticity

## FOOTNOTES

The authors would like to thank Ty Hedrick for the snap shots of *Muscidifurax raptor* and *Nasonia vitripennis* used in this paper. The work was funded by Miller's Burroughs Wellcome Fund Career Award at the Scientific Interface.

## References

**Alben, S.**(

**Alben, S.**(

**Alben, S., Shelley, M. and Zhang, J.**(

**Alben, S., Shelley, M. and Zhang, J.**(

**Bennett, L.**(

**Berman, G. J. and Wang, Z. J.**(

**Birch, J. M., Dickson, W. B. and Dickinson, M. H.**(

**Chang, J. W. and Sohn, M. H.**(

**Cheer, A. Y. L. and Koehl, M. A. R.**(

**Combes, S. A. and Daniels, T. L.**(

*Manduca sexta.*

**Denny, M.**(

**Dickinson, M. H. and Götz, K. G.**(

**Dudley, R.**(

*.*Princeton: Princeton University Press.

**Ellington, C. P.**(

*Thysanoptera*).

**Ellington, C. P.**(

**Ellington, C. P.**(

**Ellington, C. P.**(

**Etnier, S. A. and Vogel, S.**(

*Narcissus: Amaryllidaceae*) flowers in wind: drag reduction and torsional flexibility.

**Fauci, L. J.**(

**Fauci, L. J. and Fogelson, A. L.**(

**Fauci, L. J. and Peskin, C. S.**(

**Grunbaum, D., Eyre, D. and Fogelson, A.**(

**Heathcote, S., Wang, Z. and Gursul, I.**(

**Ishihara, D., Horie, T. and Denda, M.**(

**Kim, D. K., Kim, H. I., Han, J. H. and Kwon, K. J.**(

**Koehl, M. A. R.**(

**Lauder, G. V., Madden, P. G. A., Mittal, R., Dong, H. and Bozkurttas, M.**(

**Lehmann, F. O. and Pick, S.**(

**Lehmann, F. O., Sane, S. P. and Dickinson, M.**(

**Lighthill, M. J.**(

**Marden, J. H.**(

**Maxworthy, T.**(

**McQueen, M. C. and Peskin, C. S.**(

**McQueen, M. C. and Peskin, C. S.**(

**McQueen, D. M. and Peskin, C. S.**(

**Miller, L. A. and Peskin, C. S.**(

**Mittal, R.**(

**Norberg, R. A.**(

*Alucita pentadactyla*L. and

*Orneodes hexadactyla*L. (Microlepidoptera).

**Peskin, C. S.**(

*,*

**Sane, S. P.**(

**Spedding, G. R. and Maxworthy, T.**(

**Sun, M. and Xin, Y.**(

**Sunada, S., Kawachi, K., Watanabe, I. and Azuma, A.**(

**Sunada, S., Takashima, H., Hattori, T., Yasuda, K. and Kawachi,K.**(

**Tanaka, S.**(

**Vanella, M., Fitzgerald, T., Preidikman, S., Balaras, E. and Balachandran, B.**(

**Vogel, S.**(

*Drosophila*. II. Variations in stroke parameters and wing contour.

**Vogel, S.**(

**Wang, Z. J.**(

**Wang, Z. J.**(

**Wang, Z. J.**(

**Weis-Fogh, T.**(

**Weis-Fogh, T.**(

**Wu, J. C.**(

**Zhu, L. and Peskin, C. S.**(