基于Level set方法的微觀水驅油模擬分析*

高亞軍 姜漢橋 王碩亮 李俊鍵 劉傳斌 常元昊 王依誠

(1. 中國石油大學(北京)石油工程教育部重點實驗室 北京 102249； 2. 中國地質大學(北京)能源學院 北京 100083)

基于Level set方法的微觀水驅油模擬分析*

高亞軍1,2姜漢橋1王碩亮2李俊鍵1劉傳斌1常元昊1王依誠1

(1. 中國石油大學(北京)石油工程教育部重點實驗室 北京 102249； 2. 中國地質大學(北京)能源學院 北京 100083)

高亞軍,姜漢橋,王碩亮，等.基于Level set方法的微觀水驅油模擬分析[J].中國海上油氣，2016,28(6)：59-65.

Gao Yajun,Jiang Hanqiao,Wang Shuoliang,et al.Simulation analysis of microscopic water-oil displacement based on Level set method[J].China Offshore Oil and Gas,2016,28(6):59-65.

通過引入Level set方法，構建兩相界面運移數學表征方法，結合Navier-Stokes方程建立兩相微觀滲流數學模型，借助有限元方法進行求解，并利用微觀玻璃刻蝕模型驅替實驗驗證數值模擬結果。結果表明，基于該滲流數學模型的數值模擬結果與實驗結果的吻合度很高。利用Navier-Stokes方程并結合Level set方法實現了兩相流體驅替的微觀數值模擬，從而為研究微觀兩相滲流特征及驅油機理提供了一種新方法。在不等徑并聯喉道單元模型基礎上，進行了不同驅替條件下的數值模擬計算，分析了發生孔隙尺度兩相竄流的臨界喉道直徑比的變化特征，結果表明驅替壓差及油水兩相黏度比越大，越容易發生微觀竄流，更易導致在較小喉道中形成微觀殘余油。本文研究結果為微觀水驅油機理的進一步研究奠定了基礎。

數值模擬；Level set方法；微觀水驅油；微觀殘余油

目前許多學者在微觀水驅油兩相流動特征方面做了大量的實驗研究，實驗方法主要有三維平板模型[1]、玻璃刻蝕物理模擬等，但這些實驗模擬成本高，操作較為繁瑣困難，并且由于實驗監測技術和圖像獲取精度有限，導致實驗數據處理誤差較大。近年來，微觀滲流的數值模擬方法迅速發展，已形成網絡模擬和玻爾茲曼法兩大類，傳統模型中孔隙結構已用柱狀喉道和互相連接的球狀孔隙代替。Fatt[2]最早用孔隙網絡模型研究多孔介質二維網格中的流體動靜態問題，隨后網絡模型被不斷改進。盡管孔隙網絡模型可以很好地計算相對滲透率和殘余油飽和度[3-4]，但是孔喉結構太理想化，無法解釋流體動態機理[5]。玻爾茲曼法可以很好模擬巖石孔隙的兩相滲流，定性符合較好[6-7]，但也無法刻畫孔隙內部油水兩相界面動態運移特征。

在微觀水驅油研究中，油水兩相流動特征及剩余油分布狀態均是通過孔隙內油水兩相界面體現出來的。在Navier-Stokes方程的基礎上，計算流體動力學方法在兩相流領域的模擬應用逐漸發展起來[8]，所形成的兩相流體相界面動態特征的研究方法主要有界面追蹤類和界面捕獲類，其中以Level set方法為代表的界面捕獲類方法應用較多[9]。Level set方法是由Osher和Sethian提出[10]，后經Sussman等和Olsson等的發展改進[11-13]，保持了計算過程中的體積守恒。Level set方法在孔隙尺度油水兩相流領域應用較少,僅用來模擬研究孔隙尺度下水驅和CO2驅采收率的比較及孔喉中流體相界面的動態驅替特征[14-15],而且現有文獻均是直接建立孔隙模型進行數值模擬，模擬的實際精確度有待實驗驗證。

本文利用巖心切片掃描圖像構建真實孔隙拓撲分布結構，建立二維微觀孔隙模型，結合Navier-Stokes方程，用Level set方法進行數值模擬，模擬結果與可視化微觀玻璃刻蝕實驗結果達到高度吻合，證明了模擬結果的可靠性，克服了上述常規實驗方法和模擬方法存在的諸多問題。通過此模型，模擬計算殘余油富集孔道直徑及孔道內流速和壓力的分布特征，并在不等徑并聯孔道單元基礎上分析在大孔道存在的情況下小孔道內原油被驅動的孔徑條件，為微觀孔隙驅油機理的進一步研究奠定了基礎。

1 水驅實驗模型

如圖1所示，在光化學玻璃上雕蝕出兩種不同孔隙結構的模型A和B，均為中性潤濕，其中模型A孔喉大于模型B，模型A的微觀孔喉非均質性也強于模型B，且模型A中部存在明顯的優勢滲流通道。模型A的孔隙度為30.96%，尺寸為17.45 mm×14.31 mm；模型B的孔隙度為28.91%，尺寸為15.28 mm×12.13 mm。用另一塊玻璃進行粘結密封，制作成仿真巖心微觀孔隙模型，然后進行水驅油實驗。實驗步驟如下：

1) 安裝固定玻璃刻蝕模型；

2) 將模型A和B抽真空飽和原油,對應黏度分別為50 mPa·s和20 mPa·s ，并采集原油穩定后的圖像；

3) 在室溫下啟動定壓泵從模型左端注水(實驗用水密度為1 g/cm3，黏度為1 mPa·s)，進行水驅油過程，直到出口端基本不出油，用顯微放大系統和錄像機錄取實驗過程的圖像，并用計算機圖象處理系統對錄像和圖片進行處理。

圖1 飽和不同原油的玻璃刻蝕模型

2 數學模型

建立二維幾何模型，在孔隙尺度兩相流動中用Navier-Stokes方程來描述微觀兩相流體的動量傳輸特性，用Level set算法來刻畫在不同驅替時刻兩相界面動態變化，二者聯立建立微觀兩相滲流數學模型。其假設條件為：①孔隙內水油兩相流體均為不可壓縮流體，流動狀態為層流；②忽略流體重力；③注入端各孔隙入口壓力相等且恒定不變；④忽略玻璃刻蝕模型的厚度，流體在孔隙中屬于二維流動。

2.1 幾何模型

將玻璃刻蝕孔隙模型A和B在電鏡下的掃描圖進行二值化處理后，在AutoCAD軟件中轉化為孔隙輪廓圖，保持模型尺寸不變。如圖2所示，封閉曲線表示巖石顆粒，剩余聯通區域為孔隙喉道；注入端均在右端，出口端均在左端。

圖2 微觀孔隙幾何模型

2.2 Level set方法

在層流不可壓縮非混相兩相流中，用Level set方程中的φ來表示油水兩相流體各自的體積分數。φ=50%時的水平集輪廓線表示兩相流體相界面。對于油水兩相時，油相中φ值為0，水相中φ值為100%，水驅過程中φ值可以追蹤油水兩相界面。因此，兩相流體的界面運移方程可以表示為[16]

(1)

式(1)中：φ為油水兩相界面輪廓線；γ為方程求解中的重新初始化參數；t為兩相作用時間；u為流速；ε為人為設置的兩相界面厚度，一般小于計算模型中網格剖分的最小單元。

借助于Level set函數φ，流體的物性可用下述方程表示：

(2)

式(2)中：ρo和ρw分別為油和水的密度；μo和μw分別為油和水的黏度。

2.3 控制方程

Navier-Stokes方程可以用于描述不可壓縮流體質量和動量的傳輸特性。在固定的歐拉坐標系中，不可壓縮且考慮界面張力的兩相流動Navier-Stokes方程和連續性方程描述如下[12]：

(3)

式(3)中：ρ為密度，kg/m3；μ為動態黏度，Ns/m2；u為速度，m/s；p為入口壓力，Pa；I為單位矩陣；Fst為氣液界面張力，N/m。

兩相流體的界面張力可以用如下連續界面力學模型中的形式來表達[17]：

Fst=σknΓδsm

(4)

(5)

(6)

(7)

式(4)～(7)中：σ是界面張力系數；k是油水相界面曲率；nΓ是單位法向量；δsm是狄克拉函數[15]。

2．4 初始和邊界條件

如圖2所示，模型A右端有2個孔隙入口，左端有2個孔隙出口；模型B右端有7個孔隙入口，左端有9個孔隙出口。進出口邊界壓力在模型求解時設置，左右端除進出口外的邊界及頂部和底部邊界均為封閉邊界，無流體流動。初始界面設置為t=0時水相和油相的接觸面。出口和入口初始狀態方程式為

(8)

(9)

孔隙壁面條件均為中性潤濕，垂直壁面方向流速為0，方程式如下：

u·nwall=0

(10)

速度在壁面方向的法向分量為零，但由于壁面的潤濕性，速度在壁面方向的切線分量不為零，有一個極小的速度滑移。這個極小的速度分量可用如下的附加摩擦力來表示：

(11)

式(11)中：Ffr為摩擦力；μ為動力黏度；β為滑移長度，取值等于邊界部位剖分網格的尺寸。

3 模型驗證及模擬結果分析

基于上述數學模型，對這兩個微觀孔隙結構不同的模型進行求解，并進行數值模擬和實驗結果的對比，進而分析喉道內流體流動特征。

3.1 數值模擬與實驗結果對比

數學模型通過有限元方法進行求解。幾何模型用自由三角形網格差分，網格剖分數分別為93 121和46 906個。初始狀態孔隙中均飽和油，入口端注入水，水和油的密度、黏度及注入壓力如表1所示。為了保證計算過程的收斂性及節省計算時間，入口壓力與實驗壓力保持一致。孔隙內部初始壓力均為零。將數值模擬結果與實驗結果進行對比，以驗證模擬結果的可靠性。

如圖3、4所示，在不同含水飽和度條件下將模型A和B的驅替實驗所錄取的圖像與數值模擬結果進行對比，發現二者驅替相界面吻合度非常高，孔道內油水前緣相界面越尖銳，滲流阻力越小；相界面曲率越小，孔隙驅替阻力越大。

表1 油水兩相微觀滲流平面模型相關參數

圖3 模型A實驗與數值模擬的結果對比

圖4 模型B實驗與數值模擬的結果對比

記錄玻璃刻蝕實驗條件下各時刻的含水飽和度，與數值模擬得到的含水飽和度在相同歸一化時間下進行對比，發現二者的吻合度也很高(圖5)。左端見水時，模型A與模型B的含水飽和度分別為31.7%和66.8%，與圖5中曲線拐點處的含水飽和度吻合；見水后，采收率曲線變得平緩，上升幅度緩慢，水相突破后驅替阻力減小，水基本沿著主要流動通道運移，波及范圍幾乎不再擴大，驅油效率不再增加。模型A存在較明顯竄流通道，且原油黏度較高，指進現象嚴重，水沿竄流通道迅速突破，模擬得出的最終驅油效率僅為41.2%；模型B孔喉分布較模型A均勻，波及效率較高，模擬得出的最終驅油效率為72.8%。

圖5 實驗與數值模擬的含水飽和度對比

3.2 喉道內流體流動特征分析

對模型B的驅替過程分析表明，在保持驅替壓差不變的情況下，當含水飽和度達到72.8%時，流動達到穩定狀態，孔隙基本不出油。取9處富含殘余油的喉道，其孔隙及流體特性參數值如表2所示。喉道殘余油的形成主要是由于驅替相水在垂直于喉道界面的方向動力不足，對孔喉內流體的驅動速度遠小于其他喉道的平均流速，當小喉道被其它喉道率先突破的水相截斷時驅動速度基本為零，喉道殘余油形成。

由于所建立的孔隙模型壁面相對于實驗模型較為理想光滑，模擬得到的驅替速度大于實驗條件下的驅替速度。模型B在數值模擬中驅替壓力梯度為內部最大瞬時流速4.407 m/s，出口最大速度1.093 m/s，出口平均速度0.289 6 m/s。從模型B整個孔隙模型驅油效果來看，殘余油富集喉道內最大流速遠小于出口流速，流速隨喉道內驅替壓力梯度的增大而增大(圖6a)，但與喉道直徑無明顯關系(圖6b)，因此，僅從單個孔道大小無法衡量驅油阻力大小及分析殘余油的形成原因。

表2 模型B富含殘余油喉道特性參數

圖6 模型B喉道內流速與壓力梯度及孔徑的關系

連接2個孔隙的喉道有多個，且可視為不等徑并聯喉道單元(圖7)。各喉道內的驅替速度受驅油動力和油水流動阻力(即驅替壓差、油水黏度比及并聯各孔徑比值)的相互制約。在表2中，孔隙4和5、8和9屬于同一并聯喉道系統。在同一并聯喉道單元中，水相在流速較快的大喉道優先突破，在流速較慢的小喉道形成封堵，從而降低了喉道首尾壓差，導致殘余油的形成，這種現象簡稱為微觀喉道竄流。分析認為，并聯喉道直徑的比值對孔道是否發生微觀竄流影響很大。

圖7 模型 B中并聯喉道單元示意圖

為節省模型計算時間并保證數值計算的收斂性，將圖7所示的模型孔隙中各個并聯喉道單元中喉道直徑依次近似取值為10、20、30、40、50、80、100 μm，模擬計算得出驅替壓差為5 kPa時相應直徑喉道內的最大流速與直徑為10 μm的喉道內的流速比分別為1、3、53、154、235、500、660倍。由此可見，微觀孔隙結構中，孔道的直徑對孔道內兩相驅替流速的影響很大。

在水驅過程中，水相隨喉道直徑從大到小依次突破。在模擬驅替壓差為5 kPa時并聯喉道水驅過程中，當30 μm直徑喉道發生水相突破后，直徑為10 μm及以下的喉道油水相界面基本不再向前推移，此時發生微觀竄流的臨界喉道直徑比為3，臨界流速比為53。

在改變驅替壓差和油水黏度比的條件下對并聯喉道單元進行模擬計算，得到不同條件下微觀竄流臨界喉道直徑比(圖8)。由圖8a可以看出，當驅替壓差一定時，油水黏度比越大，發生竄流的臨界喉道直徑比越小，即油水黏度比越大，越易發生微觀竄流。由圖8b可以看出，在油水黏度比一定時，驅替壓差越大，發生竄流的臨界喉道直徑比越小，即驅替壓差越大，越易發生微觀竄流。

圖8 模型B中油水黏度比及驅替壓差對微觀竄流臨界喉道直徑比的影響

由以上分析可以發現，微觀條件與宏觀條件下的竄流規律吻合。在孔隙尺度微觀水驅油過程中，當驅替壓差恒定時，直徑較大的喉道中流體速度增加速率遠大于直徑較小的喉道，并且增加油水黏度比和驅替壓差均不利于增加直徑較小喉道中液體流速的相對增加量，因此小喉道更容易被優先突破大喉道的液體封堵，更易導致微觀竄流的發生，從而在相對較小的喉道中形成殘余油。

4 結論

1) 基于Navier-Stokes方程，結合Level set數學方法，進行了微觀二維孔隙尺度水驅油模擬，模擬結果與實驗得到的結果吻合度很高，為微觀水驅油機理研究提供了新的技術思路。

2) 考慮整個模型中孔道內驅替流體流速與孔徑相關性差，從不等徑并聯喉道單元的角度分析了微觀竄流臨界喉道直徑比，并對不同直徑喉道單元在不同驅替條件下的流動特征進行了模擬，結果表明驅替壓差及油水黏度比越大，微觀不等徑并聯喉道中發生竄流的可能性越大，更易導致在相對較小的喉道中形成微觀殘余油。

[1] 姜漢橋,宋亮,張賢松,等. 基于核磁共振的正韻律厚油層高含水期挖潛室內實驗[J].中國海上油氣，2014,26(6):40-43. Jiang Hanqiao,Song Liang,Zhang Xiansong,et al.Laboratory NMR experiments on tapping the production potential of positive rhythmic and thick oil reservoirs in high water-cut stage[J].China Offshore Oil and Gas,2014,26(6):40-43.

[2] FATT I.The network model of porous media I:capillary pressure characteristics[J].Trans AIME,1956,207:144-159.

[3] BLUNT M J,JACKSON M D,PIRI M,et al.Detailed physics,predictive capabilities and macroscopic consequences for pore-network models of multiphase flow[J].Advances in Water Resources,2002,25(8-12):1069-1089.

[4] LOVOLL G,MEHEUST Y,MALOY J K,et al.Competition of gravity,capillary and viscous forces during drainage in a two-dimensional porous medium:a pore scale study[J].Energy,2005,30(6):861-872.

[5] RAMSTAD T,IDOWU N,NARDI C,et al.Relative permeability calculations from two phase flow simulations directly on digital images of porous rocks[J].Transport in Porous Media,2012,94(2):487-504..

[6] RICHARD D.Oriented percolation in two dimensions[J].The Annuls Probability,1999 ,12(4):999-1040.

[7] 許友生,劉慈群,俞慧丹.多孔介質中兩相驅離的格子Boltzmann模型新研究[J].應用數學和力學,2002，23(4)： 353-358. Xu Yousheng,Liu Ciqun,Yu Huidan.New studying of Lattice-Boltzmann method for two-phase driven in porous media[J].Applied Mathematics and Mechanics,2002，23(4)： 353-358.

[8] 張超,張智,曾春珉,等.潿洲11-4油田含CO2氣井油管柱腐蝕分析[J].中國海上油氣,2015,27(4):122-125.DOI:10.11935/j.issn.1673-1506.2015.04.018. Zhang Chao,Zhang Zhi,Zeng Chunmin,et al.Analysis on tubing corrosion for gas wells with CO2in WZ 11-4 oilfield[J].China Offshore Oil and Gas，2015,27(4):122-125.DOI:10.11935/j.issn.1673-1506.2015.04.018.

[9] 張媛,王碩亮,張賢松.微圓管實驗壓力誤差校正方法探討[J].中國海上油氣,2016,28(2):94-98.DOI:10.11935/j.issn.1673-1506.2016.02.012. Zhang Yuan,Wang Shuoliang,Zhang Xiansong.Discussion on method for correcting the pressure error in micro tube experiment[J].China Offshore Oil and Gas,2016,28(2):94-98.DOI:10.11935/j.issn.1673-1506.2016.02.012.

[10] OSHER S,SETHIAN J A.Fronts propagating with curvature-dependent speed:algorithms based on Hamilton-Jacobi formulations[J].Journal of Computational Physics,1988,79(1):12-49.

[11] SUSSMAN M,ALMGREN A S,BELL J B,et al.An adaptive level set approach for incompressible two-phase flows[J].Journal of Computational Physics,1999,148(1):81-124.

[12] OLSSON E,KREISS G.A conservative level set method for two phase flow[J].Journal of Computational Physics,2005,210(1):225-246.

[13] OLSSON E,KREISS G,ZAHEDI S.A conservative level set method for two phase flow II[J].Journal of Computational Physics,2007,225(1):785-807.

[14] GUNDE A,BERA B,MITRA S K.Investigation of water and CO2(carbon dioxide) flooding using micro-CT (micro-computed tomography) images of Berea sandstone core using finite element simulations[J].Energy,2010,35(12):5209-5216.

[15] ZHU Qianlin,ZHOU Qianlong,LI Xiaochun.Numerical simulation of displacement characteristics of CO2injected in pore-scale porous media[J].Journal of Rock Mechanics and Geotechnical Engineering,2016,8(1):87-92.

[16] COMSOL Inc.COMSOL multiphysics[R].Stockholm:COMSOL Inc,2014.

[17] SHEPEL S V,SMITH B L.New finite-element/finite-volume Level set formulation for modelling two-phase incompressible flows[J].Journal of Computational Physics,2006,218(2)：479-494.

(編輯：張喜林)

Simulation analysis of microscopic water-oil displacement based on Level set method

Gao Yajun1,2Jiang Hanqiao1Wang Shuoliang2Li Junjian1Liu Chuanbin1Chang Yuanhao1Wang Yicheng1

(1.MOEKeyLaboratoryofPetroleumEngineering,ChinaUniversityofPetroleum,Beijing102249,China; 2.SchoolofEnergyResources,ChinaUniversityofGeosciences,Beijing100083,China)

A new microscopic two-phase seepage flow mathematic model is established and solved with finite element method, in which the movement characteristics of the two-phase interface are characterized with Level set method and Navier-Stokes equation. The numerical simulation is verified by the experiment conducted with micro etching glass model. Results show that the numerical simulation results agree well with that of physical experiment. Navier-Stokes equation combined with Level set method implements the micro numerical simulation of two-phase flow, which provides a new method to study micro two-phase seepage characteristics and oil displacement mechanisms. The simulation is carried out under different displacement conditions based on the unequal-diameter parallel pore throat unit model, and the variation characteristics of critical throat diameter ratio under pore-scale two-phase channeling conditions are analyzed. Results show that the larger the displacement pressure and viscosity ratio of the oil-water two-phase are, the easier the micro-channeling occurs, which leads to the formation of microscopic residual oil in the smaller throat. The work lays a foundation for further study of microscopic oil displacement mechanisms.

numerical simulation; Level set method; microscopic displacement; microscopic residual oil

1673-1506(2016)06-0059-07

10.11935/j.issn.1673-1506.2016.06.010

高亞軍，男，在讀碩士研究生，從事油氣田開發研究。地址：北京市昌平區府學路18號中國石油大學(北京)研修樓A座401室(郵編：102249)。E-mail：gaoyajuncup@163.com。

王碩亮，男，講師，2011年畢業于中國石油大學(北京)，獲博士學位，主要從事油氣田開發方面的研究。地址：北京市海淀區學院路29號(郵編：100083)。E-mail:wangshuoliang@cugb.edu.cn。

TE312

A

2016-06-28 改回日期：2016-07-22

*國家重點基礎研究發展計劃“致密油高效開發油藏工程理論與方法研究(編號：2015CB250905)”部分研究成果。