Research on Dynamic Modeling of the Supercritical Carbon Dioxide Power Cycle

: The supercritical carbon dioxide (SCO 2 ) Brayton cycle, as a substitute for the steam cycle, can be widely used in a variety of power generation scenarios. However, most of the existing SCO 2 cycle studies are restricted to basic thermodynamics research, parameter optimizations, system design in different application ﬁelds, and even economic analysis. Considering the load variability and control ﬂexibility of the power generation system, the dynamic performance research of the SCO 2 cycle is also crucial, but the work done is still limited. Based on the previous studies, Simulink software is used in this paper to develop a dynamic model of the 20 MW-SCO 2 recompression cycle, which speciﬁcally includes component models that can independently realize physical functions and an overall closed-loop cycle model. A series of comparative calculation are carried out to verify the models and the results are very positive. The SCO 2 recompression power system is built with the developed models and the dynamic model runs stably with a maximum error of 0.56%. Finally, the simulation of the dynamic switching conditions of the 20 MW-SCO 2 recompression cycle are performed and the analysis results supply instructive suggestions for the system operation and control.


Introduction
The supercritical carbon dioxide (SCO 2 ) cycle can be traced back to the 1960s, after which research on the cycle was temporarily shelved due to technical limitations [1]. As the fourth-generation nuclear reactor developed rapidly in the 21st century, the operating temperature now reaches as high as 500-900 • C. Considering that the material reliability of the ultra-supercritical steam cycle cannot be guaranteed [2], the SCO 2 Brayton cycle as an alternative cycle has attracted widespread attention from scholars. Additionally, comparative studies on the SCO 2 cycle, steam cycle, and helium cycle have also been carried out [3,4]. A large number of results show that the SCO 2 cycle is more suitable for higher temperature nuclear reactors due to the simple structure and high efficiency, and has the potential to reduce costs.
The SCO 2 cycle also has the advantages of a compact structure and good heat source matching performance, and can be widely used in nuclear energy [5][6][7], waste heat recovery [8][9][10], solar energy [11][12][13], and thermal power [14][15][16]. A lot of work has been done on the basic thermodynamic research of the SCO 2 cycle, configuration and parameter optimizations, system design in different application fields, and even economic analysis and optimization.
Sarkar [17] performed second-law analysis of the SCO 2 cycle, and the results showed that the minimum operating temperature has a greater impact on the optimal pressure ratio and cycle efficiency than the maximum operating temperature, and the irreversibility of heat exchangers is significantly higher than that of turbomachinery. AHN et al. [18] found that the recompression cycle showed the highest efficiency by comparing different cycle layouts. Moreover, the recompression cycle can effectively avoid the recuperator pinch point, which is more likely to occur in the simple recuperation cycle [19]. Ma et al. [20] analyzed the SCO 2 cycle with an intercooling system of the main compressor, and studied the influence of the pressure ratio and the pressure distribution of the two compression stages on the performance of the system. The results show that the coupled main compressor intercooling system can improve the cycle efficiency, save system cost, and improve cycle stability. Liu et al. [14] studied the effects of the reheat stages on the cycle efficiency. Li et al. [21] conceptually designed the power system of the SCO 2 cycle coupled with a small lead-cooled fast reactor, and analyzed the thermodynamic and economic performance. In addition, some scholars have conceptually designed the configurations of SCO 2 coal-fired power plants with different loads, and proposed a series of SCO 2 boiler and cycle-side temperature coupling methods to rationally utilize the flue gas heat [15,22,23]. Zhu et al. [24] used the weighted quality method to evaluate the cost of an SCO 2 boiler, and then evaluated the cost of SCO 2 coal-fired power generation, which shows the economic advantages of the SCO 2 cycle in the thermal power field. At the same time, a large quantity of related research has also been carried out in the fields of solar energy and waste heat recovery [25,26].
The above research has mainly focused on the steady state calculation. Since the SCO 2 cycle is used in the power generation system, attention needs to be focused on the performance of a variable load and flexibility. Thus, research on the performance and simulation of the SCO 2 cycle has also attracted attention from all over the world [27]. Dyreby [28] studied the method of SCO 2 Brayton cycles numerical modelling. Moisseytsev and Sienicki [29] simulated the performance characteristics of the 100 kW simple recuperated Brayton cycle configuration by developing The Plant Dynamics Code (PDC). Wu P. et al. [30] and Deng T. et al. [31] also developed the transient analysis code for the SCO 2 Brayton cycle.
As a supplement to the previous research, this paper studies the dynamic performance of the SCO 2 recompression cycle of a 20 MW sodium-cooled reactor with Simulink software.

System Configuration and Modeling Descriptions
The power generation system described in this paper adopts the SCO 2 recompression cycle, and the configuration is shown in Figure 1. The overall system includes a turbine (T), main compressor (MC), re-compressor (RC), high-temperature and low-temperature recuperators (HTRs and LTRs), intermediate heat exchanger (IHE), pre-cooler (PC), surge tank (ST), and heater. As shown in Figure 1, the high-temperature and high-pressure CO 2 from the IHE enters the turbine to do work and generates electricity. The exhausted CO 2 from turbine flows into the HTR and LTR and releases heat. The low-pressure fluid at the LTR outlet is then split into two streams: one stream is cooled in the PC and then compressed in the MC, and the other stream is directly compressed in the RC. Two streams converge at the high-pressure side outlet of LTR and are sent to the HTR and IHE for heating.
In order to meet the flexible control requirements, the MC, RC, and turbine adopt a split-axis layout. At the same time, the ST is arranged between the PC and the MC to facilitate the control of the MC inlet pressure and buffer the pressure change. The ST is connected with a CO 2 storage tank, and the MC inlet pressure is maintained above the critical point by controlling the amount of CO 2 in the ST.
This paper focuses on exploring the dynamic characteristics of the SCO 2 recompression cycle, and the mathematical model of each component is developed to reflect the dominant characteristics of the cycle. The model of each component could be able to realize independent physical functions and be mathematically independent. At the same time, there are clear boundaries and data interfaces between the modules. The internal parameters of the models can be changed according to actual needs, and the impact of each parameter on the overall circulatory system can be analyzed. Based on the above requirements, this paper develops the dynamic model of the SCO 2 system based on the Simulink platform, which is also convenient for system control strategy research.

Heat Exchanger Model
The heat exchangers involved in the SCO 2 power cycle, including HTR, LTR, PC, and IHE, adopt printed circuit heat exchangers (PCHEs), which have great advantages in dynamic modeling due to the regular layout structure. The flow and heat transfer model of PCHEs are based on the three conservation equations of mass, momentum, and energy, which can be simplified to one-dimensional models. The structure of PCHE is shown in Figure 2a. In this section, PCHE is divided into several parts, and the flow and heat transfer models are established for each part, as shown in Figure 2b. Assuming that the fluid in each part is incompressible, the steady-state model is adopted to describe the flow process. In each part, the lumped parameter method is applied to the calculation and the fluids and wall are treated as zero-dimensional nodes. The heat transfer processes between the cold fluid and the wall, the wall, and the hot fluid are considered, while the heat transfer in the direction of the fluid flow is ignored to ensure the accuracy and simplification of the calculation process. The specific calculation process is as follows [32,33].
Hot fluid: Cold fluid: Wall: The heat transfer coefficient calculation methods of different heat exchangers are selected based on the fluid type and Reynolds number (Re). For example, the Gnielinski formula is applied to calculate the Nusselt number (Nu) of the cold and hot side fluids in HTR and LTR, as shown in Equations (4) and (5). The cold and hot side of HTR and LTR are both CO 2 , so the Nusselt number (Nu) can be calculated by the Gnielinski formula. Thus, the heat transfer coefficient can be obtained. The specific heat capacity of the tube wall is determined by the material type: where Pr is the Prandtl number of the fluid. The pressure drop of each part can be determined by Equation (6). The pressure drop is proportional to the square of the flow rate and inversely proportional to the fluid density: where γ is the effective loss coefficient, which is determined by the pipeline characteristics and can be calculated with the parameters of the design condition.

Compressor Model
The dynamic response time of different processes and equipment in the SCO 2 power cycle is significantly different. The devices with a long response time adopt dynamic-state models, such as recuperators, coolers, etc. On the contrary, the devices and thermal processes with a short response time adopt steady-state models, such as turbine, compressors, working fluid separation, and mixing processes.
The compressor outlet parameters are determined by the inlet parameters and isentropic efficiency. Taking the main compressor (MC), for example, the specific calculation is shown as follows [34]: where h MC,in is the MC inlet-specific enthalpy, kJ/kg; h MC,out is the MC outlet-specific enthalpy, kJ/kg; h MC,out,is is the isentropic-specific enthalpy of the MC outlet, kJ/kg; and η MC,is is the isentropic efficiency of the MC. The outlet pressure and isentropic efficiency of the compressor are provided by the performance prediction curves. The curves describe the relationship between the outlet pressure and efficiency with the mass flow rate and rotating speed. The main compressor performance prediction curves are shown in Figure 3. The re-compressor performance prediction curves are shown in Figure 4, which are obtained through the simulation results based on the compressor design value [35].  It can be seen from the recompression configuration in Figure 1 that the low-pressure side outlet working fluid of the LTR is divided into two streams: one stream is cooled in the PC and then compressed in the MC, while the other stream is directly compressed in the RC. The split ratio (SR) is used to express the share of flow entering the recompressor, as shown in Equation (11): Therefore, the power consumption of the main compressor is: The power consumption of the recompressor is: where m is the mass flow rate of the circulation system, kg/s; W MC is the MC power consumption, kW; h RC,in is the RC inlet-specific enthalpy, kJ/kg; h RC,out is the RC outletspecific enthalpy, kJ/kg; and W RC is the RC power consumption, kW.

Turbine Model
Similar to the compressor model, the turbine outlet parameters are also determined by the inlet parameters and isentropic efficiency. The specific calculation is as follows [29]: s tur,in = s(t tur,in , p tur,in ) (15) h tur,out,is = h(s tur,in , p tur,out ) where h tur,in is the turbine inlet-specific enthalpy, kJ/kg; h tur,out is the turbine outlet-specific enthalpy, kJ/kg; h tur,out is the turbine outlet isentropic-specific enthalpy, kJ/kg; and η tur is the isentropic efficiency of the turbine. The turbine output power is calculated as follows: where m tur is the turbine flow rate, kg/s; and W tur is turbine output power, kW. The turbine performance prediction curves are shown in Figure 5 [35].

Surge Tank (ST) Model
In this section, considering the rapid response of pressure, the volume in the entire loop composed of HTR hot end-LTR hot end-PC hot end-ST is lumped into the surge tank to calculate the pressure change. Additionally, it is considered that the working fluid flow rate before and after the HTR hot end, LTR hot end, and PC hot end remains unchanged.
According to the mass conservation equation, the density change of the working fluid in the surge tank is as follows: where ρ ST is the working fluid density in the surge tank, kg/m 3 ; m HTR,h is the working fluid flow rate at the HTR hot end, kg/s; m CFP is the flow rate from the CO 2 filling pump (CFP) into the ST, kg/s; m PRV is the flow rate of the pressure-reducing valve (PRV), kg/s; V HTR is the working fluid volume at the HTR hot end, m 3 ; V LTR is the working fluid volume at the LTR hot end, m 3 ; V PC is the working fluid volume at the PC hot end, m 3 ; and V ST is the ST volume, m 3 .
According to the energy conservation equation, the specific enthalpy change of the working fluid in the surge tank is shown as below [36]: where h ST is the specific enthalpy of the CO 2 in ST, kJ/kg; h PC is the specific enthalpy of the working fluid at the PC hot end, kJ/kg; and h CFP is the specific enthalpy of the CO 2 at the CFP outlet, kJ/kg. According to the density ρ ST and specific enthalpy h ST , the temperature T ST and pressure P ST of the working fluid in surge tank can be obtained by consulting the REFPROP software.

Working Fluid Mixing Process
In the recompression cycle system, there is a node where the working fluid at the LTR cold end is mixed with that at the outlet of RC. According to the equations of mass conservation and energy conservation, the parameters of the mixed working fluid can be calculated as follows [37]: where h HTR,cin is the specific enthalpy of the mixed working fluid, kJ/kg; h LTR,cout is the specific enthalpy of the working fluid at the outlet of the LTR cold end, kJ/kg; and h RC,out is the specific enthalpy of the working fluid at the RC outlet, kJ/kg. According to the pressure and specific enthalpy, the temperature of the mixed working fluid can be obtained by consulting the REFPROP software.

Model Verification
A series of comparative calculations were carried out to verify the models. Firstly, the main component models were verified with the design values to make sure that the models are built correctly. All component models were connected based on the CO 2 flow sequence and system structure, and constitute the overall dynamic model of the power generation system. Then, based on the system design value, we carried out the system steady-state simulation.

PCHE Model Verification
The structural parameters of the PCHE are shown in Table 1. In the Simulink model, the PCHE was divided into five nodes, and the model was built with the Simulink graphical modeling tool based on the above-mentioned unsteady heat transfer relationship. We obtained the calculation results with the structural parameters and the thermal parameters based on the Simulink model. According to the comparison with the design value, the model calculation accuracy is reliable. As shown in Table 2, the outlet temperature of the hot side error is as small as 0.065%, while that of the cold side is larger, with an error of 0.25%. The results prove that the model has high accuracy. We also performed the disturbance analysis, and the results were compared with other model data to verify the dynamic simulation capacity of the model. At 300 s, the inlet temperature of the hot fluid increased and decreased by 30 • C, respectively. The characteristics of the outlet temperature of the cold and hot fluids were analyzed and compared with the literature [38]. The comparison results are shown in Figure 6: The authors of [38] carried out the PCHE thermodynamic analysis and performance prediction of a 1000 MW SCO 2 coal-fired power plant with CFD mesh simulation. Although the background and method were different from this paper, the model establishment was similar, being based on the basic conservation equations. As can be seen from the comparison, when affected by the disturbance of the step change in the hot fluid inlet temperature, the changes in the outlet temperature of the cold and hot fluids show similar dynamic characteristics: (1) When the inlet temperature of the hot fluid is subject to a step change, the outlet temperature of the cold fluid responds faster. At the same time, the outlet temperature of the cold fluid varies widely, which is equivalent to the change of the disturbance temperature; (2) The outlet temperature of the hot fluid is less affected by the inlet temperature and changes slowly after being disturbed. The temperature change is smaller as well.

Compressor Model Verification
According to the mathematical calculation formula and performance curve of the compressor in the steady state, the main compressor and recompressor models were built in Simulink. The calculation of the performance curve in the model adopted the interpolation calculation of the two-dimensional Look up table. The steady state calculation results were compared with the design value to verify the model accuracy. As shown in Table 3, the main compressor outlet temperature error is as large as 1.1%, which is considered acceptable.

Turbine Model Verification
For the turbine, the model build process is similar to the compressor model. The design value was also applied to the verification of the model and the comparison results are shown in Table 4. As can be seen, the largest error is the turbine power with a value of 0.74%, which is totally acceptable. Thus, it can be considered that the model is reliable.

System Verification
The closed-loop system model was debugged, and the initial operating parameters were calculated. With the problems of the algebraic loop and initial values solved, the SCO 2 recompression power system dynamic model can be run stably. The Simulink configuration diagram is shown in Figure 7. Similar to the system configuration and modeling description in Section 2, the Simulink configuration also mainly includes the turbine (T), main compressor (MC), recompressor (RC), high-temperature and low-temperature recuperators  The heat duty was input to the system to simulate the heater in this paper at a value of 64.3 MW. The SCO 2 in the IHE absorbs the heat and flows into the turbine to do work. After the heat release in the HTR and LTR, with a given split ratio (SR) at 0.3852, about 38.52% of the SCO 2 flows into the RC to be compressed directly while the remaining is cooled in the PC and compressed in the MC. Then, the two streams converge at the outlet of the high-pressure side of LTR and continually flow into the HTR and IHE. Thus, an integral process is completed. The recompression model was verified under design condition. The input parameters of the design condition are shown in Table 5, and the model results and error analysis are shown in Table 6. As shown in Table 6, the node number corresponds to the note number in Figure 1 in Section 2. Specifically, node 1 and node 2 are the turbine inlet and outlet, and node 3, node 4 , and node 5-6 refer to the flow process from the LTR inlet passing through PC to the MC outlet. The split SCO 2 to the RC outlet is noted by node 4 and node 5 . The two streams concentrate in the HTR inlet noted as node 8 and then flow into IHE at node 9. It can be seen from the table that the maximum error appears on the HTR inlet temperature at a value of 0.56% with a temperature difference of 2.2 • C. It is easy to note that the MC inlet temperature is kept at 304.4 • C, which is a benefit of the control of the ST. The temperature of the IHE inlet and the turbine is lower by about 1 • C, while the SCO 2 temperature from the inlet of LTR to the inlet of the HTR high-pressure side is generally higher than the design value. In the preliminary analysis, the split ratio may contribute to the error and the strong coupling of the closed-loop system would also amplify the errors of the heat transfer coefficient and physical property. However, the maximum error is only 0.56%, showing that the calculation results of this model can be considered reliable.

Experiment Verification Plan
Considering the study of the CO 2 Brayton circle is not mature enough, the detailed experimental data of a real power plant are very hard to obtain to support the experiment verification. Actually, our laboratory is preparing to carry out an experiment on the supercritical carbon dioxide system. The experimental loop is based on a simple regenerative cycle, and the schematic diagram is shown in Figure 8. The experiment loop includes a turbine, compressor, printed circuit heat exchanger (PCHE), intermediate heat exchanger (IHE), cooler, surge tank (ST), and heater. The design output power scale of the cycle is 1 MW and the circle efficiency is as large as 20%. The designed maximum temperature is 500 • C at the outlet of IHE, and the maximum pressure is 20 MW at the compressor outlet. At present, the experiment has entered the equipment manufacture stage, and the experiment will be carried out within one year. At this time, we will validate the developed equipment models and circle system based on the experimental results. In fact, the system model we developed will also serve the calculation of the experiment. Then, we will first validate the steady state operation system. Additionally, a series of transient-condition verification will also be carried out, including the working load increase and decrease process. Especially, the working state at a very low power state will also be performed, which will be our main focus.

System Dynamic Simulation of the Switching Working Condition
This section mainly studies the 20 MW-SCO 2 system dynamic characteristics in the process of switching between various operating conditions of 100%, 80%, 60%, 40%, and 20% load. The thermal parameters of each key node of the system under non-design operating conditions were estimated according to the performance curve of the equipment. The system mass flow rate was adjusted by changing the opening of the control valve in front of the turbine, so as to achieve the purpose of a variable load. Each working load had the corresponding compressor speed, heat duty, cooling water mass flow rate, and opening degree of the control valve, as shown in Table 7. In the simulation, the system initially operated stably under 20% operating conditions. When increasing the working load to the 40% load condition, the control valve opening degree was unchanged, and the compressor rotating speed, heat duty, and cooling water mass flow increased in a ramp manner. Then, the compressor and cooler maintained at the design condition and the control valve opening degree were increased continually to meet the load increase demand as shown in Figure 9. The system dynamic characteristics of the load-increasing process are shown in Figure 10. As can be seen, this operating mode can switch operating conditions by adjusting the amount of turbine distribution, while the compressor operating environment remains unchanged. This method can avoid major changes in the compressor operating environment during variable operating conditions and is beneficial to the stable operation of the system. At the same time, this method has a fast rate of changing working conditions, which can meet the demand of a fast changing load. However, this method has certain limitations for a variable load: when the load demand is below 40%, the compressor speed and cooling water mass flow experience great changes. If switching working conditions at this time, the system parameters will fluctuate greatly, and the stability time will be long. Furthermore, the compressor power consumption has to be maintained at a relatively high level, so the cycle efficiency is very low when operating at a low load.
The load-decreasing process is the opposite of the load-increasing process as shown in Figure 11, and the system dynamic characteristics are shown in Figure 12. Assuming that the system initially operates stably under 100% operating conditions, the compressor and cooler operating conditions remained unchanged at first when reducing the load. The turbine power reduced with the decrease of the turbine gas flow rate and heat duty. When the load was below 40%, reducing the compressor rotating speed, and the cooling water mass flow rate at the same time caused fluctuation of the outlet temperature in the heat exchanger and displayed a longer stable time.

Conclusions
In this study, Simulink software was applied to develop a dynamic model of the 20 MW-SCO 2 recompression cycle modularly. Corresponding component models were developed. A series of comparative calculations were carried out to verify the models. Additionally, the switching condition test of the 20 MW-SCO 2 recompression cycle was conducted. The main conclusions of the study are as follows:

•
The models of a 20 MW-SCO 2 recompression cycle were developed based on Simulink software, including heat exchanger models, compressor models, turbine models, and surge tank models.

•
The developed main component models were verified with the design values. The comparison results show that calculation errors are acceptable. The disturbance analysis results of the PCHE model were compared with PCHE in 1000 MW S-CO 2 coal-fired power plant data to verify the dynamic simulation capacity of the model. Additionally, the prediction trends were consistent.

•
The SCO 2 recompression power system was built with the developed models and the dynamic model ran stably and the maximum error was 0.56% compared to the design value.

•
The system dynamic characteristics were studied in the process of switching between various operating conditions of 100%, 80%, 60%, 40%, and 20% load. The analysis is instructive for the system operation and control.
Based on the analytical investigations, the results are very promising in that the developed model has the ability to evaluate the thermal-hydraulic characteristics of the SCO 2 Brayton circulation system, but further experimental data are still needed to validate the code.