高溫高壓氣體在低含水率介質(zhì)中遷移的數(shù)值模擬

， ，

(西北核技術(shù)研究所，西安 710024)

1 引 言

探索高溫高壓氣體在多孔介質(zhì)中的運(yùn)動(dòng)規(guī)律具有廣泛的應(yīng)用背景，可用于地下爆炸實(shí)驗(yàn)中的氣體輸運(yùn)[1]和溫室氣體減排中的二氧化碳地質(zhì)封存[2]。實(shí)際的地介質(zhì)中可能含有地下水，因此氣體在地介質(zhì)中的運(yùn)動(dòng)本質(zhì)上是一個(gè)氣液兩相滲流過程，屬于滲流力學(xué)的研究前沿問題。

Nitao[3]開發(fā)的NUFT程序可用來模擬污染物在地介質(zhì)中的運(yùn)移，如Carrigan等[4]利用USNT模塊對(duì)地下核爆炸后氙的放射性同位素的遷移行為進(jìn)行了模擬；Pruess等[5]開發(fā)的TOUGH2系列程序是目前研究多相多組分非等溫滲流的主流軟件之一，文獻(xiàn)[6]利用該程序?qū)?nèi)華達(dá)試驗(yàn)場(chǎng)的尤卡山進(jìn)行了試驗(yàn)場(chǎng)尺度的研究；Flemisch等[7]開發(fā)的DuMux是一個(gè)模擬多孔介質(zhì)內(nèi)流體運(yùn)移過程的開源軟件，涵蓋眾多模型，如三相三組分非等溫滲流模型[8]等。

姜廣輝[9]對(duì)裂隙巖體中氣液兩相的滲流特性進(jìn)行了實(shí)驗(yàn)測(cè)量和理論分析，研究了裂隙中氣液兩相的滲流機(jī)制；蔣蘭蘭[2]設(shè)計(jì)并搭建了核磁共振成像和X射線CT成像實(shí)驗(yàn)平臺(tái)，進(jìn)行了CO2排泄和吸吮的初步實(shí)驗(yàn)；鞠斌山等[10]建立了多組分非等溫CO2混相驅(qū)動(dòng)模型，考慮了對(duì)流、彌散、吸附以及儲(chǔ)層破壞現(xiàn)象；何增等[11]建立了氣液兩相多組分非等溫滲流的數(shù)學(xué)模型，并與氣體滲流模型進(jìn)行了比較。

對(duì)比可以看出，國(guó)外在該領(lǐng)域已形成眾多成熟的軟件和模擬器，而國(guó)內(nèi)在數(shù)學(xué)模型和數(shù)值模擬方面的進(jìn)展尚不充分。本文將針對(duì)文獻(xiàn)[11]的數(shù)學(xué)模型發(fā)展適宜的數(shù)值離散方法，并對(duì)高溫高壓氣體在低含水率地介質(zhì)中的輸運(yùn)過程進(jìn)行數(shù)值模擬，主要研究多孔介質(zhì)的孔隙度和滲透率對(duì)氣體遷移的影響。

2 數(shù)學(xué)模型

氣液兩相多組分非等溫滲流的數(shù)學(xué)模型，參考文獻(xiàn)[11]。兩相為液相(l)和氣相(g)，用下標(biāo)α表示；組分為水(w)、空氣(a)和示蹤物(t)，用上標(biāo)κ表示。考慮的物理機(jī)制有對(duì)流、彌散、熱傳導(dǎo)、熱輻射、氣體壓縮性以及水在氣液兩相間的平衡。

2.1.1 連續(xù)性方程

以組分為研究對(duì)象建立連續(xù)性方程。水的連續(xù)性方程為

(1)

空氣和示蹤物的連續(xù)性方程為

(κ=a,t) (2)

2.1.2 運(yùn)動(dòng)方程

多孔介質(zhì)中流體的運(yùn)動(dòng)方程為達(dá)西定律[12]，即

(3)

2.1.3 能量守恒方程

考慮對(duì)流傳熱、彌散傳熱以及熱傳導(dǎo)，建立能量守恒方程為

(4)

2.2.1 氣體的狀態(tài)方程

對(duì)氣體采用理想氣體狀態(tài)方程。氣體的總壓等于各組分氣體的分壓之和，即

(5)

(6)

本文使用國(guó)際水及水蒸氣性質(zhì)協(xié)會(huì)推薦的水的飽和蒸汽壓公式[16]。該公式的誤差小于0.1%，但形式較為復(fù)雜，指數(shù)項(xiàng)包含6項(xiàng)，其中3項(xiàng)為整數(shù)次冪，3項(xiàng)為半整數(shù)次冪。

2.2.2 相對(duì)滲透率

相對(duì)滲透率kr α一般簡(jiǎn)化為有效飽和度Se的函數(shù)，如常用的相對(duì)滲透率表達(dá)式[13,17]為

kr l=Sne,kr g=(1-Se)n

(7)

式中n=3可用于模擬考慮毛細(xì)壓力的情形[17]。有效飽和度Se的定義為

(8)

式中Sα r為α相的殘余飽和度。

2.2.3 毛細(xì)壓強(qiáng)

多孔介質(zhì)中氣液兩相的壓強(qiáng)在分界面處并不連續(xù)，其差值稱為毛細(xì)壓強(qiáng)pc，即

pc=pg-pl

(9)

毛細(xì)壓強(qiáng)一般也簡(jiǎn)化為有效飽和度的函數(shù)，本文使用Leverett[18]給出的關(guān)系式

2.120(1-Se)2+1.263(1-Se)3]

(10)

式中σ為氣液兩相間的表面張力系數(shù)。

3 數(shù)值方法

氣液兩相非等溫滲流的數(shù)學(xué)模型呈現(xiàn)出高度非線性的特征，主要表現(xiàn)在非達(dá)西流動(dòng)、水蒸氣和液態(tài)水的平衡以及復(fù)雜的本構(gòu)關(guān)系(如水的飽和蒸汽壓公式)等方面，因此合理的離散方案對(duì)數(shù)值求解至關(guān)重要。

在多相非等溫滲流中，可能會(huì)發(fā)生相的出現(xiàn)與消失。處于不同的相時(shí)，完整描述該系統(tǒng)所需要的獨(dú)立變量也不同，這組變量稱作首要變量，其他量稱作次要變量。根據(jù)吉布斯相率，一個(gè)多相多組分系統(tǒng)的自由度數(shù)F為[8]

F=C-P+2+(P-1)=C+1

(11)

表1總結(jié)了三種相態(tài)的首要變量和相出現(xiàn)的條件，圖1給出了完整的相態(tài)轉(zhuǎn)變準(zhǔn)則?？梢钥闯觯琍S1和PS2、PS1和PS3之間可直接轉(zhuǎn)變，而PS2和PS3的轉(zhuǎn)變則需要經(jīng)歷PS1的中間相態(tài)，表現(xiàn)為需要更多的迭代求解。此外，水的連續(xù)性方程在PS3時(shí)為橢圓型方程，而在其他相態(tài)時(shí)為拋物型方程，為設(shè)計(jì)統(tǒng)一的算法帶來額外的復(fù)雜性。本文采取最低飽和度法簡(jiǎn)化問題[19],即假設(shè)氣相總是存在，并設(shè)定最低飽和度(如Smin=10-6)避免氣液共存時(shí)相的消失。也就是說，利用相態(tài)PS1在Sg=Smin的極限情形來近似PS3，從而減少相態(tài)轉(zhuǎn)變的判斷。此外，通過分析可知，無論處于相態(tài)PS1還是PS2，未知數(shù)個(gè)數(shù)與方程數(shù)目都相等，即該問題是封閉的。

空間離散主要解決通量的計(jì)算問題。觀察控制方程可以發(fā)現(xiàn),運(yùn)動(dòng)方程與時(shí)間無關(guān)，使用中心差分離散；連續(xù)性方程和能量方程均具有守恒方程的形式，使用有限體積法離散，其中，對(duì)對(duì)流項(xiàng)采用迎風(fēng)格式，對(duì)擴(kuò)散項(xiàng)采用中心差分格式。計(jì)算通量還需要處理參數(shù)在界面的等效值[6,11]。

設(shè)Δtn=tn +1-tn，并利用fn+1,m表示物理量f在時(shí)間步長(zhǎng)Δtn后的第m步迭代值，迭代初值

表1 相態(tài)類型、首要變量以及相出現(xiàn)的條件

Tab.1 Phase states,corresponding primaryvariables and conditions for phase appearance

PhasestatePresentphasesPrimary variablesAppearance of phaseLiquidGasPS1l,gSl,Xwg,pg,T——PS2gXwg,Xag,pg,TTpwsat—PS3lpg,T—pwsat>pg

圖1 相態(tài)轉(zhuǎn)變準(zhǔn)則

Fig.1 Criteria of phase transition

fn +1,0=fn。記控制方程的通量為Fnβ，下標(biāo)β可取w，a，t和e，分別表示水、空氣以及示蹤氣體的連續(xù)性方程和能量守恒方程。計(jì)算步驟如下。

(1)m=0。

(4) 求解水的連續(xù)性方程

(12)

(5) 將Sn+1,m+1g代入空氣和示蹤氣體的連續(xù)性方程,

(κ=a,t)(13)

(6) 更新氣相狀態(tài)

(14)

(κ=a,t) (15)

(7) 將飽和度和氣相狀態(tài)代入能量守恒方程

(16)

求解得到Tn +1,m +1，其中Et為單元的總能量，在已知密度和飽和度時(shí)僅依賴于溫度，如在比熱容為常數(shù)時(shí),Et與溫度呈線性關(guān)系。

(8) 將溫度Tn +1,m +1代入水蒸氣分壓方程式(6)，更新水蒸氣密度為

(17)

(9) 修正氣相狀態(tài)

(18)

(κ=w,a,t)(19)

(10) 計(jì)算流體壓強(qiáng)

(20)

(21)

(11)m=m+1，重復(fù)步驟(4～10)，直至達(dá)到收斂標(biāo)準(zhǔn)，最后將收斂的數(shù)值解更新為tn +1時(shí)刻的數(shù)值， 即fn+1=fn+1,m+1。

4 計(jì)算結(jié)果與討論

采用第3節(jié)的模型和算法，編制了二維柱對(duì)稱坐標(biāo)系下求解氣液兩相非等溫滲流的計(jì)算程序。文獻(xiàn)[11]對(duì)數(shù)學(xué)模型的有效性進(jìn)行了檢驗(yàn)。本節(jié)研究高溫高壓氣體在低含水率地介質(zhì)中的輸運(yùn)問題。

求解問題如下，研究域是半徑為50 m的球體，中心處是半徑為5 m的高溫高壓氣體區(qū)，稱作空腔。在初始時(shí)刻，空腔內(nèi)充滿了水蒸氣、空氣和一氧化碳，質(zhì)量分別為802 kg，1356 kg和1727 kg，溫度為2000 K，壓強(qiáng)為4.86 MPa；研究域的其余位置與外界環(huán)境達(dá)到平衡，外界環(huán)境不含一氧化碳，氣體壓強(qiáng)為1 atm，溫度為20 ℃。主要計(jì)算參數(shù)列入表2。

圖2 Picard迭代法在一個(gè)循環(huán)內(nèi)部的完整流程

Fig.2 Complete flow chart of the Picard iteration method

表2 主要計(jì)算參數(shù)

Tab.2 Main parameters used in the calculation

K10-13 m2μl1.00×10-3 Pa·sSgr0?0.01μκg1.80×10-5 Pa·sε0.9ρl998.2 kg/m3Slr0.1

(22)

結(jié)果表明，氣體在低含水率介質(zhì)和等效孔隙度的干燥介質(zhì)內(nèi)的遷移行為基本一致。

圖4展示了孔隙度為0.010且滲透率為10-13m2條件下，空腔內(nèi)溫度和液態(tài)水質(zhì)量的時(shí)間歷程。由于空腔內(nèi)的初始溫度高于水的臨界溫度647.1 K，因此水只能以水蒸氣的形式存在。隨著空腔內(nèi)的溫度和壓強(qiáng)降低，在842 s時(shí)，空腔溫度下降到115.3 ℃，水蒸氣分壓為1.713×105Pa，等于該溫度下的飽和蒸汽壓，液態(tài)水開始出現(xiàn)。需要注意的是，該算例中水的相態(tài)轉(zhuǎn)變溫度大于100.0 ℃，這是由于水蒸氣的分壓大于一個(gè)標(biāo)準(zhǔn)大氣壓。在842 s時(shí)刻，空腔溫度下降的趨勢(shì)減緩，這是因?yàn)樗魵庖夯瘯?huì)釋放大量的潛熱。隨著溫度降低，液態(tài)水逐漸增多，但同時(shí)也在壓差作用下向周圍多孔介質(zhì)中滲透。在約104s后，空腔溫度基本保持不變，液態(tài)水的質(zhì)量由于向外滲出而逐漸減少。

圖3t=3600 s時(shí)刻示蹤氣體質(zhì)量分?jǐn)?shù)的分布(不同含水率)

Fig.3 Distribution of mass fractions of gas tracer att=3600 s，with different water contents

圖4 空腔內(nèi)溫度和液態(tài)水質(zhì)量的時(shí)間歷程

Fig.4 Time history of temperature and mass of liquid water in the cavity

4.3.1 孔隙度

圖5展示了不同孔隙度下氣相壓強(qiáng)和示蹤氣體份額在600 s時(shí)刻的分布，滲透率固定為10-13m2?？梢钥闯?，孔隙度越小，滲流發(fā)展得越快。這是因?yàn)榭紫抖仍叫〉慕橘|(zhì)內(nèi)氣壓越大，向外推進(jìn)的速度也越大，最終造成上述現(xiàn)象。也可以理解為孔隙度對(duì)空腔內(nèi)氣體演化的影響很小，如圖6所示，為了容納大致等量的氣體，孔隙度小的介質(zhì)所需的體積更多，即氣體運(yùn)動(dòng)的距離越遠(yuǎn)。

圖5t=600 s時(shí)刻氣相壓強(qiáng)和示蹤氣體份額的空間分布(不同孔隙度)

Fig.5 Distribution of gas pressure and mass fraction of gas tracer att=600 s，with different porosities

圖6 氣相壓強(qiáng)和示蹤氣體份額的時(shí)間歷程(不同孔隙度)

Fig.6 Time history of gas pressure and mass fraction of gas tracer,with different porosities

(23)

表3t=600 s時(shí)刻示蹤氣體遷移距離的模擬值與估算值的比較(不同孔隙度)

Tab.3 Comparisons between the simulated and evaluated values of the transport distance of gas tracer att=600 s，with different porosities

?0.0050.0100.0150.020Rc/m29.84(29.72)23.6320.63(20.68)18.70(18.81)

圖6展示了空腔內(nèi)和R=5.25 m處氣相壓強(qiáng)和示蹤氣體份額的時(shí)間歷程?？紫抖葘?duì)空腔內(nèi)氣體演化的影響很小，但在空腔外則不同。由于孔隙度小的介質(zhì)容納氣體的空間小，因此氣相壓強(qiáng)更大。來自空腔的氣體將與孔隙內(nèi)的氣體混合，對(duì)于孔隙度小的介質(zhì)，混合過程將更為迅速，即示蹤氣體份額上升更快。

4.3.2 滲透率

圖7展示了不同滲透率下，氣相壓強(qiáng)和示蹤氣體份額在600 s時(shí)刻的分布，孔隙度固定為0.010。根據(jù)達(dá)西定律，滲透率越高，氣體速度越快，流動(dòng)影響的范圍也越大，如K=10-12m2對(duì)應(yīng)的氣相壓強(qiáng)已接近平衡。

圖7t=600 s時(shí)刻氣相壓強(qiáng)和示蹤氣體份額的空間分布(不同滲透率)

Fig.7 Distribution of gas pressure and mass fraction of gas tracer att=600 s，with different permeabilities

圖8 氣相壓強(qiáng)和示蹤氣體份額的時(shí)間歷程(不同滲透率)

Fig.8 Time history of gas pressure and mass fraction of gas tracer,with different permeabilities

表4列出了示蹤氣體在兩個(gè)時(shí)刻的遷移距離。通過分析數(shù)據(jù)，發(fā)現(xiàn)滲透率K與遷移距離Rc的5次方近似成正比，即

(24)

表4括號(hào)內(nèi)給出了由式(24)估算的Rc，以K=10-13m2的數(shù)據(jù)為基準(zhǔn)。除了在邊界處差異較大外，估算值與計(jì)算值均吻合較好，誤差在4%以內(nèi)。

圖8展示了空腔內(nèi)和R=5.25 m處氣相壓強(qiáng)和示蹤氣體份額的時(shí)間歷程。與孔隙度不同，滲透率對(duì)空腔和介質(zhì)中的物理量均有明顯的影響。滲透率越高，空腔內(nèi)氣相壓強(qiáng)下降越快，示蹤氣體份額上升越快，介質(zhì)中的氣相壓強(qiáng)和示蹤氣體份額上升越快，這顯然是由滲流速度更快引起的。

表4 不同時(shí)刻示蹤氣體遷移距離的模擬值與估算值的比較(不同滲透率)

Tab.4 Comparisons between the simulated and evaluated values of the transport distance of gas tracer at different times,with different porosities

K/m21.0×10-143.0×10-141.0×10-13t=600 s14.42(14.91)18.23(18.57)23.63t=3600 s21.39(21.84)27.41(27.20)34.61K/m23.0×10-131.0×10-12t=600 s29.54(29.44)36.67(37.45)t=3600 s41.58(43.11)45.04(54.85)

5 結(jié) 論

本文簡(jiǎn)單介紹了兩相多組分非等溫滲流的數(shù)學(xué)模型，并建立了適宜的離散方法對(duì)其進(jìn)行數(shù)值求解。通過分析控制方程性質(zhì)和相態(tài)轉(zhuǎn)換準(zhǔn)則，采用最低飽和度法降低算法的復(fù)雜性。使用有限體積格式進(jìn)行空間離散，并針對(duì)控制方程組高度非線性的特點(diǎn)，設(shè)計(jì)了一套專門的Picard迭代法進(jìn)行時(shí)間離散，解決了氣液兩相滲流的強(qiáng)耦合問題。

利用該算法對(duì)高溫高壓氣體在低含水率介質(zhì)中的遷移行為進(jìn)行了模擬。證明了利用等效孔隙度的干燥介質(zhì)近似低含水率介質(zhì)的有效性，并分析了空腔內(nèi)的氣液相變轉(zhuǎn)換過程。在此基礎(chǔ)上，初步得到了孔隙度和滲透率對(duì)氣體輸運(yùn)的影響規(guī)律，著重研究了示蹤氣體遷移距離。計(jì)算結(jié)果表明，在相同的滲透率下，遷移距離隨孔隙度減小而增加；在相同的孔隙度下，遷移距離隨滲透率增加而增加。本文的模型和方法也可用于高含水率介質(zhì)中氣體遷移行為的研究。

:

[1] 喬登江.地下核爆炸現(xiàn)象學(xué)概論[M].北京：國(guó)防工業(yè)出版社,2002.(QIAO Deng-jiang.IntroductiontoUndergroundNuclearExplosionPhenomenology[M].Beijing：National Defense Industry Press,2002.(in Chinese))

[2] 蔣蘭蘭.CO2地質(zhì)封存多孔介質(zhì)內(nèi)氣液兩相滲流特性研究 [D].大連理工大學(xué),2014.(JIANG Lan-lan.Study on Gas/Liquid Flow Characteristic in Porous Media During CO2Geological Storage [D].Dalian University of Technology,2014.(in Chinese))

[3] Nitao J J.Reference Manual for the NUFT Flow and Transport Code,Version 3.0 [R].Report UCRL-MA-130651-REV-1,Lawrence Livermore National Laboratory,Livenmore,CA，2000.

[4] Carrigan C R,Sun Y W,Hunter S L,et al.Delayed signatures of underground nuclear explosions [J].ScientificReports,2016,6:23032.

[5] Pruess K,Oldenburg C M,Moridis G J.TOUGH2 User’s Guide,Version 2 [R].Report LBNL-43134,Lawrence Berkeley National Laboratory,Berkeley,CA，1999.

[6] Wu Y S,Pruess K.Numerical simulation of non-isothermal multiphase tracer transport in heterogeneous fractured porous media [J].AdvancesinWaterResources,2000,23(7)：699-723.

[7] Flemisch B,Darcis M,Erbertseder K,et al.DuMux：DUNE for multi-{phase,component,scale,phy-sics,...} flow and transport in porous media [J].AdvancesinWaterResources,2011,34(9)：1102-1112.

[8] Class H,Helmig R,Bastian P.Numerical simulation of non-isothermal multiphase multicomponent processes in porous media.:1.An efficient solution technique [J].AdvancesinWaterResources,2002,25(5):533-550.

[9] 姜廣輝.裂隙氣液兩相滲流特性試驗(yàn)研究 [D].中國(guó)礦業(yè)大學(xué),2014.(JIANG Guang-hui.Experimental Research on Gas-Water Two Phase Flow Characteristics in Rock Fractures[D].China University of Mining and Technology,2014.(in Chinese))

[10] Ju B S,Wu Y S,Qin J S,et al.Modeling CO2miscible flooding for enhanced oil recovery [J].PetroleumScience,2012,9(2):192-198.

[11] 何 增,田 宙,王鐵良.氣液兩相多組分非等溫滲流的數(shù)學(xué)模型 [J].現(xiàn)代應(yīng)用物理,2016,7(2):021001.(HE Zeng,TIAN Zhou,WANG Tie -liang.Mathematical model for non-isothermal liquid-gas two -phase multi-component seepage flow in porous media [J].ModernAppliedPhysics,2016,7(2):021001.(in Chinese))

[12] 孔祥言.高等滲流力學(xué)(第二版)[M].合肥：中國(guó)科學(xué)技術(shù)大學(xué)出版社,2010.(KONG Xiang-yan.AdvancedMechanicsofFluidsinPorousMedia(SecondEdition)[M].Hefei：University of Science and Tech-nology Press,2010.(in Chinese))

[13] Jambhekar V A.Forchheimer Porous-Media Flow Models：Numerical Investigation and Comparison with Experimental Data [D].Universit?t Stuttgart,2011.

[14] 王鐵良,曹 淵,張建鑫.地下爆炸空腔壓力和溫度歷程數(shù)值模擬 [J].計(jì)算物理,2011,28(5)：713-718.(WANG Tie -liang,CAO Yuan,ZHANG Jian-xin.Numerical simulation of cavity pressure and temperature in underground detonation [J].ChineseJournalofComputationalPhysics,2011,28(5)：713-718.(in Chinese))

[15] 楊世銘,陶文銓.傳熱學(xué)(第三版)[M].北京：高等教育出版社,1998.(YANG Shi-ming,TAO Wen-quan.HeatTransfer(ThirdEdition)[M].Beijing：Higher Education Press,1998.(in Chinese))

[16] Wagner W,Pruв A.The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use [J].JournalofPhysicalandChemicalReferenceData,2002,31(2)：387-535.

[17] Chueh C C,Bangerth W,Djilali N.Integrated adaptive numerical methods for transient two-phase flow in heterogeneous porous media [A].Seventh International Conference on Computational Fluid Dynamics[C].Hawaii,2012.

[18] Leverett M C.Capillary behavior in porous solids [J].TransactionsoftheAIME,1941,142(1)：152-169.

[19] Abriola L M,Pinder G F.A multiphase approach to the modeling of porous media contamination by organic compounds：2.Numerical simulation [J].WaterResourcesResearch,1985,21(1)：19-26.