Keywords

1 Introduction

Two-phase and multi-phase flows are important and common phenomena in petroleum industry, where oil, gas and water are often produced and transported together. In particular, engineers and researchers in reservoir engineering study drainage problems arising during the development and production of oil and gas reservoirs so as to obtain a high economic recovery, by developing, conducting, and interpolating the simulation of subsurface flows of reservoir fluids, including water, hydrocarbon, \(\mathrm{CO}_2\), \(\mathrm H_2\mathrm{S}\) for example in porous geological formation. Field-scale (Darcy-scale) simulation has conventionally and routinely used for this purpose [2, 4, 12, 16, 17]. A number of parameters like relative permeability and capillary pressure are taken as given functions in Darcy-scale simulation [14]. To study these parameters as well as to obtain deep understanding of porous media flow and transport, researchers develop and utilize pore-scale simulation of two-phase [1, 3, 5,6,7,8,9, 11, 13], which has been shown to be a great research tool to understand the complex hydrodynamic and phase behaviors of the systems.

2 Mathematical Modeling Framework

2.1 Mathematical Model of Bulk Properties

In this paper we consider the Van der Waals equation of state (EOS) and the Peng-Robinson EOS [15] to model the bulk properties of our fluid system. We note that the Van der Waals EOS is one of the simplest equations of state that allow phase splitting. On the other hand, the Peng-Robinson EOS is the most popular equation of state (EOS) to model and compute the fluid equilibrium property of hydrocarbon fluid and other petroleum fluids, and it is widely used in reservoir engineering and oil industries. Some of the material in this section can be found in many textbooks, but we list them here for completeness of information only.

The PVT-form of the Van der Waals EOS appears as

$$\begin{aligned} \left( p + \frac{a}{v^2}\right) \left( v-b\right) = RT. \end{aligned}$$
(1)

At the critical condition, we have \(\left( \frac{\partial p}{\partial v}\right) _T=0\) and \(\left( \frac{\partial ^2 p}{\partial v^2}\right) _T=0\). These two conditions together with the above Van der Waals EOS lead to \(v_c = 3b\), \(p_c = \frac{a}{27b^2}\), \(T_c = \frac{8a}{27Rb}\), and \(Z_c = \frac{3}{8}\). we substitute the definitions of the reduced pressure, reduced molar volume and reduced temperature with the values of critical properties (in terms of a and b) into the above Van der Waals EOS to obtain the following reduced form of the Van der Waals equation

$$\begin{aligned} \left( p_r + \frac{3}{v_r^2}\right) \left( 3v_r-1\right) = 8T_r. \end{aligned}$$
(2)

From the PVT-form of the Van der Waals EOS, we can derive the bulk Helmholtz free energy for the Van der Waals fluid as \(F(T,V,\mathbf{N}) = F_\mathrm{ideal}(T,V,\mathbf{N})\,+\,F_\mathrm{excess}(T,V,\mathbf{N}),\) where the ideal gas contribution (note that \(N=\sum _i N_i\)) is \(F_\mathrm{ideal}(T,V,\mathbf{N}) = R T \sum _{i=1}^M N_i \ln \left( \frac{N_i}{V} \right) + N C_\mathrm{intg}(T),\) with one choice of \(C_\mathrm{intg}\) as \(C_\mathrm{intg}(T) = -RT \left( \ln (\frac{n_Q}{N_A}) + 1 \right) \) and \(n_Q = \left( \frac{2\pi m k_B T}{h^2} \right) ^{3/2}.\) The excess part of Helmholtz free energy is: \(F_\mathrm{excess}(T,V,N) = -N R T \ln \left( 1 - \frac{N b}{V} \right) - \frac{N^2 a}{V}.\)

We define the Helmholtz free energy density f as the Helmholtz free energy per unit volume of fluid. It is clear that

$$\begin{aligned} f(T,\mathbf{n}) = f_\mathrm{ideal}(T,\mathbf{n}) + f_\mathrm{excess}(T,\mathbf{n}), \end{aligned}$$
(3)
$$\begin{aligned} f_\mathrm{ideal}(T,\mathbf{n}) = R T \sum _{i=1}^M n_i \ln \left( n_i \right) + n C_\mathrm{intg}(T), \end{aligned}$$
(4)
$$\begin{aligned} f_\mathrm{excess}(T,n) = -n R T \ln \left( 1 - bn \right) - a n^2. \end{aligned}$$
(5)

With the Peng-Robinson EOS, the bulk Helmholtz free energy density \(f(T, \mathbf {n})\) of a bulk fluid is determined by \(f(T,\mathbf{n}) = f_\mathrm{ideal}(T,\mathbf{n}) + f_\mathrm{excess}(T,\mathbf{n})\) and

$$\begin{aligned} f^{\text {ideal}}(T, \mathbf {n})=RT\sum _{i=1}^{M}n_{i}\left( \ln n_{i}-1\right) , \end{aligned}$$
(6)
$$\begin{aligned} f^{\text {excess}}(T, {n})=-nRT\ln \left( 1-bn\right) +\frac{a(T)n}{2\sqrt{2}b}\ln \left( \frac{1+(1-\sqrt{2})bn}{1+(1+\sqrt{2})bn}\right) . \end{aligned}$$
(7)

The two parameters a and b can be computed as follows. For a mixture, these parameters can be calculated from the ones of the pure fluids by mixing rules:

$$\begin{aligned} a(T)=\sum _{i=1}^{M}\sum _{j=1}^{M}y_{i}y_{j}a_{i}a_{j}^{1/2}(1-k_{ij}), \quad b=\sum _{i=1}^{M}y_{i}b_{i}, \end{aligned}$$
(8)

where \(y_{i}=n_{i}/n\) is the mole fraction of component i, and \(a_{i}\) and \(b_{i}\) are the Peng-Robinson parameters for pure-substance component i. We often use experimental data to fit the binary interaction coefficient \(k_{ij}\) of Peng-Robinson. For convenience, \(k_{ij}\) is usually assumed to be constant for a fixed species pair.

Even though the pure-substance Peng-Robinson parameters \(a_{i}\) and \(b_{i}\) can also be fit by using experimental data, they can also be computed from the critical properties of the species:

$$\begin{aligned} a_{i}= & {} a_{i}\left( T\right) =0.45724\frac{R^{2}T_{c_{i}}^{2}}{P_{c_{i}}}\left( 1+m_{i}\left( 1-\sqrt{\frac{T}{T_{c_{i}}}}\right) \right) ^{2},\\ b_{i}= & {} 0.07780\frac{RT_{c_{i}}}{P_{c_{i}}}. \end{aligned}$$

As intrinsic properties of the species, the critical temperature \(T_{c_{i}}\) and critical pressure \(P_{c_{i}}\) of a pure substance are available for most substances encountered in engineering practice. In the above formula for \(a_{i}\), we need also to specify the parameter \(m_{i}\) for modeling the influence of temperature on \(a_{i}\). It was suggested that one may correlates the parameter \(m_{i}\) experimentally to the accentric parameter \(\omega _{i}\) of the species by the following equations:

$$\begin{aligned} m_{i}= & {} 0.37464+1.54226\omega _{i}-0.26992\omega _{i}^{2}, \quad \omega _{i}\le 0.49,\\ m_{i}= & {} 0.379642+1.485030\omega _{i}-0.164423\omega _{i}^{2}+0.016666\omega _{i}^{3}, \quad \omega _{i}>0.49. \end{aligned}$$

The accentric parameter can be fit by using experimental data, but if we lack data, we can also compute it by using critical temperature \(T_{c_{i}}\), critical pressure \(P_{c_{i}}\) and the normal boiling point \(T_{b_{i}}\):

$$\begin{aligned} \omega = \frac{3}{7}\left( \frac{\log _{10}\left( \frac{P_{c_{i}}}{14.695\text { PSI}}\right) }{\frac{T_{c_{i}}}{T_{b_{i}}}-1}\right) -1 = \frac{3}{7}\left( \frac{\log _{10}\left( \frac{P_{c_{i}}}{1\text { atm}}\right) }{\frac{T_{c_{i}}}{T_{b_{i}}}-1}\right) -1. \end{aligned}$$

Based on the fundamental relation on thermodynamic variables, the pressure of homogeneous fluids p and the Helmholtz free energy \(f(\mathbf {n})\) can be linked in the following way

$$\begin{aligned} p&=p(\mathbf {n},T)=-\left( \frac{\partial F\left( \mathbf {n},T,\varOmega \right) }{\partial V}\right) _{T,\mathbf {N}} =-\left( \frac{\partial \left( f\left( \frac{\mathbf {N}}{V},T\right) V\right) }{\partial V}\right) _{T,\mathbf {N}}\\&=-f-V\sum _{i=1}^{M}\left( \frac{\partial f}{\partial n_{i}}\right) _{T,n_{1},\cdots ,n_{i-1},n_{i+1},\cdots n_{M}}\left( \frac{\partial \frac{N_{i}}{V}}{\partial V}\right) _{N_{i}}\\&=\sum _{i=1}^{M}n_{i}\left( \frac{\partial f}{\partial n_{i}}\right) _{T,n_{1},\cdots ,n_{i-1},n_{i+1},\cdots n_{M}}-f =\sum _{i=1}^{M}n_{i}\mu _{i}-f. \end{aligned}$$

Substitution of the Peng-Robinson expression of f leads to

$$\begin{aligned} p =\frac{nRT}{1-bn}-\frac{n^{2}a(T)}{1+2bn-b^{2}n^{2}} =\frac{RT}{v-b}-\frac{a(T)}{v(v+b)+b(v-b)}. \end{aligned}$$

2.2 Modeling the Two-Phase Systems with Interfaces

The total Helmholtz energy \(F^\mathrm{tot}\) has two contributions, one from the homogeneous bulk fluid, and another one from the interface between the two phases:

$$\begin{aligned} F^\mathrm{tot} \left( \mathbf {n}\right)= & {} F_\mathrm{bulk}\left( \mathbf {n};\, T,\varOmega \right) + F_\mathrm{interface} \left( \mathbf {n};\, T, A_I\right) \\= & {} \int _{\varOmega \setminus A_I}f(\mathbf {n};\, T)d\mathbf {x} + \int _{A_I} \sigma (\mathbf {[n]},\mathbf {\{n\}};\, T)ds, \end{aligned}$$

where \(\sigma \) is the interfacial tension (or the interfacial Helmholtz energy per unit area), which is a function of the jump and average of \(\mathbf {n}\) across the interface. In the paper, for convenience, we assume that \(\sigma \) is a given constant. With this assumption, we have

$$F^\mathrm{tot} = f\left( \mathbf {n}^L\right) V^L + f\left( \mathbf {n}^G\right) V^G + \sigma A_I,$$

where \(\mathbf {n}^L\) and \( V^L\) are the molar density and volume of liquid phase, respectively. Meanings of \(\mathbf {n}^G\) and \( V^G\) are similar to these.

We impose the total volume \(V^\mathrm{tot}\) and the total moles \(\mathbf {N}^\mathrm{tot}\) such as \(V^\mathrm{tot} = V^L + V^G\) and \(N_i^\mathrm{tot} = n_i^L V^L + n_i^G V^G\), \(i=1,2,\cdots ,M.\) The interface \(A_I\) is also a function of \(V^G\); for example, if there is only one single bubble, then \(A_I = 4\pi r^2\) while \(V^G = \frac{4\pi r^3}{3}\). Considering \(\mathbf {n}^G\) and \(V^G\) as the primary variables, we can write the total Helmholtz energy \(F^\mathrm{tot}\) as

$$F^\mathrm{tot}(\mathbf {n}^G, V^G) = f\left( \frac{\mathbf {N}^\mathrm{tot} -\mathbf {n}^G V^G}{V^\mathrm{tot}-V^G} \right) (V^\mathrm{tot}-V^G) + f\left( \mathbf {n}^G\right) V^G + \sigma A_I(V^G).$$

When \(F^\mathrm{tot}\) achieves its minimum and the function \(F^\mathrm{tot}\) is smooth, we have

$$\frac{\partial F^\mathrm{tot}(\mathbf {n}^G, V^G)}{\partial n_i^G} = 0, \quad \mathrm{and} \quad \frac{\partial F^\mathrm{tot}(\mathbf {n}^G, V^G)}{\partial V^G} = 0,$$

which implies respectively \(\mu (\mathbf {n}^G) = \mu (\mathbf {n}^L),\) and \(p^G - p^L = \sigma \frac{dA_I}{dV^G}.\) If there is only one single bubble, then \(A_I = 4\pi r^2\) and \(V^G = \frac{4\pi r^3}{3}\); the above condition simplifies to

$$p^G - p^L = \frac{2\sigma }{r},$$

which is known as the Young-Laplace equation.

3 Numerical Methods

To find the minimum of the \(F^\mathrm{tot}(\mathbf {n}^G, V^G)\), we design the following ordinary differential equation (ODE) system

$$\begin{aligned} \frac{\partial N_i^G}{\partial t}= & {} -\frac{\partial N_i^L}{\partial t} = k_{N_i}N_i^\mathrm{tot}(\mu (\mathbf {n}^L) - \mu (\mathbf {n}^G)),\quad i=1,2,\cdots ,M, \\ \frac{\partial V^G}{\partial t}= & {} -\frac{\partial V^L}{\partial t} = k_{V}V^\mathrm{tot}(p^G(\mathbf {n}^G) - p^L(\mathbf {n}^L) -\sigma \frac{dA_I}{dV^G}). \end{aligned}$$

We solve the above ODE using the following explicit Euler method, which is conditionally stable provided that the time step is less than a certain value.

$$\begin{aligned} \frac{N_i^{G,k+1}-N_i^{G,k}}{t^{k+1}-t^k}= & {} k_{N_i}N_i^\mathrm{tot}(\mu (\mathbf {n}^{L,k}) - \mu (\mathbf {n}^{G,k})),\quad i=1,2,\cdots ,M, \\ \frac{V^{G,k+1}-V^{G,k}}{t^{k+1}-t^k}= & {} k_{V}V^\mathrm{tot}(p^G(\mathbf {n}^{G,k}) - p^L(\mathbf {n}^{L,k}) -\sigma \frac{dA_I}{dV^G}), \end{aligned}$$

where \(\mathbf {n}^{G,k}\) and \(\mathbf {n}^{L,k}\) can be determined from \(V^{G,k}\), \(N_i^{G,k}\), \(i=1,2,\cdots ,M.\)

Existence, uniqueness, and energy-decay property of the solution to the above ODE equation system as well as the existence, uniqueness, and the conditional energy-decay property of the numerical solution defined above can be proved using techniques similar to the ones used in our previous work [10].

4 Numerical Examples

4.1 Effect of Capillary Pressure on the Van der Waals Fluid

We first consider the single-component two-phase fluid system modeled by the Van der Waals EOS. If we write the single-component Van der Waals EOS using reduced temperature, reduced pressure and reduced molar volume, we then may obtain an universal dimensionless Eq. (2) for the EOS. That is, all single-component two-phase fluids behave in the same way after certain linear transformation. Since we will report results in reduced quantities, it does not matter the parameters a and b we choose, but in the implementation, the choice of \(k_{N}\) and \(k_{V}\) might depend on the specific values of a and b. Without loss of generality, we choose \(a=3\) and \(b=\frac{1}{3}\); in this way, \(p_c=1\) and \(v_c=1\) and thus \(p_r = p\) and \(v_r = v.\) The units of these quantities can be any fixed units, as long as consistent units are used, such as the SI unit system.

Fig. 1.
figure 1

Phase diagrams of the single-component Van der Waals fluid: (a) reduced molar volumes of gas and liquid phases at equilibrium as a function of reduced temperature; (b) reduced boiling temperature as a function of reduced pressure

Fig. 2.
figure 2

The effect of capillary pressure on the saturation pressures of the liquid phase and the vapor (gas) bubble modeled by the single-component Van der Waals EOS: (a) \(T_r = 0.7\) (b) \(T_r = 0.8\) (c) \(T_r = 0.9\) (d) \(T_r = 0.95\)

Fig. 3.
figure 3

The effect of capillary pressure on the molar densities of the liquid phase and the vapor (gas) bubble modeled by the single-component Van der Waals EOS: (a) \(T_r = 0.7\) (b) \(T_r = 0.8\) (c) \(T_r = 0.9\) (d) \(T_r = 0.95\)

Fig. 4.
figure 4

The effect of capillary pressure on the saturation pressures of the liquid droplet and the vapor (gas) phase modeled by the single-component Van der Waals EOS: (a) \(T_r = 0.7\) (b) \(T_r = 0.8\) (c) \(T_r = 0.9\) (d) \(T_r = 0.95\)

Fig. 5.
figure 5

The effect of capillary pressure on the molar densities of the liquid droplet and the vapor (gas) phase modeled by the single-component Van der Waals EOS: (a) \(T_r = 0.7\) (b) \(T_r = 0.8\) (c) \(T_r = 0.9\) (d) \(T_r = 0.95\)

Fig. 6.
figure 6

The effect of capillary pressure on phase behaviors of a vapor (gas) bubble immersed in its saturated liquid phase modeled by the single-component Peng-Robinson EOS: (a) liquid and vapor (gas) pressures (b) molar densities of two phases

Fig. 7.
figure 7

The effect of capillary pressure on phase behaviors of a liquid droplet immersed in its vapor (gas) phase modeled by the single-component Peng-Robinson EOS: (a) liquid and vapor (gas) pressures (b) molar densities of two phases

For comparison and verification, we first let \(\sigma =0\) in our model, and we calculate the liquid-vapor phase behavior of the Van der Waals fluid without capillarity. We generate phase diagrams numerically and plot them in Figs. 1(a) and (b). Figure 1(a) is the volume-temperature phase envelope while Fig. 1(b) displays reduced boiling temperature as a function of reduced pressure.

We then consider two cases with capillary pressure. The first case is a single bubble of radius r immersed in the liquid. In this case, \(V^G = \frac{4\pi r^3}{3}\) and \(dA_I/dV^G = 2/r.\) For this single-component two-phase fluid system, our modeling ODE reduces to

$$\begin{aligned} \frac{\partial N^G}{\partial t}= & {} k_{N}N^\mathrm{tot}(\mu ({n}^L) - \mu ({n}^G)), \\ \frac{\partial r}{\partial t}= & {} \frac{k_{V}V}{4\pi r^2}^\mathrm{tot}(p^G({n}^G) - p^L({n}^L) -\frac{2\sigma }{r} ). \end{aligned}$$

The values of \(k_N\) and \(k_V\) are manually tuned for one typical simulation, and are then fixed for all other runs. In all numerical examples in this subsection, we choose \(k_N = 0.05\) and \(k_V = 0.02\). We use the unit time step \(\varDelta t = 1\) for all numerical runs in this paper.

The effect of capillary pressure on the saturation pressure of the liquid phase and the vapor (gas) phase is provided in Figs. 2(a)–(d) under various condition of reduced temperatures and capillary pressures. The horizontal axis of these plots are \(\frac{\sigma }{r p_c}\). The ratio of interfacial tension to the radius of the gas bubble \(\frac{\sigma }{r}\) is the influencing factor, we divide this ratio by the critical pressure of the fluid \(p_c\) so that we get a dimensionless quantity.

We also compare the our numerical prediction of vapor pressure with the one calculated by the Kelvin equation. The Kelvin equation can be derived by using the two equilibrium conditions \(dp^G - dp^L = d(2\sigma /r)\) and \(d\mu ^G=d\mu ^L\), and the two fundamental relations \(d\mu ^G=v^Gdp^G\) and \(d\mu ^L=v^Ldp^L\), which yield \(dp^G = \frac{d(\frac{2\sigma }{r})}{(1\,-\, \frac{v^G}{v^L})}.\) By assuming the ideal gas law for the \(v^G\) and assuming a constant molar density for \(v^L\), we can integrate the above differential equation to obtain the following Kelvin equation

$$p^G = p^{\infty } \exp (-\frac{2\sigma v^L}{rRT}),$$

which quickly reveals that the vapor pressure decreases with increasing interface curvature and increasing interfacial tension.

Figures 2(a)–(d) indicates that the vapor pressure predicted by the Kelvin equation agree well with our numerical simulation when the temperature is not too low. When the temperature is low, the vapor pressure predicted by the Kelvin equation has pronounced derivation, because the ideal gas law is no loner a good approximation at low temperature. From Figs. 2(a)–(d), we also see that the saturation liquid pressure departs from its zero-capillary value much more significantly than the saturation vapor pressure.

We observe the number of time steps required for convergence varies with the reduced temperature. When \(T_r = 0.7\), the number of time steps required for convergence is about 10 to 30 depending on the value of \(\frac{\sigma }{r p_c}\), while when \(T_r = 0.95\), the number of time steps required for convergence is can be more than 2000.

In Figs. 3(a)–(d), we show the variation of molar densities of the saturated liquid phase and saturated vapor phase as influenced by capillarity. It is evident that the molar densities of both phases decrease with increasing interface curvature and increasing interfacial tension. This is likely because the decrease of saturation pressure causes the decrease of molar densities.

The second case is a single liquid droplet of radius r immersed in the vapor phase. In this case, \(V^G = V^\mathrm{tot} - \frac{4\pi r^3}{3}\) and \(dA_I/dV^G = -2/r.\) For this single-component two-phase fluid system, our modeling ODE reduces to

$$\begin{aligned} \frac{\partial N^G}{\partial t}= & {} k_{N}N^\mathrm{tot}(\mu ({n}^L) - \mu ({n}^G)), \\ \frac{\partial r}{\partial t}= & {} \frac{k_{V}V}{4\pi r^2}^\mathrm{tot}( p^L({n}^L) - p^G({n}^G) -\frac{2\sigma }{r} ). \end{aligned}$$

Figures 4(a)–(d) provide the trend of the saturation liquid pressure of a single liquid droplet and the surrounding saturation vapor pressure as influenced by capillarity. Unlike the vapor bubble case, both the liquid pressure and the saturation vapor pressure increase here with increasing interface curvature and increasing interfacial tension. The trend of the saturation vapor pressure is also compared with the ones predicted by the Kelvin equation, which does a good job only when the temperature is not too low. In Figs. 5(a)–(d), we show the variation of molar densities of the saturated liquid phase and saturated vapor phase as influenced by capillarity. Unlike the vapor bubble case, the molar densities of both phases increase with increasing interface curvature and increasing interfacial tension.

Table 1. Relevant data of isobutane (excerpted from Table 3.1 of [4], Page 141)

4.2 Effect of Capillary Pressure on the Peng-Robinson Fluid

The Peng-Robinson EOS, unlike the single-component Van der Waals EOS, cannot transformed into a same reduced from by using reduced temperature, reduced pressure and reduced molar volume, because a third parameter, the accentric parameter appears in the Peng-Robinson EOS. We thus consider a specific example below with physical units. In the example below, we consider the species of isobutane (\(n\text {C}_{4}\)) at the temperature of 350 K. For convenience of readers, we provide in Table 1 its critical properties and its normal boiling point. The values of \(k_N\) and \(k_V\) are manually tuned for one typical simulation, and are then fixed for all other runs. In all runs in this subsection, we choose \(k_N = 1\times 10^{-3}\, \text {mol}\,\text {J}^{-1}\,\text {s}^{-1}\) and \(k_V = 1\times 10^{-9}\,\text {Pa}^{-1}\,\text {s}^{-1}\). We again use the unit time step \(\varDelta t = 1\,\text {s}.\)

Figure 6(a) displays the effect of capillary pressure on the saturation pressures of both phases when a single gas bubble is in equilibrium on its liquid phase. The corresponding variation of molar densities of both phases are plotted in Fig. 6(b). The trend of the saturation vapor pressure is also compared with the ones predicted by the Kelvin equation, which does a reasonable job in the condition simulated in this example. Clearly, like the case modeled by the Van der Waals EOS, both the liquid pressure and the saturation vapor pressure decrease here with increasing interface curvature and increasing interfacial tension, which also leads a decrease of molar density of both phases.

Figure 7(a) and (b) show the variation of the saturation pressures and molar densities of both phases when a single liquid droplet is in equilibrium on its vapor phase under various conditions of interface curvature and increasing interfacial tension. The trend seems similar to the case predicted by the single-component Van der Waals EOS qualitatively.

5 Conclusion

In this paper, we formulate a framework to model a two-phase fluid with a (sharp) interface using Van der Waals EOS and Peng-Robinson EOS. Our model is able to predict the affect of capillary pressure on phase behaviors. We propose simple but conditionally energy stable numerical schemes to solve the minimization problem of total Helmholtz energy in given total volume, temperature, and total moles. We consider bubbles in liquid as well as droplets in vapor in our numerical examples. Our numerical agree well with the prediction of Kelvin’s equation when the temperature is not too low. Due to the limitation of the space, we provide only single-component examples. Our model and schemes, however, are expected to work well with multi-component systems, which will be investigated and reported in a separate paper elsewhere.