Method SPH for numerical solution of filtering problems: including capillary forces in fluid flow simulation

This paper explores the theory of filtration within porous media, focusing on multicomponent mixtures’ dynamics and their numerical modeling. Key equations governing filtration, including Darcy’s law and mass conservation, are introduced. The concept of smoothed particle hydrodynamics (SPH) is presented as a powerful numerical technique for approximating these equations. SPH’s integral interpolation approach accommodates complex fluid-solid interactions and is applied to one-dimensional and two-dimensional test cases. Results demonstrate SPH’s effectiveness in capturing intricate dynamics and aligning with traditional methodologies. This study contributes to understanding fluid dynamics in porous media and underscores SPH’s potential for accurate simulations in diverse applications, from hydrocarbon reservoirs to environmental systems. The first section in your paper.


Introduction
In the realm of fluid dynamics and porous media, understanding and accurately modeling the intricate dynamics of multicomponent mixtures in porous media is of paramount importance.The filtration of mixtures in porous media constitutes a critical aspect of various fields, from hydrocarbon reservoir engineering to environmental remediation.This paper delves into the foundational equations and methodologies underpinning the theory of filtration within porous media, focusing on the core concepts of the smoothed particle hydrodynamics (SPH) method.
Filtration, in this context, signifies the intricate flow of multicomponent mixtures through porous bodies, where each component retains distinct physical and chemical attributes.Such phenomena are ubiquitous in nature, spanning scenarios as diverse as hydrocarbon migration in geological formations and groundwater flow through porous aquifers.By harnessing the principles of fluid flow, thermodynamics, and material behavior, researchers and engineers endeavor to comprehend and predict the intricate interplay of pressure, temperature, and saturation across different components within these porous systems.
The fundamental equations governing such processes are paramount in constructing reliable numerical models that simulate real-world scenarios.In this pursuit, the theory of filtration introduces key concepts like porosity, permeability, and relative phase permeabilities, all of which play crucial roles in the dynamics of mixtures within porous media.Darcy's law, an analogue to Fourier's law in heat conduction theory, serves as the cornerstone in linking filtration velocity with pressure gradients, setting the stage for numerical approximations.

Basic Equations of Mixture Filtration Theory in Porous Media
This article delves into the numerical modeling of multicomponent medium dynamics, building upon the concepts from [1].The premise is that the continuous medium comprises pores, cracks, and faults with dimensions significantly smaller than the overall physical size of the problem.The porosity of the medium quantifies its porous nature by relating pore volume to total volume: Where is  porosity coefficient,   is the pore volume and  is the total volume of a given element of the medium Filtration refers to the flow of multicomponent mixtures through porous structures.Here, a component denotes a segment of the mixture with distinct physical and chemical properties, maintaining its phase during simulation.
Key attributes of the system's motion encompass pressure, temperature, and component saturation.The saturation   of the k-th of component's pore space represents the fraction of pore volume occupied by that component in a unit volume: The sum of all component saturations in a volume element sums up to 1.Each component has a minimum saturation known as residual saturation (denoted as the subscript ).
Component flow is constrained within a saturation range.Effective saturations are introduced to operate within this range.The effective saturation of the k-th component can be expressed as: Fundamental to filtration theory is Darcy's law, relating filtration velocity and pressure drop.Similar to Fourier's law in heat conduction, Darcy's law governs filtration.The Darcy equation for the k-th component, neglecting gravity, is: Here, K -is absolute permeability, independent of the properties of the mixture;   -is the coefficient of dynamic viscosity of the k-th component   =   ();   =   (  ) -is the relative phase permeability of the k-th component (may depend on the saturation of other components);   -is a pressure of the k-th component.Darcy's law has applicability limits, as detailed in [2].
Another cornerstone in filtration theory is mass conservation.The mass conservation equation for the k-th component reads: Here   is the density of the k-th component.Substituting (1) into expression (2), we derive the core equation of filtration theory: Equation ( 3) forms the basis for approximation through the SPH method.

Equations of State for Liquids
This article addresses modeling two-component mixture-water and oil-accounting for their compressibility.Subscripts w, n signify water and light oil (naphtha), respectively.Initial values of saturation, temperature, and pressure are provided for each component.Density is calculated using component-specific equations of state dependent on pressure and temperature.Water and oil exhibit weak compressibility, linked linearly to temperature and pressure drop: Throughout modeling, physicochemical parameters like relative phase permeabilities and dynamic viscosities guide the process.Stone's approximation [2] determines relative phase permeabilities, represented analytically.Water's relative phase permeability is expressed as: Dynamic viscosities depend solely on temperature.For water, the dynamic viscosity follows: Hence, equations ( 3) and ( 4) form a closed system, solved for pressures, densities, and component saturations.permeability can be very low or even close to zero due to capillary forces, especially at low pressures.

2.2.Capillary forces
Capillary forces play a significant role in the behavior of liquids in porous media, particularly in very narrow pores and capillaries.
Capillary forces arise from the interaction between the liquid's surface and the porous structure of the material.In narrow pores and capillaries, such as in clayey rocks or thin fissures, capillary forces can dominate over pressure and inertia forces.At low pressures, the liquid can be "trapped" in these narrow pores due to strong capillary effects, resulting in very low permeability.
This effect is particularly important when considering low-permeability rocks, such as clay or sandstones, where pores can be extremely narrow and possess high capillary forces.Thus, permeability can be significantly reduced at low pressures due to the dominance of capillary forces.

3.Approximating Equations Using Smoothed Particle Hydrodynamics (SPH)
Historically, the smoothed particle method (SPH) was introduced for astrophysical problems like galaxy scattering.Since the 1990s, it has found extensive use in simulating hydrodynamic problems [4], involving high velocities and deformations in continuous media.Presently, SPH is applied across diverse physical fields.In [5], an algorithm for numerically solving heat conduction problems via SPH laid the foundation for utilizing the SPH algorithm in filtration theoryc [7].This stems from the mathematical similarity between the differential equations of filtration and heat conduction.
The smoothed particle hydrodynamics (SPH) method [3] revolves around smoothed integral interpolation of scalar or vector quantities in space.Smoothed interpolation for a scalar value is expressed as an integral: Here, ( −  < , ℎ) is an even interpolation kernel function distributed around a point r with a characteristic size h.The integral of the smoothing function across space sums to 1: As h approaches zero, W(r,h) tends to a delta function.Figure 2 depicts the graphical representations of W(r,h) (solid line) and its derivative W'(r,h) (dashed line).
The problem area is divided into N elements with coordinates   , where i = 1, 2, ..., N.
Each element is referred to as a smoothed particle.Shifting from integration to summation in equation ( 5) and replacing the volume element with the volume of a smoothed particle, we derive the fundamental formula for approximating an arbitrary function ( % ) at a specific point in space: Similarly, the derivative of function ( % ) with respect to   takes the form: Equations ( 6) and ( 7) serve as cornerstones in smoothed particle method theory.

3.1.Smoothed Particle Hydrodynamics for Laplace's Equation
The right-hand side of equation ( 3) must be adapted for subsequent approximation through the smoothed particle method.By denoting  = − ) * , the right-hand side is transformed identically: The transformed identity (8) involves a sum of three Laplace differential operators, to be numerically approximated using SPH.For the Laplace operator's numerical approximation, scalar A(r) s expanded into a Taylor series at adjacent points

𝐴(𝑟
Ignoring terms of fourth and higher orders, both sides of the equation are multiplied by the expression:   ( ?−  % ) ⋅  ?( ?−  % , ℎ) | ?−  % | A Continuing with the process and integrating over   , we arrive at an integral representation for the Laplace operator: By switching from integration to summation and replacing the volume element with the volume of the smoothed particle   : Substituting this approximation for Laplace operators, after simple algebraic simplification, the identity transforms into the form: Thus, equation (3) takes its final difference form: Note that the first derivative of the smoothing function   is used on the right-hand side of the difference filtration equation.

Time Integration
For numerical calculation of thermodynamic parameters in subsequent time steps, an explicit difference scheme is employed.By denoting   as the right-hand side of equation ( 3), we derive the following differential equation: The first-order accuracy difference scheme for calculating quantities at the next time step for this equation is as follows: Here, the superscript n denotes the quantity's value at time layer  =  ( , and the time integration step  is chosen based on the stability condition of the difference scheme when solving the Laplace equation  = ℎ A , where c<1.
Consequently, within each smoothed particle in space, a system of seven equations ( 3) and ( 4) is solved for seven unknown quantities.Solving this system at each point, new pressure, density, and component saturations are calculated in the next time step  =  (I' .Notably, the system of seven equations can be simplified to solving a cubic equation for pressure, assuming the incompressibility of water.Similarly, assuming the incompressibility of both water and oil, the system reduces to solving a linear equation for pressure.

Boundary Conditions
In the analyzed problems, both first-and second-kind boundary conditions are utilized to define pressure at the computational domain's boundaries.Pressure can be explicitly set or determined from boundary flow conditions.Saturations can also be explicitly set at the domain's edge.

Test Calculations
Following the outlined difference scheme, a prototype computer program was developed in the Swift language (chosen due to the Mac OS operating system) for conducting test calculations.
Let's consider a two-dimensional separation problem of a mixture containing water and oil.The objective is to separate the mixture into distinct oil and water fractions.The rectangular domain is filled with porous rocks of two types.The permeability of the first substance is constant and independent of pressure and fluids.The permeability of the second substance for water also remains unaffected by pressure.However, the permeability of the second substance for oil varies with pressure, as illustrated in Figure 1.
The mixture consisted of water and light oil.The oil density was lower than that of water.Initial values for dynamic parameters-saturation, temperature, and pressure-were given for each component.The medium was treated as isotropic throughout the volume, with weakly compressible liquid phases and uniform temperature and pressure across all mixture components.Gravity was absent.
Test calculations utilized the following values for parameters in equations ( 4      Result of calculations with radius M2=1.9 is shown in Figure 6.We see mixed water and oil flow and the same pressure in the fragment M2.Oil permeability in M2 is not zero at high pressures.

Conclusion
Through extensive investigations, a differential method utilizing smoothed particles has been developed and validated for solving the challenges posed by three-dimensional filtration of water, oil, and gas within heterogeneous porous media.Building upon this algorithm, a prototype software package has been crafted to model multidimensional filtration problems.Basic equations governing the theory of water, gas, and oil mixture filtration within porous and fractured media have been meticulously described and mathematically adapted, enabling their numerical approximation through smoothed particle theory's differential equations.
A notable advantage of the proposed method is its utilization of the smoothed function (interpolation kernel) and its first derivative for solving Laplace's differential equation.This reduction in differentiation complexity to the first derivative markedly enhances the numerical solution's quality while diminishing computational time.Nonetheless, the method's complexity arises when applying an implicit approach for differentiating the obtained equations over the temporal parameter.Addressing the inversion of the computational matrix for an irregular node set, a prerequisite for solving equations through an implicit method, necessitates substantial computational resources.
The isolated structure of the smoothed particle system permits the incorporation of new physical parameters and mathematical models.This adaptability extends to modifying equations of state for subsurface mixtures and algorithms for anisotropic behavior in porous media.This flexibility marks SPH as a versatile tool for addressing diverse challenges in complex filtration scenarios, ranging from hydrocarbon reservoir management to environmental studies.
In a two-dimensional test problem of a water and oil mixture, the objective was to achieve separation into distinct oil and water fractions.The test results indicate that increasing the area size leads to flow blockage and elevated oil pressure beyond the critical point.As permeability becomes non-zero, oil begins to flow into the hole.This suggests that the system's behavior is influenced by the interplay between area size, pressure, and permeability, impacting the flow dynamics of the mixture.

Figure 1 .
Figure 1.Oil permeability reduced at low pressures.

Figure 2 .
Figure 2. The smoothed function  and its derivative  < .

Figure 3
Figure3shows the initial geometry and bound condition.The rectangular shape M1 is filled with the first substance.In the center of the rectangular domain, there is a circle M2 filled with the second substance.Within the circle, there is a low-pressure outlet H where the liquid flow is discharged.

Figure 4
Figure4shows the initial geometry of SPH particles.There are 25 points in horizontal and 10 points in vertical direction.There are solid bound conditions in top and bottom lines.

Figure 5 .
Figure 5. Pressure profiles with radius 1 m

Figure 6 .
Figure 6.Pressure profiles with radius 1.9 m

Figure 7
Figure 7 shows the dependence of oil flow on the radius of the M2.A large size of the area M2 leads to flow blockage and an increase in pressure in the oil beyond the critical point.The permeability of the medium M2 becomes non-zero, and oil starts flowing into the hole.The critical size of M2 area is approximately r=1.4 m in this modelling geometry.

Figure 7 .
Figure 7. Dependence of oil flow on the radius of the M2.