CN108509703B - Gas reservoir state parameter while-drilling numerical inversion analysis method - Google Patents
Gas reservoir state parameter while-drilling numerical inversion analysis method Download PDFInfo
- Publication number
- CN108509703B CN108509703B CN201810238082.1A CN201810238082A CN108509703B CN 108509703 B CN108509703 B CN 108509703B CN 201810238082 A CN201810238082 A CN 201810238082A CN 108509703 B CN108509703 B CN 108509703B
- Authority
- CN
- China
- Prior art keywords
- gas
- grid
- phase
- flow
- liquid phase
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000005553 drilling Methods 0.000 title claims abstract description 101
- 238000004458 analytical method Methods 0.000 title claims abstract description 36
- 239000007788 liquid Substances 0.000 claims abstract description 64
- 239000012530 fluid Substances 0.000 claims abstract description 63
- 230000005514 two-phase flow Effects 0.000 claims abstract description 45
- 238000000034 method Methods 0.000 claims abstract description 29
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 9
- 238000004090 dissolution Methods 0.000 claims abstract description 9
- 239000007789 gas Substances 0.000 claims description 320
- 239000007791 liquid phase Substances 0.000 claims description 139
- 239000012071 phase Substances 0.000 claims description 120
- 230000006870 function Effects 0.000 claims description 83
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 claims description 62
- 239000003345 natural gas Substances 0.000 claims description 31
- 238000004364 calculation method Methods 0.000 claims description 26
- 239000010779 crude oil Substances 0.000 claims description 23
- 238000004088 simulation Methods 0.000 claims description 22
- 238000004134 energy conservation Methods 0.000 claims description 20
- 238000013508 migration Methods 0.000 claims description 19
- 230000005012 migration Effects 0.000 claims description 19
- 238000012546 transfer Methods 0.000 claims description 19
- 238000011144 upstream manufacturing Methods 0.000 claims description 19
- 230000015572 biosynthetic process Effects 0.000 claims description 18
- 230000009545 invasion Effects 0.000 claims description 17
- 239000003921 oil Substances 0.000 claims description 16
- 238000009792 diffusion process Methods 0.000 claims description 15
- 239000008398 formation water Substances 0.000 claims description 9
- 239000000126 substance Substances 0.000 claims description 9
- 239000000203 mixture Substances 0.000 claims description 7
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 7
- 238000005516 engineering process Methods 0.000 claims description 6
- 230000008569 process Effects 0.000 claims description 6
- 229920006395 saturated elastomer Polymers 0.000 claims description 6
- 239000004568 cement Substances 0.000 claims description 4
- 230000008859 change Effects 0.000 claims description 4
- 239000006185 dispersion Substances 0.000 claims description 4
- 238000000638 solvent extraction Methods 0.000 claims description 4
- CKRLIWFOVCLXTP-UHFFFAOYSA-N 4-phenyl-1-propyl-3,6-dihydro-2h-pyridine Chemical compound C1N(CCC)CCC(C=2C=CC=CC=2)=C1 CKRLIWFOVCLXTP-UHFFFAOYSA-N 0.000 claims description 3
- 238000003491 array Methods 0.000 claims description 3
- 230000004323 axial length Effects 0.000 claims description 3
- 238000013500 data storage Methods 0.000 claims description 3
- 230000008571 general function Effects 0.000 claims description 3
- 238000012905 input function Methods 0.000 claims description 3
- 239000000463 material Substances 0.000 claims description 3
- 230000007246 mechanism Effects 0.000 claims description 3
- 230000035515 penetration Effects 0.000 claims description 3
- 230000035699 permeability Effects 0.000 claims description 3
- 230000011218 segmentation Effects 0.000 claims description 3
- 230000005540 biological transmission Effects 0.000 claims description 2
- 230000010354 integration Effects 0.000 claims description 2
- 230000008901 benefit Effects 0.000 abstract description 5
- 238000011161 development Methods 0.000 abstract description 3
- 238000005259 measurement Methods 0.000 abstract description 3
- 238000012360 testing method Methods 0.000 description 21
- 238000010276 construction Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 4
- 238000012544 monitoring process Methods 0.000 description 4
- 230000000704 physical effect Effects 0.000 description 3
- 238000006243 chemical reaction Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000018109 developmental process Effects 0.000 description 2
- 230000002035 prolonged effect Effects 0.000 description 2
- 229910000831 Steel Inorganic materials 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000001556 precipitation Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 239000010959 steel Substances 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention discloses a gas reservoir state parameter while-drilling numerical inversion analysis method, which comprises the following steps: establishing a control equation; solving the domain grid division and node arrangement; discretizing a control equation; gas reservoir parameter value inversion algorithm; compiling a gas reservoir state parameter analysis simulator; the boundary condition is set with the initial condition. The invention has the advantages that: the gas reservoir state parameters are accurately analyzed in real time while drilling, and the gas well development benefits are improved. A plurality of control equations are integrated, the influence of seepage of a gas reservoir to a shaft, gas-liquid two-phase flow and heat exchange in the shaft, dissolution of gas phase in drilling fluid and temperature and pressure on gas phase thermophysical properties is comprehensively considered, a gas reservoir state parameter along a well depth profile can be obtained through inversion according to a ground gas measurement parameter, and the gas reservoir state parameter is finer; the well kick degree and the well kick risk can be accurately identified, and the method has important significance for preventing the well blowout risk and ensuring the safety of the drilling operation.
Description
Technical Field
The invention relates to the technical field of gas well testing, in particular to a gas reservoir state parameter while-drilling numerical inversion analysis method.
Background
The gas well testing operation is a common technical means for directly acquiring gas reservoir parameters in the oil and gas well exploration and development industry. Gas well testing typically has two methods: the method comprises a gas well completion test and a midway test, wherein the completion test is to pay test attention after the gas well is drilled, to extract natural gas, to analyze gas flow and pressure reaction, to analyze the physical properties of an extracted gas sample, and to obtain gas layer state parameters. And the midway test is carried out in the gas layer drilling process, a drill column is used for replacing a test pipe column, the test operation is executed, and the description of the gas layer state parameters is realized. The main difference between the two methods is that the pipe columns used for testing are different and the time nodes are different, but the main technical characteristics are similar.
The main disadvantages of the gas well testing method for obtaining the gas layer state parameters are as follows:
1. the well construction period is prolonged, and the well construction cost is increased. The gas well test operation needs to be carried out in a drill stopping state, and a test pipe column needs to be put in and taken out, so that the time consumption is long. In addition, as the drilling operation is in a static state when the test operation is carried out, complex accidents such as underground drill sticking and the like are easily caused, the well construction period is further prolonged, and the well construction cost is increased.
2. The test parameters are the overall mean value of the gas reservoir, and the profile of the gas reservoir state parameters changing along with the well depth is difficult to obtain. As the well completion test and the midway test operation take long time and have high cost, the test operation is only carried out at a special layer, the test operation can only make the integral reaction of the state parameters of the drilled gas layer section, and the difference of the gas reservoir along with the depth is difficult to reflect.
The gas logging is a ground test technology of natural gas while drilling, and gas carried by drilling fluid returned from the bottom of a well is obtained by utilizing a degasser arranged in front of a vibrating screen, and the gas is subjected to component and content monitoring and cataloguing, so that the property of formation fluid is judged according to the component and the content, and a reservoir stratum is indirectly evaluated and explained.
The basis of reservoir evaluation and interpretation by the gas logging technology is ground returned gas, and the natural gas undergoes complex processes of dissolution, precipitation and gas-liquid two-phase migration in the process of transporting from a gas reservoir to the ground, and the technology mainly has two defects:
1. time lag property. The natural gas is detected from the drilling fluid where the natural gas invades into the drilling fluid to the ground, and large time lag exists, so that the gas reservoir parameters obtained from the ground are difficult to accurately correspond to the gas reservoir depth;
2. the analysis result has large error. The natural gas quantity obtained on the ground is the difference between the gas quantity produced by the gas reservoir and the dissolved and detained natural gas quantity in the shaft, and different gases have solubility difference, so that the ground trapped gas quantity and gas components have larger difference with the invasion quantity and components of the gas reservoir at the bottom of the shaft, and the components and content analysis is directly carried out according to the natural gas obtained on the ground, so that the error is larger.
Disclosure of Invention
Aiming at the defects of the prior art, the invention provides a gas reservoir state parameter while-drilling numerical inversion analysis method, which can effectively solve the problems in the prior art.
In order to realize the purpose, the technical scheme adopted by the invention is as follows:
a gas reservoir state parameter while-drilling numerical inversion analysis method comprises the following steps:
1) establishing a control equation, comprising the following steps:
11) establishing a single liquid phase shaft pressure field model; the formula is as follows:
in the formula:
is total pressure gradient, Pa/m;
Zvis vertical depth, m;
calculating the friction resistance pressure gradient of the liquid drilling fluid by adopting a Herschel-Bulkley rheological model;
12) the method for establishing the single-liquid-phase shaft heat exchange model comprises the following steps:
121) establishing drilling fluid model in drill string
122) Building drill stem pipe body model
123) Establishing an annulus drilling fluid model
124) Establishing a stratigraphic model
125) Establishing a well wall and stratum interface model
13) The method for establishing the gas-liquid two-phase flow heat exchange model comprises the following steps:
131) establishing a conservation of mass equation comprises:
establishing a drilling fluid mass conservation equation:
establishing a free gas mass conservation equation:
establishing a dissolved gas mass conservation equation:
132) establishing a conservation of momentum equation as follows:
133) establishing an energy conservation equation comprising:
establishing an energy conservation equation in the drill string:
establishing an annular energy conservation equation:
establishing a wellbore/formation interface energy conservation equation:
wherein f (t) is a dimensionless instantaneous heat loss function;
14) establishing a gas reservoir seepage model, comprising the following steps:
141) establishing a low-pressure gas non-Darcy unstable seepage model:
142) establishing a high-pressure gas non-Darcy unstable seepage relation model:
15) establishing a natural gas PVT model, wherein the 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) establishing a natural gas viscosity model, wherein the 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
μ＝1.500143567ρg+9.84848797×10-2Tr+1.407003797×10-3Mw-8.257944449×10-1
(16)
17) establishing a dissolution rate model, wherein the formula is as follows:
in the formula:
to take into account the saturated solubility of the three main components;
cbulbis the average solubility;
k is the mass transfer coefficient, in laminar flow: k is 0.026Re0.8Sc1/3D/(dout-dinner) And when in turbulent flow:
18) establishing a gas-liquid two-phase flow velocity model, wherein the formula is as follows:
vg＝C0vm+v∞ (18)
in the formula:
C0-a velocity profile factor; the ratio of the velocity at the center of the two-phase flow to the average flow velocity depends on the flow pattern, which includes: bubbly flow, slug flow, churn flow, annular flow;
vm-average velocity of the gas-liquid two-phase flow mixture;
v∞the relative liquid phase slip velocity of the gas phase depends on the flow pattern;
2) solving the division and node arrangement of the domain grid, and arranging the positions of temperature and pressure speed nodes by using a staggered grid technology, wherein the temperature nodes are positioned in the center of the grid, and the pressure and speed items are arranged at the boundary of the grid, so that discrete variables are ensured to strictly meet the continuity of a control equation and solution variables;
the axial coordinate points to the well bottom from the well head, the number of the axial grid is gradually increased from the well head to the well bottom, the number of the well head node is 1, the well bottom is added with a virtual node N +1 with the grid axial length of 0, the grid axial positions of the same node number in each discrete area are the same, and the change part of the flow area is arranged at the grid boundary;
in the oil gas diffusion migration simulation process, in a new time step length, the gas phase migration distance is not allowed to exceed the grid length, otherwise, the problem of calculation non-convergence is caused, namely the step length of the new time step is satisfied:
3) and (3) discretizing a control equation, comprising the following steps:
31) discretization of conservation of mass equation
Discrete mass conservation equations 7, 8 and 9 result:
liquid phase:
free gas phase:
dissolving a gas phase:
tn-1time-of-day trellis parameter and upstream i +1 node tnThe time parameter is known, so the solution of the mass conservation discrete equation is substantially tnSolving mass transmission items of liquid phase, dissolved gas and free gas at the downstream boundary of the time grid; the annular air cells are divided into three grid cells according to the components of the internal fluid phase:
firstly, mixing a gas phase and a liquid phase to form a grid unit;
a gas phase front edge grid, wherein one part of the grid is a gas phase and a liquid phase, and the other part of the grid is a pure liquid phase;
③ pure liquid phase grid; the grid on the upper part of the grid on which the gas phase front edge is positioned is a pure liquid phase grid;
311) gas-liquid two-phase grid downstream boundary mass transport
The mass of three substances flowing out of the boundary of the gas-liquid two-phase grid can be calculated according to the time step, the volume gas content of the last time step, the density of the last time step and the speed of the current time;
liquid phase:
free gas:
dissolved gas:
312) pure liquid phase mesh downstream boundary mass transport
Only the liquid phase which does not contain dissolved gas flows out of the boundary of the pure liquid phase grid and can be obtained by calculation according to the time step length and the speed of the current time;
liquid phase:
free gas:
dissolved gas:
313) gas phase leading edge mesh downstream boundary mass transport
tn-1Time tnThe gas phase velocity at the downstream boundary of the time is
A liquid phase velocity of
The pure liquid phase does not contain dissolved gas;
tn-1the length of the gas-liquid two-phase part at the time node i is LgPure liquid phase fraction length Δ Zi-LgThe volume gas fraction of the gas-liquid two-phase portion is
Density of liquid phase of
Density of gas phase of
Pure liquid phase density of rhol；
Passing of gas phase front
Then reaches the downstream boundary of the grid i, and the pure liquid phase velocity in the delta t moment is the grid gas-liquid two-phase average velocity, namely
①ΔtnAt ≦ Δ t, i.e.
ΔtnThe time period gas phase front does not exceed the downstream boundary, in which case only pure liquid phase is discharged, the mass discharged being:
liquid phase:
dissolved gas:
free gas:
②Δtn>at Δ t
ΔtnThe gas phase front exceeds the downstream boundary in the time interval, and the mass of the discharged substances is determined according to whether the discharged substances contain dissolved gas or not;
a.
Δtnthe volume of liquid phase discharged in the time period is not more than tn-1The pure liquid phase volume in the grid does not discharge dissolved gas at the moment;
liquid phase:
dissolved gas:
free gas:
b.
Δtnthe volume of liquid phase discharged in time period is more than tn-1The pure liquid phase volume in the time grid and the discharged liquid phase contain dissolved gas;
liquid phase:
dissolved gas:
free gas:
in the formula
Density of liquid phase containing dissolved gas at time of i grid tn-1, pLDensity of liquid phase without dissolved gas;
314) quality determination of grid dissolved gas
A dissolved gas mass calculation region in a delta t time period is from the upstream boundary of the grid to the gas phase front edge, the dissolution rate is calculated according to the material composition and the flow characteristics from the upstream boundary of the grid to the gas phase front edge region of the grid at the n-1 moment, and the dissolved gas generated by the gas invaded by the stratum in the time period is ignored;
the mass of dissolved gas in the stationary grid after Δ t time is:
the equation is a scalar nonlinear equation and is solved by iteration
And then, obtaining the grid dissolved gas mass increment in the delta t time period as follows:
32) discrete equation of conservation of momentum
The annular drilling fluid flows upwards, the direction of an axis coordinate is increased from top to bottom, the directions of the axis coordinate and the axis coordinate are opposite, and a velocity term is positive during dispersion, so that a momentum conservation equation in a discrete format is obtained:
in the grid containing the gas phase front edge, respectively calculating friction resistance pressure drop of a single liquid phase area and a gas-liquid two-phase area, and obtaining the total pressure drop of the grid through weighting summation;
33) energy conservation equation discretization
And (3) applying a specific enthalpy derivative thermodynamic theory to eliminate the temperature of a shaft/stratum interface to obtain an annulus node temperature calculation discrete equation represented by thermal resistance:
4) the gas reservoir parameter value inversion algorithm comprises the following steps:
41) partitioning a grid
Dividing grids according to a well structure and a drilling tool structure, and setting node arrays pNode [ N ] and pNodeNew [ N ]; storing last time point grid data in the pNode, and storing new time point data in the pNodeNew;
42) calculating grid vertical coordinates, and updating the pNode node array;
calculating the vertical coordinate of the center of the grid and the vertical length of the grid according to the well track;
43) analyzing the temperature field according to the single liquid phase flow and updating the pNode node array
44) Calculating the grid flow parameters and pressure according to the single liquid phase; the method comprises the steps of assigning 0 to the mass of a liquid phase, the outflow mass flow, the density of the liquid phase, the flow velocity of the liquid phase, the pressure gradient, the pressure, the gas content of a cross section, 0 to the flow velocity of a gas phase, 0 to the outflow mass of the gas phase, 0 to the mass of dissolved gas and 0 to the outflow mass of the dissolved gas, and updating a pNode node array;
45) determining gas-invaded node number M
Traversing the grid, comparing the axial position of the grid node with the depth of the gas layer, and determining the number M of the gas invasion node;
46) initial condition application of diffusion migration simulation
Diffusion migration simulation time t is 0, gas phase leading edge grid index
Gas phase length in gas phase front grid
Point gas invasion treatment;
synchronizing the grid variable pNodeNew [ N ] into a value of pNode [ N ];
47) determining a time step dt according to the grid length and the grid gas phase speed;
48) the time advances, t-t + dt,
49) gas cut upstream node pressure initialization
The gas cut upstream node pressure adjustment value pDiff is 0,
410) determining the pressure of the gas invasion upstream node,
411) according to the gas layer parameters, drilling parameters and
calculating gas invasion
412) Computing grid parameters based on grid type
Traversing the grids from bottom to top, determining the grid type according to the grid index number i, and calculating grid flow parameters;
413) Comparing the pressure of the well mouth node with the atmospheric pressure, and calculating the relative error epsilonpDetermining a bottom boundary pressure adjustment value
414) If epsilonp>10-3Go to step 410), repeat until εp<10-3The grid parameter calculation at the new time point is finished;
415) analyzing a new temperature field according to the grid flow parameters and the pressure at the new moment;
416) updating all annulus mesh data pNode to pNodeNew;
417) if the time reaches the analysis time, the analysis is finished, otherwise, the step 47) is carried out repeatedly.
Further, a gas reservoir state parameter analysis simulator is compiled by applying a gas reservoir state parameter while-drilling numerical inversion interpretation method, and comprises the following functional modules:
ReadKickData
ReadKickData is a data input function of an analysis simulator, and completes the input of all data of a simulated well, including a well body structure, a drilling tool structure, a well track, the vertical temperature distribution of a stratum and reservoir parameters: temperature, pressure, porosity, permeability, depth, supply radius, thermophysical parameters of crude oil and natural gas, thermophysical parameters of cement, formation, drill string, etc., drilling parameters: discharge capacity, rate of penetration, torque, drill bit nozzle, simulation parameters: axial grid dimension, simulation period, whether crushed gas is considered or not, and whether a heat source is considered or not;
KickFlowBehavior
the KickFlowBehavior function is a general function module of the simulator, and other functions are assembled according to a program flow chart to finish the functions of drilling oil gas diffusion and migration simulation analysis and data storage;
TemperatureAnalysis
performing temperature field analysis on the drilling shaft by the TemperatureAnalysis function according to the grid flow parameters, and updating the temperature calculation result into a grid node array pNodeNew;
KickSimulating
performing well drilling oil gas diffusion and migration simulation analysis according to the grid node temperature and the boundary conditions of the bottom hole and the wellhead by the KickSimuling function, and updating the calculation result of the flow parameter into a grid node array pNodeNew;
GridGeneration
the GridGeneration function carries out axial segmentation and grid division on a temperature and pressure field solving domain according to a well structure and a drill string structure, and stores axial geometric information of grid nodes and geometric and medium information of a radial heat exchange object;
NormalTPNode
the normalTPnode function analyzes grid flow parameters according to a grid flow parameter calculation method that an inlet and an outlet are both gas-liquid two-phase flow;
GasFrontNode
the GasFrontNode function analyzes grid flow parameters according to a gas phase leading edge grid algorithm and updates a gas phase leading edge grid index and gas phase zone length data;
VerticalCoordinate
calculating the vertical coordinate of the center of the grid and the vertical length of the grid according to the well track; the function requires calling a DirectionParaCal function to calculate the vertical depth according to the well depth and the well inclination angle;
TOriginGeneration
the torrginggeneration function interpolates the vertical temperature distribution data of the formation to generate the original temperature at the vertical depth of the node;
Ini
applying initial conditions to solving intra-domain grid node variables by using an Ini function, and assigning initial values;
HeatResistanceSA
calculating heat exchange resistance between drilling fluid in the drill bit and annular drilling fluid; the function is realized by calling HPipe and HANnu functions to calculate the forced convection heat transfer coefficient of the liquid phase in the drill main and the gas-liquid two-phase flow in the annular space;
HeatResistanceAE
calculating heat exchange resistance between the annular gas-liquid two-phase flow drilling fluid and the stratum; the function is realized by calling a HANnu function to calculate the forced convection heat transfer coefficient of gas-liquid two-phase flow in the annulus and an instantaneous formation heat loss function to calculate the heat transfer thermal resistance from the shaft to the formation;
HPipe
the function of the HPipe function is to calculate the pipe flow forced convection heat transfer coefficient of the drilling fluid in the drill bit and the inner wall of the drill column;
HAnnu
the function of the HANnu function is to calculate the forced convection heat transfer coefficient of the annular flow of the annular gas-liquid two-phase flow drilling fluid and the well wall;
PPTP
calculating a total integration function of thermodynamic properties of the medium at a given temperature and pressure according to the type of the medium, wherein crude oil, water and natural gas are required to be used according to specific types of the medium; calling a thermophysical property calculation function of a corresponding medium to complete the function of the medium;
DPDZ_TwoPhase
calculating the flow parameters and pressure drop gradient of the gas-liquid two-phase flow according to the gas-liquid two-phase flow mechanism model, wherein the flow parameters of the gas-liquid two-phase flow comprise: gas phase velocity, liquid phase velocity, cross section gas content, gas phase density and flow pattern;
TInDrillstem
the TInDrillstem function iteratively calculates the node temperature of the drilling fluid in the drill bit according to a single liquid phase energy conservation equation;
TAnnulus
the TANNILUUS function iteratively calculates the node temperature of the annular drilling fluid according to the energy conservation equation of the gas-liquid two-phase flow;
calculating the constant pressure specific heat of the natural gas by using CpGas;
CpOil, calculating the constant pressure specific heat of the crude oil;
CpWater, calculating the constant pressure specific heat of formation water;
DenGas, calculating the natural gas density;
DenOilSatured, calculating the saturated crude oil density;
DenOilUnsated, calculating the density of the unsaturated crude oil;
RsStanding, calculating the dissolved gas-oil ratio of the crude oil;
surface tension, calculating the oil-water interfacial tension;
thcon gas, calculating the natural gas thermal conductivity;
thcon oil, calculating crude oil thermal conductivity;
ThconWater, calculating formation water thermal conductivity;
ViscGas, calculating natural gas viscosity;
ViscOil, calculate crude oil viscosity;
ViscWater, calculate the formation water viscosity.
Further, setting boundary conditions and initial conditions of the gas reservoir state parameter while-drilling numerical inversion interpretation method:
boundary condition
At the drill string inlet:
Qm＝C (45-1)
T1(z＝0)＝TE (45-2)
wherein T isEIs the drilling fluid inlet temperature, QmThe drilling fluid discharge capacity;
② the annular outlet:
p＝patm (46)
(iii) the junction of the drill string at the bottom of the well and the annulus, the temperature of the drilling fluid is the same, there are
Initial conditions
The initial condition of the temperature field is the calculated wellbore temperature under the condition of single liquid phase, namely:
Ti(z,0)＝Ti,st(z) (48)
the initial condition of the pressure field is calculated wellbore pressure under the condition of single liquid phase, namely:
pi(z,0)＝pi,st(z) (49)
and thirdly, when gas invasion starts, no gas phase exists in each mesh of the annulus:
αi(z,0)＝0 (50)。
compared with the prior art, the invention has the advantages that:
1. the gas reservoir state parameters under the condition of uninterrupted drilling operation are accurately analyzed in real time while drilling, and the gas well development benefit is improved. The method integrates 18 control equations, comprehensively considers the seepage of the gas reservoir to the shaft, the gas-liquid two-phase flow and heat exchange in the shaft, the dissolution of the gas phase in the drilling fluid and the influence of temperature and pressure on the thermophysical properties of the gas phase, and can obtain the gas reservoir state parameter along the well depth profile by inversion according to the ground gas measurement parameter by applying the method, so that the gas reservoir state parameter is finer.
2. The well kick risk is identified in time, the physical property and flow state parameter distribution of gas phase along the shaft can be obtained according to the ground gas measurement parameters, and the well bottom gas invasion amount can be obtained, so that the well kick degree and the well kick risk can be identified accurately, and the method has important significance for preventing the well blowout risk and ensuring the safety of drilling operation.
Drawings
FIG. 1 is a schematic diagram of mesh partitioning according to an embodiment of the present invention;
FIG. 2 is a schematic diagram of the positions of the nodes i and i +1 grid parameters according to the embodiment of the present invention;
FIG. 3 is a schematic view of gas phase front grid mass transport according to an embodiment of the present invention;
fig. 4 is a block diagram of a gassreservoir analysis simulator according to an embodiment of the present invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is described in further detail below with reference to the accompanying drawings and examples.
A gas reservoir state parameter while-drilling numerical inversion analysis method comprises the following parts:
controlling equation
1) Single liquid phase shaft pressure field model
In the formula:
is total pressure gradient, Pa/m;
Zvis vertical depth, m;
And calculating the friction resistance pressure gradient of the liquid drilling fluid by adopting a Herschel-Bulkley rheological model.
2) Single liquid phase shaft heat exchange model
Drilling fluid in drill string
Drill string pipe body
Annular drilling fluid
Stratum
Well wall and formation interface
3) Gas-liquid two-phase flow heat exchange model
Equation of conservation of mass
Drilling fluid mass conservation equation:
free gas mass conservation equation:
solution gas mass conservation equation:
equation of conservation of momentum
Equation of conservation of energy
Conservation of energy in the drill string equation:
annulus energy conservation equation:
wellbore/formation interface energy conservation equation:
where f (t) is a dimensionless instantaneous heat loss function.
4) Gas reservoir seepage model
Low-pressure gas non-darcy unstable seepage model:
high-pressure gas non-Darcy unstable seepage relation model:
5) natural gas PVT model:
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
μ＝1.500143567ρg+9.84848797×10-2Tr+1.407003797×10-3Mw-8.257944449×10-1 (16)
7) dissolution rate model
In the formula:
to take into account the saturated solubility of the three main components;
cbulbis the average solubility;
k is the mass transfer coefficientAnd when laminar flow is carried out: k is 0.026Re0.8Sc13D/(dout-dinner) And when in turbulent flow:
8) gas-liquid two-phase flow velocity model
vg＝C0vm+v∞ (18)
In the formula:
C0-velocity profile factor. The ratio of the velocity at the center of the two-phase flow to the average flow velocity depends on the flow pattern (bubbly, slug, churn, annular);
vm-average velocity of the gas-liquid two-phase flow mixture;
v∞the relative liquid phase slip velocity of the gas phase depends on the flow pattern.
Second, solving the domain grid division and node arrangement
The invention uses the finite volume method to disperse the mass conservation equation and the energy conservation equation, and uses the staggered grid technology to arrange the positions of the temperature node and the pressure speed node, wherein the temperature node is positioned in the center of the grid, and the pressure and speed items are arranged at the boundary of the grid, thereby ensuring that the discrete variables strictly meet the continuity of the control equation and the solved variables.
The axial coordinate points to the well bottom from the well head, the number of the axial grid is gradually increased from the well head to the well bottom, the number of the well head node is 1, the well bottom is added with a virtual node N +1 with the grid axial length being 0, the grid axial positions of the same node number in each discrete area are the same, and the change part of the flow area is arranged at the grid boundary. The meshing is shown in fig. 1.
In the oil gas diffusion migration simulation process, in a new time step length, the gas phase migration distance is not allowed to exceed the grid length, otherwise, the problem of calculation non-convergence is caused, namely the step length of the new time step is satisfied:
three, control equation discretization
Taking the annular node i as an example, the control equation is subjected to time and space dispersion. The control equation relates to physical quantities: p, T, pm、ρg、vm、vg、α、
f1. In the annular space, the drilling fluid flows from bottom to top, the node i parameter is influenced by the node i +1 parameter, and the annular space node i and the node i +1 grid parameter are arranged as shown in fig. 2. The pressure may be considered to be linearly distributed in the grid. And performing control equation discretization according to the control equation.
1) Discretization of conservation of mass equation
Discrete mass conservation equations (7), (8), and (9) yield:
liquid phase:
free gas phase:
dissolving a gas phase:
tn-1time of day trellis parameter and upstream (i +1) node tnThe time parameter is known, so the solution of the mass conservation discrete equation is substantially tnAnd solving mass transfer terms of liquid phase, dissolved gas and free gas at the downstream boundary of the time grid. The annular air cells are divided into three grid cells according to the components of the internal fluid phase:
firstly, mixing a gas phase and a liquid phase to form a grid unit;
a gas phase front edge grid, wherein one part of the grid is a gas phase and a liquid phase, and the other part of the grid is a pure liquid phase;
③ pure liquid phase grid. The grid above the grid on which the gas phase front edge is located is a pure liquid phase grid.
(1) Gas-liquid two-phase grid downstream boundary mass transport
The mass of three substances flowing out of the boundary of the gas-liquid two-phase grid can be calculated according to the time step, the volume gas content of the last time step, the density of the last time step and the speed of the current time.
Liquid phase:
free gas:
dissolved gas:
(2) pure liquid phase mesh downstream boundary mass transport
Only the liquid phase without dissolved gas flows out of the boundary of the pure liquid phase grid, and the liquid phase can be calculated according to the time step and the speed of the current time.
Liquid phase:
free gas:
dissolved gas:
(3) gas phase leading edge mesh downstream boundary mass transport
A schematic diagram of mass transport of the gas phase front grid liquid phase, dissolved gas, and free gas is shown in fig. 3.
tn-1Time tnThe gas phase velocity at the downstream boundary of the time is
A liquid phase velocity of
The pure liquid phase contains no dissolved gas.
tn-1The length of the gas-liquid two-phase part at the time node i is LgPure liquid phase fraction length of (Δ Z)i-Lg) The volume gas fraction of the gas-liquid two-phase portion is
Density of liquid phase of
Density of gas phase of
Pure liquid phase density of rhol。
Passing of gas phase front
Then reaches the downstream boundary of the grid i, and the pure liquid phase velocity in the delta t moment is the grid gas-liquid two-phase average velocity, namely
①ΔtnAt ≦ Δ t, i.e.
ΔtnThe time period gas phase front does not exceed the downstream boundary, in which case only pure liquid phase is discharged, the mass discharged being:
liquid phase:
dissolved gas:
free gas:
②Δtn>at Δ t
ΔtnThe time period gas phase front exceeds the downstream boundary, and the mass of the discharged substance is determined according to whether the discharged substance contains dissolved gas or not.
a.
ΔtnThe volume of liquid phase discharged in the time period is not more than tn-1The pure liquid phase volume in the grid at the moment does not discharge the dissolved gas.
Liquid phase:
dissolved gas:
free gas:
b.
Δtnthe volume of liquid phase discharged in time period is more than tn-1The volume of pure liquid phase in the time grid, the exhaust liquid phase, contains dissolved gas.
Liquid phase:
dissolved gas:
free gas:
in the formula
Density of liquid phase containing dissolved gas at time of i grid tn-1, pLIs the density of the liquid phase without dissolved gas.
(4) Quality determination of grid dissolved gas
And (3) obtaining a dissolved gas mass obtaining area from the upstream boundary of the grid to the gas phase front edge in the delta t time period, obtaining the dissolution rate according to the material composition and flow characteristics from the upstream boundary of the grid to the gas phase front edge area of the grid at the n-1 moment, and neglecting the dissolved gas generated by the gas invaded from the stratum in the time period.
The mass of dissolved gas in the stationary grid after Δ t time is:
the equation is a scalar nonlinear equation and is solved by iteration
And then, obtaining the grid dissolved gas mass increment in the delta t time period as follows:
2) discrete equation of conservation of momentum
The annular drilling fluid flows upwards, the direction of an axis coordinate is increased from top to bottom, the directions of the axis coordinate and the axis coordinate are opposite, and a velocity term is positive during dispersion, so that a momentum conservation equation in a discrete format is obtained:
in the grid containing the gas phase front edge, the friction resistance pressure drop is respectively calculated for a single liquid phase area and a gas-liquid phase area, and the total pressure drop of the grid is obtained through weighting summation.
3) Energy conservation equation discretization
And (3) applying a specific enthalpy derivative thermodynamic theory to eliminate the temperature of a shaft/stratum interface to obtain an annulus node temperature calculation discrete equation represented by thermal resistance:
in the formula (I), the compound is shown in the specification,
the total thermal resistance of the annulus drilling fluid to the formation.
Gas phase leading edge mesh tracking
The following variables will be used for the gas phase front grid tracking task.
Ig-a grid comprising a gas front;
Lgthe position of the gas front relative to the grid IgThe upstream edge of (a).
By moving the leading edge and the downstream velocity, these variables are updated as follows:
if it is not
Then
For a mesh containing a gas front, it is assumed that both free and dissolved gases are uniformly distributed behind the gas front, otherwise it is assumed that the gases are uniformly distributed throughout the mesh. The velocity of the drilling fluid will change as the gas front passes a grid boundary in a time step, since the velocity of the mixture is considered constant in a time step. So that the time step is split, with v before the gas phase front passesm＝vMAnd followed by vm＝(vM-αvg)/(1-α)。
Fifthly, boundary conditions and initial conditions
1) Boundary condition
At the drill string inlet:
Qm＝C (45-1)
T1(z＝0)＝TE (45-2)
wherein T isEIs the drilling fluid inlet temperature, QmIs the drilling fluid displacement.
② the annular outlet:
p＝patm (46)
(iii) the junction of the drill string at the bottom of the well and the annulus, the temperature of the drilling fluid is the same, there are
2) Initial conditions
(1) The initial conditions of the temperature field are the wellbore temperatures that have been calculated under single liquid phase conditions, i.e.:
Ti(z,0)＝Ti,st(z) (48)
(2) the initial conditions of the pressure field are the wellbore pressures that have been calculated under single liquid phase conditions, i.e.:
pi(z,0)＝pi,st(z) (49)
(3) when gas invasion begins, no gas phase exists in each grid of the annulus:
αi(z,0)＝0 (50)
sixth, algorithm
When the numerical inversion analysis of the gas reservoir parameters is carried out, firstly, the gas reservoir parameters are assumed, the ground monitoring parameters are calculated according to the following algorithm, if the error between the calculated value and the monitoring value is large, the gas reservoir parameters are adjusted, the numerical simulation is repeated until the ground calculated value and the monitoring value meet the convergence condition, and the real gas reservoir state parameters are obtained. The numerical simulation algorithm assuming the gas reservoir parameters is as follows.
Because the temperature, the pressure field, the temperature and the pressure data of the gas reservoir and the like of the annular gas-liquid two-phase flow are mutually dependent, the annular gas-liquid two-phase flow needs to be advanced according to the time step from the initial value until the analysis time is reached. Each time step needs to start from the value of the last time point, and each convergence solution is obtained through iterative calculation.
(1) Partitioning a grid
And dividing grids according to the well structure and the drilling tool structure, and setting node arrays pNode [ N ] and pNodeNew [ N ]. The pNode stores last time point grid data, and the pNodeNew stores new time point data.
(2) Calculating grid vertical coordinates, and updating the pNode node array;
the vertical coordinate of the center of the grid and the vertical length of the grid are calculated from the wellbore trajectory.
(3) Analyzing the temperature field according to the single liquid phase flow and updating the pNode node array
(4) And calculating the grid flow parameters and pressure according to the single liquid phase. The method comprises the steps of assigning 0 to the mass of a liquid phase, the outflow mass flow, the density of the liquid phase, the flow velocity of the liquid phase, the pressure gradient, the pressure, the gas content of a cross section, 0 to the flow velocity of a gas phase, 0 to the outflow mass of the gas phase, 0 to the mass of dissolved gas and 0 to the outflow mass of the dissolved gas, and updating the pNode node array.
(5) Determining gas-invaded node number M
And traversing the grid, comparing the axial position of the grid node with the depth of the gas layer, and determining the number M of the gas invasion node.
(6) Initial condition application of diffusion migration simulation
Diffusion migration simulation time t is 0, gas phase leading edge grid index
Gas phase length in gas phase front grid
(Point gas invasion treatment).
The grid variable pNodeNew [ N ] is synchronized to the value of pNode [ N ].
(7) Determining a time step dt according to the grid length and the grid gas phase speed;
(8) the time advances, t-t + dt,
(9) gas cut upstream node pressure initialization
The gas cut upstream node pressure adjustment value pDiff is 0,
(10) determining the pressure of the gas invasion upstream node,
(11) according to the gas layer parameters, drilling parameters and
calculating gas invasion
(12) Computing grid parameters based on grid type
And traversing the grids from bottom to top, determining the grid type according to the grid index number i, and calculating the grid flow parameters.
(13) Comparing the pressure of the well mouth node with the atmospheric pressure, and calculating the relative error epsilonpDetermining a bottom boundary pressure adjustment value
(14) If epsilonp>10-3Turning to the step (10), and repeatedly executing until epsilonp<10-3The grid parameter calculation at the new time point is finished;
(15) analyzing a new temperature field according to the grid flow parameters and the pressure at the new moment;
(16) updating all annulus mesh data pNode to pNodeNew;
(17) if the time reaches the analysis time, the analysis is finished, otherwise, the step (7) is carried out repeatedly.
Seventh, program module
The gas reservoir state parameter analysis simulator GasRevoruirAnalyzis compiled by the method and consists of 63 functions, and the main module structure is shown in FIG. 4. Its main functional blocks (functions) are as follows.
(1)ReadKickData
The ReadKickData is a data input function of an analysis simulator, completes the input of all data of a simulated well, and mainly comprises a well body structure, a drilling tool structure, a well track, the vertical temperature distribution of a stratum and reservoir parameters: temperature, pressure, porosity, permeability, depth, supply radius, thermophysical parameters of crude oil and natural gas, thermophysical parameters of cement, formation, drill string, etc., drilling parameters: discharge capacity, rate of penetration, torque, drill bit nozzle, simulation parameters: axial grid dimensions, simulation period, whether crushed gas is considered, whether heat source is considered.
(2)KickFlowBehavior
The KickFlowBehavior function is a general function module of the simulator, and other functions are assembled according to a program flow chart to finish the functions of drilling oil gas diffusion and migration simulation analysis and data storage.
(3)TemperatureAnalysis
And the TemperatureAnalysis function analyzes the temperature field of the well drilling shaft according to the grid flow parameters, and updates the temperature calculation result to a grid node array pNodeNew.
(4)KickSimulating
And (3) performing well drilling oil gas diffusion and migration simulation analysis according to the grid node temperature and the boundary conditions of the bottom hole and the well head by using the KickSimuling function, and updating the calculation result of the flow parameter into a grid node array pNodeNew.
(5)GridGeneration
The GridGeneration function carries out axial segmentation and grid division on a temperature and pressure field solving domain according to a well body structure and a drill string structure, and stores axial geometric information of grid nodes and geometric and medium information of a radial heat exchange object (media in a grid control body comprise 6 medium types such as single-liquid-phase drilling fluid, gas-liquid two-phase drilling fluid, steel, cement, stratum, packer fluid and the like, and each medium has corresponding thermodynamic physical properties).
(6)NormalTPNode
And analyzing the grid flow parameters by the normalTPnode function according to a grid flow parameter calculation method that the inlet and the outlet are both gas-liquid two-phase flows.
(7)GasFrontNode
The GasFrontNode function analyzes the grid flow parameters according to the gas phase leading edge grid algorithm and updates the gas phase leading edge grid index and the gas phase zone length data.
(8)VerticalCoordinate
The vertical coordinate of the center of the grid and the vertical length of the grid are calculated from the wellbore trajectory. The function calls for the DirectionParaCal function to calculate the sag from the well depth and the angle of the well.
(9)TOriginGeneration
The torrginggeneration function interpolates from the vertical temperature distribution data of the formation to produce the original temperature at the vertical depth of the node.
(10)Ini
And applying initial conditions to the variables of the grid nodes in the solution domain by the Ini function, and assigning initial values.
(11)HeatResistanceSA
And calculating the heat exchange thermal resistance between the drilling fluid in the drill bit main and the annular drilling fluid. The function is realized by calling HPipe and HANnu functions to calculate the forced convection heat transfer coefficient of the liquid phase in the drill main and the gas-liquid two-phase flow in the annular space.
(12)HeatResistanceAE
And calculating the heat exchange thermal resistance between the annular gas-liquid two-phase flow drilling fluid and the stratum. The function is realized by calling a HANnu function to calculate the forced convection heat transfer coefficient of gas-liquid two-phase flow in the annulus and an instantaneous formation heat loss function to calculate the heat transfer thermal resistance from the shaft to the formation.
(13)HPipe
The function of the HPipe function is to calculate the pipe flow forced convection heat transfer coefficient of the drilling fluid in the drill main and the inner wall of the drill string.
(14)HAnnu
The function of the HANnu function is to calculate the forced convection heat transfer coefficient of the annular flow of the annular gas-liquid two-phase flow drilling fluid and the well wall.
(15)PPTP
The method comprises the steps of calculating a total integrated function of thermodynamic properties of a medium at a given temperature and pressure according to the type of the medium, and calling a thermophysical property calculation function of the corresponding medium according to the specific type of the medium (crude oil, water and natural gas) to complete the function of the medium.
(16)DPDZ_TwoPhase
And calculating the flow parameters (gas phase velocity, liquid phase velocity, section gas content, gas phase density and flow pattern) and pressure drop gradient of the gas-liquid two-phase flow according to the gas-liquid two-phase flow mechanism model.
(17)TInDrillstem
And the TInDrillstem function iteratively calculates the drilling fluid node temperature in the drill bit according to a single liquid phase energy conservation equation.
(18)TAnnulus
And the TANNILUUS function iteratively calculates the node temperature of the annular drilling fluid according to the energy conservation equation of the gas-liquid two-phase flow.
(19) CpGas. And calculating the specific heat of the natural gas at constant pressure.
(20) CpOil. And calculating the constant pressure specific heat of the crude oil.
(21) CpWater. And calculating the constant pressure specific heat of the formation water.
(22) DenGas. And calculating the natural air density.
(23) DenOilSatured. The saturated crude density was calculated.
(24) DenOilUnsated. And calculating the density of the unsaturated crude oil.
(25) RsStanding. And calculating the dissolved gas-oil ratio of the crude oil.
(26) Surface tension. And calculating the oil-water interfacial tension.
(27) ThconGas. And calculating the thermal conductivity of the natural gas.
(28) Thcon oil. And calculating the thermal conductivity of the crude oil.
(29) Thcon water. And calculating the formation water thermal conductivity.
(30) ViscGas. The natural gas viscosity is calculated.
(31) ViscOil. The crude oil viscosity was calculated.
(32) ViscWater. The formation water viscosity is calculated.
It will be appreciated by those of ordinary skill in the art that the examples described herein are intended to assist the reader in understanding the manner in which the invention is practiced, and it is to be understood that the scope of the invention is not limited to such specifically recited statements and examples. Those skilled in the art can make various other specific changes and combinations based on the teachings of the present invention without departing from the spirit of the invention, and these changes and combinations are within the scope of the invention.
Claims (3)
1. A gas reservoir state parameter while-drilling numerical inversion analysis method is characterized by comprising the following steps:
1) establishing a control equation, comprising the following steps:
11) establishing a single liquid phase shaft pressure field model; the formula is as follows:
in the formula:
is total pressure gradient, Pa/m;
Zvis vertical depth, m;
calculating the friction resistance pressure gradient of the liquid drilling fluid by adopting a Herschel-Bulkley rheological model;
12) the method for establishing the single-liquid-phase shaft heat exchange model comprises the following steps:
121) establishing drilling fluid model in drill string
122) Building drill stem pipe body model
123) Establishing an annulus drilling fluid model
124) Establishing a stratigraphic model
125) Establishing a well wall and stratum interface model
13) The method for establishing the gas-liquid two-phase flow heat exchange model comprises the following steps:
131) establishing a conservation of mass equation comprises:
establishing a drilling fluid mass conservation equation:
establishing a free gas mass conservation equation:
establishing a dissolved gas mass conservation equation:
132) establishing a conservation of momentum equation as follows:
133) establishing an energy conservation equation comprising:
establishing an energy conservation equation in the drill string:
establishing an annular energy conservation equation:
establishing a wellbore/formation interface energy conservation equation:
wherein f (t) is a dimensionless instantaneous heat loss function;
14) establishing a gas reservoir seepage model, comprising the following steps:
141) establishing a low-pressure gas non-Darcy unstable seepage model:
142) establishing a high-pressure gas non-Darcy unstable seepage relation model:
15) establishing a natural gas PVT model, wherein the 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) establishing a natural gas viscosity model, wherein the 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
μ＝1.500143567ρg+9.84848797×10-2Tr+1.407003797×10-3Mw-8.257944449×10-1 (16)
17) establishing a dissolution rate model, wherein the formula is as follows:
in the formula:
to take into account the saturated solubility of the three main components;
cbulbis the average solubility;
k is the mass transfer coefficient, in laminar flow: k is 0.026Re0.8Sc1/3D/(dout-dinner) And when in turbulent flow:
18) establishing a gas-liquid two-phase flow velocity model, wherein the formula is as follows:
vg＝C0vm+v∞ (18)
in the formula:
C0-a velocity profile factor; the ratio of the velocity at the center of the two-phase flow to the average flow velocity depends on the flow pattern, which includes: bubbly flow, slug flow, churn flow, annular flow;
vm-average velocity of the gas-liquid two-phase flow mixture;
v∞the relative liquid phase slip velocity of the gas phase depends on the flow pattern;
2) solving the division and node arrangement of the domain grid, and arranging the positions of temperature and pressure speed nodes by using a staggered grid technology, wherein the temperature nodes are positioned in the center of the grid, and the pressure and speed items are arranged at the boundary of the grid, so that discrete variables are ensured to strictly meet the continuity of a control equation and solution variables;
the axial coordinate points to the well bottom from the well head, the number of the axial grid is gradually increased from the well head to the well bottom, the number of the well head node is 1, the well bottom is added with a virtual node N +1 with the grid axial length of 0, the grid axial positions of the same node number in each discrete area are the same, and the change part of the flow area is arranged at the grid boundary;
in the oil gas diffusion migration simulation process, in a new time step length, the gas phase migration distance is not allowed to exceed the grid length, otherwise, the problem of calculation non-convergence is caused, namely the step length of the new time step is satisfied:
3) and (3) discretizing a control equation, comprising the following steps:
31) discretization of conservation of mass equation
Discrete mass conservation equations 7, 8 and 9 result:
liquid phase:
free gas phase:
dissolving a gas phase:
tn-1time-of-day trellis parameter and upstream i +1 node tnThe time parameter is known, so the solution of the mass conservation discrete equation is substantially tnSolving mass transmission items of liquid phase, dissolved gas and free gas at the downstream boundary of the time grid; the annular air cells are divided into three grid cells according to the components of the internal fluid phase:
firstly, mixing a gas phase and a liquid phase to form a grid unit;
a gas phase front edge grid, wherein one part of the grid is a gas phase and a liquid phase, and the other part of the grid is a pure liquid phase;
③ pure liquid phase grid; the grid on the upper part of the grid on which the gas phase front edge is positioned is a pure liquid phase grid;
311) gas-liquid two-phase grid downstream boundary mass transport
The mass of three substances flowing out of the boundary of the gas-liquid two-phase grid can be calculated according to the time step, the volume gas content of the last time step, the density of the last time step and the speed of the current time;
liquid phase:
free gas:
dissolved gas:
312) pure liquid phase mesh downstream boundary mass transport
Only the liquid phase which does not contain dissolved gas flows out of the boundary of the pure liquid phase grid and can be obtained by calculation according to the time step length and the speed of the current time;
liquid phase:
free gas:
dissolved gas:
313) gas phase leading edge mesh downstream boundary mass transport
tn-1Time tnThe gas phase velocity at the downstream boundary of the time is
A liquid phase velocity of
The pure liquid phase does not contain dissolved gas;
tn-1the length of the gas-liquid two-phase part at the time node i is LgPure liquid phase fraction length Δ Zi-LgThe volume gas fraction of the gas-liquid two-phase portion is
Density of liquid phase of
Density of gas phase of
Pure liquid phase density of rhol；
Passing of gas phase front
Then reaches the downstream boundary of the grid i, and the pure liquid phase velocity in the delta t moment is the grid gas-liquid two-phase average velocity, namely
①ΔtnAt ≦ Δ t, i.e.
ΔtnThe time period gas phase front does not exceed the downstream boundary, in which case only pure liquid phase is discharged, the mass discharged being:
liquid phase:
dissolved gas:
free gas:
②Δtn>at Δ t
ΔtnThe gas phase front exceeds the downstream boundary in the time interval, and the mass of the discharged substances is determined according to whether the discharged substances contain dissolved gas or not;
a.
Δtnthe volume of liquid phase discharged in the time period is not more than tn-1The pure liquid phase volume in the grid does not discharge dissolved gas at the moment;
liquid phase:
dissolved gas:
free gas:
b.
Δtnthe volume of liquid phase discharged in time period is more than tn-1The pure liquid phase volume in the time grid and the discharged liquid phase contain dissolved gas;
liquid phase:
dissolved gas:
free gas:
in the formula
Density of liquid phase containing dissolved gas at time of i grid tn-1, pLDensity of liquid phase without dissolved gas;
314) quality determination of grid dissolved gas
A dissolved gas mass calculation region in a delta t time period is from the upstream boundary of the grid to the gas phase front edge, the dissolution rate is calculated according to the material composition and the flow characteristics from the upstream boundary of the grid to the gas phase front edge region of the grid at the n-1 moment, and the dissolved gas generated by the gas invaded by the stratum in the time period is ignored;
the mass of dissolved gas in the stationary grid after Δ t time is:
the equation is a scalar nonlinear equation and is solved by iteration
And then, obtaining the grid dissolved gas mass increment in the delta t time period as follows:
32) discrete equation of conservation of momentum
The annular drilling fluid flows upwards, the direction of an axis coordinate is increased from top to bottom, the directions of the axis coordinate and the axis coordinate are opposite, and a velocity term is positive during dispersion, so that a momentum conservation equation in a discrete format is obtained:
in the grid containing the gas phase front edge, respectively calculating friction resistance pressure drop of a single liquid phase area and a gas-liquid two-phase area, and obtaining the total pressure drop of the grid through weighting summation;
33) energy conservation equation discretization
And (3) applying a specific enthalpy derivative thermodynamic theory to eliminate the temperature of a shaft/stratum interface to obtain an annulus node temperature calculation discrete equation represented by thermal resistance:
4) the gas reservoir parameter value inversion algorithm comprises the following steps:
41) partitioning a grid
Dividing grids according to a well structure and a drilling tool structure, and setting node arrays pNode [ N ] and pNodeNew [ N ]; storing last time point grid data in the pNode, and storing new time point data in the pNodeNew;
42) calculating grid vertical coordinates, and updating the pNode node array;
calculating the vertical coordinate of the center of the grid and the vertical length of the grid according to the well track;
43) analyzing the temperature field according to the single liquid phase flow and updating the pNode node array
44) Calculating the grid flow parameters and pressure according to the single liquid phase; the method comprises the steps of assigning 0 to the mass of a liquid phase, the outflow mass flow, the density of the liquid phase, the flow velocity of the liquid phase, the pressure gradient, the pressure, the gas content of a cross section, 0 to the flow velocity of a gas phase, 0 to the outflow mass of the gas phase, 0 to the mass of dissolved gas and 0 to the outflow mass of the dissolved gas, and updating a pNode node array;
45) determining gas-invaded node number M
Traversing the grid, comparing the axial position of the grid node with the depth of the gas layer, and determining the number M of the gas invasion node;
46) initial condition application of diffusion migration simulation
Diffusion migration simulation time t is 0, gas phase leading edge grid index
Gas phase length in gas phase front grid
Point gas invasion treatment;
synchronizing the grid variable pNodeNew [ N ] into a value of pNode [ N ];
47) determining a time step dt according to the grid length and the grid gas phase speed;
48) the time advances, t-t + dt,
49) gas cut upstream node pressure initialization
The gas cut upstream node pressure adjustment value pDiff is 0,
410) determining the pressure of the gas invasion upstream node,
411) according to the gas layer parameters, drilling parameters and
calculating gas invasion
412) Computing grid parameters based on grid type
Traversing the grids from bottom to top, determining the grid type according to the grid index number i, and calculating grid flow parameters;
413) Comparing the pressure of the well mouth node with the atmospheric pressure, and calculating the relative error epsilonpDetermining a bottom boundary pressure adjustment value
414) If epsilonp>10-3Go to step 410), repeat until εp<10-3The grid parameter calculation at the new time point is finished;
415) analyzing a new temperature field according to the grid flow parameters and the pressure at the new moment;
416) updating all annulus mesh data pNode to pNodeNew;
417) if the time reaches the analysis time, the analysis is finished, otherwise, the step 47) is carried out repeatedly.
2. The method for numerical inversion analysis while drilling of the gas reservoir state parameters as recited in claim 1, wherein: the gas reservoir state parameter analysis simulator is compiled by applying a gas reservoir state parameter while-drilling numerical inversion interpretation method, and comprises the following functional modules:
ReadKickData
ReadKickData is a data input function of an analysis simulator, and completes the input of all data of a simulated well, including a well body structure, a drilling tool structure, a well track, the vertical temperature distribution of a stratum and reservoir parameters: temperature, pressure, porosity, permeability, depth, supply radius, thermophysical parameters of crude oil and natural gas, thermophysical parameters of cement, formation, drill string, etc., drilling parameters: discharge capacity, rate of penetration, torque, drill bit nozzle, simulation parameters: axial grid dimension, simulation period, whether crushed gas is considered or not, and whether a heat source is considered or not;
KickFlowBehavior
the KickFlowBehavior function is a general function module of the simulator, and other functions are assembled according to a program flow chart to finish the functions of drilling oil gas diffusion and migration simulation analysis and data storage;
TemperatureAnalysis
performing temperature field analysis on the drilling shaft by the TemperatureAnalysis function according to the grid flow parameters, and updating the temperature calculation result into a grid node array pNodeNew;
KickSimulating
performing well drilling oil gas diffusion and migration simulation analysis according to the grid node temperature and the boundary conditions of the bottom hole and the wellhead by the KickSimuling function, and updating the calculation result of the flow parameter into a grid node array pNodeNew;
GridGeneration
the GridGeneration function carries out axial segmentation and grid division on a temperature and pressure field solving domain according to a well structure and a drill string structure, and stores axial geometric information of grid nodes and geometric and medium information of a radial heat exchange object;
NormalTPNode
the normalTPnode function analyzes grid flow parameters according to a grid flow parameter calculation method that an inlet and an outlet are both gas-liquid two-phase flow;
GasFrontNode
the GasFrontNode function analyzes grid flow parameters according to a gas phase leading edge grid algorithm and updates a gas phase leading edge grid index and gas phase zone length data;
VerticalCoordinate
calculating the vertical coordinate of the center of the grid and the vertical length of the grid according to the well track; the function requires calling a DirectionParaCal function to calculate the vertical depth according to the well depth and the well inclination angle;
TOriginGeneration
the torrginggeneration function interpolates the vertical temperature distribution data of the formation to generate the original temperature at the vertical depth of the node;
Ini
applying initial conditions to solving intra-domain grid node variables by using an Ini function, and assigning initial values;
HeatResistanceSA
calculating heat exchange resistance between drilling fluid in the drill bit and annular drilling fluid; the function is realized by calling HPipe and HANnu functions to calculate the forced convection heat transfer coefficient of the liquid phase in the drill main and the gas-liquid two-phase flow in the annular space;
HeatResistanceAE
calculating heat exchange resistance between the annular gas-liquid two-phase flow drilling fluid and the stratum; the function is realized by calling a HANnu function to calculate the forced convection heat transfer coefficient of gas-liquid two-phase flow in the annulus and an instantaneous formation heat loss function to calculate the heat transfer thermal resistance from the shaft to the formation;
HPipe
the function of the HPipe function is to calculate the pipe flow forced convection heat transfer coefficient of the drilling fluid in the drill bit and the inner wall of the drill column;
HAnnu
the function of the HANnu function is to calculate the forced convection heat transfer coefficient of the annular flow of the annular gas-liquid two-phase flow drilling fluid and the well wall;
PPTP
calculating a total integration function of thermodynamic properties of the medium at a given temperature and pressure according to the type of the medium, wherein crude oil, water and natural gas are required to be used according to specific types of the medium; calling a thermophysical property calculation function of a corresponding medium to complete the function of the medium;
DPDZ_TwoPhase
calculating the flow parameters and pressure drop gradient of the gas-liquid two-phase flow according to the gas-liquid two-phase flow mechanism model, wherein the flow parameters of the gas-liquid two-phase flow comprise: gas phase velocity, liquid phase velocity, cross section gas content, gas phase density and flow pattern;
TInDrillstem
the TInDrillstem function iteratively calculates the node temperature of the drilling fluid in the drill bit according to a single liquid phase energy conservation equation;
TAnnulus
the TANNILUUS function iteratively calculates the node temperature of the annular drilling fluid according to the energy conservation equation of the gas-liquid two-phase flow;
calculating the constant pressure specific heat of the natural gas by using CpGas;
CpOil, calculating the constant pressure specific heat of the crude oil;
CpWater, calculating the constant pressure specific heat of formation water;
DenGas, calculating the natural gas density;
DenOilSatured, calculating the saturated crude oil density;
DenOilUnsated, calculating the density of the unsaturated crude oil;
RsStanding, calculating the dissolved gas-oil ratio of the crude oil;
surface tension, calculating the oil-water interfacial tension;
thcon gas, calculating the natural gas thermal conductivity;
thcon oil, calculating crude oil thermal conductivity;
ThconWater, calculating formation water thermal conductivity;
ViscGas, calculating natural gas viscosity;
ViscOil, calculate crude oil viscosity;
ViscWater, calculate the formation water viscosity.
3. The method for numerical inversion analysis while drilling of the gas reservoir state parameters as recited in claim 2, wherein: setting boundary conditions and initial conditions of the gas reservoir state parameter while-drilling numerical inversion interpretation method:
boundary condition
At the drill string inlet:
Qm＝C (45-1)
T1(z＝0)＝TE (45-2)
wherein T isEIs the drilling fluid inlet temperature, QmThe drilling fluid discharge capacity;
② the annular outlet:
p＝patm (46)
(iii) the junction of the drill string at the bottom of the well and the annulus, the temperature of the drilling fluid is the same, there are
Initial conditions
The initial condition of the temperature field is the calculated shaft temperature under the condition of single liquid phase,
namely:
Ti(z,0)＝Ti,st(z) (48)
the initial condition of the pressure field is the calculated shaft pressure under the condition of single liquid phase,
namely:
pi(z,0)＝pi,st(z) (49)
and thirdly, when gas invasion starts, no gas phase exists in each mesh of the annulus:
αi(z,0)＝0 (50)。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810238082.1A CN108509703B (en) | 2018-03-22 | 2018-03-22 | Gas reservoir state parameter while-drilling numerical inversion analysis method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810238082.1A CN108509703B (en) | 2018-03-22 | 2018-03-22 | Gas reservoir state parameter while-drilling numerical inversion analysis method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108509703A CN108509703A (en) | 2018-09-07 |
| CN108509703B true CN108509703B (en) | 2022-01-28 |
Family
ID=63378041
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810238082.1A Expired - Fee Related CN108509703B (en) | 2018-03-22 | 2018-03-22 | Gas reservoir state parameter while-drilling numerical inversion analysis method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108509703B (en) |
Families Citing this family (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110424954B (en) * | 2019-08-09 | 2020-05-05 | 西南石油大学 | Annular transient water hammer model based on wellbore annular transient multiphase flow flowing characteristics |
| CN112818503B (en) * | 2019-11-15 | 2022-11-01 | 中国石油天然气股份有限公司 | Method for determining moving speed of temporary plugging ball |
| CN111058796B (en) * | 2019-11-25 | 2021-11-09 | 西南石油大学 | Method for improving shale gas well oil layer casing well cementation quality |
| CN111415031B (en) * | 2020-02-19 | 2023-05-09 | 中石油煤层气有限责任公司 | Method for predicting productivity of coal-bed gas well |
| CN111523248B (en) * | 2020-05-12 | 2024-05-28 | 国电新能源技术研究院有限公司 | Modeling method for dynamic mechanism model of coal-fired power plant |
| CN112347675B (en) * | 2020-10-13 | 2022-08-19 | 中国石油大学(华东) | Method for cooperatively regulating and controlling reservoir natural gas hydrate phase state by drilling fluid additive and temperature and pressure field |
| CN113006769B (en) * | 2021-03-17 | 2022-07-26 | 中国石油大学(华东) | Intelligent well killing method and device for complex pressure system stratum |
| CN114482923B (en) * | 2022-01-25 | 2023-06-30 | 中国石油大学(华东) | Drilling fluid circulation heat exchange control method and system considering material phase change |
| CN116341423B (en) * | 2023-05-30 | 2023-08-01 | 西南石油大学 | Calculation method of oil-water two-phase flow sliding speed model |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CA2803315A1 (en) * | 2010-07-29 | 2012-02-02 | Exxonmobil Upstream Research Company | Methods and systems for machine-learning based simulation of flow |
| CN102943620A (en) * | 2012-08-27 | 2013-02-27 | 中国石油大学(华东) | Pressure-controlled drilling method based on drilling annulus wellbore multi-phase flow computing |
| CN103226641A (en) * | 2013-05-10 | 2013-07-31 | 中国石油大学(华东) | Coupling calculation method of deepwater gas-liquid two-phase flow circulating temperature and pressure |
| CN104153759A (en) * | 2014-07-30 | 2014-11-19 | 中国石油集团钻井工程技术研究院 | Gas-liquid two-phase flow simulating and calculating method for pressure-control well drilling |
| CN107762498A (en) * | 2017-09-27 | 2018-03-06 | 中国地质调查局油气资源调查中心 | A kind of pressure analysis method in the area of tight gas reservoir straight well volume fracturing two |
-
2018
- 2018-03-22 CN CN201810238082.1A patent/CN108509703B/en not_active Expired - Fee Related
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CA2803315A1 (en) * | 2010-07-29 | 2012-02-02 | Exxonmobil Upstream Research Company | Methods and systems for machine-learning based simulation of flow |
| CN102943620A (en) * | 2012-08-27 | 2013-02-27 | 中国石油大学(华东) | Pressure-controlled drilling method based on drilling annulus wellbore multi-phase flow computing |
| CN103226641A (en) * | 2013-05-10 | 2013-07-31 | 中国石油大学(华东) | Coupling calculation method of deepwater gas-liquid two-phase flow circulating temperature and pressure |
| CN104153759A (en) * | 2014-07-30 | 2014-11-19 | 中国石油集团钻井工程技术研究院 | Gas-liquid two-phase flow simulating and calculating method for pressure-control well drilling |
| CN107762498A (en) * | 2017-09-27 | 2018-03-06 | 中国地质调查局油气资源调查中心 | A kind of pressure analysis method in the area of tight gas reservoir straight well volume fracturing two |
Non-Patent Citations (5)
| Title |
|---|
| Coupled modeling circulating temperature and pressure of gas–liquid two phase flow in deep water wells;xuncheng song, zhichuan guan;《Journal of Petroleum Science and Engineering》;20120615;第124-131页 * |
| 多相流全瞬态温度压力场耦合模型求解及分析;何淼等;《石油钻探技术》;20150325(第02期);第25-32页 * |
| 庄茁 ; 柳占立 ; 王永亮.页岩油气高效开发中的基础理论与关键力学问题.《力学季刊》.2015, * |
| 气液两相流循环温度和压力预测耦合模型;宋洵成等;《石油钻采工艺》;20121120(第06期);第5-9页 * |
| 田冷 ; 肖聪 ; 刘明进 ; 顾岱鸿.体积压裂水平井的页岩气产能预测新方法.《大庆石油地质与开发》.2015, * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108509703A (en) | 2018-09-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108509703B (en) | Gas reservoir state parameter while-drilling numerical inversion analysis method | |
| RU2531696C2 (en) | Device and method for simulation of well structure and operating performances | |
| Yuan et al. | A novel integrated workflow for evaluation, optimization, and production predication in shale plays | |
| Yang et al. | Prediction of wellbore and formation temperatures during circulation and shut-in stages under kick conditions | |
| Behmanesh et al. | Treatment of rate-transient analysis during boundary-dominated flow | |
| CN113006769B (en) | Intelligent well killing method and device for complex pressure system stratum | |
| Modaresghazani et al. | Investigation of the relative permeabilities in two-phase flow of heavy oil/water and three-phase flow of heavy oil/water/gas systems | |
| Zhang et al. | Efficient flow rate profiling for multiphase flow in horizontal wells using downhole temperature measurement | |
| Da Silva et al. | A review of coupled dynamic well-reservoir simulation | |
| Shi et al. | Production forecasting of gas condensate well considering fluid phase behavior in the reservoir and wellbore | |
| Kortukov et al. | Fiber optic measurements as real time PLT with new transient interpretation | |
| Velasco et al. | New production performance and prediction tool for unconventional reservoirs | |
| Hasan et al. | Modeling changing storage during a shut-in test | |
| US4508169A (en) | Method for determining connate water saturation and salinity in reservoirs | |
| Takbiri-Borujeni et al. | Upscaling the steam-assisted-gravity-drainage model for heterogeneous reservoirs | |
| Zhang et al. | The optimal model of water control completion based on source function and network model | |
| Mahmud et al. | Investigating relative permeability measurements using unsteady-state core flooding method | |
| Ismail et al. | Determination of the best-fit multiphase flow correlation for high water-cut wells using prosper | |
| Wang et al. | Development of a multiphase variable mass flow model and its application in well control for variable gradient managed pressure drilling | |
| Swanson et al. | Experimental measurement and modeling of gas solubility in invert emulsion drilling fluids explains surface observations during kicks | |
| Taghavi et al. | The Impact of Autonomous Inflow Control Valve on Improved Oil Recovery in a Thin-Oil-Rim Reservoir | |
| Queipo et al. | A model for the integrated optimization of oil production systems | |
| Tokita et al. | Development of the MULFEWS multi-feed wellbore simulator | |
| Ali et al. | A novel interpretation approach for production logging | |
| Miao et al. | A Dynamic Plunger Lift Model for Shale Gas Wells |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| CB03 | Change of inventor or designer information |
Inventor after: Song Xuancheng Inventor after: Wang Haoyan Inventor after: Peng Jie Inventor after: Liu Yongwang Inventor after: Guan Zhichuan Inventor before: Song Xuancheng Inventor before: Wang Haoyan Inventor before: Peng Jie Inventor before: Liu Yongwang Inventor before: Guan Zhichuan |
|
| CB03 | Change of inventor or designer information | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20220128 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |