压裂井的三维精细表征与单井模拟分析

刘辉 廖新维 李则骅 Wei Yu

目前压裂直井或压裂水平井的单井模拟主要是基于二维或拟三维网格进行表征，但这2类网格无法精细表征裂缝、井筒的三维空间形态。为此，提出利用计算机辅助设计（CAD）技术来进行压裂井的三维几何建模方法，形成了3D精细半结构网格剖分方法和流程，以及高效数值求解方法，将该方法应用于现场一口压裂直井，并对其进行了模拟分析。分析结果表明：所构建的3D压裂直井、压裂水平井的内边界井模型，简化和统一了不同压裂井模型的井产量计算方程，该方法同样适用于更复杂的3D定向压裂井的井产量计算；水力裂缝的几何形态、尺寸和小层物性是空间压力分布的重要影响因素。所得结论可为提高复杂压裂井的三维精细表征能力和单井模拟精度提供技术支撑。

压裂井；精细表征；3D半结构网格；单井模拟；多层缝网；压力分布

3D Fine Characterization and Single-Well Simulation

Analysis of Fractured Wells

The single-well simulation for fractured vertical or horizontal wells is mainly conducted based on two-dimensional or quasi-three-dimensional grids，which，however，cannot accurately characterize the 3D spatial morphology of fractures and wellbores.In this paper，a CAD-based 3D geometric modeling method of fractured wells，a grid node constraint algorithm adapting to fracture morphology and spatial flow characteristics，and a grid-scale optimization method were proposed.Then，semi-structured grids were used to build a 3D inner boundary model of fractured well，simplifying and unifying the production calculation equations for fractured vertical and horizontal wells.Based on the differences in flow characteristics of matrix and fractures，a coupled two-point and multi-point flux approximation method was used for numerical solution.Finally，a single-well simulation was conducted on a fractured vertical well in a multilayer reservoir with complex fracture networks.The results show that the differences in physical properties of layers have a significant impact on the pressure distribution of hydraulic fractures in the layers，and the geometric shape and size of hydraulic fractures as well as the physical parameters of the layers are important factors affecting the spatial distribution of pressure.The study results are of practical significance for improving the 3D fine characterization capacity and single-well simulation accuracy of complex fractured wells.

fractured well；fine characterization；3D semi-structured grid；single-well simulation；multilayer fracture networks reservoir；spatial distribution of pressure

0 引 言

以页岩/致密储层为代表的非常规储层是提高产量的主力军。在进行大型压裂后，这类储层通常具有强非均质性的缝网，精确的裂缝表征对指导后续的开发工作具有重要意义。目前主要的裂缝表征模型有3类［1］：连续介质、双重介质和离散介质模型。连续介质模型是将裂缝渗透率平均到基质储层中，整个储层的渗透率是基质和裂缝渗透率的平均值，无法准确地刻画复杂缝网的几何拓扑和离散性。双重介质模型将储层划分为基质和裂缝2套独立的系统，能较为准确地反映裂缝的影响，但无法准确表征具有强非均质性和复杂几何拓扑结构的缝网，一般仅适用于裂缝发育、分布均匀且连通性好的裂缝性储层。离散介质模型通过显式表征裂缝的方式，可以更加精确地描述单个裂缝的物理特性以及具有复杂几何拓扑結构的缝网，近年来被广泛应用于复杂裂缝性储层的数值模拟。

离散介质模型在物理域中不考虑裂缝的开度，将N维空间中的裂缝降为N-1维，从而使复杂缝网在数值计算中能被高效、精确处理。截至目前，众多学者提出了各种不同的离散介质模型［2-7］，这些模型主要可以分为2大类：嵌入式离散介质模型（embedded discrete fracture model，EDFM）和离散基质-裂缝介质模型（discrete fracture-matrix model，DFM）。EDFM是将裂缝处理成点源或面源再嵌入到基质网格中，而DFM是通过基质网格的单元面来表征裂缝。EDFM通常是基于结构正交化网格，求解效率高但精度较低；而DFM主要是用非结构化网格进行表征，求解效率较低但精度更高。因此，EDFM一般应用于缝网储层的整体渗流特征的模拟，而DFM一般应用于需要精细刻画具有复杂几何拓扑结构的缝网储层来精细模拟缝网的局部渗流特征。

非常规储层经过大型压裂处理后，在井周围会形成由人工裂缝和天然裂缝相互交织组成的強非均质性的复杂缝网，需要采用精细开发策略来提高局部区块的产量水平。单井模拟分析方法是对单井范围或井组范围的储层进行精细建模和模拟，相比油藏尺度的数值模拟，能更精细地刻画近井裂缝特征和模拟其对产量的影响，对于提高强非均质性的非常规储层精细开发效果具有重要意义。

总结现有的相关文献研究［8-18］，针对压裂井的单井模拟研究方法，无论是直井还是水平井主要还是基于二维（2D）或拟三维（2.5D）模型，主要原因是建立三维（3D）精细几何模型和3D精细网格模型的难度大。这些2D和2.5D模型无法精细表征三维空间中的任意倾角、高度的裂缝，以及裂缝与井筒复杂的相交方式，因此采用此模型既无法准确求解真实压裂井的空间渗流特征值，也不能满足3D复杂压裂井单井模拟的高精度需求。综上，本文提出了利用计算机辅助设计（CAD）技术来实现压裂井的三维几何建模新方法，并形成了3D精细半结构网格剖分方法和流程，以及高效数值求解方法。最后将所形成的方法应用于现场多层复杂缝网储层中的一口压裂直井，对其井底瞬态压力和多层复杂缝网空间渗流特征进行了分析。

1 渗流模型

1.1 物理模型

图1为3D压裂直井、3D压裂水平井的物理模型示意图。这里采用DFM方法来表征裂缝，在3D空间中，裂缝被降维表征成不规则几何多边形，并能以任意倾角与井筒相交。

1.2 数学模型

基质系统为：

裂缝系统为：

式中：pm和pf分别为基质和裂缝的压力，Pa；cm为基质和流体的综合压缩系数，1/Pa；cf为裂缝和流体的综合压缩系数，1/Pa；t为时间，s；Km和Kf分别为基质和裂缝的渗透率，μm2；μ为流体黏度，Pa·s；φm和φf分别为基质和裂缝的孔隙度，无量纲；af为裂缝开度，m；n为相邻网格界面的单位法向量；q为基质与裂缝之间的交换流量，m/s；Ω为积分域中的单元体；s为单元体Ω的表面积，m2；V为单元体Ω的体积，m3。

与源汇模式的Peaceman井模型相比，基于半结构网格表征的压裂井内边界井模型不需要对压裂直井和压裂水平井的井模型进行分类讨论，可统一为同一种井产量计算方程。3D压裂直井内边界井模型产量计算示意图如图2所示。

在图2压裂井模型中，井产量Qw为基质和裂缝流向井底的流量之和，则可得计算方程：

式中：x为流体流动位置与井筒的距离，m;w为积分域；Sw为井筒截面积，m2。

等效井半径与真实井半径和表皮系数之间的关系为：

rwe=rw2.7-S（4）

式中：rwe和rw分别为等效井半径和真实井半径，m；S为表皮系数，无量纲。

耦合井筒储集效应的井产量Qp的计算公式为：

式中：Cw为井筒储集系数，m3/Pa；pw为井底瞬态流压，Pa。

2 3D精细半结构网格表征模型

网格质量会直接影响计算准确性和效率。这里提出能适应裂缝几何形态和空间渗流分布特征的3D半结构网格节点约束算法、网格剖分方法以及网格尺度优化方法。

2.1 DFM半结构网格节点约束算法

DFM半结构网格节点约束算法是主要考虑径向和线性分布的混合网格节点约束算法（见图3），能精确地表征裂缝几何特征、流态特征和井筒特征。

2.1.1 径向分布网格节点约束算法

假设流体不可压缩，单个径向网格体积为V，则平均网格压力为：

由裘比公式可得到径向网格区域的压力p与r的关系为：

式中：k为径向网格区域的网格渗透率，m2。

将式（8）代入式（7）中，引入网格尺寸增长因子α，可得：

基于式（9），可得到径向网格区域内的所有节点的坐标。

2.1.2 线性分布网格节点约束算法

在直井长井筒、水平井或者裂缝两侧需要使用线性分布的网格来适应其线性流动特征。沿着线性加密轨迹使用均匀线性网格节分分布，在线性加密轨迹两侧使用渐进线性网格节点分布。

均匀线性分布网格节点坐标的控制方程为：

式中：x1、y1和x2、y2分别为线性加密轨迹两个端点的坐标值；id为线性轨迹的分割序号，且id=0，1，…，nd-1；nd为沿着线性轨迹的分割数目。

渐进线性分布网格节点lid的坐标控制方程为：

lid=ηn-1l0（11）

式中：l0为轨迹的最内层网格节点与轨迹的距离;n为沿渐进线方向加密次数;η为渐进线性分布的网格尺度增长因子。

2.1.3 网格尺寸优化方法

网格尺寸优化是确保计算精度和效率的关键步骤。这里采用有限差分来计算模型的网格几何尺度的空间梯度，并基于空间梯度的值来优化整个模型的网格尺寸大小。

定义函数f（x，y，z）为：

式中：d（x，y，z）为点（x，y，z）到指定子区域的最近距离函数；κ为控制网格细化敏感度正常数。

函数f（x，y，z）在x、y和z这3个方向上的梯度范数为：

再定义一个与梯度有关的网格尺寸函数hx，y，z为：

hx，y，z=hmax-ε·‖fx，y，z‖（14）

式中：hmax是最大网格尺寸，ε为控制梯度与网格尺寸之间敏感度正常数。

使用上面的定义，当点接近指定子区域时，距离函数d（x，y，z）的值较小，则f（x，y，z）的值较大，反之f（x，y，z）的值逐渐减小。这样，可以基于f（x，y，z）的几何尺度空间梯度值来控制网格尺寸，确保在关键区域（接近子区域）获得更细的网格，在非关键区域（远离子区域）获得粗网格，减少网格数量。

2.2 3D半结构网格系统

3D半结构网格剖分前的关键一步是通過CAD技术中的布尔运算［19-20］（如交集、并集和差集，如图4所示）来组合、修改或裁剪原始几何形状，明确复杂的空间几何关系，从而建立精确的3D几何模型。

压裂直井和压裂水平井的半结构网格剖分示意图如图5所示。明确了裂缝、基质、井筒之间的空间几何关系后，基于前面介绍的网格节点约束算法和网格尺寸控制方法，在模型的各个区域确定适当的网格节点类型和密度。考虑到裂缝的几何特征、与井筒的相交特征以及流态分布特征，先在关键区域（例如井壁、裂缝与井筒交接处以及近井区域流态变化剧烈的区域）放置密度较大的线性或径向分布的网格节点，然后在远井区域放置较稀疏的网格节点。这些节点可以确保在后续网格生成过程中能高效、准确地捕捉到各区域几何和物理特征。基于上述情况放置的网格节点，使用Delaunay三角网格剖分算法［21］来生成3D四面体或六面体网格。在生成网格后，需要重点检查不同区域之间的网格尺寸和密度的相接情况，要确保不同的网格区域之间平滑过渡，避免形成大的尺寸跳跃，以确保数值稳定性和精度。图5b、图5g分别为3D压裂直井（见图5a）和3D压裂水平井（见图5f）的CAD几何模型。图5c、图5d和图5h、图5i分别为压裂直井和压裂水平井关键区域的半结构网格剖分示意图。图5e、图5j分别为压裂直井和压裂水平井的全局网格剖分示意图。

3 模型求解

现有的裂缝模型普遍采用单一的两点通量近似方法（TPFA）或多点通量近似方法（MPFA）来进行求解。这样虽能明显降低实现模型求解的难度，但单一的求解方法各自存在一定的局限性，即对于具有复杂几何形状或高度各向异性的裂缝模型，TPFA可能导致显著的误差，而MPFA虽计算精度较高，但会显著增加计算成本。为此，考虑基质和裂缝的不同渗流特征，在同一个裂缝模型使用混合离散求解方法，即在各向异性强的基质系统采用MPFA方法，在流动路径较明确、离散网格较规则的高导流裂缝系统中采用TPFA方法。

在同一求解系统中，耦合MPFA-TPFA的混合离散求解方法能在不损失求解精度的前提下，显著提高复杂缝网模型的求解效率。模型求解示意图如图6所示。

3.1 基质系统

采用MPFA方法来计算基质系统传导率。图6a是一个网格交互区内的4个相邻六面网格体的俯视图，网格体的相交线汇合交于点O。i为网格体编号，j为相邻网格的交界面，i+1表示顺时针方向上网格i的下一个网格，j+1表示顺时针方向上交界面j的下一个交界面。Pi为网格i的平均压力，Pi，j为网格i和网格i-1交界面的中心点的压力，nj为该交界面j的单位法向量，Aj为交界面的面积。vi，j为网格i中心点到交界面j的面中心的向量，vi，j⊥v′i，j且vi，j=v′i，j。

假设网格内压力线性连续分布，压力的势梯度为SymbolQC@P，可得：

Pi+1，j+1=Pi+vi，jP（15）

非零向量v′i，j和非零向量v′i，j-1线性无关，可作为网格i中的一组基向量。假设未知数α和β，可得：

P=αv′i，j+βv′i，j-1（16）

假设变量φ=vi，jv′i，j-1=vi，j-1v′i，j，联立式（15）和式（16），可计算从网格i通过界面j+1的流量ui，j+1为：

定义变量ξ+i，j和变量ξ-i，j分别为：

基于式（17）和式（18），界面j+1的两个相邻网格i和i+1通过该界面流量分别为：

由界面上的连续性条件可得：

ui，j+1=ui+1，j+1（20）

定义P=Pi，Pi+1，Pi+2，Pi-1，P′=（Pi，j，Pi+1，j+1，Pi+2，j+2，Pi-1，j-1），令A=（ai，j），B=（bi，j），C=（ci，j），D=（di，j），且

联立式（19）和式（20）整理得到：

AP+BP′=CP+DP′（22）

定义Q=（q1，q2，…，qN）T为交互区内所有相邻网格的交界面的流量q组成的矩阵;Tm为相邻网格交互区内所有交界面的传导率组成的系数矩阵，可得：

Q=TmP（23）

联立式（22）和式（23），则基质区域的传导率矩阵为：

3.2 裂缝系统

从图6b可知，在裂缝系统中采用TPFA离散求解方法计算通过两个相邻网格单元的界面f处流量为：

式中：PE和PE′分别为网格单元E和E′的网格压力;Tf为界面f的传导系数。

裂缝系统中相邻网格单元界面处的半传导率Tfracf，e的计算公式为：

式中：e为该裂缝面的边（在2D中表示节点）;af为裂缝开度;Af为界面f的面积;cf，e为连接裂缝边e的中心点到界面f中心点的位置向量;nf，e为界面f的边e的单位外法向量。

基于调和平均计算方法可以计算裂缝-裂缝的传导系数Tfrac-frac f和裂缝-基质之间的传导系数Tmatr-frac f分别为：

Tfrac-frac f=［（Tfrac f，E）-1+（Tfrac f，E′）-1）］-1（27）

Tmatr-frac f=［（Tfrac f，E）-1+（Tmatrf，E′）-1］-1（28）

对于多个裂缝（≥2）相交在同一条边e（在2D中，e表示节点）时，忽略相交位置处的体积流量，则可得：

式中：qinter-frace→f为从裂缝流向相交位置处的流量，Pe为多个裂缝相交位置处的压力。

最后，联立式（1）、式（2）、式（5）、式（24）、式（27）、式（28）和式（29），采用全隐式离散方法和牛顿-辛普森迭代方法进行数值求解。

4 实例应用

4.1 目标区块地质特征

目标区块F位于鄂尔多斯盆地的X油藏中。目标区块F地质信息和3D半结构网格划分如图7所示。根据区块F的沉积相和测井信息（见图7a），可将该区块细分为F-1层和F-2层。F-1层和F-2层的平均有效厚度分别为65和85 m。区块F的碎屑成分中，长石质量分数为32.1%，石英为28.6%，岩屑为26.9%，平均渗透率为0.58 mD，平均孔隙度为8.5%，F-2小层物性要比F-1小层更好。近年的油田测压数据显示，平均地层压力为18.41 MPa。

4.2 目标区块裂缝分布特征

通过岩心观察、测井、各种应力测试方法证明区块F裂缝发育，其裂缝多以共轭“X”形节里和大角度斜交裂缝为主。目标区块的裂缝分布玫瑰图（见图7b）显示，在北东和北西方向上存在一组“X”形裂缝带，其中北东方向的裂缝发育强度超过北西方向，并且该“X”形裂缝带穿透F-1与F-2小层。

测试井P为区块F的一口采油直井，该井射穿区块F的F-1小层和F-2小层。自2002年投产后，在不同垂深位置处进行过多次分层水力压裂。该井于2018年压裂完后投产并投入压力计进行不关井压力测试。井P位于“X”裂缝带的北方位置。

4.3 3D多层复杂缝网单井模型

在对区块F的储层特征、裂缝分布特征以及井P压裂情况认识的基础上，构建了区块F的地质模型（见图7c）和3D多层复杂缝网模型的半结构网格表征模型（见图7d、图7e），并结合基础地质参数，预估了“X”形裂缝带和随机缝网的属性参数，最后建立了井P的3D多层单井模拟模型。通过拟合井P的井底测试压力曲线，经数值反演得到一套3D缝网参数，如表1所示。

4.4 井P单井模拟结果分析

井P的单井模拟结果分析如图8所示。图8a为井P的试井曲线拟合结果与流动阶段划分情况。根据其压力导数曲线变化特征，可划分为5个阶段，分别为：井储和表皮效应、早期3D径向流、裂缝线性流、拟径向流和拟边界效应。

图8b为生产时间t=10 h时的井筒和径向加密区域的压力分布特征。从图8b可以看出，其在垂向上具有明显的差异性，下半部分的压力明显低于其上半部分。图8c为生产时间t=1 000 h时的水力裂缝的压力分布。由于F-2小层物性较好，流体能快速通过基质流向高导流水力裂缝，进而流向井筒，导致F-2小层的水力裂缝的平均压力明显高于F-1小层中的水力裂缝。图8d为生产时间t=10 h近井区域3D空间内的压力传播特征。从图8d可直观地看出，压力首先是沿着水力裂缝表面快速突进，这是压力波在近井区域扩散的主方向。其空间压力分布特征受水力裂缝几何形状、几何尺寸和倾角的影响很大。图8e为生产时间t=500 h时的F-2小层的“X”形裂缝带的压力分布特征。从图8e可见，在“X”形裂缝带的交点位置附近，北西方向裂缝带的压力明显要低于北东方向的裂缝带。图8f为生产时间t=1 000 h时的井P模拟区域内的空间压力分布特征。从图8f可以看出，由于裂缝带的阻挡作用，北东方向的裂缝带在东南区域还存在明显的高压区域，表明该区域的储量动用量很少，这可为下一步优化新井的井位提供指导。

5 结 论

（1）建立了一种新的压裂直井、压裂水平井的3D高质量半结构网格节点约束算法、网格剖分以及尺度优化方法。在基质和裂缝求解系统中采用混合离散求解方法，能有效提高含大量裂缝节点的裂缝模型的计算精度与效率的平衡比。

（2）基于精细半结构网格来表征井壁的方法，构建了3D压裂直井、压裂水平井的内边界井模型，简化和统一了不同压裂井模型的井产量计算方程，该方法同样适用于更复杂的3D定向压裂井的产量计算。

（3）实例应用分析结果表明，压裂井的早期空间压力传播特征与小层物性、裂缝的几何尺寸和与井的相交角度有很大關系，在垂向上，其井壁压力也存在明显的非均质性。此外，裂缝带的走向也会对储层压力分布产生重要影响。

［1］ BERRE I，DOSTER F，KEILEGAVLEN E.Flow in fractured porous media：A review of conceptual models and discretization approaches［J］.Transport in Porous Media，2019，130（1）：215-236.

［2］ LEE S H，LOUGH M，JENSEN C.Hierarchical modeling of flow in naturally fractured formations with multiple length scales［J］.Water resources research，2001，37（3）：443-455.

［3］ KARIMI-FARD M，DURLOFSKY L J，AZIZ K.An efficient discrete-fracture model applicable for general-purpose reservoir simulators［J］.SPE journal，2004，9（2）：227-236.

［4］ 姚军，王子胜，张允，等.天然裂缝性油藏的离散裂缝网络数值模拟方法［J］.石油学报，2010，31（2）：284-288.

YAO J，WANG Z S，ZHANG Y，HUANG Z Q.Numerical simulation method of discrete fracture network for naturally fractured reservoirs［J］.Editorial office of ACTA PETROLEI SINICA，2010，31（2）：284-288.

［5］ DONG Z，LI W，LEI G，et al.Embedded discrete fracture modeling as a method to upscale permeability for fractured reservoirs［J］.Energies，2019，12（5）：812.

［6］ 张烈辉，刘沙，雍锐，等.基于EDFM的致密油藏分段压裂水平井数值模拟［J］.西南石油大学学报（自然科学版），2019，41（4）：1-11.

ZHANG L H，LIU S，YONG R，et al.EDFM-based Numerical Simulation of Horizontal Wells with Multi-stage Hydraulic Fracturing in Tight Reservoirs［J］.Journal of Southwest Petroleum University （Science & Technology Edition），2019，41（4）：1-11.

［7］ LIU H，ZHAO X，TANG X，et al.A Discrete fracture-matrix model for pressure transient analysis in multistage fractured horizontal wells with discretely distributed natural fractures［J］.Journal of Petroleum Science and Engineering，2020，192：107275.

［8］ 苗芷芃，吴涛，张荣军，等.致密油体积压裂复杂缝网形成规律數值模拟研究［J］.石油机械，2022，50（12）：96-102.

MIAO Z P，WU T，ZHANG R J，QU L.Numerical simulation on formation law of complex fracture network by volume fracturing in tight oil reservoir［J］.China Petroleum Machinery，2022，50（12）：96-102.

［9］ 孙贺东，欧阳伟平，朱松柏，等.多尺度离散裂缝性致密砂岩气藏数值试井新方法：以塔里木盆地克拉苏气田为例［J］.天然气工业，2022，42（7）：55-64.

SUN H D，OUYANG W P，ZHU S B，et al.A new numerical well test method of multi-scale discrete fractured tight sandstone gas reservoirs and its application in the Kelasu Gas Field of the Tarim Basin［J］.Natural Gas Industry，2022，42（7）：55-64.

［10］ 赵玉龙，黄鑫，张烈辉，等.基于嵌入式离散裂缝模型优化的海陆过渡相页岩气压裂水平井数值模拟［J］.天然气工业，2023，43（4）：116-126.

ZHAO Y L，HUANG X，ZHANG L H，et al.Numerical simulation of fractured horizontal wells in transitional shale gas reservoirs based on embedded discrete fracture model optimization［J］.Natural Gas Industry，2023，43（4）：116-126.

［11］ LIU H，LIAO X，LIE K A，et al.An Integrated Modeling Framework for Simulating Complex Transient Flow in Fractured Reservoirs with 3D High-Quality Grids［C］∥SPE Reservoir Simulation Conference.OnePetro，2023［2023-06-20].

［12］ 范明涛，李社坤，李军，等.页岩气体积压裂井筒组合体力学新模型研究［J］.石油机械，2019，47（7）：63-70.

FAN M T，LI S K，LI J，et al.New mechanical model for wellbore assembly of shale gas well under stimulated reservoir volume fracturing［J］.China Petroleum Machinery，2019，47（7）：63-70.

［13］ WAN Y，LIU Y，WU N.Numerical pressure transient analysis for unfilled-caved carbonate reservoirs based on Stokes-Darcy coupled theory［J］.Journal of Petroleum Science and Engineering，2020，190：107085.

［14］ LIU H，LIAO X，TANG X，et al.A well test model based on embedded discrete-fracture method for pressure-transient analysis of fractured wells with complex fracture networks［J］.Journal of Petroleum Science and Engineering，2021，196：108042.

［15］ 朱大伟，胡永乐，崔明月，等.局部网格加密与嵌入式离散裂缝模型耦合预测压裂改造井产能［J］.石油勘探与开发，2020，47（2）：341-348.

ZHU D W，HU Y L，CUI M Y，et al.Productivity simulation of hydraulically fractured wells based on hybrid local grid refinement and embedded discrete fracture model［J］.Petroleum Exploration and Development，2020，47（2）：341-348.

［16］ QIN J，CHENG S，LI P，et al.Interference well-test model for vertical well with double-segment fracture in a multi-well system［J］.Journal of Petroleum Science and Engineering，2019，183：106412.

［17］ WANG H，KOU Z，GUO J，et al.A semi-analytical model for the transient pressure behaviors of a multiple fractured well in a coal seam gas reservoir［J］.Journal of Petroleum Science and Engineering，2021，198：108159.

［18］ 張守江，乔红军，张永飞，等.多翼裂缝垂直井的瞬态压力分析［J］.西南石油大学学报（自然科学版），2021，43（5）：137-146.

ZHANG S J，QIAO H J，ZHANG Y F，et al.Transient Pressure Analysis of Vertical Well with Multi-wing Fractures［J］.Journal of Southwest Petroleum University （Science & Technology Edition），2021，43（5）：137-146.

［19］ QIN Y，LUO Z，WEN L，et al.Research and application of Boolean operation for triangular mesh model of underground space engineering-Boolean operation for triangular mesh model［J］.Energy Science & Engineering，2019，7（4）：1154-1165.

［20］ GEUZAINE C，REMACLE J F.Gmsh：a 3D finite element mesh generator with built-in pre-and post-processing facilities［J］.International journal for numerical methods in engineering，2009，79（11）：1309-1331.

［21］ HYMAN J D，GABLE C W，PAINTER S L，et al.Conforming Delaunay triangulation of stochastically generated three dimensional discrete fracture networks：a feature rejection algorithm for meshing strategy［J］.SIAM Journal on Scientific Computing，2014，36（4）：A1871-A1894.