Study of Supercritical CO₂ Pipeline Flow Leaks: Effect of Equation of State, Impurity, and Outlet Diameter

1. Introduction

The past decade has witnessed a significant expansion of $C{O}_2$ pipeline networks globally, driven by the growing emphasis on carbon capture, utilization, and storage (CCUS) initiatives, application in enhanced oil recovery (EOR), and to reduce anthropogenic $C{O}_2$ emissions and meet global de-carbonization targets ([1,2]). Recently, significant progress has been made in CCUS technologies ([3,4,5]) and it has been found that the systems perform efficiently when the fluid ($C{O}_2$) is in supercritical state ($sC{O}_2$). The accurate modeling of supercritical carbon-dioxide $sC{O}_2$ transport thus becomes increasingly important for the design, safety assessment, cost reduction, and operational reliability of pipeline and storage systems in CCUS ([6,7,8]). One of the critical question in the safety assessment of $sC{O}_2$ transport systems is the mass flow rate during depressurization or leakage, which is strongly influenced by decompression wave speed and mixture impurities. Non-equilibrium phase change retarded by nucleation and vapor-front dynamics is found to have significant effects, but only during the leakage early stage (few milliseconds) [9,10,11]. Consequently, metastability phenomena (i.e. non-equilibrium thermodynamics) are neglected in this article, as we intend to perform CFD simulations over significantly longer periods of time. Early work by Mahgerefteh et al. ([12]) employed 1-D models to investigate the influence of friction, heat transfer, and stream impurities on decompression wave speeds in $sC{O}_2$ pipelines. Subsequent studies by Woolley et al. ([6,13,14]) combined experimental and numerical approaches to analyze multi-component $C{O}_2$ mixtures, providing foundational data for safety guidelines. The recent reviews [15,16,17,18] discuss the works and challenges in accurate modeling of the decompression wave, the thermophysical properties, phase equilibria, and dispersion of $C{O}_2$ and its mixtures. Shuaiwei et al. [19] showed that isentropic assumption during depressurization predicts the decompression wave speed for $C{O}_2$ and its mixtures quite accurately. Similarly, Magen et al. [20] assuming isentropic phase change provided formulation for predicting onset of nucleation by estimating the minimum depressurization rate near the spinodal limit. Xiaoqiang et al. [21] using their experimental data correlated non-dimensional depressurization rate and degree of superheat. Liao et al. [22] analyzed the critical mass flux and decompression wave speed for $C{O}_2$ and its mixtures through 1-D analysis using different models and provided a correlation for the mass flow rate. Cao et al. [23] focused on the flow field on cross section of pipes and describe the stratification of phases during pressure leak. Hansen et al. [24,25] analyzed the problem on a vertical pipe duct and estimated reference data for 1-D phase transition model. Bhuvanker et al. [26] studied the effect of heat transfer and gravity in $C{O}_2$ blowout. Yin et al. [27] studied the depressurization phenomenon numerically using a homogeneous equilibrium model (HEM) and a homogeneous relaxation model (HRM). Log et al. [28] and Wang et al. [29] proposed correlations for the HRM relaxation time based on initial entropy. Log et al. [30] proposed a model to include the effect of the nucleation and bubble growth into the source for mass transfer. Munkejord et al. [31] and Yu et al. [32] studied the effect of impurities for $C{O}_2-N_2$ mixture. Zhu et al. [33] analyzed the effect of water hammering on phase transition for the pipeline transport. Lio et al [34] proposed a neural network model to predict critical mass flow. Yu et al [35] studied the depressurization phenomenon during multistage throttling. While the previous research discussed were conducted for pressure puncture, many research were also done on orifices leakage. Ding et al. [36] proposed a pressure correction for estimating corrected leakage rate. Zhang et al. [37] analyzed the near field characteristics for $sC{O}_2$ during pressure leak and proposed an empirical correlation of saturation pressure $P_{sat}$ for numerical modeling. Yu et al. [38] studied the temperature distribution and $C{O}_2$ dispersion for a pipe leak buried in soil. Similarly, in [39] they analyzed the problem numerically using leakage and seepage diffusion model. Chen et. al. [40] analyzed the effect of orifice diameter, similarly Hu et al. [41] analyzed the effect of orifice diameter for underwater pipeline leakage.

The numerical modeling of operations close to the supercritical state poses a number of challenges. Firstly, the non-linear behavior of thermodynamic and transport properties ([17,32]) near the critical point. Further complicating factor is the presence of impurities such as $H_{2}O,N_2,H_2,C{H}_4,O_2,H_{2}S,S{O}_x$, and $N{O}_x$, which significantly alter mixture thermodynamics and transport properties. Pure $C{O}_2$ can be accurately modeled using multi-parameter equations of state (EoS) like Span and Wagner ([42]), but it is computationally expensive. Nazeri et al. ([43]) demonstrated that incorporating a volume correction term into the Soave-Redlich-Kwong (SRK) EoS enhances density predictions for $C{O}_2$-rich mixtures. Similarly, Demetriades et al. ([44]) derived a pressure-explicit EoS for $C{O}_2-N_2-O_2-H_2$ systems, which outperformed the GERG-2008 model ([45]) in certain regimes. Petropoulou et al. ([46]) further compared various EoS with Universal Mixing Rules (UMR) for $C{O}_2/C{H}_4$ mixtures, finding that the Peng-Robinson (PR) UMR EoS provided the best vapor-liquid equilibrium (VLE) predictions. Similarly, Li and Yan [47] analyzed five different cubic EoS for VLE calculations of $C{O}_2$ and $C{O}_2$ mixtures. Diamantonis et al. [48] compared SAFT and PC-SAFT [49] with cubic EoS for VLE modeling of $C{O}_2$ mixtures. Avendano et al. [50] has presented a parameter set for the SAFT-$\gamma$ model [51] for $C{O}_2$. Accurate prediction of transport properties (e.g., viscosity, thermal conductivity) remains another key hurdle. While empirical correlations like those of Bahadori and Vuthaluru ([52]) and Nazeri et al. ([53]) have improved viscosity estimates for $C{O}_2$-rich mixtures, thermal conductivity modeling relies heavily on correlations such as those by Jarrahian and Heidaryan ([54]) and Rostami et al. ([55]). Cross-property interactions, such as those studied by Hellmann et al. ([56]) using the classical trajectory method, further complicate the picture. A persistent challenge is the scarcity of high-quality experimental data ([57]), which limits the validation and refinement of these models. The nonlinearities introduced by real-fluid EoS, and phase change phenomenon increases the numerical stiffness of the system as well. This severely affects the stability and convergence of solvers, consequently raising the computational cost. Modeling of formation of dry ice, which is very common for the case of $sC{O}_2$ depressurization is another challenge. An accurate EoS model of solid $C{O}_2$ is needed. Trusler ([58,59]) and Jager and Span ([60]) proposed a Helmholtz and Gibbs energy based model for solid $C{O}_2$, respectively. These models have been validated for various systems ([61,62,63]) and found to be reasonably accurate with experiment. Maltby ([64]) reviewed various experimental data for various $C{O}_2$ systems and found out the deviation is minimum with PR EoS considering overall range of various conditions. Bhatia et al. ([65]) found that CPA outperforms PR for the case of condensation of $C{O}_2$ through convergent-divergent nozzle.

Despite the substantial body of work summarized above, several important limitations remain. Most existing investigations rely predominantly on one–dimensional homogeneous equilibrium or relaxation models to predict decompression wave speeds and critical mass flux. While these approaches provide valuable theoretical insight and are computationally efficient, they inherently neglect multidimensional flow structures, spatial phase stratification, and local non-uniformity that arise during realistic pipeline leakage. Furthermore, many studies emphasize either thermodynamic modeling or empirical correlations independently, with limited integration into fully coupled CFD frameworks capable of resolving transient flow–thermodynamics interactions over engineering time scales. In addition, although the influence of impurities has been recognized, a systematic numerical investigations incorporating real-fluid equations of state together with transport-property variations remain scarce. The combined effects of non-linear thermophysical behavior near the critical region, phase change, and stiff source terms pose significant challenges for stable and efficient numerical solvers, which has limited the application of high-fidelity CFD to practical large-scale problems. Consequently, there is still a lack of robust, computationally tractable models that can accurately predict depressurization dynamics and mass discharge rates for realistic multi-component $C{O}_2$ mixtures over extended transient periods relevant to safety assessment and pipeline design. The present work aims to address these gaps by developing and applying a CFD-based framework that incorporates real-fluid thermodynamics, mixture-dependent properties, and transient multiphase effects to systematically evaluate decompression behavior and leakage characteristics of $C{O}_2$ pipelines under practical operating conditions. In the present work, a high-fidelity Real Fluid Method (RFM) solver implemented on CONVERGE CFD v4.1 is used. The framework is validated for real-gas evaporating and non-evaporating sprays [66,67,68]. The thermophysical properties are calculated by a robust thermodynamic solver developed at IFP Energies Nouvelle [69]. Model validation leverages high-resolution experimental pressure–temperature data from Munkejord et al. [31,70,71], providing a rigorous benchmark for assessing accuracy. The selected binary systems are $C{O}_2-N_2$, $C{O}_2-C{H}_4$, and $C{O}_2-Ar$ representing the major non-condensable impurities relevant to current and emerging $sC{O}_2$ capture, transport, and storage applications [72]. The key contributions of this work are:

1.

A critical assessment of the influence of the equation of state on the prediction of thermophysical properties for supercritical $C{O}_2$ and its mixtures,

2.

$C{O}_2$ rich mixtures depressurization analysis of the flow rate, and pressure depression and void fraction profiles with different outlet diameters and geometries,

3.

An evaluation of the effect of impurities by simulating binary mixture on the decompression wave propagation and the resulting critical mass flow rate.

The following sections detail the numerical methodology (Section 2) and computational setup and mesh sensitivity study (Section 3). Results and discussions (Section 4) include the validation of CPA EoS and VLE thermodynamic equilibrium. Effect of impurities is studied in detail for various parameters of concern in pipeline industry. Section (Section 4) also compare simulated pressure profile, gas fraction variations, and thermodynamic path against experimental benchmarks. Effect of impurities on decompression wave speed, mass flow rate and Mach number are also discussed. Finally, the conclusion discusses key findings and highlights the broader implications of impurity-sensitive modeling for the design and transient safety assessment of industrial-scale $sC{O}_2$ transport systems.

2. Numerical Methodology

2.1. Governing Equations

The mixture is represented as a single effective fluid in which both phases share identical velocity, pressure, and temperature fields. Accordingly, the formulation reduces to four governing equations [73,74,75,76]:

$${\displaystyle \frac{\partial \rho }{\partial t}}+{\displaystyle \frac{\partial \rho u_i}{\partial x_i}}=0$$

(1)

$${\displaystyle \frac{\partial \rho u_i}{\partial t}}+{\displaystyle \frac{\partial \rho u_i{u}_j}{\partial x_j}}=-{\displaystyle \frac{\partial P}{\partial x_i}}+{\displaystyle \frac{\partial }{\partial x_j}}\left({\tau }_{ij}+{\tau }_{ij}^{SGS}\right)$$

(2)

$${\displaystyle \frac{\partial \rho e}{\partial t}}+{\displaystyle \frac{\partial \rho e{u}_j}{\partial x_j}}=-P{\displaystyle \frac{\partial u_j}{\partial x_j}}+\left({\tau }_{ij}+{\tau }_{ij}^{SGS}\right){\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial }{\partial x_j}}\left(q_j+q_j^{SGS}\right)$$

(3)

$${\displaystyle \frac{\partial \rho Y_k}{\partial t}}+{\displaystyle \frac{\partial \rho Y_k{u}_j}{\partial x_j}}=-{\displaystyle \frac{\partial }{\partial x_j}}\left(J_{kj}+J_{kj}^{SGS}\right)$$

(4)

Here, $\rho$, u, P, T, and e denote the mixture density, velocity, pressure, temperature, and specific internal energy. The molecular viscous stress tensor is given by

$${\tau }_{ij}=\mu \left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}-{\displaystyle \frac{2}{3}}{\displaystyle \frac{\partial u_k}{\partial x_k}}{\delta }_{ij}\right)$$

(5)

with ${\delta }_{ij}$ the Kronecker delta.

The heat flux and diffusive flux of species are expressed as

$$q_j=-\lambda {\displaystyle \frac{\partial T}{\partial x_j}}+\sum _m{J}_{mj}{h}_m,J_{mj}=-\rho D_m{\displaystyle \frac{\partial Y_m}{\partial x_j}}$$

(6)

where $h_m$ is the species enthalpy, $D_m$ the molecular diffusivity, and $Y_m$ the mass fraction of species m. The transport coefficients $\mu$ and $\lambda$ are obtained from the correlations of Chung et al. [77].

The sub-grid scale (SGS) stresses and fluxes are modeled by analogy with molecular terms. The SGS stress tensor is

$${\tau }_{ij}^{SGS}={\displaystyle \frac{{\mu }^{SGS}}{\mu }}{\tau }_{ij}$$

(7)

while the turbulent contributions to heat and mass transfer are

$$q_j^{SGS}=-{\lambda }_t{\displaystyle \frac{\partial T}{\partial x_j}}+\sum _m{J}_{mj}^{SGS}{h}_m,J_{mj}^{SGS}=-\rho D_{mt}{\displaystyle \frac{\partial Y_m}{\partial x_j}}$$

(8)

The turbulent transport properties are modeled through

$$D_{mt}={\displaystyle \frac{{\mu }_{SGS}}{\rho S{c}_t}},{\lambda }_t={\displaystyle \frac{C_p{\mu }_{SGS}}{P{r}_t}}$$

(9)

with $S{c}_t=0.7$, $P{r}_t=0.9$, and $C_p$ the specific heat at constant pressure. The eddy viscosity ${\mu }_{SGS}$ follows the $\sigma$-model [78],

$${\mu }_{SGS}=\rho {(C\Delta )}^2\delta$$

(10)

where $C=1.5$, $\Delta =V^{1/3}$, V being the grid-cell volume, and $\delta ={\displaystyle \frac{{\sigma }_3({\sigma }_1-{\sigma }_2)({\sigma }_2-{\sigma }_3)}{{\sigma }_1^2}}$ using the singular values ${\sigma }_i$ of the resolved velocity-gradient tensor. In this work, cross-coupling effects such as Soret and Dufour are neglected.

2.2. Thermodynamic Modeling

The Real-Fluid Method (RFM) [73,74,75,76] is based on tabulated thermodynamic data generated using the IFPEN-Carnot library [69]. The library applies PT-flash calculations for mixtures and PH-flash for pure fluids, with different cubic and association EoS (PR [79], CPA [80], PC-SAFT [49]). Flashing involves a stability test (tangent-plane distance criterion) followed by phase-split determination.

The lookup tables, discretized in $(T,P,Y_k)$ for mixtures (or $(h,P)$ for pure fluids), provide equilibrium density, internal energy, vapor fraction, heat capacity, sound speed, and transport properties. Property queries are interpolated via Inverse Distance Weighting (IDW) [81]. This methodology has been successfully applied to binary [66,68,82] and ternary [67] mixtures and is implemented in CONVERGE CFD [83].

The tables are used during the runtime for:

• Property evaluation as function $f(T,P,Y_k)$,
• Reverse lookup $T=f(e,P,Y_k)$ after solving for internal energy.

Because of the non-linear behavior of real-fluid EoS, the standard pressure equation must be reformulated. A pressure-density ratio is defined locally at each cell and time step to incorporate the effects of real fluid compressibility as

$${\varphi }^{*}={\left({\displaystyle \frac{\partial \rho }{\partial p}}\right)}_{T,Y}$$

(11)

Within the PISO loop, the pressure equation reads:

$${\displaystyle \frac{{\partial }^2(P^{**}-P^{*})}{\partial x_i\partial x_i}}-{\displaystyle \frac{{\varphi }^{*}(P^{**}-P^{*})}{d{t}^2}}={\displaystyle \frac{({\rho }^{Eo{S}^{*}}-{\rho }^n)}{d{t}^2}}+{\displaystyle \frac{1}{dt}}\left({\displaystyle \frac{\partial {\rho }^{*}{u}_i^{**}}{\partial x_i}}-S\right)$$

(12)

Here ** and * refer to successive PISO iterations, $EoS$ designate the density obtained from the table, and $n$ indicates the current time step. Although no explicit phase-change source terms are included in the transport equations under equilibrium, mass and energy transfer are consistently handled via property lookup and reverse-lookup procedures. A similar procedure is used for single component $C{O}_2$ simulations where the table inputs are ($h,P)$. In this case, h is evaluated as ($e+P/\rho$) and the temperature is among the output in the table.

2.3. Numerical Schemes

The spatial discretization is first-order accurate with a first-order Euler upwind scheme. Time integration is performed using a first-order explicit Euler forward scheme. This treatment enhances numerical robustness during the highly transient depressurization phase while maintaining computational efficiency. A dynamically adaptive time-stepping strategy is adopted. During the initial stage of depressurization, where the flow is strongly shock-dominated and characterized by steep pressure gradients, the time step is on the order of $\mathcal{O}\left(10\mathrm{ns}\right)$. As the flow transitions beyond the choked regime and pressure gradients relax, the time step is gradually increased up to $\mathcal{O}\left(10\mathrm{ms}\right)$. Numerical stability is primarily governed by shock propagation and rapid pressure-wave dynamics. The maximum allowable time step is controlled using a Courant–Friedrichs–Lewy (CFL) criterion based on the local Mach number. The maximum CFL number attained at maximum $dt$ is 2.0 for laboratory-scale pipeline simulations and 8.0 for industrial-scale pipeline configurations. The larger permissible CFL number in industrial-scale cases is attributed to the longer characteristic length scales and comparatively smoother spatial gradients following the initial decompression front propagation. All simulations are performed on AMD EPYC Genoa 9534 processors operating at 2.45 GHz on the ORION super-computing cluster. The number of cores used are 128 cores for lab scale pipeline flow and 256 cores for the industrial scale pipeline flow. Thermodynamic closure within the RFM framework is provided by the VTPR equation of state. Thermodynamic properties are evaluated using pre-tabulated data with a resolution of

$$\Delta T=1\mathrm{K},\Delta P=0.1\mathrm{bar},\Delta Y=\mathrm{Impurity}\mathrm{mass}\mathrm{fraction}.$$

(13)

The selected tabulation resolution is consistent with prior validation studies ([65,68,74,75,76]), ensuring an appropriate balance between interpolation accuracy and computational cost for multi-component CO₂ mixture depressurization simulations.

3. Computational Setup

3.1. Setup Description

The computational domain is shown in Figure 1. It is a cylindrical pipe of diameter $D=40.8$ mm and length $L=61.7$ m, relevant to the experimental setup of Munkejord et al. [31,70]. The domain has two boundaries. First, at left side is the OUTLET plane which is at the location of the disk rupture (simulating the leak), and the rest of the boundary are defined as WALL. Hence, the boundary condition are Dirichlet condition for pressure and Neumann for the rest of the depending variables such as temperature and Velocity at outlet. The outlet pressure is defined by the atmospheric condition with total pressure as 1 atm. At the Wall, no slip condition is imposed. The temperature is assumed to remain constant at the wall throughout the simulation i.e. $T_{wall}=T_0=298K$. The pipe is initially filled with supercritical $C{O}_2$, as the initial condition inside the pipeline is defined by the fluid thermodynamic state ($T_0=298K)$ and $(P_0=12.27MPa$) which is greater than the critical pressure ($P_c=7.3773MPa$) of pure $C{O}_2$. The disk ruptures at time $t=0$ and $sC{O}_2$ inside is exposed to atmospheric condition ($P_{amb}=0.1MPa,T_{amb}=277K$) causing a two-phase boiling shock including a sudden expansion wave with a violent phase change. The initial and atmospheric conditions remains the same for all the cases. The test matrix is designed to systematically assess the influence of the equation of state (EoS), impurities, outlet geometry and diameter. The test matrices for the present study are shown in Table 1, Table 2 and Table 3. In addition to these test cases, a reference case C0 of pure $C{O}_2$ is also studied, the initial and atmospheric conditions are similar i.e., $(P_0=12.22MPa,T_0=297.6K)$ and $(P_{amb}=0.1MPa,T_{amb}=281K)$ [70]. The geometry for the cases C10-C13 are taken from the reference [71] and are shown in Figure 1(b). The new outlet plane for these cases are at a distance $l=50mm$ from the previous outlet plane and the boundary area is varied through d, the diameter of the orifice or nozzle. For the case of orifice a rectangular boundary is chosen instead of circular keeping the geometry and outlet boundary net area identical. Fluid impurities and initial conditions are selected to provide detailed informations relevant to IFPEN’s current projects.

3.2. Mesh Sensitivity

The mesh sensitivity is performed for the Case C1 in Table 1. A uniform mesh is used for all the cases in the present work, with fixed 4 cells in the cross-section, i.e. ${\Delta }_x={\Delta }_y=20.4$ mm, and refined in axial direction making the setup closest to 1D-setup. The mesh sensitivity is analyzed with three different mesh resolutions i.e. $L/{\Delta }_z=4000,8000,16000$ which corresponds to ${\Delta }_z\approx 16,8,4$ mm and approximately $14.8,29.6,59.2$ thousands of cells. The results are analyzed through a monitor point at the distance $x=80$ mm, as shown in Figure 1, which corresponds to the normalized distance $x/D\approx 2$. The results from grid sensitivity are depicted in Figure 2. Apart from grid sensitivity the result also validates the numerical model with experimental data. The error percentage in the pressure and gas volume fraction profile are less than 10% before the fluctuation in properties begins due to interactions with the reflected pressure waves. The computational time taken for the different mesh resolutions i.e. $L/{\Delta }_z=4000,8000,16000$ are $3.24,6.29,20.7$ hrs, respectively. Hence, the mesh resolution with $L/{\Delta }_z=8000$ is found to be adequate for the present study.

4. Results and Discussions

Understanding the depressurization of supercritical $sC{O}_2$ requires capturing the interplay of thermodynamic and flow phenomena. The results are analyzed through pressure, temperature, and gas volume fraction data, the three independent thermodynamic property essential to define the thermodynamic state of the binary system. These data are also used to analyze the thermodynamic path of the system during depressurization. The decompression wave speed ($u_w$) represents the propagation velocity of a given pressure level P along the pipeline during depressurization. It is determined by tracking the arrival time of a specific pressure value at two axial locations. The two monitor points of choice are at $x_1=1.6$ m and $x_2=61.28$ m, giving a separation distance of $\delta x=x_2-x_1=59.68$ m. For a selected pressure level P, the corresponding arrival times $t_1\left(P\right)$ and $t_2\left(P\right)$ are identified from the pressure–time histories at $x_1$ and $x_2$, respectively. These times are defined as the instants at which the local pressure signal first reaches the given value P during depressurization. Linear interpolation between discrete data points is used to determine the arrival times. The decompression wave speed ($u_w$) is then computed as

$$u_w\left(P\right)={\displaystyle \frac{\delta x}{t_2\left(P\right)-t_1\left(P\right)}}.$$

(14)

Thus, $u_w\left(P\right)$ represents the propagation speed of an isobaric pressure level P along the pipeline. Apart from these, mass flow rate and Mach number are also studied as these informations are crucial for the design analysis of $sC{O}_2$ dispersion modeling. Mass flow rate and Mach number are obtained at the outflow boundary.

4.1. Heat Transfer and Boundary Conditions

The role of heat transfer models is significant in $sC{O}_2$ pipeline transport. The present section compares two boundary conditions i.e. adiabatic and isothermal conditions at the wall. For isothermal, the wall is assumed to be at the initial temperature i.e. $T_{wall}=T_0=298K$, throughout the observational time, i.e. $7s$. From Figure 3(a),(b) and (d), the system is found to follow the same thermodynamic path, pressure and temperature profile in the initial phase (up to $1.5s$) respectively, where the flow features are dominated by the initial shock wave expansion, phase change and Joule-Thompson cooling effect. Later on for the two-phase state fluid and especially when the fluid becomes gaseous, the heat transfer model starts showing major differences, as depicted in Figure 3(a) and Figure 3(c). It can be seen that a simplistic model of heat transfer that maintains the isothermal condition brought the depressurization pressure (Figure 3(b)) and temperature (Figure 3(c)) profiles close to experimental data. At adiabatic condition there is no source of heat intake into the fluid, hence the fluid pressure rapidly decays to atmospheric pressure and stays at minimum temperature obtained through Joule-Thompson cooling. Although, isothermal boundary condition is a better choice than adiabatic wall, this study shows that wall heat transfer needs to be computed accurately to capture the transition of the decompression wave speed as shown in Figure 3(c) for the various phases and the minimum temperature attained by the system as shown in Figure 3(d). However, the isothermal condition is found to follow well the behavior of the experimental data. Hence, the isothermal wall is the chosen boundary condition for further study in this work.

4.2. Effect of Equation of State

The cases C1–C3 in Table 1 are investigated to assess the effect of the equation of state (EoS) on the depressurization behavior of supercritical $sC{O}_2$. Three widely used EoS, namely CPA, Peng–Robinson (PR), and PC-SAFT, are considered in this study. These EoS yield slight variations in prediction of density, speed of sound, and phase equilibrium under near-critical and two-phase conditions. As a result, they influence key flow characteristics during depressurization, including pressure wave propagation, mass discharge rate, temperature evolution, and the onset and extent of phase change. The shape of the phase envelopes as computed from the three equations is shown in Figure 4. CPA and PC-SAFT show a larger critical temperature and pressure for the case of $C{O}_2+N_2(1.8\%mol)$. This also shifts the initiation point of vaporization at the bubble line (Figure 4(a)), resulting in the following order for saturation pressure: $Psa{t}_{PR}<Psa{t}_{CPA}<Psa{t}_{PC-SAFT}$ at a given temperature, while for the saturation temperature the order is: $Tsa{t}_{PR}>Tsa{t}_{CPA}>Tsa{t}_{PC-SAFT}$ at a given pressure. Because of different bubble lines of the EoS , the flow enters the two-phase flow state somewhat later for PR in comparison to CPA and PC-SAFT (see Figure 4(a)). Hence, a delayed chocked flow pressure and duration is observed for PR in comparison to CPA and PC-SAFT, as can seen at about $1s$ in the pressure and temperature evolution plotted in Figure 4(b,d). For the current monitor point, the flow enters to the gaseous phase at about $2s$, after two interactions with reflected pressure waves from the closed end of the pipe, which have led to more or less pressure oscillations in the two-phase flow region. The depressurization rate is different at each stage of interaction, and so is the gas fraction. The flow leaves the phase envelope, or enters the gaseous state at almost identical thermodynamic state for all EoS. Figure 4(c) shows the comparison of decompression wave speed ($u_w$) for various EoS, as it is highly dependent on speed of sound predicted by the EoS for a given phase and needed for ductile fracture studies, for instance. It is observed that the CPA model predicts the decompression wave speed and the transition pressure well, while the PC-SAFT over predicts the decompression wave speed and the PR over predicts the transition pressure. Hence, CPA is the chosen EoS for further studies conducted in this work.

4.3. Effect of Impurity

4.3.1. Comparison with Pure $C{O}_2$

This section compares the cases between pure $C{O}_2$ and the $C{O}_2+N_2(1.8\%mol)$ mixture, i.e. cases C0 and C1 in Table 2. It should be noted that a PH-Flash algorithm is used in the IFPEN-Carnot library [69] to determine the phase transition and equilibrium condition for pure $C{O}_2$ through the generation of an appropriate table of thermodynamic properties, as discussed in Section (Section 2.2). Figure 5 shows the comparison of thermodynamic state, pressure and gas fraction profiles, and decompression wave speed for pure and rich $C{O}_2$ mixtures. The saturation point at which the vaporization starts is well predicted for both pure $C{O}_2$ and $C{O}_2+N_2(1.8\%mol)$ mixture. As shown in Figure 5(a), the depressurization rate is slightly higher for pure $C{O}_2$ in comparison to $C{O}_2$ with impurities, as the steepness and the saturation point is lower for pure $C{O}_2$. The gas volume fraction at chocked flow state is also higher for pure $C{O}_2$ in comparison to rich $C{O}_2$ mixtures (Figure 5(b)). In this Figure, the two reflected pressure waves from the pipeline end can be clearly seen, without spurious oscillations induced by the impurity for the rich $C{O}_2$ mixture. In addition, the transition pressure for decompression wave speed ($u_w$) for pure $C{O}_2$ is smaller than with impurities, as depicted from the phase envelope as well. The thermodynamic state of entering the gaseous state is close for pure and rich $C{O}_2$ mixtures. However, the saturation pressure is slightly lower for pure $C{O}_2$, as expected (Figure 5(c)). Finally, one may note that the present simulation predicts larger chocked flow state in comparison to the experiments [31,70], which could be because of the effect of venting and pipeline divergence at the outlet.

4.3.2. Comparison Between Different Impurities $N_2,C{H}_4,Ar$

As discussed in the introduction Section 1, the selected binary systems are $C{O}_2-N_2$, $C{O}_2-C{H}_4$, and $C{O}_2-Ar$ representing the major non-condensable impurities relevant to current and emerging $sC{O}_2$ capture, transport, and storage applications [72]. The depressurization behavior of multi-component fluids is governed by a combination of thermodynamic, caloric, and transport properties, including mixture molar mass, critical properties, acentric factor, phase equilibrium characteristics, speed of sound, Joule–Thomson coefficient, and heat capacity. Figure 6 shows the comparison of thermodynamic state, pressure, gas fraction, and decompression wave speed for cases C4, C6, and C8. Since the speed of sound $C_s$ and the critical compressibility factor $Z_c$ are close to each other for all the mixtures at the given initial condition which is very close to critical point for supercritical flows, the depressurization rate and the pressure profile are also close to each other until the vaporization begins, as shown in Figure 6(a). Indeed the point of vaporization obviously differs for each case, but the thermodynamic state paths inside the phase envelope collapse again before the dew point and the gaseous phase, as shown also in Figure 6(a). The transition pressure of the decompression wave speed decreases in the same order $N_2>Ar>C{H}_4$ than the bubble saturation point in Figure 6(c). Similarly, the decompression wave speed for the gaseous phase in the Figure 6(c) also increases in the same order $N_2>Ar>C{H}_4$. Finally, changing the impurity does not induces very different results near the critical point during depressurization, although few differences in the pressure and temperature evolutions mainly during the two-phase stage, as shown in Figure 6(b,d).

4.3.3. Role of Impurity Mass Fraction

The comparison of Figure 6 and Figure 7 shows the role of the mole fraction on depressurization thermodynamic path, pressure, gas-fraction and decompression wave speed. The increase in mole fraction of impurities shifts the phase envelope upwards, accordingly shifting the onset pressure of vaporization, chocked flow pressure and gas fraction. From Figure 6(c) and Figure 7(c), for low composition mixtures, the decompression wave speed in the dense phase at $t=0$ is closer to the pure $C{O}_2$ and remains the same for all the mixtures as $410m/s$. However, for higher impurities concentration, the effect of impurities on wave speed becomes visible. The decompression wave speed for dense phase decreases as: $C{H}_4(427m/s)<Ar(410m/s)<N_2(400m/s)$. It should be noted that the minimum temperature value for the RFM tabular properties for the case of $C{O}_2+C{H}_4$ mixtures is $190K$. Hence, beyond this limit the sharp decrease in temperature is observed which is omitted from Figure 6(d) and Figure 7(d) as those are numerical artifacts induced by extrapolation of data from the out of the range of the RFM thermodynamic table. However, as shown in the pressure and gas fraction plot (Figure 6(b) and Figure 7(b)), the chosen range covers the entire two-phase transition region.

4.4. Effect of Outlet Geometry and Diameter

Cases C10-C13 in Table 3 analyze the effect of outlet geometry and diameter. The minimum cross-section area in the flow circuit determines the chocking mass flow rate [71]. The results presented in this section follows a similar trend as reported in the experimental work of Hammer et al. [71]. However, a slight difference in the chocked flow state is expected as the present simulation are with $N_2$(1.8% by mol) as an impurity. Figure 8 and Figure 9 illustrate that changing mass flow rate by changing the outlet geometry and diameter brings changes in the dense phase of the fluid, however as the vaporization begins the system follows the same thermodynamic path. The decrease in mass flow rate, reduces the rate of depressurization as shown by the pressure profile. This tends the thermodynamic path to shift right ward, representing higher density for a reduced mass flow rate. The decompression wave speed does not have a clear demarcation for the dense phase and gaseous phase for these cases as the fluid remains in the two phase state for the simulated time. Also prolonged chocked flow state, brings difficulty in estimating instantaneous time for a given pressure level, and hence the difficulty in calculating accurate decompression wave speed. Hence, the decompression wave speed plot is omitted for the present section. It is also noteworthy that the simulations become numerically more unstable under choked flow conditions, thereby requiring a smaller time step $\left(dt\right)$ to maintain stability.

4.5. Mass-Flow Rate and Mach Number

The primary objective of this work is to develop a database of mass flow rate, pressure, and gas-phase properties that can serve as inputs for dispersion models predicting the spread of $C{O}_2$ mixtures during pressure leakage events. In current industrial practice, discharge predictions largely rely on steady or quasi-steady formulations and empirical choked-flow correlations. These approaches depend strongly on the mixture speed of sound and its variation during phase change. Therefore, the present study reports both mass flow rate and Mach number for comparison. All reported quantities in Figure 10 correspond to mass-averaged values evaluated at the rupture/outflow plane. Density and velocity are coupled to each other, variation is always in both to maintain continuity. Initially the density is higher representing the presence of liquid, then reducing to air density. Accordingly the velocity varies, the variation in pressure, temperature in previous plots and mass fraction in Figure 10 and Figure 11 represent the sudden change in gas-fraction caused during two instances of the simulation: when released from the chocked flow state and when phase changes from two-phase to gas phase. The initial chocked flow condition is the duration at which the fluid exits at the maximum mass flow rate. This maximum mass flow rate for each impurity is found to be less dependent on mixture composition (i.e. mol % 3.6-5.4) and majorly the function of initial (supercritical) and outflow conditions only. However, in comparison of the different impurities this maximum mass flow rate decreases as $N_2>Ar>C{H}_4$. The flow enters the two-phase state at sub-sonic conditions but close to Mach 1 (i.e. $M_{en}\simeq 0.8$), while close to the gaseous state the fluid Mach number rises rapidly to balance the decrease in density and maintain the continuity. The maximum Mach number is found to be close to 1.7 (i.e. $M_{out}\simeq 1.7$), and is attained in the two-phase region before to transition to the gaseous phase where the Mach number decreases abruptly due to the higher speed of sound in the gaseous state. Similarly, Figure 11 reports the mass flow rate and mach number for various outlet diameter and geometries in the present study. The mass flow rate for the nozzle is found to be closer to the experimental data, but not for the orifice. This could be because of the sharp change in flow field for the case of orifice in comparison to the nozzle. The effect of vena-contracta and coefficient of discharge is severe for the case of orifice and need to be considered and modeled more accurately.

5. $sC{O}_2$ Pipeline Transport: Industrial Scale Simulation

For real life applications, the pipeline transport expands to the order of km instead of m. For such configuration, a 50km long pipeline depressurization case has been studied in [84]. A simulation of this scale generates increased complexity, particularly in terms of computational costs. Wall heat transfer model, pressure loss along the pipeline, venting tubes, turbulence, gravity effects cannot be neglected anymore.

5.1. Setup of the 50km Pipeline with Venting

A rectangular geometry is chosen instead of a circular cylinder for simplicity for the case of 50km long pipe depressurization, indeed keeping the effective area identical to the actual circular pipeline. From few preliminary simulations, it is observed that the venting configuration plays a significant role in the depressurization dynamics, as the venting pipe diameter ($d=7\mathrm{inch}$) is substantially smaller than the main pipeline diameter ($d=24\mathrm{inch}$). This pronounced area contraction governs the discharge characteristics and strongly influences the pressure decay rate. Therefore, two different venting configurations i.e. horizontal and vertical venting have been investigated taking advantage of the 3D-based finite volume method used in this work. A detailed description of both configurations are shown in Figure 12. The boundary conditions are similar to the previous short pipelines setup, i.e. Dirichlet for pressure at the outlet, isothermal wall and no slip condition at the rest of the wall. Boundary conditions along-with initial condition are also shown in Figure 12.

5.2. Pressure and Temperature Profiles

It can be seen from Figure 13(a) and Figure 13(b) that the vertical venting configuration outperforms the horizontal one to great extent, implying not only that the venting area is important but also the role of gravity for such configuration. The pressure profiles tends to follow the similar trend as the experimental for both the configurations, with vertical venting being more close to the experimental data. Some deviations begin when the simulation is close to 2 hours in physical time scale, as the effects of heat transfer and pipe loss are expected to be predominant after such long duration. Hence, the future work is to analyze the problem with a coupled conjugate heat transfer model along with skelitic pipeline configuration. The horizontal venting configuration requires approximately 2 days of wall-clock time on 256 cores to simulate 2 hours of physical time. In contrast, the vertical venting configuration requires nearly 8 days under identical computational resources and simulation duration. The extended computational time is primarily due to the severe time-step restriction imposed by the physics of the problem. The allowable time step is limited by stability constraints, which are directly linked to the smallest cell volume in the domain. Larger time steps would require larger cell volumes; however, this is not feasible in the present case. In the vertical venting configuration, resolving gravity-driven stratification and slip effects in the venting section requires to maintain a minimum of 10 cells along the vertical direction. This refinement significantly reduces the minimum cell volume in the domain, thereby imposing a much smaller stable time step and consequently increasing the total computational time.

6. Conclusions

The present work shows that the real-fluid (RFM) equilibrium framework is capable to model supercritical $C{O}_2$ pipeline transport problems with good accuracy. A good match was found between simulation results and experimental data. First, the isothermal heat transfer model and the CPA-EoS have been proven to be best suitable for $sC{O}_2$ fluid simulations. Next, the role of impurities has been investigated. Although the impurities not only shift thermodynamic boundaries but also alter the kinetics of phase transitions, it is revealed that the amount of impurities is less predominant near the critical condition operations, as the depressurization rate, wave speed and pressure profiles are found to be similar. Besides, the outlet geometry tends to impact the thermodynamic path to a great extent. It is found that for the nozzle geometries where there is a smooth change in the cross-section area, the simulation predictions were quite close to the experiment. However, for the case of the orifices where there is a sharp change in the cross-section area, the simulation under predict the leakage mass flow rate, implying a necessary computationally refined discretization for the outlet geometry. Finally, the role of venting area and gravity is found to be predominant for the 50km $sC{O}_2$ pipeline depressurization, to capture the initial physics of depressurization. The conjugate heat transfer (CHT) model and modeling of losses in pipes at the wall level become essential for subsequent stages when considering buried pipelines.