Forward flight of birds revisited. Part 2: short-term dynamic stability and trim

Thrust generation by flapping is accompanied by alternating pitching moment. On the down-stroke, it pitches the bird down when the wings are above its centre of gravity and up when they are below; on the up-stroke, the directions reverse. Because the thrust depends not only on the flapping characteristics but also on the angle of attack of the bird's body, interaction between the flapping and body motions may incite a resonance that is similar to the one that causes the swinging of a swing. In fact, it is shown that the equation governing the motion of the bird's body in flapping flight resembles the equation governing the motion of a pendulum with periodically changing length. Large flapping amplitude, low flapping frequency, and excessive tilt of the flapping plane may incite the resonance; coordinated fore–aft motion, that uses the lift to cancel out the moment generated by the thrust, suppresses it. It is probably incited by the tumbler pigeon in its remarkable display of aerobatics. The fore–aft motion that cancels the pitching moment makes the wing tip draw a figure of eight relative to the bird's body when the wings are un-swept, and a ring when the wings are swept back and fold during the upstroke.


Summary
Thrust generation by flapping is accompanied by alternating pitching moment. On the down-stroke, it pitches the bird down when the wings are above its centre of gravity and up when they are below; on the up-stroke, the directions reverse. Because the thrust depends not only on the flapping characteristics but also on the angle of attack of the bird's body, interaction between the flapping and body motions may incite a resonance that is similar to the one that causes the swinging of a swing. In fact, it is shown that the equation governing the motion of the bird's body in flapping flight resembles the equation governing the motion of a pendulum with periodically changing length. Large flapping amplitude, low flapping frequency, and excessive tilt of the flapping plane may incite the resonance; coordinated foreaft motion, that uses the lift to cancel out the moment generated by the thrust, suppresses it. It is probably incited by the tumbler pigeon in its remarkable display of aerobatics. The fore-aft motion that cancels the pitching moment makes the wing tip draw a figure of eight relative to the bird's body when the wings are un-swept, and a ring when the wings are swept back and fold during the upstroke.          (positive for nose up); vertical translation of the bird's body (and hence of the twist axis) is h. Following equations (3.1) and (5.10) of Part 1, it is assumed that the twist varies linearly along the span; moreover, it follows the flapping rateφ with where the over-dot marks a derivative with respect to time, y ∈ (−1, 1) is the span-wise coordinate, ε ∈ (0, 1) is a certain proportionality coefficient, and α g0 is the wing twist at the shoulder that remains constant throughout the flapping cycle. The flapping-rate-proportional twist, embodied in the second term on the right, was shown in Part 1 to play a central role in making the flapping propulsion aerodynamically efficient.
In order to comply with the restrictions of the aerodynamic theory developed in Part 1, it is assumed that the wing has an elliptical plan-form, with the chord length prescribed by It is also assumed that A −1 , α g0 , τ ,ḣ, φ,φ, λ andλ are small when compared with unity. It was shown in Part 1 that under this assumption, α g0 and τ are equivalent. Accordingly, without a loss of generality, it can be agreed that when not flapping, the bird is trimmed with τ = 0, and the necessary angle of attack is furnished by α g0 .

Equations of motion
Intuitively, the longitudinal dynamics of a bird in flapping flight should be characterized by three distinctive time scales, representing the short-period mode, the phugoid mode, and the flapping itself. Short-period mode is manifested in pitch and heave oscillations with practically no changes in airspeed; phugoid mode is manifested in pitch, heave and airspeed oscillations together ( [1] and [5, pp. 167-169]). The longest of the three is the time scale of the phugoid mode, a few seconds (if this mode is stable, its period is estimated as π √ 2v/g). The shortest of the three is probably the flapping period, a few tenths of a second (appendix D). With an order of magnitude separating the shortest and the longest time scales, one may expect that the interaction between the two is small, and hence the flight velocity can be practically assumed constant over the flapping period [1].
With constant flight velocity and nominally horizontal flight, the longitudinal dynamics is governed by the pair of equations and where M is the pitching moment about the centre of gravity, L is the lift and the over-bar denotes the respective quantity in the adjoint non-flapping flight ( §5.1 in Part 1). In the context of this paper, it is formally defined as the flight at the same velocity and with the same α g0 , but with This definition is consistent with the definition (5.5) of Part 1 as long as the flapping amplitude is sufficiently small. It is implicitly assumed that in the adjoint flight,L exactly offsets weight: and that the bird trims out with sweep angleλ: Based on the aerodynamic model of Part 1, it is shown in appendix A that L and M can be written in the following symbolic forms: where L ,α , L ,φ , M ,λα , M ,λφ , M ,φα 2 , M ,αφφ , M ,φφ 2 , M ,Dφ ,M ,α , M ,τ andM ,φ are the respective partial derivatives (with respect to the variable(s) following the comma), andD w is the drag of the wing in the adjoint flight; the derivatives that depend onλ have been marked by the over-bars. Explicit expressions for these derivatives, which are based on the model of Part 1, can be found in appendix A. In principle, however, most (if not all) of these derivative could have been extracted from experimental measurements, as in Taylor & Zbikowski [2] or from numerical simulations, as in Wu & Sun [6]. Substituting (4.7) in (4.1) and (4.6) in (4.2) yields the equations of motion in explicit form: In this particular form, they will be used in §8; two additional forms will be needed for § §5-7. To derive the first, we will need the reduced heave velocity, to derive the second, we will need the angle of attack of the bird's body, of (4.9), and its combination with (4.8), of (4.11) with (4.9), and the combination of its corollary, (which follows by (4.11)) with (4.8), of partial derivatives appearing in the second and the third terms on the left-hand sides of (4.13) and (4.16) (they appear multiplied by I y ) will be identified with the natural frequency and the damping of the short-period mode [5, p. 175] in the adjoint non-flapping flight. We assume that bothω 2 s and g s are real and positive-if they were not, the short-period mode would have been unstable, and hence incompatible with the assumptions underlying this short-term stability analysis. There are many arguments that can be brought in support of this statement. Suffice it to say that obtaining an unstable short-period mode contradicts the initial assumption made on its characteristic time scale, and variations in airspeed can no longer be ignored.

Flapping resonance
We would like to demonstrate that without a coordinated fore-aft sweeping motion of the wings, flapping flight could be dynamically unstable. To this end, we set λ =λ (5.1) throughout the flapping cycle and assume that the response in heave is sufficiently small to justify the neglect of the two nonlinear terms on the left-hand side of (4.13). The solution of this linearized version of (4.13) is a combination of homogeneous and particular solutions. Since the right-hand side of this equation does not vanish at all times during the flapping period, if we could demonstrate that its homogeneous part admits a diverging solution, the initial proposition would have been established. The details follow. For the sake of simplicity, we assume that the flapping is harmonic with amplitude φ 0 and (reduced) angular frequency ω: and together with (4.17), (4.18), (5.1)-(5.6) in the homogeneous part of (4.13), allows bringing it into the form: the nonlinear terms that have been neglected in due course are ψ 2 1 (ω/α g0 )(dη/dt) 2 sint and δ 2 1 (ω/α g0 )(dη 2 /dt) sint. Equation ( where b's depend on the initial conditions and μ is Floquet's exponent. b's satisfy the recurrence relation μ makes the determinant of the infinite, 5-diagonal matrix, formed on the coefficients of b's in (5.9), vanish: Equation (5.9) is obtained by substituting (5.8) in (5.7); without (5.10), the relations in (5.9) would have admitted only a trivial solution.
By mapping the solutions of (5.10) over possible combinations of γ s , κ s , ψ 1 , ψ 2 , δ 1 and δ 2 one can identify those combinations for which Im μ > −γ s ; that is, one can identify the conditions for which the solutions of (5.7) are bounded. Details can be found in appendix B; the outcome is shown in figure 1. Judging by this figure, a rational description of those conditions is hardly possible. Nonetheless, a few  conclusions can be made based on explicit form of the parameters in (5.7): and Here δ 2 /ψ 2 and δ 1 /ψ 1 are only weakly dependent on flight conditions (through the square root of the flapping frequency). For most species of the birds compiled in appendix D, δ 1 /ψ 1 ∈ (0.5, 1) at cruise; moreover, assuming ε = 0.5, δ 2 /ψ 2 ∈ (0.1, 0.2). Three combinations of δ 1 /ψ 1 and δ 2 /ψ 2 from these ranges are shown in the right three columns of figure 1.
M ,τ is contributed mainly by the tail. Since the tail is closed at cruise, one may expect this derivative to be small. Consequently, the value of γ s is determined largely by L ,α /2mω. For most of the species of birds compiled in appendix D, L ,α /2mω ∈ (0.05, 0.3). Shaded areas in figure 1 reflect values of γ s in this range.
The value of κ s depends on the stability margin and the flapping frequency. The stability margin of a bird is probably small-otherwise large pitching up moment would have been required to trim the bird in the adjoint flight [8]. Small stability margin implies small ω s , and hence one may expect that the flapping frequency will be higher than the frequency of the short-period mode in the adjoint flight. Three cases with κ s < 1 are shown in the first three rows in figure 1.
Based on the first three rows in the right three columns of figure 1, one can state that with the exception of a few singular cases, a bird is unstable if where Ψ 1 and Ψ 2 are, in general, intricate functions of κ s , γ s , δ 1 /ψ 1 and δ 2 /ψ 2 . Towards the following discussion, they will be replaced by their representative values over these nine figures, say, Ψ 1 = 1 and Ψ 2 = √ 0.1. Using (5.12), the first criterion in (5.16) translates into where f * = ωv/2π s is the dimensional flapping frequency. Whichever criterion applies, a combination of small radius of gyration, slow flapping and large flapping amplitude may lead to resonance. The data on the radius of gyration among birds is scarce. Exploiting the interpretation of the radius of gyration as half the length of a dumbbell having the same mass and inertia as the body represented by it, we estimate the radius of gyration of a bird to be comparable with half the chord at the wing's root, c 0 = 8/π A-see (3.3). In other words, when writing and

Inciting the resonance
We release now the assumption (5.1) that inhibited the wing sweep, and temporarily assume that the sweeping motion is prescribed by where k ψ 2 and k ψ 1 are certain parameters. This assumption will be released in the next section. Combined with flapping, the fore-aft motion prescribed by (6.1) makes the wing tip draw an oblique figure of eight (the first term makes it a figure of eight; the second term tilts the flapping plane forwards). Replacing (5.1) with (6.1), and repeating the steps leading from (4.13) to (5.7), yields a similar equation, only in which and Although it is plausible that some combination of k ψ 2 and k ψ 1 will make the bird stable, it is certain that sufficiently large k ψ 2 or k ψ 1 will make it unstable. This conjecture follows by observing the stability regions in the left column of figure 1 (as k ψ 2 and k ψ 1 increase, δ 1 /ψ 1 and δ 2 /ψ 2 tend to zero). It is plausible that the sweeping-motion-incited resonance is exploited by the tumbler pigeon (Columba livia) in its repetitive somersaulting. Returning to (5.7) one may note that substitutiont =t + π changes sign with the terms involving ψ 1 and δ 1 . Consequently, the stability analysis of the preceding section applies for positive and negative values of ψ 2 1 alike (only its absolute value counts), and hence the resonance can be incited with large negative values of k ψ 1 as well as with large positive values. In other words, it can be incited by tilting the flapping plane forwards and backwards alike.
All birds tilt the flapping plane backwards during the transition from forward to hovering flights [4]. In view of the above, it can incite the resonance. Since the tilt of the flapping plane is dictated by the performance requirements-the thrust is needed for lift-it is plausible that during these stages of flight the tail replaces the wing as the primary active control. In fact, the tail has been observed to open up with decreasing flight speed [4].

Trim at zero angle of attack
In general, there is infinite number of flight strategies that can be realized using active control. Here, we consider the most obvious two: keeping the angle of attack zero throughout the flapping cycle, and keeping the pitch angle zero throughout the cycle. Starting with the first strategy, the sweeping motion that makes α = 0 (7.1) is the one that makes the right-hand side of (4.16) vanish: by (4.14) and (7.1). By interpretation, the numerator on the right-hand side of (7.2)-if written as a single fraction-is the pitching moment at α = 0 and λ =λ, unmitigated by the fore-aft motion of the wing. To within a sign, the denominator in (7.2) is twice the flapping moment acting on the right wing at α = 0, (M x ) α=0 ; this conjecture follows from (A 7) and (A 8) of appendix A, and (4.26) of Part 1. In turn, can be associated with the action of lift, (L) α=0 , at the centre of pressure of the right wing, The expression on the left of (7.5) is just a definition of the centre of pressure (twice the moment generated by a single wing over twice its lift); the expression on the right follows from (4.4), (4.5), (4.24), (4.26) and (4.31) of Part 1. This interpretation of (7.2) elucidates the balancing action: the lift is moved fore-and-aft to counteract the moment created by the thrust and by the drag. A few examples of the wing tip trajectories can be found in figure 3. Individual contributions of the six terms on the right-hand side of (7.2) can be found in figure 4. The fundamental shape drawn by the wing tip is a deformed figure of eight. The three terms that are responsible for this shape are the first and the last two terms on the right-hand side of (7.2), those involving the combinationsḣφ, φφ and φφ 2 . All three can be associated with the pitching moment created by the thrust. Loosely speaking, thrust is generated by the rotation of the lift vector-it rotates forwards when the wing moves down (either due to flapping or due to heave) and backwards when the wing moves up. Since heave is in phase with the flapping angle (the bird moves up when the wings are down-see (7.3)), the lift rotation due to heave pitches the bird down during the entire flapping cycle. To counteract it, the wing tips draw a crescent shape, moving forward towards the ends of the upstroke and the down-stroke (figure 4a). At the same time, the lift rotation due to flapping creates an alternating moment, pitching the bird down at the beginnings of the up-stroke and the down-stroke, and pitching it up towards their respective ends. To counteract it, the wing tips draw a figure of eight, moving backwards during both the up-stroke and the down-stroke (figure 4g).
Modifications to the basic figure of eight come from the remaining three terms. The first one involves φ. This term can be associated with the pitching moment needed to coordinate the pitch with heave in keeping the angle of attack constant. To this end, the pitch angle should be proportional to the heave velocity, and hence to the flapping angle. In turn, it requires the pitching moment to be proportional to the pitch acceleration, and hence to the flapping angle acceleration. It brings the wings forwards (relative to their nominal positions) at the beginning of the down-stroke, and backwards at its end; the down-stroke and the up-stroke are different because the lift is different ( figure 4b).
The    Conditions are those of case 11 in appendix C. During the upstroke (π/2 < ωt < 3π/2), the unmitigated pitching moment, represented by N, is small, but the lift generated by the wing is insufficient to balance it.
centre of gravity (assuming M bt ,α small, it also makes the bird near neutrally stable-see (A 26)). When the bird starts flapping, the wing's centre of pressure moves outwards (and hence backwards) on the downstroke and inwards (and hence forwards) on the up-stroke-see (7.5). Consequently, the sweep has to be smaller on the way down than on the way up. Moreover, if the wings fold on the way up, the centre of pressure moves further inwards (and hence forwards), magnifying the effect. In fact, it can make the wing tip draw a non-intersecting ring rather than a figure of eight [4].
When the angle of attack due to flapping is comparable with the average angle of attack, α g0 -either because of high flapping rate or because of small twist or because of small α g0 -the denominator in (7.2) may vanish at certain times during the up-stroke (figures 3g,h and 5). During these events, it will be impossible to balance the pitching moment with fore-aft adjustment of the wing. A bird has two options here. One is to use the tail for control; in fact, the tail often opens up during vigorous flapping. The other is to do nothing. Since in those cases where the denominator in (7.2) can vanish, the pitching moment is small during the entire upstroke (note the range π/2 < ωt < 3π/2 in figure 5), a momentary imbalance should be of no consequence.

Trimming at zero pitch angle
in (4.8) yields the sweeping angle required to keep the pitch attitude during flight: To obtain (8.2) in this particular form, we have used the identity  A few examples of wing tip trajectories can be found in figure 6. They look similar to those shown in figure 3, but since there is no need to adjust the pitch angle to heave, less fore-aft motion on the down-stroke is required.

Concluding remarks
Without active control, an interaction between flapping, pitching and heave can cause a resonance that is similar to the one that causes the swinging of a swing. It may affect all birds that have small inertia in pitch when flapping slowly and with large amplitude, irrespective of their stability margin. This resonance can be suppressed by active control-either with fore-aft sweeping motion of the wing or with up-and-down deflection of the tail. It is most probably incited by a tumbler pigeon in its repetitive somersaulting. It is possibly incited by all birds during the transition from forward to hovering flight, and suppressed by the tail (in this flight regime most birds open their tails). It is suppressed in forward flight (where the tail is closed) by the sweeping motion of the wing. Our obtaining wing tip trajectories, that resemble those observed in birds, by simply enslaving the sweeping motion of the wing to keep either the angle of attack or the pitch attitude supports this conjecture.
Three key elements made this analysis possible. One is the aerodynamic model that was developed in Part 1. It allowed introduction of the aerodynamic derivatives (appendix A), which were instrumental in keeping the length of the associated equations in check. The second element was the theory of the third-order differential equation governing the short-term dynamics of the bird that was developed in appendix B. It allowed obtaining the stability boundaries without running time-consuming numerical solutions of this equation. The third element was in restricting the range of the six parameters of that equation based on physical data (appendix D). It would have been practically impossible to draw any definitive conclusions with six unrestricted parameters.

Appendix A. Lift and pitching moment
In this appendix, we derive explicit expressions for the lift and the pitching moment acting on a bird in flapping flight. The main contributors to both parameters are the wing(s), the body and the tail.
The combined body and tail contribution to the lift is expected to be small. In fact, exploiting the slender body theory, it should be proportional to the ratio of the tail span squared and the wing(s) area. Neglecting this contribution a priori, L is furnished by the conjunction of (4.4), (4.24) and (5.11) of Part 1 (the last equation applies because of (3.2)). We write it here in the symbolic form: In (A 1), as well as in all subsequent occurrences, a comma in the subscript stands for a partial derivative with respect to the variable(s) following it. We divide the pitching moment about the centre of gravity, into four constituents: The last one comes from the body and the tail; the first three come from the wing. Among these three, M w is the pitching moment about the quarter-chord point of the wing's root, contributed by the pressure loads acting on the wing-it is furnished by (4.27) in Part 1 (in which it was denoted M y ); Lx cg shifts this moment to the centre of gravity; M d is the pitching moment contributed by the parasite drag. Based on (4.27) of Part 1, and similar to (A 1), the expression for M w is recast as where the respective coefficients (partial derivatives), and are identified by matching between (A 6) and the original expression in Part 1. I 11 and I 13 are actually standard integrals, the former equals −1/3 and the latter equals −1/5. K 1 , K 4 and K 5 are functions of the aspect ratio found in (4.28) and (4.29) in Part 1. They are approximated with where k 1 ≈ 1.29, k 4 ≈ 2.29 and k 5 ≈ 6.12-see The aerodynamic model of Part 1 does not provide the contributions of the wing to the rate derivatives, L ,τ and M ,τ . For high-aspect ratio (almost) un-swept wings, both contributions are invariably small [5, p. 137]. In dynamic stability analysis, L ,τ is customarily neglected altogether. Towards the analysis carried out in this paper, the missing contribution to M ,τ can be easily compensated for by increasing respective contribution of the tail (see below)-it will not be needed, however.
Assuming that all wing sections have the same parasite drag coefficient, D w 0 , the wing's drag contribution to the pitching moment, M d , follows (3.3) by quadrature. For future use, it is put into the form: where M ,τ is the damping in pitch contributed mainly by the tail, M bt ,α is the stiffness in pitch contributed by the tail and the body alike, M bt 0 is the static moment contributed mainly by the tail, and M bt ,φ reflects the reaction of the tail to the periodic downwash of the wing. None of these contributions will be explicitly required in the text; in fact, since the tail is normally closed in forward flight, one may expect M bt -and, consequently, all its partial derivatives-to be small for all but long-and-wide-bodied birds.
Combining (A 1), (A 5), (A 6), (A 14) and (A 15) together, yields where, by interpretation, α = τ −ḣ is the angle of attack of the bird's body. Substituting  The following four combinations of partial derivatives: naturally arise when grouping similar terms in (4.13); in view of (A 21), they can be interpreted as partial derivatives at constant L. Additional combination, arise when groupingφ terms in (4.16); in view of (4.15), it can be interpreted as the respective derivative at constant α.

Appendix B. Solutions of (5.7)
Consider the third-order differential equation with two time-dependent coefficients: where b's depend on the initial conditions and μ is Floquet's exponent [7, p. 414]. The factor i (μ + iγ s + n), introduced in (B 2), is not a necessity, but it does make the following equations somewhat simpler. b's satisfy the recurrence relations whereκ 2 s = κ 2 s − γ 2 s ; they follow by substitution of (B 2) in (B 1). These relations admit a non-trivial solution only if the determinant H of the infinite 5-diagonal matrix, formed on the coefficients of b's in (5.9), vanishes. We will write this condition as It determines μ. We seek an analytical solution of (B 5). To this end, we form an auxiliary matrix, Consequently, where B 1 = 1 + 2 cot(πκ s )C 1 and B 2 = 2 cot(πκ s ) cot(iπγ s )C 2 replace C 1 and C 2 as the pair of independent parameters. In general, B 1 and B 2 can be established by fitting (B 11) to H at two points, say, μ = 0 and μ = μ 1 :  Exploiting standard trigonometric identities, it can be rewritten as The solution was examined att = 400 and classified as either diverging or converging. All solutions behaved as predicted.

Appendix D. Physical data
The data in this appendix were collected from the earlier studies ( [4] and [11][12][13][14]); they are reproduced in the first six columns of table 3. In particular, m * is the mass of a bird in kilograms, 2s and 2S are the span and the wing area in metres and square metres, respectively, A is the aspect ratio, f * is the wing-beat frequency at cruise in cycles per second and v is the observed velocity at cruise in metres per second. u * is the velocity at which the drag in the adjoint non-flapping flight is minimal; it was estimated with u * = (m * g) 1/2 (ρS) −1/2 (π AD 0 ) −1/4 , where D 0 is the parasite drag coefficient [15]. The numbers shown in the table were based on D 0 = 0.015 and standard density; one can verify that u * provides a fair approximation for v. u * replaced v when the latter was not included in the original data. In the following columns, ω = 2π f * s/v is the reduced frequency, m = m * /ρSs is the reduced mass; L ,α has been computed with (A 2). The next four columns reflect (5.21), (5.22), (5.14) and (5.15), respectively; ψ 2 and δ 2 /ψ 2 have been estimated with ε = 0.5. The choice of 0.015 for D 0 was somewhat arbitrary, as the actual value of the drag coefficient depends on the shape of the bird's body and its size relative to the wing; moreover, D 0 changes with the Reynolds number. Preliminary design tools found in reference [16] suggest that a tailless airplane, featuring an elliptical wing with aspect ratio of eight, and a double-ogive body, two root-chords long and a half rootchord in diameter, should have D 0 between 0.015 and 0.018 at chord-based Reynolds numbers between      200 000 and 80 000. This estimate is based on the assumption that the boundary layer on its surface is turbulent. It drops by half if the boundary layer can be assumed laminar.