Aerodynamic Ground Effect in Fruitfly Sized Insect Takeoff

Aerodynamic ground effect in flapping-wing insect flight is of importance to comparative morphologies and of interest to the micro-air-vehicle (MAV) community. Recent studies, however, show apparently contradictory results of either some significant extra lift or power savings, or zero ground effect. Here we present a numerical study of fruitfly sized insect takeoff with a specific focus on the significance of leg thrust and wing kinematics. Flapping-wing takeoff is studied using numerical modelling and high performance computing. The aerodynamic forces are calculated using a three-dimensional Navier–Stokes solver based on a pseudo-spectral method with volume penalization. It is coupled with a flight dynamics solver that accounts for the body weight, inertia and the leg thrust, while only having two degrees of freedom: the vertical and the longitudinal horizontal displacement. The natural voluntary takeoff of a fruitfly is considered as reference. The parameters of the model are then varied to explore possible effects of interaction between the flapping-wing model and the ground plane. These modified takeoffs include cases with decreased leg thrust parameter, and/or with periodic wing kinematics, constant body pitch angle. The results show that the ground effect during natural voluntary takeoff is negligible. In the modified takeoffs, when the rate of climb is slow, the difference in the aerodynamic forces due to the interaction with the ground is up to 6%. Surprisingly, depending on the kinematics, the difference is either positive or negative, in contrast to the intuition based on the helicopter theory, which suggests positive excess lift. This effect is attributed to unsteady wing-wake interactions. A similar effect is found during hovering.


Introduction
The aerodynamic forces of an air vehicle or an animal may be affected by the ground proximity. This phenomenon, known as the ground effect, has been extensively studied for aircraft [1] and rotorcraft [2]. Although the effect varies depending on many design parameters, the general trend is an increase in lift and pitching moment, and a decrease in drag. The effect decays as the distance from the ground increases, and vanishes at a distance slightly larger than the characteristic length of the vehicle. For example, for a hovering helicopter, the excess thrust vanishes if the distance to the ground exceeds 1.25 times the diameter of the main rotor [2].
Rayner [3] proposed a fixed wing lifting line theory for forward flight of birds, bats and insects. His analysis suggested that flight in ground effect provides performance improvements, if the flight speed is not too low. However, this theory could not be applied to hovering or slow forward flight at very low height, since it neglected flapping motion. Normal hovering in ground effect was considered by Gao and Lu [4]. They carried out two-dimensional numerical simulations of hovering and identified three regimes: force enhancement, force reduction, and force recovery, depending on the distance from the ground. Liu et al. [5] considered clap-andfling near the ground and found force enhancement at all distances. A three-dimensional numerical simulation of fruitfly hovering was carried out by Maeda and Liu [6]. An increase in lift and a reduction in power was found. A significant vertical force was generated on the insect's body due to the 'fountain effect'. Energetic savings have also been reported for a hummingbird hovering in ground effect [7].
Several studies considered pitching-plunging foils near a solid wall or a free surface [8][9][10]. This configuration is relevant to fish swimming as well as forward flapping flight. The ground effect mainly consists in enhanced propulsive force. However, it also generates a non-zero vertical force due to asymmetry.
The main motivation for this study comes from the fact that the ground proximity is natural for takeoff and landing. These manoeuvres, unlike hovering or forward flight, are characterized by gradual change of distance to the ground. The 'dynamic' ground effect in these circumstances may be different from the 'static' effect at a constant distance [11]. This difference may be even larger for flapping wings than for fixed wings, because animals vary their wing kinematics during takeoff.
So far, the ground effect during takeoff has been assessed for very few insects only. It was found negligible for butterflies (Pieris rapae [12], Papilio xuthus [13]), a dronefly (Eristarlis tenax) [14], and a fruitfly (Drosophila virilis) [15], but significant for a beetle (Trypoxylus dichotomus) [16]. The disparity can be attributed to significant differences in the size, morphology and kinematics of these insects. Thus, our work is motivated by the apparently contradictory conclusions on the significance of the ground effect that could be found in the animal flight literature. It is important to identify the parameters that make the ground effect strong or negligible.
In the present study, we consider a numerical model having the morphology of a fruitfly, with variable wing kinematics and leg parameters. Our objective is to determine if the ground effect can be significant for this model, and which conditions can lead to it. We thus explore the parameter space of the model and perform numerous numerical simulations using FLUSI [17], which is an open source software available on https://github.com/pseudospectators/ FLUSI/tree/plos_one_ground_effect. First, for completeness, we revisit the voluntary takeoff of a fruitfly analyzed in [15]. The main difference with respect to [15] is the use of a flight dynamics solver. We then compare takeoffs with modified parameters of the leg thrust model and wing kinematics. Finally, we consider hovering as a limiting case of very slow takeoff.
The paper is organized as follows. In section, we describe our computational approach and the takeoff parameters used in this study. The results are presented in section Leg thrust, first for a natural voluntary takeoff, then for modified takeoffs and for hovering flight. The main conclusions are summarized in section Conclusions.

Morphology and kinematics
In this work, we consider a fruitfly having mass m = 1.2 mg and wing length R = 2.83 mm, which are the values reported by Chen and Sun [15]. The body is modelled as a rigid solid, and the wings are modelled as rigid flat plates. This approximation is accurate for Drosophila during voluntary takeoff [15,18], though it occasionally fails during fast escape manoeuvres [18]. The wing contour used in this study is shown in Fig 1a. It is adapted from [15]. Its mean chord length is equal to c = 0.85 mm. The body is generated by sweeping a circular section of variable radius along a curvilinear centreline (an arc). The body has approximately the same dimensions as in [15]. The side view of the body is shown in Fig 1b. Even though the yaw and roll angles can eventually become large during takeoff, there is no significant trend for all takeoffs. Hence, to simplify the problem, we assume bilateral symmetry. Therefore, the body orientation is fully defined by the pitch angle β between the body and the horizontal axis, see Fig 1c. The wing kinematics is described by three angles: ϕ, α and θ, measured with respect to the stroke plane, as shown in Fig 1d. The positional angle ϕ defines the motion of the wing tip projection on the stroke plane. The deviation (elevation) angle θ defines the deviation of the wing tip from the stroke plane. The feathering angle α defines the rotation about the longitudinal axis of the wing, and it is related to the geometrical angle of attack (AoA) as α = 90°− AoA during downstroke and as α = 90°+ AoA during upstroke. It is convenient to refer to an 'anatomical' stroke plane angle η, i.e., to assume that the inclination of the stroke plane against the body axis is held at a constant angle for any motion of the body. Since the main focus of this study is the ground effect, it is important to ensure that the time evolution of the distance to the ground is consistent with the forces acting on the insect. For this reason, in our computations, unlike in [15], the position of the insect is dynamically computed as opposed to be prescribed. We compute the position of the body point of reference (x c , z c ), see Fig 1b, from Newton's 2nd Law, where (F ax , F az ) is the aerodynamic force, (F ℓx , F ℓz ) is the leg thrust, subscripts x and z correspond, respectively, to the horizontal and vertical components, m is the insect's mass and g is the gravitational acceleration. Eq (1) are integrated using the adaptive second order Adams-Bashforth scheme [19], simultaneously with the incompressible Navier-Stokes equations. We defined the positive z direction to be upwards and the positive x direction to be forwards (see Fig 1c).

Aerodynamics
The aerodynamic forces F ax and F az are obtained by solving the three-dimensional incompressible Navier-Stokes equations. The no-slip boundary condition at the body and wings surfaces is imposed using the volume penalization method [20], and the penalized equations are solved using a classical Fourier pseudo-spectral method. More details about the solver and the generic insect model, including a numerical validation case of fruitfly hovering, can be found in [17]. Numerical validation of the ground plane modelling using the volume penalization method is described in S1 Appendix. The computational domain in the present study is a rectangular box with sides L x , L y and L z . Suitable values of L x , L y and L z , in terms of accuracy and computational efficiency, depend on the motion of the insect within the domain. Therefore, different values are used in different simulations, as described later in the text. The domain is discretized using a uniform Cartesian grid. Periodic boundary conditions are applied on all sides of the domain, as required by the Fourier discretization. Vorticity sponge boundary conditions are imposed at the left, right, rear and front sides of the domain, as explained in [17,21], in order to minimize the effect of the finite domain size. The ground surface is modelled as a solid layer at the bottom of the domain, which by periodicity also imposes the no-slip on the top of the fluid domain. We have carried out numerical experiments to ensure that, in the numerical simulations presented in this paper, the domain size is sufficiently large, i.e., its further increase does not change significantly the forces. The dimensions that we chose are also comparable with the size of the mineral oil tank used in the experiments with a mechanical model [22].

Leg thrust
The model of the leg thrust employed in the present study is a slight modification of the compression spring model proposed in [12]. We assume that takeoff begins from rest and starts at time t = t ℓ , which can be estimated from the initiation of the legs motion in the video sequences shown in [15]. The two components of the force are given by The magnitude of the leg force is assumed to depend on the vertical component of the leg extension z = z c (t) − z c (t ℓ ) only. The force is supposed to be distributed between the three pairs of legs such that its change with horizontal displacement can be neglected, where L ℓ is the maximum leg extension length, i.e., the difference between the values of z c when the legs are fully extended at takeoff and when the insect is at rest. When the legs are fully extended, z = L ℓ , the legs lose contact with the ground and the force drops to zero. This length is estimated using video sequences in [15] to be equal to L ℓ = 1.24 mm. The spring stiffness K ℓ varies in time: it increases from K À ' before takeoff to K þ ' after takeoff. The initial value K À ' ¼ mg=L ' ensures that the insect is in equilibrium before takeoff, when the aerodynamic force is zero. The final value K þ ' is a parameter of the model that controls the maximum leg thrust. Its value can be estimated from the climb velocity at the beginning of takeoff, shown in, e.g., [15]. It may also be estimated from jumps of wingless flies [23,24] for a slightly different fruitfly, D. melanogaster. We assume the time evolution of K ℓ of the form The transition time τ ℓ can be equal to zero, in which case the leg force increases impulsively at the beginning of takeoff. However, measurements of the leg force [23] suggest a gradual increase which can be accounted for by setting τ ℓ to a value larger than zero. The value τ ℓ = 1.3ms results in the gradient dF ℓ /dt consistent with the experimental data shown in [23]. The direction ϕ ℓ also changes in time. Before takeoff, when the insect is at rest, the force is applied only in the vertical direction, i.e., À ' ¼ 90 . During takeoff, the horizontal component is non-zero, in general. We assume a time evolution of the form The values of the leg thrust model parameters used in our numerical simulations are given in Table 1.

Summary of the numerical simulations
The starting point for our study is the voluntary takeoff, as it is shown in section Voluntary takeoff (in agreement with [15]) that the ground effect is very small in that case. It is much smaller than during hovering (cf. [6]). We conjecture that this difference is due to the large takeoff vertical velocity, which is mainly the result of the leg thrust. To test this hypothesis, in section Slow takeoff, we discuss a situation in which the legs produce less force and the insect takes off slower. The ground effect becomes significant. The vertical force increases during the first two wingbeats due to the ground effect, but slightly decreases later on. We then carry out a parametric study using periodic wing kinematics in section Takeoffs with simplified kinematics, and find an even stronger adverse ground effect. Finally, in section Ground effect in hovering flight, we find similar trends during the first wingbeats in hovering flight, which can is considered as a limiting case of takeoff with zero rate of climb. Tables 1 and 2 summarize the parameters of the different cases considered in the present study. Datasets for the 'voluntary' and 'simplified' cases can be downloaded from [25].

Voluntary takeoff
In this section, we consider voluntary takeoff of a fruitfly with the parameters as in the first lines in Tables 1 and 2. This case shows some important general features of fruitfly takeoff such as the first wingbeat cycles beginning while the legs extend. Therefore it is likely that, despite some variability in voluntary takeoffs, the ground effect in general remains of the same order of magnitude in natural circumstances. The values of the body and wing angles are taken from one of the cases documented in [15]. However, the wing motion in [15] is not exactly symmetric. Therefore, the time series of ϕ, α and θ that we use for both wings correspond to the left wing data shown in [15]. Fig 2 presents the time evolution of the wing positional angle ϕ(t), the feathering angle α(t), the elevation angle θ(t) and the body pitch angle β(t), which are prescribed in our numerical simulations. The angle between the horizontal plane and the stroke plane, η − β, is also shown for reference.
Even though the wing motion is not exactly periodic, it is useful to introduce the wing beat frequency. When calculated using the average wing beat cycle period over the five cycles shown in Fig 2, it is equal to f = 169 Hz. Similarly, the average wing beat amplitude is equal to F = 134°, and the characteristic wing tip velocity is U = 2FRf = 2.23 m/s. The kinematic viscosity of air, equal to ν = 1.45 Á 10 −5 m 2 /s yields the Reynolds number Re = Uc/ν = 131. Note that U and Re do not account for the forward speed of the body.
The computational domain size is equal to L x = L y = 5R, L z = 8R, where R is the wing length. The influence of the domain size in the vertical direction is discussed in S2 Appendix. The number of grid points in each direction, respectively, is N x = N y = 640 and N z = 1280. The penalization parameter is ε = 2.5 Á 10 −4 (for details see, e.g., [21]). The aerodynamic ground effect is evaluated by comparing two numerical simulations with two different values of the initial distance from the body point of reference to the ground: z c (0) = 0.38R and 2R, which we denote 'in ground effect' (IGE) and 'out of ground effect' (OGE), respectively. The first case corresponds to a takeoff from a flat ground surface, with z c (0) being consistent with the data in [15]. In the second case, the leg model behaves as during takeoff from the ground, but the aerodynamic interaction between the insect and the ground is weak because of the large distance. This case may be interpreted as takeoff from a perch that provides enough support for the legs but has a small surface, such that the aerodynamic interactions are negligible. With the distance equal to 2R or greater, the ground effect is negligible during hovering [13]. The circulation of the wake vortices is mainly determined by the integral aerodynamic force, therefore it is not larger during takeoff than during hovering, and the spatial rate of decay of the induced velocity is the same. Hence, the ground effect with the distance equal to 2R is likely to be negligible during takeoff. The influence of the ground on the shape of the vortices is only visible during the 2nd wingbeat and later on. This difference is localized to the vicinity of the ground plane. Since the insect is relatively far from the ground by that time, this difference is unlikely to have any influence on the aerodynamic forces. Fig 3a shows the fruitfly model and the wake, IGE and OGE, at 4 subsequent time instants. The vortices created by the wings and the body are identified as the volume of fluid enclosed by the iso-surfaces of the Q-criterion. At t = 0, the air is at rest. The insect body is almost horizontal. The wings are in a pre-takeoff position from which they begin the first downstroke after t = 4.1 ms. The time t = 9.2 ms corresponds to the first reversal from downstroke to upstroke. Because of the small body pitch angle β, the stroke plane is effectively vertical. In addition, the wing tip speed during the first downstroke is smaller than during all subsequent strokes. Therefore, the vertical aerodynamic force is small, but the body lifts noticeably because of the leg thrust. The time t = 12.8 ms corresponds to the second upstroke. At this point, the distance from the body point of reference to the ground z c is already larger than the wing length R. Therefore, the aerodynamic interference with the ground is expected to be very small. Note that the kinematics during the first two wingbeat cycles are a transient. After that, the time evolution of the wing angles approaches a periodic regime and the stroke plane becomes less inclined with respect to the ground, see   The displacement of the body point of reference is shown in Fig 3b. It presents the evolution of the vertical component z(t) = z c (t) − z c (0) and the horizontal component ξ(t) = x c (t) − x c (0) over time for the cases IGE and OGE. The curves overlap. In both cases, at the end of the 4th wingbeat cycle, t = 28.4 ms, the insect gains 8.7 mm of altitude and propels 1.6 mm forward. These numbers are consistent with the trajectories shown in [15]. The displacements IGE and OGE differ by less than 1%. Therefore, the ground effect on z and ξ is indeed negligible. Fig 3c shows the two components of the leg force. At t = 0, the vertical component of the leg force is equal to the weight and the horizontal component is zero. The jump is triggered at t ℓ = 4.2 ms. At time t ℓ + τ ℓ = 5.5 ms, both components reach their peaks. After that the force decreases and vanishes at t = 9.3 ms, when the legs lose contact with the ground. Note that the leg thrust can, in principle, be different for the takeoffs IGE and OGE, because the leg model depends on the aerodynamic force via z c (t). However, for the voluntary takeoff considered here, there is no influence of the ground effect.
The horizontal and the vertical components of the aerodynamic force are shown in Fig 3d  and 3e, respectively. Over the first four wingbeat cycles, the wingbeat averaged aerodynamic forces are significantly lower than the weight. This can be explained by the large initial rate of climb due to the leg thrust, which cannot be supported by the wings. Even during the fourth wingbeat, the wing force is equal to 29% of the weight. The vertical acceleration is therefore negative after the legs lose contact with the ground, and the rate of climb slowly decreases. The ground effect is, again, negligible. Even during the first wingbeat cycle, when the wings approach the ground surface, the difference in the instantaneous vertical force between IGE and OGE is at most 0.0005 mN, i.e., about 4% of the weight. The wingbeat cycle averaged forces differ by less than 1% of the weight. Fig 3f displays the time evolution of the aerodynamic power, when operating IGE and OGE. Note that, in this study, we do not consider the inertial power because the wings have the same kinematics in both cases, IGE and OGE. Therefore, the inertial power is the same. The aerodynamic power is the aerodynamic component of the power required to actuate the wings, In Eq (6), M l and M r are the aerodynamic moments of the left and of the right wing, respectively, relative to the corresponding pivot point. O l and O r are the angular velocities of the wings and O b is the angular velocity of the body. All vectors are taken in the laboratory frame of reference. P is positive if power is consumed. We find that it is positive during most part of the takeoff (see Fig 3f). Only at the reversals during the first two cycles, when the body velocity is still small, P is slightly negative. During the 2nd wingbeat, the mean body-mass specific aerodynamic power is equal to P Ã b ¼ P ave =m ¼ 21 W/kg. Assuming that the muscles contribute to 30% of the body mass, the mean muscle-mass specific aerodynamic power is equal to P Ã m ¼ P ave =ð0:3mÞ ¼ 69 W/kg. The relative difference in the cycle averaged values between IGE and OGE is less than 0.5%.
We conclude that the ground effect is unimportant for the voluntary takeoff, a result which is in agreement with [15]. This is mainly a consequence of rapid acceleration during the first wingbeat cycle, when the legs produce a large vertical force. The main question of the next section is whether this scenario changes if the takeoff is slower and the insect remains near the ground for a longer time. The rate of climb at the beginning of takeoff is controlled by the leg model stiffness coefficient K þ ' , and the horizontal velocity is controlled by the leg angle þ ' .

Slow takeoff
This section describes a modified takeoff with the leg thrust coefficient decreased to K þ ' ¼ 0:041 N/m (see the second line in Table 1). Smaller K þ ' results in less leg thrust and slower climb, compared to the natural voluntary takeoff. Therefore we refer to this case as a 'slow takeoff'. Any further decrease of K þ ' would result in a very close approach of the insect to the ground surface and, ultimately, to a collision which would require special treatment. Modelling of such collisions could be an interesting topic for future research (see [27] for a review of structural modelling of insect wings, including impact modelling). For a fruitfly, collisions between the wings and the ground may be undesirable because of the large wingbeat frequency and light wing structure. Therefore, K þ l ¼ 0:041 N/m is an interesting limiting case. In the present numerical simulations, the computational domain size in x and y directions, the discretization grid step size and the penalization parameter are the same as in the previous section. The domain size in the vertical direction z is reduced to 6R because the insect gains much less altitude by the end of the simulation. Fig 4a shows the displacement of the body point of reference. The rate of climb is about one third of its original value and the insect only gains 3.9 mm by the end of the 4th wingbeat cycle. This is just slightly larger than the wing length R (2.83 mm). The displacement is slightly larger for IGE than for OGE in both directions, horizontal and vertical. The time evolution of leg thrust is given in Fig 4b. There is no visible difference between the two cases. The peak of the vertical force is equal to 0.051 mN, which is about four times less than in the original voluntary takeoff discussed in section Voluntary takeoff.
The time evolution of the instantaneous aerodynamic force, shown in Fig 4c and 4d, is qualitatively similar to the voluntary takeoff case considered previously. The difference between the cases OGE and IGE is negligible for the horizontal force (Fig 4c), but for the vertical force it reaches values as large as 0.0027 mN, i.e., 23% of the weight (Fig 4d). Fig 5 shows the difference between the wingbeat averaged forces in the cases IGE and OGE, normalized by the weight. The vertical force difference is shown in Fig 5b. During the 1st wingbeat, the ground effect makes the total vertical force increase by almost 6% of the weight (red line). However, during the 2nd wingbeat, the extra force decreases to only 2% of the weight. During the 3rd, the 4th and the 5th wingbeats, the difference between the vertical forces IGE and OGE is very small and negative. The increase of the vertical force during the first wingbeat is mainly due to the wings (blue line). The extra force acting on the body is only about 1% of the weight. However, later on, the contribution of the body becomes important because it remains positive, whereas for the wings it becomes negative. The horizontal force difference, shown in Fig 5a, is positive, i.e., the propulsive force increases due to the ground effect by about 2% of the weight, for all wingbeats. The contribution of the body is up to 1% of the weight. For reference, Fig 5 also shows the force differences during the voluntary takeoff. They are all smaller than 1%.
The aerodynamic power, in the cases IGE and OGE, is compared in Fig 6. The maximum difference is of about 3% in magnitude for the slow takeoff, but less than 1% for the voluntary takeoff. Considering the slow takeoff, the insect consumes more power when operating in ground effect (IGE) during the first two wingbeat cycles, but less power during the subsequent cycles. Overall, we find that the differences in the power are small.

Takeoffs with simplified kinematics
In the previous sections we noticed that the ground effect depends on the takeoff kinematics. We are mainly interested in the effects that might be generally applicable to a fruitfly sized insect takeoff. Therefore, in this section we consider parametric studies. They are performed using simplified periodic wing kinematics. The time evolution of the three wing angles over one wingbeat period, obtained by periodization of the last wingbeat in [15], is shown in Fig 7a. The wingbeat frequency is equal to f = 210 Hz. This takeoff mode can be relevant to MAVs, for which the wing kinematics and the body angle during takeoff do not change as much as for the fruitfly (see, e.g., [28]). The leg strength parameter K þ ' is varied, resulting in a variation of the takeoff rate of climb V t.o. . The body angle is constant and equal to β = 46.3°, the anatomical stroke plane angle is equal to η = 32°. In these computations, we use L x = L y = 4R, L z = 6R, N x = N y = 512 and N z = 768, corresponding to more than 200 million grid points. The penalization parameter is equal to ε = 2.5 Á 10 −4 .
Smaller K þ l implies smaller rate of climb (Fig 8c) which leads to a more significant ground effect (Fig 8d). A striking feature of Fig 8b is a significant decrease of the vertical force during the 4th, 5th and 6th wingbeats, by up to 6%. The horizontal force varies slightly, by about 1% (see Fig 8a).
The decrease in vertical force on the wings found in the IGE cases compared to OGE (Fig 8b) is an adverse ground effect which may be the result of complex wing-wake-ground interactions that depend on wing and body kinematics. Adverse ground effects have been reported previously for fixed-wing aircraft [11], but they consist in increased drag together with increased lift. However, for a two-dimensional ellipse with normal hovering kinematics [4], the mean vertical force decreases when the height from the ground is between approximately 1.5D and 4D, where D is the chord length of the ellipse. For flapping wings, an adverse ground effect was found by Quinn et al. [9]. They considered an airfoil undergoing pitch oscillations in a closed-loop water channel with prescribed free-stream velocity. Such a configuration represents a section of a bird wing in forward flight or a fish fin. Even if the pitching motion was symmetric, the proximity of the ground broke the symmetry of the flow. Thus, the airfoil produced non-zero lift. The lift was positive if the distance to the ground was less than 40% of the chord length, but it became negative at larger distances, such that the lift force pulled the airfoil towards the ground. Nevertheless, the extra propulsive force due to the ground effect was positive in all cases. Note, however, that the flows considered in [4] and [9] are effectively two-dimensional.  [15]. It is used in section Takeoffs with simplified kinematics. Also, in section Ground effect in hovering flight, this kinematics is referred to as 'P1' ('P' for 'Periodic'). The cycle begins from the upstroke. (b) Data adapted from [6]. In section Ground effect in hovering flight, it is referred to as 'P2'. The cycle begins from the downstroke. In the present work, importantly, we find an adverse ground effect in a three-dimensional configuration, which has not been previously recognized. The wingbeat cycle averaged vertical force of the wings in the ground effect is slightly larger during the first two cycles, but, as the insect flies away from the ground, the vertical force in the case IGE becomes less than that in the case OGE. Fig 9 shows a comparison of the wake at the end of the 5th wingbeat cycle (t = 23.78 ms) in the 4 different takeoffs IGE with different values of K þ l . It corresponds to the maximum decrease of the vertical force. There are noticeable differences between the vortices when the takeoff is fast and when it is slow. The part of the wake that approaches the ground deforms when it impinges on the ground. It then rolls up in a pair of vortex rings. Similar "ground vortices" are known in the context of helicopter rotor aerodynamics. In each of the 4 cases shown in Fig 9, they have different strength and position with respect to the wings. Therefore, they induce the downwash of different strength.
A detailed view of the flow near the wings is presented in Fig 10. It shows the pressure and the vorticity magnitude during the 6th wingbeat for two IGE cases with different K þ l . Five time instants are visualized in five rows, respectively. The left column shows the pressure distribution over the surface of the insect, as well as over a semisphere of radius 0.9R centred the body reference point, for K þ l ¼ 0:0095 N/m. Fig 10a is at t = 24.73 ms, during upstroke. The dark blue area near the leading edge that expands towards the wing tip is the trace of the leadingedge vortex (LEV), similar to the one discussed in [29]. The LEV at the downstroke is evident in Fig 10d at t = 27.58 ms. The pressure distributions during the reversals are more complex (Fig 10b, 10c and 10e). The middle column of Fig 10 shows the pressure iso-contours sampled on the surface of the semisphere. Two cases are compared, K þ l ¼ 0:0095 N/m and 0.0430 N/m. This choice of K þ l corresponds to the largest difference in the vertical force (see Fig 8b). For each of the cases, two isolines are drawn, p/ρR 2 f 2 = −2 and −5. There is a large difference between the contours for different K þ l in the far wake (wake far from the wings). However, in the near wake of a wing including the LEV, the difference is much smaller, which is consistent with the force changing by only a few per cent.
The iso-contours of the vorticity magnitude are compared in the right column. Here again, despite significant differences in the far wake, the LEV contours virtually overlap for K þ l ¼ 0:0095 N/m and 0.0430 N/m. This shows that the ground effect has almost no influence on the dynamics of the LEV.

Ground effect in hovering flight
In this section, we simplify the kinematics even further. We consider hovering with the insect body being fixed. The flight dynamics solver is not used in this case. The distance from the body centre to the ground is equal to 0.48R (where R is the wing length) for hovering in ground effect (IGE) and 2.4R for hovering out of ground effect (OGE). The body pitch angle and the anatomical stroke plane angle are both constant and equal to 55°, such that the stroke plane is horizontal. The wing kinematics is the same as in the previous section, see Fig 7a. We denote it  Aerodynamic Ground Effect in Fruitfly Sized Insect Takeoff as 'P1' kinematics. The wingbeat frequency is equal to f = 218 Hz. The first wingbeat starts from the upstroke, as done in [15].
In these numerical simulations we are interested in the long-time evolution of the aerodynamic forces, which after the initial transient eventually reach a periodic state. Most of the results known from the helicopter rotor theory [2] are obtained in reference to the periodic state, while the takeoffs considered in the previous sections of this paper (the slow takeoffs, in particular) last only for a few wingbeats. Therefore, it is instructive to consider the time evolution of the aerodynamic forces during hovering from t = 0 until the time when the periodic state is reached.
Since the time span of the numerical simulations presented in this section is large, it is necessary to increase the domain size in the horizontal directions. We set L x × L y × L z = 8R × 8R × 4R, where z is the vertical direction. The number of grid points is N x × N y × N z = 864 × 864 × 432. The penalization parameter is equal to ε = 2.5 Á 10 −4 .
The quantity of interest is the ratio of the wingbeat averaged forces, IGE to OGE: F z,ave,IGE / F z,ave,OGE . This quantity is shown in Fig 11a. The red solid line with "+" symbols corresponds to the total vertical force ratio. As already noticed in [6], the vertical force during hovering in ground effect, F z,IGE , reaches its periodic state significantly later than during hovering out of ground effect, F z,OGE . Therefore, the ratio of their wingbeat averages, F z,ave,IGE /F z,ave,OGE converges slowly with the number of wingbeats. It oscillates between 102% and 107%. After 27 wingbeats it reaches 106%.
A similar comparison for the force generated by the wings is shown with a red dot-dashed line in Fig 11a. It was calculated by integration of the distributed forces over the wings only, in the same numerical simulations. Therefore, the aerodynamic interaction between the body and the wings is included. This force ratio drops from 100.6% to 92.5% during the first 6 wingbeats, oscillates and then increases to 96.3%.
The time evolution of the wake vortices generated by the insect is shown in Fig 12. The visualized time instants correspond to the end of the 1st, 2nd, . . ., 25th downstroke. There are significant differences between the positions of the vortices over the first four time instants. The first wingbeat generates very strong vortex rings, that collide with the ground. Then they rebound during the second wingbeat, and parts of them moving upwards are still visible during the third wingbeat. The downwash produced by these vortices influences the nearer wake dynamics and it is likely to be responsible for the decrease of the vertical force during the first few wingbeats. After the 10th wingbeat, the wake approaches its quasi-periodic state. There are almost no visible differences between the visualizations at the end of the 20th wingbeat and at the end of the 25th wingbeat.
The pair of numerical simulations (cases IGE and OGE) that we have discussed in the above paragraphs leads to the following conclusions.
1. Over the first 27 wingbeats, F z,ave,IGE /F z,ave,OGE varies within about 5% for the total force and 9% for the wings force.
2. The wing generates less vertical force in the case IGE than in the case OGE (adverse ground effect).
3. The body makes an important contribution to the total vertical force when operating IGE, which results in the excess total vertical force (positive ground effect).
These conclusions are, of course, only valid for the particular wing shape and kinematics used in the simulation. Periodic flapping is only an approximation to the real insect wing motion, which varies from one wingbeat to another, and depends on many different conditions. To determine the effect of all existing fruitfly wing kinematics is beyond the reach of our numerical simulations. However, it is useful to compare a few different cases.
We carried out numerical simulations with the wing kinematics used in [6] (abbreviated as 'P2' in the figures). Note that, in this case, the first wingbeat begins from the downstroke, as shown in Fig 7b. The wingbeat frequency is the same, f = 218 Hz. The results of these numerical simulations are shown in Fig 11 with green lines. They are qualitatively similar to the previously shown 'P1' case, but the values are systematically larger. F z,ave,IGE /F z,ave,OGE reaches 111% for the total force and 100% for the wings, such that there is no adverse ground effect after the periodic state is established.
The adverse ground effect is rarely encountered in the aircraft or rotorcraft aerodynamics literature. However, in the context of flapping wings, it is not unusual. In the two-dimensional numerical simulations [4], a U-shape profile of the force ratio F z,ave,IGE /F z,ave,OGE versus h/c was found, where h is the distance from the wing centre to the ground and c is the wing chord. F z,ave,IGE /F z,ave,OGE was greater than 100% for h/c < 1.5, but less than 100% for h/c > 1.5, and the minimum ratio was of about 54%.
In the three-dimensional model considered in the present paper, it is possible to partially reduce the three-dimensional effects by decreasing the wing beat amplitude. One can then expect the adverse ground effect to be amplified. Indeed, this is what we find by decreasing the wingbeat amplitude by a factor of 2 (by rescaling the positional angle shown in Fig 7 such that the amplitude of the positional angle is halved but the mean positional angle is unchanged, and the wingbeat frequency remains unchanged). The results are shown in Fig 11 with blue lines. Now we have F z,ave,IGE /F z,ave,OGE < 100% at intermediate times, for wings force and for the total force. The final periodic state produces a very slight excess of the total vertical force (less than 1%). Note that, in this reduced-amplitude case, the total force ratio decreases more than the wing force ratio. This indicates, not surprisingly, that the fountain effect becomes weaker when the wingbeat amplitude is reduced. Similar computations with the 'P2' kinematics show the same trend with an even larger decrease of F z,ave,IGE /F z,ave,OGE .
The wingbeat averaged aerodynamic power ratio P ave,IGE /P ave,OGE is shown in Fig 11b. Its variation is smaller than the variation of the force, and the computations suggest that its longtime limit is between 97% and 99%, in all cases that we have considered. The shape of the time evolution profiles of the power ratio is approximately similar to the time profiles of the wings vertical force ratio. This means that a local decrease of the vertical force ratio is accompanied by a decrease of the power ratio. Therefore, if the kinematics of the wings operating in ground effect is adjusted such that P ave,IGE /P ave,OGE = 100% at any time, the force ratio F z,ave,IGE /F z,ave, OGE is likely to increase. Among other factors, the feathering angle is very likely to change passively, when in ground effect, due to compliance of the wing [30][31][32]. Such effects would need further investigation.

Conclusions
The aerodynamic ground effect in fruitfly sized insect takeoff has been studied numerically using high performance computing. The three-dimensional incompressible Navier-Stokes equations were solved using a pseudo-spectral method with volume penalization [33,34] using the FLUSI open source code [17], in order to obtain the flow field and the aerodynamic forces acting on the insect. The takeoff trajectories were calculated using a simple flight dynamics solver that accounts for the body weight, inertia, and the legs thrust. A series of computations has been carried out to explore the parametric space of the model. A natural voluntary takeoff of a fruitfly, modified takeoffs with different kinematics and leg model parameters, and hovering flights have been compared.
We found that the ground effect during the natural voluntary takeoff is negligible. The wingbeat averaged forces only differ by less than 1% of the weight. The aerodynamic power differs by less than 0.5%.
In the modified takeoffs, we decreased the leg strength. As a consequence, the rate of climb decreased and the ground effect became significant. Surprisingly, the vertical force did not always increase. It even dropped in some of the cases that we considered. This is an unsteady effect related to the vortex rings bouncing off the ground surface.
To better understand the mechanism of the adverse ground effect, we considered hovering near a flat ground surface, being the limiting case of zero rate of climb. In that case, the fountain effect produced a large upward force on the insect's body. The net ground effect was therefore positive. However, the aerodynamic force acting on the wings in ground effect was sometimes less than when the wings operate out of ground effect. The most significant decrease was observed during the first 15 wingbeats. Note that this is a much longer time period than a typical takeoff. At long time hovering, the effect was either positive or negative, depending on the wings kinematics.
The parameter space in the takeoff problem is very large. In the present study, we focused on the legs thrust and wing kinematics. However, the aerodynamic ground effect may also be sensitive to the Reynolds number, because the structure of the wake at high Re is significantly different from that at low Re. Since high Reynolds number computations are costly, they are beyond the scope of the present study, but it is an important question for future research.
Supporting Information S1 Appendix. Numerical validation of the ground plane modelling using the volume penalization method. See file S1_Appendix.pdf. (PDF) S2 Appendix. Influence of the domain size in the vertical direction. See file S2_Appendix. pdf. (PDF)