四维动态地应力建模方法及其在寿阳区块煤层气开发中的应用*

张滨海 陈峥嵘 艾传志 李莹莹 彭成勇 朱海燕 唐煊赫

（1.中海油研究总院有限责任公司 北京 100028； 2.西南石油大学 油气藏地质及开发工程国家重点实验室 四川成都 610500）

我国煤层气地质储量丰富，以沁水盆地为煤层气主要产区，但该区域煤层气开发单井产量长期处于较低水平。沁水盆地北部的寿阳区块主要产气层分别为3＃、9＃和15＃煤层，均采用水力压裂开发，2014年开发至今大部分开发井单井产气量不足200 m3／d，而产水量多高于10 m3／d。为了评估前期开发过程中压裂施工的合理性，同时为调整开发方案与二次压裂设计提供依据，需要对生产过程中地应力变化情况及当前地应力状态进行分析。

一般情况下，三维静态模型能够较为准确地反映当前地应力的状态，因此，目前常规油气藏开发过程中的储层应力状态分析多使用三维静态地应力模型。然而，一旦地层呈现强非均质性、储层岩石呈现各向异性、生产或注入过程中地应力场发生应力偏转甚至产生大量断层，若仍使用静态地应力分析方法，将必然出现分析误差。不仅如此，由于应力状态变化易出现非线性和不规律等特征，静态地应力分析方法将无法预测强非均质性储层在生产过程中的地应力变化情况，无法准确获取当前地应力状态。因此，针对上述情况，需要利用油藏模拟器和地质力学模拟器进行耦合，建立考虑时间效应的四维地应力分析预测模型。

Biot［1］最早在Terzaghi的一维流动-应力耦合理论中率先提出了三维地应力模型，而后Geertsma［2］提出了孔隙和岩石体积变化理论，并讨论了地应力变化对岩石弹性和孔隙体积的影响。从20世纪90年代开始，由于水力压裂的推广应用、油气藏长效开采过程中需要考虑压实沉降等原因，在油藏模拟过程中考虑地应力变化得出了大量的研究成果，先后出现了大量渗流-应力耦合模型［3-9］。其中，有限元方法较其他方法如有限差分、离散元等具有更好的适应性，先后发展出了Teufel等［10-11］和Chen等［13-14］提出的单孔隙和双孔隙渗流-应力耦合模型，Cuisiat等［15］、Gutierrex等［16-17］提出的全耦合有限元模型，Settari等［18］、Tran等［19］为解决全耦合收敛性差而提出的部分耦合有限元模型，Onaisi等［20］、Samier等［21-22］提出的显示积分耦合模型等。在此基础上，Hatchell等［23］、Herwanger等［24］、Onaisi等［25］发展出了考虑生产或注入过程时间效应的四维动态地应力模型，并且不断发展出各向异性地层分析和渗透率应力敏感性分析等应用。近年来，商业软件逐渐在渗流-应力耦合中发挥出强大的优势，其中应用较广、认可度较高的耦合系统主要有以CodeBright［26］、COMSOL［27］等为代表的有限差分法全耦合，以TOUGH和FLAC3D耦合［28］、ECLIPSE与ABAQUS耦合［29］等为代表的有限元法部分耦合，以ECLIPSE与VISAGE耦合［30］、ATHOS与ABAQUS耦合［31］为代表的有限元法单向耦合。然而，上述基于商业软件的多场耦合方法存在2个方面不足：一是耦合计算过程中多假设理想产量，并未基于煤层气的真实排采情况，无法反映实际生产过程中煤层气藏地应力的变化；二是未考虑储层各向异性和非均质性影响。

因此，本文基于有限差分算法油藏数值模拟器CMG和非线性有限元模拟器ABAQUS，提出了渗流-应力耦合四维动态地应力建模分析方法，并结合沁水盆地寿阳煤层气区块地质储层及岩石力学特征，以寿阳区块15＃煤层为研究对象，建立寿阳南燕竹A2区块地应力模型，分析该区块15＃煤层排水采气过程中地应力的变化情况，以期为煤层气排采制度调整、加密井井壁稳定性分析、重复压裂优化设计等提供基础数据。

1 四维地应力耦合建模方法

本文提出的四维地应力耦合建模方法，其流程是：首先在地质建模软件中进行地质建模及相关处理；然后利用CMG作为油藏模拟器，选择ABAQUS作为考虑时间效应的四维动态地应力有限元模拟器，利用自编程接口程序使得CMG和ABAQUS能够进行耦合计算，从而实现以CMG＋ABAQUS为平台进行四维动态地应力联合建模。

1.1 地质及油藏模型准备

在进行油藏模型和地应力模型耦合之前，首先需要地质建模，即建立含有地层孔渗参数（孔隙度、渗透率、孔隙压力、沉积相、岩性等）和岩石力学参数（弹性模量、泊松比等）的地质模型；然后将地质模型导入到油藏模拟器CMG中建立油藏模型，并根据目标地层的生产或注入数据及油藏渗流特性进行油藏动态模拟；最后分别将三维静态油藏模型和三维动态孔隙压力场以通用文件格式输出，即完成地质及油藏模型准备工作。

1.2 几何建模与网格处理

CMG油藏模型通常使用六面体（笛卡尔）网格，因此CMG中油藏模型的几何节点和网格可直接转换为有限元网格。利用油藏网格直接建立有限元网格的方法已经得到了应用，但该方法存在致命的缺陷，即当模型中存在断层或岩性尖灭时，必须通过降低网格尺寸来合理描述断层或岩性尖灭等地质特征，这势必使得网格密度迅速增加或出现畸形网格，最终导致计算难以收敛。此外，为了精确分析井筒附近的应力变化情况，必须要增大网格密度，而油藏模拟器无法提供井筒附近网格加密功能，因此就必须增大整体网格的密度，从而增大了计算收敛的难度。图1a为油藏模型网格及直接导入ABAQUS建立有限元模型中断层处理方法，由于断层处的不连续，将会导致ABAQUS计算难以收敛；图2a为岩性尖灭处的处理方法，油藏模型中的通常做法为不断加密网格，直至无限逼近尖灭，但这一处理方法也会造成计算量急剧增大或难以收敛。

图1 断层处网格处理方法Fig.1 Grid treatment of faults

图2 岩性尖灭处网格处理方法Fig.2 Grid treatment of pinch-outs

为了解决上述收敛问题，本文采用四面体网格进行网格划分处理，如图1b、2b所示，即利用四面体网格的几何特性，使得在断层处各节点能够共享，在岩性尖灭处用数个网格就能逼近尖灭。由于油藏模型网格和有限元模型网格在几何形状上的差异，因此无法直接将油藏网格模型直接转换为有限元网格模型。针对这一情况，本文利用自编程接口程序实现2种几何模型的转换：首先利用接口程序识别CMG导出文件中的不同层位的多个几何模型，并且转换为Pro-E模型格式文件，然后导入到ABAQUS中建立地层整体几何模型，同时利用扫描方式在各个层面之间建立十节点的四面体网格模型，最后根据CMG中导入的井眼轨迹点，对井筒附近网格进行加密处理。需要特别注意的是，由于四面体网格本身在划分时比六面体网格更容易形成畸形网格，因此可能需要采取多次网格质量检查和重划分的方式完善模型网格。

1.3 地应力有限元模型属性赋值

由于油藏模型网格和有限元模型网格的差异，使得油藏属性难以从网格中心计算点直接赋值到有限元网格的积分点中。针对这一情况，本文在接口程序中开发了属性转换功能，实现了对有限元地应力模型进行属性赋值。该功能利用邻近点法搜索离油藏网格中心计算点距离最接近的地应力模型有限元网格积分点，然后将属性赋值到该积分点中。该方法的优点在于不需要进行额外的运算，仅通过不同点的距离对比即可快速赋值，大大节约了运算成本。

1.4 ABAQUS计算模型完善

完成几何模型建立及属性赋值后，为方便后期模型的调整与修正，本文采用ABAQUS Input文件快速建模，并严格按照ABAQUS Input文件格式进行处理。同时，为了提高建模及分析效率，模型几何参数（包括节点、单元、节点集和单元集等）均以独立文件的形式存储和调用。

1）材料修正与完善。对于CMG，储层的岩石力学特性参数往往只能描述为各向同性弹性材料，而对于大多数油气藏，其储层岩石多表现为一定程度或一定方向上的各向异性，因此在ABAQUS中需要通过修改合理的弹性材料模型进行描述，同时根据已有的材料来适当地描述塑性材料。

2）接触控制。对于各个层位之间的接触或是断层接触在应力状态改变而产生的滑移，必须通过Slip Tolerance控制在相对小并且合理的范围内。

3）初始条件及应力初始化。将CMG三维静态模型转化到ABAQUS中的各项孔渗参数、初始孔隙压力以及地应力预估梯度作为初始条件，在ABAQUS进行初始化地应力平衡，得到初始化地应力，然后根据单井地应力分析结果进行校正。

4）外边界条件。在施加模型外边界条件时，本文采用的是对模型底面和侧面进行位移和旋转约束。

5）内边界条件。在每一时间步中，以变化的孔隙压力作为内边界条件。

2 寿阳区块动态地应力分析

2.1 区域及地质概况

寿阳南燕竹区位于山西省中部沁水煤田北端，地处太行山脉西麓太原东山背斜之东南翼，地势西高东低，北高南低，海拔标高为980.00～1342.10 m；但相对高差不大，一般在100～200 m。本区构造简单，断层稀少，煤层气资源丰富，含气饱和度较好，具有广阔的商业开发前景。目前，需要针对寿阳南燕竹A2区块加密井、开发井二次压裂等进行开发方案设计。

南燕竹A2区块主产气层位为3＃、9＃和15＃煤层，其中3＃和9＃煤层厚度小，孔渗条件差，因此15＃煤层为重点开发层位，其储量丰度和含气面积均具有明显优势。15＃煤层位于太原组下部，其K2下灰岩层位顶板含水层充水，底板以泥岩、砂质泥岩为主，局部为细砂岩和炭质泥岩。南燕竹A2区块目前共钻井150口，其中压裂井127口，47口排井中产气井仅23口。

目前A2区块主要存在前期压裂完成投产后产气不足和大量产水的问题，主要表现为：①相似地质条件下，压裂排量越大，产水量越大；②压裂强度过大，易导致压裂裂缝过长，形成T字缝，造成煤层沟通含水层；③压裂段煤层薄，15＃煤层部分厚度不足2 m，气井大量产水。因此，需要对该区块在生产过程中的地应力变化情况进行研究，同时准确评估该区块当前地应力状态，为后期开发方案调整及重复压裂设计提供建议。

2.2 模型的建立及应力初始化验证

为了获得寿阳南燕竹A2区块当前地应力状态，并对地应力的变化趋势进行分析，利用本文提出的四维地应力建模方法，以Petrel中导出的地质模型为原始数据，在CMG中建立了A2区块三维油藏模型，并根据该区块47口在排井的生产动态数据完成了油藏数值分析，得到该区块不同时间节点的孔隙压力变化情况，同时利用接口程序耦合建立了A2区块有限元地应力模型。为了避免模型边缘不规则导致的收敛性问题，有限元模型包含了A2区块全部和部分A1区块。图3为寿阳南燕竹A2区块地应力分析有限元几何模型的网格划分示意图。

图3 寿阳南燕竹A2区块地应力分析有限元几何模型的网格划分示意图Fig.3 Grid division diagram of FE geomechanical model of NYZ-A2 block in Shouyang

利用接口程序对有限元模型进行属性赋值，并且结合多口井的地漏实验结果在ABAQUS地应力平衡步中进行地应力初始化反演，得到A2区块的三维初始地应力场。图4为A2区块15＃煤层初始最小水平主应力σh反演结果，可以看出初始最小水平主应力整体呈现由东北方向向西南方向逐渐增大的趋势，与该区块15＃煤层埋深的走势相似。

图4 寿阳南燕竹A2区块15＃煤层最小水平主应力初始化反演结果Fig.4 Minimum horizontal principal stress initialization of NYZ-A2 15＃coalbed in Shouyang

为了验证利用本文四维地应力建模分析方法反演得到的A2区块15＃煤层最小水平主应力σh、最大水平主应力σH结果的准确性，选取A2区块SYNY-101井的初始地应力在一定深度范围内的反演结果，并与利用现场测试（测井资料＋地破试验）得到的地应力纵向剖面进行对比，结果如图5所示。由图5可以看出，A2区块SYNY-101井井深450 m以上地层为走滑断层机制，450 m以下地层为正断层机制；地应力反演结果与现场测试结果整体吻合度较高，误差小于10%，说明该地应力建模方法及初始地应力反演结果合理、可靠。

图5 寿阳南燕竹A2区块SYNY-101井地应力现场测试分析结果与模拟结果对比Fig.5 Comparison between stress field test and simulation result of Well SYNY-101 in NYZ-A2 block in Shouyang

2.3 动态地应力计算结果分析

分别筛选A2区块地质条件较好的9口在排井进行井筒储层段附近地应力的变化情况分析，其中SYNY-125、SYNY-136、SYNY-150、SYNY-161、SYNY-173、SYNY-185、SYNY-187等7口井产气量相对较高、产水量较低，SYNY-135、SYNY-190等2口井无产气且产水量较高。

1）单井地应力变化过程对比分析。

首先选取两相邻开发井SYNY-136井和SYNY-135井，两口开发井储层埋深相近，地层孔渗条件相似，但排采情况完全不同。SYNY-136井排采情况良好，230 d内的产气量大于300 m3／d，产水量小于5 m3／d；而SYNY-135井排采状况较差，230 d内的产气量始终为0，且早已完成了排水降压阶段，而产水量大于10 m3／d。

图6、7分别为SYNY-136井井筒最大水平主应力、最小水平主应力随孔隙压力的变化情况，该井井筒附近储层沉降情况如图8所示。由图6、7可知，随着孔隙压力的降低，该井井筒水平主应力均上升，其中最大水平主应力变化量ΔσH＝0.197 6 MPa，变化率为1.34%；最小水平主应力变化量Δσh＝0.222 5 MPa，变化率2.04%。由图8可以看出，孔隙压力降低0.317 2 MPa，该井储层段沉降达到0.113 4 m。

图9、10分别为SYNY-135井井筒最大水平主应力、最小水平主应力随孔隙压力的变化情况，该井井筒附近储层沉降情况如图11所示。由图9、10可知，随着孔隙压力的降低，该井井筒水平主应力均上升，其中最大水平主应力变化量ΔσH＝0.162 9 MPa，变化率为1.07%；最小水平主应力变化量Δσh＝0.169 5 MPa，变化率1.46%。由图11可以看出，孔隙压力降低0.334 75 MPa，该井储层段沉降达到0.030 3 m。

对比SYNY-136井和SYNY-135井的地应力及沉降分析结果可知，对于南燕竹A2区块，两相邻井孔隙压力下降情况接近，说明该区块储层连通性较好；同时由于SYNY-136井处于良好的排采阶段，而SYNY-135井在230 d内始终未产气，因此SYNY-136井储层段地应力变化幅度要明显大于SYNY-135井，这一特征也同样反映在储层沉降幅度的差异上。

图6 寿阳南燕竹A2区块SYNY-136井最大水平主应力随孔隙压力变化Fig.6 The maximum horizontal principal stress change with pore pressure of Well SYNY-136 in NYZ-A2 block in Shouyang

图7 寿阳南燕竹A2区块SYNY-136井最小水平主应力随孔隙压力变化Fig.7 The minimum horizontal principal stress change with pore pressure of Well SYNY-136 in NYZ-A2 block in Shouyang

图8 寿阳南燕竹A2区块SYNY-136井储层沉降情况Fig.8 The subsidence of Well SYNY-136 in NYZ-A2 block in Shouyang

图9 寿阳南燕竹A2区块SYNY-135井最大水平主应力随孔隙压力变化Fig.9 The maximum horizontal principal stress change with pore pressure of Well SYNY-135 in NYZ-A2 block in Shouyang

图10 寿阳南燕竹A2区块SYNY-135井最小水平主应力随孔隙压力变化Fig.10 The minimum horizontal principal stress change with pore pressure of Well SYNY-135 in NYZ-A2 block in Shouyang

图11 寿阳南燕竹A2区块SYNY-135井储层沉降情况Fig.11 The subsidence of Well SYNY-135 in NYZ-A2 block in Shouyang

2）地应力变化强弱对比及重复压裂井优选。

将SYNY-125、SYNY-135、SYNY-136、SYNY-150、SYNY-161、SYNY-173、SYNY-185、SYNY-187和SYNY-190等9口井的最小水平主应力变化量与最大水平主应力变化量之差、最小水平主应力变化率与最大水平主应力变化率之差分别绘制在图12中。对于寿阳区块，进行重复压裂的主要目的是为了压开新的裂缝，这就需要水平应力发生一定程度的偏转，使之与原地应力间具有一定的夹角。当目标井最大水平主应力的增长量或增长率大于等于最小水平主应力时，该井水平地应力不会产生足够的偏转，可暂时将这类井排除，需要优选最小水平主应力增长率与最大水平主应力增长率较大的井。由图12可以看出，Δσh／σh－ΔσH／σH的大小顺序为：SYNY-136＞SYNY-161＞SYNY-187＞SYNY-125＞SYNY-135＞SYNY-185＞SYNY-150＞SYNY-190＞SYNY-173。由此可见，SYNY-136井地应力偏转程度最高，二次压裂时容易产生新裂缝；而SYNY-173井地应力偏转程度最低，二次压裂时不易产生新裂缝。

图12 寿阳南燕竹A2区块9口井地应力变化结果Fig.12 Stress change results of nine wells in NYZ-A2 block in Shouyang

3 结论

1）利用CMG和ABAQUS建立了渗流-应力耦合四维动态地应力模型，并用该模型反演得到的寿阳南燕竹A2区块初始地应力场与现场实测数据吻合度较高，表明该四维地应力建模分析方法可用于煤层气开采过程中的动态地应力分析。

2）通过对寿阳南燕竹A2区块单井地应力变化对比分析结果可以得到，相邻两井之间地应力变化趋势相近，表明该区块15＃煤层储层连通性较好，但产量或孔隙压力的变化分布强弱对地应力变化和储层沉降仍有一定影响。

3）A2区块9口在排井地应力变化强弱对比结果表明，SYNY-136井地应力偏转程度最高，二次压裂时容易产生新裂缝，而SYNY-173井地应力偏转程度最低，二次压裂时不易产生新裂缝。上述认识可为煤层气排采制度调整、加密井井壁稳定性分析、二次压裂优化设计等提供参考。