页岩气藏由于超低的孔隙度和渗透率，常规开发方式无法实现商业开发，只有通过水平井+大规模体积压裂的方式才能实现效益开采^[1]。大量的微地震监测结果表明，水平井大规模压裂后形成的人工缝网异常复杂，常规数值模拟方法在进行压后数值模拟时，只能人为的设定压裂裂缝为简单的双翼对称缝，或者在人工裂缝周围通过局部网格加密的方式圈定部分网格近似模拟储层改造体积^[2]，这些处理方法都掩盖了地层真实裂缝形态对气井产能的影响。因此，有必要开展页岩压后复杂缝网下全耦合数值模拟研究。

影响水力压裂缝网形态的参数众多，其中最主要的包括地应力、天然裂缝发育情况、水平井完井设计以及压裂泵注参数，而页岩气地质工程一体化压后数值模拟是目前唯一能够综合考虑上述因素影响的新技术。地质工程一体化技术是近几年在页岩气开发中新兴的一种理念，最早较为系统明确提出这个理念是在中国南方海相页岩气开发实践中^[3]。与北美页岩气相比，中国四川盆地海相页岩气特别是南川区块盆缘复杂构造带常压页岩气开发区在地质构造、储层物性、地表地貌、工程工艺、配套设施等方面存在诸多的挑战，地层压力系数低、单井产能低对常压页岩气的效益开发提出了更高的要求。历经几年的发展，地质工程一体化的理念逐渐被业界认同并积极推广，是实现中国页岩气跨越式发展的关键途径^[4-7]。

近年来，非常规储层研究技术（如储层精细地质建模、天然裂缝建模）的进步也为地质工程一体化的实施提供了技术支撑^[8]。页岩气地质工程一体化另一个关键结合点是储层改造。储层改造效果如何评价、水力裂缝如何展布、水力裂缝和天然裂缝如何相互作用等问题是地质研究和工程作业都非常关注的问题。对于水力裂缝的模拟，传统的方法是采用加密网格或根据微地震定性描述，但难以真实描述复杂水力裂缝网络的三维展布。水力压裂时多簇裂缝之间存在应力干扰从而导致裂缝延伸方向发生改变，中外学者对此进行了一些探讨^[9-11]，研究开发了基于 DDM 的水力压裂缝网数值模拟器，分析了多簇裂缝同步延伸的影响因素。其他学者在此基础上考虑了流体沿井筒的摩阻和射孔孔眼处的压降，提出了修正的 DDM 模型^[12]，同时众多研究也考虑了天然裂缝的存在对水力裂缝延伸的影响^[13-14]。要分析压后复杂缝网对页岩气井产能的影响还需要对储层改造后的气藏开展数值模拟研究，气体在超低孔超低渗透页岩中的流动具有多重机制，包括黏性流、滑脱流动、分子扩散、表面扩散和吸附解吸等^[15]，众多学者对此进行了大量的研究^[16-19]。基于地质工程一体化理念，笔者首先建立目标井区三维地质模型和地质力学模型，然后利用所建立的裂缝扩展模型结合矿场实际泵注历史数据模拟压后复杂缝网几何形态，采用非结构化网格高分辨率数值模拟技术对压后气藏进行数值模拟研究，考虑应力敏感效应优化页岩气分段压裂水平井生产制度及水平段长度。

1 数学模型的建立

1.1 水力裂缝扩展模型

页岩气水平井分段压裂时同一压裂段内通常有4~6簇射孔，随着密切割技术的应用，簇间距离越来越小，多簇裂缝同步扩展时的应力阴影效应不可忽略。首先建立分段多簇压裂水平井物理模型，水平井筒方向与储层水平最小主应力方向平行，此时水力裂缝垂直于井筒方向延伸。数学模型包含3部分。

1.1.1 基岩变形模型

储层岩石在高水压作用下裂缝起裂扩展，属于典型的位移不连续问题。针对此类问题 CROUCH 提出了位移不连续法的求解算法，其基本思想是建立不连续面上应力与位移不连续量之间的关系^[20]。基于此，将裂缝划分为有限单元，裂缝面应力与位移不连续量关系可表示为：

1.1.2 裂缝内压裂液流动

压裂液在裂缝内包括沿裂缝长度和裂缝宽度2 个方向的流动，由于裂缝宽度（毫米级别）相较于裂缝长度（米级别）非常小，因此仅考虑在裂缝长度方向的流动，且该流动模式可看作是泊肃叶流动，其运动方程可表示为：

压裂液为胶液类幂律型流体时，其运动方程为：

假设压裂液不可压缩，裂缝内压裂液流动满足质量守恒方程，可表示为：

假设压裂液滤失满足 Carter 模型，则滤失速度可表示为：

边界条件包括给定裂缝入口端和尖端处截面流量，可表示为：

压裂液在裂缝内的流动问题可由上述控制方程和边界条件求解，但由于边界条件只给定了流量边界条件，解存在不确定性。为保证解的唯一性，引入全局质量守恒方程，其可表示为：

1.1.3 井筒中压裂液流动

水平井多簇压裂时，裂缝间扩展存在应力干扰作用，从井口泵入的压裂液并非均匀地劈分进入各簇裂缝，因此需要对压裂液在井筒中的流动进行模拟。

从流量上，井筒内总注入速率等于每条裂缝注入速率之和，可表示为

从压力上，流体压降包括流体沿水平井筒的摩阻压降和流入射孔的射孔摩阻压降，可表示为：

其中：

1.2 页岩压后流动数学模型

页岩储层改造后，水力裂缝沟通激活了大量的天然裂缝，形成复杂的缝网改造区域，因此页岩气藏压后可分为 2 个区域：储层未改造区域和储层改造区域。储层未改造区域内气体仅在页岩基质孔隙中流动，为单孔隙介质，储层改造区域内同时存在孔隙、天然裂缝和水力裂缝，因此采用双重介质模型进行模拟。引入裂缝系统应力敏感效应，以表征不同生产制度下人工裂缝系统和天然裂缝系统导流能力的变化。

储层未改造区域数学模型考虑页岩气黏性流、 Knudsen 扩散、表面扩散以及气体吸附解吸等复杂渗流机理^[15]，其数学模型为：

其中：

储层改造区域天然裂缝内数学模型为：

水力裂缝内数学模型为：

模型初始化条件：

定义气藏外边界为 Γ₁，生产井内边界为 Γ₂，天然裂缝系统与人工裂缝边界为 Γ₃，SRV 边界为 Γ₄，则边界条件为：

2 模型求解

2.1 裂缝扩展模型求解

裂缝扩展模型的求解包括 2 部分：岩体位移不连续量和流体压力，两场是相互耦合的，采用牛顿迭代法进行求解。选择流体压力 p_i（i为 1~N）、井筒根部流体压力 p₀、裂缝分流量 Q_k（k 为 1~2m）以及时间步长 dt作为迭代变量，通过构建相同数量的方程组来求解。

给定流体压力时，岩体位移不连续量可根据如下方程求得：

采用有限体积方法对压裂液在裂缝内流动方程进行离散，对（5）式进行积分可得：

根据压裂液流动方程，流量可由流压进行表示，即：

（21）式中裂缝宽度可由岩体垂向位移不连续量表达，流速可由流压表达，因此该方程未知量只有流体压力。联合（8），（9），（10）和（20）式可获得与迭代变量相同数量的线性方程，根据牛顿迭代法可求得问题的解。

2.2 流动模型求解

离散裂缝模型可利用标准伽辽金方法推导有限元的积分弱形式，可将方程离散为：

利用上一时间步的压力值求解系数 C _m，C _f 和 D _m，D _f ，继而可求得当前时刻压力值。

3 实例应用

3.1 三维地质和力学模型建立

要开展水力裂缝扩展模拟及压后数值模拟研究首先要建立地质模型。DP2井区是目前南川区块常压页岩气的主要开发区，开发层系为五峰组—龙马溪组①—⑤小层，储层厚度为 31 m 左右，地层压力系数为 1.2，为典型的常压页岩气藏，埋深为 2 500~3 000 m，最大水平主应力为65~68 MPa，最小水平主应力为 59~62 MPa，水平应力差异系数为 0.09。利用三维地震反演结果控制趋势，基于标志层及旋回标定对比，开展水平井小层归位及井周三维构造建模，其中目的层段采用 0.5 m 垂向网格精度，向上、向下网格逐渐加粗（图1a）。水平井测井解释反映储层参数变化受垂向趋势主导，基于水平井等时小层归位，融合测井垂向趋势和地震反演横向趋势控制三维展布，采用趋势建模方法结合变程分析进行高斯模拟从而建立三维属性模型。基于 DP2导眼井FMI成像测井以及地震属性解释的天然裂缝属性，通过分析优选属性，约束建立三维天然裂缝DFN模型（图1b），用于后期分析水力裂缝与天然裂缝的交互关系。利用 DP2 导眼井测井解释的地应力数据以及岩心岩石力学实验资料，首先建立一维地质力学模型，刻画页岩储层纵向应力非均质性，然后结合三维网格及属性体，计算整个区块精细三维地质力学参数，建立储层地质力学模型（图1c）。

3.2 水力压裂缝网模拟

在三维地质力学模型基础上，基于水力裂缝扩展计算模型开展 DP2HF 井压裂缝网模拟研究。由于模型引入了天然裂缝，因此要考虑天然裂缝对水力裂缝转向扩展的影响，其影响因素包括二者之间的夹角（逼近角）、最大最小水平主应力两向应力差值、天然裂缝摩擦系数和内聚力等，国际上有多种裂缝贯穿准则，采用目前国际上主流的 Gu 模型^[13] 进行求解。

首先利用 DP2 井区地质模型和现场实际压裂施工参数作为模型输入，模拟了单簇压裂时裂缝扩展规律。由不同地质参数和压裂液类型对裂缝扩展的影响（图2）可见，岩石杨氏模量越大，或者泊松比越小时，相同注入时间内裂缝半长越大，但在扩展方向上的裂缝宽度越小，容易形成“长窄型”的水力裂缝。当流体流性指数、稠度系数越大，即压裂液表观黏度越大时，压裂液在裂缝内流动时黏性能量耗散越大，相同压裂时间获得的裂缝半长有所下降，同时裂缝内流体压力越大，对裂缝壁面作用力也越大，进而促使裂缝宽度增大，后期利于支撑剂进入裂缝深处，提升支撑效果，改善裂缝导流能力。

由水平井多簇压裂时缝网扩展模拟结果（图3） 可知，水力裂缝延伸方向主要沿最大水平主应力方向进行扩展（模型中水平井筒方向垂直于最大水平主应力方向）。多簇压裂时，簇间裂缝的应力阴影效应导致中间裂缝宽度小于两边裂缝宽度，且中间裂缝长度略小于两边裂缝长度。同时由于天然裂缝的存在，未能贯穿天然裂缝的水力裂缝发生转向，沿着天然裂缝方向延伸，增加了缝网的复杂程度。当天然裂缝分布密度越大、与水力裂缝夹角越小、水平地应力差越小时，缝网改造复杂程度越高。同时由于压裂液向天然裂缝内的滤失作用，部分水力裂缝贯穿天然裂缝后继续扩展延伸能力下降，导致水力裂缝长度受限。

利用 DP2HF 井现场实际分段分簇设计和压裂施工参数开展压裂缝网模拟研究，该井共分27段进行压裂，由裂缝长度和高度模拟结果（图4）可以看出，不同段间水力裂缝的长度和高度有较大的差异，模拟平均改造裂缝长度为206 m，支撑裂缝长度为189 m，平均改造裂缝高度为31 m，有效支撑裂缝高度为4 m。水力裂缝在横向延伸方向上支撑效果较好，有效支撑裂缝长度与水力裂缝长度匹配度较高，在纵向上有效支撑裂缝高度远低于水力裂缝高度，还有较大的提高空间。

3.3 压后数值模拟

页岩储层压后存在基质孔隙、天然裂缝和水力裂缝 3 种流动空间，在不同尺度空间内页岩气的运移机制不同，利用建立的考虑复杂渗流机理和应力敏感效应的页岩气藏数值模型开展 DP2HF 井压后数值模拟研究。

3.3.1 生产制度优化

页岩气从压裂水平井中产出涉及到非常复杂的渗流机理和作用机制，而这些机理和机制又互相影响互相制约，很难准确区分各自对页岩气井产量的贡献，这也导致页岩气井合理生产制度的制定变得更加困难。例如初期放大压差生产增加压裂液返排有助于加快气体从基岩到裂缝的运移速度、促进吸附气的解吸释放，从而增加初期产量；但是大压差生产同时又会因为应力敏感效应引起裂缝闭合，导致裂缝连通性及导流能力下降，进而导致产量快速下降。地质工程一体化技术将储层、地应力、水力压裂缝网模拟、复杂渗流机理等耦合在一起，可以研究压裂缝网扩展过程中、试气返排以及衰竭生产过程中，储层应力、渗透率和缝网导流能力等参数的变化，综合考虑各种机理和机制对生产及单井最终可采储量的影响，从而优化合理的生产制度。以 DP2井区地质模型为基础，建立动态全耦合的数值模拟模型，新模型保留了原始模型各种储层物性参数、岩石力学参数、天然裂缝分布以及初始三维应力场。针对南川区块页岩气井试气产气量变化范围，设定7种不同的生产制度，初始日产气量从 4×10⁴ m³ /d 到 30×10⁴ m³ /d，研究不同生产制度下的单井稳产期及累积产气量。由 DP2HF 井不同生产制度下产气量预测曲线（图5）可知，页岩气井初期产气量越低，可以获得更长的稳产时间及累积产气量。初期配产越高，可以在早期快速获得较高的累积产气量，随着生产时间增加，低配产工作制度下的累积产气量逐渐超过高配产制度下的累积产气量，而最终的预测可采储量也更大（预测15 a）。

当然，按照低配产生产制度虽然可以获得更高的最终可采储量，但未必是最具经济效益的生产制度。按照研究区 2 000 m 水平段长单井钻采投资为 5 500万元、气价为1 420元/10³ m³、页岩气操作成本为 283元/10³ m³、增值税税率为 9%来计算不同生产制度下的净现值。由不同生产制度下单井稳产时间、15 a累积产气量和净现值（表1）可以看出，初期配产 4×10⁴ m³ /d 虽然最终可采储量最高，但是由于早期累积产气量过低，导致最终净现值也最低。经济效益最好的生产制度为初期配产 6×10⁴ m³ /d，超过该产气量后，随着初期配产的增加，经济效益逐渐下降。

3.3.2 水平段长优化

分别设置了 7 口不同水平段长的水平井，水平段长从1 250 m到2 750 m。采用研究区目前常用的压裂砂量、压裂液量等施工参数开展压裂缝网模拟研究，在此基础上开展压后数值模拟研究，预测生产15 a后的累积产气量。不同水平段长钻采投资不同，单井合理配产不同。按照 2 000 m 水平段长单井投资5 500万元、合理配产6×10⁴ m³ /d折算其他水平段长下的投资和合理配产，结合 15 a累积产气量计算不同水平段长下的内部收益率。由表2可以看出，随着水平段长的增加，内部收益率先增加后下降。当水平段长达到 1 500 m 时，内部收益率能够达到8%，水平段长增加到2 000 m时，内部收益率达到最高点，随后逐渐下降。当水平段长超过2 500 m 后，内部收益率降至8%以下，因此推荐合理水平段长为2 000 m。

4 结论

基于地质工程一体化的理念，建立了水力裂缝扩展模型和考虑复杂渗流机理的页岩气藏压后数值模型，该模型考虑了压裂时水力裂缝间应力阴影效应以及页岩气藏复杂的渗流机理和作用机制，模拟结果更加接近现场实际。

通过单簇压裂模拟明确了储层参数及流体性质对水力裂缝扩展的影响，DP2HF井压裂模拟结果显示水力压裂可形成复杂缝网结构，水力裂缝在横向延伸方向上有效支撑较好，垂向上有效支撑裂缝高度较低，有进一步改进提产的空间。

压后数值模拟显示，由于应力敏感效应单井初期配产高低对最终可采储量影响很大，综合经济评价 DP2HF 井合理生产制度为 6×10⁴ m³ /d，合理水平段长为2 000 m。

符号解释

C_f ——裂缝单元刚度矩阵；

C_l ——流体滤失系数，m/s^1/2；

C_m——基岩单元刚度矩阵；

C_nn，C_ns，C_sn，C_ss——弹性系数矩阵分量；

d_p，i ——裂缝i射孔直径，m；

D_f ——裂缝单元压缩系数矩阵；

D_m——基岩单元压缩系数矩阵；

D_n，D_s——裂缝面法向和切向位移不连续量，m；

G——三维修正系数^[21]；

i，j——裂缝单元编号；

k_f，app——天然裂缝系统视渗透率，m²；

k_F-stress——考虑应力敏感性的水力裂缝系统渗透率，m²；

k_m，app——页岩基质视渗透率，m² ；

K——稠度系数，Pa·sⁿ；

l_i，l_j ——裂缝单元i和j的长度，m；

L——裂缝长度，m；

L_i ——裂缝i的长度，m；

下标 m，f，F——基质孔隙系统、天然裂缝系统、水力裂缝系统；

m——裂缝条数；

M——甲烷气体摩尔质量，kg/mol；

n——流性指数；

n_p，i ——裂缝i射孔数；

N——裂缝单元个数；

p——流动系统孔隙压力，Pa；

p₀——井筒根部流体压力，Pa；

p_cf，i ——裂缝i井筒摩阻压降，Pa；

p_i ——气藏初始压力，Pa；

p_i，p_j——裂缝单元i和j的流体压力，Pa；

p_L——Langmuir压力，Pa；

pⁱ _net ——流体净压力，Pa；

p_pf，i ——裂缝i射孔摩阻压降，Pa；

p_w——井底流压，Pa；

p_w，i——裂缝i入口处流体压力，Pa；

q——截面流量，m² /s；

q_i——裂缝i入口端截面流量，m² /s；

q_ij——离散后裂缝单元（i，j）的流量，m² /s；

q_l ——流体滤失速度，m/s；

Q_i ——裂缝i的注入速率，m² /s；

Q_k——裂缝k的分流量，m² /s；

Q_p——基岩-裂缝窜流量，kg/s；

Q_T——总注入速率，m² /s；

R——通用气体常数，J/（K·mol），取值为8.314；

s——沿裂缝扩展方向的位移，m；

t——注入时间，s；

T——气藏温度，K；

V——体积，m³；

V_L——Langmuir体积，m³ /kg；

V_std——标准状况下摩尔体积，m³ /mol；

w——裂缝宽度，m；

x，y，z——三维空间点位置，m；

Z——气体压缩因子；

a——压后裂缝宽度，m；

α^* ——形状因子，1/m²；

μ——气体黏度，Pa·s；

μ_l ——压裂液黏度，Pa·s；

ρ——气体密度，kg/m³；

ρ_s ——页岩基岩密度，kg/m³；

ρ_sl——泥浆密度，kg/m³；

ρ_w——水密度，kg/m³；

σ_n，σ_s——裂缝面法向和切向应力，Pa；

τ₀——流体滤失开始时间，s；

ϕ——孔隙度，%。

参考文献