相間力對微孔通道內兩相驅替滲流影響規律的LBM研究

張一夫，梁 冰，張遵國

（1.遼寧工程技術大學 安全科學與工程學院，遼寧 葫蘆島 125105；2.礦山熱動力災害與防治教育部重點實驗室,遼寧 葫蘆島 125105；3.遼寧工程技術大學，遼寧 阜新 123000）

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

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

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

圖1 D2Q9離散速度模型

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

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

在多組分偽勢模型中，對多孔介質中兩相介質分別采用不同的分布函數進行描述。考慮外力項且碰撞項采用線性化（LBGK近似）的第k種組分粒子分布函數f的演化方程（格子-Boltzmann方程）為：

其表達式為：

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

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

混合流體的速度為：

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

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

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

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

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

式中：s為固相標識函數，其值取為1為固相，取0為流體相；Gkw為控制流體和固體表面相互作用力強度的參數，通過調節Gkw值，可以得到不同的表面潤濕性條件，其值取負代表為潤濕，取正代表為非潤濕。

在D2Q9模型中，Gkw值為：

當有重力存在時，可用下式表示：

式中：g為重力加速度。

當所有外力都給定時，混合流體的壓強p可由系統的狀態給定：

1.2 復雜曲線邊界處理

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

圖2 曲線邊界處節點類型示意圖

粒子碰撞僅發生在流體節點自身所在的位置，而遷移過程需沿著網格線方向以格子速度將相應的分布函數傳遞到相鄰節點。曲線邊界處粒子的遷移過程如圖3。

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

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

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

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

圖4 Bouzidi格式

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

利用二階線性插值格式獲得固體邊界附近流體節點的分布函數f值：

當q<1/2時：

當q≥1/2時：

1.3 數值方法驗證

1.3.1 液滴靜態接觸角的驗證

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

圖5 壁面潤濕性及穩態壁面接觸角模擬結果

從圖5的模擬結果來看：流固作用強度參數Gw數值的大小決定了液滴的穩定形態。隨著Gw由0.5變化到－0.5，壁面由疏水性轉變為親水性，液滴與壁面間的靜態接觸角θ由大逐漸變小，且呈近似線性規律變化（圖6）。從而驗證了偽勢模型對模擬兩相流固耦合作用的敏感性和準確性。

圖6 壁面接觸角θ與流固作用強度參數Gw之間的關系

1.3.2 毛細力液橋現象的驗證

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

圖7 不同流固作用強度參數Gw時，液橋或液滴形態的演變過程（粗線為最終穩定狀態）

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

2 相間力對多孔介質內的兩相驅替規律的影響

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

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

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

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

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

另外，從兩相界面的總體驅進速度對比來看（圖11），隨著煤體潤濕性的增強（Gw減小），兩相界面的總體驅進速度也得到了提高。

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

圖11 Gw取不同值時，兩相界面的總體驅進速度對比

3 結論

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

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

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