非飽和花崗巖殘積土水-氣兩相驅替過程數(shù)值模擬

蔡沛辰，闕 云，李 顯

（福州大學土木工程學院，福建 福州 350108）

非飽和土廣泛存在于自然界中[1]。與飽和土不同，非飽和土是一種三相土，固相（土體基質）、液相（水溶液）和氣相（空氣）同時存在。事實上，水流在非飽和土體中的入滲過程實質上是水在下滲過程中驅替空氣的兩相流問題。孔隙細觀尺度下的兩相流研究具有十分廣泛的工程背景，石油開采、廢棄物處理、地下水有機污染及泄水建筑中的摻氣水流等都與細觀多相滲流機理有著密切的聯(lián)系[2?4]。傳統(tǒng)的兩相流研究多采用試驗方法，如玻璃刻蝕、薄片平面驅替等二維方法[5]，巖心驅替物理試驗、宏觀驅替CT 掃描等三維方法[6?7]，但試驗方法研究成本較高，操作較繁瑣困難，且精度一般無法滿足工程需要，更不利于開展重復性研究。因此，為克服試驗方法的不足，細觀尺度下的數(shù)值模擬方法得到了迅速發(fā)展和應用。

數(shù)值模擬方法利用計算機刻畫細觀孔隙結構，模擬水驅氣的整個動態(tài)過程，具有重復、方便和可視化等優(yōu)點。目前，細觀尺度下多相流的數(shù)值模擬方法主要包括：格子Boltzmann、孔隙網(wǎng)絡模型、相場等方法。如：毛歡[8]基于格子Boltzmann 方法，利用Shan-Chen 模型在孔隙尺度上對兩相不混溶驅替過程進行了模擬；Qren 等[9]基于孔隙網(wǎng)絡模型方法，通過構建大量砂巖孔隙網(wǎng)絡模型，進行了多相滲流模擬，并通過模擬結果預測了多相流的相對滲透率；陳民鋒等[10]建立油水兩相三維孔隙網(wǎng)絡模型，模擬了儲層巖心的初次油驅和二次水驅過程；朱光譜等[11]基于相場方法，對孔隙尺度下表面活性劑兩相流進行了模擬研究。上述相關研究已取得豐碩成果，但多集中于巖石多相流的研究，且研究對象與實際孔隙結構存在一定差異；應用較廣泛的孔隙網(wǎng)絡模型方法無法直觀展現(xiàn)任一時刻各質點流速變化情況及無法呈現(xiàn)相界面移動狀態(tài)。Level Set 方法最早由Osher 等[12]提出，現(xiàn)已用于兩相流的數(shù)值模擬領域[13?15]，但受限于孔隙結構理想化或油氣/水兩相驅替，且未涉及土體細觀孔隙內的水驅氣兩相模擬。鑒于此，本文選用花崗巖殘積土作為研究對象，對其進行工業(yè)CT 掃描，建立二維細觀孔隙模型，借助COMSOL Multiphyscis 有限元軟件，采用Level Set 方法進行數(shù)值模擬，刻畫土壤內部各孔隙區(qū)域水-氣兩相驅替的動態(tài)可視化過程，并分析驅替過程中的驅替速度、阻力及效率的規(guī)律。

1 CT 掃描及細觀模型構建

1.1 CT 掃描及圖像處理

試驗原狀土選自福州市某地山坡，現(xiàn)場取樣如圖1所示。為保證取樣過程盡量不破壞原狀土孔隙結構，首先將取樣位置周圍土壤挖出，保留核心土，之后將中間預留的土體上部削成與模型盒相當?shù)拇笮。賹⒛Ｐ秃新蛳聣喝耄来沃貜拖鲏海敝镣馏w到達模型盒底后再向下超削約5 cm 深度，以便將土柱從根部切斷，最終獲取尺寸為15 cm×15 cm×40 cm 的土柱試樣。

圖1 現(xiàn)場取樣Fig.1 Field sampling

為保證原狀土體孔隙CT 掃描圖像的真實性，對其進行如下處理：首先，將獲取的土柱試樣用保鮮防水材料封裝，防止土體內部水分過分蒸發(fā)；其次，將封裝后的試樣放置于定制的木質模板箱中，防止搬運途中的擾動對其孔隙結構造成較大影響；最后，運至實驗室后，拆除模板和封裝材料，進行工業(yè)CT 掃描試驗。

CT 掃描試驗設備為C450KV 高能量工業(yè)CT，工作電壓為450 kV，電流為63 mA，掃描最低分辨率為0.15 mm。對獲取的圖像，通過Image J 中的Threshold功能進行二值化處理，形成只包含黑白兩色的圖像，再利用中值濾波對其進行降噪，除去圖像中孤立的噪點，最終得到一系列15 mm×15 mm 的二維掃描切片(圖2）。考慮到計算時間的因素，從原始圖像中選取孔隙連通性較好的一部分區(qū)域作為模擬對象，并將其等分為4 部分，局部放大如圖3所示，模型a—d 尺寸均為3.0 mm×1.6 mm。

圖2 二維掃描切片F(xiàn)ig.2 2D scanning slice

圖3 模型a—d 二值化圖像Fig.3 Model a—d binarized image

1.2 計算幾何模型

基于上述步驟，類似處理其他圖像以獲得不同孔隙幾何模型。最后，將這些圖像以數(shù)組的形式存儲在MATLAB 中，借助COMSOL-MATLAB 接口[16]將它們轉換為計算域，并作為計算幾何模型導入COMSOL中進行仿真模擬。以圖3 中模型b 為例，構建的計算幾何模型如圖4（a）所示。此外，為清楚查看網(wǎng)格剖分質量分布情況，將圖4（a）中紅色圈出部位放大，得到圖4（b），可以看出：模型b 的網(wǎng)格質量大多處于0.8 左右，質量分布較為均勻，也為后續(xù)數(shù)值模擬結果的精確度奠定了基礎。

圖4 模型b 的計算幾何模型及網(wǎng)格質量分布圖Fig.4 Computational geometric model and mesh quality distribution map of model b

2 水-氣兩相驅替仿真模擬

2.1 控制方程

2.1.1 納維爾·斯托克斯方程

采用不可壓縮的納維爾·斯托克斯（N-S）方程描述流體流動[17]，控制方程如下：

式中：ρ—流體密度/（kg·m?3）；

u—流速/（m·s?1）；

p—壓力/Pa；

I—單位矩陣；

μ—流體動力黏度/（Pa·s）；

F—體積力/（N·m?3）；

g—重力加速度/（m·s?2）。

2.1.2 界面張力

界面張力Fst由下式定義：

式中：σ—表面張力系數(shù)；

n—界面單位法向量；

δ—狄拉克函數(shù)；

Φ—水-氣兩相界面等值線。

2.1.3 Level Set 方法

“Level Set”接口通過跟蹤水平集函數(shù)的等值線來確定流體界面，等值線Φ= 0.5 決定界面的位置。控制Φ的傳遞和重新初始化的方程為：

式中：γ—重新初始化參數(shù)/（m·s?1）；

ε—界面厚度控制參數(shù)/m。

由于水平集函數(shù)是一個平滑階躍函數(shù)，因此可通過下式確定全局密度和動力黏度：

式中：ρw—水相密度/（kg·m?3）；

μw—水相動力黏度/（Pa·s）；

ρg—氣相密度/（kg·m?3）；

μg—氣相動力黏度/（Pa·s）。

2.2 材料屬性

采用水作為驅替相，空氣作為被驅替相，具體材料屬性見表1，其中界面張力的設定見文獻[18]。

表1 材料屬性Table 1 Material properties

2.3 邊界及初始條件

水-氣兩相驅替模擬中，依據(jù)構建的幾何模型，進行邊界條件設定：上下邊界為進出口，左右邊界為不透水層，壁為無滑移條件，初始界面設置為t= 0 時刻水相和氣相的接觸面。以模型b 為例(圖5)，其他模型的設定與模型b 類似。其中，Pin為入口水壓，Pout為出口水壓。

圖5 邊界條件設定Fig.5 Boundary condition setting

3 結果與討論

3.1 數(shù)值模型驗證

為驗證文中數(shù)值模型的正確性，采用體積守恒方法[19]驗證幾何模型中水-氣兩相驅替結果是否正確。由于多孔介質孔隙區(qū)域體積一定，故每一時刻氣相與水相的總體積是恒定不變的。對于二維多孔介質孔隙模型而言，應采用孔隙區(qū)域面積來驗證，即每一時刻氣相與水相所占的孔隙面積是恒定不變的，模型a—d 孔隙總面積依次為8.56×10?7，1.83×10?6，1.75×10?6，8.87×10?7m2。

t= 0 時刻，氣相存在于初始界面以下區(qū)域，水相存在于初始界面以上區(qū)域。水相從上邊界進入孔隙區(qū)域，以初始界面為界限，通過水壓力驅替氣相，到下邊界流出，并在孔隙區(qū)域內保持恒定的水壓差。其中，考慮到運算效率及后期驅替可觀察度等因素，以入口水壓1.1 kPa、出口水壓0.1 kPa 為例進行驅替運算仿真。圖6 為模型a—d 水相、氣相占孔隙區(qū)域面積隨時間變化的曲線。從圖中可以觀察到：4 個模型中，水-氣兩相在孔隙區(qū)域中的面積之和恒定不變，且都等于各自模型的孔隙總面積，符合守恒定律，也驗證了該數(shù)值模型的正確性。

圖6 水、氣相占孔隙區(qū)域的面積隨時間變化曲線Fig.6 The area curve of water and gas in the pore area with time

3.2 水-氣兩相驅替

3.2.1 動態(tài)驅替過程分析

采用Level Set 方法模擬不同時刻水驅氣的動態(tài)可視化過程，圖7 為模型a—d 不同時刻下水-氣兩相驅替過程示意圖，其中綠色為氣相，藍色為液相。

從圖7（a）可以發(fā)現(xiàn)：流體會優(yōu)先選擇較大孔隙進行驅替，之后才會選擇較窄孔隙，且在成圓度較高的孔隙處易出現(xiàn)“繞流”現(xiàn)象。計算知此處成圓度高達0.83，成圓度表征孔隙形狀趨近于圓的程度，定義為：4π 倍的孔隙面積與孔隙周長平方的比值[20]。由圖7（b）可看出：驅替流體存在“回流”現(xiàn)象，且部分死角孔隙流體不能進入，其原因為兩種流體均不可壓縮且不相溶，水相驅替使其形成了只含有空氣的封閉空間。由圖7（c）可知：筆直且較大的孔道，驅替速度相對較高，如通道1 明顯快于2 和3。再結合圖7（d）可見：在孔道相互交錯區(qū)域，流體并不直接相融合，而是經過一段時間后才匯合。

綜上，對于細觀尺度水-氣兩相驅替模擬，Level Set 方法能很好地捕捉兩種不混溶流體之間的界面位置。水-氣兩相驅替過程存在大孔隙優(yōu)先流特征，且“繞流”現(xiàn)象一般易于出現(xiàn)在孔隙成圓度較高處。受水-氣相參數(shù)限制，部分死角孔隙形成了只含氣相的封閉空間。筆直且較寬的孔道，驅替速度相對較高，孔道相互交錯區(qū)域，多條流體需經過一段時間后才能融合。

3.2.2 驅替速度可視化分析

運用后處理模塊中的速度可視化功能，以模型b 為例，繪制不同時刻水-氣兩相驅替過程速度可視化圖(圖8)。可以看出：多孔介質模型眾多孔隙中僅存在少數(shù)組成的主要滲流通道，且主要通道普遍有著孔隙半徑大、孔道筆直的特征。由各不同時刻速度分布云圖可知，流速較大處多出現(xiàn)在孔道相對較窄或相交匯處，且在驅替壓差不變的條件下為5×10?4，1×10?3，1.5×10?3，2×10?3s 時，最大流速分別為0.878，0.984，1.03，1.05 m/s，即隨驅替時間延長，流速有逐漸增大趨勢。上述現(xiàn)象表明：水-氣兩相驅替速度受孔道迂回程度控制，且流體優(yōu)先選擇主要通道進行滲流，存在明顯的“優(yōu)勢通道”。驅替速度與驅替時間以先快后緩的特征呈正相關變化，最大增長率為10.77%，最小僅為1.90%。

圖9 為不同孔道橫截面速度分布圖，其中孔道a、b 位置如圖7（b）所標注。從圖中可以看出：孔道b 中速度最大位置出現(xiàn)在中心位置x=3.745 mm 處，而孔道a 中速度最大位置在4.815 mm 處，并不在中心位置4.775 mm，存在右偏現(xiàn)象，與前人研究結果有所不同[13]。分析原因是：孔道a 所處位置的上方流體流入的方向為右上到左下，相比孔道b，流體會偏向孔道右側通過。此外，由圖9（a）可知驅替持續(xù)時間越長，驅替速度越大，與圖8 相對應，但孔道b 流速最大出現(xiàn)在5×10?4s 處，到1×10?3s時流速反而有所降低，不符合前文所得結論。原因分析如下：孔道b 在1×10?3s 時，出現(xiàn)了“回流”現(xiàn)象，如圖7（b）所示，使得通過孔道b 的流量驟降，進而使其速度減小，降低幅度為21.62%。

圖7 不同時刻水-氣兩相驅替過程動態(tài)示意圖Fig.7 Dynamic schematic diagram of the water-gas two-phase displacement process at different times

圖8 不同時刻水-氣兩相驅替過程速度可視化圖Fig.8 Visualized diagram of water-gas two-phase displacement process speed at different times

圖9 不同孔道橫截面速度分布圖Fig.9 Velocity distribution diagram of different tunnel cross sections

3.2.3 驅替阻力場分析

圖10 描述了不同時刻水-氣兩相驅替阻力分布情況。可以發(fā)現(xiàn)：驅替阻力隨驅替時間的延長而增加，且越靠近孔壁區(qū)域，紅色程度越大，尤其是狹窄區(qū)域，表明驅替阻力最大出現(xiàn)在孔壁處，且孔道越窄，阻力越大。

圖10 不同時刻水-氣兩相驅替阻力分布情況Fig.10 Resistance distribution of water-gas two-phase displacement at different times

3.2.4 驅替效率分析

水-氣兩相驅替效率Et定義為：t時刻下，孔隙區(qū)域內水相飽和度與水-氣兩相飽和度的比值。

式中：Swt—t時刻水相飽和度；

Sgt—t時刻氣相飽和度。

為探究不同驅替壓差對驅替過程的效率影響情況，以0.1 kPa 為初始驅替壓差，分三階段加壓，每階段加壓0.3 kPa，進行驅替模擬，并得結果如圖11（a）所示。從圖中可以發(fā)現(xiàn)：0.1～0.4 kPa、0.4～0.7 kPa、0.7～1.0 kPa 驅替效率分別增加了25.49%、6.25%和1.47%，說明驅替效率隨驅替壓差的增大而增大，且第一階段加壓增速效果更加顯著。

圖11（b）為相同驅替壓差1 kPa 下不同孔隙結構

圖11 不同情況下驅替率曲線Fig.11 Displacement rate curve under different conditions

模型驅替效率的對比曲線。計算可知，模型a—d 孔隙率分別為26.59%、38.14%、36.43%和18.48%。由圖可知：初期，孔隙率最小的模型d，驅替效率最為顯著；之后，模型c 急劇上升，而模型a、b、d 增速明顯減緩。結合圖7 分析原因：模型c 雖孔隙率較大，但孔隙結構并不復雜，其孔隙均質性較好，無明顯狹窄、不連通的孔隙。除此之外，根據(jù)文獻[21]可知，驅替效率還受到孔隙結構形態(tài)、孔徑、連通率等因素影響。

上述分析表明：驅替效率與驅替壓差成正比關系，且初期加壓增速效果顯著，可達25.49%，而后期僅為1.47%。驅替初期，相比孔隙率大的模型而言，孔隙較小模型的驅替效率較大，而隨驅替時間持續(xù)，驅替效率還將受孔隙結構特征的影響。

3.3 與其他學者兩相驅替結果對比

表2 列出了本文與其他學者對兩相流驅替研究結果的對比。從表2 可以得出：

表2 兩相流驅替研究結果對比Table 2 Comparison of research results of two-phase flow displacement

（1）目前學者對兩相流的研究多集中于巖石領域，少有提及土體的兩相驅替研究。

（2）巖石類的驅替過程一般都可得出優(yōu)勢流和指近現(xiàn)象，而本文還發(fā)現(xiàn)其存在“回流”和“繞流”現(xiàn)象。

（3）前人對驅替速度的研究多集中于速度大小分布上，本文除上述研究，還研究了孔道位置與速度的關系，得出孔道中心流速受孔隙結構影響，并非一直最大，且“回流”和“繞流”現(xiàn)象會使流速驟降。

（4）驅替阻力和驅替效率前人研究較少，提及的都只限于分析，而未做深入研究，本文一方面通過分析水-氣兩相驅替阻力場分布情況，得到驅替阻力最大出現(xiàn)在孔壁處，孔道越窄，阻力越大；另一方面，研究了壓差、孔隙結構與驅替效率的關系，得出驅替效率與壓差成正比，且初期加壓效果顯著。

4 結論

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

（2）驅替速度主要受孔道迂回程度控制，筆直、較寬孔道，驅替速度相對較高，存在明顯的“優(yōu)勢通道”，且與時間以先急后緩的特征呈正相關變化，最大增速率為10.77%，最小僅為1.90%。孔道橫截面速度大小分布與孔隙結構有關，“回流”和“繞流”現(xiàn)象會使驅替速度驟降，降低幅度可達21.62%。

（3）驅替阻力隨驅替時間的延長而增加，且最大出現(xiàn)在孔壁處，孔道越窄，阻力越大。

（4）驅替效率與驅替壓差成正比關系，初期加壓增速效果顯著，可達25.49%，后期僅為1.47%。驅替初期，相比孔隙率大的模型而言，孔隙較小模型的驅替效率較大，而隨驅替時間持續(xù)，驅替效率還將受孔隙結構特征的影響。

（5）天然情況下，水-氣兩相驅替過程中的驅替效率不僅受土體孔隙率的影響，還受孔隙連通率、孔喉大小等定量化參數(shù)的影響，而本文只涉及到不同模型間的孔隙率與驅替效率的關系。未來將基于編程法在非飽和土體中提取各種孔隙結構參數(shù)，進行定量化表征，并綜合、多因素考慮其與驅替效率之間的關系。