ABSTRACT

Two-dimensional computational fluid dynamics (CFD) is applied to better understand the effects of wing cross-sectional morphology on flow field and force production. This study investigates the influence of wing cross-section on insect scale flapping flight performance, for the first time, using a morphologically representative model of a bee (Bombus pensylvanicus) wing. The bee wing cross-section was determined using a micro-computed tomography scanner. The results of the bee wing are compared with flat and elliptical cross-sections, representative of those used in modern literature, to determine the impact of profile variation on aerodynamic performance. The flow field surrounding each cross-section and the resulting forces are resolved using CFD for a flight speed range of 1 to 5 m/s. A significant variation in vortex formation is found when comparing the ellipse and flat plate with the true bee wing. During the upstroke, the bee and approximate wing cross-sections have a much shorter wake structure than the flat plate or ellipse. During the downstroke, the flat plate and elliptical cross-sections generate a single leading edge vortex, while the approximate and bee wings generate numerous, smaller structures that are shed throughout the stroke. Comparing the instantaneous aerodynamic forces on the wing, the ellipse and flat plate sections deviate progressively with velocity from the true bee wing. Based on the present findings, a simplified cross-section of an insect wing can misrepresent the flow field and force production. We present the first aerodynamic study using a true insect wing cross-section and show that the wing corrugation increases the leading edge vortex formation frequency for a given set of kinematics.

INTRODUCTION

Biological flapping flight is an extraordinarily complex phenomenon, which occurs across multiple phylums of the animal kingdom, with wide disparities in the balance between aerodynamic performance and efficiency (Ellington, 1985; Walker et al., 2014; Mao and Gang, 2003; Miyake, 2008). The powered flight of insects is inherently unique for its combination of mechanical simplicity and high flapping frequency, resulting in a complicated vortical flow field. These flight characteristics result in agile and maneuverable flight capabilities superior to fixed wing flight such as aircraft (Taylor et al., 2003). Bees are of particular interest because of the utilization of humuli to attach their front and hind wings together during flight, causing the front and hind wings to move as one body (Basibuyuk and Quicke, 1997). Additionally, a bee can carry an additional 80% of its own body weight for miles (Winston, 1991). Bee aerodynamics is uniquely applicable to the future direction of micro-air vehicle (MAV) research as another possible solution to low velocity, low Reynolds number, large payload flight (Ma et al., 2012; Yan et al., 2001; Bai et al., 2007).

To support the creation of more effective MAV systems, it is critical that the current understanding of insect flapping flight is expanded. The first modern explorations into the fluid dynamics of insect flight were experiments utilizing particle image velocimetry (PIV) and other visualization methods (Ellington, 1984; Casey et al., 1985; Mazaheri and Ebrahimi, 2010; Carr and Ringuette, 2014). An experiment by Morse and Liburdy (2009) at low Reynolds number characterized the vortex field around a static wing at high angles of attack for flow over a flat plate. Insect flight is progressively becoming a focus area for fluid and structural researchers alike as computational capabilities improve. The dual improvement of numeric modeling methodologies and computing power is allowing for progressively more accurate and refined computational models of flapping flight (Hu and Xiao, 2014; Nakata and Liu, 2012; Aono et al., 2008). The advances in computational modeling and resulting aerodynamic understanding of complex flight phenomenon is directly aiding the development of robotic MAV systems (Wood et al., 2013). There have been multiple attempts to create low-order models approximating more complex computational fluid dynamics (CFD) results, which could be implemented into various control schemes (Amiralaei et al., 2012; Singh et al., 2004; Taha et al., 2014; Bayandor et al., 2013; Dadashi et al., 2016).

Numerous experimental and computational studies have investigated insect flapping flight (for a thorough review of these topics see Ansari et al., 2006; Saha and Celata, 2011; Shyy et al., 2010). These investigations focused on one of four insect flight strategies: fruit flies (Drosophila), which utilize a single pair of wings to flap at high frequency; dragonflies (Sympetrum), which move four wings independently at high frequency; hawkmoths (Sphingidae), which utilize four wings flapping in synchrony at low frequency; and bees (Apoidea), which utilize four wings flapping at high frequency (Miyake, 2008; Shyy et al., 2010).

Three-dimensional (3D) CFD investigations into insect flight have been performed by numerous groups. Ramamurti and Sandberg (2002) computationally investigated the unsteady flow over a fruit fly and successfully compared it to the experimental results quantified by Dickinson et al. (1999). Additional 3D CFD investigations of a fruit fly have been performed by Walker (2002) utilizing an unsteady blade element model, and Aono et al. (Aono and Liu, 2008; Aono et al., 2008) utilizing an incompressible, unsteady form of the Navier-Stokes equations. The hovering and forward flight of a wide range of insects, including bees, have been analyzed three-dimensionally for both force production and aerodynamic efficiency by Du and Sun (2008, 2012). Du and Sun also reported that a simplified form of the naturally occurring wing corrugation had minimal impact when compared to a flat plate undergoing the same kinematics (Du and Sun, 2012). In addition to the previously discussed 3D simulations, numerous two-dimensional (2D), unsteady computational analyses have been performed by various groups. Wang (Wang, 2000; Wang et al., 2004) studied hoverfly hovering, Lu and Liao (2006) explored the effects of vertical undulation of a rigid ellipse, Ishihara et al. (2009) evaluated the characteristics of passive pitching of a hoverfly undergoing hovering and Kolomenskiy et al. (2011) computationally modeled the differences between 2D and 3D for a generalized flapping system. In addition to their respective insights, each group found good agreement with 3D computational and experimental models of various flapping flight phenomenon.

At present, of the papers directly dealing with computational or experimental studies on insect flight mechanics published since 1998, five papers implement a triangular corrugated cross-section (Meng and Sun, 2011; Du and Sun, 2012; Khurana and Chahl, 2014; New et al., 2014), six papers have used a NACA0012 airfoil (Amiralaei et al., 2009; Chandar and Damodaran, 2008, 2010; Amiralaei et al., 2010; Tian et al., 2014), 10 use an elliptic cross-section (Wang, 2000; Mao and Gang, 2003; Lu and Liao, 2006; Berman and Wang, 2007; Kim and Choi, 2007; Hamdani and Naqvi, 2011; Amiralaei et al., 2011; Zhu and Zhou, 2014; Hu and Xiao, 2014; Du and Sun, 2015), and 23 utilize a flat plate (Liu et al., 1998; Ramamurti and Sandberg, 2002; Sun and Lan, 2004; Wang et al., 2004; Wu and Shao, 2004; Saputra et al., 2006; Huang et al., 2007; Zuo et al., 2007; Aono et al., 2008; Aono and Liu, 2008; Liu, 2009; Ishihara et al., 2009; Senda et al., 2010; Yin and Luo, 2010; Sudhakar and Vengadesan, 2010; Kolomenskiy et al., 2011; Lee et al., 2008, 2011; Nakata et al., 2011; Nakata and Liu, 2012; Ishihara et al., 2014; Shen and Sun, 2015; Erzincanli and Sahin, 2015). Of these papers, seven (Zuo et al., 2007; Meng and Sun, 2011; Luo and Sun, 2005; Mao and Gang, 2003; Berman and Wang, 2007; Liu, 2009; Wu and Sun, 2005) investigate the aerodynamics associated with bee aerodynamics in either hovering or forward flight. Although the aerodynamic effects of a basic corrugation have been determined to be negligible (Meng and Sun, 2011; Meng et al., 2011), the effects of more extreme corrugation associated with biologically accurate wing cross-sections remain unexplored.

To this end, the aerodynamic effects of cross-sectional geometry for a bee over the natural velocity range of bee flight is studied numerically using the unsteady, incompressible form of the Navier-Stokes equations in two dimensions. The goals of the present study are twofold. The first is to determine if misrepresentation of wing cross-section in insect flight adversely influences observed force and flow phenomenon. The second is to determine how the results of commonly used cross-sections compare with those using morphologically accurate venation. Transient and average coefficients of lift and drag are compared in depth to improve the understanding of the possible aerodynamic differences due to wing cross-section. The wing kinematics are based on high speed videography by Ellington and Dudley (1990) and are used as a basis for equations of motion to approximate the insect wing path. This work is the first computational, morphologically accurate analysis of a bee wing cross-section with direct comparisons to common cross-sectional geometries presented in the literature. The underlying fluid mechanics caused by the wing cross-sectional geometry can manifest in extreme variations in thrust, lift, and vortex development for different cross-sections undergoing the same kinematics. The results presented will determine whether the complexity of a bee wing cross-section can be simplified for further modeling and will establish the necessity of morphologically accurate models of wing geometry to predict representative aerodynamic forces for bee flight, and further, may yield additional passive methods to control vortex formation during flight that have not been studied previously.

RESULTS

Flow characteristics

The flow fields at τ=11.25 (τ=0.25 in Fig. S2) and 11.68 (between τ=0.5 and τ=0.75 in Fig. S2) are shown in Figs 1 and 2, respectively, to gain insight into the underlying aerodynamics that manifests as differences in lift and drag. The parameter tau is normalized by the flap period. Parts A and B of both figures show the pressure and vorticity distributions with streamlines, respectively, in a portion of the flow field directly around the wing. The wing cross-sections are organized in descending order from most (bee cross-section) to least geometrically complex (flat plate), as shown in Fig. 3.

Fig. 1.

Pressure and vorticity contours at τ=11.25 for the four wing cross-sections at 5 m/s. Pressure (left column) and streamlines with vorticity (right column) contours at τ=11.25 are shown in the (A) bee, (B) approximated, (C) ellipse and (D) flat plate cross-sections.

Fig. 1.

Pressure and vorticity contours at τ=11.25 for the four wing cross-sections at 5 m/s. Pressure (left column) and streamlines with vorticity (right column) contours at τ=11.25 are shown in the (A) bee, (B) approximated, (C) ellipse and (D) flat plate cross-sections.

Fig. 2.

Pressure and vorticity contours at τ=11.68 for the four wing cross-sections at 5 m/s. Pressure (left column) and streamlines with vorticity (right column) contours at τ=11.68 are shown in the (A) bee, (B) approximated, (C) ellipse and (D) flat plate cross-sections.

Fig. 2.

Pressure and vorticity contours at τ=11.68 for the four wing cross-sections at 5 m/s. Pressure (left column) and streamlines with vorticity (right column) contours at τ=11.68 are shown in the (A) bee, (B) approximated, (C) ellipse and (D) flat plate cross-sections.

Fig. 3.

The cross-sections used in the present analysis. (A) Bee, (B) approximated, (C) ellipse and (D) flat plate.

Fig. 3.

The cross-sections used in the present analysis. (A) Bee, (B) approximated, (C) ellipse and (D) flat plate.

Beginning with Fig. 1 (τ=11.25), a leading edge vortex (LEV) is already in the process of shedding from the hind wing section (lowest pressures), while the front bee wing section shows the initial formation of an LEV for both the bee and approximated wing cross-sections. The trailing edge vortex (TEV) is almost fully detached from the bee cross-section, visible in the location of the low pressure region on the upper portion of the wing, with distinct lack of vorticity in the same region. The LEV on the approximated wing section is still in the process of detaching from the wing. The bee cross-section sheds the hind wing LEV earlier in the stroke than the approximate section. The postponement of vortical shedding increases the area and pressure gradient of the vortex for the approximate cross-section, increasing lift generation. Both the ellipse and flat plate cross-sections are still in the process of forming an LEV on the front of the wing. The underside of the ellipse wing exhibits a small low pressure region near the trailing edge, transitioning to a high pressure region towards the leading edge. The transition from high to low pressure is similar to what is seen in the bee and approximated sections, where the pressure transitions from high to low to high, with the high pressure regions caused by stagnation points on the cross-section. Conversely, the flat plate section has low pressure on the upper and high pressure on the lower surfaces of the wing. These variations in pressure are caused by differences in the vortex generation and shedding capabilities of the cross-sectional geometry of the wing.

Examining Fig. 2, the flow field surrounding the four cross-sections are compared at τ=11.68. The bee and approximate cross-sections both exhibit multiple vortices shed from both the leading and trailing edges. The bee and approximate wing cross-sections have two attached vortex structures (LEV and TEV) on the upper portion of the wing. The bee cross-section exhibits a second set of vortices located further downstream which are not present in the approximate section. The ellipse and flat plate sections show progressively differing flow structures. The ellipse cross-section exhibits attached TEV and LEV similar to those seen in the approximated and bee sections, with a higher vorticity magnitude. The ellipse cross-section only shows a single remnant LEV entrained in the flow. Although the number of vortex structures immediately surrounding the ellipse are similar to those seen by the approximate and bee cross-sections, the overall locations of the vortices are different. The difference in vortex locations relative to the wing itself drastically changes the distribution of pressures along the wing, moving the low pressure region associated with the LEV towards the cross-sectional center. Comparing the bee, approximate and ellipse cross-sections with the flat plate, there is a distinct reduction observed in flow complexity with the flat plate, which generates a large single LEV.

DISCUSSION

The goal of the present study is to establish the aerodynamic differences among wing cross-sections by comparing biologically accurate wing sections with cross-sections commonly used in the computational literature. In order to establish the effects of biologically accurate wing cross-section on bee flight aerodynamics, the wing kinematics outlined in section are applied to all four cross-sections shown in Fig. 3. The induced fluid mechanics are compared using vorticity and pressure contours with the resulting forces compared from both an instantaneous and time averaged perspective.

Owing to the identical nature of the kinematics, differences in the flow field and resulting forces will be considered to be caused by the cross-sectional variation. Throughout the discussion, τ will be used to define the phase of the wing flap stroke. The discussion will refer to a given phase after a number of flap cycles as τ=11.25, where the solution is 0.25 τ after 11 strokes. The relative position of τ is shown in Fig. S2.

Variations in lift coefficient

Fig. 4 presents the CL data for four flap cycles comparing the flat plate, ellipse, approximate and bee wing cross-sections at τ at V=1−5 m/s. At all speeds, the bee and approximate wings exhibit two distinct peaks in CL, which correlate with LEV shedding at approximately τ=11.75 (the midpoint of the upstroke) and TEV shedding at τ=11.75 (midpoint of the downstroke). At 1 m/s, all cross-sections produce a very similar transient CL. At 2 m/s the flow field of the flat plate begins to differ from the other three cross-sections such that the LEV structure does not shed at the stroke transition, causing a more gradual drop in lift production. The difference in aerodynamic forces become more notable at 3 m/s where the flat plate and, to a lesser degree, the ellipse do not exhibit the reduction in CL at τ=11.5. As V increases to 5 m/s, the elliptical wing profile predicts similar CL to the bee and approximate wings. Meanwhile the flat plate profile continues to produce very high lift, more than double the peak force experienced by the bee cross-section. The trends at 4 m/s remain consistent with those observed at 3 m/s, with a progressively different shedding pattern as wing cross-sectional geometry differs from the biological model. At 5 m/s, CL for the bee wing and the approximate wing remains in good agreement. The peak CL is substantially higher for the flat plate and the minimum CL is significantly lower for both the ellipse and flat plate cross-sections. The difference in lift generation indicates that a simplified wing cross-sections will not necessarily yield biologically representative CL.

Fig. 4.

The time history of CL for the bee, approximated, elliptic and plate cross sections are shown over four flap periods. Each graph is associated with a different flow velocity beginning at 1 m/s at top left and ending at 5 m/s at the bottom.

Fig. 4.

The time history of CL for the bee, approximated, elliptic and plate cross sections are shown over four flap periods. Each graph is associated with a different flow velocity beginning at 1 m/s at top left and ending at 5 m/s at the bottom.

Comparing time-averaged CL, in Table 1, all cross-sections except for the elliptic one exhibit a significant increase in CL at 2 m/s. After 2 m/s, the CL for the bee, approximate and ellipse cross-sections decrease steadily, with the exception of 5 m/s for the approximate cross-section, which increases to 0.1. The flat plate exhibits a consistent CL from 2 to 5 m/s, while the other cross-sections display a progressive decrease. The approximate cross-section CL is consistently closest to the bee results, with a maximum difference of 0.03.

Table 1.

Average and RMSDCLandCDvalues for each cross-section at velocities from 1 to 5 m/s. RMSD values, given in parentheses below the average values, were calculated relative to the bee cross-sectional results

Average and RMSD CL and CD values for each cross-section at velocities from 1 to 5 m/s. RMSD values, given in parentheses below the average values, were calculated relative to the bee cross-sectional results
Average and RMSD CL and CD values for each cross-section at velocities from 1 to 5 m/s. RMSD values, given in parentheses below the average values, were calculated relative to the bee cross-sectional results

The root mean squared difference (RMSD) will be used to quantify the transient fluctuations from the biologically accurate cross-section present in the instantaneous data. Investigating both time-averaged force results and RMSD allows for the quantification of both instantaneous and average error.

The CL root mean square difference for each of the cross-sections is compared to the bee, where the is calculated as
formula
(1)
and n is the number of time steps analyzed over 1600 time steps from τ=3–15. The is proportional to the instantaneous difference in CL from the bee cross-section; the larger the value, the more disparate the instantaneous forces. Table 1 shows the of each of the cross-sections relative to the bee. At 1 m/s the ellipse RMSDCL was the closest to the bee. Once the velocity increases to 2 m/s, the ellipse and flat plate average CL begin to differ significantly from the bee results while the remains consistent. Increasing velocity from 2 m/s to 3 m/s, the of the approximate cross-section decreases by ∼50% while the flat plate remains consistently high. For velocities between 3 and 5 m/s, the for each of the cross-sections remains the same with the approximate section consistently exhibiting the closest similarity, followed by the ellipse and flat plate cross-sections.

Variations in drag coefficient

Fig. 5 presents the CD data for four flap cycles comparing each wing cross-section versus τ at V=1−5 m/s. The CD values at 1 m/s agree well between all four cross-sections. Each cross-section exhibits an initial peak at the start of the upstroke, before dropping to the minimum CD at τ=11.2. After τ=11.2, CD increases steadily to the maximum around τ=11.7. After the maximum at τ=11.7, there is another drop in CD at τ=11.9. From τ=11.9 to τ=12, CD increases again, peaking immediately after transitioning back to the upstroke.

Fig. 5.

The time history of CD for the bee, approximated, elliptic and plate cross sections are shown over four flap periods. Each graph is associated with a different flow velocity beginning at 1 m/s at top left and ending at 5 m/s at the bottom.

Fig. 5.

The time history of CD for the bee, approximated, elliptic and plate cross sections are shown over four flap periods. Each graph is associated with a different flow velocity beginning at 1 m/s at top left and ending at 5 m/s at the bottom.

At 3 m/s, the CD of the ellipse and flat plate cross-sections begin to shift from the approximate and bee wings. All four cross-sections still follow the same CD trends discussed with 1 m/s, with a number of minor variations. The peak at τ=11 is significantly higher for both the flat plate and ellipse sections, which agree well with one another, unlike the CD for the approximate and true bee wing. The minimum CD at τ=0.2 is significantly lower for the flat plate with a phase offset of 0.1τ, while the minima for the ellipse, approximate and flat plate agree well in both phase and magnitude.

Comparing the time-averaged CD values in Table 1, the differences in average CD are small between the four cross-sections for all of the flight velocities. The airfoils develop thrust at 1 and 2 m/s, but experience a net drag at higher velocities. The approximate section is closest to the bee CD for all simulations except 1 m/s. At 2 m/s the flat plate cross-section transitions from generating thrust to net drag production slightly sooner than the other cross-sections. The CD root mean square difference for each of the cross-sections is compared to the bee, where the is calculated as
formula
(2)
The for each cross-section relative to the bee is shown in Table 1. The ellipse section had the lowest at 1 m/s. Increasing the velocity to 2 m/s, the approximate section increases in accuracy while the elliptical and flat plate decrease, counter to what was seen with . From 3 to 5 m/s, the remains relatively consistent for the approximate wing and decreases consistently for the ellipse and flat plate. At 5 m/s, the approximate cross-section remains the closest to the bee, with the ellipse and flat plate being roughly two and three times the approximate , respectively.

The similarity in the vorticity and pressure contours manifest in the similarities in CL and CD as previously discussed, both average and RMSD. Evidenced by both instantaneous forces and the flow field in Fig. 1, the wing cross-sectional complexity coincides with distinctive changes in vortical shedding, pressure distribution over the wing, and force production at τ=11.25. Conversely, while the four cross-sections exhibit very similar instantaneous CL and CD at τ=11.68, it is clear that there can still be large disparities in the respective flow fields evidenced by the extreme variation in pressure and vorticity distribution at that instant.

The findings presented demonstrate the effects of wing profile on both the immediate flow field and the resulting aerodynamics. The approximated wing profile that represents a streamlined bee cross-section produced the most consistent similarity with the results of the bee cross-section across all metrics. The ellipse profile was the second best at representing the aerodynamics associated with the bee cross-section, capturing the underlying physics more consistently than the flat plate. However at low velocities, the profile thickness of the ellipse changed the underlying vortex dynamics. The flat plate generally captured the aerodynamics associated with the bee-profile the least, with a large LEV consistently forming and not detaching until well into the downstroke. The large LEV formation and shedding experienced by both the flat plate and ellipse, caused by the transition from the up to the downstroke (τ=0.5), does not fully shed from the wing until τ=0.6-0.7. This is earlier than the vortex structure shedding for the bee and approximate cross-sections, which is much closer to τ=0.5. The large LEV structure creates a larger net lift, consistently increasing both the average and . Finally, the present study shows that neglecting the intricacies of a highly complex wing geometry can produce results that, while physical, do not properly capture the surrounding flow field or resulting forces experienced by the biological system. The difference in performance of the bee versus the approximated, ellipse and flat plate wing cross-sections highlights the importance of high resolution modeling of insect wings, both to properly capture the underlying physics using CFD, and to maximize the effectiveness of insect-based MAV applications.

Summary

The influence of wing cross-section on aerodynamic performance was investigated using flight velocities and flapping kinematics for a bee. Four different representations of an insect wing cross-section were compared using the cross-section from a micro-computed tomography scan of an actual bee wing and an approximation of that bee wing cross-section. Two other cross-sections were also compared based on the most prevalent in the literature: an ellipse, the second most prevalent, and a flat plate cross-section, the most prevalent. Each wing cross-section was simulated using the same mesh density, flight speed, and kinematics. Each wing profile was compared using the surrounding flow field, instantaneous CL and CD, time-averaged CL and CD, and RMSD of CL and CD for each cross-section to gain a more quantitative representation of the instantaneous variation in CD and CL.

Qualitatively, these kinematic motions exhibited similar vortex patterns among the four cross-sections, generating LEV-TEV pairs. Vortex formation differed significantly in the number and timing of shedding for each cross-section causing distinctly different resulting forces. The bee and approximate cross-sections were generally comparable both qualitatively and quantitatively to one another. Vortex shedding complexity was shown to increase with wing cross-section complexity, with the flat plate and elliptical cross-sections often forming a single LEV structure over a flap period while the more morphologically accurate model shedded two distinct LEV structures for the same kinematics. With regards to CD, net thrust was produced for all wings at low velocities before transitioning to net drag production at higher velocities. The present study has shown that cross-sectional morphology can significantly influence the occurrence and frequency of LEV formation, changing observable flow structures and aerodynamic performance. This change in vortex formation and shedding frequency augments the effects of dynamic stall, a major source of lift production for biological flight. Based on the present findings, the influence of wing cross-section on aerodynamic performance cannot be assumed to be consistently negligible at the insect scale for future computational analysis and MAV applications, and should be evaluated on a case-by-case basis, with better understanding of cross-sectional effects yielding an additional passive mechanism for LEV formation control.

MATERIALS AND METHODS

Governing equations

The commercial software ANSYS Fluent (v. 15) is used to solve the velocity and pressure fields around a dynamic 2D wing cross-section. The unsteady, incompressible Navier-Stokes equations are employed, yielding the continuity equation
formula
(3)
where t is time, is the velocity vector and ρ is the fluid density. Neglecting gravitational effects, the momentum equations are
formula
(4)
where P is the pressure and is the fluid stress tensor.
The segregated pressure-based Navier-Stokes (PBNS) solver is used to solve the incompressible flow. The SIMPLE algorithm is used to resolve the pressure-velocity coupling by enforcing mass conservation using a relationship between element pressure and flux corrections (ANSYS Incorporated, 2014). The gradients are discretized using the least squares cell based (LSCB) method, momentum equations are discretized using the second-order upwind scheme, and time is discretized using a first-order implicit method. The absolute convergence criteria for the solution variables are set to 10–7. The time step is determined using a modified version of the Courant-Fredrichs-Levy (CFL) number to suit the implementation of kinematic motion:
formula
(5)
where Δt is the time step, ΔX is the smallest cell size, and VxVy are the respective kinematic velocities in the x- and y- directions. The maximum CFL value used in the present simulations is 1 to ensure a stable calculation.

Dynamic remeshing equations

An unstructured finite-volume mesh is generated around the cross-section and is updated with each time step to allow for proper mesh motion, using a combination of Laplacian smoothing and remeshing functions to retain mesh density around the wing at all points in time during the flapping motion. An unstructured, triangular 2D mesh is created around the wing profile (an example of the mesh immediately around the bee wing is shown in Fig. S1). Air enters the computational domain along the right boundary with a uniform velocity specified. The upper and lower boundaries of the fluid domain use a slip condition. The left side of the fluid domain is specified as ambient pressure (0 gauge). It should be noted that the entire domain is 88 mm×44 mm (22c×11c).

The diffusive smoothing function for the dynamic meshing is defined as
formula
(6)
where the diffusion coefficient γ defines the degree to which the boundary motion propagates through the surrounding fluid mesh. γ is calculated using the cell distance from the deforming boundary d:
formula
(7)
Solving Eqn 6 for , the mesh distribution at the next time step is solved using , Δt and the present mesh distribution :
formula
(8)
The diffusion parameter (a) was set to 1.75 for the study, and typically ranges between 0 and 2, where a=1 produces uniform diffusion and values greater than one cause regions further from the boundary to deform, while retaining resolution in the immediate vicinity of the moving boundary.
A remeshing function is implemented to resolve the extreme boundary motions present in flapping flight. A secondary background mesh is used to retain the consistency of the fluid data during remeshing. The remeshing function uses the nodal distance from the nearest boundary dmin and the most remote node from the boundaries dmax to normalize the boundary distance db:
formula
(9)
db is then used to determine the cell size at the location of interest (sizei)
formula
(10)
where Γ is the sizing function factor, defined by the size function variation α and the size function rate β. Two different equations for Γ are implemented depending on the sign of α:
formula
(11)
Combining Eqns 10 and 11 yields the allowable size of a cell in the fluid volume. Using the remeshing functions with diffusion smoothing creates a robust, highly refined mesh.

Grid resolution study

The overall lift and drag coefficients for a wing profile are defined as
formula
(12)
formula
(13)
where Fy=L is the lift force and Fx=TD is the force due to the difference between thrust T and drag D, A is the wing planform area, and Vmag is the maximum translational velocity magnitude of the wing. A positive CL means that the wing is producing lift while a positive CD represents the generation of thrust. Conversely, a negative CD indicates a drag-dominated horizontal force.

The grid convergence index (GCI) methodology is used to determine the numerical accuracy of the solution using the difference between progressively refined meshes (Celik et al., 2008), where the number of cells N is used to quantify change in resolution. The grid resolution study is performed using a static mesh but with similar cell distributions as would be expected in the dynamic mesh simulation. Due to the nature of utilizing a static wing, Vmag is replaced with the inlet velocity V for calculating Δt using CFL. The grid resolution study examined steady-state lift and drag coefficients. The calculated order of accuracy of the simulations was found to be third order with a CD at 506k yielding the maximum numeric error, 0.8%. Due to the variable nature of the mesh density when utilizing a remeshing scheme, the minimum resolution remained higher than 506k cells for all simulations presented herein (the full GCI results are given in Table S1).

2D validation

The CFD methodology is validated using digital particle image velocimetry (DPIV) experimental data reported by Wang et al. (2004) for flow over a hoverfly wing. The construction of the wing was such that the planform was the same shape as the hoverfly but with a uniform thickness. Thus, from a side view, the projected area is a simple rectangular cross-section. The 2D equivalent to the 3D geometry is analyzed as a flat plate with infinite depth in the CFD model. In the experiment, a robotic arm with three rotational degrees of freedom (DOF) and a flat Plexiglas wing with a hoverfly (Drosophila) planform was actuated in still mineral oil for Reynolds numbers between 50 and 200. The kinematics used in the in the validation are based on a fluid density of 880 kg/m3 and a kinematic viscosity of 1.15×10–4 m2/s. The kinematics used in the simulation are f=0.25 Hz, r=0.1625 m, θ1a=23.5°, θ2a=0°, θ3a=45°, φ1=π/2, φ2=0 and c=0.02385 m. All boundaries in the computational domain were modeled as no-slip walls to replicate the experimental tank. A schematic of the plate geometry orientation and kinematics are shown in Fig. 6. Fig. 6 and the remaining discussion utilize nondimensionalized time τ=ft to describe the location of the kinematics. For example, τ=4.25 coincides with the position τ=0.25 in Fig. 6 after four complete flap cycles. At τ=10 the flap cycle repeats, passing through the τ=0.125 location again.

Fig. 6.

Schematic of flat plate moving with validation kinematics over a full flap cycle. The upstroke is shown at the top with an arrow signifying motion to the left. The motion reverses for the downstroke (at the bottom) with the wing translating to the right.

Fig. 6.

Schematic of flat plate moving with validation kinematics over a full flap cycle. The upstroke is shown at the top with an arrow signifying motion to the left. The motion reverses for the downstroke (at the bottom) with the wing translating to the right.

Fig. 7 presents the DPIV data reported by Wang et al. (2004) (left) and the vorticity of the 2D simulations (right). The data at τ=4.2 (top), τ=4.5 (middle) and τ=4.7 (bottom) are compared and show good qualitative agreement between the simulations and the experiment with a very similar vortex shedding pattern. At τ=4.2, the wing is moving towards the left while rotating clockwise, causing a counter-clockwise vortex to form at the (lower) trailing edge and a clockwise vortex to form on the (upper) leading edge, which is captured in both the simulation and experiment. At τ=4.5, the wing starts moving to the right and is nearly vertical. The clockwise LEV that formed at τ=4.2 grows along the right side of the plate while the counter-clockwise TEV almost detaches. At τ=4.7 the wing moves further right, rotating clockwise while doing so. The motion causes a counter-clockwise vortex to form on the left edge of the plate and the LEV from τ=4.5 completely attaches.

Fig. 7.

Comparison of vorticities for the experimental and current simulations. Qualitative comparison of vorticities for the (A) DPIV experimental (Wang et al., 2004) and (B) current simulations, where red signifies counter-clockwise and blue clockwise vortical rotation. The (C) CL and (D) CD versus τ for the present simulation methodology and the computational and experimental data of Wang et al. (2004) are shown.

Fig. 7.

Comparison of vorticities for the experimental and current simulations. Qualitative comparison of vorticities for the (A) DPIV experimental (Wang et al., 2004) and (B) current simulations, where red signifies counter-clockwise and blue clockwise vortical rotation. The (C) CL and (D) CD versus τ for the present simulation methodology and the computational and experimental data of Wang et al. (2004) are shown.

In addition to the DPIV measurements, a 2D force sensor was attached to the base of the wing to capture the transient CL and CD of the wing (Wang et al., 2004). The transient simulation results shown in Fig. 7 exhibit close similarities to both the experimental and computational results reported by Wang et al. (2004). The largest differences occur at the stroke transition (e.g. τ=0.55 and 1.05). The variations at stroke transition are due to 3D effects in the experiment that are not captured in a 2D simulation. Quantitatively, the predicted time-averaged CL and CD values of 0.82 and 1.33 agree well with the experimental values given in Wang et al. (2004) of 0.86 and 1.34, respectively. The validation adds confidence to the 2D computational modeling to further pursue the wing profile analysis.

Cross-sectional geometries

The cross-sections shown in Fig. 3 are those used in the remainder of the analysis. In descending order, the wings cross-sections are from a real bee, an approximation of the bee profile, an ellipse and a flat plate. The profiles in Fig. 3 are shown with the leading edge at the right and trailing edge at the left, with a chord length of c=4 mm.

The bee wing cross-section (Fig. 3A) used in the analysis is from a worker bee of the species Bombus pensylvanicus. The bee had died with wings fully extended, preserving the natural humuli connection between the front and hind wings. A Skyscan 1172 high resolution micro-computed tomography scanner with a resolution of 2 μm was used to generate a 3D model of the wing. The cross-section shown in Fig. 3 is taken at the midpoint of the wing between the wingtip and joint, capturing the fore and hind wing sections in addition to the humuli interface.

The approximate wing profile (3b) is based on the bee wing cross-section using the original point cloud of the bee wing, with manually smoothed vein structures. The removal of the vein structures and smoothing of the overall profile will help establish whether the irregular corrugation of the real wing dramatically affects the vorticity development. The ellipse (3c) and flat plate cross-sections (3d) are geometric representations commonly implemented in the literature instead of the biologically accurate insect cross-section. The elliptical cross-section has a maximum thickness of 0.125c, where thicknesses between 0.1c and 0.3c (Berman and Wang, 2007; Wu and Sun, 2005) were used in the bee-related literature. The flat plate has a uniform thickness of 0.0125c, where in the literature thicknesses from 0 to 0.03c (Zuo et al., 2007; Meng and Sun, 2011) are used.

Wing kinematics

The bee wing path is modeled as three angular velocities around the base of the wing. The rotational axis used in the present study are yaw (θ1), roll (θ2) and wing pitch (θ3). The angular velocities θ1, θ2 and θ3 are used to replicate the wing path during insect flight as a 2D simulation.

Kinematic data by Ellington and Dudley (Ellington, 1984; Ellington and Dudley, 1990) follow the path of a bee wingtip in forward flight at a variety of velocities using a high speed camera (Ellington and Dudley, 1990). The results are approximated as a change in magnitude of the θ2 component based on flight speed where θ2A is equal to 10°, 15°, 20°, 24° and 26° for flight speeds from 1 to 5 m/s, respectively. The frequency f and angular amplitudes θ1A and θ3A are assumed constant at 150 Hz, 110° and 24°, respectively. The Reynolds number Re is used to nondimensionalize the flow, where Re=ρ(Vkin+V)c/μ. Due to the nature of the problem, two different velocities are possible: the flow over the wing driven by the kinematics (Vkin), and the flow over the wing by the freestream velocity (V). For the flight speeds under consideration, the Reynolds number varies from 987-1672, and is considered laminar flow.

The 3D angular amplitudes (θ1A, θ2A, θ3A) must be translated into a continuous two-dimensional motion. Of the three amplitudes, θ3 remains in angular form, as it is the angle of incidence of the wing. The approximated sinusoidal equation for the angular position is
formula
(14)
where t is the current time, ω is the angular frequency, and φ is the phase shift of the wing rotation based on the high speed videography performed by Ellington and Dudley (1990). The remaining two rotations, θ1 and θ2, are applied as translational motions in the x-and y-directions in 2D space. The translational equations are
formula
(15)
formula
(16)
where r° is the distance from the wing-base to the midpoint (r°=6.61 mm). The time derivatives of Eqns 14 and 15 are taken to determine the velocity of the wing relative to the bee in flight.

These 2D kinematic equations describe the representative up-and-down strokes shown in Fig. S2, moving from left to right. Starting from 0° angle of incidence (τ=0), the wing begins its upstroke, decreasing pitch angular velocity while increasing translational velocity until reaching 24° pitch at the upstroke midpoint (τ=0.25). For the latter half of the upstroke, the angular velocity increases while the translational velocity decreases until θ3 returns to 0°. At the top of the upstroke (τ=0.5) the wing translational velocities are zero, while the angular velocity is at the peak. The downstroke reverses the previous kinematics, returning to the original position at τ=1. The process is repeated for each beat of the wing.

Domain dependency study

Solution dependency on computational domain size is explored for insect flight at a flight velocity of 5 m/s for domains of 14c×11c, 22c ×11c and 33c×11c. The bee wing cross-section is used for all analyses. Only the horizontal domain length is varied in the analysis. The vertical dimension remains constant because it was shown that a domain height of 10c was sufficient in previous studies (Kolomenskiy et al., 2011; Ishihara et al., 2014; Du and Sun, 2015). The instantaneous vorticity in the wing wake is shown in Fig. S3 and resulting aerodynamic forces (CL and CD) are shown in Fig. S4. Based on the aerodynamic forces and observed wake structures between the three computational domains, a domain of 22c ×11c is used for the analysis.

Acknowledgements

We thank Professor Rolf Mueller and Professor Jake Socha for the use of the Skyscan 1172; Todd Estep for his assistance operating the scanner; and the members of the CRASH and CREST laboratories for their support, feedback and attention.

Footnotes

Author contributions

Conceptualization: J.F., F.B., J.B.; Methodology: J.F.; Resources: J.F.; Writing - original draft: J.F., F.B., J.B.; Writing - review & editing: J.F., F.B., J.B.; Visualization: J.F.; Supervision: F.B., J.B.; Project administration: F.B., J.B.

Funding

This research received no specific grant from any funding agency in the public, commercial or not-for-profit sectors.

References

Amiralaei
,
M. R.
,
Alighanbari
,
H.
and
Hashemi
,
S. M.
(
2009
).
Influence of unsteady parameters on the aerodynamics of a low reynolds number pitching airfoil
.
ASME 2009 International Mechanical Engineering Congress & Exposition
, Vol.
9
,
Pts A, B and C
,
323
-
329
.
Amiralaei
,
M. R.
,
Alighanbari
,
H.
and
Hashemi
,
S. M.
(
2010
).
An investigation into the effects of unsteady parameters on the aerodynamics of a low Reynolds number pitching airfoil
.
J. Fluid Struct.
26
,
979
-
993
.
Amiralaei
,
M. R.
,
Alighanbari
,
H.
and
Hashemi
,
S. M.
(
2011
).
Flow field characteristics study of a flapping airfoil using computational fluid dynamics
.
J. Fluid Struct.
27
,
1068
-
1085
.
Amiralaei
,
M. R.
,
Partovibakhsh
,
M.
and
Alighanbari
,
H.
(
2012
).
Application of Adaptive Network-Based Fuzzy Inference System (Anfis) in aerodynamics prediction of low-reynolds-number flapping motion
.
ASME 2010 International Mechanical Engineering Congress and Exposition
, Vol.
7
,
Pts A and B
,
307
-
312
.
Ansari
,
S. A.
,
Żbikowski
,
R.
and
Knowles
,
K.
(
2006
).
Aerodynamic modelling of insect-like flapping flight for micro air vehicles
.
Prog. Aerosp. Sci.
42
,
129
-
172
.
ANSYS Incorporated
(
2014
).
Fluent 15.0 User's Guide, Lebanon, NH
.
Aono
,
H.
and
Liu
,
H.
(
2008
).
A numerical study of hovering aerodynamics in flapping insect flight
.
In Bio-mechanisms of Swimming and Flying (ed. N. Kato and S. Kamimura), pp. 179-191. Tokyo: Springer
.
Aono
,
H.
,
Liang
,
F.
and
Liu
,
H.
(
2008
).
Near- and far-field aerodynamics in insect hovering flight: an integrated computational study
.
J. Exp. Biol.
211
,
239
-
257
.
Bai
,
P.
,
Cui
,
E.
,
Li
,
F.
,
Zhou
,
W.
and
Chen
,
B. Y.
(
2007
).
A new bionic MAV's flapping motion based on fruit fly hovering at low Reynolds number
.
Acta Mech. Sin.
23
,
485
-
493
.
Basibuyuk
,
H. H.
and
Quicke
,
D. L. J.
(
1997
).
Hamuli in the Hymenoptera (Insecta) and their phylogenetic implications
.
J. Nat. Hist.
31
,
1563
-
1585
.
Bayandor
,
J.
,
Bledt
,
G.
,
Dadashi
,
S.
,
Kurdila
,
A.
,
Murphy
,
I.
and
Lei
,
Y.
(
2013
).
Adaptive control for bioinspired flapping wing robots
.
American Control Conference (ACC), 2013
.
Berman
,
G. J.
and
Wang
,
Z. J.
(
2007
).
Energy-minimizing kinematics in hovering insect flight
.
J. Fluid Mech.
582
,
153
-
168
.
Carr
,
Z. R.
and
Ringuette
,
M. J.
(
2014
).
Flow structure of low-aspect-ratio rotating wings from dye visualization
.
AIAA J.
52
,
1081
-
1086
.
Casey
,
T. M.
,
May
,
M. L.
and
Morgan
,
K. R.
(
1985
).
Flight energetics of euglossine bees in relation to morphology and wing stroke frequency
.
J. Exp. Biol.
116
,
271
-
289
.
Celik
,
I. B.
,
Ghia
,
U.
,
Roache
,
P. J.
,
Freitas
,
C. J.
,
Coleman
,
H.
and
Raad
,
P. E.
(
2008
).
Procedure for estimation and reporting of uncertainty due to discretization in CFD applications
.
ASME J. Fluids Eng.
130
,
078001
.
Chandar
,
D. D. J.
and
Damodaran
,
M.
(
2008
).
Computational study of unsteady low-reynolds-number airfoil aerodynamics using moving overlapping meshes
.
AIAA J.
46
,
429
-
438
.
Chandar
,
D. D. J.
and
Damodaran
,
M.
(
2010
).
Computation of unsteady low reynolds number free-flight aerodynamics of flapping wings
.
J. Aircr.
47
,
141
-
150
.
Dadashi
,
S.
,
Feaster
,
J.
,
Bayandor
,
J.
,
Battaglia
,
F.
and
Kurdila
,
A. J.
(
2016
).
Identification and adaptive control of history dependent unsteady aerodynamics for a flapping insect wing
.
Nonlinear Dyn.
85
,
1405
-
1421
.
Dickinson
,
M. H.
,
Lehmann
,
F.-O.
and
Sane
,
S. P.
(
1999
).
Wing rotation and the aerodynamic basis of insect flight
.
Science
284
,
1954
-
1960
.
Du
,
G.
and
Sun
,
M.
(
2008
).
Effects of unsteady deformation of flapping wing on its aerodynamic forces
.
Appl. Math. Mech.
29
,
731
-
743
.
Du
,
G.
and
Sun
,
M.
(
2012
).
Aerodynamic effects of corrugation and deformation in flapping wings of hovering hoverflies
.
J. Theor. Biol.
300
,
19
-
28
.
Du
,
L.
and
Sun
,
X.
(
2015
).
Effect of flapping frequency on aerodynamics of wing in freely hovering flight
.
Comput. Fluids
117
,
79
-
87
.
Ellington
,
C. P.
(
1984
).
The aerodynamics of hovering insect flight. III. Kinematics
.
Philos. Trans. R. Soc. B Biol. Sci.
305
,
41
-
78
.
Ellington
,
C. P.
(
1985
).
Power and efficiency of insect flight muscle
.
J. Exp. Biol.
115
,
293
-
304
.
Ellington
,
C. P.
and
Dudley
,
R.
(
1990
).
Mechanics of forward flight in bumblebees. II: Quasi-steady lift and power requirements
.
J. Exp. Biol.
148
,
19
-
52
.
Erzincanli
,
B.
and
Sahin
,
M.
(
2015
).
The numerical simulation of the wing kinematics effects on near wake topology and aerodynamic performance in hovering Drosophila flight
.
Comput. Fluids
122
,
90
-
110
.
Hamdani
,
H. R.
and
Naqvi
,
A.
(
2011
).
A study on the mechanism of high-lift generation by an insect wing in unsteady motion at small Reynolds number
.
Int. J. Numer. Methods Fluids
67
,
581
-
598
.
Hu
,
J. X.
and
Xiao
,
Q.
(
2014
).
Three-dimensional effects on the translational locomotion of a passive heaving wing
.
J. Fluid Struct.
46
,
77
-
88
.
Huang
,
H.
,
Sun
,
M.
and
Huang
,
H.
(
2007
).
Dragonfly forewing-hindwing interaction at various flight speeds and wing phasing
.
AIAA J.
45
,
508
-
511
.
Ishihara
,
D.
,
Horie
,
T.
and
Denda
,
M.
(
2009
).
A two-dimensional computational study on the fluid–structure interaction cause of wing pitch changes in dipteran flapping flight
.
J. Exp. Biol.
212
,
1
-
10
.
Ishihara
,
D.
,
Horie
,
T.
and
Niho
,
T.
(
2014
).
An experimental and three-dimensional computational study on the aerodynamic contribution to the passive pitching motion of flapping wings in hovering flies
.
Bioinspir. Biomim.
9
,
046009
.
Khurana
,
M.
and
Chahl
,
J.
(
2014
).
Optimization of the leading edge segment of a corrugated wing
.
Bioinspir. Biomimet. Biorepl.
9055
,
905517
.
Kim
,
D.
and
Choi
,
H.
(
2007
).
Two-dimensional mechanism of hovering flight by single flapping wing
.
J. Mech. Sci. Technol.
21
,
207
-
221
.
Kolomenskiy
,
D.
,
Moffatt
,
H. K.
,
Farge
,
M.
and
Schneider
,
K.
(
2011
).
Two- and three-dimensional numerical simulations of the clap–fling–sweep of hovering insects
.
J. Fluid Struct.
27
,
784
-
791
.
Lee
,
J.-S.
,
Kim
,
J.-H.
and
Kim
,
C.
(
2008
).
Numerical study on the unsteady-force-generation mechanism of insect flapping motion
.
AIAA J.
46
,
1835
-
1848
.
Lee
,
K.-B.
,
Kim
,
J.-H.
and
Kim
,
C.
(
2011
).
Aerodynamic effects of structural flexibility in two-dimensional insect flapping flight
.
J. Aircr.
48
,
894
-
909
.
Liu
,
H.
(
2009
).
Integrated modeling of insect flight: from morphology, kinematics to aerodynamics
.
J. Comput. Phys.
228
,
439
-
459
.
Liu
,
H.
,
Ellington
,
C. P.
,
Kawachi
,
K.
,
Berg
,
C.
,
den
, and
Willmott
,
A. P.
(
1998
).
A computational fluid dynamic study of hawkmoth hovering
.
J. Exp. Biol.
201
,
461
-
477
.
Lu
,
X.-Y.
and
Liao
,
Q.
(
2006
).
Dynamic responses of a two-dimensional flapping foil motion
.
Phys. Fluid.
18
,
098104
.
Luo
,
G.
and
Sun
,
M.
(
2005
).
Effects of corrugation and wing planform on the aerodynamic force production of sweeping model insect wings
.
Acta Mech. Sin.
21
,
531
-
541
.
Ma
,
K. Y.
,
Felton
,
S. M.
and
Wood
,
R. J.
(
2012
).
Design, fabrication, and modeling of the split actuator microrobotic bee
.
Intelligent Robots and Systems (IROS), 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems
1133
-
1140
.
Mao
,
S.
and
Gang
,
D. U.
(
2003
).
Lift and power requirements of hovering insect flight
.
Acta Mech. Sin.
19
,
458
-
469
.
Mazaheri
,
K.
and
Ebrahimi
,
A.
(
2010
).
Experimental investigation of the effect of chordwise flexibility on the aerodynamics of flapping wings in hovering flight
.
J. Fluid Struct.
26
,
544
-
558
.
Meng
,
X.
and
Sun
,
M.
(
2011
).
Aerodynamic effects of corrugation in flapping insect wings in forward flight
.
J. Bionic Eng.
8
,
140
-
150
.
Meng
,
X. G.
,
Xu
,
L.
and
Sun
,
M.
(
2011
).
Aerodynamic effects of corrugation in flapping insect wings in hovering flight
.
J. Exp. Biol.
214
,
432
-
444
.
Miyake
,
H.
(
2008
).
Micro-structural Approach to Developing the Resonance Model of the Indirect Flight Mechanism. In Bio-mechanisms of Swimming and Flying (ed. N. Kato and S. Kamimura), pp. 129-139. Tokyo: Springer
.
Morse
,
D. R.
and
Liburdy
,
J. A.
(
2009
).
Vortex dynamics and shedding of a low aspect ratio, flat wing at low Reynolds numbers and high angles of attack
.
ASME J. Fluids Eng.
131
,
051202
.
Nakata
,
T.
and
Liu
,
H.
(
2012
).
A fluid-structure interaction model of insect flight with flexible wings
.
J. Comput. Phys.
231
,
1822
-
1847
.
Nakata
,
T.
,
Liu
,
H.
,
Tanaka
,
Y.
,
Nishihashi
,
N.
,
Wang
,
X.
and
Sato
,
A.
(
2011
).
Aerodynamics of a bio-inspired flexible flapping-wing micro air vehicle
.
Bioinspir. Biomim.
6
,
045002
.
New
,
T. H.
,
Chan
,
Y. X.
,
Koh
,
G. C.
,
Hoang
,
M. C.
and
Shi
,
S.
(
2014
).
Effects of corrugated aerofoil surface features on flow-separation control
.
AIAA J.
52
,
206
-
211
.
Ramamurti
,
R.
and
Sandberg
,
W. C.
(
2002
).
A three-dimensional computational study of the aerodynamic mechanisms of insect flight
.
J. Exp. Biol.
205
,
1507
-
1518
.
Saha
,
S. K.
and
Celata
,
G. P.
(
2011
).
Advances in modelling of biomimetic fluid flow at different scales
.
Nanoscale Res. Lett.
6
,
344
.
Saputra
,
Byun
,
D.
,
Byun
,
Y. H.
and
Park
,
H. C.
(
2006
).
Experimental and numerical study on flapping wing kinematics and aerodynamics of coleoptera
.
Exp. Mech. Nano Biotechnol.
326-328
,
175
-
178
.
Senda
,
K.
,
Hirai
,
N.
,
Iima
,
M.
,
Yokoyama
,
N.
,
Obara
,
T.
,
Kobayakawa
,
T.
and
Yokoi
,
K.
(
2010
).
Analysis of flapping flight of butterfly based on experiments and numerical simulations
.
World Automation Congress (WAC), 2010
1
-
5
.
Shen
,
C.
and
Sun
,
M.
(
2015
).
Power requirements of vertical flight in the dronefly
.
J. Bionic Eng.
12
,
227
-
237
.
Shyy
,
W.
,
Aono
,
H.
,
Chimakurthi
,
S. K.
,
Trizila
,
P.
,
Kang
,
C.-K.
,
Cesnik
,
C. E. S.
and
Liu
,
H.
(
2010
).
Recent progress in flapping wing aerodynamics and aeroelasticity
.
Prog. Aerosp. Sci.
46
,
284
-
327
.
Singh
,
S. N.
,
Simha
,
A.
and
Mittal
,
R.
(
2004
).
Biorobotic AUV maneuvering by pectoral fins: Inverse control design based on CFD parameterization
.
IEEE J. Ocean. Eng.
29
,
777
-
785
.
Sudhakar
,
Y.
and
Vengadesan
,
S.
(
2010
).
Flight force production by flapping insect wings in inclined stroke plane kinematics
.
Comput. Fluids
39
,
683
-
695
.
Sun
,
M.
and
Lan
,
S. L.
(
2004
).
A computational study of the aerodynamic forces and power requirements of dragonfly (Aeschna juncea) hovering
.
J. Exp. Biol.
207
,
1887
-
1901
.
Taha
,
H. E.
,
Hajj
,
M. R.
and
Beran
,
P. S.
(
2014
).
State-space representation of the unsteady aerodynamics of flapping flight
.
Aerosp. Sci. Technol.
34
,
1
-
11
.
Taylor
,
G. K.
,
Nudds
,
R. L.
and
Thomas
,
A. L. R.
(
2003
).
Flying and swimming animals cruise at a Strouhal number tuned for high power efficiency
.
Nature
425
,
707
-
711
.
Tian
,
F.-B.
,
Dai
,
H.
,
Luo
,
H. X.
,
Doyle
,
J. F.
and
Rousseau
,
B.
(
2014
).
Fluid-structure interaction involving large deformations: 3D simulations and applications to biological systems
.
J. Comput. Phys.
258
,
451
-
469
.
Walker
,
J. A.
(
2002
).
Rotational lift: something different or more of the same?
J. Exp. Biol.
205
,
3783
-
3792
.
Walker
,
S. M.
,
Schwyn
,
D. A.
,
Mokso
,
R.
,
Wicklein
,
M.
,
Müller
,
T.
,
Doube
,
M.
,
Stampanoni
,
M.
,
Krapp
,
H. G.
and
Taylor
,
G. K.
(
2014
).
In vivo time-resolved microtomography reveals the mechanics of the blowfly flight motor
.
PLoS Biol.
12
,
e1001823
.
Wang
,
Z. J.
(
2000
).
Two dimensional mechanism for insect hovering
.
Phys. Rev. Lett.
85
,
2216
.
Wang
,
Z. J.
,
Birch
,
J. M.
and
Dickinson
,
M. H.
(
2004
).
Unsteady forces and flows in low Reynolds number hovering flight: two-dimensional computations vs robotic wing experiments
.
J. Exp. Biol.
207
,
449
-
460
.
Winston
,
M. L.
(
1991
).
The Biology of the Honey Bee
. Cambridge, Massachusetts:
Harvard University Press
.
Wood
,
R.
,
Nagpal
,
R.
and
Wei
,
G.-Y.
(
2013
).
Flight of the Robobees
.
Sci. Am.
308
,
60
-
65
.
Wu
,
J.-S.
and
Shao
,
Y.-L.
(
2004
).
Simulation of lid-driven cavity flows by parallel lattice Boltzmann method using multi-relaxation-time scheme
.
Int. J. Numer. Methods Fluids
46
,
921
-
937
.
Wu
,
J. H.
and
Sun
,
M.
(
2005
).
The influence of the wake of a flapping wing on the production of aerodynamic forces
.
Acta Mech. Sin.
21
,
411
-
418
.
Yan
,
J.
,
Wood
,
R. J.
,
Avadhanula
,
S.
,
Sitti
,
M.
and
Fearing
,
R. S.
(
2001
).
Towards flapping wing control for a micromechanical flying insect
. In
Robotics and Automation, 2001. Proceedings 2001 ICRA. IEEE International Conference on Robotics and Automation
. Vol.
4
,
3901
-
3908
.
Yin
,
B.
and
Luo
,
H.
(
2010
).
Effect of wing inertia on hovering performance of flexible flapping wings
.
Phys. Fluids
22
.
Zhu
,
J. Y.
and
Zhou
,
C. Y.
(
2014
).
Aerodynamic performance of a two-dimensional flapping wing in asymmetric stroke
.
J. Aerosp. Eng.
228
,
641
-
651
.
Zuo
,
D.
,
Peng
,
S.
,
Chen
,
W.
and
Zhang
,
W.
(
2007
).
Numerical simulation of flapping-wing insect hovering flight at unsteady flow
.
Int. J. Numer. Methods Fluids
53
,
1801
-
1817
.

Competing interests

The authors declare no competing or financial interests.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.

Supplementary information