趙俊琪，錢琛庚，王成，孫遠翔

（1.北京理工大學(xué) 爆炸科學(xué)與技術(shù)國家重點實驗室，北京 100081；2.中國航空工業(yè)空氣動力研究院高速高雷諾數(shù)氣動力航空科技重點實驗室，遼寧，沈陽 110034）

氫能是一種清潔、制取來源豐富、燃燒熱值高的二次能源，但是氫氣的物理化學(xué)性質(zhì)決定了氫氣在制取、運輸、儲存、加注和使用等過程均存在泄漏和爆炸的風(fēng)險，使得公眾對于氫能安全問題充滿疑慮，研究加氫站內(nèi)的氫能安全事故，包括高壓氫氣泄漏擴散、氣云爆炸等，對于制定加氫站相關(guān)安全規(guī)范和標(biāo)準(zhǔn)、進行風(fēng)險評估及事故調(diào)查具有重要意義[1].高壓氫氣的危險性較高，實驗開展較為困難、可重復(fù)性差，已有的研究成果大多局限于較低的壓力或者在小尺度空間內(nèi)[2-4].而理論模型，如用于預(yù)測氣體擴散的Gauss 煙羽模型、Sutton 模型等及用于評估爆炸破壞作用的TNT 當(dāng)量模型、TNO 多能模型等，由于假設(shè)較多，僅能對簡單的事故過程進行再現(xiàn)，無法給出詳細的流場信息，也不能直觀地展現(xiàn)動態(tài)變化.因此對于加氫站氫能安全事故的大尺度空間、時間研究，基于計算流體力學(xué)（computational fluid dynamics, CFD）的數(shù)值模擬成為一種主流趨勢.目前常用的CFD 商業(yè)軟件包括Fluent、FLACS、CFX 等，國內(nèi)在工程尺度上有關(guān)分析事故發(fā)展規(guī)律、預(yù)測事故結(jié)果的研究幾乎基于以上軟件展開.但是CFD 軟件也存在缺點，國外商業(yè)軟件開發(fā)較早，數(shù)值方法精度最高為二階，依賴國外商業(yè)軟件難以實現(xiàn)大規(guī)模問題的高精度模擬仿真，因此自主研發(fā)復(fù)雜環(huán)境下可燃氣體泄漏-擴散-爆炸大規(guī)模高精度仿真軟件具有重要的工程應(yīng)用價值.

對于自然通風(fēng)情況下，加氫站高壓氫氣泄漏擴散及氣云爆炸的數(shù)值模擬，需要在計算域邊界施加風(fēng)，即亞聲速入流邊界條件，同時計算域邊界可能會接收到火焰面或多組分氣流，即亞聲速出流邊界條件.邊界條件通常被視為一個非關(guān)鍵的部分而被簡化處理，例如采用零階、一階外推方法.由于忽略了亞聲速流動的迎風(fēng)性問題，這種數(shù)值邊界條件通常不能等效于真實的物理邊界條件，導(dǎo)致高精度格式在大尺度、長時間模擬中產(chǎn)生的數(shù)值波與邊界相處作用而產(chǎn)生非物理振蕩[5]，邊界處的非物理反射還將向上游傳播及其他的邊界產(chǎn)生強數(shù)值耦合.最初POINSOT 等[5]提出了N-S 特征邊界條件（Navier-Stokes characteristic boundary conditions, NSCBC），通過對一維特征波進行分解，對進入、離開的特征波進行控制，可以顯著改善亞聲速入流、出流在邊界處的反射問題.后續(xù)，BAUM 等[6]通過將POINSOT 等[5]開發(fā)的一維特征波動關(guān)系的邊界公式擴展到用真實熱力學(xué)模型及詳細化學(xué)反應(yīng)模型描述的真實氣體混合物，給出了一種基于湍流反應(yīng)流的準(zhǔn)確、穩(wěn)定的NSCBC.目前廣泛應(yīng)用的NSCBC 應(yīng)用研究多數(shù)考慮變量在邊界上的時間推進，以計算邊界處的通量.例如，KIM 等[7]提出了在廣義坐標(biāo)系下高階格式的NSCBC 擴展保守形式，SUTHERLAND 等[8]主要研究了理想多組分氣體流動的邊界條件，DAVILLER 等[9]討論了入射特征波波幅的指定方法.在本文計算背景中，數(shù)值解設(shè)置在網(wǎng)格單元中心，計算域邊界位于兩個節(jié)點之間，如圖1 所示，若要計算流過邊界的數(shù)值通量，需要通過計算域內(nèi)節(jié)點值以及域外虛擬點值進行推導(dǎo)表示，無法直接使用傳統(tǒng)的NSCBC.目前關(guān)于使用虛擬點的NSCBC 研究較少，在研究中也未全面考慮附加項.例如，GROSS 等[10]針對單組份氣體推導(dǎo)了一種使用虛擬點的NSCBC，通過描述虛擬點內(nèi)的流動變量，使其產(chǎn)生的通量梯度滿足非反射條件，PIROZZOLI 等[11]建立了可壓縮流動的NSCBC，通過引入計算域外第一個虛擬點的流動信息參考值，使用一階偏心差分近似估計入射波特征波波幅，MOTHEAU 等[12]針對無化學(xué)反應(yīng)流動的亞聲速流動問題提出了一種基于NSCBC 給虛擬點賦值的方法.

圖1 通量離散方法Fig.1 Discretization procedure of the flux

因此，本文針對多組分反應(yīng)流，從特征變量的一維局部無黏方程出發(fā)，推導(dǎo)多組分反應(yīng)流Navier-Stokes 方程的特征波幅表達式，引入基于特征空間的反應(yīng)源項、橫向項，分析邊界處特征波的運動方向，使用計算域內(nèi)信息計算離開計算域邊界的特征波幅，使用目標(biāo)值松弛法計算進入計算域邊界的特征波幅.通過邊界處特征波幅計算原始變量梯度，結(jié)合高精度插值計算虛擬點處的守恒變量，發(fā)展基于虛擬點的反應(yīng)流N-S 特征無反射邊界條件方法.通過該程序可以更高效、精確、穩(wěn)定地實現(xiàn)通風(fēng)條件下大尺度加氫站氫氣泄漏擴散爆炸全過程的數(shù)值模擬，以給出環(huán)境風(fēng)對可燃氣云形成的影響規(guī)律以及可能出現(xiàn)的氣云爆炸事故后果.

1 基于虛擬點的Navier-Stokes 特征無反射邊界條件方法

N-S 特征邊界條件的核心思想是，在邊界處不進行強賦值以求解守恒變量，而是根據(jù)特征波的方向控制波的傳播，以更緩和的方式求解.在特征分析的過程中，為可以同時應(yīng)用于非反應(yīng)流和反應(yīng)流，引入了基于特征空間的化學(xué)反應(yīng)源項，為解決高維問題的模擬，引入了基于特征空間的橫向項.本文程序基于單元中心型均勻網(wǎng)格的有限差分方法展開，使用的插值模板為5 點模板，為求解計算域邊界通量梯度，需要在計算域外引入3 個虛擬點，而虛擬點的賦值方法與真實物理域中的流場信息有關(guān)，對于亞聲速入流、出流邊界應(yīng)用N-S 特征邊界條件方法給虛擬點賦值.

1.1 控制方程

對于包含Ns個組分，可能會產(chǎn)生Nr個化學(xué)反應(yīng)的氣體混合物，在笛卡爾坐標(biāo)系下，其非定常3D 可壓縮Navier–Stokes 控制方程組含有Ne=5+Ns-1個方程，表述如下

式中：U為守恒向量，與其對應(yīng)的原始變量向量為

F(U)為對流通量項；Fv(U)為黏性通量項；S為化學(xué)反應(yīng)源項.其中化學(xué)反應(yīng)速率的計算考慮一步化學(xué)反應(yīng)，采用Arrhenius 公式得到；下標(biāo)v表示黏性/擴散項；下標(biāo) α 、 β為x、y、z3 個方向的空間索引；下標(biāo)i為組分索引，此處對應(yīng)1 到Ns-1個組分； δαβ為克羅內(nèi)克符號.在控制方程中， ρ為流體質(zhì)量密度；u=(u,v,w)為速度向量；uα、uβ為速度分量；E為比總能量；Yi為組分i的質(zhì)量分數(shù)；p為壓力， τβα為黏性應(yīng)力張量；qβ為熱通量分量；Vi,β為組分i的質(zhì)量擴散速度分量； ω˙i為組分i的化學(xué)反應(yīng)生成速率.

本文基于多組分反應(yīng)流Navier-Stokes 控制方程組，采用大渦模擬，使用5 階WENO 格式離散對流通量項、守恒型4 階中心差分格式離散擴散通量項，使用3 階RK TVD 進行時間推進，給出高精度有限差分程序用于模擬.

1.2 數(shù)值方法

對于N-S 控制方程組，僅考慮一維問題并忽略擴散通量項，并將另外2 個方向上的空間偏導(dǎo)項合并簡寫為橫向項 T.以下為x方向上的推導(dǎo)，另外2個方向的結(jié)果完全對稱.

上式為基于守恒變量的一維控制方程組，給出x方向上守恒變量和原始變量之間的關(guān)系矩陣J以及對流通量項和原始變量之間的關(guān)系矩陣Bα，它們之間關(guān)系為

因此可改寫出基于原始變量的一維控制方程為

將非守恒雅可比矩陣F進行特征分解F=，得到特征值的對角矩陣 Λx以及對應(yīng)的左特征矩陣Vx、右特征矩陣.其中，特征值矩陣中的元素特征值λ在物理上表示特征波的傳播速度.定義特征波波幅Lx為

若忽略黏性項和源項，并以特征形式來重新描述控制方程，附加橫向項和反應(yīng)源項為

其中，基于特征變量的橫向項表示為

由于忽略體積力，只有對應(yīng)壓力和組分的源項不為0，并分別表示為sp、sYi.基于特征變量的反應(yīng)源項表示為

Lα的簡單物理含義可以描述為，除去基于特征變量的橫向項、反應(yīng)源項，為特征變量隨時間變化的相反數(shù).Lα為無化學(xué)反應(yīng)的N-S 方程組所計算出的特征波波幅，由POINSOT 等[5]給出

將橫向項、反應(yīng)源項擴展至特征波幅Lα，稱為，因此可表示出邊界處原始變量的空間梯度與的關(guān)系為

x方向左側(cè)邊界x1/2有最左側(cè)節(jié)點n=1 ，sgn=1，x方向右側(cè)邊界xnx+1/2有最右側(cè)節(jié)點n=nx，sgn=-1.通過梯度插值可得到虛擬點的值為

計算左、右邊界外虛擬點的值時，分別對應(yīng)使用計算域內(nèi)x1、xnx處的原始變量梯度.

1.3 特征波幅

1.3.1 亞音速無反射出流邊界條件

指定無窮遠處的靜壓p∞，當(dāng)計算域流出邊界的局部壓力與p∞不同時，僅部分無反射并將出口平均壓力松弛到p∞.考慮一個三維笛卡爾坐標(biāo)系下x方向邊界的流出，對于該邊界，L～φ(當(dāng)n=1 時 φ=5、sgn=1，當(dāng)n=nx時 φ=1、sgn=-1)進入計算域，而其他4+Ns-1條特征波動離開計算域.對于L～φ，施加了一個物理邊界條件，無窮遠處的靜壓，將波動振幅與壓差p-p∞關(guān)聯(lián)，并擴展包括源項、橫向項為

式中：K為壓力松弛系數(shù)； η為橫向項貢獻系數(shù).根據(jù)RUDY 和STRIKWERDA[13]的建議，壓力松弛系數(shù)可以表示為

式中：Ma為邊界處流動的平均馬赫數(shù)；l為計算域的特征尺寸，在本文中取為網(wǎng)格尺寸；c為平均當(dāng)?shù)芈曀伲?σ為人工系數(shù)，與 β的取值一樣高度依賴于流動特性.當(dāng) σ較小時，平均壓力可能出現(xiàn)漂移，而 σ越大時，邊界上的反射強度會更大.針對不同的計算情況，σ應(yīng)該取不同的值，主要依靠經(jīng)驗選取.SELLE 等[14]提出，為避免平均壓力漂移或較大反射，系數(shù)范圍可以取為1 ＜σ ＜π.

1.3.2 亞音速無反射入流邊界條件

對于亞音速無反射入流，有4+Ns-1條特征波進入計算域，僅有一條特征波Lφ離開計算域.直接施加目標(biāo)入流速度可能會導(dǎo)致聲波反射并與入流邊界相互作用，為了在邊界上維持想要施加的目標(biāo)速度（ut,vt,wt）、溫度（Tt）和質(zhì)量分數(shù)（Yit），可以在入流邊界處通過對目標(biāo)值進行松弛，來控制特征波幅并反映部分非反射的特性,

式中：松弛系數(shù)σm(m=1,2,···,5+Ns-1)用于控制施加到計算域的波動振幅阻尼； ηm用于控制橫向項的貢獻量.

1.4 典型算例驗證

1.4.1 非反應(yīng)流

該案例為在x方向均勻流場上疊加的單個渦旋，為評估邊界條件方法性能的常用案例，此外在非反應(yīng)流的案例中不引入基于特征空間的化學(xué)反應(yīng)源項.正方形計算域長度為L=13 mm，網(wǎng)格尺寸為Δx=Δy=0.1mm，初始時為常溫Tref=300 K、常壓pref=101.325kPa 下的氮氣.初始流函數(shù)的定義式為

式中：(x0,y0)為渦旋中心；Γ=0.11m2/s和Rv=0.1L=1.3mm 分別為渦旋強度和半徑.因此速度場定義為

式中：u0=100 m/s為初始流速.初始壓力場表示為

4 個計算域邊界均分別使用一階外推方法、亞聲速無反射出流邊界條件方法，無反射邊界條件中參數(shù)σ=0.25，橫向項參數(shù)為η=0.3.

圖2 明顯表現(xiàn)出來2 種方法的差異性，當(dāng)渦旋開始與邊界相互作用時，一階外推方法使渦旋在計算域邊界變形，通過邊界的過程中渦旋變形扭曲更加嚴重，開始發(fā)散，并在邊界處產(chǎn)生強烈的壓力反射返回計算域，渦旋在離開計算域時非物理反射依然存在于計算域中.而使用虛擬點的NSCBC 使渦旋在離開計算域邊界過程中幾乎沒有變形，也達到了幾何沒有反射波回到計算域的效果.可見亞聲速無反射邊界條件可以良好地抑制高維效應(yīng)帶來的壓力、聲波反射.

圖2 t=4, 26, 40, 54 μs 速度等值線和歸一化壓力場對比結(jié)果Fig.2 Velocity contours and normalized pressure field at t=4, 26, 40, 54 μs

1.4.2 反應(yīng)流

該案例計算域長度為L=3.0 mm，網(wǎng)格尺寸為Δx=Δy=12.5μm.初始時計算域有球心（1.7 mm, 1.7 mm），半徑1 mm 的點火源，為當(dāng)量比為1 的氫氣-空氣混合物的完全反應(yīng)產(chǎn)物，混合物溫度為2 250 K，計算域其他部分均為均勻預(yù)混的當(dāng)量比為1 的氫氣-空氣混合物，混合物溫度為300 K，全場壓力為101.325 kPa.3 個方向上每個邊界均分別使用一階外推方法、亞聲速無反射出流邊界條件方法，在亞聲速無反射出流邊界條件中參數(shù)σ=0.25.

圖3 顯示了二維預(yù)混火焰流出計算域不同時刻溫度等值面對比結(jié)果.非反射的邊界條件方法保持穩(wěn)定，而一階外推給出了不符合真實火焰流動的結(jié)果.在2 μs 時，火焰還未到達計算域邊界，兩種方法給出的結(jié)果一致；在7 μs 時，火焰面接近計算域邊界，一階外推對火焰具有拉伸效果，而無反射條件對于火焰?zhèn)鞑o影響；在15 μs 時，一階外推的溫度場受到邊界的反射而產(chǎn)生反轉(zhuǎn)；在28 μs 時，無反射條件的火焰面已經(jīng)接觸到右側(cè)和上側(cè)的計算域邊界，溫度場結(jié)構(gòu)保持穩(wěn)定不變，而一階外推導(dǎo)致溫度場被破壞，產(chǎn)生較大的逆溫度梯度，氣流在邊界的流動受到阻滯，產(chǎn)生了非物理的結(jié)果.

圖3 t=2, 7, 15, 28 μs 溫度等值線對比結(jié)果Fig.3 Temperature contours at t=2, 7, 15, 28 μs

2 加氫站氫氣泄漏—擴散—爆炸全鏈條數(shù)值模擬

2.1 計算模型

對某真實加氫站進行合理簡化，建立幾何模型.以加氫站的西南角為原點，建立三維直角坐標(biāo)系，向東為X軸正方向，向北為Y軸正方向，垂直向上為Z軸正方向，將計算域取為51 m×51 m×10 m，幾何模型如圖4 所示.選擇發(fā)生頻率更高的小規(guī)模連續(xù)泄漏事件作為數(shù)值模擬分析的對象，設(shè)定可能發(fā)生高壓氫氣泄漏的地方包括長管拖車A末端面板的泄氣軟管和加注機B的加注軟管.其中，均考慮配件或管道的100%損壞破裂.表1 為數(shù)值模擬考慮的泄漏場景[15 - 16].

表1 泄漏場景參數(shù)Tab.1 Parameters of release scenario

圖4 加氫站幾何模型Fig.4 Simplified geometric model of hydrogen station

對于高壓氫氣的泄漏擴散，除了內(nèi)在因素如泄漏孔徑、泄漏壓力及泄漏高度等，外部因素也對氫氣可燃氣云的形成有較大影響.其中，風(fēng)向會影響泄漏擴散的主要方向，風(fēng)速會影響擴散速度，因此泄漏氣體的擴散過程與環(huán)境風(fēng)況密切相關(guān).參考當(dāng)?shù)氐恼鎸崥夂驐l件[17]，考慮采用平均氣溫、氣壓，將模擬環(huán)境參數(shù)列于表2.

表2 環(huán)境參數(shù)Tab.2 Parameters of environment conditions

假設(shè)所有設(shè)備以及地面在事故過程中不發(fā)生破損，采用絕熱無滑移固壁邊界條件.通風(fēng)邊界采用NSCBC 亞聲速入流邊界條件，其他方向的計算域邊界均采用NSCBC 亞聲速出流邊界條件.將長管拖車的泄漏方向考慮為水平向西朝向附近的防火墻，泄漏位置為(7.4, 25.5, 1.6) m，將加注機的泄漏方向考慮為垂直向上朝向罩棚，泄漏位置為 (29.1, 32.5, 2.5) m.由于不同臨氫設(shè)備A、B的泄漏口大小不同，分別建立不同網(wǎng)格，采用Molkov 虛噴管模型對高壓氫氣射流進行簡化模擬[18]，在泄漏口處給出等效的入流條件，因此計算域網(wǎng)格尺寸放大至分別為90、140 mm，計算域網(wǎng)格總量分別為3 568 萬、947.8 萬，采用256核并行計算.

2.2 高壓氫氣泄漏擴散

高壓氫氣泄漏擴散的氣流運動十分復(fù)雜，在泄漏口附近，氣流主要受動量主導(dǎo)，沿泄漏口軸線方向氫氣濃度高速上升，而遠場區(qū)域，周圍障礙物形狀以及通風(fēng)條件等多種環(huán)境因素造成的湍流對氫氣擴散影響更大[19].環(huán)境風(fēng)對可燃氣云形成可以產(chǎn)生兩種截然不同的影響，既可以削弱自然擴散、抑制氫氣與空氣的混合并縮短擴散距離，也可以推動可燃氣云流動、增強湍流、稀釋氫氣并增大擴散距離.因此，泄漏形成的高速氣流在障礙物和環(huán)境風(fēng)共同作用下形成許多小尺度渦旋，不僅具有橫向脈動，還存在與氣流主方向相反的運動方向，使氣流動量、質(zhì)量傳遞速率急速上升.給出10 種工況下可燃氣云危險距離、可燃氣云體積的變化趨勢如圖5 所示，其中考慮人體高度，危險距離為1.5 m 高度水平方向估計.可看出水平方向危險距離總體滿足隨風(fēng)速增大而增大的趨勢，但B2在低風(fēng)速時危險距離反而比無風(fēng)時更低，這是環(huán)境風(fēng)對自然擴散的抑制造成的.對于泄漏口A/B，不同風(fēng)向時可燃氣云體積的變化趨勢相反，可見風(fēng)向?qū)τ诃h(huán)境風(fēng)稀釋可燃氣云的效果起到極大影響.

圖5 t=8 s 時所有工況可燃氣云z=1.5 m 的水平方向危險距離及氣云體積對比Fig.5 Horizontal dangerous distance at z=1.5 m and gas cloud volume of all simulation conditions at t=8 s

選取危險距離及可燃氣云體積綜合最大的工況B5進行分析，圖6 為不同時刻的可燃氣云濃度云圖.泄漏初期，射流由動量控制，可燃氣云分布不受環(huán)境風(fēng)影響，豎直向上的射流與罩棚頂部相互作用，射流速度方向改變.泄漏后期，可燃氣云在罩棚下大量聚集，隨著風(fēng)速增加，可燃氣云擴散路徑被改變，Y方向上的擴散距離小幅度減少，X方向上的擴散距離大幅度增加.東風(fēng)使可燃氣云擴散到罩棚外，傾向從壓縮機和長管拖車之間的通道運動.

圖6 不同時刻B5 可燃氣云輪廓Fig.6 Flammable gas cloud of B5 at different time

圖7 為不同時刻氫氣體積分數(shù)YH2的水平切片（z=1.5 m）及垂直切片（y=32.5 m）.由于風(fēng)速高，且下風(fēng)向存在復(fù)雜障礙物，從流線可看出各個障礙物之間存在大氣擾動并產(chǎn)生渦旋.風(fēng)場使得氣流基本朝向–X方向流動，受到障礙物的影響，主要分布于障礙物之間的通道中，造成了壓縮機附近大量氫氣聚集，因此潛在危險更大.后續(xù)因為泄漏氫氣總量不斷增加，罩棚下的可燃氣云濃度升高，可燃氣云在空中被延長，由于湍流擴散能力強，濃度較低呈現(xiàn)出不連續(xù)的碎片狀，最高到達7.18 m.

圖7 不同時刻可燃氣云z=1.5 m 水平方向切片、x=29.1 m 垂直方向切片的體積分數(shù)分布對比Fig.7 Volume fractions contour of horizontal slice at z=1.5 m, vertical slice at x=29.1 m at different time

總的來說，環(huán)境風(fēng)速對泄漏擴散的影響較為復(fù)雜.更高的風(fēng)速將抑制可燃氣云在某些方向的自然擴散，如垂直風(fēng)向的水平擴散以及豎直向上的運動等.但是風(fēng)速越高，大氣湍流越強，空氣的稀釋能力增強，環(huán)境風(fēng)對可燃氣云的推動作用也越強，可燃氣云加速擴散，在下風(fēng)向同一位置處的氫氣濃度越低.由于這兩種作用相互矛盾，導(dǎo)致風(fēng)速對可燃氣云擴散分布的影響規(guī)律呈現(xiàn)非線性，在較高風(fēng)速下第二種作用會占主導(dǎo)因素，而其余情況下則需要結(jié)合風(fēng)向、障礙物分布綜合考慮.

2.3 非均勻氣云爆炸

相對于常見的可燃氣云爆炸研究，本文由真實場景下泄漏擴散模擬得到的可燃氣云由于復(fù)雜障礙物分布和環(huán)境風(fēng)影響，爆炸極限范圍內(nèi)的氣云在空間上分布高度不連續(xù)、不均勻，在泄漏口附近空氣夾帶較少，不容易燃燒，在遠離泄漏口的地方存在復(fù)雜湍流流動以及不均勻的初始速度分布，因此其爆炸結(jié)果無法用常見的經(jīng)驗規(guī)律來預(yù)測.使用上述數(shù)值方法開展加氫站內(nèi)大規(guī)模、非規(guī)則、非均勻氣云爆炸模擬.以B5的可燃氣云分布為初始條件，設(shè)置弱點火源位置為(22.0, 31.0, 1.5) m，并設(shè)置一系列壓力監(jiān)測點(21.2, 26.4, 1.5)、(39.6, 41, 1.5)、(25.2, 52.8, 1.5)、(29.1, 32.6, 5.5) m 用于分析各設(shè)備所接收到的超壓.

t=40 ms 時壓力波已傳播至計算域邊界，圖8 為該時刻B5工況下氣云爆炸模擬結(jié)果.工況B5在風(fēng)向上存在較多障礙物，可燃氣云在環(huán)境風(fēng)與復(fù)雜障礙物群的耦合作用下產(chǎn)生了劇烈的湍流，且障礙物附近存在大氣渦旋，使得火焰面形狀變得不規(guī)則，從而增大火焰面積，導(dǎo)致火焰?zhèn)鞑ニ俣燃眲∩仙⑦M一步促進湍流的發(fā)展，同時，當(dāng)量濃度區(qū)域由于化學(xué)反應(yīng)更加劇烈，可在切片中看到相對高溫的“帶狀”區(qū)域[20-21].在40 ms 時，火焰的高溫輻射范圍不再擴展，但火焰面褶皺較多.氣云的尺寸對爆炸強度有顯著影響，氣云初始尺寸越大，火焰加速距離越長，對于前方未燃氣體的壓縮程度增加，最終正向反饋激勵了超壓上升速率，在爆心處可產(chǎn)生約600 kPa 的超壓.40 ms 之后，壓力波傳播至控制室，防火墻的存在顯著減弱了壓力波的危害，有防火墻遮蔽的控制室前比無防火墻時的最大超壓低約70 kPa.考慮使人重傷的溫度為450 K，則火焰產(chǎn)生的危險范圍約為29 m×16 m.考慮使人輕傷、重傷、死亡的超壓為20、50、400 kPa，則壓力波產(chǎn)生的輕傷、重傷和死亡危險范圍約為51 m×51 m、48 m×46 m、14 m×8 m.

圖8 t=40 ms 時氣云爆炸火焰面、壓力波陣面Fig.8 Flame and pressure wave front of gas cloud explosion at t=40 ms

為精確分析壓力波的危害，根據(jù)超壓-沖量破壞準(zhǔn)則對各設(shè)施進行危害評估.圖9 為工況B5測點P1～P4的超壓時程曲線.考慮到B5總體破壞范圍更大，選取了距離爆心更遠的監(jiān)測點，而這組監(jiān)測點滿足隨爆心越遠爆炸超壓越低、超壓變化速率更慢的趨勢.同一高度處，P1最先接收到超壓，由于處于障礙物之間，超壓峰值遠高于其他監(jiān)測點，達到620 kPa，第一次正壓作用時長為17.9 ms，沖量是第二次正壓較弱，主要破壞為一次正壓帶來.由于P2、P3和爆距相差不大，而又距離爆心較遠，氣云初始不規(guī)則分布帶來的差別開始減弱，超壓曲線相似，超壓峰值約為60 kPa，正壓作用時間分別為27.1 ms 和15 ms.P4位于泄漏口正上方，罩棚的約束導(dǎo)致了此處可燃氣云大量聚集，火焰?zhèn)鞑サ皆撎幒蟪瑝杭s為200 kPa，正壓持續(xù)了約22 ms，將在更短時間內(nèi)造成更大范圍內(nèi)的嚴重破壞.根據(jù)超壓-沖量準(zhǔn)則[22]，監(jiān)測點P1～P4對應(yīng)的設(shè)備破壞等級均為6，對于加氫站有著毀滅性的全場破壞.

圖9 工況B5 監(jiān)測點P1～P4 超壓時程曲線Fig.9 Evolution of overpressure at P1～P4 of B5

3 結(jié) 論

本文基于高精度有限差分程序，發(fā)展了基于虛擬點的無反射特征邊界條件方法，實現(xiàn)了大尺度加氫站氫氣泄漏—點火—爆炸完整過程的數(shù)值模擬，并分析了可燃氣云形成以及氣云爆炸的影響因素，完成了加氫站氫氣爆炸事故的定量危害評估.具體結(jié)論如下：

①傳統(tǒng)的Navier-Stokes 特征邊界條件方法考慮邊界上的時間梯度以求解邊界處通量，不適用于使用虛擬點的中心型網(wǎng)格，針對此問題，將一維無黏特征方程在邊界處離散，獲得并使用邊界處的空間梯度對虛擬點進行賦值，構(gòu)造了使用虛擬點的N-S特征邊界條件方法.通過測試二維渦旋和二維火焰面作用于計算域邊界的案例，證實了相對與常用的一階外推方法，改進的方法可以更好地抑制非耗散高精度格式下非反應(yīng)流反應(yīng)流在計算域邊界產(chǎn)生波的虛假反射，降低非物理振蕩.通過對目標(biāo)入流速度或出流壓力進行松弛，可以有效應(yīng)用于通風(fēng)條件的施加.

②針對某真實加氫站在多種通風(fēng)條件下的泄漏擴散過程進行模擬分析發(fā)現(xiàn)：當(dāng)風(fēng)向與無風(fēng)時可燃氣云運動方向一致時，環(huán)境風(fēng)加快可燃氣云在射流軸向上的運動速度，加快可燃氣云的稀釋，抑制徑向上的自然擴散，使可燃氣云變得細長；當(dāng)風(fēng)向與無風(fēng)時可燃氣云運動方向垂直相交時，可燃氣云失速，轉(zhuǎn)變?yōu)檠仫L(fēng)向運動；風(fēng)速越高時，可燃氣云在垂直風(fēng)向的水平擴散受到抑制越大，同時沿風(fēng)向運送加快，兩種作用相互矛盾，時間越長第二種作用更顯著；環(huán)境風(fēng)與復(fù)雜障礙物群耦合形成的空氣流動極大增強了可燃氣云的湍流強度，打亂了濃度的分層分布，擴散范圍將增加至無風(fēng)時兩倍.

③一般認為較高的風(fēng)速可以快速稀釋可燃氣云，但當(dāng)通風(fēng)不足出現(xiàn)點火源時，可燃氣云尺寸仍處于增長階段，尺寸大、高長比大，內(nèi)部又具有復(fù)雜的湍流流動與沿風(fēng)向的速度場，點火后火焰面、壓力波陣面迅速傳播，溫度輻射范圍在短時間內(nèi)突增至可燃氣云初始分布的1/2 區(qū)域，沖擊波在障礙物的作用下反射、疊加，受到循環(huán)激勵效應(yīng)，最高達到600 kPa的超高水平，將對加氫站全場造成毀滅性打擊.