0 引言
虽然研究人员通过序列切片成像法［5］、激光扫描共聚焦显微镜法［7］和X射线衍射CT扫描法［2］等岩心建模方法，将碳酸盐岩岩石碎屑颗粒之间的孔隙与岩石内部因溶蚀作用形成的溶孔简单叠加生成孔隙—溶孔型双重介质模型，但这些方法建立的岩心尺度与真实岩心尺度相比均太小［8］。核磁共振实验扫描的范围可以达到岩心尺度，但受限于实验设备精度的影响，其成像图像只能窥探孔喉大概图像，不能精细表征孔隙—溶孔型碳酸盐岩岩石特征［9］。上述几种建模方法在一定程度上对特征孔隙尺寸跨度大、连通性复杂的碳酸盐岩储层局部特征描述仍不够精确，包含准确溶孔特征的模型尺寸常常为微米—毫米级别，与真实更大尺度的厘米—米级别的岩石模型仍存在差异。卷积神经网络［10-11］具有良好的泛化能力，同时在图像特征提取和建模上具有良好的效果，并且可以结合孔隙网络模型，建立孔喉特征更准确的大尺度孔隙—溶孔双重介质模型。
1 孔隙—溶孔型碳酸盐岩建模方法
1.1 碳酸盐岩孔隙与溶孔发育特征的定量表征
表1 岩心样品L16基本物性参数Table 1 Table of properties for the core sample L16 |
| 长度 /cm | 直径 /cm | 平均喉道长度 （l）/μm | 配位数 （z） | 孔隙度/% | 渗透率/（10-3 µm2） |
|---|---|---|---|---|---|
| 实验/模拟 | 实验/模拟 | ||||
| 5.01 | 2.53 | 168 | 4.8 | 7.38/7.20 | 1.98/1.75 |
采用分辨率5 μm的设备对该样品进行CT扫描测试，岩心重构后得到如图3所示的三维图像。其中图3（a）为岩石全貌；图3（b）为在整个岩心中提取一小部分用于研究分析；图3（c）为提取的孔喉空间；图3（d）为建立的孔隙与喉道的球棍孔隙网络模型，红色球为孔隙，白色棍棒为喉道，模型区域大小为700 μm×700 μm×700 μm。图3表明，岩石L16的溶孔团主要是由多个相关的大尺寸孔隙与溶孔组合形成的。为定量表征溶孔的发育特征，采用溶孔相关长度L c来表征平均溶孔团的大小，其定义式为：L c=R c/L，其中：R c为溶孔团半径，cm；L为岩心长度，cm；L c越大，发育的溶孔团越大。统计后得到岩心L16的溶孔相关长度为0.02。
1.2 碳酸盐岩孔隙网络模型建模方法
本文研究采用三维简单立方体SC网络模型［14］为基础构建碳酸盐岩孔隙—溶孔模型。模型总结节点数为N，每个节点代表一个孔隙，孔隙之间通过半径大小各异的喉道相连，喉道长度为l；将每个节点坐标值加上一个范围于［-0.5l，0.5l］的随机数，使每个节点的坐标在半径为0.5l的球型区域内随机移动，生成一个孔喉位置随机无序分布的网络模型。通过图2所示的孔喉半径分布曲线为每个半径赋值，形成与孔喉半径相关的三维稳定随机场张量数据体，更详细的赋值步骤可见BERNABÉ等［20］的研究。配位数z表示与中心孔隙相连的喉道数量，在SC网络中，每个孔隙最多可与6个相邻孔隙相连，即最大配位数z max=6，通过将（1－z/z max）✕N数量的喉道半径设为0去除部分不连通的喉道，即可生成配位数为z的不规则孔隙网络模型，l与z取表1中的值。
在上述构建的无序网络模型的基础上，需对孔喉分布进行重组，从而进一步反映碳酸盐岩孔隙与溶孔的空间分布特征。采用地质统计学中的协方差函数构建卷积核函数［21］：E（h）=exp（-2h/L c），式中h为卷积核中心点到周围空间各个点之间的距离。利用卷积核滑过三维稳定随机场张量数据体的每一个元素进行卷积计算，得到一个新的包含了溶孔空间分布信息的三维张量数据体（溶孔也就是许多大孔隙通过卷积核函数汇聚到一块的），再将该数据体赋值到无序网络模型中，即得到如图4所示的具有不同L c的孔隙—溶孔双重介质孔隙网络模型（模型尺寸为5 cm×2.5 cm×2.5 cm，N =300×150×150，不同颜色表示不同半径的喉道）。图中蓝色区域包裹成团状分布的团块即为溶孔发育带，其余紫色空间可视为基质与微孔隙。与图1对比可以看到，利用L c=0.02建立孔隙网络模型［图4（a）］与实际岩样溶孔尺寸与空间分布具有较好的相似性，其孔隙度和渗透率由BERNABÉ等［22］提出的方法计算，计算值见表1，与实验数据非常接近，这表明本文提出的双重介质孔隙网络建模方法是有效的。图4（b）与图4（c）是通过设置不同的相关长度构建的，用以研究L c对气水相对渗透率曲线的影响。
2 孔隙网络模型的气水两相渗流模拟技术
对于气驱水过程，模型初始完全饱和黏度为μ w=1 mPa·s的湿相流体水，模拟开始后，黏度为μ g=0.017 mPa·s的非湿相氮气从入口端以一定的速度注入，保持模型边界出口端的压力恒定，使注入流体由入口端中心沿孔喉向出口边界流动。任意孔隙i和相邻孔隙j之间的喉道半径为rij，长度为lij。气水两相流动模拟存在以下假设条件：①所有的流体的压降只发生于孔隙之间的喉道中；②喉道内两相流体间只存在一个气水界面；③孔隙空间中只有活塞式驱替（图5）。当气水两相同时存在时，两相弯液面间的毛管压力pcij （Pa）使用杨拉普拉斯方程求解：
式中：γ为界面张力（设置为70 mN/m）；θ为水相接触角（设置为130°）。单个喉道中的体积流量qij （m³/s）用扩展Hagen-Poiseulle方程［18］表示：
式中：gij 为喉道的传导率，m³/（Pa·s）；pi 与pj 分别为孔隙i和孔隙j的压力，Pa；B g为气体体积系数，无量纲；Z和Z sc分别为地下和地面气体偏差因子，无量纲；T和 T sc分别为地下和地面温度，K； 为地下气体压力，Pa；p sc为地面大气压，Pa；μeff 是单个喉道中两相流体的有效黏度，通过下式计算：
式中：xij 是与气水两相凹液面位置有关的无量纲数（0≤xij ≤1），即凹液面所在位置横坐标除以整个喉道的长度（图5）。当喉道中只有单相流体时，pcij =0，μeff =μ g或μ w。
由于气体的强压缩性，流体入口压力无法瞬间从入口端到达出口端，流体的入口和出口流量不能瞬间达到平衡，孔隙与孔隙之间喉道中的流量qij 满足非稳态渗流压力扩散方程［23］：
式中：C t为气水综合压缩系数，Pa-1。
根据式（6） 遍历所有孔隙，通过泰勒展开以及Crank-Nicolson隐式有限差分法技术，可以构建线性方程组或矩阵方程并求解［24］，从而进行动态模拟，计算模型的气、水出口流量。重新计算并更新压力场、体积系数和饱和度，进行下一步计算。
3 模拟结果及讨论
3.1 模拟方法验证
图6 岩心样品L16实验与模拟气水相对渗透率曲线对比Fig.6 Comparison between the simulated and measured gas-water relative permeability for the sample L16 |
表2 不同溶孔特征条件下相对渗透率曲线特征值Table 2 Characteristic values of relative permeability acurves under different dissolved pore characteristics |
| S wi | K rg | K m | S gm | M | |
|---|---|---|---|---|---|
| 实验L c=0.02 | 0.42 | 0.48 | 0.08 | 0.29 | 0.10~0.58 |
| 模拟L c=0.02 | 0.45 | 0.47 | 0.06 | 0.30 | 0.11~0.55 |
| 模拟L c=0.20 | 0.36 | 0.52 | 0.07 | 0.44 | 0.11~0.64 |
| 模拟L c=0.35 | 0.24 | 0.33 | 0.08 | 0.49 | 0.10~0.76 |
3.2 溶孔相关长度L c对气驱水相对渗透率曲线的影响
3.3 注气速度对气驱水的影响
当气体突破出口端时，便很难再驱动剩余孔隙空间中的水，波及效率也很难再提高［26］，因此本文研究选用L c=0.02的模型重点分析了气体突破前不同注入速度对气驱水过程的影响。不同注气速度情况下气水分布情况、注气流量与压力数据存在差异，如图8所示，图中彩色部分是驱替过程中气相占据的不同孔喉半径的空间，空白区域完全含水，其中气水流量曲线横坐标的归一化时间为模拟气驱过程的实际时间与总时间的比值。图8表明，岩心气驱水过程中存在2种注气压力变化模式：①在恒定低速注气时，初始时刻气体憋压现象［17］会导致注气压力先上升后下降，束缚水饱和度较大，这是因为孔径差异引起各孔喉间毛管压力差异，孔喉两端压差需要不断改变以满足克服毛管压力，从而反复波动，最大注气压力为0.11 MPa，平均注气速度为0.002 5 mL/min，气相突破时的饱和度为0.07；②高速驱替时，两相流体快速流动会引起黏性阻力显著增加，需要较高的初始入口压力以克服毛管力和黏性力，从而满足流动条件，随气体向前推进，含水饱和度降低，水相的阻力也降低，初始入口压力就会表现出快速下降且几乎不波动，最大注气压力为0.72 MPa，平均注气速度为0.173 3 mL/min，气相突破时的饱和度为0.14。
前人［26］的实验研究已经认识到，不同的驱替速度（或压力）对非混相驱替过程的注入压力、毛管阻力与黏性阻力之间的动态平衡有重要影响，使驱替结果产生不同程度的差异。相较于常规实验，该研究提出的模拟方法可以直观地展示不同注气速度情况下气水两相流动的差异：驱替速度（或压力）较小时，渗流过程中主要受到毛管力的影响，气体会倾向于寻找毛管阻力最小的路径流动直至突破，从而驱替效率较低。增大注气速度，有助于气相克服毛管阻力和流体黏性阻力，气相流动能力与路径增加，波及效率增加。
4 结论
（1）基于卷积算法建立的孔隙—溶孔双重介质模型能够较真实地反映实际碳酸盐岩岩石的溶孔发育特征；考虑气体压缩性的孔隙网络气水两相非稳态渗流模拟方法能够有效地模拟气水两相渗流过程。
（2）碳酸盐岩储层中溶孔团越发育，气水两相共渗区越大，气水同产的时期越长；受毛管力与黏性力的影响，不同注气速度情况下气水空间分布情况、注气流量与注气压力存在差异。
甘公网安备 62010202000678号
