高壓天然氣非恒定速率泄漏擴散數值模擬研究*

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

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

0 引言

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

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

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

1 數學模型及參數設置

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

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

(1)

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

FLACS數值模型經過多次改進與驗證，其結果具有可靠性和有效性。

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

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

1)學者提出使用假設的低壓泄漏源代替實際高壓射流源。等效噴嘴模型基礎原理源自Birch等提出的虛噴嘴模型，該模型可為氣體泄漏擴散規律研究提供可靠氣云分布數據。等效噴嘴模型假設：①氣體以質量流量m從泄漏孔流出，從高壓儲罐內部到實際泄漏孔出口為等熵膨脹過程；②實際泄漏孔處氣體流速為當地聲速；③氣體由實際泄漏孔流向虛擬泄漏孔過程滿足質量守恒和動量守恒方程；④氣體在虛擬泄漏孔處溫度與壓力已恢復到環境大氣壓與環境溫度。

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

(2)

(3)

(4)

式中：P1為實際泄漏孔處壓力，Pa；P0為高壓儲罐內壓力，Pa；γ為比熱比；T0為儲罐內溫度，K；T1為實際泄漏孔處溫度，K；Rn為氣體常數，Rn=519.625 J·kg-1·K-1；ρ1為實際泄漏孔處氣體密度，kg/m3；v1為實際泄漏孔處氣體速度，m/s；c1為當地聲速，m/s。

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

ρ1v1A1=ρ2v2A2

(5)

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

(6)

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

(7)

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

2)過程模型

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

(8)

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

(9)

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

(10)

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

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

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

(11)

由罐內氣體質量與罐內壓力(pi′)變化關系及氣流質量流量(qm′)得式(12)：

(12)

令

(13)

(14)

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

(15)

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

(16)

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

(17)

(18)

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

1.3 幾何模型及參數設置

1)幾何模型與網格尺寸

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

圖1 幾何及網格模型Fig.1 Geometric and mesh models

2)泄漏速率計算及參數設置

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

(19)

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

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

圖2 質量流量隨時間變化關系Fig.2 Change of mass flow rate with time

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

表1 泄漏工況Table 1 Leakage conditions

2 泄漏模擬結果分析

2.1 風速及泄漏時間對可燃氣云分布的影響

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

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

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

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

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

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

2.2 風速及泄漏時間對真實氣體等效化學計量氣云體積的影響

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

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

(20)

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

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

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

3 結論

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

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

3)風速對可燃氣云分布和體積影響較明顯；相同泄漏速率下，可燃氣云分布范圍和體積隨風速增大而縮小，風速變大會增大風向對氣云擴散的影響。