AP04_1.vp 1 Introduction Due to agile flight of birds and insects, flapping wing pro- pulsion has already been recognized to be more efficient than conventional propellers for very small scale vehicles with wing spans of 15 cm or less, so-called micro-air vehicles (MAVs). Since the primary mission for MAVs is surveillance, they are desired to have good maneuverability and sustained flights with flight speeds of 30 to 60 kph. A current interest in the research and development community is to find the most en- ergy efficient airfoil adaptation and wing motion technologies capable of providing the required aerodynamic performance for MAV flight. Recent experimental and computational studies investi- gated the propulsive characteristics flapping airfoils, and shed some light on the relationship among the produced thrust, the amplitude and frequency of the flapping oscilla- tions, and the flow speed. Water tunnel flow visualization experiments on flapping airfoils conducted by Lai and Platzer [1] and Jones et al [2] provide a considerable amount of infor- mation on the wake characteristics of thrust producing flap- ping airfoils. In their experiments, Anderson et al [3] observed that the phase angle between pitch and plunge os- cillations plays a significant role in maximizing the propulsive efficiency. Navier-Stokes computations have been performed by Tuncer et al [4, 5, 6] and by Isogai et al [7, 8] to explore the effect of flow separation on thrust and propulsive efficiency of a single flapping airfoil in combined pitch and plunge oscil- lations. The experimental and numerical studies by Jones et al [9, 10, 11] and Platzer and Jones [12] on flapping- -wing propellers points at the gap between numerical results and the actual flight conditions in high frequency motions, and the limitation or enhancement of the performance of flapping airfoils by the onset of dynamic stall. Jones and Platzer [11] recently demonstrated a radio-controlled micro air vehicle propelled by flapping wings in biplane configura- tion (Fig. 1). The computational and experimental findings show that thrust generation and propulsive efficiency of flapping air- foils are closely connected to the flapping motion and flow parameters, such as the unsteady flapping velocity, frequency and amplitude of the pitch and plunge motions, the phase shift between them, and the air speed. It is apparent that to maximize the thrust and/or propulsive efficiency of a flapping airfoil, an optimization of all the above parameters is needed. In the earlier study [13], a gradient based numerical optimi- zation method has been applied to a flapping airfoil to maxi- mize its thrust. The preliminary results with a limited number of optimization variables compared well with the parametric studies performed earlier. In this study the optimization method developed earlier is extended to accommodate an objective function, which is a linear combination of the thrust and the propulsive efficiency of a flapping airfoil (Fig. 2). The optimization variables are taken to be the pitch and plunge amplitudes, h0 and �0, and the phase shift between the pitch and plunge motions, �. The flowfield around flapping airfoils are discretized using over- set grids. Unsteady flow solutions required to evaluate the gradients of the objective function by perturbation of the © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 31 Czech Technical University in Prague Acta Polytechnica Vol. 44 No. 1/2004 Optimization of Flapping Airfoils for Maximum Thrust and Propulsive Efficiency I. H. Tuncer, M. Kay A numerical optimization algorithm based on the steepest decent along the variation of the optimization function is implemented for maximizing the thrust and/or propulsive efficiency of a single flapping airfoil. Unsteady, low speed laminar flows are computed using a Navier-Stokes solver on moving overset grids. The flapping motion of the airfoil is described by a combined sinusoidal plunge and pitching motion. Optimization parameters are taken to be the amplitudes of the plunge and pitching motions, and the phase shift between them. Computations are performed in parallel in a work station cluster. The numerical simulations show that high thrust values may be obtained at the expense of reduced efficiency. For high efficiency in thrust generation, the induced angle of attack of the airfoil is reduced and large scale vortex formations at the leading edge are prevented. Keywords: optimization, flapping airfoils, unsteady aerodynamics, moving overset grids. Fig. 1: MAV with flapping wings [11] Fig. 2: Flapping motion of an airfoil in combined plunge and pitch optimization variables are computed in parallel in a com- puter cluster. 2 Numerical method The unsteady viscous flowfield around a flapping airfoil is computed by solving the Navier-Stokes equations on moving overset grids. The flow variables at the intergrid boundaries are interpolated from the donor grid. Computations on each subgrid are performed in parallel. PVM message passing library routines are used in the parallel solution algorithm [14]. The computed results are analyzed in terms of average thrust coefficient and propulsive efficiency values, and instan- taneous particle traces. The computational domain is discretized with overset grids. C-type grid around the airfoil is overset onto a Carte- sian background grid (Fig. 3). The flapping motion of the air- foil is imposed by moving the airfoil and the grid around it on the background grid. The flapping motion of the airfoil in combined plunge, h, and pitch, �, is specified by � �h h t� � 0 cos � , � �� �� � �0 cos � �t , where the angular frequency � is given in terms of the reduced frequency, k c U� �� . The pitching motion is about the mid-chord location. 2.1 Intergrid boundary conditions At the intergrid boundaries formed by the overset grids, the conservative flow variables are interpolated in each time- step of the unsteady solution. Intergrid boundary points are first localized in a triangular stencil in the donor grid by a di- rectional search algorithm. The localization process provides the interpolation weights to interpolate the flow variables within the triangular stencil [14]. 2.2 Optimization Optimization process is based on following the direction of the steepest ascent of an objective function, O. The direc- tion of the steepest ascent is given by the gradient vector of the objective function: � ��O O V O V nv v v� � � � � � �1 1 2 2 �, where Vn’s are the optimization variables of the objective function. The objective function is taken as a linear combination of the average thrust, Ct, and the propulsive efficiency, �, over a flap- ping period. � � 0 sets the objective function to a normalized thrust coefficient: O C C C C t t t t ( , ) ( )� � � � � � � � � � � � � 1 � �D D C T C tt d t t T � � � 1 d � � � � � � � C U T p t t S t t T 1 V Ad d where the denominator in the efficiency expression accounts for the average work required to maintain the flapping mo- tion. � denotes the optimization stepsize. The components of the gradient vector is then evaluated numerically by computing the objective function for a small perturbation of the optimization variables one at a time. It should be noted that the evaluation of these vector compo- nents requires an unsteady flow computation over a few periods of the flapping motion until a periodic behavior is reached. Once the unit vector in the ascent direction is evalu- ated by D � � � O O , the step size �V D� is to be determined. Reference [15] suggests that the value of � should be based on the Hessian of the objective function, which involves the sec- ond derivatives of the objective function with respect to opti- mization variables. The exact computation of the Hessian is expensive and the cost is proportional to the number of optimization parameters. In this work an approximation is made [13], and the step size is evaluated as follows: � � � � � �� O O D where � � � � � � �� � � O O O e n n n � � � � � D v v 1 1 . � ���O � D is then evaluated numerically. 2.3 Parallel Computation A coarse parallel algorithm based on domain decomposi- tion is implemented in a master-worker paradigm [16]. The overset grid system is decomposed into its subgrids first, and the solution on each subgrid is assigned to a separate proces- sor in the computer cluster. In addition, the background grid may also be partitioned to improve the static load balanc- ing. Intergrid boundary conditions are exchanged among subgrid processes. PVM (version 3.4.4) library routines are used for inter-process communication. In the optimization process, unsteady flow computations with perturbed optimi- zation variables, which are required to determine the gradient vector, are also evaluated in parallel. Computations are performed in a cluster of computers with dual Pentium-III processors and Linux operating system. 32 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No. 1/2004 Czech Technical University in Prague Fig. 3: Overset grid system 3 Results and discussion In this study the flapping motion of a NACA00l2 airfoil in a combined plunge and pitch is investigated. The reduced frequency of the oscillatory motion is fixed at k � 1. The optimization variables are the plunge and pitch amplitudes, (h0, �0) and the phase shift between plunge and pitch motions (�). All the flows are assumed to be laminar, and computed at Re � 1�104 and M � 0.1. Table 1 summarizes the optimization cases studied, and the initial values of the optimization vari- ables. The parallel computations with 8 processors take about 20–30 hours of wall clock time for a typical optimization case. In Case 1, where � � 0, the average thrust coefficient is maximized. The instantaneous variation of the unsteady drag (negative thrust) coefficient along a few optimization steps is shown in Fig. 4. As the optimization variables are incre- mented along the optimization steps, unsteady computations are carried out for a few periods of the flapping motion until a periodic behavior is obtained. The variation of the average thrust coefficient and the propulsive efficiency with respect to the optimization variables are given in Fig. 5. As shown, as the optimization variables are incremented along the gradi- ent of the objective function, the average thrust coefficient increases gradually, and a maximum value of 1.41 is reached © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 33 Czech Technical University in Prague Acta Polytechnica Vol. 44 No. 1/2004 Case � h0 �� � 1 0.0 0.5 5 30 2 0.5 0.5 5 30 3 1.0 0.5 5 30 4 0.0 0.5 25 60 5 0.0 1.0 5 60 6 0.0 1.0 25 90 7 1.0 0.5 25 60 8 1.0 1.0 5 60 9 1.0 1.0 25 90 Table 1: Optimization cases and starting conditions Fig. 4: Cd history along optimization steps, Case 1 Fig. 5: Maximization of thrust coefficient (� � 0), Case 1 Fig. 6: Maximization of propulsive efficiency and thrust coefficient (� � 0.5), Case 2 at h0 � 1.60, �0 � 23.5 and � � 103.4 deg. The corresponding propulsive efficiency is 28.3 %. Optimization steps for Cases 2 and 3 are shown in Figs. 6 and 7, respectively. In case 2, where � � 0.5, the average thrust 34 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No. 1/2004 Czech Technical University in Prague Fig. 7: Maximization of propulsive efficiency (� � 1.0), Case 3 Fig. 8: Instantaneous particle traces at the instant of maximum thrust along the optimization steps for Cases 1 and 2 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 35 Czech Technical University in Prague Acta Polytechnica Vol. 44 No. 1/2004 Fig. 9: Instantaneous particle traces along a period of the optimized flapping motion for Case 1 Fig. 10: Maximization of thrust coefficient (� � 0) with various starting conditions, Cases 1, 4–6 Fig. 11: Maximization of propulsive efficiency (� � 1) with various starting conditions, Cases 3, 7–9 and the propulsive efficiency have equal weights in the objec- tive function. As a results, the efficiency is improved at the expense of average thrust. It is observed that the higher effi- ciency is achieved at a lower plunge amplitude and a higher pitch amplitude. The phase shift slightly drops to 97.8 deg. In Case 3 the propulsive efficiency is maximized at the low pitch and plunge amplitudes with the corresponding very low thrust coefficient. It is apparent that the propulsive efficiency and the thrust production of flapping airfoils are inversely proportional. The unsteady flow fields along the optimization steps are investigated with particle traces. The particles are emitted along a straight line in the vicinity of the leading edge of the airfoil, and convected in the flow field with the local velocity. The line from which the particles are emitted follows the lead- ing edge of the airfoil to capture the leading edge vortex formations in more detail. In Fig. 8, the instantaneous parti- cle traces at the instant of maximum thrust (minimum drag) in a flapping period are given along the optimization steps of Cases 1 and 3. It is observed that in Case 1, the leading edge vortex formation is promoted along the optimization steps. The maximum instantaneous thrust occurs at about the mean amplitude location as the leading edge vortex forms, just before the suction field at the leading edge collapses as the leading edge vortex develops stronger. Whereas, in Case 3, the leading edge vortex formation is prevented along the optimization steps, which incidently maximizes the propul- sive efficiency. The unsteady flow becomes more streamlined with the motion of the airfoil. Fig. 8 shows the optimized flow field for maximum thrust in Case 1. The flowfield is observed to be highly vortical with strong leading edge vortices during the upstroke and the downstroke. The flow field is periodic, and antisymmetric along the upstroke and the downstroke. Next, the optimiza- tion space is searched for other possible local maximums of the objective function for Cases 1 and 3. It is implemented by initiating the optimization process from various initial con- ditions as given in Table 1. All the results of the optimization cases are also given in Table 2. The initial conditions and the optimized states at the end of the optimization processes are shown in Figs. 10 and 11 for � � 0 and � � 1, respectively. Fig. 11 reveals that all the optimization cases for � � 0 con- verge about the same value of the objective function, which is the thrust coefficient, and of the optimization variables. It suggests that the global maximum of the objective function may have been found. On the other hand, the optimization processes for � � 1 provides different optimum states for h0 36 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No. 1/2004 Czech Technical University in Prague Fig. 12: Optimized flapping motions Case h0 �0 � Ct � [%] 1 1.60 23.5 103.4 1.41 28.3 2 1.36 29.6 97.8 1.08 44.1 3 0.45 15.4 82.4 0.08 58.5 4 1.73 23.8 100.7 1.44 25.4 5 1.52 26.9 87.2 1.27 33.4 6 1.55 28.6 94.9 1.45 35.9 7 0.57 21.0 86.7 0.13 63.8 8 0.60 22.8 86.1 0.13 64.8 9 0.83 35.6 86.5 0.18 67.5 Table 2: Optimization results and �0, and about the same � values. It appears that a high flapping efficiency may be achieved for a range of h0 and �0 values, such that ao increases as h0 does. The optimum flapping motions for Cases 1–3 and 9 are shown in Fig. 12. It is clearly observed that the plunge ampli- tude plays a significant role in thrust generation. It also appears that in order to improve the efficiency, the plunge amplitude is to be reduced, and the pitch amplitude is to be increased. In addition, the phase shift between the plunge and pitch motions, which is about 90 deg for all the cases, reduces the effective angle of attack at the mid plunge loca- tion, where the plunge velocity is maximum. The variations of the effective angle of attack the airfoil sees along a flapping period are given in Fig. 13 for Cases 1–3. The period starts at 0 deg, which corresponds to the h � �h0 position of the airfoil. In agreement with the previous observations, for higher thrust production, as in Cases 1 and 2, a flapping airfoil stays at large effective angles of attack for a large fraction of the flapping period. For an efficient flapping as in Case 3, the effective angle of attack at the mid-plunge location (� � 90, 270 deg, h � 0) is set about 0 deg. Whereas, in Cases 1 and 3 the maximum effective angle of attack occurs around mid-plunge locations. 4 Concluding remarks A gradient based numerical optimization is successful- ly applied to the thrust generation and propulsive efficiency of an airfoil flapping in a combined plunge and pitch. The optimization of thrust generation and propulsive effi- ciency together is achieved with a weighted and normalized objective function. The parallel implementation of the optimization algorithm is shown to be quite robust. Thrust generation of a flapping airfoil is maximized at large plunge amplitudes as large leading edge vortices form and shed into the wake. The airfoil stays at a large effective angle of attack during the most of the flapping period. On the other hand, the propulsive efficiency of flapping airfoils may be in- creased by reducing the plunge amplitude and the effective angle of attack, and consequently by preventing the forma- tion of leading edge vortices. Further research is in progress to implement the present optimization method to the thrust generation of flapping airfoils in a biplane configuration. Acknowledgment The partial support provided by NATO Research and Technology Organization under the project TX01–02 is acknowledged. References [1] Lai J. C. S., Platzer M. F.: The Jet Characteristics of a Plunging Airfoil. 36th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, January 1998. [2] Jones K. D., Dohring C. M., Platzer, M. F.: “An Ex- perimental and Computational Investigation of the Knoller-Betz Effect.” AIAA Journal,Vol. 36 (1998), No. 7, p. 1240–1246. [3] Anderson J. M., Streitlien K., Barrett D. S., Triantafyllou M. S.: “Oscillating Foils of High Propulsive Efficiency.” Journal of Fluid Mechanics, Vol. 360 (1998), p. 41–72. [4] Tuncer I. H, Platzer M. F.: “Thrust Generation due to Airfoil Flapping.” AIAA Journal, Vol. 34 (1995), No. 2, p. 324–331. [5] Tuncer I. H., Lai J., Ortiz M. A., Platzer M. F.: Unsteady Aerodynamics of Stationary / Flapping Airfoil Combination in Tandem. AIAA Paper 97-0659, 1997. [6] Tuncer I. H., Platzer M. F.: “Computational Study of Flapping Airfoil Aerodynamics.” AIAA Journal of Aircraft, Vol. 35 (2000), No. 4, p. 554–560. [7] Isogai K., Shinmoto Y., Watanabe Y.: “Effects of Dynamic Stall on Propulsive Efficiency and Thrust of a Flap- ping Airfoil.” AIAA Journal, Vol. 37 (2000), No. 10, p. 1145–1151. [8] Isogai K., Shinmoto Y.: Study on Aerodynamic Mechanism of Hovering Insects. AIAA Paper No. 2001–2470, June 2001. [9] Jones K. D., Duggan S. J., Platzer M. F.: Flapping-Wing Propulsion for a Micro Air Vehicle. AIAA Paper No. 2001–0126, 39th AIAA Aerospace Sciences Meeting, Reno, Nevada, January 2001. [10] Jones K. D., Castro B. M., Mahmoud O., Pollard S. J., Platzer M. F., Neef M. F., Hummel D.: A Collobirative Nu- merical and Experimental Investigation of Flapping-Wing Propulsion. AIAA Paper No. 2002–0706, January 2002. [11] Jones K. D., Platzer M. F.: Experimental Investigation of the Aerodynamic Characteristics of Flapping-Wing Micro Air vehicles. AIAA Paper No. 2003–0418, January 2003. [12] Platzer M. F., Jones K. D.: The Unsteady Aerodynamics of Flapping-Foil Propellers. 9th International Symposium on Unsteady Aerodynamics, Aeroacoustics and Aero- elasticity of Turbomachines, cole Centrale de Lyon, Lyon (France), September 4–8, 2000. [13] Tuncer I. H., Kaya M.: Optimization of Flapping Airfoils for Maximum Thrust. AIAA Paper No. 2003–0420, January 2003. [14] Tuncer I. H: “A 2-D Unsteady Navier-Stokes Solution Method with Moving Overset Grids.” AIAA Journal, Vol. 35 (1997), No. 3, p. 471–476. [15] Kuruvila G., Ta’asan S., Salas M. D.: Airfoil Optimization by the One-Shot Method. AGARD Report No. 803, 1994. © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 37 Czech Technical University in Prague Acta Polytechnica Vol. 44 No. 1/2004 Fig. 13: Effective angle of attack [16] Tuncer I. H.: Parallel Computation of Multi-Passage Cascade Flows with Overset Grids. Selected Papers from Parallel Computational Fluid Dynamics Workshop, Gulcat U., Emerson D. R. (eds.), Istanbul Technical Uni- versity, Istanbul, Turkey, 1999, p. 81–89. Assoc. Prof. Ismail H.Tuncer e-mail: tuncer@ae.metu.edu.tr Mustafa Kaya, Graduate Research Assistant Department of Aerospace Engineering Middle East Technical University 06531 Ankara, Turkey 38 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No. 1/2004 Czech Technical University in Prague