## ABSTRACT

This paper responds to research into the aerodynamics of flapping wings and to the problem of the lack of an adequate method which accommodates large-scale trailing vortices. A comparative review is provided of prevailing aerodynamic methods, highlighting their respective limitations as well as strengths. The main advantages of an unsteady aerodynamic panel method are then introduced and illustrated by modelling the flapping wings of a tethered sphingid moth and comparing the results with those generated using a quasi-steady method. The improved correlations of the aerodynamic forces and the resultant graphics clearly demonstrate the advantages of the unsteady panel method (namely, its ability to detail the trailing wake and to include dynamic effects in a distributed manner).

## Introduction

The growing body of research on the aerodynamic energy requirements necessary for insect flight (Weis-Fogh, 1956; Dudley, 1991; Ellington, 1991; Sunada *et al*. 1993*a*) utilizes six prevailing methods: momentum, blade-element, hybrid momentum (or vortex), lifting-line, two-dimensional thin airfoil and lifting-surface (or vortex lattice). Unfortunately, extant analyses of the aerodynamics of flapping wings are unduly constrained for reasons related either to the inherent limitations or to the inappropriate application of prevailing methods. In response to such constraints, this paper highlights and illustrates the advantages of a distinct type of lifting-surface method known as an unsteady aerodynamic panel method. An earlier draft of this paper was presented as a poster at the Society for Experimental Biology Symposium ‘Biological Fluid Dynamics’, University of Leeds, England, July 1994.

## The analysis of aerodynamic forces on flapping wings: prevailing methods

### Momentum, blade-element and hybrid momentum methods

The concern of biologists and zoologists with the energetics of animal locomotion has led to a focus on momentum, blade-element or ‘hybrid’ momentum methods (i.e. methods that combine elements of both simple momentum and blade-element theory), which enable the prediction of energy usage as a function of forward speed or animal size and type.

The momentum (or momentum-jet) method relies on modelling the beating plane of a flapping wing as an actuator disk which accelerates the surrounding air and thus imparts a change in momentum to it, generating a thrust which balances the weight and drag of the host body. The method is thus only able to determine gross values of the aerodynamic forces and power requirements given the value of the downwind farfield velocity. It has also been used to determine the induced velocity at the wings. In this method, no reference is made to the fact that the wings are beating or to the airfoil properties of the wings themselves. As Spedding (1992) points out, the method cannot reflect changes in wing area, aspect ratio, wing-beat frequency, section geometry or kinematics.

To establish some understanding of the aerodynamic forces on flapping wings, a simplification is made by assuming that the instantaneous forces developed by a flapping wing correspond to those in steady motion at the same instantaneous velocity and attitude; that is, the quasi-steady assumption. The usual quasi-steady aerodynamic treatment is based on blade-element theory, in which the wings are divided into a number of chordwise transverse sections for which the forces can be calculated. The relative velocity of each wing section is calculated from the sum of flapping, forward and induced velocities. The induced velocity is calculated from the momentum-jet or actuator disk model as developed by Ellington (1980).

Weis-Fogh (1956) and Jensen (1956) have found that the forces produced by locusts in forward flight are explicable using quasi-steady assumptions. Ellington (1980), however, queried the validity of Weis-Fogh’s assumptions and developed a vortex theory of flight by combining blade-element and momentum theories using a pulsed actuator disc to mimic the periodic beating of the flapping wings. This ‘hybrid momentum’ method has also been developed by Rayner (1979*a,b,c*) with a much more detailed analysis of the wake, as well as by Spedding *et al*. (1984) and Spedding (1986, 1987). In both Ellington’s and Rayner’s methods, only mean values of lift and power requirements of flapping wings are determined. Parameters such as animal weight, wing length, wing flapping angle and frequency do, however, enter into the calculations.

Another ‘hybrid momentum’ method is the ‘local circulation method’ developed by Kawachi (1981) in the analysis of helicopter rotors and wind turbines and applied by Azuma *et al*. (1985) to the forward flight of dragonflies. The loading on a wing is defined by combining the approach of the blade-element analysis with a more realistic and complete analysis of the modifying effects of the unsteady wake (Azuma *et al*. 1985; Azuma and Watanabe, 1988). The wake is defined by the path of the tip and trailing edge of the wing, and an iterative procedure is invoked to balance the effect of the wake and circulation distribution of the wings. The wing planform is approximated by the action of a series of elliptically loaded airfoils of diminishing size. Nonlinear, empirical lift coefficient (*C*_{L}) *versus* angle of attack (*α*) curves are used to compute the actual forces and moments on the wings for each blade element. This method is similar to the lifting-line method (to be discussed below) and is thus subject to the same limitations usually associated with such a method (i.e. if the flow is significantly three-dimensional, the validity of using two-dimensional lift coefficient data is questionable).

Biologists and applied mathematicians have also focused on the derivation of the energy usage of beating wings by analyzing particulars of the wake. Rayner (1979*a,b,c*) regards all the force developed by beating wings to be generated during the downstroke only (for forward flight) and neglects any loading on the upstroke. Thus, a series of discrete vortices is generated. With assumptions made regarding the core size of the vortex rings, the induced flow at the wings and the mean power expenditure can be approximated. Brodsky (1991) also hypothesized the nature of the trailing wake, assuming it to be like a crooked ‘ladder’ with the dominant circulation evident in the downstroke.

After careful examination of bird flight, Spedding *et al*. (1984) characterized a ‘concertina’ or ‘roller coaster’ wake. As pointed out by Lighthill (1990), this concertina wake takes the form of a pair of vortices of equal, constant circulation but with the distance between them varying periodically. Because the circulation is constant, the wide wake (shed during the downstroke) exerts a greater momentum than the narrow wake with a smaller momentum (shed during the upstroke).

Investigations of the energetics of animal locomotion have found that predictions of power expenditure do not correlate fully with in-flight measurements of metabolic rate (Spedding, 1992). To address this problem and clarify the breakdown of insect energetics, Wilkin and Williams (1993) have measured total in-flight forces (aerodynamic plus inertial) on a tethered sphingid moth and compare their results with those from a quasi-steady, blade-element analysis. Their comparison, however, is only qualitative and thus a comprehensive analysis of both the aerodynamic and inertial forces (and hence the energy expended during each wing-stroke cycle accounting for the complex kinematics of the wings) has yet to be developed.

### Lifting-line method

To obtain better temporal estimates of the thrust, lift, power and propulsive efficiency of flapping wings, a methodology more comprehensive than the momentum-jet or wake analysis methods is required. One such approach is the lifting-line method, which also attempts to delineate the dominant dimensionless parameters that describe the thrust performance of flapping wings. In their investigation of the aerodynamic loads and propulsive efficiencies of rigid, non-twisting flapping wings, engineers such as Betteridge and Archer (1974) implement an unsteady lifting-line method that uses assumed time-dependent pressure modes, where the root axis of the wings can be set at an angle to the free stream. Archer *et al*. (1979) use a similar procedure, allowing the wing to be torsionally flexible. Their analysis, however, is limited to the in-phase twist response of a flapping wing. The relationship between the aerodynamic and structural forces necessary to produce the wing twist is, therefore, not addressed. The induced velocity at the lifting line is developed from a quasi-steady model of the wake and excludes unsteady features of the vortex wake.

Phlips *et al*. (1981) also implement a lifting-line approach in their analysis of the aerodynamics of bird flight. The wake is modelled in two parts: a near wake, which is a vortex sheet formed from the streamwise vortices traced out by the wings, and a far wake consisting of discrete transverse and streamwise vortices. Their model, however, neglects the convection of the wake.

Ahmadi and Widnall (1986) developed a low-frequency unsteady lifting-line method for a harmonically oscillating wing of large aspect ratio, using matched asymptotic expansions. Guermond and Sellier (1991) extended this method by reducing the limitation on reduced frequency and also by applying it to swept wings. Both approaches are based on linearized aerodynamic theory and, as such, are restricted to small-amplitude transverse oscillations of wings.

While the use of this method has established that the trailing wakes generated by flapping insect or bird wings are important in understanding the forces on the wings and has facilitated hypotheses of the general structure of the trailing wake, the level of modelling remains crude. Prevailing models of the wake that use the lifting-line method tend to be three-dimensional adaptations of two-dimensional theory and consequently are unable to predict in detail the dynamic behavior of the trailing wakes due to the motion of the wings when undergoing large displacements and deformations. The accuracy of the assessment of the influence of the trailing vortices on the aerodynamic forces acting on the wing is thus compromised. In general, then, the lifting-line method has enabled a multidisciplinary array of researchers to model the aerodynamics of flapping flight and established the importance of general vortex effects. The assumptions of small-amplitude motions and a high-aspect-ratio wing inherent to a low-frequency, unsteady lifting-line method are, however, unduly restrictive for the purposes of modelling the aerodynamics of flapping insect wings.

### A two-dimensional method using thin airfoil theory

To analyze flapping wings, DeLaurier (1993) and DeLaurier and Harris (1993) have adapted a method widely used by engineers in the aerodynamic analysis of helicopter blades. The method is essentially a two-dimensional strip implementation of the Theodorsen function *C*(*k*) for sinusoidally heaving and pitching oscillating wing sections. DeLaurier uses a modified Theodorsen function, as developed by Scherer (1968), to account for wings of finite aspect ratio. This theory, although convenient to implement, is strictly only applicable to oscillations of a small order of magnitude. The method has been used successfully in the stability analysis of helicopter rotors, where the deflections of the rotors are of a small magnitude (Friedmann, 1983, 1987).

### Lifting-surface or vortex lattice method

In the lifting-line method, the importance of the trailing wake and the necessity of twisting of the wing in order to produce thrust are confirmed, but detailed geometric (i.e. spatial) and kinematic effects of the wing are suppressed. The lifting-surface method, in contrast, enables the wake and the wing to be represented by a lattice of vortex elements, thus permitting a more detailed representation of the wake and wing geometry.

Lan (1979), following the work of Albano and Rodden (1969), has developed a vortex lattice approach to the modelling of oscillating flat-plate wings. The propulsive efficiency and thrust for some swept and rectangular planforms are calculated for varying phase angles between the pitching and heaving motions. The method is applied to the study of tandem wings and it is shown that tandem wings can produce high thrust with high efficiency if the pitching is in advance of the flapping and the hindwing leads the forewing with some optimum phase angle. The motion of the wings is, however, limited to small-amplitude harmonic motion.

Sunada *et al*. (1993*a,b*) have also developed a vortex lattice approach to the modelling of flat plates. The method is applied to both the analysis of splitting triangular plates and the take-off of a butterfly. The wake, however, is artificially prescribed and the method relies on experimentally developed ‘shape factors’ to account for dynamic effects.

### Prevailing aerodynamic methods: conclusions

Overall, current theory and research on flapping flight which utilizes prevailing methods establishes that flapping flight is an aeroelastic phenomenon which is inherently dynamic, is characterized by rapid reversals in stroke direction and in wing rotation which result in gross movements of and between lifting surfaces, and produces the necessary aerodynamic forces for flight in a highly efficient manner. The prevailing methods variously establish that any attempt to model the aerodynamic forces on flapping wings must accommodate both trailing vortex effects and wing force resolution in a detailed manner.

A general problem, therefore, with existing methods is that while some can detail vortex effects and others can accommodate wing force resolution, not one of the methods reviewed above is capable of detailing both. For example, the hybrid method has no detailed wake or detailed force resolution, the lifting-line method has no detailed wake resolution, is valid only for small displacements and has no detailed force resolution, and the prevailing lifting-surface method has no detailed free-wake analysis. Hence, such complaints as ‘presently, a vortex theory based on a simple wake model for forward flight of insects is lacking’ are still valid (Spedding, 1992, p. 82). Moreover, existing studies which have attempted to compare the in-flight measurements of metabolic rates with predictions of power expenditure have not provided the necessary comprehensive analysis of the aerodynamic forces involved. Given the need to model the relevant aerodynamic forces on flapping wings, and the disadvantages of prevailing aerodynamic methods, the present study advances a type of lifting-surface method known as an unsteady aerodynamic panel method.

## The advantages of an unsteady aerodynamic panel method

The unsteady aerodynamic panel method is a classical boundary element method which relies on developing a distribution of source and doublet ‘singularities’ on wing and body surfaces, and doublet ‘singularities’ to represent the wake. To date, the panel method has been used, almost exclusively, to analyze the aerodynamic forces on aircraft (Ashley and Landahl, 1985; Ashby *et al*. 1988; Katz and Plotkin, 1991). The unsteady panel method is based on potential theory which assumes non-viscous flow (see Katz and Plotkin, 1991). It is valid for Reynold’s numbers of the order of 10^{4} and above (Reynold’s number *Re =cU/𝒱*, where *c* is a nominal chord length, *U* is based on the flapping frequency, a nominal radius and the free-stream velocity and *𝒱* is the kinematic viscosity of air). Engineers have established that the advantages of the panel method are that it accommodates the detailing of the trailing wake, includes dynamic effects and includes those effects in a distributed manner. Moreover, the panel method is also capable of accommodating flexibility and interference effects. Such advantages render the panel method more useful than other prevailing methods when analyzing the aerodynamic forces on flapping wings.

While, ideally, a comprehensive demonstration of the advantages of a panel method over other prevailing methods would encompass comparing results obtained from the application of each, the scope of such a project would be enormous. Instead, to demonstrate the method’s advantages, this paper provides a comparative example using the results of Wilkin and Williams (1993), who compared ‘derived’ experimental aerodynamic forces on a tethered sphingid moth with those predicted by quasi-steady theory. They measured total forces (aerodynamic plus inertial) and obtained ‘derived’ experimental aerodynamic forces by subtracting estimates of the inertial forces from the total measured forces. It was noted above that, while the experimental data of Wilkin and Williams (1993) are useful, the correspondence they obtain between theory and measurement is only qualitative. Such results are not surprising, given that quasi-steady theory does not accommodate unsteady fluid-flow effects (such as the development and subsequent influence of free vortices) and that such effects are integral to flapping-wing flight (see Ellington, 1984; Brodsky, 1991; Spedding, 1992; Sunada *et al*. 1993*b*).

Wilkin and Williams (1993) provide wing-angle-to-plane trajectories, moth wing planform, the wind-tunnel wind velocity (3.36 m s^{-1}) and derived aerodynamic vertical and horizontal forces acting on a tethered sphingid moth. It is important to note that the latter forces are derived by subtracting estimates of the inertial forces from measured data that include both the aerodynamic and inertial forces. Drawing on their data, the present study sets out the assumptions made in modelling the in-flight aerodynamic forces of a tethered sphingid moth, details the necessary coordinate systems used, their interrelationships and kinematic relationships, and outlines a relevant potential flow model. The boundary conditions are then developed by determining the trajectories of the wings and their angular rates of change from experimental wing-angle-to-plane data. The discretization procedure is then delineated and the velocity components, pressures, loads and trailing vortex roll-up are calculated. Comparative results are then presented and general conclusions drawn.

## Modelling the flapping wings of a tethered sphingid moth using an unsteady aerodynamic panel method

Given that Reynold’s numbers (*Re*) for the aerodynamic flows of the sphingid moth are of the order of 10^{4}, this *Re* range suggests that the airflows are inertially dominated and, as such, viscous effects should not be directly important in most of the flow field. Viscous effects are thus confined to regions near the body, the wing and the wake that is shed behind the body or the wing. In contrast to high-speed flows of aeronautical interest (*Re* of the order of 10^{7}), the flow in the immediate vicinity behind the insect is considered to be laminar, only beginning to dissipate or diffuse after some distance. The diffusion of the flow is therefore assumed to have a negligible effect on the aerodynamic forces. In an analysis of the wake behind a sphingid moth, A. P. Willmott (personal communication) showed that it was clearly discernible for at least 1 wing-beat period. An initial step, therefore, is to construct an inviscid, three-dimensional, potential unsteady solver which amounts to solving Laplace’s equation in space, with impermeability conditions imposed on the moving wing surfaces. A current panel method computer program that meets the requirements of solving Laplace’s equation, as well as detailing the wake, and which is used in this investigation, is PMARC (Panel Method Ames Research Center) (Katz and Plotkin, 1991; Ashby *et al*. 1988). While PMARC does not include dynamic effects, it can be modified to do so by developing the theory as outlined below.

To model the aerodynamic forces on the flapping wings of a tethered sphingid moth using an unsteady panel method, it is therefore necessary to take the following steps: (1) develop a potential flow model, (2) discretize the wings and the wake, (3) compute the velocity components, pressures and loads acting on the wings and (4) develop the boundary conditions. A summary of these steps is found below; for a more detailed derivation of the aerodynamic forces, see Smith (1995) and Katz and Plotkin (1991).

### Potential flow model

*S*

_{B}enclosed in a volume

*V*, with the outer boundary

*S*

_{∞}and a wake model

*S*

_{W}, where

*S*

_{W}models a surface across which a discontinuity in the velocity potential or the velocity may occur. The boundary conditions apply to

*S*

_{B}and

*S*

_{∞}, respectively.

_{i}for solid bodies

*S*

_{B}and a far-field velocity potential Ф

_{∞}), the equation for the velocity potential Ф at any point

*P*in the flow field becomes (Katz and Plotkin, 1991): where

*r*is the distance from the point of interest in the flow field to a singularity on a velocity potential boundary surface,

*S*is the velocity potential boundary area, and is the normal unit vector on the inside surface of the body, wing and far-field velocity potential boundaries. This equation provides the value of Ф(

*P*) in terms of Ф and ∂ Ф/ ∂

*n*on the boundaries.

*n*indicates differentiation in the normal direction at the velocity potential surface,

*μ*is called the doublet strength and

*σ*the source strength, and letting Ф

_{∞}= Ф

_{i}(since the value of the potential inside

*S*

_{B}may be any arbitrary constant) and

*ϕ*

_{p}= Ф − Ф

_{∞}, then, on the inside surface of the body: The problem is thus reduced to finding the appropriate distribution of source (∂ Ф/ ∂

*n*) and doublet strengths (Ф) (or ‘singularities’) over the surface of the body of interest. The normal velocity at the boundary is satisfied directly by a source strength distribution: where (the kinematic velocity of the boundary surface) is a consequence of the boundary conditions to be developed below. Once the surface distribution of sources is thus determined, equation 6 provides a means to determine the unknown doublet distribution: For a given set of boundary conditions, the solution is not unique and a Kutta condition must be applied to the wakes leaving the body (wing) (Katz and Plotkin, 1991).

### Discretization, computation of velocity components, pressures and loads

*σ*) are known, and the strengths of all previous wake panels (

*μ*

_{W}) are also known from calculations for previous time steps, both may be transferred to the right-hand side of the equation and in matrix form may be written as: where

*B*

_{JK},

*C*

_{JK}and

*C*

_{JL}represent the velocity potential influence coefficients per unit singularity strength for body panel

*K*or wake panel

*L*acting on the control point of panel

*J*. The wing panelling arrangement for the combined wing model is shown in Fig. 1. A more detailed discretization does not add significantly to the numerical accuracy.

*Q*at a collocation point

*k*is the sum of the kinematic velocity plus the perturbation velocity: where

*U*(

*t*),

*V*(

*t*) and

*W*(

*t*) are the reference velocities in the

*a*

_{1},

*a*

_{2}and

*a*

_{3}directions, respectively,

*l*

_{k},

*m*

_{k}and

*n*

_{k}are the local panel coordinate directions and

*q*is the local perturbation velocity.

*c*

_{p}may be computed from Bernoulli’s equation (Katz and Plotkin, 1991): where

*Q*and

*p*are the local fluid velocity and pressures, ∂

*ϕ*

_{p}/ ∂

*t*=− ∂

*μ*/ ∂

*t*since Ф

_{i}is a constant,

*p*

_{ref}is the far-field reference pressure,

*ρ*is the density of air and

*V*

_{ref}is the magnitude of the translational velocity of the origin

*O*

_{b}.

### Description of coordinate systems

To develop the governing equations for the kinematics, and hence the boundary conditions, the definition of coordinate vector bases sets is required. The chosen bases are the body axes {*b*} and the wing axes {*a*} (Fig. 2).

*Body axes* {*b*}

Origin *O*_{b} is fixed at some reference point in the body of the moth (not necessarily the centre of gravity). The axes are defined by *b*_{1}, *b*_{2} and *b*_{3}, with the *b*_{1}-axis positive aft (parallel to the free-stream wind velocity), the *b*_{2}-axis positive starboard, and the *b*_{3}-axis positive vertically up (Fig. 2).

*Wing axes* {*a*}

Origin *O*_{a} is fixed at the wing root, rotates with the wing and is used to define ‘local’ wing rotations (azimuth, flap and pitch). The axes are defined by *a*_{1}, *a*_{2} and *a*_{3}, with the *a*_{1}-axis positive aft, the *a*_{2}-axis positive starboard, and the *a*_{3}-axis positive vertically up (Fig. 2).

### Interrelationship of the coordinate systems

*b*} and {

*a*} bases is given by

*T*

_{ab}: where

*ψ*, is the wing azimuth angle, the horizontal angle between some reference direction (e.g.

*b*

_{1}-axis) and the projection of the

*a*

_{1}-axis on the horizontal plane (positive rotation is from north to west);

*θ*is the wing flap angle, the vertical angle between the

*a*

_{1}-axis and the horizontal plane (positive rotation is up);

*ϕ*is the wing pitch angle, the angle between the

*a*

_{1}

*a*

_{3}-plane and the vertical plane containing the

*b*

_{1}-axis; positive rotation is clockwise about the

*a*

_{1}-axis, looking forward; and

### Kinematic relationships

### Development of boundary conditions

Because the experimental wings-to-planes angles are not in Euler angle form, the experimental angles must be converted before modelling is possible. With reference to Figs 2 and 3, this conversion can be described as follows. Defining the line joining the ‘root’ (*R*) and the tip of the wing (*T*) to be the *RT*-axis, the angle between the projection of the *RT*-axis on the *b*_{3}*b*_{2}-plane and the *b*_{2}-axis is *v*. The angle between the projection of the *RT*-axis on the *b*_{1}*b*_{2}-plane and the *b*_{2}-axis is *h*. The angle *h* serves as a measure of the Euler angle *θ*. From geometric considerations:

*ϕ*. The rotation of the wing about the

*RT*-axis is obtained from measurements of the projected chord of the wing on the

*b*

_{l}

*b*

_{2}-plane and the actual wing chord at different wing-sections. This rotation angle serves as a measure of the Euler angle

*θ*at the root of the wing and also as a measure of the twist of the wing along its span. In this study, both the root value of the wing twist and an average value of the root wing twist and mid-span wing twist were investigated. Thus, the wing-angles-to-planes data are converted into the Euler angles

*ψ*,

*ϕ*and

*θ*. The derivatives of these angles with respect to time may be developed numerically. Using the kinematic relationships described above, the angular velocities

*p*

_{a},

*q*

_{a}and

*r*

_{a}are developed. The determination of the normal kinematic velocity on the wing surfaces and hence the boundary conditions then follows. A plot of the derived Euler angles as a function of the wing-beat cycle is shown in Fig. 4. The motion of the wing is described using the terms ‘supination’ and ‘pronation’, ‘downstroke’ and ‘upstroke’, ‘forwardstroke’ and ‘backstroke’ in correspondence with the wing motion between the extreme values of the Euler angles

*θ, ϕ*and

*ψ*, respectively.

## Modelling assumptions

To construct a tractable mathematical model of the in-flight aerodynamic forces of a tethered sphingid moth, the following assumptions are made (these assumptions generally allow for the neglect of flexibilities, separation effects and aerodynamic interference effects between bodies, the primary concern at the present time being the incorporation of dynamic effects and inclusion of the wake). (1) The fore and aft wings are treated as one combined wing, which is considered to be a flat rigid surface hinged by a ‘universal joint’ at a ‘root’ location. (2) The wing vortices are generated at the trailing edge only. (3) The flow behind the moth is considered laminar with the vortices having no time to dissipate under the influence of viscous effects (Grodnitsky and Morozov, 1993). (4) The rounded leading edges (veins) of the wings and the ‘corrugated’ profile of the wing surface inhibit leading-edge separation (Rees, 1975). (5) Each combined wing acts independently of the other. (6) The effect of the body is not included.

## Vortex wake

Since the wake is force-free, each wake panel moves with the local free-stream velocity. This velocity is the result of the kinematic motion and the velocity components induced by the wake and the body. A view of the wake developed over one cycle is shown in Fig. 5. Running the simulation for two wing-beat cycles did not alter the results appreciably, indicating that the wake has minimal impact after one or two wing-beat periods. Indeed, if one limits the length of the wake from a half wing-beat period to a whole wing-beat period, the results are not appreciably altered (see Figs 6, 7; Tables 1, 2). This confirms the third assumption above and also justifies the use of a potential flow model.

## Results and conclusions

It was asserted above that the advantages of the panel method over other methods are its ability to accommodate the detailing of the trailing wake, to include dynamic effects and to include such effects in a distributed manner. Fig. 5 graphically captures these advantages.

Figs 6 and 7 present calculated vertical and horizontal aerodynamic forces obtained using both an unsteady panel method (this study) and a quasi-steady method (Wilkin and Williams, 1993) and compare those forces with the ‘derived’ experimental results of the latter study. Table 1 presents average values of the aerodynamic forces developed by the two methods over the downstroke, the upstroke and the full wing-beat cycle.

Table 2 reveals that the correlations between the vertical forces obtained using the panel method are greatly improved over those of Wilkin and Williams’ (1993) quasi-steady results. Admittedly, the correlation coefficients for the horizontal forces are not improved. These forces are, however, of a lower order of magnitude than the vertical forces.

Fig. 6 also reveals that the calculated vertical force is in excess of the experimental value in both the upstroke and the downstroke for the case where the wing twist is taken as the wing twist of the root. This excess could be due to the following effects: (1) the neglect of flexibility effects in the model (i.e. the variation of twist along the span), (2) the neglect of leading-edge separation effects in the aerodynamic model or (3) interference effects of the body and wings. Another cause could be (4) the accuracy of the experimental data which Wilkin and Williams (1993) derived by subtracting their estimates of the inertial forces on the wings from measured data. In order to gain some measure of the effect of twist, another case is investigated where the wing twist is considered to be the average of the twist at the root and the twist at mid-span. It is noted in Fig. 6 that, while this case develops extreme forces much closer to the experimental values, average values of lift are, however, slightly compromised.

On closer examination of Fig. 6 and Table 1, secondary effects can be noted, due to the inclusion of the extended wake in the analysis (in this instance the wake for one wing-beat period). Secondary effects include the lowering of the extreme lift force on the downstroke to correlate better with the experimental result, as well as slight increases in lift at the beginning of the downstroke and the beginning of the upstroke.

Some general limitations of the present study are that it makes certain assumptions about the morphology of moth wings which depart somewhat from biological reality. In order to make the present analysis tractable, the fore and aft wings were combined as one and assumed to be rigid even though they are in fact flexible. The present results indicate that the experimental results lie between the two sets of rigid twist motion explored (see Figs 6 and 7, downstroke vertical force) and, as such, imply that there is a ‘mean value’ of wing twist of rigid motion that would give improved results. Moreover, were the wings modelled as flexible surfaces, it is likely that the present results would be further improved. Indeed, if the fore and aft wings were modelled independently of one another, this too would introduce more flexibility and further enhance the results.

Regarding the wake structure, it is assumed that the diffusion of the vortices behind the insect does not influence the aerodynamic forces appreciably, since the wake influence dies away as it moves further from the insect. This assumption is confirmed by comparing the results of the aerodynamic forces for differing lengths of the wake (see Figs 6, 7; Tables 1, 2).

Overall, the present analysis establishes that there are distinct advantages to the use of an unsteady panel method in that it includes both wake and distributed dynamic effects and, in comparison with other prevailing methods, better models the aerodynamic forces on rigid flapping wings. The panel method also holds great promise in its ability to accommodate other effects such as flexibility and interference.

## List of symbols

*A*rigid appendage (wing)

*a*_{1},*a*_{2},*a*_{3}basis unit vectors {

*a*}forming appendage or wing reference axes*B*_{JK}source strength potential influencing coefficient matrix of wing or body

*B*rigid host body to

*A**b*_{1},*b*_{2},*b*_{3}basis unit vectors {

*b*} forming host body reference axes*C*_{JK}doublet strength potential influencing coefficient matrix of wing or body

*C*_{JL}doublet strength potential influencing coefficient matrix of wake

*C*(*k*)Theodorsen function

*C*_{L}lift coefficient

*c*nominal chord length

*c*_{p}pressure coefficient on aerodynamic panel

*h*the angle between the projection of the

*RT*-axis on the*b*_{1}*b*_{2}-plane and the*b*_{2}-axis*J*control point of aerodynamic panel

*K*body panel

*k*collocation point

*L*wake panel

*l*_{m},*l*_{n},*l*_{k}local panel coordinate directions

*n*wing-beat frequency

normal unit vector on inside surface of body, wing and far-field velocity potential boundaries

*O*_{b}origin of body axes system

*O*_{a}origin of appendage axes system

*P*point of interest in the flow field

*p*fluid static pressure

*p*_{a}relative rotation rate of the {

*a*} axes system to the {*b*} axes system about the*a*_{1}axis*p*_{ref}fluid reference static pressure

*Q*total velocity

the velocity of the fluid

*q*_{a}relative rotation rate of the {

*a*} axes system to the {*b*} axes system about the*a*_{2}axis*q*_{l}fluid perturbation velocity in panel

*l*direction*q*_{m}fluid perturbation velocity in panel

*m*direction*q*_{n}fluid perturbation velocity in panel

*n*direction*RT*-axis line joining*R*and*T**R*root point of wing

*Re*Reynolds number =

*cU*/*v**r*distance from point of interest in the flow field to singularity on a velocity potential boundary surface

*r*_{a}relative rotation rate of the {

*a*} axes system to the {*b*} axes system about the*a*_{3}axis*S*velocity potential boundary area (=

*S*_{B}+*S*_{W}+*S*_{∞})*S*_{B}velocity potential boundary area of body

*S*_{W}velocity potential boundary area of wake

*S*_{∞}velocity potential boundary area of far field

*t*time

*T*tip point of wing

*T*_{ab}coordinate transformation relationship between {

*b*} and {*a*} axes*T*_{e ωa}kinematic relationship between and

*U*(*t*)reference velocity in the

*a*_{1}direction*U*_{∞}fluid free-stream velocity in

*x*-direction*u*_{l},*v*_{l},*z*_{l}wake velocity at each time step

*V*volume of potential flow interest

*V*(*t*)reference velocity in the

*a*_{2}direction*V*_{ref}reference free-stream velocity in

*x*-directionkinematic velocity of the boundary surface

*v*the angle between the projection of the

*RT*-axis on the*b*_{3}*b*_{2}-plane and the*b*_{2}-axis*W*(*t*)reference velocity in the

*a*_{3}direction*α*wing angle of attack

aerodynamic load on aerodynamic panel

*k*aerodynamic force at

*P*_{s}*ΔS*_{k}area of aerodynamic panel

*k**μ*doublet strength

*μ*_{W}strength of wake panels

*𝒱*kinematic viscosity of air

Euler angle rotation rates between the {

*a*} and {*b*} axes =relative angle rotation rates between the {

*a*} and {*b*} axes = [*p*a*q*a*r*a]^{T}- Ф
fluid velocity potential

- Ф
_{i}internal body fluid velocity potential

- Ф
_{∞}far-field fluid velocity potential

*ϕ*_{p}fluid velocity potential relative to Ф

_{∞}*ϕ*body

*A*(or wing) flap Euler angle*ψ*body

*A*(or wing) azimuth Euler angle- ρ
air density

*σ*source strength

*θ*body

*A*(or wing) pitch Euler angle*θ*_{exp}experimentally derived pitch Euler angle of wing

- ∇
**·**differentiation with respect to time, e.g. differentiation in the normal direction at a

boundary velocity potential surface

## ACKNOWLEDGEMENTS

The authors would like to thank the anonymous reviewers for all their helpful comments.

## References

*J. Fluid Mech*

*AIAA J*

*Aeronaut. J*

*NASA, TM*

*Aerodynamics of Wings and Bodies*

*J. exp. Biol*

*J. exp. Biol*

*Aeronaut. Q*

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

*J. exp Biol*

*Aeronaut. J*

*Aeronaut. J*

*J. exp. Biol*

*Beitrag zu Struktur und Funktion biologischer Antriebsmechanismen*. Instationare Effekte an schwingenden Tierflugeln

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

*J. exp. Biol*

*Vertica*

*Vertica*

*J. exp. Biol*

*J. Fluid Mech*

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

*NASA, TM*

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

*J. Fluid Mech*

*Aeronaut. J*

*J. Fluid Mech*

*J. exp. Biol*

*J. Fluid Mech*

*J. Fluid Mech*

*Nature*

*Experimental and Theoretical Investigation of Large Amplitude Oscillating Foil Propulsion Systems*

*Corvus monedula*in slow flight

*J. exp. Biol*

*Falco tinnunculus*in flapping flight

*J. exp. Biol*

*Adv. comp. env. Physiol*

*Columba livia*) in slow flight

*J. exp. Biol*

*J. exp. Biol*

*J. exp. Biol*

*Schistocerca gregaria*)

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

*J. Physiol. Zool*