基于粘结单元法的干热岩储层水力压裂数值模拟与参数优化

摘要

水力压裂技术是实现低渗油气及地热储层的高效开发利用的关键技术手段，为了研究干热岩型地热储层水力压裂过程中水力裂缝的扩展规律，本文使用粘结单元法（Cohesive Zone Method， CZM）研究了压裂液排量、压裂液粘度以及水平地应力差对水力裂缝形态的影响，并利用正交试验对上述压裂工艺参数的组合进行优化。结果表明：压裂液排量对水力裂缝的长度具有重要影响，而压裂液的粘度对水力裂缝的宽度具有显著影响；压裂液的排量和粘度的增加，促进了分支裂缝的萌生和扩展；水平地应力差为1 MPa时，本文所建立的模型在压裂液排量和粘度分别取0.004 m³/s和0.07 Pa·s条件下，可获得最佳的压裂改造效果；随着压裂液的排量和粘度的持续增加，当压裂液的排量和粘度分别超过0.004 m³/s和0.07 Pa·s后，继续增加压裂液的排量和粘度将导致水力裂缝的长度和宽度的减小，可见在实际压裂过程中不能盲目通过提高压裂液的排量和粘度的方式实现对压裂效果的持续改进。本文丰富了干热岩储层改造的数值模拟手段，相关研究成果有望为干热岩型地热资源开采过程中裂缝扩展行为预测和压裂工艺参数的优化提供技术支撑。

关键词

干热岩; 储层改造; 水力压裂; 裂缝扩展; 粘结单元法; 正交试验; 数值模拟

0 引言

随着全球经济快速发展，石油、天然气等能源的对外依存度逐年增加，能源危机已成为世界各国需要关注的迫切问题之一^［1-2］。为了保障国家能源安全，我国提出了构建煤、油、气、核、新能源和可再生能源多元化的能源供应体系，形成常规与非常规能源共同发展的新型勘探、开发模式^［3］。

地热能因其分布广泛、低碳环保和可持续开采等特点，现已成为世界各国能源领域学者重点关注的新能源领域。与其他清洁能源如风能、太阳能相比，季节更替、气候变化以及昼夜循环等干扰因素对地热能提取的影响较小，其发电效率可达73%^［4-5］。当前，地热资源主要包括浅层地热资源、水热型地热资源和干热岩型地热资源。。

干热岩赋存在地下3000 m以深，其内部没有或仅含有极少量热蒸汽的高温低渗透性岩体（温度>180℃）^［6-7］。干热岩分布广泛，现有资料表明我国干热岩型地热资源储量丰富，资源量为2.52×10²⁵ J，埋深3~10 km范围内的资源量折合8.56×10¹⁴ t标准煤^［8-9］。然而，干热岩主要由花岗岩、花岗闪长岩等具有极低孔隙度和渗透率的岩浆岩和变质岩组成，依靠岩石基质和热储内部的构造裂隙的渗流能力获得具有商业开采价值的高温流体。因此，提出通过储层改造的方式提高热储内部导流能力，以实现干热岩型热储内部地热能的有效提取，这种通过储层改造开采干热岩的系统工程称为增强型地热系统（Enhanced Geothermal System， EGS）^［10］。

目前，EGS工程仍处于示范性工程的预研和可行性研究阶段，地质勘查技术、深部钻完井技术、储层改造技术和高效热能提取技术是EGS实现商业化开采需要突破的关键技术问题。其中，提高干热岩型热储渗透特性的储层改造技术是亟待关键突破的重点研究领域。

水力压裂是储层改造技术的核心技术手段，各国EGS示范工程中，水力压裂技术也在不断的创新和发展。EGS建储过程中的主要压裂形式包括沿用了传统石油行业的支撑剂型压裂和低粘度压裂流体形成以剪切压裂为主的水力压裂方式^［11］。水力压裂是地应力-渗流耦合的复杂物理过程，而干热岩水力压裂则是涉及地应力-渗流-温度耦合的极为复杂的物理过程^［12-13］。目前该领域研究主要通过理论、数值模拟和室内试验的方式复现深部热储内裂缝的形成、延伸和扩展规律^［14-16］。

数值模拟研究因其可重复性高，成本较低等特点逐渐被学者关注，可用于水力压裂研究的数值模拟方法主要包括位移不连续法、离散元法、扩展有限元法和有限元法。（1）位移不连续法：Chuprakov等^［17］使用位移不连续法并结合Mohr-Coulomb准则研究构造裂隙对水力裂缝扩展行为的影响。Zhao等^［18］建立三维压裂模型，并研究多裂缝同时扩展条件下的三维裂缝形态。该方法能够模拟复杂裂缝扩展问题，但在处理非均匀介质中裂缝扩展问题时存在较大难度。（2）离散元法：Abuaisha等^［19］使用该方法研究了近井区域裂缝在地应力作用下的起裂与扩展行为。Lanru等^［20］研究了地应力作用下水力裂缝的变形以及内部压裂液的流动模式。该方法通过块体或者颗粒来代表岩体，虽然能够模拟水力裂缝沿任意路径的扩展行为，但由于块体和颗粒之间的破坏通过键能的破坏来表示，其破坏过程中的接触刚度和强度参数变化需要不断的迭代计算和接触检查，导致计算量较大。（3）扩展有限元法：Wang等^［21］利用扩展有限元法研究页岩各向异性的物性条件对水力裂缝扩展行为的影响。孙可明等^［22］研究了储层内部层理对水力裂缝扩展的影响以及裂缝形态的演化法则。王涛等^［23］利用该方法研究了三维水力裂缝的扩展模拟。但是该方法在模拟水力裂缝交叉行为过程中由于交汇处节点的畸变导致数值模拟计算结果与实际结果产生较大误差。

与上述模拟方法相比，本文采用有限元法框架内的粘结单元法（Cohesive Zone Method， CZM）对水力裂缝展开研究。该方法具有以下3方面优势：（1）适应本构模型复杂的岩石模型；（2）可模拟裂缝间的交汇行为；（3）能够水力压裂过程中岩石基质与裂缝系统之间的流-固耦合过程^［23］。许多学者使用该方法对天然裂缝影响^［24］、压裂流体特性^［25］、应力场变化^［26］等因素的影响展开研究。

本文基于自主编写的网格节点分离程序，在压裂模型内部嵌入零厚度的粘结单元，进而实现对不同压裂液参数和地应力差条件下水力裂缝起裂与扩展行为的模拟，以期获得不同压裂液参数和水平地应力差对压裂效果的影响，并最终得出适合于本文所建立模型的最佳压裂参数组合。同时，综合考虑压裂过程中主裂缝、分支裂缝的形成与扩展对模型内部缝网复杂性的影响，选取累计裂缝长度和最大裂缝宽度作为正交试验结果的评价指标，相比仅通过裂缝长度和宽度来定性描述压裂效果更为合理。

1 模型描述

1.1 基本假设

建立模型的假设条件如下：

（1）所选储层岩石为非均质分布，能够服从Weibull分布；

（2）储层岩石内部为各向同性，岩体为弹性介质，服从线弹性力学规律；

（3）岩体满足Biot固结理论和有效应力原理；

（4）岩体为单相饱和介质；

（5）水为液态，不考虑水的相变；

（6）岩体力学行为用连续介质损伤力学理论描述，满足最大拉应力准则和Mohr-Coulomb准则；

（7）忽略在压裂过程中压裂液与储层岩石发生的化学或物理作用而导致的岩石力学物性产生异变。

1.2 CZM的力学本构模型与单元内流体流动模式

1.2.1 CZM的本构模型

粘合单元法单元的损伤开始出现前需要满足的应力与应变关系如下述，单元承受的牵引应力矢量t的表达式如下所示：

$$t=\left\{\begin{array}{l}{t}_{\mathrm{n}}\\ t_{\mathrm{s}}\\ t_{\mathrm{t}}\end{array}\right\}=K_{\mathrm{\epsilon }}=\left[\begin{array}{l}{K}_{\mathrm{n}\mathrm{n}},K_{\mathrm{n}\mathrm{s}},K_{\mathrm{n}\mathrm{t}}\\ K_{\mathrm{s}\mathrm{n}},K_{\mathrm{s}\mathrm{s}},K_{\mathrm{s}\mathrm{t}}\\ K_{\mathrm{t}\mathrm{n}},K_{\mathrm{t}\mathrm{s}},K_{\mathrm{t}\mathrm{t}}\end{array}\right]\left\{\begin{array}{l}{\epsilon }_{\mathrm{n}}\\ {\epsilon }_{\mathrm{s}}\\ {\epsilon }_{\mathrm{t}}\end{array}\right\}$$

（1）

式中：t_n——单元的法向方向的牵引力；t_s、t_t——分别为两个切线方向的牵引力；${{\rm K}}_{\mathrm{n}\mathrm{n}}={{\rm E}}_{\mathrm{n}\mathrm{n}}/T_{\mathrm{0}}$、${{\rm K}}_{\mathrm{s}\mathrm{s}}={\mathrm{{\rm E}}}_{\mathrm{s}\mathrm{s}}/T_{\mathrm{0}}$、${{\rm K}}_{\mathrm{t}\mathrm{t}}={{\rm E}}_{\mathrm{t}\mathrm{t}}/T_{\mathrm{0}}$分别为牵引分离刚度；$T_{\mathrm{0}}$——CZM单元的初始厚度；$E_{\mathrm{n}\mathrm{n}}$——粘结单元的法向方向的弹性模量；${{\rm E}}_{\mathrm{s}\mathrm{s}}$、${{\rm E}}_{\mathrm{t}\mathrm{t}}$——分别为粘结单元的两个切线方向的弹性模量；${\epsilon }_{\mathrm{n}}=$${\delta }_{\mathrm{n}}/T_{\mathrm{0}}$——与牵引力相对应的法向方向的分离应变；${\epsilon }_{\mathrm{s}}={\delta }_{\mathrm{s}}/T_{\mathrm{0}}$、${\epsilon }_{\mathrm{t}}={\delta }_{\mathrm{s}}/T_{\mathrm{0}}$——分别为相对应的两个切线方向的分离应变；${\delta }_{\mathrm{n}}$——法向方向的分离位移；${\delta }_{\mathrm{s}}\mathrm{、}{\delta }_{\mathrm{t}}$——分别为两个切线方向的分离位移。

基于粘弹塑性损伤模型CZM对材料的损伤开裂行为进行模拟时，第一必须深刻了解粘结单元的几何厚度和本构厚度，因为在对本构方程中刚度和应变计算时，其结果与本构厚度的取值有直接关系。第二为了保持计算过程简便，建立的此模型中本构厚度的取值默认为1，这样可以使粘结单元的应变数值与对应方向的位移保持一致。

1.2.2 CZM单元内部流体流动模型

为了达到使模型单元内部流体保持连续流动的目的，CZM的孔压-渗流单元将损伤区内的流体流动分为法向流动和切向流动，CZM孔压-渗流单元内的流体流动进一步分为法向流动和切向流动，如图1所示。

图1  CZM孔压-渗流单元内的流体流动示意

Fig.1  Schematic diagram of fluid flow in the pore pressure‑permeability unit in CZM

1.2.2.1 CZM单元内流体的法向流动

根据达西渗流定律，CZM裂隙单元内的流体在单元表面的法向流动与其接触介质单元的体积速率有直接关系，所建模型的周围岩石基质在现实压裂过程中的滤失行为可用它来进行模拟。

流体在单元上表面和下表面法向流动的计算公式为：

$$\left\{\begin{array}{l}{q}_{\mathrm{t}}=C_{\mathrm{t}}\left(p_{\mathrm{c}\mathrm{z}\mathrm{m}}-p_{\mathrm{t}}\right)\\ q_{\mathrm{b}}=C_{\mathrm{b}}\left(p_{\mathrm{c}\mathrm{z}\mathrm{m}}-p_{\mathrm{b}}\right)\end{array}\right.$$

（2）

式中：$q_{\mathrm{t}}$——CZM裂缝单元内压裂流体沿着上单元表面的体积流速，m³/s；$q_{\mathrm{b}}$——CZM裂缝单元内压裂流体沿着下单元表面的体积流速，m³/s；$C_{\mathrm{t}}$——相应单元的上表面滤失层的滤失系数，m/（Pa∙s）；$C_{\mathrm{b}}$——相应单元的下表面滤失层的滤失系数，m/（Pa∙s）；$P_{\mathrm{t}}$——与单元上表面接触部分的岩体孔隙压力，MPa；$P_{\mathrm{b}}$——与单元下表面接触部分的岩体孔隙压力，MPa；$P_{\mathrm{c}\mathrm{z}\mathrm{m}}$——CZM裂缝单元内压裂流体的压力，MPa。

1.2.2.2 CZM单元内流体的切向流动

假设单元内部压裂流体为不可压缩的牛顿流体，基于Poiseuille法则，单元内压裂流体的切向流动服从润滑方程：

$$q=-d_{\mathrm{c}\mathrm{z}\mathrm{m}}^{\mathrm{3}}/\mathrm{12}\mu \nabla p_{\mathrm{c}\mathrm{z}\mathrm{m}}$$

（3）

式中：$q$——通过CZM裂隙单元横截面的切向体积流量，m³/s；$d_{\mathrm{c}\mathrm{z}\mathrm{m}}$——裂隙单元随压裂时间而产生的不断变化的裂缝面的张开位移，m；$\mu$——裂隙内压裂流体的粘度，Pa·s；$\nabla p_{\mathrm{c}\mathrm{z}\mathrm{m}}$——CZM单元切向方向的压裂流体压降梯度，Pa/m。

假设CZM裂隙单元内流体符合两块平行板之间的层流流动的模式进行流动，可对单元内的压降梯度$\nabla p_{\mathrm{c}\mathrm{z}\mathrm{m}}$进行计算：

$$d_{\nabla p_{\mathrm{c}\mathrm{z}\mathrm{m}}}/$$$$d_x=-\mathrm{12}\mu q/(h_{\mathrm{f}}{d}_{\mathrm{c}\mathrm{z}\mathrm{m}}^{\mathrm{3}})$$

（4）

对上式进行积分可得CZM单元内压裂流体的压降方程：

$$P_{\mathrm{0}}=P{}_{\mathrm{i}}{}^{}+\mathrm{12}\mu qx/(h_{\mathrm{f}}{d}_{\mathrm{c}\mathrm{z}\mathrm{m}}^{\mathrm{3}})$$

（5）

式中：$P_{\mathrm{0}}$——压裂流体流入CZM裂隙单元时的流体压力，MPa；$P_i$——压裂流体沿切向方向穿过CZM裂隙单元后的流体压力，MPa；$x$——CZM裂隙单元切向方向的单元长度，m；$h_{\mathrm{f}}$——裂缝的高度，m；$d_{\mathrm{c}\mathrm{z}\mathrm{m}}$——裂隙单元随压裂时间而产生的不断变化的裂缝面的张开位移，m；$q$——裂缝内的流体流量，m³/s。

综上所述，基于CZM建立的如上模型可以实现水力压裂过程中缝内流体与周围岩石介质的多场耦合过程的模拟。

1.3 模型建立

本文利用自主编写的网格节点分离程序，采用ABAQUS软件在有限元计算框架内提供的基于粘弹塑性损伤模型的CZM对干热岩水力压裂裂缝的起裂与扩展行为进行模拟。

所建立的干热岩水力压裂数值模型如图2所示，该模型中部的红色圆点为压裂流体的注入点。

图2  干热岩水力压裂数值模型

Fig.2  Numerical modeling of hydraulic fracturing in hot dry rock

1.4 模拟参数选取

国内外学者针对花岗岩力学参数进行了大量的室内与原位试验^［27-28］。此外，干热岩水力压裂数值模拟方面也进行了大量的研究，结合干热岩数值模拟采用的物理力学参数^［29-30］，本文模型参数的选取如表1所示。

表1  花岗岩水力压裂数值模拟参数取值

Table 1  Mechanical parameters of hydraulic fracturing in granite

材料	参数及单位	取值
花岗岩基质	弹性模量/GPa	40
	泊松比	0.23
	渗透率/m2	1.5$\times$10-18
	孔隙压力/MPa	0
	孔隙比	0.05
CZM单元	弹性模量/GPa	40
	法向抗拉强度tn/MPa	10
	破裂面I方向的抗剪强度ts/MPa	20
	破裂面Ⅱ方向的抗剪强度tt/MPa	20
原位地应力	垂向地应力σv/MPa	10
	最大水平地力σH/MPa	9
	最小水平地力σh/MPa	4
流体性质	压裂液粘度/（mPa·s）	1/60/150
	压裂流体比重/(N·m-3)	9800/11760/14700

2 基于正交试验的压裂参数优化

2.1 正交试验基本参数

近年来，许多研究人员基于正交试验的原理对试验中的各影响因素做分析，能够对试验方案进行优化而且有效减少试验周期和成本，提高科研效率。

在进行正交试验前，需要根据试验的研究目标选择合适的试验因素和结果评价指标。本文正交试验将压裂过程中的压裂液排量、压裂液粘度以及水平地应力差作为试验因素，并对上述试验因素在水力压裂裂缝扩展行为和缝网形态的影响进行研究。

考虑到水力压裂过程中，模型内部除了形成主水力裂缝之外，还发育有分支裂缝，因此，为了准确评价水力压裂效果，本文选择累计水力裂缝长度和最大水力裂缝宽度作为水力压裂效果的评价指标，评价指标和模型内部的裂缝复杂性呈正相关，相应的评价指标越大表明储层改造效果越好。

本文使用极差分析法对正交试验结果进行分析，研究压裂液排量、压裂液粘度和水平地应力差对水力压裂效果的影响。试验因素水平见表2。

表2  正交试验因素水平

Table 2  Factors and levels of the orthogonal experimental design

因素 水平	A	B	C
	压裂液排量/ （m3·s-1)	压裂液粘度/ （Pa·s）	水平地应力差/MPa
1	0.001	0.001	1
2	0.002	0.060	3
3	0.003	0.150	5

2.2 正交试验设计

根据上述正交试验因素水平表设计正交试验，其中，基于自主编写的子程序，对模型压裂后的裂缝单元的位移信息进行自动提取和处理，进而获得正交试验中累计裂缝长度和最大裂缝宽度的值，如表3所示。

表3  正交试验设计方案

Table 3  Scheme for the orthogonal experimental design

因素 试验号	A	误差列	B	C	模拟方案	评价指标
	压裂液排量/（m3·s-1)		压裂液粘度/Pa·s	水平地应力差/MPa		累计裂缝长度/m	最大裂缝宽度/mm
1	1	1	1	1	A1B1C1	22.482	3.301
2	1	2	2	2	A1B2C2	11.509	3.251
3	1	3	3	3	A1B3C3	10.641	3.269
4	2	1	2	3	A2B2C3	12.093	3.237
5	2	2	3	1	A2B3C1	14.583	3.280
6	2	3	1	2	A2B1C2	11.410	3.015
7	3	1	3	2	A3B3C2	11.693	3.248
8	3	2	1	3	A3B1C3	15.475	3.095
9	3	3	2	1	A3B2C1	52.904	3.656

2.3 正交试验结果极差法分析

累计裂缝长度的极差分析结果如表4所示，最大裂缝宽度的极差分析结果如表5所示。

表4  累计裂缝长度的极差分析

Table 4  Range analysis of cumulative fracture length

	A	误差列	B	C
	压裂液排量/（m3·s-1)		压裂液粘度/（Pa·s）	水平地应力差/MPa
K1	44.632	46.268	49.367	89.969
K2	38.086	41.567	76.506	34.612
K3	80.072	74.955	36.917	38.209
k1	14.877	15.423	16.456	29.989
k2	12.695	13.856	25.502	11.537
k3	26.691	24.985	12.306	12.736
极差R	13.996	11.129	13.196	18.452
因素主次顺序	C>A>B
优选方案	C1A3B2

表5  最大裂缝宽度的极差分析

Table 5  Range analysis of maximum fracture width

	A		B	C
	压裂液排量/（m3·s-1)	误差列	压裂液粘度/（Pa·s)	水平地应力差/MPa
K1	9.821	9.786	9.411	10.237
K2	9.532	9.626	10.144	9.514
K3	9.999	9.940	9.797	9.601
k1	3.274	3.262	3.137	3.412
k2	3.177	3.209	3.381	3.171
k3	3.333	3.313	3.266	3.200
极差R	0.156	0.104	0.244	0.241
因素主次顺序	B>C>A
优选方案	B2C1A3

表中K_i为正交试验方案中各试验因素所对应的水平为i的结果之和。例如，表4中，B因素所在列的K₂=11.509+12.093+52.904=76.506。同理，可求得评价指标极差分析中的其余Ki的值。

表中$k_i$=$K_i/$$S$，表示所得$K_i$试验结果的算术平均值，$S$表示正交试验各试验因素水平的出现的次数，因此，表4中B因素所在的第三列中，$k_{\mathrm{1}}$=$\mathrm{49.367}/\mathrm{3}$，同理可以计算出所有$k_i$值。

$R$为极差，表中任何一列上的极差$R=\mathrm{m}\mathrm{a}\mathrm{x}\left\{k_{\mathrm{1}},k_{\mathrm{2}},k_{\mathrm{3}}\right\}-\mathrm{m}\mathrm{i}\mathrm{n}\left\{k_{\mathrm{1}},k_{\mathrm{2}},k_{\mathrm{3}}\right\}$。因此，表4中，水平地应力差（因素C）的$R=$29.989-11.537=18.452，同理，可计算得出表4和表5其余试验因素的极差。

基于正交试验的极差分析法原理，正交试验因素水平的变化对本次评价指标的影响与计算所得的极差呈正相关。因此通过表中计算得出的极差大小对各试验因素的影响程度进行排序并确定初选方案。

本文选取累计裂缝长度和最大裂缝宽度作为正交试验的评价指标，因此，应该选择各试验因素的水平值使得各评价指标的对应的k₁、k₂、k₃取得最大的值。表4中针对累计裂缝长度的正交试验中，根据极差R的大小对各因素影响的敏感性进行排序可得：C因素（水平地应力差）>A因素（压裂液排量）>B因素（压裂液粘度），因此，可得水平地应力差对压裂裂缝的累计裂缝长度具有重要影响，其次为压裂液排量，压裂液黏度对累计裂缝长度的影响最小。同时，C因素（水平地应力差）：k₁>k₃>k₂；A因素（压裂液排量）k₃>k₁>k₂；B因素（压裂液粘度）k₂>k₁>k₃。综上所述，对于累计裂缝长度的初步优化方案为C₁A₃B₂。同理，根据表5可得针对最大裂缝宽度的正交试验初步优化方案也为C₁A₃B₂。这表明水平地应力差为1 MPa条件下，压裂液排量和黏度分别选取0.003 m³/s和0.06 Pa·s，可获得表2中各因素水平范围内的最佳压裂效果。

3 基于正交试验结果的单因素分析

上述通过正交试验获得的C₁A₃B₂仅是表2中各因素水平范围内的初步优化参数组合。为了获取本文所建立模型的最佳压裂参数组合，需要对各试验因素展开单因素分析。由于水平地应力差在实际压裂工况下属于人为不可控因素，因此，本文的单因素分析主要围绕压裂液的排量和粘度展开研究。

3.1 压裂液排量对压裂效果影响的单因素分析

基于已获得的初步优化压裂参数C₁A₃B₂，在保持水平地压力差为1 MPa，压裂液粘度为0.06 Pa·s的基础上，将压裂液排量从0.003增加至0.005 m³/s，不同压裂液排量作用下模型内部的裂缝形态见图3。

图3  不同压裂液排量条件下模拟地层内部的水力裂缝形态

Fig.3  The morphology of hydraulic fracture under different fracturing fluid flow rate within the stimulated strata

根据模型内部水力裂缝形态的分析可得，水力裂缝的最大裂缝宽度与压裂液的排量呈正相关，且随着压裂液排量的增加，模型内部逐渐出现分支裂缝。

经过对数值模拟结果的后处理可得模型内部累计裂缝长度和最大裂缝宽度在不同压裂液排量作用下的模拟结果，如表6所示。

表6  不同压裂液排量下的累计裂缝长度和最大裂缝宽度

Table 6  Cumulative fracture length and maximum fracture width under different fracturing fluid flow rate

压裂液排量/ （m3·s-1）	累计裂缝长度/m	最大裂缝宽度/mm
0.003	22.631	3.196
0.004	29.005	3.369
0.005	20.741	3.329

根据表6可得，压裂液排量在0.004 m³/s时可取得累计裂缝长度和最大裂缝宽度的峰值，分别为29.005 m和3.369 mm。因此，在水平地压力差为1 MPa，压裂液粘度为0.06 Pa·s条件下，选取压裂液排量为0.004 m³/s可获得较好的压裂效果。

各排量作用下，缝内流体压力曲线见图4。

图4  不同压裂液排量下的压裂曲线

Fig.4  Pressure curves of fracturing fluids inside hydraulic fractures at different flow rate of fracturing fluid

根据图4可得，随着压裂液排量的增加，该模拟地层的起裂压裂逐渐升高，且压裂液排量越大，压裂曲线憋压阶段的增长速率越大。

综上所述，针对本文所建立的水力压裂模型中压裂液排量的单因素分析可得最优压裂方案组合为C₁B₂A_0.004，即水平地压力差为1 MPa，压裂液粘度为0.06 Pa·s条件下，压裂液排量为0.004 m³/s可获得较好的压裂效果。

3.2 针对压裂液粘度的单因素分析

同理，保持其他压裂参数不变的条件下，将压裂液粘度从0.05 Pa·s增加至0.08 Pa·s。不同压裂液粘度作用下模型内部的裂缝形态如图5所示。

图5  不同压裂液粘度条件下模拟地层内部的

Fig.5  The morphology of hydraulic fracture under different fracturing fluid viscosity within the stimulated strata

水力裂缝形态

分析模型内部水力裂缝形态可得，随着压裂液粘度增加，模型内部水力裂缝的最大宽度和缝网的复杂性也随之增加。但是，当压裂液粘度>0.07 Pa·s后，模型内部水力裂缝的宽度出现下降的趋势。

经过对数值模拟结果的后处理可得模型内部累计裂缝长度和最大裂缝宽度在不同压裂液粘度作用下的模拟结果，如表7所示。

表7  不同压裂液粘度下的累计裂缝长度和最大裂缝宽度

Table 7  Cumulative fracture length and maximum fracture width under different fracturing fluid viscosity

压裂液粘度/ （Pa·s）	累计裂缝长度/m	最大裂缝宽度/mm
0.05	18.603	3.157
0.06	12.108	3.237
0.07	26.961	3.344
0.08	18.890	3.191

由表7可得，当压裂液排量为0.004 m³/s条件下，累计裂缝长度在粘度为0.07 Pa·s时出现了峰值，分别为26.961 m和3.344 mm。因此，在水平地应力差为1 MPa，压裂液排量为0.004 m³/s条件下，选取压裂液粘度为0.07 Pa·s可获得较好的压裂效果。

不同压裂液粘度作用下缝内流体压力曲线见图6。

图6  不同压裂液粘度下的压裂曲线

Fig.6  Fracturing curves with different viscosity of fracturing fluid

根据图6可得，随着压裂液排量的增加，该模拟地层的起裂压裂逐渐升高，压裂曲线憋压阶段的增长速率与压裂液粘度呈正相关。

综上所述，针对本文所建立的水力压裂模型中压裂液排量的单因素分析可得最优压裂方案组合为C₁A_0.004B_0.07，即水平地压力差为1 MPa条件下，选取压裂液排量和粘度分别取0.004 m³/s和0.07 Pa·s，可获得本文所建模型的最佳压裂效果，即可获得累计裂缝长度和最大裂缝宽度的最大值。

4 结论

本文利用ABAQUS软件在有限元计算框架内提供的基于粘弹塑性损伤模型的Cohesive Zone Method（CZM）对干热岩模型内部水力压裂裂缝的扩展行为进行研究，并基于正交试验的原理对压裂液排量、压裂液粘度和水平地应力差进行分析，最终得到最优压裂参数组合方案。本文得到的主要结论如下：

（1）随着压裂液排量的增加，模型内部水力裂缝的宽度和长度随之增加，其对水力裂缝长度的影响较为明显；压裂液粘度与水力裂缝的长度与宽度也呈正相关，但其对裂缝宽度的影响更为显著。此外，压裂液排量和粘度的增加均能够促进分支裂缝的产生与扩展。

（2）基于极差分析与单因素分析，最终确定该模拟地层的最终优压裂方案为C₁A_0.004B_0.07，即水平地应力差是1 MPa条件下，压裂液排量选取0.004 m³/s，压裂液粘度选取0.07 Pa·s可获得最佳的压裂效果。在此条件下花岗岩模型内部的水力裂缝的最大裂缝宽度和累计裂缝长度均达到峰值，而且模型储层内部水力裂缝网络的形态也愈加复杂，分支裂缝产生并向模型内部不断扩展延伸，水力裂缝和分支裂缝的产生和扩展将极大提高模型内部裂缝网络的复杂性，提升储层改造效果。

（3）针对压裂液排量和粘度的单因素分析过程中，我们发现持续增加排量和粘度对累计裂缝长度和最大裂缝宽度的影响不是无限制增长的，因此，可以推断在实际压裂工况下，不能盲目依靠增加压裂液排量和粘度的方式实现对压裂改造效果的持续提升，应当充分了解和掌握待压裂地层的储层物性条件和不连续结构面的分布情况，在此基础上对压裂参数进行充分的优化设计，进而提高储层的压裂改造效果。