Numerical Simulation of Liquid-Solid Coupling in Multi-Angle Fractures in Pressure-Sensitive Reservoirs

This study clarifies the seepage characteristics of complex fractured pressure-sensitive reservoirs, and addresses a common technological problem, that is the alteration of the permeability degree of the reservoir bed (known to be responsible for changes in the direction and velocity of fluid flows between wells). On the basis of a new pressuresensitive equation that considers the fracture directional pressure-sensitive effect, an oil-gas-water three-phase seepage mathematical model is introduced, which can be applied to pressure-sensitive, full-tensor permeability, ultralow-permeability reservoirs with fracture-induced anisotropy. Accordingly, numerical simulations are conducted to explore the seepage laws for ultralow-permeability reservoirs. The results show that element patterns have the highest recovery percentage under a fracture angle of 45°. Accounting for the pressure-sensitive effect produces a decrease in the recovery percentage. Several patterns are considered: inverted fivesevenand nine-spot patterns and a cross-row well pattern. Finally, two strategies are introduced to counteract the rotation of the direction of the principal permeability due to the fracture directional pressure-sensitive effect.


Introduction
Based on Warren and Root's single-phase seepage equation, in 1976, Kazemi et al. [1] developed a tool for simulating the multi-well single-phase and two-phase seepage of fractured reservoirs. In 1986, Nakornthap et al. [2] built a numerical simulation model capable of simulating heterogeneous isotropic matrices and heterogeneous anisotropy fracture systems, and characterized fracture permeability in tensor form. In 1987, Chen et al. [3] coupled the mathematical equations of shafts, and developed a simulator for the numerical simulation of the dual-medium, three-dimensional, three-phase, and multi-component thermal recovery of fractured reservoirs based on the fully implicit solution. In 1988, Sonier et al. [4] considered the effect of fractures in numerical simulation, and introduced a saturation function to simulate the dynamic changes of gravity. In 1989, Por et al. [5] developed a simulator prioritizing interactions between matrix rocks, and found that the capillary force exerted a significant effect on the recovery ratio.
Previous studies on the numerical simulation of reservoirs considering fractures mostly focused on conventional reservoirs. This is not the case with domestic studies in this field. In 1997, Wang et al. [6] simulated the repeated fracturing of low-permeability reservoirs using a single-medium model. In 2003, Jiang et al. [7] introduced the method of fracture propagation in a study on hydraulic fracturing, and discovered that this method was more effective than purely increasing fracture permeability. In 2004, Yin et al. [8] mainly explored the numerical simulation of oil recovery by imbibition, and showed that this approach was only applicable to water-wet fractured reservoir beds; moreover, the larger the capillary force in the matrix system, the higher the recovery ratio of the reservoir. In 2005, Yuan et al. [9] considered the deformation characteristics of fracture systems with changing pressure in a dual-medium numerical simulation. In 2009, Chen et al. [10] adopted a dual-medium model to simulate the production of reservoirs with fractures around injection wells, and reproduced the flow status of actual reservoir beds as far as possible by coupling multiple grid types. Based on the discrete fracture model, in 2010, Zhang et al. [11] built mathematical models for matrix systems and fracture systems containing fracture-matrix and fracture-fracture channeling terms; moreover, they solved these terms using the finite element method, obtained a generalized characteristic matrix equation, and improved discrete fracture reservoirs with controlled oil-water flow directions. Ding et al. [12] performed the sophisticated numerical simulation of late-stage reservoirs under development with spatially partitioned, temporally segmented characteristic parameters. These methods have considered the differences in fracture parameters of different regions, as well as changes of physical parameters of reservoir beds (e.g., porosity, permeability, and phase seepage) with time in reservoir development. In 2012, Song [13] created relative permeability curves considering the starting pressure gradient and the stress-sensitive effect based on physical simulation experiments; Song performed a simulation study of low-permeability reservoirs using the numerical simulation software ECLIPSE and the stress sensitivity keyword ROCKTABH, and initiated the pressure gradient keyword Threshold Pressure. Zhang et al. [14] built a two-dimensional, two-phase seepage mathematical model for fractured reservoirs that considers minimum and pseudo-starting pressure gradients, and solved them using the alternative direction implicit method. In 2013, Wang et al. [15] compiled a water-drive numerical simulation software product that simultaneously considers both starting pressure gradient and reservoir bed stress sensitivity in the numerical simulation of low-permeability reservoirs. Ji [16] explored the distribution laws of the permeability degree and fracture propagation direction by adjusting fracture system and matrix system parameters based on the dual-medium model. In 2014, Wu et al. [17] built a three-dimensional, three-phase liquid-solid coupling seepage mathematical model that considers fracture-induced anisotropy; they compiled a coupling solver applicable to ultralowpermeability reservoirs, considering the dynamic changes of reservoir permeability.
In summary, international studies rarely investigated low-permeability, compact, or other special reservoir types and their fractures. Many studies have explored low-permeability, compact, and other special fractured reservoirs, but have not sufficiently assessed the pressure-sensitive deformation characteristics induced by fracture directivity. Studying the influence of the fracture directional pressuresensitive effect on the seepage mechanism of ultralow-permeability fractured reservoir is vital. Furthermore, improving the seepage theory of ultralow-permeability fractured pressure-sensitive reservoir provides vital guiding significance for improving the development effect of reservoirs.

Mathematical Modeling of Seepage
In numerical modeling, the permeability of the reservoir bed was characterized as full-tensor permeability. The fracture permeability in the permeability tensor was calculated using a pressuresensitive equation that considers the fracture directional pressure-sensitive effect [18,19]. This equation can be used to solve the problem of waterflood development pattern simulation while considering both full anisotropy and fracture pressure-sensitive effect in different directions on oilfield sites. In addition to this primary characteristic, the following basic assumptions were adopted in the modeling: 1) A reservoir contains oil, gas, and water phases, and no consideration is given to temperature fluctuations. Oil and water are immiscible, and gas is insoluble in water. The compressibility of fluids is considered. 2) Fluid seepage obeys Darcy's law.
3) Rocks in reservoirs are an isotropic medium, subject to linear elastic deformation. 4) Pressure equilibrium is achieved instantly.

Calculation of Equivalent Permeability
In the model, rocks containing micro fractures and small-scale natural fractures were treated as an equivalent continuous medium. It was assumed that fluid pressures change simultaneously in all fractures, in which case, the changes of the effective stresses borne by matrix rocks from various sides and their volumetric deformations are symmetric. This implies that the centerline of each rod-line matrix block and the center points of each square matrix block experienced no displacement with the deformation of matrix blocks. Consequently, the region enclosed by A, B, C, and D in Fig. 1 maintained a constant volume during the deformation of matrix blocks, i.e., the constant volume boundary condition was met. For this reason, this region was adopted as the fixed volume element model for fractured mediums. Grid selection was achieved as illustrated in Fig. 1, i.e., the region enclosed by the center points A, B, C, and D of matrix blocks was assumed as one grid. At the same time, because of the presence of the stress arching effect, the entire reservoir boundary [20] could be approximated as a fixed boundary.
It was assumed that fractures develop only in the x-y plane, do not develop in the z direction, and fractures in the reservoir bed are all vertical penetrating fractures. Then, the full-tensor permeability of the fractured reservoir bed can be expressed by Eq. (1): k xx k xy k xz k yx k yy k yz k zx k zy k zz where fracture permeability is calculated from the pressure-sensitive equation, considering the fracture directional pressure-sensitive effect [18] (Eq. (2)): The components of the permeability tensor obey the following relationships: k xy = k yx ; k xz = k zx ; k yz = k zy .
The fracture permeability calculation model is compared with the previous models, its innovation lies in considering the direction of fracture pressure-sensitive effect. Could be used to calculate a single set of different direction and the direction of more groups of different fractures fracture network pressure sensitive aperture and permeability variation. The process of building and derivation of the model is to fracture as an independent object of study. Instead of conventional method used in the fracture and matrix system together as an equivalent medium to study the effect of pressure sensitive.

Mathematical Models for Seepage
Based on Section 2.1, new permeability tensor calculation model considering the effect of fracture pressure sensitive direction and the corresponding seepage flow mathematical model was established. In the model, the permeability of the reservoir bed was characterized in the form of full-tensor permeability. Single-phase and oil-gas-water three-dimensional, three-phase seepage mathematical models were built that consider fracture directional pressure-sensitive effect.

Coupling Solution Method
By summarizing the previous research results, the idea of full implicit difference iterative solution is adopted. The core idea of calculation is to calculate the permeability field of the next time step by explicit method using the pressure-sensitive equation considering the directional pressure-sensitive effect of fracture (in which the influence of stress change on fracture apertureis taken into account). An explicit method is used to calculate the permeability of each grid. Then the numerical simulation of the reservoir considering the directional pressure sensitive effect of fractures is realized by updating the permeability field. In each time step, the pressure and permeability values of the previous time step are used as the basic parameters. The calculation result of this time step is output. It is then used to calculate the permeability of the next time step. Finally, the full implicit differential iterative solution of fluid flow module is realized. Detailed solution steps are shown in the following programming process of fracture pressure-sensitive numerical simulation software module. Fig. 2 shows the comparison between the calculated results (black curve) obtained after the numerical solution of the seepage mathematical model established in this paper and the data obtained from the physical simulation test (red hollow points). It can be seen that the numerical calculation results are in good agreement with the experimental test values [18]. It shows that the mathematical model established in this paper is reliable.

Overview of the Background Reservoir
The WY (Wang Yao) reservoir of the CQ (Chang Qing) oilfield is a typical low-permeability fractured pressure-sensitive reservoir and was therefore adopted as the background reservoir. The reservoir has a mean apparent permeability of 2.29 × 10 −3 μm 2 , an apparent porosity of 13.2%, a mean original formation pressure of 9.13 MPa, a saturation pressure of 6.23 MPa, and a formation saturation pressure difference of 2.90 MPa. After waterflood development, the production energy has recovered to some extent, but the production well was subjected to water breakthrough (i.e., water injection). Because of the presence of fracture-induced anisotropy, the flow direction of injected water in the reservoir is obvious.

Basic Parameters of a Typical Element Pattern
The BLACKOIL module of the ECLIPSE software E300 simulator was adopted to build a basic seepage model for the reservoir in the Cartesian coordinate system. As the reservoir itself is not a homogeneous isotropic reservoir, in the simulated development of the reservoir, the geometric boundary of a well array would not necessarily overlap with the streamline boundary in the well array. To prevent errors caused by  Table 1 shows the mesh sensitivity analysis results (using five-point well pattern, 45°fracture, pressure sensitive effect considered as an example). It can be seen that the values of recovery and water cut after 20 years of model exploitation change with the change of grid number. When the mesh number is greater than or equal to 32,760, the recovery ratio and water cut value are relatively stable and have little change. Therefore, it is reasonable to choose 32,761 grids.

Compilation of a Fracture Pressure-Sensitive Numerical Simulation Software Module
To consider the fracture directional pressure-sensitive effect in the numerical simulation of ultralowpermeability fractured pressure-sensitive reservoir, it is necessary to apply Eq. (2) in the numerical simulation. The keyword ROCKTABH of commercial numerical simulation software ECLIPSE only yields the change laws of the overall equivalent permeability of fractured rocks with pressure. However, it cannot separately reflect the influence of the fracture pressure-sensitive effect on full-tensor permeability. Therefore, a numerical simulation software module for the pressure-sensitive full-tensor permeability of fractures with different directions was compiled to optimize both the seepage field and development method of the ultralow-permeability reservoir. As illustrated in Fig. 3, the module first reads the pressure data of the last time step, and then substitutes these data into the pressure-sensitive equation considering the fracture directional pressure-sensitive effect (Eq. (2)). Then, the permeability value of the next time step is calculated, and written into the numerical simulation file to perform operations.
The study of the model of fluid-solid coupling property for the fluid drive solid.
The number of iterations setting principles: To make sure every time step to achieve convergence, or can clearly see the main residual curve tends to level. At the same time, the time step must be greater than or equal to the minimum grid length divided by the flow velocity or rotating flow velocity. In order to ensure that each iteration is within a grid range, there will be no error in the results due to cross-grid. However, in practice, the calculation speed may increase after a period of calculation. At this time, the time step can be set larger. Finally, in order to obtain a more accurate final solution, the time step is adjusted to be smaller and a subtle calculation is performed. Comparing the convergence of the model under different time steps, it is finally determined that the time step of the fluid unit simulation is 1 month. In order to accurately compare the production effect, the comparison simulation time is set to 20 years. In actual modeling, the overlap of nodes will occur on the boundary of fluid-solid coupling due to modeling reasons. You cannot use the kneading node command to turn the overlapped nodes into one. The time step for solving the stress-sensitive permeability (solid element) is equal to the time step for solving the fluid field and fluidstructure interaction. According to the idea of fluid-solid coupling, first assume the fluid pressure distribution. Using the pressure-sensitive permeability Eq. (2) to obtain the fracture aperture and permeability, then solve the fracture flow field. Compare the current calculated pressure with the previously assumed pressure, if the difference meets a certain accuracy, then enter the next time step. Otherwise, update the pressure distribution and recalculate, and iterate to the convergence condition. When 0 < α < 1/2, the iterative method is convergent. Set the maximum number of iterations to 50.

Simulation Scheme Design and Establishment of the Production System
To fully identify the influence of the fracture directional pressure-sensitive effect on the areal pattern seepage characteristics of this ultralow-permeability reservoir, a total of 80 numerical simulation models were designed (summarized in Table 2).
By combining the actual development data of the WY reservoir with related engineering methods, the production system parameters of the above models are established as presented in the following: 1) The production well is under constant-pressure production, with a flowing bottom-hole pressure (FBHP) of 3.8 MPa.
2) The injection well is injected with water at constant pressure, with an FBHP of 24.13 MPa.
3) The term of simulated development is 20 years. 4) The simulated time step length is one month.  As shown in Fig. 3, fluid injected into element patterns always preferentially flows along the direction of the maximum principal permeability k x ; therefore, oil wells in the direction of the maximum principal permeability were affected by the injected fluid to a large extent, while those in the vertical fracture direction had low efficiency. The saturation isoline near the injection well was elliptical, and the elliptical long axis was parallel to the direction of the maximum principal permeability.
(2) Effect of fracture angle on the pressure field According to Fig. 5, the seepage resistance parallel to the direction of the maximum principal permeability was small, while that perpendicular to this direction was large. When the injectionproduction direction and the direction of the maximum principal permeability featured an included angle, after the same production time, the pressure parallel to the direction of the maximum principal permeability experienced a small loss and maintained a high level. The pressure perpendicular to this direction maintained a low level because of a large loss.

Optimization of Fracture Angle and Injection-Production Direction
Figs. 6-9 show the relationship curves between recovery percentage and water content for the same type of element patterns under the effects of fractures with different directions. Tables 3-6 provide the values of the ultimate recovery ratios and water contents after 20 years. Clearly, when the water contents of the fivespot pattern, the seven-spot pattern, the nine-spot pattern, and the cross-row well pattern increased to the same value, they had the highest recovery percentage under a fracture angle of 45°. Under fracture angles of 30°, 20°, 10°, and 0°, their recovery percentages decreased successively. For element patterns with the same fracture direction, the recovery percentage considering the pressure-sensitive effect was always lower than that without considering the pressure-sensitive effect. Note: N denotes NO (not considered); the figure immediately following N denotes the fracture angle; PS denotes pressure sensitivity; FIVE denotes the inverted five-spot pattern; SEVEN denotes the inverted seven-spot pattern; NINE denotes the inverted nine-spot pattern; and STAG denotes the cross-row well pattern. For instance, "PS0FIVE" denotes an inverted five-spot pattern model, developed with a fracture of 0°and considering fracture pressure-sensitive effect; "NPS45SEVEN" denotes an inverted seven-spot pattern model developed with a fracture of 45°and without considering fracture pressure-sensitive effect.

Optimization of Fracture Angles and Spacing Patterns
When the fracture direction remained constant, element patterns of different types and areas manifested different development effects. This is caused by differences in well pattern distance, injection-production ratio, and other parameters. The following numerical simulation models ensured the same density for various patterns and the same production system for both oil and water wells. Figs. 10a-10e show the relationship curves between recovery percentage and water content after 20 years of development. The inverted nine-spot pattern achieved the best development effect when the included angles between the direction of fracture development and the direction of the macroscopic pressure gradient were 0°, 10°, 20°, 30°, and 45°, respectively. After 20 years, the differences among the four element patterns in the ultimate water content were small; however, the inverted nine-spot pattern had a longer water-free recovery period and a higher ultimate recovery ratio.

Water content
Recovery percentage Figure 9: Relationship curve between recovery percentage and water content for the cross-row well pattern

Development Adjustment Ideas
A square inverted five-spot pattern was adopted for a case study, and the characteristics of the seepage field were analyzed under different included angles between fracture direction and well pattern direction (line direction of injection-production well. The direction from the injection well to the production well is positive). The results showed that, for the square inverted five-spot pattern, when the well pattern direction in the element pattern was neither parallel nor perpendicular to the injection-production well connecting line direction, the injection-production pattern manifested complex seepage characteristics. Analyzing the destruction and reorganization of the pattern by the permeability of the reservoir with fracture-induced anisotropy  showed that the destruction and reorganization of the injectionproduction relationship of the pattern was simultaneously affected by the intensity of anisotropy and the included angle between fracture direction and well pattern direction. Table 7 shows the correlations of the two factors with the critical failure angle of the injection-production relationship of the five-spot pattern.
Because of the presence of fracture directional pressure-sensitive effect during waterflood development, the direction of the maximum principal permeability may constantly change with changing pore fluid pressure [18]. To protect the injection-production relationship of the original pattern from destruction and improve the development effect, this study proposed to maintain a reasonable formation pressure. Furthermore, the included angle between the direction of the maximum principal permeability and well pattern direction should always be kept within a safe range to protect the original element pattern from destruction. The specific formation pressure level needs to be determined according to reservoir parameters.
The five-spot pattern developed with two fractures in the reservoir bed was taken as example in this study. The included angles between the two fractures and the direction of the macroscopic pressure gradient were 15°and 60°, respectively. The 15°fracture had an initial permeability k f1 of 11.7 × 10 −3 μm 2 , and the 60°fracture had an initial permeability k f2 of 22.9 × 10 −3 μm 2 . Other parameters are listed in the following: matrix permeability: 1.36 × 10 −3 μm 2 ; original formation pressure: 15 MPa; water injection pressure: 28 MPa; Flowing bottom-hole pressure: 9 MPa. According to the above parameters, k x /k y = 25. In this case, the equation in Table 6 can be employed to calculate the critical failure angle of the injectionproduction relationship of the element pattern (either 2.5°or 42.5°). According to the superposition principle [18], the direction of the maximum principal permeability was 46.44°at the initial state. Table 8 lists the permeability values of the two fractures after a formation pressure drop of Δp and the calculation results of the direction of the maximum principal permeability. Clearly, when formation pressure decreased by 4.3 MPa, the direction of the principal permeability rotated by 3.6°. In this case, the direction of the maximum principal permeability was 42.48°, which exceeded the critical failure angle of the injectionproduction relationship of the element pattern. Thus, to prevent the direction of the maximum principal permeability in the reservoir bed from rotating by more than 2.5°, it is necessary to control the amplitude of the formation pressure drop below 4.3 MPa. Consequently, the injection-production relationship of the original element pattern is not destroyed.
The water injected into the fractured reservoir always preferentially flows along the direction of fracture development, and it would be very difficult to recover the oil from most matrices. Moreover, the fracture directional pressure-sensitive effect would cause a rotation of the direction of the principal permeability, thus changing the flow direction of injected water [18]. Based on these characteristics, this study integrated the following development adjustment ideas: Over the course of waterflood development, the direction of the principal permeability can be adjusted by changing the formation pressure based on an understanding of the remaining oil distribution in the study area (as illustrated by the example in Table 7). This changes the flow direction of the injected water, increases the swept volume, and improves the water-drive effect. The specific range of the required formation pressure adjustment depends on reservoir parameters and the distribution of remaining oil. Table 7: Destruction and reorganization of the normal five-spot pattern by anisotropy degree and well pattern direction [46] Occurrence of destruction and reorganization

Discussion
Zhao et al. [24] revised the effective stress value of a coal mass after hydraulic fracturing by applying the principle of statistical damage mechanics; they built a model for describing the liquid-solid coupling of coal seams in hydraulic fracturing (Fig. 11). Fu [26] considered the elastic deformation of reservoir bed mediums, and established a liquid-solid coupling mathematical model of the discrete fractured-vuggy model and related numerical simulation method (Fig. 12). Hood et al. [27] performed numerical simulation of liquid-solid coupling for shales with different bedding dip angles that were collected from the Niutitang Formation in Northern Guizhou Province using RFPA 2D -Flow numerical simulation software (Fig. 13); moreover, they probed into the effects of the bedding structure on the evolution laws of compressive strength, elastic modulus, and acoustic emission signal of these shales. These scholars have conducted their studies on a small scale, without exploring the numerical simulation of single-well, macroscopic well arrays. Thus, the research findings of this study present a pragmatic advancement.  Conventional numerical simulation methods mostly equate the fracture permeability to the matrix grid. When the average pore fluid pressure changes, the fracture system and matrix system deform as the same equivalent medium. When using this kind of method, it is difficult to reflect how the pressure-sensitive permeability of fractures in different directions changes. When simulating actual reservoir production, it is difficult to reflect the change of the injected fluid flow direction caused by the directional pressure-  Figure 13: Shale rupture process and acoustic emission experiment [27] sensitive effect of fractures. As shown in Fig. 14, the range of formation pressure during the production process extends outward in a circular shape.
The finite element method is used to carry out the multi-field fluid-solid coupling calculation of the actual oil reservoir, which requires a large amount of calculation. And the fractures and matrix are treated as independent units. Although this method can separately calculate the pressure-sensitive permeability change law of the fracture, but the injected fluid will always flow preferentially along the fractures. When simulating actual reservoir production, it is also difficult to reflect the change of the injected fluid flow direction caused by the directional pressure-sensitive effect of fractures. The red color in Fig. 15 indicates the area with higher pressure, and the dark blue indicates the area with lower pressure. It can be seen that as the production process progresses, pressure always preferentially propagates along the fractures.
The numerical simulation method established in this manuscript can not only calculate the pressuresensitive deformation characteristics of fractures alone, but also reflect the change of the injected fluid flow direction caused by the directional pressure-sensitive effect of fractures in the simulation of actual reservoir production (Figs. 16 and 17). Therefore, the numerical simulation method established in this manuscript is more suitable for the simulation of actual reservoir production and development.  (1) Based on a pressure-sensitive equation that considers fracture directivity, this study built an oil-gaswater three-phase seepage mathematical model for a pressure-sensitive, full-tensor permeability, ultralow-permeability reservoir with fracture-induced anisotropy. Moreover, this model was used it to perform a numerical simulation of an ultralow-permeability fractured reservoir. The innovation of this pressure-sensitive equation lies in the fracture as an independent research object to analyze its pressure-sensitive permeability. At the same time, the directional pressuresensitive effect is also considered. (2) By utilizing this pressure-sensitive equation that considers fracture directivity, a fracture pressuresensitive numerical simulation software module was compiled to explore the seepage laws and development adjustment methods of ultralow-permeability reservoirs. As shown by the results of this study, under a fracture angle of 45°, element patterns had the highest recovery percentage. For element patterns with the same fracture direction, the recovery percentage that considers pressure-sensitive effect was always lower than the recovery percentage that was calculated without considering the pressure-sensitive effect. The inverted nine-spot pattern was more suitable for ultralow-permeability fractured pressure-sensitive reservoirs than the inverted fivespot pattern, the inverted seven-spot pattern, and the cross-row well pattern. (3) In the waterflood development of an ultralow-permeability fractured reservoir, the fracture directional pressure-sensitive effect would cause a rotation of the direction of the principal permeability, thus changing the flow direction of the injected water. With regard to this problem, this study introduced two development adjustment ideas: 1) maintaining a reasonable formation pressure, and preventing the included angle between the direction of the principal permeability and well pattern direction from reaching the critical failure angle of the symmetry element; 2) changing the formation pressure and adjusting the direction of the principal permeability, thus altering the flow direction of the injected water and improving the water-drive effect.