# An Efficient Surrogate-Based Optimization Method for BWBUG Based on Multifidelity Model and Geometric Constraint Gradients

Authors: Daiyu Zhang; Bei Zhang; Zhidong Wang; Xinyao Zhu
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Performing shape optimization of blended-wing-body underwater glider (BWBUG) can significantly improve its gliding performance. However, high-fidelity CFD analysis and geometric constraint calculation in traditional surrogate-based optimization methods are expensive. An efficient surrogate-based optimization method based on the multifidelity model and geometric constraint gradient information is proposed. By establishing a shape parameterized model, deriving analytical expression of geometric constraint gradient, constructing multifidelity surrogate model, the calculation times of high-fidelity CFD model and geometric constraints are reduced during the shape optimization process of BWBUG, which greatly improve the optimization efficiency. Finally, the effectiveness and efficiency of the proposed method are verified by performing the shape optimization of a BWBUG and comparing with traditional surrogate-based optimization methods.


---

## Body

## 1. Introduction

With the development of human society, the resources on land are exhausted. Because the ocean has unique advantages in resources, environment, and space, countries in the world are gradually paying attention to the development and utilization of marine resources. As a new type of underwater vehicle, Autonomous Underwater Glider (AUG) [1–4] has the advantages of long endurance, low manufacturing cost, and high reuse. It can sail in the ocean for a long time with high efficiency. At present, it has been widely used in marine environmental investigation and monitoring [5], such as sea temperature, salinity, ocean current speed, and marine environmental pollution. It can also be used to observe submarine volcanic eruptions, detect glaciers, and serve as data transmission nodes.The basic principle of the underwater glider is to convert the hydrodynamic lift generated by the wing into the driving force through adjusting its net buoyancy. Therefore, the hydrodynamic shape has an important impact on the overall performance and navigation distance of the underwater glider. The blended-wing-body underwater glider (BWBUG) [6–8] has a fuselage with a flat airfoil profile, and the hydrofoil and the fuselage are smoothly fused together. Compared with the traditional underwater glider composed of revolving body, hydrofoil, and control surface [9, 10], its outstanding feature is that it has a larger lift-drag ratio [11, 12], which leads to a longer distance in the horizontal direction when it descends the same depth in the water. In recent years, it has received more and more attention and research.To further improve the lift-drag ratio of the underwater glider, the optimization design method is widely applied in the shape design of the underwater glider. It is the top-level structure of the shape optimization design, which controls the operation of the whole optimization design and determines the development direction of the optimization. Generally, the optimization method can be divided into gradient-based optimization method [13], gradient-free optimization method [14], and surrogate-based optimization (SBO) method [15].The gradient-based optimization method mainly determines the search direction of the next iteration by the gradient information of objective function and constraint conditions. It is efficient in finding local minima for high-dimensional nonlinear problems defined by continuous smooth functions. The typical gradient-based optimization methods include the interior point method, trust region method, and sequence quadratic programming (SQP). The tests show that the computation complexity of the gradient-based optimization method is linearly proportional to the number of variables [16, 17]. However, these optimization algorithms require the user to provide the gradient information and to guarantee the function smoothness of constraints and objectives. For the optimization problems of BWBUG, gradients might not be available when the lift-drag ratio is calculated by the CFD method which is a black box for most users. Although gradients can always be approximated with finite differences, these approximations suffer from potentially significant inaccuracies (a truncation error of Oh or Oh2 when second-order). Therefore, the gradient-based optimization method is rarely found to be applied in the shape optimization design of BWBUG at this stage.Different from the gradient-based optimization methods, the gradient-free optimization method does not need gradient information. Therefore, it is useful when gradients are not available, such as when dealing with black-box functions [18]. Another major advantage of gradient-free methods is that they tend to be robust to numerical noise and discontinuity, making them easier to use than gradient-based methods. There is a wide variety of gradient-free methods which can perform a local or global search, use mathematical or heuristic criteria, and be deterministic or stochastic [19]. The gradient-free optimization methods mainly include the Nelder–Mead algorithm which is a local search algorithm based on heuristics, the generalized pattern search (GPS) which is based on mathematical criteria, and the evolutionary algorithms (such as the genetic algorithms and particle swarm optimization algorithms) which are based on the evolution of a population of designs. For the shape optimization design of BWBUG, Fu et al. [20] used NSGA2 (Nondominated Sorting Genetic Algorithm) to optimize the shape of the autonomous underwater glider. Tang et al. [21] used BESO (Bidirectional Evolutionary Structural Optimization) to optimize the shape of the underwater glider. However, the overall cost of gradient-free optimization is sensitive to the cost of the function evaluations because they require many iterations for convergence, and the number of iterations scales poorly with the number of design variables. It is found that the computational complexity of the gradient-free optimization methods increases exponentially with the number of variables [16, 17]. Therefore, the gradient-free optimization method is very time-consuming when it is applied to the shape optimization of BWBUG with a large number of variables.To reduce the computational cost in the optimization while maintaining the advantages of the gradient-free optimization method, the surrogate-based optimization method [15, 22] has been widely developed in recent years. It mainly uses the surrogate model to replace the time-consuming numerical simulation analysis and adds new sample points according to certain criteria in the optimization process to update the surrogate model cyclically. There are many methods to build the surrogate model such as polynomials, radial basis functions, and kriging. For the shape optimization design of BWBUG, Wang et al. [23] established the hydrodynamic surrogate model of flying wing underwater glider using Gaussian kernel function and optimized it by particle swarm optimization algorithm. Sun and Li et al. [24, 25] used the Kriging model to realize the hydrodynamic shape optimization of the underwater glider. Zhang and Li et al. [26, 27] applied the surrogate-based optimization method into the shape design of BWBUG and obtained an optimal shape design. These results show that the SBO method can effectively reduce the call times of computationally expensive functions in the optimization process, but one potential issue with the surrogate model is the curse of dimensionality [28]. With the increase of design variables, the CFD evaluations are significantly increasing to build the surrogate model. Considering the high-fidelity CFD evaluation is time-consuming for BWBUG, the SBO method still needs to be improved to reduce the computational cost when it is applied to the shape optimization design of BWBUG.To further improve the optimization efficiency of the surrogate model and explore the high lift-drag ratio characteristics of BWBUG, an efficient surrogate-based optimization method based on the multifidelity model and geometric constraint gradient information is proposed to reduce the calculation times of high-fidelity model and geometric constraints in the shape optimization process of BWBUG. Section2 describes the optimization model of BWBUG, and Section 3 introduces the problems existing in the traditional surrogate-based optimization method for the shape optimization design of BWBUG. Section 4 describes the proposed efficient surrogate-based method in detail. In Section 4, firstly, the parametric modeling of BWBUG is carried out, and the gradient analytical expression of geometric constraints is given. Secondly, the high-fidelity and low-fidelity CFD models of BWBUG are established, and the nested Kriging surrogate model is constructed. Finally, based on the criteria for infilling sample points, the shape optimization framework of BWBUG is built. Section 5 performs the shape optimization design of a BWBUG and compares the optimization results to demonstrate the effectiveness and efficiency of the proposed method. Section 6 summarizes the whole paper.


## 2. Optimization Model of BWBUG

For the BWBUG, its fuselage and wing are connected smoothly, and the cross section at each spanwise position is the airfoil section, which leads to a high lift-drag ratio. Figure1 describes the geometric shape of BWBUG, and the blue curves show the airfoil profiles at nine spanwise locations.Figure 1
The geometric shape of BWBUG.The main content of the shape optimization design for the BWBUG above is to carry out parametric modeling of geometric shape and to use the efficient optimization algorithm to obtain the optimal hydrodynamic performance on the premise of meeting certain geometric constraints. The specific mathematical optimization model of BWBUG can be expressed as follows:(1)Min−LDw.r.tx∈xl,xu,s.t.gx≥0,where L/D is the lift-drag ratio of the BWBUG, and the optimization objective is to maximize the lift-drag ratio. x is the shape parameterization variable with the variable range of xl,xu. gx denotes geometric constraints, which need to meet certain design requirements such as thickness and volume.


## 3. Problems Existing in the Traditional SBO Method

To solve the optimization model of BWBUG, the main idea of the traditional surrogate-based optimization method is to use the surrogate model to replace the time-consuming numerical simulation analysis in the optimization design. Its detailed optimization process can be described as follows:(1)
The Design of Experiments (DOE) is used to generate the sample points in the design space, and the response values at these points are obtained by running accurate numerical simulations. Then, the initial surrogate model is established based on these data.(2)
Based on the established surrogate model, according to certain criteria for adding new sample points, the global optimization algorithm is applied to solve the corresponding suboptimization problem. Finally, the optimal solution is obtained as the new sample points.(3)
The response values at the new sample points are calculated and the surrogate model is updated until the added sample point sequence converges to the optimal solution.For a better understanding, the whole optimization process is shown in Figure2.Figure 2
The optimization process of the traditional SBO method.Compared with the intelligent optimization method, the optimization efficiency can be greatly improved by establishing the surrogate model of lift-drag ratio and adopting the optimization process in Figure2, but the following two problems still need to be solved.Problem 1.
Figure2 mainly uses the high-fidelity CFD model to calculate the lift-drag ratio at different sample points and builds a surrogate model of the lift-drag ratio based on these data to replace the original CFD model for optimization design, which greatly reduces the calculation consumption in the optimization process. However, if all the sample data are calculated through the high-fidelity CFD model, the calculation is still very time-consuming.Problem 2.
Figure2 often uses gradient-free optimization algorithms (such as NSGA and PSO) to solve the suboptimization problems. Therefore, the calculation times of the geometric constraints have an exponential relationship with the number of optimization variables in solving the suboptimization problems. References [16–17] show that 8 optimization variables need to be calculated with thousands of geometric constraints. Even if it only takes ten seconds or a few seconds to regenerate the geometry each time, the amount of geometric constraint calculations in the whole suboptimization process is still considerable.


## 4. Efficient Surrogate-Based Optimization Method

When the surrogate model is used to optimize the shape of the BWBUG, if the high-fidelity CFD model is directly used for the hydrodynamic analysis, the calculation cost is high. In addition, the amount of geometric constraint calculation in the process of solving the suboptimization problem is large. In this section, the shape parametric model of the BWBUG is established, and the analytic expression of the geometric constraint gradient information is derived. A multifidelity CFD model composed of high-fidelity and low-fidelity CFD models is established to use the low-fidelity model to predict the results of the high-fidelity model. Finally, an efficient surrogate-based optimization method based on the multifidelity model and geometric constraint gradient information is proposed, which can greatly reduce the calculation times of the high-fidelity model and geometric constraints in the shape optimization process of BWBUG.


### 4.1. Shape Parameterization and Constraint Gradient Computation

#### 4.1.1. Shape Parametric Modeling

For the shape of BWBUG shown in Figure1, its parametric modeling can be built by interpolating all the airfoil sections with the lofting technique. Then, by modifying the profile, twist angle, and spatial position of each airfoil section in the parametric modeling, the shape of BWBUG can be deformed in the optimization design process. In the following section, the geometric parametric modeling process is described in detail, and the related parametric variables are defined.(1)
LetDK∈R2,k=0,…,m be the m data points of the initial airfoil, and arrange them in sequence with the airfoil trailing edge as the starting and ending points. Aξ is defined as the B-spline fitting curve of these m data points, which is determined by n control points Pi,i=1,…,n. By changing the value of Pi, the curve shape of Aξ can be changed. Figure 3 shows the data points of one airfoil profile, the B-spline fitting curve, and the corresponding control points.(2)
Place the airfoil sectionAξ on the xy plane and introduce the scale factors Sx and Sy in x direction and y direction, respectively. By changing the values of Sx and Sy, the chord length and thickness of the airfoil profile Aξ can be changed.(3)
Defineθ as the twist angle of airfoil section Aξ in the xy plane. Then, the rotation transformation of airfoil section Aξ can be realized by changing the value of θ.(4)
UseT to represent the placement position of airfoil section Aξ in three-dimensional space. Then, the spatial layout of airfoil section Aξ can be controlled by changing the value of T.To sum up, the shape of BWBUG can be expressed by the following formula:(2)Wξ,η=Tη+ΘηSηAξ,whereξ is the direction of airfoil profile curve and η is the spanwise direction. Sη is the 3×3 scale matrix and its diagonal element is Sx,Sy,0. Θη represents the rotation matrix with the twist angle of θ. Tη represents the position coordinate.Figure4 shows an example of the left-right symmetrical BWBUG. For ease of description, the left half is taken for display.By changing the parameter values ofPi,Sx,Sy,θ,T, the geometric deformation of the BWBUG can be realized in the optimization design process.Figure 3
Data points andAξ fitting curve of one airfoil.Figure 4
The parametric shape of BWBUG.


#### 4.1.2. Geometric Constraint Gradient Calculation

Because the BWBUG has internal loading requirements, it is often necessary to set corresponding thickness constraints at different positions. Suppose that the thicknessh between the upper and lower surfaces of the BWBUG is required to be greater than the initial thickness h0 at position x,0,z, the geometric constraint g can be expressed as follows:(3)g=Wξ1,η1−Wξ2,η2−h0≥0,where ξ1,η1 and ξ2,η2 are the local coordinates corresponding to the upper and lower surfaces of BWBUG at position x,0,z, respectively, and h=Wξ1,η1−Wξ2,η2.According to the expression formula (2) of Wξ,η, the gradient information of the geometric constraint g relative to the parameterized variables Pi,Sx,Sy,θ,T can be derived as follows:(4)dgdPi=Θη1Sη1∂Aξ1∂Pi−Θη2Sη2∂Aξ2∂Pi,dgdS=Θη1∂Sη1∂SAξ1−Θη2∂Sη2∂SAξ2,dgdθ=∂Θη1∂θSη1Aξ1−∂Θη2∂θSη2Aξ2,dgdT=∂Tη1∂T−∂Tη2∂T.


### 4.2. Multifidelity Surrogate Model

The lift-drag ratio of BWBUG is the primary objective of shape optimization, and its calculation by high-fidelity CFD model is the most time-consuming. Considering that the results of the low-fidelity CFD model can reflect the change trend of lift-drag ratio, the nested Kriging model of high-fidelity and low-fidelity models are established to significantly improve the efficiency of building the Kriging model under the condition of achieving the same accuracy.


#### 4.2.1. CFD Model with High and Low Fidelity

In this section, a multifidelity CFD model composed of the high-fidelity model and low-fidelity model is established. Among them, the difference of high and low-fidelity models is mainly reflected in the grid elements and the number of grids. For the grid generation of the left-right symmetrical BWBUG, a half shape is modeled to reduce the computation cost. The computation domain is set as a box-topology shape and the size of the domain is set to−10L,15L×0,10L×−10L,10L (L is the reference length) in length, width, and height direction, respectively. The trimmed mesher and prism layer mesher are used to generate the volume grids in the computational domain. The number of prism layers is set to 5, and the volume growth rate is set to 1.2. In addition, two control volumes are introduced to increase the local grid density near the BWBUG. The mesh number of the high-fidelity model is set about 350 W, and that of the low-fidelity model is set about 70 W. Figures 5 and 6 show the surface grids of the high-fidelity and low-fidelity model, respectively. It can be seen that the high-fidelity model has smaller grid elements and a larger number of grids, while the low-fidelity model has larger grid elements and a smaller number of grids.Figure 5
Surface grids of the high-fidelity model.Figure 6
Surface grids of the low-fidelity model.Both the high and low-fidelity models are solved using theSSTk−w turbulence model, and the same boundary conditions and solving parameters are set. The boundary conditions include the following: the left surface of the computation domain is set as the velocity inlet, the right surface is set as the pressure outlet, the back surface is set as the symmetry plane, and the other three remaining surfaces are set as the slip walls. In addition, the BWBUG surface is set to the nonslip wall. The solution parameters include the following: the maximum number of iteration steps is set to 1000 steps, and the convergence residual is set to 1e−5.Figures7 and 8 show the pressure distribution of high and low fidelity models at 1kn speed and 4° angle of attack. From the pressure distribution, it can be seen that the high-fidelity model exhibits smoother pressure distribution at the front edges of the BWBUG than the low-fidelity model.Figure 7
Pressure distribution of the high-fidelity model.Figure 8
Pressure distribution of the low-fidelity model.


#### 4.2.2. Nested Kriging Surrogate Model

The Kriging model [29, 30] is an interpolation model. The interpolation result y^x0 at the point x0 is defined as the linear weighting of the response values ys of the known samples:(5)y^x0=wTx0ys,where wx0 is the weighting coefficient.From formula (5), the lift-drag ratio at the estimated points can be obtained as long as the expression of the weighting coefficient wx0 is given. For the Kriging model, the unknown function is regarded as a specific realization of the Gaussian static random process, and the optimal weighting coefficient wx0 is calculated by both meeting the minimum mean square error and the unbiased condition. Based on [30], the analytical expression of wx0 is given by the following formula:(6)RFFT0wx0λ=rfx0,where R is called the correlation matrix, which is composed of correlation function values between all known sample points. r is called the correlation vector, which is composed of correlation function values between unknown points and all known sample points. f is the vector of basis function, which is previously defined.By solving the linear equations (6) and substituting them into formula (5), the estimated value of the Kriging model can be obtained as follows:(7)y^x0=fTx0β^+rTR−1ys−Fsβ^,where β^=FTRF−1FTR−1ys.It is assumed that the sampling points of high and low-fidelity CFD models are as follows:(8)S1=x11,…,x1nT,S2=x21,…,x2mT,where the subscripts “1” and “2” represent high-fidelity and low-fidelity CFD models, respectively, and n is often much smaller than m.Their corresponding lift-drag ratio calculation results are as follows:(9)Y1=y11,⋯,y1nT,Y2=y21,⋯,y2mT.Then, the nested Kriging model of high- and low-fidelity CFD models is established by the following steps.Step 1. Kriging model is established for the sample data of low-fidelity CFD model as follows:(10)y^lfx0=wlfTx0Y2,wherey^lf is the low-fidelity Kriging model.Step 2. The difference between high-fidelity value and low-fidelity Kriging model prediction value is calculated at the sample points of the high-fidelity CFD model, as shown in the following formula:(11)d=Y1−y^lfS1.Step 3. The Kriging model for the difference value calculated by formula (11) is established by the following formula:(12)d^x0=whfTx0d.Step 4. By summing the Kriging model established by formulas (10) and (12), the approximate results of the high-fidelity CFD model are obtained as follows:(13)y^hfx0=y^lfx0+d^x0,wherey^hf represents the approximate value of the high-fidelity CFD model.


### 4.3. Shape Optimization Framework for BWBUG

To perform the shape optimization of BWBUG efficiently, the nested Kriging model is established according to Section4.2 to replace the time-consuming CFD model to calculate the lift-drag ratio. The suboptimization problems for adding new sample points are solved efficiently by introducing the geometric constraint gradient information calculated in Section 4.1. The sample points are infilled dynamically, and the nested Kriging model is updated iteratively until the optimization converges. The detailed optimization framework based on the above idea is described as follows:(1)
The experimental design methods, such as full factorial design, fractional factorial, orthogonal array, central composite design, or Latin hypercube design [31, 32], are used to sample the parameterized design space, and the initial sample points of high-fidelity and low-fidelity models are obtained, respectively.(2)
At the corresponding sample points, the high-fidelity and low-fidelity CFD models are solved to calculate the lift-drag ratio of BWBUG, and the nested Kriging model is established according to the method in Section4.2.2.(3)
The criterion for infilling sample point is adopted and the corresponding suboptimization problem is solved to obtain the optimal solution as the new sample point. The common criteria include Minimum of Surrogate Prediction (MSP), Expected Improvement (EI), or Probability of Improvement (PI). [33, 34]. Among them, MSP is the simplest, most direct, and the earliest criterion to be used. Its principle is to infill the sample point which can make the objective function of the surrogate model reach the minimum value. For the BWBUG, MSP is selected, and it mainly infills sample points by solving the following suboptimization problems with geometric constraints:(14)Min.y^hfxw.r.tx∈xl,xu,s.t.gx≥0,wherey^hfx is the predicted value of the nested Kriging model.For the above suboptimization problem, considering that the gradient analytic expression of geometric constraintgx is given in Section 4.1.2, gradient optimization algorithms (such as the SQP method) are used to solve the suboptimization problem efficiently.(4)
At the new sample point generated by MSP, accurate CFD analysis is carried out to determine whether it converges to the optimal solution. If it converges, the optimization process ends. Otherwise, the results are added to the existing data set and the whole optimization process is repeated.For a better understanding of the process above, an overview of shape optimization algorithm for BWBUG is shown in Algorithm1, and the shape optimization framework for BWBUG is described in Figure 9.Algorithm 1:Shape optimization algorithm for BWBUG.
Inputs:n,m: number of initial samples for high- and low-fidelity CFD modelsxl,xu: lower and upper boundsτ: Minimum Expected ImprovementOutputs:x∗: best point identifiedy∗: corresponding function valueS1=samplen, S2=samplem           SampleY1=High_CFDS1, Y2=Low_CFDS2       Calculate by multifidelity CFD modelsWhilek<kmax and yhf−y∗/y∗<τdoy^lf=krigingS2,Y2              Build Kriging for high-fidelity modeld=Y1−y^lfS1                Calculate the difference valuesd^=krigingS1,d               Build Kriging for difference valuesy^hf=y^lf+d^                 Build the nested Kriging modelx∗=MSPx                 Solve the MSP problemyhf=High_CFDx∗, ylf=Low_CFDx∗     Evaluate at predicted optimumy∗=minyhf,Y1                Update best point if necessaryS1←x∗,S1, S2←x∗,S2           Add new point to training dataY1←yhf,Y1, Y2←ylf,Y2k←k+1End whileFigure 9
The shape optimization framework for BWBUG.


## 4.1. Shape Parameterization and Constraint Gradient Computation

### 4.1.1. Shape Parametric Modeling

For the shape of BWBUG shown in Figure1, its parametric modeling can be built by interpolating all the airfoil sections with the lofting technique. Then, by modifying the profile, twist angle, and spatial position of each airfoil section in the parametric modeling, the shape of BWBUG can be deformed in the optimization design process. In the following section, the geometric parametric modeling process is described in detail, and the related parametric variables are defined.(1)
LetDK∈R2,k=0,…,m be the m data points of the initial airfoil, and arrange them in sequence with the airfoil trailing edge as the starting and ending points. Aξ is defined as the B-spline fitting curve of these m data points, which is determined by n control points Pi,i=1,…,n. By changing the value of Pi, the curve shape of Aξ can be changed. Figure 3 shows the data points of one airfoil profile, the B-spline fitting curve, and the corresponding control points.(2)
Place the airfoil sectionAξ on the xy plane and introduce the scale factors Sx and Sy in x direction and y direction, respectively. By changing the values of Sx and Sy, the chord length and thickness of the airfoil profile Aξ can be changed.(3)
Defineθ as the twist angle of airfoil section Aξ in the xy plane. Then, the rotation transformation of airfoil section Aξ can be realized by changing the value of θ.(4)
UseT to represent the placement position of airfoil section Aξ in three-dimensional space. Then, the spatial layout of airfoil section Aξ can be controlled by changing the value of T.To sum up, the shape of BWBUG can be expressed by the following formula:(2)Wξ,η=Tη+ΘηSηAξ,whereξ is the direction of airfoil profile curve and η is the spanwise direction. Sη is the 3×3 scale matrix and its diagonal element is Sx,Sy,0. Θη represents the rotation matrix with the twist angle of θ. Tη represents the position coordinate.Figure4 shows an example of the left-right symmetrical BWBUG. For ease of description, the left half is taken for display.By changing the parameter values ofPi,Sx,Sy,θ,T, the geometric deformation of the BWBUG can be realized in the optimization design process.Figure 3
Data points andAξ fitting curve of one airfoil.Figure 4
The parametric shape of BWBUG.


### 4.1.2. Geometric Constraint Gradient Calculation

Because the BWBUG has internal loading requirements, it is often necessary to set corresponding thickness constraints at different positions. Suppose that the thicknessh between the upper and lower surfaces of the BWBUG is required to be greater than the initial thickness h0 at position x,0,z, the geometric constraint g can be expressed as follows:(3)g=Wξ1,η1−Wξ2,η2−h0≥0,where ξ1,η1 and ξ2,η2 are the local coordinates corresponding to the upper and lower surfaces of BWBUG at position x,0,z, respectively, and h=Wξ1,η1−Wξ2,η2.According to the expression formula (2) of Wξ,η, the gradient information of the geometric constraint g relative to the parameterized variables Pi,Sx,Sy,θ,T can be derived as follows:(4)dgdPi=Θη1Sη1∂Aξ1∂Pi−Θη2Sη2∂Aξ2∂Pi,dgdS=Θη1∂Sη1∂SAξ1−Θη2∂Sη2∂SAξ2,dgdθ=∂Θη1∂θSη1Aξ1−∂Θη2∂θSη2Aξ2,dgdT=∂Tη1∂T−∂Tη2∂T.


## 4.1.1. Shape Parametric Modeling

For the shape of BWBUG shown in Figure1, its parametric modeling can be built by interpolating all the airfoil sections with the lofting technique. Then, by modifying the profile, twist angle, and spatial position of each airfoil section in the parametric modeling, the shape of BWBUG can be deformed in the optimization design process. In the following section, the geometric parametric modeling process is described in detail, and the related parametric variables are defined.(1)
LetDK∈R2,k=0,…,m be the m data points of the initial airfoil, and arrange them in sequence with the airfoil trailing edge as the starting and ending points. Aξ is defined as the B-spline fitting curve of these m data points, which is determined by n control points Pi,i=1,…,n. By changing the value of Pi, the curve shape of Aξ can be changed. Figure 3 shows the data points of one airfoil profile, the B-spline fitting curve, and the corresponding control points.(2)
Place the airfoil sectionAξ on the xy plane and introduce the scale factors Sx and Sy in x direction and y direction, respectively. By changing the values of Sx and Sy, the chord length and thickness of the airfoil profile Aξ can be changed.(3)
Defineθ as the twist angle of airfoil section Aξ in the xy plane. Then, the rotation transformation of airfoil section Aξ can be realized by changing the value of θ.(4)
UseT to represent the placement position of airfoil section Aξ in three-dimensional space. Then, the spatial layout of airfoil section Aξ can be controlled by changing the value of T.To sum up, the shape of BWBUG can be expressed by the following formula:(2)Wξ,η=Tη+ΘηSηAξ,whereξ is the direction of airfoil profile curve and η is the spanwise direction. Sη is the 3×3 scale matrix and its diagonal element is Sx,Sy,0. Θη represents the rotation matrix with the twist angle of θ. Tη represents the position coordinate.Figure4 shows an example of the left-right symmetrical BWBUG. For ease of description, the left half is taken for display.By changing the parameter values ofPi,Sx,Sy,θ,T, the geometric deformation of the BWBUG can be realized in the optimization design process.Figure 3
Data points andAξ fitting curve of one airfoil.Figure 4
The parametric shape of BWBUG.


## 4.1.2. Geometric Constraint Gradient Calculation

Because the BWBUG has internal loading requirements, it is often necessary to set corresponding thickness constraints at different positions. Suppose that the thicknessh between the upper and lower surfaces of the BWBUG is required to be greater than the initial thickness h0 at position x,0,z, the geometric constraint g can be expressed as follows:(3)g=Wξ1,η1−Wξ2,η2−h0≥0,where ξ1,η1 and ξ2,η2 are the local coordinates corresponding to the upper and lower surfaces of BWBUG at position x,0,z, respectively, and h=Wξ1,η1−Wξ2,η2.According to the expression formula (2) of Wξ,η, the gradient information of the geometric constraint g relative to the parameterized variables Pi,Sx,Sy,θ,T can be derived as follows:(4)dgdPi=Θη1Sη1∂Aξ1∂Pi−Θη2Sη2∂Aξ2∂Pi,dgdS=Θη1∂Sη1∂SAξ1−Θη2∂Sη2∂SAξ2,dgdθ=∂Θη1∂θSη1Aξ1−∂Θη2∂θSη2Aξ2,dgdT=∂Tη1∂T−∂Tη2∂T.


## 4.2. Multifidelity Surrogate Model

The lift-drag ratio of BWBUG is the primary objective of shape optimization, and its calculation by high-fidelity CFD model is the most time-consuming. Considering that the results of the low-fidelity CFD model can reflect the change trend of lift-drag ratio, the nested Kriging model of high-fidelity and low-fidelity models are established to significantly improve the efficiency of building the Kriging model under the condition of achieving the same accuracy.


### 4.2.1. CFD Model with High and Low Fidelity

In this section, a multifidelity CFD model composed of the high-fidelity model and low-fidelity model is established. Among them, the difference of high and low-fidelity models is mainly reflected in the grid elements and the number of grids. For the grid generation of the left-right symmetrical BWBUG, a half shape is modeled to reduce the computation cost. The computation domain is set as a box-topology shape and the size of the domain is set to−10L,15L×0,10L×−10L,10L (L is the reference length) in length, width, and height direction, respectively. The trimmed mesher and prism layer mesher are used to generate the volume grids in the computational domain. The number of prism layers is set to 5, and the volume growth rate is set to 1.2. In addition, two control volumes are introduced to increase the local grid density near the BWBUG. The mesh number of the high-fidelity model is set about 350 W, and that of the low-fidelity model is set about 70 W. Figures 5 and 6 show the surface grids of the high-fidelity and low-fidelity model, respectively. It can be seen that the high-fidelity model has smaller grid elements and a larger number of grids, while the low-fidelity model has larger grid elements and a smaller number of grids.Figure 5
Surface grids of the high-fidelity model.Figure 6
Surface grids of the low-fidelity model.Both the high and low-fidelity models are solved using theSSTk−w turbulence model, and the same boundary conditions and solving parameters are set. The boundary conditions include the following: the left surface of the computation domain is set as the velocity inlet, the right surface is set as the pressure outlet, the back surface is set as the symmetry plane, and the other three remaining surfaces are set as the slip walls. In addition, the BWBUG surface is set to the nonslip wall. The solution parameters include the following: the maximum number of iteration steps is set to 1000 steps, and the convergence residual is set to 1e−5.Figures7 and 8 show the pressure distribution of high and low fidelity models at 1kn speed and 4° angle of attack. From the pressure distribution, it can be seen that the high-fidelity model exhibits smoother pressure distribution at the front edges of the BWBUG than the low-fidelity model.Figure 7
Pressure distribution of the high-fidelity model.Figure 8
Pressure distribution of the low-fidelity model.


### 4.2.2. Nested Kriging Surrogate Model

The Kriging model [29, 30] is an interpolation model. The interpolation result y^x0 at the point x0 is defined as the linear weighting of the response values ys of the known samples:(5)y^x0=wTx0ys,where wx0 is the weighting coefficient.From formula (5), the lift-drag ratio at the estimated points can be obtained as long as the expression of the weighting coefficient wx0 is given. For the Kriging model, the unknown function is regarded as a specific realization of the Gaussian static random process, and the optimal weighting coefficient wx0 is calculated by both meeting the minimum mean square error and the unbiased condition. Based on [30], the analytical expression of wx0 is given by the following formula:(6)RFFT0wx0λ=rfx0,where R is called the correlation matrix, which is composed of correlation function values between all known sample points. r is called the correlation vector, which is composed of correlation function values between unknown points and all known sample points. f is the vector of basis function, which is previously defined.By solving the linear equations (6) and substituting them into formula (5), the estimated value of the Kriging model can be obtained as follows:(7)y^x0=fTx0β^+rTR−1ys−Fsβ^,where β^=FTRF−1FTR−1ys.It is assumed that the sampling points of high and low-fidelity CFD models are as follows:(8)S1=x11,…,x1nT,S2=x21,…,x2mT,where the subscripts “1” and “2” represent high-fidelity and low-fidelity CFD models, respectively, and n is often much smaller than m.Their corresponding lift-drag ratio calculation results are as follows:(9)Y1=y11,⋯,y1nT,Y2=y21,⋯,y2mT.Then, the nested Kriging model of high- and low-fidelity CFD models is established by the following steps.Step 1. Kriging model is established for the sample data of low-fidelity CFD model as follows:(10)y^lfx0=wlfTx0Y2,wherey^lf is the low-fidelity Kriging model.Step 2. The difference between high-fidelity value and low-fidelity Kriging model prediction value is calculated at the sample points of the high-fidelity CFD model, as shown in the following formula:(11)d=Y1−y^lfS1.Step 3. The Kriging model for the difference value calculated by formula (11) is established by the following formula:(12)d^x0=whfTx0d.Step 4. By summing the Kriging model established by formulas (10) and (12), the approximate results of the high-fidelity CFD model are obtained as follows:(13)y^hfx0=y^lfx0+d^x0,wherey^hf represents the approximate value of the high-fidelity CFD model.


## 4.2.1. CFD Model with High and Low Fidelity

In this section, a multifidelity CFD model composed of the high-fidelity model and low-fidelity model is established. Among them, the difference of high and low-fidelity models is mainly reflected in the grid elements and the number of grids. For the grid generation of the left-right symmetrical BWBUG, a half shape is modeled to reduce the computation cost. The computation domain is set as a box-topology shape and the size of the domain is set to−10L,15L×0,10L×−10L,10L (L is the reference length) in length, width, and height direction, respectively. The trimmed mesher and prism layer mesher are used to generate the volume grids in the computational domain. The number of prism layers is set to 5, and the volume growth rate is set to 1.2. In addition, two control volumes are introduced to increase the local grid density near the BWBUG. The mesh number of the high-fidelity model is set about 350 W, and that of the low-fidelity model is set about 70 W. Figures 5 and 6 show the surface grids of the high-fidelity and low-fidelity model, respectively. It can be seen that the high-fidelity model has smaller grid elements and a larger number of grids, while the low-fidelity model has larger grid elements and a smaller number of grids.Figure 5
Surface grids of the high-fidelity model.Figure 6
Surface grids of the low-fidelity model.Both the high and low-fidelity models are solved using theSSTk−w turbulence model, and the same boundary conditions and solving parameters are set. The boundary conditions include the following: the left surface of the computation domain is set as the velocity inlet, the right surface is set as the pressure outlet, the back surface is set as the symmetry plane, and the other three remaining surfaces are set as the slip walls. In addition, the BWBUG surface is set to the nonslip wall. The solution parameters include the following: the maximum number of iteration steps is set to 1000 steps, and the convergence residual is set to 1e−5.Figures7 and 8 show the pressure distribution of high and low fidelity models at 1kn speed and 4° angle of attack. From the pressure distribution, it can be seen that the high-fidelity model exhibits smoother pressure distribution at the front edges of the BWBUG than the low-fidelity model.Figure 7
Pressure distribution of the high-fidelity model.Figure 8
Pressure distribution of the low-fidelity model.


## 4.2.2. Nested Kriging Surrogate Model

The Kriging model [29, 30] is an interpolation model. The interpolation result y^x0 at the point x0 is defined as the linear weighting of the response values ys of the known samples:(5)y^x0=wTx0ys,where wx0 is the weighting coefficient.From formula (5), the lift-drag ratio at the estimated points can be obtained as long as the expression of the weighting coefficient wx0 is given. For the Kriging model, the unknown function is regarded as a specific realization of the Gaussian static random process, and the optimal weighting coefficient wx0 is calculated by both meeting the minimum mean square error and the unbiased condition. Based on [30], the analytical expression of wx0 is given by the following formula:(6)RFFT0wx0λ=rfx0,where R is called the correlation matrix, which is composed of correlation function values between all known sample points. r is called the correlation vector, which is composed of correlation function values between unknown points and all known sample points. f is the vector of basis function, which is previously defined.By solving the linear equations (6) and substituting them into formula (5), the estimated value of the Kriging model can be obtained as follows:(7)y^x0=fTx0β^+rTR−1ys−Fsβ^,where β^=FTRF−1FTR−1ys.It is assumed that the sampling points of high and low-fidelity CFD models are as follows:(8)S1=x11,…,x1nT,S2=x21,…,x2mT,where the subscripts “1” and “2” represent high-fidelity and low-fidelity CFD models, respectively, and n is often much smaller than m.Their corresponding lift-drag ratio calculation results are as follows:(9)Y1=y11,⋯,y1nT,Y2=y21,⋯,y2mT.Then, the nested Kriging model of high- and low-fidelity CFD models is established by the following steps.Step 1. Kriging model is established for the sample data of low-fidelity CFD model as follows:(10)y^lfx0=wlfTx0Y2,wherey^lf is the low-fidelity Kriging model.Step 2. The difference between high-fidelity value and low-fidelity Kriging model prediction value is calculated at the sample points of the high-fidelity CFD model, as shown in the following formula:(11)d=Y1−y^lfS1.Step 3. The Kriging model for the difference value calculated by formula (11) is established by the following formula:(12)d^x0=whfTx0d.Step 4. By summing the Kriging model established by formulas (10) and (12), the approximate results of the high-fidelity CFD model are obtained as follows:(13)y^hfx0=y^lfx0+d^x0,wherey^hf represents the approximate value of the high-fidelity CFD model.


## 4.3. Shape Optimization Framework for BWBUG

To perform the shape optimization of BWBUG efficiently, the nested Kriging model is established according to Section4.2 to replace the time-consuming CFD model to calculate the lift-drag ratio. The suboptimization problems for adding new sample points are solved efficiently by introducing the geometric constraint gradient information calculated in Section 4.1. The sample points are infilled dynamically, and the nested Kriging model is updated iteratively until the optimization converges. The detailed optimization framework based on the above idea is described as follows:(1)
The experimental design methods, such as full factorial design, fractional factorial, orthogonal array, central composite design, or Latin hypercube design [31, 32], are used to sample the parameterized design space, and the initial sample points of high-fidelity and low-fidelity models are obtained, respectively.(2)
At the corresponding sample points, the high-fidelity and low-fidelity CFD models are solved to calculate the lift-drag ratio of BWBUG, and the nested Kriging model is established according to the method in Section4.2.2.(3)
The criterion for infilling sample point is adopted and the corresponding suboptimization problem is solved to obtain the optimal solution as the new sample point. The common criteria include Minimum of Surrogate Prediction (MSP), Expected Improvement (EI), or Probability of Improvement (PI). [33, 34]. Among them, MSP is the simplest, most direct, and the earliest criterion to be used. Its principle is to infill the sample point which can make the objective function of the surrogate model reach the minimum value. For the BWBUG, MSP is selected, and it mainly infills sample points by solving the following suboptimization problems with geometric constraints:(14)Min.y^hfxw.r.tx∈xl,xu,s.t.gx≥0,wherey^hfx is the predicted value of the nested Kriging model.For the above suboptimization problem, considering that the gradient analytic expression of geometric constraintgx is given in Section 4.1.2, gradient optimization algorithms (such as the SQP method) are used to solve the suboptimization problem efficiently.(4)
At the new sample point generated by MSP, accurate CFD analysis is carried out to determine whether it converges to the optimal solution. If it converges, the optimization process ends. Otherwise, the results are added to the existing data set and the whole optimization process is repeated.For a better understanding of the process above, an overview of shape optimization algorithm for BWBUG is shown in Algorithm1, and the shape optimization framework for BWBUG is described in Figure 9.Algorithm 1:Shape optimization algorithm for BWBUG.
Inputs:n,m: number of initial samples for high- and low-fidelity CFD modelsxl,xu: lower and upper boundsτ: Minimum Expected ImprovementOutputs:x∗: best point identifiedy∗: corresponding function valueS1=samplen, S2=samplem           SampleY1=High_CFDS1, Y2=Low_CFDS2       Calculate by multifidelity CFD modelsWhilek<kmax and yhf−y∗/y∗<τdoy^lf=krigingS2,Y2              Build Kriging for high-fidelity modeld=Y1−y^lfS1                Calculate the difference valuesd^=krigingS1,d               Build Kriging for difference valuesy^hf=y^lf+d^                 Build the nested Kriging modelx∗=MSPx                 Solve the MSP problemyhf=High_CFDx∗, ylf=Low_CFDx∗     Evaluate at predicted optimumy∗=minyhf,Y1                Update best point if necessaryS1←x∗,S1, S2←x∗,S2           Add new point to training dataY1←yhf,Y1, Y2←ylf,Y2k←k+1End whileFigure 9
The shape optimization framework for BWBUG.


## 5. Results and Discussion

In this paper, a left-right symmetrical BWBUG is taken as an example to perform the shape optimization design for verifying the effectiveness and efficiency of the proposed surrogate-based optimization method.


### 5.1. Optimization Problem Description

The initial shape of the left-right symmetrical BWBUG is composed of 11 airfoil sections, and all the airfoil sections are NACA0012 airfoil. In order to reduce the number of parametric variables and the amount of calculation in the optimization process, the left half shape composed of six airfoil sections is selected for optimal design.Considering that the internal structure and loading equipment of the BWBUG have been determined with the initial shape design, the chord length of each airfoil section and the placement position of each airfoil section are set as constant to ensure that the internal structure and equipment layout remain unchanged in the optimization design. The scale factors iny direction Sy and the twist angle θ of six airfoil sections are selected as optimization design variables to improve the performance of BWBUG.In addition, the thicknesses at1/4 chord position of 6 airfoil sections are set to be no less than the initial value as geometric constraints to meet the loading requirements of internal equipment. And the interval constraints of optimization variables are also introduced not only to ensure sufficient optimization exploration space but also to avoid a sharp increase in the amount of calculation caused by an excessively large interval.To sum up, the following optimization problems are obtained:(15)Min.−LDw.r.tSyi,θi,i=1,…,6s.t.0.8≤SyiSyi0≤1.2,−80≤θi≤80,hi−hi0≥0,i=1,…,6,where Syi0 is the initial value of the variable Syi and hi0 is the initial thickness value of BWBUG at 1/4 chord position of the ith airfoil section.


### 5.2. Shape Optimization and Results Analysis

The proposed optimization method is used to solve the optimization problem (15). Firstly, the Latin hypercube method is used to select 30 and 70 sample points as the initial sample points of high-fidelity and low-fidelity CFD models, respectively. Secondly, the corresponding lift-drag ratio is calculated at each sample point, and the nested Kriging model is built according to Section 4.2.2. Finally, the gradient-based optimization algorithm is used to solve the MSP suboptimization problem and the obtained optimal points are infilled dynamically to update the nested Kriging model until the optimization converges. The convergence condition is that the difference between the two optimal sample points is less than 10−4.The shape optimization was carried out in the workstation, whose processor is Intel® Core™ i7-7700K and memory size is 16 GB. The whole optimization converged after running 362 hours. In the process, 23 sample points were infilled dynamically to update the nested Kriging model. A total of 53 high-fidelity CFD analyses and 93 low-fidelity CFD analyses were performed to obtain the optimal shape for BWBUG. Figure10 shows the convergence history of the whole optimization process.Figure 10
The convergence history of shape optimization for the BWBUG.Table1 lists all the initial values from which the optimization was started, and all the optimal values obtained after the shape optimization of BWBUG. Figure 11 compares the pressure distribution on the upper surface between the initial shape and the optimal shape of BWBUG. Figure 12 compares the pressure distribution on the lower surface of the initial shape and the optimal shape.Table 1
The initial values and optimal results for BWBUG.
ParametersInitial valuesOptimal values−L/D−13.24−16.41Sy1/Sy101.01.0972Sy2/Sy201.00.8000Sy3/Sy301.01.1999Sy4/Sy401.01.2000Sy5/Sy501.00.8000Sy6/Sy601.00.8229θ10.0−5.4593θ20.0−3.0811θ30.0−7.3963θ40.0−3.5005θ50.0−3.9580θ60.0−0.0000h1−h100.06.0473h2−h200.00.0080h3−h300.04.2409h4−h400.03.6299h5−h500.01.4569h6−h600.00.0079Figure 11
Pressure distribution on the upper surface.Figure 12
Pressure distribution on the lower surface.It can be seen from Table1 that the thicknesses at 1/4 chord position of six airfoil sections are all greater than the initial values, which meets the requirements of geometric constraints. Compared with the initial shape of BWBUG, the lift-drag ratio has an increase of 23.9% from 13.24 to 16.41. Figures 11 and 12 give the reason why the lift-drag ratio increases. As Figure 11 shows, the low-pressure area on the upper surface of the optimal shape is larger than that of the upper surface of the initial shape, especially at the leading edge. While Figure 12 shows that the area of the high-pressure zone on the lower surface of the optimal shape is larger than that of the lower surface of the initial shape. Considering the total pressure distribution on the upper and lower surfaces, the optimal shape produces a higher lift force than the initial shape, which indirectly increases the high lift-drag ratio.Taken together, the proposed optimization method tries to find a better design shape of BWBUG at each iteration (Figure10) by dynamically infilling new sample points to update the nested Kriging model of the multifidelity CFD model, and finally, it obtains an optimal shape which has a better hydrodynamic distribution meeting the design requirements than the initial shape (Table 1). It is demonstrated that the proposed optimization is effective in improving the hydrodynamic performance of BWBUG.


### 5.3. Comparison with Traditional SBO Method

In this section, the traditional SBO method and the proposed optimization method are both applied to optimize the optimization problem (15), and their optimized results are compared with each other to verify the efficiency of the proposed optimization method.In order to make the optimized results statistically significant, five initial points are selected evenly within the variable bounds to start the optimization, as shown in Table2. As for the parameter settings of the proposed optimization method, they are the same as those in Section 5.2. As for the parameter settings of the traditional SBO method, the Latin hypercube method is also adopted to select 70 sample points to establish the initial surrogate model, but the high-fidelity CFD model is used to calculate the lift-drag ratio at all sample points. The remaining parameter settings are the same as those in Section 5.2.Table 2
Initial points from which the optimization is started.
CaseSy1/Sy10Sy2/Sy20Sy3/Sy30Sy4/Sy40Sy5/Sy50Sy6/Sy60θ1θ2θ3θ4θ5θ611.01.01.01.01.01.00.00.00.00.00.00.020.80.91.11.21.01.1−4.04.00.0−8.08.04.031.10.91.20.81.20.9−8.08.0−4.04.00.0−8.041.21.10.90.91.10.80.0−8.04.0−8.0−4.08.050.91.00.81.10.91.28.00.0−8.0−4.04.0−4.0In order to shorten the time to complete the whole optimization above, the ten optimizations above were performed on a massively parallel supercomputer and each optimization was carried out concurrently with the same computational resources. Table3 lists the detailed comparison data between the traditional SBO method and the proposed optimization method.Table 3
Comparison between the traditional SBO method and the proposed method.
CaseTraditional SBO methodProposed optimization methodL/DCalls to LFCalls to HFTime (h)L/DCalls to LFCalls to HFTime (h)116.400100267.116.419353191.2216.35093247.516.378848174.5316.420102272.916.439454194.5416.39098261.416.389252187.4516.37096256.116.409050181.1Avg.16.39098261.016.409151185.7∗LF denotes the low-fidelity CFD model; HF denotes the high-fidelity CFD model; Avg. represents the average value of the five cases; h represents hours.It can be seen from Table3 that the average values of L/D obtained by the traditional SBO method and the proposed optimization method are 16.39 and 16.40, respectively, which are roughly the same as a relative error of only 0.06%. And for each case of five initial points, the maximum relative error of L/D is 0.18%. All these results imply that the proposed optimization method can get almost the same optimal objective as the traditional SBO method. Considering that the traditional SBO method has been widely proved to be an effective optimization method [35], therefore the proposed optimization method is also proved to be effective indirectly.In addition, Table3 shows that the average time cost by the traditional SBO method and the proposed optimization method is 261.0 hours and 185.7 hours, respectively. Compared with the traditional SBO method, the average time cost by the proposed optimization method is reduced by 28.9%, which demonstrates that the proposed optimization method has higher optimization efficiency than the traditional SBO method. After the in-depth analysis, the calls to the low-fidelity CFD model and high-fidelity CFD model (Table 3) by the two methods reveal the reasons why the proposed optimization method is more efficient. For the traditional SBO method, just the high-fidelity CFD model is called for the calculation of L/D and the average calls are 98. While for the proposed optimization method, the average calls to the low-fidelity CFD model are 91 and the average calls to the high-fidelity CFD model are 51. The proposed optimization method calls 91 more low-fidelity CFD models and 47 less high-fidelity CFD models than the traditional SBO method. What is more, the computational time of the high-fidelity CFD model is almost 4–6 times that of the low-fidelity CFD model for the optimization problem (15). Therefore, the proposed optimization method mainly improves the optimization efficiency by using the low-fidelity CFD model to reduce the computational time of high-fidelity model, which agrees well with the hypotheses before.However, it should be noted in Table3 that the optimal values of L/D obtained from different initial points are different with a small deviation for both the proposed optimization method and the traditional SBO method. It implies that both methods cannot strictly converge to the global optimality. Although it is difficult and uneconomical for surrogate-based optimization methods to find a strict global optimal solution to complex engineering problems, it should be fully considered in our future work.


## 5.1. Optimization Problem Description

The initial shape of the left-right symmetrical BWBUG is composed of 11 airfoil sections, and all the airfoil sections are NACA0012 airfoil. In order to reduce the number of parametric variables and the amount of calculation in the optimization process, the left half shape composed of six airfoil sections is selected for optimal design.Considering that the internal structure and loading equipment of the BWBUG have been determined with the initial shape design, the chord length of each airfoil section and the placement position of each airfoil section are set as constant to ensure that the internal structure and equipment layout remain unchanged in the optimization design. The scale factors iny direction Sy and the twist angle θ of six airfoil sections are selected as optimization design variables to improve the performance of BWBUG.In addition, the thicknesses at1/4 chord position of 6 airfoil sections are set to be no less than the initial value as geometric constraints to meet the loading requirements of internal equipment. And the interval constraints of optimization variables are also introduced not only to ensure sufficient optimization exploration space but also to avoid a sharp increase in the amount of calculation caused by an excessively large interval.To sum up, the following optimization problems are obtained:(15)Min.−LDw.r.tSyi,θi,i=1,…,6s.t.0.8≤SyiSyi0≤1.2,−80≤θi≤80,hi−hi0≥0,i=1,…,6,where Syi0 is the initial value of the variable Syi and hi0 is the initial thickness value of BWBUG at 1/4 chord position of the ith airfoil section.


## 5.2. Shape Optimization and Results Analysis

The proposed optimization method is used to solve the optimization problem (15). Firstly, the Latin hypercube method is used to select 30 and 70 sample points as the initial sample points of high-fidelity and low-fidelity CFD models, respectively. Secondly, the corresponding lift-drag ratio is calculated at each sample point, and the nested Kriging model is built according to Section 4.2.2. Finally, the gradient-based optimization algorithm is used to solve the MSP suboptimization problem and the obtained optimal points are infilled dynamically to update the nested Kriging model until the optimization converges. The convergence condition is that the difference between the two optimal sample points is less than 10−4.The shape optimization was carried out in the workstation, whose processor is Intel® Core™ i7-7700K and memory size is 16 GB. The whole optimization converged after running 362 hours. In the process, 23 sample points were infilled dynamically to update the nested Kriging model. A total of 53 high-fidelity CFD analyses and 93 low-fidelity CFD analyses were performed to obtain the optimal shape for BWBUG. Figure10 shows the convergence history of the whole optimization process.Figure 10
The convergence history of shape optimization for the BWBUG.Table1 lists all the initial values from which the optimization was started, and all the optimal values obtained after the shape optimization of BWBUG. Figure 11 compares the pressure distribution on the upper surface between the initial shape and the optimal shape of BWBUG. Figure 12 compares the pressure distribution on the lower surface of the initial shape and the optimal shape.Table 1
The initial values and optimal results for BWBUG.
ParametersInitial valuesOptimal values−L/D−13.24−16.41Sy1/Sy101.01.0972Sy2/Sy201.00.8000Sy3/Sy301.01.1999Sy4/Sy401.01.2000Sy5/Sy501.00.8000Sy6/Sy601.00.8229θ10.0−5.4593θ20.0−3.0811θ30.0−7.3963θ40.0−3.5005θ50.0−3.9580θ60.0−0.0000h1−h100.06.0473h2−h200.00.0080h3−h300.04.2409h4−h400.03.6299h5−h500.01.4569h6−h600.00.0079Figure 11
Pressure distribution on the upper surface.Figure 12
Pressure distribution on the lower surface.It can be seen from Table1 that the thicknesses at 1/4 chord position of six airfoil sections are all greater than the initial values, which meets the requirements of geometric constraints. Compared with the initial shape of BWBUG, the lift-drag ratio has an increase of 23.9% from 13.24 to 16.41. Figures 11 and 12 give the reason why the lift-drag ratio increases. As Figure 11 shows, the low-pressure area on the upper surface of the optimal shape is larger than that of the upper surface of the initial shape, especially at the leading edge. While Figure 12 shows that the area of the high-pressure zone on the lower surface of the optimal shape is larger than that of the lower surface of the initial shape. Considering the total pressure distribution on the upper and lower surfaces, the optimal shape produces a higher lift force than the initial shape, which indirectly increases the high lift-drag ratio.Taken together, the proposed optimization method tries to find a better design shape of BWBUG at each iteration (Figure10) by dynamically infilling new sample points to update the nested Kriging model of the multifidelity CFD model, and finally, it obtains an optimal shape which has a better hydrodynamic distribution meeting the design requirements than the initial shape (Table 1). It is demonstrated that the proposed optimization is effective in improving the hydrodynamic performance of BWBUG.


## 5.3. Comparison with Traditional SBO Method

In this section, the traditional SBO method and the proposed optimization method are both applied to optimize the optimization problem (15), and their optimized results are compared with each other to verify the efficiency of the proposed optimization method.In order to make the optimized results statistically significant, five initial points are selected evenly within the variable bounds to start the optimization, as shown in Table2. As for the parameter settings of the proposed optimization method, they are the same as those in Section 5.2. As for the parameter settings of the traditional SBO method, the Latin hypercube method is also adopted to select 70 sample points to establish the initial surrogate model, but the high-fidelity CFD model is used to calculate the lift-drag ratio at all sample points. The remaining parameter settings are the same as those in Section 5.2.Table 2
Initial points from which the optimization is started.
CaseSy1/Sy10Sy2/Sy20Sy3/Sy30Sy4/Sy40Sy5/Sy50Sy6/Sy60θ1θ2θ3θ4θ5θ611.01.01.01.01.01.00.00.00.00.00.00.020.80.91.11.21.01.1−4.04.00.0−8.08.04.031.10.91.20.81.20.9−8.08.0−4.04.00.0−8.041.21.10.90.91.10.80.0−8.04.0−8.0−4.08.050.91.00.81.10.91.28.00.0−8.0−4.04.0−4.0In order to shorten the time to complete the whole optimization above, the ten optimizations above were performed on a massively parallel supercomputer and each optimization was carried out concurrently with the same computational resources. Table3 lists the detailed comparison data between the traditional SBO method and the proposed optimization method.Table 3
Comparison between the traditional SBO method and the proposed method.
CaseTraditional SBO methodProposed optimization methodL/DCalls to LFCalls to HFTime (h)L/DCalls to LFCalls to HFTime (h)116.400100267.116.419353191.2216.35093247.516.378848174.5316.420102272.916.439454194.5416.39098261.416.389252187.4516.37096256.116.409050181.1Avg.16.39098261.016.409151185.7∗LF denotes the low-fidelity CFD model; HF denotes the high-fidelity CFD model; Avg. represents the average value of the five cases; h represents hours.It can be seen from Table3 that the average values of L/D obtained by the traditional SBO method and the proposed optimization method are 16.39 and 16.40, respectively, which are roughly the same as a relative error of only 0.06%. And for each case of five initial points, the maximum relative error of L/D is 0.18%. All these results imply that the proposed optimization method can get almost the same optimal objective as the traditional SBO method. Considering that the traditional SBO method has been widely proved to be an effective optimization method [35], therefore the proposed optimization method is also proved to be effective indirectly.In addition, Table3 shows that the average time cost by the traditional SBO method and the proposed optimization method is 261.0 hours and 185.7 hours, respectively. Compared with the traditional SBO method, the average time cost by the proposed optimization method is reduced by 28.9%, which demonstrates that the proposed optimization method has higher optimization efficiency than the traditional SBO method. After the in-depth analysis, the calls to the low-fidelity CFD model and high-fidelity CFD model (Table 3) by the two methods reveal the reasons why the proposed optimization method is more efficient. For the traditional SBO method, just the high-fidelity CFD model is called for the calculation of L/D and the average calls are 98. While for the proposed optimization method, the average calls to the low-fidelity CFD model are 91 and the average calls to the high-fidelity CFD model are 51. The proposed optimization method calls 91 more low-fidelity CFD models and 47 less high-fidelity CFD models than the traditional SBO method. What is more, the computational time of the high-fidelity CFD model is almost 4–6 times that of the low-fidelity CFD model for the optimization problem (15). Therefore, the proposed optimization method mainly improves the optimization efficiency by using the low-fidelity CFD model to reduce the computational time of high-fidelity model, which agrees well with the hypotheses before.However, it should be noted in Table3 that the optimal values of L/D obtained from different initial points are different with a small deviation for both the proposed optimization method and the traditional SBO method. It implies that both methods cannot strictly converge to the global optimality. Although it is difficult and uneconomical for surrogate-based optimization methods to find a strict global optimal solution to complex engineering problems, it should be fully considered in our future work.


## 6. Conclusions

When the traditional SBO method is used for the shape optimization design of BWBUG, there are some problems such as time consumption of high-fidelity CFD analysis and large amount of calculation of geometric constraints in the suboptimization process. In this paper, a shape parametric model of the BWBUG is established, the analytic expression of geometric constraint gradient is derived, and the multifidelity CFD surrogate model is constructed. Finally, an efficient surrogate-based optimization method based on the multifidelity model and geometric constraint gradient information is proposed, which reduces the calculation times of the high-fidelity model and geometric constraints in the shape optimization process, and greatly improves the optimization efficiency.The shape of one left-right symmetrical BWBUG was optimized to verify the effectiveness and efficiency of the proposed optimization method. The results show that the lift-drag ratio of the optimal shape obtained by the proposed method is 23.9% higher than that of the initial shape, and the optimization efficiency of the proposed method is 28.9% higher than that of the traditional SBO method, which proves that the proposed method has a good application prospect in the shape optimization design of the BWBUG.In addition, the main theories of the proposed optimization method, such as the parametric modeling method, the gradient calculation method of geometric constraints, and the multifidelity surrogate model method, are not just suitable for the shape of BWBUG, but also for that of the aircraft, ships, cars, and so on. Therefore, the proposed optimization method has a broad applications prospect in the shape optimization design of aircraft, ships, cars, and so on.



---

*Source: 6939863-2021-08-24.xml*