1 Introduction↩︎

According to the IEA [1], carbon capture and storage (CCS) is expected to play an important role in addressing the challenge of reducing greenhouse gas emissions. In CCS, is captured from industrial sources and transported to storage sites, e.g., depleted gas fields, where it is injected underground for permanent storage. A major share of this transport is expected to be carried out via pipelines. The captured from industries is generally not pure and can contain varying levels of impurities such as methane, nitrogen, amines, sulfur oxides, and water etc. [2]. These impurities can significantly affect phase behavior, thermophysical properties, and consequently, the design and safe operation of pipeline systems.

One of the extreme events in pipeline operations is depressurization, which can be required following a shut-in operation [3]–[5]. In such events, the temperatures and pressures can reach extreme values with the risk of exceeding the material operational limit of the pipe wall. Accurately predicting this behaviour is essential for safe pipeline design. In this paper, we investigate the depressurization from a pipeline carrying with impurities.

Some experimental studies have investigated decompression behaviour of –rich mixtures, providing valuable reference data for model validation. Drescher et al. [6] provided experimental data for mixtures. Cosham et al. [7] presented 14 shock-tube experiments, including tests with pure as well as containing impurities, to study decompression behaviour. Botros et al. [8] reported experimental data for mixtures to determine decompression wave speeds.

Two-phase flow of fluid mixtures is described by the Navier–Stokes equations. For general two-phase flow models, the reader is referred to Stewart et al. [9], and for transport specifically, to the review by Munkejord et al. [10] and the references therein. For pipeline applications, cross-sectional averaging is typically applied to obtain one-dimensional two-fluid models [11]. Such averaging may render the system of equations ill-posed. To avoid these issues and focus more on the integration between the fluid dynamics and thermodynamics, we employ the Homogeneous Equilibrium Model (HEM), which is hyperbolic and guarantees well-posedness. In HEM, the fluid is treated as a homogeneous mixture of two well-mixed phases, with the Navier-Stokes equations formulated directly in terms of mixture density, velocity, and energy. For an inviscid flow in a frictionless pipe, this formulation reduces to Euler’s equations, where all conserved quantities are expressed in terms of mixture properties.

In addition to the fluid dynamics, the phase behaviour of the real fluids plays a crucial role and is governed by thermodynamic principles. Depending on the operating conditions, transport in pipelines can either be single-phase or two-phase. To describe two-phase flow accurately, a real-gas equation of state (EOS) is required. Often, the Mie-Grüneisen stiffened gas EOS is used; see [12]–[16]. However, to accurately describe the real fluid behavior, more advanced EOSs are required. These advanced EOSs are defined in terms of the Helmholtz free energy function: given the temperature, volume, and composition of a single-phase, they allow calculation of all thermodynamic properties. It is important to note that such EOSs are defined at the single-phase level; thermodynamic properties cannot be directly evaluated from overall mixture quantities in a two-phase system without first performing a phase-split calculation. In a fluid dynamic simulation, the mixture properties, namely the volume (mesh size), velocity, and energy, are known at each time step. From these quantities, one intends to compute the mixture pressure. Naturally, this leads to a phase-split calculation, commonly referred to in thermodynamic literature as a flash calculation.

A substantial body of literature exists on flash calculations. The majority focuses on the PTN-flash (pressure–temperature specified) problem, (see [17]–[20]. The second most documented case is the VTN-flash (volume–temperature specified), (see [21]–[24]). Both PTN and VTN are isothermal flash problems. The literature, however, is rather scarce on the non-isothermal flashes, for instance, the PHN-flash (pressure–enthalpy specified), the UVN-flash (internal energy–volume specified), or the SVN-flash (entropy–volume specified). For fluid dynamical simulations, UVN-flash provides a natural formulation, since the internal energy, specific volume, and molar composition are available in each computational cell. Before carrying out the flash calculation, it is essential to know whether the fluid mixture is in a single-phase or two-phase condition. This is determined via stability analysis. If stability analysis predicts the existence of a two-phase state, the flash calculation is performed. Moreover, the outcome of the stability analysis often provides an excellent initial guess for the subsequent flash calculation.

Some researchers have investigated the numerical simulation of –transport in pipelines. Elshahomi et al. [25] presented numerical simulations of decompression-wave speeds using ANSYS software. Over the past several years, SINTEF Energy Research in Norway has been a leading contributor in this area. They have developed a dedicated flow loop for experiments and have published a steady stream of work on transport; see [3], [10], [26]. Their work covers both numerical modeling of transport with and without impurities and detailed studies of thermodynamic behaviour. Most studies focus either on the fluid dynamics; see [3], [26]–[28] and use an existing thermodynamic library as-is, or they focus on the thermodynamic behavior [29]–[31].

In this paper, we present a comprehensive treatment of two-phase multicomponent transport of fluid mixtures in pipelines that addresses both the fluid dynamics aspects as well as the thermodynamic aspects with a UVN-flash framework. In particular, we focus on the robustness of the flash problem for challenging inputs from the fluid dynamics solver. To address these challenges, we propose a reformulation of the UVN-flash problem through an alternative set of variables. This transformation improves the robustness of phase-equilibrium calculations. The key contribution of this paper is the development and testing of the flash routine tailored for dynamic pipeline transport models. The proposed methodology is tested using both tank and pipeline depressurization cases.

The paper is structured as follows. We begin with the governing equations for HEM model and the dynamics of Tank depressurization in 2. In 3, we discuss the thermodynamic aspects including the calculation of the thermodynamic properties and the choice of the EOS. Stability analysis and flash calculations for phase-split problems are discussed in 4 and 5, respectively. The time discretization method is addressed in 6. Finally, 7 presents the results and compares them with literature data, followed by the conclusions in 8.

2 Governing Equations↩︎

This section describes the governing equations for the fluid flow.

2.1 Fluid flow equations in pipeline↩︎

The Navier-Stokes(NS) equations describe the dynamics of fluid flow. In principle, one could solve the 3D equations; however, for practical applications where a pipeline could be many kilometers long, it becomes computationally prohibitive. For these situations, a more efficient approach is to use a cross-sectionally averaged 1D equations [32]. For a fully dispersed two-phase fluid flow, the governing equations can be written in terms of average mixture quantities, leading to the Homogeneous Equilibrium Model (HEM). For inviscid flow through a horizontal frictionless pipe, neglecting heat transfer effects, the HEM can be expressed as:

\[\label{Eqns:NS} \partial_{t} {\mathcal{U}}+ \partial_{x} {\mathcal{F}}({\mathcal{U}}) = \mathcal{S}({\mathcal{U}}),\tag{1}\] where \({\mathcal{U}}(x,t)\) is the vector of conserved variables, \({\mathcal{F}}({\mathcal{U}})\) the flux vector and \(\mathcal{S}({\mathcal{U}})\) the source terms. These terms have the following form: \[\label{vector95eqns} {\mathcal{U}}= \begin{bmatrix} {\rho }\\ {\rho u }\\ {\rho E} \end{bmatrix}, \quad {\mathcal{F}}({\mathcal{U}}) = \begin{bmatrix} {\rho u}\\ {\rho u^{2} + p }\\ {(\rho E + p) u} \end{bmatrix} \quad \mathcal{S}({\mathcal{U}}) = \begin{bmatrix} {0}\\ {0}\\ {0} \end{bmatrix}.\tag{2}\] Here, \(\rho, u, E, p\) denote the mixture mass density, mixture velocity, specific total energy, and the pressure, respectively. The HEM model assumes instantaneous thermodynamic equilibrium and that both phases travel at the same speed. The total specific energy \(E\) can be expressed as \[\begin{align} E = e + \frac{1}{2} u^2, \end{align}\] where \(e\) is the specific internal energy. The mixture properties can be written in terms of the phasic properties as follows \[\begin{align} \rho &:= \alpha \rho_{g} + (1 - \alpha) \rho_{l}, \tag{3} \\ \rho e &:= \alpha \rho_{g} e_{g}+ (1 - \alpha) \rho_{l} e_{l}, \tag{4} \end{align}\] where \(\alpha\) denotes the volume fraction of the gas phase, \(\rho_{g}\) and \(\rho_{l}\) denote the mass density of the gas and liquid phase, respectively, \(e_{g}\) and \(e_{l}\) denote the specific internal energy (mass basis) of the corresponding phases.

The system of equations 1 is closed by incorporating an equation of state(EOS), which relates pressure to the other thermodynamic quantities. The role of the EOS in closing the system will be discussed in detail in 3. Furthermore, for the problem to be well-posed, appropriate initial and boundary conditions must be specified; these aspects will be discussed in 7.

2.2 Spatial Discretization↩︎

The HEM model described in the 2.1 is solved using the Finite Volume Method (FVM). As a first step, we discretize the system 1 in space, which yields the following semi-discrete form: \[\begin{align} \label{eqn:semiDiscrete} \frac{\mathrm{d}\mathbf{U}_i}{\mathrm{d}t} = -\frac{1}{\triangle x} (\Hat{{\mathcal{F}}}_{i+\frac{1}{2}} - \Hat{{\mathcal{F}}}_{i-\frac{1}{2}}), \qquad i=1\ldots N, \end{align}\tag{5}\] where \(\mathbf{U}_i(t) \in \mathbb{R}^{3} \approx {\mathcal{U}}(x_{i},t)\), represents the cell-averaged conserved variables in cell \(i\), \({\triangle x}\) is the spatial grid spacing, \(\Hat{{\mathcal{F}}}_{i\pm\frac{1}{2}}\) denote the numerical fluxes at the cell interfaces \(i\pm\frac{1}{2}\) and \(N\) is the total number of finite volume cells. System 5 for the full computational domain can be expressed compactly in vectorial notation as: \[\label{eq:pipe:vector95form} \frac{\mathrm{d}\mathbf{U}(t)}{\mathrm{d}t} = \mathbf{f}_{\scriptscriptstyle\text{Pipe}}(\mathbf{U}(t), \mathbf{V}(t)),\tag{6}\] where \(\mathbf{U}(t) \in \mathbb{R}^{3N}\) is the global vector of conserved quantities for the entire domain(pipe), defined as: \[\begin{align} \mathbf{U}(t) = [\mathbf{U}_{1}(t), \ldots, \mathbf{U}_{i}(t), \ldots, \mathbf{U}_{N}(t)]^T, \qquad \mathbf{U}_{i} = [\rho_i, (\rho u)_i, (\rho E)_i]^T, \end{align}\] Here, \(\mathbf{U}(t)\) is composed of \(N\) individual cell vectors \(\mathbf{U}_i(t)\), and \(\rho_i, (\rho u)_i\) and \((\rho E)_i\) denote the mass density, momentum, and the total energy in the \(i^{th}\) cell. The vector \(\mathbf{V}(t)\) contains the non-conservative (algebraic) variables associated with the thermodynamic state in each cell. Its precise definition will be provided in  6.
The right-hand side \(\mathbf{f}_{\scriptscriptstyle\text{Pipe}}\) depends on both the conservative variables \(\mathbf{U}(t)\) and the associated non-conservative variables \(\mathbf{V}(t)\), and is defined as \[\begin{align} \mathbf{f}_{\scriptscriptstyle\text{Pipe}}(\mathbf{U}(t), \mathbf{V}(t)) := [\mathbf{f}_{\scriptscriptstyle\text{Pipe, 1}}(t), \ldots, \mathbf{f}_{\scriptscriptstyle\text{Pipe, i}}(t), \ldots, \mathbf{f}_{\scriptscriptstyle\text{Pipe, N}}(t)]^T, \end{align}\] where each local \(\mathbf{f}_{\scriptscriptstyle\text{Pipe, i}}\) is expressed as \[\begin{align} \mathbf{f}_{\scriptscriptstyle\text{Pipe, i}}(\mathbf{U}_{i-r:i+r}, \mathbf{V}_{i-r:i+r}) = -\frac{1}{\Delta x} \left( \Hat{{\mathcal{F}}}_{i+\frac{1}{2}} - \Hat{{\mathcal{F}}}_{i-\frac{1}{2}} \right). \end{align}\] where \(r\) denotes the stencil radius (half-width), \(\mathbf{U}_{i-r:i+r} := (\mathbf{U}_{i-r}, \ldots, \mathbf{U}_{i+r})\), and similarly for \(\mathbf{V}_{i-r:i+r}\). The numerical fluxes \(\Hat{{\mathcal{F}}}_{i\pm\frac{1}{2}}\) are computed using data from a stencil of width \(2r + 1\), and hence depend on the conservative variables \(\mathbf{U}_j\) and the corresponding non-conservative variables \(\mathbf{V}_j\) for \(j = i - r, \ldots, i + r\). For a first-order scheme, we put \(r = 1\). To compute the inter-cell numerical fluxes \(\Hat{{\mathcal{F}}}_{i\pm\frac{1}{2}}\), we employ the HLLC scheme; (see 10).

2.3 Tank depressurization↩︎

The tank model can be regarded as a pipeline discretised with a single computational cell. Many researchers have previously considered the dynamic simulation of the tank in the context of UVN-flash[33]–[36]. Furthermore, this model has also been employed by [37]–[39] in the context of pipeflows. The governing equations describing the evolution of the molar composition and internal energy of the tank are given by:

\[\begin{align} &\frac{\mathrm{d}N_i}{\mathrm{d}t} = \dot{N}_{i}^{\scriptscriptstyle\text{in}} - \dot{N}_{i}^{\scriptscriptstyle\text{out}} \quad \text{for } i = 1,\dots,n \tag{7}, \\ &\frac{\mathrm{d}U}{\mathrm{d}t} = \dot{H}^{\scriptscriptstyle\text{in}} - \dot{H}^{\scriptscriptstyle\text{out}} + \dot{Q}, \tag{8} \end{align}\] where \(n\) denotes the total number of components in the mixture. The variable \(N_{i}\) represents the moles of component \(i\) in the tank, while \(U\) is the total internal energy. The terms \(\dot{N}_{i}^{\scriptscriptstyle\text{in}}\) and \(\dot{N}_{i}^{\scriptscriptstyle\text{out}}\) denote the molar flow rates of component \(i\) in the input and output streams, respectively. Likewise, \(\dot{H}^{\scriptscriptstyle\text{in}}\) and \(\dot{H}^{\scriptscriptstyle\text{out}}\) denote the enthalpy flow rates associated with the input and output streams, and \(\dot{Q}\) is the rate of heat supplied to the tank.

The equations 7 —8 can be written in compact vector form as: \[\label{eq:tank:vector95form} \frac{\mathrm{d} \mathbf{U}}{\mathrm{d} t} = \mathbf{f}_{\scriptscriptstyle\text{Tank}}(\mathbf{U}(t), \mathbf{V}(t)),\tag{9}\] where the state vector \(\mathbf{U}(t) \in \mathbb{R}^{n+1}\) is defined as \[\begin{align} \mathbf{\mathbf{U}}(t) = \begin{bmatrix} N_1(t) \\ \vdots \\ N_n(t) \\ U(t) \end{bmatrix}, \end{align}\] and the vector \(\mathbf{V}(t)\) contains the variables associated with the thermodynamic state. This will be revisited in 6. The right-hand side \(\mathbf{f}_{\scriptscriptstyle\text{Tank}}(\mathbf{U}, \mathbf{V})\) is given by \[\begin{align} \label{eq:tank95rhs} \mathbf{f}_{\scriptscriptstyle\text{Tank}}(\mathbf{U}, \mathbf{V}) = \begin{bmatrix} \displaystyle \dot{N}_{1}^{\mathrm{in}} - \dot{N}_{1}^{\mathrm{out}} \\ \vdots \\ \displaystyle \dot{N}_{n}^{\mathrm{in}} - \dot{N}_{n}^{\mathrm{out}} \\ \displaystyle \dot{H}^{\mathrm{in}} - \dot{H}^{\mathrm{out}} + \dot{Q} \end{bmatrix}. \end{align}\tag{10}\]

3 Thermodynamical aspects↩︎

3.1 Definitions↩︎

For the sake of clarity, we define the following thermodynamic constructs under specified constraints of total internal energy \((U)\), volume \((V)\), and mole numbers \(\mathbf{N}\).

Reference Phase: The reference phase (\(\star\)) represents a hypothetical single-phase state characterized by the total internal energy \(U^{\star}\), volume \(V^{\star}\) and total mole numbers \(\mathbf{N}^{\star} = ( N_1^{\star}, \dots, N_n^{\star})\).

Trial Phase: The trial phase is an incipient(very small amount) phase introduced to assess the thermodynamic stability of a system. The system is perturbed by adding this phase and if the total entropy increases, the system is considered unstable as a single phase, and phase separation is favorable.

Stability Analysis: Stability analysis determines whether a given fluid at specified \(U, V\) and \(\mathbf{N}\) will remain stable as a single phase or separate into two phases, i.e., into a gas-liquid mixture. If the single-phase state is unstable, this analysis provides estimates of its temperature, molar concentration (moles per unit volume), and molar internal energy density. These estimates are used to obtain initial guesses for the subsequent flash calculation. Please see the 4 for more details.

Flash: A flash calculation determines the equilibrium phase split of a multicomponent mixture. Under specified thermodynamic constraints (e.g., UVN), it computes the distribution of extensive quantities among the coexisting phases, i.e., the enthalpy, internal energy, entropy, and volume of each phase, and the common intensive quantities, e.g., temperature, pressure, and component-wise chemical potential, which satisfy the conditions of thermal, mechanical, and chemical equilibrium.

3.2 Equation of State↩︎

To close the HEM model and the tank model, we need a constitutive relation to compute the pressure in the pipe model and the enthalpy flow rate in the tank model. This constitutive model is based on the thermodynamic properties of the fluid mixture and is commonly referred to as the Equation of State(EOS). The EOS for the fluid mixture is specified in terms of any two thermodynamic variables, such as pressure and temperature or volume and temperature, etc., in addition to the mole numbers of the components in the mixture. Thus, for an \(n\)-component mixture, \(n + 2\) variables are required to fully specify the thermodynamic state of the system. Based on the choice of the two input thermodynamic variables, the EOS can be classified into two categories: pressure-based and volume-based. Typically, the real fluid EOS is defined in terms of the Helmholtz free energy \(A\) as a function of volume \(V\), temperature \(T\), and the molar composition vector \(\mathbf{N}\): \[A = A(T, V, \mathbf{N}), \label{eqn:helmholtz95free95energy}\tag{11}\] where \(\mathbf{N}:= \{ N_1, \dots, N_n\}\) and \(N_i\) denote the number of moles of the \(i^{th}\) component in the fluid mixture. We use Peng–Robinson EOS [40] in this work. The various thermodynamic properties can be computed directly from the Helmholtz energy as below.

\[\begin{align} p &= -\left(\frac{\partial A}{\partial V}\right)_{T, \mathbf{N}} \\ S &= -\left(\frac{\partial A}{\partial T}\right)_{V, \mathbf{N}} \\ U &= A + T S \\ H &= A + T S + p V \\ C_v &= \left(\frac{\partial U}{\partial T}\right)_{V, \mathbf{N}} \\ \mu_i &= \left(\frac{\partial A}{\partial N_i}\right)_{V, \{N_j\}_{j \ne i}} \end{align}\] where \(p, S, U, H, C_v\) and \(\mu_i\) denote the pressure, entropy, internal energy, enthalpy, heat capacity at constant volume, and the chemical potential of the \(i^{th}\) component. Here, it is understood that all derivatives and functions are evaluated as functions of \(T\), \(V\), and \(\mathbf{N}\).

An EOS is defined for individual homogeneous phases and must be applied to phasic variables, namely, the temperature, volume, and molar composition of each phase. It cannot be applied directly to mixture-averaged quantities in a multiphase system. However, in practical scenarios, only the extensive properties such as the total volume and the overall molar composition are known, and the distribution of the total volume and the mole numbers among the individual phases is not available. This necessitates a phase-split computation whereby one determines the distribution of the extensive quantities among the coexisting phases and the common intensive quantities (more on this later). The phase-split calculations are commonly referred to in the thermodynamic literature as the Flash calculations.

3.3 Choice of Flash↩︎

Different types of flash problems arise depending upon which two thermodynamic variables are specified. The most common flash problem addressed in the literature is the PTN-Flash [17], in which the pressure, temperature, and the overall molar composition are specified. The objective of the PTN-flash is to find the equilibrium phase split and compute the extensive properties of each phase, such as the internal energy, enthalpy, and volume etc. Another important flash problem is the VTN-Flash, where the total volume of the mixture is specified along with the temperature and overall composition, see for example [41]. Both PTN and VTN-flashes are commonly referred to as isothermal flashes since the temperature is fixed.

In closed systems, for example, the energy balance calculation in process design, the relevant specification often includes the total internal energy \(U\), total volume \(V\), and overall molar composition \(\mathbf{N}\). This non-isothermal specification corresponds to the so-called UVN-flash problem. In our tank model and the HEM model, we typically have access to the total volume (fixed for the tank, cell volume for the pipe), total internal energy, and molar composition at each time step of the simulation. Thus, UVN-Flash formulation provides a natural approach to determine equilibrium phase-split and, consequently, the pressure in such scenarios.

4 Stability Analysis↩︎

In this section, we briefly review the stability analysis procedure. For a detailed account on this topic, we refer the interested reader to the existing literature [18], [24], [42], [43].

4.1 UVN stability↩︎

To assess the thermodynamic stability of a single-phase system with fixed total internal energy \(U^\star\), volume \(V^\star\), and mole numbers \(\mathbf{N}^\star= (N_1^\star, \dots, N_n^\star)\), we consider a hypothetical phase separation into two phases: a bulk phase \(\beta\) and an incipient (infinitesimal) phase \(\varepsilon\). The total entropy of the resulting two-phase system is given by:

\[\begin{align} S_{\text{two-phase}} = S(U^{\scriptscriptstyle\beta}, V^{\scriptscriptstyle\beta}, \mathbf{N}^{\scriptscriptstyle\beta}) + S(U^{\varepsilon}, V^{\varepsilon}, \mathbf{N}^{\varepsilon}) \end{align}\] subject to conservation constraints: \[\begin{align} U^\star&= U^{\scriptscriptstyle\beta} + U^{\varepsilon}, \\ V^\star&= V^{\scriptscriptstyle\beta} + V^{\varepsilon}, \\ N_i^\star&= N_i^{\scriptscriptstyle\beta} + N_{i}^{\varepsilon}, \qquad i=1 \dots n. \end{align}\] The entropy of the reference phase is: \[\begin{align} S_{\text{ref}} = S(U^\star, V^\star, \mathbf{N}^\star). \end{align}\] The bulk phase entropy can be computed by performing a first-order Taylor expansion of the entropy around the reference state \((U^\star, V^\star, \mathbf{N}^\star)\), treating the incipient phase contribution as an infinitesimal perturbation: \[\begin{align} S^{\scriptscriptstyle\beta} &:= S(U^\star- U^{\varepsilon}, V^\star- V^{\varepsilon}, N_1^\star- N_{1}^{\varepsilon}, \dots, N_n^\star- N_{n}^{\varepsilon}), \\ &= S(U^\star, V^\star, \mathbf{N}^\star) - U^{\varepsilon} \left(\frac{\partial S}{\partial U}\right)_{V, \mathbf{N}} - V^{\varepsilon} \left(\frac{\partial S}{\partial V}\right)_{U, \mathbf{N}} - \sum_{i=1}^n N_{i}^{\varepsilon} \left(\frac{\partial S}{\partial N_i}\right)_{U, V, \mathbf{N}}. \end{align}\] This leads to \[\begin{align} S_{\text{two-phase}} = S(U^\star, V^\star, \mathbf{N}^\star) - U^{\varepsilon} \left(\frac{\partial S}{\partial U}\right)_{V, \mathbf{N}} - V^{\varepsilon} \left(\frac{\partial S}{\partial V}\right)_{U, \mathbf{N}} - \sum_{i=1}^n N_{i}^{\varepsilon} \left(\frac{\partial S}{\partial N_i}\right)_{U, V, \mathbf{N}} + S(U^{\varepsilon}, V^{\varepsilon}, \mathbf{N}^{\varepsilon}). \end{align}\] The entropy difference \(\Delta S = S_{\text{two-phase}} - S_{\text{ref}}\) is then: \[\begin{align} \label{eq:deltaS951} \Delta S = - U^{\varepsilon} \left(\frac{\partial S}{\partial U}\right)_{V, \mathbf{N}}(U^\star, V^\star, \mathbf{N}^\star), - V^{\varepsilon} \left(\frac{\partial S}{\partial V}\right)_{U, \mathbf{N}}(U^\star, V^\star, \mathbf{N}^\star) \nonumber \\ - \sum_{i=1}^n N_{i}^{\varepsilon} \left(\frac{\partial S}{\partial N_i}\right)_{U, V, \mathbf{N}}(U^\star, V^\star, \mathbf{N}^\star) + S(U^{\varepsilon}, V^{\varepsilon}, \mathbf{N}^{\varepsilon}). \end{align}\tag{12}\] Using the standard thermodynamic identities \[\begin{align} \label{eq:thermo95ids1} \left(\frac{\partial S}{\partial U}\right)_{V, \mathbf{N}} = \frac{1}{T}, \quad \left(\frac{\partial S}{\partial V}\right)_{U, \mathbf{N}} = \frac{P}{T}, \quad \left(\frac{\partial S}{\partial N_i}\right)_{U, V, \mathbf{N}} = -\frac{\mu_i}{T}, \end{align}\tag{13}\] into 12 , we get \[\begin{align} \label{eq:deltaS952} \Delta S = - \frac{U^{\varepsilon}}{T^{\star}} - \frac{P^{\star}}{T^{\star}} V^{\varepsilon} + \sum_{i=1}^n N_{i}^{\varepsilon} \frac{\mu_i^{\star}}{T^{\star}} + S(U^{\varepsilon}, V^{\varepsilon}, \mathbf{N}^{\varepsilon}). \end{align}\tag{14}\] Furthermore, since entropy is homogeneous of degree one, Euler’s theorem gives the incipient phase entropy as: \[\begin{align} \label{eq:euler95homogeneous} S(U^{\varepsilon}, V^{\varepsilon}, \mathbf{N}^{\varepsilon}) &= \frac{U^{\varepsilon}}{T^{\varepsilon}} + \frac{P^{\varepsilon}}{T^{\varepsilon}} V^{\varepsilon} - \sum_{i=1}^{n} \frac{\mu_{i}^{\varepsilon}}{T^{\varepsilon}} N_{i}^{\varepsilon}. \end{align}\tag{15}\] Substituting 15 into 14 yields: \[\begin{align} \Delta S = \left(\frac{1}{T^{\varepsilon}} - \frac{1}{T^{\star}}\right)U^{\varepsilon} + \left(\frac{P^{\varepsilon}}{T^{\varepsilon}} - \frac{P^{\star}}{T^{\star}} \right)V^{\varepsilon} - \sum_{i=1}^n N_{i}^{\varepsilon} \left(\frac{\mu_{i}^{\varepsilon}}{T^{\varepsilon}} - \frac{\mu_i^{\star}}{T^{\star}}\right). \end{align}\] If \(\Delta S > 0\), the single-phase state is unstable, and the mixture will split into two phases. Dividing by \(V^{\varepsilon}\), we obtain the famous tangent plane distance (TPD) function

\[\begin{align} \label{eqn:tpd1} D(c_1^{\varepsilon}, \dots, c_n^{\varepsilon}, u^{\varepsilon}) &= \frac{\Delta S}{V^{\varepsilon}}, \nonumber\\ &= u^{\varepsilon} \left(\frac{1}{T^{\varepsilon}} - \frac{1}{T^{\star}}\right) + \left(\frac{P^{\varepsilon}}{T^{\varepsilon}} - \frac{P^{\star}}{T^{\star}} \right) - \sum_{i=1}^n c_{i}^{\varepsilon} \left(\frac{\mu_{i}^{\varepsilon}}{T^{\varepsilon}} - \frac{\mu_i^{\star}}{T^{\star}}\right), \end{align}\tag{16}\] where \(u^{\varepsilon} := U^{\varepsilon}/V^{\varepsilon}\) is the molar internal energy density and \(c_{i}^{\varepsilon} := N_{i}^{\varepsilon}/V^{\varepsilon}\) is molar concentration of \(i^{\text{th}}\) component in the incipient trial phase. Note that, the independent variables are \(u^{\varepsilon}\) and \(c_1^{\varepsilon}, \dots, c_n^{\varepsilon}\). The intensive properties of the incipient phase temperature \(T^{\varepsilon}\), pressure \(P^{\varepsilon}\), and chemical potentials \(\mu_i^{\varepsilon}\) can be obtained as follows. First, the temperature is determined by inverting the internal energy relation: \[u^{\varepsilon} = U(T^{\varepsilon}, 1.0, c_1^{\varepsilon}, \dots, c_n^{\varepsilon}).\] Once \(T^{\varepsilon}\) is known, the pressure can be computed as: \[P^{\varepsilon} = P(T^{\varepsilon}, 1.0, c_1^{\varepsilon}, \dots, c_n^{\varepsilon}),\] and the chemical potentials follow similarly: \[\mu_i^{\varepsilon} = \mu_i(T^{\varepsilon}, 1.0, c_1^{\varepsilon}, \dots, c_n^{\varepsilon}), \qquad i = 1, \dots, n.\]

The TPD function \(D\) can be used for stability testing for a system with a specified UVN. If \(D \le 0\) for all admissible \(\{c_{1}^{\varepsilon}, \dots, c_{n}^{\varepsilon}, u^{\varepsilon} \}\), then the system is stable as single-phase. To find out a state \(\{c_{1}^{\varepsilon}, \dots, c_{n}^{\varepsilon}, u^{\varepsilon} \}\) for which \(D > 0\), one seeks to find the local maxima of the function \(D\). Differentiating \(D\) with respect to its independent variables \(u^{\varepsilon}\) and \(\{c_1^{\varepsilon}, \dots, c_n^{\varepsilon}\}\) and applying the Gibbs—Duhem relation yields (for details, see [44], [45])

\[\begin{align} \frac{1}{T^{\varepsilon}} - \frac{1}{T^{\star}} &= 0, \tag{17}\\ \frac{1}{T^{\star}}(\mu_i^{\varepsilon} - \mu_i^{\star}) &= 0, \qquad i = 1, \dots, n. \tag{18} \end{align}\] 17 implies that the temperature of the incipient phase is equal to that of the reference phase temperature; 16 simplifies to \[\begin{align} \label{eqn:tpd2} D = \left(\frac{P^{\varepsilon}}{T^{\varepsilon}} - \frac{P^{\star}}{T^{\star}} \right) - \sum_{i=1}^n c_{i}^{\varepsilon} \left(\frac{\mu_{i}^{\varepsilon}}{T^{\varepsilon}} - \frac{\mu_i^{\star}}{T^{\star}}\right) \end{align}\tag{19}\]

Interestingly, this is identical to the TPD function obtained in VTN-stability analysis [22]. Therefore, in principle, one can solve a VTN-stability problem instead of a UVN-stability problem. The system of \(n\) equations 18 , can be solved for the incipient phase concentrations \(c_i^{\varepsilon},\; i = 1, \dots, n\) and if the resulting TPD for this set of concentrations is greater than zero, then the system is thermodynamically unstable and will split into multiple phases. Note that the temperature of the incipient phase, \(T^{\varepsilon}\), is the same as the reference phase temperature \(T^{\star}\). Furthermore, a good initial guess is critical to obtaining convergence of the stability analysis test, which is discussed next.

4.2 Initial guess for stability analysis↩︎

We employ three different types of initial guesses for the stability analysis test:

• Simplex-based approach, as discussed by [42], [43].

• Saturation pressure-based approach, as discussed by [18], [45].

• Gaussian random perturbations added to the hypothetical single-phase concentrations.

These three sets of initial conditions provide a diverse and well-distributed coverage of the solution space, which improves the robustness of the stability analysis. By exploring multiple starting points, we increase the likelihood of detecting incipient phase formation and obtaining a reliable estimate for the concentrations \(c_1^{\varepsilon}, \dots, c_n^{\varepsilon}\) and molar internal energy density \(u^{\varepsilon}\) of the trial phase. A practical consideration in implementing the stability analysis is to skip the trivial solutions, i.e., solutions where the computed trial phase molar concentration is identical to the input molar concentration.

4.3 Initial guess from stability results↩︎

A good initial guess is crucial for achieving convergence in flash calculations. Stability analysis provides the trial phase concentration and internal energy density. To initiate the flash calculations, we also need the volumetric phase split. Following Kumar et al. [43], the trial phase volume is initially set to half of the total volume and then iteratively halved until a phase split is found that yields a two-phase entropy greater than that of the corresponding single-phase state. The energy density and molar concentration can then be scaled by this trial volume to obtain initial estimates of the trial phase internal energy and mole numbers, and the corresponding bulk-phase quantities are determined by subtracting from the total quantities.

5 Reformulation of the UVN-flash for fluid dynamics↩︎

If stability analysis indicates that the mixture will split into two phases, flash is performed to determine the equilibrium phase split. At a given time step in a computational cell, the available fluid variables are the total mass density \(\rho\), the specific internal energy \(e\), and the overall mixture composition expressed as mole fractions \(\boldsymbol{z}:= \{z_1, \dots, z_n\}\). From these quantities, the total internal energy, \(U^{\scriptscriptstyle\text{tot}}\) and the total number of moles of \(i^{\scriptstyle\text{th}}\) component, \(N^{\scriptscriptstyle\text{tot}}_i\) can be computed as \[\begin{align} U^{\scriptscriptstyle\text{tot}} &= \rho e V, \\ N^{\scriptscriptstyle\text{tot}}_i &= z_i \frac{\rho V}{\sum_j z_j M_{j,w}}, \qquad i = 1,\dots,n, \end{align}\] where \(V\) is the cell volume and \(M_{j,w}\) denotes the molecular weight of component \(j\). The flow diagram for the flash calculation is provided in 1. In the following subsections, we present the UVN-flash formulation tailored for fluid-dynamical applications.

5.1 UVN-flash in TVN-Space↩︎

For a closed system whose total internal energy, volume, and mass (i.e., the number of moles) are fixed, the entropy tends to be at a maximum at equilibrium. Let the system exist in two phases, say gas and liquid. The total entropy of the system, \(S^{\scriptscriptstyle\text{tot}}\) can be expressed as

\[\label{Flash:entropy95uvn} S^{\scriptscriptstyle\text{tot}} = S(U^{g}, V^{g}, \mathbf{N}^{g}) + S(U^{l}, V^{l}, \mathbf{N}^{l}).\tag{20}\] The liquid phase quantities can be expressed in terms of the total and the gas phase quantities, i.e. \(U^{l} = U^{\scriptscriptstyle\text{tot}} - U^{g}, V^{l} = V - V^{g}\) and \(N_1^{l} = N_1^{\scriptscriptstyle\text{tot}} - N_1^{g}, \dots, N_n^{l} = N_n^{\scriptscriptstyle\text{tot}} - N_n^{g}\). 20 can be rewritten as \[\label{Flash:entropy95uvn2} S^{\scriptscriptstyle\text{tot}} = S(U^{g}, V^{g}, \mathbf{N}^{g}) + S(U^{\scriptscriptstyle\text{tot}} - U^{g}, V - V^{g}, \mathbf{N}^{\xi}),\tag{21}\] where \[\begin{align} \label{eqn:moles95remaining95phase} \mathbf{N}^{\xi} := \{N_1^{\scriptscriptstyle\text{tot}} - N_1^{g}, \dots, N_n^{\scriptscriptstyle\text{tot}} - N_n^{g}\}. \end{align}\tag{22}\] The entropy function 21 can be maximized, using the gas–phase variables as optimization variables, to determine the equilibrium phase split. This approach has been employed by Castier [46] and Smejkal et al. [42]. However, the EOS allows us to compute the entropy and other thermodynamic quantities as a function of \(T, V\) and \(\mathbf{N}\). Thus, at every sub-iteration of the maximization, the temperature must be obtained by solving \[\begin{align} \label{eqn:internal95energy95tvn} U = U(T, V, \mathbf{N}), \end{align}\tag{23}\] for a given internal energy \(U\), volume \(V\), and composition \(\mathbf{N}\). In a pipeline or a compositional reservoir simulation, where thousands or even millions of such computations need to be performed, these additional Newton iterations associated with this inversion step for determining the temperature, can become a bottleneck. This motivates the reformulation of the entropy maximization problem directly in terms of the natural EOS variables \(T, V\) and \(\mathbf{N}\). Using the thermodynamic relation \(U = A + TS\), the entropy for a given internal energy \(U^{\scriptscriptstyle\text{tot}}\) can be expressed as \[\begin{align} \label{SQ} S^{\scriptscriptstyle\text{TVN}} = \frac{U^{\scriptscriptstyle\text{tot}} - A^{\scriptscriptstyle\text{TVN}}}{T}, \end{align}\tag{24}\] where the superscript \(\scriptstyle\text{TVN}\) indicates that the function arguments are \(T, V, N_1, \dots, N_n\) and \[A^{\scriptscriptstyle\text{TVN}} := \sum_{\substack{k \in \{g, l\}}} A(T, V^{k}, \mathbf{N}^{k}).\]

When the UVN flash is performed in TVN space, the total volume \(V\) and mole numbers \(N_i^{\scriptscriptstyle\text{tot}}\) can be poorly scaled: \(V\) may range from millimeters to Kilometers and \(N_i^{\scriptscriptstyle\text{tot}}\) can vary several from nearly zero to thousands, depending on the composition. To improve numerical stability, we instead employ the gas-phase component densities \(\boldsymbol{\rho}^{g} := \{\rho_1^{g}, \dots, \rho_n^{g}\}\) together with the gas-phase volumetric fraction \(\alpha^{g}\in [0,1]\) as primary variables. These variables remain naturally well-scaled across a wide range of compositions and phase splits, thereby improving the robustness of the flash calculations.
The gas-phase volume is related to the total volume by \[\begin{align} V^{g} = \alpha^{g}V, \end{align}\] and the mole numbers in each phase are obtained from the species densities according to \[\begin{align} N_i^{g} &= \frac{\rho_i^{g} V^{g}}{M_{i,w}}, \qquad N_i^{l} = N_i - N_i^{g}, \qquad i=1,\dots,n, \end{align}\] with the inverse relation \[\begin{align} \label{rho95i95from95N95V} \rho_i^{g} = \frac{M_{i,w} N_i^{g}}{V^{g}}. \end{align}\tag{25}\] With these variables, the entropy function can be expressed as \[\begin{align} \label{SQRhoAlpha} S^{\scriptscriptstyle{T\alpha \rho}}= \frac{U^{\scriptscriptstyle\text{tot}} - A^{\scriptscriptstyle{T\alpha \rho}}}{T}, \end{align}\tag{26}\] where the superscript \(\scriptstyle{T\alpha\rho}\) indicates that the function arguments are \(T, \alpha^{g}, \rho_1^{g}, \dots, \rho_n^{g}\) and \[A^{\scriptscriptstyle{T\alpha \rho}} := A^{g}+A^{l}.\] The phasic Helmholtz energies are defined by \[\begin{align} A^{g} := A(T, \alpha^{g}, \rho_1^{g}, \dots, \rho_n^{g}), \quad A^{l} := A(T, \alpha^{l}, \rho_1^{l}, \dots, \rho_n^{l}), \end{align}\] where \(\rho_i^{l} = (N_i - N_i^{g}) M_{i, w}/(V - V^{g})\) and \(\alpha^{l}= 1 - \alpha^{g}\). The liquid-phase densities \(\rho_i^{l}\) can be expressed in terms of gas-phase variables as

\[\begin{align} \label{rhoL95rhoG} \rho_i^{l}(\rho_i^{g}, \alpha^{g}) = \frac{M_{i,w} N_i}{(1-\alpha^{g}) V} - \frac{\alpha^{g}}{1-\alpha^{g}} \, \rho_i^{g}, \quad i=1,\dots,n. \end{align}\tag{27}\] Thus, \(A^{l}\) can be expressed entirely as a function of \((T, \alpha^{g}, \boldsymbol{\rho}^{g})\): \[A^{l} = A\Bigl(T, \alpha^{l}(\alpha^{g}), \rho_1^{l}(\rho_1^{g}, \alpha^{g}), \dots, \rho_n^{l}(\rho_n^{g}, \alpha^{g})\Bigr).\]

From thermodynamics, we know that the entropy should be maximized under the prescribed \(UV\mathbf{N}\) constraints, irrespective of the particular choice of variables used to represent the entropy function. Consequently, the critical points of the entropy defined in 26 must coincide with the thermodynamic equilibrium for given \(UV\mathbf{N}\) constraints, as will be verified in the next subsection. The critical points of the entropy function are obtained as: \[\begin{align} \label{eqn:QFunction95stationarity95points} \frac{\partial S^{\scriptscriptstyle{T\alpha \rho}}}{\partial T} = 0, \qquad \frac{\partial S^{\scriptscriptstyle{T\alpha \rho}}}{\partial \alpha^{g}} = 0, \qquad \frac{\partial S^{\scriptscriptstyle{T\alpha \rho}}}{\partial \rho_i^{g}} = 0 \quad (i=1,\dots,n), \end{align}\tag{28}\] which yields system of \(n+2\) non-linear equations.

5.2 Thermodynamic consistency check↩︎

In this section, we show that the critical points of the entropy function as defined by 28 , coincide with the thermodynamic conditions for phase equilibrium. Specifically, it suffices to show that these critical points correspond to equality of pressure, equality of chemical potential for each component across the phases, a common temperature, and recovery of the prescribed internal energy constraint. The condition of equal temperature is trivially satisfied, as a common temperature is imposed in the formulation from the outset. We recall that \[\begin{align} \label{eqn:prelim95partials} \rho_i^{g} = \frac{N_i^{g} M_{i, w}}{V^{g}} = \frac{N_i^{g} M_{i, w}}{\alpha^{g} V} \implies \frac{\partial \rho_i^{g}}{\partial N_i^{g}} = \frac{M_{i, w}}{V^{g}}, \quad \frac{\partial \rho_i^{g}}{\partial \alpha^{g}} = -\frac{1}{\alpha^{g}} \frac{N_i^{g} M_{i, w}}{\alpha^{g} V} = -\frac{\rho_i^{g}}{\alpha^{g}} . \end{align}\tag{29}\] Furthermore, the inverse relations are given by \[\begin{align} \label{eqn:rho95partials} \frac{\partial N_i^{g}}{\partial \rho_i^{g}} = \frac{V^{g}}{M_{i, w}}, \quad \frac{\partial \alpha^{g}}{\partial \rho_i^{g}} = -\frac{\alpha^{g}}{\rho_i^{g}}, \quad \frac{\partial V^{g}}{\partial \rho_i^{g}} = -\frac{V^{g}}{\rho_i^{g}}. \end{align}\tag{30}\] 1) The stationarity condition with respect to the mass densities reads: \[\begin{align} 0 &= \left(\frac{\partial S^{\scriptscriptstyle{T\alpha \rho}}(T, \alpha^{g}, \rho_1^{g},\dots,\rho_n^{g})}{\partial \rho_i^{g}}\right)_{T,\alpha^{g},\{\rho_j^{g}\}_{j\ne i}} \qquad\nonumber\\ &= -\frac{1}{T}\, \left(\frac{\partial A^{\scriptscriptstyle{T\alpha \rho}}(T,\alpha^{g},\rho_1^{g},\dots,\rho_n^{l})}{\partial \rho_i^{g}}\right)_{T,\alpha^{g},\{\rho_j^{g}\}_{j\ne i}} \label{eqn:partial95Q95partial95rho95i}. \end{align}\tag{31}\] We therefore focus on the derivative of the total Helmholtz free energy with respect to \(\rho_i^{g}\). Using \(A^{\scriptscriptstyle{T\alpha \rho}}=A^{g}+A^{l}\), we obtain \[\begin{align} \label{eqn:partial95A95tot95rho95gas95tmp} \left(\frac{\partial A^{\scriptscriptstyle{T\alpha \rho}}}{\partial \rho_i^{g}}\right)_{T,\alpha^{g},\{\rho_j^{g}\}_{j\ne i}} &= \left(\frac{\partial A^{g}}{\partial \rho_i^{g}}\right)_{T,\alpha^{g},\{\rho_j^{g}\}_{j\ne i}} + \left(\frac{\partial A^{l}}{\partial \rho_i^{g}}\right)_{T,\alpha^{g},\{\rho_j^{g}\}_{j\ne i}}. \end{align}\tag{32}\] As \(\rho_i^{g}\) is a function of \(N_i^{g}\) and \(V^{g}\) (see 29 ); using the chain rule gives \[\begin{align} \left(\frac{\partial A^{g}}{\partial \rho_i^{g}}\right)_{} &= \frac{\partial V^{g}}{\partial \rho_i^{g}} \left(\frac{\partial A^{g}}{\partial V^{g}}\right)_{T,\{N_k^{g}\}} + \sum_k \frac{\partial N_k^{g}}{\partial \rho_i^{g}} \left(\frac{\partial A^{g}}{\partial N_k^{g}}\right)_{T,V^{g},\{N_{m}^{g}\}_{m\ne k}}, \tag{33} \\ \left(\frac{\partial A^{l}}{\partial \rho_i^{g}}\right)_{} &= \frac{\partial V^{l}}{\partial \rho_i^{g}} \left(\frac{\partial A^{l}}{\partial V^{l}}\right)_{T,\{N_k^{l}\}} + \sum_k \frac{\partial N_k^{l}}{\partial \rho_i^{g}} \left(\frac{\partial A^{l}}{\partial N_k^{l}}\right)_{T,V^{l},\{N_{m}^{l}\}_{m\ne k}}. \tag{34} \end{align}\] In equations 33 and 34 , only terms inside the summation where \(k=i\) contribute. Using the chemical potential definition, \(\mu_i^{g}=\left(\frac{\partial A^{g}}{\partial N_i^{g}}\right)_{T,V^{g},\{N_{j}^{g}\}_{j\ne i}}\) and rearranging, we get

\[\begin{align} \label{eqn:partial95A95rho95gas95tmp} \left(\frac{\partial A^{\scriptscriptstyle{T\alpha \rho}}}{\partial \rho_i^{g}}\right)_{} &= \mu_i^{g}\,\frac{\partial N_i^{g}}{\partial \rho_i^{g}} + \mu_i^{l}\,\frac{\partial N_i^{l}}{\partial \rho_i^{g}} + \frac{\partial V^{g}}{\partial \rho_i^{g}}\underbrace{\left(\frac{\partial A^{g}}{\partial V^{g}}\right)_{T, \{N_k^{g}\}}}_{-p^{g}} + \frac{\partial V^{l}}{\partial \rho_i^{g}}\underbrace{\left(\frac{\partial A^{l}}{\partial V^{l}}\right)_{T, \{N_k^{l}\}}}_{-p^{l}}. \end{align}\tag{35}\] Since \(V^{l}=V - V^{g}\), it follows that \(\frac{\partial V^{l}}{\partial \rho_i^{g}} = -\frac{\partial V^{g}}{\partial \rho_i^{g}}\). From 29 , we have \[\frac{\partial V^{g}}{\partial \rho_i^{g}} = -V^{g}/\rho_i^{g}, \quad \frac{\partial N_i^{g}}{\partial \rho_i^{g}}=V^{g}/M_{i,w}.\] Moreover, using mass balance relation, \(N_i^{l}=N_i-N_i^{g}\), we obtain \(\frac{\partial N_i^{l}}{\partial \rho_i^{g}}=-V^{g}/M_{i,w}\). Substituting these expressions in 35 and using 32 gives \[\begin{align} \left(\frac{\partial A^{\scriptscriptstyle{T\alpha \rho}}}{\partial \rho_i^{g}}\right)_{} &= \mu_i^{g}\frac{V^{g}}{M_{i,w}} - \mu_i^{l}\frac{V^{g}}{M_{i,w}} + \left(p^{g} - p^{l}\right) \frac{V^{g}}{\rho_i^{g}}, \nonumber\\ &= \left(\mu_i^{g}-\mu_i^{l}\right)\frac{V^{g}}{M_{i,w}} + \left(p^{g} - p^{l}\right) \frac{V^{g}}{\rho_i^{g}}. \end{align}\] The stationarity condition 31 then yields: \[\begin{align} \label{eqn:partial95Q95partial95rho95i952} 0 = (\mu_i^{g}-\mu_i^{l})\frac{\rho_i^{g}}{M_{i,w}} + (p^{g} - p^{l}), \qquad \forall i \in \{1, \dots, n\}. \end{align}\tag{36}\] Note that 36 is system of \(n\) equations, which in expanded form reads: \[\label{eqn:partial95A95partial95density95final} \begin{align} (\mu_1^{g}-\mu_1^{l})\frac{\rho_1^{g}}{M_{1,w}} &+ (p^{g} - p^{l}) = 0, \\ &\vdots \\ (\mu_n^{g}-\mu_n^{l})\frac{\rho_n^{g}}{M_{n,w}} &+ (p^{g} - p^{l}) = 0. \end{align}\tag{37}\] 2) The stationarity condition with respect to the gas phase fraction \(\alpha^{g}\) reads: \[\begin{align} \left(\frac{\partial A^{\scriptscriptstyle{T\alpha \rho}}(T, \alpha^{g}, \rho_1^{g}, \dots, \rho_n^{g})}{\partial \alpha^{g}}\right)_{T, \rho_1^{g}, \dots, \rho_n^{g}} &= 0. \end{align}\] Substitute \(A^{\scriptscriptstyle{T\alpha \rho}} = A^{g} + A^{l}\) and differentiate with respect to \(\alpha^{g}\), we get \[\begin{align} \frac{\partial A^{\scriptscriptstyle{T\alpha \rho}}}{\partial \alpha^{g}} &= \frac{\partial A^{g}}{\partial \alpha^{g}} + \frac{\partial A^{l}}{\partial \alpha^{g}}. \end{align}\] Since \(\alpha^{g}= V^{g} / V\), the partial derivative with respect to \(\alpha^{g}\) is \[\begin{align} \frac{\partial A^{g}}{\partial \alpha^{g}} = \frac{\partial A^{g}}{\partial V^{g}} \frac{\partial V^{g}}{\partial \alpha^{g}} = \frac{1}{V} \frac{\partial A^{g}}{\partial V^{g}}. \end{align}\] For the liquid phase, using \(V^{l} = V - V^{g}\), we obtain \[\begin{align} \frac{\partial A^{l}}{\partial \alpha^{g}} &= \frac{\partial A^{l}}{\partial V^{l}} \frac{\partial V^{l}}{\partial \alpha^{g}} = -\frac{1}{V} \frac{\partial A^{l}}{\partial V^{l}} . \end{align}\] Combining the two and recalling \(\frac{\partial A}{\partial V} = -p\), it follows that \[\begin{align} \frac{\partial A^{\scriptscriptstyle{T\alpha \rho}}}{\partial \alpha^{g}} &= \frac{p^{l} - p^{g}}{V}. \end{align}\] At the stationarity point \(\frac{\partial A^{\scriptscriptstyle{T\alpha \rho}}}{\partial \alpha^{g}} = 0\), which implies \[\begin{align} \label{eqn:mech95eqm} \boxed{ p^{g} = p^{l}}, \end{align}\tag{38}\] corresponding to the condition of mechanical equilibrium.

Substituting \(p^{g} = p^{l}\) in 37 yields \[\begin{align} \label{eqn:chem95eqm} \boxed{ \mu_i^{g} = \mu_i^{l}, \quad \forall \in \{1, \dots, n\}}, \end{align}\tag{39}\] which is the condition of chemical equilibrium.

3) Stationarity with respect to temperature \(T\) reads: \[\begin{align} \left(\frac{\partial S^{\scriptscriptstyle{T\alpha \rho}}(T, \alpha^{g}, \rho_1^{g}, \dots, \rho_n^{g})}{\partial T}\right)_{\alpha^{g}, \rho_1^{g}, \dots, \rho_n^{g}} &= 0. \end{align}\] Recall that \(S^{\scriptscriptstyle{T\alpha \rho}}= (U^{\scriptscriptstyle\text{tot}} - A^{\scriptscriptstyle{T\alpha \rho}})/T\). Differentiating with respect to \(T\) gives: \[\begin{align} \label{eqn:partial95S95partial95T95tmp} \frac{\partial S^{\scriptscriptstyle{T\alpha \rho}}}{\partial T} &= -\frac{U^\scriptscriptstyle\text{tot}}{T^2} - \frac{1}{T} \frac{\partial A^{\scriptscriptstyle{T\alpha \rho}}}{\partial T} + \frac{A^{\scriptscriptstyle{T\alpha \rho}}}{T^2}. \end{align}\tag{40}\] Substituting \(A^{\scriptscriptstyle{T\alpha \rho}} = A^{g} + A^{l}\) and differentiating with respect to \(T\) yields: \[\begin{align} \frac{\partial A^{\scriptscriptstyle{T\alpha \rho}}}{\partial T} &= \frac{\partial A^{g}}{\partial T} + \frac{\partial A^{l}}{\partial T}. \end{align}\] Using thermodynamic identity \(\frac{\partial A}{\partial T} = -S,\) we obtain \[\begin{align} \frac{\partial A^{\scriptscriptstyle{T\alpha \rho}}}{\partial T} &= -S^{g} - S^{l}. \end{align}\] Substituting this in 40 , we obtain \[\begin{align} \frac{\partial S^{\scriptscriptstyle{T\alpha \rho}}}{\partial T} &= -\frac{U^\scriptscriptstyle\text{tot}}{T^2} + \frac{S^{g} + S^{l}}{T} + \frac{A^{g} + A^{l}}{T^2}, \\ &= \frac{A^{g} + TS^{g} + A^{l} + TS^{l} - U^\scriptscriptstyle\text{tot}}{T^2}. \end{align}\] Applying the relation \(U = A + TS\) for each phase yields \[\begin{align} \frac{\partial S^{\scriptscriptstyle{T\alpha \rho}}}{\partial T} = \frac{U^{g} + U^{l} - U^\scriptscriptstyle\text{tot}}{T^2}. \end{align}\] The stationary condition \(\frac{\partial S^{\scriptscriptstyle{T\alpha \rho}}}{\partial T} = 0\) therefore gives \[\begin{align} \label{eqn:constraint95recovery} \boxed{ U^\scriptscriptstyle\text{tot}= U^{g} + U^{l} }, \end{align}\tag{41}\] which is the desired constraint on the internal energy. Equations 38 , 39 , and 41 together constitute the conditions of thermodynamic equilibrium under UVN specifications.

6 Temporal Discretization↩︎

We consider the differential-algebraic system (DAE) that arises from the pipe flow equations 6 and tank model 9 along with the flash problem (28 ). With the state variables partitioned into differential (conservative) variables \(\mathbf{U}\) and algebraic (non-conservative) variables \(\mathbf{V}\), the resulting DAE system can be expressed as \[\begin{align} \frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t} &= \mathbf{f}(\mathbf{U},\mathbf{V}), \tag{42} \\ 0 &= g(\mathbf{U}, \mathbf{V}) \tag{43}. \end{align}\] The precise form of \(\mathbf{U}\) and \(\mathbf{V}\) for the pipe and tank model is discussed below.

Pipe flow↩︎

For spatially discretized pipe flow, we define: \[\begin{align} &\mathbf{U}= [\mathbf{U}_{1}, \ldots, \mathbf{U}_{N}]^T, \qquad \mathbf{U}_{i} = [\rho_i, (\rho u)_i, (\rho E)_i]^T, \\ &\mathbf{V}= [\mathbf{V}_{1}, \ldots, \mathbf{V}_{N}]^T, \qquad \mathbf{V}_{i} = [\rho^{g}_{k,i}, \alpha^{g}_{i}, T_{i}]^T, \end{align}\] where \(\rho^{g}_{k}\), \(\alpha^{g}\), and \(T\) are the partial mass density of component \(k\) in gas phase, volumetric phase fraction of gas, and temperature, respectively. The subscript \(i\) refers to the the \(i^{\text{th}}\) cell. A first-order spatial discretization of the conservation laws yields the semi-discrete form: \[\mathbf{f}_i(\mathbf{U}, \mathbf{V}) := -\frac{1}{\Delta x} \left( \Hat{{\mathcal{F}}}_{i+\frac{1}{2}}(\mathbf{U}_{i+1}, \mathbf{U}_i, \mathbf{V}_{i+1}, \mathbf{V}_i) - \Hat{{\mathcal{F}}}_{i-\frac{1}{2}}(\mathbf{U}_i, \mathbf{U}_{i-1}, \mathbf{V}_i, \mathbf{V}_{i-1}) \right),\] where \(\Hat{{\mathcal{F}}}_{i \pm \frac{1}{2}}\) denotes the numerical flux function, and \(\mathbf{f}_i\) represents the discretized right-hand side of the system 1 .

Tank model↩︎

For the tank, the differential variables are defined as: \[\mathbf{U}= [N_1, \dots, N_n, U]^T,\] where \(N_i\) denotes total number of moles of component \(i\) in the mixture and \(U\) is the total internal energy. The algebraic variables \(\mathbf{V}\) are defined as: \[\mathbf{V}= [N_1^{g}, \dots, N_n^{g}, V^{g}, T]^T,\] where \(N_i^{g}\) is the number of moles of component \(i\) in the gas phase, \(V^{g}\) is the volume occupied by the gas phase, and \(T\) is the temperature. The function \(\mathbf{f}\) for the tank is defined by 10 .

Time integration scheme

To integrate the DAE system 42 in time, we employ a half-explicit forward Euler scheme [47]. This method consists of an explicit update for the differential variables, followed by a nonlinear solve to enforce the algebraic constraints. First, the differential update is computed as: \[\mathbf{U}^{n+1} = \mathbf{U}^n + \Delta t \, \mathbf{f}(\mathbf{U}^n, \mathbf{V}^n),\] followed by the solution of the nonlinear constraint: \[0 = g(\mathbf{U}^{n+1}, \mathbf{V}^{n+1}),\] to determine \(\mathbf{V}^{n+1}\). Solving this constraint corresponds to a flash calculation as formulated by system of equations 28 , in 5.

7 Numerical Experiments↩︎

In this section, we discuss the results. First, we discuss the results for the tank model, followed by the results for the pipeline depressurization. We use PR-EOS for all the results; see 12 for details. The diagonal elements of the binary interaction parameters (BIP) matrix of PR are set to zero, and the off-diagonal terms are specified in the tables 1 and 3.

7.1 Tank model↩︎

We validate our approach on the two problems from Castier [33]. The specifications of the problem are summarized in 1

7.1.1 Problem 1: Light components↩︎

In this problem, we consider a mixture of methane() and hydrogen sulfide (). The fluid in the tank is in a single-phase state at \(t=0s\). There is one input stream feeding the tank, which is in a two-phase state. Since the stream conditions are given in terms of temperature and pressure, we need to perform a PTN-flash [17] only once at the beginning of the simulation to determine the phase split and hence, the enthalpy of individual phases. Total enthalpy is obtained by adding up these individual phasic contributions. The results are shown in 2. An excellent match with the results reported in the Castier [33] is obtained. Initially, there is a sharp drop in the temperature reaching a minimum temperature (\(\approx 260K\)) around \(680s\) after which the temperature starts to rise and there is an onset of the second phase around \(1792s\). We see a change in the temperature slope at this point. Pressure, on the other hand, keeps on increasing linearly as more mass is being added to the tank.

7.1.2 Problem 2: loading↩︎

In this example, we have a mixture of medium-heavy hydrocarbons. Initially, there is a very small amount of in the tank. The tank is being fed with . For a given pressure, temperature, and feed composition, we can solve the PR–EOS for volume and thereby obtain the enthalpy flow rate. Initially, the fluid in tank is in two-phase state. Throughout the process, the heat is removed at a constant rate of 100 J/s. 3 shows the results obtained in the current work and compares them with those by Casiter [33]. Again, a very good match can be observed. Until around \(341s\), there are two phases in the tank; beyond this point, only a single phase remains. The temperature then becomes relatively constant, but the pressure exhibits a sharp increase, indicating that the fluid has transitioned to a liquid state, which behaves almost as an incompressible fluid.

To summarize, for both tank depressurization problems, our approach was able to handle the transition from single-phase to two-phase and vice versa without any difficulty. We now turn to a more challenging scenario and apply the approach to pipeline depressurization.

7.2 Pipeline Depressurization↩︎

Having validated our flash solver for transient simulation for the tank model, we now apply it to pipeline transport of a two-phase multicomponent fluid. We consider six different mixtures: 4 two-component mixtures; a five-component mixture; and a single-component fluid. Such numerical experiments have been previously considered in previous studies; see references [27], [38], [39] for single-component fluids and Munkejord et al. [3] for multicomponent mixtures.

In Munkejord et al. [3], an outflow boundary condition is prescribed at the open end(right) of the pipe. In contrast, in our simulations, we avoid imposing this boundary condition directly. Instead, we extend the computational domain to twice the original pipe length and formulate the problem as a shock tube with the initial discontinuity located at the midpoint. The reason for adopting this approach is that the flow conditions at the open end will be choked, and hence, the information from outside the pipe cannot propagate upstream into the computational domain. The simulation is terminated before the fastest traveling wave reaches the end of the pipe.

In all simulations, a grid of \(800\) cells was used (i.e., \(400\) cells per side), whereas Munkejord et al. [3] used \(2400\) cells. Furthermore, we employed a CFL value of \(0.9\), while Munkejord et al. [3] used a CFL of \(0.85\). Due to the presence of initial discontinuities, an initial time-step of \(1 \times 10^{-12}\) is used to begin the simulation. Subsequent time steps are determined based on the fastest traveling waves:

\[\begin{align} \label{simulation46dt} \Delta t = \frac{\Delta x}{\max_{i} (|u_i \pm a_i|)}, \end{align}\tag{44}\] where subscript \(i\) corresponds to the \(i^{\mathrm{th}}\) cell and \(a_i\) represents the speed of sound. The speed of sound is computed using Wood’s formula; further details are provided in 11.

The initial conditions are provided in 2 and the mixture composition is given in 3. To obtain the density and the number of moles, the compressibility factor \(Z\) is first computed. The total number of moles \(N\) is then calculated using the relation \[Z = \frac{PV}{NRT},\] where \(V\) is the total volume of the computational cell, \(R\) is the universal gas constant. The mole vector is subsequently obtained by multiplying \(N\) by the composition vector specified in the 3.

7.3 Validation↩︎

The spatial convergence is discussed in 9. Here, we focus on the validation with the results reported by Munkejord et al. [3] for –rich mixtures, with specifications provided in 2 and 3. It is important to highlight the differences between the setup considered here and that of Munkejord et al. [3] as these distinctions may contribute to the discrepancies in the results. The main differences are: (i) we consider a shock-tube setup, whereas Munkejord et al.imposed a boundary condition, and (ii) the Peng–Robinson (PR) coefficients may differ, since they are not explicitly reported in Munkejord et al. Despite the differences in modeling choices, our results show good agreement with the literature, as can be seen in 4, where the HEM trajectory in the pressure–temperature (PT) plane is shown for the binary and the five-component mixtures. Notably, since we did not take into account the heat transfer and friction, the simulation path in the left section of the pipe corresponds to an isentropic process.

5 presents the results for pressure and temperature along the pipe length for a binary mixture of and at \(t=0.1s\). Around \(x \approx 120m\), a slight change in the slope in the temperature profile is observed as highlighted by an orange rectangle. This phenomenon is intrinsic to multicomponent mixtures and absent in single-component systems, as illustrated in the 10 where the evaporation wave ends abruptly into the contact wave. This difference can be explained by considering the topology of the phase diagram: for pure substances, the saturation points form a single curve in PT-space; see 6 (a), where each pressure–temperature point can correspond to multiple values of phase fraction. For multicomponent mixtures, however, the two-phase region forms an envelope bounded by bubble and dew lines; see 6 (b). Within this envelope, the iso-quality lines (where phase fraction remains constant) span the envelope from the bubble curve (vapor fraction \(0\)) to the dew curve (vapor fraction \(1\)). Consequently, each pressure–temperature condition inside the envelope corresponds to a unique value of vapor fraction.

7 (a) presents the velocity profile. Note that the fluid velocity changes across the contact wave, unlike the single-phase case. This has also been observed in [39], [48]. 7 (b) shows the spatial variation of the speed of sound along the pipeline. The speed of sound is an important component for calculations of wave speeds in approximate Riemann solvers (HLLC in this case); see 10. A sharp drop in sound speed can be observed in the two-phase region when compared to the single-phase regions. This phenomenon is attributed to the increased compressibility of the fluid when both liquid and gaseous phases are present, and the gas bubbles dampen pressure disturbances.

The results for the five-component mixture are presented in 8. The overall wave structure is very similar to that of the binary case. The slight differences in wave amplitudes and slopes (highlighted with an orange rectangle) can be observed and arise due to differences in the initial conditions, composition, and the EOS coefficients. 9 (a) provides a magnified view of the sloping region. Notably, as is shown in 8 (b), the pressure does not change across this region. Moreover, as illustrated in 9 (b), this region corresponds to a transition from the two-phase to the single-phase regime. Finally, 10 presents the spatial pressure temperature (10 (a)) and (10 (b)) profiles for the single-component case.

To summarize, the proposed unified framework has been tested on multiple benchmark problems of varying complexity. It consistently shows a very good agreement with results from the literature, thereby demonstrating its robustness. It is important to emphasize that stability analysis and flash calculations are computationally demanding, and their corresponding routines are therefore natural candidates for future optimization.

7.4 Wave structure↩︎

For all three test cases, the solution exhibits a wave structure consistent with the typical wave patterns of a two-phase Riemann problem:

1. Left-moving rarefaction wave: This smooth wave corresponds to the expansion of the fluid in the pipeline, leading to a drop in pressure and temperature across the wave. It is a genuinely non-linear wave in the context of the Riemann Problem [49].

2. Evaporation wave: This is also a smooth wave, across which liquid boiling occurs. Notably, this wave is absent in the Riemann problem for a single-phase case. Furthermore, for

   1. Pure component: The evaporation wave ends abruptly and transitions into the single-phase region with a contact discontinuity.

   2. Multicomponent: The evaporation wave exhibits a sloping region that corresponds to varying values of vapor phase fractions in the phase envelope in \(PT\) space. Thereafter, it transitions into the single-phase region with a contact discontinuity.

3. Right-moving contact discontinuity: This wave follows the shock wave and marks material separation. No fluid mixing occurs across this wave. Unlike in single-phase fluids, it is accompanied by a change in fluid velocity.

4. Right-moving shock wave: This wave corresponds to fluid compression, is accompanied by entropy production, and results in an increase in temperature and pressure.

8 Conclusion↩︎

In this paper, we have developed a unified framework that integrates fluid dynamics and thermodynamics for simulating the transport of multicomponent fluid in two-phase conditions. The –rich mixtures relevant to CCS applications were considered as test cases. We adopted the homogeneous equilibrium model (HEM) to describe the flow, and employed a Helmholtz free energy-based equation of state to model the thermodynamic properties. The key contribution is the development and testing of a tailored UVN-flash framework for dynamic pipeline transport models. The reformulation of the UVN-flash relies on the introduction of a new set of variables, namely, \(\rho^{g}_1, \dots, \rho^{g}_n, \alpha, T\). This choice aligns seamlessly with the inputs from the fluid dynamics solver, yielding better-scaled variables when compared to the standard UVN/TVN variables. Furthermore, we demonstrated that the critical points of the entropy function, when expressed in these variables, correspond to the classical thermodynamic equilibrium conditions. Additionally, we discussed a formulation of stability analysis for the UVN-flash which allows us to reliably detect single- and two-phase states. This is a crucial component for robust multiphase transient pipeline simulations.

The methodology was tested on a set of depressurization scenarios, including two tank problems (a binary and a five-component mixture) as well as pipeline cases for pure , four binary mixtures, and one five-component mixture. Numerical experiments demonstrate that the proposed framework provides a consistent and accurate description of the coupled thermodynamics and fluid dynamics problem, and is capable of capturing the extreme conditions relevant for safe pipeline design.

A flowchart of the algorithm and details of the temporal and spatial discretization of the resulting DAE system were also provided. Since thermodynamics computations are computationally intensive, exploring strategies for speeding up these calculations is an interesting research direction. Another promising research direction could be to investigate the application of SVN-flash in pipeline depressurization, as it could offer an efficient and accurate method for modeling the thermodynamics of isentropic flow.

CRediT author statement↩︎

Pardeep Kumar: Conceptualization, Methodology, Software, Validation, Analysis, Writing – Original Draft

Patricio I. Rosen Esquivel: Project Administration, Funding Acquisition, Conceptualization, Supervision, Writing - Review & Editing

Declaration of Generative AI and AI-assisted technologies in the writing process↩︎

During the preparation of this work the authors used GitHub Copilot in order to propose wordings and mathematical typesetting. After using this tool/service, the authors reviewed and edited the content as needed. The authors take full responsibility for the content of the publication.

Declaration of competing interest↩︎

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments↩︎

This research was generously supported by Shell Projects and Technology, and we deeply appreciate their invaluable contribution.

9 Spatial Convergence↩︎

To evaluate the spatial convergence in the pipeline simulations, we consider the binary mixture case. A high-resolution run on a mesh of 6400 cells is used as the reference solution. The \(L_1\) error in temperature is reported in 11, where the results are plotted on a log-log scale. We observe increased accuracy as the mesh is refined. The observed convergence rate lies between first and second-order accuracy.

10 Spatial discretization of HEM model↩︎

To compute the numerical flux \(\hat{{\mathcal{F}}}_{i+\frac{1}{2}}\), we employ the HLLC (Harten—Lax—van Leer—Contact) scheme which is an approximate Riemann solver [50]. This method exploits the wave structure of the Riemann problem to compute fluxes at the interfaces between adjacent computational cells. For single-phase compressible flow, the Riemann solution typically consists of three distinct wave types: rarefaction waves, shock waves, and a contact discontinuity. In the context of two-phase flow, however, an additional discontinuity, referred to as the evaporation wave, may arise due to phase transition phenomena. The HLLC scheme approximates the fluxes associated with each of these waves by enforcing the Rankine-Hugoniot condition. A complete treatment of the HLLC scheme for Euler’s equations can be found in Toro [50]. Below, we present the final expressions for the numerical fluxes at the interface \(i+\frac{1}{2}\):

\[\begin{align} \hat{{\mathcal{F}}}_{i+1/2} &= \begin{cases} F_L, & \text{if } 0 \leq S_L, \\ F^{\star}_{L}, & \text{if } S_L \leq 0 \leq S^{\star}, \\ F^{\star}_{R}, & \text{if } S^{\star} \leq 0 \leq S_R, \\ F_R, & \text{if } 0 \geq S_R . \end{cases} \\ F^{\star}_{L} &= F_L + S_L (U^{\star}_{L} - U_L), \\ F^{\star}_{R} &= F_R + S_R (U^{\star}_{R} - U_R), \\ S_L &= \min(u_L - a_L, u_R - a_R), \\ S_R &= \max(u_L + a_L, u_R + a_R), \\ S^{\star} &= \frac{p_R - p_L + \rho_L u_L (S_L - u_L) - \rho_R u_R (S_R - u_R)}{\rho_L (S_L - u_L) - \rho_R (S_R - u_R)}, \\ U^{\star}_{K} &= \rho_K \left( \frac{S_K - u_K}{S_K - S^{\star}} \right) \begin{bmatrix} 1 \\ S^{\star} \\ E_K + (S^{\star} - u_K) \left [ S^{\star} + \frac{p_K}{\rho_K(S_K - u_K)}\right] \end{bmatrix}. \end{align}\] Here, \(K = L\)(left state) or \(K = R\)(right state), \(U_K \approx {\mathcal{U}}(x_{K},t)\), \(F_{K} = {\mathcal{F}}(U_{K}), a_{K}\) denotes the local speed of sound computed using the equation of state [51]. The wave speeds \(S_L\) and \(S_R\) represent the velocities of the left and right-propagating waves, respectively, and \(S^{\star}\) is the speed of the contact wave. The accuracy of the HLLC scheme is strongly influenced by the quality of the wave speeds(\(S_L\) and \(S_R\)) estimates. Here, we have presented one such method for estimating these wave speeds. For a broader discussion of wave speed approximation, the interested reader is referred to Toro [50].

11 Wood’s Speed of Sound for an \(n\)-Component, Two-Phase Mixture↩︎

Consider an \(n\)-component mixture in mechanical equilibrium between a gas (\(g\)) and a liquid (\(l\)) phase. Let the vapor volume fraction be \(\alpha^{g}\in [0,1]\). Denote the sound speeds by \(a_g, a_l\) and mass densities by \(\rho_g, \rho_l\) of the vapor and liquid phase, respectively. Wood’s mixture relation [52] then reads \[\frac{1}{\rho_m a_m^2} = \frac{\alpha^{g}}{\rho_ga_g^2} + \frac{1-\alpha^{g}}{\rho_la_l^2}, \qquad \rho_m = \alpha^{g}\rho_g+ (1-\alpha^{g})\rho_l, \label{eq:woods-mixture}\tag{45}\] where \(\rho_m\) is the overall mixture density.

Phase densities from \(N\) and \(V\)↩︎

Let \(N_i^{g}\) and \(N_i^{l}\) be the moles of component \(i\) in vapor and liquid, respectively, and \(M_{w,i}\) the molar mass. With phase volumes \(V^{g}\) and \(V^{l}\), the phasic component mass densities can be written as: \[\begin{align} \rho_i^{g} &= \frac{N_i^{g} M_{w,i}}{V^{g}}, \quad \rho_g= \sum_{i=1}^n \rho_i^{g}, \\ \rho_i^{l} &= \frac{N_i^{l} M_{w,i}}{V^{l}}, \quad \rho_l= \sum_{i=1}^n \rho_i^{l}. \end{align}\]

Phasic speed of sound↩︎

Using the thermodynamic identity for the isentropic speed of sound for each phase \[a^2 = \left( \frac{\partial p}{\partial \rho} \right)_{S} = \frac{V}{\rho} \left[ -\frac{\partial p}{\partial V} + \frac{T}{C_v} \left( \frac{\partial p}{\partial T} \right)^2 \right], \label{eq:speed-of-sound}\tag{46}\] where \(C_v\) is the isochoric heat-capacity and \(\rho\) is the phasic density (\(\rho_{g/l}\)).

12 Peng–Robinson Equation of State↩︎

We employ the Peng–Robinson equation of state (EOS) [42], which is formulated as follows: \[P(T, V, N_1, \ldots, N_n) = \frac{NRT}{V - B} - \frac{a(T)N^2}{V^2 + 2BV - B^2}, \label{eq:peng95robinson}\tag{47}\] where \(T\) is the temperature, \(V\) is the volume, \(N_i\) represents the number of moles of component \(i\) in the system, \(R\) is the universal gas constant and \(N\) is the total number of moles in the system. The parameter \(a(T)\) characterizes the attractive intermolecular forces, while the effective co-volume \(B\) is defined in terms of the parameter \(b\), which characterizes the repulsive interactions, as:

\[\begin{equation} \tag{48} B = b N, \end{equation} \begin{equation} a = \sum_{i=1}^{n} \sum_{j=1}^{n} x_i x_j a_{ij}, \tag{49} \end{equation} \begin{equation} a_{ij} = (1 - \delta_{ij}) \sqrt{a_i a_j}, \tag{50} \end{equation} \begin{equation} a_i(T) = 0.45724 \frac{R^2 T_{\text{crit},i}^2}{P_{\text{crit},i}} \left[1 + m_i \left(1 - \sqrt{T_{r,i}} \right)\right]^2, \tag{51} \end{equation} \begin{equation} b = \sum_{i=1}^{n} x_i b_i, \tag{52} \end{equation} \begin{equation} b_i = 0.0778 \frac{R T_{\text{crit},i}}{P_{\text{crit},i}}, \tag{53} \end{equation}\] where \(x_i = N_i/N\) is the mole fraction of component \(i\), \(T_{\text{crit},i}\), \(P_{\text{crit},i}\) and \(T_{r, i} = T/T_{\text{crit},i}\) are the critical temperature, critical pressure and the reduced temperature of component \(i\), and \(\delta_{ij}\) is the binary interaction parameter between component \(i\) and \(j\). The parameter \(m_i\) accounts for the acentric factor \(\omega_i\) as: \[m_i = \begin{cases} 0.37464 + 1.54226 \omega_i - 0.26992 \omega_i^2, & \omega_i < 0.5, \\ 0.3796 + 1.485 \omega_i - 0.1644 \omega_i^2 + 0.01667 \omega_i^3, & \omega_i \geq 0.5. \end{cases} \label{eq:m95i}\tag{54}\]
The residual internal energy, \(U\), in the context of the Peng-Robinson EOS is expressed as follows. \[\begin{align} U(T, V, N_1, \ldots, N_n) &= N \frac{T \partial_T(a) - a}{2 \sqrt{2}b} \ln \left[ \frac{V + \delta_1 B}{V + \delta_2 B} \right] \nonumber \\ &- NR(T - T_0) + \sum_{i=1}^{n} N_i \int_{T_0}^{T} c_{p,i}^{\mathrm{ig}}(\xi) \, d\xi + N u_0, \label{eq:internal95energy} \end{align}\tag{55}\] where \(\partial_T(a)\) is the temperature derivative of \(a(T)\), \(T_0\) is a reference temperature, \(\alpha_{ik}\) are empirical constants, \(\delta_1 = 1 + \sqrt{2}\) and \(\delta_2 = 1-\sqrt{2}\). The residual entropy, \(S\), is given as \[\begin{align} S(T, V, N_1, \ldots, N_n) &= NR \ln \left[ \frac{V - B}{V} \right] + N\frac{\partial_T(a)}{2 \sqrt{2}b} \ln \left[ \frac{V + \delta_1 B}{V + \delta_2 B} \right] \nonumber \\ &+ R \sum_{i=1}^{n} N_i \ln \frac{V P_0}{N_i RT} + \sum_{i=1}^{n} N_i \int_{T_0}^{T} \frac{c_{p,i}^{\text{ig}}(\xi)}{\xi} d\xi, \label{eq:entropy} \end{align}\tag{56}\] where \(c_{p,i}^{\text{ig}}(T)\) is the ideal gas heat capacity of component \(i\) and \(P_0\) is a reference pressure. The heat capacity \(c_{p,i}^{\text{ig}}(T)\) can be written as: \[c_{p,i}^{\text{ig}}(T) = \sum_{k=0}^{3} \alpha_{ik} T^k. \label{eq:heat95capacity}\tag{57}\] Now, we can simplify the integral in 55 as \[\int_{T_0}^{T} c_{p,i}^{\mathrm{ig}}(\xi) \, d\xi = \sum_{k=0}^{3} \alpha_{ik} \frac{T^{k+1} - T_0^{k+1}}{k+1},\] and the integral in 56 as \[\int_{T_0}^{T} \frac{c_{p,i}^{\text{ig}}(\xi)}{\xi} d\xi, = \alpha_{i0} \ln \left( \frac{T}{T_0} \right) + \sum_{k=1}^{3} \alpha_{ik} \frac{T^{k} - T_0^{k}}{k}.\] The coefficients \(\alpha_0, \alpha_1, \alpha_2, \alpha_3\) for the fluids considered in this work are listed in Table 4, while the parameters of the Peng–Robinson equation of state are summarized in Table 5. It is important to note that the arguments of logarithmic terms must remain positive in Equations 55 and 56 . If this condition is violated, the current step should be rejected or appropriately truncated to maintain physical consistency. The reference state is specified at \(T_0 = 298.15\,\mathrm{K}\) and \(P_0 = 1\,\mathrm{bar}\), where the molar internal energy is defined as \[u_0 = u(T_0, P_0) = h(T_0, P_0) - RT_0 = -RT_0 = -2478.95687512\,\mathrm{J\,mol^{-1}}.\] This definition ensures that the molar enthalpy of the ideal gas at the reference conditions is zero [42], i.e., \(h(T_0, P_0) = 0\). Furthermore, the molar entropy of each pure component as an ideal gas is also set to zero at this state, \(s_i^{\text{ideal}}(T_0, P_0) = 0\).

References↩︎