不同壓實(shí)度花崗巖殘積土的滲流模擬

王志兵，談勛勛，王 遠(yuǎn)，孫 廣，劉金明

(1.桂林理工大學(xué)土木與建筑工程學(xué)院，桂林 541004；2.廣西巖土力學(xué)與工程重點(diǎn)實(shí)驗(yàn)室，桂林 541004)

花崗巖殘積土作為一種路基填料，當(dāng)其壓實(shí)度不足時(shí)，會(huì)導(dǎo)致路基災(zāi)害或病害的發(fā)生。如2015年12月20日深圳光明新區(qū)，廢土填埋場(chǎng)由于風(fēng)化花崗巖壓實(shí)度不足導(dǎo)致滑坡[1]。以花崗巖殘積土作為路堤邊坡的廣佛肇高速公路多處出現(xiàn)淺層滑移、水毀、崩塌[2]。眾多的路基病害與填料的壓實(shí)度有緊密的關(guān)聯(lián)。因此，研究不同壓實(shí)度對(duì)花崗巖殘積土滲透特性的影響不僅對(duì)實(shí)際具有應(yīng)用價(jià)值，而且對(duì)巖土領(lǐng)域有理論意義。以往研究人員對(duì)土體的認(rèn)識(shí)更多局限在宏觀尺度上[3]，未從土體內(nèi)部的孔隙分布、孔喉分布等微觀結(jié)構(gòu)特征考慮，導(dǎo)致無(wú)法研究滲流的微觀影響因素[4]。

為了研究壓實(shí)度對(duì)多孔介質(zhì)微觀滲流特性的影響，第一步是獲得其孔隙空間結(jié)構(gòu)分布特征。目前，獲得多孔介質(zhì)孔隙結(jié)構(gòu)特征主流的研究方法有：壓汞法(mercury intrusion porosimetry,MIP)[5]、氮?dú)馕椒?nitrogen adsorption method,NAI)[6]、光學(xué)顯微鏡或掃描電鏡法[7]、CT掃描圖像法[8]。壓汞法會(huì)產(chǎn)生“墨水瓶”效應(yīng)且破壞土體的結(jié)構(gòu)。氮?dú)馕椒ㄖ荒軠y(cè)量微小的孔隙。光學(xué)顯微鏡或掃描電鏡法只能定性表征土體的局部孔隙結(jié)構(gòu)。相比與其他方法，CT掃描圖像法可以快速無(wú)損的獲取樣品孔隙空間結(jié)構(gòu)參數(shù)，定量表征樣品的孔隙空間結(jié)構(gòu)，建立三維可視化模型[9-12]。

以往表征流體流動(dòng)規(guī)律是建立在以達(dá)西定律為基礎(chǔ)的宏觀滲流方程上，由于CT掃描圖像法的出現(xiàn)，通過(guò)CT掃描圖像法提取孔隙網(wǎng)絡(luò)模型，在獲取樣品孔隙空間結(jié)構(gòu)的同時(shí)研究微米級(jí)流體輸運(yùn)特征[13-25]。由于土體孔隙結(jié)構(gòu)復(fù)雜和計(jì)算機(jī)性能的限制，導(dǎo)致直接運(yùn)用孔隙網(wǎng)絡(luò)模型進(jìn)行滲流模擬較為困難。中外學(xué)者多以流體力學(xué)的傳統(tǒng)方法CFD及格子玻爾茲曼方法(lattice boltzmann method,LBM)模擬數(shù)字巖心模型的流體輸運(yùn)特征，建立微米級(jí)滲流模擬，促進(jìn)微觀滲流機(jī)理的研究。

現(xiàn)以桂東南地區(qū)的花崗巖殘積土為例，通過(guò)Micro-CT掃描與三維可視化軟件Avizo結(jié)合，建立孔隙網(wǎng)絡(luò)模型。目前，此項(xiàng)對(duì)接技術(shù)主要用于碳酸鹽巖、砂巖、頁(yè)巖、煤、天然氣水合物及多孔碳的研究，對(duì)于流體在花崗巖殘積土的研究極少。運(yùn)用Avizo相關(guān)算法，提取花崗巖殘積土的表單元體及孔隙空間結(jié)構(gòu)，并將表單元體劃分網(wǎng)格后導(dǎo)入多場(chǎng)耦合軟件Comsol進(jìn)行花崗巖殘積土的滲流模擬，通過(guò)計(jì)算單相流體流經(jīng)該數(shù)字巖心模型的絕對(duì)滲透率，與實(shí)測(cè)值對(duì)比，檢驗(yàn)該模型準(zhǔn)確性、代表性，并研究壓實(shí)度對(duì)花崗巖殘積土滲流特性的影響。

1 基于高分辨率CT三維重建

1.1 土樣壓實(shí)度選取及制備

由室內(nèi)擊實(shí)試驗(yàn)得到桂東南地區(qū)花崗巖殘積土的最佳含水率18.9%，最大干密度為1.65 g/cm3[26]，則壓實(shí)度90%、100%、110%分別對(duì)應(yīng)的干密度為1.485、1.65、1.815 g/cm3；試樣按三個(gè)干密度值在含水率為18.9%時(shí)按壓樣法制樣。然后制作環(huán)刀樣和CT樣分別進(jìn)行滲透試驗(yàn)和CT掃描試驗(yàn)。CT樣為15 mm×15 mm×12 mm的立方體，并用聚乙烯膠盒包裝。

1.2 CT掃描

Micro-CT原理是根據(jù)花崗巖殘積土中不同密度的成分對(duì)X射線吸收系數(shù)不同，以達(dá)到區(qū)分孔隙和基質(zhì)的目的。試驗(yàn)依托三維X射線顯微鏡Xradia 510 Versa高分辨率CT系統(tǒng)，對(duì)桂東南地區(qū)不同壓實(shí)度的花崗巖殘積土進(jìn)行CT掃描，為了減小對(duì)巖心的損壞，從15 mm×15 mm×12 mm樣品中選取1 010×1 000×1 024體素的部分進(jìn)行掃描。高分辨率Micro-CT系統(tǒng)及樣品如圖1所示。得到一組二維切片(1 010張，1 000×1 024像素)，每張切片的灰度值為216，壓實(shí)度為90%、100%、110%的花崗巖殘積土的體素分別為5.81、5.81、6.2 μm。在樣品的切片中，黑色表示花崗巖殘積土孔隙，灰色表示骨架基質(zhì)。樣品的切片如圖2所示。

圖1 Micro-CT儀器及土樣Fig.1 Micro-CT and samples

1.3 圖像處理

1.3.1 濾波處理

CT掃描獲得的切片過(guò)程難免會(huì)產(chǎn)生各種系統(tǒng)噪聲，不僅降低圖像質(zhì)量，還直接影響后續(xù)圖像處理和孔隙結(jié)構(gòu)的定量分析。濾波處理可以剔除信號(hào)中特定頻率的波段，來(lái)提高圖像質(zhì)量。朱洪林[9]、呂邦民[21]使用中值濾波分別對(duì)低滲砂巖、多孔碳的灰度圖像進(jìn)行處理，何凱凱使用非局部均值濾波對(duì)煤的灰度圖像進(jìn)行處理。針對(duì)均值、中值、高斯濾波算法的優(yōu)缺點(diǎn)，且處理效果側(cè)重點(diǎn)不同，選用單一濾波算法來(lái)處理灰度圖像通常達(dá)不到預(yù)期的效果，花崗巖殘積土經(jīng)過(guò)中值濾波和非局部均值濾波的處理效果如圖2所示，根據(jù)處理效果選用非局部均值濾波(non-local mean)對(duì)灰度圖像進(jìn)行處理，該濾波結(jié)合高斯濾波和均值濾波的特點(diǎn)。花崗巖殘積土的灰度圖像經(jīng)過(guò)非局部均值濾波處理之后，孔隙和花崗巖殘積土骨架基質(zhì)的邊界更容易區(qū)分。

1.3.2 閾值分割

閾值分割是重建后的模型能否準(zhǔn)確描述實(shí)際物理模型的基礎(chǔ)，其本質(zhì)是通過(guò)圖像的灰度直方圖獲得分割的閾值。利用孔隙和基質(zhì)之間灰度值的差異來(lái)確定兩者的灰度閾值，然后對(duì)花崗巖殘積土進(jìn)行三維重建，從而提取出孔隙的相關(guān)信息。受CT掃描分辨率及花崗巖殘積土中組分復(fù)雜的限制，導(dǎo)致孔隙和基質(zhì)閾值選取更加復(fù)雜。閾值選取非常關(guān)鍵，如果選取一個(gè)合適的閾值，就能將孔隙準(zhǔn)確的分割開(kāi)來(lái)。呂邦民[21]使用OTSU(最大類間方差法 )算法對(duì)多孔碳的孔隙與基質(zhì)進(jìn)行分割。針對(duì)本次試驗(yàn)得到的CT圖像，并結(jié)合花崗巖殘積土孔隙結(jié)空間結(jié)構(gòu)特征，進(jìn)行多種閾值分割方法比選，最終選擇OTSU算法。

圖2 原始切片及濾波處理Fig.2 Original slice and filter processing

1.3.3 表單元體提取

解決計(jì)算機(jī)內(nèi)存容量、運(yùn)行速度與數(shù)字巖心尺寸大小間矛盾是通過(guò)選取合適的代表性體積單元(representative elementary volume，REV)。小于REV尺度取得的土體物理特性波動(dòng)明顯，而大于REV尺度取得的土體物理特性趨于穩(wěn)定。因此，選擇大小合適的REV不僅能代表土體的物理特性，而且也能滿足計(jì)算機(jī)的條件。當(dāng)REV尺寸與孔隙度的變化規(guī)律趨于穩(wěn)定時(shí)，REV尺寸即為對(duì)應(yīng)的最小單元體的邊長(zhǎng)。通過(guò)分析REV尺寸與孔隙度的變化規(guī)律，發(fā)現(xiàn)壓實(shí)度為90%、100%、110%花崗巖殘積土的REV尺寸分別大于220×220×220、200×200×200、200×200×200體素時(shí)，土樣的孔隙度基本不隨單元體尺寸變化，如圖3所示。因此，壓實(shí)度為90%、100%、110%花崗巖殘積土選取表單元體尺寸分別為220×220×220、200×200×200、200×200×200體素，則對(duì)應(yīng)的表單元體大小分別為1 278 μm×1 278 μm×1 278 μm、1 162 μm×1 162 μm×1 162 μm、1 240 μm×1 240 μm×1 240 μm。

2 微觀孔隙結(jié)構(gòu)模型構(gòu)建

三維可視化是將花崗巖殘積土孔隙與基質(zhì)的分布結(jié)構(gòu)以更加直觀的方式呈現(xiàn)。構(gòu)建的孔隙網(wǎng)絡(luò)模型較為真實(shí)地還原花崗巖殘積土的孔隙分布及連通性。其過(guò)程是將分割后的圖像進(jìn)行三維重構(gòu)?；谶x定的閾值和表單元體的尺寸，使用直體積渲染模塊，完成花崗巖殘積土三維重建，如圖4所示。

圖3 REV尺寸與孔隙度的關(guān)系Fig.3 The relationship between REV size and porosity

圖4 花崗巖殘積土的三維重構(gòu)Fig.4 Three-dimensional reconstruction of granite residual soil

3 滲流模擬

3.1 網(wǎng)格剖分、優(yōu)化

數(shù)字巖心模型擁有高質(zhì)量的四面體網(wǎng)格是利用有限元法進(jìn)行滲流模擬的重要條件。由于花崗巖殘積土孔隙空間結(jié)構(gòu)復(fù)雜，直接對(duì)重建后的數(shù)字巖心模型進(jìn)行網(wǎng)格劃分極其困難，而且導(dǎo)入的模型無(wú)法進(jìn)行滲流模擬。因此，采用Mimics軟件，多次修補(bǔ)破面和優(yōu)化數(shù)字巖心曲面模型，其中包括修復(fù)尖角、缺口，消除小孔、重合邊、交叉、縫隙、倒角等，然后生成高質(zhì)量的網(wǎng)格。最終將STL文件導(dǎo)入Comsol軟件進(jìn)行滲流模擬。

3.2 控制方程與邊界條件

將實(shí)際孔隙結(jié)構(gòu)的數(shù)字巖心模型作為模擬基礎(chǔ)，并展開(kāi)滲流模擬。其方法有兩類：計(jì)算流體力學(xué)的傳統(tǒng)方法(CFD)及格子玻爾茲曼Boltzman方法(LBM)。CFD是對(duì)介質(zhì)模型進(jìn)行質(zhì)量、動(dòng)量、能量守恒方程求解，并對(duì)控制方程及模型離散。LBM是一種基于分子運(yùn)動(dòng)理論模擬流體流動(dòng)的方法，流體被看成有質(zhì)量沒(méi)體積的微粒，通過(guò)模擬微粒在網(wǎng)格中的碰撞和移動(dòng)，構(gòu)建簡(jiǎn)化的運(yùn)動(dòng)論模型來(lái)反映孔隙空間中流體滲流規(guī)律。但LBM方法計(jì)算量大、耗時(shí)長(zhǎng)。因此采用CFD方法中不可壓縮Navier-Stokes方程模塊來(lái)進(jìn)行孔隙空間結(jié)構(gòu)的流體流動(dòng)模擬。

(1)

(2)

(3)

式中：ρ為流體密度，kg/m3；t為時(shí)間，s；u為流速，m/s；Δ為拉普拉斯算子；e為內(nèi)能，J；σ為應(yīng)力張量；q為熱通量。

流體屬性按常溫下水的參數(shù)賦值，對(duì)花崗巖殘積土的模型分別進(jìn)行X、Y、Z三個(gè)方向滲流模擬，設(shè)置兩個(gè)相對(duì)立面作為入口與出口邊界，4個(gè)側(cè)面設(shè)為自由滑移壁面，其余壁面視為無(wú)滑移壁面。

為了簡(jiǎn)化模擬與分析過(guò)程，作如下假設(shè)：①水只在土體孔隙內(nèi)流動(dòng)，不會(huì)滲透到花崗巖殘積土骨架基質(zhì)；②水為不可壓縮牛頓流體，且在孔隙流動(dòng)的過(guò)程中溫度保持不變；③水是連續(xù)流動(dòng)的；④水只受到重力和壓力的影響。

3.3 絕對(duì)滲透率模擬

在Comsol軟件數(shù)值模擬結(jié)果中，對(duì)建立好的模型進(jìn)行計(jì)算，對(duì)花崗巖殘積土在X、Y、Z方向的出口或入口邊界水的流動(dòng)速度進(jìn)行積分，得到體積流量，再根據(jù)達(dá)西定律獲得其絕對(duì)滲透率：

(4)

式(4)中：Q為體積流量，m3/s；A為土體滲流截面積，m2；L為土體滲流長(zhǎng)度，m；μ為流體黏度，Pa·s；ΔP為壓差，Pa；K為絕對(duì)滲透率，m2。

由表1、圖5～圖7可知，壓實(shí)度為90%、100%花崗巖殘積土在X、Y、Z方向的滲透率差異較小，但壓實(shí)度為110%花崗巖殘積土在3個(gè)方向的滲透率差異較大，X方向的滲透率是Z方向的9.27倍，其原因可能是當(dāng)壓實(shí)度增加至110%時(shí)，分層壓實(shí)對(duì)Z方向上的滲透率影響較大；當(dāng)壓實(shí)度增加至110%，花崗巖殘積土的滲透率改變2個(gè)數(shù)量級(jí)，其原因可能是當(dāng)壓實(shí)度增加至110%時(shí)，花崗巖殘積土的孔隙形態(tài)呈扁平狀，孔隙直徑減小，土顆粒與其周?chē)芤旱奈锢砘瘜W(xué)作用會(huì)產(chǎn)生結(jié)合水膜，其對(duì)土顆粒會(huì)產(chǎn)生吸附作用，由于距離不同，故吸附力的強(qiáng)弱會(huì)有所不同。距離越近，結(jié)合水膜越穩(wěn)定，在滲透中基本不發(fā)生運(yùn)動(dòng)。因此，結(jié)合水膜的存在會(huì)削弱土體的滲透性且從壓力云圖可以看出壓實(shí)度為110%花崗巖殘積土相對(duì)于壓實(shí)度為90%、100%花崗巖殘積土在3個(gè)方向的孔隙空間結(jié)構(gòu)發(fā)生較大的變化，分布不均勻，進(jìn)出口面積相差較大。壓實(shí)度為90%、100%、110%花崗巖殘積土在X、Y、Z方向上滲透率模擬值分別是實(shí)測(cè)滲透率0.89、0.53、1.29、1.74、1.71、0.62、2.02、0.88、0.22倍；其中壓實(shí)度為90%、100%、110%花崗巖殘積土滲透率模擬值平均值分別是實(shí)測(cè)滲透率0.9、1.35、0.77倍。表明孔隙結(jié)構(gòu)的數(shù)字巖心模型具有準(zhǔn)確性、代表性。存在差異主要是因?yàn)榛◢弾r殘積土在不同方向上存在一定的非均質(zhì)性且Micro-CT的分辨率有限，無(wú)法識(shí)別小于分辨率的孔隙，從而會(huì)影響滲透率的計(jì)算。

表1 花崗巖殘積土滲透率模擬結(jié)果Table 1 Permeability simulation results of granite residual soil

3.4 壓力場(chǎng)分布

通過(guò)Comsol后處理得到滲流模擬的速度場(chǎng)分布。雖然X、Y、Z方向的壓強(qiáng)差相等，但由于孔隙直徑、曲折度、連通性及形狀不同，導(dǎo)致3個(gè)方向上的壓力場(chǎng)分布存在一定的差異，如圖5～圖7所示。具體表現(xiàn)如下：壓力隨著水的流動(dòng)方向逐漸減小,最大壓力出現(xiàn)在入口附近狹窄孔隙處；當(dāng)孔隙通道突然變小時(shí)，壓力變化快，壓力表現(xiàn)為先降低后升高，主要是孔隙通道急劇減小和變彎所致。

圖5 壓實(shí)度90%花崗巖殘積土的壓力云圖Fig.5 Pressure cloud map of 90% compaction granite residual soil

圖6 壓實(shí)度100%花崗巖殘積土的壓力云圖Fig.6 Pressure cloud map of 100% compaction granite residual soil

圖7 壓實(shí)度110%花崗巖殘積土的壓力云圖Fig.7 Pressure cloud map of 110% compaction granite residual soil

3.5 速度場(chǎng)分布

由于連通性及壓力的分布導(dǎo)致X、Y、Z方向上的速度場(chǎng)分布存在一定的差異，隨著壓實(shí)度的增加，花崗巖殘積土的連通性變差，如圖8～圖10所示。雖然所有孔隙均已連通，但是有部分孔隙中沒(méi)有水流過(guò)，主要是這部分孔隙小、距入口較遠(yuǎn)且滲流通道較長(zhǎng)。平均速度隨著壓力的降低而增加，當(dāng)孔隙通道突然變小時(shí)，水的流速突然升高；選取X、Y、Z方向截面的速度圖，發(fā)現(xiàn)對(duì)于單個(gè)孔隙而言，沿著孔隙中心向外壁的方向，水的滲流速度逐漸減小。對(duì)整個(gè)連通孔隙而言，水的最大流動(dòng)速度在連通性好的中心區(qū)域。

如表2、表3所示，在600、800 μm截面上，隨著壓實(shí)度的增加，花崗巖殘積土的平均滲流速度和最大流速逐漸減小，僅壓實(shí)度為100%花崗巖殘積土在Y軸800 μm截面的流速反常；隨著壓實(shí)度的增加，在u≥2 m/s、1 m/s≤u＜2 m/s、0.5 m/s＜u＜1 m/s區(qū)間花崗巖殘積土的流速占比大致呈現(xiàn)逐漸減小，在≤0.5 m/s這個(gè)區(qū)間花崗巖殘積土的流速占比逐漸增大。

圖8 壓實(shí)度90%各截面流體流動(dòng)速度分布Fig.8 Fluid flow velocity distribution diagram of each section at 90% compaction degree

圖9 壓實(shí)度100%各截面流體流動(dòng)速度分布Fig.9 Fluid flow velocity distribution diagram of each section at 100% compaction degree

圖10 壓實(shí)度110%各截面流體流動(dòng)速度分布Fig.10 Fluid flow velocity distribution diagram of each section at 110% compaction degree

4 結(jié)論

以桂東南地區(qū)不同壓實(shí)度的花崗巖殘積土為研究對(duì)象，基于CT掃描試驗(yàn)與數(shù)字圖像處理技術(shù)，構(gòu)建其三維孔隙網(wǎng)絡(luò)模型，并基于CFD研究水在土體孔隙中的滲流規(guī)律，得到以下結(jié)論。

(1)CFD作為一種微觀滲流模擬方法，不僅可以對(duì)不同壓實(shí)度的花崗巖殘積土數(shù)字巖心模型的滲流過(guò)程進(jìn)行定量表征，還能夠直觀描述流體在孔隙內(nèi)流動(dòng)的動(dòng)態(tài)過(guò)程。且獲得花崗巖殘積土的絕對(duì)滲透率，與實(shí)測(cè)滲透率進(jìn)行驗(yàn)證，能為花崗巖殘積土滲流模擬提供可靠的模型。此次研究拓寬了CT掃描技術(shù)在巖土工程領(lǐng)域的應(yīng)用。

表3 各方向上截面各速度區(qū)間的占比及最大值Table 3 The proportion and maximum value of each velocity section in the cross section in all directions

(2)花崗巖殘積土的滲流過(guò)程中，在連通性好的孔隙中心區(qū)域流速最大，流速沿著孔隙中心向外壁的方向，流體的滲速逐漸減小；在連通性差的孔隙區(qū)域流速接近為0。

(3)隨著壓實(shí)程度的增加，在X、Y、Z方向上其截面的孔隙發(fā)生較大改變，花崗巖殘積土的連通性變差，且截面的平均滲流速度減小，這表明壓實(shí)度對(duì)花崗巖殘積土的滲流特性有重要影響。