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 (45%) 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.

The continuity equation is  
·(ρ·U)=0
(1)
The momentum equation is  
·(ρ·U·U)·(μeffective·U)=P
(2)
and the energy equation is  
·(ρ·U·h)+·(ρ·U·U22)·(ρ·αeffective·h)=0
(3)
where  
μeffective=μ+μt
(4)
 
αeffective=α+αt
(5)
 
Prt=μtρ·α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 CO2 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 (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 α, which is not defined in the two-phase region  
·(ρ·U·h)+·(ρ·U·U22)·(ρ·αeffective·h)=0
(7)
 
·(ρ·U·h)+·(ρ·U·U22)·(ρ·αt·h)=·(κ·T)
(8)

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.

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. (8) and (9). The energy equation in Eq. (7) solves enthalpy. Then, quality is determined from the enthalpy value using Eq. (9), which is then used to determine other properties  
vHEM=x·vg+(1x)·vl
(9)
 
hHEM=x·hg+(1x)·hl
(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. (10) and (11). As the high Re (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·μg+(1x)·μl
(11)
 
κHEM=x·κg+(1x)·κ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.

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 475kg/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.

In standard k-epsilon model  
DρkDt=·((μtσk+μ)k)+ρPk23ρk·Uρε
(13)
 
DρεDt=·((μtσε+μ)ε)+Cε1ρPkεk(23Cε1+Cε3)ρε·UCε2ρε2k
(14)
 
μt=Cμρk2ε
(15)
 
αt=μtρ·Prt
(16)
with coefficients provided in Table 1.
In the k-omega SST model  
DρkDt=·((αkμt+μ)k)+ρPk23ρk·Uρβωk
(17)
 
DρωDt=·((αωμt+μ)ω)+γρPkωk23ργω·Uρβω2ρ(F1)F
(18)
 
μt=ρkω
(19)
 
αt=μtρ·Prt
(20)
where  
F=2αω2(k·ω)ω,αk=Fαk1+(1F)αk2,αω=Fαω1+(1F)αω2,γ=Fγ1+(1F)γ2,β=Fβ1+(1F)β2
with coefficients provided in Table 2.

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 475kg/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, 498kg/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 Ts diagram. All of the inlet conditions correspond to a supercritical state; one of the inlet conditions (pressure is 7.7 MPa, density is 111kg/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, 498kg/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, 498kg/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.

In the flow regime map from Weisman [42,43], as shown in Fig. 19, when the velocities for both liquid and gas are high (50100m/s), a dispersed flow pattern is the result. In the following, we show that the velocities for both liquid and gas are very high (50100m/s), and a dispersed flow pattern is suitable to describe the flow. The ideal nozzle isentropic model (Eqs. (20)–(23)) combined with Moody’s slip ratio model (Eq. (24)) is used to prove this argument qualitatively. Void fraction is calculated based on Moody’s slip ratio model using Eq. (25). The total mass flux should be consistent by Eq. (26)  
hinlet=enthalpy(Pinlet,sinlet)
(21)
 
houtlet=enthalpy(Poutlet,sinlet)
(22)
 
ρoutlet=density(Poutlet,sinlet)
(23)
 
hinlet=houtlet+uoutlet2
(24)
 
S=(ρlρg)1/3=ugul
(25)
 
α=11+(1xx)(ρgρl)S
(26)
 
αρgug+(1α)ρlul=ρu
(27)

Let us consider an example case: sCO2 discharges from 9 MPa and 338kg/m3 to 6 MPa, with calculated velocities for the two phases that are ug=104m/s and ul=68m/s. With the correction factors provided by Weisman [42,43], the corrected velocities are ug/ϕ1=2205m/s and ul/ϕ2=68m/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 (25) is used to calculate the void fraction, Eq. (27) represents Fauske’s slip ratio model, and Eq. (28) is applied to calculate the averaged density of two phases  
S=(ρlρg)1/2=ugul
(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.

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 498kg/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. 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. 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

    Latin Symbols
     
  • Cp =

    specific heat, J/(kg·K)

  •  
  • h =

    specific enthalpy, J/kg

  •  
  • k =

    turbulent kinetic energy, J/kg

  •  
  • m˙ =

    mass flow rate, kg/s

  •  
  • P =

    pressure, Pa

  •  
  • Pk =

    turbulent kinetic energy production rate, J/(kg·s)

  •  
  • 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·K)

  •  
  • T =

    temperature, °C

  •  
  • t =

    time, s

  •  
  • U =

    velocity (vector), m/s

  •  
  • u =

    velocity, m/s

  •  
  • x =

    quality

    • Greek Symbols
       
    • α =

      thermal diffusivity, m2/s, or void fraction

    •  
    • αt =

      turbulence diffusivity, m2/s

    •  
    • ε =

      turbulent kinetic energy dissipation rate, J/(kgs)

    •  
    • κ =

      thermal conductivity, W/(mK)

    •  
    • υ =

      specific volume, m3/kg

    •  
    • μ =

      dynamic viscosity, Pa s

    •  
    • μt =

      turbulence viscosity, Pa s

    •  
    • ρ =

      density, kg/m3

    •  
    • ϕ1 =

      Weisman correction factor for gas phase

    •  
    • ϕ2 =

      Weisman correction factor for liquid phase

    •  
    • ω =

      turbulence frequency, 1/s

      • Subscripts
         
      • 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

        • Abbreviations
           
        • CFD =

          computational fluid dynamic

        •  
        • CO2 =

          carbon dioxide

        •  
        • DC =

          discharge coefficient

        •  
        • EXP =

          experiment

        •  
        • HEM =

          homogeneous equilibrium model

        •  
        • IS =

          isentropic model

        •  
        • OF =

          OpenFOAM

        •  
        • RANS =

          Reynolds-averaged Navier–Stokes equation

        •  
        • sCO2 =

          supercritical carbon dioxide

References

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
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
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
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
.
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 CO2 Critical Flow: Effects of L/D, and Surface Roughness
,”
Nucl. Eng. Des.
,
239
(
5
), pp.
949
955
.10.1016/j.nucengdes.2008.10.031
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
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
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
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
.
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
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
22.
Zhang
,
C.
, and
Yang
,
L.
,
2005
, “
Modeling of Supercritical CO2 Flow Through Short Tube Orifices
,”
J. Fluids Eng.
,
127
(
6
), pp.
1194
1198
.10.1115/1.2060738
23.
Yang
,
L.
, and
Zhang
,
C.
,
2009
, “
Modified Neural Network Correlation of Refrigerant Mass Flow Rates Through Adiabatic Capillary, and Short Tubes: Extension to CO2 Transcritical Flow
,”
Int. J. Refrig.
,
32
(
6
), pp.
1293
1301
.10.1016/j.ijrefrig.2009.03.005
24.
Yadav
,
A. K.
,
Ram Gopal
,
M.
, and
Bhattacharyya
,
S.
,
2012
, “
Computational Fluid Dynamic Analysis of a Supercritical CO2 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
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 CO2 in Printed Circuit Heat Exchangers
,” Supercritical CO2 Power Cycle Symposium,
Boulder, CO
.
27.
Suo-anttila
,
A. J.
, and
Wright
,
S. A.
,
2011
, “
Computational Fluid Dynamics Code for Supercritical Fluids
,” Supercritical CO2 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
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
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
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
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
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
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
.
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
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
45.
Fauske
,
H. K.
,
1961
,
Contribution to the Theory of Two-Phase, One-Component Critical Flow
,
Chicago, IL
.