Abstract
In the process of reservoir exploitation, the range of rock stress field changes is much larger than the reservoir area where the seepage occurs; so the amount of calculation of the existing seepage stress full coupling method is often very huge. Based on the theory of seepage mechanics and rock mechanics, a coupling analysis model of reservoir multiphase seepage and stress is established, and a numerical solution is established by using finite element and finite difference methods. The evolution law of the seepage field and stress field and the change of rock mechanics parameters can be studied, with emphasis on the readjustment of rock stress distribution and its deformation characteristics under the influence of the seepage field. At the same time, to improve the calculation efficiency of numerical simulation, considering the difference between the seepage calculation area and stress calculation area, the finite element method is improved. The percolation area of the reservoir is calculated with a fine grid to obtain a more accurate distribution of underground fluid. The coarse grid is used to calculate the stress calculation area and to reduce the calculation time. The mechanical equilibrium equation in the fully coupled theory is discretized on the coarse grid by the finite element method. The mass conservation equation is discretized on the fine grid by the finite volume method. The numerical simulation of Terzaghi’s one-dimensional consolidation problem and Mandel’s two-dimensional consolidation problem shows that the calculation results of this method are in good agreement with the analytical solution. Through the numerical calculation of a two-dimensional single-phase flow single well production problem and a three-dimensional two-phase flow five-point well pattern production problem, the influence of the seepage stress coupling effect in reservoir numerical simulation is analyzed.
1. Introduction
The thermal fluid-solid coupling theory of geotechnical media is the product of the mutual penetration and intersection of seepage mechanics, rock mechanics, and heat transfer. The research in this field has involved the fields of seepage and control of water conservancy and hydropower projects, reservoir-induced earthquakes, deep burial treatment of nuclear waste, coal and gas outbursts, stability of rock slopes and dam foundations, development of oil and natural gas and heat energy, utilization of groundwater resources, and so on [1, 2].
In the process of oil and gas production, with the continuous production of oil and gas, the pore pressure of the reservoir is reduced and the solid stress is redistributed, which leads to the deformation of the rock framework of the reservoir and changes the physical parameters of the reservoir, especially the porosity, permeability, and pore compressibility, which in turn affect the flow of reservoir fluid in the pore space. Therefore, the porosity, permeability, and rock deformation ability of oil and gas reservoirs are directly related to oil and gas recovery, which should be studied. In the process of drilling, oil and gas wells are soaked by fluids inside and outside the well, which directly affects the stability of the good wall. Similarly, in the mining process, due to the flow and erosion of fluid, the properties and stress of the rock skeleton around the shaft wall will change, destroying the skeleton, which leads to a large amount of sand production. Therefore, the coupling between fluid and rock must be considered in the wellbore stability analysis and sand production analysis. In the process of reservoir exploitation, due to the extraction of a large number of fluids (oil, gas, and water), the pressure of the oil (gas and water) layer is reduced, resulting in the deformation, compaction, and settlement of the overlying strata, which will bring serious consequences, such as borehole collapse, casing deformation, and damage, which is particularly important for inclined wells. In the process of water injection or polymer injection displacement production, the high pressure of the displacement fluid not only provides the driving force for the flow of oil and gas but also expands the pore space and improves the permeability of the reservoir to achieve the purpose of increasing production. Therefore, studying the fluid-solid coupling in the process of displacement is a problem that must be solved to improve oil and gas recovery [3].
As for the coupling problem of seepage stress in porous media, a more systematic and comprehensive theoretical and numerical calculation method has been established in the field of geotechnical engineering. Terzaghi [4] first introduced the concept of effective stress and proposed a one-dimensional consolidation model of rock and soil. Biot [5, 6] established a perfect three-dimensional consolidation theory and wave propagation theory with solid-liquid corresponding force balance equation and continuity equation as the main governing equations. In terms of numerical calculation, the seepage stress coupling problem of porous media is mostly studied by finite element method [7–15], such as the discrete u-p method and u-UW method of the Biot equation finite element method proposed by Liu et al. [16]. The field of geotechnical engineering pays more attention to the accurate calculation of mesh geotechnical stress and deformation. The finite element method is very suitable for the seepage stress coupling problem. The field of oil and gas reservoir engineering focuses more on the accurate calculation of the velocity and flow of oil-gas-water fluid between reservoir grids. Because the finite volume method is not limited by the orthogonal grid of the finite difference method in terms of geometric flexibility, at the same time, its physical meaning is clear, it has strict local conservation, and it can easily and accurately deal with the discontinuity of flow, so at present, the finite volume method is mostly used for numerical simulation of oil and gas reservoirs [17–20]. Akai and Tamura [21] proposed a mixed finite element finite difference method to solve the Biot equation of porous media. Subsequently, this method was developed and extended to the numerical calculation of geotechnical consolidation and saturated soil liquefaction [21–31]. On this basis, Yuan and Xiaowei [26] proposed and studied a mixed finite element finite volume method for seepage stress coupling problems.
However, different from ordinary geotechnical engineering problems, the spatial scale and time scale of reservoir numerical simulation are often large, and the seepage calculation area and stress calculation area do not coincide in the seepage stress coupling study [27]. In the actual reservoir exploitation, the underground fluid seepage occurs in the reservoir range, and the reservoir area can be selected as the seepage calculation and stress calculation area. The change of reservoir rock stress field in the production process is not limited to the scope of the reservoir but also affects the upper overburden and lower stratum of the reservoir. Therefore, the stress calculation area should generally be much larger than the reservoir range where the seepage calculation area is located. If the stress calculation area is limited to the reservoir seepage calculation area, the result of the seepage stress coupling calculation will produce a certain error [28–31]. Because reservoir numerical simulation pays more attention to the accuracy of seepage flow velocity and flow calculation, the grid in the reservoir should ensure a certain degree of fineness. In the original finite element finite volume hybrid method, the same mesh is used when the finite element method discretizes the stress balance equation and the finite volume method discretizes the mass conservation equation [26]. Applying it to the reservoir seepage stress coupling calculation will result in a very large number of grids. To improve the calculation efficiency of the reservoir seepage stress full coupling problem, this paper proposes an improved finite element finite volume hybrid method. This method uses the finite element method to discretize the stress balance equation on the coarse grid and the finite volume method to discretize the mass conservation equation on the fine grid, which not only ensures the calculation accuracy of the reservoir seepage field but also reduces the calculation burden caused by the solution of the stress field.
Therefore, the geological body is regarded as the aggregation of solid skeleton and pore fluid. The distribution characteristics and laws of oil and gas show that the migration and storage of pore fluid in porous media are closely related to the structure of the solid skeleton and rock properties. Therefore, the relationship between the solid skeleton system and the pore fluid system has become an important research content. Because the structure of a solid skeleton is controlled by tectonic stress, the study of the correlation between oil and gas migration and the stress field during tectonic activity provides a new perspective for the study of reservoir forming theory.
2. Basic Model of Porous Media
2.1. Displacement Mechanism of Tectonic Stress Field
According to the concept principle of effective stress (1), the tensor expression of effective stress in porous media is
Similarly, for the effective average stress in the fluid-solid system, there can also be the following relationship
The formula shows that the effective stress is inversely proportional to the fluid pressure without changing the total stress. The stress of pore fluid and solid skeleton depends on the characteristics of the medium itself, such as porosity, density, and plugging conditions. Assuming that the pore fluid is incompressible, under the condition of water permeability, the medium is squeezed by tectonic forces, and the skeleton deforms, resulting in a change of porosity. The pore fluid is squeezed out, and the skeleton flows to the area with low tectonic pressure or good permeability conditions. At this time, the stress borne by the skeleton is relatively larger. Under the condition of complete impermeability, the pore fluid is blocked in the medium, which increases the pore fluid pressure, and the tectonic stress will be mainly borne by the pore fluid. Therefore, the pore fluid pressure in porous media caused by tectonic stress can be expressed as [14]
c can be given according to medium parameters such as porosity and density. In addition to the average tectonic stress, the magnitude and direction of the maximum principal pressure stress will also affect the migration of pore fluid. It will make faults or fractures in the same direction in a relatively stretched state and make the vertical direction become a flow channel or make faults or structures in the vertical direction close, blocking the migration of oil and gas along the direction of the maximum principal pressure. At the same time, the flow intensity depends on the difference between the maximum principal pressure and the minimum principal pressure.
When the external structural force acts on the geological body, due to the uneven distribution of large and small structures and lithological differences in the rock mass, the differential distribution of structural stress leads to the deformation of the geological body, which changes the porosity and pore hydraulic pressure of rock masses with different lithology, resulting in the average stress difference, principal stress difference, or potential difference, driving the migration of oil and gas from high potential areas to low potential areas. When encountering suitable traps, it is possible to accumulate and form reservoirs.
2.2. Fluid-Solid Coupling Analysis of Porous Media
Under the condition of plane strain, the equilibrium equations of porous media under external force are
Among them,
According to the energy principle, the fluid potential can be simply expressed as
The strength of the force field driving oil and gas migration can be obtained from the fluid potential gradient
The fluid migration velocity can be calculated by Darcy’s law where .
The above formula shows that the seepage velocity of fluid is directly proportional to the pore pressure gradient. Through the fluid potential equation and Darcy’s law, the fluid potential field and seepage velocity of porous media caused by tectonic pressure can be obtained, and the dominant orientation and enrichment area of oil and gas migration can be analyzed.
3. Application of Fluid-Solid Coupling Theory in Reservoir
3.1. Reservoir Multiphase Fluid Seepage Equation
(1)Reservoir gas composition equation is (2)Reservoir oil composition equation is (3)Reservoir water component equation is where represents the mass of oil (steam and water) produced or injected into the formation per unit volume per unit time, kg/m3•s. are oil, gas, and water saturation, dimensionless. is the reservoir depth, . is the absolute permeability. , , and are the relative permeability of gas, oil, and water. , , and are volume coefficient of underground gas, oil, and water, dimensionless.
3.2. Deformation Field Equation of Reservoir Rock Mass
Due to the complexity of the specific geological environment, the selection of the constitutive model of reservoir rock mass is very complex, especially in the process of finite element calculation, the selection of different models may cause great differences. A large number of documents have compared the calculation results under different models. In this paper, the commonly used elastic-plastic model and the Drucker Prager yield criterion are used to obtain the mathematical equation of the deformation field of the reservoir rock mass. Equilibrium differential equation considering effective stress
Geometric equation of rock mass deformation
Elastoplastic constitutive equation
In Drucker Prager yield criterion , is the first invariant of the stress tensor, and is the second invariant of the stress tensor,
φ is the internal friction angle, and C is the cohesion. Then, there are
is the yield function vector, using the correlation flow rule.
3.3. Coupling Relationship between Fluid and Reservoir Rock Mass
Fluid not only has the mechanical effect of pore pressure on rock mass but also changes the original constitutive relationship of rock mass and changes the physical and mechanical properties of a rock mass. At the same time, due to the effect of external stress, the deformation of rock mass is caused, which changes the parameters of the seepage field, such as porosity and permeability of rock mass, and will inevitably lead to the corresponding change of seepage field.
The influence of the stress field on the seepage field is mainly reflected by porosity and permeability. A large number of experiments show that the relationship between porosity and permeability and effective pressure can be expressed in the following form: where is the porosity under the reference pressure. k is the permeability under the reference pressure, . and , they are the experimental parameters corresponding to the calculation of rock mass.
The influence of fluid seepage on the mechanical properties of rock mass can be expressed by the following formulas: where is the peak strength of rock mass, MPa. is the residual strength of rock mass, MPa. , , and are the corresponding coefficient, dimensionless. , , and are the peak strength, residual strength, and elastic modulus under reference pressure, MPa.
4. Analysis of Seepage Field
4.1. Definite Solution Condition of Seepage Field
(1)Constant pressure boundary means that the pressure at every point on the casing at the reservoir boundary is known or the bottom hole pressure of the well is known The above formula represents the given function of pressure at a point on the outer boundary G at time .(2)Constant flow boundary this boundary condition is the boundary value of the derivative of the given unknown quantity, that is, the derivative of the unknown quantity at the reservoir boundary or the well production is known where is the derivative of the pressure with respect to the outer normal direction of the boundary, and n is the normal direction. is a known function on the specified boundary.
4.2. Initial Condition
To solve the unsteady seepage problem, we also need to give the pressure and saturation distribution in the reservoir at the initial time, which is expressed as
4.3. Definite Solution Conditions of Rock Mass Deformation Field
The calculated deformation is mostly the first kind of boundary condition, that is, the surface force of rock and soil skeleton is known where is the directional derivative of the boundary. is the distribution function of surface force.
4.4. Seepage Stress Coupling Solution
The governing equations of the seepage stress coupling problem are composed of the equilibrium equation of geotechnical materials and the seepage equation of pore fluid. For the element body of the seepage stress coupling problem, the differential equilibrium equations of solid skeleton and pore fluid can be established, respectively, and the overall equilibrium equation can be obtained by adding the two equilibrium equations where is the average density of porous media. Superscript “s” refers to the solid skeleton and “f” refers to pore fluid. is porosity. is the displacement of solid skeleton; is the volume force; is the total stress of the unit body, which . The stress of solid skeleton and pore fluid, respectively. is the effective stress. is pore fluid pressure. is Kronecker symbol.
According to Terzaghi’s effective stress principle, when the Biot coefficient is taken as 1, there is the following relationship:
The seepage equation can be established according to the mass conservation equation of the pore fluid of the unit body where is the relative velocity of pore fluid relative to a solid skeleton. is the strain of a solid skeleton.
Assuming that the acceleration of the fluid relative to the solid skeleton is very small [16], that is , the relative velocity of the fluid relative to the solid skeleton can be expressed as where is permeability and is liquid viscosity.
The above formula can be regarded as Darcy’s law expression considering the deformation of solid skeleton. Compared with the original expression of Darcy’s law, the above formula adds a coupling term . The existence of this term makes the equilibrium equation and seepage equation have to be solved simultaneously, that is, a fully coupled calculation method is required.
5. Improvement of Finite Element Method
For reservoir exploitation, the influence range of the stress field is different from that of the seepage field. The solution of the seepage equation is limited to the reservoir area, and the calculation of the equilibrium equation should be carried out in a larger area including the whole reservoir, upper overburden, and lower stratum [27, 28]. In order to meet the requirements of numerical simulation of reservoir seepage stress problem, in view of the different characteristics of the calculation areas for the solution of equilibrium equation and seepage equation, the finite element method grid of discrete equilibrium equation and the finite volume method grid of discrete seepage equation is divided into grids with different degrees of fineness in the scope of the reservoir. As shown in Figure 1, Mesh 5-6-7-8 is a fine mesh element of the finite volume method. Mesh 1-2-3-4 is a coarse mesh element of the finite element method, and the finite volume method element is obtained by further subdividing the finite element method element. In order to facilitate the implementation, the fine mesh of the finite volume method can be automatically divided by connecting the M bisectors on the opposite sides of the finite element method mesh. Thus, one finite element coarse mesh contains M2 finite volume fine meshes.
According to the improved finite element method [26], for the balance equation (23), the finite element method is used to discretize the space domain on the coarse grid, and the following can be obtained: where is the node displacement increment of coarse mesh element. is the stiffness matrix. is the volume strain stiffness matrix,
is the average excess pore pressure of coarse mesh elements, where is the excess pore pressure of fine mesh elements. is the load item, is the increment of effective stress relative to the initial state.
For the seepage equation, the finite volume method is used to discretize the space domain on the fine grid, and the following can be obtained: where is a matrix composed of finite element normal shape function and the node displacement vector of fine mesh is obtained by interpolating the node displacement vector of coarse mesh. is a coupling term. is the flow terms related to the cell shape, where and are the projection length of the adjacent edges of two adjacent cells in the horizontal and vertical directions, the projection length of the centroid line of two adjacent cells in the horizontal and vertical directions, and are the length of the centroid line of two adjacent cells, and is the harmonic average of the permeability of two adjacent cells. is the excess pore pressure of the former unit, is the excess pore pressure of adjacent element .
6. Simulation Case Analysis
6.1. Terzaghi One-Dimensional Consolidation Problem
The improved coupling method is used to calculate Terzaghi’s one-dimensional elastic consolidation problem to verify the accuracy of this algorithm. The calculation model used is shown in Figure 2. The width of the calculation model is 5 m, and the height is 10 m. The top of the model is a drainage boundary; the other boundaries are not drained; the bottom is a zero displacement boundary; and the side displacement is constrained in the horizontal direction. When the vertically uniformly distributed load is applied on the top of the model, excess pore pressure will be generated inside the model. With the gradual dissipation of the excess pore pressure in the soil, the effective stress and vertical settlement displacement in the whole calculation area will gradually increase. When the excess pore pressure is completely dissipated, the settlement will reach the maximum.
For this one-dimensional consolidation problem, this method adopts different fineness densities for the finite element mesh of stress balance calculation and the finite volume mesh of pore fluid seepage calculation, as shown in Figure 3. The coarse mesh on the left is the mesh used for finite element method discretization, and the fine mesh on the right is the mesh used for finite volume method discretization. The parameters used in the calculation example are shown in Table 1.
Figures 4 and 5 , respectively, show the comparison between the numerical results and theoretical solutions of vertical displacement and pore pressure. It can be seen that the scattered data can be well-fitted with the theoretical solution, and the numerical value and change trend are relatively consistent. The calculation results of this method are in good agreement with the theoretical solution of Terzaghi’s one-dimensional consolidation problem, which verifies the accuracy of this algorithm.
6.2. Mandel’s Two-Dimensional Consolidation Problem
Mandel’s two-dimensional consolidation problem considers the porous medium sandwiched between two smooth (nonfrictional resistance), impermeable, and infinite rigid plates. Both sides are permeable boundaries, and the rest are impermeable boundaries. Apply 2F load on two plates, respectively, and consider the settlement of porous media. As shown in Figure 6, due to the symmetry of the model and the symmetry of the boundary conditions, only 1/4 of the model needs to be taken as the calculation area for research, and the calculation parameters used are shown in Table 2. Figures 7 and 8 show the numerical results and theoretical solutions of vertical and horizontal displacement of Mandel’s two-dimensional problem. It can be seen that although there are some errors between the numerical calculation results and the theoretical solution at the moment of applying the load, with the advance of time, the numerical calculation results are in good agreement with the theoretical solution, which further verifies the accuracy of this algorithm.
6.3. Calculation of Two-Dimensional Uncoupled Single-Phase Fluid Well
The model of the two-dimensional uncoupled single-phase fluid well calculation problem is shown in Figure 9. The calculation area is 2000 m long and 1500 m deep and is divided into six layers of rock strata from the surface down. The reservoir is located 1000 m underground with a thickness of 50 m, belonging to an anticline trap. The production well is located in the middle of the reservoir and is exploited at constant pressure with a bottom hole pressure of 3.5 MPa. Because the fluid flow in the nonreservoir area is not considered, only the finite element mesh is given, and the finer finite volume method mesh is used for the reservoir area. The specific parameters of each layer in the model are shown in Table 3.
The model used in the calculation is shown in Figure 10. Figures 11 and 12 show the numerical simulation results of daily production and cumulative production of production wells, respectively. In the figures, is the daily output, is the cumulative output as of that day, and is the cumulative output of the 800th day without considering the seepage stress coupling. It can be seen that after considering the seepage stress coupling effect, the daily output of the reservoir will decline faster in the early stage of production than when the seepage stress coupling effect is not considered and will enter the stable production stage earlier. The cumulative production curve is also significantly lower than that without considering the seepage stress coupling effect.
6.4. Fluid Density Calculation of Three-Dimensional Two-Phase Coupled Five-Point Well Pattern
Using the improved finite element method in this paper, the seepage stress coupling numerical simulation of a three-dimensional two-phase flow five-point well pattern production problem is carried out. Due to the symmetry of the five-point well pattern production problem model, only 1/4 of its area needs to be modeled and calculated. The whole calculation area is 1000 m long, 1000 m wide, and 1200 m deep. The stratum is divided into four layers, of which the reservoir is 800 m underground and 100 m thick. The initial porosity of the reservoir rock is 0.2, the initial permeability is 0.3 Darcy, the initial pore pressure is 12 MPa, and the initial oil saturation is 0.8. The 1/4 area of the five-point well pattern problem is exploited in the way of one injection and one production. Injection wells and production wells are arranged at two diagonally opposite corners of the reservoir area. The injection wells are injected at a constant pressure of 12 MPa, and the production wells are produced at a constant pressure of 9 MPa. A total of 20,000 days of production simulation has been carried out.
The influence of the seepage stress coupling effect on the calculation results of the oil production rate and water production rate is shown in Figure 11. Figure 12 shows the influence of the seepage stress coupling effect on the calculation results of cumulative oil production and cumulative water production. After considering the seepage stress coupling effect, the daily production and cumulative production of oil and water in production wells are significantly lower than that without seepage stress coupling. It can be seen that after considering the seepage stress coupling effect, the reservoir grid permeability decreases slightly near the injection well but greatly near the production well. At 20000 days of production, the permeability near the reservoir injection well decreases to 95% of the initial value, and the permeability near the production well decreases to 86% of the initial value. With the progress of production, the unit permeability in the whole reservoir is gradually decreasing, which is the reason why the daily production and cumulative production ratio in Figures 13 and 14 are significantly reduced when the seepage stress coupling effect is not considered.
7. Conclusion
Based on the theory of seepage mechanics and rock mechanics, this paper establishes a coupling analysis model of multiphase seepage and stress in the reservoir and uses finite element and finite difference methods to establish a numerical solution, which can study the evolution law of seepage field and stress field and the change of rock mechanics parameters and focuses on the readjustment of rock stress distribution and its deformation characteristics under the influence of seepage field. (1)Based on the theory of continuum mechanics, the mechanism of tectonic stress displacement is analyzed, and the rationality of taking the average stress in porous media as the coupling between the solid skeleton and pore fluid is expounded(2)The accuracy and effectiveness of this method are verified by using the finite element numerical simulation method. Examples of the Terzaghi one-dimensional consolidation problem and Mandel two-dimensional consolidation problem with theoretical analytical solutions are calculated, respectively, and the numerical results are in good agreement with the analytical solutions(3)The calculation of two-dimensional single-phase flow single well constant pressure production and three-dimensional two-phase flow five-point well pattern production shows that the daily production and cumulative production of production wells with seepage stress coupling are significantly lower than those without coupling. When considering the coupling of seepage and stress, the unit permeability in the whole reservoir is gradually decreasing with the progress of production. The permeability decreases slightly near the injection well and greatly near the production well
Data Availability
All data included in this study are available upon request by contact with the corresponding author.
Conflicts of Interest
The authors declare no conflict of interest.
Acknowledgments
We gratefully acknowledge the research project of the Exploration and Development Research Institute of CNPC Changqing branch (2020-132) financial support. We gratefully acknowledge the School of Geosciences, Yangtze University.