Multi-Objective Aerodynamic Design Optimisation Method of Fuel Cell Centrifugal Impeller Using Modified NSGA-II Algorithm

: This paper presents a modiﬁed NSGA-II algorithm based on the spatial density (SD) operator, combined with computer graphics-based surface parameterisation methods and computational ﬂuid dynamics (CFD) simulations. This was done to optimise the multi-objective aerodynamic design of a centrifugal impeller for a 100-kW vehicle-mounted fuel cell and improve the multi-conditions aerodynamic performance of the centrifugal impeller of the vehicle-mounted fuel cell (FC). The optimisation objectives are to maximise the isentropic efﬁciency of the rated and common operating conditions. The optimisation results showed that the efﬁciency of rated working conditions had an increase of 1.29%, mass ﬂow increase of 8.8%, pressure ratio increase of 0.74% and comprehensive margin increase of 6.2%. The efﬁciency of common working conditions had an increase of 1.2%, mass ﬂow increase of 9.1%, pressure ratio increase of 0.24% and comprehensive margin increase of 10%. The optimisation effect is obvious under the premise of satisfying the constraints, which proves the optimisation method’s engineering effectiveness and provides technical support and methodological research for the multi-objective aerodynamic design optimisation of centrifugal impellers for vehicle-mounted FCs. reduced. The decrease of the positive angle of attack also increases the ﬂow area, which increases the mass ﬂow after optimisation. The static surface pressure of the suction surface of the main blade and splitter blade is signiﬁcantly reduced, increasing the relative Mach number in the C, D and E regions. Meanwhile, the reduction of the area of the low-speed region reduces the reverse pressure gradient and improves the ﬂow performance of the downstream outlet F. Further, it can be seen from the comparison diagram of the static pressure distribution that the blade load increases after optimisation, which is conducive to the increase of work.


Introduction
New energy vehicles have been vigorously developed in recent years due to the problems of environmental pollution and oil shortages. Among them, fuel cell vehicles (FCVs) have a long-range, fast energy replenishment and high efficiency, which have more potential and are more attractive than lithium battery vehicles and traditional internal combustion engine vehicles. However, the high cost also limits their development. The fuel cell system (FCS) mainly consists of an electric stack, air supply system, hydrogen supply system, and water and thermal management systems. First, the electric motor drives a centrifugal impeller to pressurise the air. Then, the pressurised and humidified air is sent to the hydrogen fuel reactor for a chemical reaction. The final output of electricity provides energy for the electric motor and the whole vehicle. The compressor is the core component of the air supply system, accounting for 16.89% of the overall vehicle cost, second only to the fuel stack in the fuel cell (FC). Improving the aerodynamic performance and efficiency of the compressor can effectively save the overall cost.
Common types of air compressors are screw, scroll, roots, and centrifugal. A large mass flow is easy to achieve for the screw, roots and centrifugal compressors, and high pressure can be easily realized for the screw and centrifugal compressors [1]. The scroll, screw, and slide compressors are all in a type of positive displacement, which can boost the pressure through the shrinkage of inner chambers. In contrast, centrifugal compressors take in the gas and elevate the kinetic energy through a rotating impeller at high speed to boost the pressure of air. Centrifugal compressors have high efficiency, a compact structure 1.
The working environment is different. The FC reactor requires oil-free air, while the centrifugal impeller used in ICE is driven by oil-lubricated bearings, which can cause reactor contamination. 2.
The rated rotational speed is different. The ICE can make the centrifugal impeller reach more than 200,000 revolutions per minute (rpm) using high-temperature exhaust gas. Meanwhile, the exhaust gas of the FC cannot drive the centrifugal impeller to the required speed. The centrifugal impeller can only be driven by oil-free electric bearings, and its maximum speed can only reach 140,000 rpm. Further, the rated speed of 80,000 rpm-110,000 rpm is considered reasonable.

3.
It has a different working condition range. The specific speed of the centrifugal impeller in the FC is usually smaller than that of the ICE [5]. Further, their operating conditions are different. The operating range of the FC centrifugal impeller is narrow and usually runs along the surge boundary, while the operating range of the ICE centrifugal impeller is larger. If the ICE centrifugal impeller is applied to the FC system, there is a high risk of surge or blockage. In summary, FC centrifugal impellers need to be specifically designed and optimised for their own use.
Since the total pressure ratio of the FC air compressor can reach 0.2-0.4 Mpa, the driving power consumption can reach 25% of the output power of FCS [1,6]. Further, it cannot recover enough energy from the low-temperature exhaust gas to drive the highspeed motor like the ICE. Therefore, improving the aerodynamic performance of the centrifugal impeller can effectively reduce the system power consumption. According to the current study's results [7], the improvement of impeller efficiency can improve the material utilisation and comprehensive performance in FC to achieve further cost reduction. This work can be achieved through the design optimisation of the aerodynamic performance of the centrifugal impeller.
In FC gas supply systems, the higher the air pressure, the better the cost reduction. Studies have shown that the platinum content can be reduced from 0.36 g to 0.2 g/kW when the inlet pressure is increased from 1.5 to 3.0 atm. Second, the higher the air pressure, the higher the system's output power [8]. Higher air pressure is beneficial to increase the energy density of the stack, reducing the size and weight of the FC system and saving the manufacturing cost. Third, higher air pressure is beneficial to the water balance [4]. However, too high a stack inlet pressure consumes additional parasitic power, which reduces the efficiency of the whole machine. Therefore, when optimising the design of the aerodynamic performance of the FC centrifugal impeller, an appropriate increase in the pressure ratio is beneficial for the overall cost reduction. As well, the mass flow has a positive correlation with the FCS. The output power of the FCS will increase with the increase of the mass flow, which is beneficial to the efficiency of the FCS.
The traditional centrifugal impeller design method relies on the iterative design from 1D to 3D by experts. Since the 1990s, optimisation design methods have been introduced into the field of impeller aerodynamic design. They have been widely studied for their advantages of short design cycles, lack of dependence on expert experience, and breakthroughs from traditional design limitations [9,10]. Li et al. [9] presented a 3D multiobjective aerodynamic optimisation method by the multi-objective differential evolutionary algorithm, 3D blade parameterisation method, self-adaptive integration, and RANS. The optimisation results show that the isentropic efficiency and the total pressure ratio were increased by 3.06% and 1.26%, respectively. Ekradi et al. [10] used a genetic algorithm (GA), a 3D blade parameterisation method, a CFD solver and an artificial neural network, presenting a method for 3D optimisation of a transonic centrifugal impeller. The isentropic efficiency increased by 0.97%, the mass flow increased by 0.65%, and the total pressure ratio increased by 0.74%.
The optimisation algorithm is one of the core aspects of the centrifugal impeller aerodynamic design optimisation method, which determines the optimisation-seeking capability and convergence speed in the established design space. Optimisation algorithms for turbomachinery are generally classified into gradient methods, agent model methods (e.g., neural networks) and meta-heuristic algorithms (e.g., genetic algorithms). The gradient method is more efficient in the case of continuous and single-peaked variable space, but it can often only obtain the optimal local solution of the problem, not suitable for the optimisation search of complex multi-peaked problems. The agent model method can better reduce the computational cost, but it is very difficult to establish an agent model with sufficient accuracy since the dimensionality of decision variables increases. At the same time, the method is more difficult to deal with multi-objective multi-constraint problems. Although the meta-heuristic algorithms is computationally expensive, it has been widely used in impeller mechanical design optimisation due to it can obtain global optimisation solutions and has strong robustness [11]. It should be emphasised that the meta-heuristic algorithm is particularly suitable for the optimisation search of multi-objective multi-constraint problems due to its inherent concurrency property.
In recent years, meta-heuristic algorithms have been widely used to solve multiobjective optimization problems, and effectively applied to engineering practice, which are classified into two dominant classes: swarm intelligence and evolutionary. Swarm intelligence algorithms imitate the intelligence of swarms, flocks and herds of creatures while evolutionary algorithms imitate the mechanisms of evolution in nature. Mirjalili et al. [12] proposed a novel optimization algorithm named multi-objective salp swarm algorithm (MSSA) for solving optimization problems multiple objectives by imitate swarm behaviour of salps in oceans. MSSA had been tested on mathematical optimization functions and the results show that it is able to refine the initial random solutions effectively and convergence and is successfully applied to the multi-objective optimization problem of maximizing the efficiency and minimizing the cavitation of the propeller to achieve the optimal blade shape. Cao et al. [13] proposed a bi-swarm particle swarm optimizer with novel neighbourhood topology strategy (BPSO-NT) to solve the optimization of intermodal transportation planning in the field of logistics. Ten well-known benchmark functions are utilized to test the performance of BPSO-NT and the results show that the proposed BPSO-NT is superior to others PSO variants algorithm. In addition, BPSO-NT has achieved good results in solving practical engineering problems. Nadeem Shaukat et al. [14] presented an intelligent technique using modified genetic algorithms, combined with SuperMC code, to find the optimum loading pattern for the nuclear fuel, which is favourable to maximize effective multiplication factor and neutron flux densities inside the reactor core under design limits. In 2002, Deb et al. [15] proposed a modified version of NSGA, the well-known NSGA-II algorithm, which has a superior comprehensive performance compared with other evolutionary algorithms and is widely used in optimizing the centrifugal impellers. Duan et al. [16] used NSGA-II to optimise the conflict between pressure head and hydraulic efficiency. The tests proved that the optimised mini-pump increased the pressure head by about 24% and the hydraulic efficiency by 4.75% compared with the original design prototype, and the flow field analysis showed that the optimised model could alleviate the impeller wake and improve the pressure distribution and flow direction. Chen et al. [17] used NAGS-II to optimise the design of the automotive torque converter impeller. The results show that the zero-order polynomial of the Gaussian model is the best prediction for the stall torque ratio. They proposed an effective solution for the design of the hydraulic torque converter impeller geometry.
However, the NSGA-II algorithm is excessively focused on convergence, leading to controversy over its performance in the population diversity. Wang et al. [18] pointed out that the defect of NSGA-II is the selecting principle of crowding distance which is based on the cubic distance to the adjacent points. For instance, in the ZDT problems, with continuous evolution, more and more points would cluster to a small range. The diversity of individuals is lost which leads to a local optimum. Huang et al. [19] indicated that the crowding distance operator in NSGA-II can guarantee the diversity along the non-dominated front, but it may lose the horizontal diversity when wiping off excessive individuals. Due to the wide range of FC operating conditions and the higher requirements for aerodynamic performance under different operating conditions, the multi-objective problems using the classical NSGA-II algorithm easily fall into the local optimal solution.
In order to improve the population diversity of solution sets for multiple operating conditions performance of FC centrifugal impeller, a modified NSGA-II optimisation algorithm (i.e., DNSGA) is proposed, which has a novel classification method for nondominated ranking and crowding distance operator (i.e., SD operator with weights).

Algorithm Framework
There are two aspects of shortcomings for the NSGA-II algorithm. On the one hand, it is easy to fall into local optimal solutions and into premature convergence due to its focus on convergence. On the other hand, when some individuals have large differences in decision variables while their objective function values are approximately equal, the adoption of NSGA-II for the crowdedness ranking may remove potentially excellent individuals, which is unfavourable to the diversity of the population. The DNSGA algorithm is proposed to remedy the shortcomings, which essentially increases the optional range of individual decision space by improving the non-dominated ranking and spatial density ranking. This is conducive to the inheritance of excellent individuals to the next generation, thus improving the diversity of individuals in the population.
The next generation of DNSGA is selected in three main steps. First, the parents and offspring are divided into three groups (M1, M2 and M3) according to the non-dominance result. The number of individuals in M1 does not exceed 50% of the population size. The number of individuals in M2 is not smaller than the population size. The remaining individuals are divided into M3. Second, M3 is deleted. M1 and M2 are combined into group Q. Finally, the individuals with the lowest spatial density in M2 are sequentially deleted according to the SD until the number of individuals in Q is equal to the population size.
The main workflow of the algorithm is shown in Figure 1, and the specific steps are as follows: (1) Randomly generate an initial population P 0 to start the evolution; (2) Perform the binary tournament selection, simulated binary crossover and polynomial variation on P 0 to generate new offspring P * g ; (3) Evaluate the fitness function so that multi-objective values are obtained for each individual; (4) Sort all individuals in P g and P * g into F 1 , F 2 , F 3 . . . . . . F n based on non-dominance. (5) Divide F 1 , F 2 , F 3 . . . . . . F n into the three groups of M1, M2 and M3; (6) Generate the next population P g+1 according to the spatial density operator ordering; (7) Return to step 2 until the end of the maximum number of iterations.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 5 of 22 (6) Generate the next population 1 g P according to the spatial density operator ordering; (7) Return to step 2 until the end of the maximum number of iterations.

Spatial Density Ordering
The spatial density expressions for the two individuals of   , ,..., n Y x x x  with n-dimensional decision variables is expressed as Formula (1).
Empirically, the influence of the decision variables of different dimensions on the aerodynamic performance of the centrifugal impeller varies widely. Therefore, different weights are assigned to different decision variables. The mathematical expression of spatial density with weights is Equation (2).
where i w is the weight of the ith dimension variable.
In this paper, a spatial density operator ranking with weights is used to remove the redundant individuals in M2, and the procedure is as follows: (1) M1 and M2 are combined into group Q.
(2) The two individuals with the lowest spatial density in Q are found, and at least one belongs to M2. (3) If one individual belongs to M1 and the other belongs to M2. The individual belonging to M2 is directly removed from Q. If both individuals belong to M2, the one who has a lower spatial density with other individuals in Q is deleted. (4) Return (2) and (3) in a loop until the total number of individuals in Q reaches the population size.

Algorithm Validation
Five classical multi-objective algorithms, including VEGA, SPEA2, EMOEA, NSGAII and DNSGA, are run 30 times independently under two classical test functions (Table 1) from the ZDT problems proposed by Ziztler et al. [20] to test the performance of the presented DNSGA.
Empirically, the influence of the decision variables of different dimensions on the aerodynamic performance of the centrifugal impeller varies widely. Therefore, different weights are assigned to different decision variables. The mathematical expression of spatial density with weights is Equation (2).
where w i is the weight of the ith dimension variable. In this paper, a spatial density operator ranking with weights is used to remove the redundant individuals in M2, and the procedure is as follows: (1) M1 and M2 are combined into group Q.
(2) The two individuals with the lowest spatial density in Q are found, and at least one belongs to M2. (3) If one individual belongs to M1 and the other belongs to M2. The individual belonging to M2 is directly removed from Q. If both individuals belong to M2, the one who has a lower spatial density with other individuals in Q is deleted. (4) Return (2) and (3) in a loop until the total number of individuals in Q reaches the population size.

Algorithm Validation
Five classical multi-objective algorithms, including VEGA, SPEA2, EMOEA, NSGAII and DNSGA, are run 30 times independently under two classical test functions (Table 1) from the ZDT problems proposed by Ziztler et al. [20] to test the performance of the presented DNSGA.
The hypervolume ratio (HVR) and generational distance (GD) are recorded, and statistical results are shown in Table 2 and the box plots are shown in Figures 2 and 3. On the one hand, the HVR metric measures the volume covered by the obtained PF for problems where all objectives are to be minimised. A larger value of the metric indicates better diversity. On the other hand, the GD metric indicates how far the tested PF is from the real PF. A low GD value reflects good convergence performance.

Multi-Objective Aerodynamic Design Optimisation of Centrifugal Impeller
Based on the DNSGA algorithm, a new multi-objective optimisation method for the centrifugal impeller is constructed by combining the MDOF method and the CFD technique. This optimisation method is used to carry out multi-objective optimisation for the centrifugal air compressor of the vehicle-mounted FC to obtain a global optimisation scheme at different rotational speeds.

Numerical Method
CFD techniques have been widely used to study the internal flow of centrifugal impellers [21][22][23][24]. In this paper, the steady-state flow field of a centrifugal impeller is obtained by solving the three-dimensional (3D) steady-state Reynolds averaged Navier-Stokes equation using NUMECA's Fine Turbo solver. The turbulence model adopts the S-A model, which is a simple one-equation model that can capture turbulence in the endwall constrained case. Meanwhile, the S-A model is insensitive to the numerical errors

Multi-Objective Aerodynamic Design Optimisation of Centrifugal Impeller
Based on the DNSGA algorithm, a new multi-objective optimisation method for the centrifugal impeller is constructed by combining the MDOF method and the CFD technique. This optimisation method is used to carry out multi-objective optimisation for the centrifugal air compressor of the vehicle-mounted FC to obtain a global optimisation scheme at different rotational speeds.

Numerical Method
CFD techniques have been widely used to study the internal flow of centrifugal impellers [21][22][23][24]. In this paper, the steady-state flow field of a centrifugal impeller is obtained by solving the three-dimensional (3D) steady-state Reynolds averaged Navier-Stokes equation using NUMECA's Fine Turbo solver. The turbulence model adopts the S-A model, which is a simple one-equation model that can capture turbulence in the endwall constrained case. Meanwhile, the S-A model is insensitive to the numerical errors As shown in Table 2, Figures 2 and 3, for test functions the ZDT2 and the ZDT3, the mean values of the HVR metrics for DNSGA are significantly higher than that of other classical multi-objective algorithms, and the standard deviations of the HVR metrics for DNSGA are also relatively lower than other multi-objective algorithms, indicating that the population diversity performance of DNSGA is superior to others evolution algorithms.
For the ZDT2, the mean value and standard deviation of the GD metrics for DNSGA do not differ much from the other algorithms. However, for the ZDT3, the mean value of the GD metrics for DNSGA lower than VEGA and higher than other algorithms, and the standard deviation of the GD metrics for DNSGA do not differ much from the other algorithms.
To solve the defects of NAGAII algorithm, the proposed DNSGA algorithm trades off the population diversity and the convergence. By the comparison of HVR and GD metric, although the convergence of DNSGA is slightly inferior to that of other classical multiobjective optimization algorithms, DNSGA is superior to other algorithms in population diversity. For the multi-objective aerodynamic optimization problem of complex operating conditions for FC, we put more emphasis on the performance of population diversity.

Multi-Objective Aerodynamic Design Optimisation of Centrifugal Impeller
Based on the DNSGA algorithm, a new multi-objective optimisation method for the centrifugal impeller is constructed by combining the MDOF method and the CFD technique. This optimisation method is used to carry out multi-objective optimisation for the centrifugal air compressor of the vehicle-mounted FC to obtain a global optimisation scheme at different rotational speeds.

Numerical Method
CFD techniques have been widely used to study the internal flow of centrifugal impellers [21][22][23][24]. In this paper, the steady-state flow field of a centrifugal impeller is obtained by solving the three-dimensional (3D) steady-state Reynolds averaged Navier-Stokes equation using NUMECA's Fine Turbo solver. The turbulence model adopts the S-A model, which is a simple one-equation model that can capture turbulence in the end-wall constrained case. Meanwhile, the S-A model is insensitive to the numerical errors caused by mesh roughness and has been widely used in impeller machinery. The fourth-order explicit Runge-Kutta model is used for the time marching. A finite volume central difference scheme with second-order and fourth-order artificial viscosity terms is used to control the pseudo numerical oscillation in spatial discretisation. The convergence is accelerated by using multigrid, local time step, and implicit residual. The blade surface grid adopts HOH topology, the blade tip and hub gap grid adopt 4OH topology, and the thickness of the first grid layer near the wall is 1 × 10 −6 m, Y + ≤ 1.

Mesh Generation and Independence Verification
The grid independence of a single channel of the impeller at a rated speed is verified to ensure the grid quality of flow field calculation. The total grid number of the main blade and splitter blade is divided into 0.3 million, 0.66 million, 1.03 million and 1.5 million, respectively. The results are shown in Figure 4. When the grid number reaches 1.03 million, the grid independence requirement can be satisfied. Therefore, the subsequent optimisation process selects 1.03 million as the grid template for calculation. The boundary conditions for the numerical calculation are set as follows: the total compressor inlet temperature is 293 K, the total pressure is 101,325 Pa, the inlet direction is axial, and the outlet is the average static pressure. The convergence point before the first divergence point is the near-surge point by gradually increasing the backpressure and advancing the calculation from the blockage point to the near-surge point. The blade surface and end wall are set as no-slip boundary conditions. caused by mesh roughness and has been widely used in impeller machinery. The fourthorder explicit Runge-Kutta model is used for the time marching. A finite volume central difference scheme with second-order and fourth-order artificial viscosity terms is used to control the pseudo numerical oscillation in spatial discretisation. The convergence is accelerated by using multigrid, local time step, and implicit residual. The blade surface grid adopts HOH topology, the blade tip and hub gap grid adopt 4OH topology, and the thickness of the first grid layer near the wall is

Mesh Generation and Independence Verification
The grid independence of a single channel of the impeller at a rated speed is verified to ensure the grid quality of flow field calculation. The total grid number of the main blade and splitter blade is divided into 0.3 million, 0.66 million, 1.03 million and 1.5 million, respectively. The results are shown in Figure 4. When the grid number reaches 1.03 million, the grid independence requirement can be satisfied. Therefore, the subsequent optimisation process selects 1.03 million as the grid template for calculation. The boundary conditions for the numerical calculation are set as follows: the total compressor inlet temperature is 293 K, the total pressure is 101,325 Pa, the inlet direction is axial, and the outlet is the average static pressure. The convergence point before the first divergence point is the near-surge point by gradually increasing the backpressure and advancing the calculation from the blockage point to the near-surge point. The blade surface and end wall are set as no-slip boundary conditions.

Numerical Method Validation
The validity of the numerical simulation was verified using the Krain impeller as the test model [25,26], which was verified by our team [5,27]. The numerical simulation results matched well with the Krain experimental data concerning the trend of the performance curves, indicating that the numerical method adopted in this paper can provide an accurate qualitative analysis of the aerodynamic performance changes before and after optimisation.

Optimisation Object
In this paper, the performance of FCVs single-stage high-speed centrifugal impeller with a rated output of 100 kW is optimised, and its CAD model is shown in Figure 5, and the aerodynamic and geometric parameters are shown in Table 3.

Numerical Method Validation
The validity of the numerical simulation was verified using the Krain impeller as the test model [25,26], which was verified by our team [5,27]. The numerical simulation results matched well with the Krain experimental data concerning the trend of the performance curves, indicating that the numerical method adopted in this paper can provide an accurate qualitative analysis of the aerodynamic performance changes before and after optimisation.

Optimisation Object
In this paper, the performance of FCVs single-stage high-speed centrifugal impeller with a rated output of 100 kW is optimised, and its CAD model is shown in Figure 5, and the aerodynamic and geometric parameters are shown in Table 3.

Optimisation Objectives and Constraints
Centrifugal impeller optimised design conditions are usually selected at rated speed for geometric parameter improvement. Further, non-design conditions are often selected at other mass flow at the same speed, which can meet the requirements for industrial compressors working in a single environment. However, the working speed range of the vehicle-mounted centrifugal impeller is wide, and the performance of the rated speed can be improved by optimising the rated speed only. However, the performance of the common operating conditions cannot be guaranteed. The isentropic efficiency of the aerodynamic performance of the centrifugal impeller of FC at the rated operating point and common operating point is optimised by multi-objective optimisation in this paper to improve them at the same time. Take 100,000 rpm, mass flow of 118.32 g/s and pressure ratio of 2.7 as the rated working condition. The selection of the common condition point is based on the U.S. Department of Energy (DOE) technical requirements for the same class of automotive FC power systems. A total of 70% of the rated speed of the centrifugal impeller, 77.36 g/s mass flow and 1.7 total pressure ratio are set as the common operating point. The original performance curve under two working conditions is shown in Figure 6.

Optimisation Objectives and Constraints
Centrifugal impeller optimised design conditions are usually selected at rated speed for geometric parameter improvement. Further, non-design conditions are often selected at other mass flow at the same speed, which can meet the requirements for industrial compressors working in a single environment. However, the working speed range of the vehicle-mounted centrifugal impeller is wide, and the performance of the rated speed can be improved by optimising the rated speed only. However, the performance of the common operating conditions cannot be guaranteed. The isentropic efficiency of the aerodynamic performance of the centrifugal impeller of FC at the rated operating point and common operating point is optimised by multi-objective optimisation in this paper to improve them at the same time. Take 100,000 rpm, mass flow of 118.32 g/s and pressure ratio of 2.7 as the rated working condition. The selection of the common condition point is based on the U.S. Department of Energy (DOE) technical requirements for the same class of automotive FC power systems. A total of 70% of the rated speed of the centrifugal impeller, 77.36 g/s mass flow and 1.7 total pressure ratio are set as the common operating point. The original performance curve under two working conditions is shown in Figure 6. According to the above description, maximising the isentropic efficiency at two operating points is the optimisation objective, with the optimised mass flow not decreasing and the optimised total pressure ratio in-between 1 and 1.01 times the original total pressure ratio as the constraint. The mathematical expression is as follows: Objective function: max

Parameterisation Method
The development of impeller surface parameterisation methods has gone through from half-blade full parameterisation [28,29] to full-blade surface parameterisation [30], and then to the latest MDOF [27,31]. These methods are based on computer graphics by changing the blade suction surface and pressure surface control points to generate new blades directly. MDOF surface parametrisation has obvious advantages over other parametrisation methods in the process of blade surface generation, with the advantages of a According to the above description, maximising the isentropic efficiency at two operating points is the optimisation objective, with the optimised mass flow not decreasing and the optimised total pressure ratio in-between 1 and 1.01 times the original total pressure ratio as the constraint. The mathematical expression is as follows: Objective function : Constraint function: E f f d_70k is the isentropic efficiency at the design point of 70,000 rpm, E f f d_100k is the isentropic efficiency at the design point of 100,000 rpm, and Mass d_ori_70k and Mass o_70k are the mass flow at the original impeller design point and the optimised impeller design point at 70,000 rpm, respectively. Mass d_ori_100k and Mass o_100k are the mass flow at the original impeller design point and the optimised impeller design point at 100,000 rpm, respectively. T pr d_ori_70k and T pr o_70k are the total pressure ratios at the original impeller design point and the optimised impeller design point at 70,000 rpm, respectively. T pr d_ori_100k and T pr o_100k are the total pressure ratios at the original impeller design point and the optimised impeller design point at 100,000 rpm, respectively. x i is the optimisation variable, b x U i is the upper bound of the optimisation variable, and x L i is the lower bound of the optimisation variable.

Parameterisation Method
The development of impeller surface parameterisation methods has gone through from half-blade full parameterisation [28,29] to full-blade surface parameterisation [30], and then to the latest MDOF [27,31]. These methods are based on computer graphics by changing the blade suction surface and pressure surface control points to generate new blades directly. MDOF surface parametrisation has obvious advantages over other parametrisation methods in the process of blade surface generation, with the advantages of a low number of optimisation variables, smooth blade, flexible reshaping, and guaranteed blade strength. In this paper, the MDOF parameterisation method is used to parameterise the surfaces of the main blade and the splitter blade separately using two 6 × 3 order Bezier surfaces. As shown in Figure 7, each surface has seven control vertexes in the ξ direction (0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.0) and four points in the η direction (0, 0.4, 0.7, 1.0). It is necessary to keep ξ 1 and ξ 2 fixed to ensure the first-order continuity of the leading edge. At the same time, the variation of the purple control vertex of the suction surface is consistent with the variation of the green control vertex of the corresponding pressure surface during the reshaping process to ensure that the blade thickness does not get thinner. Therefore, the values of the normal coordinates of the green control vertexes on the pressure surface are the optimisation variables, and the total number is only (5 × 4) × 2 = 40. By trial and error, the variable ranges for control vertexes are set as shown in Table 4. low number of optimisation variables, smooth blade, flexible reshaping, and guaranteed blade strength. In this paper, the MDOF parameterisation method is used to parameterise the surfaces of the main blade and the splitter blade separately using two 6 × 3 order Bezier surfaces. As shown in Figure 7, each surface has seven control vertexes in the ξ direction (0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.0) and four points in the η direction (0, 0.4, 0.7, 1.0). It is necessary to keep ξ1 and ξ2 fixed to ensure the first-order continuity of the leading edge. At the same time, the variation of the purple control vertex of the suction surface is consistent with the variation of the green control vertex of the corresponding pressure surface during the reshaping process to ensure that the blade thickness does not get thinner. Therefore, the values of the normal coordinates of the green control vertexes on the pressure surface are the optimisation variables, and the total number is only (5 × 4) × 2 = 40. By trial and error, the variable ranges for control vertexes are set as shown in Table 4.

Optimisation Process and Algorithm Parameter Setting
The detailed steps of the optimisation process for this experiment are shown in Figure 8, and the optimisation parameters are set as shown in Table 5.
In Step 1, the population individuals were initialised, and the fitness values of the parent individuals were obtained using the MDOF parameterisation method and CFD numerical calculations. The more detailed processes of the CFD implemented into the optimisation model are as follows: Firstly, the mesh template is obtained through the mesh

Optimisation Process and Algorithm Parameter Setting
The detailed steps of the optimisation process for this experiment are shown in Figure 8, and the optimisation parameters are set as shown in Table 5.
The parameters of the weight coefficients of the spatial operator decision variables are set as shown in Table 6. The setting value of the weighting coefficient is based on engineering experience. The main blade of the centrifugal impeller, the middle of the blade, the trailing edge, the root and the tip of the blade have a greater impact on the aerodynamic performance, so these positions are set to a larger weighting coefficient.
In Step 4, the process returned to Step 2 and continued until the number of iterations was satisfied.    In Step 1, the population individuals were initialised, and the fitness values of the parent individuals were obtained using the MDOF parameterisation method and CFD numerical calculations. The more detailed processes of the CFD implemented into the optimisation model are as follows: Firstly, the mesh template is obtained through the mesh division and independence verification of the original impeller; Secondly, the new impeller geometry is generated by using MDOF parameterization and optimization algorithm, and then the geometric mesh of the new impeller can be obtained through the original mesh template division. Finally, the NUMECA-FINE module is called to perform aerodynamic calculations on the new mesh to achieve the aerodynamic performance of the centrifugal impeller, including flow rate, total pressure ratio and efficiency. In addition, the mass flow and total pressure ratio should satisfy the constraints, and the efficiency is searched for the maximum value in the optimization process.
In Step 2, the offspring were generated based on selection, crossover, and variation. The fitness values of the calculated offspring individuals were obtained using the MDOF parameterisation method and CFD numerical calculations.
In Step 3, the parents and offspring were merged, and the next generation of parent individuals was selected by non-dominated sorting and spatial density operator sorting. The parameters of the weight coefficients of the spatial operator decision variables are set as shown in Table 6. The setting value of the weighting coefficient is based on engineering experience. The main blade of the centrifugal impeller, the middle of the blade, the trailing edge, the root and the tip of the blade have a greater impact on the aerodynamic performance, so these positions are set to a larger weighting coefficient. In Step 4, the process returned to Step 2 and continued until the number of iterations was satisfied.

Comparison of Blade Geometry
Based on engineering experience and aerodynamic requirements for centrifugal impellers, the best individual is selected from the PF solution (i.e., the optimal geometry blade). Figure 9 shows the geometric deformation of the blade before and after optimisation. Figure 10 shows the deformation contour of the blade surface before and after optimisation. Figure 11 shows the comparison of the hub, middle and tip geometry before and after optimisation.

Comparison of Blade Geometry
Based on engineering experience and aerodynamic requirements for centrifugal impellers, the best individual is selected from the PF solution (i.e., the optimal geometry blade). Figure 9 shows the geometric deformation of the blade before and after optimisation. Figure 10 shows the deformation contour of the blade surface before and after optimisation. Figure 11 shows the comparison of the hub, middle and tip geometry before and after optimisation.

Comparison of Blade Geometry
Based on engineering experience and aerodynamic requirements for centrifugal impellers, the best individual is selected from the PF solution (i.e., the optimal geometry blade). Figure 9 shows the geometric deformation of the blade before and after optimisation. Figure 10 shows the deformation contour of the blade surface before and after optimisation. Figure 11 shows the comparison of the hub, middle and tip geometry before and after optimisation.   For the main blade along the span-wise direction, the deformations at the tip and hub are larger, while the middle section is essentially unchanged. Along the chord length direction, the change is smaller at the leading edge. The tip and middle section at the middle position bend toward the suction surface, and the tip bends toward the pressure surface. At the trailing edge, the tip bends toward the suction surface, the middle section basically remains unchanged, and the hub bends toward the pressure surface.
For the splitter blade along the span-wise direction, the deformations at the tip and hub are larger, while the middle section position changes less. Along the chord length direction, the leading edge is basically unchanged. The middle position of the tip bends to the suction surface, the middle section is basically unchanged, and the hub bends to the pressure surface. The tip of the trailing edge bends to the suction surface, the middle section is unchanged, and the hub bends to the pressure surface.
In general, the lower and middle regions of the flow path between the pressure surface of the main blade and the suction surface of the splitter blade become wider, and the top becomes narrower after optimisation. Meanwhile, the lower and middle regions of the flow path between the suction surface of the main blade and the pressure surface of the splitter blade become narrower, and the top becomes wider. The change of blade geometry and flow channel area will inevitably cause aerodynamic performance and flow field structure changes. Tables 7 and 8 show the comparison of the aerodynamic performance before and after the optimisation at the design point of rated condition and the design point of a common condition, respectively. Further, Figure 12 shows the comparison of the performance curves before and after the optimisation of the rated condition and common condition. After optimisation, the efficiency at the design point of the rated operating conditions increased by 1.29%, the mass flow increased by 8.8%, the total pressure ratio increased by 0.74%, and the comprehensive margin increased by 6.2%. The efficiency of the common operating conditions improved by 1.2%, the mass flow improved by 9.1%, the pressure ratio improved by 0.24%, and the comprehensive margin improved by 10%. The Figure 11. Comparison of the hub, middle and tip geometry before and after optimisation.

Comparison of Multi-Conditions Aerodynamic Performance before and after Optimisation
For the main blade along the span-wise direction, the deformations at the tip and hub are larger, while the middle section is essentially unchanged. Along the chord length direction, the change is smaller at the leading edge. The tip and middle section at the middle position bend toward the suction surface, and the tip bends toward the pressure surface. At the trailing edge, the tip bends toward the suction surface, the middle section basically remains unchanged, and the hub bends toward the pressure surface.
For the splitter blade along the span-wise direction, the deformations at the tip and hub are larger, while the middle section position changes less. Along the chord length direction, the leading edge is basically unchanged. The middle position of the tip bends to the suction surface, the middle section is basically unchanged, and the hub bends to the pressure surface. The tip of the trailing edge bends to the suction surface, the middle section is unchanged, and the hub bends to the pressure surface.
In general, the lower and middle regions of the flow path between the pressure surface of the main blade and the suction surface of the splitter blade become wider, and the top becomes narrower after optimisation. Meanwhile, the lower and middle regions of the flow path between the suction surface of the main blade and the pressure surface of the splitter blade become narrower, and the top becomes wider. The change of blade geometry and flow channel area will inevitably cause aerodynamic performance and flow field structure changes. Tables 7 and 8 show the comparison of the aerodynamic performance before and after the optimisation at the design point of rated condition and the design point of a common condition, respectively. Further, Figure 12 shows the comparison of the performance curves before and after the optimisation of the rated condition and common condition. After optimisation, the efficiency at the design point of the rated operating conditions increased by 1.29%, the mass flow increased by 8.8%, the total pressure ratio increased by 0.74%, and the comprehensive margin increased by 6.2%. The efficiency of the common operating conditions improved by 1.2%, the mass flow improved by 9.1%, the pressure ratio improved by 0.24%, and the comprehensive margin improved by 10%. The centrifugal impeller aerodynamic performance has been significantly improved, and the performance under different working conditions can be improved simultaneously.  centrifugal impeller aerodynamic performance has been significantly improved, and the performance under different working conditions can be improved simultaneously.    Figure 13 shows the entropy contour of the S3 section before and after the optimisation of the rated operating condition. After optimisation, the high entropy areas in channels I and II, A, B, C and D are reduced, and the losses are significantly reduced. The reason for the decrease of inlet loss is that the shock wave at the blade inlet is weakened due to the better matching of the inlet attack angle. Figure 14 shows the entropy contour of the S3 section before and after the optimisation of the common operating condition. The entropy value of A, B, C, D, E and F is obviously reduced after optimisation, given that the incidence angle and blade geometry are better matched, making the acceleration section of the blade inlet shorter and acceleration slower. This is beneficial to reduce the loss.  Figure 13 shows the entropy contour of the S3 section before and after the optimisation of the rated operating condition. After optimisation, the high entropy areas in channels I and II, A, B, C and D are reduced, and the losses are significantly reduced. The reason for the decrease of inlet loss is that the shock wave at the blade inlet is weakened due to the better matching of the inlet attack angle. Figure 14 shows the entropy contour of the S3 section before and after the optimisation of the common operating condition. The entropy value of A, B, C, D, E and F is obviously reduced after optimisation, given that the incidence angle and blade geometry are better matched, making the acceleration section of the blade inlet shorter and acceleration slower. This is beneficial to reduce the loss. As can be seen from Figures 10-12, the top of the flow channel between the pressure surface of the main blade and the suction surface of the splitter blade becomes narrower. Further, the back-bending angle is increased after optimisation, which reduces the pressure expansion angle and reduces the return loss. This is the main reason for the reduction of the top clearance and loss at the outlet of the rated and common operating conditions and the increase of efficiency. This is the main reason for the reduction of tip and outlet losses.

Flow Field Analysis
Appl. Sci. 2021, 11, x FOR PEER REVIEW 16 of 22 As can be seen from Figures 10-12, the top of the flow channel between the pressure surface of the main blade and the suction surface of the splitter blade becomes narrower. Further, the back-bending angle is increased after optimisation, which reduces the pressure expansion angle and reduces the return loss. This is the main reason for the reduction of the top clearance and loss at the outlet of the rated and common operating conditions and the increase of efficiency. This is the main reason for the reduction of tip and outlet losses.
(a) original (b) optimised Figure 13. Entropy contour of the S3 cross-section before and after optimisation (rated operating condition).
(a) original (b) optimised Figure 14. Entropy contour of the S3 cross-section before and after optimisation (common operating condition). Figure 15 shows the isentropic efficiency changes along the span-wise at the downstream outlet for different operating conditions before and after optimisation. It can be seen that the most significant efficiency improvement is roughly in the area of 80% in the span-wise direction. As such, the flow field structure at the 80% height (H) of the impeller was chosen as the focus of the analysis.  As can be seen from Figures 10-12, the top of the flow channel between the pressure surface of the main blade and the suction surface of the splitter blade becomes narrower. Further, the back-bending angle is increased after optimisation, which reduces the pressure expansion angle and reduces the return loss. This is the main reason for the reduction of the top clearance and loss at the outlet of the rated and common operating conditions and the increase of efficiency. This is the main reason for the reduction of tip and outlet losses.
(a) original (b) optimised Figure 13. Entropy contour of the S3 cross-section before and after optimisation (rated operating condition).
(a) original (b) optimised Figure 14. Entropy contour of the S3 cross-section before and after optimisation (common operating condition). Figure 15 shows the isentropic efficiency changes along the span-wise at the downstream outlet for different operating conditions before and after optimisation. It can be seen that the most significant efficiency improvement is roughly in the area of 80% in the span-wise direction. As such, the flow field structure at the 80% height (H) of the impeller was chosen as the focus of the analysis.  Figure 15 shows the isentropic efficiency changes along the span-wise at the downstream outlet for different operating conditions before and after optimisation. It can be seen that the most significant efficiency improvement is roughly in the area of 80% in the span-wise direction. As such, the flow field structure at the 80% height (H) of the impeller was chosen as the focus of the analysis. Figure 16 shows the static pressure distribution on the blade surface at the 80% height of the impeller before and after the optimisation of the rated operating condition. Figure 17 shows the relative Mach number contour on the B2B surface at the 80% height of the impeller before and after the optimisation of the rated operating condition. After optimisation, the positive angle of the attack at the inlet is reduced, the matching of airflow is improved, the area and maximum value of the high relative Mach number of the front edge A and B of the main blade suction are reduced, and the shock loss at the inlet is reduced. At the same time, the loss of the interaction between the shockwave and the boundary layer is reduced, which makes the flow separation in the latter half of the suction surface move backwards. Further, the separation area is reduced, and the separation loss is reduced. The decrease of the positive angle of attack also increases the flow area, which increases the mass flow after optimisation. The static surface pressure of the suction surface of the main blade and splitter blade is significantly reduced, increasing the relative Mach number in the C, D and E regions. Meanwhile, the reduction of the area of the low-speed region reduces the reverse pressure gradient and improves the flow performance of the downstream outlet F. Further, it can be seen from the comparison diagram of the static pressure distribution that the blade load increases after optimisation, which is conducive to the increase of work.  Figure 16 shows the static pressure distribution on the blade surface at the 80% height of the impeller before and after the optimisation of the rated operating condition. Figure  17 shows the relative Mach number contour on the B2B surface at the 80% height of the impeller before and after the optimisation of the rated operating condition. After optimisation, the positive angle of the attack at the inlet is reduced, the matching of airflow is improved, the area and maximum value of the high relative Mach number of the front edge A and B of the main blade suction are reduced, and the shock loss at the inlet is reduced. At the same time, the loss of the interaction between the shockwave and the boundary layer is reduced, which makes the flow separation in the latter half of the suction surface move backwards. Further, the separation area is reduced, and the separation loss is reduced. The decrease of the positive angle of attack also increases the flow area, which increases the mass flow after optimisation. The static surface pressure of the suction surface of the main blade and splitter blade is significantly reduced, increasing the relative Mach number in the C, D and E regions. Meanwhile, the reduction of the area of the lowspeed region reduces the reverse pressure gradient and improves the flow performance of the downstream outlet F. Further, it can be seen from the comparison diagram of the static pressure distribution that the blade load increases after optimisation, which is conducive to the increase of work.   Figure 16 shows the static pressure distribution on the blade surface at the 80% height of the impeller before and after the optimisation of the rated operating condition. Figure  17 shows the relative Mach number contour on the B2B surface at the 80% height of the impeller before and after the optimisation of the rated operating condition. After optimisation, the positive angle of the attack at the inlet is reduced, the matching of airflow is improved, the area and maximum value of the high relative Mach number of the front edge A and B of the main blade suction are reduced, and the shock loss at the inlet is reduced. At the same time, the loss of the interaction between the shockwave and the boundary layer is reduced, which makes the flow separation in the latter half of the suction surface move backwards. Further, the separation area is reduced, and the separation loss is reduced. The decrease of the positive angle of attack also increases the flow area, which increases the mass flow after optimisation. The static surface pressure of the suction surface of the main blade and splitter blade is significantly reduced, increasing the relative Mach number in the C, D and E regions. Meanwhile, the reduction of the area of the lowspeed region reduces the reverse pressure gradient and improves the flow performance of the downstream outlet F. Further, it can be seen from the comparison diagram of the static pressure distribution that the blade load increases after optimisation, which is conducive to the increase of work.    Figure 18 shows the static pressure distribution on the blade surface at the 80% height of the impeller before and after optimisation of the common operating condition. Figure  19 shows the relative Mach number contour on the B2B surface at the 80% height of the impeller before and after optimisation of the common operating condition. After optimisation, the matching of inlet flow angle is improved, the acceleration section is shortened, the area and maximum value of the high relative Mach number of suction leading-edge A and B are reduced, and the shock loss and separation loss at the inlet is reduced. Moreover, the decrease of the positive angle of attack also increases the flow area, which increases the mass flow after optimisation. The static pressure on the suction surface of the main blade and splitter blade is significantly reduced, increasing the relative Mach number in the areas of C, D and E, effectively reducing the low-speed area of the channel, reducing the reverse pressure gradient, delaying the separation position, and improving the aerodynamic performance of the downstream outlet F.   Figure 18 shows the static pressure distribution on the blade surface at the 80% height of the impeller before and after optimisation of the common operating condition. Figure 19 shows the relative Mach number contour on the B2B surface at the 80% height of the impeller before and after optimisation of the common operating condition. After optimisation, the matching of inlet flow angle is improved, the acceleration section is shortened, the area and maximum value of the high relative Mach number of suction leading-edge A and B are reduced, and the shock loss and separation loss at the inlet is reduced. Moreover, the decrease of the positive angle of attack also increases the flow area, which increases the mass flow after optimisation. The static pressure on the suction surface of the main blade and splitter blade is significantly reduced, increasing the relative Mach number in the areas of C, D and E, effectively reducing the low-speed area of the channel, reducing the reverse pressure gradient, delaying the separation position, and improving the aerodynamic performance of the downstream outlet F.  Figure 18 shows the static pressure distribution on the blade surface at the 80% height of the impeller before and after optimisation of the common operating condition. Figure  19 shows the relative Mach number contour on the B2B surface at the 80% height of the impeller before and after optimisation of the common operating condition. After optimisation, the matching of inlet flow angle is improved, the acceleration section is shortened, the area and maximum value of the high relative Mach number of suction leading-edge A and B are reduced, and the shock loss and separation loss at the inlet is reduced. Moreover, the decrease of the positive angle of attack also increases the flow area, which increases the mass flow after optimisation. The static pressure on the suction surface of the main blade and splitter blade is significantly reduced, increasing the relative Mach number in the areas of C, D and E, effectively reducing the low-speed area of the channel, reducing the reverse pressure gradient, delaying the separation position, and improving the aerodynamic performance of the downstream outlet F.  Mass d_ori_70k mass flow at the original impeller design point at 70,000 rpm Mass o_70k mass flow at the optimised impeller design point at 70,000 rpm Mass d_ori_100k mass flow at the original impeller design point at 100,000 rpm Mass o_100k mass flow at the optimised impeller design point at 100,000 rpm T pr d_ori_70k total pressure ratios at the original impeller design point at 70,000 rpm T pr o_70k total pressure ratios at the optimised impeller design point at 70,000 rpm T pr d_ori_100k total pressure ratios at the original impeller design point at 100,000 rpm T pr o_100k total pressure ratios at the optimised impeller design point at 100,000 rpm