Abstract
We consider the affect of capillary pressure on the Van der Waals fuid and on the Peng-Robinson fluid by minimizing total Helmholtz energy in given total volume, temperature, and total moles. We propose simple but conditionally energy stable numerical schemes, and we provide interesting numerical examples. We compare our numerical results with the prediction of Kelvin’s equation, indicating that Kelvin’s equation works well only when the temperature is not too low.
This work was supported by funding from King Abdullah University of Science and Technology (KAUST) through grants BAS/1/1351-01.
You have full access to this open access chapter, Download conference paper PDF
Similar content being viewed by others
Keywords
- Interfacial tension
- The Kelvin equation
- Phase diagrams
- Peng-Robinson equation of state
- Van der Waals equation of state
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
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
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
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
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:
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:
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:
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}}\):
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
Substitution of the Peng-Robinson expression of f leads to
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:
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
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
When \(F^\mathrm{tot}\) achieves its minimum and the function \(F^\mathrm{tot}\) is smooth, we have
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
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
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.
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.
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
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
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
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.
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.
References
Breure, B., Peters, C.: Modeling of the surface tension of pure components and mixtures using the density gradient theory combined with a theoretically derived influence parameter correlation. Fluid Phase Equilib. 334, 189–196 (2012)
Dawson, C., Sun, S., Wheeler, M.F.: Compatible algorithms for coupled flow and transport. Comput. Methods Appl. Mech. Eng. 193, 2565–2580 (2004)
Fan, X., Kou, J., Qiao, Z., Sun, S.: A component-wise convex splitting scheme for diffuse interface models with Van der Waals and Peng-Robinson equations of state. SIAM J. Sci. Comput. 39(1), B1–B28 (2017)
Firoozabadi, A.: Thermodynamics of Hydrocarbon Reservoirs. McGraw-Hill, New York (1999)
Kou, J., Sun, S.: Multi-scale diffuse interface modeling of multi-component two-phase flow with partial miscibility. J. Comput. Phys. 318, 349–372 (2016)
Kou, J., Sun, S.: Entropy stable modeling of non-isothermal multi-component diffuse-interface two-phase flows with realistic equations of state. Comput. Methods Appl. Mech. Eng. 341, 221–248 (2018)
Kou, J., Sun, S.: A stable algorithm for calculating phase equilibria with capillarity at specified moles, volume and temperature using a dynamic model. Fluid Phase Equilib. 456, 7–24 (2018)
Kou, J., Sun, S.: Thermodynamically consistent modeling and simulation of multi-component two-phase flow model with partial miscibility. Comput. Methods Appl. Mech. Eng. 331, 623–649 (2018)
Kou, J., Sun, S.: Thermodynamically consistent simulation of nonisothermal diffuse-interface two-phase flow with Peng-Robinson equation of state. J. Comput. Phys. 371, 581–605 (2018)
Kou, J., Sun, S., Wang, X.: An energy stable evolutional method for simulating two-phase equilibria of multi-component fluids at constant moles, volume and temperature. Comput. Geosci. 20(1), 283–295 (2016)
Kou, J., Sun, S., Wang, X.: Linearly decoupled energy-stable numerical methods for multi-component two-phase compressible flow. SIAM J. Numer. Anal. 56(6), 3219–3248 (2018)
Lake, L.W.: Enhanced Oil Recovery. Prentice Hall, Englewood Cliffs, New Jersey (1989)
Li, Y., Kou, J., Sun, S.: Numerical modeling of isothermal compositional grading by convex splitting methods. J. Nat. Gas Sci. Eng. 43, 207–221 (2017)
Moortgat, J., Sun, S., Firoozabadi, A.: Compositional modeling of three-phase flow with gravity using higher-order finite element methods. Water Resour. Res. 47, W05511 (2011)
Peng, D.Y., Robinson, D.: A new two-constant equation of state. Ind. Eng. Chem. Fundam. 15(1), 59–64 (1976)
Sun, S., Wheeler, M.F.: Symmetric and nonsymmetric discontinuous Galerkin methods for reactive transport in porous media. SIAM J. Numer. Anal. 43(1), 195–219 (2005)
Sun, S., Wheeler, M.F.: Local problem-based a \(posteriori\) error estimators for discontinuous Galerkin approximations of reactive transport. Comput. Geosci. 11(2), 87–101 (2007)
Author information
Authors and Affiliations
Corresponding author
Editor information
Editors and Affiliations
Rights and permissions
Copyright information
© 2019 Springer Nature Switzerland AG
About this paper
Cite this paper
Sun, S. (2019). Energy Stable Simulation of Two-Phase Equilibria with Capillarity. In: Rodrigues, J., et al. Computational Science – ICCS 2019. ICCS 2019. Lecture Notes in Computer Science(), vol 11539. Springer, Cham. https://doi.org/10.1007/978-3-030-22747-0_40
Download citation
DOI: https://doi.org/10.1007/978-3-030-22747-0_40
Published:
Publisher Name: Springer, Cham
Print ISBN: 978-3-030-22746-3
Online ISBN: 978-3-030-22747-0
eBook Packages: Computer ScienceComputer Science (R0)