CN112031756B - 一种页岩气藏压裂井组生产动态数值模拟方法 - Google Patents

Abstract

本发明公开了一种页岩气藏压裂井组生产动态数值模拟方法，包括以下步骤：页岩储层和压裂水平井三维地质体生成及非结构四面体网格离散；建立页岩气藏多级压裂水平井综合渗流数学模型；基于控制体积有限元法构建综合渗流数学模型的全隐式数值模型；建立基于四面体网格离散的全隐式多级压裂水平井数值井模型，并嵌入全隐式数值模型，构建完整的多系统同步联立求解矩阵；对多系统同步联立求解矩阵迭代求解，并分析朗格缪尔体积等对页岩气藏压裂井组生产动态的影响。本发明弥补了现有页岩气藏二维模拟器的不足，综合考虑页岩储层物性层间非均质性和水力压裂裂缝纵向延伸扩展特征，实现了三维深层页岩储层中立体开发压裂井组生产动态的准确模拟和预测。

Description

一种页岩气藏压裂井组生产动态数值模拟方法

技术领域

本发明涉及页岩气藏的勘探与开发技术领域，具体涉一种页岩气藏压裂井组生产动态数值模拟方法。

背景技术

目前，我国已初步形成本土化的页岩气开发技术，可支撑3500m以浅的页岩气资源开发利用。但是，对于埋深超过3500m的深层页岩的规模效益开发尚面临诸多待解的难题。页岩气藏开发实践表明：长井段水平井体积压裂技术是获得页岩气藏工业产能的关键。然而，由于深层页岩层间物性差异大，仍然套用浅层页岩针对主力产层的平面压裂水平井井组开采方式，将导致纵向动用程度低，单井最终可采储量较低。同时，深层页岩孔隙结构特征更加复杂，页岩基质有机质和无机矿物表面气体附存状态差异大，采用传统的单重/双重连续介质模型难以准确表征深层页岩多尺度、强非线性流动规律。近年来，基于正交网格和有限差分方法的嵌入式离散裂缝模拟器、基于非结构三角形和有限元方法的离散裂缝模拟器被广泛应用于页岩气藏压裂水平井生产动态的模拟预测。但是，上述数值模拟方法均普遍采用二维或拟三维网格对页岩储层及压裂水平井进行空间网格离散，未考虑纵向流动特征和压裂裂缝纵向延伸扩展特征，难以实现对深层页岩立体开发压裂井组生产动态的模拟预测。

发明内容

鉴于上述原因，本发明的目的是提供一种页岩气藏压裂井组生产动态数值模拟方法，本方法采用四面体网格对深层页岩立体开发压裂井组三维地质体和三维压裂裂缝进行空间离散；综合考虑页岩基质滑脱、克努森扩散、吸附解吸、天然裂缝应力敏感和压裂裂缝线性流动特征，基于多重连续介质-离散裂缝耦合模型构建页岩气藏多尺度流动综合渗流模型；利用控制体积有限元法建立全隐式数值模型，并获得数值解；应用数值模拟方法对朗格缪尔体积、压裂裂缝渗透率和压裂井开采模式对压裂井组生产动态的影响进行了敏感性分析。本发明的技术方案为：

一种页岩气藏压裂井组生产动态数值模拟方法，包括以下步骤：

S1：页岩储层和压裂水平井三维地质体生成及非结构四面体网格离散；

S2：基于页岩气藏孔隙结构特征和多重连续介质-离散裂缝耦合模型，结合初始条件和边界条件，建立页岩气藏多级压裂水平井综合渗流数学模型；

S3：基于控制体积有限元法，构建综合渗流数学模型的全隐式数值模型；

S4：对内边界条件进行处理，建立基于四面体网格离散的全隐式多级压裂水平井数值井模型，并嵌入步骤S3的全隐式数值模型，构建完整的多系统同步联立求解矩阵；

S5：对构建的多系统同步联立求解矩阵进行迭代求解，并分析朗格缪尔体积、压裂裂缝渗透率和压裂井开采模式对页岩气藏压裂井组生产动态的影响。

进一步，所述步骤S1还包括以下步骤：

S11：根据研究工区的实际地质情况和水平井井轨迹及压裂裂缝分布特征，生成地质体；

S12：按照开源四面体网格剖分软件Tetgen的数据结构要求，对地质体数据进行编辑导入，生成网格离散模型。

进一步，所述步骤S2还包括如下步骤：

S21：基于页岩气藏孔隙结构特征和多重连续介质-离散裂缝耦合模型，将页岩多尺度储集空间划分为：微纳尺度有机孔、介观无机宏孔、天然裂缝和大尺度水力压裂裂缝系统；构建耦合微纳孔隙滑脱、克努森扩散、有机质表面吸附解吸、天然缝网应力敏感效应以及各储集空间流量传输和水力压裂裂缝线性流动的综合渗流方程为：

其中，

式中符号定义：O、I、f、F分别代表有机质、无机质、天然裂缝和水力压裂裂缝系统；k_a为考虑滑脱和克努森扩散的有机质微纳孔隙表观渗透率，mD；k_m为页岩基质气测渗透率，mD；k_fe为考虑应力敏感影响的页岩天然裂缝渗透率，mD；k_fi为页岩天然裂缝原始渗透率，mD；d_f为天然裂缝应力敏感系数，MPa^-1；β为滑脱系数；k_F为水力压裂裂缝平均渗透率，mD；C_gO为有机质孔隙中气体压缩系数，MPa^-1；D_k为克努森扩散系数，m²/s；μ_gO为有机质微纳孔隙中气体粘度，mPa·s；μ_gl为无机质介观孔隙中气体粘度，mPa·s；μ_gf为天然裂缝系统中气体粘度，mPa·s；μ_gF为压裂裂缝系统中气体粘度，mPa·s；B_gO为有机质微纳孔隙中气体体积系数；B_gl为无机质介观孔隙中气体体积系数；B_gf为天然裂缝系统中气体体积系数；B_gF为压裂裂缝系统中气体体积系数；q_des为标况下的吸附解吸量，方/天；V_L为朗格缪尔体积，cm³/g；p_L为朗格缪尔压力，MPa；ρ_s为页岩岩石密度，g/cm³；q_g,lO为有机质和无机质系统间气体窜流流量，方/天；q_g,fl为天然裂缝系统和无机质系统间气体窜流流量，方/天；q_g,Ff为压裂裂缝系统和天然裂缝系统间气体窜流流量，方/天；α_(·)为窜流形状因子，m^-2；φ_m为页岩基质孔隙度；φ_f为天然裂缝系统孔隙度；φ_F为压裂裂缝系统孔隙度；f_O为有机质在页岩基质中所占体积比；f_l为无机质在页岩基质中所占体积比；p_gO为有机质孔隙内气体压力，MPa；p_gl为无机质孔隙内气体压力，MPa；p_gf为天然裂缝内气体压力，MPa；p_gF为水力压裂裂缝内气体压力，MPa；p_gf0为天然裂缝内原始气体压力，MPa；q_gsct为标况下的地面产气量，方/天。

S22：设置初始条件和外边界条件

初始条件：

p_gf＝p_gf0 (5)

p_gl＝p_gO＝p_gm0 (6)

p_gF＝p_gF0 (7)

式中符号定义：p_gm0为页岩基质原始压力，MPa；p_gF0为水力压裂裂缝系统原始压力，MPa。

外边界条件：

进一步，所述步骤S3还包括如下步骤：

S31：选取一个包含压裂裂缝和有机质-无机质-天然裂缝连续体的四面体单元，构建单元特性矩阵，定义单元内不同系统压力试探解为：

式中符号定义：

为l系统中单元网格的压力试探解；N_lv为l系统中单元网格的插值函数；v＝1，2，3，4代表四面体的四个顶点，P_glv为l系统中单元网格的各顶点压力值。

基于Galerkin加权余量法和高斯定律，将式(9)带入到步骤S2中建立的综合渗流数学模型中，构建积分弱形式，则对流项变化为：

其中，渗透率张量计算公式为：

式中符号定义：V_1l为l系统中四面体单元网格的体积，m³；k_lx表示l系统中x方向渗透率，mD；k_ly表示l系统中y方向渗透率，mD；k_lz表示l系统中z方向渗透率，mD。

将式(10)进一步展开得到流入四面体网格任意顶点i的流量计算式为：

其中

式中符号定义：b_i和b_v为四面体网格顶点在x方向坐标系数，m；c_i和c_v为四面体网格顶点在y方向坐标系数，m；d_i和d_v为四面体网格顶点在z方向坐标系数，m；p_il为l系统中四面体网格顶点i处的压力值，MPa；p_vl为l系统中四面体网格顶点v处的压力值，MPa。

对四面体网格单元的四个顶点构建如式(12)的流量计算式，从而得到四面体网格单元传导率矩阵格式为：

进一步利用控制体积有限元法，建立不同流量传输特性矩阵，包括，

窜流项单元特性矩阵为：

吸附解吸项单元特性矩阵为：

时间项单元特性矩阵为：

同时，基于离散裂缝模型原理，将三维水力压裂裂缝“降维”处理为宽度为w_F的二维平面，采用非结构三角形网格对裂缝面进行空间离散，并通过式(17)将水力压裂缝线性流动方程式(4)嵌入到多重连续介质渗流空间中；基于“窜流平衡”理论，在天然裂缝与压裂裂缝的相交面上，满足天然裂缝网格压力与压裂裂缝网格压力相等，从而约去压裂裂缝网格压力未知变量，并抵消掉天然裂缝与压裂裂缝之间的窜流量q_g,Ff：

式中符号定义：PDEs代表偏微分方程组；w_F为压裂裂缝宽度，m；V_t为天然裂缝和压裂裂缝网格单元总体积，m³；V_f为天然裂缝网格单元体积，m³；Ω_F为二维压裂裂缝单元面积，m²；

采用全隐式计算格式，则下一时步的各系统压力值可表示为：

式中符号定义：δ为算子，表示第k次迭代到k+1次迭代后的压力变化；n为上一时步；n+1为下一时步。

将式(13)～式(18)，代入步骤S21中综合渗流方程，建立多系统全隐式单元特性矩阵为：

对离散区域内的每一个四面体网格，建立类似的单元矩阵；假设存在N个网格节点，通过整体叠加，形成关于有机质、无机质和裂缝系统的3N×3N大矩阵，即综合渗流数学模型的全隐式数值模型：

[k]_3N*3N[δX]_3N*1＝[R]_3N*1 (22)

式中符号定义：K为系数矩阵；δX为未知变量矩阵；R为余量矩阵。

进一步，所述步骤S4还包括如下步骤：

S41：对于式(4)中的标况下的压裂井组地面产气量q_gsct，采用修正的拟稳态计算公式展开得到：

其中：

式中符号定义：N_well为压裂水平井井数；N_Fi为i压裂井的水力压裂裂缝总条数；p_bhi为i压裂井井底流压，MPa；p_ave,ij为i压裂井第j条压裂裂缝井源汇点所在网格块的平均压力，MPa；Pl_g,ij为i压裂井第j条压裂裂缝采气指数，方/天/MPa；r_wi为i压裂井的井筒半径，m；r_0，ij为i压裂井第j条压裂裂缝井源汇点的等效井半径，m；

S42：对于压裂裂缝井等效井半径r_0，ij，将气体向压裂井的流动视为沿裂缝面向厚度为裂缝宽度w_F，ij的等效直井的拟稳态径向流动，则各裂缝面等效井半径可表示为：

式中符号定义：A_ij为i压裂井第j条压裂裂缝井源汇点二维裂缝面的面积，m²。

S43：将式(18)带入式(23)中，得到全隐式压裂水平井的数值井计算格式为：

S44：将式(25)带入到式(22)中，构建完整的页岩储层和压裂井组多系统同步联立求解矩阵。

进一步，所述步骤S5还包括如下步骤：

S51：由有限元基本原理，对于封闭外边界条件，添加0向量到式(22)系数矩阵K中；

S52：采用共轭梯度法对完整的页岩储层和压裂井组多系统同步联立矩阵进行求解，通过牛顿-拉夫逊迭代方法，获得有机质、无机质以及裂缝系统中各离散网格节点在一个时步下的压力变化值：δp_gf1,…δp_gfN；δp_gl1,…δp_gl1N；δp_gO1,…δp_gO1N；进而由式(18)获得下一时步的值；

S53：输出各时步下的离散网格节点压力值，以及定井底流压生产条件下的压裂井组日产气量和累产气量。

本发明有益效果如下：1.基于多重连续介质-离散裂缝耦合模型，建立了页岩气藏压裂井组综合渗流模型，实现了对微纳尺度有机孔、介观无机宏孔、天然裂缝和大尺度水力压裂裂缝系统多尺度孔隙空间中耦合微纳孔隙滑脱、克努森扩散、有机质表面吸附解吸、天然缝网应力敏感效应以及各系统间流量传输和水力压裂裂缝线性流的强非线性流动规律的联合表征；2.弥补了现有页岩气藏二维模拟器的不足，综合考虑页岩储层物性层间非均质性和水力压裂裂缝纵向延伸扩展特征，实现了三维深层页岩储层中立体开发压裂井组生产动态的准确模拟和预测。

附图说明

图1是页岩气藏压裂井组生产动态数值模拟流程图；

图2是页岩气藏压裂水平井井组三维地质体及网格离散示意图；

图3是四面体网格单元及多重连续介质微元体示意图；

图4是压裂水平井拟稳态井模型示意图；

图5是朗格缪尔体积对生产动态的影响图；

图6是压裂裂缝渗透率对生产动态的影响图；

图7是压裂井开采模式对生产动态的影响图。

具体实施方式

下面将结合本发明实施例，对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

一种页岩气藏压裂井组生产动态数值模拟方法，如图1所示，包括以下步骤：

S1：页岩储层和压裂水平井三维地质体生成及非结构四面体网格离散；

S2：基于页岩气藏孔隙结构特征和多重连续介质-离散裂缝耦合模型，结合初始条件和边界条件，建立页岩气藏多级压裂水平井综合渗流数学模型；

S3：基于控制体积有限元法，构建综合渗流数学模型的全隐式数值模型；

S4：对内边界条件进行处理，建立基于四面体网格离散的全隐式多级压裂水平井数值井模型，并嵌入步骤S3的全隐式数值模型，构建完整的多系统同步联立求解矩阵；

S5：对构建的多系统同步联立求解矩阵进行迭代求解，并分析朗格缪尔体积、压裂裂缝渗透率和压裂井开采模式对页岩气藏压裂井组生产动态的影响。

所述步骤S1还包括以下步骤：

S11：根据研究工区的实际地质情况和水平井井轨迹及压裂裂缝分布特征，生成地质体；

S12：按照开源四面体网格剖分软件Tetgen的数据结构要求，对地质体数据进行编辑导入，生成网格离散模型，如图2。

所述步骤S2还包括如下步骤：

S21：基于页岩气藏孔隙结构特征和多重连续介质-离散裂缝耦合模型，将页岩多尺度储集空间划分为：微纳尺度有机孔、介观无机宏孔、天然裂缝和大尺度水力压裂裂缝系统；构建耦合微纳孔隙滑脱、克努森扩散、有机质表面吸附解吸、天然缝网应力敏感效应以及各储集空间流量传输和水力压裂裂缝线性流的综合渗流方程为：

其中，

式中符号定义：O、I、f、F分别代表有机质、无机质、天然裂缝和水力压裂裂缝系统；k_a为考虑滑脱和克努森扩散的有机质微纳孔隙表观渗透率，mD；k_m为页岩基质气测渗透率，mD；k_fe为考虑应力敏感影响的页岩天然裂缝渗透率，mD；k_fi为页岩天然裂缝原始渗透率，mD；d_f为天然裂缝应力敏感系数，MPa^-1；β为滑脱系数；k_F为水力压裂裂缝平均渗透率，mD；C_gO为有机质孔隙中气体压缩系数，MPa^-1；D_k为克努森扩散系数，m²/s；μ_gO为有机质微纳孔隙中气体粘度，mPa·s；μ_gl为无机质介观孔隙中气体粘度，mPa·s；μ_gf为天然裂缝系统中气体粘度，mPa·s；μ_gF为压裂裂缝系统中气体粘度，mPa·s；B_gO为有机质微纳孔隙中气体体积系数；B_gl为无机质介观孔隙中气体体积系数；B_gf为天然裂缝系统中气体体积系数；B_gF为压裂裂缝系统中气体体积系数；q_des为标况下的吸附解吸量，方/天；V_L为朗格缪尔体积，cm³/g；p_L为朗格缪尔压力，MPa；ρ_s为页岩岩石密度，g/cm³；q_g,lO为有机质和无机质系统间气体窜流流量，方/天；q_g,fl为天然裂缝系统和无机质系统间气体窜流流量，方/天；q_g,Ff为压裂裂缝系统和天然裂缝系统间气体窜流流量，方/天；α_(·)为窜流形状因子，m^-2；φ_m为页岩基质孔隙度；φ_f为天然裂缝系统孔隙度；φ_F为压裂裂缝系统孔隙度；f_O为有机质在页岩基质中所占体积比；f_l为无机质在页岩基质中所占体积比；P_gO为有机质孔隙内气体压力，MPa；P_gl为无机质孔隙内气体压力，MPa；P_gf为天然裂缝内气体压力，MPa；P_gF为水力压裂裂缝内气体压力，MPa；P_gf0为天然裂缝内原始气体压力，MPa；q_gsct为标况下的地面产气量，方/天。

S22：设置初始条件和外边界条件

初始条件：

p_gf＝p_gf0 (5)

p_gl＝p_gO＝p_gm0 (6)

p_gF＝p_gF0 (7)

式中符号定义：P_gm0为页岩基质原始压力，MPa；P_gF0为水力压裂裂缝系统原始压力，MPa；

外边界条件：

进一步，所述步骤S3还包括如下步骤：

S31：选取一个包含压裂裂缝和有机质-无机质-天然裂缝连续体的四面体单元，如图3，构建单元特性矩阵，定义单元内不同系统压力试探解为：

式中符号定义：

为l系统中单元网格的压力试探解；N_lv为l系统中单元网格的插值函数，v＝1，2，3，4代表四面体的四个顶点，P_glv为l系统中单元网格的各顶点压力值。

基于Galerkin加权余量法和高斯定律，将压力试探解带入到步骤S2建立的综合渗流数学模型中，构建积分弱形式，则对流项变化为：

其中，渗透率张量计算公式为：

式中符号定义：V_1l为l系统中四面体单元网格的体积，m³；k_lx表示l系统中x方向渗透率，mD；k_ly表示l系统中y方向渗透率，mD；k_lz表示l系统中z方向渗透率，mD；

将式(10)进一步展开得到流入四面体网格任意顶点i的流量计算式为：

其中

式中符号定义：b_i和b_v为四面体网格顶点在x方向坐标系数，m；c_i和c_v为四面体网格顶点在y方向坐标系数，m；d_i和d_v为四面体网格顶点在z方向坐标系数，m；p_il为l系统中四面体网格顶点i处的压力值，MPa；p_vl为l系统中四面体网格顶点v处的压力值，MPa。

对四面体网格单元的四个顶点构建如式(12)的流量计算式，从而得到四面体网格单元传导率矩阵格式为：

进一步利用控制体积有限元法，建立不同流量传输特性矩阵，包括，

窜流项单元特性矩阵为：

吸附解吸项单元特性矩阵为：

时间项单元特性矩阵为：

同时，基于离散裂缝模型原理，将三维水力压裂裂缝“降维”处理为宽度为w_F的二维平面，采用非结构三角形网格对裂缝面进行空间离散，并通过式(17)将水力压裂缝线性流动方程式(4)嵌入到多重连续介质渗流空间中；基于“窜流平衡”理论，在天然裂缝与压裂裂缝的相交面上，满足天然裂缝网格压力与压裂裂缝网格压力相等，从而约去压裂裂缝网格压力未知变量，并抵消掉天然裂缝与压裂裂缝之间的窜流量q_g,Ff：

式中符号定义：PDEs代表偏微分方程组；w_F为压裂裂缝宽度，m；V_t为天然裂缝和压裂裂缝网格单元总体积，m³；V_f为天然裂缝网格单元体积，m³；Ω_F为二维压裂裂缝单元面积，m²；

采用全隐式计算格式，则下一时步的各系统压力值可表示为：

式中符号定义：δ为算子，表示第k次迭代到k+1次迭代后的压力变化；n为上一时步；n+1为下一时步。

将式(13)～式(18)，代入步骤S21中综合渗流方程，建立多系统全隐式单元特性矩阵为：

对离散区域内的每一个四面体网格，建立类似的单元矩阵；假设存在N个网格节点，通过整体叠加，形成关于有机质、无机质和裂缝系统的3N×3N大矩阵，即综合渗流数学模型的全隐式数值模型：

[k]_3N*3N[δX]_3N*1＝[R]_3N*1 (22)

式中符号定义：K为系数矩阵；δX为未知变量矩阵；R为余量矩阵。

进一步，所述步骤S4还包括如下步骤：

S41：对于式(4)中的标况下的压裂井组地面产气量q_gsct，采用修正的拟稳态计算公式展开得到：

其中：

式中符号定义：N_well为压裂水平井井数；N_Fi为i压裂井的水力压裂裂缝总条数；p_bhi为i压裂井井底流压，MPa；p_ave,ij为i压裂井第j条压裂裂缝井源汇点所在网格块的平均压力，MPa；PI_g,ij为i压裂井第j条压裂裂缝采气指数，方/天/MPa；r_wi为i压裂井的井筒半径，m；r_0，ij为i压裂井第j条压裂裂缝井源汇点的等效井半径，m；

S42：对于压裂裂缝井等效井半径r_0,ij，如图4，将气体向压裂井的流动视为沿裂缝面向厚度为裂缝宽度w_F,ij的等效直井的拟稳态径向流动，则各裂缝面等效井半径可表示为：

式中符号定义：A_ij为i压裂井第j条压裂裂缝井源汇点二维裂缝面的面积，m²。

S43：将式(18)带入式(23)中，得到全隐式压裂水平井的数值井计算格式为：

S44：将式(25)带入到式(22)中，构建完整的页岩储层和压裂井组多系统同步联立求解矩阵。

进一步，所述步骤S5还包括如下步骤：

S51：由有限元基本原理，对于封闭外边界条件，添加0向量到式(22)系数矩阵K中；

S52：采用共轭梯度法对完整的页岩储层和压裂井组多系统同步联立矩阵进行求解，通过牛顿-拉夫逊迭代方法，获得有机质、无机质以及裂缝系统中各离散网格节点在一个时步下的压力变化值：δp_gf1,…δp_gfN；δp_gl1,…δp_gl1N；δp_gO1,…δp_gO1N；进而由式(18)获得下一时步的值；

S53：输出各时步下的离散网格节点压力值，以及定井底流压生产条件下的压裂井组日产气量和累产气量。

图5是朗格缪尔体积对生产动态的影响图；图6是压裂裂缝渗透率对生产动态的影响图；图7是压裂井开采模式对生产动态的影响图。分析图5-7可知：①朗格缪尔体积越大，页岩基质吸附能力越强，随着压力降低解吸出的吸附气量越多，压裂井组日产气量和累产气量均增加；②压裂裂缝渗透率越大，气体在压裂裂缝中流动阻力越小，气井早期日产气量越大，累产气量越大；③同单压裂井开采模式进行对比分析，对于深层页岩气藏，采用立体压裂井组模式开采能够更有效地提高页岩纵向动用程度，有利于提高累产气量。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (3)

1.一种页岩气藏压裂井组生产动态数值模拟方法，其特征在于，包括以下步骤：

S1：页岩储层和压裂水平井三维地质体生成及非结构四面体网格离散；

S2：基于页岩气藏孔隙结构特征和多重连续介质-离散裂缝耦合模型，结合初始条件和边界条件，建立页岩气藏多级压裂水平井综合渗流数学模型；

S21：基于页岩气藏孔隙结构特征和多重连续介质-离散裂缝耦合模型，将页岩多尺度储集空间划分为：微纳尺度有机孔、介观无机宏孔、天然裂缝和大尺度水力压裂裂缝系统；构建耦合微纳孔隙滑脱、克努森扩散、有机质表面吸附解吸、天然缝网应力敏感效应以及各储集空间流量传输和水力压裂裂缝线性流动的综合渗流方程为：

其中，

式中符号定义：O、I、f、F分别代表有机质、无机质、天然裂缝和水力压裂裂缝系统；k_a为考虑滑脱和克努森扩散的有机质微纳孔隙表观渗透率，mD；k_m为页岩基质气测渗透率，mD；k_fe为考虑应力敏感影响的页岩天然裂缝渗透率，mD；k_fi为页岩天然裂缝原始渗透率，mD；d_f为天然裂缝应力敏感系数，MPa^-1；β为滑脱系数；k_F为水力压裂裂缝平均渗透率，mD；C_gO为有机质孔隙中气体压缩系数，MPa^-1；D_k为克努森扩散系数，m²/s；μ_gO为有机质微纳孔隙中气体粘度，mPa·s；μ_gI为无机质介观孔隙中气体粘度，mPa·s；μ_gf为天然裂缝系统中气体粘度，mPa·s；μ_gF为压裂裂缝系统中气体粘度，mPa·s；B_gO为有机质微纳孔隙中气体体积系数；B_gI为无机质介观孔隙中气体体积系数；B_gf为天然裂缝系统中气体体积系数；B_gF为压裂裂缝系统中气体体积系数；q_des为标况下的吸附解吸量，1/d；V_L为朗格缪尔体积，cm³/g；p_L为朗格缪尔压力，MPa；ρ_s为页岩岩石密度，g/cm³；q_g,IO为有机质和无机质系统间气体窜流流量，1/d；q_g,fI为天然裂缝系统和无机质系统间气体窜流流量，1/d；q_g,Ff为压裂裂缝系统和天然裂缝系统间气体窜流流量，1/d；α_(·)为窜流形状因子，m^-2；φ_m为页岩基质孔隙度；φ_f为天然裂缝系统孔隙度；φ_F为压裂裂缝系统孔隙度；f_O为有机质在页岩基质中所占体积比；f_I为无机质在页岩基质中所占体积比；P_gO为有机质孔隙内气体压力，MPa；P_gI为无机质孔隙内气体压力，MPa；P_gf为天然裂缝内气体压力，MPa；P_gF为水力压裂裂缝内气体压力，MPa；P_gf0为天然裂缝内原始气体压力，MPa；q_gsct为标况下的地面产气量，1/d；

S22：设置初始条件和外边界条件

初始条件：

p_gf＝p_gf0 (5)

p_gI＝p_gO＝p_gm0 (6)

p_gF＝p_gF0 (7)

式中符号定义：p_gm0为页岩基质原始压力，MPa；p_gF0为水力压裂裂缝系统原始压力，MPa；外边界条件：

S3：基于控制体积有限元法，构建综合渗流数学模型的全隐式数值模型；

S31：选取一个包含压裂裂缝和有机质-无机质-天然裂缝连续体的四面体单元，构建单元特性矩阵，定义单元内不同系统压力试探解为：

式中符号定义：

为l系统中单元网格的压力试探解；N_lv为l系统中单元网格的插值函数，v＝1，2，3，4代表四面体的四个顶点；P_glv为l系统中单元网格的各顶点压力值；

基于Galerkin加权余量法和高斯定律，将式(9)带入到步骤S2建立的综合渗流数学模型中，构建积分弱形式，则对流项变化为：

其中，渗透率张量计算公式为：

式中符号定义：V_1l为l系统中四面体单元网格的体积，m³；k_lx表示l系统中x方向渗透率，mD；k_ly表示l系统中y方向渗透率，mD；k_lz表示l系统中z方向渗透率，mD；

将式(10)进一步展开得到流入四面体网格任意顶点i的流量计算式为：

其中

式中符号定义：b_i和b_v为四面体网格顶点在x方向坐标系数，m；c_i和c_v为四面体网格顶点在y方向坐标系数，m；d_i和d_v为四面体网格顶点在z方向坐标系数，m；p_il为l系统中四面体网格顶点i处的压力值，MPa；p_vl为l系统中四面体网格顶点v处的压力值，MPa；

对四面体网格单元的四个顶点构建类似的流量计算式，从而得到四面体网格单元传导率矩阵格式为：

利用控制体积有限元法，建立不同流量传输特性矩阵，包括，

窜流项单元特性矩阵为：

其中，

吸附解吸项单元特性矩阵为：

时间项单元特性矩阵为：

同时，基于离散裂缝模型原理，将三维水力压裂裂缝“降维”处理为宽度为w_F的二维平面，采用非结构三角形网格对裂缝面进行空间离散，并通过式(17)将水力压裂缝线性流动方程式(4)嵌入到多重连续介质渗流空间中；基于“窜流平衡”理论，在天然裂缝与压裂裂缝的相交面上，满足天然裂缝网格压力与压裂裂缝网格压力相等，从而约去压裂裂缝网格压力未知变量，并抵消掉天然裂缝与压裂裂缝之间的窜流量q_g,Ff：

式中符号定义：PDEs代表偏微分方程组；w_F为压裂裂缝宽度，m；V_t为天然裂缝和压裂裂缝网格单元总体积，m³；V_f为天然裂缝网格单元体积，m³；Ω_F为二维压裂裂缝单元面积，m²；

采用全隐式计算格式，则下一时步的各系统压力值表示为：

式中符号定义：δ为算子，表示第k次迭代到k+1次迭代后的压力变化；n为上一时步；n+1为下一时步；

将式(13)～式(18)，代入步骤S21中综合渗流方程，建立多系统全隐式单元特性矩阵为：

对离散区域内的每一个四面体网格，建立多系统全隐式单元特性矩阵；假设存在N个网格节点，通过整体叠加，形成关于有机质、无机质和裂缝系统的3N×3N大矩阵，即综合渗流数学模型的全隐式数值模型：

[K]_3N*3N[δX]_3N*1＝[R]_3N*1 (22)

式中符号定义：K为系数矩阵；δX为未知变量矩阵；R为余量矩阵；

S4：对内边界条件进行处理，建立基于四面体网格离散的全隐式多级压裂水平井数值井模型，并嵌入步骤S3的全隐式数值模型，构建完整的多系统同步联立求解矩阵；

S41：对于式(4)中的标况下的压裂井组地面产气量q_gsct，采用修正的拟稳态计算公式展开得到：

其中：

式中符号定义：N_well为压裂水平井井数；N_Fi为i压裂井的水力压裂裂缝总条数；p_bhi为i压裂井井底流压，MPa；p_ave，ij为i压裂井第j条压裂裂缝井源汇点所在网格块的平均压力，MPa；PI_g,ij为i压裂井第j条压裂裂缝采气指数，方/天/MPa；r_wi为i压裂井的井筒半径，m；r_0,ij为i压裂井第j条压裂裂缝井源汇点的等效井半径，m；

S42：对于压裂裂缝井等效井半径r_0,ij，将气体向压裂井的流动视为沿裂缝面向厚度为裂缝宽度w_F，ij的等效直井的拟稳态径向流动，则各裂缝面等效井半径表示为：

式中符号定义：A_ij为i压裂井第j条压裂裂缝井源汇点二维裂缝面的面积，m²；

S43：将式(18)带入式(23)中，得到全隐式压裂水平井的数值井计算格式为：

S44：将式(25)带入到式(22)中，构建完整的页岩储层和压裂井组多系统同步联立求解矩阵；

S5：对构建的多系统同步联立求解矩阵进行迭代求解，并分析朗格缪尔体积、压裂裂缝渗透率和压裂井开采模式对页岩气藏压裂井组生产动态的影响。

2.根据权利要求1所述的一种页岩气藏压裂井组生产动态数值模拟方法，其特征在于，所述步骤S1还包括以下步骤：

S11：根据研究工区的实际地质情况和水平井井轨迹及压裂裂缝分布特征，生成地质体；

S12：按照开源四面体网格剖分软件Tetgen的数据结构要求，对地质体数据进行编辑导入，生成网格离散模型。

3.根据权利要求1所述的一种页岩气藏压裂井组生产动态数值模拟方法，其特征在于，所述步骤S5还包括以下步骤：

S51：由有限元基本原理，对于封闭外边界条件，添加0向量到式(22)系数矩阵K中；

S52：采用共轭梯度法对完整的页岩储层和压裂井组多系统同步联立矩阵进行求解，通过牛顿-拉夫逊迭代方法，获得有机质、无机质以及裂缝系统中各离散网格节点在一个时步下的压力变化值：δp_gf1,…δp_gfN；δp_gI1,…δp_gI1N；δp_gO1,…δp_gO1N；进而由式(18)获得下一时步的值；

S53：输出各时步下的离散网格节点压力值，以及定井底流压生产条件下的压裂井组日产气量和累产气量。

Images