非常規儲層孔隙網絡兩相流動模擬研究

王秀坤,劉海成,吳忠維,崔傳智

(1.中國石油大學(北京)非常規油氣科學技術研究院，北京 102249；2.中國石化勝利油田分公司勘探開發研究院，東營 257015；3.中國石油大學(華東)石油工程學院，青島 266580)

油氣藏多相流動模擬的關鍵在于表征其毛管力與相對滲透率曲線，傳統油氣藏主要依靠室內巖心實驗測定獲取，一般耗時較長，尤其對于非常規儲層，如致密、頁巖儲層，其實驗往往難以開展?；跀底謳r心孔隙網絡流動模擬技術，在孔隙尺度上高效模擬非常規儲層內油、氣、水的多相流動過程，并計算毛管力與相對滲透率曲線，是行之有效的研究方法之一。

借助微納米CT(computed tomography)、FIB-SEM(focused ion beam scanning electron microscope)等實驗手段，可獲得巖石微觀孔隙空間的三維數據體，二值化分割后即可獲得數字巖心[1]?？紫冻叨攘鲃幽M方法分為直接模擬和孔隙網絡模擬，前者直接在數字巖心的孔隙空間的體素點上開展數值模擬計算[2-3]，計算成本高。而后者則一般利用最大球體法[4]、分水嶺算法[5]等算法抽提獲得表征其孔喉幾何形狀和拓撲空間的孔隙網絡，然后據此開展多相流動模擬計算[6]。通常采用的準靜態孔隙網絡流動模擬方法只考慮毛細管力的作用，可快速模擬排驅和吮吸過程中油水的三維空間展布變化，繼而獲取宏觀的流動物性參數，如絕對滲透率、相對滲透率曲線和毛細管力曲線等。起初，孔隙網絡兩相流動模擬法主要是在挪威國家石油公司研究中心的Oren等[7]和英國帝國理工大學Blunt[8]研究基礎之上而創建的。Valvatne等[9]于2004年開發了業界公認的一套的兩相流動模擬軟件，實現了對初次排驅和二次吮吸過程的模擬，可得到相對滲透率、毛細管力曲線和電阻率變化規律等宏觀參數，極大推動了孔隙網絡流動模擬的研究工作。但 van Dijke等[10]指出，Blunt軟件中對油膜坍塌機制的考慮不合理，并利用能量守恒提出了更精確的油膜存在條件模型。Helland等[11]指出了Blunt軟件中斜三角形截面傳導系數的誤差，并提出了更貼近孔喉特征的四角星形截面模型。Ryazanov[12]據此四星形模型，建立了常規儲層的孔隙網絡兩相流動模擬方法，而對非常規儲層沒有涉及。Petrovskyy等[13]完善了Ryazanov[12]孔隙網絡多相流動的模擬計算速度，提出了更高效的非潤濕相圈閉檢測算法。近年來，中國學者也對孔隙網絡流動模擬方法開展了與國際接軌的相關性研究，取得了一定的特色成果[14-16]。

目前，基于孔隙網絡模型的非常規儲層兩相流動模擬研究較少。為此，基于四角星形截面的處理方式，結合高性能圖論算法和開源可視化表征軟件Paraview[17]，在考慮活塞式驅替、卡斷填充、協同式填充、帶油膜充填等孔喉微觀填充機理基礎上，建立了孔隙網絡兩相流動模擬方法；通過與實驗結果對比驗證了該方法的有效性和準確性；并運用該孔隙網絡兩相流動模擬方法對吉木薩爾頁巖儲層的復雜油水兩相流動機理進行了探討。

1 孔隙網絡兩相流動的模擬方法

1.1 基本幾何關系及排驅過程模擬

當孔隙網絡模型中每個孔和喉的截面采用四角星形時(圖1)，通過幾何關系可得孔或喉的截面積A，其計算公式為

圖1 四角星形截面幾何關系示意圖Fig.1 Schematic of the four-pointed crossing sectional area

A=2r2(1+cotγ)

(1)

式(1)中：A為孔或喉的截面積，m2；r為孔或喉內切圓半徑，m；γ為四角星的半角，(°)。

(2)

(3)

式中：G為截面的形狀因子，無量綱；P為周長，m。

對于初次排驅過程，即油藏成藏過程，儲層首先被水充滿，烴源巖生油后通過生烴壓力的作用進入儲層，油相克服毛細管力驅替水相，最終形成油藏。排驅和后續吮吸過程，均認為巖石為水濕。對于這個過程的孔隙網絡模擬，毛細管力逐步增大，喉道起到主要的約束作用，任一喉道入口毛細管力大小可表示為式(4)[12]。當毛細管力超過式(4)時，油相進入此喉道以及與之連接的孔隙。

(4)

每個孔或喉的含水飽和度Sw可通過截面積占比獲得，可表示為

(5)

1.2 水濕油藏水驅油過程模擬

對于水濕油藏水驅油過程，即二次吮吸過程，毛細管力逐步降低，水相逐步填充孔和喉。喉道主要有活塞式填充和卡斷填充兩種微觀填充機理，活塞式喉道填充發生時需要其所連接的孔隙有一個已被水相填充，否則只能發生卡斷填充。喉道卡斷填充的臨界毛細管力為[12]

(6)

需要說明的是，當孔隙被水相填充后，與之相連接的喉道也會隨之被填充，孔隙制約著喉道活塞式充填的發生，這與排驅過程中喉道制約著孔隙被油相侵入的道理是相對應的。孔隙一般為協同式填充，其臨界毛細管力為[9]

(7)

式(7)中：

對于排驅和吮吸的過程，需要考慮潤濕滯后現象[18]。排驅過程中采用后退角，而吮吸過程采用前進角，且前進角大于后退角。初次排驅結束之后，標志著油藏的形成和最大的毛細管力。二次吮吸過程發生時，毛細管力逐步降低，此時角隅的水相區域膨脹，接觸角從后退角逐步增大，直至變為前進角，含水飽和度緩慢增加，但僅以角隅和水膜的形式存在，對應的含水飽和度計算公式見式(5)。一般在模擬計算時會直接給定前進角和后退角的值，水濕巖心，典型的前進角為50°～60°，后退角0°～10°。

式(6)、式(7)就構成了孔隙網絡準靜態兩相流動模擬的主要微觀填充機理，初次排驅較簡單，但對于二次吮吸過程，多個微觀機理相互之間是競爭關系，需要按照公式中所計算的各類填充臨界毛細管力的大小進行排序，較大的臨界毛管力對應的填充機理優先發生。但由于孔隙喉道輸入巨大，每一步填充一個孔或喉是不現實的，為獲取更高的效率，每一個毛細管力變化內，有很多孔和喉均有填充發生，但填充發生的早晚影響著協同式孔隙填充發生過程中n的值，為此提出：首先，依據式(6)確定所有卡斷填充發生的喉道并令其填充，繼而，依據式(7)判斷的協同式填充是否發生，但對于n=0時，也需要額外按照n=1判斷孔隙內油是否為死油，最后，與已被水填充的孔隙的喉道也令其均被水填充。

此外，一個重要的問題是，油相作為非潤濕相會以油團的形式被圈閉在孔隙網絡中，形成殘余油，油團可能是一個孔[式(7)，n為0時]，也可能是若干個孔隙喉道的連通體，此過程需要高效算法實時監測不可動油團的形成[13]。通過采用先進的圖論算法，借助NetworkX[19]高效實現了對孔隙網絡連通關系的動態監測。其主要思想為：①增加 Outlet 節點與出口處所有孔節點連接；②找到與Outlet連接的所有含油的孔節點，即為可動油；③所有含油孔節點去掉可動油孔節點即為不可動的死油團。

1.3 油濕及混合潤濕油藏水驅油過程模擬

當儲層為油濕時，水驅油過程中，毛細管力為阻力，此過程類似于初次排驅過程，區別在于存在水膜，主要的微觀充填機理有帶油膜活塞式充填和不帶油膜活塞式充填。

(8)

(9)

此外，對于混合潤濕巖心水驅油過程的模擬。認為孔隙半徑較大的孔隙為油濕孔，與之連接的喉道也為油濕，給定油濕孔隙體積所占比，標記油濕孔隙和喉道。對于水濕的孔隙和喉道所組成的孔隙網絡按照2.2節開展模擬計算，此時保持油濕孔隙網絡均為存在水膜的油相充填，并按照整體的孔隙網絡計算相對滲透率曲線和毛細管力，繼而按照上述算法對油濕的孔隙網絡開展模擬計算，最終獲得整體的相對滲透率曲線和毛細管力曲線。

2 宏觀物性參數計算方法

首先開展單相流的模擬計算，獲取巖心的絕對滲透率。具體地，計算每個孔隙和喉道的傳導系數T，計算公式為

(10)

式(10)中：T為傳導系數，m3；L為巖芯長度，m。

對于孔隙i與孔隙j之間的傳導系數可按照調和平均獲得，可表示為

(11)

式(11)中：Ti為孔隙i的傳導系數，m3；Tj為孔隙j的傳導系數，m3；Tij為喉道的傳導系數，m3。

孔隙尺度模擬均假設流體不可壓縮，因此對于任意一個孔隙與周圍所連接孔隙之間的流入和流出量之和為零，以此可建立以孔隙壓力為未知量的線性方程組?；趫D論的處理方式，利用孔隙喉道之間的連通關系，建立鄰接矩陣，鄰接矩陣的系數即為傳導系數。添加鄰接矩陣的對角線為每一行矩陣累加和的負數，并且根據邊界孔隙編號，將邊界孔隙對應的矩陣進行處理，設定入口和出口壓力值。將所建立的矩陣變為稀疏矩陣，利用大規模稀疏矩陣的求解算法，獲得每一個孔隙的壓力。繼而根據壓力場和傳導系數，獲取出口端面的流量q，利用達西定律獲得絕對滲透率，即

(12)

式(12)中：q為流量，m3/s；μ為液體黏度，Pa·s；pin、pout分別為入口和出口處的壓力，Pa。

對于兩相流過程，對于每個時間步，根據飽和度和截面內油水的分布形式，確定計算每個孔隙和喉道的油相和水相的傳導系數[20]

(13)

To=T(1-Sw)2

(14)

式(13)中：Tw為水相的傳導系數，m3；To為油相的傳導系數，m3；C(G,θ)為潤濕相額外附加流動阻力的修正系數，具體值參考文獻[21]。

最后按照與單相流計算相同的方式，計算油相和水相的壓力場，獲取油相滲透率ko、水相滲透率kw，即可獲得油、水相對滲透率

而對應的整體含水飽和度可通過孔和喉的體積加權獲取。

3 實例應用

3.1 Berea砂巖巖心

kr為相對滲透率圖2 油驅水(初次排驅)過程相對滲透率曲線Fig.2 Relative permeabilities of oil injection (primary drainage) process

讀入已獲得的經典Berea砂巖的孔隙網絡StateOil格式數據[9]，按照上述計算方法，對于油驅水(成藏)、水濕巖心水驅油以及油濕巖心水驅油過程進行了模擬。本代碼效率較高，普通筆記本電腦模擬計算時間少于5 min，得到圖2、圖3所示的模擬計算結果，與室內實驗測量結果[22]吻合率高，證實了所建立孔隙網絡流動模擬方法的有效性和準確性。由圖3可知，水濕巖心相對滲透率曲線殘余油飽和度高，且水相相對滲透率端點值小于0.2；而油濕巖心相對滲透率曲線殘余油飽和度低，且對應的水相端點值很高，接近0.8。這些相對滲透率曲線特征是與前人的研究是一致的[23]。

圖3 巖心水驅油過程的相對滲透率曲線Fig.3 Relative permeabilities of water flooding process

在此研究基礎之上，對接開源軟件Paraview實現了兩相流動的三維可視化表征。圖4為水濕和油濕情況下水驅剩余油分布，其中孔隙以球體的形式表征，而喉道以棍棒的形式表征，其顏色表征著流體的飽和度，紅色為油相，藍色為水相。圖4可證實水濕巖心剩余油主要以孤立油團的形式存在，且主要分布于孔喉比較大的孔隙內，而油濕巖心水驅剩余油主要分布于小孔喉的區域，而大孔動用程度高，與水濕情況剛好相反。

紅色為油相；藍色為水相圖4 水驅后剩余油分布的三維可視化表征Fig.4 Three dimensional visualization of water-flooded rock

3.2 吉木薩爾頁巖

在Berea砂巖兩相流動模擬結果與實驗驗證的基礎上，利用FIB-SEM實驗實現了對吉木薩爾頁巖的三維孔隙空間數字化表征。其分辨率為20 nm×20 nm×20 nm，數據空間大小為500×500×500。應用快速傅立葉變換算法，去除了三維圖像中的噪聲偽影，繼而通過設置閾值灰度值對孔隙空間進行圖像分割，并采用最大球體算法提取頁巖的三維孔隙網絡。具體的頁巖三維孔隙空間和孔隙網絡的Paraview可視化如圖5所示。

圖5 吉木薩爾頁巖孔隙空間的三維可視化及孔隙網絡模型可視化表征Fig.5 Three dimensional visualization of pore space and pore network of Jimsar shale

吉木薩爾頁巖孔隙網絡共有13 419個孔隙和31 393個喉道，連通孔隙度為16.35%，與實驗室內氣測孔隙度19%接近?？紫堵实膿p失是由于小于20 nm的孔和喉道無法識別造成的。孔隙和喉道的平均半徑分別為29.75 nm和19.13 nm。盡管有幾個相當大的孔和喉，但整體平均孔喉比為2.29。應用所提出的孔隙網絡方法，計算出x、y和z方向的達西滲透率為0.004～0.006 mD，這也是頁巖油藏在地層條件下的合理值。頁巖儲層通常認為是混合潤濕，此處采用1/2水濕和1/2油濕的混合潤濕條件，將較大/較小部分的孔隙和喉道設置為油濕，其余為水濕情況，模擬計算得到的油水相對滲透率曲線如圖6所示。

圖6 吉木薩爾頁巖混合潤濕條件下油水兩相相對滲透率曲線Fig.6 Oil-Water relative permeability curves of mixed-wetted Jimsar shale

由圖6可知，兩種混合潤濕性情況均表明殘余油飽和度較高(即接近40%)，且兩相流區域狹窄，這意味著頁巖油的采收率遠低于常規高滲透油藏。對于較大孔隙油濕情況，水相相對滲透率極低，且端點值小于0.01，這意味著頁巖基質中的水幾乎不會流動，這與致密儲層“滲透率圈閉”現象是一致的[24-26]。對于較小孔隙油濕情況，隨著含水飽和度的增加，水相相對滲透率增加得更快，端點值在0.3左右。兩者的差異表明，較大的孔隙具有主要的流動能力，當較大的孔隙為油濕時，水首先充滿水濕的較小孔隙，然后充滿油的較大孔隙被卡斷形成死有團。然而，當較大的孔隙為水濕時，水首先填充孔隙，然后延伸至較小的油濕孔隙。油濕孔隙的分布對油水流動機制起著重要作用，未來的研究需要建立更大尺度的孔隙網絡模型和更精細的復雜潤濕性表征，使其更能代表頁巖地層，這是未來工作的研究方向。

4 結論

基于實際數字巖心抽提的孔隙網絡模型，利用四角星形截面流體微觀賦存機理，結合高性能圖論算法和開源可視化表征軟件，考慮了活塞式驅替、卡斷填充、協同式填充等孔喉微觀填充機理，建立了孔隙網絡兩相流動模擬方法。實現了對孔隙網絡內非潤濕相連通關系的動態監測，模擬了初次排驅和水驅油過程，基于經典Berea砂巖的孔隙網絡，計算得到了不同潤濕性條件下油水的相對滲透率曲線，對比實驗結果驗證了模擬結果的有效性和準確性，可視化表征結果方便的表征了不同潤濕條件下剩余油的分布規律。得出如下主要結論。

(1)證實了水濕巖心剩余油主要以孤立油團的形式存在，且主要分布于孔喉比較大的孔隙內，而油濕巖心水驅剩余油主要分布在小孔喉的區域，而大孔動用程度高，與水濕情況剛好相反。

(2)對吉木薩爾頁巖儲層油水兩相流動進行了初步探討，結果表明頁巖兩相流區域狹窄，殘余油飽和度高，水相存在“滲透率圈閉”現象，而且明確了油濕孔隙的分布對油水流動機制起的重要作用。