黏土膨胀储层伤害数值模拟研究

黄波， 徐建平， 蒋官澄，3， 王巧智， 苏延辉， 白相双

（1.中海油能源发展股份有限公司工程技术分公司， 天津 300450；2.中国石油大学（北京）石油工程学院， 北京 102249；

3.中国石油大学（北京）油气资源与探测国家重点实验室，北京102249；4.中国石油吉林油田钻井工艺研究院， 吉林松原138000）

0 引言

钻井、 完井、 注水等过程向储层中引入水相，往往会导致黏土水化膨胀储层伤害。黏土膨胀储层伤害是一类重要的储层伤害类型[1-3]。从矿场实践来看，目前中海油各区块的一些钻、完井，以及注水井不同程度地出现了黏土膨胀伤害，因此研究黏土膨胀的时间、空间性质具有较高的现实意义。从黏土膨胀机理来看，目前主要研究的有2类，一类是渗透水化[1，4-5]，即孔隙中水分子通过固-液界面向岩石固相内部渗透扩散，导致岩石内部黏土矿物吸水膨胀，损害介质孔隙度和渗透率；另一类是表面水化[6]，黏土矿物表面含有强亲水性的羟基组和含氧官能团，在液体溶液体系中，溶液中的水分子会被这些亲水基团吸附到颗粒表面形成一层有序水分子边界层，即水化膜[6-7]。如果不考虑表面水化，则等价于认为这层水化膜的厚度为0。从储层伤害建模研究的角度，一般是考虑渗透水化[1]。而对于表面水化，往往用于钻井液造浆[8-9]、高泥化煤泥水高效沉降澄清处理工艺[7]等，其共同特点是黏土固相分散到溶液体系中，是液体主导的系统，而从黏土膨胀储层伤害的角度上看，则是溶液体系在孔隙固相介质中渗流和扩散，是固体主导的系统，这是不一样的物理过程，因此在该模型的建立中主要考虑渗透水化[1]。

从研究方法上看，黏土膨胀储层伤害的研究主要集中在实验和数值模拟2方面。实验方面主要包括黏土稳定剂的研发[10]、井壁稳定性的研究[11]、膨胀特性实验研究[12]等；数值模拟方面较流行的有采用分子动力模拟技术来研究黏土矿物晶体层间结构[13]等。总体来看，黏土膨胀储层伤害的宏观时空模拟研究不多，但对于矿场实践，需要对黏土膨胀伤害的空间分布特征和其随时间变化特征有一定的定量或半定量认识。在黏土膨胀储层伤害的时空模拟中，膨胀的限度和速率是2个核心问题。可以比较容易地建立水分在储层孔隙中运动的数学模型，但水分和黏土接触之后黏土的膨胀行为的量化描述却是相对更难的一件事。笔者建立了黏土膨胀储层伤害的时空模型。当然，任何模型都是有简化和适用条件的，这里也不例外。

1 模型的建立

1.1 局部膨胀建模

黏土膨胀损害是比较复杂的过程。但其基本的物理过程大致分为3个阶段，①水相进入储层内部孔隙空间渗流；②水相和孔隙岩石壁面接触，开始向固相内部扩散，进入固相的水分导致固相内的黏土成分开始膨胀；③驱动水相进入固相的动力和阻力平衡，水分无法再进入介质固相内部，膨胀趋于结束。第1个阶段比较容易描述，只需要建立常规水相对流扩散的渗流方程即可。而后2者是黏土膨胀储层伤害模型的关键与核心。只要建立了第2阶段的模型，第3阶段便只是第2阶段的渐进情况。

首先关注水分子通过固-液界面向固相介质内部扩散的过程。出于简便考虑，假设水分子通过固-液界面由孔隙中的液相向岩石壁面内部扩散的过程服从菲克扩散。先设几个重要的变量：岩石固相中的初始含水体积分数c0、岩石固相中某一时刻的含水体积分数c、孔隙中某一时刻的含水体积分数c1、介质初始孔隙度φ0、岩石孔隙中的束缚水饱和度 Swc。则显而易见的是 c1（t=0）=Swc。在某个微观孔隙壁面的水侵过程可以由图1来表达[1]。

图1 储层多孔介质中水相渗流及固-液界面水分子渗透示意图

如图1所示，取指向固相内部为z方向，亦即水分扩散的方向。根据菲克扩散律可以建立水分子在壁面向内部扩散的方程[1，14]：

式中，D为水分子向固相内部扩散的扩散系数。这个方程的初始条件即为c（t=0）=c0。边界条件则为中S为累积进入固相的水体积分数，则为水吸速率。系数k即为著名的膜交换系数（Film Mass Transfer Coefficient）[1，15]， 它和溶液中离子相互作用及矿化度等有关[15]。在这样的边界条件和初始条件下， 可以求解方程（1）， 对原方程和边界、 初始条件同时进行拉普拉斯变换得到该系统在拉普拉斯变量p域上的解。

对上式进行拉普拉斯逆变换可以得出待求解c（z，t）。通过查表[14]或者一般性围道积分逆变换公式[16]即可得出 c（z，t），

这个结果与文献[1]一致。在式 （4）推导过程中用到了如下的数学事实，即：

式（5）的推导可以通过简单的数学操作得出。基于式（5）有2个必要的结果，

由式（6）就可以推出式（4）。在这里，D刻画了水分向固相内部运动的动力；而k则刻画了水分向固相内部运动的阻力。据文献[17]可知，黏土膨胀导致的孔隙度变化率可以由如下方程来表征。

其中，λ为黏土膨胀系数。λ满足下述关系[17]：

其中，PI为岩石塑性系数，无量纲，若PI＜（1～2），则为脆性岩石；2＜PI＜6，塑脆性岩石；PI＞6则为塑性岩石。Cc为岩石中的黏土质量分数。k′为经验参数。式（7）表明，只要＞0，黏土就会膨胀，使得，即孔隙度下降。决定正负的就是c1，c0的相对大小。若c1＞c0，则＞0，说明孔隙中水分含量大于固相中水分含量， 将向固相内扩散；反之， 若 c1＜c0， 按照式（4）、 式（7）将给出＜0，这意味着孔隙度增加，似乎说明出现某种岩石失水而收缩的情形，实际上笔者认为这是不合实际的。假设在人为钻完井或注水改变储层之前，孔隙和基质间的水迁移早已处于平衡状态的话，＜0意味着某种从孔隙中向外抽出了水分，打破平衡之后岩石基质向孔隙中失水的情形。而在研究的这些过程中都是在向孔隙中补充水。所以在模型中对于加上限制， 当 c1＜c0时， 令=0，即没有孔隙中的水相向固相内扩散，这是显而易见的，水分子不会从低浓度往高浓度自发运动。对于钻完井和注水过程，外界引入的水相是增加了孔隙中的水分含量，而使得c1＞c0。因此有下式。

该公式意味着，外界引入的水使黏土膨胀。如果膨胀，则其膨胀的速率和限度由驱动水分子运动的动力D和阻力k决定。例如， 在k=1.7×10-6m/s，D=1×10-7m2/s， c1=0.08， c0=0.05 时， 黏土吸水速率随时间变化趋势如图2所示。需要注意的是，严格意义上讲c1是场而不是标量，所以图2表示的应是某个空间点处的水吸速率随时间的变化趋势。

图2 黏土吸水速率随时间的变化趋势

由图2可知， 其吸水速率从初始时刻的5×10-8m/s，到后期最终衰减为0。其特点是早期剧烈衰减，中晚期衰减趋势平缓。这就决定了黏土膨胀这种储层伤害类型具有早期迅速增大（如果能发生）， 中后期增大趋势逐渐变缓趋于停止的特征。膨胀总量几乎由早期贡献， 但伤害的总时间可以延伸较长的时期。当然， 参数值不同会导致此曲线的一些差异。但其趋势规律是有广泛意义的。目前中海油各区块已确认发生黏土膨胀的井中， 确实出现了和这种趋势一致的行为。尤其是注水井， 因为注水井可以注很长时间， 可跨越黏土膨胀损害整个生命周期。比如绥中36-1油田M13注水井， 从2014年6月新井投注至2015年10月试井测表皮注水时间约480 d，实测总表皮为8.38， 黏土膨胀伤害占该井伤害量的近一半。该井到2015年1月开始出现欠注， 注水压力上升到10 MPa，从2015年1月到4月较短时期内配注量从600 m3/d快速下降到340 m3/d， 此后直到约第480 d试井发现储层污染期间，配注量自340 m3/d缓慢下降到200～300 m3/d之间。这种趋势基本说明，黏土膨胀伤害大致具有这种早期剧烈增加，中晚期缓慢增加直至停止的特征。但由于现场难以获取定量的单因素损害数据，这种“验证”不完全严谨。而对于钻完井过程， 其特征时间为几天到几十天， 主要位于黏土膨胀早期， 处在损害快速上升的阶段。比如LF13-1油田10a井， 钻井耗时90 d， 完钻测试达30的表皮系数， 据估计该井的黏土膨胀损害也占到约30%～40%的比例，在较短的时期内就上升到了可观的损害程度。这时如果积极采取防膨措施，可以基本杜绝这种损害的发生。

1.2 黏土膨胀伤害时空模拟

由于人为因素导致外来水相进入储层渗流， 使得c1成为一个时空分布的变量。所以模型的另一部分是对c1的描述。外来水相在储层中的运动满足对流扩散律，不难写出c1满足的方程：

其中，K0为渗透率，mD，μ为流体黏度，P为压力，Pa。另有压力传导方程为：

其中，Ctotal为岩石流体综合压缩系数。当得出了孔隙度的变化情况之后，可以得出无因次渗透率和孔隙度之间的经验关系如下[1，18]。

φ（r，t），Kd（r，t）表明它们是随时间和空间变化的场函数。这里给出一个算例。模拟的是注水井。算例边界条件为 c1（r=rw， t）=（表明注入点处孔隙完全被水充满）；初始条件为为束缚水饱和度。模型的必要输入参数列在表1中。

表1 模拟所需参数取值列表

模拟井周时，需要采用柱坐标系，即有下式。

对空间采用中心差分、对时间隐式差分，可以较容易地求解此系统。之所以选择模拟注水井是因为注水井注水时间一般很长，可以跨越黏土膨胀损害从产生到逐渐结束的整个时间段。实际上，损害的严重程度跟其他参数也密切相关。不同的注水井损害严重程度不同，同时有的钻井过程，虽然时间比注水过程短得多，但它造成的损害程度可能比注水过程更大。比如前面提到的LF13-1油田10a井，钻井90 d导致的储层伤害就比绥中36-1油田M13注水井注水480 d导致的储层伤害还更严重。

2 模拟结果

基于模型（13）和表1的输入参数，可以模拟出这个算例的黏土膨胀储层伤害导致的无因次渗透率空间分布的变化情况，包括表皮系数的变化趋势，见图3和图4。这里的表皮大小只是基于算例而言的。若用户采用此模型需要输入当地参数值。

图3 黏土膨胀伤害模型无因次渗透率空间分布随时间变化的趋势（基于算例）

图4 黏土膨胀伤害模型表皮系数随时间的变化趋势（基于算例）

从图3和图4可知，黏土膨胀伤害特点是某空间点的损害率不一定特别大，但其空间扩展的范围可以较远，是不可忽略的损害因素。同时，从其表皮系数的增加趋势也可以看出，早期快速增长，而中晚期缓慢增加直至不变。这样就可以大致了解一定输入条件下黏土膨胀伤害的速率和限度，对实际产生一定的指导意义。

最后，针对模拟结果的实验验证问题，目前确实有很大的困难。难点主要有2个方面， 一是该实验验证应是在矿场尺度， 需要组织大量的人力物力来完成， 目前的研究工作允许的是室内岩心实验，但是岩心实验无法代表矿场尺度的伤害空间分布情况及其随时间变化特征， 因此对该模型的验证作用不大；二是来自于问题本身的困难， 该研究是对黏土膨胀储层伤害单独进行的建模研究， 而现场试井测试给出的表皮系数是近井地带伤害的综合反映，难以定量分解某种损害类型占总表皮的比例， 针对该定量化的实验研究当前不太现实。事实上， 目前在现场遇到的储层伤害诊断问题中， 如何科学地将总表皮分解到各损害类型是亟待解决的难题。该模型只是笔者尝试解决该问题的一部分尝试。近期，开发的储层伤害诊断软件“砂岩储层伤害预测诊断及储层保护系统” 在现场进行了试应用， 将软件计算表皮结果同实测表皮对比以测试软件模型的准确性。总共收集了中海油20口井的资料， 其中10口井包含试井表皮系数， 10口井不包含表皮系数。包含表皮系数的井中， 软件计算表皮和实际测试表皮吻合度在90.53%～98.64%之间。不包含表皮系数的井中， 软件对主要损害类型的预测以及严重程度排序符合现场描述。这20口井中， 黏土膨胀储层伤害较为广泛。并且软件中黏土膨胀模块的算法正是采用的本文模型。因此， 从初期应用的效果上说，该模型具有一定的实际意义和较高的准确性， 否则必然会对整体准确度造成不利影响。必须承认的是，这并非最好的验证方式， 但却是目前能做到的最好的验证方式。在未来的工作中， 应该探索更好更直接的矿场实验手段来验证单个模型的准确度。

3 结论

1.模拟了黏土膨胀储层伤害的时空变化特点。其主要特征由2部分组成，①从空间局部看，膨胀吸水作用早期衰减剧烈，中晚期缓慢衰减直至为0，这导致损害量主要由早期贡献，中晚期贡献量少却可以持续很长时间。因此在作业早期采用适当的防膨措施可以控制黏土膨胀储层伤害。②从损害效果看，空间某点的损害率往往不会特别高，但其空间扩展范围较广。这也是为什么有不少井虽然有较明显的黏土膨胀损害，实施强注仍能保证一定的配注量。从储层保护的角度，黏土膨胀应该早预防早控制，才能保证钻完井及注采效果。

2.局限性是考虑了经典的菲克扩散，实际上近年来反常扩散领域的研究表明，水分子在地层多孔固相中的扩散行为很可能会偏离菲克扩散，呈现反常扩散[20]。不过限于篇幅和目的，没有深入探讨该扩散机制，在将来的研究中，可以将反常扩散纳入考虑之中，可以预计的是，在那种情形下其数学结构要复杂得多。

3.局限性之二是该黏土膨胀本质上模拟的是渗透水化，也就是说该模型中处理边界条件z=0的隐藏假设是固-液界面膜厚度为0。近年来有些研究表明，所谓表面水化也可能较重要，这意味着界面膜厚度不为0，而是由界面物理化学性质决定的一个有限值。这个有限膜厚度在建立的模型中会造成较大数学上的困难，因为膜内部的膨胀性质以及其对水吸的额外影响难以量化。所以，这里采用了这一简化，在一定条件下足以对矿场实际产生指导意义。

4.此类动力学模型是储层伤害模拟和诊断研究的趋势，未来的储层伤害模型应该能给出伤害的时空演变规律，模拟储层伤害的动态变化规律，而目前多数模型都不能给出伤害的场特征和时间性质。