多孔介質(zhì)曲折度對(duì)膨潤(rùn)土襯墊滲透性能的影響

侯娟 滕宇陽(yáng) 李昊 劉磊

摘要：膨潤(rùn)土的顆粒膨脹所引起的孔隙率及流徑曲折度的改變是影響膨潤(rùn)土襯墊（Geo? synthetic Clay Liner，GCL）防滲性能的關(guān)鍵.通過(guò) COMSOL 建立模擬膨潤(rùn)土中滲流的數(shù)值模型，研究了膨潤(rùn)土顆粒膨脹過(guò)程對(duì) GCL 防滲性能的影響. COMSOL模型研究結(jié)果表明，流體通過(guò)膨潤(rùn)土粒間孔隙時(shí)，邊界通道的最大速度明顯小于非邊界通道速度，這與泊肅葉公式的基本結(jié)論一致，初步驗(yàn)證了模型的準(zhǔn)確性.同時(shí)，微觀粒子速度分布規(guī)律與宏觀流體場(chǎng)的速度分布規(guī)律基本一致，因此可以從微觀的角度研究GCL 的宏觀滲透規(guī)律.利用 COMSOL 中的粒子軌跡監(jiān)測(cè)方法，分析了顆粒膨脹對(duì)膨潤(rùn)土孔隙率和流徑曲折度的影響.研究結(jié)果表明，流徑的曲折度隨孔隙率的增加而減小，且呈指數(shù)函數(shù)關(guān)系.基于 GCL 試樣的多孔介質(zhì)毛管模型和修正泊肅葉公式，考慮 GCL 孔隙率與流徑曲折度之間的指數(shù)函數(shù)關(guān)系，提出了能綜合反映孔隙率與流徑曲折度影響的 GCL滲透系數(shù)理論預(yù)測(cè)模型，并將理論預(yù)測(cè)結(jié)果與已有文獻(xiàn)試驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比，兩者的比值介于1/5～5，驗(yàn)證了模型的準(zhǔn)確性.

關(guān)鍵詞：膨潤(rùn)土襯墊;滲透系數(shù);COMSOL;曲折度;理論模型

中圖分類號(hào)：TU411.4? 文獻(xiàn)標(biāo)志碼：A

Influence of Porosity and Tortuosity on Hydraulic Conductive of Geosynthetic Clay Liner

HOU Juan1，2，3?，TENG Yuyang1，LI Hao1，LIU Lei3

（1. School of Mechanics and Engineering Sciences，Shanghai University，Shanghai 200444，China;

2. School of Engineering，University of Virginia，Charlottesville，VA 22904，USA;

3. State Key Laboratory of Geomechanics and Geotechnical Engineering，Institute ofRock and Soil Mechanics，Chinese Academy of Sciences，Wuhan 430071，China）

Abstract：The change of porosity and the tortuosity induced by the swelling progress of Geosynthetic Clay Liner （GCL） influence the hydraulic conductive obviously. The model simulating the seepage in bentonite is established by COMSOL to investigate the effect of the expansion process of bentonite particles on the hydraulic conductive of GCL.The COMSOL simulation results showed that the maximum velocity of the boundary channel is less than that of the non-boundary channel， when the liquid passes through the pores between the bentonite particles， which is consistent with the Poiseuille equation and initially verifies the accuracy of the model. The microscopic particle velocity distri ? bution law is consistent with the velocity distribution law of macroscopic fluid field. Therefore， the macroscopic hy? draulic law can be investigated in the microscopic aspect. Meanwhile， the particle expansion on the porosity and the tortuosity of bentonite is analyzed by using the particle trajectory monitoring method in COMSOL. The results showed that the tortuosity decreases with the increase of porosity， and it is exponentially functional. Based on the capillary model of porous media and the Poiseuille Formula， the theoretical prediction model of the hydraulic conductive of GCL， which can reflect the influence of porosity and the tortuosity， is proposed， considering the exponential function relationship between the porosity and the tortuosity. The theoretical calculation prediction model for predicting the hydraulic conductive of GCL is proposed. The ratio of the theoretical and experiments is between 1/5 and 5， which confirms the accuracy of the theoretical model of GCL.

Key words：Geosynthetic Clay Liner（ GCL）;hydraulic conductivity;COMSOL;tortuosity;theoretical model

膨潤(rùn)土襯墊（Geosynthetic Clay Liner，GCL）是由兩層土工織物間夾裹一層薄的膨潤(rùn)土構(gòu)成的，近幾十年來(lái)廣泛應(yīng)用于垃圾填埋場(chǎng)的頂部封場(chǎng)及底部襯墊系統(tǒng)，具有防滲性能好、厚度小以及自愈力強(qiáng)等優(yōu)點(diǎn).目前對(duì)于 GCL滲透系數(shù)的研究，主要集中在宏、微觀試驗(yàn)研究方面[1-2]，Shackelford 等[3]對(duì) GCL 進(jìn)行了大量滲透試驗(yàn)，發(fā)現(xiàn)影響 GCL 滲透性的因素有孔隙比、膨潤(rùn)土級(jí)配、厚度等.劉志彬等[4]對(duì)膨潤(rùn)土內(nèi)部微孔隙進(jìn)行了定量分析，并與掃描電鏡圖像進(jìn)行了對(duì)比.王寶等[5]通過(guò)試驗(yàn)發(fā)現(xiàn)孔隙比與滲透系數(shù)存在良好的線性關(guān)系，該結(jié)論與 Kang 等[6]用去離子水滲透 GCL結(jié)果相同. Mesri和 Olson[7]發(fā)現(xiàn)，滲透系數(shù)和孔隙比在對(duì)數(shù)坐標(biāo)系下存在線性關(guān)系. Frataloc? chi[8]研究了應(yīng)力對(duì) GCL 滲透系數(shù)的影響，發(fā)現(xiàn)滲透系數(shù)的對(duì)數(shù)與孔隙比之間存在較好的線性關(guān)系，得到與 Petrov 等[9]相似的結(jié)論.何俊等[10]通過(guò)微觀結(jié)構(gòu)分析，推導(dǎo)了 GCL 的膨潤(rùn)土飽和滲透系數(shù)與孔隙比之間的計(jì)算公式.

然而，僅用孔隙率往往存在不能充分體現(xiàn)土體內(nèi)部粒間孔隙特征等問(wèn)題，也不能反映顆粒尺寸本身對(duì)滲透系數(shù)的影響[11].自Kozeny[12]提出流徑曲折度概念以來(lái)，Carman[13]考慮孔隙通道的不規(guī)則性，在 Poiseuille 定律的基礎(chǔ)上，結(jié)合Kozeny公式得到了曲折度與滲透系數(shù)之間的表達(dá)式. Nooruddin 和 Hos? sain[14]根據(jù)曲折度與土顆粒之間的關(guān)系，改進(jìn)了Kozeny-Carman 模型. Tsang[15]研究了滲流路徑的曲折度對(duì)流場(chǎng)的影響. Murata 等[16]基于孔隙介質(zhì)的滲透系數(shù)Kozeny-Carman 方程，推導(dǎo)得到了修正的流量關(guān)系表達(dá)式.此外，還有眾多學(xué)者采用試驗(yàn)研究[17]、理論分析[10，18]以及數(shù)值模擬[4]等手段研究了多孔介質(zhì)孔隙流徑的曲折度與孔隙率之間的關(guān)系.

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

1模型建立和參數(shù)確定

COMSOL 以有限元為基礎(chǔ)，通過(guò)求解偏微分方程（單場(chǎng)）和方程組（多場(chǎng)）實(shí)現(xiàn)多場(chǎng)耦合的仿真模擬[25] .建模的過(guò)程一般分為模型假設(shè)、確定物理場(chǎng)控制方程、幾何模型建立、邊界與初始條件確定、網(wǎng)格劃分以及最終求解6個(gè)部分.

1.1模型假設(shè)

參考張志紅等[26]的研究，本文基本假設(shè)為：多孔介質(zhì)內(nèi)部為均一幾何孔隙;過(guò)程等溫;滲流液均質(zhì)各向同性;忽略毛細(xì)管吸力作用;流體不可壓縮且為單相穩(wěn)定滲流.

1.2物理場(chǎng)控制方程1）流體物理場(chǎng)

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

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

2）粒子追蹤場(chǎng)

采用了歐拉-拉格朗日法對(duì)運(yùn)動(dòng)顆粒進(jìn)行追蹤，同時(shí)采用牛頓第二定律描述運(yùn)動(dòng)過(guò)程.歐拉-拉格朗日法采用式（2）記錄每個(gè)粒子的不同狀態(tài).

式中：mp為顆粒質(zhì)量;vp為顆粒速度;F 為顆粒受力.

1.3模型驗(yàn)證

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

1.4幾何模型建立

1.3節(jié)驗(yàn)證了 COMSOL 中的蠕動(dòng)流場(chǎng)可對(duì) GCL 的層流運(yùn)動(dòng)進(jìn)行宏觀模擬.由于在孔隙率不變時(shí)，等效滲透系數(shù)與孔隙尺寸的平方成正比[30-31]，考慮到有限元計(jì)算的時(shí)效性，采用幾何相似性理論，將驗(yàn)證模型尺寸進(jìn)行適當(dāng)擴(kuò)大，進(jìn)行下列計(jì)算.同時(shí)，也有研究發(fā)現(xiàn)顆粒的幾何外形[32]對(duì)其滲透系數(shù)的影響并不顯著，因此，本模型構(gòu)建邊長(zhǎng)為10 mm 的大正方形表示 GCL試樣，交錯(cuò)布置9×9邊長(zhǎng)為1 mm 的小正方形表示膨潤(rùn)土顆粒，在孔隙率 n=0～1范圍內(nèi)進(jìn)行優(yōu)化計(jì)算.

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

1.5 邊界與初始條件確定

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

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

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

1.6 網(wǎng)格劃分

考慮時(shí)效性，參考 Zimmerman[25]的研究，采用三角形網(wǎng)格進(jìn)行劃分計(jì)算，并分別對(duì) n=0～1之間不同孔隙率的多孔介質(zhì)模型進(jìn)行模擬.

值得一提的是，在 GCL 滲透實(shí)驗(yàn)中，測(cè)得 GCL 試樣的初始孔隙率一般為0.6[5，37].但當(dāng)垃圾填埋場(chǎng)開(kāi)始服役時(shí)，GCL 襯墊在吸水膨脹與上覆荷載作用下，孔隙率的數(shù)值將明顯減小[38-39].同時(shí)，Seiphoori等[40] 對(duì) MX-80膨潤(rùn)土試樣進(jìn)行膨脹試驗(yàn)后表明，不同壓實(shí)狀態(tài)下膨潤(rùn)土的孔隙率介于0.34～0.58之間，因此，以下主要以孔隙率 n=0.41的工況為例進(jìn)行分析計(jì)算，具體網(wǎng)格劃分如圖4所示.

2仿真結(jié)果與分析

2.1 蠕動(dòng)流計(jì)算結(jié)果

基于 COMSOL 軟件自身的求解器的特性，一般認(rèn)為穩(wěn)態(tài)求解可用以驗(yàn)證瞬時(shí)求解的結(jié)果是否趨于穩(wěn)定，而瞬時(shí)求解有利于分析宏觀試驗(yàn)現(xiàn)象隨時(shí)間的變化趨勢(shì)，因此，將蠕動(dòng)流場(chǎng)和粒子追蹤模塊進(jìn)行耦合瞬時(shí)計(jì)算.

首先，采用穩(wěn)態(tài)對(duì)蠕動(dòng)流場(chǎng)進(jìn)行模擬.入口邊界選取充分發(fā)展的流動(dòng)命令以減少邊界效應(yīng)，即當(dāng)流體到達(dá)左邊界（入口）時(shí)，水流趨于平穩(wěn)，這樣可以有效地減少因?yàn)檫吔缂訅憾鴮?dǎo)致的入口處水流速度不均勻的現(xiàn)象.圖5所示為模擬試樣中流體的速度圖與壓力矢量圖.在速度矢量圖（圖5（a））中，可觀測(cè)到 A、B、D 點(diǎn)處路徑的最大速度小于 C 點(diǎn)處最大速度.分析原因是 C點(diǎn)的水流速度源自上、下兩側(cè)豎向通道，類似于通道半徑在此處增加.因此，根據(jù)泊肅葉公式，通道半徑增加會(huì)導(dǎo)致水流速度增加，即通道個(gè)數(shù)直接影響流速.同時(shí)，由壓力矢量圖（圖5（b））可見(jiàn)，粒間孔隙壓力從左到右（出口）出現(xiàn)均勻遞減的趨勢(shì).此外，圖5（b）左邊界下部壓力分布（虛線）比上部壓力分布（實(shí)線）更緊密，主要原因是幾何模型孔隙通道布置不均勻，導(dǎo)致壓力在上、下邊界上出現(xiàn)不均勻分布.

2.2 流動(dòng)粒子追蹤

在瞬時(shí)狀態(tài)下，將蠕動(dòng)流場(chǎng)和粒子追蹤模塊進(jìn)行耦合計(jì)算.模擬過(guò)程中發(fā)現(xiàn)粒子追蹤場(chǎng)在20000μs 時(shí)基本達(dá)到穩(wěn)定流的速度狀態(tài)，因此，僅監(jiān)測(cè)流動(dòng)粒子在 t=0～20000μs 內(nèi)的位置坐標(biāo).圖6給出了部分時(shí)刻流動(dòng)粒子軌跡圖.

對(duì)比圖5（a）和圖6（a）可知，微觀流動(dòng)粒子分布呈現(xiàn)與宏觀流體場(chǎng)分布類似的規(guī)律.圖6（a）中，E、F 點(diǎn)路徑內(nèi)的最大速度小于 G 點(diǎn)的最大速度（實(shí)線所示），分析其原因是 G點(diǎn)等效通道半徑增加，導(dǎo)致速度增加;在相同的原因下，導(dǎo)致 H點(diǎn)的最大速度小于 I 點(diǎn)最大速度（虛線所示）.因此，圖6（a）中宏觀蠕動(dòng)流場(chǎng)是細(xì)觀粒子流動(dòng)場(chǎng)作用的結(jié)果.

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

2.3 GCL 滲透系數(shù)理論預(yù)測(cè)模型

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

同時(shí)，由于土體的滲流路徑 L 具有一定的曲折性. Collins[41]采用曲折度τ來(lái)反映多孔介質(zhì)中平均有效流動(dòng)路徑長(zhǎng)度 Li （m）與流體流動(dòng)方向的直線路徑 L （m）的比值，如公式（3）所示.不少學(xué)者[42-47]對(duì)曲折度τ進(jìn)行進(jìn)一步研究后發(fā)現(xiàn)，孔隙率與曲折度的關(guān)系可通過(guò)數(shù)學(xué)表達(dá)式進(jìn)行表述，如表3所示.

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

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

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

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

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

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

值得一提的是，盡管式（7）中的曲折度τ是基于多孔介質(zhì)本身性質(zhì)的參數(shù)，其難于用試驗(yàn)手段直接測(cè)得，但可以通過(guò)孔隙率 n 計(jì)算獲得（見(jiàn)表3）.

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

由擬合模擬的一系列工況所得結(jié)果（圖8中實(shí)線標(biāo)示），可得到孔隙率與曲折度之間的數(shù)學(xué)表達(dá)式為：

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

為驗(yàn)證式（9）的準(zhǔn)確性，本文采用大量已有試驗(yàn)結(jié)果（Chai 等[50]、Jo 等[51-52]、Petrov 和 Rowe[53]、Chen 等[54]、Setz等[55]、Bradshaw 和 Benson[56]以及Seiphoor[57]）進(jìn)行計(jì)算比較，結(jié)果如圖9所示.以 Chai 等[50]數(shù)據(jù)為例，具體的計(jì)算過(guò)程如下.

由式（9）可得此時(shí) GCL 的預(yù)測(cè)滲透系數(shù) k：

本文采用對(duì)數(shù)坐標(biāo)系下滲透系數(shù)偏差分析的表示方法[58-61]（如圖9所示）.圖中，橫坐標(biāo)為預(yù)測(cè)滲透系數(shù)值，縱坐標(biāo)為試驗(yàn)滲透系數(shù)值.從圖中可看出，預(yù)測(cè)值與試驗(yàn)值的范圍為0.23～3.86.按照行業(yè)標(biāo)準(zhǔn)《鈉基膨潤(rùn)土復(fù)合防水襯墊》（FZ/T 64036—2103），滲透系數(shù)值的誤差常在一個(gè)數(shù)量級(jí)內(nèi)進(jìn)行表征，如規(guī)定滲透系數(shù)≤5×10-9～1×10-10 cm/s.同時(shí)，對(duì)數(shù)坐標(biāo)系下滲透系數(shù)的偏差一般介于1/10～10[37，58-61].本文 GCL 的理論預(yù)測(cè)模型計(jì)算結(jié)果與試驗(yàn)值的比在1/5～5，小于一個(gè)數(shù)量級(jí)，說(shuō)明該 GCL 理論預(yù)測(cè)模型可以較為準(zhǔn)確地用于預(yù)測(cè) GCL 的滲透系數(shù)，可以指導(dǎo)工程實(shí)際.

3結(jié)論與展望

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

1）當(dāng)入口邊界條件為均勻壓力時(shí)，在蠕動(dòng)流場(chǎng)中，流體通道半徑的增加會(huì)導(dǎo)致流體速度增加，同時(shí)，流體壓力從入口到出口處呈現(xiàn)均勻遞減的趨勢(shì).

2） GCL試樣中，微觀流動(dòng)粒子分布與宏觀流體場(chǎng)分布規(guī)律基本一致，且大部分流動(dòng)粒子在整體時(shí)間的前1/3段完成運(yùn)移.

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

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

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

參考文獻(xiàn)

[1]? BOUAZZA A. Geosynthetic clay liners[J]. Geotextiles and Geomembranes，2002，20（1）：3-17.

[2]? ROWE R K，MUKUNOKI T，LI M H，et al. Effect of freeze-thawon the permeation of arctic diesel through a GCL[J]. Journal of ASTM International，2004，1（2）：11741.

[3]? SHACKELFORD C D，BENSON C H，KATSUMI T，et al. Evaluating the hydraulic conductivity of GCLs permeated with non- standard liquids[J]. Geotextiles and Geomembranes，2000，18（2/3/4）：133-161.

[4] 劉志彬，施斌，王寶軍.改性膨脹土微觀孔隙定量研究[J].巖土工程學(xué)報(bào)，2004，26（4）：526-530.

LIU Z B，SHI B，WANG B J. Quantitative research on micropores of modified expansive soils[J]. Chinese Journal of Geotechnical Engineering，2004，26（4）：526-530.（In Chinese）

[5] 王寶，董興玲.不同有效應(yīng)力下礦山滲濾液對(duì)土工合成黏土襯墊滲透特性影響的試驗(yàn)研究[J].巖土力學(xué)，2017，38（5）：1350-1358.

WANG B，DONG X L. Hydraulic conductivity of mine leachate through geosynthetic clay liners under different effective stresses [J]. Rock and Soil Mechanics，2017，38（5）：1350-1358.（ In Chinese）

[6]? KANG J B，SHACKELFORD C D. Consolidation of a geosynthetic clay liner under isotropic states of stress [J]. Journal of Geotechnical and Geoenvironmental Engineering，2010，136（1）：253-259.

[7]? MESRI G，OLSON R E. Mechanisms controlling the permeabilityof clays[J]. Clays and Clay Minerals，1971，19（3）：151-158.

[8]? FRATALOCCHI E. Hydraulic perfomance of GCLS with differentstress? history [ C]//Intemational? Workshop “Hydra-physica- mechanies of Landfills”. Grenoble，F(xiàn)rance：University of Joseph Fourier，2005：1-4.

[9]? PETROV R J，ROWE R K，QUIGLEY R M. Selected factors influencing GCL hydraulic conductivity[J]. Journal of Geotechni? cal and Geoenvironmental Engineering，1997，123（8）：683-695.

[10]何俊，施建勇.膨潤(rùn)土中飽和滲透系數(shù)的計(jì)算[J].巖石力學(xué)與工程學(xué)報(bào)，2007，26（ S2）：3920-3925.

HE J，SHI J Y. Calculation of saturated permeability coefficient of bentonite[J]. Chinese Journal of Rock Mechanics and Engi? neering，2007，26（ S2）：3920-3925.（In Chinese）

[11] PITTMAN E D. Relationship of porosity and permeability to various parameters derived from mercury injection capillary pressure curves for sand-stones [J]. AAPG Bulletin ， 1992，76（2）：191-198.

[12] KOZENY J. Uber kapillareleitung des wassersimboden[J]. Sitzungsberichte der Wiener Akademie der Wissenschaften，1927，136（2a）：271-306.

[13] CARMAN P C. Fluid flow through granular beds[J]. ChemicalEngineering Research and Design，1997，75：S32-S48.

[14] NOORUDDIN H A，HOSSAIN M E. Modified Kozeny-Carmencorrelation for enhanced hydraulic flow unit characterization[J]. Journal of Petroleum Science and Engineering，2011，80（1）：107-115.

[15] TSANG Y W. The effect of tortuosity on fluid flow through a singlefracture [J]. Water Resources Research ，1984，20（9）：1209-1215.

[16] MURATA S，MITSUISHI H，SAITO T. Characterization of fracture permeability by using a fractal model[ C]// SPE Asia Pacific Oil and Gas Conference and Exhibition. Melbourne，Australia： SPE，2002：1-8.

[17] DIAS R，TEIXEIRA J A，MOTA M，et al. Tortuosity variation in alow density binary particulate bed[J]. Separation & PurificationTechnology，2006，51（2）：180-184.

[18] XU P，YU B M. Developing a new form of permeability andKozeny-Carman constant for homogeneous porous media by means of fractal geometry[J]. Advances in Water Resources，2008，31（1）：74-81.

[19] CARMAN P C. Flow of gases through porous media[M]. London：Butterworths Scientific Publications，1956：182.

[20] SHAN H Y，LAI Y J. Effect of hydrating liquid on the hydraulicproperties of geosynthetic clay liners [J]. Geotextiles and Geo? membranes，2002，20（1）：19-38.

[21] ROWE R K. Geosynthetic clay liners：Perceptions and misconceptions [J].? Geotextiles and Geomembranes，2020，48（2）：137-156.

[22]郭建春，莊園，劉超.考慮非達(dá)西效應(yīng)的酸蝕裂縫流場(chǎng)數(shù)值模擬[J].巖土力學(xué)，2015，36（11）：3315-3321.

GUO J C，ZHUANG Y，LIU C. Numerical simulation of flow field of acid etched fractures considering non-Darcy effect[J]. Rock and Soil Mechanics，2015，36（11）：3315-3321.（In Chinese）

[23]侯會(huì)明，胡大偉，周輝，等.考慮細(xì)觀等效熱學(xué)參數(shù)的高放廢物處置庫(kù)圍巖應(yīng)力-滲流-溫度耦合模擬方法[J].巖土力學(xué)，2019，40（9）：3625-3634.

HOU H M，HU D W，ZHOU H，et al. Thermo-hydro-mechanical coupling simulation method of surrounding rock in high-level ra? dioactive waste repository considering effective meso-thermal pa ? rameters [J]. Rock and Soil Mechanics，2019，40（9）：3625-3634.（In Chinese）

[24] KOVACS G. Seepage Hydraulics[ M]. Amsterdam：Elsevier Scientific Publishing Company，1981：71.

[25] ZIMMERMAN W B J. COMSOL Multiphysics 有限元法多物理建模與分析[M].北京：人民交通出版社，2007：89-113.

ZIMMERMAN W B J. COMSOL Multiphysics modeling with finite element methods [ M]. Beijing： China Communications Press，2007：89-113.（In Chinese）

[26]張志紅，師玉敏，朱敏.黏土墊層水力–力學(xué)–化學(xué)耦合模型研究[J].巖土工程學(xué)報(bào)，2016，38（7）：1283-1290.

ZHANG Z H，SHI Y M，ZHU M. Coupled hydro-mechanical- chemical model for clay liner[J]. Chinese Journal of Geotechni? cal Engineering，2016，38（7）：1283-1290.（In Chinese）

[27]金磊，曾亞武，程濤，等.基于 LBM-IMB-DEM 分析巖土多孔介質(zhì)的滲透淤堵特性[J].巖土工程學(xué)報(bào)，2021，43（5）：909-917.

JIN L，ZENG Y W，CHENG T，et al. Seepage clogging character? istics of rock and soil porous media using LBM-IMB-DEM simula? tion method[J]. Chinese Journal of Geotechnical Engineering，2021，43（5）：909-917.（In Chinese）

[28]梁越，陳鵬飛，林加定，等.基于透明土技術(shù)的多孔介質(zhì)孔隙流動(dòng)特性研究[J].巖土工程學(xué)報(bào)，2019，41（7）：1361-1366.

LIANG Y，CHEN P F，LIN J D，et al. Pore flow characteristics of porous media based on transparent soil technology[J]. Chinese Journal of Geotechnical Engineering，2019，41（7）：1361-1366.（In Chinese）

[29] AHN H S，JO H Y. Influence of exchangeable cations on hydraulic conductivity of compacted bentonite [J]. Applied Clay Sci? ence，2009，44（1/2）：144-150.

[30]古成中，吳新躍.有限元網(wǎng)格劃分及發(fā)展趨勢(shì)[J].計(jì)算機(jī)科學(xué)與探索，2008，2（3）：248-259.

GU C Z，WU X Y. A review of FEM and trend of development [J]. Journal of Frontiers of Computer Science & Technology，2008，2（3）：248-259.（In Chinese）

[31]李守巨，劉迎曦，孫偉，等.多孔材料孔隙尺寸對(duì)滲透系數(shù)影響的數(shù)值模擬[J].力學(xué)與實(shí)踐，2010，32（4）：12-17.

LI S J，LIU Y X，SUN W，et al. Numerical simulation of relation ? ship between hydraulic conductivity and porosity size for porous geotechnical materials[J]. Mechanics in Engineering，2010，32（4）：12-17.（In Chinese）

[32]周力沛，唐曉武，程冠初，等.雙尺度漸進(jìn)展開(kāi)計(jì)算黏土滲透率影響因素研究[J].巖土工程學(xué)報(bào)，2018，40（7）：1205-1211.

ZHOU L P，TANG X W，CHENG G C，et al. Influencing factors for calculating clay permeability using asymptotic expansion method[J]. Chinese Journal of Geotechnical Engineering，2018，40（7）：1205-1211.（In Chinese）

[33] LAIRD D A. Influence of layer charge on swelling of smectites[J]. Applied Clay Science，2006，34（1/2/3/4）：74-87.

[34] LIU L C. Prediction of swelling pressures of different types ofbentonite in dilute solutions[J]. Colloids and Surfaces A：Physi? cochemical and Engineering Aspects，2013，434：303-318.

[35]項(xiàng)國(guó)圣.含鹽環(huán)境中膨潤(rùn)土的膨脹理論及膨脹衰減機(jī)理研究[ D].上海：上海交通大學(xué)，2007：15-19.

XIANG G S. Theory of swelling propertied and mechanism on swell attenuation of bentonite in solution[ D]. Shanghai：Shanghai Jiaotong University，2007：15-19.（In Chinese）

[36] SARABIAN T，RAYHANI M T. Hydration of geosynthetic clayliners from clay subsoil under simulated field conditions [J]. Waste Management，2013，33（1）：67-73.

[37]王寶，陳彬，竇桐桐，等.水泥固化體淋濾液對(duì) GCL 防滲性能的影響[J].巖土工程學(xué)報(bào)，2019，41（2）：390-396.

WANG B，CHEN B，DOU T T，et al. Influences of stabilization/ solidification product leachates on hydraulic performance of geo ? synthetic clay liners[J]. Chinese Journal of Geotechnical Engi? neering，2019，41（2）：390-396.（In Chinese）

[38] SUMMER M E. Interdisciplinary aspects of soil science. Handbook of Soil Science. Section G[ M]. Boca Raton USA：CRC Press，2000：8.

[39]何俊.基于微觀分析的填埋場(chǎng) GCL 中運(yùn)移規(guī)律的理論研究[ D].南京：河海大學(xué)，2006：45-76.

HE J. Theoretical study of transport through GCL in landfill based on microscopic analysis [ D]. Nanjing：Hohai University，2006：45-76.（In Chinese）

[40] SEIPHOORI A，LALOUI L，F(xiàn)ERRARI A，et al. Water retentionand swelling behaviour of granular bentonites for application in Geosynthetic Clay Liner （ GCL） systems [J]. Soils and Founda? tions，2016，56（3）：449-459.

[41] COLLINS R E. Flow of fluids through porous materials [ M]. Reinhold，New York：Reinhold Publishing Corp，1961：250-256.

[42]孫軍杰，田文通，劉琨，等.基于泊肅葉定律的土體滲透系數(shù)估算模型[J].巖石力學(xué)與工程學(xué)報(bào)，2016，35（1）：150-161.?? SUN J J，TIAN W T，LIU K，et al. Estimation model of soil perme ? ability coefficient based on Poiseuille′s law[J]. Chinese Journal of Rock Mechanics and Engineering，2016，35（1）：150-161.（ In Chinese）

[43] IVERSEN N，J?RGENSEN B B. Diffusion coefficients of sulfateand methane in marine sediments：Influence of porosity[J].Geo? chimica et Cosmochimica Acta，1993，57（3）：571-578.

[44] PISANI L. Simple expression for the tortuosity of porous media[J]. Transport in Porous Media，2011，88（2）：193-203.

[45] KOPONEN A，KATAJA M，TIMONEN J. Tortuous flow in porousmedia[J]. Physical Review E，1996，54（1）：406-410.

[46] ATTIA A M. Effects of petrophysical rock properties on tortuosityfactor[J]. Journal of Petroleum Science and Engineering，2005，48（3/4）：185-198.

[47] KHABBAZI A E，HINEBAUGH J，BAZYLAK A. Determining theimpact of rectangular grain aspect ratio on tortuosity-porosity cor? relations of two-dimensional stochastically generated porous media [J]. Science Bulletin，2016，61（8）：601-611.

[48] YU B M，CHENG P.? A fractal permeability model for bidispersed porous media [J]. International Journal of Heat and Mass Transfer，2002，45（14）：2983-2993.

[49] XIAO B Q，TU X，REN W，et al. Modeling for hydraulic permeability and Kozeny-Carman constant of porous nanofibers using a fractal approach[J]. Fractals，2015，23（3）：1550029.

[50] CHAI J C，PRONGMANEE N. Barrier properties of a geosynthetic clay liner using polymerized sodium bentonite[J].Geotex? tiles and Geomembranes，2020，48（3）：392-399.

[51] JO H Y，KATSUMI T，BENSON C H，et al. Hydraulic conductivity and swelling of nonprehydrated GCLs permeated with single- species salt solutions[J]. Journal of Geotechnical and Geoenvi? ronmental Engineering，2001，127（7）：557-567.

[52] JO H Y，BENSON C H，SHACKELFORD C D，et al. Long-termhydraulic conductivity of a geosynthetic clay liner permeated with inorganic salt solutions[J]. Journal of Geotechnical and Geoenvironmental Engineering，2005，131（4）：405-417.

[53] PETROV R J，ROWE R K. Geosynthetic clay liner （ GCL）-chemical compatibility by hydraulic conductivity testing and fac ? tors impacting its performance[J]. Canadian Geotechnical Jour? nal，1997，34（6）：863-885.

[54] CHEN J，F(xiàn)ANG Y G，GU R G，et al. Study on pore size effect oflow permeability clay seepage [J]. Arabian Journal of Geosci? ences，2019，12（7）：1-10.

[55] SETZ M C，TIAN K，BENSON C H，et al. Effect of ammonium onthe hydraulic conductivity of geosynthetic clay liners[J].Geotex? tiles and Geomembranes，2017，45（6）：665-673.

[56] BRADSHAW S L，BENSON C H. Effect of municipal solid wasteleachate on hydraulic conductivity and exchange complex of geo ? synthetic clay liners[J]. Journal of Geotechnical and Geoenviron ? mental Engineering，2014，140（4）：04013038.

[57] SEIPHOORI A，F(xiàn)ERRARI A，LALOUI L. Water retention behav?iour and microstructural evolution of MX-80 bentonite during wet? ting and drying cycles[J].Géotechnique，2014，64（9）：721-734.

[58] LI Q，CHEN J N，BENSON C H，et al. Hydraulic conductivity ofbentonite-polymer composite geosynthetic clay liners permeated with bauxite liquor[J]. Geotextiles and Geomembranes，2021，49（2）：420-429.

[59] BENSON C H，CHEN J N，EDIL T B，et al. Hydraulic conductiv?ity of compacted soil liners permeated with coal combustion prod ? uct leachates[J]. Journal of Geotechnical and Geoenvironmental Engineering，2018，144（4）：04018011.

[60] CHEN J N，SALIHOGLU H，BENSON C H，et al. Hydraulic con?ductivity of bentonite – polymer composite geosynthetic clay liners permeated with coal combustion product leachates[J]. Journal of Geotechnical and Geoenvironmental Engineering，2019，145（9）：04019038.

[61] SHACKELFORD C D，SEVICK G W，EYKHOLT G R. Hydraulicconductivity of geosynthetic clay liners to tailings impoundment so ? lutions [J].? Geotextiles and Geomembranes，2010，28（2）：149-162.