密閉多艙室內柱形氣云爆炸過程二維數值模擬

徐維錚，吳衛國，況 正

(1. 武漢理工大學 高性能艦船技術教育部重點實驗室，湖北 武漢 430063；2. 武漢理工大學 交通學院船舶、海洋與結構工程系，湖北 武漢 430063)

0 引 言

工業可燃性氣云爆炸事故大多是由弱點火點燃可燃氣體引起的，其傳播形式大多是以亞音速傳播的爆燃波。爆燃波的傳播過程很復雜，受環境條件和物理因素的影響極大，當氣體在傳播過程中由于障礙物阻礙發生湍流加強效應或是在初始時刻由高能量直接起爆，其爆炸模式將由爆燃轉為爆轟或直接以爆轟模式發生爆炸。氣體發生爆轟以后將在空氣中形成強度較大的帶有負壓區的空氣沖擊波，對工作人員和周邊結構設施將造成較大的損害。針對空中可燃氣云、云霧爆炸模擬，近年來徐勝利等學者進行過相應的數值研究[1 - 4]，他們將可燃氣云、云霧等效的簡化為一團高壓氣體，而艙室內部可燃氣云爆炸研究較少。

爆炸屬于強間斷問題，其數值模擬對激波捕捉格式提出較嚴格要求。目前，WENO格式作為一種典型的高精度激波捕捉格式，對流場內的激波間斷具有較高的分辨率，適于求解包含激波、膨脹波以及接觸間斷等復雜結構的流場。LIU等[5]于1994年提出了 WENO（weighted essentially non-oscillation scheme）格式，之后，SHU 等[6 - 7]發展了該格式并擴展了其應用。

考慮到WENO格式高精度、穩定性較好的優勢，基于FORTRAN平臺，采用三階WENO有限差分格式，自主開發了密閉多艙室內爆炸波二維數值計算程序。利用所開發的程序開展了多艙室內柱形氣云爆炸波傳播過程及不同艙室內部爆炸載荷特性研究。

1 控制方程

含激波流場采用可壓縮歐拉方程進行描述，其具體形式如下：

在方程（1）的每個方向上均可以看成是一個雙曲守恒律方程：

例如，針對x方向，方程（5）的數值離散形式為：

2 數值方法

2.1 三階 WENO-JS 格式

三階WENO格式的數值離散和推導過程如下[6]：

利用上述2種二階通量的凸組合計算最終具有三階精度的數值通量，即

其中，本文參數ε取值為10-6。光滑因子的表達式如下：

2.2 時間項數值離散

采用三階TVD-RK法對歐拉方程（1）中的時間項進行數值離散，具體的離散格式如下[6]：

3 試驗驗證

為了檢驗所開發數值程序的可靠性和正確性，選取文獻[8]開展的平面激波繞方腔流動試驗進行對比分析。圖1給出計算域的尺寸和網格劃分示意圖（本文計算中，網格尺寸為1 mm）。

在標準狀況下，平面激波馬赫數為1.3，根據沖擊波物理可以確定波后物理參數。圖2給出該算例初始條件。

圖3給出不同時刻數值模擬紋影圖與試驗的對比，圖中不同時刻的波系結構和臺階處繞射形成的漩渦結構均證明了文中程序開發的可靠性和正確性，為進行多艙室內爆炸模擬奠定了基礎。

圖1 計算區域尺寸及網格劃分Fig.1 Size of the computational zone

圖2 激波繞方腔流動初始條件Fig.2 Initial condition of the planar shock wave over a square cavity

圖3 數值與文獻[8]中試驗結果對比Fig.3 Comparisons between numerical and experimental results[8]

4 多艙室內爆炸過程數值計算

本節主要采用驗證的程序對密閉多艙室內平面柱形高壓氣云爆炸過程進行數值模擬研究，主要探討爆炸波傳播路徑及不同艙室內部爆炸載荷特性，為后期的結構抗爆設計及毀傷評估提供一定的參考價值。

4.1 艙室尺寸及測點布置

多艙室有3個艙室組成，分別標記為艙室1、艙室2、艙室3，艙室1與艙室2及艙室3之間通過方形連接導管相連通，每個艙室的具體尺寸參見 圖4。每個艙室的壁面中點處均設置2個壓力測點（分別編號為 NO1，NO2，NO3，NO4，NO5，NO6）對爆炸超壓載荷時間歷程進行輸出（見圖4（a））。壁面邊界條件設置為剛性邊界（不考慮艙室結構的變形），均勻網格尺寸取為1 mm（見圖4（b））。

圖4 多艙室布置及網格劃分Fig.4 Size of multi-cabin and mesh distribution

4.2 初始條件及工況設置

柱形均勻高壓氣云的具體參數：半徑為20 mm，密度為 1.99 kg/m3，壓力為 13.8×105Pa。為了探討爆炸發生位置對爆炸過程的影響規律，本文選取2種典型爆炸工況：工況1表示爆炸發生在艙室1內部，如圖5（a）所示，工況2表示爆炸發生在艙室3內部，如圖5（b）所示。均勻高壓氣云均設置在艙室中間（見圖5中的淺色區域），周圍深色區域表示空氣域，空氣密度為 1.225 kg/m3，壓力為 1.0×105Pa。

4.3 工況一爆炸過程分析

本節主要根據數值模擬結果分析工況1爆炸后爆炸波傳播過程及不同艙室爆炸載荷分布特性。

4.3.1 爆炸波傳播過程分析

為了探討工況1氣云爆炸后爆炸波在多艙室內部的傳播過程，這里給出不同時刻爆炸壓力云圖和數值紋影圖（見圖6）。根據圖6可知工況1多艙室內爆炸過程爆炸波的傳播路徑：爆炸波首先在艙室1中進行自由膨脹運動（圖6（a））爆炸波初次到達艙室壁面時發生正規則反射，壁面處形成局部高壓區域（圖6（b）），之后爆炸波沿著連接導管進入艙室2和艙室3內部（圖6（c）），并緊接著在艙室2和艙室3上下壁面處進行系列正規則斜反射和馬赫反射（圖6（d）和圖6（e）），當爆炸波傳播到艙室2和艙室3的最右端壁面處時再次形成復雜的壁面反射波（圖6（f））。由于多艙室處于封閉空間，因此可以預知爆炸波經過上述在不同艙室之間的傳播過程，爆炸波強度將不斷衰弱，爆炸氣體的動能逐漸轉化為氣體的內能，最終艙室內部將趨近于一個準靜態過程。

4.3.2 爆炸載荷特性分析

為了探討工況1氣云爆炸后不同艙室內部爆炸載荷分布特性，這里給出不同艙室內部壁面測點超壓時間歷程曲線和沖量歷程時間曲線。針對工況1，由于艙室2和艙室3上下對稱，因此，這里只需要討論艙室1和艙室2內部壁面測點的爆炸載荷特性。圖7給出艙室1壁面測點NO1，NO2爆炸載荷時間歷程，從圖7可知艙室內爆炸載荷主要呈現峰值不斷衰減的沖擊波載荷和準靜態壓力。雖然測點NO1，NO2與爆源處于近似相同的距離，但是測點NO1超壓載荷及沖量明顯大于測點NO2，通過分析發現這主要是由于測點NO2靠近艙室2和艙室3，當爆炸波通過連接導管傳播到艙室2和艙室3內部后將明顯減弱測點NO2處的沖擊波強度。

圖5 初始條件及爆炸工況Fig.5 Initial condition and explosion cases

圖6 工況 1 不同時刻爆炸波傳播過程Fig.6 Blast wave propagation process at different time for explosion case one

圖7 艙室 1 壁面測點 NO1，NO2 爆炸載荷時間歷程（工況1）Fig.7 Blast load time histories of gauging point NO1, NO2 in cabin 1 for explosion case one

圖8 給出艙室2壁面測點NO3，NO4爆炸載荷時間歷程，從圖8可知在爆炸初期，測點NO3超壓載荷及沖量稍微大于測點NO4，而在爆炸后期，測點NO3超壓載荷及沖量演變為小于測點NO4，通過分析發現這主要是由于爆炸初期測點NO3相較測點NO4更靠近爆源，因此其接觸到的初始沖擊波強度較大。而在爆炸后期由于測點NO4與連接導管正對，在后續的爆炸波傳播過程中，其將會受到爆炸波的連續直接沖擊，反而造成其超壓載荷、沖量時間歷程稍微高于測點NO3。

圖9給出艙室1，艙室2壁面測點NO1，NO2，NO3，NO4爆炸載荷時間歷程，從圖9可知艙室1內部測點NO1，NO2爆炸載荷強度明顯強于艙室2內部測點NO3，NO4，通過分析發現這主要是由于爆炸發生在艙室1內部，因此其主要能量集中在艙室1內部。同時根據圖9（b）中的沖量歷程曲線還得知，不同艙室內部最終的準靜態壓力相同，如果將沖量積分曲線開始變為直線段的轉折點作為準靜態壓力開始形成的時間，則在該爆炸工況下不同艙室內部測點到達準靜態壓力的時間基本相同（近似為0.5 ms）。

圖8 艙室 2 壁面測點 NO3，NO4 爆炸載荷時間歷程（工況1）Fig.8 Blast load time histories of gauging point NO3, NO4 in cabin 2 for explosion case one

圖9 艙室 1 及艙室 2 壁面測點 NO1，NO2，NO3，NO4 爆炸載荷時間歷程（工況1）Fig.9 Blast load time histories of gauging point NO1, NO2, NO3,NO4 in cabin 1, 2 for explosion case one

圖10 工況 2 不同時刻爆炸波傳播過程Fig.10 Blast wave propagation process at different time for explosion case two

4.4 工況 2 爆炸過程分析

本節主要根據數值模擬結果分析工況2爆炸后爆炸波傳播過程及不同艙室爆炸載荷分布特性。

4.4.1 爆炸波傳播過程分析

為了探討工況2氣云爆炸后爆炸波在多艙室內部的傳播過程，這里給出不同時刻爆炸壓力云圖和數值紋影圖（見圖10）。根據圖10可知工況2多艙室內爆炸過程爆炸波的傳播路徑：爆炸波首先在長方形艙室3中進行自由膨脹運動，初次到達艙室上下壁面時發生正規則反射，壁面處形成局部高壓區域（圖10（a）），之后爆炸波沿著連接導管進入艙室1內部并以球形方式進行擴張（圖10（b）和圖10（c）），在擴張過程中由于艙室1上壁面的限制，形成局部馬赫反射波（圖10（d）），緊接著，爆炸波傳播到艙室2內部并形成復雜的壁面反射波（圖10（e）和圖10（f））。

圖11 艙室 1 壁面測點 NO1，NO2 爆炸載荷時間歷程（工況2）Fig.11 Blast load time histories of gauging point NO1, NO2 in cabin 1 for explosion case two

圖12 艙室 2 壁面測點 NO3，NO4 爆炸載荷時間歷程（工況2）Fig.12 Blast load time histories of gauging point NO3, NO4 in cabin 2 for explosion case two

4.4.2 爆炸載荷特性分析

為了探討工況2氣云爆炸后不同艙室內部爆炸載荷分布特性，這里給出不同艙室內部壁面測點超壓時間歷程曲線和沖量歷程時間曲線。圖11給出艙室1壁面測點NO1，NO2爆炸載荷時間歷程，從圖11可知測點NO1超壓載荷及沖量稍微大于測點NO2，通過分析發現這主要是由于測點NO1靠近艙室3的連接導管，當爆炸波通過連接導管傳播到艙室1內部時，測點NO1將遭受到更強的爆炸波沖擊。

圖12給出艙室2壁面測點NO3，NO4爆炸載荷時間歷程，從圖12可知測點NO3，NO4雖然超壓載荷時間歷程差異較大，但是沖量載荷歷程基本相同，通過分析發現這主要是由于測點NO3，NO4位于艙室2內部，相較于艙室1中的爆源處于遠場區域，當爆炸波通過連接導管擴張到艙室1再傳播到艙室2內部后，由于艙室2容積也較艙室1偏小，傳入的爆炸波沖擊很快使得艙室2內部達到準靜態過程。

圖13給出艙室3壁面測點NO5，NO6爆炸載荷時間歷程，從圖13可知測點NO6初始超壓載荷明顯強于測點NO5，后期強度衰減速度快于測點NO5，這主要是由于艙室3為長方形艙室，測點6距離爆源較近，因此初期沖擊波載荷強度大，由于測點6距離連接導管近，因此其后期超壓衰減速率更快。艙室3容積較小，雖然測點NO6距離連接導管較近，對該處的沖量歷程有一定的影響，但是測點6距離爆源較近，最終使得NO5，NO6沖量歷程曲線基本重合。

圖14給出艙室1～艙室3壁面測點NO1～NO6爆炸載荷時間歷程，從圖14可知艙室3內部測點NO5，NO6爆炸載荷強度明顯強于艙室1內部測點NO1，NO2，通過分析發現這主要是由于爆炸發生在艙室3內部，因此其主要能量集中在艙室3內部。艙室1內部測點NO1，NO2爆炸載荷強度稍微強于艙室2內部測點NO3，NO4，通過分析發現這主要是由于NO1，NO2相較測點NO3，NO4距離爆炸源更近一些，沖擊波強度隨著傳播距離的增大，強度逐漸衰弱。同時根據圖14（b）中的沖量歷程曲線還得知，盡管爆炸波達到各測點的時間存在差異性，但是該爆炸工況下不同艙室內部測點到達準靜態壓力的時間基本相同（近似為0.5 ms），這一結論與4.3.2小節中的分析結果相同。因此，通過這2個工況的爆炸分析，可以初步得出如下結論：對于特定的封閉艙室，準靜態壓力形成時間主要與艙室容積和爆源能量相關。其實，從準靜態壓力形成的物理過程來分析也很容易闡述這種關系，由于準靜態壓力全場均勻，因此其在全場空間內必然應同時到達。

圖13 艙室 3 壁面測點 NO5，NO6 爆炸載荷時間歷程（工況2）Fig.13 Blast load time histories of gauging point NO5, NO6 in cabin 3 for explosion case two

圖14 艙室 1～艙室 3 壁面測點 NO1，NO2，NO3，NO4，NO5，NO6爆炸載荷時間歷程（工況2）Fig.14 Blast load time histories of gauging point NO1, NO2, NO3,NO4, NO5, NO6 in cabin 1,2,3 for explosion case two

圖15 工況1及工況2壁面測點NO1爆炸載荷時間歷程對比圖Fig.15 Comparisons of blast load time histories of gauging point NO1between explosion case one and two

4.5 工況 1 與工況 2 對比分析

為了進一步分析爆源位置對多艙室內部不同測點處爆炸載荷特性的影響規律，圖15～圖20分別給出2種工況下測點NO1～NO6的超壓載荷和沖量時間歷程曲線圖。

根據圖15和圖16可知，針對艙室1內部測點NO1，NO2，工況1相較工況2具有更強的超壓載荷和更大的沖量，這主要是由于工況1中爆炸發生在艙室1中。

根據圖17和圖18可知，針對艙室2內部測點NO3，NO4，2種工況超壓載荷歷程差異較大，但是沖量載荷歷程差別較小，這主要是由于艙室2中的測點相較2種工況中的爆源距離較大，屬于遠場的范圍，沖量的歷程主要由準靜態壓力決定。

根據圖19和圖20可知，針對艙室3內部測點NO5，NO6，工況2相較工況1具有更強的超壓載荷和更大的沖量，這主要是由于工況2中爆炸發生在艙室3中。

圖16 工況1及工況2壁面測點NO2爆炸載荷時間歷程對比圖Fig.16 Comparisons of blast load time histories of gauging point NO2 between explosion case one and two

圖17 工況1及工況2壁面測點NO3爆炸載荷時間歷程對比圖Fig.17 Comparisons of blast load time histories of gauging point NO3 between explosion case one and two

圖18 工況1及工況2壁面測點NO4爆炸載荷時間歷程對比圖Fig.18 Comparisons of blast load time histories of gauging point NO4 between explosion case one and two

圖19 工況1及工況2壁面測點NO5爆炸載荷時間歷程對比圖Fig.19 Comparisons of blast load time histories of gauging point NO5 between explosion case one and two

圖20 工況1及工況2壁面測點NO5爆炸載荷時間歷程對比圖Fig.20 Comparisons of blast load time histories of gauging point NO5 between explosion case one and two

5 結 語

采用三階WENO有限差分格式，自主開發了密閉多艙室內爆炸波二維數值計算程序。利用所開發的程序開展了多艙室內柱形氣云爆炸波傳播過程及不同艙室內部爆炸載荷特性研究。通過本文的研究主要得到以下幾點結論：

1）多艙室內爆炸過程中，爆炸所在艙室內部爆炸載荷強度明顯強于相鄰艙室，尤其需要重點關注爆炸所在艙室通過連接導管對直接相鄰艙室爆炸載荷的影響；爆源位置對不同艙室內部超壓載荷時間歷程影響較大，而對處于遠場的艙室內沖量時間歷程影響較小。

2）對于特定的封閉艙室，準靜態壓力峰值及準靜態壓力形成時間均全場均勻，準靜態壓力形成時間主要與艙室容積和爆源能量相關。

3）本文研究可為工程抗爆結構設計及爆炸毀傷評估提供一定的參考和指導。