On the role of tail in stability and energetic cost of bird flapping flight
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 flightFlight 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 – gsin 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 + gcos 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)
Figure 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).Full size imageThe 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 modelThe 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} = qhat{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 (2pi) that saturates at an effective angle of attack (alpha _{r}) of (pm 15^{circ }).Drag production by body and wingThe 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 modelSince 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} More
