氣液兩相多組分非等溫滲流的數(shù)學(xué)模型

何　增，田　宙，王鐵良

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

?

氣液兩相多組分非等溫滲流的數(shù)學(xué)模型

何增，田宙，王鐵良

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

摘要：通過分析氣液兩相在多孔介質(zhì)中的相互作用過程，建立了一套描述兩相多組分非等溫滲流的數(shù)學(xué)模型。所建立的控制方程中，既包括對(duì)流、彌散和源匯作用，也包括非達(dá)西流動(dòng)、水氣轉(zhuǎn)變以及輻射傳熱等高溫高壓環(huán)境下的特殊物理過程。提供了一套完整的使控制方程封閉的本構(gòu)關(guān)系，并推導(dǎo)了比能和比焓的計(jì)算公式。利用本文模型對(duì)地下爆炸氣體的遷移過程進(jìn)行了數(shù)值模擬。結(jié)果表明：與氣相滲流模型相比，該模型可以更合理地描述氣體在地質(zhì)介質(zhì)中的輸運(yùn)行為。

關(guān)鍵詞：兩相滲流；非等溫流動(dòng)；數(shù)學(xué)模型；本構(gòu)關(guān)系；氣體輸運(yùn)

兩相非等溫滲流技術(shù)在環(huán)境保護(hù)、生物科學(xué)、能源工業(yè)等領(lǐng)域廣泛應(yīng)用[1-2]。兩相滲流的研究以滲流力學(xué)為基礎(chǔ)，與巖石力學(xué)、表面物理、傳熱學(xué)等學(xué)科交叉，一直是滲流領(lǐng)域的研究熱點(diǎn)。國(guó)內(nèi)外學(xué)者開展了大量的研究工作，并建立了許多理論模型。例如，美國(guó)勞倫斯·伯克利國(guó)家實(shí)驗(yàn)室Pruess等開發(fā)的TOUGH2系列程序[3-4]可用于模擬孔隙-裂隙介質(zhì)中多相多組分的非等溫流動(dòng)，主要應(yīng)用于地?zé)岽鎯?chǔ)、核廢料處置、環(huán)境評(píng)估及修復(fù)等領(lǐng)域，如文獻(xiàn)[5]利用該程序?qū)?nèi)華達(dá)試驗(yàn)場(chǎng)尤卡山中的多相示蹤劑輸運(yùn)問題進(jìn)行了試驗(yàn)場(chǎng)尺度的研究，Class等建立了低壓下適用于10～300℃的三相三組分非等溫滲流模型[6]，并進(jìn)行了驗(yàn)證與確認(rèn)。此外，常見的程序還有美國(guó)勞倫斯·利弗莫爾國(guó)家實(shí)驗(yàn)室Nitao等研制的NUFT程序[7-8]、斯圖加特大學(xué)Helmig等發(fā)布的MUFTE-UG程序[9]等。國(guó)內(nèi)劉昌軍推導(dǎo)了不可壓縮水氣兩相滲流的微分方程，并對(duì)降水入滲下飽和-非飽和水氣兩相滲流作了深入探討[10]。劉曉麗等基于巖體滲流水力學(xué)和多相滲流力學(xué)理論，建立了不可壓縮水氣兩相滲流與雙重介質(zhì)變形的流固耦合數(shù)學(xué)模型[11]。黃朝琴等建立了不可壓縮兩相滲流的離散裂隙型模型[12]，研究了裂隙性油藏的注水開發(fā)過程。

對(duì)比國(guó)內(nèi)外相關(guān)研究可以看出，國(guó)內(nèi)研究所使用的模型較為簡(jiǎn)單，未考慮氣體的壓縮性和非等溫效應(yīng)。本文綜合考慮氣液在多孔介質(zhì)中的相互作用過程，建立了描述氣液兩相多組分非等溫滲流的數(shù)學(xué)模型，包括控制方程和本構(gòu)關(guān)系兩部分，并以此為基礎(chǔ)研究了地下爆炸后氣體的遷移行為。

1氣液兩相多組分非等溫滲流的數(shù)學(xué)模型

兩相為液相(l)和氣相(g)，用α表示相態(tài)；組分為水(w)、空氣(a)和示蹤物(t，地下爆炸的示蹤物一般為一氧化碳)，用κ表示組分。本文將空氣簡(jiǎn)化為單一組分的氣體。控制方程包括質(zhì)量守恒方程、運(yùn)動(dòng)方程和能量守恒方程，為了使方程封閉，還應(yīng)該補(bǔ)充相應(yīng)的本構(gòu)關(guān)系。本文模型可用來描述任意多組分氣體，這里以三組分為例進(jìn)行說明。

模型的適用條件為：1) 空氣和示蹤物均不溶于液態(tài)水；2) 在氣液兩相共存時(shí)，水蒸氣為飽和蒸汽；3) 氣體組分均滿足理想氣體狀態(tài)方程；4) 固體骨架和液態(tài)水不可壓縮；5) 不考慮毛細(xì)壓力滯后效應(yīng)。

1.1控制方程

以組分為研究對(duì)象可以很方便地得到質(zhì)量守恒方程，而不必考慮組分在相間轉(zhuǎn)化的細(xì)節(jié)。相作為一個(gè)整體是一起運(yùn)動(dòng)的，故運(yùn)動(dòng)方程需要分相建立。采取局部熱平衡假設(shè)構(gòu)建能量方程：在一個(gè)代表性單元內(nèi)，流體和骨架具有相同的溫度。質(zhì)量和能量轉(zhuǎn)移的途徑在液相中有對(duì)流遷移與源匯作用，在氣相中還包含水動(dòng)力彌散效應(yīng)。

1.1.1質(zhì)量守恒方程(連續(xù)性方程)

水在氣相和液相共存，而空氣和示蹤物只存在于氣相。因此，水的連續(xù)性方程為

(1)

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

(2)

其中，φ為多孔介質(zhì)孔隙度；Sl、Sg分別為液相、氣相的飽和度；ρl、ρg為液相、氣相的密度；Xg(κ)為組分κ在氣相中的質(zhì)量分?jǐn)?shù)；vl、vg分別為液相、氣相的滲流速度；Dg(κ)為組分κ在氣相中的水動(dòng)力彌散系數(shù)；q(κ)為組分κ的質(zhì)量源匯項(xiàng)。

對(duì)彌散通量，存在如下制約關(guān)系：

(3)

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

根據(jù)達(dá)西定律[1]，多孔介質(zhì)中液相和氣相的運(yùn)動(dòng)方程為

(4)

其中，K為固有滲透率張量；kr,α為α相的相對(duì)滲透率；μα為α相的動(dòng)力黏度；pα為α相的壓強(qiáng)；g為重力加速度矢量。

在流速較大時(shí)，許多學(xué)者對(duì)達(dá)西定律提出了修正方案[2]。對(duì)于多相非達(dá)西流動(dòng)，F(xiàn)orchheimer方程的形式為[14]

(5)

其中，η為固有通過率；ηr,α為α相的相對(duì)通過率；β為Forchheimer系數(shù)，其與固有通過率的關(guān)系為β=1/η。

式(5)右端的第3項(xiàng)表示慣性力的貢獻(xiàn)。為了表征流動(dòng)偏離達(dá)西定律的程度，定義第3項(xiàng)與第2項(xiàng)之比為Forchheimer數(shù)[15]

(6)

在應(yīng)用時(shí)，首先由達(dá)西定律求出速度vα，再由式(6)計(jì)算Foα。將Foα與其判斷標(biāo)準(zhǔn)Foc進(jìn)行比較，如果Foα>Foc，則用式(5)重新求vα。

1.1.3能量守恒方程

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

(7)

其中，uα為α相的比能；hα為α相的比焓；hg(κ)為處于氣相的組分κ的比焓；λp為多孔介質(zhì)的表觀熱導(dǎo)率；qh為能量源匯項(xiàng)；ρs和us分別表示固體骨架的密度和比能，下標(biāo)s表示固體骨架。

研究表明，在溫度變化范圍較大時(shí)，忽略輻射傳熱將引起較大誤差[16]。輻射傳熱等于輻射熱流密度φ與面積的乘積[17]，φ的表達(dá)式為

(8)

其中，T為溫度；為介質(zhì)黑度，大部分非金屬材料的黑度在0.85～0.95之間，本文取=0.9；C0為黑體輻射系數(shù)。

1.2本構(gòu)關(guān)系

為使上述控制方程封閉，還需要提供相應(yīng)的本構(gòu)關(guān)系。本節(jié)將給出具體的本構(gòu)關(guān)系，重點(diǎn)對(duì)比能和比焓的計(jì)算公式進(jìn)行推導(dǎo)。需要說明的是，部分參數(shù)(如相對(duì)滲透率)的計(jì)算有多種經(jīng)驗(yàn)公式，本文只給出一種。

1.2.1約束性條件

氣相和液相的飽和度滿足

(9)

氣相中各組分的質(zhì)量分?jǐn)?shù)滿足

(10)

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

采用理想氣體狀態(tài)方程描述氣體。氣體的總壓力等于各組分氣體的分壓之和：

(11)

其中，R(κ)為組分κ的氣體常數(shù)。在氣液兩相共存時(shí)，假設(shè)水蒸氣為飽和蒸汽，即水蒸氣的分壓等于相同溫度下的飽和蒸汽壓。

(12)

本文使用國(guó)際水及水蒸氣性質(zhì)協(xié)會(huì)(TheInternationalAssociationforPropertiesofWaterandSteam,IAPWS)推薦的水的飽和蒸汽壓公式[18]。

1.2.3水動(dòng)力彌散系數(shù)

水動(dòng)力彌散系數(shù)包含機(jī)械彌散和分子擴(kuò)散兩項(xiàng)，其分量計(jì)算公式為[5,19]

(13)

其中，αT、αL分別為介質(zhì)的橫向、縱向彌散度；τ為介質(zhì)的彎曲因子；Dd(κ)為組分κ在氣相中的分子擴(kuò)散系數(shù)。

利用Wilke方法[13]，通過雙組分?jǐn)U散系數(shù)得到多組分?jǐn)U散系數(shù)Dd(w)和Dd(t)為

(14)

其中，Dd(w-a)、Dd(w-t)和Dd(a-t)分別為水蒸氣-空氣、水蒸氣-示蹤物及空氣-示蹤物的擴(kuò)散系數(shù)；xg(w)、xg(a)、xg(t)分別為組分w、a、t的摩爾分?jǐn)?shù)。在缺乏觀測(cè)資料時(shí)，雙組分?jǐn)U散系數(shù)一般采取Wilke-Lee方法進(jìn)行估計(jì)[20]?？諝獾膹浬⑾禂?shù)Dg(a)可根據(jù)限制條件式(3)得到。

Millington和Quirk提出的彎曲因子的表達(dá)式為[21]

(15)

1.2.4相對(duì)滲透率與相對(duì)通過率

相對(duì)滲透率kr,α描述的是多相滲流中有效流動(dòng)面積減小產(chǎn)生的影響，一般將其簡(jiǎn)化為有效飽和度Se的函數(shù)，如常用的相對(duì)滲透率表達(dá)式為[14,22]

(16)

本文考慮毛細(xì)壓力的情形[22]，取n=3。有效飽和度Se的定義為

(17)

其中，Sres,l和Sres,g分別為l相和g相流體的殘余飽和度。類似地，相對(duì)通過率的常用表達(dá)式為[14]

(18)

文獻(xiàn)[14]推薦m的取值可為3，5，6，本文取m=3。

1.2.5動(dòng)力黏度

混合氣體的黏度與其成分有關(guān)。本文采用Wilke方法[20]計(jì)算N組分氣體的黏度

(19)

其中，φij的計(jì)算公式為

其中，M(i)表示組分i的分子量。

1.2.6毛細(xì)壓力

氣液兩相的壓差稱為毛細(xì)壓力pc，即

pc=pg-pl

(20)

通常將毛細(xì)壓力簡(jiǎn)化為有效飽和度的函數(shù)，如Leverett給出的關(guān)系式為[23]

(21)

其中，σ為氣液兩相間的表面張力。

1.2.7比能和比焓

IAPWS協(xié)會(huì)選定水三相點(diǎn)的液相水作為參考狀態(tài)。三相點(diǎn)液相水的參數(shù)為[18]

(22)

三相點(diǎn)水蒸氣的比焓為[18]

(23)

本文采納上述規(guī)定，推導(dǎo)比能和比焓的計(jì)算公式。

假設(shè)液態(tài)水不可壓縮，其比能為

(24)

其中，cl為液態(tài)水的比熱容；Δul=cl·Tt。式(24)最后的近似條件是在所研究的溫度范圍內(nèi)，比熱容近似為常數(shù)。液態(tài)水的比焓為

(25)

水蒸氣的比焓為

(26)

其中，Δhg=cp,g(w)·Tt-hg,t(w)。水蒸氣的比能為

(27)

其中，ug,t(w)=hg,t(w)-RwTt為水蒸氣在三相點(diǎn)的比能；Δug=cv,g(w)·Tt-ug,t(w)=Δhg。式(27)的推導(dǎo)利用了理想氣體比熱容的性質(zhì)，即cp,g(κ)-cv,g(κ)=R(κ)。

根據(jù)在同一溫度下兩種物質(zhì)的比能之差與參考點(diǎn)無關(guān)的性質(zhì)，可得氣相組分κ的比能為

(28)

(29)

固體骨架的比能為

(30)

在得到各組分的比能和比焓后，混合氣體的比能和比焓分別為

(31)

在常溫下，氣體組分κ的比熱容為

(32)

其中，γ為比熱比。對(duì)于水蒸氣，γ(w)=4/3；對(duì)于空氣，γ(a)=7/5。示蹤物的比熱比與其性質(zhì)有關(guān)，如γ(CO)=7/5。

1.2.8熱導(dǎo)率

平行模型[13]將多孔材料等效為相互平行的圓柱形孔隙通道，假設(shè)熱流方向與多孔孔隙通道平行，則多孔介質(zhì)的表觀熱導(dǎo)率為

(33)

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

采用上節(jié)模型，編制了二維柱對(duì)稱坐標(biāo)系下求解氣液兩相非等溫滲流的計(jì)算程序。下面利用本文模型研究地下爆炸氣體的輸運(yùn)問題，并將計(jì)算結(jié)果與文獻(xiàn)[24]建立的氣相滲流模型的結(jié)果進(jìn)行比較，對(duì)兩種模型所得結(jié)果的異同性進(jìn)行解釋。

在處理水蒸氣-液態(tài)水的相變時(shí)，氣相滲流模型[24]假設(shè)水蒸氣在[Tb-ΔT,Tb+ΔT]間均勻液化，其中，Tb=373.15 K為水在1.01×105Pa的沸點(diǎn)，ΔT一般取為10 K。該模型沒有給出相變引起的質(zhì)量和能量源匯項(xiàng)，也未更新水蒸氣和空氣的份額變化。

2.1計(jì)算條件

假定在半徑為7 m的空腔中，初始時(shí)刻充滿了水蒸氣、空氣和一氧化碳，質(zhì)量分別為2.40×103， 3.65×103， 5.18×103kg， 氣體均勻分布， 溫度為2 000 K。研究域(含空腔)是半徑為86 m的球體，初始時(shí)刻與外界環(huán)境達(dá)到平衡。外界環(huán)境不含一氧化碳，氣體密度為1.02 kg·m-3，溫度為293.15 K。計(jì)算參數(shù)如表1所列。

表1　計(jì)算參數(shù)

2.2數(shù)值方法

圖1為網(wǎng)格剖分示意圖，計(jì)算區(qū)域?yàn)榍蛐?。任意點(diǎn)(θi,rj)實(shí)際代表圖中所示的虛線圓圈。坐標(biāo)原點(diǎn)位于空腔中心，圖中的紅圈代表空腔壁面，黑色網(wǎng)格代表多孔介質(zhì)。研究表明，空腔區(qū)域可等效為孔隙率為1、滲透率為相鄰多孔介質(zhì)的滲透率與孔隙度之比的多孔介質(zhì)[25]。采用交錯(cuò)網(wǎng)格的方式存儲(chǔ)物理量，即只有速度位于網(wǎng)格界面，其他物理量均位于網(wǎng)格中心。

圖1網(wǎng)格剖分示意圖Fig.1Schematic diagram of the grid

對(duì)運(yùn)動(dòng)方程，使用中心差分法離散，其中，對(duì)相對(duì)滲透率與黏度之比采用迎風(fēng)格式；對(duì)質(zhì)量和能量守恒方程使用有限體積法離散，其中，對(duì)對(duì)流項(xiàng)采用迎風(fēng)格式，對(duì)擴(kuò)散項(xiàng)采用中心差分格式。數(shù)值計(jì)算中還需要處理參數(shù)在界面處的等效值。例如，根據(jù)文獻(xiàn)[19]，當(dāng)運(yùn)動(dòng)方向與界面平行(或垂直)時(shí)，界面處的等效滲透率為相鄰邊界滲透率的加權(quán)平均值(或調(diào)和平均值)。

2.3結(jié)果與分析

2.3.1600s時(shí)的物理量分布

圖2給出了600s時(shí)本文模型和氣相滲流模型[24]的計(jì)算結(jié)果。可以看出，氣相壓強(qiáng)、氣相密度和滲流速度的分布基本吻合，因?yàn)榇藭r(shí)空腔內(nèi)尚未發(fā)生氣液轉(zhuǎn)變，未更新組分份額的影響也不太顯著。示蹤氣體的質(zhì)量分?jǐn)?shù)在空腔內(nèi)吻合是由于對(duì)流占主導(dǎo)地位，即氣體組分幾乎是以相同的比例運(yùn)動(dòng)的。滲透到多孔介質(zhì)內(nèi)的高溫氣體迅速降低到常溫，根據(jù)式(12)，水蒸氣的密度將降低、質(zhì)量分?jǐn)?shù)將相對(duì)減小，相應(yīng)的空氣和示蹤氣體的質(zhì)量分?jǐn)?shù)將增加。隨著時(shí)間的推進(jìn)，液態(tài)水開始出現(xiàn)，氣相份額的變化逐漸凸顯，兩種模型的差異越來越大，如圖3所示。

2.3.2空腔內(nèi)和r=7.5m處物理量的時(shí)間歷程

圖4和圖5分別為空腔內(nèi)和r=7.5m處物理量的時(shí)間歷程。圖4表明，在1.25×103s后，兩種模型的計(jì)算結(jié)果才有顯著差異，此時(shí)本文模型開始出現(xiàn)液態(tài)水，對(duì)應(yīng)的溫度為393.15K，水蒸氣分壓為1.99×105Pa，而氣相滲流模型[24]的相變發(fā)生在1.42×103s。液化釋放的潛熱將減緩溫度下降速度，故本文給出的空腔溫度更高。氣相滲流模型的相變過程非?？?，使得空腔壁面氣體流速迅速降低(見圖3(d))，故得到的氣相密度更大。此外，根據(jù)本文模型，在相態(tài)發(fā)生變化后，水蒸氣的質(zhì)量分?jǐn)?shù)減小，空氣和示蹤氣體的質(zhì)量分?jǐn)?shù)增加。圖5表明，在r=7.5m處，兩種模型給出的結(jié)果基本吻合。例如，溫度的差異始終在3K之內(nèi)，這與空腔的溫度變化范圍相比可以忽略。

(a)Pressureofthegasphase(b)Densityofthegasphase

(c)Massfractionofgastracer(d)Darcy’svelocityintheradialdirection

圖2t=600s時(shí)的各物理量分布

Fig.2Distributionofphysicalquantitiesatt=600s

(a)Pressureofthegasphase(b)Densityofthegasphase

(c)Massfractionofthegastracer(d)Darcy’svelocityintheradialdirection

圖3t=3 600s時(shí)的各物理量分布

Fig.3Distributionofphysicalquantitiesatt=3 600s

(a)Pressureofthegasphase(b)Densityofthegasphase

(c)Massfractionofthegastracer(d)Temperature

圖4空腔內(nèi)各物理量的時(shí)間歷程

Fig.4Timehistoryofphysicalquantitiesinthecavity

(a)Pressureofthegasphase(b)Densityofthegasphase

(c)Massfractionofthegastracer(d)Temperature

圖5r=7.5m處各物理量的時(shí)間歷程

Fig.5Timehistoryofphysicalquantitiesatr=7.5m

3小結(jié)

通過分析氣液兩相在多孔介質(zhì)中的相互作用過程，建立了描述兩相多組分非等溫滲流的數(shù)學(xué)模型。其中，針對(duì)高溫高壓環(huán)境，考慮了非達(dá)西流動(dòng)、水氣轉(zhuǎn)變以及輻射傳熱等特殊過程，且提供了一套完整的本構(gòu)關(guān)系使控制方程封閉。應(yīng)用該模型模擬了爆炸氣體的遷移行為，并與氣相滲流模型[24]的計(jì)算結(jié)果進(jìn)行了比較。對(duì)比結(jié)果表明，本文模型可以更合理地描述空腔內(nèi)的氣液相變過程和空腔外的氣體輸運(yùn)過程。

參考文獻(xiàn)

[1]孔祥言. 高等滲流力學(xué) [M]. 2版. 合肥: 中國(guó)科學(xué)技術(shù)大學(xué)出版社, 2010. (KONGXiang-yan.AdvancedMechanicsofFluidsinPorousMedia[M]. 2nded.Hefei:UniversityofScienceandTechnologyPress, 2010.)

[2]劉偉, 范愛武, 黃曉明. 多孔介質(zhì)傳熱質(zhì)理論與應(yīng)用 [M]. 北京: 科學(xué)出版社, 2006. (LIUWei,FANAi-wu,HUANGXiao-ming.TheoryandApplicationofHeatandMassTransferinPorousMedia[M].Beijing:SciencePress, 2006.

[3]PRUESSK.TOUGH2:Ageneral-purposenumericalsimulatorformultiphasefluidandheatflow[R].LBL-29400,LawrenceBerkeleyNationalLaboratory,UniversityofCalifornia, 1991.

[4]PRUESSK,OLDENBURGCM,MORIDISGJ.TOUGH2user’sguide[R]. 2nded.LBL-43134,LawrenceBerkeleyNationalLaboratory,UniversityofCalifornia, 1998.

[5]WUYS,PRUESSK.Numericalsimulationofnon-isothermalmultiphasetracertransportinheterogeneousfracturedporousmedia[J].AdvancesinWaterResources, 2000, 23(7): 699-723.

[6]CLASSH,HELMIGR,BASTIANP.Numericalsimulationofnon-isothermalmultiphasemulticomponentprocessesinporousmedia: 1.Anefficientsolutiontechnique[J].AdvancesinWaterResources, 2002, 25(5): 533-550.

[7]NITAOJJ.User’smanualfortheUSNTmoduleoftheNUFTcode[R]. 2nded.UCRL-MA-130653,LawrenceLivermoreNationalLaboratory, 1998.

[8]NITAOJJ.ReferencemanualfortheNUFTflowandtransportcode[R]. 2nded.UCRL-MA-130651-REV-1,LawrenceLivermoreNationalLaboratory, 2000.

[9]HELMIGR,BASTIANP,CLASSH,etal.ArchitectureofthemodularprogramsystemMUFTE-UGforsimulatingmultiphaseflowandtransportprocessesinheterogeneousporousmedia[J].MathematischeGeologie, 1998, 2(64): 123-131.

[10]劉昌軍. 非飽和水氣兩相滲流數(shù)值模擬研究及應(yīng)用 [D]. 南京: 河海大學(xué), 2005. (LIUChang-jun.Studyonnumericalsimulationofunsaturatedseepagebasedonwater-gastwo-phaseflowtheoryanditsapplication[D].Nanjing:HehaiUniversity, 2005.)

[11]劉曉麗, 梁冰, 王思敬, 等. 水氣二相滲流與雙重介質(zhì)變形的流固耦合數(shù)學(xué)模型 [J]. 水利學(xué)報(bào), 2005, 36(4): 405-412. (LIUXiao-li,LIANGBing,WANGSi-jing,etal.Fluid-solidcoupledmathematicalmodelforwater-airtwo-phaseinfiltrationanddeformationofdualporousmedia[J].JournalofHydraulicEngineering, 2005, 36(4): 405-412.)

[12]黃朝琴, 姚軍, 王月英, 等. 裂隙介質(zhì)兩相滲流數(shù)值模擬研究 [J]. 中國(guó)科學(xué): 技術(shù)科學(xué), 2011, 41(9): 1 240-1 248. (HUANGZhao-qin,YAOJun,WANGYue-ying,etal.Numericalstudyontwo-phaseflowthroughfracturedporousmedia[J].SciChinaTechSci, 2011, 41(9): 1 240-1 248.)

[13]FALTAR,PRUESSK,JAVANDELI,etal.Numericalmodelingofsteaminjectionfortheremovalofnonaqueousphaseliquidsfromthesubsurface: 1.Numericalformulation[J].WaterResourcesResearch, 1992, 28(2): 433-449.

[14]JAMBHEKARVA.Forchheimerporous-mediaflowmodels:Numericalinvestigationandcomparisonwithexperimentaldata[D].Stuttgart:Universit?tStuttgart, 2011.

[15]ZENGZW,GRIGGR.Acriterionfornon-Darcyflowinporousmedia[J].TransportinPorousMedia, 2006, 63(1): 57-69.

[16]王鐵良, 曹淵, 張建鑫. 地下爆炸空腔壓力和溫度歷程數(shù)值模擬 [J]. 計(jì)算物理, 2011, 28(5): 713-718. (WANGTie-liang,CAOYuan,ZHANGJian-xin.Numericalsimulationofcavitypressureandtemperatureinundergrounddetonation[J].ChineseJournalofComputationalPhysics, 2011, 28(5): 713-718.)

[17]楊世銘, 陶文銓. 傳熱學(xué) [M]. 3版. 北京: 高等教育出版社, 1998. (YANGShi-ming,TAOWen-quan.HeatTransfer[M]. 3rded.Beijing:HigherEducationPress, 1998.)

[18]WAGNERW,PRUSSA.TheIAPWSformulation1995forthethermodynamicpropertiesofordinarywatersubstanceforgeneralandscientificuse[J].JPhysChemRefData, 2002, 31(2): 387-535.

[19]王洪濤. 多孔介質(zhì)污染物遷移動(dòng)力學(xué) [M]. 北京: 高等教育出版社, 2008. (WANGHong-tao.DynamicsofFluidFlowandContaminantTransportinPorousMedia[M].Beijing:HigherEducationPress, 2008.)

[20]REIDRC,POLINGBE,PRAUSNITZJM.ThePropertiesofGasesandLiquids[M]. 4thed.NewYork:McGrawHill, 1987.

[21]MILLINGTONRC,QUIRKJP.Permeabilityofporoussolids[J].TransactionsoftheFaradaySociety, 1961, 57: 1 200-1 207.

[22]CHUEHC,BANGERTHW,DJILALIN.Integratedadaptivenumericalmethodsfortransienttwo-phaseflowinheterogeneousporousmedia[C]// 7thInternationalConferenceonComputationalFluidDynamics,Hawaii, 2012.

[23]LEVERETTMC.Capillarybehaviorinporoussolids[J].TransAIME, 1941, 142(1): 152-169.

[24]王鐵良, 張建鑫, 曹淵. 地下爆炸氣體輸運(yùn)的一種雙孔隙度數(shù)學(xué)模型 [J]. 應(yīng)用力學(xué)學(xué)報(bào), 2011, 28(6): 565-569. (WANGTie-liang,ZHANGJian-xin,CAOYuan.Adualporositymathematicalmodelofgastransportationonundergrounddetonation[J].ChineseJournalofAppliedMechanics, 2011, 28(6): 565-569.)

[25]王鐵良, 曹淵, 鄺飛虹, 等. 氣體在包含空洞的巖體中滲透的數(shù)學(xué)模型簡(jiǎn)化 [J]. 現(xiàn)代應(yīng)用物理, 2014, 5(1): 45-50. (WANGTie-liang,CAOYuan,KUANGFei-hong,etal.Simplifiedmathematicalmodelforgaspermeatingthroughhollowrockmass[J].ModernAppliedPhysics, 2014, 5(1): 45-50.)

收稿日期：2016-01-20；修回日期：2016-03-04

作者簡(jiǎn)介：何增(1991- )，男，陜西西安人，研究實(shí)習(xí)員，碩士研究生，主要從事計(jì)算流體力學(xué)研究。 E-mail:hezeng@nint.ac.cn

中圖分類號(hào)：O35

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：2095-6223(2016)021001(9)

Mathematical Model of Non-Isothermal Liquid-Gas Two-Phase Multi-Component Seepage Flow in Porous Media

HE Zeng,TIAN Zhou,WANG Tie-liang

(Northwest Institute of Nuclear Technology,Xi’an710024,China)

Abstract:A set of mathematical model is developed for non-isothermal two-phase multi-component seepage flow by analyzing the interactions of gas and liquid in porous media. Special processes are considered in high temperature and high pressure environment, such as non-Darcy flow, water transition and heat radiation, the convection, dispersion, source, and sink terms. Moreover, a complete set of constitutive relations is proposed to close the governing equations, among which the formulas for calculating the specific internal energy and the specific enthalpy are derived. Finally, this model is used to simulate the migration of explosive gas after underground detonation. The result suggests that, compared with the gas phase seepage model, this model can provide a more reasonable description for the gas transport in geological medium.

Key words:two-phase seepage flow；non-isothermal flow；mathematical model；constitutive relations；gas transport