Shale gas reservoir multi-section fractured horizontal well prediction method considering stress sensitivity

Abstract

The invention relates to the technical field of shale gas reservoir development, in particular to a shale gas reservoir multi-section fractured horizontal well prediction method considering stress sensitivity. The method is characterized by comprising the following steps: acquiring shale gas reservoir related data, including physical property parameters of a reservoir and shale gas and production parameters of a gas well; calculating the relationship among the reservoir permeability, the shale gas physical property and the reservoir pressure, and performing fitting; dividing grids according to crack types, and determining crack positions through certainty and a random modeling method; calculating parameters such as crack section length, equivalent distance and conductivity; constructing a relationship among shale matrix grid pressure, fracture grid pressure and conductivity, and forming a numerical model; the physical property parameters are substituted, and the gas well yield and the cumulative yield under different working systems are solved through a full implicit method; and predicting a yield decline rule and an accumulated yield change rule of the shale gas well according to the curve graph. The productivity prediction method which is accurate, rapid and wide in application range is provided for efficient development of shale gas reservoirs.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of shale gas reservoir development, in particular to a shale gas reservoir multi-stage fracturing horizontal well prediction method considering stress sensitivity. BACKGROUND

[0002] In recent years, China has discovered shale gas reservoirs with great commercial development potential in Sichuan, Shaanxi and Yunnan, among which Fuling shale gas reservoir is a typical representative. The successful development of shale gas fields is of great significance for promoting China's energy structure adjustment, alleviating the supply pressure of natural gas market in the central and eastern regions, and ensuring national energy security.

[0003] On the one hand, shale gas reservoirs usually have poor porosity and permeability, and usually need to use hydraulic fracturing to improve shale gas production capacity. At the same time, the natural fractures of shale gas reservoirs that have been put into commercial development are generally developed, resulting in complex fracture distribution characteristics of shale gas reservoirs. On the other hand, a large number of laboratory experiments and field production example analyses show that shale matrix usually has strong stress sensitivity effect, causing the permeability of shale reservoir to change with the change of pore pressure; at the same time, the decrease of shale pore pressure will cause the desorption of adsorbed gas, and the above two reactions will cause the control equation of shale matrix to have strong nonlinear characteristics. At present, the nonlinear characteristics caused by stress sensitivity are usually handled by using small perturbation method to obtain its approximate solution, and the approximate solution obtained by using this method usually has certain error, and the production solution or pressure solution obtained is usually a semi-analytical solution or an approximate solution, which makes its application range narrow. Secondly, the adsorption / desorption effect of shale gas will also cause the control equation to have strong nonlinearity, and the existing technology usually uses the definition of dimensionless adsorption / desorption constant to represent the adsorption / desorption effect of shale gas, so as to reduce the nonlinearity of the control equation. The semi-analytical solution obtained by the above method cannot accurately and intuitively reflect the influence of shale gas adsorption / desorption effect on production decline and cumulative production due to the existence of quadratic term.

[0004] Furthermore, the coupling of natural fractures and artificial fractures will produce complex fracture networks, and the existing technology usually uses fracture wall flow normalization and pressure normalization method when dealing with the coupling of complex fracture networks and shale matrix, but for complex fracture networks, the fracture distribution characteristics are not consistent, which leads to the inability to accurately simulate the pressure loss and flow loss near the fracture end. Finally, in the actual production process of shale gas wells, mixed working systems such as pressure control production or pressure release production are usually used, and the existing technology usually cannot have production decline and cumulative production under different working systems of shale gas wells.

[0005] In general, the prior art has certain limitations in simulating the production decline and cumulative production of multi-stage fractured horizontal wells in shale gas reservoirs. At present, there is a lack of an accurate, fast and widely applicable multi-stage fractured horizontal well productivity prediction method that takes into account the stress sensitivity of shale gas reservoirs. SUMMARY

[0006] The present application aims to overcome the deficiencies of the prior art and provides a stress-sensitive shale gas reservoir multi-stage fractured horizontal well prediction method, comprising the following steps:

[0007] Obtain the physical property parameter data of the shale gas reservoir and shale gas and the working system data of the gas well being studied;

[0008] Based on the physical property data of the shale gas reservoir, calculate the relationship and graph between the shale gas reservoir permeability, shale gas physical property parameters and shale gas reservoir pressure, and use the polynomial fitting method to fit the relationship between the shale gas physical property parameters and the shale gas reservoir pressure;

[0009] Divide the grid according to the fracture type in the shale reservoir, and use the deterministic modeling method and the stochastic modeling method to determine the position of the fracture in the shale matrix;

[0010] According to the distribution of the matrix grid and the position of the fracture in the matrix grid, respectively calculate the embedded fracture segment length, equivalent distance and conductivity parameter in the shale matrix grid;

[0011] Respectively construct the corresponding relationship between the shale matrix grid pressure, fracture grid pressure and conductivity to form a numerical model of the stress-sensitive shale gas reservoir multi-stage fractured horizontal well;

[0012] Substitute the obtained physical property parameters into the stress-sensitive shale gas reservoir multi-stage fractured horizontal well numerical model, and use the fully implicit method to respectively solve the gas well production decline and cumulative production under the two working systems of pressure control production and pressure release production, and draw the relationship curve between the gas well production decline and cumulative production and the production time;

[0013] According to the stress-sensitive shale gas reservoir multi-stage fractured horizontal well production decline and cumulative production-time relationship curve drawn, predict the production decline law and cumulative production change law of the shale gas well being studied.

[0014] The relationship between the shale gas reservoir permeability and the reservoir pressure is calculated by the following formula:

[0015]

[0016] In the formula: K mi is the initial permeability of the shale gas reservoir, mD; γ is the stress sensitivity, MPa -1 ; P iP0 is the original formation pressure, MPa; p is the shale gas reservoir pressure, MPa.

[0017] The shale gas properties include shale gas compressibility factor, shale gas compressibility coefficient, shale gas volume coefficient, shale gas viscosity, shale gas adsorption / desorption compressibility coefficient;

[0018] The shale gas compressibility factor is obtained by solving the following two Z g expressions by iteration method, and then substituting the obtained y value into any one Z g equation for calculation:

[0019]

[0020]

[0021] In the formula: Z is the apparent corresponding pressure, dimensionless; Z g is the shale gas compressibility factor, dimensionless; t is the apparent corresponding temperature reciprocal, dimensionless; y is the shale gas corresponding density, dimensionless;

[0022] The relationship between shale gas compressibility coefficient and reservoir pressure is calculated by the following formula:

[0023]

[0024] In the formula: Z is the apparent corresponding pressure, dimensionless; Z g is the shale gas compressibility factor, dimensionless; The normalized temperature is the relative temperature, which is the ratio of the current temperature to a certain reference temperature, dimensionless;

[0025] The relationship between shale gas volume coefficient and reservoir pressure is calculated by the following formula:

[0026]

[0027] In the formula: p sc and T sc are the pressure and temperature under standard conditions, MPa, K; p and T are the shale gas reservoir pressure and temperature, MPa, K; Z g is the shale gas compressibility factor, dimensionless;

[0028] The relationship between shale gas viscosity and reservoir pressure is calculated by using the modified Dempsey model:

[0029]

[0030] In the formula: μ g is the shale gas viscosity, Pa·s;

[0031] μ1 = (1.709 x 10 -5 -2.062 x 10 -6 γ g )(1.8T + 32) + 8.118 x 10 -3 -6.15 x 10 -3 log(γ g )

[0032] a0 = -2.4621182 a1 = 2.970547414 a2 = -0.286264054 a3 = 0.008054205

[0033] a4 = 2.80860949 a5 = -3.49803305 a6 = 0.36037302 a7 = -0.01044324

[0034] a8 = -0.793385648 a9 = 1.39643306 a 10 = -0.149144925 a 11 = 0.004410155

[0035] a 12 = 0.083938718 a 13 = -0.186408848 a 14 = 0.020336788 a 15 = -0.000609579

[0036] The relationship between the shale gas adsorption / desorption compressibility and the reservoir pressure is calculated by the following formula:

[0037]

[0038] wherein C d is the adsorption / desorption compressibility, dimensionless; p s is the shale density, kg / m 3 ; V L is the liquid volume interacting with the shale gas, m 3 ; p L is the pressure of the liquid, Pa; p is the current shale gas pressure, Pa.

[0039] The deterministic modeling method determines the location of artificial fractures in the shale matrix, and the stochastic modeling method determines the location of natural fractures in the shale matrix, and the natural fracture stochastic modeling method is as follows:

[0040] Natural fracture starting point:

[0041] X nf1= D xnf + (x e - 2 * D xnf ) * rand (1); Y hf1 = D ynf + (y e - 2 * D ynf ) * rand (1)

[0042] Natural fracture end point:

[0043] X nf2 = X nf1 + sign (rand (1) - 0.5) * d xnf ; Y hf2 = Y hf1 + sign (rand (1) - 0.5) * d ynf

[0044] Where:

[0045]

[0046] In the formula, rand (x) is a random function that generates a random seed, and sign (x) is a sign function.

[0047] The fracture segment length embedded in the shale matrix grid is calculated by the following formula:

[0048]

[0049] In the formula, (x s1 , y s1 ) and (x s2 , y s2 ) are the coordinates of the two intersection points of the fracture and the matrix grid line, respectively.

[0050] The equivalent distance is divided into three categories: a. The equivalent distance from the matrix grid to the fracture; b. The equivalent distance between adjacent numbered fracture segments in the same fracture; c. The equivalent distance of intersecting fracture segments;

[0051] a. The equivalent distance from the matrix grid to the fracture is calculated using the following formula:

[0052]

[0053] In the formula, S is the cross-sectional area of the matrix grid into which the fracture segment is embedded, x n is the normal distance from the matrix grid to the fracture;

[0054] b. The equivalent distance between adjacent numbered fracture segments in the same fracture is calculated using the following formula:

[0055]

[0056] In the formula,

[0057] This represents the equivalent distance between the first segment of a crack with adjacent numbers in the same crack and the length of its embedment into the matrix mesh. The equivalent distance representing the length of another adjacent crack segment embedded into the matrix mesh within the same crack;

[0058] l f1 The length of the first adjacent numbered segment of a crack within the same crack embedded into the matrix mesh;

[0059] l f2 The length of another adjacent-numbered segment of the same crack embedded into the matrix grid;

[0060] l f1i After discretizing the first crack segment embedded in the matrix mesh at equal intervals, the distance from the midpoint of the i-th discrete segment to the intersection of adjacent crack segments of the same crack segment (calculated using the two-point distance formula) is N. f1 part;

[0061] l f2i After discretizing another crack segment embedded in the matrix mesh at equal intervals, the distance from the midpoint of the i-th discrete segment to the intersection of adjacent crack segments of the same crack segment (calculated using the two-point distance formula) is N. f2 part;

[0062] The number of discrete segments is based on the shorter of the two crack segments. The shorter crack has 100 discrete segments, and the longer crack has the following number of discrete segments:

[0063]

[0064] c. The equivalent distance between intersecting crack segments is calculated using the following formula:

[0065]

[0066] In the formula,

[0067] It is the equivalent distance from the i-th discrete segment in a certain intersecting crack segment to the intersection point;

[0068] It is the equivalent distance from the i-th small discrete segment in another intersecting crack segment to the intersection point;

[0069] N f1 d is the number of discrete segments, where the length of the first segment of the intersecting crack is embedded in the matrix, at equal intervals. f1i Let be the equivalent distance from the i-th discrete segment in a given intersecting crack segment to the intersection point;

[0070] N f2 d is the number of discrete segments, where the length of the other segment of the intersecting crack is embedded in the matrix, at equal intervals. f2i This is the equivalent distance from the i-th small discrete segment in another intersecting crack segment to the intersection point;

[0071] The number of discrete segments is based on the shorter crack among the intersecting cracks. The shorter crack has 200 discrete segments, and the longer crack has the following number of discrete segments:

[0072]

[0073] The conductivity is classified into four categories: a. conductivity between matrix grids; b. conductivity between matrix grids and fracture grids; c. conductivity of adjacent fracture segments within the same fracture; d. conductivity of intersecting fracture segments.

[0074] a. Matrix mesh and its conductivity:

[0075] Grid in xx direction:

[0076]

[0077] YY direction grid:

[0078]

[0079] In the formula, dx, dy, and dz are the mesh sizes in the x, y, and z directions, respectively;

[0080] b. Conductivity of matrix mesh and fracture mesh:

[0081]

[0082] In the above formula, d m-f k is the length of the crack embedded into the matrix mesh. m and k f These are the permeability of the matrix and the cracks, respectively.

[0083] c. Conductivity of adjacent crack segments with the same crack number:

[0084]

[0085] d. Conductivity of intersecting fracture segments:

[0086]

[0087] In the formula k f1 and k f2 These represent the permeability of the intersecting fractures; w f1 and w f2 These represent the widths of the intersecting cracks.

[0088] The relationship between the pressure and conductivity of the shale matrix grid is calculated using the following equation:

[0089]

[0090] The relationship between crack mesh pressure and conductivity is calculated using the following equation:

[0091]

[0092] By combining the above equations relating shale matrix grid pressure and conductivity, fracture grid pressure and conductivity, and the Peaceman well equation, a numerical model of a multi-stage fractured horizontal well in a stress-sensitive shale gas reservoir can be obtained.

[0093] The methods for the two working systems of controlled pressure production and depressurized production are calculated using the following formula:

[0094] The production conditions for depressurization production are: p w =p wf

[0095] The production conditions for pressure-controlled production are: p w =p wmax -D t *n

[0096] In the two formulas above, p wf The bottom-hole flowing pressure during depressurization production; p wmax To control the maximum flowing pressure at the bottom of the well during pressure production; D t The pressure level is denoted by n; n is the nth time step.

[0097] The convergence condition for the fully implicit method is as follows:

[0098] Matrix equation:

[0099]

[0100] Crack equation:

[0101]

[0102] Compared with the prior art, the present invention has the following advantages:

[0103] This invention utilizes numerical simulation methods and embedded discrete fracture models to establish pressure-release and pressure-controlled production models for multi-stage fractured horizontal wells in stress-sensitive shale gas reservoirs. It analyzes the impact of pressure-release and pressure-controlled production on the production decline and cumulative production of multi-stage fractured horizontal wells using the controlled variable method, establishing the production decline law of shale gas reservoirs under different operating conditions. By adopting pressure-release and pressure-controlled production as common development modes for shale gas fractured horizontal wells, it overcomes the limitation that the rapid production decline of shale gas wells makes it impossible to directly adopt constant production, thus conforming to the actual development mode of multi-stage fractured horizontal wells in shale gas reservoirs. By comprehensively considering the production decline and cumulative production characteristics of multi-stage fractured horizontal wells in shale gas reservoirs under pressure-release and pressure-controlled production, it establishes the production decline and cumulative production variation laws of multi-stage fractured horizontal wells in stress-sensitive shale gas reservoirs under different development modes. This approach comprehensively considers the characteristics of shale gas reservoirs and development modes, engineering practice and theoretical application, overcoming the shortcomings of existing technical solutions that can only analyze production decline and cumulative production variation laws under constant production and constant pressure development modes. This invention improves the accuracy of production capacity prediction for multi-stage fractured horizontal wells in stress-sensitive shale gas reservoirs under different development modes, and can meet the shale gas production needs of the study area.

[0104] This invention, based on oil and gas reservoir numerical simulation methods and embedded discrete fracture models, establishes the production decline and cumulative production variation laws of multi-stage fractured horizontal wells in stress-sensitive shale gas reservoirs, enabling the prediction of production capacity variation laws and optimization of hydraulic fracture parameters for this type of gas well. The matrix grid is partitioned using oil and gas reservoir numerical simulation methods. The fracture distribution characteristics of shale gas reservoirs in the study area are selected, and an embedded discrete fracture model is constructed using embedded discrete fractures. Simultaneously, a fully implicit method is used to solve and predict the production of newly commissioned multi-stage fractured horizontal wells in the study area's shale gas reservoirs. This establishes a method for optimizing hydraulic fracturing parameters for this type of gas reservoir. This approach avoids the overly idealistic problem of the Blasingame chart (BLASINGAME TA, MCCRAY TL, LEE W J. Decline curve analysis for variable pressure drop / variable flowrate systems [C]. SPE 21513, 1991.) established by the semi-analytical method in existing technologies, and also solves the shortcomings of poor convergence and narrow applicability of production capacity prediction models established by numerical or analytical methods in existing technologies. The production capacity predicted by the multi-stage fracturing horizontal wells in stress-sensitive shale gas reservoirs using the different operating systems of pressure control production and pressure release production proposed in this invention has a high degree of agreement with the actual production data of shale gas reservoirs. In the example, the production decline data of one multi-stage fracturing horizontal well is compared with the prediction results of this invention, and the relative error is less than 6.1%, which shows high prediction accuracy and can basically meet the actual needs of shale gas reservoirs in the field. Attached Figure Description

[0105] Figure 1 A flowchart illustrating the method for predicting multi-stage fracturing horizontal wells in stress-sensitive shale gas reservoirs, as described in this invention.

[0106] Figure 2 A flowchart illustrating an embodiment of the method for predicting multi-stage fracturing horizontal wells in stress-sensitive shale gas reservoirs according to the present invention;

[0107] Figure 3 This is a distribution map of crack morphology in the matrix mesh;

[0108] Figure 4 A graph showing the relationship between shale gas compressibility factor and pressure;

[0109] Figure 5 This is a graph showing the relationship between the compressibility coefficient of shale gas and pressure.

[0110] Figure 6 A graph showing the results of diminishing returns and cumulative returns under different development models;

[0111] Figure 7 This is a graph showing the impact of depressurization level on output decline and cumulative output under the depressurization production mode.

[0112] Figure 8 This is a graph showing the impact of pressure control level on output decline and cumulative output under pressure-controlled production mode.

[0113] Figure 9 The graph shows the effects of stress sensitivity on production decline and cumulative production.

[0114] Figure 10 The graph shows the impact of crack distribution characteristics on production decline and cumulative production.

[0115] Figure 11 This is a comparison chart of the production decline results of well SGFHW01 in the study block with the actual production data. Detailed Implementation

[0116] The specific embodiments of the present invention will be described in detail below. The present invention can be implemented in many other ways different from those described herein, and those skilled in the art can make similar modifications without departing from the spirit of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed.

[0117] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. The terminology used is for describing particular embodiments only and is not intended to limit the invention.

[0118] A method for predicting the productivity of multi-stage fractured horizontal wells in stress-sensitive shale gas reservoirs includes the following steps:

[0119] S101: Obtain physical property data of shale gas reservoir and shale gas, as well as operating system data of the gas wells under study.

[0120] Shale gas reservoir physical properties: density ρ of shale s Original formation pressure P i Initial permeability K of shale matrix mi Shale porosity φ m The porosity compressibility coefficient c of shale m , wellbore radius r w Gas reservoir thickness h, gas reservoir extent (x) e ,y e The stress sensitivity coefficient γ of the gas reservoir, the temperature T of the shale gas reservoir, and the permeability k of the fracture system (natural and artificial fractures) are all considered. f Porosity φ fm Crack width w f and the compressibility coefficient c of the crack f ;

[0121] Physical properties of shale gas: Langmuir adsorption pressure P L Langmuir adsorption volume V L The compressibility factor Z of shale gas g shale gas volume factor B g The compressibility coefficient C of shale gas g Viscosity μ of shale gas g Shale gas adsorption / desorption compressibility coefficient C d ;

[0122] Work system data: Gas well work system p w It is divided into a pressure release production system and a pressure control production system.

[0123] S102: Based on the physical property data of shale gas reservoirs, calculate the relationship and relationship chart between shale gas reservoir permeability, shale gas physical property parameters and shale gas reservoir pressure, and use the polynomial fitting method to fit the relationship between shale gas physical property parameters and shale gas reservoir pressure.

[0124] Based on the reservoir physical property data of shale gas reservoirs, the relationship between reservoir permeability and reservoir pressure p is calculated using the following formula:

[0125]

[0126] In the formula: K mi γ represents the initial permeability of the shale gas reservoir, in mD; γ is the stress-sensitive coefficient, in MPa. -1 ;P i p represents the original formation pressure, in MPa; p represents the shale gas reservoir pressure, in MPa.

[0127] Based on the physical properties of shale gas, the relationship between shale gas compressibility factor, shale gas compressibility coefficient, shale gas volume coefficient, shale gas viscosity, shale gas adsorption / desorption compressibility coefficient, and reservoir pressure is calculated.

[0128] Shale gas compressibility factor is obtained by combining the following two Z... g The expression is solved iteratively to obtain the y-value, and then the obtained y-value is substituted into any Z-value. g Calculate using the equation:

[0129]

[0130] In the formula: The corresponding pressure is dimensionless; Z g is the shale gas compressibility factor, dimensionless; t is the reciprocal of the apparent corresponding temperature, dimensionless; y is the corresponding density of shale gas, dimensionless.

[0131] The compressibility coefficient of shale gas is calculated using the following formula:

[0132]

[0133] In the formula: The corresponding pressure is dimensionless; Z g is the shale gas compressibility factor, dimensionless; This represents the normalized temperature, which is a relative temperature. It represents the ratio of the current temperature to a certain reference temperature and is dimensionless.

[0134] The shale gas volume factor is calculated using the following formula:

[0135]

[0136] In the formula: p sc and T sc The pressure and temperature under standard conditions are represented by p and T, respectively, in MPa and K; the pressure and temperature of the shale gas reservoir are represented by p and T, respectively, in MPa and K; Z g is the shale gas compressibility factor, dimensionless;

[0137] Shale gas viscosity was calculated using a modified Dempsey model (Stockman FD, Dempsey JR, Preston F. Practical application of a two-dimensional numerical model for gas reservoir studies[J]. Journal of Petroleum Technology, 1967, 19(9): 1127-1136.).

[0138]

[0139] Where: μ g Shale gas viscosity, Pa·s;

[0140] μ1=(1.709×10 -5 -2.062×10 -6 γ g (1.8T+32)+8.118×10 -3 -6.15×10 -3 log(γ g (6)

[0141] a0=-2.4621182 a1=2.970547414 a2=-0.286264054 a3=0.008054205

[0142] a4=2.80860949 a5=-3.49803305 a6=0.36037302 a7=-0.01044324

[0143] a8=-0.793385648 a9=1.39643306 a 10 =-0.149144925 a 11 =0.004410155

[0144] a 12 =0.083938718 a 13 =-0.186408848 a 14 =0.020336788 a 15 = -0.000609579

[0145] The shale gas adsorption / desorption compressibility coefficient is calculated using the following formula:

[0146]

[0147] Where: C d ρ is the adsorption / desorption compressibility coefficient, dimensionless; s The density of shale is kg / m³. 3 V L m is the volume of liquid interacting with shale gas. 3 ;p L p is the pressure of the liquid, in Pa; p is the current shale gas pressure, in Pa.

[0148] Based on the graph showing the relationship between shale gas physical properties and shale gas reservoir pressure, a polynomial fitting method was used to fit the relationship between the above physical properties and shale gas reservoir pressure, so as to calculate the physical properties corresponding to the shale gas reservoir pressure at any given time. Among them, the relationship between shale gas compressibility factor and shale gas reservoir pressure is fitted according to equation (2), the relationship between shale gas compressibility coefficient and shale gas reservoir pressure is fitted according to equation (3), the relationship between shale gas volume coefficient and shale gas reservoir pressure is fitted according to equation (4), the relationship between shale gas viscosity and shale gas reservoir pressure is fitted according to equation (6), and the relationship between shale gas adsorption / desorption compressibility coefficient and shale gas reservoir pressure is fitted according to equation (7).

[0149] S103: Grid the shale reservoir according to the fracture type, and use deterministic and stochastic modeling methods to determine the location of fractures in the shale matrix.

[0150] Based on the reservoir characteristics of the study block, a block-centered grid was used for the shale matrix to precisely define the initiation, end, length, and orientation of fractures within the shale matrix. For the distribution characteristics of natural fractures, a stochastic modeling method was used to randomly assign the initiation and end coordinates of natural fractures, thereby determining their distribution characteristics, including length and orientation. Meanwhile, the characteristics of artificial fractures were determined using a deterministic modeling method based on the on-site hydraulic fracturing construction plan to identify the hydraulic fracture characteristics of the study area.

[0151] The coordinates of the starting point of the natural crack: (X) nf1 =[x s1 ,x s2 ...,x sn ],Y nf1 =[y s1 ,y s2 ...,y sn ])

[0152] The endpoint coordinates of the natural crack: (X) nf2 =[x t1 ,x t2 ...,x tn ],Y nf2 =[y t1 ,y t2 ...,y tn ])

[0153] Coordinates of the starting point of the artificial crack: (X) hf1 =[x hs1 ,x hs2 ...,x hsn ],Y hf1 =[y hs1 ,y hs2 ...,y hsn ])

[0154] The endpoint coordinates of the artificial crack: (X) hf2 =[x ht1 ,x ht2 ...,x htn ],Y hf2 =[y ht1 ,y ht2 ...,y htn ])

[0155] The starting and ending points of the random modeling method for natural cracks are as follows:

[0156] starting point:

[0157] X nf1 =D xnf +(x e -2*D xnf )*rand(1); Y hf1 =D ynf +(y e -2*D ynf )*rand(1)

[0158] end:

[0159] X nf2 =X nf1 +sign(rand(1)-0.5)*d xnf ;Y hf2 =Y hf1 +sign(rand(1)-0.5)*d ynf (8)

[0160] in:

[0161]

[0162] In the formula, rand(x) is a random function that generates a random seed, and sign(x) is a sign function.

[0163] S104: Based on the distribution of the matrix grid and the location of the fractures in the matrix grid, calculate the fracture segment length, equivalent distance, and conductivity parameters embedded in the shale matrix grid.

[0164] Based on the principles of numerical simulation of oil and gas reservoirs, the shale matrix system is sequentially numbered into grids, and the fracture system is further numbered sequentially according to the order of artificial fractures and natural fractures.

[0165] The length of the fracture segment embedded in the shale matrix grid is calculated using the following formula:

[0166]

[0167] In the formula, (xs1 ,y s1 ) and (x s2 ,y s2 ) are the coordinates of the two intersection points of the crack and the matrix grid line, respectively.

[0168] Based on the embedding relationship between the crack and the matrix, and the connection relationship between cracks, the equivalent distance is divided into three categories: a. the equivalent distance between the matrix mesh and the crack; b. the equivalent distance between adjacent numbered crack segments in the same crack; c. the equivalent distance between intersecting crack segments.

[0169] a. The equivalent distance from the matrix mesh to the crack is calculated using the following formula:

[0170]

[0171] In the above formula, S is the cross-sectional area of ​​the crack segment embedded in the matrix mesh, and x n It is the normal distance from the matrix mesh to the crack;

[0172] b. The equivalent distance between adjacent numbered crack segments within the same crack is calculated using the following formula:

[0173]

[0174] In the formula:

[0175] This represents the equivalent distance between the first segment of a crack with adjacent numbers in the same crack and the length of its embedment into the matrix mesh. The equivalent distance representing the length of another adjacent crack segment embedded into the matrix mesh within the same crack;

[0176] l f1 The length of the first adjacent numbered segment of a crack within the same crack embedded into the matrix mesh;

[0177] l f2 The length of another adjacent-numbered segment of the same crack embedded into the matrix grid;

[0178] l f1i After discretizing the first crack segment embedded in the matrix mesh at equal intervals, the distance from the midpoint of the i-th discrete segment to the intersection of adjacent crack segments of the same crack segment (calculated using the two-point distance formula) is N. f1 part;

[0179] l f2i After discretizing another crack segment embedded in the matrix mesh at equal intervals, the distance from the midpoint of the i-th discrete segment to the intersection of adjacent crack segments of the same crack segment (calculated using the two-point distance formula) is N. f2 part;

[0180] The number of discrete segments is based on the shorter of the two crack segments. The shorter crack has 100 discrete segments, and the longer crack has the following number of discrete segments:

[0181]

[0182] c. The equivalent distance between intersecting crack segments is calculated using the following formula:

[0183]

[0184] In the formula,

[0185] It is the equivalent distance from the i-th discrete segment in a certain intersecting crack segment to the intersection point;

[0186] It is the equivalent distance from the i-th small discrete segment in another intersecting crack segment to the intersection point;

[0187] N f1 d is the number of discrete segments, where the length of the first segment of the intersecting crack is embedded in the matrix, at equal intervals. f1i Let be the equivalent distance from the i-th discrete segment in a given intersecting crack segment to the intersection point;

[0188] N f2 d is the number of discrete segments, where the length of the other segment of the intersecting crack is embedded in the matrix, at equal intervals. f2i This is the equivalent distance from the i-th small discrete segment in another intersecting crack segment to the intersection point;

[0189] The number of discrete segments is based on the shorter crack among the intersecting cracks. The shorter crack has 200 discrete segments, and the longer crack has the following number of discrete segments:

[0190]

[0191] Based on the embedding relationship between the fracture and the matrix, and the connection relationship between fractures, the conductivity is divided into four categories: a. conductivity between matrix meshes; b. conductivity between matrix meshes and fracture meshes; c. conductivity of adjacent fracture segments within the same fracture; d. conductivity of intersecting fracture segments.

[0192] a. Matrix mesh and its conductivity:

[0193] Grid in xx direction:

[0194]

[0195] YY direction grid:

[0196]

[0197] In the formula, dx, dy, and dz are the mesh sizes in the x, y, and z directions, respectively;

[0198] b. Conductivity of matrix mesh and fracture mesh:

[0199]

[0200] In the above formula, d m-f k is the length of the crack embedded into the matrix mesh. m and k f These represent the permeability of the matrix and the cracks, respectively.

[0201] c. Conductivity of adjacent crack segments with the same crack number:

[0202]

[0203] d. Conductivity of intersecting fracture segments:

[0204]

[0205] In the formula k f1 and k f2 These represent the permeability of the intersecting fractures; w f1 and w f2 These represent the widths of the intersecting cracks.

[0206] S105: Construct the correspondence between shale matrix grid pressure, fracture grid pressure and conductivity respectively, and form a numerical model of multi-stage fractured horizontal wells in shale gas reservoirs that takes into account stress sensitivity.

[0207] The relationship between shale matrix grid pressure and conductivity is calculated using the following formula:

[0208]

[0209] The relationship between crack mesh pressure and conductivity is calculated using the following formula:

[0210]

[0211] By combining equations (19), (20), and the Peaceman well equation, a numerical model of a multi-stage fractured horizontal well in a shale gas reservoir considering stress sensitivity can be obtained.

[0212] S106: Substitute the obtained physical property parameters into the numerical model of a multi-stage fractured horizontal well in a stress-sensitive shale gas reservoir, and use a fully implicit method to solve the well production decline and cumulative production under two working systems: controlled pressure production and depressurized production. Then, plot the relationship curves between well production decline and cumulative production and production time.

[0213] Based on the relationship between grid pressure and conductivity of the matrix system and fracture system, the matrix equation set of multi-stage fractured horizontal wells in stress-sensitive shale gas reservoirs is obtained. The gas well production and cumulative production under two working systems, namely pressure release production and pressure control production, are calculated using a fully implicit method.

[0214] The method for characterizing the two working systems of controlled pressure production and depressurized production is calculated using the following formula:

[0215] The production conditions for depressurization production are: p w =p wf

[0216] The production conditions for pressure-controlled production are: p w =p wmax -D t *n

[0217] In the two formulas above, p wf The bottom-hole flowing pressure during depressurization production; p wmax To control the maximum flowing pressure at the bottom of the well during pressure production; D t The pressure level is denoted by n; n is the nth time step.

[0218] When using a fully implicit method to solve the problem, the convergence condition is:

[0219] Matrix equation:

[0220]

[0221] Crack equation:

[0222]

[0223] S107: Based on the plotted curves showing the relationship between production decline and cumulative production and production time of multi-stage fractured horizontal wells in stress-sensitive shale gas reservoirs, predict the production decline pattern and cumulative production change pattern of the studied shale gas wells.

[0224] We selected newly commissioned shale gas wells in the study area that were developed using natural energy, and analyzed the impact of fluid properties and fracture characteristics on the production and cumulative production of the shale gas wells. Based on the current development status of shale gas wells in the study area, we evaluated the initial production capacity of the shale gas wells and predicted the production decline and cumulative production variation patterns of the shale gas wells.

[0225] Example:

[0226] This embodiment specifically includes the following steps:

[0227] S101: Obtain physical property data of shale gas reservoir and shale gas, as well as operating system data of the gas wells under study.

[0228] Actual data from the research block were collected and organized. For the field implementation case, the number of shale gas wells with declining production to be predicted was one, a multi-stage fractured horizontal well named SGFHW01. The horizontal well in the area where the implementation case was located had 5 fractured stages, 20 natural fractures in the near-wellbore formation, and a stress sensitivity coefficient of 0.001359 MPa. -1 The total production time is T = 2000 days; multi-stage fracturing horizontal wells are either pressure-release production or pressure-controlled production, where the pressure-release level is p wf =10.3MPa, pressure control level is p wf =20MPa-0.05*n.

[0229] The initial formation pressure of the shale gas well is Pi = 26.2 MPa, the formation temperature is T = 333.15 K, and the shale density is 2.5 t / m³. 3 The pore compressibility coefficient of the shale matrix is ​​0.000435 MPa. -1 The shale reservoir thickness is 21.9 m; the matrix permeability and porosity are 0.0005 mD and 0.05, respectively; the permeability of artificial fractures and natural fractures are 27400 mD and 1000 mD, respectively; the porosity is 0.75 and 0.55, respectively; both are 0.001 m wide; the relative density of shale gas is 0.65; the Langmuir adsorption pressure and adsorption volume of shale gas are 18 m³ / s. 3 / t and 5MPa.

[0230] S102: Based on the physical property data of shale gas reservoirs, calculate the relationship and relationship chart between shale gas reservoir permeability, shale gas physical property parameters and shale gas reservoir pressure, and use the polynomial fitting method to fit the relationship between shale gas physical property parameters and shale gas reservoir pressure.

[0231] Based on the physical properties of shale gas reservoirs, the relationship between reservoir stress sensitivity and pressure is calculated using formula (1), and the relationship between fluid stress sensitivity and pressure is calculated using formulas (2)-(6). Figure 4 The compression factor shown is as follows: Figure 5 The compressibility factor, volume factor, fluid viscosity, and shale gas adsorption / desorption compressibility factor are shown; the compressibility factor, compressibility factor, volume factor, and shale gas viscosity within the fracture are solved using the same method.

[0232] S103: Grid the shale reservoir according to the fracture type, and use deterministic and stochastic modeling methods to determine the location of fractures in the shale matrix.

[0233] Based on the hydraulic pressure distribution observed in the field implementation case, the discrete number of the shale matrix grid was set to 127*153, and the grid size was 10*10. The distribution of artificial and natural fractures within the shale matrix grid was determined using formulas (8) and (9), such as... Figure 3 As shown; in Figure 3 In the diagram, red lines represent artificial fractures, blue lines represent natural fractures, and brownish-red lines represent horizontal wellbore. In this case, the half-length of the artificial fracture is 101m, and the maximum half-length of the natural fracture is 52m.

[0234] S104: Based on the distribution of the matrix grid and the location of the fractures in the matrix grid, calculate the fracture segment length, equivalent distance, and conductivity parameters embedded in the shale matrix grid.

[0235] The shale matrix grid is numbered sequentially. After the matrix grid is numbered, the fracture grid is numbered in the order of hydraulic fractures first and natural fractures second. After the matrix grid and fracture grid are numbered, the fracture segment embedded in the matrix grid is determined, the sequence number of the fracture segment is saved, and the length of the fracture segment embedded in the matrix grid is calculated using formula (10) and saved.

[0236] Based on the embedding of the crack system in the matrix mesh, the location of the crack intersection point and the number of the two intersecting cracks are found, and the mesh number of the crack segment corresponding to the intersection point is also found. In this case, the coordinates of the crack intersection point, the number of the intersecting crack and the corresponding crack segment mesh number are shown in Table 1.

[0237] Table 1 Summary of Crack Intersections (This Patent)

[0238]

[0239] Based on the matrix-crack embedding relationship, the connection and intersection relationship between cracks, the equivalent distance between the matrix and crack, the equivalent distance between adjacent crack segments within the same crack, and the equivalent distance between intersecting cracks are calculated using formulas (11)-(15), respectively.

[0240] Based on the above three classifications of equivalent distances and the relationship between adjacent matrix grids, the conductivity of adjacent matrix grids, the conductivity of matrix grid-fracture grid, the conductivity of adjacent fracture segments within the same fracture, and the conductivity of intersecting fractures are calculated using formulas (16)-(19), respectively.

[0241] S105: Construct the correspondence between shale matrix grid pressure, fracture grid pressure and conductivity respectively, and form a numerical model of multi-stage fractured horizontal wells in shale gas reservoirs that takes into account stress sensitivity.

[0242] Based on the relationship between matrix mesh, fracture mesh pressure and conductivity, a matrix equation system of conductivity, matrix mesh pressure and fracture mesh pressure is constructed using formulas (20)-(21).

[0243] S106: Substitute the obtained physical property parameters into the numerical model of a multi-stage fractured horizontal well in a stress-sensitive shale gas reservoir, and use a fully implicit method to solve for the well production decline and cumulative production under both controlled and depressurized production regimes, and plot the relationship curves between well production decline and cumulative production and production time.

[0244] The equations in S105 are solved using a fully implicit method. The pressures of the matrix grid and fracture grid obtained are substituted into step S102 to recalculate the physical properties of the shale reservoir and fluid until the convergence conditions (22) and (23) are met, at which point the calculation ends.

[0245] Calculate the production and cumulative production of the shale gas well during pressure release and pressure control production in step S101, and plot the production decline and cumulative production graphs, as follows. Figure 6 As shown;

[0246] The effects of different pressure control levels and different pressure release levels on production decline and cumulative production are compared separately, as follows: Figure 7 and Figure 8 As shown, where:

[0247] Different pressure control levels: P wf =20MPa-D t *n;D t =0.02MPa,0.04MPa,0.06MPa,0.08MPa

[0248] Different pressure relief levels: P wf =5MPa, 10MPa, 15MPa, 20MPa

[0249] The effects of different stress sensitivity levels in shale reservoirs and different distribution characteristics of natural fractures in shale gas reservoirs on production decline and cumulative production were calculated separately, as follows: Figure 9 and Figure 10 As shown, where:

[0250] Shale reservoirs are sensitive to stress at different levels:

[0251] γ = 10 -2 MPa -1 5×10 -2 MPa -1 10 -1 MPa -1 0MPa -1

[0252] Different distribution characteristics of natural fractures in shale gas reservoirs: n f =30,60,90,120; n f Number of natural cracks

[0253] S107: Based on the plotted curves showing the relationship between production decline and cumulative production and production time of multi-stage fractured horizontal wells in stress-sensitive shale gas reservoirs, predict the production decline pattern and cumulative production change pattern of the studied shale gas wells.

[0254] Substitute the parameters from steps S101 and S103 back into step S109 to calculate the production decline, and compare the result with the actual production data of the study block. The comparison result is as follows: Figure 11 As shown. From Figure 11 As can be seen, the predicted decrease in output in this patented solution is basically consistent with the actual situation.

Claims

1. A method for predicting multi-fractured horizontal wells in shale gas reservoirs considering stress sensitivity, characterized in that, The method comprises the following steps: acquiring physical parameter data of a shale gas reservoir and shale gas and working system data of a studied gas well; based on the physical data of the shale gas reservoir, calculating the relationship and graph between the permeability of the shale gas reservoir, the physical parameters of the shale gas and the pressure of the shale gas reservoir, and fitting the relationship between the physical parameters of the shale gas and the pressure of the shale gas reservoir by using a polynomial fitting method; dividing a grid according to the fracture type in the shale reservoir, and determining the position of the fracture in the shale matrix by using a deterministic modeling method and a random modeling method; according to the distribution of the matrix grid and the position of the fracture in the matrix grid, respectively calculating the length of the fracture segment embedded in the matrix grid of the shale, the equivalent distance and the conductivity parameter; respectively constructing the corresponding relationship between the matrix grid pressure, the fracture grid pressure and the conductivity, and forming a numerical model of a multi-fractured horizontal well of a stress-sensitive shale gas reservoir; putting the acquired physical parameters into the numerical model of the multi-fractured horizontal well of the stress-sensitive shale gas reservoir, and respectively solving the gas well production decline and cumulative production under two working systems of pressure control production and pressure release production by using a fully implicit method, and drawing a relationship curve between the gas well production decline and cumulative production and production time; according to the drawn relationship curve between the production decline and cumulative production of the multi-fractured horizontal well of the stress-sensitive shale gas reservoir and the production time, predicting the production decline law and the cumulative production change law of the studied shale gas well.

2. The method for predicting stress-sensitive shale gas reservoir multi-stage fractured horizontal wells according to claim 1, wherein, The relationship between the permeability of the shale gas reservoir and the reservoir pressure is calculated by the following formula: where: K mi is the initial permeability of the shale gas reservoir, mD; γ is the stress sensitivity, MPa -1 ; P i is the initial formation pressure, MPa; p is the pressure of the shale gas reservoir, MPa.

3. The method for predicting stress-sensitive shale gas reservoir multi-fractured horizontal wells according to claim 1, characterized in that, The shale gas physical properties include a shale gas compression factor, a shale gas compressibility, a shale gas volume coefficient, shale gas viscosity and a shale gas adsorption / desorption compressibility; The shale gas compression factor is solved by iteration method through the following two Z g expressions, and then the solved y value is substituted into any one Z g equation for calculation: wherein: Z is the dimensionless compressibility factor of the shale gas; t is the inverse of the corresponding temperature, dimensionless; y is the corresponding density of the shale gas, dimensionless. g Z is the dimensionless compressibility factor of the shale gas; t is the inverse of the corresponding temperature, dimensionless; y is the corresponding density of the shale gas, dimensionless. The relationship between the shale gas compressibility and the reservoir pressure is calculated by the following formula: wherein: Z is a dimensionless number representing the corresponding pressure; g Z is a dimensionless number representing the shale gas compressibility factor; Tn is a normalized temperature, which is a relative temperature, representing the ratio of the current temperature to a reference temperature, and is dimensionless. The relationship between the shale gas volume coefficient and the reservoir pressure is calculated by the following formula: where: p sc and T sc are pressure and temperature, respectively, at standard conditions, MPa, K; p and T are pressure and temperature, respectively, of the shale gas reservoir, MPa, K; Z g is the shale gas compressibility factor, dimensionless; The relationship between the shale gas viscosity and the reservoir pressure is calculated by using a modified Dempsey model: where: μ g is the shale gas viscosity, Pa-s; μ1 = (1.709 x 10 -5 -2.062 x 10 -6 γ g )(1.8T + 32) + 8.118 x 10 -3 -6.15 x 10 -3 log(γ g ) a0=-2.4621182a1=2.970547414a2=-0.286264054a3=0.008054205 a4=2.80860949a5=-3.49803305a6=0.36037302a7=-0.01044324 a8 = -0.793385648 a9 = 1.39643306 a 10 = -0.149144925 a 11 = 0.004410155 a 12 = 0.083938718a 13 = -0.186408848a 14 = 0.020336788a 15 = -0.000609579 The relationship between the shale gas adsorption / desorption compressibility and the reservoir pressure is calculated by the following formula: where C d is the adsorption / desorption compressibility, dimensionless; p s is the pressure, Pa; p is the current shale gas pressure, Pa. 3 ; V L is the volume of the liquid interacting with the shale gas, m 3 ; p L is the pressure of the liquid, Pa; p is the current shale gas pressure, Pa.

4. The method for predicting stress-sensitive shale gas reservoir multi-fractured horizontal wells according to claim 1, characterized in that, The deterministic modeling method determines the position of the artificial fracture in the shale matrix, and the random modeling method determines the position of the natural fracture in the shale matrix, and the natural fracture random modeling method is as follows: a natural fracture starting point: X nf1 = D xnf + (x e - 2 * D xnf ) * rand(1); Y hf1 = D ynf + (y e - 2 * D ynf ) * rand(1) a natural fracture ending point: X nf2 = X nf1 + sign(rand(l)-0.5)*d xnf ; Y hf2 = Y hf1 + sign(rand(l)-0.5)*d ynf wherein rand(x) is a random function generating a random seed, and sign(x) is a sign function. The length of the fracture segment embedded in the matrix grid of the shale is calculated by the following formula:

5. The method for predicting stress-sensitive shale gas reservoir multi-fractured horizontal wells according to claim 1, characterized in that, The equivalent distance is divided into three categories: a. the equivalent distance between the matrix grid and the fracture; b. the equivalent distance between adjacent numbered fracture segments in the same fracture; and c. the equivalent distance between intersecting fracture segments; where (x s1 ,y s1 ) and (x s2 ,y s2 ) are the coordinates of the two intersection points of the fracture and the matrix grid lines, respectively.

6. The method for predicting stress-sensitive shale gas reservoir multi-fractured horizontal wells according to claim 1, characterized in that, a. the equivalent distance between the matrix grid and the fracture is calculated by the following formula: b. the equivalent distance between adjacent numbered fracture segments in the same fracture is calculated by the following formula: S is the cross-sectional area of the fracture segment embedded in the matrix grid, x n is the normal distance from the matrix grid to the fracture; wherein ​ represents the equivalent distance of embedding of the first segment of the fracture of adjacent number in the same fracture into the length of the matrix grid; represents the equivalent distance of embedding of the other segment of the fracture of adjacent number in the same fracture into the length of the matrix grid; l f1 the length of the first segment of the fracture of the same fracture with the adjacent number is embedded into the matrix grid; l f2 the length of the other segment of the fracture of the same fracture with the adjacent number is embedded into the matrix grid; l f1i After the first segment of the fracture is discretized equidistantly and embedded into the matrix grid, the distance between the midpoint of the ith discrete segment and the intersection point of the adjacent fracture segments of the same fracture segment (calculated by the two-point distance formula), which is the N f1 segment. l f2i After equidistant discretization for another fracture segment embedded in the matrix grid, the distance from the midpoint of the ith discrete segment to the intersection point of the adjacent fracture segments of the same fracture segment (calculated by the two-point distance formula), which is N f2 segment; The discrete number of segments is based on the shorter of the two intersecting fractures, with the shorter fracture having a discrete number of 100 segments and the longer fracture having a discrete number of: The equivalent distance of the intersecting fracture segments is calculated using the following equation: where, is the equivalent distance of the i-th small discrete segment in a certain intersecting fracture segment to the intersection point; is the equivalent distance of the i-th small discrete segment in another intersecting fracture segment to the intersection point; N f1 the number of discrete segments for equidistant discretization of the length of the fracture embedded into the matrix for the first segment of the intersecting fracture, d f1i the equivalent distance of the i-th small discrete segment in a certain intersecting fracture segment to the intersection point; N f2 the number of discrete segments for equidistant discretization of the length of the fracture embedded in the matrix for another segment of the intersecting fracture, d f2i the equivalent distance of the i-th small discrete segment in another intersecting fracture segment to the intersection point; The discrete number of segments is based on the shorter of the two intersecting fractures, with the shorter fracture having a discrete number of 200 segments and the longer fracture having a discrete number of:

7. The method for predicting stress-sensitive shale gas reservoir multi-fractured horizontal wells according to claim 1, characterized in that, The conductivities are divided into four categories: a. matrix grid to matrix grid conductivities; b. matrix grid to fracture grid conductivities; c. fracture grid to fracture grid conductivities; and d. intersecting fracture segment conductivities. a. matrix grid to matrix grid conductivities: x-x direction grid: y-y direction grid: where dx, dy, and dz are the grid sizes in the x, y, and z directions, respectively. b. matrix grid to fracture grid conductivities: c. fracture grid to fracture grid conductivities: d in the above equation m-f is the length of the fracture embedded in the matrix grid, k m and k f are the permeability of the matrix and fracture, respectively; d. intersecting fracture segment conductivities: The relationship between the shale matrix grid pressure and the conductivities is calculated using the following equation: where k f1 and k f2 are the permeabilities of the intersecting fractures; w f1 and w f2 are the widths of the intersecting fractures.

8. The method for predicting stress-sensitive shale gas reservoir multi-fractured horizontal wells according to claim 1, characterized in that, The relationship between the fracture grid pressure and the conductivities is calculated using the following equation: The numerical model for a multi-fractured horizontal well in a stress-sensitive shale gas reservoir is obtained by simultaneously solving the above equations for the shale matrix grid pressure and the conductivities, the fracture grid pressure and the conductivities, and the Peaceman well equation. The method for the two work systems of pressure-controlled production and pressure-released production is calculated using the following equation:

9. The method for predicting stress-sensitive shale gas reservoir multi-fractured horizontal wells according to claim 1, wherein, The convergence condition for the fully implicit method is: The production conditions for the pressure relief production are: p w = p wf The production conditions for the controlled pressure production are: p w = p wmax - D t *n p wf is the bottom hole flowing pressure during open pressure production; p wmax is the maximum bottom hole flowing pressure during pressure controlled production; D t is the pressure control level; n is the nth time step.

10. The method for predicting stress-sensitive shale gas reservoir multi-fractured horizontal wells according to claim 1, wherein, Matrix equation: Fracture equation: ​

Images