CN104895550A - Tight gas fracturing horizontal well numerical value well testing model building and solving method - Google Patents

Abstract

The invention provides a tight gas fracturing horizontal well numerical value well testing model building and solving method. The method includes the following steps: step one, generating two-dimensional geologic body and a three-dimensional geologic body of a tight gas reservoir fracturing horizontal well; step two, performing grid discretization on the generated two-dimensional geologic body and the generate three-dimensional geologic body of the tight gas reservoir fracturing horizontal well; step three, calculating a pressure-difference-free seepage model of a horizontal well shaft; step four, building a coupling model, solving the built coupling model, and enabling acquired solutions to generate a well testing theoretical curve; step five, fitting the theoretical curve acquired in step four with a measured curve to acquire parameters of well testing explanation. The method has the advantages of high calculation speed, good curve fitting effect and accurate explanation result.

Description

Method for solving set up by a kind of tight gas pressure break horizontal well numerical well testing model

Technical field

The present invention relates to a kind of tight gas pressure break horizontal well numerical well testing model and set up the method solved, belong to petroleum industry Oil/gas Well well testing field.

Background technology

Tight gas is as one of three large Unconventional forage, and stock number is enriched, and potentiality to be exploited is large.Tight gas reservoir has the features such as hypotonic, low pressure, low abundance, and gas well natural production ability is low, needs after reservoir reconstruction measure, just have commercial mining and is worth.Hydraulic fracturing technology and horizontal well technology improve the effective ways of tight gas reservoir production capacity.The exploitation of current tight gas reservoir generally adopts multistage pressure break horizontal well technology.

The well test analysis of tight gas reservoir pressure break horizontal well is the important means obtaining the rear fracture parameters of pressure and reservoir parameter, is also effective ways seepage flow mechanism being carried out to directly checking.Due to the complexity of tight gas reservoir multistage pressure break Horizontal Well Flow mechanism and well type, domestic at present also do not have the special WELL TEST INTERPRETATION MODEL for tight gas reservoir pressure break horizontal well, the multistage pressure break horizontal well analytic modell analytical model that main employing routine business software Saphir and EPS software provide makes an explanation analysis, and this badly influences the correct explanation of tight gas pressure break horizontal well well test data.

The main Problems existing of current tight gas pressure break horizontal well well test analysis:

1) Reservoir Seepage mechanism aspect, the WELL TEST INTERPRETATION MODEL adopted at present does not consider the non linear fluid flow through porous medium mechanism such as stress sensitive, free-boundary problem of tight gas reservoir.Conventional well test model based on Darcy linear seepage flow mechanism is not suitable for tight gas reservoir.If directly adopt conventional pressure break horizontal well test model to make an explanation analysiss to tight gas reservoir pressure break horizontal well well test data, will bring difficulty to the matching of well test analysis, also may there is very big error in the result that matching obtains.

2) reservoir aspect, conventional well testing analytic modell analytical model hypothesis reservoir is uniform dielectric, the non-homogeneity of reservoir cannot be considered, but actual reservoir has obvious non-homogeneity, reservoir will certainly produce material impact to tight gas seepage flow and bottom pressure response, does not therefore consider that reservoir heterogeneity also can bring certain influence to well test curve match and explanation results at present.

3) fractue spacing aspect, conventional pressure break horizontal well well testing analytic modell analytical model cannot process the situation of crack non-equidistance at present, and the position of the pressure break of actual conditions is nearly all unequal-interval, and this is the no small error to practical application band also.

4) pit shaft multiphase flow and well track aspect, nearly all do not consider the impact of pit shaft multiphase flow and well track at present in WELL TEST INTERPRETATION MODEL.The exploitation of tight gas is usually along with the output of water, and causing in pit shaft is biphase gas and liquid flow, and gas-liquid two-phase fails to be convened for lack of a quorum and strengthens the flow resistance of downhole well fluid, causes pressure reduction larger in pit shaft, thus has influence on well test analysis result.In addition, in current horizontal well well testing test pressure gauge usually under enter the position of 10 ~ 20m more than kickoff point (KOP), distance net horizontal section has the distance of more than 500m.And in well test analysis, usually the pressure measured by pressure gauge is used as the pressure of horizontal wellbore section, this will certainly cause the error of well test analysis.The best method addressed this problem is exactly set up the pressure break horizontal well test model considering pit shaft multiphase flow and real well track.

Summary of the invention

In order to overcome the shortcoming of above-mentioned prior art, the object of the present invention is to provide that a kind of computational speed is fast, curve good, explanation results accurately tight gas pressure break horizontal well numerical well testing model set up the method solved.

In order to achieve the above object, the technical scheme that the present invention takes is: method for solving set up by a kind of tight gas pressure break horizontal well numerical well testing model, comprises the following steps:

Step one: the plastid two-dimensionally of tight gas reservoir pressure break horizontal well and the generation of three-dimensional geologic;

Step 2: to the tight gas reservoir pressure break horizontal well generated two-dimensionally plastid and three-dimensional geologic to carry out grid discrete;

Step 3: pit shaft calculates without the flow model in porous media of pressure reduction;

Step 4: set up coupling model, and the coupling model set up is solved, and the solution obtained is generated Well Testing Theory curve;

Step 5: say that the theoretical curve that obtains in step 4 and measured curve carry out matching, obtain the parameter of well test analysis.

Generate plastid two-dimensionally and the three-dimensional geologic of tight gas reservoir pressure break horizontal well in described step one, concrete steps are as follows:

1) external boundary of geologic body residing for pressure break horizontal well, pit shaft inner boundary, crack and recombination region, then by arranging inner and outer boundary and the crack attribute concrete size of plastid and shape definitely, drawing and setting up plastid two-dimensionally;

2) according to the two-dimentional physique of foundation and the position of well track and reservoir up-and-down boundary, solid Boolean calculation is utilized to generate three-dimensional geologic.

The concrete steps that plastid two-dimensionally in described step 2 and three-dimensional geologic carry out grid discrete are as follows:

First, Netgen open source software is bundled into merit and compiles, build running environment;

Then, the inner and outer boundary in plastid two-dimensionally and three-dimensional geologic is belonged to and forms two-dimensional grid file and three-dimensional grid file according to the requirement of Netgen grid file form respectively, then it is discrete to carry out grid according to the grid discrete step that Netgen is arranged.

Described horizontal wellbore without the flow model in porous media calculating concrete grammar of pressure reduction is:

1) stratum and the fisstured flow equation of stress sensitive is considered

Stratum filtration equation:

$\begin{array}{c}{K}_{\mathrm{xD}}\frac{{\∂}^2{p}_{\mathrm{DR}}}{{\∂x}_D^2}+K_{\mathrm{yD}}\frac{{\∂}^2{p}_{\mathrm{DR}}}{{\∂y}_D^2}+K_{\mathrm{zD}}\frac{{\∂}^2{p}_{\mathrm{DR}}}{{\∂z}_D^2}-K_{\mathrm{xD}}{\mathrm{\γ}}_D{\left(\frac{\∂p_{\mathrm{DR}}}{{\∂x}_D}\right)}^2\\ -K_{\mathrm{yD}}{\mathrm{\γ}}_D{\left(\frac{{\∂p}_{\mathrm{DR}}}{{\∂y}_D}\right)}^2-K_{\mathrm{zD}}{\mathrm{\γ}}_D{\left(\frac{{\∂p}_{\mathrm{DR}}}{{\∂z}_D}\right)}^2=\frac{e^{{\mathrm{\γ}}_D{p}_{\mathrm{DR}}}}{K_{\mathrm{lD}}}\frac{{\∂p}_{\mathrm{DR}}}{{\∂t}_D}\end{array}---\left(1\right)$

Fisstured flow equation:

$\frac{{\∂}^2{p}_{\mathrm{Df}}}{{\∂r}_D^2}+K_{\mathrm{zD}}\frac{{\∂}^2{p}_{\mathrm{Df}}}{{\∂z}_D^2}=\frac{1}{K_{\mathrm{fD}}}\frac{{\∂p}_{\mathrm{Df}}}{{\∂t}_D}---\left(2\right)$

Primary condition:

p _D(x,y,z,0)＝0 (3)

Internal boundary condition:

$\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{A}_j{K}_{\mathrm{jD}}{e}^{-{\mathrm{\γ}}_D{p}_{\mathrm{jD}}}{\left(\frac{{\∂p}_{\mathrm{jD}}}{{\∂n}^{\′}}\right)|}_{{\mathrm{\Γ}}_{\mathrm{in}}}=2\mathrm{\π}{h}_D(1-C_{\mathrm{DL}}\frac{{\mathrm{dp}}_{\mathrm{wD}}}{{\mathrm{dt}}_D})---\left(4\right)$

$p_{\mathrm{wD}}\left(t_D\right)=p_{\mathrm{jD}}-S_t\frac{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{A}_{\mathrm{jD}}{K}_{\mathrm{jD}}{e}^{-{\mathrm{\γ}}_D{p}_{\mathrm{jD}}}}{2\mathrm{\π}{h}_D}\left(\frac{{\∂p}_{\mathrm{jD}}}{\∂n^{\′}}\right)-{\mathrm{MP}}_{\mathrm{jD}}---\left(5\right)$

Outer Boundary Conditions:

Close: $\frac{\∂p_D}{\∂n^{\′}}{|}_{{\mathrm{\Γ}}_{\mathrm{out}}}=0---\left(6\right)$

Level pressure: $p_D{|}_{{\mathrm{\Γ}}_{\mathrm{out}}}=0---\left(7\right)$ Symbol implication in formula:

P _dRfor the dimensionless pressure of subterranean formation zone; p _dffor the dimensionless pressure of crack area; t _dfor nondimensional time; C _dLfor dimensionless wellbore storage constant; K _xDfor x direction dimensionless permeability; K _yDfor y direction dimensionless permeability; K _zDfor z direction dimensionless permeability; K _fDfor dimensionless fracture permeabgility; γ _dfor dimensionless permeability modules; p _wDfor the dimensionless pressure of Article 1 crack and pit shaft point of intersection; MP _jDfor j point dimensionless pressure and p _wDbetween difference; A _jfor inner boundary triangle dimensionless area; h _dfor dimensionless reservoir thickness; S _tfor wellbore skin coefficient;

2) equation solution

First introduce conversion, by nonlinear filtration equation linearisation, then adopt Finite Element Method to solve, transformation for mula is:

$p_D=-\frac{1}{{\mathrm{\γ}}_D}\ln (1-{\mathrm{\γ}}_D\mathrm{\η})---\left(8\right)$

Adopt mixed finite element method to be solved in the stratum after conversion and fisstured flow equations simultaneousness, the finite element equation of stratum and Fracture System is decomposed into the finite element equation of subterranean formation zone (on the right of formula 9 Section 1) and represents the finite element equation (on the right of formula 9 Section 2) of Fracture System.

${\∫\∫\∫}_{\mathrm{\Ω}}\mathrm{FEQd\Ω}={\∫\∫\∫}_{{\mathrm{\Ω}}_m}\mathrm{FEQd}{\mathrm{\Ω}}_m+w_f\·{\∫\∫}_{\stackrel{\‾}{{\mathrm{\Ω}}_f}}\mathrm{FEQd}\stackrel{\‾}{{\mathrm{\Ω}}_f}---\left(9\right)$

A. subterranean formation zone three-dimensional finite element equation is:

$\begin{array}{c}{\mathrm{VK}}_{\mathrm{LD}}(K_{\mathrm{xD}}{b}_i^2+K_{\mathrm{yD}}{c}_i^2+K_{\mathrm{zD}}{d}_i^2+\frac{1}{10K_{\mathrm{LD}}\mathrm{\Δ}{t}_D}){\mathrm{\η}}_i^{e,n+1}+{\mathrm{VK}}_{\mathrm{LD}}\left(\begin{array}{c}{K}_{\mathrm{xD}}{b}_i{b}_j+K_{\mathrm{yD}}{c}_i{c}_j+K_{\mathrm{zD}}{d}_i{d}_j\\ +\frac{1}{20K_{\mathrm{LD}}\mathrm{\Δ}{t}_D}\end{array}\right){\mathrm{\η}}_j^{e,n+1}\\ +{\mathrm{VK}}_{\mathrm{LD}}(K_{\mathrm{xD}}{b}_i{b}_k+K_{\mathrm{yD}}{c}_i{c}_k+K_{\mathrm{zD}}{d}_i{d}_k+\frac{1}{20K_{\mathrm{LD}}\mathrm{\Δ}{t}_D}){\mathrm{\η}}_k^{e,n+1}+{\mathrm{VK}}_L\left(\begin{array}{c}{K}_{\mathrm{xD}}{b}_i{b}_m+K_{\mathrm{yD}}{c}_i{c}_m+\\ K_{\mathrm{zD}}{d}_i{d}_m+\frac{1}{20K_{\mathrm{LD}}\mathrm{\Δ}{t}_D}\end{array}\right){\mathrm{\η}}_m^{e,n+1}\\ -\frac{{\mathrm{AK}}_{\mathrm{LD}}}{6}\frac{\∂{\mathrm{\η}}_i^{e,n+1}}{\∂n^{\′}}-\frac{{\mathrm{AK}}_{\mathrm{LD}}}{12}\frac{\∂{\mathrm{\η}}_{j/k/m}^{e,n+1}}{\∂n^{\′}}-\frac{{\mathrm{AK}}_{\mathrm{LD}}}{12}\frac{\∂{\mathrm{\η}}_{k/m/j}^{e,n+1}}{\∂n^{\′}}\\ =\frac{V}{10\mathrm{\Δ}{t}_D}{\mathrm{\η}}_i^{e,n}+\frac{V}{20\mathrm{\Δ}{t}_D}{\mathrm{\η}}_j^{e,n}+\frac{V}{20\mathrm{\Δ}{t}_D}{\mathrm{\η}}_k^{e,n}+\frac{V}{20\mathrm{\Δ}{t}_D}{\mathrm{\η}}_m^{e,n}\end{array}---\left(10\right)$

$\frac{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{A}_{\mathrm{jD}}{K}_{\mathrm{jD}}\left(\frac{\∂{\mathrm{\η}}_{\mathrm{jD}}}{\∂n^{\′}}\right){|}_{{\mathrm{\Γ}}_{\mathrm{in}}}}{2\mathrm{\π}{h}_D}+\frac{C_{\mathrm{DL}}}{\mathrm{\Δ}{t}_D}{\mathrm{\η}}_w^{e,n+1}=1+\frac{C_{\mathrm{DL}}}{\mathrm{\Δ}{t}_D}{\mathrm{\η}}_w^{e,n}---\left(11\right)$

${\mathrm{\η}}_w\left(t_D\right)=\mathrm{\η}{|}_{{\mathrm{\Γ}}_{\mathrm{in}}}-S_t\frac{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{A}_{\mathrm{jD}}{K}_{\mathrm{jD}}}{2\mathrm{\π}{h}_D}\left(\frac{\∂{\mathrm{\η}}_j}{\∂n^{\′}}\right)---\left(12\right)$

B. fracture surface two dimensional finite element equation is:

$\begin{array}{c}A(K_{\mathrm{rD}}{b}_i^2+K_{\mathrm{zD}}{c}_i^2+\frac{1}{6\mathrm{\Δ}{t}_D}){\mathrm{\η}}_i^{e,n+1}+A(K_{\mathrm{rD}}{b}_i{b}_j+K_{\mathrm{zD}}{c}_i{c}_j+\frac{1}{12\mathrm{\Δ}{t}_D}){\mathrm{\η}}_j^{e,n+1}\\ +A(K_{\mathrm{rD}}{b}_i{b}_k+K_{\mathrm{zD}}{c}_i{c}_k+\frac{1}{12\mathrm{\Δ}{t}_D}){\mathrm{\η}}_k^{e,n+1}-\frac{L}{3}\frac{\∂{\mathrm{\η}}_i^{e,n+1}}{\∂n^{\′}}-\frac{L}{6}\frac{\∂{\mathrm{\η}}_{j,k}^{e,n+1}}{\∂n^{\′}}\\ =\frac{A}{6\mathrm{\Δ}{t}_D}{\mathrm{\η}}_i^{e,n}+\frac{A}{12\mathrm{\Δ}{t}_D}{\mathrm{\η}}_j^{e,n}+\frac{A}{12\mathrm{\Δ}{t}_D}{\mathrm{\η}}_k^{e,n}\end{array}---\left(13\right)$

$\frac{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{L}_{\mathrm{jD}}{K}_{\mathrm{jD}}\left(\frac{\∂{\mathrm{\η}}_{\mathrm{jD}}}{\∂n^{\′}}\right)}{2\mathrm{\π}}+\frac{C_{\mathrm{DL}}}{\mathrm{\Δ}{t}_D}{\mathrm{\η}}_w^{e,n+1}=1+\frac{C_{\mathrm{DL}}}{\mathrm{\Δ}{t}_D}{\mathrm{\η}}_w^{e,n}---\left(14\right)$

${\mathrm{\η}}_w\left(t_D\right)={\mathrm{\η}}_j-S_t\frac{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{L}_{\mathrm{jD}}{K}_{\mathrm{LD}}}{2\mathrm{\π}{h}_D}\left(\frac{{\∂\mathrm{\η}}_j}{{\∂n}^{\′}}\right)---\left(15\right)$

By finite element equation (10) ~ (15) simultaneous composition system stiffness matrix, the SuperLU numerical solution device of parallelization is utilized to solve large linear systems, pressure field distribution and the inner boundary normal pressure gradient of whole reservoir can be obtained, then calculate each crack production flow thus:

$Q_{\mathrm{fi}}=Q_{\mathrm{sc}}\frac{w_{\mathrm{fD}}{L}_{\mathrm{iD}}{K}_{\mathrm{fD}}{e}^{-{\mathrm{\γ}}_D{p}_{\mathrm{iD}}}}{2\mathrm{\π}{h}_D}\left(\frac{{\∂p}_{\mathrm{iD}}}{{\∂n}^{\′}}\right){|}_{{\mathrm{\Γ}}_{\mathrm{in}}}---\left(16\right)$

Symbol implication in formula:

η is transformation to linearity parameter; η _wfor the conversion parameter corresponding to dimensionless bottom pressure value; w _ffor crack width, m; w _fDfor dimensionless crack width; L _jDfor crack inner boundary unit wires length; V is tetrahedron volume; B, c, d are finite element coefficient; I, j, k, m are finite element tetrahedron four summit sequence numbers; Q _fibe the flow of the i-th crack, m ³/ d; Q _scfor the flow of gas well under mark condition, m ³/ d.

Set up coupling model in described step 4, and the coupling model set up is solved carry out according to following steps:

A, according to the wellbore pressure p that step 4 calculates _wDwith crack flow Q _fi, adopt pit shaft multiphase flow design formulas to calculate, obtain the dimensionless pressure reduction MP of pit shaft each point _iD, concrete formula is as follows:

The basic equation that pit shaft multiphase flow calculates is:

$-\frac{{\mathrm{dp}}_w}{\mathrm{dZ}}=\frac{({\mathrm{\ρ}}_l{H}_L+{\mathrm{\ρ}}_g(1-H_L))\sin \mathrm{\θ}+\frac{{\mathrm{\λGv}}_m}{2\mathrm{AD}}}{1-\left\{\right[{\mathrm{\ρ}}_l{H}_L+{\mathrm{\ρ}}_g(1-H_L)\left]v_m{v}_{\mathrm{sg}}\right\}/p}---\left(17\right)$

Symbol implication in formula:

ρ _lfor fluid density, kg/m ³; ρ _gfor gas density, kg/m ³; G is gas-liquid mixture mass flow, kg/s; v _mfor mixture flowing velocity, m/s; v _sgfor gas superficial flow velocity, m/s; A is pit shaft oil pipe sectional area, m ²; D is pipe aperture, m.

Wherein liquid holdup H _ladopt Beggs-Brill method to calculate with coefficient of frictional resistance λ, obtain wellbore pressure gradient according to formula (17), then to add up the pressure difference obtained between each point and shaft bottom standard point according to pit shaft inner boundary relative distance:

${\mathrm{MP}}_i=\underset{j=0}{\overset{i}{\mathrm{\Σ}}}{(-\frac{{\mathrm{dp}}_w}{\mathrm{dZ}})}_i{\mathrm{\ΔZ}}_j---\left(18\right)$

By each point pressure reduction MP obtained _icarry out nondimensionalization, each point dimensionless pressure reduction MP can be obtained _iD.

B, the inner boundary each point pressure reduction MP obtained in steps A _iD, be brought in flow model in porous media and calculate wellbore pressure p ' _wDwith crack flow Q ' _fi, the coupling condition of pit shaft multiphase flow model and flow model in porous media is as follows:

On stratum and pit shaft interface:

$p_{\mathrm{DRi}}{|}_{{\mathrm{\Γ}}_{\mathrm{in}}}=p_{\mathrm{wD}}+{\mathrm{MP}}_{\mathrm{iD}}---\left(19\right)$

On crack and pit shaft interface:

$p_{\mathrm{Dfj}}{|}_{{\mathrm{\Γ}}_{\mathrm{in}}}=p_{\mathrm{wD}}+{\mathrm{MP}}_{\mathrm{jD}}---\left(20\right)$

C, the p ' that front and back iteration step calculates _wDand p _wDsubtract each other, when poor absolute value is less than ε, then continue the calculating of future time step, when the absolute value of both differences is more than or equal to ε, then the current p ' obtained between the two _wDand Q ' _finew wellbore pressure pWD and crack flow Q is obtained after being averaged with the value of back _fi, and be brought in steps A and carry out iteration, until the p ' obtained _wDwith p _wDthe absolute value carrying out the difference of subtracting each other is less than ε, and the general value of ε is 10 ^-4;

D, the p ' that step C determines _wDwith the Q ' of correspondence _ficarry out record, if now time step k< total time walks n, then continue the calculating of future time step, according to current wellbore pressure p _wD, crack flow Q _fi, from steps A, calculate the p ' of a new round _wDand Q ' _fi;

E, obtains p ' institute in n calculating in step D _wDand Q ' _fivalue generate Well Testing Theory curve.

Described step 5 mainly theoretical curve and measured curve contrasts matching in log-log graph, semilogarithmic plot and full historical pressures curve map, according to the fitting degree of curve, the adjustment of fitting parameter can be carried out, finally make theoretical curve and measured curve can both be coincide preferably in log-log graph, semilogarithmic plot and full historical pressures curve map, curve just completes; After curve completes, namely can obtain the parameter of well test analysis, comprise fracture parameters and reservoir parameter.

The present invention adopts above technical scheme, has the following advantages, and has that computational speed is fast, curve good, explanation results is accurate.

Accompanying drawing explanation

Fig. 1 consider the pressure break horizontal well three-dimensional geologic of well track and reservoir heterogeneity and grid discrete;

The flow chart that Fig. 2 coupling model is set up and solved.

Detailed description of the invention

Below in conjunction with drawings and Examples, the present invention is described in detail.

Embodiment 1

A method for solving set up by tight gas pressure break horizontal well numerical well testing model as shown in Figure 2, comprises the following steps:

Step one: the plastid two-dimensionally of tight gas reservoir pressure break horizontal well and the generation of three-dimensional geologic; Mainly

1) external boundary of geologic body residing for pressure break horizontal well, pit shaft inner boundary, crack and recombination region, then by arranging inner and outer boundary and the crack attribute concrete size of plastid and shape definitely, drawing and setting up plastid two-dimensionally;

2) according to the two-dimentional physique of foundation and the position of well track and reservoir up-and-down boundary, solid Boolean calculation is utilized to generate three-dimensional geologic.

Step 2: as shown in Figure 1 to the tight gas reservoir pressure break horizontal well generated two-dimensionally plastid and three-dimensional geologic to carry out grid discrete; First, Netgen open source software is bundled into merit and compiles, build running environment;

Then, the inner and outer boundary in plastid two-dimensionally and three-dimensional geologic is belonged to and forms two-dimensional grid file and three-dimensional grid file according to the requirement of Netgen grid file form respectively, then it is discrete to carry out grid according to the grid discrete step that Netgen is arranged.

Step 3: horizontal wellbore calculates without the flow model in porous media of pressure reduction;

1) stratum and the fisstured flow equation of stress sensitive is considered

Stratum filtration equation:

$\begin{array}{c}{K}_{\mathrm{xD}}\frac{{\&PartialD;}^2{p}_{\mathrm{DR}}}{{\&PartialD;x}_D^2}+K_{\mathrm{yD}}\frac{{\&PartialD;}^2{p}_{\mathrm{DR}}}{{\&PartialD;y}_D^2}+K_{\mathrm{zD}}\frac{{\&PartialD;}^2{p}_{\mathrm{DR}}}{{\&PartialD;z}_D^2}-K_{\mathrm{xD}}{\mathrm{\&gamma;}}_D{\left(\frac{\&PartialD;p_{\mathrm{DR}}}{{\&PartialD;x}_D}\right)}^2\\ -K_{\mathrm{yD}}{\mathrm{\&gamma;}}_D{\left(\frac{{\&PartialD;p}_{\mathrm{DR}}}{{\&PartialD;y}_D}\right)}^2-K_{\mathrm{zD}}{\mathrm{\&gamma;}}_D{\left(\frac{{\&PartialD;p}_{\mathrm{DR}}}{{\&PartialD;z}_D}\right)}^2=\frac{e^{{\mathrm{\&gamma;}}_D{p}_{\mathrm{DR}}}}{K_{\mathrm{lD}}}\frac{{\&PartialD;p}_{\mathrm{DR}}}{{\&PartialD;t}_D}\end{array}---\left(1\right)$

Fisstured flow equation:

$\frac{{\&PartialD;}^2{p}_{\mathrm{Df}}}{{\&PartialD;r}_D^2}+K_{\mathrm{zD}}\frac{{\&PartialD;}^2{p}_{\mathrm{Df}}}{{\&PartialD;z}_D^2}=\frac{1}{K_{\mathrm{fD}}}\frac{{\&PartialD;p}_{\mathrm{Df}}}{{\&PartialD;t}_D}---\left(2\right)$

Primary condition:

p _D(x,y,z,0)＝0 (3)

Internal boundary condition:

$\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{A}_j{K}_{\mathrm{jD}}{e}^{-{\mathrm{\&gamma;}}_D{p}_{\mathrm{jD}}}{\left(\frac{{\&PartialD;p}_{\mathrm{jD}}}{{\&PartialD;n}^{\&prime;}}\right)|}_{{\mathrm{\&Gamma;}}_{\mathrm{in}}}=2\mathrm{\&pi;}{h}_D(1-C_{\mathrm{DL}}\frac{{\mathrm{dp}}_{\mathrm{wD}}}{{\mathrm{dt}}_D})---\left(4\right)$

$p_{\mathrm{wD}}\left(t_D\right)=p_{\mathrm{jD}}-S_t\frac{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{A}_{\mathrm{jD}}{K}_{\mathrm{jD}}{e}^{-{\mathrm{\&gamma;}}_D{p}_{\mathrm{jD}}}}{2\mathrm{\&pi;}{h}_D}\left(\frac{{\&PartialD;p}_{\mathrm{jD}}}{\&PartialD;n^{\&prime;}}\right)-{\mathrm{MP}}_{\mathrm{jD}}---\left(5\right)$

Outer Boundary Conditions:

Close: $\frac{\&PartialD;p_D}{\&PartialD;n^{\&prime;}}{|}_{{\mathrm{\&Gamma;}}_{\mathrm{out}}}=0---\left(6\right)$

Level pressure: $p_D{|}_{{\mathrm{\&Gamma;}}_{\mathrm{out}}}=0---\left(7\right)$ Symbol implication in formula:

P _dRfor the dimensionless pressure of subterranean formation zone; p _dffor the dimensionless pressure of crack area; t _dfor nondimensional time; C _dLfor dimensionless wellbore storage constant; K _xDfor x direction dimensionless permeability; K _yDfor y direction dimensionless permeability; K _zDfor z direction dimensionless permeability; K _fDfor dimensionless fracture permeabgility; γ _dfor dimensionless permeability modules; p _wDfor the dimensionless pressure of Article 1 crack and pit shaft point of intersection; MP _jDfor j point dimensionless pressure and p _wDbetween difference; A _jfor inner boundary triangle dimensionless area; h _dfor dimensionless reservoir thickness; S _tfor wellbore skin coefficient;

2) equation solution

First introduce conversion, by nonlinear filtration equation linearisation, then adopt Finite Element Method to solve, transformation for mula is:

$p_D=-\frac{1}{{\mathrm{\&gamma;}}_D}\ln (1-{\mathrm{\&gamma;}}_D\mathrm{\&eta;})---\left(8\right)$

Adopt mixed finite element method to be solved in the stratum after conversion and fisstured flow equations simultaneousness, the finite element equation of stratum and Fracture System is decomposed into the finite element equation of subterranean formation zone (on the right of formula 9 Section 1) and represents the finite element equation (on the right of formula 9 Section 2) of Fracture System.

${\&Integral;\&Integral;\&Integral;}_{\mathrm{\&Omega;}}\mathrm{FEQd\&Omega;}={\&Integral;\&Integral;\&Integral;}_{{\mathrm{\&Omega;}}_m}\mathrm{FEQd}{\mathrm{\&Omega;}}_m+w_f\&CenterDot;{\&Integral;\&Integral;}_{\stackrel{\&OverBar;}{{\mathrm{\&Omega;}}_f}}\mathrm{FEQd}\stackrel{\&OverBar;}{{\mathrm{\&Omega;}}_f}---\left(9\right)$

A. subterranean formation zone three-dimensional finite element equation is:

$\begin{array}{c}{\mathrm{VK}}_{\mathrm{LD}}(K_{\mathrm{xD}}{b}_i^2+K_{\mathrm{yD}}{c}_i^2+K_{\mathrm{zD}}{d}_i^2+\frac{1}{10K_{\mathrm{LD}}\mathrm{\&Delta;}{t}_D}){\mathrm{\&eta;}}_i^{e,n+1}+{\mathrm{VK}}_{\mathrm{LD}}\left(\begin{array}{c}{K}_{\mathrm{xD}}{b}_i{b}_j+K_{\mathrm{yD}}{c}_i{c}_j+K_{\mathrm{zD}}{d}_i{d}_j\\ +\frac{1}{20K_{\mathrm{LD}}\mathrm{\&Delta;}{t}_D}\end{array}\right){\mathrm{\&eta;}}_j^{e,n+1}\\ +{\mathrm{VK}}_{\mathrm{LD}}(K_{\mathrm{xD}}{b}_i{b}_k+K_{\mathrm{yD}}{c}_i{c}_k+K_{\mathrm{zD}}{d}_i{d}_k+\frac{1}{20K_{\mathrm{LD}}\mathrm{\&Delta;}{t}_D}){\mathrm{\&eta;}}_k^{e,n+1}+{\mathrm{VK}}_L\left(\begin{array}{c}{K}_{\mathrm{xD}}{b}_i{b}_m+K_{\mathrm{yD}}{c}_i{c}_m+\\ K_{\mathrm{zD}}{d}_i{d}_m+\frac{1}{20K_{\mathrm{LD}}\mathrm{\&Delta;}{t}_D}\end{array}\right){\mathrm{\&eta;}}_m^{e,n+1}\\ -\frac{{\mathrm{AK}}_{\mathrm{LD}}}{6}\frac{\&PartialD;{\mathrm{\&eta;}}_i^{e,n+1}}{\&PartialD;n^{\&prime;}}-\frac{{\mathrm{AK}}_{\mathrm{LD}}}{12}\frac{\&PartialD;{\mathrm{\&eta;}}_{j/k/m}^{e,n+1}}{\&PartialD;n^{\&prime;}}-\frac{{\mathrm{AK}}_{\mathrm{LD}}}{12}\frac{\&PartialD;{\mathrm{\&eta;}}_{k/m/j}^{e,n+1}}{\&PartialD;n^{\&prime;}}\\ =\frac{V}{10\mathrm{\&Delta;}{t}_D}{\mathrm{\&eta;}}_i^{e,n}+\frac{V}{20\mathrm{\&Delta;}{t}_D}{\mathrm{\&eta;}}_j^{e,n}+\frac{V}{20\mathrm{\&Delta;}{t}_D}{\mathrm{\&eta;}}_k^{e,n}+\frac{V}{20\mathrm{\&Delta;}{t}_D}{\mathrm{\&eta;}}_m^{e,n}\end{array}---\left(10\right)$

$\frac{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{A}_{\mathrm{jD}}{K}_{\mathrm{jD}}\left(\frac{\&PartialD;{\mathrm{\&eta;}}_{\mathrm{jD}}}{\&PartialD;n^{\&prime;}}\right){|}_{{\mathrm{\&Gamma;}}_{\mathrm{in}}}}{2\mathrm{\&pi;}{h}_D}+\frac{C_{\mathrm{DL}}}{\mathrm{\&Delta;}{t}_D}{\mathrm{\&eta;}}_w^{e,n+1}=1+\frac{C_{\mathrm{DL}}}{\mathrm{\&Delta;}{t}_D}{\mathrm{\&eta;}}_w^{e,n}---\left(11\right)$

${\mathrm{\&eta;}}_w\left(t_D\right)=\mathrm{\&eta;}{|}_{{\mathrm{\&Gamma;}}_{\mathrm{in}}}-S_t\frac{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{A}_{\mathrm{jD}}{K}_{\mathrm{jD}}}{2\mathrm{\&pi;}{h}_D}\left(\frac{\&PartialD;{\mathrm{\&eta;}}_j}{\&PartialD;n^{\&prime;}}\right)---\left(12\right)$

B. fracture surface two dimensional finite element equation is:

$\begin{array}{c}A(K_{\mathrm{rD}}{b}_i^2+K_{\mathrm{zD}}{c}_i^2+\frac{1}{6\mathrm{\&Delta;}{t}_D}){\mathrm{\&eta;}}_i^{e,n+1}+A(K_{\mathrm{rD}}{b}_i{b}_j+K_{\mathrm{zD}}{c}_i{c}_j+\frac{1}{12\mathrm{\&Delta;}{t}_D}){\mathrm{\&eta;}}_j^{e,n+1}\\ +A(K_{\mathrm{rD}}{b}_i{b}_k+K_{\mathrm{zD}}{c}_i{c}_k+\frac{1}{12\mathrm{\&Delta;}{t}_D}){\mathrm{\&eta;}}_k^{e,n+1}-\frac{L}{3}\frac{\&PartialD;{\mathrm{\&eta;}}_i^{e,n+1}}{\&PartialD;n^{\&prime;}}-\frac{L}{6}\frac{\&PartialD;{\mathrm{\&eta;}}_{j,k}^{e,n+1}}{\&PartialD;n^{\&prime;}}\\ =\frac{A}{6\mathrm{\&Delta;}{t}_D}{\mathrm{\&eta;}}_i^{e,n}+\frac{A}{12\mathrm{\&Delta;}{t}_D}{\mathrm{\&eta;}}_j^{e,n}+\frac{A}{12\mathrm{\&Delta;}{t}_D}{\mathrm{\&eta;}}_k^{e,n}\end{array}---\left(13\right)$

$\frac{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{L}_{\mathrm{jD}}{K}_{\mathrm{jD}}\left(\frac{\&PartialD;{\mathrm{\&eta;}}_{\mathrm{jD}}}{\&PartialD;n^{\&prime;}}\right)}{2\mathrm{\&pi;}}+\frac{C_{\mathrm{DL}}}{\mathrm{\&Delta;}{t}_D}{\mathrm{\&eta;}}_w^{e,n+1}=1+\frac{C_{\mathrm{DL}}}{\mathrm{\&Delta;}{t}_D}{\mathrm{\&eta;}}_w^{e,n}---\left(14\right)$

${\mathrm{\&eta;}}_w\left(t_D\right)={\mathrm{\&eta;}}_j-S_t\frac{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{L}_{\mathrm{jD}}{K}_{\mathrm{LD}}}{2\mathrm{\&pi;}{h}_D}\left(\frac{{\&PartialD;\mathrm{\&eta;}}_j}{{\&PartialD;n}^{\&prime;}}\right)---\left(15\right)$

By finite element equation (10) ~ (15) simultaneous composition system stiffness matrix, the SuperLU numerical solution device of parallelization is utilized to solve large linear systems, pressure field distribution and the inner boundary normal pressure gradient of whole reservoir can be obtained, then calculate each crack production flow thus:

$Q_{\mathrm{fi}}=Q_{\mathrm{sc}}\frac{w_{\mathrm{fD}}{L}_{\mathrm{iD}}{K}_{\mathrm{fD}}{e}^{-{\mathrm{\&gamma;}}_D{p}_{\mathrm{iD}}}}{2\mathrm{\&pi;}{h}_D}\left(\frac{{\&PartialD;p}_{\mathrm{iD}}}{{\&PartialD;n}^{\&prime;}}\right){|}_{{\mathrm{\&Gamma;}}_{\mathrm{in}}}---\left(16\right)$

Symbol implication in formula:

η is transformation to linearity parameter; η _wfor the conversion parameter corresponding to dimensionless bottom pressure value; w _ffor crack width, m; w _fDfor dimensionless crack width; L _jDfor crack inner boundary unit wires length; V is tetrahedron volume; B, c, d are finite element coefficient; I, j, k, m are finite element tetrahedron four summit sequence numbers; Q _fibe the flow of the i-th crack, m ³/ d; Q _scfor the flow of gas well under mark condition, m ³/ d.

Step 4: set up coupling model, and the coupling model set up is solved, and the solution obtained is generated Well Testing Theory curve; And the coupling model set up is solved to carry out according to following steps:

A, according to the wellbore pressure p that step 4 calculates _wDwith crack flow Q _fi, adopt pit shaft multiphase flow design formulas to calculate, obtain the dimensionless pressure reduction MP of pit shaft each point _iD, concrete formula is as follows:

The basic equation that pit shaft multiphase flow calculates is:

$-\frac{{\mathrm{dp}}_w}{\mathrm{dZ}}=\frac{({\mathrm{\&rho;}}_l{H}_L+{\mathrm{\&rho;}}_g(1-H_L))\sin \mathrm{\&theta;}+\frac{{\mathrm{\&lambda;Gv}}_m}{2\mathrm{AD}}}{1-\left\{\right[{\mathrm{\&rho;}}_l{H}_L+{\mathrm{\&rho;}}_g(1-H_L)\left]v_m{v}_{\mathrm{sg}}\right\}/p}---\left(17\right)$

Symbol implication in formula:

ρ _lfor fluid density, kg/m ³; ρ _gfor gas density, kg/m ³; G is gas-liquid mixture mass flow, kg/s; v _mfor mixture flowing velocity, m/s; v _sgfor gas superficial flow velocity, m/s; A is pit shaft oil pipe sectional area, m ²; D is pipe aperture, m.

Wherein liquid holdup H _ladopt Beggs-Brill method to calculate with coefficient of frictional resistance λ, obtain wellbore pressure gradient according to formula (17), then to add up the pressure difference obtained between each point and shaft bottom standard point according to pit shaft inner boundary relative distance:

${\mathrm{MP}}_i=\underset{j=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-\frac{{\mathrm{dp}}_w}{\mathrm{dZ}})}_i{\mathrm{\&Delta;Z}}_j---\left(18\right)$

By each point pressure reduction MP obtained _icarry out nondimensionalization, each point dimensionless pressure reduction MP can be obtained _iD.

B, the inner boundary each point pressure reduction MP obtained in steps A _iD, be brought in flow model in porous media and calculate wellbore pressure p ' _wDwith crack flow Q ' _fi, the coupling condition of pit shaft multiphase flow model and flow model in porous media is as follows:

On stratum and pit shaft interface:

$p_{\mathrm{DRi}}{|}_{{\mathrm{\&Gamma;}}_{\mathrm{in}}}=p_{\mathrm{wD}}+{\mathrm{MP}}_{\mathrm{iD}}---\left(19\right)$

On crack and pit shaft interface:

$p_{\mathrm{Dfj}}{|}_{{\mathrm{\&Gamma;}}_{\mathrm{in}}}=p_{\mathrm{wD}}+{\mathrm{MP}}_{\mathrm{jD}}---\left(20\right)$

C, the p ' that front and back iteration step calculates _wDand p _wDsubtract each other, when poor absolute value is less than ε, then continue the calculating of future time step, when the absolute value of both differences is more than or equal to ε, then the current p ' obtained between the two _wDand Q ' _finew wellbore pressure p is obtained after being averaged with the value of back _wDwith crack flow Q _fi, and be brought in steps A and carry out iteration, until the p ' obtained _wDwith p _wDthe absolute value carrying out the difference of subtracting each other is less than ε, and the general value of ε is 10 ^-4;

D, the p ' that step C determines _wDwith the Q ' of correspondence _ficarry out record, if now time step k< total time walks n, then continue the calculating of future time step, according to current wellbore pressure p _wD, crack flow Q _fi, from steps A, calculate the p ' of a new round _wDand Q ' _fi;

E, obtains p ' institute in n calculating in step D _wDand Q ' _fivalue generate Well Testing Theory curve.Step 5: say that the theoretical curve that obtains in step 4 and measured curve carry out matching, obtain the parameter of well test analysis;

Described step 5 mainly theoretical curve and measured curve contrasts matching in log-log graph, semilogarithmic plot and full historical pressures curve map, according to the fitting degree of curve, the adjustment of fitting parameter can be carried out, finally make theoretical curve and measured curve can both be coincide preferably in log-log graph, semilogarithmic plot and full historical pressures curve map, curve just completes; After curve completes, namely can obtain the parameter of well test analysis, comprise fracture parameters and reservoir parameter.

Embodiment 2

1. realize tight gas reservoir pressure break horizontal well geologic body to generate fast

First set up plastid two-dimensionally, then be converted to three-dimensional geologic by plastid two-dimensionally.By writing the code that software geometric figure is drawn, realize the external boundary drawing geologic body residing for pressure break horizontal well, the functions such as pit shaft inner boundary, crack and recombination region, again by arranging inner and outer boundary and the crack attribute concrete size of plastid and shape definitely, plastid two-dimensionally can be set up fast in this way.By plastid two-dimensionally, in conjunction with the position of well track and reservoir up-and-down boundary, utilize solid Boolean calculation directly can generate three-dimensional geologic, adopt OpenCasCade instrument to carry out 3-D view display, three-dimensional geologic as shown in Figure 1.

2. realize plastid and the automatic discrete functionality of three-dimensional geological volume mesh two-dimensionally

First Netgen open source software is bundled into merit to compile, build running environment, then the inner and outer boundary attribute in plastid is two-dimensionally formed two-dimensional grid file according to the requirement of Netgen grid file form, then it is discrete to carry out grid according to the grid discrete step that Netgen is arranged.The departure process of three-dimensional geologic is similar with the departure process of plastid two-dimensionally, unlike three-dimensional grid discrete first need, three-dimensional geologic is exported as common format (as STEP form), then it is discrete to adopt Netgen to carry out grid to the file of this output format.The discrete time of plastid is shorter two-dimensionally, and required time is generally between 10s ~ 30s, and the discrete time of three-dimensional geologic is longer, and need determine according to the discrete density of grid, usual discrete time is no more than 5 minutes.The mesh node that the discrete rear automatic display of grid is discrete and grid line, as shown in Figure 1.

3. set up WELL TEST INTERPRETATION MODEL, the type of partitioning model parameter

According to the type of flow involved by pressure break horizontal well, set up the WELL TEST INTERPRETATION MODEL considering various factors coupling.Parameter involved in model is divided into known parameters and unknown parameter (parameter to be explained), and known parameters should arrange different input interfaces from unknown parameter, in order to avoid obscure.Known parameters inputs according to actual conditions, and unknown parameter can be determined after curve completes.

4. the code of compiling model numerical solution, implementation model rapid solving

After having set up WELL TEST INTERPRETATION MODEL, according to model feature, the numerical solution algorithm designed a model, the computer code of compiling model numerical solution, model solution flow chart as shown in Figure 2.First input basic data, set up model of geological structure body, it is discrete model of geological structure body to be carried out grid.Initial step carries out flow model in porous media calculating according to pit shaft without the mode of pressure reduction, calculate initial flowing bottomhole pressure (FBHP) and each crack flow distribution, pit shaft multiphase flow computation model is utilized to calculate the Pressure difference distribution of pit shaft according to this result, pressure reduction carries out flow model in porous media calculating thus again, obtain new bottom pressure and crack flow distribution, contrast with result of calculation before the iteration judging this time step and whether terminate.If both differences are less than ε in a small amount, carry out the calculating of next time step, otherwise then continue iteration, until carry out the calculating of future time step again after stable.Utilize the method finally to complete the calculating of all setting-up times step, namely complete solving coupling model.

5. carry out measured curve and theoretical curve matching

After model solution completes, automatically generate Well Testing Theory curve, theoretical curve and measured curve are contrasted matching in log-log graph, semilogarithmic plot and full historical pressures curve map.Log-log graph is in front view, and semilogarithmic plot and full historical pressures curve map are in auxiliary view district, according to the fitting degree of curve, can carry out the adjustment of fitting parameter.Finally make theoretical curve and measured curve can both be coincide preferably in log-log graph, semilogarithmic plot and full historical pressures curve map, curve just completes, after curve completes, namely can obtain the parameter of well test analysis, comprise fracture parameters and reservoir parameter.

Embodiment 3

1) tight gas reservoir pressure break horizontal well two and three dimensions geologic body generates and display methods fast;

First set up plastid two-dimensionally, then be converted to three-dimensional geologic by plastid two-dimensionally.The Heterogeneous Characteristics of reservoir is set according to Horizontal Well Log Interpretation result, comprise the distribution of degree of porosity, gas saturation and original permeability, different regions is divided in geologic body, each Regional Representative isotropic body, each isotropic body has different degree of porosity, water saturation and original permeability, is made up of the non-homogeneity of whole geological system many different isotropic bodies.Isotropic body number is more, and parameter differences is larger, then the non-homogeneity of geologic body is stronger; By plastid two-dimensionally, then in conjunction with the position of well track and reservoir up-and-down boundary, directly can generate three-dimensional geologic; In order to carry out the display of three-dimensional geologic better, carry out the comprehensive display of 3-D view based on OpenCasCade instrument, displaying contents comprises three-dimensional geologic, the discrete figure of grid, calculates cloud atlas etc.

2) realization that tight gas reservoir pressure break horizontal well grid is automatically discrete;

The geologic body of generation is converted to STEP formatted file, IGES formatted file or BREP formatted file, then Netgen Open-Source Tools is adopted to realize the automatic discrete functionality of grid of plastid and three-dimensional geologic two-dimensionally, two-dimensional grid is triangular mesh, and three-dimensional grid is tetrahedral grid.For the feature of pressure break horizontal well well testing problem, the quantity by arranging crack and the discrete node of pit shaft in grid departure process carrys out the grid near automatic infill well cylinder and crack, carrys out the discrete quality of control mesh by arranging grid discrete parameter.Grid is more intensive, and grid amount is more, and grid discrete time is longer, and calculate required time also longer, and the too small meeting of grid amount has influence on computational accuracy, therefore actual conditions should select moderate mesh-density.Two-dimensional grid amount is between 10,000 ~ 30,000 under normal circumstances, and three-dimensional grid amount is between 100,000 ~ 500,000.

3) foundation of tight gas reservoir pressure break horizontal well stratum filtration and pit shaft multiphase flow coupling model;

Tight gas reservoir pressure break horizontal well is divided into three flow regions: pit shaft, crack and stratum.The flow mechanism of zones of different is different, and pit shaft is biphase gas and liquid flow, and crack is high speed non-darcy flow, and stratum is the non linear fluid flow through porous medium with stress sensitive.First based on tight gas pressure break Horizontal Well Flow feature, crack district is seen as Thief zone high speed non-darcy flow region, because crack width is very little, the flowing of crack district fluid is considered as two-dimensional flow, sets up a kind of pressure break horizontal well two Dimension Numerical Value well test model and three-dimensional numerical value well test model; Two dimensional model is mainly used in horizontal wellbore with the situation of barrel forms completion, and threedimensional model is used for horizontal wellbore with the situation of bore hole mode completion; Pit shaft biphase gas and liquid flow is considered on the basis of this flow model in porous media again, adopts the flow resistance computational methods of Beggs-Brill method establishment pit shaft biphase gas and liquid flow; Flow model in porous media is coupled with pressure continuous print mode on pit shaft contact surface according to stratum with pit shaft multiphase flow model.

4) fast solution method of tight gas reservoir pressure break horizontal well coupling model.

Adopt mixed finite element method to carry out numerical solution to model, flow model in porous media is carried out being coupled iterative computation with pit shaft gas-liquid two-phase flow calculation model.According to the feature of well testing flow model in porous media, be log series model by Node configuration computing time, arrange 10 ~ 20 calculation levels in each logarithm period, the initial time of calculating and termination time can adjust according to actual conditions.The condition that each time step iteration terminates is that result of calculation not only meets flow model in porous media but also meet pit shaft multiphase flow computation model.The system of linear equations of SuperLU solver to two dimensional model is utilized to solve, the system of linear equations of SuperLU solver to threedimensional model of parallelization is utilized to solve, the Thread Count that parallel computation adopts is determined according to the grid amount calculated, and parallel line number of passes is set to 6 ~ 10 and can meets computation requirement under normal circumstances.

5. carry out measured curve and theoretical curve matching

After model solution completes, automatically generate Well Testing Theory curve, theoretical curve and measured curve are contrasted matching in log-log graph, semilogarithmic plot and full historical pressures curve map.Log-log graph is in front view, and semilogarithmic plot and full historical pressures curve map are in auxiliary view district, according to the fitting degree of curve, can carry out the adjustment of fitting parameter.Finally make theoretical curve and measured curve can both be coincide preferably in log-log graph, semilogarithmic plot and full historical pressures curve map, curve just completes, after curve completes, namely can obtain the parameter of well test analysis, comprise fracture parameters and reservoir parameter.

More than exemplifying is only illustrate of the present invention, does not form the restriction to protection scope of the present invention, everyly all belongs within protection scope of the present invention with the same or analogous design of the present invention.

Claims (6)

1. a method for solving set up by tight gas pressure break horizontal well numerical well testing model, it is characterized in that, comprises the following steps:

Step one: the plastid two-dimensionally of tight gas reservoir pressure break horizontal well and the generation of three-dimensional geologic;

Step 2: to the tight gas reservoir pressure break horizontal well generated two-dimensionally plastid and three-dimensional geologic to carry out grid discrete;

Step 3: pit shaft calculates without the flow model in porous media of pressure reduction;

Step 4: set up coupling model, and the coupling model set up is solved, and the solution obtained is generated Well Testing Theory curve;

Step 5: say that the theoretical curve that obtains in step 4 and measured curve carry out matching, obtain the parameter of well test analysis.

2. method for solving set up by a kind of tight gas pressure break horizontal well numerical well testing model according to claim 1, and it is characterized in that, generate plastid two-dimensionally and the three-dimensional geologic of tight gas reservoir pressure break horizontal well in described step one, concrete steps are as follows:

1) external boundary of geologic body residing for pressure break horizontal well, pit shaft inner boundary, crack and recombination region, then by arranging inner and outer boundary and the crack attribute concrete size of plastid and shape definitely, drawing and setting up plastid two-dimensionally;

2) according to the two-dimentional physique of foundation and the position of well track and reservoir up-and-down boundary, solid Boolean calculation is utilized to generate three-dimensional geologic.

3. method for solving set up by a kind of tight gas pressure break horizontal well numerical well testing model according to claim 1, and it is characterized in that, the concrete steps that the plastid two-dimensionally in described step 2 and three-dimensional geologic carry out grid discrete are as follows:

First, Netgen open source software is bundled into merit and compiles, build running environment;

Then, the inner and outer boundary in plastid two-dimensionally and three-dimensional geologic is belonged to and forms two-dimensional grid file and three-dimensional grid file according to the requirement of Netgen grid file form respectively, then it is discrete to carry out grid according to the grid discrete step that Netgen is arranged.

4. method for solving set up by a kind of tight gas pressure break horizontal well numerical well testing model according to claim 1, it is characterized in that, described horizontal wellbore without the flow model in porous media calculating concrete grammar of pressure reduction is:

1) stratum and the fisstured flow equation of stress sensitive is considered

Stratum filtration equation:

$\begin{array}{c}{K}_{\mathrm{xD}}\frac{{\&PartialD;}^2{p}_{\mathrm{DR}}}{{\&PartialD;x}_D^2}+K_{\mathrm{yD}}\frac{{\&PartialD;}^2{p}_{\mathrm{DR}}}{{\&PartialD;y}_D^2}+K_{\mathrm{zD}}\frac{{\&PartialD;}^2{p}_{\mathrm{DR}}}{{\&PartialD;z}_D^2}-K_{\mathrm{xD}}{\mathrm{\&gamma;}}_D{\left(\frac{{\&PartialD;p}_{\mathrm{DR}}}{{\&PartialD;x}_D}\right)}^2\\ -K_{\mathrm{yD}}{\mathrm{\&gamma;}}_D{\left(\frac{{\&PartialD;p}_{\mathrm{DR}}}{{\&PartialD;y}_D}\right)}^2-K_{\mathrm{zD}}{\mathrm{\&gamma;}}_D{\left(\frac{{\&PartialD;p}_{\mathrm{DR}}}{{\&PartialD;z}_D}\right)}^2=\frac{e^{{\mathrm{\&gamma;}}_D{p}_{\mathrm{DR}}}}{K_{\mathrm{lD}}}\frac{{\&PartialD;p}_{\mathrm{DR}}}{{\&PartialD;t}_D}\end{array}---\left(1\right)$

Fisstured flow equation:

$\frac{{\&PartialD;}^2{p}_{\mathrm{Df}}}{{\&PartialD;r}_D^2}+K_{\mathrm{zD}}\frac{{\&PartialD;}^2{p}_{\mathrm{Df}}}{{\&PartialD;z}_D^2}=\frac{1}{K_{\mathrm{fD}}}\frac{{\&PartialD;p}_{\mathrm{Df}}}{{\&PartialD;t}_D}---\left(2\right)$

Primary condition:

p _D(x,y,z,0)＝0 (3)

Internal boundary condition:

$\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{A}_j{K}_{\mathrm{jD}}{e}^{-{\mathrm{\&gamma;}}_D{p}_{\mathrm{jD}}}{\left(\frac{{\&PartialD;p}_{\mathrm{jD}}}{{\&PartialD;n}^{\&prime;}}\right)|}_{{\mathrm{\&Gamma;}}_{\mathrm{in}}}={2\mathrm{\&pi;h}}_D(1-C_{\mathrm{DL}}\frac{{\mathrm{dp}}_{\mathrm{wD}}}{{\mathrm{dt}}_D})---\left(4\right)$

$p_{\mathrm{wD}}\left(t_D\right)=p_{\mathrm{jD}}-S_t\frac{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{A}_{\mathrm{jD}}{K}_{\mathrm{jD}}{e}^{-{\mathrm{\&gamma;}}_D{p}_{j}D}}{{2\mathrm{\&pi;h}}_D}\left(\frac{{\&PartialD;o}_{\mathrm{jD}}}{{\&PartialD;n}^{\&prime;}}\right)-{\mathrm{MP}}_{\mathrm{jD}}---\left(5\right)$

Outer Boundary Conditions:

Close: ${\frac{{\&PartialD;p}_D}{{\&PartialD;n}^{\&prime;}}|}_{{\mathrm{\&Gamma;}}_{\mathrm{out}}}=0---\left(6\right)$

Level pressure: ${p_D|}_{{\mathrm{\&Gamma;}}_{\mathrm{out}}}=0---\left(7\right)$

Symbol implication in formula:

P _dRfor the dimensionless pressure of subterranean formation zone; p _dffor the dimensionless pressure of crack area; t _dfor nondimensional time; C _dLfor dimensionless wellbore storage constant; K _xDfor x direction dimensionless permeability; K _yDfor y direction dimensionless permeability; K _zDfor z direction dimensionless permeability; K _fDfor dimensionless fracture permeabgility; γ _dfor dimensionless permeability modules; p _wDfor the dimensionless pressure of Article 1 crack and pit shaft point of intersection; MP _jDfor j point dimensionless pressure and p _wDbetween difference; A _jfor inner boundary triangle dimensionless area; h _dfor dimensionless reservoir thickness; S _tfor wellbore skin coefficient;

2) equation solution

First introduce conversion, by nonlinear filtration equation linearisation, then adopt Finite Element Method to solve, transformation for mula is:

$p_D=-\frac{1}{{\mathrm{\&gamma;}}_D}\ln (1-{\mathrm{\&gamma;}}_n\mathrm{\&eta;})---\left(8\right)$

Adopt mixed finite element method to be solved in the stratum after conversion and fisstured flow equations simultaneousness, the finite element equation of stratum and Fracture System is decomposed into the finite element equation of subterranean formation zone (on the right of formula 9 Section 1) and represents the finite element equation (on the right of formula 9 Section 2) of Fracture System.

${\&Integral;\&Integral;\&Integral;}_{\mathrm{\&Omega;}}\mathrm{FEQd\&Omega;}={\&Integral;\&Integral;\&Integral;}_{{\mathrm{\&Omega;}}_m}{\mathrm{FEQd\&Omega;}}_m+w_f\&CenterDot;{\&Integral;\&Integral;\&Integral;}_{\stackrel{\&OverBar;}{{\mathrm{\&Omega;}}_f}}\mathrm{FEQd}\stackrel{\&OverBar;}{{\mathrm{\&Omega;}}_f}---\left(9\right)$

A. subterranean formation zone three-dimensional finite element equation is:

$\begin{array}{c}{\mathrm{VK}}_{\mathrm{LD}}(K_{\mathrm{xD}}{b}_i^2+K_{\mathrm{yD}}{c}_i^2+K_{\mathrm{zD}}{d}_i^2+\frac{1}{{10K}_{\mathrm{LD}}{\mathrm{\&Delta;t}}_D}){\mathrm{\&eta;}}_i^{e,n+1}+{\mathrm{VK}}_{\mathrm{LD}}\left(\begin{array}{c}{K}_{\mathrm{xD}}{b}_i{b}_j+K_{\mathrm{yD}}{c}_i{c}_j+K_{\mathrm{zD}}{d}_i{d}_j\\ +\frac{1}{20K_{\mathrm{LD}}+{\mathrm{\&Delta;t}}_D}\end{array}\right){\mathrm{\&eta;}}_j^{e,n+1}\\ +{\mathrm{VK}}_{\mathrm{LD}}(K_{\mathrm{xD}}{b}_i{b}_k+K_{\mathrm{yD}}{c}_i{c}_k+K_{\mathrm{zD}}{d}_i{d}_k+\frac{1}{{20K}_{\mathrm{LD}}{\mathrm{\&Delta;t}}_D}){\mathrm{\&eta;}}_k^{e,n+1}+{\mathrm{VK}}_L\left(\begin{array}{c}{K}_{\mathrm{xD}}{b}_i{b}_m+K_{\mathrm{yD}}{c}_i{c}_m+\\ K_{\mathrm{zD}}{d}_i{d}_m+\frac{1}{{20K}_{\mathrm{LD}}{\mathrm{\&Delta;t}}_D}\end{array}\right){\mathrm{\&eta;}}_m^{e,n+1}\\ -\frac{{\mathrm{AK}}_{\mathrm{LD}}}{6}\frac{{\&PartialD;\mathrm{\&eta;}}_i^{e,n+1}}{{\&PartialD;n}^{\&prime;}}-\frac{{\mathrm{AK}}_{\mathrm{LD}}}{12}\frac{{\&PartialD;\mathrm{\&eta;}}_{j/k/m}^{e,n+1}}{{\&PartialD;n}^{\&prime;}}-\frac{{\mathrm{AK}}_{\mathrm{LD}}}{12}\frac{{\&PartialD;\mathrm{\&eta;}}_{k/m/j}^{e,n+1}}{{\&PartialD;n}^{\&prime;}}\\ =\frac{V}{10{\mathrm{\&Delta;t}}_D}{\mathrm{\&eta;}}_i^{e,n}+\frac{V}{{20\mathrm{\&Delta;t}}_D}{\mathrm{\&eta;}}_j^{e,n}+\frac{V}{{20\mathrm{\&Delta;t}}_D}{\mathrm{\&eta;}}_k^{e,n}+\frac{V}{20{\mathrm{\&Delta;t}}_D}{\mathrm{\&eta;}}_m^{e,n}\end{array}---\left(10\right)$

$\frac{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{A}_{\mathrm{jD}}{K}_{\mathrm{jD}}{\left(\frac{{\&PartialD;\mathrm{\&eta;}}_{\mathrm{jD}}}{{\&PartialD;n}^{\&prime;}}\right)|}_{{\mathrm{\&Gamma;}}_m}}{{2\mathrm{\&pi;h}}_D}+\frac{C_{\mathrm{DL}}}{{\mathrm{\&Delta;t}}_D}{\mathrm{\&eta;}}_w^{e,n+1}=1+\frac{C_{\mathrm{DL}}}{{\mathrm{\&Delta;t}}_D}{\mathrm{\&eta;}}_w^{e,n}---\left(11\right)$

${\mathrm{\&eta;}}_w\left(t_D\right)={\mathrm{\&eta;}|}_{{\mathrm{\&Gamma;}}_{\mathrm{in}}}-S_t\frac{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{SA}}_{\mathrm{jD}}{K}_{\mathrm{jD}}}{{2\mathrm{\&pi;h}}_D}\left(\frac{{\&PartialD;\mathrm{\&eta;}}_j}{{\&PartialD;n}^{\&prime;}}\right)---\left(12\right)$

B. fracture surface two dimensional finite element equation is:

$\begin{array}{c}A(K_{\mathrm{rD}}{b}_i^2+K_{\mathrm{zD}}{c}_i^2+\frac{1}{6{\mathrm{\&Delta;t}}_D}){\mathrm{\&eta;}}_i^{e,n+1}+A(K_{\mathrm{rD}}{b}_i{b}_j+K_{\mathrm{zD}}{c}_i{c}_j+\frac{1}{12{\mathrm{\&Delta;t}}_D}){\mathrm{\&eta;}}_j^{e,n+1}\\ +A(K_{\mathrm{rD}}{b}_i{b}_k+K_{\mathrm{zD}}{c}_i{c}_k+\frac{1}{12{\mathrm{\&Delta;t}}_D}){\mathrm{\&eta;}}_k^{e,n+1}-\frac{L}{3}\frac{{\&PartialD;\mathrm{\&eta;}}_i^{e,n+1}}{{\&PartialD;n}^{\&prime;}}-\frac{L}{6}\frac{{\&PartialD;\mathrm{\&eta;}}_{j,k}^{e,n+1}}{{\&PartialD;n}^{\&prime;}}\\ =\frac{A}{{6\mathrm{\&Delta;t}}_D}{\mathrm{\&eta;}}_i^{e,n}+\frac{A}{{12\mathrm{\&Delta;t}}_D}{\mathrm{\&eta;}}_j^{e,n}+\frac{A}{12{\mathrm{\&Delta;t}}_D}{\mathrm{\&eta;}}_k^{e,n}\end{array}---\left(13\right)$

$\frac{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{L}_{\mathrm{jD}}{K}_{\mathrm{jD}}\left(\frac{{\&PartialD;\mathrm{\&eta;}}_{\mathrm{jD}}}{{\&PartialD;n}^{\&prime;}}\right)}{2\mathrm{\&pi;}}+\frac{C_{\mathrm{DL}}}{{\mathrm{\&Delta;t}}_D}{\mathrm{\&eta;}}_w^{e,n+1}=1+\frac{C_{\mathrm{DL}}}{{\mathrm{\&Delta;t}}_D}{\mathrm{\&eta;}}_w^{e,n}---\left(14\right)$

${\mathrm{\&eta;}}_w\left(t_D\right)={\mathrm{\&eta;}}_j-S_t\frac{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{L}_{\mathrm{jD}}{K}_{\mathrm{jD}}}{{2\mathrm{\&pi;h}}_D}\left(\frac{{\&PartialD;\mathrm{\&eta;}}_j}{{\&PartialD;n}^{\&prime;}}\right)---\left(15\right)$

By finite element equation (10) ~ (15) simultaneous composition system stiffness matrix, the SuperLU numerical solution device of parallelization is utilized to solve large linear systems, pressure field distribution and the inner boundary normal pressure gradient of whole reservoir can be obtained, then calculate each crack production flow thus:

$Q_{\mathrm{fi}}=Q_{\mathrm{sc}}\frac{w_{\mathrm{fD}}{L}_{\mathrm{iD}}{K}_{\mathrm{fD}}{e}^{{-\mathrm{\&gamma;}}_D{p}_{\mathrm{iD}}}}{2{\mathrm{\&pi;h}}_D}{\left(\frac{{\&PartialD;p}_{\mathrm{iD}}}{{\&PartialD;n}^{\&prime;}}\right)|}_{{\mathrm{\&Gamma;}}_{\mathrm{in}}}---\left(16\right)$

Symbol implication in formula:

η is transformation to linearity parameter; η _wfor the conversion parameter corresponding to dimensionless bottom pressure value; w _ffor crack width, m; w _fDfor dimensionless crack width; L _jDfor crack inner boundary unit wires length; V is tetrahedron volume; B, c, d are finite element coefficient; I, j, k, m are finite element tetrahedron four summit sequence numbers; Q _fibe the flow of the i-th crack, m ³/ d; Q _scfor the flow of gas well under mark condition, m ³/ d.

5. method for solving set up by a kind of tight gas pressure break horizontal well numerical well testing model according to claim 1, it is characterized in that, sets up coupling model in described step 4, and solves the coupling model set up and carry out according to following steps:

A, according to the wellbore pressure p that step 4 calculates _wDwith crack flow Q _fi, adopt pit shaft multiphase flow design formulas to calculate, obtain the dimensionless pressure reduction MP of pit shaft each point _iD, concrete formula is as follows:

The basic equation that pit shaft multiphase flow calculates is:

$-\frac{{\mathrm{dp}}_w}{\mathrm{dZ}}=\frac{({\mathrm{\&rho;}}_l{H}_L+{\mathrm{\&rho;}}_g(1-H_L))\sin \mathrm{\&theta;}+\frac{{\mathrm{\&lambda;Gv}}_m}{2\mathrm{AD}}}{1-\left\{\right[{\mathrm{\&rho;}}_l{H}_L+{\mathrm{\&rho;}}_g(1-H_L)\left]v_m{v}_{\mathrm{sg}}\right\}/p}---\left(17\right)$

Symbol implication in formula:

ρ _lfor fluid density, kg/m ³; ρ _gfor gas density, kg/m ³; G is gas-liquid mixture mass flow, kg/s; v _mfor mixture flowing velocity, m/s; v _sgfor gas superficial flow velocity, m/s; A is pit shaft oil pipe sectional area, m ²; D is pipe aperture, m.

Liquid holdup H _ladopt Beggs-Brill method to calculate with coefficient of frictional resistance λ, obtain wellbore pressure gradient according to formula (17), then to add up the pressure difference obtained between each point and shaft bottom standard point according to pit shaft inner boundary relative distance:

${\mathrm{MP}}_i=\underset{j=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-\frac{{\mathrm{dp}}_w}{\mathrm{dZ}})}_i{\mathrm{\&Delta;Z}}_j---\left(18\right)$

By each point pressure reduction MP obtained _icarry out nondimensionalization, each point dimensionless pressure reduction MP can be obtained _iD.

B, the inner boundary each point pressure reduction MP obtained in steps A _iD, be brought in flow model in porous media and calculate wellbore pressure p ' _wDwith crack flow Q ' _fi, the coupling condition of pit shaft multiphase flow model and flow model in porous media is as follows:

On stratum and pit shaft interface:

${p_{\mathrm{DRi}}|}_{{\mathrm{\&Gamma;}}_{\mathrm{in}}}=p_{\mathrm{wD}}+{\mathrm{MP}}_{\mathrm{iD}}---\left(19\right)$

On crack and pit shaft interface:

${p_{\mathrm{Dfj}}|}_{{\mathrm{\&Gamma;}}_{\mathrm{in}}}=p_{\mathrm{wD}}+{\mathrm{MP}}_{\mathrm{jD}}---\left(20\right)$

C, the p ' that front and back iteration step calculates _wDand p _wDsubtract each other, when poor absolute value is less than ε, then continue the calculating of future time step, when the absolute value of both differences is more than or equal to ε, then the current p ' obtained between the two _wDand Q ' _finew wellbore pressure p is obtained after being averaged with the value of back _wDwith crack flow Q _fi, and be brought in steps A and carry out iteration, until the p ' obtained _wDwith p _wDthe absolute value carrying out the difference of subtracting each other is less than ε, and the general value of ε is 10 ^-4;

D, the p ' that step C determines _wDwith the Q ' of correspondence _ficarry out record, if now time step k< total time walks n, then continue the calculating of future time step, according to current wellbore pressure p _wD, crack flow Q _fi, from steps A, calculate the p ' of a new round _wDand Q ' _fi;

E, obtains p ' institute in n calculating in step D _wDand Q ' _fivalue generate Well Testing Theory curve.

6. method for solving set up by a kind of tight gas pressure break horizontal well numerical well testing model according to claim 1, it is characterized in that, described step 5 mainly theoretical curve and measured curve contrasts matching in log-log graph, semilogarithmic plot and full historical pressures curve map, according to the fitting degree of curve, the adjustment of fitting parameter can be carried out, finally make theoretical curve and measured curve can both be coincide preferably in log-log graph, semilogarithmic plot and full historical pressures curve map, curve just completes; After curve completes, namely can obtain the parameter of well test analysis, comprise fracture parameters and reservoir parameter.