蔡沛辰， 闕 云， 李 顯

(福州大學(xué) 土木工程學(xué)院， 福建 福州 350108)

1 研究背景

土體孔隙的滲流問題在壩體工程、邊坡防護(hù)、基坑開挖等領(lǐng)域具有非常重要的研究意義[1-3]，由于土體滲流涉及水-氣兩相耦合問題，使得求解過程變得較為棘手[4]。而傳統(tǒng)的滲流研究多集中于宏觀尺度，用試驗(yàn)手段測得滲透系數(shù)來表征土體的滲流特性。但一方面這種研究無法直觀準(zhǔn)確地描述孔隙結(jié)構(gòu)特征、流體滲流路徑以及孔隙結(jié)構(gòu)與滲流特性之間的內(nèi)在聯(lián)系[5]，另一方面，試驗(yàn)方法耗時、耗力、耗資，且不利于重復(fù)性開展。因此，為克服上述試驗(yàn)方法的不足，一些細(xì)觀尺度下的數(shù)值仿真方法得到了廣泛應(yīng)用并迅速發(fā)展。

目前，細(xì)觀尺度下數(shù)值仿真手段主要包括以下兩大類：

(1)基于計算機(jī)編程的仿真手段。例如苗杰[6]借助Poreflow半代碼編程模擬器，對煤巖進(jìn)行了水-氣兩相流模擬，并將模擬結(jié)果與實(shí)測結(jié)果進(jìn)行了對比驗(yàn)證；徐宗恒等[7]采用格子玻爾茲曼方法純代碼編程，對土體孔隙進(jìn)行了整個滲流過程的定量化研究。

(2)基于專業(yè)有限元軟件的模擬手段，例如李輝等[8]從土石壩滲流角度出發(fā)，借助COMSOL軟件建立模型，對不同水位下的壩體浸潤線、滲透水壓力以及壩體內(nèi)比降等相關(guān)滲流參數(shù)進(jìn)行了研究；戚藍(lán)等[9]、湛文濤等[10]和邢皓楓等[11]從降雨強(qiáng)度角度出發(fā)，分別采用ABAQUS和SEEP軟件分析了降雨類型和降雨強(qiáng)度對邊坡滑移特征時空分布規(guī)律的影響，并詳細(xì)探討了降雨強(qiáng)度和歷時對路基土邊坡滲流場的影響情況。上述研究雖然已取得了豐碩的成果，但計算機(jī)編程仿真方法要求研究者不僅需具備較強(qiáng)的專業(yè)知識，且需熟悉相應(yīng)的編程語言，適用對象面相對較窄，而應(yīng)用最為廣泛的專業(yè)軟件模擬方法雖對程序語言能力要求較低，但還是較側(cè)重于針對專業(yè)的科研工作者使用，且重復(fù)性操作較為繁瑣、耗時。

鑒于此，本文以原狀花崗巖殘積土為研究對象，對其進(jìn)行工業(yè)CT掃描，建立反映真實(shí)土體孔隙結(jié)構(gòu)的幾何模型，采用虛擬仿真APP開發(fā)手段封裝土體細(xì)觀尺度下水-氣兩相滲流數(shù)值仿真平臺，對孔隙內(nèi)水-氣兩相流動過程進(jìn)行動態(tài)可視化研究，并驗(yàn)證了結(jié)果的真實(shí)性，以期探尋一種高效、快捷的研究手段，使科研工作者面對友好的自定義界面，從而提高科研效率。

2 CT掃描及土體細(xì)觀模型構(gòu)建

試驗(yàn)原狀土選自福建省福州市某地山坡，現(xiàn)場取樣如圖1所示。取樣點(diǎn)選取植被茂密的位置，取樣前先清理表面腐殖層，清理面積為100 cm×100 cm，厚度為20 cm，最終獲取尺寸為15 cm×15 cm×40 cm的土柱試樣。測試土樣的基本物理參數(shù)，結(jié)果如表1所示，其中，滲透率定義為：在一定壓差下，土體本身允許流體通過的能力，由于花崗巖殘積土滲透性較小，滲透率測定一般采用變水頭試驗(yàn)[12]。

表1 所選試樣的基本物理參數(shù)

圖1 試驗(yàn)原狀土現(xiàn)場取樣[5]

對獲取的原狀土采用英華公司C450KV高能量工業(yè)CT掃描儀進(jìn)行XZ面掃描，CT掃描工作電壓為450 kV，工作電流為63 mA，掃描最低分辨率為0.15 mm，之后得到一系列能真實(shí)體現(xiàn)土壤孔隙分布的CT掃描切片，如圖2所示。

通過MATLAB軟件中的imread和im2bw等函數(shù)功能，對CT掃描圖像進(jìn)行二值化處理，形成只包含黑白兩色的圖像，再基于小波變換對二值化圖像進(jìn)行降噪處理，除去圖像中孤立的噪點(diǎn)，最后選取中間位置的切片，截取其中孔隙連通性較好的區(qū)域作為模擬對象，如圖3所示，截取的區(qū)域大小為3 mm×1.65 mm。將圖像以數(shù)組的形式存儲在MATLAB中，借助COMSOL-MATLAB接口[13]將它們轉(zhuǎn)換為計算域，并作為計算幾何模型導(dǎo)入COMSOL中進(jìn)行仿真模擬，如圖4所示。

圖2 試驗(yàn)土樣二維掃描切片[5] 圖3 本文所采用的土樣二值化掃描圖像 圖4 進(jìn)行仿真模擬的土樣幾何模型

3 虛擬仿真平臺及APP開發(fā)

3.1 虛擬仿真平臺

COMSOL Multiphysics被譽(yù)為“第一款真正的任意多物理場直接耦合分析軟件”，是一款具備完善的理論基礎(chǔ)、集功能性、靈活性和實(shí)用性于一體，并可通過內(nèi)部預(yù)設(shè)的專業(yè)求解模塊方便用戶進(jìn)行應(yīng)用和進(jìn)一步開發(fā)拓展的軟件[14]。其中，COMSOL APP是通過高度專業(yè)化的用戶界面與COMSOL模型進(jìn)行交互的一種直觀而有效的研究方式。

3.2 仿真APP構(gòu)建過程

下面以滲流場研究中經(jīng)典的水-氣兩相流問題為例，詳細(xì)敘述如何通過基于COMSOL Multiphysics構(gòu)建的虛擬仿真平臺封裝搭建獨(dú)立運(yùn)行的仿真APP，便于科研工作者省時、省力地研究土體孔隙中的滲流規(guī)律[15]。

3.2.1 數(shù)值模擬平臺建立

(1)幾何模型導(dǎo)入。將建立的計算模型通過COMSOL-MATLAB接口轉(zhuǎn)換為計算域(圖4)，并作為計算幾何模型導(dǎo)入COMSOL中進(jìn)行仿真模擬。

(2)分別采用N-S(Navier-Stokes)方程和Level Set方法描述流體流動狀態(tài)和水-氣兩相界面的瞬時動態(tài)變化情況[16]。其中Level Set方法接口通過跟蹤水平集函數(shù)的等值線Φ來確定流體界面，Φ=0表示流體為氣相，Φ=1表示流體為水相，一般取Φ=0.5的等值線作為相界面?？刂痞档膫鬟f和重新初始化的方程為：

(1)

式中：Φ為水-氣兩相界面等值線函數(shù)；t為滲流時間，s；u為流速，m/s；γ為重新初始化參數(shù)，m/s；ε為界面厚度控制參數(shù)，m。

為使Level Set方程在計算過程中具有較高的穩(wěn)定性，公式(1)中的ε通常可取值為ε=hc/2，其中hc為界面區(qū)域的網(wǎng)格大小。此外，可通過下式確定全局密度ρ(kg/m3)和動力黏度μ(Pa·s)：

ρ=ρw+(ρg-ρw)Φ

(2)

μ=μw+(μg-μw)Φ

(3)

式中：ρw、ρg分別為水相和氣相的密度，kg/m3；μw、μg分別為水相和氣相的動力黏度，Pa·s。

計算步驟如下：

第1步，將構(gòu)建的多孔介質(zhì)幾何模型導(dǎo)入到COMSOL中，設(shè)置相應(yīng)的尺寸參數(shù)，并用多邊形矢量功能設(shè)置兩相滲流的初始界面位置，將多孔介質(zhì)模型分為兩部分。

第2步，對多孔介質(zhì)幾何模型兩部分分別添加研究者所需要的仿真材料屬性，驅(qū)替相添加為“水”，被驅(qū)替相添加為“空氣”。

第3步，采用水平集模塊和層流模塊進(jìn)行多物理場耦合計算，驅(qū)替相與被驅(qū)替相本構(gòu)關(guān)系均為牛頓流，且在動量方程中包含表面張力，邊界條件為不透水層，邊壁為無滑移條件。水相從上邊界進(jìn)入孔隙區(qū)域，以初始界面為界線，通過水壓力驅(qū)替氣相，到下邊界流出，整個過程中孔隙區(qū)域內(nèi)保持恒定的水壓差。

第4步，使用軟件自帶的“自由三角形網(wǎng)格”功能對模型進(jìn)行網(wǎng)格劃分，類型和尺寸設(shè)置為流體動力學(xué)和極細(xì)化(若劃分不成功，即由于多孔介質(zhì)模型孔隙結(jié)構(gòu)較為復(fù)雜或均質(zhì)性差，需要對網(wǎng)格調(diào)整為粗化)。

第5步，使用相初始化求解器和瞬態(tài)求解器共同進(jìn)行求解，求解器的時間步長為range(0～0.2 s,每10-5s計算1次)，點(diǎn)擊計算，經(jīng)過計算后可得到兩相滲流的動態(tài)可視化過程以及速度場、壓力場等示意圖。

(3)交互搭建APP。使用COMSOL中APP開發(fā)器功能對模型建立過程進(jìn)行交互APP封裝搭建。自主設(shè)計封裝完成的APP啟動界面如圖5所示。該系統(tǒng)主要包括“參數(shù)設(shè)置區(qū)”“邊界設(shè)定區(qū)”“功能操作區(qū)”和“結(jié)果展示區(qū)”等，詳細(xì)如圖6所示。其中，“參數(shù)設(shè)置區(qū)”和“邊界設(shè)定區(qū)”主要包含“模型導(dǎo)入”“材料屬性設(shè)定”“邊界條件選定”“壓差參數(shù)輸入”等；“功能操作區(qū)”主要包含“構(gòu)建模型”“網(wǎng)格剖分”“計算”“創(chuàng)建報告”和“幫助”等按鈕；“結(jié)果展示區(qū)”主要為模型可視化展示以及計算后得到兩相滲流過程、速度場等結(jié)果示意圖的展示。

圖5 自主設(shè)計封裝完成的APP啟動界面

3.2.2 仿真APP封裝

(1)主窗口設(shè)定。在主窗口下菜單欄中添加3個子項(xiàng)，并分別設(shè)置為“保存”“另存為”和“退出”，快捷方式鏈接為“CTRL+S”“CTRL+SHIFT+S”和“CTRL+E”。

(2)建立相應(yīng)功能的表單，并在主表單中完成APP運(yùn)行界面的整體布局規(guī)劃。在表單窗口建立4個表單，其中第1個表單作為主表單用于APP窗口的整體搭建布局設(shè)計，為防止表單過多，將“參數(shù)設(shè)置區(qū)”“功能操作區(qū)”和“結(jié)果展示區(qū)”布置在其上；其他3個表單分別為集合表單、信息表單和APP主要功能簡介表單，將信息表單和APP主要功能簡介表單及主表單復(fù)合到幾個表單中，以主表單為默認(rèn)窗口展示。

(3)程序編寫，添加相應(yīng)功能。對于所需輸入的參數(shù)及仿真過程在之前已完成，現(xiàn)可以直接鏈接到APP中，并為參數(shù)設(shè)置區(qū)、功能操作區(qū)和結(jié)果展示區(qū)添加相應(yīng)功能。

“參數(shù)設(shè)置區(qū)”的功能實(shí)現(xiàn)選擇全局定義命令下的所需參數(shù)進(jìn)行鏈接，可通過添加Method并輸入程序語句實(shí)現(xiàn)參數(shù)初始化功能；“邊界設(shè)定區(qū)”需要插入“選擇輸入”，并以“顯示”的方式與之前建模過程進(jìn)行鏈接，顯示于主表單的圖像窗口；“功能操作區(qū)”的功能實(shí)現(xiàn)選擇模型ROOT命令下的相關(guān)操作；“結(jié)果展示區(qū)”的功能實(shí)現(xiàn)選擇結(jié)果命令下的相關(guān)可視化顯示，且可在工具欄命令下添加圖像放大、縮小、打印等基本功能。

(4)APP搭建封裝完成。至此，所有按鈕及相關(guān)功能已封裝完成，通過compiler功能對其編譯處理，并可根據(jù)自身需求設(shè)置相應(yīng)的APP圖標(biāo)和啟動界面，最后輸出可供Windows、Linux和macOS三大系統(tǒng)使用的EXE格式的APP文件。

圖6 虛擬仿真APP運(yùn)行界面

4 結(jié)果與分析

基于CT掃描技術(shù)構(gòu)建的土體孔隙細(xì)觀模型，以水-氣兩相滲流為例(水相為驅(qū)替相，氣相為初始相)，借助虛擬仿真APP開發(fā)技術(shù)，將復(fù)雜的理論和物理場封裝起來，模擬不同時刻水-氣兩相滲流的動態(tài)可視化過程。為使得模擬結(jié)果更好地貼近于真實(shí)情況，土體孔隙滲流方向設(shè)為沿深度方向進(jìn)行，并設(shè)置模型初始壓差以Pin=110 Pa入滲、Pout=10 Pa流出，詳細(xì)邊界條件設(shè)置如圖7所示，其他相關(guān)參數(shù)如表2所示。

圖7 土體孔隙水-氣兩相滲流邊界條件設(shè)定

表2 土體孔隙水-氣兩相滲流計算參數(shù)

4.1 仿真數(shù)值模型驗(yàn)證

為驗(yàn)證所構(gòu)建數(shù)值模型的正確性，采用文獻(xiàn)[18]的方法，驗(yàn)證幾何模型中滲流是否符合Darcy定律。水流通過壓力進(jìn)行驅(qū)動，從上邊界進(jìn)入孔隙區(qū)域，到下邊界流出，并在孔隙區(qū)域內(nèi)保持恒定的水壓差100 Pa。整個孔隙區(qū)域XZ面的水流流速云圖如圖8所示。由圖8可以發(fā)現(xiàn)，在滲流路徑上存在一些高速區(qū)，經(jīng)計算得出整個孔隙區(qū)域的平均流速為0.041 4 m/s，平均體積流量為8.19×10-8m3/s。

圖8 土體孔隙區(qū)域XZ面模型流速云圖

通過改變模型邊界的水壓差計算相應(yīng)的體積流量，對其進(jìn)行參數(shù)化研究，圖9為滲流體積流量隨水壓差的變化趨勢。由圖9可看出，滲流體積流量隨施加的水壓差基本呈線性變化，符合Darcy定律，也驗(yàn)證了該計算模型的正確性。

圖9 土體孔隙水-氣兩相滲流體積流量隨水壓差的變化趨勢

4.2 水-氣兩相滲流過程

圖10為土體孔隙水-氣兩相滲流過程動態(tài)示意圖。由圖10(a)、10(b)可看出，水體會優(yōu)先選擇較大孔隙進(jìn)行滲流，之后才會選擇較窄孔隙。由圖10(c)可發(fā)現(xiàn)，滲流相在成圓度較高的孔隙處易出現(xiàn)“繞流”現(xiàn)象，通過計算得出此處成圓度達(dá)0.76，成圓度表征孔隙形狀趨近于圓的程度，定義為4π倍的孔隙面積與孔隙周長平方的比值[19]。再對比觀察圖10(c)、10(d)可以發(fā)現(xiàn)，滲流水體還存在“回流”現(xiàn)象，且部分死角孔隙水流不能進(jìn)入，其原因是孔隙中的水、氣均不可壓縮且不相溶，使孔隙中形成了只含有空氣的封閉空間。

圖10 土體孔隙水-氣兩相滲流過程動態(tài)示意圖

綜上所述，通過虛擬仿真APP平臺對于細(xì)觀尺度水-氣兩相滲流進(jìn)行模擬，能很好地顯示兩種不混溶流體之間的界面位置；水-氣兩相滲流過程存在大孔隙優(yōu)先流特征，且“繞流”現(xiàn)象一般易于出現(xiàn)在孔隙成圓度較高處；受水、氣相參數(shù)限制，部分死角孔隙形成了只含氣相的封閉空間。

4.3 水-氣兩相滲流速度場

在水-氣兩相滲流細(xì)觀速度場的基礎(chǔ)上，運(yùn)用COMSOL后處理模塊中的粒子追蹤方法繪制不同時刻流速場下的粒子追蹤圖像，如圖11所示(粒子起點(diǎn)控制參數(shù)表示為X：range(0,0.006,3)，Z：1.65，單位：mm)，圖中粒子軌跡線條的顏色與圖例相對應(yīng)，表示流速的大??；線條密度表示通過孔隙的水流流量大小。

對比圖11(a)和 11(b)可知，滲流路徑1-1′和2-2′相比其他路徑孔道的彎曲程度更低，滲流速度也相應(yīng)較大。再對比觀察圖11(c)和11(d)可以發(fā)現(xiàn)，t=4×10-3s時刻孔隙中的流量明顯大于t=8×10-4s時刻，同時水-氣兩相滲流存在“繞流”現(xiàn)象，出現(xiàn)“繞流”現(xiàn)象部位的速度相比其周圍速度更小，且這種現(xiàn)象隨著滲流時間的增加而更加顯著，如圖11(d)中圓圈標(biāo)出的“繞流”部位的平均“繞流”流速僅為0.005 8 m/s，而其周圍的滲流平均流速為0.087 m/s，兩者相差14倍。

圖11 土體孔隙水-氣兩相滲流速度場不同時刻的粒子追蹤圖像

上述現(xiàn)象表明，水-氣兩相滲流流速大小受孔道彎曲程度控制，且流體優(yōu)先選擇連通性較好的孔道進(jìn)行滲流。水-氣兩相滲流中發(fā)生“繞流”部位的流速遠(yuǎn)小于其周圍的滲流流速。

4.4 水-氣兩相滲流壓力場

圖12為土體孔隙水-氣兩相滲流不同時刻壓力場分布示意圖。由圖12可知，在初始階段，孔隙所受的壓力沿兩相滲流方向逐漸減小，且孔徑越小則壓力越大。但隨著滲流歷時的持續(xù)增加，孔隙所受到的壓力分布范圍沿滲流方向逐漸擴(kuò)大，且孔徑大、連通性好的孔道的壓力值越大。此外，由于土體內(nèi)部的孔隙結(jié)構(gòu)復(fù)雜，造成了壓力場分布的不均勻性。

圖12 土體孔隙水-氣兩相滲流不同時刻壓力場分布示意圖

4.5 討 論

表3列出了本文與其他學(xué)者關(guān)于土(巖)體孔隙水-氣(油)兩相流的研究結(jié)果對比。由表3可以得出以下結(jié)論：

(1)目前國內(nèi)外研究學(xué)者對兩相流的研究多集中于巖石和理論模型領(lǐng)域，少有涉及真實(shí)原狀土體的兩相滲流研究。

(2)巖石類的水-氣兩相流研究一般均可得出優(yōu)勢流和指近現(xiàn)象，而除此之外，本文還發(fā)現(xiàn)土體孔隙中存在“回流”和“繞流”現(xiàn)象。

(3)本文與前人研究角度不同，不局限于直接研究土體兩相滲流場的整體流速大小，而是研究了“繞流”區(qū)域內(nèi)流速分布與周圍區(qū)域流速的大小關(guān)系。此外，還研究了兩相滲流動態(tài)過程中壓力場的分布變化情況，發(fā)現(xiàn)土體孔隙兩相滲流壓力場分布存在不均勻性的特點(diǎn)。

表3 本文與其他學(xué)者關(guān)于土(巖)體孔隙水-氣(油)兩相流的研究結(jié)果對比

5 結(jié) 論

本文以原狀花崗巖殘積土水-氣兩相滲流研究為例，借助數(shù)值模型構(gòu)建虛擬仿真APP，將復(fù)雜的理論和物理場進(jìn)行封裝，并驗(yàn)證了虛擬仿真APP開發(fā)在土體孔隙細(xì)觀兩相滲流研究中的可行性。得出以下結(jié)論：

(1)水-氣兩相滲流過程存在大孔隙優(yōu)先流特征，且“繞流”現(xiàn)象一般易出現(xiàn)于孔隙成圓度較高的部位。受水、氣相參數(shù)的限制，部分死角孔隙形成了只含氣相的封閉空間。

(2)兩相滲流流速主要受到孔道迂回程度的影響，孔徑大、彎曲程度低的孔隙中滲流速度相對較大，存在明顯的“優(yōu)勢通道”。發(fā)生“繞流”現(xiàn)象部位的流速遠(yuǎn)小于其周圍的滲流流速。

(3)土體孔隙兩相滲流壓力場分布存在不均勻性。在滲流初期，孔隙所受壓力沿滲流方向逐漸減小，且孔徑越小則壓力越大。但隨著滲流歷時的增加，孔隙壓力沿滲流方向增大，且孔徑大、連通性好的孔隙所受壓力越大。

(4)土體孔隙兩相流是一個復(fù)雜的動態(tài)過程，本研究中模擬了土體孔隙中水-氣兩相滲流過程，關(guān)注點(diǎn)為土體孔隙結(jié)構(gòu)特征對兩相流的影響，但未考慮水、氣相互溶合性及基質(zhì)域滲流的問題，故將此作為今后的研究方向。本文研究成果能夠?yàn)樵瓲钔馏w細(xì)觀滲流機(jī)理的定量化研究提供一定幫助，所采用的研究思路對其他類似巖土體的滲流研究也具有一定的借鑒意義。