The increasing demand for oil products from limited oil reserves has strongly motivated the research in the area of reservoir simulation to predict and optimize the production facilities. Out of the variety of contemporary problems related to petroleum engineering, secondary and tertiary (or enhanced) oil recovery involving injection of fluids and chemical reactions have become important production processes today. Simulations of oil recovery consist in solving numerically the differential equations that govern the multiphase flow in a porous medium, which have no analytical solution, taking into account all geological, physical, chemical, and mechanical phenomena occurring during operation. The analysis of all these phenomena can be performed, up to some extent, by means of laboratory experiments and/or field tests at a small scale. While such experiments and tests are in general expensive and difficult to conduct, the use of mathematical models aimed at reproducing such experiments on rock core samples are becoming increasingly more important. In addition to more traditional numerical techniques, such as finite-difference, finite-element, and discrete-element methods, Smoothed Particle Hydrodynamics (SPH) is becoming a potentially promising and powerful alternative technique. Although many tasks and challenges still remain, the great success of SPH developments and applications in engineering and science has recently motivated its application to the multiphase flow of immiscible and miscible fluids in heterogeneous porous media. In this paper, we briefly review the fundamentals of the SPH method for multiphase flows along with its most recent developments and applications to flow in oil reservoir rocks.

Keywords: particle methods, numerical modeling, multiphase flow, oil reservoir rock

SPH, smoothed particle hydrodynamics; SOR, secondary oil recovery; EOR, enhanced oil recovery; CBF, continuum boundary force; PF, pairwise force

In the petroleum industry the implementation of new technologies for optimizing the efficiency of oil extraction is of primary importance. In general, oil production is divided into three main stages: primary, secondary (SOR), and tertiary or enhanced oil recovery (EOR). In particular, SOR involves the injection of water and gas to push the oil through the well and drive it to the surface, while EOR entails changing the properties of hydrocarbons through the injection of steam, chemicals (i.e., polymers), and gas (like, for example, natural gas, nitrogen, or carbon dioxide) in order to enhance the ability of trapped oil to flow through the reservoir. After primary and secondary recovery more than 70% of the oil remains trapped in the porous rock, while using EOR the efficiency of extraction can increase from 30 to 75%.¹ An accurate determination of the physical characteristics of oil reservoirs is still beyond the reach at the moment, and so statistical tools together with complex laboratory experiments and field tests are used to extrapolate these properties. In order to reduce the model uncertainties, oil companies are addressing much effort to the research and development of new technologies. Among these, the use of numerical simulations has become a useful tool to complement traditional type experiments for a large number of phenomena occurring during oil extraction operations.^2–6

Numerical models aimed at solving oil extraction problems involve mainly the solution of the nonlinear fluid equations for multiphase and multicomponent flows through a heterogeneous medium (porous rock).^3,3–9 When describing the flow at a microscopic level (or pore scale), the Navier-Stokes equations coupled to the mass conservation law for multiphase flow must be solved, while at a macroscopic level Darcy´s law is employed in place of the Navier-Stokes equations to describe the flow behavior. At the microscopic scale, the fluid motion is mainly controlled by the effects of interfacial tension, viscosity, and capillary pressures. In compositional flows, chemical components between coexisting phases may react, modifying the fluid phase properties. In addition, the interactions between the multicomponent fluid phases and the porous rock may also alter the wettability. Other complexities come from the heterogeneity and anisotropy of the grain structure especially when flow simulations using digital images of micro-tomographic data of real rock core samples are required to numerically reproduce the experimental set-up and results from laboratory displacement tests.¹⁰

In general, most mathematical models relying on a continuum approach involve developing conceptual models.^11–16 They incorporate the geometrical information of the fracture-matrix system, define the conservation laws for the fracture-matrix domain, and solve a number of discrete nonlinear algebraic and constitutive equations, which express relations and constraints of physical processes, variables, and parameters as functions of primary unknowns.

In this paper we present a short overview of recent applications and developments of SPH to simulations of oil reservoir rocks addressed principally to reproduce laboratory displacement tests on rock core samples.

The governing equations for multiphase flow can be written in terms of the continuity and Navier-Stokes equations for each phase α as

$\frac{\text{d}{\rho }_{\alpha }}{\text{d}t}= -{\rho }_{\alpha }\left(\nabla .v_{\alpha }\right)$ ………………… (1)

${\rho }_{\alpha }\frac{\text{d}{v}_{\alpha }}{\text{d}t}= -\nabla p_{\alpha }+\nabla .T_{\alpha }+ {\rho }_{\alpha } g$ ………………………. (2)

where ${\rho }_{\alpha }$ is the fluid density, $v_{\alpha }$ is the fluid velocity vector, $p_{\alpha }$ is the pressure, $T_{\alpha }=\left[{\mu }_{\alpha }\left(\nabla v_{\alpha }+ \nabla v_{\alpha }^t\right)\right]$ is the viscous stress tensor, ${\mu }_{\alpha }$ is the dynamic viscosity, g is the gravitational acceleration, and d⁄dt  is the total or material derivative. A constitutive equation of state (EOS) for each phase, $p_{\alpha }=f\left({\rho }_{\alpha }\right)$ is commonly employed to close the system of equations (1) and (2). The solution of these equations must be subjected to no-slip conditions at the fluid-solid boundary, i.e., the normal and tangential fluid velocity for each fluid phase must vanish at the solid surface.^10,17

In SPH applications of multiphase flow in a porous medium, the nonlinearities of these equations at the fluid-fluid interfaces are modeled using the Young-Laplace boundary condition for pressure and velocity, while at the fluid-fluid-solid interfaces the Young formula is used to prescribe the contact angle at the contact line.^18-20

For chemically reacting fluid phases, terms accounting for the interfacial transfer rate of mass $I_{\alpha }$ and momentum $M_{\alpha }$ must be added as sources on the right-side of Equations (1) and (2), respectively [7]. On the other hand, if heat transfer effects are important, Equations (1) and (2) must be complemented with a transport equation for the internal energy or enthalpy. However, a simplification of the above model can be obtained in the framework of continuum mixture theory.²¹ While the transport phenomena are mathematically described by the basic principles of conservation for each phase separately and by appropriate kinematic and dynamic conditions at the interfaces, this difficulty can be more easily covered by mixture theory by making use of the volume fraction ${\phi }_{\alpha }$ occupied by phase α, which is defined as a scalar function of position and time such that $0 \le {\phi }_{\alpha }\le 1$ . Expressing the continuity equation in terms of the volume fraction (or phase saturation) and solving Equation (2) for the mixture velocity reduces the number of variables necessary to provide a description of the multiphase flow compared to full multiphase models. In addition, there is no need to define and track the interfaces between coexisting fluid phases. In spite of these simplifications, there is at present a lack of SPH implementations relying on the mixture model approach.

The method of SPH was originally developed for applications to astrophysical flows^23,24 and since then it has been applied to a broad range of areas related to fluid and solid mechanics.²⁴ In SPH, a fluid is represented by a set of interpolation points (or particles), which are then advected with the fluid velocity in a Lagrangian sense.

A continuous field function $f\left(x\right)$ is approximated by a two-step procedure: (1) a kernel or weighted interpolation estimate of f(x), say $\langle f\left(x\right)\rangle$ defined by the integral

$\langle f\left(x\right)\rangle ={{\displaystyle \int }}^{\text{}}f\left(x\text{'}\right)W\left(x-x^{\prime },h\right)\text{d}{x}^{\prime }$ …………………… (3)

where $W\left(x-x\text{'},h\right)$ is referred to as the smoothing (or kernel) function and h is the so-called smoothing length specifying the width of the kernel; and (2) a particle approximation where the kernel estimate of the function at the position of particle i, say $f_i=\langle f\left(x_i\right)\rangle ,$ is evaluated by replacing the integral in Equation (3) by a summation over all neighbors of particle i lying within the kernel support of radius $\le kh$ (where k is some real number greater or equal to one), that is

$f_i={\displaystyle \sum }_{j=1}^n\frac{f_j}{n_j}W\left(x_i-x_j, h\right)$ ………………………. (4)

Where $n_j={\rho }_j/m_j$ is the particle number density, m_j is the mass, and ${\rho }_j$ is the density of particle $j$ at position $x_j$ , while the ratio $1/n_j=m_j/{\rho }_j$ is used to represent the volume occupied by particle j in order to ensure consistency with the kernel approximation (3). In passing, we note that consistency of the kernel approximation demands that the volume integral of W must be equal to unity and that the limit of W when $h\to 0$ must tend to the Dirac-δ distribution. Similar expressions to (3) and (4) follow for the gradient of f with< $f\to \nabla f$ . One important realization when multiple fluid phases with significantly different densities are present is to use n_j rather than the ratio $m_j/{\rho }_j$ in Equation (4).^19,25 This will effectively remove density discrepancies and artificial surface tension effects that would otherwise appear at phase interfaces due to the jump in density between coexisting phases.

The SPH representation of the particle number density of particle i within phase a is readily obtained from Equation (4) as

$n_{\alpha i}={\displaystyle \sum }_{j=1}^n{W}_{ij}$ ………….. (5)

where $W_{ij}=W\left(x_i-x_j, h\right)$ and the sum is over all n neighbors within the compact support of the kernel.

This form guarantees exact mass conservation. The most commonly employed SPH representation of the Navier-Stokes equations (2) for multiphase flow is obtained using the pressure gradient and viscous dissipation particle approximations derived by Morris et al.²⁶ which in terms of the particle number density is

$m_{\alpha i}\frac{\text{d}{v}_{\alpha i}}{\text{d}t}=-{\displaystyle \sum }_{j=1}^n\left(\frac{p_{\alpha i}+p_{\alpha j}}{n_{\alpha i}{n}_{\alpha j}}\right){\nabla }_i{W}_{ij}+{\displaystyle \sum }_{j=1}^n\frac{4{\mu }_{\alpha i}{\mu }_{\alpha j}}{{\mu }_{\alpha i}+{\mu }_{\alpha j}}\frac{\left(v_{\alpha i}-v_{\alpha j}\right)}{n_{\alpha i}{n}_{\alpha j}}\frac{x_{ij}.{\nabla }_i{W}_{ij}}{x_{ij}\text{'}{x}_{ij}}+m_{\alpha i}g+F_{\alpha i}$ ………………… (6)

where $x_{ij}=x_i-x_j$ and the added force $F_{\alpha i}$ accounts for the interfacial tension on the fluid-fluid interface and/or the wetting and non-wetting behaviors at the fluid-fluid-solid interface. In general, for immiscible phases these forces are modeled in SPH using two different approaches to impose the Young-Laplace boundary condition at the fluid-fluid interface and the Young condition (contact angle) at the fluid-fluid-solid interface. One approach is based on the so-called Continuum Boundary Force (CBF) scheme ,²⁷ where these forces can be calculated directly from the boundary conditions. The other approach is based on the assumption that these forces are pair wise molecular-like forces.^18-20 In this approach the force $F_{\alpha i}$ >when the field particle i belongs to phase α is defined as

$F_{\alpha i}$ with $F_{ij}=-F_{\alpha \beta }\left(x_{ij}\right)\frac{x_{ij}}{x_{ij}}$ ………………… (7)

where $F_{ij}$ is a pair wise interaction force acting between particle i of phase α and particle j of phase β, where β can be a fluid or a solid phase.²⁰ To mimic the interaction forces at a molecular level, $F_{\alpha \beta }$ is defined to be short-range repulsive for $x_{ij}\le x<kh$ and long-range attractive for $x<x_{ij}\le kh$ where $x_{ij}^2=x_{ij}.x_{ij}$ and x is some distance from the position of particle i within the kernel support. If the wall is wetted by fluid α early SPH simulations using this approach employed a Lennard-Jones force to define $F_{\alpha \beta }$ in Equation (7), while the non-wetting behavior was simulated defining $F_{\alpha \beta }=-C/x_{ij}$ where C is a coefficient characterizing the non-wetting intensity of the fluid.¹⁸ However, a drawback of this scheme is that for both types of forces the wetting and non-wetting intensities are specified by coefficients that must be defined to produce the desired behaviors. In other words, they act as free parameters that must be set arbitrarily. A more sophisticated model relying on a closed-form relation between the force $F_{\alpha \beta }$ in Equation (7) and the macroscopic properties of the multiphase system was recently derived by Tartakovsky and Panchenko,²⁰ the so-called Pairwise Force SPH (PF-SPH), which is based on Gibbs thermodynamics theory.

As a final remark, the treatment of no-slip boundary conditions has evolved to the point where it is now possible to simulate experimental tests on real-world rock core samples. Such boundary conditions use a state specific form of the particle number density to approximate a fluid particle’s proximity to a wall boundary.¹⁰ This method produces accurate no-slip conditions for highly irregular wall shapes, thereby allowing flow simulation and characterization using digital images of X-ray micro-tomographic data of rock samples.

Two-dimensional pore-scale PF-SPH models of immiscible flow in heterogeneous porous media were first generated by using randomly located, non-intersecting circular grains of different sizes.¹⁹ In these models, pore-scale anisotropy was also investigated by creating a small fracture in the middle of the porous medium. The results for these rather idealized rock geometries were found to be in good agreement with analytical solutions and experimental observations. In particular, the amount of non-wetting fluid that is not displaced by injection of the wetting fluid into the porous medium was found to depend on the flow direction relative to that of the micro-fracture. Such simulations have demonstrated the ability of SPH to account for the effects of pore-scale heterogeneity on concentration evolution and estimate transport properties of the porous medium such as mass transfer and dispersion coefficients.

SPH models of oil displacement in cavity-fracture structures were also performed using Lennard-Jones repulsive forces for the wetting phase and simple repulsions of the form $F_{\alpha \beta }=-C/x_{ij}$ for the non-wetting phase.¹⁸ For vertically oriented fractures connecting with a circular cavity, water was seen to flow in the form of plug flow displacing the oil with an efficiency of oil recovery close to 100% when the wall is water-wet. For horizontal fractures, water was seen to rapidly sink down the cavity with the water-oil interface gradually rising until reaching the height of the connecting fractures. In both cases, the numerical simulations were found to compare well with corresponding experiments. For models in which the cavity connects to a top horizontal fracture, these simulations have shown that the water flows in a continuous way for water-wet wall and that there exists a critical width of the vertical fracture beyond which the oil in the cavity can be driven out. In this case, as the oil density and the inflow velocity is lowered the critical width becomes correspondingly smaller for water-wet wall.

Simulations of three-phase flow in a porous medium were recently performed using a novel formulation of the PF method.²⁰ Flow simulations were carried out for two liquids and a gas in a two-dimensional heterogeneous porous medium. As is typical for water/oil/gas systems, the liquids were assumed to be 1000 times denser and 100 times more viscous than the gas phase. The effects of wettability and surface tensions between the phases were studied for different scenarios showing acceptably good agreement with experimental tests. While most existing grid-based multiphase models are limited to no more than two phases, these PF-SPH calculations possibly represent one of the very few cases that can readily deal with three-phase flows in a porous medium.

As a next step, future SPH simulations at the pore scale should point toward three-dimensional models and the analysis of flow through real model geometries derived from X-ray micro-tomographic data in order to reproduce true experimental tests on rock core samples. The development of new no-slip boundary conditions for low Reynolds number flows that can now handle the degree of complexity of rock samples is bringing SPH to become relatively close to these goals. An important aspect of this type of simulations is that they will allow direct prediction of the rock properties such as the absolute and relative permeability, which are critical to the characterization of oil and gas reservoirs. Numerical permeability measurements are also key factors for managing and enhancing resource recovery.

In this paper, we have presented a brief overview of the most recent developments and applications of the SPH method to multiphase flow in oil reservoir rocks. Many tasks and challenges still remain that will be key to acceptance of SPH by the industry as a reliable alternative to more traditional reservoir modeling tools. While SPH presents true advantages over other techniques when adding the complexity of multiple fluid phases and rock heterogeneity and anisotropy, future developments and applications must include simulations of chemically reacting compositional (multiphase and multicomponent) flows as well as a unified Lagrangian formulation for both elastic wall and immiscible fluids, thus allowing deformation of the solid matrix due to changes of the inner pressure. As a further prospect, flow simulations from digital images of real rocks would demand a multi-scale approach to solve for small pores and fractures with enough accuracy and efficiency. A promising way is to use adaptive kernels based on density estimations as the ones proposed by Silverman,²⁸ which have been successfully employed to resolve free surfaces, fluid-fluid interfaces, and sharp discontinuities accurately.^29,30

None.

The author declares no conflict of interest.

1. A Firoozabadi. Recovery mechanisms in fractured reservoirs and field performance. Journal of Canadian Petroleum Technology. 2000;39(11):13–17.
2. JR Franchi. Principles of Applied Reservoir Simulation. UK: Elsevier; 2006.
3. Z Chen, G Huan, Y Ma (2006) Computational Methods for Multiphase Flows in Porous Media. Society for Industrial and Applied Mathematics. Philadelphia, USA, 2006. 531 p.
4. M Cancelliere, D Viberti, F Verga. A step forward to closing the loop between static and dynamic reservoir modeling. Oil & Gas Science and Technology–Revue d’IFP Energies Nouvelles. 2014;69(7):1201–1225.
5. KA Lie, BT Mallison. Mathematical models for oil reservoir simulation. In: B Engquist, editor. Encyclopedia of Applied and Computational Mathematics. Springer–Verlag. Germany, 2015. p. 850–856.
6. P Druetta, P Tesi, C De Persis, et al. Methods in oil recovery processes and reservoir simulation. Advances in Chemical Engineering and Science. 2016;6:399–435.
7. MB Allen. Numerical modelling of multiphase flow in porous media. Advances in Water Research. 1985;8:162–187.
8. PM Adler, H Brenner. Multiphase flow in porous media. Annual Review of Fluid Mechanics. 1988;20:35–59.
9. J Bear. Dynamics of Fluids in Porous Media. Dover, New York, USA, 1988. 764 p.
10. DW Holmes, JR Williams, P Tilke, et al. Characterizing flow in oil reservoir rock using SPH:Absolute permeability. Computational Particle Mechanics. 2016;3(2):141–154.
11. M Bai, D Elsworth, JC Roegiers. Multiporosity/multipermeability approach to the simulation of naturally fractured reservoirs. Water Resources Research. 1993;29(6):1621–1633.
12. S Stothoff, D Or. A discrete–fracture boundary integral approach to simulating coupled energy and moisture transport in a fractured porous medium. In: B Faybishenko, PA Witherspoon, SM Benson. editor. Dynamics of Fluids in Fractured Rocks, AGU Geophysical Monograph, American Geophysical Union. Washington, USA, 2000;122:269–279.
13. Z Chen. Homogenization and simulation for compositional flow in naturally fractured reservoirs. Journal of Mathematical Analysis and Applications. 2007;326(1):12–32.
14. YS Wu, G Qin. A generalized numerical approach for modeling multiphase flow and transport in fractured porous media. Communications in Computational Physics. 2009;6(1):85–108.
15. P Lemonnier, B Bourbiaux. Simulation of naturally fractured reservoirs. State of the art. Part 1. Physical mechanisms and simulator formulation. Oil & Gas Science and Technology –Revue d IFP Energies Nouvelles. 2010;65(2):239–262.
16. P Lemonnier, B Bourbiaux. Simulation of naturally fractured reservoirs. State of the art. Part 2. Matrix–fracture transfers and typical features of numerical studies. Oil & Gas Science and Technology–Revue d’IFP Energies Nouvelles. 2010;65(2):263–286.
17. DW Holmes, JR Williams, P Tilke. Smooth particle hydrodynamics simulations of low Reynolds number flows through porous media. International Journal for Numerical and Analytical Methods in Geomechanics. 2011;35(4):419–437.
18. G Zhou, Z Chen, W Ge, et al. SPH simulation of oil displacement in cavity–fracture structures. Chemical Engineering Science. 2010;65:3363–3371.
19. AM Tartakovsky, P Meakin.
20. Pore scale modeling of immiscible and miscible fluid flows using smoothed particle hydrodynamics. Advances in Water Resources. 2006;29:1464–1478.
21. AM Tartakovsky, A Panchenko. Pair wise force smoothed particle hydrodynamics model for multiphase flow: Surface tension and contact line dynamics. Journal of Computational Physics. 2016;305:1119–1146.
22. M Manninen, V Taivassalo, S Kallio. On the Mixture Model for Multiphase Flow. VTT Publications 288. Finland, 1996. 67 p.
23. RA Gingold, JJ Monaghan. Smoothed particle hydrodynamics:Theory and application to non–spherical stars. Monthly Notices of theRoyal Astronomical Society. 1977;181(3):375–389.
24.  LB Lucy. A numerical approach to the testing of the fission hypothesis. Astronomical Journal. 1977;82(12):1013–1024.
25. GR Liu, MB Liu. Smoothed Particle Hydrodynamics: A Mesh free Particle Method. Singapore: World Scientific Publishing; 2003. 449 p.
26. XY Hu, NA Adams. A multi–phase SPH method for macroscopic and mesoscopic flows. Journal of Computational Physics. 2006;213:844–861.
27. JP Morris, PJ Fox, Y Zhu. Modeling low Reynolds number incompressible flows using SPH. Journal of Computational Physics. 1997;136(1):214–226.
28. S Adami, X Hu, N Adams. A new surface–tension formulation for multi–phase SPH using a reproducing divergence approximation. Journal of Computational Physics. 2010;229(13):5011–5021.
29. BW Silverman. Density Estimation for Statistics and Data Analysis. Chapman & Hall/CRC. USA, 1986. 176 p.
30. L Di G Sigalotti, H López. Adaptive kernel estimation and SPH tensile instability. Computers & Mathematics with Applications. 2008;55:23–50.
31. L Di G Sigalotti, H López, A Donoso, et al. A shock–capturing SPH scheme based on adaptive kernel estimation. Journal of Computational Physics. 2006;212:124–149.