An efficient geometric constraint handling method for surrogate-based aerodynamic shape optimization

Abstract

Handling a large number of geometric constraints brings a big challenge to the surrogate-based aerodynamic shape optimization (ASO) driven by computational fluid dynamics (CFD). It is not feasible to calculate the geometric constraint functions directly during the sub-optimization of a surrogate-based optimization, as the geometric constraint functions are to be evaluated thousands of times for each updating cycle and the total cost of a number of cycles can be prohibitive. This article proposes an efficient method of handling geometric constraints within the framework of a surrogate-based optimization to address this problem. The core idea is to use the Kreisselmeier-Steinhauser (KS) method to aggregate all geometric constraints into one that can be approximated by a cheap surrogate model, in order to avoid the large computational cost associated with tremendous calculation of geometric constraint values. The proposed method is verified by an analytical test case. Then, the proposed method is demonstrated and compared with the methods of building surrogate models of all geometric constraints and calculating all geometric constraints directly during each sub-optimization by drag minimizations of NACA0012 air foil and ONERA M6 wing in transonic flows. To investigate the ability of the proposed method for handling various geometric constraints, drag minimization of CRM wing in viscous transonic flow driven by CFD is performed. Results show that the proposed method can dramatically improve the optimization efficiency of ASO with the number of geometric constraints ranging from 15 to 1429 and the number of types of geometric constraints up to 3, which offers great potential for handling a larger number and more types of geometric constraints.

KEYWORDS:

• Surrogate-based optimization
• aerodynamic shape optimization
• constraint handling method
• constrained optimization

Abbreviations

SBO	=	surrogate-based optimization
ASO	=	aerodynamic shape optimization
CFD	=	computational fluid dynamics
DoE	=	design of experiment
Ks	=	Kreisselmeier-Steinhauser
LHS	=	Latin hypercube sampling
ISC	=	infill-sampling criteria
EI	=	expected improvement
MSP	=	minimizing surrogate prediction
LCB	=	lower confidence bounding
MSE	=	mean squared error
PI	=	probability of improvement

1. Introduction

Surrogate-based optimization (SBO) has become one of the most attractive methods for aerodynamic shape optimization (ASO) driven by computationally expensive computational fluid dynamics (CFD) simulations. In a SBO, the expensive cost functions are approximated by cheap-to-evaluate surrogate models to drive the adding and evaluation of new sample point(s) towards the global optimum (Queipo et al., 2005; Forrester & Keane, 2009; Haftka et al., 2016). A typical process of a SBO includes two stages (Parr et al., 2012; Han, 2016a; Liu et al., 2019): (1) using the sample points generated by a design of experiment (DoE) method, initial surrogate models are constructed for the objective and constraint functions, respectively; (2) new samples are suggested by sub-optimization defined by certain infill-sampling criteria, then their responses obtained by a solver (such as CFD solver) are used to update the surrogate models. It employs much less function evaluations to search the global optimum compared with a genetic algorithm, especially when expensive high-fidelity CFD or multidisciplinary analysis are employed (Viana et al., 2014; Han et al., 2013; Lu et al., 2021; Park et al., 2020). SBO is free of gradient information and could be used to solve the problems with high nonlinear and multi-modal design space (Haftka et al., 2016; Han, 2016b). Therefore, SBO gained much attention in the area of engineering design optimization. Despite the popularity of SBO, it is still suffering from the difficulty associated with ‘curse of dimensionality', which means that the computational cost grows exponentially as the number of design variables increases. To improve the efficiency of a SBO, some advanced technologies have been proposed, such as evaluating the responses of the samples in parallel (Liu et al., 2017), variable-fidelity surrogate modeling (Courrier et al., 2016; Han et al., 2013; Han et al., 2020; Ha et al., 2014; Jo et al., 2016; Leifsson & Koziel, 2015, 2016; Zhang et al., 2018), surrogate modeling using cheap gradients (Laurenceau et al., 2012; Han et al., 2017; Mortished et al., 2016; Bouhlel & Martins, 2019), and applying dimension reduction in the optimization (Li et al., 2019; Bouhlel et al., 2018).

However, handling a large number of geometric constraints presents a big challenge to the surrogate-based aerodynamic shape optimization. First, it is not feasible to calculate the geometric constraint functions directly during the sub-optimization performed by a global optimization algorithm such as genetic algorithm within the framework of a SBO, in which these functions defined by massive surface meshes and parameterization methods such as free form deformation (FFD) are evaluated thousands of times for each updating cycle, the total cost of a number of cycles can be prohibitive. Second, the method of building surrogate models for all geometric constraints during the optimization will suffer from difficulties associated with the prohibitive computational cost caused by training the large number of hyper-parameters in a one by one manner. Finally, it is improper to replace all constraints with the most violated one such as the maximum thickness constraint, since the maximum thickness t_max may changes during the optimization, which leads the constraint function to be discontinuous, and may result in an increased number of iterations of the optimization and even premature convergence. At present, a large number of geometric constraints have been widely applied in gradient-based aerodynamic shape optimization of different aerodynamic configurations, such as the one-stage turbine (Xu et al., 2015), regional Jet (Bons et al., 2018), large civil aircraft wing (Lei et al., 2019; Kenway & Martins, 2016; Chen et al., 2016), strut-braced wing (Secco & Martins, 2019), 10 MW wind turbine (Madsen et al., 2019), blended-wing-body aircraft (Reist & Zingg, 2017; Lyu et al., 2017), as well as fairing (Brelje et al., 2020). Therefore, it is necessary to investigate an efficient geometric constraint handling method for SBO to solve ASO problems involving a large number of geometric constraint.

In the field of structural optimization, the KS method proposed by Kreisselmeier and Steinhauser (1979) provided a way to handle large-scale constraints, since it can make good use of the KS function to aggregate a large number of constraints into one. The KS method is developed by Raspanti et al. (2000) to handle inequality constraints in a nonlinear programming problem. Akguen et al. (2001) used it to improve the efficiency of solving the sensitivities of large-scale stress and buckling constraints in gradient-based structural optimizations. Then, the KS method was widely used in the gradient-based structural optimization (Paris et al., 2010; Lee & Martins, 2012) and aerostructural optimization (Martins, 2002; Martins et al., 2002; Hwang et al., 2014). It is found that the KS function with a small aggregation parameter ρ is overly conservative. It becomes more accurate as ρ increases, but may cause numerical difficulties and oscillations in a gradient-based optimization due to the ill-conditioned Hessian matrix. Poon and Martins (2005, 2007) addressed this problem by an adaptive KS method, more accurate and robust gradient-based optimization results were obtained compared with the KS method. The adaptive KS method was further used for solving the gradient-based aerostructural optimization (Kenway et al., 2014; Kenway & Martins, 2014) and aeroservoelastic optimization (Haghighat et al., 2012). However, it is mainly used to solve the numerical difficulties and oscillations in gradient-based optimizations, but cannot provide a higher approximate accuracy to the original constraints than the KS function with a larger value of ρ, theoretically. Different KS functions were introduced into the surrogate-based aerostructural optimization by Zhang et al. (2019). Results show that the KS method can handle large-scale structural stress constraints by aggregating all constraints to a continuous one with accurate approximation of the original constraint bound. However, the KS method has not been applied to handle strict geometric constraints of an aerodynamic shape optimization within the framework of a SBO, which motivates the study of this work.

The objective of this article is to investigate an efficient method to solve CFD-driven aerodynamic shape optimization problems involving a large number of geometric constraints within the framework of a SBO, in order to effectively improve the capability of SBO in engineering design. The core idea is to use the KS method to aggregate all geometric constraints into one that can be approximated by a cheap surrogate model, and to calculate the all geometric constraint values during the sub-optimization by using the surrogate model of the aggregated constraint, which leads SBO to avoid the large computational cost associated with building surrogate models for all geometric constraints and calculating all geometric constraint values tremendously in each iteration of optimization. This is nearly a blank area according to literature survey, which inspires the study of this article.

This article continues in section 2 to give a detailed description of flow simulation method used in this article. In section 3, surrogate-based optimization for constrained aerodynamic shape optimization is described. In section 4, the proposed method for handling geometric constraints in SBO is dedicated. In section 5, first, the G9 analytical test case is used to verify the proposed method. Second, the drag minimizations of the NACA0012 airfoil and ONERA M6 wing with a small number of geometric constraints are used to compare different geometric constraint handling methods. Moreover, the ONERA M6 wing optimization example with 1000 geometric constraints is used to further study the ability of handling a large number of geometric constraints by the proposed method. Finally, to investigate the ability of the proposed method for handling various geometric constraints, drag minimization of NACA CRM wing in viscous transonic flow driven by CFD is performed.

2. Flow simulation method

In this article, an in-house code ‘PMNS3D' is used as the CFD solver to perform 2D and 3D CFD simulations. It is a finite-volume, cell-centered multi-block solver for solving the Euler and Reynolds-averaged Navier–Stokes (RANS) equations. The Jameson-Schmidt-Turkel (JST) scheme is used for the spatial discretization. The Lower–Upper Symmetric Gauss Seidel (LU-SGS) scheme is implemented for time advancing with multi-grid technique to accelerate the iteration converging to steady states. The turbulence viscous coefficient is calculated by Spalart-Allmaras (SA) one-equation turbulence model. The in-house CFD code has been validated and used to the optimization of various airfoils (Han et al., 2019; Zhang et al., 2018) and wings (Liu et al., 2017; Han et al., 2020).

2.1. Governing equation of flow and discretization scheme

The Navier–Stokes equations in integral form for the control volume V and surface element dS can be defined as (1) $\underset{\mathrm{\Omega }}{\int \int \int }\frac{\mathrm{\partial }Q}{\mathrm{\partial }t}dV+\int \underset{\mathrm{\partial }\mathrm{\Omega }}{\int }\mathit{F}\bullet \mathit{n}dS=\int \underset{\mathrm{\partial }\mathrm{\Omega }}{\int }{\mathit{F}}_{\mathit{v}}\bullet \mathit{n}dS\text{,}$(1) where $\mathbf{\text{Q}}\text{ = }[\rho ,\rho u,\rho v,\rho w,\rho E{]}^T$ denotes the conservation variables, ${\rho }_a$ is the density. $\mathbf{\text{q}}=(u,v,w{)}^T$is the Cartesian velocity vector. $E$ is the specific total energy. n is the surface normal of control volume boundary $\mathrm{\partial }\Omega$. $\mathit{F},{\mathit{F}}_{\mathit{v}}$ denote the inviscid and viscous flux vector, respectively. Their detailed form can be found in reference (Han, 2021). For the Euler equations, the viscous flux vector of Eq. (1)(1) $\underset{\mathrm{\Omega }}{\int \int \int }\frac{\mathrm{\partial }Q}{\mathrm{\partial }t}dV+\int \underset{\mathrm{\partial }\mathrm{\Omega }}{\int }\mathit{F}\bullet \mathit{n}dS=\int \underset{\mathrm{\partial }\mathrm{\Omega }}{\int }{\mathit{F}}_{\mathit{v}}\bullet \mathit{n}dS\text{,}$(1) is neglected.

By using the cell-centered finite volume method for each grid cell i, the spatial discretization of Eq. (1)(1) $\underset{\mathrm{\Omega }}{\int \int \int }\frac{\mathrm{\partial }Q}{\mathrm{\partial }t}dV+\int \underset{\mathrm{\partial }\mathrm{\Omega }}{\int }\mathit{F}\bullet \mathit{n}dS=\int \underset{\mathrm{\partial }\mathrm{\Omega }}{\int }{\mathit{F}}_{\mathit{v}}\bullet \mathit{n}dS\text{,}$(1) can be expressed as (2) $\frac{\mathrm{\partial }({\mathit{Q}}_i{V}_i)}{\mathrm{\partial }t}+{\overline{\mathit{R}}}_i=\mathbf{\text{0}}\text{,}$(2) where $V_i$ is the volume of current grid cell i. ${\overline{\mathit{R}}}_i={{\mathit{R}}_i}^c-{{\mathit{R}}_i}^v-{\mathit{D}}_i$ is the total residual of Eq. (2)(2) $\frac{\mathrm{\partial }({\mathit{Q}}_i{V}_i)}{\mathrm{\partial }t}+{\overline{\mathit{R}}}_i=\mathbf{\text{0}}\text{,}$(2) , ${\mathit{R}}_i,{{\mathit{R}}_i}^v$are the residual values of the inviscid flux and the viscous flux respectively, and ${\mathit{D}}_i$ is the artificial dissipation term added to prevent non-physical oscillations.

Implicit solution algorithms have been implemented based on the LU-SGS method. The first-order backward difference is performed on the time derivative term, and the viscous flux is treated explicitly, Eq. (2)(2) $\frac{\mathrm{\partial }({\mathit{Q}}_i{V}_i)}{\mathrm{\partial }t}+{\overline{\mathit{R}}}_i=\mathbf{\text{0}}\text{,}$(2) can be written into the following implicit form (3) $\frac{V_i({{\mathit{Q}}_i}^{k+1}-{{\mathit{Q}}_i}^k)}{\Delta t}+{\overline{\mathit{R}}}_i^{{}^{k+1}}=\mathbf{\text{0}}\text{,}$(3) where ${{\overline{\mathit{R}}}_i}^{{}^{k+1}}$can be written as (4) ${{\overline{\mathit{R}}}_i}^{{}^{k+1}}={\mathit{R}}_i^{c,k+1}-{\mathit{R}}_i^{v,k}-{\mathit{D}}_i^k\text{,}$(4) where the superscript k represents the number of the time step. Furthermore, the JST scheme is used for the spatial discretization and the SA turbulence model is used to compute turbulence viscosity coefficient. More details can be found in reference (Han, 2021).

2.2. Boundary condition

2.2.1. Wall boundary condition

The no-slip boundary condition and the adiabatic wall boundary condition are used as the boundary conditions for RANS-based CFD simulations. Since the wall is stationary, $\mathbf{\text{q}}=(u,v,w{)}^T\text{ = }\mathbf{\text{0}}$ at the wall. For Euler based CFD simulations, as the viscosity is neglected, the normal component of the relative velocity is ${\mathbf{\text{q}}}_b\text{ = }(u_b,v_b,w_b{)}^T\text{ = }\mathbf{\text{0}}$. The no-slip boundary condition can be defined as (5) $(\mathbf{\text{q}}-{\mathbf{\text{q}}}_b)\cdot \mathit{n}\text{ = }\mathbf{\text{0}}.$(5) The adiabatic wall boundary condition assumes that there is no heat exchange between the wall and the fluid, therefore, the gradient of temperature (6) $\frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}=0$(6)

2.2.2. Symmetric boundary condition

The symmetry boundary condition defines a mirror face that reflects all the flow distribution to the other side. The are no flows and scalar flux across the boundary. By using this boundary condition, the domain can essentially be halved to reduce the time to achieve a solution.

2.2.3. Far-field boundary condition

In this paper, the local one-dimensional Riemannian invariants are used to deal with the far-field boundary. the local one-dimensional Riemannian invariants are defined as (7) $\begin{array}{rl}{R}^{-}& =q_n-\frac{2a}{\gamma -1}\\ R^{+}& =q_n+\frac{2a}{\gamma -1},\end{array}$(7) where $q_n$ denotes the normal component of the velocity at the far-field boundary, a is the sound velocity.

3. Surrogate-Based optimization for constrained aerodynamic shape optimization driven by CFD

3.1. Framework of surrogate-based constrained optimization

The purpose of this article is to explore an efficient method that can be used to handle a large number of geometric constraints of an aerodynamic shape optimization problem in the framework of a SBO. The generic constrained single-objective optimization problem can be defined as (8) $\begin{array}{l}\text{min}\text{. }f(\mathit{x})\\ \text{s}\text{.t}\text{. }{c}_j(\mathit{x})\le 0,j=1,\cdot \cdot \cdot ,N_c,\\ {\mathit{x}}_{\text{l}ow}\le \mathit{x}\le {\mathit{x}}_{up}\end{array}$(8) where f (x) is the objective function, x is the vector of design variables, $c_j(\mathit{x}),j=1,\cdot \cdot \cdot ,N_c$ are the constraint functions, $N_c$ denotes the number of the constraints. ${\mathit{x}}_{\text{l}ow}$ and ${\mathit{x}}_{up}$ denote the lower and upper boundary of the design space, respectively.

An in-house surrogate-based optimizer, SurroOpt, is applied to address the optimization problem mentioned above. It has been tested with numerous benchmark numerical examples and practical engineering problems (Han, 2016b; Liu et al., 2017; Han et al., 2020). Figure 1 sketches the framework of SurroOpt. The optimization process could be described as

1. Latin hypercube sampling (LHS) is selected to generate initial sample points. Then, the responses of the objective, aerodynamic constraints (e.g. C_L, C_M) are evaluated in parallel by the in-house CFD code, and geometric constraints are evaluated by another in-house codes.

2. Based on the initial sample points obtained by LHS and corresponding responses, initial surrogate models of the objective and constraint functions mentioned above are built and tuned. In this article, kriging surrogate model is used.

3. The ISC in conjunction with expected improvement (EI), mean squared error (MSE), minimizing surrogate prediction (MSP), probability of improvement (PI) and lower confidence bounding (LCB) are used to obtain five new sample points. The new samples are selected sequentially via sub-optimization defined by certain ISC via the genetic algorithm. Their responses are evaluated by the same in-house CFD code in parallel later again.

4. The new sample points and corresponding responses are added to the sampled database and the kriging surrogate models are updated.

5. Steps 3 and 4 are repeated until one of the termination conditions is satisfied.

Figure 1. Framework of surrogate-based optimizer, SurroOpt (Han, 2016b).

3.2. Infill-sampling criteria (ISC) and constraint handling methods

As we described in previous section, new feasible samples are selected by ISC. This step is called the sub-optimization in SBO, which turns out to be a constrained problem that involves the maximization or minimization of a function defined by the ISC. In SurroOpt, ISC involving EI, MSP, PI, LCB and MSE are combined to generate five new sample points. Note that each infill criterion can only be used to select one point. Since each infill criterion has its advantages and disadvantages, the simultaneous implementation of all criteria can complement each other. According to these ISC, the constraints handling methods can be divided into two categories.

3.2.1. MSP and LCB infill criteria with constraint handling

MSP considers adding a new sample at the current optimum. It is not guaranteed that MSP can converge to the global optimum due to the lack of exploration. The sub-optimization problem of MSP is (9) $\begin{array}{l}\text{min }\hat{f}(\mathit{x})\\ \text{s}\text{.t}\text{. }{\hat{c}}_i(\mathit{x})\le 0,i=1,\cdots ,N_C\end{array}$(9) where ‘${}^{\wedge }$' means the surrogate model of the corresponding item of Eq. (8)(8) $\begin{array}{l}\text{min}\text{. }f(\mathit{x})\\ \text{s}\text{.t}\text{. }{c}_j(\mathit{x})\le 0,j=1,\cdot \cdot \cdot ,N_c,\\ {\mathit{x}}_{\text{l}ow}\le \mathit{x}\le {\mathit{x}}_{up}\end{array}$(8) .

LCB is developed by taking the uncertainty of surrogate models into account, the following sub-optimization problem would be solved (10) $\begin{array}{l}\text{min LCB(}\mathit{x}\text{)}=\hat{f}\text{(}\mathit{x}\text{)}-A\hat{s}\text{(}\mathit{x}\text{)}\\ \text{s}\text{.t}\text{. }{\hat{c}}_i\text{(}\mathit{x}\text{)}\le \text{0, }i=1,\dots ,N_C,\end{array}$(10) where A denotes a user-defined parameter. $\hat{s}\text{(}\mathit{x}\text{)}$is the estimated standard deviation of the surrogate model prediction.

For these two ISC, in order to address a constrained problem, the global optimization algorithm GA is used. The constrained problem is transformed into an unconstrained one by adding a one pass penalty function to the model prediction (Deb, 2000), (11) $J(\mathbf{\text{x}})=\{\begin{array}{l}\hat{j}(\mathit{x})\\ \text{if}{\hat{c}}_i(\mathit{x})\le 0\mathrm{\forall }i=1,2,\cdots ,N_{\text{C}}\\ {\hat{j}}_{max}+\sum _{i=1}^{N_{\text{C}}}max({\hat{c}}_i(\mathit{x}),0)\text{otherwise}\end{array},$(11) where $J(\mathbf{\text{x}})$ is the objective function value of the worst feasible solution in the population, $\hat{j}$ is the objective function of ISC for Eq. (9)(9) $\begin{array}{l}\text{min }\hat{f}(\mathit{x})\\ \text{s}\text{.t}\text{. }{\hat{c}}_i(\mathit{x})\le 0,i=1,\cdots ,N_C\end{array}$(9) and Eq. (10)(10) $\begin{array}{l}\text{min LCB(}\mathit{x}\text{)}=\hat{f}\text{(}\mathit{x}\text{)}-A\hat{s}\text{(}\mathit{x}\text{)}\\ \text{s}\text{.t}\text{. }{\hat{c}}_i\text{(}\mathit{x}\text{)}\le \text{0, }i=1,\dots ,N_C,\end{array}$(10) , $max(\cdot )$ returns the largest value. ${\hat{c}}_i(\mathbf{\text{x}})$ is the normalized constraint model.

3.2.2. PI, EI and MSE infill criteria with constraint handling

Using the statistical feature of kriging model, both the exploitation and exploration of surrogate model can be accounted to select new sample points by estimating the regions with a high probability of improvement. The probability of improvement function is defined as (12) $P[I(\mathbf{\text{x}})]\text{ = }\mathrm{\Phi }\left(\frac{f_{min}\text{ - }\hat{f}(\mathbf{\text{x}})}{s(\mathbf{\text{x}})}\right).$(12) Correspondingly, the sub-optimization of PI is given by (13) $\begin{array}{l}\text{max }P(I(\mathbf{\text{x}}))\\ \text{s}\text{.t}\text{. }{\hat{c}}_i(\mathbf{\text{x}})\le 0,i=1,\cdots ,N_C\end{array}$(13) PI does not indicate the magnitude of improvement. It only identifies areas where some improvement may be expected. The magnitude of improvement is exposed in the concept of EI, (14) $\begin{array}{rl}& E[I(\mathbf{\text{x}})]\\ & =\{\begin{array}{cl}(f_{min}-\hat{f})\mathrm{\Phi }\left(\frac{f_{min}-\hat{f}}{s}\right)+s\varphi (\frac{f_{min}-\hat{f}}{s})& \mathrm{if}\mathrm{s}>0\\ 0& \text{ if s = 0}\end{array}.\end{array}$(14) Thus, the sub-optimization to be solved is (15) $\begin{array}{l}\text{max }E(I(\mathit{x}))\\ \text{s}\text{.t}\text{. }{\hat{c}}_i(\mathit{x})\le 0,i=1,\cdots ,N_C.\end{array}$(15) In addition, the new update samples can be added at the location with maximum MSE, to improve the global accuracy of surrogate model, (16) $\begin{array}{l}\text{max }MSE(\hat{f}(\mathit{x}))\\ \text{s}\text{.t}\text{. }{\hat{c}}_i(\mathit{x})\le 0,i=1,\cdots ,N_C.\end{array}$(16) The constrained optimization method using penalty function mentioned in section 3.2.1 can also be used for PI, EI and MSE method. Making use of the statistical feature of kriging model, a probabilistic approach is proposed (Schonlau, 1997). Taking the EI criterion as an example, the probabilistic approach is implemented by using a product of the expected improvement of the objective function and the PoF calculated from the constraint functions. The PoF is calculated in the same manner as the probability of improvement, which is given by (17) $P[C_i(\mathbf{\text{x}})\le 0]=\mathrm{\Phi }\left(-\frac{{\hat{c}}_i(\mathit{x})}{s_{c,i}(\mathit{x})}\right),$(17) where ${\hat{c}}_i(\mathit{x})$ is the surrogate model of the constraint function, $s_{c,i}$ is the root mean square error (RMSE) of ${\hat{c}}_i(\mathit{x})$. Then a new sample can be obtained by solving the following optimization problem: (18) $\text{max }J(\mathit{x})\cdot \prod _{i=1}^{N_C}P[C_i(\mathit{x})\le 0],$(18) where $J(\mathit{x})$ is $P(I(\mathit{x}))$, $E(I(\mathit{x}))$, $MSE(\hat{f}(\mathit{x}))$ respectively for PI, EI and MSE infill criteria. Parr et al. (2010) found that the probabilistic method of constraint handling performs better than the penalty method on a number of problems.

4. Methods of handling geometric constraints

In ASO, a large number of geometric constraints are needed due to the structure or assembly requirements, e.g. thickness constraints on the wing or cross-sectional area constraints on fuselage/nacelle, which becomes a great challenge for SBO. Generally, the geometric parameter (e.g. thickness or cross-sectional area) should be no less than a certain value at required locations, (19) $\begin{array}{ll}\text{min}\text{. }& f(\mathit{x})\\ \text{w}\text{.r}\text{.t}\text{.}& {\mathit{x}}_{low}\le \mathit{x}\le {\mathit{x}}_{up}\\ \text{s}\text{.t}\text{. }& c_j(\mathit{x})\le 0,j=1,\cdot \cdot \cdot ,N_c\\ & g_i(\mathit{x})=g{p}_{i,0}(\mathit{x})-g{p}_i(\mathit{x})\le 0,i=1,\cdots ,N_t\end{array}$(19) where $c_i(\mathit{x})$ denotes the remaining constraints, such as the lift and moment constraints, and $N_c$ is the number of $c_i(\mathit{x})$. $g_i(\mathit{x})$ denotes the geometric constraint. $g{p}_{0,i}$,$g{p}_i$ stand for the geometric parameter of the baseline aerodynamic shape and the new shape, respectively. $N_t$ is the number of geometric constraints.

4.1. Evaluating all geometric constraints directly (EGCD)

Solving geometric constraints directly during the sub-optimization of a SBO can provide true values of the geometric constraints for the optimization. Then, Eq. (19)(19) $\begin{array}{ll}\text{min}\text{. }& f(\mathit{x})\\ \text{w}\text{.r}\text{.t}\text{.}& {\mathit{x}}_{low}\le \mathit{x}\le {\mathit{x}}_{up}\\ \text{s}\text{.t}\text{. }& c_j(\mathit{x})\le 0,j=1,\cdot \cdot \cdot ,N_c\\ & g_i(\mathit{x})=g{p}_{i,0}(\mathit{x})-g{p}_i(\mathit{x})\le 0,i=1,\cdots ,N_t\end{array}$(19) becomes (20) $\begin{array}{ll}min.& \hat{f}(\mathit{x})\\ \text{w}\text{.r}\text{.t}\text{.}& {\mathit{x}}_{low}\le \mathit{x}\le {\mathit{x}}_{up}\\ \text{s}\text{.t}\text{.}& {\hat{c}}_j(\mathit{x})\le 0,j=1,\cdot \cdot \cdot ,N_c\\ & g_i(\mathit{x})\le 0,i=1,\cdots ,N_t,\end{array}$(20) where $g_i(\mathit{x})$ is the real geometric constraint value, which is calculated directly in the sub-optimization. The values of $f(\mathit{x})$ and $c(\mathit{x})$ are still obtained by predicting their surrogate models.

4.2. Building surrogate models for all geometric constraints (BSMGC)

This method is to build surrogate models for the objective and constraints in a one-by-one manner, Eq. (19)(19) $\begin{array}{ll}\text{min}\text{. }& f(\mathit{x})\\ \text{w}\text{.r}\text{.t}\text{.}& {\mathit{x}}_{low}\le \mathit{x}\le {\mathit{x}}_{up}\\ \text{s}\text{.t}\text{. }& c_j(\mathit{x})\le 0,j=1,\cdot \cdot \cdot ,N_c\\ & g_i(\mathit{x})=g{p}_{i,0}(\mathit{x})-g{p}_i(\mathit{x})\le 0,i=1,\cdots ,N_t\end{array}$(19) becomes (21) $\begin{array}{ll}min.& \hat{f}(\mathit{x})\\ \text{w}\text{.r}\text{.t}\text{.}& {\mathit{x}}_{low}\le \mathit{x}\le {\mathit{x}}_{up}\\ \text{s}\text{.t}\text{.}& {\hat{c}}_j(\mathit{x})\le 0,j=1,\cdot \cdot \cdot ,N_c\\ & {\hat{g}}_i(\mathit{x})\le 0,i=1,\cdots ,N_t,\end{array}$(21) where ‘${}^{\wedge }$' stands for the surrogate model of the corresponding item of Eq. (19)(19) $\begin{array}{ll}\text{min}\text{. }& f(\mathit{x})\\ \text{w}\text{.r}\text{.t}\text{.}& {\mathit{x}}_{low}\le \mathit{x}\le {\mathit{x}}_{up}\\ \text{s}\text{.t}\text{. }& c_j(\mathit{x})\le 0,j=1,\cdot \cdot \cdot ,N_c\\ & g_i(\mathit{x})=g{p}_{i,0}(\mathit{x})-g{p}_i(\mathit{x})\le 0,i=1,\cdots ,N_t\end{array}$(19) .

The probability of feasibility of a sample point, satisfying all constraints, is the product of the PoF of each constraint, which is calculated using the surrogate models of the constraint functions. Then, the new updates can be obtained via solving the sub-optimization problem defined by the ISC.

4.3. The KS method

If the surrogate model of the aggregated constraint is built, the optimization problem becomes (22) $\begin{array}{ll}min& \hat{f}(\mathit{x})\\ \text{s}\text{.t}\text{.}& {\hat{g}}_{KS}(\mathit{x})\le 0\\ & \hat{\mathit{c}}(\mathit{x})\le 0,i=1,\cdots ,N_c,\\ & {\mathit{x}}_{\text{l}ow}\le \mathit{x}\le {\mathit{x}}_{up}\end{array}$(22) where $g_{KS}$ denotes the aggregated constraint of original constraints $g_i(i=1,\cdots ,N_t)$ in Eq. (19)(19) $\begin{array}{ll}\text{min}\text{. }& f(\mathit{x})\\ \text{w}\text{.r}\text{.t}\text{.}& {\mathit{x}}_{low}\le \mathit{x}\le {\mathit{x}}_{up}\\ \text{s}\text{.t}\text{. }& c_j(\mathit{x})\le 0,j=1,\cdot \cdot \cdot ,N_c\\ & g_i(\mathit{x})=g{p}_{i,0}(\mathit{x})-g{p}_i(\mathit{x})\le 0,i=1,\cdots ,N_t\end{array}$(19) . The subscript denotes the KS function, which could aggregate the constraints into one, is defined as (23) $\begin{array}{r}KS[\mathit{g}(\mathit{x})]=g_{max}(\mathit{x})+\frac{1}{\rho }\ln \left[\sum _{j=1}^n{e}^{\rho (g_j(\mathit{x})-g_{max}(\mathit{x}))}\right],\end{array}$(23) where ρ is the aggregation parameter. g_max(x) denotes the maximum value among all constraints at the current sample x, which is defined as (24) $g_{max}(\mathit{x})\text{ = }max({\mathit{g}}_i(\mathit{x})),i=1,\dots ,N_t.$(24) The range of the KS function at a sample point x is (25) $g_{max}(\mathit{x})<KS[\mathit{g}(\mathit{x})]<g_{max}(\mathit{x})+\frac{\ln n}{\rho }.$(25) From Eq. (25)(25) $g_{max}(\mathit{x})<KS[\mathit{g}(\mathit{x})]<g_{max}(\mathit{x})+\frac{\ln n}{\rho }.$(25) , one can observe that, (1) the upper bound value of the KS function decrease with the increase of ρ, and when ρ gets close to infinity, the value of the KS function just becomes equivalent to g_max(x); (2) If KS [g(x)] ≤ 0 is satisfied, g_max(x) must be less than zero, which may lead to a conservative evaluation of the original active constraints. Therefore, the feasible region defined by the KS function can be a subset of the original feasible region.

In the following single-variable example, we use the KS function to aggregate the inequality constraints. Consider the following convex inequality constraints: (26) $\begin{array}{l}{g}_1(x)=-0.5(x\text{ + }1)\le 0\\ g_2(x)=0.5x^2-1\le 0\\ g_3(x)=e^x-2\le 0.\end{array}$(26) The expression of the aggregation constraint is (27) $\begin{array}{rl}& KS[g_1(x),g_2(x),g_3(x)]\\ & =g_{max}(x)+\frac{1}{\rho }\ln \left[\sum _{j=1}^{\text{3}}{e}^{\rho (g_j(\mathit{x})-g_{max}(\mathit{x}))}\right],\end{array}$(27)

Figure 2 shows that the aggregated constraint changes with different values of ρ. Two numerical corners (corner 1 and 2) are located at where the constraints intersect in the feasible region, and the details of these two numerical corners on the right side of Figure 2. It can be observed clearly that the value of the KS function approaches the maximum constraint g_max(x) as ρ increases. If ρ is large enough, as ρ = 1000, the aggregated constraint is almost coinciding with the original active constraints.

Figure 2. Schematics of the aggregated constraint changing with different values of ρ.

4.4. Discussion of computational time of sub-optimization by using different constraint handling methods

To evaluate the contribution of the KS constraint method to the optimization efficiency, we compare the computational time of three constraints handling methods mentioned above by a 40-design-variable wing ASO case. The objective is to minimize the drag of the wing, 100 thickness constraints and another aerodynamic constraint are used to limit the wing internal volume and lift coefficient. To obtain new samples that satisfy all constraints, performing the sub-optimization is essential. In our case, GA is used for sub-optimization, with the maximum evolutional generation set as 500 and the population size set as 300. The computational cost of once sub-optimization performed by GA in each iteration of SBO can be expressed as (28) $t_{subopt}=n_{ME}\cdot n_{PS}\cdot (t^{\ast }\text{ + }{n}_{SM}\cdot t^{\mathrm{\prime }}),$(28) where $n_{ME},n_{PS}$ denote the maximum evolutional generation and the population size set for GA respectively. Their product $n_{ME}\cdot n_{PS}$ denotes the number of iterations of sub-optimization. $t^{\ast }$ is the computational cost of calculating all real geometric constraint values for once by generating a new aerodynamic shape. $n_{SM}$ denotes the number of surrogate models, and $t^{\mathrm{\prime }}$ is the computational cost of evaluating the objective or a constraint value for once by using corresponding surrogate models. $(t^{\ast }\text{ + }{n}_{SM}\cdot t^{\mathrm{\prime }})$ denotes the computational cost of evaluating the values of the objective and all constraints for once time.

The computational time of sub-optimization for optimizations by using different constraint handling methods at the 1000th sample are evaluated according to Eq. (28)(28) $t_{subopt}=n_{ME}\cdot n_{PS}\cdot (t^{\ast }\text{ + }{n}_{SM}\cdot t^{\mathrm{\prime }}),$(28) and list in Table 1. All above calculations are performed with Intel Xeon CPU e5-2690 at 2.6 GHz. From Table 1 one can see that, (1) for the optimization via EGCD method, $t^{\ast }\text{ = }0.35(s)$,$t^{\mathrm{\prime }}\text{ = }1.4\times {10}^{-3}(s)$. Although calculating 100 real geometric constraint values for once only takes 0.35s, the computational time of the sub-optimization via EGCD can be achieved according to Eq. (28)(28) $t_{subopt}=n_{ME}\cdot n_{PS}\cdot (t^{\ast }\text{ + }{n}_{SM}\cdot t^{\mathrm{\prime }}),$(28) , i.e. $t_{subopt}$ = 500 × 300 × (0.35 + 2 × 0.0014) =52,920 (s) = 14.70 (h); (2) For the optimization via BSMGC method, it takes about 26.5 minutes to train a surrogate model with 1000 samples. However, training 102 surrogate models using BSMGC costs up to 45.73 hours, and so many surrogate models result in 5.95 hours for once sub-optimization, i.e. $t_{subopt}$ = 500 × 300 ×(0 + 102 × 0.0014) =  21,420 (s) = 5.95 (h); (3) For the KS method, all geometric constraints are aggregated into one, so there are 3 surrogate models in total, and $t_{subopt}$ = 500 × 300 × (0 + 3 × 0.0014) =  630 (s)≈0.18 (h). In contrast, large benefit can be obtained by the proposed method. It uses KS function to aggregate all geometric constraints into one, and then the cheap surrogate model of the aggregated constraint is used to search for new shapes, which can avoid the large computational cost associated with building surrogate models for all geometric constraints and tremendous calculation of real geometric constraint values caused by sub-optimization. It is obvious that only the proposed method can be applied to the wing design instead of two other methods.

Table 1. Time cost of the wing design with different constraint handling methods at the 1000th sample.

	No. of surrogate models	Total No. of	Time cost of training	Time cost of sub-optimization (hours)
Methods	of geometric constraints	surrogate models	surrogate models (hours)	$t^{\ast }$(s)	$t^{\mathrm{\prime }}$(s)	Total cost (h)
EGCD	0	2	0.90	0.35	1.4 × 10^−3	14.70
BSMGC	100	102	45.73	0	1.4 × 10^−3	5.95
The KS method (This paper)	1	3	1.35	0	1.4 × 10^−3	0.18

EGCD: evaluating all geometric constraints directly.

BSMGC: building surrogate models for geometric constraints.

5. Examples and results

In this section, the proposed method is verified by G9 test case. Then, it is demonstrated by three typical CFD-driven aerodynamic shape optimization examples. First, the drag minimizations of the NACA0012 airfoil and ONERA M6 wing with a small number of geometric constraints are used to compare different geometric constraint handling methods. Second, the ONERA M6 wing optimization with 1000 geometric constraints is used to further study the ability of handling a large number of geometric constraints by the proposed method. Finally, drag minimization of NACA CRM wing in viscous transonic flow is performed to investigate the ability of the proposed method for handling various geometric constraints.

The optimizations are performed using the advanced TH-1A computing system at National Supercomputer Center (NSCC) in Tianjin. Each computing node in this cluster has two 2.6 GHz Intel Xeon CPUs (e5-2690 V4) and 128 GB RAM in total. Please note that, direct comparisons between SBO and other optimization methods are beyond the scope of this article, since we mainly concern about investigating a method that can efficiently handle geometric constraints within the framework of a SBO.

5.1. G9 test case

5.1.1. Problem statement

This is a benchmark case for strongly constrained global optimization. The mathematical model is (29) $\begin{array}{rl}Min.f(\mathbf{\text{x}})& =(x_1-10{)}^2+5(x_2-12{)}^2+{x_3}^4\\ & +3(x_4-11{)}^2+10{x_5}^6+7{x_6}^2\\ & +{x_7}^4-4x_6{x}_7-10x_6-8x_7\\ s.t.g_1& =-127+2{x_1}^2+3{x_2}^4+x_3+4{x_4}^2+5x_5\le 0\\ g_2& =-282+7x_1+3x_2+10{x_3}^2+x_4-x_5\le 0\\ g_3& =-196+23x_1+{x_2}^2+6{x_6}^2-8x_7\le 0\\ g_4& =4{x_1}^2+{x_2}^2-3x_1{x}_2+2{x_3}^2+5x_6-11x_7\le 0\\ {}^{}{x}_i& \in [-10,10]\text{, }i\text{ = }1\text{,}..\text{,}7.\end{array}$(29) This case contains 7 design variables and 4 constraints, in which g₁ and g₄ are active constraints. It is very difficult to optimize since the feasible region is only 0.5% of the design space and the order of the objective function value ranges from 10² to 10⁷. Although the real optimum is unknown, the optimization community has observed a best solution of $f(x^{\ast })=680.6300573$.

5.1.2. Aggregation parameter setting

As mentioned in the section 4, the approximation accuracy of the KS function, which influences the optimal results to a large extent, significantly depends on the aggregation parameter ρ. Therefore, it is necessary to find an appropriate aggregation parameter of the KS function for this case.

We select two sample points from the sampled datasets generated by DoE randomly. One is feasible that satisfy all constraint, while the other is infeasible that does not satisfy any constraint. Since the aggregated constraint KS (x, ρ) can provide a conservative approximation to the boundary of original constraints, we define the relative error as (30) $\text{RE}=\text{KS}(\text{x},\rho )-{\text{g}}_{\text{max}}(\text{x}),$(30) where g_max(x) is the maximum value of geometric constraint functions. Figure 3 shows the relative errors (REs) of both the feasible and infeasible samples get close to zero with the increase of ρ. The values of g_max and relative errors of the selected feasible and infeasible designs are shown in Table 2. One can see that the setting ρ = 1000 is reasonable for this case and the aggregation constraints are accurate enough.

Figure 3. Schematics of the relative errors changing with the increase of ρ.

Table 2. The maximum constraint g_max and relative errors for the selected feasible and infeasible design.

		Relative error
Description	gmax	ρ=50	ρ=100	ρ=1000
Feasible design	−0.11505	6.67 × 10^−3	1.45 × 10^−3	8.77 × 10^−12
Infeasible design	0.02825	1.72 × 10^−2	4.87 × 10^−3	1.18 × 10^−9

5.1.3. Results of optimizations

LHS method is used to generate 30 initial sample points. The ISC in conjunction with EI, MSP are applied to add 2 new sample points sequentially at each iteration of a SBO. In order to obtain new samples that satisfy geometric constraints, GA is used for sub-optimization. The maximum evolutional generation is set as 200 and the population size is set as 300. Since G9 case is difficult to optimize, the acceptable maximum number of sample points is set as 1200. To eliminate the randomness in the process of optimization, the optimization using each method of handling geometric constraints is repeated 5 times.

Figure 4 shows a typical convergence history of an optimization, in which the constant parameter ρ is set as 1000. The feasible and infeasible points are denoted by the blue squares and red triangles, respectively. The convergence history of the objective function is shown in Figure 5. The result of the optimization is listed in Table 3. It is observed that the value of objective function decreases as the increase of ρ. The best objective functional of 680.66319 is obtained when ρ = 1000 due to the accurate aggregated constraint, which is very close to the reference value of 680.63006. Besides, all constraints are satisfied strictly.

Figure 4. The initial DoE and infill-sampling process of a typical SBO.

Figure 5. Convergence histories of the optimizations.

Table 3. Results of G9 test case with the KS method.

Aggregation	Number		Surrogate prediction	Constraint values
method	of runs	Optimum	error at optimum (%)	Cons. 1	Cons.2	Cons. 3	Cons. 4
KS (ρ=50)	5	684.0511942	0.114414	−1.2661	−255.2084	−146.9092	−0.1057
KS (ρ=100)	5	682.6945800	0.008746	−0.0458	−252.8843	−142.7912	−0.5369
KS (ρ=1000)	5	680.6631931	0.027330	−1.4629	−260.2386	−162.4911	−9.5231

5.2. Aerodynamic shape optimization of NACA0012 airfoil in inviscid transonic flow

5.2.1. Problem statement

The design point of NACA0012 airfoil is that the angle of attack AoA = 0°and the freestream Mach number Ma = 0.85 in inviscid flow. The optimization problem is defined as (31) $\begin{array}{l}Min.C_d\\ s.t.C_l=0\\ \begin{array}{cc}y\ge y_{\text{baseline}}& \mathrm{\forall }x\in [0,1]\end{array},\end{array}$(31) where $C_l=0$ means the lift coefficient of the airfoil should be equal to zero in the whole optimization process. $y\ge y_{baseline}$ means the local thicknesses should be no less than the ones of the baseline geometry, at all the positions from leading edge to the trailing edge.

The NACA0012 airfoil is modified slightly, which leads to a zero-thickness trailing edge: (32) $\begin{array}{rl}{y}_{baseline}& =\pm 0.6(0.2969\sqrt{x}-0.1260x-0.3516x^2\\ & +0.2843x^3-0.1036x^4).\end{array}$(32) We use the free-form deformation method to parameterize the airfoil (see Figure 6). There are 19 control points on the upper and lower surface of the FFD control volume respectively. To ensure the airfoil shape to be symmetric, the displacements of the active control points on the upper surface are forced to be the same as their counterparts on the lower surface. The design variables are the y-coordinate of the active control points. As a result, there are 19 variables practically. In addition, as shown in Figure 6, 15 thickness constraints are distributed on NACA0012 airfoil and can be probed uniformly along the chord. The design space is shown in Figure 7. Please note that the optimization is solved in the fixed design space to avoid the failure of optimization due to the emergence of odd geometries.

Figure 6. Design variables and thickness constraints on NACA0012 airfoil.

Figure 7. Illustration of the design space.

5.2.2. Grid convergence study

The computational grid of NACA0012 airfoil with 51,200 cells (L2 grid) is generated by using our in-house O-type grid generator ‘2d-Grid'. To determine the resolution accuracy of this grid, a grid convergence study is performed. All grids (see Figure 8) with deficient size are calculated by our in-house CFD solver ‘PMNS3D'. Figure 9 shows that the pressure distributions computed using all grids are in good agreement with the reference data (Batina 1991), which turns out that the in-house CFD solver has sufficient accuracy. The simulation results are listed in Table 4, one can see that the L2 grid is of sufficient accuracy, the difference in the drag coefficient for the L2 grid and the finest grid (L0 grid) is only about one drag count.

Figure 8. O-type grids of NACA0012 airfoil for the grid convergence study.

Figure 9. Comparison of pressure distributions computed using all grids with reference data (Poole et al., 2015).

Table 4. Mesh-convergence study for the baseline NACA0012airfoil.

Grid level	Grid size	Cd (counts)	Cl	Cm
L0	204,800	480.20	0	0
L1	115,200	480.57	0	0
L2	51,200	481.27	0	0
L3	28,800	483.97	0	0
L4	12,800	488.24	0	0

5.2.3. Aggregation parameter setting

We select two sample points from the sampled datasets generated by DoE randomly. One is feasible, and the other is infeasible (see Figure 10) in the design space. Figure 11 shows how the REs of both the feasible and infeasible samples get close to zero with the increase of ρ. The values of g_max and REs of the selected feasible and infeasible designs are listed in Table 5. One can see that the setting ρ = 1000 is reasonable for our optimization problem and the aggregation constraints are accurate enough.

Figure 10. The selected feasible and infeasible designs compared with baseline airfoil.

Figure 11. Schematics of the relative errors changing with the increase of ρ.

Table 5. The maximum constraint g_max and relative errors of the selected feasible and infeasible designs.

Description	gmax	Relative error
Feasible design	0.11505	3.42 × 10^−7
Infeasible design	−0.02825	1.25 × 10^−9

Figure 12. The initial DoE and infill-sampling process of a typical SBO.

5.2.4. Results of optimizations

Initial surrogate models are built based on 50 sample points generated by the LHS method. The ISC in conjunction with EI, MSE, MSP, PI and LCB are applied to add 5 new sample points sequentially at each iteration. In order to obtain new samples that satisfy geometric constraints, GA is used for the sub-optimization. The maximum evolutional generation is set as 100 and the population size is set as 150. To compare different methods of handling geometric constraints, we set the acceptable maximum number of sample points to 600. Besides, to eliminate the randomness in the process of optimization, the optimization using each method of handling geometric constraints is repeated 5 times.

Figure 12 shows a typical convergence history of an optimization, in which the geometric constraints are calculated directly. The feasible and infeasible points are denoted by the blue squares and red triangles, respectively. The average convergence histories of the objective over time are shown in Figure 13. The results of optimizations are listed in Table 6. The time costs of NACA0012 airfoil designs using different constraint handling methods within the framework of a SBO are shown in Table 7. From Figure 13, Tables 6 and 7, one can see that,

1. The using of BSMGC results in the most time-consuming optimization, since the large number of surrogate models are built and used to predict new samples during the optimization. On one hand, the cost of training surrogate models for the objective and all geometric constraints becomes higher and higher as the number of sample points increases. On the other hand, predicting the values of the objective and constraints with so many surrogate models result in the high-cost sub-optimization. In addition, the optimization takes 51.55 hours, which is too expensive for a simple airfoil aerodynamic shape optimization driven by Euler-based CFD simulations. Therefore, BSMGC is only suitable for solving optimization problems with a very small number of geometric constraints.

2. The minimum C_d is obtained by EGCD, because the real values of geometric constraints are calculated directly and applied for the sub-optimization to select new samples. However, if it is expensive to solve the geometric constraint values, the application of this method will be limited due to the unacceptable cost of the sub-optimization solved by GA.

3. As ρ increases, the approximate accuracy of the KS aggregated constraint to the original geometric constraints improves. As a result, the average C_d decreases as ρ changes from 50 to 1000. When ρ = 1000, the accuracy of the KS function is high enough, both the minimum and the average C_d obtained by the KS method are very close to those obtained by EGCD.

4. The KS method reduces the number of surrogate models by aggregating all geometric constraints into one, and the aggregated constraint is used to select new samples that satisfy all geometric constraints during the sub-optimization, which makes it more efficient than the other methods. The average total time cost of the optimization via the KS method is 11.57 hours, which is only 22.44% and 35.11% of that via BSMGC and EGCD respectively.

Figure 13. Average convergence histories of aerodynamic shape optimization of NACA0012 airfoil.

Table 6. Results of NACA0012 airfoil design with different constraint handling methods.

		Cd (counts)			Surrogate prediction error
Methods	Number of runs	minimum	Maximum	average	at optimum (%)
BSMGC	5	66.84	71.32	68.82	0.981
EGCD	5	60.82	65.84	63.84	0.345
KS (ρ=50)	5	85.53	89.39	89.15	0.651
KS (ρ=100)	5	69.68	73.74	74.18	0.416
KS (ρ=1000)	5	62.87	65.57	64.88	0.985

BSMGC: building surrogate models for geometric constraints.

EGCD: evaluating all geometric constraints directly.

KS: the KS method.

Table 7. Time cost of NACA0012 airfoil design within the framework of a SBO (hours).

	Number of surrogate	Time cost of training	Time cost of	Total time cost
Methods	models	surrogate models	sub-optimization	of SBO
BSMGC	16	32.23	9.84	51.55
EGCD	1	0.41	23.06	32.95
KS	2	1.03	1.07	11.57

BSMGC: building surrogate models for geometric constraints.

EGCD: evaluating all geometric constraints directly.

KS: the KS method.

Figure 14. Comparison of NACA0012 airfoil and optimal airfoils.

Figure 15. Comparison of surface pressure distributions of NACA0012 airfoil with optimal airfoils.

The airfoil shapes and pressure distributions of the baseline airfoil and the best designs obtained from repeated optimizations are compared and sketched in Figures 14 and 15, respectively. One can see that the optimal airfoils become thicker at the leading and trailing edges compared with the baseline airfoil, the shock waves are significantly weakened. The comparison of Mach number contours of baseline and optimal airfoil achieved via the KS method with ρ = 1000 is shown in Figure 16, it is clearly shown that the strongly shock wave on the baseline airfoil surface is broken into an attached shock wave and a detached shock wave. Since the attached shock wave is much weaker, the shock wave drag is greatly reduced.

Figure 16. Comparison of Mach number contours of baseline and optimal airfoil.

The thicknesses of the optimal airfoils are compared with the baseline airfoil, as listed in Table 8. One can see that the thickness constraints marked in red are active constraints, as well as all constraints are satisfied strictly because the thicknesses of these optimal shapes are bigger than those of the baseline shape.

Table 8. Thicknesses of optimal airfoils compared with the baseline airfoil.

	Dimensionless
Index of	thicknesses of	Dimensionless thicknesses of optimal airfoils
thicknesses	baseline airfoil	BSMGC	EGCD	KS (ρ=50)	KS (ρ=100)	KS (ρ=1000)
Thick1	0.02434	0.04697	0.04557	0.04684	0.04670	0.04577
Thick2	0.08383	0.09402	0.09391	0.09527	0.09384	0.09322
Thick3	0.10569	0.10815	0.10727	0.10971	0.10827	0.10775
Thick4	0.11621	0.11629	0.11622	0.11893	0.11716	0.11660
Thick5	0.11989	0.12032	0.12005	0.12359	0.12151	0.12141
Thick6	0.11879	0.12110	0.12114	0.12459	0.12255	0.12346
Thick7	0.11412	0.11949	0.12125	0.12249	0.12106	0.12311
Thick8	0.10666	0.11723	0.12052	0.11882	0.11834	0.12106
Thick9	0.09699	0.11573	0.11862	0.11537	0.11588	0.11850
Thick10	0.08548	0.11457	0.11580	0.11287	0.11442	0.11574
Thick11	0.07239	0.11199	0.11230	0.11069	0.11279	0.11193
Thick12	0.05787	0.10682	0.10750	0.10696	0.10833	0.10630
Thick13	0.04196	0.09925	0.10062	0.09977	0.10011	0.09878
Thick14	0.02463	0.08715	0.08928	0.08882	0.08926	0.08849
Thick15	0.00574	0.03981	0.04264	0.04489	0.04377	0.04483

5.3. Aerodynamic shape optimization of ONERA M6 wing in inviscid transonic flow

5.3.1. Problem statement

In this example, drag minimization of ONERA M6 wing is performed in a fixed design space. This optimization problem subjects to lift and thickness constraints, the design condition is that the angle of attack AoA = 3.06°and the freestream Mach number Ma = 0.8395 in inviscid transonic flow. We apply the FFD method to parameterize the wing with a control volume which contains five control sections. Each control section contains five control points at the upper and lower surface respectively, as shown in Figure 17. Note that the control points (shown in black color) at the trailing edge are also fixed in this case, which results in 40 design variables in total. The drag minimization problem is formulated as (33) $\begin{array}{ll}Min.& C_D\\ s.t.& C_L\ge C_{L,0}\\ & Thicknes{s}_i\ge Thicknes{s}_{i,0}\text{,}i=1,\dots ,n.\end{array}$(33) The optimization problem is defined to subject to 20 (4 spanwise × 5 chordwise) full thickness constraints (see Figure 18).

Figure 17. FFD parameterization for ONERA M6 wing (40 design variables).

Figure 18. Geometric constraints enforced on the M6 wing (20 constraints).

5.3.2. Grid convergence study

The structured C–H grid with 1,245,816 cells (L2 grid) is generated by our in-house 3D grid generator ‘3D_CH_grid' for ONERA M6 wing. A grid convergence study is performed to determine the resolution accuracy of this grid. All grids (see Figure 19) are calculated by our in-house CFD solver ‘PMNS3D’. Figure 20 shows that the pressure distributions computed using all grids are in good agreement with the reference data (Batina, 1991), which turns out that the in-house 3D CFD solver has sufficient accuracy. The simulation results are listed in Table 9, one can see that the difference in the drag coefficient for the L2 grid and the finest grid (L0 grid) is only 1.2 drag counts, which turns out that the L2 grid is of sufficient accuracy.

Figure 19. C-H grids of ONERA M6 wing for the grid convergence study.

Figure 20. Comparison of pressure distributions computed using all grids with reference data (Batina, 1991) at two typical cross-sections.

Table 9. Mesh-convergence study for the baseline ONERA M6 wing.

Grid level	Grid size	CD (counts)	CL	CM
L0	2,056,392	119.66	0.2910	−0.2186
L1	1,617,408	119.92	0.2908	−0.2184
L2	1,245,816	120.84	0.2905	−0.2182
L3	479,232	124.70	0.2886	−0.2168
L4	202,176	132.47	0.2862	−0.2152

5.3.3. Aggregation parameter setting

The feasible and infeasible sample points are selected from the sampled datasets generated by DoE randomly. Figure 21 shows how REs get close to zero with the increase of ρ. Table 10 lists detail information about g_max and REs of the selected sample points. It is obvious that the aggregated constraints with ρ = 2000 is accurate enough for this optimization case.

Figure 21. Schematics of the relative errors changing with the increase of ρ.

Table 10. The maximum constraint g_max and relative errors for the selected feasible and infeasible designs.

Description	gmax	Relative error
Feasible design	5.03 × 10^−3	2.17 × 10^−6
Infeasible design	−0.30883	7.21 × 10^−11

5.3.4. Results of optimizations

For this case, 100 initial sample points are generated to build the initial surrogate models for the objective and constraint functions. EI, MSE, MSP, PI and LCB are used to select 5 new samples sequentially, the responses are solved in parallel to update the surrogate models. GA is used for sub-optimization to obtain new samples that satisfy such large number of geometric constraints, the maximum evolutional generation is set as 500 and the population size is set as 300. The maximum number of sample points is set as 750 and the optimization is repeated for 3 times to eliminate the randomness in the process of optimization.

A typical optimization convergence history is shown in Figure 22, in which the geometric constraints are handled by the KS method. The average convergence history of the objective over the same time are shown in Figure 23. The pressure distributions of the baseline and optimal shapes obtained from repeated optimizations are compared and sketched in Figure 24. The results of optimizations are listed in Table 11. The time costs of ONERA M6 wing designs using different constraint handling methods within the framework of a SBO are shown in Table 12. From Figures 23 and 24, Tables 11 and 12, one can see that, compared with the other geometric constraint handling methods, the KS method makes SBO achieve the preset maximum number of sample points first after 130 iterations, and the optimal wing shape with the minimum drag coefficient is obtained by the KS method within the same run time. It is obvious that without such an efficient geometric constraint handling method, performing a wing shape optimization that subjects to so many geometric constraints would be very intractable.

Figure 22. The initial DoE and infill-sampling process of a typical SBO.

Figure 23. Average convergence histories of wing aerodynamic shape optimization.

Figure 24. Comparison of pressure distributions and geometric shapes of the baseline and optimal wings.

Table 11. Results of M6 wing optimization with 20 geometric constraints.

	Number of	CD (counts)			Surrogate prediction
Methods	independent runs	Minimum	Maximum	Average	error at optimum
BSMGC	3	96.37	100.78	98.17	0.0576
EGCD	3	91.98	93.80	93.08	0.1687
KS	3	85.78	87.84	86.80	0.0019

BSMGC: building surrogate models for geometric constraints.

EGCD: evaluating all geometric constraints directly.

KS: the KS method.

Table 12. Time cost of M6 wing design within about 210 h.

			Time cost of	Time cost of
	Number of	Number of	training surrogate	sub-optimization
Methods	surrogate models	SBO iterations	models (hours)	(hours)
BSMGC	22	71	68.96	50.90
EGCD	2	26	0.57	175.86
KS	3	130	33.74	5.37

BSMGC: building surrogate models for geometric constraints.

EGCD: evaluating all geometric constraints directly.

KS: the KS method.

The thicknesses of the optimal wings are compared with the baseline, as listed in Table 13. We can see that most thickness constraints are active constraints. Moreover, all geometric constraints are satisfied strictly.

Table 13. Thicknesses of optimal wings compared with thebaseline.

		Dimensionless
Index of	Dimensionless	thicknesses of optimal wings
thicknesses	thicknesses of baseline	BSMGC	EGCD	KS
Thick1	0.03935	0.04016	0.04796	0.05493
Thick2	0.07489	0.07940	0.07979	0.08669
Thick3	0.07561	0.07781	0.07671	0.08455
Thick4	0.05133	0.05406	0.05374	0.05932
Thick5	0.01263	0.01317	0.01310	0.01356
Thick6	0.03564	0.04044	0.04387	0.03875
Thick7	0.06425	0.06613	0.07003	0.06744
Thick8	0.06490	0.06517	0.06519	0.06886
Thick9	0.04465	0.04705	0.04483	0.04876
Thick10	0.01245	0.01306	0.01266	0.01308
Thick11	0.03145	0.03146	0.03446	0.03151
Thick12	0.05361	0.05396	0.05812	0.05376
Thick13	0.05419	0.05600	0.05969	0.05625
Thick14	0.03797	0.04024	0.04181	0.04133
Thick15	0.01226	0.01271	0.01279	0.01302
Thick16	0.02683	0.02935	0.02914	0.02751
Thick17	0.04297	0.04907	0.05116	0.05185
Thick18	0.04346	0.04827	0.05212	0.05377
Thick19	0.03127	0.03289	0.03776	0.03875
Thick20	0.01199	0.01208	0.01362	0.01377

5.3.5. Further study of the KS method for handling a large number of geometric constraints

In order to explore the ability for handling large-scale geometric constraints of the KS method, the optimization of M6 wing is defined to subject to 1000 (20 panwise × 50 chordwise) full thickness constraints, which is shown in Figure 25.

Figure 25. 1000 thickness constraints (20 spanwise × 50 chordwise) for M6 wing design.

Two sample points are selected from the sampled datasets generated by DoE randomly. One is feasible and the other one is infeasible. Figure 26 shows how relative errors get close to zero with the increase of ρ. Table 14 lists the values of g_max and relative errors of the selected sample points. One can see that the aggregated constraints with ρ = 2000 is accurate enough for this optimization case.

Figure 26. Schematics of the relative errors changing with the increase of ρ.

Table 14. The maximum constraint g_max and relative errors of the selected feasible and infeasible designs.

Description	gmax	RE
Feasible design	−0.01334	7.03 × 10^−4
Infeasible design	0.40298	3.83 × 10^−9

For this case, the maximum number of sample points is set as 1000 to get a convergent optimization and the optimization is repeated for 3 times. The other optimization parameters for SBO are set the same as those mentioned in section 5.3.4. A typical optimization convergence history is shown in Figure 27. The average convergence history is shown in Figure 28. It can be seen that all optimizations are converged by using the KS method. The pressure distribution of the optimal wing shape is compared with the baseline, as shown in Figure 29. One can see that for the optimal wing shape, the shock is weakened significantly, as well as all constraints are satisfied.

Figure 27. The initial DoE and infill-sampling process of a typical SBO.

Figure 28. Average convergence history of wing aerodynamic shape optimization.

The optimization results are listed in Table 15. We can see that, compared with the baseline, large drag reduction is achieved. Please note that, since the number of the geometric constraints is much more than 20, the minimum objective value is slightly larger than that obtained via solving the problem with 20 geometric constraints by using the KS method. We analyze the results of the 20-constraint-problem and find that although all the geometric constraints are strictly satisfied, the thicknesses are smaller than the original shape at many locations that are not constrained. It is obvious that the KS method is an efficient geometric constraint handling method and can be used to solving aircraft design problems with large number of constraints.

Table 15. Optimization results of ONREA M6 wing optimization with 1000 geometricconstraints.

Number of geometric	Number of	CD (counts)			Surrogate prediction
constraints	independent runs	Minimum	Maximum	Average	error at optimum
1000	3	89.96	91.08	90.54	0.0017

5.4. Aerodynamic shape optimization of NASA CRM wing in viscous transonic flow

5.4.1. Problem statement

In this example, drag minimization of NASA CRM wing is performed in a fixed design space. The design condition is that the lift coefficient C_L= 0.5, the freestream Mach number Ma = 0.85, and the Reynolds number Re = 5 × 10⁶ in viscous transonic flow. The reference point of NASA CRM wing is at 25% mean aerodynamic chord, which corresponds to the position of the center of gravity. All coordinates are scaled by the mean aerodynamic chord (7.00532 m). The resulting reference length is 1.0, the half span is 3.758151 and the reference area is 3.407014. The coordinate for the reference point is at (x, y, z) = (1.2077, 0.007669, 0). The FFD method is used to parametrize the CRM wing. As shown in Figure 30, there are 5 control sections on the FFD volume of CRM wing, each control section has 18 control points, resulting 90 control points in total. In order to ensure that the locations of the leading edge and the trailing edge of the wing does not change, the two control points at the leading edge and trailing edge of each control section share one design variable to move the same distance along the opposite direction. We also choose y-direction displacements of the other 70 control points and the angle of attack as design variables. So there are 81 wing shape variables in total. The mathematical model of the optimization is (34) $\begin{array}{ll}Min.& C_D\\ s.t.& C_L=0.5\\ & C_{MZ}\ge -0.17\\ & t_i-0.85t_{i,0}\ge 0,i=0,\dots ,1400\\ & are{a}_j\ge are{a}_{j,0},j=0,\dots ,28\\ & v\ge v_0\\ & \Delta y_{\text{LE}}^{\text{upper}}\text{ = - }\Delta y_{\text{LE}}^{\text{lower}}\\ & \Delta y_{\text{TE}}^{\text{upper}}\text{ = - }\Delta y_{\text{TE}}^{\text{lower}}.\end{array}$(34) Three types of geometric constraints including thickness constraints, cross-sectional area constraints and the volume constraint are enforced on the wing (see Figure 31). Meanwhile, there are 1400 thickness constraints imposed at 50 chordwise and 28 spanwise locations of the wing, they are enforced to be no less than 85% of ones of the baseline wing at each location. The cross-sectional areas of these 28 wing spanwise locations are enforced to be no less than those of the baseline wing, as well as the volume is enforced to be no less than that of the baseline wing. As a result, there are 1429 geometric constraints in total. Another 2 aerodynamic constraints are that the lift coefficient of the wing is required to be equal to 0.5, and the pitch moment of the wing C_MZ is enforced to be no less than −0.17. The mathematical model of the optimization is

Figure 29. Comparison of pressure distribution of the baseline and the optimal wing.

Figure 30. FFD parameterization for NASA CRM wing (81 design variables).

Figure 31. Geometric constraints enforced on the NASA CRM.

5.4.2. Grid convergence study

The multi-block CFD structured grid with 2,442,240 cells (L1 grid) is generated by ANAYS ICEM CFD for NASA CRM wing. A grid convergence study is performed to determine the resolution accuracy of this grid. All grids (see Figure 32) are calculated by our in-house CFD solver ‘PMNS3D’. Figure 33 shows that the pressure distributions computed using all grids are in good agreement with the reference data (Kenway & Martins, 2016), which turns out that the in-house CFD solver has sufficient accuracy. The simulation results are listed in Table 16, one can see that the difference in the drag coefficient for the L1 grid and the finest grid (L0 grid) within 1 drag count, which turns out that the L1 grid has sufficient accuracy.

Figure 32. O-H grids of NACA CRM wing for the grid convergence study.

Figure 33. Comparison of pressure distributions computed using all grids with reference data (Kenway & Martins, 2016) at typical cross-sections.

Table 16. Mesh-convergence study for the baseline NASA CRM wing.

Grid level	Y^+	Grid size	CD (counts)	CL	CM
L0	2/3	6,101,760	205.23	0.5000	−0.1744
L1	4/9	2,442,240	205.90	0.5000	−0.1749
L2	1.0	1,243,264	213.07	0.5000	−0.1729
L3	1.2	306,336	246.38	0.5000	−0.1735

5.4.3. Aggregation parameter setting

The feasible and infeasible sample points are selected from the sampled datasets generated by DoE randomly. Figure 34 shows how REs get close to zero with the increase of ρ. Table 17 lists detail information about g_max and REs of the selected sample points. It is obvious that the aggregated constraints with ρ = 5000 is accurate enough for this optimization case.

Figure 34. Schematics of the relative errors changing with the increase of ρ.

Table 17. The maximum constraint g_max and relative errors for the selected feasible and infeasible designs.

Description	gmax	Relative error
Feasible design	8.96 × 10^−3	5.25 × 10^−6
Infeasible design	−1.42 × 10^−2	4.30 × 10^−10

5.4.4. Results of optimizations

For this case, 100 initial sample points are generated to build the initial surrogate models for the objective and constraint functions. EI and MSP are used to select 2 new samples sequentially, the responses are solved in parallel to update the surrogate models. GA is used for sub-optimization to obtain new samples that satisfy such large number of geometric constraints, the maximum evolutional generation is set as 500 and the population size is set as 300. The maximum number of sample points is set as 700 and the optimization is repeated for 2 times to eliminate the randomness in the process of optimization.

A typical optimization convergence history is shown in Figure 35, in which the geometric constraints are handled by the KS method. The average convergence history of the objective over the same time are shown in Figure 36. The pressure distributions of the baseline and optimal shapes obtained from repeated optimizations are compared and sketched in Figure 37. The results of optimizations are listed in Table 18. From Figure 37 and Table 18, one can see that, compared with the baseline, the shock is weakened significantly for the optimal wing shape, the drag coefficient drops from 205.9 counts to 197.8 counts, as well as all constraints are satisfied. It is obvious that the KS method is an efficient geometric constraint handling method and can be used to solving aircraft design problems with a large number and multi-types of geometric constraints.

Figure 35. The initial DoE and infill-sampling process of a typical SBO.

Figure 36. Average convergence histories of wing aerodynamic shape optimization.

Figure 37. Comparison of pressure distributions and geometric shapes of the baseline and optimal wing.

Table 18. Results of NASA CRM wing optimization with 3 types of geometric constraints.

Number of geometric	Number of	CD (counts)			Surrogate prediction
constraints	independent runs	Minimum	Maximum	Average	error at optimum
1429	2	197.80	197.87	197.84	0.0023

5. Conclusions

In this article, the constraint handling methods for solving CFD-driven aerodynamic shape optimization problems with geometric constraints were investigated within the framework of a SBO, in order to improve the capability of SBO for the engineering design optimizations. First, G9 numerical test case was performed to verify whether the KS method can handle the constraints effectively. Second, a comparative study of different methods of handling geometric constraints was presented via the aerodynamic shape optimizations of NACA0012 airfoil and M6 wing in transonic flows. Moreover, the M6 wing aerodynamic shape optimization with up to 1000 geometric constraints was conducted to explore the ability of the KS method for handling a large number of geometric constraints. Finally, drag minimization of NACA CRM wing in viscous transonic flow is performed to investigate the ability of the proposed method for handling various geometric constraints. According to the current test cases, following conclusions can be drawn as

1. The cost of SBO via building surrogate models for all geometric constraints (BSMGC) is very expensive, which induces a long design cycle. As the number of constraints increases, the cost becomes higher and higher, even prohibitive. Therefore, it is only suitable for solving optimization problems with very small number of constraints.

2. Evaluating all geometric constraints directly (EGCD) is not efficient, since for each of sub-optimizations using genetic algorithm (GA) the geometric modeling of aerodynamic configuration has to be conducted thousands of times or even more, and the total cost associated with a number of updating cycles can be extremely large. Therefore, the method EGCD will be limited due to aerodynamic shape optimization of simple configuration and with small number of updating cycles.

3. The proposed method of this article can be used to handle constraints effectively within the framework of a SBO by aggregating all inequality constraints into one. As ρ increases, the approximate accuracy of the KS constraint to the original geometric constraints improves. As a result, the KS method can dramatically improve the optimization efficiency of CFD-driven aerodynamic shape optimizations with the number of geometric constraints in the range from 15 to 1000 and the number of types of geometric constraints up to 3, which improves the usability of SBO in engineering design optimization and offers great potential for handling a larger number and more types of geometric constraints.

4. The parameter ρ needs to be selected according to a specific example, because the accuracy of the KS function may vary from a sample point to another and it is sensitive to the number of the geometric constraints. The method to obtain the optimal parameter ρ should be investigated to make the KS function has high accuracy for any sample.

The future work beyond the scope of this article will focus on improving the method of achieving the optimal parameter ρ to make the KS function has high accuracy for any sample, as well as the optimization of the complex aerodynamic shape such as wing-body-tail configuration with a large number and many types of geometric constraints.

Acknowledgement

This research is sponsored by the National Natural Science Foundation of China (11772261 and 11972305), as well as the Shaanxi Science Fund for Distinguished Young Scholars (2020JC-13) and Natural Science Fund of Shaanxi Province (2020JM-127). The authors would like to thank Dr. Yu Zhang for her linguistic assistance during the preparation of this article.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by National Natural Science Foundation of China: [Grant Number 11772261,11972305]; Shaanxi Science Fund for Distinguished Young Scholars: [Grant Number 2020JC-13]; Natural Science Fund of Shaanxi Province: [Grant Number 2020JM-127].