AU2020103953A4

AU2020103953A4 - Method and system for predicting production of fractured horizontal well in shale gas reservoir

Info

Publication number: AU2020103953A4
Application number: AU2020103953A
Authority: AU
Inventors: Shuyong Hu; Wenhai HUANG; Tingting Qiu; Daqian RAO; Jiayi Zhang
Original assignee: Southwest Petroleum University
Current assignee: Southwest Petroleum University
Priority date: 2020-12-08
Filing date: 2020-12-08
Publication date: 2021-02-11
Anticipated expiration: 2028-12-08

Abstract

The present disclosure relates to a method and system for predicting the production of a fractured horizontal well in a shale gas reservoir, and relates to the technical field of fractured horizontal wells. Considering the different diffusion modes of the matrix in different zones, the present disclosure uses Fick's First Law to describe the pseudo-steady-state diffusion of the matrix in the fracture network zone, Fick's Second Law to describe the unsteady-state diffusion of the matrix in the pure matrix zone, and Darcy's Law to describe the seepage in the fracture network. The present disclosure predicts the production of the fractured horizontal well in the shale gas reservoir under the conditions of matrix-microfracture coupling and hydraulically created fracture-microfracture coupling. The present disclosure improves the prediction accuracy of shale gas well production, and more accurately describes the actual flow law of the shale gas reservoir.

Description

METHOD AND SYSTEM FOR PREDICTING PRODUCTION OF FRACTURED HORIZONTAL WELL IN SHALE GAS RESERVOIR TECHNICAL FIELD

The present disclosure relates to the technical field of fractured horizontal wells, in particular to a method and system for predicting the production of a fractured horizontal well in a shale gas reservoir.

BACKGROUND

At present, many studies at home and abroad are based on the apparent permeability model that considers the multiple flow mechanisms of matrix pores to study the unsteady-state flow or pseudo-steady-state flow in the matrix and the Darcy flow or non-Darcy flow in the fracture and establish the unsteady-state or pseudo-steady-state productivity model of matrix-fracture coupling. In the nano-scale pores of the matrix, gas diffusion based on the concentration difference is an important component of gas flow. Shale gas diffusion can be described in two different ways: pseudo-steady-state diffusion and unsteady-state diffusion. In the hydraulically fractured zone and the pure matrix zone, the size of the matrix block is very different, and the diffusion mode of gas in the matrix is also different. Therefore, it is difficult for the model to describe the true gas flow by considering only the pseudo-steady-state diffusion or unsteady-state diffusion of the matrix, which leads to the low accuracy of the model in predicting the production of the shale gas well.

Any discussion of documents, acts, materials, devices, articles or the like which has been included in the present specification is not to be taken as an admission that any or all of these matters form part of the prior art base or were common general knowledge in the field relevant to the present disclosure as it existed before the priority date of each of the appended claims.

SUMMARY

Throughout this specification the word "comprise", or variations such as "comprises" or "comprising", will be understood to imply the inclusion of a stated element, integer or step, or group of elements, integers or steps, but not the exclusion of any other element, integer or step, or group of elements, integers or steps.

Some embodiments of the present disclosure aim to provide a method and system for predicting the production of a fractured horizontal well in a shale gas reservoir, so as to improve the prediction accuracy of the production of a shale gas well.

According to the present disclosure, there is provided a method for predicting the production of a fractured horizontal well in a shale gas reservoir includes:

dividing a fractured horizontal well to be predicted in a shale gas reservoir into five seepage zones according to a matrix block and a fracture network after hydraulic fracturing, where the five seepage zones include: hydraulically fractured zone I, fracture network zone II, pure matrix zone III, pure matrix zone IV and pure matrix zone V;

obtaining a zone I seepage differential equation for hydraulically created fracture zone I, a zone II seepage differential equation and a zone II diffusion equation for fracture network zone II, a zone III diffusion equation for pure matrix zone III, a zone IV diffusion equation for pure matrix zone IV and a zone V diffusion equation for pure matrix zone V;

obtaining a preset dimensionless transform relationship;

solving the zone V diffusion equation by the dimensionless transform relationship and Laplace transform to obtain a solution of the zone V diffusion equation;

solving the zone IV diffusion equation by the dimensionless transform relationship and Laplace transform to obtain a solution of the zone IV diffusion equation;

solving the zone III diffusion equation by the dimensionless transform relationship, Laplace transform and the solution of the zone V diffusion equation to obtain a solution of the zone III diffusion equation;

solving the zone II seepage differential equation by the dimensionless transform relationship, Laplace transform, the zone II diffusion equation, the solution of the zone IV diffusion equation and the solution of the zone III diffusion equation to obtain a solution of the zone II seepage differential equation;

solving the zone I seepage differential equation by the dimensionless transform relationship, Laplace transform and the solution of the zone II seepage differential equation to obtain a solution of the zone I seepage differential equation;

obtaining a first preset condition;

using the solution of the zone I seepage differential equation to obtain a bottom hole pseudo-pressure solution according to the first preset condition;

obtaining a dimensionless production solution by Duhamel's Principle according to the bottom hole pseudo-pressure solution; and predicting the production of the fractured horizontal well in the shale gas reservoir by using a Stehfest numerical inversion method according to the dimensionless production solution.

According to the present disclosure, there is further provided a system for predicting the production of a fractured horizontal well in a shale gas reservoir includes:

a seepage zone dividing module, configured to divide a fractured horizontal well to be predicted in a shale gas reservoir into five seepage zones according to a matrix block and a fracture network after hydraulic fracturing, where the five seepage zones include: hydraulically fractured zone I, fracture network zone II, pure matrix zone III, pure matrix zone IV and pure matrix zone V;

a seepage zone equation obtaining module, configured to obtain a zone I seepage differential equation for hydraulically created fracture zone I, a zone II seepage differential equation and a zone II diffusion equation for fracture network zone II, a zone III diffusion equation for pure matrix zone III, a zone IV diffusion equation for pure matrix zone IV and a zone V diffusion equation for pure matrix zone V;

a dimensionless transform relationship obtaining module, configured to obtain a preset dimensionless transform relationship;

a zone V diffusion equation solving module, configured to solve the zone V diffusion equation by the dimensionless transform relationship and Laplace transform to obtain a solution of the zone V diffusion equation;

a zone IV diffusion equation solving module, configured to solve the zone IV diffusion equation by the dimensionless transform relationship and Laplace transform to obtain a solution of the zone IV diffusion equation;

a zone III diffusion equation solving module, configured to solve the zone III diffusion equation by the dimensionless transform relationship, Laplace transform and the solution of the zone V diffusion equation to obtain a solution of the zone III diffusion equation;

a zone II seepage differential equation solving module, configured to solve the zone II seepage differential equation by the dimensionless transform relationship, Laplace transform, the zone II diffusion equation, the solution of the zone IV diffusion equation and the solution of the zone III diffusion equation, to obtain a solution of the zoneII seepage differential equation; a zone I seepage differential equation solving module, configured to solve the zone I seepage differential equation by the dimensionless transform relationship, Laplace transform and the solution of the zone II seepage differential equation to obtain a solution of the zone I seepage differential equation; a first preset condition obtaining module, configured to obtain a first preset condition; a bottom hole pseudo-pressure solution obtaining module, configured to use the solution of the zone I seepage differential equation to obtain a bottom hole pseudo-pressure solution according to the first preset condition; a dimensionless production solution obtaining module, configured to obtain a dimensionless production solution by Duhamel's Principle according to the bottom hole pseudo-pressure solution; and a production predicting module, configured to predict the production of the fractured horizontal well in the shale gas reservoir by using a Stehfest numerical inversion method according to the dimensionless production solution.

According to the specific embodiments of the present disclosure, the present disclosure has the following technical effects.

The present disclosure provides a method and system for predicting the production of a fractured horizontal well in a shale gas reservoir. Considering the different diffusion modes of the matrix in different zones, the present disclosure uses Fick's First Law to describe the pseudo-steady-state diffusion of the matrix in the fracture network zone, Fick's Second Law to describe the unsteady-state diffusion of the matrix in the pure matrix zone, and Darcy's Law to describe the seepage in the fracture network. The present disclosure predicts the production of the fractured horizontal well in the shale gas reservoir under the conditions of matrix-microfracture coupling and hydraulically created fracture-microfracture coupling. The present disclosure improves the prediction accuracy of shale gas well production, more accurately describes the actual flow law of the shale gas reservoir, provides a new method for production decline analysis and prediction of the shale gas well, and provides theoretical basis and production suggestions for the beneficial development of the shale gas reservoir.

BRIEF DESCRIPTION OF DRAWINGS

To describe the technical solutions in the embodiments of the disclosure or in the prior art more clearly, the accompanying drawings required for the embodiments are briefly described below. Apparently, the accompanying drawings in the following description show merely some embodiments of the present disclosure, and a person of ordinary skill in the art may still derive other accompanying drawings from these accompanying drawings without creative efforts.

FIG. 1 is a flowchart of a method for predicting the production of a fractured horizontal well in a shale gas reservoir according to an embodiment of the present disclosure.

FIG. 2 is a simplified schematic diagram of a seepage pattern of the fractured horizontal well in the shale gas reservoir according to the embodiment of the present disclosure.

FIG. 3 is a schematic diagram of seepage zones of the fractured horizontal well in the shale gas reservoir according to the embodiment of the present disclosure.

FIG. 4 is a schematic diagram showing an effect of the size of a simulated reservoir volume (SRV) zone on the production according to the embodiment of the present disclosure.

FIG. 5 is a schematic diagram showing an effect of a wellbore storage coefficient on the production according to the embodiment of the present disclosure.

FIG. 6 is a schematic diagram showing an effect of a skin coefficient on the production according to the embodiment of the present disclosure.

FIG. 7 is a schematic diagram showing an effect of an inter-porosity flow coefficient on the production according to the embodiment of the present disclosure.

FIG. 8 is a schematic diagram showing an effect of fracture conductivity on the production according to the embodiment of the present disclosure.

FIG. 9 is a schematic diagram showing an effect absorbed gas and free gas in a matrix on the production according to the embodiment of the present disclosure.

DETAILED DESCRIPTION

The technical solutions in the embodiments of the present disclosure are clearly and completely described below with reference to the accompanying drawings in the embodiments of the present disclosure. Apparently, the described embodiments are merely a part rather than all of the embodiments of the present disclosure. All other embodiments obtained by a person of ordinary skill in the art based on the embodiments of the present disclosure without creative efforts should fall within the protection scope of the present disclosure.

An objective of the present disclosure is to provide a method and system for predicting the production of a fractured horizontal well in a shale gas reservoir, so as to improve the prediction accuracy of the production of a shale gas well.

To make the foregoing objective, features, and advantages of the present disclosure clearer and more comprehensible, the present disclosure is further described in detail below with reference to the accompanying drawings and specific embodiments.

This embodiment provides a method for predicting the production of a fractured horizontal well in a shale gas reservoir. FIG. 1 shows a flowchart of the method for predicting the production of a fractured horizontal well in a shale gas reservoir according to the embodiment of the present disclosure. As shown in FIG. 1, the method for predicting the production of a fractured horizontal well in a shale gas reservoir includes:

Step 101: Divide a fractured horizontal well to be predicted in a shale gas reservoir into five seepage zones according to a matrix block and a fracture network after hydraulic fracturing, where the five seepage zones include: hydraulically fractured zone I, fracture network zone II, pure matrix zone III, pure matrix zone IV and pure matrix zone V.

The shale reservoir has very low permeability, and it is believed that the effective discharge boundary after hydraulic fracturing is equal or close to the length of a hydraulically created fracture. Therefore, the fractured horizontal well in the shale gas reservoir can be simulated as a horizontal well with a rectangular drainage area. The rectangular drainage area of the horizontal well is composed of a fracture network divided by a matrix block. The pure matrix zone without being hydraulically fractured also contributes to the production of shale gas and cannot be ignored. The seepage pattern of the fractured horizontal well in the shale gas reservoir is shown in FIG. 2. A single-stage treatment zone consists of multiple single-cluster treatment zones. In FIG. 2, there are three single-cluster treatment zones in a single-stage treatment zone, and the single-cluster treatment zone is composed of a hydraulically fractured network and an untreated zone. A quarter of a fracture network (a single-cluster treatment zone) in FIG. 2 is taken as the research object to establish a five-linear-flow mathematical model. The five-linear-flow mathematical model includes five seepage zones: hydraulically fractured zone (zone I), fracture network zone (zone II), and pure matrix zones (zone III, zone IV and zone V), as shown in FIG. 3. The gas flows into an artificial fracture from a formation, and then flows into a wellbore along the fracture. In FIG. 3, the one-way arrow indicates a mass transfer direction (flow direction and diffusion direction) of the gas. In FIG. 3, the horizontal axis y indicates the direction of the horizontal wellbore; the vertical axis x indicates the direction of the hydraulically created fracture. ye indicates the width of the research object, which is one-half of the cluster spacing. yl indicates the width of the fracture network zone (i.e. stimulated reservoir volume (SRV) zone). In FIG. 3, zone II is the SRV zone. xe is a distance from the boundary of the untreated zone to the horizontal wellbore (which can be obtained from microseismic monitoring). XF is the length of the hydraulically created fracture in zone I. WF is the width of the hydraulically created fracture, unit: m. The length of the hydraulically created fracture in zone I is equal to the length of those in zones II and III.

Assumptions are made for the five-linear flow mathematical model. (1) The reservoir is a closed reservoir of equal thickness and is always in an isothermal state during the mining process. (2) The gas flow in the fracture obeys the Darcy flow law, and is single-phase gas flow, ignoring the influence of gravity and capillary force. (3) The hydraulically created fractures are perpendicular to the wellbore, symmetrically and evenly distributed, with the same properties and characteristics. (4) The fracture network zone (SRV zone) has a fracture-pore medium, and the outside zones, that is, the zones out of the SRV zone, specifically, zone III, zone IV and zone V, all have a single porous medium with the same physical properties. (5) The gas reservoir is in equilibrium before exploitation, and the adsorbed gas and free shale gas are also in dynamic equilibrium.

Step 102: Obtain a zone I seepage differential equation for hydraulically created fracture zone I, a zone II seepage differential equation and a zone II diffusion equation for fracture network zone II, a zone III diffusion equation for pure matrix zone III, a zone IV diffusion equation for pure matrix zone IV and a zone V diffusion equation for pure matrix zone V.

Step 102 may specifically include:

After hydraulic fracturing, the main hydraulically created fracture in hydraulically fractured zone I is perpendicular to the horizontal wellbore, the fracture height is equal to the height of the reservoir, and the seepage direction is perpendicular to the wellbore direction. The fluid (gas source) in the hydraulically created fracture zone comes from the SRV zone, and there is no fluid supply by itself, that is, the source term q=0. Substituting Darcy's formula into the material balance equation leads to the zone I differential equation:

(OFPlF)-V(p1FkIFVp1F)=O

at p (1)

In Eq. (1), 1F represents a fracture porosity in zone I; 1F represents a fracture gas density in zone I; t represents time;V represents a gradient operator; kF represents fracture permeability in zone I; represents a gas viscosity of the fractured horizontal well in the shale gas reservoir; VPlF represents a gradient operator of fracture pressure in zone I.

The shale matrix diffusion process is a process from unsteady-state diffusion to pseudo-steady-state diffusion. The matrix size in the SRV zone is relatively small (mostly on the order of meters in diameter), and the pressure drop generated by diffusion can quickly spread to the inside of the bedrock and enter a pseudo-steady state. Therefore, Fick's First Law is used to describe the diffusion process from the matrix in fracture network zone II to the fracture, to obtain the zone II diffusion equation:

8V m *.Dm(V-Vm) (2) at In the equation, Vm represents the total concentration of shale gas inside the matrix block, including the concentration of adsorbed phase and free phase; 0~ represents a matrix block shape factor; Dm represents a coefficient of shale gas diffusion inside the matrix block; VE represents an apparent concentration of shale gas when the gas supply from the matrix block to the fracture network reaches equilibrium.

Vm Vm= c PmsPm+(1-m) c PmVL Vm is expressed by: se P +PL (3)

VE is expressed by: VE Z S CZ,, sc +pf Gf P('m) VL (4) ZTP Pf+PL

In the equations, 0. represents a matrix porosity; Z represents a gas compressibility under standard conditions; Tsc represents a temperature under standard conditions, 273.15k; Z represents a gas compressibility under gas reservoir conditions; T represents an actual temperature; Psc represents a standard pressure, 1.01x1 5 Pa; Pm represents a matrix pressure; VL represents a Langmuir volume; PL represents a Langmuir pressure; and f represents a fracture pressure.

In fracture network zone II, the matrix supplies gas to the fracture, that is, the source term is q2m. According to Darcy's Law and the principle of material balance, the seepage differential equation of the fracture in zone II is:

2f2f Pgsc (5) 2m at p

In the equation, 2f represents a fracture porosity in zone II; 2f represents a fracture gas density in zone II; k2f represents a fracture permeability in zone II; represents a gas viscosity of the fractured horizontal well in the shale gas reservoir; VP 2 represents a gradient operator of fracture pressure in zone I; Pg"c represents a gas density under standard conditions; q2m represents the amount of gas supplied to the fracture by the matrix.

The fracturing treatment has little effect on pure matrix zone III, pure matrix zone IV and pure matrix zone V, so they are regarded as pure matrix zones. The matrix gas in zone III, zone IV and zone V is supplied to the fracture network in a diffusion manner. Assuming that Vm (x, y, z, t) is the gas concentration at a certain point in the matrix zone at a certain moment, x, y and z respectively represent the x direction, y direction and z direction in a three-dimensional (3D) Cartesian coordinate system of the fractured horizontal well in the shale gas reservoir. In this embodiment, there is no gas mass transfer in the z direction, which is quantitative and will not be considered in model establishment and solution. Due to the large size of the outer matrix zones, the diffusion in the matrix zone is considered to be unsteady-state diffusion. Vm (x, y, z, t) is a function of time and space. Therefore, Fick's Second Law is used to describe the diffusion process from the matrix to the fracture in zones III, IV and V, and the zone III diffusion equation, the zone IV diffusion equation and the zone V diffusion equation are uniformly written as: 8V 2 Vm (6) mD at 2 In the equation, v is defined as the square of the gradient operator, meaning the second-order partial derivative of Vm to x and y.

Step 103: Obtain a preset dimensionless transform relationship. For the convenience of solving, dimensionless variables are defined, and the preset dimensionless transform relationships include:

V rkfhT,

( Dimensionless pseudo-pressure: P s ; dimensionless time:

D XD YD WFD WF XF ; dimensionless distance: Xf Xi; dimensionless

7F 7F F 7FD conductivity factor: 7/ , where: f cfC p#CF dimensionless kF W D FD _F F DmD conductivity: kfxF dimensionless diffusion coefficient: storage

co= D. A - F .D coefficient: Z ,q,, inter-porosity flow coefficient: 7; dimensionless

concentration: ED =VE i VmD =m i

In the above equations, V'i represents a dimensionless fracture pseudo-pressure of an i-th zone, where i=1 represents zone I, i=2 represents zone II, i=3 represents zone III, i=4 represents zone IV, and i=5 represents V zone. There is no subscript i in the SRV model (zone II model). 7 k represents the circumference ratio. I represents the fracture permeability of zone II, unit: m2. h represents the thickness of the formation, unit: m. qsc represents the surface production of the 3 horizontal well, unit: m /s. O represents the pseudo-pressure of the entire gas reservoir at an initial moment. Vl represents the pseudo-pressure of the i-th zone at a certain moment. tD

represents a dimensionless time. 1I represents a dimensionless conductivity factor. XF

represents the length of the hydraulically created fracture, unit: m. XD represents a dimensionless variable x, where x represents a variable in the x direction in the 3D coordinate system. xf represents the fracture length of zone II. YD represents a dimensionless variable y, where y represents a variable in the y direction in the 3D coordinate system. WFD represents a dimensionless width of the hydraulically created fracture. WF represents the width of the hydraulically created fracture, unit: m. FD represents a dimensionless conductivity factor of the fracture in zone I. 7F represents the conductivity factor of the main fracture in zone I. 7I represents the fracture conductivity factor of zone II. / represents the porosity of the fracture network. represents thcmrsiiiyofatrznI.kF the compressibility of fracture zone II. F represents the permeability of zone I. OF represents the porosity of the hydraulically created fracture. CF represents the compressibility of zone I, unit: I/Pa. FCD represents a dimensionless conductivity. DmD represents a dimensionless diffusion coefficient. D. represents the diffusion coefficient, unit: m 2 /s. o represents the storage coefficient. X represents the inter-porosity flow coefficient. ED represents a dimensionless apparent concentration of shale gas when the gas supply from the matrix block to the fracture network reaches equilibrium. Vi represents an initial concentration of shale gas inside the matrix block. VmD represents a dimensionless total concentration of shale gas inside the matrix block. The subscript D represents dimensionless. The subscript f represents the fracture network. represents porosity. The subscript m represents the matrix. The subscript F represents the hydraulically created fracture. The subscript 0 represents an initial state. Vim represents a concentration in different zones, i=1, 2, 3, 4, 5; Vim represents a shale gas concentration in zone I, including the concentration of adsorbed phase and free phase.

The initial concentration of shale gas inside the matrix block is expressed as:

#mZcT p( ) PlVL I ZTpSC pA + PL In the equation, A represents an original formation pressure.

Step 104: Solve the zone V diffusion equation by the dimensionless transform relationship and Laplace transform to obtain a solution of the zone V diffusion equation.

Step 104 may specifically include:

Dimensionless transform is performed on the zone V diffusion equation by the dimensionless transform relationship to obtain a zone V dimensionless diffusion equation.

Laplace transform is performed on the zone V dimensionless diffusion equation to obtain a zone V dimensionless diffusion equation in the Laplace space. The diffusion direction of zone V is the -x direction of the Cartesian coordinate system, then the zone V dimensionless diffusion equation is subjected to Laplace transform based on the dimensionless time tD to obtain a zone V dimensionless diffusion equation in the Laplace space:

a2v 65mD = ___V

2 5mD D D5mD (7)

In the equation, 5mD represents a dimensionless matrix gas concentration of zone V in the Laplace space; XD represents a dimensionless variable x; s is a Laplace variable, defined during

the Laplace transform; D5mD represents a dimensionless matrix diffusion coefficient of zone V; mD represents a dimensionless matrix gas concentration in zone V.

A boundary condition of pure matrix zone V is obtained. The zone V boundary condition

includes an outer boundary condition 05Dand an inner boundary condition 8 xD OD D

5m 3mDX where, l represents a dimensionless matrix gas concentration in

zone III in the Laplace space; XeD represents a dimensionless boundary length in the x direction;

V3mD represents a dimensionless matrix gas concentration in zone III.

The zone V boundary condition is used to solve the zone V dimensionless diffusion equation in the Laplace space to obtain a solution of the zone V diffusion equation. The solution of the zone V diffusion equation (the zone V matrix diffusion equation) is:

cosh[ s /D 5mD(eD XD)] V5mD 3m.D D cosh [s /D 5 mD X - 1)]

Step 105: Solve the zone IV diffusion equation by the dimensionless transform relationship and Laplace transform to obtain a solution of the zone IV diffusion equation.

Step 105 may specifically include:

Dimensionless transform is performed on the zone IV diffusion equation by the dimensionless transform relationship to obtain a zone IV dimensionless diffusion equation.

Laplace transform is performed on the zone IV dimensionless diffusion equation to obtain a zone IV dimensionless diffusion equation in the Laplace space.

A boundary condition of pure matrix zone IV is obtained. The zone IV boundary condition

aV4 mD 0

includes an outer boundary condition VD XD andaninnerboundarycondition

4mD 'XD- 21D 'XD 1 , where, 4mD represents a dimensionless matrix gas concentration in zone

IV in the Laplace space; V/1 represents a dimensionless fracture gas concentration in zone II in the Laplace space.

The zone IV boundary condition is used to solve the zone IV dimensionless diffusion equation in the Laplace space to obtain a solution of the zone IV diffusion equation. The solution of the zone IV diffusion equation (the zone IV matrix diffusion equation) is:

V4 mD v = V~j cosh [s / D 4mD (XeD XD)] 2fD -(XeD

In the equation, V 4 mD represents a dimensionless matrix gas concentration in zone IV in the

Laplace space; D4mD represents a dimensionless matrix gas diffusion coefficient of zone IV;

2fD represents a dimensionless gas concentration in the fracture network in zone II.

Step 106: Solve the zone III diffusion equation by the dimensionless transform relationship, Laplace transform and the solution of the zone V diffusion equation to obtain a solution of the zone III diffusion equation. Since the diffusion direction from the zone III matrix to the zone II fracture network is perpendicular to the diffusion direction from the zone V matrix to the zone III matrix, the diffusion of the zone III matrix is two-dimensional diffusion.

Step 106 may specifically include:

Dimensionless transform is performed on the zone III diffusion equation by the dimensionless transform relationship to obtain a zone III dimensionless diffusion equation.

Laplace transform is performed on the zone III dimensionless diffusion equation to obtain a zone III dimensionless diffusion equation in the Laplace space. The zone III dimensionless diffusion equation in the Laplace space is:

a2 V3mD DD 5mD av5mD _ I1 v3mD aYD 2 D3mD XD D1 D3mD tD

In the equation, YD represents a dimensionless variable y; D3mD represents a dimensionless matrix gas diffusion coefficient of zone III.

Derivative finding is performed on the solution of the zone V diffusion equation. The derivative finding on the solution of the zone V matrix diffusion equation, i.e. Eq. (8), is as follows:

VKrD VinD s /D5 , tanh [s/ D (xeD )] (11

To facilitate the solution below, a first function fs(s) is defined as:

f(S)= D mD D3mD 5 s /D5mD tanhFS/D 5 mD (XD )]± D3mD (12).

A boundary condition of pure matrix zone III is obtained. The zone III boundary condition

includes an outer boundary condition VmD 2fD D andannnerboundary YDYmD DDYD

condition 3D= 0, where, YlD represents a dimensionless width of zone II; YeD OY D =YeD

represents a dimensionless boundary length in the y direction.

According to the solution of the zone V diffusion equation after derivative finding, the zone III boundary condition is used to solve the zone III dimensionless diffusion equation in the Laplace space to obtain a solution of the zone III diffusion equation. By substituting the solution of the zone V diffusion equation (Eq. (11)) after derivative finding and the first functionfs(s) (Eq. (12)) into the zone III dimensionless diffusion equation (Eq. (10)) in the Laplace space, Eq. (10) is reduced to:

82V3mD -(13)

YD 2 3(s3mD

Solving Eq. (13) with the zone III boundary condition yields a solution of the zone III diffusion equation (the zone III matrix diffusion equation) as follows:

3 2fD cosh f 3 (S)(YeD YD)] (14) YD YlD coshL 3(S)(YeD YD)j

Step 107: Solve the zone II seepage differential equation by the dimensionless transform relationship, Laplace transform, the zone II diffusion equation, the solution of the zone IV diffusion equation and the solution of the zone III diffusion equation to obtain a solution of the zone II seepage differential equation.

Step 107 may specifically include:

Dimensionless transform is performed on the zone II diffusion equation by the dimensionless transform relationship to obtain a zone II dimensionless diffusion equation. Using the 2 A= F20-D dimensionless transform relationships ED E i mD m i and f to nondimensionalize the diffusion process from the matrix of zone II to the fracture (Eq. (2)) yields a-vmD V = A(VED VmD) a zone II dimensionless diffusion equation: D , where VmD represents a dimensionless total concentration of shale gas inside the matrix block.

Laplace transform is performed on the zone II dimensionless diffusion equation to obtain a zone II dimensionless diffusion equation in the Laplace space. A pseudo-steady-state diffusion equation is used to describe the difussion process from the matrix block to the fracture, and Laplace transform is performed on the zone II dimensionless diffusion equation to obtain a zone II dimensionless diffusion equation in the Laplace space:

VmD A ED (15)

In the equation, VmD represents a dimensionless matrix gas concentration in the Laplace space; VED represents a dimensionless gas concentration in the lapalce space when the matrix and the fracture are in equilibrium.

Dimensionless transform is performed on a pressure function of fracture network zone II by the dimensionless transform relationship to obtain a dimensionless shale gas concentration when the gas supply from the matrix block to fracture network zone II reaches equilibrium. The shale gas concentration VE when the gas supply from the matrix block to the fracture network reaches equilibrium is a function of the fracture network pressure p, which produces the following equation according to the definition of dimensionless variable:

- mZscqsc piZi pLVL pscqscT pZ D Zckh 2pI (pf +pL)(pi+pL)rckhT, 2 pi V (16)

In the equation, k represents a fracture permeability; 4 represents a gas viscosity in the initial state; Zi represents a gas deviation coefficient in the initial state; '/fD represents a dimensionless pseudo-pressure of the fracture system in the Laplace space.

Substituting the dimensionless shale gas concentration into the zone II dimensionless diffusion equation in the Laplace space leads to a matrix gas concentration of fracture network zone II. For the convenience of calculation, an adsorption/desorption index 01 and a free gas index 02 are defined to characterize the influence of the adsorbed gas and free gas in the matrix system on the gas supply to the fracture network. The adsorption/desorption index 01 and the free gas index 02 are expressed as:

01- mscqsc 2i (17) Z;ckh 2pi

02=(1-m) PLVL(18) ( pf +pL)( pi+±pL)z khT 2pi

Reducing Eq. (15) by Eqs. (17), (18) and (16) yields the matrix gas concentration in fracture network zone II:

v A(01+02) VmD +fD (19)

A pseudo-steady-state diffusion seepage differential equation of fracture network zone II is

obtained according to fracture network zone II's matrix gas concentration, gas flow mechanism

and seepage differential equation. The diffusion of the matrix block in zone II is

pseudo-steady-state diffusion, so the gas flow mechanism of fracture network zone II is:

aV q2m 2m , where q22 is the amount of gas diffused from the matrix (or the gas supply 8t amount from the matrix to the fracture), and V2. is the gas concentration of the matrix in zone

IV2 =VmD'

By substituting the zone II matrix gas concentration equation (19) and the gas flow mechanism q2. of fracture network zone II to the zone II seepage differential equation, the pseudo-steady-state diffusion seepage differential equation of fracture network zone II is obtained.

Dimensionless transform is performed on the pseudo-steady-state diffusion seepage

differential equation by the dimensionless transform relationship to obtain a differential equation

of fracture network zone II. The gas of the fracture network in zone II comes from the matrix in

zone IV and the matrix block in the fracture network of zoneII. The defined dimensionless

D,, D, rhD variables, i.e. dimensionless diffusion coefficient 'i , storage coefficient

2

inter-porosity flow coefficient and dimensionless concentration VED E are

used to nondimensionalize the pseudo-steady-state diffusion seepage differential equation,

reducing the pseudo-steady-state diffusion seepage differential equation to the differential

equation of fracture network zone II.

Laplace transform is performed on the zone II differential equation to obtain a zone II differential equation after the Laplace transform:

2 fD4mD ±2 + f2 0 2 2 2 fD 2mD (20) YD aXD XDD 2mD

In the equation, V2fD represents a dimensionless pseudo-pressure of the zone II fracture in

2rh D2 = 2hD 2m the Laplace space; G2 is the defined storage coefficient,

represents a diffusion coefficient of zone II; D 2 mD represents a dimensionless diffusion coefficient of zone II; 2mD represents a dimensionless gas concentration of the zone II matrix in the Laplace space; V2D represents a dimensionless pseudo-pressure in the zone II fracture.

Derivative finding is performed on the solution of the zone IV diffusion equation. Because the second term on the left side of Eq. (20) includes the expression of the zone IV matrix diffusion equation after derivative finding, the solution of Eq. (20) requires the derivative finding of the zone IV matrix diffusion equation. The derivative finding on the zone IV matrix diffusion equation is as follows:

8 V~ / ah ~i~ 4m 2fD s / D4mD tanh s / D4mD D(21) D XD=D

Using the solution of the zone IV diffusion equation after derivative founding to reduce the Laplace-transformed zone II differential equation leads to a reduced zone II differential equation (that is, a reduced Laplace-transformed zone II differential equation). To facilitate the solution, a second functionfp(s) is defined:

#2 2f,,q pZ, slD ta[s/D nD XD-1)]+S+ So>2 AU(012m+022.) (22) f2 (S 2 Zrckh 2p, 4 D

In the equation, #2j represents a porosity of the fracture network;km represents an 0 inter-porosity flow coefficient of the zone II matrix; 12m represents an adsorption/desorption index of the zone II matrix; 0 2 2m represents a free gas index of zone II.

Substituting the solution of the zone IV diffusion equation (Eq. (21)) after derivative finding and the second functionf2(s) (Eq. (22)) into the Laplace-transformed zone II differential equation leads to a reduced zone II differential equation:

2 V2 fD 0(23) 2 .YD

In the equation, (2f represents a dimensionless pseudo-pressure in the zone II fracture in the Laplace space.

The reduced zone II differential equation is solved to obtain a general solution of the reduced zone II differential equation. The general solution of the seepage differential equation of fracture network zone II is:

V 2 fD(YD) = A cosh[ f 2 (s)(yD YD)]+ Bsinh f 2 (s)(yD YD)] (24)

In the equation, A and B represent undetermined coefficients for the general solution of the seepage differential equation of fracture network zone II.

Derivative finding is performed on the solution of the zone III diffusion equation. The derivative finding on the solution of the zone III matrix diffusion equation is as follows:

3_2fD YDY J3 (S) tanh[ J3(S) (YeD Y1D) (25) (25) aOYD O2DfD YD

A second preset condition is obtained, where the second preset condition is yD=ylD.

Derivative finding is performed on the general solution in the second preset condition. When YD=YD, A and B are expressed by:

=__A O2fD -Bf 2 (s) YD LYD A aYD D

Substituting the expressions of A and B into Eq. (24) yields the general solution in the second preset condition:

V 2 fD (YD /2fDYDYID cosh[ 2 (S)(YD YID)] +a 2 fD sinh[ 2 (YD -YID)]/ J 2 (S)

(26)

Derivative finding is performed on Eq. (26).

According to the solution of the zone III diffusion equation after derivative finding and a relationship between the fracture gas concentration and pseudo-pressure, the general solution in the second preset condition after derivative finding is used to obtain an outer boundary condition of fracture network zone II. Considering the relationship between the fracture gas concentration and pseudo-pressure:

V2JDYD D ElID = rcTkhqS02fC 2pi, 2JD Y D=YlD (27)

the solutions of the zone III diffusion equation after derivative finding, namely Eq. (25) and Eq. (27), are substituted into the general solution in the second preset condition after derivative finding, to obtain the outer boundary condition of the zone II seepage differential equation:

aY2fD 02f2 q~c p Z aVIJD V=2flZD YD 0)3 3(s) tanh[j3(S)(YeD D] (28) ayD ID YDYD K kh 2p

In the equation, C3 is the defined storage coefficient, Dsc 3m , and D3m represents a diffusion coefficient of zone III.

To facilitate the solution, a third function z3(s) is defined:

z 2 f Zscqsc p,Z, Z3(S=CV (Zrckh tanh f 3 (s)(yD 2p, )ah[J S YD]J1D 3) S 2()(9

Reducing Eq. (28) by the third function z3(s) yields the outer boundary condition of fracture

az 3 sf 2 ()fDf YD YID aYD network zone II: YD YID (30).

An inner boundary condition of fracture network zone II is obtained. The pressure at the

interface of zone II and zone I is equal, so the inner boundary condition of zone II is:

Y2D D=WFD2 1FD D=WFD2 (31). In the equation,1FD represents a dimensionless

pseudo-pressure in the main fracture of zone I in the Laplace space.

The zone II outer boundary condition, the zone II inner boundary condition and the general solution in the second preset condition are used to obtain a solution of the zone II seepage differential equation. Substituting the zone II outer boundary condition (Eq. (30)) into the general solution (Eq. (26)) in the second preset condition leads to:

'2./D V2,D cosh(YD)((yD 1D) S) YD-YD=ll YM] -(S)(32) osh[J2 3 D sinh[ 2 (S) (YD-YD)]

Eq. (32) produces:

cosh[V2s (yD 1 (~ih(D 1 D Vf~fD (YD)I= coshDf_1D D 2()y YID)]-Z 3 (S) silh f 2(S)(YD -YID)]}

(33)

Substituting the zone II inner boundary condition (Eq. (31)) into Eq. (33) leads to:

V D -FD1 = fDL-FD1

/2JD YD D cosh[J 2 (WFD 2-YD)] -Z 3 (s)sinh[ 2 (WF D 1/2-YlD)]} (34)

Eq. (34) produces:

YD D FD /{D coshf2 (W/2D -yD)] Z 3(s)sinh f2 (WFD /2YID)] (35)

To facilitate the solution, a fourth function h2 (s) and a fifth function c 2 (s,YD) redefined:

h2 (s)= cosh f 2 (s)(WFD /2-yl)]z 3 (s)sinh f 2 (s)(WFD/2-YlD)j (36)

cosh f2 (s)(yD D) Z 3 (s)sinh[ f 2 (s)(yD YD)] c2 (SY) (37) 2'D h2 (s)

Reducing Eq. (35) by Eqs. (36) and (37) yields a solution of the zoneII seepage differential equation:

V2JD D (Y)C2 (SYD) V 1FDYDWFD/2 (38)

Step 108: Solve the zone I seepage differential equation by the dimensionless transform relationship, Laplace transform and the solution of the zone II seepage differential equation to obtain a solution of the zone I seepage differential equation.

Step 108 may specifically include:

Dimensionless transform is performed on the zone I seepage differential equation by the dimensionless transform relationship to obtain a zone I dimensionless seepage differential equation. According to the equation of state of an ideal gas and the definition of pseudo-pressure, Eq. (1) may be transformed into:

2 F 2 F __CF a /F + __ (39) ax 2 ay 2 kF at

In the equation, VF represents a pseudo-pressure in the main fracture of zone I.

According to the dimensionless definition, Eq. (39) is reduced into a zone I dimensionless seepage differential equation:

a2VFD a2 VFD __ C9FD +.D 1FD (40) qX D 2 2 9YD 77FD &tD

In the equation, VFD represents a dimensionless pseudo-pressure of the zone I fracture, and 7 7FD represents a dimensionless conductivity factor of the zone I fracture.

A calculus operation is performed on the zone I dimensionless seepage differential equation, and a zone I dimensionless seepage differential equation is obtained after the calculus operation. Along the seepage direction of the fracture network, integration is performed from the interface of the fracture network and the hydraulically created fracture to the half width of the hydraulically created fracture. By performing the calculus operation on the zone I dimensionless seepage differential equation (Eq. (40)), Eq. (40) is reduced into a zone I dimensionless seepage differential equation after the calculus operation:

2 FD 2 0 fFD 1 ' fFD 7 2 FD D FD 7FD D 2

A continuous relationship of a gas flow flux at the interface between fracture network zone II and the hydraulically created fracture is obtained. It is considered that the gas flow flux at the interface between the fracture network and the hydraulic fracture is continuous, that is, there is a continuous relationship of the gas flow flux:

kFhaVf = kfh (42). ay y=W/ 2 ay y=W1 2

In the equation, /F represents a pseudo-pressure of the hydraulically created fracture, and Vf represents a pseudo-pressure of the zone II fracture.

According to the continuous relationship of the gas flow flux and the dimensionless transform relationship, dimensionless transform is performed on the zone I dimensionless seepage differential equation after calculus operation, to obtain a reduced zone I seepage differential equation. According to the continuous relationship of the gas flow flux (Eq. 42) and the dimensionless transform relationship, dimensionless transform is performed on the zone I dimensionless seepage differential equation (Eq. 41) after calculus operation, to obtain a reduced zone I seepage differential equation:

8 VFD 2 'D//FD VfFD DXD 2 FCD 8 YD 8 _ FD 7FD tD D2 (43)

Laplace transform is performed on the reduced zone I seepage differential equation, and a zone I seepage differential equation is obtained after the Laplace transform. Laplace transform is performed on the seepage differential equation (Eq. 43) of the hydraulically created fracture based on the dimensionless time tD, and a zone Iseepage differential equation is obtained after the Laplace transform:

a2 /FD 2 aytf __S

FCD OYD _W 7F ~~ FD OXD2 YD 2 (44)

In the equation, yfFD represents a dimensionless pseudo-pressure of the hydraulically created fracture in the Laplace space.

Derivative finding is performed on the solution of the zoneII seepage differential equation. The solution of the zone II seepage differential equation after derivative finding is:

__2_D Ffs)sinh[j7(s(WFD 12-yD3)]-z3(s)jj cosh[j(s)(WD 12-yD)] alWyD D/ =%2l2h 2(s) C9YD"_F/2 k()(45)

To facilitate the solution, a six function F2(s) is defined:

F() =/2(s)sinh7f2(s)(WFD 12-YlD)J-z3 (S2'7(cSCoshf(s(WFD 1YD)J (46) h2(s)

Reducing Eq. (45) by the sixth function F2(s) (Eq. (46)) leads to:

1FD F2(s) (47) YWD2 aD YD=WFD12 YD=WFD/

A boundary condition of hydraulically fractured zone I is obtained. The zone I boundary condition is defined by the characteristics of hydraulically fractured zone I. The outer boundary of zone I is a distal end of the fracture, which is assumed to be a no-flow boundary; the inner boundary of zone I is the interface between the hydraulically created fracture and the wellbore. aV I FD 0

The zone I boundary condition includes an outer boundary condition: aXD XD andan VFD _

8 XD sFCD F inner boundary condition: XD CD where CD s a dimensionles conductivity.

According to the solution of the zone II seepage differential equation after derivative finding, the zone I boundary condition is used to solve the zone I seepage differential equation after Laplace transform, to obtain a solution of the zone I seepage differential equation. The solution (Eq. (47)) of the zone II seepage differential equation after derivative finding is substituted into the zone I seepage differential equation (Eq. (44)) after Laplace transform. In order to facilitate the solution, a seventh function g2(s) is defined:

s 2 92(S) - F 2 (s) 771FD FCD (48)

In the equation, rFD represents a dimensionless conductivity factor of the hydraulically created fracture.

Substituting Eq. (47) into Eq. (44) and reducing Eq. (44) by the defined seventh function g2(s) lead to: ax 62k1Y1FD

9D (49)

Solving Eq. (49) with the zone I boundary condition yields a solution of the zone I seepage differentialequation:

VIFT- cosh g 2 (s)(1-xD)] 1FD sFCD g2 (s) sinh(Vg2 (s)) (50)

Step 109: Obtain a first preset condition, where the first preset condition is XD=O.

Step 110: Use the solution of the zone I seepage differential equation to obtain a bottom hole pseudo-pressure solution according to the first preset condition.

Step 110 may specifically include:

Substituting the first preset condition into the solution of the zone I seepage differential equation leads to a bottom hole pseudo-pressure solution. The first preset condition XD=0 is substituted into the solution (Eq. (50)) of the zone I seepage differential equation to obtain a bottom hole pseudo-pressure solution. When XD=0, the bottom hole pseudo-presure solution is:

VwD _[ 1fFD '1D=O sFD g2(s)tanh g2(s) (51)

In the equation, Vf1 vD represents a dimensionless bottom hole pseudo-pressure in the Laplace space.

Step 111: Obtain a dimensionless production solution by Duhamel's Principle according to the bottom hole pseudo-pressure solution. According to Duhamel's Principle, the bottom hole pseudo-pressure which considers the skin effect and reservoir effect in the Laplace space is:

sywD +Sc PwD (52) s+CDS2VfwD + Sc)

In the equation, pvD represents a dimensionless bottom hole pseudo-pressure considering the skin effect and the storage effect in the Laplace space; Sc is a skin coefficient, and CD is a dimensionless storage coefficient.

According to the research resultsError! Reference source not found.of Van Everdingen and Hurst, in the Laplace space, a dimensionless bottom hole pseudo-pressure under constant production conditions and a dimensionless production under constant pressure conditions have the following conversion relationship:

_ 1 qiwD 2 (53) S pwD

In the equation, qw represents a dimensionless production considering the skin effect and the storage effect in the Laplace space under constant pressure conditions.

Substituting Eq. (51) into Eq. (52) yields the dimensionless bottom hole pseudo-pressure considering the skin effect and the storage effect under constant production conditions, and substituting the conversion relationship (Eq. (53)) and the dimensionless bottom hole pseudo-pressure (Eq. (52)) yields the dimensionless production solution in the Laplace space.

Step 112: Predict the production of the fractured horizontal well in the shale gas reservoir by using a Stehfest numerical inversion method according to the dimensionless production solution.

Step 112 may specifically include: The Stehfest numerical inversion method is used to numerically solve the dimensionless production solution to obtain the predicted production of the fractured horizontal well in the shale gas reservoir.

This embodiment further provides a system for predicting the production of a fractured horizontal well in a shale gas reservoir. The system for predicting the production of a fractured horizontal well in a shale gas reservoir includes: a seepage zone dividing module, a seepage zone equation obtaining module, a dimensionless transform relationship obtaining module, a zone V diffusion equation solving module, a zone IV diffusion equation solving module, a zone III diffusion equation solving module, a zone II seepage differential equation solving module, a zone I seepage differential equation solving module, a first preset condition obtaining module, a bottom hole pseudo-pressure solution obtaining module, a dimensionless production solution obtaining module and a production predicting module.

The seepage zone dividing module is configured to divide a fractured horizontal well to be predicted in a shale gas reservoir into five seepage zones according to a matrix block and a fracture network after hydraulic fracturing, where the five seepage zones include: hydraulically fractured zone I, fracture network zone II, pure matrix zone III, pure matrix zone IV and pure matrix zone V.

The seepage zone equation obtaining module is configured to obtain a zone I seepage differential equation for hydraulically created fracture zone I, a zone II seepage differential equation and a zone II diffusion equation for fracture network zone II, a zone III diffusion equation for pure matrix zone III, a zone IV diffusion equation for pure matrix zone IV and a zone V diffusion equation for pure matrix zone V.

The dimensionless transform relationship obtaining module is configured to obtain a preset dimensionless transform relationship.

The zone V diffusion equation solving module is configured to solve the zone V diffusion equation by the dimensionless transform relationship and Laplace transform to obtain a solution of the zone V diffusion equation.

The zone V diffusion equation solving module may specifically include a zone V dimensionless transform unit, a zone V Laplace transform unit, a zone V boundary condition obtaining unit and a zone V dimensionless diffusion equation solving unit. The zone V dimensionless transform unit is configured to perform dimensionless transform on the zone V diffusion equation by the dimensionless transform relationship to obtain a zone V dimensionless diffusion equation.

The zone V Laplace transform unit is configured to perform Laplace transform on the zone V dimensionless diffusion equation to obtain a zone V dimensionless diffusion equation in a Laplace space.

The zone V boundary condition obtaining unit is configured to obtain a boundary condition of pure matrix zone V.

The zone V dimensionless diffusion equation solving unit is configured to use the zone V boundary condition to solve the zone V dimensionless diffusion equation in the Laplace space to obtain a solution of the zone V diffusion equation.

The zone IV diffusion equation solving module is configured to solve the zone IV diffusion equation by the dimensionless transform relationship and Laplace transform to obtain a solution

of the zone IV diffusion equation.

The zone IV diffusion equation solving module may specifically include a zone IV dimensionless transform unit, a zone IV Laplace transform unit, a zone IV boundary condition obtaining unit and a zone IV dimensionless diffusion equation solving unit. The zone IV dimensionless transform unit is configured to perform dimensionless transform on the zone IV diffusion equation by the dimensionless transform relationship to obtain a zone IV dimensionless diffusion equation.

The zone IV Laplace transform unit is configured to perform Laplace transform on the zone IV dimensionless diffusion equation to obtain a zone IV dimensionless diffusion equation in the Laplace space.

The zone IV boundary condition obtaining unit is configured to obtain a boundary condition of pure matrix zone IV.

The zone IV dimensionless diffusion equation solving unit is configured to use the zone IV boundary condition to solve the zone IV dimensionless diffusion equation in the Laplace space to obtain a solution of the zone IV diffusion equation.

The zone III diffusion equation solving module is configured to solve the zone III diffusion equation by the dimensionless transform relationship, Laplace transform and the solution of the zone V diffusion equation to obtain a solution of the zone III diffusion equation.

The zone III diffusion equation solving module may specifically include a zone III dimensionless transform unit, a zone III Laplace transform unit, a zone V diffusion equation solution derivative finding unit, a zone III boundary condition obtaining unit and a zone III dimensionless diffusion equation solving unit. The zone III dimensionless transform unit is configured to perform dimensionless transform on the zone III diffusion equation by the dimensionless transform relationship to obtain a zone III dimensionless diffusion equation.

The zone III Laplace transform unit is configured to perform Laplace transform on the zone III dimensionless diffusion equation to obtain a zone III dimensionless diffusion equation in the Laplace space.

The zone V diffusion equation solution derivative finding unit is configured to perform derivative finding on the solution of the zone V diffusion equation.

The zone III boundary condition obtaining unit is configured to obtain a boundary condition of pure matrix zone III.

The zone III dimensionless diffusion equation solving unit is configured to use the zone III boundary condition to solve the zone III dimensionless diffusion equation in the Laplace space according to the solution of the zone V diffusion equation after derivative finding, to obtain a solution of the zone III diffusion equation.

The zone II seepage differential equation solving module is configured to solve the zone II seepage differential equation by the dimensionless transform relationship, Laplace transform, the zone II diffusion equation, the solution of the zone IV diffusion equation and the solution of the zone III diffusion equation, to obtain a solution of the zoneII seepage differential equation.

The zone II diffusion equation solving module may specifically include a zone II dimensionless transform unit, a zone II first Laplace transform unit, a dimensionless shale gas concentration solving unit, a zone II matrix gas concentration solving unit, a pseudo-steady-state diffusion seepage differential equation obtaining unit, a zone II differential equation obtaining unit, a zone II second Laplace transform unit, a zone IV diffusion equation solution derivative finding unit, a zone II differential equation reducing unit, a zone II differential equation general solution obtaining unit, a zone III diffusion equation solution derivative finding unit, a second preset condition obtaining unit, a second preset condition general solution derivative finding unit, a zone II outer boundary condition obtaining unit, a zoneII inner boundary condition obtaining unit and a zone II seepage differential equation solving unit. The zoneII dimensionless transform unit is configured to perform dimensionless transform on the zone II diffusion equation by the dimensionless transform relationship to obtain a zone II dimensionless diffusion equation.

The zone II first Laplace transform unit is configured to perform Laplace transform on the zone II dimensionless diffusion equation to obtain a zone II dimensionless diffusion equation in the Laplace space.

The dimensionless shale gas concentration solving unit is configured to perform dimensionless transform on a pressure function of fracture network zone II by the dimensionless transform relationship to obtain a dimensionless shale gas concentration when the gas supply from the matrix block to fracture network zone II reaches equilibrium.

The zone II matrix gas concentration solving unit is configured to substitute the dimensionless shale gas concentration into the zone II dimensionless diffusion equation in the Laplace space to obtain a matrix gas concentration of fracture network zone II.

The pseudo-steady-state diffusion seepage differential equation obtaining unit is configured to obtain a pseudo-steady-state diffusion seepage differential equation of fracture network zone II according to fracture network zone II's matrix gas concentration, gas flow mechanism and seepage differential equation.

The zone II differential equation obtaining unit is configured to perform dimensionless transform on the pseudo-steady-state diffusion seepage differential equation by the dimensionless transform relationship to obtain a differential equation of fracture network zone II.

The zone II second Laplace transform unit is configured to perform Laplace transform on the zone II differential equation to obtain a zone II differential equation after the Laplace transform.

The zone IV diffusion equation solution derivative finding unit is configured to perform derivative finding on the solution of the zone IV diffusion equation.

The zone II differential equation reducing unit is configured to use the solution of the zone IV diffusion equation after derivative founding to reduce the Laplace-transformed zone II differential equation to obtain a reduced zone II differential equation.

The zone II differential equation general solution obtaining unit is configured to solve the reduced zone II differential equation to obtain a general solution of the reduced zone II differential equation.

The zone III diffusion equation solution derivative finding unit is configured to perform derivative finding on the solution of the zone III diffusion equation.

The second preset condition obtaining unit is configured to obtain a second preset condition.

The second preset condition general solution derivative finding unit is configured to perform derivative finding on a general solution in the second preset condition.

The zone II outer boundary condition obtaining unit is configured to use the general solution in the second preset condition after the derivative finding to obtain an outer boundary condition of fracture network zone II, according to the solution of the zone III diffusion equation after derivative finding and a relationship between a fracture gas concentration and a pseudo-pressure.

The zone II inner boundary condition obtaining unit is configured to obtain an inner boundary condition of fracture network zone II.

The zone II seepage differential equation solving unit is configured to use the zone II outer boundary condition, the zone II inner boundary condition and the general solution in the second preset condition to obtain a solution of the zoneII seepage differential equation.

The zone I seepage differential equation solving module is configured to solve the zone I seepage differential equation by the dimensionless transform relationship, Laplace transform and the solution of the zone II seepage differential equation to obtain a solution of the zone I seepage differential equation.

The zone I diffusion equation solving module may specifically include a zone I dimensionless transform unit, a calculus operation performing unit, a gas flow flux continuous relationship obtaining unit, a zone I dimensionless seepage differential equation reducing unit, a zone I Laplace transform unit, a zone II seepage differential equation solution derivative finding unit, a zone I boundary condition obtaining unit and a zone I seepage differential equation solving unit. The zone I dimensionless transform unit is configured to perform dimensionless transform on the zone I seepage differential equation by the dimensionless transform relationship to obtain a zone I dimensionless seepage differential equation.

The calculus operation performing unit is configured to perform a calculus operation on the zone I dimensionless seepage differential equation to obtain a zone I dimensionless seepage differential equation after the calculus operation.

The gas flow flux continuous relationship obtaining unit is configured to obtain a continuous relationship of a gas flow flux at an interface between fracture network zone II and the hydraulically created fracture.

The zone I dimensionless seepage differential equation reducing unit is configured to perform dimensionless transform on the zone I dimensionless seepage differential equation after calculus operation according to the continuous relationship of the gas flow flux and the dimensionless transform relationship, to obtain a reduced zone I seepage differential equation.

The zone I Laplace transform unit is configured to perform Laplace transform on the reduced zone I seepage differential equation to obtain a zone I seepage differential equation after the Laplace transform.

The zone II seepage differential equation solution derivative finding unit is configured to perform derivative finding on the solution of the zoneII seepage differential equation.

The zone I boundary condition obtaining unit is configured to obtain a boundary condition of hydraulically fractured zone I.

The zone I seepage differential equation solving unit is configured to use the zone I boundary condition to solve the zone I seepage differential equation after the Laplace transform according to the solution of the zone II seepage differential equation after derivative finding, to obtain a solution of the zone I seepage differential equation.

The first preset condition obtaining module is configured to obtain a first preset condition.

The bottom hole pseudo-pressure solution obtaining module is configured to use the solution of the zone I seepage differential equation to obtain a bottom hole pseudo-pressure solution according to the first preset condition.

The bottom hole pseudo-pressure solution obtaining module may specifically include a bottom hole pseudo-pressure solution obtaining unit, configured to substitute the first preset condition into the solution of the zone I seepage differential equation to obtain a bottom hole pseudo-pressure solution.

The dimensionless production solution obtaining module is configured to obtain a dimensionless production solution by Duhamel's Principle according to the bottom hole pseudo-pressure solution.

The production predicting module is configured to predict the production of the fractured horizontal well in the shale gas reservoir by using a Stehfest numerical inversion method according to the dimensionless production solution.

The production predicting module may specifically include a production predicting unit, configured to use the Stehfest numerical inversion method to numerically solve the dimensionless production solution to obtain the predicted production of the fractured horizontal well in the shale gas reservoir.

The method and system for predicting the production of a fractured horizontal well in a shale gas reservoir in the present disclosure are used to carry out sensitivity analysis of production decline. The present disclosure takes the physical parameters and engineering parameters of the gas reservoir as input, sets the time required for prediction, obtains the predicted production, and plots the predicted production and the corresponding dimensionless time on the coordinates to obtain a typical curve of production decline in real space. The physical parameters and engineering parameters of the gas reservoir include reservoir parameters, gas reservoir temperature, original formation pressure, gas density, gas viscosity and gas compressibility. The reservoir parameters are shown in Table 1. By changing the parametric values, the sensitivity of factors affecting production decline is analyzed.

Table 1 Reservoir parameters

Reservoir parameter Value Reservoir parameter Value

Zone II matrix diffusion coefficient D2m Zone II size xexyexh (m) 50x30x70 5x 10-5 2 (m /s)

Zone III matrix diffusion coefficient D4m Model size xexyexh (m) 100x100x70 5x 10-6 (m2 /s)

Width of hydraulically created fracture WF Zone IV matrix diffusion coefficient D 4 m 0.005 5x 10-6 (m) (m2/s)

Half length of hydraulically created Zone V matrix diffusion coefficient Dm 50 5x10-6 fracture (m) (m2/s)

Permeability of hydraulically created Matrix porosity(%) 3 500 fracture (mD)

SRV zone fracture network permeability 0.1 Isotherm adsorption volume V (sm/mlt) 3 (mD)

SRV zone matrix diffusion coefficient D2m 5x 10-5 Isothermal adsorption pressure pL (MPa) 4 (m2/s)

Dimensionless wellbore storage 0.001 Skin coefficient Sc 0.1 coefficient Sc

The parameters in Table 1 are substituted into the method and system for predicting the production of the fractured horizontal well in the shale gas reservoir. Only by changing the value of the width of the SRV zone, the dimensionless production solution corresponding to the width of the SRV zone is obtained. They are provided in FIGS. 4-9. The legends of FIGS. 4-9 show qD

and the derivative of qD at different yl. As shown in FIG. 4, the width of the SRV zone does not affect the overall shape of the production decline curve, and it mainly affects the bilinear flow of the fracture and the matrix as well as the linear flow of the matrix. A larger value of the width of the SRV zone leads to farther propagation of the fracture, a longer duration of the bilinear flow of the fracture and the matrix, and a longer time for the gas well to maintain high production. The vertical axis in FIG. 4 shows the dimensionless production and the derivative thereof, where qD represents the dimensionless production, and qD represents the derivative of the dimensionless production (that is, the rate of change of production); yi represents the width of the SRV zone.

As shown in FIG. 5, the wellbore storage coefficient mainly affects the initial stage of production. It is not easy to observe in actual production due to the short duration. The wellbore storage effect hardly affects the production stage following the linear flow of the fracture. In FIG. , CD represents the dimensionless storage coefficient, and the horizontal axis represents the dimensionless time tD.

The skin coefficient Sc represents the degree of pollution near the wellbore. A larger skin coefficient indicates more serious pollution and a larger flow resistance, which means a lower production of the gas well under the same pressure difference. As shown in FIG. 6, the skin coefficient does not affect the gas flow during the storage effect stage of the wellbore. A larger skin coefficient leads to a greater depression of the dimensionless production derivative curve. The skin coefficient has a more obvious effect on the linear flow of the fracture in the early stage. After the bilinear flow of the fracture and the matrix is formed, this effect gradually decreases until it almost disappears.

The inter-porosity flow coefficient reflects the difference in the physical properties of the matrix and the fracture. As shown in FIG. 7, a larger inter-porosity flow coefficient k indicates a greater difference, which makes it easier for the fluid in the matrix system to flow to the fracture network, and leads to an earlier start of the bilinear flow, a longer duration and a higher production.

The fracture conductivity reflects the seepage capacity of the hydraulically created fracture. As shown in FIG. 8, a greater fracture permeability leads to a lower flow resistance for the fluid in the fracture, and indicates a stronger fracture conductivity, an earlier start of the linear flow of the fracture, and a higher production of the gas well. The fracture conductivity hardly affects the duration of each production decline stage. In FIG. 8, FCD represents the fracture conductivity factor.

For the dual-medium SRV zone, the contribution of adsorbed gas and free gas in the matrix system to the fracture network is considered, and the adsorption/desorption index 01 and free gas index 02 are introduced to characterize the influence of the adsorbed gas and free gas in the matrix system on the gas supply to the fracture network. As shown in FIG. 9, when the gas supply of the matrix system to the fracture network is ignored, that is, when the adsorption/desorption index 01 and the free gas index 02 are both zero, the duration of the bilinear flow stage of the hydraulically created fracture and the fracture network in the treatment zone is very short. The Langmuir model is a pressure-based desorption model. When the pressure decreases, the adsorbed gas desorbs rapidly, so the contribution of the adsorbed gas to the fracture network in the matrix system is greater than that of the free gas.

The sensitivity analysis of production decline shows that: a greater width of the SRV zone leads to a longer time for the gas well to maintain a high production; a greater storage coefficient of the wellbore leads to a higher production in the initial stage; the skin coefficient has a more obvious effect on the linear flow in the early fracture; a larger inter-porosity flow coefficient indicates an earlier start of the bilinear flow, a longer duration and a higher production; a stronger fracture conductivity indicates an earlier start of the linear flow of the fracture and a higher production of the gas well; the contribution of the adsorbed gas to the fracture network in the matrix system is greater than that of free gas.

The present disclosure divides the fractured horizontal well in the shale gas reservoir into a fracture network zone (SRV zone) with a fracture-pore medium and pure matrix zones with a single porous medium. Considering the different diffusion modes of the matrix in different zones, the present disclosure uses Fick's First Law to describe the pseudo-steady-state diffusion of the matrix in the fracture network zone, Fick's Second Law to describe the unsteady-state diffusion of the matrix in the pure matrix zone, and Darcy's Law to describe the seepage in the fracture network. The present disclosure establishes a model for predicting the production of the fractured horizontal well in the shale gas reservoir under the conditions of matrix-microfracture coupling and hydraulically created fracture-microfracture coupling (as shown in steps 102 to 112). The present disclosure more accurately describes the actual flow law of the shale gas reservoir, provides a new method for production decline analysis and prediction of the shale gas well, and provides theoretical basis and production suggestions for the beneficial development of the shale gas reservoir. The matrix-microfracture coupling is reflected in the gas concentration difference between the matrix and the fracture, and the matrix supplies gas to the fracture by diffusion.

Each embodiment of the present specification is described in a progressive manner, each embodiment focuses on the difference from other embodiments, and the same and similar parts between the embodiments may refer to each other. For a system disclosed in the embodiments, since the system corresponds to the method disclosed in the embodiments, the description is relatively simple, and reference can be made to the method description.

In this specification, several specific embodiments are used for illustration of the principles and implementations of the present disclosure. The description of the foregoing embodiments is used to help illustrate the method of the present disclosure and the core ideas thereof. In addition, those of ordinary skill in the art can make various modifications in terms of specific implementations and scope of application in accordance with the ideas of the present disclosure. In conclusion, the content of this specification should not be construed as a limitation to the present disclosure.

Claims

THE CLAIMS DEFINING THE INVENTION ARE AS FOLLOWS:

1. A method for predicting the production of a fractured horizontal well in a shale gas reservoir, comprising:

dividing a fractured horizontal well to be predicted in a shale gas reservoir into five seepage zones according to a matrix block and a fracture network after hydraulic fracturing, wherein the five seepage zones comprise: hydraulically fractured zone I, fracture network zone II, pure matrix zone III, pure matrix zone IV and pure matrix zone V;

obtaining a preset dimensionless transform relationship;

obtaining a first preset condition;

using the solution of the zone I seepage differential equation to obtain a bottom hole pseudo-pressure solution according to the first preset condition; obtaining a dimensionless production solution by Duhamel's Principle according to the bottom hole pseudo-pressure solution; and predicting the production of the fractured horizontal well in the shale gas reservoir by using a Stehfest numerical inversion method according to the dimensionless production solution.

2. The method for predicting the production of a fractured horizontal well in a shale gas reservoir according to claim 1, wherein the solving the zone V diffusion equation by the dimensionless transform relationship and Laplace transform to obtain a solution of the zone V diffusion equation specifically comprises:

performing dimensionless transform on the zone V diffusion equation by the dimensionless transform relationship to obtain a zone V dimensionless diffusion equation;

performing Laplace transform on the zone V dimensionless diffusion equation to obtain a zone V dimensionless diffusion equation in a Laplace space;

obtaining a boundary condition of pure matrix zone V; and

using the zone V boundary condition to solve the zone V dimensionless diffusion equation in the Laplace space to obtain a solution of the zone V diffusion equation.

3. The method for predicting the production of a fractured horizontal well in a shale gas reservoir according to claim 1, wherein the solving the zone IV diffusion equation by the dimensionless transform relationship and Laplace transform to obtain a solution of the zone IV diffusion equation specifically comprises:

performing dimensionless transform on the zone IV diffusion equation by the dimensionless transform relationship to obtain a zone IV dimensionless diffusion equation;

performing Laplace transform on the zone IV dimensionless diffusion equation to obtain a zone IV dimensionless diffusion equation in the Laplace space;

obtaining a boundary condition of pure matrix zone IV; and

using the zone IV boundary condition to solve the zone IV dimensionless diffusion equation in the Laplace space to obtain a solution of the zone IV diffusion equation;

wherein the solving the zone III diffusion equation by the dimensionless transform relationship, Laplace transform and the solution of the zone V diffusion equation to obtain a solution of the zone III diffusion equation specifically comprises: performing dimensionless transform on the zone III diffusion equation by the dimensionless transform relationship to obtain a zone III dimensionless diffusion equation; performing Laplace transform on the zone III dimensionless diffusion equation to obtain a zone III dimensionless diffusion equation in the Laplace space; performing derivative finding on the solution of the zone V diffusion equation; obtaining a boundary condition of pure matrix zone III; and using the zone III boundary condition to solve the zone III dimensionless diffusion equation in the Laplace space according to the solution of the zone V diffusion equation after derivative finding, to obtain a solution of the zone III diffusion equation; wherein the solving the zone II seepage differential equation by the dimensionless transform relationship, Laplace transform, the zone II diffusion equation, the solution of the zone IV diffusion equation and the solution of the zone III diffusion equation to obtain a solution of the zone II seepage differential equation specifically comprises: performing dimensionless transform on the zone II diffusion equation by the dimensionless transform relationship to obtain a zone II dimensionless diffusion equation; performing Laplace transform on the zone II dimensionless diffusion equation to obtain a zone II dimensionless diffusion equation in the Laplace space; performing dimensionless transform on a pressure function of fracture network zone II by the dimensionless transform relationship to obtain a dimensionless shale gas concentration when the gas supply from the matrix block to fracture network zone II reaches equilibrium; substituting the dimensionless shale gas concentration into the zone II dimensionless diffusion equation in the Laplace space to obtain a matrix gas concentration of fracture network zone II; obtaining a pseudo-steady-state diffusion seepage differential equation of fracture network zone II according to fracture network zone II's matrix gas concentration, gas flow mechanism and seepage differential equation; performing dimensionless transform on the pseudo-steady-state diffusion seepage differential equation by the dimensionless transform relationship to obtain a differential equation of fracture network zone II; performing Laplace transform on the zone II differential equation to obtain a zone II differential equation after the Laplace transform; performing derivative finding on the solution of the zone IV diffusion equation; using the solution of the zone IV diffusion equation after derivative founding to reduce the Laplace-transformed zone II differential equation to obtain a reduced zone II differential equation; solving the reduced zone II differential equation to obtain a general solution of the reduced zone II differential equation; performing derivative finding on the solution of the zone III diffusion equation; obtaining a second preset condition; performing derivative finding on the general solution in the second preset condition; using the general solution in the second preset condition after the derivative finding to obtain an outer boundary condition of fracture network zone II, according to the solution of the zone III diffusion equation after derivative finding and a relationship between a fracture gas concentration and a pseudo-pressure; obtaining an inner boundary condition of fracture network zone II; and using the zone II outer boundary condition, the zone II inner boundary condition and the general solution in the second preset condition to obtain a solution of the zone II seepage differential equation; wherein the solving the zone I seepage differential equation by the dimensionless transform relationship, Laplace transform and the solution of the zone II seepage differential equation to obtain a solution of the zone I seepage differential equation specifically comprises: performing dimensionless transform on the zone I seepage differential equation by the dimensionless transform relationship to obtain a zone I dimensionless seepage differential equation; performing a calculus operation on the zone I dimensionless seepage differential equation to obtain a zone I dimensionless seepage differential equation after the calculus operation; obtaining a continuous relationship of a gas flow flux at an interface between fracture network zone II and the hydraulically created fracture; performing dimensionless transform on the zone I dimensionless seepage differential equation after calculus operation according to the continuous relationship of the gas flow flux and the dimensionless transform relationship, to obtain a reduced zone I seepage differential equation; performing Laplace transform on the reduced zone I seepage differential equation to obtain a zone I seepage differential equation after the Laplace transform; performing derivative finding on the solution of the zoneII seepage differential equation; obtaining a boundary condition of hydraulically fractured zone I; and using the zone I boundary condition to solve the zone I seepage differential equation after the Laplace transform according to the solution of the zone II seepage differential equation after derivative finding, to obtain a solution of the zone I seepage differential equation; wherein the using the solution of the zone I seepage differential equation to obtain a bottom hole pseudo-pressure solution according to the first preset condition specifically comprises: substituting the first preset condition into the solution of the zone I seepage differential equation to obtain a bottom hole pseudo-pressure solution.

4. A system for predicting the production of a fractured horizontal well in a shale gas reservoir, comprising:

a seepage zone dividing module, configured to divide a fractured horizontal well to be predicted in a shale gas reservoir into five seepage zones according to a matrix block and a fracture network after hydraulic fracturing, wherein the five seepage zones comprise: hydraulically fractured zone I, fracture network zone II, pure matrix zone III, pure matrix zone IV and pure matrix zone V;

a dimensionless transform relationship obtaining module, configured to obtain a preset dimensionless transform relationship; a zone V diffusion equation solving module, configured to solve the zone V diffusion equation by the dimensionless transform relationship and Laplace transform to obtain a solution of the zone V diffusion equation; a zone IV diffusion equation solving module, configured to solve the zone IV diffusion equation by the dimensionless transform relationship and Laplace transform to obtain a solution of the zone IV diffusion equation; a zone III diffusion equation solving module, configured to solve the zone III diffusion equation by the dimensionless transform relationship, Laplace transform and the solution of the zone V diffusion equation to obtain a solution of the zone III diffusion equation; a zone II seepage differential equation solving module, configured to solve the zone II seepage differential equation by the dimensionless transform relationship, Laplace transform, the zone II diffusion equation, the solution of the zone IV diffusion equation and the solution of the zone III diffusion equation, to obtain a solution of the zoneII seepage differential equation; a zone I seepage differential equation solving module, configured to solve the zone I seepage differential equation by the dimensionless transform relationship, Laplace transform and the solution of the zone II seepage differential equation to obtain a solution of the zone I seepage differential equation; a first preset condition obtaining module, configured to obtain a first preset condition; a bottom hole pseudo-pressure solution obtaining module, configured to use the solution of the zone I seepage differential equation to obtain a bottom hole pseudo-pressure solution according to the first preset condition; a dimensionless production solution obtaining module, configured to obtain a dimensionless production solution by Duhamel's Principle according to the bottom hole pseudo-pressure solution; and a production predicting module, configured to predict the production of the fractured horizontal well in the shale gas reservoir by using a Stehfest numerical inversion method according to the dimensionless production solution.

5. The system for predicting the production of a fractured horizontal well in a shale gas reservoir according to claim 4, wherein the zone V diffusion equation solving module specifically comprises:

a zone V dimensionless transform unit, configured to perform dimensionless transform on the zone V diffusion equation by the dimensionless transform relationship to obtain a zone V dimensionless diffusion equation;

a zone V Laplace transform unit, configured to perform Laplace transform on the zone V dimensionless diffusion equation to obtain a zone V dimensionless diffusion equation in a Laplace space;

a zone V boundary condition obtaining unit, configured to obtain a boundary condition of pure matrix zone V; and

a zone V dimensionless diffusion equation solving unit, configured to use the zone V boundary condition to solve the zone V dimensionless diffusion equation in the Laplace space to obtain a solution of the zone V diffusion equation;

wherein the zone IV diffusion equation solving module specifically comprises:

a zone IV dimensionless transform unit, configured to perform dimensionless transform on the zone IV diffusion equation by the dimensionless transform relationship to obtain a zone IV dimensionless diffusion equation;

a zone IV Laplace transform unit, configured to perform Laplace transform on the zone IV dimensionless diffusion equation to obtain a zone IV dimensionless diffusion equation in the Laplace space;

a zone IV boundary condition obtaining unit, configured to obtain a boundary condition of pure matrix zone IV; and

a zone IV dimensionless diffusion equation solving unit, configured to use the zone IV boundary condition to solve the zone IV dimensionless diffusion equation in the Laplace space to obtain a solution of the zone IV diffusion equation.