基于QSGS 法3D 重構(gòu)土體滲流場的LBM 數(shù)值模擬

闕云 ，邱婷，蔡沛辰，馬宏巖，謝秀棟，薛斌

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

在基坑、邊坡、堤壩、隧道等工程項(xiàng)目中，滲流作用往往會(huì)引起土體或結(jié)構(gòu)的變形失穩(wěn)，嚴(yán)重威脅工程安全［1-6］，明晰土體內(nèi)部滲流機(jī)理是預(yù)防和治理結(jié)構(gòu)滲流破壞的關(guān)鍵一環(huán).目前，國內(nèi)外研究者多通過宏觀滲流試驗(yàn)、數(shù)值模擬方法研究土體的滲流特性，但傳統(tǒng)滲流試驗(yàn)方法不可避免地會(huì)對原有結(jié)構(gòu)產(chǎn)生擾動(dòng)，試驗(yàn)數(shù)據(jù)將存在一定測量誤差，且反映的是宏觀滲流表現(xiàn)，無法從根本上揭示多孔介質(zhì)孔隙內(nèi)部的滲流過程，以及孔隙參數(shù)對滲流特性的影響機(jī)理.因此，從細(xì)觀尺度出發(fā)探究土體滲流特性尤為必要.

當(dāng)前細(xì)觀滲流研究主要集中在兩方面：重構(gòu)多孔介質(zhì)模型和細(xì)觀滲流數(shù)值模擬.1）多孔介質(zhì)模型重構(gòu)方法.其主要包括球體沉降法［7］、硬球Monte-Carlo法［8］、分?jǐn)?shù)布朗運(yùn)動(dòng)法［9］、隨機(jī)生長方法［10-11］、模擬退火法（Simulate Anneal Arithmetic，SAA）［12］和CT掃描重構(gòu)技術(shù).采用前三種方法構(gòu)建多孔介質(zhì)模型相對復(fù)雜［13］，故常采用后三種方法.SAA法重構(gòu)模型與真實(shí)土體孔隙結(jié)構(gòu)較為接近［14］，隨機(jī)法具有高效便捷且節(jié)約計(jì)算機(jī)內(nèi)存資源的優(yōu)勢，但上述方法存在生成孔隙結(jié)構(gòu)的尺寸和位置難以控制的缺點(diǎn).為解決此問題，Wang 等［11］提出了四參數(shù)隨機(jī)生長法（Quartet Structure Generation Set，QSGS）.研究土體孔隙參數(shù)對滲透特性影響情況時(shí)采用QSGS 法比CT 掃描法更為方便、快捷，CT 掃描法更適合研究真實(shí)土體孔隙結(jié)構(gòu)內(nèi)部的滲流特性.2）細(xì)觀滲流數(shù)值模擬方法.其主要包括基于連續(xù)介質(zhì)模型的CFD方法［15］、基于分子動(dòng)力學(xué)模型的格子Boltzmann 方法（Lattice Boltzmann Method，LBM）［16］.LBM 方法具有可準(zhǔn)確求解流體流動(dòng)Navier-Stokes 方程以及處理復(fù)雜幾何邊界的優(yōu)勢，廣泛應(yīng)用于多相流和多組分流動(dòng)中，如多相流密度對比模擬、非混相和部分混相驅(qū)替過程模擬、兩相滲流模擬等［17-18］.

由于QSGS 能與LBM 方法靈活對接實(shí)現(xiàn)土體細(xì)觀滲流場仿真模擬，諸多研究者已采用QSGS-LBM法研究了重構(gòu)土體微細(xì)觀滲流問題，如周瀟等［19］基于QSGS 法重構(gòu)土體微觀結(jié)構(gòu)，通過LBM 方法建立飽和土體滲流模型，探討了微觀土體結(jié)構(gòu)中水的滲流變化規(guī)律.申林方等［20］對傳統(tǒng)四參數(shù)隨機(jī)生長法進(jìn)行改進(jìn)，并采用LBM 方法研究了恒定流速下重構(gòu)土體的細(xì)觀滲流場.Jun 等［21］基于QSGS-LBM 方法分析了蒸壓加氣混凝土砌塊的平面滲流場特性.蔡沛辰等［22］采用QSGS 法構(gòu)建土體模型，采用LBM 方法研究重構(gòu)土在不同條件下的細(xì)觀滲流機(jī)理.李滔等［23］通過QSGS-LBM 方法分析了多孔介質(zhì)滲透率與孔隙尺度各向異性、孔隙分布非均質(zhì)性的關(guān)系.

雖然QSGS-LBM 方法在細(xì)觀滲流仿真中取得了較為豐碩的成果，但仍存在以下不足：一方面，傳統(tǒng)QSGS 法重構(gòu)土體模型存在土顆粒團(tuán)大小、形狀較為均一的缺陷，并不符合真實(shí)土體細(xì)觀結(jié)構(gòu)中土顆粒團(tuán)分布情況.另一方面，LBM滲流場仿真中模型尺寸大小與計(jì)算機(jī)內(nèi)存容量及運(yùn)行速度之間存在矛盾，而最優(yōu)模型尺寸的選取是解決問題的關(guān)鍵.在細(xì)觀滲流研究領(lǐng)域，關(guān)于孔隙參數(shù)、孔隙結(jié)構(gòu)各向異性等因素對3D滲流場特性影響機(jī)理的研究相對欠缺.

鑒于此，本文采用傳統(tǒng)與改進(jìn)QSGS法等效重構(gòu)各向同性/各向異性3D土體孔隙模型，綜合孔隙連通性和運(yùn)算時(shí)間因素選取最優(yōu)仿真模型，并結(jié)合LBM進(jìn)行滲流場數(shù)值模擬，直觀展現(xiàn)各孔隙區(qū)域內(nèi)滲流速度及流線分布情況，分析模型尺寸、孔隙結(jié)構(gòu)、孔隙率、土體顆粒大小及模型各向異性等因素下重構(gòu)土體滲流場的細(xì)觀滲流特性.

1 土體細(xì)觀結(jié)構(gòu)表征及構(gòu)建

1.1 細(xì)觀結(jié)構(gòu)表征

作為一種多孔介質(zhì)材料，土體通常可劃分為土體顆粒（固相）和孔隙區(qū)域，其空間分布函數(shù)如下［20］：

式中：G（x）為隨機(jī)變量，反映孔隙分布情況，且G（x）的平均值為.

1.2 四參數(shù)隨機(jī)生長法

四參數(shù)隨機(jī)生長法［11，13］可通過控制分布概率pc、生長概率pd、概率密度pirs（在i方向上第r相在第s相上的生長概率）和孔隙率n四個(gè)參數(shù)重構(gòu)土體細(xì)觀孔隙結(jié)構(gòu).令固相（土體顆粒）作為生長相，孔隙作為非生長相，則3D多孔介質(zhì)模型的重構(gòu)流程如下：

1）在3D 空間中，令生長相以一定的分布概率pc隨機(jī)布設(shè)，pc必須小于設(shè)定的孔隙率n.

2）在3D 空間中，令分布的生長相單元以一定的生長概率pd朝19個(gè)運(yùn)動(dòng)方向生長，如圖1所示.

3）重復(fù)上述步驟，當(dāng)生長相達(dá)到設(shè)定的孔隙率n時(shí)終止生長，3D多孔介質(zhì)模型重構(gòu)完畢.

控制分布概率pc=0.01（大顆粒），孔隙率n=0.36、0.50 和0.64，各向生長概率為一定值0.01（假設(shè)土體為各向同性），經(jīng)過生長，得到不同孔隙率下60×60×60格點(diǎn)大小的3D模型，如圖2所示.

圖2 不同孔隙率的 QSGS 模型（以pc=0.01為例）Fig.2 QSGS model with different porosity（taking pc=0.01 as an example）

1.3 改進(jìn)四參數(shù)隨機(jī)生長法

為更好地研究孔隙分布各向異性對3D 多孔介質(zhì)滲透能力的影響，采用兩個(gè)依次遞進(jìn)的生長過程來構(gòu)建3D 多孔介質(zhì)模型，當(dāng)各方向的生長速度不同時(shí)，即可獲得各向異性3D模型，詳細(xì)步驟如下：

步驟1：大團(tuán)顆粒生成較大的固相，固相生長核概率為pt，生長概率為pd，其中d代表大團(tuán)顆粒的19個(gè)生長方向，孔隙率為n1（pt

步驟2：小顆粒團(tuán)生成較小的固相，固相生長核概率為ps，生長概率為pk，其中k代表小顆粒團(tuán)的19個(gè)生長方向，且ps遠(yuǎn)小于pt.

兩個(gè)步驟相互結(jié)合，直至達(dá)到設(shè)定的孔隙率，最終生成具有各向異性的3D重構(gòu)土體細(xì)觀模型.

各向異性的參數(shù)可通過原狀土體CT 切片中孔隙的典型結(jié)構(gòu)和幾何形態(tài)特征確定，若同時(shí)調(diào)整各個(gè)方向的pd可生成與CT 掃描模型良好吻合的各向異性多孔介質(zhì)模型.本文模型通過設(shè)定孔隙率n=0.597（真實(shí)土體孔隙率平均值），生長概率p1～7=0.01、p8～11=0.2、p12～18=0.01、p19=0.005，所生成的3D重構(gòu)土體細(xì)觀模型如圖3所示，圖中灰色區(qū)域?yàn)榭紫?、黑色區(qū)域?yàn)榛|(zhì).

圖3 考慮各向異性的 3D 重構(gòu)模型Fig.3 Three-dimension anisotropic soil reconstruction model

2 格子 Boltzmann 模型

2.1 格子 Boltzmann 方程

格子Boltzmann 方程是離散形式的Boltzmann-BGK 方程［24］，其可從細(xì)觀層面仿真模擬多孔介質(zhì)中的流體滲流過程.不含外力項(xiàng)時(shí)，粒子分布函數(shù)的演化方程可表示為：

式中：Fα（ω，t）為格點(diǎn)ω位置處在t時(shí)刻沿α方向的粒子分布函數(shù)；eα為離散速度；δt為離散時(shí)間；τ為弛豫時(shí)間；Feqα(ω，t)為局部平衡態(tài)粒子分布函數(shù).

2.2 格子 Boltzmann 基本模型

格子Boltzmann 模型通常由三部分組成，即離散速度模型（格子）、平衡態(tài)分布函數(shù)和分布函數(shù)的演化方程［25］.本文選擇的模型是3D 層面最為常用的D3Q19模型，如圖4所示.

圖4 D3Q19 模型Fig.4 D3Q19 model

D3Q19 模型一共有如圖4 所示的19 個(gè)運(yùn)動(dòng)方向，離散速度eα滿足式（3）.

式中：c=δx/δt，其中δx為網(wǎng)格步長，δt為時(shí)間步長，均取1.

D3Q19模型的權(quán)系數(shù)wα如下：

DdQm模型的局部平衡態(tài)粒子分布函數(shù)可表示為：

式中：cs為格子聲速，取值為；ρ為密度；wα為D3Q19模型的權(quán)系數(shù)；u為宏觀速度.

同時(shí)基于LBE 基本模型，采用Chapman-Enskog展開方法［26］推導(dǎo)所對應(yīng)的Navier-Stokes 方程.密度ρ、宏觀速度u、壓力差p及運(yùn)動(dòng)黏度系數(shù)μ與弛豫時(shí)間τ（無量綱）之間關(guān)系如下：

2.3 真實(shí)單位與格子單位轉(zhuǎn)換

針對實(shí)際問題的數(shù)值仿真，常用的編程思路有兩種：程序代碼中物理量直接采用實(shí)際物理單位；程序代碼中物理參數(shù)都做無量綱化處理［24］.本文采用的是后者，即無量綱化的格子單位，參考Succi 著作［27］中的方法對LBM 單位換算部分進(jìn)行闡述.模型中基本參數(shù)（格子單位）：長度L、密度ρ、時(shí)間t、壓力差p和運(yùn)動(dòng)黏度系數(shù)μ.模型對應(yīng)的參數(shù)（物理單位）：長度為L′、密度為ρ′、時(shí)間為t′、壓力差為p′、運(yùn)動(dòng)黏度系數(shù)為μ′.為實(shí)現(xiàn)上述兩者的轉(zhuǎn)換，需引入部分參考量，即參考長度Lr、參考密度ρr和參考速度ur［24］.定義為：

式中：c′和c分別為物理單位和格子單位下的聲速.

對于一個(gè)特定的問題，模擬中的長度L、密度ρ、聲速c和運(yùn)動(dòng)黏度系數(shù)μ是已知的；實(shí)際物理量也可以通過相關(guān)公式獲得.故而參考密度ρr、參考速度ur可以確定，但L′與Lr仍無法確定.鑒于此，補(bǔ)充如下關(guān)系式：

此外，時(shí)間t、壓力p與t′、p′之間的轉(zhuǎn)換，可基于下式進(jìn)行求解：

至此，格子與實(shí)際物理單位間的換算完成［24］.通常，可采用下式將格子單位下的δx=δy=δz=1，δt=1 及c2=1/3，轉(zhuǎn)換成物理單位.

2.4 邊界條件處理

本文采用標(biāo)準(zhǔn)反彈格式［28］模擬滲流土體顆粒與流體間無滑移的滲流行為，其基本思想為：流體節(jié)點(diǎn)從（α-1，t）入射到（α，t）的粒子分布函數(shù)Fα1，遷移碰撞后沿原路徑返回，據(jù)此獲得節(jié)點(diǎn)（α，t）處的粒子分布函數(shù)Fα0，其余節(jié)點(diǎn)類似，粒子分布函數(shù)Fα0可表示為：

式中：α1 為指向壁面的方向；α0 為α1 的相反方向；xb為固體壁面處的格點(diǎn)；xf為流體格點(diǎn)，xf=xb-eαδt.

2.5 收斂性判別

參考文獻(xiàn)［28］中收斂判別方法判斷結(jié)果是否收斂.給定一個(gè)小量ε=10-6，若Error＜ε，則計(jì)算結(jié)果收斂，否則返回計(jì)算.3D判別式如下：

2.6 LBM算例驗(yàn)證

通過Poiseuille 流驗(yàn)證自編LBM 程序的正確性.選取60×60×100 的網(wǎng)格區(qū)域作為驗(yàn)證計(jì)算的3D 模型，具體計(jì)算參數(shù)如表1所示.

表1 LBM驗(yàn)證算例的計(jì)算參數(shù)表（格子單位）Tab.1 Calculation parameter of LBM verification example（grid unit ）

本文所有計(jì)算參數(shù)單位均為格子單位［20，30］，邊界條件處理同2.4節(jié).

對比自編LBM 程序計(jì)算結(jié)果與Poiseuille 流理論值，如圖5 所示.從圖5 中可知，LBM 計(jì)算值與Poiseuille流理論值基本吻合，流速最大誤差僅為1.47%.

圖5 Poiseuille 理論值與 LBM 計(jì)算值對比（格子單位）Fig.5 Comparison of Poiseuille theoretical value and LBM calculated value（grid unit ）

3 結(jié)果與討論

基于QSGS重構(gòu)技術(shù)和LBM對不同類型的3D土體模型進(jìn)行滲流仿真模擬，設(shè)定模型初始?jí)翰睿簆in=1.000 6 流入，pout=1.000 0 流出，邊界條件詳細(xì)設(shè)定如圖6所示，其他參數(shù)設(shè)置參見2.6節(jié)算例.

圖6 計(jì)算模型邊界條件Fig.6 Boundary conditions of the calculation of model

3.1 最優(yōu)仿真模型尺寸確定方法

本文選用pc=0.01，n=0.64，尺寸為30、40、50、60、80、100、120、150 和200 的立方體模型，通過分析模型尺寸大小對滲流穩(wěn)定時(shí)間和孔隙連通情況的影響來確定數(shù)值仿真最優(yōu)模型尺寸的大小.滲流穩(wěn)定狀態(tài)判斷標(biāo)準(zhǔn)［19］：當(dāng)格點(diǎn)離散速度不再隨時(shí)步增加而發(fā)生變化時(shí)，可認(rèn)為滲流已經(jīng)達(dá)到穩(wěn)定狀態(tài)，圖7為滲流穩(wěn)定后不同典型尺寸模型的速度場圖像.

圖7 不同典型尺寸模型滲流穩(wěn)定后速度場切片F(xiàn)ig.7 Velocity field slices after seepage stability of different typical size models

采用Matlab 軟件編程，統(tǒng)計(jì)3D 模型的孔隙參數(shù)，統(tǒng)計(jì)結(jié)果如表2所示.從表2中可看出：隨模型尺寸擴(kuò)大，孔隙連通程度明顯增大，模型尺寸擴(kuò)大至100×100×100 格點(diǎn)大小后，繼續(xù)擴(kuò)大模型尺寸時(shí)發(fā)現(xiàn)連通孔隙率nc增加并不明顯，再結(jié)合圖7（b）分析可知，隨模型尺寸增大滲流穩(wěn)定所需時(shí)間呈直線形增大.鑒于此，綜合孔隙連通性和運(yùn)算時(shí)間因素，本文最終選取100×100×100格點(diǎn)大小的多孔介質(zhì)模型作為最優(yōu)仿真單元模型，用于后續(xù)研究.該研究方法可為其他類似3D巖土體最優(yōu)仿真單元選取提供一定的借鑒.

表2 不同尺寸模型孔隙參數(shù)表（格子單位）Tab.2 Pore parameters of different size models（grid units ）

圖8 為典型模型尺寸對滲流速度和穩(wěn)定時(shí)間的影響情況.從圖8 中可發(fā)現(xiàn)：模型尺寸越大，格點(diǎn)速度u先隨之增大后減小，而穩(wěn)定所需時(shí)間t始終呈增長趨勢，統(tǒng)計(jì)知尺寸30、40、50、60、80、100、120、150和200的模型所需時(shí)間分別為212 min、657 min、853 min、1 091 min、1 957 min、2 661 min、4 582 min、6 283 min和103 043 min（運(yùn)算所采用計(jì)算機(jī)配置為 Intel Core i5-10400@3.60GHz 8G RAM，64 bit operating system，Windows10）.分析原因是：由表2 結(jié)論可知，模型尺寸越大，孔隙間的連通率隨之增大，連通孔隙增多，即速度也將增大，但當(dāng)模型尺寸增大到一定界限時(shí)，連通孔隙率nc將會(huì)保持在一個(gè)穩(wěn)定值左右；同時(shí)隨模型尺寸擴(kuò)大，模型的孔隙結(jié)構(gòu)也更加復(fù)雜，流體在孔隙結(jié)構(gòu)中來回碰撞、遷移的次數(shù)隨之增加，這將使得滲流達(dá)到穩(wěn)定所需的時(shí)間增加，故速度在一定程度上會(huì)降低.

圖8 3D模型尺寸對滲流時(shí)間和速度的影響Fig.8 Influence of 3D model size on seepage time and velocity

3.2 孔隙結(jié)構(gòu)對滲流場的影響

以pc=0.01，n=0.64，100×100×100格點(diǎn)大小的模型為例進(jìn)行滲流模擬，探討孔隙結(jié)構(gòu)特征對滲流過程的影響.

圖9（a）～（d）分別列出了滲流穩(wěn)定后不同切面的滲流速度及流線分布示意圖.從圖9（a）～圖9（c）中可明顯地發(fā)現(xiàn)：流體粒子自入口到出口，整個(gè)過程中格點(diǎn)速度呈逐漸減小趨勢.同時(shí)結(jié)合三維模型孔隙結(jié)構(gòu)特征分析可知，土體的孔隙結(jié)構(gòu)特征對流速影響較大.如：格子坐標(biāo)（1，12，46）處，格點(diǎn)速度為2.12×10-3，遠(yuǎn)高于其他孔隙區(qū)域速度.再觀察圖9（d）流線圖可知，QSGS 重構(gòu)的三維孔隙模型待滲流穩(wěn)定后，速度分布較為集中，且存在x向主滲流通道.

圖9 3D滲流模擬計(jì)算結(jié)果（pc=0.01，n=0.64）Fig.9 Results of three-dimensional seepage simulation（pc=0.01，n=0.64）

3.3 孔隙率對滲流場的影響

以pc=0.01，n=0.36、0.50、0.64、0.78，100×100×100 格點(diǎn)大小的模型為例，分析孔隙率變化情況下，滲流穩(wěn)定后滲流速度及滲透率變化情況.

圖10（a）～圖10（d）分別列出了不同孔隙率的三維模型滲流穩(wěn)定后典型切面的滲流速度分布示意圖.從圖10 中可看出：隨著模型孔隙率逐漸增大，滲流流體在孔隙中的分布范圍更加廣泛且均勻，但不論模型孔隙率大小，速度最大值都出現(xiàn)在孔隙較大處.

圖10 不同孔隙率的3D滲流模擬結(jié)果（pc=0.01）Fig.10 Results of 3D seepage simulation with different porosity（pc=0.01）

圖11 為不同孔隙率的模型沿x向截面的平均滲流速度變化曲線.從圖11 中不難發(fā)現(xiàn)，總體上孔隙率越大，平均滲流速度就越大，但臨近出口位置處（x=80～100）變化較為雜亂，無規(guī)律可循.此外，為探究孔隙率大小與滲透率k之間的關(guān)系，繪制了孔隙率n=0.36～0.78 與滲透率的關(guān)系，并對其進(jìn)行擬合，得到兩者的關(guān)系表達(dá)式，如圖12 所示.從圖12 中可發(fā)現(xiàn)：計(jì)算滲透率隨模型孔隙率的增大逐漸增大，且近似呈線性關(guān)系，擬合表達(dá)式為y=0.148x-0.012 2，相關(guān)系數(shù)為0.938 4.

圖12 孔隙率與滲透率的擬合關(guān)系Fig.12 Fitting relationship between porosity and permeability

3.4 土顆粒大小對滲流場的影響

以pc=0.01（大顆粒）、0.05（中顆粒）、0.1（小顆粒），n=0.64，100×100×100格點(diǎn)大小的3D模型為例，分析土顆粒大小變化情況下，滲流場速度分布情況.

圖13 為不同土顆粒大小的滲流模擬結(jié)果.從圖13（a）～圖13（c）中可以看出：當(dāng)分布概率為pc=0.01 時(shí)，滲流場穩(wěn)定后的流速主要分布范圍為0～2×10-3；當(dāng)分布概率為pc=0.05 時(shí)，滲流場穩(wěn)定后的流速主要分布范圍為0～1×10-3；當(dāng)分布概率為pc=0.1時(shí)，滲流場穩(wěn)定后的流速主要分布范圍為0～1×10-3.并且隨著土顆粒粒徑減小，對應(yīng)的速度場分布更為均勻，但是速度呈減小的趨勢.從圖13（d）～圖13（f）中可以看出：當(dāng)土顆粒粒徑越小時(shí)，滲流場流線越密集，相應(yīng)的滲流孔道也越細(xì)窄，滲流流體粒子在土體孔隙中的分布越集中.

圖13 不同土顆粒大小的3D滲流模擬結(jié)果Fig.13 Results of three-dimensional seepage simulation with different soil particle sizes

為分析土體顆粒大小與速度間定量化關(guān)系，繪制不同土顆粒大小截面出口位置滲流速度分布曲線，如圖14 所示.由圖14 可知：當(dāng)分布概率pc=0.01時(shí)，出口位置滲流速度波動(dòng)在0.25×10-3～0.45×10-3區(qū)間，最大波動(dòng)幅度可達(dá)1.18×10-3；當(dāng)分布概率pc=0.05時(shí)，出口位置滲流速度在0.25×10-3～0.40×10-3間波動(dòng)，最大波動(dòng)幅度可達(dá)0.89×10-3；而當(dāng)分布概率pc=0.1 時(shí)，出口位置滲流速度波動(dòng)在0.25×10-3～ 0.35×10-3區(qū)間，最大波動(dòng)幅度為0.88×10-3.分析上述原因是，3D重構(gòu)土體顆粒越小，土體內(nèi)部孔隙分布越均勻，流體滲流過程中的流速分布也更加集中.由上可見，土體中的土顆粒大小對滲流特性具有重要的作用.

圖14 不同土顆粒大小出口滲流速度分布曲線Fig.14 Distribution curve of outlet seepage velocity with different soil particle size

圖15 為不同土顆粒大小與滲透率的關(guān)系曲線.由圖15可知，孔隙率n相同情況下，土體滲透率隨土顆粒尺寸減小而減小，且分布概率pc=0.05 時(shí)，滲透率降低特別顯著，降低幅度達(dá)到42.86%，當(dāng)pc＞0.05時(shí)，滲透率基本穩(wěn)定在0.078上下.

圖15 不同土顆粒大小與滲透率的關(guān)系曲線Fig.15 Relationship curves between soil particle size and permeability

3.5 各向異性對滲流場的影響

對1.3 節(jié)各向異性3D 重構(gòu)模型施加如圖6 所示的邊界條件，并借助LBM進(jìn)行滲流場仿真計(jì)算，模擬結(jié)果如圖16 和圖17 所示.從圖16 中可發(fā)現(xiàn)：各向異性模型孔隙滲流速度均小于各向同性的滲流速度，且整個(gè)過程中格點(diǎn)速度呈逐漸減小趨勢.分析原因是：各向異性孔隙模型重構(gòu)過程中，考慮因素更多，重構(gòu)結(jié)果更加貼近于真實(shí)的土體孔隙結(jié)構(gòu)，孔隙結(jié)根據(jù)能量守恒定律可知，相比各向同性孔隙模型滲構(gòu)復(fù)雜無序，這將對流體滲流過程造成一定的阻礙，滲流速度勢必會(huì)更小.

圖17 3D各向異性模型出口邊界滲流速度分布情況Fig.17 Seepage velocity distribution at exit boundary of 3D anisotropic soil reconstruction model

再觀察圖16（d）可知，各向異性重構(gòu)孔隙模型待滲流穩(wěn)定后，流線分布相互交錯(cuò)，分布均勻.此外，對比各向同性模型流線分布可以明顯發(fā)現(xiàn)，各向異性模型的流線較為稀疏，不集中，且x向主滲流通道不明顯.

圖17 為各向異性3D 孔隙模型出口邊界滲流速度分布情況.從圖17 中可發(fā)現(xiàn)：考慮各向異性的模型出口邊界滲流速度波動(dòng)幅度較大、變化趨勢規(guī)律性不明顯；出口邊界平均滲流速度為0.24×10-3，且存在多處格點(diǎn)速度為0 的區(qū)域（圖17 中格點(diǎn)20～ 30 范圍）.分析其主要原因可能是：相較于各向同性模型，各向異性模型重構(gòu)所得孔隙分布復(fù)雜無序，且孔隙空間結(jié)構(gòu)不均一，具有非均勻性.

綜上，3D各向同/異性模型滲流場的格點(diǎn)速度均存自滲流入口到出口呈逐漸減小的趨勢.但各向異性重構(gòu)孔隙模型滲流場流線稀疏、分布相互交錯(cuò)，在x方向的主滲流通道不明顯，且出口各個(gè)位置的流速變化無明顯規(guī)律，平均流速為0.24×10-3，存在多處流速為0 的區(qū)域.而各向同性重構(gòu)孔隙模型速度分布較為集中，且存在明顯的x向主滲流通道，最大流速出現(xiàn)在孔隙較大處，可達(dá)2.12×10-3.

此外，對比文獻(xiàn)［31］中提出的2D 各向異性模型，其在孔隙連通性好的孔道中形成了明顯的主滲流通道，而本文中3D 各向異性模型在x向主滲流通道不明顯，分析其主要原因可能是：模型尺寸、孔隙率、初始?jí)翰?、各向異性等參?shù)設(shè)置不同，并且3D 模型相較于2D 模型在空間結(jié)構(gòu)上更加復(fù)雜、無序，因此，對流體滲流過程造成更大阻礙.

4 對比分析

周瀟等［19］、申林方等［20］、Jun等［21］、蔡沛辰等［22］、李滔等［23］基于LBM模擬細(xì)微觀尺度下流體在多孔介質(zhì)中滲流過程，表3列出了本文與其他學(xué)者關(guān)于多孔介質(zhì)重構(gòu)模型微細(xì)觀滲流研究結(jié)果對比情況.

表3 多孔介質(zhì)滲流研究結(jié)果對比Tab.3 Comparison of research results of porous media seepage

從表3中可以發(fā)現(xiàn)：

1）研究者對多孔介質(zhì)滲流研究一般都直接給出了研究對象的尺寸，而對于最優(yōu)仿真模型尺寸選取鮮有提及.本文綜合孔隙連通性和滲流穩(wěn)定時(shí)間等因素選取3D 最優(yōu)仿真模型尺寸為100×100×100 格點(diǎn)大小，可見3D 模型尺寸比2D 模型［22］尺寸小1/3左右.

2）目前滲流研究分析對象多為速度場、孔隙率和計(jì)算滲透率，本文則更加深入地研究了3D 孔隙結(jié)構(gòu)與速度間的關(guān)系、土顆粒大小對滲流場影響等，發(fā)現(xiàn)多孔介質(zhì)滲流過程中存在明顯的“指進(jìn)”和優(yōu)先流效應(yīng)；土顆粒越小，速度場分布越均勻，但速度更小.在2D 模型中最大格點(diǎn)速度為3.24×10-2，而3D 模型孔隙滲流最大格點(diǎn)速度僅為2.12×10-3，可見3D 滲流速度明顯小于2D 滲流速度，分析原因可能是：3D 模型孔隙相比于2D 模型孔隙結(jié)構(gòu)更為復(fù)雜，對滲流過程造成一定程度阻礙.

3）在土體細(xì)微觀滲流研究領(lǐng)域，僅有少部分學(xué)者在重構(gòu)土體模型時(shí)考慮各向異性，并且結(jié)論多為定性表述.本文基于3D 各向異性模型分析了不同位置處滲流速度大小分布情況，結(jié)論為定量化表征，并發(fā)現(xiàn)滲流速度u在出口邊界處波動(dòng)幅度較大，平均滲流速度u=0.24×10-3，且存在多處u=0的區(qū)域.

4）重構(gòu)細(xì)觀多孔介質(zhì)模型，已被廣泛用于研究多孔介質(zhì)的孔隙結(jié)構(gòu)及預(yù)測其水力學(xué)特征等方面，但模型能否真實(shí)反映實(shí)際對象的孔隙結(jié)構(gòu)，目前尚無定論［32］.本文在前人基礎(chǔ)上，采用改進(jìn)QSGS 法生成了具有各向異性的3D 重構(gòu)土體孔隙結(jié)構(gòu)模型，更接近真實(shí)原狀土體CT掃描模型，具有較好的適用性.

5 結(jié)論

1）3D 模型隨尺寸擴(kuò)大，孔隙連通程度明顯增大，最大增長率為19.23%，出現(xiàn)在100×100×100格點(diǎn)大小處，繼續(xù)擴(kuò)大模型尺寸發(fā)現(xiàn)連通孔隙率nc增加并不明顯.

2）流體粒子易在孔隙連通性好、孔徑大的區(qū)域形成主滲流通道，且存在“指進(jìn)效應(yīng)”，孔道中間流速最大，可達(dá)格點(diǎn)速度0.002 1，越靠近孔壁滲流速度越小，接近0.

3）相同土顆粒大小下，隨模型孔隙率逐漸增大，滲流流體在孔隙中的分布范圍愈加廣泛，平均滲流速度也越大.3D 模型計(jì)算滲透率隨模型孔隙率的增大逐漸增大，且近似呈線性關(guān)系，擬合表達(dá)式為y=0.148x-0.012 2，相關(guān)系數(shù)為0.938 4.

4）孔隙率n相同的情況下，3D 土體滲透率隨著土顆粒尺寸減小而減小，且分布概率pc=0～0.05 時(shí)滲透率顯著降低，降低幅度達(dá)42.86%，當(dāng)pc＞0.05時(shí)，滲透率基本穩(wěn)定在0.078上下.

5）考慮各向異性模型滲流場的格點(diǎn)速度自滲流入口到出口呈逐漸減小的趨勢，且在出口邊界處的滲流速度波動(dòng)幅度較大、變化趨勢規(guī)律性不明顯，平均滲流速度為0.24×10-3，同時(shí)存在多處格點(diǎn)速度為0的區(qū)域.此外，相較于各向同性模型流線分布，各向異性模型的流線相互交錯(cuò)，較為稀疏、不集中，且主滲流通道不明顯.