A computational study on the influence of insect wing geometry on bee flight mechanics

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.

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.

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.

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.

Fig. 4 presents the C_L 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 C_L, 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 C_L. 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 C_L at τ=11.5. As V_∞ increases to 5 m/s, the elliptical wing profile predicts similar C_L 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, C_L for the bee wing and the approximate wing remains in good agreement. The peak C_L is substantially higher for the flat plate and the minimum C_L 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 C_L.

Comparing time-averaged C_L, in Table 1, all cross-sections except for the elliptic one exhibit a significant increase in C_L at 2 m/s. After 2 m/s, the C_L 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 C_L from 2 to 5 m/s, while the other cross-sections display a progressive decrease. The approximate cross-section C_L is consistently closest to the bee results, with a maximum difference of 0.03.

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.

Fig. 5 presents the C_D data for four flap cycles comparing each wing cross-section versus τ at V_∞=1−5 m/s. The C_D 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 C_D at τ=11.2. After τ=11.2, C_D increases steadily to the maximum around τ=11.7. After the maximum at τ=11.7, there is another drop in C_D at τ=11.9. From τ=11.9 to τ=12, C_D increases again, peaking immediately after transitioning back to the upstroke.

At 3 m/s, the C_D 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 C_D 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 C_D for the approximate and true bee wing. The minimum C_D 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.

The similarity in the vorticity and pressure contours manifest in the similarities in C_L and C_D 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 C_L and C_D 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.

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 C_L and C_D, time-averaged C_L and C_D, and RMSD of C_L and C_D for each cross-section to gain a more quantitative representation of the instantaneous variation in C_D and C_L.

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 C_D, 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.

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 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, V_mag 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 C_D 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).

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/m³ and a kinematic viscosity of 1.15×10^–4 m²/s. The kinematics used in the simulation are f=0.25 Hz, r=0.1625 m, θ_1a=23.5°, θ_2a=0°, θ_3a=45°, φ₁=π/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. 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.

In addition to the DPIV measurements, a 2D force sensor was attached to the base of the wing to capture the transient C_L and C_D 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 C_L and C_D 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.

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.

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 (θ₁), roll (θ₂) and wing pitch (θ₃). The angular velocities θ₁, θ₂ and θ₃ 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 θ₂ 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=ρ(V_kin+V_∞)c/μ. Due to the nature of the problem, two different velocities are possible: the flow over the wing driven by the kinematics (V_kin), 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.

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 θ₃ 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.

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 (C_L and C_D) 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.

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.