Supercritical $CO2$ ($sCO2$) is a promising working fluid for future high-efficiency power conversion cycles. In order to develop these cycles, it is necessary to understand supercritical and two-phase $CO2$ flow. This paper presents a methodology for the computational fluid dynamic (CFD) simulation of $sCO2$ flowing through a restriction under a wide range of flow conditions. Under an accidental situation, such as a pipe break, the inventory of $sCO2$ leaks out through a small restriction. In this research, we use circular and annular orifices to mimic the behavior of restrictions. As the atmospheric pressure is much smaller than the operating pressure in the pipe, a two-phase choked flow will happen. Such behavior is considered in the simulation. The homogeneous equilibrium model (HEM) is employed to model the two-phase state. To correctly simulate the behavior of the power cycle under this accidental scenario, the inventory leakage rate should be calculated precisely. Therefore, at the current state, this study only focuses on the prediction of mass flow rate through orifices.

## Introduction and Motivation

A supercritical fluid is defined as a fluid in a state where pressure and temperature are above its critical point [1]. For $CO2$, the critical point is 7.38 MPa at 31.1°C. For water, the critical point is 22.06 MPa at 373.95°C. When a fluid is supercritical, it does not exhibit a saturation point where two phases can be distinguished from one another. However, supercritical fluid does exhibit a pseudocritical point, around which there is a dramatic change in its thermodynamic and transport properties, e.g., density and specific heat as seen in Fig. 1.

Supercritical water has been considered for future nuclear power plants for many years [2–6]. More recently, other fluids like $C2$ and helium have also been considered as candidates [7–10]. The main advantage of using supercritical fluid in a power cycle is that the appearance of a two-phase condition is avoided in the major components as shown in Figs. 2 and 3. The single-phase situation can simplify the design of components such as the compressor, turbine, and heat exchanger. With $sCO2$, the cycle can achieve high efficiency ($\u223c45%$) by utilizing the properties’ change near the critical point with relatively low turbine inlet temperature [7,9].

In an accident scenario of $sCO2$ power cycle, such as a pipe break, the inventory of $sCO2$ leaks out of the system. In order to simulate the behavior of the system under such situation, it is very critical to determine the leakage rate of the inventory. From the traditional approach, the ideal nozzle isentropic model, which is discussed with more detail in Section 5.3, combined with a discharge coefficient is used to evaluate the leakage rate. However, to the authors’ knowledge, the discharge coefficient is significantly sensitive to the geometry of the restriction. As a result, the computational fluid dynamic (CFD) method is used to consider the influence of geometry. In this research, we combine the properties of $sCO2$ and the open-source CFD code OpenFOAM to calculate the leakage rate. The homogeneous equilibrium model (HEM) is used to calculate two-phase properties, as two phases may appear in some situation.

## Previous Work

Facilities for supercritical fluid blowdown experiments were built to simulate pipe break scenarios. Experiments of this type using $sCO2$ were performed by Mignot [11], Chen [12], and Liu [13]. Experiments using supercritical water in choked flow were conducted by Chen [14] and Muftuoglu [15]. In UW-Madison, Rodarte [16], Edlebeck [17], and Wolf [18] constructed and improved the experiment facility to measure $sCO2$ flow rate through orifices. Experimental data from this facility are used in this work to compare with the simulation results.

Wallis [19] provides a thorough summary of modeling techniques for two-phase critical flow. Hesson and Peck [20] conducted experiments on $CO2$ two-phase choked flow and developed a correlation to predict their experimental data. However, their experiment and model had only been tested for saturated liquid, saturated vapor, and two-phase region. Simoneau and Hendrick [21] used a homogeneous model to develop general charts that can be used to calculate choked two-phase flow for a simple cryogenic liquid. None of the papers referenced show results in the supercritical region. Zhang and Yang [22,23] conducted tests in which the inlet condition was $CO2$ in the supercritical region and developed a model to predict critical mass flow rate using a one-dimensional model. Their model has about 10% error compared to their experimental data.

Commercial CFD codes have been used to simulate supercritical fluids. Yadav [24] used ANSYS-FLUENT [26] to simulate $sCO2$ in a natural circulation loop. ANSYS-FLUENT has internally implemented the property routine REFPROP [25] into its code. However, users can only use REFPROP in single-phase or supercritical states with ANSYS-FLUENT. Van Abel [26] also used ANSYS-FLUENT to simulate heat transfer in printed circuit heat exchangers for $sCO2$. With a proper turbulence model and careful geometric description, his simulation shows very good agreement with the experiment. Suo-Anttila [27] used a commercial CFD code called C3D and implemented an interface with REFPROP to compute properties. Qiu [28,29] used the homogeneous equilibrium assumption to study the condensation process for regular and retrograde fluid. His simulation shows very good agreement with the experimental data. All the research discussed above is limited to a single phase and the supercritical region, and none of their research is extended to the two-phase domain. Fairweather [30] developed a CFD model for $sCO2$, which is capable of simulating two-phase conditions. They used the HEM in their model and simulated a pipe break problem in order to study an accident scenario. Their CFD model was compared to their experiment data and the agreement was good; however, only the temperature distribution data were provided and not the leakage rate.

## Numerical Methodology

As fluids discharge from a supercritical state, two phases may appear. In order to simulate two-phase flow, the CFD code cannot use temperature and pressure as the only primitive variables to determine other properties. This requirement eliminates some commercial CFD codes for solving this problem. Most commercial CFD software products do not provide access to the source code, which would be required to make changes in the underlying equations and property data. For this reason, we decided to use the open-source CFD code OpenFOAM [31] for this research. Part of the code is modified to handle $CO2$ in the supercritical and two-phase state. The solver used in this research is designed for the steady-state Reynolds-averaged Navier–Stokes (RANS) equations; a quick review of the steady-state RANS equations is helpful for further discussion.

### Solver.

The *rhoSimplecFoam* [31] solver in OpenFOAM is used as the prototype for the modified solver. This solver uses the *SIMPLEC* [32] algorithm to solve Reynolds-averaged Navier–Stokes equations for steady-state compressible flow. As the SIMPLEC algorithm is a standard algorithm for CFD, readers can find many references that have discussed this algorithm in detail. However, the original solver is designed using an ideal gas assumption. This study modified the solver to let it handle $CO2$ in supercritical and two-phase states.

### Modification of Energy Equation.

*rhoSimplecFoam*and uses the same algorithm for pressure-velocity coupling; however, some modifications are made to the energy equation. Equation (6) shows the energy equation in the unmodified solver. In the modified solver, the energy equation is changed to the form of Eq. (7) to avoid using $\alpha $, which is not defined in the two-phase region

### Real Fluid Model.

REFPROP [25] was initially used to provide the $CO2$ properties for the simulation. However, it was found that the use of REFPROP causes simulations to run extremely slow. Fluid properties are currently calculated using the FIT software library published by Northland Numerics [33], which provides an interpolated representation of properties; however, the underlying property data are still obtained from REFPROP. In this way, the FIT code calculates fluid properties much faster than REFPROP with acceptable error. Figure 4 shows a comparison of the FIT and REFPROP property calculation for $CO2$-specific heat. Each plot is based on 1 million points, distributed evenly from 240 to 1500 K and from 1 kPa to 100 MPa. As we can see, the maximum relative difference for specific heat is less than 1e-5. Comparisons with other properties show the same agreement.

### Homogeneous Equilibrium Model.

### Geometry and Boundary Conditions.

Since both circular and annular orifices are axisymmetric, a two-dimensional axisymmetric geometry method is used to save computational time. OpenFOAM, unlike ANSYS-FLUENT, does not have the ability to solve 2D problems in cylindrical coordinates because every term of Navier–Stokes equations in OpenFOAM is designed for the Cartesian system [31,34]. Therefore, the OpenFOAM community uses a pseudoaxisymmetric method for 2D axisymmetric problems [31,34]. The computational domains for circular and annular orifices are shown in Figs. 5 and 6. However, readers should keep in mind that these figures are not created in scale.

There are several boundary conditions used in simulations as shown in Figs. 5 and 6. The “wall” boundary condition means nonsplit and nonpenetration for velocity, and zero gradient for pressure. Due to the adiabatic assumption in this process, this study assumes zero heat flux at the wall, which applies zero gradient of temperature. This is achieved by the combination of zero gradient for both pressure and enthalpy. The “wedge” boundary condition is recommended by the OpenFOAM User Guide [31] for pseudoaxisymmetric geometry when there is only one level of cells in the corresponding direction. It assumes the adjacent cell is the cell itself. Therefore, the “wedge” boundary condition is actually a specialized periodic boundary condition, which is designed for pseudoaxisymmetric geometry. The “empty” boundary condition makes the cell faces on that patch have a zero area, which means that there are no fluxes, including momentum flux and energy flux. IThe OpenFOAM User Guide [31] recommends to use the “empty” boundary condition at the axis of a pseudoaxisymmetric geometry. The cells on the axis degrade to wedge-shaped blocks as the vertices on the axis collapse to one. However, other cells are still hexahedral.

### Meshing Sensitivity.

We use the OpenFOAM default meshing utility *blockMesh* [31] to create meshes. The annular orifice geometry with different meshing numbers was simulated and compared to demonstrate meshing sensitivity of leakage rate. The dimension of the annular orifice is shown in Fig. 6. The inlet condition is 10 MPa at $475\u2009\u2009kg/m3$. The standard k-epsilon turbulence model is used. More discussion on the choice of turbulence models is presented in Section 3.7. This study is only concerned about the prediction of mass flow rate through orifices. Figure 7 shows the choked mass flow rate prediction for different meshing numbers. From Fig. 7, it is sufficient to conclude that the meshing is refined enough for this research. The meshing is maintained at the highest level as shown in Fig. 7.

### Turbulence Model.

In order to study the turbulence effect on simulation, several turbulence models should be applied and the results compared. The solver employed herein is designed for the RANS turbulence model; as a result, only RANS turbulence models are considered. He [35] and Yoo [36] discussed the choice of turbulence models for CFD simulation for supercritical fluid. However, their discussion is based on relatively low Re situations, where heat transfer is a major concern. In this research, all walls are assumed to be adiabatic, and only the mass flow rate through orifices is measured in the experiment and extracted from the simulation. Moreover, the Re number is about 200,000, which is much larger than the situation discussed in [35] and [36]. It is not clear theoretically which turbulence model is best to use for our research problem; the standard $k$-epsilon model [37], which is designed for high Re situations, is apparently our best choice. The $k$-omega SST model [38] is also used to compare the to the results of the standard $k$-epsilon model. The equations and coefficients used in OpenFOAM for these two models are listed below.

The comparison is based on an annular orifice with dimensions as shown in Fig. 6; the data are shown in Fig. 8. The two previously introduced turbulence models are used in the simulation under the same inlet condition of 10 MPa at $475\u2009\u2009kg/m3$. Figure 8 shows that the leakage rate predictions from these two turbulence models are very close. It seems that different turbulence models do not predict different leakage rates. At the current stage, a more detailed analysis of turbulence modeling is not performed, but will be performed in future work. Therefore, currently, the standard $k$-epsilon turbulence model is used for all the simulations herein.

There have been many studies on the constant $Prt$ assumption in supercritical fluid flow [35,39–41]. It is believed that the current constant $Prt$ assumption is not valid in supercritical fluid flow, and the value of $Prt$ deviates from unity. However, in [41], the recommended value of $Prt$ does not significantly deviates from unity. As a result, three different constant $Prt$ values are tested to see its influence on leakage rate. These values are 0.9, 1.0, and 1.1. Figure 9 shows the results of different $Prt$ at the same inlet condition (7.7 MPa, $498\u2009\u2009kg/m3$) for a circular orifice. As seen in Fig. 9, different $Prt$ numbers give the same prediction of leakage rate through a circular orifice. Thus, the assumption of $Prt=1.0$ is used in all simulations in this research. In this problem, the heat transfer due to turbulence is much smaller than the energy transport by flow. Therefore, different $Prt$ numbers do not change the leakage rate significantly.

## Test Facility

The test setup that provided the experiment data is depicted in Fig. 10. The setup was used to measure the leakage rate and pressure drop for circular and annular orifices under $sCO2$ conditions. A single-stage compressor made by Hydro-Pac is used to raise the pressure of $CO2$ above the critical pressure (7.38 MPa). The outlet of the compressor is connected to a buffer tank, which is a large-volume tank used to minimize the pressure fluctuations due to the compressor cycling, and enhances the stability of the whole system. After the compressor, the flow passes through precooler and preheater stages. The precooler and preheater are automatically controlled to cool or heat the fluid to the target temperature. Before the test section, the flow goes through a flow meter where both the density and mass flow rate are recorded. Thermocouples and pressure transducers are placed before and after the test section to measure the upstream and downstream conditions, corresponding to states 1 and 2 in Fig. 11. After the flow passes through the test section, it enters a reservoir tank to help attain a stable downstream condition. The flow then returns to the inlet condition of the compressor. A description of the conditions typically viewed at various points around the loop during a standard experiment can be seen in Fig. 11. This figure shows a temperature-specific entropy diagram with points labeled with respect to Fig. 10. The uncertainty analysis of this test facility is discussed with great detail in Refs. [16,44,45]. Readers who want to know more about this test facility can refer to these references.

## Results and Discussion

The CFD simulation method discussed in the previous section was used to simulate the flow of $sCO2$ through circular and annular orifices. The same geometries were tested in the experimental facility to provide validation data. After the data are presented, the ensuing discussion explores the working region of the proposed methodology.

### Circular Orifice.

This section summarizes the simulation data for the circular orifice. The circular orifice has a diameter of 1.014 mm and a length of 3.2 mm. Twelve inlet conditions have been simulated. Figure 12 shows the inlet conditions on a $T\u2013s$ diagram. All of the inlet conditions correspond to a supercritical state; one of the inlet conditions (pressure is 7.7 MPa, density is $111\u2009\u2009kg/m3$) is far away from the two-phase dome while the rest are located just above the two-phase dome. As a result, the flow exhibits flashing when the outlet condition is in the two-phase region.

Figure 13 shows one set of data for one inlet condition (7.7 MPa, $498\u2009\u2009kg/m3$, the point with arrow in Fig. 12); data from other inlet conditions show similar agreement with the experimental measurement. In Fig. 13, besides experiment and simulation data, the mass flow rate prediction from the ideal nozzle isentropic model is also presented for comparison (with the notation “IS”). The isentropic model is introduced in more detail in Section 5.3 as it is used as a theoretical model for discussion. In Fig. 13, the mass flow rate prediction from the simulation is within a few percent compared to the experimental measurements. Figure 14 presents the phase distribution for inlet condition of 7.7 MPa, $498\u2009\u2009kg/m3$, and outlet of 5.2 MPa. The trajectory line of this discharge process is also shown in Fig. 21. However, in Fig. 14, the quality is defined as 1.0 when fluid is in the single-phase state, including the supercritical phase. When fluid is in the two-phase state, the quality returns to its original definition. In Fig. 15, we present the simulation data versus experiment data for all the tested inlet conditions in Fig. 12. Figure 15 shows that the maximum leakage rate difference between simulation and experiment is less than 5%.

### Annular Orifice.

Simulations and experiments are also performed on three annular orifices with different lengths. Annular orifices with small, medium, and large length-to-diameter ratios are simulated. The results presented show that as the length of the annular orifice increases, friction becomes more important in the simulation. As a result, the effect of roughness on the simulation is investigated at the same time. Table 3 gives the geometric and inlet conditions of these annular orifices. The roughness effect on the circular orifice is negligible as its length-to-diameter ratio is small.

In Figs. 16 and 17, the number after the “OF” legend represents the roughness applied in the simulation. As can be seen, roughness plays an important role for the long annular orifices, but it is less significant in the short annular orifices. The best estimate for the roughness is approximately 5 μm. The measured roughness of these annular orifices is about 4.6 μm, which is very close to the estimated value. All the data from simulations with roughness of 5 μm are compared with experimental measurements in Fig. 18. Like the results discussed in the previous section, the maximum difference is about 5%.

### Discussion.

From the data presented, it is clear that HEM works very well for the study problem. In this section, we conduct some qualitative analysis to explain the applicability of HEM to this problem. First, we identify the flow regime. Theoretically, HEM should work in the regime of dispersed flow, i.e., the two phases are well mixed. The goal is to show the flow is dispersed in this problem. We then discuss the applicable working region of the proposed methodology with some simple theoretical models.

Let us consider an example case: $sCO2$ discharges from 9 MPa and $338\u2009\u2009kg/m3$ to 6 MPa, with calculated velocities for the two phases that are $ug=104\u2009\u2009m/s$ and $ul=68\u2009\u2009m/s$. With the correction factors provided by Weisman [42,43], the corrected velocities are $ug/\varphi 1=2205\u2009\u2009m/s$ and $ul/\varphi 2=68\u2009\u2009m/s$. When we apply the result in Fig. 19, a dispersed flow pattern results. This analysis is qualitative, but the magnitude should be similar for real situations.

In the previous discussion, we qualitatively prove that the dispersed flow pattern applies to this problem, which proves the applicability of HEM. In this part, the working region of HEM is discussed. In this study, only the mass flow rate through orifices is measured and compared; as a result, the averaged properties, such as density, viscosity, conductivity, will influence the result. However, under the current high Re situation, the flow is highly turbulent, which means turbulent viscosity and conductivity are much larger than their thermal counterparts. Therefore, the averaging method for viscosity and conductivity is less important. As a result, the density is used as the only criteria to discuss the working region of HEM in this session.

Here, we investigate the density prediction as a function of quality at different pressure values. Figure 20 uses Eq. (28) to calculate the averaged density and compares it with HEM. Figure 20 shows that at a pressure of 7 MPa, the density calculations from other models are similar to HEM. However, at 4 MPa, the results from different models deviate more.

Based on the 10% criteria mentioned previously, for each slip ratio model we can draw a characteristic line in the two-phase domain that divides itself into two subdomains. In one subdomain, the density difference is more than 10%, and in the other, it is within 10%. The characteristic lines from these two slip ratio models, Moody’s and Fauske’s, are presented in Fig. 21. The characteristic line in Fig. 21 means that between this line and the saturation line, the density calculation from HEM is within 10% of other models (depending on what model is used). Thus, in this region, HEM should work pretty well. However, when the thermodynamic state is below this characteristic line, then the density calculated from HEM deviates from other models by more than 10%, which means HEM works poorly. Readers should keep in mind that, specific to our problem, once the flow gets choked, further reducing the downstream pressure will not affect the mass flow rate. That is to say, when the flow chokes, the thermodynamic state at the exit of the restriction remains at a fixed state. As long as this exit state is within the working region of HEM in Fig. 21, HEM should predict the mass flow rate quite well.

However, in our experiment and simulation, we found that for some inlet conditions that are close to the critical point, the choking exit state is actually below the characteristic line. In Fig. 21, there is a trajectory line representing the corresponding outlet and choking conditions for an inlet condition of 7.7 MPa and $498\u2009\u2009kg/m3$. As shown in Fig. 21, when the flow gets choked, the exit condition is actually below these two characteristic lines. That is why we emphasize that the analysis above is qualitative, not quantitative, and the transition of the HEM working region is smooth, not sharp.

From the preceding discussion, we can build the following qualitative concept for the working region of the proposed methodology:

- 1.
The flow pattern should be dispersed flow. However, the flow map used here is obtained from air-water flow in tubes, and it is not quite applicable to this problem. Currently, there is a lack of both experimental and numerical data for this specific situation. Future work includes finding a better way to determine the flow pattern for two-phase flow near the critical point of the fluid.

- 2.
The density calculation from HEM should be close to what is calculated from other more complicated models. However, in the analysis above, only the Moody’s and Fauske’s slip ratio models are applied, and the slip ratio is not known exactly. As a result, Fig. 21 only presents the HEM working region qualitatively, not quantitatively.

## Conclusion and Future Work

In order to simulate the $sCO2$ Brayton cycle under a pipe break scenario, it is important to calculate the leakage rate of the inventory. The leakage rate of $sCO2$ is sensitive to the upstream condition and restriction geometry. To take these effects into account, a method to simulate $CO2$ flow in supercritical and two-phase regions was developed with OpenFOAM. An interface was implemented in OpenFOAM to link a computationally efficient external property code FIT. The HEM assumption was made for two-phase states. This methodology was used to simulate $sCO2$ flow through circular and annular orifices. A detailed comparison between experimental and simulation results shows differences of less than 5% in the mass flow rate and proves the credibility of this proposed methodology.

A discussion of why HEM works well for this research is presented. We provide two arguments to explain this observation. First, the dispersed flow pattern applies to this problem; however, it needs additional work to be further proved. Second, there is a region in the two-phase dome where the density calculation from HEM is similar to other models. As soon as the thermodynamic state is within this region, HEM should work properly. However, this analysis is qualitative not quantitative. Future work will try to determine the working region of HEM with other more complicated models.

## Acknowledgment

The authors would like to thank the US DOE NEUP grant for supporting this research.

## Nomenclature

### Latin Symbols

- $Cp$ =
specific heat, $J/(kg\xb7K)$

- $h$ =
specific enthalpy, $J/kg$

- $k$ =
turbulent kinetic energy, $J/kg$

- $m\u02d9$ =
mass flow rate, $kg/s$

- $P$ =
pressure, Pa

- $Pk$ =
turbulent kinetic energy production rate, $J/(kg\xb7s)$

- $Prt$ =
turbulent Prandtl number

- PR =
pressure ratio, outlet pressure over inlet pressure, $Poutlet/Pinlet$

- Re =
Reynolds number

- $S$ =
slip ratio, gas velocity over liquid velocity, $ug/ul$

- $s$ =
specific entropy, $J/(kg\xb7K)$

- $T$ =
temperature, °C

- $t$ =
time, s

- $U$ =
velocity (vector), $m/s$

- $u$ =
velocity, $m/s$

- $x$ =
quality

- $\alpha $ =
thermal diffusivity, $m2/s$, or void fraction

- $\alpha t$ =
turbulence diffusivity, $m2/s$

- $\epsilon $ =
turbulent kinetic energy dissipation rate, $J/(kg\u2009s)$

- $\kappa $ =
thermal conductivity, $W/(m\u2009K)$

- $\upsilon $ =
specific volume, $m3/kg$

- $\mu $ =
dynamic viscosity, Pa s

- $\mu t$ =
turbulence viscosity, Pa s

- $\rho $ =
density, $kg/m3$

- $\varphi 1$ =
Weisman correction factor for gas phase

- $\varphi 2$ =
Weisman correction factor for liquid phase

- $\omega $ =
turbulence frequency, $1/s$

- choked =
choked condition

- effective =
combination of thermal and turbulence properties

- EXP =
data obtained from experiment

- $g$ =
gas

- HEM =
property from HEM model

- inlet =
condition at inlet

- $l$ =
liquid

- outlet =
condition at outlet

- OF
data obtained from simulation

- slip =
property from slip ratio model