We present the first integrative computational fluid dynamics (CFD) study of near- and far-field aerodynamics in insect hovering flight using a biology-inspired, dynamic flight simulator. This simulator, which has been built to encompass multiple mechanisms and principles related to insect flight, is capable of flying' an insect on the basis of realistic wing–body morphologies and kinematics. Our CFD study integrates near-and far-field wake dynamics and shows the detailed three-dimensional (3D)near- and far-field vortex flows: a horseshoe-shaped vortex is generated and wraps around the wing in the early down- and upstroke; subsequently, the horseshoe-shaped vortex grows into a doughnut-shaped vortex ring, with an intense jet-stream present in its core, forming the downwash; and eventually,the doughnut-shaped vortex rings of the wing pair break up into two circular vortex rings in the wake. The computed aerodynamic forces show reasonable agreement with experimental results in terms of both the mean force (vertical,horizontal and sideslip forces) and the time course over one stroke cycle(lift and drag forces). A large amount of lift force (approximately 62% of total lift force generated over a full wingbeat cycle) is generated during the upstroke, most likely due to the presence of intensive and stable,leading-edge vortices (LEVs) and wing tip vortices (TVs); and correspondingly,a much stronger downwash is observed compared to the downstroke. We also estimated hovering energetics based on the computed aerodynamic and inertial torques, and powers.

Flapping-flying insects employ unsteady aerodynamic mechanisms to keep them afloat, and there have been many studies on this topic(Ellington, 1984a; Ellington, 1984b; Ellington, 1984c; Ellington, 1984d; Ellington et al., 1996; van den Berg and Ellington,1997a; van den Berg and Ellington, 1997b; Dickinson et al., 1999; Liu,2002; Liu, 2005; Sane, 2003; Lehmann, 2004a; Lehmann, 2004b; Wang, 2005). A general conclusion from these studies is that insects obtain enough lift force to support their weight through the sophisticated vortices generated by the flapping wings. Among the many mechanisms involved in insect flight, delayed stall (Ellington et al.,1996), which is featured by prolonged attachment of a leading-edge vortex (LEV) on a wing, has been widely recognized as an important unsteady aerodynamic mechanism contributing to the enhancement of lift force generation in flapping-flying insects. However, the delayed stall has been deduced based exclusively on experiments on the hawkmoth Manduca sexta, which is one of the largest insects (wing span 5 cm). The question still remains as to whether the delayed stall presents and plays similar roles in smaller insects such as a fruit fly (wing span 0.2 cm). Birch et al.'s dynamically scaled robotic wing model experiments (Birch et al., 2004) offered some indirect evidence towards answering the question. Their experimental results suggested that the transport of vorticity from the leading edge to the wake, which permits prolonged vortex attachment and enhanced force production, might take different forms at different Reynolds numbers (Re; i.e. in insects of different sizes).

Besides the aforementioned LEV, other unsteady mechanisms might contribute to force production in insect flight. Dickinson et al. suggested that rotational circulation and wake capture increased aerodynamic force during the rotational phase of wing motion (Dickinson et al., 1999). The correct angular difference between two counter-lateral wings during dorsal stroke reversal (clap-and-fling') was found to increase total lift force by up to 17%, which indicated the possible role of wing–wing interaction in insect flight(Lehmann et al., 2005). The flow around a real hawkmoth in forward flight was visualized using digital particle image velocimetry (DPIV) (Bomphrey et al., 2005; Bomphrey et al.,2006). The visualization results described an apparent vortex structure across the insect thorax at the end of the downstroke, suggesting a potential effect of wing–body interaction on aerodynamic force generation.

Overall the above-mentioned studies, together with other relevant studies(Birch and Dickinson, 2001; Birch and Dickinson, 2003; Lehmann, 2002; Fry et al., 2003; Fry et al., 2005; Wang et al., 2004), have made significant contributions to understanding of the different aspects of the aerodynamic mechanisms involved in insect flight. However, most of the previous studies have been focused exclusively on the near-field flow and its correlation with aerodynamic force production. Therefore the challenging problem to quantify the near- and far-field flow structures and to correlate them to force-production still remains.

Both near- and far-field flows are strongly three-dimensional (3D) in space and unsteady in time. A reasonable investigation of the interaction between them entails obtaining sufficient information on the unsteady 3D flows. Obtaining such a large amount of information might be extremely difficult experimentally; at this point, computational fluid dynamics (CFD)-based simulation may provide a relatively more effective method. In fact, during the last decade CFD has been widely applied to studies concerning insect flight(Liu and Kawachi, 1998; Liu et al., 1998; Wang, 2000; Ramamurti and Sandberg, 2002; Sun and Tang, 2002a; Sun and Tang, 2002b; Wang et al., 2004; Miller and Peskin, 2004; Miller and Peskin, 2005; Liu, 2005). These simulation-based studies had varying degrees of success in dealing with the specific issues under consideration, but most of them were still either limited to two-dimensional (2D) computations, or lacking in quantitative evaluation of the effect of wing–wing and wing–body interactions.

In this study, we present the first integrative CFD study of the unsteady 3D near- and far-field vortex wake dynamics in a hovering fruit fly and their relation to lift force generation using a biology-inspired dynamic flight simulator (see Movies 1 and 2 in supplementary material; H.L., manuscript in preparation). To reproduce a hovering fruit fly on a computer using the simulator, both the wing–body morphological and kinematic models were built based faithfully on measurements from a real fruit fly, Drosophila melanogaster. The large amount of information on flow field and aerodynamic force offered by the computation enabled us to perform in-depth analysis of the correlation of vortex flow structure with force generation,the interaction between near- and far-field flows, and the effects of wing–wing and wing–body interactions. We also simulated a hovering hawkmoth and, by comparing the computed results for the two insects, we elucidate the marked dependence of the spanwise flow and the delayed stall on Re.

### A biology-inspired dynamic flight simulator

This study employs a biology-inspired, dynamic flight simulator, consisting of in-house software developed recently for quantitatively investigating the aerodynamics of a flapping-flying insect (see Movies 1 and 2 in supplementary material; H.L., manuscript in preparation). This simulator is able to fly' an insect with realistic body–wing morphologies and flapping-wing and body kinematics, and to evaluate unsteady aerodynamics including detailed vortex flow fields and flying energetics involving aerodynamic and inertial torques and powers. Details of the simulator can be found in the Appendix.

### Morphological and kinematic model

We constructed a realistic wing–body morphological model(Fig. 1) for the fruit fly Drosophila melanogaster, body length 2.78 mm and wing length 2.39 mm. We traced the borderlines of the wings and the body from pictures of a fruit fly taken in two perpendicular views, and then reconstructed the shape of the fruit fly on a computer based on the borderlines, by assuming that the cross sections of both the body and the wings are in an ellipse shape. We assumed a uniform wing thickness (1.2% of the mean chord) for the two wings, which resembles the wing geometry of a real fruit fly. Finally, curve smoothing was performed at the leading and trailing edges, and at the wing tip. Note that to avoid the attachment of the wing on the body surface we added a virtual portion of the wing length (approximately /32, where is the mean wing chord length) at the wing base, which could largely improve the numerical convergence but seldom affect the results in the hovering flight. More details can be found in the Appendix (morphological modeling).

Fig. 1.

A morphological model of a fruit fly Drosophila melanogaster. (A)A fruit fly with a computational model superimposed on the right half(http://www.tmd.ac.jp/artsci/biol/textlife/fruitfly.jpg). The fruit fly has a body length of 2.78 mm, a wing length of 2.39 mm (mean wing chord length =0.78 mm), and an aspect ratio of 3.06. (B) A multi-block grid system of the two wings and body of the fruit fly (wing, 45×45×31; body, 45×47×95) with a distance between the body surface and the outer boundary of 20.

Fig. 1.

A morphological model of a fruit fly Drosophila melanogaster. (A)A fruit fly with a computational model superimposed on the right half(http://www.tmd.ac.jp/artsci/biol/textlife/fruitfly.jpg). The fruit fly has a body length of 2.78 mm, a wing length of 2.39 mm (mean wing chord length =0.78 mm), and an aspect ratio of 3.06. (B) A multi-block grid system of the two wings and body of the fruit fly (wing, 45×45×31; body, 45×47×95) with a distance between the body surface and the outer boundary of 20.

The wing–body kinematic model was established based on the measurements of a hovering fruit fly, Drosophila melanogaster(Fry et al., 2005)(Fig. 2). Reynolds number was defined by Re=c̄Uref/ν,where ν is the kinematic viscosity of air (1.5×10–5m2 s–1). According to the measured data(=0.78 mm, R=2.39 mm,Φ=2.44 rad (140°), f=218 s–1) of the fruit fly (see List of symbols and abbreviations for definitions) in this study, Re was estimated to be 134 with a reduced frequency K=0.212. More details can be found in the Appendix (kinematic modeling).

### Far-field flow: vortex wake structures and downwash

The absolute iso-vorticity surfaces around the hovering fruit fly and the velocity vectors at the plane perpendicular to the stroke plane (see Fig. 2B) are illustrated in Fig. 3Ai,ii,Bi,ii,Ci,ii,Di) at four typical moments over a flapping cycle (see Fig. 2Ca,b,e,g). Note that the absolute iso-vorticity surfaces around the hovering fruit fly without the body are also illustrated in Fig. 3Aiv,v,Biv,v,Civ,v,Diii to evaluate the effect of the body. In addition, to understand the relationship between the flow fields and the intensity of the downwash, the downward speed contours and the velocity vectors in the stroke plane are depicted in Fig. 3Aiii,Biii,Ciii,Dii and Fig. 4.

Fig. 2.

Schematic diagram of the computational system of a fruit fly Drosophila melanogaster. (A) The local wingbase-fixed (x, y, z) and the global earth-fixed (X, Y, Z)coordinate systems. The origin O′ of the wingbase-fixed coordinate system lies at the wing base, with the x-axis normal to the stroke plane [the yz plane as defined by Ellington(Ellington, 1984b)], the y-axis vertical to the body axis and z-direction parallel to the stroke plane. (B) The wing kinematics are described by the positional angle φ, the feathering angle (angle of attack of the wing) α, the elevation angle θ, and the stroke plane angle β; the link to the earth-fixed frame of reference comes through the body angle χ. We assume a body angle χ of 45° and a stroke plane angle β of 0°(Fry et al., 2005). (C)Instantaneous positional angle φ, feathering angle α, and elevation angle θ of the fruit fly wing over one complete flapping cycle. Green solid, orange broken and blue dash-dot lines represent the positional angleφ, the feathering angle α and the elevation angle θ,respectively. Red points a–g: (a) mid pronation, (b) early downstroke,(c) mid downstroke, (d) late downstroke, (e) early upstroke, (f) mid upstroke and (g) late upstroke. T, dimensionless period of one flapping cycle.

Fig. 2.

Schematic diagram of the computational system of a fruit fly Drosophila melanogaster. (A) The local wingbase-fixed (x, y, z) and the global earth-fixed (X, Y, Z)coordinate systems. The origin O′ of the wingbase-fixed coordinate system lies at the wing base, with the x-axis normal to the stroke plane [the yz plane as defined by Ellington(Ellington, 1984b)], the y-axis vertical to the body axis and z-direction parallel to the stroke plane. (B) The wing kinematics are described by the positional angle φ, the feathering angle (angle of attack of the wing) α, the elevation angle θ, and the stroke plane angle β; the link to the earth-fixed frame of reference comes through the body angle χ. We assume a body angle χ of 45° and a stroke plane angle β of 0°(Fry et al., 2005). (C)Instantaneous positional angle φ, feathering angle α, and elevation angle θ of the fruit fly wing over one complete flapping cycle. Green solid, orange broken and blue dash-dot lines represent the positional angleφ, the feathering angle α and the elevation angle θ,respectively. Red points a–g: (a) mid pronation, (b) early downstroke,(c) mid downstroke, (d) late downstroke, (e) early upstroke, (f) mid upstroke and (g) late upstroke. T, dimensionless period of one flapping cycle.

Fig. 3.

Far-field flow structures around a hovering fruit fly. Absolute iso-vorticity surfaces and velocity vectors with a body (Ai,ii,Bi,ii,Ci,ii,Di)and without a body (Aiv,v,Biv,v,Civ,v,Diii) at four instances of (a), (b), (e)and (g), respectively, as illustrated in Fig. 2C. (Aiii,Biii,Ciii,Dii)The corresponding velocity vectors and contours of vertical velocity components in the stroke plane (red broken circle).

Fig. 3.

Far-field flow structures around a hovering fruit fly. Absolute iso-vorticity surfaces and velocity vectors with a body (Ai,ii,Bi,ii,Ci,ii,Di)and without a body (Aiv,v,Biv,v,Civ,v,Diii) at four instances of (a), (b), (e)and (g), respectively, as illustrated in Fig. 2C. (Aiii,Biii,Ciii,Dii)The corresponding velocity vectors and contours of vertical velocity components in the stroke plane (red broken circle).

As seen in Fig. 3Ai,ii, a ring-shaped vortex wake structure, which results from the vortices generated and detached from the wings and body during the preceding upstroke, is evident during pronation (see Fig. 2Ca). This wake structure is the most prominent feature over a flapping cycle. Since the starting vortex (the detached TEV) balances the circulation of the bound vortex, the air within the vortex ring is given downward momentum, the downwash, which is observed flowing downward through the center of the vortex ring as a jet(Fig. 3Ai,ii).

In the first half of the downstroke (see Fig. 2Cb, Fig. 3Bi,ii), a horseshoe-shaped vortex is observed, which wraps around the wing and comprises a leading-edge vortex (LEV), a wing tip vortex (TV) and a trailing-edge vortex(TEV) (or starting vortex). A similar horseshoe-shaped vortex was also observed in an experimental DPIV study of an impulsively started, dynamically scaled, flapping wing (Poelma et al.,2006). During the downstroke, the LEV and the TEV grow steadily in size and expand towards the wing base. Eventually the TEV begins to detach from the trailing edge and joins the TV. Subsequently, the LEV, the TV and the shed TEV in toto form a doughnut-shaped vortex ring for one wing, and hence a pair of vortex rings for the wing-pair. During most of the downstroke,the doughnut-shaped vortex ring pair has an intense, downward jet-flow through the doughnut' hole, which forms the downwash during the downstroke(Fig. 3Bi,ii,iii). In the second half of the downstroke, the TV gradually enlarges, subsequently the LEV and the TV weaken and detach from the upper surface, and the doughnut-shaped vortex rings of the wing pair break up into two downward circular vortex rings, forming the far-field wake below the fly.

During the upstroke (see Fig. 2Ce,g), as illustrated in Fig. 3Ci,ii,Di, the wing also generates a ring-shaped vortex wake structure that resembles, but is more structured than, the wake during the downstroke. The exact shape of the wake depends on how well the vortices shed from the wing merge with those vortices still attached to the wing(Fig. 3Ai,ii,Ci,ii,Di). The upstroke's downwash also resembles that of the downstroke.

Fig. 4.

Downward velocity maps at two cross sections around a hovering fruit fly at mid down- and upstroke. (A) Schematic diagram of the two cross sections around the hovering fruit fly, showing the position of velocity maps at the two cross sections at mid down- and upstroke (B–E). (B,C) Velocity maps at the plane parallel to the X-axis (a–a′); (D,E) Velocity maps at the plane parallel to the Y-axis (b–b′). Color of vectors indicates magnitude of downward velocity.

Fig. 4.

Downward velocity maps at two cross sections around a hovering fruit fly at mid down- and upstroke. (A) Schematic diagram of the two cross sections around the hovering fruit fly, showing the position of velocity maps at the two cross sections at mid down- and upstroke (B–E). (B,C) Velocity maps at the plane parallel to the X-axis (a–a′); (D,E) Velocity maps at the plane parallel to the Y-axis (b–b′). Color of vectors indicates magnitude of downward velocity.

Fig. 4 illustrates that the downward flow (downwash) velocity maps at two cross sections around the hovering fruit fly during mid down- and up-stroke: one is parallel to the X-axis and the other to the Y-axis. The computed results demonstrate that, throughout a flapping cycle, despite the existence of downward flows beneath the body, the strong downwash is only present in a limited area swept by the flapping wing, which shows a maximum speed of approximately 0.8Uref over one flapping cycle and is largely attenuated within 1.5 times of the body length of the fruit fly far behind the wing (see Fig. 3Aiii,Biii,Ciii,Diii and Fig. 4).

Fig. 5.

Streamlines around the wing at mid down- and upstroke. Streamlines are color-coded to indicate the absolute flow speed. (A) The LEV and the TV at mid downstroke (see Fig. 2Cc). (B)The LEV and the TV at mid upstroke (see Fig. 2Cf). LEV, leading-edge vortex; TV, wing tip vortex (for further explanation, see text).

Fig. 5.

Streamlines around the wing at mid down- and upstroke. Streamlines are color-coded to indicate the absolute flow speed. (A) The LEV and the TV at mid downstroke (see Fig. 2Cc). (B)The LEV and the TV at mid upstroke (see Fig. 2Cf). LEV, leading-edge vortex; TV, wing tip vortex (for further explanation, see text).

### Near-field flow: the leading-edge and the wing tip vortices

#### Downstroke

During the early downstroke, the flapping wings generate a horseshoe-shaped vortex, which initially is limited to the wing tip region. This vortex comprises a LEV and a TEV connecting to a TV. In the first half of the downstroke, the LEV and the TEV first develop and stretch from the wing tip towards the wing base. Both vortices mostly show a 2D structure with a very weak axial flow. Subsequently the TEV (or a starting vortex) detaches from the trailing edge but remains joined to the TV, forming the doughnut-shaped vortex ring. Note that the LEV begins to grow from the wing tip towards the wing base, possibly as a result of the onset of the TV. As the positional angle of wing approaches zero, the LEV continues to grow while connecting to the TV(Fig. 5A). In the second half of the downstroke, the LEV initially remains attached to the wing surface,with no evidence of breaking down or shedding. At the end of the downstroke and the start of supination, the LEV breaks down at approximately 70–80%of wing length from the wing base. As the wing rotates and decelerates, the LEV is transformed from a 2D into a 3D structure(Fig. 6A). At the same time,the TV becomes increasingly unstable and gradually grows towards the wing base and eventually it covers almost half the wing.

#### Supination

Supination of the wing is completed when its positional angle changes from–36.7° to –55.0°, and the angle of attack from–32.0° to 50.2°. During early supination, the downstroke LEV and TV are shed due to the rotation of the wing and are washed downward. Moreover,because the wing translational speed becomes very slow(Fig. 6A), a stopping vortex forms around the leading edge. Immediately after the wing reversal, a horseshoe-shaped vortex forms and wraps around the wing. When the wing translates upwards, rotates and increases its angle of attack, the TEV is shed and washed downward, while the LEV develops and grows gradually along the leading edge from the wing tip to the wing base. In the second half of supination, the wing intercepts its own wake.

#### Upstroke

At the beginning of the upstroke immediately after the wing reversal, a horseshoe-shaped vortex is visible wrapping around the wing and consisting of a TV, a LEV and a TEV. Similar to the downstroke, the TEV subsequently detaches from the trailing edge but connects to the TV; and eventually a doughnut-shaped vortex ring is observed. As the LEV continues to grow, a spanwise flow becomes indiscernible in its core, but the LEV is mostly in a 2D structure (Fig. 7Bii,iii). The spanwise flow is not, however, observed in the core of the LEV but along the rearmost half of the wing. As the positional angle of the wing approaches zero at mid upstroke, a large LEV and TV become visible(Fig. 5B) and a corresponding negative pressure region is observed on the upper surface of the wing. The upstroke LEV continues to grow without breaking down or shedding throughout the translational phase of the upstroke. When the wing starts pronation, the LEV and the TV become unstable and a stopping vortex is generated, and together form a complicated 3D vortex structure(Fig. 6B).

#### Pronation

During early pronation, the LEV and the TV are unstable but remain attached to the wing. At approximately 70% of the wing length from the wing base on the upper wing surface, the TV pulls the LEV upward away from the wing surface(Fig. 6B). The doughnut-shaped vortex ring detaches from the wing and merges with the downward wake. During late pronation, the motion of fruit fly wing resembles the fling motion employed by many small insects (Weis-Fogh,1973), which generates two vortices between the left and right wings. We can observe this process in our simulation in the form of a TV and the initial horseshoe-shaped vortex. After the pronation, while the wing translates downwards and rotates, and its angle of attack increases, the starting vortex is shed off from the trailing edge.

Fig. 6.

Streamlines around a fruit fly at early supination and pronation. Streamlines around a hovering fruit fly at early supination (A) and at early pronation (B). A color map represents pressure contours on the wing and body surfaces: red is high and blue is low; streamlines are released at the leading edge from the wing surfaces, color of streamlines indicates the absolute flow speed. LEV, leading-edge vortex; TV, wing tip vortex; SV, stopping vortex (for further explanation, see text).

Fig. 6.

Streamlines around a fruit fly at early supination and pronation. Streamlines around a hovering fruit fly at early supination (A) and at early pronation (B). A color map represents pressure contours on the wing and body surfaces: red is high and blue is low; streamlines are released at the leading edge from the wing surfaces, color of streamlines indicates the absolute flow speed. LEV, leading-edge vortex; TV, wing tip vortex; SV, stopping vortex (for further explanation, see text).

#### Re dependence on the role of flapping wing motion in the LEV formation

To identify the Re dependence on the role of translational and rotational motion of the flapping wing in the formation of the LEV, we compared the computed velocity vector maps and the spanwise flow contours at the plane transecting the wing at 60% span at Re=6300 (hawkmoth) and Re=134 (fruit fly) (Fig. 8Aii,Bii). The computed results suggest that the translational motion of the wings (rotation of positional angle of the wings) dominates the development of the LEV, while wing rotation about the leading edge (angle of attack of wings) benefits the stability of the LEV. The role of wing rotation is more evident at the low Re (134) than at the high Re(6300). At the high Re (6300), an axial flow at the core of the LEV is much more pronounced, and the LEV forms a 3D structure near the leading edge; by contrast, at the low Re (134), only a weak axial flow is detected, indicating that the vortical structure near the leading edge might be mostly in a 2D structure. Yet at such a low Re, a remarkable spanwise flow is predicted along the rearmost half of the wing at mid down-and upstroke (Fig. 7Aii,iii,Bii,iii).

### Evaluation of hovering energetics

#### Forces

The mean aerodynamic forces are quantified by the following equations (Eqn 1, 2, 3) using the mean force coefficients over a flapping cycle:
$\ \mathrm{Vertical}{\ }\mathrm{force}{\ }(\mathrm{lift})=0.5{\rho}U_{\mathrm{ref}}^{2}S_{\mathrm{w}}C_{\mathrm{L}},$
(1)
$\ \mathrm{Horizontal}{\ }\mathrm{force}{\ }(\mathrm{drag}{\ }\mathrm{and}{\ }\mathrm{thrust})=0.5{\rho}U_{\mathrm{ref}}^{2}S_{\mathrm{w}}C_{\mathrm{D}},$
(2)
$\ \mathrm{Sideslip}{\ }\mathrm{force}=0.5{\rho}U_{\mathrm{ref}}^{2}S_{\mathrm{w}}C_{\mathrm{S}},$
(3)
where ρ is the density of the air (1.23 kg m–3), Uref is the reference velocity, Sw is the planform area of the wing and CL,CD, CS are the three non-dimensional mean force coefficients of the vertical, horizontal and sideslip forces,respectively. According to Eqn 1, 2, 3, with the values assigned to the parameters defined in the CFD simulation (Uref=2.54 m s–1), the mean lift force is calculated to be 9.60×10–6 N, which exceeds the weight of the fruit fly(9.41×10–6 N) by 2%; the horizontal and sideslip forces computed are less than 4% of the vertical force. These force predictions agree very well with the situation expected for hovering flight and therefore indirectly validate our simulations.

While these mean force coefficients provide a validation of our time-averaged results, the time courses of computed vertical (lift) and horizontal (drag and thrust) forces over one flapping cycle can be validated against the experimental data (Fry et al.,2005) (Fig. 9). This comparison confirms the match between time-averaged computational and experimental data. However, there are discrepancies between the time courses of the lift force. The main discrepancy occurs in the phase of the periodic force signature, especially in the early down- and upstroke. These phase-shift differences might be caused by different locations of the wing's axis of rotation which, in the computations, is assumed to lie at a quarter chord length from the leading edge. For horizontal force (drag and thrust), we obtain satisfactory agreements between the computation and the experiment, for both the mean and instantaneous magnitudes. The lift force peaks twice, at mid downstroke and mid upstroke (Fig. 9A). The two peaks are most likely due to the existence of a large negative pressure area on the upper wing surface, induced by the LEV and the TV (Fig. 5). The drag force is produced exclusively during the downstroke while the thrust force is produced only during the upstroke (Fig. 9B).

#### Torques

The time courses of the computed aerodynamic and inertial torques of the wing are plotted in Fig. 10A–C for rolling, yawing and pitching, respectively. The maximums of the aerodynamic rolling torque (ART) and the aerodynamic yawing torque (AYT) are approximately 3.0×10–10 Nm, which is only 1% of that of the aerodynamic pitching torque (APT)(2.0×10–8 Nm). Obviously, the APT dominates during most of the flapping cycle, except in the early down- and upstroke, where the APT and the inertial pitching torque (IPT) become comparable. Hence, it is the APT that controls the pitching motion and the IPT may work as a secondary controller that adjusts the total pitching torques (TPT). The IPT might also play an important role in insect steering maneuvers(Wang et al., 2005).

Fig. 7.

Velocity vectors and contours of spanwise flow velocity at mid down- and upstroke. A wing–body computational model of a fruit fly (Ai) at mid downstroke, and (Bi) at mid upstroke. Velocity vectors at two cross sections at 50% and 60% of the wing length, (Aii) and (Aiii) at mid downstroke, and(Bii) and (Biii) at mid upstroke. A color map denotes the spanwise flow velocity; a negative magnitude of the spanwise velocity points to a flow from the wing base to the wing tip, and vice versa. L.E., leading edge;T.E., trailing edge.

Fig. 7.

Velocity vectors and contours of spanwise flow velocity at mid down- and upstroke. A wing–body computational model of a fruit fly (Ai) at mid downstroke, and (Bi) at mid upstroke. Velocity vectors at two cross sections at 50% and 60% of the wing length, (Aii) and (Aiii) at mid downstroke, and(Bii) and (Biii) at mid upstroke. A color map denotes the spanwise flow velocity; a negative magnitude of the spanwise velocity points to a flow from the wing base to the wing tip, and vice versa. L.E., leading edge;T.E., trailing edge.

#### Powers

Based on computed instantaneous aerodynamic forces and wing velocities, we calculate muscle-mass-specific inertial, aerodynamic and mechanical powers(Fig. 11). For comparison, the data measured by Fry et al. (Fry et al.,2005) are also given in Fig. 11. The muscle-mass-specific inertial power Piner denotes the power consumed to accelerate the mass of the wing (Eqn A14). The time course of the computed inertial power is in satisfactory qualitative, but not quantitative agreement with the experimental results(Fry et al., 2005), which is thought to be due to the difference in the wing shape between the computational and experimental models (Fry et al., 2005). By integrating inertial power over a flapping cycle, we find that an average 56.9 W kg–1 is needed to accelerate the wing. Note that wing deceleration is assumed to accrue no cost and, at the same time, there is no elastic power storage. The muscle-mass-specific aerodynamic power Paero is the power necessary to overcome air resistance (Eqn A15). The comparison between computations and measurements shows good agreement except for a slight phase lag between the computational and the experimental data (Fry et al.,2005). The computed mean aerodynamic power (89.3 W kg–1) is very close to the experimental result of Fry et al.(Fry et al., 2005). Maximum aerodynamic power is reached at each mid stroke when the flapping wing speed and the aerodynamic forces approach maximum values. The muscle-mass-specific total mechanical power Ptotal is the power required to move the wing. We simply calculate Ptotal using Eqn A16 as the sum of the muscle-mass-specific aerodynamic and inertial powers. It can be seen that the curve of Ptotal turns to be negative in the second half of the upstroke because more power is required to overcome aerodynamic forces when the wing decelerates (Fig. 11C).

Fig. 8.

Comparison of near-field flow between a fruit fly and a hawkmoth at mid downstroke. A wing-body computational model of the hawkmoth (Re=6300, Uref=5.05 m s–1, cm=1.83 cm) (Ai) and the fruit fly (Re=134, Uref=2.54 m s–1, cm=0.78 mm) (Bi), with the LEVs visualized by streamlines; (Aii, Bii) the corresponding velocity vectors at the cross sections at 60% of the wing length. A color map denotes the spanwise flow velocity; a negative magnitude of the spanwise velocity points to a flow from the wing base to the wing tip, and vice versa. L.E., leading edge;T.E., trailing edge.

Fig. 8.

Comparison of near-field flow between a fruit fly and a hawkmoth at mid downstroke. A wing-body computational model of the hawkmoth (Re=6300, Uref=5.05 m s–1, cm=1.83 cm) (Ai) and the fruit fly (Re=134, Uref=2.54 m s–1, cm=0.78 mm) (Bi), with the LEVs visualized by streamlines; (Aii, Bii) the corresponding velocity vectors at the cross sections at 60% of the wing length. A color map denotes the spanwise flow velocity; a negative magnitude of the spanwise velocity points to a flow from the wing base to the wing tip, and vice versa. L.E., leading edge;T.E., trailing edge.

### The leading-edge vortex

Our CFD analysis indicates that the LEV of fruit flies differs in shape from that of hawkmoths. The main difference is in the position on the wing where the LEV is attached and in the intensity of the axial flow in the LEV's core. Birch and Dickinson's experimental results(Birch and Dickinson, 2001)indicate that the LEV of fruit flies exhibit a rather stable vortex structure without separation during most of the down- and upstroke, and show no evidence of axial flow in the vortex core. These observations are in stark contrast to those obtained with hawkmoths, whose LEV detaches from the wing surface at approximately 75% of the wing length and whose LEV shows a strong axial flow in the core (van den Berg and Ellington,1997a). The computed results confirm the LEV features described by Birch and Dickinson (Birch and Dickinson,2001), who speculated that the difference between hawkmoths and fruit flies can be attributed to effects of size or Re (100–250 for fruit flies, >6000 for hawkmoths). Birch and Dickinson further suggested that considering that a typical adult insect is closer in size to a fruit fly and that therefore their observations of flow patterns are more likely to be typical, the insects are more likely to prolong the attachment of the LEV due to the attenuating effect of the downwash. They explained the absence of the axial flow in the LEV core of fruit flies by pointing out that the pressure gradients within the vortex core might be too small to drive a substantial axial flow on the smaller wing(Birch and Dickinson,2001).

Fig. 9.

Time courses of the vertical (lift; A) and the horizontal (drag and thrust;B) forces over a flapping cycle. Blue, red and yellow lines represent the measurements of the upper (Exp_u), average (Exp_a) and lower (Exp_l) values obtained by Fry et al. (Fry et al.,2005), respectively; broken line is the computed result (Com_a). T, dimensionless period of one flapping cycle.

Fig. 9.

Time courses of the vertical (lift; A) and the horizontal (drag and thrust;B) forces over a flapping cycle. Blue, red and yellow lines represent the measurements of the upper (Exp_u), average (Exp_a) and lower (Exp_l) values obtained by Fry et al. (Fry et al.,2005), respectively; broken line is the computed result (Com_a). T, dimensionless period of one flapping cycle.

To test this hypothesis, Fig. 12 illustrates the pressure gradient contours on the wing of a fruit fly and a hawkmoth. Obviously, the computed pressure gradients on the wing of the hawkmoth are much larger than those of the fruit fly. This observation suggests that the pressure gradient very likely causes the strong axial flow at the LEV core and enhances the LEV's stability as long as Reynolds numbers are high enough – in the order of several thousands(Ellington et al., 1996; van den Berg and Ellington,1997a; Liu and Kawachi,1998; Liu et al.,1998). At a low Re of 100–250, e.g. for the fruit fly, the flapping wing cannot create a LEV strong enough to generate a steep pressure gradient at the vortex core. Nevertheless, our results indicate that a fruit fly can produce a stable LEV during most of the down- and upstroke. While the LEV of the hovering hawkmoth breaks down roughly at mid downstroke(van den Berg and Ellington,1997a; Liu et al.,1998), the LEV of the hovering fruit fly remains attached to the leading edge and grows stably throughout the downstroke, eventually breaking down during the subsequent supination (Fig. 6A). At low Re values, the LEV on a flapping or rotary wing may translate over a longer distance before growing sufficiently to break down. Still, it is valid to ask which forces or mechanisms are responsible for enhancing LEV stability. As observed in other studies on the flapping wing of the fruit fly (Birch and Dickinson,2001; Birch et al.,2004), we also find a fairly strong spanwise flow outside that of the LEV, but in the rearmost half of the wing, during most of the down- and upstroke. While we show that the pressure gradient is not the cause(Fig. 12), centrifugal and Coriolis forces just might be sufficient to create the spanwise flow outside the LEV's core in the rearmost half of the wing that is indicted by our CFD analysis.

Fig. 10.

Time courses of aerodynamic and inertial torques of the wing over a flapping cycle. Aerodynamic and inertial torques of the wing as shown in Fig. 2B. (A) Rolling torques:aerodynamic rolling torque (blue; Tr_aero_t), inertial rolling torque (red;Tr_iner_t) and total rolling torque (green; Tr_total). (B) Yawing torques:aerodynamic yawing torque (blue; Ty_aero_t), inertial yawing torque (red;Ty_iner_t) and total yawing torque (green; Ty_total). (C) Pitching torques:aerodynamic pitching torque (blue; Tp_aero_t), inertial pitching torque (red;Tp_iner_t) and total pitching torque (green; Tp_total). T,dimensionless period of one flapping cycle.

Fig. 10.

Time courses of aerodynamic and inertial torques of the wing over a flapping cycle. Aerodynamic and inertial torques of the wing as shown in Fig. 2B. (A) Rolling torques:aerodynamic rolling torque (blue; Tr_aero_t), inertial rolling torque (red;Tr_iner_t) and total rolling torque (green; Tr_total). (B) Yawing torques:aerodynamic yawing torque (blue; Ty_aero_t), inertial yawing torque (red;Ty_iner_t) and total yawing torque (green; Ty_total). (C) Pitching torques:aerodynamic pitching torque (blue; Tp_aero_t), inertial pitching torque (red;Tp_iner_t) and total pitching torque (green; Tp_total). T,dimensionless period of one flapping cycle.

The computed results of our analysis indicate that the lift-boosting mechanism of the LEV observed in hovering hawkmoths(Ellington et al., 1996) is also likely to enhance lift force in the hovering fruit fly, because the delayed stall also strengthens the fruit fly LEV and hence augments lift force. However, we find significant discrepancies in lift force production during the down- and upstrokes: the fruit fly produces 62% of the total lift force during the upstroke and 38% during the downstroke. In contrast,hawkmoths produce approximately 30–40% of the total lift force during the upstroke and 60–70% during the downstroke(Liu and Kawachi, 1998; Liu et al., 1998; Aono and Liu, 2006), and hummingbirds produce approximately 25% of the total lift force during the upstroke and 75% during the downstroke(Warrick et al., 2005). This implies that the size differences in wing and body kinematics play a key role in the aerodynamic force generation of hovering insects and birds.

Fig. 11.

Time courses of muscle-mass-specific powers over a flapping cycle. (A)Inertial powers: solid line (Exp) and broken line (Com) express the experimental data (Fry et al.,2005) and the computational data, respectively. (B) Aerodynamic powers: upper (blue; Exp_u), average (red; Exp_a) and lower (yellow; Exp_l)data from Fry et al. (Fry et al.,2005; green broken line is the computed data (Com_a). (C) Total mechanical powers: blue solid line (Exp) and orange broken line (Com) express the experimental data (Fry et al.,2005) and the computational data. T, dimensionless period of one flapping cycle.

Fig. 11.

Time courses of muscle-mass-specific powers over a flapping cycle. (A)Inertial powers: solid line (Exp) and broken line (Com) express the experimental data (Fry et al.,2005) and the computational data, respectively. (B) Aerodynamic powers: upper (blue; Exp_u), average (red; Exp_a) and lower (yellow; Exp_l)data from Fry et al. (Fry et al.,2005; green broken line is the computed data (Com_a). (C) Total mechanical powers: blue solid line (Exp) and orange broken line (Com) express the experimental data (Fry et al.,2005) and the computational data. T, dimensionless period of one flapping cycle.

### The vortex structure and the downwash

During most of the down- and upstroke, as shown in Fig. 3, the doughnut-shaped vortex rings of the wing pair are observed, which eventually detach from the wing and body during the subsequent supination and pronation, breaking up into two circular vortex rings, with strong downward flow through the core (hole)of the vortex ring. These phenomena, of a vortex ring and its breaking-up,were also predicted in an experimental study based an analysis of the vortex wake structures around a robotic hawkmoth model(van den Berg and Ellington,1997b). Note that the root vortex as observed by van den Berg and Ellington is absent here at the wing base during most of the downstroke, very likely because of the existence of the insect body. van den Berg and Ellington also observed a vortex wake ring with an intensive downwash through its centre in the wing wake during the downstroke, and called it a dumbbell-shaped vortex structure (van den Berg and Ellington,1997b), but they did not quantify the formation, development and break-up of the vortex ring because of technical limitations with their smoke-rake flow visualization. Nevertheless, they predicted that the dumbbell-shaped vortex structure should eventually break up into two single circular vortex rings rather than merge together into a single one. The computed vortex wake structure during the downstroke of the fruit fly analyzed here also differs from the descriptions of the shed vortex wake observed in previous studies of insect flight(Ellington, 1984d; Grodnitsky and Morozov, 1993; Dickinson and Götz, 1996; Willmott et al., 1997; Birch and Dickinson, 2001).

In our fruit fly and hawkmoth simulations, we also find this vortex ring pair (we call them a doughnut-shaped vortex ring pair). It has an intense downwash or jet-stream at its core, and eventually breaks up into two single circular vortex rings during the subsequent supination and pronation. A follow-up study on the details of a lattice wake structure and dynamics of the downwash will be conducted in the near future. Moreover, our results also reveal that the vortex wake structure is more compact than the wakes of other insects and mechanical insect models reported by many previous studies(Ellington, 1984d; Grodnitsky and Morozov, 1993; Dickinson and Götz, 1996). However, whether the compact wake structure could benefit the flight efficiency is not clear in the scope of this study.

Although the ring-shaped vortex wake structure is observed over a complete flapping cycle, the wake formed during the upstroke is more structured and stronger than the wake formed during the downstroke. Corresponding to the intensity of the vortex, the strongest downwash is present on the far-body side of the ring-shaped vortex core (Figs 3 and 4). The relationship between the downwash generated by the flapping wings and instantaneous aerodynamic force production is illustrated in Figs 4 and 9. At mid down- and upstroke,high-lift force is produced by the flapping wings (see Fig. 9), and then the remarkable downwash in the stroke plane is also observed. In attention to the intensity of the downwash, the downwash at mid upstroke is stronger than that at mid downstroke. However, more detailed quantitative study is required to understand how the downwash is linked to the aerodynamic force generation in flapping-flying insects.

Fig. 12.

Pressure gradient contours on the wings of a fruit fly (A) and a hawkmoth(B) at mid downstroke. The white arrows represent the core and direction of the spanwise pressure gradient on the wing surfaces.

Fig. 12.

Pressure gradient contours on the wings of a fruit fly (A) and a hawkmoth(B) at mid downstroke. The white arrows represent the core and direction of the spanwise pressure gradient on the wing surfaces.

For the fruit fly, the very weak axial flow in the LEV core implies that the LEV does not feed significant amounts of vorticity into the TV, as is the case for the hawkmoth. In other words, the LEV discharges little energy into the TV and hence into the total vortex wake structure. However, adding flow energy locally by prolonging the LEV's attachment to the wing does not necessarily increase global flight efficiency because of wake effects, such as an increased downwash. Based on the results of robotic experiments, Birch and Dickinson suggested that the downwash had a potent inhibitory effect on the magnitude of force production (Birch and Dickinson, 2001). This is confirmed by our computed results as shown in Fig. 13, where time course of the lift force is plotted for three complete stroke cycles. The mean lift force generated during the second cycle is reduced by 22.5% compared with the first cycle, but almost no difference is observed between the second and the third cycle. Considering that insects begin to flap in static air, the drop in lift force between the first and second stroke cycle can be attributed to the inhibitory effect of the downwash, which is fully established at the end of the first cycle.

Fig. 13.

Time course of vertical (lift) force during the first three flapping cycles. The mean lift force generated during the second cycle is reduced significantly by 22.5% compared with the first cycle, but almost no difference is observed between the second and third cycles.

Fig. 13.

Time course of vertical (lift) force during the first three flapping cycles. The mean lift force generated during the second cycle is reduced significantly by 22.5% compared with the first cycle, but almost no difference is observed between the second and third cycles.

### The effect of the body on the wake dynamics of hovering flight

Fig. 3Aiv,v,Biv,v,Civ,v,Diiishows that removing the body changes the vortex wake structure only slightly. Specifically, the trace of the shed TEV seems to expand slightly now that the barrier formed by the body is removed. This slight change in the flow field is also confirmed by its effect, albeit small, on aerodynamic forces shown in Fig. 14. Albeit fluctuations of up to 100% in instantaneous lift and sideslip force coefficients are observed at particular moments, the magnitudes of the mean lift and sideslip force coefficients with versus without the presence of the body differ by less than 2%. In short, the effect of the body on hovering aerodynamics is very small and hence it can be neglected. Nevertheless, in other flight modes, such as during forward flight and quick turn, the body-effect-related instantaneous forces or the inertial torques of the body may play important roles in the flight control.

Fig. 14.

Comparison of aerodynamic force coefficients between a fruit fly with and without the body over a flapping cycle. While the mean aerodynamic force coefficients over one flapping cycle are similar, time courses of aerodynamic force coefficients differ between the simulation with body' and without body'. Solid lines, lift coefficients with the body (orange) and without the body (brown); broken lines, drag coefficients with the body (dark blue) and without the body (sky blue); dotted lines, sideslip force coefficients with the body (pink) and without the body (purple). T, dimensionless period of one flapping cycle.

Fig. 14.

Comparison of aerodynamic force coefficients between a fruit fly with and without the body over a flapping cycle. While the mean aerodynamic force coefficients over one flapping cycle are similar, time courses of aerodynamic force coefficients differ between the simulation with body' and without body'. Solid lines, lift coefficients with the body (orange) and without the body (brown); broken lines, drag coefficients with the body (dark blue) and without the body (sky blue); dotted lines, sideslip force coefficients with the body (pink) and without the body (purple). T, dimensionless period of one flapping cycle.

### A biology-inspired dynamic flight simulator

In the following we give a brief description of the methodology of the simulator with a specific focus on three relevant modeling steps: (1)morphological modeling, (2) kinematic modeling, and (3) multi-block- and overset-grid-based computational fluid dynamic (CFD) modeling.

### Morphological modeling

We obtain a morphological model of an object in four steps. First, we image the object digitally; secondly, we segment the image to extract the object's shape as a wire frame and/or skeleton model; thirdly, we smooth and curve/surface-fit the wire frame to construct a morphological model; and finally, we render the surface and/or volume to reconstruct the object and then decompose the object in the computational domain to generate a grid. We further develop an efficient computer-aided method that unifies a morphological and kinematic modeling of 3D flyers(Liu, 2002; Liu, 2005) (H.L., manuscript in preparation).

### Kinematic modeling

A flying insect comprises a body and flapping-wing kinematics(Fig. 2). The body movements are quantified by the body angle χ (inclination of the body relative to horizontal plane), and the stroke plane angle β (plane in which the wing flaps). The wingbase-fixed coordinate system illustrated in Fig. 2A,B has its origin at the wing base, with the x-axis normal to the stroke plane, the y-axis perpendicular to the body axis, and the z-axis parallel to the stroke plane. The flapping-wing kinematics consist of three basic motions within the stroke plane: (1) flapping about the x-axis in the wingbase-fixed coordinate system, described by the positional angleφ; (2) rotation of the wing about the z-axis, described by the elevation angle θ; and (3) rotation (feathering) of the wing about the y-axis by varying the angle of attack α. Here, a general definition of the positional angle, the elevation angle and the angle of attack are given in degrees using the first three Fourier terms(Fig. 2C):
$\ {\varphi}(t)={{\sum}_{n=0}^{3}}[{\varphi}_{\mathrm{cn}}\mathrm{cos}(nKt)+{\varphi}_{\mathrm{sn}}\mathrm{sin}(nKt)],$
(A1)
$\ {\theta}(t)={{\sum}_{n=0}^{3}}[{\theta}_{\mathrm{cn}}\mathrm{cos}(nKt)+{\theta}_{\mathrm{sn}}\mathrm{sin}(nKt)],$
(A2)
$\ {\alpha}(t)={{\sum}_{n=0}^{3}}[{\alpha}_{\mathrm{cn}}\mathrm{cos}(nKt)+{\alpha}_{\mathrm{sn}}\mathrm{sin}(nKt)].$
(A3)
Note that t is dimensionless time and parameter K is the reduced frequency defined by 2πfc̄/2Uref,where f is flapping wing frequency, is the mean wing chord length(reference length) and Uref is the reference velocity at the wingtip defined by 2ΦRf, where Φ is the wing beat amplitude and R is the wing length. The Fourier coefficientsϕ cn, ϕsn, θcnsn, αcn and αsn are determined accordingly where n is integer varying from 0 to 3.

#### Regridding for a flapping wing and a moving body

The 3D movements of flapping wings and a body cause large wing deformations and 6 d.f. (degrees of freedom) displacements of the body. Modeling such movements requires an efficient and robust grid generator that fits the instantaneously deforming wing surface as well as the moved body and other boundaries. To model the 3D movements of a flapping wing(Fig. 2), we employed a previously described method (Liu and Kawachi, 1998; Liu et al.,1998; Liu, 2005)that uses the initial grid and the wing kinematics to analytically regenerate the wing-fitted grid, while minimizing additional computational requirements. The method is implemented in three steps: (1) rotating grids in the whole wing-fitted sub-domain (the grid) according to the rigid' feathering motions of the wing; (2) rotating the feathering-based grids in the whole sub-domain according to the rigid' flapping motion; and (3) rotating the feathering- and flapping-based grids in the whole sub-domain according to the elevation motion.

### Multi-blocked, overset grid method

Modeling the shape of an insect with two or four wings is a challenging problem for CFD simulations. The wings are not only undergoing large-scale movements relative to the body, but also flap rapidly, requiring us to model highly unsteady vortical flows about multiple and moving bodies. To achieve this, we develop a multi-blocked, overset grid method and incorporate it into an in-house CFD solver (Liu,2005; Aono and Liu,2006; H.L., manuscript in preparation). This grid method uses three individual structured grid systems, one for the body and the one for each wing. Each grid system is made to fit the object (body or wing), moving and deforming with the object. A trilinear interpolation technique ensures the communication of velocities and pressures among overlapping grids(Liu, 2005; Aono and Liu, 2006; H.L.,manuscript in preparation).

As shown in Fig. 1B, three grids are generated for the body and the two wings of the fruit fly. The wing grid comprises 45×45×31 cells with the outer boundary 2 mean chord lengths away from the wing surface; the grids for both wings are identical copies, using the relation of geometrical symmetry of the two wings about the body axis. The body grid is much larger because it is used as a background grid to envelope the two wing grids for the interpolation; it comprises 45×47×95 cells, and the grid is approximately 20 mean chord lengths wide (measured as the distance between outer boundary and body surface).

### Solutions to the Navier-Stokes equations

The governing equations are the three-dimensional, incompressible unsteady Navier-Stokes (NS) equations, written in a strong conservative form for momentum and mass, and non-dimensionalized in an integral form, such that:
$\ {{\int}}_{V(t)}\left(\frac{{\partial}\mathbf{q}}{{\partial}{\tau}}\right)\mathrm{d}V+\frac{{\partial}}{{\partial}t}{{\int}}_{V(t)}\mathbf{Q}\mathrm{d}V+{{\int}_{S(t)}}(\mathbf{f}-\mathbf{Qu}_{\mathrm{g}}){\cdot}\mathbf{n}\mathrm{dS}=0.$
(A4)
where the term f=(F+Fv, G+Gv, H+Hv) represents the net flux across the cell surfaces. The sub-terms are defined as:
$\mathbf{Q}=\left[\begin{array}{c}u\\v\\w\\0\end{array}\right],\mathbf{q}=\left[\begin{array}{c}u\\v\\w\\p\end{array}\right],\mathbf{F}=\left[\begin{array}{c}u^{2}+p\\uv\\uw\\{\gamma}u\end{array}\right],\mathbf{G}=\left[\begin{array}{c}vu\\v^{2}+p\\vw\\{\gamma}v\end{array}\right],\mathbf{H}=\left[\begin{array}{c}wu\\wv\\w^{2}+p\\{\gamma}w\end{array}\right],\mathbf{F}_{\mathrm{v}}=-\frac{1}{Re}\left[\begin{array}{c}2u_{\mathrm{x}}\\u_{\mathrm{y}}+u_{\mathrm{x}}\\u_{\mathrm{z}}+w_{\mathrm{x}}\\0\end{array}\right],\mathbf{G}_{\mathrm{v}}=-\frac{1}{Re}\left[\begin{array}{c}v_{\mathrm{x}}+u_{\mathrm{y}}\\2v_{\mathrm{y}}\\v_{\mathrm{z}}+w_{\mathrm{y}}\\0\end{array}\right],\mathbf{H}_{\mathrm{v}}=-\frac{1}{Re}\left[\begin{array}{c}w_{\mathrm{z}}+u_{\mathrm{x}}\\w_{\mathrm{y}}+u_{\mathrm{z}}\\2w_{\mathrm{z}}\\0\end{array}\right].$
In the preceding equations, γ is the pseudo-compressibility coefficient; p is pressure; u, v, w are the x, y, z velocity components in the Cartesian coordinate system; t denotes physical time, τ is pseudo time; Re is the Reynolds number. Note that the term q associated with pseudo time is designed for an inner-iteration at each physical time step, and will vanish when the divergence of velocity is driven to zero so as to satisfy the equation of continuity. Time-dependent solutions of the incompressible NS equations are formulated in an ALE(Arbitrary Lagrangian–Eulerian) manner with the FVM (Finite Volume Method) and are performed in a time-marching manner with a pseudo-compressibility method; we enforce conservation of mass and momentum in both time and space. More details can be found elsewhere(Liu and Kawachi, 1998).

### Boundary conditions

As shown in Fig. 1, the solutions to the NS equations with a multi-blocked, overset grid for a flapping-flying insect require appropriate boundary conditions for the overlapping zones among the different single grid block, the moving walls of the wing and the body, and the far-field outside boundary.

For individual grid blocks and for the total wing and total body we use the fortified solutions to the NS equations by adding a forcing term with communication of a vector q* to offer the boundary conditions for velocity and pressure in the overlapping zones of the two grids(H.L., manuscript in preparation). The fortified equations are solved inside the computational domain, except for holes and the single grid boundary. In the case of a fruit fly, each time step requires that we solve the fortified NS equations three times, once for each grid cell block. On the body surface,the no-slip condition is applied to calculate the velocity components. To account for dynamic effects due to the accelerations of the oscillating body(moving and/or deforming body surface), pressure divergence at the surface stencils is derived from the local momentum equation, such that
$\ (u,v,w)=(u_{\mathrm{body}},v_{\mathrm{body}},w_{\mathrm{body}}),$
(A5)
$\ {\partial}p{/}{\partial}\mathbf{n}=-\mathrm{a}_{0}{\cdot}\mathbf{n},$
(A6)
where the velocity (ubody, vbody,wbody) and the acceleration (a0) on the solid wall are evaluated and updated using the renewed grids on the body surface at each time step.

For the background grid of the insect body we need to define appropriate boundary conditions at the outside boundary(Fig. 1B). Consider that, when an insect hovers or flies forward at a velocity Vf, the boundary conditions for the velocity and the pressure may be given such as:(1) at upstream V(u, v, w)=Vf while pressure p is set to zero; (2) at downstream zero-gradient condition is taken for both velocity and pressure, i.e. ∂(u, v, w,p)/∂n=0, where n is the unit outward normal vector at the outside boundary.

### Evaluation of forces, torques and powers

Given the wing kinematics, according to the Newton's second law, the non-dimensional inertial forces acting upon the wing may be calculated as the sum of the inertial forces on all wing cells, such that:
$\ \mathbf{F}_{\mathrm{iner}}^{{\ast}}=-{{\sum}_{i}}\frac{\mathrm{d}[{\rho}_{\mathrm{w}}{/}{\rho}{\Delta}v_{\mathrm{i}^{\mathrm{i}}}\mathbf{v}^{{\ast}}(t)]}{\mathrm{d}t},$
(A7)
where wing mass density ρw is assumed to be constant throughout the wing, Δvi is the volume of the cell constructed by eight grid points on the wing, and
$$\mathrm{v}_{\mathrm{i}}^{{\ast}}$$
(t) is the computed wing velocity at each cell center at time t. The aerodynamic forces,
$$\mathbf{F}_{\mathrm{aero}}^{{\ast}}$$
of lift and drag coefficients acting on the wing surface can be calculated from the pressure and stresses along its surface based on the solutions to the NS equations. The resultant lift and thrust forces are calculated first in the local wingbase-fixed coordinate system (x, y, z), and then transformed into the earth-based coordinate system (X, Y, Z)accounting for the stroke plane angle (Fig. 2B), yielding vertical, horizontal and sideslip forces. Both the aerodynamic and inertial forces are non-dimensionalized with the reference velocity Uref, the reference length , the air density ρ and the planform area of the wing Sw, such that:
$\ \mathbf{F}_{\mathrm{iner}}^{{\ast}}=\frac{\mathbf{F}_{\mathrm{iner}}}{0.5{\rho}U_{\mathrm{ref}}^{2}S_{\mathrm{w}}},$
(A8)
$\ \mathbf{F}_{\mathrm{aero}}^{{\ast}}=\frac{\mathbf{F}_{\mathrm{aero}}}{0.5{\rho}U_{\mathrm{ref}}^{2}S_{\mathrm{w}}}.$
(A9)
The non-dimensionalized inertial and aerodynamic torques of the wing are calculated as the sum of the cross product of each force and the positional vector at each cell center of the wing about the origin of body (O″),such that:
$\ \mathbf{T}_{\mathrm{iner}}^{{\ast}}={{\sum}_{i}}(\mathbf{r}_{\mathrm{i}}^{{\ast}}{\times}\mathbf{F}_{\mathrm{iner},\mathrm{i}}^{{\ast}}),$
(A10)
$\ \mathbf{T}_{\mathrm{aero}}^{{\ast}}={{\sum}_{i}}(\mathbf{r}_{\mathrm{i}}^{{\ast}}{\times}\mathbf{F}_{\mathrm{aero},\mathrm{i}}^{{\ast}}),$
(A11)
where
$$\mathbf{r}_{\mathrm{i}}^{{\ast}}$$
denotes the non-dimensionalized positional vector of the cell center on the wing surface,and
$$\mathbf{F}_{\mathrm{aero},\mathrm{i}}^{{\ast}}$$
and
$$\mathbf{F}_{\mathrm{iner},\mathrm{i}}^{{\ast}}$$
represent the aerodynamic and inertial forces at each cell center, respectively.
The aerodynamic and inertial powers are calculated as the scalar products of the velocity and the aerodynamic and inertial forces of the wing as:
$\ P_{\mathrm{iner}}^{{\ast}}={{\sum}_{i}}(\mathbf{F}_{\mathrm{iner},\mathrm{i}}^{{\ast}}{\cdot}\mathbf{v}_{\mathrm{i}}^{{\ast}}),$
(A12)
$\ P_{\mathrm{aero}}^{{\ast}}={{\sum}_{i}}(\mathbf{F}_{\mathrm{aero},\mathrm{i}}^{{\ast}}{\cdot}\mathbf{v}_{\mathrm{i}}^{{\ast}}),$
(A13)
where
$$\mathbf{v}_{\mathrm{i}}^{{\ast}}$$
is the computed wing velocity at the cell i. Furthermore, the muscle-mass-specific aerodynamic and inertial powers can be calculated as:
$\ P_{\mathrm{iner}}=P_{\mathrm{iner}}^{{\ast}}{/}M_{\mathrm{m}},$
(A14)
$\ P_{\mathrm{aero}}=P_{\mathrm{aero}}^{{\ast}}{/}M_{\mathrm{m}}$
(A15)
where the mass of flight muscle, Mm is assumed to contribute to 30% of the total body mass (Lehman and Dickinson, 1997). The muscle-mass-specific total mechanical power Ptotal is simply the sum of the muscle-mass-specific aerodynamic and inertial powers,such that:
$\ P_{\mathrm{total}}=P_{\mathrm{iner}}+P_{\mathrm{aero}}.$
(A16)

### Verification and validation

A variety of benchmark tests for verification and validation have been undertaken to assess the reliability of the single-grid NS solver for unsteady flows about a moving and deforming body(Liu and Kawachi, 1998; Liu et al., 1998; Liu, 2002). Verification and validation for the present multi-block- and overset grid-based in-house NS solver were further conducted through an extensive study of unsteady flows past a single, rowing-feathering fin with a single BFC (body-fitted coordinate) grid (81×31×31) and with a two-block grid consisting of a single grid (case1: 51×31×25 and case2: 31×25×13)fitted to the fin and a cubic background grid (81×31×31), at a Reynolds number Re of 1.597×104 and the reduced frequency K of 3.0 (Liu and Kato,2004 see Fig. 5). The computed results with the two-block grid match very well those of the single-grid, and the computed time course of three force coefficients. Cx, Cy and Cz,show reasonable agreement with the measurements(Liu and Kato, 2004) (see Fig. 6).

Verification of the present integrated, computational system is performed with a specific focus on its self-consistency in terms of grid refinement and time step effect. A grid sensitivity analysis demonstrates that the wing grid(33×35×19) and the body grid (33×35×35) achieved reasonably accurate solutions for a hovering fruit fly. As illustrated in Fig. A1A, a maximum difference in lift and drag coefficients is obtained within 5% among this set of grids and a set of fine grids (wing grid: 45×45×31, body grid:45×43×95) and finest grids (wing grid: 45×45×31, body grid: 57×55×121). The effect of time increment on the force generation is further investigated using two time steps of 0.005 and 0.0025,and the computed results show almost no difference between the two cases;hence a physical time step of 0.005 is used throughout the simulations(Fig. A1B). An extended study of validation of the hovering fruit fly is further discussed in the Results by comparing the vertical (lift) and horizontal (drag and thrust) forces as well as the powers with the experimental data(Fry et al., 2005).

Fig. A1.

Effects of grid and time step on aerodynamic force production. (A) Grid sensitive analysis for aerodynamic force coefficients over a flapping cycle. Three grid systems are used: c, a coarse grid system (wing:33×35×19, body: 33×35×35); fn, a fine grid system(wing: 45×45×31, body: 45×47×95); and ft, a finest grid system (wing: 45×45×31, body: 57×55×121). (B)Time step sensitive analysis for aerodynamic force coefficients over a flapping cycle. Two time steps are used: t1, time step dt=0.005; t2,time step dt=0.0025.

Fig. A1.

Effects of grid and time step on aerodynamic force production. (A) Grid sensitive analysis for aerodynamic force coefficients over a flapping cycle. Three grid systems are used: c, a coarse grid system (wing:33×35×19, body: 33×35×35); fn, a fine grid system(wing: 45×45×31, body: 45×47×95); and ft, a finest grid system (wing: 45×45×31, body: 57×55×121). (B)Time step sensitive analysis for aerodynamic force coefficients over a flapping cycle. Two time steps are used: t1, time step dt=0.005; t2,time step dt=0.0025.

List of symbols and abbreviations

• a0

acceleration (or deceleration) of the insect wing and body

•
• ALE

arbitrary Lagrangian–Eulerian

•
• APT

aerodynamic pitching torque

•
• ART

aerodynamic rolling torque

•
• AYT

aerodynamic yawing torque

•
• BFC

body-fitted coordinate

•
• mean wing chord length (reference length)

•
• CFD

computational fluid dynamics

•
• Cx, Cy,Cz

dimensionless force coefficients

•
• CD

dimensionless coefficient of horizontal force

•
• CL

dimensionless coefficient of vertical force

•
• CS

dimensionless coefficients of sideslip force

•
• DPIV

digital particle image velocimetry

•
• dt

time increments

•
• f

flapping wing frequency

•
• f

net flux across cell surface

•
• FVM

finite volume method

•
• Faero

aerodynamic force

•
• Finer

inertial force

•
• F*aero

dimensionless aerodynamic force

•
• F*iner

dimensionless inertial force

•
• F*aero,i

dimensionless aerodynamic force of the cell (i)

•
• F*iner,i

dimensionless inertial force of the cell (i)

•
• i

cell index

•
• IPT

inertial pitching torque

•
• K=fc̄/2Uref

reduced frequency

•
• LEV

•
• Mm

mass of flight muscle

•
• n=(nx,ny,nz)

unit outward normal vector

•
• O

origin of earth-fixed Cartesian coordinates

•
• O′

origin of wingbase-fixed Cartesian coordinates

•
• O″

origin of insect body

•
• p

pressure

•
• P*aero

dimensionless aerodynamic power

•
• P*iner

dimensionless inertial power

•
• Paero

muscle-mass-specific aerodynamic power

•
• Piner

muscle-mass-specific inertial power

•
• Ptotal=Paero+Piner

muscle-mass-specific total mechanical power

•
• q

flux vector with respect to pseudo-compressibility

•
• q*

communication vector in overlapping zones of the two grids

•
• r*i

non-dimensionalized positional vector of the cell

•
• R

wing length

•
• Re

Reynolds number

•
• S(t)

surface of the control volume

•
• Sw

planform area of a wing

•
• t

dimensionless time

•
• T

dimensionless period of one flapping cycle

•
• TV

wing tip vortex

•
• TEV

trailing-edge vortex

•
• TPT

total pitching torque

•
• T*aero

dimensionless aerodynamic torque

•
• T*iner

dimensionless inertial torque

•
• ubody,vbody,wbody

velocity

•
• ug

velocity of a moving cell

•
• Uref=2RfΦ

reference velocity at wing tip

•
• ΔvI

volume of the cell(i)

•
• V(t)

volume of control volume

•
• Vf

velocity of the insect body

•
• v*I

dimensionless velocity of the cell (i)

•
• x, y, z

wingbase-fixed Cartesian coordinates

•
• X, Y, Z

earth-fixed Cartesian coordinates

•
• α

feathering angle (or angle of attack of the wing)

•
• β

stroke plane angle

•
• γ

pseudo-compressibility coefficient

•
• χ

body angle

•
• Φ

wing beat amplitude

•
• θ

elevation angle

•
• φ

positional angle

•
• φ cn, φ sn, θcn, θsn, αcn, αsn

coefficients of the kinematic data

•
• ρ

air density

•
• ρw

mass density of the wing

•
• τ

pseudo-time

•
• ν

kinematic viscosity of air

We thank the anonymous referees for their valuable comments and suggestions on the manuscript. This work was partially supported by a PRESTO (Precursory Research for Embryonic Science and Technology) program of the Japan Science and Technology Agency (JST), and Grant-in-Aid for Scientific Research Nos. 18656056 and 18100002, Japan Society for the Promotion of Science (JSPS). The simulations were performed in a supercomputer [Super Combined Cluster (RSCC)],RIKEN (The Institute of Physical and Chemical Research), Japan.

Aono, H. and Liu, H. (
2006
). Vortical structure and aerodynamics of hawkmoth hovering.
J. Biomech. Sci. Eng.
1
,
234
-245.
Birch, J. M. and Dickinson, M. H. (
2001
). Spanwise flow and the attachment of the leading-edge vortex on insect wings.
Nature
412
,
729
-733.
Birch, J. M. and Dickinson, M. H. (
2003
). The influence of wing–wake interactions on the production of aerodynamic forces in flapping flight.
J. Exp. Biol.
206
,
2257
-2272.
Birch, J. M., Dickinson, W. B. and Dickinson, M. H.(
2004
). Force production and flow structure of the leading edge vortex on flapping wings at high and low Reynolds numbers.
J. Exp. Biol.
207
,
1063
-1072.
Bomphrey, R. J., Lawson, N. J., Taylor, G. K. and Thomas, A. L. R. (
2005
). The aerodynamics of Manduca sexta:digital particle image velocimetry analysis of the leading-edge vortex.
J. Exp. Biol.
208
,
1079
-1094.
Bomphrey, R. J., Lawson, N. J., Taylor, G. K. and Thomas, A. L. R. (
2006
). Application of digital particle image velocimetry to insect aerodynamics: measurement of the leading-edge vortex and near wake of a hawkmoth.
Exp. Fluids
40
,
546
-554.
Dickinson, M. H. and Götz, K. G. (
1996
). The wake dynamics and flight forces of the fruit fly, Drosophila melanogaster.
J. Exp. Biol.
199
,
2085
-2104.
Dickinson, M. H., Lehmann, F.-O. and Sane, S.(
1999
). Wing rotation and the aerodynamic basis of insect flight.
Science
284
,
1954
-1960.
Ellington, C. P. (
1984a
). The aerodynamics of insect flight. I. The quasi-steady analysis.
Philos. Trans. R. Soc. Lond. B Biol. Sci.
305
,
1
-15.
Ellington, C. P. (
1984b
). The aerodynamics of insect flight. III. The kinematics.
Philos. Trans. R. Soc. Lond. B Biol. Sci.
305
,
41
-78.
Ellington, C. P. (
1984c
). The aerodynamics of insect flight. IV. Aerodynamic mechanisms.
Philos. Trans. R. Soc. Lond. B Biol. Sci.
305
,
79
-113.
Ellington, C. P. (
1984d
). The aerodynamics of insect flight. V. A vortex theory.
Philos. Trans. R. Soc. Lond. B Biol. Sci.
305
,
115
-144.
Ellington, C. P., van den Berg, C., Willmott, A. P. and Thomas,A. L. R. (
1996
). Leading-edge vortices in insect flight.
Nature
384
,
626
-630.
Fry, S. N., Sayaman, R. and Dickinson, M. H.(
2003
). The aerodynamics of free-flight maneuvers in Drosophila.
Science
300
,
495
-498.
Fry, S. N., Sayaman, R. and Dickinson, M. H.(
2005
). The aerodynamics of hovering flight in Drosophila.
J. Exp. Biol.
208
,
2303
-2318.
Grodnitsky, D. L. and Morozov, P. P. (
1993
). Vortex formation during tethered flight of functionally and morphologically two-winged insects, including evolutionary considerations on insect flight.
J. Exp. Biol.
182
,
11
-40.
Lehmann, F.-O. (
2002
). The constraints of body size on aerodynamics and energetic flying fruit flies: an integrative view.
Zoology
105
,
287
-295.
Lehmann, F.-O. (
2004a
). The mechanisms of lift enhancement in insect flight.
Naturewissenschaften
91
,
10
-122.
Lehmann, F.-O. (
2004b
). Aerial locomotion in flies and robots: kinematic control and aerodynamics of oscillating wings.
Arthropod Struct. Dev.
33
,
331
-345.
Lehmann, F.-O. and Dickinson, M. H. (
1997
). The changes in power requirements and muscle efficiency during elevated force production in the fruit fly Drosophila.
J. Exp. Biol.
200
,
1133
-1143.
Lehmann, F.-O., Sane, S. P. and Dickinson, M. H.(
2005
). The aerodynamic effects of wing–wing interaction in flapping insect wings.
J. Exp. Biol.
208
,
3075
-3092.
Liu, H. (
2002
). Computational biological fluid dynamics: digitizing and visualizing animal swimming and flying.
Integr. Comp. Biol.
42
,
1050
-1059.
Liu, H. (
2005
). Simulation-based biological fluid dynamics in animal locomotion.
ASME Appl. Mech. Rev.
58
,
269
-282.
Liu, H. and Kato, N. (
2004
). A numerical study of unsteady hydrodynamics of a mechanical pectoral fin.
J. Bionics Eng.
1
,
108
-120.
Liu, H. and Kawachi, K. (
1998
). A numerical study of insect flight.
J. Comput. Phys.
146
,
124
-156.
Liu, H., Ellington, C. P., Kawachi, K., van den Berg, C. and Willmott, A. P. (
1998
). A computational fluid dynamic study of hawkmoth hovering.
J. Exp. Biol.
201
,
461
-477.
Miller, L. A. and Peskin, C. S. (
2004
). When vortices stick: an aerodynamic transition in tiny insect flight.
J. Exp. Biol.
207
,
3073
-3088.
Miller, L. A. and Peskin, C. S. (
2005
). A computational fluid dynamics of clap and fling' in the smallest insects.
J. Exp. Biol.
208
,
195
-212.
Poelma, C., Dickson, W. B. and Dickinson, M. H.(
2006
). Time-resolved reconstruction of the full velocity field around a dynamically scaled flapping wing.
Exp. Fluids
41
,
213
-225.
Ramamurti, R. and Sandberg, W. C. (
2002
). A three-dimensional computational study of the aerodynamic mechanisms of insect flight.
J. Exp. Biol.
205
,
1507
-1518.
Sane, S. (
2003
). The aerodynamics of insect flight.
J. Exp. Biol.
206
,
4191
-4208.
Sun, M. and Tang, J. (
2002a
). Lift and power requirements of hovering flight in Drosophila virlis.
J. Exp. Biol.
205
,
2413
-2427.
Sun, M. and Tang, J. (
2002b
). Unsteady aerodynamic force generation by a model fruit fly wing in flapping motion.
J. Exp. Biol.
205
,
55
-70.
van den Berg, C. and Ellington, C. P. (
1997a
). The three-dimensional leading-edge vortex of a hovering' model hawkmoth.
Philos. Trans. R. Soc. Lond. B Biol. Sci.
352
,
329
-340.
van den Berg, C. and Ellington, C. P. (
1997b
). The vortex wake of a `hovering' model hawkmoth.
Philos. Trans. R. Soc. Lond. B Biol. Sci.
352
,
317
-328.
Wang, H., Inada, Y. and Liu, H. (
2005
). A numerical analysis of inertial torques in the steering maneuvers of hovering Drosophila.
JSME Int. J. Ser. C
48
,
499
-512.
Wang, Z. J. (
2000
). Two dimensional mechanism for insect hovering.
Phys. Rev. Lett.
85
,
2216
-2219.
Wang, Z. J. (
2005
). Dissecting insect flight.
Annu. Rev. Fluid Mech.
37
,
183
-210.
Wang, Z. J., Birch, J. M. and Dickinson, M. H.(
2004
). Unsteady forces and flows in low Reynolds number hovering flight: two-dimensional computations vs robotic wing experiments.
J. Exp. Biol.
207
,
449
-460.
Warrick, D. R., Tobalske, B. W. and Powers, D. R.(
2005
). Aerodynamics of the hovering hummingbird.
Nature
435
,
1094
-1097.
Weis-Fogh, T. (
1973
). Quick estimates of flight fitness in hovering animals, including novel mechanisms for lift production.
J. Exp. Biol.
59
,
169
-230.
Willmott, A. P., Ellington, C. P. and Thomas, A. L. R.(
1997
). Flow visualization and unsteady aerodynamics in the flight of the hawkmoth, Manduca sexta.
Philos. Trans. R. Soc. Lond. B Biol. Sci.
352
,
303
-316.