基于Realizable k-ε 湍流模型的氨氣泄漏數(shù)值模擬研究*

孫恩吉，李紅果，王　敏

(1. 中國安全生產(chǎn)科學研究院，北京 100012；2. 國家安全生產(chǎn)監(jiān)督管理總局信息研究院，北京 100012)

0　引言

隨著我國工業(yè)化生產(chǎn)和存儲的快速發(fā)展，現(xiàn)代化工藝的液氨制冷冷庫不斷增多，冷庫無論從規(guī)模上、設備上和自動化程度上都有了顯著提高。然而，我國仍存在很多老舊冷庫，這些冷庫由于設備陳舊、安全管理落后導致生產(chǎn)安全事故時有發(fā)生，安全生產(chǎn)形勢依舊嚴峻。2012年10月22日，湖北省洪湖市德炎水產(chǎn)公司發(fā)生氨氣泄漏事故[1]，造成400多人氨氣中毒搶救。2009年8月5日，內(nèi)蒙古赤峰制藥集團氨水配制車間液氨卸車過程中，發(fā)生氨氣泄漏[2]，致使202人接受治療。液氨制冷站制冷系統(tǒng)所需液氨儲存在儲罐及輸液管道內(nèi)，液氨輸運管道、閥門或儲罐容器壁等處意外破損、爆裂都會導致液氨大量泄漏，液氨泄漏后會迅速蒸發(fā)并在站內(nèi)迅速擴散[3]，對站內(nèi)及附近人員的人身安全將造成致命的危害。

采用計算流體動力學軟件來數(shù)值模擬泄漏擴散的過程，研究氣象條件、地表情況、泄漏源等對氣體擴散的影響是氣體泄漏擴散研究的重要方法之一。熊立春等[4]引入時間疊加的高斯液氨泄漏擴散模擬及人員疏散；J. Labovsky等[5]對Fladis試驗過程進行數(shù)值模擬，使得CFD動態(tài)模擬得到成功應用；Galeev等[6-7]模擬了氨氣擴散的過程，研究了連續(xù)釋放過程中風速和障礙物對擴散分布范圍的影響；S.M. Tauseef等[8]計算了在障礙物情況下重氣的擴散方式；鄭茂輝等[9-10]著重分析了建筑物的擾動和不同風速對氨氣濃度的影響；王志鵬[11]采用了高斯模型對液氨儲罐泄漏事故進行了研究，提出了泄漏后的快速處置方法；于加收[12]對不同風向條件下氨的泄漏擴散規(guī)律進行了研究；吳其榮等[13]模擬了火電廠區(qū)的氨泄漏及其擴散狀況。以往數(shù)值模擬更多的關(guān)注液氨儲罐事故與自然風影響下的區(qū)域擴散規(guī)律，針對室內(nèi)不同的泄漏點、不同的排風條件下氨氣的運移規(guī)律的研究相對較少。

因此，提出了基于Realizablek-ε的氨氣泄漏限元數(shù)值模擬分析方法，計算了液氨泄漏質(zhì)量，分析了增加排風口、不同液氨泄漏口及不同液氨泄漏量的氨氣擴散規(guī)律及濃度變化情況，對于設置應急處置裝備、采取合理措施具有指導意義。

1　模型的建立

1.1　氨氣泄漏物理模型

氨氣泄漏物理模型包括液氨儲罐，排風口A，B，進風口，空氣壓縮機3臺。液氨制冷車間尺寸為20 m×9 m。液氨儲罐為高5 m，直徑2.5 m的圓形罐體。壓縮機尺寸為1.5 m×1 m。進、排風口直徑尺寸均為1 m，分別位于左側(cè)的上下方及右側(cè)上方；通風口的尺寸為1 m，模型具體布置如圖1所示。

網(wǎng)格劃分是流體數(shù)值模擬的重要基礎，在模型離散化的過程中，模型被劃分成若干離散量元素。運用六邊形和非結(jié)構(gòu)化的網(wǎng)格劃分，將模擬區(qū)域劃分成3部分不同網(wǎng)格密度的計算區(qū)域，其中，泄漏口上方和中間泄漏區(qū)域采用最密集的四邊形網(wǎng)格，兩側(cè)區(qū)域網(wǎng)格為較疏的四邊形網(wǎng)格，網(wǎng)格包括大約500 000六邊形網(wǎng)格元素和400 000四邊形網(wǎng)格元素。以下關(guān)于液氨泄漏蒸發(fā)后的氨氣擴散過程模擬，選取其上部排風口為壓力出口,下方通風進口為速度入口邊界條件。

1.2　氨氣泄漏量的計算

在常壓下，液氨在空氣中會極速膨脹，大量氣化，并擴散到較大的空間范圍。液氨制冷車間正常運轉(zhuǎn)時，液氨儲罐泄漏發(fā)生在液相區(qū)域，在常溫、常壓下將以氣態(tài)形式存在，假設泄漏路徑較短，不會形成汽化核心而使部分液氨在泄漏路徑上汽化形成閃蒸兩相流，同時，假定液氨泄漏時全部被蒸發(fā)，而不考慮形成液池，其泄漏速率用伯努利方程計算：

(1)

式中:Q0為液氨泄漏速率，kg/s；Cd為泄漏系數(shù)；A為泄漏口橫截面積，m2；Pt為貯氨器內(nèi)壓，Pa；P0為大氣壓，Pa；h為泄漏點距液面高度，m；圓形泄漏口，其Cd取值為0.65，泄漏口直徑為0.01 m；240 K時，液氨密度ρt為681 kg/m3；儲罐內(nèi)壓Pt為950 kPa；大氣壓力P0為101 kPa；h取值為0.4 m；g取值為9.8 m/s2；液氨泄漏質(zhì)量初始流量為1.645 kg/s。

泄漏的液氨會發(fā)生閃蒸，泄漏時直接蒸發(fā)的液體所占百分比F可按下式計算：

(2)

式中：CP為液體的定壓比熱，J/kg·K；T為泄漏前液體的溫度，K；T0為液體在常壓下的沸點，K；H為液體的汽化熱，J/kg。根據(jù)經(jīng)驗，當F>0.2時，一般不會形成液池；當F<0.2時，F(xiàn)與蒸發(fā)液體之比，有線性關(guān)系，即當F=0時，沒有液氨被蒸發(fā)；當F=0.1時，有50%的液體被蒸發(fā)。

以下數(shù)值模擬主要針對制冷車間內(nèi)液氨儲罐短時間(10 s，30 s和60 s)泄漏的氨氣運移規(guī)律，根據(jù)制冷車間內(nèi)部的實際布置和工況，可作如下基本假設：

1)假設液氨儲罐環(huán)境溫度與儲罐初始溫度相同。

2)液氨儲罐泄漏口處的氨氣泄漏速度恒定。

3)液氨儲罐泄漏過程中氨氣作為理想氣體考慮。

4)泄漏過程不發(fā)生任何相變反應和化學反應。

5)進風口處的風速恒定。

6)泄漏過程中，溫度恒定且與外界無熱量交換。

1.3　Realizable k-ε模型

Realizablek-ε模型由Launder and Spalding提出后，出現(xiàn)了RNGk-ε模型及Realizablek-ε模型，被廣泛應用于工業(yè)流場及熱場的交換模擬中。Realizablek-ε模型為湍流粘性增加了1個限制條件公式，并為耗散率增加了傳輸方程，對旋轉(zhuǎn)流動、流動分離及復雜二次流都有更好的實現(xiàn)，對平板和圓柱射流發(fā)散比率有更精準的預測。利用Realizablek-ε模型在對氨氣泄漏過程中的射流現(xiàn)象的數(shù)值分析更能體現(xiàn)出優(yōu)越性，其表達式如下：

(3)

(4)

2　模擬結(jié)果分析

2.1　增加排風口的濃度分布

根據(jù)《呼吸防護用品的選擇、使用與維護》(GB/T18664-2002)中的規(guī)定，空氣中氨氣立即威脅生命或健康的濃度(IDLH)為360 mg/m3。選擇此濃度為臨界點，換算成體積摩爾濃度為2.12×10-5kmol/m3。模擬將排風口B設置為壓力出口，其他條件不變，增加排風口A時，計算泄漏時間為10 s，30 s，60 s時的氨氣濃度變化，結(jié)果如圖2所示。

圖2　氨泄漏10 s, 30 s, 60 s時氨氣IDLH區(qū)域(增加排風口A)Fig.2　Ammonia leak 10 s, 30 s, and 60 s when ammonia IDLH area (increased exhaust port A)

由圖2分析可以得出，增設排風口A，氨氣泄漏10 s后，氨氣由泄漏口垂直排入泄漏區(qū)，因其比空氣輕，上升至一定高度并開始在頂層飄逸。由于排風口A側(cè)下方為進風口，氨氣仍偏向排風口B側(cè)進行擴散，并擴散至泄漏區(qū)的頂部，在其到達頂部后，受到頂部阻力的影響，開始向下并沿四周空間擴散。在氨氣擴散和上升過程中，其高濃度區(qū)域不連續(xù)，并出現(xiàn)“擾流”現(xiàn)象。氨氣泄漏30 s后，泄漏氨氣向排風口B側(cè)擴散的更快速，B處氨氣濃度也相對較高，泄漏至頂部的氨氣沿頂部擴散的范圍和區(qū)域也有所增加，并且在頂部方向向排風口A側(cè)擴展。氨氣泄漏60 s后，泄漏氨氣氣流開始向排風口A側(cè)擴散，進風口對擴散的影響也逐漸加大。氨氣在到達B處頂部后，開始沿頂部及右側(cè)墻壁向下進行擴散。

2.2　不同泄漏口的濃度分布

將泄漏口位置設置在儲罐右上部，邊界條件設置為速度入口，泄漏口泄漏量為1.645 kg/s，將原泄漏口設置為固壁邊界，其他條件不變，模擬泄漏口位置不同、泄漏時間為10 s，30 s，60 s時的氨氣濃度變化情況，結(jié)果如圖3所示。

圖3　氨泄漏10 s, 30 s, 60 s時氨氣IDLH區(qū)域(改變泄漏口位置)Fig.3　Ammonia leak 10 s, 30 s, and 60 s ammonia IDLH area (changing leak port position)

改變泄漏位置，氨氣泄漏10 s后，泄漏氨氣氣流向排風口B側(cè)方向擴散，由于泄漏時間較短、泄漏速度大，氨氣并未向上擴散。氨氣泄漏30 s后，泄漏氣流開始向排風口B側(cè)擴散，受到第2臺壓縮機的阻擋，氨氣向廠房頂部擴散，儲罐與壓縮機間的氨氣濃度開始增大，氨氣在到達排風口B處濃度也開始逐步增加。氨氣泄漏60 s后，進風口的風速對泄漏擴散的影響增大。由于儲罐與壓縮機間無法進行氣體擴散，其間的氨氣濃度泄漏接近泄漏口處的氨氣濃度，氨氣氣流向排風口B側(cè)擴散趨勢明顯。

2.3　不同泄漏量的濃度分布

將泄漏量由1.645 kg/s分別增加至1.81 kg/s，1.99 kg/s，其他條件不變時，模擬泄漏區(qū)泄漏量變化，以及泄漏時間為10 s，30 s，60 s時的氨氣濃度變化情況。泄漏時間為10 s，泄漏量為1.81 kg/s，1.99 kg/s時氨氣濃度變化對比如圖4所示。

圖4　氨泄漏10 s時氨氣IDLH區(qū)域Fig.4　Ammonia leak 10 s ammonia IDLH area in different Leakage(1.81 kg/s upper panel, 1.99 kg/s lower panel)

泄漏量增加、泄漏時間10 s時，泄漏區(qū)氨氣最高摩爾濃度由0.397 kmol/m3增加至0.399 kmol/m3；泄漏量為1.99 kg/s時，氨氣濃度更容易上升至頂部，同時也主要集中在排風口B側(cè)。泄漏時間為30 s，泄漏量為1.81 kg/s，1.99 kg/s時氨氣濃度變化對比如圖5所示。

圖5　氨泄漏30 s時氨氣IDLH區(qū)域Fig.5　Ammonia leak 30 s ammonia IDLH area in different Leakage(1.81 kg/s upper panel, 1.99 kg/s lower panel)

泄漏區(qū)氨氣最高摩爾濃度由0.420 kmol/m3增加至0.426 kmol/m3，泄漏口豎直方向氨氣高濃度區(qū)域較多，排風口B處的氨氣濃度升高。泄漏時間為60 s，泄漏量為1.81 kg/s，1.99 kg/s時氨氣濃度變化對比如圖6所示。泄漏氨氣擴散到達上部頂板后，沿頂板向四周運動。

圖6　氨泄漏60 s時氨氣IDLH區(qū)域Fig.6　Ammonia leak 60 s ammonia IDLH area in different Leakage(1.81 kg/s upper panel, 1.99 kg/s lower panel)

3　結(jié)論

1)通過數(shù)值模擬分析發(fā)現(xiàn)，在增加排風口的情況下，氨氣泄漏60 s后，與進風口同側(cè)的排風口氨氣濃度增加幅度小于1%，氨氣主要朝對流方向的排風口飄逸。在不同泄漏口模擬條件下，由于氨氣氣流會受到廠房內(nèi)不同設備設施的阻擋，60 s后進風口的風速對氨氣泄漏擴散影響增大。不同氨氣泄漏量的情況下，通過1或2個湍流模型可產(chǎn)生更為合理的數(shù)值模擬結(jié)果。

2)數(shù)值模擬結(jié)果的重要性是所選擇的模型需要考慮在多相局部穩(wěn)定性的影響和介質(zhì)以及渦流擴散的非各向同性效應。氨氣的分子量小于空氣的平均分子量，對流排風口位置對于泄漏氨氣濃度影響較大，增加對流端的排風口有利于氨氣氣云的快速擴散和稀釋。氨氣泄漏量的增加使得泄漏口垂直方向氨氣濃度增加顯著，同時，空氣流動速度增加導致氨氣擴散趨于穩(wěn)定，能有效減小泄漏源造成的災害危險范圍。

3)基于Realizablek-ε湍流模型的氨氣泄漏有限元數(shù)值模擬可以有效的分析和預測氨氣泄漏后運移規(guī)律，該方法可為液氨制冷企業(yè)應急救援預案的制定、氨泄漏事故預防及應急處置提供參考。

[1] 張兵.冷庫建筑防火設計分析——以寶源豐禽業(yè)有限公司特別重大火災事故為例[J].中國安全生產(chǎn)科學技術(shù), 2013, 9(10):157-167.

ZHANG Bing. Analysis on fire safety design of cold storage building based on the lesson of fire disaster in BAO YUAN FENG Poultry Co. Ltd [J]. Journal of Safety Science and Technology, 2013, 9(10):157-167.

[2] 莊學強.大型液化天然氣儲罐泄漏擴散數(shù)值模擬[D] .武漢：武漢理工大學，2012.

[3] 趙明,陳求穩(wěn).城市重大危險源區(qū)域風險評價方法研究[J].中國安全生產(chǎn)科學技術(shù), 2014, 10(9):158-164.

ZHAO Ming, CHEN Qiuwen. Study on method of regional risk assessment for urban major hazard [J]. Journal of Safety Science and Technology, 2014, 10(9):158-164.

[4] 熊立春,陳建宏,石東平.引入時間疊加的高斯液氨泄漏擴散模擬及人員疏散研究[J].中國安全生產(chǎn)科學技術(shù), 2015, 11(11):76-82.

XIONG Lichun, CHEN Jianhong, SHI Dongping. Simulation on liquid ammonia leakage by Gaussian model introducing temporal superposition and personnel evacuation [J]. Journal of Safety Science and Technology, 2015, 11(11):76-82.

[5] J. Labovsky, L. Jelemensky. Verification of CFD pollution dispersion modeling based on experimental data [J]. Journal of Loss Prevention in the Process Industries, 2011, 24(2): 166-177.

[6] A. D. Galeev, E. V. Starovoytova, S.I. Ponikarov. Numerical simulation of the consequences of liquefied ammonia instantaneous release using FLUENT software [J]. Process Safety and Environmental Protection, 2013, 91:191-201.

[7] A. D. Galeev, A. A. Salin, S.I. Ponikarov. Consequence analysis of aqueous ammonia spill using computational fluid dynamics [J]. Journal of Loss Prevention in the Process Industries, 2013, 26: 628-638.

[8] S.M. Tauseef, D. Rashtchian. S.A. Abbasi. CFD-based simulation of dense gas dispersion in presence of obstacles [J]. Journal of Loss Prevention in the Process Industries, 2011, 24(4): 371-376.

[9] 鄭茂輝，金敏，許建明.城市建筑群環(huán)境有毒有害氣體擴散數(shù)值模擬[J].同濟大學學報自然科學版，2013, 41(1):48-52.

ZHENG Maohui, JIN Min, XU Jianming. Numerical simulation of hazardous gas dispersion around buildings in urban environment [J].Journal of Tongji University (Natural Science) 2013, 41(1):48-52.

[10] 鄭茂輝，郭月容.城市居住區(qū)毒氣擴散數(shù)值模擬[J].環(huán)境保護科學，2012, 38(1):1-4.

ZHENG Maohui, GUO Yuerong. Numerical simulation of toxic gas dispersion in urban residential areas [J].Environmental Protection Science, 2012, 38(1):1-4.

[11] 王志鵬.氨氣泄漏數(shù)值模擬及應急響應措施研究[D]. 北京；中國地質(zhì)大學，2013.

[12] 于加收.開放空間液氨泄漏擴散規(guī)律及人員疏散的研究[D].北京：中國地質(zhì)大學(北京)，2014.

[13] 吳其榮，杜云貴，喻江濤，等. 火電廠氨泄漏擴散特性模擬研究[J].環(huán)境科學與工程.2014, 37(2): 155-159.

WU Qirong，DU Yungui，YU Jiangtao, et al. Simulation of characteristics of ammonia leak diffusion in thermal power plant [J].Environmental Science & Technology, 2014, 37(2):155-159.