陳小凡 唐 潮 杜志敏 湯連東 魏嘉寶 馬 旭

1.“油氣藏地質及開發(fā)工程”國家重點實驗室·西南石油大學 2.中國石化勝利油田分公司河口采油廠3.中國石油長慶油田分公司第五采氣廠

0 引言

合理考慮不同類型介質中頁巖氣的流動機理對于頁巖氣數(shù)值模擬具有重要意義[1]?，F(xiàn)有的頁巖氣數(shù)值模擬模型主要包括雙重介質、多重離散介質和等效介質模型，其中雙重介質模型使用較為廣泛。Kucuk和Sawyer[2]首次研究了基于雙孔模型的頁巖氣儲層壓力變化。隨后，Bumb和McKee[3]通過對Langmuir等溫吸附方程添加額外的吸附系數(shù)來研究吸附—解吸行為對不穩(wěn)定流動瞬態(tài)行為的影響。然而，以上研究工作都忽略了納米及微米尺度的擴散過程。Carlson和Mercer[4]通過在雙孔模型中引入解吸擴散項研究頁巖氣直井壓力的變化，由于該模型沒有考慮滑脫效應，在生產(chǎn)初期對頁巖氣產(chǎn)能預測結果偏小。Swami等[5]建立了一個考慮Knudsen擴散、滑脫以及吸附—解吸過程的雙重介質模型，并通過實驗室數(shù)據(jù)進行了驗證。

雖然雙重介質模型在各類商業(yè)軟件中廣泛使用，由于其固有缺陷，即使進行全網(wǎng)格或局部加密，仍然對多尺度以及大面積分布的縫網(wǎng)系統(tǒng)適應性較差。首先，傳統(tǒng)雙重介質模型均沒有考慮頁巖中多尺度流動機理，不能真實反映頁巖氣井生產(chǎn)實際[6-8]。其次，由于頁巖儲層中存在復雜的多尺度縫網(wǎng)，造成數(shù)值模擬工作中滲透率場非均質性嚴重，如果要對局部裂縫區(qū)域進行網(wǎng)格剖分，勢必產(chǎn)生大量的極小化網(wǎng)格，導致計算收斂性差甚至無法收斂[9-10]。同時，基質—裂縫界面處流度的劇烈變化也產(chǎn)生了階躍現(xiàn)象，造成的非物理振蕩會對計算精度產(chǎn)生影響[11-12]。因此，雙重介質模型應用在頁巖氣模擬中最大的限制在于難以選擇合適的特征體積單元，由于裂縫的多尺度性，裂縫性儲層是否存在特征體積單元或特征體積單元大小的確定仍存在較大爭議。因此，部分學者提出采用多重介質模型來對頁巖氣井產(chǎn)能進行模擬研究[13-14]?；诙嘀亟橘|的概念，Schepers[15]、Dehghanpour等[16]分別建立了將吸附—解吸過程與基質中流動耦合的達西流動模型。Wu等[17]和張烈輝等[18]將裂縫劃分為天然微裂縫和人工壓裂縫，建立了考慮裂縫應力敏感和滑脫效應的致密裂縫性油藏多重介質模型，氣體在基質中的流動考慮滑脫效應，在裂縫中的流動為高速非達西流動，并比較了多重介質模型與雙重介質模型的差別。Aboaba和Cheng[19]利用線性流模型描述了頁巖氣儲層中壓裂水平井的典型產(chǎn)能曲線及產(chǎn)量變化規(guī)律，但沒有考慮吸附和擴散的影響。Wang等[20]采用嵌入式離散裂縫描述大尺度裂縫、連續(xù)模型描述中小尺度裂縫的方法建立了一種模擬頁巖多尺度裂縫中的流體流動模型。方文超等[21]考慮致密儲集層的可壓縮性和基質流體的非線性滲流，建立了基于二維體積壓裂的多尺度滲流裂縫模型。

在前人研究的基礎上，筆者將基質—天然微裂縫視為雙重介質、人工壓裂縫視為離散介質，考慮基質孔隙中頁巖氣的吸附—解吸效應，建立了描述頁巖氣復雜流動過程的數(shù)學模型，基于有限體積方法離散得到頁巖氣三維滲流數(shù)值模型，利用降維方法將壓裂縫處理成二維平面實體，井筒處理為一維線元實體，再利用順序求解方法進行求解，進而模擬頁巖氣多段壓裂水平井生產(chǎn)動態(tài)和儲層壓力分布變化。所取得的研究成果可用于頁巖氣儲層體積壓裂設計以及多段壓裂水平井的生產(chǎn)動態(tài)預測。

1 數(shù)學模型

1.1 模型假設

圖1為頁巖儲層中一口多段壓裂水平井網(wǎng)格模型示意圖，整個模擬區(qū)域沿x軸方向長度為1 200 m、y軸方向寬度為800 m、z軸方向高度為100 m。水平井段長度為700 m，多段壓裂水平井空間配置如圖1-a所示，人工壓裂縫主縫半長100 m，縫高60 m；Ⅰ級分支裂縫縫長70 m，縫高60 ；Ⅱ級分支裂縫縫長50 m，縫高60 m。人工壓裂縫用二維面元實體表示，水平井段用一維線元實體表示。為簡化模型，作出如下假設：①氣藏為矩形封閉頁巖氣藏，流動為單相氣體等溫滲流，忽略重力、氣體滑脫、Knudsen擴散及應力敏感影響；②水平井位于氣藏中心位置，人工壓裂主裂縫垂直于井筒且關于井筒對稱，氣體只能通過人工壓裂縫流入井筒；③氣體在人工壓裂縫、天然微裂縫中的滲流服從達西定律，基質中頁巖氣吸附—解吸過程采用Langmuir等溫吸附方程描述；④氣井定壓生產(chǎn)，忽略水平井水平段壓降。

1.2 控制方程

1.2.1 控制方程推導

根據(jù)真實氣體狀態(tài)方程，頁巖氣密度計算式為：

式中ρgi（i=m、f、hf）表示基質、天然微裂縫和人工壓裂縫中氣體密度，kg/m3；Mg表示氣體摩爾質量，g/mol；pi（i=m、f、hf）表示基質、天然微裂縫和人工壓裂縫壓力，MPa；Z表示氣體偏差因子，無量綱；

圖1 網(wǎng)格模型示意圖

注：點C、F為兩個相鄰控制體的形心；f為界面（綠色背景）；b為界面f的邊界（紅色線）；dhf為人工裂縫無因次開度；dCf為控制體形心指向界面形心的位置向量，m；Sf為界面f上指向外側的法向向量，m2R表示氣體普適常數(shù)，8 314 J/(kmol·K)；T表示氣藏溫度，K。

依據(jù)物質守恒定律，基質中頁巖氣連續(xù)性方程為：

式中φm表示基質孔隙度，無因次；t表示時間，s；Km表示基質滲透率，mD；μg表示氣體黏度，mPa·s；qm表示基質解吸氣率，kg/（m3·s）；qmf表示基質向微裂縫竄流率，kg/（m3·s）。

式（2）左邊第1項為單位體積基質微元內流體質量變化，第2項為通過微元表面流出的流體質量通量，第3項為基質解吸氣量，第4項為基質向裂縫竄流量。在此假設微孔隙內為單相氣，由頁巖基質表面的吸附氣和孔隙內游離氣組成，根據(jù)Langmuir等溫吸附方程得到基質中解吸氣量計算式為：

式中Vm表示單位體積基質吸附氣量，m3；Vstd表示標準狀況下氣體摩爾體積，m3/mol；VL表示Langmuir體積，m3/kg；pL表示Langmuir壓力，MPa。

而頁巖氣在微裂縫中的連續(xù)性方程為：式中φf表示天然微裂縫孔隙度，無因次；Kf表示天然微裂縫滲透率，mD；qhf表示微裂縫向人工壓裂縫竄流率，kg/（m3·s）。

同理，在水力壓裂縫中的連續(xù)性方程為：

式中φhf表示人工裂縫孔隙度，無因次；Khf表示人工壓裂縫滲透率，mD；qwell表示水平井產(chǎn)氣率，kg/（m3·s）。

1.2.2 初始條件

式中pinit表示原始地層壓力，MPa。

1.2.3 邊界條件

令Γout表示求解域外邊界，Γin表示內邊界。假設模型外邊界為封閉邊界，生產(chǎn)井定壓生產(chǎn)，則邊界條件為：

2 模型離散與求解

2.1 求解域離散

由于四面體網(wǎng)格幾何中心與形心重合，為簡化計算步驟，本文以四面體網(wǎng)格和Delaunay三角形網(wǎng)格分別對基質—微裂縫和人工壓裂縫進行網(wǎng)格剖分。其他類型網(wǎng)格需要重新計算控制體形心，在此不再詳述。首先將基質—微裂縫系統(tǒng)用四面體網(wǎng)格表示，考慮為連續(xù)介質（圖1-b）；人工壓裂縫降維處理成二維平面，由四面體網(wǎng)格間Delaunay三角形界面表示（圖1-c），對裂縫降維是提高多尺度模擬計算收斂性的關鍵處理方法，若在網(wǎng)格系統(tǒng)內將裂縫考慮為三維，必然會在開度很小的裂縫空間內進行四面體網(wǎng)格剖分，進而形成大量的極小化網(wǎng)格，導致后續(xù)求解過程無法繼續(xù)[22]。Juanes等[23]研究表明，將裂縫考慮成二維比將其考慮成三維的計算收斂性顯著提高?？刂企w單元由網(wǎng)格剖分所得到的四面體確定（圖1-c），該方法可保證控制體單元互不重疊地覆蓋整個研究區(qū)域。然后，將整個求解區(qū)域Γd劃分為兩個子域Γm－f和Γhf，分別代表基質—微裂縫區(qū)域和人工壓裂縫區(qū)域，根據(jù)連續(xù)介質場理論，整個模擬區(qū)域流動控制方程（Flow Governing Equation，簡稱為FGE）的積分形式為：

2.2 空間離散

本節(jié)以基質系統(tǒng)流動方程為例說明在四面體網(wǎng)格中的求解方法。微裂縫與人工壓裂縫求解方法與之類似。首先將式（1）、（3）代入式（2），并在控制體VC范圍內進行積分，即

式中VC表示以C為形心的控制體；dV表示控制體微元體積，m3。

式（9）說明，對于任意控制體，在某一時間段內基質系統(tǒng)氣體質量流量的變化等于流入流出基質氣量與基質解吸氣量以及基質—微裂縫竄流量之和，因此保證了該模型在任意局部網(wǎng)格仍然遵守質量守恒定律。

式（9）中不穩(wěn)定流動項可以近似表達為：

根據(jù)散度定理，將體積積分轉化為面積分，則式（9）中對流項可以寫為：

進一步對控制體VC的面積分寫作對控制體各界面上通量面積分之和的形式為：

式中SC表示構成控制體VC的所有界面，以圖1-c為例，控制體VC由4個三角形界面構成。

進一步由梯形積分公式可知，在界面f上的面積分可以寫為：

式中Sf表示界面f上指向外側的法向向量，有

因此，式（9）中對流項可以寫為：

式中λm表示基質中氣體流度，m2/（mPa·s）。同理，得到源匯項有限體積計算式為：

式（15）等號右邊第1項為控制體內解吸氣量，第2項為基質—微裂縫竄流量。由于解吸氣量是時間的函數(shù)，則根據(jù)式（3）可以將基質系統(tǒng)流動方程的有限體積離散格式改寫為：

同理可以得到微裂縫及人工壓裂縫中流動方程離散格式為：

式中λf和λhf表示天然微裂縫和人工壓裂縫中氣體流度，m2/（mPa·s）； SCF表示控制體VC與大尺度裂縫共用的界面；b表示構成界面f的所有邊界；lb表示垂直邊界b的外法向向量，其中 l是邊界b的長度（二維界面厚度為1）。

2.3 順序求解方法

所謂順序求解方法，就是在每個時間步上，首先求解出某個變量，再代入其他變量表達式進行迭代求解。此方法保證了在每個時間步上計算量小于整體求解方法。假設當前時間步為k，則所有與人工壓裂縫壓力及基質壓力有關的變量都采用（k+1）時間步的值進行隱式求解，則式（16）、（18）可以寫為：

注意到式（19）、（20）中pf使用的是k時間步的值，是一個已知數(shù)。因此兩式中各含有一個未知變量phf和pm。因此可以使用Newton-Raphson方法迭代求解。將求解得到的帶入微裂縫流動方程對裂縫壓力進行顯式求解，得

2.4 壓力梯度的近似計算

顯然，為求解式（19）、（20）、（21），需要計算控制體內壓力梯度的變化，在此給出一種基于格林高斯定理的方法，該方法相對簡單且適合各種幾何結構的網(wǎng)格（結構/非結構，正交/非正交）[24-25]。體積為VC、形心為C的控制體內平均壓力梯度計算式為：

根據(jù)散度定理將體積積分轉化為面積分，有

對于離散界面，式（23）可以寫為：

然后，將沿控制體表面的面積分用積分中值定理近似表示為表面形心處的插值乘以表面向量的形式，即

由式（25）可知，為計算控制體內壓力梯度，首先需要知道控制體表面向量Sf，并通過控制體形心處壓力值pC和相鄰控制體形心處壓力值pF插值得到界面f上的壓力，其插值格式為：

其中gC和gF為考慮控制體體積VC和VF的權重值，計算公式為：

考慮到相鄰的兩個控制體（網(wǎng)格編號分別為j和k），表面向量Sf不能同時向外，因此對于特定網(wǎng)格定義表面向量方向由控制體網(wǎng)格編號確定，始終由編號較小的控制體指向編號較大的控制體，則式（25）改寫為：

2.5 非正交網(wǎng)格中表面向量的處理

由于本文中采用的非結構網(wǎng)格，其網(wǎng)格形態(tài)往往是非正交的。即控制體表面向量與連接控制體形心的向量不共線，也就是說控制體表面法向向量與流體速度（壓力梯度）方向不一致。這種情況下，因為存在一個垂直于向量的分量，控制體的壓力梯度不能寫成f(pC, pF)形式[26-29]。

定義e為沿著形心C、F連線所確定方向的單位矢量，m，則沿著e方向的壓力梯度可以表示為：

式中rC、rF表示形心C和F的位置向量，m；dCF表示相鄰控制體形心C與F間的距離，m。

因此，為了實現(xiàn)非正交網(wǎng)格中通量的線性化，表面向量Sf應該寫作兩個向量Ef、Tf之和，即

其中向量Ef沿著形心C、F連線方向，m2；Tf為Sf的非正交分量，m2，則通過界面f的通量可以表示為：

對于向量Ef、Tf的計算通常有以下3種方法（表1）。

上述3種分解方法所得到的結果均滿足式（7）中擴散通量計算需要，其主要差別主要體現(xiàn)在非正交網(wǎng)格計算過程中的準確性和穩(wěn)定性。根據(jù)前人研究發(fā)現(xiàn)，即使在高度非正交的網(wǎng)格中超松弛修正方法仍能保證通量計算的穩(wěn)定性[24-25,30-31]。

表1 非正交網(wǎng)格中控制體表面向量分解方法

2.6 模型驗證

2.6.1 模型對比驗證

為了驗證有限體積方法的正確性，建立矩形封閉氣藏模型，在基質和微裂縫中僅考慮黏性流的影響，忽略吸附—解吸機制，并認為人工壓裂縫為無限導流縫，模型參數(shù)如表2所示。將模擬結果與油氣藏商業(yè)數(shù)值模擬軟件Eclipse（以下簡稱軟件Eclipse）計算的結果進行對比，如圖2所示，本文模型與軟件Eclipse計算的水平井產(chǎn)氣量基本一致，說明本文采用的數(shù)值計算方法正確、可行。

如圖2所示，在早期段本文數(shù)值解與商業(yè)軟件計算結果存在一定誤差，原因在于壓裂縫附近網(wǎng)格精度不夠導致（最大單元步長15 m，最小單元步長3 m，最大單元遞增率1.2，單元曲率因子0.5，狹窄區(qū)域分辨率0.8），可以通過對壓裂縫位置處的網(wǎng)格進行加密處理以提高計算精度。

2.6.2 求解方法對比

基于前述模型，通過不同求解方法（順序求解方法和全隱式求解方法）進行計算，結果顯示差距主要在生產(chǎn)初期，順序求解方法計算得到的開井瞬時產(chǎn)氣量為195×104m3/d，而全隱式求解方法計算得到的開井瞬時產(chǎn)氣量為167×104m3/d，但隨著計算推進，兩種方法計算的產(chǎn)氣量迅速趨于一致（圖3），進一步驗證了模型的正確性。

表2 多段壓裂水平井數(shù)值模擬參數(shù)表[32]

圖2 本文數(shù)值解與軟件Eclipse計算結果對比圖

3 實例模擬

3.1 預設參數(shù)

根據(jù)本文參考文獻[1]，建立頁巖氣多段壓裂水平井模型，原始地層壓力30 MPa，地層溫度為343.15 K，基質孔隙度為8%，基質滲透率為0.003 mD，天然微裂縫孔隙度為0.3%，天然微裂縫滲透 率 為 0.5 mD，Langmuir體 積 為 4×10－3m3/kg，Langmuir壓力為5 MPa，頁巖密度為2 600 kg/m3，甲烷分子摩爾質量為16 g/mol，標況下頁巖氣摩爾體積為 0.022 4 m3/mol，氣體黏度為 1.85×10－2mPa·s，井底流壓為5 MPa。通過軟件Eclipse和CMG計算發(fā)現(xiàn)，在上述模型壓力變化過程中氣體偏差因子變化范圍介于0.908 9～0.982 8，利用多項式擬合將偏差因子改寫成壓力函數(shù)的形式。多段壓裂水平井空間配置如圖1-a所示，人工壓裂縫主縫導流能力為15 D·cm，Ⅰ級分支裂縫導流能力為10 D·cm，Ⅱ級分支裂縫導流能力為7 D·cm。

3.2 模擬結果分析

3.2.1 Langmuir體積的影響

基于物理模型參數(shù)，模擬頁巖氣多段壓裂水平井生產(chǎn)情況，如圖4所示，在前5年儲層壓力下降區(qū)域主要集中在體積改造區(qū)域，而未壓裂區(qū)域內儲層壓力變化幅度較小，表明該階段產(chǎn)出的氣體主要來自體積改造區(qū)域內的游離氣及解吸氣。

圖5為不同Langmuir體積所對應的平均地層壓力、解吸氣量、產(chǎn)氣量以及解吸氣貢獻比變化曲線。可以看出，解吸氣使地層壓力上升，但作用有限，對產(chǎn)氣量影響不大；隨著生產(chǎn)時間延長，解吸氣量絕對值逐漸減少，但是在產(chǎn)氣量中所占的比例逐漸上升。

圖3 不同求解方法計算結果對比圖

圖4 不同生產(chǎn)時間地層壓力分布圖

圖5 不同Langmuir體積下水平井生產(chǎn)動態(tài)預測曲線圖

3.2.2 裂縫形態(tài)的影響

通過改變壓裂段數(shù)進一步研究裂縫形態(tài)對頁巖氣水平井產(chǎn)氣量的影響，模擬結果顯示，生產(chǎn)5年時，3段壓裂模式下的壓降區(qū)域范圍與4段壓裂模式下的壓降區(qū)域范圍較接近（圖6）；4段壓裂模式下投產(chǎn)初期產(chǎn)氣量（64×104m3/d）為3段壓裂模式下初期產(chǎn)氣量（47×104m3/d）的1.37倍，但4段壓裂模式下產(chǎn)氣量遞減速度更快，5年后兩種壓裂模式下的水平井產(chǎn)氣量基本一致，且該階段的累產(chǎn)氣量曲線平行（圖7）。由此可見，壓裂段數(shù)繼續(xù)增加對于增加累產(chǎn)氣量效果已不顯著，需要通過增加裂縫半長來擴大儲層改造體積。確定合理的壓裂段數(shù)，同時獲得較長的壓裂縫長，是頁巖氣水平井增產(chǎn)改造措施的核心內容。

圖6 不同壓裂段數(shù)水平井生產(chǎn)5年時地層壓力分布圖

圖7 不同壓裂段數(shù)下頁巖氣水平井生產(chǎn)曲線圖

4 結論

1）采用所建立的模型與軟件Eclipse計算的多段壓裂水平井產(chǎn)氣量基本一致，說明該數(shù)值計算方法正確、可行；通過順序求解方法和全隱式求解方法的對比，進一步驗證模型的正確性。

2）解吸氣對地層壓力有補充作用，但作用有限，對產(chǎn)氣量影響不大，隨著生產(chǎn)時間的延長解吸氣量在產(chǎn)氣量中所占比例逐漸上升。

3）確定合理的壓裂段數(shù)，同時獲得較長的壓裂縫長，是頁巖氣水平井增產(chǎn)改造措施的核心內容。