A kind of gas reservoir state parameter is with boring numerical inversion analysis method
Technical field
The present invention relates to gas well technical field of measurement and test, more particularly to a kind of gas reservoir state parameter is analyzed with numerical inversion is bored
Method.
Background technology
Gas well test jobs are the technological means that oil/gas well exploration and development industry commonly directly acquires gas reservoir parameter.Gas
Usually there are two types of methods for well test:Gas well completion is tested and drillstem test, and well completing test is that tripping in test is closed after gas well finishing drilling
Note produces natural gas, analyzes gas flow and stress reaction, the physical property of analysis extraction gas sample, to obtain gas-bearing formation state ginseng
Number.Drillstem test is carried out in gas-bearing formation drilling process, is replaced test string with drill string, is executed test jobs, realizes gas
The description of layer state parameter.The two main distinction is the tubing string difference and timing node difference for test, but major technique is special
Point is similar.
The major defect that gas well test method obtains gas-bearing formation state parameter is shown as:
1. extending well construction period, increase shaft building cost.Gas well test jobs needs are carried out in the case where stopping brill state, and are needed
Tripping in and taking-up test string, time-consuming.Further, since drillng operation remains static when test jobs carry out, easily draw
The complex accidents such as underground bit freezing are sent out, further extend well construction period, increase shaft building cost.
2. test parameter is gas-bearing formation population mean, it is difficult to obtain gas reservoir state parameter and change section with well depth.Due to completion
Time-consuming, of high cost with drillstem test operation for test, therefore test jobs are only carried out in special layers position, and test jobs are only capable of making
The W-response of finishing drilling gas-bearing formation section state parameter, it is difficult to reflect otherness of the gas reservoir with depth.
Gas detection logging is a kind of with boring natural gas ground test technology, is obtained using the degasser before vibrating screen is placed in
The gas entrained by drilling fluid returned from shaft bottom, carries out it monitoring of component and content and edits and records, and judges ground laminar flow accordingly
Volume property is evaluated and is explained to reservoir indirectly.
Gas logging technology is that outlet is returned on ground to the foundation that evaluating reservoir is explained, natural gas is from gas reservoir to the defeated of ground
Fortune during experience dissolving, be precipitated and gas-liquid two-phase migrate complex process, the technology mainly have there are two disadvantage:
1. time-lag effect.From bore meet gas reservoir, natural gas intrusion drilling fluid detect natural gas to ground, there are it is larger when
Between hysteresis quality, be difficult to gas reservoir depth so as to cause the gas reservoir parameter that ground obtains accurate corresponding;
2. errors of analytical results is big.The amount of natural gas that ground obtains is gas reservoir output tolerance and dissolving in pit shaft and delay
The difference of amount of natural gas, and different gas, there are dissolubility difference, old place face capture tolerance and gas component are invaded with shaft bottom gas reservoir
Entering amount and component, there are larger differences, and the natural gas that direct basis ground obtains carries out component and content analysis, error are larger.
Invention content
The present invention in view of the drawbacks of the prior art, provides a kind of gas reservoir state parameter with boring numerical inversion analysis method,
It can effectively solve the problem that the above-mentioned problems of the prior art.
In order to realize the above goal of the invention, the technical solution adopted by the present invention is as follows:
A kind of gas reservoir state parameter includes the following steps with numerical inversion analysis method is bored:
1) governing equation is established, is included the following steps:
11) single liquid phase wellbore pressure field model is established;Formula is as follows:
In formula:For stagnation pressure force gradient, Pa/m;
ZvFor vertical depth, m;
For single liquid phase drilling fluid friction pressure gradient, Pa/m;
Liquid drilling fluid friction pressure gradient is calculated using Herschel-Bulkley rheological models;
12) single liquid phase pit shaft heat exchange models are established, are included the following steps:
121) drilling fluid model in drill string is established
122) drill string tube body model is established
123) annular space drilling fluid model is established
124) stratigraphic model is established
125) borehole wall and bed boundary model are established
13) biphase gas and liquid flow heat exchange models are established, are included the following steps:
131) establishing mass-conservation equation includes:
Establish drilling fluid mass-conservation equation:
Establish free gas mass-conservation equation:
Establish solution gas mass-conservation equation:
132) momentum conservation equation, such as following formula are established:
133) energy conservation equation is established, including:
Establish drill string self-energy conservation equation:
Establish annular space energy conservation equation:
Establish pit shaft/stratum interface energy conservation equation:
F (t) is the instantaneous heat-loss function of zero dimension in formula.
14) gas reservoir flow model in porous media is established, is included the following steps:
141) the non-darcy unstable flow through porous media model of low-pressure gas is established:
142) the non-darcy transient seepage flow relational model of high pressure gas is established:
15) natural gas PVT models are established, formula is as follows:
Wherein:
ρr=0.27pr/(Z·Tr), pr=p/pc,Tr=T/Tc, Z=pv/ (RT), ρ=1/v
Other parameters are as follows:
A1=0.3265, A2=-1.0700, A3=-0.5339, A4=0.01569
A5=-0.05165, A6=0.5475, A7=-0.7361, A8=0.1844
A9=0.1056, A10=0.6134, A11=0.7210;
Step 16) establishes Natural Gas Viscosity model, and formula is as follows:
R1:0 < ρg≤0.26g/cm3,0.7≤Tr≤3.4,16≤Mw≤72.2
μ=5.6563271 × 10-2ρg+4.9374602×10-3Tr+4.1949307×10-5Mw+2.93978342×10-3
R2:0.26 < ρg≤0.46g/cm3,0.88≤Tr≤3.4,16≤Mw≤72.2
μ=1.3707401 × 10-1ρg+2.72296913×10-3Tr+7.64990184×10-6Mw+4.050623771×
10-3
R3:0.46 < ρg≤0.595g/cm3,0.645≤Tr≤1.5,21.6≤Mw≤119
μ=4.8834347 × 10-1ρg+4.437225271×10-2Tr+6.5756117×10-4Mw-2.002744453×
10-1
R4:0.595 < ρg≤0.76g/cm3,0.5≤Tr≤1.4,35≤Mw≤119
The ρ of μ=1.500143567g+9.84848797×10-2Tr+1.407003797×10-3Mw-8.257944449×10-1(16)
17) rate of dissolution model is established, formula is as follows:
In formula:To consider the saturation solubility of three kinds of key components;
cbulbFor average solubility;
K is mass transport coefficient, when laminar flow:K=0.026Re0.8Sc1/3D/(dout-dinner), when turbulent:
18) biphase gas and liquid flow flow-rate profile is established, formula is as follows:
vg=C0vm+v∞ (18)
In formula:
C0--- velocity profile coefficient.The ratio of speed and mean flow rate at two-phase flow center depends on flow pattern (blister
Stream, slug flow, agitation stream, annular flow);
vm--- biphase gas and liquid flow mixture average speed;
v∞--- the opposite liquid phase sliding velocity of gas phase depends on flow pattern.
2) domain mesh generation and inserting knot are solved, using staggered-mesh technical arrangement temperature and pressure velocity node position
It sets, wherein temperature nodes are located at grid element center, and pressure and speed term are arranged at net boundary, to ensure that discrete variable is tight
Lattice meet governing equation and solve the continuous of variable.
Axial coordinate is directed toward shaft bottom from well head, and axial grid number is gradually increased from well head to shaft bottom, well head node serial number
It is 1, shaft bottom increases the dummy node N+1 that a grid axial length is 0, the grid axis of each zone of dispersion same node point number
It is identical to position, it is arranged at net boundary at area of passage change;
During oil-gas diffusion transported simulation, in new time step, gas phase migration distance does not allow more than Gridding length,
It is no that it will cause calculating, convergence problem, i.e., the step-length of new time step should not meet:
3) governing equation is discrete, includes the following steps:
31) mass-conservation equation is discrete
Discrete mass conservation equation 7, equation 8 and equation 9 obtain:
Liquid phase:
Free gas phase:
Dissolve gas phase:
tn-1Moment mesh parameter and upstream (i+1) node tnMoment parameter is it is known that event, conservation of mass discrete equation solve
Substantially tnThe mass transport item of liquid phase, solution gas and free gas is sought at moment grid downstream boundary;Cell is pressed in annular space
It is divided into three kinds of grid cells according to internal flow phase composition:
1. gas-liquid two-phase hybrid grid unit;
2. gas phase leading edge grid, a part is gas-liquid two-phase in the grid, and a part is pure liquid phase;
3. pure liquid phase grid.Grid on grid top where gas phase leading edge is pure liquid phase grid.
311) gas-liquid two-phase grid downstream boundary mass transport
The three kinds of material masses flowed out at gas-liquid two-phase net boundary, can be according to time step, a upper time step
The respective speed of gas volume fraction, the respective density of a upper time step and current time is calculated.
Liquid phase:
Free gas:
Solution gas:
312) pure liquid phase grid downstream boundary mass transport
At pure liquid phase net boundary flow out only without solution gas liquid phase, can according to time step and it is current when
Between speed be calculated.
Liquid phase:
Free gas:
Solution gas:
313) gas phase leading edge grid downstream boundary mass transport
tn-1Moment is to tnGas phase velocity is at moment downstream boundaryLiquid velocity isWithout dissolving in pure liquid phase
Gas.
tn-1The a length of L in gas-liquid two-phase part in moment node ig, pure liquid phase part length is (△ Zi-Lg), gas-liquid two-phase portion
Point gas volume fraction beDensity of liquid phase isDensity of gas phase isPure density of liquid phase is ρl。
Gas phase leading edge is passed through
It reaches at grid i downstream boundaries afterwards, pure liquid velocity is grid gas-liquid two-phase average speed in △ t moments, i.e.,
①ΔtnWhen≤Δ t, i.e.,
ΔtnPeriod gas phase leading edge is less than downstream boundary, the only pure liquid phase being discharged in this case, the quality of discharge
For:
Liquid phase:
Solution gas:
Free gas:
②Δtn>When Δ t
ΔtnPeriod gas phase leading edge is more than downstream boundary, whether contains solution gas in two kinds of situation according to the substance of discharge
Determine the quality for excluding substance.
a.
ΔtnThe liquid phase volume of period discharge is not more than tn-1Pure liquid phase volume in moment grid, is not discharged solution gas.
Liquid phase:
Solution gas:
Free gas:
b.
ΔtnThe liquid phase volume of period discharge is more than tn-1Pure liquid phase volume in moment grid, be discharged liquid phase in include
Solution gas.
Liquid phase:
Solution gas:
Free gas:
In formulaFor density of liquid phase of the i grid tn-1 moment containing solution gas, ρLFor the density of liquid phase without solution gas.
314) grid solution gas quality is sought
Δ t period solution gas quality seek region be grid upstream boundary to gas phase leading edge, rate of dissolution is according to n-1
The substance swum over on moment grid in grid gas phase front edge area is constituted and flow performance is sought, and is ignored in the period
The solution gas that the gas of stratum intrusion generates.
Solution gas quality after the Δ t times in static grid is:
The equation is scalar nonlinear equation, is iteratively solved outAfterwards, then grid solution gas matter in the Δ t periods is obtained
Measuring incrementss is:
32) momentum conservation equation is discrete
Annular space drilling fluid flows up, and axial coordinate direction increases from top to bottom, and the two direction is on the contrary, by speed when discrete
Item takes positive value, obtains the momentum conservation equation of discrete scheme:
Containing in the up-front grid of gas phase, friction pressure drop is sought respectively to single liquid phase region and Gas-liquid phase region, by adding
Power sum the grid total pressure drop.
33) energy conservation equation is discrete
Using specific enthalpy difference quotient thermodynamic argument, pit shaft/bed boundary temperature is eliminated, the annular space section indicated with thermal resistance is obtained
Point temperature computation discrete equation:
In formula,For the entire thermal resistance of annular space drilling fluid to stratum.
4) gas reservoir parameter values inversion algorithm, includes the following steps:
41) grid division
According to casing programme and drilling tool structure grid division, setting node array pNode [N] and pNodeNew [N].
Upper time dot grid data are stored in pNode, pNodeNew stores new time point data.
42) the vertical coordinate of grid is calculated, pNode node arrays are updated;
The vertical coordinate and grid vertical length of grid element center are calculated according to well track.
43) it is flowed according to single liquid phase and carries out temperature field analysis, update pNode node arrays
44) grid flow parameter and calculation of pressure are carried out simultaneously according to single liquid phase.Including liquid phase quality, outflow quality, outflow
Mass flow, density of liquid phase, liquid phase flow rate, barometric gradient, pressure, void fraction assign 0, gas phase flow velocity and assign the outflow of 0, gas phase
Quality assigns 0, solution gas quality and assigns 0, outflow solution gas quality tax 0, updates pNode node arrays.
45) gas cut node serial number M is determined
Grid is traversed, the gentle layer depth of grid node axial position is compared, determines gas cut node serial number M.
46) diffusion And Movement simulates primary condition application
Diffusion And Movement simulated time t=0, gas phase forward position grid indexGas phase length in the grid of gas phase forward position(processing of point gas cut).
Grid variable pNodeNew [N] is synchronous for pNode [N] numerical value.
47) time step dt is determined according to Gridding length and grid gas phase velocity;
48) time stepping method, t=t+dt,
49) gas cut upstream node pressure assigns initial value
Gas cut upstream node pressure adjusted value pDiff=0,
410) gas cut upstream node pressure is determined,
411) according to gas reservoir factor, drilling parameter andCalculate gas cut amount
412) mesh parameter is calculated according to trellis-type
Grid is traversed from bottom to top, and trellis-type is determined according to grid index i, calculates grid flow parameter.
It is flowed for single liquid phase in grid
Gas phase forward position grid, updateWith
Gas-liquid two-phase stream trellis
413) compare well head node pressure and atmospheric pressure, calculate relative error εp, determine bottom boundary pressure adjusted value
εp=| pDiff |/101000
If 414) εp>10-3, step 410) is gone to, is repeated, until εp<10-3, new time point mesh parameter, which calculates, to be tied
Beam;
415) according to new moment grid flow parameter and the new temperature field of Pressure Analysis;
416) it is pNodeNew to update whole annular space grid data pNode;
If 417) time reaches analysis time, analysis terminates, otherwise, goes to step 47) and repeat.
Further, using gas reservoir state parameter mould is analyzed with numerical inversion means of interpretation establishment gas reservoir state parameter is bored
Quasi- device, including function module are as follows:
ReadKickData
ReadKickData is the data input function of analysis mode device, completes the input of simulation all data of well, mainly
Including casing programme, drilling tool structure, well track, stratum vertical temperature distribution, reservoir parameter:Temperature, porosity, is oozed at pressure
Saturating rate, depth, drainage radius, former oil and gas thermal physical property parameter, the thermal physical property parameter of cement, stratum, drill string etc., drilling
Parameter:Discharge capacity, rate of penetration, torque, bit nozzle, analog parameter:Whether axial mesh scale, simulation cycle consider to be crushed
Whether gas considers heat source;
KickFlowBehavior
KickFlowBehavior functions are the general function modules of simulator, and other functions are assembled according to program flow diagram,
Complete drilling well oil-gas diffusion and transported simulation analysis and data storage function;
TemperatureAnalysis
TemperatureAnalysis functions carry out bored shaft temperature field analysis according to grid flow parameter, by temperature
In result of calculation update to grid node array pNodeNew;
KickSimulating
KickSimulating functions thoroughly do away with grid node temperature and shaft bottom and carry out drilling well oil gas expansion with well head boundary condition
It dissipates and transported simulation is analyzed, and will be in the update to grid node array pNodeNew of flow parameter result of calculation;
GridGeneration
GridGeneration functions solve domain to temperature and pressure field according to casing programme and drill column structure and carry out axially
The geometry and medium information of segmentation and mesh generation, the axial geological information of save mesh node and radial heat exchange object;
NormalTPNode
NormalTPNode functions are all the grid flow parameter computational methods of biphase gas and liquid flow according to entrance and exit,
Analyze grid flow parameter;
GasFrontNode
GasFrontNode functions analyze grid flow parameter, and update gas phase forward position according to gas phase forward position trellis algorithm
Grid index and gas phase zone length data;
VerticalCoordinate
The vertical coordinate and grid vertical length of grid element center are calculated according to well track.Function performance requires to call
DirectionParaCal functions calculate vertical depth according to well depth and hole angle;
TOriginGeneration
TOriginGeneration functions generate node vertical depth according to the vertical temperature distribution data interpolation on stratum
The original temperature at place;
Ini
Ini function pairs solve grid node variable application primary condition in domain, assign initial value;
HeatResistanceSA
Calculate the heat exchanged thermoresistance for boring drilling fluid and annular space drilling fluid in main.Realize the function performance need call HPipe and
HAnnu functions calculate liquid phase and are boring main interior and forced-convection heat transfer coefficient of the biphase gas and liquid flow in annular space;
HeatResistanceAE
Calculate the heat exchanged thermoresistance between annular space biphase gas and liquid flow drilling fluid and stratum.Realize that the function performance needs to call
HAnnu functions calculate forced-convection heat transfer coefficient and instantaneous stratum heat-loss function of the biphase gas and liquid flow in annular space and calculate well
Cylinder arrives the heat exchanged thermoresistance on stratum;
HPipe
The function of HPipe functions is to calculate the pipe stream forced-convection heat transfer coefficient for boring drilling fluid and drill string inner wall in main;
HAnnu
The function of HAnnu functions is to calculate the Annular cutting forced-convection heat transfer of annular space biphase gas and liquid flow drilling fluid and the borehole wall
Coefficient;
PPTP
According to total aggregation function of macroscopic property of the media type calculation medium under given temperature and pressure, need
According to specific media type crude oil, water, natural gas;The hot Calculation of Physical Properties function of respective media is called to complete its function;
DPDZ_TwoPhase
Biphase gas and liquid flow flow parameter and pressure drop gradient, biphase gas and liquid flow stream are calculated according to biphase gas and liquid flow mechanism model
Dynamic parameter includes:Gas phase velocity, liquid velocity, void fraction, density of gas phase, flow pattern;
TInDrillstem
TInDrillstem functions are iterated to calculate according to single liquid phase energy conservation equation bores drilling fluid node temperature in main;
TAnnulus
TAnnulus functions iterate to calculate annular space drilling fluid node temperature according to biphase gas and liquid flow energy conservation equation;
CpGas calculates natural gas specific heat at constant pressure;
CpOil calculates crude oil specific heat at constant pressure;
CpWater calculates water flooding specific heat at constant pressure;
DenGas calculates natural gas density;
DenOilSaturated calculates saturation oil density;
DenOilUnsaturated calculates unsaturated oil density;
RsStanding calculates crude oil dissolved gas oil ratio;
SurfaceTension calculates oil water interfacial tension;
ThconGas calculates natural gas thermal conductivity;
ThconOil calculates crude oil thermal conductivity;
ThconWater calculates water flooding thermal conductivity;
ViscGas calculates Natural Gas Viscosity;
ViscOil calculates viscosity of crude;
ViscWater calculates water flooding viscosity.
Further, gas reservoir state parameter is arranged with the boundary condition for boring numerical inversion means of interpretation with primary condition:
Boundary condition
1. drill string inlet:
Qm=C (45-1)
T1(z=0)=TE (45-2)
Wherein TEFor drilling fluid inlet temperature, QmFor drilling fluid displacement.
2. annular space exit:
P=patm (46)
3. shaft bottom drill string and annular space junction, drilling fluid temperature is identical, has
Primary condition
1. the primary condition in temperature field is the temperature in wellbore calculated under single liquid-phase condition,
I.e.:
Ti(z, 0)=Ti,st(z) (48)
2. the primary condition of pressure field is the wellbore pressure calculated under single liquid-phase condition,
I.e.:
pi(z, 0)=pi,st(z) (49)
3. gas cut starts, without gas phase in each grid of annular space:
αi(z, 0)=0 (50).
Compared with prior art the advantage of the invention is that:
1. the gas reservoir state parameter in the case of drillng operation is uninterrupted improves gas well exploitation effect with accurate analysis in real time is bored
Benefit.Be integrated with 18 governing equations, considered gas liquid two-phase flow in gas reservoir to the seepage flow, pit shaft of pit shaft and heat exchange,
It is anti-can to survey parameter using this method according to ground gas for dissolved, temperature and pressure influence to gas phase hot physical property of the gas phase in drilling fluid
Drill to obtain gas reservoir state parameter along well depth section, gas reservoir state parameter is finer.
2. well kick risk identifies in time, physical property and flow regime of the parameter acquiring gas phase along pit shaft can be surveyed according to ground gas
Parameter distribution, and obtain shaft bottom gas cut amount, to accurately identify well kick degree and well kick risk, to take precautions against blowout risk and
Ensure that drillng operation is of great significance safely.
Description of the drawings
Fig. 1 is mesh generation schematic diagram of the embodiment of the present invention;
Fig. 2 is node i of the embodiment of the present invention and i+1 mesh parameter position views;
Fig. 3 is gas phase of embodiment of the present invention forward position grid mass transportation schematic diagram;
Fig. 4 is GasReservoirAnalysis emulation module structure charts of the embodiment of the present invention.
Specific implementation mode
To make the objectives, technical solutions, and advantages of the present invention more comprehensible, develop simultaneously embodiment referring to the drawings,
The present invention is described in further details.
A kind of gas reservoir state parameter is with boring numerical inversion analysis method, including with lower part:
One, governing equations
1) single liquid phase wellbore pressure field model
In formula:For stagnation pressure force gradient, Pa/m;
ZvFor vertical depth, m;
For single liquid phase drilling fluid friction pressure gradient, Pa/m.
Liquid drilling fluid friction pressure gradient is calculated using Herschel-Bulkley rheological models.
2) single liquid phase pit shaft heat exchange models
1. drilling fluid in drill string
2. drill string tube body
3. annular space drilling fluid
4. stratum
5. the borehole wall and bed boundary
3) biphase gas and liquid flow heat exchange models
1. mass-conservation equation
Drilling fluid mass-conservation equation:
Free gas mass-conservation equation:
Solution gas mass-conservation equation:
2. momentum conservation equation
3. energy conservation equation
Drill string self-energy conservation equation:
Annular space energy conservation equation:
Pit shaft/stratum interface energy conservation equation:
F (t) is the instantaneous heat-loss function of zero dimension in formula.
4) gas reservoir flow model in porous media
The non-darcy unstable flow through porous media model of low-pressure gas:
The non-darcy transient seepage flow relational model of high pressure gas:
5) natural gas PVT models:
Wherein:
ρr=0.27pr/(Z·Tr), pr=p/pc,Tr=T/Tc, Z=pv/ (RT), ρ=1/v
Other parameters are as follows:
A1=0.3265, A2=-1.0700, A3=-0.5339, A4=0.01569
A5=-0.05165, A6=0.5475, A7=-0.7361, A8=0.1844
A9=0.1056, A10=0.6134, A11=0.7210
6) Natural Gas Viscosity model:
R1:0 < ρg≤0.26g/cm3,0.7≤Tr≤3.4,16≤Mw≤72.2
μ=5.6563271 × 10-2ρg+4.9374602×10-3Tr+4.1949307×10-5Mw+2.93978342×10-3
R2:0.26 < ρg≤0.46g/cm3,0.88≤Tr≤3.4,16≤Mw≤72.2
μ=1.3707401 × 10-1ρg+2.72296913×10-3Tr+7.64990184×10-6Mw+4.050623771×
10-3
R3:0.46 < ρg≤0.595g/cm3,0.645≤Tr≤1.5,21.6≤Mw≤119
μ=4.8834347 × 10-1ρg+4.437225271×10-2Tr+6.5756117×10-4Mw-2.002744453×
10-1
R4:0.595 < ρg≤0.76g/cm3,0.5≤Tr≤1.4,35≤Mw≤119
The ρ of μ=1.500143567g+9.84848797×10-2Tr+1.407003797×10-3Mw-8.257944449×10-1
(16)
7) rate of dissolution model
In formula:To consider the saturation solubility of three kinds of key components;
cbulbFor average solubility;
K is mass transport coefficient, when laminar flow:K=0.026Re0.8Sc1/3D/(dout-dinner), when turbulent:
8) biphase gas and liquid flow flow-rate profile
vg=C0vm+v∞ (18)
In formula:
C0--- velocity profile coefficient.The ratio of speed and mean flow rate at two-phase flow center depends on flow pattern (blister
Stream, slug flow, agitation stream, annular flow);
vm--- biphase gas and liquid flow mixture average speed;
v∞--- the opposite liquid phase sliding velocity of gas phase depends on flow pattern.
Two, domain mesh generation and inserting knot are solved
Present invention application Finite Volume Method is discrete to mass-conservation equation and energy conservation equation progress, using staggeredly
Grid arranges that temperature and pressure velocity node position, wherein temperature nodes are located at grid element center, pressure and speed term arrangement
At net boundary, to ensure discrete variable strictly meet governing equation and solve variable it is continuous.
Axial coordinate is directed toward shaft bottom from well head, and axial grid number is gradually increased from well head to shaft bottom, well head node serial number
It is 1, shaft bottom increases the dummy node N+1 that a grid axial length is 0, the grid axis of each zone of dispersion same node point number
It is identical to position, it is arranged at net boundary at area of passage change.Mesh generation is as shown in Figure 1.
During oil-gas diffusion transported simulation, in new time step, gas phase migration distance does not allow more than Gridding length,
It is no that it will cause calculating, convergence problem, i.e., the step-length of new time step should not meet:
Three, governing equation is discrete
By taking annular space node i as an example, time and spatial spreading are carried out to governing equation.The physical quantity that governing equation is related to
For:p、T、ρm、ρg、vm、vg、α、f1.In annular space, the bottom-up flowing of drilling fluid, node i parameter is by node i+1
Parameter influences, and annular space node i and+1 mesh parameter of node i arrangement are as shown in Figure 2.Pressure can be considered linear distribution within a grid.
It is discrete that governing equation is carried out accordingly.
1) mass-conservation equation is discrete
Discrete mass conservation equation (7), (8) and (9) obtains:
Liquid phase:
Free gas phase:
Dissolve gas phase:
tn-1Moment mesh parameter and upstream (i+1) node tnMoment parameter is it is known that event, conservation of mass discrete equation solve
Substantially tnThe mass transport item of liquid phase, solution gas and free gas is sought at moment grid downstream boundary.Cell is pressed in annular space
It is divided into three kinds of grid cells according to internal flow phase composition:
1. gas-liquid two-phase hybrid grid unit;
2. gas phase leading edge grid, a part is gas-liquid two-phase in the grid, and a part is pure liquid phase;
3. pure liquid phase grid.Grid on grid top where gas phase leading edge is pure liquid phase grid.
(1) gas-liquid two-phase grid downstream boundary mass transport
The three kinds of material masses flowed out at gas-liquid two-phase net boundary, can be according to time step, a upper time step
The respective speed of gas volume fraction, the respective density of a upper time step and current time is calculated.
Liquid phase:
Free gas:
Solution gas:
(2) pure liquid phase grid downstream boundary mass transport
At pure liquid phase net boundary flow out only without solution gas liquid phase, can according to time step and it is current when
Between speed be calculated.
Liquid phase:
Free gas:
Solution gas:
(3) gas phase leading edge grid downstream boundary mass transport
Gas phase forward position grid liquid phase, solution gas and free gas mass transport schematic diagram are as shown in Figure 3.
tn-1Moment is to tnGas phase velocity is at moment downstream boundaryLiquid velocity isWithout dissolving in pure liquid phase
Gas.
tn-1The a length of L in gas-liquid two-phase part in moment node ig, pure liquid phase part length is (△ Zi-Lg), gas-liquid two-phase portion
Point gas volume fraction beDensity of liquid phase isDensity of gas phase isPure density of liquid phase is ρl。
Gas phase leading edge is passed through
It reaches at grid i downstream boundaries afterwards, pure liquid velocity is grid gas-liquid two-phase average speed in △ t moments, i.e.,
①ΔtnWhen≤Δ t, i.e.,
ΔtnPeriod gas phase leading edge is less than downstream boundary, the only pure liquid phase being discharged in this case, the quality of discharge
For:
Liquid phase:
Solution gas:
Free gas:
②Δtn>When Δ t
ΔtnPeriod gas phase leading edge is more than downstream boundary, whether contains solution gas in two kinds of situation according to the substance of discharge
Determine the quality for excluding substance.
a.
ΔtnThe liquid phase volume of period discharge is not more than tn-1Pure liquid phase volume in moment grid, is not discharged solution gas.
Liquid phase:
Solution gas:
Free gas:
b.
ΔtnThe liquid phase volume of period discharge is more than tn-1Pure liquid phase volume in moment grid, be discharged liquid phase in include
Solution gas.
Liquid phase:
Solution gas:
Free gas:
In formulaFor density of liquid phase of the i grid tn-1 moment containing solution gas, ρLFor the density of liquid phase without solution gas.
(4) grid solution gas quality is sought
Δ t period solution gas quality seek region be grid upstream boundary to gas phase leading edge, rate of dissolution is according to n-1
The substance swum over on moment grid in grid gas phase front edge area is constituted and flow performance is sought, and is ignored in the period
The solution gas that the gas of stratum intrusion generates.
Solution gas quality after the Δ t times in static grid is:
The equation is scalar nonlinear equation, is iteratively solved outAfterwards, then grid solution gas matter in the Δ t periods is obtained
Measuring incrementss is:
2) momentum conservation equation is discrete
Annular space drilling fluid flows up, and axial coordinate direction increases from top to bottom, and the two direction is on the contrary, by speed when discrete
Item takes positive value, obtains the momentum conservation equation of discrete scheme:
Containing in the up-front grid of gas phase, friction pressure drop is sought respectively to single liquid phase region and Gas-liquid phase region, by adding
Power sum the grid total pressure drop.
3) energy conservation equation is discrete
Using specific enthalpy difference quotient thermodynamic argument, pit shaft/bed boundary temperature is eliminated, the annular space section indicated with thermal resistance is obtained
Point temperature computation discrete equation:
In formula,For the entire thermal resistance of annular space drilling fluid to stratum.
Four, gas phase leading edge mesh tracing
Following variables will be used for gas phase leading edge mesh tracing task.
Ig--- include the up-front grid of gas;
Lg--- the up-front position of gas, relative to grid IgUpstream edge.
It is updated according to the following formula by mobile forward position and velocity of downstream, these variables:
IfThen
For comprising the up-front grid of gas, it is assumed that free gas and dissolved gas are all generally evenly distributed in gas
After leading edge, otherwise assume that gas is evenly distributed in entire grid.When gas leading edge passes through one in a time step
When net boundary, the speed of drilling fluid will change, since the speed of the mixture in a time step is considered normal
Number, so the speed of drilling fluid will change.Time step is division in this way, has v before gas phase leading edge passes throughm=vM,
And there is v laterm=(vM-αvg)/(1-α)。
Five, boundary condition and primary condition
1) boundary condition
1. drill string inlet:
Qm=C (45-1)
T1(z=0)=TE (45-2)
Wherein TEFor drilling fluid inlet temperature, QmFor drilling fluid displacement.
2. annular space exit:
P=patm (46)
3. shaft bottom drill string and annular space junction, drilling fluid temperature is identical, has
2) primary condition
(1) primary condition in temperature field is the temperature in wellbore calculated under single liquid-phase condition, i.e.,:
Ti(z, 0)=Ti,st(z) (48)
(2) primary condition of pressure field is the wellbore pressure calculated under single liquid-phase condition, i.e.,:
pi(z, 0)=pi,st(z) (49)
(3) when gas cut starts, without gas phase in each grid of annular space:
αi(z, 0)=0 (50)
Six, algorithm
When gas reservoir parameter values back analysis, it is first assumed that gas reservoir parameter calculates ground monitoring ginseng according to following algorithms
Number adjusts gas reservoir parameter if calculated value and monitor value error are larger, repeats numerical simulation, until ground calculated value and monitor value
Meet the condition of convergence, both true gas reservoir state parameter.It is assumed that the numerical simulation algorithm under gas reservoir parameter is as follows.
Since annular space biphase gas and liquid flow temperature, pressure field, gas reservoir temperature and pressure data etc. interdepend, therefore, it needs from initial value
It sets out, is promoted according to time step, until reaching analysis time.Each time step is needed from upper time point value, warp
Iterative calculation is crossed, until obtaining each convergence solution.
(1) grid division
According to casing programme and drilling tool structure grid division, setting node array pNode [N] and pNodeNew [N].
Upper time dot grid data are stored in pNode, pNodeNew stores new time point data.
(2) the vertical coordinate of grid is calculated, pNode node arrays are updated;
The vertical coordinate and grid vertical length of grid element center are calculated according to well track.
(3) it is flowed according to single liquid phase and carries out temperature field analysis, update pNode node arrays
(4) grid flow parameter and calculation of pressure are carried out simultaneously according to single liquid phase.Including liquid phase quality, outflow quality, outflow
Mass flow, density of liquid phase, liquid phase flow rate, barometric gradient, pressure, void fraction assign 0, gas phase flow velocity and assign the outflow of 0, gas phase
Quality assigns 0, solution gas quality and assigns 0, outflow solution gas quality tax 0, updates pNode node arrays.
(5) gas cut node serial number M is determined
Grid is traversed, the gentle layer depth of grid node axial position is compared, determines gas cut node serial number M.
(6) diffusion And Movement simulates primary condition application
Diffusion And Movement simulated time t=0, gas phase forward position grid indexGas phase length in the grid of gas phase forward position(processing of point gas cut).
Grid variable pNodeNew [N] is synchronous for pNode [N] numerical value.
(7) time step dt is determined according to Gridding length and grid gas phase velocity;
(8) time stepping method, t=t+dt,
(9) gas cut upstream node pressure assigns initial value
Gas cut upstream node pressure adjusted value pDiff=0,
(10) gas cut upstream node pressure is determined,
(11) according to gas reservoir factor, drilling parameter andCalculate gas cut amount
(12) mesh parameter is calculated according to trellis-type
Grid is traversed from bottom to top, and trellis-type is determined according to grid index i, calculates grid flow parameter.
It is flowed for single liquid phase in grid
Gas phase forward position grid, updateWith
Gas-liquid two-phase stream trellis
(13) compare well head node pressure and atmospheric pressure, calculate relative error εp, determine bottom boundary pressure adjusted value
εp=| pDiff |/101000
(14) if εp>10-3, step (10) is gone to, is repeated, until εp<10-3, new time point mesh parameter, which calculates, to be tied
Beam;
(15) according to new moment grid flow parameter and the new temperature field of Pressure Analysis;
(16) it is pNodeNew to update whole annular space grid data pNode;
(17) if the time reaches analysis time, analysis terminates, otherwise, goes to step (7) and repeat.
Seven, program module
Gas reservoir state parameter analysis mode device GasReservoirAnalysis is worked out using this method, by 63 functions
It constitutes, main module structure is as shown in Figure 4.Its main functional modules (function) is as follows.
(1)ReadKickData
ReadKickData is the data input function of analysis mode device, completes the input of simulation all data of well, mainly
Including casing programme, drilling tool structure, well track, stratum vertical temperature distribution, reservoir parameter:Temperature, porosity, is oozed at pressure
Saturating rate, depth, drainage radius, former oil and gas thermal physical property parameter, the thermal physical property parameter of cement, stratum, drill string etc., drilling
Parameter:Discharge capacity, rate of penetration, torque, bit nozzle, analog parameter:Whether axial mesh scale, simulation cycle consider to be crushed
Whether gas considers heat source.
(2)KickFlowBehavior
KickFlowBehavior functions are the general function modules of simulator, and other functions are assembled according to program flow diagram,
Complete drilling well oil-gas diffusion and transported simulation analysis and data storage function.
(3)TemperatureAnalysis
TemperatureAnalysis functions carry out bored shaft temperature field analysis according to grid flow parameter, by temperature
In result of calculation update to grid node array pNodeNew.
(4)KickSimulating
KickSimulating functions thoroughly do away with grid node temperature and shaft bottom and carry out drilling well oil gas expansion with well head boundary condition
It dissipates and transported simulation is analyzed, and will be in the update to grid node array pNodeNew of flow parameter result of calculation.
(5)GridGeneration
GridGeneration functions solve domain to temperature and pressure field according to casing programme and drill column structure and carry out axially
The geometry and medium information (grid of segmentation and mesh generation, the axial geological information of save mesh node and radial heat exchange object
Medium in control volume has:6 kinds of media such as single liquid phase drilling fluid, gas-liquid two-phase drilling fluid, steel, cement, stratum, packer fluids
Type, each medium have corresponding thermodynamics physical property).
(6)NormalTPNode
NormalTPNode functions are all the grid flow parameter computational methods of biphase gas and liquid flow according to entrance and exit,
Analyze grid flow parameter.
(7)GasFrontNode
GasFrontNode functions analyze grid flow parameter, and update gas phase forward position according to gas phase forward position trellis algorithm
Grid index and gas phase zone length data.
(8)VerticalCoordinate
The vertical coordinate and grid vertical length of grid element center are calculated according to well track.Function performance requires to call
DirectionParaCal functions calculate vertical depth according to well depth and hole angle.
(9)TOriginGeneration
TOriginGeneration functions generate node vertical depth according to the vertical temperature distribution data interpolation on stratum
The original temperature at place.
(10)Ini
Ini function pairs solve grid node variable application primary condition in domain, assign initial value.
(11)HeatResistanceSA
Calculate the heat exchanged thermoresistance for boring drilling fluid and annular space drilling fluid in main.Realize the function performance need call HPipe and
HAnnu functions calculate liquid phase and are boring main interior and forced-convection heat transfer coefficient of the biphase gas and liquid flow in annular space.
(12)HeatResistanceAE
Calculate the heat exchanged thermoresistance between annular space biphase gas and liquid flow drilling fluid and stratum.Realize that the function performance needs to call
HAnnu functions calculate forced-convection heat transfer coefficient and instantaneous stratum heat-loss function of the biphase gas and liquid flow in annular space and calculate well
Cylinder arrives the heat exchanged thermoresistance on stratum.
(13)HPipe
The function of HPipe functions is to calculate the pipe stream forced-convection heat transfer coefficient for boring drilling fluid and drill string inner wall in main.
(14)HAnnu
The function of HAnnu functions is to calculate the Annular cutting forced-convection heat transfer of annular space biphase gas and liquid flow drilling fluid and the borehole wall
Coefficient.
(15)PPTP
According to total aggregation function of macroscopic property of the media type calculation medium under given temperature and pressure, need
The hot Calculation of Physical Properties function of respective media is called to complete its function according to specific media type (crude oil, water, natural gas).
(16)DPDZ_TwoPhase
Biphase gas and liquid flow flow parameter (gas phase velocity, liquid velocity, section are calculated according to biphase gas and liquid flow mechanism model
Void fraction, density of gas phase, flow pattern) and pressure drop gradient.
(17)TInDrillstem
TInDrillstem functions are iterated to calculate according to single liquid phase energy conservation equation bores drilling fluid node temperature in main.
(18)TAnnulus
TAnnulus functions iterate to calculate annular space drilling fluid node temperature according to biphase gas and liquid flow energy conservation equation.
(19)CpGas.Calculate natural gas specific heat at constant pressure.
(20)CpOil.Calculate crude oil specific heat at constant pressure.
(21)CpWater.Calculate water flooding specific heat at constant pressure.
(22)DenGas.Calculate natural gas density.
(23)DenOilSaturated.Calculate saturation oil density.
(24)DenOilUnsaturated.Calculate unsaturated oil density.
(25)RsStanding.Calculate crude oil dissolved gas oil ratio.
(26)SurfaceTension.Calculate oil water interfacial tension.
(27)ThconGas.Calculate natural gas thermal conductivity.
(28)ThconOil.Calculate crude oil thermal conductivity.
(29)ThconWater.Calculate water flooding thermal conductivity.
(30)ViscGas.Calculate Natural Gas Viscosity.
(31)ViscOil.Calculate viscosity of crude.
(32)ViscWater.Calculate water flooding viscosity.
Those of ordinary skill in the art will understand that the embodiments described herein, which is to help reader, understands this
The implementation of invention, it should be understood that protection scope of the present invention is not limited to such specific embodiments and embodiments.This
The those of ordinary skill in field can make according to the technical disclosures disclosed by the invention various does not depart from of the invention essence
Various other specific variations and combinations, these variations and combinations are still within the scope of the present invention.