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

（福州大學 土木工程學院，福建 福州，350108）

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

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

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

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

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

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

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

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

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

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

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

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

2）在3D 空間中，令分布的生長相單元以一定的生長概率pd朝19個運動方向生長，如圖1所示.

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

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

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

1.3 改進四參數(shù)隨機生長法

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

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

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

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

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

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

2 格子 Boltzmann 模型

2.1 格子 Boltzmann 方程

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

式中：Fα（ω，t）為格點ω位置處在t時刻沿α方向的粒子分布函數(shù)；eα為離散速度；δt為離散時間；τ為弛豫時間；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 個運動方向，離散速度eα滿足式（3）.

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

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

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

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

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

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

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

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

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

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

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

2.4 邊界條件處理

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

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

2.5 收斂性判別

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

2.6 LBM算例驗證

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

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

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

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

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

3 結(jié)果與討論

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

圖6 計算模型邊界條件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ù)值仿真最優(yōu)模型尺寸的大小.滲流穩(wěn)定狀態(tài)判斷標準［19］：當格點離散速度不再隨時步增加而發(fā)生變化時，可認為滲流已經(jīng)達到穩(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)計3D 模型的孔隙參數(shù)，統(tǒng)計結(jié)果如表2所示.從表2中可看出：隨模型尺寸擴大，孔隙連通程度明顯增大，模型尺寸擴大至100×100×100 格點大小后，繼續(xù)擴大模型尺寸時發(fā)現(xiàn)連通孔隙率nc增加并不明顯，再結(jié)合圖7（b）分析可知，隨模型尺寸增大滲流穩(wěn)定所需時間呈直線形增大.鑒于此，綜合孔隙連通性和運算時間因素，本文最終選取100×100×100格點大小的多孔介質(zhì)模型作為最優(yōu)仿真單元模型，用于后續(xù)研究.該研究方法可為其他類似3D巖土體最優(yōu)仿真單元選取提供一定的借鑒.

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

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

圖8 3D模型尺寸對滲流時間和速度的影響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格點大小的模型為例進行滲流模擬，探討孔隙結(jié)構(gòu)特征對滲流過程的影響.

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

圖9 3D滲流模擬計算結(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 格點大小的模型為例，分析孔隙率變化情況下，滲流穩(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)系，并對其進行擬合，得到兩者的關(guān)系表達式，如圖12 所示.從圖12 中可發(fā)現(xiàn)：計算滲透率隨模型孔隙率的增大逐漸增大，且近似呈線性關(guān)系，擬合表達式為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格點大小的3D模型為例，分析土顆粒大小變化情況下，滲流場速度分布情況.

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

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

為分析土體顆粒大小與速度間定量化關(guān)系，繪制不同土顆粒大小截面出口位置滲流速度分布曲線，如圖14 所示.由圖14 可知：當分布概率pc=0.01時，出口位置滲流速度波動在0.25×10-3～0.45×10-3區(qū)間，最大波動幅度可達1.18×10-3；當分布概率pc=0.05時，出口位置滲流速度在0.25×10-3～0.40×10-3間波動，最大波動幅度可達0.89×10-3；而當分布概率pc=0.1 時，出口位置滲流速度波動在0.25×10-3～ 0.35×10-3區(qū)間，最大波動幅度為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 時，滲透率降低特別顯著，降低幅度達到42.86%，當pc＞0.05時，滲透率基本穩(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進行滲流場仿真計算，模擬結(jié)果如圖16 和圖17 所示.從圖16 中可發(fā)現(xiàn)：各向異性模型孔隙滲流速度均小于各向同性的滲流速度，且整個過程中格點速度呈逐漸減小趨勢.分析原因是：各向異性孔隙模型重構(gòu)過程中，考慮因素更多，重構(gòu)結(jié)果更加貼近于真實的土體孔隙結(jié)構(gòu)，孔隙結(jié)根據(jù)能量守恒定律可知，相比各向同性孔隙模型滲構(gòu)復雜無序，這將對流體滲流過程造成一定的阻礙，滲流速度勢必會更小.

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

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

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

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

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

4 對比分析

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

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

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