基于數字巖心與格子Boltzmann方法的致密砂巖自發滲吸模擬研究

汪勇，孫業恒，梁棟，蔡建超

1 中國石化勝利油田分公司勘探開發研究院，東營 257015

2 中國石化勝利油田分公司科技管理部，東營 257015

3 中國石油大學油氣資源與探測國家重點實驗室，北京 102249

*通信作者, wangyong733.slyt@sinopec.com

0 引言

致密油氣藏廣泛發育于我國各大含油氣盆地，其資源量約為常規油氣的3倍[1]。中國致密油藏的勘探與開發起步較晚，但發展速度較快，目前已在鄂爾多斯、松遼、準噶爾、渤海灣等多個含油氣盆地發現致密油勘探開發有利區[2-3]。隨著常規石油資源量的減少，低滲—致密油藏的重要性日益突出，提高致密油藏開發技術水平[4]，對我國石油工業持續穩定發展具有重要戰略意義。

在水力壓裂和注水開發過程中，自發滲吸是低滲—致密油藏基質—裂縫間產生油水置換現象的重要機制[5-6]，基質中的油氣得以進入高滲裂縫網絡從而提高油氣產能。致密油藏中自發滲吸是一個復雜過程，受多種因素影響：巖石微觀結構、流體特性及邊界條件等[7-8]。不同條件下，影響滲吸規律的主導因素不盡相同。精確刻畫滲吸規律是提高油藏采收率的關鍵基礎問題之一。

近年來，數值模擬方法已被廣泛應用于油藏開發滲流特性的研究[9-10]，而基于數字巖心模型的孔隙尺度流動模擬為研究熱點[11-12]。格子Boltzmann方法(Lattice Boltzmann Method，LBM)[13]是一種通過流體粒子的碰撞和遷移來描述流體流動的介觀方法。相比于傳統計算流體力學，該方法易于處理復雜幾何邊界和實現并行計算，常用于微觀滲流模擬，或直接計算巖石宏觀物性參數等[14-15]。

本文基于勝利油田東營凹陷金家—樊家鼻狀構造帶致密砂巖微觀孔隙結構的數字巖心模型，應用基于LBM的多相流模型開展孔隙尺度自發滲吸模擬研究。研究聚焦于巖心微觀孔隙結構與潤濕性對自發滲吸過程中流體界面變化及采出程度的影響，有助于明確致密油藏水力壓裂開發中自發滲吸效應的特征及探明其影響因素。

1 致密砂巖數字巖心模型構建

巖心樣品來自勝利油田樊154區塊，巖性為致密砂巖，埋藏深度2506.91～3205.01 m不等。由于原始樣品具有一定的含油性，所以首先進行洗油處理。之后在洗完油的標準巖心上鉆取了直徑為3 mm的樣品，用于CT掃描構建數字巖心模型。CT掃描的圖像數據來自于VerXRM500型微米CT顯微鏡，光源電壓30～160 kV，最高分辨率 0.5 μm。

共有三塊巖心進行了微米CT掃描，掃描電壓為110 kV，掃描分辨率為1 μm。掃描結束后，掃描出的二維圖像由重構軟件重構成三維數據體。圖1a、b、c為此次研究區編號為1、2和3巖心的二維掃描原始灰度圖像。從掃描圖像可以看出，1號樣品孔隙發育較好，孔隙連通性較高，孤立孔隙較少，孔隙數量較高；2號樣品孔隙呈裂縫狀分布，同樣具有較好的孔隙連通性，孔隙數量較高。3號樣品孔隙發育較差，孔隙連通性不高，孤立孔隙較多，孔隙數量較少。為了提升計算效率與排除孤立孔隙對滲吸模擬的影響，三塊巖心的大小為1503像素的子體積模型進行一系列包括提取、濾波、二值化等處理手段從原始CT掃描數據中提取出用于自發滲吸模擬(圖1d、e、f，孔隙度依次為0.186，0.139，0.103)。從圖中可以看出，提取出的子體積模型很好的表征了三塊巖心的孔隙結構，1號巖心孔隙較多，連通性較好；2號巖心孔隙呈片狀發育，連通性較好；3號巖心連通孔隙較少，分布非均質性較強。

2 孔隙尺度自發滲吸數值模擬

此次研究中主要應用基于LBM的顏色梯度模型開展致密儲層數字巖心模型自發滲吸模擬計算。下面簡要介紹控制流體輸運的LBM基礎理論。

LBM利用離散空間網格格點上具有不同速度方向的流體粒子的碰撞和遷移來表征流體流動。一個格點某一離散速度方向的流體粒子占有量即為分布函數fi(x)，其在固定格子上的運動方程如下所示

式中，Ωi為碰撞算子，決定著演化方程表征的流動物理及刻畫流動問題的能力[13]。流體流動過程中宏觀物理量，如密度和速度可分別由下式計算得到

圖1 編號為1、2、3號巖心的二維CT掃描灰度圖像與三維子體積圖像Fig. 1 2D CT scanning grayscale images and 3D sub-volume images of cores 1, 2 and 3

LBM的介觀特征為從底層構建考慮多相流體微觀相互作用的宏觀流動模型提供了基礎理論框架[16]。目前已發展了多種多相流模型，其中顏色梯度模型(Colour Gradient Model)由于其具有界面張力可獨立調節的特性和較小的界面厚度等特點，被廣泛用于多孔介質多相流體流動研究[17]。該模型以紅藍兩色標記非混溶流體，通過在碰撞項中添加微擾動以引入表面張力，即在該模型中表面張力被視為壓力局部各向異性。碰撞后非混溶流體粒子通過重新標色驅使流體流向相同顏色的流體區域，達到相分離的目的[16,18]。

由于致密油儲層孔隙結構的復雜性及毛細流動過程中較小的毛管數(較低的流速)，為了提高基于三維數字巖心模型滲吸模擬的精度和穩定性，模擬中利用了優化擾動項的顏色梯度模型及多松弛(MRT)碰撞項[19]。優化的顏色梯度模型中包括三個時空演化方程。其中一個為全局分布函數，用于控制壓力和速度場的演化，另外兩個演化方程用于模擬兩相界面的演化。優化的顏色梯度如下[20]：

兩個獨立的演化方程用于計算紅藍流體密度場(ρr和ρb)對流，其中紅色流體的演化方程和平衡分布函數分別為

自發滲吸模擬開始前潤濕流體(藍色)位于入口端20層格子緩沖層內，非潤濕流體(紅色)飽和整個孔隙空間和出口端20層格子緩沖層。在整個模擬過程中兩相流體不受外力作用，完全依靠毛管力驅動潤濕流體滲入致密油儲層孔隙空間，因此計算域六面均采用周期邊界條件。與滲吸方向垂直的兩個側面，在遷移前需要與緩沖層內流體序參數進行匹配。固體壁面采用反彈邊界以適應復雜的孔隙空間結構。由于本次研究在孔隙尺度考察孔隙結構特征和潤濕性對自發滲吸的影響，并且忽略重力的影響，因此在模擬中潤濕和非潤濕流體的黏度均設置為0.01，密度設置為1.0(格子單位)。模擬中不同的潤濕角是通過設定固體壁面序參數而實現。在不同孔隙結構自發滲吸模擬中設置的潤濕角與不同巖心樣品的實驗測取值進行匹配(三塊巖心實驗測量的潤濕角均在30°～40°范圍內)。

3 模擬結果及分析

3.1 微觀孔隙結構對自發滲吸的影響

為了直觀對比不同微觀孔隙結構自發滲吸驅替過程中的滲吸前緣演化特征，圖2以一定時間間隔展示了自發滲吸過程中潤濕和非潤濕流體分布。結果顯示不同孔隙結構內潤濕和非潤濕流體相間界面演化規律整體相似。這是由于在數值模擬中三種孔隙結構的潤濕角設置值較為接近，且均為強潤濕性。根據孔隙尺寸的相對大小，潤濕流體的流動方式可明顯分為兩類：在大孔隙內潤濕流體優先沿孔隙壁面或角隅流動，該現象造成大孔隙內兩相界面明顯滯后。當大孔隙內的潤濕流體流經喉道處時，由于“卡斷”效應，造成非潤濕流體被捕獲在孔隙中央。尺寸相對小的孔隙內潤濕流體快速填充整個孔隙，當大孔隙周圍有很多小孔隙存在時，大小孔隙內潤濕流體填充速率的差異引起“繞流”效應，造成大量的非潤濕流體殘留在大孔隙中。1號巖心模型的孔隙結構尺寸大小和連通性相對均勻，非潤濕流體主要因“卡斷”效應而滯留；2號巖心模型孔隙結構片狀發育明顯及連通性相對較好，非潤濕流體主要因“繞流”效應而滯留；3號巖心模型孔隙結構也呈現出片狀發育，但由于其連通性差，該模型內非潤濕流體主要因“卡斷”效應而滯留。

圖3為不同巖心采出程度隨時間的變化。1號巖心孔隙結構的滲吸速率一直變化，計算終止時刻仍未達到穩定狀態，滲吸過程仍在緩慢發生。2號和3號巖心孔隙結構滲吸開始時刻滲吸速率十分迅速，達到拐點后2號巖心后期滲吸速率幾乎沒有明顯變化，而3號巖心滲吸速率減緩，滲吸過程仍在持續，直到達到下一個拐點，滲吸速率基本達到穩定。三種類型巖心孔隙結構自發滲吸的最終采出程度分別為0.60，0.70，0.51。2號巖心孔隙結構內的最終采出程度最大，這也可以從計算停止后兩相流體的賦存狀態定性判斷出(圖2)。而較大的最終采出程度與該巖心微觀孔隙結構的片狀發育及較好的連通性直接相關。

圖4更為直觀的展示了不同孔隙結構模型滲吸過程中滲吸速率隨時間的變化規律。滲吸初始時刻，3號巖心孔隙結構的滲吸速率最大，1號巖心孔隙結構的滲吸速率最小。滲吸初始階段，三種孔隙結構的自發滲吸速率均迅速減小，隨后以一定的速率持續緩慢變化。在0.25×106～1.0×106計算時間內，2號巖心的滲吸速率仍然保持較大的速率，而在此時間段內其它兩塊巖心孔隙結構的自發滲吸速率相對較小，并以較小的波動緩慢減小。

3.2 巖石潤濕性對自發滲吸的影響

為深入認識不同潤濕性對致密砂巖自發滲吸過程中孔隙內流體動態填充和整體采收率的影響，以2號巖心的微觀孔隙結構為基礎，在不同潤濕條件下開展自發滲吸模擬計算。模擬中不同潤濕角設置值分別為20.0°，37.3°，60.0°，80.0°。

圖2 編號為1、2、3號巖心子體積模型內自發滲吸過程中兩相界面演化(藍色為潤濕流體，紅色為非潤濕流體，巖石骨架不顯示)Fig. 2 Evolution of two-phase interface during spontaneous imbibition in the sub-volume model of cores 1, 2 and 3 (blue indicates the wetting fluid, red indicates the non-wetting fluid, and the rock skeleton is not shown)

圖3 三塊巖心自發滲吸采出程度隨時間變化關系Fig. 3 The recovery degree of spontaneous imbibition of three cores varies with time

圖4 三塊巖心自發滲吸速率隨時間變化關系Fig. 4 The spontaneous imbibition rates of three cores varies with time

圖5以一定時間間隔展示了不同潤濕條件下自發滲吸過程中的兩相分布。結果顯示，潤濕性對滲吸前緣的形態及兩相流體的空間分布有較大影響。強潤濕條件下(潤濕角≤37.3°)，潤濕流體主要以“角流”形式優先侵入孔隙角隅。滲吸前緣整體雜亂、分散，主終端界面滯后明顯。隨著潤濕角增大，以“角流”形式的流動明顯減小，兩相界面的形態也相應變得規則、緊湊。當潤濕角進一步增大到80°時，潤濕流體流動能力明顯減弱，僅能侵入入口端附近的部分孔隙。

圖5 2號樣品不同潤濕條件下孔隙內滲吸前緣演化(藍色-潤濕流體，紅色-非潤濕流體，巖石骨架不顯示)Fig. 5 Evolution of imbibition interface in the pores under different wetting conditions of core 2 (blue indicates the wetting fluid, red indicates the non-wetting fluid, and the rock skeleton is not shown)

圖6 滲吸初期不同潤濕條件下非潤濕流體分布(紅色-非潤濕流體，灰色-巖石骨架，潤濕流體不顯示)Fig. 6 Distribution of non-wetting fluid under different wetting conditions in the early stage of imbibition (red indicates the non-wetting fluid, gray indicates the rock skeleton, and the wetting fluid is not shown)

滲吸初期，強、弱潤濕條件下(潤濕角≤60°)模型入口端均有一定體積的非潤濕相液滴析出(圖6)，定性說明了在該潤濕條件下可以發生逆向滲吸現象。但不同潤濕條件下逆向滲吸發生程度和發生位置明顯不同。潤濕性越強(潤濕角越小)，入口端面上尺寸較小和較大孔隙均可作為逆向滲吸發生位置。隨著潤濕角增大，逆向滲吸現象僅能通過尺寸較大的孔隙發生，當潤濕性接近中性潤濕時，基本不發生逆向滲吸現象。

不同潤濕條件下滲吸采出程度隨時間的變化如圖7所示。可以看出，同一微觀孔隙結構模型的自發滲吸驅替效果顯著受潤濕強度的影響。潤濕性接近中性潤濕時，自發滲吸現象不明顯。當潤濕角減小時，滲吸驅替效率明顯提高，但潤濕角由37.3°進一步減小到20°時，滲吸采收率曲線提升效果并不顯著。滲吸采收率曲線的這種差異主要源于不同潤濕條件下孔隙尺度滲吸前緣不同的演化特征。中性潤濕條件下潤濕流體侵入程度低并很快達到穩定，滲吸驅替效果最差。弱潤濕(潤濕角為60°)和強潤濕(潤濕角≤37.3°)條件下，滲吸前緣主要沿著片狀孔隙和小孔隙形成的優勢路徑向多孔介質內部流動。由于同一模型內該優勢通道較為相似，相應地這三種潤濕條件下的最終采收程度較為接近。此外，由于潤濕角越小，毛管力越大，滲吸速率越快，因此強潤濕條件下的滲吸采收率曲線高于弱潤濕條件下的采收率曲線。當微觀孔隙內兩相流體達到穩定時，強、弱潤濕條件下模型從入口到出口范圍內均有大量非潤濕流體殘留。其中，單一孔隙內的非潤濕相液滴主要由于“卡斷”效應被捕獲，而位于模型中部較大體積的“團簇”狀非潤濕相主要由于“繞流”效應被捕獲。計算終止時刻不同潤濕條件下對應的最終采出程度如表1所示。潤濕角越小，自發滲吸最終采出程度越高。

圖7 2號巖心不同潤濕條件下采出程度隨時間的變化Fig. 7 The recovery degree of core 2 under different wetting conditions varies with time

圖8 2號巖心不同潤濕條件下自發滲吸速率隨時間的變化Fig. 8 The spontaneous imbibition rates of core 2 under different wetting conditions varies with time

表1 2號巖心不同潤濕條件下滲吸驅油效果Table 1 The oil displacement effect by imbibition of core 2 under different wetting conditions

圖8直觀展示了不同潤濕條件下滲吸速率隨時間的變化。滲吸初始階段，潤濕角越小，滲吸速率越大。強(潤濕角≤37.3°)、弱(潤濕角為 60°)和中性潤濕(潤濕角為80°)條件下的滲吸速率差異顯著，而強潤濕條件下(潤濕角≤37.3°)的滲吸速率的大小和變化趨勢較為接近。該變化趨勢表明，在強潤濕條件下進一步減小潤濕角并不能顯著提高滲吸速率。不同潤濕條件下自發滲吸持續時間顯著不同。強潤濕條件下滲吸速率快，潤濕流體快速突破并達到穩定狀態；中性潤濕條件下滲吸速率最慢，潤濕流體侵入程度低，滲吸持續時間最短；而弱潤濕條件下滲吸速率居中，滲吸持續時間最長。

4 結論

本研究通過致密砂巖的數字巖心模型進行了自發滲吸模擬計算，重點分析了孔隙結構特征和潤濕條件對自發滲吸的影響規律。結論如下：

(1)致密砂巖儲層的微觀孔隙結構特征顯著影響自發滲吸驅替機理和采收規律。形態特征以片狀發育且連通性較好的孔隙結構中滲吸速率快且非潤濕流體主要以“繞流”形式捕集，最終采出程度高；孔隙尺寸細小且連通性較好的孔隙結構內滲吸速率穩定，無較大波動，滲吸現象持續時間長，非潤濕流體主要以“卡斷”方式捕集，最終采出程度一般；形態特征以片狀發育但連通性較差的孔隙結構滲吸速率波動顯著，非潤濕流體主要以“卡斷”方式捕集，最終采出程度低。

(2)潤濕性顯著影響自發滲吸過程中滲吸前緣演化。潤濕角越小，潤濕流體以“角流形式”優先侵入孔隙角隅，兩相界面雜亂、分散，主終端液面滯后明顯，滲吸前緣后非潤濕相滯留明顯。

(3)潤濕性對滲吸初期的逆向滲吸及之后的順向滲吸過程影響顯著。逆向滲吸的發生程度和發生位置與潤濕性直接相關，強潤濕條件下，逆向滲吸發生程度高和發生位置多。順向滲吸過程中，潤濕性越強，滲吸作用越明顯，滲吸速率越快，滲吸采出程度越高。