李文欣 黃啟玉 孫旭 張靖 孫超超

1中國石油大學（北京）油氣管道輸送安全國家工程實驗室·石油工程教育部重點實驗室

2中國石化工程建設有限公司

近年來，隨著外浮頂儲油罐不斷朝大型化發(fā)展，由密封圈密封不嚴或失效導致的油氣泄漏、起火爆炸事故逐年增加，不僅造成能源浪費，還會引起環(huán)境污染甚至人員傷亡。CHANG 和LIN[1]對國內(nèi)外529 起油罐火災事故進行統(tǒng)計表明，在雷擊引起的浮頂油罐火災事故中，因雷擊密封圈引起的約占83%。因此，有必要對大型外浮頂儲油罐密封圈油氣擴散規(guī)律進行研究，以便正確預測油氣分布規(guī)律、預估起火爆炸風險及范圍，為外浮頂罐火災事故的預防和救援提供理論支持。

在過去幾十年里，人們開展了多種工業(yè)場景下的氣體擴散規(guī)律研究。MAROTZKE、HUBER[2]、劉國梁[3]、吳晉湘[4]和施志榮[5]等開展了大量的風洞實驗，研究可燃氣體擴散規(guī)律和燃爆范圍。其中MAROTZKE 等搭建平地、斜坡等實驗場景，測取重氣繞各種不同類型障礙物的擴散結果；HUBER等考慮了長方體建筑物對流場結構及氣體擴散的影響，得出建筑物尾流區(qū)氣體濃度分布的數(shù)學模型；劉國梁等建立直流式風洞實驗系統(tǒng)，研究圍墻和樹對重氣擴散特征的影響；吳晉湘等研究了液化石油氣在有限空間內(nèi)的擴散過程，實驗結果表明，可燃區(qū)域隨泄漏時間推移而逐漸向上方發(fā)展；施志榮運用示蹤技術，在油氣儲運安全綜合實驗平臺上模擬事故現(xiàn)場的有害氣體的擴散規(guī)律。宋賢生[6]、王建[7]、趙剛[8]等采用雷諾時均（RANS）方程與湍流模型相結合的方式模擬流場，研究氣體濃度分布和障礙物對擴散過程的影響。其中宋賢生等利用CFX軟件模擬罐區(qū)油氣擴散過程，研究群罐和防火堤對油氣濃度分布的影響；王建等通過建立大型球罐區(qū)可燃氣體泄漏擴散和燃爆的數(shù)值模型，確定可燃范圍并對燃爆強度進行評估；趙剛等建立海上油氣鉆采平臺的三維模型，分析平臺結構、設備布局、風速對天然氣濃度分布的影響。竺柏康[9]、趙晨露[10]、郝慶芳[11]等采用小型試驗、數(shù)值模擬方法研究浮頂罐的油氣泄漏問題。其中竺柏康等建立浮頂罐試驗模型，發(fā)現(xiàn)風速比溫度對密封圈油氣濃度的不均勻分布影響更大；趙晨露、郝慶芳等采用RANS 方法與湍流模型相結合的方式模擬流場，研究不同風速、不同液位、不同浮盤泄漏孔隙位置的油氣擴散規(guī)律。

已有儲罐油氣擴散數(shù)值研究中，均采用了RANS 與湍流模型相結合的求解方法。由于RANS方法僅能求解流場時均量，不能準確模擬流場旋渦結構和油氣非定常耗散過程。為克服RANS 方法的缺點，采用大渦模擬（LES）方法，該方法可直接求解流場大尺度旋渦結構，更加準確描述油氣隨流場的輸運和耗散過程。GOUSSEAU 等[12]對比了RANS 和LES 的計算精度，發(fā)現(xiàn)隨著障礙物回流區(qū)對泄漏源影響的增大，RANS 模型對流通量的預測結果相較于LES 模型會產(chǎn)生較大誤差。本文采用LES 方法，能更加準確描述大型外浮頂罐繞流流場及密封圈泄漏油氣耗散過程。

1 控制方程及數(shù)值解法

根據(jù)大渦模擬的基本思想[13]，對不可壓縮流體的N-S 方程作濾波處理，可得大渦模擬控制方程，即

式中：上劃線為濾波后的場變量；下標i、j=1～3，為x、y、z三個方向的分量；ui、uj為氣流速度分量，m/s；xi、xj為坐標分量，m；t為時間，s；ρ為氣流密度，kg/m3；p為壓強，Pa；ν為氣流運動黏度，m2/s；τij為亞格子應力，m2/s2。

在控制方程中，和pˉ為待求量。由于式（3）中項不是和的顯式表示，τij也是未知量，因此方程組（1）～（3）不封閉，需要通過亞格子尺度模型建立τij與、的聯(lián)系。基于BOUSSINESQ 假設，亞格子應力可采用下式計算

式中：νSGS為亞格子渦黏系數(shù)，m2/s；為可解尺度應變率張量，s-1；τkk為亞格子應力的各向同性部分，m2/s2；δij為Kronecker 符號。

式（1）、（2）、（4）、（5）即為外浮頂罐繞流風場的控制方程。

在流過儲罐的氣流影響下，從密封圈泄漏出的油氣將向周圍環(huán)境擴散。由于油氣組分與空氣的濃度差較小，且環(huán)境溫差不大，忽略浮力的影響，對不可壓縮湍流中的組分輸運方程進行過濾，可得油氣擴散控制方程，即

式中：c為油氣濃度，kg/m3；λ為質量耗散系數(shù)，m2/s；為亞格子組分輸運過程，可通過建立與式（4）、（5）類似的亞格子模擬描述。

選用FLUENT 提供的DSM 模型（動態(tài)Smargorinsky-Lilly 模型）求解外浮頂罐油氣擴散問題。DSM 模型以SM 模型（Smargorinsky-Lilly 模型）為基礎，通過多次過濾，把湍流局部結構信息引入到亞格子應力中，可在計算過程中調整模型系數(shù)。該模型確定的亞格子渦擴散系數(shù)隨流場和標量場的不同而發(fā)生變化，可以克服SM 模型中將湍流普朗特數(shù)視作常數(shù)的缺點。大渦模擬求解過程具體設置如下：控制方程離散采用三維雙精度基于壓力的分離式求解器；壓力-速度耦合采用PISO 算法；空間離散中梯度項離散采用Least Squares Cell Based 格式，壓力項離散采用Standard 格式，動量方程離散采用Bounded Central Differencing 格式，油氣對流擴散項離散采用Second Order Upwind 格式，瞬態(tài)項求解采用Bounded Second Order Implicit 格式；壓力項求解亞松弛因子設置為0.01；動量求解亞松弛因子設置為0.05；甲烷求解亞松弛因子設置為1；無量綱時間步長設置為Δt*=Δt×u∞/(1.5H)=0.062（Δt為時間步長，s；u∞為來流速度，m/s；H為儲罐高度，m）；每個時間步長內(nèi)的最大迭代數(shù)設置為200；連續(xù)性方程、油氣濃度和x、y、z三個方向速度的殘差收斂準則均為10-5。

2 大渦模擬解法驗證

大型外浮頂罐的罐容大且儲存介質易燃易爆，難以通過現(xiàn)場試驗獲取強風作用下的密封圈油氣擴散數(shù)據(jù)。浮頂油罐的油氣擴散研究歸根結底是確定風場繞流障礙物時對油氣的裹挾、運移規(guī)律，因此基于HUBER 試驗，驗證本文大渦模擬解法的準確性。根據(jù)HUBER 等的風洞實驗建立數(shù)值計算模型，數(shù)值計算模型與實驗模型的尺寸比為1∶1（圖1）。在計算域為26H×14H×7H（H為長方體建筑物的高度，0.25 m）的流場中放置尺寸為H×H×2H的長方體建筑物（長邊垂直于來流方向）。建筑物迎風面距流場上游進口5H，建筑物背風側距流場下游出口20H，中心距流場兩側邊界7H，建筑物頂部距流場頂部6H。緊貼建筑物后墻面有1 個高為1.5H、直徑為0.042H的煙囪，煙囪頂部垂直向上排放泄漏氣體（甲烷和空氣的混合物，甲烷體積分數(shù)為1%），氣體泄漏速度為We=1.5u∞。流場進口自由來流對應的雷諾數(shù)為Re=6.4×104（基于H）。邊界條件設置如下：流場進口采用速度進口，流場出口采用自由出流邊界條件，流場兩側采用對稱邊界條件，流場頂部采用滑移邊界條件，煙囪頂部設置速度入口邊界條件，流場底部、建筑物表面和煙囪側面均采用無滑移邊界條件。如圖1 所示，在建筑物下游5H和10H處共設置6 條垂直于來流直線并檢測其上濃度變化。其中，LINE1、LINE4 垂直于地面，LINE2、LINE3、LINE5 和LINE6 平行于地面。LINE2 和LINE5 距地面1.5H，LINE3 和LINE6 在地面上。

圖1 氣體擴散基準模型Fig.1 Benchmark model of gas diffusion

采用高質量六面體結構化網(wǎng)格對流場進行網(wǎng)格剖分，所得網(wǎng)格包括4 565 210 個節(jié)點，如圖2 所示。根據(jù)HUBER 等的實驗，各計算參數(shù)設置如下：參考壓力為101 325 Pa；參考溫度為300 K；啟動組分輸運模型，混合物為空氣和甲烷（混合物之間不發(fā)生化學反應），空氣的密度為1.225 kg/m3、動力黏度為1.789 4×10-5Pa·s，甲烷的密度為0.667 9 kg/m3、動力黏度為1.087×10-5Pa·s。濃度場求解穩(wěn)定后，將t*（計算時間，無量綱）=600～900 的時均甲烷濃度無量綱化，與風洞實驗數(shù)據(jù)進行對比，檢驗大渦模擬計算結果的準確性。定義濃度系數(shù)K為

式中：χ為甲烷濃度，kg/m3；A為建筑物的截面積，m2；u∞為對應參考高度處的自由來流速度，m/s；Qe為甲烷的泄漏速率，kg/s。

圖2 氣體擴散基準模型的計算網(wǎng)格Fis.2 Computational grid for gas diffusion benchmark model

如圖3 所示，大渦模擬計算所得6 條檢測線上的濃度系數(shù)分布與HUBER 等的實驗結果吻合良好，驗證了本文油氣擴散模型和大渦模擬解法的準確性。

圖3 大渦模擬計算的濃度系數(shù)與風洞實驗對比Fig.3 Comparison of concentration coefficients calculated by LES and wind tunnel experiments

3 大型外浮頂罐密封圈油氣擴散

3.1 計算模型

采用驗證后的氣體擴散大渦模擬解法模擬10×104m3大型外浮頂罐的油氣擴散過程。大型外浮頂罐模型的流場計算域和密封圈泄漏區(qū)域劃分如圖4 所示，其中D為儲罐直徑，H為儲罐的罐體高度（H=0.275D），H*為儲罐最大安全液位（H*=0.25D），h為液位高度（即浮盤所在位置）。流場計算域范圍為15D×5D×4H，采用結構化網(wǎng)格對計算域進行空間離散。流場進口距儲罐中心5D，采用速度入口邊界條件；出口位于儲罐下游10D處，選用自由出流邊界條件；兩側邊界分別距離儲罐中心2.5D，施加對稱邊界條件；頂部邊界距離罐頂3H，采用滑移壁面邊界條件；流場底部和儲罐壁面選用無滑移邊界條件。

圖4 大型外浮頂罐模型的流場計算域和密封圈泄漏區(qū)域劃分Fig.4 Flow field calculation domain and seal leakage area division of large external floating-roof tank model

外浮頂罐的密封圈在日常操作、風吹日曬等影響因素下，會出現(xiàn)不同程度的磨損和破壞，密封圈的油氣泄漏位置和泄漏組分因罐而異，難以通過一個模型反映出所有儲罐的泄漏場景。為反映不同的泄漏情況，本文的數(shù)值模型將環(huán)形密封圈均勻分為4 份，選用最易揮發(fā)的輕烴組分CH4作為泄漏物，對應設置4 種泄漏位置。如圖4c、圖4d 所示，上風側為 -45°＜θ＜45° 的紅色區(qū)域；側風向為45°＜θ＜135° 的綠色區(qū)域和225°＜θ＜315° 的黃色區(qū)域；下風側為135°＜θ＜225°的藍色區(qū)域；整個環(huán)向為0°＜θ＜360°的整個環(huán)形密封圈區(qū)域（θ為泄漏點與迎風子午線的夾角）。儲罐蒸發(fā)損耗與油品性質、風速、所在地的大氣壓及密封形式等因素有關，現(xiàn)場數(shù)據(jù)不易獲取。因此，本文將某油田5×104m3原油外浮頂罐的年油氣損耗量折算為泄漏速度[14]，設置泄漏區(qū)域的速度入口邊界條件，以甲烷作為泄漏介質，假設其以u=0.03u∞的速度沿豎直方向從泄漏口逸出。對于現(xiàn)場實測數(shù)據(jù)已知的情況，可將實測泄漏速度與油氣組分代入，采用與本文相同計算流程求解。來流雷諾數(shù)固定在Re=1.64 × 106（基于儲罐直徑D），計算液位h=25%H*、50%H*、75%H*、100%H*和泄漏位置為上風側、側風向、下風側和整個環(huán)向等共16種工況。無量綱時間步長Δt*=Δt×u∞/D=0.037 5，待CH4濃度波動穩(wěn)定后，將t*=30～70 的CH4體積分數(shù)進行時均統(tǒng)計。

3.2 環(huán)向油氣濃度分布

圖5 給出了不同情況下密封圈上方2.5%H*高度處的CH4體積分數(shù)分布。如圖5a 所示，當油氣從上風側泄漏時，甲烷濃度沿環(huán)向緩慢升高，達到峰值后迅速減?。浑S液位升高，濃度峰值沿環(huán)向后移，但泄漏擴散范圍變化不大。如圖5b 所示，當油氣從側風向泄漏時，甲烷濃度沿兩端逐漸減小，且45°＜θ＜135°范圍內(nèi)的上風側甲烷濃度略高于下風側。液位高度為h=25%H*和h=50%H*時，在泄漏口上風側0°＜θ＜45°范圍內(nèi)充滿了較低濃度的甲烷，而液位高度為h=75%H*和h=100%H*時上風側的甲烷濃度幾乎為0。說明低液位時泄漏氣體在浮盤上方空間的擴散更顯著。如圖5c 所示，當油氣從下風側泄漏時，泄漏的甲烷聚集在泄漏口附近，不會向上游擴散；隨液位升高，濃度零值沿環(huán)向略有前移。如圖5d 所示，當油氣從整個環(huán)向泄漏時，液位較低時，上風側甲烷濃度高于下風側，濃度最高值位于θ=45°附近；滿液位時，甲烷氣體沿整個環(huán)向分布均勻。

為理解泄漏氣體在環(huán)向的分布規(guī)律，圖6 給出了不同液位高度的浮盤表面瞬態(tài)流線圖。如圖6a、圖6b 和圖6c 所示，儲罐液位較低時，浮盤上方的大部分氣流撞擊儲罐內(nèi)壁的后側，沿反方向運動，使下風側的油氣向上風側輸運，上風側泄漏的油氣仍聚集于此，而不進入其他區(qū)域。當油氣泄漏出現(xiàn)在側風向和下風側時，油氣則會沿流線向上風側擴散，小部分氣流撞擊儲罐內(nèi)壁前側再次產(chǎn)生回流，沿環(huán)向向兩側流動，將油氣向兩側輸運，因此，迎風子午線處（θ=0°）的油氣濃度不是最高值。這兩股氣流交匯于θ=45° 附近并產(chǎn)生旋渦，因此，當油氣從整個環(huán)向均勻泄漏時，油氣濃度最高值約位于θ=45° 。如圖6d 所示，儲罐接近滿液位時，氣流產(chǎn)生多股分流，在浮盤上方形成多處旋渦，因而此時的環(huán)向油氣濃度分布較均勻。

圖5 不同液位高度下泄漏口上方環(huán)向油氣濃度分布（ t*=30～70）Fig.5 Annular oil vapor concentration distribution above the leakage port at different liquid level heights（t*=30～70）

圖6 不同液位高度下浮盤表面瞬態(tài)流線（ t*=70）Fig.6 Transient flow line of floating roof surface at different liquid level heights（t*=70）

3.3 流向油氣濃度分布

圖7、圖8 分別給出了環(huán)向泄漏時z=0 剖面上的CH4體積分數(shù)分布和瞬態(tài)流場。如圖7a、圖8a 所示，當儲罐位于低液位h=25%H*時，在浮盤上方形成一個大尺度旋渦，渦量最大。在該旋渦的輸運作用下，從密封圈泄漏出的油氣可以充滿整個空腔。常用渦量來描述旋渦運動的大小和方向，流場的中央量為流體速度的旋度，ω*為無量綱渦量。如圖7b、圖7c、圖8b、圖8c 所示，隨著儲罐液位從h=25%H*升高到h=50%H*和h=75%H*，空腔大尺度旋渦分裂成多個小尺度旋渦，渦量顯著降低，旋渦對密封圈泄漏油氣的輸運作用減弱，因此，空腔內(nèi)的油氣擴散范圍逐漸減小，油氣主要聚集在密封圈附近。如圖7d、圖8d 所示，當儲罐液位為h=100%H*時，浮盤上方空腔不再被旋渦充滿，僅有浮盤上風側被旋渦覆蓋，由密封圈泄漏的油氣極易被輸運到主流中，因此，在浮盤表面只分布有少量油氣。

3.4 可燃區(qū)分析

甲烷的爆炸極限為5%～15%（體積分數(shù)），低于爆炸下限時不燃燒，高于爆炸上限后安靜燃燒。基于濃度計算結果，繪制爆炸下限CL=0.05 和爆炸上限CH=0.15 對應的時均濃度等值面圖，進一步分析儲罐浮盤上方空腔內(nèi)的三維油氣分布。如圖9 所示，離泄漏口較近的紅色曲面為爆炸上限等值面，包裹紅色曲面的藍色曲面為爆炸下限等值面。對爆炸下限的等值面進行面積積分，并除以相應的密封圈泄漏面積，可得單位泄漏面積對應的可燃區(qū)，用SL（無量綱）表示。如圖9、圖10 所示，低液位時（h=25%H*），浮盤上方大尺度旋渦的輸運作用使可燃區(qū)的范圍最大；當泄漏口位于下風側時，SL最高可達55。隨著液位增加，空腔內(nèi)旋渦輸運強度減弱，可燃區(qū)不斷減小。滿液位時（h=100%H*）的可燃區(qū)面積遠小于其他三個液位，當泄漏口位于下風側時SL達到最低值15。從圖9、圖10 還可以看出，當液位較低時，可燃區(qū)的面積差異較大，從下風側泄漏出來的油氣對應的SL最大，側風向次之，上風側最小。儲罐接近滿液位時，可燃區(qū)的面積大致相同，從側風向泄漏出來的油氣對應的SL最大，上風側次之，下風側最小。

圖7 不同液位高度下油氣瞬時濃度場（ t*=70）Fig.7 Instantaneous concentration fields of oil vapor at different liquid level heights（t*=70）

圖8 不同液位高度下儲罐繞流的瞬時渦量場和流線（t*=70）Fig.8 Instantaneous vorticity field and flow line of the flow around the rank at different liquid level heights（t*=70）

圖9 密封圈泄漏的甲烷時均濃度等值面（CL=0.05，CH=0.15）Fig.9 Isosurface of the mean concentration of methane leaked from the seal（CL=0.05，CH=0.15）

圖10 單位泄漏面積對應的可燃區(qū)Fig.10 Combustible area per unit leakage area

4 結論

利用大渦模擬高精度計算方法模擬10×104m3大型外浮頂罐密封圈泄漏油氣的擴散過程，得到不同液位儲罐瞬態(tài)繞流風場的旋渦分布和演變過程，重點分析泄漏位置和液位對罐頂油氣分布和可燃區(qū)的影響。模擬結果如下：

（1）當液位較低時，氣流在儲罐浮盤上方空腔內(nèi)形成大尺度旋渦，在旋渦的輸運作用下，泄漏油氣從下風側向上風側擴散，上風側油氣濃度高于下風側；隨著液位升高，浮盤上方空腔內(nèi)大尺度旋渦分裂為小尺度旋渦，輸運作用降低，油氣聚集在泄漏口附近。

（2）當油氣沿整個環(huán)向密封圈均勻泄漏時，油氣濃度最高值位于迎風面左右兩側45°附近，應及時維護和檢修這一區(qū)域的可燃液體監(jiān)測裝置。

（3）隨著液位升高，泄漏油氣在浮盤上的環(huán)向擴散作用減弱，主要聚集在泄漏口附近，可燃區(qū)減小；滿液位時，浮盤只有在靠近迎風面的部分區(qū)域被旋渦覆蓋，同時氣流將部分泄漏油氣帶離浮盤，可燃區(qū)達到最低值。

（4）當油氣分別從上風側、側風向和下風側泄漏時，低液位儲罐的可燃區(qū)范圍由大到小依次為下風側、側風向、上風側，滿液位儲罐的可燃區(qū)范圍由大到小依次為側風向、上風側、下風側，低液位儲罐的可燃區(qū)范圍比高液位儲罐高1.3～1.7 倍。