发明内容

本发明针对上述技术问题，提供了非均质储层水平井砂液产出剖面动态演化模拟与预测方法，实现非均质储层水平井生产段砂液产出剖面随生产时间的动态演化过程模拟和产出剖面预测。最终将对水平井出砂和出水剖面的判断识别由初始静态提升到动态层面，为解决长井段水平井砂水协同控制难题提供有效支撑，对于水平井长期开发方式动态优化、单井生产制度实时优化都具有普遍参考价值和科学意义。

专业术语解释：入流剖面包括“综合入流系数剖面”、“入流速度剖面”、“出砂强度剖面”，“综合入流系数剖面”、“入流速度剖面”、“出砂强度剖面”分别为全井段综合入流系数、入流速度和出砂强度在该时刻的集合。

本发明解决技术问题的技术方案为：

非均质储层水平井砂液产出剖面动态演化模拟与预测方法，包括以下步骤：

S1：根据水平井井身结构和水平生产段长度，对水平井水平生产段进行网格划分；根据初始声波时差、密度和中子测井资料计算得到储层岩石内聚强度、孔隙度、渗透率的非均质分布剖面；根据生产层位测试物性数据进行孔隙度和渗透率总体校正；

所述步骤S1中，对水平井长生产井段以ΔH（建议取0.2~0.5m）为间隔进行分段网格划分，共分为N段，用序号i表示任一分段的序号，1≤i≤N，如图1所示；

根据声波时差、密度测井资料计算储层岩石内聚强度的公式如式（1）所示：

(1)

式中，为第j个测井深度岩石内聚强度，MPa；σ_c为岩石单轴抗压强度，MPa；ρ_r为地层岩石密度，kg/m³；∆t_h为横波时差，μs/m；∆t_v为纵波时差，μs/m；

根据中子测井资料计算储层孔隙度的公式如式（2)、(3）所示：

(2)

(3)

式中，为中子孔隙度，无量纲；为岩石骨架密度，g/cm³;为地层水密度，g/cm³;为泥质含量，无量纲；为泥质密度，g/cm³;为计算得到初始孔隙度，无量纲；为密度孔隙度，无量纲；一般砂岩储层分别取2.65、1、2.42。

计算储层渗透率的计算公式如(4)所示：

(4)

式中，k为渗透率，mD；S_wi为束缚水饱和度，无量纲。

定义孔隙度校正系数、渗透率校正系数为：储层油藏工程标定的储层平均孔隙度、渗透率与根据测井资料计算得到的非均匀剖面的孔隙度、渗透率的平均值的比值，用(5)、(6)计算：

(5)

(6)

式中分别为孔隙度校正系数、渗透率校正系数，无量纲；分别为储层油藏工程标定的储层平均孔隙度和渗透率，单位分别为小数（无量纲）和mD；分别为根据测井资料和式（3）、(4）计算得到的第j个测井深度的初始孔隙度和渗透率，单位分别为小数（无量纲）和mD；M为测井数据数量，无量纲。

根据油藏工程标定的储层平均孔隙度和渗透率，校正根据测井数据得到的孔隙度、渗透率的非均匀分布数据，如式（7)、(8）所示：

(7)

(8)

式中，为校正后第j个测井深度的初始孔隙度，无量纲；为校正后第j个测井深度的初始渗透率，mD；

定义为i井段平均孔隙度，无量纲，其计算方法为i井段范围内全部的算术平均值。

S2：根据水平井产量、渗透率、生产段长度和生产压差，拟合计算平均拟采油指数、拟采气指数和拟采水指数；按照生产段均匀入流情形计算水平段入流强度分布；根据水平段入流强度分布和地层静压及趾端井底流压，计算水平段流动压降及各网格段井筒压力；

根据本发明优选的，所述步骤S2中，定义拟采油/水/气指数分别如式（9)、(10)、(11）所示并用其计算：

(9)

(10)

(11)

式中，为平均拟采油指数，；为水平井产油量，m³/d;k₀为储层油藏工程标定的储层平均渗透率，mD；L_h为水平井生产段长度，m；为生产压差，MPa；

J_XW为平均拟采水指数，m²/（N·s）；Q_w为水平井产水量，m³/d;

J_xg为平均拟采气指数，m²/（N·s）；Q_g为水平井产气量，m³/d。

按照全井均匀入流剖面初始化井筒入流强度，用于计算井筒压力分布，水平段各井段的实际井底流压的计算公式如式（12）所示：

(12)

式中，P_wf(i)为第i井段实际井底流压，MPa，当i=1时，，为标定井底流压；为第i井段单位长度流动压力损失，MPa/m；根据流动摩阻压降公式计算得到，第i井段实际井筒流量为该井段至趾端全部网格入流量之和。

S3：使用地层静压减去井筒压力得到各网格段的储层入流生产压差；根据拟采油指数、拟采气指数和拟采水指数计算网格段的油、气、水入流强度，形成初始入流剖面；

根据本发明优选的，所述步骤S3中，水平段各井段的实际生产压差的计算公式如式（13）所示：

(13)

式中，P_r—地层静压，MPa；—第i点的实际生产压差，MPa。

各分井段的储层岩石骨架体积计算公式如式（14）所示：

(14)

式中，TV_(i)为第i井段储层岩石骨架体积，m³；R_w为井筒半径，m；R_m为地层出砂半径，m。

流体入流强度计算公式如式（15)、(16)、(17)、(18）所示：

(15)

(16)

(17)

(18)

式中，Q_o(i)为i段处油入流强度，m³/(d·m）；Q_w（i）为i段处水入流强度，m³/（d·m）；Q_g（i）为i段处气体入流强度，m³/（d·m）；Q_l(i)为i段处流体入流强度，m³/（d·m）；k_（i）为i井段内的全部k_0j的平均值；

总流速计算公式如式（19）所示：

(19)

V_f(i)为流体入流流速，m/min；D_w为井筒直径，m。

S4：根据各网格段流体入流流速和岩石内聚强度计算出砂速率和含砂浓度，形成出砂强度非均匀分布剖面；计算综合入流系数得到综合入流非均匀剖面；

根据本发明优选的，所述步骤S4中，计算全井段的各段处的V_f(i)与各测井深度的S_0j，再通过全井段计算平均值得到平均流体入流流速V_fa和平均内聚强度S_0a，所述的V_fa、S_0a分别为全井段的V_f(i)和S_0j的算术平均值。

拟出砂指数计算公式如式（20）所示：

(20)

式中，J_xs为拟出砂指数，m/s；L_qs为出砂强度，L/（d·m）；U_f为流体黏度，mPa.s。

根据入流条件计算出砂速率和含砂浓度，如式（21）所示：

(21)

Q_s（i）为i井段处出砂速率，L/（d·m）；S_0（i）为第i井段的平均内聚强度，为第i井段内的全部S_0j的平均值；

含砂浓度计算公式如式（22）所示：

(22)

C_s（i）为i井段处含砂浓度，无量纲。

综合入流系数计算公式如式（23）所示：

(23)

式中，Z_（i）为综合入流系数，无量纲；Q_sa为全井段平均出砂速率，L/（d·m），计算方式为对各井段的计算平均值。

S5：绘制初始投产状态的水平段储层岩石内聚强度、孔隙度、渗透率、油气水入流强度、出砂强度、液体含砂率、综合入流指数等指标的非均匀分布剖面图；

根据本发明优选的，所述步骤S5中，以水平井井深为横坐标，分别以水平井各井段储层岩石内聚强度、孔隙度、渗透率、油气水入流强度、出砂强度、液体含砂率、综合入流指数等指标为纵坐标，以图2为例绘制非均匀分布剖面折线图。

S6：设定生产时间步长∆t。在初始时刻（t=0）各指标非均匀入流剖面基础上，计算∆t时间后的各网格段上的各入流指标及其分布。首先计算∆t时间段内各网格段的出砂量，根据井段各网格段的出砂量计算储层孔隙度变化及变化后的孔隙度、渗透率、内聚强度；

根据本发明优选的，所述步骤S6中，∆t时间内i井段网格的累积出砂量计算公式如式（24）所示：

(24)

式中，ΔV_s（i）为∆t时间内i井段累积出砂量，m³。

根据出砂量计算∆t时刻末孔隙度，计算公式如式（25）所示：

(25)

式中，为i井段∆t时刻末近井地层孔隙度，无量纲；为i井段∆t时刻初近井地层孔隙度，无量纲。在第一次迭代中，等于步骤S1计算的，第一次以后的每一次迭代循环中，每个∆t时刻初的等于上一个∆t时刻末的的计算结果；

孔隙度、渗透率、内聚强度变化幅度参数计算公式如式（26)、(27)、(28）所示：

(26)

(27)

(28)

为∆t时刻末i井段孔隙度变化率幅度参数，无量纲；为∆t时刻末i井段渗透率变化率幅度参数，无量纲；为∆t时刻末i井段内聚强度变化率幅度参数，无量纲。

计算渗透率变化，按照孔隙度变化比折算，计算公式如式（29）所示：

(29)

式中，k_（i,t1）为i井段在∆t时刻末近井地层渗透率，mD；k_（i,t0）为i井段在∆t时刻初近井地层渗透率，mD，在第一次迭代中，k_（i,t0）等于步骤S3计算的k_（i），第一次以后的每一次迭代循环中，每个∆t时刻初的等于上一个∆t时刻末的k_（i,t1）的计算结果；

计算内聚强度变化，如式（30）所示：

(30)

式中，S_{0（i，t1）}为i井段在∆t时刻末近井地层内聚强度，MPa；S_{0（i，t0）}为i井段在∆t时刻初近井地层内聚强度，MPa，在第一次迭代中，S_{0（i，t0）}等于步骤S4计算的S_0（i），第一次以后的每一次迭代循环中，每个∆t时刻初的S_{0（i，t0）}等于上一个∆t时刻末的S_{0（i，t1）}的计算结果；

S7：根据初始时刻（t=0）入流剖面重新计算生产段井筒各网格段压力分布，使用地层静压减去井筒压力得到各网格段的储层入流生产压差；根据∆t时刻的孔隙度、渗透率计算新的网格段的油、气、水入流强度，形成新的入流剖面；

根据本发明优选的，所述步骤S7中，采用初始时刻（t=0）时刻各井段井筒入流强度计算井筒压力分布。水平段各点的实际井底流压的计算公式如式（31）所示

(31)

式中，P_{wf（i，t1）}为Δt时刻末第i井段实际井底流压，MPa；P_{wf（i-1，t1）}为Δt时刻末第i-1井段实际井底流压，Mpa；ΔP_{k（i，t1）}为Δt时刻末第i井段单位长度流动压力损失，MPa/m。ΔP_{k（i，t1）}为根据初始时刻（t=0）入流量采用流动摩阻压降公式计算得到，第i井段实际井筒流量为该井段至趾端全部网格入流量之和；

水平段各井段的实际生产压差的计算公式如式（32）所示：

(32)

式中，ΔP_（i，t1）为Δt时刻末第i井段的实际生产压差，MPa。

∆t时刻流体产出速度计算公式如式（33)、(34)、(35)、(36）所示：

式中，Q_o(i,t1)为Δt时刻末i井段处油入流强度，m³/（d·m）；Q_w(i,t1)为Δt时刻末i井段处水入流强度，m³/（d·m）；Q_g(i,t1)为Δt时刻末i井段处气体入流强度，m³/（d·m）；Q_l(i,t1)为Δt时刻末i井段处流体入流强度，m³/（d·m）；k_(i,t1)为Δt时刻末i井段处各渗透率平均值。

总流速计算公式如式（37）所示：

(37)

V_f(i,t1)为Δt时刻末i井段处流体入流流速，m/min；Da为i井段处井筒直径，m。

S8：根据∆t时刻新的入流强度、内聚强度计算新的各网格段的出砂量、含砂浓度和综合入流指数，得到∆t时刻全部入流指标新的非均匀分布剖面图；

根据本发明优选的，所述步骤S8中，计算全井段平均流体入流流速和内聚强度，得到Δt时刻末平均流体入流流速V_fa(t1)和平均内聚强度S_0a(t1)，V_fa(t1)和S_0a(t1)分别为全部井段的V_f(i,t1)和S_0(i,t1)的算术平均值。

出砂指数计算公式如式（38）所示：

(38)

式中，J_xs(t1)为∆t时刻末出砂指数。L_qs为出砂强度，L/（d·m）；U_f为流体黏度，mPa.s。

根据入流条件计算出砂速率和含砂浓度，如式（39）所示：

(39)

Q_s(i,t1)为∆t时刻末i井段处出砂速率，L/（d·m）。

含砂浓度计算公式如式（40）所示：

(40)

C_s(i,t1)为Δt时刻末i井段处含砂浓度，无量纲。

综合入流系数计算公式如式（41）所示：

(41)

式中，Z_(i,t1)为Δt时刻末i井段处综合入流系数，无量纲；Q_sa(t1)为Δt时刻全井段平均出砂速率，L/（d·m），对全井段的Q_s(i,t1)取平均值。

S9：计算局部入流系数（Partial Inflow Factor，简称PIF系数）。

根据本发明优选的，所述步骤S9中，PIF系数计算公式如式（42)、(43）所示：

式中，PIF_s为出砂强度局部入流系数，无量纲；PIF_l为产液强度局部入流系数，无量纲；Q_smax(t1)为Δt时刻末全井段最高出砂速率，L/（d·m），V_fmax(t1)为Δt时刻末全井段最高流体入流流速，m/min，对全井各点处的Q_s(i,t1)和V_f(i,t1)取最高值。

S10：生产时间递增∆t，从S6开始继续下一个时间步长的迭代计算，循环实现砂液产出剖面的动态演化，并实现任意时刻t的各种入流剖面预测。

本发明的有益效果：

1. 易出砂储层水平井控砂和控水是保证高效生产的关键技术，由于缺少非均质储层水平井砂液产出剖面动态演化模拟与预测方法，导致其控砂控水措施缺乏目的性。利用本发明提供的方法，可以简易快速模拟砂液产出剖面演化过程，预测给定生产条件和生产时间下的地层流体（油、气、水）和地层砂产出剖面，为水平井控砂和控水管柱设计提供重点产出部位，即重点控制井段和位置，提供控砂控水的目的性和综合效果。

2. 利用本发明提供的方法，可以快速识别非均质储层水平井的高速入流位置，判断防砂控水完井筛管的高冲蚀发生位置，指导优化筛管结构，在高冲蚀风险位置使用高抗冲蚀的筛管，提高完井筛管的综合抗冲蚀能力。另外，利用本发明提供的局部入流系数（PIF系数）可以快速计算实际局部高速入流位置的最高流速，用于冲蚀损坏速率预测及预防。

3. 本发明提供的方法将对水平井出砂和出水剖面的判断识别由静态提升到动态层面，为解决长井段水平井砂水协同控制难题提供有效支撑；相比于静态的砂水剖面识别，防砂效率大幅度提高。并可为智能控砂控水装置提供参数设计依据，随地层生产变化进行更为合理的出砂出水干预，以最低的成本来避免未来可能出现的油井大量出砂/出水情况，避免了出砂出水后的补救作业以及轻微出砂井段的过度防砂。对于油气储层长期开发方式动态优化、单井生产制度实时优化、提高采收率方式动态调整，具有普遍参考价值和科学意义。