In this section, we introduce flapping flight dynamics and describe the bird model used in our computational framework. Furthermore, we describe how such a dynamical model is used in order to identify steady and level flapping flight regimes, study their stability, and assess their energetic performance.
Equations of motion modelling flapping flight
Flight dynamics is restricted to the longitudinal plane and thus the bird main body is captured as a rigid-body with three degrees of freedom, i.e. two in translation and one in rotation. This model preserves symmetry with respect to this plane, without any lateral force and moment. The aerodynamic model of the wing relies on the theory of quasi-steady lifting line23. Additionally, the present work does not account for the inertial forces due to the acceleration of the wing, and thus also neglecting the so-called inertial power. This inertial power was shown to be negligible in fast forward flight conditions, in comparison to the other contributions to actuation power24, and is thus systematically neglected in similar work10,11,25 since wing inertia is neglected.
The body is thus modelled with a mass \(m_b\) and a rotational inertia \(I_{yy}\) about its center of mass. The equations of motion are expressed in the body frame \(G(x’, y’, z’)\) with unit vectors \((\hat{\textbf{e}}_{x’}, \hat{\textbf{e}}_{y’}, \hat{\textbf{e}}_{z’})\), and an origin located at the center of mass, as pictured in Fig. 1a. The state space vector is thus
$$\begin{aligned} \textbf{x} = \{u, w, q, \theta \} \end{aligned}$$
where u and w are the body velocities along the \(x’-\) and \(z’-\)axis and \(\theta\) and q are the pitch angle and its time derivative about the \(y’-\)axis, respectively. Consequently, the equations of motion read11,13,26
$$\begin{aligned} \begin{aligned} \dot{u}&= -qw – g\sin \theta + \frac{1}{m_b}\big ( {F_{x’}(\textbf{x}(t), t)} \\&\quad + {F_{x’, t}(\textbf{x}(t), t)} + D (\textbf{x}(t), t) \big )\\ \dot{w}&= qu + g\cos \theta +\frac{1}{m_b} \big ( F_{z’}(\textbf{x}(t), t) + F_{z’, t}(\textbf{x}(t), t) \big ) \\ \dot{q}&=\frac{1}{I_{yy}} \big ( M_{y’}(\textbf{x}(t), t) + M_{y’, t}(\textbf{x}(t), t) \big )\\ \dot{\theta }&= q \end{aligned} \end{aligned}$$
(1)
(a) Bird model for describing the flight dynamics in the longitudinal plane. The state variables are expressed with respect to the moving body-frame located at the flier’s center of mass \(G(x’,z’)\). These state variables are the component of forward flight velocity, u, the velocity component of local vertical velocity, w, the orientation of this body-centered moving frame with respect to the fixed frame, \(\theta\) and its angular velocity, q. A second frame \(O(x’_{w}, z’_{w})\) is used to compute the position of the wing, relative to the body. The wings (dark gray) and the tail (red) are the surfaces of application of aerodynamic forces. (b) Top view of the bird model. The left wing emphasizes a cartoon model of the skeleton. The shoulder joint s connects the wing to the body via three rotational degrees of freedom (RDoF), the elbow joint e connects the arm with the forearm via one RDoF and the wrist joint w connects the forearm to the hand via two RDoF. Each feather is attached to a bone via two additional RDoF, except the most distal one ”1” which is rigidly aligned with the hand. The right wing further emphasizes the lifting line (red) which is computed as a function of the wing morphing. The aerodynamic forces generated on the wing are computed on the discretized elements \(P_{i}\). The tail is modelled as a triangular shape with fixed chord \(c_{t}\) and maximum width \(b_{t}\) that can be morphed as a function of its opening angle \(\beta\). (c) Wing element i between two wing profiles, identifying a plane \(\Sigma\) containing the lifting line (red). (d) Cross section of the wing element, containing the chord point \(\mathbf {P_i}\) where the velocities are computed (Color figure online).
The forcing terms in Eq. (1) are the aerodynamic forces and moments applied to the wing (namely \(F_{x’}\), \(F_{z’}\), and \(M_{y’}\) ) and to the tail (\(F_{x’, t}\), \(F_{z’, t}\), and \(M_{y’, t}\)). The whole drag is captured by an extra force D that sums contributions due to the body \(D_{b}\), the skin friction of the wing (wing profile) \(D_{p,w}\), and the skin friction of the tail (tail profile) \(D_{p,t}\). These terms are described in detail in the next sections. Importantly, we accounted for the drag acting purely along \(x’\) direction, after proving that the projection of the drag forces along \(z’\)-axis is between two and three orders of magnitude smaller with respect to the aerodynamic forces produced by two other main lifting surfaces. This assumptions holds for the fast forward flight regime that are subject of our study, but such components of drag along \(z’\) axis should be accounted for other flight situations.
Wing model
The bird has two wings. Each wing is a rigid poly-articulated body, comprising the bird arm, forearm and hand, as pictured in Fig. 1b. Each segment is actuated by a joint to induce wing morphing. We refer to13,15 for a complete description of this wing kinematic model.
Each joint is kinematically driven to follow a sinusoidal trajectory specified as:
$$\begin{aligned} q_{i}(t) = q_{0,i}(t) + A_{i} \sin (\omega t + \phi _{0,i}) \end{aligned}$$
(2)
with \(\omega = 2 \pi f\) and f being the flapping frequency which is identical for each joint, \(q_{0,i}\) being the mean angle over a period (or offset), \(A_i\) the amplitude, and \(\phi _{0,i}\) the relative phase of joint i. A complete wingbeat cycle is therefore described through a set of 19 kinematic parameters, including the frequency f.
We assume that the wing trajectory is rigidly constrained, and therefore we do not need to explicitly solve the wing dynamics. Under this assumption, the motion generation does not require the computation of joint torques. The model further embeds seven feathers of length \(l_{ki}\) in each wing. The feathers in the model have to be considered a representative sample of the real wing feathers. They thus have a limited biological relevance; their number is chosen so as to interpolate the planform satisfactorily and to smoothly capture the morphing generated by the bone movements. These feathers are attached to their respective wing bones via two rotational degrees of freedom allowing them to pitch and spread in the spanwise direction. These two degrees of freedom are again kinematically driven by relationships that depend on the angle between the wing segments13. This makes the feathers spreading and folding smoothly through the wingbeat cycle. In sum, the kinematic model of the wing yields the position of its bones and feathers at every time step. This provides a certain wing morphing from which the wing envelope (leading edge and trailing edge) can be computed (see Fig. 1b). From the wing envelope, the aerodynamic chord and the lifting line are computed. The lifting line is the line passing through the quarter of chord, which is itself defined as the segment connecting the leading edge to the trailing edge and orthogonal to the lifting line (Fig. 1b). This extraction algorithm is explained in detail in15.
In order to calculate the aerodynamic forces, the angle of attack of the wing profile has to be evaluated. Each wing element defines a plane containing the lifting line and the aerodynamic chord as pictured in Fig. 1c. The orientation of the plane is identified by the orthogonal unit vectors \((\hat{\textbf{e}}_n, \hat{\textbf{e}}_t, \hat{\textbf{e}}_b)\), where \(\hat{\textbf{e}}_n\) is the vector perpendicular to the plane and \(\hat{\textbf{e}}_t\) is the tangent to the lifting line. To compute the effective angle of attack, the velocity perceived by the wing profile is computed as the sum of the velocities due to the body and wing motion, and the velocity induced by the wake. The first contribution, \(\textbf{U}\), accounts for
$$\begin{aligned} \textbf{U} = \textbf{U}_{\infty } – \textbf{v}_{kin} – \textbf{v}_{q}\end{aligned}$$
where \(\textbf{U}_{\infty } = u \hat{\textbf{e}}_{x’} + w \hat{\textbf{e}}_{z’}\) is the actual flight velocity, \(\textbf{v}_{kin}\) is the relative velocity of the wing due to its motion, and \(\textbf{v}_{q}\) is the component induced by the angular velocity of the body q and calculated as
$$\begin{aligned} \textbf{v}_{q} = q\hat{\textbf{e}}_{y’} \wedge (\textbf{P}_{i} – \textbf{G})\end{aligned}$$
This velocity vector \(\textbf{U}\) defines the angle \(\alpha\), as pictured in Fig. 1d.
The second contribution is due to the induced velocity field by the wake, i.e. the downwash velocity \(w_{d}\), and acting along the normal unit vector \(-w_{d}\hat{\textbf{e}}_n\). The resulting effective angle of attack, \(\alpha _{r}\), is thus
$$\begin{aligned} \alpha _{r} = \alpha – \frac{w_{d}}{|\textbf{U}|}\end{aligned}$$
The downwash velocity \(w_d\) is computed according to the Biot-Savart law23, assuming the wake being shed backwards in the form of straight and infinitely long vortex filaments at each time step of the simulation13,15. This quasi-steady approximation is justified a posteriori by ensuring that our reduced frequency, inversely proportional to the unknown airspeed, never exceeds the value of 0.2, below which the effects of time-dependent wake shapes on wing circulation are negligible (e.g. see discussion in27). Once the downwash is evaluated, it is possible to evaluate the circulation, and consequently the aerodynamic force and moment acting at the element \(P_i\), i.e. \(F_{x’, i}(\textbf{x}(t), t), F_{z’, i}(\textbf{x}(t), t), M_{y’, i}(\textbf{x}(t), t)\), as explained in detail in13. We use the thin airfoil theory for the estimation of the lift coefficient, with a slope of \(2\pi\) that saturates at an effective angle of attack \(\alpha _{r}\) of \(\pm 15^{\circ }\).
Drag production by body and wing
The main body and the wings induce drag that should be accounted for in a model aiming at characterizing energetic performance. Body-induced drag is named parasitic because the body itself does not contribute to lift generation, and only induces skin friction and pressure drag around its envelope28. The total body drag is
$$\begin{aligned} D_{b} = \frac{1}{2}\rho C_{d, b} S_{b}|\textbf{U}_{\infty }|^{2} \end{aligned}$$
(3)
where \(\rho\) is the air density. We implemented the model described by Maybury28 to compute the body drag coefficient \(C_{d, b}\). This depends on the morphology of the bird and the Reynolds number Re according to
$$\begin{aligned} C_{d,b} = 66.6m_{b}^{-0.511}FR_{t}^{0.9015}S_{b}^{1.063}Re^{-0.197} \end{aligned}$$
(4)
with \(S_{b}\) and \(FR_{t}\) are respectively the frontal area of the body and the fitness ratio of the bird, and both of them can be estimated from other allometric formulas i.e.28,29.
$$\begin{aligned} S_{b}= & {} 0.00813m_{b}^{2/3} \end{aligned}$$
(5)
$$\begin{aligned} FR_{t}= & {} 6.0799m_{b}^{0.1523} \end{aligned}$$
(6)
The Reynolds number \(Re = \rho |\textbf{U}_{\infty }| \overline{c} / \mu\) is calculated with the reference length of the mean aerodynamic chord \(\overline{c}\), with \(\mu\) being the dynamic viscosity. This model is found to be suitable for Reynolds number in the range \(10^{4}-10^{5}\)28. Another source of drag is the profile drag due to friction between the air and the feathers on the wings. It is the sum of the profile drag at each section along the wingspan, i.e.
$$\begin{aligned} D_{p,w} = \frac{1}{2} \rho C_{d, pro} \sum _{i=1}^{n} c_{i}|\textbf{U}_{r,i}|^{2} ds_{i} \end{aligned}$$
(7)
with \(c_{i}\) the chord length, \(ds_{i}\) the length of the lifting line element along the wingspan, and \(\textbf{U}_{r,i}\) the velocity at the wing section i accounting for the body velocity, the kinematics velocity of the wing and the downwash velocity (Fig. 1c,d). We used a value of profile drag of \(C_{d, pro} = 0.02\) and this is assumed to be constant over the wingspan and throughout the flapping cycle30. In reality, due to the wing motion, this value should be gait dependent. However, the aforementioned assumption has been largely used in previous works31,32.
Tail model
Since the span of the tail is of the same magnitude as its aerodynamic chord, here the lifting line approach cannot be used23. Therefore, the tail is modelled according to the slender delta wing theory, as a triangular planform33. Its morphology is illustrated in Fig. 1b and characterized by the opening angle \(\beta\) and the chord \(c_t\). This latter parameter is kept constant, thus the tail span is controlled via \(\beta\) from the trigonometrical relationship
$$\begin{aligned} b_{t} = 2c_{t}\tan \frac{\beta }{2}\end{aligned}$$
The main limitation of this framework is the low range of angles of attack (\(\alpha _{tail}<5^{^\circ }\)) within which it provides accurate results34. This limitation is valid in our context of fast forward flight, where the bird flight is straight, horizontal, and the forward velocity u is much larger than the vertical one w. The velocity component acting on the tail-like surface is
$$\begin{aligned} \textbf{U}_{t}(t) = \textbf{U}_{\infty } + \textbf{v}_{ind}^{w\rightarrow t} + \textbf{v}_{ind,b} \end{aligned}$$
(8)
where the term \(\textbf{v}_{ind}^{w\rightarrow t}\) is the velocity acting on the tail, induced by vortex filament shed by the wing. This velocity is calculated according to Biot-Savart law23, and
$$\begin{aligned} \textbf{v}_{ind,b} = q \hat{\textbf{e}}_{y’} \wedge (\textbf{G}-\textbf{N})\end{aligned}$$
is the velocity induced by the body rotation, with \((\textbf{G}-\textbf{N})\) the vector between the center of mass of the body (\(\textbf{G}\)) and the point of application of the forces on the tail (\(\textbf{N}\)) taken at two third along the tail chord, as illustrated in Fig. 1b. Furthermore, our tail spread takes \(45^{\circ }\) as upper value, i.e. a maximum aspect ratio of 1.65 for our triangular tail. The tail tilt relative to the body is constrained to \(0^{\circ }\), so that the tail pitch is the same as the body’s. Furthermore, as the tail kinematic velocity stays far lower than in the wing sections, the stall zone is avoided.
The forces generated by the tail are computed according to33
$$\begin{aligned} \begin{aligned} F_{x’, t}&=\left( \frac{\pi }{4}\rho \alpha _{t} |\textbf{U}_{t}|^{2} b_{t}^{2} \right) \cdot \hat{\textbf{e}}_{x’} \\ F_{z’, t}&= \left( \frac{\pi }{4}\rho \alpha _{t} |\textbf{U}_{t}|^{2} b_{t}^{2} \right) \cdot \hat{\textbf{e}}_{z’} \end{aligned} \end{aligned}$$
(9)
These forces are applied at point \(\textbf{N}\). In addition, adding this tail-like surface introduces another source of drag that needs to be accounted for. This parasitic drag contribution is, according to33
$$\begin{aligned} D_{p, t}=\frac{1}{2}\rho |\textbf{U}_{t}|^{2} S_{t} C_{D, f} \end{aligned}$$
(10)
where \(S_{t}\) is the tail planar surface and \(C_{D, f}\) the dimensionless friction coefficient.
Numerical framework
Due to the fact that Eq. (1) are driven by time-varying forces and periodic actuation, they do not converge to steady equilibrium points11,35. Instead, the steady condition corresponds to a limit cycle of this model, and its stability can be assessed via dedicated tools. The numerical identification of such a periodic orbit and the characterization of its stability are challenged by the fact that this orbit is unknown a priori. In previous work, we developed a method to achieve both at the same time via a multiple-shooting algorithm. This framework has been released as a Python toolbox called multiflap13. It takes as input a set of ordinary differential equations (ODE) like (1), i.e. of the form
$$\begin{aligned} \dot{\textbf{x}} = \textbf{v}(\textbf{x}, t, \nu ) \end{aligned}$$
(11)
with \(\textbf{x}\) being the state variables, \(\textbf{v}\) a vector field describing the system dynamics, t the time and \(\nu\) a set of configuration parameters. A periodic orbit is thus a particular solution of this set of ODE such that
$$\begin{aligned} \textbf{x}(t, \nu ) = \textbf{x}(t+T, \nu )\end{aligned}$$
with T being the cycle period. Such a periodic solution, defines a steady flight regime.
The stability of such closed orbits can also be assessed through the long-term response to a perturbation: if the perturbed trajectory converges back to the orbit, then it is stable, and vice-versa. We perform this stability analysis via the Floquet theory12. Stability of the set of equations is governed by the eigenvalues of the so-called Floquet matrix \(\mathbb {J}\), also known as Floquet multipliers, \(\Lambda _{i}\). This Floquet matrix maps perturbations within an infinitesimal sphere around a point of the limit cycle \((\mathbf {x_{0}}, t_{0})\) into an ellipsoid after a time T equal to the period of the orbit. Stretching or contracting ratios of the principal axes of this first order transformation are governed by the Floquet Multipliers. Floquet multipliers have the property of being invariant along the limit cycle, whereas the Floquet matrix and its eigenvectors depend on it. Concretely, the Floquet matrix \(\mathbb {J}\) can be calculated as the solution of the variational equation:
$$\begin{aligned} \begin{aligned} \frac{d\mathbb {J}}{dt}(\textbf{x}_0) \Big |_{t_0}^{t}&= \mathbb {A}(\textbf{x}, t) \mathbb {J}(\textbf{x}_0) \Big |_{t_0}^{t}\\ \mathbb {J} (\textbf{x}_0)\Big |_{t_0}^{t_0}&= \mathbb {I} \end{aligned} \end{aligned}$$
(12)
where the matrix
$$\begin{aligned} \mathbb {A}(\textbf{x}, t) = \nabla {\textbf{v}(\textbf{x}, t)} | _{\textbf{x}=\mathbf {x^{*}}} \end{aligned}$$
(13)
is called the Stability Matrix36 and is T-periodic on the limit cycle. If the absolute values of all Floquet multipliers \(\Lambda _{i}\) are smaller than one, the corresponding periodic orbit is stable. Conversely, if the absolute value of at least one multiplier is larger than one, the corresponding orbit is unstable and the perturbation spirals out of the limit cycle along the corresponding eigendirection(s). This framework provides another important feature, namely the stretching/contracting rate per unit of time, or Floquet exponent, \(\lambda _{i}\)36,37
$$\begin{aligned} \lambda _{i} = \frac{1}{T} \ln |\Lambda _{i}| \end{aligned}$$
(14)
Application to steady and level flight
In the present study, we leverage the aforementioned multiple-shooting algorithm, for seeking steady flight regimes within the following parametric space
$$\begin{aligned} \nu = (\beta , \psi _{s,z}, A_{s, x}, q_{w,y}) \end{aligned}$$
(15)
where \(\beta\) is the tail opening angle, \(\psi _{s,z}\) is the shoulder sweep offset, \(A_{s, x}\) is the wingbeat amplitude, and \(q_{w,y}\) is the mean rotation angle of the wing profiles of the forearm about the axis y, see Fig. 2. The other parameters defining the wing kinematics are kept fixed to values similar to those reported in13.
(a) Front view of the bird model. The wingbeat amplitude \(A_{s,x}\) about the \(x’\)-axis is the main kinematic parameter governing the wing amplitude of movement. (b) Top view of the bird model. The shoulder sweep offset \(\psi _{s,z}\) captures the average angle of the arm bone with respect to \(x’_{w}\) over the period. The angle \(\beta\) captures the magnitude of tail opening. (c) Section of the wing. It highlights the wing rotation angle governed by the wrist degree of freedom \(q_{w,y}\).
Previous studies have shown that these four parameters decisively govern the flight regime in bird flapping and gliding modes. The tail opening \(\beta\) and shoulder sweep offset \(\psi _{s,z}\) influence flight stability, since these are the parameters having a paramount influence on the generation of pitching moment. Then, the shoulder wingbeat amplitude \(A_{s,x}\) has a direct impact on thrust production and therefore on airspeed and power consumption38. The last parameter \(q_{w,y}\) modulates the generation of lift13,15.
On top of seeking for steady flight regimes, it is important to identify those corresponding to level flight, i.e. with the bird flying at a constant altitude. This level flight condition thus corresponds to an average mean vertical velocity being equal to zero over the period, in a fixed reference frame. Figure 1 shows that this instantaneous vertical velocity can be computed as
$$\begin{aligned} W_{ff} = -\dot{Z} = u\sin {\theta } – w\cos {\theta } \end{aligned}$$
(16)
Concretely, satisfying the level flight condition isolates a three-dimensional manifold within the four-dimensional parametric space of Eq. (15). Finding this manifold is done by searching the value of the parameter \(q_{w,y}\) that corresponds to level flight for all possible combinations of the three other parameters \(\beta\), \(\psi _{s,z}\) and \(A_{s,x}\). In other words, we report here only the limit cycles that belong to the manifold satisfying
$$\begin{aligned} \left| \overline{{W}_{ff}}({\beta }, {\psi _{s,z}}, {A}_{s, x},q_{w,y}) \right| < \epsilon \end{aligned}$$
(17)
with \(\epsilon = 5 \times {10}^{-3}\) ms\(^{-1}\), which corresponds to a maximum vertical deviation of 1 mm per flapping cycle.
Power consumption and cost of transport
Each limit cycle corresponds to a particular flapping gait with its own mechanical power consumption and the corresponding cost of transport.
Since inertial power for accelerating and decelerating a wing is neglected, the actuation power produced by the wing joints is exactly equal to the power transferred by this wing to the environment, i.e.
$$\begin{aligned} P_{act}(t) = \sum _{i=1}^{n}\textbf{F}_{aero,i}(t) \cdot (-\textbf{v}_{kin, i}(t)) \end{aligned}$$
(18)
where \(\textbf{v}_{kin, i}(t)\) is the velocity of the lifting line computed at the discretized point \(P_i\) and time t, and \({F}_{aero,i}(t)\) is the corresponding aerodynamic load on the wing element i, computed by the quasi-steady lifting line model. The mean power consumption over one flapping period is thus
$$\begin{aligned} \overline{P}_{act} = \frac{1}{T}\int _{0}^{T} P_{act}(t) dt \end{aligned}$$
(19)
Another important metric to assess locomotion performance is the so-called Cost of Transport (CoT), i.e. a dimensionless ratio equal to the mechanical work produced by the actuators to transport a unit of body weight across a unit of distance39. Here, it is thus defined as
$$\begin{aligned} CoT = \frac{\overline{P}_{act}}{mg\overline{|\textbf{U}_{\infty }|}} \end{aligned}$$
(20)
with \(\overline{|\textbf{U}_{\infty }|}\) being the magnitude of the flight speed of the corresponding limit cycle, averaged over one period.
Parametrization of bird morphology and wing kinematics
The above framework and concepts are here embodied in a model of the northern bald ibis (Geronticus eremita). This species does not only display features at the origin of possible conflict between energetic and stability demands (continuous flapping propulsion and long migratory flights), but it is often studied by biologists because of its endangered status. The length of bones and feathers are set up accordingly. To the best of our knowledge, the precise wingbeat kinematics of this particular bird have not been reported in the literature. Consequently, we follow the same approach as in13, consisting in scaling up the detailed kinematic pattern of another bird40 to the morphology of ours. This scaling process, and the parameter tuning are extensively detailed in the Supplementary Material S1. The morphological and kinematic parameters used to describe the bird are reported in Table S1 and the resulting flapping gait is illustrated in Figs. S1, S2 and S3.
This study is performed within the following parametric space: tail opening \(\beta \in [0^{\circ },45^{\circ }]\), wingbeat amplitude \(A_{s,x} \in [29^{\circ },45^{\circ }]\) and sweep offset \(\psi _{s,z} \in [9^{\circ },15^{\circ }]\). The parametric space is meshed with an uniform grid spaced along \(\psi _{s,z}\) and \(A_{s,x}\) with a step size of \(0.5^{\circ }\), and a step size of \(1^{\circ }\) along \(\beta\). This resulted in 19,734 possible flight configurations. In the results, we report all solutions satisfying Eq. (17), with the addition of two exclusion criteria. First, we excluded limit cycles that do not correspond to fast forward flight. In41, this corresponds to flying modes 4 and 5 and requires the body pitch angle to stay close to the horizon tail configuration. Concretely, we excluded from the results limit cycles corresponding to a mean body pitch angle larger than \(6^{\circ }\) in absolute value. Second, we excluded limit cycles corresponding to biologically incompatible kinematics. This was implemented by excluding limit cycles with a mean rotation angle of the wing \(q_{w,y}\) larger than \(12^{\circ }\) in absolute value. Indeed, remembering that the related amplitude of this joint was fixed to \(30^{\circ }\) (Table S1 in the Supplementary Material), this criteria excluded solutions corresponding to geometrical rotation of the forearm larger than \(\pm 42^{\circ }\) in absolute value, which we considered to be not physiologically consistent.
