多孔介質曲折度對膨潤土襯墊滲透性能的影響

侯娟，滕宇陽，李昊，劉磊

（1.上海大學力學與工程科學學院，上海 200444；2.弗吉尼亞大學工學院，弗吉尼亞州夏洛茨維爾22904；3.中國科學院武漢巖土力學研究所巖土力學與工程國家重點實驗室，湖北武漢 430071）

膨潤土襯墊（Geosynthetic Clay Liner，GCL）是由兩層土工織物間夾裹一層薄的膨潤土構成的，近幾十年來廣泛應用于垃圾填埋場的頂部封場及底部襯墊系統，具有防滲性能好、厚度小以及自愈力強等優點.目前對于GCL 滲透系數的研究，主要集中在宏、微觀試驗研究方面［1-2］，Shackelford 等［3］對GCL 進行了大量滲透試驗，發現影響GCL 滲透性的因素有孔隙比、膨潤土級配、厚度等.劉志彬等［4］對膨潤土內部微孔隙進行了定量分析，并與掃描電鏡圖像進行了對比.王寶等［5］通過試驗發現孔隙比與滲透系數存在良好的線性關系，該結論與Kang 等［6］用去離子水滲透GCL 結果相同.Mesri 和Olson［7］發現，滲透系數和孔隙比在對數坐標系下存在線性關系.Fratalocchi［8］研究了應力對GCL 滲透系數的影響，發現滲透系數的對數與孔隙比之間存在較好的線性關系，得到與Petrov 等［9］相似的結論.何俊等［10］通過微觀結構分析，推導了GCL 的膨潤土飽和滲透系數與孔隙比之間的計算公式.

然而，僅用孔隙率往往存在不能充分體現土體內部粒間孔隙特征等問題，也不能反映顆粒尺寸本身對滲透系數的影響［11］.自Kozeny［12］提出流徑曲折度概念以來，Carman［13］考慮孔隙通道的不規則性，在Poiseuille 定律的基礎上，結合Kozeny 公式得到了曲折度與滲透系數之間的表達式.Nooruddin 和Hossain［14］根據曲折度與土顆粒之間的關系，改進了Kozeny-Carman 模型.Tsang［15］研究了滲流路徑的曲折度對流場的影響.Murata 等［16］基于孔隙介質的滲透系數Kozeny-Carman 方程，推導得到了修正的流量關系表達式.此外，還有眾多學者采用試驗研究［17］、理論分析［10，18］以及數值模擬［4］等手段研究了多孔介質孔隙流徑的曲折度與孔隙率之間的關系.

GCL 中起主要防滲作用的膨潤土，是一種顆粒遇水會膨脹的特殊多孔介質.然而，盡管目前已有大量研究表明GCL中膨潤土孔隙率與滲透系數呈正相關的關系［19］，但是，如何從細觀角度考慮膨潤土顆粒尺寸及孔隙特征對GCL 滲透性能的影響，成為進一步研究和解釋膨潤土滲透機理的難點［9，20］.GCL 中曲折度是影響滲透系數的關鍵因素之一［21］，但對其膨脹過程中量化分析的研究比較鮮見，目前尤其缺乏對孔隙滲流通道的曲折度與孔隙率之間關系的研究.隨著計算機技術的高速發展，數值模擬成為目前研究多孔介質曲折度的主要手段之一［22］.其中，COMSOL 在分析多物理場耦合以及在多孔介質中流體流動等復雜問題方面具有明顯的優勢［23］，因此，本文采用COMSOL 軟件，構建了包含孔隙率與曲折度等細觀變量的多孔介質GCL結構模型，從細觀角度，研究土體的顆粒尺寸及孔隙特征等對GCL滲透性能的影響.首先，通過構建流體物理場和粒子追蹤場，模擬分析蠕動流場的穩態速度圖和壓力矢量圖，得到宏觀流場分布規律；然后，通過對比宏觀蠕動流場與微觀流動（水）粒子分布規律，探究宏觀流場分布與微觀流動粒子之間的關系；同時，通過實時追蹤液體粒子的運移過程，得到液體在多孔介質內流徑并通過數學變換得到曲折度與流徑曲折度的表達式；最后，采用毛管模型描述粒間孔隙組成的土骨架［24］，并基于此提出能夠綜合考慮孔隙率與流徑曲折度的GCL滲透系數理論預測模型，驗證模型的準確性.

1 模型建立和參數確定

COMSOL 以有限元為基礎，通過求解偏微分方程（單場）和方程組（多場）實現多場耦合的仿真模擬［25］.建模的過程一般分為模型假設、確定物理場控制方程、幾何模型建立、邊界與初始條件確定、網格劃分以及最終求解6個部分.

1.1 模型假設

參考張志紅等［26］的研究，本文基本假設為：多孔介質內部為均一幾何孔隙；過程等溫；滲流液均質各向同性；忽略毛細管吸力作用；流體不可壓縮且為單相穩定滲流.

1.2 物理場控制方程

1）流體物理場

由于GCL 的滲透系數較小，在膨潤土粒間孔隙通道中，形成斯托克斯流（Stokes Flow），由Navier-Stokes連續性方程得蠕動流方程：

式中：v為動力黏度；ρ為液體密度；P為流體壓力；u為流體的速度矢量；?為拉普拉斯算符.

2）粒子追蹤場

采用了歐拉-拉格朗日法對運動顆粒進行追蹤，同時采用牛頓第二定律描述運動過程.歐拉-拉格朗日法采用式（2）記錄每個粒子的不同狀態.

式中：mp為顆粒質量；vp為顆粒速度；F為顆粒受力.

1.3 模型驗證

首先采用宏觀試驗對COMSOL 數值模型進行驗證.依據Petrov 等［9］試驗結果，建立COMSOL 數值模型，并采用穩態對蠕動流場進行模擬，得到滲透系數分布圖如圖1 所示.計算得到滲流出口的滲透系數k為1.57×10-11m/s，與試驗滲透系數比值為1.121.同時，分別建立其他試驗結果的COMSOL 數值模型，得到的滲透系數與試驗對比匯總如表1所示.由表1可知，該COMSOL 數值模型能夠較好地模擬GCL 在去離子水滲透時的工況.其次采用宏觀流場對COMSOL 數值模型進行驗證.在本文的仿真模擬過程中，宏觀流場（圖1 所示）的規律與金磊等［27］的研究結果類似，即在孔隙通道半徑大處流速增加明顯.同時，本文流場速度隨孔隙通道半徑增大而增大，其結果與梁越等［28］的研究結論一致.

圖1 滲透系數分布圖（n=0.4）Fig.1 Distribution of the hydraulic conductive（n=0.4）

表1 模型準確性對比分析Tab.1 Comparative analysis of the model accuracy

1.4 幾何模型建立

1.3 節驗證了COMSOL 中的蠕動流場可對GCL的層流運動進行宏觀模擬.由于在孔隙率不變時，等效滲透系數與孔隙尺寸的平方成正比［30-31］，考慮到有限元計算的時效性，采用幾何相似性理論，將驗證模型尺寸進行適當擴大，進行下列計算.同時，也有研究發現顆粒的幾何外形［32］對其滲透系數的影響并不顯著，因此，本模型構建邊長為10 mm 的大正方形表示GCL試樣，交錯布置9×9邊長為1 mm的小正方形表示膨潤土顆粒，在孔隙率n=0～1范圍內進行優化計算.

已有學者［33-35］研究發現，在膨潤土膨脹過程中的晶體膨脹階段，是膨潤土固相顆粒的膨脹間距不斷增大的膨脹過程，即顆粒幾何尺寸擴大而孔隙率減小的過程.同時采用正方形的放大命令，確保模型膨脹前后幾何外形一致，也可保證模型達到所需的最小孔隙率，實現模擬孔隙率n=0～1 的全區域變化過程.因此通過COMSOL 內逐漸放大小正方形的命令模擬顆粒膨脹過程，從而控制幾何模型的孔隙率逐漸變小，部分工況如圖2 所示.選用蠕動流模塊和流動粒子追蹤模塊進行瞬時計算，通過實時監測流動粒子的位置坐標，最終確定液體在多孔介質內流通的滲流路徑.

圖2 不同孔隙率下的幾何模型Fig.2 The geometric models for different porosities

1.5 邊界與初始條件確定

蠕動流模塊邊界條件設定為小正方形（顆粒）外圍及大正方形（膨潤土試樣）上下邊界均為封閉邊界，流體受壓力P作用下，從左邊界流入，右邊界單向流出.

流動粒子追蹤模塊邊界條件為粒子均勻排布在左邊界，由蠕動流中的速度場控制粒子運動至右側壁，最后附著在右側壁上后其速度為0.上下邊界、小正方形外圍邊界均為保持粒子動量守恒的粒子反彈壁，使得粒子在接觸到上述邊界后立即反彈，速度大小不變，不發生動能損失.整體COMSOL模型邊界條件如圖3所示.

圖3 GCL簡化幾何2D模型Fig.3 The simplified geometric 2D model of the GCL

工程實際中，垃圾填埋場中GCL 襯墊上部滲濾液收集系統及第一層垃圾體上覆荷載一般約為14 kPa［36-37］.因此，模型取入口邊界壓力為P=14 kPa，該條件同時滿足蠕動流場雷諾數小于1 的要求.模型參數詳見表2.

表2 模型參數確定Tab.2 Setting parameters of the model

1.6 網格劃分

考慮時效性，參考Zimmerman［25］的研究，采用三角形網格進行劃分計算，并分別對n=0～1 之間不同孔隙率的多孔介質模型進行模擬.

值得一提的是，在GCL滲透實驗中，測得GCL試樣的初始孔隙率一般為0.6［5，37］.但當垃圾填埋場開始服役時，GCL 襯墊在吸水膨脹與上覆荷載作用下，孔隙率的數值將明顯減?。?8-39］.同時，Seiphoori 等［40］對MX-80 膨潤土試樣進行膨脹試驗后表明，不同壓實狀態下膨潤土的孔隙率介于0.34～0.58 之間，因此，以下主要以孔隙率n=0.41 的工況為例進行分析計算，具體網格劃分如圖4所示.

圖4 網格劃分Fig.4 Mesh division

2 仿真結果與分析

2.1 蠕動流計算結果

基于COMSOL 軟件自身的求解器的特性，一般認為穩態求解可用以驗證瞬時求解的結果是否趨于穩定，而瞬時求解有利于分析宏觀試驗現象隨時間的變化趨勢，因此，將蠕動流場和粒子追蹤模塊進行耦合瞬時計算.

首先，采用穩態對蠕動流場進行模擬.入口邊界選取充分發展的流動命令以減少邊界效應，即當流體到達左邊界（入口）時，水流趨于平穩，這樣可以有效地減少因為邊界加壓而導致的入口處水流速度不均勻的現象.圖5 所示為模擬試樣中流體的速度圖與壓力矢量圖.在速度矢量圖（圖5（a））中，可觀測到A、B、D 點處路徑的最大速度小于C 點處最大速度.分析原因是C 點的水流速度源自上、下兩側豎向通道，類似于通道半徑在此處增加.因此，根據泊肅葉公式，通道半徑增加會導致水流速度增加，即通道個數直接影響流速.同時，由壓力矢量圖（圖5（b））可見，粒間孔隙壓力從左到右（出口）出現均勻遞減的趨勢.此外，圖5（b）左邊界下部壓力分布（虛線）比上部壓力分布（實線）更緊密，主要原因是幾何模型孔隙通道布置不均勻，導致壓力在上、下邊界上出現不均勻分布.

2.2 流動粒子追蹤

在瞬時狀態下，將蠕動流場和粒子追蹤模塊進行耦合計算.模擬過程中發現粒子追蹤場在20 000 μs時基本達到穩定流的速度狀態，因此，僅監測流動粒子在t=0～20 000 μs 內的位置坐標.圖6 給出了部分時刻流動粒子軌跡圖.

對比圖5（a）和圖6（a）可知，微觀流動粒子分布呈現與宏觀流體場分布類似的規律.圖6（a）中，E、F點路徑內的最大速度小于G 點的最大速度（實線所示），分析其原因是G 點等效通道半徑增加，導致速度增加；在相同的原因下，導致H 點的最大速度小于I 點最大速度（虛線所示）.因此，圖6（a）中宏觀蠕動流場是細觀粒子流動場作用的結果.

圖5 試樣中流體速度圖與壓力矢量圖Fig.5 Fluid velocity and pressure vectors in the specimen

但同時對比發現，圖6（a）中蠕動流場平均速度為1.53 m/s，因此，如果不考慮曲折度影響，粒子按直線規律運移的話，應該在t=6 535 μs 時到達右邊界.然而，如圖6（d）所示，實際粒子運移中，當t=6 535 μs時，其僅運移到距右邊界一定位置處（圖6（d）中虛線所示），而在t=8 000 μs（約為整體運移時間的前三分之一）時，大部分粒子才到達右邊界，少數粒子甚至需要比這更長的時間（如圖6（e）所示）.這是由于流徑具有一定的曲折度，使得粒子實際運移路徑大于直線距離.此外，觀察圖6（b）～（e）（粒子運移過程）可發現，在流體速度控制下，流動粒子先從左邊界均勻分布，通過粒間孔隙通道，通道較窄，速度增加，隨后發生分流，在上、下側通道內速度減緩；繼續通過粒間孔隙通道，通道收窄，速度再次增加.依次反復后，流動粒子呈不均勻狀黏附在出口處.

圖6 流動粒子軌跡圖Fig.6 The trajectory of the particle flowing

2.3 GCL滲透系數理論預測模型

GCL 襯墊中起主要防滲作用的膨潤土是一種多孔介質.Kozeny［12］和Carman［13］提出了預測多孔介質滲透率的半經驗公式（Kozeny-Carman（K-C）方程），該式適用于多孔介質滲透系數較低的情況.首先，Nooruddin和Hossain［14］在K-C方程基礎上，將土體通道簡化為毛管模型［41］（圖7）.長度為L、截面面積為A的多孔介質土體（圖7（a）），此時的固相顆粒為不規則分布，由孔隙組成的滲流通道也不均勻分布，均質土體的固相顆粒一般為土顆粒；采用等流量條件下滲流通道的等效毛管化，將土顆粒規則排布后，形成土顆粒以外的孔隙通道，等效后為直徑2R的通道（圖7（b））；但由于土體的實際滲流路徑如圖7（c）中實線所示，即在土顆粒與土顆粒之間形成曲折的路徑；假設該路徑中孔隙通道直徑為2r（r≤R），如圖7（d）所示，以此形成等效毛管模型來表示液體在多孔介質中的移動規律.

圖7 滲透通道等效毛管化分析示意圖Fig.7 Schematic of equivalent capillary of permeation channels

同時，由于土體的滲流路徑L具有一定的曲折性.Collins［41］采用曲折度τ來反映多孔介質中平均有效流動路徑長度Li（m）與流體流動方向的直線路徑L（m）的比值，如公式（3）所示.不少學者［42-47］對曲折度τ進行進一步研究后發現，孔隙率與曲折度的關系可通過數學表達式進行表述，如表3所示.

表3 孔隙率與曲折度關系Tab.3 Relationship between porosity and tortuosity

由于膨潤土具有自相似性，因此可將土顆粒通過單位體積法簡化為單元體.同時，將簡化后的單元體表示為毛管模型，以便通過滲流路徑獲得滲流的等效理論模型.本文等效模型作如下假設：假設在等效滲透通道模型中，總共有M根半徑為ri、長度分別為Li的毛管束.結合修正的哈根-泊肅葉［48-49］方程可知，單位時間內通過M根毛管的滲流速度q（m3/s）為：

式中：μ為滲濾液的動力黏度，Pa·s；ρw為液體的密度，kg/m3；g為重力加速度，m/s2；π=3.141 59；ri為在M根彎曲毛管中，第i根毛管的通道半徑，m；Li為第i根毛管的滲流長度，m.

將式（3）代入式（4），得通過土體的總流量Q（m3/s）為：

又由達西定律可知，水流通量J（m/s）與水力梯度ΔH/L成正比：

式中：k為滲透系數，m/s；ΔH為水力梯度；L為滲流路徑，m.

由于水流通量J是總流量Q與土體截面面積A（m2）的比值（J=Q/A），將式（6）代入式（5），換算可得滲透系數k為：

值得一提的是，盡管式（7）中的曲折度τ是基于多孔介質本身性質的參數，其難于用試驗手段直接測得，但可以通過孔隙率n計算獲得（見表3）.

因此，本文孔隙率n從0 到1 的一系列COMSOL模型，采用流動粒子追蹤方法獲得孔隙率全范圍內流動粒子的運動軌跡，并統計計算相應的滲流路徑以及曲折度，模型選取的試驗驗證參數如表4 所示.由此計算得到不同孔隙率與曲折度之間的關系（圖8中散點標示），可以看出，曲折度隨著孔隙率的增加而減小.

表4 試驗驗證參數［50］Tab.4 The experimental parameters verification

圖8 孔隙率與曲折度擬合圖Fig.8 The fitting curve of porosity and tortuosity

由擬合模擬的一系列工況所得結果（圖8 中實線標示），可得到孔隙率與曲折度之間的數學表達式為：

將式（8）代入式（7）中，即可得到能綜合反映孔隙率與流徑曲折度影響的GCL 滲透系數理論預測模型：

為驗證式（9）的準確性，本文采用大量已有試驗結果（Chai 等［50］、Jo 等［51-52］、Petrov 和Rowe［53］、Chen等［54］、Setz等［55］、Bradshaw和Benson［56］以及Seiphoor［57］）進行計算比較，結果如圖9 所示.以Chai 等［50］數據為例，具體的計算過程如下.

圖9 GCL預測滲透系數與試驗滲透系數Fig.9 Predicted and measured hydraulic conductivities of the GCL

由式（9）可得此時GCL的預測滲透系數k：

此時，所得預測結果和試驗結果比值為：

本文采用對數坐標系下滲透系數偏差分析的表示方法［58-61］（如圖9 所示）.圖中，橫坐標為預測滲透系數值，縱坐標為試驗滲透系數值.從圖中可看出，預測值與試驗值的范圍為0.23～3.86.按照行業標準《鈉基膨潤土復合防水襯墊》（FZ/T 64036—2103），滲透系數值的誤差常在一個數量級內進行表征，如規定滲透系數≤5×10-9～1×10-10cm/s.同時，對數坐標系下滲透系數的偏差一般介于1/10～10［37，58-61］.本文GCL 的理論預測模型計算結果與試驗值的比在1/5～5，小于一個數量級，說明該GCL 理論預測模型可以較為準確地用于預測GCL 的滲透系數，可以指導工程實際.

3 結論與展望

本文通過COMSOL 模擬了GCL在膨脹變形過程中，曲折度隨孔隙率的變化而變化的規律，得到多孔介質滲流流徑的細觀模型，研究了孔隙率和流徑曲折度對GCL滲透系數的影響，主要結論如下：

1）當入口邊界條件為均勻壓力時，在蠕動流場中，流體通道半徑的增加會導致流體速度增加，同時，流體壓力從入口到出口處呈現均勻遞減的趨勢.

2）GCL 試樣中，微觀流動粒子分布與宏觀流體場分布規律基本一致，且大部分流動粒子在整體時間的前1/3段完成運移.

3）GCL 試樣中，流徑的曲折度隨孔隙率的增加而減小，且呈指數函數關系.

4）提出的能綜合反映孔隙率與流徑曲折度影響的GCL 滲透系數理論預測模型，可以較為準確地預測試驗結果，兩者比值介于1/5～5.

需要注意的是，本文膨潤土曲折度與孔隙率之間的關系主要是基于模擬均一幾何尺寸和均一幾何孔隙表征膨潤土試樣得到的.因此，研究結果可以從概念上定性地分析曲折度隨膨潤土顆粒膨脹的變化規律，同時反映曲折度對GCL 滲透系數的影響.但是，后續還需要對大量不規則試樣展開進一步的深入和系統的研究，才有可能將研究結果推廣到定量分析以及實際應用中.