Nomenclature w, 　state vector; T, 　time period;
f_ci, 　inviscid flux vector; t, 　time;
F_ci, 　Cartesian inviscid flux vector; t_f , final time;
f_vi, 　viscous flux vector; t_n, 　nondimensional time;
F_vi, 　Cartesian viscous flux vector; V , cell volume;
u_i, 　velocity components of the fluid; Pr, 　Prandtl number;
u_bi, 　velocity components of the boundary; δ_ij , Kronecker delta;
k, 　coefficient of thermal conductivity; i_var, the sequence number of design variables;
E, 　total energy; K_ij , 　transformation metrics;
p, 　pressure; C_L, 　lift coefficient;
σ_ij , 　viscous stresses; C_D, 　drag coefficient;
ξ_i, 　computational coordinates; I, 　objective function;
R, 　residual; Ψ, Lagrange multiplier;
S, 　shape function; G, 　local gradient. 1 Introduction

With the development of computational fluid dynamics,the theory of airfoil optimization has been maturing over the past few decades^[1].Nowadays,main airfoil optimization approaches include direct search methods,stochastic methods,and gradient-based methods.Among all these methods,the adjoint method,which was proposed by Jameson^{[2, 3]}in the 1980s,is most widely adopted.Since the computational expense incurred in the calculation of the gradient is independent of the number of design variables,the adjoint approach can reduce the cost of the computation.The adjoint method can reduce the cost of the computation since the computational expense incurred in the calculation of the complete gradient is effectively independent of the number of design variables.During the last two decades,this method has been developed rapidly.Reuther and Jameson^[4],Jameson et al.^[5],and Reuther et al.^{[6, 7]}carried out studies investigating the design of airfoils,wings,and the whole aircraft using adjoint methods.Anderson and Venkatakrishnan^[8]and Bueno-Orovio et al.^[9]applied this approach on unstructured grids.

However,most of the studies focused on the analysis of steady flows.While many flow fields which we are interested in,such as the flow fields around the rotorcraft,turbo machinery,and flapping flight,are unsteady in nature.At present,there are only few publications that have addressed the optimization of unsteady airfoils.In recent years,some researchers used the adjoint method in the design of unsteady airfoils.Most of the works were based on the Euler equations considering the difficulties encountered in dealing with the viscosity term.Nadarajah and Jameson^{[10, 11]}derived the time-accurate adjoint equations with the Euler equations.Srinath and Mittal^[12]presented a new method that combined the finite element method and the continuous adjoint method to optimize plunging airfoils.

Traditionally,the optimization approach which has been adopted for the unsteady optimization of airfoil profiles only requires a small extension of a steady flow design code.In the method studied before,the steady flow solutions are computed for a various number of cases by varying freestream conditions to indicate the unsteady flow solutions.Compared with this approach,the method presented in this paper is a full unsteady optimization method.In the present work,a method based on the continuous adjoint method is formulated for the optimization of unsteady airfoils in viscous flows.A basic framework of time-accurate unsteady airfoils optimization using the time-averaged aerodynamic coefficients as the objective functions is presented,and the time-accurate continuous adjoint equation and its boundary conditions are derived.Then,we solve both the flow field governing equations and the continuous adjoint equations numerically using the finite volume method (FVM) formulation.The multi-grid method is adopted during the calculations to accelerate the convergence.The numerical results are validated by comparison with the direct numerical simulate (DNS) data.Finally,the plunging NACA0012 airfoil is adopted to verify the optimization method.The Hicks-Henne bump functions are used to modify the shape of the airfoil.Two different objective coefficients,i.e.,the maximize time-averaged lift coefficient and the minimize time-averaged drag coefficient,are used to obtain optimal shapes in unsteady viscous flows.

2 Governing equations 2.1 Arbitrary Lagrangian-Eulerian (ALE) formulations of Navier-Stokes equations

In order to simulate the movement of the airfoils,the ALE formulations of Navier-Stokes equations are adopted as the governing equations,i.e.,

where the state vector w,the inviscid flux vector f_ci,and the viscous flux vector f_vi are described by

In the above definitions,u_i and u_bi are the Cartesian velocity components of the fluid and the boundary,respectively.E is the total energy,and k is the coefficient of thermal conductivity.The pressure is determined by the equation of state

The viscous stresses are written as

where μ is the sum of the laminar coefficients of viscosity and the eddy viscosity,and λ is the second coefficient of viscosity.In this paper,the flow field is treated as laminar. 2.2 Numerical simulation methods

The methodology which is used here is the FVM.When using a discretization on a bodyconforming structured mesh,we adopt a transformation to the computational coordinates ξ_i,which is defined by the matrices as follows:

Then,the governing equations can be written in the computational space as where

V is the cell volume,and F_ci=S_ijf_cj,F_vi=S_ijf_vj,and S_ij=JK_cj⁻¹ .

As mentioned above,the unsteady flow field is solved numerically using the FVM formulation.The governing equations are spatially discretized using the Jameson-Schmidt-Turkel (JST) scheme.Time integration is handled by a second-order accurate dual-time stepping approach.The multi-grid method is adopted during the calculations to accelerate the convergence.

In order to validate the numerical method that we have adopted,the unsteady viscous flow field around the plunging NACA0012 airfoil is calculated.Simulations are finished on a structured "C-mesh" counting 240×121 cells.The computational domain extends roughly 20 chord lengths downstream and 10 chord lengths upstream,as seen in Fig. 1.The mesh is subject to the rigid body motion and moves with the airfoil.The Reynolds number is 1 850,the Mach number is 0.2,and the Strouhal number is 0.6.Figure 2 contains a comparison of the lift and drag coefficients versus time between the computational results and the DNS results^[13].The computational results agree well with the DNS results,which means that this numerical method can obtain accurate unsteady airfoil aerodynamic parameters effectively.

3 Time-accurate continuous adjoint method 3.1 Continuous unsteady adjoint equations

The cost function in this problem is defined as the function of the flow field variable,w,and the physical location of the boundary,which can be represented by the function,S.The concrete form of the objective function is given as follows:

Here,a change in S results in a change of the cost function.

Then,introduce the governing equations as a constraint

We have

Combining the items that contain the flow filed information yields

Now,we have the unsteady continuous adjoint equation as

The boundary condition is The variation of the objective function is

3.2 Numerical simulation methods

Comparing (2) with (8),we find that the forms of the continuous adjoint equation and the flow field governing equation are almost consistent.Therefore,we can use the same numerical scheme which is used to the unsteady flow field to solve the adjoint equation.It is important to note that the unsteady adjoint equation requires integration in reverse time.Therefore,we have to make a minor adjustment in the computational code.

4 Design optimization of plunging airfoils 4.1 Implementation of optimization procedure

The unsteady design process in this paper is shown in Fig. 3.First,calculate the unsteady flow with the numerical methods mentioned above and save the flow solution at each time step during the calculation.Second,solve the adjoint equations with the flow variables obtained from the first step.Then,calculate the gradient information using the numerical solutions.Then,modify the airfoil shape in the direction of improvement using a simple descent method.The whole design process is repeated until the optimal design results meet the requirements.

4.2 Minimization of time-averaged drag coefficient

An NACA0012 airfoil is selected as the initial airfoil.Here,the computations are performed at a free stream where the Mach number is 0.2,and the Reynolds number is 1 850.The Strouhal number of the plunging airfoil is 0.2.There are 38 points up and down along the surface of the airfoil chosen as the design variables.The Hicks-Henne bump functions are used,which can be added to the original airfoil geometry to modify the shape.The time-averaged drag coefficient is chosen as the objective function.Meanwhile,the time-averaged lift coefficient is chosen as the constraint condition.Figure 4 presents the gradients of the objective function at the end of the first calculate cycle.After 25 design cycles,the solution converges.Figure 5 presents the optimization result of the minimized time-averaged drag.According to the solution,we find that the optimized airfoil overall thickness decreases,and the drag coefficient reduces from 0.063 1 to 0.049 4,while the lift coefficient is kept at 0.

4.3 Maximization of time-averaged lift coefficient

In this case,the time-averaged lift coefficient is chosen as the objective function,while the time-averaged drag coefficient is chosen as the constraint condition.After 29 design cycles,the solution converges.Figure 6 presents the optimization result of maximization of the timeaveraged lift coefficient.According to the calculation results,we find that the optimized airfoil overall camber increases,and the lift coefficient increases from 0 to 0.440 89,while the drag coefficient maintains at 0.06.

5 Conclusions

This article presents a complete formulation of the continuous unsteady viscous adjoint approach to aerodynamic optimal design with both theoretical derivation and case validation.

A significant reduction in the time-averaged drag coefficient is achieved for the plunging NACA0012 airfoil while maintaining the time-averaged lift coefficient.The time-averaged lift coefficient increases a lot with the fixed time-averaged drag coefficient.These results demonstrate that this method can effectively solve the problem of unsteady airfoil optimization design as an ideal optimization approach with low cost and high precision.