Influence of Thickness Variation on the Flapping Performance of Symmetric NACA Airfoils in Plunging Motion

In order to investigate the impact of airfoil thickness on flapping performance, the unsteady flow fields of a family of airfoils from an NACA0002 airfoil to an NACA0020 airfoil in a pure plunging motion and a series of altered NACA0012 airfoils in a pure plunging motion were simulated using computational fluid dynamics techniques. The “class function/shape function transformation“ parametric method was employed to decide the coordinates of these altered NACA0012 airfoils. Under specified plunging kinematics, it is observed that the increase of an airfoil thickness can reduce the leading edge vortex LEV in strength and delay the LEV shedding. The increase of the maximum thickness can enhance the time-averaged thrust coefficient and the propulsive efficiency without lift reduction. As the maximum thickness location moves towards the leading edge, the airfoil obtains a larger time-averaged thrust coefficient and a higher propulsive efficiency without changing the lift coefficient.


Introduction
Since the Micro Air Vehicle MAV was generally defined by the Defense Advanced Research Projects Agency DARPA in 1997 1 , the Flapping Wing MAV FWMAV has been receiving more and more attention from military and civilian application domains.There is therefore an increasing interest to understand the aerodynamics of the flapping wing by experimental and numerical methods 2 .The first researchers who observed the unsteady flow dynamic characteristics of a flapping wing are Knoller 3 and Betz 4 , and in the middle of 1930s, von Kármán and Burgers gave a theoretical explanation for the different patterns of a largescale drag-indicative wake and a thrust-indicative wake 5 .It reinterested fluid scientists and biologists about two decades ago and now is a very active research area.Ellington gave a very comprehensive description of the insect hovering aerodynamics and unsteady aerodynamic effects were highlighted in these good series papers in 1984 6 .Anderson et al. showed that oscillating foil could have a very high propulsive efficiency, as high as 87%, under specific conditions by water tunnel experiments 7 .Dickinson et al. demonstrated three distinct mechanisms, delayed stall, rotational circulation, and wake capture in enhanced aerodynamic performance of insects using a robotic fly apparatus 8 .To investigate the flow field and effects of flapping parameters on the thrust generation and the propulsive efficiency numerically, an unsteady panel method 9 , and Navier-Stokes equations computations 10-15 have been employed during past decade, especially the latter are used more and more widely with benefits of continuous improvements of computers and Computational Fluid Dynamics CFD techniques.In addition, Shyy has a very good review of geometric similarity, scaling laws of birds, bats, insects, and artificial flying vehicles in 1999 16  Nevertheless, the influence of airfoil geometry on its flapping performance has not been explored enough.Most of researchers employed a flat plate 19 , an ellipse 20-22 , or specified NACA series airfoils as their study objects 7, 11, 14, 23 .However, Jones and Platzer showed that the airfoil thickness has very little effect on the propulsive efficiency for pure plunging motion with infinitesimal amplitude 23 .Lentink and Gerritsma concluded that a thin cambered airfoil outperformed a thick airfoil with respect to thrust coefficient and propulsive efficiency at very low Reynolds number Re O 150 24 .Ansari et al. studied a number of synthetic planform shapes in hovering flight and concluded that increasing aspect ratio, wing length, and wing area all enhance lift, albeit at different rates 25 .
The shape of wing section and other geometry parameters of birds, bats and insects are decided by the evolution of nature, but the parameters of wing sections airfoils of manmade FWMAVs could be selected by the designers.It is therefore useful to explore the effects of wing shape on its flapping performance.The present study mainly focuses on the influence of thickness of 2D NACA series symmetric airfoils with pure plunging motion which is the most popular style of current artificial flapping wing flying vehicles.Not only the influence of the maximum thickness magnitude but also the influence of the maximum thickness location on the flapping performance of the airfoil was discussed.
This paper is organized mainly into five parts.Following this introduction is the definition of the flapping model, several important parameters involved in the flapping model and how to calculate the time-averaged thrust coefficient and the propulsive efficiency are described.In Section 3, the governing equations for the 2-dimensional 2D unsteady flow field and the discritization of the computational domain are described.The computational method is validated by the published data for a rigid NACA0014 airfoil with plunging motion.After that, the influence of the magnitude of the maximum thickness on the timeaveraged thrust coefficient, the propulsive efficiency, and the lift coefficient time history are shown and analyzed in Section 4. Furthermore, the effects of locations of the maximum thickness based on a family of altered NACA0012 airfoils are described in Section 5, the "Class function/Shape function transformation" parametric method was employed to define the coordinates of these altered airfoils.Finally, a brief summary is concluded at the end of the paper.

Flapping Model
As the evolution of academic research on flapping wings 2 , some practical flapping wing micro air vehicles have been fabricated during past years.An unconventional example, shown in Figure 1, was proposed by Jones and Platzer in 2006 26 .In this design, the backward flapping biwing with a 25 cm span generates the thrust and the forward stationary wing provides lift.The present work is closely related to this model.
The real flapping model shown in Figure 1 is a 3-dimensional problem.However, it can be simplified to be a 2-dimensional flapping airfoil model with pure plunging motion.Under this simplification, the small-scale translation in flight direction and the trivial pitching motion along the leading edge are both ignored.Correspondingly, the wing tip effects and the unsteady flow along the wing span during flapping are also neglected.The flapping motion, that is, the pure plunging motion, is described using a harmonic function as follows: where y t stands for the plunging motion, c is the chord length, h is the dimensionless plunging amplitude with respect to the chord length, ω is the angular frequency in rad/s.One dimensionless parameter, called reduced frequency, is often used to describe the flapping frequency.The reduced frequency, k, is defined by where f is the flapping frequency in Hertz, u is the flow velocity of far field.This flapping model is illustrated in Figure 2.
In addition, there are two other important similarity numbers, Reynolds number Re and Strouhal number St .The former is employed to evaluate the flow viscosity, and the latter is used to make a comparison between the flapping speed and the far field flow speed.where ν is the kinematic viscosity of concerned fluid, A is the wake width and can be estimated using the peak-to-peak excursion of the trailing edge, or more simply by twice the plunging amplitude.
In classical aerodynamics, the lift coefficient C L , the drag coefficient C D , and the pitch moment coefficient C M can be defined as follows: where F x and F y are the components of resulting aerodynamics force along horizontal parallel with u direction and vertical normal to u direction directions, ρ is the flow density, S is the reference area and equals to c in value for 2D problem, L is the reference length and equals to c here.Then, the time-averaged thrust coefficient C T and the power input coefficient C P in one flapping cycle is calculated by 2.5 and 2.6 .Correspondingly, the propulsive efficiency is defined using 2.7 where T are the period in seconds and T 1/f, ẏ t is the first-order time derivations of y t .
For simplicity, we also use C T to stand for C T in the next sections.

Governing Equations
The commercial CFD solver FLUENT v6.3.26 was employed to simulate the unsteady flow fields around moving airfoils with predefined motions.The two-dimensional timedependent Navier-Stokes equations were solved using the finite volume method, assuming incompressible laminar flow.The mass and momentum equations were solved in a fixed inertial reference frame incorporating a dynamic mesh.The dimensionless mass and momentum conservation equations are given by 22 where u and p are the dimensionless velocity and dimensionless pressure, respectively.

Grid Scheme
The hybrid mesh which is shown schematically in Figure 3 was employed to analyze the flow field.The computational domain is divided into two distinct zones: moving zone and remeshing zone.The moving zone consists of C-type structured quadrilateral mesh, and the remeshing zone unstructured triangular mesh.The airfoil is located in the center of the computational domain, and has a no-slip wall boundary condition applied.The spacial scale of each zone and corresponding boundary condition are also shown in Figure 3.The whole moving zone mesh, including the interfaces between these two zones, moves with the airfoil together according to the predefined airfoil motion.This means remeshing only occurs at a distance of 20 to 45 reference length away from the airfoil body, which ensures that the flow simulation around the airfoil is a little affected by the moving mesh.The C-type mesh in the very near neighborhood of the airfoil is shown in Figure 4, with 201 nodes on each single airfoil surface.The hybrid mesh was generated in GAMBIT v2.3.16.

Validations
A grid sensitivity study was carried out first to evaluate the independence of the numerical solution on the mesh size.Some specified unsteady flow fields around a rigid NACA0014 airfoil with pure plunging motion were computed under conditions of k 2.0, h 0.4, u 34.7 m/s Ma 0.1 , c 0.064 m, and Re 1.0 × 10 4 .The sizes of involved grids were 201 × 201 nodes 201 along every single airfoil surface, 201 in vertical direction with the thickness of the first layer grid around the airfoil of 0.0002c, 401 × 201 nodes with the first layer thickness equals to 0.0002c, 201 × 101 nodes with the first layer thickness of 0.0005c, 201 × 101 nodes with the first layer thickness of 0.0001c, 201 × 101 nodes with a first layer thickness of 0.00005c, and 201 × 101 nodes with the first layer thickness of 0.0002c.Each grid scheme was simulated in 7 periods, and both of the lift time history and drag time history showed smooth periodic properties after 3 or 4 cycles.There were 500 time steps in every plunging cycle.Time histories of drag coefficients and lift coefficients based on these six grid   schemes match very well, as shown in Figures 5 and 6.However, the 201 × 101 grid with the fist layer thickness of 0.0002c was employed for the next simulations.
Furthermore, to validate the accuracy of the present hybrid mesh, simulations based on the 201 × 101 size grid were performed with 1000 and 500 time steps in one plunging period under conditions of k 2.0, h 0.4, u 34.7 m/s Ma 0.1 , c 0.064 m, and Re 1.0 × 10 4 , respectively.The coupling between the pressure and the velocity was achieved by means of a SIMPLE algorithm.Meanwhile, the discretizations of pressure and momentum terms were the second-order scheme and the second-order Upwind scheme.The accuracy was set to be double-precision and the time discretization was the first-order implicit, which is one of the most straightforward methods in FLUENT for dynamic mesh module 27 .The time variation of the plunging position and the time histories of C D in one period are shown in Figure 7 with the comparisons with results obtained by Tuncer and Kaya 11 and Miao and Ho 14 .These four results are in good agreement, though different mesh schemes were employed in these studies.Furthermore, the time history of C L with 10 times plunging position is shown in Figure 8.  size grid with first layer thickness of 0.0002c, and 500 time steps were employed for all of the next simulations.
The three wiggles of the lift coefficient near the peak value, marked as 1 , 2 , and 3 in Figure 8, are caused by the coupling of the formation of the leading edge vortex LEV and the interaction of the airfoil with the LEV shed during the previous half period PLEV .The pressure contour of the airfoil neighborhood at wiggle 1 is shown in the left of Figure 9.At this time instant, the airfoil is meeting the PLEV while the LEV is just forming.The low pressure region under the airfoil reduces the lift coefficient, and the low pressure region above the airfoil enhances the lift coefficient.The former overwhelmingly dominates in the coupling effects, which results in the decreasing of the lift coefficient.From wiggle 1 to the wiggle 2 , the low pressure region caused by the LEV is increasing with the PLEV moving to the wake, and the LEV can balance the PLEV at wiggle 2 at which pressure contour is shown in the middle of Figure 9. From the wiggle 2 to wiggle 3 , the low pressure region caused by the LEV continues to increase and results in increasing of the lift coefficient.Finally, the LEV is shedding and the lift coefficient achieves its maximum value at wiggle 3 at which pressure contour is shown in the right of Figure 9.It should be mentioned that these phenomena depend on the kinematics of the airfoil.If one of the kinematics, for example, the plunging amplitude, is changed, the phenomena may appear in another way.

Influence of Maximum Thickness
This section mainly focuses on the effects of the magnitudes of the maximum thickness.
The popular NACA 4-digit series airfoils are employed, and the conventional definition of symmetrical 4-digit NACA airfoils can be described by the following 28 where x is the position along the chord from 0.0 to c, y is the half thickness at a given value of x, t is the maximum thickness as a fraction of the chord.This definition results in the location of the maximum thick of at 30% 0.3c as illustrated in Figure 10.
Computations for pure plunging motion of a family of airfoils from a NACA0002 airfoil to a NACA0020 airfoil were performed under the conditions of k 2.0, h 0.4, u 10.0 m/s, c 0.1 m, and Re 1.0 × 10 4 , where these flapping kinematics were modified from the Tuncer and Kaya's paper 11 and Miao and Ho's paper 14 with the same reduced frequency, the same Reynolds number and the same Strouhal number.The relationships of the time-averaged thrust coefficient and the propulsive efficiency with the ratio of the maximum thickness to the chord are shown in Figure 11.It is observed that the time-averaged thrust coefficient and the propulsive efficiency increase almost linearly with the maximum thickness ratio up to a value of about 0.12, though the rates of increasing speeds are a little different.After that, the effect of the maximum thickness becomes less intensive, but is still not negligible.The interesting thing is that the lift histories are very similar for these airfoils as shown in Figure 12.Pressure contours of the NACA0006 airfoil and the NACA0020 airfoil at n 0 T time instant, n 1/4 T time instant, n 1/2 T time instant, and n 3/4 T time instant are shown in Figures 13 and 14, respectively.It is observed that the increase of an airfoil thickness can reduce the LEV strength and delay the LEV shedding.However, the increase of the thickness also induces the change of the airfoil geometry.This change results in the change of the direction of the pressure exerted on the airfoil, as the pressure must be perpendicular to the surfaces of the airfoil.The increase of the time-averaged thrust coefficient and the propulsive efficiency without lift reduction are induced by the coupling effects of the interaction between the airfoil and the LEV and the change of the airfoil geometry.

Influence of Maximum Thickness Location
In order to check the effect of the maximum thickness location Loc on the flapping performance, the "Class function/Shape function Transformation CST " parametric geometry representation method is employed to parameterize the 2D airfoil.The universal CST method is defined as follows 29 where C x/c is the class function, S x/c is the shape function, and Δz TE is the thickness of trailing edge.S x/c can be described using Bernstein polynomials of any specified order n, and in most of CFD simulations, Δz TE is set to be 0.0.When the order of Bernstein polynomials BPO is 2, 5.1 can be written as 5.2 for round nose airfoil, where A i i 0, 1, 2 are the parameters to be determined, and A 0 2R LE /c, where R LE is the radius of leading edge.
A i can be determined by the least square fitting method based on the original NACA0012 coordinates obtained by 4.1 .As shown in Figure 15, the coordinates defined by CST method match very well with the original coordinates, corresponding A i are shown in Table 1.The leading edge radius and the value of the maximum thickness were kept constant in the next series of simulations.Keeping the value of the first coefficient, A 0 , constant keeps the leading edge radius constant.The analytic equation to specify the location of maximum thickness for a constant maximum thickness can be developed as follows.Denote y/c and x/c as y and x for simplicity, then the thickness of the airfoil, Th x , is described as follows: Suppose that Th x gets its maximum value at x 0 , then dTh x 0 dx 0 0 5.4 must be satisfied.Substitute 5.3 to 5.4 , is obtained.Meanwhile, should be satisfied.After x 0 is specified, A 1 and A 2 can be obtained numerically using 5.5 and 5.6 .The values of A i for different x 0 /c are shown in Table 1, and the geometries of these five airfoils are shown in Figure 16.
Computations for the unsteady flow fields with pure plunging motion of a series of altered NACA0012 airfoils with Loc 30.0% to Loc 50.0% were performed under the conditions of k 2.0, h 0.

Conclusions
The unsteady flow fields for a pure plunging motion of a family of airfoils from an NACA0002 airfoil to a NACA0020 airfoil, and a series of altered NACA0012 airfoils were analyzed using CFD techniques.Some interesting results can be concluded here.Under specified conditions, k 2.0, h 0. airfoil with Loc 50% decrease by 17.8% and 11.9%, respectively, compared to ones of the original NACA0012 airfoil with Loc 30%.Though the numerical simulations were just performed under specified conditions and a single reduced frequency, we hope these conclusions shed a little light to the explanation of the difference between the speed of flat-plat-like wing flying animals, insects and thick-airfoil-like wing flying animals, birds.It is not sure if a similar phenomenon appears in pure pitching motion, combination of plunging motion and pitching motion, and a complete 3D wing.These are recommendations for future research.
, and Sane gave a very detailed review of aerodynamics characteristics of insect flight in 2003 17 .Very recently, Shyy et al. reviewed the recent progress in flapping wing aeroelasticity, both of the spanwise and chordwise flexibility in flapping wings were summarized 18 .

Figure 2 :
Figure 2: Illustration of the flapping model in pure plunging motion.

Figure 3 :
Figure 3: Hybrid mesh topology with boundary conditions.

Figure 4 :
Figure 4: C-type grid very near the airfoil.
schemes match very well, as shown in Figures5 and 6.However, the 201 × 101 grid with the fist layer thickness of 0.0002c was employed for the next simulations.Furthermore, to validate the accuracy of the present hybrid mesh, simulations based on the 201 × 101 size grid were performed with 1000 and 500 time steps in one plunging period under conditions of k 2.0, h 0.4, u 34.7 m/s Ma 0.1 , c 0.064 m, and Re 1.0 × 10 4 , respectively.The coupling between the pressure and the velocity was achieved by means of a SIMPLE algorithm.Meanwhile, the discretizations of pressure and momentum terms were the second-order scheme and the second-order Upwind scheme.The accuracy was set to be double-precision and the time discretization was the first-order implicit, which is one of the most straightforward methods in FLUENT for dynamic mesh module 27 .The time variation of the plunging position and the time histories of C D in one period are shown in Figure7with the comparisons with results obtained by Tuncer and Kaya 11 and Miao and Ho 14 .These four results are in good agreement, though different mesh schemes were employed in these studies.Furthermore, the time history of C L with 10 times plunging position is shown in Figure8.Figures7 and 8 also show that 500 time steps in one cycle are good enough to get the complete details.Based on above-mentioned validations, 201 × 101

−Figure 5 :Figure 6 :
schemes match very well, as shown in Figures5 and 6.However, the 201 × 101 grid with the fist layer thickness of 0.0002c was employed for the next simulations.Furthermore, to validate the accuracy of the present hybrid mesh, simulations based on the 201 × 101 size grid were performed with 1000 and 500 time steps in one plunging period under conditions of k 2.0, h 0.4, u 34.7 m/s Ma 0.1 , c 0.064 m, and Re 1.0 × 10 4 , respectively.The coupling between the pressure and the velocity was achieved by means of a SIMPLE algorithm.Meanwhile, the discretizations of pressure and momentum terms were the second-order scheme and the second-order Upwind scheme.The accuracy was set to be double-precision and the time discretization was the first-order implicit, which is one of the most straightforward methods in FLUENT for dynamic mesh module 27 .The time variation of the plunging position and the time histories of C D in one period are shown in Figure7with the comparisons with results obtained by Tuncer and Kaya 11 and Miao and Ho 14 .These four results are in good agreement, though different mesh schemes were employed in these studies.Furthermore, the time history of C L with 10 times plunging position is shown in Figure8.Figures7 and 8 also show that 500 time steps in one cycle are good enough to get the complete details.Based on above-mentioned validations, 201 × 101

Figure 12 :
Figure 12: Time histories of C L .
4, u 10.0 m/s, c 0.1 m, and Re 1.0 × 10 4 .The relationships of the time-averaged thrust coefficient and the propulsive efficiency with the ratio of the maximum thickness location to chord are shown in Figure 17.It is observed that the timeaveraged thrust coefficient and the propulsive efficiency almost decrease linearly with the ratio of location of the maximum thickness.Similar to the former results, the lift coefficient

Figure 17 :Figure 18 :
Figure 17: Trends of C T and η versus Loc.
history is very similar to each other as illustration in Figure18.Pressure contours of the original NACA0012 airfoil and the altered NACA0020 airfoil with Loc 50.0% at n 0 T time instant, n 1/4 T time instant, n 1/2 T time instant, and n 3/4 T time instant are shown in Figures19 and 20, respectively.It is observed that the pressure contours near these two airfoils are very similar to each other.The difference of the flapping performance for these two airfoils is caused mainly by the geometry, because the pressure exerted on an airfoil is perpendicular to the surfaces of the airfoil.

Table 1 :
Values of A i for different x 0 /c.