基于Level Set有限元方法的微觀水驅(qū)油數(shù)值模擬

高亞軍，姜漢橋，王碩亮，劉傳斌，常元昊,王依誠

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

?

基于Level Set有限元方法的微觀水驅(qū)油數(shù)值模擬

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

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

在微觀水驅(qū)油實驗的基礎上，建立了二維微觀孔隙模型；引入Level Set數(shù)學方法，結(jié)合N-S方程,建立了微觀兩相滲流數(shù)學模型，借助有限元方法，進行水驅(qū)油兩相數(shù)值模擬，研究微觀水驅(qū)油動態(tài)特征。研究了孔喉非均質(zhì)性、潤濕非均質(zhì)性、油水兩相黏度以及驅(qū)替壓差對驅(qū)油效果的影響；通過對模擬結(jié)果中孔隙內(nèi)流體的波及效率、兩相界面的運移速度以及優(yōu)勢滲流通道的分析發(fā)現(xiàn)，Level Set方法能很好地處理各因素下油水兩相界面的拓撲變化，為微觀水驅(qū)油機理的研究提供了一種新的技術方法。

水驅(qū)實驗；數(shù)值模擬；微觀驅(qū)替特征

許多學者在微觀水驅(qū)油機理和兩相流動特征的研究方面做了大量的實驗，實驗方法主要有巖心薄片平面驅(qū)替和玻璃刻蝕物理模擬等二維方法[1-3]，以及驅(qū)替巖心CT掃描技術的三維方法[4-5]，但以上方法實驗模擬仿真成本高，操作較繁瑣和困難,并且實驗監(jiān)測技術和圖像獲取精度有限，導致實驗數(shù)據(jù)誤差較大。數(shù)值模擬方法不需要實際巖心參與，利用計算機對微觀巖心孔隙模型進行刻畫，模擬實際水驅(qū)油整個過程的演化規(guī)律，研究油水兩相動態(tài)變化特征，具有便捷和可視化等優(yōu)點。目前應用較多的是孔隙網(wǎng)絡模型[6-7]，但此方法一方面無法更進一步清晰呈現(xiàn)孔道內(nèi)流體兩相界面運移狀態(tài)，另一方面無法直觀獲取各流體質(zhì)點在任何時刻的壓力和流速變化。

針對以上方法所存在的問題，本文提出基于Level Set(水平集)方法的微觀水驅(qū)油兩相數(shù)值模擬，對孔隙內(nèi)流動特征進行可視化研究。水平集方法是由Osher和Sethian[8]提出的，最初主要應用于智能控制、圖像處理等方面，近年來經(jīng)國外學者對求解方法的不斷改進[9-10]，逐漸被用于氣液或液液兩相流動的數(shù)值模擬研究領域中。有學者用此方法做過簡單模型的油水兩相模擬[11-12]，但局限于理想孔隙結(jié)構(gòu)或油氣兩相驅(qū)替，均未與實驗進行對比驗證，且未涉及到微觀孔隙內(nèi)的油水兩相模擬。為此筆者根據(jù)巖心CT掃描圖像，進一步處理后，建立二維微觀孔隙模型，結(jié)合N-S方程，引入Level Set方法來對油水兩相界面的動態(tài)運移進行追蹤，并用有限元方法求解。

1　油水兩相驅(qū)數(shù)學模型

1.1Level Set方法

假設地層中巖石孔隙內(nèi)水油兩相流體均為不可壓縮流體，流動狀態(tài)為層流。在層流兩相流中，水平集界面自動生成界面運移的方程，在Level Set方程中，當φ=0時的水平集輪廓線代表兩相流體界面。對于此油水兩相，在油中φ=0，水中φ=1。Level Set方程可以被看作是表示水油兩相流中水的體積分數(shù)，因此，兩相流體的界面運移方程可以用下式表示：

(1)

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

借助于Level Set函數(shù) ，流體的物性可用下述方程表示，即：

(2)

式中：ρo、ρw分別為油和水的密度；μo、μw分別為油、水的黏度。經(jīng)過這樣定義后，流體的密度和黏度就可以在兩相接觸面處平滑過渡，不會形成突變，保證計算的穩(wěn)定收斂。

1.2控制方程

(3)

式中：ρ為密度，kg/cm3；μ為動態(tài)黏度，Pa·s；u為速度，m/s；I為單位矩陣；p為入口壓力，Pa；Fst為油水界面表面張力，N/m 。

表面張力項可以通過下式計算[13]：

Fst=·T

(4)

T=σ(I-(nnT))δ

(5)

式中：n為界面法向量；σ為表面張力系數(shù)；δ表示狄克拉函數(shù)，在兩相界面以外區(qū)域此函數(shù)值均為零。函數(shù)被近似為：

(6)

界面法向量可表示為：

(7)

當用有限元方法求解方程時，先使方程變形，然后在計算區(qū)域內(nèi)分部積分。層流兩相流中，為了表示壁面與流體的相互作用力，在計算域上求積分可以用加上一個邊界積分的形式表示如下：

(8)

式中：θ為壁面接觸角，當使用無滑移邊界條件即非潤濕條件時，則無此邊界項，因為在邊界上時test(u)=0，不能設置接觸角。然而，如果允許流體細微的邊界滑移(可以忽略不計)，就可以設置接觸角。上式方程中給出了潤濕壁面所增加的邊界條件項，因此可以設置兩相的壁面接觸角。

1.3邊界條件

入口、出口邊界上的壓力由模型設置，兩側(cè)邊界為壁面封閉邊界，初始界面設置為t=0時水相和油相的接觸面，壁面條件為可潤濕壁面，垂直壁面方向流速為零。即：

(1)入口邊界：p=p0,n·μ(u+(u)T)=0

(2)出口邊界：p=p1,n·μ(u+(u)T)=0

式中：β是滑移長度，其中，β值的大小為邊界壁面部位剖分網(wǎng)格的尺寸大小；Frf為壁面作用力。

2　水驅(qū)油仿真模擬

油水兩相驅(qū)仿真模擬研究表明，構(gòu)建仿真模型是研究驅(qū)替機理的重要環(huán)節(jié)[13-14]。本文用CT掃描巖心切片處理后建立的玻璃刻蝕模型進行仿真模擬。模型的構(gòu)建需要經(jīng)過多次反復試驗、驗證、修繕，才能最終建立可行的模型。本文水驅(qū)油特征動態(tài)模擬的基本步驟如下：

(1)對實際儲層巖心進行CT掃描，經(jīng)圖像處理后，建立微觀孔喉模型。

(2)根據(jù)流體運動方程，結(jié)合Level Set數(shù)學方法，建立數(shù)學模型。

(3)對實際物理模型的水驅(qū)油實驗中壓力、黏滯力、溫度、潤濕性等因素所反映出的兩相動態(tài)特征進行分析，建立數(shù)學模型的邊界控制條件。

(4)改變仿真模擬中兩相流體參數(shù)以及邊界控制條件，根據(jù)整個仿真圖像中油水運移狀態(tài)，研究不同條件下水驅(qū)機理及剩余油的控制因素。

根據(jù)上述步驟，首先需要對實際所建立的微觀孔喉模型進行水驅(qū)實驗，為此筆者選取根據(jù)實際巖心CT掃描孔喉特征所建立的微觀非均質(zhì)性不同的玻璃刻蝕模型A和B進行驅(qū)替實驗，對所建立的微觀兩相驅(qū)仿真模型進行對比驗證。

2.1水驅(qū)油實驗

實驗中采用光化學刻蝕的玻璃模型A和B，模型尺寸分別為17.45 mm×14.31 mm和15.28 mm×12.13 mm。為了能夠較好地體現(xiàn)出油水兩相在驅(qū)替過程中相界面的動態(tài)變化特征，實驗所用原油黏度較大，飽和原油黏度分別為50 mPa·s和20 mPa·s，孔隙度分別為為30.96%和28.91%。實驗用鹽水密度為1.05 g/cm3，黏度為1 mPa·s。實驗步驟如下：

(1)安裝固定玻璃刻蝕模型，將模型抽空飽和油，然后再進行水驅(qū)油；

(2)按實驗方案中壓力進行水驅(qū)油實驗，用顯微放大系統(tǒng)和錄像機錄取實驗過程的圖像，并用計算機圖象處理系統(tǒng)對錄像和圖片進行處理；

(3)對整個驅(qū)替實驗的錄像進行分析，統(tǒng)計不同時刻的驅(qū)油效率，分析微觀水驅(qū)相界面特征變化并解釋剩余油的形成。

2.2數(shù)值模擬與實驗的對比分析

針對玻璃刻蝕模型，將CT掃描的孔隙結(jié)構(gòu)經(jīng)過二值化轉(zhuǎn)換處理可得到二維孔隙輪廓模型，對孔隙區(qū)域進行自由三角網(wǎng)格剖分，建立相關有限元模型進行求解。實驗和仿真模擬均從右端注入，左端采出，各參數(shù)如表1所示。

表1　數(shù)值模擬模型參數(shù)設置

數(shù)值模擬條件較為理想，而實驗條件下影響因素較多，相同含水飽和度下，實際驅(qū)替時間與數(shù)值模擬的驅(qū)替時間不相等。因此，將各階段記錄時間與數(shù)值模擬的時間在(0,1)作歸一化處理進行對比分析。記錄玻璃刻蝕實驗條件下各時刻的的含水飽和度，與數(shù)值模擬得到的含水飽和度在相同歸一化時間下進行對比，發(fā)現(xiàn)吻合度很高(圖1)。剛注入時刻，由于數(shù)值模擬不存在實驗條件下注入速度的穩(wěn)定過程所帶來的誤差，因此初始時刻的擬合存在一定的相對誤差(27.2%)，隨著驅(qū)替的進行，相對誤差越來越小(5.5%)。

圖1　實驗與數(shù)值模擬的含水飽和度對比

將驅(qū)替實驗A和B所錄取的圖像與數(shù)值模擬結(jié)果進行等含水飽和度對比發(fā)現(xiàn)(圖2)，驅(qū)替前緣兩相界面吻合度很高，其中，紅色為油相，藍色為水相。A左端見水時的含水飽和度為31.7%，與圖1中模型A曲線拐點處的含水飽和度非常接近。相對于模型B，模型A的非均質(zhì)性較嚴重，左端見水后，水相突破后流速通道阻力減小，注入水很難波及到其它孔隙，驅(qū)油效率降低，最終采收率為41.2%。模型B左端初始見水的含水飽和度為66.8%(圖2)，同樣與圖1中曲線拐點基本一致，見水后驅(qū)油效率基本不再增加，整個驅(qū)替過程中，油水兩相界面推進較均勻，波及范圍大，最終采收率達到72.8%。由于模型A的原油黏度遠大于B，因此驅(qū)替過程中，水相在進入細小孔喉時的阻力非常大，模型A的指進現(xiàn)象較為嚴重。

在與實驗的驅(qū)替過程對比之下，發(fā)現(xiàn)本文所建

圖2　實驗與數(shù)值模擬的結(jié)果對比

立的數(shù)學有限元仿真模擬方法能夠較好地刻畫水驅(qū)油整體過程的相界面動態(tài)分布，因此可以對其它特征進行進一步研究。

3　微觀水驅(qū)油特征

3.1微觀孔喉非均質(zhì)性

由于孔隙介質(zhì)存在微觀非均質(zhì)性，不同半徑孔喉的毛管力差異較大，驅(qū)替相所受阻力也不同，此時繞流和指進為主要的微觀驅(qū)油機理。對整個二維平面模型來看，水相由入口端進入，沿著一條或數(shù)條阻力較小、孔徑相對較大的孔道向前突進并首先到達出口，即微觀指進現(xiàn)象。從實驗和仿真模擬中均可以看到當，當模型出口端一旦見水，注入水大部分便會只沿著突破通道流動，模型A從左端初期見水時采收率31.7%到最終穩(wěn)定采收率41.2%，驅(qū)油效率增加幅度很小(圖1)。模型B的注入水推進較為均勻，指進現(xiàn)象較弱，水驅(qū)效果較好(圖2b)。對比模型A和B，可以發(fā)現(xiàn)A的孔喉非均質(zhì)性比B嚴重，模型A中部存在大孔道，模型B孔喉分布相對較為均質(zhì)，且B的孔喉比A小，雖然驅(qū)替需要更大的壓力，但是壓力在平行驅(qū)替方向孔喉內(nèi)分布較均勻，油水前緣均勻推進，波及效率較高，驅(qū)油效果較好。

3.2潤濕非均質(zhì)性

孔隙內(nèi)油水兩相的流動狀態(tài)與壁面潤濕性緊密相關。為了體現(xiàn)潤濕性對水驅(qū)油時流動孔道的選擇特點，在本文所建立的仿真模型中，以壁面接觸角表示潤濕角。設置部分孔隙水濕接觸角為60°(圖3，b1)，其余條件不變，在相同驅(qū)替時間下，將水濕條件下的模擬結(jié)果與孔隙壁面為中性潤濕條件下的模擬結(jié)果進行對比(圖3)，以同一個壓力系統(tǒng)相同驅(qū)替條件的a2和b2區(qū)域為參考基準，發(fā)現(xiàn)在水濕條件下，孔隙中油水兩相界面運移更快(圖3，b1)。部分壁面水濕加快了水相的指進，導致大孔隙繞流現(xiàn)象更嚴重；大孔隙油水界面的迅速推移，導致小孔隙內(nèi)憋壓不足，油水界面停止向前推移(圖3，b2)。經(jīng)過以上分析，可以發(fā)現(xiàn)本文所建立的微觀仿真數(shù)學模型能夠較好地處理壁面潤濕性所引起的油水兩相運移動態(tài)變化特征。

3.3油水黏度比

原油的黏度在流動過程中反映為原油內(nèi)部的摩擦阻力，決定兩相與壁面的摩擦力大小，從而影響兩相驅(qū)替的流動速度和驅(qū)油特征。取同一驅(qū)替時刻不同黏度下模擬結(jié)果(圖4)可見，在相同的驅(qū)替時間下，當油水黏度比較小時，驅(qū)替阻力更小，兩相界面運移更快。油水黏度比越大，水相進入小孔隙的阻力越大，水沿著大孔道的指進現(xiàn)象越來越明顯。由圖4可看出，level set方法能夠較好地處理油水兩相黏度變化所引起的相界面的拓撲變化。

圖3　模型A不同潤濕角下模擬結(jié)果

3.4驅(qū)替壓力

高壓驅(qū)替時，水沿大孔道中心迅速突破，指進明顯，截斷細小孔道，使油相從連續(xù)相變?yōu)榉沁B續(xù)相，形成殘余油，但高壓驅(qū)替可將部分可動剩余油驅(qū)替出來(圖5和圖6)。選取其中典型局部區(qū)域放大進行詳細分析如下：

(1)在油水前緣到達同等距離時，對比區(qū)域a和區(qū)域b中驅(qū)替現(xiàn)象(圖5和圖6)，發(fā)現(xiàn)高壓力驅(qū)替下的區(qū)域b中，油的體積分數(shù)達到50%左右，油水過渡帶很長，指進現(xiàn)象明顯，驅(qū)替效率明顯低于壓力較小時的驅(qū)替結(jié)果。

圖4　模型A不同油相黏度下模擬結(jié)果

圖5　P=50 000 Pa時模型模擬結(jié)果

圖6　P=100 000 Pa時模型B模擬結(jié)果

(2)對比圖中c、d區(qū)域，發(fā)現(xiàn)驅(qū)替壓差過大(d區(qū)域)，即水相推進速度過快時，水主要沿著大孔道中心流動，繞過連通性不好、孔喉非均質(zhì)嚴重的區(qū)域，將孔隙角落的原油封堵死，形成角狀殘余油。

(3)對于一些小孔道(e區(qū)域)，毛管力較大，當驅(qū)替壓力小于原油的流動阻力時，不能驅(qū)替小孔喉內(nèi)的原油，形成可動殘余油，殘余油主要以塞狀或柱狀形式存在。

3.5流速可視化分析

在數(shù)值模擬流速示意圖中，可以清楚地看出存在主要流速通道(圖7)。對整個驅(qū)替過程中的壓力和流速分析表明，水通過優(yōu)勢通道到達出口端時，整個流域的流動阻力降低，其它小孔隙壓力相對降低，波及效率降低，有利于剩余油的形成。流體在孔隙通道中心的流速最大，越靠近孔隙壁面，流速越小。那些微小孔隙、封閉空隙中，油水幾乎不參與流動，流速較小，孔隙中壓力降落較快，坡度較陡；而對于那些大孔隙高流速的空隙中，壓力在流動方向上降落較慢。

圖7　模型B流速示意圖

沿y軸方向?qū)ψ⑷攵说?個孔隙編號依次為y1，y2，y3，y4，y5，y6，y7，統(tǒng)計同一驅(qū)替時間步長下入口端的驅(qū)替流速，做出入口端整個縱剖面在固定壓力下的注入速度，可以清楚看到在7個孔隙入口處流速大小(Vy2>Vy5>Vy4>Vy7>Vy1>Vy6>Vy3)，入口速度與圖7中的主要流速通道的分布相對應。但由于模型的絕對理想化，在與實際驅(qū)替規(guī)律不變的情況下，模擬驅(qū)替流速大于實際的驅(qū)替流速。在孔隙y5入口中可以發(fā)現(xiàn)，孔喉中心流速最大，越靠近孔隙壁面越小(圖8)。其中在孔道中心處可以看到流速并非最大，呈下凹現(xiàn)象，主要是因為y5孔喉入口內(nèi)部中心巖石顆粒阻擋形成分流所致。在驅(qū)替壓力梯度不變的情況下，隨著驅(qū)替時間的進行，孔隙中流速逐漸增大。主要是因為在微尺度驅(qū)替下，慣性力可以忽略不計，黏滯阻力占主導地位，當?shù)宛ざ纫后w驅(qū)替高黏度液體時，黏滯阻力逐漸減小[15]，因此驅(qū)替速率逐漸增大。

圖8　單孔隙入口流速隨時間變化規(guī)律

4　結(jié)論

(1)引入Level Set數(shù)學方法，進行微觀二維平面水驅(qū)油仿真模擬，數(shù)值模擬的結(jié)果與實驗得到的結(jié)果吻合度很高，為微觀水驅(qū)油的特征研究提供了新的可視化技術手段。

(2)Level Set方法在處理油水兩相黏度、潤濕性以及驅(qū)替壓力等引起的相界面拓撲變化時，具有較高的仿真度，可用此方法來進行更為復雜的微觀流體流動機理研究。

(3)Level Set數(shù)值模擬方法具備常規(guī)物理仿真和數(shù)值仿真所不具備的優(yōu)點，能夠清楚研究孔隙內(nèi)任何部位流體的壓力、流速等變化特征。油水兩相流在孔道中心流速最大，越靠近孔隙壁面流速越小；孔徑越小，壓力降落越快；在驅(qū)替壓力梯度和油水黏度比恒定時(油相黏度大于水相)，孔道內(nèi)油水兩相驅(qū)替速度逐漸加快。

[1]高輝，孫衛(wèi)，路勇，等.低滲透砂巖儲層油水微觀滲流通道與驅(qū)替特征實驗研究——以鄂爾多斯盆地延長組為例[J].油氣地質(zhì)與采收率，2011，18(1)：58-62.

[2]李洪璽，劉全穩(wěn)，何家雄，等.物理模擬研究微觀剩余油分布[J].新疆石油地質(zhì)，2006，27(3)：351-353.

[3]劉太勛，徐懷民.扇三角洲儲層微觀剩余油分布模擬試驗[J].中國石油大學學報：自然科學版，2011，35(4)：20-26.

[4]呂偉峰，冷振鵬，張祖波，等.應用CT掃描技術研究低滲透巖心水驅(qū)油機理[J].油氣地質(zhì)與采收率，2013，20(2)：87-90.

[5]曹永娜.CT掃描技術在微觀驅(qū)替實驗及剩余油分析中的應用[J].CT理論與應用研究，2015，24(1)：47-56.

[6]高慧梅，姜漢橋，陳民鋒，等.儲集層微觀參數(shù)對油水相對滲透率影響的微觀模擬研究[J].石油勘探與開發(fā)，2006，33(6)：734-737.

[7]陳民鋒，姜漢橋.基于孔隙網(wǎng)絡模型的微觀水驅(qū)油驅(qū)替特征變化規(guī)律研究[J].石油天然氣學報，2006，28(5)：91-95.

[8]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.

[9]Sussman M.A sharp interface method for incompressible two-phase flows[J].Journal of Computational Physics, 2007,221(2):469-505.

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

[11]王琳琳，田輝，李國君.基于Level Set方法對油水和氣水兩相界面的數(shù)值模擬[J].應用力學學報，2010，27(2)：298-302.

[12]Qianlin Zhu,Qianlong Zhou,Xiaochun Li.Numerical simulation of displacement characteristics of CO2injected in pore-scale porous media，2016,81:87-92.

[13]孫煥泉，孫國，程會明，等.勝坨油田特高含水期剩余油分布仿真模型[J].石油勘探與開發(fā)，2002，29(3)：66-68.

[14]徐守余，朱連章，王德軍.微觀剩余油動態(tài)演化仿真模型研究[J].石油學報，2005，26(2)：69-72.

[15]A.R Ershov,Z M Zorin,V D Sobolev, et al.Displacement of silicone oil from the hydrophobic surface by aqueous trisiloxane solutions[J].Colloid Journal，2001,63:290-295.

編輯：李金華

1673-8217(2016)05-0091-06

2016-05-01

高亞軍，1991年生，2015年畢業(yè)于中國地質(zhì)大學(北京)石油工程專業(yè)，在讀碩士研究生，從事油氣田開發(fā)方向研究。基金項目：國家重點技術研究發(fā)展計劃(973計劃)項目“致密油高效開發(fā)油藏工程理論與方法研究”(2015CB250905)。

TE311

A