## SUMMARY

A finite element flow solver was employed to compute unsteady flow past a three-dimensional *Drosophila* wing undergoing flapping motion. The computed thrust and drag forces agreed well with results from a previous experimental study. A grid-refinement study was performed to validate the computational results, and a grid-independent solution was achieved. The effect of phasing between the translational and rotational motions was studied by varying the rotational motion prior to the stroke reversal. It was observed that, when the wing rotation is advanced with respect to the stroke reversal,the peak in the thrust forces is higher than when the wing rotation is in phase with the stroke reversal and that the peak thrust is reduced further when the wing rotation is delayed. As suggested by previous authors, we observe that the rotational mechanism is important and that the combined translational and rotational mechanisms are necessary to describe accurately the force time histories and unsteady aerodynamics of flapping wings.

## Introduction

Flapping foil propulsion has received considerable attention in the past few years as an alternative to the propeller. This mode of propulsion, which involves no body undulation, has many applications, such as propulsion of submersibles, maneuvering and flow control, that are of interest to the hydrodynamic community and unconventional aerodynamics of micro aerial vehicles (MAVs) and the study of aircraft flutter that are of interest to the aerodynamic community.

Flapping foil propulsion is also important in the area of biofluid dynamics for the study of propulsion in insects, birds and certain aquatic animals. Flying animals generate lift and thrust as a consequence of the interaction between the flapping motions of the wings and the surrounding air. These animals also perform maneuvers involving rapid plunging and pitching motions. Conventional steady-state theories do not predict sufficient forces to meet those required for flight (Ellington,1984). Therefore, we need to understand the unsteady aerodynamics of flapping wings undergoing highly three-dimensional motions with widely varying geometries.

Experimental work on two-dimensional flapping foils has been carried out by Anderson (1996) and Freymuth(1999). Computational studies have been performed by Jones and Platzer(1997) and Ramamurti and Sandberg (2001). While two-dimensional wing section investigations can yield useful insights on the coupled pitching and heaving dynamics, nothing can be learned concerning the influence of spanwise flow. It is therefore essential to carry out computations for actual three-dimensional insect wings. Ramamurti et al.(1996) computed three-dimensional unsteady flow past moving and deforming geometries in a simulation of a swimming tuna with caudal fin oscillation.

The three-dimensional wing strokes of insects can be divided into two translational phases and two rotational phases. During the translational phases, the upstroke and the downstroke, the wings move through the air with high angles of attack, and during the rotational phases, pronation and supination, the wings rotate rapidly and reverse direction. Dickinson et al.(1999) studied the effects of translational and rotational mechanisms of the wing in *Drosophila melanogaster*. They directly measured the forces produced by flapping wings and explained the aerodynamics of insect flight by interactions between three unsteady flow mechanisms. The `delayed' stall mechanism is a translational mechanism which, in two dimensions, produces high lift in the initial phases of translation until eventual flow separation; in three dimensions, the spanwise flow effectively prevents stall. Rotational circulation and wake capture are rotational mechanisms that depend mainly on the pronation and supination of the wing during stroke reversal. Walker and Westneat (1997) have studied experimentally the kinematics of fin motion in a fish, the bird wrasse *Gomphosus varius*. Liu and Kawachi(1998) numerically investigated the flow over a hovering hawkmoth, *Manduca sexta*. They reported the presence of a spiral leadingedge vortex to which they attributed the lift enhancement mechanism. They validated their results by comparing the computational streak lines found in a two-dimensional hovering airfoil with the experimental smoke visualization, but they did not directly compare instantaneous forces.

Here, we extend the two-dimensional pitching and heaving airfoil computations to three dimensions. This study will address the role of the rotational motion in detail. Also, the role of the leading edge vortex and the interaction between the axial flow and this leading edge vortex are investigated. The primary objectives were (i) to validate the three-dimensional computations by comparing the forces with the experimental results of Dickinson et al.(1999) and performing gridrefinement studies, to verify the hypothesis of Dickinson et al.(1999) that rotational mechanisms of the wing form the basis by which the insect modulates the magnitude and direction of the forces during flight and (ii) to provide data on the forces, moments and power required for the development of a robotic fly.

To this end, computations are performed for various phase angles between the rotation and translation motions, and the time history of the unsteady forces is compared with the experimental results. The flow solver we employ is a finiteelement-based incompressible flow solver based on simple, low-order elements. The simple elements enable the flow solver to be as fast as possible, reducing the overheads in building element matrices, residual vectors, etc. The governing equations are written in Arbitrary Lagrangian Eulerian form, which enables flow with moving bodies to be simulated. The details of the flow solver, the rigid body motion and adaptive remeshing are given by Ramamurti et al.(1995) and are summarized below.

## Materials and methods

### The incompressible flow solver

*p*denotes pressure, and

**v**-

_{a}=v*w*, the advective velocity vector, where

**v**is flow velocity and

**w**is mesh velocity and the material derivative is with respect to the mesh velocity

**w**. Both the pressure

*p*and the stress tensor σ have been normalized by the (constant) density ρ and are discretized in time

*t*using an implicit time-stepping procedure. It is important for the flow solver to be able to capture the unsteadiness of a flow field. The present flow solver is built as time-accurate from the onset, allowing local time stepping as an option. The resulting expressions are subsequently discretized in space using a Galerkin procedure with linear tetrahedral elements. To be as fast as possible, the overheads in building element matrices, residual vectors, etc., should be kept to a minimum. This requirement is met by employing simple, low-order elements that have all the variables (

**v**,

*p*) at the same location. The resulting matrix systems are solved iteratively using a preconditioned conjugate gradient algorithm (PCG), as described by Martin and Löhner(1992). The flow solver has been successfully evaluated for both two-dimensional and three-dimensional laminar and turbulent flow problems by Ramamurti and Löhner(1992) and Ramamurti et al.(1994).

To carry out computations of the flow about oscillating and deforming geometries, one needs to describe grid motion on a moving surface, i.e. to couple the moving surface grid to the volume grid. The volume grid in the proximity of the moving surface is then remeshed to eliminate badly distorted elements. The velocity of the mesh is obtained in a manner so as to reduce this distortion. A detailed description of the various mesh movement algorithms is given in Ramamurti et al.(1994). In that study,smoothing of the coordinates was employed for the mesh movement with a specified number of layers of elements that move rigidly with the wing. In two-dimensional studies (Ramamurti and Sandberg, 2001), the grid showed that the elements at the edge of the rigid layers were very distorted after one cycle of oscillation. This is due to a residual mesh velocity that is present as a result of the non-convergence of the mesh velocity field. This will appear whether a spring-analogy or a Laplacian-based smoothing is used.

^{n+1}were obtained as a weighted average of the original grid point location at time

*t*=0(x

^{0}) and the location of the point as if it moved rigidly with the body(x

_{rigid}

^{n+1}):

*f*) is a simple linear function based on the distance from the center of rotation

*r*, and is given by:

The values of *r*_{min} and *r*_{max} used in this study are 2.0 and 10.0, respectively. To reduce the computational effort,the experimental arrangement of Dickinson et al.(1999) was approximated by introducing a symmetry plane. Because of the proximity of the wing at the beginning of the downstroke and the rotation of the wing during the pronation phase, the normal component mesh velocity of the points on the symmetry plane can become non-zero. This would result in the points being pulled away from the symmetry plane. To avoid this problem, the points on the symmetry plane are allowed to glide along this plane. A similar technique has been employed for the simulation of torpedo launch from a submarine by Ramamurti et al.(1995) where the gap between the launch tube and the torpedo was small.

## Results and discussion

*Drosophila melanogaster*is shown in Fig. 1A. The coordinate system (

*x,y,z*) is fixed to the body with the

*x*coordinate normal to the stroke plane. During the translational phases (upstroke and downstroke), the wing moves from close to the

*y*axis through an angleϕ, the wingbeat amplitude. The flapping wing configuration used in the flow simulations is shown in Fig. 1B. This is based on the experimental arrangement of Dickinson et al. (1999). Fig. 1B shows the position of the wing at three different times during the flapping cycle. The coordinate system (

*x′,y′,z′*) is fixed to the wing, and the wing rotates about the

*z′*axis throughout the cycle. The planform of the wing is derived from the

*Drosophila*wing and is 25 cm long and 3.2 mm thick. The experimental apparatus consisted of two wings immersed in a tank of mineral oil. The viscosity of the oil, the length of the wing and the frequency of the flapping motion were chosen to match the Reynolds number (

*Re*) of a typical

*Drosophila melanogaster*,approximately 136. The

*Re*for the present calculations is defined on the basis of the mean chord of the wing

*U*

_{t}(ignoring the forward velocity), as follows:

*U*

_{t}=2ϕ

*nR, R*is the wing length (25 cm),

*AR*is the aspect ratio of the wing,

*n*is the frequency of flapping motion, ϕ is the wingbeat amplitude (peak to peak, in rad) and v is the kinematic viscosity (115 cSt=115×10

^{-6}m

^{2}s

^{-1}).

The kinematics of the wing motion is obtained from the experiments of Dickinson et al. (1999). The angles of rotation about the *x* axis (roll) and the *z*′axis (pitch) are shown in Fig. 2A. The motion of the wing is prescribed using these two angles. Fig. 2B shows the translational velocity of the wing tip and the rotational (angular) velocity of the wing. Three different phases between the translational and rotational motions were used. In the `advanced' case, wing rotation precedes stroke reversal by 8 % of the wingbeat cycle; `symmetrical' wing motion is where the wing rotation occurs symmetrically with respect to stroke reversal; in `delayed' wing motion, rotation is delayed by 8 % with respect to stroke reversal. Wingbeat amplitude is 160°, flapping frequency is 0.145 Hz and the angle of attack at midstroke is approximately 40° during both upstroke and downstroke.

### Symmetrical case

The flow solver described here is employed to compute the flow past the *Drosophila* wing undergoing translation and rotation. First, an inviscid solution was obtained using a grid consisting of 178 219 points and 965 877 tetrahedral elements. An initial steady-state solution was obtained in 1500 time steps. The unsteady solution using the prescribed kinematics(Fig. 2) is then obtained. The surface pressure on the wing is integrated to obtain the forces on the wing along the three axes (*F*_{x}, *F*_{y}, *F*_{z}). The thrust *T* and the drag *D* forces are then computed as *T*=-*F*_{x} and *D*=√(*F*_{y}^{2}+*F*_{z}^{2}),respectively. These forces are compared with those obtained from the experiments of Dickinson et al.(1999).

*r*

_{2}

^{2}(

*s*) is the second moment of the dimensionless area of the wing (0.40). The variation of these forces during the translational phase of the wing is also predicted correctly,but the magnitude of the thrust force during the downstroke is higher than that of Dickinson et al.(1999). The kinematics is symmetrical between the up- and downstroke, so the resultant force should also be symmetrical. That this is not the case may be due the mechanical play in the experimental arrangement, as suggested by M. H. Dickinson (personal communication). To understand the different mechanisms occurring during the wingbeat cycle, we can divide the cycle into two rotational and two translational phases. The rotational phase near the beginning of the downstroke (pronation) occurs between time

*t*0 and

*t*3(Fig. 3A). Thrust decreases between

*t*0 and

*t*1, then increases until

*t*2. This behavior can be explained by a rotational mechanism. The wing continuously rotates in a counterclockwise direction producing a circulation pointing nearly along the +

*y*direction. Between

*t*0 and

*t*1,the wing is translating in the -

*z*direction, resulting in a force pointing in the -

*x*direction, thus producing a peak in thrust at

*t*0. If a rotational mechanism alone were present, the thrust should continue to decrease until

*t*3; in fact, the thrust force increases between

*t*1 and

*t*2. This happens after the wing changes direction at the start of each half-stroke. Dickinson et al.(1999) attribute this increase in thrust to a wake-capture mechanism in which the wing passes through the shed vorticity of the previous stroke.

The position of the wing at the beginning of the downstroke is shown in Fig. 4A. The chord in this case of symmetrical rotation is aligned with the *x* axis at this instant. We found a separation bubble attached to the leading edge during the interval *t*0-*t*1. Particles were released from a rake of rectangular grid of points in a plane 0.8 mm away from the bottom surface of the wing. Using the instantaneous velocity field, the positions of these particles were obtained by integrating the velocity at these rake points in both the positive and negative velocity directions until the length of the traces exceeded a specified length or the particles ended on a solid boundary or exited the computational domain. These instantaneous particle traces are colored according to the magnitude of velocity at that location. Fig. 5 shows the leading edge vortex with vorticity oriented in the +*y* direction. This leading edge vortex was created at the end of the preceding upstroke. This vortex is located below the wing and is rotating in the counterclockwise direction,which can be seen from the velocity vectors shown in Fig. 6A. A possible explanation for the increase in thrust between *t*1 and *t*2 is that the wing moving through this wake benefits from the shed vorticity. As the wing moves through this vortex during the downstroke, it produces a stagnation region at the bottom of the wing near the *z*′ axis(Fig. 5), resulting in an increase in thrust.

At the beginning of the downstroke (*t*1), the flow separates at the leading edge and reattaches on the bottom of the wing as shown in Fig. 6A. As the wing continues to move down, between *t*1 and *t*2, the separation point of this bubble moves back along the wing chord(Fig. 6). During the interval *t*2-*t*3, we observe a trailing edge separation bubble forming(see Fig. 7A-C). A similar separation region forms at the wing tip, as can be seen in Fig. 7D-F. Fig. 4B shows the position of the wing at *t*=12.5 s. Thrust production reaches a local minimum around this instant (Fig. 3). In Fig. 7A,D, a large recirculation region can be seen in the wake of the wing. This separated flow from the wing tip and trailing edge will result in a higher pressure on the upper surface of the wing and, hence, a reduction in the thrust. During the interval *t*1-*t*3, the magnitude of the translational acceleration of the wing decreases while that of the angular acceleration increases (symmetrical phase, Fig. 8).

Between *t*2 and *t*3, the magnitude of the translational acceleration is large enough to overcome the rotational effect and, when the angular acceleration becomes large enough, the rotational mechanism dominates,resulting in the observed reduction in thrust until *t*3. Between *t*3 and *t*4, the translational effect should result in a constant thrust because the translational acceleration is almost constant during this period. The rotational effect produces an increase in thrust between *t*3 and *t*4, with a plateau in the middle, which occurs when the trailing edge vortex is shed. Similar trends are observed during the supination phase prior to the beginning of the upstroke(*t*4-*t*5) and at the beginning of the upstroke(*t*5-*t*7).

Fig. 9 shows the instantaneous traces or streaklines of particles released 3.0 cm above or 3.5 cm below the wing in a plane parallel to the wing and near the leading edge. In Fig. 9A, a wing tip vortex can be seen, but no leading edge vortex is visible above the wing surface. A leading edge vortex spinning in the counterclockwise direction is found below the wing surface (Fig. 9B). Particle traces near the mid stroke are shown in Fig. 10A. At this instant, the wing rotation axis *z*′ is aligned with the body coordinate *z* (see Fig. 1B). Here,we can see the beginnings of the leading edge vortex on the upper surface of the wing that is also shown by the velocity vectors in the *xy* plane at *z*=20 cm (Fig. 10B). Fig. 10Cshows the pressure contours on the upper surface of the wing. The pressure is non-dimensionalized with respect to the dynamic head,½ρ *U*_{t}^{2}, where ρ is the density of the mineral oil (880 kg m^{-3}) used in the experiments of Dickinson et al. (1999). A region of constant pressure is present from the root of the wing up to approximately 60%of the span and extending to the trailing edge. Fig. 11 shows the particle traces prior to the end of the downstroke. Here, the leading edge vortex is clearly seen on the upper surface of the wing.

### Grid-refinement study

To assess the sensitivity of our computational results, we carried out a grid-refinement study. The resolution of the grid in the vicinity of the wing was doubled. The computations were carried out using a grid consisting of approximately 238×10^{3} points and 1.3×10^{6}tetrahedral elements. The time step was halved for this computation. The computed thrust forces are shown in Fig. 12. It can be seen that the agreement between the two analyses is very good; even the coarse grid produces adequate resolution.

### Viscous effect

To the study the effects of viscosity, a laminar viscous computation was carried out, for *Re*=120. Because of the lack of information on the transition to and the presence of uncertainties in turbulence modeling,laminar flow was assumed first. Fig. 13 shows the time history of the thrust and drag forces for the inviscid and the viscous cases. The finer mesh employed in the grid-refinement study (see above) was used for this computation, and the mesh size near the leading and trailing edge of the wing was approximately 0.02. It is clear that viscous effects are minimal, and the thrust and drag forces are dominated by translational and rotational mechanisms. Hence, inviscid computations were carried out to study the effect of phasing between the translational and rotational motions.

### Effect of phasing

Fig. 14 compares the forces produced when the rotational motion precedes stroke reversal (`advanced' case)with those of Dickinson et al.(1999). Again, the agreement with the experimental results is good. In this case, the peak in the thrust force is achieved prior to the beginning of the downstroke at *t*=10.76 s and is approximately 0.56N compared with a value of 0.47N for the symmetrical case (see Fig. 3A). This can be explained by the rotational mechanisms discussed above. The rotational effect diminishes prior to the beginning of the downstroke,producing a negative thrust of 0.2N. The thrust then increases until *t*=11.98s. During this period, the wing moves through the wake created during the upstroke, as in the symmetrical case, resulting in a high pressure on the bottom of the wing. The velocity vectors near the leading edge are shown in Fig. 15. During this period, both the translational and rotational accelerations are in phase(Fig. 8). The peak thrust is approximately 0.48N compared with a value of 0.28N for the symmetrical case. Thereafter, the combined effect of rotational and translational motions produces reduced thrust until a second peak arises due to the rotational motion at *t*=14.23s, prior to the beginning of the upstroke. The mean thrust force for one wingbeat cycle is approximately 0.312N, and the mean thrust coefficient

In the `delayed' case, where wing rotation is delayed with respect to stroke reversal, the rotational motion does not produce any thrust prior to the beginning of the downstroke (Fig. 16A). The mean thrust force for one wingbeat cycle is approximately 0.206N and the mean thrust coefficient

*t*=12.05 s. In this case, the high pressure on the bottom of the wing together with the orientation of the wing cause a reduction in thrust. Subsequently, the combined translational and rotational mechanisms result in an increase in thrust. At

*t*=12.8 s,we observe a plateau region in the thrust(Fig. 16A). During this period, the presence of a trailing edge vortex on the upper surface(Fig. 17C,D) increases the pressure on the upper surface of the wing, resulting in a temporary loss of thrust; when this vortex leaves the trailing edge, thrust increases again.

Fig. 18 shows the magnitude of velocity in the *xz* plane at y=10 cm at the beginning of the downstroke for the three cases of wing motion in the wake created by the wing. Velocities are greatest for the advanced case and smallest for delayed rotation. The wing moving through the higher-velocity fluid therefore produces an additional thrust in the advanced rotation case, whereas the wing for the delayed case intercepts the flow at an angle that produces negative thrust. Similar velocity fields can be seen in the particle image velocimetry data of Dickinson et al. (1999).

### Mechanical aspects of the flapping wing

The results of the present study were then used to derive the input forces,moments, power requirement and efficiency for the creation of a robotic fly. First, the forces on the wing were integrated to a particular spanwise location from the root. The forces were obtained by marking this location and then computing the force contribution of all the surface elements on the wing up to this location. This location was then tracked using the prescribed rigid body motion. The elements contributing to the force up to this location are recomputed as the surface mesh is regenerated due to the motion. Three different spanwise locations were chosen: quarter span, half span and three-quarter span of the wing. Fig. 19 shows their force contributions compared with the total contribution of the wing and half the total thrust force generated by the wing for the symmetrical rotation case. It is clear from this figure that nearly half the thrust is generated by the outer 25% of the wing.

Moments about the wing root in the wing coordinate system(*x′,y′,z′*) were then computed from the moments about the fixed body coordinates and the prescribed kinematics of the wing(Fig. 20). The moment about the wing rotation axis *z′* (*M*_{z′}) is nearly zero throughout the cycle, implying that there is no torsional load on the system. The variation of moment about the *x′* (=*x*)axis (*M*_{x′}) is anti-symmetrical because only one wing is considered for the moment computation.

*P*

_{in}is computed by:

**F**is the force vector and W

_{wing}is the velocity of the surface of the wing. Fig. 21shows the instantaneous power requirement for one wing for the three phases of wing rotation. The mean power required is 0.024 W for the symmetrical case,0.039 W for the advanced case and 0.024 W for the delayed case. The mean thrust for the symmetrical case is 0.318 N; values for the advanced and delayed cases are 0.312 N and 0.206 N, respectively.

## Acknowledgements

This work was supported by the Office of Naval Research through the Tactical Electronic Warfare Division Micro Air Vehicles Program of the Naval Research Laboratory. The authors would like to thank Professor Michael Dickinson and Mr Sanjay Sane for providing the experimental kinematics and data for comparison and for very useful discussions and Professor Rainald Löhner of the George Mason University for his support throughout the course of this work. The computations carried out for this work were supported in part by a grant of HPC time from the DoD HPC centers, ARL MSRC SGI-O2K and NRL SGI-O2K.

## References

**Anderson, J. M.**(

**Dickinson, M. H., Lehmann, F.-O. and Sane, S.**(

**Ellington, C. P.**(

**Freymuth, P.**(

**Jones, K. D. and Platzer, M. F.**(

**Liu, H. and Kawachi, K.**(

**Martin, D. and Löhner, R.**(

**Meirovitch, L.**(

**Ramamurti, R. and Löhner, R.**(

**Ramamurti, R., Löhner, R. and Sandberg, W. C.**(

**Ramamurti, R., Löhner, R. and Sandberg, W. C.**(

**Ramamurti, R. and Sandberg, W. C.**(

**Ramamurti, R., Sandberg, W. C. and Löhner, R.**(

**Walker, J. A. and Westneat, M. W.**(

*Gomphosus varius*(Labridae).