A two-dimensional computational study on the fluid–structure interaction cause of wing pitch changes in dipteran flapping flight

Flying insects depend on unsteady fluid dynamic effects to generate lift in their flapping flights. The majority of the experimental(Ellington et al., 1996; Dickinson et al., 1999; Usherwood and Ellington, 2002; Birch and Dickinson, 2003) and numerical (Liu et al., 1998; Wang, 2000; Sun and Tang, 2002; Rammamurti and Sandberg, 2002; Wang et al.,2004; Luo and Sun,2005; Miller and Peskin,2005) investigations have been performed with the rigid wing models that are given combined active flapping and pitching motions. These investigations have found that the unsteady fluid dynamic effects are significantly affected by the kinematical characteristics of the pitch motions such as the phase lag between pitching and flapping(Dickinson et al., 1999). On the other hand, many observations have reported the flexibility of the insect wings during flapping flights with various modes of deformation. The most significant among them is the high torsional flexibility in Diptera which is concentrated on the narrow wing basal and short root regions. The consequences of this are twofold: (1) it is hard to transmit the active torsion applied by the muscle to the outer wing and (2) the wing can twist easily to provide a passive change in the angle of attack(Ennos, 1987) or passive pitch motion.

In general, the passive pitch motion of the wings that solves the equation of motion balances three forces: wing inertial, elastic and fluid forces. According to Ennos (Ennos,1988b), the wing inertial force accounts for much of the wing rotation at the stroke reversals in dipteran flapping flight. Although the wing inertial force might be dominant at the stroke reversals, where the acceleration is large, it is not enough to explain the passive pitch motion during the parts of the stroke other than the wing reversals. Investigation of the passive pitch motion during the entire stroke should be performed by taking into account all three of these forces. From the view point of the fluid surrounding the wings, these forces exerted on the wings act back on the fluid, producing a cycle of interaction between the fluid and the wings,termed the fluid–structure interaction (FSI), through the equilibrium on their interface.

The aim of this paper was to reveal the details of the passive pitch motion due to wing torsional flexibility and its effects on lift generation by using(a) the non-linear FSI finite element method to analyze the precise motion of the wing passive pitching and the surrounding fluid flow, (b) characterization of the insect flight using a FSI similarity law, (c) the lumped torsional flexibility model as a simplified dipteran wing and (d) the analytical wing model.

The finite element method (Ishihara and Yoshimura, 2005; Ishihara et al., 2008) was used to solve the motion of the wing model which undergoes active sinusoidal flapping motion and passive pitching in the two-dimensional fluid. The numerical simulations are guided by the FSI similarity law (D.I., M.D. and T.H., manuscript submitted), which was used to correctly incorporate data from the selected real insect, the crane fly Tipula obsolete (Ellington,1984a; Ellington,1984b), and the robotic fly(Wang et al., 2004) into our wing model. The lumped torsional flexibility model consists of a plate with unit extent in the z-direction and an attached spring. The former models the wing cross-section averaged over the wing span. The latter models the concentrated torsional flexibility of the wing and allows the wing to twist around its torsional axis. We also introduced the analytical wing model,a single degree of freedom mass–spring–dashpot system, to explain characteristics of the passive pitching motion in the simulation.

The passive pitching motion of our wing model simulates well the real insect's pitching motion. The model wing keeps the high attack angle during its translation and rotates at the stroke reversals. The resulting pitch angle is approximately equal to that of the selected insect, the crane fly. It is especially important that the model wing begins to twist before it changes its flapping direction. Such advanced pitch motion is necessary in order for the insect flight to increase the lift force(Dickinson et al., 1999) and is widely observed in insect flight. The analytical wing model explains the advanced pitch motion of the passive pitching as well as the torsion wave observed in the dipteran wing (Ennos,1988b). The lift force generated by such passive pitching almost meets the force required to support the weight of the crane fly, but it is not quite enough. This could be attributed to the loosely attached leading edge vortex on the wing due to the long wing chord travel of the crane fly for the two-dimensional simulation. Indeed the lift of the fictitious insect flight,whose Reynolds number and stroke–wing chord ratio are much closer to those of the fruit fly model (Wang et al.,2004), shows a 35% increase compared with that of the crane fly flight due to the more tightly attached leading edge vortex on the wing.

To analyze wing motion, the surrounding fluid motion and their interaction,we used the finite element FSI analysis method(Ishihara and Yoshimura, 2005; Ishihara et al., 2008) based on the monolithic algorithm (Zhang and Hisada, 2001; Rugonyi and Bathe, 2001), which simultaneously solves the following equations.

• Navier–Stokes equations for the incompressible viscous fluid:

  \[\ {\rho}^{\mathrm{f}}\frac{{\partial}v_{i}^{\mathrm{f}}}{{\partial}t}+{\rho}^{\mathrm{f}}v_{j}^{\mathrm{f}}\frac{{\partial}v_{i}^{\mathrm{f}}}{{\partial}x_{j}}=\frac{{\partial}{\sigma}_{ji}^{\mathrm{f}}}{{\partial}x_{j}}\]

  (1a)

  and

  \[\ \frac{{\partial}v_{i}^{\mathrm{f}}}{{\partial}x_{i}}=0,\]

  (1b)

  where ρ, v_i and σ_ij are the mass density, velocity and stress with the superscript f indicating the fluid. The fluid is assumed to be Newtonian. The stress

  \({\sigma}_{ij}^{\mathrm{f}}=p{\delta}_{ij}+{\mu}({\partial}v_{i}{/}{\partial}x_{j}+{\partial}v_{j}{/}{\partial}x_{i})\)

  is used, where δ_ij is the Kronecker delta, μ is the fluid viscosity and p is the fluid pressure. Body forces acting on the fluid are ignored for simplicity.
• Equation of motion for the elastic body:

  \[\ {\rho}^{\mathrm{s}}\frac{\mathrm{d}v_{i}^{\mathrm{s}}}{\mathrm{d}t}=\frac{{\partial}{\sigma}_{ji}^{\mathrm{s}}}{{\partial}x_{j}},\]

  (2)

  where the superscript s refers to the structural quantity. We ignore the body forces acting on the elastic body. Let us consider the equation of motion of a flapping insect wing. Assuming that the insect is hovering perfectly, the net upward aerodynamic force acting on the wing will be equal in magnitude and opposite in direction to the gravitational force on the insect. Note that the wing weight is 2–4% of the body weight in the crane fly according to Ellington (Ellington, 1984a),whose data are used in this study. Thus the gravitational force on the wing can be ignored and the body force term is dropped in Eqn 2 describing the wing motion. In this study, the second Piola–Kirchhoff stress and the Green–Lagrange strain tensors are used. The latter is defined by

  \(e_{ij}^{\mathrm{s}}=1{/}2[{\partial}u_{i}{/}{\partial}x_{j}+{\partial}u_{j}{/}{\partial}x_{i}+({\partial}u_{k}{/}{\partial}x_{i})({\partial}u_{k}{/}{\partial}x_{j})]\)

  ⁠. We assume that the increments of these tensors are related by Hooke's law and the material non-linearity observed in biomaterials is ignored. This assumption is justified by the torsional test for the dipteran wing by Ennos(Ennos, 1988a). The test result for the wing without immobilization of its basal and root regions has shown that the relationship between the applied torque and the resulting angular displacement is approximately linear in pronation or supination.
• Geometrical compatibility and equilibrium conditions on the fluid–structure interface:

  \[\ v_{i}^{\mathrm{f}}=v_{i}^{\mathrm{s}},\]

  (3a)

  \[\ {\sigma}_{ij}^{\mathrm{f}}n_{j}^{\mathrm{f}}+{\sigma}_{ij}^{\mathrm{s}}n_{j}^{\mathrm{s}}=0,\]

  (3b)

  where

  \(n_{i}^{\mathrm{f}}\)

  and

  \(n_{i}^{\mathrm{s}}\)

  denote the outward unit normal vectors to the fluid and structure, respectively.

The arbitrary Lagrangian–Eulerian (ALE) method(Hughes et al., 1981) is used in order to take into account the deformable fluid–structure interface. For the elastic body, the total Lagrangian formulation is used in order to take into account the geometrical non-linearity due to the large deformation. The finite elements used for the fluid are the stabilized continuous linear velocity and pressure elements (Tezduyar et al., 1992), while those for the elastic plate are the mixed interpolation of the tensor component elements(Bathe and Dvorkin, 1985; Noguchi and Hisada, 1993). The details of the algorithm of the FSI method used in this study as well as its verification for the basic FSI problems are given by Ishihara and Yoshimura(Ishihara and Yoshimura,2005).

In this study, we introduced a dynamic similarity law for the FSI (D.I.,M.D. and T.H., manuscript submitted) using the dimensional analysis of the equations of motion of the FSI system, Eqns 1a, 2 and 3b. Let U, V, L, f and P be the characteristic displacement, velocity, length, frequency and pressure, respectively. The similarity law consists of six non-dimensional numbers: Re=ρ^fVL/μ (the Reynolds number,the ratio between the inertial force due to the convective acceleration and the viscous force), St=fL/V (the Strouhal number,the ratio between the inertial force due to the Eulerian time derivative acceleration and the inertial force due to the convective acceleration), R_s=ρ^sL²f²/E(the ratio between the inertial force due to the Lagrangian time derivative acceleration and the elastic force), R_M=ρ^f/ρ^s (the mass number,the ratio between the fluid inertial force due to the Eulerian time derivative acceleration and the structural inertial force due to the Lagrangian time derivative acceleration), R_Is=fμ/E(the ratio between the fluid viscous and elastic forces) and R_Ip=ρ^fVLf/E (the ratio between the fluid dynamic pressure and the elastic force), where the characteristic pressure P is evaluated by the magnitude of the dynamic pressure ρ^fV². Only four of them are independent due to the relationships R_Ip/R_Is=Re and R_M=R_s×St×R_Ip. The numbers R_s, R_Is and R_Ip are new as far as we know.

Two systems that are geometrically similar are also dynamically similar if the non-dimensional numbers introduced above are equal. The FSI similarity law was verified through the numerical examples of (a) the forced vibration of an elastic cantilever plate in the quiescent fluid (the problem analyzed in this paper in the context of flapping flight), and (b) the flow-induced vibration of an elastic cantilever plate situated in the wake of a rectangular prism in the uniform flow (D.I., M.D. and T.H., manuscript submitted). In each example,five dynamically similar systems with different fluid and/or reference length were analyzed by our finite element FSI method. In order to verify the validity of the present law, all the governing equations were cast in the dimensional form deliberately. The dynamic similarity of these systems was successfully verified, with the relative error of the results being less than 1% in all cases.

In the dimensional analysis of the equations of the motion of the FSI system, we assumed the linear strain tensor e^s_ij=1/2(∂u_i/∂x_j+∂u_j/∂x_i)and a linear relationship (Hooke's law) between the stress and strain tensors. We also used an alternative approach to the dimensional analysis of the general FSI system including, but not limited to, the linear elastic equation. We defined the following fundamental variables of the FSI system on which the system motion is dependent: L, U, V, P, ρ^s,ρ ^f, E (Young's modulus) and μ. Applying the standard procedure of the dimensional analysis to these variables, we found the same non-dimensional numbers as before. The former analysis ensures that our law is the general framework of the dynamic similarity law. The latter analysis ensures that our law can be applied to the FSI system with large deformation of the elastic body.

The numerical simulations are guided by the FSI similarity law, which was used to correctly incorporate data from a selected real insect, the crane fly(Ellington, 1984a; Ellington, 1984b), and a robotic fly (Wang et al.,2004) into our wing model. It is very interesting how these parameters vary for a variety of flying insects and their flight performances change as each number is varied. In the numerical experiments of this study,comparison of the leading edge vortex behaviors between the crane fly and fictitious insect flight deals to some extent with these questions. The leading edge vortex becomes more tightly attached to the wing as Reand St are decreased while R_M and R_Is are preserved. In future work, we will tackle these interesting questions.

Most of the wing torsional flexibility in Diptera is concentrated on the narrow basal and short root regions of the wings as shown in Fig. 1A, and allows the wing to twist around its torsional axis in a span-wise direction as shown in Fig. 2(Ennos, 1987; Ennos, 1988a). Based on these facts, (1) we modeled the wing two-dimensionally by a rectangular plate, with chord length c and a unit extent in the out-of-plane direction(z-direction), which is the span-wise average of the varying cross-section, and (2) we modeled the torsional flexibility of the wing as an elastic spring, i.e. the torsional flexibility of the wing is lumped into the spring and its torsional stiffness is characterized by the spring constant(the lumped torsional flexibility model). The concept of the lumped torsional flexibility model is shown in Fig. 1B. Initially, the wing is at rest with the angular displacement or pitch angle a=0 in the two-dimensional fluid and then the sinusoidal flapping motion in the x-direction(Eqn 4a) is applied at its leading edge, where it is supported by a roller. The flapping motion is approximated by the sinusoidal function in this study because the actual kinematics of flies and bees are similar to this simple motion(Ellington, 1984b). We used a continuum plate model as shown in Fig. 1C. This is because the plate spring is easily discretized by using finite elements. This model consists of two parts: the upper narrow part of length l_s and the rest with low and high bending rigidities to simulate the spring (s) and the rigid body (r), respectively. The thickness and Young's modulus of the former are h_s and E_s and those for the latter are h_r and infinity, respectively. The infinite Young's modulus is approximated by setting a value far larger than E_s. The mass density is set to a uniform value ρ^s over the entire plate. Few data concerning the wing torsional flexibility are available. Ennos(Ennos, 1988a) has investigated the wing torsional stiffness of some Diptera such as crane, drone and blow flies. The Reynolds number of crane fly flight (Re≈290) is near to that of the robotic fly flight (Re≈100)(Wang et al., 2004), while the Reynolds numbers of the flights of drone and blow flies are over 1000. We selected the crane fly in this study because the applicability of our numerical method to the flapping wing has been verified in the flow regime of low Reynolds number using the robotic fly data (see Appendix). The high torsional flexibility region around the wing base of the crane fly is schematically denoted by the arrows in Fig. 1A. The chord length of this region is approximately 10–30%of the maximum chord length that occurs at mid-span. Since the upper narrow part of the continuum plate model corresponds to the high torsional flexibility region, we define the length l_s of the continuum plate model as 0.2c. Under this setting, we obtain the following non-dimensional numbers in the FSI similarity law.

The parameters of the continuum plate model are determined so that the non-dimensional numbers (Eqns 5 and 6) agree with those for the real insect. The real insect selected was the crane fly(Ellington, 1984a; Ellington, 1984b), with span length of one wing R=1.27 cm, total area of the wing pair S=0.59 cm², stroke angle φ=123 deg., flapping frequency f=45.5 Hz and mass of the wing pair m_w=0.000245 g. The total mass of the insect is m=0.0114 g. The averaged chord length c is given by S/(2R)=0.23 cm. Let ρ^s=1.2 g cm^–3 be the wing mass density(Jensen and Weis-Fogh, 1962; Wainwright et al., 1982); the averaged thickness is then h=m_w/(ρ^sS)=0.00069 cm,which is identified as the thickness h_r of the hypothetical rigid part of the real insect. The stroke amplitude A₀ of this flapping is estimated by the arc length traveled by the chord at the mid-span and is given by 0.5R(φ/180)π=1.36 cm. The value of the Young's modulus E_s is 6.1 GPa(Wainwright et al., 1982) and those for the fluid mass density ρ^f and viscosity μ areρ ^f=0.0012 g cm^–3 (air) and μ=0.00018 g(cms)^–1 (air), respectively. The non-dimensional numbers St, Re, R_Is, R_M andγ _r are determined to be 0.054, 290,1.4×10^–13, 1.0×10^–3 and 340,respectively, using the above properties. On the other hand, the numberγ _s is determined based on the expectation that it is of the same order as γ_r; among the numerically tested values ofγ _s in the region of γ_r=340, we have found that γ_s=900 maximizes the averaged lift. To test the validity of this selection of γ_s=900, consider the bending of the continuum plate spring with its upper end fixed and the moment Mapplied at the lower end. The slope of the plate at its lower end, using the Euler–Bernoulli beam assumption, is given by a=Ml_s/(E_sI), where I=h_s³/12 is the second moment of area for the plate spring. Thus the spring constant of the plate is given by k_s=E_sI/l_s=E_sh_s³/(12l_s),which is set to the macroscopic torsional stiffness G_w of the insect wing. Notice that the selected value ofγ _s=c/h_s=900 gives k_s=1.9 g cm² (s²rad)^–1, which lies in the range of values for G_w=0.8–15.4 g cm² (s²rad)^–1 (Ennos,1988a) obtained for the crane fly.

The torsional axis of our model is located within the upper 20% of the wing chord and that of the real insect is located near the leading edge. We investigate the effects of the location of the rotational axis in the Appendix, where we show that they can be neglected with slight perturbation in the force amplitude and phase shift.

In order to explain the characteristics of the passive pitching of the continuum plate model, we introduced a single degree of freedom mass–spring–dashpot system of the wing(Fig. 3) for the theoretical analysis. The wing is modeled by a rigid rod of length c and its angular displacement a(t) is assumed to be small enough that it can be measured by the relative x-displacementδ(t)=u(t)–U(t) of the wing center according to a(t)=δ(t)/(c/2) (assumption of linearization), where U(t) and u(t) are the displacement of the wing base and center, respectively. We assume that the wing mass is concentrated on the wing center and, to be consistent with the scheme of linearization, the resultant force on the wing is assumed to act on the wing center in the x-direction.

Fig. 4 shows the dependence of the phase shift b on the frequency ratio f/f_n for three representative values of the damping ratios ζ=0.5 (under damping), 1 (critical damping) and 2 (over damping). Regardless of the value of ζ, the phase shift b is (a)positive (advanced phase shift) for f/f_n<1,(b) zero (no phase shift) for f/f_n=1 and (c)negative (delayed phase shift) for f/f_n>1,respectively. Note that the amount of the advanced phase shift increases as the value of f/f_n (<1) gets smaller.

Once the chord length c and the fluid material propertiesρ ^f and μ are specified in our continuum plate model, all the other properties, the flapping frequency f, stroke amplitude A₀, plate thicknesses h_r and h_s, plate material properties ρ^s and E_s, are determined to match the non-dimensional numbers of the FSI similarity law Re, St, R_M, R_Is, γ_r and γ_s obtained from the real insect. We used c=6.7 cm and the material properties of mineral oil ρ^f=0.88 g cm^–3 and μ=0.715 g(cm s)^–1 in the scaled computation, following the dynamically scaled model of Dickinson and colleagues(Dickinson et al., 1999). Since our plate is flexible at the upper narrow part, unlike the rigid wing model of Dickinson and colleagues (Dickinson et al.,1999), not only the fluid motion but also the plate motion and their interaction are dynamically scaled to the real insect flapping flight. Fig. 5 shows the schematic view of the analysis domain and the finite element mesh employed, in which the no-slip condition is assumed on the outer boundary of the fluid domain and the wing chord. A convergence test where the fluid grid is refined is provided in the Appendix. The center of the wing is located at the center of the fluid domain initially. We first apply the sinusoidal flapping motion Eqn 4a in the x-direction at the leading edge of the flexible wing to measure the induced passive pitching in terms of the angular displacement or pitch angle a(t) of the trailing edge of the wing. The angular displacement a(t) is measured from the rest position of the wing as shown in Fig. 1C. The anti-clockwise direction is positive. Note that no angular displacement a(t) is applied actively. Next, we consider a rigid wing model and apply both the sinusoidal flapping motion Eqn 4a and the measured angular displacement a(t) actively to calculate the fluid forces. As shown in Fig. 6, the fluid forces obtained in the two cases are in fine agreement. Since the ability of our method to simulate the two-dimensional wing motion under the combined active flapping and pitching for the rigid wing has been verified in the Appendix, this agreement serves as a benchmark to validate our method for the flexible wing that induces passive pitching under active flapping.

Figs 7 and 8 summarize the passive pitching motion of our continuum plate model wing. Fig. 7 also shows the generated lift history. After the initial transient phase, the angular displacement or pitch angle a oscillates regularly after around four cycles as shown in Fig. 7. The initial transient phase can be explained by the equation of motion(Eqn 9) of the analytical wing model. The analytical solution of Eqn 9 consists of the free and forced vibration modes initially. However, the free vibration mode disappears due to the damping effect after multiple strokes and the forced vibration mode remains. Note that the analytical solution of the forced vibration mode is given by Eqn 10. The same phenomenon occurs in the pitching motion of the continuum plate model wing. Following the stabilization of the wing motion, the lift coefficient also oscillates regularly.

The passive pitching motion shown in Fig. 8 seems to simulate well the real insect's pitching motion. The model wing keeps the high attack angle during its translation and rotates at the stroke reversals. In addition, our wing model gives a maximum pitch angle of 50–60 deg., which lies in the range of values for the crane fly, 45–65 deg.(Ellington, 1984b). Note that the pitch angle is a non-dimensional variable. It is especially important that the model wing begins to twist before it changes its flapping direction as shown in Figs 7 and 8. Such advanced pitch motion is important for the insect flight to increase the lift force by intercepting the wake of the previous stroke (wake capture)(Dickinson et al., 1999) and is widely observed in insect flight.

The advanced pitch motion observed in Figs 7 and 8 can be explained by the analytical wing model as follows. The time history of the pitch angle a(t) takes a form similar to a sinusoidal function with the same frequency f as flapping in the regular oscillation after around four cycles. Thus, the time history of pitching in the steady state is approximately expressed by Eqn 4bwith the positive phase shift b. The fact that our analytical model,excited by the sinusoidal force in Eqn 9 and subject to damping, attains the sinusoidal oscillation(Eqn 10) for the relative displacement δ(t)=ca(t)/2 in the steady state provides convincing support for the presence of the sinusoidal passive pitching a(t) in our continuum plate model. As shown in Fig. 4, regardless of the value of ζ, the phase shift b is (a) positive (advanced phase shift)for f/f_n<1, (b) zero (no phase shift) for f/f_n=1 and (c) negative (delayed phase shift) for f/f_n>1, respectively. For our model wing the spring constant k_s′, given by Eqn 8 with k_s=1.9 g cm² (s²rad)^–1 and c=0.23 cm, is equal to 144 g s^–2, while the wing mass m_w=ρ^sc(0.2h_s+0.8h)=0.000166 g. Thus the natural frequency calculated by Eqn 12 is f_n=148 Hz which gives the ratio f/f_n=0.31<1, where f=45.5 Hz. This shows that our continuum plate model wing satisfies the condition f/f_n<1 for advanced pitching.

The analytical wing model also provides an insight into the tip to base torsion wave observed in dipteran flight(Ennos, 1988b). As shown in Fig. 4, the value of the advanced phase shift b increases as the value of f/f_n (<1) gets smaller. Since the natural frequency f_n given by Eqn 12 increases when the wing mass m_w is reduced, the phase shift is bigger for the lighter wing mass. Consequently, when we move from the base to the tip of the wing, the wing mass becomes lighter(Ennos, 1989) and the phase shift becomes larger. Such a phase shift gradient along the span-wise direction appears as the torsion wave in the three-dimensional wing. We need to be careful, however, with the limitations of our continuum plate model of the wing for which the chord length is obtained as an average over the span-wise direction and the motion represented is an average over the span. Nevertheless, our simple analytical model provides an insight into the tip to base torsion wave observed in dipteran flight(Ennos, 1988b).

We found that the lift force generated by passive pitching almost meets the required force to support the weight of the crane fly. Fig. 7B shows the time history of the lift coefficient C_L (normalized by 0.5ρ^fcV_max²) in the wing motion with passive pitching. Its mean value C_L,mean=0.46 gives the mean combined lift force f_L,mean=6.8×10^–5N for both wings of the crane fly, which is comparable to but smaller than the weight of the crane fly w=11×10^–5N. This could be attributed to the loosely attached leading edge vortex on the wing due to the long wing chord travel of the crane fly for the two-dimensional simulation. In the two-dimensional simulation with Re=75–115 and the stroke–wing chord ratio A₀/c=2.8–4.8 for the fruit fly, the wing reverses before the leading edge vortex has time to separate during each stroke. Thus in these cases the additional mechanism is not required to stabilize the leading edge vortex and the lift generated provides a good approximation of the corresponding lift in the three-dimensional experiment (Wang et al.,2004). In our two-dimensional simulation with Re=290 and A₀/c=5.9 (St=0.054) for the crane fly,the leading edge vortex appears to be slightly separated from the wing as shown in Fig. 9. Due to the larger wing travel length or the larger stroke–wing chord ratio and the smaller Strouhal number than those of Wang and colleagues(Wang et al., 2004), the leading edge vortex separates from the wing before the wing reverses in our two-dimensional simulation with no additional stabilization mechanism such as the span-wise flow effect. If we adopt a fictitious insect with Re=200 and A₀/c=3.2 (St=0.1),whose characteristics are closer to those of the model of Wang and colleagues(Wang et al., 2004) than to those of the crane fly, the leading edge vortex is more tightly attached to the wing as shown in Fig. 10. Consequently, the lift coefficient of the fictitious insect is higher than that of the crane fly as shown in Fig. 11 with an approximately 35% increase in the mean lift coefficient(C_L,mean=0.62) compared with that of the crane fly. Pitching in this case is also advanced, as in the case of the crane fly.

It is important to remember here that it is hard to transmit the active torsion applied by the muscle to the outer wing due to the wing torsional flexibility in Diptera (Ennos,1987), while the wing pitch is adequately controlled to produce the lift to enable the insect to hover. Our wing model meets these criteria and explains the mechanism of insect flight with passive pitching. It seems to be natural for insects to select such a passive mechanism to minimize the effort required to move their wings. This passive mechanism might be useful as the design principle for developing micro-air vehicles, i.e. engineers can avoid the need to develop the active mechanism to drive and control the wing pitch adequately.

In the benchmark problem described above, where comparison with the results of Wang and colleagues (Wang et al.,2004) was made, the position of the axis of rotation is located at the center of the wing. In contrast, the position of the rotational axis of our lumped flexibility model is located within the upper 20% of the wing,while that of the real insects is located near the leading edge. Thus the effects of the location of the axis of rotation need to be investigated. We applied U(t) and a(t) at four positions;center (O), top (A), upper 1/10 (B) and upper 1/3 (C). Fig. A3 shows the time history of C_L and C_D for these positions, and Fig. A4 shows the relationship between the rotational axis positions and the second peak value of the force coefficients from 3.5 to 4 cycles in Fig. A3. As shown in Fig. A3, the time histories of C_L and C_D for the four positions are similar. As shown in Fig. A4, the force coefficients vary linearly. C_L varies from 1.46 to 1.65 and C_D varies from 2.26 to 2.71 as the rotational axis position varies from positions C (upper 30% of the wing) to A (top of the wing). The rotational axis of our model is located within the upper 20% of the wing chord and that of real insects is located near the leading edge. As a consequence, the effect of the location of the rotational axis is negligible as long as the rotational axis is located near the top of the wing chord or the leading edge.

A convergence test in which the fluid grid is refined is also provided here. Fig. A5 shows the three fluid grids (meshes O, A and B) employed. Mesh O (the number of nodes and elements is 3417 and 6600, respectively) is used to give the results in Figs A2, A3 and A4. The refined meshes A (the number of nodes and elements is 3905 and 7560, respectively) and B (the number of nodes and elements is 4697 and 9120, respectively) have 1.4 and 2 times finer space resolutions, respectively, around the wing. As shown in Fig. A6, the force coefficients given by using these meshes show very good agreement. As a consequence, we employed mesh O throughout this study.