深部煤储层水力压裂改造效应及其能量控因

Transformation effect and energy control mechanism of hydraulic fracturing in deep coal reservoirs

摘要

水力压裂过程涉及能量输入、释放和耗散，探究压裂过程中能量演化机理及其改造效应，可为深部煤储层压裂工艺优化提供理论依据和技术支撑.以鄂尔多斯盆地东缘延川南区块2^#煤层为研究对象，在优化真三轴水力压裂物理模拟方法基础上，结合数值模拟及数学建模，分析了深部煤层水力压裂应力分布及裂隙起裂⁃扩展特征，揭示了输入能量演化机理及其对压裂效果的控制作用.结果表明：钻孔内部填充NaCl粉末、抗拉拔植筋胶加固、试样外侧橡胶垫片缓冲的真三轴水力压裂物理模拟优化方法降低了试样应力加载破碎概率及端部效应，满足深部煤岩水力压裂物理模拟的需要.随排量增大，造缝效果呈现出由小规模、低缝宽复杂缝网向大规模、中缝宽复杂缝网，再向超大规模、高缝宽简单裂缝的转变，排量9 mL/min可形成规模和裂隙复杂度均较高的裂隙体系，破裂压力点是裂隙带扩展由缓至快的转折点.综合考虑地应力和水压力影响，构建了水力压裂应力分布模型，发现环向应力由压应力转化为拉应力是压裂裂隙起裂的驱动力，对于深部正断型应力状态，裂隙不仅形成于垂向应力和最大水平主应力所在的平面，而且优先从平行于垂向应力方向的井筒附近起裂.水力压裂物理模拟裂隙起裂的能量来自于水压增大阶段（阶段Ⅱ）积累的输入能，其由压强能和动能组成，压强能是能量主体.排量与输入能、耗散能呈二次函数关系，与能量输入速率呈对数函数关系，增大排量能够提高能量输入速率，但能量输入效率降低.裂隙率与能量输入速率之间呈高斯分布，当阶段Ⅱ能量输入速率介于0.70~1.15 J/s时，物理模拟压裂改造效果最好.

Abstract

The hydraulic fracturing process involves energy input， release， and dissipation. Investigating the energy evolution mechanisms and stimulation effects during fracturing processes can provide theoretical foundations and technical support for the optimization of fracturing techniques in deep coal reservoirs. Taking the No. 2 coal seam in the southern Yanchuan block in the eastern margin of the Ordos Basin as the research object， by optimizing the physical simulation method of true triaxial hydraulic fracturing， combined with numerical simulation and mathematical modeling， the stress distribution and fracture initiation and propagation characteristics of hydraulic fracturing in deep coalbeds are analyzed， revealing the evolution mechanism of input energy and its control effect on the fracturing effect. The results show the following： The optimization method for the physical simulation of true triaxial hydraulic fracturing， in which the borehole is filled with NaCl powder， reinforced with anti⁃pull⁃resistant glue， and buffered with rubber gaskets on the outer side of the sample， has good reliability， which greatly reduces the probability and end effect of the sample breaking during stress loading and meets the needs of the physical simulation of hydraulic fracturing in deep coal. As the injection rate increases， the fracture creation effect changes from a small⁃scale， low⁃fracture width complex fracture network to a large⁃scale， medium⁃fracture width complex fracture network and then to ultralarge scale， high⁃fracture width simple fractures. A injection rate of 9 mL/min can form a fracture system with a high scale and fracture complexity. The fracture pressure point is the turning point from slow to fast expansion of the fracture zone. Comprehensively considering the effects of in situ stress and water pressure， a model of hydraulic fracturing stress distribution is constructed. The calculation revealed that the transformation of circumferential stress from compressive stress to tensile stress is the driving force for the initiation of fracturing fractures. In the deep normal faulting stress regime， fractures not only form in the plane where the vertical stress and the maximum horizontal principal stress are located but also preferentially initiate near the wellbore parallel to the vertical stress direction. The energy of hydraulic fracturing physical simulation of fracture initiation comes from the accumulated input energy during the water pressure increase stage （stage Ⅱ）. It consists of pressure energy and kinetic energy， and pressure energy is the main body of energy. The injection rate exhibits a quadratic relationship with both the input energy and the dissipated energy， while it follows a logarithmic relationship with the energy input rate. Increasing injection rate can increase the energy input rate， but the energy input efficiency decreases. There is a Gaussian distribution between the fracture rate and the energy input rate. When the energy input rate in Stage Ⅱ is between 0.70 and 1.15 J/s， the effect of physical simulation fracturing is the best.

Graphical abstract

关键词

深部煤层气 / 水力压裂 / 物理模拟 / 数值模拟 / 能量演化 / 改造效应

Key words

deep coalbed methane / hydraulic fracturing / physical simulation / numerical simulation / energy evolution / transformation effect

深部煤层气是我国煤层气产业发展的重要组成部分^［1‑2］.目前，已在鄂尔多斯盆地、准噶尔盆地、四川盆地等地实施了深部煤层气勘探开发工程，并在多个区块取得了重大产气突破，部分区块已实现整体效益开发，使我国煤层气产业发展进入快车道^［3‑5］.但需指出的是，早期深部煤层气勘探开发的突破是在缺乏地质认识和资料支撑的基础上，通过现场压裂试验、技术攻关取得的，深部煤层气地质认识整体落后于现场开发工程的需要^［6］.与此同时，深部新区块开发实践也表明，深部煤储层改造效果存在明显的区块差异性，深部煤储层压裂改造效果的好坏是深部地质条件和工程因素协同作用的结果，在某个区块取得较好改造效果的压裂技术并不能简单地迁移到其他区块.与浅部相比，深部煤层气勘探开发成本更高，只有深化深部煤储层地质理论及储层改造机理的研究认识，才能解决深部煤储层压裂改造过程中的关键科学问题，进而更好地指导现场煤层气压裂工程.

针对深部煤储层水力压裂，前人开展了大量研究工作，发现深部煤储层存在高温、高地应力、高储层压力、低孔隙度、低渗透率、高游离气占比、高饱和度、强塑性、高地层能量等地质特征^［7‑9］，其在压裂时存在破裂压力高、裂隙扩展困难、裂隙易闭合、支撑剂易嵌入等问题^{［1， 10‑11］}.基于这些认识，前人不断创新现场压裂设计理念，主要思路是通过高排量、大规模注入、组合支撑剂及压裂液优化来应对上述问题^{［12‑13］}.现场压裂设计理念的创新取得了较好的应用效果，促使深部煤层气井产气量不断提高，日产气量达1万m³的井大量涌现，部分煤层气井突破10万m³/d^［14］.但是，现场试验的认识更多是基于数据归纳和经验总结，并未从压裂破坏的细观力学机制和内在能量机理方面揭示深部煤层压裂破坏过程中的基础理论问题，这是某一区块效果较好的压裂工程参数在其他区块适用性变差的主要原因.针对现场试验的局限性，前人通过水力压裂物理模拟和数值模拟，研究了压裂破坏的细观力学机制，发现压裂液增压、煤基质软化、水楔效应等是压裂液渗滤及孔隙压力重分布的结果^［15］，裂隙起裂过程受煤岩力学性质、地应力状态、压裂工程参数综合制约^{［16‑17］}.但是，传统物理模拟方法用水泥包裹煤样以制备真三轴水力压裂物理模拟试样^{［18‑19］}，这不仅存在制样周期长、养护过程复杂等不足，水泥还会影响煤样压裂力学性质.此外，数值模拟研究集中于裂缝扩展影响因素敏感性分析、流固耦合模型构建、压裂参数优化等方面^{［20‑21］}，关于深部煤储层水力压裂物理模拟及数值模拟的协同优化研究鲜见.

近年来，随着深部煤储层压裂认识的逐渐提高，前人提出了地质⁃工程精益压裂理念，主要是在保障压裂有效性的同时，降低施工过程中的资源消耗和环境影响，目的是寻找压裂改造效果和开发经济性的最佳结合点，促使深部煤层气产业经济高效发展^［22］.但是，目前缺乏有效表征压裂有效性和开发经济性的量化指标.

水力压裂过程中，压裂液是水压传递的载体，更是能量传递的介质^{［23‑24］}.水力压裂将具有高弹性能的压裂液注入煤层，致使煤岩变形直至破裂，促使煤储层⁃煤层气⁃压裂液耦合的能量系统发生复杂变化.水力压裂过程中能量的输入、演化机理是影响储层改造效果的根本原因，也是量化表征压裂有效性与开发经济性的桥梁和纽带.但前人关于煤储层能量的研究成果主要集中在煤层气成藏过程及机理^{［25‑26］}、煤层气地质选区评价^［27］和煤层气排采过程中排采制度优化^［28］等方面，对水力压裂过程中的能量变化规律及其控制作用的研究鲜有涉及.

为此，本文以鄂尔多斯盆地东南缘延川南区块深部煤储层为研究对象，基于储层能量的研究思路，结合真三轴水力压裂物理模拟、数值模拟及CT扫描等技术手段，以压裂排量、液量、能量为研究切入点，致力于阐明深部煤储层压裂破坏细观力学机制，揭示压裂过程中能量的演化特征及其改造效应，既可以为该区深部煤层气压裂工程实施提供科学依据，也为其他深部煤层气区块压裂参数量化表征提供新的思路.

1 深部煤储层水力压裂物理模拟

1.1 水力压裂物理模拟装置

水力压裂物理模拟装置由计算机、应力控制柜、真三轴水力压裂物理模拟试验机、样品室、超高压恒速恒压泵、中间容器、柱塞计量泵、预加载系统、水压传感器、电热管等组件组成.柱塞计量泵受控制柜控制，用于向试验机的活塞压头泵注液压油.试验机内部的活塞压头可以分别施加和调控最大水平主应力（σ_H）、最小水平主应力（σ_h）和垂向应力（σ_V），三向应力最大可加载至100 MPa.试验机内部布置有电热管，通过计算机程序设定并给样品室加温，加温上限为60 ℃.样品室尺寸灵活，可以开展边长分别为100，200，300 mm立方体试样的压裂物理模拟.超高压恒速恒流泵具有恒压、恒流2种模式，水压上限为150 MPa，最大流量为15 mL/min.中间容器材质为哈氏合金，可实现多种溶剂压裂.计算机连接应力控制柜、水压传感器、电热管和超高压恒速恒压泵，可通过配套软件调控三向应力、温度和恒速恒压泵的注入参数，测试并记录三向应力和实时水压数据（图1）.

1.2 物理模拟试样制备新方法

水力压裂物理模拟试样取自延川南区块华宁煤矿（简称HN）2^#煤层工作面，压裂样品制备流程如图2所示，具体步骤如下：1）选取表面裂隙发育程度低的大块煤样，用石材切割机沿与层理平行/正交方向切割尺寸为100 mm×100 mm×100 mm的立方体块样.2）用磨光机初步打磨煤样表面，并依次用磨料粒径75 μm和10 μm的砂纸将煤样表面打磨光滑.3）用⌀12 mm钻机沿煤样某一平行层理面表面的中心向煤样体心钻长度55 mm的孔，并用铜刷除去钻孔内部的煤粉.4）用螺纹车刀在外径10 mm、内径8 mm的钢管外侧加工长度30 mm的螺纹以增大摩擦力，在内侧加工长度30 mm的螺纹，用于与接头对接.5）在钻孔底部填充高度10 mm的NaCl粉末，向钻孔内部注入抗拉拔植筋胶，用生料带包裹外径10 mm、内径8 mm、长度45 mm的不锈钢管外侧有螺纹、内侧无螺纹的一端的端首，并将该端沿顺时针方向旋转插入钻孔，至不锈钢管完全没入钻孔且螺纹顶端与柱样端面处于同一平面，静置48 h.6）对煤样开展CT扫描，观察植筋胶填充固结情况和NaCl粉末状态，若NaCl粉末均位于钻孔底部且呈细颗粒状、NaCl粉末之间主要为孔隙、钢管外壁与煤样钻孔之间充满植筋胶，表明钻孔加固效果良好.7）用细长钢针刺破不锈钢管端面的生料带，并向钻孔内部注入水以溶解NaCl粉末.

按上述步骤制作水力压裂物理模拟试样，相当于还原了固井过程，并在柱样中心预留了直径10 mm、高度10 mm的裸眼段.CT扫描结果如图2所示，可以看出植筋胶在钻孔和钢管间的填充效果良好，裸眼段NaCl粉末呈颗粒状且富含孔隙，表明试样达到设计要求.

1.3 压裂模拟方案及流程优化

水力压裂物理模拟压裂排量计算方法如下^［29］：

\frac{M_{\mathrm{L}}{M}_{\mathrm{W}}{M}_{\mathrm{Z}}}{M_{\mathrm{Q}}{M}_{\mathrm{T}}}=\mathrm{1},（1）

式中：M_L为裂隙带长度相似比例系数；M_W为裂隙带宽度相似比例系数；M_Z为裂隙带高度相似比例系数；M_Q为压裂排量相似比例系数；M_T为压裂时间相似比例系数.以M_L为例，其计算公式如下：

M_{\mathrm{L}}=\frac{L_{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}}}{L_{\mathrm{f}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}},（2）

式中：L_model为物理模拟试验最大缝长，m；L_field为现场压裂最大缝长，m.

延川南区块常规压裂排量最小为6 m³/min，微地震监测裂缝带长度、宽度、高度依次为250，100，10 m（图3）.水力压裂物理模拟过程中，σ_V>σ_H>σ_h，因此裂缝缝长方向主要平行于σ_V，压裂物理模拟时主裂缝沿轴向贯穿煤样，则缝高为100 mm，按照现场实际压裂的裂缝带长度、宽度、高度之间的比例关系，得物理模拟理论裂缝带长度为2 500 mm，裂缝带宽度为1 000 mm.前人研究发现，水力压裂起裂时间与排量负相关^{［30‑31］}.课题组采集延川南区块煤样开展水力压裂物理模拟参数敏感性分析时发现，起裂时间与排量之间呈负指数关系（图4）.研究区现场压裂数据表明，常规排量压裂起裂时间多介于20~30 min.对于实验室物理模拟，受样品尺寸制约，试样在达到破裂压力后迅速起裂破坏，因此起裂前阶段是本研究的重点.为了对比不同排量压裂时水压曲线的动态变化特征，兼顾试验安全性及设备稳定性，本研究在参数设置时延长物理模拟压裂时间，对M_T赋值为2.

将上述参数代入式（1）和（2），得现场6 m³/min压裂排量对应物理模拟的排量为3 mL/min.此外，延川南区块近年开展了大规模、超大规模水力压裂，开发效果进一步提高.因此，本研究将水力压裂物理模拟压裂排量（Q）设置为3，6，9，12，15 mL/min，分别模拟现场6，12，18，24，30 m³/min排量压裂.研究区2^#煤层埋深介于800~1 600 m之间，平均埋深为1 280 m^［3］，以该煤层平均埋深对应地应力及地温为水力压裂物理模拟应力及温度（T）参数设定的依据，具体为：σ_V=32 MPa，σ_H=26 MPa，σ_h=20 MPa，T=45 ℃.详细压裂方案见表1.

水力压裂物理模拟流程如下：

1）将切割后的立方体煤样置于CT系统载物台上，开展第一次CT扫描，获取样品初始裂隙.

2）按照前文步骤加工试样，并用保鲜膜包裹，置于事先布置好垫块的样品室内.在样品室垫块与煤样接触的6个面之间分别放置尺寸为100 mm×100 mm×5 mm的橡胶垫片，起缓冲作用，防止刚性垫块直接将应力作用在煤上使煤样碎裂.其中，在一个垫片中心钻直径为15 mm的孔洞，用于穿过钢管.装样时，煤样层理面平行于地面，钻孔方向平行于σ_h（图5a）.

3）装样完成后，将样品室吊装至试验机，并将煤样中的钢管与压裂管线相连.封闭试验机，打开电热管，设定T=45 ℃.

4）预加载系统以0.01 MPa/s的速度将三向应力加载至1 MPa，并稳定20 min.

5）设定σ_V=32 MPa，σ_H=26 MPa，σ_h=20 MPa，以0.02 MPa/s恒定应力增速施加三向应力，达到设定应力后稳定30 min，使煤样充分压实.

6）超高压恒速恒压泵以恒流模式向中间容器注入纯净水，驱动中间容器内部的活性水压裂液沿压裂管线进入煤样，通过水压传感器实时监测水压变化（图5b）.本研究压裂液流量分别为3，6，9，12，15 mL/min.

7）达到破裂压力后，继续向煤样内部注入压裂液，直至注入液量与滤失液量达到动态平衡，关闭超高压恒速恒压泵，压裂结束.

8）按照与步骤2）~5）相反的顺序打开试验机，并取出压裂试样，将试样用保鲜膜和透明胶紧紧包裹.

9）对压裂后的煤样再次进行CT扫描，并与初始CT扫描结果进行对比.

1.4 水力压裂物理模拟结果

1.4.1 压裂物理模拟曲线及特征参数

深部真三轴水力压裂物理模拟的水压⁃时间曲线表现出分阶段演化特征：1）自开始压裂至水压由0 MPa转增的临界点为阶段Ⅰ，压裂液充注管线和井筒，未对煤样造成实质影响；2）自水压由0 MPa开始增大的临界点至破裂压力点为阶段Ⅱ，阶段起点表明井筒中已充满压裂液，水压随压裂液注入而增大，阶段终点表明煤样形成宏观裂隙；3）破裂压力点之后为阶段Ⅲ，绝大多数样品水压波动下降，少部分样品在某一压力水平保持相对稳定（图6a）.

由于物理模拟试样尺寸较小，裂缝起裂后迅速扩展至试样边界.因此，阶段Ⅰ、Ⅱ的水压⁃时间曲线具有更好的可对比性.分析注液量、起裂时间、破裂压力等参数发现，不同试样在阶段Ⅰ的注液量均介于79~80 mL，说明管线差异对压裂过程的影响可以忽略不计；阶段Ⅱ注液量介于51.50~88.00 mL，随排量增大，呈先增大后减小趋势，在9 mL/min时达到峰值；起裂前总注液量的变化特征与阶段Ⅱ相似，其差别主要受阶段Ⅱ注液量影响.排量与起裂所需时间呈负指数关系，表明大排量压裂可以大幅缩短压裂起裂时间.破裂压力介于33.66~41.62 MPa，随排量增大，破裂压力先升高后降低，当Q=9 mL/min时达到最大值，这与阶段Ⅱ注液量与排量的关系相似.破裂压力与阶段Ⅱ注液量线性正相关，说明在有限空间内，大液量能够积聚更高的压力（图6b）.综合分析认为，排量、阶段Ⅱ注液量直接影响煤样压裂起裂过程，地应力状态及煤样物质组成的非均质性通过控制煤样物理力学性质，间接影响煤样起裂过程.

1.4.2 压裂物理模拟应力演化特征

在物理模拟基础上，分析水平井压裂时井筒附近应力分布特征.首先，沿井筒径向取截面，该截面在压裂过程中的受力情况如图7所示，在忽略煤岩自重的前提下，煤基质受地应力及水压（σ_w）共同影响，其中，地应力指σ_V及σ_H.为便于构建模型，假设储层是均质、各向同性且为弹性的.

根据力的叠加原理，水力压裂应力分布受地应力引起的静应力和压裂液造成的动应力的叠加影响.前人研究发现，当储层渗透率低于1 mD时，可以忽略压裂液渗流、滤失对流体压力传递规律的影响^［32］.对于地应力，井筒周围静应力分布基于下式求解^［33］：

\left\{\begin{array}{l}{\sigma }_{\rho }=\left(\frac{{\sigma }_{\mathrm{H}}+{\sigma }_{\mathrm{V}}}{\mathrm{2}}\right)\left(\mathrm{1}-\frac{r^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}\right)+\left(\frac{{\sigma }_{\mathrm{H}}-{\sigma }_{\mathrm{V}}}{\mathrm{2}}\right)\left(\mathrm{1}-\frac{r^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}\right)\left(\mathrm{1}-\frac{\mathrm{3}{r}^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}\right)\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{2}\phi ,\\ {\sigma }_{\phi }=\left(\frac{{\sigma }_{\mathrm{H}}+{\sigma }_{\mathrm{V}}}{\mathrm{2}}\right)\left(\mathrm{1}+\frac{r^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}\right)-\left(\frac{{\sigma }_{\mathrm{H}}-{\sigma }_{\mathrm{V}}}{\mathrm{2}}\right)\left(\mathrm{1}+\frac{\mathrm{3}{r}^{\mathrm{4}}}{{\rho }^{\mathrm{4}}}\right)\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{2}\phi ,\\ {\tau }_{\rho \phi }=\left(\frac{{\sigma }_{\mathrm{H}}-{\sigma }_{\mathrm{V}}}{\mathrm{2}}\right)\left(\mathrm{1}-\frac{r^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}\right)\left(\mathrm{1}+\frac{\mathrm{3}{r}^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}\right)\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{2}\phi ,\end{array}\right.（3）

式中：σ_ρ 为径向应力，MPa；σ_φ 为环向应力，MPa；τ_ρφ 为切应力，MPa；r为井筒半径，m；ρ为计算点到井筒中心的距离，m；φ为极角，（°）.

对于σ_w，将煤样简化为外径为R，内径为r的圆筒（图8a），进而求解σ_w影响下应力分量的分布及变化（图8b）.在不考虑地应力的影响下，圆筒外侧不受力的约束，即围压（σ₃）为0 MPa，内部应力随压裂液注入不断增大.应力函数（Φ）以极径（ρ）和极角（φ）为自变量，假设煤样切面均质，则煤样受力状态呈轴对称状，应力函数仅为ρ的函数.此时应力分量表达式可写为

\left\{\begin{array}{l}{\sigma }_{\rho }=({\sigma }_x{)}_{\phi =\mathrm{0}}=\frac{\mathrm{1}}{\rho }\frac{\partial \Phi }{\partial \rho }+\frac{\mathrm{1}}{{\rho }^{\mathrm{2}}}\frac{{\partial }^{\mathrm{2}}\Phi }{\partial {\phi }^{\mathrm{2}}}=\frac{\mathrm{1}}{\rho }\frac{\mathrm{d}\Phi }{\mathrm{d}\rho },\\ {\sigma }_{\phi }=({\sigma }_y{)}_{\phi =\mathrm{0}}=\frac{{\partial }^{\mathrm{2}}\Phi }{\partial {\rho }^{\mathrm{2}}}=\frac{{\mathrm{d}}^{\mathrm{2}}\Phi }{\mathrm{d}{\rho }^{\mathrm{2}}},\\ {\tau }_{\rho \phi }=({\tau }_{xy}{)}_{\phi =\mathrm{0}}=-\frac{\partial }{\partial \rho }·\left(\frac{\mathrm{1}}{\rho }\frac{\partial \Phi }{\partial \phi }\right)=\mathrm{0}.\end{array}\right.（4）

为求解式（4），基于弹性力学基本理论的变形协调方程如下：

{\left(\frac{{\partial }^{\mathrm{2}}}{\partial {\rho }^{\mathrm{2}}}+\frac{\mathrm{1}}{\rho }·\frac{\partial }{\partial \rho }+\frac{\mathrm{1}}{{\rho }^{\mathrm{2}}}·\frac{{\partial }^{\mathrm{2}}}{\partial {\phi }^{\mathrm{2}}}\right)}^{\mathrm{2}}Φ=0.（5）

而轴对称条件下式（5）可简化为

\begin{array}{c}{\mathrm{\Delta }}^{\mathrm{2}}\left(\frac{{\mathrm{d}}^{\mathrm{2}}\Phi }{\mathrm{d}{\rho }^{\mathrm{2}}}+\frac{\mathrm{1}}{\rho }·\frac{\mathrm{d}\Phi }{\mathrm{d}\rho }\right)=\\ \left(\frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{\rho }^{\mathrm{2}}}+\frac{\mathrm{1}}{\rho }·\frac{\mathrm{d}}{\mathrm{d}\rho }\right)\left(\frac{{\mathrm{d}}^{\mathrm{2}}\Phi }{\mathrm{d}{\rho }^{\mathrm{2}}}+\frac{\mathrm{1}}{\rho }·\frac{\mathrm{d}\Phi }{\mathrm{d}\rho }\right)=\mathrm{0},\end{array}（6）

式中：Δ²为拉普拉斯算子，轴对称条件下可写为

{\mathrm{\Delta }}^{\mathrm{2}}=\left(\frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{\rho }^{\mathrm{2}}}+\frac{\mathrm{1}}{\rho }·\frac{\mathrm{d}}{\mathrm{d}\rho }\right)=\frac{\mathrm{1}}{\rho }·\frac{\mathrm{d}}{\mathrm{d}\rho }\left(\rho \frac{\mathrm{d}}{\mathrm{d}\rho }\right).（7）

将式（7）代入式（6），可得轴对称条件下极坐标形式的变形协调方程：

\frac{\mathrm{1}}{\rho }·\frac{\mathrm{d}}{\mathrm{d}\rho }\left\{\rho \frac{\mathrm{d}}{\mathrm{d}\rho }\left[\frac{\mathrm{1}}{\rho }·\frac{\mathrm{d}}{\mathrm{d}\rho }\left(\rho \frac{\mathrm{d}\Phi }{\mathrm{d}\rho }\right)\right]\right\}=\mathrm{0}.（8）

对式（8）积分4次，得到Φ的通解：

\Phi =A\mathrm{l}\mathrm{n}\rho +B{\rho }^{\mathrm{2}}\mathrm{l}\mathrm{n}\rho +C{\rho }^{\mathrm{2}}+D,（9）

式中：A，B，C，D为待定系数.

将式（9）代入式（4），得轴对称条件下应力分量的一般性解答：

\left\{\begin{array}{l}{\sigma }_{\rho }=\frac{A}{{\rho }^{\mathrm{2}}}+B(\mathrm{1}+\mathrm{2}\mathrm{l}\mathrm{n}\rho )+\mathrm{2}C,\\ {\sigma }_{\phi }=-\frac{A}{{\rho }^{\mathrm{2}}}+B(\mathrm{3}+\mathrm{2}\mathrm{l}\mathrm{n}\rho )+\mathrm{2}C,\\ {\tau }_{\rho \phi }=\mathrm{0}.\end{array}\right.（10）

而圆筒内、外表面的应力边界条件为

\left\{\begin{array}{l}({\tau }_{\rho \phi }{)}_{\rho =r}=\mathrm{0},\\ ({\tau }_{\rho \phi }{)}_{\rho =R}=\mathrm{0},\\ ({\sigma }_{\rho }{)}_{\rho =r}={\sigma }_{\mathrm{w}},\\ ({\sigma }_{\rho }{)}_{\rho =R}={\sigma }_{\mathrm{3}}=\mathrm{0}.\end{array}\right.（11）

式中　 σ₃为围压，MPa.

由式（10）可知，式（11）中前2个条件自然满足，后2个条件可联立式（10）、式（11）得到：

\left\{\begin{array}{l}\frac{A}{r^{\mathrm{2}}}+B(\mathrm{1}+\mathrm{2}\mathrm{l}\mathrm{n}r)+\mathrm{2}C={\sigma }_{\mathrm{w}},\\ \frac{A}{R^{\mathrm{2}}}+B(\mathrm{1}+\mathrm{2}\mathrm{l}\mathrm{n}R)+\mathrm{2}C=\mathrm{0}.\end{array}\right.（12）

进一步考察位移单值条件，在轴对称应力状态下，相同ρ值对应的环向位移相等，结合位移分量表达式，并排除多值项，可知B=0.将B=0代入式（12），可知：

\left\{\begin{array}{l}A=\frac{r^{\mathrm{2}}{R}^{\mathrm{2}}{\sigma }_{\mathrm{w}}}{R^{\mathrm{2}}-r^{\mathrm{2}}},\\ C=-\frac{{\sigma }_{\mathrm{w}}{r}^{\mathrm{2}}}{\mathrm{2}(R^{\mathrm{2}}-r^{\mathrm{2}})}.\end{array}\right.（13）

将式（13）代入式（10），整理得：

\left\{\begin{array}{l}{\sigma }_{\rho }=-\frac{(R^{\mathrm{2}}/{\rho }^{\mathrm{2}})-\mathrm{1}}{(R^{\mathrm{2}}/r^{\mathrm{2}})-\mathrm{1}}·{\sigma }_{\mathrm{w}},\\ {\sigma }_{\phi }=\frac{(R^{\mathrm{2}}/{\rho }^{\mathrm{2}})+\mathrm{1}}{(R^{\mathrm{2}}/r^{\mathrm{2}})-\mathrm{1}}·{\sigma }_{\mathrm{w}}.\end{array}\right.（14）

地层条件下，模型外径R趋近于无穷大，即R→∞，可得到具有圆孔的无限大薄板，或具有圆形孔道的无限大弹性体，可将式（14）改写为

\left\{\begin{array}{l}{\sigma }_{\rho }=-\frac{r^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}{\sigma }_{\mathrm{w}},\\ {\sigma }_{\phi }=\frac{r^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}{\sigma }_{\mathrm{w}}.\end{array}\right.（15）

将式（3）与式（15）合并，可得煤基质在地应力和σ_w共同影响下的应力分量表达式：

\left\{\begin{array}{l}{\sigma }_{\rho }=\left(\frac{{\sigma }_{\mathrm{H}}+{\sigma }_{\mathrm{V}}}{\mathrm{2}}\right)\left(\mathrm{1}-\frac{r^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}\right)+\left(\frac{{\sigma }_{\mathrm{H}}-{\sigma }_{\mathrm{V}}}{\mathrm{2}}\right)\left(\mathrm{1}-\frac{r^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}\right)\left(\mathrm{1}-\frac{\mathrm{3}{r}^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}\right)\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{2}\phi +\frac{r^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}{\sigma }_{\mathrm{w}},\\ {\sigma }_{\phi }=\left(\frac{{\sigma }_{\mathrm{H}}+{\sigma }_{\mathrm{V}}}{\mathrm{2}}\right)\left(\mathrm{1}+\frac{r^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}\right)-\left(\frac{{\sigma }_{\mathrm{H}}-{\sigma }_{\mathrm{V}}}{\mathrm{2}}\right)\left(\mathrm{1}+\frac{\mathrm{3}{r}^{\mathrm{4}}}{{\rho }^{\mathrm{4}}}\right)\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{2}\phi -\frac{r^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}{\sigma }_{\mathrm{w}},\\ {\tau }_{\rho \phi }=\left(\frac{{\sigma }_{\mathrm{H}}-{\sigma }_{\mathrm{V}}}{\mathrm{2}}\right)\left(\mathrm{1}-\frac{r^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}\right)\left(\mathrm{1}+\frac{\mathrm{3}{r}^{\mathrm{2}}}{{\rho }^{\mathrm{2}}}\right)\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{2}\phi .\end{array}\right.（16）

基于研究区2^#煤层地应力状态和压裂数据，设定σ_H=26 MPa，σ_V=32 MPa，σ_w最大值为55 MPa，分析σ_ρ 、σ_φ 随σ_w和ρ的变化规律.计算过程中，定义压应力为正、拉应力为负.

x轴正方向，φ=0°，当σ_w=0 MPa时，σ_ρ 、σ_φ 均为压应力，σ_ρ 在近井端小于远井端，近井端为0 MPa，远井端接近σ_H，且随相对半径（r/ρ）增大，σ_ρ 的增速变小；近井端σ_φ 大于远井端，意味着近井端存在应力集中现象，随r/ρ增大，σ_φ 的降速变小，在远井端的最小值接近σ_V.表明自近井端至远井端，靠近井筒部位σ_ρ，σ_φ 变化更为显著.对比距井筒不同位置处σ_ρ，σ_φ 随σ_w的变化，发现σ_ρ 与σ_w线性正相关，近井端增速显著大于远井端，表明σ_w对近井端σ_ρ 的影响最大.σ_φ 与σ_w线性负相关，近井端降速显著大于远井端，表明σ_w对近井端σ_φ 的影响也最大.σ_w增大过程中，井筒附近σ_ρ 迅速增大、σ_φ 迅速减小且有由压转拉趋势，若σ_w继续增大，可能在σ_φ 转为拉应力并超过试样抗拉强度后促使试样发生拉伸破坏（图9a）.

y轴正方向，φ=90°，当σ_w=0 MPa时，随r/ρ增大，σ_ρ 和σ_φ 的变化与φ为0°时趋势一致，仅应力分量的数值存在差异，表现为：与φ=0°相比，σ_ρ 在φ=90°方向远井端的稳定值更大，近井端差异较小；而σ_φ 在近井端和远井端的稳定值更小.随σ_w增大，σ_ρ 始终处于压应力状态，且与σ_w正相关，σ_φ 则由低σ_w时的压应力状态转换为高σ_w时的拉应力状态，这是近井端压裂形成拉张裂隙的关键.r/ρ固定时，σ_ρ 与σ_w线性正相关，σ_φ 则线性负相关，距离井筒越近，斜率越大（图9b）.对比图9a和9b发现，在σ_w增大过程中，y轴正方向σ_ρ 的增大量和σ_φ 的减小量与x轴正方向相等，表明σ_w作用下σ_φ 和σ_ρ 在井筒各方位均匀变化.叠加地应力状态对σ_φ 的影响，σ_φ 沿y轴正方向率先由压转拉.

井筒附近σ_ρ 的分布特征如图10所示，当σ_w=0 MPa时，σ_ρ 在近、远井端均处于压应力状态，由近及远，σ_ρ 逐渐增大，在水平方向趋近于σ_H，垂直方向趋近于σ_V.对比发现，水平方向和垂直方向的σ_ρ 均随σ_w增大而增大，且水平方向始终小于垂直方向，表明地应力状态对σ_ρ 分布的影响显著.

σ_φ 的分布特征如图11所示，当σ_w=0 MPa时，井筒附近存在应力集中现象，且沿σ_H方向更为显著，σ_φ 最大可达70 MPa，接近σ_H的3倍.随σ_w增大，σ_φ 的应力集中程度逐渐减弱.当σ_w超过46 MPa后，σ_φ 率先沿σ_V方向由压转拉，是压裂形成拉张裂隙的优势方位.与σ_ρ 和σ_φ 不同，τ_ρφ 不随水压变化，仅受地应力控制.由式（16）可知，无论是沿σ_H方向还是σ_V方向，τ_ρφ 始终等于0 MPa.

1.4.3 压裂裂隙起裂‑扩展特征

水力压裂物理模拟前后煤样CT扫描及三维重构结果如图12所示，压裂前煤样结构致密，仅发育少量规模较小的显微裂隙，表明不同试样的差异较小，可以尽可能减弱天然裂隙、结构差异等对物理模拟结果的影响，有助于分析排量对压裂裂缝扩展的控制作用.

HN3⁃1的Q=3 mL/min，压裂后沿垂直于σ_h，σ_H和σ_V方向形成了3组宏观裂隙，大多数裂隙未完全贯穿煤样，主要在井筒附近形成复杂缝网，仅在试样上部形成沿天然裂隙扩展的裂隙.HN1⁃1的Q=6 mL/min，压裂后在井筒附近形成了较复杂的裂缝体系，包括2条“X”形相交、方向平行于σ_V的主裂隙.HN1⁃2的Q=9 mL/min，压裂形成了2组互相垂直且贯穿煤样的宏观裂隙，2组主裂隙分别沿平行于σ_V和σ_H方向扩展.与HN3⁃1和HN1⁃1相比，HN1⁃2的裂隙规模更大，且在主裂隙周围形成了大量分支裂隙.HN1⁃4的Q=12 mL/min，与前述3个煤样相比，压裂裂缝较为简单，主裂隙呈“T”形，优势方向垂直或低角度斜交于σ_h.HN2⁃1的Q=15 mL/min，沿垂直于σ_h方向形成了1条贯穿煤样的简单裂缝，裂隙的复杂程度最低，裂隙面平直光滑（图12）.

裂隙重构结果表明，天然裂隙对深部煤储层水力压裂的影响较小，随着排量增大，地应力对人工裂隙方向及复杂程度的影响逐渐增大.在相同地应力状态下，随排量增大，压裂裂缝由小规模复杂裂隙向大规模复杂裂隙和大规模简单裂隙转变，且在Q=9 mL/min时形成规模和裂隙复杂度均较高的裂隙体系.结合前文应力分量演化特征，可知对于深部正断型应力状态，人工裂隙不仅优先形成于σ_H和σ_V所在的平面，且首先沿平行于σ_V方向起裂.

在三维重构基础上，计算煤样不同切片的裂隙面积分数和裂隙率，发现所有煤样均在井筒附近达到裂隙面积分数的最大值，表明井筒附近裂隙最复杂，与橡胶垫片接触部位裂隙发育程度低，垫片起到了保护煤样的效果.随排量增大，裂隙率先增大后减小，一方面说明过小的压裂排量会限制压裂裂隙向远端扩展，影响储层改造效果；另一方面说明压裂规模并非越大越好，对于特定埋深和地质条件下的煤储层，存在最优排量.在本研究中，当排量介于6~12 mL/min之间时，裂隙率较高，且在Q=9 mL/min时达到最大值，是给定应力、温度条件下的最优排量（图13）.

随排量增大，裂隙率先增大后减小的地质控制机理在于：压裂排量较小时，缝内净压力较低，压裂液优先进入近井地带的结构弱面，并通过滤失进入煤基质.由于压裂液流速较低，所携能量不足以促使裂缝网络向远井端扩展，最终形成局限于近井地带的复杂裂缝网络.与整个煤样相比，裂缝网络的规模较小，因此整个煤样的裂隙率较低.随着排量增大，压裂液携带能量大幅增加，缝内净压力显著提高，驱动人工裂隙沿多个方向向远井端扩展，在此过程中，部分天然裂隙被捕获并得以扩展，最终形成大范围裂缝网络.与整个煤样相比，裂缝网络规模显著增大，因此整个煤样的裂隙率提高.当排量进一步增大，缝内净压力极高，压裂液携带的巨大能量在短时间内沿着阻力最小的方向扩展，天然裂隙未被有效捕获、扩展，最终形成平行于σ_V的单一主裂缝.由于缺乏分支缝，因此超大排量压裂时的裂隙率反而最小.

与文献［34］的CT扫描结果相比，本研究煤样的裂隙率较低，原因在于水力压裂物理模拟煤样的尺寸较大，CT扫描的最低分辨率是67.69 μm.因此，本文计算裂隙率中不包括尺度低于67.69 μm的孔隙和极小的微裂隙.由于水力压裂形成裂隙的尺寸往往大于天然微裂隙，所以CT扫描分辨率不会对裂隙率变化的整体趋势造成太大影响.

研究还发现，在本研究温压条件下，无论排量大小，均在钻孔底部形成密集的微裂隙，附近煤基质呈层状或团块状从钻孔内壁剥离，本文将此区域称为微裂隙密集区；当水压达到破裂压力后，宏观裂隙在钻孔底部微裂隙密集区形成，并向试样外侧扩展，将此区域称为主裂隙扩展区（图14）.分析认为，微裂隙密集区的成因包括以下三方面：一是煤的抗拉强度低，σ_φ 随σ_w增大而增大，当σ_φ 达到自由面附近煤基质抗拉强度后，煤基质开裂，并在基质单元差异应变作用下剥落；二是煤具有强非均质性，随σ_w增大，煤基质沿径向发生差异性位移，导致钻孔附近剪应力集中，当其大于煤基质单元抗剪强度时，自由面附近煤岩隆起、开裂，进而沿径向发展为层状剥落；三是煤在高温、高应力条件下呈弹塑性，近井端塑性区能耗增大，压裂液在煤塑性变形的同时渗入塑性区基质，改变孔隙压力分布，导致微裂隙宽度增大但延伸遇阻.在达到破裂压力之前，井筒附近优先形成密集的微裂缝，而后在水压、地应力、密集微裂隙、层理等共同影响下形成宏观裂隙.

对比不同排量压裂效果，发现与小排量压裂相比，大排量压裂在缩短起裂时间的同时，有效避免了压裂液在近井附近复杂破碎带的滤失，使更多的压裂液向远井端流动，可增大压裂改造范围，形成大规模裂隙.但当压裂排量过大时，容易形成延伸距离远、难以对煤岩整体增渗的简单裂隙.总体而言，在σ_V=32 MPa、σ_H=26 MPa、σ_h=20 MPa的深部正断型应力状态下，随排量增大，煤样压裂改造效果并非越来越好，而是先变好后变差.当排量低于9 mL/min时，往往沿不同方向形成多条主裂隙，其受地应力状态的影响较小；当排量大于9 mL/min时，仅形成1~2条主裂隙，且主裂隙的方向受地应力状态的影响显著；在本试验模拟研究中，当压裂排量为9 mL/min（对应现场18 m³/min）时，取得的压裂改造效果最好.

2 深部煤储层水力压裂数值模拟

真三轴水力压裂物理模拟具有结果可靠、直观等优点，但也存在成本高昂、可重复性差的局限性.前人研究已经证实，数值模拟同样是一种切实可行的研究手段.真实破裂过程分析（RFPA）兼具有限元和统计损伤理论的优点，适于解决岩石损伤破坏问题^［35］.本文用RFPA^3D⁃Flow版开展水力压裂数值模拟.模拟过程参照水力压裂物理模拟试验参数，以压裂排量、液量及改造效果为研究对象，研究裂隙起裂⁃扩展过程，并与物理模拟结果对比验证.

2.1 水力压裂数值模拟方法

RFPA基本原理可概括为以下3条：1）将复杂宏观非线性问题转化为简单细观线性问题；2）将复杂非连续介质问题转化为简单连续介质力学问题；3）通过不同类型细观单元变形、破坏的积累反映宏观行为的演化^［36］.

对于每个计算单元，当单元的应力、应变达到设定阈值后，该单元开始累计损伤，以弹性模量（E）为例，其损伤累计过程可用下式^［37］描述：

E_{\mathrm{I}}=(\mathrm{1}-D_{\mathrm{m}})E,（17）

式中：E_I为损伤后的弹性模量，GPa；D_m为损伤变量.

对于拉伸状态，D_m计算公式如下：

（18）

式中：η为单元残余抗拉强度系数；ε为拉伸应变，%；ε_t0为单元初始损伤拉伸应变，%；ε_tu 为单元极限拉伸应变，%.

对于压缩过程，采用Mohr⁃Coulomb准则判别单元所处应力状态为抗压或剪切状态：

{\mathrm{\sigma }}_{\mathrm{1}}-{\mathrm{\sigma }}_{\mathrm{3}}\frac{\mathrm{1}+\mathrm{s}\mathrm{i}\mathrm{n}\phi }{\mathrm{1}-\mathrm{s}\mathrm{i}\mathrm{n}\phi }\ge {\sigma }_{\mathrm{c}},（19）

式中：σ₁为轴压，MPa；σ_c为抗压强度，MPa.

此时，D_m可由下式计算：

（20）

式中：ξ为单元残余抗压强度系数；ε_c0为单元极限压缩应变，%.

煤在细观尺度上主要发生拉伸破坏，其次为剪切破坏.因此，RFPA将拉伸破坏准则作为第一损伤阈值准则，Mohr⁃Coulomb准则作为第二损伤阈值准则^［36］.

数值模拟过程中，认为单元尺度上的水力特征呈各向同性，基于多孔介质Darcy定律，流体流动控制方程^［38］如下：

\mathrm{\Delta }\left[\frac{k{\rho }_{\mathrm{a}\mathrm{c}}g}{\mu }(\mathrm{\Delta }{P}_{\mathrm{w}}-{\rho }_{\mathrm{a}\mathrm{c}}g\mathrm{\Delta }z)\right]=S\frac{P_{\mathrm{w}}}{\partial t}-F\frac{\partial {\epsilon }_{\mathrm{v}}}{\partial t},（21）

式中：k为渗透系数，m/s；ρ_ac为流体密度，kg/m³；g为重力加速度，m/s²；μ为流体黏度，Pa·s；P_w为流体压强，MPa；z为高度，m；S为储水率，m^-1；F为孔隙水压力系数；t为时间，s；ε_v为单元体的体应变，%.

单元的变形控制方程^［37］为

G{u}_{i,j}+\frac{G}{\mathrm{1}-\mathrm{2}\nu }{u}_{j,i}-F{P}_w+f_i=\mathrm{0},（22）

式中：G为剪切模量，MPa；u_i，j 和u_j，i 为位移分量，%；ν为泊松比；f_i 为体积力分量，MPa.

完整煤岩和损伤煤岩的应力⁃渗流耦合方程分别见式（23）和（24）^［37］：

K_{\mathrm{e}}=K_{\mathrm{0}}{\mathrm{e}}^{-\chi \left(\frac{{\sigma }_{\mathrm{V}}+{\sigma }_{\mathrm{H}}+{\sigma }_{\mathrm{h}}}{\mathrm{3}}-\alpha P_{\mathrm{w}}\right)},（23）

K_{\mathrm{d}}=\frac{\sqrt[\mathrm{3}]{V_{\mathrm{i}}^{\mathrm{2}}}{\rho }_{\mathrm{a}\mathrm{c}}g}{\mathrm{108}\mu }{\epsilon }_{\mathrm{v}}^{\mathrm{2}},（24）

式中：K_e为未损伤单元的渗透率，m²；K₀为单元初始渗透率，m²；χ为耦合系数；α为Biot系数；K_d为损伤单元的渗透率，m²；V_i为单元体积，m³.

2.2 数值模拟参数标定

RFPA水力压裂数值模拟用到的参数包括平均细观弹性模量、平均细观抗压强度、压拉比、泊松比、内摩擦角、均质度系数、残余强度系数、初始渗透率、孔隙度、地应力、压裂参数等.其中，压拉比由煤样抗压、抗拉强度计算得到，抗压强度测试参照GB/T 23561.7—2009《煤和岩石物理力学性质测定方法 第7部分：单轴抗压强度测定及软化系数计算方法》及GB/T 23561.9—2009《煤和岩石物理力学性质测定方法 第9部分：煤和岩石三轴强度及变形参数测定方法》；抗拉强度测试参照GB/T 23561.10—2010《煤和岩石物理力学性质测定方法 第10部分：煤和岩石抗拉强度测定方法》，测试得到煤样压拉比为13.30，三轴压缩内摩擦角为19.88°；初始渗透率、孔隙度测试参照SY/T 6385—2016《覆压下岩石孔隙度和渗透率的测定方法》，测试结果分别为2.3×10^⁃16 m²，4.52%.此外，部分参数为煤岩的细观力学参数，无法通过试验直接测得，但其近似值可由试验结果计算得到，并在试算基础上调整优化.Weibull分布是描述材料缺陷和强度统计分布的经典模型，其中，均质度系数（m）是反映煤岩细观单元力学参数非均质性大小的参数，m越大，材料越均质；m越小，材料非均质性越强.因此，对m赋值并用式（25）和式（26）试算平均细观弹性模量（E_m）和平均细观单轴抗压强度（σ_m）^［39］：

E_{\mathrm{m}}=\frac{E}{\mathrm{0.1412}\mathrm{l}\mathrm{n}m+\mathrm{0.6476}}\mathrm{ }(\mathrm{1}.\mathrm{2}\le m\le \mathrm{10}),（25）

{\sigma }_{\mathrm{m}}=\frac{{\sigma }_{\mathrm{c}}}{\mathrm{0.2602}\mathrm{l}\mathrm{n}m+\mathrm{0.0233}}\mathrm{ }(\mathrm{1}.\mathrm{2}\le m\le \mathrm{50}).（26）

为了验证m赋值的科学性及合理性，获取E_m及σ_m后，利用RFPA软件构建单轴压缩数值模型，钻取2个尺寸⌀50 mm×100 mm的标准柱塞煤样并开展平行的单轴压缩试验，用2个煤样的单轴压缩试验结果检验、标定数值模拟结果，调整参数m，直至数值模拟试样的峰值强度、应力⁃应变曲线与试验结果高度契合，得到可以反映煤岩宏观力学性质的m、E_m和σ_m（图15）.借助RFPA软件构建与水力压裂物理模拟煤样尺寸相同的数值模型，对数值模型的参数进行赋值（表2），开展真三轴水力压裂数值模拟.

2.3 压裂数值模拟结果

在水力压裂物理模拟基础上，利用数值模拟方法探究排量和液量对裂隙起裂、扩展的影响.数值模拟排量与物理模拟一一对应，分别是3，6，9，12，15 mL/min.数值模拟过程中，每计算步代表5 s，为了确保不同模拟样品注入压裂液量相同，5组模拟分别计算150，75，50，38，30步，使压裂液量相同，探究排量对压裂效果的影响.

当Q=3 mL/min时，沿多个方向形成复杂裂缝网络，从x轴、y轴中间切片可以看出，裂隙带集中于试样中下部，从z轴中间切片可以看出试样中部裂隙最发育；该试样形成了扩展范围较有限、缝宽较小的复杂裂隙带，绝大多数部位缝宽小于1 mm（图16a）.当Q=6 mL/min时，裂隙带的规模与3 mL/min试样相当，但形状存在差异.该试样裂隙带形状呈以x轴中间切面为“赤道板”的纺锤体形.x轴中间切片的面积大于y轴和z轴中间切片，且x轴中间切片附近裂隙的缝宽更大，表明主裂隙垂直于σ_h.与3 mL/min压裂条件相比，6 mL/min压裂形成主裂隙的缝宽更大（图16b）.当Q=9 mL/min时，裂隙带的规模与3和6 mL/min试样相当.从x轴、y轴、z轴中间切面可以看出，不同方向裂隙带的发育程度相当，且沿3个方向均形成了缝宽较大的主裂隙，表明9 mL/min压裂试样形成的裂缝网络具有更强的导流能力，压裂效果更好（图16c）.当Q=12 mL/min时，未形成复杂裂缝网络，而是沿近似垂直于σ_h方向形成了缝宽更大的简单裂隙.从x轴中间切片可以看出，裂隙自试样中心沿平行于σ_V和σ_H方向扩展；自中心向外，缝宽逐渐减小（图16d）.当Q=15 mL/min时，同样形成了近似垂直于σ_h的简单裂隙（图16e）.与12 mL/min压裂裂隙带相比，15 mL/min压裂形成裂隙带沿σ_V和σ_H方向的长度略小于前者，局部区域的裂隙缝宽略大于前者，表明同样是形成简单裂隙，在液量相同情况下，大排量形成裂隙的缝宽更大.

总结水力压裂数值模拟裂隙带分布特征，发现在压裂液量固定条件下，随排量增大，压裂裂隙的复杂程度先增大后降低.具体表现在，低排量压裂形成缝宽较小的复杂裂缝网络；随排量增大，形成缝宽较大的复杂裂缝网络；当排量进一步增大，形成缝宽显著增大的简单裂隙.在本文试验条件下，水力压裂物理模拟和数值模拟得到了相似的规律，即均在排量9 mL/min时达到最佳改造效果，当排量低于9 mL/min时形成复杂裂隙带，当排量大于9 mL/min时形成简单裂隙.

选取压裂效果最好的排量9 mL/min的数值模拟试样，分析压裂液量对裂隙带扩展的影响，结果如图17所示.图中，蓝色线条为水压⁃时间步曲线，红色团簇分别为第1（0.75 mL）、第5（3.75 mL）、第10（7.50 mL）、第15（11.25 mL）、第20（15.00 mL）、第25（18.75 mL）、第30（22.50 mL）、第35（26.25 mL）、第40（30.00 mL）、第45（33.75 mL）、第50（37.50 mL）步的裂隙带形貌.可以看出，煤样在第30步达到破裂压力，30步以前，裂隙集中于井筒附近，规模非常有限，且扩展缓慢，表明煤样在达到破裂压力前仅在井筒周边形成微裂隙.当水压达到破裂压力后，裂隙沿多个方向迅速向外扩展，裂隙带范围扩大.对比不同时间步裂隙带的形貌特征，发现在达到破裂压力前，裂隙带的规模随液量增大不发生显著变化，达到破裂压力后，裂隙带的规模随液量增大迅速增大.水压增大过程和裂隙带的扩展特征表明，破裂前的时间步占总时间步的2/3，为破裂后裂隙带迅速扩展积累了大量压强能，是压裂裂隙起裂扩展最主要的动力和能量来源.

3 深部煤储层水力压裂的能量控因

3.1 能量计算方法

以压裂液为载体的外部能量的持续输入是煤在水力压裂过程中损伤并破坏的本质，人工裂缝的形成是煤岩压裂过程中能量不可逆耗散的外在表现^［40］.水力压裂物理模拟试样内部不含气体，储层压力可视为0 MPa，因此计算总能量时可忽略气体弹性能和水体弹性能的影响，仅考虑煤基质在地应力影响下储存的煤基块弹性能（U_e）.

由于煤储层在原位条件下处于三向应力状态，结合广义胡克定律，可得3个主应力方向应变分量的表达式：

\left\{\begin{array}{l}{\epsilon }_x=\frac{\mathrm{1}}{E}[{\sigma }_x-\mathrm{\nu }({\sigma }_y+{\sigma }_z)],\\ {\epsilon }_y=\frac{\mathrm{1}}{E}[{\sigma }_y-\mathrm{\nu }({\sigma }_x+{\sigma }_z)],\\ {\epsilon }_z=\frac{\mathrm{1}}{E}[{\sigma }_z-\mathrm{\nu }({\sigma }_x+{\sigma }_y)],\end{array}\right.（27）

式中：σ_x，σ_y，σ_z 分别为x，y，z方向的主应力，MPa；ε_x，ε_y，ε_z 分别为x，y，z方向的应变，%.

实际地应力条件下，分别用σ_H，σ_h，σ_V替代σ_x，σ_y，σ_z，用ε_H，ε_h，ε_V代替ε_x，ε_y，ε_z .U_e可由煤样的应力及应变共同表征：

U_e=\frac{\mathrm{1}}{\mathrm{2}}σ_H{\epsilon }_{\mathrm{H}}^{\mathrm{e}}+\frac{\mathrm{1}}{\mathrm{2}}σ_h{\epsilon }_{\mathrm{h}}^{\mathrm{e}}+\frac{\mathrm{1}}{\mathrm{2}}σ_V{\epsilon }_{\mathrm{V}}^{\mathrm{e}}，（28）

式中：{\epsilon }_{\mathrm{H}}^{\mathrm{e}}，{\epsilon }_{\mathrm{h}}^{\mathrm{e}}，{\epsilon }_{\mathrm{V}}^{\mathrm{e}}分别为σ_H，σ_h，σ_V方向的弹性应变，%.

将式（27）与式（28）相结合，得到U_e的表达式：

\begin{array}{l}{U}_{\mathrm{e}}\approx \frac{\mathrm{1}}{\mathrm{2}E}[{\sigma }_{\mathrm{V}}^{\mathrm{2}}+{\sigma }_{\mathrm{H}}^{\mathrm{2}}+{\sigma }_{\mathrm{h}}^{\mathrm{2}}-\\ \mathrm{2}\nu ({\sigma }_{\mathrm{V}}{\sigma }_{\mathrm{H}}+{\sigma }_{\mathrm{H}}{\sigma }_{\mathrm{h}}+{\sigma }_{\mathrm{V}}{\sigma }_{\mathrm{h}})]\end{array} . （29）

对于水力压裂物理模拟，受试样尺寸影响，裂缝形成时间短暂，研究压裂起裂特征更有意义，因此本文着重分析水力压裂阶段Ⅱ的能量演化特征及其对起裂破坏过程的控制作用.本研究中，假设压裂液为不可压缩流体，具有恒定的黏度和密度.对于物理模拟试样而言，从压裂开始至结束，压裂液携带的能量是唯一的能量来源.由于水力压裂物理模拟管线较短，活性水压裂液的黏度低，且压裂管线表面光滑，摩阻小，可忽略管线中的能量损失.因此，压裂输入煤样的能量由两部分组成，分别是压强能（U_ac）和动能（U_v）.U_ac、U_v和输入能（U）可分别由式（30）、式（31）和式（32）计算得到：

U_{\mathrm{a}\mathrm{c}}={\int }_{t_{\mathrm{1}}}^{t_{\mathrm{2}}}\frac{{\sigma }_{\mathrm{w}}}{{\rho }_{\mathrm{a}\mathrm{c}}}{m}_{\mathrm{a}\mathrm{c}}\mathrm{d}t,（30）

U_{\mathrm{v}}={\int }_{t_{\mathrm{1}}}^{t_{\mathrm{2}}}\frac{m_{\mathrm{a}\mathrm{c}}{v}^{\mathrm{2}}}{\mathrm{2}}\mathrm{d}t,（31）

U=U_{\mathrm{a}\mathrm{c}}+U_{\mathrm{v}}.（32）

式中：m_ac为单位时间注入压裂液质量，kg；v为液体运移的速度，m/s.

输入能中，形成人工裂隙耗散的表面能对压裂效果影响显著，可由式（33）计算：

U_{\mathrm{f}\mathrm{r}\mathrm{a}\mathrm{c}}=G_{\mathrm{I}\mathrm{C}}{S}_{\mathrm{f}},（33）

式中：U_frac为表面能，J；G_IC为基于能量释放率的断裂韧性，J/m²；S_f为裂隙总表面积，通过Avizo提取水力压裂物理模拟煤样的三维裂隙得到，m².

其中，G_IC可由下式计算：

G_{\mathrm{I}\mathrm{C}}=\frac{{K^{\mathrm{2}}}_{\mathrm{I}\mathrm{C}}}{E\mathrm{\text{'}}},（34）

式中：E'为有效弹性模量，MPa；K_IC为Ⅰ型应力强度因子，MPa·m^1/2，可由下式求得：

K_{\mathrm{I}\mathrm{C}}=P_{\mathrm{t}}\sqrt[]{\mathrm{\pi }a},（35）

式中：P_t为抗拉强度，通过巴西劈裂试验测得，MPa；a为裂隙长度，等于圆盘煤样半径，m.

水力压裂过程一般被视为平面应变问题，则E'由下式得到：

E\text{'}=\frac{E}{\mathrm{1}-{\nu }^{\mathrm{2}}}.（36）

3.2 能量演化特征及其对压裂效果的控制作用

将水力压裂物理模拟阶段Ⅱ的试验结果代入式（30）和式（31），发现不同排量压裂煤样的U_v介于1.16×10^-7~1.95×10^-6 J，且两者线性正相关（图18a）；U_ac介于287.84~562.82 J，随排量增大，U_ac呈先增大后减小趋势（图18b）.对比U_v与U_ac发现，前者比后者小8~9个数量级，即在水力压裂物理模拟过程中，绝大部分U由U_ac提供，U_v的贡献很小，可忽略不计，因此可将U_ac近似视为U.

能量输入煤样后，部分U在人工裂隙形成过程中耗散，属于水力压裂的“有效能量”.将煤样的P_t、E、ν等力学参数与水力压裂物理模拟试样的裂隙表面积提取结果代入式（33）~（36），计算发现，裂隙扩展耗散的U_frac介于85.85~368.96 J，随排量增大，U_frac呈先增大后减小趋势（图18c）.U_ac、U_frac与排量之间均呈二阶多项式关系，且变化趋势相似，无论是从输入总能量角度，还是人工裂隙扩展耗散的表面能角度，均表明存在较优排量区间，将压裂排量设定在该区间范围内有助于形成大规模体积缝网.在此基础上，结合裂隙率与排量的关系，可知对于本研究模拟条件，排量6~12 mL/min属于较优排量区间，9 mL/min时压裂改造效果最理想（图13）.

对比不同排量压裂时U与改造效果之间的差异，发现在压裂阶段Ⅱ中，排量为3 和15 mL/min煤样的U分别为290.50和287.84 J，两者排量差异巨大，但U的差异很小（图18b）.结合不同排量压裂物理模拟得到的三维裂隙重构结果（图12）、裂隙分布和裂隙率变化特征（图13）及数值模拟得到的裂隙带分布特征（图16），可知同样的U可能取得不同的压裂改造效果，即能量随压裂液输入煤样的过程和输入能本身同样重要.此外，将三轴压缩试验得到相应地应力状态下煤样的弹性模量和泊松比代入式（29），并按照样品尺寸计算得到压裂物理模拟试样所储存的煤基块弹性能为41.32 J.将U与U_e对比，发现前者是后者的6.97~13.62倍，这表明，水力压裂起裂前输入煤层的能量远大于煤基质在原位条件下储存的U_e.

为了对比不同排量压裂能量输入过程的差异，定义“能量输入速率”为压裂过程中目标时段内单位时间的能量输入量，在本文中指阶段Ⅱ输入煤样的总能量与该阶段总时间的比值，其表达式为

u=\frac{U_{\mathrm{Ⅱ}}}{t_{\mathrm{Ⅱ}}}，（37）

式中：u为能量输入速率，J/s；U_Ⅱ 为阶段Ⅱ输入能，J；t_Ⅱ 为阶段Ⅱ总时间，s.

研究发现，阶段Ⅱ排量越大，能量输入速率越快，两者呈对数函数关系，表明随排量增大，能量输入速率的增速降低，换言之，增大排量能够提高能量输入速率，但能量输入的效率在降低（图19a）.进一步对比不同排量压裂阶段Ⅱ能量输入速率增量的变化，发现中低排量下单位排量增量（3 mL/min）对应能量输入速率的增量整体大于中高排量，这也是对大排量压裂能量输入效率降低的佐证（图19b）.上述分析表明，增大压裂排量可以提高阶段Ⅱ能量输入速率，有助于煤样内部能量在短时间内迅速积累和耗散，实际压裂过程中可以促使压裂裂缝突破微小的天然弱面，更容易在地应力状态控制下形成大规模裂隙.当压裂排量进一步增大，能量输入速率也随之增大，在其超过某一阈值后，煤样迅速破坏并形成延伸较远的简单裂隙，难以形成复杂的高效导流通道.

进一步分析裂隙率与阶段Ⅱ能量输入速率的关系，发现随能量输入速率增大，裂隙率呈先增大后减小趋势，两者呈高斯分布（图19c）.结合不同排量水力压裂煤样CT扫描的裂隙分布特征，煤样裂隙率在压裂排量为9 mL/min时达到最大值.将这一排量代入图19a中能量输入速率与排量的拟合关系式，得到使本文水力压裂物理模拟压裂效果由好转差的能量输入速率的阈值为1.01 J/s.因此，基于数据拟合结果，若以峰值裂隙率的90%作为评价压裂效果优劣的依据，则对于本研究边长100 mm的立方体煤样，当阶段Ⅱ能量输入速率介于0.70~1.15 J/s时，压裂改造效果最好.

4 结　论

1）钻孔内部填充NaCl粉末、抗拉拔植筋胶加固、试样外侧橡胶垫片缓冲的真三轴水力压裂物理模拟优化方法具有制样过程简单、钢管固定可靠、井筒密封性良好、缓冲层有效等优点，解决了深部三向应力施加过程中煤样易破碎的问题，可以取代在煤样外表面浇筑水泥以保护煤样的常规方法.

2）构建了深部煤储层水力压裂应力分布模型，随煤基质所受水压σ_w增大，井筒附近环向应力σ_φ 由压应力转为拉应力，是裂隙起裂的主因.径向应力σ_ρ 随σ_w增大而增大，水平方向小于垂直方向；σ_φ 在井筒附近出现应力集中，并随σ_w增大而减小，且水平方向大于垂直方向.深部正断型应力状态下，随σ_w增大，水力压裂裂隙不仅优先形成于最大水平主应力σ_H和垂向应力σ_V所在的平面，且首先沿平行于σ_V方向起裂.

3）对于深部正断型应力状态，排量3~15 mL/min范围内，随排量增大，依次形成小规模复杂缝网⁃大规模复杂缝网⁃超大规模简单裂缝.自井筒向外，裂隙带分为微裂隙密集区和主裂隙扩展区，前者煤体从钻孔内壁剥离与煤岩弹塑性力学性质、非均质性、基质单元差异应变等相关；后者宏观裂隙自微裂隙密集区向外扩展，以拉张裂隙为主、剪切裂隙为辅.

4）水力压裂物理模拟裂隙起裂的动力来自于阶段Ⅱ积累的输入能U，主要由压强能U_ac贡献.裂隙率、U_ac、表面能U_frac与排量之间均呈二阶多项式关系，排量6~12 mL/min属于较优压裂区间，9 mL/min时最佳.能量输入速率与排量呈对数关系，与裂隙率呈高斯分布，阶段Ⅱ能量输入速率介于0.70~1.15 J/s时，压裂改造效果最好.

参考文献

[1] 秦　勇. 中国深部煤层气地质研究进展［J］. 石油学报， 2023，44（11）： 1791⁃1811.

[2] QIN Yong. Progress on geological research of deep coalbed methane in China［J］. Acta Petrolei Sinica， 2023，44（11）： 1791⁃1811.

[3] 桑树勋， 郑司建， 刘世奇， 等. 煤系气及深部煤层气高效勘探开发若干研究进展［J］. 中国矿业大学学报， 2025，54（1）： 1⁃25.

[4] SANG Shuxun， ZHENG Sijian， LIU Shiqi， et al. Research advances in efficient exploration and development of coal measure gases and deep coallbed methane［J］. Journal of China University of Mining & Technology， 2025，54（1）： 1⁃25.

[5] 陈贞龙， 王运海， 刘　晓， 等. 延川南深部煤层气开发关键技术与地质工程一体化实践［J］. 煤田地质与勘探， 2025，53（1）： 142⁃151.

[6] CHEN Zhenlong， WANG Yunhai， LIU Xiao， et al. Critical technologies and geology⁃engineering integration practices for deep CBM production in the Yanchuannan CBM field［J］. Coal Geology & Exploration， 2025，53（1）： 142⁃151.

[7] 杨兆彪， 李存磊， 郭巧珍， 等. 新疆准噶尔盆地白家海凸起深部煤层气不同赋存态分配规律［J］. 中国矿业大学学报， 2024，54（1）： 127⁃137.

[8] YANG Zhaobiao， LI Cunlei， GUO Qiaozhen， et al. Distribution patterns of various occurrence states of deep coalbed methane in the Baijiahai Uplift， Junggar Basin， Xinjiang［J］. Journal of China University of Mining & Technology， 2024，54（1）： 127⁃137.

[9] 房大志， 程泽虎， 李佳欣. 渝东南地区超深层煤层气高效压裂技术及精细排采制度研究与实践：以NY1井为例［J］. 煤田地质与勘探， 2022，50（5）： 50⁃56.

[10] FANG Dazhi， CHENG Zehu， LI Jiaxin. Eefficient fracturing technology and fine drainage system of ultra⁃deep coalbed methane in southeast Chongqing： A case study of NY1 well［J］. Coal Geology & Exploration， 2022，50（5）： 50⁃56.

[11] 徐凤银， 熊先钺， 侯　伟， 等. 深部煤层气产业升级与“八个一体化”体系的建立［J］. 石油学报， 2025，46（2）： 289⁃305.

[12] XU Fengyin， XIONG Xianyue， HOU Wei， et al. Upgrading of deep coalbed methane industry and establishment of the “Eight⁃in⁃One” system［J］. Acta Petrolei Sinica， 2025，46（2）： 289⁃305.

[13] 方志明， 王润东， 杨晨龙. 深部煤层气开发的机遇与挑战［J］. 中国矿业大学学报， 2025，54（1）： 34⁃51.

[14] FANG Zhiming， WANG Rundong， YANG Chenlong. Opportunities and challenges in deep coalbed methane development［J］. Journal of China University of Mining & Technology， 2025，54（1）： 34⁃51.

[15] 许　浩， 汤达祯， 陶　树， 等. 深、浅部煤层气地质条件差异性及其形成机制［J］. 煤田地质与勘探， 2024，52（2）： 1⁃7.

[16] XU Hao， TANG Dazhen， TAO Shu， et al. Differences in geological conditions of deep and shallow coalbed methane and their formation mechanisms［J］. Coal Geology & Exploration， 2024，52（2）： 1⁃7.

[17] WANG X L， PAN J N， WANG H C， et al. Transformation effect of multi⁃scale pores induced by true triaxial supercritical carbon dioxide fracturing in high⁃rank coal［J］. Gas Science and Engineering， 2025，139： 205644.

[18] LI S， QIN Y， TANG D Z， et al. A comprehensive review of deep coalbed methane and recent developments in China［J］. International Journal of Coal Geology， 2023，279： 104369.

[19] 邱宇超， 梁卫国， 李　静， 等. 非均质弹塑性煤体水压致裂裂纹形态研究［J］. 煤炭学报， 2022，47（10）： 3668⁃3679.

[20] QIU Yuchao， LIANG Weiguo， LI Jing， et al. Study on fracture morphology of hydraulic fracturing in heterogeneous elastoplastic coal［J］. Journal of China Coal Society， 2022，47（10）： 3668⁃3679.

[21] 刘长华， 陈国青， 高　宇， 等. 沁水盆地深部软煤煤层气开发实践［J］. 矿业安全与环保， 2024，51（2）： 141⁃145.

[22] LIU Changhua， CHEN Guoqing， GAO Yu， et al. Development practice of deep soft coal CBM in Qinshui Basin［J］. Mining Safety & Environmental Protection， 2024，51（2）： 141⁃145.

[23] ZHOU D， LYU Z Q， CAO Y X， et al. Fractal evolution characteristics of pore structure in coal⁃acidified stimulation［J］. Fractal and Fractional， 2025，9（2）： 62.

[24] 徐凤银， 聂志宏， 孙　伟， 等. 鄂尔多斯盆地东缘深部煤层气高效开发理论技术体系［J］. 煤炭学报， 2024，49（1）： 528⁃544.

[25] XU Fengyin， NIE Zhihong， SUN Wei， et al. Theoretical and technological system for Highly efficient development of deep coalbed methane in the eastern edge of Erdos Basin［J］. Journal of China Coal Society， 2024，49（1）： 528⁃544.

[26] YAN X， YU H T. Numerical simulation of hydraulic fracturing with consideration of the pore pressure distribution based on the unified pipe⁃interface element model［J］. Engineering Fracture Mechanics， 2022，275： 108836.

[27] 房孝杰， 吴财芳， 申　建， 等. 深部煤储层温⁃压耦合力学响应机制及裂隙扩展规律：以延川南区块为例［J/OL］. 煤炭学报， 1⁃17 ［2025⁃10⁃16］.

[28] FANG Xiaojie， WU Caifang， SHEN Jian， et al. Mechanical response mechanism and fracture propagation law of deep coal reservoirs under temperature⁃pressure coupling conditions： In a case of the southern Yanchuan block［J/OL］. Journal of China Coal Society， 1⁃17 ［2025⁃10⁃16］.

[29] ZOU Y S， ZHAO Z W， ZHANG S C， et al. Adjustment of fracture network geometry during temporary plugging and diverting fracturing in deep coalbed methane reservoirs： An experimental study［J］. Physics of Fluids， 2024，36（12）： 126628.

[30] WANG X L， PAN J N， WANG K， et al. Fracture variation in high⁃rank coal induced by hydraulic fracturing using X⁃ray computer tomography and digital volume correlation［J］. International Journal of Coal Geology， 2022，252： 103942.

[31] QIN J X， LIU W， WEI Y L， et al. Characteristics of hydraulic crack propagation in adjacent coal seams for coalbed methane production［J］. Energy & Fuels， 2025，39（6）： 3069⁃3086.

[32] HAN L， MISCH D， SHEN J， et al. Influence of viscous and capillary forces on residual water in anthracite fracture networks［J］. Natural Resources Research， 2023，32（2）： 603⁃617.

[33] CHEN B， LI S， TANG D Z， et al. Numerical simulation study on hydraulic fracture propagation of multi⁃cluster fracturing of horizontal well in deep fractured coal seams［J］. Engineering Fracture Mechanics， 2025，318： 110983.

[34] 徐凤银， 甄怀宾， 李曙光， 等. 深部煤层气储层改造技术迭代升级历史与发展方向：以鄂尔多斯盆地东缘大吉区块为例［J］. 煤炭科学技术， 2025，53（3）： 1⁃18.

[35] XU Fengyin， ZHEN Huaibin， LI Shuguang， et al. History and development direction of iterative upgrading of deep coalbed methane reservoir reconstruction technology： Taking the Daji Block in the eastern margin of the Ordos Basin as an example［J］. Coal Science & Technology， 2025，53（3）： 1⁃18.

[36] FANG X J， WU C F， ZHANG H W， et al. Stress distribution properties and deformation⁃fracture mechanisms in hydraulic fracturing of coal［J］. Fuel， 2023，351： 129049.

[37] 陈世达， 汤达祯， 侯　伟， 等. 深部煤层气地质条件特殊性与储层工程响应［J］. 石油学报， 2023，44（11）： 1993⁃2006.

[38] CHEN Shida， TANG Dazhen， HOU Wei， et al. Geological particularity and reservoir engineering response of deep coalbed methane［J］. Acta Petrolei Sinica， 2023，44（11）： 1993⁃2006.

[39] 吴财芳， 秦　勇， 傅雪海， 等. 煤基块弹性能及其与地质控制因素之间的关系［J］. 中国矿业大学学报， 2005，34（5）： 636⁃639.

[40] WU Caifang， QIN Yong， FU Xuehai， et al. Coal matrix flexibility energy and the relation with geological controlling factors［J］. Journal of China University of Mining & Technology， 2005，34（5）： 636⁃639.

[41] 姜　玮， 吴财芳. 织纳煤田煤储层弹性能及其对有利选区的控制作用［J］. 煤炭学报， 2011，36（10）： 1674⁃1678.

[42] JIANG Wei， WU Caifang. Flexibility energy of coal⁃bed gas reservoir in Zhina coalfield and their control function on the favorable areas［J］. Journal of China Coal Society， 2011，36（10）： 1674⁃1678.

[43] 吴财芳， 秦　勇. 煤储层弹性能及其控藏效应：以沁水盆地为例［J］. 地学前缘， 2012，19（2）： 248⁃255.

[44] WU Caifang， QIN Yong. Flexibility energies of coal⁃bed gas reservoir and the controlling function on coal⁃bed gas reservoir formation： A case study from Qinshui Basin［J］. Earth Science Frontiers， 2012，19（2）： 248⁃255.

[45] LIU X L， JIA T R， WEI J P， et al. Modeling and evolution characteristics of coal reservoir energy during gas and water production［J］. Journal of Natural Gas Science and Engineering， 2021，96： 104329.

[46] 柳贡慧， 庞　飞， 陈治喜. 水力压裂模拟实验中的相似准则［J］. 石油大学学报（自然科学版）， 2000，24（5）： 45⁃48.

[47] LIU Gonghui， PANG Fei， CHEN Zhixi. Development of scaling laws for hydraulic fracture simulation tests［J］. Journal of the University of Petroleum， China， 2000，24（5）： 45⁃48.

[48] TIAN Z Z， XIONG Z， WEI Y H， et al. Numerical simulation of fracture propagation in deep coal seam reservoirs ［J］. Energy Science & Engineering， 2023，11（10）： 3559⁃3574.

[49] ZHAO Z H， WU T Y， GUO J C， et al. Study on the spatial and temporal evolution law of large⁃size fracture propagation in deep coalbed based on acoustic emission technology［J］. Rock Mechanics and Rock Engineering， 2025，58（3）： 3799⁃3814.

[50] 朱　葛. 变流量水力压裂井底动压及储层动应力与裂缝扩展规律的研究［D］. 秦皇岛： 燕山大学， 2022.

[51] ZHU Ge. Study on bottom hole dynamic pressure， reservoir dynamic stress and fracture propagation law of variable flow hydraulic fracturing［D］. Qinhuangdao： Yanshan University， 2022.

[52] 蔡美峰， 何满潮， 刘东燕. 岩石力学与工程（第二版）［M］.北京： 科学出版社， 2013.

[53] CAI Meifeng， HE Manchao， LIU Dongyan. Rock mechanics and engineering （second edition） ［M］. Beijing： Science Press， 2013.

[54] 屈　晶， 申　建， 韩　磊， 等. 基于CT图像的高阶煤不同宏观煤岩组分裂隙差异发育规律［J］. 天然气工业， 2022，42（6）： 76⁃86.

[55] QU Jing， SHEN Jian， HAN Lei， et al. Characteristics of fractures in different macro⁃coal components in high⁃rank coal based on CT images［J］. Natural Gas Industry， 2022，42（6）： 76⁃86.

[56] 朱万成， 唐春安， 杨天鸿， 等. 岩石破裂过程分析（RFPA^2D）系统的细观单元本构关系及验证［J］. 岩石力学与工程学报， 2003，23（1）： 24⁃29.

[57] ZHU Wancheng， TANG Chun̓an， YANG Tianhong， et al. Constitutive relationship of mesoscopic elements used in RFPA^2D and its validations［J］. Chinese Journal of Rock Mechanics and Engineering， 2003，23（1）： 24⁃29.

[58] 夏英杰， 李连崇， 唐春安， 等. 储层砂岩破坏特征与脆性指数相关性影响的试验及数值研究［J］. 岩石力学与工程学报， 2016，36（1）： 10⁃28.

[59] XIA Yingjie， LI Lianchong， TANG Chun̓an， et al. Experiment and numerical research on failure characteristic and brittleness index for reservoir sandstone［J］. Chinese Journal of Rock Mechanics and Engineering， 2016，36（1）： 10⁃28.

[60] 李志超. 低渗储层水力裂缝扩展特征的数值模拟研究［D］. 大连： 大连理工大学， 2019.

[61] LI Zhichao. Numerical research on the propagation behavior of hydraulic fractures in the low⁃permeability reservoirs［D］. Dalian： Dalian University of Technology， 2019.

[62] 朱　颖. 油页岩层理对其力学特性及裂缝起裂与扩展的影响研究［D］. 长春： 吉林大学， 2022.

[63] ZHU Ying. Study on the influence of oil shale bedding planes on its mechanical properties， fracture initiation and extension［D］. Changchun： Jilin University， 2022.

[64] 李连崇， 唐春安， 梁正召， 等. 煤层底板陷落柱活化突水过程的数值模拟［J］. 采矿与安全工程学报， 2009，26（2）： 158⁃162.

[65] LI Lianchong， TANG Chun̓an， LIANG Zhengzhao， et al. Numerical simulation on water inrush process due to activation of collapse columns in coal seam floor［J］. Journal of Mining & Safety Engineering， 2009，26（2）： 158⁃162.

[66] 谢和平， 鞠　杨， 黎立云， 等. 岩体变形破坏过程的能量机制［J］. 岩石力学与工程学报， 2008，27（9）： 1729⁃1740.

[67] XIE Heping， JU Yang， LI Liyun， et al. Energy mechanism of deformation and failure of rock masses［J］. Chinese Journal of Rock Mechanics and Engineering， 2008，27（9）： 1729⁃1740.

基金资助

中央高校基本科研业务费专项资金项目(2025QN1072)