申林方，李騰風(fēng)，王志良，李 澤，王鵬宇

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

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

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

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

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

1 壁面滑移邊界

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

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

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

2 格子Boltzmann 模型

2.1 格子Boltzmann 方程

2.2 邊界處理

2.3 單位轉(zhuǎn)換

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

2.4 收斂判斷依據(jù)

2.5 算法驗證

2.5.1 蒸汽中的懸浮液滴模擬

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

根據(jù)Laplace 定律，穩(wěn)定后液滴的內(nèi)外壓差與其半徑成反比，可表示為[24]：

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

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

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

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

圖5 液滴內(nèi)外壓差與半徑間的關(guān)系Fig. 5 Relationship between pressure difference and droplet radius

2.5.2 壁面接觸角模擬

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

表2 接觸角與壁面密度對應(yīng)關(guān)系Table 2 Relationship between contact angle and wall density

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

3 數(shù)值分析研究

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

圖7 計算模型Fig. 7 Computational model

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

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

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

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

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

巖體微裂隙隙寬變化改變了固液界面表面積與流體體積之比，從而進一步影響到壁面浸潤性與流體滲流間的關(guān)系。為此，針對重力作用下接觸角分別為57°、84°、114°和147°的四種裂隙，研究了裂隙寬度從0.08 mm 增加至0.32 mm 時，巖體微裂隙的滲流規(guī)律，其平均滲流流速與無滑移邊界相對誤差的變化規(guī)律，如圖8 所示。由圖可知，在壁面浸潤性相同的情況下，隨著裂隙隙寬的增加，其滲流流速與無滑移邊界的相對誤差逐漸減小。說明壁面浸潤性對微裂隙滲流特性的影響，隨隙寬的增加而逐漸減弱。此外，壁面的親/疏水性愈強，微裂隙隙寬對其滲流特性的影響愈發(fā)顯著。接觸角為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，研究了壓力梯度對微裂隙滲流特性的影響。裂隙平均流速與壓力梯度的對應(yīng)關(guān)系如圖9 所示，由圖可知，裂隙的平均滲流流速隨壓力梯度的增加而增大，且兩者呈線性關(guān)系，這一規(guī)律與無滑移邊界Poiseuille 流的理論公式相一致。同時，裂隙壁面疏水能力越強，其直線斜率越大，說明強疏水壁面對微裂隙滲流具有顯著的促進作用。對于親水或弱疏水壁面，其斜率變化并不顯著。

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

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

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

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

4 結(jié)論

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

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

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

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

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