A three-dimensional computational fluid dynamics simulation is performed for a ruby-throated hummingbird (Archilochus colubris) in hovering flight. Realistic wing kinematics are adopted in the numerical model by reconstructing the wing motion from high-speed imaging data of the bird. Lift history and the three-dimensional flow pattern around the wing in full stroke cycles are captured in the simulation. Significant asymmetry is observed for lift production within a stroke cycle. In particular, the downstroke generates about 2.5 times as much vertical force as the upstroke, a result that confirms the estimate based on the measurement of the circulation in a previous experimental study. Associated with lift production is the similar power imbalance between the two half strokes. Further analysis shows that in addition to the angle of attack, wing velocity and surface area, drag-based force and wing–wake interaction also contribute significantly to the lift asymmetry. Though the wing–wake interaction could be beneficial for lift enhancement, the isolated stroke simulation shows that this benefit is buried by other opposing effects, e.g. presence of downwash. The leading-edge vortex is stable during the downstroke but may shed during the upstroke. Finally, the full-body simulation result shows that the effects of wing–wing interaction and wing–body interaction are small.
Unlike birds of larger sizes, hummingbirds can perform sustained hovering in addition to regular cruise flight and manoeuvres. Many studies have been done to characterize the kinematics, physiology and aerodynamics of the hummingbird wing [1–4], and they were summarized in the work of Tobalske et al. . In general, hummingbirds use similar aerodynamics to those of insects, e.g. presence of a leading-edge vortex (LEV) over the wing surface [6,7], for lift production. However, differences between hummingbird and insect aerodynamics are conceivable as the anatomy and physiology of the hummingbird wing are distinct from those of the insect wing. For example, recent evidence shows that hummingbirds can achieve the inversion of the angle of attack through active wing rotation at the wrist . This actuation mechanism is different from that of insects whose wing inversion can be realized through combined muscle activation at the wing root and the passive deformation of the wing surface . The implication of this difference on the lift and power efficiency of hummingbirds is still unclear. In order to better understand aerodynamics of hummingbirds, their lift and flow characteristics are needed. Unfortunately, such data are so far very limited. To set the context for this work, we briefly summarize several recent studies on the force production and flow behaviour in hummingbird flight.
Altshuler et al.  used a dissected hummingbird wing and tested lift production of the wing revolving in one direction. By comparing with wing models of increasing realism, i.e. those with sharpened leading edges and with substantial camber, they found that the real hummingbird wing generates more lift, and their result suggests that some geometric details such as the presence of camber tend to increase lift. Using particle image velocimetry (PIV), Warrick et al. [4,10] studied the flight of rufous hummingbirds. They were able to measure the flow in the two-dimensional planes that are perpendicular to the wing axis during the entire stroke cycle. Based on the PIV data, they visualized the LEV and calculated the circulation at different spanwise locations. Interestingly, their result shows that the average bound circulation during the downstroke is 2.1 times of that during the upstroke . Assuming that the conventional aerofoil theory holds, that is, lift is linearly proportional to the bound circulation, the authors suggest that the lift production possesses the same amount of asymmetry. They further proposed the possible mechanisms that may have contributed to such lift characteristics. For example, the wing velocity and the angle of attack during the downstroke are greater than those during the upstroke. Other variables they suggested include longer wing span and formation of a positive camber during the downstroke. In another PIV study, Altshuler et al.  measured the wake flow of the wings and tail of hovering Anna's hummingbirds, and they proposed a vortex-ring model for the wake generated by the wings. Later, Wolf et al.  conducted further PIV study of the same hummingbird species, and from the strength of the shed vortices, they also concluded that lift production is highly asymmetric between the two half strokes.
Despite these previous efforts, there exists no direct study on the detailed force characteristics and the three-dimensional flow pattern of the hummingbird wing in hovering flight. As a useful tool, computational fluid dynamics (CFD) has been applied in many previous works to study aerodynamics of flapping wings, including both rigid and flexible wing models [13–16]. Here, we are motivated to perform a CFD study to quantify the force histories in a stroke cycle and to investigate any particular force production mechanisms used by the hummingbird. The main questions we would like to answer through this work include: (i) what are the characteristics of the force history, and what are the underlying mechanisms for the possible downstroke–upstroke asymmetry? (ii) What is the three-dimensional wake pattern like, and how may it be associated with the force characteristics? (iii) How much aerodynamic power does the hummingbird have to spend on hovering and what is the efficiency? (iv) Does the hummingbird use some of the mechanisms that insects use for lift enhancement, such as wake capture and wing–wing interaction?
2. Material and method
2.1. Experiment and reconstruction of the wing kinematics
The hummingbird, a female ruby-throated (Archilochus colubris) with a body mass of 3.41 g, is used as the subject in this study. High-speed filming experiment was conducted to record the wing motion of the bird. In the experiment, the bird was trained to fly in a 0.4 × 0.4 × 0.5 m3 netted chamber and was recorded 1000 frames per second with a 1/5000th shutter by four high-speed cameras: two Phantom v. 7.1 (Vision Research Inc., Wayne, NJ, USA), one Photron SA-3 and one Photron 1024 PCI (Photron USA Inc., San Diego, CA, USA). Each flapping cycle contains about 25 frames. The bird was labelled prior to the experiment using 1-mm diameter dots of non-toxic white paint, as shown in figure 1a. The experimental set-up is described in detail in Hedrick et al. . The nine markers numbered in figure 1a and located on the outline of the left wing are used in this study. These markers include five points on the leading edge, one at the wingtip and three on the trailing edge. To avoid blind spots, the cameras were positioned with one directly behind the bird in the same horizontal plane, two with an elevated oblique and slightly rear view and one with a ventral view of the bird (see the electronic supplementary material, figure S1). After the videos are taken, a custom MATLAB program  was used to automatically track the markers frame by frame and to extract their three-dimensional coordinates. A principal components analysis has been done to verify that these nine points are sufficient to characterize the wing motion.
To reconstruct the wing geometry and motion, spline interpolation is used to connect the outline of the wing at each instantaneous time frame. Then, both the leading edge and the trailing edge are evenly discretized by 41 nodes each. The wing chord is approximated with straight segments which have rounded ends and an effective thickness 7% of the average chord length.
A triangular mesh is then generated to discretize the wing surface, which is assumed to be smooth. Corrugations caused by the feathers are ignored as their effect on the laminar boundary layer is expected to be small at the current Reynolds number. Discussions on the effect of feathers at higher Reynolds numbers can be found in a recent experimental study . A single wing consists of 1129 elements and 615 Lagrangian nodes. To increase the time resolution for the small-step solution of the simulation, the trajectory of each mesh node is also refined by spline interpolation in time. Eight cycles of wing kinematics are reconstructed from the imaging data and are used for the simulation (see the electronic supplementary material, movie). Note that dynamic deformations of the wing such as spanwise bending and twisting have been included in the reconstructed kinematics (see the electronic supplementary material, figure S2) and thus their aerodynamic consequences will be incorporated in the simulation results.
As seen in figure 1b, the entire wing surface exhibits a twist along the wing axis, and the twist angle changes dynamically in a stroke cycle due to the pitching motion of the wing. To define the wing posture and the time-varying angle of attack, we select three points on the wing: the wing tip, the leading-edge point and the trailing-edge point of the mid-chord. These three points form a triangle approximating posturing of the distal half of the wing surface, as indicated in figure 1c. The chord angle, αc, is defined as the instantaneous acute angle between the plane spanned by this triangle and the horizontal plane. This angle will be used to measure orientation and pitch rotation of the distal wing surface. The angle of attack, α, is defined as the instantaneous angle between the tip velocity vector and the triangle.
2.2. Simulation set-up and model validation
The air is assumed to be governed by the viscous incompressible Navier–Stokes equation. The equation is solved by a second-order accurate immersed-boundary method  that is able to handle large displacement of the moving boundaries (see numerical method in the electronic supplementary material). A fixed, non-uniform, single-block Cartesian grid is employed to discretize the domain (figure 2a). The rectangular domain is 20 × 20 × 18 cm3. For the single-wing simulation, 330 × 250 × 210 (17 million) points are used for the baseline simulation. A coarser mesh with 232 × 180 × 140 (6 million) points and a finer mesh with 420 × 310 × 240 (31 million) points are also used in the single-wing case to verify grid convergence. All three meshes have maximum resolution around the wing, which is 0.05, 0.033 and 0.025 cm in all three directions, for the coarser, baseline and finer mesh, respectively. The two extra simulations are run for two cycles, and they produce a maximum 3% difference from the baseline mesh in the mean and RMS values of the vertical force. The full-body simulation employs 336 × 408 × 216 (30 million) points, and the resolution around the body and wings is the same as in the baseline case for the single wing.
The numerical method has been previously validated for flapping-wing simulations against both experimental and simulation data in Dai et al. , where a fruit fly model and an impulsively started plate were studied. To further validate the model in this work, we compare the flow field with that obtained from the PIV experiment by Warrick et al. . Note that the rufous hummingbird (Selasphorus rufus) was used in the experiment, whereas the ruby-throated hummingbird (A. colubris) is used in this study. However, these two species are very similar to each other in terms of the morphological data and wing kinematics. Table 1 lists some of the key parameters of the current hummingbird along with those from Warrick et al. , including the body mass, M, the flapping frequency, f, the wing length, R, the wing span, b, average chord length, c, the wing area, S and the wingbeat amplitude, Φ. It can be seen that all the parameters in this study fall well within the ranges in the experiment. We also converted the angle of attack and the chord angle of the present hummingbird using the definitions in the experimental study, and the result of comparison is generally consistent (in Tobalske et al. , the chord angle is 14° ± 7° for downstroke and 31° ± 4° for upstroke; in Warrick et al. , the angle of attack is 36° ± 12° and 26° ± 13°. In this study, the chord angle is 16° for downstroke and 48° for upstroke; and the angle of attack is 33° and 24°. All angles are measured according to their definitions). The Reynolds number of the flow is set to be . This non-dimensional number represents the ratio between the fluid inertia and the viscous effect.
Figure 3b shows a typical spanwise slice of the instantaneous flow during mid-downstroke at 70% wingspan from the wing root. Note that the experimental data are shown for the slice at approximately 80% wingspan, or 4 cm from the wing root. A discussion on the choice of the spanwise location is deferred to the end of this section. It can be seen that in both cases, a strong shear layer exists on the dorsal surface of the wing and is generally attached to the wing surface. In the experiment, the shear layer on the ventral side of the wing is not visible due to the shadow effect. Both figures show that a large clockwise vortex is located in the wake of the wing and is about one chord length away from the trailing edge, though the strength of vortex is weaker in the simulation. Overall downwash is created in both cases, which corresponds to lift production. There are also other visible differences between the two plots. In particular, the vortices in the experiment appear to be multiple blobs above the wing surface, whereas in the simulation a continuous vortex sheet is formed and is slightly separated from the wing near the leading edge. Comparison of other time frames also displays similarities in general flow patterns but considerable differences in flow details (see the electronic supplementary material, figures S3–S6). We point out that variations in the wing kinematics of bird individuals may have led to discrepancies in the flow field observed here. In addition, some of these differences are likely caused by low resolution in the experiment where around 17 points per centimetre were used for the velocity field. In the simulation, 30 points per centimetre in the baseline grid and 40 points per centimetre in the finest grid are employed around the wing. Furthermore, the two grids displayed a consistent form of shear layers.
We further compare the bound circulation around the wing chord with the data from the experiment. Figure 4 shows the phase-averaged circulation, Γ, defined as , along a circular path that encloses the wing chord. The diameter of the circle is 10% greater than the chord length. Increasing this diameter by 20% only changes the maximum circulation by 5%. In Warrick et al. , the phase-averaged circulation is shown at 80% wingspan for the entire stroke cycle. However, their results also show that the spanwise location of the maximum circulation varies largely among the bird individuals, although in general the maximum happens between 40 and 80% of wingspan. In this study, we found that the maximum bound circulation takes place between 50 and 70% of wingspan. Therefore, we plot Γ for 50, 70 and 80% wingspan locations and compare them with the experiment data. For the same reason, in the validation of the flow field, we chose to use the slice at 70% of wingspan.
Figure 4 shows that the present circulation at 50% wingspan matches the best with the experimental data. At both 70 and 80% wingspan, the circulation has a significant drop after the mid-downstroke. In the experimental result, the ratio of the downstroke and upstroke circulations is 2.1 ± 0.1 in magnitude. In the simulation, this ratio is 2.2, 2.3 and 2.0 for 50, 70 and 80% wingspan, respectively.
3. Results and discussion
We first report the forces, power and efficiency of the single-wing simulation and then discuss the characteristics of the forces and flow field. In the end, we will also discuss the full-body simulation.
3.1. Force, power and efficiency
The global coordinate system is shown in figure 2a, where X, Y and Z denote the forward, spanwise and vertical direction, respectively. The resultant force components, FX, FY and FZ, are normalized by the fluid density, ρ, the average wing area, S and the average tip velocity, , according to 3.1 where CX, CY and CZ are the force coefficients and is the coefficient of the second moment of area of the wing surface about the axis passing through the wing base point and parallel to the wing. In this study, S = 5.68 cm2 and m s−1 are averaged from the reconstructed wing motion. The air density is chosen to be 1.23 kg m−3. From these data, the reduced frequency of the wing as defined by is 0.16.
Figure 5 shows the time courses of the force coefficients and power coefficient. Note that the cycle-to-cycle variations seen in this figure are due to the non-periodic features in the wing kinematics. The aerodynamic power here is calculated by directly integrating the dot product between the wing velocity and the aerodynamic loading over the entire wing surface. The power coefficient is defined by normalizing the power by , where is the dimensionless third moment of the area of the wing. From the result, the average vertical force coefficient is which corresponds to 3.12 g of total weight support provided by two wings. The total lift is about 91% of the weight of the bird. The remaining lift could be provided by the wing camber [3,4], which is not incorporated into the current model.
The most striking feature of the vertical force is that the downstroke produces clearly much higher lift than the upstroke. The data show that CZ averaged during the downstroke is 2.5 times of that during the upstroke, which is generally consistent with the lift estimated based on the circulation in the experiments [4,10]. Note that the ratio of the bound circulation between the downstroke and upstroke is 2.1 ± 0.1 in Warrick et al. . Another observation in figure 5 is that the forces and power contain a significant dip during the upstroke. This dip corresponds to the LEV shedding from the wing, which will be discussed in §3.4.
The averaged forward force coefficient is , which is much smaller than . The average spanwise force coefficient is . These forces can be cancelled out for the real bird when taking into consideration of two-wing symmetry (for the Y-direction), tail motion and possibly the bird–feeder interaction in the imaging experiment (the latter two for the X-direction).
The power coefficient in figure 5 also exhibits similar asymmetry as the vertical force coefficient. Further calculation shows that the downstroke requires 2.8 times as much power as the upstroke. The averaged power coefficient throughout the cycles is . Defined as the ratio between the lift coefficient and the power coefficient, the aerodynamic efficiency of the wing is thus . Using the dimensional values of ρ, and S, and the body mass, we obtain the mass-specific power of the bird, which is 55 W kg−1. Altshuler et al.  estimated the power consumption of the hummingbirds using the empirically derived drag coefficient measured from a revolving hummingbird wing. For the hummingbirds flying at elevation below 1000 m (body mass ranging from 2.5 to 9 g), the mass-specific power for hovering was estimated to be between 23 and 33 W kg−1 in their work, which is about half of the current result.
Chai & Dudley  reported the oxygen consumption and therefore metabolic power input of ruby-throated hummingbirds to be around 260 W kg−1. Thus, our aerodynamic power output implies a muscle efficiency of 21%. Vertebrate muscle efficiency can reach slightly less than 30%, but hummingbirds are expected to be slightly less efficient because of adaptations for maintaining continuous high mass-specific power output and due to the unmeasured cost of accelerating the wing mass during each half stroke.
The overall muscle efficiency of 21% found here is substantially greater than that reported in earlier studies  that use simpler models to predict aerodynamic power requirements and report efficiencies of around 10%. However, other recent Navier–Stokes simulations of hovering animal flight have also reported higher power requirements than predicted  and that revolving wing experiments do not necessarily reproduce the same flow conditions and thus force coefficients as flapping wings .
3.2. Circulation and wing rotation
As shown in figure 4, the bound circulation around the wing chord is consistent with the measurement of Warrick et al. . Furthermore, the circulation is sustained through the wing reversal. For example, during the downstroke, circulation around the translating wing is developed, and towards the end of downstroke and throughout supination, the circulation does not vanish but remains the same sign, e.g. positive or counterclockwise from the right side view. Similarly, the circulation developed during the upstroke translation remains negative throughout pronation, as shown in figure 4. The lingering circulation is caused by the pitching rotation of the wing around its own axis . Unlike a spinning cylinder in a uniform flow, this circulation cannot always be used for lift production (e.g. when the translational speed is zero or the wing surface is vertical and thus has zero projected area on the horizontal). Therefore, the vertical force as shown in figure 5a is still nearly zero at wing reversals.
To better see the phase relationship between the lift production and the wing motion, we plot in figure 6 the vertical force coefficient, the translational velocity of the wing as represented by the tip velocity, Vtip, the chord angle αc and also the pitching velocity represented by . Fewer cycles are plotted henceforth to show details within a cycle, although statistics are taken from all cycles available. From this figure, we may see additional pitching effect other than pronation and supination: during mid-downstroke, there is a positive peak in and this peak also roughly corresponds to the maximum translational speed of the wing. Such backward pitching rotation would increase the circulation and, along with the wing translation, help to enhance lift production during the downstroke. On the other hand, during mid-upstroke, the magnitude of the negative peak in is much lower. This difference could have increased the force asymmetry between the downstroke and upstroke, as will be discussed in detail next.
3.3. Asymmetric lift production
3.3.1. Force asymmetry
Figure 5a shows that lift production is highly asymmetric, with the downstroke generating much greater weight support than the upstroke. The average vertical force provided by the downstroke is 0.022 N and by the upstroke is 0.0090 N. Thus, the ratio of asymmetry is 2.5. Table 2 further lists the lift coefficient, the power coefficient and the lift-to-power ratio separately for the downstroke and upstroke. It can be seen that the downstroke produces more lift, but it is also more power-consuming. By rescaling the lift and power using the respective wing velocity and surface area of each half stroke to obtain and , we see that the lift-to-power coefficient is similar for the downstroke and upstroke. Thus, despite that their aerodynamics are quite different, the two half strokes still have similar efficiency.
In Warrick et al. , the force asymmetry between the upstroke and downstroke was attributed to several mechanisms, including the wing velocity, angle of attack, wing area and camber. Except that the camber effect is not included in this study, all the other mechanisms have been observed in the simulation, as will be discussed next. In addition, we found that other mechanisms also have contributed to the asymmetry, which include the drag-based vertical force, wing–wake interaction and pitching rotation. The effect of pitching rotation has been discussed in §3.2. So we will focus on the other effects.
First, table 2 provides the comparison of a few key kinematic parameters between the downstroke and upstroke, including the average tip speed, angle of attack and wing area. It can be seen that the ratio of the average wing area between the downstroke and upstroke is only 1.09, and the ratio of the average tip speed is only 1.13. The ratio between the velocity squares is 1.21 only. That is, the combination of the wing area and velocity is much less than the ratio of 2.5 in the force asymmetry. Therefore, some other mechanisms must be significant in leading to the large imbalance of two half strokes.
3.3.2. Drag-based vertical force
First, we consider the effect of deviation, i.e. the non-reciprocal path of the wing in a stroke cycle. Observing the wing motion from the side view, we note that the wing tip traces a roughly elliptical path whose long axis has a small angle with respect to the horizontal plane. This deviation from the mean stroke plane is shown in figure 7a by plotting the cycle-averaged trajectory of the right wing tip in the XZ-plane. In the figure, the mean stroke plane is tilted forward by approximately 12° with respect to the horizontal. This observation motivates us to decompose the forces generated by the wing into the aerodynamic lift, i.e. the force perpendicular to the wing translation, and the aerodynamic drag, i.e. the force opposite to the wing translation, as illustrated in figure 7a.
To do this analysis, we define the instantaneous stroke plane as the plane spanned by the instantaneous tip velocity vector and the wing axis. The instantaneous stroke plane angle, β, is the angle between this plane and the horizontal plane. Both the instantaneous and cycle-averaged values of β are plotted in figure 7b, along with the instantaneous and cycle-averaged angle of attack α. Note that these two angles are defined in the three-dimensional space and are shown in the two-dimensional plot in figure 7a for illustration purposes only. It can be seen that after the pronation, β is around −50° and then drops in magnitude during more than half of the downstroke. During early downstroke, the angle of attack is large and drops from 80° to 39°. The two angle histories indicate that during early downstroke, the wing is pressing downward while sweeping forward. Towards the end of downstroke, β becomes positive, but its magnitude is less than 25° before supination. In comparison, during the upstroke β is around 10° and only varies slightly.
We define the resultant force normal to the instantaneous stroke plane as lift, FL, and the force opposite to the direction of the instantaneous tip velocity as drag, FD. Figure 8 shows the normalized lift and drag by , CL and CD, and also their projections in the vertical direction, CL,Z and CD,Z, for two cycles. In figure 8a, CL and CD correlate with each other and have similar magnitude. The average data are listed in table 3 separately for the downstroke and upstroke. Figure 8b shows that during the downstroke, drag has a significant positive contribution to the vertical force during the first half of the downstroke and has only a small negative contribution during the second half of the downstroke. During the upstroke, drag has a mostly negative contribution, and the magnitude is small. On average, the drag-based vertical force, CD,Z, is 0.63 or 24% of the total vertical force CZ during the downstroke, and it is −0.12 or 12% of CZ during the upstroke. As CD,Z of the downstroke is 61% of CZ of the upstroke and CD,Z of the upstroke is small, we can conclude that drag contributes 0.61 out of the asymmetry ratio 2.5 in the vertical force.
Figure 8a,b also show that after excluding the drag-based vertical force, the lift coefficient, CL, is still asymmetric between the downstroke and upstroke, and so is its vertical component, CL,Z. On average, the downstroke-to-upstroke ratio in CL is 1.80.
3.3.3. Wing speed and angle of attack
As pointed out by Warrick et al. , the differences in the translational speed and angle of attack between the downstroke and upstroke may have been a major effect for the lift asymmetry. To test this hypothesis, we designed a revolving-wing model for the current hummingbird. In this model, a rigid wing with a flat surface is created by projecting the actual wing during a mid-downstroke onto a plane (so the spanwise twist is eliminated), and the modified wing accelerates from the stationary position to a maximum velocity and then continues to revolve at that velocity (see the electronic supplementary material, figure S7). The stroke plane is horizontal, and the angle of attack is kept constant throughout the entire process. Two cases are simulated in this test. In the first case, the wing tip follows the translation history of the wing tip in an actual downstroke chosen from one typical cycle, from 0 to the maximum velocity 15 m s−1 within the time period 0.2T, and the angle of attack is α = 41°. In the second case, the wing tip follows the translation history of an actual upstroke of the same cycle, from 0 to the maximum velocity 12 m s−1 within the time period 0.15T, and the angle of attack is α = 28°. The air properties (density and viscosity) remain the same in this set-up.
The results show that the ratio of the lift during steady translation is 1.57 between the downstroke revolving wing and the upstroke revolving wing. Thus, the combined effect of translation and angle of attack is confirmed. However, it should be noted that comparing the revolving wing and flapping wing, the transient histories of lift display considerable differences, as seen in figure 9. This result suggests that the rotational motion of the flapping wing during the acceleration phase is still important.
3.3.4. Wing–wake interaction
Wing–wake interaction is a unique feature of flapping wings. In the previous study of the aerodynamics of the fruit fly, Dickinson et al.  suggest that the wing–wake interaction enhances lift production and is able to generate a peak force at the beginning of a half stroke if the angle of attack is reversed timely (which is the case for advanced pitching and symmetric pitching). It would be interesting to see to what extent a similar effect exists in hummingbird flight, and also whether this effect influences the downstroke and upstroke differently.
First, the lift graph in figure 5 shows that there is no clear peak in CZ in the beginning of either downstroke or upstroke. To investigate the presence of the wing–wake interaction, in figure 10 we visualize the flow in a XZ-plane shortly after the wing reversal by plotting the velocity vectors tangent to the plane. In figure 10a where a typical downstroke is shown, the wing moves somewhat downward and translates at a lower elevation, and this allows the wing to capture the opposite flow produced by the preceding upstroke. Note that the opposite flow also travels downward due to the overall downwash. On the other hand, in figure 10b where a typical upstroke is shown, the wing moves somewhat upward and translates at a higher elevation, and thus it misses the opposite flow produced by the preceding downstroke. Therefore, qualitatively speaking, the downstroke benefits more from the wing–wake interaction than the upstroke, although the interaction does not generate a separate lift peak because of its timing.
To further the investigation, we simulate each half stroke in separate runs with otherwise identical wing kinematics. The start and end of the simulation are based on the observation of the wing positions at pronation and supination. Thus, the effect of wing–wake interaction is excluded in such isolated wing strokes. One issue to bear in mind is that in the isolated strokes, the wing does not encounter a mean downwash as it does in the continuous strokes. The downwash reduces the effective angle of attack and thus weakens lift production.
Figure 11 shows the lift coefficient, CZ, of the isolated strokes along with the data for the continuous strokes. In the first downstroke, the two simulations produce identical results and thus are not shown. For the other strokes, notable differences can be seen between the two simulations. For downstrokes, lift produced by the isolated strokes is close to that produced by the corresponding continuous strokes, whereas for the upstrokes, the isolated strokes produce greater lift than the continuous strokes. On average, the ratio of lift between the continuous and the isolated strokes is 93.2% for downstroke and 83.1% for upstroke. This result suggests that for the present hummingbird, the lift-enhancing effect of the wing–wake interaction does not exceed the mitigating effects of other possible mechanisms present, e.g. the downwash. On the other hand, the wing–wake interaction does affect the lift asymmetry, as the downstroke-to-upstroke ratio in the vertical force is reduced to 2.2 for the isolated strokes.
Finally, it should be noted that the upstroke–downstroke force asymmetry was also observed in the hovering flight of some insects such as the hawkmoth  and fruit fly , though for the fruit fly the upstroke produces greater vertical force. It may be possible that some of the effects discussed in this study have led to the observed asymmetry. For example, from the tip trajectory of those insect wings and the force history provided in the earlier studies [24,27], one can see a similar correspondence between the downward wing translation and the large lift production, which is likely caused by a similar drag-based effect.
3.4. Three-dimensional vortex structures
Figure 12 shows a few selected snapshots of the three-dimensional flow field in a stroke cycle, which is identified by plotting an isosurface of the imaginary part of the complex eigenvalue of the instantaneous velocity gradient tensor . This technique allows one to identify regions where rotation dominates over strain.
An LEV is developed in the beginning of the downstroke, and this LEV grows stronger and remains stably attached to the wing during most of the downstroke. During wing translation, the LEV, the tip vortex (TV) and the shed trailing-edge vortex (TEV) are connected end to end, forming a vortex loop, within which the air moves downward (figure 12a). Towards the end of downstroke, the wing rotates rapidly along its own axis, and the LEV is divided into two branches, known as dual LEV , as seen in figure 12a. Corresponding to the stable LEV, there is no clear lift drop throughout the downstroke translation. At the end of downstroke, the LEV starts to shed from the wing as seen in figure 12b. Another feature of the downstroke is that the wing catches the vortex loop produced by the preceding upstroke and disrupts this loop through the wing–wake interaction.
During the upstroke, an LEV is also formed in the beginning (figure 12c), but the distal portion of this LEV is pinched off during mid-upstroke, as shown in figure 12d. Correspondingly, there is a visible dip in the vertical force around the same time of the upstroke shown in figure 5a. Later, the LEV will be formed again and will also form branches like dual LEVs. As discussed earlier, during the upstroke the wing misses the wake produced by the preceding downstroke. As a result, the vortices generated by the downstroke are better preserved in the wake.
Figure 12 also shows that the wake contains many slender-shaped vortices. These vortices are formed mainly due to break-up of the TV and TEV at the current high Reynolds number. This flow behaviour is consistent with the result of a recent work  that demonstrated a similar phenomenon of vortex break-up at Re = 1500. To further confirm the accuracy of these vortices, we have compared the simulations from the baseline and the finest mesh as described in §2.2, and the results show that the general characteristics of the vortices are consistent. In the regions far away from the wing, the isolated vortices likely contain artificial effects due to reduced resolution there.
3.5. Full-body simulation
A full-body model of the hovering hummingbird is also created by using symmetric kinematics for the left and right wings. The body of the bird is approximated by a sequence of ellipses with different sizes and aspect ratios. The bird model is run in an extended domain in the Y-direction. The typical flow field is shown in figure 13 for mid-downstroke and shortly after supination. From the vortex structures in the flow, we note that LEV and the TV during the downstroke are similar to those in the single-wing simulation. However, during supination, the two wings are near each other (the included angle is about 30°). The flows around the two wings are close enough to interact. In particular, when the wings move away from each other, the vortices generated from each wing grow and collide with one another. The interaction is complex and leads to further break-up of the vortices (see a movie in the electronic supplementary material). Other than that, the major vortex structures, such as the LEV and TV, are similar to those seen in the single-wing simulation.
Despite the effect of the wing–wing interaction on the three-dimensional vortex structures, the lift production is not significantly affected. Figure 14 provides a comparison of the lift coefficient between the full-body and the single-wing simulations. It can be seen that the forces from the two simulations are very close to each other. This result suggests that the wing–wing interaction and the wing–body interaction do not play an important role in lift production of the hummingbird.
A three-dimensional simulation was performed for a hovering hummingbird with the realistic wing motion reconstructed from imaging data. The simulation captures the lift and power characteristics in a stroke cycle and also details of the flow field. Our result confirms and provides specific data for the lift asymmetry that was previously suggested based on the measurement of the circulation around the wing. Furthermore, we quantitatively analysed the sources of the lift asymmetry and pointed out the mechanisms that lead to the asymmetry. Summarizing the results, the downstroke produces 150% higher vertical force than the upstroke. Among the factors, the wing area contributes 10% greater force, the drag-based effect contributes 60%, the wing–wake interaction contributes 30% and the remaining 50% can be attributed to the combined wing speed, angle of attack and wing rotation.
This research was supported by the NSF under CBET-0954381 to H.L. and IOS-0920358 to T.L.H. The computing resources were provided by the NSF XSEDE and the Vanderbilt ACCRE.
- Received May 21, 2014.
- Accepted June 18, 2014.
- © 2014 The Author(s) Published by the Royal Society. All rights reserved.