Full wellbore pressure calculation method, device and computer readable storage medium for fracturing horizontal well

Abstract

The invention discloses a full wellbore pressure calculation method, device and computer-readable storage medium for fracturing horizontal wells, includes steps: collect basic parameters, establish fully implicit numerical model of formation flow; solve the full implicit numerical model under an inner boundary condition of a constant gas production rate; calculate wellbore pressure change of each fracturing section of horizontal well; calculate the bottom-hole pressure of each fracture initiation point; repeat steps until the variables converge, the bottom-hole flow pressure and the bottom-hole pressure at each fracture initiation point are obtained at this time step; the wellbore pressure and wellhead casing pressure are calculated at this time step. The method is for predicting the variation law of full wellbore pressure and wellhead casing pressure by using production data in the absence of bottom-hole pressure test data, and has practical value for accurate prediction of production performance of gas reservoir fracturing horizontal wells.

Claims

1. A method for calculating a full wellbore pressure of a fracturing horizontal well, wherein the method comprises the following steps: S1, collect the wellbore structural parameters, fluid physical properties, reservoir parameters and historical production data at different time steps of the gas well, and establish a formation flow fully implicit numerical model based on an embedded discrete fracture model; S2, calculate the formation flow fully implicit numerical model based on the embedded discrete fracture model under an inner boundary condition of a constant gas production rate at a certain time step, to obtain a gas production, liquid production, grid block pressure and initial bottom-hole flow pressure of each fracture; S3, establish a pressure drop model considering the gravity loss, frictional resistance loss and fracture convergence loss by combining the gas production, liquid production, and the pressure of the grid block where the fracture initiation point is located, to obtain a wellbore pressure variation of each fracturing section of the horizontal well; S4, calculate the bottom-hole pressure of each fracture initiation point from the initial bottom-hole flow pressure and the wellbore pressure change of each fracturing section; S5, perform the formation flow fully implicit numerical model based on the embedded discrete fracture model in step S2 by combining the bottom-hole flow pressure and the bottom-hole pressure of each fracture initiation point, and perform iterative cycles of the steps S2-S4 until the variables converge, to obtain the bottom-hole flow pressure and bottom-hole pressure at each fracture initiation point under this time step; S6, bring the bottom-hole flow pressure and the total gas production and liquid production at the time step into the pressure drop model considering gravity loss and frictional resistance loss in the inclined section and vertical well section, to obtain wellbore pressure and wellhead casing pressure of each section through the iterative cycle at each section; S7, for each time step, repeat steps S2-S6 to obtain the full wellbore pressure at different production times.

2. The full wellbore pressure calculation method of fracturing horizontal well according to claim 1, wherein establishment of the formation flow fully implicit numerical model based on the embedded discrete fracture model described in the step S1 comprises the following steps: S11, establish a formation flow mathematical model based on the embedded discrete fracture model, in which the dual medium model is used to describe the flow exchange between the matrix and micro-fractures, and the embedded discrete fracture model is used to describe the flow exchange between the matrix and artificial fractures, based on mass conservation equation and Darcy's law to establish seepage differential equation, the formation flow mathematical model of the embedded discrete fracture model is as follows: matrix system seepage equation: $\begin{array}{cc}\nabla .Math.\left(\frac{k_m{k}_{\mathrm{rlm}}}{{\mu }_{\mathrm{lm}}{B}_{\mathrm{lm}}}\nabla p_{\mathrm{lm}}\right)-q_{l,\mathrm{fm}}=\frac{\partial }{\partial t}\left(\frac{{\varphi }_m{S}_{\mathrm{lm}}}{B_{\mathrm{lm}}}\right)& \left(1\right)\end{array}$ micro-fracture system seepage equation: $\begin{array}{cc}\nabla .Math.\left(\frac{k_f{k}_{r/f}}{{\mu }_{\mathrm{lf}}{B}_{\mathrm{lf}}}\nabla p_{\mathrm{lf}}\right)+q_{l,\mathrm{fm}}-q_{l,\mathrm{Ff}}=\frac{\partial }{\partial t}\left(\frac{{\varphi }_f{S}_{\mathrm{lf}}}{B_{\mathrm{lf}}}\right)& \left(2\right)\end{array}$ artificial fracture system seepage equation: $\begin{array}{cc}.Math.q_{l,F}+.Math.q_{l,\mathrm{FF}}+.Math.q_{l,\mathrm{Ff}}+q_{l,\mathrm{well}}=\frac{\partial }{\partial t}\left(\frac{{\varphi }_T{S}_{\mathrm{lF}}}{B_{\mathrm{lF}}}\right)& \left(3\right)\end{array}$ supplementary equation is as follows: $\begin{array}{cc}{q}_{l,\mathrm{fm}}=\frac{\alpha k_m{k}_{\mathrm{rlm}}}{{\mu }_{\mathrm{lm}}{\beta }_{\mathrm{lm}}}\left(p_{\mathrm{lm}}-p_{\mathrm{lf}}\right)& \left(4\right)\end{array}$ $q_{l,F}=G_F\frac{1}{{\mu }_{\mathrm{lF}}{\beta }_{\mathrm{lF}}}\left(p_{\mathrm{lF}}-p_{\mathrm{lF}}\right)$ $q_{l,\mathrm{FF}}=G_{\mathrm{FF}}\frac{1}{{\mu }_{\mathrm{lF}}{\beta }_{\mathrm{lF}}}\left(p_{\mathrm{lm}}-p_{\mathrm{lf}}\right)$ $q_{l,\mathrm{Ff}}=G_{\mathrm{Ff}}\frac{1}{{\mu }_{\mathrm{lf}}{\beta }_{\mathrm{lf}}}\left(p_{\mathrm{lf}}-p_{\mathrm{lF}}\right)$ $q_{l,\mathrm{well}}=G_{\mathrm{well}}\frac{1}{{\mu }_{\mathrm{lf}}{\beta }_{\mathrm{lf}}}\left(p_{\mathrm{lF}}-p_{\mathrm{wf}}\right)$ wherein k.sub.m and k.sub.f represent the permeability of matrix, and micro-fracture, respectively, the unit is D; l=g or w, where g represents gas phase and w represents liquid phase; k.sub.rlm and k.sub.rlf represent the relative permeability of matrix, and micro-fracture, the unit is D; μ.sub.lm, μ.sub.lf, μ.sub.lF represent the fluid viscosity of matrix, micro-fracture and artificial fracture system, the unit is Pa.Math.s; B.sub.lm, B.sub.lf, B.sub.lF represent the fluid volume coefficients of the matrix, micro-fracture, and artificial fracture systems, respectively, dimensionless; φ.sub.m, φ.sub.f, φ.sub.F represent the porosity of the matrix, micro-fractures, and artificial fracture systems, respectively, dimensionless; S.sub.lm, S.sub.lf, S.sub.lF represent saturation, dimensionless; t represents production time, the unit is day; p.sub.lm, p.sub.lf, p.sub.lF represent the two-phase pressure of matrix, micro-fracture and artificial fracture system, the unit is MPa; q.sub.l,fm represents the flow exchange between matrix and micro-fractures based on the dual medium model, the unit is 1/day, wherein 1=g,w; α represents the shape factor between matrix and micro-fractures, the unit is m.sup.−2; q.sub.l,F represents the flow exchange between the grids after an artificial fracture is divided into multiple grids by the matrix grid, the unit is 1/day; q.sub.l,FF represents the flow exchange between different artificial fracture grids, the unit is 1/day; q.sub.l,Ff represents the flow exchange between the artificial fracture grid and the matrix grid, the unit is 1/day; q.sub.l,well represents the flow exchange between the artificial fracture grid and the connected well grid, the unit is 1/day; G.sub.F represents geometric parameters for the connection pairs of grids in the same fracture, dimensionless; G.sub.FF represents the geometric parameter of the connection pairs of intersected fracture grids between different fractures, dimensionless; G.sub.Ff is the geometric parameter of the connection pair between the matrix grid and intersected fracture grid, dimensionless; G.sub.well represents geometric parameters for connected pairs between fracture grid and intersected well grid, dimensionless; S12, use a three-dimensional cubic grid to perform a space grid discretization for the reservoir, and based on the embedded discrete fracture model, the fracturing fractures are embedded in the three-dimensional cubic grid in a two-dimensional plane; S13, using a finite difference method to establish the formation flow fully implicit numerical model based on the embedded discrete fracture model, comprising steps: for the seepage equation (1) of the matrix system, the flow exchange relationship between the matrix grids on the left is written in the flow form: $\begin{array}{cc}.Math.q_{l,\mathrm{mm}}-q_{l,\mathrm{fm}}=\frac{\partial }{\partial t}\left(\frac{{\varphi }_m{S}_{\mathrm{lm}}}{B_{\mathrm{lm}}}\right)& \left(5\right)\end{array}$ in the formula: q.sub.l,mm is the flow exchange between connected matrix grids, the unit is 1/day; the finite difference method is used to derive and expand the flow exchange term on the left of the above seepage equation, if there is a flow exchange between i and j grids, the fully implicit iterative form of the flow exchange between the two grids is: $\begin{array}{cc}{q}_{\mathrm{lij},\mathrm{mm}}^{n+1}=a_{\mathrm{lmi}}\delta p_{\mathrm{lmi}}+a_{\mathrm{lmj}}\delta p_{\mathrm{lmj}}+b_{\mathrm{lmi}}\delta S_{\mathrm{lmwi}}+b_{\mathrm{imj}}\delta S_{\mathrm{lmwj}}+d_{\mathrm{imij}}& \left(6\right)\end{array}$ in the formula: n+1 representing the next time step; δp.sub.lmi, δp.sub.lmj representing the iterative pressure variables in the i-th and j-th grids of the matrix system; a.sub.lmi, a.sub.lmj representing coefficients of iterative pressure variables in the i-th and j-th grids of the matrix system; δS.sub.lmi, δS.sub.lmj representing the iterative saturation variable in the i-th and j-th grids of the matrix system; b.sub.lmi, b.sub.lmj representing coefficients the iterative saturation variables in the i-th and j-th grids of the matrix system; d.sub.lmij representing the known quantity of flow exchange between the i-th and j-th grids of the matrix system; the mass accumulation term on the right of the above seepage equation also derives and expands the unknowns, the coefficients of the unknowns are added to the coefficient matrix, and the known quantities are brought into the known quantity matrix to obtain:
mm×δ.sub.m=d.sub.m  (7) in the formula: mm representing the coefficient matrix in the matrix flow equation; δ.sub.m representing the unknown quantity matrix in the matrix flow equation; d.sub.m representing the matrix of known quantities in the matrix flow equation; in the same way obtaining the coefficient matrix ff, FF, ww in the micro-fracture, artificial fracture and wellbore flow equation; the flow exchange coefficient matrix mf between the matrix and the micro-fracture system; the flow exchange coefficient matrix fF between the micro-fracture and the artificial fracture; flow exchange coefficient matrix Fw between artificial fracture and wellbore; known quantity matrix d.sub.f, d.sub.F, d.sub.w in micro-fracture, artificial fracture, and wellbore flow equation; unknown matrix δf, δF, δw in micro-fractures, artificial fractures, and wellbore flow equations; establish a fully implicit solver matrix based on the finite difference method: $\begin{array}{cc}\left[\begin{array}{cccc}\mathrm{mm}& \mathrm{mf}& 0& 0\\ \mathrm{mf}& ff& \mathrm{jF}& 0\\ 0& \mathrm{jF}& \mathrm{FF}& \mathrm{Fw}\\ 0& 0& \mathrm{Fw}& \mathrm{ww}\end{array}\right]\left[\begin{array}{c}\delta m\\ \delta f\\ \delta F\\ \delta w\end{array}\right]=\left[\begin{array}{c}{d}_m\\ d_f\\ d_F\\ d_w\end{array}\right]& \left(8\right)\end{array}$ equation (8) is solved to obtain the gas production and liquid production of each fracture, the grid block pressure where each fracture is located, and the initial bottom-hole flow pressure.

3. The full wellbore pressure calculation method of fracturing horizontal well according to claim 2, it is characterized in that, the step S3 comprises the following steps: S31, calculate apparent flow rates of gas and liquid at the initiation point of each fracture from the gas production and liquid production of each fracture, and combine the apparent flow rates of gas and liquid into the Mukherjee-Brill liquid holdup formula to calculate a liquid holdup at each fracture initiation point: $$\begin{gathered} \begin{array}{cc}{v}_{sg}=q_{sc}{B}_g/A \\ {v}_{\mathrm{sl}}=q_w{B}_w/A& \left(9\right)\end{array} \end{gathered}$$ wherein v.sub.sl represents the apparent flow rate of the liquid phase, the unit is m/s; v.sub.sg representing the apparent flow rate of the gas phase, the unit is m/s; q.sub.sc represents the fracture gas production, the unit m.sup.3/s; q.sub.w represents the fracture fluid production, the unit is m.sup.3/s; B.sub.g represents the volume coefficient of the gas phase, dimensionless; B.sub.w represents the liquid phase volume coefficient, dimensionless; A represents the cross-sectional area of the wellbore, the unit is m.sup.2; $\begin{array}{cc}{H}_l=\exp \left[\left(c_1+c_2\sin \theta +c_3{\sin }^2\theta +c_4{N}_l^2\right)\frac{N_{gv}^{c_5}}{N_{\mathrm{lv}}^{c_6}}\right]& \left(10\right)\end{array}$ in the formula: H.sub.l represents the liquid holdup, dimensionless; c.sub.1=−0.380113;c.sub.2=−0.129875;c.sub.3=−0.119788;c.sub.4=2.343227;c.sub.5=0.475686; c.sub.6=0.288657; θ represents the angle between the wellbore and the horizontal direction, the unit is °; the N.sub.lv, N.sub.gv, N.sub.l in the above formula are as follows: $\begin{array}{cc}\mathrm{Liquid}\mathrm{velocity}\mathrm{number}\text{:}{N}_{lv}={v_{sl}\left(\frac{{\rho }_l}{G\delta }\right)}^{1/4}& \left(11\right)\end{array}$ $\mathrm{Gas}\mathrm{velocity}\mathrm{number}\text{:}{N}_{gv}={v_{sg}\left(\frac{{\rho }_l}{G\delta }\right)}^{1/4}$ $\mathrm{Liquid}\mathrm{viscosity}\mathrm{number}\text{:}{N}_l={{\mu }_l\left(\frac{G}{{\rho }_l{\delta }^3}\right)}^{1/4}$ in the formula: ρ.sub.l represents the liquid density, the unit is kg/m.sup.3; δ represents the gas-liquid interfacial tension, the unit is N/m; μ.sub.l repesents the liquid viscosity, the unit is Pa.Math.s; S32, calculate the density ρ.sub.m of the mixture, friction coefficient f.sub.m from the liquid holdup at each fracture initiation point and the grid block pressure, and the pressure gradient equation considering gravity loss, frictional resistance loss and fracture convergence loss is established as follows: $\begin{array}{cc}\frac{\mathrm{dp}}{\mathrm{dL}}={\rho }_{m}G\sin \theta +f_m\frac{{\rho }_m{v}_m^2}{2D}+{\rho }_m{v}_m\frac{{\mathrm{dv}}_m}{\mathrm{dL}}& \left(12\right)\end{array}$ wherein dp/dL represents the pressure gradient, the unit is MPa/m; dL represents the length of the horizontal section between adjacent fractures, the unit is m; ρ.sub.m represents the mixed density, the unit is kg/m.sup.3; G represents the acceleration of gravity, the unit is m/s.sup.2; θ represents the angle between the wellbore and the horizontal direction; f.sub.m represents the two-phase friction coefficient, dimensionless; v.sub.m represents the mixing speed, the unit is m/s; D represents the diameter of the pipe, the unit is m; S33, combine the pressure gradient equation (12) and the distance ΔL.sub.i from the i-th fracture initiation point to the bottom-hole, to obtain the wellbore pressure changes in each fracturing section of the horizontal well as follows: $\begin{array}{cc}\Delta p_{\mathrm{wi}}=\Delta L_i.Math.\frac{\mathrm{dp}}{\mathrm{dL}}\left(i=2,3,.Math.N_F\right)& \left(13\right)\end{array}$ in the formula: NF represents the total number of fractures; ΔL.sub.i represents the distance from the i-th fracture initiation point to the bottom-hole of the well, the unit is m; Δpwi represents the horizontal wellbore pressure change from the bottom-hole to the i-th fracture, the unit is MPa.

4. A device for calculating the full wellbore pressure of a fracturing horizontal well, it is characterized in that, comprising: processor, and acquisition module configured to acquire the initial data for calculating the full wellbore pressure of the fracturing horizontal well; storage module having stored thereon a program operable on the processor for calculating the full wellbore pressure of a fracturing horizontal well, the program for calculating the full wellbore pressure of a fracturing horizontal well being processed by the processor implements the steps of the method according to any one of claims 1-3; output module configured to output calculation results.

5. A non-transitory tangible computer readable storage medium, in which program code executable by a processor is stored, it is characterized in that, the computer readable storage medium comprises a plurality of instructions configured to cause the processor to perform the method for calculating the full wellbore pressure of a fractured horizontal well of any one of claims 1-3.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is a schematic diagram of the discrete spatialgrid of the gas well based on EDFM in the embodiment of the present invention;

(2) FIG. 2 is reservoir pressure distribution at different production times of the gas well in the embodiment of the present invention;

(3) FIG. 3 is the prediction curve of the full wellbore pressure distribution at different production times of the gas well according to the embodiment of the present invention;

(4) FIG. 4 is a flow chart of the method for calculating the full wellbore pressure of a fracturing horizontal well according to an embodiment of the present invention; FIG. 4A is a flow chart of establishment of a formation flow fully implicit numerical model according to an embodiment of the present invention; and FIG. 4B is a flow chart of obtaining a wellbore pressure variation of each fracturing section of the horizontal well according to an embodiment of the present invention.

(5) FIG. 5 is a device for calculating the full wellbore pressure of a fracturing horizontal well for the embodiment of the present invention.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

(6) In order to have a clearer understanding of the technical features, purposes and beneficial effects of the present invention, the technical solutions of the present invention are now described in detail below, but should not be construed as limiting the scope of implementation of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by persons of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

Example

(7) To calculate the full wellbore pressure of fracturing horizontal wells, the flow chart of this calculation method can be referred to FIG. 4, including the following steps:

(8) S1, collect wellbore structural parameters, fluid physical properties, reservoir parameters and production history data at different time steps, and establish a fully implicit numerical model of formation flow based on embedded discrete fracture model (EDFM);

(9) This example takes an actual shale gas well as an example, and the specific basic data are as follows:

(10) TABLE-US-00001 TABLE 1 Basic data of an actual shale gas well Numerical Numerical Parameter value Parameter value Outer zone 1400*400*50 Formation 40 volume(m.sup.3) pressure(MPa) Length of main 200 Daily gas 15 fracture(m) production(10,000 cubic meters/day) Horizontal section 1000 Main fracture 500 length(m) permeability(mD) The number of 11 Micro-fracture 0.005 main fracture permeability(mD) Matrix 0.0001 Matrix porosity 5% permeability(mD) Stress Sensitivity 0.1 Micro-fracture 1% coefficient(MP.sup.−1) porosity Sleeve size(m) 0.1214 Shape factor(m.sup.−2) 0.01

(11) The establishment of the formation flow fully implicit numerical model based on the embedded discrete fracture model (EDFM) includes the following steps:

(12) S11, establish a formation flow mathematical model based on the embedded discrete fracture model, use the dual medium model to describe the flow exchange between the matrix and micro-fractures, and use the embedded discrete fracture model to describe the flow exchange between the matrix and artificial fractures, based on the mass conservation equation, the seepage differential equation is established with Darcy's law, and the formation flow mathematical model of the embedded discrete fracture model is as follows:

(13) Matrix system seepage equation:

(14) $\begin{array}{cc}\nabla .Math.\left(\frac{k_m{k}_{rlm}}{{\mu }_{lm}{B}_{lm}}\nabla p_{lm}\right)-q_{l,\mathrm{fm}}=\frac{\partial }{\partial t}\left(\frac{{\varphi }_m{S}_{lm}}{B_{lm}}\right)& \left(1\right)\end{array}$

(15) The seepage equation of the micro-fracture system:

(16) $\begin{array}{cc}\nabla .Math.\left(\frac{k_f{k}_{r/f}}{{\mu }_{lf}{B}_{lf}}\nabla p_{lf}\right)+.Math.q_{l,\mathrm{fm}}-q_{l,\mathrm{Ff}}=\frac{\partial }{\partial t}\left(\frac{{\varphi }_f{S}_{\mathrm{lf}}}{B_{\mathrm{lf}}}\right)& \left(2\right)\end{array}$

(17) Seepage equation of artificial fracture system:

(18) $\begin{array}{cc}.Math.q_{l,F}+.Math.q_{l,\mathrm{FF}}+.Math.q_{l,\mathrm{Ff}}+q_{l,\mathrm{well}}=\frac{\partial }{\partial t}\left(\frac{{\varphi }_F{S}_{lF}}{B_{lF}}\right)& \left(3\right)\end{array}$

(19) The supplementary equation is as follows:

(20) $\begin{array}{cc}{q}_{l,\mathrm{fm}}=\frac{\alpha k_m{k}_{\mathrm{rlm}}}{{\mu }_{\mathrm{lm}}{B}_{\mathrm{lm}}}\left(p_{lm}-p_{\mathrm{lf}}\right)& \left(4\right)\end{array}$$q_{l,F}=G_F\frac{1}{{\mu }_{\mathrm{lF}}{B}_{\mathrm{lF}}}\left(p_{\mathrm{lf}}-p_{\mathrm{lF}}\right)$$q_{l,\mathrm{FF}}=G_{\mathrm{FF}}\frac{1}{{\mu }_{\mathrm{lF}}{B}_{\mathrm{lF}}}\left(p_{\mathrm{lm}}-p_{\mathrm{lf}}\right)$$q_{l,\mathrm{Ff}}=G_{\mathrm{Ff}}\frac{1}{{\mu }_{\mathrm{lf}}{B}_{\mathrm{lf}}}\left(p_{\mathrm{lf}}-p_{\mathrm{lF}}\right)$$q_{l,\mathrm{well}}=G_{\mathrm{well}}\frac{1}{{\mu }_{\mathrm{lF}}{B}_{\mathrm{lF}}}\left(p_{\mathrm{lF}}-p_{\mathrm{wf}}\right)$

(21) Wherein k.sub.m and k.sub.f represent the permeability of matrix, micro-fracture, and artificial fracture system, respectively, the unit is D; k.sub.rlm and k.sub.rlf represent the relative permeability of matrix, micro-fracture and artificial fracture system, the unit is D, wherein 1=g,w; μ.sub.lm, μ.sub.lf, μ.sub.lF represent the fluid viscosity of matrix, micro-fracture and artificial fracture system, the unit is Pa.Math.s, wherein 1=g,w; B.sub.lm, B.sub.lf, B.sub.lF represent the fluid volume coefficients of the matrix, micro-fracture, and artificial fracture systems, respectively, dimensionless, wherein 1=g,w; φ.sub.m, φ.sub.f, φ.sub.F represent the porosity of the matrix, micro-fractures, and artificial fracture systems, respectively, dimensionless; S.sub.lm, S.sub.lf, S.sub.lF represent saturation, dimensionless, wherein 1=g,w; t represents production time, the unit is day; p.sub.lm, p.sub.lf, p.sub.lF represent the two-phase pressure of matrix, micro-fracture and artificial fracture system, the unit is MPa, wherein 1=g,w; r.sub.l represents the phase density, the unit is kg/m.sup.3, wherein 1=g,w; G represents the acceleration of gravity, the unit is m/s.sup.2; z represents vertical depth, the unit is m; q.sub.l,fm represents the flow exchange between matrix and micro-fractures based on the dual medium model, the unit is 1/day, wherein 1=g,w; α represents the shape factor between matrix and micro-fractures, the unit is m.sup.−2; q.sub.l,F represents the flow exchange between the grids after an artificial fracture is divided into multiple grids by the matrix grid, the unit is 1/day; q.sub.l,FF represents the flow exchange between different artificial fracture grids, the unit is 1/day; q.sub.l,Ff represents the flow exchange between the artificial fracture grid and the matrix grid, the unit is 1/day; q.sub.l,well represents the flow exchange between the artificial fracture grid and the connected well grid, the unit is 1/day; G.sub.F represents geometric parameters for the connection pairs of grids in the same fracture; G.sub.FF represents the geometric parameter of the connection pairs of intersected fracture grids between different fractures; G.sub.Ff is the geometric parameter of the connection pair between the matrix grid and intersected fracture grid; G.sub.well represents geometric parameters for connected pairs between fracture grid and intersected well grid.

(22) S12, use a three-dimensional cubic grid to discretize the reservoir in space, based on EDFM, the fractures are embedded in the three-dimensional cubic grid in a two-dimensional plane.

(23) S13, use the finite difference method to establish the formation flow fully implicit numerical model based on the embedded discrete fracture model (EDFM), and solve the formation flow fully implicit numerical model; include the following steps:

(24) For the seepage equation (1) of the matrix system, the flow exchange relationship between the matrix grids on the left can be written in the form of flow:

(25) $\begin{array}{cc}.Math.q_{,\mathrm{lmm}}-q_{l,\mathrm{fm}}=\frac{\partial }{\partial t}\left(\frac{{\varphi }_m{S}_{lm}}{B_{lm}}\right)& \left(5\right)\end{array}$

(26) In the formula: q.sub.l,mm represents the flow exchange between connected matrix grids, the unit is 1/day.

(27) The finite difference method is used to derivative expansion of the flow exchange term Σq.sub.i,mm−Σq.sub.j,mm on the left side of the above seepage equation, if there is flow exchange between i and j grids, the fully implicit iterative form of the flow exchange between the two grids is.

(28) $\begin{array}{cc}{q}_{\mathrm{lij},\mathrm{mm}}^{n+1}=a_{\mathrm{lmi}}\delta p_{\mathrm{lmi}}+a_{\mathrm{lmj}}\delta p_{\mathrm{lmj}}+b_{\mathrm{lmi}}\delta S_{\mathrm{lmwi}}+b_{\mathrm{imj}}\delta S_{\mathrm{lmwj}}+d_{\mathrm{imij}}& \left(6\right)\end{array}$

(29) In the formula: n+1 represents the next time step; v represents the v-th iteration step; δp.sub.lmi, δp.sub.lmj represents the iterative pressure variables in the i-th and j-th grids of the matrix system; a.sub.lmi, a.sub.lmj representing coefficients of iterative pressure variables in the i-th and j-th grids of the matrix system; δS.sub.lmi, δS.sub.lmj represents the iterative saturation variable in the i-th and j-th grids of the matrix system; b.sub.lmi, b.sub.lmj representing coefficients the iterative saturation variables in the i-th and j-th grids of the matrix system; d.sub.lmij is the known quantity of flow exchange between the i-th and j-th grids of the matrix system.

(30) The mass accumulation term

(31) $\frac{\partial }{\partial t}\left(\frac{{\varphi }_m{S}_{lm}}{B_{lm}}\right)$
on the right side of the above seepage equation also derives and expands the unknowns, the coefficients of the unknowns are added to the coefficient matrix, and the known quantities are brought into the known quantity matrix to obtain:
mm×δ.sub.m=d.sub.m  (7)

(32) In the formula: mm represents the coefficient matrix in the matrix flow equation; δ.sub.m represents the unknown quantity matrix in the matrix flow equation; d.sub.m is the matrix of known quantities in the matrix flow equation. Therefore, the fully implicit solution matrix based on the finite difference method for micro-fractures, artificial fractures, the flow exchange terms of each system and the internal flow of the wellbore can be deduced as follows:

(33) $\begin{array}{cc}\left[\begin{array}{cccc}\mathrm{mm}& \mathrm{mf}& 0& 0\\ \mathrm{mf}& \mathrm{ff}& \mathrm{fF}& 0\\ 0& \mathrm{fF}& \mathrm{FF}& \mathrm{Fw}\\ 0& 0& \mathrm{Fw}& \mathrm{ww}\end{array}\right]\left[\begin{array}{c}\delta m\\ \delta f\\ \delta F\\ \delta w\end{array}\right]=\left[\begin{array}{c}{d}_m\\ d_f\\ d_F\\ d_w\end{array}\right]& \left(8\right)\end{array}$

(34) Wherein ff, FF, and ww represent the coefficient matrices in the microfractures, artificial fractures, and wellbore flow equations, respectively; mf, fm represent the flow exchange coefficient matrices between the matrix and the microfracture system; fF, Ff represent the microfractures and artificial fractures; Fw represents the flow exchange coefficient matrix between the artificial fracture and the wellbore; d.sub.f, d.sub.F, d.sub.w represent the known quantity matrices in the microfracture, artificial fracture, and wellbore flow equations, respectively; δf, δF, and δw represent the unknown quantity matrices in the microfracture, artificial fracture, and wellbore flow equations, respectively.

(35) By solving Equation (8), the gas production and liquid production of each fracture, the grid block pressure where each fracture is located, and the initial bottom-hole flow pressure can be obtained.

(36) S2, solve the formation flow fully implicit numerical model based on the embedded discrete fracture model (EDFM) under an inner boundary condition of a constant gas production rate at a certain time step, and obtain the gas production, liquid production, grid block pressure and initial bottom-hole flow of each fracture; when this step is performed, in the first iteration step, the pressure drop change in the horizontal section is not considered; in the subsequent iteration steps, the pressure drop change in the horizontal section is considered.

(37) S3, bring the gas production, liquid production, and the pressure of the grid block where the fracture initiation point is located into the pressure drop model considering the gravity loss, frictional resistance loss and fracture convergence loss, and obtain the wellbore pressure change of each fracturing section of the horizontal well. Specifically include the following steps:

(38) S31, calculate the apparent flow rate of gas and liquid at the initiation point of each fracture from the gas production and liquid production of each fracture, and bring it into the Mukherjee-Brill (M-B) liquid holdup formula to calculate the liquid holdup at each fracture initiation point:

(39) $$\begin{gathered} \begin{array}{cc}{v}_{sg}=q_{sc}{B}_g/A \\ {v}_{s1}=q_w{B}_w/A& \left(9\right)\end{array} \end{gathered}$$

(40) Wherein v.sub.sl represents the apparent flow rate of the liquid phase, the unit is m/s; v.sub.sg represents the apparent flow rate of the gas phase, the unit is m/s; q.sub.sc represents the fracture gas production, the unit m.sup.3/s; q.sub.w represents the fracture fluid production, the unit is m.sup.3/s; B.sub.g represents the volume coefficient of the gas phase, dimensionless; B.sub.w represents the liquid phase volume coefficient, dimensionless; A represents the cross-sectional area of the wellbore, the unit is m.sup.2.

(41) 0$\begin{array}{cc}{H}_l=\exp \left[\left(c_1+c_2\sin \theta +c_3{\sin }^2\theta +c_4{N}_l^2\right)\frac{N_{gv}^{c_5}}{N_{\mathrm{lv}}^{c_6}}\right]& \left(10\right)\end{array}$

(42) In the formula: H.sub.l represents the liquid holdup, dimensionless; c.sub.1=−0.380113;c.sub.2=−0.129875;c.sub.3=−0.119788;c.sub.4=2.343227;c.sub.5=0.475686; c.sub.6=0.288657; θ represents the angle between the wellbore and the horizontal direction, the unit is degree (°).

(43) The N.sub.lv, N.sub.gv, N.sub.l in the above formula are as follows:

(44) $\begin{array}{cc}\mathrm{Liquid}\mathrm{velocity}\mathrm{number}\text{:}{N}_{lv}={v_{sl}\left(\frac{{\rho }_l}{G\delta }\right)}^{1/4}& \left(11\right)\end{array}$$\mathrm{Gas}\mathrm{velocity}\mathrm{number}\text{:}{N}_{gv}={v_{sg}\left(\frac{{\rho }_l}{G\delta }\right)}^{1/4}$$\mathrm{Liquid}\mathrm{viscosity}\mathrm{number}\text{:}{N}_l={{\mu }_l\left(\frac{G}{{\rho }_l{\delta }^3}\right)}^{1/4}$

(45) In the formula: ρl reprents the liquid density, the unit is kg/m.sup.3; δ represents the gas-liquid interfacial tension, the unit is N/m; μ.sub.l reprents the liquid viscosity, the unit is Pa.Math.s.

(46) S32, calculate the density ρ.sub.m of the mixture, friction coefficient f.sub.m from the liquid holdup at each fracture initiation point and the grid block pressure, and the pressure gradient equation considering gravity loss, frictional resistance loss and fracture convergence loss is established as follows:

(47) $\begin{array}{cc}\frac{\mathrm{dp}}{\mathrm{dL}}={\rho }_{m}G\sin \theta +f_m\frac{{\rho }_m{v}_m^2}{2D}+{\rho }_m{v}_m\frac{d{v}_m}{\mathrm{dL}}& \left(12\right)\end{array}$

(48) In the formula: (dp/dL) represents the pressure gradient, the unit is MPa/m; dL represents the length of the horizontal section between adjacent fractures, the unit is m; ρ.sub.m represents the mixed density, the unit is kg/m.sup.3; G represents the acceleration of gravity, the unit is m/s.sup.2; θ represents the angle between the wellbore and the horizontal direction; f.sub.m represents the two-phase friction coefficient, dimensionless; v.sub.m represents the mixing speed, the unit is m/s; D represents the diameter of the pipe, the unit is m.

(49) S33, combine the pressure gradient equation (12) and the distance ΔL.sub.i from the i-th fracture initiation point to the bottom-hole (i. e, the first fracture initiation point), the wellbore pressure changes in each fracturing section of the horizontal well are obtained as follows:

(50) $\begin{array}{cc}\Delta p_{\mathrm{wi}}=\Delta L_i.Math.\frac{\mathrm{dp}}{\mathrm{dL}}\left(i=2,3,.Math.N_F\right)& \left(13\right)\end{array}$

(51) In the formula: N.sub.F represents the total number of the fractures; ΔL.sub.i represents the distance from the i-th fracture initiation point to the bottom-hole of the well, the unit is m; Δp.sub.wi represents the horizontal wellbore pressure change from the bottom-hole to the i-th fracture, the unit is MPa.

(52) S4, calculate the bottom-hole pressure of each fracture initiation point from the initial bottom-hole flow pressure and the wellbore pressure change of each fracturing section;

(53) S5, bring the bottom-hole flow pressure and the bottom-hole pressure of each fracture initiation point into the formation flow fully implicit numerical model based on the embedded discrete fracture model (EDFM) in step S2, and iteratively cycle steps S2-S4 until the iterative variables converge, obtain the bottom-hole flow pressure and bottom-hole pressure of each fracture initiation point at this time step;

(54) S6, bring the bottom-hole flow pressure and the total gas production and liquid production at the time step into the pressure drop model considering gravity loss and frictional resistance loss in the inclined section and vertical well section, and obtain the wellbore pressure and wellhead casing pressure of each section through the iterative cycle section by section.

(55) S7, for each time step, repeat steps S2-S6 to obtain the full wellbore pressure under different production times and draw a pressure distribution curve. For specific results, see FIG. 2 and FIG. 3.

(56) Referring to FIG. 5, a device for calculating the full wellbore pressure of a fracturing horizontal well, includes a processor, an acquisition module for acquiring initial data for calculating the full wellbore pressure of a fracturing horizontal well, an output module and a storage module for outputting calculation results; a program for calculating the full wellbore pressure of a fracturing horizontal well is stored on the storage module, which can be run on the processor, and the program for calculating the full wellbore pressure of a fracturing horizontal well is executed by the processor to implement the above method step.

(57) A computer readable storage medium having stored therein program code executable by a processor, the computer readable storage medium including a plurality of instructions configured to cause the processor to perform the described fracturing horizontal well process. wellbore pressure calculation method.

(58) The above is not intended to limit the present invention in any form. Although the present invention has been disclosed through the above-mentioned embodiments, it is not intended to limit the present invention. Any person skilled in the art, within the scope of the technical solution of the present invention, When the technical contents disclosed above can be used to make some changes or modifications to equivalent embodiments with equivalent changes, any simple modifications or equivalents to the above embodiments according to the technical essence of the present invention do not depart from the content of the technical solution of the present invention. Changes and modifications still fall within the scope of the technical solutions of the present invention.