Numerical optimization of flapping foil kinematics using high－order spectral difference method

YU Meilin，WANG Zhijian，HU Hui

（1.University of Maryland，Baltimore County，Baltimore，MD 21250，USA；2.TheUniversity of Kansas，Lawrence，KS 66045，USA；3.Iowa State University，Ames，IA 50011，USA）

Nomenclature

0 Introduction

Unsteady flapping-wing aerodynamics has attracted considerable research attention recently due to the increased interestin designing Micro Air Vehicles(MAVs)and NanoAir Vehicles(NAVs)for a variety of missions，including reconnaissance and surveillance，targeting and sensing.Comprehensive reviews［1-8］are available for the latest progress in this area，but there are still many open problems needed to be tackled before our current understanding on unsteady flapping flight can be systematically integrated into design of aerodynamically efficient MAVs and NAVs.Numerical optimization is such a powerful tool to address these issues.On one hand，numerical optimization can provide usefulgeometric and kinematic information for the design of MAVs and NAVs;on the other hand，the optimization outputs(e.g.，flow fields)can further enhance our understanding on unsteady flapping flight.One example to demonstrate the power of numerical optimization on both roles is the work by Berman and Wang［9］.In this work，they reported the energy-minimizing kinematics of several hovering insects.Interestingly，by examining the optimization results，they also explained why in the optimal kinematics the stroke maintains a constant leading edge throughout a flapping period even when a symmetry between the leading and trailing edges of the wing is assumed.

In the last three decades，there has been a continuous research interest in optimizing the flapping wing geometry and/or kinematics to achieve the optimal aerodynamic performance in terms of lift，lift to drag ratio，thrust，propulsive efficiency，power consumption，etc.，with different level of fidelity.Early in 1980s，Jones［10］had calculated the optimum load distribution along the wing for a given wing-root bending moment during the pure flapping motion based on classic aerodynamic theories.They found that the optimum loading was able to generate thrust more efficiently compared with an elliptic loading，which was suggested to be optimal for steady flight.After that，numerical optimization based on low fidelity aerodynamic models is carried out to optimize the hovering or forward flapping flight.Most of these studies adopted blade element momentum(BEM)theories［9，11-14］or vortex lattice methods［15-19］to model aerodynamic loadings.A trade-off of these low fidelity aerodynamic models is the low computational cost.Although these researches have made substantial progress in optimal flapping flight，low fidelity aerodynamic models may fail to provide appropriate level of complexity of the flow information when viscous effects became non-negligible［17］or the 3D nonlinear effects dominate the flow fields［20］.

The growing sophistication of computational hardware and simulation tools now allows for the direct coupling between optimization algorithms and high-fidelity Navier-Stokes solvers for the purpose of flapping motion optimization.Tuncer and Kaya［21］optimized the weighted sum consisting of thrust and propulsive efficiency of a flapping airfoil in combined plunging and pitching motion by coupling a gradient-based optimization method with a 2D Navier-Stokes solver supporting moving overset grids.Willis et al.［20］exploited a collection of computational tools with multiple fidelity levels in the design of effective flapping wing vehicles. In 2009， Pesavento and Wang［22］further confirmed that 2D flapping flight can be aerodynamically more efficient than the optimal steady motion.Recently，Culbreth et al.［23］presented several optimization results for flapping airfoils and wings obtained from an optimization framework constructed by a gradientbased optimization algorithm and a low-dissipation kinetic energy preserving(KEP)finite volume Navier-Stokes solver.They found that the maximum propulsive efficiency appeared to occur right before the incipience of leading edge separation.Jones and Yamaleev［24］performed shape and kinematics optimization of flapping wings for maximum propulsive efficiency using the 3D unsteady Reynolds averaged Navier-Stokes simulation.It is found that the combined optimization of wing shape and kinematics can significantly improve the performance of thrust and propulsive efficiency.

In the present study，we focus on maximizing endurance time for long-duration cruising flight.Therefore，the propulsive efficiency is chosen as the objective function.Optimization studies［11，13，16，19，21］for forward flight have been carried out with the vortex lattice method or low-order(≤2)CFD methods.Considering that flapping wing aerodynamics features vortex-dominated flow，we propose to take a high-order accurate CFD method based cost-effective high-fidelity optimization approach to provide further insights into flapping motions with maximum propulsive efficiency.At the first stage，a gradient-based optimization algorithm is coupled with a high-order spectral difference(SD)Navier-Stokes solver to search for the optimal kinematics on a coarse mesh.Then a high-order SD Navier-Stokes solver is used to simulate the flow fields with the optimized flapping kinematics on a fine mesh.Note that grid and polynomial order refinement studies are performed to make sure that the optimal kinematics obtained on the coarse mesh is still optimal on the fine mesh.Results from the refinement study will be presented in Section 2.Simulations on fine meshes can better resolve the detailed vortex structures associated with the optimal kinematics and help understand flow physics.

The remainder of this paper is organized as follows.In Section 1，the optimization framework including the numerical method is introduced.Numerical results are displayed and discussed in Section 2. Section 3 briefly concludes the present work and points out possible future work.

1 Optimization framework

1.1 Governing equations

Numerical simulations are performed with an unsteady compressible Navier-Stokes solver using dynamic unstructured grids based high-order spectral difference(SD)method developed in Ref.［25］.The 2D unsteady compressible Navier-Stokes equations in conservation form read，

where Q=(ρ，ρu，ρv，E)Tare the conserved variables，and F，G are the total fluxes including both the inviscid and viscous flux vectors，i.e.，F=Fi－Fvand G=Gi－Gv，which takes the following form

In Eq.(2)，ρ is the fluid density，u and v are the Cartesian velocity components，p is the pressure，and E is the total energy，μ is dynamic viscosity，Cpis the specific heat at constant pressure，Pr is the Prandtl number，and T is the temperature.The stress tensors in Eq.(2)take the following form

On assuming that the perfect gas law is obeyed，the pressure is related to the total initial energy by E=)with constantγ，which closes the solution system.

To achieve an efficient implementation，a time-dependent coordinate transformation from the physical domain(t，x，y)to the computational domain(τ，ξ，η)，as shown in Fig.1(a)，is applied to Eq.(1).And we obtain

Herein，τ=t and(ξ，η)∈［－1，1］2，are the local coordinates in the computational domain.In the transformation shown above，the Jacobian matrix J takes the following form

Note that the grid velocity vg=(xτ，yτ)is related with(ξt，ηt)by

Fig.1 Introduction of the spectral difference method

1.2 Space discretization

The SD method is used for the space discretization.In the SD method，two sets of points are given，namely the solution and flux points，as shown in Fig.1(b).Conservative variables are defined at the solution points(SPs)，and then interpolated to flux points to calculate local fluxes.In the present study，the solution points are chosen as the Chebyshev-Gauss quadrature points and the flux points are selected to be the Legendre-Gauss points with end points as-1 and 1.Then using Lagrange polynomials，we reconstruct all the fluxes at the flux points.Note that this reconstruction is continuous within a standard element，but discontinuous on the cell interfaces. To ensure conservation， the AUSM+-up approximate Riemann solver［26］is used to reconstruct a common inviscid flux on the interface.The reconstruction of the common viscous flux is based on the average of the‘left’and ‘right’fluxes.Readers can refer to［27］for more detailed implementation procedures.

1.3 Simulation parameters

In this study，the kinematics of the airfoil is specified as follows:

Plunging motion:y=h sin(ωt);

Pitching motion:θ=A sin(ωt+φ)

where y and θ are the plunging displacement and pitching angle of the airfoil respectively，h and A are the plunging and pitching amplitude respectively，ω is the angular frequency of the oscillation and φ is the phase lag between the plunging and pitching motions.

The aerodynamic parameters are specified as follows.The Reynolds number(Re)based on the airfoil chord length and the free stream velocity is defined as Re=ρU∞C/μ，well within the insect flight regime.The Strouhal number Str based on the plunging amplitude h is defined as Str=2fh/U∞，where f=ω/(2π)is the oscillating frequency.The reduced frequency k based on the airfoil chord length C is defined as k=2πfC/U∞.The free stream Mach number(Ma)in this study is fixed at Ma=0.1.The AUSM+-up approximate Riemann solver for all speeds［26］is used to resolve the low Mach number flow.

1.4 Optimization procedure

For simplicity，we consider an unconstrained optimization problem of the form

where ηpoweris the propulsive efficiency，〈CT〉is the time averaged thrust coefficient，〈Cpower〉is the time averaged power coefficient of the oscillating airfoil，and x denotes the optimization parameters.For 2D simulations，the time averaged thrust coefficient〈CT〉is defined asThrust(t)/(0.5ρ∞C)d t，where C is the chord length.The rule-of-thumb definition of time averaged power coefficient〈Cpower〉in numerical simulations iswhere θ is the pitch angle of the airfoil.The unsteady flow solver is used to evaluate the propulsive efficiency

The Broyden-Fletcher-Goldfarb-Shanno(BFGS)method［28］is adopted for the optimization.It is a quasi-Newton algorithm and can achieve a superlinear rate of convergence.The BFGS algorithm is briefly reviewed here for completeness［28］.

Denote the function to be minimized by f(f=ηpowerin our case).Givean educated guess of the optimization parameters as x0.Let the convergence tolerance be ε＞0，s0=x1－x0，and y0=▽f1－▽f2.Then，the initial inverse Hessian approximation H0can be specified as below:

Compute line search direction

Set xk+1=xk+αkpk，where αkis the line search step length，which satisfies the Wolfe condition，

Update sk， yk， and Hk+1using the following formulas:

Typical values for c1and c2in the Wolfe condition suggested in Ref.［28］are c1=1×10-4and c2=0.9.In the present optimization，a relatively large line search step length is used，namely αk=1.0.If the Wolfe condition is not satisfied，then a contraction ratio 0.1 is multiplied to αkto limit the step size.The convergence criteria ε for‖▽η‖2is set as 10-5.Based on our previous studies［29］，although sometimes ‖▽η‖2does not converge to the desired error level(e.g.，‖▽η‖2oscillates around 10-2)，yet the values of the optimization parameters converge almost surely towards the local extrema.Under such circumstances，another criterion is augmented to judge convergence.It states that if the differences of the absolute value of η in successive optimization steps are smaller than 10-3，then the optimization process will stop and the current optimization parameters are deemed as optimal.

2 Optimization results and discussions

Optimization variables are flapping kinematics parameters，including the phase shift between plunge and pitch motions φ，the pitch and plunge amplitudes A，h and the oscillation angular frequency ω.Since the aerodynamic performance can vary significantly with airfoils thickness［27，30］， three NACA 4-digit airfoils，namely，NACA0006，NACA0012 and NACA0030，in combined plunge and pitch motion are used to test the effects of airfoil thickness on optimal flapping kinematics.The effects of Reynolds number and pitching pivot location on optimal flapping kinematics are also presented in this section.

2.1 Numerical verification

First of all，the optimization framework is verified using the NACA0012 airfoil.The optimization parameter is the phase shift between plunge and pitch motions φ.Initially，φ is set as 70°.Other parameters are specified as Str=0.3，k=1.0，A=20°and Re=1，200.Both hrefinement(grid refinement)and p-refinement(polynomial degree)studies are carried out to determine the suitable grid and scheme accuracy.The results are displayed in Table 1.For the coarse mesh，34 elements are placed on the airfoil surface and 23 elements are assigned in the radial direction from the airfoil.For the medium mesh，the element number in the radial direction is kept the same as that for the coarse mesh，while the element number in the circumferential direction(along the airfoil surface)is 70 elements.For the fine mesh，the element number in the circumferential direction is 100，and that in the radial direction is 35.The minimum cell size normalized by the airfoil chord length in the radial direction is 0.001 for all cases，which has been demonstrated small enough to resolve boundary layer flow features for Reynolds numbers of order 104when P2reconstruction is used［25］.From Table 1，it is found that the results(i.e.，η and φ)using P3reconstruction on the coarse mesh agree well with those using P2reconstruction on both the coarse and medium meshes.The results using P1reconstruction are converging towards those values.To save computational cost，P2reconstruction on the coarse mesh is used in all optimization simulations.A thorough verification of the current Navier-Stokes flow solver for unsteady flapping wing aerodynamics simulations can be found in Ref.［25，27］.

Table 1 Optimization results for the oscillating NACA0012 airfoil from h-refinement and p-refinement studies at Str=0.3，k=1.0，A=20°and Re=1，200

The convergence histories of the L2norm of ▽η and φ using P2reconstruction on the coarse grid are shown in Fig.2.We observe that the optimal propulsive efficiency is achieved after only six optimization iterations.From the convergence path，it is clear that the propulsive efficiency increases about 5% as φ increases from 70°to 87.3°.

Fig.2 Convergences of the L2 norm of ▽η and the phase lag φ

2.2 Optimization results

Since the gradient-based optimization algorithm is used to search for the local extrema，initial guess needs to be sufficiently close to the optimal point to ensure a solution.Thus，the initial guess must be based on the flow physics of flapping flight.In this paper，initial values of the optimization variables are chosen based on the previous study［27］.After optimization，in order to elucidate the detailed vortical structures under the optimal airfoil oscillation，flow fields associated with the optimized kinematics are simulated using a 3rdorder accurate SD scheme on the fine mesh.

2.2.1 Effects of airfoil thickness

The optimization results of flapping kinematics for NACA0006，NACA0012 and NACA0030 airfoils are presented in Table 2.The Reynolds number is chosen as 1，200，and the pitch pivot is located at C/3 from the leading edge of the airfoil.From this table，it is found that the optimized flapping kinematics and propulsive efficiency for NACA0006 and NACA0012 are very similar.For both airfoils，the propulsive efficiency is larger than 50%.The maximum propulsive efficiency for the oscillating NACA0030 airfoil is about 32%，and the optimized flapping kinematics is very different from that of NACA0006 and NACA0012.These results indicate that airfoil thickness has large influence on the optimal flapping kinematics，and relatively thin airfoils have more potential to acquire superior aerodynamic performance.This trend agrees well with the observation from previous studies［31，27］.

Table 2 Optimization results for a series of NACA 4-digit airfoils with Re=1，200 and the pitch pivot located at C/3 from the leading edge of the airfoil

The vorticity fields around the oscillating NACA0006，NACA0012 and NACA0030 airfoilswith the optimized kinematics at four phases，namely，0，π/2，π and 3π/2，are displayed in Fig.3 ～ Fig.5，respectively.It is found that large scale leading edge separation only occurs when the airfoil starts to reverse its plunging direction.During this short oscillation reversion period，the vortices are quickly shed from the airfoil and the flow then attaches to the airfoil surface again.During most of the oscillation cycle，no leading edge separation shows up.These flow features may be critical for maintaining high propulsive efficiency，which has been reported in［21，23］.This will be further discussed in the Section 2.3.We also observe that the flow fields for NACA0006 and NACA0012 airfoils are very similar.The separation region on the NACA0030 airfoil at phasesπ/2 and 3π/2 is much larger than that on the thin airfoils.This will affect the aerodynamic force generation.

Fig.3 Vorticity fields with optimized kinematics at different phases for the NACA0006 airfoil

Fig.4 Vorticity fields with optimized kinematics at different phases for the NACA0012 airfoil

Fig.5 Vorticity fields with optimized kinematics at different phases for the NACA0030 airfoil

The time histories of the thrust and lift force coefficients for all three airfoils are compared in Fig.6.From the results it is observed that both the thrust and lift force histories of NACA0006 and NACA0012 airfoils are very similar，while that of the NACA0030 airfoil exhibits larger peaks.Especially for the lift force history，the peak lift value on NACA0030 near the phases π/2 and 3π/2 is more than twice of that on both NACA0006 and NACA0012.This corresponds well with the flow fields discussed previously.

Furthermore，these phenomena explain why the optimized propulsive efficiencies for NACA0006 and NACA0012 are much larger than that for NACA0030 as shown in Table 2.Specifically，according to the definition of the propulsive efficiency η，if the thrust coefficient〈CT〉is similar，then the propulsive efficiency is determined by the input power〈Cpower〉.If we assume that the mass and moment of inertia of the airfoil can be neglected(i.e.，no inertia)，then the input power is actually the work done by the lift force and the aerodynamic moment per unit time.It has been examined that the work done by the aerodynamic moment is much smaller than that done by the lift force.Therefore，it is clear that if the thrust coefficient〈CT〉is almost the same，the lift force in fact dominates the propulsive efficiency.Small lift generation is conducive to maintain high propulsive efficiency.The large lift generation on the NACA0030 airfoil results in inferior propulsion performance.

Fig.6 Time histories of the thrust and lift force coefficients with optimized kinematics for a series of NACA 4-digit airfoils with Re=1，200 and the pitch pivot located at C/3 from the leading edge of the airfoil

2.2.2 Effects of Reynolds number

The optimization results of flapping kinematics for NACA0012 under different Reynolds numbers are presented in Table 3.The pitch pivot is located at from the leading edge of the airfoil.Several observations are concluded from this table.First of all，larger propulsive efficiency is obtained as the Reynolds number increases.According to Anderson et al.［32］，the propulsive efficiency of an oscillating NACA0012 airfoil can reach as high as 87% in the Strouhal number range 0.3 ～0.4 at Re=40，000.In this study，the maximum propulsive efficiency of the NACA0012 airfoil at Re=1，200 is achieved at Str=0.42 with a value of 53.4%.The relatively low propulsive efficiency is probably due to the relatively low Reynolds number adopted in the study.Secondly，for the case with Re=300，the optimal Strouhal number is 0.78.This value is much larger than the optimal Strouhal number range 0.25 ～ 0.35 reported in［33］.This is due to that the analysis in［33］is based on high Reynolds number flow.Current results show similarity with the literature reported for low Reynolds number range(e.g.，100)［34］.

Table 3 Optimization results for a NACA0012 airfoil with different Reynolds numbers

The vorticity fields around the oscillating NACA0012 airfoil with the optimized kinematics at the first two phases，namely，0 and π/2，for both Reynolds numbers 300 and 600 are shown in Fig.7.Note that since anti-symmetric vortex patterns will show up for phases π and 3π/2，as illustrated in Fig.3 ～ Fig.5，only the vorticity fields at the first two phases are presented.Similarly as in the cases for different thicknesses，no leading edge separation shows up during most of the oscillation cycle;large scale leading edge separation only occurs when the airfoil starts to reverse its plunging direction.

Fig.7 Vorticity fields with optimized kinematics at phases 0 and π/2 for different Reynolds numbers

The time histories of the thrust and lift force coefficients for the three Reynolds numbers are compared in Fig.8.It is found that the aerodynamic force results of the cases with Re=600 and Re=1，2000 are similar，and depart from those of the case with Re=300.Besides，the detailed vortical structures at Re=300 are quite different from those at Re=600 and Re=1，200 as shown in Fig.4 and Fig.7.This is primarily due to the distinct differences presented in the optimal kinematics as shown in Table 3.

Fig.8 Time histories of the thrust and lift force coefficients with optimized kinematics for a NACA0012 airfoil with different Reynolds numbers(the pitch pivot is located at C/3 from the leading edge of the airfoil)

2.2.3 Effects of pivot location

The optimization results of flapping kinematics with different pitching pivot locations are presented in Table 4.Three frequently adopted pivot locations，in previous experimental and numerical studies［4］，namely，C/4，C/3 and C/2，are used in this section.The NACA0012 airfoil is adopted in this section，and the Reynolds number is fixed at 1，200.The vorticity fields around the oscillating NACA0012 airfoil with the optimized kinematics at the first two phases，namely，0 and π/2，for both pitching pivot locations at C/4 and C/2 are presented in Fig.9.The time histories of the thrust and lift force coefficients for the three different pitching pivot locations are compared in Fig.10.From Table 4 and Fig.10，it is found that the optimized flapping kinematics and the aerodynamic performance for all these cases are very similar.From Fig.9，it is observed that although the aredynamic performance of the cases with different pitching pivot locations is similar，the detailed vortical structures with pitching pivot located at C/4 are apparently different from those with pitching pivot located at C/3 and C/2.Based on all the results，we conclude that the locations of the three frequently adopted pitching pivots pitching pivot do not significantly affect the kinematics optimization.

Fig.9 Vorticity fields with optimized kinematics at phases 0 andπ/2 for different pitch pivot locations

Fig.10 Time histories of the thrust and lift force coefficients with optimized kinematics for a NACA0012 airfoil with different pitch pivot locations(the Reynolds numbers is fixed at 1，200)

Table 4 Optimization results for a NACA0012 airfoil with different pitch pivot locations(the Reynolds numbers is fixed at 1，200)

2.3 Discussion of the optimized kinematics

In this section，aerodynamic performance and the corresponding flow fields from the optimized and non-optimized cases around an oscillating NACA0006 airfoil are compared.From Table 5，it is clear that although the time averaged thrust coefficient is almost the same for these two cases，the propulsive efficiency for the optimized case is about 2.5 times of that for the non-optimized case.The time histories of the thrust and lift force coefficients with optimized and non-optimized kinematics for a NACA0006 airfoil with Re=1，200 and the pitch pivot located at C/3 are presented in Fig.11.The vorticity fields with the non-optimized kinematics parameters at four phases are displayed in Fig.12.It is obvious that large leading edge separation appears during the entire oscillation cycle.Also observed from Table 5，there exists apparent difference on the pitching amplitude between the optimized and non-optimized cases.This infers that the large leading edge separation from the non-optimized case could be connected with the variation of the effective angle of attack(AOA).

The definition of the effective AOA for the combined plunge and pitch motion is given as

assuming that the pitch speed is negligible.Sketch of the effective AOA during the oscillation is shown in Fig.13.The variation of the effective AOA for both optimized and non-optimized cases of NACA0006 is displayed in Fig.14(a).It is obvious that for the optimized case，the peak absolute value of effective AOA is well controlled below 15°.This feature can help prevent the appearanceof large leading edge separation.However for the nonoptimized case，the peak absolute value of effective AOA can reach about 40°and during the entire oscillation cycle，the effective AOA remains at large value.Then a large leading edge separation appears during the entire oscillation cycle for the non-optimized case.

Table 5 Comparison of kinematics parameters and aerodynamic performance of a NACA0006 airfoil with optimized and non-optimized kinematics parameters(the Reynolds number is fixed at Re=1，200，and the pitch pivot is located at C/3 from the leading edge of the airfoil)

Fig.11 Time histories of the thrust and lift force coefficients with optimized and non-optimized kinematics for a NACA0006 airfoil with Re=1，200 and the pitch pivot located at C/3 from the leading edge of the airfoil

Fig.12 The vorticity fields with non-optimized kinematics at different phases for the NACA0006 airfoil at Re=1，200(the pitch pivot is located at C/3 from the leading edge of the airfoil.Other parameters are tabulated in Table 5)

Fig.13 Sketch of the effective angle of attack during the combined plunge and pitch motion

Fig.14 Effective angle of attack

The variation of the effective AOA for different airfoils，Reynolds numbers and pivot locations with optimized kinematics parameters is shown in Fig.13(b)～ Fig.13(d).For the airfoil thickness tests，it is clear that the variation trend for both NACA0006 and NACA0012 is quite similar，but different from that for NACA0030.It is also found that the peak absolute value of effective AOA of NACA0030 is about 30°，which is twice that of the other two airfoils.Thus the lift generation on NACA0030 can be much larger than that on NACA0006 and NACA0012.Based on the analyses in Section 2.2.1，the lower optimized propulsive efficiency of NACA0030 can be attributed to the large variation of the effective AOA.For both the Reynolds number and pitching pivot location tests，it is found that the peak absolute values of effective AOA are well controlled below 20°for all these cases.

3 Conclusions and future work

A high-order SD method based numerical optimization approach is developed to optimize the kinematics of flapping airfoils and provide further insights into flapping motions with maximum propulsive efficiency.As a first step，a gradient-based optimization algorithm is coupled with a high-order SD Navier-Stokes solver to search for the optimal kinematics for a certain airfoil undergoing combined plunging and pitching motion on a coarse mesh. Specifically， the quasi-Newton optimization algorithm BFGS is used to speed up the convergence.Then the high-order SD scheme is used to simulate the detailed vortex structures associated with the optimal kinematics of the flapping flight on a fine mesh.

The proposed optimization procedure isused to optimize several flapping kinematics，including the oscillation frequency，plunging and pitching amplitudes and the phase lag between the plunging and pitching motion.The effects of airfoil thickness，Reynolds number，and pivot location on the optimization output are investigated.It is found that the maximum propulsive efficiency for NACA0006 and NACA0012 airfoils is more than 50% after four kinematics parameters are optimized.However，the optimized propulsive efficiency for a NACA0030 airfoil can only reach about 30%.As Reynolds numbers increase，the propulsive efficiency becomes larger.Itis also found that locations of pitching pivot do not significantly affect the kinematics optimization.Based on the comparison between optimized and non-optimized cases for the NACA0006 airfoil，it is found that the variation of effective AOA plays a key role in achieving high propulsive efficiency.Specifically，under optimized conditions，the variation of effective AOA is well controlled to avoid the occurrence of large leading edge separation during most of the oscillation cycle.

We note that this work is a first attempt to incorporate high-order accurate numerical schemes into the optimization framework for unsteady vortex-dominated flow.In future work，more sophisticated flapping flight conditions will be considered，such as including the constraints of lift and power.In order to extend the current work for efficient 3D optimization of both kinematics and geometry，an adjoint-based numerical optimization framework needs to be established.These issues are reserved for future publications.

［1］HO S，NASSEF H，PORNSINSIRIRAK N，et al.Unsteady aerodynamics and flow control for flapping wing flyers［J］.Progress in Aerospace Sciences，2003，39(8):635-681.

［2］PLATZER M F，JONES K D，YOUNG J，et al.Flapping-wing aerodynamics:progress and challenges［J］.AIAA Journal，2008，46(9):2136-2149.

［3］ROZHDESTVENSKY K V，RYZHOV V A.Aerohydrodynamics of flapping-wing propulsors［J］.Progress in Aerospace Sciences，2003，39(8):585-633.

［4］SHYY W，AONO H，CHIMAKURTHI S，et al.Recent progress in flapping wing aerodynamics and aeroelasticity［J］.Progress in Aerospace Sciences，2010，48(7):284-327.

［5］SHYY W，BERG M，LJUNGQVIST D.Flapping and flexible wings for biological and micro air vehicles［J］.Progress in Aerospace Sciences，1999，355:455-505.

［6］SHYY W，LIAN Y，TANG J，et al.Aerodynamics of low Reynolds number flyers［M］.New York:Cambridge Univ.Press，2008.

［7］TRIANTAFLYLLOU M S，TECHET A H，HOVER F S.Review of experimental work in biomimetic foils［J］.IEEE Journal of Oceanic Engineering，2004，293:585-594.

［8］WANG Z J.Dissecting insect flight［J］.Annual Review of Fluid Mechanics，2005，371:183-210.

［9］BERMAN G J，WANG Z J.Energy-minimizing kinematics in hovering insect flight［J］.J.Fluid Mech.，2007，582:153-168.

［10］JONES R T.Wing flapping with minimum energy［R］.NASA Technical Memorandum 81174，1980.

［11］HAMDAOUI M，MOURET J B，DONCIEUX S，et al，Optimization of kinematics for birds and UAVs using evolutionary algorithms［J］.World Academy of Science，Engineering and Technology，2008，47:181-192.

［12］STANFORD B，KURDI M，BERAN P，et al.Shape，structure，and kinematic parameterization of a power-optimal hovering wing［R］.In:51st AIAA/ASME/ASCE/AHS/ASC Structures，Structural Dynamics， and Materials Conference. Orlando， Florida，2010.AIAA-2010-2963.

［13］GOGULAPATI A，FRIEDMANN P P，MARTINS J R R A.Optimization of the kinematics of a flapping wing MAV in hover for enhanced performance［R］.In:54th AIAA/ASME/ASCE/AHS/ASC Structures，Structural Dynamics，and Materials Conference.Boston，Massachusetts，2013.2013-1646.

［14］WANG Q，GOOSEN H，KEULEN F V.Optimal hovering kinematics with respect to various flapping-wing shapes［C］.In:International Micro Air Vehicle Conference and Flight Competition(IMAV2013).Toulouse France，2013.

［15］HALL K C，HALL S R.Minimum induced power requirements for flapping flight［J］.Journal of Fluid Mechanics，1996，323:285-315.

［16］ITO K.Optimization of flapping wing motion［C］.In:ICAS 2002 CONGRESS，2002.

［17］STANFORD B K，BERAN P S.Analytical sensitivity analysis of an unsteady vortex-lattice method for flapping-wing optimization［J］.Journal of Aircraft，2010，472:647-662.

［18］GHOMMEM M，HAJJ M R，WATSON L T，et al.Deterministic global optimization of flapping wing motion for micro air vehicles［C］.In:2010 AIAA ATIO/ISSMO Conference，2010.

［19］GHOMMEM M， COLLIER N， NIEMI A H，et al.On the shape optimization of flapping wings and their performance analysis［J］.Aerospace Science and Technology，2014，32:274-292.

［20］WILLIS D J，PERSSON P O，ISRAELI E R，et al.Multifidelity approaches for the computational analysis and design of effective flapping wing vehicles［C］.In:46th AIAA Aerospace Sciences Meeting and Exhibit.Reno，Nevada，2008.AIAA-2008-518.

［21］TUNCER I H，KAYA M.Optimization of flapping airfoils for maximum thrust and propulsive efficiency［J］.AIAA Journal，2005，4311:2329-2336.

［22］PESAVENTO U，WANG Z J.Flapping wing flight can save aerodynamic power compared to steady flight［J］.Physical Review Letters，2009，103:118102-4.

［23］CULBRETH M，ALLANEAU Y，JAMESON A.High-fidelity optimization of flapping airfoils and wings［R］.In:29th AIAA Applied Aerodynamics Conference.Honolulu，Hawaii，2011.AIAA-2011-3521.

［24］JONES M，YAMALEEV N K.Adjoint-based shape and kinematics optimization of flapping wing propulsive efficiency［R］.In 43rd AIAA Fluid Dynamics Conference.San Diego，CA，2013.AIAA-2013-2472.

［25］YU M L，WANG Z J，HU H.A high-order spectral difference method for unstructured dynamic grids［J］.Computer ＆ Fluids，2011，48:84-97.

［26］LIOU M S.A sequel to AUSM，Part II:AUSM+-up for all speeds［J］.Journal of Computational Physics，2006，214:137-170.

［27］YU M L，WANG Z J，HU H.High fidelity numerical simulation of airfoil thickness and kinematics effects on flapping airfoil propulsion［J］.Journal of Fluids and Structures，2013，42:166-186.

［28］NOCEDAL J，WRIGHT S.Numerical optimization，Springer，2000.

［29］YU M L，WANG Z J，HU H.High-fidelity optimization offlapping airfoils for maximum propulsive efficiency［R］.In:51st AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition，Grapevine(Dallas/Ft.Worth Region).Texas，2013.AIAA-2013-0085

［30］YU M L，WANG Z J，HU H.Airfoil thickness effects on the thrust generation of plunging airfoils［J］.Journal of Aircraft，2012，495:1434-1439.

［31］ASHRAF M A，YOUNG J，LAI J C S.Reynolds number，thickness and camber effects on flapping airfoil propulsion［J］.Journal of Fluids and Structures，2011，27:145-160.

［32］ANDERSON J M，STREITLIEN K，BARRETT D S，etal.Oscillating foils of high propulsive efficiency［J］.Journal of Fluid Mechamics，1998，360:41-72.

［33］TRIANTAFLYLLOU G S，TRIANTAFLYLLOU M S，GROSENBAUGH M A.Optimal thrust development in oscillating foils with application to fish propulsion［J］.Journal of Fluids and Structures，1993，7:205-224.

［34］ZHANG L.Numerical studies of hydrodynamics of fish locomotion and schooling by a vortex particle method［R］.ProQuest，2008.［35］YOUNG J，LAI J C S.Oscillation frequency and amplitude effects on the wake of a plunging airfoil［J］.AIAA Journal，2004，42:2042-2052.