# Numerical Study of the Aerodynamics of a Hovering Hummingbird’s Wing with Dynamic Morphing

Authors: Yuanbo Dong; Bifeng Song; Dong Xue; Yuanbo Li; Wenqing Yang
Journal: International Journal of Aerospace Engineering
 (2022)

---

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


---

## Body

## 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 201×401 points for each became available.Figure 1
(a) Amazilia hummingbird sucking nectar provided by a feeder, (b) the geometry surface of the wing at the moment, and (c) a corresponding discrete grid of the captured wing.
(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 6×11 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.Figure 2
The comparison between the experimental model and the reconstructed model at two random time instants: (a) group 1 and (b) group 2. The gray in each group represents the experimental model, while the green represents the reconstructed model.
(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.Table 1
Basic morphological and kinematic parameters of the hummingbird wing.
ParameterNomenclatureValueMaximumMinimumFrequencyf28.8Hz——AmplitudeΦ103°——Stroke plane angleβ11.7°——Wing lengthR69.3mm71.2mm66.8mmWing surface areaS1365mm21478mm21250mm2Mean chord lengthc¯m19.5mm——Figure 3
Kinematics of the right wing of an Amazilia hummingbird in a stroke period. (a) The image on the left shows the wing kinematics during the downstroke, while (b) the image on the right represents that during the upstroke.
(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 102−104. 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.
(1)∇⋅V=0,(2)∂ρairV∂t+∇⋅ρairVV=−∇p+∇⋅τ−ρairg,where V is air velocity, p is the pressure, τ is the stress vector, and g 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 250mm and 49mm, 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 (c¯) of 77 mm, an area of the plane of 0.0193 m2, and a second moment of area radius (R2) 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:
(3)ϕ=60°cos2πt∗,(4)ψ=90°−αmsinπt∗Δτ,0≤t∗<Δτ2,αm,Δτ2≤t∗<1−Δτ2,90°−αmcosπt∗−1‐Δτ/2Δτ,1−Δτ2≤t∗<1+Δτ2,90°+αm,1+Δτ2≤t∗<1−Δτ2,90°+αmcosπt∗−1−Δτ/2Δτ,1−Δτ2≤t∗<1,where ϕ and ψ represent the instantaneous harmonic flapping angle and trapezoidal rotation angle, respectively. The normalized time t∗ is defined by t∗=t/T, in which t is the time instant and T is the flapping period. αm 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 45° and 0.3, respectively. We adopted the definition of Re=lUref/ν, in which l is the distance travelled by a point at the R2 position of the wing per flapping cycle and ν is the kinematic viscosity of the fluid medium. The numerical study was conducted at Re=4000, 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.Figure 4
Kinematics of (a) the harmonic model and (b) the trapezoidal model used in the flapping and rotating motions.
(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 D=40c¯, 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(5)∇⋅γ∇u⟶=0,where u⟶ is the mesh displace velocity. The diffusion coefficient γ as a function of the boundary distance is of the form
(6)γ=1dα,where d 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 u⟶ 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 (Δt). After this boundary condition was determined, Equation (5) was solved through the biconjugate stabilized (BCGSTAB) method with finite element discretization and u⟶ was obtained directly at each mesh node. The node positions were then updated according to
(7)x⟶new=x⟶old+u⇀Δt,where x⇀new and x⇀old 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.Figure 5
Comparisons of (a) the lift coefficient and (b) the aerodynamic power coefficient with the previous computational and experimental results of Lua et al. [16].
(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(8)Re=ρfLrefUrefμf=2ρfc¯mΦfRμf=2ρairc¯mΦfRμair=9367.4,where ρf is the fluid density, μf is the dynamic viscosity of a fluid, Lref is reference length (mean chord length c¯m), Uref is reference velocity (average wingtip velocity Uref=2ΦfR=7.1758m⋅s−1), and Φ is the peak-to-peak amplitude of the stroke angle. Air density ρair is 1.205kg⋅m−3, and its viscosity μair is 1.81×10−5Pa⋅s. 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.Figure 6
The mesh system used in the current study: (a) the component mesh, (b) the surface mesh of the hummingbird wing model, (c) the background mesh, and (d) the total mesh system.
(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(9)CF=F1/2ρairUref2S,where F is the vertical (X-axis) or horizontal force (Z-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 T/200 and T/400 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 T/200 was acceptable.Figure 7
Effects of mesh density and time step size on the computed (a) horizontal force coefficient and (b) vertical force coefficient. Solid, long-dashed, dash-dot, and dash-dot-dot lines correspond to 0.7, 1.3, 2.5, and 1.3 million tetrahedral elements, respectively. A time stepT/200 was selected for solid, long-dashed, and dash-dot lines, while for dash-dot-dot lines, the time step was T/400.
(a)(b)


## 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 201×401 points for each became available.Figure 1
(a) Amazilia hummingbird sucking nectar provided by a feeder, (b) the geometry surface of the wing at the moment, and (c) a corresponding discrete grid of the captured wing.
(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 6×11 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.Figure 2
The comparison between the experimental model and the reconstructed model at two random time instants: (a) group 1 and (b) group 2. The gray in each group represents the experimental model, while the green represents the reconstructed model.
(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.Table 1
Basic morphological and kinematic parameters of the hummingbird wing.
ParameterNomenclatureValueMaximumMinimumFrequencyf28.8Hz——AmplitudeΦ103°——Stroke plane angleβ11.7°——Wing lengthR69.3mm71.2mm66.8mmWing surface areaS1365mm21478mm21250mm2Mean chord lengthc¯m19.5mm——Figure 3
Kinematics of the right wing of an Amazilia hummingbird in a stroke period. (a) The image on the left shows the wing kinematics during the downstroke, while (b) the image on the right represents that during the upstroke.
(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 102−104. 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.
(1)∇⋅V=0,(2)∂ρairV∂t+∇⋅ρairVV=−∇p+∇⋅τ−ρairg,where V is air velocity, p is the pressure, τ is the stress vector, and g 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 250mm and 49mm, 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 (c¯) of 77 mm, an area of the plane of 0.0193 m2, and a second moment of area radius (R2) 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:
(3)ϕ=60°cos2πt∗,(4)ψ=90°−αmsinπt∗Δτ,0≤t∗<Δτ2,αm,Δτ2≤t∗<1−Δτ2,90°−αmcosπt∗−1‐Δτ/2Δτ,1−Δτ2≤t∗<1+Δτ2,90°+αm,1+Δτ2≤t∗<1−Δτ2,90°+αmcosπt∗−1−Δτ/2Δτ,1−Δτ2≤t∗<1,where ϕ and ψ represent the instantaneous harmonic flapping angle and trapezoidal rotation angle, respectively. The normalized time t∗ is defined by t∗=t/T, in which t is the time instant and T is the flapping period. αm 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 45° and 0.3, respectively. We adopted the definition of Re=lUref/ν, in which l is the distance travelled by a point at the R2 position of the wing per flapping cycle and ν is the kinematic viscosity of the fluid medium. The numerical study was conducted at Re=4000, 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.Figure 4
Kinematics of (a) the harmonic model and (b) the trapezoidal model used in the flapping and rotating motions.
(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 D=40c¯, 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(5)∇⋅γ∇u⟶=0,where u⟶ is the mesh displace velocity. The diffusion coefficient γ as a function of the boundary distance is of the form
(6)γ=1dα,where d 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 u⟶ 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 (Δt). After this boundary condition was determined, Equation (5) was solved through the biconjugate stabilized (BCGSTAB) method with finite element discretization and u⟶ was obtained directly at each mesh node. The node positions were then updated according to
(7)x⟶new=x⟶old+u⇀Δt,where x⇀new and x⇀old 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.Figure 5
Comparisons of (a) the lift coefficient and (b) the aerodynamic power coefficient with the previous computational and experimental results of Lua et al. [16].
(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(8)Re=ρfLrefUrefμf=2ρfc¯mΦfRμf=2ρairc¯mΦfRμair=9367.4,where ρf is the fluid density, μf is the dynamic viscosity of a fluid, Lref is reference length (mean chord length c¯m), Uref is reference velocity (average wingtip velocity Uref=2ΦfR=7.1758m⋅s−1), and Φ is the peak-to-peak amplitude of the stroke angle. Air density ρair is 1.205kg⋅m−3, and its viscosity μair is 1.81×10−5Pa⋅s. 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.Figure 6
The mesh system used in the current study: (a) the component mesh, (b) the surface mesh of the hummingbird wing model, (c) the background mesh, and (d) the total mesh system.
(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(9)CF=F1/2ρairUref2S,where F is the vertical (X-axis) or horizontal force (Z-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 T/200 and T/400 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 T/200 was acceptable.Figure 7
Effects of mesh density and time step size on the computed (a) horizontal force coefficient and (b) vertical force coefficient. Solid, long-dashed, dash-dot, and dash-dot-dot lines correspond to 0.7, 1.3, 2.5, and 1.3 million tetrahedral elements, respectively. A time stepT/200 was selected for solid, long-dashed, and dash-dot lines, while for dash-dot-dot lines, the time step was T/400.
(a)(b)


## 3. Results and Discussion

### 3.1. Force, Power, and Efficiency

As shown in Figure6, X, Y, and Z 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, FX and FZ, are normalized by air density, the wing area, and the average wingtip velocity according to equation (9). Therefore, CX and CZ represent the horizontal force coefficient and the vertical force coefficient, respectively. The aerodynamic power CP 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
(10)P=−∬pdA+τdA⋅V=−∬pn+τ⋅VdA,where p is the static pressure, τ is the wall shear stress vector, V is velocity vector, A is a surface vector, n is the unit normal vector of the surface, and A is the magnitude of the surface vector A. The power coefficient CP is defined by normalizing the power P by 1/2ρairUref3S.Note that there exist slight cycle-to-cycle variations as you can see in Figure8. 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:
(11)CX¯=1Δt∫tt+ΔtCXdt,(12)CZ¯=1Δt∫tt+ΔtCZdt,(13)CP¯=1Δt∫tt+ΔtCPdt,where t is the start time instant and Δt 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 ηL=CZ¯/CP¯.Figure 8
The two normalized force components (a)CX, (b) CZ in the single-wing simulation, and (c) the aerodynamic power coefficient CP normalized by 1/2ρairUref3S.
(a)(b)(c)Figure 9
The histories ofCX, CZ, and CP in one stroke cycle (obtained by averaging the curves over the five computational periods in Figure 8).
(a)(b)(c)Table 2
Time-averaged values of the horizontal force, vertical force, and power coefficients during the different stroke phases and the corresponding aerodynamic efficiency.
CX¯CZ¯CP¯ηLA stroke cycle0.160.810.691.17Downstroke0.661.261.181.07Upstroke0.340.360.201.80The ratio1.943.505.900.60Figure9shows 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 wasCZ¯=0.81, corresponding to 13.7g of the total weight support provided by the two wings. The time-averaged horizontal force coefficient was CX¯=0.16 and much smaller than the value of CZ¯. A positive CX¯ 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 nondimensionalQ criterion [19]. The positive value is indicative of vortexes in the Q-criterion formulation. Here, the Q value was normalized according to
(14)Q∗=QUref/c¯m2.All three-dimensional vortex structures could be clearly shown by adopting an appropriate value ofQ∗. 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 Q∗=1.0.Figure 10
Three-dimensional vortex structures on the dorsal side of the wing during a stroke cycle. ULEV, URV, UTEV, and UTV are abbreviations for the upstroke leading edge vortex, the upstroke root vortex, the upstroke trailing edge vortex, and the upstroke tip vortex, respectively. Similar labels (DLEV, DRV, DTEV, and DTV) were used to describe the vortex produced during the downstroke.Figure 11
Three-dimensional vortex structures on the ventral side of the wing during a stroke cycle.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. TheQ∗ 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.Figure 12
The evolution of vortex rings during a stroke cycle.


### 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 Cp (the p here is lowercase, while the P in the identifier of the power coefficient is uppercase) was obtained by normalizing the pressure difference according to
(15)Cp=Δp1/2ρairUref2=p−pref1/2ρairUref2,where p is the static pressure, pref is reference pressure (1.01325 Pa), and Δp 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.Figure 13
The contours of pressure coefficient distribution and the streamlines on a plane near the middle of the initial wingspan during the downstroke, proceeding left to right and top to bottom (a)–(f).
(a)(b)(c)(d)(e)(f)Figure 14
The contours of pressure coefficient distribution and the streamlines on a plane near the middle of the initial wingspan during the upstroke, proceeding left to right and top to bottom (a)–(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 Figure13 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 7.49m⋅s‐1. That value is close to the value of 7.58m⋅s‐1 obtained from the Amazilia hummingbird wing.Table 3
A comparison of the wing parameters and simulation results.
An Amazilia hummingbirdA ruby-throated hummingbirdDifferenceWing length69.30mm44.70mm35.5%Mean chord19.50mm11.70mm40.0%Wing area1365mm2586mm258.4%Stroke amplitude103°114°-11.0%Frequency28.8Hz42.0Hz-45.8%Average wingtip velocity7.58m⋅s‐17.49m⋅s‐11.19%Vertical force coefficient0.810.783.70%Power coefficient0.690.647.20%Aerodynamic efficiency1.171.77-51.28%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 (Table3). 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 Figure13), 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.


## 3.1. Force, Power, and Efficiency

As shown in Figure6, X, Y, and Z 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, FX and FZ, are normalized by air density, the wing area, and the average wingtip velocity according to equation (9). Therefore, CX and CZ represent the horizontal force coefficient and the vertical force coefficient, respectively. The aerodynamic power CP 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
(10)P=−∬pdA+τdA⋅V=−∬pn+τ⋅VdA,where p is the static pressure, τ is the wall shear stress vector, V is velocity vector, A is a surface vector, n is the unit normal vector of the surface, and A is the magnitude of the surface vector A. The power coefficient CP is defined by normalizing the power P by 1/2ρairUref3S.Note that there exist slight cycle-to-cycle variations as you can see in Figure8. 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:
(11)CX¯=1Δt∫tt+ΔtCXdt,(12)CZ¯=1Δt∫tt+ΔtCZdt,(13)CP¯=1Δt∫tt+ΔtCPdt,where t is the start time instant and Δt 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 ηL=CZ¯/CP¯.Figure 8
The two normalized force components (a)CX, (b) CZ in the single-wing simulation, and (c) the aerodynamic power coefficient CP normalized by 1/2ρairUref3S.
(a)(b)(c)Figure 9
The histories ofCX, CZ, and CP in one stroke cycle (obtained by averaging the curves over the five computational periods in Figure 8).
(a)(b)(c)Table 2
Time-averaged values of the horizontal force, vertical force, and power coefficients during the different stroke phases and the corresponding aerodynamic efficiency.
CX¯CZ¯CP¯ηLA stroke cycle0.160.810.691.17Downstroke0.661.261.181.07Upstroke0.340.360.201.80The ratio1.943.505.900.60Figure9shows 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 wasCZ¯=0.81, corresponding to 13.7g of the total weight support provided by the two wings. The time-averaged horizontal force coefficient was CX¯=0.16 and much smaller than the value of CZ¯. A positive CX¯ 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 nondimensionalQ criterion [19]. The positive value is indicative of vortexes in the Q-criterion formulation. Here, the Q value was normalized according to
(14)Q∗=QUref/c¯m2.All three-dimensional vortex structures could be clearly shown by adopting an appropriate value ofQ∗. 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 Q∗=1.0.Figure 10
Three-dimensional vortex structures on the dorsal side of the wing during a stroke cycle. ULEV, URV, UTEV, and UTV are abbreviations for the upstroke leading edge vortex, the upstroke root vortex, the upstroke trailing edge vortex, and the upstroke tip vortex, respectively. Similar labels (DLEV, DRV, DTEV, and DTV) were used to describe the vortex produced during the downstroke.Figure 11
Three-dimensional vortex structures on the ventral side of the wing during a stroke cycle.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. TheQ∗ 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.Figure 12
The evolution of vortex rings during a stroke cycle.


## 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 Cp (the p here is lowercase, while the P in the identifier of the power coefficient is uppercase) was obtained by normalizing the pressure difference according to
(15)Cp=Δp1/2ρairUref2=p−pref1/2ρairUref2,where p is the static pressure, pref is reference pressure (1.01325 Pa), and Δp 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.Figure 13
The contours of pressure coefficient distribution and the streamlines on a plane near the middle of the initial wingspan during the downstroke, proceeding left to right and top to bottom (a)–(f).
(a)(b)(c)(d)(e)(f)Figure 14
The contours of pressure coefficient distribution and the streamlines on a plane near the middle of the initial wingspan during the upstroke, proceeding left to right and top to bottom (a)–(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 Figure13 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 7.49m⋅s‐1. That value is close to the value of 7.58m⋅s‐1 obtained from the Amazilia hummingbird wing.Table 3
A comparison of the wing parameters and simulation results.
An Amazilia hummingbirdA ruby-throated hummingbirdDifferenceWing length69.30mm44.70mm35.5%Mean chord19.50mm11.70mm40.0%Wing area1365mm2586mm258.4%Stroke amplitude103°114°-11.0%Frequency28.8Hz42.0Hz-45.8%Average wingtip velocity7.58m⋅s‐17.49m⋅s‐11.19%Vertical force coefficient0.810.783.70%Power coefficient0.690.647.20%Aerodynamic efficiency1.171.77-51.28%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 (Table3). 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 Figure13), 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.



---

*Source: 5433184-2022-02-21.xml*