SUMMARY
In this study, the passive pitching due to wing torsional flexibility and its lift generation in dipteran flight were investigated using (a) the nonlinear finite element method for the fluid–structure interaction, which analyzes the precise motions of the passive pitching of the wing interacting with the surrounding fluid flow, (b) the fluid–structure interaction similarity law, which characterizes insect flight, (c) the lumped torsional flexibility model as a simplified dipteran wing, and (d) the analytical wing model, which explains the characteristics of the passive pitching motion in the simulation. Given sinusoidal flapping with a frequency below the natural frequency of the wing torsion, the resulting passive pitching in the steady state, under fluid damping, is approximately sinusoidal with the advanced phase shift. We demonstrate that the generated lift can support the weight of some Diptera.
 dipteran flight
 dynamic similarity law
 finite element method
 fluid–structure interaction
 wing torsional flexibility
INTRODUCTION
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 nonlinear 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 twodimensional 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 zdirection and an attached spring. The former models the wing crosssection 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 twodimensional 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.
MATERIALS AND METHODS
Nonlinear finite element FSI analysis
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: (1a) and (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 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: (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 . We assume that the increments of these tensors are related by Hooke's law and the material nonlinearity 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: (3a) (3b) where and 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 nonlinearity 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).
We have selected the twodimensional wing model similar to that used by Wang and colleagues (Wang et al., 2004) as a benchmark to test the validity of our method. Miller and Peskin (Miller and Peskin, 2005) have used a similar wing model to demonstrate the validity of their numerical technique. The wing is modeled by a thin rigid elastic plate of thickness h with chord length c; in our elastic model a large Young's modulus is assigned to simulate the rigid body behavior. The model is twodimensional having a unit extent in the zdirection (outofplane direction) without variation in this direction. Following Wang et al. (Wang et al., 2004), an xdisplacement U(t) and an angular displacement around the zdirection a(t): (4a) (4b) are actively applied to the center of the wing, where A_{0},β and b are the stroke amplitude, the amplitude of the pitching angle of attack and the phase shift, respectively. The common frequency f of the flapping and pitching is set to produce the Reynolds number Re=ρ^{f}V_{max}c/μ=75, whereρ ^{f} is the fluid mass density, μ is the fluid viscosity and V_{max}=πfA_{0} is the maximum wing velocity. The details of this test are described in the Appendix, where the time histories of force coefficients given by our method agree well with the results of Wang and colleagues (Wang et al., 2004).
FSI similarity law
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 nondimensional numbers: Re=ρ^{f}VL/μ (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}=ρ^{s}L^{2}f^{2}/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}=ρ^{f}VLf/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 ρ^{f}V^{2}. 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 nondimensional 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 flowinduced 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 nondimensional 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 Re and St are decreased while R_{M} and R_{Is} are preserved. In future work, we will tackle these interesting questions.
Lumped torsional flexibility model
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 spanwise direction as shown in Fig. 2 (Ennos, 1987; Ennos, 1988a). Based on these facts, (1) we modeled the wing twodimensionally by a rectangular plate, with chord length c and a unit extent in the outofplane direction (zdirection), which is the spanwise average of the varying crosssection, 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 twodimensional fluid and then the sinusoidal flapping motion in the xdirection (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 midspan. 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 nondimensional numbers in the FSI similarity law.
Dynamic numbers: (5a) (5b) (5c) (5d)
Geometric numbers: (6a) (6b)
The parameters of the continuum plate model are determined so that the nondimensional 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^{2}, 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 WeisFogh, 1962; Wainwright et al., 1982); the averaged thickness is then h=m_{w}/(ρ^{s}S)=0.00069 cm, which is identified as the thickness h_{r} of the hypothetical rigid part of the real insect. The stroke amplitude A_{0} of this flapping is estimated by the arc length traveled by the chord at the midspan 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 nondimensional 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 M applied 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_{s}I), where I=h_{s}^{3}/12 is the second moment of area for the plate spring. Thus the spring constant of the plate is given by k_{s}=E_{s}I/l_{s}=E_{s}h_{s}^{3}/(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^{2} (s^{2} rad)^{–1}, which lies in the range of values for G_{w}=0.8–15.4 g cm^{2} (s^{2} 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.
Analytical wing model
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 xdisplacementδ (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 xdirection.
The linearized motion of the concentrated wing mass is given by: (7) where the first, second and third terms represent the wing inertial force f_{i}=m_{w}d^{2}u/dt^{2}, the fluid damping force f_{c}=c_{f}(du/dt–dU/dt)= c_{f}dδ(t)/dt and the restoring force f_{s}=k_{s}′[u(t)–U(t)]=k_{s}′δ(t) due to the spring. The coefficient k_{s}′ is the spring constant for the relative xdisplacementδ (t)=u(t)–U(t) and is obtained by linearizing the moment and angular displacement relationship M=k_{s}a according to a=δ/(c/2) and M=f_{s}(c/2), with the result f_{s}=k_{s}′δ, where: (8) In addition, c_{f} is the damping force coefficient due to the fluid and the subtraction dU/dt is the fluid advection effect induced by the wing global flapping motion.
This equation of motion is reduced to the following equation: (9) where F_{0}=m_{w}U_{0}(2πf)^{2} is the amplitude of the periodic exciting force F_{0} given in terms of the amplitude U_{0} (=A_{0}/2) and the frequency f of the spring base displacement U(t) given by Eqn 4a. The analytical solution of Eqn 9 in the steady state is given by: (10) with: (11a) (11b) (11c) where (12) and (13) are the natural frequency and the damping ratio, respectively.
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.
RESULTS AND DISCUSSION
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_{0}, plate thicknesses h_{r} and h_{s}, plate material properties ρ^{s} and E_{s}, are determined to match the nondimensional 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 noslip 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 xdirection 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 anticlockwise 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 twodimensional 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.
Passive pitching motion
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 nondimensional 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 4b with 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^{2} (s^{2} rad)^{–1} and c=0.23 cm, is equal to 144 g s^{–2}, while the wing mass m_{w}=ρ^{s}c(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 spanwise direction appears as the torsion wave in the threedimensional 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 spanwise 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).
Lift generated by wing motion including passive pitching
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ρ^{f}cV_{max}^{2}) 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^{–5}N for both wings of the crane fly, which is comparable to but smaller than the weight of the crane fly w=11×10^{–5}N. 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 twodimensional simulation. In the twodimensional simulation with Re=75–115 and the stroke–wing chord ratio A_{0}/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 threedimensional experiment (Wang et al., 2004). In our twodimensional simulation with Re=290 and A_{0}/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 twodimensional simulation with no additional stabilization mechanism such as the spanwise flow effect. If we adopt a fictitious insect with Re=200 and A_{0}/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 microair vehicles, i.e. engineers can avoid the need to develop the active mechanism to drive and control the wing pitch adequately.
APPENDIX
Testing the validity of our FSI analysis method
We have selected the twodimensional wing model similar to that used by Wang and colleagues (Wang et al., 2004) as a benchmark to test the capability of our method. The wing is modeled by a thin rigid elastic plate of thickness h with chord length c; in our elastic model a large Young's modulus is assigned to simulate the rigid body behavior. The model is twodimensional having the unit extent in the zdirection (outofplane direction) without variation in this direction. Following Wang et al. (Wang et al., 2004), an xdisplacement U(t) and an angular displacement around the zdirection a(t): (A1a) (A1b) are actively applied to the center of the wing as shown in Fig. A1, where A_{0}=2.8c, β=π/4 and b=0 are the stroke amplitude, the amplitude of the pitching angle of attack and the phase shift, respectively. The common frequency f of the flapping and pitching is set to produce the Reynolds number Re=ρ^{f}V_{max}c/μ=75, whereρ ^{f} is the fluid mass density, μ is the fluid viscosity and V_{max}=πfA_{0} is the maximum wing velocity. The wing starts its downstroke with the initial angular displacement a(0)=a_{0}=0. Notice that this wing motion agrees with that used by Wang and colleagues (Wang et al., 2004) with symmetrical pitching. Miller and Peskin (Miller and Peskin, 2005) have used a similar wing model with symmetrical pitching to demonstrate the validity of their numerical technique. Fig. A2 gives the time history of lift coefficient C_{L} and drag coefficient C_{D}, normalized by 0.5ρ^{f}cV_{max}^{2}, which shows very good agreement with the results of Wang and colleagues (Wang et al., 2004). As, in both simulations, the leading edge vortex is attached to the wing during the wing stroke, the lift coefficient agrees well with that obtained from the threedimensional experiment even though the threedimensional axial flow is absent.
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.
FOOTNOTES

D.I. is grateful for support from a GrantinAid from Japan Ministry of Education, Culture, Sports, Science and Technology (No. 17760088).