Smoothed particle hydrodynamics method for numerical solution of multidimensional filtering problems

The paper is devoted to solving the filtering problem for a mixture of water, gas, and oil in a homogeneous porous medium. The basic equations of filtration theory are converted into a special form for the numerical approximation with the smoothed particle hydrodynamics (SPH) method. A numerical difference scheme is constructed on the basis of SPH method. An algorithm for setting the boundary conditions is proposed and a number of isothermal 1D and 2D tests on the filtering process simulation for a mixture of water, oil, and gas.


1.Basic equations of filtration theory for mixtures in a porous medium
The problem of numerically simulating the dynamics of a multicomponent medium based on paper [1] is considered. It is assumed that the continuous medium is a system of pores, fractures, and faults, which typical dimensions are small in comparison with the physical size of the problem. The ratio of the pore volume to the total volume is the medium porosity metric:

= /
where is the porosity coefficient, is the pore volume, is the total volume of a given element of the medium. The multicomponent mixture flow in porous media is called "filtration". Here, a "component" is a part of the mixture with its own physical and chemical properties and with no changes in its phase state during the simulation. The main parameters of a moving system are the pressure, temperature, and saturation of each component. Saturation of the k-th component void space is the pore volume fraction occupied by this component in a volume element: where n is the number of components in mixture. The summarised saturation of all components in a volume element equals 1: There is a saturation minimum for each component. These values are called "residual saturations" (subscript ). The simultaneous flow of components takes place in a certain range of saturation values only. For the operation in this range of saturation values, the notion of effective saturations of the components The basic filtration law (Darcy rule) establishes a relationship between the flow velocity vector and the pressure drop that causes the filtration flow. The flow velocity, for each component is written (regardless of the gravity force) in the following form: where K is the porous medium characteristic called "absolute permeability", which is independent of the mixture properties; (4) Equation (4) will be further used for approximation by the SPH method.

2.Equations of state for liquid and gas
The paper poses the problem of simulating a three-component mixture of gas, water, and oil with regard to their compressibility. For convenience, the following index notations are further used for these components: for water, n for light oil (naphtha), and g for gas. The initial values of saturation, temperature, and pressure are specified for each component. Densities are calculated from the corresponding equations of state of the components depending on the pressure and temperature. The water and oil components are considered to be weakly compressible and linearly dependent on the temperature and pressure gradients: where = , and 0 is the known density value of the k-th component, corresponding to pressure 0 and temperature 0 . As for gas, the equation of state of an ideal gas is assumed to be valid: In the simulation the physicochemical parameters of the mixture, such as relative phase permeabilities and dynamic viscosities of the components are used. Thus, equations (4), (5), and (6) form a closed system which is solved with respect to pressures, densities, and saturations of the components.

3.Equation approximation with SPH method
The smoothed particle hydrodynamics (SPH) method [2] is based on smoothed integral interpolation of any scalar, or vector quantity specified in space. Smoothed interpolation ( ) of some scalar quantity ( ) is expressed as an integral: where ( − ′ , ℎ) is an even function also called an "interpolation kernel". Figure 1 graphically represents function W(r, h) (solid line) and its derivative, W'(r, h) (dashed line). The computational domain in the problem is partitioned into N elements with coordinates , where = 1, . Each element is hereinafter referred to as a smoothed particle. Passing from integration to summation in (7) and replacing volume element ′ by the volume of smoothed particle, we obtain the basic formula for the difference approximation of an arbitrary function ( ) in a given point : Accordingly, the derivative of function ( ) with respect to has the following form: Equations (8) and (9) The right-hand side of equation (10) contains the sum of three Laplace differential operators which should be approximated with SPH method.

4.Time integration
To calculate thermodynamic parameters, an explicit difference scheme is used in the next time step. Let the right-hand side of equation (4) be we obtain the following differential equation: The first-order difference scheme for the calculation of values at the next time layer for equation (14) looks like +1 +1 = + ⋅

5.Test computations
The initial conditions for tests are taken from [3]. In all problems, the rock porosity, density, and absolute permeability are considered invariable in the whole volume, the residual saturation values are assumed to equal zero, the computational domain size is 1 m × 1 m × 1 m. One-dimensional and twodimensional problems are solved in Cartesian coordinates. The mixture consists of gas, water, and light oil. For each component, the initial values of all dynamic parameterssaturation, temperature, and pressure -are considered to be specified. The medium is isotropic in the simulated volume. It is assumed that liquid phases are weakly compressible, gas is ideal, the temperature and pressure are the same for all mixture components. The temperature is considered invariable in the course of simulations. There is no gravity. Figure 2 graphically represents computation results for a 1D test from paper [3]. These computations were carried out using an explicit difference scheme on a regular and space-uniform computational grid of 50 nodes.  The results obtained are visually the same as those in paper [3].

Two-dimensional test on the gas injection under pressure
The physicochemical parameters and initial conditions for the computational domain in the two-dimensional test are the same as in the one-dimensional test described above.

6.Conclusion
In the numerical simulation of a real physical flow process (for example, oil production in an underground oil field) it is important to construct a three-dimensional computational grid in which the cell size in the oil well vicinity equals the typical size of the wellhead diameter (15-30 cm) and the cell size