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

蔡沛辰，闕 云，李 顯

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

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

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

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

1.1 CT 掃描及圖像處理

試驗(yàn)原狀土選自福州市某地山坡，現(xiàn)場取樣如圖1所示。為保證取樣過程盡量不破壞原狀土孔隙結(jié)構(gòu)，首先將取樣位置周圍土壤挖出，保留核心土，之后將中間預(yù)留的土體上部削成與模型盒相當(dāng)?shù)拇笮?，再將模型盒慢慢向下壓入，依次重?fù)削壓，直至土體到達(dá)模型盒底后再向下超削約5 cm 深度，以便將土柱從根部切斷，最終獲取尺寸為15 cm×15 cm×40 cm 的土柱試樣。

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

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

CT 掃描試驗(yàn)設(shè)備為C450KV 高能量工業(yè)CT，工作電壓為450 kV，電流為63 mA，掃描最低分辨率為0.15 mm。對獲取的圖像，通過Image J 中的Threshold功能進(jìn)行二值化處理，形成只包含黑白兩色的圖像，再利用中值濾波對其進(jìn)行降噪，除去圖像中孤立的噪點(diǎn)，最終得到一系列15 mm×15 mm 的二維掃描切片(圖2）。考慮到計(jì)算時(shí)間的因素，從原始圖像中選取孔隙連通性較好的一部分區(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 計(jì)算幾何模型

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

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

2 水-氣兩相驅(qū)替仿真模擬

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 決定界面的位置?？刂痞档膫鬟f和重新初始化的方程為：

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

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

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

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

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

ρg—?dú)庀嗝芏?（kg·m?3）；

μg—?dú)庀鄤恿︷ざ?（Pa·s）。

2.2 材料屬性

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

表1 材料屬性Table 1 Material properties

2.3 邊界及初始條件

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

圖5 邊界條件設(shè)定Fig.5 Boundary condition setting

3 結(jié)果與討論

3.1 數(shù)值模型驗(yàn)證

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

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

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

3.2 水-氣兩相驅(qū)替

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

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

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

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

3.2.2 驅(qū)替速度可視化分析

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

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

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

圖8 不同時(shí)刻水-氣兩相驅(qū)替過程速度可視化圖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 驅(qū)替阻力場分析

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

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

3.2.4 驅(qū)替效率分析

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

式中：Swt—t時(shí)刻水相飽和度；

Sgt—t時(shí)刻氣相飽和度。

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

圖11（b）為相同驅(qū)替壓差1 kPa 下不同孔隙結(jié)構(gòu)

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

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

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

3.3 與其他學(xué)者兩相驅(qū)替結(jié)果對比

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

表2 兩相流驅(qū)替研究結(jié)果對比Table 2 Comparison of research results of two-phase flow displacement

（1）目前學(xué)者對兩相流的研究多集中于巖石領(lǐng)域，少有提及土體的兩相驅(qū)替研究。

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

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

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

4 結(jié)論

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

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

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

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

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