氣液非混相驅替過程中的卡斷機理及模擬研究1)

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

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

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

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

引言

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

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

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

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

1 單組分多相偽勢模型

1.1 經典BGK 碰撞模型

在格子玻爾茲曼方法中,采用離散的密度分布函數表征流體的運動規律,標準BGK 碰撞過程可以表示為[36]

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

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

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

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

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

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

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

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

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

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

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

1.3 狀態方程

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

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

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

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

1.4 外力項處理格式

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

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

1.5 邊界條件

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

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

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

2 模型驗證

2.1 熱力學一致性驗證

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

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

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

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

2.2 界面張力驗證

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

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

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

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

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

2.3 靜態接觸角驗證

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

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

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

圖4 不同靜態接觸角模擬Fig.4 Simulation of different static contact angles

2.4 角隅液相滯留特征

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

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

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

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

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

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

2.5 網格獨立性檢驗

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

圖7 網格獨立性驗證Fig.7 Grid independence verification

3 動態卡斷機理及影響因素分析

3.1 靜態卡斷準則評價

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

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

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

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

式中,Lp為整個孔喉系統的長度,L為喉道長度.基于此,可以得到卡斷靜態準則為R*≤0.53.

建立301 × 51 的格子空間,在120 ＜x＜ 180 及10 ＜y＜ 40 的位置處設置寬度為30 lu 的喉道(L*=0.2,R*=0.6).x方向的邊界設置與2.5 節中一致,在y=0 及y=50 及喉道壁面設置為固體壁面,潤濕角設置為0°,并采用反彈邊界.模擬的兩相流體性質為:氣相密度ρ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 化學數據庫).

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

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

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

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

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

3.2 毛管數的影響

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

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

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

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

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

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

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

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

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

3.3 孔喉長度比的影響

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

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

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

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

圖14 不同L*條件下氣液流動狀態與驅替壓差的關系Fig.14 The relationship between gas-liquid flow state and displacement pressure difference with different L*

3.4 孔喉寬度比的影響

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

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

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

圖16 不同R*條件下氣液流動狀態與驅替壓差的關系Fig.16 The relationship between gas-liquid flow state and displacement pressure difference with different R*

4 結論

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

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

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

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

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