碳酸盐岩多尺度三维酸蚀蚓孔立体延伸动态模拟

薛 衡,赵立强,刘平礼,崔明月,蒋卫东,梁 冲,叶颉枭,徐 彬

(1.西南石油大学 油气藏地质及开发工程国家重点实验室，四川 成都 610500； 2.中国石油 勘探开发研究院 廊坊分院，河北 廊坊 065007； 3.中国石油 西南油气田分公司 工程技术研究院，四川 成都 610017)



碳酸盐岩多尺度三维酸蚀蚓孔立体延伸动态模拟

薛 衡1,赵立强1,刘平礼1,崔明月2,蒋卫东2,梁 冲2,叶颉枭3,徐 彬1

(1.西南石油大学 油气藏地质及开发工程国家重点实验室，四川 成都 610500； 2.中国石油 勘探开发研究院 廊坊分院，河北 廊坊 065007； 3.中国石油 西南油气田分公司 工程技术研究院，四川 成都 610017)

在双重尺度模型基础上建立了多尺度三维酸蚀蚓孔延伸模型并阐述了其数值求解方法。同时，利用地质建模方法建立了孔隙空间关联分布模型，成功解决了孔隙空间突变问题，使模拟结果更加真实可靠。通过模拟发现，注入速度是影响蚓孔形态的主要因素,随着注入速度的增大，溶蚀结构可分为面溶蚀、锥形结构、主蚓孔、分支结构和均匀溶蚀5类。当扩散速度与对流速度相当时会产生主蚓孔通道，越靠近注入端口的蚓孔越粗，而越接近出口端的蚓孔越细。主蚓孔在延伸过程中大致会经历以下4个阶段：蚓孔竞争发育、形成优势通道、主蚓孔突破岩心和蚓孔内径扩宽。其次，孔隙空间联通性是蚓孔扩展路径的决定因素，蚓孔最终沿初始高孔渗分布轨迹突破岩心。最后，通过现场算例验证了该模型及方法用于模拟碳酸盐岩储层酸化、酸压中蚓孔扩展的可行性与正确性，并为该研究领域的精细化模拟预测提供了新的手段。

孔隙介质；酸蚀蚓孔；三维模拟；数值模型；碳酸盐岩；提高采收率

大量实验研究表明蚓孔根据溶蚀形态的不同可以分为五种类型：端面溶蚀、锥形溶蚀、主蚓孔、分支型蚓孔和均匀溶蚀[1]。主蚓孔溶蚀对应酸化效率最高，此时用最少的酸液便可以形成有效穿透伤害带的流动通道，达到增产最大化的目的[2]。然而，蚓孔的形成包括一系列复杂的物理化学反应，其中储层压力、温度和物性等不可控因素都会对该物理化学行为产生直接影响。因此在碳酸盐岩储层酸化改造中不仅要研究可控参数(排量与酸液的最优组合等)对蚓孔延伸的影响。同时，应兼顾考虑各种不可控参数(地层应力、温度和非均质性等)对蚓孔延伸的影响。McDuff[3]等人开展了大尺度径向渗流条件下的酸蚀蚓孔实验研究，与以往线性流动岩心实验不同，该实验研究更符合现场实际，但实验成功率较低。由于实验往往更加注重机理研究，通过实验获得的参数往往很难甚至不能直接应用于现场。因此，为了解决实际工程问题，多种蚓孔预测数学模型相继提出。其中，Buijse[4]和Gdanski[5]提出的管道模型假设管状蚓孔已存在，这与实际情况符合程度有限。而Fredd和Fogler[6]、Hoefner和Fogler[7]提出的网络模型，其模拟结果与实验结果往往差异较大。与前两种数学模型相比，Panga[8]等提出的双重尺度模型更好地描述了线性流条件下酸-岩反应动力学及传质控制等机理。2007年，Kalia和Balakotaiah[9]在Panga[8]等基础上推导了更符合现场应用的径向流(r-θ方向)双重尺度模型。

双重尺度模型在碳酸盐岩储层增产改造研究领域中得到了广泛运用与延伸。例如Liu P[10]等, Liu M[11]等，Ratnakar[12]等和Bulgakova[13]等在其基础上开展了不同变粘酸体系的转向酸化规律研究；Ziauddin[14]等，Kalia[9,15]等，Liu M[16-17]等和Izgec[18]等在其基础上研究了溶洞及裂缝型碳酸盐岩介质中的酸-岩流动反应机理。然而，国内外研究目前主要集中在孔隙型介质的二维模拟上，全三维模拟鲜有报道。由于现有模型未考虑各物性参数在三维空间中的变化规律，使得其模拟结果与实际结果存在较大误差。比如蚓孔突破体积的模拟结果往往比实验结果大，蚓孔的模拟形态与实际形态差异显著等。因此本文重点研究了多尺度三维蚓孔扩展模型建模，并阐述了其数值求解方法；最后开展了三维酸蚀蚓孔立体延伸规律研究和现场算例优化。

1 多尺度三维蚓孔扩展模型建模

根据Panga[8]等提出的双重尺度模型，推导出多尺度三维蚓孔扩展模型，假设条件如下：①酸化处理对象为孔隙介质，忽略裂缝及孔洞；②孔隙度的非均质性具有空间相关性；③岩石矿物成分为单一灰岩介质。

1.1 用于描述酸-岩流动反应的达西尺度方程

描述流体渗流的连续性方程

(1)

(2)

式中：U为速度，m/s；K为岩心渗透率，10-3μm2；μeff为有效粘度，mPa·s；p为压力，MPa；Φ为孔隙度，%；t为时间，s；为梯度算符。

描述酸液传质、扩散、反应的物质平衡方程

(3)

描述孔隙内部H+传质、反应的局部平衡方程

(4)

描述孔隙度变化的瞬态方程

(5)

式中：Cf为孔隙内部酸液质量浓度，kmol/m3；Cs为孔隙壁面酸液质量浓度，kmol/m3；Dei为酸液在i方向的有效扩散系数，m2/s；i为下标，x，y和z；kc为传质速度，m/s；R(Cs)为单步不可逆反应的溶蚀速度，R(Cs)=ksCs，m·kmol/(s·m3)；ks为反应速度，m/s；α为酸液的溶蚀能力，kg/kmol；ρs为岩石密度，kg/m3；av为比表面积，m2/m3。

1.2 用于描述岩石物性变化的微观尺度方程

岩石溶解导致孔隙介质中的各物理性质发生改变，Garman-Kozeny模型给出了孔隙度、渗透率、比表面和孔道半径间的微观变化关系。

(6)

式中：Φ0为岩心原始孔隙度，%；K0为岩心原始渗透率，10-3μm2；β为由实验获取的指数，无因次；rp为孔隙半径，m；rp 0为原始孔隙半径，m；a0为原始比表面积，m2/m3。

Panga[8]等给出了传质速度kc与舍伍德数间的关系公式和不同方向上的扩散系数方程。

(7)

(8)

式中：Sh为舍伍德数，无因次；Sh∞为渐进舍伍德数，无因次；m为孔隙长度与孔隙直径的比值，无因次；Rep为孔隙尺度的雷诺数，无因次；Sc为施密特数，无因次；Dm为分子有效扩散系数，m2/s；λi为取决于孔隙结构的常数，无因次；αos为取决于孔隙结构的常数，无因次(αos=1.0表示球体填充床)。

1.3 初始条件和边界条件

t=0时的初始条件

(9)

t>0时的压力边界

(10)

t>0时的浓度边界

(11)

2 数值模型求解方法

首先我们利用地质建模软件建立了孔隙度具有空间关联性的三维随机分布模型(图1)。x，y和z方向上的网格数量分别为60×30×30，尺寸为5 cm×2.5 cm×2.5 cm，岩心平均孔隙度和平均渗透率分别为0.2×10-3μm2和10×10-3μm2。由于全三维精细化的流化耦合数值求解是非常复杂的，计算所需内存及耗时非常大，我们利用隐式差分法离散了公式(2)和公式(3)，并采用优化并行算法成功进行编程求解。

公式(2)离散后的连续性方程

(12)

其中

公式(3)离散后的反应方程

(13)

其中，方程组(13)中的系数此处从略。

每个时间步长中空间各节点的物理性质(如压力、浓度、孔隙度、渗透率等)需要结合初始条件和边界条件[公式(9)—(11)]进行迭代求解，其详细的求解步骤如下：①通过地质建模软件建立初始孔隙空间分布模型，并进行网格划分；②通过公式(1)-(2)计算n+1时刻的压力场和流速场; ③通过公式(7)-(8)和(14)计算当前时刻各空间节点处的反应速度，传质速度和H+扩散速度；④通过公式(3)-(4)计算n+1时刻的酸浓度场；⑤通过公式(5)-(6)依次计算出n+1时刻的孔隙度、渗透率、比表面、孔喉半径，并更新各空间节点的属性；⑥判断：如果蚓孔突破岩心则计算完毕，否则重复步骤②。

图1 空间关联随机孔隙分布模型Fig.1 Spatial correlation porosity distribution model

3 三维酸蚀蚓孔立体延伸模拟及分析

温度和酸液浓度是影响酸-岩反应速度的重要因素，用伊拉克米桑油田的碳酸盐岩(灰岩成分>92%)进行了动力学实验研究，并通过实验数据拟合得到了反应速度与浓度、温度间的关系公式：

(14)

式中：R是气体常数，8.314 J/(mol·K)；T是环境温度，K。其他重要物性参数和注入参数见表 1。

以图1给出的三维孔隙空间分布模型作为初始条件，将表1中的物性参数和注入参数代入数学模型中进行计算，当入口端注入压力与初始注入压力的比值(下文统称为IPR(Injection Pressure Ratio))降为1%时停止计算[8]。图2是不同时刻岩心内部蚓孔空间结构。在驱替作用下，流体在会优先进入高孔隙度区域，因此注酸初期在岩心端面附近先形成了竞争发育的蚓孔(图2a)。当注酸90 min时，岩心中逐渐形成优势通道(图2b)。注酸中后期，大部分酸液进入优势通道内使其持续增长，最终形成穿透整个岩心的主蚓孔通道(图2c)。另外，越靠近注入端口的蚓孔越粗，而越接近出口端的蚓孔越细，此时最大主蚓孔直径约为0.66 cm，最小主蚓孔直径为0.25 cm。持续注酸直到IPR等于0.01，主蚓孔通道进一步变粗，而其它部位孔隙结构无显著变化(图2d)，这是因为在蚓孔穿透岩心后形成无阻流动通道，大部分鲜酸都将流经主蚓孔通道并持续对蚓孔壁面进行冲刷刻蚀，最终使蚓孔直径扩宽。

表1 物性参数和注入参数

Omer[18]等和姚奕明[19]等均发现蚓孔趋于沿大孔道分布的地方扩展。而通过三维蚓孔模拟分析，我们进一步验证了该结论的正确性(图3)。图3a是部分孔隙初始分布剖面，这里虚线圈出了孔隙度较大的位置，并用于对比蚓孔的扩展路径。由于酸液在驱替作用下会优先进入渗流阻力较小的孔隙，因此这些大孔隙分布的位置会更先被溶蚀，最终形成图 3b中的孔隙空间分布剖面。

Damkohler数是描述扩散速度与对流速度关系的变量，也是用于判断酸岩反应溶蚀形态的重要参数。

图2 注酸过程中三维蚓孔形态Fig.2 Wormhole patterns in different stages of acidizing a.45 min,IPR=0.75;b.90 min,IPR=0.34;c.145 min,IPR=0.02;d.299 min,IPR=0.01

图3 酸化前后沿x方向孔隙分布剖面对比Fig.3 Comparison of porosity profiles in x coordinate before and after acidizinga.酸化前;b. 酸化后

图4 不同排量下的孔隙空间溶蚀结构和注入速度与突破孔隙体积(PVBT)关系曲线Fig.4 Wormhole patterns at different injection rates and curve of PVBT vs. injection rate a.Qinj=0.02 mL/min; b.Qinj=0.07 mL/min;c.Qinj=0.2 mL/min;d.Qinj=5 mL/min;e.Qinj=50 mL/min;f.注入速度与突破孔隙体积关系曲线

在低注入排量下，由于对流速度较慢使得Damkohler数非常大，扩散速度在整个反应过程起主导，此时会形成面溶蚀结构(图4a)。适当提高注入排量，此时扩散速度略微占主导，则会形成图4b的锥形溶蚀结构。当扩散速度与对流速度之间达到相当时，则会形成主蚓孔通道，此刻的酸化效率最高(图4c)。而在大排量下，由于对流速度相对更快，岩心更容易形成多分枝结构或均匀溶蚀结构(图4d,e)。而在锥形溶蚀结构、主蚓孔结构和多分枝结构下蚓孔延伸轨迹更易受到孔隙空间结构的影响。从图4f可以看出，该岩心最优注入排量约为0.5 mL/min，对应突破孔隙体积为3.9 PVBT。

4 现场算例分析

以1口垂直井进行单层酸化改造为例，目标层段厚度为10 m,储层初始平均渗透率为10×10-3μm2,初始孔隙度为0.2,井筒半径为0.12 m,设计解堵半径1.5 m。其他参数取典型值，见表1。其中C0=4.4 kmol/m3对应浓度为15%的盐酸体系。为了便于现场实施酸化设计方案，模拟排量范围取0.1～5 m3/min。以排量为横坐标，酸液突破解堵半径(rd=1.5 m)时对应的每米用酸量为纵坐标绘制关系曲线(图5)。

从图5可见，Qinj为0.5 m3/min时对应用酸强度最小，为2.45 m3/m。同时，不论增大还是降低排量，都将使用酸强度增加，从而增加施工成本。其次，排量在0.3～1.5 m3/min内，排量大小对用酸强度影响程度有限。因此在情况许可条件下应尽量提高排量(但不应超过1.5 m3/min)，从而减少作业时间，降低施工安全风险。图6是最优排量下的酸蚀蚓孔扩展形态，从图中可见，井筒附近形成了多条具有高导流能力的主蚓孔通道，局部位置蚓孔长度达1.5 m，有效突破伤害半径，实现了碳酸盐岩储层酸化解堵的目的。

图5 排量与用酸强度关系曲线Fig.5 Injection strength vs. injection rate

图6 最优注入排量下的近井带三维溶蚀结构Fig.6 Wormholes near wellbore under optimal injection rate

5 结论与建议

本文在双重尺度模型基础上，推导出了多尺度三维蚓孔扩展模型，并系统研究了孔隙联通性、注酸速度和酸液用量对酸蚀蚓孔形态的影响。通过现场实例计算表明，该模型及方法能够正确用于指导碳酸盐岩储层的酸化优化设计，减少施工盲目性。对于裂缝型储层的酸压改造涉及的关键技术难题[20](如酸蚀蚓孔发育、蚓孔/天然裂缝滤失、裂缝壁面酸-岩非均匀刻蚀等)均可结合本文相关内容开展进一步研究工作。

[1] Fredd C N,Miller M J.Validation of carbonate matrix stimulation models[J].SPE,2000:39-52.

[2] Huang T,Hill A D,Schechter R S.Reactionrate and fluid loss:The keys to wormhole initiation and propagation in carbonate acidizing[J].Spe Journal,2000,5(3):287-292.

[3] Mcduff D,Jackson S,Shuchart C E,et al.Understanding wormholes in carbonates:Unprecedented experimental scale and 3D visualization[J].Journal of Petroleum Technology,2010,62(10):78-81.

[4] Buijse M A.Understanding wormholing mechanisms can improve acid treatments in carbonate formations[J].Spe Production & Facilities,2000,15(3):168-175.

[5] Gdanski R.A fundamentally new model of acid wormholing in carbonates[C]//Spe European Formation Damage Conference.Society of Petroleum Engineers，1999.

[6] Fredd C N,Fogler H S.Influence of transport and reaction on wormhole formation in porous media[J].Aiche Journal,1998,44(9):1933-1949.

[7] Hoefner M L,Fogler H S.Pore evolution and channel formation during flow and reaction in porous media[J].Aiche Journal,1988,34(1):45-54.

[8] Panga M K R,Ziauddin M,Balakotaiah V.Two-scale continuum model for simulation of wormholes in carbonate acidization[J].Aiche Journal,2005,51(51):3231-3248.

[9] Kalia N,Balakotaiah V.Modeling and analysis of wormhole formation in reactive dissolution of carbonate rocks[J].Chemical Engineering Science,2007,62(4):919-928.

[10] Liu P,Xue H,Zhao L Q,et al.Analysis and simulation of rheological behavior and diverting mechanism of In Situ Self-Diverting acid[J].Journal of Petroleum Science & Engineering,2015,132:39-52.

[11] Liu M,Zhang S,Mou J,et al.Diverting mechanism of viscoelastic surfactant-based self-diverting acid and its simulation[J].Journal of Petroleum Science & Engineering,2013,105(3):91-99.

[12] Ratnakar R R,Kalia N,Balakotaiah V.Modeling,analysis and simulation of wormhole formation in carbonate rocks with in situ cross-linked acids[J].Chemical Engineering Science,2013,90(10):179-199.

[13] Bulgakova G T,Kharisov R Y,Pestrikov A V,et al.Experimental study of a viscoelastic surfactant-based in Situ Self-Diverting acid system:Results and interpretation[J].Energy & Fuels,2013,28(3):1674-1685.

[14] Ziauddin M,Bize E.The effect of pore scale heterogeneities on carbonate stimulation treatments[J].Projectinnovation,2007,44(4):1045-1050.

[15] Kalia N,Balakotaiah V.Effect of medium heterogeneities on reactive dissolution of carbonates[J].Chemical Engineering Science,2009,64(2):376-390.

[16] Liu M,Zhang S,Mou J.Effect of normally distributed porosities on dissolution pattern in carbonate acidizing[J].Journal of Petroleum Science & Engineering,2012,s 94-95(5):28-39.

[17] Liu M,Zhang S,Mou J,et al.Wormhole propagation behavior under reservoir condition in carbonate acidizing[J].Transport in Porous Media,2012,96(1):203-220.

[18] Izgec O,Zhu D,Hill A D.Numerical and experimental investigation of acid wormholing during acidization of vuggy carbonate rocks[J].Journal of Petroleum Science & Engineering,2010,74(1):51-66.

[19] 姚奕明,李双明,蔡鑫,等.孔隙空间分布规律对酸蚀蚓孔扩展规律的影响[J].科学技术与工程,2014 (14):33-37.

Yao Yiming,Li Shuangming,Cai Xin,et al.Effects of spatial distribution of pores on wormhole propagation[J].Science Technology & Engineering,2014 (14):33-37.

[20] 赵立强,刘飞,王佩珊,等.复杂水力裂缝网络延伸规律研究进展[J].石油与天然气地质,2014(4):562-569.

ZhaoLiqiang,Liu Fei,Wang Peishan,et al.A review of creation and propagation of complex hydraulic fracture network[J].Oil & Gas Geology,2014,35(4):562-569.

(编辑 张亚雄)

Dynamic simulation of 3-D multiple-scale wormhole propagation in carbonate rocks

Xue Heng1,Zhao Liqiang1,Liu Pingli1,Cui Mingyue2,Jiang Weidong2,Liang Chong2,Ye Jiexiao3,Xu Bin1

(1.StateKeyLaboratoryofOil&GasReservoirGeologyandExploitationEngineering,SouthwestPetroleumUniversity,Chengdu,Sichuan610500,China;2.LangfangBranchofResearchInstituteofPetroleumExploration&Development,PetroChina,Langfang,Hebei065007,China; 3.EngineeringTechnologyResearchInstituteofSouthwestOil&GasFieldCompany,PetroChina,Chengdu,Sichuan610017,China)

Three dimensional (3-D) simulation of wormhole propagation is crucial for predicting stimulation effects of acidizing or acid fracturing in carbonate reservoirs.This paper presented a 3-D multiple-scale wormhole propagation model built based on the double-scale model,and its solving details.Meanwhile,geological modeling method was used to establish the geology model of porosity correlation distribution in 3-D space,which deals with the saltation of porosity successfully and makes the simulation results more realistic and reliable.Through the simulations,we found that the injection rate is the key factor affecting wormhole patterns.With the increase of the injection rate,five dissolving patterns may be observed,namely face dissolution,conical wormhole,dominant wormhole,ramified wormhole and uniform dissolution.When the convection and dispersion are comparable,the dominant wormhole may form with a thick root and a thin tip.The do-minant wormhole propagation can be divided into the following successive four stages,including competition stage,dominant stage,breaking stage and broadening stage.Most importantly,we found that the wormhole propagation pathway is highly correlated with the spatial correlation of porosity,since the wormhole breaks the core through those large pores with well connectivity.Finally,the field simulation results verified that the models and methods proposed in this paper are feasible and reliable for wormhole simulation of acidizing and acid fracturing in carbonate reservoirs.Those works could provide theoretical support for fine simulation and prediction in the relevant research areas.

porous medium,acidizing wormhole,3D simulation,numerical model,carbonate rock,EOR

0253-9985(2016)05-0792-06

10.11743/ogg20160520

2015-09-09;

2016-08-20。

薛衡(1988—)，男,博士研究生，压裂酸化。E-mail:xuehengbbc@gmail.com。

国家自然科学基金项目(51474182)。

TE357

A