Mathematical modeling of elastic soil acid treatment

The study is devoted to an initial boundary value problem describing the process of cleaning the bottomhole zone of an oil well with acid solution. It is assumed that the soil frame is an elastic solid body and the pore space has double porosity. The physical process is described at the microscopic level by the Stokes equations for the liquid component, the diffusion-convection equation for the concentrations of acid and chemical reaction products, and the Lame equations for the solid frame. Due to soil dissolution, the pore space has an unknown (free) boundary. The lattice Boltzmann method (LBM), the immersed boundary method (IBM) and the finite element method (FEM) are used for computer simulation. Finally, discretization of equations and results of numerical solutions are presented.


Introduction
The bottomhole cleaning process (see figure 1) is very important to increase oil production. Oil reservoirs are complex geologically heterogeneous objects. Heterogeneity means that the medium properties change in space. The analysis of wells and cores shows that the geological properties (porosity, permeability, etc.) of oil reservoirs are heterogeneous even within one field. Often, insufficient consideration of heterogeneity at the planning stage becomes evident too late, when the acid solution pumped into the soil through injection wells is far from the intended place. In addition, the concentration of injected acid and injection modes play an important role. Therefore, understanding the dynamics of fluids in heterogeneous porous media and the mechanism of rock dissolution by acids is of fundamental importance for the effective management of oil production. This can be achieved by creating a hydrodynamic simulator of the bottomhole zone based on the corresponding mathematical model, which allows to optimize the considered technological process.
There is a wide range of mathematical models describing the process of rock dissolution at the macroscopic level (see [1][2][3], and references therein), where typical dimensions are meters or tens of meters. At the macroscopic level, at each point of the continuous medium, there is both a solid frame (rock) and a liquid in the pores of this frame. All such models are based on the same principle. Fluid dynamics are usually described by the Darcy system of filtration equations or some of its modification. Equations describing the migration of acid and chemical reaction products are simply postulated and are modified diffusion-convection equations for the corresponding concentrations. The main thing in these postulates is the form of coefficients in equations, because the main mechanism of the physical process is concentrated on the free boundary between the pore space and the solid frame. Moreover, during the process, the geometry of the pore space changes. All these fundamentally important changes occur at the microscopic level, corresponding to the average size of pores or cracks in rocks, while any of the proposed macroscopic models operates on completely different scales.
The exact formulation of rock dissolution problem in one porosity medium for an absolutely rigid solid body at the microscopic level was first presented in the monograph [4], but mathematically rigorous results on the existence of any solution were absent there.
This paper deals with a more complex problem, where pores of size ε = l/L and cracks of size δ = ε ν (0 < ν < 1) are taken into account in elastic soil. Here l is the characteristic pore size and L is the characteristic size of the considered physical region.
Pores are modeled by the system of cylindrical channels of radius ε r p (x, t), 0 < r p < 1/2, where its axes are parallel to the axes coordinates located at a distance ε from each other. Cracks are modeled by sphere of radius δ r c (x, t), 0 < r c < 1/2, located at a distance δ from each other. The motion of viscous fluid in an elastic periodic skeleton with double porosity was previously considered by A.M. Meirmanov [5].
Here the physical process at the microscopic level is described by the laws of continuum mechanics. These are the Stokes equations for the velocity v ε and the pressure p ε of a viscous incompressible fluid in domain Ω ε f (t) with dimensionless viscosity κ µ , Lame equations for displacement w ε of elastic solid frame Ω ε s (t), diffusion-convection equation for the acid concentration c ε and the conditions on a strong discontinuity (the required free boundary) Γ ε (t) "elastic frame-pore space" [6]. There are no inertial terms in the Stokes and Lame equations, as it has a sufficiently long duration in time (the fluid filtration rate is several meters per year) The problem is supplemented by the boundary conditions on the free boundary separating the pore space and the soil frame, which follow from the conservation laws of classical mechanics in its integral form, boundary conditions at given boundaries, initial conditions, and an additional condition on the free boundary following from the laws of theoretical chemistry [7] and allowing to define this boundary. Here d ε n is the velocity Γ ε (t) in the direction of unit normal n to the boundary Γ ε (t), c ε is the concentration of acid in the liquid, κ is a given constant.

Materials and methods
The physical process is considered in a bounded domain Ω (see figure 2) with piecewise Lipschitz boundary S,S = ∂Ω, The domain Ω ε f (t) ⊂ Ω describes the region occupied by liquid, and Ω ε s (t) ⊂ Ω denotes elastic frame (see figure 2). The boundary S 0 is impermeable to fluid in the pore space, the boundary S 1 models the injection well. In dimensionless variables the dynamics of acid solution in an unknown domain Ω ε f (t) (fluid domain) is described by the system of differential equations for the velocity v ε and the pressure p ε of continuous medium, p a is atmospheric pressure, c ε is the reagent concentration, I is a unit matrix.
In elastic frame Ω ε s (t) the displacements of continuous medium w ε are described by the Lame There are boundary condition (1) and conditions on the free boundary expressing the laws of conservation of mass [8].
In (2)-(9) v ε n = v · n is the normal component of the fluid velocity vector v, n is the unit normal vector to the free boundary Γ ε (t), outer with respect to the region Ω ε f (t), κ c , κ i , β and β i (i = 1, ..., m) are given positive constants, v ε and p ε are the velocity and pressure in the liquid, c ε i are concentrations of chemical reactions products. Condition (6) simplifies the boundary condition (8) On the given boundaries, consisting of the injection well S 1 , the lateral boundary S 2 and the impermeable boundary S 0 , the following conditions are satisfied In (8)-(16) n is the unit normal vector to the boundaries S 0 , The problem is completed by initial conditions For simplicity, assume p 0 , p 1 and c 1 are given positive constants. Note that the problem (1), (6)-(8), (10)- (14), (17) is solved independently of defining functions c i (i = 1, 2, ..., m). The latter are determined after finding the functions {v ε , p ε , c ε , r}.
All functions of the form Φ(x, t, y, z) are considered 1-periodic in variables y and z respectively, where (x, t) ∈ Ω and The structure of the pore space Y f (r p ) is given by the desired function r p (x, t), figure 3a). The structure of the crack space is given by the required function r c (x, t), 0 < r c (x, t) < 1/2 as Z f (r c ) = {z ∈ Z : z 2 figure 3b). With a fixed ε > 0 domain Ω ε f (t) filled with fluid and elastic frame Ω ε s (t) are determined by the characteristic function φ(x, t, y, z) as follows n (x, t) we denote velocity of the free boundary in the direction of the normal n(x, t) to the surface Γ ε (t) at the point x ∈ Γ ε (t).

Results and discussions
In this study, the lattice Boltzmann method (LBM) was used to solve the system of equations (1)-(4). The method describes the behavior of particles collection and is based on the Boltzmann equation, which describes the gas dynamics on the mesoscopic scale [9] ∂f ∂t where f is the particle density at the point x of the time t, Θ is the collision operator that determines the local distribution after the collision of particles. Using LBM and equation (2) with respect to time, rewrite our equations for velocities as follows In equations (22)-(25), ∆x and ∆t are steps in space and time, f i is the velocity distribution function (see figure 4), f eq i is the equilibrium distribution function, and τ is the dimensionless relaxation time  Here 2 is the number of spatial dimensions (x and y), and 9 is the number of speeds . After these two steps are completed, the distribution functions f i are calculated for the current time step, and then the fluid velocity is found from Similar transformations are carried out for displacements of the elastic soil frame described by the system of equations (5).
The product flux of chemical reactions can be described using the convection-diffusion equation [10] ∂C j ∂t where C j is the concentration of the products of chemical reactions j, κ c j is the dimensionless diffusion coefficient determined for each reaction product. It is worth to note the absence of the term responsible for the chemical reaction, in view of the presence of the boundary conditions (8) and (9) at the "liquid-solid" interface.
The LB-analogue of the equation (28) will be as follows where Solving equations (29)-(31), find the concentration in equation (28) Since the time scale of fluid flow is much smaller than other processes, the velocity reaches steady state during the selected time step. The equilibrium distribution functions are calculated using the equation (31), then the collision step (see figure 5a) is performed according to the following equation The concentration distribution functions near the "liquid-solid" interface (see figure 5b) will have the following form Since the time scale of fluid flow is much smaller than other processes, the velocity reaches steady state during the selected time step. The equilibrium distribution functions are calculated using the equation (31), then the collision step is performed according to the following equation The geometry of elastic solid frame is updated using the rules detailed in [11]. Before performing the propagation stage, to apply the rock dissolution effect, first update the distribution function and then calculate the propagation stage After that, solving (32) obtain the concentrations of the products of chemical reactions and proceed to the next time step.
As stated earlier, the multiple relaxation in time (MRT) scheme [12] is used to improve the stability of calculations. This scheme uses the transformation matrix M . Thus, the collision operator will be where S = diag(0, s e , s , 0, s q , 0, s q , s v , s v ) is a diagonal matrix. The parameters s v and the relaxation time τ are related as To improve accuracy, it is extremely important to choose the correct relaxation parameter values. Thus the establishment of so-called viscous relationship between s q and s v [13] is used to reduce the nonphysical behavior of fluid near the free boundary in the course of numerical solution For s e and s s e = s = s v .

Immersed boundary method (IBM)
The immersed boundary method is used to track the free boundary behavior. The main idea of the IB method is that the flow field is influenced by the force distribution created by the free boundary of the solid body. IBM uses independent meshes for fluid and solid sampling. The fluid is discretized by a set of Eulerian points, which are in fact fixed regular Cartesian lattice points, and the body boundary is discretized by a set of Lagrangian markers. Let us note the importance of coefficient β in the condition on the free boundary. It affects the dissolution rate of elastic solid frame. Figure 8 shows the concentration field of the reagent with β = [1; 10] at fixed time moment t=2.7 * 10 5 .  Finally, figure 9 shows the position of the free boundary at fixed time moment t=1 * 10 6 . The gray borders denote the geometry of elastic soil frame before acid exposure.

Conclusion
In this study, a correct mathematical model for cleaning the bottomhole zone of oil wells was derived based on the fundamental laws of continuum mechanics and theoretical chemistry, which accurately describe the physical processes under consideration at the microscopic level (on the scale of pores and cracks). The problem of acid treatment of a non-periodic elastic solid frame with double porosity was considered. Numerical experiments of the problem are carried out. The influence of acid concentration in the injection well and coefficients in the condition on the common boundary "liquid-elastic solid body" is investigated. On each temporarily layer, the border was adjusted twice. First, as a result of fluid action on an elastic body, then as a result of acid dissolution. As a result, an increase in the products of the chemical reaction can be noted over time at a given value of the concentration of the reagent at the inlet boundary. Also revealed the influence of the coefficient β in the condition on the free boundary. With its increase, the dissolution process proceeds faster.