李昂, 于浩波, 謝斌, 黃波, 馮耀榮, 宮敬, 郝鵬飛

(1.中國石油大學(北京)油氣管道輸送安全國家工程實驗室, 北京 102249; 2.大慶油田采油工程研究院, 黑龍江 大慶 163453; 3.中國石油新疆油田工程技術(shù)研究院, 新疆 克拉瑪依 834000; 4.清華大學航天航空學院, 北京 100084)

0 引 言

致密油藏由于其儲量可觀,開發(fā)潛力大,20世紀中期在國外已經(jīng)受到重視并開采,取得了巨大成功,如竇宏恩等[1]介紹的美國巴肯油藏。在微尺度下研究致密油儲層的孔隙結(jié)構(gòu)及流體在其中的傳熱傳質(zhì)對于儲層的高效開發(fā)及環(huán)保領(lǐng)域有著重要意義[2-3]。

致密儲層巖心屬于結(jié)構(gòu)復(fù)雜的多孔介質(zhì),研究流體在其中的流動首先要獲取其內(nèi)部的三維空間結(jié)構(gòu),然后才能進行流動模擬。中國學者對國外研究成果進行了總結(jié)對比后均指出,使用X射線對樣品進行CT掃描建立數(shù)字巖心是一種準確而有效的獲取巖心內(nèi)部空間結(jié)構(gòu)的手段[4-6]。

對孔隙中的微觀流動進行數(shù)值模擬是近年國內(nèi)外學者研究的熱點,對于其模擬方法有基于孔隙網(wǎng)絡(luò)模型法(PNM)[7]、格子Boltzmann法(LBM)[8-11]及計算流體力學(CFD)方法[12-18]。PNM方法,在簡化模型的基礎(chǔ)上進行流動模擬,得出的結(jié)果很可能與實際狀況產(chǎn)生較大誤差[12]。LBM方法,雖然能夠處理復(fù)雜的邊界條件,編程簡單,但是運算量巨大,計算耗時長[13]。運用CFD方法是孔隙流動模擬的發(fā)展趨勢。當前的研究成果主要集中在孔隙中的流動上,還很少有涉及到其中傳熱問題的研究,而Wu M等[19]與Yang L等[20]指出研究其中的傳熱狀態(tài)也很有必要,而且其研究對象主要是孔隙度較高、滲透率較大的巖心模型,比如Berea、Fontainebleau砂巖模型以及Sandpack填充模型等。關(guān)于低滲透率的致密巖心(滲透率<1.0×10-3μm2)的研究還很少。

本文選取低滲透率致密巖心作為研究樣品,結(jié)合Micro-CT掃描與后續(xù)圖像處理及網(wǎng)格劃分技術(shù),使用基于有限體積法(FVM)的CFD軟件ANSYS FLUENT計算巖心的速度與溫度分布狀況,為研究巖心類多孔介質(zhì)中的流動與傳熱問題提供一種新方法。

1 巖心三維模型的建立

CT掃描是利用穿透性好的X光拍攝巖心不同角度的透射圖像,再利用層析成像算法反演出每個成像點的透射率值,該值大小直接反映出該點的物質(zhì)密度;以白色表示高密度物質(zhì),以黑色表示孔隙,形成灰度圖像。

本文中使用ZEISS公司生產(chǎn)的Versa510對樣品進行成像,該設(shè)備最高分辨率為300 nm。巖心樣品來自新疆油田吉木薩爾31號井,為白云質(zhì)巖,在樣品中鉆取直徑約1 mm圓柱進行掃描實驗,可以獲得分辨率為916×994的二維灰度圖像970張[見圖1(a)],空間分辨率為1 μm/像素。在二維灰度圖像中可以清楚地看到,黑色的部分為孔隙(低密度),灰白色的部分為骨架(高密度)。通過二維圖像的疊加可以得到樣品的三維灰度圖像。將所有圖像用于重建數(shù)字巖心對計算機的性能有非常高的要求,本研究中使用Mostaghimi P等計算出的數(shù)字巖心的代表元體積(REV)來確定重構(gòu)區(qū)域尺寸[15],即從圖像中選取體積為200 μm×200 μm×200 μm的區(qū)域作為研究對象。圖像處理在三維可視化軟件AVIZO中進行。

圖1 CT掃描實驗說明

使用Micro-CT掃描獲得的灰度圖像存在各類噪聲,會對結(jié)果產(chǎn)生影響。進行圖像處理之前首先需要對灰度圖像進行降噪。研究中使用降噪效果比較好的中值濾波法對原始圖像進行降噪;之后對二維灰度圖像進行閾值分割,根據(jù)灰度值的差異將孔隙空間與巖石骨架區(qū)分開,圖1(b)為閾值分割后二維切片圖。將閾值分割后的圖像從下往上疊加,再使用Marching Cubes算法產(chǎn)生三維孔隙表面,去除孤立孔隙后通過提取孔隙結(jié)構(gòu)的中心線建立起表征孔隙空間簡化結(jié)構(gòu)的網(wǎng)絡(luò)模型,三維重構(gòu)結(jié)果以及孔隙網(wǎng)絡(luò)模型見圖2。

巖心的三維重構(gòu)結(jié)果可以直觀地看出孔隙的分布狀況,但是要深入認識孔隙結(jié)構(gòu)的特點還需要對其特征參數(shù)進行統(tǒng)計分析。首先在重建模型中使用的快速分水嶺算法將每個孔隙進行標記,所有孔隙的體積和等效直徑便可以統(tǒng)計出來,再除以重構(gòu)區(qū)域體積便可以算出重建模型孔隙度大小為8.34%。與實測孔隙度8.73%比對發(fā)現(xiàn),計算孔隙度稍小于實測孔隙度,可能原因為在圖像處理中將一部分小的孤立孔隙去除所導(dǎo)致,巖心的非均質(zhì)性也可能會造成計算值與實測值存在一些誤差[12]?？紫逗淼赖闹睆郊绑w積的分布見圖3。

圖2 數(shù)字巖心構(gòu)建結(jié)果

圖3 孔隙喉道直徑與體積分布

2 計算過程

圖4 網(wǎng)格及邊界條件

CFD計算方法主要包括有限差分法(FDM)、有限元法(FEM)以及有限體積法(FVM)。其中,FVM是最為普遍的一種數(shù)值計算方法,計算結(jié)果準確且適用于復(fù)雜幾何結(jié)構(gòu),通用性較廣的CFD軟件中大多都采用該方法[21],然而該種方法在孔隙級的流動模擬中使用卻很少。為此,需搭建重構(gòu)模型與基于FVM的計算軟件ANSYS FLUENT對接的橋梁。

首先對幾何模型進行網(wǎng)格劃分,該步驟在AVIZO軟件中的Generate Tera Grid模塊中進行。在產(chǎn)生網(wǎng)格之前,必須先對重構(gòu)數(shù)字巖心進行表面處理,去除孔隙空間表面交叉的部分并保持幾何模型的閉合性,使之能夠產(chǎn)生網(wǎng)格,經(jīng)表面處理后產(chǎn)生的網(wǎng)格數(shù)約為537萬個;在網(wǎng)格產(chǎn)生后還需要對產(chǎn)生的非結(jié)構(gòu)化四面體網(wǎng)格進行優(yōu)化,以免在后續(xù)迭代計算過程中出現(xiàn)不收斂或者運算超出浮點數(shù)。研究以X、Y、Z這3個坐標軸為流動方向進行計算,不同流動方向進出口邊界條件亦不同,需要分別設(shè)置,如圖4所示(以X軸流動方向為例),其中圖4(b)中的紫色平面和圖4(c)中的金黃色平面分別為流體的進口和出口。計算時進口壓力分別設(shè)置為10 Pa,出口壓力設(shè)置為0 Pa,流體設(shè)置為液態(tài)水,進口溫度設(shè)為300 K,壁面溫度設(shè)置為310 K,邊界條件如圖4(d)(以X軸流動方向為例)。由于流體在孔隙中的流速較低,故在計算前假設(shè)其中的流動為層流狀態(tài),即低雷諾數(shù)Stokes流,且壁面無滑移,其流動控制方程為

·V=0

(1)

μ2V-p=0

(2)

式中,V為流體在多孔介質(zhì)中的速度矢量;μ為流體動力黏度;為拉普拉斯算子;p為多孔介質(zhì)中流體的壓力。溫度場則通過求解能量方程得出。

經(jīng)過求解器的迭代計算,在一定壓差下水通過巖心孔隙空間的體積流量的數(shù)值解便可以求解出來,并使用達西定律計算出絕對滲透率。在達西定律中,絕對滲透率視為常數(shù),只與孔隙的幾何形狀有關(guān),是巖石本身的性質(zhì),與通過的流體性質(zhì)無關(guān)。達西定律的表達形式為

(3)

式中,Q為通過多孔介質(zhì)的流體的體積流量,μm3/s;S為流體流動方向上樣品的截面積,μm2;K為絕對滲透率,μm2;μ為流體的動力黏度,Pa·s;Δp為多孔介質(zhì)兩側(cè)壓差,Pa;L為流動方向上樣品的長度,μm;Q/S為通過多孔介質(zhì)的滲流速度,即達西速度。

圖5 三維滲流流線圖

3 計算結(jié)果及分析

3.1 絕對滲透率計算結(jié)果

重構(gòu)區(qū)域內(nèi)的3個不同方向的絕對滲透率數(shù)值計算結(jié)果與實測滲透率見表1。從滲透率實測值可知實驗巖心屬于白云質(zhì)巖中滲透率較高的樣品;同時每個流動方向的滲透率都互不相同,且理論值與實驗結(jié)果相比存在一些誤差,杜新龍[12]提出這種誤差是受到巖心的非均質(zhì)性以及掃描的分辨率等影響所致,但理論計算結(jié)果與實驗值均在同一數(shù)量級并且已較為接近,對于孔隙級別的流動來說該誤差在可以接受的范圍內(nèi)。與用有限元軟件COMSOL來計算的方法相比,雖然使用COMSOL計算得到的結(jié)果非常接近,但是使用ANSYS FLUENT進行迭代計算時收斂速度更快,在同一臺計算機上達到相同計算精度的時間只需要COMSOL的1/5左右,在孔隙級別流動模擬計算中使用ANSYS FLUENT進行計算更有優(yōu)勢。

表1 巖心滲透率計算結(jié)果

3.2 滲流流場分析

將計算所得數(shù)據(jù)進行后處理作出的三維流線圖見圖5所示。從流線圖5可以直觀地觀測到流體在孔隙內(nèi)部的分布狀況,有流線穿過的區(qū)域表明該區(qū)域有流體流動,流線顏色靠近紅色表示該處流速越快,越靠近藍色則表示流速越慢。

之前算出X、Y、Z這3個流動方向絕對滲透率值互不相同的原因可以從流線分布狀況來說明。若將孔隙內(nèi)部流動通道等效為毛細管,則根據(jù)Poiseuille定律可知,在進出口端壓差、通道長度與流體性質(zhì)相同的情況下,通過該通道的流量與通道的尺寸相關(guān)。

(4)

式中,Q為流量;R為圓管半徑;Δp為兩端壓差;η為液體粘滯系數(shù);l為圓管長度。通過三維滲流流線圖可以看出,各流動方向上流體都會匯集于一段孔隙通道中(見圖5),即存在匯流通道,也可稱之為流動控制喉道。在每段匯流通道中等距離取10個與流線相垂直的截面,并對該截面進行積分計算,可以得出X、Y、Z流動方向的匯流通道平均等效半徑Rx、Ry、Rz分別為10.46、16.14、21.37 μm,等效半徑比Rx∶Ry∶Rz=1∶1.39∶1.84,Rx4∶Ry4∶Rz4=1∶5.67∶17.4,而3個流動方向絕對滲透率之比Kx∶Ky∶Kz=1∶6.3∶16.3,與Poiseuille定律所表現(xiàn)的趨勢基本吻合,可以認為不同流動方向上絕對滲透率產(chǎn)生差異的原因是匯流通道的尺寸大小不同。

圖6 X流動方向不同截面溫度分布

3.3 溫度場分析

通過流動方向上不同截面孔隙內(nèi)溫度云圖可表示溫度分布情況。以X流動方向為例,圖6顯示了距進口平面0、4、8、12、16 μm截面孔隙內(nèi)流體溫度分布狀況,各個截面的最大溫差為9.7、6.1、3.1、1.4、0.4 ℃,可以看出沿著該流動方向流體被迅速加熱至壁面溫度附近;同時通過設(shè)置將進口壓力設(shè)置為100 Pa和500 Pa再計算后發(fā)現(xiàn),上述各截面溫度分布情況變化非常小,最大溫差值基本不變。表明孔隙內(nèi)部有著非常高的流動傳熱效率,主要是由于孔隙尺度小且流體流速低導(dǎo)致熱混合充分。

4 結(jié) 論

(1) 使用Micro-CT技術(shù)對吉木薩爾致密油儲層巖心樣品進行掃描獲取了巖心三維灰度圖像,結(jié)合圖像處理技術(shù),解決了網(wǎng)格劃分及優(yōu)化的難題,成功將幾何模型導(dǎo)入基于FVM的CFD軟件ANSYS FLUENT中進行流動與傳熱耦合計算,與有限元軟件COMSOL比,收斂速度大大加快,為研究孔隙中的流動與傳熱現(xiàn)象提供了一種新的途徑。

(2) 通過對巖心內(nèi)部滲流流場進行分析,發(fā)現(xiàn)在REV大小的重構(gòu)區(qū)域內(nèi)部流場存在一段特殊的匯流通道,巖心絕對滲透率大小與匯流通道的平均等效半徑有關(guān),其規(guī)律與Poiseuille定律所表現(xiàn)的趨勢基本吻合。

(3) 通過計算在不同進口壓力下溫度場分布,獲取了不同截面最大溫度差;發(fā)現(xiàn)沿著流動方向上流體被迅速加溫至壁面溫度附近,說明孔隙級流動中有著很高的傳熱效率。

參考文獻:

[1] 竇宏恩, 馬世英. 巴肯致密油藏開發(fā)對我國開發(fā)超低滲透油藏的啟示 [J]. 石油鉆采工藝, 2012, 34(2): 120-124.

[2] OVAYSI S, PIRI M. Direct Pore-level Modeling of Incompressible Fluid Flow in Porous Media [J]. Journal of Computational Physics, 2010, 229(19): 7456-7476.

[3] BLUNT M J, BIJELJIC B, HU D, et al. Pore-scale Imaging and Modeling [J]. Advances in Water Resources, 2013, 51: 197-216.

[4] 黃豐. 多孔介質(zhì)模型的三維重構(gòu)研究 [D]. 合肥: 中國科學技術(shù)大學, 2007.

[5] 趙秀才. 數(shù)字巖心及孔隙網(wǎng)絡(luò)模型重構(gòu)方法研究 [D]. 北京: 中國石油大學, 2009.

[6] 劉學鋒, 張偉偉, 孫建孟. 三維數(shù)字巖心建模方法綜述 [J]. 地球物理學進展, 2013, 28(6): 3066-3072.

[7] DONG H, BLUNT M J. Pore-network Extraction from Micro-computerized-tomography Images [J]. Physical Review E, 2009, 80(2): 1957-1974.

[8] CHEN S, DOOLEN G D. Lattice Boltzmann Method for Fluid Flows [J]. Annual Review of Fluid Mechanics, 2003, 30(1): 329-364.

[9] OKABE H, BLUNT M J. Prediction of Permeability for Porous Media Reconstructed using Multiple-point Statistics [J]. Physical Review E, 2005, 70(6): 66-135.

[10] BOEK E S, Venturoli M. Lattice-Boltzmann Studies of Fluid Flow in Porous Media with Realistic Rock Geometries [J]. Computers & Mathematics with Applications, 2010, 59(7): 2305-2314.

[11] COON E T, PORTER M L, KANG Q. Taxila LBM: A Parallel, Modular Lattice Boltzmann Framework for Simulating Pore-scale Flow in Porous Media [J]. Computational Geosciences, 2013, 18(1): 1-11.

[12] 杜新龍. 低滲透砂巖油層微流動機理研究 [D]. 成都: 西南石油大學, 2012.

[13] 朱洪林. 低滲砂巖儲層孔隙結(jié)構(gòu)表征及應(yīng)用研究 [D]. 成都: 西南石油大學, 2014.

[14] LIU Z, WU H. Pore-scale Study on Flow and Heat Transfer in 3D Reconstructed Porous Media Using Micro-tomography Images [J]. Applied Thermal Engineering, 2016, 100: 602-610.

[15] MOSTAGHIMI P, BLUNT M J, BIJELJIC B. Computations of Absolute Permeability on Micro-CT Images [J]. Mathematical Geosciences, 2013, 45(1): 103-125.

[16] ZHANG S, MAESTRA F D, COMBARET N, et al. The Analysis and Simulation of Rock Properties Using FIBSEM and Virtual Material Studio [M]. NAFEMS World Congress, 2011.

[17] 劉向君, 朱洪林, 梁利喜. 基于微CT技術(shù)的砂巖數(shù)字巖石物理實驗 [J]. 地球物理學報, 2014, 57(4): 1133-1140.

[18] SIENA M, HYMAN J D, RIVA M, et al. Direct Numerical Simulation of Fully Saturated Flow in Natural Porous Media at the Pore Scale: a Comparison of Three Computational Systems [J]. Computational Geosciences, 2015, 19(2): 423-437.

[19] WU M, LI M, XU C, et al. The Impact of Concrete Structure on the Thermal Performance of the Dual-media Thermocline Thermal Storage Tank Using Concrete as the Solid Medium [J]. Applied Energy, 2014, 113(1): 1363-1371.

[20] YANG L, SHEN H. Effects of the Porous Media Distribution on the Performance Improvement for Isothermal Chamber [J]. Applied Thermal Engineering, 2015, 86(1): 301-308.

[21] 陶文銓. 數(shù)值傳熱學 [M]. 2版. 西安: 西安交通大學出版社, 2001.