Electronic Supplemental Material to “Non-linear flight dynamics and stability ofhovering model insects”

Bin Liang and Mao Sun *

Ministry-of-Education Key Laboratory of Fluid Mechanics,Beijing University of Aeronautics and Astronautics, China(*)

S1Equations of motion of an insect

Equations of motion for a model insect with a body and N flapping rigid wings have been derived by Sun etal. [1] (also see [2]). Three frames of reference have been used in the derivation (Fig. S1): Frame (xE, yE, zE) is an inertial frame (fixed on the earth);xE is in the horizontal direction and zE in the vertical direction.Frame (xb, yb, zb) is a frame fixed on the insect body with its origin at the centre of mass of the wingless body;xband zb are in the longitudinal symmetrical plane and yb points to the right side of the insect. Frame (xw, yw, zw) is a frame fixed on an insect wing with its origin at the root of the wing; xw is in the spanwise direction, pointing to the wing tip, and zw in the chordwise direction, pointing to the trailing edge of the wing. With respect to the inertial frame, the centreof mass of the wingless body moves at velocity vcg, the body rotates at angular velocity ωbd(angular velocity of the wing relative to the body is prescribed). Let FA be the total aerodynamic force and MA the total aerodynamicmoment about the center of mass of the body, m the total mass of the insect (body mass plus wing mass), mwg the mass of a wing, Ibd the matrix of moments and products of inertial of the body, Iwg the moments and products of inertial of a wing, Rh the vector from the body center of mass to the root of a wing, Rwg the vector from wing root to the center of mass of a wing, and gthe gravitational acceleration and t the time. The equations of motionderived by Sun et al. [1]can be written as:

(S1a)

(S1b)

where

(S2a)

(S2b)

(S2c)

(S2d)

(S2e)

where bvcg, bωbd, bFA, bMA and bIbd, represent the xb, yb and zb components of vcg and ωbd, FA, MA, Ibd, respectively; N is the number of wings; Iwg is the moments and products of inertial of a wing with respect to the wing frame (,,);wIwg represents thexw, yw and zw component of Iwg; and are matrix of direction cosines (determined by the flapping angles), representing the transformation from wing frame to body frame and body frame to wing frame, respectively; symbol “~” is defined as follows: let be a vector, then denotes the following matrix (see [3]):

(S3)

InEq. S1,a1represents the inertial force of the wings due to the rotation and angular acceleration of the body and a2 represents the inertial forces due to wing acceleration and the gyroscopic forces due to the wing and body rotations; b1 and b2represent the corresponding moments;represents the time-varying moments and products of inertia of the flapping wings. The second term on the left-hand-side of Eq. S1b is the moment about the center of mass of the body due to the weight of the wings.

Let u, v and w denote the xb, yband zb components of vcg, respectively, and p, qand r denote the xb, yband zb componentsof ωbd, respectively. Let XA, YA and ZA denote the xb, yband zb components of FA, respectively, and LA, MAand NA denote the xb, yband zb components of MA, respectively. Let XI, YI and ZI denote the xb, yband zb components of (a1+b1), respectively, and LI, MI and NIdenote the xb, yband zb components of (a2+b2), respectively(XI,andZIare called as wing-inertial forces and MI is called as wing-inertial moment). Let MXG, MYG and MZGdenote thexb, yband zb components of the moment produced by the weight of the wings. The xb, yband zb components of Eq. S1 are:

(S4a)

(S4c)

(S4d) (S4e)

(S4f)

where andare thepitch and roll angles of the body, respectively, andIx,bw, Iy,bw,etc., are the elements ofof (Ibd+c):

(Ibd+c)=

As seen above, some terms in EqsS4a-f are dependent on Euler angles and  and equations that govern the Euler angle rates are needed, and they are [10]:

(S4g)

(S4h)

(S4i)

It is noted that when  is close to +90° (or -90°), some terms in the right side of EqsS4g and i become infinity. This problem is avoided by using quaternions to replace the Euler angles in EqsS4g-i.

u, vw, p, q, r, ,and , which determine the motion of the flying system, are called as state variables.

S2 Solution method of the Navier-Stokes equations

To obtain the aerodynamic forces and moments, in principle we need to compute the flows around the wings and the body. But near hovering, the aerodynamic forces and moments of the body are negligibly small compared to those of the wings because the velocity of the body motion is very small [4]. Furthermore, the interaction between wing and body has been shown to be negligibly small: the aerodynamic force in the case with the body/wing interaction is less than 2.5% different from that without body/wing interaction [5]. Therefore, in the present CFD model, the body is neglected.

Eq.2.2 is solved over moving overset gridsbecause there are relative motion between the left and right wings.The computational method is the same as thatused bySun and Yu [6] and only an outline of the method is given here. The algorithm is based on the method of artificial compressibility. It was first developed by Rogers et al. [7]for single-grid, and then extended by Rogers and Pulliam [8] to overset grids.The time derivatives of the momentum equations are differenced using a second-order, three-point backward difference formula. To solve the time discretized momentum equations for a divergence free velocity at a new time level, a pseudo-time level is introduced into the equations and a pseudo-time derivative of pressure divided by an artificial compressibility constant is introduced into the continuity equation. The resulting system of equations are iterated in pseudo-time until the pseudo-time derivative of pressure approaches zero, thus, the divergence of the velocity at the new time level approaches zero. The derivatives of the viscous fluxes in the momentum equation are approximated using second-order central differences. For the derivatives of convective fluxes, upwind differencing based on the flux-difference splitting technique is used. A third-order upwind differencingis used at the interior points and a second-order upwind differencingis used at points next to boundaries. With overset grids (Fig. S2), for each wing there is a body-fitted curvilinear grid, which extends a relatively short distance from the wing, and in addition, there is a background Cartesian grid, which extends to the far-field boundary of the domain. The solution method for single-grid is applied to each of these grids; data are interpolated from one grid to another at the inter-grid boundary points.Details of this algorithm can be found in Refs. Refs. [7, 8].For far field boundary conditions, at the inflow boundary, the velocity components are specified as free-stream conditions while pressure is extrapolated from the interior; at the outflow boundary, pressure is set equal to the free-stream static pressure and the velocity is extrapolated from the interior. On the wing surfaces, impermeable wall and no-slip boundary conditions are applied, and the pressure on the boundary is obtained through the normal component of the momentum equation written in the moving coordinate system.

In the present study, the background grid hasonly traslational motion and it moves at the velocity of the center of mass of the insect. The wing grids are generated using a Poisson solver which is based on the work of Hilgenstock [9]; they are of O-H type grids. The background Cartesian grid is generated algebraically.

Before proceeding to compute the flows for obtaining the solutions of the equilibrium flight and the disturbance motions, we conducted grid resolution test to ensure that the flow calculation is grid independent. For a clear description of the time courses of the forces, we express the time during a cycle as a non-dimensional parameter, , such that =0 at the start of a downstroke, and =1 at the end of the subsequentupstroke.Fig. S3 shows the computed lift, drag and pitch moment coefficients of the wings of the hawkmoth at hovering flight; results calculated with three grid-systems are plotted. In all the grid-systems, the outer boundary of the wing-grid is set at about2.6c from the wing surface and that of the background-grid at20cfrom the wings. For grid-system 1, the wing grid has dimensions 47*34*30 around the wing, in the normal directionand in the spanwise direction, respectively (first layer grid thickness was 0.0015c), and the background grid has dimensions63*63*63 in thex, y andz directions, respectively. For grid-system 2, the corresponding grid dimensions are 73*50*45 and 97*97*97 (first layer grid thickness was 0.001c). For grid-system 3, the corresponding grid dimensions are109*75*67 and 145*145*145(first layer grid thickness was 0.00067c).For all the three grid-systems, grid points of the background grid concentrate in the near field of the wings where its grid density is approximately the same as that of the outer part of the wing-grid. It is observed that the first grid refinement produces only small changes in the forces and moments, and the second grid refinement produces almost no changes.Calculations were also conducted using a larger computational domain. The domain was enlarged by adding more grid points to the outside of the background grid of grid-system 2. The calculated results showed that there was no need to put the outer boundary further than that of grid-system 2. From the above discussion, it is concluded that grid-system 2 is proper for the model dronefly. Similar grid test was conducted for the model dronefly and it was concluded that a grid-system with a 71*43*58 wing-grid and a 105*105*105background-grid was was proper. The effect of time step value was considered and it was found that a numerical solution effectively independent of the time step was achieved if (is the non-dimensional time, defined as =tc/U); therefore, , is used in the present study.

S3Integration method of coupling the equations of motion with the Navier-Stokes equations and method for obtaining periodic solution

S3.1 The integration method of coupling the equations of motion with the Navier-Stokes equations

The equations of motion (Eq.2.1) and the Navier-Stokes equations (Eq.2.2) must be coupled in the solution process because the aerodynamic forces and moments in the equations of motion must be obtained from the solution of the Navier-Stokes equations and the boundary condition of Navier-Stokes equations must be obtained from the motion of the insect determined by the equations of motion. The integration process is as follows. Suppose that at time tnthe motion of the insect body (u, v,w, p, q, r, , and ) are known (the wing motion with respect to the body is prescribed), then the boundary and initial conditions of the Navier-Stokes equations would be known, and also suppose that the flow beforetnis known. Then the Navier-Stokes equations can be solved to provide the aerodynamic forces and moments at tnand Eq.2.1 is marched to the next time station tn+1, using the Euler predictor-corrector method (which has second order time accuracy).The process is repeated in the next step, and so on (see below for how the calculation starts).

S3.2 The method of finding periodic solution

For hovering flight, the mean velocity and mean pitch rate are zero, and the mean pitch angle is a constant. Therefore, we look for periodic solutions of Eq. 2.1 with zero mean horizontal velocity, vertical velocity and pitch rate. If the periodic solutions are stable, a forward integration of the equations would converge to the solutions. But linear stability analysis showed that hovering and forward flight were unstable [10]. Therefore, the shooting method is used to find the periodical solutions. In the shooting method, the initial conditions are iterated, so that they coincide with the terminal conditions. Let x represent the vector of the state variables [uvwpqr]T and f(x, t) represent the vector function on the right side of the system of Eq. 2.1. The equations of motion and the periodical conditions can be written as:

(S5a)

(S5b)

where T is the period of the solution (in the present study, T is the wingbeat period and is known). A vector parameter, s, is introduced, and it is chosen as. Eq. S5 can be written as

(S6a)

(S6b)

where g(s) is a vector function representing the dependence of the equations on s,. s is iterated using the Newton’s method until the boundary conditions, i.e. , are satisfied to an acceptableaccuracy (see [2] for details of the shooting method).

S3.3How to determine the initial flow field and the solution process

In the shooting method, the equations of motion are repeatedly integrated form t=0 to t=T with different initial state variables (u, v, w, p, q, r, , and ), until the state variables at t=T equal to those at t=0. During each of the iterations, the Navier-Stokes equations is solved only in one wingbeat period (from t=0 to t=T). Note that flow is unsteady and the solution of unsteady Navier-Stokes att=0 depends on the flow field before t=0, specifically, on the flow field at the time step before t=0 (i.e. at ; is the time step). Wu etal. showed that the periodic wake of a flapping wing is established in less than 3 cycles after the wing starts to beat, and that when the body oscillations were small, one could use the flow field established by the insect without body oscillation to provide the required initial flow field data [2].With these assumptions, the solution process of obtaining the periodic solutions is as follows:

1. Give the initial values of the state variables of the body, i.e. u, w, q and θ at (v, p, r and  are always zero at hover flight). In the first iteration, they are given as: u=0, w=0,q=0 andθ= ( is the mean body angle, known from measured data). In each of the later iterations, they are updated by the Newton’s method. For convenience, we write the initial values as , , and .

2. Start the wingbeat at , with the body fixed (u, w and q=0 andθ=), and at the same time start to compute the flow with the Navier-stokes equations. From to , u, w, q and θchange smoothly from zero to their initial values, , , and , respectively.

3. At , start the computation of coupling the equations of motion and the Navier-Stokes equations (the flow field at the time step before , which is needed for solving the Navier-Stokes equations at , is provided by the flow solution in Step 2), and obtain the solution in the period from to .

4. Compare the values of u, w, q and θat with those at , i.e. , , and . If they are respectively equal, the periodic solution is obtained and the calculation finished. If not, go to the next step.

5. Compute the Jacobian matrix in the Newton’s method and obtain the new initial state variables. With the new initial values, go to Step 2 for next iteration.

S3.4Correcting the wing kinematic parameters and the calculation flow diagram

In the above described solution process, the wings are supposed to flap with kinematic parameters of hover flight. Experimentally measured data are to be used. These are necessarily small errors in the measurement. As shown in Wu etal., when there are small errors in the wing kinematic parameters, a periodic solution not representing a hover flight, but a flight with some small speed is obtained [2]. For example, if the angle of attack used is a little larger than the real one, the computed periodic solution gives an upward flight of small speed. To overcome this problem, we do not use experimental data to prescribe αd, αu and , but determine them in the solution process, by requiring that they take values such that hovering flight is ensured. The reason for only three parameters being determined in this way is that there are three conditions to be met in hovering flight [2]: zero mean horizontal velocity, zero mean vertical velocity and zero mean pitch rate. The reason for choosing αd, αu and to be determined in the solution process is that experimental data for αd and αu have relatively large error [11] and aerodynamic forces and moments are very sensitive to the variations inαd, αu and .

The steps of obtaining the hovering flight periodic solution are summarized in a flow diagram shown in Fig. S4.

S4 Periodic solutionsof hovering flight

The equations of motion coupled with the Navier-Stokes equations are solved to obtain the periodic solution that represents the hovering flight. The numerical solutions, u(t), w(t), q(t) and θ(t),for the hovering model dronefly and model hawkmoth are shown in Fig. S5a and b, respectively (since hovering is symmetrical flight, v, p, r,  and  are zero). As discussed in Section 2 and Section S2, , and, are determined in the solution process using the conditions of equilibrium flight: mean horizontal and vertical velocities of the body center of mass and mean pitch rate of body are zero. It was found that for the hovering model hawkmoth, when=47.9 deg, =47.1 deg and =9.6 deg, periodic solutions with zero means in u, w and q were obtained, and for the dronefly, the corresponding , and values are 34.0, 34.5 and 3.0 deg., respectively.For the hawkmoth, the peak-to-peak value of the x component of the velocity of the center of mass is about 0.12U (Fig. S5a) and the peak-to-peak value of the pitch angle is about 5 deg. For the dronefly (Fig. S5b), the corresponding values are 0.025Uand 0.3 deg. The body oscillations in the hawkmoth are much larger than that of the dronefly. This is expected because the wingbeat frequency of the hawkmoth is lower and wing mass (relative to the body mass) is much larger than that of the dronefly.