Aerodynamic Optimization of Transonic Rotor Using Radial Basis Function Based Deformation and Data-Driven Differential Evolution Optimizer

Abstract

The complicated flow conditions and massive design parameters bring two main difficulties to the aerodynamic optimization of axial compressors: expensive evaluations and numerous optimization variables. To address these challenges, this paper establishes a novel fast aerodynamic optimization platform for axial compressors, consisting of a radial basic function (RBF)-based blade parameterization method, a data-driven differential evolution optimizer, and a computational fluid dynamic (CFD) solver. As a versatile interpolation method, RBF is used as the shape parameterization and deformation technique to reduce optimization variables. Aiming to acquire competitive solutions in limited steps, a data-driven evolution optimizer is developed, named the pre-screen surrogate model assistant differential evolution (pre-SADE) optimizer. Different from most surrogate model-assisted evolutionary algorithms, surrogate models in pre-SADE are used to screen the samples, rather than directly estimate them, in each generation to reduce expensive evaluations. The polynomial regression model, Kriging model, and RBF model are integrated in the surrogate model to improve the accuracy. To further save optimization time, the optimizer also integrates parallel task management programs. The aerodynamic optimization of a transonic rotor (NASA Rotor 37) is performed as the validation of the platform. A differential evolution (DE) optimizer and another surrogate model-assisted algorithm, committee-based active learning for surrogate model assisted particle swarm optimization (CAL-SAPSO), are introduced for the comparison runs. After optimization, the adiabatic efficiency, total pressure ratio, and surge margin are, respectively, increased by 1.47%, 1.0%, and 0.79% compared to the initial rotor. In the same limited steps, pre-SADE gets a 0.57% and 0.51% higher rotor adiabatic efficiency than DE and CAL-SAPSO, respectively. With the help of parallel techniques, pre-SADE and DE save half the optimization time compared to CAL-SAPSO. The results verify the effectiveness and the rapidity of the fast aerodynamic optimization platform.

1. Introduction

The rapid development of an advanced aero-engine proposes high requirements to the compressor. As one of the core components of the engine, an advanced compressor needs to meet the following three requirements: high single-stage pressure ratio, high adiabatic efficiency, and wide stability surge margin. These requirements bring challenges to the aerodynamic design of the compressor, especially the blade shaping. In the traditional design process, the blade shaping consists of the quasi-3-dimensional (Q3D) design and blade stacking process, followed with the 3-dimensional (3-D) optimization process. Generally, 100–150 design variables are needed to describe a single blade shape [1,2]. As for the multi-stage compressor, the number rises to 1000–2000 [2]. The 3-D viscous computational fluid dynamic (CFD) programs, which need massive computing resources, are widely used in aerodynamic optimization problems. Thus, two main difficulties need to be considered in the 3-D aerodynamic optimization of the axial compressor: numerous optimization variables and expensive evaluations. To overcome these barriers, 3-D blade parameterization methods and data-driven optimization algorithms were used and developed in the past few decades.

The 3-D parameterization methods can help reduce the number of design variables. To comprehensively describe the blade deformation, a parameterization method needs to cover the following characteristics: (1) have a wide range of deformations, (2) be robust to unique geometry, (3) guarantee the continuity and smoothness of the geometry, (4) be easy to implement in each deformation, and (5) reduce the total number of parameters. Traditional geometry parameterization methods of the compressor blade mainly focus on the 2-D blade sections using B-spline or Bezier spline. In the optimization process, the geometry deformation usually performs at the 2-D blade section before the stacking process. Combining with effective optimization algorithms and evaluation programs, the Q3D aerodynamic optimization methods were developed in the past few decades [3,4,5]. However, because a single blade is usually stacked by several sections, the parameterization of a 2-D section may still introduce too many variables in the optimization of a multi-stage axial compressor. Thus, full 3-D parameterization techniques, such as free-form-deformation (FFD) and radial basis functions (RBF), are applied to the compressor blade. These techniques can directly describe the 3-D geometry and perform deformation within a few parameters. FFD was first proposed by Sederberg and Parry [6], and widely used in external and internal flow aerodynamic optimization [7,8,9]. However, the control volumes need to have corresponding vertex points in three dimensions, and the number of control points in each layer must be equal, which leads to an increase in the number of control points. As the derivative method of the FFD, direct manipulation of free-form deformation (DFFD) can directly control the points on the blade and subsequently compute the deformations of the rest points. The control volume of the DFFD method can be placed arbitrarily in space to further reduce the parameterization points. Hang and Chen [10] performed the aerodynamic optimization of a 9-stage axial compressor with DFFD parameterization and got a significant performance improvement. RBF interpolation is another parameterization and deformation method which only needs to define the control points rather than the control volumes [11]. With the ability of direct and precise mesh deformation, RBF interpolation methods are widely used in the aerodynamic optimization of industry structures [12,13]. In the field of external flow, Biancolini et al. [14] adopted RBF to optimize the sail trim. Poirier and Nadarajah [15], as well as Jakobsson and Amoignon [16], implemented the aerodynamic optimization of the Onera M6 wing using RBF interpolation. For the internal flow of turbomachinery, Tang et al. [17] used RBF interpolation to perform deformation at the generated computation mesh to optimize the transonic rotor (NASA Rotor 67), aiming to save the mesh generation time in each evaluation. Khalfallah and Ghenaiet [18] combined the RBF parameterization method with sequential sampling to optimize the centrifugal impeller. However, in the turbomachinery optimization loop, it is usually difficult for the deformed meshes to meet the CFD calculation requirements, and it is difficult for the constraints to correspond to the parameterized points. Moreover, the parameterization process itself will be a heavy calculation burden due to the cloud of mesh points. Thus, Gagliardi and Giannakoglou [19] directly performed deformation on the wings. The RBF interpolation method is barely used in the direct deformation of a blade of an axial compressor, especially in the transonic blade. Considering the calculation resources and robustness of the methods, this paper uses the directly performed RBF interpolation method to implement the blade deformation. Adjoint methods provide another way to solve high-dimensional problems in the aerodynamic optimization of axial compressors [20]. The number of optimization variables is basically independent of the optimization time. However, the adjoint method needs to be supported by an open-source CFD solver to obtain the adjoint matrix, which brings difficulties in the combination with commercial CFD software. Furthermore, the adjoint method is a gradient-based algorithm, which is fast but easily falls into local optima. Thus, this paper mainly focuses on gradient-free algorithms.

Function evaluations (FE) are time-consuming in most 3-D aerodynamic optimizations of compressors. In recent years, many data-driven optimization algorithms were proposed due to the fast development of machine learning techniques. The surrogate model, which can quickly estimate the value of the objective function, is widely used in the amelioration of traditional evolutionary algorithms [21,22]. These algorithms are called surrogate model-assisted evolutionary algorithms (SAEAs). SAEAs first randomly generate a fully evaluated training database. Based on the database, a surrogate model is established to substitute the exact FEs in the next optimization loops. Then the algorithm will select limited samples for exact evaluating to append the database and improve the accuracy of the estimation model. The selections of the samples are based on unique strategies. Supported by a trained surrogate model, SAEA can save considerable optimization time. Zhen et al. [23] combined the surrogate model with the differential evolution algorithm to solve expensive optimization problems. Wang et al. [24] proposed committee-based active learning for surrogate-assisted particle swarm optimization (CAL-SAPSO), which can get competitive solutions while significantly reducing the number of exact FEs (100–300). However, these algorithms have limitations in optimization problems, which have narrow global optimums and many local optimums, which are blamed on the total substitution of the exact FEs. In each loop, the optimizer only gets FEs from the trained surrogate model. If the randomly initialized database fails to cover the global optimum, the optimizer may have difficulties in getting global optimum solutions and easily fall into the local optimum. Specific strategies, such as additional random coefficients, are introduced to strengthen the exploration. To overcome the challenges, this paper proposes a data-driven evolution optimization algorithm named the pre-screen surrogate model assistant differential evolution (pre-SADE) optimizer. The surrogate model is trained and used to pre-screen the FEs to save total optimization times. The polynomial regression (PR) model, Kriging model, and RBF model are used to construct the surrogate model. This algorithm is executed in the framework of a differential evolution optimizer, which is widely used in the industry and academic optimization problems. Different from most SAEAs, the samples are pre-screened and directly evaluated in each loop of pre-SADE. More exact FEs also help to avoid the algorithm falling into the local optimum. The algorithm is also upgraded to support parallel operations to save significant running time.

This paper first establishes a novel aerodynamic optimization platform for axial compressors, which is based on the RBF interpolation parameterization technique and pre-SADE optimizer. A transonic rotor (NASA Rotor 37) is optimized to verify the effectiveness of the platform. The remainder of this paper is organized as follows: Section 2 presents the key components of the optimization platform, including the deductions of the RBF interpolation method, as well as the diagrams and the estimations of pre-SADE algorithm. The detailed study case and the optimization settings are provided in Section 3. Section 4 gives the optimization results, detailed analyses of the optimized geometry, and the complex flow conditions. Finally, several conclusions are drawn in Section 5.

2. Methods

2.1. Parameterization Method

2.1.1. RBF Shape Parametrization Technique for Compressors

Buhmann [11] summarized and deduced RBF kernel and interpolation methods. Benefiting from the versatility and robustness, RBF interpolation was widely used in the real-time freeform shaping editing of graphics [25,26], as well as the mesh deformation [27]. The parameterization method is based on a defined map, $\mathscr{H}\left(\mathit{x}\right):{\mathit{R}}^{\mathit{n}}\to {\mathit{R}}^{\mathit{n}}$. The radial basis kernel is a real-valued function:

$$\phi \left(\mathit{x}\right)=\phi (\Vert \mathit{x}-{\mathit{x}}_i\Vert )$$

(1)

where $\Vert \mathit{x}-{\mathit{x}}_{\mathit{i}}\Vert$ denotes the Euclidean distance $r$ from point $\mathit{x}$ to point ${\mathit{x}}_{\mathit{i}}$. ${\mathit{x}}_{\mathit{i}}$ is the location of the RBF source node. The formulation shows that the value of the RBF kernel function only depends on the distances to the source nodes. To implement the translation and rotation deformations, a linear polynomial $p\left(\mathit{x}\right)$ is introduced to the formulation. Thus, the basic formulation of the map takes the form:

$$\mathscr{H}\left(\mathit{x}\right)=p\left(\mathit{x}\right)+{\displaystyle \sum }_{i=1}^N{\mathit{w}}_{\mathit{i}}{\phi }_{\mathit{i}}\left(\mathit{x}\right)$$

(2)

where $\mathscr{H}\left(\mathit{x}\right)$ is the interpolated function. The index $i$ identifies the source nodes of the RBF, ${\mathit{w}}_{\mathit{i}}$ is the weight of each kernel function. $N$ denotes the total number of control points. ${\mathit{\gamma }}_{\mathit{i}}$ indicates the polynomial coefficients of the linear polynomial $p\left(\mathit{x}\right)$.

To execute the deformation, the displacement ${\mathit{\delta }}_{\mathit{i}}$ of the source nodes (control points) needs to be interpolated to the origin scatters. ${\mathit{w}}_{\mathit{i}}$ needs to be solved through the equation of the displacement at source nodes:

$$\mathscr{H}\left(x_i\right)={\mathit{\delta }}_{\mathit{i}}$$

(3)

To close the equations, new requirements need to be introduced because the polynomial term has no contributions at source points.

$${\displaystyle \sum }_{i=1}^N{\mathit{w}}_{\mathit{i}}q\left(\mathit{x}\right)=0$$

(4)

The degrees of all polynomials $q\left(\mathit{x}\right)$ are equal or less than $p\left(\mathit{x}\right)$.

Thus, the coefficient ${\mathit{\gamma }}_{\mathit{i}}$ and the weight ${\mathit{w}}_{\mathit{i}}$ will be computed by solving the linear system.

$$\left[\begin{array}{cc}\mathit{M}& \mathit{P}\\ {\mathit{P}}^T& 0\end{array}\right]\left(\begin{array}{c}\mathit{W}\\ \mathit{\gamma }\end{array}\right)=\left(\begin{array}{c}\Delta \\ \mathbf{0}\end{array}\right)$$

(5)

where $\mathit{M}$ is the radial interpolation matrix of the source nodes,

$$\mathit{M}=\left[\begin{array}{ccc}{\phi }_1\left({\mathit{x}}_1\right)& \cdots & {\phi }_1\left({\mathit{x}}_N\right)\\ ⋮& \ddots & ⋮\\ {\phi }_N\left({\mathit{x}}_1\right)& \cdots & {\phi }_N\left({\mathit{x}}_N\right)\end{array}\right]$$

(6)

$\mathit{P}$ is the matrix introduced to close the linear system. In this 3D interpolation method, $\mathit{P}$ takes the form,

$$\mathit{P}=\left[\begin{array}{cc}1& {\mathit{x}}_1\\ \begin{array}{c}1\\ \begin{array}{c}⋮\\ 1\end{array}\end{array}& \begin{array}{c}{\mathit{x}}_2\\ \begin{array}{c}⋮\\ {\mathit{x}}_N\end{array}\end{array}\end{array} \begin{array}{c}{\mathit{y}}_1\\ \begin{array}{c}{\mathit{y}}_2\\ \begin{array}{c}⋮\\ {\mathit{y}}_N\end{array}\end{array}\end{array} \begin{array}{c}{\mathit{z}}_1\\ \begin{array}{c}{\mathit{z}}_2\\ \begin{array}{c}⋮\\ {\mathit{z}}_N\end{array}\end{array}\end{array}\right]$$

(7)

$\mathit{W}$, $\mathit{\gamma }$, $\Delta$ are formed by following matrix forms, respectively.

$$\mathit{W}=\left[\begin{array}{c}{\mathit{w}}_1^T\\ {\mathit{w}}_2^T\\ \begin{array}{c}⋮\\ {\mathit{w}}_N^T\end{array}\end{array}\right], \mathit{\gamma }=\left[\begin{array}{c}{\mathit{\gamma }}_1^T\\ {\mathit{\gamma }}_2^T\\ \begin{array}{c}⋮\\ {\mathit{\gamma }}_4^T\end{array}\end{array}\right], \Delta =\left[\begin{array}{c}{\mathit{\delta }}_1^T\\ {\mathit{\delta }}_2^T\\ \begin{array}{c}⋮\\ {\mathit{\delta }}_N^T\end{array}\end{array}\right]$$

(8)

Thus, the Equation (5) is a $\left(N+4\right)\times \left(N+4\right)$ linear system. By solving the linear system, $\mathit{W}$ and $\mathit{\gamma }$ are calculated and the interpolations are deducted.

2.1.2. The Parameterization and Deformation of Rotor 37 Based on RBF Interpolation

The implementation of RBF interpolation in this paper is based on PyGeM [28], which is a Python library that consists of several popular parameterization and deformation techniques. In the RBF interpolation method, the inputs consist of the original geometry, control points and their movements, the radial basis kernel functions, and the affect radius. The affect radius describes the influence area of the kernel functions. The values of the kernel function will be zero if the distance is larger than the affect radius, which leads to non-deformation outside the influence area in each source node. The value of the affect radius depends on the specific interpolation problems. Figure 1 shows the flow chart of the RBF interpolation method and the deformation samples of a single blade. The interpolation method first reads the original geometry. Then the number of control points and the appropriate locations are decided, mainly based on the input geometry, and the affect radius needs to be subsequently determined to balance the local deformation with the global smoothness of the geometry. With the determinations of the control points and affect radius, the parameterization process is finished. In the deformation process, the movements of the control points are introduced and interpolated to each point of the original geometry by solving the linear system. After updating the locations of all the points, the new geometry is generated. The local deformations (such as the precise control of the mean camber line) and the global deformations (such as the bow, sweep, and twist operations of the blade) can be intuitively controlled through the RBF interpolation method.

2.2. Optimization Algorithm

2.2.1. Pre-Screen Surrogate Model Assisted Differential Evolution Optimizer

Aiming to save the total predicting times, most SAEAs totally substitute the expensive evaluating program in the estimations of samples. However, these strategies may make the algorithms easily fall into the local optimum, causing the total substitution to not cover the global optimum. Moreover, to further save the optimization time, these algorithms may only provide one sample to evaluate at the same time, which cannot make full use of computing resources. With the fast development of computing resources, parallel calculating methods play an important role in optimization problems. With parallel techniques, several FEs can be evaluated at the same time, and the calculation of a single expensive FE can be significantly accelerated. Thus, as an improvement, pre-SADE slightly increases the number of exact FEs to overcome the limitations of the exploration without increasing the total optimization time.

Emmerich [29] proposed the metamodel-assisted evolution strategies in single- and multi-objective optimization. The FEs are pre-screened by the metamodel in each loop to save total optimization time. With the inspiration of pre-screen strategies and the integrated surrogate model, this paper proposes the pre-screen surrogate model-assisted differential evolution optimizer (pre-SADE). Figure 2 shows the diagram of the pre-SADE.

Pre-SADE runs in the framework of a differential evolution (DE) optimizer, which was proposed and details were deducted by Storn and Price [30]. The main operations of a DE consist of initialization, mutation, crossover, and selection. Apart from the initialization, new FEs are generated by mutation and crossover processes in each generation. The pre-screen strategy is used in the calculation of FEs. In the first generation, randomly initialized FEs need to be exactly evaluated. In each generation, the candidates, which need to be fully calculated, are put into the task pool. The management program (MP) automatically distributes the computation resources for parallel calculation. Every local parallel node calculates one candidate from the task pool independently and puts the result into the result pool. Each node contains several parallel cores. The evaluating process finishes until the task pool is empty. Based on the fully evaluated samples of the first generation, the global surrogate model (GSM) is established.

The framework of pre-SADE is based on Python, and a surrogate model toolbox [31] based on MATLAB, is used to implement the global surrogate model. The GSM consists of a PR model, Kriging model, and RBF model. The PR model ($s_{PR}$) is a quadratic polynomial model. The Kriging model ($s_{KR}$) is a simple Gaussian kernel Kriging model. The RBF model ($s_{RBF}$) is a RBF network with one hidden layer. The output ($s_{GSM}$) of the surrogate model is the weighted summary of the three models, which takes the form:

$$s_{GSM}\left(\mathit{x}\right)=w_1{s}_{PR}\left(\mathit{x}\right)+w_2{s}_{KR}\left(\mathit{x}\right)+w_3{s}_{RBF}\left(\mathit{x}\right)$$

(9)

where $w_i$, as the weight of the $i$ (1 ≤$i$ ≤ 3) surrogate model, is defined by

$$w_i=0.5-\frac{e_i}{2\left(e_1+e_2+e_3\right)}$$

(10)

where $e_i$ denotes the root mean square error of the $i \left(1\le i\le 3\right)$ surrogate model [32].

In the following generations, all the newly generated FEs are input to the GSM and quickly obtain estimations. Then the FEs are ranked by the estimated value. The exact evaluation ratio ($e_r)$ denotes the ratio of exact FEs to all FEs in each generation. The former $e_r$ samples are calculated and added to the training database to update the GSM. The evaluations of the rest samples are set as the estimated value. Finally, all FEs of the new generation are calculated and returned to the next loop of the DE optimizer. The algorithm terminates once reaching the max generation. The pseudo-code of the algorithm is provided in Appendix A.

2.2.2. The Verification of Pre-SADE

To examine the performance of pre-SADE, several benchmark optimization problems are introduced in this paper. The genetic algorithm (GA), DE, and CAL-SAPSO are employed as the comparison algorithm. The CAL-SAPSO needs 11$d$ steps to get acceptable solutions [33], where $d$ is the dimension of the optimization variables. Similarly, pre-SADE aims to get competitive solutions with hundreds of exact function evaluations. With the help of parallel techniques, the max step number is set to 15–20$d$ to save the total optimization time. The detailed optimization settings of the benchmark functions are shown in

Table 1. The settings of test problems.

Problem	$\mathit{d}$	$\mathbf{Search} \mathbf{Space} \left(\mathit{d}\right)$	Optimum	Max Step	Notes
Ackley	10, 15, 20	[−32.768, 32.768]	0	200, 200, 300	Multi-Modal
Griewank	10, 15, 20	[−600, 600]	0	200, 200, 300	Multi-Modal
Levy	10, 15, 20	[−10, 10]	0	200, 200, 300	Multi-Modal

. The Ackley function, Griewank function, and Levy function [34] are widely used for testing optimization algorithms. The Ackley function is characterized by a nearly flat outer region, and a large hole at the center. The Griewank function and Levy function have many widespread local minima, which are regularly distributed.

Table 2. The optimization results of the optimized value of CAL-SAPSO, DE, GA, and pre-SADE.

Problem	d	CAL-SAPSO	DE	GA	Pre-SADE
Ackley	10	19.426 ± 0.896	18.247 ± 0.448	18.454 ± 0.551	6.537 ± 0.392
Ackley	15	19.930 ± 0.580	18.923 ± 0.491	19.315 ± 0.815	9.374 ± 0.295
Ackley	20	19.454 ± 0.588	18.674 ± 0.715	18.471 ± 0.520	9.349 ± 0.489
Griewank	10	1.027 ± 0.299	56.924 ± 3.416	61.901 ± 2.825	2.233 ± 0.221
Griewank	15	0.771 ± 0.217	86.339 ± 6.432	83.464 ± 6.184	8.229 ± 0.513
Griewank	20	0.975 ± 0.164	192.354 ± 16.720	127.471 ± 12.472	9.575 ± 0.982
Levy	10	0.091 ± 0.003	17.035 ± 1.057	15.677 ± 1.144	0.007 ± 0.003
Levy	15	0.814 ± 0.074	46.932 ± 1.976	37.126 ± 1.725	0.144 ± 0.014
Levy	20	0.868 ± 0.009	58.179 ± 1.787	51.347 ± 1.830	0.790 ± 0.068

shows the results of the tests. Each test runs three times to involve the randomness in the testing processes. The mean values and standard deviations of objective functions are shown in Table 2. With the limited steps of exact function evaluations, pre-SADE obtains the best performance in the Ackley function, where other methods perform poorly. The optimization history of Ackley is shown in Figure 3. Pre-SADE shows significant improvements compared to the other three methods. In the Griewank problem, CAL-SAPSO obtains the best results, but acceptable solutions are also obtained by pre-SADE. Moreover, as shown in Figure 4, CAL-SAPSO gets the best result through a steep fall in every converge history. The steep fall may accelerate the converge process, but it also makes the algorithm fall more easily into the local optimum. Once the surrogate model fails to cover the global optimum solution, such as in Ackley, CAL-SAPSO may be stuck in the local optimum. In contrast, pre-SADE shows a continuous downward trend in every optimization history. Figure 5 shows the optimization history of the Levy function, and pre-SADE obtains the best performance. Compared with DE and GA, pre-SADE also acquires considerable accelerations in the converging process. These benchmark results show that pre-SADE can find competitive solutions in the search space within 20$d$ steps. Furthermore, pre-SADE can parallelly calculate the FEs in each generation, which saves considerable calculating time compared with CAL-SAPSO.

The Ackley function is a multi-modal optimization problem with a narrow and deep valley near the global optimum. The surrogate models have difficulties in getting enough information near the peak within limited exact samples. Thus, the total substitution algorithm may perform poorly in this type of problem. DE and GA may perform well in the Ackley problem in more generations, because each FE is exactly calculated, including the narrow deep valley. To avoid losing important information, surrogate models are only used to screen the FEs in pre-SADE, rather than directly calculate the value. With the same number of exact FEs, pre-SADE performs more generations than the DE and GA to get better solutions. In the aerodynamic optimization problems of compressors, the global optimum areas are usually very narrow due to the complicated and sensitive flow conditions. Thus, pre-SADE is used as the optimization method in the fast aerodynamic platform of axial compressors.

3. Rotor 37 Optimization

3.1. Case Studied

NASA Rotor 37 [35], a transonic compressor rotor developed by NASA in the 1970s, which has comprehensive experiment data, is widely used in aerodynamic optimization research [36,37]. Figure 6 shows the representation of Rotor 37. As a single row compressor with a high pressure ratio, Rotor 37 has typical transonic flow characteristics, such as strong shock wave boundary layer interactions, end wall separation, and tip-leakage vortex. These complex flow conditions bring challenges in aerodynamic optimization. To describe the geometry and implement the deformation in detail, the parameterization methods, such as FFD were introduced in Rotor 37 aerodynamic optimization [38]. This paper uses the RBF interpolation method to parameterize Rotor 37.

3.2. Optimization Settings

The aerodynamic optimization of Rotor 37 in this research mainly consists of three key approaches: (1) the RBF interpolation method provides the deformation of geometry, (2) pre-SADE is used as the optimizer, and (3) the 3-D computational fluid dynamic predicting program calculates the function evaluations of the samples.

As shown in Figure 7, 18 control points (CPs) are selected and distributed in the tip, mid, and hub surface of the rotor blade. Each surface contains six CPs. Two CPs are located at the leading and trailing edges. Four CPs are located at the position of 0.3 and 0.6 in the chord length of the suction surface and the pressure surface, respectively. To keep the smoothness of the geometry, the distance of the CP to the blade is set to 0.4 cm, and the affect radius of the RBF interpolation method is set to 1.0 mm. The type of RBF kernel function is selected as the Gaussian spline, $e^{-r^2}$in this method. As the blade heights are usually coupled with the flow path, the z direction deformations are not introduced in the optimization. In multi-stage compressors, the locations of the stages and the distances between the rows are generally determined in the framework. To make this optimization more applicable, the CPs of the hub surface are fixed to guarantee the axial location of the row. In this case, the movement of CPs in the mid and tip surfaces are selected as the optimization variables, and each point contains two direction movements (X and Y direction). Thus, there are 24 independent variables in this optimization problem. The ranges of the movements are set as [−1 cm, 1 cm].

NUMECA, as a widely used commercial turbomachinery CFD software, is introduced as the 3-D viscous steady predicting program in this optimization. The Spalart–Allmaras model is used as the turbulence model. The time integration is the fourth-order Runge–Kutta scheme. The inlet boundary conditions are provided according to the experimental stations. The inlet total pressure is 101,325 Pa and the inlet total temperature is 288.15 K. The back static pressure is selected as the outlet boundary condition. The Autogrid5 module is used for the mesh generation of each sample. The rotor tip clearance is set to 0.356 mm. After the verification of mesh independence, the radial station number of the computing meshes is set to 57, which contains 17 stations in the tip clearance. The cell width of the mesh near wall is set to 0.001 mm to guarantee the near-wall $y^{+}$< 5 [39]. The mesh topology of the rotor is O4H. The 11 parallel processes are adapted in each 3-D calculation. Figure 8a shows the illustration of 3-D mesh.

Figure 8b,c shows the comparisons of the simulation and experiment results about Rotor 37. At the design point, the adiabatic efficiency and total pressure ratio of the initial rotor are 86.17% and 2.008, respectively. The efficiency is, on average, 2% lower than the experimental data, and the total pressure ratio agreement is reasonable, which is consistent with the previous study [40]. Generally, considering the overall trends of the efficiency and pressure ratio, the simulation results have good agreement with the experimental data.

The basic parameters of pre-SADE depend on the dimension of optimization variables, the calculation cost of a single sample, as well as the number of parallel nodes. The F factor, which describes the amplification of the differential variation in DE, is set to 0.5. The recombination probability of DE is set to 0.7. The max number of generations is set to 10, and the population of each loop is set to 100. To balance the accuracy and the training cost, the max volume of the training database of the integrated surrogate model is set to 300. The first generation of pre-SADE needs to be exactly evaluated, and the exact evaluation ratio ($e_r)$ in the next generations is set to 25%. Thus, 325 samples need to be exactly evaluated in the optimization process. Parallel techniques are introduced to accelerate the optimization. Five parallel computation nodes, which can calculate five samples at the same time, are used in the prediction of exact FEs. To further reduce the calculation time, the CFD meshes of the samples are generated at the same time (25 parallel nodes) before the calculation process. The optimization objective is set to maximize the adiabatic efficiency near the design point, which is calculated under the outlet back pressure of 130,000 Pa. The constraints of the optimization are set as follows:

$$\{\begin{array}{c}\left|\frac{mas{s}_{opt}-mas{s}_{ini}}{mas{s}_{ini}}\right|\le 1.0\%\\ \left|\frac{TP{R}_{opt}-TP{R}_{ini}}{TP{R}_{ini}}\right|\le 1.0\%\end{array}$$

(11)

where $mas{s}_{opt}$ and $mas{s}_{ini}$ denote the flow rate of the optimized and initial rotors, respectively. $TP{R}_{opt}$ and $TP{R}_{ini}$, respectively, refer to the total pressure ratio of the optimized and initial rotors.

Thus, the optimization statements are given in the following:

• Variables: 24. (The movements of the X, Y directions of the selected 12 control points).
• Ranges: −1 cm to 1 cm.
• Objective: Maximize the adiabatic efficiency at the design point.
• Constrains: The changes of flow rate and total pressure ratio are lower than 1.0%.
• Predicting software: NUMECA (Spalart–Allmaras turbulence model).

4. Results and Discussion

4.1. Optimization Results

In this paper, AMD-3990X (64 core, 2.9 GHz) is used to execute the optimization. In the optimization process, 325 samples are exactly calculated, and the total run time is 2.007 h. CAL-SAPSO and DE are selected to implement the same optimization problem as comparisons. The optimization results are shown in

Table 3. The optimization results of pre-SADE, CAL-SAPSO, and DE.

Method	Exact Calculated Steps	Optimization Results	Increment	Time (h)	Parallel Nodes
Pre-SADE	325	87.64%	1.47%	2.0	5
CAL-SAPSO	325	87.13%	0.96%	5.1	5 (initialization) + 1
DE	325	87.07%	0.90%	1.9	5
Pre-SADE	825	87.94%	1.77%	6.0	5
CAL-SAPSO	825	87.42%	1.25%	25.6	5 (initialization) + 1
DE	825	87.54%	1.37%	5.2	5

. In the limited number of steps, the adiabatic efficiency at the design point optimized by pre-SADE is 1.47% higher than the initial rotor. DE and CAL-SAPSO get a 0.90% and 0.96% improvement, respectively. Figure 9a shows the predicting accuracy of the surrogate model (SM). The predicting accuracy is defined as the ratio of the number of samples, which are both in the SM-screened assembles and the true sorted assembles, to the total true calculated number. The predicting accuracy of SM is tested based on the run of DE, which exactly calculates all the samples in each generation. As the database updates, the predicting accuracy keeps higher than 60% in most of the generations, meaning that the SM can effectively screen the samples before exact calculations. Moreover, few inappropriately screened samples can provide the exploration to the algorithm. Figure 9b shows the optimization history of the three optimizers in 325 steps. To further tap the potential of pre-SADE, the number of generations is increased to 30. The max step number is 825, subsequently. DE and CAL-SAPSO are also executed as comparisons. The results are shown in Table 3 and the converge history is shown in Figure 9c. At the design point, pre-SADE, DE, and CAL-SAPSO get a 1.77%, 1.25%, and 1.37% increase in adiabatic efficiency, respectively. In the 325 step runs, pre-SADE and DE finished the optimization in 2.0 h benefiting from the parallel improvement, compared to 5.1 h of CAL-SAPSO, which also used parallel systems in the initialization process. In the 825 step runs, pre-SADE and DE finished in 6.0 h compared to 25.6 h of CAL-SAPSO.

Pre-SADE shows significant improvements compared to DE and CAL-SAPSO, which verifies that the fast optimization platform based on the RBF interpolation method and pre-SADE can implement the compressor aerodynamic optimization of Rotor 37 in limited steps with competitive solutions. Thus, the fast optimization platform can finish the aerodynamic optimization process of Rotor 37 in hours, which makes the optimization or design process of the axial compressor more flexible.

4.2. The Analysis of 3-D Deformations and Flow Conditions

Figure 10 shows the comparisons of the optimized and the initial blade. The 3-D view and the 2-D section view of a 90% and 50% span are provided, respectively. Through the 2-D blade section view, the maximum thickness of each section is not changed that can keep the structural characteristic of the rotor. The main deformation occurs at the shape of the mean camber line. In the 50% span section, the optimized blade presents a nearly S-shape along the camber line. An S-shape blade can generate pre-compression waves to weaken the shock wave, and subsequently reduce the flow loss. In the 90% span section, the suction surface is slightly concave near the leading edge. The detailed analyses of the performances and flow conditions are shown in the following section.

The performances of the initial and optimized rotors at the design point are shown in

Table 4. The overall performance of optimized and initial rotors.

Run	Adiabatic Efficiency	Surge Margin	Total Pressure Ratio	Flow Rate (kg/s)
Initial	86.17%	13.2%	2.008	20.99
Optimized	87.64%	14.2%	2.024	21.15
Increment	+1.47%	1.0%	+0.79%	+0.76%

. The surge margin is defined as follows:

$$SM=\left(\frac{TP{R}_c/mas{s}_c}{TP{R}_d/mas{s}_d}-1\right)\times 100\%$$

(12)

where $TP{R}_c$ and $mas{s}_c$ denote the total pressure ratio and flow rate of the choke point, respectively. $TP{R}_d$ and $mas{s}_d$, respectively, refer to the total pressure ratio and flow rate at the design point.

As the optimization objective, the adiabatic efficiency is increased by 1.47%. The surge margin is increased by 1.0%, which means that the optimization process did not affect the aerodynamic stability of the rotor. The optimized total pressure ratio and flow rate of the design point are, respectively, increased by 0.79% and 0.76% compared to the initial performance, which are consistent with the constraints of the optimization settings.

Figure 11 shows the comparisons of the performance characteristic between optimized and initial rotors at the design point. At the range of flow rate, the adiabatic efficiency gains significant improvement through optimization, and the total pressure ratio is consistent with initial performance after the optimization.

In the actual airflow, more flow loss usually leads to higher entropy. The contours of entropy distributions in the blade channel are shown in Figure 12. Two channel sections are selected to show the entropy comparison in detail. In the upstream of the shock waves, channel Section 1 is placed near the 0.3 camber length. The detail entropy distribution shows that the optimized entropy is slightly decreased at the tip clearance. Channel Section 2, located near the 0.6 camber length, lies in the downstream of the shock waves. After optimization, the entropy near the suction surface is significantly reduced, which indicates that the initial rotor suffers more flow loss than the optimized geometry near the suction surface in the downstream of the shockwave.

To further analyze the reasons for the entropy decrease in the suction surface, the detailed flow conditions are shown in Figure 13. The limiting streamline describes the 3-D flow near the blade. The contours of the relative Mach number (RMa) are shown through the blade to blade (B2B) view, where the shape of the shock waves and the location of the flow separation can be directly observed.

The strong shock waves in the blade flow channels cause the large separation and the complex shock wave boundary layer interaction, which reflects in the separation line at the limiting streamline on the suction surface. Compared to the initial rotor, the separation line of the optimized rotor obviously moves downstream at the tip and mid-height of the blade. In the 90% span and 50% span of the B2B view, the distributions of the RMa also indicate that the flow separation near the trailing edge is notably reduced through the optimization.

The limiting streamline shows that the reverse flow near the 70% span is reduced, and the flow reattachment occurs near the trailing edge. In the B2B view, the airflow also shows obvious reattachment after the slight separation at the downstream of the shock wave. The reattachment can reduce the flow separation near the tip region. Near the 50% span, there exists a strong reverse flow in the initial rotor, which leads to the large flow separation and flow blockage. After optimization, the reverse flow area is released, and the jammed flow is dispersed to the mid-height of the blade.

At the 50% span section of the B2B view, the max Mach number of the initial blade is 1.54, and the Mach number is reduced to 1.43 after optimization. The shock waves of the initial flow are nearly normal. The strong normal shock waves bring large flow loss in the blade channels. After optimization, the normal shock waves become oblique shock waves, which can reduce the flow loss of the shock waves in the blade channels.

In summary, the decrease in the flow separation, the reduction in the reverse flow area, and the weakness of the shock waves all result in the increase in the adiabatic efficiency in the overall performance.

5. Conclusions

This paper first combines the direct RBF interpolation with the pre-screen data-driven optimization algorithm to establish a fast aerodynamic optimization platform of axial compressors. In the platform, the RBF interpolation method provides the parameterization and the deformation techniques for the blade to reduce optimization variables. Pre-SADE, which is proposed and tested in this paper, can accomplish the optimization in a limited budget of expensive function evaluations. Rotor 37, which is very sensitive to the deformations of the blade due to the complex flow conditions, is selected and optimized as the verification sample of the platform. The main findings are listed below:

• Different from most surrogate model-assisted algorithms, pre-SADE screens the samples through the integrated surrogate model to save total optimization time. By avoiding directly estimating the samples, pre-SADE can reduce the dependence on the accuracy of the surrogate model. Pre-SADE, CAL-SAPSO, DE, and GA are used to run the benchmark tests. Under the limitation of max steps, pre-SADE shows significant improvements in the Ackley function and gets competitive solutions in the Griewank function in the benchmark tests of 10, 15, and 20 dimensions. The results of the benchmark tests illustrate the effectiveness of pre-SADE in the optimization with a limited budget of exact function evaluations;
• The aerodynamic optimization of Rotor 37, based on the presented platform, modified the blade shape mainly in the 50% span to the tip of the blade. The optimization leads to the decrease in large flow separation at the suction surface, the reduction in the reverse flow area, as well as the weakening in normal shock waves. The improvements in the flow conditions contribute to the promotion of the overall performance at the design point. The adiabatic efficiency, total pressure ratio, mass flow rate, and surge margin of the final optimized compressor are 87.64%, 2.024, 20.40 kg/s, and 14.0%. Compared to the initial compressor, the adiabatic efficiency and the surge margin have, respectively, increased by 1.47% and 1.0% under the optimization constraints;
• The limitation to the number of exact function evaluations is set to 325. Pre-SADE gets a 0.51% and 0.57% higher adiabatic efficiency than CAL-SAPSO and DE, respectively. Pre-SADE can parallel computing though an original integrated parallel management program. With five parallel nodes, pre-SADE and DE finish the optimization process in 2 h, compared to the 5.1 h of CAL-SAPSO. The validation optimization results show that this platform has the ability to quickly implement the aerodynamic optimization of axial compressors.
• Through the RBF interpolation method, the global deformations (bow, sweep, and twist) and local deformations (S-shape) of the Rotor 37 blade can be implemented by the movements of several control points (18 in this paper). The optimization results show that direct RBF interpolation has good performance as the parameterization method in the aerodynamic optimization of axial compressors.