張晟庭 * 李 靖 *, 陳掌星 *, 張 濤 ** 吳克柳 * 馮 東 * 畢劍飛 * 李相方 *

* (中國石油大學(xué) (北京) 油氣資源與探測國家重點實驗室,北京 102249)

? (卡爾加里大學(xué)化學(xué)與石油工程系,加拿大阿爾伯塔 T2N1N4)

** (西南石油大學(xué)油氣藏地質(zhì)及開發(fā)工程國家重點實驗室,成都 610500)

引言

復(fù)雜多孔介質(zhì)中的非混相驅(qū)替在兩相流動中占據(jù)重要地位,并在氣驅(qū),氣水交替及泡沫驅(qū)等提高油氣采收率領(lǐng)域具有重要研究意義[1-5].在氣驅(qū)提高油氣采收率過程中,連續(xù)氣相注入儲層后,將會形成復(fù)雜的氣液非混相驅(qū)替模式;注入儲層中的氣相一般為非潤濕相,在非潤濕相驅(qū)替潤濕相的過程中,孔喉壁面上將會滯留一層液膜,液膜的存在可能會造成驅(qū)替過程中發(fā)生卡斷甚至水鎖現(xiàn)象.因此,在孔隙尺度研究氣液非混相驅(qū)替過程中的卡斷機理有利于人們認識微觀孔喉中的氣液流動狀態(tài).

研究人員通過理論研究,實驗及數(shù)值模擬等方法,對氣液兩相驅(qū)替過程中的卡斷機理展開了大量研究工作.在理論方面,Roof[6]基于孔喉毛管壓力不平衡這一機制,建立了預(yù)測卡斷模型的靜態(tài)準(zhǔn)則,認為當(dāng)最小喉道半徑小于孔隙半徑的一半時,潤濕相將沿著孔隙壁面發(fā)生回流,也即導(dǎo)致非潤濕相發(fā)生卡斷.Gauglitz 等[7]給出了收縮圓柱形孔喉系統(tǒng)的卡斷準(zhǔn)則,認為毛管數(shù)Ca對卡斷發(fā)生的時間有一定影響,并且對于不同寬度比的孔喉系統(tǒng)影響規(guī)律不同.Ransohoff 等[8]給出了方形喉道截面的以及三角形喉道截面的靜態(tài)卡斷準(zhǔn)則,并認為在非圓形截面的孔喉中更容易發(fā)生卡斷現(xiàn)象.然而,上述卡斷準(zhǔn)則并沒有考慮動態(tài)參數(shù)對卡斷影響,Tsai 和Miksis[9]通過數(shù)值方法發(fā)現(xiàn)毛管數(shù)過大或者過小時,氣泡不會發(fā)生卡斷現(xiàn)象.Deng 等[10-11]擴展了圓形及非圓形縮頸孔喉系統(tǒng)的靜態(tài)卡斷準(zhǔn)則,發(fā)現(xiàn)之前的靜態(tài)準(zhǔn)則對卡斷的預(yù)測過于保守,同時發(fā)現(xiàn)即使喉道系統(tǒng)達到了卡斷發(fā)生的靜態(tài)準(zhǔn)則,高毛管數(shù)也會抑制卡斷的發(fā)生.在實驗方面,隨著微流控模型的發(fā)展,大量學(xué)者采用微流控模型對兩相流動的卡斷現(xiàn)象進行可視化實驗,Tian 等[12]利用2D 微模型研究了氣水流動過程中的卡斷現(xiàn)象,并驗證了Roof 等提出的靜態(tài)準(zhǔn)則.Cha 等[13]利用2D 矩形截面微模型研究了氣泡卡斷機理,并與提出的卡斷理論模型得到了較好的一致性.Tetteh 等[14]利用微流控模型觀察到了低礦化度水驅(qū)過程中的油相卡斷.Wu 等[15]利用實驗研究非縮頸孔喉的氣泡卡斷現(xiàn)象,并且認為喉道長度及喉道寬度對卡斷有一定影響.

相比于實驗,數(shù)值模擬方法可以更好的幫助理解及預(yù)測孔隙尺度氣液兩相流動狀態(tài)及卡斷現(xiàn)象.在孔隙尺度上模擬多相流的方法有孔隙網(wǎng)絡(luò)模型(PNM)[16]、光滑粒子動力學(xué)方法(SPH)[17]、密度泛函流體動力學(xué)方法(DFH)[18]、流體體積方法(VOF)[19]以及格子玻爾茲曼方法[20]等.部分學(xué)者采用VOF方法研究了縮頸孔喉系統(tǒng)中的卡斷現(xiàn)象,Raeini 等[21]利用VOF 方法研究了星形喉道截面的卡斷現(xiàn)象,并考慮了動態(tài)因素的影響.Starnoni 和Pokrajac[22]采用VOF 方法考慮接觸角及黏度比對卡斷影響,研究表明,喉道截面圓度越高發(fā)生卡斷的臨界接觸角越小,黏度比的增大也會降低發(fā)生卡斷的臨界接觸角.Zhang 等[23]采用VOF 方法研究了方形截面喉道的卡斷現(xiàn)象,并認為高毛管數(shù)及高黏度比可以有效抑制卡斷的發(fā)生,Cha 等[13]利用VOF 方法模擬了矩形截面的卡斷現(xiàn)象,并給出了卡斷發(fā)生的孔喉寬度比條件.然而,VOF 方法很難準(zhǔn)確計算界面的法向和曲率,且界面重構(gòu)方法復(fù)雜,向高維推廣困難[24].格子玻爾茲曼方法是介于微觀分子動力學(xué)模擬方法和宏觀傳統(tǒng)數(shù)值模擬方法之間的一種介觀數(shù)值模擬方法,與VOF 方法等界面跟蹤或界面捕捉方法相比,格子玻爾茲曼方法中的氣液界面可以自動產(chǎn)生,演化與遷移,無須采用任何界面跟蹤或界面捕捉技術(shù)[25],并能直接處理流體與流體間的作用力及流體與固壁間的作用力[26].目前多相格子玻爾茲曼模型可以分為四類,包括顏色梯度模型[27]、偽勢模型[28]、自由能模型[29]以及相場模型[30].張磊等[31]基于顏色梯度模型在孔喉系統(tǒng)中模擬非混相驅(qū)替過程,并認為多孔介質(zhì)的非均質(zhì)性越強,流體越容易發(fā)生卡斷.趙玉龍等[32]基于顏色梯度模型,模擬地層高溫高壓條件下致密氣驅(qū)替地層水的流動過程,并發(fā)現(xiàn)在多孔介質(zhì)中大量連通的微小通道被卡斷水占據(jù).Alpak等[33]采用相場模型模擬3D 縮頸孔隙中的兩相驅(qū)替過程,觀察到了卡斷現(xiàn)象并與Roof[6]靜態(tài)準(zhǔn)則保持較好的一致性,但是模擬兩相流體的黏度相差不大,不能反應(yīng)氣液兩相的黏度比.Wei 等[34]采用偽勢模型模擬單個氣泡通過喉道過程,并觀察到了卡斷現(xiàn)象.Zhang等[35]采用改進的雙組分偽勢模型模擬氣驅(qū)水過程,成功在喉道壁面捕捉到了水膜的存在,并認為微納尺度氣水兩相流動過程中壁面液膜的厚度不可忽略.

但是,目前對卡斷現(xiàn)象的研究往往建立在縮頸漸變孔喉系統(tǒng)中,弱化了壁面液膜厚度的影響,同時宏觀VOF 方法不能充分表征壁面和液膜相互作用.為此,本文在原始的偽勢格子玻爾茲曼方法的基礎(chǔ)上,改進了流體之間作用力格式并添加流固作用力以捕捉氣液兩相驅(qū)替過程中喉道壁面上的液膜變化規(guī)律及其對氣液兩相流動狀態(tài)的影響.并進一步研究了驅(qū)替壓差(毛管數(shù)Ca),喉道長度及寬度等因素對液膜厚度及流動狀態(tài)的影響規(guī)律.本文的研究內(nèi)容為表征復(fù)雜多孔介質(zhì)的卡斷機理提供了理論及模擬基礎(chǔ).

1 單組分多相偽勢模型

1.1 經(jīng)典BGK 碰撞模型

在格子玻爾茲曼方法中,采用離散的密度分布函數(shù)表征流體的運動規(guī)律,標(biāo)準(zhǔn)BGK 碰撞過程可以表示為[36]

式中,fi(x,t)為i方向密度分布函數(shù),Δx和Δt分別表示網(wǎng)格步長和時間步長,τ表示松弛時間,其取值與流體的運動黏度ν相關(guān),表示為

ei為各方向的離散速度,對于D2Q9 模型表示為[37]

1.2 流體之間作用力及流固作用力

原始偽勢模型中流體之間的作用力定義為[28]

式中,G為格林函數(shù),ψ(x,t)為有效質(zhì)量,其定義為ψ(ρ)=ρ0[1-exp(-ρ/ρ0)],ρ0=1 為參考密度.進一步,考慮第一層粒子之間的作用力,上述作用力格式可以簡化為

式中,g為流體之間作用強度的控制參數(shù),w(|ei|2)為權(quán)重系數(shù),|ei|2=1,w(|ei|2)=1/3;|ei|2=2,w(|ei|2)=1/12.采用原始的偽勢模型的作用力格式會在兩相界面處存在較大的虛假速度.因此,本文采用Gong 等[38]提出的作用力改進格式以降低虛假速度,具體表示為

進一步考慮最近及次近的節(jié)點之間作用力,Chen 等[39]提出了上述作用力的簡化格式

式中,β為調(diào)節(jié)參數(shù),通過改變β大小可以改善模型的熱力學(xué)一致性并降低虛假速度.若僅考慮第一層節(jié)點的作用力,可將式(9)進一步擴展為[40]

式中,γ1和γ2為作用力格式的離散系數(shù),其中γ1=1/3,γ2=1/12.上述離散作用力格式也稱為E4 作用力格式.讀者可以參考文獻[40]的研究,將作用力擴展至E6 或E8 作用力格式獲得更低的虛假速度及更精確的界面張力.

流體與壁面的作用力類似于流體之間的作用力格式,其定義如下[41]

式中,s(x+eiΔt,t)為判斷流體節(jié)點和壁面的標(biāo)識函數(shù),如果x+eiΔt是固體節(jié)點其值為1,x+eiΔt是流體節(jié)點則其值為0.ρw為壁面密度,通過調(diào)節(jié)ρw的大小可以模擬不同的接觸角.

1.3 狀態(tài)方程

原始的偽勢模型的狀態(tài)方程可以通過泰勒展開獲得[28],具體表示為

當(dāng)采用上述的狀態(tài)方程時,會導(dǎo)致嚴(yán)重的熱力學(xué)不一致性,同時模擬的最大密度比和黏度比均在10 以內(nèi).文獻[42]提出將真實的流體狀態(tài)方程與偽勢模型結(jié)合,大大提高了偽勢模型的熱力學(xué)一致性.本文采用常用的Redlich-Kwong (RK)狀態(tài)方程表征氣液兩相,RK 狀態(tài)方程的定義為

式中,pEOS為實際流體的壓力,R為氣體常數(shù),T為系統(tǒng)的溫度,ρ為實際流體的密度,a,b為臨界參數(shù),其中a=0.427 48R2Tc2.5/pc,b=0.086 62RTc/pc.Huang等[43]給出了偽勢模型中的參數(shù)取值,a=2/49,b=2/21,R=1,臨界溫度Tc=0.196 1,臨界壓力pc=0.178 4,臨界密度ρc=2.988 7.在添加RK 狀態(tài)方程后可以采用下述方法計算偽勢

值得注意的是,當(dāng)采用上述格式計算偽勢時,g將會在計算壓力p的時候消掉,不再具有實際物理意義,僅作為保證根號內(nèi)的符號為正的作用,故本文的模擬中取g=-1.

1.4 外力項處理格式

在添加了流體-流體作用力及流固作用力后,本文采用Kupershtokh 等[44]提出的精確差分方法(EDM)格式處理外力項,其定義如下

式中,Ftotal表示流體之間作用力與流固作用力之和,EDM 格式能夠適應(yīng)較寬的溫度范圍,是目前最常見的外力處理格式之一.

1.5 邊界條件

格子玻爾茲曼方法中的邊界條件對數(shù)值計算精度,計算穩(wěn)定性以及計算效率具有很大影響[45].常用于模擬兩相流動邊界條件包括周期邊界以及無滑移反彈邊界等.其中,周期邊界指流體粒子從一側(cè)邊界離開流場時,在下一時間步會從流場的另一側(cè)流入流場,周期邊界可以保證整個系統(tǒng)的質(zhì)量和動量守恒.傳統(tǒng)的周期邊界采用增加虛擬流體節(jié)點(x0,xN+1)的方法實現(xiàn),可以表示為[46]

對于靜止的固體邊界,常用的處理方法就是對邊界上的粒子進行反彈處理,反彈邊界可以嚴(yán)格保證邊界上速度無滑移,反彈邊界可以表示為[46]

式中,ybottom和ytop分別為位于底部及頂部的固體邊界.

2 模型驗證

2.1 熱力學(xué)一致性驗證

熱力學(xué)不一致性一直是偽勢模型中的主要問題之一,將直接影響氣液密度比及界面張力的計算精度.因此,需要首先評價模型的熱力學(xué)一致性.考慮到氣液彎液面引起的毛管壓力會對氣液密度有一定影響,采用假想的平直界面氣液模型來研究氣液共存密度,并與麥克斯韋理論計算結(jié)果進行對比來評價模型的熱力學(xué)一致性.

首先,建立31 × 201 的格子空間,在x和y方向均設(shè)置為周期邊界,只改變溫度,以模擬不同溫度下的氣液共存密度,采用Huang 等[47]密度初始化方法提高數(shù)值穩(wěn)定性,該方法可以表示為

式中,W為相界面寬度,一般為2～ 5 個格子,ρv和ρl分別為麥克斯韋理論的氣相和液相密度.采用上述密度初始化方法之后,在50 ≤y≤150 的位置為液相,其他位置為氣相,氣液界面為直線,可以忽略毛管壓力作用[47].對比β=1 和β=1.125 及β=1.25 三種不同參數(shù)條件下的氣液共存密度和麥克斯韋理論計算的氣液共存密度,結(jié)果如圖1 所示,β=1.125 時,LB 模型與麥克斯韋理論所得的氣液共存密度基本一致,據(jù)此采用β=1.125.

圖1 氣液共存密度曲線對比Fig.1 Comparison of gas-liquid coexistence density curves

2.2 界面張力驗證

在驗證熱力學(xué)一致性之后,通過模擬靜態(tài)圓形液滴來計算界面張力.首先,構(gòu)建3 種不同網(wǎng)格尺寸(101 × 101,151 × 151,201 × 201)的格子空間,模擬靜態(tài)圓形液滴密度分布并進行網(wǎng)格獨立性檢驗.采用Huang 等[48]提出的密度初始化方法

式中,R0為液滴初始化半徑,(x0,y0)為液滴初始圓心位置,模擬條件選取為:T=0.8Tc,ρv=0.342,ρl=6.60,四周采用周期邊界,其余參數(shù)與2.1 節(jié)中保持一致.為了便于比較,將不同網(wǎng)格尺寸中液滴的直徑與網(wǎng)格尺寸的比值保持一致,模擬時間步選為50 000步,足以使模擬結(jié)果達到穩(wěn)定.模擬結(jié)果如圖2 所示,不同網(wǎng)格尺寸下中心線處(x=Nx/2)對應(yīng)的氣相及液相密度分布幾乎一致,說明模擬結(jié)果具有網(wǎng)格獨立性.模擬的密度分布結(jié)果僅在兩相過渡區(qū)域存在微小偏差,當(dāng)網(wǎng)格密度超過151 × 151 時,密度分布的模擬結(jié)果幾乎一致,足以滿足模擬精度要求,為了節(jié)約計算資源量,選取網(wǎng)格尺寸為151 × 151 的格子空間進行界面張力模擬驗證.

圖2 中心線處密度分布Fig.2 Density distribution of center line

進一步模擬T=0.7Tc,T=0.75Tc以及T=0.8Tc三種不同溫度條件下圓形液滴的靜態(tài)行為,模擬達到平衡后根據(jù)等式(14)計算實際的氣液壓差,利用Young-Laplace 方程ΔP=Pin-Pout=σ/R確定界面張力[49].本文認為氣液的分界面處的密度為(ρv+ρl)/2,并以此為根據(jù)確定液滴直徑大小,通過改變液滴的不同直徑,獲得一系列的壓差和液滴半徑的對應(yīng)關(guān)系.如圖3 所示,壓差與半徑的倒數(shù)呈現(xiàn)良好的線性關(guān)系(R2均大于0.999),通過線性擬合得到不同溫度下氣液界面張力大小分別為0.463 lu,0.361 lu,0.264 lu (其中l(wèi)u 代表格子單位).

圖3 界面張力驗證Fig.3 Verification of interfacial tension

2.3 靜態(tài)接觸角驗證

本節(jié)在考慮流體-流體的作用力的基礎(chǔ)上,添加流體與固體壁面之間的作用力,模擬液滴與壁面間的不同接觸角.Chen 等[39]給出了一種液滴在壁面上的接觸角計算方法

式中,θ為液滴的潤濕角,ξ1為液滴與固體壁面接觸線的長度,ξ2為液滴的高度.為了模擬驗證不同的潤濕角,建立151 × 151 的格子區(qū)域,在x方向設(shè)置周期邊界,在y=0 及y=Ny設(shè)置為固體壁面,并采用反彈邊界[50],模擬溫度T=0.8Tc,密度初始化方法采用Li 等[51]提出的密度初始化方法,其具體格式如下

式中,R0為液滴初始化半徑,本文取15 lu,x0,y0為液滴初始化圓心位置,本文設(shè)置(x0,y0)=(75,135).模擬結(jié)果如圖4 所示,通過改變不同的壁面密度ρw可以模擬0°～ 180°的潤濕角.

圖4 不同靜態(tài)接觸角模擬Fig.4 Simulation of different static contact angles

2.4 角隅液相滯留特征

在方形截面孔隙及喉道中,氣相驅(qū)替液相之后,液相將會在孔隙及喉道的角隅處滯留,液相的滯留特征(曲率半徑)對非圓形孔隙中的卡斷機理有一定的影響[52].為此,本文采用提出的偽勢模型模擬液相在角隅的曲率半徑與液相飽和度Sl的關(guān)系,驗證模型的準(zhǔn)確性.建立151 × 151 的格子空間,模擬溫度T=0.8Tc,四周均為固體壁面并采用反彈邊界,潤濕角設(shè)置為0°,其余參數(shù)與2.3 節(jié)中一致.在文獻[51]基礎(chǔ)上提出一種新的初始化密度方法以表征液相在角隅中的不同曲率半徑,具體格式為

式中,(x0,y0)=(75,75) 為氣泡初始圓心坐標(biāo)位置,R1為初始氣泡半徑,將R1的取值范圍設(shè)為75～之間即可表征角隅液的不同曲率半徑(0～ 75).如圖5 所示,所提模型可以捕捉到壁面的液膜及液相在角隅處的曲率半徑.壁面上靜態(tài)液膜厚度的表征方法可以參考文獻[53-57],本文主要聚焦評價液相在角隅處的不同曲率直徑d及飽和度的Sl關(guān)系.

圖5 角隅液相滯留特征Fig.5 Corner liquid retention characteristics

通過將液相在角隅的曲率半徑與液相飽和度模擬結(jié)果與Li 等[58]提出的方形孔隙角隅液的曲率半徑與液相飽和度的關(guān)系相對比來驗證模型的準(zhǔn)確性,若不考慮壁面液膜厚度的影響,理論模型可以簡化為

式中,ε為無量綱曲率直徑,ε=d/(L1-2h),L1為方形孔隙的長度,h為液膜厚度,d為角隅液的曲率直徑.圖6 給出了模擬結(jié)果與理論結(jié)果的對比,可以看出,模擬結(jié)果與理論結(jié)果基本一致,說明所提模型可以準(zhǔn)確評價液相在角隅的賦存規(guī)律以及確定角隅液的曲率半徑.

圖6 角隅液曲率直徑與飽和度的關(guān)系Fig.6 Relationship between corner liquid curvature diameter and saturation

2.5 網(wǎng)格獨立性檢驗

通過模擬不同網(wǎng)格尺寸條件下單管內(nèi)氣相驅(qū)替液相過程(彎液面頂點位置隨時間變化),對模型進行網(wǎng)格獨立性檢驗.建立4 種不同的網(wǎng)格尺寸(241 ×41,301 × 51,361 × 61,421 × 71),將y方向的底部及頂部設(shè)置為固體邊界并采用無滑移反彈格式.在入口和出口分別設(shè)置為5～10 格子的氣相緩沖區(qū)域,其余位置設(shè)置為液相,在每個格子節(jié)點添加外力驅(qū)動流動,在流動方向(x方向)采用周期邊界,當(dāng)液相被驅(qū)替至出口的氣相緩沖區(qū)域中時,液相密度在下一個時間步被自動替換為氣相密度,實現(xiàn)緩沖區(qū)域內(nèi)液相“消失”,從而實現(xiàn)驅(qū)替過程[35].模擬條件選取為:T=0.8Tc,ρv=0.342,ρl=6.60,將模擬的驅(qū)替壓差Δp以及模擬時間步長與取點時間間隔的比值保持一致.模擬結(jié)果如圖7 所示,對于不同的網(wǎng)格尺寸,模擬結(jié)果幾乎一致.因此,在后續(xù)驅(qū)替的模擬中,選取的模擬條件網(wǎng)格尺寸為301 × 51.

圖7 網(wǎng)格獨立性驗證Fig.7 Grid independence verification

3 動態(tài)卡斷機理及影響因素分析

3.1 靜態(tài)卡斷準(zhǔn)則評價

Ransohoff 等[8]基于Roof 提出的卡斷判斷機制,(即喉道毛管力Pc-neck大于前緣毛管力Pc-front),提出了方形孔喉截面的靜態(tài)卡斷準(zhǔn)則.其核心假設(shè)是在方形截面孔喉系統(tǒng)中,液相(潤濕相)在角隅處以角流為主要流動方式,如圖8 所示,忽略壁面液膜的影響,縮頸方形截面孔喉系統(tǒng)的卡斷判斷準(zhǔn)則可以表示為

圖8 氣液兩相流動示意圖Fig.8 Schematic diagram of gas-liquid two-phase flow

式中,Rc為孔隙曲率半徑,Rt為喉道曲率半徑,Rzt為喉道橫向曲率半徑,Cm為無量綱曲率,與孔隙截面形狀有關(guān).基于最小表面能模型,Ransohoff 等[8]給出了方形孔隙中無量綱曲率Cm=1.89.當(dāng)孔隙收縮性較平緩時1/Rzt可以忽略,式(27)進一步簡化為Rc＞1.89Rt.

Ransohoff 等[8]提出的靜態(tài)準(zhǔn)則一般適用于驅(qū)替壓差較小的準(zhǔn)靜態(tài)情況.為了驗證該靜態(tài)準(zhǔn)則的適用性,本文采用改進的偽勢模型進行模擬驗證.為了方便將模擬結(jié)果與Ransohoff 靜態(tài)準(zhǔn)則進行對比,首先定義無量綱孔喉長度比L*及孔喉寬度比R*表征孔隙結(jié)構(gòu)的變化,L*和R*可以表示為

式中,Lp為整個孔喉系統(tǒng)的長度,L為喉道長度.基于此,可以得到卡斷靜態(tài)準(zhǔn)則為R*≤0.53.

建立301 × 51 的格子空間,在120 ＜x＜ 180 及10 ＜y＜ 40 的位置處設(shè)置寬度為30 lu 的喉道(L*=0.2,R*=0.6).x方向的邊界設(shè)置與2.5 節(jié)中一致,在y=0 及y=50 及喉道壁面設(shè)置為固體壁面,潤濕角設(shè)置為0°,并采用反彈邊界.模擬的兩相流體性質(zhì)為:氣相密度ρv=0.53 lu、液相密度ρl=6.08 lu、氣相黏度μv=0.088 lu、液相黏度μl=1.013 lu、兩相界面張力為σ=0.178 lu.選取的模擬密度比(ρr=11.5)與15 MPa,350 K 溫度條件下的水(980 kg/m3)與甲烷(90 kg/m3)的密度比(ρr=11)接近(NIST 化學(xué)數(shù)據(jù)庫).

為了更清楚的表示驅(qū)替過程,筆者將整個驅(qū)替過程分為三個階段(以時間為分界線).第一階段,氣相在左端的孔隙中的驅(qū)替階段;第二階段,氣相侵入喉道的驅(qū)替階段;第三階段,氣相突破喉道進入右端孔隙的驅(qū)替階段.在這三個階段中,選取氣相飽和度Sg以及滯留液相飽和度Srl作為量化表征液相變化的參數(shù).氣相飽和度Sg是指在驅(qū)替過程中氣相體積分?jǐn)?shù)與整個區(qū)域體積分?jǐn)?shù)之比;滯留液相飽和度Srl是指在驅(qū)替過程中,以兩相彎液面的頂點為分界線(垂直線),將整個區(qū)域分為兩部分,分界線以內(nèi)的液相體積分?jǐn)?shù)與整個區(qū)域體積分?jǐn)?shù)的比值.換句話說,Sg+Srl等效為石油工程領(lǐng)域中活塞式驅(qū)替(界面為平直界面)過程中氣驅(qū)波及的全部面積.捕捉了驅(qū)替過程中的Sg以及Sg+Srl隨時間的變化規(guī)律以及三種階段某一時刻下的兩相流體的密度分布(圖9標(biāo)注).第一階段(t＜ 9000),氣相在左端的孔隙中的驅(qū)替階段,驅(qū)替過程中形成相對穩(wěn)定的彎液面,Srl保持較低的數(shù)值;第二階段(9000≤t≤13 000),氣相侵入喉道的驅(qū)替階段,此時在孔隙角隅處及喉道壁面上存在滯留的液相,隨著模擬時間步的增大Srl逐漸增大;第三階段(t＞ 13 000),氣相突破喉道進入右端孔隙的驅(qū)替階段,此時由于壁面的強潤濕性,氣相無法與壁面接觸,導(dǎo)致右端孔隙壁面上滯留大量液相,Srl迅速增大,并且孔隙右側(cè)滯留的液相(潤濕相)隨時間逐漸發(fā)生回流到喉道中,最終卡斷相界面.此外,所提模型孔喉條件為(R*=0.6),并不滿足卡斷發(fā)生的靜態(tài)準(zhǔn)則條件(R*≤ 0.53),但是通過模擬可以觀察到卡斷現(xiàn)象.這是因為在驅(qū)替的第二和第三階段,氣相無法與喉道壁面以及右端孔隙壁面接觸,導(dǎo)致氣相突破喉道后的前緣曲率半徑始終小于孔隙半徑,使得Cm取值偏大,從而導(dǎo)致基于角流推導(dǎo)的靜態(tài)準(zhǔn)則預(yù)測結(jié)果偏保守.

圖9 不同時間步下的Sg 和Srl 變化Fig.9 Sg and Srl changes at different time steps

進一步計算了驅(qū)替過程中的第三階段內(nèi)的Pc-neck及Pc-front數(shù)值變化.如圖10 所示,隨著時間步的增大Pc-neck不斷增大,Pc-front不斷減小,這是因為Srl增大使得Pc-neck不斷增大,氣相突破喉道之后前緣曲率半徑不斷增大導(dǎo)致Pc-front不斷減小,最終Pc-neck＞Pc-front,使得滯留的液相發(fā)生回流從而卡斷相界面.

圖10 不同時間步下Pc-neck 及Pc-front 隨時間變化Fig.10 Pc-neck and Pc-front changes at different time steps

上述結(jié)果與之前靜態(tài)準(zhǔn)則判斷發(fā)生卡斷的條件一致,說明Roof[6]基于毛管壓力不平衡的卡斷判斷機制在突變孔喉系統(tǒng)中同樣適用.同時,上述2D 模擬結(jié)果與Wu 等[15]在3D 孔喉結(jié)構(gòu)中氣泡卡斷實驗結(jié)果基本一致,說明模擬結(jié)果對實際的3D 情況仍然有效.然而,大量研究表明,動態(tài)因素對卡斷現(xiàn)象有較大影響[11],但是目前的理論模型都是準(zhǔn)靜態(tài)幾何準(zhǔn)則,無法考慮動態(tài)因素的影響.為此,利用提出的偽勢模型,采用數(shù)值模擬方法評價毛管數(shù)Ca,L*及R*對氣液兩相驅(qū)替中相界面卡斷的影響規(guī)律.

3.2 毛管數(shù)的影響

毛管數(shù)Ca是表征兩相流動的重要參數(shù),其定義為

式中,ρn,νn分別代表非潤濕相的密度及運動黏度,在本文的模擬中與氣相的密度及運動黏度相同,uave為兩相驅(qū)替的平均速度[9].

通過改變驅(qū)替壓差大小進而改變毛管數(shù),研究毛管數(shù)對于卡斷的影響.建立301 × 51 的格子空間,在100 ＜x＜ 200 及10 ＜y＜ 40 的位置處設(shè)置長度為100 lu,寬度為30 lu 的喉道(L*=0.33,R*=0.6),驅(qū)替壓差Δp=0.069,其余條件均與3.1 中相同.每隔500 時間步,捕捉相界面前緣位置,通過相界面位置與時間關(guān)系曲線獲得兩相驅(qū)替平均速度(uave=Δs/Δt).需要指出的是,Δs和Δt表示發(fā)生卡斷現(xiàn)象之前的位移和時間.圖11(a)給出了Ca=4.8 × 10-3時兩相驅(qū)替過程中的Sg以及Sg+Srl隨時間的變化規(guī)律以及三個階段某一時刻下的兩相流體的密度分布,第一階段的模擬結(jié)果與3.1 節(jié)中一致;當(dāng)氣相進入喉道之后,孔隙角隅及喉道壁面始終存在滯留的液相,Srl逐漸增大.由于驅(qū)替壓差較小,氣相不能突破喉道,最終氣液相界面在喉道內(nèi)保持平衡不再延伸,Srl保持不變(Srl=0.11).造成上述現(xiàn)象的原因是氣液界面的毛管壓力以及固液界面的作用力形成較大的流動阻力,這種現(xiàn)象稱為毛管“釘扎”效應(yīng)[59].因此,發(fā)生卡斷的首要條件即為驅(qū)替壓差可以克服毛管“釘扎”效應(yīng),形成有效驅(qū)替.

增大驅(qū)替壓差至Δp=0.078,圖11(b)給出了Ca=5.6 × 10-3時兩相驅(qū)替過程的Sg以及Sg+Srl隨時間的變化規(guī)律以及三個階段某一時刻下的兩相流體的密度分布.驅(qū)替的第一階段與前文模擬結(jié)果一致;在氣相進入喉道之后,孔隙角隅及喉道壁面始終存在滯留的液相,Srl逐漸增大,且滯留在喉道壁面上的液相厚度呈非均勻分布.隨著模擬時間步的增大,滯留的液相逐漸回流并卡斷相界面.值得一提的是,發(fā)生卡斷時間屬于驅(qū)替的第二階段,說明在非漸變孔喉系統(tǒng)中卡斷可以在氣相突破喉道之前發(fā)生.

繼續(xù)增大驅(qū)替壓差至Δp=0.105,其余條件保持不變,圖11(c)給出了Ca=8.0 × 10-3時兩相驅(qū)替過程中的Sg以及Sg+Srl隨時間的變化規(guī)律以及三個階段某一時刻下的兩相流體的密度分布.此時毛管數(shù)相對較大,相比于之前的模擬結(jié)果,在第二階段,喉道壁面上滯留的液相飽和度Srl較小;隨著模擬的時間步的增大,氣相不斷向前流動,在突破喉道前并未發(fā)生卡斷,直至氣相到達出口端,隨著模擬時間步的繼續(xù)增大,滯留的液相發(fā)生回流并形成卡斷.并且隨著毛管數(shù)的增大,發(fā)生卡斷的位置向喉道出口端移動,這與Wei 等[60]采用3D 孔喉結(jié)構(gòu)研究不同毛管數(shù)條件下氣泡卡斷的實驗結(jié)果一致.

圖11 不同時間步下的Sg 和Srl 變化Fig.11 Sg and Srl changes at different time steps

進一步增大驅(qū)替壓差至Δp=0.114,其余條件保持不變,圖11(d)給出了Ca=8.6 × 10-3時兩相驅(qū)替過程中的Sg以及Sg+Srl隨時間的變化規(guī)律以及三個階段某一時刻下的兩相流體的密度分布.此時驅(qū)替壓差很大,相比于之前的模擬結(jié)果Srl較小;隨著時間步的增大,氣相從喉道中突破,并最終到達出口端.在氣相到達出口端以后,隨著模擬時間步的增大,并未發(fā)生卡斷現(xiàn)象.這說明對于固定的孔喉結(jié)構(gòu),存在一個臨界毛管數(shù),當(dāng)系統(tǒng)的毛管數(shù)大于臨界毛管數(shù)時,將會抑制滯留的液相在喉道內(nèi)發(fā)生回流,從而抑制卡斷現(xiàn)象的發(fā)生.

圖12 給出了50 000 時間步內(nèi),上述四種條件下相界面前緣頂點位置隨時間的變化.發(fā)生卡斷之后,把回退到喉道中的彎液面頂點作為實際的相界面位置,這樣當(dāng)相界面位置突然變化時,即可判斷為發(fā)生卡斷.可以看出,氣相為連續(xù)相時,驅(qū)替過程中可以發(fā)生連續(xù)卡斷,并且每次發(fā)生卡斷的位置幾乎不會變,該過程類似于液滴/氣泡的形成過程[61-62].同時,對于特定的孔喉系統(tǒng),只有在滿足一定的毛管數(shù)范圍內(nèi)才會發(fā)生卡斷,較大的毛管數(shù)會抑制驅(qū)替過程中滯留液相的回流從而抑制卡斷的發(fā)生,較小的毛管數(shù)無法克服毛管“釘扎”效應(yīng)形成有效驅(qū)替.模擬結(jié)果與文獻[9]在漸變孔喉中研究氣泡在壓力梯度作用下的卡斷行為所得出的結(jié)論一致.

圖12 相界面位置隨時間的變化Fig.12 Change of phase interface position with time

3.3 孔喉長度比的影響

孔喉結(jié)構(gòu)對卡斷同樣具有重要影響[9-10],為此本節(jié)利用提出的偽勢模型研究喉道長度對卡斷的影響.建立301 × 51 的格子空間,在10 ＜y＜ 40 的位置處設(shè)置相同寬度(R*=0.6),不同長度的喉道來研究不同L*對卡斷的影響.對于不同的喉道長度,驅(qū)替壓差比毛管數(shù)更能體現(xiàn)動態(tài)因素對卡斷的影響,因此采用驅(qū)替壓差代替毛管數(shù)表征動態(tài)參數(shù)的影響.

在模擬過程中,改變驅(qū)替壓差從0.036 至0.123,改變L*從0.08 至0.4 來確定卡斷發(fā)生的臨界條件;這里僅用Δp=0.084,L*=0.3 和Δp=0.084,L*=0.2 兩種模擬條件為例,展示不同孔喉長度比的驅(qū)替效果.圖13 給出了Δp=0.084,L*=0.3 及L*=0.2 條件下,不同時間步的驅(qū)替過程.隨著喉道長度的減小,流動阻力減小,兩相驅(qū)替平均速度增大,驅(qū)替過程中滯留在孔隙角隅及喉道壁面處的液相飽和度Srl較低.值得一提的是,設(shè)置的孔喉寬度比R*=0.6,不滿足靜態(tài)準(zhǔn)則預(yù)測發(fā)生卡斷的條件(R*≤0.53),但是模擬結(jié)果可以觀察到卡斷現(xiàn)象,說明即使不滿足靜態(tài)準(zhǔn)則預(yù)測發(fā)生卡斷的孔喉寬度比,足夠長的喉道長度也會促進卡斷的發(fā)生.因此,在相同驅(qū)替壓差的作用下,不同喉道長度的孔喉系統(tǒng)會呈現(xiàn)不同的流動狀態(tài).

圖13 兩相驅(qū)替過程對比Fig.13 Comparison of two-phase displacement processes

為了確定不同孔喉長度比的卡斷發(fā)生條件,開展了一系列模擬,得到不同孔喉長度比下發(fā)生卡斷的臨界驅(qū)替壓差,如圖14 所示.對于固定的喉道寬度的孔喉系統(tǒng),存在臨界孔喉長度比(L*=0.08),使得無論驅(qū)替壓差如何變化,氣液兩相驅(qū)替過程中不會出現(xiàn)卡斷現(xiàn)象,這與Deng 等[10]提出的短波長的漸變孔喉中卡斷會被抑制的結(jié)論一致.隨著喉道長度增大,不同的驅(qū)替壓差會使得氣液驅(qū)替過程中出現(xiàn)卡斷,連續(xù)流動以及無效驅(qū)替三種流動狀態(tài).當(dāng)驅(qū)替壓差小于臨界驅(qū)替壓差下界時,無法克服毛管“釘扎”效應(yīng),形成無效驅(qū)替;當(dāng)驅(qū)替壓差大于臨界驅(qū)替壓差上界時,氣相將保持連續(xù)流動的狀態(tài),抑制卡斷的發(fā)生.

圖14 不同L*條件下氣液流動狀態(tài)與驅(qū)替壓差的關(guān)系Fig.14 The relationship between gas-liquid flow state and displacement pressure difference with different L*

3.4 孔喉寬度比的影響

除了喉道長度,喉道的寬度對于卡斷現(xiàn)象也有一定影響,本節(jié)通過固定L*改變R*,研究不同喉道寬度對卡斷的影響.建立301 × 51 的格子空間,在125 ＜x＜ 175 設(shè)置長度為50 lu 的喉道(L*=0.167),其余參數(shù)設(shè)置與3.3 節(jié)相同.在模擬過程中,改變驅(qū)替壓差從0.033 至0.1,改變R*大小從0.52 至0.68 來確定卡斷臨界條件;這里僅用Δp=0.084,R*=0.52 和Δp=0.084,R*=0.6 兩組條件為例,展示不同孔喉寬度比對氣液兩相驅(qū)替過程中的流動狀態(tài)的影響.圖15 給出了Δp=0.084,R*=0.52 及R*=0.6 兩種條件下的驅(qū)替過程.隨著喉道寬度的減小,流動阻力顯著增大,兩相驅(qū)替平均速度減小,驅(qū)替過程中滯留在孔隙角隅及喉道壁面處的液相飽和度Srl較高.值得一提的是,R*=0.52 滿足靜態(tài)準(zhǔn)則預(yù)測發(fā)生卡斷的條件(R*≤ 0.53),模擬結(jié)果同樣觀察到了卡斷現(xiàn)象,驗證了模擬結(jié)果的準(zhǔn)確性.對于R*=0.6 的情況下,驅(qū)替過程中始終不會發(fā)生卡斷現(xiàn)象,這也說明了喉道寬度增大將會抑制卡斷的發(fā)生.

圖15 兩相驅(qū)替過程對比Fig.15 Comparison of two-phase displacement processes

為了確定不同孔喉寬度比對卡斷現(xiàn)象的影響,開展了一系列模擬,得到了不同孔喉寬度比下發(fā)生卡斷的臨界驅(qū)替壓差.如圖16 所示,與3.3 節(jié)類似,對于固定喉道長度的孔喉系統(tǒng),氣液兩相在不同驅(qū)替壓差的作用下形成氣相連續(xù)流動,卡斷及無效驅(qū)替3 類流動狀態(tài).隨著喉道寬度的增加,發(fā)生卡斷的驅(qū)替壓差范圍逐漸減小,當(dāng)孔喉寬度比增大至R*=0.68 時,無論施加多大的驅(qū)替壓差,喉道內(nèi)都不會發(fā)生卡斷現(xiàn)象.值得一提的是,對于R*=0.52 的孔喉結(jié)構(gòu)(滿足靜態(tài)準(zhǔn)則預(yù)測發(fā)生卡斷的條件),當(dāng)驅(qū)替壓差大于0.1 時,驅(qū)替過程中將不會發(fā)生卡斷現(xiàn)象,這說明即使達到靜態(tài)準(zhǔn)則所預(yù)測的卡斷條件,較大的驅(qū)壓差也可以抑制卡斷的發(fā)生.同時,基于等式(27)預(yù)測臨界孔喉寬度比為(R*=0.53),相比于通過模擬得出的臨界孔喉寬度比(R*=0.68)降低28.3%.

圖16 不同R*條件下氣液流動狀態(tài)與驅(qū)替壓差的關(guān)系Fig.16 The relationship between gas-liquid flow state and displacement pressure difference with different R*

4 結(jié)論

本文在原始單組分偽勢格子玻爾茲曼模型基礎(chǔ)上改進流體-流體作用力格式,建立了改進的偽勢模型,并在此基礎(chǔ)上模擬了孔喉系統(tǒng)中的氣液兩相在不同驅(qū)替條件下的流動狀態(tài),得出以下結(jié)論.

(1) 明確了非漸變孔喉結(jié)構(gòu)卡斷機理,即孔喉毛管壓力不平衡使得驅(qū)替過程中滯留的液相(潤濕相)隨時間逐漸發(fā)生回流是造成相界面卡斷的根本原因.但是,模擬過程中在孔隙角隅及喉道壁面觀察到大量滯留的液相,這導(dǎo)致基于角流的假設(shè)思想建立的靜態(tài)卡斷判斷準(zhǔn)則將會高估喉道右側(cè)氣泡的曲率半徑從而低估卡斷發(fā)生的條件.

(2) 揭示了動態(tài)因素驅(qū)替壓差(毛管數(shù))對氣液兩相驅(qū)替的影響規(guī)律.在非漸變孔喉系統(tǒng)中,只有 驅(qū)替壓差(毛管數(shù))大小在一定范圍時,喉道內(nèi)才會發(fā)生卡斷;超過該毛管數(shù)上界,即使?jié)M足靜態(tài)準(zhǔn)則條件,卡斷也會被抑制;低于該下界,無法完成驅(qū)替;同時毛管數(shù)的增大使得發(fā)生卡斷的位置向喉道出口端移動.

(3) 揭示了孔喉長度比對氣液兩相驅(qū)替的影響規(guī)律.對于固定寬度的孔喉系統(tǒng),即使不滿足靜態(tài)卡斷發(fā)生的孔喉寬度比(R*≤ 0.53),足夠長的喉道也可以促進卡斷的發(fā)生,隨著喉道長度增大,發(fā)生卡斷的驅(qū)替壓差范圍越大.并且存在一個臨界喉道長度,使得無論驅(qū)替壓差如何變化,喉道內(nèi)都不會發(fā)生卡斷,對于本文的模型,臨界孔喉長度比L*=0.08.

(4) 揭示了孔喉寬度比對氣液兩相驅(qū)替的影響規(guī)律.對于固定長度的孔喉系統(tǒng),喉道寬度越大,發(fā)生卡斷的驅(qū)替壓差范圍越小;并且存在一個臨界喉道寬度,使得無論驅(qū)替壓差如何變化,喉道內(nèi)都不會發(fā)生卡斷.對于本文的模型,臨界孔喉寬度比R*=0.68.若采用靜態(tài)準(zhǔn)則預(yù)測該模型發(fā)生卡斷的條件(R*=0.53),將會低估28.3%.