壓裂井的三維精細表征與單井模擬分析

劉輝 廖新維 李則驊 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所示。明確了裂縫、基質、井筒之間的空間幾何關系后，基于前面介紹的網格節點約束算法和網格尺寸控制方法，在模型的各個區域確定適當的網格節點類型和密度?？紤]到裂縫的幾何特征、與井筒的相交特征以及流態分布特征，先在關鍵區域（例如井壁、裂縫與井筒交接處以及近井區域流態變化劇烈的區域）放置密度較大的線性或徑向分布的網格節點，然后在遠井區域放置較稀疏的網格節點。這些節點可以確保在后續網格生成過程中能高效、準確地捕捉到各區域幾何和物理特征?；谏鲜銮闆r放置的網格節點，使用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.