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

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

0 引言

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

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

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

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

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

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

(1)

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

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

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

1.2 等效噴嘴及泄漏過程模型

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

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

(2)

(3)

(4)

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

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

ρ1v1A1=ρ2v2A2

(5)

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

(6)

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

(7)

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

2)過程模型

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

(8)

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

(9)

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

(10)

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

由式(2)～(4)可得，實際泄漏孔處壓力、溫度、密度、質(zhì)量流量等參數(shù)。因為C>0,γ>1，所以儲罐內(nèi)壓力和溫度會隨時間增長而下降。

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

(11)

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

(12)

令

(13)

(14)

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

(15)

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

(16)

泄漏孔壓力為大氣壓力，即P1=Pa；泄漏孔處溫度(T1)和質(zhì)量流量(m)計算如式(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)格尺寸

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

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

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

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

(19)

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

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

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

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

表1 泄漏工況Table 1 Leakage conditions

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

2.1 風(fēng)速及泄漏時間對可燃氣云分布的影響

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

圖3 可燃氣云分布3維視圖(T=200 s)Fig.3 3D view of combustible gas cloud distribution(T=200 s)

當(dāng)T分別為10，100，200 s時，對應(yīng)風(fēng)速條件下離地面1.5 m處XY平面可燃氣云分布如圖4所示。為研究風(fēng)速對可燃氣云分布的動態(tài)作用，分析同等高度下，不同時間對應(yīng)不同風(fēng)速條件下可燃氣云分布最長距離，統(tǒng)計結(jié)果見表2。

圖4 不同時間對應(yīng)不同風(fēng)速條件下XY平面可燃氣云分布Fig.4 Distribution of combustible gas cloud on XY plane at different time and different wind speeds

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

表2 不同時間對應(yīng)風(fēng)速條件下可燃氣云分布最長距離統(tǒng)計Table 2 Statistics on longest distance of combustible gas cloud distribution at different time and different wind speeds

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

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

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

(20)

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

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

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

3 結(jié)論

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

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

3)風(fēng)速對可燃氣云分布和體積影響較明顯；相同泄漏速率下，可燃氣云分布范圍和體積隨風(fēng)速增大而縮小，風(fēng)速變大會增大風(fēng)向?qū)庠茢U散的影響。