A surrogate-model based multidisciplinary shape optimization method with application to a 2D subsonic airfoil

Computers & Fluids 36 (2007) 520–529 www.elsevier.com/locate/compfluid

A surrogate-model based multidisciplinary shape optimization method with application to a 2D subsonic airfoil J.-C. Jouhaud a

a,*

, P. Sagaut b, M. Montagnac a, J. Laurenceau

a

CERFACS, European Center For Research and Formation in Advanced Computations, 42 Avenue Coriolis, 31057 Toulouse Cedex 01, France b Laboratoire de Mode´lisation en Me´canique, Universite´ Pierre et Marie Curie, Paris VI, 4 Place Jussieu, 75252 Paris Cedex 5, France Received 6 February 2005; received in revised form 1 March 2006; accepted 4 April 2006 Available online 22 August 2006

Abstract A surrogate-model based shape optimization method is presented and applied to the case of the multidisciplinary shape optimization of a 2D NACA subsonic airfoil. The cost function is designed so that both the far-field radiated noise and the aerodynamic forces are controlled. The surrogate model is based on the Kriging optimal interpolation technique. In order to increase the efficiency of the method, a dynamic Kriging method is developed, which can be interpreted as an Adaptive Mesh Refinement method in the shape optimization parameters.  2006 Elsevier Ltd. All rights reserved.

1. Introduction Shape optimization is the final purpose of most engineering works dealing with Computational Fluid Dynamics [1,2]. Despite a huge amount of work has been devoted to the design of efficient optimization techniques, the definition of a robust, general-purpose reliable technique for multidisciplinary shape optimization is still a research topic. Several questions must be addressed to build such a methodology. First, the cost function to be minimized is a priori not globally convex in the general case, leading to the use of non-gradient-based techniques able to capture the existence of multiple minima and to identify the global minimum. The systematic exploration of the space of solution spanned by the shape optimization parameters is not realistic in practical cases, because it requires an enormous number of evaluation of the solution and leads to an amazingly large computational cost. This cost can be reduced using some smart techniques such as genetic algorithms [3,5,6,8,9] or pattern identification algorithms, but these

*

Corresponding author. Tel.: +33 561 19 30 51; fax: +33 561 19 30 00. E-mail address: [email protected] (J.-C. Jouhaud).

0045-7930/$ - see front matter  2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2006.04.001

methods also induced a large computational cost since they require an initial exploration of the space of solutions associated with the optimization problem. Therefore, it is important to reduce the computational cost while using a method well suited for globally non-convex cost function. A key idea is to parameterize the space of possible solutions via a simple, computationally unexpensive model, and to use this model to generate inputs for the optimization algorithm. Such a model is often referred to as the response surface of the system to be optimized, leading to the definition of a so-called surrogate-model based optimization methodology [4,7]. Another expected feature of an optimization methodology is generality: different users will be interested in different problems, and the method must be well suited for multipurpose, multidisciplinary optimization. It is advocated here that such an optimization tool should be an external module coupled to the CFD tool, but should not be developed inside the CFD solver for the sake of generality and flexibility. The aim of the present paper is to assess the potential use of the Kriging method to generate a reliable response surface for multidisciplinary shape optimization of a twodimensional wing in a turbulent flow. The cost function under consideration here includes both aerodynamic and

J.-C. Jouhaud et al. / Computers & Fluids 36 (2007) 520–529

aeroacoustic criteria. The use of the Kriging method for surrogate-model based shape optimization was considered in [1,2] in the case of a low-Reynolds, incompressible flow past a trailing edge using Large-Eddy Simulation. We address here the problem of the high-Reynolds number, low-subsonic flow around a non-symmetric NACA profile, the turbulent flow being computed using a RANS-type approach. Since only the mean flow is computed, a trailing-edge noise metric is used to evaluate the radiated noise. In order to improve the accuracy of the response surface, a dynamic multiresolution technique in the optimization parameter space is developed, which makes it possible to minimize the computational cost of the response surface building step. The article is organized as follows. Governing equations, including turbulence modeling and main features of the numerical method, are presented in Section 2. The Kriging method and its present dynamic multilevel implementation are detailed in Section 3. The selected test case and the numerical results are discussed in Section 4.

521

where l0 is the dynamic viscosity at the reference temperature T0 and the constant Cs equals to 110.3 K. With a constant Prandtl number, the heat conductivity can be written lC as K T ¼ Prp with Cp the specific heat at constant pressure and Pr = 0.72 for air. For a caloric perfect gas, the state equation is given by p = qRT where the gas constant R is equal to 287 (J/kg K) for air. 2.2. Wilcox k  w turbulence model The k  x model equations [11] are given by the system oqk ~ Þ  r:ððl þ rk lt ÞrkÞ ¼ P k  bk qxk; þ r:ðqk U ot oqx ~ Þ  r:ððl þ rx lt ÞrxÞ ¼ ax x P k  bx qx2 þ r:ðqxU ot k ð6Þ ~ , the eddy viscosity with the production term P k ¼ s:rU qk lt ¼ x , the closure coefficients rk = 0.5, rx = 0.5, bk = 0.09, bx = 0.075, ax = 0.5, k is the turbulent kinetic energy and x is the specific turbulent dissipation.

2. Governing equations and flow solver 2.3. Flow solver—elsA software 2.1. Navier–Stokes equations for compressible fluids The governing equations are the 3D Navier–Stokes equations which describe the conservation of mass, momentum and energy of a viscous fluid flow. Using Cartesian coordinates, these equations can be expressed in a conservative form as follows: oW þ r:G ¼ 0 ot

ð1Þ

The state vector W and the flux G = Ge  Gv decomposed in an inviscid and a viscous part are given by T

~ ; qE ; W ¼ ½q; qU ¼

~ ; qU ~U ~ þ p I; U ~ ðqE þ pÞT ; Ge ¼ ½qU ¼ ¼

ð2Þ

T

~ ~ Gv ¼ ½0; s ; s U q

~ the velocity, p the pressure and E where q is the density, U the total energy. For a Newtonian fluid, the shear stress ¼ tensor s is given by ¼

T

¼

~ þ ðrU ~ Þ Þ þ kr  U ~I s ¼ lðrU

ð3Þ

with l the dynamic viscosity and k the second coefficient of viscosity. The Stokes assumption reduces the Lame´’s relation to 2l + 3k = 0. The heat flux q is given by Fourier’s law ~ q ¼ K T rT

ð4Þ

with T the temperature and KT the thermal conductivity coefficient. The dynamic viscosity l is given by the Sutherland’s formulae  32 T T 0 þ Cs l ¼ l0 ð5Þ T 0 T þ Cs

2.3.1. Numerical methods The elsA software [10] solves the 3D compressible Navier–Stokes equations using a cell-centered finite-volume method. Integrating Eq. (1) over a domain X and applying Green’s divergence theorem yield the following integral form: Z I oU dV þ G:~ n dS ¼ 0 ð7Þ X ot oX with ~ n the outward normal of the boundary oX of the control volume X. The separated time/space discretization process leads to the following delta form: ADU nþ1 ¼ 

Dt RðU n Þ jXj

ð8Þ

where the residual R comes from the space discretization and depends on the conservative variable field Un. The Jacobian matrix A comes from the implicit time discretization and DUn+1 = Un+1  Un corresponds to the field correction also called the time increment. In this paper, a standard multigrid [13] method combined with a local time stepping is used in order to accelerate the convergence to steady solutions. The classical second order central scheme of Jameson et al. [12] is used for spatial discretization and a LU-SSOR implicit method [14] for solving the time integration system (8). 3. Response surface construction using the Kriging method 3.1. Introduction to the Kriging method The term ‘‘Kriging’’ denotes a family of interpolation methods where weighting coefficients are chosen to minimize

522

J.-C. Jouhaud et al. / Computers & Fluids 36 (2007) 520–529

the variance of the error [15]. First applied in geological analysis, it has been extended to many fields of application, including agriculture, human geography [16], epidemiology [18], biostatistics or archeology and computer experiments [23]. The Kriging method is a linear interpolation method that predicts values at unknown locations from data observed at known locations (control points). The estimation at any location is obtained as a weighted average of the known control point values. Since the estimator is linear, the error variance can be expressed using covariances. So, the covariogram [17] of the function to be estimated is used to define the weighting coefficients that determine the contribution of each control point to the prediction of new values at unsampled locations in the space spanned by the uncertain parameters. The covariogram represents the spatial correlation of the data. It is desirable that the estimator has the unbiasedness property that is, the expectation of the predicted value is equal to the expectation of the true value. This property is enforced thanks to requirements on the weighting coefficients. Numerous assumptions to determine the weighting coefficients give rise to many types of methods such as simple, ordinary and universal Kriging methods. The difference between these three models lies in the assumption about the mean value of the variable under study: simple Kriging requires a known constant mean value while ordinary Kriging assumes a constant but unknown mean value in the searching neighbourhood. Thus, these two approaches model a spatial surface as deviations from a constant mean, where the deviations are spatially correlated (application to steady problems). Within the universal Kriging framework, the local mean is approximated as a sum of low-order polynomials defined as functions of the spatial coordinates. This type of model is appropriate when there are some heterogeneities in the expectation of surface function (application to unsteady problems). In the following, only Ordinary Kriging (OK) will be considered. As a matter of fact, simple Kriging implies that functions to be estimated have a known mean, which is not true in the present application. And universal Kriging is not justified in the present work, since only steady flow simulations are concerned in the present work. The OK method enforces that the sum of all coefficients is one to give an unbiased estimator. Therefore, in deterministic simulations, the OK interpolator is said to be exact since the predicted values at observed input values are exactly equal to the observed (simulated) output values. Moreover, it assumes a stationary model so that statistical properties (such as mean, variance, covariance, etc.) do not depend on exact spatial locations, so the mean and variance of a variable at one location are equal to the mean and variance at another location. The expectation of the function is constant all over the domain. Moreover, the correlation between any two locations depends only on the vector separating them, and not on

their exact locations, leading to the definition of an isotropic model. 3.2. Mathematical formulation In the following, the parameter space or the study region is noted S  R2 . As explained in the previous section, the principle of the OK method is to estimate the value attributed to a surface function f at a location xp 2 S, where the true value is obviously unknown, using a linear function of data values FT = {f(x1), . . . , f(xi), . . . , f(xn)} at control points XT = {x1, . . ., xi, . . . , xn}: n X ci ðxp Þf ðxi Þ ¼ CT ðxp ÞF ðX Þ ð9Þ f^ ðxp Þ ¼ i¼1

where f^ is the linear estimator function of f, xp a vector in the two-dimensional parameter space S and {ci(xp)} the weighting functions. The region S contains n control points with data values {f(xi)}i=1, . . ., n. The main issue remains to determine the weighting functions with the isotropic stationary model assumption and the unbiasedness constraints. The assumption of an isotropic stationary model implies that the covariance C of the surface function f of two locations (xi, xj) is given by a function depending only on the distance between these two locations. The covariance C is defined by the following formulation C½f ðxi Þ; f ðxj Þ ¼ E½ðf ðxi Þ  Eðf ðxi ÞÞðf ðxj Þ  Eðf ðxj ÞÞ ¼ Cðjxi  xj jÞ Covariances are usually represented covariance matrix 0 r2 Cðkx1  x2 kÞ B Cðkx  x kÞ r2 2 1 B B C¼B .. .. @ . . Cðkxn  x1 kÞ

Cðkxn  x2 kÞ

ð10Þ as a matrix called the 1    Cðkx1  xn kÞ    Cðkx2  xn kÞ C C C .. .. C A . . 

r2 ¼ Cð0Þ ð11Þ

where r2 is the variance of the set of control points. Similarly, for plain explanation, the covariance vector ~ c is introduced 1 0 Cðkxp  x1 kÞ B Cðkxp  x2 kÞ C C B C ~ ð12Þ cðxp Þ ¼ B .. C B A @ . Cðkxp  xn kÞ The estimated value f^ ðxp Þ will most likely differs from the actual value f(xp) and the difference is called the estimation error ep ¼ f^ ðxp Þ  f ðxp Þ

ð13Þ

To obtain an optimal function estimator, the weighting functions ci(x), i = 1, . . . , n are computed so as to minimize

J.-C. Jouhaud et al. / Computers & Fluids 36 (2007) 520–529

the variance of the error or equivalently the mean square error (MSE) of f^ ðxp Þ o ½MSEðf^ ðxp ÞÞ ¼ 0 oCðxp Þ

ð14Þ

where MSE is defined by the following formula: MSEðf^ ðxp ÞÞ ¼ E½ðf^ ðxp Þ  f ðxp ÞÞ2 

i¼1 2

þr 2

 lnðcorÞ 9h 2 R=h ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 52 þ 52

ci ðxp ÞCðjxp  xi jÞ

i¼1

¼ CT ðxp ÞCCðxp Þ þ r2  2CT ðxp Þ~ cðxp Þ

ð15Þ

As a consequence, the solution of Eq. (14) is Cðxp Þ ¼ C 1~ cðxp Þ

ð16Þ

In the ordinary Kriging method, the stationary model assumption gives the following constraint: n X CT ðxp ÞI ¼ ci ðxp Þ ¼ 1 8 xp 2 S ð17Þ i¼1

To enforce the unbiasedness constraint about the summation of the weighting coefficients, a Lagrange multiplier k(xp) is introduced. The original problem (14) with the constraint (17) is then rewritten as o MSEf^ ðxp Þ þ CT ðxp ÞIkðxp Þ ¼ 0 oCðxp Þ

ð18Þ

leading to Cþ ðxp Þ ¼ C 1 cþ ðxp Þ þ~

ð19Þ

where  Cþ ðxp Þ ¼

Cðxp Þ kðxp Þ



 ;

Cþ ¼

C

1

1

0



 ;

~ cþ ðxÞ ¼

~ cðxÞ



1 ð20Þ

The estimator of the function f(x) is then given by cT ðxp ÞC 1 F ðX Þ f^ ðxp Þ ¼ ~

ð21Þ

In ordinary Kriging, the covariogram function is arbitrarily chosen from typical theoretical functions, or estimated from the observed data. Here, to estimate the surface function, exponential covariogram is used CðhÞ ¼ r2 expðhhÞ

ð22Þ

where h is the Euclidian distance between two locations of S. The parameter h governs the strength of the covariance between two locations. As the model is isotropic, it is not direction dependant. First, assuming that all control points are correlated, we have 8h 2 R 8ði; jÞ 2 n2 ;

ð24Þ

In our case,

j¼1 n X

The correlation is a decreasing function of the distance. Considering our grid distribution of control points, the minimum occurs on the diagonal of the grid (Fig. 4) for ði; jÞ ¼ ð~i; ~jÞ. Thus, by fixing the strength of the correlation cor between the two most distant control points, we obtain 9h 2 R= expðhkx~i  x~j kÞ ¼ cor

¼ E½f^ 2 ðxp Þ þ E½f 2 ðxp Þ  2E½f ðxp Þf^ ðxp Þ n X n X ¼ cj ðxp Þci ðxp ÞCðjxi  xj jÞ

523

9cor 20; 1½= expðhkxi  xj kÞ P cor ð23Þ

ð25Þ

Finally, according to our parameters and cost function the correlation factor chosen is cor = 0.2, giving h ’ 0.03. The model covariogram used is   3 2 h ð26Þ CðhÞ ¼ r exp  100 The linear system (19) is solved using the mathematical library of LAPACK (Linear Algebra PACKage) [19]. 3.3. Automatic Kriging computational suite In this work, the Kriging method has been implemented in an Automatic Kriging Computational Suite (AKCS) which is coupled with the elsA solver. Each step of the process consists in five stages: (1) Definition of points to be computed: (a) Sampling points in the uncertain parameter space are chosen. (b) Automatic mesh generation around the profile for each sampling point. (2) CFD computations with elsA: (a) Realization of the simulations for each sampling point in the uncertain parameter space. (3) Data processing: (a) Computations of the objective values at sampling points with the function to be interpolated. (b) Creation of data files for Kriging method. (4) Kriging method: (a) Reconstruction of the surface function. (b) Computation of the mean square error of Kriging method. (5) Automatic refinement: (a) Computation of the surface function gradients. (b) Refinement—Automatic determination of new points. The keystone of this Dynamic Kriging Computational Suite is provided by an Automatic Refinement Tool [22] which determines automatically new sampling points in regions where the surface function exhibits large gradients (see Figs. 4–6). This tool is based on a function called sensor which corresponds to an estimation of the surface

524

J.-C. Jouhaud et al. / Computers & Fluids 36 (2007) 520–529

function gradient using the a coefficient for each sampling point: a¼

! kgradðf Þk ! maxðkgradðf ÞkÞ

After the first computation on the initial mesh (coarse response surface), a is calculated for all cells. If its value is greater than a critical a0 value, the cells are flagged to be refined. Here, we have chosen a0 = 0.15 which was found to be the optimal choice to avoid the use of too numerous points in the vicinity of large gradient regions. Then, the smallest block that contains all flagged cells is determined. If the ratio of the number of flagged cells over the total number of cells in this block is greater than a defined value (usually around 80%), then it is created, otherwise the block is divided by two along the largest topological direction. The process is repeated for the two blocks until convergence. At last, coarse flagged cells are here linearly subdivided by a factor of 2. This form-recognizing process is called Grouping/Clustering Algorithm and it has been fully described by Quirk [21]. It is worth noting that the present implementation is non-intrusive, since it does not require any modification of the basic CFD tool: the coupling between the Kriging tool and the CFD solver is performed via data file transfers. Another important feature is the capability of using a multiresolution approach in the uncertain parameter space to minimize the error in the response surface interpolation. 4. Surrogate-model based multidisciplinary optimization of NACAXX16 airfoils 4.1. Airfoil shape modeling To illustrate the surrogate-model based optimization procedure, the problem of the multidisciplinary shape optimization of a 4-digit NACA airfoil is considered. Let us first recall the nomenclature of the four-digit NACA airfoil and discuss the choice of the optimization parameters. In a 4-digit NACA airfoil of chord c, the first digit m is the value of the maximum camber (expressed in percent of the chord), the second digit is the position p of the maximum camber (measured from the leading edge in tenths of the chord), and the last two digits denote the maximum thickness t of the airfoil in percent of the chord. Here, digits are replaced by real values and the thickness percentage is set equal to 16. The choice of the optimal shape can be performed using several approaches. A first one consists in defining control points on the shape surface and to search for their optimal displacement. This method usually leads to the definition of an optimization problem with a very high dimensionality. A second solution is to parameterize the shape using a very low dimensional model. It is chosen here to use the second way, which is best suited for the definition of a sur-

rogate model, and to use m and p as optimization parameters, leading to the definition of the two-dimensional optimization problem. The airfoil mean camber line is given by Eq. (27). m y c ¼ 2 ð2px  x2 Þ; x 2 ½0; p½; p ð27Þ m 2 ðð1  2pÞ þ 2px  x Þ; x 2 ½p; c yc ¼ ð1  pÞ2 The thickness distribution above and below the mean line is defined as pffiffiffi t ð0:2969 x  0:126x  0:3516x2 yt ¼  0:2 þ 0:2843x3  0:1015x4 Þ ð28Þ The coordinates for the airfoil upper surface (xu, yu) and lower surface (xl, yl) are given by xu ¼ x  y t sin h;

xl ¼ x þ y t sin h;

y u ¼ y c þ y t cos h;

tan h ¼

dhy c ; dhx

ð29Þ

y l ¼ y c  y t cos h

4.2. Multidisciplinary cost function definition It is chosen to use a multidisciplinary cost function in the present to show the efficiency of the method and its flexibility. The present design problem is associated with a multidisciplinary mono-objective cost function f ðxÞ : R2 ! R which is written as follows: ( Find xopt 2 S  R2 ; ð30Þ Such that jf ðxopt Þj P jf ðxÞj 8x 2 S with x = (m, p) the vector of design parameters and S is the subspace of permitted values. The objective is to reach a fixed lift coefficient C oL and to minimize both the drag coefficient and the intensity I of the far-field radiated noise measured at a control position. These coefficients are combined to maximize the cost function f ðxÞ ¼

1 ðC L 

C oL Þ2

þ C 2D þ I 2

ð31Þ

The quantities CL, CD and I vary in intervals whose bounds depend on the range of the optimization parameters. The different parts appearing in the cost function are made non-dimensional so that each of them ranges from 0 to 1. Computation of the mean drag and lift is straightforward using aerodynamic data retrieved from the RANS computation. An additional problem arises when the radiated noise is addressed, since acoustic sources are intrinsically tied to unsteady fluctuations. As a consequence, a model for the radiated noise must be used to evaluate the acoustic intensity. It is worth recalling that on a clean subsonic airfoil like the one considered in the present study the main acoustic source is the trailing edge noise, which originates in the scattering of acoustic waves generated

J.-C. Jouhaud et al. / Computers & Fluids 36 (2007) 520–529

x leading edge

7

θ

ψ

525

679.0 675.4 674.8 668.5 613.0 547.0 481.0 415.0 349.0 283.0 217.0 151.0 85.0 19.0

6

β

5

p

y H

4

z 3

receiver Fig. 1. Directivity angles used in the definition of the noise metric.

2

2

3

4

5

6

7

m

due to the passage of turbulent boundary layers over the trailing edge. Therefore, the trailing edge noise metric developed in [20] is retained in the present study. The far-field noise intensity per unit volume of acoustic sources is approximated as q1 5 1 Dðh; wÞ u l cos3 ðbÞ ð32Þ 2p3 a21 0 0 H2  where Dðh; /Þ ¼ 2 sin2 h2 sinðwÞ is the directivity term, and h and w the polar and azimuthal directivity angles, respectively. u0 and l0 are the characteristic velocity scale and length scale for turbulence. H is the distance to the receiver and b is the trailing edge sweep angle. q1 and a1 are respectively the free-stream density and sound speed. These different variables are displayed in Fig. 1. The characteristic turbulent velocity u0 is chosen to be a function of the turbulent kinetic energy k. pffiffiffi u0 ¼ maxð k Þ ð33Þ

Fig. 3. Isolines of the cost function at the first AKCS stage of the dynamic Kriging procedure. 7

I

6

p

5

4

3

2

The characteristic turbulent length scale l0 is modeled as u0 ð34Þ l0 ¼ x0

2

3

4

5

m

Fig. 4. Sample points in the parameter space at the first stage.

6

(7,7)

(2,2)

(7,2)

p

(2,7)

4

2

2

3

4

6

m

5

6

Fig. 2. NACA(2,2,16), NACA(2,7,16), NACA(7,2,16), NACA(7,7,16) profiles.

7

7

526

J.-C. Jouhaud et al. / Computers & Fluids 36 (2007) 520–529

with x0 the turbulence frequency observed at the location where k is maximum. The variables k and x are obtained from the turbulence model equations used in the fully turbulent RANS computations (without transition imposition). In the present study, a single control point defined by w = 90 and b = 0 is used to evaluate the cost function. 4.3. Flow parameters and numerical results The Mach number is set equal Ma = 0.2, the angle of incidence is taken equal to 2 and the Reynolds number based on the chord is Re = 106. The computational mesh 7

is made of two H-grids with 137 points in the streamwise direction and 65 in the normal direction of the airfoil. The subspace of authorized values for the optimization parameters is defined as S = [2, 7] · [2, 7]. The airfoil shapes associated to extreme values are displayed in Fig. 2. In the following, the figures will show the dimensionalized parameters. The CL and CD coefficients are respectively multiplied by 103 and 104 to emphasize the changes. To initialize the dynamic Kriging procedure, the optimization parameter subspace S is discretized with 6 sampling points uniformly distributed in each direction, leading to 36 evaluations of the cost function. Fig. 3 shows the isolines of the response surface for the cost function f reconstructed on a 500 · 500 grid. It is observed that the cost function exhibits two possible maxima. 7

6

700.1 647.8 595.5 543.2 490.8 438.5 386.2 333.9 281.6 229.2 176.9 124.6 72.3 20.0

6

p

5 5

p

4

4

3 3

2 2

3

4

m

5

6

7 2

Fig. 5. Sample points in the parameter space at the second stage—one level of refinement.

2

3

4

5

6

7

m Fig. 7. Isolines of the cost function at the second AKCS stage of the dynamic Kriging procedure.

7 7

6 670.0 619.9 569.8 519.8 469.7 419.6 369.5 319.4 269.4 219.3 169.2 119.1 69.1 19.0

6

5

p

p

5

4 4

3 3

2 2

3

4

m

5

6

7

Fig. 6. Sample points in the parameter space at the third stage—two levels of refinement.

2

2

3

4

m

5

6

7

Fig. 8. Isolines of the cost function at the third AKCS stage of the dynamic Kriging procedure.

J.-C. Jouhaud et al. / Computers & Fluids 36 (2007) 520–529

the last refinement level, which leads to the definition of the most accurate response surface, reveals the existence of a significant number of local maxima (Fig. 7). This is understood looking at the response surface of the lift coefficient (Fig. 9), the drag coefficient (Fig. 10) and the modeled noise intensity (Fig. 11), which are relatively smooth functions with different gradients. Their non-linear combination to obtain a multidisciplinary cost function

The dynamic refinement procedure yields the definition of two refinement levels. The whole optimization process required a total number of 257 CFD simulations, 36 coming from the initial step, 74 for the first level of refinement and 147 point for the last level. The final sampling set is displayed in Fig. 6, while the evolution of the response surface as a function of the grid refinement level is displayed in Figs. 4–6. It is observed that

.7

846 .9

780 .0

818

5

21 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1

991.0 962.2 933.4 904.6 875.7 846.9 818.1 780.0 760.4 731.6 702.8 674.0 645.2 616.3 587.5 558.7 529.9 501.0 472.2 443.4 414.6

21 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1

258.8 253.8 248.9 243.9 238.9 234.0 229.0 224.1 219.1 214.2 209.2 204.2 199.3 194.3 189.4 184.4 179.4 174.5 169.5 164.6 159.6

.1

4

87 5

p

.6

.1 .4

702.8

3

760

6 73 1.

616

2

8 18

8 702 .

501.0

2

90 4

.0

.6

529.9

645.2

587.5

674.0

616.3

558.7

472.2

3

780

731

501.0

443.4

414.6

4

5.7 87 .9 846 .1

81 8 .4

.0 2 645 .

529.9

5

760

674

5 587. 558.7

472.2 443.4

.3

.8

6 16

702

501.0

6

33.4

80.0

31.6

9 529.

7

527

6

7

m Fig. 9. Isolines of the lift coefficient CL (·103).

22 1 4.

20 9.2

.3

.3

204

199

19 4

.4

.4

18 4

18 9

.5

p

8 3.

9 8. 4

1 74

.2

179.4

164.6

199.3

194.3

189.4

184.4

174.5

169.5

159.6

179

2

23

.4

5 16 9.

159.6

3

23

179

6 16 4.

5

4

25

19

4.4

9.4

18

18

.5

6

2 4. 21 9.1 9.2 2 . 21 2 0 2 04 9.3 19 .3 4

74

.5

6 159.

169

7

2 2

3 3

4 4

5 5

6 6

m Fig. 10. Isolines of the drag coefficient CD (·104).

7 7

J.-C. Jouhaud et al. / Computers & Fluids 36 (2007) 520–529 13

9

9

15

13

8 7

190

7 9 10 11

11

10

7 8

10

12

13

14 11

8

9

6

5

11

5

4

685

10.9 17 18 0

16

14

170

2

6

20

12

12 9 5 11 6 8 3

4

7

12

528

p

13

8

13

12

4

10

11

9

7 12

3

21 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1

6.66E-02 5.09E-02 3.88E-02 2.96E-02 2.26E-02 1.32E-02 7.68E-03 2.60E-03 1.74E-03 1.55E-03 1.37E-03 1.22E-03 1.10E-03 1.03E-03 9.58E-04 8.99E-04 8.46E-04 6.75E-04 5.15E-04 3.93E-04 3.00E-04

6 8

10

5

2 2

3

4

5

6

7

m Fig. 11. Isolines of the acoustic intensity coefficient I.

makes several local maxima arising. This reconstruction of individual response surfaces for each element of the full multidisciplinary cost function also illustrates another advantage of the surrogate-based method: once these elementary response surfaces computed, new optimization problems based on cost functions defined as non-linear functions of the same base elements (drag, lift, radiated noise at control position) can be solved without additional use of the Navier–Stokes solver. Once the Kriging-based surrogate-model is built, the optimal parameters are obtained via a systematic exploration of the response surface. Since the present problem is associated with a two-dimensional problem, this method is efficient and the use of advanced methods is not necessary. We now address the problem of defining a robust shape optimization procedure based on the surrogate model. Here, the term robust means that the selected extremum retained as being the optimal solution (mopt, popt) must be such that a small variation of the design parameters, i.e. considering (m, p) 2 [mopt  e, mopt + e] · [popt  e, popt + e] (where e is an error tolerance) must not result in a dramatic change in the cost function value. A simple method to achieve that is to use the filtered response surface f , which is defined as Z  f ðxÞ ¼ Gðe; jx  yjÞf ðyÞ dy ð35Þ S

where G(e, jx  yj) is the regularization kernel with a characteristic cutoff scale e. The use of this smoothing function removes the narrow maxima, which are not good candidates for robust optimization.

0.4

0.2

0

-0.2

-0.4 0

0.2

0.4

0.6

0.8

1

Fig. 12. NACA(mrobust, probust, 16) profile (solid line) and symmetric NACA(0,0,16) profile (dashed line).

Using a top-hat filter and e = 0.05, one obtains the following robust optimal values mrobust = 5.8 and probust = 4.075. The filtered response surface is displayed in Fig. 8, while the associated wing profile is shown in Fig. 12. 5. Conclusions A surrogate-model based shape optimization method was proposed, which relies on the Kriging interpolation technique. The Kriging method is used to build the response surface of the cost function to be maximized. In order to improve the efficiency of the method, the number of required solution of the Navier–Stokes equations is decreased thanks to an Adaptive Mesh Refinement type

J.-C. Jouhaud et al. / Computers & Fluids 36 (2007) 520–529

technique in the optimization parameter space. The case of the multidisciplinary optimization of a 2D airfoil was used to assess the proposed method. The method used in the present paper requires the use of a significant number of sampling points, and must be modified to handle 3D configurations with a large number of optimization parameters. The issues associated with large-dimension optimization spaces and sampling set optimization are beyond the scope of the present paper. Screening methods to decrease the number of effective parameters and methods for generation of optimal sampling sets were therefore not considered, and will be investigated in future works. References [1] Marsden AL, Wang M, Dennis JE, Moin P. Optimal aeroacoustic shape design using the surrogate management framework. Optim Eng 2004;5(2):235–62. Special Issue on Surrogate Optimization. [2] Marsden AL, Wang M, Koumoutsakos P. Optimal aeroacoustic shape design using approximation modeling. Annual Research Briefs, Center for Turbulence Research, 2002, p. 201–13. [3] El-Beltagy MA, Keane AJ. Evolutionary optimization for computationally expensive problems using Gaussian processes. In: Arabnia H, editor. Proc. int. conf. artificial intelligence IC-Al’2001. Las Vegas: CSREA Press; 2001. p. 708–14. [4] Booker AJ, Dennis Jr JE, Frank PD, Serafini DB, Torczon V, Trosset MW. A rigourous framework for optimization of expensive functions by surrogates. Struct Optim 1999;17(1):1–13. [5] Buche D, Schraudolph NN, Koumoutsakos P. Accelerating evolutionary algorithms with Gaussian process fitness function models. IEEE Trans Systems, Man Cybern. Part C 2005;35(2):183–94. [6] Simpson TW, Mauery TM, Korte JJ, Mistree F. Comparaison of response surface and Kriging models for multidisciplinary design optimization. In: 7th AIAA/USAF/NASA/ISSMO symposium on multidisciplinary analysis and optimization, St. Louis, USA, AIAA Paper 98–4758, 1998. [7] Torczon V, Trosset MW. Using approximations to accelerate engineering design optimization, ICASE Report No. 98–33, Tech. Rep., NASA Langley Research Center Hampton, VA 23681–2199, 1998.

529

[8] Jin Y, Olhofer M, Sendhoff B. A framework for evolutionary optimization with approximate fitness functions. IEEE Trans Evol Comput 2002;6(5):481–94. [9] Jin Y. A comprehensive survey of fitness approximation in evolutionary computation. Soft Comput 2005;9(1):3–12. [10] Cambier L, Gazaix M. elsA: an efficient object-oriented solution to CFD complexity. In: 40th AIAA aerospace science meeting and exhibit, Reno, AIAA Paper 2002-0108, 2002. [11] Wilcox DC. Simulation of transition with a two-equation turbulence model. AIAA J 1994;32(2):247–55. [12] Jameson A, Schmidt W, Turkel E. Numerical solutions of the Euler equations by finite volume methods using Runge–Kutta timestepping schemes, AIAA Paper 81-1259, 1981. [13] Jameson A. Steady state solutions of the Euler equations for transonic flow by a multigrid method. Advances in scientific comp. Academic Press; 1982. p. 37–70. [14] Yoon S, Jameson A. An LU-SSOR scheme for the Euler and Navier– Stokes equations, AIAA Paper 87-0600, 1987. [15] Krige DG. A statistical approach to some basic mine valuations problems on the witwatersrand. J Chem, Metal Mining Soc South Africa 1951;52:119–39. [16] Oliver MA, Webster R. Kriging: a method of interpolation for geographical information system. Int J Geograph Informat Systems 1990;4(3):313–32. [17] Journel AG, Huijbregts CJ. Mining geostatistics. Academic Press; 1978. [18] Christakos G, Hristopulos DT. Spatiotemporal environmental health modeling. Kluwer; 1988. [19] Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J, et al. LAPACK users’ guide. 3rd ed. Society for Industrial and Applied Mathematics; 1999. [20] Hosder S, Schetz JA, Grossman B, Mason WH. Airframe noise modeling appropriate for multidisciplinary design and optimization. In: 42nd AIAA aerospace science meeting and exhibit, Reno, Nevada, January 5–8, 2004. [21] Quirk JJ. An adaptive grid algorithm for computational shock hydrodynamics. Ph.D. Thesis, Cranfield Institute of Technology, College of Aeronautics, 1991. [22] Jouhaud J-C, Montagnac M, Tourrette L. A multigrid adaptive mesh refinement strategy for 3D aerodynamics design. Int J Numer Meth Fluids 2005;47:367–85. [23] Sacks J, Welch WJ, Mitchell TJ, Wynn HP. Design and analysis of computer experiments. Statist Sci 1989;4:409–35.