超臨界CO2
--水兩相流與CO2毛細捕獲：微觀孔隙模型實驗與數值模擬研究1)

胡 冉 陳益峰 萬嘉敏 周創兵

?(武漢大學水資源與水電工程科學國家重點實驗室，武漢430072)?(勞倫斯伯克利國家實驗室地球能源科學部，美國伯克利94720)

超臨界CO2
--水兩相流與CO2毛細捕獲：微觀孔隙模型實驗與數值模擬研究1)

胡 冉?,?陳益峰?,2)萬嘉敏?周創兵?

?(武漢大學水資源與水電工程科學國家重點實驗室，武漢430072)?(勞倫斯伯克利國家實驗室地球能源科學部，美國伯克利94720)

CO2毛細捕獲機制是CO2地質封存中的關鍵科學問題，然而有關孔隙尺度下(微米極)超臨界CO2毛細捕獲的研究較少.采用高壓流體--顯微鏡--微觀模型實驗裝置，開展超臨界CO2條件(8.5MPa,45?C)下CO2驅替水(排水)和水驅替CO2(吸濕)實驗，采用高分辨率照相機采集CO2水兩相流運動圖像，并借助光學顯微鏡直接觀測孔隙尺度下CO2毛細捕獲特征.同時，采用計算流體動力學方法對實驗過程進行三維數值模擬.數值模擬不僅反映了實驗過程中兩相流驅替鋒面的推進過程，還刻畫了孔隙尺度下被捕獲的CO2液滴/團簇三維空間形態特征.最后，基于數值模擬給出了CO2初始飽和度與殘余飽和度曲線，即毛細捕獲曲線，并對比分析了3種毛細捕獲曲線預測模型(即Jurauld模型、Land模型和Spiteri模型)的優劣.分析表明，Jurauld模型的描述能力稍優于Land模型，Spiteri模型的描述能力較弱.由于Land模型只需單個參數，且參數具有明確的物理意義，因此在實際工程中，建議優先采用Land模型.

CO2地質封存,微觀模型,兩相流,數值模擬,毛細捕獲

引言

大量研究表明，化石燃料燃燒產生的CO2超量排放是全球氣候變暖的主要原因.目前，CO2地質封存(geological carbon sequestration,GCS)被廣泛認為是緩解CO2引起溫室效應的有效途徑之一.該方法將CO2注入地下以達到長期穩定、安全封存的目的.CO2地質封存場所包括衰竭油氣藏、不可開采煤層和地下咸水含水層[1].世界氣候變化委員會研究報告指出[2]，全世界范圍內CO2在咸水含水層中的封存量可占2005年—2050年總排放量的20%～50%.我國鹽水盆地眾多，CO2地質封存潛力巨大[34].當CO2注入地下至少800m深的咸水層中，CO2以超臨界狀態滲入巖石孔隙中，孔隙水被CO2驅替；當CO2停止注入時，地下水回流，原先注入到孔隙中的CO2一部分被水的回流驅替，而殘余部分由于毛細力和孔隙結構的雙重作用，被捕獲在孔隙中.這種捕獲機制被稱為毛細捕獲(capillary trapping或residual trapping)[5].在GCS中，毛細捕獲意義尤其重要：一方面，由于毛細捕獲主要依靠巖石孔隙間的毛細力捕獲CO2，對巖石的完整性要求不高，從而提高了地質封存的安全性[6]；另一方面，CO2以團簇形態位于孔隙中，CO2、水、礦物三者的接觸，促進了溶解和化學反應的進行，從而對CO2超長期礦物捕獲具有重要影響.

孔隙介質中CO2毛細捕獲量(飽和度)受控于CO2水兩相流過程，與介質的濕潤性[710]、孔隙結構[1012]、CO2注入速率[1314]與初始飽和度等因素密切相關.CO2初始飽和度與殘余飽和度的關系稱為毛細捕獲曲線，該曲線具有重要的工程意義.國內外研究人員對CO2毛細捕獲機制開展了大量實驗與數值模擬研究.在實驗方面，目前的研究手段主要包括：(1)基于CT(computed tomography)掃描的巖芯驅替或玻璃珠驅替實驗[1517];(2)基于照相機/顯微鏡觀測的微觀模型實驗[1820].Al-Raoush[8]研究了介質濕潤性對毛細捕獲量及其分布的影響，結果表明介質的親水性愈強，毛細捕獲量愈大，被捕獲的液滴尺寸愈大.Chaudhary等[10]研究介質濕潤性與孔隙孔喉比對CO2毛細捕獲的影響，實驗結果表明：親水性愈強，孔喉比愈大，CO2的捕獲量愈高.Andrew等[21]針對石灰巖開展了CO2水排驅實驗，實驗結果表明，被毛細捕獲的CO2液滴/團簇，其尺寸呈指數分布特征.武愛兵等[22]通過超臨界CO2(sc CO2)排驅飽水砂巖巖芯實驗，研究了鹽水濃度對CO2毛細量的影響，并指出兩相流的排驅形式可分為活塞式、攜帶式和溶解式.Niu等[17]研究了溫度、壓力和咸水濃度對砂巖孔隙中CO2殘余總量與分布特征的影響.在上述實驗研究中，盡管巖芯試樣能夠較精確反映巖石的三維孔隙結構，然而由于CT掃描精度有限，較難獲取微米尺度下CO2毛細捕獲的細部特征.

另一方面，基于激光刻蝕技術制備而成的微觀孔隙模型，可以采用顯微鏡和照相機直接觀測微米尺度下的兩相流動狀態，能夠采集到比CT掃描更高分辨率的流體運動圖像[2325].研究人員基于微觀模型實驗研究了驅替速率與流體黏滯系數[1820]、介質濕潤性[26]和孔隙結構[27]對驅替模式和毛細捕獲的影響.Lenormand等[20]和Zhang等[28]研究了兩種互不混溶流體驅替過程中的黏性指進、毛細指進與穩定驅替模式，表明驅替模式受控于兩個特征參數，即毛細數和黏度比.Cottin等[29]則研究了微觀模型濕潤性對指流形態的影響，表明當介質由親水性過渡為中等親水性時，指流形態向穩定的方向發展.但以上研究較少關注超臨界CO2條件下CO2的毛細捕獲特征.

除實驗研究外，數值模擬也在孔隙尺度下的兩相流運動規律與毛細捕獲機制研究方面發揮了重要作用.常用的數值模擬方法有孔隙網絡模型方法[30]、格子玻爾茲曼方法[31]、光滑粒子流體動力學方法[32]和計算流體動力學方法[3334].其中，前3種方法所需的模型參數較多，且模型參數不獨立，因此在計算之前，需要對參數進行標定[35].計算流體動力學方法的優點在于模型參數較少，無需標定，但對計算資源的要求較高.

本文基于高壓流體--顯微鏡--微觀模型觀測系統，研究超臨界CO2條件(8.5MPa、45?C)下CO2水兩相流運動規律以及孔隙尺度下 CO2毛細捕獲特征，采用直接求解納維--斯托克斯方程模擬兩相流運動過程，并與實驗結果進行對比研究.最后，基于數值模擬得到的CO2初始飽和度與殘余飽和度曲線(毛細捕獲曲線)，對比分析了3種殘余飽和度預測模型的優劣.

1 實驗裝置與實驗流程

1.1 高壓--顯微鏡--微觀模型實驗裝置

采用如圖1(a)所示的實驗裝置開展超臨界CO2水驅排實驗，該裝置主要由微觀模型觀測子系統、流量控制子系統、溫度控制子系統和CO2水反應器子系統四大部分組成.

(1)微觀模型觀測子系統由微觀模型安裝平臺、照相機(Carl Zeiss,AxioCam MRc5)和倒置的光學顯微鏡(Carl Zeiss,Observer Z1.m)組成.微觀模型長20mm，寬10mm，孔隙喉道寬度為50μm，孔隙深度為 40μm，孔隙總體積 (pore volume)為 1.97μL，如圖1(b).微觀模型由特定的螺絲安裝在平臺上.照相機采集微觀模型全范圍(20mm×10mm)內的流體流動狀態圖像(1幀/秒)，顯微鏡采集孔隙尺度(微米級)下的流體流動圖像.照相機和顯微鏡通過數據線與電腦連接，采集的圖像實時傳入電腦.

(2)流量控制子系統由高精度高壓流量泵(Teledyne ISCO,500HP×2和65HP)組成.如圖1(a)所示，泵A以恒定流量向微觀模型入口注入超臨界CO2，泵B以恒定流量注入水，泵C連接微觀模型出口，確保流體系統為恒定壓力(8.5MPa).

(3)CO2水反應器子系統由攪拌反應器 (Parr,model 4560)、泵A和CO2氣缸組成.盡管水和CO2是兩種互不混溶流體，但在微觀孔隙模型中，由于孔隙體積極小(1.97μL)，無論是少量CO2溶解于水(形成不穩定的H2CO3)，還是微量水溶解于CO2，都會對兩相流運動特性產生顯著影響[24].為了消除溶解對兩相流的影響，在泵A恒壓(8.5MPa)條件下，攪拌棒勻速轉動，確保攪拌器上層CO2與下層水充分反應，以達到兩相平衡(充分飽和).

(4)溫度控制子系統由隔熱室、溫度傳感控制器和加熱燈組成.(1)、(2)和(3)三個子系統位于密閉的隔熱室中，室內由加熱燈提供熱量，室內的溫度傳感器(Omega,Type T)與室外的溫度控制器(Cole Parmer,EW-89000-10)相連，以確保隔熱室內的溫度保持在 (45±1)?C.

圖1 高壓流體--顯微鏡--微觀模型實驗裝置Fig.1 High-pressure fluids-microsco y-micromodel system

1.2 實驗流程

實驗過程包括CO2、水兩種流體的制備、微觀模型清洗、CO2與水驅排、CO2與水兩相流平衡4個階段.試驗重復3次，具體過程如下：

(1)CO2與水兩相流體制備.在本研究中，不考慮CO2與水的相互溶解，因此實驗中的CO2和水須達到相平衡狀態.為反映地層咸水濃度，采用脫離子水和高純度NaCl(Sigma-Aldrich,ACS reagent grade)配置0.01mol/L的NaCl溶液.將該溶液注入攪拌器，并將CO2由氣缸注入泵A.泵A與攪拌器連接，將流體壓力和溫度分別提至8.5MPa和45?C.攪拌棒勻速轉動24h，確保CO2與水兩者充分飽和.

(2)微觀模型清洗.采用脫離子水以500μL/min的流量清洗微觀模型，持續10min，然后放置在烘箱內在120?C條件下烘24h.

(3)CO2與水驅排.常壓下，采用泵C向微觀模型注入脫離子水，確保微觀模型達到飽水狀態，然后將壓力提升至8.5MPa.打開泵C與微觀模型管路，將制備好的水(0.01mol/L NaCl)以200μL/min流量注入微觀模型，持續5min.然后，打開泵A與微觀模型管路，CO2與水在注入口A處相接觸.此后將實驗靜置3h，以確保CO2與水達到熱平衡狀態.3h后，泵A以192μL/min的流量向微觀模型中注入超臨界CO2以驅替水，該過程為排水過程(drainage).照相機以1幀/秒的速率采集圖像.由于微觀模型孔隙體積較小 (1.97μL)，該排水過程僅持續 2s便可達到穩定狀態.排水結束后，關閉泵A與微觀模型管路，打開泵B與微觀模型管路，以同樣的流量(192μL/min)注入水，即吸濕過程(imbibition).該過程僅持續3s便可達到穩定狀態

(4)CO2與水兩相流平衡.驅排過程結束后，關閉微觀模型入口與出口閥門.根據照相機采集的圖像，5min后微觀模型內流體便可達到平衡狀態.此時，采用顯微鏡采集微觀模型中被毛細捕獲的CO2液滴/團簇的形態與分布.

本實驗裝置是基于文獻[23]改進而成的，有關實驗細節，可參考文獻[23].

2 數值模擬

2.1 控制方程

微觀模型中，CO2與水兩相驅替物理過程可由納維--斯托克斯方程描述，即流體動量守恒方程

與質量守恒方程

式中,u為流體速度,p為流體壓力，μ為流體黏滯系數，fs表征了由CO2與水界面引起的毛細力對流體動量的貢獻.fs可表達為

式中,σ為界面張力,κ為曲率,n為界面指向濕潤相(水相)的單位法向量，δΓT為Dirac delta函數，界面處為1，其他處為0.

式(1)涉及到兩相流運動，在各相內部，式(1)可分別描述各個相的流體運動.對于兩相流體界面，目前廣泛采用流體體積法(volume of flui，VOF)刻畫界面的演化[35].定義無量綱變量，即流體體積分數α，來表征兩相界面的空間位置(如圖2所示)

圖2 流體體積法示意圖Fig.2 Illustration of the volume of flui method

在流體體積方法中，體積分數的演化可由如下方程描述

式中，ur為人為構造的流速，以提升界面處的數值穩定性與魯棒性[36].

根據式(4)，式(1)中的流體黏度與密度可表達為流體體積分數的函數

式中，μw和μCO2分別為水和CO2的黏滯系數，ρw和ρCO2分別為水和CO2的密度.

考慮到界面法向量 n=?α，界面曲率κ=?·(?α/‖?α‖).將 n,κ與α的關系式與式(3)代入式(1)，并聯立式(2)和式(5)，可得到CO2水兩相流運動的控制方程

2.2 模擬方法

OpenFOAM[36]是一個得到國內外廣泛應用的計算流體動力學開源軟件包，其開源特性使研究人員可根據問題求解的需要更改底層代碼.本文在Open-FOAM平臺上采用有限體積法對式(7)進行了數值實現.微觀模型的計算網格如圖3所示，共剖分六面體單元1962924個，節點2712905個，單元的尺寸控制在6～9μm.流體參數如表1所示.

圖3 微觀模型有限體積網格Fig.3 Finite volume mesh for the micromodel

表1 超臨界CO2與水兩相流參數Table 1 The parameters for the supercritical CO2-water two-phase flui fl w

需要指出的是，本文微觀模型壁面為親水性.在8.5MPa和45?C的高溫高壓條件下，基于不同區域的顯微鏡照片，測量了10組CO2與水在微觀孔隙模型壁面上的接觸角，平均值為22?.在OpenFOAM中，接觸角對兩相流界面的影響由下式表征[35]

式中,θ為壁面接觸角，ns和nt分別為壁面的單位法向量和切向量.

根據實驗條件，數值模擬的初始條件為：在模型區域內，u=0，α=1，p=8.5MPa.邊界條件為：(1)在排水階段(drainage)，在注入口A，u=0.32m/s(根據流量折算)，?p=0，α =0；在出口，?u/?x=?u/?y= ?u/?z=0，p=8.5MPa，?α =0，(2)在吸濕階段(imbibition)，在注入口B，u=0.32m/s(根據流量折算)，?p=0，α =1；在出口，?u/?x= ?u/?y=?u/?z=0，p=8.5MPa，?α =0，(3)在平衡階段(注入口與出口閥門關閉)，在注入口B和出口，u=0m/s，p=8.5MPa，?α =0.

數值模擬采用并行計算，以確保在可接受的時間內獲得解答.基于國家天河II號超級計算機，采用96個CPU核心并行模擬CO2水排水、吸濕與平衡過程，共耗機時23040核·小時.對于3.3節毛細捕獲曲線的數值模擬，采用768個CPU核心進行數值模擬，共耗機時115200核·小時.

3 結果分析

如1.2節所述，對實驗過程重復3次，獲得3組實驗圖像.對于CO2的分布特征與CO2液滴/團簇形態，采用第1組實驗圖像與數值模擬進行對比分析；對于CO2的飽和度，采用3組實驗數據的平均值與數值模擬結果進行對比.采用ImageJ軟件[39]對實驗圖形進行后處理，提取CO2飽和度的量值.

3.1 超臨界CO2水驅排實驗與數值模擬對比

圖4給出了兩相不混溶流體驅替過程中呈現出的3種流動特征，即黏性指進、毛細指進與穩定流.這3種流動特征受控于兩個特征參數，即毛細數(Ca)和黏度比(M).毛細數反映了兩相流所受到的黏滯力與毛細力之比的量值.Ca和M的定義如下

式中,ui和μi分別為驅替流體的特征流速與黏滯系數，μ2為被驅替流體的黏滯系數.由于實驗條件不一致，Lenormand等(虛線)[20]和Zhang等(實線)[28]給出的結果有一定差異，如圖4所示.

圖4 兩種互不相容流體流動狀態相圖Fig.4 The phase diagram for displacement pattern within the two-phase immiscible flui fl w

將表1所示參數代入式(8)和式(9)可得：在排水階段(drainage，CO2驅替水)，毛細數與黏度比分別為，Ca=2.24×10-5，M=0.038；在吸濕階段(imbibition，水驅替CO2)，Ca=5.87×10-4，M=26.24.由圖4可知，在本實驗的排水階段(drainage)，兩相流驅替過程既有毛細指進也有黏性指進.

圖5給出了排水0.1s時，微觀模型中CO2分布的數值模擬結果.由于CO2注入口位于左上方，而恒定壓力出口位于右下方，因此CO2指進向右下方傾斜.由圖可知，CO2驅替水過程呈現強烈的指進現象.根據指進的寬度與長度，可初步判定：在左上方，CO2鋒面推進為毛細指進，而在左下方，CO2鋒面前移為黏性指進.圖5所示的局部放大圖給出了CO2驅替鋒面處，CO2與水界面的三維形態.對于本文研究的親水性介質而言，界面向濕潤相(水相)凸起.

圖5 排水過程(drainage)0.1 s時刻數值模擬結果Fig.5 The numerical results for the two-phase flui fl w after 0.1 s of drainage

圖6給出了排水1s時，CO2分布的實驗觀測與數值模擬對比結果.由圖可知，在排水1 s時刻，驅替相 (CO2)已到達出口處，驅替穿透.考慮到排水實驗注入口與出口的位置，在微觀模型左下方與右上方局部區域，CO2的局部飽和度較小.圖6(a)的CO2飽和度為51.1%，三組實驗的CO2飽和度平均值為(51.0±3.9)%.CO2飽和度的數值模擬結果為53.9%(圖6(b)).排水結束后，CO2飽和度的實驗結果為(58.3±2.1)%，數值模擬結果為61.3%，兩者基本一致.

圖6 排水過程(drainage)1 s時刻實驗(由照相機采集)與數值模擬對比((a)實驗結果,(b)數值模擬結果)Fig.6 Comparison between numerical results and the DSLR captured images after 1 s of drainage((a)experimental result,(b)numerical simulation)

排水過程持續2s后，水從微觀模型左下方以恒定流量注入，以驅替排水過程中注入的CO2.圖7給出了吸濕0.1s時，CO2相分布特征的數值模擬結果.由圖可知，在水相注入口附近，CO2局部飽和度趨于零，表明在此區域，水以“活塞式”方式驅替CO2.圖7中的局部放大圖表明，在吸濕過程中，部分CO2以液滴形式存在于孔隙中.

吸濕1s時，CO2分布的實驗觀測與數值模擬對比結果如圖8所示.由圖可知，在吸濕1s時刻，實驗觀測結果表明CO2以液滴/團簇形式存在于單個或多個孔隙中.CO2局部飽和度在微觀模型左端和右端較高，而在中間偏左側較低.數值結果表明,在微觀模型左側 (水相注入口附近區域)，水以 “活塞式”驅替CO2，導致CO2的局部飽和度趨于零.而在此區域，實驗觀測的CO2飽和度較大.其原因是在排水階段(drainage)，有一定體積的超臨界CO2進入左下方入口處(圖8(a)).當吸濕開始時(imbibition)，在排水階段進入的超臨界CO2被水從入口區域驅替至微觀孔隙中，從而導致左下方附近區域的CO2飽和度過大.由此可見，對于本文采用的微觀流體模型，由于流體入口段不均勻，使得微觀孔隙模型中的左邊界不能成為均勻的過流斷面，從而在一定程度上影響入口端附近流場的分布規律以及CO2或水的局部飽和度.

圖7 吸濕過程(imbibition)0.1s時刻數值模擬結果Fig.7 The numerical results for the two-phase flui fl w after 0.1s of imbibition

圖 8(a)的 CO2飽和度為 29.4%，三次重復實驗的平均值為(27.7±2.3)%，數值模擬結果為21.9%(圖8(b)).由此可見，在吸濕1s時刻，數值模擬與實驗結果有一定差異.造成偏差的原因如下：(1)如前所述，在排水階段進入入口段的CO2，在吸濕初期(imbibition)被水驅替至微觀孔隙中，導致CO2飽和度的實驗值偏大；(2)由于制造工藝，微觀模型表面有微小的(約0.5μm)起伏度(粗糙度)，而在數值模擬中，幾何模型沒有反映這種起伏度.考慮到微觀模型喉道尺寸為50μm，深度為40μm，在吸濕階段，毛細數增大，黏滯力對流體運動的影響增大.因此，起伏度對流體運動的影響在吸濕階段被放大，從而導致數值模擬與實驗結果的偏差增大；(3)在本文數值模擬中，接觸角取常數，而忽略了動態接觸角的影響.當毛細數大于10-4時，流體與孔隙壁面的前沿和后沿接觸角具有一定差異，接觸角的這種動態效應對流體界面的推進具有一定影響[4041].

圖8 吸濕過程(imbibition)1s時刻實驗(由照相機采集)與數值模擬對比((a)實驗結果,(b)數值模擬結果)Fig.8 Comparison between numerical results and the DSLR captured images after 1s of imbibition((a)experimental result,(b)numerical simulation)

3.2 CO2毛細捕獲特征的實驗與數值模擬對比

水驅替CO2持續3s后，吸濕過程結束.關閉微觀模型注入口與出口.由于此刻注入口與出口端的壓力差為零，微觀模型中的流體流速將逐漸遞減至零，流體壓力將變為恒定壓力值(8.5MPa).由照相機與顯微鏡采集的圖像可知，閥門關閉5min后，兩相流達到平衡狀態.圖9給出了平衡狀態下被捕獲的CO2液滴/團簇形態與位置.由圖可知，在本文研究的親水性介質中，被捕獲的 CO2以液滴或團簇的形式存在于孔隙中，CO2液滴或團簇呈球狀(或近似球狀)，該實驗觀測結果與以往的實驗結果一致[8,10].這些液滴或團簇位于單個孔隙(pore body)中(圖9(d),圖9(f)和圖9(g))，或位于兩個相鄰的孔隙中(圖9(e))，相鄰的孔隙由孔喉連接.數值模擬結果較好地反映了被捕獲CO2液滴/團簇的這種分布特征.CO2毛細捕獲的數值模擬結果如圖9(h)～圖9(n)所示.被捕獲的CO2液滴/團簇部分位于單個孔隙中(圖9(h)、圖9(j)和圖9(k))，部分位于兩個相鄰的孔隙中(圖9(i)).由圖9(h)～圖9(k)的三維視圖可知，對于親水性介質，CO2液滴的三維形態為凸狀體，這與基于CT掃描的實驗觀測結果一致[8].

如圖 9(b)所示，CO2殘余飽和度為 10.3%，三組實驗的平均值為 (9.2±1.1)%，數值模擬結果為7.9%(圖9(m)).由于CO2殘余飽和度與驅替流速(或毛細數)密切相關，研究表明[4243]，當毛細數大于1×10-6時，一般服從Snwr～λlg(Ca)的規律，其中λ為常數.毛細數越小，兩種流體的毛細力對兩相流運動的影響越顯著.本實驗中，在吸濕階段(imbibition)，毛細數為Ca=5.87×10-4，M ＞1，毛細力對兩相流的影響不顯著，驅替模式屬于毛細指進和穩定驅替的過渡區(圖4).

圖9 超臨界CO2毛細捕獲實驗結果(a)～(g)與數值模擬結果(h)～(n).(a)和(c)為顯微鏡低倍物鏡采集的被捕獲CO2液滴/團簇分布；(b)為照相機采集的被捕獲CO2液滴/團簇整體分布；(d)～(g)為顯微鏡高倍物鏡采集的單個孔隙中被捕獲CO2團簇分布形態.(m)為被捕獲的CO2液滴/團簇整體分布；(l)和(n)為被捕獲的CO2團簇在多個孔隙內的分布；(h)～(k)為被捕獲的CO2團簇在單個孔隙中的形態Fig.9 The observed(a)～ (g)and simulated(h)～ (n)shape and distribution of trapped CO2droplets/clusters：(a)and(c)the distribution of trapped CO2droplets/clusters at multiple pores;(b)the distribution of trapped CO2droplets/clusters at the full scale of micromodel;(d)～(g)the distribution of trapped CO2droplets/clusters at a single pore.(m)the distribution of trapped CO2droplet/cluster at the full scale of micromodel;(l)and(n)the distribution of trapped CO2droplets/clusters at multiple pores;(h)～(k)the distribution of trapped CO2droplets/clusters at a single pore

受微觀孔隙模型實驗周期的限制，在本實驗中，驅替流速設定為192μL/min，沒有采取相對較小的驅替流速.此外，微觀模型可承受最大值為10MPa的內外壓力差.在以往的微觀孔隙模型實驗研究中[24,44]，一般將微觀模型置入透明壓力盒中.微觀孔隙模型內部加壓時，壓力盒(微觀孔隙外部)也同時加壓，以確保微觀孔隙內外壓差不太大.本研究未采用壓力盒，因而實驗裝置簡單，無需染色劑便可直接觀測微觀模型中的兩相流運動[23].然而，由于沒有壓力盒抵消壓差，在內外8.5MPa的壓差條件下，若實驗持續較長時間(6h)，微觀孔隙模型存在“爆裂”的危險.由于小流量驅替實驗將顯著延長實驗的持續時間，因此為了減低實驗風險，本研究采用較大的驅替流速，以保證在相對較短的時間內獲得兩相驅替的全過程與毛細捕獲的實驗結果.在以后的研究中，將改進實驗裝置，在保證實驗安全可靠的前提下，采用更小的驅替流速，研究兩相流的運動規律.

3.3 毛細捕獲曲線實驗/模擬結果與預測模型對比

在CO2地質封存和石油開采中，孔隙介質的非濕潤相初始飽和度(Snwi)與殘余飽和度(Snwr)曲線對于工程應用尤其重要，該曲線稱為毛細捕獲曲線.Land[45]提出了單參數預測模型，描述兩者之間的關系.Land模型表達如下

式中,Snwi和Snwr分別為初始飽和度和殘余飽和度，C為Land模型參數.Snwi為排水過程結束或吸濕過程開始時，介質的非濕潤相飽和度.Land模型參數可通過介質的最大飽和度和最大的殘余飽和度計算得到

式中,Snwr,max和Snw,max分別為介質非濕潤相最大殘余飽和度和最大飽和度.

Jerauld[46]在Land模型的基礎上，提出了兩參數模型，以反映介質濕潤性的非均質性對介質殘余飽和度的影響，表達如下

式中,b為模型參數.對比式(11)與式(13)可知，當b=0時，Jerauld模型直接退化為Land模型.

Spiteri等[47]提出了能夠反映不同濕潤性對介質殘余飽和度影響的預測模型

式中,α和β為模型參數.

對于上述3種殘余飽和度預測模型，采用數值模擬與實驗結果進行對比研究.在本文實驗中，當排水過程結束時，CO2分布不再改變，因此可認為CO2達到了初始飽和度的最大值.由于本文開展的排水、吸濕實驗在5s內便分別達到穩定狀態，較難獲得初始飽和度和殘余飽和度的其他量值.為獲取本文微觀模型不同的Snwi與Snwr量值，采用數值方法模擬不同持續時間(小于2s)的排水過程，得到初始飽和度，繼而模擬吸濕過程，得到對應的殘余飽和度.

圖10給出了微觀模型初始飽和度與殘余飽和度的實驗、數值模擬與預測模型擬合結果.表2列出了基于數值模擬結果，各個模型的最優擬合參數與均方偏差.由表2可知，Jerauld模型與Land模型與數值模擬結果的均方偏差分別為0.0066和0.01，兩者相差不大，而Spiteri模型的均方偏差為0.0190.因此，Jerauld模型與Land模型與數值模擬結果吻合較好，但Jerauld模型優于Land模型，Spiteri模型的描述能力最差.

表2 各模型的最優擬合系數與均方偏差Table 2 The fittin parameters and the root-mean-square-error(RMSE)for the three models

圖10 超臨界CO2毛細捕獲曲線數值模擬結果與不同模型對比Fig.10 The comparison of the capillary trapping curve between the numerical simulation/experimental data and three di ff erent models

由圖10可知，當初始飽和度小于7%時，毛細捕獲曲線與1：1直線重合，表明當注入CO2的體積分數小于0.07時，全部CO2都被捕獲于介質中.在此階段，Jerauld模型吻合效果明顯優于Land模型與Spiteri模型.而當初始飽和度大于30%時，殘余飽和度基本保持不變.Jerauld模型與Land模型都能反映這種趨勢，而Spiteri模型曲線在初始飽和度大于45%后，開始下降.圖10也給出了初始飽和度與殘余飽和度的實驗結果及其誤差圖.對于實驗結果：Snwi=(58.3±2.1)%，Snwr=(9.2±1.1)%；對于數值模擬：Snwi=61.3%，Snwr=7.9%，表明無論是初始飽和度還是殘余飽和度，數值模擬結果均小于實驗結果.數值模擬與實驗結果的偏差主要受3.1節所述的3個方面原因的影響.Raeini等[33]針對LV60砂巖開展了CO2毛細捕獲的實驗與數值模擬對比研究，表明在初始飽和度大于0.5的區域，CO2毛細捕獲的數值模擬結果與實驗結果的絕對偏差為2%～4%.在本文對比研究中，實驗與模擬的絕對偏差為1.3%.

4 結論

本文采用高壓流體--顯微鏡--微觀模型實驗裝置，開展了超臨界CO2條件下CO2水驅排實驗，研究了孔隙尺度下CO2毛細捕獲特征，同時基于開源軟件包OpenFOAM對CO2水兩相流實驗過程進行了數值模擬，最后基于數值模擬與實驗結果對比分析了3種毛細捕獲模型的優劣.主要結論如下：

(1)超臨界CO2驅替水呈現出強烈的指進現象，這種指進既包含黏性指進，也包含毛細指進.數值模擬結果較好地反映了這種指進現象，同時也很好地刻畫了指進鋒面的三維形態.

(2)對于本文研究的親水性介質，在實驗觀測中，被捕獲 CO2液滴/團簇位于單個或者多個孔隙中.數值模擬不僅反映了這一特征，還刻畫了CO2液滴/團簇的三維空間形態.

(3)本文對比分析了3種毛細捕獲曲線預測模型的描述能力.結果表明：盡管Jerauld模型描述能力最強，但比Land模型多一個參數，且Land模型參數物理意義明確.因此，在實際工程中，建議優先采用Land模型.

本文針對CO2地質封存中毛細捕獲機制的微觀實驗和數值模擬研究具有重要意義.然而，由于驅替流速較大，本文的實驗與數值模擬結果只能反映場地CO2注入井附近的超臨界CO2毛細捕獲特征.此外，在實際GCS工程中，介質濕潤性與孔隙結構對CO2毛細捕獲具有決定性影響.今后將改進現有實驗平臺，開展極小驅替流速下的實驗研究.

1 Zhang W,Li Y,Xu T,et al.Long-term variations of CO2trapped in di ff erent mechanisms in deep saline formations：a case study of the Songliao Basin,China.International Journal of Greenhouse Gas Control,2009,3(2)：161-180

2 Metz B,Davidson O,De Coninck H,et al.IPCC special report on carbon dioxide capture and storage.Prepared by Working Group III of the Intergovernmental Panel on Climate Change.IPCC,Cambridge,United Kingdom and New York,USA：Cambridge University Press,2005

3李小春,劉延鋒,白冰等.中國深部咸水含水層CO2儲存優先區域選擇.巖石力學與工程學報,2006,25(5)：963-968(Li Xiaochun,Liu Yanfeng,Bai Bing,et al.Ranking and screening of CO2Saline aquifer storage zones in China.Chinese Journal of Rock Mechanics and Engineering,2006,25(5)：963-968(in Chinese))

4李小春,方志明,魏寧等.我國CO2捕集與封存的技術路線探討.巖土力學,2009,30(9)：2674-2678(Li Xiaochun,Fang Zhiming,Wei Ning,et al.Discussion on technical roadmap of CO2capture and storage in China.Rock and Soil Mechanics,2009,30(9)：2674-2678(in Chinese))

5 Suekane T,Nobuso T,Hirai S,et al.Geological storage of carbon dioxide by residual gas and solubility trapping.International Journal of Greenhouse Gas Control,2008,2(1)：58-64

6 Al Mansoori SK,Itsekiri E,Iglauer S,et al.Measurements of nonwetting phase trapping applied to carbon dioxide storage.International Journal of Greenhouse Gas Control,2010,4(2)：283-288

7 Tanino Y,Blunt M.Laboratory investigation of capillary trapping under mixed-wet conditions. Water Resources Research,2013,49(7)：4311-4319

8 Al-Raoush RI.Impact of wettability on pore-scale characteristics of residual nonaqueous phase liquids.Environmental Science&Technology,2009,43(14)：4796-4801

9 Iglauer S,Pentland C,Busch A.CO2wettability of seal and reservoir rocks and the implications for carbon geo-sequestration.Water Resources Research,2015,51(11)：729-774

10 Chaudhary K,Bayani Cardenas M,Wolfe WW,et al.Pore-scale trapping of supercritical CO2and the role of grain wettability and shape.Geophysical Research Letters,2013,40(15)：3878-3882

11 Tanino Y,Blunt MJ.Capillary trapping in sandstones and carbonates：Dependence on pore structure.Water Resources Research,2012,48(8)：W08525

12 Moura M,Fiorentino EA,M?l?y K,et al.Impact of sample geometry on the measurement of pressure-saturation curves：Experiments and simulations.Water Resources Research,2015,51(12)：8900-8926

13 Kimbrel EH,Herring AL,Armstrong RT,et al.Experimental characterization of nonwetting phase trapping and implications for geologic CO2sequestration.International Journal of Greenhouse Gas Control,2015,42：1-15

14 Chatzis I,Kuntamukkula M,Morrow N.E ff ect of capillary number on the microstructure of residual oil in strongly water-wet sandstones.SPE Reservoir Engineering,1988,3(3)：902-912

15 El-Maghraby RM,Blunt MJ.Residual CO2trapping in Indiana limestone.Environmental Science&Technology,2012,47(1)：227-233

16 Andrew M,Bijeljic B,Blunt MJ.Pore-scale imaging of trapped supercritical carbon dioxide in sandstones and carbonates.International Journal of Greenhouse Gas Control,2014,22：1-14

17 Niu B,Al-Menhali A,Krevor SC.The impact of reservoir conditions on the residual trapping of carbon dioxide in Berea sandstone.Water Resources Research,2015,51(4)：2009-2029

18 Zhang C,Oostrom M,Grate JW,et al.Liquid CO2displacement of water in a dual-permeability pore network micromodel.Environmental Science&Technology,2011,45(17)：7581-7588

19 Wang Y,Zhang C,Wei N,et al.Experimental study of crossover from capillary to viscous fingerin for supercritical CO2–water displacement in a homogeneous pore network.Environmental Science&Technology,2012,47(1)：212-218

20 Lenormand R,Touboul E,Zarcone C.Numerical models and experiments on immiscible displacements in porous media.Journal Of Fluid Mechanics,1988,189：165-187

21 Andrew M,Bijeljic B,Blunt MJ.Pore-scale imaging of geological carbon dioxide storage under in situ conditions.Geophysical Research Letters,2013,40(15)：3915-3918

22武愛兵,李銥,常春等.不同成分鹽水驅CO2.現代地質,2014,28(5)：1061-1067(Wu Aibing,Li Yi,Chang Chun,et al.The residual gas saturation of di ff erent components of saline floodin CO2.Geoscience,2014,28(5)：1061-1067(in Chinese))

23 Kim Y,Wan J,Kneafsey TJ,et al.Dewetting of silica surfaces upon reactions with supercritical CO2and brine：pore-scale studies in micromodels.Environmental Science&Technology,2012,46(8)：4228-4235

24 Chang C,Zhou Q,Kneafsey TJ,et al.Pore-scale supercritical CO2dissolution and mass transfer under imbibition conditions.Advancesin Water Resources,2016,92：142-158

25 Chang C,Zhou Q,Oostrom M,et al.Pore-scale supercritical CO2dissolution and mass transfer under drainage conditions.Advances in Water Resources,2017,100：14-25

26 Zhao B,MacMinn CW,Juanes R.Wettability control on multiphase fl w in patterned microfluidics Proceedings of the National Academy of Sciences,2016,113(37)：10251-10256

27 Xu W,Ok JT,Xiao F,et al.E ff ect of pore geometry and interfacial tension on water-oil displacement efficiency in oil-wet microfluidi porous media analogs.Physics of Fluids,2014,26(10)：093102

28 Zhang C,Oostrom M,Wietsma TW,et al.Influenc of viscous and capillary forces on immiscible flui displacement：Pore-scale experimental study in a water-wet micromodel demonstrating viscous and capillary fingering Energy&Fuels,2011,25(8)：3493-3505

29 Cottin C,Bodiguel H,Colin A.Influenc of wetting conditions on drainage in porous media：A microfluidi study.Physical Review E,2011,84(2)：026311

30 Raeesi B,Piri M.The e ff ects of wettability and trapping on relationships between interfacial area,capillary pressure and saturation in porous media：A pore-scale network modeling approach.Journal of Hydrology,2009,376(3)：337-352

31 Liu H,Ju Y,Wang N,et al.Lattice Boltzmann modeling of contact angle and its hysteresis in two-phase fl w with large viscosity di ff erence.Physical Review E,2015,92(3)：033306

32 Bandara U,Tartakovsky AM,Oostrom M,et al.Smoothed particle hydrodynamics pore-scale simulations of unstable immiscible fl w in porous media.Advances in Water Resources,2013,62：356-369

33 Raeini AQ,Bijeljic B,Blunt MJ.Modelling capillary trapping usingfinite-olumesimulationoftwo-phasefl wdirectlyonmicro-CT images.Advances in Water Resources,2015,83：102-110

34 Ferrari A,Jimenez-Martinez J,Borgne TL,et al.Challenges in modeling unstable two-phase fl w experiments in porous micromodels.Water Resources Research,2015,51(3)：1381-1400

35 Ferrari A,Lunati I.Direct numerical simulations of interface dynamics to link capillary pressure and total surface energy.Advances in Water Resources,2013,57：19-31

36 Greenshields C.OpenFOAM User Guide.Version 2.4.0.Open-FOAM Foundation Ltd,2015

37 Batzle M,Wang Z.Seismic properties of pore fluids Geophysics,1992,57(11)：1396-408

38 Wang S,Tokunaga TK.Capillary pressure-saturation relations for supercritical CO2and brine in limestone/dolomite sands：Implications for geologic carbon sequestration in carbonate reservoirs.Environmental Science&Technology,2015,49(13)：7208-7217

39 Abr`amo ffMD,Magalh?aes PJ,Ram SJ.Image processing with Image J.Biophotonics International,2004,11(8)：36-42

40 Roman S,Soulaine C,AlSaud MA,et al.Particle velocimetry analysis of immiscible two-phase fl w in micromodels.Advances in Water Resources,2015,95：199-211

41 Horgue P,Augier F,Duru P,et al.Experimental and numerical study of two-phase fl ws in arrays of cylinders.Chemical Engineering Science,2013,102(15)：335-345

42 Geistlinger H,Ataei-Dadavi I,Mohammadian S,et al.The impact of pore structure and surface roughness on capillary trapping for 2-D and 3-D porous media：Comparison with percolation theory.Water Resources Research,2015,51(11)：9094-9111

43 Mohammadian S,Geistlinger H,Vogel HJ.Quantificatio of gasphase trapping within the capillary fringe using computed microtomography.Vadose Zone Journal,2015,14(5)：

44 Jimenez-Mart′?nez J,Porter ML,Hyman JD,et al.Mixing in a threephase system：Enhanced production of oil-wet reservoirs by CO2injection.Geophysical Research Letter,2016,43(1)：196-205

45 Land CS.Calculation of imbibition relative permeability for twoand three-phase fl w from rock properties.Society of Petroleum Engineers Journal,1968,8(2)：149-156

46 JerauldG.Generalthree-phaserelativepermeabilitymodelforPrudhoe Bay.SPE Reservoir Engineering,1997,12(04)：255-263

47 Spiteri EJ,Juanes R,Blunt MJ,et al.A new model of trapping and relative permeability hysteresis for all wettability characteristics.SPE Journal,2008,13(3)：277-288

SUPERCRITICAL CO2WATER DISPLACEMENTS AND CO2CAPILLARY TRAPPING:MICROMODEL EXPERIMENT AND NUMERICAL SIMULATION1)

Hu Ran?,?Chen Yifeng?,2)Wan Jiamin?Zhou Chuangbing??(State Key Laboratory of Water Resources and Hydropower Engineering Science，Wuhan University，Wuhan 430072，China)?(Lawrence Berkeley National Laboratory，Energy Geosciences Division，Berkeley CA 94720，USA)

The CO2capillary trapping is an important scientifi issue in geological carbon sequestration,but few researches focus on the trapping mechanism at pore scale under supercritical CO2condition.In this study,based on the high-pressure fluids-microsco y-micromodel experimental system,we performed drainage experiment,i.e.supercritical CO2displacing water,and imbibition experiment,i.e.water displacing CO2,under the conditions of 45?C and 8.5MPa.The DSLR camera was used to capture pictures of CO2-water two-phase immiscible fl w and the microscopy was used to capture the capillary trapping behavior for the supercritical CO2at the pore scale.The computational flui dynamic method was adopted to simulate the two-phase flui fl w processes.The numerical results are generally in agree-ment with the experimental observations,and further provide three-dimensional geometries on the interface during the drainage-imbibition processes and the trapped supercritical CO2droplet/cluster.Finally,the capillary trapping curve,i.e.the relationship between the initial CO2saturation and the residual saturation,was obtained from the numerical results,and we made an assessment of the three capillary trapping models,i.e.Land’s,Jurauld’s and Spiteri’s trapping models.A comparison of the models performance indicates that Jurauld’s model behaves slightly better than Land’s model,whereas Spiteri’s model behaves poorly.However,given that Land’s model only contains one parameter of clear physical meaning,it is recommended for practical use.

geological carbon sequestration,micromodel,two-phase fl w,numerical simulation,capillary trapping

O363.2,V211.1+7

：A

10.6052/0459-1879-16-237

2016–08–26 收稿，2017–02–14 錄用,2017–02–15 網絡版發表.

1)國家自然科學基金(51409198,51579188)及中國博士后科學基金(2015T80833)資助項目.

2)陳益峰,教授,主要研究方向：巖土介質多場耦合與滲流控制理論.E-mail：csyfchen@whu.edu.cn

胡冉,陳益峰,萬嘉敏,周創兵.超臨界CO2--水兩相流與CO2毛細捕獲：微觀孔隙模型實驗與數值模擬研究.力學學報,2017,

49(3)：638-648

Hu Ran,Chen Yifeng,Wan Jiamin,Zhou Chuangbing.Supercritical CO2water displacements and CO2capillary trapping：micromodel experiment and numerical simulation.Chinese Journal of Theoretical and Applied Mechanics,2017,49(3)：638-648