## ABSTRACT

A computational fluid dynamic (CFD) modelling approach is used to study the unsteady aerodynamics of the flapping wing of a hovering hawkmoth. We use the geometry of a *Manduca sexta*-based robotic wing to define the shape of a three-dimensional ‘virtual’ wing model and ‘hover’ this wing, mimicking accurately the three-dimensional movements of the wing of a hovering hawkmoth. Our CFD analysis has established an overall understanding of the viscous and unsteady flow around the flapping wing and of the time course of instantaneous force production, which reveals that hovering flight is dominated by the unsteady aerodynamics of both the instantaneous dynamics and also the past history of the wing.

A coherent leading-edge vortex with axial flow was detected during translational motions of both the up- and downstrokes. The attached leading-edge vortex causes a negative pressure region and, hence, is responsible for enhancing lift production. The axial flow, which is derived from the spanwise pressure gradient, stabilises the vortex and gives it a characteristic spiral conical shape.

The leading-edge vortex created during previous translational motion remains attached during the rotational motions of pronation and supination. This vortex, however, is substantially deformed due to coupling between the translational and rotational motions, develops into a complex structure, and is eventually shed before the subsequent translational motion.

Estimation of the forces during one complete flapping cycle shows that lift is produced mainly during the downstroke and the latter half of the upstroke, with little force generated during pronation and supination. The stroke plane angle that satisfies the horizontal force balance of hovering is 23.6°, which shows excellent agreement with observed angles of approximately 20–25°. The time-averaged vertical force is 40% greater than that needed to support the weight of the hawkmoth.

## INTRODUCTION

Flapping flight is employed by most flying insects (approximately 750 000 species), birds and bats (nearly 8000 species). Animal flight covers a wide range of Reynolds numbers, from approximately 10 for the smallest insects with fringed wings to approximately 10^{4} for large insects such as certain beetles, moths and butterflies, and up to approximately 10^{5} for most birds and bats.

The question of how flapping animals fly can be translated from the aerodynamic point of view to the basic problem of the interaction between the wings and the surrounding air. The way the wings are put in motion, the characteristics of this motion and the main features of the wing–airflow interaction form the essence of the flight mechanism. Flying animals, in general, differ from conventional aircraft in beating their wings rapidly to hover and dart, in performing manoeuvres with rapid accelerations and decelerations, and in quick plunging and pitching motions. Therefore, we need to understand the aerodynamics of the flapping wings undergoing such highly three-dimensional and unsteady flapping motions, i.e. of the complicated three-dimensional, unsteady and viscous flow around wings with a variable geometry.

The goal of understanding the production of aerodynamic forces can be approached from either of two directions: by characterising the structure of the resultant wake and/or the near-field flow or by studying the time course of force generation by the flapping wings. The approach that visualises the resultant wake, combined with energy conservation assumptions, allows the force production of the wings to be calculated from the flows past them. The key point here is how accurately to characterise and model the unsteady wake pattern behind the flapping wings, which affects the advance ratio, the strength and the location of the wake and, hence, the forces on the body. Rayner (1979) proposed a vortex theory of animal flight based on the wake being modelled as a series of vortex rings periodically shed into the wake. However, flow visualisation of the wake behind a jackdaw *Corvus monedula* (Spedding, 1986) shows that this model is not quite realistic; the momentum of the vortex rings produced less than 35% of the momentum necessary to support the weight of the bird. Maxworthy (1979) succeeded in his pioneering experimental study in visualising the near-field flow around a rigid three-dimensional model of the ‘fling’ motion (Weis-Fogh, 1972) and found a leading-edge vortex with spanwise flow. Using smoke-visualisation of *Manduca sexta* (Willmott and Ellington, 1997) and of a hawkmoth-based robotic wing (Van den Berg and Ellington, 1997), this leading-edge vortex and the axial (spanwise) flow at the core of the vortex were clearly observed and confirmed to be essential for enhancing lift in insect flight.

The second approach, obtaining the time course of force generation by flapping wings, has also been followed both experimentally and numerically. There are relatively few measurements of the instantaneous forces produced by the three-dimensional flapping wings of animals. Jensen (1956) measured the mean lift of isolated locust wings. Cloupeau (1979) measured instantaneous lift in tethered locusts, and concluded that unsteady mechanisms must generate extra lift forces. Wilkin (1990) also succeeded in measuring the fluctuating lift forces on tethered locusts flying in a wind tunnel, finding that the peak lift was approximately twice the quasi-steady value. Wilkin (1993) and Wilkin and Williams (1993) produced further estimates of the flight forces of both locusts and moths, determining aerodynamic forces by subtracting estimates of the inertial forces from the measured forces. Dickinson (1996) used flow visualisation and instantaneous force measurements of tethered fruit flies (*Drosophila melanogaster*) to study the dynamics of force generation during flight, but could not obtain a cohesive picture of their relationship mainly because of the time delay involved in the force measurements.

Numerical studies of the unsteady flow around three-dimensional flapping wings are even fewer, largely because of the difficulties in modelling the complicated geometry and in simulating the three-dimensional movements of wings. One popular method that is available for computing the unsteady flow about arbitrary bodies is the potential flow panel method; a review of the various panel methods is given by Katz (1991). More recently, an extension of such a three-dimensional panel methods was developed by Vest (1996), in which the motion of the body is prescribed and each wing is assumed to be rigid and to rotate about a common axis. By using a time-stepping solution procedure, the shape of the trailing wake is determined; all vorticity is assumed to be confined to the boundary layer on the surface of the body and in a thin vortex wake shed behind the body. Smith (1996) also developed an unsteady aerodynamic three-dimensional panel method which can describe the kinematics of a flapping wing relatively realistically and can also account for wing flexibility using a finite element method in analysing the dynamic structure.

These panel methods are generally coupled with the discrete vortex method, and thus can predict the wake as well as the unsteady force generation. However, by ignoring the viscosity of fluid, they fail to predict reasonably the near-field flow, e.g. a large attached leading-edge vortex that has been recognised recently (Ellington *et al*. 1996) to play a key role in lift generation for insect flight at Reynolds numbers (*Re*) around 10^{3}. At such low *Re* values, flows around flapping wings behave unsteadily, usually demonstrating occurrence/break-down of vortices, separation and reattachment quite different from the behaviour of a fluid at high Reynolds numbers.

As a long-standing goal of studying the aerodynamics of flapping flight, the best approach is to solve the complete viscous and unsteady near-field flows around the flapping wings of animals. Solutions of the full Navier–Stokes equations for the three-dimensional, unsteady flow fields are challenging, and for the flow around variable-geometry bodies this problem is complicated even further. To understand the unsteady aerodynamics induced by a flapping wing, we have developed a computational fluid dynamic (CFD) model (Liu *et al*. 1995, 1996; Liu, 1997) designed to study both forward flight and hovering of a three-dimensional virtual wing which can have arbitrary variable geometry and can undergo realistic three-dimensional motions, accurately mimicking the movements of insect/bird wings. This CFD model is based on a robust and effective solver of the three-dimensional, incompressible, unsteady Navier–Stokes equations (Liu *et al*. 1995, 1996), which incorporates recent techniques. In the present study, by defining the three-dimensional geometry of a *Manduca*-based robotic wing, and by hovering the object with the three-dimensional movements of the wing of a hovering hawkmoth, we address an overall understanding of viscous and unsteady flows around a three-dimensional flapping wing during the whole process of translational and rotational motions, as well as the corresponding time course of instantaneous force generation.

## MATERIALS AND METHODS

### Resolution of the Navier–Stokes equations

The three-dimensional incompressible unsteady Navier– Stokes equations in the conservative form of momentum and mass are nondimensionalized in a local wingbase-fixed coordinate system (see Fig. 1A), written in an integral form as:

**q**=[

*u, v, w, p*]

^{T}in which

*u, v*and

*w*are three components of velocity, and

*p*is pressure. The term

**F**

_{i}expresses flow flux coming entering and leaving the cell boundary surfaces. The term

*V*(

*t*) denotes the volume of the cell (

*i*,

*j*,

*k*), a hexahedron constructed from eight grid points, and

*S*(

*t*) are its six boundary surfaces with unit outward normal vectors of

**n**. The influence of the moving grid is contributed by the velocity of the moving cell

**u**

_{g}. The term τ denotes the computational time. The governing equations are built and solved in the local wingbase-fixed coordinate system illustrated in Fig. 1A. A moving grid system is introduced to make grids fit the moving and deforming wing at each physical time step, and therefore to ensure sufficient grid density in resolving the viscous and unsteady flows around the wing surface.

The governing equations are discretized using the finite volume method (FVM) and are solved in a time-marching manner using the pseudo-compressibility technique. The time-dependent solution to the governing equations is realised by introducing an inner iteration to allow the divergence of velocity to vanish at each physical time step, and the time derivatives in the momentum equations are differenced using a first-order, two-point, backward-difference implicit formula (Liu *et al*. 1995).

In developing boundary conditions on the body surface, the no-slip condition is used for the velocity components: i.e. the velocity of a fluid particle is equal to the velocity of the wing at that point. This means that the wing has a smooth surface, although the roughness of the real hawkmoth wing may result in some discrepancy. For the dynamic effect due to rapid acceleration/deceleration of the wing, pressure divergence at the surface stencils is used for the pressure condition derived from the local momentum equation. At open boundary stencils, to satisfy the conservation of momentum, zero-gradient conditions are utilised for velocity components, which means that the local fluid can enter and leave the domain freely. Pressure is unchangeable at the open boundary. Details relating to the algorithm will be published (H. Liu, in preparation).

### Kinematics of the flapping wing

_{F}) and the angle of attack of the hindwing (α

_{H}). The positional angle, elevation angle and the angles of attack of the fore- and hindwing are defined (in rad) using the first three Fourier terms in Willmott and Ellington (1997) as illustrated in Fig. 2:

*t*is dimensionless and parameter

*K*is the reduced frequency defined by 2π

*fc*

_{m}/2

*U*

_{ref}, where

*f*is flapping frequency,

*c*

_{m}is the mean wing chord length (also used as reference length) and

*U*

_{ref}is a reference velocity at the radius of the second moment of wing area [

*r*

_{2}(s): as a fraction of wing length

*R*] (Willmott and Ellington, 1997). The coefficients ϕ

_{cn}, ϕ

_{sn}, θ

_{cn}, θ

_{sn}, α

_{Hcn}, α

_{Fsn}, α

_{Hcn}and α

_{Hsn}are determined from the empirical kinematic data as in Willmott and Ellington (1997) where

*n*is an integer varying from 0 to 3. The reference velocity equals , where is the mean angular velocity of the wing. is equal to twice the wingbeat amplitude (Φ, in rad) multiplied by the wingbeat frequency (

*f*, in s

^{−1}). The Reynolds number

*Re*is given by: where ν is the kinematic viscosity of air, 1.5×10

^{−5}m

^{2}s

^{−1}. For the hawkmoth,

*c*

_{m}=1.83 cm,

*r*

_{2}(s)=0.52,

*R*=4.83 cm, Φ=2.0 rad and

*f*=26.1 s

^{−1}, so

*Re*is between approximately 3000 and 4000. The reduced frequency parameter

*K*was calculated as 0.37, resulting in a dimensionless computational period of 8.5. The dimensionless parameters

*Re*and

*K*describing the flow were the same for the real moth, the flapper and the computational case.

### Geometry of the virtual wing and domain discretization

We base our geometry of the object on a tracing of the outline of a fore- and hindwing of a male *M. sexta*. The flapper wing shape was based on the same tracing (Van den Berg and Ellington, 1997), whose aspect ratio is 5.3. A uniform thickness was used, but with elliptic smoothing at the leading and trailing edges as well as at the tip. On the basis of the kinematics of the hawkmoth, in which the two wings do not clap but maintain a certain distance between them, we model a single virtual wing to mimic the hovering mode employed by insects that do not take advantage of the clap- and-fling mechanism.

Since the main spars of insect wings radiate from the bases of the wings, we constructed our virtual wing in a similar manner. Also, since a hawkmoth wing couple consists of a connected forewing and hindwing, a two-block method was developed in which each block of grids corresponded to one of the wings. The two parts were joined by a co-grid line at the co-axis between them (Fig. 3A).

To facilitate resolution of the viscous flow at the leading edge, trailing edge and wing tip as well as to save computing time, an O-O type grid topology was introduced (Fig. 3B). In this scheme, the whole physical domain is a semi-ellipsoid. A grid system with 72 000 (79×29×29) nodes was generated around the wing. Grids were clustered to the wing surface, in particular to the leading edge, the trailing edge and the wing tip. To avoid the unstable reflection of the solution at the open boundary, the distance far from the wing was six times the mean chord length.

To model the moving/deforming boundary during movements of the wing, a special method was developed to regenerate grids explicitly at each time step, which minimizes the computing time in regridding. To deal with large deformations of the wing, regridding is implemented in two steps: (1) rotation of the grids in the whole domain according to the ‘rigid’ rotational motion of the wing, i.e. flapping, elevation and feathering, which keeps the qualities of the original grids; followed by (2) nonlinear twisting of the grids to incorporate flexing and camber, and stretching the grid lines to fit the deformed wing completely. In the movements of the flapper (Van den Berg and Ellington, 1997), camber is formed as a result of the difference in angle of attack between the forewing and hindwing. This camber causes shear in the wing membrane, which is particularly pronounced during supination and pronation, and which is relieved by a small spanwise sliding motion of the hindwing. We realise this deformation numerically by first rotating the grids according to the angle of attack of the flexing forewing, and then flexing the hindwing to account for the difference between the angles of attack. To mimic the realistic wing, this flexing of the hindwing is achieved by moving grid points vertically (H. Liu, in preparation).

### Evaluation of forces

The Navier–Stokes equations, once solved, describe the near-field flows and the resultant wake in terms of the velocity components and pressure at discretized grid points or cell centres. The lift and drag forces acting on the wing surface can be evaluated from the pressure and stresses along its surface. The resultant lift and thrust forces are calculated first in the local wingbase-fixed coordinate system (*x, y, z*), as illustrated in Fig. 1A, and then transformed into the global space-fixed coordinate system (*X, Y, Z*) by accounting for the stroke plane angle.

*x*-,

*y*- and

*z*-directions (Fig. 1A); this can minimise the numerical error arising from integrating the projection of pressure and stresses. All the force-related quantities are nondimensionalized with respect to the reference velocity

*U*

_{ref}, the reference length

*c*

_{m}and the platform area of the wing

*S*

_{wing}such that: where

*L*is lift,

*T*is thrust,

*C*

_{L}and

*C*

_{T}are the lift coefficient and the thrust coefficient, respectively, and ρ is the air density. To assess the vertical force component that supports the insect’s weight, we need to transform the lift and thrust into the global coordinate system. Considering the inclination angle (β) of the stroke plane with respect to the horizontal plane, we can obtain the vertical (

*C*

_{v}) and horizontal (

*C*

_{h}) force coefficients respectively as and The force balance of hovering requires the mean vertical force to equal the weight, and the mean horizontal force to be zero. In addition to inserting the observed stroke plane angle β into equations 9 and 10, therefore, we can also calculate β such that: where and are the mean values of

*C*

_{L}and

*C*

_{T}, respectively.

## RESULTS AND DISCUSSION

### Validation

A series of validation studies have been undertaken and have shown that this computational model is robust and capable of reasonably estimating unsteady fluid–structure interactions in biomechanical problems (Liu *et al*. 1995, 1996; H. Liu, in preparation). To establish further the validity of the model in simulating the aerodynamics of hawkmoth hovering, we present below a comparison between our computed flows and the experimental results of Van den Berg and Ellington (1997). Details of this comparison will be published (H. Liu, in preparation).

Van den Berg and Ellington (1997) used video and 35 mm photographs to record the flow around the flapper using smoke released from the leading edge. The visualised flows were given at five spanwise cross sections at distances of 0.25*R*, 0.5*R*, 0.63*R*, 0.75*R* and 0.87*R* from the base, respectively. The core positions of the leading-edge vortices and the reattachment points are shown in Fig. 4. The flow at three positional angles ϕ during the downstroke with values of 30°, 0° and −36° is also illustrated in Fig. 4. For comparison, the computed flows were visualised at the same three positions, and the instantaneous streamlines are plotted in top and side views relative to the stroke plane (Fig. 5). Pressure contours over the wing are coloured: red indicates high pressure, while blue represents low pressure. There is good qualitative agreement between our CFD analysis (Fig. 5) and the experimental results (Fig. 4) in terms of three-dimensional leading-edge vortex flow. The CFD results show clear axial flow at the core of the leading-edge vortex during the middle and late downstroke, giving the vortex a spiral conical shape. Also, the negative pressure region over the wing surface near the leading edge matches quite well the location of the measured cores of the leading-edge vortices. Additionally, a spanwise velocity gradient can be seen in Fig. 6. The shed vortex from the early half of the downstroke can still be seen above the wing, which is in good agreement with the smoke-visualised flow in Fig. 4B (the red-coloured vortex). A strong tip vortex with large velocity magnitudes (Fig. 6) is detected at the position where the leading-edge vortex breaks down and is shed from the wing tip.

In contrast, the computed leading-edge vortex is slightly smaller than the experimental one (Figs 4, 5). In a study of the visualised airflow around a tethered hawkmoth in a wind tunnel (Willmott and Ellington, 1997), the leading-edge vortex was also observed to be small at low near-hovering airspeeds. Note that streaklines, not streamlines, were observed in the flapper and hawkmoth experiments.

We had no measured time course of instantaneous forces to validate the computed force-related quantities. With the good qualitative agreement of the flows between the simulation and the measurements, however, we think that the present CFD model can give a reasonable understanding of the complicated unsteady fluid phenomena in insect flight. In the next section, we address an overall understanding of how the viscous and unsteady flows near the wing behave and how the instantaneous forces are produced during the course of a complete flapping cycle.

### CFD-visualised flows during a flapping cycle

#### During downstroke

As illustrated in Fig. 5, a large leading-edge vortex is readily produced during most of the downstroke. During the early downstroke (Fig. 5A), a coherent leading-edge vortex (*LEV1*) stretches from the base to approximately 60–75% of the wing length. *LEV1* is quite two-dimensional in structure, although some axial flow is visible, with the swirls closely attached to the leading edge of the upper surface and with the core enlarging towards the wing tip.

As the wing moves to a horizontal position (Fig. 5B), *LEV1* becomes a large, spiralling vortex with strong axial flow at the core, towards the wing tip. The average angle of attack relative to the local flow at this instant is approximately 50°, ignoring the mean induced velocity. This is much higher than the angle, approximately 20–30°, at which stall is initiated for a stationary hawkmoth wing (Willmott and Ellington, 1997). *LEV1* is very probably created by dynamic stall. Note that *LEV1* is very three-dimensional in structure, consisting of a conical spiral, enlarging as it is swept along the wing. Under *LEV1*, there is a region of negative pressure (coloured blue) on the wing surface, which corresponds to the diameter and location of the vortex and which will augment the lift force (Fig. 5B).

The axial (spanwise) flow is probably closely related to the spanwise pressure gradient resulting from the spanwise velocity gradient along the flapping wing. The axial flow is accelerated by the spanwise pressure gradient to approximately 50–75% of the wing length, where a negative pressure island is observed (Fig. 5B). The low-pressure region is quite three-dimensional in structure where the spanwise pressure gradient reverses. This corresponds to a large iso-speed area, where the axial flow may reach its peak (Fig. 6) and then starts to decelerate due to a *reverse* pressure gradient (Fig. 5B). At approximately 75% of the wing length, the axial flow is reduced to zero (Fig. 5B). These results strongly support those of Van den Berg and Ellington (1997), who found that the largest axial flow velocities were comparable with the mean wing tip velocity and were reached at approximately 50–60% of the wing length.

Immediately after the middle of the downstroke, the core of *LEV1* breaks down at 75% of the wing length, and the tip region separates from the wing early in the second half of the downstroke. We think that this breakdown and separation of *LEV1* into the tip vortex is because of the instability of the flow with the reverse pressure gradient combined with the onset of wing deceleration after the middle of the downstroke. *LEV1* does not show a significant change from the base to the 75% position.

Towards the end of the downstroke, a new leading-edge vortex (*LEV2*) is created between the wing tip and the broken-down position (Fig. 5C). *LEV2* also shows a spiral axial flow, but directed towards the base, and it connects to the original *LEV1*. This *reverse* axial flow is apparently created by the *reverse* pressure gradient originally generated by the tip vortex. *LEV2* is strongly affected by the tip vortex and is relatively unstable; the reverse pressure gradient meeting the spanwise gradient leads to a weak axial flow.

As the wing approaches the end of the downstroke and starts supinatory rotation, *LEV2* begins to be pushed off the leading edge because of the deceleration of the wing, and the breakdown point of *LEV1* moves in to 50% of the wing length (Fig. 5C). Note that *LEV2* looks like a free vortex but is not shed, maintaining a constant vorticity and staying connected to the tip vortex because of the axial flow.

### During supination

When the positional angle of the wing reaches its minimum (−45°), i.e. at the bottom of the stroke (Fig. 7), the angles of attack of the fore- and hindwing have risen to 25° and 10° (Fig. 2), respectively, resulting in a cambered profile. The *LEV*s produced during the downstroke and the tip vortex appearing during early supination are not shed yet but, with a slight deformation, turn into a ‘hook-shaped vortex’ (Fig. 7B). The hook-shaped vortex stays attached to the leading edge and wraps around the wing tip, producing an island at the centre of the ‘hook’ where all the reattachment (and/or attachment) flows meet. Hence, there is a stagnation point at the centre, resulting in a positive pressure region. Additionally, a shear-layer vortex (*SLV*, a vortex created by the steep gradient of the shear layer close to the wing surface) is observed on the lower surface of the wing, running from the base to the wing tip, causing a region of negative pressure (Fig. 7C).

During the latter half of supination, when the wing translates upwards, rotates, flexes and increases its angle of attack, the hook-shaped vortex is strongly deformed and very quickly collapses beyond the mid-span joint of the two *LEV*s. *LEV1* is pushed onto the upper surface of the wing with a straight line of axial flow, but is virtually separated at the mid-span and connects to a large flux of flow going over the upper surface and shedding from the trailing edge. This shedding over the trailing edge during late supination is also observed in the flow-visualisation experiments of a tethered peacock butterfly *Inachis io* (Brodsky, 1991). *LEV2* rolls over the leading edge onto the lower surface, connects to the *SLV* and is eventually shed from the trailing edge.

### During upstroke

Immediately after shedding of the hook-shaped vortex at the end of supination, a very small leading-edge vortex appears at the wing tip at a positional angle of −36° (Fig. 8A). It is not separated, like a tip vortex. The vortex is quite two-dimensional in structure, although there is some spiral axial flow inside. We did not observe a separate starting vortex in the early upstroke. The leading-edge vortex is limited to a narrow region close to the wing tip and barely affects force generation at this instant (see below).

By the time that the wing reaches the middle position (Fig. 8B), the vortex extends from the wing tip to the base. The vortex still has a two-dimensional structure and is very small, probably because of the small angle of attack with respect to the local flow. The vortex is tightly attached to the leading edge but merely leads to a very small negative pressure region at the wing tip. The flow follows both the upper and lower surfaces smoothly. The axial flow at the core of the vortex initially travels from the wing tip to the base, but reverses direction during the latter half of the upstroke.

After the middle of the upstroke, the leading-edge vortex grows rapidly and hence enlarges the negative pressure region. When the wing reaches a positional angle of 51° (Fig. 8C), the vortex dimensions are comparable with those observed at early downstroke. Just as during the later part of the downstroke, a reverse pressure gradient and a reverse axial flow run from the break-down point of the upstroke leading-edge vortex to the wing tip. The vortex breaks down at approximately 60–70% of the wing length, but there is no significant shedding of the tip vortex, which is quite different from the downstroke. We believe that, as in the downstroke, the position of breakdown results from the balance between the spanwise pressure gradient and the reverse pressure gradient at the wing tip.

### During pronation

During pronation, the wing has a small velocity as the positional angle reaches its peak (65°) and the angles of attack of the fore- and hindwing approach zero (Fig. 2). The near-field flow around the wing is visualised in top and side views in Fig. 9. The upstroke *LEV* remains attached to the leading edge of the lower wing surface (Fig. 9A).

A trailing-edge vortex (*TEV*) is also detected (Fig. 9B,C) below the lower surface, which is probably due to wing rotation during the first half of pronation. This vortex runs from the base to the tip, is larger than the *LEV* and lies below it. Note that the *LEV* spirals to the wing tip, crosses the lower surface of the wing and eventually connects to the *TEV*. This leads to a large area of flow attachment and hence a positive pressure region (Fig. 9B). Furthermore, a shear-layer vortex (*SLV*) (Fig. 9C) considerably larger than that at supination covers the upper surface and leads to a large region of negative pressure (Fig. 9C).

As the wing continues to translate and rotate, changing its angle of attack, the separated region of the *SLV* gradually grows and moves towards the trailing edge, where it is eventually shed into the wake. As the *TEV* is shed into the wake, the *LEV* breaks down below the wing with a large flux of flow at the midwing. At the end of pronation, when the wing is about to start the downstroke, air flows smoothly over the lower surface but leaves a large messy flow from the shed SLV of the upper surface. This implies that all the vortices observed during pronation are shed before the downstroke and thus that, in order to enhance lift generation, the wing must produce a new leading-edge vortex on the downstroke.

### Estimation of force generation

To link the CFD-visualized flows to the time course of instantaneous forces, the vortex system is illustrated in Fig. 10. The time variation curves (Fig. 11) of the vertical and horizontal forces are given in the global space-fixed coordinate system (see Fig. 1A). The mean horizontal force must be zero in hovering, and substitution of the computed forces into equation 11 suggests that a stroke plane angle (β) of 23.6° is required to satisfy the force balance. This value is in excellent agreement with the range of approximately 20–25° reported by Willmott and Ellington (1997) for hovering *Manduca sexta* and it has been used in equations 9 and 10.

Comparing Fig. 10 with Fig. 11, we see that the shedding of all the vortices at the end of pronation apparently leads to a small decrease in vertical force. The first half of the downstroke shows a steep increase in lift, and the vertical force coefficient rises almost linearly to a peak of approximately 0.8. This is probably due to delayed stall, and the increasing size, and hence strength, of the leading-edge vortex (*LEV1* in Fig. 5A,B) produces an increasing negative pressure region. In contrast, the lift shows relatively non-linear features during the latter half of the downstroke (Fig. 11) but maintains a very high value. We believe that these forces are predominantly produced by the two leading-edge vortices *LEV1* and *LEV2* (Fig. 5C), which are also created by delayed stall. *LEV2*, with reverse axial flow originating from the wingtip vortex with strong non-linearity, rapidly increases its strength and hence its fractional contribution to lift production. Therefore, the aerodynamics of the flapping wing, particularly during the latter half of the downstroke, is strongly dependent upon the unsteady three-dimensional movements of the wing.

During the upstroke, the lift coefficient shows a moderate increase, from a slightly negative value immediately after supination to a peak of approximately 0.4 just prior to pronation (Fig. 11). A negative horizontal force, or thrust, is produced throughout the upstroke, which balances the backwardly directed horizontal force of the downstroke. During the first half of the upstroke, therefore, the dominant forces are very likely to be drag-based (both pressure and friction drag). In the late upstroke (Fig. 8C), we see a *LEV* with axial flow developing quite quickly, and this non-linearises the force variation (Fig. 11) as pronation is approached.

To investigate the extent to which the translational motion is responsible for force production during one complete flapping cycle, we calculated the mean vertical and horizontal forces during the downstroke and upstroke by integrating the instantaneous forces. Mean vertical and horizontal force coefficients are respectively 0.46 and −0.31 during the downstroke, and 0.14 and 0.31 during the upstroke. Note that most of the vertical force, nearly 77% of the total in one cycle, is produced during the downstroke, but 23% is still produced during the upstroke. The backward thrust force during the downstroke is completely cancelled out by a forward thrust force during the upstroke because of the stroke plane angle calculated using equation 11.

The complex vortex structure during the upstroke, with the *LEV*, the *TEV* and the *SLV* (Fig. 9), complicates force production during pronation. On the lower wing surface, the interaction between the *LEV* and the *TEV* results in a large region of positive pressure (Fig. 9C); on the upper surface, we see a large region of negative pressure at the *SLV* core. Note that the lift force remains positive during pronation (Fig. 11). With the calculated stroke plane angle of 23.6°, we see a large instantaneous vertical force coefficient of approximately 0.35 immediately after the highest wing position, which is even larger than the mean value of 0.3 during one complete cycle (Fig. 11). Because pronation is so brief, however, it contributes only slightly to lift production. Analysis of the aerodynamics during pronation suggests that the different three-dimensional movements (Fig. 2) of the wing during supination are of great importance in avoiding negative lift production. The asymmetric kinematics of supination creates a complex ‘hook-shaped’ vortex on the upper surface (Fig. 7B) with a small *SLV* on the lower surface (Fig. 7C), which is quite different from that seen during pronation (Fig. 9). Therefore, little negative lift and thrust forces are produced during supination (Fig. 11).

*U*

_{ref}=2.6 m s

^{−1}and the surface area

*S*

_{wing}=8.8 cm

^{2}of a real

*Manduca sexta*wing (Willmott and Ellington, 1997), the vertical force

*F*

_{v}is equal to: where the density of air ρ is taken as 1.225×10

^{−6}kg m

^{−3}. The total vertical force is therefore equal to 11.05 mN for one wing. The mean force required per wing is half the weight of the hawkmoth (mass 1.58 g) or approximately 7.8 mN. Hence, the calculated mean force in a complete flapping cycle is 1.41 times that required to support the weight of the hawkmoth. The insect body may interact with the flow, but it is very unlikely that this would result in an almost 41% reduction of the vertical force.

### Leading-edge vortices during up- and downstroke

Large leading-edge vortices that enhance the circulatory lift and that are likely to feature prominently in the delayed stall of insect wings are commonly observed in two-dimensional model experiments at appropriate Reynolds numbers (e.g. Newman 1977; Saharon and Luttges, 1989; Dickinson and Götz, 1993). However, these large two-dimensional leading-edge vortex mechanisms can only give approximate explanations of insect flight because they lack the axial flow generated by the more sophisticated three-dimensional mechanism that couples both three-dimensional motions and the geometry of the wing.

A leading-edge vortex is observed (Figs 4B, 5B) and a large lift is detected (Fig. 11) during the downstroke, but how is the leading-edge vortex created? Using the instantaneous velocity based on the kinematic data (Fig. 2), we first calculated the angles of attack with respect to the local flow at the span position *r*_{2}(*s*)*R* and plotted them in Fig. 12. These values do not include the effect of the downwash from the induced velocity field. The angle of attack is roughly flat during most of the downstroke with a value of approximately 50°, which would be reduced by approximately 15° by the downwash (Willmott and Ellington, 1997). It is during this period that the large vortices (*LEV1* and *LEV2*) are closely attached to the leading edge on the upper wing surface. Compared with the stall angle of attack (approximately 20–30°) of a hawkmoth wing in steady flow (Willmott and Ellington, 1997), this indicates that the leading-edge vortices are created by delayed stall. The constancy of the angle of attack further implies that it might be optimal for creating the leading-edge vortex (Fig. 5) and hence the large lift production (Fig. 11).

Why is the leading-edge vortex spiral? The pressure gradient due to the velocity gradient and/or spanwise centrifugal acceleration may well have been responsible for the observed axial flow. However, these conditions also apply to helicopter rotors and wind turbine blades, but large-scale spanwise flows have not been observed for them (De Vries, 1983). So why is the spiral leading-edge vortex not formed for rotary wings whereas it is clearly present for flapping wings?

On the basis of our analysis given above, we think that there may be two essential conditions for the existence of a three-dimensional vortex with axial flow: vorticity and a pressure gradient at the core of the vortex. Flapping insect wings at a large angle of attack can produce a strong leading-edge vortex very probably by delayed stall. At the low Reynolds numbers at which insects fly, this *LEV*, when it occurs, would usually maintain chordwise flow (flow from the leading-edge to the trailing-edge) without shedding, which would further increase its strength. The steep spanwise pressure gradient, which is probably due to the velocity gradient, can produce the strong axial flow that concentrates the *LEV*, preventing it from accumulating into a large vortex that would be unstable, as observed in previous two-dimensional experiments relevant to insect flight (Edwards, 1982; Dickinson, 1996).

In helicopter rotors and wind turbine blades, the first condition is probably not satisfied. The high Reynolds numbers (up to approximately 10^{6}) together with a constant small angle of attack will not result in a strong and stable vortex; even if a vortex were generated during dynamic stall at high angles of attack, it would not be able to resist the very strong chordwise flow and would immediately collapse. Thus, even though the second condition may be fulfilled, we observe only a small radial velocity component (Chee, 1990).

How is the reverse pressure gradient at the wing tip (Fig. 5B) produced? The large angle of attack of the wing usually leads to a steep velocity gradient at the wing tip and thus produces a pressure gradient between the lower and upper surfaces (Fig. 5B). As a result, as is often seen in flows around a delta wing, a wing tip vortex is observed. Differences in velocities at the wing tip during the downstroke therefore lead to a reverse pressure gradient against the spanwise gradient produced by the velocity gradient (Figs 5B, 6). Note that both flow and pressure distribution at the wing tip are very much three-dimensional in structure and therefore are strongly dependent upon both the local geometry and the wing motions. Combining with the spanwise pressure gradient, the reverse pressure gradient decelerates the axial flow to zero at approximately the 75% wing length (Fig. 5A,B). During the latter half of the downstroke, the wing decelerates and hence the velocity gradient decreases, leading to a reduced spanwise pressure gradient. In addition, with increasing positional angle, the velocity gradient over the wing tip between the upper and lower surfaces may become steeper, which could result in an increase in the reverse pressure gradient. As a result, the reverse pressure gradient becomes increasingly dominant and the point of breakdown is continuously pushed towards the base of the wing. At the end of the downstroke, it is located near the middle of the wing (Fig. 5C). This gives a good explanation of what happens in the dark region from the mid-span to the wing tip as seen in the smoke-visualised flow in Fig. 4C.

### Aerodynamics during pronation and supination

The presence of the *LEV* during pronation as well as supination poses a question: how can a *LEV* created during previous translational motion remain attached to the wing and not be shed? This is a general question that can be applied to the flapping flight of insects and other animals that do not take advantage of the clap- and-fling mechanism (Weis-Fogh, 1972). As described in previous sections, we believe that this is because of the appropriate coupling of unsteady motions involving the deceleration, the elevation and the rotational motion, which *in toto* delays shedding of the vortex during pronation and supination (Fig. 10).

Additionally, the appearance of the shear-layer vortex (*SLV*) can be explained by the fact that the rotational motion and deceleration of the wing lead to the separation point of the *LEV* rolling over the leading edge and attaching onto the upper surface, which means that the leading edge cuts into the vortex. Additionally, with the help of the elevation motion, flow within the separated region on the upper surface of the wing would be accelerated, which could also enlarge this region of suction and lead to the production of the large shear-layer vortex. Note that the vortex shows a three-dimensional conical structure. Therefore, the region of suction leads to the production of a large negative pressure region covering most of the upper surface of the wing which, together with the positive pressure region on the lower surface, may contribute to some extent to lift production during pronation, depending upon the orientation of the wing.

Through our CFD analysis, we believe that the rotational flow structures are apparently dependent on the *history* of the aerodynamics, i.e. the leading-edge vortex created during previous translational motion, and also on the complex three-dimensional motions of the wing. The complicated vortex system observed during pronation and supination could not occur in the absence of the upstroke and the downstroke. The rotational motion seems to capture and deform the existing vortex, but not create a new leading-edge vortex. For insects taking advantage of the clap- and-fling mechanism, the clap at the top of the stroke would probably shed the leading-edge vortex generated during the upstroke, and the quick fling of the wings would produce a large new leading-edge vortex on the upper surface that would enhance the circulation and lift on the downstroke. Although many insects could use the circulation created by a fling/peel motion, most seem to shun the mechanism in spite of the aerodynamic benefits. If nothing else, there must be significant mechanical wear and damage to the wings caused by repeated clapping (Ellington, 1995). For most insects that use a near fling or near peel, where the wings remain separated by a fraction of the chord length and do not actually touch, however, it is obvious that they receive aerodynamic benefit from the previous translational motion.

## Acknowledgement

This research was supported by the Japan Science and Technology Corporation and by the British Council. Thanks go to Dr Shigeru Sunada for valuable discussions and suggestions.

## References

*Inachis io*L. (Lepidoptera, Nymphalidae) and some aspects of insect flight evolution

*J. exp. Biol*

*Proceedings of the 46th Annual Forum of the American Helicopter Society*

*J. exp. Biol*

*Ann, Rev, Fluid Mech*

*Drosophila melanogaster*

*J. exp. biol*

*J. exp. Biol*

*J. Fluid Mech*

*Biological Fluid Dynamics*

*Nature*

*Expts Fluids*

*J. exp. Biol*

*Comp. Phys. Comms*

*Phil. Trans. R. Soc. Lond. B*

*Low-speed Aerodynamics: From Wing Theory to Panel Methods*

*FED*

*Numerical Developments in CFD*. ASME

*J. exp. Biol*

*J. Fluid Mech*

*Aeschna*dragonfly

*Scale Effects in Animal Locomotion*

*J. Fluid Mech*

*AIAA Paper*no

*AIAA J*

*Corvus monedula*) in slow flight

*J. exp. Biol*

*J. exp. Biol*

*Phil. Trans. R. Soc. Lond. B*

*AIAA J*

*Drosophila*

*J. exp. Biol*

*Schistocerca gregaria*(Orthoptera: Acrididae), flying in a wind tunnel

*J. Kansas ent. Soc*

*J. exp. Biol*

*Physiol. Zool*

*Manduca sexta*. II. Aerodynamic consequences of kinematic and morphological variation

*J. exp. Biol*