CN104866681A - Temperature pressure numerical simulation method in closing process of high-temperature high-pressure oil gas inclined shaft - Google Patents

Abstract

The invention belongs to the technical field of management of oil gas reservoir exploitation engineering and discloses a temperature pressure numerical simulation method in a closing process of a high-temperature high-pressure oil gas inclined shaft, which is used for precise prediction on temperature and pressure distribution of a shaft flow. The simulation method comprises the following steps: A, establishing a shaft-closing differential equation coupling model; B, adopting a limited differential algorithm to solve the coupling model. According to the invention, the temperature and pressure distribution of the shaft flow is precisely predicted, design level of the oil gas exploitation equipment is greatly improved and the exploitation of the oil gas reservoir is facilitated.

Description

Temperature, pressure method for numerical simulation in High Temperature High Pressure oil gas inclined shaft closing well process

Technical field

The invention belongs to development of oil and gas reservoir administrative skill field, be specially temperature, pressure method for numerical simulation in a kind of High Temperature High Pressure oil gas inclined shaft closing well process, it relates in High Temperature High Pressure Oil/gas Well inclined shaft closing well process, temperature, the coupling pressure Analysis on Mechanism of oil gas two-phase transient state stream, Mathematical Models, numerical simulation algorithm design etc.

Background technology

In Oil-Gas Well Engineering, closing well is important and the control problem of complexity.May blowout out of control be caused when the stratum, bottom of well forms the too low or overflow of declining, pass well operations must be implemented.The operation of oil well closing well will cause the change of wellhead pressure, and this may constitute a serious threat to wellhead equipment, and then causes serious consequence.

In closing well process, understand fully that the state of uncertain stream is very necessary for the effective control measure of enforcement.Specify temperature, pressure distribution, especially well head maximum pressure ensures safety for the suitable wellhead equipment of selection and has great importance.Therefore, the raising that accurately predicting greatly can promote oil gas production equipment design level is carried out to the temperature and pressure distribution of well stream.

Inclined shaft two-phase flow pressure and temp is coupled: before Modling model, be first described below basic assumption:

(1) the hydrocarbon flow flow direction is in the wellbore simply, and the heat between two-phase and the conversion of energy are ignored.

(2) air-flow and oil stream have same pressure, temperature and flow velocity at arbitrary xsect of pit shaft.

(3) friction force between gentle of pasta is ignored.

(4) vertical formation temperature is by based on known geothermic gradient linear distribution.

In closing well process, under liquid flow occurs in the condition of turbulent flow appearance.As shown in Figure 1, pit shaft pitch angle is θ to example slant well system, and cross-sectional area is A, and hydraulic diameter is d, and length is Z.The direction that the distance of coordinate flows along oil pipe.In closing well process, the hoistway of fluid becomes Transient Flow from steady flow, the state of last relative quiescent.

Material balance model: take shaft bottom as initial point, direction is forward straight up.A clock infinitesimal and length infinitesimal is respectively dt and dz.To every one phase flow service property (quality) conservation theorem.In time dt, the quality flowing through dz is the quality variable caused due to tubing string distortion and the compression of fluid is $\frac{\∂}{\∂t}\left({\mathrm{\ρ}}_k{H}_k\mathrm{Adz}\right)\mathrm{dt}.$

${\mathrm{\ρ}}_k{\mathrm{\υ}}_k{H}_k\mathrm{Adt}-[{\mathrm{\ρ}}_k{\mathrm{\υ}}_k{H}_{k}A+\frac{\∂}{\∂z}\left({\mathrm{\ρ}}_k{\mathrm{\υ}}_k{H}_{k}A\right)]\mathrm{dt}=\frac{\∂}{\∂t}\left({\mathrm{\ρ}}_k{H}_k\mathrm{Adz}\right)\mathrm{dt}---\left(1\right)$

Wherein ρ represents density, and v represents speed, and H express liquid is detained, and subscript represents gas phase or liquid phase.Above-mentioned formula (1) be can be exchanged into:

$\frac{\∂\left({\mathrm{\ρ}}_k{H}_k\right)}{\∂t}+\frac{\∂\left({\mathrm{\ρ}}_k{H}_k{v}_k\right)}{\∂z}=0---\left(2\right)$

So, liquid phase is had

$\frac{\∂\left({\mathrm{\ρ}}_l{H}_l\right)}{\∂t}+\frac{\∂\left({\mathrm{\ρ}}_l{H}_l{v}_l\right)}{\∂z}=0---\left(3\right)$

Gas phase is had

$\frac{\∂\left({\mathrm{\ρ}}_g{H}_g\right)}{\∂t}+\frac{\∂\left({\mathrm{\ρ}}_g{H}_g{v}_g\right)}{\∂z}=0---\left(4\right)$

Momentum balance: in time dt, the momentum flowing into dz is the momentum flowing out dz is the power acted on dz comprises the pressure on xsect gravity ρ _kg cos θ Adz, friction force τ _kbs _kbdz, interfacial stress τ _kjs _kjdz.In addition, the momentum increment in time dt is had by the law of conservation of momentum

$\frac{\∂\left({\mathrm{\ρ}}_k{H}_k{\mathrm{\υ}}_k\right)}{\∂t}A+\frac{\∂\left({\mathrm{\ρ}}_k{H}_k{\mathrm{\υ}}_k^2\right)}{\∂z}A+\frac{\∂\left(H_{k}P\right)}{\∂z}A=-{\mathrm{\ρ}}_{k}g\cos \mathrm{\θA}-{\mathrm{\τ}}_{\mathrm{kb}}{S}_{\mathrm{kb}}+{\mathrm{\τ}}_{\mathrm{kj}}{S}_{\mathrm{kj}}---\left(5\right)$

Be equivalent to

${\mathrm{\ρ}}_k{H}_k\frac{\∂{\mathrm{\υ}}_k}{\∂t}+H_k{\mathrm{\υ}}_k\frac{\∂{\mathrm{\ρ}}_k}{\∂t}+H_k\frac{\∂P}{\∂z}+H_k{\mathrm{\υ}}_k^2\frac{\∂{\mathrm{\ρ}}_k}{\∂z}+2{\mathrm{\ρ}}_k{H}_k{\mathrm{\υ}}_k\frac{\∂{\mathrm{\υ}}_k}{\∂z}+{\mathrm{\ρ}}_{k}g\cos \mathrm{\θ}+\frac{{\mathrm{\τ}}_{\mathrm{kb}}{S}_{\mathrm{kb}}}{A}-\frac{{\mathrm{\τ}}_{\mathrm{kj}}{S}_{\mathrm{kj}}}{A}=0---\left(6\right)$

Write as liquid phase respectively

${\mathrm{\ρ}}_l{H}_l\frac{\∂{\mathrm{\υ}}_l}{\∂t}+H_l{\mathrm{\υ}}_l\frac{{\∂\mathrm{\ρ}}_l}{\∂t}{+H}_l\frac{\∂P}{\∂z}+H_l{\mathrm{\υ}}_l^2\frac{{\∂\mathrm{\ρ}}_l}{\∂z}+{2\mathrm{\ρ}}_l{H}_l{\mathrm{\υ}}_l\frac{{\∂\mathrm{\υ}}_l}{\∂z}+{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ}+\frac{{\mathrm{\τ}}_{\mathrm{lb}}{S}_{\mathrm{lb}}}{A}-\frac{{\mathrm{\τ}}_{\lg }{S}_{\mathrm{lb}}}{A}=0---\left(7\right)$

And gas phase

${\mathrm{\ρ}}_g{H}_g\frac{\∂{\mathrm{\υ}}_g}{\∂t}+H_g{\mathrm{\υ}}_g\frac{{\∂\mathrm{\ρ}}_g}{\∂t}{+H}_g\frac{\∂P}{\∂z}+H_g{\mathrm{\υ}}_g^2\frac{{\∂\mathrm{\ρ}}_g}{\∂z}+{2\mathrm{\ρ}}_g{H}_g{\mathrm{\υ}}_g\frac{{\∂\mathrm{\υ}}_g}{\∂z}+{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ}+\frac{{\mathrm{\τ}}_{\mathrm{gb}}{S}_{\mathrm{gb}}}{A}-\frac{{\mathrm{\τ}}_{\mathrm{gl}}{S}_{\mathrm{gl}}}{A}=0---\left(8\right)$

Here, τ _lbs _lband τ _gbs _gbrepresent liquid level gentle friction stree respectively.τ _lgs _lgthe shear stress between two-phase.Bubble flow surround by liquid phase, i.e. τ _gbs _gb=0.Meanwhile, the shear stress on interface between two-phase and mass exchange are ignored, then τ _lgs _lg=τ _gls _gl=0.Now, (7) and (8) can be rewritten as

${\mathrm{\ρ}}_l{H}_l\frac{\∂{\mathrm{\υ}}_l}{\∂t}+H_l{\mathrm{\υ}}_l\frac{{\∂\mathrm{\ρ}}_l}{\∂t}{+H}_l\frac{\∂P}{\∂z}+H_l{\mathrm{\υ}}_l^2\frac{{\∂\mathrm{\ρ}}_l}{\∂z}+{2\mathrm{\ρ}}_l{H}_l{\mathrm{\υ}}_l\frac{{\∂\mathrm{\υ}}_l}{\∂z}+{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ}+\frac{{\mathrm{\τ}}_{\mathrm{lb}}{S}_{\mathrm{lb}}}{A}=0---\left(9\right)$

And

${\mathrm{\ρ}}_g{H}_g\frac{\∂{\mathrm{\υ}}_g}{\∂t}+H_g{\mathrm{\υ}}_g\frac{{\∂\mathrm{\ρ}}_g}{\∂t}{+H}_g\frac{\∂P}{\∂z}+H_g{\mathrm{\υ}}_g^2\frac{{\∂\mathrm{\ρ}}_g}{\∂z}+{2\mathrm{\ρ}}_g{H}_g{\mathrm{\υ}}_g\frac{{\∂\mathrm{\υ}}_g}{\∂z}+{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ}=0---\left(10\right)$

Energy equilibrium: as shown in Figure 2, in time dt, the energy flowing into dz is infinitesimal dz the energy flowed out is $m_k{C}_{\mathrm{pk}}T(z+\mathrm{dz})+\frac{1}{2}{m}_k{\mathrm{\υ}}_k^2(z+\mathrm{dz})+m_{k}g\cos \mathrm{\θ}(z+\mathrm{dz});$ The heat of conduct radiation is Q; In addition, in time dt, the increment of energy is wherein C _pkfor specific heat capacity.Had by principle of conservation of energy:

$\begin{array}{c}{m}_k{C}_{\mathrm{pk}}T\left(z\right)+\frac{1}{2}{m}_k{\mathrm{\υ}}_k^2\left(z\right)+m_{k}g\cos \mathrm{\θ}\left(z\right)-m_k{C}_{\mathrm{pk}}T(z+\mathrm{dz})-\frac{1}{2}{m}_k{\mathrm{\υ}}_k^2(z+\mathrm{dz})-m_{k}g\cos \mathrm{\θ}(z+\mathrm{dz})-Q\\ =\frac{\∂}{\∂t}\left({\mathrm{\ρ}}_k{H}_k{C}_{\mathrm{pk}}T\right)\mathrm{Adz}\end{array}---\left(11\right)$

Have further

$\frac{\∂}{\∂t}\left({\mathrm{\ρ}}_k{H}_k{C}_{\mathrm{pk}}T\right)A+\frac{\∂}{\∂z}{m}_k{C}_{\mathrm{pk}}T+\frac{1}{2}\frac{\∂}{\∂z}{m}_k{\mathrm{\υ}}_k^2+\frac{\∂}{\∂z}{m}_{k}g\cos \mathrm{\θ}+Q=0---\left(12\right)$

In dt, the heat radiation conduction between fluid and stratum is 2 π r _tou _to(T-T _h) H _kdz.Heat radiation conduction between stratum is [2 π K _e(T _h-T _e) H _k/ f (t _d)] dz.So total heat radiation conduction is

$Q=H_k\frac{2\mathrm{\π}{r}_{\mathrm{to}}{U}_{\mathrm{to}}{k}_e(T-T_e)}{r_{\mathrm{to}}{U}_{\mathrm{to}}f\left(t_D\right)+K_e}\mathrm{dz}---\left(13\right)$

In addition,

m _k＝ρ _kH _kυ _kA (14)

Make a=(2 π r _tou _tok _e)/A (r _tou _tof (t _d)+K _e); By (12), (13) and (14), we obtain gas phase energy conservation equation

${\mathrm{\ρ}}_g{C}_{\mathrm{pg}}\frac{\∂T}{\∂t}+{\mathrm{\ρ}}_g{C}_{\mathrm{pg}}{\mathrm{\υ}}_g\frac{\∂T}{\∂z}+C_{\mathrm{pg}}T\frac{\∂{\mathrm{\ρ}}_g}{\∂t}+({\mathrm{\ρ}}_g{C}_{\mathrm{pg}}T+\frac{3}{2}{\mathrm{\ρ}}_g{\mathrm{\υ}}_g^2+{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ})\frac{\∂{\mathrm{\υ}}_g}{\∂z}+(\frac{1}{2}{\mathrm{\υ}}_g^3+g{\mathrm{\υ}}_g\cos \mathrm{\θ})\frac{\∂{\mathrm{\ρ}}_g}{\∂z}+a(T-T_e)=0---\left(15\right)$

Similarly, liquid phase energy conservation equation can be obtained

${\mathrm{\ρ}}_1{C}_{\mathrm{pl}}\frac{\∂T}{\∂t}+{\mathrm{\ρ}}_l{C}_{\mathrm{pl}}{\mathrm{\υ}}_l\frac{\∂T}{\∂z}+({\mathrm{\ρ}}_l{C}_{\mathrm{pl}}T+\frac{3}{2}{\mathrm{\ρ}}_l{\mathrm{\υ}}_l^2+{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ})\frac{\∂{\mathrm{\υ}}_l}{\∂z}+a(T-T_e)=0---\left(16\right)$

Stratum heat transfer: stratum heat conduction model is

$\frac{\∂T_e}{\∂t}=\frac{{\mathrm{\λ}}_e}{C_{P_e}{\mathrm{\ρ}}_e}(\frac{{\∂}^2{T}_e}{\∂r^2}+\frac{1}{r}\frac{\∂T_e}{\∂r})---\left(\text{17}\right)$

Starting condition: T (0)=T ₀, 0≤z≤h,

Boundary condition: $\frac{\∂T_e}{\∂r}=0,\mathrm{if\; r}\→\∞;\mathrm{dQ}=-2\mathrm{\π}{r}_{\mathrm{cem}}{\mathrm{\λ}}_{\mathrm{cem}}\frac{\∂T_e}{\∂r}\mathrm{dz}{|}_{r=r_{\mathrm{cem}}}$

Carry out nondimensionalization process, r _d=r/r _cemand formula (17) is equivalent to

$\frac{\∂T_e}{\∂t_D}=\frac{{\mathrm{\λ}}_e}{C_{P_e}{\mathrm{\ρ}}_e}(\frac{{\∂}^2{T}_e}{\∂r_D^2}+\frac{1}{r_D}\frac{\∂T_e}{\∂r_D})---\left(\text{18}\right)$

Boundary condition can be converted into

$\frac{\∂T_e}{\∂r_D}{|}_{r_D=1}=-\frac{\mathrm{dQ}}{\mathrm{dz}}{\left(2\mathrm{\π}{\mathrm{\λ}}_e\right)}^{-1},\frac{\∂T_e}{\∂r_D}{|}_{r_D\→\∞}=0.---\left(19\right)$

Summary of the invention

Technical matters to be solved by this invention is: propose temperature, pressure method for numerical simulation in a kind of High Temperature High Pressure oil gas inclined shaft closing well process, carries out accurately predicting to the temperature and pressure distribution of well stream.

The technical solution adopted for the present invention to solve the technical problems is: temperature, pressure method for numerical simulation in High Temperature High Pressure oil gas inclined shaft closing well process, comprises the following steps:

A. closing well differential equation coupling model is set up;

B. finite differential algorithm is adopted to solve coupling model.

Further, the differential equation of closing well described in steps A coupling model is:

$\left\{\begin{array}{c}\frac{\∂\left({\mathrm{\ρ}}_l{H}_l\right)}{\∂t}+\frac{\∂\left({\mathrm{\ρ}}_l{H}_l{\mathrm{\υ}}_l\right)}{\∂z}=0\\ \frac{\∂\left({\mathrm{\ρ}}_g{H}_g\right)}{\∂t}+\frac{\∂\left({\mathrm{\ρ}}_g{H}_g{\mathrm{\υ}}_g\right)}{\∂z}=0\\ {\mathrm{\ρ}}_l{H}_l\frac{\∂{\mathrm{\υ}}_l}{\∂t}+H_l{\mathrm{\υ}}_l\frac{\∂{\mathrm{\ρ}}_l}{\∂t}+H_l\frac{\∂P}{\∂z}+H_l{\mathrm{\υ}}_l^2\frac{\∂{\mathrm{\ρ}}_l}{\∂z}+2{\mathrm{\ρ}}_l{H}_l{\mathrm{\υ}}_l\frac{\∂{\mathrm{\υ}}_l}{\∂z}+{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ}+\frac{{\mathrm{\τ}}_{\mathrm{lb}}{S}_{\mathrm{lb}}}{A}=0\\ {\mathrm{\ρ}}_g{H}_g\frac{\∂{\mathrm{\υ}}_g}{\∂t}+H_g{\mathrm{\υ}}_g\frac{\∂{\mathrm{\ρ}}_g}{\text{\∂t}}+H_g\frac{\∂P}{\∂z}+H_g{\mathrm{\υ}}_g^2\frac{\∂{\mathrm{\ρ}}_g}{\∂z}+2{\mathrm{\ρ}}_g{H}_g{\mathrm{\υ}}_g\frac{\∂{\mathrm{\υ}}_g}{\∂z}+{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ}=0\\ {\mathrm{\ρ}}_l{C}_{\mathrm{pl}}\frac{\∂T}{\∂t}+{\mathrm{\ρ}}_l{C}_{\mathrm{pl}}{\mathrm{\υ}}_l\frac{\∂T}{\∂z}+({\mathrm{\ρ}}_l{C}_{\mathrm{pl}}T+\frac{3}{2}{\mathrm{\ρ}}_l{\mathrm{\υ}}_l^2+{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ})\frac{\∂{\mathrm{\υ}}_l}{\∂z}+a(T-T_e)=0\\ {\mathrm{\ρ}}_g{C}_{\mathrm{pg}}\frac{\∂T}{\∂t}+{\mathrm{\ρ}}_g{C}_{\mathrm{pg}}{\mathrm{\υ}}_g\frac{\∂T}{\∂z}+C_{\mathrm{pg}}T\frac{\∂{\mathrm{\ρ}}_g}{\∂t}+({\mathrm{\ρ}}_g{C}_{\mathrm{pg}}T+\frac{3}{2}{\mathrm{\ρ}}_g{\mathrm{\υ}}_g^2+{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ})\frac{\∂{\mathrm{\υ}}_g}{\∂z}+(\frac{1}{2}{\mathrm{\υ}}_g^3+g{\mathrm{\υ}}_g\cos \mathrm{\θ})\frac{\∂{\mathrm{\ρ}}_g}{\∂z}+a(T-T_e)=0\\ \frac{\∂T_e}{\∂t_D}=\frac{{\mathrm{\λ}}_e}{C_{P_e}{\mathrm{\ρ}}_e}(\frac{{\∂}^2{T}_e}{\∂r_D^2}+\frac{1}{r_D}\frac{\∂T_e}{\∂r_D}),\frac{\∂T_e}{\∂r_D}{|}_{r_D=1}=-\frac{\mathrm{dQ}}{\mathrm{dz}}{\left(2\mathrm{\π}{\mathrm{\λ}}_e\right)}^{-1},\frac{\∂T_e}{\∂r_D}{|}_{r_D\→\∞}=0\\ H_l+H_g=1,C_{\mathrm{pl}}=\frac{{\mathrm{\ρ}}_o{H}_o{C}_{\mathrm{po}}+{\mathrm{\ρ}}_w{H}_w{C}_{\mathrm{pw}}}{{\mathrm{\ρ}}_l{H}_l}\\ P(z,0)=P_0\left(z\right),T(z,0)=T_0\left(z\right),\mathrm{\υ}(z,0)={\mathrm{\υ}}_0\left(z\right),T_e(r_D,0)=T_e^0\left(r_D\right),0\≤0\≤Z\\ P(0,t)={\stackrel{\‾}{P}}_0\left(t\right),T(0,t)={\stackrel{\‾}{T}}_0\left(t\right),\mathrm{\υ}(0,t)={\stackrel{\‾}{\mathrm{\υ}}}_0\left(t\right),T_e(0,t)=T_e^0\left(t\right)\end{array}\right.---\left(20\right)$

Based on above-mentioned coupling model, (P, T) forecast of distribution model adopts following matrix representation:

$A\frac{\∂\mathrm{\ψ}}{\∂t}+B\frac{\∂\mathrm{\ψ}}{\∂z}=C---\left(21\right)$

Wherein A and B is flow variables matrix, and C is the vector comprising all algebraic terms, and ψ is solution vector; Can obtain

$\begin{array}{cc}A=\left[\begin{array}{cc}0& {\mathrm{\ρ}}_g{C}_{\mathrm{pg}}\\ 0& {\mathrm{\ρ}}_l{C}_{\mathrm{pl}}\end{array}\right]& B=\left[\begin{array}{cc}{H}_g& {\mathrm{\ρ}}_g{C}_{\mathrm{pg}}{\mathrm{\υ}}_g\\ H_l& {\mathrm{\ρ}}_l{C}_{\mathrm{pl}}{\mathrm{\υ}}_l\end{array}\right]\end{array}$

$C=\left[\begin{array}{c}-{\mathrm{\ρ}}_g{H}_g\frac{\∂{\mathrm{\υ}}_g}{\∂t}-(2{\mathrm{\ρ}}_g{H}_g{\mathrm{\υ}}_g+{\mathrm{\ρ}}_g{C}_{\mathrm{pg}}T+\frac{3}{2}{\mathrm{\ρ}}_g{\mathrm{\υ}}_g^2+{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ})\frac{\∂{\mathrm{\υ}}_g}{\∂z}-(H_g{\mathrm{\υ}}_g+C_{\mathrm{pg}}T)\frac{\∂{\mathrm{\ρ}}_g}{\∂t}\\ -(H_g{\mathrm{\υ}}_g^2+\frac{1}{2}{\mathrm{\υ}}_g^3+g{\mathrm{\υ}}_g\cos \mathrm{\θ})\frac{\∂{\mathrm{\ρ}}_g}{\∂z}-{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ}-a(T-T_e)\\ -{\mathrm{\ρ}}_l{H}_l\frac{\∂{\mathrm{\υ}}_l}{\∂t}-(2{\mathrm{\ρ}}_l{H}_l{\mathrm{\υ}}_l+{\mathrm{\ρ}}_l{C}_{\mathrm{pl}}T+\frac{3}{2}{\mathrm{\ρ}}_l{\mathrm{\υ}}_l^2+{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ})\frac{\∂{\mathrm{\υ}}_l}{\∂z}-H_l{\mathrm{\υ}}_l\frac{\∂{\mathrm{\ρ}}_l}{\∂t}-H_l{\mathrm{\υ}}_l^2\frac{\∂{\mathrm{\ρ}}_l}{\∂z}\\ -{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ}-a(T-T_e)-\frac{{\mathrm{\τ}}_{\mathrm{lb}}{S}_{\mathrm{lb}}}{A}\end{array}\right]$

Variable (υ, ρ) parametric equation is as follows:

$\left\{\begin{array}{c}{\mathrm{\ρ}}_l\frac{\∂{\mathrm{\υ}}_l}{\∂t}-{\mathrm{\ρ}}_g\frac{\∂{\mathrm{\υ}}_g}{\∂t}+{\mathrm{\υ}}_l\frac{\∂{\mathrm{\ρ}}_l}{\∂t}-{\mathrm{\υ}}_g\frac{\∂{\mathrm{\ρ}}_g}{\∂t}+{\mathrm{\υ}}_l^2\frac{\∂{\mathrm{\ρ}}_l}{\∂z}-{\mathrm{\υ}}_g^2\frac{\∂{\mathrm{\ρ}}_g}{\∂z}+2{\mathrm{\ρ}}_l{\mathrm{\υ}}_l\frac{\∂{\mathrm{\υ}}_l}{\∂z}-2{\mathrm{\ρ}}_g{\mathrm{\υ}}_g\frac{\∂{\mathrm{\υ}}_g}{\∂z}+\frac{{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ}}{H_l}-\frac{{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ}}{H_g}+\frac{{\mathrm{\τ}}_{\mathrm{lb}}{S}_{\mathrm{lb}}}{{\mathrm{AH}}_l}=0\\ H_l\frac{\∂{\mathrm{\ρ}}_l}{\∂t}+{\mathrm{\ρ}}_l{H}_l\frac{\∂{\mathrm{\υ}}_l}{\∂z}+{\mathrm{\υ}}_l{H}_l\frac{\∂{\mathrm{\ρ}}_l}{\∂z}=0\\ H_g\frac{\∂{\mathrm{\ρ}}_g}{\∂t}+{\mathrm{\ρ}}_g{H}_g\frac{\∂{\mathrm{\υ}}_g}{\∂z}+{\mathrm{\υ}}_g{H}_g\frac{\∂{\mathrm{\ρ}}_g}{\∂z}=0\\ \frac{\∂T_e}{\∂t_D}=\frac{{\mathrm{\λ}}_e}{C_{P_e}{\mathrm{\ρ}}_e}(\frac{{\∂}^2{T}_e}{\∂r_D^2}+\frac{1}{r_D}\frac{\∂T_e}{\∂r_D}),\frac{\∂T_e}{\∂r_D}{|}_{r_D=1}=-\frac{\mathrm{dQ}}{\mathrm{dz}}{\left(2\mathrm{\π}{\mathrm{\λ}}_e\right)}^{-1},\frac{\∂T_e}{\∂r_D}{|}_{r_D\→\∞}=0\\ H_l+H_g=1,C_{\mathrm{pl}}=\frac{{\mathrm{\ρ}}_o{H}_o{C}_{\mathrm{po}}+{\mathrm{\ρ}}_w{H}_w{C}_{\mathrm{pw}}}{{\mathrm{\ρ}}_l{H}_l}\\ P(z,0)=P_0\left(z\right),T(z,0)=T_0\left(z\right),\mathrm{\υ}(z,0)={\mathrm{\υ}}_0\left(z\right),T_e(r_D,0)=T_e^0\left(r_D\right),0\≤z\≤Z\\ P(0,t)={\stackrel{\‾}{P}}_0\left(t\right),T(0,t)={\stackrel{\‾}{T}}_0\left(t\right),\mathrm{\υ}(0,t)={\stackrel{\‾}{\mathrm{\υ}}}_0\left(t\right),T_e(0,t)=T_e^0\left(t\right)\end{array}\right.---\left(22\right).$

Further, in step B, adopt finite differential algorithm to solve coupling model and comprise:

B1. finite differential framework is set up;

B2. design conditions parameter;

B3. based on finite differential framework, coupling model is solved.

Further, set up finite differential framework described in step B1 to comprise:

By F (z _j, t _k) be designated as wherein, z _j=jh, j=1,2 ..., m, t _k=k τ, k=1,2 ..., n

At point (z _j+1/2, t _k) place pair carry out discretize respectively with U, can obtain:

$\left\{\begin{array}{c}{\left(\frac{\∂U}{\∂t}\right)}_{j+1/2}^{k+1}=\frac{U_{j+1}^{k+1}+U_j^{k+1}-U_{j+1}^k-U_j^k}{2\mathrm{\τ}}\\ {\left(\frac{\∂U}{\∂z}\right)}_{j+1/2}^{k+1}=\frac{U_{j+1}^{k+1}-U_j^{k+1}}{h}\\ U_{j+1/2}^{k+1}=\frac{U_{j+1}^{k+1}+U_j^{k+1}}{2}\end{array}\right.---\left(23\right).$

Further, in step B2, described design conditions parameter comprises:

(a) heat diffusion equation:

$f\left(t_D\right)=\left\{\begin{array}{c}1.128\sqrt{t_D}(1-0.3\sqrt{t_D})t_D\≤1.5\\ (0.4063+0.5\ln t_D)(1+0.6/t_D)t_D>1.5\end{array}\right.$

Wherein t _d=at/r ² _wh;

B () pit shaft diverse location is to the heat-conduction coefficient U of second contact surface _to:

$\frac{1}{U_{\mathrm{to}}}=r_{\mathrm{ti}}\frac{1}{{\mathrm{\λ}}_{\mathrm{ins}}}\ln \left(\frac{r_{\mathrm{ci}}}{r_{\mathrm{to}}}\right)+\frac{1}{h_c+h_r}+r_{\mathrm{ti}}\frac{1}{{\mathrm{\λ}}_{\mathrm{cem}}}\ln \left(\frac{r_{\mathrm{cem}}}{r_{\mathrm{co}}}\right);$

Friction force between (c) fluid and tube wall:

${\mathrm{\τ}}_{\mathrm{lb}}{S}_{\mathrm{lb}}/A=f{\mathrm{\ρ}}_l{\mathrm{\υ}}_l^2/\left(2d\right);$

(d) friction factor f:

$1/\sqrt{f}=1.14-21g(e/d+21.25/{\mathrm{Re}}^{0.9})$

E () gas-liquid mixed ratio of specific heat is held:

C _Pg＝1697.5107P ^0.0661T ^0.0776，C _Pl＝4.2KJ/(Kg·℃)。

Further, in step B3, described solving coupling model based on finite differential framework comprises:

B31: the step-length of setting-up time and the degree of depth, obtains the cant angle theta of every bit _j=θ _j-1+ (θ _k-θ _k-1) Δ s _j/ Δ s _k, wherein j represents calculating infinitesimal, s _krepresent to tilt for θ _k; Time measure the degree of depth;

B32: specify suitable finite differential grid;

B33: given starting condition and boundary condition;

B34: make T=T _k, solve the T that following system of equations obtains time t _e:

$\left\{\begin{array}{c}\frac{\∂T_e}{\∂t_D}=\frac{{\mathrm{\λ}}_e}{C_{P_e}{\mathrm{\ρ}}_e}\frac{{\∂}^2{T}_e}{\∂r_D^2}+\frac{{\mathrm{\λ}}_e}{r_D{C}_{P_e}{\mathrm{\ρ}}_e}\frac{\∂T_e}{\∂r_D}\\ \frac{\∂T_e}{\∂r_D}{|}_{r_D=1}=-\mathrm{\α}(T-T_e){\left(2\mathrm{\π}{\mathrm{\λ}}_e\right)}^{-1}\\ \frac{\∂T_e}{\∂r_D}{|}_{r_D\→\∞}=0\\ T_e\left(0\right)=T_e^0\end{array}\right.---\left(24\right)$

Order for the temperature of time t, degree of depth z, i=1,2 ..., t _s, j=1,2 ..., N, wherein t _stime and last radial point is represented respectively with N, it is the initial temperature on stratum; Finite differencing method is applied to formula (24) to obtain

τ and represent respectively between time interval and radial portion; Canonical form can be converted into

Differential method is used to carry out sliding-model control to boundary condition; To r _d=1, have

$\frac{\∂T_e}{\∂r_D}{|}_{r_D=1}-\frac{\mathrm{\α}}{2\mathrm{\π}{\mathrm{\λ}}_e}T{|}_{r_D=1}=-\frac{\mathrm{\α}{T}_k}{2\mathrm{\π}{\mathrm{\λ}}_e}---\left(27\right)$

So

To r _d=N, has

$T_{e,N}^{i+1}-T_{e,N-1}^{i+1}=0---\left(29\right)$

Integrate (26), (28) and (29), the calculably symbolic solution of layer temperature, can obtain T _ediscretize matrix $\left[T_{e,j}^i\right];$

B35: use 3-diagonal matrix Algorithm for Solving coupling model;

B36: along time dimension, performs B34 and B35, until the time terminates.

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

Accompanying drawing explanation

Fig. 1 is pit shaft schematic diagram;

Fig. 2 is heat radiation schematic diagram;

Fig. 3 is the discretize schematic diagram in finite solution region;

Fig. 4 (a) for pressure in closing well process in time with change in depth curve map, Fig. 4 (b) is wellhead pressure changing trend diagram in time;

Fig. 5 (a) for temperature in closing well process in time with change in depth curve map, Fig. 4 (b) is wellhead temperature changing trend diagram in time;

Fig. 6 (a) for gas velocity in closing well process in time with change in depth curve map, Fig. 6 (b) is well head gas velocity changing trend diagram in time;

Fig. 7 (a) for liquid velocity in closing well process in time with change in depth curve map, Fig. 7 (b) is well head liquid velocity changing trend diagram in time;

Fig. 8 (a) for gas density in closing well process in time with change in depth curve map, Fig. 8 (b) is well head gas volume density changing trend diagram in time;

Fig. 9 (a) is the pressure distribution schematic diagram under different output; Fig. 9 (b) is the Temperature Distribution schematic diagram under different output;

Figure 10 (a) is for resulting fluid speed is with change of production schematic diagram; Figure 10 (b) is for resulting fluid density is with change of production schematic diagram.

Embodiment

The present invention is intended to propose temperature, pressure method for numerical simulation in a kind of High Temperature High Pressure oil gas inclined shaft closing well process, carries out accurately predicting to the temperature and pressure distribution of well stream.

On concrete enforcement, in the High Temperature High Pressure oil gas inclined shaft closing well process in the present invention, temperature, pressure method for numerical simulation comprises following implementation:

Closing well differential equation coupling model: by (3), (4), (9), (10), (15) and (16) can obtain the derivative scalar coupling in closing well process:

$\left\{\begin{array}{c}\frac{\∂\left({\mathrm{\ρ}}_l{H}_l\right)}{\∂t}+\frac{\∂\left({\mathrm{\ρ}}_l{H}_l{\mathrm{\υ}}_l\right)}{\∂z}=0\\ \frac{\∂\left({\mathrm{\ρ}}_g{H}_g\right)}{\∂t}+\frac{\∂\left({\mathrm{\ρ}}_g{H}_g{\mathrm{\υ}}_g\right)}{\∂z}=0\\ {\mathrm{\ρ}}_l{H}_l\frac{\∂{\mathrm{\υ}}_l}{\∂t}+H_l{\mathrm{\υ}}_l\frac{\∂{\mathrm{\ρ}}_l}{\∂t}+H_l\frac{\∂P}{\∂z}+H_l{\mathrm{\υ}}_l^2\frac{\∂{\mathrm{\ρ}}_l}{\∂z}+2{\mathrm{\ρ}}_l{H}_l{\mathrm{\υ}}_l\frac{\∂{\mathrm{\υ}}_l}{\∂z}+{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ}+\frac{{\mathrm{\τ}}_{\mathrm{lb}}{S}_{\mathrm{lb}}}{A}=0\\ {\mathrm{\ρ}}_g{H}_g\frac{\∂{\mathrm{\υ}}_g}{\∂t}+H_g{\mathrm{\υ}}_g\frac{\∂{\mathrm{\ρ}}_g}{\text{\∂t}}+H_g\frac{\∂P}{\∂z}+H_g{\mathrm{\υ}}_g^2\frac{\∂{\mathrm{\ρ}}_g}{\∂z}+2{\mathrm{\ρ}}_g{H}_g{\mathrm{\υ}}_g\frac{\∂{\mathrm{\υ}}_g}{\∂z}+{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ}=0\\ {\mathrm{\ρ}}_l{C}_{\mathrm{pl}}\frac{\∂T}{\∂t}+{\mathrm{\ρ}}_l{C}_{\mathrm{pl}}{\mathrm{\υ}}_l\frac{\∂T}{\∂z}+({\mathrm{\ρ}}_l{C}_{\mathrm{pl}}T+\frac{3}{2}{\mathrm{\ρ}}_l{\mathrm{\υ}}_l^2+{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ})\frac{\∂{\mathrm{\υ}}_l}{\∂z}+a(T-T_e)=0\\ {\mathrm{\ρ}}_g{C}_{\mathrm{pg}}\frac{\∂T}{\∂t}+{\mathrm{\ρ}}_g{C}_{\mathrm{pg}}{\mathrm{\υ}}_g\frac{\∂T}{\∂z}+C_{\mathrm{pg}}T\frac{\∂{\mathrm{\ρ}}_g}{\∂t}+({\mathrm{\ρ}}_g{C}_{\mathrm{pg}}T+\frac{3}{2}{\mathrm{\ρ}}_g{\mathrm{\υ}}_g^2+{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ})\frac{\∂{\mathrm{\υ}}_g}{\∂z}+(\frac{1}{2}{\mathrm{\υ}}_g^3+g{\mathrm{\υ}}_g\cos \mathrm{\θ})\frac{\∂{\mathrm{\ρ}}_g}{\∂z}+a(T-T_e)=0\\ \frac{\∂T_e}{\∂t_D}=\frac{{\mathrm{\λ}}_e}{C_{P_e}{\mathrm{\ρ}}_e}(\frac{{\∂}^2{T}_e}{\∂r_D^2}+\frac{1}{r_D}\frac{\∂T_e}{\∂r_D}),\frac{\∂T_e}{\∂r_D}{|}_{r_D=1}=-\frac{\mathrm{dQ}}{\mathrm{dz}}{\left(2\mathrm{\π}{\mathrm{\λ}}_e\right)}^{-1},\frac{\∂T_e}{\∂r_D}{|}_{r_D\→\∞}=0\\ H_l+H_g=1,C_{\mathrm{pl}}=\frac{{\mathrm{\ρ}}_o{H}_o{C}_{\mathrm{po}}+{\mathrm{\ρ}}_w{H}_w{C}_{\mathrm{pw}}}{{\mathrm{\ρ}}_l{H}_l}\\ P(z,0)=P_0\left(z\right),T(z,0)=T_0\left(z\right),\mathrm{\υ}(z,0)={\mathrm{\υ}}_0\left(z\right),T_e(r_D,0)=T_e^0\left(r_D\right),0\≤0\≤Z\\ P(0,t)={\stackrel{\‾}{P}}_0\left(t\right),T(0,t)={\stackrel{\‾}{T}}_0\left(t\right),\mathrm{\υ}(0,t)={\stackrel{\‾}{\mathrm{\υ}}}_0\left(t\right),T_e(0,t)=T_e^0\left(t\right)\end{array}\right.---\left(20\right)$

Based on above-mentioned model, (P, T) forecast of distribution model can use following matrix representation:

$A\frac{\∂\mathrm{\ψ}}{\∂t}+B\frac{\∂\mathrm{\ψ}}{\∂z}=C---\left(21\right)$

Wherein A and B is flow variables matrix, and C is the vector comprising all algebraic terms, and ψ is solution vector; According to these agreements, can obtain

$\begin{array}{cc}A=\left[\begin{array}{cc}0& {\mathrm{\ρ}}_g{C}_{\mathrm{pg}}\\ 0& {\mathrm{\ρ}}_l{C}_{\mathrm{pl}}\end{array}\right]& B=\left[\begin{array}{cc}{H}_g& {\mathrm{\ρ}}_g{C}_{\mathrm{pg}}{\mathrm{\υ}}_g\\ H_l& {\mathrm{\ρ}}_l{C}_{\mathrm{pl}}{\mathrm{\υ}}_l\end{array}\right]\end{array}$

$C=\left[\begin{array}{c}-{\mathrm{\ρ}}_g{H}_g\frac{\∂{\mathrm{\υ}}_g}{\∂t}-(2{\mathrm{\ρ}}_g{H}_g{\mathrm{\υ}}_g+{\mathrm{\ρ}}_g{C}_{\mathrm{pg}}T+\frac{3}{2}{\mathrm{\ρ}}_g{\mathrm{\υ}}_g^2+{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ})\frac{\∂{\mathrm{\υ}}_g}{\∂z}-(H_g{\mathrm{\υ}}_g+C_{\mathrm{pg}}T)\frac{\∂{\mathrm{\ρ}}_g}{\∂t}\\ -(H_g{\mathrm{\υ}}_g^2+\frac{1}{2}{\mathrm{\υ}}_g^3+g{\mathrm{\υ}}_g\cos \mathrm{\θ})\frac{\∂{\mathrm{\ρ}}_g}{\∂z}-{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ}-a(T-T_e)\\ -{\mathrm{\ρ}}_l{H}_l\frac{\∂{\mathrm{\υ}}_l}{\∂t}-(2{\mathrm{\ρ}}_l{H}_l{\mathrm{\υ}}_l+{\mathrm{\ρ}}_l{C}_{\mathrm{pl}}T+\frac{3}{2}{\mathrm{\ρ}}_l{\mathrm{\υ}}_l^2+{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ})\frac{\∂{\mathrm{\υ}}_l}{\∂z}-H_l{\mathrm{\υ}}_l\frac{\∂{\mathrm{\ρ}}_l}{\∂t}-H_l{\mathrm{\υ}}_l^2\frac{\∂{\mathrm{\ρ}}_l}{\∂z}\\ -{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ}-a(T-T_e)-\frac{{\mathrm{\τ}}_{\mathrm{lb}}{S}_{\mathrm{lb}}}{A}\end{array}\right]$

Variable (υ, ρ) parametric equation is as follows:

$\left\{\begin{array}{c}{\mathrm{\ρ}}_l\frac{\∂{\mathrm{\υ}}_l}{\∂t}-{\mathrm{\ρ}}_g\frac{\∂{\mathrm{\υ}}_g}{\∂t}+{\mathrm{\υ}}_l\frac{\∂{\mathrm{\ρ}}_l}{\∂t}-{\mathrm{\υ}}_g\frac{\∂{\mathrm{\ρ}}_g}{\∂t}+{\mathrm{\υ}}_l^2\frac{\∂{\mathrm{\ρ}}_l}{\∂z}-{\mathrm{\υ}}_g^2\frac{\∂{\mathrm{\ρ}}_g}{\∂z}+2{\mathrm{\ρ}}_l{\mathrm{\υ}}_l\frac{\∂{\mathrm{\υ}}_l}{\∂z}-2{\mathrm{\ρ}}_g{\mathrm{\υ}}_g\frac{\∂{\mathrm{\υ}}_g}{\∂z}+\frac{{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ}}{H_l}-\frac{{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ}}{H_g}+\frac{{\mathrm{\τ}}_{\mathrm{lb}}{S}_{\mathrm{lb}}}{{\mathrm{AH}}_l}=0\\ H_l\frac{\∂{\mathrm{\ρ}}_l}{\∂t}+{\mathrm{\ρ}}_l{H}_l\frac{\∂{\mathrm{\υ}}_l}{\∂z}+{\mathrm{\υ}}_l{H}_l\frac{\∂{\mathrm{\ρ}}_l}{\∂z}=0\\ H_g\frac{\∂{\mathrm{\ρ}}_g}{\∂t}+{\mathrm{\ρ}}_g{H}_g\frac{\∂{\mathrm{\υ}}_g}{\∂z}+{\mathrm{\υ}}_g{H}_g\frac{\∂{\mathrm{\ρ}}_g}{\∂z}=0\\ \frac{\∂T_e}{\∂t_D}=\frac{{\mathrm{\λ}}_e}{C_{P_e}{\mathrm{\ρ}}_e}(\frac{{\∂}^2{T}_e}{\∂r_D^2}+\frac{1}{r_D}\frac{\∂T_e}{\∂r_D}),\frac{\∂T_e}{\∂r_D}{|}_{r_D=1}=-\frac{\mathrm{dQ}}{\mathrm{dz}}{\left(2\mathrm{\π}{\mathrm{\λ}}_e\right)}^{-1},\frac{\∂T_e}{\∂r_D}{|}_{r_D\→\∞}=0\\ H_l+H_g=1,C_{\mathrm{pl}}=\frac{{\mathrm{\ρ}}_o{H}_o{C}_{\mathrm{po}}+{\mathrm{\ρ}}_w{H}_w{C}_{\mathrm{pw}}}{{\mathrm{\ρ}}_l{H}_l}\\ P(z,0)=P_0\left(z\right),T(z,0)=T_0\left(z\right),\mathrm{\υ}(z,0)={\mathrm{\υ}}_0\left(z\right),T_e(r_D,0)=T_e^0\left(r_D\right),0\≤z\≤Z\\ P(0,t)={\stackrel{\‾}{P}}_0\left(t\right),T(0,t)={\stackrel{\‾}{T}}_0\left(t\right),\mathrm{\υ}(0,t)={\stackrel{\‾}{\mathrm{\υ}}}_0\left(t\right),T_e(0,t)=T_e^0\left(t\right)\end{array}\right.---\left(22\right)$

Finite differential algorithm design; For the sake of simplicity, some length is divided into by pit shaft to be the part of h; The deformation of every part and pipe thickness, in pipe, the density of outer fluid is relevant; Length is divided into be the decile of τ the test duration; Model calculates from moment 0 and shaft bottom.Temperature, pressure, speed and density etc. are calculated piecemeal.

(1) finite differential framework: be generally difficult to the analytic solution obtaining model (21), available values analogy method obtains approximate solution.Be applicable to planning problem due to finite differencing method and can obtain separating comparatively accurately, therefore adopting this method.The first step of finite differential determines differential structrue, replaces partial derivative by difference quotient.Thus partial differential equation is converted into algebraic equation.So solve partial differential equations to become Solving Algebraic Equation group, limited differential structrue can be built based on Taylor expansion.

For (21), three class formulas need discretize, and U.By Z-t plane, use horizontal line ot

Split, i.e. z _j=jh, j=1,2 ..., m, the longitudinal axis carries out similar segmentation, has t _k=k τ, k=1,2 ..., n, see Fig. 3.

For simplicity, by F (z _j, t _k) be designated as at point (z _j+1/2, t _k) place is to above-mentioned carry out discretize respectively with U, three differential equations can be obtained:

$\left\{\begin{array}{c}{\left(\frac{\∂U}{\∂t}\right)}_{j+1/2}^{k+1}=\frac{U_{j+1}^{k+1}+U_j^{k+1}-U_{j+1}^k-U_j^k}{2\mathrm{\τ}}\\ {\left(\frac{\∂U}{\∂z}\right)}_{j+1/2}^{k+1}=\frac{U_{j+1}^{k+1}-U_j^{k+1}}{h}\\ U_{j+1/2}^{k+1}=\frac{U_{j+1}^{k+1}+U_j^{k+1}}{2}\end{array}\right.---\left(23\right)$

(2) parameter calculates:

(a) heat diffusion equation

$f\left(t_D\right)=\left\{\begin{array}{c}1.128\sqrt{t_D}(1-0.3\sqrt{t_D})t_D\≤1.5\\ (0.4063+0.5\ln t_D)(1+0.6/t_D)t_D>1.5\end{array}\right.---\left(24\right)$

Wherein $t_D=\mathrm{\αt}/r_{\mathrm{wh}}^2.$

B () pit shaft diverse location is to the heat-conduction coefficient U of second contact surface _to

$\frac{1}{U_{\mathrm{to}}}=r_{\mathrm{ti}}\frac{1}{{\mathrm{\λ}}_{\mathrm{ins}}}\ln \left(\frac{r_{\mathrm{ci}}}{r_{\mathrm{to}}}\right)+\frac{1}{h_c+h_r}+r_{\mathrm{ti}}\frac{1}{{\mathrm{\λ}}_{\mathrm{cem}}}\ln \left(\frac{r_{\mathrm{cem}}}{r_{\mathrm{co}}}\right)$

Friction force between (c) fluid and tube wall

${\mathrm{\τ}}_{\mathrm{lb}}{S}_{\mathrm{lb}}/A=f{\mathrm{\ρ}}_l{\mathrm{\υ}}_l^2/\left(2d\right)$

(d) friction factor f

$1/\sqrt{f}=1.14-21g(e/d+21.25/{\mathrm{Re}}^{0.9})$

E () gas-liquid mixed ratio of specific heat is held

C _Pg＝1697.5107P ^0.0661T ^0.0776，C _Pl＝4.2KJ/(Kg·°C)

(3) algorithm steps design:

After differential framework is determined, the parameters such as pressure, temperature, speed and density can be calculated piecemeal.Steady production state is as the original state of closing well.

Step 1: the step-length of setting-up time and the degree of depth.Obtain the cant angle theta of every bit _j=θ _j-1+ (θ _k-θ _k-1) Δ s _j/ Δ s _k,

Wherein j represents calculating infinitesimal, s _krepresent to tilt for θ _ktime measure the degree of depth.

Step 2: specify suitable finite differential grid.

Step 3: given starting condition and boundary condition.

Step 4: make T=T _k.Solve the T that following system of equations obtains time t _e:

$\left\{\begin{array}{c}\frac{\∂T_e}{\∂t_D}=\frac{{\mathrm{\λ}}_e}{C_{P_e}{\mathrm{\ρ}}_e}\frac{{\∂}^2{T}_e}{\∂r_D^2}+\frac{{\mathrm{\λ}}_e}{r_D{C}_{P_e}{\mathrm{\ρ}}_e}\frac{\∂T_e}{\∂r_D}\\ \frac{\∂T_e}{\∂r_D}{|}_{r_D=1}=-\mathrm{\α}(T-T_e){\left(2\mathrm{\π}{\mathrm{\λ}}_e\right)}^{-1}\\ \frac{\∂T_e}{\∂r_D}{|}_{r_D\→\∞}=0\\ T_e\left(0\right)=T_e^0\end{array}\right.---\left(24\right)$

Order for the temperature of time t, degree of depth z, i=1,2 ..., t _s, j=1,2 ..., N, wherein t _stime and last radial point is represented respectively with N. it is the initial temperature on stratum.Finite differencing method is applied to (24) to obtain

τ and represent respectively between time interval and radial portion.Canonical form can be converted into

Differential method is used to carry out sliding-model control to boundary condition.To r _d=1, have

$\frac{\∂T_e}{\∂r_D}{|}_{r_D=1}-\frac{\mathrm{\α}}{2\mathrm{\π}{\mathrm{\λ}}_e}T{|}_{r_D=1}=-\frac{\mathrm{\α}{T}_k}{2\mathrm{\π}{\mathrm{\λ}}_e}---\left(27\right)$

So

To r _d=N, has

$T_{e,N}^{i+1}-T_{e,N-1}^{i+1}=0---\left(29\right)$

Integrate (26), (28) and (29), the calculably symbolic solution of layer temperature, can obtain T _ediscretize matrix $\left[T_{e,j}^i\right].$

Step 5: use 3-diagonal matrix Algorithm for Solving coupling model.

Step 6: along time dimension, performs step 4 and step 5, until the time terminates.

Embodiment:

For Y gas well at HTHP, utilize the above model set up, the parameters such as the wellhead pressure temperature in closing well process are calculated.As described in above-mentioned model analysis and solution procedure, well section is split into some perforation unit from bottom, the time also will do similar division, then calculates according to above-mentioned calculation procedure.

Correlation parameter and measurement data: well depth 3700m, crude production rate is 14m every day ³, key pressure 4.968MPa, thermal conductive coefficient 2.06, temperature at the bottom of ground 16 DEG C, underground temperature gradient 0.02 DEG C/m, factor of porosity 0.2, the roughness 0.000015 of well inwall, total compressibility 0.03, gas liquid ratio 100000:1, length section 10m, time period 10mimutes; All the other parameters are in table 1 and table 2.

Table 1

Table 2

Result of calculation is analyzed: by numerical simulation, can obtain a series of analysis result.The pressure and temp distribution of pit shaft in main concern closing well process.From Fig. 4 a, when the time is fixed, pressure increases along with the degree of depth and increases; When the degree of depth is fixed, pressure increases along with the time and increases, and reaches steady state (SS) in the some time; Along with the degree of depth increases, pressure amplification reduces.In closing well process, the fluid on stratum enters pit shaft, causes the corresponding increase of the pressure in shaft bottom.After closing well, the fluid in pit shaft is close to and stops, and the pressure wave of the generation that bottomhole wellbore pressure rises is transmitted to well head, and then causes wellhead pressure to increase.The variation tendency of wellhead pressure is shown in Fig. 4 b, and pressure is Rapid Variable Design in early days, a period of time laggard enter stationary state.

From Fig. 5 a, when the time is fixed, temperature increases along with the degree of depth and increases.When the degree of depth is fixed, temperature increases along with the time and reduces.This is because, along with the increase of time after closing well, there is the heat conduction of downhole well fluid to surrounding formations, so temperature progressively declines.Reach steady state (SS) in the some time, fluid temperature (F.T.) reaches Formation Depth according to linear distribution in each degree of depth.The variation tendency of wellhead temperature is shown in Fig. 5 b, and temperature declines in early days fast, enters closely layer temperature after a period of time.

Variable parameter is analyzed: in numerical simulation, and speed and density are regarded variable element.From Fig. 6 (a), when the time is fixed, gas velocity increases along with the degree of depth and reduces; When the degree of depth is fixed, gas velocity increases along with the time and reduces.In addition, from Fig. 6 (b), along with passage of time, well head gas velocity changes in very little scope.From Fig. 7 (a), when the time is fixed, liquid velocity increases along with the degree of depth and reduces; When the degree of depth is fixed, liquid velocity increases along with the time and reduces.As can be seen from Fig. 7 (b), along with passage of time, well head liquid velocity changes in very little scope; From Fig. 8 (a), when the time is fixed, gas density increases along with the degree of depth and slowly increases; When the degree of depth is fixed, gas velocity increases along with the time and increases.As can be seen from Fig. 8 (b), along with passage of time, well head gas volume density variation range is very little.

Sensitivity analysis: except the improper closing well caused due to improper operation, normal closing well is generally because the underproduction causes the reduction forming energy, when formation can return to certain level, needs to drive a well oil recovery.Therefore, crude production rate has material impact to parameters such as temperature, pressures.According to field data, three kinds of yield levels are used for sensitivity analysis: 6m ³/ day, 14m ³/ day, and 24m ³/ day.For studying the difference of different output, by numerical simulation algorithm, obtain series of results, the test duration is 30minutes after closing well.

Fig. 9 (a) shows the pressure distribution under different output: in the identical degree of depth, and pressure increases along with the increase of output.In fact, reduce law by the output of pit shaft, being formed can be the key factor affecting oil well output, reflects the deliverability of oil reservoir to oil well.Mineralization pressure is larger, and oil well output is larger.Hydrodynamic pressure is with the corresponding increase of increase of oil yield.When output is from 6m ³/ day is increased to 24m ³/ day, wellhead pressure is increased to 31.4Mpa by 22.7Mpa.But bottomhole wellbore pressure change is little, and other factors existing and affect output are described.

Fig. 9 (b) shows the Temperature Distribution under different output: in the identical degree of depth, and temperature raises along with the increase of output.This is because along with output increase, the corresponding increase of pressure, and then temperature raises because forming the increase of energy.When output is from 6m ³/ day is increased to 24m ³/ day, wellhead temperature is increased to 162.75 DEG C by 160.045 DEG C.

As shown in Figure 10 (a), in the identical degree of depth, fluid velocity increases along with output and increases.The maximal rate of well head oil is the maximal rate 17.5m/s of 14.4m/s, gas.If fluid velocity is excessive, friction loss increases, and this can bring infringement to the borehole wall.Therefore, from engineering viewpoint, fluid velocity should control in certain scope.

As shown in Figure 10 (b), fluid density increases along with crude production rate and reduces.When output is from 6m ³/ day is increased to 24m ³/ day, resulting fluid density is from 88.06kg/m ³become 81.58kg/m ³.

Claims (6)

1. temperature, pressure method for numerical simulation in High Temperature High Pressure oil gas inclined shaft closing well process, is characterized in that, comprise the following steps:

A. closing well differential equation coupling model is set up;

B. finite differential algorithm is adopted to solve coupling model.

2. temperature, pressure method for numerical simulation in High Temperature High Pressure oil gas inclined shaft closing well process as claimed in claim 1, it is characterized in that, the differential equation of closing well described in steps A coupling model is:

$\left\{\begin{array}{c}\frac{\∂\left({\mathrm{\ρ}}_l{H}_l\right)}{\∂t}+\frac{\∂\left({\mathrm{\ρ}}_l{H}_l{\mathrm{\υ}}_l\right)}{\∂z}=0\\ \frac{\∂\left({\mathrm{\ρ}}_g{H}_g\right)}{\∂t}+\frac{\∂\left({\mathrm{\ρ}}_g{H}_g{\mathrm{\υ}}_g\right)}{\∂z}=0\\ {\mathrm{\ρ}}_l{H}_l\frac{{\∂\mathrm{\υ}}_l}{\∂t}+H_l{\mathrm{\υ}}_l\frac{\∂{\mathrm{\ρ}}_l}{\∂t}+H_l\frac{\∂P}{\∂z}+H_l{\mathrm{\υ}}_l^2+\frac{\∂{\mathrm{\ρ}}_l}{\∂z}+2{\mathrm{\ρ}}_l{H}_l{\mathrm{\υ}}_l\frac{\∂{\mathrm{\υ}}_l}{\∂z}+{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ}+\frac{{\mathrm{\τ}}_{\mathrm{lb}}{S}_{\mathrm{lb}}}{A}=0\\ {\mathrm{\ρ}}_g{H}_g\frac{{\∂\mathrm{\υ}}_g}{\∂t}+H_g{\mathrm{\υ}}_g\frac{\∂{\mathrm{\ρ}}_g}{\∂t}+H_g\frac{\∂P}{\∂z}+H_g{\mathrm{\υ}}_g^2\frac{{\∂\mathrm{\ρ}}_g}{\∂z}+2{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ}=0\\ {\mathrm{\ρ}}_l{C}_{\mathrm{pl}}\frac{\∂T}{\∂t}+{\mathrm{\ρ}}_l{C}_{\mathrm{pl}}{\mathrm{\υ}}_1\frac{\∂T}{\∂z}+({\mathrm{\ρ}}_l{C}_{\mathrm{pl}}T+\frac{3}{2}{\mathrm{\ρ}}_l{\mathrm{\υ}}_l^2+{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ})\frac{{\∂\mathrm{\υ}}_l}{\∂z}+a\left({T-T}_e\right)=0\\ {\mathrm{\ρ}}_g{C}_{\mathrm{pg}}\frac{\∂T}{\∂t}+{\mathrm{\ρ}}_g{C}_{\mathrm{pg}}{\mathrm{\υ}}_g\frac{\∂T}{\∂z}+C_{\mathrm{pg}}T\frac{{\∂\mathrm{\ρ}}_g}{\∂t}+({\mathrm{\ρ}}_g{C}_{\mathrm{pg}}T+\frac{3}{2}{\mathrm{\ρ}}_g{\mathrm{\υ}}_g^2+{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ})\frac{{\∂\mathrm{\υ}}_g}{\∂z}+(\frac{1}{2}{\mathrm{\υ}}_g^3+{\mathrm{g\υ}}_g\cos \mathrm{\θ})\frac{{\∂\mathrm{\ρ}}_g}{\∂z}+a\left({T-T}_e\right)=0\\ \frac{\∂T_e}{\∂t_D}=\frac{{\mathrm{\λ}}_e}{C_{P_e{\mathrm{\ρ}}_e}}(\frac{{\∂}^2{T}_e}{\∂r_D^2}+\frac{1}{r_D}\frac{\∂T_e}{\∂r_D}),\frac{\∂T_e}{\∂r_D}{|}_{r_D=1}=-\frac{\mathrm{dQ}}{\mathrm{dz}}{\left(2{\mathrm{\π\λ}}_e\right)}^{-1},\frac{\∂T_e}{\∂r_D}{|}_{r_D\→\∞}=0\\ H_l+H_g=1,C_{\mathrm{pl}}=\frac{{\mathrm{\ρ}}_o{H}_o{C}_{\mathrm{po}}+{\mathrm{\ρ}}_w+H_w{C}_{\mathrm{pw}}}{{\mathrm{\ρ}}_l{H}_l}\\ P(z,0)=P_0\left(z\right),T(z,0)=T_0\left(z\right),\mathrm{\υ}(z,0)={\mathrm{\υ}}_0\left(z\right),T_e(r_D,0)=T_e^0\left(r_D\right),0\≤z\≤Z\\ P(0,t)={\stackrel{\‾}{P}}_0\left(t\right),T(0,t)={\stackrel{\‾}{T}}_0\left(t\right),\mathrm{\υ}(0,t)={\stackrel{\‾}{\mathrm{\υ}}}_0\left(t\right),T_e(0,t)=T_e^0\left(t\right)\end{array}\right.---\left(20\right)$

Based on above-mentioned coupling model, (P, T) forecast of distribution model adopts following matrix representation:

$A\frac{\∂\mathrm{\ψ}}{\∂t}+B\frac{\∂\mathrm{\ψ}}{\∂z}=C---\left(21\right)$

Wherein A and B is flow variables matrix, and C is the vector comprising all algebraic terms, and ψ is solution vector; Can obtain

$A=\left[\begin{array}{cc}0& {\mathrm{\ρ}}_g{C}_{\mathrm{pg}}\\ 0& {\mathrm{\ρ}}_l{C}_{\mathrm{pl}}\end{array}\right]$ $B=\left[\begin{array}{cc}{H}_g& {\mathrm{\ρ}}_g{C}_{\mathrm{pg}}{\mathrm{\υ}}_g\\ H_l& {\mathrm{\ρ}}_l{C}_{\mathrm{pl}}{\mathrm{\υ}}_1\end{array}\right]$

$C=\left[\begin{array}{c}-{\mathrm{\ρ}}_g{H}_g\frac{{\∂\mathrm{\υ}}_g}{\∂t}-(2{\mathrm{\ρ}}_g{H}_g{\mathrm{\υ}}_g+{\mathrm{\ρ}}_g{C}_{\mathrm{pg}}T+\frac{3}{2}{\mathrm{\ρ}}_g{\mathrm{\υ}}_g^2+{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ})\frac{\∂{\mathrm{\υ}}_g}{\∂z}-(H_g{\mathrm{\υ}}_g+C_{\mathrm{pg}}T)\frac{{\∂\mathrm{\ρ}}_g}{\∂t}\\ -(H_g{\mathrm{\υ}}_g^2+\frac{1}{2}{\mathrm{\υ}}_g^3+g{\mathrm{\υ}}_g\cos \mathrm{\θ})\frac{\∂{\mathrm{\ρ}}_g}{\∂z}-{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ}-a(T-T_e)\\ -{\mathrm{\ρ}}_l{H}_l\frac{{\∂\mathrm{\υ}}_l}{\∂t}-(2{\mathrm{\ρ}}_l{H}_l{\mathrm{\υ}}_l+{\mathrm{\ρ}}_l{C}_{\mathrm{pl}}T+\frac{3}{2}{\mathrm{\ρ}}_l{\mathrm{\υ}}_l^2+{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ})\frac{\∂{\mathrm{\υ}}_l}{\∂z}-H_l{\mathrm{\υ}}_l\frac{\∂{\mathrm{\ρ}}_l}{\∂t}-H_l{\mathrm{\υ}}_1^2\frac{\∂{\mathrm{\ρ}}_l}{\∂z}\\ -{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ}-a(T-T_e)-\frac{{\mathrm{\τ}}_{\mathrm{lb}}{S}_{\mathrm{lb}}}{A}\end{array}\right]$

Variable (υ, ρ) parametric equation is as follows:

$\left\{\begin{array}{c}{\mathrm{\ρ}}_l\frac{\∂{\mathrm{\υ}}_l}{\∂t}-{\∂}_g\frac{\∂{\mathrm{\υ}}_g}{\∂t}+{\mathrm{\υ}}_l\frac{\∂{\mathrm{\υ}}_l}{\∂t}-{\mathrm{\υ}}_g\frac{\∂{\mathrm{\ρ}}_g}{\∂t}+{\mathrm{\υ}}_l^2\frac{\∂{\mathrm{\ρ}}_l}{\∂z}-{\mathrm{\υ}}_g^2\frac{\∂{\mathrm{\ρ}}_g}{\∂z}+2{\mathrm{\ρ}}_l{\mathrm{\υ}}_l\frac{\∂{\mathrm{\υ}}_l}{\∂z}-2{\mathrm{\ρ}}_g{\mathrm{\υ}}_g\frac{\∂{\mathrm{\υ}}_g}{\∂z}+\frac{{\mathrm{\ρ}}_{l}g\cos \mathrm{\θ}}{H_l}-\frac{{\mathrm{\ρ}}_{g}g\cos \mathrm{\θ}}{H_g}+\frac{{\mathrm{\τ}}_{\mathrm{lb}}{S}_{\mathrm{lb}}}{{\mathrm{AH}}_l}=0\\ H_l\frac{\∂{\mathrm{\ρ}}_l}{\∂t}+{\mathrm{\ρ}}_l{H}_l\frac{\∂{\mathrm{\υ}}_l}{\∂z}+{\mathrm{\υ}}_l{H}_l\frac{\∂{\mathrm{\ρ}}_l}{\∂z}=0\\ H_g\frac{\∂{\mathrm{\ρ}}_g}{\∂t}+{\mathrm{\ρ}}_g{H}_g\frac{\∂{\mathrm{\υ}}_g}{\∂z}+{\mathrm{\υ}}_g{H}_g\frac{\∂{\mathrm{\ρ}}_g}{\∂z}=0\\ \frac{\∂T_e}{\∂t_D}=\frac{{\mathrm{\λ}}_e}{C_{\mathrm{Pe}}{\mathrm{\ρ}}_e}(\frac{{\∂}^2{T}_e}{\∂r_D^2}+\frac{1}{r_D}\frac{\∂T_e}{\∂r_D}),\frac{\∂T_e}{\∂r_D}{|}_{r_D=1}=-\frac{\mathrm{dQ}}{\mathrm{dz}}{\left(2\mathrm{\π}{\mathrm{\λ}}_e\right)}^{-1},\frac{\∂T_e}{\∂r_D}{|}_{r_D\→\∞}=0\\ H_l=H_g=1,C_{\mathrm{pl}}=\frac{{\mathrm{\ρ}}_a{H}_a{C}_{\mathrm{pa}}+{\mathrm{\ρ}}_w{H}_w{C}_{\mathrm{pw}}}{{\mathrm{\ρ}}_l{H}_l}\\ P(z,0)=P_0\left(z\right),T(z,0)=T_0\left(z\right),\mathrm{\υ}(z,0)={\mathrm{\υ}}_0\left(z\right),T_e(r_D,0)=T_e^0\left(r_D\right)0,\≤z\≤Z\\ P(0,t)={\stackrel{\‾}{P}}_0\left(t\right),T(0,t)={\stackrel{\‾}{T}}_0\left(t\right),\mathrm{\υ}(0,t)={\stackrel{\‾}{\mathrm{\υ}}}_0\left(t\right),T_e(0,t)=T_e^0\left(t\right)\end{array}\right.---\left(22\right).$

3. temperature, pressure method for numerical simulation in High Temperature High Pressure oil gas inclined shaft closing well process as claimed in claim 1, is characterized in that, in step B, adopts finite differential algorithm to solve coupling model and comprises:

B1. finite differential framework is set up;

B2. design conditions parameter;

B3. based on finite differential framework, coupling model is solved.

4. temperature, pressure method for numerical simulation in High Temperature High Pressure oil gas inclined shaft closing well process as claimed in claim 3, is characterized in that, set up finite differential framework and comprise described in step B1:

By F (z _j, t _k) be designated as wherein, z _j=jh, j=1,2 ..., m, t _k=k τ, k=1,2 ..., n

At point (Z _j+1/2, t _k) place pair carry out discretize respectively with U, can obtain:

$\left\{\begin{array}{c}{\left(\frac{\∂U}{\∂t}\right)}_{j+1/2}^{k+1}=\frac{U_{j+1}^{k+1}+U_j^{k+1}-U_{j+1}^k-U_j^k}{2\mathrm{\τ}}\\ {\left(\frac{\∂U}{\∂z}\right)}_{j+1}^{k+1}=\frac{U_{j+1}^{k+1}-U_j^{k+1}}{h}\\ U_{j+1/2}^{k+1}=\frac{U_{j+1}^{k+1}+U_j^{k+1}}{2}\end{array}\right.---\left(23\right).$

5. temperature, pressure method for numerical simulation in High Temperature High Pressure oil gas inclined shaft closing well process as claimed in claim 3, it is characterized in that, in step B2, described design conditions parameter comprises:

(a) heat diffusion equation:

$f\left(t_D\right)=\left\{\begin{array}{c}1.128\sqrt{t_D}(1-0.3\sqrt{t_D})t_D\≤1.5\\ (0.4063+0.5\ln t_D)(1+0.6/t_D)t_D>1.5\end{array}\right.$

Wherein t _d=at/r ² _wh;

B () pit shaft diverse location is to the heat-conduction coefficient U of second contact surface _to:

$\frac{1}{U_{\mathrm{to}}}=r_{\mathrm{ti}}\frac{1}{{\mathrm{\λ}}_{\mathrm{ins}}}\ln \left(\frac{r_{\mathrm{ci}}}{r_{\mathrm{to}}}\right)+\frac{1}{h_c+h_r}+r_{\mathrm{ti}}\frac{1}{{\mathrm{\λ}}_{\mathrm{cem}}}\ln \left(\frac{r_{\mathrm{cem}}}{r_{\mathrm{co}}}\right);$

Friction force between (c) fluid and tube wall:

${\mathrm{\τ}}_{\mathrm{lb}}{S}_{\mathrm{lb}}/A=f{\mathrm{\ρ}}_l{\mathrm{\υ}}_l^2/\left(2d\right);$

(d) friction factor f:

$1/\sqrt{f}=1.14-21g(e/d+21.25/{\mathrm{Re}}^{0.9})$

E () gas-liquid mixed ratio of specific heat is held:

6. temperature, pressure method for numerical simulation in High Temperature High Pressure oil gas inclined shaft closing well process as claimed in claim 3, is characterized in that, in step B3, described solving coupling model based on finite differential framework comprises:

B31: the step-length of setting-up time and the degree of depth, obtains the cant angle theta of every bit _j=θ _j-1+ (θ _k-θ _k-1) Δ s _j/ Δ s _k,

Wherein j represents calculating infinitesimal, s _krepresent to tilt for θ _k, time measure the degree of depth;

B32: specify suitable finite differential grid;

B33: given starting condition and boundary condition;

B34: make T=T _k, solve the T that following system of equations obtains time t _e:

$\left\{\begin{array}{c}\frac{\∂T_e}{\∂t_D}=\frac{{\mathrm{\λ}}_e}{C_{P_e{\mathrm{\ρ}}_e}}\frac{{\∂}^2{T}_e}{{\∂r}_D^2}+\frac{{\mathrm{\λ}}_e}{r_D{C}_{P_e{\mathrm{\ρ}}_e}}\frac{{\∂T}_e}{\∂r_D}\\ {\frac{\∂T_e}{\∂r_D}|}_{r_D=1}=-\mathrm{\α}(T-T_e){\left(2\mathrm{\π}{\mathrm{\λ}}_e\right)}^{-1}\\ {\frac{\∂T_e}{{\∂r}_D}|}_{r_D\→\∞}=0\\ T_e\left(0\right)=T_e^0\end{array}\right.---\left(24\right)$

Order for the temperature of time t, degree of depth z, i=1,2 ..., t _s, j=1,2 ..., N, wherein t _stime and last radial point is represented respectively with N, it is the initial temperature on stratum; Finite differencing method is applied to formula (24) to obtain

τ and represent respectively between time interval and radial portion; Canonical form can be converted into

Differential method is used to carry out sliding-model control to boundary condition; To r _d=1, have

${\frac{\∂T_e}{\∂r_D}|}_{r_D=1}-\frac{\mathrm{\α}}{2\mathrm{\π}{\mathrm{\λ}}_e}T{|}_{r_D=1}=-\frac{\mathrm{\α}{T}_k}{2\mathrm{\π}{\mathrm{\λ}}_e}---\left(27\right)$

So

To r _d=N, has

$T_{e,N}^{i+1}-T_{e,N-1}^{i+1}=0---\left(29\right)$

Integrate (26), (28) and (29), the calculably symbolic solution of layer temperature, can obtain T _ediscretize matrix $\left[T_{e,j}^i\right];$

B35: use 3-diagonal matrix Algorithm for Solving coupling model;

B36: along time dimension, performs B34 and B35, until the time terminates.