胡之牮，李樹新，王建君，周鴻，趙玉龍，張烈輝

（1.西南石油大學油氣藏地質(zhì)及開發(fā)工程全國重點實驗室，四川 成都 610500；2.中國石油煤層氣有限責任公司，北京 100028；3.中國石油浙江油田分公司，浙江 杭州 310023；4.中國石油西南油氣田分公司蜀南氣礦，四川 瀘州 646001）

隨著非常規(guī)能源在世界能源中占據(jù)越來越重要的地位，頁巖氣作為“清潔非常規(guī)能源”的開發(fā)成為世界能源開發(fā)的研究重點。中國是頁巖氣大國，擁有豐富的頁巖氣資源量和廣闊的開采前景，并且已經(jīng)在四川盆地及其鄰區(qū)實現(xiàn)了海相頁巖氣資源的有效開發(fā)[1?2]。頁巖儲層具有超低孔、超低滲的特點，普通的開采方式不能充分將其開采，必須通過大型水力壓裂技術(shù)來實現(xiàn)工業(yè)化開采[3]，因此，頁巖氣的產(chǎn)能研究對頁巖氣藏的有效開發(fā)有著重要意義。

頁巖氣主要以游離態(tài)和吸附態(tài)方式賦存在納米孔隙中，并且在多尺度流動空間內(nèi)存在許多復雜流動機理。目前，主流的研究方法是根據(jù)克努森數(shù)劃分流動形態(tài)[4?5]，運用線性加權(quán)[6]或分子碰撞頻率加權(quán)[7?8]的耦合方式建立綜合考慮滑脫、解吸和擴散[9]的表觀滲透率模型，但還很少有研究者將頁巖氣藏孔隙分為無機質(zhì)孔隙和有機質(zhì)孔隙進行研究。國內(nèi)外預測頁巖氣藏壓裂水平井產(chǎn)能主要有解析法、半解析法和數(shù)值法。解析法[10?11]大多都考慮儲層為多重介質(zhì)，并劃分幾種流動區(qū)域，利用拉式變換求得產(chǎn)量的解析解；半解析法[12?13]通過離散裂縫網(wǎng)格，構(gòu)造壓降矩陣，求解線性方程組進行產(chǎn)能預測；數(shù)值法[14?15]利用有限元、有限差分等方法求解滲流偏微分方程，從而模擬生產(chǎn)動態(tài)。但水力壓裂形成的裂縫不僅形狀復雜，并且呈現(xiàn)“不等間距、不等裂縫半長[16]”的非均勻分布，研究者們的產(chǎn)能計算方法沒有全面綜合地考慮到這些，造成計算結(jié)果不準確。

研究基于頁巖氣的多重運移機制理論，將儲層分為無機質(zhì)孔隙和有機質(zhì)孔隙，建立耦合了2種孔隙的頁巖表觀滲透率模型（圖1）。在此基礎上，考慮壓裂裂縫的復雜形狀和非均勻分布，利用源匯思想和疊加原理建立頁巖氣藏壓裂復雜裂縫水平井產(chǎn)能預測模型，分析了頁巖氣微觀滲流、裂縫形狀和裂縫非均勻分布對產(chǎn)能的影響。

1 模型建立

1.1 頁巖氣表觀滲透率模型

頁巖氣儲層多為微米—納米級孔隙，孔隙結(jié)構(gòu)復雜、形狀多樣[17?18]。運用迂曲毛管束的分形理論[19]來描述頁巖儲層孔隙，將頁巖儲層分為有機質(zhì)孔隙和無機質(zhì)孔隙，全面考慮黏性流、克努森擴散和表面擴散的流動機理，推導頁巖表觀滲透率。

無機質(zhì)孔隙中存在黏性流[20]與克努森擴散?？紤]稀薄氣體效應修正黏性流公式，再引入加權(quán)因子表征黏性流和克努森擴散的貢獻率[8]，得到無機質(zhì)單根毛管截面的總質(zhì)量通量為：

式中：Jin為無機質(zhì)單根毛管截面的總質(zhì)量通量，單位mol/（m2·s）；εV為黏性流貢獻率；α為稀薄氣體效應系數(shù)；Kn為克努森數(shù)；D為孔隙直徑，單位m；p為地層壓力，單位MPa；T為地層溫度，單位K；Z為偏差因子；R為氣體常數(shù)，單位J/（mol·K）；Mg為氣體分子摩爾質(zhì)量，單位kg/mol；μg為氣體黏度，單位mPa·s；εN為克努森擴散貢獻率；Cg為氣體壓縮系數(shù)，單位1/MPa；Δp為毛管束兩端壓力差，單位MPa；L為迂曲毛細管長度，單位m。

式（1）乘以截面面積得到單根毛管質(zhì)量流量，再對單根毛管沿最小和最大孔徑進行積分，利用達西定律計算無機質(zhì)孔隙滲透率為：

式中：Kin為無機質(zhì)孔隙滲透率，單位10?3μm2；φ為孔隙度；Dp為孔隙面積分形維數(shù)；Dt為迂曲度分形維數(shù)；Dmax為最大孔隙直徑，單位m；Dmin為最小孔隙直徑，單位m。

有機質(zhì)孔隙表面存在吸附氣，不僅會減小有效孔徑，并且隨著儲層壓力降低，吸附氣發(fā)生的表面擴散還會對滲透率做出貢獻。氣體分子覆蓋率、有效孔徑和有效孔隙度[21]計算公式分別為：

式（3）—式（5）中：θ為氣體分子覆蓋率；pL為Langmuir壓力，單位MPa；De為有效孔徑，單位m；dm為氣體分子直徑，單位m；φe為有效孔隙度；Dave為平均孔隙直徑，單位m。

利用考慮吸附氣厚度的有效參數(shù)計算3 種復雜流動的流量，采用壁聯(lián)擴散[22]的方式進行耦合，得到有機質(zhì)孔隙滲透率為：

式中：Kor為有機質(zhì)孔隙滲透率，單位10?3μm2；Le為有機質(zhì)迂曲毛細管長度，單位m；DB為表面擴散系數(shù)，單位m2/s；NA為阿伏伽德羅常數(shù)，單位1/mol；下標e 代表有機質(zhì)有效參數(shù)。

通過對頁巖進行電鏡掃描，確定有機質(zhì)占比ε即可加權(quán)得出頁巖表觀滲透率為：式中：Kapp為頁巖表觀滲透率，單位10?3μm2；ε為頁巖儲層有機質(zhì)占比。

1.2 壓裂水平井產(chǎn)能模型

模型設定條件如下：①水平井筒位于氣藏中心，水平井長度方向為x軸，裂縫長度方向為y軸，裂縫高度方向為z軸；②裂縫垂直于水平井筒，縫寬和縫高沿裂縫長度方向減小，并且沿水平井長度方向呈現(xiàn)“不等間距、不等半長”的分布；③氣藏初始壓力均勻分布；④氣藏內(nèi)流動為單相氣體流動且忽略重力的影響；⑤氣體僅通過裂縫流入井筒并且不考慮井筒壓降。

定義氣體擬壓力[23]為：

式中：m（p）為氣體擬壓力，單位MPa；pi為原始地層壓力，單位MPa；μgi為初始氣體黏度，單位mPa·s；Zi為初始偏差因子。

根據(jù)高等滲流力學理論[24]，封閉地層點匯擬壓力形式為：

式中：m（pi）為氣體初始擬壓力，單位MPa；q（t）為點匯流量，單位m3/s；Ct為綜合壓縮系數(shù)，單位1/MPa；S(x,t)S(y,t)S(z,t)為x、y、z方向基本匯函數(shù)。

其中3個方向基本匯函數(shù)分別為：

式（10）—式（12）中：xe、ye、ze分別為儲層長度、寬度、厚度，單位m；xw、yw、zw為點匯坐標，單位m；x、y、z為氣藏任意一點坐標，單位m；t為時間，單位s。

將式（8）代入式（9），引入平均物性參數(shù)來近似等效氣體黏度和偏差因子的變化，可以得到真實氣體狀態(tài)下的點匯壓力解為：

由于地應力和壓裂施工等因素的影響，水力壓裂不會形成理想的矩形裂縫，而是沿著裂縫長度方向會變窄、變矮，常見的復雜形狀裂縫如圖2。

圖2 復雜形狀裂縫示意圖Fig.2 Schematic diagram of complex-shaped fracture

水平井上有n條裂縫，每條裂縫的上下兩翼沿中心線各離散成n份，每條裂縫共離散成2n份，整個縫網(wǎng)共有2n×n個裂縫網(wǎng)格，每個裂縫網(wǎng)格中心點可視為一個點匯（圖3）。

圖3 單條裂縫離散網(wǎng)格圖Fig.3 Grid discretization of single fracture

根據(jù)點匯函數(shù)，初始時刻各邊界并無流體流動，裂縫上翼微元對氣藏中任意一點產(chǎn)生的壓降為：

式中：p（x,y,z,t）為地層任意一點的壓力，單位MPa；為第j條裂縫上翼第i個微元流量，單位m3/s；上標（a）代表上翼方程，如果換成（b）代表下翼方程，上下兩翼方程具有鏡像性，因此，只給出上翼方程，下同。

定義氣藏壓降系數(shù)為：

式中：R(i,j)為氣藏壓降系數(shù)，單位MPa/m3；i為裂縫序號；j為單翼裂縫網(wǎng)格序號。

線性方程組共有2n×n個方程，把方程組分成n份，每一份的前n個方程代表這條裂縫的上翼，后n個方程代表這條裂縫的下翼。上翼總方程組如下：

式中：為第j條裂縫上翼第k個微元壓力，單位MPa；o、k、l、m為序號。

由達西定律可知，第j條裂縫上翼任意裂縫微元i到井筒的壓降為：

式中：pwf為井底流壓，單位MPa；psc為標準狀況下壓力，單位MPa；Tsc為標準狀況下溫度，單位K；Kf為裂縫滲透率，單位10?3μm2；第j條裂縫上翼長度，單位m；f函數(shù)和g函數(shù)分別為縫寬和縫高發(fā)育函數(shù)。

2 模型求解

聯(lián)立公式（15）和（16）得到耦合流動方程組，方程組由壓力平方差、壓降矩陣和裂縫微元流量構(gòu)成。方程如下：

式（18）—式（20）中：nt為時間步數(shù)；Δt為時間步長；為第nt?1個時間步數(shù)的壓力平方差矩陣；qntΔt為第nt個時間步數(shù)的裂縫流量為第nt?1個時間步數(shù)的地層壓力，單位MPa，當nt=1時，pntΔt=pi。

A為縫網(wǎng)系統(tǒng)總壓降矩陣，由氣藏壓降矩陣和裂縫壓降矩陣組成：

氣藏壓降矩陣如下：

裂縫壓降矩陣如下，O為零矩陣，代表第j條裂縫的壓降矩陣：

其中表示第j條裂縫上翼的裂縫壓降矩陣：

第j條裂縫上翼裂縫壓降矩陣系數(shù)為：

所求的方程為線性方程組，運用高斯—賽德爾迭代法即可求解。模型求解流程如圖4。

圖4 模型求解流程圖Fig.4 Model solving flow chart

3 模型驗證與分析

3.1 模型驗證

為驗證模型準確性，分別采用Petrel 數(shù)值模擬軟件和本文模型計算壓裂水平井產(chǎn)能，計算參數(shù)見表1。

表1 頁巖氣藏基礎參數(shù)Table 1 Basic parameters of shale gas reservoir

不考慮裂縫的復雜形狀和非均勻分布，模型退化為等間距壓裂水平井產(chǎn)能模型，模型計算結(jié)果與數(shù)值模擬軟件計算結(jié)果基本一致（誤差小于5.0%），驗證了模型的準確性（圖5）。

圖5 等間距壓裂水平井產(chǎn)能模型計算結(jié)果對比Fig.5 Comparison of calculated results for production capacity models of equally spaced fractured horizontal wells

3.2 參數(shù)敏感性分析

3.2.1 微觀滲流機理

頁巖氣藏與常規(guī)氣藏最大的不同是頁巖氣在微納米孔隙中存在多重運移機制。采用定壓生產(chǎn)進行模擬，由圖6a可知，1 500 d后不考慮微觀滲流的累計產(chǎn)氣量為0.98×108m3，考慮微觀滲流的累計產(chǎn)氣量達到了1.03×108m3，兩者相差4.85%；而在生產(chǎn)初期（即為剛開始生產(chǎn)時，t=1 d）考慮微觀滲流與不考慮微觀滲流的日產(chǎn)氣量差距達到了20.3 %。將式（3）變換成如下式子：

圖6 微觀滲流對頁巖氣產(chǎn)量的影響Fig.6 Influence of micro seepage on shale gas production

由式（28）結(jié)果可知：當?shù)貙訅毫^高時，氣體分子覆蓋率值較大，有機質(zhì)孔隙壁面充滿著吸附態(tài)氣體分子；隨著地層壓力下降，氣體分子覆蓋率值變小，吸附態(tài)氣體分子不斷地從壁面解吸出來。在生產(chǎn)初期，地層壓力下降較快，此時，地層開始大量出現(xiàn)解吸氣，當生產(chǎn)到300 d 時，地層壓力下降開始變得緩慢，吸附氣的解吸變得緩慢，因此，微觀滲流對頁巖氣井產(chǎn)量的影響不可忽略，需要構(gòu)建全面考慮氣體多重運移機制的頁巖表觀滲透率模型評價頁巖氣井的產(chǎn)能（圖6b）。

3.2.2 裂縫形狀

在其他參數(shù)不變的情況下，裂縫形狀的變化顯著影響水平井產(chǎn)量，所有復雜形狀裂縫產(chǎn)量都低于理想矩形裂縫，產(chǎn)量由大到小依次為：矩形裂縫、橢形裂縫、楔形裂縫、星形裂縫；模擬生產(chǎn)1 500 d后，與理想矩形裂縫產(chǎn)量相比，橢形裂縫、楔形裂縫、星形裂縫的累產(chǎn)氣量分別減少了2.04%、3.87%、6.63%；初始日產(chǎn)氣量減幅分別達到了4.92 %、9.63 %和17.01%。綜上所述，考慮裂縫為理想矩形裂縫不僅與實際情況不符，還高估了其產(chǎn)能（圖7）。

圖7 裂縫形狀對頁巖氣產(chǎn)量的影響Fig.7 Influence of fracture shape on shale gas production

3.2.3 裂縫非均勻分布

由于壓裂的成本與裂縫長度呈正相關，考慮裂縫總數(shù)與長度一定，需分析裂縫橫向非均勻分布和縱向非均勻分布對產(chǎn)能的影響，因此，建立裂縫縱、橫向非均勻分布模型（圖8）。

圖8 裂縫縱、橫向非均勻分布模型Fig.8 Models for non-uniform distribution of fractures in both longitudinal and transverse directions

由模擬結(jié)果（圖9）可知：在裂縫總數(shù)和長度不變時，橫向非均勻分布模式的產(chǎn)量由大到小依次為：外密內(nèi)疏型、均勻型、外疏內(nèi)密型。1 500 d后，外密內(nèi)疏型、均勻型、外疏內(nèi)密型的累產(chǎn)氣量分別為1.00×108m3、0.98×108m3、0.95×108m3，相較于均勻型，外密內(nèi)疏型的累產(chǎn)氣量增加了2.04 %，外疏內(nèi)密型的累產(chǎn)氣量減少了5.26 %。這是由于外疏內(nèi)密型會造成中間裂縫的控制區(qū)域重疊，造成浪費?？v向非均勻分布模式中，交錯型的累產(chǎn)氣量最高，為1.01×108m3；紡錘型獲得的累產(chǎn)氣量次之，為0.99×108m3；啞鈴型的產(chǎn)量最低，為0.97×108m3。相較于均勻型，交錯型的累產(chǎn)氣量增加了3.06%，紡錘型的累產(chǎn)氣量增加了1.02 %，啞鈴型的累產(chǎn)氣量減少了1.03 %。分析可得，相較紡錘型布局和啞鈴型布局，交錯型布局的裂縫控制面積最大。綜上所述，在頁巖氣藏的壓裂施工中選擇外密內(nèi)疏?交錯型的布縫方式可以提高水平井的產(chǎn)能。

圖9 裂縫非均勻分布對頁巖氣產(chǎn)量的影響Fig.9 Influence of non-uniform distribution of fractures on shale gas production

3.3 實例應用

選取長寧頁巖氣藏一口水平井進行生產(chǎn)數(shù)據(jù)擬合，基礎計算參數(shù)見表2。

表2 長寧頁巖氣井基礎參數(shù)Table 2 Basic parameters of Changning shale gas well

全面考慮頁巖氣的微觀滲流機理和壓裂裂縫的復雜拓展，模擬生產(chǎn)時間為20 a。由圖10 可以看出，本研究模型計算結(jié)果與生產(chǎn)數(shù)據(jù)接近，20 a后的累產(chǎn)氣量為1.33×108m3，運用文獻[25]中模型預測該井的EUR（最終可采儲量）為1.36×108m3，兩者相差2.21%，再次證明模型可靠。

圖10 生產(chǎn)數(shù)據(jù)擬合圖Fig.10 Production data fitting diagram

4 結(jié)論

1）基于氣體多重運移機制，建立耦合了無機質(zhì)孔隙和有機質(zhì)孔隙的頁巖儲層表觀滲透率，并運用實空間點匯函數(shù)和疊加原理，建立了非均勻分布的復雜形狀裂縫水平井產(chǎn)能預測模型，該模型還考慮了裂縫的非均勻分布，可以描述裂縫的復雜形狀，更符合現(xiàn)場實際情況。

2）考慮微觀滲流的頁巖氣井累產(chǎn)氣量約為不考慮微觀滲流頁巖氣井累產(chǎn)氣量的1.05 倍，因此，在對頁巖氣井進行產(chǎn)能評價時需要建立全面考慮頁巖氣微觀滲流的頁巖表觀滲透率模型。

3）地應力和壓裂施工等因素的存在使得壓裂裂縫不會呈現(xiàn)理想矩形形狀，常見的復雜形狀裂縫有橢形裂縫、楔形裂縫和星形裂縫，它們的產(chǎn)能依次降低，考慮壓裂裂縫為矩形裂縫會高估水平井產(chǎn)能。

4）裂縫的非均勻分布會影響水平井的產(chǎn)能，橫向非均勻中外密內(nèi)疏型得到的產(chǎn)量最多，縱向非均勻分布模式中交錯型得到的產(chǎn)量最多。因此，在壓裂施工中，采用外密內(nèi)疏?交錯型的布局有利于提高水平井的產(chǎn)量。