阅读说明：本技术 非均质储层水平井砂液产出剖面动态演化模拟与预测方法 (Heterogeneous reservoir horizontal well sand liquid production profile dynamic evolution simulation and prediction method ) 是由 董长银 王力智 甘凌云 皇凡生 王卫阳 陈德春 周博 宋雅君 陈琛 于 2021-10-08 设计创作，主要内容包括：本发明公开了非均质储层水平井砂液产出剖面动态演化模拟与预测方法,涉及油气开采工程技术领域,利用本发明提供的方法,可以简易快速模拟砂液产出剖面演化过程,预测给定生产条件和生产时间下的地层流体和地层砂产出剖面,为水平井控砂和控水管柱设计提供重点产出部位,即重点控制井段和位置,提供控砂控水的目的性和综合效果。利用本发明提供的方法,可以快速识别非均质储层水平井的高速入流位置,判断防砂控水完井筛管的高冲蚀发生位置,指导优化筛管结构,在高冲蚀风险位置使用高抗冲蚀的筛管,提高完井筛管的综合抗冲蚀能力。另外,利用本发明提供的局部入流系数可以快速计算实际局部高速入流位置的最高流速,用于冲蚀损坏速率预测及预防。(The invention discloses a method for simulating and predicting the dynamic evolution of a sand liquid production profile of a horizontal well of a heterogeneous reservoir, which relates to the technical field of oil and gas exploitation engineering. By utilizing the method provided by the invention, the high-speed inflow position of the horizontal well of the heterogeneous reservoir can be quickly identified, the high erosion occurrence position of the sand-control water-control well completion sieve tube is judged, the optimization of the sieve tube structure is guided, the sieve tube with high erosion resistance is used at the high erosion risk position, and the comprehensive erosion resistance of the well completion sieve tube is improved. In addition, the local inflow coefficient provided by the invention can be used for quickly calculating the highest flow speed of the actual local high-speed inflow position for erosion damage rate prediction and prevention.)

1. The method for simulating and predicting the dynamic evolution of the sand-liquid production profile of the horizontal well of the heterogeneous reservoir is characterized by comprising the following steps of:

s1: according to the well body structure and the length of the horizontal production section of the horizontal well, carrying out grid division on the horizontal production section of the horizontal well:

carrying out sectional grid division on the long production well section of the horizontal well at intervals of delta H, wherein the horizontal well is divided into N sections in total, the serial number of any one section is represented by a serial number i, and i is more than or equal to 1 and less than or equal to N;

calculating to obtain a heterogeneous distribution profile of the cohesive strength, the porosity and the permeability of the rock according to the initial acoustic time difference, the density and neutron logging data:

the formula for calculating the cohesive strength of the reservoir rock according to the acoustic time difference and the density logging data is shown as the formula (1):

(1)

in the formula (I), the compound is shown in the specification,the cohesive strength of the rock at the jth logging depth is MPa; sigma_cThe uniaxial compressive strength of the rock is MPa; rho_rIs the density of stratum rock in kg/m³；∆t_hIs the transverse wave time difference, mu s/m; Δ t_vIs the longitudinal wave time difference, mu s/m;

the formula for calculating the neutron porosity of the reservoir according to the neutron logging data is shown in the formulas (2) and (3):

(2)

(3)

in the formula (I), the compound is shown in the specification,is neutron porosity, dimensionless;is the density of the rock skeleton in g/cm³;Is the density of the water of the stratum in g/cm³;The mud content is dimensionless;is the density of mud in g/cm³;The initial porosity is obtained by calculation, and the method is dimensionless;density porosity, dimensionless;

the calculation formula for calculating the reservoir permeability is shown in (4):

(4)

wherein k is permeability, mD; s_wiDimensionless for irreducible water saturation;

and (3) performing overall correction on porosity and permeability according to the physical property data of the production horizon test:

defining the ratio of the porosity correction coefficient and the permeability correction coefficient to the average porosity and permeability of the reservoir calibrated by reservoir oil deposit engineering and the average porosity and permeability of the non-uniform profile calculated according to logging data, and calculating by using the following steps (5) and (6):

(5)

(6)

in the formulaThe porosity correction coefficient and the permeability correction coefficient are respectively dimensionless;the average porosity and permeability of the reservoir are respectively calibrated for reservoir oil reservoir engineering, and the units are respectively decimal (dimensionless) and mD;respectively calculating the initial porosity and permeability of the jth logging depth point obtained according to the logging information and the logging information (3) and (4), wherein the units are decimal dimensionless and mD; m is the logging data quantity and is dimensionless;

correcting non-uniform distribution data of porosity and permeability obtained according to logging data according to reservoir average porosity and permeability calibrated by reservoir engineering, as shown in (7) and (8):

(7)

(8)

in the formula (I), the compound is shown in the specification,to be corrected toInitial porosity for j log depths, dimensionless;initial permeability, mD, for corrected jth log depth;

definition ofIs the average porosity of the i well section, is dimensionless, and is calculated by the method that the i well section rangesAn arithmetic mean;

S2：

fitting and calculating an average oil-production-simulated index, a gas-production-simulated index and a water-production-simulated index according to the yield, the permeability, the length of the production section and the production pressure difference of the horizontal well; calculating the flow pressure drop of the horizontal section and the wellbore pressure of each grid section according to the inflow intensity distribution of the horizontal section, the static pressure of the stratum and the bottom hole flow pressure of the toe end;

the average pseudo oil recovery/water/gas index is calculated using equations (9), (10), (11), respectively:

(9)

(10)

(11)

in the formula (I), the compound is shown in the specification,in order to average the pseudo-oil-production index,；m is the oil production of the horizontal well³K0 is the reservoir average permeability, mD, calibrated by reservoir oil deposit engineering; l is_hM is the length of a horizontal well production section;production differential pressure, MPa;

J_XWis an average pseudo-water sampling index, m²/（N·s）；Q_wM is the water yield of the horizontal well³/d;

J_xgIs an average pseudo-gas production index, m²/（N·s）；Q_gFor gas production of horizontal wells, m³/d；

Initializing the inflow intensity of the shaft according to the uniform inflow profile of the whole well, and using the intensity to calculate the pressure distribution of the shaft, wherein the calculation formula of the actual bottom hole flow pressure of each well section of the horizontal section is shown as the formula (12):

(12)

in the formula, P_wf(i)The actual bottom hole flowing pressure of the ith well section is MPa, when i =1,to calibrate the bottom hole flow pressure;the flow pressure drop per unit length of the ith well section is MPa/m;calculating according to a flow friction resistance pressure drop formula;

S3：

subtracting the actual wellbore pressure of each well section by using the formation static pressure to obtain the reservoir inflow production pressure difference of each grid section; calculating the inflow intensity of oil, gas and water of the grid section according to the oil-production-simulated index, the gas-production-simulated index and the water-production-simulated index to form an initial inflow profile:

the calculation formula of the actual production pressure difference of each well section of the horizontal section is shown as the formula (13):

(13)

in the formula (I), the compound is shown in the specification,-actual differential production pressure, MPa, at point i; p_r-formation static pressure, MPa;

the volume calculation formula of the reservoir rock framework of each well section is shown as the formula (14):

(14)

in the formula, TV_(i)Is the volume of the rock skeleton of the reservoir at the ith well section, m³；R_wIs the wellbore radius, m; r_mThe formation sand production radius, m;

the fluid inflow intensity calculation formula is shown in formulas (15), (16), (17) and (18):

(15)

(16)

(17)

(18)

in the formula, Q_o(i)Oil inflow intensity at i-section, m³/(d·m）；Q_w（i）Water inflow intensity at section i, m³/（d·m）；Q_g（i）Intensity of gas inflow at section i, m³/（d·m）；Q_l(i)Intensity of fluid inflow at segment i, m³/（d·m）；k_（i）All k in i well section_0jAverage value of (d);

the total flow rate calculation is shown in equation (19):

(19)

V_f(i)the fluid inflow velocity is m/min; d_wIs the wellbore diameter, m;

S4：

calculating the sand production rate and the sand concentration according to the inflow velocity of the fluid of each grid section and the cohesive strength of the rock to form a sand production strength non-uniform distribution profile; calculating a comprehensive inflow coefficient to obtain a comprehensive inflow non-uniform section:

calculating V at each section of the whole wellbore section_f(i)S with each log depth_0jAnd then calculating the average value through the whole well section to obtain the average fluid inflow velocity V_faAnd average cohesive strength S_0aSaid V is_fa、S_0aEach of V of the whole well section_f(i)And S_0jThe arithmetic mean of (a):

the calculation formula of the simulated sand index is shown as the formula (20):

(20)

in the formula, J_xsTo draw up a sand index, m/s; l is_qsThe sand production strength is L/(d.m); u shape_fFluid viscosity, mpa.s;

calculating the sand production rate and the sand concentration according to the inflow conditions, as shown in the formula (21):

(21)

Q_s（i）the sand production rate at the well section i is L/(d.m); s_0（i）Average cohesive strength of the ith interval, total S in the ith interval_0jAverage value of (d);

the sand content calculation formula is shown as formula (22):

(22)

C_s（i）the sand concentration at the well section i is dimensionless;

the integrated inflow coefficient calculation formula is shown in equation (23):

(23)

in the formula, Z_（i）The comprehensive inflow coefficient is dimensionless; q_saThe average sand production rate for the whole interval, L/(d.m), is Q for each interval_s（i）Calculating an average value;

s5: drawing a non-uniform distribution profile diagram of indexes such as horizontal section reservoir rock cohesive strength, porosity, permeability, oil-gas-water inflow strength, sand production strength, liquid sand content, comprehensive inflow index and the like in the production state at the initial time (t = 0): drawing a non-uniform distribution profile broken line graph by taking the well depth of the horizontal well as a horizontal coordinate and taking indexes such as the cohesive strength, porosity, permeability, oil-gas-water inflow strength, sand production strength, liquid sand content, comprehensive inflow index and the like of reservoir rock of each well section of the horizontal well as a vertical coordinate;

s6: setting the production time step Δ t, calculating the inflow indexes and distribution thereof on each grid section after the time t on the basis of the non-uniform inflow profile of each index at the initial time (t = 0), and calculating the porosity change of the reservoir and the porosity, permeability and cohesive strength after the change according to the sand output of each grid section of the well section:

the sand output of each grid section in the time t period is calculated firstly, the accumulated sand output of the i well section grid in the time t is calculated according to the formula (24):

(24)

in the formula,. DELTA.V_s（i）The i well section accumulated sand yield m within t time³；

Calculating the porosity at the end of the time when the sand is discharged according to the sand discharge amount, wherein the calculation formula is shown as the formula (25):

(25)

in the formula (I), the compound is shown in the specification,the porosity of the near-well stratum at the end of the time t of the i well section is dimensionless;the porosity of the near well stratum at the time of the i well section Δ t is dimensionless; in the first iteration of the process,equal to that calculated in step S1；

The calculation formula of the parameters of the porosity, permeability and cohesive strength variation amplitude is shown in the formulas (26), (27) and (28):

(26)

(27)

(28)

the i well section porosity change rate amplitude parameter at the end of t time is dimensionless;the i well section permeability change rate amplitude parameter at the end of the t time is dimensionless;the cohesive strength change rate amplitude parameter of the i well section at the end of the t time is dimensionless;

and (3) calculating the permeability change, and according to the conversion of the porosity change ratio, calculating the formula as shown in formula (29):

(29)

in the formula, k_（i,t1）The near-well stratum permeability at the end of the time at the i well section, mD; k is a radical of_（i,t0）For the initial near-well permeability at the time of Δ t for the i-well section, mD, in the first iteration, k_（i,t0）K equal to that calculated in step S3_（i）；

Calculating the change of cohesive strength as shown in formula (30):

(30)

in the formula, S_{0（i，t1）}The cohesive strength of the near-well stratum at the end of the time of i well section is MPa; s_{0（i，t0）}At the beginning of the i well sectionWell formation cohesive strength, MPa, in the first iteration, S_{0（i，t0）}Is equal to S calculated in step S4_0（i）；

S7: recalculating pressure distribution of each grid section of a well bore of the production section according to an inflow profile at the initial time (t = 0), and subtracting the pressure of the well bore by using formation static pressure to obtain reservoir inflow production differential pressure of each grid section; calculating the inflow intensity of oil, gas and water of a new grid section according to the porosity and permeability at the time t to form a new inflow profile;

calculating the pressure distribution of the shaft by adopting the inflow intensity of the shaft of each shaft section at the initial moment (t = 0), wherein the calculation formula of the actual bottom hole flow pressure of each point of the horizontal section is shown as the formula (31):

(31)

in the formula, P_{wf（i，t1）}The actual bottom hole flowing pressure of the ith well section at the end of delta t is MPa; p_{wf（i-1，t1）}The actual bottom hole flow pressure of the i-1 th well section at the end of delta t is Mpa; delta P_{k（i，t1）}The flow pressure loss of the ith well section in unit length at the end of delta t is MPa/m; delta P_{k（i，t1）}The flow rate is calculated by adopting a flow friction pressure drop formula according to the inflow rate at the initial time (t = 0), and the actual wellbore flow rate of the ith wellbore section is the sum of the inflow rates of all grids from the wellbore section to the toe end;

the calculation formula of the actual production pressure difference of each well section of the horizontal section is shown as the formula (32):

(32)

in the formula,. DELTA.P_（i，t1）The actual production pressure difference of the ith well section at the end of delta t is MPa;

the fluid output speed calculation formula at the time of Δ t is shown as the formulas (33), (34), (35) and (36):

in the formula, Q_o(i,t1)The oil inflow intensity at i well section at the end of time delta t, m³/（d·m）；Q_w(i,t1)Water inflow intensity at i well section at end of time Δ t, m³/（d·m）；Q_g(i,t1)The gas inflow intensity at i well section at the end of time Δ t, m³/（d·m）；Q_l(i,t1)The fluid inflow intensity at i well section at the end of time at, m³/（d·m）；k_(i,t1)The average value of each permeability at the i well section at the end of the delta t moment;

the total flow rate calculation is shown in equation (37):

(37)

V_f(i,t1)the fluid inflow velocity at the i well section at the end of the time delta t is m/min; da is the diameter of a shaft at the well section i, and m;

s8: calculating the sand output, sand content and comprehensive inflow index of each new grid section according to the new inflow strength and cohesion strength at the time t, and obtaining a new non-uniform distribution profile of all inflow indexes at the time t:

calculating the average inflow velocity and cohesive strength of the whole well section to obtain the average inflow velocity V of the fluid at the end of the time delta t_fa(t1)And average cohesive strength S_0a(t1)，V_fa(t1)And S_0a(t1)V of all well sections respectively_f(i,t1)And S_0(i,t1)The arithmetic mean of (a);

the sand production index calculation formula is shown as the formula (38):

(38)

in the formula, J_xs(t1)At time t, the sand index_qsThe sand production strength is L/(d.m); u shape_fFluid viscosity, mpa.s;

calculating the sand production rate and the sand concentration according to the inflow conditions, as shown in the formula (39):

(39)

Q_s(i,t1)the sand production rate at the i well section at the time t is L/(d.m);

the sand content calculation formula is shown as the formula (40):

(40)

C_s(i,t1)the sand concentration at the end i well section at the time t is dimensionless;

the integrated inflow coefficient calculation formula is shown as equation (41):

(41)

in the formula, Z_(i,t1)The comprehensive inflow coefficient at the i well section at the end of the delta t time is dimensionless; q_sa(t1)Average sand production rate of whole well section at delta t time, L/(d.m), Q for whole well section_s(i,t1)Taking an average value;

s9: calculating local inflow coefficients

The calculation formula is shown in equations (42) and (43):

in the formula, PIF_sThe local inflow coefficient of the sand production strength is dimensionless; PIF_lThe local inflow coefficient is the liquid production intensity, and is dimensionless; q_smax(t1)The highest sand production rate of the whole well section at the end of time delta t, L/(d.m), V_fmax(t1)The maximum fluid inflow velocity of the whole well section at the end of the time delta t, m/min, and Q at each point of the whole well_s(i,t1)And V_f(i,t1)Taking the highest value;

s10: the production time is increased by Δ t, the iterative calculation of the next time step Δ t is continued from S6, the dynamic evolution of the sand liquid output profile is circularly realized, and various inflow profile predictions at any time t are realized.

2. The heterogeneous reservoir horizontal well sand fluid production profile dynamic evolution simulation and prediction method as claimed in claim 1, wherein Δ H is in a range of 0.2-0.5 m.

3. The heterogeneous reservoir horizontal well sand fluid production profile dynamic evolution simulation and prediction method of claim 1,respectively taking 2.65, 1 and 2.42.

Technical Field

The invention relates to the technical field of oil and gas exploitation engineering, in particular to a method for simulating and predicting dynamic evolution of a sand liquid production profile of a horizontal well of a heterogeneous reservoir.

Background

The horizontal well technology is an important technical means for developing petroleum and natural gas, and has the advantages of long production well section (hundreds of meters to thousands of meters), large drainage area, high oil extraction speed and the like. For a horizontal well with loose sandstone and an easily-produced oil and gas reservoir, the exploitation of the horizontal well faces two troublesome problems of sand production and water production. Sand control and water control of long-well-section horizontal wells are one of core technologies for ensuring normal exploitation and production of the long-well-section horizontal wells.

During horizontal well production, it is generally believed or assumed that the feed inflow profile of the reservoir to the horizontal wellbore is uniform. Whether for a high producing gas well (e.g., 200 million square production per day) or an oil well (e.g., 2000 square fluid production per day), the flow rate of fluid to the screens in the wellbore is extremely low, far from an erosive damage condition, if calculated in terms of a horizontal segment uniform inflow. However, a number of field practices have shown that partially sanded wells and horizontal wells exhibit washout damage to the completion screen and that within a limited longer interval, the washout perforation locations are unique, as shown in figure 8. It was concluded that the feed inflow profile of the reservoir to the horizontal well was non-uniform.

Due to the heterogeneity of the cementing strength, permeability and the like of the oil and gas reservoir and the influence of the reservoir on the distribution rule of the porosity and permeability of the reservoir in the sand production process, the change is caused; the change of the porosity and permeability distribution of the reservoir along the axis direction of the shaft can influence the reservoir circulation and inflow profile, and further reversely influence the sand production rule. Namely, a synergistic influence mechanism exists between the formation sand of the heterogeneous sand-prone reservoir and the fluid production profile. The above synergistic influence mechanism causes the sand-liquid inflow profile of the horizontal well of the heterogeneous reservoir to dynamically evolve along with the production time, so that the heterogeneity of sand-liquid inflow becomes stronger and stronger, as shown in fig. 9.

The sand control and water control of horizontal wells of reservoirs easy to produce sand are key technologies for ensuring high-efficiency production. Because the horizontal well production section is long in length, the fluid production section and the sand production section of the horizontal well long production section have heterogeneity caused by the difference of physical properties of reservoirs, and the heterogeneity dynamically evolves along with the production time. The accurate depiction of the sand liquid output section can make sand control and water control operation have a target, and improve the sand control and water control effect. At present, reservoir physical properties are mainly recognized by taking static data measured at the initial development stage as a main part, and a simulation and prediction method is lacked for the dynamic evolution of a sand-liquid production profile of a horizontal well of a heterogeneous reservoir, so that the optimization design of sand and water control measures of the heterogeneous reservoir is lacked with a dynamic basis.

Disclosure of Invention

Aiming at the technical problems, the invention provides a method for simulating and predicting the dynamic evolution of the sand-liquid production section of the horizontal well with the heterogeneous reservoir, and realizes the simulation of the dynamic evolution process of the sand-liquid production section of the horizontal well with the production time with the production section of the heterogeneous reservoir and the prediction of the production section. Finally, the judgment and identification of the horizontal well sand production and water production profile are promoted from the initial static state to the dynamic state, so that effective support is provided for solving the problem of sand and water cooperative control of the long-well-section horizontal well, and the method has general reference value and scientific significance for the dynamic optimization of the horizontal well long-term development mode and the real-time optimization of the single-well production system.

Interpretation of professional terms: the inflow profile comprises a comprehensive inflow coefficient profile, an inflow velocity profile, a sand production intensity profile, a comprehensive inflow coefficient profile, an inflow velocity profile and a sand production intensity profile which are respectively the set of the comprehensive inflow coefficient, the inflow velocity and the sand production intensity of the whole well section at the moment.

The technical scheme for solving the technical problem of the invention is as follows:

the method for simulating and predicting the dynamic evolution of the sand-liquid production profile of the horizontal well of the heterogeneous reservoir comprises the following steps:

s1: carrying out grid division on the horizontal production section of the horizontal well according to the well body structure of the horizontal well and the length of the horizontal production section; calculating to obtain a heterogeneous distribution profile of the cohesive strength, the porosity and the permeability of the reservoir rock according to the initial acoustic time difference, the density and the neutron logging information; performing overall correction on porosity and permeability according to the physical property data of the production horizon test;

in the step S1, carrying out sectional grid division on the horizontal well long production well section by taking delta H (0.2-0.5 m is recommended) as an interval, dividing the horizontal well long production well section into N sections in total, using a serial number i to represent the serial number of any one section, and using i to be more than or equal to 1 and less than or equal to N, as shown in figure 1;

the formula for calculating the cohesive strength of the reservoir rock according to the acoustic time difference and the density logging data is shown as the formula (1):

(1)

in the formula (I), the compound is shown in the specification,the cohesive strength of the rock at the jth logging depth is MPa; sigma_cThe uniaxial compressive strength of the rock is MPa; rho_rIs the density of stratum rock in kg/m³；∆t_hIs the transverse wave time difference, mu s/m; Δ t_vIs the longitudinal wave time difference, mu s/m;

the formula for calculating the porosity of the reservoir according to the neutron logging data is shown in the formulas (2) and (3):

(2)

(3)

in the formula (I), the compound is shown in the specification,is neutron porosity, dimensionless;is the density of the rock skeleton in g/cm³;Is the density of the water of the stratum in g/cm³;The mud content is dimensionless;is the density of mud in g/cm³;The initial porosity is obtained by calculation, and the method is dimensionless;density porosity, dimensionless; general sandstone reservoirRespectively taking 2.65, 1 and 2.42.

The calculation formula for calculating the reservoir permeability is shown in (4):

(4)

wherein k is permeability, mD; s_wiTo restrict water saturation, noneDimension.

Defining the porosity correction coefficient and the permeability correction coefficient as follows: the ratio of the average porosity and permeability of the reservoir calibrated by reservoir oil reservoir engineering to the average porosity and permeability of the non-uniform profile calculated according to the logging data is calculated by the following steps (5) and (6):

(5)

(6)

in the formulaThe porosity correction coefficient and the permeability correction coefficient are respectively dimensionless;the average porosity and permeability of the reservoir are respectively calibrated for reservoir oil reservoir engineering, and the units are respectively decimal (dimensionless) and mD;respectively calculating the initial porosity and permeability of the jth logging depth according to the logging information and the formulas (3) and (4), wherein the units are decimal (dimensionless) and mD; and M is the logging data quantity and is dimensionless.

Correcting non-uniform distribution data of porosity and permeability obtained according to logging data according to the average porosity and permeability of a reservoir layer calibrated by oil reservoir engineering, wherein the non-uniform distribution data are shown as formulas (7) and (8):

(7)

(8)

in the formula，The corrected initial porosity of the jth logging depth is dimensionless;initial permeability, mD, for corrected jth log depth;

definition ofIs the average porosity of the i well section, is dimensionless, and is calculated by the method that the i well section rangesIs calculated as the arithmetic mean of (1).

S2: fitting and calculating an average oil-production-simulated index, a gas-production-simulated index and a water-production-simulated index according to the yield, the permeability, the length of the production section and the production pressure difference of the horizontal well; calculating inflow intensity distribution of the horizontal section according to the uniform inflow situation of the production section; calculating the flow pressure drop of the horizontal section and the wellbore pressure of each grid section according to the inflow intensity distribution of the horizontal section, the static pressure of the stratum and the bottom hole flow pressure of the toe end;

preferably, in step S2, the pseudo oil recovery/water/gas index is defined and calculated as shown in formulas (9), (10), (11), respectively:

(9)

(10)

(11)

in the formula (I), the compound is shown in the specification,in order to average the pseudo-oil-production index,；m is the oil production of the horizontal well³/d;k₀Reservoir average permeability, mD, calibrated for reservoir engineering; l is_hM is the length of a horizontal well production section;production differential pressure, MPa;

J_XWis an average pseudo-water sampling index, m²/（N·s）；Q_wM is the water yield of the horizontal well³/d;

J_xgIs an average pseudo-gas production index, m²/（N·s）；Q_gFor gas production of horizontal wells, m³/d。

Initializing the inflow intensity of the shaft according to the uniform inflow profile of the whole well, and using the intensity to calculate the pressure distribution of the shaft, wherein the calculation formula of the actual bottom hole flow pressure of each well section of the horizontal section is shown as the formula (12):

(12)

in the formula, P_wf(i)The actual bottom hole flowing pressure of the ith well section is MPa, when i =1,to calibrate the bottom hole flow pressure;the flow pressure loss per unit length of the ith well section is MPa/m;and (4) calculating according to a flow friction pressure drop formula to obtain the actual wellbore flow of the ith wellbore section, wherein the actual wellbore flow of the ith wellbore section is the sum of all grid inflow flows from the wellbore section to the toe end.

S3: subtracting the wellbore pressure by using the formation static pressure to obtain the reservoir inflow production pressure difference of each grid section; calculating the inflow intensity of oil, gas and water of the grid section according to the oil-production-simulated index, the gas-production-simulated index and the water-production-simulated index to form an initial inflow profile;

preferably, in step S3, the calculation formula of the actual production pressure difference of each well section of the horizontal section is as shown in formula (13):

(13)

in the formula, P_r-formation static pressure, MPa;actual differential pressure in production, MPa, at point i.

The volume calculation formula of the reservoir rock framework of each well section is shown as the formula (14):

(14)

in the formula, TV_(i)Is the volume of the rock skeleton of the reservoir at the ith well section, m³；R_wIs the wellbore radius, m; r_mThe formation sand radius, m.

The fluid inflow intensity calculation formula is shown in formulas (15), (16), (17) and (18):

(15)

(16)

(17)

(18)

in the formula, Q_o(i)Oil inflow intensity at i-section, m³/(d·m）；Q_w（i）Water inflow intensity at section i, m³/（d·m）；Q_g（i）Intensity of gas inflow at section i, m³/（d·m）；Q_l(i)Intensity of fluid inflow at segment i, m³/（d·m）；k_（i）All k in i well section_0jAverage value of (d);

the total flow rate calculation is shown in equation (19):

(19)

V_f(i)the fluid inflow velocity is m/min; d_wIs the wellbore diameter, m.

S4: calculating the sand production rate and the sand concentration according to the inflow velocity of the fluid of each grid section and the cohesive strength of the rock to form a sand production strength non-uniform distribution profile; calculating a comprehensive inflow coefficient to obtain a comprehensive inflow non-uniform section;

preferably, in step S4, V at each section of the whole well section is calculated_f(i)S with each log depth_0jAnd then calculating the average value through the whole well section to obtain the average fluid inflow velocity V_faAnd average cohesive strength S_0aSaid V is_fa、S_0aEach of V of the whole well section_f(i)And S_0jIs calculated as the arithmetic mean of (1).

The calculation formula of the simulated sand index is shown as the formula (20):

(20)

in the formula, J_xsTo draw up a sand index, m/s; l is_qsThe sand production strength is L/(d.m); u shape_fFluid viscosity, mpa.s.

Calculating the sand production rate and the sand concentration according to the inflow conditions, as shown in the formula (21):

(21)

Q_s（i）the sand production rate at the well section i is L/(d.m); s_0（i）Average cohesive strength of the ith interval, total S in the ith interval_0jAverage value of (d);

the sand content calculation formula is shown as formula (22):

(22)

C_s（i）the sand concentration at the well section i is dimensionless.

The integrated inflow coefficient calculation formula is shown in equation (23):

(23)

in the formula, Z_（i）The comprehensive inflow coefficient is dimensionless; q_saThe average sand production rate, L/(d.m), for the entire interval is calculated for each intervalThe average value is calculated.

S5: drawing a non-uniform distribution profile diagram of indexes such as horizontal section reservoir rock cohesive strength, porosity, permeability, oil-gas-water inflow strength, sand production strength, liquid sand content, comprehensive inflow index and the like in an initial production state;

preferably, in step S5, the horizontal well depth is taken as the horizontal axis, and indexes such as the cohesive strength, porosity, permeability, inflow strength of oil, gas and water, sand production strength, liquid sand content, and comprehensive inflow index of reservoir rock at each interval of the horizontal well are taken as the vertical axis, and the non-uniform distribution profile line graph is drawn by taking fig. 2 as an example.

S6: setting the production time Δ t. And calculating the inflow indexes and the distribution thereof on each grid section after the time t on the basis of the non-uniform inflow profile of each index at the initial time (t = 0). Calculating the sand yield of each grid section in the time period of the Δ t, and calculating the porosity change of the reservoir and the porosity, permeability and cohesive strength after the change according to the sand yield of each grid section of the well section;

according to the invention, in the step S6, the accumulated sand yield of the i well section grid within the Δ t time is preferably calculated as shown in formula (24):

(24)

in the formula,. DELTA.V_s（i）The i well section accumulated sand yield m within t time³。

Calculating the porosity at the end of the time when the sand is discharged according to the sand discharge amount, wherein the calculation formula is shown as the formula (25):

(25)

in the formula (I), the compound is shown in the specification,the porosity of the near-well stratum at the end of the time t of the i well section is dimensionless;the porosity of the near well stratum at the time of i well section Δ t is dimensionless. In the first iteration of the process,equal to that calculated in step S1In each iteration cycle after the first time, each is at the first timeAt the end of one Δ tThe calculation result of (2);

the calculation formula of the parameters of the porosity, permeability and cohesive strength variation amplitude is shown in the formulas (26), (27) and (28):

(26)

(27)

(28)

the i well section porosity change rate amplitude parameter at the end of t time is dimensionless;the i well section permeability change rate amplitude parameter at the end of the t time is dimensionless;and the cohesive strength change rate amplitude parameter of the i well section at the end of the t time is dimensionless.

And (3) calculating the permeability change, and according to the conversion of the porosity change ratio, calculating the formula as shown in formula (29):

(29)

in the formula, k_（i,t1）The near-well stratum permeability at the end of the time at the i well section, mD; k is a radical of_（i,t0）For the initial near-well permeability at the time of Δ t for the i-well section, mD, in the first iteration, k_（i,t0）K equal to that calculated in step S3_（i）In each iteration cycle after the first time, each is at the first timeEqual to one at end_（i,t1）The calculation result of (2);

calculating the change of cohesive strength as shown in formula (30):

(30)

in the formula, S_{0（i，t1）}The cohesive strength of the near-well stratum at the end of the time of i well section is MPa; s_{0（i，t0）}The cohesive strength of the near-well stratum at the time of the arrival of the i well section, MPa, S in the first iteration_{0（i，t0）}Is equal to S calculated in step S4_0（i）In each iteration cycle after the first time, S at each time t_{0（i，t0）}Equal to S at the end of one_{0（i，t1）}The calculation result of (2);

s7: recalculating pressure distribution of each grid section of a well bore of the production section according to an inflow profile at the initial time (t = 0), and subtracting the pressure of the well bore by using formation static pressure to obtain reservoir inflow production differential pressure of each grid section; calculating the inflow intensity of oil, gas and water of a new grid section according to the porosity and permeability at the time t to form a new inflow profile;

preferably, in step S7, the wellbore pressure distribution is calculated by using the wellbore inflow intensity of each well section at the initial time (t = 0). The calculation formula of the actual bottom hole flowing pressure of each point of the horizontal section is shown as the formula (31)

(31)

In the formula, P_{wf（i，t1）}The actual bottom hole flowing pressure of the ith well section at the end of delta t is MPa; p_{wf（i-1，t1）}The actual bottom hole flow pressure of the i-1 th well section at the end of delta t is Mpa; delta P_{k（i，t1）}The flow pressure loss per unit length of the ith well section at the end of time delta t is MPa/m. Delta P_{k（i，t1）}In order to obtain the flow friction pressure drop formula calculation according to the inflow at the initial time (t = 0),the actual wellbore flow of the ith well section is the sum of all grid inflow rates from the well section to the toe end;

the calculation formula of the actual production pressure difference of each well section of the horizontal section is shown as the formula (32):

(32)

in the formula,. DELTA.P_（i，t1）The actual production pressure difference of the ith well section at the end of the time delta t is MPa.

The fluid output speed calculation formula at the time of Δ t is shown as the formulas (33), (34), (35) and (36):

in the formula, Q_o(i,t1)The oil inflow intensity at i well section at the end of time delta t, m³/（d·m）；Q_w(i,t1)Water inflow intensity at i well section at end of time Δ t, m³/（d·m）；Q_g(i,t1)The gas inflow intensity at i well section at the end of time Δ t, m³/（d·m）；Q_l(i,t1)The fluid inflow intensity at i well section at the end of time at, m³/（d·m）；k_(i,t1)The average value of each permeability at the i well section at the end of the time delta t is shown.

The total flow rate calculation is shown in equation (37):

(37)

V_f(i,t1)the fluid inflow velocity at the i well section at the end of the time delta t is m/min; da is the diameter of the well shaft at the well section i and m.

S8: calculating the sand output, sand content and comprehensive inflow index of each new grid section according to the new inflow strength and cohesion strength at the time t to obtain a new non-uniform distribution profile of all inflow indexes at the time t;

preferably, in step S8, the average inflow velocity and the cohesive strength of the whole wellbore section are calculated to obtain the end of the Δ t timeAverage fluid inflow velocity V_fa(t1)And average cohesive strength S_0a(t1)，V_fa(t1)And S_0a(t1)V of all well sections respectively_f(i,t1)And S_0(i,t1)Is calculated as the arithmetic mean of (1).

The sand production index calculation formula is shown as the formula (38):

(38)

in the formula, J_xs(t1)The sand index at the time of t. L is_qsThe sand production strength is L/(d.m); u shape_fFluid viscosity, mpa.s.

Calculating the sand production rate and the sand concentration according to the inflow conditions, as shown in the formula (39):

(39)

Q_s(i,t1)the sand production rate at the i well section at the time t is L/(d.m).

The sand content calculation formula is shown as the formula (40):

(40)

C_s(i,t1)the sand concentration at the i well section at the end of the delta t time is dimensionless.

The integrated inflow coefficient calculation formula is shown as equation (41):

(41)

in the formula, Z_(i,t1)The comprehensive inflow coefficient at the i well section at the end of the delta t time is dimensionless; q_sa(t1)Average sand production rate of whole well section at delta t time, L/(d.m), Q for whole well section_s(i,t1)And (6) taking an average value.

S9: a local Inflow coefficient (PIF coefficient for short) is calculated.

Preferably, in step S9, the PIF coefficient calculation formula is represented by equations (42) and (43):

in the formula, PIF_sThe local inflow coefficient of the sand production strength is dimensionless; PIF_lThe local inflow coefficient is the liquid production intensity, and is dimensionless; q_smax(t1)The highest sand production rate of the whole well section at the end of time delta t, L/(d.m), V_fmax(t1)The maximum fluid inflow velocity of the whole well section at the end of the time delta t, m/min, and Q at each point of the whole well_s(i,t1)And V_f(i,t1)The highest value is taken.

S10: and the production time is increased by Δ t, the iterative calculation of the next time step is continued from S6, the dynamic evolution of the sand liquid production profile is realized circularly, and various inflow profile predictions at any time t are realized.

The invention has the beneficial effects that:

1. the sand control and water control of horizontal wells of reservoirs with easy sand production are key technologies for ensuring high-efficiency production, and due to the lack of a simulation and prediction method for the dynamic evolution of the sand liquid production profile of horizontal wells of heterogeneous reservoirs, the sand control and water control measures of the horizontal wells of heterogeneous reservoirs are lack of purpose. By utilizing the method provided by the invention, the evolution process of the sand liquid production profile can be simply and rapidly simulated, the stratum fluids (oil, gas and water) and the stratum sand production profile under given production conditions and production time can be predicted, and key production positions are provided for the design of the horizontal well sand control and water control tubular columns, namely key control well sections and positions, so that the purposes and comprehensive effects of sand control and water control are provided.

2. By utilizing the method provided by the invention, the high-speed inflow position of the horizontal well of the heterogeneous reservoir can be quickly identified, the high erosion occurrence position of the sand-control water-control well completion sieve tube is judged, the optimization of the sieve tube structure is guided, the sieve tube with high erosion resistance is used at the high erosion risk position, and the comprehensive erosion resistance of the well completion sieve tube is improved. In addition, the local inflow coefficient (PIF coefficient) provided by the invention can be used for quickly calculating the highest flow speed of an actual local high-speed inflow position for erosion damage rate prediction and prevention.

3. The method provided by the invention promotes the judgment and identification of the sand production and water production profile of the horizontal well from static state to dynamic state, and provides effective support for solving the problem of sand-water cooperative control of the horizontal well at a long well section; compared with static sand-water profile identification, the sand control efficiency is greatly improved. And parameter design basis can be provided for the intelligent sand and water control device, more reasonable sand and water outlet intervention can be carried out along with the production change of the stratum, the possible condition of a large amount of sand/water outlet of the oil well in the future can be avoided at the lowest cost, and the remedial operation after sand and water outlet and the excessive sand prevention of a slight sand outlet well section can be avoided. The method has general reference value and scientific significance for dynamic optimization of a long-term development mode of an oil and gas reservoir, real-time optimization of a single-well production system and dynamic adjustment of a recovery efficiency improving mode.

Drawings

FIG. 1 is a schematic diagram of a method for calculating evolution of a dynamic profile produced by heterogeneous sand liquid;

FIG. 2 is a heterogeneous cross-sectional view of the initial porosity of a horizontal well;

FIG. 3 is a schematic diagram of the dynamic evolution distribution of the porosity of a horizontal well;

FIG. 4 is a schematic diagram of the dynamic evolution distribution of permeability of a horizontal well;

FIG. 5 is a schematic diagram of dynamic evolution distribution of a comprehensive inflow coefficient of a horizontal well;

FIG. 6 is a schematic diagram of the dynamic evolution distribution of inflow velocity of a horizontal well;

FIG. 7 is a schematic diagram of the dynamic evolution distribution of sand production intensity of a horizontal well;

FIG. 8 is a schematic diagram of non-uniform inflow of a horizontal well;

FIG. 9 is a schematic diagram of the dynamic evolution of a horizontal well sand fluid production profile of a heterogeneous reservoir.

Detailed Description

In order to clearly explain the technical features of the present invention, the following detailed description of the present invention is provided with reference to the accompanying drawings. The following disclosure provides many different embodiments, or examples, for implementing different features of the invention. To simplify the disclosure of the present invention, the components and arrangements of specific examples are described below. Furthermore, the present invention may repeat reference numerals and/or letters in the various examples. This repetition is for the purpose of simplicity and clarity and does not in itself dictate a relationship between the various embodiments and/or configurations discussed. It should be noted that the components illustrated in the figures are not necessarily drawn to scale. Descriptions of well-known components and processing techniques and procedures are omitted so as to not unnecessarily limit the invention.

The simulation and prediction method for the dynamic evolution of the sand liquid production profile of the horizontal well with the heterogeneous reservoir is applied to the horizontal section of the horizontal well with a certain oil reservoir in the victory oil field, the horizontal section of the horizontal well is 685m long, the water content is 95%, the liquid production amount is about 97t/D, the sand content of the fluid volume is 0.01%, the viscosity of degassed crude oil is 400mPa.s, the average permeability of the stratum is 0.4D, and the porosity is 33%. By using the method, the porosity, the permeability, the comprehensive inflow index, the inflow speed and the sand production intensity of a horizontal segment are predicted and obtained by taking 5 years as the total calculation time and taking 30 days as the single calculation step length, and the dynamic evolution distribution schematic diagrams are shown in fig. 3 to 7.

According to the dynamic evolution prediction result of the output section, stratum particles are peeled off along with the production process, and the porosity and the permeability are integrally and gradually increased and gradually accelerated. The pore permeability change is obvious at 950-1100 m, and the increase range of the porosity and the permeability has obvious positive correlation with the initial physical property. Along with the continuous production process, the inflow speed and the sand production strength are gradually enhanced, areas with high sand production and high permeability are at 950 + 1100 m, 1480 m and 1580 m, and the sand production and the water production have obvious enhancement tendency. And because of the difference between the formation physical properties and the production pressure difference, the overlapping degree of the initial sand production and water production profiles is lower, and the initial sand production and water production profiles gradually tend to overlap along with the production continuation, which shows that the sand production and the water production converge with each other. According to the preliminary simulation result, the inflow speed and the sand output intensity of the high-speed inflow point are increased by 4-6 times, and the sand output intensity of the area is higher than that of other areas by more than 8 times, so that the area is a key sand prevention and erosion prevention object.

Although the embodiments of the present invention have been described with reference to the accompanying drawings, the scope of the present invention is not limited thereto, and various modifications and variations which do not require inventive efforts and which are made by those skilled in the art are within the scope of the present invention.