任豪， 張翔， 雷真

(云南大學(xué)建筑與規(guī)劃學(xué)院， 昆明 650504)

廣泛分布于云南的紅黏土是由碳酸鹽類巖石在高溫潮濕環(huán)境下風(fēng)化形成的，其性狀不同于常規(guī)黏土，是一種具有高含水率、低壓縮性、高液限、低密度、大孔隙比的高塑性黏土，同時(shí)具有優(yōu)良的抗?jié)B性，是較好的筑壩材料。針對(duì)紅黏土在壩體防滲中的應(yīng)用已有較多研究，符鋒等[1]研究了紅黏土作為防滲心墻土在土石壩的應(yīng)用；張燕燕等[2]和許樹春等[3]通過調(diào)節(jié)紅黏土含水率并開展碾壓實(shí)驗(yàn)來研究提高紅黏土防滲心墻施工質(zhì)量的方法。云南省內(nèi)也有較多水利工程采用了這種土料作為防滲材料，如曲靖花山水庫、新平黃草壩水庫等。

云南紅黏土中游離氧化鐵膠結(jié)作用水穩(wěn)定性好，膠結(jié)體在水中不易分散，故抗?jié)B性能較好。常規(guī)的變水頭實(shí)驗(yàn)測(cè)得的滲透系數(shù)一般在1×10-6～1×10-7cm/s。

土的宏觀滲透特性決定于其微觀結(jié)構(gòu)，特別是微觀尺度下的孔隙網(wǎng)絡(luò)形態(tài)。為進(jìn)一步了解土的滲透機(jī)理，近年來很多學(xué)者采用了不同的方法對(duì)該問題開展了研究工作。探究土的微觀結(jié)構(gòu)和微觀滲流的方法主要包括直接掃描建模或數(shù)學(xué)建模以及數(shù)學(xué)分析方法，其中直接掃描法可以獲得真實(shí)的微觀圖像，例如，王超等[4]使用CT(computed tomography)技術(shù)構(gòu)建了黃土巖的孔隙網(wǎng)絡(luò)并分析其不同方向的滲透率；Mukunoki等[5]使用CT掃描構(gòu)建了土體的三維孔隙結(jié)構(gòu)并分析其水力學(xué)性質(zhì)及侵蝕機(jī)理，但受信噪比、分辨率等因素的影響，掃描尺度有限制；陳建平等[6]使用SEM(scanning electron microscope)技術(shù)定量分析了軟土的微觀結(jié)構(gòu)參數(shù)；Ural[7]在研究黏性土的微觀結(jié)構(gòu)時(shí)采用了SEM技術(shù)對(duì)土樣進(jìn)行了三維重構(gòu)，雖然提高了精度，但也相應(yīng)地增加了時(shí)間和成本。除了直接掃描技術(shù)，數(shù)學(xué)建模方法也是重要手段：周力沛等[8]和Sun等[9]用雙尺度方法計(jì)算黏土滲透率，分別選擇圓形和正方形單元作為表征單元體(representative elementary volume， REV)模型，通過探討REV模型的特征來描述土體的滲透特性；胡其志等[10]究黏性土的滲透系數(shù)時(shí)，將土顆粒假定為大小相等、連續(xù)分布的球體，并對(duì)太沙基滲透公式修正后計(jì)算滲透率。這些數(shù)學(xué)建模方法為研究土的微觀結(jié)構(gòu)和滲流提供了新的思路，但在模型構(gòu)建時(shí)都采用了孔隙形態(tài)簡(jiǎn)化，無法體現(xiàn)出孔隙介質(zhì)結(jié)構(gòu)的隨機(jī)性?；谝陨峡紤]，探求一種不依賴于孔隙幾何形態(tài)固定化與簡(jiǎn)單化的孔隙介質(zhì)建模方法勢(shì)在必行。

現(xiàn)提出一種土體微觀結(jié)構(gòu)建模和微觀滲流的數(shù)值分析方法。通過二相化方法將連續(xù)的高斯隨機(jī)場(chǎng)轉(zhuǎn)化為二相場(chǎng)來代表孔隙介質(zhì)的基質(zhì)和孔隙，構(gòu)建出滿足紅黏土真實(shí)孔隙率和孔徑分布的三維孔隙結(jié)構(gòu)模型，該方法能夠擺脫傳統(tǒng)的基于簡(jiǎn)單單元的孔隙介質(zhì)建模。完成建模后，將數(shù)學(xué)形態(tài)學(xué)圖像分析方法計(jì)算機(jī)實(shí)現(xiàn)后用于模型特征分析，得到與滲流有關(guān)的參數(shù)，結(jié)合Kozeny-Carman方程，計(jì)算滲透率，并與實(shí)驗(yàn)數(shù)據(jù)作對(duì)比。通過以上研究流程，可以從微觀孔隙網(wǎng)絡(luò)的角度探討紅黏土的滲透特征。

1 紅黏土孔隙結(jié)構(gòu)與三維建模

1.1 云南紅黏土孔隙結(jié)構(gòu)參數(shù)

研究選取云南曲靖花山水庫紅黏土土樣進(jìn)行孔隙結(jié)構(gòu)參數(shù)研究與孔隙結(jié)構(gòu)建模，如圖1所示。通過壓實(shí)實(shí)驗(yàn)獲得紅黏土壓實(shí)樣，壓實(shí)樣的最優(yōu)含水量約為34%，最大干密度1.35 g/cm3，孔隙比1.03。通過壓汞法測(cè)得壓實(shí)樣孔徑分布于0.01～0.2 μm，如圖2所示。

圖1 云南紅黏土和試樣Fig.1 Yunnan red clay and compacted specimen

圖2 云南紅黏土累積孔徑分布曲線Fig.2 Cumulative pore size distribution curve of Yunnan red clay

根據(jù)孔隙是否連通，可以將孔隙分為開孔隙和閉孔隙，一般來說，閉孔隙在壓實(shí)黏土中占比很少，故可認(rèn)為壓汞法所得孔隙特征近似為真實(shí)孔隙特征[11]。

1.2 紅黏土孔隙結(jié)構(gòu)三維建模

紅黏土三維孔隙結(jié)構(gòu)模型的構(gòu)建采用高斯隨機(jī)場(chǎng)閾值化方法來完成。通過對(duì)空間內(nèi)連續(xù)分布的高斯隨機(jī)場(chǎng)設(shè)定閾值，將連續(xù)場(chǎng)轉(zhuǎn)化為二相場(chǎng)來模擬孔隙介質(zhì)的基質(zhì)和孔隙網(wǎng)絡(luò)。

具體來說：對(duì)于一個(gè)定義在范圍F的一維連續(xù)函數(shù)f(x)，如將函數(shù)的值域設(shè)定一個(gè)閾值n后，函數(shù)可分為兩部分，如圖3所示。將小于閾值的部分定義為一相Qs，其余部分為另一相，即完成一維連續(xù)函數(shù)的二相化處理，所選擇的閾值n決定了二相的比例。此外，函數(shù)越波動(dòng)起伏，二相化后兩相的尺寸越狹窄。

圖3 連續(xù)函數(shù)的二相化Fig.3 Binarization of continuous function

將一維函數(shù)的二相化處理方式擴(kuò)展到三維空間函數(shù)。對(duì)連續(xù)場(chǎng)設(shè)定閾值后轉(zhuǎn)換為二相場(chǎng)，可用于代表孔隙介質(zhì)的二相，如圖4所示，圖4(a)中橘色代表隨機(jī)場(chǎng)三維空間內(nèi)的高取值區(qū)域，藍(lán)色代表低取值區(qū)域，橘色與藍(lán)色相間反映了隨機(jī)場(chǎng)的波動(dòng)情況，圖4(b)中灰色部分代表孔隙對(duì)于高斯隨機(jī)場(chǎng)(期望為0，方差為1)而言，其相關(guān)長(zhǎng)度LC控制著場(chǎng)的形態(tài)起伏[12]隨機(jī)場(chǎng)的相關(guān)長(zhǎng)度越小，具有相關(guān)性的點(diǎn)的分布范圍就越小，場(chǎng)的波動(dòng)也越大，二相化后孔隙的尺寸越??；反之相關(guān)長(zhǎng)度越大，場(chǎng)波動(dòng)越平緩，二相化后孔隙的尺寸也就越大。與此同時(shí)，所設(shè)定的閾值n控制著二相比例，即控制孔隙率的大小，其關(guān)系通過模型分析，可用反誤差函數(shù)表達(dá)，即

圖4 高斯隨機(jī)場(chǎng)二相化Fig.4 Binarization of Gaussian random field

(1)

式(1)中：φ為孔隙率。

通過調(diào)節(jié)這兩個(gè)參數(shù)，可以將高斯隨機(jī)場(chǎng)二相化以后生成不同形態(tài)的孔隙結(jié)構(gòu)模型。

對(duì)紅黏土壓實(shí)樣而言，其孔徑分布在0.01～0.2 μm，最大與最小孔隙尺寸相差20倍。高斯隨機(jī)場(chǎng)的相關(guān)長(zhǎng)度決定了波動(dòng)程度，也決定了二相化之后孔徑的平均大小和有限的分布范圍[1]為了使構(gòu)建的孔隙模型能夠與紅黏土壓實(shí)樣的真實(shí)孔徑分布相吻合，需要將若干個(gè)具有不同相關(guān)長(zhǎng)度的二相場(chǎng)疊加起來以擴(kuò)大模型中的孔徑分布范圍，如圖5所示。同時(shí),每個(gè)場(chǎng)的閾值選擇決定了孔隙率，需參考累積孔徑分布曲線來測(cè)試調(diào)整。

圖5 不同相關(guān)長(zhǎng)度二相場(chǎng)疊加Fig.5 Union of two-phase fields with different LC

本文中建模選取邊長(zhǎng)為0.9 μm的立方體，網(wǎng)格劃分為300×300×300，最小單元尺寸為0.003 μm，此選擇可以呈現(xiàn)紅黏土的全部孔隙部分(0.01～0.2 μm)根據(jù)圖2，設(shè)定8個(gè)不同相關(guān)長(zhǎng)度的高斯隨機(jī)場(chǎng)，設(shè)定孔隙率后計(jì)算閾值，隨機(jī)場(chǎng)閾值化后計(jì)算二相場(chǎng)的孔隙率，結(jié)果如表1所示。

表1 建模相關(guān)參數(shù)Table 1 Modeling related parameters

考慮到多場(chǎng)疊加時(shí)，小孔隙可能會(huì)被大孔隙吞并，因此在根據(jù)孔徑分布曲線設(shè)定各場(chǎng)孔隙率時(shí)，小孔隙二相場(chǎng)的孔隙率會(huì)適當(dāng)增加。將以上八個(gè)二相場(chǎng)聯(lián)合，得到紅黏土壓實(shí)樣三維孔隙結(jié)構(gòu)模型，模型總孔隙率為55.1%，如圖6所示。

圖6 紅黏土孔隙結(jié)構(gòu)模型Fig.6 Pore structure model of red clay

1.3 孔徑分布驗(yàn)證

在模型孔隙率基本滿足(壓實(shí)樣孔隙比1.03)的前提下，借助數(shù)學(xué)形態(tài)學(xué)方法對(duì)模型孔徑分布進(jìn)行驗(yàn)證。數(shù)學(xué)形態(tài)學(xué)是一種基于拓?fù)鋵W(xué)、格論、集合論的幾何結(jié)構(gòu)分析及處理的理論，用于處理二值圖像。圖像的一些幾何和拓?fù)涮卣?，如尺寸、測(cè)地線距離、形狀、連通性等，可通過形態(tài)學(xué)方法進(jìn)行分析[13-15]，需用到數(shù)學(xué)形態(tài)學(xué)開運(yùn)算，即設(shè)定一個(gè)結(jié)構(gòu)元素，將孔隙網(wǎng)絡(luò)和結(jié)構(gòu)元素進(jìn)行比較，小于結(jié)構(gòu)元素尺寸的孔隙部分將被剔除，統(tǒng)計(jì)剩余的孔隙率[16]。不斷增大結(jié)構(gòu)元素的尺寸，即可對(duì)模型孔隙網(wǎng)絡(luò)尺寸進(jìn)行篩分，如圖7所示。如對(duì)紅黏土模型進(jìn)行結(jié)構(gòu)元素為21 nm的開運(yùn)算，剩余孔隙率為44.0%；進(jìn)行結(jié)構(gòu)元素為51 nm的開運(yùn)算，剩余孔隙率為17.6%，開運(yùn)算結(jié)果如圖8所示。

圖7 形態(tài)學(xué)開運(yùn)算示意圖Fig.7 Schematic diagram of morphological opening

圖8 紅黏土孔隙模型開運(yùn)算Fig.8 Morphological opening on red clay pore model

模型孔徑分布與紅黏土壓實(shí)樣實(shí)際孔徑分布對(duì)比如圖9所示，規(guī)律基本吻合。通過上述方法，在保證孔隙特征參數(shù)(孔隙率和孔徑分布)的前提下，呈現(xiàn)出了隨機(jī)孔隙網(wǎng)絡(luò)形態(tài)的紅黏土模型。

圖9 累積孔徑分布曲線對(duì)比Fig.9 Comparison of cumulative pore size distribution curves

2 云南紅黏土滲透性能與數(shù)值模擬

2.1 云南紅黏土滲透性

通過變水頭滲透實(shí)驗(yàn)測(cè)定云南紅黏土擊實(shí)樣滲透率，測(cè)得結(jié)果如表2所示。

表2 云南紅黏土滲透實(shí)驗(yàn)結(jié)果Table 2 Permeability test results of Yunnan red clay

2.2 紅黏土滲透數(shù)值模擬

根據(jù)Kozeny-Carman方程，孔隙介質(zhì)的滲透率估算方法為

(2)

式(2)中:K為滲透率；φ為孔隙率；C為Kozeny-Carman常數(shù)；S為連通孔隙比表面積。其中C可以通過路徑迂曲度τ來計(jì)算[17],即

C=2τ2

(3)

基于2.1節(jié)所構(gòu)建的紅黏土壓實(shí)樣三維孔隙結(jié)構(gòu)模型，可通過數(shù)學(xué)形態(tài)學(xué)方法對(duì)圖像進(jìn)行分析，計(jì)算模型連通孔隙的迂曲度τ和比表面積S，進(jìn)而采用Kozeny-Carman方程對(duì)滲透系數(shù)進(jìn)行數(shù)值模擬。具體采用的研究方法如下：

2.2.1 迂曲度計(jì)算

迂曲度τ反映流體滲流路徑的曲折程度，表達(dá)式為

(4)

如圖10所示，孔隙網(wǎng)絡(luò)形態(tài)復(fù)雜，滲流發(fā)生在連通孔隙中。在數(shù)值分析時(shí)，需剔除模型中的非連通孔隙，該步驟通過數(shù)學(xué)形態(tài)學(xué)中的測(cè)地重建圖像處理方法來實(shí)現(xiàn)[17-18]。

圖10 迂曲度示意圖Fig.10 Diagram of tortuosity

如圖11所示，方法具體描述為：將孔隙模型的整個(gè)孔隙網(wǎng)絡(luò)圖像設(shè)為一個(gè)集合，稱為mask(紅色部分)，將與滲流方向的上游邊界面接觸的孔隙部分設(shè)為另一個(gè)集合，稱為marker(藍(lán)色部分)。完成初始設(shè)定后，沿著marker的邊界對(duì)marker進(jìn)行等尺寸膨脹，尺寸決定于所選擇的結(jié)構(gòu)元素的大小。膨脹后與mask取交集，則marker會(huì)沿著開口孔隙的路徑往模型內(nèi)部延伸一段距離(取決于膨脹的大小)。隨后重復(fù)“膨脹+相交”步驟直至交集后沒有新增孔隙，則與上游邊界面相連的開孔隙可被提取出來。

圖11 測(cè)地重建示意圖Fig.11 Schematic diagram of geodesic reconstruction

同理，將模型下游邊界面上的孔設(shè)為marker,重復(fù)測(cè)地重建，可得到與下游面相連的孔隙。將上、下游測(cè)地重建的結(jié)果取交集后可獲得模型上下游面間的連通孔隙。

提取的紅黏土模型滲流孔隙網(wǎng)絡(luò)如圖12所示(流向：深色端流向淺色端)。根據(jù)模型邊長(zhǎng)以及“膨脹+交集”迭代運(yùn)算的次數(shù)可獲得滲流路徑長(zhǎng)度和迂曲度的取值范圍，結(jié)果如表3所示。

圖12 紅黏土模型滲流孔隙網(wǎng)絡(luò)Fig.12 Seepage pore network of red clay model

表3 紅黏土模型滲流路徑分析結(jié)果Table 3 Analysis results of seepage paths of red clay model

2.2.2 比表面積與滲透系數(shù)的計(jì)算

比表面積定義為單位體積模型內(nèi)孔隙的總表面積，由于閉孔隙并不參與滲流過程，在計(jì)算表面積時(shí)只考慮連通孔隙表面積。

由于孔隙模型構(gòu)建于0.9 μm邊長(zhǎng)、300×300×300網(wǎng)格劃分的立方體內(nèi)，將模型圖像放大，可以發(fā)現(xiàn)不同形態(tài)的孔隙表面均呈階梯狀，如圖13(a)所示。如將基質(zhì)的單元格表示為0，孔隙的單元格表示為1，如圖13(b)所示,通過統(tǒng)計(jì)兩種類型單元格交界面的個(gè)數(shù)，再乘以每個(gè)交界面的面積(即單元格一個(gè)表明的面積)，可計(jì)算出孔隙的表面積。將此思路程序?qū)崿F(xiàn)后應(yīng)用于3.2.1所得的連通孔隙模型，得到總表面積為2.3×107nm2，除以模型體積后，比表面積為0.03 nm2/nm3。

圖13 計(jì)算孔隙表面積方法示意圖Fig.13 Schematic diagram of pore surface area calculation method

將迂曲度代入式(3)，求出Kozeny-Carman常數(shù)。將各個(gè)參數(shù)代入Kozeny-Carman方程[式(2)]，可計(jì)算模型的滲透率K, 再計(jì)算滲透系數(shù)k，計(jì)算公式為

(5)

式(5)中：ρ為水的密度；g為重力加速度；η為水在20 ℃時(shí)的動(dòng)力黏滯性系數(shù)(1.01×10-6kPa·s)，計(jì)算結(jié)果如表4所示。

表4 計(jì)算模型滲透系數(shù)Table 4 Calculation of permeability coefficient of model

對(duì)比紅黏土壓實(shí)樣滲透系數(shù)計(jì)算值(0.24×10-7～0.54×10-7cm/s)與實(shí)驗(yàn)值(1.32×10-7～3.05×10-7cm/s)結(jié)果相似程度較高，驗(yàn)證了模型構(gòu)建和滲流分析方法的可行性，但需要說明和注意的是：

(1)模型網(wǎng)格劃分越密，單元體越小，一方面在建模時(shí)候能夠涵蓋的孔徑跨度更大，建模越準(zhǔn)確；另一方面，模型中階梯狀的孔隙表面越能接近光滑曲面，表面積會(huì)減小，相應(yīng)的計(jì)算滲透系數(shù)會(huì)越大。增加網(wǎng)格密度會(huì)導(dǎo)致相應(yīng)的建模和計(jì)算時(shí)間也越長(zhǎng)，對(duì)計(jì)算機(jī)計(jì)算性能的要求也越高。

(2)紅黏土中含有大量黏粒，表面帶負(fù)電，水在滲流過程中會(huì)受黏粒表面電性吸引，形成吸附水層。當(dāng)水層厚度與孔隙尺寸相當(dāng)時(shí)，會(huì)構(gòu)成結(jié)合水阻礙滲流。因此，紅黏土的滲流可能并非發(fā)生在整個(gè)連通的孔隙網(wǎng)絡(luò)中，某些較小尺寸的連通孔隙可能并非滲流的基座。而在模型計(jì)算時(shí)考慮的是整個(gè)連通孔隙網(wǎng)絡(luò)，并未對(duì)該部分孔隙專門分析，也可能導(dǎo)致計(jì)算滲透系數(shù)與真實(shí)情況有偏差。后續(xù)需針對(duì)該問題進(jìn)一步研究。

3 結(jié)論

以云南紅黏土作為筑壩材料的滲透性能為研究背景，提出了一種針對(duì)孔隙介質(zhì)三維孔隙結(jié)構(gòu)建模和滲透系數(shù)數(shù)值分析的方法，具體結(jié)論如下。

(1)模型構(gòu)建基于紅黏土壓實(shí)樣的孔隙結(jié)構(gòu)參數(shù)。通過設(shè)定閾值將連續(xù)分布的高斯隨機(jī)場(chǎng)轉(zhuǎn)化為二相場(chǎng)來表征孔隙介質(zhì)的基質(zhì)和孔隙。將多個(gè)具有不同相關(guān)長(zhǎng)度和閾值的二相場(chǎng)疊加起來，構(gòu)建出滿足紅黏土真實(shí)孔隙率和孔徑分布的三維孔隙結(jié)構(gòu)模型。該方法不同于常見的對(duì)孔隙網(wǎng)絡(luò)進(jìn)行形態(tài)簡(jiǎn)化的建模方法，能夠提供隨機(jī)的孔隙結(jié)構(gòu)。其中，模型孔徑分布的驗(yàn)算通過數(shù)學(xué)形態(tài)學(xué)圖像處理方法中的開運(yùn)算算法完成。

(2)基于所建模型，采用數(shù)學(xué)形態(tài)學(xué)中的測(cè)地重建和開運(yùn)算方法，計(jì)算了滲流路徑的相關(guān)參數(shù)，如路徑長(zhǎng)度、迂曲度、體積分?jǐn)?shù)、表面積等，并通過Kozeny-Carman方程估算了模型的滲透系數(shù)。通過與實(shí)驗(yàn)數(shù)據(jù)對(duì)比，驗(yàn)證了方法的可行性。

(3)數(shù)值分析結(jié)果表明：紅黏土的微觀孔隙網(wǎng)絡(luò)結(jié)構(gòu)決定了其低滲透率，決定了紅黏土作為筑壩材料的較好的抗?jié)B性能。

(4)文中提及的方法可為孔隙介質(zhì)的建模和滲流分析提供新的思路。