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

#### **Dongbo Shi and Yonghui Xie \***

State Key Laboratory for Strength and Vibration of Mechanical Structures, School of Energy and Power Engineering, Xi'an Jiaotong University, Xi'an 710049, China; sdb\_xjtu@126.com **\*** Correspondence: yhxie@mail.xjtu.edu.cn; Tel.: +86-029-82664443

Received: 22 February 2020; Accepted: 16 March 2020; Published: 19 March 2020

#### **Featured Application: In this research, a design-optimization method of a supercritical carbon dioxide compressor with high performance without a high speed requirement is proposed, which can improve the economy and reliability of the application system.**

**Abstract:** 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.

**Keywords:** supercritical carbon dioxide; centrifugal compressor; aerodynamic optimization design; numerical simulation

#### **1. 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. CO2 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-CO2 Brayton cycle is above the critical temperature of carbon dioxide. Additionally, S-CO2 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 CO2 can reduce compressor power consumption and improve cycle efficiency. Therefore, S-CO2 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.

At present, the SNL (Sandia National Laboratories), TIT (Tokyo Institute of Technology) and KIER (Korea Institute of Energy Research) have designed and completed the S-CO2 compressor equipment used in the experimental system. Table 1 summarizes the key parameters of the S-CO2 centrifugal compressor in the test system of the three research institutions [5–8].


**Table 1.** The key parameters of the S-CO2 centrifugal compressor.

It can be seen from Table 1 that, due to the high density characteristics of S-CO2, the centrifugal compressors designed by various research institutions have a smaller diameter, compact structure and higher speed. At present, the design of S-CO2 compressor is not perfect. The design of the diffuser, bearing, seal and other auxiliary parts is difficult. Therefore, the S-CO2 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-CO2 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-CO2 centrifugal compressor has been increasingly studied by scholars of various countries. The design of the S-CO2 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-CO2 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-CO2 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-CO2 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-CO2 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-CO2. 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-CO2 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.

#### **2. Design of the S-CO2 Compressor**

#### *2.1. Design Method*

Based on our previous design experience and methods [14,15], the thermodynamic design software of S-CO2 centrifugal compressor (S-CO2CPTD) 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:

$$\mathfrak{gl}\_{\rm tl} = \frac{1}{\mathcal{S}} (c\_{2\iota}\mu\_2 - c\_{1\iota}\mu\_1) \tag{1}$$

where *hth* (J/kg) is the energy head gained by one kilogram of fluid in the impeller blade channel, *g* (m/s2) 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.

$$h\_{\text{tot}} = h\_{\text{th}} + h\_{df} + h\_t \tag{2}$$

where *htot* (J/kg) is the total energy head, *hdf* (J/kg) is the loss of impeller resistance and *ht* (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:

$$h\_{tot} = \frac{c\_p}{A}(T\_2 - T\_1) + \frac{c\_2^2 + c\_1^2}{2g} \tag{3}$$

where *cp* (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:

$$h\_{th} = \int\_{1}^{2} \frac{dp}{\chi} + \frac{c\_2^2 - c\_1^2}{2g} + h\_{f1} \tag{4}$$

where *p* (Pa) is the pressure, γ (N/m3) is the specific weight and *hf*<sup>1</sup> (J/kg) is the flow loss of fluid flowing in the impeller.

When adopting the Bernoulli equation in the diffuser, *hth* is zero:

$$\frac{c\_2^2 - c\_3^2}{2g} = \int\_2^3 \frac{dp}{\gamma} + h\_{f2} \tag{5}$$

where *hf*<sup>2</sup> (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:

$$D = \sqrt{\frac{\dot{m}}{\rho\_2 \mu\_{2r} \pi \tau\_1 \tau\_2}}\tag{6}$$

where . *m* (kg/s) is the mass flow rate, ρ (kg/m3) is density, τ<sup>1</sup> is the blocking factor of the impeller outlet, τ<sup>2</sup> = *l*2/*D* is the relative width of the impeller outlet and *l*<sup>2</sup> (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:

$$l\_1 = \frac{(1+\tau\_4) \cdot q\_1}{\pi \tau\_3 D' u\_{1r}} \tag{7}$$

where τ<sup>4</sup> is the leakage loss coefficient of the impeller cover, *q* (m3/s) is the volume flow rate, τ<sup>3</sup> 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.

**Figure 1.** Geometric design diagram of the impeller blade.

#### *2.2. 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.


**Table 2.** Key parameters of the initial design of the compressor.

**Figure 2.** The initial design of the S-CO2 centrifugal compressor: (**a**) rotor's meridian face profile and (**b**) 3-D model.

#### *2.3. 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 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\_r(\mathbf{y}) = f\_a(\mathbf{y}) + \Delta f\_i(\mathbf{y}) \tag{8}$$

where Δ*fi*(**y**) is the deviation between the thermal design *fa*(**y**) and CFD accurate calculation results. The calculation deviation satisfies Gauss distribution at each point in the whole design space:

$$
\Delta f(\mathbf{y}) \sim \mathcal{N}(\mu(\mathbf{y}), \sigma^2(\mathbf{y})) \tag{9}
$$

The mean μ(**y**) and variance σ(**y**) of Gaussian distribution are:

$$\begin{cases} \mu(\mathbf{y}) = \mu + \mathbf{R}\_{\mathbf{y}D} \mathbf{R}\_{DD}^{-1} (\Delta f \mathbf{D} - 1 \mu) \\ \sigma(\mathbf{y}) = \sigma^2 (1 - \mathbf{R}\_{\mathbf{y}D} \mathbf{R}^{-1} \_{DD} \mathbf{R}\_{D\mathbf{y}}) \end{cases} \tag{10}$$

where *D* is the training set, **Ry***<sup>D</sup>* and **R***D***<sup>y</sup>** are the covariance vector of the calculated condition and the new sampling point, **R***DD* is the covariance matrix and Δ*fD* 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:

$$\begin{cases} \widecheck{f}\_r(\mathbf{y}) = f\_a(\mathbf{y}) + \mu(\mathbf{y}) + \xi \sigma(\mathbf{y})\\ \widecheck{f}\_r(\mathbf{y}) = f\_b(\mathbf{y}) + \mu(\mathbf{y}) - \xi \sigma(\mathbf{y}) \end{cases} \tag{11}$$

where *<sup>f</sup> <sup>r</sup>*(**y**) is the lower limit estimation, *f <sup>r</sup>*(**y**) is the upper limit estimation and ξ is the confidence constant.

**Figure 3.** The design-optimization method based on the Gauss process regression.

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.

**Figure 4.** The optimization process.

#### *2.4. 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.

**Figure 5.** The impeller of the optimal design of the S-CO2 centrifugal compressor: (**a**) rotor's meridian face profile and (**b**) 3-D model of the impeller.

**Figure 6.** The diffuser of the optimal design of the S-CO2 centrifugal compressor: (**a**) blade profile on the B2B cross section and (**b**) 3-D model of the diffuser.

#### **3. Numerical Methods**

#### *3.1. 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–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 m2/s2) and epsilon (30,000 m2/s3) 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 <sup>×</sup> 10-6 and the iterative deviation less than 0.1% of the outlet temperature and impeller blade torque were considered as the convergence conditions.

#### *3.2. Mesh and Grid Independence*

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<sup>+</sup> near the wall was about 1, which met the calculation requirements of the SST turbulence model.

**Figure 7.** The grid schematic diagram of the impeller–diffuser fluid domain.

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.

**Figure 8.** The grid independence.

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.

$$
\Delta h\_s = c\_p T\_{\text{inlet}} \left[ 1 - \left( \frac{p\_{\text{out}}}{p\_{\text{inlet}}} \right)^{\frac{\kappa - 1}{\kappa}} \right] \tag{12}
$$

$$\eta\_{\rm fs} = \frac{P}{\dot{m}(\Delta \text{l}\_{\rm s} + c\_{\rm inlet}^2/2)}\tag{13}$$

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.

#### *3.3. 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 9.** The numerical validation.

#### **4. Results and Discussion**

#### *4.1. The Initial Design*

Figure 10 shows the velocity vector and partial enlargement of the 50% blade height section of the compressor initial design.

**Figure 10.** The velocity vector distribution of the 50% blade height section of the 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 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.

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 low-pressure area, and the working fluid would enter the two-phase region. This is consistent with the conclusions of other scholars [23,24].

**Figure 11.** The pressure distribution on the surface of impeller blades.

**Figure 12.** The cross critical region of the blade surface.

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-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.


**Table 5.** Performance analysis of the initial design compressor.

#### *4.2. The Optimal Design*

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 13.** The velocity vector distribution of the 50% blade height section for the optimal design of the compressor.

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.

**Figure 14.** Second flow velocity on S3 sections of the optimal compressor impeller: (**a**) 30% axial chord length; (**b**) 40% axial chord length and (**c**) 90% axial chord length.

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.

**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.** The pressure distribution of the blade surface.

**Figure 17.** The pressure curves of the blade surface with different blade high cross sections.

**Figure 18.** The cross critical region of the blade surface.

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.

**Figure 19.** The entropy distribution of different sections: (**a**) different blade height sections; (**b**) impeller outlet and (**c**) diffuser outlet.

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.

**Figure 20.** The compressor mass flow rate-total pressure ratio curve.

**Figure 21.** The compressor mass flow rate-total static efficiency curve.

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.


**Table 6.** Comparison of compressor performances.

#### **5. Conclusions**

In this paper, a thermodynamic design software of the S-CO2 centrifugal compressor (S-CO2CPTD) was developed based on the 1-D thermal design method. The 150 kW S-CO2 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:


**Author Contributions:** Conceptualization, D.S. and Y.X.; investigation, D.S.; methodology, D.S.; resources, Y.X.; software, D.S. and Y.X.; supervision, Y.X.; validation Y.X.; writing—original draft preparation, D.S.; writing—review and editing, D.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Design and Performance Analysis of a Supercritical Carbon Dioxide Heat Exchanger**

#### **Han Seo 1,\*, Jae Eun Cha 1, Jaemin Kim 2, Injin Sah <sup>1</sup> and Yong-Wan Kim <sup>1</sup>**


Received: 29 May 2020; Accepted: 23 June 2020; Published: 30 June 2020

**Abstract:** This paper presents a preliminary design and performance analysis of a supercritical CO2 (SCO2) heat exchanger for an SCO2 power generation system. The purpose of designing a SCO2 heat exchanger is to provide a high-temperature and high-pressure heat exchange core technology for advanced SCO2 power generation systems. The target outlet temperature and pressure for the SCO2 heat exchanger were 600 ◦C and 200 bar, respectively. A tubular type with a staggered tube bundle was selected as the SCO2 heat exchanger, and liquefied petroleum gas (LPG) and air were selected as heat sources. The design of the heat exchanger was based on the material selection and available tube specification. Preliminary performance evaluation of the SCO2 heat exchanger was conducted using an in-house code, and three-dimensional flow and thermal stress analysis were performed to verify the tube's integrity. The simulation results showed that the tubular type heat exchanger can endure high-temperature and high-pressure conditions under an SCO2 environment.

**Keywords:** supercritical CO2; heat exchanger; flow analysis; thermal stress analysis

#### **1. Introduction**

The supercritical carbon dioxide (SCO2) Brayton cycle has been considered one of the most promising alternatives to existing power generation systems, such as the steam Rankine and gas Brayton cycles. The steam Rankine cycle can achieve high thermal efficiency due to the low pumping power, but the overall size of the system is large because the steam density of the low pressure side is lower than the atmospheric pressure. The turbine inlet temperature (TIT) of the gas Brayton cycle is higher than in the steam Rankine cycle. It can achieve a high thermal power output, but material integrity and high compression work problems still remain. The SCO2 cycle has the advantages of both the steam Rankine and gas Brayton cycles because of its fluid characteristics. First of all, the supercritical region of the SCO2 is readily accessible (Tc = 31.1 ◦C, Pc = 73.8 bar); thus, the system can be controlled easily compared to other critical state fluids. Because the SCO2 cycle operates near the critical point, the fluid reflects both liquid and gas properties. The compression work of the SCO2 consumes less than the conventional gas Brayton cycle. In addition, the density difference between the hot and cold sides are small. This means that the overall size of the SCO2 power generation system can be minimized compared to the conventional power generation systems. Based on the advantages, the SCO2 cycle can be applied to various technologies, such as concentrating power, fossil fuel, geothermal, nuclear, ship-board propulsion, waste heat recovery, etc. Therefore, the development of the SCO2 power cycle has been extensively studied.

Studies on the SCO2 cycle have taken place since the 1960s, but its development did not occur due to these technological advances [1,2]. In the 2000s, Dostal et al. [3] studied the SCO2 cycle as the next power generation system for Generation IV nuclear reactors because the TIT for advanced nuclear reactors is about 550 ◦C, a system operating condition which is between the steam Rankine and gas Brayton cycles. Because the SCO2 cycle can be designed in a wide range of the TIT, the SCO2 system has received more attention as the next power generation system. In addition, the SCO2 cycle can also be used in various heat sources, such as fossil fuel, concentrating solar power, shipboard propulsion, waste heat recovery, and geothermal.

After validating the SCO2 cycle as a new power generation system, numerous studies have been conducted to verify the cycle. Sandia National Laboratories designed a recompression configuration, and various experimental facilities have been developed, such as turbine and compressor performance characterizations [4–6]. The Naval Nuclear Laboratory designed and built an integrated system test (IST) for the SCO2 Brayton cycle. The IST is a simple recuperated Brayton cycle for variable turbomachinery tests. The system demonstrated that the SCO2 Brayton cycle was controlled throughout the entire system operation, but inherent problems related to the SCO2 Brayton cycle were identified [7–9]. Ecogen Power Systems designed a SCO2 power cycle (EPS-100) for exhaust heat recovery applications. The EPS-100 is the first commercial-scale SCO2 system, and it has a 7MWe scale with multiple stages of recuperation and extraction from the heat source. The validation test of the EPS-100 was completed, and commercialization took place [10,11]. The SunShot program, which develops the SCO2 Brayton cycle for a concentrating solar power (CSP) system, was initiated with a simple recuperated cycle [12]. The target demonstration of the SCO2 system is 10 MWe with a 50% net thermal efficiency.

To demonstrate the major components in the SCO2 power generation system, the SCO2's heat exchanger and turbo-expander have been tested in a 1 MWe test loop [13,14]. A tubular-type heat exchanger was selected as the primary heater, which was connected with a commercial natural gas-fired combustor [13]. The SunShot program is the first MW-scale SCO2 power cycle demonstration in which the TIT is higher than 700 ◦C. A demonstration of the SCO2 power cycle was performed by considering the unique characteristics of the CSP system [14].

In addition, the US Department of Energy (DOE) designated SCO2 research as a cross-cutting technology and supercritical transformational electric power (STEP) program with a collaboration between fossil, nuclear, and renewable energy [15]. The program has focused on designing, constructing, commissioning, and operating a 10-MWe SCO2 pilot plant test facility. The detailed design of the facility and equipment is now proceeding. Fabrication and construction of a pilot test facility with a simple-cycle test will be finished at the end of 2020, and facility operation and testing with the recompression cycle is scheduled until 2022 [16,17].

In Korea, various kinds of experimental facilities have been developed and investigated for the SCO2 power cycle [18–21]. Recently, the development and testing of the SCO2 power cycle for a waste heat recovery system started in the Korea Atomic Energy Research Institute (KAERI). The target of the project is to develop a prototype SCO2 power generation system for a waste heat recovery system. The project aims at developing SCO2 core technology, such as turbomachinery and heat exchangers. A simple recuperated cycle was selected. The major components of the prototype are an SCO2 compressor and turbine, precooler, recuperator, and waste heat recovery heat exchanger. The target of the TIT is about 430 ◦C. To secure the TIT, it is important to manage the high-temperature and high-pressure heat exchange technologies under a SCO2 environment.

In this research, the SCO2 heat exchanger that will be used in the SCO2 power generation prototype is designed. The maximum target of the SCO2 outlet temperature and pressure is 200 bar and 600 ◦C, respectively. Based on the design conditions of the SCO2 heat exchanger, heat exchanger type selection, material selection, and tube specification based on the commercial availability were conducted. Preliminary performance analysis of the SCO2 heat exchanger was conducted using an in-house heat exchanger code, and the flow and thermal stress analysis of the SCO2 heat exchanger were performed using commercially available computational fluid dynamic (CFD) codes.

#### **2. Design Considerations of Supercritical CO2 Heat Exchanger**

#### *2.1. Operating Condition of SCO2 Heat Exchanger*

Table 1 lists the operating conditions of the SCO2 heat exchanger. The mass flow rate of the SCO2 is 1 kg/s, and the inlet and outlet temperatures of the SCO2 are 300 ◦C and 600 ◦C, respectively. The outlet pressure of the SCO2 is 200 bar. Based on the operating condition of the SCO2 power cycle, the pressure drop of the heat exchanger was limited to 1.5 bar. The heat duty of the SCO2 heat exchanger was 380kW. A combustion system was selected as a heat source for the SCO2. Flue gas is composed of liquefied petroleum gas (LPG) and air. The inlet temperature of the flue gas was set at 800 ◦C. The inlet temperature is based on the maximum allowable stress for the heat exchanger tube. The inlet mass flow rate of the flue gas was calculated as 0.8497 kg/s.


**Table 1.** Operating condition of the SCO2 heat exchanger.

#### *2.2. Heat Exchanger Type Selection*

Heat exchanger types can be classified based on the number of working fluids, compactness, flow arrangements, and heat transfer mechanisms. Tubular, plate type, extended surface, and printed circuit heat exchanger types are typical heat exchangers used in industrial areas. Among the heat exchangers, the tubular heat exchanger is popular due to its flexibility: the core shape can be easily changed by the tube diameter, length, and arrangements. In addition, tubular heat exchangers are usually used in high-temperature and high-pressure conditions. A plate type heat exchanger consists of two flow membranes, and a number of plates are compressed or welded with a gasket. Therefore, it is not appropriate to use it in extreme operating conditions due to the possibility of leakage. Compared to the tubular and plate type heat exchangers, a higher effectiveness can be achieved by using extended surface heat exchangers. However, a high pressure drop can appear on extended surface heat exchangers. For compact size heat exchangers, printed circuit heat exchangers (PCHE) have been widely studied. The volume of a PCHE can be minimized up to 1/30 compared to conventional shell-and-tube heat exchangers with the same heat duty [22]. However, maintenance and inspection of a PCHE are difficult because these heat exchangers are manufactured by a diffusion bonding process. In addition, there are limitations in material selection for diffusion bonding processes.

The type of SCO2 heat exchanger can be determined by the desired operating condition. Because the target operating condition of the SCO2 heat exchanger is at a high temperature and high pressure (600 ◦C and 200 bar), the heat exchanger should endure high thermal stress and thermal shock. In addition, the maintenance and inspection of the heat exchanger should be easy. The tubular type heat exchanger has a low pressure drop and offers the least-risk design for the thermal shock resistance, and it has modest effectiveness. Because the flue gas was considered as the heat source, the pressure drop on the hot side should also be minimized. The flexible design of the tubular heat exchanger can offer a pressure drop on the flue gas side. Therefore, the tubular type was selected as the SCO2 heat exchanger.

#### *2.3. Heat Exchanger Material Selection*

Heat exchanger material selection is based on a combination of cost, moderate properties under the operating condition, fabricability, and availability. In addition, candidate materials are required to have good corrosion, oxidation, carburization, and brittleness resistance under the SCO2 condition. Because the operating condition of the SCO2 heat exchanger is at a high temperature and high pressure, it is important to determine an appropriate material that can endure extreme operating conditions. Based on the operating condition, the maximum material surface temperature was assumed to be 650 ◦C. The maximum allowable stress values were considered as criteria for the heat exchanger material selection [23]. According to the maximum allowable stress values for candidate materials in the SCO2 heat exchanger, S31042, S34709, and S34710 were selected because they have good corrosion and carburization resistance in SCO2 environments. Based on the experimental results of the corrosion and carburization resistance, experience at similar operating conditions, cost, and availability, S34709 was finally selected as the SCO2 heat exchanger material.

#### **3. Design of Supercritical CO2 Heat Exchanger**

Figure 1 illustrates a schematic of the SCO2 heat exchanger tube bundle. The tube specification was based on commercial availability and cost. The available tube diameter and the thickness of the S34709 material were 21.7 mm and 4.9 mm, respectively. The basic configuration of the SCO2 heat exchanger had a rectangular duct fed by a flue gas. The heat transfer from the hot flue gas to the cold SCO2 occurred in the rectangular duct. A staggered tube array with a counter-crossflow arrangement was considered. The tube length pitch (SL) and tube height pitch (ST) were 35 mm and 60 mm, respectively. The selection of the tube length pitch was based on the minimum thickness for pipe bends for induction and incrementing bending [24]. The straight line of the heat exchanger tube was 800 mm. The height and length of the tube were 471.7 mm and 2086.7 mm, respectively. Figure 2 shows the heat exchanger nozzle, header, and tube supporting structure. The total length of the SCO2 heat exchanger was estimated as 5132 mm, and the lengths of the SCO2 heat exchanger combustor and SCO2 heat exchanger chamber were 1577 mm and 2765 mm, respectively. As a concept for the tube-supporting structure, rectangular plates made by the welding method were considered.

Preliminary simulation of the SCO2 heat exchanger was performed using an in-house heat exchanger analysis code. Figure 3 shows a program flow chart for the analysis of the SCO2 heat exchanger. An effective number of the transfer unit method was used for the simulation code. The thermo-properties of the flue gas and SCO2 were obtained from the NIST REFPROP Database 23, Version 9.1 [25]. The fluid properties, such as Reynolds number, heat transfer coefficient, and friction factor, were calculated based on the thermal properties of the working fluids. Then, the heat transfer characteristics, such as overall heat transfer coefficient, number of transfer units, and effectiveness, were computed to obtain the outlet temperature of the working fluid. The pressure drop of each fluid was then obtained when the calculation of the outlet temperature was converged. For the shell side, the heat transfer correlation proposed by Zukauskas [26] was used. The pressure drop correlation for the staggered tube banks was employed [27]. For the tube side, the heat transfer correlation proposed by Gnielinski [28] was used, and the pressure drop was calculated by considering the entrance, momentum, core friction, and exit effects [29]. The validation of the developed code was verified with other commercially available heat exchanger code. The simulation results showed that the outlet temperature (600.6 ◦C) can be obtained with the present design considerations of the SCO2 heat exchanger. In addition, the pressure drop on the SCO2 side satisfied the design constraints (<1.5 bar). However, there is a limitation in using the in-house heat exchanger code because it only represents the outlet conditions. This means that it is difficult to find local heat transfer characteristics along the heat exchanger tubes as well as thermal stress along the tubes. Therefore, the analysis of thermal characteristics was conducted using commercial three-dimensional CFD codes.

**Figure 1.** Schematic of the SCO2 heat exchanger tube bundle.

**Figure 2.** SCO2 heat exchanger header, nozzle, and tube-supporting structure.

**Figure 3.** Heat exchanger performance analysis flow chart.

#### **4. CFD Analysis of a Supercritical CO2 Heat Exchanger**

Three-dimensional commercially available CFD codes were used to analyze the flow and thermal stress characteristics of the SCO2 heat exchanger. Flow and thermal stress analysis were performed. Based on the temperature and heat transfer coefficient distributions obtained from the flow analysis, thermal stress analysis was conducted to evaluate the tube's integrity at the SCO2 heat exchanger operating condition.

#### *4.1. Flow Analysis*

The thermal characteristics of the SCO2 heat exchanger were analyzed using a commercial CFD code, CFX. The overall performance simulation of the SCO2 heat exchanger was difficult due to the computing power. Therefore, a design of the scaled SCO2 heat exchanger was created for the CFD analysis. Figure 4 shows the scaled SCO2 heat exchanger: Figure 4a indicates the overall three-dimensional geometry, and Figure 4b illustrates the staggered tube array used in the flow CFD analysis. In the SCO2 heat exchange chamber, three heat exchanger tubes with a staggered tube array were positioned. The center of the heat exchanger tube reflected the full-scale tube specification, while the half-scale tube specification was considered in the bottom and top of the tubes. The numbers of nodes and elements for the analysis were 9,555,606 and 10,698,900, respectively. Figure 5 shows the mesh distributions on the scaled SCO2 heat exchanger. Figure 5a illustrates the overall mesh distributions in the scaled SCO2 heat exchanger, Figure 5b indicates the mesh structure around the straight tube position, and Figure 5c shows the mesh formation near the tube bending location. Precise cells near the wall surface were considered to keep the turbulence effect. Several grid layers were applied around the heat exchanger tubes to consider the wall effect on the working fluids.

For the turbulence model, a two-equation turbulence model of RNG k-ε was used: k is the turbulence kinetic energy, which is defined by the fluctuation in velocity, and ε is the turbulence eddy dissipation. The RNG k-ε model was improved from the standard k-ε model, and it was derived using the renormalization group theory. The RNG k-ε model has an additional term, which considers the eddy dissipation, average shear stress, and swirl effect. These features provide a more accurate turbulence model compared to the standard k-ε model. The transport equations for turbulence generation and dissipation are the same as those for the standard k-ε model, but the model constants differ [30].

**Figure 4.** Scaled SCO2 heat exchanger; (**a**) overall CFD analysis structure; (**b**) detailed view around heat exchanger tubes.

**Figure 5.** Mesh distributions on the scaled SCO2 heat exchanger: (**a**) overall mesh distributions, (**b**) mesh distributions around the straight tube position, (**c**) mesh formation near the tube bending area.

For the continuity equation:

$$\frac{\partial \rho}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j} (\rho \mathbf{U}\_j) = 0 \tag{1}$$

For the momentum equation:

$$\frac{\partial \rho \mathcal{U} I\_i}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j} \Big(\rho \mathcal{U} I\_i \mathcal{U}\_j\Big) = -\frac{\partial p'}{\partial \mathbf{x}\_i} + \frac{\partial}{\partial \mathbf{x}\_j} \Big(\mu\_{eff} \Big(\frac{\partial \mathcal{U}\_i}{\partial \mathbf{x}\_j} + \frac{\partial \mathcal{U}\_j}{\partial \mathbf{x}\_i}\Big)\Big) + S\_M \tag{2}$$

For the transport equation for turbulence dissipation:

$$\frac{\partial(\rho\varepsilon)}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_{j}}(\rho \mathrm{Id}\_{i}\varepsilon) = \frac{\partial}{\partial \mathbf{x}\_{j}} \bigg[ \left( \mu + \frac{\mu\_{t}}{\sigma\_{\varepsilon \mathrm{RNG}}} \right) \frac{\partial \varepsilon}{\partial \mathbf{x}\_{j}} \bigg] + \frac{\varepsilon}{k} (\mathbb{C}\_{\varepsilon \mathrm{I RNG}} P\_{\mathrm{k}} - \mathbb{C}\_{\varepsilon \mathrm{ZNG}} p \varepsilon + \mathbb{C}\_{\varepsilon \mathrm{I RNG}} P\_{\varepsilon \mathrm{b}}) \tag{3}$$

where ρ is the density, p' is the modified pressure, μ is the viscosity, SM is the sum of body forces and Pk is the shear production of turbulence, while σεRNG, Cε1RNG constant, and C<sup>ε</sup>2RNG constant are the RNG k-ε constants.

The boundary conditions were based on the actual operation conditions of the SCO2 heat exchanger. For the actual SCO2 heat exchanger, the number of heat exchanger tubes was 16, as shown in Figure 1. Uniform mass flow distributions on each heat exchanger tube were assumed: the mass flow rate of the SCO2 in each tube was 0.0625 kg/s. For the flue gas in the scaled SCO2 heat exchanger, the mass flow rate was calculated as 0.1062 kg/s. The fluid properties were implemented in the CFD code with a real gas properties table format, using NIST REFPROP Database 23, Version 9.1 [25]. For the convergence criteria, the SCO2 outlet temperature was monitored in each step as well as solution imbalances (mass, momentum, and energy) less than 1%.

Figure 6 shows the CFD results for the center position of the SCO2 heat exchanger: Figure 6a,b show the temperature and velocity distributions, respectively. High-velocity regions were focused on the left and right sides of the SCO2 heat exchange chamber because there was empty space for installing the SCO2 heat exchanger and tube displacement margin due to the thermal stress. However, the in-house code did not consider the empty space at the corner of the SCO2 heat exchange chamber. The empty space could have resulted in different simulation results between the in-house code and the CFD analysis. However, the outlet temperature of the SCO2 showed a similar performance (602 ◦C) compared to the in-house code result. The pressure drop of the SCO2 flow path was calculated as 0.374 bar, which is lower than the design constraint. Further pressure drops should be considered in the heat exchanger headers and nozzles.

**Figure 6.** *Cont.*

**Figure 6.** CFD results at the centerline of the scaled SCO2 heat exchanger: (**a**) temperature and (**b**) velocity distributions.

The temperature and heat transfer coefficient distributions on the heat exchanger are important because these parameters can influence the tube's integrity. Figure 7 shows the temperature and heat transfer coefficient distributions on the inner and outer heat exchanger tubes. Figure 7a,b are the results of the tube's outer surface and Figure 7c,d are the results of the tube's inner surface. The tube's maximum inner and outer surface temperatures were 647 ◦C and 637 ◦C, respectively. The tube's maximum surface temperature was lower than the temperature assumption value (650 ◦C) in the heat exchanger material selection. For the tube's inner area, the maximum heat transfer coefficient was about 4000 W/m2K, which was located at the tube's bending location. For the outer tube area, the maximum heat transfer coefficient was located near the outlet position of the flue gas and was 70 W/m2K. Based on the temperature and heat transfer coefficient distributions on the heat exchanger tube, thermal stress analysis was performed.

**Figure 7.** *Cont.*

**Figure 7.** CFD results around the heat exchanger tube: (**a**) temperature distributions on the outer tube surface, (**b**) heat transfer coefficient distributions on the outer tube surface, (**c**) temperature distributions on the inner tube surface, (**d**) heat transfer coefficient distributions on the inner tube surface.

#### *4.2. Thermal Stress Analysis*

To show the integrity of the SCO2 heat exchanger tube, thermal stress analysis was conducted using a commercial CFD code, ABAQUS [31]. High-temperature regions of the heat exchanger tube located near the SCO2 outlet region (flue inlet region) were analyzed. Therefore, 10 straight lines with 10 bending flow paths were modeled for the thermal stress simulation. With consideration of the symmetrical structure of the heat exchanger tube, a half scale of the heat exchanger tube was analyzed. 1,114,135 elements were used for the thermo-mechanical analysis. Three-dimensional continuum element DCC3D8 was used to obtain the temperature distribution, and the stress analysis was carried out using the C3D8 continuum solid element.

In order to obtain temperature distributions for the heat exchanger tube, temperature and heat transfer coefficient profiles obtained from the thermal analysis were used as the input for thermal stress analysis. Because the pressure drop in the heat exchanger and heat exchange chamber was negligible compared to the operating condition, the pressures in the SCO2 heat exchanger tube and the chamber were assumed to be 200 bar and atmospheric pressure, respectively. Figure 8 illustrates the temperature distribution in the scaled heat exchanger tube. A solid temperature of the heat exchanger tube is closer to the SCO2 temperature because the SCO2 heat transfer coefficient is higher than that of the flue gas. The evaluation of thermal stress analysis was based on the allowable stress value of S34709. The allowable stress value at a temperature of 650 ◦C is 539 bar, which is based on ASME Sec. II [23]. The finite element stress analysis was performed with three cases separately: (1) thermal loading, (2) pressure loading, and (3) thermal and pressure loading.

**Figure 8.** Heat exchanger tube temperature distributions used in the thermal stress analysis.

Figure 9 shows the tube displacement due to the temperature distribution in the heat exchanger tube. Figure 9a,b represent the tube displacement shape and the principal strain distributions, respectively. The high-temperature area near the exit of the SCO2 had a large tube expansion, and it gradually decreased. The principal strain distributions showed similar distributions compared to the tube displacement. The maximum principal strain value was discovered to be 0.017. If the deformation of the heat exchanger tube is not considered, thermal buckling due to excessive pressure stress can occur. In the present SCO2 heat exchanger design, empty space in the chamber was considered, which maintained the integrity of the tube even when the tube displacement appeared.

**Figure 9.** Tube displacements on the tube: (**a**) tube displacement shape, (**b**) the maximum principal strain distributions.

Von-Mises stress distributions are shown in Figure 10 along the scaled heat exchanger tube. Figure 10a,b are the results of the thermal and pressure loading cases, respectively. In the case of thermal loading, the maximum stress was found near the tube bending area. On the other hand, the maximum stress was located at the tube bending area for the pressure loading condition. The maximum local von-Mises stress values were 35.5 bar and 516 bar for the thermal and pressure loading cases, respectively. It was confirmed that these stress values were lower than the allowable stress value. Figure 10c is the result of thermal stress analysis considering both the thermal and pressure loading cases. The stress distributions in the tube were similar to the combination of the thermal and pressure loading cases. The maximum von-Mises value was calculated as 523 bar, which is lower than the allowable stress value considered in the present study. The stress evaluation loaded on the heat exchanger tube was conducted based on the ASME Sec. VIII. The membrane stress of the heat exchanger tube was 287 bar, while the allowable stress value was 539 bar. The sum of the membrane stress and the bending stress was 378 bar, which is lower than the constraint value (1.5 × allowable stress value = 808 bar). Therefore, the stress state of the SCO2 heat exchanger satisfied the ASME criteria.

**Figure 10.** Von-Mises stress distributions: (**a**) thermal loading, (**b**) pressure loading, (**c**) thermal and pressure loading cases.

#### **5. Conclusions**

This study focused on the design of an SCO2 heat exchanger for obtaining high-temperature and high-pressure heat exchange technologies under an SCO2 environment. A tubular type heat exchanger was selected because it has high durability in extreme conditions, such as having low pressure losses in both the hot and cold sides. The heat exchanger material selection was conducted based on the maximum allowable stress, corrosion resistance, cost, and availability. A staggered tube array with a counter-cross flow arrangement was determined and the overall size of the SCO2 heat exchanger was based on the tube bending criteria and the results of in-house heat exchanger performance code. Commercially available three-dimensional CFD codes were then used to analyze the flow and thermal characteristics of the SCO2 heat exchanger. The temperature and heat transfer coefficient distributions on the SCO2 heat exchanger were analyzed. Then, thermal stress analysis was conducted based on the obtained flow analysis results. The stressed state of the SCO2 heat exchanger was evaluated based on the ASME procedure. The membrane stress, bending stress, and local stress were lower than the allowable stress. The results indicate that the stress of the present heat exchanger satisfied the ASME criteria. Based on the design of the SCO2 heat exchanger, the manufacturing process can be performed.

**Author Contributions:** Conceptualization, H.S. and J.E.C.; Methodology, H.S. and J.E.C.; Software, H.S., J.K., and Y.-W.K.; Validation, H.S., J.K., and Y.-W.K.; Formal analysis, H.S., J.E.C., and I.S.; Investigation, H.S., J.K., and I.S.; Resources, H.S., Data Curation, H.S., J.K., and Y.-W.K.; Writing-Original Draft, H.S.; Writing-Review & Editing, J.E.C.; Visualization, H.S.; Supervision, H.S.; Project administration, H.S.; Funding acquisition, J.E.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was financially supported by the institute of Civil Military Technology Cooperation funded by the Defense Acquisition Program Administration and Ministry of Trade, Industry and Energy of Korean government under grant No. 17-CM-EN-04.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**

#### **Nomenclature**


#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Performance Analysis of the Supercritical Carbon Dioxide Re-compression Brayton Cycle**

#### **Mohammad Saad Salim, Muhammad Saeed** † **and Man-Hoe Kim \***

School of Mechanical Engineering & IEDT, Kyungpook National University, Daegu 41566, Korea; msaadsalim@hotmail.com (M.S.S.); saeed.aarib@gmail.com (M.S.)

**\*** Correspondence: kimmh@asme.org; Tel.: +82-53-950-5576

† Present Address: Mechanical Engineering Department, Khalifa University of Science and Technology, SAN Campus, P.O. Box 2533, Abu Dhabi, UAE.

Received: 30 December 2019; Accepted: 5 February 2020; Published: 7 February 2020

**Abstract:** This paper presents performance analysis results on supercritical carbon dioxide (*sCO*2) re-compression Brayton cycle. Monthly exergy destruction analysis was conducted to find the effects of different ambient and water temperatures on the performance of the system. The results reveal that the gas cooler is the major source of exergy destruction in the system. The total exergy destruction has the lowest value of 390.1 kW when the compressor inlet temperature is near the critical point (at 35 ◦C) and the compressor outlet pressure is comparatively low (24 MPa). The optimum mass fraction (*x*) and efficiency of the cycle increase with turbine inlet temperature. The highest efficiency of 49% is obtained at the mass fraction of *x* = 0.74 and turbine inlet temperature of 700 ◦C. For predicting the cost of the system, the total heat transfer area coefficient (*UATotal*) and size parameter (*SP*) are used. The *UATotal* value has the maximum for the split mass fraction of 0.74 corresponding to the maximum value of thermal efficiency. The *SP* value for the turbine is 0.212 dm at the turbine inlet temperature of 700 ◦C and it increases with increasing turbine inlet temperature. However the *SP* values of the main compressor and re-compressor increase with increasing compressor inlet temperature.

**Keywords:** re-compression Brayton cycle; carbon dioxide; supercritical; thermodynamic; exergy; cycle simulation; design point analysis

#### **1. Introduction**

The main cause of pollution is the combustion of fossil fuels to create energy for heavy industrialization and urbanization. Fossil fuel reserves are diminishing due to this process; thus, a big demand for power generation from green energy sources at high efficiency has been created. Global warming is another big concern, as has been pointed out by the United Nations Framework Convention on Climate Change (UNFCC). The proposal from the conference was to undertake efforts so that the rise in global average temperature increase could be limited to well below 2 ◦C above preindustrial levels. Due to the use of fossil fuels and harmful working fluids, significantly harmful effects on the environment are causing problems such as global warming and acid rain. The effects of pollution tend to bring unpredictable changes in the global climate, as has been asserted by the Intergovernmental Panel on Climate Change (IPCC), and rising sea levels are making large parts of the Earth uninhabitable [1]. Thus, green sources of energy are the need of the hour to solve these issues, which has led to research being conducted on different forms of green energy, such as biogas [2,3], geothermal energy [4,5], energy from human excreta [6] and solar energy [7]. In order to cope with the aforementioned global climate challenges, carbon dioxide (*CO*2)-based power systems present an environmentally friendly option and are capable of providing power at high efficiency.

The Rankine cycle and the air-standard Brayton cycle are well-known thermodynamic cycles. The benefit of the Rankine cycle is that high efficiency can be achieved because the pump consumes a very small amount of work, since compression is carried out when the working fluid is in the liquid state [8]. The advantage of the Brayton cycle is that the turbine inlet temperature is high; thus, it can achieve high efficiency, but the disadvantage is that the work consumed by the compressor is very large. Due to this, the air-standard Brayton cycle's efficiency is not significantly higher than that of the steam-based Rankine cycle. The primary advantage of the supercritical *CO*<sup>2</sup> (*sCO*2) Brayton cycle is that the positive points of the steam-based Rankine cycle and air-standard Brayton cycle are both combined. The turbine inlet temperature in the *sCO*<sup>2</sup> Brayton cycle is high. Moreover, since the compressor operates near the *CO*<sup>2</sup> critical point at very high pressure at which the density is significantly high and the compressibility factor is small, the work that is consumed by the compressor is significantly low. The *sCO*<sup>2</sup> Brayton cycle operates above the critical point, so the need for condensing the system is removed and the system has a simple layout. Since the *sCO*<sup>2</sup> re-compression Brayton cycle has very high operating pressures compared to the steam-based Rankine cycle, the size of the *sCO*2-based power system's components is considerably smaller [9].

Currently, extensive research is being conducted on *sCO*2-based power systems, and the *sCO*<sup>2</sup> Brayton cycle can be found in a variety of arrangements in the literature, such as in reheated and intercooled re-compression layouts [10]. Crespi et al. [11] reviewed the different single and combined layouts of *sCO*<sup>2</sup> Brayton cycle power systems with efficiencies of 40%–50% and 50%–60%, respectively, while Saeed and Kim [8] analyzed a re-compression *sCO*<sup>2</sup> Brayton cycle power system with an integrated turbine design and optimization algorithm. In their study, they proposed that the cycle performs best when the inlet temperature of the compressor is set near the *CO*<sup>2</sup> critical temperature (i.e., 32–37 ◦C) and the compressor inlet pressure is set slightly above the *CO*<sup>2</sup> critical pressure (i.e., 7.8–8.1 MPa) along with a moderate pressure ratio (i.e., 2.9–3.1). Saeed et al. [12] carried out a design optimization and performance analysis of the *sCO*<sup>2</sup> re-compression Brayton cycle. They developed detailed mathematical models of the cycle components and simulation codes for the turbine, compressor and heat exchanger. These codes were used to analyze the performance of the cycle under the design conditions as well as off-design conditions.

To analyze the performances of thermal power systems, researchers have used different performance parameters, including thermal efficiency (η*th*) and exergy efficiency (η*ex*). Moreover, to indicate the cost and size of the thermodynamic system, parameters such as the total heat transfer area coefficient (*UATotal*) and size parameter (*SPTotal*) have also been used to indicate the heat exchanger and turbomachinery sizes, respectively [13–15]. Patel et al. [16] studied the optimization of a waste-heat-based organic Rankine cycle (ORC)-powered cascaded vapor compression-absorption refrigeration system. They used the log mean temperature difference (LMTD) to determine the *UA* value of each heat exchanger in the system. The purpose was to minimize the *UA* value in order to minimize the area required for heat exchange, and thus minimize the cost of the heat exchangers in the system.

This paper presents the various benefits of *sCO*<sup>2</sup> power systems. A detailed investigation has been conducted for the performance of the system with regard to key performance parameters such as η*th*, η*ex*, *UA* and size parameter (*SP*) of the turbomachinery. The system analysis also has been performed on the basis of the changing ambient temperature (*T*0) and water temperature (*Tw*,*in*,) values to signify how the system's exergetic performance changed on a monthly basis. To the author's best knowledge, monthly/seasonal analysis using for *sCO*2-based power systems are not available in the literature. Moreover, the performance of the system at different turbine inlet temperatures (*T*7), compressor inlet temperatures (*T*1) and compressor outlet pressures (*P*2) is also presented to signify the effect of each variable on various performance parameters of the system, and to indicate which values are best for use as the cycle's design points.

#### **2. Methodology**

#### *2.1. Cycle Processes*

Figure 1 shows a cycle layout and temperature-entropy (T-s) diagram of the system considered in the study. As shown in Figure 1a, the system consists of a turbine, a primary heat exchanger, high and low-temperature recuperators, the main compressor, a re-compression compressor and a gas cooler. In comparison with the recuperated cycle, the re-compression Brayton cycle includes an added intermediate compressor, an additional recuperator and split/mixing flow values. A fraction of mass is taken from the mainstream and fed to the re-compression compressor bypassing the low-temperature recuperator (LTR). This fraction of flow enters the mainstream again increases its temperature before entering into the high-temperature recuperator (HTR). This arrangement in turns increases the thermal efficiency of the cycle (42% to 50% for cycles operating with the lowest and highest temperatures, 37 ◦C and 700 ◦C respectively [10]) by reducing heat rejection in the pre-cooler.

**Figure 1.** Schematic of re-compression *sCO*<sup>2</sup> Brayton cycle.

The stream exiting the re-compressor in State 5 and the stream exiting the recuperator in State 3 are combined in State 4. This stream then passes the HTR where its temperature further increases to State 6. After this, the heat addition process takes place at the main heat exchanger and the stream reaches *T7* in State 7. In the turbine, the expansion process takes place until State 8. From State 8 to State 9, the HTR recuperates heat, and then, from State 9 to State 10, the LTR recuperates heat. After this, cooling takes place in the gas-cooler in State Point 1. The gas cooling process from State 10 to State 1 is used to transfer heat to the coolant (water) which is input to the gas cooler at the monthly temperature value (see Figure 2) and has to be heated to 40 ◦C for domestic uses, such as floor heating.

#### *2.2. Design Parameters*

The cycle simulations were developed using MATLAB [17], and the thermodynamic properties of *CO*<sup>2</sup> at each state point in the cycle were calculated NIST's REFPROP [18]. The design parameters and governing equations for defining the operation of the cycle are listed in Table 1.

**Figure 2.** Ambient and water temperature data for the city of Busan [19].



For the monthly exergy analysis presented in this study, ambient temperature (*T0*) and water (coolant) inlet temperatures (*Tw,in*) in the gas cooler were used from the data for the city of Busan [19] as shown in Figure 2.

#### *2.3. Energy and Exergy Analysis*

To simplify the analysis, the following assumptions are made:


The following set of governing equations represent the energy analysis for the *sCO*<sup>2</sup> cycle:

$$\mathcal{W}\_T = m\_{CO\_2} (h\_7 - h\_8)\_\prime \tag{1}$$

$$\mathcal{W}\_{\rm MC} = m\_{\rm CC} \mathbf{x} (h\_2 - h\_1),\tag{2}$$

$$\mathcal{W}\_{\rm RC} = m\_{\rm CO\_2} (1 - \mathbf{x}) (h\_5 - h\_{10}),\tag{3}$$

$$Q\_{\rm in} = m\_{\rm CO\_2} (h\_7 - h\_6),\tag{4}$$

$$Q\_{\mathbb{GC}} = m\_{\mathbb{GC}2} \mathbf{x} (h\_{10} - h\_1),\tag{5}$$

$$\mathcal{W}\_{net} = \mathcal{W}\_T - \mathcal{W}\_{\text{MC}} - \mathcal{W}\_{\text{RC}\_{\prime}} \tag{6}$$

$$
\eta\_{\rm tl} = \frac{\mathcal{W}\_{\rm net}}{Q\_{\rm in}}.\tag{7}
$$

Thermal efficiency (η*th*) is a parameter that takes into account how much of the energy input from the heat source is converted into the net shaft work by the turbine, but it does not reflect the irreversibilities that are involved in the *sCO*<sup>2</sup> re-compression Brayton cycle. The exergy analysis of the system is significantly useful, since energy is conserved but exergy is destroyed by the irreversibility [20]. The second law (exergy) efficiency of the system is an indicator of the maximum theoretical work that can be generated as the system is brought to equilibrium with the environment.

The specific exergy of flow at any state in the system is given by

$$
\sigma = h - h\_0 - T\_0(s - s\_0). \tag{8}
$$

By applying exergy balance over each component, the exergy destruction in the different components of the system can be defined as follows.

**Primary heat exchanger:**

$$I\_{HX} = m\_{\text{CO}\_2}(\mathbf{e}\_6 - \mathbf{e}\_7) + m\_s(\mathbf{e}\_{\text{si}} - \mathbf{e}\_{\text{so}}),\tag{9}$$

where *ms* is the mass flow rate of the heat transfer fluid.

**Main compressor:**

$$I\_{\rm MC} = m\_{\rm CO\_2} \mathbf{x}(e\_1 - e\_2) + W\_{\rm MC}.\tag{10}$$

**Re-compressor:**

$$I\_{R\mathbb{C}} = m\_{CO\_2} (1 - \mathfrak{x}) (\mathfrak{e}\_{10} - \mathfrak{e}\_5) + \mathcal{W}\_{RC}.\tag{11}$$

**Gas cooler:**

$$I\_{\rm GC} = m\_{\rm CO\_2} \mathbf{x} (\mathbf{e}\_{10} - \mathbf{e}\_1) + m\_{\rm av} (\mathbf{e}\_{\rm air} - \mathbf{e}\_{\rm uvout}) . \tag{12}$$

**Turbine:**

$$I\_T = m\_{\rm CO\_2} (\mathbf{e}\_7 - \mathbf{e}\_8) - W\_T. \tag{13}$$

**High-temperature recuperator:**

$$I\_{HTR} = m\_{CO\_2} (\mathbf{e}\_8 - \mathbf{e}\_9 + \mathbf{e}\_4 - \mathbf{e}\_6). \tag{14}$$

**Low-temperature recuperator:**

$$I\_{\rm LTR} = m\_{\rm CO\_2} [\varepsilon \mathfrak{g} - \varepsilon\_{10} + \mathfrak{x}(\varepsilon\_2 - \varepsilon\_3)].\tag{15}$$

The total exergy destruction in the system is given by

$$I\_{Total} = I\_{HX} + I\_{M\underline{C}} + I\_{R\underline{C}} + I\_{G\underline{C}} + I\_T + I\_{HTR} + I\_{LTR} \tag{16}$$

The heat rejection from the system in the gas cooling process, which can be utilized for district heating purposes, is determined from the following equation:

$$Q\_{out} = m\_w(h\_{nout} - h\_{win}).\tag{17}$$

By balancing the exergy throughout the whole system, the exergy input to the system is

$$E\_{\rm in} = W\_{\rm net} + I\_{Total} + \left(\frac{T\_0}{T\_{w\,\,avg}} - 1\right) Q\_{\rm out}.\tag{18}$$

The exergy efficiency of the system is given by

$$
\eta\_{\rm ex} = \frac{\mathcal{W}\_{\rm net}}{E\_{\rm in}}.\tag{19}
$$

The total irreversibility ratio (*IR*) of the system can be written as

$$IR = \frac{I\_{\text{total}}}{E\_{in}}.\tag{20}$$

#### *2.4. Total Heat Transfer Area Coe*ffi*cient and Total Size Parameter Value (SPTotal)*

*UATotal* and *SPTotal* are useful parameters for estimating the cost of the system. *UATotal* indicates the total area required for the heat exchangers in the *sCO*<sup>2</sup> re-compression Brayton cycle, and thus the cost associated with the investment and maintenance of the heat exchangers. The hypothesis is that differences between the heat transfer coefficients in the *sCO*<sup>2</sup> cycle heat exchangers are not significant. The total heat transfer area in the heat exchangers increases as a result of increasing *UATotal*, thereby increasing the investment and maintenance costs of the heat exchangers, and thus the system [21]. A smaller value of *UATotal* is desirable for a cost-effective design.

For a heat exchanger in general, we have

$$
\Delta IL = \frac{\mathcal{Q}}{\Delta T\_M} \,\tag{21}
$$

$$
\Delta T\_M = \frac{\Delta T\_{\text{max}} - \Delta T\_{\text{min}}}{\ln \left( \frac{\Delta T\_{\text{max}}}{\Delta T\_{\text{min}}} \right)} \,\tag{22}
$$

where Δ*TM* is logarithmic mean temperature difference and Δ*Tmax* and Δ*Tmin* are the maximum and minimum temperature differences, respectively, at the two ends of the heat exchanger. However, in the case of supercritical carbon dioxide Brayton cycle (*sCO*2-*BC*), properties of the working fluid change swiftly, and the definition of the Δ*TM* is not applicable in this case. In order to cope with the situation, the length of the heat exchanger was divided into N number of segments. The length of each segment was kept small enough that the variation within a particular small segment can be ignored, as shown in Figure 3. The value of the number of segments (N) depends upon the operation region of the heat exchanger. For the pre-cooler and LTR, the numbers of segments will be large, as these two heat exchangers operate close to the critical point, and properties of *sCO*<sup>2</sup> change at high rates. For HTR, in contract with LTR and pre-cooler, but operating away from the critical point where variation in the properties is small, a lesser number of segments would be required. Further details on the model can be found in the previous studies [8,12,22,23].

Based on the discretized model definition of the logarithmic mean temperature difference (Δ*TM*) *i* , the *i th* segment is given by the Equation (23)

$$\begin{aligned} \left(\Theta\_{2}\right)^{\mathrm{i}} &= \left(\mathrm{T}^{\mathrm{h}}\right)^{\mathrm{i}} - \left(\mathrm{T}^{\mathrm{c}}\right)^{\mathrm{i}};\\ \left(\Theta\_{1}\right)^{\mathrm{i}} &= \left(\mathrm{T}^{\mathrm{h}}\right)^{\mathrm{i}+1} - \left(\mathrm{T}^{\mathrm{c}}\right)^{\mathrm{i}+1};\\ \left(\Delta\mathrm{T}\_{\mathrm{M}}\right)^{\mathrm{i}} &= \frac{\max\{\left(\left(\bullet\_{1}\right)^{\mathrm{i}}, \left(\bullet\_{2}\right)^{\mathrm{i}}\right) - \min\{\left(\bullet\_{1}\right)^{\mathrm{i}}, \left(\bullet\_{2}\right)^{\mathrm{i}}\}}{\log\left(\frac{\max\left(\left(\bullet\_{1}\right)^{\mathrm{i}}, \left(\bullet\_{2}\right)^{\mathrm{i}}\right)}{\min\left(\left(\bullet\_{1}\right)^{\mathrm{i}}, \left(\bullet\_{2}\right)^{\mathrm{i}}\right)}\right)} \end{aligned} \tag{23}$$

Thus, UA value through the *i th* segment could be computed as given by the following equation.

$$(\ellIA)^i = \frac{(dq)i}{\left(\Delta T\_M\right)^i}.\tag{24}$$

And UA value associated with the whole heat exchanger can be calculated using the following relation.

$$
\natural LA = \sum\_{i=1}^{N} (\natural LA)^{i}.\tag{25}
$$

Later, combining the *UA* values for the *sCO*<sup>2</sup> system, we obtain

$$
\Pi A\_{Total} = \Pi A\_{HX} + \Pi A\_{HTR} + \Pi A\_{LTR} + \Pi A\_{GC}.\tag{26}
$$

*SP* indicates the size of the turbomachinery in the system, and a high *SPTotal* means comparatively large sizes of the turbine and compressors. Consequently, the investment and maintenance costs of the turbomachinery will increase, which, in turn, will increase the total cost of the system. Because of this, a small *SPTotal* is desirable.

**Figure 3.** Discretized model of the heat exchanger to capture the effect of variation in the properties of sCO2.

Macchi and Perdichizzi [14] used the turbine *SP* to determine the expander size, which for the *sCO*<sup>2</sup> re-compression Brayton cycle, is given by

$$SP\_T = \frac{\sqrt{V\_{8s}}}{\sqrt[4]{\Delta h\_s}} = \frac{\sqrt{m\_{CO\_2}v\_{8s}}}{\sqrt[4]{h\_7 - h\_{8s}}},\tag{27}$$

where *V8s* is the isentropic value of the volume flow rate of the working fluid at the turbine exit, and Δ*hs* is the isentropic specific enthalpy drop.

The main compressor and re-compression compressor *SP* values can be respectively defined as

$$SP\_{\rm MC} = \frac{\sqrt{m\_{\rm CO\_2} \pi v\_1}}{\sqrt[4]{h\_{2s} - h\_1}},\tag{28}$$

$$SP\_{RC} = \frac{\sqrt{m\_{CO\_2}(1-\chi)v\_{10}}}{\sqrt[4]{h\_{5s} - h\_{10}}}.\tag{29}$$

For the system under study,

$$SP\_{Total} = SP\_T + SP\_{MC} + SP\_{RC} \,. \tag{30}$$

Hence,

$$(SP)\_{Total} = \frac{\sqrt{m\_{CO\_2}v\_{8s}}}{\sqrt[4]{h\_7 - h\_{8s}}} + \frac{\sqrt{m\_{CO\_2}xv\_1}}{\sqrt[4]{h\_{2s} - h\_1}} + \frac{\sqrt{m\_{CO\_2}(1-x)v\_{10}}}{\sqrt[4]{h\_{5s} - h\_{10}}}.\tag{31}$$

#### *2.5. Properties of sCO*<sup>2</sup>

The selection of the working fluid has an important role in optimizing the performance of the system. Efficiency, cost and environmental aspects are key performance metrics that need to be considered before selecting a working fluid. Amongst the desired characteristics are an ozone depletion potential (ODP) of 0 and a low global warming potential (GWP) value. Table 2 lists the properties of *CO*2, which is the working fluid used for the system in the study. Moreover, *CO*<sup>2</sup> is a non-toxic, cheap and easily available substance. A point of interest is that *CO*<sup>2</sup> has a critical temperature that is very close to the ambient temperature; thus, it is relatively easy for it to be used in the supercritical state, and it can be matched to a number of power cycles. The properties of *CO*<sup>2</sup> show large variations near their critical points, which is highly advantageous for the compression process. With a small increase in temperature, the energy of the fluid increases significantly during the compression process, thereby making it highly efficient [10].

**Table 2.** Typical properties of *CO*2.


#### **3. Results and Discussion**

#### *3.1. Thermal E*ffi*ciency (*η*th)*

Figure 4 shows thermal efficiency, η*th* of the *sCO*<sup>2</sup> system. As expected, the increment of *T7* increases the efficiency of the system, meaning that more power output can be derived for the same heat input. The optimal value for the split mass fraction (*x*) is the one for which the thermal efficiency shows the maximum value for each *T*<sup>7</sup> case. Moreover, as *T*<sup>7</sup> increases, the optimal value of *x* slightly shifts to the right. The best efficiency is at *x* = 0.74 and *T*<sup>7</sup> = 700 ◦C.

**Figure 4.** The thermal efficiency of the *sCO*<sup>2</sup> cycle at different turbine inlet temperatures (*T*7).

#### *3.2. Exergy E*ffi*ciency (*η*ex) and Irreversibility Ratio (IR)*

Figure 5 presents the *IR* and the η*ex* values for the *T*<sup>7</sup> of 700 ◦C. The optimal value for *x* is the one for which η*ex* shows the maximum value. The *IR* is the minimum, and at the same time, η*ex* is the maximum at *x* = 0.74.

**Figure 5.** Exergy efficiency (η*ex*) and irreversibility ratio (*IR*) of the *sCO*<sup>2</sup> cycle at a turbine inlet temperature (*T*7) of 700 ◦C.

#### *3.3. Exergy Destruction in Di*ff*erent Components*

The total exergy destruction of the *sCO*<sup>2</sup> re-compression Brayton cycle increased with an increase in *T*7, as shown in Figure 6. One of the reasons for this is the higher entropy generation at higher *T*<sup>7</sup> values causing comparatively higher exergy destruction in the turbine. Figure 7 presents the exergy destruction in different components of the *sCO*<sup>2</sup> cycle at *T*<sup>7</sup> = 700 ◦C and *x* = 0.74. It can be observed that the gas cooler shows the maximum exergy destruction because a lot of exergy is destroyed when the *sCO*<sup>2</sup> is cooled down from State 10 to State 1. Figure 8 shows that when the outlet pressure of the compressor is high, the exergy destruction is also comparatively high. Moreover, in the summer months, the exergy destruction was high because of high ambient and water temperatures. Figure 9 shows that as the compressor inlet temperature is kept to a lower value (near the critical point), the exergy destruction is minimized. The minimum exergy destruction of 390 kW was observed in January at a compressor inlet temperature of 35 ◦C and a compressor outlet pressure of 24 MPa.

**Figure 6.** Total exergy destruction vs. turbine inlet temperature (*T*7).

**Figure 7.** Exergy destruction in different components of the cycle at *T*<sup>7</sup> = 700 ◦C.

**Figure 8.** Total exergy destruction at different compressor outlet pressures (P2).

**Figure 9.** Total exergy destruction at different compressor inlet temperatures (*T*1).

#### *3.4. Heat Transfer Area Coe*ffi*cient (UA) and Turbomachinery Size (SP) Values*

*UA* and *SP* values are the parameters for the estimation of the cost of the system. *UA* value indicates the area required for heat exchangers, and thus indicates their cost. As the *UA* value increases, the total heat transfer area in the heat exchangers also increases, thereby increasing the investment and maintenance costs of the heat exchangers, and thus the system [21]. *SP* values indicate the size of the turbomachinery in the system; thus, a lower value is desirable.

Figure 10 illustrates the relationship of *UA* with the different heat exchangers in the *sCO*<sup>2</sup> re-compression Brayton cycle with different *T*<sup>7</sup> values. It should be noted that the *UA* value for the primary heat exchanger was the highest and increased at an increasing rate, whereas the *UA* value for the HTR showed a linear increase with respect to *T*7. The net result is that for a high *T*<sup>7</sup> value, a large heat exchange area is required when considering all of the heat exchangers in the system. The gas cooler showed the lowest *UA* value and was not sensitive to *T*7.

**Figure 10.** Heat transfer area coefficient (*UA*) vs. turbine inlet temperature (*T*7).

Figure 11 shows the effect of *T*<sup>1</sup> on the heat exchanger sizes, as indicated by the *UA* value for each heat exchanger. It should be noted that the *UA* values of the primary heat exchanger, gas cooler and LTR all decrease with increasing *T*1, whereas that of HTR shows only a slight increase. As shown in Figure 12, the net result was a decrease in the *UATotal* value with increasing *T*1, which implies that the total area required for heat exchange decreases, and smaller heat exchanger sizes will be required for higher compressor inlet temperatures (which are much higher than the critical temperature value).

**Figure 11.** Heat transfer area coefficient (*UA*) vs. compressor inlet temperature (*T*1).

**Figure 12.** Total heat transfer area coefficient (*UATotal*) vs. compressor inlet temperature (*T*1).

Figure 13 shows the relationship between *UATotal* and *x* in the *sCO*<sup>2</sup> re-compression Brayton cycle. The plot shows a maximum *UATotal* value at *x* = 0.74, which is the same point at which η*th* was maximum. This implies that the size, and thus, the cost of the heat exchangers in the system, will need to be high to achieve the maximum efficiency.

**Figure 13.** Total heat transfer area coefficient (*UATotal*) vs. split mass fraction (*x*) at *T*<sup>7</sup> = 700 ◦C.

Figures 14 and 15 show the *SP* values for the turbine, main compressor and re-compressor in the *sCO*<sup>2</sup> cycle. The *SP* of the turbine increases at higher *T*<sup>7</sup> values, implying that a large turbine should be used for high *T*<sup>7</sup> values. Moreover, the *SP* values of the main compressor and re-compressor are much less sensitive to *T*7, meaning that the same compressor size can be used when *T*<sup>7</sup> is either low or high. The main compressor has a higher *SP* value than the re-compressor. Figure 15 shows the variation in *SP* of the different turbomachinery components with *T*1. The *SP* of the turbine does not depend on *T*1, as can be seen by the straight line. However, the *SP* values of the main compressor and re-compressor increase with an increase in *T*1. Moreover, the increase in the *SP* value for the main compressor is more pronounced than that of the re-compressor. Hence, relatively large compressor sizes will be required at a high *T*1.

**Figure 14.** Size parameter (*SP*) vs. turbine inlet temperature (*T*7).

**Figure 15.** Size parameter (*SP*) vs. compressor inlet temperature (*T*1).

#### **4. Conclusions**

The present study investigated seasonal performance analysis of s*CO*2-based re-compression Brayton power system. From the analysis, the following conclusions were drawn:


#### **5. Future Work**


**Author Contributions:** M.S.S. did the simulation analysis and drafted the manuscript. M.S. did heat exchanger analysis and edited the manuscript. M.-H.K. supervised the research and edited the manuscript. The authors read and approved the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**


#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
