相間力對微孔通道內(nèi)兩相驅(qū)替滲流影響規(guī)律的LBM研究

張一夫，梁 冰，張遵國

（1.遼寧工程技術(shù)大學(xué) 安全科學(xué)與工程學(xué)院，遼寧 葫蘆島 125105；2.礦山熱動力災(zāi)害與防治教育部重點(diǎn)實(shí)驗(yàn)室,遼寧 葫蘆島 125105；3.遼寧工程技術(shù)大學(xué)，遼寧 阜新 123000）

孔隙介質(zhì)尺度下的細(xì)觀多相滲流的研究具有十分廣泛的工程背景[1-2]，如煤層注水、石油天然氣開采、核廢料處理、有機(jī)物對土壤的污染范圍等都與細(xì)觀多相滲流的特點(diǎn)和機(jī)理有密切關(guān)系。對于孔隙尺度的研究，因其涉及到復(fù)雜的物理力學(xué)機(jī)制，對其內(nèi)部現(xiàn)象的描述較為困難。既要考慮到不同流體之間的相互作用，還要考慮到流體與孔隙骨架之間的流固耦合，特別是當(dāng)滲流通道達(dá)到了微孔尺度，許多宏觀尺度下的作用力，如重力、慣性力等不再是影響流動特點(diǎn)的主要因素，而像表面張力、毛細(xì)力及孔隙介質(zhì)表面的潤濕機(jī)理等作用更加凸顯，成為影響兩相滲流規(guī)律的主要因素。對兩相流動的數(shù)值研究，傳統(tǒng)的方法是依托Navier-stokes方程并結(jié)合對界面追蹤方程的求解來追蹤兩相界面運(yùn)動，常見的方法有流體體積函數(shù)法（VOF）[3]、水平集方法（Level set）[4]、Front tracking[5]和 phase-field[6]等。這些方法通常是用來描述少數(shù)大的界面運(yùn)動，但在孔隙尺度下，由于物理模型的復(fù)雜性，往往很難刻畫出流動的細(xì)節(jié)和大量細(xì)小、分散界面的準(zhǔn)確位置。格子Bohzmann方法(LBM)是近二十幾年發(fā)展起來的一種基于分子動力學(xué)的模擬方法[7]，具有清晰的物理背景。由于其微觀本質(zhì)和介觀特點(diǎn)，可以方便地描述出流動中不同相間的相互作用，特別是邊界條件易于實(shí)施的特點(diǎn)，在模擬復(fù)雜流動問題上具有天然的優(yōu)勢。近年來提出的基于動理學(xué)理論的多組分流動LBM模型主要有顏色模型[8]、偽勢模型[9]、自由能模型[10]等，其中偽勢模型因其能方便地描述相間作用力、自動追蹤相界面且計(jì)算程序簡單易行，在多相流領(lǐng)域得到充分發(fā)展、應(yīng)用最為廣泛。利用LBM方法在孔隙尺度上對微孔介質(zhì)內(nèi)的流動進(jìn)行數(shù)值模擬，擬通過改變相間力的大小和作用機(jī)制，找到反映微孔通道內(nèi)兩相滲流的運(yùn)動規(guī)律和特點(diǎn)。

1 基于多組分偽勢模型的LBM方法

取經(jīng)典的二維九速正方格子-玻爾茲曼模型，簡稱D2Q9模型，對格子速度eα進(jìn)行離散。采用單松弛時間的LBM多組分偽勢模型對多孔介質(zhì)內(nèi)的兩相驅(qū)替過程進(jìn)行研究。D2Q9離散速度模型如圖1。

圖1 D2Q9離散速度模型

式中：α為9個格子速度分量的方向；c為格子速度。

1.1 考慮相間力的多組分偽勢模型

在多組分偽勢模型中，對多孔介質(zhì)中兩相介質(zhì)分別采用不同的分布函數(shù)進(jìn)行描述?？紤]外力項(xiàng)且碰撞項(xiàng)采用線性化（LBGK近似）的第k種組分粒子分布函數(shù)f的演化方程（格子-Boltzmann方程）為：

其表達(dá)式為：

式中：u、ρ分別為混合流體的宏觀速度和密度；cs為格子聲速，取；wα為權(quán)重系數(shù)，在D2Q9模型中其值分別為w0=4/9，w1-4=1/9和w5-8=1/36。

每種組分流體的宏觀密度和速度由下式確定：

混合流體的速度為：

作用在第k種組分上的總作用力Fk包括3個部分：流體組分之間的相互作用力、流體組分與固體壁面之間的作用力以及流體組分所受的重力。

偽勢模型通過建立組分流體間的相互作用勢，得到來自于相鄰組分流體的作用力：

式中：ψk為第k種組分“有效密度”；Gkk為格林（Green）函數(shù)。

式中：g為第k種組分和第種組分之間相互作用的強(qiáng)度。

當(dāng)存在固-液界面時，為了考慮固體表面潤濕性對流體流動的影響時，固體表面對第k種組分流體的作用力可寫成：

式中：s為固相標(biāo)識函數(shù)，其值取為1為固相，取0為流體相；Gkw為控制流體和固體表面相互作用力強(qiáng)度的參數(shù)，通過調(diào)節(jié)Gkw值，可以得到不同的表面潤濕性條件，其值取負(fù)代表為潤濕，取正代表為非潤濕。

在D2Q9模型中，Gkw值為：

當(dāng)有重力存在時，可用下式表示：

式中：g為重力加速度。

當(dāng)所有外力都給定時，混合流體的壓強(qiáng)p可由系統(tǒng)的狀態(tài)給定：

1.2 復(fù)雜曲線邊界處理

當(dāng)研究孔隙尺度下的多相滲流問題時，必須要考慮固體骨架的具體形態(tài)對流動的影響。由于LBM方法是通過碰撞和遷移過程來描述流動粒子運(yùn)動，因此對于流體粒子而言，在固體骨架附近的粒子的遷移和碰撞過程結(jié)果的尤為重要。曲線邊界處節(jié)點(diǎn)類型示意圖如圖2。

圖2 曲線邊界處節(jié)點(diǎn)類型示意圖

粒子碰撞僅發(fā)生在流體節(jié)點(diǎn)自身所在的位置，而遷移過程需沿著網(wǎng)格線方向以格子速度將相應(yīng)的分布函數(shù)傳遞到相鄰節(jié)點(diǎn)。曲線邊界處粒子的遷移過程如圖3。

圖3 曲線邊界處粒子的遷移過程

圖3中xf位置的分布函數(shù)在碰撞后需沿著1到8 個格子方向分別以 e1、e2、e3、...、e8的速度進(jìn)行遷移。以發(fā)生在5方向xf點(diǎn)與xs點(diǎn)之間的遷移過程為例，由邊界節(jié)點(diǎn)xw上的分布函數(shù)（xw，t）運(yùn)動到流體節(jié)點(diǎn)xf的表達(dá)式為：

因此，對于直角正交網(wǎng)格下的LBM模型，在每一個碰撞-遷移循環(huán)中，遷移前均需利用插值獲得與物理邊界相交的網(wǎng)格線上的分布函數(shù)值（xw，t），插值結(jié)果即為此處的曲線邊界條件，據(jù)此插值分布函數(shù)值（xw，t）進(jìn)入式（12）的遷移過程。

為了獲得從邊界進(jìn)入流場內(nèi)部的分布函數(shù)值，采用Bouzidi等[11]提出的一種結(jié)合反彈格式和空間插值的邊界處理方法。Bouzidi格式如圖4。

圖4 Bouzidi格式

該方法通過其引入表征邊界位置與相鄰固體節(jié)點(diǎn)、流體節(jié)點(diǎn)間相對距離的線性插值因子q，從而可以更加精準(zhǔn)地描述固體邊界的位置。其表達(dá)式為：

利用二階線性插值格式獲得固體邊界附近流體節(jié)點(diǎn)的分布函數(shù)f值：

當(dāng)q<1/2時：

當(dāng)q≥1/2時：

1.3 數(shù)值方法驗(yàn)證

1.3.1 液滴靜態(tài)接觸角的驗(yàn)證

為了驗(yàn)證流固作用參數(shù)Gw對壁面潤濕性情況以及液滴靜態(tài)接觸角的影響，建立了液滴-固壁物理模型。將初始半徑為20的液滴置于100×100的二維格子系統(tǒng)中，液滴初始圓心置于（50,0），上下邊界為固體壁面，采用半步長反彈邊界條件，左右為周期性邊界條件。2 種組分流體的參數(shù)分別為 ρ1，in＝ρ2，out＝1.0，v1，in＝v2，out＝0.16，組分間作用力參數(shù) Gf為 2.4。壁面潤濕性及穩(wěn)態(tài)壁面接觸角模擬結(jié)果如圖5。

圖5 壁面潤濕性及穩(wěn)態(tài)壁面接觸角模擬結(jié)果

從圖5的模擬結(jié)果來看：流固作用強(qiáng)度參數(shù)Gw數(shù)值的大小決定了液滴的穩(wěn)定形態(tài)。隨著Gw由0.5變化到－0.5，壁面由疏水性轉(zhuǎn)變?yōu)橛H水性，液滴與壁面間的靜態(tài)接觸角θ由大逐漸變小，且呈近似線性規(guī)律變化（圖6）。從而驗(yàn)證了偽勢模型對模擬兩相流固耦合作用的敏感性和準(zhǔn)確性。

圖6 壁面接觸角θ與流固作用強(qiáng)度參數(shù)Gw之間的關(guān)系

1.3.2 毛細(xì)力液橋現(xiàn)象的驗(yàn)證

為進(jìn)一步驗(yàn)證流固作用力及兩相間表面張力模型施加的正確性，建立不同流固作用參數(shù)Gw下液橋的穩(wěn)定性模型。將初始寬度為10的液柱置于100×100的二維格子系統(tǒng)的中心位置，上下邊界為固體壁面，采用半步長反彈邊界條件，左右為周期性邊界條件。2 種組分流體的參數(shù)分別設(shè)置為，ρ1，in＝ρ2，out＝1.0，v1，in＝v2，out＝0.16，組分間作用力參數(shù) Gf為 2.7。不同流固作用強(qiáng)度參數(shù)Gw時，液橋或液滴形態(tài)的演變過程如圖7。

圖7 不同流固作用強(qiáng)度參數(shù)Gw時，液橋或液滴形態(tài)的演變過程（粗線為最終穩(wěn)定狀態(tài)）

從圖7的模擬結(jié)果來看，流固作用參數(shù)強(qiáng)度Gw對穩(wěn)態(tài)液柱形態(tài)有顯著影響，相間作用力及流固作用力是決定能否形成穩(wěn)定液橋的主要影響因素。當(dāng)相間作用力（相間力作用參數(shù)Gf）不變時，流固作用力過大（圖 7（a））或者過小（圖 7（d））都很難形成穩(wěn)定的液橋，只有在合適的流固作用力下（圖7（b）和圖7（c））才能形成穩(wěn)定的液橋。

2 相間力對多孔介質(zhì)內(nèi)的兩相驅(qū)替規(guī)律的影響

為探討相間力對孔隙介質(zhì)兩相滲流的影響，建立500×100個格子單位的規(guī)則孔隙結(jié)構(gòu)，上下邊界為固體壁面，采用半步長反彈邊界條件，左右Zuo-He邊界條件，入口流速為0.1。2種組分流體的參數(shù)分別設(shè)置為 ρ1，in＝ρ2，out＝1.0，v1，in＝v2，out＝0.16，組分間作用力參數(shù)Gf為2.1。分別取不同的Gw值，對多孔介質(zhì)內(nèi)兩相驅(qū)替過程的進(jìn)行了數(shù)值模擬，得到了兩相界面演變規(guī)律（圖8）。

從圖8中可以看出，隨著兩相流體的驅(qū)替過程向右推進(jìn)到孔隙介質(zhì)區(qū)域時，相界面開始發(fā)生扭曲、破碎和變形。隨著時間的推移，原本連續(xù)的相界面變成了大量細(xì)小、分散的界面。當(dāng)孔隙介質(zhì)表面設(shè)置為中性潤濕性界面（Gw取0）時，驅(qū)替流體前緣界面向前推進(jìn)時，后續(xù)流體由于受孔隙介質(zhì)骨架的阻礙，推進(jìn)較慢，導(dǎo)致兩相流體的相界面被拉得很長、黏性指進(jìn)現(xiàn)象明顯、驅(qū)替效率較低。

當(dāng)孔隙介質(zhì)表面的親水性逐漸增強(qiáng)時（Gw=-0.5和 Gw=-1.0）（圖 9 和圖 10），隨著驅(qū)替流體相界面的前緣向前推進(jìn)，后續(xù)驅(qū)替流體依靠相間力和流固作用力，加快其在孔隙骨架表面上的潤濕速度，驅(qū)替流體的整體推進(jìn)速度較快、黏性指進(jìn)現(xiàn)象減弱、驅(qū)替效果提高。

圖9 Gw取-0.5時兩相界面的時空演變規(guī)律

圖10 Gw取-1.0時兩相界面的時空演變規(guī)律

另外，從兩相界面的總體驅(qū)進(jìn)速度對比來看（圖11），隨著煤體潤濕性的增強(qiáng)（Gw減?。?，兩相界面的總體驅(qū)進(jìn)速度也得到了提高。

由此可見，相間力對微孔通道內(nèi)兩相滲流驅(qū)替過程影響較為明顯，隨著孔隙介質(zhì)表面潤濕性的增強(qiáng)，兩相驅(qū)替的總體效果增強(qiáng)，驅(qū)替效果較好。因此，為了提高煤層的注水效果，可以從改善煤體表面潤濕性角度出發(fā)，通過注入液體改性來提高煤體表面的潤濕性，特別是低滲透煤層，由于煤體滲透性差，孔隙盲孔、角隅較多，容易形成封閉氣，通過改性能使在在孔隙通道內(nèi)兩相驅(qū)替過程中表面張力、毛細(xì)力的作用更加凸顯，明顯改善注水效果和驅(qū)替速度。

圖11 Gw取不同值時，兩相界面的總體驅(qū)進(jìn)速度對比

3 結(jié)論

1）基于偽勢模型的多相流格子波爾滋蔓方法可以準(zhǔn)確地反映出相間力對微孔通道內(nèi)兩相滲流軀替過程力學(xué)特性的影響規(guī)律。

2）從模擬結(jié)果的分析來看，兩相驅(qū)替的總體效果與相間力有較大關(guān)系，多孔骨架的親水性越強(qiáng)，兩相滲流的總體驅(qū)替效果就越好，從而有利于驅(qū)替速度的提高，反之亦然。

3）可以通過改進(jìn)注入液體在煤體孔隙表面的潤濕性，減少低滲透煤層注水過程中封閉氣體的形成，從而有效提高注水效果。