Aerodynamic Optimization Design of a 150 kW High Performance Supercritical Carbon Dioxide Centrifugal Compressor without a High Speed Requirement

Supercritical carbon dioxide (S-CO2) Brayton cycle technology has the advantages of excellent energy density and heat transfer. The compressor is the most critical and complex component of the cycle. Especially, in order to make the system more reliable and economical, the design method of a high efficiency compressor without a high speed requirement is particularly important. In this paper, thermodynamic design software of a S-CO2 centrifugal compressor is developed. It is used to design the 150 kW grade S-CO2 compressor at the speed of 40,000 rpm. The performance of the initial design is carried out by a 3-D aerodynamic analysis. The aerodynamic optimization includes three aspects: numerical calculation, design software and the flow part geometry parameters. The aerodynamic performance and the off-design performance of the optimal design are obtained. The results show that the total static efficiency of the compressor is 79.54%. The total pressure ratio is up to 1.9. The performance is excellent, and it can operate normally within the mass flow rate range of 5.97 kg/s to 11.05 kg/s. This research provides an intelligent and efficient design method for S-CO2 centrifugal compressors with a low flow rate and low speed, but high pressure ratio.


Introduction
Energy, environment and development are three major themes we are facing today. The extensive use of fossil energy has posed a great threat to the living space of mankind. CO 2 is a low-cost fluid with a low critical point (31.1 • C and 7.38 MPa). It is non-toxic and non-combustible and has great thermal stability, physical properties and safety [1]. The working temperature of the S-CO 2 Brayton cycle is above the critical temperature of carbon dioxide. Additionally, S-CO 2 has good transitivity and fast mobility, and its density is close to that of liquid. So, it can make the pressure of the fluid high [2,3]. Meanwhile, the unique properties of microsupercritical CO 2 can reduce compressor power consumption and improve cycle efficiency. Therefore, S-CO 2 Brayton cycle technology has been adopted in thermal power, nuclear power, solar power generation and so on [4]. As a key component, compressor works near the critical point of carbon dioxide. The local condensation is easy to occur, so the design of compressor is very difficult. It can be seen from Table 1 that, due to the high density characteristics of S-CO 2 , the centrifugal compressors designed by various research institutions have a smaller diameter, compact structure and higher speed. At present, the design of S-CO 2 compressor is not perfect. The design of the diffuser, bearing, seal and other auxiliary parts is difficult. Therefore, the S-CO 2 compressor is usually not able to reach the design speed in the actual operation. The experimental speed is far below the design value. For instance, in the experiment of TIT, the speed of turbine compressor sets is only 55,000 r/min [7]. In the experiment of KIER, the speed of turbine compressor sets is 30,000 r/min [8]. Moreover, the internal flow field of the compressor is far from the design point. It will lead to greater efficiency and power changes. For example, when the system of TIT is running, the centrifugal compressor efficiency is only 48%, and the circulating power generation capacity is only 0.11 kW [9]. Therefore, a series of problems such as aerodynamic design, actual operation rule, bearing and gas seal structure design, diffuser matching problem and material property changes of the S-CO 2 centrifugal compressor need further study. In particular, it is necessary to study the design method of a high efficiency compressor at low design speed.
In recent years, the mechanism of the S-CO 2 centrifugal compressor has been increasingly studied by scholars of various countries. The design of the S-CO 2 centrifugal compressor with a low speed coefficient is studied by Lettieri et al. of MIT (Massachusetts Institute of Technology) [10]. It is found that the use of the vane diffuser can improve the efficiency of the compressor in this case. Monje et al. [11] studied the design method of the compressor in the S-CO 2 Brayton cycle system, and the key parameter selection of the one-dimensional design program and multidimensional design method are introduced. Budinis et al. [12] have carried on the design and the analysis to the control system of the S-CO 2 compressor, and studied the variable operating curves and the surge control methods in detail. Shao et al. [13] of Dalian University of Technology introduced the concept of 'condensate margin' in the design process of the S-CO 2 centrifugal compressor to evaluate the working fluid state of the impeller inlet.
In summary, the compressor works in the microsupercritical point. Thus, the change of S-CO 2 physical property depends on the inlet conditions to some extent. The change of density with the pressure is also different from the ideal gas. Therefore, the existing scientific principles and the formulas related to avoiding surge may no longer be applicable to this system. At present, the compressor aerodynamic design lacks the empirical parameter and the range of estimated parameter corresponding to S-CO 2 . This makes the design more difficult. In addition, the design cycle is greatly increased, and designers often need a lot of experience. Especially, high speed can improve the efficiency of compressor. However, it will greatly improve the design difficulty and cost of the high-speed motor. It will also affect the reliability of bearing and system. Based on the research status and difficulties mentioned above, the thermal design of a 150 kW high performance S-CO 2 centrifugal compressor without a high speed requirement is carried out using the design software. The aerodynamic analysis, aerodynamic optimization and off-design performance analysis are also carried out.

Design Method
Based on our previous design experience and methods [14,15], the thermodynamic design software of S-CO 2 centrifugal compressor (S-CO 2 CPTD) was developed. The software was adopted in the compressor thermodynamic design of this research. The software has five main functions: data input, core calculation, data output, drawing and secondary exploration. The software is designed by Visual Basic 6.0 and Intel Fortran 2017.
Based on 1-D flow theory, the following main equations were used in the core calculation.
(1) Euler Equation According to the energy transformation and conservation law, the mechanical energy transferred to rotor blades is converted into fluid energy. Thus, the energy gained by one kilogram of fluid is: where h th (J/kg) is the energy head gained by one kilogram of fluid in the impeller blade channel, g (m/s 2 ) is the acceleration of gravity, c (m/s) is the absolute speed and u (m/s) is the circumferential speed. The subscript 1 represents the impeller inlet, the subscript 2 represents the impeller outlet and the subscript u represents the circumferential component.
(2) Energy Equation The total power of per kilogram of fluid consumed in centrifugal compressor stages is considered to consist of three parts: the work of the impeller on the fluid, the loss of internal leakage and the loss of impeller resistance. Based on the energy equation, the work of the impeller on fluid is converted to the fluid energy.
where h tot (J/kg) is the total energy head, h df (J/kg) is the loss of impeller resistance and h t (J/kg) is the loss of internal leakage. In Formula (3), the theoretical energy head is converted to the working substance in the form of mechanical energy. The loss of internal leakage and the loss of impeller resistance are converted to fluid in the form of heat. Therefore, the total energy equation can be written as: where c p (J/(kg·K)) is the specific heat capacity at constant pressure, A is the thermal equivalent of work (J/cal) and T (K) is the temperature. (

3) Bernoulli Equation
For the whole compressor stage, the work of the impeller on the fluid is converted to the following three parts: (1) Improving the static pressure energy of fluid. (2) Improving the kinetic energy of fluid, but in general, the kinetic energy increases little, and is often negligible. (3) Overcoming the flow loss of working fluids in stages.
By adopting the Bernoulli equation in the impeller, the following equation can be obtained: where h f2 (J/kg) represents the flow loss of fluid flowing in the diffuser. The subscript 2 represents the diffuser inlet and 3 represents the diffuser outlet.
(4) Key geometric design Equations The thermal design is mainly to calculate the flow in the centrifugal compressor impeller to determine the geometric parameters. According to the design conditions and estimated parameters, the impeller outer diameter can be obtained: where . m (kg/s) is the mass flow rate, ρ (kg/m 3 ) is density, τ 1 is the blocking factor of the impeller outlet, τ 2 = l 2 /D is the relative width of the impeller outlet and l 2 (mm) is the outlet blade height of impeller. The subscript r represents the radial component.
The inlet blade height of the impeller can be obtained from the following equation: where τ 4 is the leakage loss coefficient of the impeller cover, q (m 3 /s) is the volume flow rate, τ 3 is the blocking factor of the impeller inlet and D' (mm) is the outer diameter of impeller inlet section. In this study, the radial linear design method [16] was adopted for the impeller blade, and the cylindrical parabola geometric design method was adopted for the blade profile, as shown in Figure 1. In Figure 1, m is the thickness of the blade outlet, n is the thickness of the blade inlet and z' is the axial length of blade at different radius positions. The detailed methods and specific equations are given by Reference 16. The diffuser is in the form of the airfoil blade, and the optimal design is realized by optimizing the profile and the geometry angle of the inlet and outlet.
Appl. Sci. 2020, 10, x FOR PEER REVIEW  4 of 21 where p (Pa) is the pressure, γ (N/m 3 ) is the specific weight and hf1 (J/kg) is the flow loss of fluid flowing in the impeller.
When adopting the Bernoulli equation in the diffuser, hth is zero: where hf2 (J/kg) represents the flow loss of fluid flowing in the diffuser. The subscript 2 represents the diffuser inlet and 3 represents the diffuser outlet.
(4) Key geometric design Equations The thermal design is mainly to calculate the flow in the centrifugal compressor impeller to determine the geometric parameters. According to the design conditions and estimated parameters, the impeller outer diameter can be obtained: where m (kg/s) is the mass flow rate, ρ (kg/m 3 ) is density, τ1 is the blocking factor of the impeller outlet, τ2=l2/D is the relative width of the impeller outlet and l2 is the outlet blade height of impeller. The subscript r represents the radial component.
The inlet blade height of the impeller can be obtained from the following equation: where τ4 is the leakage loss coefficient of the impeller cover, q (m 3 /s) is the volume flow rate, τ3 is the blocking factor of the impeller inlet and D' (mm) is the outer diameter of impeller inlet section.
In this study, the radial linear design method [16] was adopted for the impeller blade, and the cylindrical parabola geometric design method was adopted for the blade profile, as shown in Figure  1. In Figure 1, m is the thickness of the blade outlet, n is the thickness of the blade inlet and z' is the axial length of blade at different radius positions. The detailed methods and specific equations are given by Reference 16. The diffuser is in the form of the airfoil blade, and the optimal design is realized by optimizing the profile and the geometry angle of the inlet and outlet.

The Intial Design
The key parameters of the initial design of compressor are shown in Table 2. The design speed was 40,000 rpm. There were 12 main blades and 12 splitter blades. The rotor's meridian face profile and 3-D model are shown in Figure 2. On the basis of the initial design, the Case A (9 main blades + 9 splitter blades) and the Case B (12 main blades) were designed to analyze the influence of blade form and number on compressor performance. The geometric parameters and blade profiles of the initial design, Case A and Case B were the same.

The Intial Design
The key parameters of the initial design of compressor are shown in Table 2. The design speed was 40,000 rpm. There were 12 main blades and 12 splitter blades. The rotor's meridian face profile and 3-D model are shown in Figure 2. On the basis of the initial design, the Case A (9 main blades + 9 splitter blades) and the Case B (12 main blades) were designed to analyze the influence of blade form and number on compressor performance. The geometric parameters and blade profiles of the initial design, Case A and Case B were the same.

Optimization Method
At present, the empirical coefficients and estimated parameters matching the thermal design of the S-CO2 centrifugal compressor have not been disclosed. The size of the flow passage obtained by the 1-D flow method does not necessarily conform to the actual flow requirements. At the same time, the design and optimization period of the traditional design method is long. Therefore, based on the 1-D flow theory, a fast thermal design method of S-CO2 compressor was established in this study. Combined with the high-precision three-dimensional aerodynamic analysis method, a designoptimization method based on Gauss process regression was proposed, as shown in Figure 3.

Optimization Method
At present, the empirical coefficients and estimated parameters matching the thermal design of the S-CO 2 centrifugal compressor have not been disclosed. The size of the flow passage obtained by the 1-D flow method does not necessarily conform to the actual flow requirements. At the same time, the design and optimization period of the traditional design method is long. Therefore, based on the 1-D flow theory, a fast thermal design method of S-CO 2 compressor was established in this study. Combined with the high-precision three-dimensional aerodynamic analysis method, a design-optimization method based on Gauss process regression was proposed, as shown in Figure 3.
According to the thermal design results and aerodynamic design results, and based on the Gauss process regression assumption, the true isentropic efficiency of an unknown design condition y can be estimated:f where ∆ f i (y) is the deviation between the thermal design f a (y) and CFD accurate calculation results. The calculation deviation satisfies Gauss distribution at each point in the whole design space: Appl. Sci. 2020, 10, 2093 6 of 20 The mean µ(y) and variance σ(y) of Gaussian distribution are: where D is the training set, R yD and R Dy are the covariance vector of the calculated condition and the new sampling point, R DD is the covariance matrix and ∆ f D is the deviation between the thermal design results and CFD accurate calculation results of all samples in the training set. Therefore, the real isentropic efficiency (aerodynamic analysis result) of an unknown design condition y is estimated as follows: where f r (y) is the lower limit estimation, f r (y) is the upper limit estimation and ξ is the confidence constant.  According to the thermal design results and aerodynamic design results, and based on the Gauss process regression assumption, the true isentropic efficiency of an unknown design condition y can be estimated: where () i f  y is the deviation between the thermal design () a f y and CFD accurate calculation results.
The calculation deviation satisfies Gauss distribution at each point in the whole design space: The mean ()  y and variance ()  y of Gaussian distribution are: where D is the training set, Therefore, the real isentropic efficiency (aerodynamic analysis result) of an unknown design condition y is estimated as follows: This method has the characteristics of having fast speed and being self-adaptive. The detailed optimization theory and method can refer to Reference [17]. Key geometric parameters were used as optimization variables, including: the relative width of the impeller outlet, inlet and outlet geometric angle, leading edge sweep angle, etc. Generally, the better design of a given blade profile can be obtained by only 30 steps of iteration (i.e., the number of the aerodynamic analysis). The optimization time is only 1-2 days. According to the better design results, three aspects are optimized: numerical calculation, compressor design software and flow passage geometry parameters. The optimization process is shown in Figure 4. The numerical calculation is improved by changing the numerical calculation method and optimizing the physical property database. By changing the control parameters of numerical calculation, the stability and speed of convergence can be improved. The interpolation density of physical database is optimized to improve the calculation speed on the premise of ensuring the calculation accuracy. The compressor design software is modified by adjusting the experience coefficient, the loss model and the design method through the thermal design and Gauss process regression method. Meanwhile, the flow passage geometric parameters are optimized, such as the blade profile, geometric angle and wrap angle, the axial length of the impeller, etc. In addition, the matching of diffuser and impeller should be considered. Iterating the optimal process until the optimal design is achieved. calculation method and optimizing the physical property database. By changing the control parameters of numerical calculation, the stability and speed of convergence can be improved. The interpolation density of physical database is optimized to improve the calculation speed on the premise of ensuring the calculation accuracy. The compressor design software is modified by adjusting the experience coefficient, the loss model and the design method through the thermal design and Gauss process regression method. Meanwhile, the flow passage geometric parameters are optimized, such as the blade profile, geometric angle and wrap angle, the axial length of the impeller, etc. In addition, the matching of diffuser and impeller should be considered. Iterating the optimal process until the optimal design is achieved.

The Optimal Design
The key parameters of the optimal design of the compressor are shown in Table 3. The design speed is 40,000 rpm. There were 12 main blades and 17 diffuser blades. Figure 5 and Figure 6 respectively show the profile and 3-D model of the impeller and diffuser in the optimal design. Table 3. Key parameters of the optimal design of the compressor.

Key Parameters Value Unit
Thermodynamic parameters

The Optimal Design
The key parameters of the optimal design of the compressor are shown in Table 3. The design speed is 40,000 rpm. There were 12 main blades and 17 diffuser blades. Figures 5 and 6 respectively show the profile and 3-D model of the impeller and diffuser in the optimal design. Table 3. Key parameters of the optimal design of the compressor.

Key Parameters Value Unit
Thermodynamic parameters

Boundary Conditons
In this research, NUMECA-FineTurbo was used to solve the three-dimensional N-S equations. The single impeller-diffuser flow passage was used as the calculation model. S-CO2 was used as the working substance. According to the literature [18] and related data (NIST database), the thermal physical data of microsupercritical and microsubcritical CO2 were integrated. The thermal physical property continuous function was constructed by the Kriging surrogate model. The points near the critical point were locally encrypted because of the sharp change of physical properties near the critical point. The dynamic database of thermal physical properties of S-CO2 was established for invocation, which can ensure the accuracy and convergence of calculation. The pressure range of working medium was set to 1-20 MPa, and the temperature range was set to 273.15-600 K. The fluid model included the parameters of all the phases of CO2 in the above pressure and temperature range. The Shear Stress Transport (extended wall function) turbulence model was selected. Some other scholars [19][20][21] have also used the turbulence model in the analysis of centrifugal compressor, and obtained reasonable results. For the impeller inlet, the total temperature was 309.15 K and the total pressure was 7.8 MPa. k (5 m 2 /s 2 ) and epsilon (30,000 m 2 /s 3 ) were selected as turbulent quantities. For the diffuser outlet, the flow rate was set as 6.66 kg/s, and the influence of backflow was taken into consideration. The impeller fluid domain had a rotational speed of 40,000 rpm around the Z axis. The cross section between the impeller outlet and diffuser inlet was set as a coupling interface with

Boundary Conditons
In this research, NUMECA-FineTurbo was used to solve the three-dimensional N-S equations. The single impeller-diffuser flow passage was used as the calculation model. S-CO2 was used as the working substance. According to the literature [18] and related data (NIST database), the thermal physical data of microsupercritical and microsubcritical CO2 were integrated. The thermal physical property continuous function was constructed by the Kriging surrogate model. The points near the critical point were locally encrypted because of the sharp change of physical properties near the critical point. The dynamic database of thermal physical properties of S-CO2 was established for invocation, which can ensure the accuracy and convergence of calculation. The pressure range of working medium was set to 1-20 MPa, and the temperature range was set to 273.15-600 K. The fluid model included the parameters of all the phases of CO2 in the above pressure and temperature range. The Shear Stress Transport (extended wall function) turbulence model was selected. Some other scholars [19][20][21] have also used the turbulence model in the analysis of centrifugal compressor, and obtained reasonable results. For the impeller inlet, the total temperature was 309.15 K and the total pressure was 7.8 MPa. k (5 m 2 /s 2 ) and epsilon (30,000 m 2 /s 3 ) were selected as turbulent quantities. For the diffuser outlet, the flow rate was set as 6.66 kg/s, and the influence of backflow was taken into consideration. The impeller fluid domain had a rotational speed of 40,000 rpm around the Z axis. The cross section between the impeller outlet and diffuser inlet was set as a coupling interface with

Boundary Conditons
In this research, NUMECA-FineTurbo was used to solve the three-dimensional N-S equations. The single impeller-diffuser flow passage was used as the calculation model. S-CO 2 was used as the working substance. According to the literature [18] and related data (NIST database), the thermal physical data of microsupercritical and microsubcritical CO 2 were integrated. The thermal physical property continuous function was constructed by the Kriging surrogate model. The points near the critical point were locally encrypted because of the sharp change of physical properties near the critical point. The dynamic database of thermal physical properties of S-CO 2 was established for invocation, which can ensure the accuracy and convergence of calculation. The pressure range of working medium was set to 1-20 MPa, and the temperature range was set to 273.15-600 K. The fluid model included the parameters of all the phases of CO 2 in the above pressure and temperature range. The Shear Stress Transport (extended wall function) turbulence model was selected. Some other scholars [19][20][21] have also used the turbulence model in the analysis of centrifugal compressor, and obtained reasonable results. For the impeller inlet, the total temperature was 309.15 K and the total pressure was 7.8 MPa. k (5 m 2 /s 2 ) and epsilon (30,000 m 2 /s 3 ) were selected as turbulent quantities. For the diffuser outlet, the flow rate was set as 6.66 kg/s, and the influence of backflow was taken into consideration. The impeller fluid domain had a rotational speed of 40,000 rpm around the Z axis. The cross section between the impeller outlet and diffuser inlet was set as a coupling interface with conservative coupling by the pitchwise row mixed model. The diffuser wall was set as absolutely static, and the impeller wall was set as relatively static. The shroud and hub were set as adiabatic, and the conditions of non-slip flow were satisfied. The residual less than 1 × 10 -6 and the iterative deviation less than 0.1% of the outlet temperature and impeller blade torque were considered as the convergence conditions. Figure 7 shows the grid schematic diagram of the impeller-diffuser fluid domain. In this paper, NUMECA-AutoGrid5 was used to mesh. The H type mesh was used in the inlet and outlet extension sections. For higher mesh quality, the O type mesh was used for blade meshing. The grids of the tip clearance and near the wall were refined to obtain accurate flow parameters. The value of Y + near the wall was about 1, which met the calculation requirements of the SST turbulence model.

Mesh and Grid Independence
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 21 conservative coupling by the pitchwise row mixed model. The diffuser wall was set as absolutely static, and the impeller wall was set as relatively static. The shroud and hub were set as adiabatic, and the conditions of non-slip flow were satisfied. The residual less than 1 × 10 -6 and the iterative deviation less than 0.1% of the outlet temperature and impeller blade torque were considered as the convergence conditions. Figure 7 shows the grid schematic diagram of the impeller-diffuser fluid domain. In this paper, NUMECA-AutoGrid5 was used to mesh. The H type mesh was used in the inlet and outlet extension sections. For higher mesh quality, the O type mesh was used for blade meshing. The grids of the tip clearance and near the wall were refined to obtain accurate flow parameters. The value of Y + near the wall was about 1, which met the calculation requirements of the SST turbulence model. In order to effectively utilize computing resources, the grid independence verification was carried out. The optimal design of the compressor was taken as an example. The grid independence is shown in Figure 8. With the increase of the grid number, the accuracy of the calculation model also increased. More details of flow loss could be captured in the numerical calculation, so the efficiency was gradually reduced. When the grid number increased from 800,000 to 1,600,000, the relative change of total static efficiency was less than 0.2%. So, when the grid number was more than 800,000, the calculation model could meet the demand of the calculation accuracy.  In order to effectively utilize computing resources, the grid independence verification was carried out. The optimal design of the compressor was taken as an example. The grid independence is shown in Figure 8. With the increase of the grid number, the accuracy of the calculation model also increased. More details of flow loss could be captured in the numerical calculation, so the efficiency was gradually reduced. When the grid number increased from 800,000 to 1,600,000, the relative change of total static efficiency was less than 0.2%. So, when the grid number was more than 800,000, the calculation model could meet the demand of the calculation accuracy.

Mesh and Grid Independence
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 21 conservative coupling by the pitchwise row mixed model. The diffuser wall was set as absolutely static, and the impeller wall was set as relatively static. The shroud and hub were set as adiabatic, and the conditions of non-slip flow were satisfied. The residual less than 1 × 10 -6 and the iterative deviation less than 0.1% of the outlet temperature and impeller blade torque were considered as the convergence conditions. Figure 7 shows the grid schematic diagram of the impeller-diffuser fluid domain. In this paper, NUMECA-AutoGrid5 was used to mesh. The H type mesh was used in the inlet and outlet extension sections. For higher mesh quality, the O type mesh was used for blade meshing. The grids of the tip clearance and near the wall were refined to obtain accurate flow parameters. The value of Y + near the wall was about 1, which met the calculation requirements of the SST turbulence model. In order to effectively utilize computing resources, the grid independence verification was carried out. The optimal design of the compressor was taken as an example. The grid independence is shown in Figure 8. With the increase of the grid number, the accuracy of the calculation model also increased. More details of flow loss could be captured in the numerical calculation, so the efficiency was gradually reduced. When the grid number increased from 800,000 to 1,600,000, the relative change of total static efficiency was less than 0.2%. So, when the grid number was more than 800,000, the calculation model could meet the demand of the calculation accuracy.  The total static efficiency was determined by the output power, mass flow rate, isentropic enthalpy rise and inlet velocity. The isentropic enthalpy drop and the total static efficiency are shown in the Formulas (12) and (13), respectively.

Mesh and Grid Independence
where κ is the adiabatic exponent and P (W) is the power. The subscript inlet represents the compressor inlet and the subscript outlet represents the compressor outlet.

Numerical Validation
In order to verify the accuracy of the numerical method, the single-stage compressor model of SNL was adopted. The key geometric parameters of the model are summarized in Table 4. The case of the total pressure of 7.69 MPa, total temperature of 305.95 K and rotation speed of 50,000 rpm was analyzed. The verification results were compared with the numerical results from Rinaldi et al. [22] and the experimental results from Wright et al. [5], as shown in Figure 9. The calculated efficiency curve was quite close to that from Rinaldi et al., and shows the same trend. However, there was a certain error between the setting of boundary conditions in the numerical calculation and the actual conditions in the experiment. Therefore, the difference between the numerical calculation and experiment was also acceptable. In general, the numerical method in this study was accurate.   Figure 10 shows the velocity vector and partial enlargement of the 50% blade height section of  Figure 10 shows the velocity vector and partial enlargement of the 50% blade height section of the compressor initial design.  Figure 10 shows the velocity vector and partial enlargement of the 50% blade height section of the compressor initial design. It can be seen from the diagram that there was no obvious flow separation phenomenon inside the impeller. However, the flow inside the impeller was not uniform. When the fluid flowed along the main blade to 50% of the chord length (the leading edge (LE) of the splitter blade), the flow velocity of the suction side (SS) increased rapidly, while the fluid on the pressure side (PS) still flowed at a lower velocity. This led to a larger pressure gradient between the PS and SS. Thus, the working fluid on the PS would leak through the tip clearance to the SS. It would lead to a significant reduction in compressor performance inevitably. Figure 11 shows the pressure distribution of the impeller blades surface in detail. Generally speaking, the PS pressure was greater than SS pressure. The pressure distribution on the splitter blade surface was basically consistent with that on the corresponding area of the main blade. Under the action of centrifugal force, the pressure of CO2 increased gradually along the blade profile. There was a larger pressure gradient at the trailing edge (TE). There was no reverse pressure region in the It can be seen from the diagram that there was no obvious flow separation phenomenon inside the impeller. However, the flow inside the impeller was not uniform. When the fluid flowed along the main blade to 50% of the chord length (the leading edge (LE) of the splitter blade), the flow velocity of the suction side (SS) increased rapidly, while the fluid on the pressure side (PS) still flowed at a lower velocity. This led to a larger pressure gradient between the PS and SS. Thus, the working fluid on the PS would leak through the tip clearance to the SS. It would lead to a significant reduction in compressor performance inevitably. Figure 11 shows the pressure distribution of the impeller blades surface in detail. Generally speaking, the PS pressure was greater than SS pressure. The pressure distribution on the splitter blade surface was basically consistent with that on the corresponding area of the main blade. Under the action of centrifugal force, the pressure of CO 2 increased gradually along the blade profile. There was a larger pressure gradient at the trailing edge (TE). There was no reverse pressure region in the pressure distribution. That means there was no flow separation phenomenon on the blade surface. There was a small low-pressure region at the LE of the impeller blade. It was caused by the local acceleration of the flow impact.

The Initial Design
When the pressure in the low-pressure region of the LE is lower than the critical point, the transcritical phenomenon will occur. To some extent, it will affect the performance of the compressor. In order to display the cross critical region of the blade surface more clearly, the maximum pressure of the contour was set to 7.38 MPa (CO 2 critical pressure), as shown in Figure 12. The phenomenon was mainly caused by two factors. Firstly, the working fluid entered the compressor along the axial direction, and there was a flow acceleration phenomenon at the rotor LE. The curvature change of the molded lines at the rotor LE was the largest. Thus, the gradient of the velocity increase was correspondingly the largest. At the corner position of the LE on the top of the blade, the accelerating fluid on both sides of the blade would cause an 'ejection' effect. The gradient of the velocity increase was the largest at this position, so the fluid velocity increased sharply. It would produce an obviously low pressure and low temperature region. CO 2 entered the two-phase region and was far below the critical state. So, it would result with the possibility of 'condensation'. Secondly, the phenomenon of an 'off vortex at the LE on the top of the blade' near this position would cause an obviously low-pressure area, and the working fluid would enter the two-phase region. This is consistent with the conclusions of other scholars [23,24].
Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 21 pressure distribution. That means there was no flow separation phenomenon on the blade surface.
There was a small low-pressure region at the LE of the impeller blade. It was caused by the local acceleration of the flow impact. Figure 11. The pressure distribution on the surface of impeller blades.
When the pressure in the low-pressure region of the LE is lower than the critical point, the transcritical phenomenon will occur. To some extent, it will affect the performance of the compressor. In order to display the cross critical region of the blade surface more clearly, the maximum pressure of the contour was set to 7.38 MPa (CO2 critical pressure), as shown in Figure 12. The phenomenon was mainly caused by two factors. Firstly, the working fluid entered the compressor along the axial direction, and there was a flow acceleration phenomenon at the rotor LE. The curvature change of the molded lines at the rotor LE was the largest. Thus, the gradient of the velocity increase was correspondingly the largest. At the corner position of the LE on the top of the blade, the accelerating fluid on both sides of the blade would cause an 'ejection' effect. The gradient of the velocity increase was the largest at this position, so the fluid velocity increased sharply. It would produce an obviously low pressure and low temperature region. CO2 entered the two-phase region and was far below the critical state. So, it would result with the possibility of 'condensation'. Secondly, the phenomenon of an 'off vortex at the LE on the top of the blade' near this position would cause an obviously lowpressure area, and the working fluid would enter the two-phase region. This is consistent with the conclusions of other scholars [23,24].   When the pressure in the low-pressure region of the LE is lower than the critical point, the transcritical phenomenon will occur. To some extent, it will affect the performance of the compressor. In order to display the cross critical region of the blade surface more clearly, the maximum pressure of the contour was set to 7.38 MPa (CO2 critical pressure), as shown in Figure 12. The phenomenon was mainly caused by two factors. Firstly, the working fluid entered the compressor along the axial direction, and there was a flow acceleration phenomenon at the rotor LE. The curvature change of the molded lines at the rotor LE was the largest. Thus, the gradient of the velocity increase was correspondingly the largest. At the corner position of the LE on the top of the blade, the accelerating fluid on both sides of the blade would cause an 'ejection' effect. The gradient of the velocity increase was the largest at this position, so the fluid velocity increased sharply. It would produce an obviously low pressure and low temperature region. CO2 entered the two-phase region and was far below the critical state. So, it would result with the possibility of 'condensation'. Secondly, the phenomenon of an 'off vortex at the LE on the top of the blade' near this position would cause an obviously lowpressure area, and the working fluid would enter the two-phase region. This is consistent with the conclusions of other scholars [23,24].  On the basis of the initial design of the compressor, this paper also contrasted and analyzed Case A and Case B. The calculated performance parameters are shown in Table 5. As can be seen from the table, compared with the initial design, although the total pressure ratio of Case B decreased from 1.67 to 1.65, the total static efficiency increased from 50.3% to 64.56% greatly. The pressure ratio of Case A was 1.67 and the total static efficiency was 63.36%. It was slightly lower than the efficiency of Case B. If the blade number is adjusted, the performance of the compressor may be equivalent to that of Case B, and may even be slightly higher. Compared with the numerical results, it was found that the addition of the splitter blade had little influence on the streamline, pressure and temperature distribution in the whole fluid domain. This is because the S-CO 2 compressor had a higher speed and a more compact structure than the centrifugal compressor with a conventional working medium. Therefore, it is only necessary to set the appropriate number of main blades to restrict the flow at the impeller outlet and make the outlet flow angle meet the requirements of the design value. Additionally, the compressor design with the splitter blade structure had a smaller flow capacity and larger torque. This would lead to a significant decrease in the compressor efficiency. Besides, the addition of splitter blade structure would greatly increase manufacturing difficulty and manufacturing cost. Therefore, considering the factors such as performance, machining and economy, the splitter blade structure was not considered in the following research. Overall, the compressor performance of initial design was quite different from the design value, thus the three-dimensional aerodynamic optimization is needed.  Figure 13 shows the velocity vector distribution of the 50% blade height section for the compressor optimal design, and the partial enlargement of the LE and TE of rotor and diffuser blades. It can be seen from the diagram that there was no serious flow separation phenomenon and backflow phenomenon in the impeller and diffuser. The flow field was uniform. The local acceleration phenomenon was inevitable only in the LE of the impeller blade and diffuser blade. In addition, there was a small vortex produced by the impeller blade blunt TE. distribution in the whole fluid domain. This is because the S-CO2 compressor had a higher speed and a more compact structure than the centrifugal compressor with a conventional working medium. Therefore, it is only necessary to set the appropriate number of main blades to restrict the flow at the impeller outlet and make the outlet flow angle meet the requirements of the design value. Additionally, the compressor design with the splitter blade structure had a smaller flow capacity and larger torque. This would lead to a significant decrease in the compressor efficiency. Besides, the addition of splitter blade structure would greatly increase manufacturing difficulty and manufacturing cost. Therefore, considering the factors such as performance, machining and economy, the splitter blade structure was not considered in the following research. Overall, the compressor performance of initial design was quite different from the design value, thus the three-dimensional aerodynamic optimization is needed.  Figure 13 shows the velocity vector distribution of the 50% blade height section for the compressor optimal design, and the partial enlargement of the LE and TE of rotor and diffuser blades. It can be seen from the diagram that there was no serious flow separation phenomenon and backflow phenomenon in the impeller and diffuser. The flow field was uniform. The local acceleration phenomenon was inevitable only in the LE of the impeller blade and diffuser blade. In addition, there was a small vortex produced by the impeller blade blunt TE.  Figure 14 shows the secondary flow velocity vector with different axial chord lengths on the S3 cross section. The vortex structure and development process of the compressor impeller can be analyzed by secondary flow velocity. For the horseshoe vortex, its intensity and influence range were obviously weaker than that of the radial inflow turbine. Only horseshoe vortex bifurcation on the PS (HVp) was found near the LE of the rotor blade, while horseshoe vortex bifurcation on the SS (HVs) was not found. Under the pressure gradient, the boundary layer was rolled up and the horseshoe vortex bifurcation was formed on the PS. As the fluid flows downstream, the influence range increased continuously. It finally merged with the passage vortex under the dissipation of the adverse pressure gradient. Additionally, the difference between the passage vortex and other vortex structures was not obvious. Figure 14a shows that there was an up passage vortex (PVu) and down passage vortex (PVd) with the same rotational direction at the tip and root of the blade on the PS, respectively. However, the intensity was small and the influence range was narrow. The up passage vortexes were generated mainly by the interaction of three parts. The first part was the scraping flow produced by the relative motion of the upper end wall due to the high-speed rotation of the impeller. The second part was the leakage flow from the PS to the SS in the tip clearance. Another part was the effect of transverse pressure gradient. As the fluid flowed along the radial direction, the scraping effect caused by the high-speed flow of the impeller became stronger, and the vortices on the PS expanded gradually. The influence range of vortex in the down passage gradually decreased. Its position was gradually squeezed into the SS and eventually dissipated and disappeared, as shown in Figure 14b. Near the TE of the rotor blade, the pressure on the SS was basically the same as that on the PS. The leakage flow and the transverse pressure gradient disappear, and the vortex system was dissipated and disappeared, as shown in Figure 14c.

The Optimal Design
analyzed by secondary flow velocity. For the horseshoe vortex, its intensity and influence range were obviously weaker than that of the radial inflow turbine. Only horseshoe vortex bifurcation on the PS (HVp) was found near the LE of the rotor blade, while horseshoe vortex bifurcation on the SS(HVs) was not found. Under the pressure gradient, the boundary layer was rolled up and the horseshoe vortex bifurcation was formed on the PS. As the fluid flows downstream, the influence range increased continuously. It finally merged with the passage vortex under the dissipation of the adverse pressure gradient. Additionally, the difference between the passage vortex and other vortex structures was not obvious. Figure 14a shows that there was an up passage vortex (PVu) and down passage vortex (PVd) with the same rotational direction at the tip and root of the blade on the PS, respectively. However, the intensity was small and the influence range was narrow. The up passage vortexes were generated mainly by the interaction of three parts. The first part was the scraping flow produced by the relative motion of the upper end wall due to the high-speed rotation of the impeller. The second part was the leakage flow from the PS to the SS in the tip clearance. Another part was the effect of transverse pressure gradient. As the fluid flowed along the radial direction, the scraping effect caused by the high-speed flow of the impeller became stronger, and the vortices on the PS expanded gradually. The influence range of vortex in the down passage gradually decreased. Its position was gradually squeezed into the SS and eventually dissipated and disappeared, as shown in Figure 14b. Near the TE of the rotor blade, the pressure on the SS was basically the same as that on the PS. The leakage flow and the transverse pressure gradient disappear, and the vortex system was dissipated and disappeared, as shown in Figure 14c.  Figure 15 shows the pressure and temperature distribution of the 50% blade height section for the compressor optimal design. Generally speaking, the pressure and temperature of the PS were greater than that of the SS. Under the action of centrifugal force, the CO 2 pressure and temperature from the LE to the TE of the rotor gradually increased, and there was a larger gradient at the TE. The diffusing action of diffuser maximized the outlet pressure and outlet temperature. The outlet static pressure reached 14.25 MPa. Figure 15 shows the pressure and temperature distribution of the 50% blade height section for the compressor optimal design. Generally speaking, the pressure and temperature of the PS were greater than that of the SS. Under the action of centrifugal force, the CO2 pressure and temperature from the LE to the TE of the rotor gradually increased, and there was a larger gradient at the TE. The diffusing action of diffuser maximized the outlet pressure and outlet temperature. The outlet static pressure reached 14.25 MPa.
(a) (b) Figure 15. The pressure and temperature distribution of the 50% blade height section for the optimal design of the compressor: (a) pressure distribution and (b) temperature distribution. Figure 16 shows the surface pressure distributions of the impeller blades and diffuser blades of the compressor optimal design. The blade surface pressure curves with different blade height sections are shown in Figure 17. In general, the pressure increased gradually from the inlet to the outlet, the pressure gradient distribution was uniform. The maximum pressure was at the diffuser outlet. Due to the decrease of the CO2 leakage caused by the tip clearance, the surface pressure of rotor was basically the same along the blade height direction from the position of the 50% chord length to the TE. The cross critical region of the blade surface is shown in Figure 18. Similarly, the fluid in the blade tip clearance at the impeller blade LE was 'ejected'. Meanwhile, the fluid at the LE of the impeller blade would be affected by the phenomenon of an 'off vortex at the LE on the top of blade'. The simultaneous action of the two phenomena led to obviously low-pressure-temperature regions at the impeller blade LE, and the carbon dioxide entered the two-phase region. The fluid state at the corner position of the LE on the top of the blade was far below the critical state. This might cause 'condensation'. The flow in the transcritical region was very complex and the physical properties changed dramatically. This would cause great loss and affect the flow condition in the impeller. By optimizing the geometry of the flow passage, the influence range of the transcritical phenomenon in the optimal design was obviously weakened. It was mainly concentrated on the SS of the LE.  Figure 16 shows the surface pressure distributions of the impeller blades and diffuser blades of the compressor optimal design. The blade surface pressure curves with different blade height sections are shown in Figure 17. In general, the pressure increased gradually from the inlet to the outlet, the pressure gradient distribution was uniform. The maximum pressure was at the diffuser outlet. Due to the decrease of the CO 2 leakage caused by the tip clearance, the surface pressure of rotor was basically the same along the blade height direction from the position of the 50% chord length to the TE. The cross critical region of the blade surface is shown in Figure 18. Similarly, the fluid in the blade tip clearance at the impeller blade LE was 'ejected'. Meanwhile, the fluid at the LE of the impeller blade would be affected by the phenomenon of an 'off vortex at the LE on the top of blade'. The simultaneous action of the two phenomena led to obviously low-pressure-temperature regions at the impeller blade LE, and the carbon dioxide entered the two-phase region. The fluid state at the corner position of the LE on the top of the blade was far below the critical state. This might cause 'condensation'. The flow in the transcritical region was very complex and the physical properties changed dramatically. This would cause great loss and affect the flow condition in the impeller. By optimizing the geometry of the flow passage, the influence range of the transcritical phenomenon in the optimal design was obviously weakened. It was mainly concentrated on the SS of the LE.     In order to study the loss of flow in the compressor more intuitively, the entropy distribution of different blade height sections, the impeller outlet and diffuser outlet are given in Figure 19. It can be found that the flow loss inside the compressor was more and more serious from the hub to the shroud. It is consistent with the gradual deterioration of the flow situation mentioned above. The    In order to study the loss of flow in the compressor more intuitively, the entropy distribution of different blade height sections, the impeller outlet and diffuser outlet are given in Figure 19. It can be found that the flow loss inside the compressor was more and more serious from the hub to the shroud. It is consistent with the gradual deterioration of the flow situation mentioned above. The In order to study the loss of flow in the compressor more intuitively, the entropy distribution of different blade height sections, the impeller outlet and diffuser outlet are given in Figure 19. It can be found that the flow loss inside the compressor was more and more serious from the hub to the shroud. It is consistent with the gradual deterioration of the flow situation mentioned above. The entropy of the 10% blade height section had little change. The working medium was pressurized in the impeller, which is similar to the isentropic process. In the diffuser, due to the influence of the viscous boundary layer and the wake, an increase in entropy appeared on the blade surface and the extension of the outlet. At the 50% blade height section, the entropy of blade TE increased greatly under the influence of the wake. The flow loss in this region was increased. A large range of higher entropy increases occurred at the 90% blade height section near the blade top. This area is close to the wall, and the working fluid was strongly sheared by the impeller at high speed. Therefore, there was a large viscous dissipation loss. At the same time, this area was affected by the tip clearance layer. Various vortex structures were blended here, and the flow was very complicated. The entropy at the impeller outlet was evenly distributed in the circumferential direction and gradually increased along the blade height. This is because the area near the shroud was greatly affected by the tip clearance layer and the cutting of the wall surface. The entropy at the diffuser outlet was basically unchanged. There was only a slight increase in entropy on the suction side. This means that there were no backflow and secondary flow in this area. In summary, the optimal design had small flow loss and good aerodynamic performance. entropy increases occurred at the 90% blade height section near the blade top. This area is close to the wall, and the working fluid was strongly sheared by the impeller at high speed. Therefore, there was a large viscous dissipation loss. At the same time, this area was affected by the tip clearance layer. Various vortex structures were blended here, and the flow was very complicated. The entropy at the impeller outlet was evenly distributed in the circumferential direction and gradually increased along the blade height. This is because the area near the shroud was greatly affected by the tip clearance layer and the cutting of the wall surface. The entropy at the diffuser outlet was basically unchanged. There was only a slight increase in entropy on the suction side. This means that there were no backflow and secondary flow in this area. In summary, the optimal design had small flow loss and good aerodynamic performance. In the off-design analysis, eight kinds of outlet pressure boundary conditions were set up. The mass flow rate-total pressure ratio curve and mass flow rate-total static efficiency curve were obtained, as shown in Figure 20 and Figure 21. In the off-design analysis, eight kinds of outlet pressure boundary conditions were set up. The mass flow rate-total pressure ratio curve and mass flow rate-total static efficiency curve were obtained, as shown in Figures 20 and 21.
It can be seen from the figure that as the mass flow rate increased, the total pressure ratio and total static efficiency showed a decreasing trend. It is worth noting that the compressor performance was higher under the small flow condition deviating from the design condition. When the mass flow rate was 5.97 kg/s, the compressor had the highest pressure ratio and total static efficiency, 1.90% and 80% respectively. When the mass flow rate was in the range of 5.97-9.52 kg/s, the compressor performance curve changed gently. At this time, the compressor had a higher pressure ratio and efficiency. Under the condition of a mass flow rate less than 5.97 kg/s, the pressure ratio was too large, which led to the increase of outlet pressure and further caused the backflow and compressor surge. The sensitivity of mass flow rate to compressor performance was greatly increased while the mass flow rate was greater than 9.52 kg/s. The total pressure ratio and the total static efficiency decreased sharply with the increase of mass flow rate, and the blocking phenomenon was becoming more and more serious. As the mass flow rate continued to increase to 11.05 kg/s, the flow at a throat in the compressor channel reached a critical state. At this time, the flow rate was at its maximum. No matter how much the compressor back pressure was lowered, the flow rate would no longer increase. The compressor would enter the blocking condition.  It can be seen from the figure that as the mass flow rate increased, the total pressure ratio and total static efficiency showed a decreasing trend. It is worth noting that the compressor performance was higher under the small flow condition deviating from the design condition. When the mass flow rate was 5.97 kg/s, the compressor had the highest pressure ratio and total static efficiency, 1.90 and 80% respectively. When the mass flow rate was in the range of 5.97-9.52 kg/s, the compressor performance curve changed gently. At this time, the compressor had a higher pressure ratio and efficiency. Under the condition of a mass flow rate less than 5.97 kg/s, the pressure ratio was too large, which led to the increase of outlet pressure and further caused the backflow and compressor surge. The sensitivity of mass flow rate to compressor performance was greatly increased while the mass flow rate was greater than 9.52 kg/s. The total pressure ratio and the total static efficiency decreased sharply with the increase of mass flow rate, and the blocking phenomenon was becoming more and more serious. As the mass flow rate continued to increase to 11.05 kg/s, the flow at a throat in the compressor channel reached a critical state. At this time, the flow rate was at its maximum. No matter how much the compressor back pressure was lowered, the flow rate would no longer increase. The compressor would enter the blocking condition.  It can be seen from the figure that as the mass flow rate increased, the total pressure ratio and total static efficiency showed a decreasing trend. It is worth noting that the compressor performance was higher under the small flow condition deviating from the design condition. When the mass flow rate was 5.97 kg/s, the compressor had the highest pressure ratio and total static efficiency, 1.90 and 80% respectively. When the mass flow rate was in the range of 5.97-9.52 kg/s, the compressor performance curve changed gently. At this time, the compressor had a higher pressure ratio and efficiency. Under the condition of a mass flow rate less than 5.97 kg/s, the pressure ratio was too large, which led to the increase of outlet pressure and further caused the backflow and compressor surge. The sensitivity of mass flow rate to compressor performance was greatly increased while the mass flow rate was greater than 9.52 kg/s. The total pressure ratio and the total static efficiency decreased sharply with the increase of mass flow rate, and the blocking phenomenon was becoming more and more serious. As the mass flow rate continued to increase to 11.05 kg/s, the flow at a throat in the compressor channel reached a critical state. At this time, the flow rate was at its maximum. No matter how much the compressor back pressure was lowered, the flow rate would no longer increase. The compressor would enter the blocking condition. The performance of the compressor optimal design in this paper was compared with the most advanced SNL centrifugal compressor in the current public literature, as shown in Table 6. Compared with the compressor of SNL, the compressor speed of the optimal design was 40,000 rpm, and the speed was almost reduced to 50%. It is worth noting that in this paper, the design speed was greatly reduced. This means that it had more testability, lower motor cost, simpler system composition and other advantages. Meanwhile, the total pressure ratio was 1.9, bigger than 1.8 of SNL. Besides, in terms of total static efficiency, it was greatly exceeded, which was 13.16% higher than the compressor of SNL. In general, the compressor designed in this study had a low speed requirement and strong comprehensive performance.

Conclusions
In this paper, a thermodynamic design software of the S-CO 2 centrifugal compressor (S-CO 2 CPTD) was developed based on the 1-D thermal design method. The 150 kW S-CO 2 centrifugal compressor at the speed of 40,000 rpm was designed by using the developed software. The 3-D aerodynamic analysis and aerodynamic optimization was carried out in the compressor initial design. The concrete conclusions are as follows: