#### Abstract

Multistage axial compressor is the key component of aeroengine and gas turbine to realize energy conversion. In order to avoid the “curse of dimensionality” problem in the global optimization process of AL-31F four-stage low-pressure compressor under multiple working conditions, an optimization method based on phased parameterization strategy is proposed. The method uses the idea of “exploration before exploitation” for reference and divides the optimization process into two phases. In the first phase, the traditional parametric modification method based on stacking line is adopted; in the second phase, the full-blade surface parametric modification method with significant low-dimensional characteristics is adopted. Based on the improved artificial bee colony algorithm, a multitask concurrent optimization system is built on the supercomputing platform, and the engineering optimization solution is obtained within 91 hours. The optimization results are as follows: under the condition of meeting the constraints, the adiabatic efficiency is increased by 0.3% and the surge margin is 4.0% at the design speed; the adiabatic efficiency is increased by 0.8% and the surge margin is 2.3% at the off-design speed. These results verify the usefulness and reliability of the optimization method in the field of aerodynamic optimization of a multistage axial flow compressor.

#### 1. Introduction

A multistage axial flow compressor is the core component of aeroengines and gas turbines and plays an important role in military and civil fields. To better realize the “three high” performance (high efficiency, high pressure ratio, and high surge margin) in the design of the compressors, the optimization design method is usually carried out after expert design. Although computing power has greatly advanced in recent decades, the aerodynamic optimization design of multistage axial flow compressors still requires excessive calculation costs. The basic reason is that the high-dimensional, expensive, , and black-box (HEB) characteristics of the aerodynamic optimization of multistage axial flow compressors easily suffer from the “curse of dimensionality” [1], and it is difficult to obtain an optimized solution in the limited time range of engineering.

In view of the above HEB problems, researchers have proposed many solutions, which can be divided into three categories: the surrogate model method [2–5], the adjoint method [6–8], and the dimensionality reduction method [9–11]. The surrogate model method uses a simple approximate function to replace the original complex PDE numerical solution, which greatly reduces the optimization time. The drawback of this method is that the number of samples needed to train the surrogate model increases sharply with the increase of optimization variables. Therefore, the surrogate model method is not suitable for the optimization design of a multistage axial flow compressor with many design variables. The adjoint optimization method is first applied in the field of external aerodynamic optimization, and then it is introduced into the field of internal flow by Jemenson [12, 13]. The characteristic of this method is that the number of iterations required for optimization is independent of the size of optimization control parameters, which is a considerable advantage for HEB problems. However, this method depends on the initial design and is a local optimization algorithm. To find the global optimized solution, many different initial blades must be designed first, which increases the complexity and difficulty of the optimization and makes it more difficult to deal with multiobjective and multiconstraint problems. The dimensionality reduction method is an important method used to reduce the dimensionality of an optimization problem for a multistage axial flow compressor from different perspectives to reduce the design space scale. It can be divided into parameterization mapping dimensionality reduction [14, 15], spatial decoupling [12], stage-by-stage optimization [16], and variable selection [17, 18]. References [19, 20] note that the parameterization mapping dimensionality reduction method is one of the most promising solutions to the HEB problem, and its core concept is to reduce design variables on the premise of preserving the optimal solution of high-dimensional space as much as possible. The surface parameterization method is an effective parameterized dimensionality reduction method that is first applied to the optimization design of an external flow airfoil [21–24] and later introduced to internal flow [25–27]. Because the essence of this method is to constrain the radial parameters so that the sections are no longer independent changes, thus greatly reducing the design variables, the radial smoothness of the profile is realized. However, the problem of this method is that it only exerts an effect on half of the blade surface. At the same time, the leading edge blade angle, which has an important influence on the development of air flow, cannot be adjusted, which limits the space for performance improvements. In this paper, a full-blade surface parameterization method is used in the optimization process. Compared with the semiblade surface parameterization [25–27], this method has excellent modification flexibility. It can keep the low-dimensional characteristics while taking into account the changes of suction, pressure surfaces, and leading edge metal angle.

The advantage of the automatic optimization method is that it is general and suitable for all kinds of optimization problems; it also ignores the particularity of different optimization tasks, so it does not make full use of the specific experience and knowledge of experts for design tasks. This “man-in-the-loop” concept is proposed to solve this problem [28, 29]. The core concept is to introduce the specific experience of the designer for the design object into the optimization process to avoid excessive invalid exploration of the design space during the optimization process. In the 1990s, when the optimization method was in its initial phases, “man-in-the-loop” was ignored with the gradual development of the optimization algorithm. With the increasing scale of optimization problems and the prominence of HEB problems, the concept has gained the attention of researchers at home and abroad and has begun to play a role in the field of aerospace engineering. The “man-in-the-loop” concept improves the efficiency of optimization design in two ways: (1) The first is introducing the long-term experience accumulated by experts into the optimization process in the form of an objective function or constraint. For the transonic wing design of narrow body airliners, Zhang [30] mathematically described the weak shockwave shape as a constraint condition and introduced it into the optimization process according to the design experience. During the optimization process of the evolutionary algorithm, individuals with poor matching between the pressure distribution and the expected shape can be eliminated quickly, and the optimization can be conducted more efficiently for the final direction of synthesis. (2) The second is optimization process monitoring and intervention. In evolutionary algorithm optimization, population regulation [31, 32], adaptive optimization algorithm development [33–35], and global and local algorithm combination [36–38] are carried out. Han et al. [39] proposed a hybrid design method combining inverse design and optimization design with a highly efficient global surrogate optimization framework that has been developed over many years, which allows designers to directly “intervene” in the numerical optimization process according to their design experience and their understanding of the flow mechanism; they designed three low-drag and high-lift double-blunt nose wings for different sections of a blade root, which is a good realization of their design goal.

Applying the “man-in-the-loop” concept from the perspective of the parameterization method is an important way to reduce the design space and the optimization time. Trepanier et al. [40] used a parameterization method that combined global scope and local scope to optimize compressor blades and introduced the design experience in the process, which achieved their optimization goal. Guangwen et al. [41] used the recursive algorithm of the Bezier curve to optimize the compressor blade profile by a multilevel parameterization method, which accelerated the convergence process. Su et al. [42] used the three-dimensional CST parameterization method to carry out two-level parameterization optimization of aircraft geometry, combining global parameterization with local parameterization, which significantly reduced the calculation consumption while ensuring the space search ability. Although the multilevel parameterization method can reduce the optimization time to a certain extent, it does not address the dimensionality reduction problem of multistage axial flow compressor aerodynamic optimization due to the lack of complete utilization of expert experience in the optimization loop.

To solve the HEB problem in high-fidelity aerodynamic optimization of a multistage axial flow compressor, this paper proposes an optimization strategy based on phased parameterization under the guidance of the “man-in-the-loop” concept: in the first stage, the “mode search” of the stacked line modification is adopted and combined with the blade-row-linkage strategy to address the “repeated stage flow” phenomenon of a multistage axial flow compressor; the design space is explored by a large step under the optimal Latin square test method. In the second stage, based on the analysis of the optimal solution in the first stage, a full-blade surface parameterization method is adopted for local fine optimization.

The structure of this paper is as follows: Section 2 introduces the specific connotation of the phased parameterization optimization strategy; Section 3 introduces the different parameterization methods of the two phases; Section 4 introduces the optimization algorithm needed for the optimization task and the multitask concurrent optimization system built on the supercomputing platform; Section 5 adopts the proposed method to perform the optimization for a four-stage low-pressure compressor under multiple working conditions; Section 6 presents some conclusions.

#### 2. Phased Parameterization Optimization Strategy for Multistage Axial Compressors

The purpose of engineering optimization is to find the suboptimal solution in a design space in a limited time range. Therefore, the original intention of using an optimization strategy in engineering optimization tasks is to sacrifice certain global optimization accuracies to improve the speed of optimization.

“Balance of global exploration and local exploitation” is an important optimization concept that has been widely studied and discussed in the fields of evolutionary algorithms and surrogate models, but insufficient attention has been paid to the importance of this concept in parameterization. In fact, the parameterization method directly determines the structure of a design space, while the optimization algorithm only looks for the optimized solution in the design space with a certain structure. Therefore, for the whole optimization task, one must first determine the structure of design space, and a reasonable design space structure also needs to apply the concept of “balancing global exploration and local exploitation” to achieve a better outcome.

The shape of stacking line determines the three-dimensional configuration of the blade, reflects the relative spatial position of each two-dimensional radial profile, and has an important impact on the aerodynamic performance of the blade. It is also the last design step in the traditional blade modeling process. Besides, the change of blade body determines the shape of each two-dimensional blade profile. Theoretically, the effects of three-dimensional configuration and two-dimensional profile geometry on blade performance are coupled and need to be considered as a whole. However, in order to further reduce the dimension, it is acceptable to sacrifice a certain design space and decouple the modification of stack line and blade body. The proposed optimization strategy divides the optimization process into two phases: each phase adopts different parameterization methods, and the design optimization process introduces the designer’s experience to better achieve the balance between global exploration and local exploitation.

Because of the “repeated stage flow” phenomenon in a multistage axial flow compressor, that is, because of the effect of radial mixing and turbulent diffusion, the radial distribution of the outlet velocity of each stage is the same as that of the inlet velocity of this stage. To be precise, not all stages of the low-pressure compressor are repetitive stages. Since there are shock waves in the first stage of the low-pressure compressor and no shock waves in the later stages, the “repeated stage phenomenon” is more obvious in the second, third, and fourth stages. However, for further dimensionality reduction, this paper still links all rotor blades and stator blades separately. The “repeated stage phenomenon” corresponds to the “repeated stage design”; i.e., the same rotor and stator blade geometry can be used for different stages. For the design optimization performed in this paper, the blade geometry of different stages is modified in the same way. This paper proposes a blade-row-linkage strategy; that is, during the first optimization phase, the rotor blade and the stator blade are modified in the same way, respectively, and various possible optimization modes of the blade are explored to lay a good foundation for the second phase of local exploitation. The blade-row-linkage strategy can greatly reduce the dimension of multistage axial compressor optimization task and effectively promote the solution of the HEB problem.

The phased parameterization optimization strategy is explained as follows:(a)In the first phase, the parameterization method based on the stacking line modification and the blade-row-linkage strategy is used to perform the synchronous change of the sweep, bend, and stagger angle of all rotor and stator blades, respectively, to achieve a large range of blade modifications and a large step of space exploration. At the same time, the optimal Latin square test method is used to obtain the whole space information with as few samples as possible and explore various possible optimized blade modes. This phase can effectively reduce the optimization task dimension of the multistage axial compressor.(b)In the second phase, according to expert experience, the optimized blade of the first phase is selected, and the flow field is analyzed. In view of the part of the flow field that is relatively bad, the corresponding blade control variables are selected, and the second phase optimization is carried out by using the full-blade surface parameterization method and the improved artificial bee colony (IABC) algorithm to realize local exploitation under the existing blade mode.

The above two-phase parameterization method and optimization algorithm selection, the use of the blade-row-linkage strategy, and the analysis of the flow field optimization in the first phase have all introduced expert experience, which makes the “curse of dimensionality” problem of multistage axial flow compressor aerodynamic optimization well solved.

#### 3. Parameterization Methods

The optimization method proposed in this paper involves two different parameterization methods: one is based on stacking line, and the other is based on the full-blade surface parameterization.

##### 3.1. Parameterization Method Based on Stacking Line

The two-dimensional blade profile of each section obtained in the modeling process needs to be accumulated radially to generate three-dimensional blades. For a compressor, the center of gravity stacking method is generally adopted. By changing the shape of stacking line, the bending, sweeping, and stagger angle of the blade can be changed. As shown in Figure 1(a), the bent blade can be formed by moving the two-dimensional blade profile along the direction perpendicular to the chordwise; as shown in Figure 1(b), the swept blade can be formed by moving the two-dimensional blade along the chordwise direction; as shown in Figure 1(c), the change in the blade stagger angle can be realized by rotating the two-dimensional blade profile along the center of gravity.

**(a)**

**(b)**

**(c)**

The blade parameterization method based on the stacking line needs to select several control points of the shape of stacking line firstly and then change the shape of stacking line smoothly through these control points, so as to obtain the required corresponding moving distance or rotation angle of each radial section, so as to obtain a new blade geometry.

##### 3.2. Full-Blade Surface Parameterization Method

The traditional blade parameterization method is based on the traditional modeling method, which needs to change some two-dimensional radial section geometry. Generally speaking, more than 60 parameters are needed for single-row blade. In contrast, the surface parameterization method has a significant advantage of low dimension, while the newly developed full-blade surface parameterization method has good modification flexibility, while ensuring the low dimension characteristics, and can take into account the changes of suction surface, pressure surface, and leading and trailing edge blade angle. Figure 2(a) shows a schematic diagram of the full-blade surface parameterization principle. It can be seen that the method is realized by the following steps:(a)The original blade is opened radially along the leading edge to form a whole surface(b)The coordinates of each geometric point on the whole surface are parameterized by chord length, which is converted to the (*ξ*, *η*) coordinate system(c)The Bezier surface is generated(d)A corresponding value on Bezier surface is found for each data point of blade surface in the (*ξ*, *η*) coordinate system.(e)This value is superimposed on the circumferential direction of the corresponding points on the original blade surface to form a new blade.

**(a)**

**(b)**

The chord length parameterization formula in the above step (b) is (1) and (2), and the Bezier surface construction formula in step (c) is (3)–(5):where ; refers to the number of data points of a section; refers to the total number of sections; refers to the sum of the chord lengths of each segment of the *j*th section in the direction; refers to the sum of the chord lengths of the *i*th section line in the *ξ* direction; is the *m*th chord length of the *j*th section line in the *η* direction; and is the *n*th chord length of the *i*th section line in the *ξ* direction. refers to the change value of each point on the perturbed surface; is (*m* + 1) × (*n* + 1) control points of Bezier surface; are the Bernstein basis function determined by equation (4); and *u* are two independent variables in the computational domain of the perturbed surface, and its variation range is [0, 1]; and is the number of combinations determined by equation (5).

Figure 2(b) shows the distribution of control points in the full-blade surface parameterization method. The red points are the free active points; i.e., each red point can be changed independently. The yellow points are the limited active points; i.e., four rows of yellow points are changed synchronously to ensure the geometric continuity of the leading edge and increase the degree of freedom of geometric change of leading edge. It can also be seen from the figure that the distribution of points in the middle region and trailing edge area is sparse, while the distribution near the leading edge is dense because the blade performance is sensitive to the shape of the leading edge. The distribution law of the *η* direction can be adjusted according to the actual demand of engineering design. Generally, it is more suitable to arrange 7 points in the *ξ* direction and 4 points in the *η* direction because under the action mechanism of the Bezier surface, more points mean higher dimensions and increase the probability of mutual interference between the action ranges of different control points. This parameterization method can reduce the control parameters of a single row of blades to 16, which can well realize the dimensionality reduction of optimization tasks.

Compared with the traditional blade parameterization method, the full-blade surface parameterization has obvious advantages: (1) the blade surface modification is realized by surface superposition, and there is no need to carry out surface fitting through back calculation of control points, so there is no problem of fitting accuracy, and it has good construction convenience; (2) due to the high-order continuity of Bezier surface, the smoothness of the parameterized blade surface is not lower than that of the original blade surface; (3) because the essence of the surface parameterization method is to constrain the control points in both *ξ* and *η* directions, the method has significant low-dimensional characteristics.

#### 4. Optimization Algorithm and Optimization Platform

Generally, the design space of a multistage axial flow compressor has multiple peak characteristics. To avoid falling into a local extremum and find a better global solution, ABC algorithm [43] is used to perform the design optimization task. The principle of IABC method is shown in Figure 3. It essentially uses initialization, exploration of employed bees, onlooker bees, and detector bees to recognize global and local information. As an improvement of the standard artificial bee colony algorithm (for a detailed introduction of the ABC algorithm, refer to [44]), it is mainly implemented from three aspects: (1) an initialization method, (2) a food source exploration mechanism of employed bees, and (3) a food source exploitation mechanism of onlooker bees. The detailed explanation is as follows:(a)Initialization mode: by using the optimal Latin square test method instead of the random initialization method of the ABC algorithm, more comprehensive design space information can be obtained with fewer samples.(b)Food source exploration mechanism of employed bees: in IABC algorithm, the food source adopts formula (7) to replace the exploration mechanism in formula (6) in the ABC algorithm, making full use of the global and local optimal food source information.(c)Food source exploitation mechanism of onlooker bees: the exploitation mechanism of onlooker bees in IABC algorithm uses formula (8) instead of formula (6), making full use of the information of the local neighborhood optimal food source, realizing more faster convergence speed and better global optimization accuracy.where refers to the *j*th component of the position of the *i*th bee’s newly explored food source, ; , , and are all random numbers between ; and is the most suitable location of food source in the overall information. refers to the best food source location in the neighborhood, and refers to the location of new onlooker bees to explore.

The neighborhood is determined by the Chebyshev distance, which is defined as where is the Chebyshev distance between and .

Compared with GA and ABC algorithms, IABC algorithm has better global optimization accuracy and convergence speed, and more details can be found in [44].

To further reduce the optimization time, a supercomputing platform with strong computing power is considered. Because the evolutionary algorithm does not need the gradient information of the objective function during the optimization process, it can achieve multitask concurrency; therefore, this paper builds a multitask concurrency optimization system of a multistage axial flow compressor based on ABC algorithm and supercomputing platform. The optimization system can reduce the optimization time to approximately one-tenth of the original.

Figure 4 shows the schematic diagram of the optimization system, which consists of a user’s computer, logged node, and calculating node. The process of initialization and food source exploration of employed bees and onlooker bees in IABC algorithm is the colony behavior of bees. When these colony behaviors occur, many engineering folders are generated by scripts on the user’s computer. Different geometric files are generated in each folder by IABC algorithm. They are transmitted to the logged node through an online transmission function, recorded, and then distributed to the calculating node for calculation (grid generation, flow field calculation, and postprocessing). After the calculating node obtains the results, it feeds back to the logged node, makes records, feeds back to the user’s computer through the synchronization function, and updates the geometric files in different folders. According to this iteration, after the feedback of *N* calculated results, it generates *N* new task files until it jumps out of the cycle and obtains the optimization results.

#### 5. Aerodynamic Optimization of AL-31F Four-Stage Low-Pressure Compressor under Multiple Working Conditions

The proposed optimization method is used to perform the aerodynamic optimization for AL-31F four-stage low-pressure compressor under multiple working conditions. AL-31F is a fourth-generation aeroengine produced in Russia. Its compression parts are composed of a four-stage low-pressure compressor and a nine-stage high-pressure compressor. Among them, the four-stage low-pressure compressor has 10 rows of blades, including an inlet guide vane and an outlet cascade, as shown in Figure 5. The design pressure ratio of the four-stage low-pressure compressor is 3.6, the designed adiabatic efficiency is 81%, and the surge margin is 14.9%.

##### 5.1. Numerical Methods

Checking the numerical method is the necessary process before optimization. This paper uses the Fine_Turbo module of the commercial software NUMECA V11.0 to calculate the steady flow field. The grid generation adopts the Autogrid5 module to build the 4HO grid and adopts a gapless model for all blade rows. The fourth-order Runge–Kutta time discrete scheme and the central difference scheme are used to discretize the time and space; the auxiliary artificial viscosity technique is used to eliminate the calculation oscillation; and the local time step, multigrid, and implicit residual smoothing techniques are used to accelerate the convergence. The compressor performance calculated in this paper is typical of end-wall constrained flow, and there is no significant flow scale variation in this case. The S-A model is a simple one-equation model for this end-wall constrained case and can capture the turbulence in this case. Meanwhile, the S-A model is less sensitive to numerical errors caused by mesh roughness and has been widely used in turbomachinery. Therefore, the S-A model is used as the turbulence model in this paper.

The boundary conditions are set as follows: the designed rotation speed of the four-stage low-pressure compressor is 10,200 r/min. The total inlet temperature and total pressure are set as 293 K and 101325 Pa, respectively, and the inlet air flow is along the axial direction. No slip model is adopted for the wall boundary. Given the back pressure at the outlet boundary, the whole compressor characteristic line is calculated under different back pressure conditions. From the design point, the blocking point can be approached gradually by reducing the back pressure; from the design point, the surge point can be approached gradually by increasing the back pressure.

The flow field calculation method is checked by comparing the experimental data [47] with the simulation data of a four-stage low-pressure compressor. As seen from Figure 6, the black and red lines represent the experimental and simulated compressor performance data, respectively. It is clear that the simulated performance has no flow margin. The reason why there is no flow margin in the numerical simulation results is the phenomenon where various errors (errors in the turbulence model, numerical errors, etc.) are continuously amplified and accumulated in the numerical calculation. Thus, the characteristic line cannot really be calculated to the surge point, but only at the near-stall point. This calculated surge margin is certainly not accurate enough. The difference between experiment and numerical simulation under the small flow rate conditions should be attributed to the steady calculation adopted by this research. Besides, the flow field calculation only involves a single blade-passage, not the full blade-passage, may be another reason. However, since the goal of this paper is to optimize the blade, it is only necessary to ensure the consistency of the numerical calculation method used before and after optimization. It can be considered that making the flow field of the near-stall point smoother can expand the surge margin. The maximum pressure ratio calculated by numerical calculation is slightly larger than the experimental data by 3.2%, and the maximum efficiency is 1.5% larger than the experimental value. At the design speed, the pressure ratio range and efficiency range of simulation value are slightly increased compared with the experimental value, but the flow margin of simulation is 4.7% smaller than the experimental value. In general, the trend of calculated value and experimental value under different back pressures is the same, and the coincidence is good. It can be considered that the numerical simulation method is effective and reliable.

**(a)**

**(b)**

In addition to checking the flow field calculation method, grid independence verification is also needed. Figures 7(a) and 7(b) show the 3D blade grid and its details and the B2B mesh of the hub of R1, and Figure 7(c) shows the grid independence verification. The numbers of single-row blade grids of the four sets of grids are 250,000, 500,000, 930,000, and 1,250,000. The near-wall condition index *y*^{+} is a very important index for ensuring the mesh quality, and the value of *y*^{+} for the four sets of meshes is kept less than 3, as shown in Figure 7(c). Figure 7(d) shows that when the numbers of grids in a single row of blades is close to 1 million, the increase in the number of grids has a negligible impact on the performance. It has been noted in [45, 46] that coarse mesh optimization does not affect the trend of blade performance changing with geometry, so the optimization time can be further reduced by using coarse grids for optimization and fine grids for verification. In this paper, the first set of grids (250,000 in a row) is used for optimization, and the third grid (930,000 in a row) is used for numerical simulation of the flow field and aerodynamic performance. The use of this strategy can reduce the optimization time by half.

**(a)**

**(b)**

**(c)**

**(d)**

##### 5.2. Application of Optimization Method

The optimization method explained in Figure 8 is used to improve the aerodynamic performance of AL-31F four-stage low-pressure compressor. The original blade geometry is the input of the optimization cycle. In the first phase, the blade-row-linkage strategy is adopted and combined with the parameterization method based on the stacking line, and the optimal Latin square test method is adopted to explore the design space. Then, according to the designer’s experience, the optimized blade in the first phase is selected, and the flow field is analyzed to determine the optimization variables in the second phase. In the second phase, the full-blade surface parameterization method and ABC algorithm are used to find the optimized solution in the second phase. Finally, according to the designer’s experience, the final optimized blade is determined.

###### 5.2.1. Optimization in the First Phase

*(1) Optimization Objectives and Constraints*. This optimization is a multiple condition optimization conducted at the design speed (10,200 r/min) and at a reduced speed of 0.8 × the design speed (8160 r/min). The optimization objective is to improve the efficiency of the design point, working point, and surge margin at the two speeds. According to the author’s optimization experience, selecting the near-stall condition as the optimization condition (the design speed 340,000 Pa back pressure, the off-design speed 240,000 Pa back pressure) can achieve an optimized solution that can not only improve the efficiency at design point or working point but also improve the surge margin when only one flow field calculation is carried out in each optimization iteration. The objective function setting of the actual optimization process is shown as follows:

The constraints are shown as follows:where and are the isentropic efficiencies of the optimization operating point at the design speed and off-design speed. and represent the calculated efficiencies of the operating point at the design speed and off-design speed. and are their respective weighting coefficients. This optimization hopes to improve the performance at the design speed and off-design speed in a balanced way, so and are set to 1, and are the total pressure ratio in the optimization process and the referenced total pressure ratio at the design speed, and mean the total pressure ratio at the optimization operating point and off-design speed, and are the new flow rate and referenced flow at the design speed, and and are the flow rate at the off-design speed. “minimum” is a minimum value given by human, which aims to eliminate the solution that does not meet the constraints in the optimization process of IABC algorithm. is the design variable; is the lower limit and is the upper limit of the optimization variable.

The formulas for , , and SM (surge margin) are as follows:where are the total pressure ratios at the outlet and inlet, respectively; is the specific heat ratio; are the total pressure ratios at the outlet and inlet, respectively; SM refers to the surge margin; and refer to the total pressure ratios at the surge point and design speed, respectively.

*(2) Optimization Strategy and Control Parameters Setting*. In the first phase of optimization, based on the “repeated stage design” concept, all rotor blades and all stator blades are linked separately, and the parameterization method based on the stacking line is used for blade modification.

The optimization control parameters are set as follows: using the method of stacking line parameterization to set the bending, sweeping, and stagger angle rotating, the rotor blade and the stator blade are linked, and there are total 18 design variables in three radial stations (the rotor blade is 20%, 60%, and 100% of the radial height, and the stator blade is 0%, 50%, and 100% of the radial height). After repeated attempts, the range of sweeping and bending of rotor and stator blades is [−2, 2] mm, and the change range of stagger angle rotating of the rotor and stator blades is [−1, 1] rad. The experiment times of the optimal Latin square are 150. Each optimization task needs 27 cores in parallel, and each iteration needs to calculate the optimization operating point at the design speed and off-design speed, with a total of 300 flow field calculations. Due to supercomputing constraints, 10 tasks can be assigned at a time in parallel (without queuing). The time of each optimization task is 1 h, and the total optimization time is 35 h.

*(3) Comparative Analysis of the Optimization Results and Flow Field in the First Phase*. Table 1 shows the performance comparison in the first phase between the optimized blade and the original blade at the design speed and off-design speed. Under the condition of satisfying the constraints, the optimized efficiency at design point is increased by 0.3%, while the surge margin is 5.7%; at the off-design speed, the optimized efficiency at work point is increased by 0.66%, and the surge margin is 1.5%.

In the first phase, the surge margin at the off-design speed is reduced, so in the second phase we focus on improving the surge margin at the off-design speed. This paper analyzes the flow field near the surge point at the off-design speed to determine the blade row with large loss and optimize it. Figure 9 shows the comparison of static pressure distribution at the hub, middle, and tip sections after the first phase optimization of the near-surge point (260,000 Pa back pressure) at the off-design speed. Note that the static pressure here is the dimensionless static pressure obtained after dividing the absolute value of the static pressure by 101,325. All following static pressure distributions are treated in this way. Figure 10 shows the relative Mach number distribution at three sections (hub, middle, and tip) after the first phase optimization of the near-surge point at the off-design speed. Figure 9 indicates that each blade row inlet at each section has a positive incidence angle, where the middle section of R2 and the middle and tip sections of R1 have the largest positive incidence angle. With the increase of the back pressure, the positive incidence angle will further increase, which will cause the air flow to separate in a large range at the back of the blade, resulting in surge. From Figure 10, we can see that a large low-speed separation area appears in the middle and tip channel of R1, thus increasing the loss there. In the middle section of R2, there is a large range of low-speed areas in the channel, where there will be a large loss. It can be seen from the above analysis that the control points in the middle and tip of R1 and R2 should be selected as the optimization variables in the second phase.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

###### 5.2.2. Optimization in the Second Phase

*(1) Optimization Objectives and Constraints in the Second Phase*. To ensure that the second phase can effectively improve the surge margin at the off-design speed on the premise that the performance at the design speed is still The purpose of optimization is as follows: to improve the surge margin at the off-design speed in the second phase, and at the same time, to ensure that the efficiency at the design speed is not lower than the original blade. In order to achieve this goal, the optimization objectives are set as follows:

Establish the constraint condition by where refers to the efficiency of working point (220,000 Pa back pressure) at the off-design speed and refers to the efficiency of near-surge point (260,000 Pa back pressure) at the off-design speed. andare the calculated efficiencies of the working point and near-surge point (260,000 Pa back pressure), respectively, at the off-design speed. and are the weighting coefficients of two operating points and both are set to 1. and are the new total pressure ratio in the optimization process and the original total pressure ratio at the design speed, and are the total pressure ratio of optimized operating point at the off-design speed, and are the new flow rate and original flow rate at the design speed, and and are the new flow rate and the original flow rate at the off-design speed. refers to the efficiency in the second phase and refers to the efficiency of the original blade at the operating point of 340,000 Pa back pressure at the design speed. This constraint ensures that the efficiency at the design speed is larger than that of the referenced blade but can be lower than that after the first phase.

*(2) The Relevant Optimization Parameter Settings in the Second Phase*. According to the flow field analysis of the optimized blade in the first phase, R1 and R2 are selected as the modified blades in the second phase. To ensure that the flow rate change in the optimization process will not exceed the constraints, the leading edge control points of R1 and R2 are fixed to ensure that the leading edge metal angle remains unchanged and the control variables of blade body are selected. The rest of the blade rows are fixed as the optimized blades obtained in the first phase, and the distribution of the control variables of R1 and R2 is shown in Figure 11.

As can be seen in Figure 7, the second phase adopts the IABC algorithm to find the global optimal solution. The engineering optimization pursues the engineering optimized solution rather than the theoretical global optimal solution. Therefore, the colony size and number of generations are generally determined according to the size of the design variables, and the upper limit of the number of iterations is used as a marker to exit the cycle.

In the second phase, the full-blade surface parameterization method is used with 12 optimization variables. After repeated attempts, the range of R1 and R2 blade control parameters is set as [−3.0, 3.0] mm. In IABC algorithm, the bee colony size is 100, and the number of iterations is 3. Approximately 300 calculations are performed on the supercomputing platform; each optimization task needs 27 cores in parallel and can be concurrent with up to 10 tasks at a time (when the task does not queue on the supercomputing platform); and the optimization time of each optimization task is 1.5 hours (each task needs to calculate three operating points, 340,000 Pa back pressure at the design speed, 220,000 Pa back pressure and 260,000 Pa back pressure at the off-design speed), with a total time of 56 h.

###### 5.2.3. Optimization Results and Analysis

Table 2 shows the performance comparison before and after multicondition optimization in the second phase. Compared with the referenced blade, the optimized efficiency in the second phase increases by 0.3% at the design speed, and the surge margin increases by 4.0%; the efficiency increases by 0.8% at the off-design speed, and the surge margin increases by 2.3%. Compared with the optimized solution in the first phase, the surge margin in the second stage decreased by 1.7%; the efficiency at working point and the surge margin increased by 0.1% and 3.8% at the off-design speed. In the second phase, although the efficiency at the design speed is slightly reduced, the surge margin at the off-design speed is greatly improved, and the four optimization objectives under multiple working conditions are realized. Therefore, the selection of this geometry as the final optimized solution meets the design requirements.

Figures 12 and 13 show the comparison of the aerodynamic performance after two phases at the two speeds, respectively. It can be seen that the flow rate change range of the original blade and the optimized blade at the design speed is almost 0, and the optimized pressure ratio range is larger than the original one; from Figure 13, it can be seen that the efficiency of the above range of the pressure ratio at the working point is larger than that of the original one at 0.8 reduced speed, thus improving the surge margin of compressor at 0.8 reduced speed of the four-stage low-pressure compressor.

**(a)**

**(b)**

**(a)**

**(b)**

The geometry of typical compressor blades before and after optimization is selected for comparison, as shown in Figure 14. It can be seen that, compared with the original blade, the variation trend of optimized R1 and R2 is basically the same. In the hub region, the optimized blade is slightly forward swept; in the middle region, the optimized blade camber angle and leading edge metal angle are slightly reduced; in the tip region, the leading edge metal angle of optimized blade is significantly reduced, and the camber angle of blade profile is slightly increased. Because all blades except R1 and R2 are fixed in the second phase, the geometric changes of R3 and R4 are the same as those of R1 and R2 in the first phase: the stagger angle of hub blade profile decreases, the stagger angle of middle blade profile decreases slightly and presents the characteristics of reverse bending, and the stagger angle of tip section increases. Optimized S5 represents the change of all stator blades: the hub profile sweeps back, the middle profile bends slightly, the tip profile sweeps forward, and the stagger angle increases slightly.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

**(k)**

**(l)**

To further analyze the reasons for the performance improvement of the optimized blade at the design speed, a typical flow field before and after optimization is selected for comparative analysis. It can be seen from Figure 15(a) that the supersonic Mach number in the channel of optimized R1 decreases slightly, the shock wave intensity decreases, and the corresponding shock wave loss and the loss of shock wave interaction with the boundary layer decrease. The decrease of shock wave intensity in this region is due to the decrease of leading edge metal angle, which leads to the decrease of positive incidence angle of inlet flow. At the same time, due to the slight reduction of the blade camber angle, the reverse pressure gradient in the middle section of R1 decreases, causing the low-speed region near the trailing edge to be slightly reduced and the corresponding flow separation loss to be reduced. Similarly, there is the flow field of R2. From the comparison of S5 flow field before and after optimization in Figure 15(b), it can be found that, because the stagger angle of optimized S5 increases slightly, the positive incidence angle of the incoming flow decreases, and the aerodynamic load of the whole blade profile decreases, so that the low-speed region at the trailing edge is obviously reduced, and the corresponding flow separation loss is reduced. All the optimized stator blades have similar flow improvement to the optimized S5.

**(a)**

**(b)**

As shown in Figure 16, the limiting streamline and static pressure distribution of R2 and R3 before and after optimization at the design speed are compared. Because the leading edge metal angle of the optimized R2 blade at the tip is reduced, the positive incidence angle of the incoming flow is reduced, so that the shock wave position of the blade profile is slightly pushed back on the suction surface near the tip, as shown in Figure 16(a), and the corresponding flow separation loss is reduced. In Figure 16(b), due to the decrease of the leading edge metal angle and the increase of the negative incidence angle of the incoming flow, the position of the shock wave on the suction surface of optimized R3 is pushed back, and the separation area behind the shock wave is reduced, which is shown in the streamline diagram; that is, the vortex area is reduced, thus reducing the corresponding flow separation loss.

**(a)**

**(b)**

The comparison of typical flow field before and after optimization of working point (220000 Pa back pressure) at the off-design speed is shown in Figure 17, and the most significant flow field change is at the tip region. The results show that the Mach number on the suction surface decreases obviously with the change of the tip profile of R1, and the corresponding loss of the shock wave and the loss of interaction between shock wave and boundary layer are reduced, similar to R2, R3, and R4.

It can be seen from the comparison of the limiting streamline diagram in Figure 18 that optimized R1 significantly shortens the length of the radial separation line at the off-design speed and reduces the size of the separation bubble, thus reducing the corresponding separation loss. In addition, the optimized R1 also reduces the open separation zone near the tip, causing the decrease of the corresponding air loss, while the radial separation of other rotor blades is not obvious.

**(a)**

**(b)**

Figure 19 shows the comparison of entropy distribution in the channel before and after optimization at the off-design speed. Due to the decrease of the stagger angle in the hub area of optimized blade, the positive incidence angle of airflow at the off-design speed decreases, the airflow in the main flow area is smoother, the low entropy area increases, and the corresponding airflow loss decreases.

The following is an analysis of the causes for the increase in the surge margin at the off-design speed. Since the real surge point is difficult to obtain accurately by numerical calculation, we believe that the change in flow field at the near-stall point can predict the change in the surge point. Figure 20 compares the relative Mach number distribution of the near-stall point of the blade-tip section at the off-design speed before and after optimization, from which we can see that, compared to the optimized blade, the original tandem blades have a large range of low-speed regions. The optimized blade has a significant reduction in the low-speed region and enhanced circulation capacity, while the corresponding airflow loss is reduced.

**(a)**

**(b)**

Figure 21 compares the static pressure distributions at the tips of the first two rows of rotor blades before and after optimization. It can be seen that the shock wave of the optimized R1 pushes back and the shock wave intensity is slightly weakened. The optimized R2 has a significant reduction in shock wave intensity, thus significantly reducing the loss of shock wave at this location and the loss of interaction between the shock wave and the boundary layer.

From the limiting streamline comparison of the suction surface of the rotor blades shown in Figure 22, it can be seen that the vortex at the tip of the optimized R1 is significantly weakened, such that separation loss at this location is reduced. At the same time, the forward-flow region of the optimized blades increases, and the separated airflow achieves reattachment and forms separation bubbles at 80% to 90% of the radial area after separation. Meanwhile, the original blade suction surface has a wide range of open separation, so the optimized R1 effectively reduces the separation loss. Similar to the situation at the tip of optimized R1, the vortex at the tips of the optimized R2 and R3 is significantly weakened. The separation line at the hub region near the trailing edge of the optimized R2 is significantly shortened. The improvement of flow in these areas reduces the corresponding losses and leads to an increase in the flow capacity of the optimized blades.

**(a)**

**(b)**

**(c)**

In summary, after two phases of optimization, the performance of the four-stage low-pressure compressor at the design speed and off-design speed has been improved. The performance improvement at the design speed is mainly due to the reduction of the Mach number at the middle and tip regions of R1 and R2 and the reduction of the separation area in the exit cascade. However, due to the limited degree, the efficiency improvement is limited. The improvement of the performance at the off-design speed is mainly due to the weakening of the intensity of shock wave at the tip section of the R1 and R2 and the reduction of the radial secondary flow separation of the optimized R1. The improvement in the surge margin of the optimized blades is mainly due to reductions in the separation areas of the suction surfaces of the first three rows of rotor blades and a significant reduction in the low-speed zone of the tandem blades.

#### 6. Conclusion

Based on the phased parametric optimization strategy and IABC algorithm, a multitask concurrent optimization system is built on the supercomputing platform to perform the optimization for the AL-31F four-stage low-pressure axial flow compressor under multiple working conditions, and the following conclusions are obtained:(1)The phased parameterization strategy can greatly reduce the dimension of multistage axial compressor aerodynamic optimization and shorten the optimization time. In the first phase, the parameterization method based on the stacking line and “blade-row-linkage” is used to search the mode and explore the design space by a large step. In the second phase, refined exploitation is adopted, and good results are achieved. At the same time, the optimization strategy is mainly embodied in the selection of the parameterization method and optimization algorithm, the selection of optimization control variables, and the introduction of expert experience after the first phase optimization, which avoids invalid exploration in the optimization process.(2)The multitask concurrent optimization system based on IABC algorithm and supercomputing platform is useful optimization system for multistage compressors, and the calculation efficiency can be increased by more than 10 times. The total optimization time of the two phases is 91 hours, which is within the acceptable time range in engineering.(3)The proposed optimization method is adequate for the multiple condition optimization of an AL-31F four-stage low-pressure compressor. Under the constraint conditions, the adiabatic efficiency at the design speed is increased by 0.3%, and the surge margin is 4.0%; the adiabatic efficiency and the surge margin at the off-design speed are increased by 0.8% and 2.3%. This method has wide application value in the field of the engineering optimization of multistage axial compressors.(4)The performance improvement of AL-31F at the design speed is mainly due to the reduction of Mach number in the middle and tip sections of R1 and R2, as well as the reduction of separation loss at the exit of outlet cascade; the performance improvement at the off-design speed is mainly due to the reduction of the loss of shock wave, caused by the decrease of Mach number in the tip sections of R1 and R2, and the significant reduction of secondary flow on the suction surface of R1. The improvement in the surge margin is mainly due to decreases in the separation areas of the suction surfaces of the first three rows of blades and a significant decrease in the low-speed area of the tandem blades.

#### Abbreviations

HEB: | High-dimensional, expensive, and black-box |

PDE: | Partial differential equations |

CST: | Class-shape function transformation |

IABC: | Improved artificial bee colony |

ABC: | Artificial bee colony |

GA: | Genetic algorithm |

3D: | Three-dimensional |

Ori: | Original |

Opt: | Optimized |

: | Chordwise direction |

: | Spanwise direction |

mm: | Millimeter |

rad: | Radian |

S-A: | Spalart‐Allmaras |

TPR: | Total pressure ratio |

eff: | Efficiency |

SM: | Surge margin. |

#### Data Availability

The calculated data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.