基于SLAB的重氣泄漏擴散模擬分析

何娟霞，劉 磊，周琪勇，李昂昂，馬 野，周冬梅

(1.廣西大學資源環境與材料學院，南寧 530004；2.華南理工大學材料科學與工程學院，廣州 510641)

汽油、乙醇、液氯等具有燃爆性、毒性、腐蝕性的危險化學品一旦發生泄漏，將會形成重質氣云。由于重質氣云密度比空氣大，在重力作用下，將會沉降到地面，并沿風向緩慢擴散。一旦發生火災爆炸、中毒窒息事故，將給周邊人群生命安全、大氣環境及水資源安全帶來嚴重危害[1-3]。

近年來眾多學者對重質氣云擴散分析采用SLAB模型。姜傳勝等[4]通過SLAB模擬與風洞實驗的對比，發現結果基本一致，并剖析了重氣連續擴散的特點；黃江平等[5]通過改進SLAB模型中風速廓線的計算方法，得到適用于模擬城市重氣擴散的SLAB_URBAN模型；朱紅亞等[6]通過對SLAB模型預測值與試驗觀測值的比較，發現數據吻合較好，并在SLAB模型基礎上發展了多源重氣泄漏擴散模型；劉昕等[7]將SLAB模型用于氨泄漏擴散后下風向濃度分布的模擬；熊發等[8]對SLAB模型大氣穩定度及地表粗糙度兩個參數進行了優化；Li等[9]考慮氣象環境尤其是風向、風速的時變性，提出了SLAB改進模型SLABi。眾多學者為SLAB后續研究奠定了基礎。

目前基于SLAB模擬下風向三維空間濃度分布的研究較少，在SLAB模型的基礎上，結合應急需求，從地平面濃度分布、敏感點濃度隨時間變化、濃度曲面隨高度變化3個方面對擴散數據進行分析。以便快速、有效預測事故空間影響區域、減少人員傷亡，為事故應急、環境保護提供科學依據。

1 重氣云團的形成與判定

1.1 重氣云團的形成

重氣云團(heavy gas, dense gas)是指分子量大于空氣的氣體形成的氣云或氣體本身分子量雖比空氣小，但由于從高壓、低溫環境泄漏等原因形成的比空氣重的氣云[10-12]。重質氣云的形成過程主要有①高壓、低溫液化氣體泄漏到大氣中，處于過熱狀態，迅速氣化形成夾雜有液滴的重質氣云，如液氯泄漏擴散；②常壓液態物質泄漏后，在地面集聚形成液池，進而揮發產生重質氣云，如汽油揮發擴散；③氣態物質直接在空氣中擴散形成重質氣云，如硫化氫井噴[13-14]。

1.2 重氣云團的判定

重氣云團一般用理查森數(Richardson number，Ri)判定[15-16]。理查森數原是研究大氣及海洋邊界層的重要無量綱參數，后被推廣到研究氣象學中的對流及大尺度運動中[17]。當初始理查森數Ri0大于臨界值R′，則認為是重氣云團，用重氣擴散模型對其擴散進行模擬；反之，則認為是非重氣云團。臨界值R′近似取50。Ri0計算公式為[12]

(1)

(2)

式中：g0′=(ρ0-ρa)g/ρa，m/s2；g為重力加速度，m/s2；ρ0為氣體初始密度，kg/m3；ρa為大氣密度，kg/m3；V0′為氣體連續泄漏時所形成的云羽的初始體積通量，m3/s；u為環境風速，m/s；D為泄漏源的特征水平尺寸，m；u*為摩擦速率，m/s；V0為氣體瞬時泄漏時所形成云團的初始體積，m3。

2 重氣擴散模型

2.1 概述

重氣擴散模型主要有唯象模型、箱模型、淺層模型、計算流體力學(computational fluid dynamics, CFD)模型等[18-19]，其中淺層模型結合了CFD模型與一維模型的優點，可以基本保持擴散控制方程的完備性，又能節約計算資源，提高計算效率，模擬結果吻合度較好[13]。

2.2 SLAB模型

SLAB模型由淺層模型發展而來，于1990年美國勞倫斯利佛莫爾國家實驗室(Lawrence Livermore National Laboratory，LLNL)開發。SLAB模型考慮到空氣卷吸、云團加熱、羽流上升、云團曲流等對重氣擴散的影響，通過求解質量、能量、動量、組分守恒等平衡方程，得到相應地點、時間的計算濃度。如表1所示，SLAB模型可對4種泄漏類型、3種氣云擴散類型進行分析，擴散類型取決于泄漏源持續泄漏時長(continuous source duration，TSD)。

表1 泄漏類型及氣云擴散類型Table 1 Types of leakage and gas diffusion

SLAB手冊[20](User’s Manual for SLAB)對SLAB模型做了詳細的說明，包括數學建模過程、SLAB程序的使用指南、SLAB程序使用案例，如圖1所示。數學建模過程包括穩態羽流擴散、瞬時煙團擴散、子模型及時間平均濃度求解。其中子模型包括環境流速剖面模型、卷吸率、熱與動量描述、熱力學模型、羽流上升模型。圖2所示為SLAB模型計算流程，首先確定泄漏源類型，然后選擇擴散模型，進而根據擴散方程計算空間平均濃度，最后得到時間平均濃度。

圖1 SLAB手冊架構Fig.1 Structure of SLAB manual

圖2 SLAB模型計算流程Fig.2 Calculation flow of SLAB

3 案例分析

3.1 泄漏場景介紹

以南寧某紙業公司氯氣鋼瓶整體破裂瞬時泄漏擴散為例。氯(chlorine)具有刺激性、氧化性、毒性，可對大氣、水體產生嚴重危害。美國工業衛生協會(American Industrial Hygiene Association，AIHA)制定的緊急反應計劃指南(Emergency Response Planning Guidelines，ERPGs)[21]將氯濃度水平劃分為3個等級，如表2所示。本案例取3、9、58、300 mg/m34個濃度標準，分別作為輕危區、中危區、重危區、致死區4個應急區域的劃分界限[15, 19]。

該公司周圍地形平坦，無高大建筑，廠內員工500人，周邊環境風險受體如表3所示。液氯鋼瓶容積為715 m3，最高工作壓力為1.1 MPa，正常工作溫度為-12 ℃，采用單層材質，最大儲存量為1 t。

表2 氯氣應急濃度標準Table 2 Emergency concentrations of chlorine

表3 泄漏點周邊環境風險受體Table 3 Risk receptors near leakage point

3.2 瞬時煙團模型控制方程

氯氣鋼瓶破裂瞬時泄漏擴散符合SLAB中的瞬時煙團模型。該模型根據質量、動量、能量、組分守恒方程，用空氣卷吸的概念衡量湍流對云團擴散的影響，將云團質心的擴散時間作為獨立變量。得到t時刻下風向(x,y,z)點的體積濃度值C(x,y,z,t)的計算式為

C(x,y,z,t)=4BxByhC(t)C1(x-

Xc,bx,βx)C1(y,by,βy)C2(z,Zc,σ)

(3)

t時刻的云團體積平均濃度C(t)為

(4)

水平濃度輪廓函數C1(y,by,βy)為

(5)

(6)

附加濃度輪廓函數C1(x-Xc,bx,βx)的對應法則與水平濃度輪廓函數C1(y,by,βy)相同，只需將對應的參數替換即可。垂直濃度輪廓函數C2(z,Zc,σ)為

(7)

(8)

式中：Bx、By分別為云團下風向與側風向半寬，m；bx、βx為云團下風向半寬參數；by、βy為云團側風向半寬參數；h為云團高度，m；m(t)為t時刻的云團質量濃度；Ma、Ms分別為干空氣、泄漏源物質的分子量；Xc為云團質心坐標，m；Zc為高度參數；erf為誤差函數。

3.3 參數確定

廣西壯族自治區南寧市屬于亞熱帶季風氣候，全年主導風向為東風、南風。表3所述泄漏點北側人口最為密集，考慮事故后果的嚴重程度，案例設定事發時風向為南風。SLAB模擬需要輸入的參數及其取值情況如表4所示。該環境條件下云團理查森數約為4 794，屬重氣云團。

3.4 模擬結果及分析

濃度模擬結果都是時間平均值。所謂某一點的時間平均濃度，是指在時間區間(tpk-TAV/2，tpk+TAV/2)內的濃度平均值，如式(9)所示。

(9)

式(9)中：Ctav為時間平均濃度，mg/m3；TAV為濃度平均時間，TAV不應大于云團持續時間TCD，s；tpk為(x,y,z)點濃度達到最大值時的時間，s；C(x,y,z,t)為(x,y,z)點在時間點t的濃度，mg/m3，見式(3)。

3.4.1 濃度分布曲線

根據已有數據點，調用MATLAB中griddata命令的cubic插值法，繪制下風向半邊濃度云圖，如圖3所示。下風向地面處(z=0 m)氯氣濃度分布情況如圖4所示，圖中由外到內4個閉合曲線代表由低到高的4個濃度標準。

結果顯示：廠內大部分位于致死區，濃度在300 mg/m3以上，且事故發生后氯氣云團隨風流擴散較快，廠區最遠點(下風向300 m處)濃度到達最大值只需85 s；A村處于事故點下風向位置，約50%面積位于重危區內，受氯氣云團影響較大，且居民較多，是本次事故的重點敏感區域；B村距事故地點1 200 m，但由于偏離了下風向位置，只有較少的面積受到影響，承受的風險遠小于A村；C、D、E村由于遠離下風向位置，幾乎不受事故的影響。分析可知，重氣云團對下風向區域影響最大，各濃度值影響區域都近似呈現狹長的橢圓形，廠區內人員及A村居民向垂直于下風向的方向逃生，可較快脫離危險區域。

3.4.2 濃度-時間曲線

A村中心點(下風向1 500 m處)的濃度-時間曲線如圖5所示。當泄漏發生后，150 s內空氣中氯氣濃度低于0.000 8 mg/m3；150 ～ 200 s期間，氯氣濃度逐漸緩慢升高；在200 ～ 285 s區間，空氣中氯氣濃度急速上升并達到最大值61.1 mg/m3；此后不斷下降，約在450 s左右降至3 mg/m3以下，對人體的傷害逐漸減輕。A村中心點暴露在50～60 mg/m3濃度區間內的時間約為60 s，能刺激人的呼吸道，引起咳嗽，可對健康造成一定的損害，沒有嚴重的生命威脅。

表4 SLAB模擬參數匯總Table 4 Summary of SLAB input parameters

圖3 半邊濃度云圖Fig.3 Half of concentration contour map

圖4 下風向地面處氯氣濃度平面分布Fig.4 Downwind concentration distribution of chlorine near the ground

圖5 A村中心點濃度-時間變化曲線Fig.5 Concentration-time curve in the A village

3.4.3 三維濃度曲面

圖6所示為各濃度三維散點輪廓圖。圖7所示為各濃度曲面側面輪廓線。可見濃度曲面呈狹長橢圓拱形。近距離處由于云團中氯氣的擴散，使濃度曲面的影響高度與影響寬度擴張；遠距離處受空氣卷吸與湍流影響，導致云團稀釋、濃度曲面收縮。各濃度最大影響高度分別為118、82、43、23 m。由于負浮力(重力)效應，加之泄漏源溫度較低引起的密度差，使云水平擴散距離遠大于縱向爬升高度，濃度曲面呈現扁平型，且近地處濃度曲面隨高度變化坡度較大，遠地處曲面較平緩。

圖6 各濃度三維散點輪廓圖Fig.6 Three-dimensional scatter diagram of all concentrations

圖7 各濃度曲面側面輪廓線Fig.7 Side profile of all concentrations

4 結論

利用SLAB重氣擴散模型，對某企業液氯鋼瓶瞬時泄漏擴散事故進行模擬，得到以下結論。

(1)得到泄漏源下風向三維濃度分布，對該類泄漏事故預防提供了科學依據，可為事故應急疏散、環境保護提供參考。

(2)根據計算結果，1 min左右致死區即可覆蓋整個廠區，致死區域寬度覆蓋144 m，空間距地高度最高達23 m，下風向最遠距離1 230 m。A村中心點到達濃度最大值的時間在285 s左右，顯著快于風速到達該處的時間，該點暴露在50 ～ 60 mg/m3濃度區間內的時間約為60 s。

(3)各濃度的三維輪廓曲面呈現狹長、扁平的橢圓拱形，且濃度不同，曲面高度、長度、寬度也不同。輕危、中危、重危、致死區的高度分別為118、82、43、23 m。