CN105184061A - Numerical simulation method of temperature and pressure distribution of gas producing well - Google Patents

Abstract

The invention relates to the technical field of oil and gas reservoir development engineering management, and discloses a numerical simulation method of temperature and pressure distribution of a gas producing well, which is used for accurately simulating the distribution conditions of temperature and pressure parameters of the gas producing well to provide guidance for the production. The method comprises the following steps: A, establishing a temperature and pressure coupling differential equation model; B, estimating related parameters of differential equations; C, using a difference method for algorithm design; D, constructing simulation software based on a C# development platform; and E, inputting initial data in the simulation software to obtain a simulation result of the temperature and pressure distribution. The numerical simulation method disclosed by the invention can be used for accurately predicting the temperature and pressure distribution of the gas producing well to greatly improve the design level of oil and gas mining equipment, which is beneficial to the development oil and gas reservoirs.

Description

The method for numerical simulation of gas-producing well temperature, pressure distribution

Technical field

The present invention relates to development of oil and gas reservoir administrative skill field, be specially the method for numerical simulation of gas-producing well temperature, pressure distribution, relate to the temperature, pressure coupling mechanism analysis of gas-producing well, Mathematical Models, method for numerical simulation design etc.

Background technology

Produce gas well usually along with high temperature, high flow rate, high pressure, this is due to reasons such as friction, pit shaft deformation, heat transfer.Along with the evolution of production environment, comprise deep water and high-temperature pressure condition, optimize in the design of petroleum engineering facility, hydrate prevention, the performance analysis of producing well, be necessary to carry out accurate prediction to temperature, pressure.

When liquid produces at first from a region, its temperature may think the same in shaft bottom.This is invalid to gas.If Joule-Thompson effect is appropriately considered, gasinlet temperature available stratigraphic temperature is estimated.Therefore, the temperature of wellbore bottom, can be reliably estimated.But along with the rising of fluid, its temperature is significantly higher than the temperature of surrounding formation, this is because formation temperature reduces along with the degree of depth and declines.

When stratum and fluid have temperature difference, heat transfer phenomenon will be there is.In any degree of depth, formation temperature is not only relevant with radiation length, also relevant with the production time.When reaching steady flow, turbulent flow ensure that in certain degree of depth, and the temperature of fluid is constant.So the thermal loss in fluid reduced along with the time, and depend on the different thermal resistance of pit shaft red heat liquid and surrounding formation.

Complete system by fluid, the annular gap containing low-pressure air, sleeve pipe, setting and stratum, as shown in Figure 1:

Catheter diameter is r _ti, external diameter is r _t0, the internal diameter of sleeve pipe is r _ci, external diameter is r _c0, heat is conducted by the air in ring.Radiation and convection current also occur.When body is heated, the conduction velocity of radiation energy depends on the temperature of body.The conduction of the radiation energy between conduit and sleeve pipe, depends on that interface sends and absorbs the feature of heat.In a lot of situation, the gap between conduit and hole is by shutoff.Because the conduction of cement may, lower than the stratum of periphery, calculate by pit shaft stage, progressively upwards.Using the intersection point of the air intake opening of gas well as coordinate axis, pit shaft is downwards the positive dirction of coordinate axis.The pressure deduction that Fig. 2 shows pit shaft infinitesimal is analysed, and P is hydrodynamic pressure, and v is gas velocity, and l is the degree of depth, and dv is the speed increment on dl, and dp is the pressure increment on dl, and θ is the pitch angle of pit shaft.

Summary of the invention

Technical matters to be solved by this invention is: the method for numerical simulation proposing the distribution of a kind of gas-producing well temperature, pressure, carries out accurate simulation to the temperature, pressure parameter distribution situation of gas-producing well, thus provides guidance foundation for producing.

The present invention solves the problems of the technologies described above adopted scheme: the method for numerical simulation of gas-producing well temperature, pressure distribution, comprises the following steps:

A, set up temperature, pressure coupled-differential equations model;

B, estimation differential equation group correlation parameter;

C, use method of difference carry out algorithm design;

D, based on C# development platform build simulation softward;

E, in simulation softward, input primary data, obtain temperature, pressure distribution simulation result.

Further, in steps A, described temperature, pressure coupled-differential equations model of setting up specifically comprises:

Model construction condition based on hypothesis obtains energy conservation equation, mass-conservation equation, momentum conservation equation respectively:

Wherein, energy conservation equation:

$q+\frac{d}{dl}(h+v^2/2-gl\sin \mathrm{\θ})=0---\left(1\right)$

Mass-conservation equation:

$\mathrm{\ρ}\frac{dv}{dl}+v\frac{d\mathrm{\ρ}}{dl}=0---\left(2\right)$

Momentum conservation equation:

$-\frac{1}{\mathrm{\ρ}}\frac{dP}{dl}+f\frac{v^2}{2r_{ti}}+gsin\mathrm{\θ}=\frac{vdv}{dl}---\left(3\right)$

On infinitesimal dl, the radiation heat from fluid to setting-bed boundary conduction is

$q=\frac{(-2\mathrm{\π})r_{ti}{U}_{ti}}{w}(T-T_h)---\left(4\right)$

The radiation heat conducted from setting-bed boundary to surrounding formation is:

$q=-\frac{2\mathrm{\π}ke}{wf\left(t_D\right)}(T_h-T_e)---\left(5\right)$

Convolution (4) and formula (5) obtain fluid and the heat conducting differential equation of surrounding formation:

$q=\frac{2\mathrm{\π}}{w}\[\frac{r_{ti}{U}_{ti}ke}{r_{ti}{U}_{ti}f\left(t_D\right)+ke}\](T-T_e)---\left(6\right)$

(1) and (6) is integrated, obtains following ordinary differential equation:

$\frac{2\mathrm{\π}}{w}\[\frac{r_{ti}{U}_{ti}ke}{r_{ti}{U}_{ti}\left(t_D\right)+ke}\](T-T_e)+\frac{d}{dl}(h+v^2/2-gl\sin \mathrm{\θ})=0---\left(7\right)$

When fluid flows in the wellbore time, because caliber change is very little, Joule – Thomson coefficient is negligible, therefore dh=C _pdT, the enthalpy change amount overall for formula (7) is:

$\frac{2\mathrm{\π}}{w}\[\frac{r_{ti}{U}_{ti}ke}{r_{ti}{U}_{ti}f\left(t_D\right)+ke}\](T-T_e)+C_p\frac{dT}{dl}+v\frac{dv}{dl}-g\sin \mathrm{\θ}=0---\left(8\right)$

The differential prescription journey of hydrodynamic pressure, temperature, density and speed in gasser can be obtained by formula (1), (2), (3), (8):

$\{\begin{array}{c}\frac{d\mathrm{\ρ}}{dl}=\frac{\frac{-8.314z\mathrm{\ρ}}{C_{p}M}\[\frac{\left(2\mathrm{\π}\right)r_{ti}{U}_{ti}ke}{w\[r_{ti}{U}_{ti}f\left(t_D\right)+ke\]}(T-T_e)-g\cos \mathrm{\θ}\]-f\frac{v^2\mathrm{\ρ}}{2r_{ti}}-g\mathrm{\ρ}\cos \mathrm{\θ}}{v^2-(\frac{8.314{\mathrm{zv}}^2}{C_{p}M}+\frac{8.314Tz}{M})}\\ \frac{dv}{dl}=-\frac{v}{\mathrm{\ρ}}\frac{d\mathrm{\ρ}}{dz}\\ \frac{dP}{dl}=\mathrm{\ρ}g\cos \mathrm{\θ}+f\frac{v^2\mathrm{\ρ}}{2r_{ti}}+v^2\frac{d\mathrm{\ρ}}{dl}\\ \frac{dT}{dl}=\frac{\[\frac{v^2}{\mathrm{\ρ}}\frac{d\mathrm{\ρ}}{dl}+g\cos \mathrm{\θ}-\frac{\left(2\mathrm{\π}\right)r_{ti}{U}_{ti}ke}{w\[r_{ti}{U}_{ti}f\left(t_D\right)+ke\]}(T-Te)\]}{C_p}\end{array}---\left(9\right)$

Use y _i(i=1,2,3,4) replace ρ, v, P, T, then system of equations can be reduced to

$\frac{{\mathrm{dy}}_i}{dl}=F_i(y_1,y_2,y_3,y_4),(i=1,2,3,4)---\left(10\right).$

Further, the model construction condition of described hypothesis comprises:

The flowing of rock gas is one way stable; The heat transfer of pit shaft is stable state; Stratum heat transfer is unstable; Conduit and sleeve pipe are concentric.

Further, in step B, estimate that differential equation group correlation parameter specifically comprises:

B1. heat-conduction coefficient is estimated:

$U_{ti}^{-1}=\frac{1}{hc+hr}+\frac{r_{ti}\ln \left(\frac{r_{cem}}{r_{co}}\right)}{k_{cem}}+\frac{r_{ti}\ln \left(\frac{r_{ci}}{r_{to}}\right)}{k_{ang}}---\left(11\right)$

B2. specific heat at constant pressure C is calculated _p:

C _p＝1243+3.14T+7.931×10 ^-4T ²-6.881×10 ^-7T ³(12)

B3. stratum heat diffusion equation:

$\left\{\begin{array}{cc}f\left(t_D\right)=1.1281\sqrt{t_D}(1-0.3\sqrt{t_D}& (t_D\≤1.5)\\ f\left(t_D\right)=(1+\frac{0.6}{t_D})\[0.4063+0.5\ln \left(t_D\right)\]& (t_D>1.5)\end{array}\right.---\left(13\right)$

$t_D=\frac{t\mathrm{\α}}{r_{wb}^2}$

B4. friction factor f is calculated:

$\frac{1}{\sqrt{f}}=1.14-2\lg (\frac{0.00001524}{r_{ti}}+\frac{21.25}{{\mathrm{Re}}^{0.9}})---\left(14\right).$

Further, in step C, described use method of difference is carried out algorithm design and is comprised solving temperature, pressure coupled-differential equations model:

Gas-producing well is divided into n interval, j (j=1,2 ..., n) be node, n is shaft bottom, and step-length is h;

The boundary condition of temperature, pressure coupled-differential equations is:

$P\[n\]=P_0,T\[n\]=T_0,\mathrm{\ρ}\[n\]=0.000001\×3484.48{\mathrm{\γ}}_g\frac{P_0}{{\mathrm{zT}}_0},v\[n\]=\frac{101000\×300000T_0}{293\×86400P_{0}A}---\left(15\right)$

In order to calculate the correlative of each node, use following iterative equation:

$y_i\[j-1\]=y_i\[j\]-\frac{h}{6}(a\[i\]+2b\[i\]+2c\[i\]+d\[i\]),(i=1,2,3,4;j=1,2,3,4...n)---\left(16\right)$

For using above equation, need to estimate a [i], b [i], c [i] and d [i]:

$\left\{\begin{array}{c}a\[i\]=F_i(y_1\[n\],y_2\[n\],y_3\[n\],y_4\[n\])\\ b\[i\]=F_i(y_1\[n\]+\frac{h}{2}a\[1\],y_2\[n\]+\frac{h}{2}a\[2\],y_3\[n\]+\frac{h}{2}a\[3\],y_4\[n\]+\frac{h}{2}a\[4\])\\ c\[i\]=F_i(y_1\[n\]+\frac{h}{2}b\[1\],y_2\[n\]+\frac{h}{2}b\[2\],y_3\[n\]+\frac{h}{2}b\[3\],y_4\[n\]+\frac{h}{2}b\[4\])\\ d\[i\]=F_i(y_1\[n\]+hc\[1\],y_2\[n\]+hc\[2\],y_3\[n\]+hc\[3\],y_4\[n\]+hc\[4\])\end{array}\right.---\left(17\right).$

Further, in step D, described simulation softward comprises subscriber interface module, exports output module, algoritic module and chart module;

Described subscriber interface module is for providing the initial information of gas well, and these information can be divided into static data, fluid data and production data:

Static data: describe static lower gas well feature;

Fluid data: the state describing fluid;

Production data: the parameter that the production phase is given;

Described input/output module is used for reading in from storer, preserving data;

Described algoritic module is used for carrying out mathematical computations according to the algorithm of design;

Described chart module is used for the data result of simulation to carry out pictorialization displaying.

The invention has the beneficial effects as follows: accurately predicting is carried out to the temperature and pressure distribution of gas-producing well, greatly promotes the raising of oil gas production equipment design level, thus be conducive to Reservoir Development.

Accompanying drawing explanation

Fig. 1 is well segment structure schematic diagram;

Fig. 2 is that in conduit, pressure deduction analyses schematic diagram;

Fig. 3 is the method for numerical simulation process flow diagram of gas-producing well temperature, pressure of the present invention distribution;

Fig. 4 is simulation software construction schematic diagram;

Fig. 5 is program circuit schematic diagram;

Fig. 6 is simulation pressure distribution curve figure out;

Fig. 7 is simulation temperature distributing curve diagram out.

Embodiment

The present invention is intended to the method for numerical simulation proposing the distribution of a kind of gas-producing well temperature, pressure, carries out accurate simulation to the temperature, pressure parameter distribution situation of gas-producing well, thus provides guidance foundation for producing.

The parameter lexical or textual analysis likely related in the present invention is as follows:

A: inner tube surface amasss, m ²;

C _p: specific heat at constant pressure, J/Kg.K;

F: friction factor;

F (T _d): stratum transient heat conduction equation of time;

G: acceleration of gravity;

H: specific enthalpy;

K _ang: air heat conduction in ring, J/m.K;

K _cem: setting heat conductivity, J/m.K;

K _e: stratum heat conductivity, J/m.K;

L: pit shaft length, l;

M: the molar average weight of gas, g/mol;

P: fluid pressure, pa;

P _pe: critical pressure, pa;

Rcem: setting external diameter, m

Re: Reynolds number;

T: production time, s;

T: fluid temperature (F.T.), K;

T _e: surrounding formation temperature, K;

U _ti: overall heat transfer coefficient, W/m.K;

V: natural gas rate, m/s;

W: total mass velocity;

Z: compressibility factor;

α: the thermal diffusivity on stratum, m ²/ s.

Do further to describe to the solution of the present invention below in conjunction with drawings and Examples:

As shown in Figure 3, the method for numerical simulation of gas-producing well temperature, pressure distribution of the present invention comprises:

A, set up temperature, pressure coupled-differential equations model;

B, estimation differential equation group correlation parameter;

C, use method of difference carry out algorithm design;

D, based on C# development platform build simulation softward;

E, in simulation softward, input primary data, obtain temperature, pressure distribution simulation result.

On concrete enforcement, in steps A, the model construction condition based on hypothesis obtains energy conservation equation, mass-conservation equation, momentum conservation equation respectively;

For the situation of realistic gas-producing well, make the following assumptions:

(1) flowing of rock gas is one way stable;

(2) heat transfer of pit shaft is considered to stable state;

(3) based on dimensionless transient heat conduct equation of time, stratum heat transfer is unstable;

(4) conduit and sleeve pipe are concentric;

Energy conservation equation:

$q+\frac{d}{dl}(h+v^2/2-gl\sin \mathrm{\θ})=0---\left(1\right)$

Mass-conservation equation:

$\mathrm{\ρ}\frac{dv}{dl}+v\frac{d\mathrm{\ρ}}{dl}=0---\left(2\right)$

Momentum conservation equation:

$-\frac{1}{\mathrm{\ρ}}\frac{dP}{dl}+f\frac{v^2}{2r_{ti}}+gsin\mathrm{\θ}=\frac{vdv}{dl}---\left(3\right)$

Coupled-differential equations model: on infinitesimal dl, the radiation heat from fluid to setting-bed boundary conduction is

$q=\frac{(-2\mathrm{\π})r_{ti}{U}_{ti}}{w}(T-T_h)---\left(4\right)$

The radiation heat conducted from setting-bed boundary to surrounding formation is

$q=\frac{2\mathrm{\π}}{w}\[\frac{r_{ti}{U}_{ti}ke}{r_{ti}{U}_{ti}f\left(t_D\right)+ke}\](T-T_e)---\left(5\right)$

By on (4) and (5), obtain fluid and the heat conducting differential equation of surrounding formation:

$q=\frac{2\mathrm{\π}}{w}\[\frac{r_{ti}{U}_{ti}ke}{r_{ti}{U}_{ti}f\left(t_D\right)+ke}\](T-T_e)---\left(6\right)$

(1) and (6) is integrated, obtains the following differential equation

$\frac{2\mathrm{\π}}{w}\[\frac{r_{ti}{U}_{ti}ke}{r_{ti}{U}_{ti}\left(t_D\right)+ke}\](T-T_e)+\frac{d}{dl}(h+v^2/2-gl\sin \mathrm{\θ})=0---\left(7\right)$

When fluid flows in the wellbore time, because caliber change is very little, Joule – Thomson coefficient is negligible.Therefore, dh=C _pdT, overall enthalpy change amount is

$\frac{2\mathrm{\π}}{w}\[\frac{r_{ti}{U}_{ti}ke}{r_{ti}{U}_{ti}f\left(t_D\right)+ke}\](T-T_e)+C_p\frac{dT}{dl}+v\frac{dv}{dl}-g\sin \mathrm{\θ}=0---\left(8\right)$

(1) (2) (3) (8) are put together, arrange, the differential prescription journey of the hydrodynamic pressure in gasser, temperature, density and speed can be obtained:

$\{\begin{array}{c}\frac{d\mathrm{\ρ}}{dl}=\frac{\frac{-8.314z\mathrm{\ρ}}{C_{p}M}\[\frac{\left(2\mathrm{\π}\right)r_{ti}{U}_{ti}ke}{w\[r_{ti}{U}_{ti}f\left(t_D\right)+ke\]}(T-T_e)-g\cos \mathrm{\θ}\]-f\frac{v^2\mathrm{\ρ}}{2r_{ti}}-g\mathrm{\ρ}\cos \mathrm{\θ}}{v^2-(\frac{8.314{\mathrm{zv}}^2}{C_{p}M}+\frac{8.314Tz}{M})}\\ \frac{dv}{dl}=-\frac{v}{\mathrm{\ρ}}\frac{d\mathrm{\ρ}}{dz}\\ \frac{dP}{dl}=\mathrm{\ρ}g\cos \mathrm{\θ}+f\frac{v^2\mathrm{\ρ}}{2r_{ti}}+v^2\frac{d\mathrm{\ρ}}{dl}\\ \frac{dT}{dl}=\frac{\[\frac{v^2}{\mathrm{\ρ}}\frac{d\mathrm{\ρ}}{dl}+g\cos \mathrm{\θ}-\frac{\left(2\mathrm{\π}\right)r_{ti}{U}_{ti}ke}{w\[r_{ti}{U}_{ti}f\left(t_D\right)+ke\]}(T-Te)\]}{C_p}\end{array}---\left(9\right)$

Use y _i(i=1,2,3,4) replace ρ, v, P, T, then system of equations can be reduced to

$\frac{{\mathrm{dy}}_i}{dl}=F_i(y_1,y_2,y_3,y_4),(i=1,2,3,4)---\left(10\right)$

Estimate described in step B that differential equation group correlation parameter comprises:

B1. heat-conduction coefficient is estimated:

$U_{ti}^{-1}=\frac{1}{hc+hr}+\frac{r_{ti}\ln \left(\frac{r_{cem}}{r_{co}}\right)}{k_{cem}}+\frac{r_{ti}\ln \left(\frac{r_{ci}}{r_{to}}\right)}{k_{ang}}---\left(11\right)$

Above formula shows in wellbore system, and between sleeve pipe/conduit annular space, material (gas, oil, water or potpourri) plays important effect in the heat conducting process of decision.

B2. specific heat at constant pressure C is calculated _p:

C _p＝1243+3.14T+7.931×10 ^-4T ²-6.881×10 ^-7T ³(12)

B3. stratum heat diffusion equation:

$\left\{\begin{array}{cc}f\left(t_D\right)=1.1281\sqrt{t_D}(1-0.3\sqrt{t_D}& (t_D\≤1.5)\\ f\left(t_D\right)=(1+\frac{0.6}{t_D})\[0.4063+0.5\ln \left(t_D\right)\]& (t_D>1.5)\end{array}\right.---\left(13\right)$

$t_D=\frac{t\mathrm{\α}}{r_{wb}^2}$

B4. friction factor f is calculated:

$\frac{1}{\sqrt{f}}=1.14-2\lg (\frac{0.00001524}{r_{ti}}+\frac{21.25}{{\mathrm{Re}}^{0.9}})---\left(14\right).$

Method of difference is used to carry out algorithm design in step C; Due to C _pbe not constant with θ, pit shaft must be divided into some length of interval, and step-length is h.Suppose that well can be divided into n interval, j (j=1,2 ..., n) be node, n is shaft bottom.According to temperature and the pressure in shaft bottom, density and the fluid velocity of corresponding rock gas can be calculated.The boundary condition of the above-mentioned differential equation is:

$P\[n\]=P_0,T\[n\]=T_0,\mathrm{\ρ}\[n\]=0.000001\×3484.48{\mathrm{\γ}}_g\frac{P_0}{{\mathrm{zT}}_0},v\[n\]=\frac{101000\×300000T_0}{293\×86400P_{0}A}$

In order to calculate the correlative of each node, use following iterative equation;

$y_i\[j-1\]=y_i\[j\]-\frac{h}{6}(a\[i\]+2b\[i\]+2c\[i\]+d\[i\]),(i=1,2,3,4;j=1,2,3,4...n)$

For using above equation, need to estimate a [i], b [i], c [i] and d [i]:

$\left\{\begin{array}{c}a\[i\]=F_i(y_1\[n\],y_2\[n\],y_3\[n\],y_4\[n\])\\ b\[i\]=F_i(y_1\[n\]+\frac{h}{2}a\[1\],y_2\[n\]+\frac{h}{2}a\[2\],y_3\[n\]+\frac{h}{2}a\[3\],y_4\[n\]+\frac{h}{2}a\[4\])\\ c\[i\]=F_i(y_1\[n\]+\frac{h}{2}b\[1\],y_2\[n\]+\frac{h}{2}b\[2\],y_3\[n\]+\frac{h}{2}b\[3\],y_4\[n\]+\frac{h}{2}b\[4\])\\ d\[i\]=F_i(y_1\[n\]+hc\[1\],y_2\[n\]+hc\[2\],y_3\[n\]+hc\[3\],y_4\[n\]+hc\[4\])\end{array}\right.$

In step D, build simulation softward based on C# development platform; As shown in Figure 4, it comprises subscriber interface module, input/output module, algoritic module and icon module to software configuration;

(1) subscriber interface module: the initial information that gas well is mainly provided, these information can be divided into static data, fluid data and production data:

(i) static data: describe static lower gas well feature, be abbreviated as S.

(ii) fluid data: the state describing fluid, Ru Shui, oil or gas, be abbreviated as F.

(iii) production data: the parameter that the production phase is given, is abbreviated as P.

(2) input/output module: this module is read in, preserved data from storer.MicrosoftSQLServer2000 is used to set up the database of pit shaft.

(3) algoritic module: this module will be applied to all mathematical computations.The temperature and pressure from shaft bottom to well head arbitrary node can be calculated.

(4) chart module: this module is by data drawing list.By CrystalReports can draw temperature, pressure distribution curve and Output rusults.

After simulation softward has designed, just can utilize and input primary data in simulation softward, obtain temperature, pressure distribution simulation result, specific algorithm flow process as shown in Figure 5.

Embodiment:

For certain a bite gas-producing well, its input data are in table 1:

Table 1: gas-producing well correlated inputs tables of data

By software simulation, obtain series of results.As Fig. 6 and Fig. 7, temperature and pressure is from shaft bottom to well head, more and more less.This is because the heat transfer of gravity, friction and pit shaft formation.Temperature is at 0 to 3000 meters, and decline a lot, but tend towards stability subsequently, temperature is then with degree of depth linear change.Fig. 6 shows, and at same well depth, gas production rate is larger, and because frictional resistance volume increase pressure is less, because flow velocity increases, temperature raises.

Claims (6)

1. the method for numerical simulation of gas-producing well temperature, pressure distribution, is characterized in that, comprise the following steps:

A, set up temperature, pressure coupled-differential equations model;

B, estimation differential equation group correlation parameter;

C, use method of difference carry out algorithm design;

D, based on C# development platform build simulation softward;

E, in simulation softward, input primary data, obtain temperature, pressure distribution simulation result.

2. the method for numerical simulation of gas-producing well temperature, pressure distribution as claimed in claim 1, it is characterized in that, in steps A, described temperature, pressure coupled-differential equations model of setting up specifically comprises:

Model construction condition based on hypothesis obtains energy conservation equation, mass-conservation equation, momentum conservation equation respectively:

Wherein, energy conservation equation:

$q+\frac{d}{dl}(h+v^2/2-gl\sin \mathrm{\θ})=0---\left(1\right)$

Mass-conservation equation:

$\mathrm{\ρ}\frac{dv}{dl}+v\frac{d\mathrm{\ρ}}{dl}=0---\left(2\right)$

Momentum conservation equation:

$-\frac{1}{\mathrm{\ρ}}\frac{dP}{dl}+f\frac{v^2}{2r_{ti}}+gsin\mathrm{\θ}=\frac{vdv}{dl}---\left(3\right)$

On infinitesimal dl, the radiation heat from fluid to setting-bed boundary conduction is

$q=\frac{(-2\mathrm{\π})r_{ti}{U}_{ti}}{w}(T-T_h)---\left(4\right)$

The radiation heat conducted from setting-bed boundary to surrounding formation is:

$q=-\frac{2\mathrm{\π}ke}{wf\left(t_D\right)}(T_h-T_e)---\left(5\right)$

Convolution (4) and formula (5) obtain fluid and the heat conducting differential equation of surrounding formation:

$q=\frac{2\mathrm{\π}}{w}\[\frac{r_{ti}{U}_{ti}ke}{r_{ti}{U}_{ti}f\left(t_D\right)+ke}\](T-T_e)---\left(6\right)$

(1) and (6) is integrated, obtains following ordinary differential equation:

$\frac{2\mathrm{\π}}{w}\[\frac{r_{ti}{U}_{ti}ke}{r_{ti}{U}_{ti}\left(t_D\right)+ke}\](T-T_e)+\frac{d}{dl}(h+v^2/2-gl\sin \mathrm{\θ})=0---\left(7\right)$

When fluid flows in the wellbore time, because caliber change is very little, Joule – Thomson coefficient is negligible, therefore dh=C _pdT, the enthalpy change amount overall for formula (7) is:

$\frac{2\mathrm{\π}}{w}\[\frac{r_{ti}{U}_{ti}ke}{r_{ti}{U}_{ti}f\left(t_D\right)+ke}\](T-T_e)+C_p\frac{dT}{dl}+v\frac{dv}{dl}-g\sin \mathrm{\θ}=0---\left(8\right)$

The differential prescription journey of hydrodynamic pressure, temperature, density and speed in gasser can be obtained by formula (1), (2), (3), (8):

$\{\begin{array}{c}\frac{d\mathrm{\ρ}}{dl}=\frac{\frac{-8.314z\mathrm{\ρ}}{C_{p}M}\[\frac{\left(2\mathrm{\π}\right)r_{ti}{U}_{ti}ke}{w\[r_{ti}{U}_{ti}f\left(t_D\right)+ke\]}(T-T_e)-g\cos \mathrm{\θ}\]-f\frac{v^2\mathrm{\ρ}}{2r_{ti}}-g\mathrm{\ρ}\cos \mathrm{\θ}}{v^2-(\frac{8.314{\mathrm{zv}}^2}{C_{p}M}+\frac{8.314Tz}{M})}\\ \frac{dv}{dl}=-\frac{v}{\mathrm{\ρ}}\frac{d\mathrm{\ρ}}{dz}\\ \frac{dP}{dl}=\mathrm{\ρ}g\cos \mathrm{\θ}+f\frac{v^2\mathrm{\ρ}}{2r_{ti}}+v^2\frac{d\mathrm{\ρ}}{dl}\\ \frac{dT}{dl}=\frac{\[\frac{v^2}{\mathrm{\ρ}}\frac{d\mathrm{\ρ}}{dl}+g\cos \mathrm{\θ}-\frac{\left(2\mathrm{\π}\right)r_{ti}{U}_{ti}ke}{w\[r_{ti}{U}_{ti}f\left(t_D\right)+ke\]}(T-Te)\]}{C_p}\end{array}---\left(9\right)$

Use y _i(i=1,2,3,4) replace ρ, v, P, T, then system of equations can be reduced to

$\frac{{\mathrm{dy}}_i}{dl}=F_i(y_1,y_2,y_3,y_4),(i=1,2,3,4)---\left(10\right).$

3. the method for numerical simulation of gas-producing well temperature, pressure distribution as claimed in claim 2, it is characterized in that, the model construction condition of described hypothesis comprises:

The flowing of rock gas is one way stable; The heat transfer of pit shaft is stable state; Stratum heat transfer is unstable; Conduit and sleeve pipe are concentric.

4. the method for numerical simulation of gas-producing well temperature, pressure distribution as claimed in claim 1, is characterized in that, in step B, estimates that differential equation group correlation parameter specifically comprises:

B1. heat-conduction coefficient is estimated:

$U_{ti}^{-1}=\frac{1}{hc+hr}+\frac{r_{ti}\ln \left(\frac{r_{cem}}{r_{co}}\right)}{k_{cem}}+\frac{r_{ti}\ln \left(\frac{r_{ci}}{r_{to}}\right)}{k_{ang}}---\left(11\right)$

B2. specific heat at constant pressure C is calculated _p:

C _p＝1243+3.14T+7.931×10 ^-4T ²-6.881×10 ^-7T ³(12)

B3. stratum heat diffusion equation:

$\left\{\begin{array}{cc}f\left(t_D\right)=1.1281\sqrt{t_D}(1-0.3\sqrt{t_D}& (t_D\≤1.5)\\ f\left(t_D\right)=(1+\frac{0.6}{t_D})\[0.4063+0.5\ln \left(t_D\right)\]& (t_D>1.5)\end{array}\right.---\left(13\right)$

$t_D=\frac{t\mathrm{\α}}{r_{wb}^2}$

B4. friction factor f is calculated:

$\frac{1}{\sqrt{f}}=1.14-2\lg (\frac{0.00001524}{r_{ti}}+\frac{21.25}{{\mathrm{Re}}^{0.9}})---\left(14\right).$

5. the method for numerical simulation of gas-producing well temperature, pressure distribution as claimed in claim 1, it is characterized in that, in step C, described use method of difference is carried out algorithm design and is comprised solving temperature, pressure coupled-differential equations model:

Gas-producing well is divided into n interval, j (j=1,2 ..., n) be node, n is shaft bottom, and step-length is h;

The boundary condition of temperature, pressure coupled-differential equations is:

$P\[n\]=P_0,T\[n\]=T_0,\mathrm{\ρ}\[n\]=0.000001\×3484.48{\mathrm{\γ}}_g\frac{P_0}{{\mathrm{zT}}_0},v\[n\]=\frac{101000\×300000T_0}{293\×86400P_{0}A}---\left(15\right)$

In order to calculate the correlative of each node, use following iterative equation:

$y_i\[j-1\]=y_i\[j\]-\frac{h}{6}(a\[i\]+2b\[i\]+2c\[i\]+d\[i\]),(i=1,2,3,4;j=1,2,3,4...n)---\left(16\right)$

For using above equation, need to estimate a [i], b [i], c [i] and d [i]:

$\left\{\begin{array}{c}a\[i\]=F_i(y_1\[n\],y_2\[n\],y_3\[n\],y_4\[n\])\\ b\[i\]=F_i(y_1\[n\]+\frac{h}{2}a\[1\],y_2\[n\]+\frac{h}{2}a\[2\],y_3\[n\]+\frac{h}{2}a\[3\],y_4\[n\]+\frac{h}{2}a\[4\])\\ c\[i\]=F_i(y_1\[n\]+\frac{h}{2}b\[1\],y_2\[n\]+\frac{h}{2}b\[2\],y_3\[n\]+\frac{h}{2}b\[3\],y_4\[n\]+\frac{h}{2}b\[4\])\\ d\[i\]=F_i(y_1\[n\]+hc\[1\],y_2\[n\]+hc\[2\],y_3\[n\]+hc\[3\],y_4\[n\]+hc\[4\])\end{array}\right.---\left(17\right).$

6. the method for numerical simulation of gas-producing well temperature, pressure distribution as claimed in claim 1, is characterized in that, in step D, described simulation softward comprises subscriber interface module, exports output module, algoritic module and chart module;

Described subscriber interface module is for providing the initial information of gas well, and these information can be divided into static data, fluid data and production data:

Static data: describe static lower gas well feature;

Fluid data: the state describing fluid;

Production data: the parameter that the production phase is given;

Described input/output module is used for reading in from storer, preserving data;

Described algoritic module is used for carrying out mathematical computations according to the algorithm of design;

Described chart module is used for the data result of simulation to carry out pictorialization displaying.