Investigation on Unsteady Flow Characteristics of a SCO2 Centrifugal Compressor

Supercritical carbon dioxide (SCO2) is a vital working fluid in the application of power units and its high density helps to achieve a compact mechanical structure. Centrifugal compressors are of vital use in various kinds of equipment. In this paper, a SCO2 centrifugal compressor of large input power and mass flow rate is designed and numerically investigated. A thorough numerical analysis of the unsteady flow field in the centrifugal compressor is performed in ANSYS-CFX. The computation adopts hexahedral mesh, finite volume method, and the RNG k-ε two-equation turbulence model. Streamlines, temperature, pressure, and Mach number distributions at different time steps in one revolution period are covered to present the unsteady effect of turbomachinery. Meanwhile, the force on a single rotor blade is monitored to investigate the frequency components of exciting force, thus providing the foundation for vibration analysis. Moreover, the torque, output power, pressure ratio, and isentropic efficiency in the steady and the unsteady time-averaged condition are calculated and compared with the design condition to measure the validity of the design. In summary, the unsteady computation result reveals that the unsteady flow characteristics are prominent in the designed compressor and the design of impeller and diffuser meet the requirement.


Introduction
Supercritical carbon dioxide (SCO 2 ) is a kind of supercritical fluid which has a critical point at around room temperature.Concretely, its critical temperature and pressure are respectively 304.2 K and 7.4 MPa [1].As a working fluid in a centrifugal compressor, SCO 2 possesses lots of merits.For example, its critical condition is easy to reach and it is an environmentally friendly working fluid, which makes SCO 2 safe and cheap for industrial use [2].The most useful advantages of SCO 2 are its high density, low viscosity, and large pressure increase with small enthalpy increase, which result in high efficiency, compact mechanical structure, and a larger pressure ratio with smaller input power in the centrifugal compressor [3][4][5].The research on SCO 2 has gradually been increasing lately.In recent years, the investigations on SCO 2 Brayton cycle or SCO 2 -based solar Rankine cycle system have been covered [6][7][8].SCO 2 cycle has been considered as the most promising and attractive candidate in concentrated solar power, coal power system, or bottoming cycle of fuel cells due to its abovementioned advantages [9][10][11].As per turbomachinery, SCO 2 is also a promising working fluid.In miniature power units, compact radial turbines and centrifugal compressors with complex geometries are designed and used [12].In large power units, axial turbines and compressors with SCO 2 as a working fluid can be obtained [13].
Centrifugal compressors of various power inputs with SCO 2 as a working fluid have drawn more and more attention lately.Both experimental and computational methods are applied.Monje B. et al. [10,14] conducted thorough research on SCO 2 systems and compared the aerodynamic performance of air and SCO 2 in conical diffusers.Budinis S. and Thornhill N.F.[15] studied the computer-based design and analysis of control systems for centrifugal compressors when the operating fluid is supercritical CO 2 .Utamura M. et al. [16] from Tokyo Institute of Technology tested the aerodynamic characteristics of a centrifugal compressor working in supercritical carbon dioxide by an experimental method.Computational fluid dynamics (CFD) has been widely used as a tool for designing, evaluating, and improving turbomachinery devices, including compressors [17,18], turbines [19,20], and pumps [21,22].Zhao H. et al. [23] carried out the numerical investigation of two-phase flow characteristics on the blade tip of a supercritical CO 2 centrifugal compressor.Behafarid F. and Podowski M.Z.[24] presented a novel modeling approach and the corresponding computer simulation of the SCO 2 flow inside a high-speed compact centrifugal compressor.Kim S.G. et al. [25] conducted a 3D numerical study for the full SCO 2 compressor geometry including diffuser and volute.The predictions were compared with measurements obtained from the constructed SCO 2 compressor test facility.It demonstrated the accuracy and the feasibility of numerical investigation.Compared to normal CO 2 , the centrifugal compressor with SCO 2 has higher efficiency, a more compact mechanical structure, and a larger pressure ratio [26].Hence, the work in this paper is of great importance for future power systems.Nevertheless, one potential disadvantage of SCO 2 lies in the phase change near the critical point, which can lead to the instability of actual operation.In this paper, the state parameters of inlet are above the critical point of CO 2 and the geometry of impeller and diffuser is well designed to avoid the subcritical state of the working fluid as much as possible [27].
In this paper, based on the real gas property of SCO 2 , a centrifugal compressor of large power and high pressure ratio, which has a large mass flow rate, has been designed.The corresponding geometrical parameters of the impeller and diffuser are obtained.Then, a three-dimensional model of the centrifugal compressor is established.With the application of the CFD method, the designed SCO 2 compressor is analyzed in steady and unsteady state.Hence, the detailed flow characteristics and operating performance are revealed.The objective of the numerical investigation is to measure the overall performance of this SCO 2 centrifugal compressor and to obtain the fundamental data for the design optimization.

Design of a Centrifugal Compressor
The state parameters of inlet and outlet were selected by taking an example from Zhao H. et al. [23].The pressure of the inlet was increased to guarantee that the SCO 2 fluid was always in supercritical condition.The design point was selected to optimize performance at the 100% power condition.The compressor design was developed using a three-step design process.First, a quick parametric screening study using meanline analysis and empiricism was conducted to synthesize the effects of various design parameters on performance [28].Once a skeletal geometry model that had a performance potential consistent with the requirements was established, the detailed geometry of the centrifugal compressor was sculpted using ANSYS Workbench as a turbomachinery geometry modeler.The development of the detailed geometry was guided by various aerodynamic analysis tools including 3D viscous CFD.

1D Calculations and Analyses
The 1D meanline model computed the overall compressor geometric size along with performance estimates.
Aside from meeting the required pressure rise and efficiency, the compressor must also operate stably over a broad range of conditions.A prime concern for compressor operability is the possibility of having localized condensed phase change of the working fluid at the inlet to the compressor.Any phase change of the working fluid SCO 2 will result in large variations in density which can adversely impact not only the performance of the compressor, but also its durability.To reduce the possibility of local acceleration of the fluid at the inlet, and possible condensed phase change, the inlet was sized to maximize the local static pressure similar to the traditional approach for sizing a pump to preclude cavitation.This produced a larger inlet area than required for minimum relative velocity at the shroud.In this case, the increase in impeller size did not have a significant effect on low-loss operating range since the inlet Mach number was relatively low.Considering all the difficulties mentioned above, the performance estimates as well as the state parameters at inlet and outlet of the designed centrifugal compressor are listed in Table 1.

Geometry of Impeller
The detailed geometric parameters of the impeller are determined in the design process.It is of significance to select a proper profile of the impeller.In this paper, multiple tests were conducted to ensure the high efficiency of the SCO 2 centrifugal compressor.To avoid the subcritical state at the inlet of the impeller, several CFD results of different impeller blade profiles were compared to select the one without transcritical phenomenon.This process guided the evolution of the impeller blade profile.We finally chose the profile of best performance after several improvements had been made.
The meridian plane of the impeller is presented in Figure 1a below, while the 3D impeller model is shown in Figure 1b.In Figure 1a, the blue arrow stands for low pressure and low temperature while the red arrow is just the reverse.
Appl.Sci.2017, 7, 310 3 of 17 reduce the possibility of local acceleration of the fluid at the inlet, and possible condensed phase change, the inlet was sized to maximize the local static pressure similar to the traditional approach for sizing a pump to preclude cavitation.This produced a larger inlet area than required for minimum relative velocity at the shroud.In this case, the increase in impeller size did not have a significant effect on low-loss operating range since the inlet Mach number was relatively low.Considering all the difficulties mentioned above, the performance estimates as well as the state parameters at inlet and outlet of the designed centrifugal compressor are listed in Table 1.18.9

Geometry of Impeller
The detailed geometric parameters of the impeller are determined in the design process.It is of significance to select a proper profile of the impeller.In this paper, multiple tests were conducted to ensure the high efficiency of the SCO2 centrifugal compressor.To avoid the subcritical state at the inlet of the impeller, several CFD results of different impeller blade profiles were compared to select the one without transcritical phenomenon.This process guided the evolution of the impeller blade profile.We finally chose the profile of best performance after several improvements had been made.
The meridian plane of the impeller is presented in Figure 1a below, while the 3D impeller model is shown in Figure 1b.In Figure 1a, the blue arrow stands for low pressure and low temperature while the red arrow is just the reverse.The impeller adopts twisted blades, which reveal preferable performance in most centrifugal compressors.The profile of hub and shroud as well as the blade lean angle distributions along hub to shroud are shown in Figure 2, respectively.The blade was modeled by stretching from hub profile to shroud profile, according to the blade lean angle distributions.The impeller adopts twisted blades, which reveal preferable performance in most centrifugal compressors.The profile of hub and shroud as well as the blade lean angle distributions along hub to shroud are shown in Figure 2, respectively.The blade was modeled by stretching from hub profile to shroud profile, according to the blade lean angle distributions.Leading edge (LE) and trailing edge (TE) of the blades adopt ellipse type and circular type, respectively.The LE elliptic ratio is 2.0.These types can effectively reduce the impact at LE and TE, thus decreasing the aerodynamic loss.
The specific geometric parameters of the impeller meridian plane and the twisted blades are given in Table 2. To guarantee the sealing performance, a shrouded impeller was adopted.

Geometry of Diffuser
The detailed geometric parameters of the diffuser were also determined in the design process.In the diffuser, the velocity of working fluid SCO2 decrease, and the pressure can further increase.Hence, a well-designed diffuser can significantly raise the efficiency of centrifugal compressors.Likewise, multiple tests have been conducted to ensure the validity of the diffuser.To avoid the circulation reflux and the vortex in the wake flow of the diffuser, several CFD results of different diffuser blade profiles were compared to select the one without apparent vortex and backflow.Finally, we chose the profile of best aerodynamic performance and least vortex after several rounds of improvement, which ensured the least flow loss in the diffuser.Leading edge (LE) and trailing edge (TE) of the blades adopt ellipse type and circular type, respectively.The LE elliptic ratio is 2.0.These types can effectively reduce the impact at LE and TE, thus decreasing the aerodynamic loss.
The specific geometric parameters of the impeller meridian plane and the twisted blades are given in Table 2. To guarantee the sealing performance, a shrouded impeller was adopted.

Geometry of Diffuser
The detailed geometric parameters of the diffuser were also determined in the design process.In the diffuser, the velocity of working fluid SCO 2 decrease, and the pressure can further increase.Hence, a well-designed diffuser can significantly raise the efficiency of centrifugal compressors.Likewise, multiple tests have been conducted to ensure the validity of the diffuser.To avoid the circulation reflux and the vortex in the wake flow of the diffuser, several CFD results of different diffuser blade profiles were compared to select the one without apparent vortex and backflow.Finally, we chose the profile of best aerodynamic performance and least vortex after several rounds of improvement, which ensured the least flow loss in the diffuser.
The meridian plane of the diffuser is presented in Figure 3a below, while the profile of the diffuser is shown in Figure 3b.In this paper, we employed a single type of straight blade (i.e., the profile of hub and shroud are identical).In Figure 3a, the blue arrow stands for low pressure and low temperature while the red arrow is just the reverse.
Appl.Sci.2017, 7, 310 5 of 17 The meridian plane of the diffuser is presented in Figure 3a below, while the profile of the diffuser is shown in Figure 3b.In this paper, we employed a single type of straight blade (i.e., the profile of hub and shroud are identical).In Figure 3a, the blue arrow stands for low pressure and low temperature while the red arrow is just the reverse.The LE and TE of the diffuser blade adopted ellipse type and circular type, respectively.The LE elliptic ratio was 2.0.These types can effectively reduce the impact at LE and TE when the working fluid SCO2 flow passes the diffuser channel.
In order to obtain a convergence solution, the outlet of diffuser was extended.Furthermore, the gap between the impeller and diffuser was 8 mm, which was also tested and filtered to reduce the impact loss towards the diffuser.The 3D diffuser model is shown in Figure 4.The specific geometric parameters of the diffuser meridian plane and the straight blades are given in Table 3.

Computational Method
The 3D viscous CFD simulation in this paper was carried out in commercially-available software ANSYS-CFX.The LE and TE of the diffuser blade adopted ellipse type and circular type, respectively.The LE elliptic ratio was 2.0.These types can effectively reduce the impact at LE and TE when the working fluid SCO 2 flow passes the diffuser channel.
In order to obtain a convergence solution, the outlet of diffuser was extended.Furthermore, the gap between the impeller and diffuser was 8 mm, which was also tested and filtered to reduce the impact loss towards the diffuser.The 3D diffuser model is shown in Figure 4.The meridian plane of the diffuser is presented in Figure 3a below, while the profile of the diffuser is shown in Figure 3b.In this paper, we employed a single type of straight blade (i.e., the profile of hub and shroud are identical).In Figure 3a, the blue arrow stands for low pressure and low temperature while the red arrow is just the reverse.The LE and TE of the diffuser blade adopted ellipse type and circular type, respectively.The LE elliptic ratio was 2.0.These types can effectively reduce the impact at LE and TE when the working fluid SCO2 flow passes the diffuser channel.
In order to obtain a convergence solution, the outlet of diffuser was extended.Furthermore, the gap between the impeller and diffuser was 8 mm, which was also tested and filtered to reduce the impact loss towards the diffuser.The 3D diffuser model is shown in Figure 4.The specific geometric parameters of the diffuser meridian plane and the straight blades are given in Table 3.

Computational Method
The 3D viscous CFD simulation in this paper was carried out in commercially-available The specific geometric parameters of the diffuser meridian plane and the straight blades are given in Table 3.

Computational Method
The 3D viscous CFD simulation in this paper was carried out in commercially-available software ANSYS-CFX.

Governing Equations
The governing equations were derived from the conservation of mass, momentum, and energy.The general form of these equations can be written as [29]: The terms in the above equation are respectively transient, convective, diffusive, and source term.Where φ is the universal variable, and φ = 1 indicates continuity equation.When φ is temperature, it represents energy equation, and when φ is set as the velocity of an arbitrary direction, Equation ( 1) corresponds to being a momentum equation.In addition, Γ is generalized diffusion coefficient and S is generalized source term.
Turbulence modeling is a process to make Reynolds averaged equations closed.Compared to the standard k-ε turbulence model, RNG k-ε is demonstrated to be more accurate by Hu et al. [30], which appears high precision.Hence, the RNG k-ε turbulence model was employed for the investigation in this paper.The transport equations of turbulent kinetic energy k and turbulent dissipation rate ε respectively perform as follows [29]: ∂ ∂t where α k , α ε , C 1ε , C 2ε , and C 3ε are empirical constants.
The abovementioned G k is the production of turbulence kinetic energy caused by averaged velocity gradient, G b is the production of turbulence kinetic energy caused by buoyancy, and the dilatation dissipation term Y M can exert the compressibility effect.G k , G b , and Y M respectively perform [31]: where β is the coefficient of thermal expansion, Pr t is the turbulence Prandtl number, and M t is the turbulence Mach number.When analyzing the unsteady flow characteristics, the transient rotor stator method in ANSYS-CFX was applied to the interface between the impeller and the diffusers.Virtual and true time step were introduced to the time discretization in order to avoid the restriction of time step length in explicit scheme and the complication transform of implicit scheme [32].
where the time level is denoted by the superscript n, q is the symbol of backward difference operator, and R is the residual.
In order to describe the practical property of the working fluid SCO 2 in the designed centrifugal compressor, a modified Soave Redlich Kwong equation of state in CFX database MATERIALS-redkw was adopted [33].It has been proven to be precise and effective by Kim S.G. et al. [25].In the database, "CO2RK" represents the working fluid carbon dioxide.With the restriction of pressure and temperature, it is able to be limited in the supercritical region, thus validly simulating the property of SCO 2 in the designed compressor [34].The explicit state equation is shown below in Equation (8).
where V is the molar volume of SCO 2 , P and T are respectively pressure and temperature of fluid, while R represents gas constant.The correction factors perform as follows: where, it is worthwhile to mention that correction factor a is related to the fluid temperature T.
In Equation ( 10), the acentric factor ω is set as 0.225 in the computational process and T r is calculated as follows: In this way, this equation of state strikes a great balance between simplicity and accuracy when performing a real fluid computational simulation [35].Likewise, the Soave Redlich Kwong state equation is considered to possess reasonable accuracy, especially near and above the critical point, thus ensuring the validity of SCO 2 simulation [36,37].
Moreover, the calculation of specific heat capacity should be considered.The specific heat at 0 gage pressure is computed by a fourth order polynomial as below, which ensures high precision:

Discretization Method and Boundary Conditions
In this paper, we utilize Reynolds Averaged Navier-Stokes (RANS) equations to dispose turbulent flow.The 3D, steady, and unsteady RANS equations were discretized using finite volume scheme.For the determination of near-wall velocities, a scalable wall function was chosen.Total energy equation including the viscous work term was used for energy conservation.A high resolution advection scheme of CFX was used.
Boundary conditions are indispensable to solve the aforementioned governing equations.The boundary conditions employed in this paper are as follows: mass flow rate of inlet (273.1 kg/s), static temperature of inlet (306.7 K), and static pressure of outlet (18.9 MPa).Meanwhile, the rotating speed was 15,000 rpm.The time step was set as 1.3333 × 10 −5 s, which indicates that 20 steps were needed to pass a diffuser passage.It has been concluded by Liu Y.W. et al. that, in general, unsteady simulation is recommended to compute at least 10 to 20 time steps per period.To obtain a more accurate result, 20 steps per passage were chosen [38].Hence, 300 steps were able to cover a whole period of impeller revolution.

Generation of Mesh
A high quality three-dimensional computational mesh was necessary for the CFD simulation in this paper.In this case, we adopted hexahedral meshes in all fluid domains.As shown in Figure 5, O-type mesh was applied around impeller blades, which greatly improved the mesh quality at the leading edge and the trailing edge of the blades.Additionally, H-type mesh was mainly employed in the flow channels, which also enhanced the mesh quality.It can also be observed that the mesh in the boundary layers of the wall was densely generated, to be exact, the hub wall, the shroud wall, and the blade surface.Hence, the y + all over the impeller domain was between 20~110, which corresponds to the requirement of the RNG k-ε turbulence model.The mesh in this investigation excluded tip clearances since a shrouded rotor was designed.As shown in Figure 6, O-type mesh was also applied around diffuser blades and H-type mesh was employed in the flow channels.The mesh in the boundary layers of the wall was densely generated, and the y + all over the diffuser domain was between 30~80, which also corresponds to the requirement of the RNG k-ε turbulence model.As explained above, in order to be within the wall function validity limits, the first mesh line was placed at the y + range of 20~110 in near the wall region.To guarantee the accuracy of the numerical simulation, a mesh independence test was conducted before the numerical analysis of unsteady flow in the designed centrifugal compressor.In the independence test, we used steady simulation for convenience.Since steady result served as the initial value for unsteady computation, it was efficient to examine the mesh validity and obtain the initialization simultaneously.In the steady simulation, the boundary conditions and the rotating speed were the same as unsteady calculation.Additionally, the frozen rotor method in ANSYS-CFX was applied to the interface between the impeller and the diffusers.
Four additional grids were generated to examine mesh convergence criteria and solution dependence on the mesh resolution.The mesh sensitivity was examined in terms of the torque of a single blade.The result of mesh independence test is presented in Table 4.As shown in Figure 6, O-type mesh was also applied around diffuser blades and H-type mesh was employed in the flow channels.The mesh in the boundary layers of the wall was densely generated, and the y + all over the diffuser domain was between 30~80, which also corresponds to the requirement of the RNG k-ε turbulence model.As shown in Figure 6, O-type mesh was also applied around diffuser blades and H-type mesh was employed in the flow channels.The mesh in the boundary layers of the wall was densely generated, and the y + all over the diffuser domain was between 30~80, which also corresponds to the requirement of the RNG k-ε turbulence model.As explained above, in order to be within the wall function validity limits, the first mesh line was placed at the y + range of 20~110 in near the wall region.To guarantee the accuracy of the numerical simulation, a mesh independence test was conducted before the numerical analysis of unsteady flow in the designed centrifugal compressor.In the independence test, we used steady simulation for convenience.Since steady result served as the initial value for unsteady computation, it was efficient to examine the mesh validity and obtain the initialization simultaneously.In the steady simulation, the boundary conditions and the rotating speed were the same as unsteady calculation.Additionally, the frozen rotor method in ANSYS-CFX was applied to the interface between the impeller and the diffusers.
Four additional grids were generated to examine mesh convergence criteria and solution dependence on the mesh resolution.The mesh sensitivity was examined in terms of the torque of a single blade.The result of mesh independence test is presented in Table 4.As explained above, in order to be within the wall function validity limits, the first mesh line was placed at the y + range of 20~110 in near the wall region.To guarantee the accuracy of the numerical simulation, a mesh independence test was conducted before the numerical analysis of unsteady flow in the designed centrifugal compressor.In the independence test, we used steady simulation for convenience.Since steady result served as the initial value for unsteady computation, it was efficient to examine the mesh validity and obtain the initialization simultaneously.In the steady simulation, the boundary conditions and the rotating speed were the same as unsteady calculation.Additionally, the frozen rotor method in ANSYS-CFX was applied to the interface between the impeller and the diffusers.
Four additional grids were generated to examine mesh convergence criteria and solution dependence on the mesh resolution.The mesh sensitivity was examined in terms of the torque of a single blade.The result of mesh independence test is presented in Table 4.It can be concluded that when the element number exceeds approximate 599.6 × 10 4 , the torque of a single blade is of small variation.Hence, based on the balance of time cost and precision, we ultimately selected sequence number 4 to conduct the unsteady numerical simulation; in other words, the total element number was 599.6 × 10 4 .Elements in rotor part and stator part were respectively about 300 × 10 4 .

Results and Discussion
In this section, the result of compressor design and 3D viscous CFD analysis are provided and discussed.
Due to the influence of rotation, the flow in turbomachinery is inherently unsteady.Hence, it is more accurate to simulate the flow characteristics in an unsteady state.With the abovementioned boundary conditions, an unsteady numerical analysis was conducted based on the steady simulation result.The steady simulation was considered converged when RMS residuals were below 10 −4 and the monitored parameters remained stable.
In the unsteady calculation process, to monitor the convergence status, we chose to extract the mass flow rate at compressor outlet in each time step.Figure 7 shows the variation of the monitored parameters with time step.It can be concluded that when the element number exceeds approximate 599.6 × 10 4 , the torque of a single blade is of small variation.Hence, based on the balance of time cost and precision, we ultimately selected sequence number 4 to conduct the unsteady numerical simulation; in other words, the total element number was 599.6 × 10 4 .Elements in rotor part and stator part were respectively about 300 × 10 4 .

Results and Discussion
In this section, the result of compressor design and 3D viscous CFD analysis are provided and discussed.
Due to the influence of rotation, the flow in turbomachinery is inherently unsteady.Hence, it is more accurate to simulate the flow characteristics in an unsteady state.With the abovementioned boundary conditions, an unsteady numerical analysis was conducted based on the steady simulation result.The steady simulation was considered converged when RMS residuals were below 10 −4 and the monitored parameters remained stable.
In the unsteady calculation process, to monitor the convergence status, we chose to extract the mass flow rate at compressor outlet in each time step.Figure 7 shows the variation of the monitored parameters with time step.As Figure 5 indicates, in steady computation, the mass flow rate of compressor outlet was close to the inlet (i.e., 273.1 kg/s).The unsteady computation exhibited apparent periodical changes after approximate 400 time steps, and the averaged mass flow rate of outlet approached 273.1 kg/s.Hence, it can be concluded that the unsteady calculation reached convergence criterion.We chose one revolution of impeller as the analysis segment, as shown in Figure 7 (i.e., 300 time steps), to run the aftermentioned analyses in this paper.As Figure 5 indicates, in steady computation, the mass flow rate of compressor outlet was close to the inlet (i.e., 273.1 kg/s).The unsteady computation exhibited apparent periodical changes after approximate 400 time steps, and the averaged mass flow rate of outlet approached 273.1 kg/s.Hence, it can be concluded that the unsteady calculation reached convergence criterion.We chose one revolution of impeller as the analysis segment, as shown in Figure 7 (i.e., 300 time steps), to run the aftermentioned analyses in this paper.

Unsteady Flow Phenomenon in the Centrifugal Compressor
The rotational condition resulted in an extremely unsteady flow state in the centrifugal compressor.In this section, streamlines, pressure, Mach number distribution at different time steps in one period are covered to analyze the unsteady flow in the designed SCO 2 centrifugal compressor.Temperature distribution on the impeller blade and diffuser blade is also exhibited to show the influence of LE and TE towards SCO 2 flow.
We chose a period of revolution after the convergence of unsteady computation to conduct the following discussion.The results of four key moments in a rotational period T are compared to illuminate the unsteady flow phenomenon.In the SCO 2 centrifugal compressor, the SCO 2 fluid first flows past the impeller passage.The input power drives the impeller to rotate.On account of rotation, the centrifugal force acts on SCO 2 fluid, resulting in the increase of pressure.Afterwards, the working fluid flows into the diffuser passage, in which its kinetic energy transforms into the static pressure energy.Hence, the static pressure further enlarges until the SCO 2 fluid flows out of diffuser.
Figure 8 gives the streamline in the whole passage of the SCO 2 centrifugal compressor.

Unsteady Flow Phenomenon in the Centrifugal Compressor
The rotational condition resulted in an extremely unsteady flow state in the centrifugal compressor.In this section, streamlines, pressure, Mach number distribution at different time steps in one period are covered to analyze the unsteady flow in the designed SCO2 centrifugal compressor.Temperature distribution on the impeller blade and diffuser blade is also exhibited to show the influence of LE and TE towards SCO2 flow.
We chose a period of revolution after the convergence of unsteady computation to conduct the following discussion.The results of four key moments in a rotational period T are compared to illuminate the unsteady flow phenomenon.In the SCO2 centrifugal compressor, the SCO2 fluid first flows past the impeller passage.The input power drives the impeller to rotate.On account of rotation, the centrifugal force acts on SCO2 fluid, resulting in the increase of pressure.Afterwards, the working fluid flows into the diffuser passage, in which its kinetic energy transforms into the static pressure energy.Hence, the static pressure further enlarges until the SCO2 fluid flows out of diffuser.
Figure 8 gives the streamline in the whole passage of the SCO2 centrifugal compressor.As shown in Figure 8, due to the influence of rotation, the diffusers are situated at different locations towards the impeller in each time step.The inlet of centrifugal compressor first serves as an intake of SCO2 fluid.In the impeller passage, the working fluid rotates along with the impeller, thus resulting in vortices in the passages.In the meantime, the relative velocity remains almost constant.Then, SCO2 enters into the diffuser passage.The absolute velocity increases with rotational speed taken into account.At the LE of diffuser blade, a velocity stagnation point appears.Afterwards, the working fluid passed through the diffuser passage with the increasing pressure and the decrease in velocity.Comparing the four moments in a period of revolution, the distribution of velocity streamlines are roughly identical while the distribution of vortices vary in different time steps of a period, respectively, as shown in Figure 8a-d.To be exact, at 0.5T and 0.75T, the diffuser blades hold back the wake flow from impeller, thus resulting in more vortices in the rotor passages.Especially at 0.75T, the diffuser blade locates at the middle of rotor passage.Due to As shown in Figure 8, due to the influence of rotation, the diffusers are situated at different locations towards the impeller in each time step.The inlet of centrifugal compressor first serves as an intake of SCO 2 fluid.In the impeller passage, the working fluid rotates along with the impeller, thus resulting in vortices in the passages.In the meantime, the relative velocity remains almost constant.Then, SCO 2 enters into the diffuser passage.The absolute velocity increases with rotational speed taken into account.At the LE of diffuser blade, a velocity stagnation point appears.Afterwards, the working fluid passed through the diffuser passage with the increasing pressure and the decrease in velocity.Comparing the four moments in a period of revolution, the distribution of velocity streamlines are roughly identical while the distribution of vortices vary in different time steps of a period, respectively, as shown in Figure 8a-d.To be exact, at 0.5T and 0.75T, the diffuser blades hold back the wake flow from impeller, thus resulting in more vortices in the rotor passages.Especially at 0.75T, the diffuser blade locates at the middle of rotor passage.Due to the obstacle influence of diffuser blade, the maximum velocity locates near the LE of diffuser blade.Hence, the streamlines at 0.75T are the most disorganized.
Figure 9 shows the pressure contours at mid-span of impeller and diffuser blades at different time steps.
the obstacle influence of diffuser blade, the maximum velocity locates near the LE of diffuser blade.Hence, the streamlines at 0.75T are the most disorganized.
Figure 9 shows the pressure contours at mid-span of impeller and diffuser blades at different time steps.As Figure 9 indicates, the pressure achieves the minimum value at the inlet of compress.On account of rotation, the centrifugal force acts on SCO2 fluid, resulting in the increase of pressure in the rotor passage.Afterwards, the working fluid flows into the diffuser passage, in which its kinetic energy transforms into the static pressure energy.Hence, the static pressure further enlarges until the SCO2 fluid flows out of diffuser.Comparing the four moments in a period of revolution, the distribution of pressure at the mid-span of blades varies.When the TE of impeller blade is close to the LE of diffuser blade, the pressure of the outlet reaches its maximum, as shown in Figure 9a.Moreover, there exists a stagnation point at the LE of diffuser blade, respectively shown in Figure 9a-d.To be exact, at 0.75T, the diffuser blade is located at the middle of rotor passage and the diffuser blades hold back the wake flow from the impeller, thus resulting in a maximum stagnation pressure of 23.29 MPa.
Figure 10   As Figure 9 indicates, the pressure achieves the minimum value at the inlet of compress.On account of rotation, the centrifugal force acts on SCO 2 fluid, resulting in the increase of pressure in the rotor passage.Afterwards, the working fluid flows into the diffuser passage, in which its kinetic energy transforms into the static pressure energy.Hence, the static pressure further enlarges until the SCO 2 fluid flows out of diffuser.Comparing the four moments in a period of revolution, the distribution of pressure at the mid-span of blades varies.When the TE of impeller blade is close to the LE of diffuser blade, the pressure of the outlet reaches its maximum, as shown in Figure 9a.Moreover, there exists a stagnation point at the LE of diffuser blade, respectively shown in Figure 9a-d.To be exact, at 0.75T, the diffuser blade is located at the middle of rotor passage and the diffuser blades hold back the wake flow from the impeller, thus resulting in a maximum stagnation pressure of 23.29 MPa.
Figure 10 provides the Mach number contours at the mid-span of the impeller and diffuser blades at different time steps.
the obstacle influence of diffuser blade, the maximum velocity locates near the LE of diffuser blade.Hence, the streamlines at 0.75T are the most disorganized.
Figure 9 shows the pressure contours at mid-span of impeller and diffuser blades at different time steps.As Figure 9 indicates, the pressure achieves the minimum value at the inlet of compress.On account of rotation, the centrifugal force acts on SCO2 fluid, resulting in the increase of pressure in the rotor passage.Afterwards, the working fluid flows into the diffuser passage, in which its kinetic energy transforms into the static pressure energy.Hence, the static pressure further enlarges until the SCO2 fluid flows out of diffuser.Comparing the four moments in a period of revolution, the distribution of pressure at the mid-span of blades varies.When the TE of impeller blade is close to the LE of diffuser blade, the pressure of the outlet reaches its maximum, as shown in Figure 9a.Moreover, there exists a stagnation point at the LE of diffuser blade, respectively shown in Figure 9a-d.To be exact, at 0.75T, the diffuser blade is located at the middle of rotor passage and the diffuser blades hold back the wake flow from the impeller, thus resulting in a maximum stagnation pressure of 23.29 MPa.
Figure 10 provides the Mach number contours at the mid-span of the impeller and diffuser blades at different time steps.The Mach number in Figure 10 was calculated in stationary reference frame.As Figure 10 indicates, the SCO2 fluid at the mid-span of the blades accelerates in the impeller passage and reaches the maximum value at the outlet of impeller.In the diffuser passage, the working fluid gradually slows down.Meanwhile, it can be demonstrated that the maximum Mach number in the unsteady computation locates near the LE of the diffuser blade in Figure 10c (i.e., 0.75T).At 0.75T, the diffuser blade is located at the middle of the rotor passage, thus resulting in a disorganized flow condition.Due to the obstacle influence of the diffuser blade, the maximum Mach number reaches 0.58.
Figure 11 gives the temperature distributions of impeller and diffuser blades of the SCO2 centrifugal compressor.As shown in Figure 11, the distribution of blade temperature varies in different time steps.The working fluid SCO2 first achieves the minimum value at inlet of the compressor.With the pressure increases in the impeller passage, the temperature on the impeller blade also increases.Afterwards, the working fluid flows into the diffuser passage, and the temperature on the diffuser blade further increases.In comparing the four moments in a period of revolution, the maximum temperature in the unsteady computation is located at the LE of diffuser blade in Figure 11c.When the diffuser blade is located at the middle of the rotor passage, the stagnation point at the LE of the The Mach number in Figure 10 was calculated in stationary reference frame.As Figure 10 indicates, the SCO 2 fluid at the mid-span of the blades accelerates in the impeller passage and reaches the maximum value at the outlet of impeller.In the diffuser passage, the working fluid gradually slows down.Meanwhile, it can be demonstrated that the maximum Mach number in the unsteady computation locates near the LE of the diffuser blade in Figure 10c (i.e., 0.75T).At 0.75T, the diffuser blade is located at the middle of the rotor passage, thus resulting in a disorganized flow condition.Due to the obstacle influence of the diffuser blade, the maximum Mach number reaches 0.58.
Figure 11 gives the temperature distributions of impeller and diffuser blades of the SCO 2 centrifugal compressor.The Mach number in Figure 10 was calculated in stationary reference frame.As Figure 10 indicates, the SCO2 fluid at the mid-span of the blades accelerates in the impeller passage and reaches the maximum value at the outlet of impeller.In the diffuser passage, the working fluid gradually slows down.Meanwhile, it can be demonstrated that the maximum Mach number in the unsteady computation locates near the LE of the diffuser blade in Figure 10c (i.e., 0.75T).At 0.75T, the diffuser blade is located at the middle of the rotor passage, thus resulting in a disorganized flow condition.Due to the obstacle influence of the diffuser blade, the maximum Mach number reaches 0.58.
Figure 11 gives the temperature distributions of impeller and diffuser blades of the SCO2 centrifugal compressor.As shown in Figure 11, the distribution of blade temperature varies in different time steps.The working fluid SCO2 first achieves the minimum value at inlet of the compressor.With the pressure increases in the impeller passage, the temperature on the impeller blade also increases.Afterwards, the working fluid flows into the diffuser passage, and the temperature on the diffuser blade further increases.In comparing the four moments in a period of revolution, the maximum temperature in the unsteady computation is located at the LE of diffuser blade in Figure 11c.When the diffuser blade is located at the middle of the rotor passage, the stagnation point at the LE of the As shown in Figure 11, the distribution of blade temperature varies in different time steps.The working fluid SCO 2 first achieves the minimum value at inlet of the compressor.With the pressure increases in the impeller passage, the temperature on the impeller blade also increases.Afterwards, the working fluid flows into the diffuser passage, and the temperature on the diffuser blade further increases.In comparing the four moments in a period of revolution, the maximum temperature in the unsteady computation is located at the LE of diffuser blade in Figure 11c.When the diffuser blade is located at the middle of the rotor passage, the stagnation point at the LE of the diffuser blade is the most remarkable.The highest temperature 335.1 K in Figure 11c also indicates maximum aerodynamic loss at 0.75T.

Unsteady Forces on a Rotor Blade
The periodic flow exciting force on a rotor blade caused by the unsteady working condition is vital when studying the vibration of the impeller and the reliability of the entire rotor.
The axial and circumferential forces on a single rotor blade in one rotation period are presented in Figure 12.In the centrifugal compressor, the rotor blade is exposed to the violent unsteady loading due to the periodically high turbulence and energy dissipation flow generating in the cascade of the diffuser's leading edge.Figure 12 reveals obvious periodic characteristic in each of the 20 time steps (i.e., every time a rotor blade passes a diffuser flow passage).
Appl.Sci.2017, 7, 310 13 of 17 diffuser blade is the most remarkable.The highest temperature 335.1 K in Figure 11c also indicates maximum aerodynamic loss at 0.75T.

Unsteady Forces on a Rotor Blade
The periodic flow exciting force on a rotor blade caused by the unsteady working condition is vital when studying the vibration of the impeller and the reliability of the entire rotor.
The axial and circumferential forces on a single rotor blade in one rotation period are presented in Figure 12.In the centrifugal compressor, the rotor blade is exposed to the violent unsteady loading due to the periodically high turbulence and energy dissipation flow generating in the cascade of the diffuser's leading edge.Figure 12 reveals obvious periodic characteristic in each of the 20 time steps (i.e., every time a rotor blade passes a diffuser flow passage).The Fast Fourier Transform (FFT) is usually used to investigate frequency components of exciting force.The spectral function is defined as below: where k is the order of exciting force.
To study the distribution of axial and circumferential forces in time domain, FFT is employed in this paper.The spectrum analyses of the circumferential and the axial forces are conducted.
Hence, the spectrum analyses of circumferential and axial exciting force on a rotor blade are shown in Figure 13, in which the horizontal axis represents the multiple frequencies of the sampled signal and the vertical axis is the corresponding amplitude.The Fast Fourier Transform (FFT) is usually used to investigate frequency components of exciting force.The spectral function is defined as below: where k is the order of exciting force.
To study the distribution of axial and circumferential forces in time domain, FFT is employed in this paper.The spectrum analyses of the circumferential and the axial forces are conducted.
Hence, the spectrum analyses of circumferential and axial exciting force on a rotor blade are shown in Figure 13, in which the horizontal axis represents the multiple frequencies of the sampled signal and the vertical axis is the corresponding amplitude.

Unsteady Forces on a Rotor Blade
The periodic flow exciting force on a rotor blade caused by the unsteady working condition is vital when studying the vibration of the impeller and the reliability of the entire rotor.
The axial and circumferential forces on a single rotor blade in one rotation period are presented in Figure 12.In the centrifugal compressor, the rotor blade is exposed to the violent unsteady loading due to the periodically high turbulence and energy dissipation flow generating in the cascade of the diffuser's leading edge.Figure 12 reveals obvious periodic characteristic in each of the 20 time steps (i.e., every time a rotor blade passes a diffuser flow passage).The Fast Fourier Transform (FFT) is usually used to investigate frequency components of exciting force.The spectral function is defined as below: where k is the order of exciting force.
To study the distribution of axial and circumferential forces in time domain, FFT is employed in this paper.The spectrum analyses of the circumferential and the axial forces are conducted.
Hence, the spectrum analyses of circumferential and axial exciting force on a rotor blade are shown in Figure 13, in which the horizontal axis represents the multiple frequencies of the sampled signal and the vertical axis is the corresponding amplitude.The largest amplitude of exciting forces lies in the frequency of 3750 Hz.Meanwhile, the exact value of circumferential force and axial force are respectively 160.2 N and 59.3 N. The exciting-vibration force factor is defined as follows: where F max is the maximum exciting force and the denominator represents the average exciting force.The maximum exciting-vibration force factors are respectively 0.06 for circumferential force and 0.08 for axial force, which can guarantee the safety of the compressor operation.The exciting force of low frequency is derived from the interference between stator and rotor.The blocking effect of the diffusers results in the largest exciting force.Additionally, other frequencies of large amplitude correspond to a certain multiple of rotational frequency.For instance, the second largest amplitude 47.6 N and 6.1 N (respectively circumferential and axial) correspond to 7500 Hz, which results from the nonuniformity at the TE of rotor blades.In addition, in the spectrum figure, the exciting force magnitude of high frequency is much smaller than that of low frequency, which suggests that the low frequency should be considered in a vibration analysis.

Power and Efficiency of the Designed Turbine
The process of analyzing the power and efficiency of the designed centrifugal compressor from the numerical simulation result is able to estimate whether the SCO 2 centrifugal compressor satisfies the design requirement.In this section, we compare the torque of steady result and unsteady time-averaged result.Meanwhile, the input power, pressure ratio, and isentropic efficiency are calculated and compared with the design condition.
The comparison of the single blade torque is performed in Figure 14.
The largest amplitude of exciting forces lies in the frequency of 3750 Hz.Meanwhile, the exact value of circumferential force and axial force are respectively 160.2 N and 59.3 N. The exciting-vibration force factor is defined as follows: where Fmax is the maximum exciting force and the denominator represents the average exciting force.
The maximum exciting-vibration force factors are respectively 0.06 for circumferential force and 0.08 for axial force, which can guarantee the safety of the compressor operation.The exciting force of low frequency is derived from the interference between stator and rotor.The blocking effect of the diffusers results in the largest exciting force.Additionally, other frequencies of large amplitude correspond to a certain multiple of rotational frequency.For instance, the second largest amplitude 47.6 N and 6.1 N (respectively circumferential and axial) correspond to 7500 Hz, which results from the nonuniformity at the TE of rotor blades.In addition, in the spectrum figure, the exciting force magnitude of high frequency is much smaller than that of low frequency, which suggests that the low frequency should be considered in a vibration analysis.

Power and Efficiency of the Designed Turbine
The process of analyzing the power and efficiency of the designed centrifugal compressor from the numerical simulation result is able to estimate whether the SCO2 centrifugal compressor satisfies the design requirement.In this section, we compare the torque of steady result and unsteady time-averaged result.Meanwhile, the input power, pressure ratio, and isentropic efficiency are calculated and compared with the design condition.
The comparison of the single blade torque is performed in Figure 14.The torque of a single blade in unsteady simulation has a periodic variation each time that the impeller rotates past a diffuser passage.The time-averaged unsteady result is 212.92N•m, which is slightly lower than the steady result 216.52 N•m.Since the unsteady computation can more accurately simulate the aerodynamic performance, the value 212.92 N•m is more reliable.
The pressure ratio, input power, and isentropic efficiency are significant parameters when measuring the performance of a designed centrifugal compressor.Isentropic efficiency is defined as follows: where h represents enthalpy, and the subscript is stands for isentropic parameter.The subscript 1 and 2 respectively represent the inlet and outlet of compressor.The pressure ratio, input power, and isentropic efficiency are significant parameters when measuring the performance of a designed centrifugal compressor.Isentropic efficiency is defined as follows: where h represents enthalpy, and the subscript is stands for isentropic parameter.The subscript 1 and 2 respectively represent the inlet and outlet of compressor.
In Table 5, we calculate the input power based on the value of torque mentioned above.The pressure ratio, input power and isentropic efficiency of design condition, steady CFD result and unsteady CFD result are compared.It can be concluded from Table 5 that the input power meets the design requirement.Nevertheless, the isentropic efficiency resulted from the CFD process is inferior to that of the design.The reason lies in the inaccurate computation of flow loss in the 1D viscous design.The profiles of the diffuser and the impeller have been ignored while the loss coefficient of leakage and impeller friction is mainly from empirical values.The 3D numerical investigation can more precisely simulate the flow process in this SCO 2 centrifugal compressor, especially in an unsteady simulation which has taken the influence of rotation into account.

Conclusions
In this paper, a SCO 2 centrifugal compressor of large input power and mass flow rate was designed and numerically investigated.A thorough numerical analysis of the unsteady flow field in the centrifugal compressor was performed in ANSYS-CFX.Streamlines, temperature, pressure, and Mach number distributions at different time steps in one revolution period were covered to present the unsteady effect of turbomachinery.Meanwhile, the force on a single rotor blade was monitored to investigate the frequency components of exciting force, thus providing foundation for vibration analysis.Moreover, the torque, output power, pressure ratio, and isentropic efficiency in the steady and the unsteady time-averaged condition were calculated and compared with the design condition to measure the validity of the design.
The preliminary design provided the key geometry parameters of impeller and diffuser.For instance, the inlet diameter of diffuser was 220 mm.In the meantime, multiple tests were conducted to ensure the validity of the impeller and the diffuser.We finally chose the geometry of best performance and least vortex after several rounds of improvements.
A thorough numerical analysis of the unsteady flow field in the centrifugal compressor was performed in ANSYS-CFX.The computation adopted hexahedral mesh, finite volume method, and the RNG k-ε two-equation turbulence model.In the unsteady computation, the monitored parameters revealed obvious periodic characteristics.Comparing the calculation result of four different time steps in a revolution, the result revealed obvious unsteady flow characteristics.To be specific, the flow at 0.75T had the worst flow characteristic while the flow at 0.25T was the best due to the different location of the diffuser towards the impeller.The largest amplitude of exciting forces was found in the frequency of 3750 Hz and the exact value of circumferential and axial forces were respectively 160.2 N and 59.3 N. The pressure ratio, input power, and isentropic efficiency mainly met the design requirement, although the efficiency was slightly lower in the viscous CFD computation owing to the inaccurate calculation of flow loss in the design.
In conclusion, the unsteady flow characteristics are prominent in the centrifugal compressor and unsteady CFD simulation is able to provide a more accurate estimation of compressor performance.On the whole, the design of the 5 MW centrifugal compressor in this paper can achieve a satisfactory aerodynamic performance.In the actual operation, the air exciting-vibration force caused by unsteady effect should be considered in the strength assessment.Also, it has been concluded that this design

Figure 7 .
Figure 7. Variation of the monitored parameters.

Figure 7 .
Figure 7. Variation of the monitored parameters.
provides the Mach number contours at the mid-span of the impeller and diffuser blades at different time steps.(a) (b)

Figure 12 .
Figure 12.Unsteady forces on a rotor blade.

Figure 12 .
Figure 12.Unsteady forces on a rotor blade.
the most remarkable.The highest temperature 335.1 K in Figure11calso indicates maximum aerodynamic loss at 0.75T.

Figure 12 .
Figure 12.Unsteady forces on a rotor blade.

Figure 14 .
Figure 14.Torque of a single blade in the SCO2 compressor.

Figure 14 .
Figure 14.Torque of a single blade in the SCO 2 compressor.

Table 1 .
Performance estimates and state parameters of the designed compressor.

Table 1 .
Performance estimates and state parameters of the designed compressor.

Table 4 .
Results of mesh independence test.

Table 4 .
Results of mesh independence test.

Table 5 .
Comparison of design and CFD results.