Passive mechanism of pitch recoil in flapping insect wings

The high torsional flexibility of insect wings allows for elastic recoil after the rotation of the wing during stroke reversal. However, the underlying mechanism of this recoil remains unclear because of the dynamic process of transitioning from the wing rotation during stroke reversal to the maintenance of a high angle of attack during the middle of each half-stroke, when the inertial, elastic, and aerodynamic effects all have a significant impact. Therefore, the interaction between the flapping wing and the surrounding air was directly simulated by simultaneously solving the incompressible Navier–Stokes equations, the equation of motion for an elastic body, and the fluid-structure interface conditions using the three-dimensional finite element method. This direct numerical simulation controlling the aerodynamic effect revealed that the recoil is the residual of the free pitch vibration induced by the flapping acceleration during stroke reversal in the transient response very close to critical damping due to the dynamic pressure resistance of the surrounding air. This understanding will enable the control of the leading-edge vortex and lift generation, the reduction of the work performed by flapping wings, and the interpretation of the underlying necessity for the kinematic characteristics of the flapping motion.


Introduction
Approximately 300 million years ago, insects became the first organisms on Earth to develop flight, an ability to which they owe much of their evolutionary success [1]. Insect flight has attracted the attention of physicists and engineers as well as biologists for more than a century because it involves many challenging subjects of great complexity. For example, mimicking insect flight is one approach being actively pursued for application in micro air vehicles [2]. The insect flight system can be decomposed into subsystems with different length and time scales, including inner locomotion, the surrounding air, and control, from anatomical, physiological, ecological, energetic, and mechanical perspectives [1]. These subsystems cooperate to create sophisticated insect behavior in flight.
The aerodynamic mechanisms of insect flight under well-characterized wing kinematics have been studied extensively [3,4,9,[50][51][52][53]. Figure 1(a) shows the schematic of wing motions. The basic pattern of the flapping motion is close to sinusoidal but with a larger acceleration or deceleration during stroke reversals and a constant velocity during the middle of each half-stroke [5][6][7]. Therefore, the time variant of the angular displacement to describe the flapping motion (flapping angular displacement) can be modeled as shown in figure 2. The basic pattern of the pitch motion consists of a wing rotation during the stroke reversal, a subsequent recoil, and the maintenance of a high angle of attack during the middle of each halfstroke [5,[7][8][9]. Figure 1(b) shows the schematic of the time variation of the wing pitch, which is drawn based on [5,9]. As shown in this figure, the wings rotate through a large angle at the end of each half-stroke during 10%-20% of the flapping period. The rotation carries the wing to a pitch angle some 10°more than that during the middle of each half-stroke. The wing then recoils to the proper value. However, the mechanisms that determine the wing motions remain unclear because they work across different subsystems.
The morphological structure of the wing base suggests the passivity of the pitch motion, i.e., the high torsional flexibility of the wing base prevents insects from transmitting the torsional moment actively applied by their own muscle to the outer wing [10,11]. Initially, a mechanical evidence was presented demonstrating the inertial cause of the wing rotation during stroke reversal using isolated subsystems [7], such as the aerodynamics under the boundary condition of an active pitch motion [12]. Recently, it was shown that the basic pattern of the pitch motion can be simulated without the aid of the active pitch motion in a fluidstructure interaction (FSI) system [13][14][15][16]. A high angle of attack can be maintained during the middle of each half-stroke by ensuring the wing's elastic and aerodynamic moments remain in equilibrium [16]. The generation of sufficient lift suggests that the pitch motion throughout the entire stroke is based on the passive mechanisms of inertia and aeroelasticity [15].
The recoil might be due to the wing's torsional elasticity but no proof was offered as described in [5]. However, the elastic effect alone is not likely to explain the recoil, since the natural frequency of the wing's torsion is a few times larger than the flapping frequency [14,22,47,48]. Most recently, a two-dimensional CFD analysis using the prescribed recoil showed that the recoil can improve the lift-to-drag ratio [49], and the image analysis revealed the time variation of the recoil in detail [9]. However, it remains unclear how the recoil is created. In this study, we revealed a c) the model wing using lumped torsional flexibility as a first approximation. Φ is the stroke angle; θ is the pitch angular displacement; j is the flapping angular displacement; N denotes the Nth cycle; n and n * are the unit vector normal to the wing's upper surface and its projection onto the stroke plane, respectively; F S and F S * are the fluid surface force acting on the wing and its projection onto the stroke plane, respectively; F L is the lift force defined as the projection of F S onto the y-axis; and F D is the drag force defined as (F S * · n * ).
passive mechanism of the recoil based on the vibration theory. As shown in figure 1(b), the recoil is the dynamic process of transitioning from the wing rotation during stroke reversal to the maintenance of a high angle of attack during the middle of each half-stroke. The inertial, elastic, and aerodynamic effects all play a significant role in this process. Because of this complexity of the recoil, its mechanism remains unclear as long as we know. Therefore, a complete FSI system was numerically analyzed for the purpose of performing a direct simulation of this process. We revealed from the present direct simulation that the recoil is the residual of the free pitch vibration caused by the flapping acceleration in the transient response very close to critical damping resulting from the dynamic pressure resistance of the surrounding air.
On the basis of this understanding of the recoil, we found the recoil effect on the lift as follows: the recoil enhances the early development of the leading edge vortex, and increases the mean lift during the early period of each half stroke. Moreover, we found that the kinematic characteristic of the flapping motion contributes to the energy efficiency without sacrificing the hovering performance due to the recoil effect invoked by them. Note that the kinematic characteristic of the flapping motion is close to sinusoidal but with a larger acceleration or deceleration during stroke reversals and a constant velocity during the middle of each half-stroke. This result of the energy efficiency would explain the underlying necessity of the kinematic characteristic of the flapping motion observed in actual insect wings.

Insect wing model
A high torsional flexibility concentrated at the base of an insect wing [10,11], which is schematically shown in figure 1(a), leads to a high angle of attack that separates the air flow around the wing and creates a leading-edge vortex that provides sufficient lift for the insect to maintain flight. However, other deformations to the wing, such as twist and camber, are relatively Time histories of (a) the stroke angular displacement j, (b) the stroke angular velocity dj/dt, and (c) the stroke angular acceleration d 2 j/dt 2 for flapping motion modeling. Φ is the stroke angle, T j is the flapping period, and t a is the acceleration time. The red lines indicate the histories for t a =T j /4, the blue lines indicate the histories for t a =T j /2.5, and the gray lines indicate the sinusoidal histories.
small [9,17,18], and flow separation is not very sensitive to such deformations [19]. The morphological structure of the wing base contributing to the pitch motion can be reduced to the macroscopic constitutive relationship between the torsional moment applied to the wing and the resultant pitch angular displacement from the perspective of the mechanics of materials. The constitutive parameters can be dynamically controlled to consider maneuvering flight [20]. A static constitutive equation, which was obtained from a torsional test of the wing [11], was used to describe hovering flight. Therefore, as a first approximation, an insect wing can be modeled using lumped torsional flexibility [14], which consists of a rigid flat-plate wing and a torsional spring, as shown in figure 1(c). This type of model [13,15,16,[20][21][22][23][24][25] and the actively controlled rigid flat-plate wing model [3,4,7,12,[26][27][28][29][30][31][32][33] have recently gained popularity.

Model parameters and governing equations
Assuming an ideal hovering state, the model wing using the lumped torsional flexibility in figure 1(c) is defined as follows. The wing base is placed at the origin of a Cartesian coordinate system. The stroke plane is set in the horizontal xz plane because it is almost horizontal in many hovering insects [34]. The flapping axis is defined as the vertical y-axis. The angular displacement of the flapping motion j (−Φ/ 2jΦ/2, where Φ is the stroke plane angle) is defined as the angle between the leading edge of the wing and the z-axis and is positive for counterclockwise rotation about the y-axis. The time variation of j is close to sinusoidal but has a larger acceleration or deceleration during the stroke reversals and a constant velocity during the middle of each half-stroke [5][6][7]. Therefore, dj/dt can be expressed as a trapezoidal function, as shown in figure 2(b), where the acceleration time t a is used to define the period of increasing or decreasing the flapping velocity in each cycle. The flapping period T j is the inverse of the flapping frequency f j . Note that the change in the value of dj/ dt averaged over each half-stroke is less than 1% when the acceleration time t a is varied. The initial value of j is set to −Φ/2. The pitch angular displacement θ is defined as the angle between the wing chord and the vertical direction and is positive for counterclockwise rotation about the wing span axis. The initial state of the wing chord is set to be at rest and oriented vertically (normal to the stroke plane).
The nondimensional lift F L is defined as the ycomponent of F S , and the nondimensional drag F D is defined as the inner product of F S * and n * (figure 1(c)), where F S is the total fluid surface force acting on the wing nondimensionalized by the dynamic pressure term ρ f A w V m 2 /2 (ρ f is the fluid density, A w is the wing area, and V m is the mean flapping velocity), F S * is the projection of F S onto the stroke plane, n is the unit vector normal to the upper surface of the wing, and n * is the projection of n onto the stroke plane. The mean lift force is dependent on the nondimensional radius r 2 of the second moment of the wing area [15,17,35]. r 2 of many insects [17,18] is approximately equivalent to that of a rectangular wing at 1/√3. Therefore, a rectangular wing model is used for the sake of simplicity. Its planar shape can be defined by the aspect ratio r A =2L w /c m , where L w is the span length of one wing and c m is the mean chord length. In this study, r A is set to 11 based on the morphological data of an actual crane fly [17]. The torsional flexibility of the model wing is defined by the torsional compliance C w [11], which is the torsional angular displacement of the wing divided by the applied torsional moment.
The behavior of the model wing can be fully described by the governing system of equations for FSI [14], which consists of the three-dimensional partial differential equations of motion for the elastic body, the three-dimensional incompressible Navier-Stokes (NS) equations for the fluid, and the compatibility and equilibrium equations for the fluid-structure interface as follows: The equilibrium equation for the elastic body to describe the motion of the deformable model wing can be expressed as where d/dt in the left-hand side is the Lagrangian time derivative, the superscript s indicates a quantity corresponding to the structure, ρ s is the mass density of the structure, v i s is the ith component of the structural velocity vector, and ij s s is the ijth component of the Cauchy stress tensor of the structure. Note that i equals 1, 2, or 3, and corresponds to the x, y, and z direction, respectively, in the given Cartesian coordinate system.
The incompressible NS equations to describe the fluid flow surrounding the wing using the arbitrary Lagrangian-Eulerian (ALE) method [36] can be expressed as The interface conditions to describe the interaction between the wing and the surrounding fluid can be expressed as where n i f and n i s are the ith components of the outward unit normal vectors on the fluid-structure interface corresponding to the fluid and the structure, respectively.

Dimensionless numbers
The four independent dimensionless numbers considered in this study (the Strouhal number St [38], the Reynolds number Re [38], the Cauchy number Ca [39], and the mass ratio M [40]) are obtained by the dimensional analysis for the governing system of equations describing FSI (1)-(3) and can be used to measure the dynamic similarity between the model wing and an actual insect wing as follows: These four dimensionless numbers can be defined with respect to the characteristic length c m ; the characteristic velocity, defined as the mean flapping velocity V m ; and the characteristic time where ρ f is the fluid mass density, μ f is the fluid viscosity, and m w is the wing mass. The definition of St is equivalent to the geometric similarity condition because V m =2r 2 ΦL w f j reduces St to 1/(r 2 r A Φ).
Because the nondimensional radius r 2 and the aspect ratio r A for the model wing are approximately equal to those for an actual insect wing, the stroke plane angle Φ in the model wing is set to 123°based on kinematic data for an actual crane fly [5]. Substituting the morphological [17], kinematic [5], and material [11] data for an actual crane fly into equations (4a)-(4d), St=0.073, Re=254, Ca=0.19, and M=16 are obtained as the dimensionless numbers describing the case that is dynamically similar to an actual wing. These values are hereafter referred to as the reference values.

Monolithic equation system
Applying finite element discretization to the total Lagrangian formulation [37] of the equation (1), the nonlinear equilibrium equation system can be obtained in matrix form. Similarly, applying finite element discretization to the equations (2a) and (2b), the nonlinear equation system can be obtained in matrix form. Applying the interface conditions (3a) and (3b) to these spatially discretized governing equations, the monolithic equation system can be obtained as follows: The equilibrium equation as where M is the mass matrix, C is the diffusive matrix, G is the divergence operator matrix, N is the convective term vector, q is the elastic internal force vector, g is the external force vector, a is the acceleration vector, v is the velocity vector, u is the displacement vector, p is the pressure vector, the subscript L indicates the lumping of the matrix, and the subscript τ indicates the transpose of the matrix. The monolithic equation system (5) is linearized using the increments of the state variables to obtain the following equations: where the pressure term and the elastic interior force term are evaluated implicitly, M * is the generalized mass matrix, Δ denotes the time increment, t denotes the current time, Δg and Δh are the residual vectors of the equilibrium equation (5a) and the incompressibility constraint (5b), respectively, G ε is come from the pressure stabilization term of the PSPG method [46], and the relations between the increments of the state variables are given as Δu=βΔt 2 Δa and Δv=γΔtΔa based on Newmark's β method.
The fluid convection and diffusion terms are evaluated explicitly to reduce the computational cost. Therefore, M * is given as L M+βΔt 2 K, while the Courant's number condition and the Diffusion number condition are imposed on the time increment Δt for the stability of the time integration. The predictormulticorrector algorithm is used for the time integration.

Projection method
The monolithic method solves the linearized monolithic equation system (6) and satisfies the interface conditions (3a) and (3b) to avoid spurious numerical power on the interface, which yields numerical instability. However, the formulation leads to an illconditioned equation system. Therefore, in this study, the projection method using the algebraic splitting [43] is used to avoid this difficulty. The present method was first proposed in [41], and modified in [42,43]. The present method is briefly described as follows: The state variables are predicted as the intermediate state variables from the equilibrium equation (5a) for the known pressure, which is linearized as where â is the intermediate acceleration or the prediction of the acceleration, and a Dˆis its increment from the known acceleration. Subtracting both sides of equation (7) from those of equation (6a) gives, after suitable rearrangement, where v is the intermediate velocity or the prediction of the velocity, and the relations between the increments of the intermediate state variables as well as those for the state variables are given based on Newmark's β method. Left multiplying both sides of equation (8) by τ G L M −1 , the following equation is obtained: is solved, then equation (9) is reduced as Since the linear convergence of the state variables is expected for the present definition of M * , v agrees with v asymptotically in the nonlinear iterations. Therefore, the second term of equation (11) will vanish asymptotically, and the incompressibility constraint for the unknown fluid velocity is satisfied. It follows from the above formulation that the monolithic equation system is split into the equilibrium equations (6a) and (7) and the PPE (10), where the intermediate velocity is used and, different from the other studies using block-LU factorization or substructuring, the Schur complement matrix is never produced. The proposed method is summarized as follows: in the nonlinear iterations, the linearized equilibrium equation for the previous pressure (7) is solved to derive the intermediate velocity, the PPE (10) is solved to determine the current pressure such that the current velocity satisfies the incompressibility constraint, and the linearized equilibrium equation for the current pressure (6a) is solved to derive the current velocity.

Mesh model of the interface
In this study, mixed interpolation of tensorial components (MITC) shell elements (four nodes in each element) [44,45] are used for a thin elastic structure, while linear equal-order-interpolation velocity-pressure elements (four nodes in each elements) [46] are used for the fluid. The same arrangements of the fluid and shell nodes can be used [54][55][56]. Therefore, a pair of two overlapping fluid nodes per shell node on the interface is used to describe the pressure discontinuity between the obverse and reverse of the shell structure. Then the interface conditions can be written in matrix form as where v, Q, and g are the nodal components of the velocity, and the internal and external force vectors, respectively, and the subscripts o and r indicate the obverse and reverse sides, respectively.

Analysis setup
The setup of the numerical analysis is identical to that of the dynamically scaled experiment in [15]. The rectangular wing in figure 3(b) flaps in a rectangular tank containing fluid in figure 3(a). In the corresponding experiment, the tank is sealed at the fluid surface to prevent wave generation. Therefore, the noslip condition is imposed on all outer boundaries of the fluid domain. The model wing consists of the leading edge, the plate spring, and the wing plate as shown in figure 3(b). The sectional shape of the model wing is symmetrical with respect to the wing chord.
The leading edge is implemented as a moving boundary, while the plate spring and the wing plate are implemented as a flexible shell structure. The wing model's longitudinal length is L w =22.5 cm and its chord length is c m =4.3 cm, which is given using the aspect ratio of the actual cranefly wing [17], the leading edge has a rectangular cross section with dimensions of 0.7 cm×0.7 cm, the dimensions of the plate spring are given as a length in the wing's longitudinal direction, l ps =22.5 cm, a length in the wing's chord direction, c ps =1.0 cm, and a thickness of t ps =0.03 cm, the dimensions of the wing plate are given as a length in the wing's longitudinal direction, l wp =22.5 cm, a length in the wing's chord direction, c wp =2.6 cm, and a thickness of t wp =0.1 cm. The stroke plane angle Φ is set at 123 deg [5], the flapping frequency f j is set at 0.521 Hz, the flapping angular displacement j is initially set at −Φ/2. In the corresponding experiment, the material properties of the plate spring are a mass density of 0.36 s ps r = g cm −3 , a Young's modulus of E ps =6.56×10 9 g (cm −1 s 2 ), and a Poisson's ratio of ν ps =0.35, the material properties of the wing plate are a mass density of 1.4 s wp r = g cm −3 , a Young's modulus of E wp =2.45×10 10 g (cm −1 s 2 ), and a Poisson's ratio of ν wp =0. 4.
The values of Re, St, and Ch using the above setup are equivalent to those in section 2.3, while the value of M is about 1% of that in section 2.3. This is because of the limitation of available solid materials. The material properties of the model wing from the corresponding experiment described in the above are used for the purpose of demonstrating the validity of the present numerical framework in this section. In contrast, the wing mass m w in section 3 is given so as to satisfy the condition of M in section 2.3. In the computer implementation, m w was divided into the leading edge mass (35%), the plate spring mass (34%), and the wing plate mass (31%) following the mass distribution along the wing chord measured in an actual dipteran insect's wing [7].
The plate spring and the wing plane are modeled by MITC shell elements (number of nodes: 121, number of elements: 100), while the fluid domain is modeled by stabilized linear equal-order-interpolation velocity-pressure elements (number of nodes: 46 911, number of elements: 254 352). The time increment Δt is set at T j /5 000. The consumed time in each analysis time step was about 4.7 s (average from 6 to 7cycles) using a multiple core processing (10 core Xeon 2.8 GHz×2CPU, 32GB memory).

Validity of the present numerical framework
The present numerical framework consists of the projection method using the algebraic splitting for the monolithic equation system described in section 2.4.2, the mesh model of the interface described in section 2.4.3, the dynamic similarity law to conduct the numerical simulation described in section 2.3, and the setup in section 2.5. Figures 4(a) and (b) show the time histories of the pitch angle θ and the lift force F L (the vertical component of the fluid force acting on the model wing) from the numerical results using the above setup. The finer mesh that had twice the division in each direction was also used to check mesh convergence. As shown in these figures, the numerical results from the mesh in figure 3 were very close to those from the finer mesh. Therefore, we use the mesh in figure 3 in the following section. Figure 4(c) shows the relationship between the acceleration time t a and the mean of F L . As shown in this figure, the numerical results were very close to the experimental results. Therefore, the present numerical framework has the enough accuracy for the present simulation.

Varying fluid viscosity
The model wing was simulated under various fluid viscosities to investigate the effect of viscosity on its behavior. To vary the fluid viscosity, Re was set to different values in the range of 64-1270, which is the typical Re range for insect flight, while the other dimensionless numbers were maintained at their reference values (St=0.073, Ca=0.19, and M=16, as described in section 2.3). The acceleration time was set to t a =T j /4. The vorticity of the leading-edge vortex was found to increase with decreasing fluid viscosity, as shown in figures 5 and 6. As a result, the mean lift increased by approximately 13% when the fluid viscosity was decreased by increasing Re to 1270, as shown in figures 7(a) and (b). However, the variation in the mean magnitude of the fluid surface force was less than approximately 4%. This is because the mean magnitude of the drag decreased by approximately 9% when the fluid viscosity was decreased by increasing Re to 1270, as shown in figures 7(a) and (b). Therefore, the variation in the mean pitch angular displacement was no more than 2%. Figure 7(c) shows the pitch histories in Re=64, 85, 127, 254, 508, 762, and 1270. As shown in this figure, they are hardly distinguishable with each other. It follows from these results that the fluid viscous force has little effect on the pitch vibration of flapping insect wings in the range of Re considered in this study.

Varying fluid density
In the case dynamically similar to an actual insect wing with the acceleration time t a =T j /4, where the nondimensional numbers were maintained at their reference values (St=0.073, Re=254, Ca=0.19, and M=16, as described in section 2.3), a characteristic peak in θ appeared after the wing rotation during stroke reversal, as shown in figure 8(a) (the case of 100% fluid density). The recoil occurred after this peak, and a high angle of attack was subsequently maintained during the middle of each half-stroke. These characteristics of the present pitch history are remarkably similar to those from the detailed image analysis of the deforming wing kinematics of actual insects during hovering flight conducted in a previous study [9]. Note that any pitching motion was not prescribed in the present model wing. Therefore, their remarkable similarity demonstrated that the pitch motion in actual insect wings can be created passively. The wing rotation can be caused by the inertial moment [7,12], and a high angle of attack can be maintained by ensuring the elastic and aerodynamic moments are in equilibrium [16]. However, the mechanism of the recoil has not yet been elucidated. The difficulty lies in the process of transitioning from the wing rotation to the maintenance of a high angle of attack, when the inertial, elastic, and aerodynamic effects all appear. Therefore, the surrounding air was rarefied to decrease the aerodynamic effect from that in the case dynamically similar to an actual insect as follows: Ca, M, and Re represent the ratios of the dynamic pressure to the elastic force [39], the mass of the wing to the additional mass due to the fluid [40], and the fluid inertial force to the fluid viscous force [38], respectively. Therefore, to rarefy the fluid, Ca and M were respectively decreased and increased at the same rate from their reference values while Re was maintained at its reference value. For example, to obtain a fluid density equal to 10% of the reference fluid density while keeping all other fluid characteristics constant, Ca was decreased from 0.19 to 0.019, M was increased from 16 to 160, and Re was maintained at its reference value of 256.
Initially, the acceleration time t a =T j /4 was considered. The pitch motion at 10% of the reference fluid density was predominantly a free vibration as shown in figure 8(a). This free pitch vibration would be caused by the inertial force due to the flapping   acceleration during stroke reversal. This inertial force acted on the wing like step loading, since the acceleration time t a is close to the natural period T n of the pitch mode of the model wing, which is about 0.24 times of T j . The center of the wing's mass was below the torsional axis. Therefore, the inertial force rotated the wing around its torsional axis. The evidence of the free pitch vibration is that the frequency of the observed vibration in 10% fluid density is very close to the natural frequency of the pitch mode of the model wing, which is about 4.1 times of f j . The continuous shift of the characteristic peak of θ in the vicinity of the stroke reversal (circles in figure 8(a)) from 100% to 10% of the reference fluid density indicates that the wing rotation during stroke reversal and the subsequent recoil are based on the free pitch vibration.
As shown in figure 8(a), the free vibration was significantly damped by the surrounding air, and was reduced to the transient response very close to critical damping, i.e., the free vibration almost dissipated until the middle of each half-stroke, where the pitch angular displacement θ converged to satisfy the equilibrium between the elastic and dynamic pressure moments. Note that this significant damping is a result of the dynamic pressure resistance of the surrounding fluid because the contribution of the fluid viscous force to the pitch vibration is very small compared to that of the fluid dynamic pressure, as discussed in section 3.1. In this vibration, the pitch angle monotonically decreased from the value of the characteristic peak after the stroke reversal to a large pitch angle during the middle of each half stroke. This phenomenon corresponds to the recoil in actual insect (see the schematic view in figure 1(b)). Therefore, the recoil can be understood as the residual of the free pitch vibration caused by the flapping acceleration in the transient response very close to critical damping resulting from the dynamic pressure resistance of the surrounding air.
It follows from this understanding that the underlying free pitch vibration is significant in cases of apparent recoil, as in the case shown in figure 8(a). In the case of the acceleration time t a =T j /2.5, which yields a smaller acceleration during stroke reversal (figure 2), the recoil does not appear, because the underlying free pitch vibration is not significant, as shown in figure 8(b). It is interesting to investigate the effect of Re on the pitch vibration in the rarefied fluid. For example, at very small Re, can the fluid viscous damping create the characteristic pitching motion instead of the dynamic pressure damping? This question will be investigated in our future work.

Varying acceleration time
As discussed in section 3.2, the recoil can be understood as the residual of the free pitch vibration caused by the flapping acceleration in the transient response very close to critical damping resulting from the dynamic pressure resistance of the surrounding air. On the basis of this understanding of the recoil, the leading-edge vortex and the lift can be controlled by varying the acceleration time t a . The acceleration time t a was thus varied for the case dynamically similar to an actual insect wing, where the nondimensional numbers were maintained at their reference values as St=0.073, Re=254, Ca=0.19, and M=16. The recoil corresponds to the slope of the pitch history that follows the characteristic peak after the stroke reversal as shown figure 1(b). Therefore, in this study, the recoil is said to be strong when this slope of the pitch history is significant. As shown in figures 9(a) and (c), the recoil strengthened with decreasing the acceleration time t a , or increasing acceleration during stroke reversal. This is because the increase in the acceleration during stroke reversal enhanced the underlying free vibration. Figures 10(a)-(e) show the time histories of the lift, the pitch angle, and the flapping speed in the acceleration time t a =T j /4.5, T j /4, T j /3.5, T j /3, and T j / 2.5, respectively, and figure 10(f) shows the relationships between t a and the lift at the start, middle, and end times of the constant flapping velocity. As shown in these figures, the lift at the middle and end times of the constant flapping velocity increased with increasing t a due to increasing the constant flapping speed. Similarly, the lift at the start time of the constant flapping velocity increased with increasing t a in the range of t a >=T j /3.5. On the other hand, the lift at the start time of the constant flapping velocity increased with decreasing t a in the range of t a <=T j /3.5, irrespective of decreasing the constant flapping speed.  The vertical dotted lines divide each half-stroke into thirds. The vertical solid lines indicate times of 6.51, 6.53, 6.54, 6.55, 6.56, 6.58, 6.60, 6.61, 6.63, 6.64, and 6.66 cycles. Note that 6.61 cycle is the start time of the constant flapping velocity in T j /4.5 and 6.64 cycle is the end time of the constant flapping velocity in T j /3.5, which are indicated by arrow. Figures 11 and 12 show the flow fields around the wing for t a =T j /4.5 and T j /3.5, respectively, during the first third of the half stroke (6.51-6.66 cycle). As shown in these figures, the leading-edge vortex in t a =T j /4.5 developed quickly compared with that in t a =T j /3.5 after starting the recoil in t a =T j /4.5 at about 6.55 cycle, while the leading-edge vortex in t a =T j /3.5 developed slowly throughout the first third of the half stroke.
The recoil induces the rotational flow around the wing so as to enhance the leading-edge vortex. The recoil in t a =T j /4.5 was stronger than that in t a =T j /3.5, as shown in figure 9(c), and the enhancement of the leading-edge vortex due to the recoil in t a =T j /4.5 was stronger than that in t a =T j /3.5. Therefore, the leading-edge vortex in t a =T j /4.5 developed quickly compared with that in t a =T j /3.5 after starting the recoil in t a =T j /4.5 at about 6.55 cycle, and the lift in t a =T j /4.5 was larger than that in t a =T j /3.5 during a large part of the first third of the half stroke, as shown in figures 9(b) and (d). This effect of the recoil on the leading-edge vortex will explain the reason why the lift at the start time of the constant flapping speed increased with decreasing t a in the range of t a <=T j /3.5, irrespective decreasing the constant flapping speed, as shown in figure 10(f).
As shown in figures 10(a) and (b), two characteristic peaks appeared in the lift histories for t a =T j /4.5 and T j /4, which can be explained as follows: in t a =T j /3.5, the lift increased during the period of the constant flapping speed, and, after this period, it decreased rapidly due to decreasing the flapping speed rapidly, as shown in figure 10(c). Therefore, the lift peak was created during the second half of each half stroke. Similarly, in t a =T j /4.5 and T j /4, the lift peak during the second half of each half stroke would be created due to decreasing the flapping speed rapidly after the period of the constant flapping speed. Moreover, the lift peak during the first half of each half stroke would be created due to the effect of the recoil on the leading-edge vortex.
The lift decreased in the middle of each half-stroke with decreasing t a , as shown in figures 9(b) and 10(f), because the translational speed decreased to maintain a stroke plane angle constant (figure 2). As shown in figure 13, however, the decrease in the mean lift did not exceed 2%, and the mean lift was nearly independent of t a . This is because the increase in the lift due to the recoil was comparable with the decrease in the translational lift ( figure 13). As with the mean lift, the mean drag was found to be nearly independent of t a .
The nondimensional work done by a flapping wing during the nth cycle of the stroke is defined as The mean W n decreases monotonically by approximately 10% when t a decreases from T j /2.5 to T j /4.5 as shown in figure 14, whereas the mean of the total lift remains approximately constant due to the recoil effect. As shown in figure 2, the flapping motion at smaller t a is similar to that for an actual insect wing. Therefore, the kinematic characteristics of the flapping motion contribute to the energy efficiency without sacrificing the hovering performance provided by the pitch motion based on the passive mechanisms.

Concluding remarks
In this study, we revealed a passive recoil mechanism by conducting a direct numerical simulation based on the dynamic similarity law for FSI. The recoil is the residual of the free pitch vibration induced by the flapping acceleration in the transient response very Figure 11. Fluid velocity and vorticity fields for t a =T j /4.5 at (a) 6.51 cycle, (b) 6.53 cycle, (c) 6.54 cycle, (d) 6.55 cycle, (e) 6.56 cycle, (f) 6.57 cycle, (d) 6.58 cycle, (e) 6.60 cycle, (f) 6.61 cycle, (g) 6.63 cycle, (h) 6.64 cycle, and (i) 6.66 cycle on a cylindrical plane whose radius divided by the wing longitudinal length is equal to r 2 of the model wing. The color and direction of each arrow represent the magnitude and direction of the fluid velocity, respectively, with the magnitude ranging from 0 (blue) to 50 cm s −1 (pink). The color contours in the background of the plots indicate the magnitude of the x-component of the fluid vorticity, ranging from −80 (blue) to 80 s −1 (pink), of which color legend is shown in figures 5 and 6.
close to critical damping due to the dynamic pressure resistance of the surrounding air. The leading-edge vortex and the lift were controlled on the basis of this understanding. As a result, the recoil effect or the Figure 12. Fluid velocity and vorticity fields for t a =T j /3.5 at (a) 6.51 cycle, (b) 6.53 cycle, (c) 6.54 cycle, (d) 6.55 cycle, (e) 6.56 cycle, (f) 6.57 cycle, (g) 6.58 cycle, (h) 6.60 cycle, (i) 6.61 cycle, (j) 6.63 cycle, (k) 6.64 cycle, and (l) 6.66 cycle on a cylindrical plane whose radius divided by the wing longitudinal length is equal to r 2 of the model wing. The color and direction of each arrow represent the magnitude and direction of the fluid velocity, respectively, with the magnitude ranging from 0 (blue) to 50 cm s −1 (pink). The color contours in the background of the plots indicate the magnitude of the x-component of the fluid vorticity, ranging from −80 (blue) to 80 s −1 (pink), of which color legend is shown in figures 5 and 6.  Relationship between t a /T j and nondimensional work. The nondimensional work was averaged from 1 to 10 cycle. enhancement of the leading-edge vortex and the increase in the lift due to the recoil were elucidated. Because the recoil of the wing is the process of transitioning from the wing rotation during stroke reversal to the maintenance of a high angle of attack during the middle of each half-stroke, its passive mechanism naturally includes the known passive mechanisms of these two processes. On the basis of this understanding of the passive mechanisms that are at play during the entire stroke, we found that the kinematical characteristics of the actual flapping motion contribute to the energy efficiency without sacrificing the hovering performance due to the recoil effect invoked by them. These characteristics can reduce the work done by flapping wings and, at the same time, enhance the recoil. The enhanced recoil increases the lift during the first half of each halfstroke so as to maintain the mean of the total lift.