非常規(guī)儲層孔隙網(wǎng)絡(luò)兩相流動模擬研究

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

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

油氣藏多相流動模擬的關(guān)鍵在于表征其毛管力與相對滲透率曲線，傳統(tǒng)油氣藏主要依靠室內(nèi)巖心實(shí)驗(yàn)測定獲取，一般耗時較長，尤其對于非常規(guī)儲層，如致密、頁巖儲層，其實(shí)驗(yàn)往往難以開展。基于數(shù)字巖心孔隙網(wǎng)絡(luò)流動模擬技術(shù)，在孔隙尺度上高效模擬非常規(guī)儲層內(nèi)油、氣、水的多相流動過程，并計(jì)算毛管力與相對滲透率曲線，是行之有效的研究方法之一。

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

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

1 孔隙網(wǎng)絡(luò)兩相流動的模擬方法

1.1 基本幾何關(guān)系及排驅(qū)過程模擬

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

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

A=2r2(1+cotγ)

(1)

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

(2)

(3)

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

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

(4)

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

(5)

1.2 水濕油藏水驅(qū)油過程模擬

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

(6)

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

(7)

式(7)中：

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

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

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

1.3 油濕及混合潤濕油藏水驅(qū)油過程模擬

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

(8)

(9)

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

2 宏觀物性參數(shù)計(jì)算方法

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

(10)

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

對于孔隙i與孔隙j之間的傳導(dǎo)系數(shù)可按照調(diào)和平均獲得，可表示為

(11)

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

孔隙尺度模擬均假設(shè)流體不可壓縮，因此對于任意一個孔隙與周圍所連接孔隙之間的流入和流出量之和為零，以此可建立以孔隙壓力為未知量的線性方程組?；趫D論的處理方式，利用孔隙喉道之間的連通關(guān)系，建立鄰接矩陣，鄰接矩陣的系數(shù)即為傳導(dǎo)系數(shù)。添加鄰接矩陣的對角線為每一行矩陣?yán)奂雍偷呢?fù)數(shù)，并且根據(jù)邊界孔隙編號，將邊界孔隙對應(yīng)的矩陣進(jìn)行處理，設(shè)定入口和出口壓力值。將所建立的矩陣變?yōu)橄∈杈仃?，利用大?guī)模稀疏矩陣的求解算法，獲得每一個孔隙的壓力。繼而根據(jù)壓力場和傳導(dǎo)系數(shù)，獲取出口端面的流量q，利用達(dá)西定律獲得絕對滲透率，即

(12)

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

對于兩相流過程，對于每個時間步，根據(jù)飽和度和截面內(nèi)油水的分布形式，確定計(jì)算每個孔隙和喉道的油相和水相的傳導(dǎo)系數(shù)[20]

(13)

To=T(1-Sw)2

(14)

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

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

而對應(yīng)的整體含水飽和度可通過孔和喉的體積加權(quán)獲取。

3 實(shí)例應(yīng)用

3.1 Berea砂巖巖心

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

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

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

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

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

3.2 吉木薩爾頁巖

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

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

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

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

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

4 結(jié)論

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

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

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