基于水动力不稳定性的油水乳化储层伤害数值模拟

易飞 ， 徐建平 ， 蒋官澄 ，3， 王巧智 ， 苏延辉 ， 高波

0 引言

钻完井过程一般会不可避免地向近井地带引入水基工作液，完井结束之后若不能有效清除这些水相，则滞留在近井地带。一旦开始生产，油相流入井筒之前可能在一定条件下与前述水相形成油水乳化液，乳化液滴堵塞渗流通道形成渗透率损害和表皮。但随着油田含水率的不断提高，单纯采出原油的情况已经很少。对于此种情况，若使用油基钻井液仍需考虑油水乳化伤害的可能性。因此，对乳化损害的研究仍具有现实意义。比如中海油渤海某油田几口井就出现了反相乳化储层损害，造成产液指数下降等问题。一些室内实验研究也表明混合了钻井液滤液和储层流体的乳液对储层可造成较明显的伤害。原油乳化实际上是一种重要的储层伤害[1]。乳化液中所含液滴能堵塞孔隙、增加黏度、降低碳氢化合物的有效流动能力，损害油气层产能。

目前关于乳化的研究主要集中于两个方面，一是对乳化物稳定性的物理化学研究[2-4]，二是研究从分相到乳化的过程所遵循的物理规律[5-8]。其中，第一类研究最多，其本质上是平衡问题，为慢过程，实验或理论研究都可以更细致地展开。而后者本质上是速率问题， 为快过程， 速率问题中动力学因素占主导因素，相关研究更少。但对于乳化储层伤害的模拟而言，后者恰恰更重要，要知道乳化损害的空间分布和随时间变化情况，就必须要模拟油水两相从分相到乳化的动态过程[9-12]。两相界面的破裂准则以及破裂后的两相界面如何形成特定粒径的乳化液滴，都对乳化损害的时空分布定量描述至关重要。

基于两相界面的水动力不稳定性理论，考虑相剪切速率和油水界面张力为乳化损害的力学因素，模拟从完成新井到开始生产过程中，乳化损害在时空的分布及变化规律。当然，正如前述，乳化物稳定性的物理化学特征，特别是油水界面的界面化学性质，是乳化能否形成的主要条件，比如乳化剂（有机处理剂等）的分子结构[13]、HLB值[14]、油相水相组成与特性等内在因素。但主要考虑乳化损害的动力学过程，其主要影响因素是外部力学性质，比如相剪切速率，虽是外部因素，但必不可少。另外，虽然仅用界面张力取代复杂的界面化学性质显然是不够的，但在目前的理论研究中，难以从上述化学性质直接过渡到从分相到乳化液滴的力学行为。而界面张力的取值却直接导致了相相对运动的水动力学特征。换言之， 在界面化学性质允许的条件下，乳化损害发生后，损害在空间的分布和变化规律是什么样的，是无法从理论上基于微观化学结构直接导出的。在这个过程中，上述微观性质对宏观运动的影响只能通过界面张力这样的宏观力学参数来表征。

1 模型建立

1.1 乳化条件和液滴半径

从力学角度看，多孔介质中，油水两相间的低界面张力和高的机械剪切力是引起乳状液形成的主要因素[5]。因此，建立的乳化损害模型必须要能够考虑油水界面张力以及流速导致的机械剪切对乳化液滴形成的影响。笔者在这里的建模主要考虑模拟乳化液滴形成的条件和形成的液滴半径。因为这两者在乳化损害的整个描述中占据重要地位，前者回答是否发生乳化，而后者回答如果发生乳化，损害程度有多大。乳化损害的重要机理即为大量乳化液滴堵塞地层孔隙喉道，造成渗透率损害。

考虑图1所示相界面模型。内相破碎前以某一相对流速穿过外相，其可以是水相也可以是油相。若内相为油相外相为水相形成乳化液；若内相为水相外相为油相则形成反相乳化液。实际中形成乳化液还是反相乳化液，取决于特定的储层与流体性质，也取决于含水率的大小。在此模型中，可以灵活地转换外相和内相以处理这2种情况（见图1）。

图1 两相界面示意图

含水率不同，形成的乳化液形式也不同。以渤海某油田某井为例，该井从2011年完钻后开始生产到2015年，含水率从不到5%逐渐上升到60%～70%，其油水乳化特征动态表现为3个阶段。第1阶段含水率小于反相乳化点，含水率增加，流压和产液下降缓慢或者平稳；第2阶段为反相乳化区，含水率接近反相乳化点，含水率相对稳定，缓慢上升，产液指数下降，产液量最低；第3阶段含水率超过反相乳化点，含水率和产液量快速上升。现场测试反相乳化点在25%～40%之间。

在两相相对运动的过程中，若相界面稳定，则局部的界面扰动将很快衰减，保持界面的稳定性。若相对流速过大或界面张力过低，局部界面扰动被放大，使得整个界面以扰动波长λ1和λ2（分别在x，y方向）被分割为若干小块，这些小块散入外相形成液滴。图1中h为内相液膜厚度。水动力学不稳定性理论的作用在这里就是确定扰动波长λ1和λ2。据文献[5]，上述两波长和一阶水力扰动最高频率～vmax之间满足下述关系：

式中，λ1和λ2分别为方向的扰动波长，μm；v～max为一阶水力扰动最高频率，μm-1。v～max满足下式[5]，

式中，φ为介质孔隙度，μ为油相黏度，k为渗透率，m/s；u为渗流速度，ρ为油相密度，g为重力加速度，σ为油水界面张力。一般有假设λ1≈λ2[5]，由式(1)和式(2)可知：

因此，λ12·h体积的小块形成液滴的半径ro满足关系，由此关系及式(3)可以得到液滴半径如下。

由式(4)可以看出，界面张力越小、相对运动速度越高（剪切越强），形成的液滴半径越小。同时，应该注意到，对于特定的h，可能计算出的r0比h还大，从物理上来说是不可能的。这种情况实际对应着不能发生乳化，等效为r0=0。亦即，若乳化发生，则必有r0＜h且不为0。基于式(4)，可以作出在不同界面张力下r0和渗流速度u之间的关系，如图2所示。所需输入参数值由表1中部分输入数据提供。

图2 乳化液滴半径在不同界面张力下与渗流速度的关系

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

由图2可以看出，在某一特定的界面张力下，只有当渗流速度超过某一临界速度后，才会有乳化液滴产生。并且随着界面张力的增大，这个临界流速也逐渐增大。图2对比了h=6 μm和h=12 μm 2种情况，同等界面张力下，大的液膜厚度可以在更低流速下产生乳化液滴。

一般而言，在假设多孔介质孔隙完全被液体充满的情况下，这里的h就约为地层孔径r。即高渗大孔道的储层较低渗小孔道的储层更容易发生乳化伤害。不添加任何表面活性剂的条件下，油水界面张力约为25 mN/m左右，这时如果相对运动速度只能达到渗流速度级别，按照这里的结果，乳化是不会发生的。但在储层中的油水界面由于各种固相微粒、化学剂及气相充当了乳化剂的作用，使其界面张力更小，那么在油层渗流速度下就有可能发生乳化。同时，若发生乳化，其局限于近井并流动剪切强的区域。

1.2 液滴堵塞模拟

空间不同位置处渗流速度不同，如果产生乳化，其液滴半径也会不同。在确定每个时刻产生的液滴是否堵塞孔隙喉道时，主要考虑尺寸拒止（Size-Exclusion[9]）,即堵塞由液滴和孔喉几何尺寸差异导致。如图3(a)所示。

图3 空间某位置产生乳化液滴堵塞示意图及堵塞概率

在空间某一位置，若某时刻产生了半径为r0（i∆t）的液滴，则此液滴要么运移到下游，要么停止在这里造成渗透率损害。设i时刻产生的液滴r0（i∆t）堵塞空间位置R（储层中某点距离井中心的距离）处的概率为 β（β=β（r0（i∆t）），则自初始时刻i=1至时刻n，空间R处累积孔径损失为无因次渗透率K可由孔径d损失和原孔径的比值度量。r0（i∆t）表示t=i∆t时形成的液滴半径。∆t为时间步长。i为迭代次数。

其中，r 即为地层平均孔径，μm；mK为经验常数，一般取2[10]。要确定β的值，显然知道它是r0的函数：大的液滴堵塞概率大，小的液滴堵塞概率小。其定量的函数关系可由储层孔径分布函数给出，如图3(b)所示。大量研究表明，地层孔径分布函数N(r)近似为log-normal 函数：

式中，μs、σs为分布特征参数，可以看出该函数已经归一化。孔径均值和标准差分别为。通常，通过实际数据能知道的是r 和SD（r）。可以通过这两式反算出μs和σs，如下。

得出分布函数的具体形式。对于某一半径为ro的液滴， 只有分布函数中小于ro的那部分孔径才会被液滴堵塞，因此，堵塞概率β即为0到ro的累积分布：

1.3 乳化损害时空模拟

除了上节所述，某一位置R处随着时间变化产生的乳化液滴有所差异，在同一时刻，不同的位置R处产生的乳化液滴和堵塞概率也是不同的。在取定界面张力的条件下，速度场的变化成为影响乳化损害的关键因素。渗流速度服从达西定律：

其中P为压力场。压力P满足下述的压力传导方程：

其中，Ctotal为岩石流体综合压缩系数，MPa-1；μ为油黏度；φ为孔隙度；k为储层渗透率。不考虑周向变化，梯度算符和拉普拉斯算子可以写为：

综合前述结果，可以建立如下所示的乳化损害模型：

式中，r0（R，i∆ t） 表示在 t=i∆ t时 R 处产生的液滴的半径，单位μm。β（r0）是产生半径为r0的液滴的概率，无单位。Kd（R， t=n∆ t） 表示在 t=i∆ t时R处的无因次渗透率，无单位。

如果模拟的情景为某完钻井开始生产，流体从开始流动到逐渐达到稳定的速度剖面分布。显然，早期近井含油饱和度高，类似2.1中的案例，形成的应是反相乳化。对空间采用中心差分、对时间隐式差分，可以较易地求解此系统。速度场的求解并不复杂，模型的难度在于每个空间点和时间点对液滴半径的判断及历史累积。所需的输入参数总结在表1中。

2 模拟结果

如图4(a)所示，速度场（为负，代表速度向井筒中心）逐渐趋于稳定分布。开井时井壁处压力梯度最大，井壁附近速度大，远处小；随着压降逐渐向储层内传播，井壁处压力梯度逐渐下降，速度下降，而远处速度逐渐增大，直到逐渐达到稳态速度分布。图4(b)为无因次渗透率的时空变化情况。由图4（a）中可以看出，随着曲线逐渐变得更平坦，距离井中心更远的流体也开始参与流动，这个速度场在扩展。由图4（b）可以看出，随着速度场的扩展，符合乳化条件的位置更多，故渗透率损害的前缘也在逐渐扩展。但由于速度场逐渐趋于稳态分布，渗透率损害的前缘不会再向前传播。值得注意的是，渗透率损害最严重的地方并不是最接近井壁的点，而是井壁和前缘之间的某点。这是因为，虽然井壁近旁流体剪切最强，最有利于乳化发生，据前面小节知，在能发生乳化的前提下，剪切越强，产生的液滴半径越小，堵塞概率相应也更小。故最严重的乳化损害发生在损害前缘和井壁之间某个剪切强度适当的位置。但总体上看，乳化损害发生在近井流动剪切强的区域。虽然乳化损害发生的区域很有限，但渗透率损害严重，导致表皮系数较高。当然，若所取的油水界面张力更大，则表皮会更小甚至无法发生乳化以致零表皮。

图4 单井开井生产渗流速度分布和乳化损害后无因次渗透率分布随时间变化情况

如图5所示，随着界面张力逐渐增大，表皮系数逐渐减小，直至没有乳化损害。此结果基于表1的算例，具体参数不同，表皮系数会有差异。

图5 模拟乳化损害表皮系数随时间变化情况

进一步地，其他输入参数，如油相黏度等，也是决定表皮趋势的重要因素。比如2.1中提到的案例，该井油质重，黏度较高，储层条件下估计达几十个毫帕秒，储层油水界面张力按5～25 mN/m估算，则按本文模型计算其导致的乳化损害，表皮系数约为2～6左右，与实际的乳化损害严重程度基本符合。但由于现场表皮测试均为总表皮，不能得到某单一损害因素的表皮系数值，故不可在此精确验证模型的准确度。将这个估算结果同算例对比可知，相比低油黏度（5 mPa·s）在0.5 mN/m的界面张力下（其他参数值基于表1数据）达到6左右的表皮系数，高油黏度下可以在更高的界面张力（5～25 mN/m）达到同等的损害水平，而在这样高的界面张力下基于水动力学不稳定性理论，低黏油体系不发生油水乳化。同理，若两者界面张力相同，则高黏油体系较低黏油体系更易产生乳化损害。这个结论同目前渤海油田出现的乳化损害案例基本是高黏油井相符。另一方面，一般高黏油为胶质、沥青含量较高的体系，含有较多可以充当天然乳化剂的物质，这也导致乳化损害更容易发生。

总体来看，乳化现象非常复杂，无论是形成乳化还是反相乳化，其产生和变化的过程是复杂的化学、物理、力学过程。研究的目的主要是从力学角度着手，研究乳化损害（如果发生）、其空间分布规律及随时间变化的规律。这个模型建立在简化和假设的基础之上，也不可能涵盖所有的乳化机理。

3 结论

1.在仅考虑油水乳化损害的力学特征的前提下，基于水动力不稳定性理论模拟了乳化损害的时空分布规律，结果表明，在考虑界面张力和流动剪切为主控因素时，乳化损害空间范围局限于近井区域，但其造成的渗透率损害在一定条件下却非常可观，导致较大的表皮系数。因此，为了避免给油井生产造成不利影响，在钻完井过程中应减少或避免水基工作液侵入储层，或者对于含水量高的储层避免使用油基钻井液。

2.不同于大部分乳化损害研究，考虑从分相界面到形成乳化的动力学过程，一定程度上更加贴近地描述了乳化储层伤害产生的特点和规律，并且据此实现乳化损害的时空模拟。

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

4.本数学模型的局限性在于主要考虑力学性的乳化形成机理，是乳化产生和演变的外部机理，实际上还有化学性机理作为其内部机理。比如可能考虑乳化剂、表面活性剂、固相微粒等对乳化条件的影响。不过，这类影响都被本模型反映在界面张力的取值大小上，在对界面张力合理的估计下，也不尝为一种简化方法。

[1]WILLIAM E FOXENBERG，SYED A ALI，MINGLIE KE，et al. Eliminate formation damage caused by high density completion fluid - crude oil emulsion[J].SPE 1998，39444.

[2]ALEXEY KABALNOV. The rmodynamic and theoretical aspects of emulsions and their stability[J].Current Opinion in Colloid & Interface Science，1998，3（3）：270-275.

[3]IVAN B IVANOV，PETER A KRALCHEVSKY.Stability of emulsions under equilibrium and dynamic conditions[J].Colloids and Surfaces A Physicochemical and Engineering Aspects，1997，128（1）：155-175.

[4]JEROME BIBETTE，FERNANDO LEAL-CALDERON.Surfactant-stabilized emulsions[J].Current Opinion in Colloid & Interface Science，1996，1（6）：746-751.

[5]R RAGHAVAN，S S MARSDEN JR. Theoretical aspects of emulsification in porous media[J].SPE 1971，308911（2）：153-161.

[6]MONK G W.Viscous energy dissipated during atomisation of a liquid[J].Journal of Applied Physics，1952，23（2）：288-288.

[7]RAJAGOPAL E S. Hydrodynamic aspects of the formation of emulsions (in rheology of emulsions )[M].Oxford, England:Pergamon Press, 1963. 15.

[8]OSTRACH S，KOESTEL A. Film instabilities in twophase flow[J]. AIChE Journal，1965，11（2）：294-303.

[9]ZHENJIANG YOU，ALEXANDER BADALYAN，et al.Size-exclusion colloidal transport in porous media -stochastic modeling and experimental study[J].SPE Journal，2013，162941 ：620-633.

[10]FARUK CIVAN. Non-isothermal permeability impairment by fines migration and deposition in porous media including dispersive transport[J].Transport in Porous Media，2010，85：233-258.

[11]宋兆辉，蓝强.乳化钻井液钻井过程中引起的地层伤害[J].国外油田工程，2008，24（5）：31-36.

SONG Zhaohui，LAN Qiang.Formation damages caused by emulsions during drilling with emulsified drilling fluids[J].Foreign Oil field Engineering， 2008， 24（5）：31-36.

[12]INGEBRET FJELDE. Formation damages caused by emulsions during drilling with emulsified drilling fluids[J]. SPE，2007，105858.

[13]范维玉，赵品晖，康剑翘，等.分子模拟技术在乳化沥青研究中的应用[J].中国石油大学学报(自然科学版 )，2014，38（6）：179-185.

FAN Weiyu，ZHAO Pinhui，KANG Jianqiao，et al.Application of molecular simulation technology to emulsified asphalt study[J].Journal of China University of Petroleum (Edition of Natural Science)，2014，38（6）：179-185.

[14]林吉生，唐培忠，郝清滟，等.稠油HLB值与含水率反相点的关联研究[J].石油炼制与化工，2015，46（8）：61-65.

LIN Jisheng，TANG Peizhong，HAO Qingyan，et al. Correlation of heavy oil HLB value and its moisture inverted point [J]. Petroleum Processing and Petrochemicals，2015，46（8）：61-65.