基于多场耦合效应的水平井压裂应力场分析

官 兵，李士斌，张立刚，陈双庆，王业强.

(1.东北石油大学石油工程学院，黑龙江大庆 163318；2.大庆油田有限责任公司第三采油厂，黑龙江大庆 163000)



基于多场耦合效应的水平井压裂应力场分析

官 兵1，李士斌1，张立刚1，陈双庆1，王业强2.

(1.东北石油大学石油工程学院，黑龙江大庆 163318；2.大庆油田有限责任公司第三采油厂，黑龙江大庆 163000)

为了有效分析水平井压裂过程中的应力场分布情况，考虑应力-渗流-温度多物理场因素的耦合作用，利用有限元数值模拟手段，开展了对水平井压裂裂缝近井应力场的模拟分析。首先通过岩石物理力学参数测定试验，得到了有限元分析的基本数据；之后引入应力-渗流-温度耦合理论，建立了水平井扰动应力场的有限元分析模型，模拟计算得到研究区的压力、温度和应力分布情况；最后，综合多场耦合效应得到近水平井筒裂缝周围的应力场变化规律，确定裂缝在平行井筒方向对其周围应力场的扰动范围主要集中在1倍的裂缝长度以内，裂缝产生的诱导应力降低了原井筒附近的压应力大小。将数值解结果与已有扰动应力数学模型计算的相应解析解进行对比，验证了有限元模型求解的合理性，为水平井压裂应力场的分析研究提供了一定的基础和参考。

水平井压裂；多场耦合；应力场；有限元；合理性

水平井压裂应力场是影响裂缝起裂和扩展行为的关键因素之一[1-2]，而压裂过程中应力场又受控于多个物理场的共同作用。水平井压裂改造过程中，井筒进入储层，打破原始地层的力学平衡，在原地应力场的约束和井筒内液柱压力的共同作用下，井筒周围产生应力集中效应，压裂裂缝产生诱导应力，二者叠加使地层岩石发生变形，从而扰动原始地应力场。同时由于注入压力与地层孔隙压力的差异造成压裂液渗滤引起附加应力，加之注入液体温度与储层温度的差异也会产生热应力作用。因而可知，水平井压裂是一个应力-渗流-温度多场耦合的过程。

早在1995年，陈勉[3]等人就基于多孔弹性理论、叠加原理，建立了考虑内压、地应力及压裂液渗流效应的大斜度井筒周围应力分布计算模型；2004年，徐严波[4]在此基础上引入地层温度产生的热应力，建立了水平井筒周围的应力分布模型；之后刘洪[5]、尹建[6]、翁定为[7]等都对人工裂缝的诱导应力场进行了研究，推导出了它的计算数学模型，其中尹建[8]、翁定为等人还利用有限元数值模拟方法分析了压裂裂缝周围的应力场变化。基于这些研究，已经建立了比较成熟的水平井井筒及压裂裂缝所产生的诱导应力场计算模型，并且提出了多场耦合作用的影响，但是对于同时考虑流-固-热耦合理论和水平井与裂缝共同形成的扰动应力场研究却较少。本文在前人研究的基础上，针对考虑应力-渗流-温度多场作用各种因素影响下的水平井压裂扰动应力场计算模型开展了进一步的研究，提出有限元数值分析方法，模拟计算了压裂过程中多物理场耦合作用下的应力场变化，并将其与解析解对比进行了合理性验证。

1 试验数据

选取大庆油田外围的扶杨油层致密砂岩组作为本文的研究对象，遵循中华人民共和国地质矿产行业标准(DZ/T0276.4—2015)《岩石物理力学性质试验规程》，挑选完整且表面没有明显裂纹的岩石制作成直径×高度规格为25 mm×50 mm的圆柱体岩心若干，记录好取样深度和对应标号。通过开展孔隙度、渗透率测试及三轴应力试验等获得了研究区储层岩石的基本物性参数和力学参数(表1)，为储层裂缝数值模拟分析提供了基础。

表1 取样岩心试验数据结果

2 多场耦合理论

多场耦合作用是指一个系统中所存在的物理场相互影响，区别于单一的物理作用，不同物理场中每一个都在其他物理场中扮演着不同的角色，始终存在耦合关系，包括其起始状态和整个过程。分析一个复杂多物理场的前提是引入耦合关系。基于国内外大量学者[9-11]对不同物理场耦合关系数学模型的分析，推导出了一组考虑多物理场完全耦合而又具有可求解性的方程组。其中，假设以拉应力为正，孔隙水压力为负，得到应力-渗流-温度耦合模型的完整方程组如下：

(1)

式中G——岩石剪切模量，Pa；ν——地层泊松比；p——地层孔隙水压力，Pa；T——地层温度，℃；T0——无应力状态下的绝对温度，℃；t——时间，s；K′——排水体积模量，且K′=2G(1+ν)/3(1-2ν)，Pa；

Ks——岩石骨架的有效体积模量，Pa；

αp——地层Biot系数，αp=1-K′/Ks，αp≤1，其值取决于地层岩石压缩性；

αT——岩体热膨胀系数，℃-1；

al——流体热膨胀系数，℃-1；

as——岩石骨架的体积膨胀系数，℃-1；

Fi——i方向的体积力，N；

ui(i=x,y,z)——i方向的位移，m；

εv——体积应变，εv=εxx+εyy+εzz；

Φ——孔隙度，%；

κ——岩石固有渗透系数，m2；

Vl——多孔介质中流体速度矢量，s-1；

βl——流体的体积模量，Pa；

ρl——流体密度，kg·m-3；

γ——γ=(2μ+3λ)β，μ和λ为拉梅常数，J/m3，β为各向同性固体的线性热膨胀系数，℃-1；

g——重力加速度，m·s-2；

z——垂直坐标，m；

(ρcp)t——充满了流体的多孔介质的比热容，J/(m3·℃)；

Kt——充满了流体的多孔介质的热传导系数，W/(m3·℃)；

qt——充满了流体的多孔介质的热源强度，W/m3。

由于应力-渗流-温度多场耦合方程组既考虑了不同物理场的相互作用又是时间的物理过程，有典型的非线性，除少数特殊情况可求出解析解外，一般需要用数值方法进行求解。基于应力-渗流-温度多场耦合理论模型的建立，本文选用专门针对偏微分方程组求解的有限元分析软件COMSOLMultiphysics对其进行求解。

3 有限元模型求解

3.1 模型建立

利用有限元数值分析方法分析水平井压裂裂缝扰动应力场，首先需要建立扰动应力场的三维几何模型。设定地层模型的长宽高尺寸分别为2000m×1000m×1000m；井筒半径为5m，水平井筒长1000m，假设水平井筒沿最小水平主应力方向延伸，裂缝与水平井筒夹角为90°；压裂人工裂缝形态设定为椭球型，其半长为100m，半高为50m，开度为5m。绘制其几何模型如图1所示。其中，直角坐标系(x,y,z)与地层岩石所受的3个主应力(σh,σH,σv)方向相对应。几何模型建立以后，需要将其转化为有限元数学模型。有限元模拟数学模型的建立首先要按照有限元数值分析所要求的数学、力学规则进行单元划分。本模型采用自由四面体单元对模型进行剖分，为了提高计算精度，在剖分井筒和裂缝单元时需要进行网格加密处理，最后得到研究区块的扰动应力场有限元模型如图2所示。

图1 扰动应力场的三维几何模型示意图Fig.1 3D geometry model of disturbed stress field

通过COMSOL有限元数值分析软件，在建立有限元模型之后，只要确定储层岩石的物理力学参数，并给予模型合理的施加载荷和边界约束，就能够对水平井压裂扰动应力场的分布进行模拟计算。

图2 扰动应力场的有限元模型示意图Fig.2 Finite element model of disturbed stress field

3.2 基本数据

假设地层岩石为多孔弹性介质，其中的流体流动满足Darcy定律，井筒与裂缝中充满单相液流，井筒与裂缝完全连通，井壁与裂缝面所受压力和温度等同于其中液体的压力和温度，注入液体比地层温度低24℃，其他给定分析水平井压裂扰动应力场分布的模型计算基本数据见表2。

3.3 求解结果

根据多场耦合理论模型及地层岩石的物理力学参数，对构建的水平井压裂裂缝有限元分析模型进行计算，得到了水平井压裂裂缝的应力-渗流-温度多场耦合的模拟计算结果。

表2 模型计算基本数据

图3是水平井压裂过程中的流场压力作用云图。如图3a所示，地层整体的压力分布在20 MPa左右，井筒和裂缝内的流体压力较高为30 MPa，而地层多孔介质具有渗透作用，井筒和裂缝周围会发生一定的液体渗滤作用(其流动方向如图中流线方向所示)，因而也会产生压力传递，从而出现压力过渡阶段(其作用方向如图中红色箭头方向，呈发散状)。截取水平井筒轴线所在的二维水平面，分析水平井筒和压裂裂缝周围的局部压力变化，如图3b所示。井壁和裂缝面处的压力最高为30 MPa，而随着沿箭头方向距离的增加，压力逐渐降低，直至接近原地层孔隙压力20 MPa。

图3 流场压力计算云图及等值线图Fig.3 Flow field pressure result and contour map

图4是水平井压裂过程中的温度变化云图。如图4a所示，地层整体温度随着深度的增加而逐渐升高，其分布范围为20～59.903℃。由于井筒内注入初始液体温度较低，假设井筒和裂缝壁面处的温度与地面温度一致为20℃。同样截取二维水平面，如图4b所示。水平井筒所在深度的原地层温度分布在44 MPa左右，地层温度高于井筒内的液体温度，从而产生传热过程，地层温度向井筒方向传递，使井筒和裂缝周围的温度逐渐升高。

图4 温度场的计算云图及等值线图Fig.4 Temperature field result and contour map

水平井筒压裂的流场压力和温度场变化都会影响其周围的应力分布，同时考虑原始地层的应力场作用，可以得到多场耦合作用下的水平井压裂过程中的地层水平主应力的分布变化情况，如图5和图6所示。其中压应力为负，拉应力为正。

图5 最小水平主应力计算云图及等值线图Fig.5 Minimum horizontal principal stress result and contour map

如图5a所示，地层整体的最小水平主应力方向没有发生明显偏转，依旧平行于x轴方向；地层整体的最小水平主应力大小集中在30 MPa左右，与原地应力场设置一致，只在裂缝周围发生了显著变化。截取水平井筒轴线所在的二维水平面，选取压裂裂缝周围的局部应力场进行分析，如图5b所示。从图中可以看出水平井筒和裂缝的出现使原地应力场的应力发生了扰动，裂缝周围应力场发生了一定的偏转，其扰动范围为-37.536～6.5907 MPa；而只有在裂缝壁面附近出现了正的拉应力，其应力值为6.6 MPa左右；随着距离裂缝壁面水平距离x轴方向的增加，其周围应力由拉应力变为压应力而逐渐增大，裂缝尖端出现较大的应力集中，为-37 MPa左右。

图6 最大水平主应力计算云图及等值线图Fig.6 Maximum horizontal principal stress result and contour map

如图6a所示，地层整体的最大水平主应力方向没有发生明显偏转，依旧平行于y轴方向；地层整体的最大水平主应力大小集中在50 MPa左右，与原地应力场设置一致，与最小水平主应力不同的是，水平井井筒和裂缝周围都发生了较为明显的变化。同样截取水平井筒轴线所在的二维水平面，分析压裂裂缝周围的局部应力场，如图6b所示。裂缝周围应力场发生了一定的偏转，而井筒和裂缝同时作用，使最大水平主应力发生扰动，其扰动范围为-97.062～-29.388 MPa；裂缝壁面附近出现了正的拉应力，与原地应力场的负压应力叠加，其值为-29 MPa左右；随着距离裂缝壁面x轴方向水平距离的增加，其周围压应力逐渐增大，裂缝尖端出现较大的应力集中，为-97 MPa左右，井筒周围的应力变化主要分布在y轴方向，其变化规律同裂缝。

通过以上对水平井压裂裂缝应力场、渗流场和温度场的分析可以看出，水平井压裂裂缝会使井壁附近的原始地层条件发生扰动，它是一个多场耦合作用的过程，因而在分析其压裂应力场扰动规律的同时，必须考虑其他物理场的作用效果。

3.4 结果分析

有限元数值分析方法能够快速高效地对多场耦合方程组进行求解，但是针对水平井分段压裂扰动应力场的计算分析，有限元分析是否正确合理还需要展开进一步研究。以大庆外围扶杨油层某水平井分段压裂的水平井作为研究对象，提取水平井筒井壁所在直线(x轴方向)上的计算应力值，绘制出扰动应力变化规律的曲线。通过将有限元数值模拟计算的数值结果与已有的成熟扰动应力数学模型[5,12-13]计算的解析结果进行对比(图7)，验证有限元模型求解的合理性。

图7 井壁处扰动应力数值解与解析解对比Fig.7 Comparison between the numerical and analytical solutions of disturbance stress in well wall

从图7中可以看出，针对水平井压裂扰动应力场分析这一问题，利用有限元数值模拟软件得到的结果曲线在局部范围内会出现一定的波动，这是由于在进行三维建模求解过程中网格的精度划分所造成的。而受限于数值模拟软件和计算机硬件配置，由网格精度所引起的计算误差无法避免，但不影响其对应力大小随距离变化规律的分析。在平行井筒方向，裂缝对其周围的应力场扰动范围主要分布在1倍的裂缝长度以内，扰动后的最小水平主应力下降速率先缓后快，当其与裂缝面x轴方向的距离接近1倍裂缝长度后，裂缝产生的扰动逐渐消失；而扰动后的最大主应力先快速降至最低，在其0.8倍裂缝半长处出现最低值，之后又缓慢上升，直至其影响程度趋于0。

通过与已有数学模型所计算的解析解结果进行对比，两种方法所得到的应力变化规律基本一致，其结果基本相符。因此，本文利用有限元数值分析方法展开对水平井压裂扰动应力分析是合理的。

4 结论

(1)通过对研究区致密砂岩进行的密度、孔隙度、渗透率基本物性测试，以及三轴应力试验，掌握了该储层岩石的物性特征：孔隙度平均值为4.70%；渗透率平均值小于0.5 mD，分布在纳达西级别；岩石杨氏模量平均值为20387.1 MPa，泊松比平均值为0.130，抗压强度平均值为213.54 MPa。

(2)基于应力-渗流-温度多场耦合方程，利用有限元数值分析软件COMSOL Multiphysics建立了水平井分段压裂扰动应力场的有限元分析模型，提出了求解多物理场的有限元数值方法，并计算了裂缝周围的压力、温度和应力分布情况。

(3)针对水平井分段压裂问题，选取相同研究对象，将有限元分析方法模拟计算得到的数值解与已有水平井压裂扰动应力数学模型计算的解析解进行对比，得到的求解结果基本一致，验证了有限元数值法研究水平井压裂扰动应力场的合理性。

[1] MANCHANDA R, SHARMA M M. Impact of completion design on fracture complexity in horizontal wells[C]. SPE159899, 2012.

[2] ROUSSEL N P, MANCHANDA R, SHARMA M M. Implications of fracturing pressure data recorded during a horizontal completion on stage spacing design[C]. SPE152631, 2012.

[3] 陈勉,陈治喜,黄荣樽.大斜度井水压裂缝起裂研究[J].石油大学学报(自然科学版),1995,19(2):30-35.

[4] 徐严波.水平井水力压裂基础理论研究[D].成都:西南石油大学,2004:1-128.

[5] 刘洪,胡永全,赵金洲,等.重复压裂气井诱导应力场模拟研究[J].岩石力学与工程学报,2004,23(23):4022-4027.

[6] 尹建,郭建春,曾凡辉.水平井分段压裂射孔间距优化方法[J].石油钻探技术,2012,40(5):67-71.

[7] 翁定为,严星明,卢拥军,等.考虑应力干扰的致密油气压裂设计及实现方法[J].中国矿业大学学报,2014,43(4):639-645.

[8] 尹建.水平井分段压裂诱导应力场研究与应用[D].成都:西南石油大学,2014:1-75..

[9] 王自明,江同文,宋文杰,等.弹性气藏中气体渗流的气-固-热耦合数学模型[J].天然气工业,2003,23(3):85-87.

[10] 贾善坡,邹臣颂,王越之,等.基于热-流-固耦合模型的石油钻井施工过程数值分析[J].岩土力学,2012,33(S2):321-328.

[11] 曹伟,康志勤,吕义清,等.基于COMSOL Multiphysics的煤层气对流热采分析[J].地下空间与工程学报,2014,10(5):1139-1145.

[12] 邵尚奇,田守嶒,李根生,等.水平井缝网压裂裂缝间距的优化[J].石油钻探技术,2014,42(1):86-90.

[13] 曾顺鹏,张国强,韩家新,等.多裂缝应力阴影效应模型及水平井分段压裂优化设计[J].天然气工业,2015,35(3):55-59.

Analysis of Stress Field of Horizontal Well FracturingBased on Multi-Field Coupling Effect

Guan Bing1, Li Shibin1, Zhang Ligang1, Chen Shuangqing1, Wang Yeqiang2

(1.School of Petroleum Engineering, Northeast Petroleum University, Daqing, Heilongjiang 163318, China;2.No.3 Oil Production Plant of Daqing Oilfield Co., Ltd., Daqing, Heilongjiang 163000, China)

In order to analyze the distribution of stress field in the process of horizontal wells fracturing effectively, the multi-field coupling effect was considered and the finite element numerical simulation analysis was applied. Firstly, through the determination test of rock physical and mechanical parameters, the basic data were obtained. Then introducing the theory of stress - seepage - temperature coupling, the finite element analysis model of the disturbance stress field in horizontal well fracturing was established. The results of pressure, temperature and stress distribution in the study area were simulated and calculated. Finally, combined with the multi-field coupling effect, the change law of the stress field around the fractures was achieved. It was found that in the parallel wellbore direction the stress disturbance range of one fracture was mainly concentrated within 1 time of its length and its induced stress made the original compressive stress around the wellbore reduced. Comparing the numerical solution results with the corresponding analytic solutions of the existing disturbed stress mathematic model, the rationality of the finite element model was verified. It can provide a certain basis and reference to study the stress field in horizontal wells fracturing.

horizontal well fracturing; multi-field coupling; stress field; finite element; rationality

东北石油大学研究生创新科研项目“致密油气藏体积压裂局部应力场改造控制机理研究”(YJSCX2015-017NEPU)、国家自然科学基金项目“致密储层井工厂压裂应力场演化机理及调控方法”(51504067)联合资助。

官兵(1990—)，女，黑龙江省双城市人，2013年毕业于东北石油大学石油工程专业，在读博士研究生，主要从事非常规压裂裂缝应力场数值模拟分析的研究工作。邮箱：sygcgb1990@163.com.

TE357

A