基于格子玻爾茲曼方法的巖心切片滲流通道研究

劉博偉，牟松茹，李彥來，張 墨，張文童

（中海石油（中國）有限公司天津分公司渤海石油研究院，天津300459）

數(shù)字巖心技術是一種依托高精度成像設備在巖心孔隙級層面建模，繼而分析儲層物性參數(shù)的技術，目前已成為一種較為成熟、可用于油氣田開發(fā)工程的分析手段。相對于傳統(tǒng)室內巖心實驗，數(shù)字巖心以其實驗周期短、實驗可重復、微觀因素可控、節(jié)約巖心等優(yōu)勢被廣泛關注及應用。數(shù)字巖心掃描分辨率已達到納米級，將油氣流動研究推向微觀孔隙流動研究的大門[1-5]。

數(shù)字巖心孔隙流動具有異常復雜的孔隙空間，孔隙和骨架的邊界極不規(guī)則，直接求解Navier-Stokes方程不太符合實際，有限元、有限體積、有限差分等方法也難以處理異常復雜邊界[6-7]。格子玻爾茲曼方法（Lattice Boltzmann Method，LBM），作為介觀尺度計算流體力學方法，具有介于微觀分子動力學模型和宏觀連續(xù)模型的介觀模型特點，因其流體相互作用易描述、復雜邊界易設置、并行計算容易等優(yōu)勢，逐漸被應用到多孔介質復雜空間流動模擬[8-13]。

本文以渤海S油田東營組儲層數(shù)字巖心切片為研究對象，通過圖像處理完成孔隙網(wǎng)絡建模，利用格子玻爾茲曼方法開展巖心切片二維孔隙流動模擬及速度場表征，在此基礎上設計實驗，分析生產壓差、孔喉優(yōu)化對滲流通道的影響，探尋合適方法提高微觀孔隙滲流波及程度。

1 數(shù)字巖心切片流動模擬

實驗采用渤海S油田東營組某砂巖巖心，巖心儲層深度為1 379.0～1 379.1 m，巖心CT掃描結果如圖1所示，其中黑色部分為孔隙，灰色部分為骨架。圖1（a）為全直徑巖心掃描結果，對應掃描分辨率為 100 μm；圖 1（b）為巖心柱塞掃描結果，對應掃描分辨率為 13 μm；圖 1（c）為柱塞子樣掃描結果，對應掃描分辨率為2 μm。

圖1 渤海S油田某儲層數(shù)字巖心掃描結果Fig.1 Digital core scanning results of a reservoir in Bohai S oilfield

根據(jù)巖心掃描切片可知，全直徑巖心由于掃描精度很低，無法有效識別巖心孔隙及裂縫；巖心柱塞掃描分辨率中等，巖心骨架邊緣灰度相近，識別存在誤差，影響巖心孔隙分析；柱塞子樣掃描分辨率較高，能夠有效識別骨架邊緣、巖心孔隙。因此，本文基于柱塞子樣掃描結果，隨機選取巖心某一切片開展孔隙網(wǎng)絡建模及切片流動模擬（見圖2），分析切片滲流情況。

圖2 數(shù)字巖心某二維切片F(xiàn)ig.2 Digital two-dimensional slice

1.1 巖心切片圖像處理

在數(shù)字巖心分層CT掃描切片過程中，受掃描設備影響不可避免地會產生一定程度的噪聲，干擾骨架邊界識別，因此需要開展灰度圖像降噪處理。常用降噪處理方式有均值濾波和中值濾波。由于中值濾波屬非線性降噪技術，能夠相對有效地保護邊緣信息，在數(shù)字巖心灰度圖像降噪處理的同時兼顧保護骨架邊界的不規(guī)則形狀。因此，選用中值濾波技術對數(shù)字巖心掃描圖像進行噪聲過濾，降噪效果見圖3。

圖3 切片局部原始及過濾后灰度圖像Fig.3 Local original and filtered grayscale image of a core plunger

數(shù)字巖心CT掃描結果為灰度圖像堆棧，而用于流動模擬的孔隙網(wǎng)絡模型僅關注骨架、孔隙兩種狀態(tài)，因此需要開展閾值分割，將灰度圖像二值化，用1代表骨架，0代表孔隙。采用單閾值處理方式，數(shù)字巖心數(shù)據(jù)體設定單一閾值，灰度值大于閾值視為骨架，小于閾值視為孔隙。選用試湊法以儲層孔隙度為依據(jù)進行閾值選取，當閾值取76時，數(shù)字巖心孔隙度與儲層實際值一致，以此為閾值進行圖像二值化，二值化處理效果見圖4。

圖4 數(shù)字巖心切片圖像二值化處理Fig.4 Digital core slice image binarization

1.2 格子玻爾茲曼方法模型

基于Y.H.Qian等[14]所提出的DdQm系列格子玻爾茲曼方法模型開展二維切片流動模擬，本文僅選用D2Q9系列模型，通過C++語言編程實現(xiàn)模型計算。

格子玻爾茲曼模型主要由三部分組成，即格子（離散速度模型）、平衡態(tài)分布函數(shù)、分布函數(shù)的演化方程[15-16]。

數(shù)字巖心切片為二維二值化圖像，因此離散速度模型選用D2Q9模型（即適用范圍為二維空間，模型包含9個離散速度方向，見圖5），離散速度配置見式（1）。

圖5 D2Q9離散速度模型Fig.5 D2Q9 discrete velocity model

其中，α=0,1,2,…,8代表D2Q9模型離散速度方向；eα代表對應α方向上的速度；c代表單位格子速度，即格子步長與時間步長之比。

平衡態(tài)分布函數(shù)選用Y.H.Qian等[14]提出的LBGK模型（由于本文所研究的切片滲流僅局限于單相流體流動，因此不考慮附加作用力項），平衡態(tài)分布函數(shù)見式（2）、（3）。

平衡態(tài)分布函數(shù)對應的演化方程見式（4）。

其中，r代表流體粒子的空間位置；t代表時間；δt代表時間步長。

根據(jù)LBGK模型理論，宏觀密度、速度定義見式（5）。

其中，ρ′和u′分別表示模型的宏觀密度、速度。

數(shù)字巖心切片流動模擬對應的孔隙網(wǎng)絡模型呈矩形（見圖6），入口端/出口端選用壓力邊界，頂面/底面及孔隙骨架邊界選用標準反彈邊界。標準反彈邊界，指在某微粒格點上，微粒按照運動軌跡運動至邊界格點，其后微粒發(fā)生反彈，反彈后沿原運動軌跡相反的方向運動；標準反彈邊界常用于描述無滑移壁面反彈，即沿壁面方向無速度梯度的反彈。壓力邊界條件，參考Q.S.Zou等[17]根據(jù)非平衡反彈思想提出的動力學壓力邊界格式假設；在孔隙網(wǎng)絡模型流入端壓力值已知，即 Pin=cs2ρin（Pin為流入面壓力，ρin為流入微粒密度）已知，而流入端速度未知,邊界處沿壁面方向速度分量為零，水平速度分量ux未知，結合離散速度模型及平衡態(tài)分布函數(shù)與宏觀變量的關系，聯(lián)立方程組即可求解邊界速度。

圖6 孔隙網(wǎng)絡模型示意Fig.6 Schematic diagram of pore network model

1.3 切片滲流過程表征

基于格子玻爾茲曼方法切片流動模擬結果，計算孔隙內流體速度場分布情況，通過對比滲流速度大小區(qū)分判斷孔隙優(yōu)勢滲流通道。圖7為切片孔隙內流體速度場隨時間變化情況。當速度場不再隨時間發(fā)生變化時，達到穩(wěn)態(tài)。本文優(yōu)勢通道研究均是在滲流達到穩(wěn)態(tài)條件下開展的，選取相鄰時間步長速度變化量小于10-3為穩(wěn)態(tài)判斷條件。

2 巖心切片優(yōu)化滲流通道分析

基于格子玻爾茲曼方法，結合達西滲流理論，開展不同生產壓差、不同孔喉尺寸條件下流體穩(wěn)態(tài)速度場分布研究，分析二維切片優(yōu)勢滲流通道的影響因素，開展優(yōu)化孔隙滲流方法探索。

數(shù)字巖心二維切片宏觀尺寸為1.00 mm×0.56 mm，格子尺寸為500×280；實驗模擬流體為氣體，宏觀密度為1.293 kg/m3。

2.1 生產壓差對滲流影響

設置入口端/出口端生產壓差分別為0.005、0.010、0.015、0.020、0.025 MPa，其他條件不變，可獲得對應穩(wěn)態(tài)速度場分布，圖8為不同壓差條件下滲流通道表征圖。

根據(jù)圖8可知，伴隨生產壓差逐漸增大，切片孔隙內的滲流速度逐漸增大，流體的滲流通道也趨于穩(wěn)定。圖8（a）中，由于生產壓差較小，滲流速度較小，但能夠辨別出滲流通道為主通道；圖8（b-e）中，滲流通道仍以主通道為主，但隨著壓差進一步增大，次級通道也有較弱的滲流。

圖7 巖心二維切片速度場Fig.7 Core two-dimensional slice velocity field diagram

圖8 不同壓差條件下滲流通道圖Fig.8 Seepage channel diagram under different pressure difference conditions

根據(jù)模擬結果可知，僅通過放大生產壓差，次級通道滲流速度較小，無法高效驅替局部微觀剩余油。

2.2 孔喉優(yōu)化對微觀滲流影響

利用數(shù)字巖心微觀因素可控的優(yōu)勢，結合二維切片優(yōu)勢滲流通道展布形態(tài)，類比現(xiàn)有的油田增產技術（儲層酸化壓裂，即擴大巖心孔縫；注采井網(wǎng)調驅，即封堵孔喉改變流場），有針對性地改造二維切片局部骨架形態(tài)，測試其對孔隙微觀滲流的影響。

假設能夠實現(xiàn)儲層巖石微觀骨架孔喉封堵，以此為基礎開展二維切片骨架改造。將原二維切片主通道入口進行孔喉封堵，結果如圖9所示。根據(jù)穩(wěn)態(tài)滲流場分布可知，孔喉封堵能夠有效改變切片滲流通道分布，難滲流區(qū)域開始形成有效滲流，達到增大微觀孔隙波及程度的目的。然而，次級通道入口孔徑較小，封堵優(yōu)勢通道也造成滲流速度變小，不利于儲層滲流。

圖9 封堵局部孔喉條件下滲流通道Fig.9 Seepage channel diagram under local pore throat conditions

針對微觀孔喉封堵后造成的滲流速度變小的問題，提出微觀壓裂酸化的假設，以此進行切片骨架孔喉擴徑改造，增大孔隙滲流速度。在封堵主通道入口的基礎上，將局部區(qū)域孔喉擴徑，增大次級通道入口的半徑，結果如圖10所示。根據(jù)穩(wěn)態(tài)滲流場分布可知，難滲流區(qū)域開始形成高效滲流，切片微觀波及程度進一步擴大，達到擴大微觀儲層改造的目的。

基于二維切片骨架改造的假設，僅針對圖10中次級通道入口進行改造，開展不同孔徑大小對孔隙滲流速度的影響。模擬結果如圖11所示，當次級通道入口孔徑相對較小時，其對滲流速度的制約作用越大；當孔徑擴大到一定程度時，滲流速度增幅逐漸減小。

圖10 局部孔喉擴徑條件下滲流通道圖Fig.10 Seepage channel diagram under local pore throat expansion condition

圖11 次要通道孔徑大小對滲流影響曲線Fig.11 Influencing curve of secondary channel aperture size on seepage flow

2.3 優(yōu)化微觀滲流影響

儲層開發(fā)初期，可逐步增大生產壓差，實現(xiàn)儲層主力滲流通道孔隙中的原油開發(fā)，但仍存在微觀局部剩余油；開發(fā)中后期，主力滲流通道驅油效率達到較高水平，此時可考慮將該類通道封堵，再輔助孔喉擴徑措施，實現(xiàn)滲流通道優(yōu)化，挖潛巖石微觀未波及孔隙的剩余油。

3 結 論

（1）基于數(shù)字巖心掃描數(shù)據(jù)結果，利用圖像處理建立孔隙網(wǎng)絡模型，通過格子玻爾茲曼方法，選用單松弛模型、壓力邊界、標準反彈邊界，實現(xiàn)巖心切片二維流動滲流通道表征。

（2）設計實驗分析生產壓差、孔喉優(yōu)化對滲流通道的影響；模擬結果表明，生產壓差能夠增強滲流速度，但無法改變滲流通道；孔喉優(yōu)化能夠有效改變滲流通道擴大孔隙波及程度。