考慮壁面浸潤性的光滑巖體微裂隙滲流特性數值模擬研究

申林方，李騰風，王志良，李 澤，王鵬宇

(昆明理工大學建筑工程學院，昆明，云南 650500)

巖體在長期的自然風化、地質構造變動及人類活動等作用下，形成了大量平均隙寬僅為0.2 mm的微裂隙[1]。這些微裂隙構成了地下水運移的主要通道，影響著巖體的整體滲透特性。而巖體微裂隙的滲流問題廣泛存在于實際工程中[2-3]，在油氣開采時人們希望提高巖體滲透率以增加開采率；而在水利、地下空間開發、核廢棄物處置等工程建設中，又希望降低其滲透性增加結構安全。因此研究巖體微裂隙的滲流機制，具有非常重要的理論和工程應用價值。

由于巖體微裂隙的隙寬微小，固液界面表面積與液體體積之比較大，使得固液壁面間的浸潤性顯著影響流體的運動特性。與宏觀裂隙相比，微裂隙具有滲透性弱、持水性強、導水性差、隱蔽性強等特點。同時，巖石通常由多種礦物組成，而不同礦物成分的浸潤性存在著顯著差異，如：文象花崗巖中石英的水接觸角約為144.0°、螢石的水接觸角約為74.6°等[4]，從而進一步加劇了巖體微裂隙滲流的復雜性。因此，有必要考慮壁面浸潤性影響對巖體微裂隙的滲流機制展開研究。

在微尺度下，流體受到裂隙壁面浸潤性的影響，導致宏觀條件下的無滑移邊界不再適用[5]。針對于微尺度滲流特性研究，國內外學者在實驗研究、理論研究及數值模擬等方面，已進行了大量的研究工作，并取得了較好的進展。在實驗研究方面：Holt 等[6]通過實驗研究發現當水通過微尺度碳納米管時，其流速比Poiseuille 流理論解高3 個數量級。Song 等[7]采用實驗方法研究了超疏水表面毫米級氣-水界面的滑移速度?；羲乇蟮萚8]基于水在超疏水和超親水微通道內的滲流實驗，指出超疏水微通道內的流動阻力顯著降低，最大可達25%?？傮w而言，采用實驗方法研究雖然可以直觀得到微尺度下的滲流規律，但其花費高、耗時長，且實驗結果往往受控于實驗方法、儀器精度等[9]。在理論研究方面：Navier 最早注意到固體壁面附近的滑移現象，并認為壁面處流體的滑移速度與其剪切速率成正比[10]。Spikes 等[11]通過引入滑移模型，建立了新的描述牛頓流體運動的數學方程。Yang 等[12]研究了流體在疏水壁面矩形微通道內的流動，并提出了一種改進方法用以確定疏水平行板微通道的滑移系數。Gennes[13]根據剪應力與摩擦阻力間的關系，推導出疏水壁面滑移長度的計算公式。對于微通道滲流的理論研究，通常其基本假定過于理想化，所得結論只能用于定性分析，與實際情況尚存在一定的偏差。隨著計算機技術的快速發展，數值計算方法在研究微通道內的流體流動方面得到了廣泛應用。Martell 等[14]基于直接數值模擬方法(DNS)，研究了超疏水表面在微通道中的減阻性能，并討論了不同粗糙度影響下的滑移速度、壁面剪應力及雷諾應力等。劉趙淼和逄燕[15]模擬了液體在圓形和梯形截面微通道內的流動，并分析了不同壓力差、微通道尺寸和表面粗糙度等因素對流動摩擦系數的影響。楊大勇和劉瑩[16]基于有限單元法，研究了微通道表面粗糙度對其電滲流特性的影響。然而，常規的宏觀數值計算方法，無法真實的反映微通道固液壁面間的作用機制，考慮壁面浸潤性影響從微觀的角度揭示其滲流特征。

基于分子動理論發展起來的Boltzmann 方法(LBM)，是一種宏觀離散、微觀連續的介觀數值計算方法。它通過考慮不同相不同組分間分子作用力的影響，描述流體與流體以及流體與固體之間的微觀作用機制，能夠彌補宏觀數值計算方法的缺陷，從微觀角度出發解釋宏觀物理現象，因此在微通道多相流體界面動力學的研究中得到了廣泛應用。付宇航等[17]基于偽勢格子Boltzmann 方法模擬了傾斜壁面浸潤性梯度驅動液滴的運動過程。黃橋高等[18]采用格子Boltzmann 方法研究了壁面浸潤性對疏水表面滑移流動及減阻特性的影響規律。Genty 等[19]基于格子Boltzmann 方法，研究了深部非飽和泥巖微裂隙內部的有效擴散系數。然而考慮壁面浸潤性影響的多相流問題分布于不同的研究領域，現有成果相對分散且各自研究的側重點也有所不同，尚不足以全面揭示巖體微裂隙的滲流機制。鑒于此，本文基于Shan-Chen偽勢模型的格子Boltzmann 方法，考慮壁面浸潤性的影響建立了光滑巖體微裂隙滲流的數值模型，結合兩個經典算例驗證了計算模型的有效性，并討論了裂隙壁面浸潤性、裂隙隙寬、壓力梯度及流體黏滯性等因素對其滲流特性的影響。

1 壁面滑移邊界

當流體在裂隙中流動時，經典的宏觀流體力學認為壁面附近流體速度為零，即無滑移邊界。隨著裂隙隙寬的減小，固液壁面接觸面積與流體體積之比隨之增大，導致壁面與流體間的相互作用不可忽略，從而促使與壁面接觸的部分流體產生切向速度，發生相對運動，稱之為滑移邊界[20]。

圖1 分別為無滑移、黏滯層、正滑移邊界的示意圖。其中，us為壁面處的流體流速，在無滑移和黏滯層邊界中us=0。當壁面強疏水時，壁面與附近流體之間存在排斥力，此時產生正滑移。反之，當壁面強親水時，壁面處流體在吸引力作用下，易形成“類固體層”，即黏滯層。

圖1 滑移邊界類型Fig. 1 Slip boundary type

2 格子Boltzmann 模型

2.1 格子Boltzmann 方程

2.2 邊界處理

2.3 單位轉換

本文將所有計算參數進行無量綱化處理，并確保無量綱化前后的流動準則數一致，將無量綱參數作為聯系的橋梁，實現物理單位與格子單位之間的單位轉換。對于本文的滲流問題，引入如下無量綱參數：

2.4 收斂判斷依據

2.5 算法驗證

2.5.1 蒸汽中的懸浮液滴模擬

基于Shan-Chen 偽勢模型的格子Boltzmann 方法，模擬了懸浮在蒸汽中液滴的演化過程。計算區域網格數為128× 128 格子，計算參數采用格子單位，假定初始狀態下蒸汽(密度ρg為0.032)中存在一個橢圓形液滴(密度ρl為0.279)，氣液相粒子間的相互作用強度G =1/3，液滴的長軸與短軸之比為1.25，如圖4 所示，計算模型四周均采用周期性邊界。在氣液兩相分子間作用力的影響下，液滴逐漸由橢圓形演化為圓形，并趨于穩定。將穩定后蒸汽與液滴密度平均值的位置作為氣液相界面，并用于計算液滴半徑。

根據Laplace 定律，穩定后液滴的內外壓差與其半徑成反比，可表示為[24]：

式中： pl為液滴內壓力； pg為蒸汽中的壓力；γgl為氣液表面張力；r 為穩定后的液滴半徑。

圖4 懸浮在蒸汽中的液滴Fig. 4 A suspended liquid droplet in vapor

本文模擬了長軸長度從12 增至40 格子數的情況下，橢圓形液滴的演化過程，并得到了穩定后液滴內外壓差與其半徑間的關系，如圖5 所示。由圖可知，圓形液滴內外壓差與其半徑的倒數呈線性關系，線性擬合曲線關系式為：

擬合曲線的相關系數為 R2=0.9986，這與Laplace 定律相符，也充分證明了計算程序的正確性。

圖5 液滴內外壓差與半徑間的關系Fig. 5 Relationship between pressure difference and droplet radius

2.5.2 壁面接觸角模擬

為了研究固體壁面的浸潤性，假定液態水在固體壁面形成液滴，基于偽勢Shan-Chen 模型的格子Boltzmann 方法，模擬了不同浸潤性壁面的接觸角。計算模型采用200× 150 的網格系統，計算參數均采用格子單位。初始狀態下，在固體壁面的中央設置一個半徑r=30 的半圓形液滴，其液相密度為ρl=0.279，液滴周圍的氣相密度為ρg=0.032，兩相粒子間的相互作用強度G =1/3，計算模型四周均采用周期性邊界。通過改變虛擬壁面密度ρs從0.06 增至0.18，使不同相間的分子作用力發生變化，從而改變了固體壁面的浸潤性。液滴在分子間作用力的影響下，其形態逐漸發展演化，穩定后得到不同浸潤性壁面的接觸角θ。經計算得到，不同浸潤性固體壁面接觸角與虛擬壁面密度的對應關系如表2 所示，其典型接觸角的幾何形態如圖6 所示。在本文后續的分析討論中，以此為基礎，通過設置不同的虛擬壁面密度，選取相應的接觸角。

表2 接觸角與壁面密度對應關系Table 2 Relationship between contact angle and wall density

圖6 不同固體壁面接觸角的幾何形態Fig. 6 Geometric shapes of different solid wall contact angles

3 數值分析研究

為了研究光滑巖體微裂隙的滲流特性，基于Shan-Chen 偽勢模型的格子Boltzmann 方法，考慮壁面浸潤性影響建立了相應的數值模型。計算模型如圖7 所示，微裂隙的L×W=1.6 mm×0.2 mm，劃分成800×100 的網格，在初始狀態下裂隙中充滿液體，且液體的運動黏度ν=1.0×10-6m2/s，裂隙滲流的驅動力為重力，系統中計算參數物理單位與格子單位間的轉換，如表1 所示?；诮⒌奈⒘严稘B流數值模型，研究了裂隙壁面浸潤性、裂隙隙寬、壓力梯度、流體黏滯性等因素對微裂隙滲流特性的影響，并將計算結果與無滑移邊界的Poiseuille 流進行了對比分析。

圖7 計算模型Fig. 7 Computational model

在壓力驅動下，無滑移邊界Poiseuille 流平均流速的計算公式為[28]：

3.1 壁面浸潤性對微裂隙滲流特性的影響

為了研究壁面浸潤性對巖體微裂隙滲流特性的影響，針對裂隙隙寬W=0.2 mm 不同接觸角的裂隙方案，進行了滲流場的數值計算，并將計算結果與無滑移邊界的解析解進行了對比，其相對誤差 δr如表3 所示。由表3 可知，對于疏水壁面，由于其對附近流體的排斥作用而產生加速效果，使得平均滲流流速大于無滑移邊界，且隨著裂隙壁面接觸角的增加(疏水程度增強)，其相對誤差逐漸增大。而親水壁面由于壁面附近流體受到吸引而產生黏滯效應，使得其平均滲流流速小于無滑移邊界，且隨著接觸角的減小(親水程度增強)，相對誤差逐漸增大。但從總體變化趨勢看，疏水壁面對微裂隙滲流的影響要比親水壁面顯著。

表3 不同接觸角微裂隙的平均滲流流速Table 3 Average velocity in a micro-fracture with different contact angles

3.2 隙寬對微裂隙滲流特性的影響

巖體微裂隙隙寬變化改變了固液界面表面積與流體體積之比，從而進一步影響到壁面浸潤性與流體滲流間的關系。為此，針對重力作用下接觸角分別為57°、84°、114°和147°的四種裂隙，研究了裂隙寬度從0.08 mm 增加至0.32 mm 時，巖體微裂隙的滲流規律，其平均滲流流速與無滑移邊界相對誤差的變化規律，如圖8 所示。由圖可知，在壁面浸潤性相同的情況下，隨著裂隙隙寬的增加，其滲流流速與無滑移邊界的相對誤差逐漸減小。說明壁面浸潤性對微裂隙滲流特性的影響，隨隙寬的增加而逐漸減弱。此外，壁面的親/疏水性愈強，微裂隙隙寬對其滲流特性的影響愈發顯著。接觸角為147°、114°、84°及57°的壁面，微裂隙隙寬從0.32 mm 降至0.08 mm 的過程中，其平均流速與無滑移邊界的相對誤差分別增加了7.42%、2.07%、1.18%及3.19%。

圖8 平均流速相對誤差隨微裂隙隙寬的變化趨勢Fig. 8 Relative error variation of average velocity with micro-fracture width

3.3 壓力梯度對微裂隙滲流特性的影響

壓力梯度以及流體與固體壁面間的作用力，共同影響著巖體微裂隙的滲流特性。為此，保持裂隙隙寬W=0.2 mm 不變，針對接觸角為147°、114°、84°及57°的壁面，將壓力梯度G*從0.6ρg 逐漸增加至1.8ρg，研究了壓力梯度對微裂隙滲流特性的影響。裂隙平均流速與壓力梯度的對應關系如圖9 所示，由圖可知，裂隙的平均滲流流速隨壓力梯度的增加而增大，且兩者呈線性關系，這一規律與無滑移邊界Poiseuille 流的理論公式相一致。同時，裂隙壁面疏水能力越強，其直線斜率越大，說明強疏水壁面對微裂隙滲流具有顯著的促進作用。對于親水或弱疏水壁面，其斜率變化并不顯著。

圖9 微裂隙平均流速與壓力梯度的關系Fig. 9 Relationship between average velocity of microfracture and pressure gradient

3.4 流體黏滯性對微裂隙滲流特性的影響

巖體微裂隙中流體的黏滯性，隨著溫度、溶質成分及濃度等因素的不同而有所差異。為此，在裂隙寬度W=0.2 mm 及壓力梯度G*=ρg不變的情況下，分析了流體運動黏度從0.6×10-6m2/s 增加至1.6×10-6m2/s 的過程中，壁面接觸角為57°、84°、114°和147°時，裂隙滲流平均流速的變化規律。微裂隙平均滲流流速與流體運動黏度間的關系，如圖10 所示。由圖可知，隨著流體運動黏度的增加，流動阻力增大，微裂隙滲流的平均流速逐漸減小，兩者呈反比例關系，這與Poiseuille 流的解析式(26)相符。當流體的運動黏度相同時，裂隙滲流的平均流速隨接觸角的增加而增大，且對于強疏水壁面其影響尤為顯著。

圖10 微裂隙平均流速與流體運動黏度的關系Fig. 10 Relationship between average velocity of microfracture and kinematic viscosity of fluid

4 結論

基于Shan -Chen 偽勢模型的格子Boltzmann 方法，考慮壁面浸潤性的影響建立了光滑巖體微裂隙滲流的數值模型，結合蒸汽中懸浮液滴和壁面接觸角的模擬，證明了該模型的準確性和有效性。最后，考慮巖體壁面浸潤性、裂隙隙寬、壓力梯度及流體黏滯性等因素影響，研究了巖體微裂隙的滲流特性，得出以下結論。

(1)疏水壁面對附近流體的排斥作用，使得微裂隙的平均滲流流速要大于無滑移邊界，且隨著疏水程度的增強，兩者間的相對誤差逐漸增大；而親水壁面對流體的黏滯作用使其產生減速效果?？傮w而言，疏水壁面對微裂隙滲流的影響比親水壁面更加顯著。

(2)壁面浸潤性對光滑巖體微裂隙滲流特性的影響，隨著隙寬的減小而逐漸增強，且壁面的親/疏水性越強，微裂隙隙寬對其滲流特性的影響也越顯著。

(3)巖體微裂隙的平均滲流流速隨壓力梯度的增加而增大，且兩者呈線性關系。同時，由于強疏水壁面對微裂隙滲流具有顯著的促進作用，導致其斜率較大。

(4)隨著流體運動黏度的增加，流動阻力增大，導致微裂隙滲流的平均流速逐漸減小，流體的運動黏度與滲流平均流速呈反比例關系。