原狀花崗巖殘積土大孔隙細(xì)觀滲流場(chǎng)的格子Boltzmann模擬

蔡沛辰， 闕云， 楊鵬飛

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

0 引言

花崗巖殘積土作為一種典型結(jié)構(gòu)性土， 具有孔隙比大、 擾動(dòng)性強(qiáng)等特點(diǎn)[1]. 除自身存在孔隙外， 受動(dòng)植物孔隙、 土壤干濕凍融循環(huán)影響， 土壤內(nèi)分布著許多大孔隙， 且大孔隙的存在導(dǎo)致水流繞過(guò)大部分基質(zhì)土壤， 快速到達(dá)土壤深層， 即產(chǎn)生大孔隙流效應(yīng)[2]， 并常伴有邊坡失穩(wěn)、 水土流失等自然災(zāi)害出現(xiàn). 因此， 在細(xì)觀尺度上揭示花崗巖殘積土大孔隙流水分遷移機(jī)理日趨緊迫.

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

鑒于此， 本研究采用原狀花崗巖殘積土作為對(duì)象， 對(duì)其進(jìn)行工業(yè)CT掃描， 建立真實(shí)反映土體孔隙結(jié)構(gòu)的計(jì)算模型. 采用格子Boltzmann方法進(jìn)行數(shù)值模擬， 直觀展現(xiàn)各大孔隙區(qū)域的滲流速度分布情況， 并分析不同時(shí)步和不同深度情況下的細(xì)觀滲流特性， 以期探究和認(rèn)識(shí)不同條件下花崗巖殘積土大孔隙的細(xì)觀滲流規(guī)律.

1 CT掃描及細(xì)觀模型構(gòu)建

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

圖1 現(xiàn)場(chǎng)取樣Fig.1 Field sampling

表1 所選試樣的基本物理參數(shù)

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

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

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

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

2 格子Boltzmann模型

2.1 格子Boltzmann方程

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

(1)

2.2 格子Boltzmann基本模型

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

圖4 D2Q9模型Fig.4 D2Q9 model

D2Q9模型共有9個(gè)運(yùn)動(dòng)方向， 速度配置如下

(2)

式中:c=δx/δt，δx、δt分別為網(wǎng)格步長(zhǎng)和時(shí)間步長(zhǎng)， 通常取為1.

平衡態(tài)分布函數(shù)為

(3)

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

(4)

同時(shí)采用文獻(xiàn)[18]中Chapman-Enskog展開(kāi)方法， 來(lái)推導(dǎo)BLE基本模型所對(duì)應(yīng)的流體力學(xué)中的Navier-Stokes方程. 模型的宏觀密度、 速度及黏滯系數(shù)ν與無(wú)量綱弛豫時(shí)間τ之間關(guān)系如下

2.3 邊界條件處理

細(xì)觀滲流模擬中， 依據(jù)構(gòu)建的計(jì)算模型， 進(jìn)行邊界條件設(shè)定： 上下邊界為流體進(jìn)出口， 采用非平衡態(tài)外推格式邊界處理方式； 左右邊界及內(nèi)部孔壁為不透水層， 采用標(biāo)準(zhǔn)反彈格式邊界處理方式， 如圖5所示.

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

(5)

(6)

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

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]在頂蓋驅(qū)動(dòng)流中的收斂判據(jù)方式來(lái)判斷計(jì)算結(jié)果是否收斂， 給定收斂標(biāo)準(zhǔn)的小量ε=10-6. 若Error<ε， 則計(jì)算結(jié)果收斂， 否則返回計(jì)算. 判別式如下：

(8)

2.5 LBM算例驗(yàn)證

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

表2 驗(yàn)證算例計(jì)算參數(shù)表(格子單位)

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

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

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

3 結(jié)果與討論

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

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

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

圖8 不同步長(zhǎng)滲流場(chǎng)速度分布圖Fig.8 Asynchronous long seepage field velocity distribution diagram

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

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

為得到土體大孔隙滲流場(chǎng)流速穩(wěn)定后不同深度格點(diǎn)位置處滲流速度沿水平方向的變化情況， 選取Z=240(模型入口邊界)、Z=160(模型1/3深度)、Z=80(模型2/3深度)和Z=0(模型出口邊界)格點(diǎn)所在深度進(jìn)行研究， 結(jié)果繪制如圖9所示. 從圖中可以看出： 水平方向格點(diǎn)處有超過(guò)50%的區(qū)域滲流速度都為0， 對(duì)照?qǐng)D3土體細(xì)觀結(jié)構(gòu)模型可以發(fā)現(xiàn)大孔隙存在區(qū)域， 格點(diǎn)流速分布呈單峰狀或多峰狀；Z=0處格點(diǎn)流速分布雜亂， 無(wú)規(guī)律可尋， 特別是圖中水平方向格點(diǎn)0～100范圍內(nèi)， 查看土體細(xì)觀結(jié)構(gòu)模型發(fā)現(xiàn)此處存在多個(gè)復(fù)雜大孔隙相關(guān)聯(lián). 總體來(lái)看Z=240(模型入口邊界)和Z=0(模型出口邊界)格點(diǎn)所在深度處滲流速度最大， 而Z=160(模型1/3深度)處格點(diǎn)速度普遍最小， 且?guī)缀跛懈顸c(diǎn)位置處的流速均小于入口邊界流速u(mài)=0.1. 分析原因是： 花崗巖殘積土孔隙比大， 且原狀土體內(nèi)部孔隙結(jié)構(gòu)復(fù)雜無(wú)序， 這將對(duì)流體滲流過(guò)程勢(shì)必造成一定程度上的阻礙， 根據(jù)能量守恒定律可知， 穩(wěn)定后流速必小于最初的初始速度.

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

4 結(jié)語(yǔ)

研究發(fā)現(xiàn)， 格子Boltzmann數(shù)值方法可有效對(duì)原狀花崗巖殘積土大孔隙細(xì)觀滲流進(jìn)行計(jì)算， 采用Poiseuille流對(duì)程序的正確性進(jìn)行驗(yàn)證， 誤差僅為0.69%.

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

最后， 研究不同深度格點(diǎn)位置處滲流速度沿水平向的變化情況， 發(fā)現(xiàn)大孔隙存在區(qū)域的格點(diǎn)流速分布呈單峰狀或多峰狀， 而多個(gè)大孔隙相交互區(qū)域內(nèi)格點(diǎn)流速分布呈散點(diǎn)狀， 且無(wú)規(guī)律可尋.