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

劉輝 廖新維 李則驊 Wei Yu

目前壓裂直井或壓裂水平井的單井模擬主要是基于二維或擬三維網(wǎng)格進行表征，但這2類網(wǎng)格無法精細(xì)表征裂縫、井筒的三維空間形態(tài)。為此，提出利用計算機輔助設(shè)計（CAD）技術(shù)來進行壓裂井的三維幾何建模方法，形成了3D精細(xì)半結(jié)構(gòu)網(wǎng)格剖分方法和流程，以及高效數(shù)值求解方法，將該方法應(yīng)用于現(xiàn)場一口壓裂直井，并對其進行了模擬分析。分析結(jié)果表明：所構(gòu)建的3D壓裂直井、壓裂水平井的內(nèi)邊界井模型，簡化和統(tǒng)一了不同壓裂井模型的井產(chǎn)量計算方程，該方法同樣適用于更復(fù)雜的3D定向壓裂井的井產(chǎn)量計算；水力裂縫的幾何形態(tài)、尺寸和小層物性是空間壓力分布的重要影響因素。所得結(jié)論可為提高復(fù)雜壓裂井的三維精細(xì)表征能力和單井模擬精度提供技術(shù)支撐。

壓裂井；精細(xì)表征；3D半結(jié)構(gòu)網(wǎng)格；單井模擬；多層縫網(wǎng)；壓力分布

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 引 言

以頁巖/致密儲層為代表的非常規(guī)儲層是提高產(chǎn)量的主力軍。在進行大型壓裂后，這類儲層通常具有強非均質(zhì)性的縫網(wǎng)，精確的裂縫表征對指導(dǎo)后續(xù)的開發(fā)工作具有重要意義。目前主要的裂縫表征模型有3類［1］：連續(xù)介質(zhì)、雙重介質(zhì)和離散介質(zhì)模型。連續(xù)介質(zhì)模型是將裂縫滲透率平均到基質(zhì)儲層中，整個儲層的滲透率是基質(zhì)和裂縫滲透率的平均值，無法準(zhǔn)確地刻畫復(fù)雜縫網(wǎng)的幾何拓?fù)浜碗x散性。雙重介質(zhì)模型將儲層劃分為基質(zhì)和裂縫2套獨立的系統(tǒng)，能較為準(zhǔn)確地反映裂縫的影響，但無法準(zhǔn)確表征具有強非均質(zhì)性和復(fù)雜幾何拓?fù)浣Y(jié)構(gòu)的縫網(wǎng)，一般僅適用于裂縫發(fā)育、分布均勻且連通性好的裂縫性儲層。離散介質(zhì)模型通過顯式表征裂縫的方式，可以更加精確地描述單個裂縫的物理特性以及具有復(fù)雜幾何拓?fù)浣Y(jié)構(gòu)的縫網(wǎng)，近年來被廣泛應(yīng)用于復(fù)雜裂縫性儲層的數(shù)值模擬。

離散介質(zhì)模型在物理域中不考慮裂縫的開度，將N維空間中的裂縫降為N-1維，從而使復(fù)雜縫網(wǎng)在數(shù)值計算中能被高效、精確處理。截至目前，眾多學(xué)者提出了各種不同的離散介質(zhì)模型［2-7］，這些模型主要可以分為2大類：嵌入式離散介質(zhì)模型（embedded discrete fracture model，EDFM）和離散基質(zhì)-裂縫介質(zhì)模型（discrete fracture-matrix model，DFM）。EDFM是將裂縫處理成點源或面源再嵌入到基質(zhì)網(wǎng)格中，而DFM是通過基質(zhì)網(wǎng)格的單元面來表征裂縫。EDFM通常是基于結(jié)構(gòu)正交化網(wǎng)格，求解效率高但精度較低；而DFM主要是用非結(jié)構(gòu)化網(wǎng)格進行表征，求解效率較低但精度更高。因此，EDFM一般應(yīng)用于縫網(wǎng)儲層的整體滲流特征的模擬，而DFM一般應(yīng)用于需要精細(xì)刻畫具有復(fù)雜幾何拓?fù)浣Y(jié)構(gòu)的縫網(wǎng)儲層來精細(xì)模擬縫網(wǎng)的局部滲流特征。

非常規(guī)儲層經(jīng)過大型壓裂處理后，在井周圍會形成由人工裂縫和天然裂縫相互交織組成的強非均質(zhì)性的復(fù)雜縫網(wǎng)，需要采用精細(xì)開發(fā)策略來提高局部區(qū)塊的產(chǎn)量水平。單井模擬分析方法是對單井范圍或井組范圍的儲層進行精細(xì)建模和模擬，相比油藏尺度的數(shù)值模擬，能更精細(xì)地刻畫近井裂縫特征和模擬其對產(chǎn)量的影響，對于提高強非均質(zhì)性的非常規(guī)儲層精細(xì)開發(fā)效果具有重要意義。

總結(jié)現(xiàn)有的相關(guān)文獻研究［8-18］，針對壓裂井的單井模擬研究方法，無論是直井還是水平井主要還是基于二維（2D）或擬三維（2.5D）模型，主要原因是建立三維（3D）精細(xì)幾何模型和3D精細(xì)網(wǎng)格模型的難度大。這些2D和2.5D模型無法精細(xì)表征三維空間中的任意傾角、高度的裂縫，以及裂縫與井筒復(fù)雜的相交方式，因此采用此模型既無法準(zhǔn)確求解真實壓裂井的空間滲流特征值，也不能滿足3D復(fù)雜壓裂井單井模擬的高精度需求。綜上，本文提出了利用計算機輔助設(shè)計（CAD）技術(shù)來實現(xiàn)壓裂井的三維幾何建模新方法，并形成了3D精細(xì)半結(jié)構(gòu)網(wǎng)格剖分方法和流程，以及高效數(shù)值求解方法。最后將所形成的方法應(yīng)用于現(xiàn)場多層復(fù)雜縫網(wǎng)儲層中的一口壓裂直井，對其井底瞬態(tài)壓力和多層復(fù)雜縫網(wǎng)空間滲流特征進行了分析。

1 滲流模型

1.1 物理模型

圖1為3D壓裂直井、3D壓裂水平井的物理模型示意圖。這里采用DFM方法來表征裂縫，在3D空間中，裂縫被降維表征成不規(guī)則幾何多邊形，并能以任意傾角與井筒相交。

1.2 數(shù)學(xué)模型

基質(zhì)系統(tǒng)為：

裂縫系統(tǒng)為：

式中：pm和pf分別為基質(zhì)和裂縫的壓力，Pa；cm為基質(zhì)和流體的綜合壓縮系數(shù)，1/Pa；cf為裂縫和流體的綜合壓縮系數(shù)，1/Pa；t為時間，s；Km和Kf分別為基質(zhì)和裂縫的滲透率，μm2；μ為流體黏度，Pa·s；φm和φf分別為基質(zhì)和裂縫的孔隙度，無量綱；af為裂縫開度，m；n為相鄰網(wǎng)格界面的單位法向量；q為基質(zhì)與裂縫之間的交換流量，m/s；Ω為積分域中的單元體；s為單元體Ω的表面積，m2；V為單元體Ω的體積，m3。

與源匯模式的Peaceman井模型相比，基于半結(jié)構(gòu)網(wǎng)格表征的壓裂井內(nèi)邊界井模型不需要對壓裂直井和壓裂水平井的井模型進行分類討論，可統(tǒng)一為同一種井產(chǎn)量計算方程。3D壓裂直井內(nèi)邊界井模型產(chǎn)量計算示意圖如圖2所示。

在圖2壓裂井模型中，井產(chǎn)量Qw為基質(zhì)和裂縫流向井底的流量之和，則可得計算方程：

式中：x為流體流動位置與井筒的距離，m;w為積分域；Sw為井筒截面積，m2。

等效井半徑與真實井半徑和表皮系數(shù)之間的關(guān)系為：

rwe=rw2.7-S（4）

式中：rwe和rw分別為等效井半徑和真實井半徑，m；S為表皮系數(shù)，無量綱。

耦合井筒儲集效應(yīng)的井產(chǎn)量Qp的計算公式為：

式中：Cw為井筒儲集系數(shù)，m3/Pa；pw為井底瞬態(tài)流壓，Pa。

2 3D精細(xì)半結(jié)構(gòu)網(wǎng)格表征模型

網(wǎng)格質(zhì)量會直接影響計算準(zhǔn)確性和效率。這里提出能適應(yīng)裂縫幾何形態(tài)和空間滲流分布特征的3D半結(jié)構(gòu)網(wǎng)格節(jié)點約束算法、網(wǎng)格剖分方法以及網(wǎng)格尺度優(yōu)化方法。

2.1 DFM半結(jié)構(gòu)網(wǎng)格節(jié)點約束算法

DFM半結(jié)構(gòu)網(wǎng)格節(jié)點約束算法是主要考慮徑向和線性分布的混合網(wǎng)格節(jié)點約束算法（見圖3），能精確地表征裂縫幾何特征、流態(tài)特征和井筒特征。

2.1.1 徑向分布網(wǎng)格節(jié)點約束算法

假設(shè)流體不可壓縮，單個徑向網(wǎng)格體積為V，則平均網(wǎng)格壓力為：

由裘比公式可得到徑向網(wǎng)格區(qū)域的壓力p與r的關(guān)系為：

式中：k為徑向網(wǎng)格區(qū)域的網(wǎng)格滲透率，m2。

將式（8）代入式（7）中，引入網(wǎng)格尺寸增長因子α，可得：

基于式（9），可得到徑向網(wǎng)格區(qū)域內(nèi)的所有節(jié)點的坐標(biāo)。

2.1.2 線性分布網(wǎng)格節(jié)點約束算法

在直井長井筒、水平井或者裂縫兩側(cè)需要使用線性分布的網(wǎng)格來適應(yīng)其線性流動特征。沿著線性加密軌跡使用均勻線性網(wǎng)格節(jié)分分布，在線性加密軌跡兩側(cè)使用漸進線性網(wǎng)格節(jié)點分布。

均勻線性分布網(wǎng)格節(jié)點坐標(biāo)的控制方程為：

式中：x1、y1和x2、y2分別為線性加密軌跡兩個端點的坐標(biāo)值；id為線性軌跡的分割序號，且id=0，1，…，nd-1；nd為沿著線性軌跡的分割數(shù)目。

漸進線性分布網(wǎng)格節(jié)點lid的坐標(biāo)控制方程為：

lid=ηn-1l0（11）

式中：l0為軌跡的最內(nèi)層網(wǎng)格節(jié)點與軌跡的距離;n為沿漸進線方向加密次數(shù);η為漸進線性分布的網(wǎng)格尺度增長因子。

2.1.3 網(wǎng)格尺寸優(yōu)化方法

網(wǎng)格尺寸優(yōu)化是確保計算精度和效率的關(guān)鍵步驟。這里采用有限差分來計算模型的網(wǎng)格幾何尺度的空間梯度，并基于空間梯度的值來優(yōu)化整個模型的網(wǎng)格尺寸大小。

定義函數(shù)f（x，y，z）為：

式中：d（x，y，z）為點（x，y，z）到指定子區(qū)域的最近距離函數(shù)；κ為控制網(wǎng)格細(xì)化敏感度正常數(shù)。

函數(shù)f（x，y，z）在x、y和z這3個方向上的梯度范數(shù)為：

再定義一個與梯度有關(guān)的網(wǎng)格尺寸函數(shù)hx，y，z為：

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

式中：hmax是最大網(wǎng)格尺寸，ε為控制梯度與網(wǎng)格尺寸之間敏感度正常數(shù)。

使用上面的定義，當(dāng)點接近指定子區(qū)域時，距離函數(shù)d（x，y，z）的值較小，則f（x，y，z）的值較大，反之f（x，y，z）的值逐漸減小。這樣，可以基于f（x，y，z）的幾何尺度空間梯度值來控制網(wǎng)格尺寸，確保在關(guān)鍵區(qū)域（接近子區(qū)域）獲得更細(xì)的網(wǎng)格，在非關(guān)鍵區(qū)域（遠(yuǎn)離子區(qū)域）獲得粗網(wǎng)格，減少網(wǎng)格數(shù)量。

2.2 3D半結(jié)構(gòu)網(wǎng)格系統(tǒng)

3D半結(jié)構(gòu)網(wǎng)格剖分前的關(guān)鍵一步是通過CAD技術(shù)中的布爾運算［19-20］（如交集、并集和差集，如圖4所示）來組合、修改或裁剪原始幾何形狀，明確復(fù)雜的空間幾何關(guān)系，從而建立精確的3D幾何模型。

壓裂直井和壓裂水平井的半結(jié)構(gòu)網(wǎng)格剖分示意圖如圖5所示。明確了裂縫、基質(zhì)、井筒之間的空間幾何關(guān)系后，基于前面介紹的網(wǎng)格節(jié)點約束算法和網(wǎng)格尺寸控制方法，在模型的各個區(qū)域確定適當(dāng)?shù)木W(wǎng)格節(jié)點類型和密度?？紤]到裂縫的幾何特征、與井筒的相交特征以及流態(tài)分布特征，先在關(guān)鍵區(qū)域（例如井壁、裂縫與井筒交接處以及近井區(qū)域流態(tài)變化劇烈的區(qū)域）放置密度較大的線性或徑向分布的網(wǎng)格節(jié)點，然后在遠(yuǎn)井區(qū)域放置較稀疏的網(wǎng)格節(jié)點。這些節(jié)點可以確保在后續(xù)網(wǎng)格生成過程中能高效、準(zhǔn)確地捕捉到各區(qū)域幾何和物理特征?；谏鲜銮闆r放置的網(wǎng)格節(jié)點，使用Delaunay三角網(wǎng)格剖分算法［21］來生成3D四面體或六面體網(wǎng)格。在生成網(wǎng)格后，需要重點檢查不同區(qū)域之間的網(wǎng)格尺寸和密度的相接情況，要確保不同的網(wǎng)格區(qū)域之間平滑過渡，避免形成大的尺寸跳躍，以確保數(shù)值穩(wěn)定性和精度。圖5b、圖5g分別為3D壓裂直井（見圖5a）和3D壓裂水平井（見圖5f）的CAD幾何模型。圖5c、圖5d和圖5h、圖5i分別為壓裂直井和壓裂水平井關(guān)鍵區(qū)域的半結(jié)構(gòu)網(wǎng)格剖分示意圖。圖5e、圖5j分別為壓裂直井和壓裂水平井的全局網(wǎng)格剖分示意圖。

3 模型求解

現(xiàn)有的裂縫模型普遍采用單一的兩點通量近似方法（TPFA）或多點通量近似方法（MPFA）來進行求解。這樣雖能明顯降低實現(xiàn)模型求解的難度，但單一的求解方法各自存在一定的局限性，即對于具有復(fù)雜幾何形狀或高度各向異性的裂縫模型，TPFA可能導(dǎo)致顯著的誤差，而MPFA雖計算精度較高，但會顯著增加計算成本。為此，考慮基質(zhì)和裂縫的不同滲流特征，在同一個裂縫模型使用混合離散求解方法，即在各向異性強的基質(zhì)系統(tǒng)采用MPFA方法，在流動路徑較明確、離散網(wǎng)格較規(guī)則的高導(dǎo)流裂縫系統(tǒng)中采用TPFA方法。

在同一求解系統(tǒng)中，耦合MPFA-TPFA的混合離散求解方法能在不損失求解精度的前提下，顯著提高復(fù)雜縫網(wǎng)模型的求解效率。模型求解示意圖如圖6所示。

3.1 基質(zhì)系統(tǒng)

采用MPFA方法來計算基質(zhì)系統(tǒng)傳導(dǎo)率。圖6a是一個網(wǎng)格交互區(qū)內(nèi)的4個相鄰六面網(wǎng)格體的俯視圖，網(wǎng)格體的相交線匯合交于點O。i為網(wǎng)格體編號，j為相鄰網(wǎng)格的交界面，i+1表示順時針方向上網(wǎng)格i的下一個網(wǎng)格，j+1表示順時針方向上交界面j的下一個交界面。Pi為網(wǎng)格i的平均壓力，Pi，j為網(wǎng)格i和網(wǎng)格i-1交界面的中心點的壓力，nj為該交界面j的單位法向量，Aj為交界面的面積。vi，j為網(wǎng)格i中心點到交界面j的面中心的向量，vi，j⊥v′i，j且vi，j=v′i，j。

假設(shè)網(wǎng)格內(nèi)壓力線性連續(xù)分布，壓力的勢梯度為SymbolQC@P，可得：

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

非零向量v′i，j和非零向量v′i，j-1線性無關(guān)，可作為網(wǎng)格i中的一組基向量。假設(shè)未知數(shù)α和β，可得：

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

假設(shè)變量φ=vi，jv′i，j-1=vi，j-1v′i，j，聯(lián)立式（15）和式（16），可計算從網(wǎng)格i通過界面j+1的流量ui，j+1為：

定義變量ξ+i，j和變量ξ-i，j分別為：

基于式（17）和式（18），界面j+1的兩個相鄰網(wǎng)格i和i+1通過該界面流量分別為：

由界面上的連續(xù)性條件可得：

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(yù)=（ai，j），B=（bi，j），C=（ci，j），D=（di，j），且

聯(lián)立式（19）和式（20）整理得到：

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

定義Q=（q1，q2，…，qN）T為交互區(qū)內(nèi)所有相鄰網(wǎng)格的交界面的流量q組成的矩陣;Tm為相鄰網(wǎng)格交互區(qū)內(nèi)所有交界面的傳導(dǎo)率組成的系數(shù)矩陣，可得：

Q=TmP（23）

聯(lián)立式（22）和式（23），則基質(zhì)區(qū)域的傳導(dǎo)率矩陣為：

3.2 裂縫系統(tǒng)

從圖6b可知，在裂縫系統(tǒng)中采用TPFA離散求解方法計算通過兩個相鄰網(wǎng)格單元的界面f處流量為：

式中：PE和PE′分別為網(wǎng)格單元E和E′的網(wǎng)格壓力;Tf為界面f的傳導(dǎo)系數(shù)。

裂縫系統(tǒng)中相鄰網(wǎng)格單元界面處的半傳導(dǎo)率Tfracf，e的計算公式為：

式中：e為該裂縫面的邊（在2D中表示節(jié)點）;af為裂縫開度;Af為界面f的面積;cf，e為連接裂縫邊e的中心點到界面f中心點的位置向量;nf，e為界面f的邊e的單位外法向量。

基于調(diào)和平均計算方法可以計算裂縫-裂縫的傳導(dǎo)系數(shù)Tfrac-frac f和裂縫-基質(zhì)之間的傳導(dǎo)系數(shù)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表示節(jié)點）時，忽略相交位置處的體積流量，則可得：

式中：qinter-frace→f為從裂縫流向相交位置處的流量，Pe為多個裂縫相交位置處的壓力。

最后，聯(lián)立式（1）、式（2）、式（5）、式（24）、式（27）、式（28）和式（29），采用全隱式離散方法和牛頓-辛普森迭代方法進行數(shù)值求解。

4 實例應(yīng)用

4.1 目標(biāo)區(qū)塊地質(zhì)特征

目標(biāo)區(qū)塊F位于鄂爾多斯盆地的X油藏中。目標(biāo)區(qū)塊F地質(zhì)信息和3D半結(jié)構(gòu)網(wǎng)格劃分如圖7所示。根據(jù)區(qū)塊F的沉積相和測井信息（見圖7a），可將該區(qū)塊細(xì)分為F-1層和F-2層。F-1層和F-2層的平均有效厚度分別為65和85 m。區(qū)塊F的碎屑成分中，長石質(zhì)量分?jǐn)?shù)為32.1%，石英為28.6%，巖屑為26.9%，平均滲透率為0.58 mD，平均孔隙度為8.5%，F(xiàn)-2小層物性要比F-1小層更好。近年的油田測壓數(shù)據(jù)顯示，平均地層壓力為18.41 MPa。

4.2 目標(biāo)區(qū)塊裂縫分布特征

通過巖心觀察、測井、各種應(yīng)力測試方法證明區(qū)塊F裂縫發(fā)育，其裂縫多以共軛“X”形節(jié)里和大角度斜交裂縫為主。目標(biāo)區(qū)塊的裂縫分布玫瑰圖（見圖7b）顯示，在北東和北西方向上存在一組“X”形裂縫帶，其中北東方向的裂縫發(fā)育強度超過北西方向，并且該“X”形裂縫帶穿透F-1與F-2小層。

測試井P為區(qū)塊F的一口采油直井，該井射穿區(qū)塊F的F-1小層和F-2小層。自2002年投產(chǎn)后，在不同垂深位置處進行過多次分層水力壓裂。該井于2018年壓裂完后投產(chǎn)并投入壓力計進行不關(guān)井壓力測試。井P位于“X”裂縫帶的北方位置。

4.3 3D多層復(fù)雜縫網(wǎng)單井模型

在對區(qū)塊F的儲層特征、裂縫分布特征以及井P壓裂情況認(rèn)識的基礎(chǔ)上，構(gòu)建了區(qū)塊F的地質(zhì)模型（見圖7c）和3D多層復(fù)雜縫網(wǎng)模型的半結(jié)構(gòu)網(wǎng)格表征模型（見圖7d、圖7e），并結(jié)合基礎(chǔ)地質(zhì)參數(shù)，預(yù)估了“X”形裂縫帶和隨機縫網(wǎng)的屬性參數(shù)，最后建立了井P的3D多層單井模擬模型。通過擬合井P的井底測試壓力曲線，經(jīng)數(shù)值反演得到一套3D縫網(wǎng)參數(shù)，如表1所示。

4.4 井P單井模擬結(jié)果分析

井P的單井模擬結(jié)果分析如圖8所示。圖8a為井P的試井曲線擬合結(jié)果與流動階段劃分情況。根據(jù)其壓力導(dǎo)數(shù)曲線變化特征，可劃分為5個階段，分別為：井儲和表皮效應(yīng)、早期3D徑向流、裂縫線性流、擬徑向流和擬邊界效應(yīng)。

圖8b為生產(chǎn)時間t=10 h時的井筒和徑向加密區(qū)域的壓力分布特征。從圖8b可以看出，其在垂向上具有明顯的差異性，下半部分的壓力明顯低于其上半部分。圖8c為生產(chǎn)時間t=1 000 h時的水力裂縫的壓力分布。由于F-2小層物性較好，流體能快速通過基質(zhì)流向高導(dǎo)流水力裂縫，進而流向井筒，導(dǎo)致F-2小層的水力裂縫的平均壓力明顯高于F-1小層中的水力裂縫。圖8d為生產(chǎn)時間t=10 h近井區(qū)域3D空間內(nèi)的壓力傳播特征。從圖8d可直觀地看出，壓力首先是沿著水力裂縫表面快速突進，這是壓力波在近井區(qū)域擴散的主方向。其空間壓力分布特征受水力裂縫幾何形狀、幾何尺寸和傾角的影響很大。圖8e為生產(chǎn)時間t=500 h時的F-2小層的“X”形裂縫帶的壓力分布特征。從圖8e可見，在“X”形裂縫帶的交點位置附近，北西方向裂縫帶的壓力明顯要低于北東方向的裂縫帶。圖8f為生產(chǎn)時間t=1 000 h時的井P模擬區(qū)域內(nèi)的空間壓力分布特征。從圖8f可以看出，由于裂縫帶的阻擋作用，北東方向的裂縫帶在東南區(qū)域還存在明顯的高壓區(qū)域，表明該區(qū)域的儲量動用量很少，這可為下一步優(yōu)化新井的井位提供指導(dǎo)。

5 結(jié) 論

（1）建立了一種新的壓裂直井、壓裂水平井的3D高質(zhì)量半結(jié)構(gòu)網(wǎng)格節(jié)點約束算法、網(wǎng)格剖分以及尺度優(yōu)化方法。在基質(zhì)和裂縫求解系統(tǒng)中采用混合離散求解方法，能有效提高含大量裂縫節(jié)點的裂縫模型的計算精度與效率的平衡比。

（2）基于精細(xì)半結(jié)構(gòu)網(wǎng)格來表征井壁的方法，構(gòu)建了3D壓裂直井、壓裂水平井的內(nèi)邊界井模型，簡化和統(tǒng)一了不同壓裂井模型的井產(chǎn)量計算方程，該方法同樣適用于更復(fù)雜的3D定向壓裂井的產(chǎn)量計算。

（3）實例應(yīng)用分析結(jié)果表明，壓裂井的早期空間壓力傳播特征與小層物性、裂縫的幾何尺寸和與井的相交角度有很大關(guān)系，在垂向上，其井壁壓力也存在明顯的非均質(zhì)性。此外，裂縫帶的走向也會對儲層壓力分布產(chǎn)生重要影響。

［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］ 姚軍，王子勝，張允，等.天然裂縫性油藏的離散裂縫網(wǎng)絡(luò)數(shù)值模擬方法［J］.石油學(xué)報，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的致密油藏分段壓裂水平井?dāng)?shù)值模擬［J］.西南石油大學(xué)學(xué)報（自然科學(xué)版），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］ 苗芷芃，吳濤，張榮軍，等.致密油體積壓裂復(fù)雜縫網(wǎng)形成規(guī)律數(shù)值模擬研究［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］ 孫賀東，歐陽偉平，朱松柏，等.多尺度離散裂縫性致密砂巖氣藏數(shù)值試井新方法：以塔里木盆地克拉蘇氣田為例［J］.天然氣工業(yè)，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］ 趙玉龍，黃鑫，張烈輝，等.基于嵌入式離散裂縫模型優(yōu)化的海陸過渡相頁巖氣壓裂水平井?dāng)?shù)值模擬［J］.天然氣工業(yè)，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］ 范明濤，李社坤，李軍，等.頁巖氣體積壓裂井筒組合體力學(xué)新模型研究［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］ 朱大偉，胡永樂，崔明月，等.局部網(wǎng)格加密與嵌入式離散裂縫模型耦合預(yù)測壓裂改造井產(chǎn)能［J］.石油勘探與開發(fā)，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］ 張守江，喬紅軍，張永飛，等.多翼裂縫垂直井的瞬態(tài)壓力分析［J］.西南石油大學(xué)學(xué)報（自然科學(xué)版），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.