高壓天然氣非恒定速率泄漏擴(kuò)散數(shù)值模擬研究*

程方明，張安邦，王 燾，羅振敏，王 濤，陳 言

(1.西安科技大學(xué) 安全科學(xué)與工程學(xué)院，陜西 西安 710054；2.西安市城市公共安全與消防救援重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710054；3.兵器工業(yè)衛(wèi)生研究所，陜西 西安 710065)

0 引言

在天然氣工業(yè)中，壓縮天然氣(Compressed Natural Gas，CNG)生產(chǎn)、儲(chǔ)存、運(yùn)輸、使用一般在高壓容器中，且工藝涉及設(shè)備多，流程復(fù)雜。高壓容器及附屬設(shè)施在使用過(guò)程中常因腐蝕或操作失誤導(dǎo)致天然氣泄漏。天然氣儲(chǔ)配廠(chǎng)儲(chǔ)存大量高壓天然氣，一旦發(fā)生泄漏，遇明火會(huì)發(fā)生火災(zāi)或爆炸，對(duì)人員安全、公共設(shè)施及周邊環(huán)境造成嚴(yán)重威脅。因此研究?jī)?chǔ)配廠(chǎng)高壓天然氣泄漏，加強(qiáng)此類(lèi)事故的預(yù)防與控制很有必要。

國(guó)內(nèi)外學(xué)者對(duì)氣體泄漏與擴(kuò)散做了大量研究：Li等[1]通過(guò)使用FLACS軟件研究安全間距對(duì)氣體擴(kuò)散的影響發(fā)現(xiàn)，安全間距降低船型和圓柱形浮式液化天然氣(FLNG)結(jié)構(gòu)中的氣云尺寸，增大圓柱形FLNG結(jié)構(gòu)氣云尺寸；董剛等[2]采用FLUNT模擬高壓管道內(nèi)天然氣泄漏與擴(kuò)散過(guò)程，研究風(fēng)速對(duì)擴(kuò)散過(guò)程的影響；Li等[3]使用CFX軟件模擬LNG燃料船發(fā)動(dòng)機(jī)室中天然氣泄漏擴(kuò)散，優(yōu)化氣體探測(cè)器布局；Birch等[4-5]首次定義虛擬噴口，通過(guò)質(zhì)量和動(dòng)量守恒定律，推導(dǎo)出虛擬噴口處各有效物性參數(shù)，進(jìn)一步完善氣體泄漏模型；Chenoweth等[6-7]基于理想氣體模型，提出適用于A(yíng)bel-Noble狀態(tài)方程的泄漏模型；Chefer等[8]對(duì)理想氣體模型和Abel-Noble狀態(tài)方程的泄漏模型預(yù)測(cè)結(jié)果進(jìn)行對(duì)比發(fā)現(xiàn)，2者差異較??；劉凱迪[9]通過(guò)建立2區(qū)模型對(duì)高壓欠膨脹H2射流進(jìn)行數(shù)值模擬研究；劉自亮等[10]修正Brich泄漏模型，并結(jié)合CFD軟件模擬氣體泄漏及擴(kuò)散過(guò)程；李雪芳等[11]提出基于理想氣體狀態(tài)方程的理論模型和基于A(yíng)bel-Noble狀態(tài)方程的數(shù)值計(jì)算模型，并分析初始罐內(nèi)壓力為34.5，70 MPa時(shí)H2泄漏問(wèn)題，得到罐內(nèi)H2滯止參數(shù)隨時(shí)間變化曲線(xiàn)及出口氣流速度與流量隨時(shí)間變化曲線(xiàn)。

現(xiàn)有研究大都基于低壓和恒定泄漏速率條件，但實(shí)際氣體泄漏速率會(huì)隨儲(chǔ)罐壓力、溫度等參數(shù)變化而變化，且大多數(shù)氣體存儲(chǔ)設(shè)備壓力均高于環(huán)境壓力。當(dāng)儲(chǔ)存壓力與環(huán)境壓力之比大于氣體臨界壓力比(vcr)時(shí)，氣體泄漏就會(huì)產(chǎn)生欠膨脹射流。發(fā)生欠膨脹射流的泄漏口處速度會(huì)達(dá)到當(dāng)?shù)芈曀?，離開(kāi)泄漏口的氣體會(huì)以超音速繼續(xù)向外膨脹，最終在射流的遠(yuǎn)場(chǎng)區(qū)完全膨脹，氣體流動(dòng)完全發(fā)展。以某天然氣儲(chǔ)配廠(chǎng)為例，使用泄漏過(guò)程模型和等效噴嘴模型，利用FLACS軟件對(duì)儲(chǔ)配廠(chǎng)存儲(chǔ)壓力1.05 MPa高壓天然氣非恒定速率泄漏擴(kuò)散過(guò)程進(jìn)行數(shù)值模擬，考察環(huán)境風(fēng)速、泄漏時(shí)間對(duì)天然氣泄漏擴(kuò)散的影響，為研究欠膨脹射流高壓天然氣泄漏事故，采取相應(yīng)預(yù)防措施提供理論依據(jù)。

1 數(shù)學(xué)模型及參數(shù)設(shè)置

1.1 氣體泄漏擴(kuò)散數(shù)值模型及有效性驗(yàn)證

對(duì)于氣體泄漏擴(kuò)散，使用標(biāo)準(zhǔn)k-ε湍流模型，采用SIMPLE壓力修正法算法，并將其推廣到處理含附加源項(xiàng)的可壓縮流動(dòng)。通過(guò)建立描述流體特性的質(zhì)量、動(dòng)量、能量以及組分守恒方程，結(jié)合邊界條件求解計(jì)算區(qū)域中超壓、濃度、溫度等變量值，湍流和化學(xué)反應(yīng)影響包含在方程當(dāng)中[12]，方程如式(1)所示：

(1)

式中：Sφ為源項(xiàng)；φ為通用求解變量，包括質(zhì)量、動(dòng)量、能量等變量；ρ為氣體密度，kg/m3；xj為j方向積分；μi為i方向上速度矢量；Γφ為擴(kuò)散系數(shù)。

FLACS數(shù)值模型經(jīng)過(guò)多次改進(jìn)與驗(yàn)證，其結(jié)果具有可靠性和有效性。

2004年，Hanna等選用Kit Fox實(shí)驗(yàn)、MUST實(shí)驗(yàn)、Prairie Grass實(shí)驗(yàn)及EMU的L型建筑物風(fēng)洞數(shù)據(jù)集[13-17]，通過(guò)空氣質(zhì)量模型性能測(cè)量方法對(duì)FLACS氣體泄漏擴(kuò)散模型進(jìn)行有效性驗(yàn)證。結(jié)果表明：FLACS的CFD模型可以很好地預(yù)測(cè)氣體泄漏擴(kuò)散情況，其結(jié)果具有可靠性和有效性。2010年，Hansen[18]利用實(shí)驗(yàn)數(shù)據(jù)(China Lake,Burro,Coyote)驗(yàn)證CFD模型，在不同大氣條件下得到較精準(zhǔn)的測(cè)量值和預(yù)測(cè)值。

1.2 等效噴嘴及泄漏過(guò)程模型

1)學(xué)者提出使用假設(shè)的低壓泄漏源代替實(shí)際高壓射流源。等效噴嘴模型基礎(chǔ)原理源自Birch等提出的虛噴嘴模型，該模型可為氣體泄漏擴(kuò)散規(guī)律研究提供可靠氣云分布數(shù)據(jù)。等效噴嘴模型假設(shè)：①氣體以質(zhì)量流量m從泄漏孔流出，從高壓儲(chǔ)罐內(nèi)部到實(shí)際泄漏孔出口為等熵膨脹過(guò)程；②實(shí)際泄漏孔處氣體流速為當(dāng)?shù)芈曀?；③氣體由實(shí)際泄漏孔流向虛擬泄漏孔過(guò)程滿(mǎn)足質(zhì)量守恒和動(dòng)量守恒方程；④氣體在虛擬泄漏孔處溫度與壓力已恢復(fù)到環(huán)境大氣壓與環(huán)境溫度。

根據(jù)可壓縮理論，確定實(shí)際泄漏孔處條件如式(2)～(4)所示：

(2)

(3)

(4)

式中：P1為實(shí)際泄漏孔處壓力，Pa；P0為高壓儲(chǔ)罐內(nèi)壓力，Pa；γ為比熱比；T0為儲(chǔ)罐內(nèi)溫度，K；T1為實(shí)際泄漏孔處溫度，K；Rn為氣體常數(shù)，Rn=519.625 J·kg-1·K-1；ρ1為實(shí)際泄漏孔處氣體密度，kg/m3；v1為實(shí)際泄漏孔處氣體速度，m/s；c1為當(dāng)?shù)芈曀?，m/s。

氣體由實(shí)際泄漏孔流向虛擬泄漏孔滿(mǎn)足質(zhì)量守恒和動(dòng)量守恒方程，如式(5)～(6)所示：

ρ1v1A1=ρ2v2A2

(5)

ρ1v12A1+A1(P1-P2)=ρ2v22A2

(6)

聯(lián)立式(5)和式(6)可得虛擬泄漏孔處氣體速度與直徑，如式(7)所示：

(7)

式中：P2為虛擬泄漏孔處壓力，該壓力等于大氣壓力，Pa；v2為虛擬泄漏孔處速度，m/s；d1為實(shí)際泄漏孔直徑，mm；d2為虛擬泄漏孔直徑，mm。

2)過(guò)程模型

過(guò)程模型可根據(jù)儲(chǔ)罐內(nèi)初始?jí)毫?、溫度、?chǔ)罐體積和泄漏孔面積計(jì)算不同時(shí)刻泄漏孔壓力、溫度等參數(shù)。高壓氣體泄漏過(guò)程分2個(gè)階段：第1階段是超臨界流動(dòng)狀態(tài)，該階段儲(chǔ)罐內(nèi)氣體壓力與環(huán)境壓力比大于該氣體臨界壓力比；第2階段是亞臨界流動(dòng)狀態(tài)，該階段開(kāi)始標(biāo)志是儲(chǔ)罐內(nèi)壓力降為該氣體臨界壓力[11]。在第1階段，儲(chǔ)罐內(nèi)氣體壓力(Pi)隨時(shí)間變化如式(8)所示：

(8)

式中：常數(shù)如式(9)所示。

(9)

儲(chǔ)罐內(nèi)溫度(Ti)隨時(shí)間的變化如式(10)所示：

(10)

式中：Pi,0為儲(chǔ)罐內(nèi)初始?jí)毫?，Pa；Ti,0為儲(chǔ)罐內(nèi)初始溫度，K；vi,0為儲(chǔ)罐內(nèi)氣體初始比體積，m3/kg。

由式(2)～(4)可得，實(shí)際泄漏孔處壓力、溫度、密度、質(zhì)量流量等參數(shù)。因?yàn)镃>0,γ>1，所以?xún)?chǔ)罐內(nèi)壓力和溫度會(huì)隨時(shí)間增長(zhǎng)而下降。

在第2階段，泄漏孔處壓力與外部環(huán)境壓力相等，且不再隨時(shí)間變化而變化。設(shè)該時(shí)間為ttrans，如式(11)所示：

(11)

由罐內(nèi)氣體質(zhì)量與罐內(nèi)壓力(pi′)變化關(guān)系及氣流質(zhì)量流量(qm′)得式(12)：

(12)

令

(13)

(14)

聯(lián)立式(13)與式(14)，式(12)可簡(jiǎn)化為式(15)：

(15)

根據(jù)壓力微分方程，最終可得儲(chǔ)罐內(nèi)氣體壓力(Pi′)和溫度(Ti′)隨時(shí)間的變化[19]，如式(16)所示：

(16)

泄漏孔壓力為大氣壓力，即P1=Pa；泄漏孔處溫度(T1)和質(zhì)量流量(m)計(jì)算如式(17)～(18)所示[19]：

(17)

(18)

式中：Pa為大氣壓力，Pa；P1為泄漏孔處壓力，Pa；A1為泄漏孔面積，m2；Rn=519.625 J·kg-1·K-1，為氣體常數(shù)。

1.3 幾何模型及參數(shù)設(shè)置

1)幾何模型與網(wǎng)格尺寸

以某天然氣儲(chǔ)配廠(chǎng)為例，對(duì)罐區(qū)高壓天然氣泄漏進(jìn)行數(shù)值模擬。該儲(chǔ)配廠(chǎng)面積240 m×170 m，儲(chǔ)配廠(chǎng)天然氣球型罐共4臺(tái)，壓力1.05 MPa，各儲(chǔ)罐最大充裝系數(shù)0.8，依據(jù)現(xiàn)場(chǎng)實(shí)際情況按1∶1大小建立幾何模型?？傆?jì)算區(qū)域擴(kuò)大為280 m×210 m×50 m，為保證計(jì)算精度，節(jié)省計(jì)算時(shí)間，對(duì)泄漏點(diǎn)周?chē)W(wǎng)格進(jìn)行加密，在核心區(qū)域(60 m×90 m×20 m)采用大小一致立方體網(wǎng)格，在細(xì)化網(wǎng)格與均一網(wǎng)格之間使用光滑(Smooth)命令處理，對(duì)核心區(qū)域外網(wǎng)格進(jìn)行延伸處理，延伸比例小于1.2。根據(jù)FLACS軟件網(wǎng)格設(shè)置指導(dǎo)[12]及網(wǎng)格獨(dú)立性分析，設(shè)置氣體泄漏核心區(qū)域立方體網(wǎng)格尺寸為1 m，細(xì)化區(qū)域網(wǎng)格大小為0.35 m，網(wǎng)格X、Y、Z方向的相鄰控制體厚度遞差均為20.29%，滿(mǎn)足FLACS軟件計(jì)算網(wǎng)格要求。最終幾何模型及網(wǎng)格劃分如圖1所示，網(wǎng)格單元約49萬(wàn)個(gè)。

圖1 幾何及網(wǎng)格模型Fig.1 Geometric and mesh models

2)泄漏速率計(jì)算及參數(shù)設(shè)置

由于天然氣儲(chǔ)罐初始?jí)毫?.05 MPa，由式(19)可得，環(huán)境壓力與儲(chǔ)罐內(nèi)滯止壓力之比小于甲烷的臨界壓力比，發(fā)生泄漏時(shí)將產(chǎn)生欠膨脹射流。

(19)

式中：甲烷比熱比γ取1.3。

為研究較嚴(yán)重泄漏場(chǎng)景，泄漏孔徑設(shè)定100 mm。由式(1)～(3)可得實(shí)際泄漏孔處初始泄漏質(zhì)量流量15.54 kg/s；由式(4)～(6)可得初始虛擬泄漏孔直徑345.86 mm；結(jié)合過(guò)程模型，在T=200 s時(shí)，質(zhì)量流量14.83 kg/s，虛擬泄漏孔直徑336.49 mm。在T=7 558 s時(shí)，氣體泄漏進(jìn)入亞臨界流動(dòng)狀態(tài)，泄漏質(zhì)量流量降為3.11 kg/s。質(zhì)量流量隨時(shí)間變化曲線(xiàn)如圖2所示。

圖2 質(zhì)量流量隨時(shí)間變化關(guān)系Fig.2 Change of mass flow rate with time

結(jié)合模型及現(xiàn)場(chǎng)實(shí)際情況可知，儲(chǔ)罐上部與下部容易發(fā)生泄漏。天然氣主要組分為甲烷，密度比空氣小，當(dāng)儲(chǔ)罐上部發(fā)生泄漏，氣體隨泄漏初始動(dòng)能與空氣作用向下風(fēng)向和更高處擴(kuò)散，且高處出現(xiàn)點(diǎn)火源的可能性較小，所以本次泄漏位置選儲(chǔ)罐底部，泄漏方向?yàn)?Y方向，泄漏工況見(jiàn)表1。參考廠(chǎng)區(qū)所在地區(qū)氣象條件，設(shè)置風(fēng)向45°。根據(jù)風(fēng)向及FLACS軟件指導(dǎo)手冊(cè)，設(shè)置X，Y，Z正邊界條件為Wind，負(fù)邊界條件為Nozzle。

表1 泄漏工況Table 1 Leakage conditions

2 泄漏模擬結(jié)果分析

2.1 風(fēng)速及泄漏時(shí)間對(duì)可燃?xì)庠品植嫉挠绊?/h3>

不同風(fēng)速條件下T=200 s時(shí)可燃?xì)庠品植?維視圖如圖3所示(甲烷體積分?jǐn)?shù)為5%～15%)。由圖3可知，氣體泄漏后在近場(chǎng)區(qū)受初始動(dòng)能主導(dǎo)，從高壓儲(chǔ)罐泄漏的天然氣，本身具有很大動(dòng)能，即使使用虛擬泄漏孔代替實(shí)際泄漏孔，在虛擬泄漏孔處初始泄漏速度也能達(dá)239.97 m/s，初始動(dòng)能447.44 kJ。在初始動(dòng)能主導(dǎo)下向泄漏方向擴(kuò)散時(shí)，隨擴(kuò)散距離增加，需克服風(fēng)力等阻力，動(dòng)能逐漸減小，在距離泄漏源較遠(yuǎn)區(qū)域，氣云受風(fēng)力和浮力作用才開(kāi)始顯現(xiàn)。通過(guò)對(duì)比不同風(fēng)速條件下氣云分布可知，自然風(fēng)速增大可有效減少可燃?xì)庠品植紖^(qū)域，減小可燃?xì)庠企w積。

圖3 可燃?xì)庠品植?維視圖(T=200 s)Fig.3 3D view of combustible gas cloud distribution(T=200 s)

當(dāng)T分別為10，100，200 s時(shí)，對(duì)應(yīng)風(fēng)速條件下離地面1.5 m處XY平面可燃?xì)庠品植既鐖D4所示。為研究風(fēng)速對(duì)可燃?xì)庠品植嫉膭?dòng)態(tài)作用，分析同等高度下，不同時(shí)間對(duì)應(yīng)不同風(fēng)速條件下可燃?xì)庠品植甲铋L(zhǎng)距離，統(tǒng)計(jì)結(jié)果見(jiàn)表2。

圖4 不同時(shí)間對(duì)應(yīng)不同風(fēng)速條件下XY平面可燃?xì)庠品植糉ig.4 Distribution of combustible gas cloud on XY plane at different time and different wind speeds

由表2可知，同等高度條件下，風(fēng)速越大，可燃?xì)庠圃赮方向擴(kuò)散距離越短；在X方向，由于風(fēng)向與氣體泄漏方向不一致，氣云隨風(fēng)速增大，擴(kuò)散最遠(yuǎn)距離先增大后減小。從時(shí)間角度分析，可燃?xì)庠圃赥為10～100 s，擴(kuò)散最長(zhǎng)距離變化明顯，但在100～200 s擴(kuò)散過(guò)程中，可燃?xì)庠谱铋L(zhǎng)擴(kuò)散距離變化不明顯，最大變化距離1.76 m。結(jié)果表明，泄漏持續(xù)一段時(shí)間后，泄漏時(shí)間增加對(duì)可燃?xì)庠品植加绊戄^小，因?yàn)閺牡?節(jié)泄漏質(zhì)量流量計(jì)算結(jié)果可知，持續(xù)泄漏200 s內(nèi)，泄漏質(zhì)量流量?jī)H變化0.71 kg/s，對(duì)泄漏擴(kuò)散影響不明顯，所以各風(fēng)速條件下高壓天然氣變速泄漏擴(kuò)散會(huì)受風(fēng)力作用及浮力作用，使氣體擴(kuò)散在一定時(shí)間內(nèi)達(dá)到動(dòng)態(tài)穩(wěn)定；在相同泄漏速率下，風(fēng)速增大可縮短氣云趨于動(dòng)態(tài)穩(wěn)定時(shí)間，風(fēng)向?qū)τ谛孤U(kuò)散的影響越來(lái)越明顯，泄漏氣體會(huì)有明顯往下風(fēng)向擴(kuò)散的趨勢(shì)。

表2 不同時(shí)間對(duì)應(yīng)風(fēng)速條件下可燃?xì)庠品植甲铋L(zhǎng)距離統(tǒng)計(jì)Table 2 Statistics on longest distance of combustible gas cloud distribution at different time and different wind speeds

2.2 風(fēng)速及泄漏時(shí)間對(duì)真實(shí)氣體等效化學(xué)計(jì)量氣云體積的影響

等效方法可在模擬次數(shù)最少情況下得到爆炸載荷的均勻分布。分析等效化學(xué)計(jì)量氣云體積，對(duì)爆炸研究意義重大。在FLACS中用于評(píng)估爆炸后果的等效化學(xué)計(jì)量云有9種，其中Q8和Q9較為常用：Q8主要用于阻塞度較高的場(chǎng)景，Q9主要用于阻塞率低，通風(fēng)良好的開(kāi)放場(chǎng)景。本文模擬場(chǎng)景為天然氣儲(chǔ)罐區(qū)，四周空曠，所以選擇Q9進(jìn)行分析。Q9表達(dá)公式如式(20)所示[20]：

Q9=∑V×BV×E/(BV×E)

(20)

式中：V為可燃?xì)怏w體積，m3；BV為層流燃燒速度，m/s；E為在空氣中恒壓燃燒所產(chǎn)生的體積膨脹。

圖5 不同風(fēng)速條件下等效化學(xué)計(jì)量云體積Fig.5 Equivalent stoichiometric cloud volume under different wind speeds

4種風(fēng)速條件下等效化學(xué)計(jì)量云體積如圖5所示。由圖5可知，風(fēng)速為1.3 m/s時(shí)，等效化學(xué)計(jì)量氣云體積在60 s后不再發(fā)生明顯變化，在200 s時(shí)達(dá)到575.87 m3；風(fēng)速為3 m/s時(shí)，等效化學(xué)計(jì)量氣云的體積在40 s之后變化不明顯，基本穩(wěn)定在330.68 m3；風(fēng)速為5.0,8.0 m/s時(shí)，在15 s之后分別穩(wěn)定在184.69，104.19 m3。結(jié)果表明，在持續(xù)泄漏的200 s內(nèi)，真實(shí)氣體等效化學(xué)計(jì)量氣云體積在氣體泄漏初始動(dòng)能、穩(wěn)定風(fēng)場(chǎng)及浮力共同作用下，一定時(shí)間后趨于穩(wěn)定；在持續(xù)泄漏的200 s內(nèi)，達(dá)到動(dòng)態(tài)穩(wěn)定狀態(tài)后，泄漏時(shí)間增加對(duì)等效化學(xué)計(jì)量氣云體積幾乎無(wú)影響。

3 結(jié)論

1)儲(chǔ)存壓力為1.05 MPa天然氣儲(chǔ)罐發(fā)生泄漏將產(chǎn)生欠膨脹射流，根據(jù)等效噴嘴模型，使用虛擬泄漏孔代替實(shí)際泄漏孔，泄漏初期仍具有447.44 kJ的高動(dòng)能，氣體泄漏擴(kuò)散在近場(chǎng)受初始動(dòng)能主導(dǎo)。

2)結(jié)合過(guò)程模型，在持續(xù)泄漏的200 s內(nèi)，氣體泄漏質(zhì)量流量?jī)H產(chǎn)生0.71 kg/s變化，對(duì)氣體泄漏擴(kuò)散影響不明顯；該時(shí)間段內(nèi)，不同風(fēng)速條件下高壓天然氣變速泄漏，并在初始動(dòng)能、穩(wěn)定風(fēng)場(chǎng)及浮力共同作用下，使可燃?xì)庠企w積及分布在一定時(shí)間內(nèi)達(dá)到動(dòng)態(tài)穩(wěn)定狀態(tài)，等效化學(xué)計(jì)量氣云體積不再發(fā)生明顯變化；風(fēng)速增大可縮短動(dòng)態(tài)穩(wěn)定時(shí)間；隨泄漏時(shí)間進(jìn)一步增加，質(zhì)量流量會(huì)有明顯變化；當(dāng)T=7 558 s時(shí)，氣體泄漏將進(jìn)入亞臨界流動(dòng)狀態(tài)，質(zhì)量流量下降為3.11 kg/s，相比初始泄漏流量減少12.43 kg/s。

3)風(fēng)速對(duì)可燃?xì)庠品植己腕w積影響較明顯；相同泄漏速率下，可燃?xì)庠品植挤秶腕w積隨風(fēng)速增大而縮小，風(fēng)速變大會(huì)增大風(fēng)向?qū)庠茢U(kuò)散的影響。