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.

Description

FIELD OF THE INVENTION

The invention relates to the technical field of numerical simulation of oil and gas reservoirs, in particular to a method, device and computer-readable storage medium for calculating a full wellbore pressure for fracturing horizontal wells.

BACKGROUND OF THE INVENTION

Unconventional shale gas reservoirs, tight gas reservoirs and other original reservoirs generally have low porosity and permeability, and low natural depletion production. In order to obtain industrial production capacity, it is necessary to generate multiple artificial fractures through long horizontal section plus hydraulic fracturing technology to improve the seepage capacity of unconventional reservoirs. In the past research, it is generally believed that the production of each artificial fracture is equal, and the pressure loss in the wellbore is ignored, and the fracturing horizontal well is regarded as having infinite conductivity. However, the actual situation is that due to the heterogeneity of the reservoir and the interference between fractures in the production process, the production of each artificial fracture is not equal. In addition, due to the confluence of reservoir fluid, the flow in the horizontal section of the wellbore belongs to variablemass multi-phase flow, which will generate frictional pressure drop, gravity pressure drop, and accelerated pressure drop caused by fluid confluence in fractures. The pressure change in the wellbore will in turn affect the fluid confluence of each fracture, resulting in the change of the productivity of each fracture. Therefore, it is necessary to study the coupled formation wellbore flow model of fracturing horizontal wells to achieve accurate prediction of fracture production, reservoir pressure variation and the pressure distribution characteristics of the full wellbore in fracturing horizontal wells.

SUMMARY OF THE INVENTION

In view of the above reasons, the purpose of the present invention is to provide a method for calculating the full wellbore pressure of a fracturing horizontal well. The technical solution of the present invention is as follows.

A method for calculating a full wellbore pressure of a fracturing horizontal well, includes the following steps:

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).

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 well bottom flow pressure of each fracture.

S3, bring the gas production volume, liquid production volume 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 variation of each fracturing section of the horizontal well.

S4, calculate the well bottom-hole pressure of each fracture initiation point from the initial well bottom-hole flow pressure and the wellbore pressure change of each fracturing section.

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 the 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.

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 each section wellbore pressure and wellhead casing pressure through the iterative cycle step by step.

S7, for each time step, repeat steps S2-S6 to obtain the full wellbore pressure at different production times.

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 for outputting calculation results, and a storage module; a program for calculating the full wellbore pressure of a fracturing horizontal well that can be run on the processor is stored thereon, and the program for calculating the full wellbore pressure of a fracturing horizontal well is executed by the processor to realize the above-mentioned method steps.

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 wellbore pressure calculation method of fracturing horizontal well.

Compared with the prior art, the present invention has the following beneficial effects:

1. based on the dual medium-embedded discrete fracture coupling model, a fully implicit numerical model for the seepage flow of fracturing horizontal wells in unconventional tight gas reservoirs and shale gas reservoirs is established, which realizes the accurate characterization of multi-medium flow laws in unconventional gas reservoirs;

2. considering the pressure drop, frictional resistance and gravity pressure drop caused by the mass-varying flow of formation fluid into the wellbore, a wellbore pressure drop model for fracturing horizontal wells in unconventional gas reservoirs with gas-liquid two-phase flow was established, and the wellbore pressure drop model for fracturing horizontal wells in unconventional gas reservoirs was established. The full wellbore pressure distribution characteristics are simulated and predicted, and in the absence of bottom-hole pressure test data, the production data can be used to predict the change law of the whole wellbore pressure and wellhead casing pressure, which saves the test cost and has important practical significance.

BRIEF DESCRIPTION OF THE DRAWINGS

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

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

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;

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.

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

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

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:

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);

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

TABLE 1
Basic data of an actual shale gas well
	Numerical		Numerical
Parameter	value	Parameter	value
Outer zone volume	1400*400*50	Formation	40
(m3)		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−1)		porosity
Sleeve size(m)	0.1214	Shape factor(m−2)	0.01

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

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:

Matrix system seepage equation:

$\begin{array}{cc}\nabla ·\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}$

The seepage equation of the micro-fracture system:

$\begin{array}{cc}\nabla ·\left(\frac{k_f{k}_{\mathrm{rlf}}}{{\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}$

Seepage equation of artificial fracture system:

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

The supplementary equation is as follows:

$\begin{array}{cc}{q}_{l,\mathrm{fm}}=\frac{\alpha k_m{k}_{\mathrm{rlm}}}{{\mu }_{\mathrm{lm}}{B}_{\mathrm{lm}}}\left(p_{\mathrm{lm}}-p_{\mathrm{lf}}\right)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)& \left(4\right)\end{array}$

Wherein k_m and k_f represent the permeability of matrix, micro-fracture, and artificial fracture system, respectively, the unit is D; k_rlm and k_rlf represent the relative permeability of matrix, micro-fracture and artificial fracture system, the unit is D, wherein l=g,w; μ_lm, μ_lf, μ_lF represent the fluid viscosity of matrix, micro-fracture and artificial fracture system, the unit is Pa·s, wherein l=g,w; B_lm, B_lf, B_lF represent the fluid volume coefficients of the matrix, micro-fracture, and artificial fracture systems, respectively, dimensionless, wherein l=g,w; μ_m, μ_f, μ_F represent the porosity of the matrix, micro-fractures, and artificial fracture systems, respectively, dimensionless; S_lm, S_lf, S_lF represent saturation, dimensionless, wherein l=g,w; t represents production time, the unit is day; p_lm, p_lf, p_lF represent the two-phase pressure of matrix, micro-fracture and artificial fracture system, the unit is MPa, wherein l=g,w; r_l represents the phase density, the unit is kg/m³, wherein l=g,w; G represents the acceleration of gravity, the unit is m/s²; z represents vertical depth, the unit is m; q_l,fm represents the flow exchange between matrix and micro-fractures based on the dual medium model, the unit is 1/day, wherein l=g,w; α represents the shape factor between matrix and micro-fractures, the unit is m⁻²; q_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_l,FF represents the flow exchange between different artificial fracture grids, the unit is 1/day; q_l,Ff represents the flow exchange between the artificial fracture grid and the matrix grid, the unit is 1/day; q_l,well represents the flow exchange between the artificial fracture grid and the connected well grid, the unit is 1/day; G_F represents geometric parameters for the connection pairs of grids in the same fracture; G_FF represents the geometric parameter of the connection pairs of intersected fracture grids between different fractures; G_Ff is the geometric parameter of the connection pair between the matrix grid and intersected fracture grid; G_well represents geometric parameters for connected pairs between fracture grid and intersected well grid.

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.

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:

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:

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

In the formula: q_l,mm represents the flow exchange between connected matrix grids, the unit is 1/day.

The finite difference method is used to derivative expansion of the flow exchange term Σq_l,mm−q_l,fm 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.

q_lij,mmⁿ⁺¹=a_lmiδp_lmi+a_lmjδp_lmj+b_lmiδS_lmwi+b_lmjδS_lmwj+d_lmij  (6)

In the formula: n+1 represents the next time step; v represents the v-th iteration step; δp_lmi, δp_lmj represents the iterative pressure variables in the i-th and j-th grids of the matrix system; a_lmi, a_lmj representing coefficients of iterative pressure variables in the i-th and j-th grids of the matrix system; δS_lmi, δS_lmj represents the iterative saturation variable in the i-th and j-th grids of the matrix system; b_lmi, b_lmj representing coefficients the iterative saturation variables in the i-th and j-th grids of the matrix system; d_lmij is the known quantity of flow exchange between the i-th and j-th grids of the matrix system.

The mass accumulation term

$\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×δm=d_m  (7)

In the formula: mm represents the coefficient matrix in the matrix flow equation; δ_m represents the unknown quantity matrix in the matrix flow equation; d_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:

$\begin{array}{cc}\left[\begin{array}{cccc}mm& mf& 0& 0\\ mf& \mathrm{ff}& \mathrm{fF}& 0\\ 0& \mathrm{fF}& FF& Fw\\ 0& 0& Fw& 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}$

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_f, d_F, d_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.

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.

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.

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:

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:

v_sg=q_scB_g/A

v_sl=q_wB_w/A  (9)

Wherein v_sl represents the apparent flow rate of the liquid phase, the unit is m/s; v_sg represents the apparent flow rate of the gas phase, the unit is m/s; q_sc represents the fracture gas production, the unit m³/s; q_w represents the fracture fluid production, the unit is m³/s; B_g represents the volume coefficient of the gas phase, dimensionless; B_w represents the liquid phase volume coefficient, dimensionless; A represents the cross-sectional area of the wellbore, the unit is m².

$\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_l represents the liquid holdup, dimensionless; c₁=−0.380113; c₂=−0.129875; c₃=−0.119788; c₄=2.343227; c₅=0.475686; c₆=0.288657; θ represents the angle between the wellbore and the horizontal direction, the unit is degree (°).

The N_lv, N_gv, N_l in the above formula are as follows:

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

In the formula: μl reprents the liquid density, the unit is kg/m³; δ represents the gas-liquid interfacial tension, the unit is N/m; μ_l reprents the liquid viscosity, the unit is Pa·s.

S32, calculate the density ρ_m of the mixture, friction coefficient f_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{dp}{dL}={\rho }_{m}G\sin \theta +f_m\frac{{\rho }_m{v}_m^2}{2D}+{\rho }_m{v}_m\frac{d{v}_m}{dL}& \left(12\right)\end{array}$

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; ρ_m represents the mixed density, the unit is kg/m³; G represents the acceleration of gravity, the unit is m/s²; θ represents the angle between the wellbore and the horizontal direction; f_m represents the two-phase friction coefficient, dimensionless; v_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_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:

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

In the formula: N_F represents the total number of the fractures; ΔL_i represents the distance from the i-th fracture initiation point to the bottom-hole of the well, the unit is m; Δp_wi represents the horizontal wellbore pressure change from the bottom-hole to the i-th fracture, the unit is MPa.

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, 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;

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.

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.

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.

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.

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.

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: ∇ · ( k m ⁢ k rlm μ l ⁢ m ⁢ B l ⁢ m ⁢ ∇ p l ⁢ m ) - q l, fm = ∂ ∂ t ( ϕ m ⁢ S l ⁢ m B l ⁢ m ) ( 1 ) ∇ · ( k f ⁢ k rlf μ lf ⁢ B lf ⁢ ∇ p lf ) + q l, fm - q l, Ff = ∂ ∂ t ( ϕ f ⁢ S lf B lf ) ( 2 ) ∑ q l, F + ∑ q l, FF + ∑ q l, Ff + q l, well = ∂ ∂ t ( ϕ F ⁢ S lF B lF ) ( 3 ) q l, fm = α ⁢ k m ⁢ k rlm μ l ⁢ m ⁢ B l ⁢ m ⁢ ( p l ⁢ m - p lf ) ⁢ q l, F = G F ⁢ 1 μ lf ⁢ B lf ⁢ ( p lf - p lF ) ⁢ q l, FF = G FF ⁢ 1 μ lF ⁢ B lF ⁢ ( p l ⁢ m - p lf ) ⁢ q l, Ff = G Ff ⁢ 1 μ lf ⁢ B lf ⁢ ( p lf - p lF ) ⁢ q l, well = G well ⁢ 1 μ lf ⁢ B lf ⁢ ( p lf - p wf ) ( 4 ) ∑ q l, m ⁢ m - q l, fm = ∂ ∂ t ( ϕ m ⁢ S l ⁢ m B l ⁢ m ) ( 5 ) [ m ⁢ m m ⁢ f 0 0 m ⁢ f ff fF 0 0 fF F ⁢ F F ⁢ w 0 0 F ⁢ w w ⁢ w ] [ δ ⁢ m δ ⁢ f δ ⁢ F δ ⁢ w ] = [ d m d f d F d w ] ( 8 ) 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.

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:

micro-fracture system seepage equation:

artificial fracture system seepage equation:

supplementary equation is as follows:

wherein km and kf 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; krlm and krlf represent the relative permeability of matrix, and micro-fracture, the unit is D; μlm, μlf, μlF represent the fluid viscosity of matrix, micro-fracture and artificial fracture system, the unit is Pa·s; Blm, Blf, BlF represent the fluid volume coefficients of the matrix, micro-fracture, and artificial fracture systems, respectively, dimensionless; φm, φf, φF represent the porosity of the matrix, micro-fractures, and artificial fracture systems, respectively, dimensionless; Slm, Slf, SlF represent saturation, dimensionless; t represents production time, the unit is day; plm, plf, plF represent the two-phase pressure of matrix, micro-fracture and artificial fracture system, the unit is MPa; ql,fm represents the flow exchange between matrix and micro-fractures based on the dual medium model, the unit is 1/day, wherein l=g,w; α represents the shape factor between matrix and micro-fractures, the unit is m−2; ql,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; ql,FF represents the flow exchange between different artificial fracture grids, the unit is 1/day; ql,Ff represents the flow exchange between the artificial fracture grid and the matrix grid, the unit is 1/day; ql,well represents the flow exchange between the artificial fracture grid and the connected well grid, the unit is 1/day; GF represents geometric parameters for the connection pairs of grids in the same fracture, dimensionless; GFF represents the geometric parameter of the connection pairs of intersected fracture grids between different fractures, dimensionless; GFf is the geometric parameter of the connection pair between the matrix grid and intersected fracture grid, dimensionless; Gwell 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:

in the formula: ql,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: qlij,mmn+1=almiδplmi+almjδplmj+blmiδSlmwi+blmjδSlmwj+dlmij  (6)

in the formula: n+1 representing the next time step; δplmi, δplmj representing the iterative pressure variables in the i-th and j-th grids of the matrix system; almi, almj representing coefficients of iterative pressure variables in the i-th and j-th grids of the matrix system; δSlmi, δSlmj representing the iterative saturation variable in the i-th and j-th grids of the matrix system; blmi, blmj representing coefficients the iterative saturation variables in the i-th and j-th grids of the matrix system; dlmij 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×δm=dm  (7)

in the formula: mm representing the coefficient matrix in the matrix flow equation; δm representing the unknown quantity matrix in the matrix flow equation; dm 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 df, dF, dw 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:

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: H l = exp [ ( c 1 + c 2 ⁢ sin ⁢ θ + c 3 ⁢ sin 2 ⁢ θ + c 4 ⁢ N l 2 ) ⁢ N g ⁢ v c 5 N lv c 6 ] ( 10 ) Liquid ⁢ velocity ⁢ number: N l ⁢ v = v s ⁢ l ( ρ l G ⁢ δ ) 1 / 4 ⁢ Gas ⁢ velocity ⁢ number: N g ⁢ v = v s ⁢ g ( ρ l G ⁢ δ ) 1 / 4 ⁢ Liquid ⁢ viscosity ⁢ number: N l = μ l ( G ρ l ⁢ δ 3 ) 1 / 4 ( 11 ) d ⁢ p d ⁢ L = ρ m ⁢ G ⁢ sin ⁢ θ + f m ⁢ ρ m ⁢ v m 2 2 ⁢ D + ρ m ⁢ v m ⁢ dv m d ⁢ L ( 12 ) Δ ⁢ p wi = Δ ⁢ L i · d ⁢ p d ⁢ L ⁢ ( i = 2, 3, … ⁢ N F ) ( 13 ) in the formula: NF represents the total number of fractures; ΔLi 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.

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: vsg=qscBg/A vsl=qwBw/A  (9)

wherein vsl represents the apparent flow rate of the liquid phase, the unit is m/s; vsg representing the apparent flow rate of the gas phase, the unit is m/s; qsc represents the fracture gas production, the unit m3/s; qw represents the fracture fluid production, the unit is m3/s; Bg represents the volume coefficient of the gas phase, dimensionless; Bw represents the liquid phase volume coefficient, dimensionless; A represents the cross-sectional area of the wellbore, the unit is m2;

in the formula: Hl represents the liquid holdup, dimensionless; c1=−0.380113; c2=−0.129875; c3=−0.119788; c4=2.343227; c5=0.475686; c6=0.288657; θ represents the angle between the wellbore and the horizontal direction, the unit is °;

the Nlv, Ngv, Nl in the above formula are as follows:

in the formula: ρl represents the liquid density, the unit is kg/m3; δ represents the gas-liquid interfacial tension, the unit is N/m; μl represents the liquid viscosity, the unit is Pa·s;

S32, calculate the density ρm of the mixture, friction coefficient fm 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:

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; ρm represents the mixed density, the unit is kg/m3; G represents the acceleration of gravity, the unit is m/s2; θ represents the angle between the wellbore and the horizontal direction; fm represents the two-phase friction coefficient, dimensionless; vm 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 ΔLi 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:

4. A device for calculating the full wellbore pressure of a fracturing horizontal well, it is characterized in that, comprising: 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;

processor, and

acquisition module configured to acquire the initial data for calculating the full wellbore pressure of the fracturing horizontal well;

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.