原狀花崗巖殘積土大孔隙細觀滲流場的格子Boltzmann模擬

蔡沛辰， 闕云， 楊鵬飛

(福州大學土木工程學院， 福建 福州 350108)

0 引言

花崗巖殘積土作為一種典型結構性土， 具有孔隙比大、 擾動性強等特點[1]. 除自身存在孔隙外， 受動植物孔隙、 土壤干濕凍融循環影響， 土壤內分布著許多大孔隙， 且大孔隙的存在導致水流繞過大部分基質土壤， 快速到達土壤深層， 即產生大孔隙流效應[2]， 并常伴有邊坡失穩、 水土流失等自然災害出現. 因此， 在細觀尺度上揭示花崗巖殘積土大孔隙流水分遷移機理日趨緊迫.

現階段大孔隙滲流的研究方法主要有試驗方法和數值模擬兩種. 在試驗研究方面， 呂捷等[3]采用降雨染色示蹤試驗， 對天然邊坡土體內大孔隙優勢流的空間響應進行了研究； Forrer等[4]借助亮高藍染色劑研究了土壤的優勢流現象. 這些試驗研究方法可以部分揭示土體滲流的一些大孔隙機理， 但這種研究無法細致準確地表征大孔隙結構特征、 滲流路徑以及孔隙結構與滲流特性之間的內在聯系. 近年來， 數值模擬在大孔隙滲流方面的研究主要采用宏觀連續模型、 計算流體力學(computational fluid mechanics， CFD)和格子Boltzmann等方法. 1) 基于宏觀連續模型方法. 馮杰[5]對存在大孔隙的土壤中水及溶質運移機制進行了研究， 并分析了大孔隙滲流的特性； Robert等[6]采用改進的MACRO模型模擬了土壤中水分優先遷移的情況， 結果表明考慮大孔隙優先流動可以更加豐富黏土土壤水態的一維描述. 2) 基于CFD方法. 汪鴻山等[7]以瀝青路面為研究對象， 研究了不同結構層孔隙率對路面滲透性能的影響. 相關研究已取得較大成果， 但宏觀連續模型方法將固相和水體均假設作為連續介質， 實際上固相是離散體， 連續介質處理方式勢必與實際孔隙結構存在一定差異， 對實踐不能起到很好的指導意義； 而最為成熟、 常用的CFD方法多應用于處理較為簡單的模型邊界條件， 且準確性與可信性難以保證. 3) 基于格子Boltzmann方法對原狀殘積土體大孔隙細觀滲流過程進行動態可視化研究. 格子Boltzmann方法最早在1989年由Mcnamara和Zanetti提出[8]， 隨后迅速在流體力學領域發展起來， 具有構思新穎、 算法簡單、 程序易于實現及并行計算等優點[9]. 學者Batany等[10]、 申林方等[11]、 崔冠哲等[12]采用此方法研究過孔隙滲流模擬， 但受限于孔隙結構理想化或考慮因素單一化， 且均未涉及到細觀尺度下原狀花崗巖殘積土大孔隙內的水分遷移模擬.

鑒于此， 本研究采用原狀花崗巖殘積土作為對象， 對其進行工業CT掃描， 建立真實反映土體孔隙結構的計算模型. 采用格子Boltzmann方法進行數值模擬， 直觀展現各大孔隙區域的滲流速度分布情況， 并分析不同時步和不同深度情況下的細觀滲流特性， 以期探究和認識不同條件下花崗巖殘積土大孔隙的細觀滲流規律.

1 CT掃描及細觀模型構建

試驗原狀土選自福州市西北一處坡地， 現場取樣如圖1所示， 選取植被茂密的位置， 取樣前先用鐵鍬清理表面腐殖層， 清理面積為100 cm×100 cm， 厚度為20 cm， 最終獲取尺寸為15 cm×15 cm×40 cm的土柱試樣， 測試土樣的基本物理參數， 結果如表1所示.

圖1 現場取樣Fig.1 Field sampling

表1 所選試樣的基本物理參數

在土壤孔隙結構特征質量方面， 工業CT掃描明顯高于醫學CT掃描[13]， 因此對獲取的原狀土采用英華公司C450KV高能量工業CT掃描儀進行XZ向掃描， CT掃描工作電壓為450 kV， 工作電流為63 mA， 掃描最低分辨率為0.15 mm， 之后得到一系列能真實體現土壤孔隙分布的CT掃描切片， 如圖2所示.

圖2 CT掃描切片 Fig.2 CT scan slice

通過Matlab軟件中的imread和im2bw等函數功能， 對CT掃描圖像進行二值化處理， 形成只包含黑白兩色的圖像， 再基于小波變換對二值化圖像進行降噪處理， 除去圖像中孤立的噪點， 最后選取中間位置的切片， 截取其中孔隙連通性較好的區域作為模擬對象， 如圖3所示， 大小為6 mm×3.2 mm.

圖3 土體細觀孔隙模型Fig.3 Soil meso-pore model

2 格子Boltzmann模型

2.1 格子Boltzmann方程

格子Boltzmann方程是Boltzmann-BGK方程的一種特殊離散形式[14]， 即可通過求解時間t、 位置ω處粒子分布函數F(ω,t)的離散 Boltzmann 方程來推導出Navier-Stokes方程， 從細觀尺度來模擬流體流動的規律[15]. 不含外力項時， 粒子分布函數F(ω,t)的演化可以通過離散的格子Boltzmann方程表示為

(1)

2.2 格子Boltzmann基本模型

格子Boltzmann模型一般由格子(離散速度模型)、 平衡態分布函數及分布函數的演化方程組成[16]. 1992年， Qian等[17]提出的DdQm系列模型是格子Boltzmann方法的基本模型， 常見的有一維模型D1Q3， 二維模型D2Q7、 D2Q9, 三維模型D3Q17、 D3Q19等， 其中d表示維度，m表示離散速度個數. 本研究選擇的是二維模型最為常用的D2Q9， 如圖4所示.

圖4 D2Q9模型Fig.4 D2Q9 model

D2Q9模型共有9個運動方向， 速度配置如下

(2)

式中:c=δx/δt，δx、δt分別為網格步長和時間步長， 通常取為1.

平衡態分布函數為

(3)

式中：cs為格子聲速；ρ為密度；eα為離散速度；u為宏觀速度；wα為權系數， D2Q9模型權系數配置如下

(4)

同時采用文獻[18]中Chapman-Enskog展開方法， 來推導BLE基本模型所對應的流體力學中的Navier-Stokes方程. 模型的宏觀密度、 速度及黏滯系數ν與無量綱弛豫時間τ之間關系如下

2.3 邊界條件處理

細觀滲流模擬中， 依據構建的計算模型， 進行邊界條件設定： 上下邊界為流體進出口， 采用非平衡態外推格式邊界處理方式； 左右邊界及內部孔壁為不透水層， 采用標準反彈格式邊界處理方式， 如圖5所示.

圖5 邊界條件處理 Figs.5 Boundary condition processing

(5)

(6)

圖5(b)為標準反彈格式邊界條件處理示意圖， 其基本思想是： 自流體節點(i-1, 2)入射到邊界及節點(i, 1)的分布函數F8， 不發生碰撞即原路返回， 由此獲得節點(i, 1)上的F6[20]， 其他節點相類似. 即分布函數可表示為

F2, 5, 6(i, 1)=F4, 7, 8(i, 1)

(7)

式中:F4, 7, 8(i, 1)可分別由F4(i, 2)、F7(i+1, 2)和F8(i-1, 2)遷移得到.

2.4 收斂性判別

本研究參考何雅玲等[14]在頂蓋驅動流中的收斂判據方式來判斷計算結果是否收斂， 給定收斂標準的小量ε=10-6. 若Error<ε， 則計算結果收斂， 否則返回計算. 判別式如下：

(8)

2.5 LBM算例驗證

為驗證本數值計算方法的正確性， 采用Poiseuille流對自編LBM程序的正確性進行驗證. 選取20×5的矩形網格區域作為驗證計算模型， 其邊界條件的處理辦法同上處理方式. 具體計算參數如表 2 所示. 其中，Nx和Ny為計算模型的長度和寬度；u為入口流體速度.在 LBM 數值計算中， 為便于計算和分析， 一

表2 驗證算例計算參數表(格子單位)

般采用無量綱的格子單位[15， 20]， 故所有參數均采用格子為單位.

采用LBM方法編程對選取的計算模型進行計算， 本文模擬結果與Poiseuille流理論值對照結果如圖6所示. 從圖中可以看出， 本文LBM模擬值與Poiseuille理論值擬合度較好， 計算知最大流速處誤差僅0.69%.

圖6 Poiseuille理論值與LBM模擬值對比Fig.6 Comparison of Poiseuille theoretical value and LBM simulation value

3 結果與討論

基于工業CT掃描技術及圖像處理方法構建的花崗巖殘積土大孔隙細觀模型， 借助格子Boltzmann方法模擬流體(水)在土體大孔隙中不同時步下的滲流情況. 為使得模擬結果更貼近于現實情況， 將土體中孔隙滲流方向沿土體的深度方向進行， 并設置模型以初始流速u=0.1入滲， 具體計算模型邊界條件設置如圖7所示. 此外， 選取430×240的網格作為計算區域， 土體模型的孔隙率φ為31.71%. 其他相關設置及計算參數均與上述驗證算例取值相同.

圖7 計算模型邊界條件Fig.7 Calculation model boundary conditions

圖8(a)～(f)分別列出了100、 300、 500、 2 000、 10 000和50 000時步的滲流場速度分布圖. 從圖中可以看出： 流體粒子在孔隙連通性好、 且孔徑較大的大孔隙區域快速流過， 呈閃電狀分布， 且流速較大， 而在部分孤立大孔隙中流體無法流通； 不同時步下滲流流體隨時步的增加， 分布范圍越加廣泛， 同時大孔隙優先流效應逐漸減小， 分析原因是所選取模型土體趨于飽和狀態. 如： 格子坐標(60, 200)處， 100、 300、 500、 2 000、 10 000和50 000時步下的滲流場格子速度分別為0.042 11、 0.053 22、 0.027 65、 0.011 37、 0.010 22、 0.010 21.

圖8 不同步長滲流場速度分布圖Fig.8 Asynchronous long seepage field velocity distribution diagram

將圖8與圖3土體細觀結構圖相對比可發現： 并非所有大孔隙都會出現優先流現象， 部分孤立或封閉的孔隙滲流速度近乎為0. 圖8(a)～(f)中分別對格子坐標范圍(Nx： 278～314，Ny： 196～240)處的大孔隙流的推進過程進行了局部放大(紅色圈內). 從中可以看出大孔隙滲流過程中存在較為明顯的指進效應， 孔道中間流速最大， 可達格子速度0.093 2， 越靠近孔壁滲流速度越小， 接近于0， 其結果較為符合Poiseuille圓管流理論.

圖8(c)～(f)分別對計算模型左側流體孔道中大孔隙區域進行了橢圓紅色圈注， 由圖中可以知道， 流體自進口端沿大孔隙存在的孔道到達出口端時約500時步， 相比其他孔道率先流出模型， 大孔隙優先流效應較為顯著.

為得到土體大孔隙滲流場流速穩定后不同深度格點位置處滲流速度沿水平方向的變化情況， 選取Z=240(模型入口邊界)、Z=160(模型1/3深度)、Z=80(模型2/3深度)和Z=0(模型出口邊界)格點所在深度進行研究， 結果繪制如圖9所示. 從圖中可以看出： 水平方向格點處有超過50%的區域滲流速度都為0， 對照圖3土體細觀結構模型可以發現大孔隙存在區域， 格點流速分布呈單峰狀或多峰狀；Z=0處格點流速分布雜亂， 無規律可尋， 特別是圖中水平方向格點0～100范圍內， 查看土體細觀結構模型發現此處存在多個復雜大孔隙相關聯. 總體來看Z=240(模型入口邊界)和Z=0(模型出口邊界)格點所在深度處滲流速度最大， 而Z=160(模型1/3深度)處格點速度普遍最小， 且幾乎所有格點位置處的流速均小于入口邊界流速u=0.1. 分析原因是： 花崗巖殘積土孔隙比大， 且原狀土體內部孔隙結構復雜無序， 這將對流體滲流過程勢必造成一定程度上的阻礙， 根據能量守恒定律可知， 穩定后流速必小于最初的初始速度.

圖9 滲流速度穩定后不同深度處格點速度分布示意圖(格子單位) Fig.9 Schematic diagram of grid velocity distribution at different depths after the seepage velocity stabilizes (grid unit)

4 結語

研究發現， 格子Boltzmann數值方法可有效對原狀花崗巖殘積土大孔隙細觀滲流進行計算， 采用Poiseuille流對程序的正確性進行驗證， 誤差僅為0.69%.

在孔隙連通性好、 孔徑較大的區域流體粒子以較高速度通過， 且呈閃電狀分布， 而在部分孤立孔隙中流體無法流通. 其中， 大孔隙中滲流存在明顯的“指進效應”， 孔道中間流速最大， 可達格子速度0.093 2， 越靠近孔壁滲流速度越小， 接近于0， 其結果較為符合Poiseuille圓管流理論.

最后， 研究不同深度格點位置處滲流速度沿水平向的變化情況， 發現大孔隙存在區域的格點流速分布呈單峰狀或多峰狀， 而多個大孔隙相交互區域內格點流速分布呈散點狀， 且無規律可尋.