氣體泄爆壓力分步計算模型及其湍流修正*

孫 松，高康華，邱艷宇，王明洋

（陸軍工程大學爆炸沖擊防災減災國家重點實驗室，江蘇 南京 210007）

可燃性氣體因其來源廣、熱能高、污染小等特點已經成為人們生產、生活中最常用的能源之一；與此同時，頻繁發生的氣體爆炸事故也給人們的生命財產安全帶來了巨大威脅。泄爆是工程實際中最常用的防護措施，研究者們對于泄爆條件下空間內的爆炸壓力進行了大量的實驗研究[1-5]，對氣體爆炸的物理過程及各影響因素有了較全面的了解。

同時，出于安全性、經濟性和實用性的考慮，研究者們以大量實驗為基礎結合氣體爆炸物理過程進行了理論模型的研究，以期能夠預測空間內氣體爆炸壓力的變化情況。Bradley 等[6]根據守恒方程、理想氣體狀態方程和等熵壓縮方程推導出關于升壓速率的偏微分方程組，成為其后理論模型發展的基礎。Han 等[7]通過對比模型計算值與實驗結果發現，壓力計算模型可以較好地描述建筑空間內的壓力變化，隨后通過模型計算探討了泄放面積、泄放壓力、屋室尺寸、氣體濃度和點火位置對爆燃壓力的影響。Ugarte 等[8]根據不同點火位置、容器形狀在一維等熵泄放模型的基礎上對火焰面積、火焰傳播速度和外部爆燃壓力進行修正，對比發現修正后的模型與實驗結果吻合更好。美國防火規范NFPA68-2013[9]中總結大量實驗數據得到不同工況下泄爆壓力及泄爆面積的簡化確定方法，為實際工程防護提供了重要指導。Sustek 等[10]使用常用經驗與半經驗公式對不同實驗工況進行計算，對比發現，由于各經驗公式中的參數很多都是基于實驗數據擬合得到，因此各公式適用范圍較窄，使用經驗公式計算不同工況可能會產生較大誤差。

現有的氣體爆燃壓力計算方法多是以大量實驗數據為基礎評估空間內的最大峰值超壓，無法準確描述空間內的壓力變化過程；數值模擬方法計算成本高，耗時長，不利于工程快速評估應用；而常用的壓力計算模型[11-13]多為關于升壓速率的偏微分方程組，其求解過程較復雜，同時這些模型大多假設泄放構件開啟后先泄放已燃氣體，再泄放未反應氣體，無法計算同時泄放已燃氣體和未反應氣體的工況。本文中，通過將爆燃過程進行微分，假設每個微小時間段內爆燃泄放過程均按照先燃燒、再泄放、最后壓力平衡的過程獨立分步進行，得到了泄爆壓力分步計算模型；同時在尺寸為2 m×1.2 m×0.6 m 的爆炸腔體一端安裝不同面積的泄壓結構進行泄爆實驗，對分步壓力計算模型進行驗證，并根據實驗結果對模型進行湍流修正。該模型清晰地反映了氣體爆燃泄放機理，通過對微分時段內的分步計算過程進行迭代，可以較好地描述泄爆條件下空間內的壓力及各狀態參數的變化過程，計算過程簡單，可為工程實際中快速評估內爆炸壓力提供參考；經過湍流加速修正的壓力計算模型可以較好地描述因泄放截面改變引起的湍流區域內的壓力變化情況。

1 理論模型

1.1 基本假設

由于爆燃過程非常復雜，涉及到燃燒學、熱動力學、流體動力學等多方面內容，因此在模型計算時假設可燃氣體密度均勻充滿密閉空間且符合理想氣體狀態方程；同時假設密閉容器為絕熱壁面，忽略熱交換及內部反應的能量耗散過程。

泄放條件下氣體爆燃的超壓變化過程受燃燒產生的壓力和泄放降低的壓力共同影響，變化過程較復雜，為簡化計算將燃燒過程和泄放過程單獨考慮，將氣體爆燃泄放過程分成3 個獨立階段：（1）只燃燒不泄放(confined combustion stage)。將此階段簡化為密閉空間內的爆燃過程，不考慮泄放對于爆燃壓力的影響。（2）泄放過程(venting stage)。第1 階段燃燒結束后室內壓力上升，高于外部環境壓力，此時燃燒產物和未燃氣體將分別發生泄放。（3）泄放后壓力平衡階段(pressure equilibrium stage)。該階段泄放過程結束，此時空間內部燃燒產物與未燃氣體壓力不同，空間內部會發生氣體流動最終達到壓力平衡。密閉空間內氣體爆燃泄放過程如圖1 所示，圖1(a)～(c)為燃燒反應階段，圖1(d)～(e)為泄放階段，圖1(f)為壓力平衡階段。

圖 1 爆燃泄放分步計算模型示意圖Fig.1 Diagrams of the sub-steps model for combustion venting

1.2 控制方程

1.2.1 燃燒階段

通常假設預混氣體以規則球形燃燒并向四周擴散，則在 Δt時間內所燃燒的體積為：

式中： VB,j為第j 個時間段內燃燒的氣體體積；S 為火焰傳播速度； Vh,j-1為j-1 時刻已燃氣體體積。

現定義燃盡率 Δμj為 Δt時間內燃燒的氣體質量與可燃氣體總質量之比，則：

式中： ρc,j-1為 j-1 時刻未燃氣體密度， V0為 爆炸空間總體積， ρ0為初始氣體密度。

又因為根據等熵壓縮方程，在絕熱過程中有：

式中： γu為未燃氣體絕熱指數， pj-1為j-1 時刻空間內的壓力，p0為空間內的初始壓力。

則 Δμj可以寫成：

由于空間內的壓力隨氣體燃燒分數的增大而升高，假設室內壓力的升高與燃燒產物生成量成線性關系[14]，則圖1 中(c)階段室內壓力為：

式中： pb為 密閉空間內氣體爆燃產生的最大壓力，該值與氣體種類和濃度有關； pn,j為j 時刻密閉燃燒階段結束時內部壓力值。

由于整個反應過程為絕熱過程且壓力 pn,j已知，則根據等熵壓縮方程：

可得圖1 中j 時刻(c)階段未燃氣體體積為：

已燃氣體體積為：

式中： Vc,j-1為 j-1 時刻未燃氣體體積； Vh,j-1為j-1 時刻已燃氣體的體積。

由質量守恒可得圖1 中(c)階段未燃氣體密度為：

式中： ρn,c,j為j 時刻燃燒階段結束時未燃氣體密度, ρc,j-1為j-1 時刻未燃氣體密度。

已燃氣體密度為：

式中： ρn,h,j為j 時刻燃燒階段結束時燃燒產物密度, ρh,j-1為j-1 時刻燃燒產物密度。

此時第j 個時間段內燃燒階段結束，開始進行泄放階段的壓力計算。

1.2.2 泄放階段

考慮實際情況，假設該階段未燃氣體與燃燒產物分別發生泄放且外界氣壓恒為標準大氣壓，則以次臨界速度（ p0/ pj-1≥βk）泄放時[6,14]，單位時間內泄放的燃燒產物和未燃氣體體積分別為：

以臨界速度（ p0/ pj-1≤βk）泄放時，單位時間內泄放的燃燒產物和未燃氣體體積分別為：

式中： Ah,j和 Ac,j分別為 Δt 時間內燃燒產物和未燃產物的泄放流動面積； γb為 燃燒產物的絕熱指數； βk為臨界比值，其計算公式為 βk=[2/(γu+1)]γu/(γu-1)。

由于發生氣體泄放，未燃氣體和燃燒產物的密度及壓力均發生變化，根據等熵壓縮方程可得圖1 中(e)階段j 時刻泄放階段結束時空間內未燃氣體和燃燒產物的壓力分別為：

j 時刻泄放階段結束時空間內未燃氣體和燃燒產物的密度分別為：

在第2 階段中假設只發生氣體的泄放流動，但空間內未燃氣體與已燃氣體所占體積不發生改變，此時第2 階段結束，開始進行空間內壓力平衡階段的計算。

1.2.3 壓力平衡階段

由等熵壓縮方程可知:

將式(19)代入式(20)并變形后得完整反應過程結束時空間內的壓力：

式中： β= pr,j/pj。

求得圖1 中(f)階段的壓力pj后，根據狀態方程與等熵壓縮方程可以求得(f)階段的氣體物理參數：

式中： ρh,j和 ρc,j分別為j 時刻整個爆燃泄放過程結束后空間內未燃氣體和燃燒產物的密度； Vh,j和 Vc,j分別為此時空間內未燃氣體和燃燒產物的體積。

此時第3 階段結束，一個微分時間段內腔體中氣體壓力、密度、體積的完整變化過程結束，開始迭代進行下一個時間步的計算。

反應后期若空間內混合氣體已全部燃盡，則燃燒過程結束，此時氣體繼續由內部流向外部空間，直至內外壓力平衡。該階段在迭代計算過程中省略圖1 中(b)、(c)過程，即該情況下完整反應流程圖變為圖1(a)、(d)、(e)、(f)。

2 實驗驗證

2.1 實驗裝置

實驗裝置由配氣系統、長方體爆炸容器、點火系統、泄爆系統與數據采集系統構成，如圖2 所示。

圖 2 實驗系統示意圖Fig.2 Schematic of the experimental system

爆炸腔體為自行設計的長方體容器，長、寬、高分別為2、1.2、0.6 m，最大可承受1 MPa 的壓力。爆炸容器一端為開敞的用于安裝泄爆結構的卡箍式法蘭盤，其尺寸為1.2 m×0.6 m，泄爆構件示意圖如圖3所示；另一端壁面中部通過法蘭結構與電阻絲點火器連接，構成點火系統，其最大點火溫度可達500 ℃。實驗中選用乙烯作為可燃氣體，通過在爆炸腔體側面等距離設置進氣孔和抽氣孔運用壓力分配法配氣，氣體濃度由流量計和壓力表控制。進氣孔內配置有細長進氣管道，管道上均勻分布有小孔，這種均布直管式進氣系統可以保證腔體內預混氣體混合較均勻。在爆炸腔體頂部裝有PCB113B26 系列壓電式高頻壓力傳感器，分別距泄放端口1.2、0.1 m。數據采集系統將采集到的電信號轉化為壓力信號并輸出，其數據采樣頻率為200 kHz。高速攝影儀置于距泄爆口側面8 m 處，拍照頻率為1 000 s-1。

圖 3 泄爆構件示意圖(單位為mm)Fig.3 Schematic diagrams of vent covers (unit in mm)

如圖3 所示，實驗中泄爆構件采用擊穿壓力均為7.5 kPa 的泄爆板和泄爆膜。泄爆板泄放面積為1.1 m×0.5 m，泄爆膜泄放面積為0.3 m×0.6 m。采用泄爆膜作為泄爆構件開展實驗時，在腔體泄爆端安裝尺寸為1.2 m×0.6 m 的隔板，并在隔板中部開設面積為0.6 m×0.3 m 的泄爆口，使用泄爆膜對其進行封閉，如圖3(a)所示；采用泄爆板作為泄爆構件開展實驗時，直接在腔體泄爆端安裝尺寸為1.2 m×0.6 m 的泄爆板，如圖3(b)所示。構件擊穿壓力通過在泄爆結構背面粘貼應變片以應變信號消失時刻腔體內對應的壓力確定，實驗時選取的可燃氣體的體積分數為4%～8%。

2.2 結果對比

使用分步壓力計算模型對實驗工況進行計算，由于實驗工況為矩形容器端部中心點火，因此需針對火焰的傳播形狀對計算模型進行修正。由高速攝影觀察到火焰由泄爆端流出時其陣面可近似認為是平面，因此為計算簡便假設火焰未接觸壁面時為球形傳播，接觸壁面時為平面傳播，則端部點火時火焰的體積表達式為：

式中：2a、2c 為爆炸腔體的長和高，x 為火焰半徑。

泄爆板和泄爆膜工況的計算值及實驗測量結果分別如圖4～5 所示。由圖4 可知：采用泄爆板作為泄爆構件時，2 個傳感器測得的壓力時程曲線基本重合，符合準靜態壓力形式[15]；此時腔體內部為單峰值壓力曲線，分析認為當構件擊穿后由于泄放面積較大使得短時間內大量氣體外泄，引起腔體內部壓力下降至環境壓力。該情況下模型計算值與實驗結果吻合較好。

圖 4 可燃氣體的體積分數不同、采用泄爆板作為泄爆構件時，泄爆壓力的計算值與實驗測量結果的比較Fig.4 Comparison of gas explosion venting pressure between calculated results by the model and experimental ones in the experiments with vent plates at different volume fractions of combustible gas

圖 5 可燃氣體的體積分數不同、采用泄爆膜作為泄爆構件時，泄爆壓力的計算值與實驗測量結果的比較Fig.5 Comparison of gas explosion venting pressures between the model and the tests with vent films at different volume fractions of combustible gas

觀察圖5 發現，當采用泄爆膜作為泄爆構件時腔體內壓力曲線均為雙峰值曲線，分析認為由于泄爆膜面積較小，構件擊穿后腔體內部氣體外流引起壓力下降產生第1 個壓力峰值；氣體外泄同時又會加劇腔體內部湍流擾動，扭曲火焰陣面，提高燃燒速率與升壓速率，引起腔體內壓力再次上升產生第2 個峰值[1]。同時由圖5 發現當乙烯體積分數較高時靠近泄爆口的壓力傳感器T1 測得的第2 個峰值壓力要大于內部傳感器T2 相應的測量值，分析原因認為由于泄爆膜開口面積較小，外流氣體流動面積的突然改變、結合容器壁的限制作用使得近泄爆口處湍流效應明顯，可燃氣體得到充分燃燒，腔體內部產生壓力梯度，泄爆口處壓力大于腔體內部壓力，且越接近最佳體積分數該現象越明顯[16-17]。此時理論模型可以較好地描述腔體內部的壓力變化情況，但與近泄爆口處壓力變化情況存在誤差。

為了進一步驗證計算模型的適用性，運用該模型計算其他研究者的實驗工況[5,7,18]，其結果如圖6所示。通過對比發現計算結果與實驗結果吻合較好，該方法具有較好的普適性。

圖 6 運用本文計算模型計算其他研究者的實驗工況[5, 7, 18]下的瀉爆壓力Fig.6 Gas explosion venting pressures calculated by the model proposed in this paper for the experimental conditions[5, 7, 18] of other researchers

3 模型修正

由于傳感器T1 測得的第2 個峰值壓力大于腔體內部相應的壓力值是因泄放截面改變引起的火焰湍流傳播引起，因此在此引入湍流加速因數λ 對火焰傳播速度進行修正，以計算靠近泄放口位置的壓力情況。湍流加速作用因泄放截面的改變引起，因此將傳播路徑中的障礙物截面與火焰傳播截面的比值，即阻塞率, 作為其增強因子的影響因素。研究者們對障礙物對火焰傳播的加速作用進行了大量的實驗研究[16,19-20]，雖然實驗環境不同，但是障礙物阻塞率對于火焰的湍流加速效果表現出一定的規律性。現根據美國防火規范NFPA68 確定因截面改變引起的湍流加速因數[9]：

式中：Aobs為障礙物截面面積，As為火焰傳播路徑截面面積， λ0為因自湍流等因素引起的湍流因數。在本文的泄爆膜實驗工況中Aobs取0.37， As取0.55， λ0取2。

則未因傳播截面改變引起湍流加速的火焰傳播速度為[21]:

近泄爆口處因截面改變產生湍流加速的火焰傳播速度為：

式中： S0為基本火焰傳播速度；p 和p0分別為空間內的壓力和初始壓力；T0和Tu分別為初始溫度和未燃氣體溫度； α 和 η為 根據實驗確定的參數，此處 α取2，η取-0.15。

現將修正后的火焰傳播速度代入壓力計算模型，則計算結果如圖7 所示。由圖7 可知，經湍流修正后的壓力計算模型得到的峰值壓力及升壓時間與傳感器T1 測得的實驗結果吻合較好，使用湍流修正因數預測因傳播截面改變引起的湍流區域的壓力值具有可行性。

圖 7 針對體積分數不同的可燃氣體，修正后的壓力計算模型得到的瀉爆壓力與實驗結果的比較Fig.7 Comparison of gas explosion venting pressures calculated by the modified model and the test data at different volume fractions of combustible gas

4 結 論

（1）通過將爆燃泄放過程進行微分，假設每個微小時間段內爆燃泄放過程均按照先燃燒、再泄放、最后壓力平衡的過程獨立分步進行，得到了泄爆壓力分步計算模型。該模型計算簡單，不僅可以計算爆燃泄放過程中的壓力變化情況，還可以得到各時刻空間內的氣體密度、體積等物理狀態參數。

（2）在尺寸為2 m×1.2 m×0.6 m 的長方體容器中充滿乙烯-空氣混合氣體進行實驗驗證，發現當采用大面積泄爆板作為泄爆構件時爆燃壓力曲線均為單峰值曲線，該條件下壓力模型計算結果與實驗結果吻合較好。

（3）當采用小面積泄爆膜作為泄爆構件時會因外流氣體流動截面改變引起火焰湍流加速，使得腔體內產生壓力梯度，近泄爆口處壓力值大于腔體內部。通過湍流加速因子對火焰速度進行修正，則修正后的壓力計算模型可以較好地描述近泄爆口處的壓力變化情況。