Design Optimization of aMultistage Axial Flow Compressor Based on Full-Blade Surface Parameterization and Phased Strategy

A new design optimization method is proposed for the problem of high-precision aerodynamic design of multistage axial compressors. The method mainly contains three aspects: full-blade surface parametrization can significantly reduce the number of control variables per blade row and increase the degrees of freedom of the leading edge blade angle compared with the traditional semiblade parametric method; secondly, the artificial bee colony algorithm improved initialization and food source exploration and exploitation mechanism to enhance the global optimization ability and convergence speed, and a distributed optimization system is built on the supercomputing platform based on this method; finally, a phased optimization strategy based on the “synchronous change in multirow blades” is proposed, and expert experience is introduced to achieve a better balance between exploration and exploitation. The optimization method is applied to the AL-31F four-stage low-pressure compressor. As a result, the adiabatic efficiency is improved by 0.67% and the surge margin is improved by 3.1% under the premise that the total pressure ratio and mass flow rate satisfy the constraints, which verifies the effectiveness and engineering practicality of the proposed optimization method in the field of multistage axial flow compressor aerodynamic optimization.


Introduction
Aeroengines and gas turbines are needed for strategic national development, and one of the core components of aeroengines and gas turbines is the multistage axial flow compressor. For the design of multistage axial flow compressors, an aerodynamic optimization design method can effectively reduce the dependence on manual experience and greatly improve the design capability. It is well known that the aerodynamic optimization design of a multistage axial flow compressor has the characteristics of typical high-dimensionality, expensive cost, and black box (HEB) [1] problem. With an increase in optimization control variables, the design space will grow exponentially and fall into the "curse-of-dimensionality." Therefore, obtaining the optimized solution within the acceptable time range of engineering is a frontier problem in the field of aero-dynamic optimization of multistage axial flow compressors, which has a very important engineering application value.
Before the year 2000, limited by the computing power at that time, aerodynamic optimization of the multistage axial flow compressor was mainly limited to one dimension and two dimensions [2][3][4][5]. With the great improvement in computing power, three-dimensional steady aerodynamic optimization of compressors has been developed rapidly, and a few three-dimensional unsteady optimization methods have even been proposed [6]. However, the cost of unsteady calculation is too large, much simplification is required, the actual calculation accuracy is limited, and there is no engineering practicability. At present, research on multistage axial flow compressors still mainly focuses on three-dimensional steady optimization. With the development of optimization regarding multidiscipline and high-fidelity problems, the HEB problem mentioned above has increasingly become the key to restrict the design of multistage axial flow compressors.
Traditional methods for solving the HEB problem of aerodynamic optimization of multistage axial flow compressors can be divided into six categories: (1) screening variables [7,8], (2) stage-by-stage optimization [9], (3) spatial decoupling [10], (4) surrogate model method [11][12][13][14][15], (5) adjoint algorithm [16][17][18], and (6) parametric dimensionality reduction [19,20]. Among them, the screening method is only suitable for examples with obvious bad flow fields but is not universal. The stage-by-stage optimization process is complex, and the coupling effect of the flow field between the front and back stages is ignored. The spatial decoupling method proposed in literature [10] is only suitable for the case of a weak radial flow but not for a blade with a low aspect ratio. The surrogate model method is difficult to train in high-dimensional problems, and the number of samples needed increases exponentially with the dimensionality, so the sample training cost is difficult to handle. Although the adjoint algorithm has the advantage that optimization is independent of the dimensionality of variables, it easily falls into local extrema [21] and has defects in handling multiconstraint multiobjective problems.
The literature [22,23] notes that one of the most promising strategies for solving the HEB problem is parametric mapping dimensionality reduction. For aerodynamic optimization of multistage axial flow compressors, the aim is to find a parametric method to reduce the design variables and achieve dimensionality reduction while retaining as much of the optimal solution from the original design space as possible. The traditional parametric method has remarkably high-dimensional characteristics because each section of the blade is modeled separately and independently of each other [24]. An effective method to constrain the control parameters of each section is the surface parametric method, which was first applied in outflow aerodynamic optimization and is an effective dimensionality reduction method [25][26][27][28]. Lee et al. compared B-spline and CST surface methods and concluded that the two parametric methods have their own advantages and disadvantages and are suitable for different situations. In 2003, the semiblade surface parametric method (blade suction surface modification by a Bezier surface) was applied in the internal flow field of compressors by Burguburu and le Pape [29], which was further developed by Wang and Zhou [30] and Cheng et al. [1]. However, the defect of this surface parametric method is that only half of the blade surface is modified, and the influence of the other half surface on the flow field is ignored. In addition, this method cannot change the blade angle at the leading edge; this defect plays an important role in the development of the flow field as it results in limiting the improvement space of the aerodynamic performance for blade geometry. In this paper, a full-blade surface parametric method is presented that regards the suction surface and the pressure surface as a whole surface and superimposes the perturbed surface on it. This method reduces the control variables and realizes an effective dimensionality reduction under the condition of ensuring the degree of freedom of the leading edge, trailing edge, and blade body.
Another important part of the optimization process is the optimization algorithm. To avoid falling into a local extremum, this paper adopts an improved artificial bee colony (IABC) algorithm with good global optimization ability [31], which does not depend on the initial solution. In addition, due to the improvements in initialization, food source exploration, and exploitation mechanisms from the standard artificial bee colony [32] (ABC) algorithm, its convergence speed is improved. Because the evolutionary algorithm has the advantage that it does not depend on the gradient information of the function, it can perform a simultaneous search of multiple individuals in the design space with the potential for multitask concurrency. Combined with a commercial supercomputing platform with powerful parallel capabilities, the IABC algorithm can greatly reduce the optimization time.
However, because evolutionary algorithms essentially require many high-fidelity performance evaluations, the above strategy is still not enough to reduce the computational cost to an acceptable range. This paper proposes a phased optimization strategy based on the "synchronous change in multirow blades," which can greatly reduce the time cost of optimization. The core of this strategy is to draw on the idea of "exploration first, exploitation later" in the optimization process and divide the entire optimization process into two phases to effectively reduce the number of iterations required for the spatial search, thereby reducing the optimization time.
The structure of this paper is as follows: Section 2 shows the optimization method; Section 3 performs the optimization case verification and analysis; Section 4 draws a few conclusions.

Optimization Method Based on Full-Blade
Surface Parameterization 2.1. Parametric Method. The traditional parametric method shapes each radial section of the blade separately; on the one hand, it makes the blade profile difficult to smooth; on the other hand, it makes the blade control parameters numerous (a single row of blades can have hundreds of parameters). To better promote the solution of the HEB problem of multistage axial compressor aerodynamic optimization, the surface parametric method with significant lowdimensional characteristics was selected. This method understands the blade as a combination of the suction surface and pressure surface rather than the accumulation of twodimensional sections. This paper proposes a full-blade surface parametric method that is an improvement of the semiblade surface parametric method [29]. Figure 1 shows the construction process of the method. The blade is opened radially from the leading edge to form an unfolded surface, and the Bezier perturbed surface is superimposed to form a new unfolded surface. Then, the new unfolded surface is stitched at the leading edge in the radial direction to generate a new blade. The parametric method can be divided into the following four steps: (1) The blade opens radially at the leading edge, forming an unfolded surface 2 International Journal of Aerospace Engineering (2) Chord length parameterization is applied. Because the three-dimensional surface of the physical space and the unit Bezier surface are to be superimposed, the coordinates of each point of the blade in the physical space need to be converted to the unit plane in the calculation space, and ξ and η are the axes of the computational space. The formula for chord length parameterization is as follows: where i ∈ ð1, n p Þ and j ∈ ð1, n s Þ; n p and n s refer to the number of data points of a section and the total number of sections, respectively; T l j refers to the sum of the chord lengths of each segment of the jth section in the η direction; T L i refers to the sum of the chord lengths of the ith section in the ξ direction; l m is the mth chord length of the jth section in the η direction; and L n is the nth chord length of the ith section in the ξ direction (3) Bezier perturbed surface generation is applied. The perturbed surface is defined by the following formula: where R ! refers to the change value of each point on the perturbed surface; P k,l is the ðm + 1Þ × ðn + 1Þ control points of the perturbed surface; B m l ðvÞ and B n k ðuÞ refer to the Bernstein basis functions calculated by equation (4); v and u refer to the independent coordinate variables of the perturbed surface, and its variation range is ½0, 1; and C n k is the combination number calculated by equation (5  The full-blade surface parametric method is a new surface parametric method that retains the high-order smoothness characteristic while improving the existing semiblade surface parametric method, which is mainly manifested in the following two points: (1) Since the full-blade surface parametric method combines the suction surface and the pressure surface as a whole, it can take into account the influence of the suction surface and pressure on the aerodynamic performance while not significantly increasing the number of control variables, and the low-dimensional characteristics are maintained well. The control parameters of a single blade row can be reduced to 16 (2) The full-blade surface parametric method increases the degree of freedom of the geometric variation at the leading edge. The semiblade surface parametric method lacks the freedom of the blade angle at the leading edge, so it has certain restrictions on the geometric exploration function. Figure 2 shows a schematic diagram of the point distribution of the full-blade parametric method. The yellow point is a free active point, and the blue point is a limited active point; that is, all the blue points need to keep changing synchronously so that the leading edge of the suction pressure can be smoothly connected 2.2. Optimization Algorithm. For the purpose of searching for the global optimal solution of the multistage axial flow compressor, an evolutionary algorithm with global search ability is adopted. Since evolutionary algorithms do not require the gradient information of function values, multiple individuals can be concurrent, which can greatly reduce the optimization time. Among various evolutionary algorithms, the ABC [32] algorithm has a better optimization accuracy and convergence speed. In this paper, the ABC algorithm is developed through an initialization mechanism, a food source exploration, and an exploitation mechanism, and its schematics are shown in Figure 3. The IABC [31] algorithm involves employed bees, onlooker bees, and detector bees. At the beginning, the employed and onlooker bees constitute the entire bee colony, and the number of bees is divided equally between these types. The food source in the schematic refers to a set of feasible solutions, and the food source concentration refers to the fitness of the feasible solutions. The detailed steps of the IABC algorithm are as follows: (1) Food source initialization. To obtain as much design space information as possible from as few sample points as possible, instead of the random initialization mechanism in the standard ABC algorithm, the optimal Latin square experiment method [33] is utilized to initialize the distribution of the food source, and the feasible solution X i ði = 1, 2, 3,⋯,NÞ of ND-dimensional vectors is generated (2) The concentration of the initial food source is calculated, the best value is recorded, and the concentra-tion is ranked according to the food source. The first half of the food source's fitness is taken as the target of the employed bees, and the rest is left to the onlooker bees (3) According to formula (1), all bees explore the vicinity of the initial source and calculate the concentration of the food source. When the concentration of the new food source is greater than that of the original one, the original food source will be replaced by the better one. If the opposite is true, the bees will continue to explore where V j i refers to the jth component of the food source location of the ith bee, j ∈ f1, 2,⋯,Dg ; k ∈ f1, 2,⋯,N/2g, i ∈ f1, 2,⋯,N/2g and k ≠ i ; ξ j i and λ j i are random numbers between ½−1, 1 ; and X j best is the best food source location in the general information of the food source (4) Onlooker bees explore the food source. The onlooker bees use the Russian roulette rule to select food sources, as shown in formula (2); that is, the food source is selected with a probability proportional to the concentration of the food sources determined by the employed bees and is explored in the vicinity according to formula (3). If the concentration of the obtained food source is higher than the concentration of the original food source, the original food source is replaced, and if not, the exploration is continued where the value of N e is N/2, referring to the number of employed bees. P refers to the probability of selecting a certain food source; f ðX i Þ refers to the concentration of the ith and f ðX m Þ refers to the mth food sources; ðX j Neib Þ best refers to the food source position that has the largest concentration in the adjacent area; and ðV j Neib Þ best refers to the location of the new food source of onlooker bees. The Chebyshev distance formula is used to calculate the best point in the neighborhood, which is shown in formula (4): where j ∈ f1, 2,⋯,Dg and dði, tÞ is the Chebyshev distance between X i and X t . The neighborhood of X i is determined by formula (5), which means that if the Chebyshev distance between X i 4 International Journal of Aerospace Engineering and X t is less than or equal to the product of r (the neighborhood radius) and md i (the average Chebyshev distance), then X t locates in the neighborhood of X i ; otherwise, it does not: where md i is the average Chebyshev distance between X i and the entire onlooker bee population, S is the neighborhood of X i , and r is the radius of the neighborhood. When r is 1, the algorithm has the best convergence.
ðX Neib Þ best is calculated by where fitðÞ refers to the individual fitness of each bee. For the purpose of testing the performance of the above algorithm, the ABC algorithm and the genetic algorithm [34] (GA) are used for comparative analysis.
The algorithm parameters of IABC and ABC are set the same, and the total number of bee colonies is 200, among which the number of employed bees is 100, the number of    Table 1. Table 1 shows that the IABC algorithm has obvious advantages for the first three benchmark functions. The minimum value obtained is approximately 10 orders of magnitude smaller than that for the other two algorithms. For the Rosenbrock function, the ABC algorithm is optimal, but the IABC algorithm is at the same level, and the disadvantage is not obvious. For the Zakharov function, the GA is optimal, but the IABC algorithm is suboptimal by less than an order of magnitude; hence, the disadvantage is acceptable. Figure 4 shows the convergence history of these benchmark functions. It can be seen that IABC has a better and faster convergence performance in the first four benchmark function tests than the other algorithms.

Optimization Strategy.
Engineering optimization often does not require an optimal solution, but only an optimized solution that is better than the original one. In this case, a certain optimization strategy is needed to shorten the optimization cost. The advanced parametric methods and optimization algorithms introduced in Sections 2 and 3 reduce the dimensionality of the optimization task for the multistage axial compressor, thus reducing the optimization calculation cost to a certain extent. However, considering the high cost of the high-fidelity numerical calculation of multistage axial compressors, it is necessary to further adopt an appropriate optimization strategy to complete the optimization task within the acceptable time range of the project. This article uses a combination of the following three optimization strategies.

Coarse Grid Optimization and Fine Grid Verification.
The literature [36,37] notes that the coarse mesh has been able to capture the near wall viscosity, so coarse and fine meshes can predict the same trend of performance variation during optimization. Although the computing power of computer is improving, it is still insufficient for the high-precision design optimization of multistage axial compressor. Therefore, it is still necessary to adopt certain optimization strategies to reduce the time consumption of single flow field calculation. "Coarse mesh optimization, fine mesh verification" is a common and effective strategy to reduce the time consumption of single flow field calculation. The flow field calculation in the optimization process uses a coarser grid, and after obtaining the optimization solution, a fine grid is used to check the calculation to obtain the optimized solution under the high-fidelity flow field calculation.

Multitask Concurrent Parallel
Computing Based on the IABC Algorithm and Supercomputing Platform. Since evolutionary algorithm optimization does not depend on the gradient information of the function, it naturally possesses multitask concurrent computing capabilities. Using China's advanced "Taihu Lake Supernatural Light" supercomputing platform, a multitask concurrent optimization system based on the IABC algorithm is built, which can reduce the optimization time considerably. Figure 5 shows a schematic diagram of the multitask concurrent system based on the IABC algorithm and the supercomputing platform. On the user's computer, the initial blade is used as the input of the program to start the optimization algorithm; initialize the food source distribution; iteratively conduct employed bee, onlooker bee, and detector bee optimization; and record the location of the best food source for each bee colony generation (i.e., the optimal solution for each generation). If the loop exit condition has not been met, then the loop is continued. In any loop iteration, the initial food source, employed bees, onlooker bees, and detector bees all adopt a concurrent mechanism, so the login node needs to be started as an intermediary. Since there is an upper limit for supercomputing concurrence, the number of concurrent bee colonies must be limited to Max_Run; that is, Max_ Run bees are concurrently processed. After counting one bee, the next new bee is started to join the parallel. The process of concurrent bees is shown on the user's computer as Max_Run. New geometry files are generated through the parametric model and then placed in the corresponding Max_Run new engineering task folder. At the login node, multitask files are generated through scripts, and then, multitask files are assigned. The function of the multitask file is to start a remote computing node to perform concurrent computing. On the computing node, each task runs independently and must complete three steps, namely, grid generation, flow field calculation, and postprocessing analysis. On the login node, the script responsible for regularly detecting the status of the computing node runs at all times. If an engineering task in the computing node ends, performance results will be written in donelist.txt, and a line of donelist.txt will be added. At this time, by regularly synchronizing donelist.txt on the user's computer and the login node, the operation status on the local user's computer is known, and a decision on whether to dispatch a new bee can be made. Eventually, the entire optimization cycle reaches the set number of iteration steps and exits the cycle to obtain the optimized blade for this case.        International Journal of Aerospace Engineering evolutionary algorithm, this paper proposes a phased optimization strategy based on the "synchronous change in multirow blades"; that is, the aerodynamic optimization process of the multistage axial compressor is divided into two phases, as Figure 6 shows. The first phase: the "synchronous change in multirow blades" strategy is adopted; that is, for the purpose of exploring the design space at a large step, each blade row of the multistage axial compressor is modified in the same way and by the same amount. After the optimization reaches a certain number of steps, which is determined according to the designer's experience, the decision of when to end the first phase of the optimization exploration is made, and the optimized blade of the first phase is obtained. The first phase focuses on the exploration of the design space.
After the first phase, the optimization process comes to a temporary end and the existing optimization results are ana-lyzed to identify the most promising blade rows and corresponding blade regions (largest loss in this area). These blade rows and blade regions are selected for the second phase of automatic optimization. In this analysis and selection process, the experience of the designer is incorporated.
The second phase: it is the independent optimization of the above selected potential blade rows. The aim is to build on the results of the first phase of exploration to enable further exploitation of the corresponding design space to obtain better performing blade shapes.
The phased optimization strategy based on the "synchronous change in multirow blades" can play a practical role in the aerodynamic optimization of multistage axial compressors for the following two reasons: (1) The optimized solution only exists in a few small local areas in the huge design space. In the evolutionary  9 International Journal of Aerospace Engineering algorithm, the change in each blade row in each optimization iteration is very small, and it is difficult to quickly explore the design space. If only the number of variables of each optimization in the algorithm is increased, then it is easy to increase the possibility of the mutual cancellation of the effects of different variables, but it is also difficult to quickly explore the space. The "synchronous change in multirow blades" strategy can increase the pace of exploration and greatly increase the probability of finding the optimized solution neighborhood within a limited number of iterations (2) After the first phase, the designer's experience is utilized to analyze the flow field so that fewer optimization variables need to be used for targeted optimization exploitation in the second phase, avoiding invalid exploration of the design space and improving the optimization efficiency

Optimization Case Verification and Analysis
3.1. Optimization Object. The optimization example used in this paper to verify the effectiveness of the above optimization method is the AL-31F four-stage axial low-pressure compressor, with a total of ten rows of blades, as shown in Figure 7. The inlet is a guide vane, and the outlet is tandem. The reason for using a tandem blade in the last stage of the stator blade is to ensure that the outlet airflow is along the axial direction and at the same time to avoid large separation. If only one stator blade is used, the blade-camber angle of that stator blade will be too large, and thus, the airflow will be prone to large separation at the back of the blade, causing great losses and even stalling and surging. Compared with the performance of international advanced low-pressure compressors, as shown in Table 2, it can be seen that the average stage pressure ratio, adiabatic efficiency, and surge margin of the AL-31F four-stage low-pressure compressors have a certain gap compared with other advanced low-pressure compressors. Therefore, there is some room for improvement in the AL-31F's aerodynamic optimization design.

Numerical Method Verification.
A CFD numerical simulation uses the Fine module of NUMECA to solve the threedimensional Reynolds average N-S equation in finite volume form. The Autogrid5 module of NUMECA is used to generate the grid, and the grid topology of the flow passage and blade is established automatically. In this test case, the model without tip clearance and the Spalart-Allmaras one-equation turbulence model are adopted. The time discretization is  The calculation domain and boundary condition setting for four-stage low-pressure compressor are shown in Figure 8. The total temperature at the inlet boundary is 293.15 K, the total pressure is 101,325 Pa, and the inlet air flow direction is axial. The boundary of the solid wall is adiabatic and slip-free. Given the back pressure at the outlet boundary, by gradually adjusting the outlet back pressure during the calculation process, it will advance towards the near stall point and the blocking point. The near stall point is the condition before calculating divergence that arises from increasing the back pressure of the design point, and the blocking point is the condition before calculating divergence that arises from reducing the back pressure of the design point. Figure 7 shows the comparison between the characteristic curve of the low-pressure compressor at the design speed obtained by numerical simulation and the experimental data. At the design speed, the maximum pressure ratio measured in the experiment is 4.05, the maximum pressure ratio obtained by simulation is 4.18, and the simulated value is 3.2% larger than the experimental value. Figure 7(b) illustrates that the simulated value of the maximum efficiency is  11 International Journal of Aerospace Engineering 1.5% larger than the experimental value. At the design speed, the calculated pressure ratio range and efficiency range are larger than the experimental values, and the flow margin is small. The simulated design speed margin is 10.2%, and the experimentally measured margin is 14.9%. However, the calculated characteristic curve of the low-pressure compressor is consistent with the trend of the experimental curve, which proves the validity of the numerical simulation method.
To make the verification results more credible, grid independence verification is performed on the optimized objects. Figure 9(a) shows the grid of the three-dimensional blade, Figure 9(b) shows the B2B grids of hub section, and Figure 9(c) shows the calculation performance of four sets of grids under a backpressure of 300,000 Pa. The mesh quality satisfies the near wall condition y + < 2. As seen from Figure 9(c), when the number of single channel grids in each row of blades reaches approximately 1 million, the changes in efficiency and pressure ratio are very small, so it is reasonable to believe that the third set of grids (on average, 930,000 per single passage of a blade) has met the grid independence requirements. However, by considering further shortening the optimization time, the "coarse grid optimization, fine grid verification" strategy is adopted, in which the first set of grids (average 250,000 per single passage of a blade) is used in the optimization process and the third set of grids is used for the flow field calculation of the optimization results. Figure 10 shows the optimization process. First, the flow field of the initial blade is calculated. If the performance does not meet the requirements, it enters into the optimization cycle. The IABC algorithm is adopted, which includes four steps: initialization, exploration of employed bees, exploration of onlooker bees, and establishment of detector bees. Each step will produce a corresponding number of new food sources, and each food source position represents a new set of control parameters. A new blade is generated under the full-blade surface parametric control method, and then, mesh generation, flow field calculation, and performance analysis are carried out for the new blade. If the requirements are met, the optimization will exit the cycle, and the optimized blade will be obtained; otherwise, the iterative optimization of the cycle will be continued. Usually, the decision regarding when to exit the cycle should depend on the designer's experience.

Optimization Settings and Result
Analysis. The three optimization strategies described in Section 4 are combined to conduct the optimization task of the four-stage lowpressure compressor in two phases. (1) Optimization Goals and Constraints in the First Phase. From the comparison of the above numerical simulation and experimental results, it can be seen that the efficiency and surge margin of the design speed of the AL-31F fourstage low-pressure compressor still have room for improvement. Therefore, the optimization objective of the first phase is mainly to maximize the efficiency of the design point, taking into account the surge margin improvement. According to the author's optimization experience, when 340,000 Pa of back pressure near the stall point is selected as the optimization condition and the adiabatic efficiency under the back pressure is taken as the objective function of the optimization process, an optimization solution that can improve the design point efficiency and surge margin can be obtained under the condition that only one flow field calculation is carried out in each optimization iteration. The flow rate constraint changes within 1%, and the total pressure ratio does  with respect to x subject to TPR − TPR ori ≥ 0 mass-mass ori mass ori ≤ 1:0% where eff refers to the adiabatic efficiency under optimized conditions; TPR and TPR ori refer to the total pressure ratio in the optimization process and the total pressure ratio of the original compressor, respectively; mass and mass ori refer to the flow rate in the optimization process and the flow rate of the original compressor, respectively; minus is the minimum value given by experts; and x i refers to the optimization variable, whose lower and upper limits are x L i and x U i , respectively. The solutions that do not satisfy the constraints are directly excluded from the IABC program, and the next generation of bees will not converge towards this food source, or the probability of convergence towards this source is almost zero. In this way, the function of constraint can be realized in an unconstrained algorithm.
The formulas for TPR, eff , and SM are as follows: where P out and P in are the total pressure ratios at the outlet and inlet, respectively; k is the specific heat ratio; T out and T in are the total temperature at the outlet and inlet, respectively; SM refers to the comprehensive surge margin; and TPR s and TPR 0 refer to the total pressure ratios at the surge point and design point, respectively.
(2) Optimization Settings in the First Phase. The full-blade surface parametric method is utilized, and all the blade rows are changed synchronously, with 16 design variables per row. The four parameters controlling the leading edge are within [-5.0, 5.0] mm, and the 12 parameters controlling the blade body geometry are within [-6.0, 6.0] mm. In the IABC optimization algorithm, the population size is 100, and the number of iterations is 8. The maximum number of mines is 3, which means that the maximum number of times the next generation of bees can explore at a particular previous generation of food source cannot exceed three times. That is, if the calculated fitness is smaller than the original one three times in a row, a new food source needs to be generated randomly again. On the supercomputing platform, approximately 850 calculations were carried out, and among them, the effective point was approximately 800 (the invalid points that did not meet the constraints were removed). Each optimization task needs 27 cores in parallel, and up to 24 tasks can be performed concurrently (when the task is not queued after supercomputing). The optimization time of each optimization task is 30 min; this process takes approximately 30 hours in total.
(3) Optimization Results and Comparative Analysis in the First Phase. Table 3 shows the performance of the fourstage low-pressure compressor at the design point with the design speed. It can be seen that under the premise of the "synchronous change in multirow blades" strategy in the first phase, the optimized blade meets the constraints of the flow and pressure ratio, and the efficiency increases by 0.48% and the surge margin 4.8%.
As shown in Figure 11, the static pressure and Mach number distribution in the hub, middle, and tip sections of the optimized solution in the first phase are presented. The losses in the R3 and R4 flow fields are not shown because they are not significant. It is illustrated that the positive incidence angles of airflow of R1 and R2 at the hub and middle sections are too large; the airflow at the leading edge is accelerated, thus increasing the shock wave intensity; and the maximum Mach number at the middle reaches 1.4, leading to a corresponding large loss at the tip section. At the same time, it can also be noted that there are two other relatively large low-speed areas: one is the tip trailing edge of S1 and the other is the outlet of S5 at the tip, thus causing corresponding loss. Therefore, these areas should be the target of the exploitation optimization in the second phase.

The Second Phase: Optimization Based on Exploitation Strategy
(1) Optimization Objectives and Constraints in the Second Phase. According to the flow field analysis of the optimized blade in the first phase, the control variables of the leading edge at the hub and the middle of R1 and R2, the leading edge at the tip of R3, the leading edge at the tip of S1, and the blade body at the middle and tip of S5 are selected as the design 13 International Journal of Aerospace Engineering variables of the second phase. The rest of the blade row is fixed to the optimized blade obtained in the first phase.
The optimization objectives, conditions, and constraints of the second phase are the same as those in the first phase.
(2) Optimization Parameter Settings in the Second Phase. After the flow field calculation and analysis of the optimized geometry obtained in the first phase, a total of five rows of blade geometry of R1, R2, R3, S1, and S5 are selected as optimized in the second phase. To reduce the optimization    International Journal of Aerospace Engineering variables for R1, R2, R3, and S1, the leading edge control parameters are optimized to reduce the positive incidence angle. To reduce the optimization variables for S1 and S5, the blade body is controlled in addition to the leading edge variables. Refer to Figure 12 for the specified control variables of each row: all the green points are leading edge control points, which are variable in the radial direction, but the four points of the same row must keep pace to avoid distortion at the leading edge; all the red points are blade-independent control points, and all the black points are fixed points.
In the second phase, the optimization parameters are set as follows: the full-blade surface parametric method is still adopted, with a total of 25 optimization variables. The 3 control parameters of the leading edge of R1 and R2 and the 2 control parameters of the leading edge at the hub of R3 all change within the range [-6.0, 0.0] mm. The variation range of the 8 control parameters of S1 and 9 control parameters of S5 is within [-6.0, 6.0] mm. The number of iterations is 10, the size of the bee colony in the IABC optimization algorithm is 100, and the maximum number of exploitations is 3. On the supercomputing platform, a total of approximately 1,000 calculations are performed, and the number of effective points is approximately 980 (without the invalid points violating the constraints). Each optimization task requires 27 cores in parallel, and a maximum of 24 tasks can be performed concurrently (under the condition that the tasks are not queued up). The optimization time of each optimization task is 30 min, and the optimization takes approximately 48 hours in the second phase.
The time consumption of the optimization method proposed in this paper is compared with the conventional method in Table 4. From Table 4, it can be seen that the proposed optimization method can greatly reduce the number of control variables in the optimization process, thus greatly reducing the times of flow field calculations required. At the same time, the optimization time can be reduced by an order of magnitude due to the construction of a supercomputing-based concurrent system. Although the optimization time of the DOE (design of experiment) method is not too long, it only has the function of global exploration and does not have the ability of local exploitation. The proposed method in this paper has the shortest optimization consumption time among several methods.
(3) Comparative Analysis of the Optimization Results in the Second Phase. Since a single-objective optimization (adiabatic efficiency) is adopted in the second phase, it cannot find a solution that is superior to the optimized solution found in the first phase in terms of adiabatic efficiency and the surge margin in the results. According to the designer's experience and needs, the point where the adiabatic efficiency is most improved and the surge margin is slightly decreased is selected as our final optimized solution. Table 5 shows a performance comparison of the optimization design point before and after optimization in the second phase. It can be seen from the table

18
International Journal of Aerospace Engineering that after two phases of optimization, the mass flow rate and pressure ratio at the design point of the optimized blade are almost unchanged, while the adiabatic efficiency increases by 0.67% and the surge margin 3.1%. Comparing the optimization results of the first phase and the second phase, it can be seen that the adiabatic efficiency of the design point in the second phase increases by 0.19%, while the surge margin decreases by 1.7%. Figure 13 shows a performance comparison between the optimized blade of the second phase and the original one. The change in the flow rate of the four-stage low-pressure compressor after optimization is almost 0. This is a common situation in multistage axial flow compressors. This situation indicates that although the back pressure is changing, there is a flow rate restriction in one or more stages in the multistage compressor, probably due to the fact that it has stalled there. We tend to expand the flow surge margin by adjusting the angle of multiple rows of stator blades in the normal operating condition of a multistage axial compressor. Our optimization process does not consider the adjustment of the angle of the stator blades, so it appears as a straight line on the pressurizer characteristic graph. The optimized blade can maintain a high flow rate even under higher back pressure and thus improve its surge margin.
Each blade row changed before and after optimization. Figure 14 shows the geometric comparison of R1, R2, R3, R4, S1, S2, and S5 before and after optimization. Among them, R4 and S2 underwent only one optimization, and the rest of the blade deformations that underwent only one optimization were similar to them. R1, R2, R3, S1, and S5 underwent two optimizations.
From Figures 14(d) and 14(f), it can be seen that the hub and tip shape of the optimized R4 becomes thinner and the blade-camber angle increases slightly; the leading edge part of the middle airfoil is basically unchanged, while the blade angle at the trailing edge increases slightly; thus, the blade-camber angle increases slightly; while the optimized S2 is the opposite of R4, after optimization, the hub and tip airfoils of S2 become thicker and the blade-camber angle becomes smaller, while the blade-camber angle of the middle airfoil decreases slightly.
From Figure 14(a), it can be seen that the hub bladecamber angle of R1 slightly increases, and the middle and tip airfoils do not change much with only a slight normal shift. From Figures 14(b) and 14(c), it can be seen that the hub blade-camber angle of the optimized R2 and R3 increases significantly and the leading edge blade angle decreases. The middle airfoil is similar to the hub, except that the degree is reduced.
For S1, the optimized blade-camber angle at the hub, middle, and tip regions was reduced, most significantly at the tip regions, as shown in Figure 14(e). For S5, the optimized blade-camber angles at the hub and the tip were significantly reduced, most significantly at the tip, and the blade-camber angle at the middle increased, as shown in Figure 14(g).     (4) Comparative Analysis in the Second Phase. The flow field before and after optimization should be compared and analyzed. However, it is difficult to highlight the details of the local flow field if all of them are displayed indiscriminately. Therefore, in order to better reveal the reason for the performance improvement of the optimized blade, typical blade rows are chosen for flow field analysis. The flow field changes in the four rotor blades are similar, and the flow field details of the typical blade row regions are selected to show and analyze the changes before and after optimization. Figure 15 shows a comparison of the relative Mach number distribution at the hub, midspan, and tip of R2, R3, and S5. Figure 16 illustrates the static pressure distribution of each row of blades at the hub, midspan, and tip. Combined with Figures 14-16, it can be seen that although the bladecamber angle of the hub and mid airfoil shapes of the optimized R2 and R3 becomes larger, the inverse pressure gradient near the trailing edge is reduced due to the change in load and slightly thinner thickness, resulting in a significant reduction in the airflow separation area at this location and a reduction in the airflow loss there. The situation of R3 is basically the same as that of R2. As seen from Figures 15(g) and 15(h), as the blade-camber angle of the optimized tandem S5 tip is reduced, the back pressure gradient near the trailing edge is reduced, and hence, the low-speed region caused by airflow mixing at the trailing edge of S5 is significantly reduced.
As a result, the flow in the channel is smoother, thereby reducing the airflow loss in the S5 outlet area.
The limit streamlines on the suction side of R1, R2, R3, and R4 are shown in Figure 17. From Figures 17(a)-17(d), it can be seen that for R1 and R2, the changes before and after optimization are not significant; for R3, there are large changes before and after optimization. The closed separation area at the tip of the optimized R3 is significantly reduced after optimization, thus reducing the corresponding loss. In addition, the separation line of the optimized R3 is shortened in the radial direction and thus reduce the radial separation area. At the same time, it can be observed that the optimized R3 will be attached again after separation, and then, a second separation line will be formed, reducing the separation area along the flow direction and thus reducing the corresponding separation loss. For R4, the optimized R4 has an additional reflux region at the tip, which increases the loss at the tip. For R4, original R4 has an additional reflux region at the tip, which increases the loss at the tip. Figure 18 shows the entropy distribution at midspan and tip section before and after optimization, from which it can be seen that the entropy value of the airflow is significantly lower after optimization (such as regions A, B, C, D, E, and F), which indicates a smoother airflow and a corresponding reduction in losses.

22
International Journal of Aerospace Engineering Based on the above analysis, it can be seen that after two phases of optimization, the R2 and R3 blade has a reduced reverse pressure gradient near the trailing edge at the hub and middle, resulting in a significant reduction in the separation area after the shock wave, and a separation bubble for reattachment is formed under these conditions (shown as optimized R3), thus reducing the corresponding loss and improving the efficiency. The reverse pressure gradient near the trailing edge of the S5 blade is also reduced due to the reduction of the blade-camber angle at the tip, narrowing the low-speed area near the trailing edge and downstream and reducing the corresponding airflow loss.

Conclusion
This paper proposes a new optimization method for highfidelity global optimization of multistage axial compressors based on the full-blade surface parameterization. The new method is applied to the optimization of the AL-31F four-stage low-pressure compressor aerodynamics. The following conclusions are drawn: (1) Compared with the traditional parametric method, the full-blade surface parametric method can effectively reduce the dimensionality and promote the solution of the HEB problem for a multistage axial flow compressor. At the same time, compared with the semiblade surface parametric method, the full-blade surface parametric method increases the degree of freedom of the blade angle change, therefore enhancing the probability of the existence of an optimal solution (2) After verifying on the benchmark function, the IABC algorithm is shown to have better optimization performance than the GA and ABC algorithms. The initialization method and exploration mechanism are applicable to the multipeak problem of multistage axial compressor optimization (3) The phased optimization strategy based on the "synchronous change in multirow blades" achieves a good balance between global exploration and local development. By introducing expert experience, the exploration of an invalid design space by optimization algorithms is avoided, and the optimization efficiency is greatly improved (4) The new optimization method based on the fullblade parameterization can be successfully applied to the AL-31F four-stage low-pressure compressor. When the flow and pressure ratio meet the constraints, the adiabatic efficiency and surge margin are increased by 0.67% and 3.1%, respectively Abbreviations HEB: High dimensionality, expensive cost, 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 mm: Millimeter rad: Radian TPR: Total pressure ratio eff: E fficiency SM: Surge margin DOE: Design of experiment Max_Run: The maximum number of tasks that can be concurrent at one time.
Nomenclature ξ: Chordwise direction η: Spanwise direction n p : Number of data points n s : Total number of sections T l j : Sum of the chord lengths of each segment of the jth section in the η direction T L i : Sum of the chord lengths of the ith section in the ξ direction l m : The mth chord length of the jth section in the η direction L n : The nth chord length of the ith section in the ξ direction The best food source location in the general information of the food source N e : Number of employed bees P: Probability of selecting a certain food source f ðX i Þ: Concentration of the ith food sources f ðX m Þ: Concentration of the mth food sources ðX j Neib Þ best : Food source position that has the largest concentration in the adjacent area ðV j Neib Þ best : The location of the new food source of onlooker bees dði, tÞ: Chebyshev distance between X i and X t md i : Average Chebyshev distance between X i and the entire onlooker bee population r: Radius of the neighborhood fitðÞ: Individual fitness of each bee P out : The total pressure ratios at the outlet P in : The total pressure ratios at the inlet T out : The total temperature at the outlet T in : The total temperature at the inlet.

Data Availability
The geometry and performance 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.