2.1 Experimental systems and techniques

Encarsia formosa Gahan were acquired from the Laboratory of Biological Invasion of Institute of Plant Protection, Chinese Academy of Agricultural Sciences, which were descendants of wild-caught E. formosa. The wasps were housed in small chambers in groups of three to five individuals and fed with 20 % sugar solution. The flight of the wasps in a transparent flight chamber ( $40\times 40\times 30~\text{mm}^{3}$ ) was filmed using three orthogonally aligned synchronized high-speed cameras (FASTCAM Mini UX100, Photron Inc., San Diego, CA, USA; 10 000 frames per second, shutter speed $20~\unicode[STIX]{x03BC}\text{s}$ , resolution $896\times 488$  pixels) mounted on an optical table (figure 1). Each camera was equipped with a 60 mm micro-Nikkor lens and 12 mm extension tube. Each camera view was backlit using a 50 W integrated red light emitting diode (LED; luminous flux, 4000 lm; wavelength, 632 nm). The light was made uniform by two lenses. We manually triggered the synchronized cameras when the insect was observed to fly steadily in the filming area of approximately $4\times 4\times 4~\text{mm}^{3}$ , which represented the intersecting field of views of the three cameras. The experiment was conducted at room temperature $25{-}27\,^{\circ }\text{C}$ and relative humidity 50–60 %.

Figure 1.

A sketch showing the experimental set-up: the flight chamber and cameras; each camera view is backlit using a integrated red light emitting diode (LED).

The method used to extract the three-dimensional (3-D) body and wing kinematics from the filmed data was a modification of that used in several previous works of our group (Liu & Sun Reference Liu and Sun2008; Mou, Liu & Sun Reference Mou, Liu and Sun2011). The original method and the modification are described in appendix A.

2.2 Methods of flow and aerodynamic-force computation

The flows around and the aerodynamic force acting on the insect were computed using the method of computational fluid dynamics. The governing equations of the flow are the incompressible Navier–Stokes equations, the dimensionless forms of which are

(2.1) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\left(\frac{1}{Re}\right)\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}$ is the non-dimensional fluid velocity, $p$ the non-dimensional fluid pressure, $\unicode[STIX]{x1D70F}$ the non-dimensional time, $\unicode[STIX]{x1D735}$ the gradient operator, $\unicode[STIX]{x1D6FB}^{2}$ the Laplacian operator and $Re$ the Reynolds number; in the non-dimensionalization, $c$ , $U$ and $c/U$ are taken as reference length, velocity and time, respectively (here $c$ is the mean chord length of wing and $U$ is the mean flapping velocity at the radius of gyration of wing).

In the present study, the linked fore- and hindwings were modelled as a single wing; but the forewing could rotate relative to the hindwing around the ‘leading edge’ of the hindwing, which happened during the ‘fling’. The wing of E. formosa consists of the membrane and the brim of marginal hairs (figure 2 a), but here the hairs that line the wing are modelled as an extension of the membranous wing, i.e. the planform of the model wing is similar to that of a wing consisting of the membrane and the brim of marginal hairs (figure 2 b). The effects and limitations of these simplifying assumptions will be discussed in a later section (§ 7). The section thickness of the model wing is 3 % of the local chord length and the leading and trailing edges are rounded. The body shape is given by fitting ellipses to the body in the video images.

Figure 2.

(a) Wing planform of an Encarsia formosa, extracted from microscope images. Blue circles indicate the boundary between fore- and hindwings. (b) The planform of the model wing used in computations. The red dot indicates the wing base and the dashed line indicates the axis of pitch rotation.

Figure 3.

Portions of computational grids.

The equations were solved over moving overset grids because there are relative movements between the left and right wings and between the body and wing. There was a body-fitted curvilinear grid for each of the wings and the body and a background Cartesian grid which extends to the far-field boundary of the domain (figure 3). In the calculation process, the wing grid was regenerated in every time step because of the spanwise bending deformation and the rotation of the forewing relative to the hindwing during the ‘fling’. The computational method is the same as that used in several previous studies of our group (Sun & Yu Reference Sun and Yu2003, Reference Sun and Yu2006; Du & Sun Reference Du and Sun2012) and only an outline of the method is given here. The algorithm is based on the method of artificial compressibility. It was first developed by Rogers, Kwak & Kiris (Reference Rogers, Kwak and Kiris1991) for a single grid, and then extended by Rogers & Pulliam (Reference Rogers and Pulliam1994) to overset grids. The time derivatives of the momentum equations are differenced using a second-order, three-point backward difference formula. To solve the time discretized momentum equations for a divergence-free velocity at a new time level, a pseudo-time level is introduced into the equations and a pseudo-time derivative of pressure divided by an artificial compressibility constant is introduced into the continuity equation. The resulting system of equations is iterated in pseudo-time until the pseudo-time derivative of pressure approaches zero; thus, the divergence of the velocity at the new time level approaches zero. The derivatives of the viscous fluxes in the momentum equation are approximated using second-order central differences. For the derivatives of convective fluxes, upwind differencing based on the flux-difference splitting technique is used. A third-order upwind differencing is used at the interior points and a second-order upwind differencing is used at points next to boundaries. With overset grids, the solution method for a single grid is applied to each of the wing grid and the background grid, and data are interpolated from one grid to another at the inter-grid boundary points using tri-linear interpolation. Details of this algorithm can be found in Rogers et al. (Reference Rogers, Kwak and Kiris1991) and Rogers & Pulliam (Reference Rogers and Pulliam1994). For the far-field boundary conditions, at the inflow boundary, the velocity components are specified according to relative velocity at the boundary while pressure is extrapolated from the interior; at the outflow boundary, pressure is set equal to the static pressure of the still air and the velocity is extrapolated from the interior. On the wing surfaces, impermeable wall and no-slip boundary conditions are applied, and the pressure on the boundary is obtained through the normal component of the momentum equation written in the moving coordinate system. The background Cartesian grid was generated algebraically. The wing grid (O–H type) and body grid (O–O type) were generated using a Poisson solver which was based on the work of Hilgenstock (Reference Hilgenstock, Sengupta, Hauster, Eiseman and Thompson1988). The grids will be further described in § 4.1, as will be the analysis of the convergence of solutions.

4.1 Code validation and grid test

The flow solver had been tested by unsteady aerodynamic forces measured by translating (Meng & Sun Reference Meng and Sun2013) and rotating (Sun & Yu Reference Sun and Yu2006) model insect wings. These tests showed that the unsteady aerodynamic forces computed by the present solver agreed well with the experimental measurements. But these tests were for cases of higher Reynolds number ( $Re$ above 200). Here we further test the solver at Reynolds numbers ( $Re\sim 10$ ) identical to the wasp flight. First, we made a conceptual test for a single, flapping wing at $Re=10$ (the flapping motion is an idealized flapping motion often used in insect flight study, given by Dickinson et al. Reference Dickinson, Lehmann and Sane1999). We computed the flow first using the single grid, and then computed the flow using the multi-block, overset grids. Since the physical problem is the same in the two computations, the results should be the same. This can be a test on the multi-block, overset solver. Figure 8(a) shows the comparison in the lift and drag coefficients ( $C_{L}$ and $C_{D}$ ) between the two computations. There is very good agreement between these two sets of results, showing that the multi-block, overset Navier–Stokes solver works well at $Re$ around 10.

Figure 8.

(a) Lift and drag coefficients of a single wing flapping at $Re=10$ computed using a single grid, compared with those computed using the multi-block, overset grids. (b) The calculated drag coefficient of a sphere at $Re=0.5{-}20$ , compared with the theoretical results and experimental data.

To further validate the low- $Re$ simulation, we computed the flows around a sphere at $Re=0.5{-}20$ , for which there exist Stokes or Oseen solutions and accurate experiment data (White Reference White1991). Figure 8(b) compares the computed drag coefficient with the theoretical and experimental data. The difference between the computed results and the data is less than 2 %.

As mentioned above, the flow solver was tested by measured data of translating and rotating wings. Recently, lift coefficient ( $C_{L}$ ) was measured for model insect wings in flapping motion at $Re=180$ (Han, Chang & Han Reference Han, Chang and Han2016) and in revolving motion with an acceleration phase, a constant-speed phase and a deceleration phase at $Re=500$ (Manar & Jones Reference Manar and Jones2014). Here we make a further test of the solver using these experiment data. As seen in figure 9, the computed $C_{L}$ agrees very well with the measured data.

Figure 9.

Comparison between the calculated and measured lift coefficients on a rectangular wing in revolving motion with an acceleration phase, a constant-speed phase and a deceleration phase ( $Re=500$ ) (a), and on a model fruitfly wing in flapping motion ( $Re=180$ ) (b).

Before proceeding to compute the flows, grid resolution tests were conducted to ensure that the flow calculation was grid independent. Three grid systems were considered. For grid system 1, the wing grid had dimensions $41\times 61\times 43$ in the normal direction, around the wing, and in the spanwise direction, respectively (first layer grid thickness was $0.0015c$ ); the body grid had dimensions $66\times 51\times 35$ along the body, in the azimuthal direction and in the normal direction, respectively; the background grid had dimensions $81\times 81\times 81$ in the  $x$ , $y$ and $z$ directions, respectively. For grid system 2, the corresponding grid dimensions were $61\times 91\times 65$ , $99\times 77\times 53$ and $121\times 121\times 121$ ( $0.001c$ ). For grid system 3, the corresponding grid dimensions were $91\times 135\times 96$ , $149\times 118\times 80$ and $181\times 181\times 181$ ( $0.00067c$ ). For all three grid systems, the outer boundary of the background grid was $20c$ from the wing root; the grid points of the background grid concentrated in the near field of the wings where its grid density was approximately the same as that of the outer part of the body grid. The typical cell size in the near wake, where vortices from a stroke may interact with the wing in the subsequent strokes, was less than $0.038c$ , $0.061c$ and $0.1c$ for grid systems 3, 2 and 1, respectively. The aerodynamic forces and vorticity patterns computed from three grid systems are shown in figure 10 (the wings started to flap in a quiescent fluid, and after about three cycles the flow became periodical). It is observed that there is almost no difference between the forces calculated by the three grid systems (figure 10 a). The first grid refinement produces a little change in the vorticity plot, and the second grid refinement produces almost no change (figure 10 b). Calculations were also conducted using a larger computational domain. The domain was enlarged by adding more grid points to the outside of the background grid of grid system 2 (the outer boundary of the background grid of grid system 2 was increased from $20c$ to  $40c$ ). The calculated results showed that there was no need to put the outer boundary further than that of grid system 2. The non-dimensional time step was 0.02 (non-dimensionalized by $c/U$ ). The effect of time step value was studied and it was found that a numerical solution effectively independent of the time step was achieved if the time step value was ${\leqslant}0.02$ . From the above discussion, it was concluded that grid system 2 and time step 0.02 were appropriate for the calculation.

Figure 10.

Grid resolution test. (a) Time courses of the vertical and horizontal forces of the wings in one cycle, calculated with three grid systems. (b) The corresponding non-dimensional spanwise vorticity contours at section $0.7R$ at $t/T=0.14$ ; the magnitude of the non-dimensional vorticity at the outer contour is 3 and the contour interval is 1.

4.2 Aerodynamic forces

In the present study, the vertical and horizontal components of the total aerodynamic force of a wing are referred to as vertical force and horizontal force, respectively. Because the wing does not translate in the horizontal plane, the vertical force is not the lift of the wing and the horizontal force is not the drag. The velocity at the radius of gyration ( $r_{2}$ ) of the wing is used to represent the velocity of the wing (figure 11 e). Lift and drag are defined as the components of the total aerodynamic force that are perpendicular and parallel to the velocity of the wing, respectively. Figure 11(a–d) gives the time histories of the computed aerodynamic forces of a wing in one cycle for EF1 (results for the other seven individuals are given in appendix B, see figure 23). Because the computed forces show great consistency for the eight individuals, we analyse the results of EF1 (figure 11) as an example. As seen from the figure, a large vertical-force peak is produced during the rowing (figure 11 a, $t/T\approx 0{-}0.2$ ) and another smaller one during the ‘fling’ (figure 11 a, $t/T\approx 0.55{-}0.7$ ); the other parts of the cycle do not produce any positive vertical force. From the data in figure 11(a), it is calculated that the rowing contributes 70 % of the total vertical force and the ‘fling’ the other 30 %. The rowing produces a forward pointing horizontal force and the ‘fling’ and the subsequent forward sweeping produces a backward pointing horizontal force (figure 11 b), and the two approximately cancel, as required for hovering and ascending flight.

Figure 11.

Aerodynamic forces. (a–d) Vertical force, horizontal force, lift and drag of a wing. (e) Diagram showing the wing position, velocity vector and force vector at various times in one stroke cycle (side view). (f) Back view of the wing in ‘fling’ motion, represented by the wing section at $0.7R$ , and the corresponding velocity vector and force vector. Blue arrow, velocity; red arrow, force.

As mentioned, the insects are in hovering or very slow ascending flight. Thus the computed mean force vector of the wings should be vertical. As seen in table 1, the mean force vector is very close to being vertical: its horizontal component ( $F_{w,h}$ ) is small and the angle between the force vector and the vertical ( $\unicode[STIX]{x1D702}$ ) is less than $6^{\circ }$ for five of the eight wasps and approximately $10^{\circ }$ for the other three. The magnitude of the mean force vector of the wings ( $F_{w}$ ) is shown in table 1. When in slow ascending flight, there is a small drag on the body (the computed body drag, $D_{b}$ , is also given in table 1). The value of ( $F_{w}$ – $D_{b}$ ) should be equal to the weight of the insect. As mentioned above, we could not measure the weight of the insects in our experiment. Weis-Fogh (Reference Weis-Fogh1973) measured the weight of ten wasps of the same species as those in our experiment (E. Formosa) and gave the average weight of the insects. Since the wasps in our experiment and those in Weis-Fogh’s experiment are the same species, they are geometrically similar. Thus we can make an indirect comparison here. For the insects in Weis-Fogh’s experiment, the average wing length was 0.62 mm, the average wing-tip speed was $1.17~\text{m}~\text{s}^{-1}$ and the average weight was $0.245~\unicode[STIX]{x03BC}\text{N}$ . For the eight insects in our experiment, the average wing length is 0.6 mm and the average wing-tip speed is $0.1075~\text{m}~\text{s}^{-1}$ . Assuming that the aerodynamic force is proportional to the wing area (or the square of the wing length) and to the square of wing-tip speed, the average weight of the insects in our experiment can be estimated as

(4.1) $$\begin{eqnarray}\left({\textstyle \frac{0.6}{0.62}}\right)^{2}\times \left({\textstyle \frac{1.075}{1.17}}\right)^{2}\times 0.245~\unicode[STIX]{x03BC}\text{N}=0.194~\unicode[STIX]{x03BC}\text{N}.\end{eqnarray}$$

The computed average $F_{w}$ – $D_{b}$ for the eight insects is $0.196~\unicode[STIX]{x03BC}\text{N}$ , in good agreement with the above estimated average weight ( $0.194~\unicode[STIX]{x03BC}\text{N}$ ).

4.3 Aerodynamic-force mechanisms

The vertical and horizontal forces of a wing come from its lift and drag. The rowing produces a large drag and a much smaller lift (figure 11 c,d). Therefore, the large vertical force by the rowing is mainly from the drag of the wing (approximately 70 % of the vertical force during the rowing period is contributed by the drag). In this period, the wing is in a transient motion (velocity increases from zero at $t/T=0$ to a large value at $t/T=0.15$ and then decreases to a small value at $t/T\approx 0.2$ ) and the angle of attack is rather large, between $46^{\circ }$ and $65^{\circ }$ (figure 11 e). The transient motion and force here are very similar to those of an oar of a boat. Let us examine how the large drag is produced.

Since the wing is in fast acceleration ( $t/T=0{-}t/T\approx 0.15$ ) and deceleration ( $t/T\approx 0.15{-}t/T\approx 0.2$ ), we first look at the role of added mass. Howe (Reference Howe1989, Reference Howe1995) made an analysis of the force exerted on a rigid body in unsteady motion in a uniform incompressible viscous fluid, and showed that the force on the body is the summation of three parts. The first part is caused by the acceleration of the body and the expression of this part is identical to the added-mass force obtained by solving the potential flow of no free vorticity, and the second and third parts are due to the free vorticity in the wake and the viscous skin friction, respectively. His results show that the added mass of a rigid body moving in an incompressible fluid can be computed by solving a separate potential flow. Our group recently developed a code for computing the added mass force of flapping wings (Liu & Sun Reference Liu and Sun2018). Using the code we calculated the added-mass force during the rowing period and the results are shown in figure 12(a,b). It is seen that the added-mass forces are relatively small. The large drag must be produced by other mechanisms.

Figure 12.

(a,b) Added-mass contributions to the vertical and horizontal forces. (c,d) The contributions of the pressure- and friction-based components to the vertical and horizontal forces.

Figure 12(c,d) shows the contributions of the pressure- and friction-based components to the vertical and horizontal forces. It is seen that the large vertical force or the large drag of the wing is mainly contributed by the pressure-based component (figure 12 c). Let us investigate how the large pressure-based force is generated. Flow field data in this period are shown in figure 13. The fluid is pushed downwards by the ‘lower’ wing surface and sucked downwards by the ‘upper’ wing surface (figure 13 c). Strong vortices are shed at the leading and trailing edges (figure 13 a,b). The pushing results in a large pressure on the lower surface of the wing and the induced velocity of the leading-edge vortex (LEV) and trailing-edge vortex (TEV) results in a large suction pressure on the upper surface of the wing (figure 13 d), giving the large drag or vertical force. Note that for the tiny insects here, both the LEV and TEV are used in enhancing the aerodynamic force production, unlike the cases for larger insects whose $Re$ is larger than 100, where LEV dominates the enhancement of aerodynamic force production. We thus see that the large pressure force that gives the large drag or vertical force is generated by the unsteady effect similar to that of starting flow of a plate started in its normal direction.

Figure 13.

Flow fields and surface pressure distributions during the rowing period. (a) Iso-vorticity surface plots (top view) at $t/T=0.09$ , 0.12, 0.15 and 0.18 (in this and other plots of iso-vorticity surfaces, the magnitude of the non-dimensional vorticity is 2). Vorticity is non-dimensionalized by $U/c$ . (b) Non-dimensional spanwise vorticity contours at section  $0.7R$ ; solid and dashed lines denote the anticlockwise and clockwise vorticity, respectively. The magnitude of the non-dimensional vorticity at the outer contour is 1 and the contour interval is 1. (c) Non-dimensional velocity vector plot at section  $0.7R$ . (d) Non-dimensional surface pressure distributions at section  $0.7R$ . Non-dimensional pressure, $C_{P}$ , is defined as $(p-p_{\infty })/0.5\unicode[STIX]{x1D70C}U^{2}$ .

To further show the unsteady effect of the transient motion on the aerodynamic forces in the rowing period, we computed the corresponding steady flows of the wing, i.e. at each time point in the rowing period, we use the instantaneous speed and angle of attack of the wing at that time to compute a steady-state flow. For example, at $t/T=0.09$ in figure 13, the wing’s tip speed and angle of attack were known; we let the wing move (revolve) at the given tip speed and angle of attack until the flow reached steady state (after the wing revolved for more than one revolution, the flow became steady), thus we obtain the corresponding steady flow and force for this time point. We do the same for other time points, $t/T=0.12$ , 0.15 and 0.18 and so on. The steady-state results are shown in figure 14. It is seen that in the steady-state cases, the vortices shed at the leading and trailing edges are much weaker and the pressure on the lower surface and the suction on the upper surface are much smaller than the corresponding ones in the transient motion (comparing results in figure 14 with those in figure 13). The reason for the LEV and TEV in the unsteady case being much stronger than those in the steady-state case is as follows. The wing in the unsteady case is similar to a plate in starting motion in its normal direction, in which vorticity shed at the edges of the plate could accumulate to a large amount behind the plate. We computed the vertical force from the steady-state case and it is 45 % less than that produced by the unsteady (rowing) motion.

Figure 14.

Flow fields and surface pressure distributions for the quasi-steady case, corresponding to those in figure 13. (Data here are non-dimensionalized using the same velocity as that used in figure 13.)

Next, let us look at the vertical-force peak during the ‘fling’. It is also due to the large drag on the wing (see figure 11 a,d,f). The production mechanism of the large drag has been explained in detail by many previous studies (Spedding & Maxworthy Reference Spedding and Maxworthy1986; Sun & Yu Reference Sun and Yu2003; Miller & Peskin Reference Miller and Peskin2005; Arora et al. Reference Arora, Gupta, Sanghi, Aono and Shyy2014). Flow field data during the ‘fling’ are shown in figure 15. As pointed out in previous studies (Spedding & Maxworthy Reference Spedding and Maxworthy1986; Sun & Yu Reference Sun and Yu2003; Miller & Peskin Reference Miller and Peskin2005), the mechanism of the vertical-force peak during the ‘fling’ is that the opening of the wing pair generates a vortex ring containing a strong downward jet (see figure 15 a,c) in a very short period, or a large fluid impulse, resulting in the large force.

Figure 15.

Flow fields during the ‘fling’ period. (a) Iso-vorticity surface plots (top view) at $t/T=0.55$ , 0.60, 0.65 and 0.70. (b) Non-dimensional spanwise vorticity contours at section  $0.7R$ . Solid and dashed lines denote the anticlockwise and clockwise vorticity, respectively. (c) Non-dimensional velocity vector plot at section  $0.7R$ .