Abstract
To deepen our understanding of the aerodynamics by which hummingbirds use flexible wings to hover efficiently during flapping flight, a three-dimensional wing model with dynamic morphing was developed according to the morphological and kinematic data of a hovering hummingbird’s wing. Navier-Stokes equations were solved on a dynamically deforming grid to study wing aerodynamics. In numerical simulations, boundary-based smoothing and overset methodologies were used in combination to update the interior nodes of cells in the computational domain, allowing those nodes to accommodate the motion of the flexible wall. This study showed that the leading edge vortex (LEV) attached to the wing was stable during the downstroke but extremely unstable and shed continuously during the upstroke. In the results of the downward stroke, the different vortices separated from the surface and formed a vortex ring. The difference is that the leading edge vortex induced a vortex ring near the root and a smaller and weaker vortex ring near the wingtip during the upstroke. A significant enhancement in aerodynamic forces was found during the downstroke, along with a large number of power consumptions. We found that the asymmetry in the time-averaged vertical force between the two half strokes was 3.5, and the value was higher than that reported earlier. Aerodynamic force coefficients and efficiency matched well with those of another hummingbird wing (the ruby-throated hummingbird). It is worth noting that hummingbirds can maintain a similar wingtip speed by flapping their wings, but different strategies are adapted to hover efficiently due to the differences in size and body weight. The results of this study help to better understand the aerodynamics of the hovering hummingbird.
1. Introduction
Hummingbirds, the smallest bird species with excellent flight capability, can perform both regular cruise flights and aerobatic flights [1]. In particular, they can perform sustained hovering. The wings of hovering hummingbirds are usually rotated around their long axis during supination (transition from downstroke to upstroke) and kept outstretched during the upstroke [2]. These behavioral characteristics are more similar to the flight characteristics of insects. Insects usually achieve the inversion of wings through muscle activation at the wing root and passive deformation of the wings [3]. Unlike insects, hummingbirds can reverse the angle of attack through active wing rotation at the wrist. And the dynamic wing morphing of the hummingbird wings is more remarkable than the passive deformations observed on smaller insect wings. Because hummingbird wings are anatomically and physiologically different from insect wings, hummingbirds generally do not fly in a similar way as insects do. Scientists and biologists have focused on the hovering flight of insects and hummingbirds for a long time in the past. Compared to tiny insects, giant hummingbirds are one of the few extensively studied vertebrate species.
Over the last decade, there have been many experimental studies on the underlying flight mechanism of hummingbirds. Previous studies primarily employed experimental measurements to examine the evolution of vortices in the flow field of hummingbird flight, which explained the unstable effect. Warrick et al. [4] measured the velocity field around a rufous hummingbird wing by conducting particle image velocimetry (PIV) experiments and discovered that a stable and conical leading edge vortex had a dominant role in the lift enhancement, especially during the downstroke. In another study [5], they found two vortex rings in the region swept by a pair of hummingbird wings. Wolf et al. [6] analyzed the strength of the shedding vortex of a hummingbird of the same species. They concluded that there existed a high degree of asymmetry in lift production based on the results of circulation. Some researchers tried to explain the mechanism of aerodynamic force generation using fluid dynamics theory by analyzing wing kinematics. Altshuler et al. [7] proposed that the flapping speed and the angle of attack during the downstroke were higher than those during the upstroke. In addition, they found a longer wingspan and positive camber during the downstroke. So, the asymmetries in force generation were attributed to motion asymmetries. However, the experimental results were not complete though we have gained much understanding of the aerodynamics of hovering hummingbirds thanks to the efforts of previous researchers. Based on the current measurement techniques, it is still challenging to obtain the details of the pressure distribution on the wing surface during the flight of hummingbirds. And it is not easy to make reliable predictions of transient aerodynamic forces using unsteady aerodynamics theory. Nafi et al. [8] estimated the aerodynamic forces on freely flying birds by analyzing the wingbeat kinematics and near-wake flow measurements using long-duration time-resolved particle image velocimetry. Lentink et al. [9] developed a method to measure aerodynamic forces on flying animals and biomimetic robots. However, they obtained their measurements by tying the animals and robots, which prevented the animals from moving freely.
Computational fluid dynamics (CFD) has been widely applied to the study of the aerodynamics of flapping wings because it has excellent advantages in predicting the transient aerodynamic forces and visualizing the flow field around flapping wings. Using simplified hummingbird wing models, some work on the aerodynamics of hovering hummingbirds captured the asymmetry of aerodynamic forces as well. Liang and Dong [10] modelled a hummingbird wing as a rigid body in their numerical simulations. They found that the amount of lift during the downstroke was 2.95 times that that during the upstroke. Song and Zhang [11] also conducted a relevant numerical simulation and neglected some details in the numerical model, including wing flexibility, camber, thickness, microstructure, and stroke plane deviation. Ignoring those that can influence the evolution of the flow structure could result in a mismatch of aerodynamic forces with those on realistic wings. However, the conclusions from previous CFD analyses might be partial due to the oversimplification of the kinematics and morphology of the hummingbird wings in those models. To produce a realistic simulation, an appealing methodology is to determine the wing motion and deformation and then use the wing shape captured as an input to the numerical solver. Using this methodology, Song et al. [12] conducted a three-dimensional numerical study on a hovering hummingbird wing (ruby-throated). As reported in their study, the asymmetry in lift production was 2.5. They also investigated possible reasons for the asymmetry from the aspects of wing kinetics and wing wake interaction. Although Song et al. [12] adopted actual kinematics in their numerical model, the wing chord was modelled approximately as straight segments without considering the cambers in cross sections in the process of reconstructing the kinematics and morphology. However, Maeda et al. [2] pointed out that one of the most prominent features of the cross sections was the temporal variation in the cambers. Altshuler et al. [13] found that realistic hummingbird wings (with sharpened leading edges and with substantial camber) produced higher lift than artificial wings, and some geometric details such as the presence of camber tended to improve the lift force. Maybe the variation in camber stems from differences in hovering strategies was adopted by hummingbirds of different species and individuals. To discover the general aerodynamic mechanisms of hovering hummingbirds, we should choose as many species-specific hummingbirds as possible and analyze their similarities and differences.
In this paper, we reconstructed a flapping wing with dynamic morphing using morphological and kinematic data from a realistic hummingbird (Amazilia amazilia). Numerical simulations were performed to quantify the aerodynamic forces and power. In the numerical simulations, we combined boundary-based smoothing and overset methodology to update the nodes of cells in the computational domain. The time-varying position of each node on the surface mesh was prescribed directly without considering complex fluid-structure interaction modelling. A self-developed program (User-Defined Function) was used to impose the wing motions and deformations. The validation of the model showed the effectiveness of this methodology. To explore some particular mechanisms of aerodynamic enhancement and better understand the aerodynamics of hummingbird hovering, we analyzed the time courses of aerodynamic forces, the characteristics of the flow field, and the evolution of vortex structures based on the simulation results, revealing the aerodynamic characteristics for efficient hovering of hummingbirds. Comparing the results of this study with those of a previous study on another hummingbird enabled a better understanding of the aerodynamics of hovering hummingbirds.
2. Materials and Methods
2.1. Morphological Modelling and Reconstruction of the Wing Kinematics
The object of study in this article is an Amazilia hummingbird in hover, and Maeda et al. [2] recorded the kinematics of its right wing in the Tama Zoological Park in Tokyo, Japan. Neither the gender nor the age of the hummingbird was known, but its age was estimated to be at least 10 years. The process of this experiment was briefly described here, and see their paper for details on that. In the experiment, the hummingbird frequently visited a nectar provided by a feeder (Figure 1(a)) and was thus able to be recorded at 2000 frames per second by four high-speed cameras. The outline and boundaries between various parts of the wing were manually traced and extracted as points with equal spacing, as shown in Figure 1(c). After the three-dimensional calibration coefficients for stereo reconstruction through a DLT (direct linear transformation) program were available, a process was repeated to find the best match between tracked curves and projected curves using camera coordinates. As a result, 17 discrete grids consisting of points for each became available.
(a)
(b)
(c)
The three-dimensional coordinates of those points in each frame could be obtained. All original data used in this study can be found in the electronic supplementary material of the reference [2]. In addition, the wing kinematics and dynamic wing morphing can be visually observed in the videos they shared. As you can see in Figure 1(c), the original points in each discrete grid make it easy to reflect the detailed features of the kinematics and morphology of various parts of the wing. However, irregular changes in some details, especially the gaps and slots caused by feather slip, might bring challenges to the dynamic update of nodes on the surface mesh and the deformation of cells around the wing in numerical simulations. Besides, large numbers of points might result in low computational efficiency when constructing numerical models through interpolation using such existing data sets. It is therefore necessary to filter some noisy points and take samples along the chordwise and spanwise. In a previous numerical study [11], nine markers located on the outline of a hummingbird wing, including five points on the leading edge, one at the wingtip, and three on the trailing edge, were used to characterize the wing motion. Tay et al. [14] investigated the possibility of shortening the digitization process by using fewer points or fewer frames per cycle and still maintaining a reasonable level of accuracy in the numerical results. Their results showed that it was possible to shorten and simplify the digitization process of a flapping wing by decreasing suitably the wing points and the number of frames per cycle. However, there were at most 16 markers in their tests, and they believed that is enough to model the motions of a flapping wing. Therefore, we finally selected sample points to reconstruct the kinematics and morphology of the wing. To model the spanwise bending and twist, there were 11 and 6 sample points with equal spacing extracted along the wingspan and chord, respectively. Thin-plate spline interpolation [15] is a robust spatial data interpolation method that is capable of approximating almost all bio-related deformations. Therefore, the spatial interpolation method, one of the radial basis function (RBF) interpolation methods, was able to be used to predict the positions of other missing points based on the three-dimensional coordinates of the sample points. A comparison between the experimental model and the reconstructed model at two random time instants in a stroke period is shown in Figure 2. Figures 2(a) and 2(b) correspond to a certain time instant in the downstroke and one in the upstroke, respectively.
(a)
(b)
With the discretization of the fluid domain finished, the three-dimensional coordinates of each node on the wing surface mesh were available. The initial wing configuration was used as a reference for temporal interpolation. We refined the trajectory of each selected sample point using RBF interpolation in order to obtain the small-step solution of the numerical simulation. Using interpolation, we were able to reconstruct a model that passed through the selected sample points accurately. As the model includes various dynamic wing morphing, their aerodynamic effects are incorporated into the simulation results. Note that the corrugations caused by the feathers and the small holes in the feather vanes were ignored, which could affect the boundary layer and thus the generation of aerodynamic forces. A study that uses kinematic data from different spanwise positions to quantify dynamic wing morphing has been published [2]. The morphological and kinematic parameters of the wing are listed in Table 1, and those values were used as reference values in numerical simulations. Figures 3(a) and 3(b) show the kinematics of the wing during the downstroke and upstroke, respectively. For the real hummingbird hovering, the wing length and wing area vary with time in a flapping cycle. The maximum value and minimum values of wing length and wing surface area in one stroke period are also listed in Table 1.
(a)
(b)
2.2. Setup of the Flow Solver and Model Validation
In the present work, computational fluid dynamics (CFD) simulations were performed using the commercial software ANSYS Fluent (Fluent 19.2 package). The approach, the finite volume method, was the same as that adopted in Reference [16] and has been shown to be valid for studying the unstable aerodynamics of flapping wings [17]. Young and Lai [18] used the commercial CFD solver to estimate the force generation of a flapping wing at Reynolds numbers ranging from 100 to 50000, and the results had no significant difference whether a laminar model or a turbulent model was adopted at the same Reynolds number. In addition, similar CFD with a laminar model was used in previous studies [10, 11]. For hummingbirds hovering in still air, the Reynolds numbers (denoted by Re) root in the order of . As a result, the flow field around hovering hummingbirds was assumed to be laminar, without considering the effects of turbulence in this realm. Consequently, at low Reynold numbers, air flow was believed to be governed by the three-dimensional viscous incompressible, unsteady continuity, and Navier-Stokes equations for laminar flow, as shown in their conservation forms. where is air velocity, is the pressure, is the stress vector, and is gravity.
Here, a rigid three-dimensional flapping wing was selected from a study [16] as the object in this validation. The wing has a semielliptical planform with semimajor and semiminor axes of and , respectively. An offset of 113 mm is present between the root and the center of the flapping axis. As a result, the wing has a mean chord () of 77 mm, an area of the plane of 0.0193 m2, and a second moment of area radius () of 229 mm. The wing with a thickness of 1.5 mm was assumed to be rigid and was confined to a horizontal stroke plane. Therefore, its motions were defined by periodic variations of the flapping angle and rotation angles. The flapping follows a simple harmonic profile (see Figure 4(a)), while the time courses of rotation angles consist of invariant sections joined by sinusoidal curves (see Figure 4(b)). The harmonic and trapezoidal motions were defined as follows: where and represent the instantaneous harmonic flapping angle and trapezoidal rotation angle, respectively. The normalized time is defined by , in which is the time instant and is the flapping period. is the midstroke angle of attack, and is the normalized duration of the rotation motion. Note that the downstrokes and upstrokes are symmetrical. A test set in series 1 was selected for comparison. In this case, the midstroke angle of attack and the normalized duration of the rotation motion correspond to and 0.3, respectively. We adopted the definition of , in which is the distance travelled by a point at the position of the wing per flapping cycle and is the kinematic viscosity of the fluid medium. The numerical study was conducted at , which falls within the Reynolds numbers of hovering hummingbirds. In this study, if the average wingtip velocity had been specified as the reference velocity, the value of Re would be 6344.
(a)
(b)
In this numerical simulation, the second-order upwind scheme and the central difference scheme were used to discretize the convective and diffusive terms, respectively. We obtained the pressure and velocity solutions using the coupled algorithm. A User-Defined Function (UDF) was used to realize the predefined flapping motion described in Equations (3) and (4). A macro–DEFINE_GRID_MOTIO–allows us to specify the three-dimensional coordinates for each node on the surface mesh to update their positions at each time step. We developed a program in the preprocessing to achieve the updates of the nodes on the rigid wall at any time instant by using the rotation transformation matric and the node coordinate at the starting time instant. The computation domain is a hemisphere with a diameter of , and it was discretized into tetrahedral fluid cells using ANSYS Meshing software. The grid consists of two parts, the background mesh representing the far-field region and the separate component mesh around the wing. A no-slip wall condition was applied on the wing surface and the outer boundary of the computational domain. There was only one simulated wing instead of a pair of identical wings, so the flat plane of the sphere was set to a symmetric boundary condition. The outer boundary of the component mesh was modified to the overset type, and an overset interface was created to connect the two fluid domains. The boundary-based smoothing method and the overset methodology were combined to make the interior nodes accommodate the movement of the wings.
For diffusion-based smoothing, the mesh motion is governed by the diffusion equation where is the mesh displace velocity. The diffusion coefficient as a function of the boundary distance is of the form where is a normalized boundary distance and is a user input parameter. The Laplace equation (Equation (5)) describes how the prescribed boundary motion diffuses into the interior of the deforming mesh. With a nonuniform diffusion coefficient, mesh nodes in regions with high diffusivity tend to move together (that is, with less relative motion). The velocity of the nodes on the boundary mesh can be calculated by the solver according to the changes of the specified node coordinate at each time step and the time step size (). After this boundary condition was determined, Equation (5) was solved through the biconjugate stabilized (BCGSTAB) method with finite element discretization and was obtained directly at each mesh node. The node positions were then updated according to where and is the new position of the node at the next time step and the old position at the current time step, respectively. Therefore, the outer boundary of the component mesh was deformed by absorbing the motion of the inner boundary and reconnected to nonconformal interfaces at each time step, which connected the cell zones by interpolating the cell data in the overlapping regions. Finally, the mesh density and the time step size were selected on the basis of comparative studies. The time history curves of lift and aerodynamic power coefficients are shown in Figure 5, and the experimental and computational results from the reference [16] were also plotted for comparison. The variation of the lift and power coefficients is in good agreement with that of the previous study. A few inconsistencies in the results might be caused by the differences arising from the mesh topology and the grid update strategy. The sliding mesh feature in this solver was used to impose the appropriate wing kinematics in another research, and there might be some unavoidable measurement errors in the experimental results.
(a)
(b)
2.3. Simulation Set-Up and Model Verification
In this study, the definition of Reynolds number is consistent with those of previous studies on the aerodynamics of hovering hummingbirds. Based on the rigid body assumption, it was calculated to be where is the fluid density, is the dynamic viscosity of a fluid, is reference length (mean chord length ), is reference velocity (average wingtip velocity ), and is the peak-to-peak amplitude of the stroke angle. Air density is , and its viscosity is . The values of the other parameters mentioned above can be found in Table 1. Equations (1) and (2) were solved to simulate the reconstructed model of the hummingbird wing. The spatial and temporal discretization schemes, the solution method of the flow governing equations, the grid update strategy, and the boundary conditions were the same as those used in the previous validation. The solver will begin to read the interpolation data for each node from an external file and store them in memory as soon as a macro–DEFINE_ON_DEMAND–is executed in the self-developed program. To control the change in position of each node on the surface mesh independently, another macro, DEFINE_GRID_MOTION, was used to assign the three-dimensional coordinates for each node based on the interpolation data at each time step. The solver automatically updated the position of the wing at each time step as the computation progressed. There existed relative motion between the nodes on the surface mesh, so the current wing was modelled as a flexible wall involving fluid-structure interaction rather than a rigid wall. As shown in Figure 6, the tetrahedral cells in the overlapping zone were guaranteed to have similar sizes, generating the matching interface. The wing thickness is 0 mm in this model, and thus, the thickness effect has been ignored.
(a)
(b)
(c)
(d)
In order to select a larger enough quantity of fluid cells and a fine enough time step for the numerical simulations, several cases with a distinctive grid density and time step were tested in this study. Those grids had the same cell topology, but the sizes of the cells were scaled at the corresponding locations to conduct the grid-independent test. After the aerodynamic forces became available, they were normalized according to where is the vertical (-axis) or horizontal force (-axis) on the wing. As shown in Figure 7, the three grids consist of 0.7, 1.3, and 2.5 million tetrahedral cells, respectively. And they correspond to the coarse mesh, the medium mesh, and the fine mesh, respectively. The good convergences of the simulation results were obtained when the number of tetrahedral meshes reached 0.7 million. If the aerodynamic forces were our sole interests, then, the grid density would be adequate. However, we selected a denser grid, the grid with 1.3 million tetrahedral cells, to pay attention to the details of the vortex structures around the wing as much as possible. We conducted simulations on this grid at the time steps of and to verify that the results were independent of the selection of time steps. Considering the restriction of computing power, we believed that the time step size of was acceptable.
(a)
(b)
3. Results and Discussion
3.1. Force, Power, and Efficiency
As shown in Figure 6, , , and represent the three orthogonal axes of the global coordinate system and correspond to the backward, spanwise, and vertical directions, respectively. The horizontal and vertical components of the net aerodynamic force, and , are normalized by air density, the wing area, and the average wingtip velocity according to equation (9). Therefore, and represent the horizontal force coefficient and the vertical force coefficient, respectively. The aerodynamic power was calculated by integrating the dot product between the velocity vector and the stress (pressure and wall shear stress vector) over all surface meshes according to where is the static pressure, is the wall shear stress vector, is velocity vector, is a surface vector, is the unit normal vector of the surface, and is the magnitude of the surface vector . The power coefficient is defined by normalizing the power by .
Note that there exist slight cycle-to-cycle variations as you can see in Figure 8. We averaged these coefficients for the last three cycles, and the final results are shown in Figure 9. The time-averaged values of the force and power coefficients are defined as follows: where is the start time instant and is a time interval. The time-averaged values of the force and power coefficients generated during the downstroke, upstroke, and whole stroke cycle are listed in Table 2, respectively. Efficiency is defined as the ratio of the time-averaged vertical force coefficient to the time-averaged power coefficient and is thus .
(a)
(b)
(c)
(a)
(b)
(c)
Figure 9shows that the vertical force was highly asymmetric. A high peak was observed in the force coefficients during the second half phase of the downstroke. The magnitudes of the force coefficients are relatively smaller at most time instants. However, it can be seen that larger vertical forces were generated during the downstroke. In addition, the wing produced drag force during the downstroke, but thrust force during the upstroke. However, the magnitude of the horizontal force coefficient was lower than that of the vertical force coefficient during the entire stroke cycle. There were multiple maxima for the vertical force coefficient during the upstroke, and the power coefficient followed a similar trend to the vertical force coefficient. The wing could produce lift and thrust during the upstroke without a large number of power consumptions.
The time-averaged vertical forces generated during the downstroke and upstroke were 5.3 g and 1.5 g, respectively, causing an asymmetry of 3.5. The value in time during a stroke cycle was , corresponding to of the total weight support provided by the two wings. The time-averaged horizontal force coefficient was and much smaller than the value of . A positive value refers to a net drag force, while a negative value refers to a net thrust force. Power consumption during the downstroke is 5.90 times more than power consumption during the upstroke, and it is 0.21 W for the entire stroke cycle. There existed an improvement in the production of aerodynamic forces, resulting in higher power consumption during the downstroke. The aerodynamic efficiency during a stroke cycle was 1.17, and the wing had a higher efficiency during the upstroke.
3.2. Three-Dimensional Vortex Structures and Wake
According to the CFD results, the vortex structures were visualized using the isosurfaces of the nondimensional criterion [19]. The positive value is indicative of vortexes in the -criterion formulation. Here, the value was normalized according to
All three-dimensional vortex structures could be clearly shown by adopting an appropriate value of . The vortex structures around the wing were analyzed to understand their underlying effects on aerodynamics. Figures 10 and 11 show a few selected snapshots of the vortex structures near the wing surface at different time instants, respectively. And they were identified by plotting the isosurface of .
At the beginning of the downstroke, the dorsal side of the wing was characterized by the formation of the leading edge vortex (DLEV). The DLEV spiraled from root to tip and grew in size. However, the leading edge vortex (ULEV), tip vortex (UTV), and root vortex (URV) on the ventral side of the wing that had been developed during the previous cycle did not separate from the surface of the wing. Subsequently, the vortex near the wingtip began to move along the chord and soon separated from the wall. As the wing accelerated during the translation phase, the leading edge vortex, the tip vortex, and the trailing edge vortex on the dorsal side of the wing gradually developed and grew in size. At 0.5 T, the leading edge vortex in the outer segment, the tip vortex, and the trailing edge vortex separated from the wing surface, and they were eventually fused into the wake. During the earlier downstroke, the root vortex (URV) still adhered to the ventral side of the wing, which caused the formation of the dorsal root vortex (DRV) to be delayed. Before a transition from the downstroke to the upstroke occurred, the root vortex (DRV) grew in a spiral pattern from the leading edge to the trailing edge. The root vortex (DRV) remained attached to the dorsal side at the beginning of the stroke and gradually degenerated until it was replaced by the wing root vortex (URV).
During the upstroke, especially at 0.65 T, the leading edge vortex, the root vortex, and the trailing edge vortex joined together to form a vortex ring on the dorsal side of the wing. At the beginning of the downstroke, no vortex rings were observed on the ventral side of the wing. That was because the trailing edge vortex had dissipated in the wake. There were also no vortex rings on the dorsal side of the wing during the upstroke. The differences in the vortex structures during the two stroke phases were that the leading edge vortex and the tip vortex were more unstable, dispersed, and smaller in size, and showed the characteristics of continuous shedding. Near the end of the upstroke, there were no significant differences in the vortex structures with the previous ones.
The trailing vortex could not be detected according to the current visualization results. The value of 1.0 was changed to a value of 0.01 to visualize the vortex structures around the wing to a greater extent. Figure 12 shows that the periodically shed vortex rings were the most distinguishing features of the wake. At first, the trailing edge vortex was stably attached on the wing surface, and a vortex ring was clearly visible. A vortex ring, consisting of the trailing edge vortex, the root vortex, the tip vortex, and the leading edge vortex, broke off from the wing surface after the stroke reversed. Due to the viscous dissipation of the surrounding flow, the vortex ring became weaker and weaker as it went farther below the hummingbird wing. A detached trailing edge vortex replaced the attached trailing edge vortex in the ring during the upstroke. There existed a smaller vortex ring at the wingtip, and it consisted of the detached leading edge vortex and tip vortex.
3.3. Pressure Distribution and Wake Topology in the Streamwise
To elucidate the effects of the development and evolution of the vortex structures on wing aerodynamics, a slice was extracted at a position near half the wingspan to plot the pressure distribution and streamlines in the flow. And the spanwise flow velocity was thus assumed to be 0 m·s-1 in this plane. The pressure coefficient (the here is lowercase, while the in the identifier of the power coefficient is uppercase) was obtained by normalizing the pressure difference according to where is the static pressure, is reference pressure (1.01325 Pa), and is the pressure difference between the static pressure and reference pressure. The pressure distribution and streamlines in this slice at various time instants are shown in Figure 13 (during the downstroke) and Figure 14 (during the upstroke). The three-dimensional wing model was also plotted in low transparency and green in each snapshot for reference.
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
At the beginning of the downstroke, the leading edge vortex could be seen in this slice. As the downstroke proceeded, the images in Figure 13 show the evolution of the leading edge vortex relatively closer to the outer segment of the wing. The leading edge vortex moved backward along the chord and finally separated from the trailing edge, which was consistent with that discussed in the last section. The leading edge vortex near the wingtip was detached from the wing surface. In contrast to the earlier sustained growth, the vertical force coefficient experienced a slight decrease for the first time. The leading edge vortex in the proximal section became larger and stronger with increasing stroke speed. And it adhered stably to the leading edge of the wing. Therefore, the vertical force coefficient was kept increasing during the subsequent downstroke. The vortex core induced the airflow to move towards the wing surface, and a low-pressure area was thus formed on the upper surface of the wing. The lower surface had a high pressure region, which caused the second peak in the vertical force coefficient. The horizontal force coefficient was also higher at this time as a result of the enhancement in aerodynamic forces. When the wing decelerated before the stroke direction changed, the high-pressure region on the lower surface was gradually shrinking. The leading edge vortex was becoming smaller and more unstable until it detached from the surface, resulting in the vertical force coefficient dropping to the valley from the peak.
In contrast, the lead4ing edge vortex was smaller in size during the upstroke. There was continuous shedding of the leading edge vortex, so it is no wonder that the minority of weight support stems from the vertical force produced during the upstroke. Before the stroke direction changed, the leading edge vortex generated in the earlier phase adhered to the leading edge. Secondly, the leading edge vortex got smaller due to the flipping of the wing. Finally, the new leading edge vortex developed on the upper surface and gradually moved toward the trailing edge. Another new leading edge vortex entered the eye during the middle phase of the upstroke. Its size increased slightly during the subsequent upstroke, and the vortex structures on the upper surface did not show dramatic changes. However, the high-pressure zone on the lower surface was initially enlarged and then diminished. The leading edge vortex and the trailing edge vortex induced strong flow on the wing surface, and the vortex structures attached to the wing surface were also affected by the detached vortex in the wake. Although the spanwise flow was not present in the current analysis, the flow and the pressure distribution around the wing were also the results of the combined action of the tip vortex and the root vortex.
3.4. A Comparison Study on Aerodynamic Characteristics
The time courses of the force and power coefficients obtained by us were compared with those of two previous studies: one on a ruby-throated hummingbird [12] and the other on a rufous hummingbird (Selasphorus rufus) [11]. We found that the vertical force coefficients had similar trends. There were two peaks during the downstroke, whereas the force coefficients showed a fluctuation of multiple peaks as the upstroke continued. Using the morphological and kinematic data provided by Reference [12], we converted those results into equivalent values according to the definitions in this document for comparison. Table 3 lists the results from the comparison study. The ruby-throated hummingbird wing has a shorter wing length, a shorter mean chord, and thus a smaller wing area, but a larger stroke amplitude at the wingtip and a higher stroke frequency bring the average wingtip velocity up to . That value is close to the value of obtained from the Amazilia hummingbird wing.
A time-averaged vertical force coefficient of 0.78 and a time-averaged power coefficient of 0.64 were obtained from the conversion of these results. And their values are very close to the time-averaged vertical force coefficient of 0.81 and the time-averaged power coefficient of 0.69 in the current model (Table 3). But the previous study demonstrated an even higher aerodynamic efficiency of 1.77. And their results also showed that the ratio of the time-averaged vertical force coefficients produced during the two half stroke phases was only 2.49 and smaller than that found in our study. Hummingbirds might maintain a similar average wingtip velocity of 7.50 m/s to hover efficiently. But different species of hummingbirds might adopt different strategies to achieve efficient hovering due to differences in size and weight.
3.5. Aerodynamic Characteristics for Efficient Hovering
The leading edge vortex developed during the two half strokes due to large magnitudes of the kinematic angle of attack across the wing sections. The spanwise twist for each wing section resulted in smaller angles of incidence in the distal sections. So the enhancement in aerodynamics was due to the delayed stall mechanism of the leading edge vortex. The twist was more significant during the upstroke, so the leading edge vortex was smaller and more unstable. The factor was one reason for the presence of the asymmetries in aerodynamic forces and power consumption. In addition, the distance between the leading edge vortex and the trailing edge vortex was constantly increased as the downstroke progressed. That led to a generation of the high first moment of vorticity proportional to the aerodynamic forces. During the upstroke, the trailing edge vortex was much smaller than the vortex core formed at the end of the downstroke. So the first moment of vorticity was lower, explaining why the aerodynamic forces had larger values during the downstroke. The vertical force peaks were absent during the stroke transition because the timing of the wing flipping was earlier than that of the change in the stroke direction. The advanced rotation also avoided the generation of large negative lift generation. Furthermore, the wing flip occurred at the wrist joint while the root was still at a positive angle of attack. A large spanwise twist gradient existed at the inboard section, and thus, the positive and negative components at the inboard and outboard section might cancel each other out, so the vertical force did not fluctuate remarkably. And our results have shown that hovering hummingbirds with an inclined stroke plane produce the most flight forces during the downstroke phase. The wake shed by the hovering hummingbird is best described as a series of stacked vortex rings.
The wing did not interact with the vortex ring during the upstroke (see Figure 13), and it had the spanwise out-of-stroke bending that allowed it to bypass the vortex ring separated at the end of the downstroke. When the wing had finished the upstroke and was ready for the subsequent downstroke, the wake had no significant influence because there was only a smaller and more unstable vortex ring at the wingtip. Therefore, it is reasonable to conclude that such an unsteady high lift mechanism as wake capture did not play a dominant role in hovering hummingbirds.
4. Conclusions
A three-dimensional numerical study was performed for a hovering hummingbird wing with remarkable dynamic morphing. We developed a program based on the User-Defined Function in the commercial software ANSYS Fluent. The method can be used to model a flexible wall that involves fluid-structure interactions. The boundary-based smoothing method and the overset methodology were combined to make the interior nodes in the computational domain accommodate the wing motion in the present numerical simulations. This method can be generalized to other fluid simulations related to biology. The simulation results have shown the characteristics of the aerodynamic forces and the power during a stroke cycle. The vertical force was highly asymmetric, and the aerodynamic power had a similar change. A significant enhancement in the vertical force and horizontal force was found during the downstroke, along with a large number of power consumptions. The thrust was generated during the upstroke to balance the drag generated during the downstroke, causing a decrease in vertical force. It was found that the asymmetry in the time-averaged vertical force during the two half strokes was 3.5. Details of the vortex structures and wake were revealed and discussed in our study. The leading edge vortex (LEV) attached to the wing was stable during the downstroke but extremely unstable and was continuously shed during the upstroke. Successive downstroke wake vortex was connected with distinct vortices detached from the surface to form vortex rings. During the upstroke, the leading edge vortex induced a vortex ring near the root and a smaller vortex ring in the outer segment of the wing. After the vortex ring separated from the wing at the end of the downstroke, it moved to the lower right of the wing. The interaction between these vortexes affects the flow around the wing and the pressure distribution on the wing surface during the entire stroke period. The aerodynamic characteristics for efficient hovering were discussed, and the reasons behind them were explained using wing kinematics. Finally, a comparing study has shown that hummingbirds of different species may maintain a similar average wingtip velocity to achieve efficient hovering.
Data Availability
All data included in this study are available upon request by contact with the corresponding author.
Conflicts of Interest
The authors declare that they have no conflict of interest.
Acknowledgments
The work was supported in part by the National Natural Science Foundation of China under Grant 11872314 and the Key R&D Program in Shanxi Province of China under Grants 2020GY-154, 2020ZDLGY06-05, and 2021ZDLGY09-10.