Unsteady aerodynamic force generation by a model fruit fly wing in flapping motion
Institute of Fluid Mechanics, Beijing University of Aeronautics and Astronautics, Beijing 100083, P.R. China
*e-mail: sunmao{at}public.fhnet.cn.net
Accepted 19 October 2001
![]() |
Summary |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Considerable lift can be produced when the majority of the wing rotation is conducted near the end of a stroke or wing rotation precedes stroke reversal (rotation advanced), and the mean lift coefficient can be more than twice the quasi-steady value. Three mechanisms are responsible for the large lift: the rapid acceleration of the wing at the beginning of a stroke, the absence of stall during the stroke and the fast pitching-up rotation of the wing near the end of the stroke.
When half the wing rotation is conducted near the end of a stroke and half at the beginning of the next stroke (symmetrical rotation), the lift at the beginning and near the end of a stroke becomes smaller because the effects of the first and third mechanisms above are reduced. The mean lift coefficient is smaller than that of the rotation-advanced case, but is still 80 % larger than the quasi-steady value. When the majority of the rotation is delayed until the beginning of the next stroke (rotation delayed), the lift at the beginning and near the end of a stroke becomes very small or even negative because the effect of the first mechanism above is cancelled and the third mechanism does not apply in this case. The mean lift coefficient is much smaller than in the other two cases.
Key words: insect, flight, fruit fly, Drosophila sp., wing, computational fluid dynamics, unsteady aerodynamics.
![]() |
Introduction |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Ellington et al. (1996) performed flow-visualization studies on a hawkmoth Manduca sexta using a mechanical model of the hawkmoth wings. They discovered that the dynamic stall vortex on the wing did not shed during the translational motion of the wing in both the up- and downstrokes because it was stabilized by a strong spanwise flow. The authors suggested that this was a new mechanism of lift enhancement, prolonging the benefit of the delayed stall for the entire stroke. This mechanism of lift enhancement was confirmed by computational fluid-dynamic analyses (Liu et al., 1998
; Lan and Sun, 2001
).
Recently, Dickinson et al. (1999) conducted force measurements on flapping robotic fruit fly wings and showed that a large lift force was maintained during the translational phases of the up- and downstrokes and brief bursts of high lift occurred during stroke reversal when the wing was rotating. They explained the large lift force during the translational phases by the delayed-stall mechanism of Ellington et al. (1996
) and suggested some possible explanations for the large lift force recorded during wing rotation. Only force measurements were conducted by Dickinson et al. (1999
); if simultaneous flow-field information had been obtained, further insights into the high lift generation process could be obtained.
In the present paper, we conduct a computational fluid-dynamic analysis on unsteady aerodynamic force generation for a model fruit fly wing in flapping motion. Unsteady aerodynamic forces and flow fields are obtained simultaneously by numerical solution of the NavierStokes equations. The high lift generation process can be explained further on the basis of force and flow-field information. The flapping motion of normal hovering flight (Fig. 1), similar to that used by Dickinson et al. (1999), is considered.
|
![]() |
Materials and methods |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
|
![]() | (1) |
where ut+=ut/U, =t*U/c, where t* is (dimensional) time,
1 is the time at which the deceleration near the end of a stroke starts and
1+
t is the time at which the acceleration at the beginning of the next stroke finishes. In the calculation, ut+ determines the azimuthal-rotational speed of the wing. Denoting the azimuthal-rotational speed as
, we have
(
)=ut/r0, where
is non-dimensional time. The angle of attack of the wing is denoted by
, which also takes a constant value except at the beginning or near the end of a stroke. Following the experimental work of Dickinson et al. (1999
),
is set as 40° in the present study. Around stroke reversal,
changes with time, and the angular velocity,
, is given by:
![]() | (2) |
where +=
c/U,
0+ is a constant,
r is the non-dimensional time at which the rotation starts and
r is the time interval over which the rotation lasts. In the time interval
r, the wing rotates from
=40° to
=140°. Therefore, when
r is specified,
0+ can be determined (around the next stroke reversal, the wing would rotate from
=140° to
=40°, so the sign of the right-hand side of equation 2 should be reversed). The Reynolds number is defined as Re=cU/
, where
is the kinematic viscosity, and Re=136 is used here, as in Dickinson et al. (1999
), as typical of a fruit fly wing. A stroke angle of 155° is assumed, approximately the same as that used by Dickinson et al. (1999
). On the basis of the stroke angle, the non-dimensional time taken for one cycle (two strokes) is approximately 10.8. In the calculation, the wing starts the flapping motion in still air. The calculation is stopped when periodicity in the aerodynamic force and flow structures is approximately reached.
The NavierStokes equations and the computational method
In the flapping motion considered in the present paper, the wing performs translational motion (azimuthal rotation) and pitching rotation. To calculate the flow around a body performing unsteady motion (such as the present flapping wing), one approach is to write and solve the governing equations in a body-fixed, non-inertial frame of reference with inertial force terms added to the equations. An advantage of this approach is that the coordinate transformation that generates a body-conforming computational grid does not need to vary with time, and the grid is generated only once. But, in this approach, some treatment is needed to handle the far-field boundary conditions when the body is rotating (velocity at the far-field boundary tends to be infinite), introducing extra terms into the equations.
Another approach is to write and solve the governing equations in an inertial frame of reference. By using a time-dependent coordinate transformation and the relationship between the inertial and non-inertial frames of reference, a body-conforming computational grid in the inertial frame of reference (which varies with time) can be obtained from a body-conforming grid in the body-fixed non-inertial frame, which needs to be generated only once. This approach does not need special treatment for the far-field boundary conditions and, moreover, since no extra terms are introduced into the equations, existing numerical methods can be applied directly to the solutions of the equations. This approach is employed here.
The non-dimensionalized three-dimensional incompressible unsteady NavierStokes equations, written in the inertial coordinate system o,x,y,z (see Fig. 1A), are as follows:
| (3) |
| (4) |
| (5) |
| (6) |
where u, v and w are three components of the non-dimensional velocity and p is the non-dimensional pressure. In the non-dimensionalization, U, a constant velocity given above, c, the mean chord length of the wing, and c/U are taken as the reference velocity, length and time, respectively. Equations 36 are transformed from the Cartesian coordinate system (x,y,z,t) to the curvilinear coordinate system (,
,
,
) using a general time-dependent coordinate transformation in the form:
![]() | (7) |
The transformed equations written in conservative form are as follows:
| (8) |
| (9) |
where J is the Jacobian of the transformation and:
![]() | (10) |
![]() | (11) |
![]() | (12) |
| (13) |
| (14) |
| (15) |
| (16) |
| (17) |
| (18) |
| (19) |
where the symbol is the gradient operator, and the velocity gradients and the metrics of the transformation were written as:
| (20) |
| (21) |
Fig. 1A also shows a wing-fixed coordinate system (o',x',y',z'). The inertial coordinates (o,x,y,z) are related to the wing-fixed coordinates (o',x',y',z') through the following relationship:
| (22) |
where a is the distance between o and o'. Using equation 22, the metrics in the inertial coordinate system, (x,
y,
z,
t), (
x,
y,
z,
t), (
x,
y,
z,
t), which are needed in equations 8 and 9, can be calculated from those in the body-fixed, non-inertial coordinate system, (
x',
y',
z'), (
x',
y',
z'), (
x',
y',
z'), which need to be calculated only once. As the wing moves, the coordinate transformation functions vary with (x,y,z,t) such that the grid system moves and always fits the wing. The wing-fixed non-inertial frame of reference (o',x',y',z') is used in the initial grid generation and also in the description of the calculated results.
Equations 8 and 9 are solved using the algorithm developed by Rogers and Kwak (1990) and Rogers et al. (1991
) and are in the same form as that solved by Rogers et al. (1991
). The algorithm is based on the method of artificial compressibility, which introduces a pseudotime derivative of pressure into the continuity equation. Time accuracy in the numerical solutions is achieved by subiterating in pseudotime for each physical time step. The algorithm uses a third-order flux-difference splitting technique for the convective terms and the second-order central difference for the viscous terms. The time derivatives in the momentum equation are differenced using a second-order, three-point, backward-difference formula. The algorithm is implicit and has second-order spatial and time accuracy. For details of the algorithm, see Rogers and Kwak (1990
) and Rogers et al. (1991
).
At the inflow boundary, the velocity components are specified as free-stream conditions while the pressure is extrapolated from the interior. At the outflow boundary, the pressure is set equal to the free-stream static pressure and the velocity is extrapolated from the interior. On the wing surface, impermeable wall and no-slip boundary conditions were applied, and the pressure on the boundary is obtained through the normal component of the momentum equation.
A body-conforming grid was generated by using a Poisson solver based on the work of Hilgenstock (1988). The grid topology used was an O-H grid. A portion of the grid used for the wing is shown in Fig. 2.
A code based on the above numerical method was written by Lan and Sun (2001), and it was verified by the analytical solution of the boundary layer flow on a flat plate and validated by comparing the calculated and measured pressure distributions on a wing. When the code is used in the present study, before making observations pertaining to the physical aspect of the flow, estimates of the accuracy of the computed solutions must be provided. In the present calculations, numerical uncertainties are mainly associated with time discretization (i.e. time step values), spatial resolution and far-field boundary location. The effects of these variables on the computed solutions will be discussed below together with the calculated results.
Evaluation of the aerodynamic forces
Once the NavierStokes equations have been numerically solved, the velocity components and pressure at discretized grid points for each time step are available. The aerodynamic force acting on the wing arises from the surface pressure and viscous stress along the wing surface. Integrating the pressure and viscous stress over the wing surface at a time step gives the total force acting on the wing at the corresponding time instant. The lift of the wing, L, is the component of the total force perpendicular to the translational velocity and is positive when it is directed upwards. The drag, D, is the component of the total force parallel to the translational velocity and is positive when directed opposite to the direction of the translational velocity of the downstroke. The lift and drag coefficients, denoted by CL and CD, respectively, are defined as follows:
| (23) |
| (24) |
where is the fluid density.
![]() |
Results and discussion |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
|
|
|
Variations in CL and CD over a flapping cycle in the present study were similar to those of Dickinson et al. (1999) (see their Figs 1D and 3A). A comparison of the calculated lift with the measured values will be given below.
As noted above, since the flow conditions in the down- and upstrokes are essentially the same, only one stroke is discussed here. At the beginning of the stroke (=16.718), CL increases to peak at approximately 1.2, and variation in CD is small. During this period, the wing accelerates towards the right and rotates clockwise by only a small angle, i.e. from
=134.4° at
=16.7 to
=140° at
=17.5 (see Fig. 4A,B); thus, its motion is mainly translational acceleration. In fast translational acceleration, as shown by Lan and Sun (2001
), a peak in CL could be produced as a result of the generation of strong vorticity layers in a short period, giving a large rate of change of fluid impulse. As can be seen in Fig. 4C, at the end of the acceleration, a new layer of positive vorticity has formed around the leading edge and upper surface of the wing, and a new layer of negative vorticity has formed on the lower surface of the wing extending beyond its trailing edge (starting vortex).
Between =18 and
=19.4, the wing translates (azimuthally rotates) with constant speed at
'=40° (where
'=180°
and
=140°); a large CL, approximately 1.3, is maintained. If the wing is not in azimuthal rotation but is moving like an aeroplane wing, as shown by Lan and Sun (2001
), stall would accur at approximately two chord lengths of travel after starting, or at
18.7, and CL would start to decrease at about this time. In the present case (Figs 4CE, 5AC), the dynamic stall vortex does not shed, stall is absent and a large CL can therefore be maintained; this is the stall-absence mechanism revealed by Ellington et al. (1996
).
Near the end of the stroke, between =19.4 and
=22.1, large values of CL and CD occur, peaking at
20.6 (Fig. 3). From
=19.4 to
=20.6, while still translating with constant speed, the wing rotates (pitches up) rapidly. This motion, fast pitching-up while translating with constant speed, results in sharp increases in CL and CD, also due to the generation of strong vorticity layers over a short period (Lan and Sun, 2001
). As seen in Fig. 4F, a new layer of positive vorticity is formed around the leading edge of the wing and a new layer of negative vorticity near the trailing edge of the wing. From
=20.8 to
=22.1, the wing is in deceleration, causing sharp decreases in CL and CD (Lan and Sun, 2001
). The above discussion shows that pitching-up rotation causes CL and CD to increase rapidly, and deceleration causes them to decrease rapidly, forming the large peaks in CL and CD in the last part of the stroke.
When the wing moves like an aeroplane wing under steady-state conditions, at the same Reynolds number and angle of attack (Re=136 and =40°), the calculated CL is 0.59. Vogel (1967
) measured the forces on a thin plate cut to the shape of a fruit fly wing in a wind tunnel, at Re=200 and
=40°, CL was approximately 0.6. In fact, from our steady-state calculations and from the measurements of Vogel (1967
), the maximum value of CL at steady-state conditions is approximately 0.6. For reference, we used this maximum CL value and the translational velocity of the wing, ut (velocity at r0), to estimate the quasi-steady lift of the wing, 0.6(0.5
ut2S), and therefore the quasi-steady lift coefficient, 0.6(ut2/U2). The quasi-steady lift coefficient is plotted in Fig. 3B for comparison. For the flapping motion described above, CL is much larger than the corresponding quasi-steady value for most of the stroke. The mean value of CL is 1.2, more than twice the quasi-steady value (0.5 in Fig. 3). As seen in Fig. 3, this large mean CL results from three causes: the CL peak at the beginning of the stroke, the large CL in the middle of the stroke and the large CL peak near the end of the stroke. The above analysis showed that the CL peak at the beginning of the stroke could be explained by the rapid acceleration of the wing, the large CL in the middle of the stroke by the stall-absence mechanism and the large CL peak near the end of the stroke by the rapid pitching-up rotation of the wing while in constant-speed translation. However, Dickinson et al. (1999
) gave somewhat different explanations for the CL peaks at the beginning and near the end of the stroke. These explanations are investigated in greater detail below.
Effects of acceleration at the beginning of a stroke and effects of the wake of the previous strokes
Fig. 6 presents results for when t is varied in equation 1 (smaller
t corresponds to a larger acceleration at the beginning of a stroke); other conditions are the same as in case A. A larger acceleration results in a larger CL peak at the beginning of the stroke, indicating that the CL peak at the beginning of a stroke is closely related to the rapid acceleration of the wing. Dickinson et al. (1999
) considered that the large CL peak at the beginning of a stroke could be explained by the mechanism of wake capture. This mechanism was revealed by Dickinson (1994
) for a two-dimensional aerofoil in a simplified flapping motion with two opposite strokes, in which the flow generated by the first stroke could increase the effective fluid velocity, and hence the lift, at the start of the next stroke. In a recent study by Sun and Hamdani (2001
), it was shown that, for a two-dimensional aerofoil, if the first stroke was very short and the dynamic stall vortex did not shed (or a vortex street similar to the von Karman vortex street did not form), the mechanism of wake capture would not exist. For a three-dimensional wing in flapping motion, Ellington et al. (1996
) showed (see also above) that the dynamic stall vortex does not shed and, therefore, that the mechanism of wake capture might not exist.
|
|
|
Effects of varying the location of the wing-rotation axis and the rotation rate
The large CL peak near the end of a stroke was explained by Dickinson et al. (1999) as analogous to the Magnus effect. When a moving cylinder or sphere spins, the friction between the fluid and the surface of the body tends to drag the fluid near the surface in the same direction as the rotational motion. Superimposed onto the usual non-spinning flow, this extra velocity contribution creates a higher-than-usual velocity (hence lower pressure) at the top of the body and a lower-than-usual velocity (hence higher pressure) at the bottom, resulting in lift. This phenomenon is called the Magnus effect. The Magnus effect is a steady-state flow phenomenon. Note that it cannot explain the large CD peak that accompanies the large CL peak (see Figs 3 and 6) (see also fig. 1D in Dickinson et al., 1999
). It was shown above that the large CL and CD peaks at the end of the stroke were due to the effects of rapid pitching-up rotation. This effect would be expected to become weaker when the location of the rotation axis is moved rearward or when the rotation rate is lower. Fig. 9 gives the results for three different locations of the rotation axis (including that in case A, where it was located at 0.2c from the leading edge of the wing); other conditions are the same as in case A. When the rotation axis is moved rearward, the peak CL (and CD) decreases greatly.
|
|
|
Fig. 12 gives examples of the vorticity plots for symmetrical rotation. At the beginning of the stroke (Fig. 12A,B), the wing performs the same translational acceleration as in case A (from =16.7 to
=18, towards the right). But since it also rotates clockwise (from
=86.5° to
=140°) and the rotation axis is near the leading edge of the wing (0.2c from the leading edge), a large part of the wing moves backwards (towards the left). As a result, the starting vortex (Fig. 12B) has less vorticity and moves less far downstream from the wing than in case A at the corresponding time (Fig. 4C). Sun and Hamdani (2001
) and Lan and Sun (2001
) showed that, if the starting vortex is weaker and moves less far downstream in a given time, less lift will be generated. This helps to explain the smaller CL at the beginning of the stroke. Near the end of the stroke, e.g. at
=20.9 (Fig. 12D), the wing has just started pitching-up rotation and a small new vorticity is generated [at the corresponding time in case A (see Fig. 4G) a strong new vorticity layer is generated]. This helps to explain the smaller CL peak near the end of the stroke.
|
|
Comparison with model-wing experiment
Dickinson et al. (1999) measured the unsteady lift of a robotic wing modelled on the fruit fly. Their lift is presented in dimensional form and, for comparison, we need to convert it into the lift coefficient. In their experiment, the fluid density
was 0.88x103 kg m3, the wing length l was 0.25 m and the wing area S was 0.0167 m2. The translational speed at the wing tip during the constant-speed translation phase of a stroke, given in fig. 3D of Dickinson et al. (1999
), is 0.235 m s1 and, therefore, the reference speed (speed at r0=0.58l) is calculated as U=0.14 m s1. From the above data,
. Using the definition of CL (equation 23), the lift in fig. 3A of Dickinson et al. (1999
) can be converted to CL (Fig. 14). The CL curve from Fig. 6B with
t=1.14 (approximately the same as that used by Dickinson et al., 1999
) is plotted for comparison. The calculated CL is clearly smaller than the measured values, and they differ by almost the same amount throughout the majority of the stroke. The model wing used in the present calculations had an aspect ratio of 2.76, which is much smaller than that used by Dickinson et al. (1999
), which can be calculated as
. This partly explain the lower CL. Despite the quantitative differences between the calculated and experimental results, the overall pattern is very similar.
|
We investigated whether the insect weight can be balanced by the mean lift calculated in the present study. From the definition of CL in equation 23:
![]() | (25) |
where is the mean lift and
is the mean lift coefficient (air density
is taken as 1.25x106 g cm4 s2). The speed at the radial location r0 was defined as the translational speed of the wing, and the translational speed during the constant-speed translation phase of a stroke was taken as the reference speed U. To calculate of the results shown in Figs 3 and 11, the mean translational speed over a stroke was 0.825U. For the hovering fruit fly of Weis-Fogh (1973
), the mean translational speed is 2
nr0=219 cm s1, and U can be calculated as U=219/0.825=265 cm s1. Inserting the values of
, U and St into equation 25, the mean lift can be written as
. For
to equal the insect weight,
should therefore be 1.96x105/(2.50x105)=0.78.
For the symmetrical rotation and rotation-advanced cases, was 0.91 and 1.20, respectively. Both exceed that needed to balance the weight. It should be noted, however, that the angle of attack at midstroke used in these calculations was 40°, which may be larger than the actual value. Calculations were performed with smaller angles of attack at midstroke (other conditions kept the same). It was shown that, when the angle of attack was 36°,
was 0.79 for symmetrical rotation and 1.02 for the rotation-advanced case; when the angle of attack was 34°,
was 0.75 and 0.98, respectively.
The above results show that, with symmetrical rotation and an angle of attack at midstroke of approximately 36°, the fruit fly can produce enough lift for hovering flight. Ellington (1984b) observed many small insects in hovering flight, including the fruit fly, and found an angle of attack of approximately 35°, so our predicted value is in very good agreement with the observations. The above results also show that, if the insect employs a larger angle of attack or changes the timing of wing rotation, much greater lift can be produced for manoeuvring or other purposes.
![]() |
List of symbols |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
c mean chord length
CD drag coefficient
CL lift coefficient
D drag
J Jacobian of the transformation between (x,y,z,t) and (,
,
,
)
l wing length
L lift
mean lift
n wingbeat frequency
o origin of the inertial frame of reference
o' origin of the non-inertial frame of reference
p non-dimensional fluid static pressure
r radial position along the wing length
r0 radius of the second moment of wing area
Re Reynolds number
S area of one wing
St area of a wing pair
t* time
u,v,w non-dimensional velocity component in the x,y,z directions, respectively
ut translational velocity of the wing
ut+ non-dimensional translational velocity of the wing
U reference velocity
x,y,z coordinates in the inertial frame of reference
x',y',z' coordinates in the non-inertial frame of reference
angle of attack
' 180°
angular velocity of pitching rotation
+ non-dimensional angular velocity of pitching rotation
r duration of pitching rotation, non-dimensional
t duration of translational acceleration, non-dimensional
kinematic viscosity
,
,
transformed coordinates
density of fluid
non-dimensional time (
=t)
1 time when translational deceleration starts, non-dimensional
r time when pitching rotation starts, non-dimensional
stroke amplitude
azimuthal rotation angle
angular velocity of azimuthal rotation
z' spanwise component of vorticity, non-dimensional
gradient operator (
/
x,
/
y,
/
x)
![]() |
Acknowledgments |
---|
![]() |
References |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Dickinson, M. H. (1994). The effects of wing rotation on unsteady aerodynamic performance at low Reynolds numbers. J. Exp. Biol. 192, 179206.
Dickinson, M. H. and Götz, K. G. (1993). Unsteady aerodynamic performance of model wings at low Reynolds numbers. J. Exp. Biol. 174, 4564.
Dickinson, M. H., Lehman, F. O. and Sane, S. P. (1999). Wing rotation and the aerodynamic basis of insect flight. Science 284, 19541960.
Ellington, C. P. (1984a). The aerodynamics of hovering insect flight. I. The quasi-steady analysis. Phil. Trans. R. Soc. Lond. B 305, 115.
Ellington, C. P. (1984b). The aerodynamics of hovering insect flight. III. Kinematics. Phil. Trans. R. Soc. Lond. B 305, 4178.
Ellington, C. P., van den Berg, C. and Willmott, A. P. (1996). Leading edge vortices in insect flight. Nature 384, 626630.
Hilgenstock, A. (1988). A fast method for the elliptic generation of three dimensional grid with full boundary control. In Numerical Grid Generation in CFM88 (ed. S. Sengupta, J. Hauser, P. R. Eiseman and J. F. Thompson), pp. 137146. Swansea UK: Pineridge Press Ltd.
Lan, S. L. and Sun, M. (2001). Aerodynamic properties of a wing performing unsteady rotational motions at low Reynolds number. Acta Mech. 149, 135147.
Liu, H., Ellington, C. P., Kawachi, K., Van Den Berg, C. and Willmott, A. P. (1998). A computational fluid dynamic study of hawkmoth hovering. J. Exp. Biol. 201, 461477.
Rogers, S. E. and Kwak, D. (1990). Upwind differencing scheme for the time-accurate incompressible NavierStokes equations. AIAA J. 28, 253262.
Rogers, S. E., Kwak, D. and Kiris, C. (1991). Numerical solution of the incompressible NavierStokes equations for steady-state and dependent problems. AIAA J. 29, 603610.
Spedding, G. R. (1992). The aerodynamics of flight. In Advances in Comparative and Enviromental Physiology, vol. II, Mechanics of Animal Locomotion (ed. R. McN. Alexander), pp. 51111. London: Springer-Verlag.
Sun, M. and Hamdani, H. (2001). High-lift generation by an airfoil performing unsteady motion at low Reynolds number. Acta Mech. Sinica 17, 97144.
Vogel, S. (1967). Flight in Drosophila. III. Aerodynamic characteristics of fly-wing and wing models. J. Exp. Biol. 46, 431443.
Weis-Fogh, T. (1973). Quick estimates of flight fitness in hovering animals, including novel mechanism for lift production. J. Exp. Biol. 59, 169230.