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; WeisFogh, 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'. WeisFogh 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 (WeisFogh, 1973). The clap and fling motion has also been reported in the greenhouse whitefly Trialeurodes vaporariorum (WeisFogh, 1975) and Thrips physapus (Ellington, 1984a). Miller and Hedrick have observed clap and fling kinematics in highspeed videos of the tiny parasitoid wasps Muscidifurax raptor and the jewel wasp Nasonia vitripennis (see supplementary material Fig. 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; WeisFogh, 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 twodimensional, twolink 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).
While flexibility will probably reduce the drag forces required for clap and fling, the wings should not be so flexible that the lift forces produced are significantly diminished. Ellington suggested that a peel mechanism with flexible wings might actually serve to augment lift forces relative to the rigidfling case (Ellington, 1984b). In the case of peel, the wings are pulled apart along the leading edges and curve along the wing chords. As the peel progresses, the separation point between the wings moves from the leading edge to the trailing edge (see Fig. 1A,B). The peeling motion of the wings has been observed in images of 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: (1) where x_{e} is the length of the exposed section of the wing, and β is the constant halfangle between the wings and 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).
Similarly, flexible clap might be thought of as a reversed peel. As the wings clap together, the angle of attack increases and the wings meet at the leading edges. The meeting point between the wings then moves from the leading edge to the trailing edge. As described by Ellington (Ellington, 1984b) and visualized by Maxworthy (Maxworthy, 1979), the clap produces a downward jet of air and the wings consequently experience an upward force (see Fig. 1C,D). This downward jet and brief augmentation of the lift force was also observed in twodimensional numerical simulations (Miller and Peskin, 2005; Sun and Xin, 2003) and physical models of clap with rigid wings (Lehmann et al., 2005). The mean lift force during clap can be roughly estimated using a `reverse peel' model. The mass of the air per unit span between the wings isρβ c^{2}. This air mass gains velocity u_{c} in the duration 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: (2) This is an inviscid approximation of the mean force, and the force generated in the viscous case is likely to be smaller because some momentum is lost to viscosity.
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 twodimensional 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).
The dimensional analysis presented here is used to determine the forces relevant to the problem of clap and fling in small insects. The Navier–Stokes equations in three dimensions can be rewritten in dimensionless form as follows: (3) (4) (5) where 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 nondimensional velocity, position, pressure and time, respectively. F′ is the nondimensional 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.
Using the dimensional analysis described by Alben (Alben, 2008; Alben, 2009), the dimensionless equations describing a twodimensional wing immersed in a fluid can be written as follows: (6) where m′ is the dimensionless mass, T is the tension, k′ is the dimensionless bending stiffness, f′ is the dimensionless force acting on the wing, is the curvature, 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: (7) (8) where ρ_{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 nondimensional 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 noslip condition.
Further details of the immersed boundary method and its discretization may be found in the Appendix.
Numerical simulations
The twodimensional 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 twodimensional 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 twowing 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 `nearlyrigid' 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.
Lift and drag forces were calculated as functions of time by summing the forces that each immersed boundary point of the model wing applied to the fluid at each time step and taking the opposite sign of that value. As done in experiments, lift and drag coefficients were filtered to remove high frequency `noise' from the vibrations of the elastic boundary. This did not change the basic shape of the graphs. The lift and drag coefficients are defined as follows: (9) (10) where 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 highspeed 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.
At the end of the upstroke and the beginning of the downstroke, the wings were placed 1/10 chord lengths apart. This distance was chosen based on the analysis of Sun and Xin who showed that the forces acting on the wings during clap and fling do not change significantly for distances of 1/10 of the wing chord and lower (Sun and Xin, 2003). The kinematics of the left wing during the clap strokes (upstrokes) are described here. The right wing (when present) was the mirror image of the left wing at all times during its motion, and the kinematics of the fling (downstroke) was symmetrical to the upstroke. The translational velocities during the clap stroke were constructed using a series of equations to describe each part of the stroke. Plots of translational and angular velocities as functions of time are shown in Fig. 3. The velocity during acceleration at the beginning of the clap stroke is given by: (11) (12) where 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.
The translational velocities during deceleration of the left wing during the clap stroke are given as: (13) where τ_{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).
The angles of attack were similarly defined using a set of equations describing the angular velocity during the rotational phase of the stroke. Letα be defined as the angle of attack of the wing relative to the horizontal plane. For all clap strokes, the wings were translated at a constant angle of attack of 45 deg. during the upstroke and rotated to 90 deg. at the end of the upstroke. For all fling strokes, the wings were rotated fromα =90 deg. to α=45 deg. at the beginning of the downstroke. After rotation, the angle of attack was held constant for the remainder of the stroke. Let θ be defined as the angle between the left wing and the positive xaxis (the origin is defined as the intersection of the wing with the xaxis at the initial time). The angular velocity of the left wing during the rotational phase at the end of the upstroke is given by: (14) (15) where ω_{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 translationalrotational 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 (WeisFogh, 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 ageold question, what is the sound of one wing clapping, we now consider a single wing moving according to the same kinematics as in the twowing 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 onewinged 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 twowinged 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 onewinged 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.
Twowinged 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). Twowinged 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. Twowinged 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). Twowinged 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 twowinged 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 twowinged 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 twowing 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
To understand the aerodynamic mechanisms of lift and drag generation, we apply the general aerodynamic theory for viscous flows presented by Wu (Wu, 1981). For the particular case of fling shown in Fig. 15A, assume two wings were rotated apart along their trailing edges and are now translating away from each other along a horizontal plane. During rotation, two large leading edge vortices (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 R_{f} be of negligible vorticity. Note that the subscript 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): (16) where ω is the absolute value of the vorticity. This equation basically states that the total lift on both wings is proportional to the difference between the magnitude of the time rate of change of the first moment of vorticity associated with the leading edge vorticity and the time rate of change of the first moment of trailing edge vorticity. Therefore, suppression of the formation of the trailing edge vortex from the motion of the peel should transiently enhance lift forces during translation.
To understand the aerodynamic mechanism of drag generation, consider the case of a wing in translation as shown in Fig. 15B,C. Both wings are started from rest and are now translating at a constant velocity. The 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 (R_{n}) and a trailing edge vortex (R_{p}) of equal strength and opposite sign were formed and remain attached to the wing. Let the rest of the fluid domain R_{f} be of negligible vorticity. The vortices formed are more diffuse for the lower Re. The total drag force acting on both wings can then be defined as follows: (17) This equation basically states that increasing the time rate of change of vorticity in the ydirection 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 quasisteady 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 singlewing translation. In addition, lift over drag ratios are lower for twowinged clap and fling than for onewinged 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
The dimensional equations of motion for a twodimensional fluid are as follows: (A1) (A2) where 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.
The dimensional interaction equations between the fluid and the boundary are given by: (A3) (A4) where 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 twodimensional 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 noslip condition. Each of these equations involves a twodimensional 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.
The immersed boundary equations are given by: (A5) (A6) (A7) (A8)
These equations describe the forces applied to the fluid by the boundary in Lagrangian coordinates. Eqn A5 describes 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 A6 describes 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.
The system of differentiointegral equations given by Eqns A1, A2, A3, A4, A5, A6, A7, A8 was solved on a square grid with periodic boundary conditions. The Navier–Stokes equations were discretized on a fixed Cartesian grid, and the wings were discretized on a moving Lagrangian array of points. Consider the discretization of the equations that describe the force applied to the fluid as a result of the deformation of the boundary (Eqns A5, A6, A7, A8). Let Δ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: (A9)
Using these definitions, the target force, f_{targ}, defined by Eqn A5, is discretized as follows: (A10)
The bending force, f_{beam}, defined by Eqn A8, is discretized as follows: (A11) where , , n_{f} is the total number of grid points in the boundary andδ _{ml} is the Kronecker symbol, which is defined as follows: (A12)
Of course, Eqn A11 can be `simplified' by making use of the definition of δ_{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 l=1, 2, n_{f}–2 and n_{f}–1. An efficient implementation of Eqn A11 which handles these cases seamlessly is to loop over m rather than l and, for each m, to compute the contribution to all of the relevant values of (f_{beam})_{l}. For more details on the discretization of the bending force, see Zhu and Peskin (Zhu and Peskin, 2002). The stretch force, f_{str}, defined by Eqn A7, is discretized as follows: (A13)
Finally, the total force applied to the fluid by the boundary is given by the following equation: (A14)
Because the Lagrangian array of boundary points does not necessarily coincide with the fixed Cartesian grid used for the computation of the fluid velocities, a smoothed approximation to the Dirac delta function (Eqns A3, A4) is used to handle the fluid–boundary interaction. This approximate delta function applies a force from the boundary to the fluid and interpolates the fluid velocity. The approximate twodimensional Dirac delta function,δ _{h}, used in these calculations, is given by the following equations: (A15) where h is the size of the spatial step in the Cartesian grid and x=(x, y).
Using this choice of δ_{h}, Eqns A3, A4 can be discretized as follows: (A16) (A17) where form is the sum over all of the discrete collection of points of the r=mΔr, and is the sum over all of the discrete points of the form x=(ih, jh), where i, j, m are integers and h is the mesh width.
The final step is to discretize the incompressible Navier–Stokes equations. The method used in this paper is implicitly defined as follows: (A18) (A19) where D^{0} is the central difference approximation to ▿. It is defined as follows: (A20) where {e_{1}, e_{2}} is the standard basis of R^{2}. The operators are forward and backward difference approximations of ∂/∂x_{α}. They are defined as follows: (A21) (A22)
Thus, the viscous term given as is a difference approximation to the Laplace operator. Finally, S_{h}(u) is a skewsymmetrical difference operator, which serves as a difference approximation of the nonlinear term u•▿u. This skewsymmetrical difference operator is defined as follows: (A23)
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^{n} and 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 ▿
 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′
 nondimensional 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′
 nondimensional 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)
 skewsymmetrical difference operator
 t
 dimensional time
 t′
 nondimensional time
 T
 tension
 U
 velocity vector in Lagrangian framework
 U
 characteristic velocity
 u
 fluid velocity
 u
 represents the velocity of the separation point
 u′
 nondimensional 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′
 nondimensional position
 x_{e}
 length of the exposed section of the wing
 Y(r, t)
 desired Eulerian position of the wing
 β
 constant halfangle 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)
 twodimensional delta function
 δ_{h}(x)
 discretization of the twodimensional delta function
 δml
 Kronecker symbol
 curvature
 v(τ)
 translational velocity at dimensionless time τ
 μ
 dynamic viscosity of the fluid
 ρ
 density of the fluid
 ρs
 mass per unit length of the wing
 sum over all of the discrete collection of points of the form r=mΔr
 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
Supplementary material available online at http://jeb.biologists.org/cgi/content/full/212/19/3076/DC1
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.