劉 沖，杜 揚，梁建軍，張培理，孟 紅

（1. 陸軍勤務學院油料系，重慶 401311；2. 陸軍勤務學院勤務指揮系，重慶 401311）

汽油作為一種工業(yè)制成品，廣泛應用于交通運輸和工業(yè)清洗行業(yè)。但是，作為一種易燃易爆易揮發(fā)的危險品，在儲存、運輸及加注的全過程中，汽油時刻有著安全威脅。為有效應對安全威脅，降低人類的生命和財產(chǎn)損失，非常有必要針對汽油蒸氣/空氣混合氣的爆炸特征開展深入研究。此前，已有大量針對不同可燃氣體爆炸特征的研究，發(fā)現(xiàn)受限空間結構和形狀特征對可燃氣爆炸超壓和火焰的變化規(guī)律有較大影響。主要影響因素包括為狹長受限空間長徑比[1-2]、半受限空間開口尺寸及位置[3-5]、受限空間障礙物形狀及尺寸[6-8]和所含彎管及分支結構特征[9-11]等。特別地，彎管及分支結構是儲油洞庫坑道及輸油管路中常見的結構形式，研究表明彎管及分支結構的存在對爆炸超壓和火焰?zhèn)鞑ニ俣扔休^強的激勵作用，并在一定程度上可與擋板式障礙物類比[9]。Zhu 等[12]發(fā)現(xiàn)含彎管結構受限空間可燃氣爆炸規(guī)律主要受彎管結構數(shù)量影響，而受彎管方向影響較小。Niu 等[13]對甲烷在輸氣管網(wǎng)中的爆炸特征進行了研究，發(fā)現(xiàn)爆炸超壓峰值及火焰?zhèn)鞑ニ俣仁苓B接平行管網(wǎng)的分支長度影響，兩分支長度不同時爆炸更為劇烈復雜。

盡管彎管和分支結構在形狀特征上有一定相似，但由于分支結構具有更復雜的通道和多個流動終點[14]，因此確有必要對含多分支結構受限空間油氣爆炸特征進行研究。Zhang 等[14]對5 種不同分支結構受限空間爆炸進行研究，發(fā)現(xiàn)油氣爆燃指數(shù)的增加與分支數(shù)量呈函數(shù)關系。Li 等[15]、杜揚等[16]通過實驗研究，發(fā)現(xiàn)分支結構可以增大汽油蒸氣/空氣混合氣爆炸超壓峰值、升壓速率及火焰?zhèn)鞑ニ俣?。杜揚等[17]對含雙側分支受限空間油氣爆炸進行研究，發(fā)現(xiàn)火焰在含有雙側分支結構的管道內(nèi)會呈現(xiàn)“浪花狀”火焰。目前，人們已經(jīng)對含分支結構受限空間油氣爆炸進行了一定的研究，但研究多集中于實驗分析和現(xiàn)象描述，對爆炸流場和火焰的耦合關系仍需深入研究。

本文中，在已有實驗研究[17]的基礎上，通過大渦模擬方法對含雙側分支受限空間油氣爆炸過程進行模擬，并將火焰形態(tài)、傳播速度和爆炸動態(tài)超壓等模擬結果與實驗進行驗證，分析含雙側分支受限空間油氣爆炸流場結特征、火焰?zhèn)鞑バ袨楹蛣討B(tài)超壓特性，指出“浪花狀”火焰形成的內(nèi)在原因，揭示流場、火焰和超壓的耦合機理。

1 數(shù)學模型

1.1 大渦模擬控制方程

大渦模擬通過對小尺度湍流脈動進行過濾，求解計算量比直接模擬方法小，但空間分辨率比雷諾平均方法高，適合于捕捉油氣爆炸的湍流特征[18]。

大渦模擬對小尺度湍流脈動進行Favre 過濾后，利用N-S 方程對大尺度湍流脈動進行直接求解。經(jīng)Favre 過濾后，大渦模擬控制方程的連續(xù)性方程、動量方程、能量方程、反應進程變量及其守恒方程為：

式中：ρ 為密度，t 為時間，ui、uj為速度分量，xi、xj為方向分量，τij為黏性應力張量分量，p 為壓力，hs為顯焓，λ 為熱導率，T 為溫度，D 為擴散系數(shù)，為歸一化化學反應速率。Yf為局部未燃燃料質(zhì)量分數(shù)，Yf0為初始未燃混合氣中燃料質(zhì)量分數(shù)；c 為反應進程變量，c=0 代表化學反應未進行，c=1 代表化學反應完全；符號～表示Favre 過濾量，—表示物理空間過濾量。

為使大渦模擬控制方程封閉，需建立亞網(wǎng)格模型。壁面自適應局部渦黏模型（WALE）能夠較好地捕捉流場從層流到湍流的轉(zhuǎn)變，且計算量較小[19]。因此，采用WALE 模型作為計算的亞網(wǎng)格模型。渦黏模型方程為：

式中：δij為單位張量分量；μS為亞網(wǎng)格黏性系數(shù)，為：

1.2 亞網(wǎng)格燃燒模型

湍流燃燒過程中，火焰鋒面厚度較小，為取得精確的結果，必須采用更精細的網(wǎng)格。但這樣會使網(wǎng)格量增大，運算成本陡增。因此，采用Zimont 模型[20]，對火焰鋒面進行增厚處理，確保在有限運算量下精確模擬火焰鋒面。在Zimont 模型中，湍流火焰速度公式為：

2 數(shù)值模擬

2.1 物理模型

圖1 實驗系統(tǒng)Fig. 1 Experimental system

以文獻[6]中的實驗數(shù)據(jù)為基礎開展數(shù)值模擬。其實驗系統(tǒng)如圖1 所示，由實驗管道、點火系統(tǒng)、壓力采集系統(tǒng)、高速攝影系統(tǒng)、油氣產(chǎn)生系統(tǒng)、油氣濃度測試系統(tǒng)、同步觸發(fā)系統(tǒng)和數(shù)據(jù)采集系統(tǒng)等部分組成。圖2 為實驗管道尺寸，數(shù)值模擬幾何模型與實驗管道尺寸一致。主管道及兩側分支管道內(nèi)部橫截面尺寸均為100 mm×100 mm，主管道左端面密閉，中心處為點火點，右端面開口率為50%。實驗前，用聚乙烯薄膜密封右端開口，但爆炸破膜壓力很小，對爆炸發(fā)展演變影響較小，因此在數(shù)值計算過程中未考慮薄膜的影響[8]。為保證計算結果的準確性，數(shù)值計算中在管道右端開口處設置了一個700 mm×500 mm×800 mm 的外部空間模擬實驗外場，以抑制爆炸回流對管道內(nèi)部的影響，如圖2 所示。沿爆炸傳播方向，在分支前、后各275 mm 的主管道處分別設置兩個壓阻型傳感器P1、P2，對應數(shù)值計算相同位置兩個壓力監(jiān)測點。實驗時，管道充裝汽油體積分數(shù)1.73%的汽油-空氣混合氣，環(huán)境壓力約為1.013×105Pa，環(huán)境溫度約為300 K。高速攝影儀記錄爆炸過程火焰發(fā)展過程，幀率為500 Hz。

2.2 網(wǎng)格劃分

圖2 實驗管道Fig. 2 Experimental pipe

大渦模擬計算過程中，對較大尺寸渦團直接求解、對較小尺寸渦團過濾的特性決定了利用大渦模擬進行數(shù)值計算沒有“網(wǎng)格無關性”，越小的網(wǎng)格尺寸和時間步長會獲得越精確的計算結果，這與雷諾平均方法不同。但隨著網(wǎng)格尺寸趨近更小，計算量將趨近于直接模擬方法，因此在網(wǎng)格劃分上要盡可能平衡計算精度和計算成本[21]。經(jīng)過調(diào)試驗證后，本文中網(wǎng)格劃分如圖3 所示，圖3(a)為計算域整體網(wǎng)格劃分，圖3(b)為分支處局部加密網(wǎng)格。因計算區(qū)域規(guī)則，內(nèi)場和外場網(wǎng)格均為六面體結構網(wǎng)格。在內(nèi)場分支結構區(qū)域加密網(wǎng)格，最小為2 mm；外場沿爆炸傳播方向網(wǎng)格逐漸增大，最大為10 mm：計算域網(wǎng)格總數(shù)約為490 萬個。

圖3 網(wǎng)格劃分Fig. 3 Computational grids

2.3 邊界條件和初始條件

為與實驗條件相對應，對含雙側分支結構受限空間充裝汽油體積分數(shù)1.73%的汽油-空氣混合氣的爆炸過程進行模擬。主管道和分支管道壁面設置為絕熱無滑移壁面且溫度為300 K，外場空間壁面設置為壓力出口且壓力為1.013×105Pa。環(huán)境溫度為300 K，流場初始速度和反應進程變量為零。在實驗管道左端面中心處設置一個半徑5 mm 的半球形區(qū)域，將其內(nèi)部反應進程變量設為1，用以模擬實驗點火過程。

2.4 計算方法

計算過程中，采用二階迎風格式對對流項進行離散，采用二階中心差分格式對擴散項進行離散，采用SIMPLE 算法求解壓力和速度耦合方程組。時間步長為10-4s，每時間步迭代20 次，保證能量殘差小于10-6，動量殘差小于10-5，質(zhì)量殘差小于10-6，反應進程變量殘差小于10-4。所有計算在Think Station P700 工作站完成，耗時約100 h。

3 結果分析

油氣爆炸關鍵特征通常包括火焰?zhèn)鞑バ袨楹统瑝鹤兓?guī)律，為驗證數(shù)值模擬結果的準確性，選取火焰?zhèn)鞑バ螒B(tài)、火焰鋒面?zhèn)鞑ニ俣群统瑝呵€三方面進行數(shù)值模擬，并和文獻[17]實驗數(shù)據(jù)對比。

3.1 火焰?zhèn)鞑バ袨?/h3>

圖4 為火焰形態(tài)的實驗和數(shù)值模擬結果。在初始階段（32 ms 前），火焰受到管道側壁的約束沿軸向拉伸，以典型的“指形”火焰?zhèn)鞑?，且此時火焰鋒面未受到分支結構的擾動作用，火焰鋒面包絡面清晰光滑，為層流狀。當火焰?zhèn)鞑ブ练种Ч艿栏浇鼤r，因壁面約束消失和分支管道尖角的障礙物作用，致使流場波動劇烈，火焰從層流狀向湍流狀轉(zhuǎn)捩，火焰鋒面褶皺彎曲?；鹧?zhèn)鞑ブ练种Ч艿溃?4～64 ms）時，火焰鋒面形成了三路傳播路徑，即火焰繼續(xù)沿主管道向末端傳播，同時在兩側分支管道中對稱傳播。在分支管道中，火焰鋒面因受到右側壁面的約束，在傳播中逐漸向左側壁面彎曲，并形成了實驗中對稱的“浪花狀”火焰。并且，隨著反應的不斷進行，狹長火焰鋒面兩側大量未燃氣被同時引燃，火焰鋒面嚴重變形并形成明顯褶皺。

圖4 實驗和數(shù)值模擬的火焰形態(tài)Fig. 4 Experimental and simulated flame structures

由圖4 可以看出，火焰在主管道中傳播時，數(shù)值模擬能夠直觀準確地反映實驗中的火焰形態(tài)變化。當火焰?zhèn)鞑ミM入分支管道時，數(shù)值模擬能夠較準確地反映實驗中“浪花狀”火焰的輪廓和鋒面位置，但在細節(jié)上和實驗仍存在一定差別。究其原因有兩點：一方面，高速攝影儀在拍攝較大尺寸物體時像素數(shù)會呈指數(shù)降低，實驗照片對火焰的細節(jié)捕捉不夠；另一方面，數(shù)值模擬采用了Zimont 模型計算火焰?zhèn)鞑ヤh面，而該模型為了更好地捕捉到火焰鋒面位置，對火焰鋒面進行了增厚處理，導致在60 和62 ms 時對分支管道中火焰鋒面末端彎曲細節(jié)表達不精準。鑒于本文的研究目的，是為了分析火焰與流場耦合規(guī)律，因此可以忽略。

在火焰形態(tài)對比基礎上，對火焰位置和傳播速度進行了對比，如圖5～6 所示。由圖5 可以發(fā)現(xiàn)，數(shù)值模擬的火焰位置能夠較好地反應火焰?zhèn)鞑ノ恢秒S時間的變化，最大誤差為14.436%，發(fā)生在42 ms處。這是因為，此時火焰鋒面?zhèn)鞑ブ练种Ч艿栏浇?，爆炸流場復雜，火焰鋒面結構出現(xiàn)嚴重潰散，擾動隨機性較大，導致誤差增大。由圖6 可以發(fā)現(xiàn)，數(shù)值模擬的火焰速度與實驗的變化規(guī)律整體一致，在爆炸初期（28 ms 前），因參與反應的油氣較少，反應速率較小，火焰?zhèn)鞑ニ俣染徛?；隨著化學反應持續(xù)進行（實驗28～40 ms，數(shù)值模擬28～45 ms），火焰鋒面向末端傳播過程中沿軸向不斷拉伸，火焰表面積增大，參與反應的單位時間油氣增多，化學反應速率增大，進而火焰?zhèn)鞑ニ俣妊杆僭龃?；當火焰鋒面?zhèn)鞑サ椒种Ч艿栏浇鼤r，初期分支管道的泄壓作用較大，火焰?zhèn)鞑ニ俣仍鲩L趨勢出現(xiàn)停滯，其后隨著分支管道內(nèi)油氣被引燃，火焰鋒面表面積呈指數(shù)增加，火焰?zhèn)鞑ニ俣壤^續(xù)增大。數(shù)值模擬誤差最大為24%，發(fā)生在45 ms 處。

圖5 實驗和數(shù)值模擬的火焰位置Fig. 5 Experimental and simulated flame front locations

圖6 實驗和數(shù)值模擬的火焰速度Fig. 6 Experimental and simulated flame speeds

3.2 動態(tài)超壓變化規(guī)律

在壓力傳感器P1 處實驗和數(shù)值模擬的動態(tài)超壓曲線，如圖7 所示。表1 為實驗和數(shù)值模擬的3 個典型超壓峰值的對比。可見，兩者超壓曲線的規(guī)律基本一致，分別在65 和68 ms 達到爆炸超壓峰值，峰值誤差為8.6%。但在數(shù)值模擬中，未考慮薄膜破裂時的破膜壓力對管道內(nèi)壓力場的影響，因此在爆炸超壓上升階段，實驗存在3 個峰值，而數(shù)值模擬只存在2 個峰值。

表1 實驗和數(shù)值模擬的典型超壓峰值Table 1 Experimental and simulated typical overpressure peaks

在不考慮破膜壓力的情況下，爆炸超壓升壓過程可以劃分為4 個階段，以數(shù)值模擬得到的曲線為例：第Ⅰ階段為緩慢升壓階段（0～20 ms），該階段從模擬點火開始至火焰拉伸呈“指形”終止，在該階段爆炸集中于點火源附近，參與反應未燃氣量少，反應釋能小，超壓隨之上升緩慢；第Ⅱ階段為超壓上升階段（20～50 ms），隨著爆炸進行，參與反應的未燃氣逐漸增多，單位時間釋能迅速增大，超壓開始明顯上升；第Ⅲ階段為超壓平衡階段（50～56 ms），主管道泄流效應和分支管道二次泄壓作用弱化了內(nèi)場超壓上升，超壓上升速率減慢；第Ⅳ階段為快速升壓階段（56～65 ms），分支管道增大了燃燒面積，反應速率迅速上升，超壓再一次快速上升，直至可燃氣消耗殆盡轉(zhuǎn)入下降階段。

在小尺度狹長空間發(fā)生油氣泄爆時，通常管道內(nèi)壓力波動以靜壓為主，動壓對超壓的影響非常小，因此空間內(nèi)部各點的超壓基本一致。而對于含雙側分支的受限空間，當壓力波傳播經(jīng)過分支附近時，內(nèi)部空間橫截面的變化必然影響壓力波的傳播，因此在分支管道前后測點的超壓一定會有所差異。為此，對比了分支管道前后測點P1 和P2 處的超壓曲線，如圖8 所示。由圖8 可以看出，P1 和P2 處的超壓曲線演變規(guī)律基本一致，僅在局部有數(shù)值大小等方面的差異，具體有3 點：（1）在緩慢升壓階段和超壓上升階段，P2 處超壓上升滯后于P1；（2）P1 處最大超壓峰值大于P2，形成時刻略滯后于P2；（3）P2 處的超壓振蕩比P1 劇烈。P1 和P2 處的超壓特性差異說明，分支管道導致主管道橫截面突擴，在分支管道附近誘導稀疏波導致壓力波衰減，分支管道對主管道起到了泄壓作用，并且P2 位置更靠近泄爆口，受泄流效應的影響也更大。

圖7 P1 處實驗和數(shù)值模擬的爆炸超壓曲線Fig. 7 Experimental and simulated overpressure curves at point P1

3.3 流場-火焰-超壓耦合

圖8 P1 和P2 處數(shù)值模擬的爆炸超壓曲線Fig. 8 Simulated overpressure curves at points P1 and P2

由文獻[22]，在爆炸發(fā)展過程中，流場變化會影響火焰?zhèn)鞑バ袨?，火焰?zhèn)鞑バ袨榈母淖円鄷Ρ鲌霎a(chǎn)生影響，火焰?zhèn)鞑ズ屯牧髁鲌鲩g存在一定的正反饋耦合機制。此外，火焰結構在一定程度上表征反應速率的快慢，進而影響爆炸超壓的波動。

圖9 為火焰?zhèn)鞑バ螒B(tài)和流場速度矢量圖，白色箭頭代表速度矢量，顏色越白速度數(shù)值越大，速度矢量所在平面為管道中軸面（y=0.05 m）。由于重點關注分支管道附近的火焰和流場形態(tài)，因此火焰?zhèn)鞑ブ练种Ч艿狼埃?0、32 ms）的圖片包含點火端至分支管道；而為了清晰顯示火焰在分支管道中傳播形成“浪花狀”的過程，這個階段（48、60、64 ms）的圖片僅包含分支管道區(qū)域。

在20 ms 時，火焰以層流狀態(tài)傳播，火焰鋒面前端附近速度矢量沿半球面指向未燃區(qū)，主管道未燃區(qū)氣體在火焰膨脹作用下沿水平方向向管道末端流動。當流動經(jīng)過分支管道時，兩側壁面對未燃氣的約束消失，未燃氣對稱地向兩側分支管道膨脹，在主管道與分支管道交界處建立了層流流場，此時分支管道內(nèi)部流場速度趨近于零。隨著反應進行，在32 ms 時火焰以“指形”傳播，火焰面積增大，主管道軸向速度場迅速增強，由于受到分支管道尖角的擾動流場在分支管道口附近向湍流轉(zhuǎn)捩，且因分支管道右側壁面的約束，在主管道與分支管道交界處形成了對稱的渦旋結構，上側旋轉(zhuǎn)方向為逆時針，下側為順時針，此外分支管道內(nèi)部流場速度逐漸增大。當反應進行到48 ms、火焰?zhèn)鞑ブ练种Ч艿揽跁r，火焰鋒面在復雜湍流流場作用下顯著變形，而早期建立的對稱渦旋結構在火焰沖擊作用下向分支管道內(nèi)部發(fā)展且旋轉(zhuǎn)方向維持不變，但渦旋結構沿分支管道軸向拉伸成為較大的橢圓形，此時火焰面積較小且受到較強的泄壓作用，導致火焰?zhèn)鞑ニ俣壬仙厔轀p緩。此后，火焰?zhèn)鞑ミM入分支管道，在漩渦結構作用下火焰鋒面沿分支管道右側壁面拉伸，而后向左側壁面彎曲形成60 ms 時的“浪花狀”火焰。與此同時，火焰鋒面前側未燃氣不斷被引燃，火焰膨脹促進分支管道內(nèi)部速度矢量指向主管道，分支管道內(nèi)未燃氣在此流場作用下進入主管道，導致主管道與分支管道交界處形成新的燃燒面，如64 ms 時。

為了便于分析火焰形態(tài)和動態(tài)超壓的耦合關系，將P1 處數(shù)值模擬的動態(tài)超壓曲線和典型時刻火焰形態(tài)結構一起呈現(xiàn)，如圖10 所示。根據(jù)前述分析，爆炸超壓變化過程可劃分為4 個升壓階段和1 個振蕩下降階段。第Ⅰ階段，在點火初期，只有點火頭周圍少量油氣參與化學反應，而后火焰以半球形傳播，參與反應的可燃氣有限，爆炸緩慢發(fā)展，超壓上升趨勢不明顯。第Ⅱ階段，爆炸進入快速升壓階段，初期火焰以“指形”結構傳播（32 ms），火焰面積隨著鋒面拉伸不斷增大，參與反應的油氣增多，超壓逐漸上升，當火焰?zhèn)鞑ブ练种Ч艿栏浇鼤r，因流場復雜導致的火焰鋒面褶皺嚴重進一步導致反應更加劇烈，超壓上升趨勢更明顯。第Ⅲ階段，火焰?zhèn)鞑ミM入分支管道形成“偽受限空間”燃燒（54、56 ms），初始火焰?zhèn)鞑ニ俣容^慢，火焰面積變化不明顯，參與反應的油氣量沒有顯著增加，而火焰鋒面靠近末端泄爆口，主管道泄壓作用增強，導致這個階段升壓速率下降。第Ⅳ階段，當火焰在分支管道中渦旋結構的作用下充分發(fā)展時，火焰鋒面不規(guī)則潰散嚴重，同時參與反應的油氣量也大幅增多，爆炸劇烈進行，升壓速率較前3 階段都大。振蕩下降階段，火焰?zhèn)鞑サ?4 ms 時，由于主管道末端的泄壓效應，超壓迅速下降，而后在壓力波來回傳播和管道內(nèi)殘余油氣持續(xù)爆炸的兩方面作用下形成超壓振蕩。

圖10 數(shù)值模擬的超壓曲線和火焰?zhèn)鞑バ螒B(tài)Fig. 10 Simulated overpressure curve and flame propagation structure

4 結 論

利用大渦模擬WALE 模型，模擬了雙側分支結構受限空間油氣泄壓爆炸的超壓特性和火焰行為，與實驗數(shù)據(jù)對比了火焰?zhèn)鞑バ螒B(tài)、火焰?zhèn)鞑ニ俣群蛣討B(tài)超壓曲線，驗證了所選取模型對含多分支結構受限空間泄壓爆炸計算的有效性。然后，運用該模型對雙側分支結構受限空間油氣泄壓爆炸的空間流場結構、火焰行為和超壓特性進行了研究。得出以下結論：

（1）火焰?zhèn)鞑ミM入分支管道前，火焰行為會影響分支管道內(nèi)部流場，在主管道和分支管道交界處會產(chǎn)生方向相反的對稱渦旋結構，隨著火焰進一步傳播，渦旋結構會向分支結構內(nèi)部發(fā)展，建立旋轉(zhuǎn)流場；

（2）當火焰?zhèn)鞑ミM入分支管道后，分支管道內(nèi)部前期已建立流場決定了火焰的形態(tài)，火焰鋒面在渦旋結構作用下形成“浪花狀”火焰，此后火焰和流場相互影響，流場向湍流轉(zhuǎn)捩，火焰鋒面潰散；

（3）爆炸超壓升壓過程可劃分為4 個階段，受到火焰鋒面面積、主管道泄流效應和分支管道二次泄壓共同作用，爆炸流場、火焰行為和動態(tài)超壓具有顯著耦合性。