Three-dimensional simulation for fast forward flight of a calliope hummingbird

We present a computational study of flapping-wing aerodynamics of a calliope hummingbird (Selasphorus calliope) during fast forward flight. Three-dimensional wing kinematics were incorporated into the model by extracting time-dependent wing position from high-speed videos of the bird flying in a wind tunnel at 8.3 m s−1. The advance ratio, i.e. the ratio between flight speed and average wingtip speed, is around one. An immersed-boundary method was used to simulate flow around the wings and bird body. The result shows that both downstroke and upstroke in a wingbeat cycle produce significant thrust for the bird to overcome drag on the body, and such thrust production comes at price of negative lift induced during upstroke. This feature might be shared with bats, while being distinct from insects and other birds, including closely related swifts.


Introduction
Hummingbirds are distinguished and extremely agile flyers among birds. They are capable of not only sustained hovering flight, but also fast forward flight and various rapid manoeuvers. Recent studies of fluid dynamics have mainly focused on hummingbirds' unique hovering capability and unsteady aerodynamics associated with the wings that move in a relatively horizontal plane [1][2][3][4][5][6].

Model configuration and simulation approach 2.1. Reconstruction of the wing kinematics
A female calliope hummingbird (Selasphorus calliope) was the study subject, whose basic morphological data are provided in table 1. The experimental study was conducted to obtain the wing kinematics at a sustained flight speed of U = 8.3 m s −1 , at which the wingbeat frequency is 45.5 Hz while the stroke plane angle between the stroke plane and the horizontal is 67.9 • . In the experiment, the bird was placed in an open-circuit, variable-speed wind tunnel with a feeder at the middle of the tunnel, and it was trained to adapt to the wind while feeding. We recorded flight kinematics of the hummingbird using three highspeed cameras distributed outside the wind tunnel: 1 Photron SA-3 (Photron USA Inc., San Diego, CA, USA) and 2 two Photron 1024 PCI, all with electronically synchronized shutter timing. Two cameras were placed dorsally to the bird, and one was placed laterally, as shown by the views in figure 1. Video recordings were made at 1000 Hz with a shutter speed of 1/10 000 s. The working section of the tunnel is 85 cm in length, square in cross section, 60 × 60 cm 2 at the inlet and increasing to 61.5 × 61.5 cm 2 at the outlet to accommodate boundary-layer thickening [28]. Maximum deviations in velocity within a cross section are less than 10% of the mean; the boundary layer is less than 1 cm thick and turbulence is 1.2%. After the videos were taken, a custom MATLAB program [29] was used to track the markers frame by frame and to extract their three-dimensional coordinates. These markers were pre-labelled on the wings using non-toxic paint, and they included five points on the leading edge, on at the wingtip and three on the trailing edge, as shown in figure 2, where a comparison of the reconstructed model with the camera view shows that the instantaneous wing position and deformation are well captured. In a similar study for the hovering hummingbird [5], a principal components analysis verifies that these points are sufficient to characterize the wing motion. The wing geometry reconstruction process is similar to that in    the previous study [5]. The wing profile was constructed using spline interpolation through the marker points, and the wing surface was then built using triangular elements within the profile. The bird body was reconstructed using the camera views of the hummingbird. In the current reconstruction, a single wing consists of 1335 elements and 718 nodes, while the body surface consists of 3560 elements and 1782 nodes. A total of 13 cycles of wingbeats during steady flight were captured, and each cycle contains approximated 22 frames. To increase the time resolution of the wing position, the trajectories of the wing mesh nodes are also refined by spline interpolation in time. Seven cycles of wingbeats were reconstructed from the imaging data and used for the simulation. Figure 3 shows a sequence of wing positions within a cycle and also the wingtip trajectory (see the electronic supplementary material for an animation).   between the chord and the flight direction. The angle of attack, α, is defined as the angle between the chord and the relative flow direction that combines both the freestream velocity and the translational velocity of the chord at the leading edge. These two angles are plotted in figure 5 for two chords and five cycles, one proximal chord at dimensionless locationr = r/R = 0.15 and one distal chord atr = 0.9, which are denoted by subscripts p and d, respectively. From these plots, we can see large differences between the proximal chord and distal chord. For the proximal chord, the chord angle ψ p and angle of attack α p are both positive during the entire cycle. For the distal chord, these angles change the sign and vary significantly. During the downstroke, ψ d is negative, i.e. the leading edge tilting downward, but α d is positive owing to fast translation of the chord. During upstroke, ψ d is positive, i.e. the leading edge tilting upward, but the α d is negative, indicating that the pressure surface and the suction surface are swapped at that moment. Wing twist can be described by the difference between the two chord angles, ψ d − ψ p , which is plotted in figure 6. It is shown that rsos.royalsocietypublishing.org R. Soc. open sci. the twist angle reaches its extreme value during mid-downstroke and mid-upstroke; however, it is more pronounced during upstroke (near 40 • ) than during downstroke (near 25 • ). These differences between the proximal section and the distal section lead to highly non-uniform pressure distribution on the wing surface as shown later.

Simulation set-up and verification
In the model, the Reynolds number, defined as Uc/ν, is set to be Re = 3000, wherec is the average chord length and ν is the kinematic viscosity. The flow is assumed to be governed by the viscous incompressible Navier-Stokes equation, which is solved by an in-house code that adopts a secondorder immersed-boundary finite-difference method. The code is able to handle large displacement of the moving boundaries [30]. A fixed, non-uniform, single-block Cartesian grid is employed to discretize the domain  Table 2. Comparison of the force coefficients for the wings and body from the two different meshes, where C Z and C T are for the vertical force and thrust of one wing, respectively, and C Z,b and C D,b are for the vertical force and drag of the body, respectively. million) points are used for the baseline simulation. A finer mesh is also used in the simulation to verify grid convergence. Both of these meshes have maximum resolution around the wing, which is 1 60 cm in all three directions for the baseline case and 1 70 cm for the refined case. The simulation was run in parallel using domain decomposition and Message Passing Interface (MPI). The time step is t = 5 µs, which leads to approximately 4400 steps per wingbeat cycle. A multigrid method was employed to accelerate convergence of the Poisson solver. A total of 96 processor cores were used for the baseline case, and 128 cores for the refined case.
The vertical force F Z and thrust F T = −F X generated by one wing are normalized by the fluid density, ρ, the flight speed, U, and the surface area of the wing according to The lift and drag on the bird body, F Z b and F D b , are normalized in the same manner. The aerodynamic power coefficient of one wing is defined as where f is the stress on the wing surface, and u is the velocity of a point on the wing in the body-fixed coordinate system. The simulation results for two wingbeat cycles from both meshes are shown in figure 7 and table 2 for comparison. In figure 7 and also other figures from herein, the shaded area indicates downstroke, while the white area indicates upstroke. These results include the time-averaged lift and thrust of one wing and also lift and drag of the bird body. From the table, we see that the maximum difference of all the forces is less than 5%. Thus, the baseline resolution is deemed satisfactory for the current study.

Aerodynamic forces
The force and power coefficients are shown in figure 8, which include both instantaneous and phase-averaged data. The cycle-averaged data are listed in table 3 for an entire cycle and for downstroke/upstroke separately. Figure 8a shows that the weight support is mostly generated during downstroke where C Z is positive. Mid-downstroke corresponds to the maximum lift production. During supination and early upstroke, the wings are still able to generate some weight support. Around mid-upstroke, vertical force becomes negative even though its amplitude is not particularly high.
On the other hand, figure 8b shows that thrust is mostly positive during both downstroke and upstroke. Furthermore, thrust has a greater peak during upstroke than during downstroke. However, the data in table 3 show that downstroke on average produces more thrust than upstroke. Figure 8c shows that the power consumption during both half wingbeats are significant. However, the power requirement is greater for downstroke, about twice as high as upstroke. This feature is similar to hovering, where downstroke power is nearly 2.8 times of upstroke power according to Song et al. [5], who studied the ruby-throated hummingbird. Using the equation P = 1 2 C P ρU 3 (2S) for the aerodynamic power, P, we have P = 94.5 mW for the calliope hummingbird, and body mass-specific power 34 W kg −1 . Thus, mass-specific power output of the hummingbird is within the range reported for larger bird species. For example, cockatiel power output ranges from 17 W kg −1 at 5 m s −1 to 47 W kg −1 at 14 m s −1 , and dove power output ranges from 31 W kg −1 at 7 m s −1 to 54 W kg −1 at 17 m s −1 [31]. Compared with the hovering ruby-throated hummingbird at 55 W kg −1 [5], forward flight in the calliope hummingbird requires less power, which is expected since hovering generally is more energy-demanding than forward  Table 3. Force production and power of upstroke, downstroke and entire cycle (i.e. average between upstroke and downstroke). flight. On the other hand, forward flight at 8.3 m s −1 is not necessarily minimum power speed for the hummingbird, as the mechanical power output of birds can be described as a U-shaped curve function of the flight speed [31][32][33].
Using the present force coefficients and equation F total = 1 2 (2C Z + C Z,b )ρU 2 S, we obtain the total vertical lift produced by the bird, which is around 94% of the bird weight. The bird body itself generates about 22.2% of body weight. This result will be discussed later. The thrust generated on the two wings together is 152% of the body drag. The imbalance of the vertical and horizontal forces could have been caused by several reasons: (i) the absence of camber in the wing model, (ii) error in digitization of the wing position, (iii) beak-feeder interaction as the bird was attempting to feed during the recording,  (iv) interaction between the bird's body and the wake of the feeder, and (v) an underestimate of drag coefficient for the body.

Force production mechanism
Overall force production of the hummingbird can be explained from the wing kinematics as viewed from the global coordinate system, i.e. the coordinate system fixed with the ambient air. Figure 9 shows the proximal and distal chord moving in the global coordinate system with their trajectories traced out. During downstroke, the angle of attack is positive for both the proximal chord and the distal chord. Therefore, both wing sections generate weight support. Since the leading edge of the distal section tilts downward, the aerodynamic lift has a forward component that leads to thrust generation during downstroke.
During upstroke, both wing sections move forward in air, even though the stroke plane angle is less than 90 • and the wings move backward with respect to the body. Nevertheless, positive thrust is generated during this half cycle by the distal section. As shown in figure 9, the angle of attack of the distal chord is negative at upstroke, and the overall force on the section points downward and forward. Figure 10 shows the pressure distribution within four selected vertical slices at mid-downstroke and mid-upstroke. It can be seen that the roles of the distal section and proximal section are different. For both downstroke and upstroke, the proximal wing has pressure surface on the ventral side and suction surface on the dorsal side. Thus, its main role is for vertical force generation. However, the distal wing flips its angle of attack between the two half cycles. Thus, positive (negative) pressure is distributed on the ventral (dorsal) side during downstroke, and the opposite is true during upstroke. This pressure differential leads to weight support during downstroke only, but thrust production during both downstroke and upstroke.

Vortex structures
Vortex structures, which are induced by the wing motion and dominate the wake, have been a focal point in the study of force production of flapping wings and fish fins. They can also be used to evaluate whether a bird adopts slow gait or fast gait [22]. As has been pointed out by previous researchers, at slow gait the trailing-edge vortices (TEVs) form rings after each downstroke, and the sequence would look like a series of smoke rings [34,35]; while at fast gait, the tip vortices (TVs) form undulating vortex tubes from the tip, and the TEVs form cylinders from the trailing edges, both being convected downstream [24,25,36].
In the current study of hummingbird flight, we used the iso-surface to show the vortex structures. The scalar quantity is defined as the maximum value of the imaginary part of the eigenvalue of the velocity gradient tensor ∇u, and it describes the strength of the local rotation of fluid [37]. It is found that the TVs are continuously shed from the wingtip, and the TEVs shed with the shape of separate cylinders. Such ladder-type vortex structures were also observed in the measurement of forward flight of swifts [25], which are closely related to hummingbirds, and like hummingbirds, do not flex their wings during flight. As pointed out by the authors [25], the ladder-type wake is formed by continuous shedding of the spanwise vortices from the wing surface and is different from an earlier speculation that for swifts and hummingbirds, a ladder-like wake would be generated by a distinct vortex in each of downstroke and upstroke [38]. Even though they both generate continuous ladder-like vortex structures, swifts do not forego weight support on upstroke like we report for hummingbirds.
Several snapshots of the flow field are shown in figure 11. These snapshots show roughly the shape of the TVs that follow the trajectory of the wingtips. In addition, vortex shedding from the trailing edge is evident. Formation of the leading-edge vortices (LEVs) during both downstroke and upstroke is visible, and the LEVs are stable for both downstroke and upstroke. The behaviour of the LEVs has to do with both the instantaneous angle of attack and pitching rotation of the wing [39]. From figure 5d, the angle of attack of the distal section keeps a maximum value around 25 • for a significant period of time, which would have caused LEV shedding and stall, if the wing simply translated without changing its pitch. However, the wings are also performing rapid pitching around their axes, as seen from variation of the chord angle plotted in figure 5b. That is, the chord angle magnitude decreases during downstroke after t/T > 0.2, and quickly increases during upstroke before mid-upstroke. Such rotational motion has been known to maintain stability of LEVs and to enhance lift production of the wings [39].

Forces on the bird body
Lift and drag on the bird body are affected by the orientation of the bird during flight. In general, the inclination angle of the body decreases with the increase of flight speed [12,23,27,40,41]. In the current study, the body angle of the hummingbird is χ b = 12 • , which is close to the angle of the rufous hummingbird at a speed of 8 m s −1 , where χ b = 11 • [12]. Figure 12 shows both the instantaneous and phase-averaged data of the forces on the body. Downstroke-, upstroke-and cycle-averaged data are listed in table 3. These results show that the lift on the body provides 22.2% of the weight support. Furthermore, lift on the body during downstroke is 1.76 times of that during upstroke. In figure 12a, lift on the body oscillates significantly during a wingbeat cycle. On the other hand, drag on the body does not vary very much in a cycle and is nearly equal on average between downstroke and upstroke, which are reflected in figure 12b and table 3.
The high percentage of lift produced by the body and the oscillations of the body lift in a cycle may be attributed to aerodynamic interaction between the wings and the body that is assumed to be stationary in the current study. To verify this possibility, we also simulated separately the same flow around the isolated body without the wings attached. The shape and orientation of the body remain the same in the test. Figure 13a shows the pressure distribution on the bird body from the isolated body simulation, which can be compared with the result from the full-body simulation shown in figure 13b and c for middownstroke and mid-upstroke, respectively. For the isolated body, even though a high-pressure zone     is established below the body and near the head, the flow passes around the body in the absence of the wings and merges behind and above the body, where the pressure is partially recovered. As a result, the overall lift by the body is small. When the wings are present, the flow from below is prevented from passing around the body by the wings. Furthermore, the wing-wing interaction mechanism, similar to that proposed by Lehmann et al. [42], apparently has played a role here. That is, when the two wings are separating from one another from pronation to mid-downstroke, they create a low-pressure zone above the bird body, as shown in figure 13b, thus leading to a net upward force. This mechanism also explains why during upstroke, the low-pressure zone above the body, as plotted in figure 13c, is significantly smaller when compared with downstroke. The present result shows that the bird body has significant contribution to the overall weight support. In comparison, previous experimental studies of insects and other birds indicated that the body lift is only a small portion of the animal weight, as shown in table 4. However, we point out that in those previous studies, the force was measured for the isolated animal body only, while in the current study, the wings are present and are in constant motion. For the isolated hummingbird body, we also observed  [48], bumblebee: [15,17], rufous hummingbird: [12], magpie, pigeon and zebra finch: [40,41], hawk moth: [49], barn swallow: [50], bat: [51]). low lift production. As shown in table 4, lift of the isolated hummingbird body is only 8% of the weight and is comparable with previous data for insects and also birds (e.g. 15-20% during flexed-wing bounds in zebra finch).

Comparison of hummingbirds with other flying animals
The advance ratio J and stroke plane angle β are two primary factors that affect the force production of flapping wings during forward flight. These two variables differ largely among animal species. Figure 14 shows a few species on the β − J map with the data directly collected or derived from various sources. It can be seen that hummingbirds largely fall within the range of the insects but also extend into the range of other birds. For all species, the stroke plane angle increases with the advance ratio, which is expected since at fast flight speed, the animals not only reduce the body angle for drag reduction, which would naturally cause the stroke plane angle to increase, but may also tilt the stroke plane further to enhance thrust generation. For the small insects like bumblebees and fruit flies, the advance ratio is usually less than one [14]. For such slow flight, lift production is predominant over thrust production. Since the back-sweeping velocity of the distal wing exceeds the forward flight speed at upstroke [15,20,33], the wingtip trajectory traced out in the global coordinate system is highly backward skewed at upstroke, which is shown in figure 15 for a bumblebee at J = 0.6. In this case, downstroke is mainly for lift production, while upstroke is mainly for thrust production. If the flight speed is further reduced, with a more skewed trajectory, upstroke may even produce lift as well. Overall, this strategy of using upstroke is also known as 'backward-flick' [18,23,31]. An exception is the fruit fly which is shown by a recent study that its upstroke uses a paddling mode to produce drag-based thrust [20].
For most birds, the advance ratio is significantly greater than one, and the stroke plane angle is close to 90 • . Thus, the wingtip trajectory in the global coordinate system becomes much less backward-skewed and the wavelength-to-amplitude ratio becomes greater as shown in figure 15 for a pigeon. In this case, upstroke is not suitable for thrust generation. Instead, the wings are either feathered with little force produced [24,25], or swept on upstroke with some lift production (e.g. pigeon) [26,27]. On the other hand, a powerful downstroke is used to produce both lift and thrust.
For the hummingbird in the current study, the advance ratio is between that of insects and other birds. The wingtip trajectory is moderately skewed as shown in figure 15. Therefore, with proper angle of attack, the wings can still produce thrust during upstroke. However, since the overall force points downward, some lift has to be sacrificed. From this figure, it can be seen that thrust can be produced when the wing speed at upstroke is comparable to or possibly even lower than the flight speed. This thrust mechanism is analogous to a sail that moves against wind and thus is termed a 'sail mode' in this work. On the other hand, downstroke of the hummingbird is similar to that of big birds, as shown in figure 15, where both lift and thrust are generated.  It should be pointed out that at slow flight speeds, force production of the hummingbird still appears close to that of insects. As shown by Tobalske et al. [12], when J is below 0.7, the stroke plane angle of the hummingbird is small and the wingtip trajectory is also highly skewed like that of insects. Similarly, some insects can also perform fast flight at J > 1, e.g. hawk moth, as seen from figure 14. It would be interesting to see whether their force production mechanism is similar to that described here for the hummingbird.
There is significant similarity between the hummingbird and bats, both being capable of hovering and forward flight at J = 1 [52]. Recent flow visualization studies of bats have suggested that aerodynamic function of upstroke changes with forward flight speed [16]. At hovering and slow speeds, bat wings are inverted during upstroke to aid weight support; and at fast speeds, downstroke produces weight support and thrust, but upstroke may generate extra thrust at cost of negative lift. These features are similar to what we have reported for the hummingbird, which is interesting given that, unlike hummingbirds, bats flex their wings during upstroke.
It is also interesting to point out the similarities and differences between the hummingbird and the swift, as both of them keep their wings extended during the entire wingbeat cycle regardless of flight speeds. This is highly unusual for birds, as most birds flex the wings during upstroke [53]. Regardless of keeping the wings extended, the combination of spanwise twist and the stroke plane angle reported here for the hummingbird is a mechanism for permitting flight over a wide range of speed. For the swift, it is believed that wing twist is an important factor leading to optimal efficiency of flapping flight [54]. On the other hand, the swift's flight style is dedicated to cruising flight, and it is not adept at hovering, and some swift species are even unable to accelerate into flight from a standing start without first dropping in altitude to gain velocity [55]. Therefore, one conclusion could be that the swift pattern of force development is optimized for more continuous weight support during cruising flight in a way that the closely related hummingbird does not achieve; while for the hummingbird, its evolutionary trajectory may have favoured the capacity for some thrust production at the expense of the more-constant weight support style in the swift [54]. This hypothesis merits testing in a broader phylogenetic context, as it is not appropriate to infer evolutionary trajectories from two-species comparisons [56].

Conclusion
Three-dimensional computer simulation has been performed to study aerodynamics of a hummingbird in fast forward flight and whose wing motion was captured by filming the bird flight in a wind tunnel. The finding places hummingbirds in an interesting position relative to insects and large birds. At a speed of 8.3 m s −1 , the advance ratio of the hummingbird is between those of typical insects and large birds, and the simulation results show that the hummingbird uses a different strategy for lift and thrust production. In particular, its power downstroke generates both weight support and thrust just like other birds, but its upstroke further enhances thrust by setting the distal wings at a proper angle of attack with respect to the oncoming air, even though such thrust enhancement comes at cost of some negative lift at upstroke. These features are thus similar to those of bats flying at J = 1. Ultimately, caution is necessary in interpreting our results as they emanate from the study of a single bird of one species. New studies in a phylogenetic context will be useful for understanding the evolutionary trajectories and selective pressures that drove the hummingbird flight style.
Ethics. All protocols associated with hummingbird care and experimentation were approved by the University of Data accessibility. Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.8ch1b.