2.1 Three-dimensional dynamic flight simulation

To simulate three-dimensional free flight with flapping wings, we solve the Newton–Euler equations for the coupled wing–body system (Chang & Wang Reference Chang and Wang2014). The insect model consists of $(n+1)$ rigid bodies, where $n$ is the number of wings on the body. Each wing is modelled as an ellipsoid connected to the body, also an ellipsoid, through a ball joint that allows for three degrees of freedom in rotation. The body kinematics are given by its position $\boldsymbol{r}^{b}$ , linear velocity $\boldsymbol{v}^{b}$ , body orientation quaternion $[q^{b}]$ and angular velocity $\unicode[STIX]{x1D74E}^{b}$ . In our current implementation we use quaternions to represent the body and wing orientations. This has the advantage of avoiding gimbal lock and simplifies the algebra. For the results the quaternions are converted to Euler angles, which are easier to understand, as they refer to the rotations about body axes.

The Newton–Euler equations governing the body dynamics are

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle m^{b}\boldsymbol{a}^{b}=m^{b}\boldsymbol{g}-\mathop{\sum }_{i=1}^{n}\boldsymbol{F}_{i}^{c}, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D644}^{b}\unicode[STIX]{x1D737}^{b}=-\unicode[STIX]{x1D74E}^{b}\times (\unicode[STIX]{x1D644}^{b}\unicode[STIX]{x1D74E}^{b})-\mathop{\sum }_{i=1}^{n}\unicode[STIX]{x1D749}_{i}^{c}-\mathop{\sum }_{i=1}^{n}\boldsymbol{r}_{i}^{b}\times \unicode[STIX]{x1D749}_{i}^{c}. & \displaystyle\end{eqnarray}$$

Similarly, the governing equations for the $i$ th wing are

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle m_{i}^{w}\boldsymbol{a}_{i}^{w}=m_{i}^{w}\boldsymbol{g}+\boldsymbol{F}_{i}^{c}+\boldsymbol{F}_{i}^{a}, & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D644}_{i}^{w}\unicode[STIX]{x1D737}_{i}^{w}=-\unicode[STIX]{x1D74E}_{i}^{w}\times (\unicode[STIX]{x1D644}_{i}^{w}\unicode[STIX]{x1D74E}_{i}^{w})+\unicode[STIX]{x1D749}_{i}^{c}+\boldsymbol{r}_{i}^{w}\times \unicode[STIX]{x1D749}_{i}^{c}+\unicode[STIX]{x1D749}_{i}^{a}. & \displaystyle\end{eqnarray}$$

Here $b$ denotes body, $w$ denotes wing, $m$ is mass, $\unicode[STIX]{x1D644}$ is the moment of inertia tensor, $\boldsymbol{a}$ is the linear acceleration, $\unicode[STIX]{x1D737}$ is the angular acceleration, $\boldsymbol{g}$ is the gravitational constant, $\unicode[STIX]{x1D74E}$ is the angular velocity, $\boldsymbol{F}^{a}$ and $\unicode[STIX]{x1D749}^{a}$ are the aerodynamic force and torque on the wing, $\boldsymbol{F}_{i}^{c}$ and $\unicode[STIX]{x1D749}_{i}^{c}$ are the internal force and torque to be determined, $\boldsymbol{r}_{i}^{b}$ is the position of the $i$ th wing root relative to the body centre of mass (COM) and $\boldsymbol{r}_{i}^{w}$ is the position of the $i$ th wing root relative to the COM of the $i$ th wing.

Because we have introduced unknown variables in the form of the internal forces and torques at the joint, $\boldsymbol{F}^{c}$ and $\unicode[STIX]{x1D749}^{c}$ , additional equations must be specified. These equations are two kinematic constraints applied at each wing joint. The first is on the angular acceleration of the wing relative to the body,

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D737}_{i}^{r}=\unicode[STIX]{x1D737}_{i}^{w}-\unicode[STIX]{x1D737}^{b},\end{eqnarray}$$

which states that the wing acceleration relative to the body must be the prescribed acceleration. The second constraint is the matching condition for the linear acceleration of the wing hinge:

(2.6) $$\begin{eqnarray}\boldsymbol{a}^{b}+\unicode[STIX]{x1D737}^{b}\times \boldsymbol{r}_{i}^{b}+\unicode[STIX]{x1D74E}^{b}\times (\unicode[STIX]{x1D74E}^{b}\times \boldsymbol{r}_{i}^{b})=\boldsymbol{a}^{w}+\unicode[STIX]{x1D737}^{w}\times \boldsymbol{r}_{i}^{w}+\unicode[STIX]{x1D74E}^{w}\times (\unicode[STIX]{x1D74E}^{w}\times \boldsymbol{r}_{i}^{w}).\end{eqnarray}$$

At each instance in time, the coupled dynamic equations (2.1)–(2.4) together with the constraint equations (2.5) and (2.6) can be cast into a linear system, $\text{A}x=\text{B}$ , where $x=[\boldsymbol{a}^{b},\unicode[STIX]{x1D737}^{b},\boldsymbol{a}_{i}^{w},\unicode[STIX]{x1D737}_{i}^{w},\boldsymbol{F}_{i}^{c},\unicode[STIX]{x1D749}_{i}^{c}]$ , and $\text{A}$ contains $m_{i}$ and $\unicode[STIX]{x1D644}_{i}$ , and the vector $\text{B}$ contains the known variables with $\unicode[STIX]{x1D74E}_{i}$ , $\boldsymbol{F}_{i}^{a}$ and $\unicode[STIX]{x1D749}_{i}^{a}$ evaluated at the current time. We solve this system of equations by inverting the matrix $\text{A}$ , using standard lower–upper (LU) decomposition, to find the body accelerations.

Once the body accelerations $\boldsymbol{a}^{b}$ and $\unicode[STIX]{x1D737}^{b}$ are obtained, the body kinematic state evolves in time according to

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{\boldsymbol{r}}^{b}=[q]^{b}(\boldsymbol{v}^{b}), & \displaystyle\end{eqnarray}$$

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{\boldsymbol{v}}^{b}=\boldsymbol{a}^{b}-\unicode[STIX]{x1D74E}^{b}\times \boldsymbol{v}^{b}, & \displaystyle\end{eqnarray}$$

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle [\dot{q}]^{b}={\textstyle \frac{1}{2}}[q]^{b}\cdot [\unicode[STIX]{x1D714}]^{b}, & \displaystyle\end{eqnarray}$$

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D74E}}^{b}=\unicode[STIX]{x1D737}^{b}. & \displaystyle\end{eqnarray}$$

Here equation (2.7) rotates the velocity vector from the body frame to the laboratory frame. This quaternion rotation of a vector can be written out in normal vector products by separating the real part, $q_{0}$ , and complex part, $\boldsymbol{q}$ , of the quaternion as $[q]=(q_{0},\boldsymbol{q})$ . Equation (2.7) becomes

(2.11) $$\begin{eqnarray}\dot{\boldsymbol{r}}^{b}=\boldsymbol{v}^{b}+2\boldsymbol{q}^{b}\times (\boldsymbol{q}^{b}\times \boldsymbol{v}^{b}+q_{0}^{b}\boldsymbol{v}^{b}).\end{eqnarray}$$

In (2.9) the operator is quaternion multiplication, which is possible after promoting $[\unicode[STIX]{x1D714}]^{b}$ to a quaternion by adding a zero real part, i.e.  $[\unicode[STIX]{x1D714}]^{b}=(0,\unicode[STIX]{x1D74E}^{b})$ . This set of nonlinear ordinary differential equations is solved using an adaptive order Adams–Bashforth–Moulton algorithm that ensures the quaternion $[q]^{b}$ remains unitary and thus a representation of the body orientation (Shampine & Gordon Reference Shampine and Gordon1975).

2.1.1 Aerodynamic force model on the wing

The dynamic equations, (2.1)–(2.4), require the specification of the aerodynamic force and torque on the wings. Our force model is a quasi-steady force approximation that takes into account the main unsteady effects, including the dynamic stall, the coupling between wing translation and rotation and the added mass term (Wang Reference Wang2005). The forces on the wings are decomposed into three parts, corresponding to the lift, the drag and the added mass, i.e.  $\boldsymbol{F}_{i}^{a}=\boldsymbol{F}_{L}+\boldsymbol{F}_{D}+\boldsymbol{F}_{A}$ . We use the blade element approximation and divide the wing into infinitesimal wing blades, $\boldsymbol{F}=\int \boldsymbol{f}(s)\,\text{d}s$ and $\unicode[STIX]{x1D749}=\int \boldsymbol{s}\times \boldsymbol{f}(s)\,\text{d}s$ , with $s$ the coordinate along the span.

The specific form of the lift and drag on each blade element is deduced from experiments and numerical simulations of a free falling plate in fluid (Pesavento & Wang Reference Pesavento and Wang2004; Andersen, Pesavento & Wang Reference Andersen, Pesavento and Wang2005). The lift force density is given by $\boldsymbol{f}_{L}(s)=-\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D71E}(s)\times \boldsymbol{v}(s)$ , with the circulation $\unicode[STIX]{x1D71E}(s)=(c(s)/2)[C_{T}(\unicode[STIX]{x1D6FC})v(s)\sin (2\unicode[STIX]{x1D6FC}(s))-C_{R}(\unicode[STIX]{x1D6FC})c(s)\unicode[STIX]{x1D714}_{x}]\hat{\boldsymbol{s}}$ . Here $\hat{\boldsymbol{s}}$ is the outward unit vector along the span, $\unicode[STIX]{x1D70C}_{f}$ is the fluid density, $C_{T}$ and $C_{R}$ are force coefficients that depend on $\unicode[STIX]{x1D6FC}$ , $c(s)$ is the local chord length, $\unicode[STIX]{x1D714}_{x}$ is the $x$ -component of the angular velocity and $v(s)$ is the magnitude of the local wing velocity. $\unicode[STIX]{x1D6FC}$ denotes the angle of attack, which is the angle between the wing chord and the velocity vector. The damping force density is given by: $\boldsymbol{f}_{D}=-\unicode[STIX]{x1D70C}_{f}k(s)\boldsymbol{v}(s)$ , with $k(s)=(1/2)c(s)v(s)(\text{A}-\text{B}\cos (2\unicode[STIX]{x1D6FC}(s)))$ and $\text{A}$ and $\text{B}$ are constants. The aerodynamic force parameters are fixed at $C_{T}=1.5$ , $C_{R}=\unicode[STIX]{x03C0}$ , $\text{A}=1.4$ , $\text{B}=1$ and $\unicode[STIX]{x1D70C}_{f}=1.293~\text{kg}~\text{m}^{-3}$ .

The added mass term represents the force needed to accelerate the wing in a potential flow. For ellipsoids, the added mass and torque have exact expressions:

(2.12) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}_{added}=\left(\begin{array}{@{}c@{}}-m_{33}v_{z}\unicode[STIX]{x1D714}_{y}\\ m_{33}v_{z}\unicode[STIX]{x1D714}_{x}\\ -m_{33}\dot{v}_{z}\end{array}\right), & \displaystyle\end{eqnarray}$$

(2.13) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D749}_{added}=\left(\begin{array}{@{}c@{}}-m_{44}\dot{\unicode[STIX]{x1D714}}_{x}+m_{55}\unicode[STIX]{x1D714}_{y}\unicode[STIX]{x1D714}_{z}-m_{33}v_{y}v_{z}\\ -m_{55}\dot{\unicode[STIX]{x1D714}}_{y}-m_{44}\unicode[STIX]{x1D714}_{x}\unicode[STIX]{x1D714}_{z}+m_{33}v_{x}v_{z}\\ (m_{44}-m_{55})\unicode[STIX]{x1D714}_{x}\unicode[STIX]{x1D714}_{y}\end{array}\right). & \displaystyle\end{eqnarray}$$

Here $v$ and $\unicode[STIX]{x1D714}$ are the velocity and angular velocity at the centre of mass of the wing. The added mass coefficients, $m_{33}$ , $m_{44}$ and $m_{55}$ , are (Tuckerman Reference Tuckerman1926)

(2.14) $$\begin{eqnarray}\displaystyle & \displaystyle m_{33}=\frac{4}{3}\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}_{f}\frac{ab^{2}}{E(k)}, & \displaystyle\end{eqnarray}$$

(2.15) $$\begin{eqnarray}\displaystyle & \displaystyle m_{44}=\frac{4}{15}\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}_{f}ab^{4}\frac{a^{2}-b^{2}}{(2a^{2}-b^{2})E(k)-b^{2}K(k)}, & \displaystyle\end{eqnarray}$$

(2.16) $$\begin{eqnarray}\displaystyle & \displaystyle m_{55}=\frac{4}{15}\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}_{f}a^{3}b^{2}\frac{a^{2}-b^{2}}{(a^{2}-2b^{2})E(k)+b^{2}K(k)}, & \displaystyle\end{eqnarray}$$

with $K(k)$ and $E(k)$ the elliptic integrals of the first and second kind respectively. $k=1-(b/a)^{2}$ , with $a$ the semi-major axis, half-span, and $b$ the semi-minor axis, half-chord.

2.2 Wing motion

The wings flap back and forth along a horizontal stroke plane following a pattern similar to that observed in flies. We use a generalized family of wing motions for the three wing angles: $\unicode[STIX]{x1D719}^{w}(t)$ , the stroke angle, $\unicode[STIX]{x1D703}^{w}(t)$ , the deviation angle and $\unicode[STIX]{x1D713}^{w}(t)$ , the wing-pitch angle (Berman & Wang Reference Berman and Wang2007),

(2.17) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}^{w}(t)=\unicode[STIX]{x1D719}_{0}+\unicode[STIX]{x1D719}_{m}\frac{\arcsin (K\sin (2\unicode[STIX]{x03C0}ft))}{\arcsin (K)}, & \displaystyle\end{eqnarray}$$

(2.18) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}^{w}(t)=\unicode[STIX]{x1D703}_{0}+\unicode[STIX]{x1D703}_{m}\cos (N\cdot 2\unicode[STIX]{x03C0}ft+\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D703}}), & \displaystyle\end{eqnarray}$$

(2.19) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}^{w}(t)=\unicode[STIX]{x1D713}_{0}+\unicode[STIX]{x1D713}_{m}\frac{\tanh (C\sin (2\unicode[STIX]{x03C0}ft+\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D713}}))}{\tanh (C)}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D719}_{0}$ , $\unicode[STIX]{x1D703}_{0}$ and $\unicode[STIX]{x1D713}_{0}$ set the mean; $\unicode[STIX]{x1D719}_{m}$ , $\unicode[STIX]{x1D703}_{m}$ and $\unicode[STIX]{x1D713}_{m}$ are amplitudes; $f$ is the wing-beat frequency; $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D703}}$ and $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D713}}$ are phase shifts; $N=1$ or $2$ , $0<K<1$ and $C>0$ are waveform control parameters. $N=1$ corresponds to one vertical oscillation per stroke, and $N=2$ corresponds to a figure-eight motion. $\unicode[STIX]{x1D719}^{w}$ becomes sinusoidal when $K$ is close to $0$ and triangular when $K$ is close to $1$ . $\unicode[STIX]{x1D713}^{w}$ becomes sinusoidal at small $C$ and a step function at large $C$ . For simplicity, the deviation $\unicode[STIX]{x1D703}$ from the main stroke plane is neglected, and the phase shift $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D713}}$ in $\unicode[STIX]{x1D713}^{w}(t)$ is such that the wing pitches in advance of the wing stroke reversal. For the calculations reported here, the base parameters are $f=250~\text{Hz}$ , $\unicode[STIX]{x1D719}_{m}=70^{\circ }$ , $\unicode[STIX]{x1D719}_{0}=0^{\circ }$ , $K=0.7$ , $\unicode[STIX]{x1D703}_{m}=\unicode[STIX]{x1D703}_{0}=0^{\circ }$ , $\unicode[STIX]{x1D713}_{m}=53^{\circ }$ , $\unicode[STIX]{x1D713}_{0}=90^{\circ }$ , $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D713}}=72.4^{\circ }$ and $C=2.4$ (Berman & Wang Reference Berman and Wang2007; Chang & Wang Reference Chang and Wang2014). When a parameter is varied in the study, it is stated in the relevant figure.

2.3 Morphological parameters

For comparison with related studies we use morphological parameters similar to those of fruit flies, Drosophila melanogaster, we start with the same set of parameters as in (Chang & Wang Reference Chang and Wang2014): body mass $m_{b}=1.1~\text{mg}$ , length $=2.4~\text{mm}$ , width $=1.2~\text{mm}$ ; wing mass $m_{w}=3.6\times 10^{-3}~\text{mg}$ , span $=2~\text{mm}$ , maximal chord $=1~\text{mm}$ and maximal thickness $=0.1~\text{mm}$ . The wings are attached $h=1~\text{mm}$ vertically up from and $0.36~\text{mm}$ to the side of the body centre of mass (COM), with the COM of the wing $s_{x}=0.87~\text{mm}$ along the span and $s_{y}=0.25~\text{mm}$ along the chord direction from the hinge location (see figure 1 d). When a parameter is varied, it is stated in the relevant figure. For theoretical interest, the hinge position is chosen to be in the range of $h\in [-2,8]~\text{mm}$ , which corresponds to 1 body length below and 4 body lengths above the centre of mass of the fly.

Figure 1.

(a) The main mode of flight stability is the pitch instability, and it is due to the coupling between the forward motion and the body pitching motion. The stick shows the body pitch axis, with the circle indicating the top of the body. In both hovering and ascending flight, the body oscillation grows and the flight becomes unstable. Only in descending flight can an upward drag force stabilize the oscillation. (b) The Poincaré section for longitudinal flight. The state variables are $(v_{x},v_{y},\unicode[STIX]{x1D703},\unicode[STIX]{x1D714})$ . The flight model maps the initial state $s(t)$ to $s(t+T)$ after one wing beat. Periodic flight corresponds to fixed points of this map. (c) The wing motion has two degrees of freedom, the back and forth stroke $\unicode[STIX]{x1D719}^{w}$ , and the rotation about the pitching axis, the wing pitch $\unicode[STIX]{x1D713}^{w}$ . Left and right symmetry ensures that the body moves along the longitudinal plane. The tilt of the body away from the vertical is $\unicode[STIX]{x1D703}_{b}$ . (d) Morphological parameters. $h$ denotes the wing hinge position (red dot above the body centre of mass (COM)). $s_{x}$ and $s_{y}$ define the location of the wing COM relative to the hinge.

2.4 Periodic orbit search and stability computation

The stability of periodic flights are analysed using Floquet analysis. For this, we first need to identify periodic flight states (figure 1 b). The left and right wings flap symmetrically leading to longitudinal flights. The state vector is four-dimensional, $\boldsymbol{s}=[v_{x},v_{z},\unicode[STIX]{x1D703}_{b},\unicode[STIX]{x1D714}_{y}]$ , where $v_{x}$ and $v_{z}$ are the horizontal and vertical velocities, and $\unicode[STIX]{x1D703}_{b}$ and $\unicode[STIX]{x1D714}_{y}$ are the body pitch angle and the pitching rate. A periodic state is found by adjusting the initial state $\boldsymbol{s}_{0}$ until it minimizes the periodicity error $\unicode[STIX]{x1D716}=\Vert M(\boldsymbol{s}_{0})-\boldsymbol{s}_{0}\Vert$ , using root-finding algorithms. Here $M$ represents the evolution of flight dynamics, $\boldsymbol{s}(t+T)=M(\boldsymbol{s}(t)),$ with $T$ the wing-beat period. The search algorithm consists of successive minimal searches, with each cycle combining the previous best periodic state with periodic states of nearby model parameters to seed a new minimum search.

The function $M$ defines a Poincaré map on our flight state space. The periodic flight states correspond to fixed points of this map. We also investigated whether other periodic orbits exist, by starting from very different initial conditions and looking for period-2 orbits. All results converge to the same periodic states as reported here, suggesting a unique solution.

The stability of flight is given by the eigenvalues of the Jacobian of $M$ . The eigenvalues of the Jacobian matrix $\unicode[STIX]{x1D63F}$ measure the growth rate of the deviation from periodic flight (Grimshaw Reference Grimshaw1991; Kuznetsov Reference Kuznetsov2013). To find these eigenvalues we first compute $\unicode[STIX]{x1D63F}$ ,

(2.20) $$\begin{eqnarray}\unicode[STIX]{x1D63F}=\left.\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{s}}\left(\frac{M(T,\boldsymbol{s})-\boldsymbol{s}}{T}\right)\right|_{s=s_{0}}.\end{eqnarray}$$

The eigenvalues $\unicode[STIX]{x1D706}$ of $\unicode[STIX]{x1D63F}$ are known as the Lyapunov exponents. If $\max (\text{Re}(\unicode[STIX]{x1D706}))<0$ the deviations decrease over time and the periodic flight is stable. Moreover, the eigenvectors of $\unicode[STIX]{x1D63F}$ will indicate which deviations lead to instability.

For ease of interpretation, we define $\unicode[STIX]{x1D6FC}=\exp (\unicode[STIX]{x1D706}T)$ , which is the multiplicative factor by which a deviation grows after a single wing beat. These $\unicode[STIX]{x1D6FC}$ values are called the multipliers of the Poincaré map associated with our model $M$ . If all $|\unicode[STIX]{x1D6FC}_{i}|<1$ , then the deviations decay and the flight is stable. However, if there is at least one $|\unicode[STIX]{x1D6FC}_{i}|>1$ then the flight is unstable.