A study on forces acting on a flapping wing

In order to study the forces acting on a flapping wing, an experimental investigation is performed in steady water flow. In this study, a SD7003 airfoil undergoes combined pitching and plunging motion which simulates the forward flight of small birds. The frequency of pitching motion is equal to the frequency of plunging motion and pitch leads the plunge by a phase angle of 90 degrees. The experiments are conducted at Reynolds numbers of 2500 ≤ Re ≤ 13700 and the vortex formation is recorded using the digital particle image velocimetry (DPIV) technique. A prediction of thrust force and efficiency is calculated from the average wake deficit of DPIV data, the near-wake vorticity patterns and time dependent velocity vectors are determined to comment on the thrust and drag indication. Direct force measurements are attempted using a Force/Torque sensor which is capable of measuring forces and moments in three axial directions.


Introduction
Mankind has wondered about flying animals since the early ages.The natural way of producing lift and propulsive force in flight is flapping and the wings let animals to produce thrust force in low Reynolds number range.Animals and insects which generate thrust in water and air has been researched for the last 10-15 years to understand the principles of the flight of these animals [1,2].Biomimetic research assists to understand the fundamental science of flight.Flapping wing aerodynamics may improve the design and flight performance of Micro Air Vehicles (MAVs) having application areas of targeting, surveillance, reconnaissance, scouting, biochemical sensing and remote observation of hazardous environments, in addition to their military use.
The wing structure and flight characteristics of insects are different from birds.Furthermore, big birds which are able to glide and soar have different wing shape deformation and flight characteristics.Large animals mostly employ pure heaving or plunging motion in the low frequency regime at larger Reynolds numbers whereas small birds and insects mostly employ combined pitching and plunging motions in low Reynolds number regime.Therefore, flapping wing MAVs are mostly inspired by insects and small birds.
In low Reynolds number flow over a stationary airfoil, a regular Karman vortex street is generated where the upper row of vortices rotate clockwise and the lower row of vortices rotate counterclockwise.Since velocity distribution in the wake results in a distinct momentum deficit, indicative of drag, the regular Karman vortex street yields a drag producing wake [3].
In their earliest studies, Knoller [4] and Betz [5] realized and explained the thrust production of a plunging airfoil and Katzmayer [6] confirmed the theory by experiments.Karman and Burgers [7] commented that a thrust force could be produced by a wake of consisting of two rows of counter-rotating vortices on an airfoil in an incompressible flow.When the oscillation of the airfoil has high amplitude and frequency, the downstream velocity distribution becomes jet-like which is indicative of thrust on the airfoil [8].
Kramer's experiments [9] on a fixed wing with varied airflow incidence via rotatable vanes showed that increasing the angle of attack of the airfoil increases the maximum lift coefficient.Carr [10] indicated that vortex formed on the airfoil during unsteady motion is associated with extra lift.
Although there are many studies in literature which focus on pure pitching or pure plunging motions of airfoils, animals fly with a combined pitching and plunging motion in nature.However, since unsteady flows considering coupled pitching and plunging motions includes a large parameter space, combined motion of foils has been subject to relatively few investigations.The experimental investigation of Anderson et al. [11] showed that optimum propulsive efficiencies were obtained within an approximate range of Strouhal numbers 0.2-0.4.This range coincides with that presented by Triantafyllou et al. [12] for observed fish and cetaceans swimming at or near their maximum observed speed.However, Young and Lai [13] [14].Rival et al. [15] investigated the leading edge vortices experimentally using Particle Image Velocimetry (PIV) technique.They tested different airfoil kinematics including asymmetric and peak-shifted plunging motions as well as different motion kinematics according to the strength of LEV, timing of LEV shedding and its convection into the wake.Rank and Ramaprian [16] investigated the flow over a NACA 0015 airfoil in light stall conditions and obtained instantaneous velocity fields by using PIV to reveal the formation, growth and eventual shedding of a weak vortex.
The velocity distribution and its derivatives on the airfoil directly affect the loads acting on the airfoil [17,18].One of the major topics is the calculation of moments and forces acting on the airfoil.Over the last two decades, the investigations focus on force prediction on two dimensional circular cylinders, rigid, oscillating or in transitional motion [19][20][21][22][23] and airfoils or flat plates [24,25].The relation of the flow structures with the loading on the body is also investigated [24] for a sinusoidally pitching airfoil.
Estimating aerodynamic forces by using Particle Image Velocimetry (PIV) images with respect to wake excess velocity is also a popular topic in the flapping wing studies.Comparing pure plunging and pure pitching motions, Rival and Tropea [26] state that the variations in lift and moment are only a function of the shed vortex positioning relative to the airfoil itself.Ol et al. [27] studied the fluid physics of a pitching and plunging rigid airfoil in nominally two-dimensional conditions using experimental, computational and theoretical approaches which include direct force measurements and phaseaveraged 2D PIV in a water tunnel, 2D Reynoldsaveraged Navier-Stokes equations and Theodorsen's formula.
Rival et al. [15] emphasized experimental difficulties to measure forces since both aerodynamic and inertial forces affect the wing, and to obtain essentially high resolution images in a wide parameter range.Bernal et al. [28] conducted an experimental investigation to measure the unsteady forces acting on a pitching and plunging airfoil using Fiber Bragg Gratings (FBG) in a water channel.In their study, the lift, drag and pitching moment are measured around a SD7003 airfoil at Reynolds number of 60 000 under steady and unsteady flow conditions.Gursul and Ho [29] studied the relationship between freestream velocity and lift force.They measured the aerodynamic forces on a stationary NACA0012 wing in a water channel with sinusoidally varying freestream.They concluded that the lift increases when the streamlines reattach downstream of the LEV on the airfoil and the lift decreases when the vortex detaches from the leading edge.Hubel and Tropea [30] performed direct force measurements, parallel and perpendicular to the oncoming flow, with time averaged PIV measurements in the wake of a flapping airfoil.They state that the vertical force coefficient for flapping wing is higher than that of the fixed wing.
In order to study the forces acting on a flapping wing, an experimental investigation is performed in steady water flow.A SD7003 airfoil undergoes combined pitching and plunging motion which simulates the forward flight of small birds.In a previous study [31], investigated cases are classified into five flow structure categories based on instantaneous and averaged vorticity patterns and velocity fields around and in the near wake of the airfoil.In this study, according to the flapping parameters, the related leading and trailing edge vortex formations and their interactions are observed to comment on how they contribute to lift and propulsion.A prediction of thrust force and efficiency is calculated from the average wake deficit of DPIV data, the nearwake vorticity patterns and time dependent velocity vectors are determined to comment on the thrust or drag indication.Direct force measurements are attempted using a Force/Torque sensor which is capable of measuring forces and moments in three axial directions.The goal of the study is to correlate the predicted forces from PIV images for each flow structure category with direct measurement of aerodynamic forces.

Summary of Significant Parameters
The airfoil has two kinematic motions which are pitching and plunging.They are sinusoidal motions and each of them has a specific equation of motion as; where h(t) is the linear plunge motion, ψ is phase angle and h amp is the plunge amplitude.Additionally, α(t) is the angular pitch motion, α 0 is the initial angle of attack and α amp is the pitch amplitude.The oscillation frequency f is kept equal for pitching and plunging motions and pitch leads the plunge by a phase angle of 90 degrees.The displacement definitions are shown in figure 1.The reduced frequency (k) characterizes the oscillation behavior and is defined as; where f is the flapping frequency; c is chord length and U ∞ is the velocity of the freestream flow.
Another dimensionless parameter is Strouhal number (St) which characterizes the vortex shedding in the wake of flapping wing.Its definition is given as; where A is the peak-to-peak displacement of the airfoil trailing edge.The amplitude ratio of pitching and plunging motions is defined as; Due to the presence of plunging motion, it is necessary to consider the vertical velocity of the airfoil and consequently its effective angle of attack.The vertical velocity is the derivative of plunging motion given in equation ( 1); Therefore the effective angle of attack becomes; 7

Experimental Setup
The experiments are carried out in the close-circuit, large scale water channel in Trisonic Research Center at the Faculty of Aeronautics and Astronautics of Istanbul Technical University.The test section of the water channel is transparent Plexiglas and includes three segments with cross-sectional dimensions of 1010mm x 790mm.The range of freestream velocity of the water is between 25.5 mm/s and 137.5 mm/s which corresponds to Reynolds number range of 2500 ≤ Re ≤ 13700.SD7003 airfoil profile is optimized for low Reynolds number flows in the range considered for micro air vehicle operations [33,34].The transparent Plexiglas SD7003 wing is manufactured in CNC milling machine and has a chord length of 10 cm with a span of 30 cm.The airfoil model is mounted in a vertical cantilevered arrangement in the water channel about its quarter chord.The connection rod connects the airfoil to a servo motor which provides the sinusoidal pitching motion via a coupling system which itself is connected to a linear table with another servo motor for plunging motion.A rectangular Plexiglas endplate is suspended from the top of the water channel above the fully submerged airfoil in order to reduce the end and free-surface effects.The experimental setup can be seen in figure 2.

Force Measurements
An ATI Mini40 F/T sensor is attached to the vertical cantilevered arrangement between the wing and pitch servo motor.The applied sensor is capable of measuring forces and moments in three axial directions with sensing range of ±40N in F x and F y directions and ±120N in F z direction and ±2Nm in T x , T y and T z .

Image Acquisition and Post-Processing
A Digital Particle Image Velocimetry (DPIV) system is used for capturing and processing the instantaneous quantitative flow images.This method is a non-intrusive, optical laser flow measurement technique to acquire time dependent velocity fields.
Prior to the experiments, the water in the channel is seeded with silver coated hollow glass spheres with a mean particle diameter of 10μm.
DynamicStudio v2.21 software (Dantec Dynamics A/S) is employed for image acquisition and interrogation.The flow is illuminated by a dual cavity Nd:Yag laser with 15 Hz repetition for each cavity.The maximum output energy of the laser is 120 mJ per pulse.A beam expander is used to form and adjust the thickness of the laser sheet which horizontally illuminates the flow field and the vertically mounted airfoil from one side of the water channel.
Two 10-bit CCD cameras are positioned underneath the water channel for image acquisition.Two different setups for cameras are used in order to focus: -on only the near-wake of the airfoil with a full vertical extent of wake (Setup-1) -on the airfoil and in its near wake at the same time (Setup-2) The data for Setup-1 is obtained in a previous study [32] where the wing is confined in between two endplates, therefore represent 2D results.Image acquisition areas of two cameras for both setups can be seen in figure 3.In post processing, two images from the two cameras are stitched before interrogation by an in-house script run in the Matlab link of Dynamic Studio software.The cross-correlation technique is used to interrogate stitched images with a window size of 64x64 pixels, 50% overlapping in each direction.

Thrust Estimation from DPIV Data
An average thrust force and efficiency estimate were calculated from the DPIV velocity data using the wake excess velocity and the momentum theorem expressed for a control volume bounded by a control surface, using the steps provided by Anderson et al. [11].The control volume is sketched in figure 4 where the freestream velocity U ∞ is across the upstream boundary bc and the velocity u(y) is across the downstream boundary ad.Assuming incompressible, steady flow with constant pressure along the control volume boundaries, the thrust on the foil T can be written in terms of the momentum flux in the x-direction as; 8 Applying mass conservation, Equation 9 simplifies to:

9
Although full field DPIV velocity data is obtained to compute the time rate of change of momentum inside the control volume, the pressure along the control volume boundaries is not known and the assumption of steady flow is incorrect since the vortex wakes are unsteady with variable pressure.In spite of this, it is shown that the approximation using momentum flux calculation is a very good measure of the force, compared to actual force measurements for flapping foils and oscillating cylinders [1].
The thrust coefficient of the foil, C f , is defined similar to drag coefficient as follows: where c is the chord and s is the span of the airfoil.
The thrust coefficient for a flapping motion can be also calculated using the area swept by the airfoil.In that case, peak to peak value of the plunge motion should be used instead of the airfoil span.Therefore, the thrust coefficient based on the area swept by the wing, C T , will be in terms of thrust coefficient C f based on the planform area of the wing; The efficiency is the ratio of the work usefully applied to the total amount of work supplied.Using Prandtl's calculations [35], the theoretical efficiency for a propeller is defined as;

Results and Discussions
In the current study, the investigated cases are those representing five flow structure categories defined in the PhD Thesis of Fenercioglu [32].The non-dimensional values of the cases representing five flow categories with a sub-category are given in table 1.The flow structures around and in the near-wake of the airfoil are classified according to instantaneous vorticity patterns around the airfoil, instantaneous nearwake vorticity patterns and averaged flow fields.The mean angle of attack, the amplitude of pitching motion and the phase angle between pitch and plunge are fixed to be a 0 =8°, a amp =8.6° and ψ=-π/2 where pitch leads the plunge motion.Pitch and plunge frequencies, f, are kept equal to each other.

Flow structure categories
In Category A1, two counter-rotating pairs of vortices are formed from the leading and trailing edges during one oscillation cycle.One pair is on the lower side of the airfoil during the upstroke and the other pair is on the upper side during the downstroke.On the other hand, only one pair of counter-rotating vortices is shed in the near-wake during a cycle of oscillation.The timing of vortex shedding is such that the wake exhibits a jet-like velocity profile.
Category A2 is a subset of category A. In terms of formation of flow structures, it is very similar to the Category A1.The differences are due to its lower amplitude of plunging motion and higher frequency of oscillation.During one oscillation cycle, two stronger counter-rotating pairs of vortices are formed from the leading and trailing edges.The major difference of flow structures in comparison with the previous subset A1 is that the roll-up radius of the vortex shedding is small and the vortices are stronger.
Category B also presents two subsets with same distinguishable differences in between depending on the amplitude of plunge motion, similar to Category A. The presence of counter rotating pairs of vortices is still observed in Category B2, however, the vortices formed at the leading edge cannot have enough time to evolve as they are swept away during the following half cycle of the motion.While a negative vortex enlarges preserving its strength and detaches from the leading edge, a positive vortex is shed in the wake at the end of the downstroke.The process does not seem to be nearly reversed during the upstroke as in the case for Category A. In terms of loading, the presence of a reverse Karman street announces thrust production.
Category C is observed only for h amp /c = 0.25.The major characteristic of flow structures in this category is the shedding of two negative vortices per cycle of motion.During the upstroke, the formation of a negative vortex at the trailing edge is observed while another negative vortex is formed at the leading edge and moves towards the trailing edge.At the mid-plunge position they are about to be connected, however this does not happen; the one formed at the trailing edge is shed right after the mid-plunge position, while the other one is shed at the end of the upstroke.However, due to the different field of view (Setup-2) chosen in this study, this could not be observed during force measurement experiments.
In Category D, the roll-up of a positive vortex at the trailing edge is considerably reduced.The formation of counter-rotating pair of vortices during the upstroke as seen in the high Strouhal number cases of previous categories is no longer evident.The near-wake patterns indicate that the flapping motion still induces a jet-like velocity profile at this low Strouhal number regime.
In Category E, a negative vorticity layer on the upper surface and a positive vorticity layer on the lower surface of the airfoil are present during most of the stages and the wake also oscillates up and down with the airfoil motion.During the upstroke the vorticity layers are thin and attached to the surface; however, the negative vorticity layer detaches from the surface during the downstroke and exhibits localized cells of vorticity concentration, revealing Karman-like vortex formations in the nearwake at the beginning of the downstroke.The near-wake vorticity patterns indicate a velocity deficit.

Velocity deficit profiles from DPIV images
The velocity deficit profiles in the near-wake of the airfoil for all categories are given in table 2 which presents the results of drag/thrust coefficient and efficiency for each case along with the deviation from the freestream velocity in the near-wake, obtained at 1.1 chords away from the trailing edge.Only Category E does not exhibit thrust production and the resulting deviation from the freestream is so low that the velocity vectors are enlarged five times and colored in light blue to point out the difference.The jet-like flow is inclined in the cases designated with Category A. Therefore, a lift production along with thrust is expected.In accordance with Equation ( 12), high efficiencies are observed in general for low Strouhal numbers, just before a C d ≥ 0 value is observed.These cases are in general designated as Category D. 01028-p.5

DPIV and Direct force measurements
For each case investigated, 200 PIV images are acquired for 25 cycles of motion.The first image of each data set is taken at the beginning of the third period of motion.8 images of vorticity patterns representing a cycle of motion are analyzed for each category along with lift variation with respect to angle of attack and effective angle of attack in a period of motion and for all duration of force measurements.For the sake of brevity, only the DPIV results for Category A (A2) and Category E will be presented herein.
Figure 5 shows vorticiy patterns of an example case (U ∞ =0.03m/s, h amp /c=0.25,f=0.6Hz) representing Category A2.The first image belongs to the beginning of motion cycle where the airfoil is at its mid-plunge position moving up. Figure 6 presents an example case (U ∞ =0.1m/s, h amp /c=0.5, f=0.1Hz) for Category E. While strong LEVs are formed and shed to result a reverse Karman street for Category A2 case, an oscillating wake up and down with the airfoil motion is evident for Category E case.Although we are interested in thrust production, its existence can easily be identified from the wake surveys.On the other hand, lift variations are mostly studied (e.g.[27], [36], etc.) for oscillatory airfoil flows as lift can be calculated through quasi-steady approximation or modeled using harmonic motions in pitch and plunge as in the Theodorsen model [37].Indeed, thrust is associated with the timing of the LEV formation and separation which leads to the measurement and estimation of the maximum attainable lift in the dynamic-stall process.EFM 2012 Lift variation in a cycle for each case is given in table 3.Each DPIV image and vertical lines on the lift variations are numbered accordingly.The added mass term is subtracted from the measured directional forces, F x and F y which are defined on a reference frame fixed to the wing; F x is along the chord line and F y is perpendicular to the chord line.Lift is then calculated considering instantaneous geometric angle of attack due to pitching motion and instantaneous effective angle of attack including the effect of plunging motion.Table 3 also presents the variation of the effective angle of attack for a cycle of motion.The range of the vertical axes used in the graphs of table 3 for lift and effective angle of attack are fixed.The lift variations for all categories are observed to be sinusoidal in a general sense.Although the effective angle of attack variations are also sinusoidal and reach their maximum at the mid-cycle, maximum lift occurs at a different time in the cycle of oscillation.The case representing Category E exhibits very low amplitude of lift.For all other flow structure categories yielding thrust, there is no major difference to be noted in between the lift variations obtained with respect to geometric or effective angle of attack.However, the variations are more likely in terms of modulations as we move from Category A to D. According to quasi-steady approximation, one would expect to observe the maximum lift at the mid-cycle as it is in Categories B2 and C. The maximum lift is shifted by −π/2 with respect to the occurrence of the maximum effective angle of attack for Categories A2 and D.

Conclusion
In order to study the forces acting on a flapping wing, an experimental investigation is performed in steady water flow where a SD7003 airfoil undergoes combined pitching and plunging motion.Cases representing five flow structure categories are investigated using DPIV technique and direct force measurements.The observations are in accordance with previous studies.Different approaches including observation of flow structures, thrust estimation using velocity profiles in the near wake and direct force measurement yield similar major conclusions.Below a certain low Strouhal number, approximately St ≤ 0.2, the flapping motion does not result in thrust production and only the related flow structure category does not exhibit jet-like velocity profile in its near-wake.

Fig. 3 .
Fig. 3. Image areas of the cameras for two setups.

Table 1 .
The experiment categories.

Table 2 .
Velocity deficit profiles (u-U ∞ ) in the near wake and force coefficients.

Table 3 .
Lift and effective angle of attack variations for a cycle of motion.