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

王志兵，談勛勛，王 遠，孫 廣，劉金明

(1.桂林理工大學土木與建筑工程學院，桂林 541004；2.廣西巖土力學與工程重點實驗室，桂林 541004)

花崗巖殘積土作為一種路基填料，當其壓實度不足時，會導致路基災害或病害的發生。如2015年12月20日深圳光明新區，廢土填埋場由于風化花崗巖壓實度不足導致滑坡[1]。以花崗巖殘積土作為路堤邊坡的廣佛肇高速公路多處出現淺層滑移、水毀、崩塌[2]。眾多的路基病害與填料的壓實度有緊密的關聯。因此，研究不同壓實度對花崗巖殘積土滲透特性的影響不僅對實際具有應用價值，而且對巖土領域有理論意義。以往研究人員對土體的認識更多局限在宏觀尺度上[3]，未從土體內部的孔隙分布、孔喉分布等微觀結構特征考慮，導致無法研究滲流的微觀影響因素[4]。

為了研究壓實度對多孔介質微觀滲流特性的影響，第一步是獲得其孔隙空間結構分布特征。目前，獲得多孔介質孔隙結構特征主流的研究方法有：壓汞法(mercury intrusion porosimetry,MIP)[5]、氮氣吸附法(nitrogen adsorption method,NAI)[6]、光學顯微鏡或掃描電鏡法[7]、CT掃描圖像法[8]。壓汞法會產生“墨水瓶”效應且破壞土體的結構。氮氣吸附法只能測量微小的孔隙。光學顯微鏡或掃描電鏡法只能定性表征土體的局部孔隙結構。相比與其他方法，CT掃描圖像法可以快速無損的獲取樣品孔隙空間結構參數，定量表征樣品的孔隙空間結構，建立三維可視化模型[9-12]。

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

現以桂東南地區的花崗巖殘積土為例，通過Micro-CT掃描與三維可視化軟件Avizo結合，建立孔隙網絡模型。目前，此項對接技術主要用于碳酸鹽巖、砂巖、頁巖、煤、天然氣水合物及多孔碳的研究，對于流體在花崗巖殘積土的研究極少。運用Avizo相關算法，提取花崗巖殘積土的表單元體及孔隙空間結構，并將表單元體劃分網格后導入多場耦合軟件Comsol進行花崗巖殘積土的滲流模擬，通過計算單相流體流經該數字巖心模型的絕對滲透率，與實測值對比，檢驗該模型準確性、代表性，并研究壓實度對花崗巖殘積土滲流特性的影響。

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

1.1 土樣壓實度選取及制備

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

1.2 CT掃描

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

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

1.3 圖像處理

1.3.1 濾波處理

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

1.3.2 閾值分割

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

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

1.3.3 表單元體提取

解決計算機內存容量、運行速度與數字巖心尺寸大小間矛盾是通過選取合適的代表性體積單元(representative elementary volume，REV)。小于REV尺度取得的土體物理特性波動明顯，而大于REV尺度取得的土體物理特性趨于穩定。因此，選擇大小合適的REV不僅能代表土體的物理特性，而且也能滿足計算機的條件。當REV尺寸與孔隙度的變化規律趨于穩定時，REV尺寸即為對應的最小單元體的邊長。通過分析REV尺寸與孔隙度的變化規律，發現壓實度為90%、100%、110%花崗巖殘積土的REV尺寸分別大于220×220×220、200×200×200、200×200×200體素時，土樣的孔隙度基本不隨單元體尺寸變化，如圖3所示。因此，壓實度為90%、100%、110%花崗巖殘積土選取表單元體尺寸分別為220×220×220、200×200×200、200×200×200體素，則對應的表單元體大小分別為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 微觀孔隙結構模型構建

三維可視化是將花崗巖殘積土孔隙與基質的分布結構以更加直觀的方式呈現。構建的孔隙網絡模型較為真實地還原花崗巖殘積土的孔隙分布及連通性。其過程是將分割后的圖像進行三維重構。基于選定的閾值和表單元體的尺寸，使用直體積渲染模塊，完成花崗巖殘積土三維重建，如圖4所示。

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

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

3 滲流模擬

3.1 網格剖分、優化

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

3.2 控制方程與邊界條件

將實際孔隙結構的數字巖心模型作為模擬基礎，并展開滲流模擬。其方法有兩類：計算流體力學的傳統方法(CFD)及格子玻爾茲曼Boltzman方法(LBM)。CFD是對介質模型進行質量、動量、能量守恒方程求解，并對控制方程及模型離散。LBM是一種基于分子運動理論模擬流體流動的方法，流體被看成有質量沒體積的微粒，通過模擬微粒在網格中的碰撞和移動，構建簡化的運動論模型來反映孔隙空間中流體滲流規律。但LBM方法計算量大、耗時長。因此采用CFD方法中不可壓縮Navier-Stokes方程模塊來進行孔隙空間結構的流體流動模擬。

(1)

(2)

(3)

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

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

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

3.3 絕對滲透率模擬

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

(4)

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

由表1、圖5～圖7可知，壓實度為90%、100%花崗巖殘積土在X、Y、Z方向的滲透率差異較小，但壓實度為110%花崗巖殘積土在3個方向的滲透率差異較大，X方向的滲透率是Z方向的9.27倍，其原因可能是當壓實度增加至110%時，分層壓實對Z方向上的滲透率影響較大；當壓實度增加至110%，花崗巖殘積土的滲透率改變2個數量級，其原因可能是當壓實度增加至110%時，花崗巖殘積土的孔隙形態呈扁平狀，孔隙直徑減小，土顆粒與其周圍溶液的物理化學作用會產生結合水膜，其對土顆粒會產生吸附作用，由于距離不同，故吸附力的強弱會有所不同。距離越近，結合水膜越穩定，在滲透中基本不發生運動。因此，結合水膜的存在會削弱土體的滲透性且從壓力云圖可以看出壓實度為110%花崗巖殘積土相對于壓實度為90%、100%花崗巖殘積土在3個方向的孔隙空間結構發生較大的變化，分布不均勻，進出口面積相差較大。壓實度為90%、100%、110%花崗巖殘積土在X、Y、Z方向上滲透率模擬值分別是實測滲透率0.89、0.53、1.29、1.74、1.71、0.62、2.02、0.88、0.22倍；其中壓實度為90%、100%、110%花崗巖殘積土滲透率模擬值平均值分別是實測滲透率0.9、1.35、0.77倍。表明孔隙結構的數字巖心模型具有準確性、代表性。存在差異主要是因為花崗巖殘積土在不同方向上存在一定的非均質性且Micro-CT的分辨率有限，無法識別小于分辨率的孔隙，從而會影響滲透率的計算。

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

3.4 壓力場分布

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

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

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

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

3.5 速度場分布

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

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

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

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

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

4 結論

以桂東南地區不同壓實度的花崗巖殘積土為研究對象，基于CT掃描試驗與數字圖像處理技術，構建其三維孔隙網絡模型，并基于CFD研究水在土體孔隙中的滲流規律，得到以下結論。

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

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

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

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