Simulation of Supercritical CO2 Flow Through Circular and Annular Orifice

Yuan, Haomin; Edlebeck, John; Wolf, Mathew; Anderson, Mark; Corradini, Michael; Klein, Sanford; Nellis, Gregory

doi:10.1115/1.4029337

Supercritical ${CO}_{2}$ (⁠ ${sCO}_{2}$ ⁠) 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 ${CO}_{2}$ flow. This paper presents a methodology for the computational fluid dynamic (CFD) simulation of ${sCO}_{2}$ flowing through a restriction under a wide range of flow conditions. Under an accidental situation, such as a pipe break, the inventory of ${sCO}_{2}$ 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 ${CO}_{2}$ ⁠, 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.

Fig. 1

View large Download slide

Property change of ${CO}_{2}$ near the pseudocritical point at pressure of 8 MPa

Supercritical water has been considered for future nuclear power plants for many years [2–6]. More recently, other fluids like $C_{2}$ 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 ${sCO}_{2}$ ⁠, the cycle can achieve high efficiency (⁠ $\sim 45 %$ ⁠) by utilizing the properties’ change near the critical point with relatively low turbine inlet temperature [7,9].

Fig. 2

View large Download slide

Schematic of ${sCO}_{2}$ simple Brayton cycle

Fig. 3

View large Download slide

$T - s$ diagram of ${sCO}_{2}$ simple Brayton cycle

In an accident scenario of ${sCO}_{2}$ power cycle, such as a pipe break, the inventory of ${sCO}_{2}$ 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 ${sCO}_{2}$ 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 ${sCO}_{2}$ 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 ${sCO}_{2}$ 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 ${CO}_{2}$ 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 ${CO}_{2}$ 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 ${sCO}_{2}$ 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 ${sCO}_{2}$ ⁠. 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 ${sCO}_{2}$ ⁠, 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 ${CO}_{2}$ 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.

The continuity equation is

\nabla \cdot (ρ \cdot U) = 0

(1)

The momentum equation is

\nabla \cdot (ρ \cdot U \cdot U) - \nabla \cdot (μ_{effective} \cdot \nabla U) = - \nabla P

(2)

and the energy equation is

\nabla \cdot (ρ \cdot U \cdot h) + \nabla \cdot (ρ \cdot U \cdot \frac{U^{2}}{2}) - \nabla \cdot (ρ \cdot α_{effective} \cdot \nabla h) = 0

(3)

where

μ_{effective} = μ + μ_{t}

(4)

α_{effective} = α + α_{t}

(5)

\Pr_{t} = \frac{μ_{t}}{ρ \cdot α_{t}}

(6)

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 ${CO}_{2}$ in supercritical and two-phase states.

Modification of Energy Equation.

The modified solver inherits most of its code from rhoSimplecFoam and uses the same algorithm for pressure-velocity coupling; however, some modifications are made to the energy equation. Equation () shows the energy equation in the unmodified solver. In the modified solver, the energy equation is changed to the form of Eq. () to avoid using

α

⁠, which is not defined in the two-phase region

\nabla \cdot (ρ \cdot U \cdot h) + \nabla \cdot (ρ \cdot U \cdot \frac{U^{2}}{2}) - \nabla \cdot (ρ \cdot α_{effective} \cdot \nabla h) = 0

(7)

\nabla \cdot (ρ \cdot U \cdot h) + \nabla \cdot (ρ \cdot U \cdot \frac{U^{2}}{2}) - \nabla \cdot (ρ \cdot α_{t} \cdot \nabla h) = \nabla \cdot (κ \cdot \nabla T)

(8)

Real Fluid Model.

REFPROP [25] was initially used to provide the ${CO}_{2}$ 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 ${CO}_{2}$ -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.

Fig. 4

View large Download slide

Comparisons of FIT and REFPROP for $C_{p}$

Homogeneous Equilibrium Model.

In the HEM, all phases are assumed to have the same velocity, temperature, and pressure. With the HEM assumption, the state of two phases can be determined by pressure and enthalpy uniquely as shown in Eqs. () and (). The energy equation in Eq. () solves enthalpy. Then, quality is determined from the enthalpy value using Eq. (), which is then used to determine other properties

v_{HEM} = x \cdot v_{g} + (1 - x) \cdot v_{l}

(9)

h_{HEM} = x \cdot h_{g} + (1 - x) \cdot h_{l}

(10)

For transport properties (e.g., viscosity and conductivity), the average value for a two-phase state can be based on area, mass, or volume. The CFD model here obtains two-phase transport properties for viscosity and conductivity based on the mass average using Eqs. () and (). As the high Re (⁠

\sim 200,000

⁠) encountered in the problems, turbulence viscosity and diffusivity are much larger than their thermal counterparts. As a consequence, different averaging methods for transport properties do not produce different leakage rate predictions

μ_{HEM} = x \cdot μ_{g} + (1 - x) \cdot μ_{l}

(11)

κ_{HEM} = x \cdot κ_{g} + (1 - x) \cdot κ_{l}

(12)

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.

Fig. 5

View large Download slide

Computational domain for circular orifice

Fig. 6

View large Download slide

Computational domain for annular orifice

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 {kg / m}^{3}$ ⁠. 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.

Fig. 7

View large Download slide

Choked mass flow rate for different meshing numbers for annular orifice

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.

In standard

k

-epsilon model

\frac{D ρ k}{D t} = \nabla \cdot ((\frac{μ_{t}}{σ_{k}} + μ) \nabla k) + ρ P_{k} - \frac{2}{3} ρ k \nabla \cdot U - ρ ε

(13)

\frac{D ρ ε}{D t} = \nabla \cdot ((\frac{μ_{t}}{σ_{ε}} + μ) \nabla ε) + C_{ε 1} \frac{ρ P_{k} ε}{k} - (\frac{2}{3} C_{ε 1} + C_{ε 3}) ρ ε \nabla \cdot U - C_{ε 2} \frac{ρ ε^{2}}{k}

(14)

μ_{t} = C_{μ} \frac{ρ k^{2}}{ε}

(15)

α_{t} = \frac{μ_{t}}{ρ \cdot \Pr_{t}}

(16)

with coefficients provided in Table 1.

Table 1

Coefficients for standard $k$ -epsilon model in OpenFOAM

$C_{μ}$	$C_{ε 1}$	$C_{ε 2}$	$C_{ε 3}$	$σ_{k}$	$σ_{ε}$	$\Pr_{t}$
0.09	1.44	1.92	$- 0.33$	1.0	1.3	1.0

In the

k

-omega SST model

\frac{D ρ k}{D t} = \nabla \cdot ((α_{k} μ_{t} + μ) \nabla k) + ρ P_{k} - \frac{2}{3} ρ k \nabla \cdot U - ρ β^{*} ω k

(17)

\frac{D ρ ω}{D t} = \nabla \cdot ((α_{ω} μ_{t} + μ) \nabla ω) + γ \frac{ρ P_{k} ω}{k} - \frac{2}{3} ρ γ ω \nabla \cdot U - ρ β ω^{2} - ρ (F - 1) F

(18)

μ_{t} = \frac{ρ k}{ω}

(19)

α_{t} = \frac{μ_{t}}{ρ \cdot \Pr_{t}}

(20)

where

F = \frac{2 α_{ω 2} (\nabla k \cdot \nabla ω)}{ω}, α_{k} = F α_{k 1} + (1 - F) α_{k 2}, α_{ω} = F α_{ω 1} + (1 - F) α_{ω 2}, γ = F γ_{1} + (1 - F) γ_{2}, β = F β_{1} + (1 - F) β_{2}

with coefficients provided in Table 2.

Table 2

Coefficients for $k$ -omega SST model in OpenFOAM

$α_{k 1}$	$α_{k 2}$	$α_{ω 1}$	$α_{ω 2}$	$β_{1}$	$β_{2}$	$β^{*}$	$γ_{1}$	$γ_{2}$	$\Pr_{t}$
0.85034	1.0	0.5	0.85616	0.075	0.0828	0.09	0.5532	0.4403	1.0

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 {kg / m}^{3}$ ⁠. 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.

Fig. 8

View large Download slide

Short annular orifice data for comparison of standard $k$ -epsilon and $k$ -omega SST

There have been many studies on the constant $\Pr_{t}$ assumption in supercritical fluid flow [35,39–41]. It is believed that the current constant $\Pr_{t}$ assumption is not valid in supercritical fluid flow, and the value of $\Pr_{t}$ deviates from unity. However, in [41], the recommended value of $\Pr_{t}$ does not significantly deviates from unity. As a result, three different constant $\Pr_{t}$ 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 $\Pr_{t}$ at the same inlet condition (7.7 MPa, $498 {kg / m}^{3}$ ⁠) for a circular orifice. As seen in Fig. 9, different $\Pr_{t}$ numbers give the same prediction of leakage rate through a circular orifice. Thus, the assumption of $\Pr_{t} = 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 $\Pr_{t}$ numbers do not change the leakage rate significantly.

Fig. 9

View large Download slide

Circular orifice data for different $\Pr_{t}$ using standard $k$ -epsilon model

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 ${sCO}_{2}$ conditions. A single-stage compressor made by Hydro-Pac is used to raise the pressure of ${CO}_{2}$ 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.

Fig. 10

View large Download slide

Schematic diagram of ${sCO}_{2}$ test facility at UW-Madison

Fig. 11

View large Download slide

$T - s$ diagram of ${sCO}_{2}$ test facility at UW-Madison

Results and Discussion

The CFD simulation method discussed in the previous section was used to simulate the flow of ${sCO}_{2}$ 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 - s$ diagram. All of the inlet conditions correspond to a supercritical state; one of the inlet conditions (pressure is 7.7 MPa, density is $111 {kg / m}^{3}$ ⁠) 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.

Fig. 12

View large Download slide

Circular orifice test inlet conditions on $T - s$ diagram

Figure 13 shows one set of data for one inlet condition (7.7 MPa, $498 {kg / m}^{3}$ ⁠, 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 {kg / m}^{3}$ ⁠, 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%.

Fig. 13

View large Download slide

Circular orifice data for inlet condition of 7.7 MPa, $498 kg / m^{3}$

Fig. 14

View large Download slide

Quality distribution of circular orifice for inlet condition of 7.7 MPa, $498 kg / m^{3}$ ⁠, outlet condition of 5.2 MPa

Fig. 15

View large Download slide

Mass flow rate comparisons for circular orifice

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.

Table 3

Geometry parameter for annular orifices

	Inner diameter (in.)	Outer diameter (in.)	Length (in.)	Length/clearance	Test conditions
Short	0.118	0.1267	0.05	11.5	(10 MPa, $475 kg / m^{3}$ ⁠)
Medium	0.118	0.1267	0.1	23	(10 MPa, $325 kg / m^{3}$ ⁠), (10 MPa, $475 kg / m^{3}$ ⁠)
Long	0.118	0.1258	1.23	315	(10 MPa, $475 kg / m^{3}$ ⁠), (11 MPa, $498 kg / m^{3}$ ⁠)

	Inner diameter (in.)	Outer diameter (in.)	Length (in.)	Length/clearance	Test conditions
Short	0.118	0.1267	0.05	11.5	(10 MPa, $475 kg / m^{3}$ ⁠)
Medium	0.118	0.1267	0.1	23	(10 MPa, $325 kg / m^{3}$ ⁠), (10 MPa, $475 kg / m^{3}$ ⁠)
Long	0.118	0.1258	1.23	315	(10 MPa, $475 kg / m^{3}$ ⁠), (11 MPa, $498 kg / m^{3}$ ⁠)

$Clearance = (outer diameter - inner diameter) / 2$ ⁠.

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

Fig. 16

View large Download slide

Short annular orifice data for inlet condition of 10 MPa, $475 kg / m^{3}$

Fig. 17

View large Download slide

Long annular orifice data for inlet condition of 10 MPa, $475 kg / m^{3}$

Fig. 18

View large Download slide

Mass flow rate comparisons for annular orifice

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.

In the flow regime map from Weisman [42,43], as shown in Fig. 19, when the velocities for both liquid and gas are high (⁠

\sim 50 - 100 m / s

⁠), a dispersed flow pattern is the result. In the following, we show that the velocities for both liquid and gas are very high (⁠

\sim 50 - 100 m / s

⁠), and a dispersed flow pattern is suitable to describe the flow. The ideal nozzle isentropic model (Eqs. ()–()) combined with Moody’s slip ratio model (Eq. ()) is used to prove this argument qualitatively. Void fraction is calculated based on Moody’s slip ratio model using Eq. (). The total mass flux should be consistent by Eq. ()

h_{inlet} = enthalpy (P_{inlet}, s_{inlet})

(21)

h_{outlet} = enthalpy (P_{outlet}, s_{inlet})

(22)

ρ_{outlet} = density (P_{outlet}, s_{inlet})

(23)

h_{inlet} = h_{outlet} + u_{outlet}^{2}

(24)

S = {(\frac{ρ_{l}}{ρ_{g}})}^{1 / 3} = \frac{u_{g}}{u_{l}}

(25)

α = \frac{1}{1 + (\frac{1 - x}{x}) (\frac{ρ_{g}}{ρ_{l}}) S}

(26)

α ρ_{g} u_{g} + (1 - α) ρ_{l} u_{l} = ρ u

(27)

Fig. 19

View large Download slide

Weisman flow regime map for vertical flow [42,43]

Let us consider an example case: ${sCO}_{2}$ discharges from 9 MPa and $338 {kg / m}^{3}$ to 6 MPa, with calculated velocities for the two phases that are $u_{g} = 104 m / s$ and $u_{l} = 68 m / s$ ⁠. With the correction factors provided by Weisman [42,43], the corrected velocities are $u_{g} / ϕ_{1} = 2205 m / s$ and $u_{l} / ϕ_{2} = 68 m / 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.

In HEM, two phases are assumed to have the same velocity, which results in a lower averaged density. In the following part, two slip ratio models, Moody’s [44] and Fauske’s [45], are used to compare the averaged density with HEM. As HEM significantly deviates from other models, its performance should also degrade. To indicate this transition, this paper has arbitrarily chosen the criteria that when the averaged density from HEM deviates more than 10% from other models, HEM works poorly. However, this transition should be smooth, not sharp. Equation () is used to calculate the void fraction, Eq. () represents Fauske’s slip ratio model, and Eq. () is applied to calculate the averaged density of two phases

S = {(\frac{ρ_{l}}{ρ_{g}})}^{1 / 2} = \frac{u_{g}}{u_{l}}

(28)

ρ = α ρ_{g} + (1 - α) ρ_{l}

(29)

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.

Fig. 20

View large Download slide

Density ratio of other models compare to HEM

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.

Fig. 21

View large Download slide

Working region for HEM

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 {kg / m}^{3}$ ⁠. 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 ${sCO}_{2}$ Brayton cycle under a pipe break scenario, it is important to calculate the leakage rate of the inventory. The leakage rate of ${sCO}_{2}$ is sensitive to the upstream condition and restriction geometry. To take these effects into account, a method to simulate ${CO}_{2}$ 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 ${sCO}_{2}$ 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

$C_{p}$ =: specific heat, $J / (kg \cdot K)$
$h$ =: specific enthalpy, $J / kg$
$k$ =: turbulent kinetic energy, $J / kg$
$\dot{m}$ =: mass flow rate, $kg / s$
$P$ =: pressure, Pa
$P_{k}$ =: turbulent kinetic energy production rate, $J / (kg \cdot s)$
$\Pr_{t}$ =: turbulent Prandtl number
PR =: pressure ratio, outlet pressure over inlet pressure, $P_{outlet} / P_{inlet}$
Re =: Reynolds number
$S$ =: slip ratio, gas velocity over liquid velocity, $u_{g} / u_{l}$
$s$ =: specific entropy, $J / (kg \cdot K)$
$T$ =: temperature, °C
$t$ =: time, s
$U$ =: velocity (vector), $m / s$
$u$ =: velocity, $m / s$
$x$ =: quality

References

1.

Klein

,

S.

, and

Nellis

,

G.

,

2012

,

Thermodynamics

,

Cambridge University Press

,

Cambridge, UK

.

2.

Dobashi

,

K.

,

Kimura

,

A.

,

Oka

,

Y.

, and

Koshizuka

,

S.

,

1998

, “

Conceptual Design of a High Temperature Power Reactor Cooled, and Moderated by Supercritical Light Water

,”

Ann. Nucl. Energy

,

25

(

8

), pp.

487

–

505

.10.1016/S0306-4549(97)00079-0

Google Scholar

Crossref

3.

Oka

,

Y.

, and

Koshizuka

,

S.

,

2001

, “

Supercritical-Pressure, Once-Through Cycle Light Water Cooled Reactor Concept

,”

J. Nucl. Sci. Technol.

,

38

(

12

), pp.

1081

–

1089

.10.1080/18811248.2001.9715139

Google Scholar

Crossref

4.

Oka

,

Y.

,

Koshizuka

,

S.

, and

Yamasaki

,

T.

,

1992

, “

Direct Cycle Light Water Reactor Operating at Supercritical Pressure

,”

J. Nucl. Sci. Technol.

29

(

6

), pp.

585

–

588

.10.1080/18811248.1992.9731568

Google Scholar

Crossref

5.

Chatharaju

,

M.

,

2011

,

Computational Study of Critical Flow Discharge in Supercritical Water Cooled Reactors

,

McMaster University

,

Hamilton, Ontario, Canada

.

6.

Bishop

,

A.

,

Efferding

,

L.

, and

Tong

,

L.

,

1962

, “

A Review of Heat Transfer, and Fluid Flow of Water in the Supercritical Region, and During “Once-Through” Operation

,” ,

Westinghouse Electric Corporation, Atomic Power Division

, Pittsburgh, PA,

106

p.

7.

Dostal

,

V.

,

2004

,

A Supercritical Carbon Dioxide Cycle for Next Generation Nuclear Reactors

,

Massachusetts Institute of Technology

,

Cambridge, MA

.

8.

Dostal

,

V.

,

Hejzlar

,

P.

, and

Driscoll

,

M. J.

,

2006

, “

The Supercritical Carbon Dioxide Power Cycle: Comparison to Other Advanced Power Cycles

,”

Nucl. Technol.

,

154

(

3

), pp.

283

–

301

.

Google Scholar

Crossref

9.

Wright

,

S. A.

,

Radel

,

R. F.

,

Vernon

,

M. E.

,

Rochau

,

G. E.

, and

Pickard

,

P. S.

,

2010

,

Operation and Analysis of a Supercritical CO2 Brayton Cycle

,

Albuquerque

,

NM and Livermore, CA

.

10.

Wright

,

S. A.

,

Vernon

,

M. E.

, and

Pickard

,

P. S.

,

2006

,

Concept Design for a High Temperature Helium Brayton Cycle with Interstage Heating and Cooling

,

Albuquerque

,

NM and Livermore, CA

.

11.

Mignot

,

G. P.

,

Anderson

,

M. H.

, and

Corradini

,

M. L.

,

2009

, “

Measurement of Supercritical

{CO}_{2}

Critical Flow: Effects of L/D, and Surface Roughness

,”

Nucl. Eng. Des.

,

239

(

5

), pp.

949

–

955

.10.1016/j.nucengdes.2008.10.031

Google Scholar

Crossref

12.

Chen

,

J. P.

,

Liu

,

J. P.

,

Chen

,

Z. J.

, and

Niu

,

Y. M.

,

2004

, “

Trans-critical R744, and Two-Phase Flow Through Short Tube Orifices

,”

Int. J. Therm. Sci.

,

43

(

6

), pp.

623

–

630

.10.1016/j.ijthermalsci.2003.10.011

Google Scholar

Crossref

13.

Liu

,

J. P.

,

Niu

,

Y. M.

,

Chen

,

J. P.

,

Chen

,

Z. J.

, and

Feng

,

X.

,

2004

, “

Experimentation, and Correlation of R744 Two-Phase Flow Through Short Tubes

,”

Exp. Therm. Fluid Sci.

,

28

(

6

), pp.

565

–

573

.10.1016/j.expthermflusci.2003.08.004

Google Scholar

Crossref

14.

Chen

,

Y.

,

Zhao

,

M.

,

Yang

,

C.

,

Bi

,

K.

,

Du

,

K.

, and

Zhang

,

S.

,

2010

, “

Critical Flow of Water Under Supercritical Pressures

,”

2010 14th International Heat Transfer Conference

, Vol.

2

,

ASME

,

Washinton, DC

, pp.

319

–

326

.

15.

Muftuoglu

,

A.

, and

Teyssedou

,

A.

,

2014

, “

Experimental Study of Abrupt Discharge of Water at Supercritical Conditions

,”

Exp. Therm. Fluid Sci.

,

55

(

May

), pp.

12

–

20

.10.1016/j.expthermflusci.2014.02.009

Google Scholar

Crossref

16.

Rodarte

,

M. A.

,

2011

,

The Development of an Experimental Test Facility to Measure Leakage through Labyrinth Seals

,

University of Wisconsin-Madison

,

Madison, WI

.

17.

Edlebeck

,

J. P.

,

2013

,

Measurements and Modeling of the Flow of Supercritical Carbon Dioxide

, Master thesis,

University of Wisconsin-Madison

,

Madison, WI

.

18.

Wolf

,

M. P.

,

2014

,

Flow of Supercritical Carbon Dioxide Through Annuli and Labyrinth Seals

, Master thesis,

University of Wisconsin-Madison

,

Madison, WI

.

19.

Wallis

,

G. B.

,

1980

, “

Critical Two-Phase Flow

,”

Int. J. Multiphase Flow

,

6

(

1

), pp.

97

–

112

.

Google Scholar

Crossref

20.

Hesson

,

J. C.

, and

Peck

,

R. E.

,

1958

, “

Flow of Two-Phase Carbon Dioxide Through Orifices

,”

AIChE J.

,

4

(

2

), pp.

207

–

210

.10.1002/(ISSN)1547-5905

Google Scholar

Crossref

21.

Simoneau

,

R. J.

, and

Hendricks

,

R. C.

,

1977

, “

Generalized Charts for Computation of Two-Phase Choked Flow of Simple Cryogenic Liquids

,”

Cryogenics

,

17

(

2

), pp.

73

–

76

.10.1016/0011-2275(77)90099-6

Google Scholar

Crossref

22.

Zhang

,

C.

, and

Yang

,

L.

,

2005

, “

Modeling of Supercritical

{CO}_{2}

Flow Through Short Tube Orifices

,”

J. Fluids Eng.

,

127

(

6

), pp.

1194

–

1198

.10.1115/1.2060738

Google Scholar

Crossref

23.

Yang

,

L.

, and

Zhang

,

C.

,

2009

, “

Modified Neural Network Correlation of Refrigerant Mass Flow Rates Through Adiabatic Capillary, and Short Tubes: Extension to

{CO}_{2}

Transcritical Flow

,”

Int. J. Refrig.

,

32

(

6

), pp.

1293

–

1301

.10.1016/j.ijrefrig.2009.03.005

Google Scholar

Crossref

24.

Yadav

,

A. K.

,

Ram Gopal

,

M.

, and

Bhattacharyya

,

S.

,

2012

, “

Computational Fluid Dynamic Analysis of a Supercritical

{CO}_{2}

Based Natural Circulation Loop With End Heat Exchangers

,”

Int. J. Adv. Eng. Sci. Appl. Math.

,

4

(

3

), pp.

119

–

126

.10.1007/s12572-012-0062-2

Google Scholar

Crossref

25.

Lemmon

,

E. W.

,

Huber

,

M. L.

, and

Mclinden

,

M. O.

,

2007

,

NIST Reference Fluid Thermodynamic and Transport Properties-REFPROP

, Boulder, CO.

26.

Van Abel

,

E.

,

Anderson

,

M.

, and

Corradini

,

M.

,

2011

, “

Numerical Investigation of Pressure Drop and Local Heat Transfer of Supercritical

{CO}_{2}

in Printed Circuit Heat Exchangers

,” Supercritical

{CO}_{2}

Power Cycle Symposium,

Boulder, CO

.

27.

Suo-anttila

,

A. J.

, and

Wright

,

S. A.

,

2011

, “

Computational Fluid Dynamics Code for Supercritical Fluids

,” Supercritical

{CO}_{2}

Power Cycle Symposium,

Boulder, CO

.

28.

Qiu

,

L.

, and

Reitz

,

R. D.

,

2014

, “

Condensation Processes in a Motoring Engine

,”

J. Supercrit. Fluids

,

90

(

June

), pp.

84

–

100

.10.1016/j.supflu.2014.03.013

Google Scholar

Crossref

29.

Qiu

,

L.

,

Wang

,

Y.

, and

Reitz

,

R. D.

,

2014

, “

On Regular, and Retrograde Condensation in Multiphase Compressible Flows

,”

Int. J. Multiphase Flow

,

64

(

1

), pp.

85

–

96

.10.1016/j.ijmultiphaseflow.2014.05.004

Google Scholar

Crossref

30.

Fairweather

,

M.

,

Falle

,

S.

,

Hebrard

,

J.

,

Jamois

,

D.

,

Proust

,

C.

,

Wareing

,

C.

, and

Woolley

,

R.

,

2012

, “

Reynolds-Averaged Navier-Stokes Modelling of the Near-Field Structure of Accidental Releases of Carbon Dioxide from Pipelines

,”

Proceedings of the 22nd European Symposium on Computer Aided Process Engineering

,

London, UK

.

31.

OpenCFD Limited

,

2011

,

OpenFOAM User Guide

, [Online]. Available: http://www.openfoam.org/docs/user/.

32.

Van Doormaal

,

J. P.

, and

Raithby

,

G. D.

,

1984

, “

Enhancements of the Simple Method for Predicting Incompressible Fluid Flows

,”

Numer. Heat Transfer.

,

7

(

2

), pp.

147

–

163

.10.1080/01495728408961817

Google Scholar

Crossref

33.

Northland Numerics

,

2012

,

FIT

, [Online]. Available:http://www.northlandnumerics.com/.

34.

OpenFOAM Foundation

,

2011

,

Open FOAM Programmer’s Guide

, [Online]. Available: http://www.foamcfd.org/Nabla/guides/ProgrammersGuide.html.

35.

He

,

S.

,

Kim

,

W. S.

, and

Bae

,

J. H.

,

2008

, “

Assessment of Performance of Turbulence Models in Predicting Supercritical Pressure Heat Transfer in a Vertical Tube

,”

Int. J. Heat Mass Transf.

,

51

(

19–20

), pp.

4659

–

4675

.10.1016/j.ijheatmasstransfer.2007.12.028

Google Scholar

Crossref

36.

Yoo

,

J. Y.

,

2013

, “

The Turbulent Flows of Supercritical Fluids with Heat Transfer

,”

Annu. Rev. Fluid Mech.

,

45

(

1

), pp.

495

–

525

.10.1146/annurev-fluid-120710-101234

Google Scholar

Crossref

37.

Jones

,

W.

, and

Launder

,

B. E.

,

1972

, “

The Prediction of Laminarization with a Two-Equation Model of Turbulence

,”

Int. J. Heat Mass Transf.

,

15

(

2

), pp.

301

–

304

.10.1016/0017-9310(72)90076-2

Google Scholar

Crossref

38.

Menter

,

F. R.

,

1993

, “

Zonal Two Equation

k - ω

Turbulence Models for Aerodynamic Flows

,”

24th Fluid Dynamics Conference

,

Orlando, Florida

.

39.

Yoo

,

J. Y.

,

2013

, “

The Turbulent Flows of Supercritical Fluids with Heat Transfer

,”

Annu. Rev. Fluid Mech.

,

45

(

1

), pp.

495

–

525

.

40.

Bae

,

Y.

,

Hong

,

S.

, and

Kim

,

Y.

,

2012

, “

Numerical Simulation of Supercritical Heat Transfer Under Severe Axial Density Gradient in a Narrow Vertical Tube

,”

The 12th International Congress on Advances in Nuclear Power Plants

,

Chicago,IL

.

41.

Jaromin

,

M.

, and

Anglart

,

H.

,

2013

, “

A Numerical Study of the Turbulent Prandtl Number Impact on Heat Transfer to Supsercritical Water Flowing Upward Under Deteriorated Conditions

,”

The 15th International Topical Meeting on Nuclear Reactor Thermalhydraulics

,

Pisa, Italy

.

42.

Weisman

,

J.

,

Duncan

,

D.

,

Gibson

,

J.

, and

Crawford

,

T.

,

1979

, “

Diameter on Two-Phase Flow Patterns in Horizontal Lines

,”

Int. J. Multiphase Flow

,

5

(

1

), pp.

437

–

462

.

Google Scholar

Crossref

43.

Weisman

,

J.

, and

Kang

,

S.

,

1981

, “

Flow Pattern Transitions in Vertical, and Upwardly Inclined Lines

,”

Int. J. Multiphase Flow

,

7

(

1

), pp.

271

–

291

.10.1016/0301-9322(81)90022-7

Google Scholar

Crossref

44.

Moody

,

F. J.

,

1965

, “

Maximum Flow Rate of a Single Component, Two-Phase Mixture

,”

J. Heat Transfer

,

87

(

1

), pp.

134

–

141

.10.1115/1.3689029

Google Scholar

Crossref

45.

Fauske

,

H. K.

,

1961

,

Contribution to the Theory of Two-Phase, One-Component Critical Flow

,

Chicago, IL

.

Simulation of Supercritical CO₂ Flow Through Circular and Annular Orifice

Introduction and Motivation

Previous Work

Numerical Methodology

Solver.

Modification of Energy Equation.

Real Fluid Model.

Homogeneous Equilibrium Model.

Geometry and Boundary Conditions.

Meshing Sensitivity.

Turbulence Model.

Test Facility

Results and Discussion

Circular Orifice.

Annular Orifice.

Discussion.

Conclusion and Future Work

Acknowledgment

Nomenclature

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

Simulation of Supercritical CO2 Flow Through Circular and Annular Orifice

Introduction and Motivation

Previous Work

Numerical Methodology

Solver.

Modification of Energy Equation.

Real Fluid Model.

Homogeneous Equilibrium Model.

Geometry and Boundary Conditions.

Meshing Sensitivity.

Turbulence Model.

Test Facility

Results and Discussion

Circular Orifice.

Annular Orifice.

Discussion.

Conclusion and Future Work

Acknowledgment

Nomenclature

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only

Simulation of Supercritical CO₂ Flow Through Circular and Annular Orifice