Abstract

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.

1. 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 large-scale 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 [1015] 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], 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].

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 [2022], 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 man-made 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 time-averaged 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.

2. 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 stands for the plunging motion, is the chord length, 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, , is defined by where is the flapping frequency in Hertz, 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 () and Strouhal number (). 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. and are described by where is the kinematic viscosity of concerned fluid, 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 , the drag coefficient , and the pitch moment coefficient can be defined as follows: where and are the components of resulting aerodynamics force along horizontal (parallel with direction) and vertical (normal to direction) directions, is the flow density, is the reference area and equals to in value for 2D problem, is the reference length and equals to here. Then, the time-averaged thrust coefficient and the power input coefficient in one flapping cycle is calculated by (2.5) and (2.6). Correspondingly, the propulsive efficiency is defined using  (2.7) where are the period in seconds and , is the first-order time derivations of . For simplicity, we also use to stand for in the next sections.

3. Approach for Computations

3.1. 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 time-dependent 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 and are the dimensionless velocity and dimensionless pressure, respectively.

3.2. 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.

3.3. 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 , ,  m/s (Ma = 0.1),  m, and . The sizes of involved grids were nodes (201 along every single airfoil surface, 201 in vertical direction) with the thickness of the first layer grid around the airfoil of 0.0002, nodes with the first layer thickness equals to 0.0002, nodes with the first layer thickness of 0.0005, nodes with the first layer thickness of 0.0001c, nodes with a first layer thickness of 0.00005, and nodes with the first layer thickness of 0.0002. 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 grid with the fist layer thickness of 0.0002 was employed for the next simulations.

Furthermore, to validate the accuracy of the present hybrid mesh, simulations based on the size grid were performed with 1000 and 500 time steps in one plunging period under conditions of , ,  m/s (Ma = 0.1),  m, and Re = , 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 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 with 10 times plunging position is shown in Figure 8. Figures 7 and 8 also show that 500 time steps in one cycle are good enough to get the complete details. Based on above-mentioned validations, size grid with first layer thickness of 0.0002, 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.

4. 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 is the position along the chord from 0.0 to , is the half thickness at a given value of , is the maximum thickness as a fraction of the chord. This definition results in the location of the maximum thick of at 30% (0.3) 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 , ,  m/s,  m, and Re = , 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 time instant, time instant, time instant, and 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.

5. 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 is the class function, is the shape function, and is the thickness of trailing edge. can be described using Bernstein polynomials of any specified order , and in most of CFD simulations, 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 are the parameters to be determined, and , where is the radius of leading edge.

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 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, , 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 and as and for simplicity, then the thickness of the airfoil, , is described as follows: Suppose that gets its maximum value at , then must be satisfied. Substitute (5.3) to (5.4), is obtained. Meanwhile, should be satisfied. After is specified, and can be obtained numerically using (5.5) and (5.6). The values of for different 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 , ,  m/s,  m, and Re = . 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 time-averaged 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 history is very similar to each other as illustration in Figure 18. Pressure contours of the original NACA0012 airfoil and the altered NACA0020 airfoil with Loc = 50.0% at time instant, time instant, time instant, and time instant are shown in Figures 19 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.

6. 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, , ,  m/s,  m, and Re = , corresponding to St = 0.2546, the increase of an airfoil thickness can reduce the leading edge vortex in strength and delay the leading edge vortex shedding. The coupling effects of the LEV and the airfoil geometry induce the increase of the time-averaged thrust coefficient and the propulsive efficiency without lift reduction. Particularly, the time-averaged thrust coefficient of a NACA0020 airfoil is about 12.5 times of the time-averaged thrust coefficient of a NACA0002 airfoil, and similar feature happens to the propulsive efficiency. Furthermore, the closer that the location of the maximum thickness is to the leading edge, the larger the time-averaged thrust coefficient and the higher the propulsive efficiency with very similar lift coefficient time history. Particularly, the time-averaged thrust coefficient and the propulsive efficiency of an altered NACA0012 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.

Acknowledgments

The Grant supported from National Science Foundation of China (no. 10972034) and National Science Foundation for Post-doctoral Scientists of China (20090460216) are greatly acknowledged. The first author would like to express his great appreciation to Dr. Brenda M. Kulfan, from the Boeing Company, for her invaluable suggestions on the CST method.