Near- and far-field aerodynamics in insect hovering flight: an integrated computational study

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.

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.

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 c̄/32, where c̄ 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).

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̄U_ref/ν,where ν is the kinematic viscosity of air (1.5×10^–5m² s^–1). According to the measured data(c̄=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).

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.

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 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.8U_ref 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).

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 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.

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).

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.

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).

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).

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).

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 P_iner 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 P_aero 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 P_total is the power required to move the wing. We simply calculate P_total using Eqn A16 as the sum of the muscle-mass-specific aerodynamic and inertial powers. It can be seen that the curve of P_total 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).

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).

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.

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.

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 Gö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 Gö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.

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. 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.

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.

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).

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.

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).

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 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 V_f, the boundary conditions for the velocity and the pressure may be given such as:(1) at upstream V(u, v, w)=V_f 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.

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×10⁴ 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. C_x, C_y and C_z,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).

List of symbols and abbreviations

 
• a₀

  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

 
• c̄

  mean wing chord length (reference length)

 
• CFD

  computational fluid dynamics

 
• C_x, C_y,C_z

  dimensionless force coefficients

 
• C_D

  dimensionless coefficient of horizontal force

 
• C_L

  dimensionless coefficient of vertical force

 
• C_S

  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

 
• F_aero

  aerodynamic force

 
• F_iner

  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=2πfc̄/2U_ref

  reduced frequency

 
• LEV

  leading-edge vortex

 
• M_m

  mass of flight muscle

 
• n=(n_x,n_y,n_z)

  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

 
• P_aero

  muscle-mass-specific aerodynamic power

 
• P_iner

  muscle-mass-specific inertial power

 
• P_total=P_aero+P_iner

  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

 
• S_w

  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

 
• u_body,v_body,w_body

  velocity

 
• u_g

  velocity of a moving cell

 
• U_ref=2RfΦ

  reference velocity at wing tip

 
• Δv_I

  volume of the cell(i)

 
• V(t)

  volume of control volume

 
• V_f

  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.