## Abstract

Animal flight performance has been studied using models developed for man-made aircraft. For an aeroplane with fixed wings, the energetic cost as a function of flight speed can be expressed in terms of weight, wing span, wing area and body area, where more details are included in proportionality coefficients. Flying animals flap their wings to produce thrust. Adopting the fixed wing flight model implicitly incorporates the effects of wing flapping in the coefficients. However, in practice, these effects have been ignored. In this paper, the effects of reciprocating wing motion on the coefficients of the fixed wing aerodynamic power model for forward flight are explicitly formulated in terms of thrust requirement, wingbeat frequency and stroke-plane angle, for optimized wingbeat amplitudes. The expressions are obtained by simulating flights over a large parameter range using an optimal vortex wake method combined with a low-level blade element method. The results imply that previously assumed acceptable values for the induced power factor might be strongly underestimated. The results also show the dependence of profile power on wing kinematics. The expressions introduced in this paper can be used to significantly improve animal flight models.

## 1. Introduction

In the past, animal flight performance has been studied in various ways, often borrowing basic principles used in the design and operation of aeroplanes and helicopters. To a large extent, this will result in good approximations, as much of the physics is very similar. However, aeroplanes use separated systems for propulsion and lift production, which results in a model for energy use that is relatively simple. Helicopters use a mechanism that combines lift and thrust production, which is more similar to flying animals, but animals use reciprocating wings, rather than rotating wings. This adds another level of complexity, where besides the wingbeat frequency also the wingbeat amplitude is of importance for optimization. This work involves a theoretical study on the effects of wingbeat frequency, amplitude and stroke-plane angle, on the aerodynamic efficiency of animal flapping flight. In addition to providing an accurate aerodynamic characterization of self-powered animal flight, the utility of models like the one presented here is fundamental to understand adaptive flight behaviours in wild animals, ranging from display, foraging, food transport and migration [1]. For example, the power–speed relationship predicts alternative ‘optimal’ flight speeds that should be selected depending on the ecological context, and observations support the notion that birds and bats are capable of selecting an appropriate context-related flight speed in broad agreement with predictions [2,3].

Aerodynamic power is commonly separated into components with different underlying mechanisms:
*P*_{ind}, is due to the active acceleration of mass flow in order to produce a force opposing weight and drag (for fixed wing aircraft, the induced power owing to thrust production is considered separately as part of engine efficiency). Profile power *P*_{pro} is due to the local drag on the wings, and body power *P*_{body} due to drag on the body. In 1968, Pennycuick [5] applied this model to estimate flight performance in birds. An actuator disc model was used to obtain induced power, whereas blade-element theory, using quasi-steady aerofoil properties, was used to compute profile power. The energetic model was then expressed in terms of weight, wingspan, body frontal area, air density and air speed as a general model for bird and bat flight. For each component, there is a modifying factor, which includes all effects not accounted for by these measurable variables. These factors have been adjusted over the years following field observations of flight speeds [6,7] and wind tunnel experiments [8]. The effects of wing kinematics should thus be implicitly incorporated in these correction factors. However, little effort has gone into systematically quantifying the effects of reciprocating wings.

The effects of wing kinematics have been studied before, and several alternative flight models have been developed. The flapping wings are simulated at different levels of detail. A blade-element approach, similar to [5], was used by Parslew & Crowther [9] to optimize wing kinematics of a pigeon for minimal power. Rayner [10] used a continuous vortex lifting line model for forward flight, accounting for unsteadiness owing to flapping and wake deformation. Phlips *et al.* [11] used a similar lifting line model, but with an active downstroke and an inactive upstroke. Hall & Hall [12] developed a method involving a sheet of continuously varying vorticity shed from the flapping wings. This generalized vortex wake method allows for a continuous transition between wake topologies, resembling those observed in wind tunnel studies on birds [13]. The output of the Hall and Hall model is the ideal distribution of circulation over the shed wake for minimum power for a prescribed cyclic wake geometry. This ideal circulation distribution relates to the ideal lift distribution over the wing, which for steady planar wings is the familiar elliptical distribution [14]. Salehipour & Willis [15] used this vortex wake method to investigate an extended parameter space for animal flight, varying reduced frequency, stroke amplitude, thrust and lift. They found that wingbeat kinematics, specifically amplitude and frequency, could be optimized for minimum power. These studies all show that wing kinematics affect the power requirements for flight.

The aim of this work is to formulate explicit expressions that describe the effects of flapping wings on the different components of aerodynamic power in forward flight (equation (1.1)). These expressions are generally applicable to flying vertebrates, i.e. birds and bats, that use a span reduction during the upstroke. Using the Hall and Hall model [12] a set of solutions was produced for different combinations of thrust ratio, reduced frequencies and stroke-plane angles, where the wingbeat amplitude was optimized for minimum induced power (§2). We then used nonlinear regressions to express the effects of wing flapping on the drag and power components as functions of those three parameters (§3). Finally, the expressions are applied to an example to illustrate the implications of the model (§4).

## 2. Wake optimization

The approach we used for computing aerodynamic power is modified from the Hall and Hall model [12]. Here, we will first briefly describe the method of finding the optimal circulation distribution for a prescribed wake geometry. This is followed by a description of the wake geometry and how it is optimized.

### (a) Aerodynamic power optimization for a given geometry

The optimal wake method optimizes the distribution of circulation on a predefined wake geometry, as shown in figure 1. This can be interpreted as a sheet of vortex rings shed at each location the wing passes through the air. The strength (circulation) of each vortex ring is contained in the vector ** Γ**. Then, an influence matrix

**K**

_{ind}is constructed, describing the induced velocity owing to each vortex ring on control points related to all other vortex rings. This matrix is computed using the Biot–Savart law for a unit strength vortex for each ring using the prescribed geometry of the wake:

*j*,

**n**

_{j}is the normal to the enclosed plane and

*A*

_{j}its area. The velocity contribution of vortex ring

*i*on the control point (centroid) in vortex ring

*j*is indicated by

**V**

_{i,j}. For the computation of

**K**

_{ind}, the wake is repeated upstream and downstream of the wake under consideration, representing succeeding and preceding wingbeats. Formally, the wake needs to be repeated to infinity both upstream and downstream, but tests for convergence indicated the largest error on the induced power for using only four repeated wakes on each side was approximately 0.5%. With matrix

**K**

_{ind}, the induced power can be computed,

*ρ*is the air density and

*f*is the wingbeat frequency. The distribution

**is still unknown, but it has to satisfy the constraint that the wake should correspond to a prescribed aerodynamic force. This force can be expressed as**

*Γ***F**=

*ρf*

**B**

**, where the vector**

*Γ***F**contains the thrust force

*T*, a side force

*F*

_{y}, and a lift force

*L*. Matrix

**B**essentially contains the vortex ring area and orientation:

**for which ∂**

*Γ**P*

_{ind}/∂

**=**

*Γ***0**. As the force constraint has also to be satisfied, the equation that is minimized is formulated as

**F**

_{R}is the required force vector to be produced and

**λ**is a Lagrange multiplier. By setting the derivatives with respect to

**and**

*Γ***λ**to zero, the solution

*ϑ*(

*t*), stroke-plane angle

*φ*, the wingspan

*b*(

*t*) and the travelled distance per wingbeat

*U*/

*f*, with

*U*being the flight speed (figure 1). The ratio between wingspan and travelled distance is found in the reduced frequency

*k*

_{f}=2

*πfb*/

*U*, which is therefore a useful parameter to characterize the wake. The solutions for thrust (

*T*) and lift (

*L*) are linearly independent. This means that the solutions can also be characterized by the thrust ratio

*T*/

*L*.

In successive work on the Hall and Hall model [16], the model was extended to include profile drag. Aerofoil viscous drag may be approximated by the equation *V* is the local velocity, *c* the local chord, *C*_{Dv,0} the zero-lift drag coefficient, *C*_{Dv,2} the lift-dependent drag coefficient and *C*_{ℓ} the sectional lift coefficient. This is referred to as viscous drag, as the coefficients are Reynolds number dependent. The lift coefficient can be expressed as *C*_{ℓ}=2*Γ*/*Vc*, which makes it possible to include the lift-dependent viscous drag in equation (2.1) by extending the influence matrix
*c*_{j} is the chord length of the wing at the span position corresponding to vortex ring *j*. By using this formulation, lift-dependent viscous drag is taken into account during the optimization of the circulation distribution. The modified circulation distribution also affects the induced power and induced drag, but because this change is due to adding viscous properties, it will be counted as lift-dependent viscous power, i.e.
*Γ*_{vis} is the optimal circulation distribution including the viscous model, and *Γ*_{ind} the solution of the inviscid model.

Several studies have measured the drag of birds [17–20], bats [21] or their wings [22–25]. However, few studies have explicitly described the profile drag polar. In the viscous Hall and Hall model [16], a value of *C*_{Dv,2}=0.012 was taken from empirical measurements on a NACA 4412 aerofoil for a Reynolds number of 6×10^{6}. This Reynolds number is too high to be representative for animal flight [4]. In the work of Salehipour & Willis [15], values for *C*_{Dv,2} were estimated based on numerical computations on a NACA 0006 aerofoil. Over a range of Reynolds numbers 10^{3}<*Re*<10^{5}, values varied between 0.01 and 0.05 depending on which angle of attack was chosen for the quadratic fit. In the computations for the current work, we use a value *C*_{Dv,2}=0.03, which is in the middle of the above range (The resulting drag polar also compares well with experimental data from [20]; see the electronic supplementary material, S1).

Zero-lift drag is independent of the distribution of circulation and therefore not taken into account in the wake optimization (equation (2.1)), but its contribution to the power consumption will be estimated from the wake geometry. The zero-lift drag coefficient strongly depends on the local Reynolds number, *Re*_{c}=*Vc*/*ν*, with *ν* the kinematic viscosity. The state of the boundary layer, i.e. the position of turbulent transition and possibly separation, will have an effect, but this is ignored in the current work. For an approximation, we assume fully laminar attached flow and use Blasius' flat plate boundary layer solution [14,26]:
**V**|=d*s*/d*t* representing the geometric speed of the wing section, where (d*s*)^{2}=(d*x*)^{2}+(d*y*)^{2}+(d*z*)^{2}.

To complete the viscous model, a chord distribution along the wingspan is required. There is a great variety of wing shapes among birds and bats, some having rather rectangular wings, whereas others have more pointed wings. In this work, we use an elliptical chord distribution
*y*′=*y*/*b* and *c*_{max} is the maximum chord such that *S* being the wing area. It should be noted that this drag model suggests potentially beneficial effects of wing shape depending on the situation, as seen when exposing all the chord dependencies in the viscous drag equation,

The amount of thrust a bird has to produce for steady flight is, apart from body drag (*D*_{body}), depending on the drag force on the wing itself, i.e. induced drag and profile drag. These forces are evaluated locally and the components along the flight direction are summed. We compute the local induced drag as the local induced power divided by the local (geometric) velocity, so that
*T*=*D*_{ind}+*D*_{v,0}+*D*_{v,2}+*D*_{body}.

### (b) Optimization of the wake geometry

The wake geometry represents the shape of the vortex sheet as it separates from the trailing edge of the wings and convects with the freestream flow. Here, it is assumed that the induced flow does not significantly displace or deform the wake, which requires that the flight speed needs to be sufficiently large relative to the induced velocity. We assume a sinusoidal flapping motion, with equal upstroke and downstroke duration. The downstroke duration is known to vary between species [27–29] and with flight speed, but is generally close to 50% [30–35]. We leave the possible effects of varying downstroke time to future investigations. We used *N*_{U}×*N*_{b}=41×40 points (flight direction and span direction) for the nodes of a triangular mesh (each triangle representing a vortex ring), with coordinates from 0 to *U*/*f* in flight direction and −*b*/2 to *b*/2 in span direction. Then, a wing retraction transform was applied to the locations of the points. Birds and bats partially retract their wings during the upstroke [36]. In this work, wing retraction is modelled as
*b*_{max} is the maximum wingspan of the animal, *θ* the amplitude of the flapping motion and *τ* the normalized time variable running from 0 at the beginning of the downstroke, to 1 at the end of the upstroke. This transform ensures a smooth gradient from gliding flight with no upstroke wing retraction to large amplitude flapping flight with maximum wing retraction. Next, the wake mesh was transformed by rotating the wing *θ* represents the maximum excursion from the centre, rather than from the peak-to-peak amplitude often reported in studies of animal flight kinematics. A third transform tilts the stroke-plane an angle *φ* around the spanwise axis. Note that the angle *φ* is measured from the plane perpendicular to the flight direction, rather than from the horizontal plane often reported in studies of animal flight kinematics [27,29,37], and that the angle is positive when the downstroke has a forward component (figure 1).

A jackdaw (*Corvus monedula*) was chosen as a bird model, having a wingspan *b*=0.60 , a wing area *S*=0.059 ^{2} and a body mass of 0.180 [38]. The inviscid solution is characterized by reduced frequency and thrust ratio only, so the specific bird species has no influence on the optimization. For the components due to zero-lift viscous drag, the solutions are non-dimensionalized as explained in §3, thereby removing the influence of bird species. Only the lift-dependent viscous drag could be affected by the choice of species, specifically as a function of the aspect ratio, as the optimal circulation distribution (equation (2.2)) depends on the absolute chord length. However, this effect was found to be negligible (electronic supplementary material, S2, showing the sensitivity of the model to changes in *C*_{Dv,2}/*c*).

For a fixed stroke-plane angle the flapping amplitude *θ* was optimized for 200 randomized combinations of thrust ratio, reduced frequency, flight speed and lift production, over ranges specified in table 1. A thrust ratio *T*/*L*=0 corresponds to non-flapping flight, whereas *T*/*L*=0.3 is about twice the maximum expected value at cruising speed as computed using the Pennycuick model for a database of 220 birds [4], (fig. 13.11, p. 370). From that same data, one could find that birds are expected to fly with a reduced frequency of 2–3 at the minimum power speed, so that the range 1<*k*_{f}<6 covers a speed range from half the minimum power speed to twice the minimum power speed (assuming a constant flapping frequency). As the computations are performed in absolute dimensions, flight speed *U* was varied over a range plausible for the model bird, independent of reduced frequency. This forces variation in absolute flapping frequency. Finally, the lift production is varied relative to the typical bird weight, so to cover scenarios in which the bird is carrying variable loads (e.g. fat stores or prey).

The wake geometry was optimized for minimal induced power, even though one would expect an animal to minimize total power. However, at lower speeds, induced power will be dominant for most birds and bats. At high speeds, when viscous power components become dominant, the additional velocity component owing to flapping will be relatively small for the expected wingbeat frequencies. With the effect on the viscous power components being small (see also figure 8) and induced power species independent, optimizing for induced power will produce a more generally applicable model. For each optimization, an initial estimate for *θ* was used (*θ*=45°). In the second iteration, the amplitude was increased by 5°, and in the third an amplitude 5° below the initial estimate was chosen. For the next iterations, a quadratic polynomial was fitted to the three trials with the lowest power, of which the analytical minimum would serve as the next estimate. The search was limited to 0<*θ*<90°, the lower limit being non-flapping flight, whereas at the upper limit, the two wings touch at maximal excursion. The convergence tolerance was set to 1°. For each converged solution, the obtained amplitude *θ* and the resulting force components (*D*_{ind}, *D*_{v,0} and *D*_{v,2}) and power components (*P*_{ind}, *P*_{v,0} and *P*_{v,2}) were stored for analysis. The optimization was repeated for different stroke-plane angles from 0° to 50° in steps of 10°, using the same 200 combinations of thrust ratio, weight support, reduced frequency and flight speed.

## 3. Relating optimal flapping flight to non-flapping flight

The objective of this study is to characterize the effects of the reciprocating wings of flying animals, which we do by relating the optimized flapping wake data to their non-flapping equivalents (i.e. the limit *T*/*L*=0).

The optimal flapping amplitude depends on the thrust ratio *T*/*L* and reduced frequency *k*_{f}. Then, the expected dependency of the drag and power components on flapping amplitude can be substituted by a dependency on *T*/*L* and *k*_{f}. All the functions used for the nonlinear regressions are based on observation of the data. This means that the relations fit the data very well, but they are not based on any formal derivations. The obtained relations should therefore be interpreted as interpolating functions applicable to forward flapping flight within the ranges of *k*_{f} and *T*/*L* specified in table 1. The approximating functions were fitted to the data using the `nlinfit()` function from the statistics toolbox for Matlab R2013a (8.1.0.604), using bisquare robust fitting. For the drag and power components the robust fitting weights of each component were combined as

In §3a, the optimized amplitude *θ* is expressed in terms of *T*/*L*, *k*_{f} and *φ*. In §3b, the derivation of the non-dimensionalized drag components is described and the approximating functions are explained, followed by a similar analysis for the power components. §3d shortly describes the effects of stroke-plane angle.

### (a) Optimized flapping amplitude

To the best of our knowledge, there is no known analytical relation between optimal flapping amplitude, thrust ratio and reduced frequency, so a function was designed to fit the data from the optimization procedure (§2b). Flapping amplitude is limited to 90°. By using a tangent transform on the flapping amplitude, the transformed data can take any positive value. When we plot *T*/*L* and reduced frequency *k*_{f} the data points are on a single surface, which corresponds to the previously made statement that the inviscid wake is characterized by these two variables. More specifically, *T*/*Lk*_{f} collapsed all data points onto a curve, which could be approximated by the function *q*, *p*_{1/2}, *p*_{1} and *p*_{4} were assumed to be quadratic polynomials of *p*_{1/2,1}, *p*_{1/2,2}, *p*_{1,2} and *p*_{4,2}), and the model was refitted. The removal of these coefficients narrowed the confidence bounds of the remaining coefficients and reduced the overall magnitude of the residuals. The final fitting coefficients are shown in table 2. This fit approaches 95% of the data points to within the set convergence tolerance of 1° for the stroke amplitude optimization. The remaining 5% of the data points are dominated by the highest stroke-plane angles (*φ*=50° has 66% of all residuals larger than 1°, *φ*=40° has 25%) and on average deviate 2.4° from the fitted model. Figure 2 shows the fitted model as a function of reduced frequency and stroke-plane angle for different thrust ratios, indicating how increased thrust requirements require higher flapping amplitudes while increasing frequency allows for smaller amplitudes. Tilting the stroke-plane backward (increasing *φ*) requires larger stroke amplitudes, which intuitively makes sense as a compensation for the reduced vertical displacement of the wing in order to maintain thrust.

### (b) Drag components and thrust

In the wake optimization, induced drag *D*_{ind}, zero-lift viscous drag *D*_{v,0} and lift-dependent viscous drag *D*_{v,2} were computed (equations (2.6)–(2.8)). Each component has a corresponding expression for non-flapping flight (limit of *T*/*L*=0). Induced drag is expressed as

The non-flapping solution for zero-lift drag is *C*_{D,v,0} from equation (2.3)). For the wing shape described by equation (2.4) expression evaluates to approximately

The non-flapping solution for lift-dependent viscous drag can be expressed as *D*^{′}_{v,2}=2*C*_{Dv,2}*L*^{2}/*ρU*^{2}*S*. Note here the similarity to induced drag
*D*/*D*′=1+*f*_{Dv}(*k*_{f})(*T*/*L*), which results in good fits for each separate stroke-plane angle. As only *f*_{D} changes with stroke-plane angle, it can be expressed as
*p*_{0} and *p*_{1} are quadratic functions of *f*_{D}(*k*_{f},*φ*) are visualized in figure 3.

The animal flaps its wings with the purpose of generating thrust, which in steady forward flight has to counter aerodynamic drag. The total drag on an animal can be expressed as the sum of the body drag and the drag on the wings, so that the thrust required for steady flight is *T*=*D*_{body}+*D*_{ind}+*D*_{v,0}+*D*_{v,2}. With the approximations introduced above, equations (3.3), (3.5) and (3.7), all being linear functions of thrust, equation (3.9) for the required thrust ratio is found

### (c) Power components

In steady flight models, the required aerodynamic power is obtained by multiplying the drag with the flight speed, *P*=*DU*. However, to find the aerodynamic power for flapping flight, it is not sufficient to multiply the required thrust found from equation (3.9) with the flight speed, because this thrust is only the wingbeat averaged force along the direction of flight. Within the wingbeat, work is also done perpendicular to this direction. Therefore, in a similar manner as with the drag components, we non-dimensionalized the power components by their non-flapping equivalents, and examined their relation to reduced frequency and thrust ratio.

The general behaviour of the non-dimensionalized power components is very similar to their drag counterparts, so that similar models can be used
*r*, which notably improved the fit for low reduced frequencies. Again, only the factor in front of the thrust ratio is a function of stroke-plane angle and reduced frequency
*p*_{0} and *p*_{1} are again modelled as quadratic polynomials of *f*_{P} are visualized in figure 3. For all components, it seems that *f*_{D} is less than *f*_{P}, i.e. that the penalty owing to flapping is more severe on the power than on the drag.

Equation (3.9) can be used in equations (3.10), (3.11) and (3.12), which combine with *P*_{body}=*D*_{body}*U* to total aerodynamic power
*P*/*P*′) is a function of reduced frequency *k*_{f}, stroke-plane angle *φ* and thrust ratio *T*/*L*. Of these, *T*/*L* is itself a function of *k*_{f}, *φ*, and the non-flapping drag components *D*′. Note that the reciprocal of equation (3.13) is equivalent to the definition of propulsive efficiency commonly used as a performance measure for the propulsion system of fixed wing aircraft, i.e.

### (d) Optimizing stroke-plane angle

Equation (3.13) is a function of the non-flapping drag components and two kinematic parameters: stroke-plane angle *φ* and the flapping frequency in the form of reduced frequency *k*_{f}. From figure 3, it can be observed that the stroke-plane angle affects lift-dependent components different from the lift-independent components. Figure 4*a* illustrates the optimal stroke-plane angle (d*f*_{D}/d*φ*=0 and d*f*_{P}/d*φ*=0) in the case where induced drag is dominant, *D*^{′}_{ind}≫*D*^{′}_{v0},*D*^{′}_{v2}, as a function of reduced frequency. For a constant flapping frequency, reduced frequency is inversely proportional to the flight speed, so the lower the speed, the more horizontal becomes the stroke-plane.

On the other hand, when the zero-lift drag component is dominant, *b* shows only a maximum for d*f*_{D}/d*φ*=0 and d*f*_{P}/d*φ*=0 within the range of 0≤*φ*≤50°. This plot implies tilting the stroke-plane forward (*φ*<0) allows for lower values of *f*_{Dv0} and *f*_{Pv0}. Even though the non-lift producing case, which this limit implies, has not been simulated, intuitively the tendency towards a forward tilted stroke-plane makes sense, considering its similarity to breaststroke swimming (i.e. drag-based thrust).

A flying bird or bat will be between these limiting cases, so at this stage, one can only argue that the actual optimal stroke-plane angle for any species will be less than that required for minimizing induced drag (dashed line in figure 4*a*). Finding the specific optimal stroke-plane angle requires evaluating the relative contributions of each non-flapping drag component, which will differ across species.

## 4. Jackdaw example

To illustrate the implications of equations (3.3)–(3.7) and (3.10)—(3.12), again the jackdaw is used [38]. We start with determining the non-flapping drag components, the sum of which determines the base level of thrust that needs to be produced, as depicted in figure 5. These drag components are then used in equation (3.13) and (3.9). For a given flight speed and frequency, the stroke-plane angle can now be optimized. In figure 6, this is done for the jackdaw flying with a wingbeat frequency of 6.4 Hz, which is a reference frequency derived from the allometric relation *f*_{ref}=*m*^{3/8}*g*^{1/2}*b*^{−23/24}*S*^{−1/3}*ρ*^{−3/8} [4,39]. The optimal stroke-plane angle gets more horizontal for low flight speeds, which corresponds with observations of slow-flying animals [27,35]. Here, it should be noted that this horizontal orientation in this model is not because it increases the angle of attack of the wings, but purely related to increasing the horizontal component of the airflow over the wings. At high flight speeds, the optimal stroke-plane angle is approximately vertical. For speeds above 17 m s^{−1}, an actual optimum (d*η*/d*φ*=0) does not exist for positive *φ*, but considering the large separation between the contours of equal efficiency, the penalty for using any other stroke-plane angle is very small.

In §2, it was also noted that the model is only applicable to fast forward flight. In helicopter theory, the lower limit is often set to twice the induced velocity in hover [40]. In the current model, this limit is in most cases implicitly enforced by the upper limit on thrust ratio. The induced velocity is related to induced drag by *v*_{i}/*V* =*D*_{ind}/*L*. From Glauert's interpolation of induced velocity for slow flight, one can derive
*v*_{ih} is the induced velocity in hover [41]. This means that at *V*/*v*_{ih}=2 the induced drag-to-lift ratio is *D*_{ind}/*L*≈0.24. When adding profile drag it is very unlikely that the required thrust ratio will remain below 0.3 at this low speed, so that the assumption of fast forward flight is implicitly justified by the upper limit on thrust ratio.

As indicated by equation (3.9), the drag will be slightly increased owing to the flapping of the wings. For a range of flapping frequencies, figure 5 shows the extent of this increase. Focusing on frequencies above the reference frequency of 6.4 Hz, it appears the deviations from the non-flapping thrust requirement are minor. The effects are more clear when looking at the thrust efficiency, which for the reference frequency is lowest around 0.85. For higher frequencies, it seems the efficiency is almost speed independent with a value around 95%. Another interesting observation from these plots is the behaviour at low speeds. The thrust efficiency seems to increase above unity, which at first might seem wrong. However, the efficiency is taken relative to the non-flapping case. In slow flight, with a tilted stroke-plane, the beating of the wing temporarily increases the airspeed over the wing, which reduces the induced drag. If then during the recovery stroke, less lift is produced, also here there will be less induced drag, compared with the non-flapping situation.

Figure 7 shows propulsive efficiency, *η*_{P}=*T*′*U*/*P*, visualized for a range of wingbeat frequencies, together with the corresponding optimized stroke-plane angles and wingbeat amplitudes. Higher wingbeat frequencies generally increase efficiency, but the relative gain in efficiency decreases with increasing frequency. Also, the relative benefit of higher frequency decreases with increasing flight speed. Figure 7 also shows how high propulsive efficiency is related to lower wingbeat amplitudes. However, here it should be kept in mind that this concerns the optimized amplitude, i.e. reducing the amplitude without changing frequency and stroke-plane angle will inevitably result in a lower efficiency.

Considering again the reference wingbeat frequency of 6.4 Hz, the jackdaw flies with a propulsive efficiency of 80% around 10 m s^{−1} increasing to over 90% for speeds above 15 m s^{−1}. In this range, the optimal wingbeat amplitude would be around 35°. At 10 m s^{−1}, the stroke-plane angle would be tilted 20° from vertical and tilted more vertical for higher speeds. The power curve corresponding to this frequency is shown in figure 7*b*. Also the non-flapping power curve is shown, illustrating how the flapping increases the power requirement. The speed dependency of the propulsive efficiency causes the speed for minimum power to increase compared with the non-flapping curve.

For reference, the power curve as computed by Flight 1.25 (Pennycuick's model) is also shown. Below the minimum power speed, Pennycuick's model is between the flapping and the non-flapping curve. At higher speeds, it predicts much lower power requirements. A notable effect of this is the higher maximum range speed (minimal *P*/*U*). Both the non-flapping and the flapping curve predict much lower maximum range speeds. This might be related to concerns raised by field observations where birds generally fly slower than expected [7]. Field observations of jackdaw flight speeds (A. Hedenström and S. Åkesson 2014, unpublished data) indicate a preferred flight speed between 11.6 and 14.3 m s^{−1} (95% confidence interval for the corrected mean observed flight speed) for level flying individual birds (electronic supplementary material, S3). For this specific bird, the Pennycuick model predicts a maximum range speed of 15.9 m s^{−1}, which is outside of this confidence interval, whereas the model we developed in this work predicts this characteristic speed at 12.1 m s^{−1}, which is well within the observed range.

As a final consideration, the effect of flapping can be examined for each power component separately. Figure 8 shows the power factors *P*_{v,0}/*P*^{′}_{v,0}) and *k*, has been a subject of debate. Early estimates of this parameter were in the range of 1.1–1.2, commonly used for helicopters. Recently, a value as low as 0.9 has been suggested, to match the Pennycuick model with field observations [7]. As can be observed from figure 8, all of these estimates might be very optimistic, as for the jackdaw reference frequency the induced power factor never goes below 1.4. Even for the extremely high wingbeat frequency of 15 Hz, the parameter only approaches the previous worst-case estimate of 1.2. In figure 8*d*–*f,* the contribution owing to flapping is shown relative to the total non-flapping power. This shows that around the minimum power speed the reference frequency of 6.4 Hz would result in roughly 30% increase of aerodynamic power solely owing to the increased induced power, demonstrating the importance of estimating this parameter correctly. The zero-lift power factor increases strongly for decreasing speed. This is a result of the normalization with the non-flapping power component that is very small, whereas the flapping wings are experiencing a higher level of viscous drag. Looking at the relative contribution owing to flapping (figure 8*e*), it follows that the effects of flapping on the zero-lift power are relatively small. This supports the arguments in §2b for optimizing induced power only. The lift-dependent viscous power factor is very similar in behaviour to the inviscid induced power factor, though slightly lower. The relative effect on the total power curve is less than half that of the induced power factor, which makes sense considering that for this bird the non-flapping lift-dependent viscous power is approximately 57% of the induced power (equation (3.6)).

## 5. Concluding remarks

We developed a model for the aerodynamic cost of vertebrate flight that incorporates the main effects of flapping wings. For a range of combinations of thrust ratio, reduced frequency and stroke-plane angle, the wingbeat amplitude *θ* was optimized for minimum power, using a numerical vortex wake model. The results were converted into an explicit general model, represented by equation (3.13). Our model provides potential explanations to divergence between the model that is commonly used to estimate flight power in flying vertebrates [4,5] and observations [7].

The Hall and Hall model [12] computes the energetically most efficient distribution of circulation for a given wake geometry. This results in an ideal solution, equivalent to that for fixed wing lifting line theory. In all subsequent computations, it is assumed that a bird or bat wing will be able to produce this ideal wake. However, it may occur that the wake prescribes a lift coefficient that is not physically attainable for the aerofoil. Additionally, when the required lift coefficient would be realistic for the aerofoil, it may still require a geometric angle that is not within the range of motions of the wing. Figure 9*a*–*c* shows the distribution of lift coefficient along the span throughout the wingbeat at three different speeds: 5, 10 and 15 m s^{−1}, estimated for the example of a jackdaw of §4 (see electronic supplementary material, S4 for details). At 5 m s^{−1}, large lift coefficients are prescribed for the body region and the inner wing, exceeding the maximum lift coefficient for many aerofoils (though unsteady mechanisms could increase the maximum lift coefficient [42]). This suggests the bird may not be able to produce this optimal load distribution and it has to resort to alternative distributions, which inevitably has to result in increased induced power. However, the model does not take into account the additional lifting tail surface, usually maximally spread in slow flight, which would reduce the required lift coefficient at the inner wing. To achieve the prescribed lift distribution along the wing at this low speed, a considerable amount of wing twist is required, as illustrated in figure 9*d*, but bird wings are flexible and wing tip reversal is not uncommon among birds at low speeds [43]. Around the minimum power speed (10 m s^{−1}), the maximum lift coefficient is below the stall value of most aerofoils and the range of required wing twist is reduced (figure 9*b*,*e*). This trend is continued with increasing speed (figure 9*c*,*f*).

On a similar note, the effect of wake roll-up has not been taken into account. Even before the wake leaves the wing, it starts to roll-up into the tip vortices, which could alter the realized lift distribution. Fixed wing aircraft frequently use wing tip devices, which are designed to improve the efficiency of the wing [44]. The separated outer primaries of bird wings have been suggested to serve a similar purpose [45–47]. In fixed wing aerodynamics, these deviations are dealt with by introducing a span efficiency factor. A similar factor should be expected for flapping flight, though its value could be much more variable compared with its fixed wing equivalent. A consideration in future models could thus be to quantify these effects in real animals and include them in the model.

## Ethics statement

This study did not involve any experiments on human or animal subjects.

## Data accessibility

All computations were performed in Matlab R2013a (8.1.0.604) (The MathWorks, Inc). The output dataset that was used for generating the explicit functions is uploaded as electronic supplementary material.

## Funding statement

The research was supported by the Swedish Research Council (A.H. 621-2009-4965/621-2012-3585 and L.C.J. 621-2013-4596) and from The Royal Physiographic Society in Lund (M.K.). This work received support from the Centre for Animal Movement Research (CAnMove) financed by a Linnaeus grant (349-2007-8690) from the Swedish Research Council and Lund University.

## Author contributions

M.K. conceived of the study, performed the computations and analysis, and drafted the manuscript.. L.C.J and A.H. contributed to the design of the study, interpreting the data and drafting the manuscript. All authors gave final approval for publication.

## Conflict of interests

We have no competing interests.

## Acknowledgements

We are grateful to Colin J. Pennycuick for development of software for obtaining and analysing bird flight speed data and to Susanne Åkesson for her help with measuring jackdaw flight speeds. We thank the referees for their constructive comments on the manuscript.

- Received December 9, 2014.
- Accepted March 10, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.