張一夫，梁 冰，張遵國(guó)

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

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

1 基于多組分偽勢(shì)模型的LBM方法

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

圖1 D2Q9離散速度模型

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

1.1 考慮相間力的多組分偽勢(shì)模型

在多組分偽勢(shì)模型中，對(duì)多孔介質(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個(gè)部分：流體組分之間的相互作用力、流體組分與固體壁面之間的作用力以及流體組分所受的重力。

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

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

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

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

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

在D2Q9模型中，Gkw值為：

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

式中：g為重力加速度。

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

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

當(dāng)研究孔隙尺度下的多相滲流問題時(shí)，必須要考慮固體骨架的具體形態(tài)對(duì)流動(dòng)的影響。由于LBM方法是通過碰撞和遷移過程來描述流動(dòng)粒子運(yùn)動(dòng)，因此對(duì)于流體粒子而言，在固體骨架附近的粒子的遷移和碰撞過程結(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 個(gè)格子方向分別以 e1、e2、e3、...、e8的速度進(jìn)行遷移。以發(fā)生在5方向xf點(diǎn)與xs點(diǎn)之間的遷移過程為例，由邊界節(jié)點(diǎn)xw上的分布函數(shù)（xw，t）運(yùn)動(dòng)到流體節(jié)點(diǎn)xf的表達(dá)式為：

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

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

圖4 Bouzidi格式

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

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

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

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

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

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

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

圖5 壁面潤(rùn)濕性及穩(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)證了偽勢(shì)模型對(duì)模擬兩相流固耦合作用的敏感性和準(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)的中心位置，上下邊界為固體壁面，采用半步長(zhǎng)反彈邊界條件，左右為周期性邊界條件。2 種組分流體的參數(shù)分別設(shè)置為，ρ1，in＝ρ2，out＝1.0，v1，in＝v2，out＝0.16，組分間作用力參數(shù) Gf為 2.7。不同流固作用強(qiáng)度參數(shù)Gw時(shí)，液橋或液滴形態(tài)的演變過程如圖7。

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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