致密砂岩气藏气水流动规律及储层干化作用机理

赵玉龙 刘香禺 张烈辉 唐洪明 熊 钰 郭晶晶 单保超

1. “油气藏地质及开发工程”国家重点实验室·西南石油大学 2. 煤燃烧国家重点实验室·华中科技大学

0 引言

致密砂岩气藏储层通常具有孔隙度与渗透率低、非均质性强、毛细管现象突出、原始含水饱和度高等特点，储层极易遭受损害，并且储层中广泛存在的地层水使气体渗流受阻，导致致密砂岩气（以下简称致密气）高效开发面临诸多难题[1-5]。在钻井、压裂等作业过程中滤失的工作液加上储层中地层水会使得井筒及其附近储层高含水，对孔喉造成堵塞，减小气体渗流空间，形成所谓的“水锁”现象[6]。而储层干化则可以有效地解决该问题[7-9]，即通过向生产井注入干化剂，使得井筒及其附近一定范围内储层中地层水和侵入的工作液迅速反应气化，让原本被水占据的通道重新成为气体渗流通道，从而提高气体的渗流能力。但是，针对储层高温高压条件下气水两相流动模拟及储层干化对气体渗流能力的影响，目前的相关研究报道还鲜见。

致密气藏中气水两相流动复杂，传统的计算流体力学方法（CFD方法）在处理孔隙尺度两相流动问题时面临诸多挑战。虽然流体体积函数法、水平集方法等可以实现对相界面的追踪和捕获[10]，但是难以从微观粒子间的作用力出发来进行描述，并且上述方法对相界面的表征需要专门进行追踪和捕获计算，费时费力。因此，学者们开始注意到格子Boltzmann方法（以下简称LBM）在孔隙尺度多组分多相流动模拟上的优越性。LBM作为一种较新的数值计算方法，已经被证明是可以用于模拟流体流动和输运现象的方法，在多孔介质流动模拟中得到了广泛应用。不同于传统的CFD方法，LBM基于微观模型和介观动理学方程，能够在不跟踪、捕捉非混溶相或者组分间界面的情况下实现界面自动分离，并且能够表征出微观粒子间的相互作用。这些优势在传统CFD方法中是很难实现的。学者们采用LBM对非常规天然气开展了广泛的流动模拟研究，但大多数研究是针对简单结构通道中的单相气体[11-15]，对复杂多孔介质中的气水两相流动模拟则鲜有报道。

为了研究储层高温高压条件下气水两相流动规律，进而揭示储层干化对致密气渗流能力的影响，笔者基于LBM模拟地层高温高压条件下致密气驱替地层水的流动过程，得到地层中束缚水的分布，然后采用激光刻蚀模型，进行储层干化实验，并借鉴该实验的可视化结果对储层干化数值模拟进行简化，在此基础上利用数值模拟手段研究了储层干化对致密气渗流能力的影响。

1 LBM基本理论

前人基于LBM在模拟多组分多相流动方面已开展了大量研究工作，形成了一系列格子Boltzmann模型[16-21]，笔者此次采用颜色梯度模型。格子Boltzmann模型主要包括格子模型、演化方程和平衡态分布函数，此次研究二维气水两相流动，采用D2Q9格子模型。Grunau等[22]对颜色梯度模型进行改进，实现了不同颜色相的分布函数独立参与演化过程。在两相流颜色梯度模型中，演化方程为：

式中下标k表示不同相（R、B相）的流体；下标i表示离散速度方向，取值为0, 1,, 8；表示格子空间 位置处k相流体粒子在时刻t处于i方向的粒子分布函数；表示碰撞算子，用于表征流体粒子之间的相互作用；表示扰动算子，用于表征表面张力的影响； 表示格子速度；δt表示格子空间时间步。在LBM中所有物理量的单位均为格子空间单位，无量纲。

式中τk表示k相流体的无因次松弛时间；上标eq表示平衡态。

τk和k相流体格子运动黏度（vk）的关系如式（3）所示[24]。在实际模拟中，可通过设置两相无因次松弛时间来模拟不同黏度比。

式中cs表示格子声速，。

平衡态分布函数计算式为：

式中ρk表示k相流体在格子空间的宏观密度；wi表示权重系数；表示格子空间宏观速度。

参数 表达式为：

参数αR、αB与两相密度比的关系式为：

式中ρR、ρB分别表示R、B相流体在格子空间中的宏观密度。

此时满足两相压力相等，有

式中pk表示k相流体压力。

可以通过调整式（6）中αk来满足预定的密度比，并且满足。

权重系数（wi）和离散格子速度的计算式分别为：

其中

式中c表示格子速度。一般取δx=δt，则c=1。

式中Ψ表示表示不同位置处两相差别的局部序参数。

Ψ计算式为：

式中Ak表示与界面张力相关的参数。

参数Si在不同的离散速度方向有不同的取值，有

在原始颜色梯度模型中，粒子分布函数在经过上述演变后，还需要重新进行标色处理，该过程涉及极大化求解问题，需要进行大量计算，费时费力。为此，不少学者提出了简化方法来重新标色粒子分布函数[26-28]，笔者采用Leclaire等[28]提出的简化方法，即

ρk、ρ计算式分别为：

2 模型验证

2.1 Laplace定律验证

Laplace定律是验证两相流动模型和代码正确性的常用定律，此次采用该定律验证上述模型在地层高温高压条件下的准确性。根据该定律，在界面张力的作用下非混相两相流体会分离并形成圆形相界面，界面内外相压差和圆形相界面的半径成反比例关系，即

式中pin表示圆形相界面内部流体圆心处压力；pout表示外部流体压力；σ表示界面张力；R表示圆形相界面半径。

考虑地层温度为400 K，通过查询美国国家标准与技术研究院（NIST）化学数据库，得到压力介于10～100 MPa时地层水和致密气密度比小于20，运动黏度比小于2，与地面条件下的结果相差甚远。在正方形流场中放置一个初始静止气泡，周围为地层水，正方形四周进出口处均设置为周期性边界条件[29]，为验证模型的适应性，设定内外相运动黏度比为1∶50，在密度比为5、15、25、35、45的条件下进行模拟测试。模拟结果表明，在不同密度比下，R2均在0.99左右，满足Laplace定律（图1）。

2.2 平直通道两相Poiseuille流动模拟验证

平直通道中两相Poiseuille流动如图2所示，图2中L、H分别表示通道长度和宽度，a、b分别表示两相流体分界面及通道壁面到通道中心线的距离，垂直于流动方向的横截面上速度剖面满足式（20），即

图1 不同密度比下Δp和

图2 平直通道中两相Poiseuille流动模拟示意图

式中Ux表示流动方向速度；表示合外力，其y方向分力为0；v1、v2分别表示“流体相1”和“流体相2”的运动黏度；ρ1、ρ2分别表示“流体相1”和“流体相2”的密度。

致密砂岩气储层岩石通常具有亲水性，因此考虑图2中“流体相1”为地层水。模拟气水两相黏度比为1∶25，合外力大小为0.000 001，上下边界采用Bounce-back反弹格式来处理，左右边界为周期性边界条件，进而模拟无限长通道中流动[29]。如图3所示，可以看出采用格子Boltzmann模型计算的两相Poiseuille流沿流动方向速度（Ux）的数值解与解析解基本一致，表明该模型可以用于地层条件下气水非混相驱替的模拟。

图3 基于LBM模型计算的两相Poiseuille流速度数值解与解析解结果对比图

3 基于数字岩心的气驱水流动模拟

3.1 数字岩心生成及流动参数设置

多孔介质数字岩心构建是孔隙尺度流动模拟的前提，其构建方法主要分为物理方法重构和数学方法生成两类[30]。常规实验手段无法直接获取岩心内部的孔隙结构信息，通常需要借助CT断层扫描仪或高倍光学显微镜得到一系列有序的高精度二维切片，然后将切片叠加从而构建三维空间上的数字岩心模型，即物理方法重构。虽然该方法可以得到分辨率较高的数字岩心模型，但是此类数字岩心切片在二维平面上不具有连通性，这显然对于二维空间上的流动模拟并不适用。笔者是为了研究气水两相流动规律，并不要求二维数字岩心的形状和真实孔隙空间完全一致。因此，采用数学方法直接生成连通的孔隙网络模型是一种可行的方法。直接生成数字岩心的方法有多种[29]，其中四参数随机生长法与LBM相结合在孔隙尺度的研究中应用广泛[31]，笔者此次也采用四参数随机生长法构建孔隙网络模型。

对于多孔介质储层，一般水平方向渗透率高于垂直方向渗透率，在构建多孔介质时应该考虑固体相在水平方向有更大的生长概率。基于前述方法构建多孔介质模型，流场x、y方向分别设置220个、150个网格，其中入口和出口边界处各预留10个网格宽度区域作为缓冲区域，剩余中间200个网格宽度的区域为多孔介质。为确保随机生成的多孔介质在流动方向上具有一定程度的连通性，设置孔隙度为50%。根据达西定律计算得到多孔介质渗透率，其计算式为：

式中Kw表示真实物理空间渗透率，m2，1 m2≈1012D；λw表示真实物理空间长度，m；λ表示λw对应格子空间中的长度。

以格子空间中单位长度对应真实物理空间中10 nm长度为例，此时的渗透率约为0.04 mD。入口端采用非平衡态外推格式边界条件模拟稳定驱替相源[32]，出口端采用对流形式的开放边界条件[33]，上下边界和固体壁面采用半步长反弹格式处理[29]。

考虑地层温度为400 K，压力为40 MPa，通过NIST化学数据库可以查到地层水和致密气（纯甲烷）的密度及黏度（表1）。模拟时将致密气相密度设置为1.00格子单位密度，地层水相密度设置为5.45格子单位密度，取致密气在格子空间中的运动黏度为0.166 667，地层水运动黏度为0.291 936，即可使得致密气和地层水在给定温度、压力下的主要物性参数之比在格子空间和真实物理空间中相等。对于两相流动还应该考虑其在固体壁面润湿性的不同所带来的影响，本文通过不混溶两相接触角设置方法[34]，模拟地层水与固体壁面的接触角分别为30°、45°、60°和75°的情况，采用多核CPU并行计算，模拟计算时间步为500 000步。

表1 储层条件下地层水和致密气的物性参数表

3.2 结果分析

计算得到不同时间步多孔介质中含水饱和度（Sw）的变化，如图4所示。可以看到整个驱替过程中含水饱和度逐渐下降，并且地层水和孔壁接触角越大，含水饱和度下降速度越快；在130 000步附近存在明显转折点，该转折点恰好与致密气驱替前缘突破多孔介质相对应。当致密气突破后，气体主要沿连通优势孔道流动，含水饱和度下降的速度将减缓。

图4 不同接触角下气驱水过程中多孔介质含水饱和度（Sw）随时间步变化曲线图

图5 气驱水过程中气水两相分布图（接触角为30e）

图5展示了地层水和孔壁接触角为30°时致密气驱替地层水在不同时间步的气水两相分布情况。图5中蓝色区域代表地层水、黄色区域代表致密气、黑色区域代表岩石骨架。可以看出，随着驱替时间步增加，中部连通的大孔道中的地层水首先被驱替。当致密气在连通的大孔道中突破多孔介质外流时，上下相邻连通的小孔道中的地层水仅有少量被驱替；继续增加驱替时间，上部连通的小孔道中致密气突破。此后，致密气主要集中在上部和中部优势通道，中部和下部连通的小孔道中驱替几乎停滞。大约437 000步后含水饱和度曲线趋于平缓，流动达到相对稳定状态，含水饱和度为44.47%，与致密气藏含水饱和度高的特点相符。当地层水和孔壁接触角为45°、60°和75°时，在模拟结束前含水饱和度仍在缓慢下降，这是由于在水湿条件下毛细管力是气驱水的阻力，润湿角越大阻力越小，此时致密气仍能驱替部分小孔道中地层水。

根据模拟结果将束缚水分布类型分为盲端孔隙水、死孔隙水、连通通道壁面吸附水膜（以下简称吸附水膜）以及卡断水4类（图6）。存在卡断水的位置将产生额外的毛细管力，阻碍驱替进行甚至无法驱替；吸附水膜则占据孔道空间，使气体流动的有效孔径减小。同时可以看到，多孔介质中大量连通的微小通道被地层水占据，从而无法实现有效驱替，存在明显的“水锁”现象，严重影响了致密气在多孔介质中的渗流能力。

图6 驱替结束后束缚水分布类型划分图

4 储层干化对致密气流动影响研究

4.1 储层干化对地层水的消耗实验研究

在开展激光刻蚀薄片驱替实验以后，研究储层干化对地层水的消耗作用。实验选取厚度为2 mm，长度、宽度分别为20 mm、10 mm的环氧树脂基材料为刻蚀试样，限于现有工艺技术水平，多孔介质区域刻蚀深度为50 μm。利用激光刻蚀技术将设计好的孔隙网络模型在微小透明面板上重构，得到刻蚀薄片模型。透明刻蚀薄片置于高倍显微镜载物台上，在显微镜上方加装摄像机对实验过程进行拍摄，同时将数据传输到电脑并实时显示（图7）。

图7 激光刻蚀薄片驱替实验装置示意图

首先，向刻蚀薄片模型中注入甲烷气，将其中的空气排出，然后采用微流量泵以0.05 mL/min的速度注入配置的地层水。为便于观察，在地层水中添加微量蓝色染剂。由于模型孔道为亲水的，在地层水驱替甲烷气达到稳定后，模型大部分孔道被地层水占据，只有少量孔道为甲烷气与地层水共存。利用微流量泵向薄片入口端以0.05 mL/min的速度注入甲烷气体，直到出口端无地层水流出。最后，将2 PV提前配置好的干化剂从出口端以同样的速度反向注入，此时干化剂与束缚水接触反应。

如图8-a所示，盲端孔隙中存在部分气泡，无法实现地层水的完全饱和，这是由界面张力和贾敏效应引起的；如图8-b所示，束缚水的存在形式与前述数值模拟结果相符，包括盲端孔隙水、死孔隙水、吸附水膜及卡断水；如图8-c所示，干化剂与束缚水接触后发生反应，产生大量气泡，吸附水膜、卡断水以及大量盲端孔隙水被消耗掉，仅剩余死孔隙水以及少量盲端孔隙水。对比图8-b、c，发现干化剂存在明显作用深度，即超过作用深度的盲端孔隙水无法被有效消耗。

4.2 储层干化对致密气流动能力影响的数值模拟研究

注入干化剂主剂为离子型碳化物超细微粉末，通过添加增溶剂使其悬浮于易挥发携带剂中，与水接触反应后消耗地层水，同时产生甲烷气体，从而达到降低含水饱和度的目的。实际上，干化剂与地层水的反应复杂，而在数值模拟中直接添加其化学反应困难；同时，采用室内实验模拟地层高温高压条件下驱替及干化难度大，无法有效获取干化后储层中地层水的分布情况。因此，借鉴激光刻蚀模型中常温常压下干化剂对地层水的消耗过程，对干化数值模拟进行简化，然后借助数值模拟手段研究储层干化对致密气渗流能力的影响。

图8 驱替实验不同阶段激光刻蚀薄片局部区域相分布图

考虑干化发生在气水界面处，随着气水界面逐渐向地层水占据的空间扩展，地层水逐渐被消耗，束缚水饱和度逐渐降低。具体简化过程如下：①首先将前述模拟得到的气驱水结束后的束缚水密度场输入密度矩阵（R），并通过各个格点位置的密度矩阵值判断其所属相态（致密气/地层水/固体骨架颗粒），致密气格点用“1”表示，地层水格点用“0”表示，固体骨架颗粒格点用“－1”表示，进而得到三相分布矩阵（S）；②遍历致密气格点，判断其是否位于气水两相界面处，如果是则赋值为“1”，其余所有格点赋值为“0”，进而得到气水两相界面分布矩阵（F）；③遍历气水两相界面格点，由这部分格点向其邻近的地层水格点进行干化拓展处理；遍历S中地层水格点，如果与气水两相界面相邻，将S中相应位置格点由“0”转换为“1”；④重复步骤②、③，直到完成T次干化拓展，相应将干化拓展深度表示为D(T)。

图9 不同干化拓展次数（T）下气水两相分布图

将前述致密气驱替地层水在模拟时间步为500 000步的气水两相分布作为干化拓展处理时的初始气水分布，干化效果如图9所示。可以看出，随着干化拓展次数（T）增加，干化对地层水的消耗越来越多，原本被地层水占据的通道逐渐被打开，成为气体流动通道；干化拓展次数越多，气水两相界面推进得越深。

为了研究干化对致密气流动的影响，仍然采用上述模型，基于不同干化拓展次数得到的气水两相分布，进行致密气驱替地层水流动模拟。模拟参数与前述3.1部分设置的参数相同，模拟时间步为200 000步。

如图10所示，储层干化后气体参与流动的区域范围明显增大，同时气相流场中的最高流速也明显提高，并且高速流动是出现在连通的大通道的狭小喉道处。图11中UT,ave表示不同干化拓展次数（T）下的沿程平均格子速度，U0,ave表示不考虑储层干化时的沿程平均格子速度，可以看到在相同x位置处，随着干化拓展次数增加，增大，说明干化可以有效提升储层多孔介质中气体渗流速度。如图11-b所示，在＜0.7范围内，干化拓展次数（T）等于50和100对应的基本一致，并且均大于T等于10时的，说明在范围内、T等于50和100情况下的储层干化效果基本相同，并且已经达到极限；＞0.7后，T等于100对应的明显大于T等于50对应的，图9-a显示该区域多孔介质中存在大片“水锁”区域，说明在此范围内干化作用越强对于解除“水锁”效果越好。

为了进一步表征储层干化对致密气渗流能力的影响，计算出不同T下致密气渗透率提升倍数（Kr）：随变化曲线图

图10 储层干化前后多孔介质气相速度场云图

图11 不同干化拓展次数（T）下UT,ave与

式中KT表示经过T次干化拓展后获得的致密气渗透率；K0表示不考虑干化获得的致密气渗透率。

如图12所示，随着T增加Kr整体呈现上升趋势。在T较小（介于0～20）时，储层干化对致密气渗流能力的提升效果显著，但T大于20后，Kr增幅变缓，致密气渗流能力的提升效果将不显著。T=20时，Kr=1.37，即干化后致密气渗透率增大了37%；而当T增至80时Kr仅为1.5，较前者（1.37）仅提升了9.49%。结合上述认识，在进行储层干化工艺的现场应用时，需确定出合理的干化强度，从而能够获得较好的经济效益。

图12 KrüT关系曲线图

5 结论

1）所建立的格子Boltzmann模型在地层高温高压条件下满足Laplace定律，采用该模型得到的两相Poiseuille流速度数值解与解析解结果基本一致，表明该模型可以用于地层条件下气水非混相驱替的模拟。

2）致密气在多孔介质连通的大孔道中优先突破，并且在突破后驱替地层水的速度显著下降；地层水与岩石壁面的接触角显著影响气水两相流动，岩石亲水性越强驱替速度越慢。

3）致密砂岩气藏中束缚水可分为吸附水膜、盲端孔隙水、死孔隙水和卡断水4类；多孔介质中大量连通的微小通道被卡断水和吸附水膜占据，存在明显的“水锁”现象，严重影响致密气在储层多孔介质中的渗流能力。

4）干化剂可与束缚水反应并且产生大量气泡，将吸附水膜、卡断水和盲端孔隙水消耗掉从而提高气体渗流能力。

5）对于由卡断水形成的“水锁”区域，增大干化强度可以有效改善气体渗流能力，整体上随干化强度增大，致密气渗透率增大，但干化强度超过一定程度后，致密气渗透率的增幅逐渐减小。