劉傳喜，方文超，秦學杰

（中國石化石油勘探開發(fā)研究院，北京 100083）

非常規(guī)油氣藏壓裂水平井進入生產(chǎn)期，壓裂縫網(wǎng)形態(tài)并不是靜態(tài)不變的，而是隨著開發(fā)過程的進行，縫內(nèi)壓力發(fā)生變化，表現(xiàn)出動態(tài)縫網(wǎng)特征［1-4］。2016年，Vermilye等［5］采用生產(chǎn)期微地震監(jiān)測手段，獲取一口壓裂水平井投產(chǎn)后由于應力和流體擾動所產(chǎn)生的微地震事件點，證實了縫網(wǎng)的動態(tài)變化特征，并提出了動態(tài)生產(chǎn)體積APV（active production volume）概念。監(jiān)測結(jié)果表明儲層中的壓裂改造范圍在生產(chǎn)期迅速縮小，2年后的APV約為初始改造體積的一半，3年后的APV幾乎消失。這與目前非常規(guī)油氣藏生產(chǎn)井初期產(chǎn)量相對較高但隨后快速遞減并維持長期低產(chǎn)的生產(chǎn)規(guī)律相吻合［6-7］。

壓裂縫網(wǎng)的動態(tài)變化對生產(chǎn)期儲層的供給能力和井的產(chǎn)量大小具有重要影響，是開發(fā)調(diào)整方案制定和實施的重要依據(jù)［8-10］。但目前尚未建立針對實際儲層中縫網(wǎng)動態(tài)變化過程的數(shù)值模擬方法，即目前用于油藏數(shù)值模擬的縫網(wǎng)模型保持為壓裂后的初始縫網(wǎng)形態(tài)，不隨儲層壓力發(fā)生變化。本研究在構(gòu)建壓裂縫網(wǎng)形變數(shù)學模型和數(shù)值模型的基礎(chǔ)上，采用傳導率修正技術(shù)將裂縫形變模型與離散裂縫數(shù)值模擬器進行耦合，研發(fā)了具有模擬動態(tài)縫網(wǎng)功能的數(shù)值模擬器，并開展實例應用。

1 裂縫形變數(shù)學模型

非常規(guī)油氣藏水平井壓裂形成縫網(wǎng)后，隨著開采進行，裂縫中流體壓力會產(chǎn)生較大波動，加之裂縫面上巖石應力的作用，導致開采過程中裂縫屬性（寬度和滲透率）發(fā)生變化。實驗表明［11-13］，壓裂裂縫寬度隨裂縫面上有效應力的變化遵循雙曲模型（圖1）。曲線的變化分為3個階段：隨著裂縫面上有效應力的增大（縫內(nèi)壓力降低），初始階段裂縫寬度呈非線性急劇下降；裂縫空間壓縮到一定程度后，支撐劑的抗壓強度發(fā)揮作用，此時繼續(xù)增大有效應力則使裂縫寬度呈緩慢擬線性下降；當裂縫空間壓縮逼近極限時，裂縫寬度對有效應力變化不敏感，此時寬度變化進入平緩段。3個階段在橫坐標上的分界點為σ1（有效應力，MPa）和σ2，其中σ2為裂縫寬度下降至極限點的極限有效應力，MPa；對應的寬度稱為殘余裂縫寬度wr，m；當有效應力大于σ2后，所建立的數(shù)學模型規(guī)定裂縫寬度不隨有效應力的變化而變化，保持為wr。

將圖1的雙曲模型用如下公式進行定量表征：

式中：wf為裂縫寬度，m；wmax為初始裂縫寬度，m；β為與裂縫力學性能相關(guān)的無因次特征參數(shù)；Ψ為裂縫剛度，MPa；σn為裂縫面上有效應力，MPa，其計算公式為：

式中：σnormal為裂縫面上的總正應力，MPa；pf為裂縫內(nèi)的流體壓力，MPa。

特征參數(shù)β的計算表達式為：

式中：θres為無因次殘余裂縫寬度系數(shù)，即殘余裂縫寬度與最大裂縫寬度的比值：

裂縫剛度Ψ與雙曲模型曲線的曲率密切相關(guān)，Ψ越大，曲率越小，意味著裂縫寬度下降的速度越小，曲線越接近線性，如圖2所示，圖中計算參數(shù)分別為σ2=20 MPa，θres=0.1。

圖2 裂縫變形雙曲模型曲線形態(tài)與裂縫剛度之間的關(guān)系Fig.2 Hyperbolic curves for fracture deformation under different fracture stiffness

在實際儲層中，沿裂縫走向的裂縫寬度分布并不是一個定值，而是呈“兩頭窄中間寬”的分布形態(tài)。為實現(xiàn)對裂縫變形更加真實的模擬，先對裂縫進行線剖分，然后采用二次函數(shù)對各個線單元的初始裂縫寬度wmax進行賦值（圖3）。

圖3 實際裂縫寬度的非均勻分布及其表征方法Fig.3 Non-uniformdistribution of fracture width along a real fracture and the characterization method of the fracture width distribution

將描述裂縫寬度隨應力變化的雙曲模型與描述裂縫寬度非均勻分布的二次函數(shù)模型進行結(jié)合，則不僅可模擬生產(chǎn)過程中裂縫寬度的變化，還可以模擬裂縫長度的變化，使裂縫變形模擬更貼近實際（圖4）。

圖4 模擬裂縫長度和寬度動態(tài)變化示意圖Fig.4 Schematic sketch of the dynamic fracture width and length modeled in the simulator

2 裂縫形變數(shù)值模擬方法

采用基于非結(jié)構(gòu)化網(wǎng)格的離散裂縫模型（DFM）開展裂縫形變模擬，離散裂縫模型在模擬過程中保持了裂縫原始分布和形態(tài)，對裂縫形變過程的模擬具有較好適應性。技術(shù)思路為：通過傳導率修正的方式，將裂縫形變數(shù)學模型與離散裂縫數(shù)值模擬模型進行耦合，形成具有模擬裂縫形變功能的離散裂縫數(shù)值模擬器。

2.1 離散裂縫模型傳導率計算

在DFM模型中，任意兩個相鄰網(wǎng)格i和j之間的傳導率可采用以下公式計算［14］：

式中：Tij為i和j之間的全傳導率，m3；Ti和Tj分別為網(wǎng)格i和網(wǎng)格j與網(wǎng)格交界面之間的半傳導率，m3。

上式即為DFM模型中常采用的兩點流動近似（TPFA）方法［15-16］。該方法分開計算不同類型介質(zhì)（基質(zhì)、裂縫）的傳導率，更易于觀察不同類型介質(zhì)中流動機理的不同。

單個網(wǎng)格與網(wǎng)格交界面之間半傳導率的計算公式為：

式中：Aij為網(wǎng)格i與j的網(wǎng)格交界面面積，m2；Ki為網(wǎng)格i滲透率，m2；d為網(wǎng)格交界面到網(wǎng)格中心點（基質(zhì)網(wǎng)格為三角形形心，裂縫網(wǎng)格為四邊形對角線交點0）垂直距離，m。對于裂縫網(wǎng)格：

如圖5所示，DFM模型的網(wǎng)格類型有兩種，分別為基質(zhì)網(wǎng)格（m）和裂縫網(wǎng)格（f），對應著兩種網(wǎng)格半傳導率Tm和Tf。網(wǎng)格的鄰接類型有3種，分別為基質(zhì)-基質(zhì)（①）、基質(zhì)-裂縫（②）和裂縫-裂縫（③），對應著3種全傳導率Tm-m，Tm-f和Tf-f。對裂縫形變進行模擬的核心方法就是利用裂縫形變數(shù)學模型對裂縫半傳導率進行修正，進而修正與裂縫網(wǎng)格相關(guān)的全傳導率，達到裂縫形變與流體流動耦合模擬的目的。

圖5 模擬器中2種網(wǎng)格和3種傳導率示意圖Fig.5 Schematic sketch of two kinds of grids and three kinds of transmissibility in the simulator

2.2 裂縫形變數(shù)值表征

將裂縫形變模型與裂縫網(wǎng)格半傳導率耦合，即將公式（1）代入公式（6）：

式中：Tf0為裂縫網(wǎng)格初始半傳導率，m3，計算式為：

式中：Ai，f表示網(wǎng)格i與裂縫網(wǎng)格f交界面面積，m2；Kf為裂縫網(wǎng)格滲透率，m2。考慮裂縫形變效應的裂縫半傳導率具有時變特征，公式（8）中σn與裂縫中的流體壓力密切相關(guān)。考慮裂縫變形效應的跨尺度耦合流動全傳導率計算仍基于TPFA方法［公式（5）］，此時，與裂縫網(wǎng)格相關(guān)的全傳導率（基質(zhì)-裂縫，裂縫-裂縫）也具備時變特征。即，基質(zhì)網(wǎng)格（m）-裂縫網(wǎng)格（f）流動全傳導率計算公式為：

裂縫網(wǎng)格（f1）-裂縫網(wǎng)格（f2）流動全傳導率計算公式為：

在模擬計算過程中，首先計算出只與網(wǎng)格幾何參數(shù)有關(guān)的裂縫網(wǎng)格初始半傳導率矩陣，而后在每一步模擬過程中對半傳導率和全傳導率進行修正。

3 動態(tài)縫網(wǎng)模擬實例應用

3.1 典型模型動態(tài)縫網(wǎng)模擬

建立致密油藏中一口多段壓裂水平井典型模型，如圖6所示，以直觀研究裂縫形變效應對油藏開發(fā)的影響。典型模型的尺寸為400 m×200 m，水平井長度設計為300 m，壓裂段數(shù)為7段，半縫長為60 m。

圖6 用于裂縫形變效應模擬的多段壓裂水平井網(wǎng)格模型Fig.6 Mesh model containing multi-stage fractures for the simulation of fracture deformation

表1總結(jié)了涉及到的模擬參數(shù)，模擬時間為365 d。裂縫應力參數(shù)設置為：殘余裂縫寬度系數(shù)θres=0.1，極限有效應力σ2=12 MPa，油藏總正應力在各網(wǎng)格間均勻分布σnormal=20 MPa。為模擬水力壓裂縫的縫寬在長度方向上的非均勻分布現(xiàn)象，對每一條裂縫的不同位置賦了不同的縫寬初值，每條裂縫的寬度呈現(xiàn)“兩頭窄中間寬”的分布規(guī)律，縫寬范圍為1～10 mm。

表1 多段壓裂水平井模型數(shù)值模擬參數(shù)Table 1 Simulation parameters of the multi-stage fractured well model

模型中設置的裂縫形變程度分為無裂縫形變、弱裂縫形變（Ψ=10 MPa）、強裂縫形變（Ψ=2 MPa）。利用模擬器計算得到的水平井產(chǎn)量和累計產(chǎn)量變化曲線如圖7所示。由圖可知，裂縫形變顯著降低了油井初產(chǎn)，其后一直保持低產(chǎn)水平，主要因為在衰竭式開采條件下，開井初期裂縫中的壓力快速下降，導致裂縫快速閉合，儲層供給能力急劇降低，從而極大降低了水平井產(chǎn)能，如圖8所示不同裂縫剛度對應的不同模擬時間點壓力分布圖。圖8表明，考慮裂縫形變效應后，裂縫的導流能力顯著降低，井周圍的壓降幅度明顯縮小。此外，對比可看出，Ψ=2 MPa時，裂縫導流能力的減弱首先發(fā)生在裂縫的兩端，且減弱程度最大。當裂縫中間處的導流能力處于相對較高的水平時，裂縫兩端的壓降程度仍然很小，這主要是因為裂縫兩端的寬度小，因此在壓力下降過程中優(yōu)先閉合。

圖7 不同裂縫剛度對應的水平井生產(chǎn)動態(tài)曲線Fig.7 Production dynamic curves obtained under different fracture stiffness

圖8 不同裂縫剛度對應的壓力分布變化Fig.8 Reservoir pressure variation under different fracture stiffness with time

圖9為不同時期縫網(wǎng)裂縫寬度分布圖，即動態(tài)縫網(wǎng)模型?？煽闯?，裂縫形變的過程由裂縫兩端向中間延伸。在生產(chǎn)初期，裂縫形態(tài)變化較劇烈，閉合速度較快；到了生產(chǎn)中后期，閉合速度變慢；模擬期末，裂縫兩端的縫寬基本接近極限寬度，這與圖1所示的裂縫寬度的雙曲變化規(guī)律相吻合。

圖9 數(shù)值模擬不同生產(chǎn)時期的動態(tài)縫網(wǎng)模型Fig.9 Dynamic fracture network models at different production periods in the numerical simulator

3.2 實際井動態(tài)縫網(wǎng)模擬

為研究油田實際井的動態(tài)縫網(wǎng)變化特征及其對井產(chǎn)量的影響，選取鄂爾多斯盆地某致密油藏EN125井開展動態(tài)縫網(wǎng)構(gòu)建和油藏數(shù)值模擬。EN125井水平段長780 m，壓裂段數(shù)7段，泵注總液量1 482 m3，總砂（陶粒）量176 m3，排量3 m3∕min，井周圍局部發(fā)育有天然裂縫。采用壓裂模擬方法預測了該井的壓裂縫網(wǎng)模型，如圖10a所示。

圖10b為EN125井離散裂縫數(shù)值模擬網(wǎng)格模型，基質(zhì)滲透率0.36×10-3μm2，基質(zhì)孔隙度8%，裂縫初始寬度范圍1～10 mm，原始油藏壓力18 MPa，流體粘度3.2 mPa·s，模型考慮了致密油藏啟動壓力梯度，實驗測試結(jié)果為0.18 MPa∕m。擬合后的壓裂裂縫形變參數(shù)：殘余裂縫寬度系數(shù)θres=0.1，裂縫剛度Ψ=2 MPa，極限有效應力σ2=12 MPa，油藏總正應力在各網(wǎng)格間均勻分布σn=18 MPa。

圖10 EN125井壓裂縫網(wǎng)模型（a）及對應的油藏數(shù)值模擬網(wǎng)格模型（b）Fig.10 Fracture network model（a）and corresponding reservoir numerical simulation mesh model（b）for Well EN125

圖11為模擬1年得到的EN125井在不同時期的動態(tài)縫網(wǎng)模型。可見，裂縫動態(tài)閉合主要發(fā)生在生產(chǎn)初期的1～2個月內(nèi)，到后期逐漸放緩。整個過程中動態(tài)縫網(wǎng)的變化規(guī)律為：主裂縫溝通的天然裂縫首先閉合，水力裂縫的寬度變??；而后水力裂縫由兩端向中間逐漸閉合。圖11表明，壓裂形成的初始縫網(wǎng)在1年后會發(fā)生顯著變化，縫網(wǎng)規(guī)模明顯變小，與文獻中的實際井監(jiān)測結(jié)果相符［5］。

圖11 EN125井不同生產(chǎn)時期的動態(tài)縫網(wǎng)模型Fig.11 Dynamic fracture network models at different production periods for Well EN125

為更直觀研究動態(tài)縫網(wǎng)對壓裂水平井產(chǎn)量的影響，采用離散裂縫數(shù)值模擬器對EN125井開展了兩種情況下的模擬。一種情況是考慮裂縫閉合，另一種情況是不考慮裂縫閉合，其他參數(shù)均一致。將兩種情況下的模擬產(chǎn)液量與該井實際產(chǎn)液量進行擬合，結(jié)果如圖12所示。

圖12表明，當不考慮裂縫閉合時，模擬產(chǎn)液量明顯高于井的實際產(chǎn)液量，特別是初期產(chǎn)量吻合率較低，因為此期間儲層中裂縫閉合對產(chǎn)量的影響并未在模擬器中體現(xiàn)。而在相同物性和流體參數(shù)條件下，當在模擬器中考慮裂縫閉合時，模擬產(chǎn)液量與實際產(chǎn)液量差別較小，特別是初期產(chǎn)量較容易擬合上。因此，動態(tài)縫網(wǎng)對致密油藏壓裂水平井的產(chǎn)量影響顯著，在數(shù)值模擬過程中考慮動態(tài)縫網(wǎng)更加符合致密油藏生產(chǎn)實際，開發(fā)方案設計精度也將更高。

4 結(jié)論

1）壓裂縫網(wǎng)是非常規(guī)油氣藏中的主要滲流通道，油氣井投產(chǎn)后，由于壓裂縫內(nèi)壓力發(fā)生變化，裂縫發(fā)生動態(tài)形變，表現(xiàn)出動態(tài)縫網(wǎng)特征，對油氣井產(chǎn)能具有重要影響。

2）建立了表征裂縫形變的雙曲數(shù)學模型及描述裂縫寬度非均勻分布的二次函數(shù)模型，可同時模擬裂縫寬度和長度的動態(tài)變化。

3）采用傳導率修正方式，將裂縫形變數(shù)學模型與離散裂縫數(shù)值模擬模型耦合，形成了具有模擬裂縫形變功能的離散裂縫數(shù)值模擬器。

4）建立了一口壓裂水平井典型模型和一口致密油藏實際壓裂水平井模型，應用模擬器進行模擬，構(gòu)建了不同生產(chǎn)時期的動態(tài)縫網(wǎng)，驗證了動態(tài)縫網(wǎng)模擬功能和可靠性。