Abstract
Hydraulic fracturing is a key technology in unconventional reservoir production, yet many simulators only consider the single-phase flow of shale gas, ignoring the two-phase flow process caused by the retained fracturing fluid in the early stage of production. In this study, a three-dimensional fluid–gas–solid coupling reservoir model is proposed, and the governing equations which involve the early injection water phenomenon and stress-sensitive characteristics of shale gas reservoirs are established. The finite element–finite difference method was used for discretisation of stress and strain equations and the equations of flow balances. Further, a sensitivity analysis was conducted to analyse fracture deformation changes in the production. Fracture characteristics under different rock mechanics coefficients were simulated, and the influence of rock mechanics parameters on productivity was further characterised. The stimulated reservoir volume zone permeability could determine the retrofitting effect, the permeability increased from 0.02 to 0.1 mD, and cumulative gas production increased from 18.08 to 26.42 million m3, thus showing an increase of 8.34 million m3, or 46%. The effect of Young’s modulus on the yield was smaller than Poisson’s ratio and the width and length of the fractures. Production was most sensitive to the length of the fractures. The length of the fracture increased from 200 to 400 m, and the cumulative gas production increased from 26.44 to 38.34 million m3, showing an increase of 11.9 million m3, or 45%. This study deepens the understanding of the production process of shale gas reservoirs and has significance for the fluid–gas–solid coupling of shale gas reservoirs.
1. Introduction
Shale gas reservoirs have vast reserves globally and have received considerable attention for decades [1, 2]. Hydraulic fracturing is a key technology which creates a wide range of fracture networks around the well of a shale gas reservoir. The permeability of fracture networks is higher than that of the shale matrix and makes shale gas flow available; therefore, the development of the fracture networks influences the gas flow rate and total recovery. It is essential to study the characteristics of complex fracture networks after hydraulic fracturing to facilitate more accurate predictions of production.
The formation of fracture networks is complex, and many factors influence their properties [3, 4]. Laboratory experiments show that the permeability of shale is highly stress-sensitive. Hideaki et al. [5]. investigated rock permeability and found that permeability was determined by the chemomechanical effects. Liu et al. [6] introduced a two-dimensional matrix-fracture flow experiment to study the stress-sensitive phenomenon and found that the decrease in flow pressure led to the compression of the rock volume and the deformation of the fracture geometry, and the porosity and permeability of the fracture decreased directly, which resulted in a decrease in flow capacity recovery [7–9]. In a field test, Karal et al. [10] used well log data to analyse the permeability of natural fractures. However, both laboratory tests and field tests have limitations. The monitoring of the in situ stress field generally only operates in a relatively small area around the test well and is expensive, making massive implementation difficult. While in the laboratory, the test environment is typically inconsistent with the actual geologic conditions; thus, certain physical processes cannot be successfully implemented. Therefore, for convenience and to facilitate higher accuracy, numerical simulation is a suitable method to study fracture permeability with fluid–solid coupling in the entire reservoir [11, 12].
During the hydraulic fracturing process, the fracturing fluid, which is retained in the fracture network, can be produced with shale gas in the early stage of production, which decreases the fracture permeability. Many cases have been reported in field development [13]. Several researchers have established a numerical simulator to solve problems related to fluid–solid coupling from different aspects; the procedure for numerical simulation includes mathematical model establishment, model discretisation methods, and the solution of model equations. Gan and Elsworth [14] introduced a continuum model which considers the permeability of fracture rock masses in discrete fracture networks considering the coupled stress and fluid flow. Further, commercial numerical reservoir simulators use the rock compression coefficient to simulate the rock block deformation. However, the rock compression coefficient is measured under laboratory conditions; thus, it fails to represent the true variations in stress changes, especially for a complex fracture network reservoir. The unique geological characteristics of shale gas reservoirs make it difficult for the compression coefficient to meet the conditions of laboratory measurements. Julio et al. [15] noticed the geomechanical effects of the fracture rock masses incorporated in the model for fluid flow in deformable fractured media and represented the evolution of the fracture rock mass permeability. Sanaee et al. [16] analysed the core flood laboratory experimental data using two different numerical methods, and both methods involved the overburden stress effects on the fracture-matrix core permeability. It is considered more effective to solve the problem by combining a numerical reservoir simulator and an in situ stress simulator to compute the effective stress distribution in the formation.
Therefore, it is essential to build a special fluid structure coupling simulator, especially for shale gas reservoirs with complex fracture networks [17]. In the numerical simulation, there are many discrete methods to cope with the complex fracture network problem, including the finite difference, finite volume, and finite element methods. Among them, the finite element method is widely used in geomechanics. Shin and Santamarina [18] applied the finite element framework to an implicit joint-continuum model to simulate the hydromechanical process in fractured rock masses. Similarly, Li et al. [4] introduced the extended finite element method to simulate the interaction in complex fracture networks and found that different working conditions, such as fracture azimuth and fluid pressure, showed different fracture morphologies. Sutopo and Sato [19] conducted a mixed finite volume element method on a flux continuous model to simulate flow in fractured reservoirs. The finite volume method is based on irregular grid simulation. Using this, Xu et al. [20] introduced an element-based finite volume method to a three-dimensional embedded discrete fracture model and improved the flexibility of fracture patterns. Special discrete methods such as the discontinuous finite element method can solve many difficult problems, such as highly permeable anisotropic reservoirs. Vu et al. [21] introduced the boundary element method to solve the two-dimensional flow problem with an anisotropic porous medium including several intersecting curved fractures. To date, the finite difference method is the first preference for commercial simulators because of its stability and convenience. However, there are few reports on the combined finite element–finite difference (FE-FD) method to deal with the coupling problem. This is likely because the realisation of FE–FD in commercial software is often based on an explicit solution, while its numerical cost is high and stability is doubtful.
This study begins by addressing this issue. Section 2 describes the establishment of the equilibrium equation, boundary conditions, and discrete method. Sections 3 and 4 consider the influence of different fracture characteristic parameters on the reservoir properties and field performance and present conclusions.
2. Mathematical Model Equations
This chapter presents a fluid–solid coupling model. The reservoir is divided into a series of regular units on a three-dimensional scale. The physics process assumes that oil and water are immiscible and that the fluid follows Darcy’s law. The capillary force and gravity were not considered, and consideration of the matrix diffusion was excluded temporarily. For the gas component, an equation for mass conservation is established in the local element: where the physical meaning of is the well production volume per unit grid volume. The physical meaning of Equation (1) is that in the local element, the sum of the mass change caused by the diffusion term (left one) and the accumulation term (left two) is zero. In the uncoupled process, the diffusion term only considers the flow of the fluid. For the fluid–solid coupling process, the diffusion term of the fluid in the porosity media needs to consider the rock volume strain caused by solid deformation, where is the fluid density, is the porosity, is the fluid saturation, and is the fluid l source.
The rate of change of the bulk modulus is the total normal displacement of solid deformation per unit time, assuming , , and are the rock displacements along the , , and axes, respectively, and the equation is as follows:
The flow velocity of the fluid in the porous media is expressed by Darcy’s law, as described above: where is the relative permeability of the gas, K is a function of the gas saturation, is the gas viscosity, and is the gas pressure. Because the capillary force and gravity are not considered, the gas pressure is equal to the fluid pressure .
The second term is the derivative of the cumulative term, over time. The derivative of the cumulative term over time is divided into three parts: density, porosity, and saturation, as follows:
For the unknown terms in the equation, the volume strain, pressure, and saturation are the first unknown terms, which can be solved directly. The remaining unknown terms, including permeability, phase permeability, density, and porosity, are the second unknown terms. These are functions of time and the first unknown terms. They must be transformed into a function of the first unknown terms. Next, we discuss how to represent the second unknown terms.
For the gas component, the change in density with pressure cannot be directly expressed by the compression coefficient. The compression coefficient is not constant or linear with pressure, showing nonlinear characteristics. It is assumed that shale gas is an ideal gas and that the density can be expressed by the equation of state of the ideal gas:
By deriving the above formula from time, we can obtain the reciprocal of density with respect to time, which is proportional to the derivative of pressure:
Permeability is important for accurately describing the fluid flow. At present, many empirical formulas have been used to simulate the change in permeability with effective stress. However, for shale gas reservoirs with complex fracture networks, the change in the stress field is more complex because of the characteristics of the fracture networks. Thus, laboratory measurement conditions cannot accurately match the underground conditions; therefore, it is difficult to directly apply the permeability formula related to the effective stress. However, another definition of permeability may also be used, in which permeability is a function of porosity. Under the condition of low permeability, the permeability is expressed as
Porosity is expressed as the fluid pressure , block modulus of elasticity, initial porosity , the Biot coefficient , and volume strain (Erdinc, 2017):
By integrating these formulas, the control equations for the displacement, pressure, and saturation are established:
For the control equation of the water component, the derivation method is the same as the gas component with one difference, i.e., the compressibility of water is much smaller than that of gas, so the compressibility of water capacity can be expressed by a constant compression coefficient, and the water density is directly proportional to the compression coefficient: where is the density of the liquid, and is the compressibility of the liquid.
To obtain the solution of the model, the aforementioned governing equation requires auxiliary conditions. In addition to the assumption that the capillary force and gravity are ignored, it is also assumed that the pressure of the gas and that of water can be uniformly expressed as fluid pressure:
According to the principle of phase equilibrium, the sum of the liquid phase saturation and gas phase saturation is equal to 1, that is,
The governing equation of water can be expressed as follows:
Next, the in situ stress governing equation for a local element is described. For this, the reservoir is assumed to be an elastic medium; therefore, in a given element of the reservoir, the sum of the derivations of the stress in all directions is zero: where is the total stress tensor, and its physical meaning is the sum of the effective stress and fluid pressure in the element. Generally, the Biot coefficient is added to describe the degree of fluid pressure on rock deformation. The effective stress of the element can be expressed as where the plus/minus sign indicates the direction of stress.
In the numerical model of the reservoir, the absolute displacement of the element has no practical application significance. Generally, the relative displacement is considered to estimate the degree of deformation of an element. It is assumed that under the initial conditions of the reservoir, a force balance is obtained between the initial in situ stress and fluid pressure in the shale gas reservoir, and the effective stress is zero. Thus, the displacement of each point in the initial state of the formation should also be zero. When the fluid pressure changes, the calculated value corresponds to the relative displacement of the element. By substituting these conditions into the above formula, the following relationships can be obtained:
Because the unknown term to be solved is displacement, these equations of stress need to be transformed into a displacement equation. After special conservation, the following equilibrium equations of displacement and fluid pressure are obtained, where denotes the elastic tensor, and is the displacement tensor:
Finally, Equations (10), (14), and (18) are combined to form a complete fluid–solid coupling equation for the shale gas reservoir. The above mathematical equations need to be discretised into numerical models, in which the standard finite element method is used to discretise the geomechanical equations and the finite difference method to discretise the fluid equations. In the finite difference method, the upstream weight method was used to determine the relative permeability to ensure the stability of the numerical simulation. The displacement variable is at the corner of the element, and the pressure and saturation are at the centre of the element. The locations of the pressure, saturation, and displacement in the cell are shown in Figure 1.
The discretisation results are as follows, where is the trial function of the finite element, is the unit volume, and is the reservoir area:
In the production of shale gas reservoirs, the conditions of production are usually well-maintained under constant wellbore pressure conditions. The Peaceman well model was used to simulate the production well in this study.
The initial condition of the reservoir was set as follows. Under the assumption of setting the fluid pressure, the initial in situ stress in the model is obtained by the following formula:
The boundary conditions for the stress and flow are listed. The boundary of the flow conforms to Neumann boundary conditions, and the shale gas reservoir is in a closed state:
For the boundary condition of the in situ stress, the normal direction of the upper boundary is influenced by the constant overburden stress. The displacement from other surfaces in the direction of the normal direction was 0. The mathematical expressions are as follows:
The complex fracture network consists of a hydraulic fracture and a natural fracture group around it. To describe the fluid flow between the fracture and rock around smoothly, the local grid refinement (LGR) scheme was used to indicate a high-permeability fracture network. The distribution of the cells around the hydraulic fractures obeys a logarithmic distribution, and the natural fractures caused by hydraulic fracturing around the hydraulic fracture were represented in the REV region. The coupled set of equations was then solved using the Newton–Raphson method.
3. Results
A numerical simulation model for the horizontal well staged fracturing development of a shale gas reservoir was established (Figure 2). The basic geological parameters of the base case are listed in Table 1. The model was 1220 m long, 400 m wide, and 25 m high, with an area of 0.488 km2. The grid cells of length direction, width direction, and height direction were , respectively, and the total number of grids was 420,000.
The model considers the fluid–structure coupling effect, and the permeability and porosity change with the change in pressure in the production process. The Durcan model was used to describe the stress-sensitive effect of the shale reservoir. The model considered the adsorption of shale gas. Curtis, a well-known scholar in the field of shale gas research in the United States, proposed in 2002 that shale gas is essentially continuously generates biochemically derived gas and thermally derived gas. He states that shale gas has a hidden accumulation mechanism with multiple types of lithological closure, a relatively short migration distance, and a general formation saturation gas content, which can exist in both a free state in natural cracks and pores and in an adsorbed state on the surface of kerogen and clay particles, and can also even dissolve into kerogen and asphaltene states [22]. In the study of adsorption phenomena, the adsorption isotherm curve is the most common way to express the adsorption performance, and its shape can well reflect the physical and chemical interactions between the adsorbent and the adsorbed substance. Therefore, in this study, the Langmuir isothermal adsorption equation was used to describe the change process of adsorbed gas with pressure in the shale gas production process. The numerical simulation was conducted using the simulation software of Schlumberger.
3.1. Degree of Rock Deformation
The degree of formation degeneration is controlled by Young’s modulus and Poisson’s ratio . In an elastomer, the higher the value of Young’s modulus, the less obvious the deformation effect is under the premise of the same stress [23, 24]. Poisson’s ratio determines the ratio of the longitudinal to transverse deformation of the elastomer. In the fluid–structure coupling process of this study, Young’s modulus determines the volume modulus of the rock and controls the volume strain of the reservoir in the formation [25, 26].
Young’s modulus and Poisson’s ratio mainly affect the fracture elongation effect, which then affects the gas well productivity. Therefore, before studying the influence of Young’s modulus and Poisson’s ratio on productivity, it is necessary to simulate the fracture generation process. The fracturing simulation software FRCPT was used to simulate the fracturing process, and the effect of fracture extension under different Young’s moduli and Poisson’s ratios under the same fracturing scale and pumping program were studied. The research conclusions are presented in Tables 2 and 3.
It is evident from the tables that the influence of Young’s modulus on the fracture elongation effect is greater than that of the Poisson’s ratio, mainly because in general, Poisson’s ratio has a smaller variation range, while Young’s modulus has a larger variation range. Under different Poisson’s ratios, the difference in fracture width variation was no more than 0.3 m, and the difference in fracture length was no more than 1.4 m. Under such width and fracture length variations, there was almost no influence on productivity. Therefore, this sensitivity analysis does not discuss the influence of Poisson’s ratio on production and instead focuses on the influence of Young’s modulus on productivity.
Based on the fracturing simulation with different Young’s moduli, a reservoir numerical simulation model was established to predict the productivity of horizontal wells with different Young’s moduli (Table 4).
The final results (Figures 3 and 4) show that with the increase in Young’s modulus, the range of SRV zones becomes smaller under the same fracturing system, and accordingly, the reservoir transformation effect becomes worse. The difference in peak gas production was not large, increasing from 34,000 to 34,500 m3/d. The cumulative gas production decreased from 25.99 to 23.27 million m3, a production decrease of 10%.
(a)
(b)
(c)
(d)
(e)
(f)
Compared with the other three parameters, Young’s modulus has less effect on the production, mainly in the process of 2 GPa to 10 GPa, the fracturing parameters, fracturing scale, and fracturing system under the same fracture of the SRV area changed little, and the extent and effect of the SRV area were the key factors influencing the shale gas reservoir; therefore, in comparison, the influence of the output of Young’s modulus on the other three parameters.
3.2. SRV Regional Permeability
Owing to the influence of hydraulic fracturing technology, a large number of natural microfractures are produced around the pressure fractures. These natural fractures improve the reservoir permeability in these areas and are the main seepage channels of shale gas reservoirs, which are important factors for improving productivity. We call these artificially fractured high-permeability SRV zones. Permeability in SRV zones is usually obtained through subsurface core sampling and logging experiments.
Under the original conditions, the permeability of the shale gas reservoir is relatively small, and commercial airflow cannot be obtained without reservoir modification. The permeability of the reformed shale gas reservoir can reach 0.1–0.01 mD in the modified area. Combined with the high-pressure characteristics of shale gas reservoirs, commercial airflow can be obtained under the conditions of permeability.
On the one hand, the permeability of the SRV zone determines the productivity of the shale gas well. On the other hand, the SRV zone itself is small; thus, the productivity will inevitably reach an inflexion point after the permeability reaches a certain height. Therefore, based on the permeability of the revamp zone, five groups of parameters of the revamp zone (0.02, 0.04, 0.06, 0.08, and 0.1 mD) were selected to study the influence of the permeability of the SRV zone on the productivity of a single well.
The results (Figures 5 and 6 and Table 5) show that the SRV zone permeability had a significant impact on the final cumulative gas production per well. When the permeability increased from 0.02 to 0.1 mD, the cumulative gas production increased from 18.08 to 26.42 million m3, a 46% increase in production. However, when the SRV permeability reached 0.06 mD, the increase in cumulative gas production decreased. This is owing to two factors: (1) the SRV area permeability could meet the corresponding area of the gas extraction; thus, permeability was no longer the bottleneck of production and (2) from the point of production at the end of pressure distribution, the gas was output from most but not all of the SRV area of low permeability, and the pressure nearly did not decrease the SRV area, which reached a considerable recovery. At that time, the increasing permeability had not increased the gas well production and only by expanding the SRV area could the single-well production capacity be increased; therefore, in the SRV transformation, the transformed effect of permeability did not affect the SRV area, but the SRV area itself expanded in width, as demonstrated in the fracture length sensitivity analysis.
(a)
(b)
(c)
(d)
(e)
(f)
3.3. Fracture Spacing
Fracture spacing is an important factor that determines the effect of reservoir development [27, 28]. Moreover, fracture spacing can change the control area of fractures and the variation law of stress field between fractures.
In certain cases, the horizontal well length determines the spacing between the fractures’ own control areas, in that, the greater the gap spacing and the fracture control area, the less the recovery efficiency of the horizontal control area, cost reduction, the 1000 m horizontal wells, based on base case, respectively, set up seam article 6, article 7, article 8 seam sewing seam, article 9, 10, the corresponding joint spacing for, in order to study the effect of joint spacing on horizontal well productivity.
The results (Figures 7 and 8 and Table 6) showed that with the decrease in fracture spacing, the number of fractures increased and the SRV zone expanded and the cumulative gas production per well increased accordingly. Fractures in the bar or the increase in the number did not quickly improve the recovery of a single well. From the point of this model, when the fracture number increased to 7, the pressure between the fractures began to decline, compared with 8 fractures, article 9 and article 10 of fractures of the pressure did not expand much or affect the scope. The inflection point on this performance in the accumulative gas production appeared between 6 to 7 fractures; unlike the increase in fracture length, the increase in fracture number was less sensitive to enhanced recovery. The cumulative gas production of 6 to 7 fractures increased from 4 m3 to 4 m3, and one fracture increased production by 13 to 2.46 million m3. However, with 8 fractures, 1 fracture increased production by 9% to 1.69 million m3, and the production capacity of a single fracture decreased.
(a)
(b)
(c)
(d)
(e)
(f)
Moreover, in the process of fracturing, the increase in fracture stages is associated with an increase in the economic cost. Therefore, in the optimisation process of fracture stages, the influence of economic factors should be fully considered, and the optimisation of fracture stages should be conducted with economic indicators as the object.
3.4. Fracture Length
The following will consider the impact of fracture length on development results. The fracture length will directly determine the width of the SRV zone, which has a significant impact on seepage flow [29, 30].
Based on base case, three groups of fracture lengths were set: 200 m (fracture half-length is 100 m), 300 m (half-length, 150 m), and 400 m (half-length, 200 m). In all three cases, the length of the SRV zone increased with the increase in fracture length, whereas the fracture width remained unchanged.
The numerical simulation results (Figures 9 and 10 and Table 7) show that the fracture length has a direct effect on the cumulative gas production of a single well. The longer the fracture, the greater the cumulative gas production of a single well, and there was no inflexion point. Therefore, the size of the SRV region is a key factor in determining the success of an SRV. When the fracture length increased from 200 to 400 m, cumulative gas production per well increased from 4 m3 to 4 m3, an increase of 45%.
(a)
(b)
(c)
(d)
4. Conclusions
(1)The effect of Young’s modulus on yield is smaller than that of the other three parameters, and under different Poisson’s ratios, the width and length of fractures vary by a relatively small range. Poisson’s ratio itself changed slightly from 0.2 to 0.32, the fracture width and length changed by 0.3 and 1.4 m, respectively, and the variation in the width and length had almost no effect on productivity(2)The fracture length has a direct effect on cumulative gas production per well, in that, the longer the fracture, the greater the cumulative gas production per well, and there was no inflexion point. The length of the fracture increased from 200 to 400 m, and the cumulative gas production increased from 26.44 to 38.34 million m3, an increase of 11.9 million m3, or 45%. Therefore, in SRV retrofitting, it is not the SRV zone permeability that determines the retrofitting effect, but the SRV zone expansion width(3)As the number of fractures increased from 6 to 10, the cumulative gas production increased from 19.26 to 26.44 million m3, an increase of 7.18 million m3, or 37%. However, in the process of fracturing, such an increase in fracture stages would mean an increase in the economic cost. Therefore, in the optimisation process of fracture stages and the influence of economic factors should be fully considered, and the optimisation of fracture stages should be conducted with economic indicators as the object
Data Availability
The research results of this paper are based on the application of formula derivation, and all the formulas and related data are included in the paper.
Conflicts of Interest
The authors declare that they have no conflicts of interest.