TOUGH-FLAC3D热流固耦合模拟煤储层水力压裂过程

袁学浩，姚艳斌，甘 泉，刘大锰，周 智[.中国地质大学(北京) 能源学院，北京 0008； .阿伯丁大学 地球科学学院 石油地质系，英国 阿伯丁 AB4 TU；.中国石油 华北油田分公司，河北 任丘 0655]

煤层气藏的开发极大依赖于裂缝对储层渗透率的改造作用[1]。由于我国绝大多数煤层的原位渗透率很低[2-3]，必须采取水力压裂技术，使之在煤储层中形成高渗的导流通道。李同林认为水力压裂在地层中产生裂缝的条件与岩石力学特性、地应力条件、压裂液性质等诸多因素有关[4]。单学军和郝艳丽对煤层压裂的研究发现[5-6]，影响裂缝形态的因素中除了受埋深影响之外，还与构造应力和先存裂缝的发育情况有关。Blanton认为水力裂缝与先存裂缝之间的逼近角和水平主应力差是判断裂缝扩展方向的两个重要参数[7]。

在压裂设计软件方面，目前主要以国外软件为主，如FracproPT，E-StimPlan，Terrfrac和Meyer。FracproPT软件系统是拟三维压裂软件，提供支撑剂和酸化压裂增产的设计、模拟优化功能。FracproPT的独特技术是它的实时数据管理和分析能力。E-StimPlan是全三维压裂设计与分析软件，它具备目前进行压裂优化设计所需的压裂设计诊断和经济优化评价功能，能够完成压裂前地层评估。Meyer是一套拟三维压裂设计分析软件，在压裂充填设计方面有一定的优势。Terrfrac它采用了二维流动方式实现了裂缝扩展的全三维模拟，主要应用地热开发、核废料处理等领域，但是它仅针对水力压裂已知方案模拟，在国外很少见到该软件设计的油田压裂实例[8-9]。

煤储层的水力压裂过程要涉及温度-流体-力学多物理场的耦合，上述软件在处理多场耦合的问题上均稍显不足。多场耦合的问题涉及各因素之间的互相作用，鉴于其多学科交叉理论的复杂性，相关方面的研究才刚刚起步[10]。目前国内外主要是通过数值模拟技术来对其进行描述。COMSOL Multiphy-sics是其中应用较为广泛的一款软件，通过求解偏微分方程组来实现真实物理现象的仿真，最突出特点是具有专业的计算模型库，内置各种常用的物理模型，但是对于模拟储层压裂后形成的不规则裂缝，还有一定局限性。

本次选择将TOUGH与FLAC3D编程搭接来模拟煤层水力压裂。TOUGH为美国劳伦斯伯克利实验室开发的多相流体和反应地球化学模拟程序[11]。FLAC3D是Itasca公司开发的商业软件，擅长进行岩体的热-流-固耦合分析[12]。Zhou L认为FLAC3D的源代码并不适用于水力压裂模拟，他发展了移动边界条件的新算法内嵌入FLAC3D中，模拟结果同拟三维的Fracpro有很好一致性[13]。Rutqvist将TOUGH和FLAC3D两个计算代码搭接起来，模拟了孔隙岩石中的多相渗流，热传导及变形问题[14]。Chen和Jeffrey的模型耦合了流体流动、裂缝扩展以及裂缝表面的摩擦，但未考虑由于温度变化引起的流体性质的改变[15]。Gan Q将TOUGH和FLAC3D搭接，定义了多级流场下的渗透率张量，研究了在多场耦合下的储层压裂过程[16-17]。

本文基于前人在沁水盆地南部煤层气地质研究成果，开展了数值模拟研究工作。研究的特色主要表现在4个方面：①模拟的原始地质模型是根据安泽区块的实际储层地质条件，并结合真实大样三轴物理模拟实验而建立的；②充分考虑了煤储层水力压裂过程中的多场耦合效应，可以体现应力场、流体场与温度场的互相作用与影响过程；③考虑了煤储层的双重孔隙结构，将煤基质与微裂隙分开表征且实现参数的互相传递；④考虑了顶底板的存在，其岩石力学性质与煤储层差异明显，使得模拟结果更加符合实际。

1 模型的建立

1.1 地质模型的确定

地质模型是根据安泽区块的煤储层原位地质条件，并结合真实大样三轴物理模拟实验[18-19](许露露等，2014)建立的。其中物理模拟的天然煤岩样品采自安泽区块周边矿区山西组3#煤层，镜质组反射率分布范围为1.89%～2.78%，属于瘦煤-贫煤。选用不同质量配比的水泥和石英砂来配制成抗压强度为22.2 MPa和31.0 MPa的模拟顶底板，分别代表泥岩和砂岩。制作完成后样品的大小为300 mm×300 mm×300 mm(图1)。根据安泽区块的实际地质参数与物理实验参数设定的模拟方案如表1所示。

1.2 数值模型的建立

根据Gan Q[16-17]的研究，结合安泽地区的实际地质情况，设定了数值模拟参数表(表2)。其中力学参数、煤层与顶底板参数、双重介质参数的取值，均来自于研究区样品实验数据或试井结果，地应力的测试采用凯撒声发射方法。值得指出的是，数值模型中所需参数是井底压力，而施工参数只有压裂液排量。本文根据该区煤层气井的压裂施工资料，将压裂液最大排量转化为井底压力，并取值为40 MPa。

程序耦合时，首先TOUGH计算储层内的水流和热流的传递，输出地层的温度、孔隙度和渗透率等参数；然后通过本构模型把这些参数传递到FLAC3D中进行应力应变计算，通过剪胀角来反映裂缝由于剪切变形导致缝宽的变化，进而修正渗透率[20-21]。煤层中模拟的两条先存裂缝如图2所示，裂缝初始宽度由缝长计算得到[22]。

图1 水力压裂试样示意图Fig.1 Hydraulic fracturing sample sketch

d=0.000 012 5lI

(1)

式中:d为裂缝宽度，m；l为缝长，m；I为指数因子，参考前人研究,I取值为0.8[16]。

2 模拟结果

2.1 逼近角对压力场的影响

图3中1#煤样与2#煤样分别模拟了先存裂缝的逼近角为30°和70°时的压力场分布。可以看出逼近角为30°裂缝比70°裂缝在相同时间内增压范围更大，说明当先存裂缝与水力裂缝之间的逼近角较小时，压裂液更容易沿着先存裂缝突进。另外在等压线展布形态上，30°裂缝在最大主应力方向上压力增加迅速，在模拟第40 s时已波及到边界，而70°裂缝沿最大主应力压力传播受阻，第40 s时压力没传到边界，反而在裂缝的两端形成了压力的相对集中。

图4为1#和2#两个煤样地层流体压力在最大主应力方向(图3中水平方向)上的变化曲线。当逼近角为30°时，地层流体的压力曲线是光滑的，且随着时间增加平缓升高。当逼近角为70°时，地层流压在最大主应力方向上存在明显的传播受阻。图4中2#煤样的裂缝位置压降存在小陡坎，随着时间增加，这种憋压造成的压力前缘会逐渐前移并消减。

2.2 逼近角对温度场的影响

压裂液泵入地层中会引起储层温度场的扰动，其中温度场与渗流场的变化存在相似情况，这是由于水流和热流的运移在TOUGH中二者物理方程形式相同[21]，均满足Fourier定律。在压裂液进入地层的初始阶段(模拟第20 s时)，逼近角为30°裂缝其温度沿裂缝方向出现“近井低、远井高”的现象，降温区域更大；而逼近角为70°裂缝整体温度均一，降温效果不明显(图5)。另外，裂缝温度随模拟时间也呈现下降趋势，逼近角为30°裂缝的温度下降速度明显快于逼近角为70°裂缝。

2.3 岩石力学性质对裂缝扩展的影响

裂缝是否会垂向穿层扩展关键看能否在缝端处产生大于岩石抗张强度的张应力。由于煤岩具有较低的杨氏模量和较高的泊松比，在高压作用下更易变形；同时煤层中广泛发育的微裂隙增加了煤层的导流能力，使得压力在煤层当中的传播快于在顶底板中，因而在3#煤样30 s时的压力分布中，观察到压力在以井为中心呈现环状分布的同时存在垂向上的3层分布。模拟第60 s时观察到在煤层靠近上下隔层的交界面处，形成两个压力集中区域(图6)。那么较大的压力一旦达到顶底板的抗拉强度便会突破顶底板，在此处造成裂缝的垂向穿层扩展。

表1 模拟实验参数设定汇总Table 1 Summary of parameter setting in the simulation experiment

表2 沁水盆地安泽区块数值模拟参数表Table 2 Numerical simulation parameters of the Anze Block in the Qinshui Basin

图2 煤层中的模拟裂缝示意图Fig.2 Sketch map displaying simulated fractures in coalbed

2.4 地应力对裂缝扩展的影响

图7中3#与4#煤样分别模拟深部(930 m)与浅部(650 m)煤层，二者的层内最大主应力等值线图均为椭圆形，长轴沿着水平最大地应力方向，说明了在该方向是裂缝起裂的首选方向。3#煤样由于埋深较深，原地应力较大，煤样中应力分布均匀稳定，煤层中更易发生塑性变形；而4#煤样的原地应力较低，更易发生脆性断裂，其内部应力分布不稳定，有较为明显的错断(图7中白线内)，且易发生于交界面附近。

图3 第40 s时流体压力分布Fig.3 Fluid pressure distribution at the 40th second

图4 流体压力在最大主应力方向上变化曲线Fig.4 Fluid pressure variation curve at σH direction

图5 第20 s时地层温度分布Fig.5 Formation temperature distribution at the 20th second

图6 3#煤样的压力分布模拟结果Fig.6 Pressure distribution simulation results of the 3# coal sample

图7 煤层内最大主应力分布模拟结果Fig.7 Simulation results of the maximum principal stress distribution in coalbed

3 讨论

3.1 裂缝在不同逼近角下扩展规律

周健等采用大尺寸真三轴实验系统，探究先存裂缝对水力裂缝扩展的影响认为二者逼近角θ越大，水力裂缝越容易直接穿过先存裂缝[23]。先存裂缝的临界破裂压力为：

pn=σHcos2θ+σhsin2θ+σt1

(2)

式中:pn为先存裂缝临界破裂压力，MPa；σt1为先存裂缝的抗拉强度，MPa；θ为逼近角， (°)。式(2)经变换可得先存裂缝的临界破裂压力为：

(3)

可知在地应力一定情况下，先存裂缝的临界破裂压力pn会随逼近角增加而增加。30°裂缝逼近角小，临界破裂压力低，压裂液更容易沿此裂缝突进，产生的温度降低也越明显；而70°裂缝逼近角大，临界破裂压力高，压力不断积累容易横穿先存裂缝，在最大主应力方向上形成新的裂缝。

数值模拟的结果与许露露等人[18]的物理实验模拟结果具有很好的一致性(表3)。1#煤样的压裂结果显示水力裂缝只沿先存裂缝扩展，且由于压裂液容易进入先存裂缝而造成压力释放，其顶底板是完好的；而2#煤样压裂后在最大主应力方向上出现新的裂缝，且由于大逼近角下的裂缝的压力封堵作用，使得2#煤样顶板被压穿。

3.2 裂缝垂向扩展穿层规律

赵海峰等[24]认为水力裂缝在垂向上扩展至储层与上下隔层交界面时，可能发生3中情况：①裂缝停止扩展；②裂缝沿地层界面转向扩展；③裂缝直接穿过地层界面进以简化为：

表3 含顶底板煤岩真三轴压裂实验结果Table 3 Results of real tri-axial fracturing experiment with roof,floor and coal rock

(4)

p2=σvcos2θ+σhsin2θ+σt

(5)

(6)

式中:p1为裂缝在界面停止扩展的临界压力，MPa；p2为裂缝沿交界面扩展的临界压力，MPa；p3为裂缝穿透顶底板的临界压力，MPa；K1c为煤层的断裂韧性，MPa·m1/2;K2c为顶底板的断裂韧性，MPa·m1/2；σt为界面抗拉强度，MPa；θ为煤层与顶底板的倾角， (°)。

p1,p2和p3的值用来预测煤层中水力裂缝扩展的模式。如果min(p1,p2,p3)=p1, 则裂缝与界面相交后停止扩展；如果min(p1,p2,p3)=p2, 则裂缝与界面相交后沿地层界面扩展, 裂缝发生转向；如果min(p1,p2,p3)=p3, 则裂缝与界面相交后穿透界面进入顶底板。

虽然3#煤样的抗压强度要大于4#煤样，理应是3#煤样的顶底板不易被压穿，但是结果却显示3#煤样的底板被压穿，而4#煤样的顶底板均保持完整(表3)，这要归因于二者的地应力条件不同。说明相比于岩石力学性质的差异，煤样所处的三轴地应力条件对裂缝是局限在顶底板中还是压穿顶底板具有更为重要的影响。

由于煤层普遍具有明显的应力敏感性[25]，3#煤样由于埋深较大，裂隙较为闭合，导致压裂液的滤失系数较小。因此相比4#煤样，3#煤样更容易在靠近顶底板的交界面处，形成压力相对集中区，当积累的缝内流体压力超过顶底板的临界破裂压力时，顶底板将被压穿。

另一个原因是3#样品具有较大的垂向应力，由于模型中煤层与顶底板之间的倾角为0°，式(5)可以简化为：

p2=σv+σt

(7)

由于煤岩的抗拉强度非常低可以忽略。因此，沿煤岩与顶底板之间交界面扩展的临界压力主要受到垂向应力的影响，由于3#煤样的垂向应力大于4#煤样，因此3#煤样中裂缝沿界面扩展的临界压力大于4#煤样，从而导致3#煤样中流体压力的积累，当流压大于其底板的临界破裂压力后，底板则被压穿。

4 模型验证

图8分别为四块煤样的压裂曲线。1#煤样破裂压力为13.6 MPa，煤层破裂后压力突降后又缓慢上升，最终压裂液排量和滤失达到平衡，压力稳定；2#煤样破裂压力为10.5 MPa，随后压裂液在煤层中均匀渗透，当压裂液渗透到顶板和煤层的胶结面时，压力开始增加，当压力增加到最大值时，顶板开裂。2#煤样大逼近角裂缝对压力起到一定封堵作用，使流压积累，因此在5.2 min就压开了顶板；而1#煤样顶底完好，是由于小逼近角的裂缝在2.6 min就被压开，减缓了流压积累。

3#煤岩在2.2 min起裂，起裂压力不明显，在6.6 min压力突降，说明压裂液进入到下部隔层，压穿底板；4#煤层的破裂压力不明显，随着压裂的进行，在5.6 min压裂液到达外边界，压力变化趋于稳定。

3#煤样由于埋深较大，裂隙闭合，压裂液的滤失系数较小，因此相同时间的泵压值均大于4#煤样。由于3#样品的垂向应力大于4#样品，导致水力裂缝沿着3#煤层与交界面的扩展的临界压力大于裂缝穿透顶底板的临界压力，最终在6.6 min压穿了底板，而4#煤样则顶底完好。

煤样物理实验与压裂曲线可以初步验证将TOUGH与FLAC3D搭接来模拟煤储层水力压裂的准确性，对煤储层双重介质和多场耦合问题考虑之后，可以对压裂过程中渗流场、应力场以及温度场的变化有更为精细的刻画。但模型在反映实际复杂地质问题上还有所欠缺，同时相比其他商业软件，如FracproPT，在提供完整压裂设计方案上还有待提高。

5 结论

1) 先存裂缝的临界破裂压力会随逼近角的增加而增加。逼近角大时，压裂液易沿先存裂缝突进，形成较大的压力场和温度场波及范围；逼近角小时，水力裂缝易横穿先存裂缝扩展，且容易压穿顶底板。

2) 由于煤岩与顶底板力学性质差异显著且煤层中裂隙发育，使得煤层导流能力强，压力在垂向上呈现三层分布；在煤层中靠近上下隔层的交界面处，会形成两个压力集中区，裂缝容易在此处突破顶底板。

3) 煤层中最大主应力等值线为椭圆形，长轴沿着水平最大地应力方向，说明地应力最大的方向是裂缝起裂的首选方向。同时浅部煤层各向地应力较小，更易在交界面处发生脆性变形。

4) 相对于顶底板的力学性质的差异，煤样的原地应力对顶底板是否被压穿作用更为明显。深部煤层由于压裂液滤失系数较小和垂向应力较大造成的流体压力积累，是顶底板被压穿的主要原因。

参 考 文 献

[1] Ferrill D A,Morris A P,Hennings P H,et al.Faulting and fracturing in shale and self-source reservoirs:Introduction[J].AAPG Bulletin,2014,98(11):2161-2164.

[2] Huang B,Liu C,Fu J,et al.Hydraulic fracturing after water pressure control blasting for increased fracturing[J].International Journal of Rock Mechanics & Mining Sciences,2011,48(6):976-983.

[3] Zhang J C.Numerical simulation of hydraulic fracturing coalbed methane reservoir.Fuel 2014; 136:57-61.

[4] 李同林.煤岩层水力压裂造缝机理分析[J].天然气工业,1997,17(tr):53-56.

Li Tonglin.Analysis of coalbed hydraulic fracturing mechanism[J].Natural Gas Industry,1997,17(tr):53-56.

[5] 单学军,张士诚,李安启,等.煤层气井压裂裂缝扩展规律分析[J].天然气工业,2005,25(1):130-132.

Shan Xuejun,Zhang Shicheng,Li Anqi,et al.Analyzing the fracture extended law of hydraulic fracturing in coalbed gas wells[J].Natural Gas Industry,2005,25(1):130-132.

[6] 郝艳丽,王河清,李玉魁.煤层气井压裂施工压力与裂缝形态简析[J].煤田地质与勘探,2001,29(3):20-22.

Hao Yanli,Wang Qinghe,Li Yukui.Simplified analysis of fracture treating pressure and fracture morphology in coalbed gas well[J].Coal Geology & Exploration,2001,29(3):20-22.

[7] Blanton T L.An experimental study of interaction between hydraulically induced and pre-existing fractures.In:SPE 10847,presented at the SPE/DOE unconventional gas recovery symposium,Pittsburgh; May 16- 18,1982.

[8] 计勇,郭大立,王东浩,等.综合应用多种压裂软件技术研究[J].重庆科技学院学报:自然科学版,2009,11(3):135-137.

Ji Yong,Guo Dali,Wang Donghao,et al.Research on application of various fracturing software technology[J].Journal of Chongqing University of Science and Technology,2009,11(3):135-137.

[9] 李勇明,郭建春,赵金洲.油藏整体压裂数值模拟软件的研制与应用[J].断块油气田,2004,11(5):37-39.

Li Yongming,Guo Jianchun,Zhao Jinzhou.Development and application of integral fracturing simulation software[J].Fault-Block Oil & Gas Field,2004,11(5):37-39.

[10] 李士斌,李磊,张立刚.清水压裂多场耦合下裂缝扩展规律数值模拟分析[J].石油化工高等学校学报,2014,27(1):42-47.

Li Shibin,Li Lei,Zhang Ligang.The numerical simulation analysis of crack propagation law under riverfrac treatment multi-field coupling[J].Journal of Petrochemical Universities,2014,27(1):42-47.

[11] Xu T,Sonnenthal E,Spycher N,et al.TOUGHREACT-a simulation program for non-isothermal multiphase reactive geochemical transport in variably saturated geologic media:Applications to geothermal injectivity and CO2 geological sequestration[J].Computers & Geosciences,2006,32(2):145-165.

[12] 于子望,张延军,张庆,等.TOUGHREACT搭接FLAC3D算法[J].吉林大学学报(地球科学版),2013(1):199-206.

Yu Ziwang,Zhang Yanjun,Zhang Qing,et al.Algorithm of TOUGHREACT links to FLAC3D[J].Journal of Jilin University(Earth Science Edition),2013(1):199-206.

[13] Zhou L,Guo J C.3D Modeling of hydraulic fracturing in tight gas reservoirs by using of FLAC3D and validation through comparison with fracPro[J].Advanced Materials Research,2012,524-527:1293-1299.

[14] Rutqvist J,Wu Y S,Tsang C F,et al.A modeling approach for analysis of coupled multiphase fluid flow,heat transfer,and deformation in fractured porous rock[J].International Journal of Rock Mechanics & Mining Sciences,2002,39(4):429-442.

[15] Jeffrey R G.Propped fractured geometry of three hydraulic fractures in Sydney basin coal seams[J],SPE 50061,1998．

[16] Gan Q,Elsworth D.Production optimization in fractured geothermal reservoirs by coupled discrete fracture network modeling[J].Geothermics,2016,62:131-142.

[17] Gan Q,Elsworth D.A continuum model for coupled stress and fluid flow in discrete fracture networks[J].Geomechanics and Geophysics for Geo-Energy and Geo-Resources,2016,2(1):43-61.

[18] 许露露.沁南盆地南部郑庄-安泽区煤储层水力裂缝扩展规律研究[D].北京:中国地质大学(北京),2014:59-77.

Xu Lulu.Mechanism of hydraulic fracture propagation of the coal reservoir in the Zhengzhuang-Anze field,southern Qinshui Basin,China[D].China University of Geosciences,Beijing,2014:59-77.

[19] 许露露,崔金榜,黄赛鹏,等.煤层气储层水力压裂裂缝扩展模型分析及应用[J].煤炭学报,2014,(10):2068-2074.

Xu Lulu,Cui Jingbang,Huang Saipeng,et al.Analysis and application of fracture propagated model by hydraulic fracturing in coal-bed methane reservoir[J].Journal of China Coal Society,2014,(10):2068-2074

[20] 张庆.增强型地热系统人工压裂机理研究及应用[D].吉林大学,2014.

Zhang Qing.Study and application of artificial fracturing mechanism of enhanced geothermal systems[D].Jilin University,2014.

[21] 于子望.多相多组分THCM耦合过程机理研究及其应用[D].吉林大学,2013.

Yu Ziwang.Research on multiphase-multi component THCM coupling mechanism and its application[D].Jilin University,2013.

[22] Olson J E.Sublinear scaling of fracture aperture versus length:An exception or the rule?[J].Journal of Geophysical Research Solid Earth,2003,108(B9):ETG 3-1.

[23] 周健,陈勉,金衍,等.裂缝性储层水力裂缝扩展机理试验研究[J].石油学报,2007,28(5):109-113.

Zhou Jian,Chen Mian,Jin Yan,et al.Experimental study on propagation mechanism of hydraulic fracture in naturally fractured reservoir[J].Acta Petrolei Sinica,2007,28(5):109-113.

[24] Zhao H,Chen M.Extending behavior of hydraulic fracture when reaching formation interface[J].Journal of Petroleum Science & Engineering,2010,74(1):26-30.

[25] Li S,Tang D Z,Pan Z J,et al.Characterization of the stress sensitivity of pores for different rank coals by nuclear magnetic resonance[J].Fuel,2013,111:746-754.