爆炸驅動燃料云團界面不穩定性發展及顆粒分散過程

申 洋， 丁 玨， 翁培奮

(上海大學上海市應用數學和力學研究所， 上海200072)

可燃氣體、液固燃料空氣混合物的爆炸性質對工業生產安全造成嚴重危害， 如金屬顆粒在工廠發生的粉塵爆炸， 但同時因其爆炸毀傷效應而在國防和軍事有巨大用途. 利用中心裝藥爆炸驅動燃料運動形成云團， 這樣的燃料空氣混合物， 稱作燃料空氣炸藥(fuel air explosive，FAE)[1]. FAE 是以揮發性液體碳氫化合物或固體粉質可燃物為燃料， 以空氣的氧氣為氧化劑組成的非均相爆炸性混合物， 具有能量高、分布爆炸性、原料易得等特點.

燃料的爆炸驅動是提高燃料空氣炸藥裝置威力的關鍵技術， 國內外學者在爆炸作用下FAE 裝置殼體破裂[2-4]、燃料界面發展[5-6]及射流形成[7-9]、液體燃料運動及破碎[10-13]、FAE云團形成的規律性[14]等方面開展了很多有意義的研究.

李磊等[15]設計了一種有上下平面約束的液體爆炸分散裝置， 應用平面激光誘導熒光技術得到了徑向膨脹液體環內誘導的熒光照片. Kandan等[16]通過激波管實驗和數值計算， 研究了沖擊波加載下燃料不穩定性的發生和發展. Shi等[17]在垂直矩形激波管中， 利用半導體激光器、信號放大電路、數字示波器、高速攝影， 以及自行設計的測量系統測量了激波驅動的界面速度;并對Richtmyer-Meshkov 不穩定性(Richtmyer-Meshkov instability， RMI)的機理進行了實驗研究. 結果表明: 界面速度隨馬赫數的增加而增大， 界面運動與混合區寬度的距離與時間呈線性關系.

Frost等[18]基于實驗和數值方法研究了爆炸驅動下顆粒到達不同位置時的速率， 分析顆粒和周圍流場動量和能量之間的交換. 廖斌等[19]采用實驗與數值計算方法對沖擊作用下液滴在環境液體中的演變過程進行了研究， 重點討論液膜破碎之前液滴的4 種演變模式和流場.

發生兩相爆轟的條件如下: 空間懸浮一定濃度的燃料， 形成燃料云團; 而爆炸驅動液體燃料擴散是液固燃料空氣混合物形成的重要過程. 以往對于爆炸驅動燃料分散的研究主要集中于云團形狀和狀態， 而對于爆炸驅動最初階段的力學性質， 由于受到實驗或試驗條件的限制，研究成果較少， 特別是涉及波、渦與燃料顆粒的相互作用機制、氣液界面的不穩定演化規律尚未清晰. 因此， 本工作開展液體燃料在爆炸沖擊波作用下分散機制的研究， 為工業生產的安全防治以及國防領域燃料空氣炸藥威力的提高提供理論指導.

1 爆炸驅動燃料云團運動的物理數學模型和數值計算方法

在工業生產、國防軍事領域， 沖擊波驅動燃料使得顆粒處于高溫、高壓和高速氣流中， 速度和溫度差導致兩相間動量和能量的傳遞; 同時， 顆粒發生蒸發效應， 使兩相間出現質量傳遞，逐漸形成由可燃氣體、燃料與空氣混合的云團. 本工作將重點研究爆炸驅動下流場和燃料云團的運動性質.

爆炸驅動裝置通常為柱對稱結構， 通過中心軸裝填的高能炸藥爆炸釋放的能量來驅動燃料(見圖1). 圖2 顯示燃料質量為12.5 kg 爆炸驅動裝置在58 ms 時刻云團的狀態.

圖1 爆炸驅動裝置Fig.1 Explosive driven device

圖2 爆炸驅動燃料云團的形成Fig.2 Fuel cloud formation by explosion-driven

1.1 燃料云團湍流運動的物理數學模型

針對爆炸驅動燃料云團的形成和運動過程， 本工作建立相間耦合的兩相流模型來描述， 即把氣體當作連續介質(流體相)， 把燃料作為離散體系(顆粒相)， 用歐拉-拉格朗日方法來描述，并考慮液體燃料的蒸發效應.

1.1.1 流體相控制方程

本工作討論的是非定常的流場， 二維多相流的控制方程為

式中：Φ 為通用變量， 分別表示 u， v， k， ε， H 和 Y 等求解變量; ρ 為密度; r 為柱坐標中的徑向; Uj為j 方向的速度; Γ 為擴散系數; 在直角坐標和柱坐標下， 上標m = 0 時用于平面流動，m = 1 時用于軸對稱流動; SΦ和SPΦ分別表示氣相本身的項和液體燃料顆粒蒸發引起的額外源項.

為了分析氣液界面不穩定的發展， 刻畫沖擊波誘導旋渦產生的過程， 本工作采用成熟的、適用于高雷諾數流場的標準k-ε 湍流模型， 數值研究氣體的湍流性質.

氣相遵守理想氣體狀態方程，

1.1.2 顆粒相控制方程

假設燃料顆粒的初始粒徑相同， 且顆粒間的相互作用可以忽略. 控制方程為

1.1.3 液體燃料蒸發的動力學方程

本工作討論的燃料為庚烷. 顆粒在波后加速， 由于對流傳熱升溫發生蒸發， 形成燃料蒸氣.根據顆粒蒸發模型， 其質量變化關系式為

式中: Sc 是施密特數; Sh 是Sherwood 數， 表示對流作用對傳質過程的影響，

B 為 Spalding 數， 也稱質傳遞數，

其中Ys表示顆粒表面的蒸汽質量分數， Y 為燃料蒸汽在氣相中的質量分數.

顆粒表面溫度變化率為

式中: Q 代表氣相傳遞給液相燃料的熱量; Nu 為努塞爾數; Ts為顆粒表面溫度; T∞為環境溫度; Lv為顆粒氣化潛熱; cp，l為液體等壓比熱容; D 為顆粒直徑; λ 為導熱系數.

1.2 燃料云團湍流運動的數值計算方法

本工作采用歐拉-拉格朗日方法， 利用壓力隱式算子分裂(pressure implicit split operator，PISO)算法和壓力耦合方程組半隱式方法(semi-implicit method for pressure linked equations，SIMPLE)相結合的PIMPLE 算法， 其中對流項用二階精度的守恒型單調迎風格式(monotonic upwind scheme for conservation laws， MUSCL)， 數值研究爆炸沖擊波的驅動以及燃料云團的形成過程.

為了驗證數值格式的精度， 針對Sod 問題開展數值驗證[20]. Riemann 類型的初始條件如下: (ρ，v，p)L=(1，0，1); (ρ，v，p)R=(0.125，0，0.100).

圖3 和4 給出了t=2 時刻壓力和密度數值解與精確解的對比. 由圖可以看出， 數值解與精確解吻合較好， 數值格式對波和間斷具有較高的分辨率.

為了刻畫云團的形成性質， 開展對爆炸驅動12.5 kg 燃料云團運動過程的數值模擬. 研究結果顯示， 液體燃料分散是云團在空間逐漸形成和展開的動力學過程， 即中心裝藥發生爆轟， 而后爆轟產物氣體作用于周圍介質形成沖擊波， 波和爆轟產物驅動燃料運動. 在驅動初期， 爆炸作用力占主導; 隨后的階段， 氣動阻力占主導， 逐漸形成燃料和空氣混合的可爆性云團. 圖5 和6 給出了云團尺寸的計算結果與實驗數據的比較. 可以看到: 在60 ms 時間內， 計算所得的云團演化尺寸與實驗數據非常接近， 即在爆炸作用力下， 云團尺寸發展較快; 當云團內的壓力接近外界環境壓力時， 云團膨脹運動變慢， 云團尺寸的變化趨緩. 計算所得的云團尺寸隨時間演化滿足指數函數規律: 半徑R = 5.026 46-5.137 63e-0.079t; 高度H = 1.774 59-1.853 92e-0.084t(t : 0 ～60 ms). 按照云霧膨脹的相似率[21]， 函數中隨時間t的增長系數(云團半徑增長系數為0.079， 云團高度增長系數為0.084)反映了中心裝藥釋放的總能量， 與爆炸驅動力密切相關.

圖3 Sod 問題的壓力曲線Fig.3 Pressure curve for Sod problem

圖4 Sod 問題的密度曲線Fig.4 Density curve for Sod problem

圖5 燃料云團半徑隨時間的變化Fig.5 Fuel cloud radius over time

圖6 燃料云團高度隨時間的變化Fig.6 Fuel cloud height over time

2 爆炸驅動下沖擊波、湍流渦與燃料之間的相互作用

在爆炸驅動最初期， 沖擊波作用于燃料， 氣液界面發生演化， 由初始較為細微的擾動逐漸發展成湍流. 同時， 波與液體燃料分界面作用發生折射， 氣液界面獲得加速， 促進不穩定性發展， 這些對燃料爆炸拋撒、燃料混合起著重要作用.

下面將著重研究和分析最初期兩相耦合的動力學性質. 圖7 給出了計算區域示意圖， 其中區域尺寸為 170 mm×170 mm(長×高). 為了簡化， 沖擊波的效應采取10 mm×170 mm(長×高)高溫高壓(T =1 000 K，P =1 519.875 kPa)區域來代替.Γ1為對稱平面，Γ2為區域的對稱軸，Γ3和Γ4為無反射邊界條件.燃料為液態庚烷，質量為4.242 mg. 在數值模擬中，為了比較和分析數值研究的網格，跟蹤流場中的點(0.06 m，0.04 m)處的參數. 當計算的網格尺寸分別為0.5 和0.8 mm 時，流場跟蹤點處壓力峰值分別為578.404 和568.707 kPa，與網格尺度1 mm(壓力峰值567.106 kPa)分別相差1.70%和0.77%. 綜合考慮計算量和誤差， 采取尺寸為1 mm 的網格系統， 數值模擬中網格數目為28 900. 圖8 給出了燃料云團的示意圖， 其中藍色箭頭代表誘導渦旋轉方向為順時針， 而紅色箭頭代表逆時針.

圖7 計算區域Fig.7 Computation field

圖8 燃料云團示意圖Fig.8 Schematic diagram of fuel cloud

2.1 兩相介質中沖擊波的性質

2.1.1 沖擊波的動力學性質

波的傳播是隨時間演化的非定常運動過程， 波后介質會發生動力學和熱力學狀態的變化.當高溫高壓氣體產生的沖擊波作用于氣相和燃料介質中， 在介質內存在能量傳遞過程. 壓縮作用也由近及遠地傳遞到介質其他部分. 在傳遞過程中， 擾動部分與未擾動部分的分界面是波的波陣面. 由于空氣和燃料介質中波傳播時的阻抗不同， 波傳播速度不同， 因此波在兩相中陣面發生彎曲(見圖9).

圖9 沖擊波在兩相介質中的傳播Fig.9 Transmission of shock waves in two-phase media

圖10 給出了在燃料介質和氣相介質中沖擊波壓力的分布情況. 可以看到: 燃料介質中壓力分布存在2 個峰值(peak1 和peak2)， 而氣相介質中只存在1 個壓力峰值. 分析原因主要是波在燃料中傳播時發生反射， 反射波朝爆源運動時， 由于波的匯聚導致壓力升高， 在爆源附近又形成一個新的波. 隨后， 這個新生成的波開始在燃料介質中向前傳播. 此外， 波掃過燃料， 燃料介質的動量和能量發生變化. 在40 μs 時刻前， 燃料云團中波的壓力峰值較空氣中更大. 由于稀疏波的影響， 壓力峰值隨時間逐漸降低. 在105 μs 時刻， 燃料云團中壓力峰值的超壓值降至0.21 MPa 左右.

圖10 兩相介質中沖擊波壓力峰值在空間的分布Fig.10 Spatial distribution of shock wave pressure peaks in two-phase media

氣相和燃料介質中波的超壓峰值隨時間的演化如圖11 所示. 選取波與燃料云團開始的21 μs 為初始作用時刻， 波即將穿過云團的105 μs 為結束時刻. 計算結果顯示: 波在氣相場中傳播時， 由于受到燃料介質的影響， 在 21～105 μs 時間段內， 壓力峰值相對穩定， 衰減較為緩慢; 然而， 在 21～105 μs 時間段內， 燃料云團介質內超壓發生快速衰減. 第一道波在 66 μs 時刻后穿過云團， 第二道波于105 μs 時刻后離開云團.

圖11 氣相和燃料介質波超壓峰值的演化Fig.11 Evolution of overpressure peaks in gas and fuel medium

沖擊波在氣相、燃料介質中超壓峰值的演化均滿足冪數律的關系，

式中: a， b， c 為系數; Z(m·kg-1/3)為比例距離， 由沖擊波傳播的距離和中心裝藥量的立方根來確定[22].

由于氣相和燃料介質中波的傳播特性不同， 波的衰減性質略有區別.

圖12 給出了在31 μs 時刻， 波陣面位置及中軸線上流場溫度、燃料顆粒溫度的分布. 可以看到， 31 μs 時刻沖擊波在云團內部傳播(沖擊波在平面中的位置用藍色線標記)， 流場溫度、顆粒溫度增加， 因此顆粒的蒸發動力學效應增強.

圖12 t=31 μs 時刻在y =0 m 軸線上氣相場溫度以及燃料蒸氣濃度的分布Fig.12 Distribution of gas field temperature and fuel vapor concentration along y =0 m axis at 31 μs

2.1.2 液體燃料的運動

燃料的運動取決于流場的氣體動力學與燃料云團的相互作用. 在321 μs 時間內， 燃料云團在x， y 方向上獲得的平均速度和加速度演化如圖13 所示. 可以看到， 在爆炸驅動最初期，沖擊波作用使得燃料獲得較大的徑向、垂向加速度; 隨后， 燃料在反射波和氣液界面湍流渦的影響下， 徑向加速度急劇下降， 運動速度減慢; 在200 μs 時刻后， 由于云團內側附近流場壓力的升高， 推動燃料介質繼續沿徑向運動， 加速度又逐漸上升. 圖 13 中， 在 112 μs 時刻， x 方向上的速度幾乎達到最大值22.8 m/s; 而y 方向上的速度在295 μs 時刻達到最大值8.97 m/s.

圖13 燃料云團的平均速度和加速度Fig.13 Average velocity and acceleration of fuel cloud

分析計算所得的平均速度和加速度曲線可知， 在氣動力作用下， 燃料介質加速度的變化存在急劇上升和下降2 個階段.

式中: ax1為常數; axmax為加速度的最大值; b 為衰減系數. 在氣動力作用下， y 方向燃料的加速度也滿足上述函數關系.

分析和整理本工作的數值計算結果， 給出加速度的指數函數關系式如下.

可見， 燃料獲得 x， y 方向加速度的增長系數分別為 49.92 和 8.26; 隨后， 燃料在 x， y 方向加速度的衰減系數降至1.07 和1.97. 增長和衰減系數反映出氣動力和大、小尺度渦對燃料運動的作用效應， 這也使得燃料介質呈現復雜的運動性質.

2.2 沖擊波誘導渦與燃料界面的相互作用

在爆炸驅動的最初期階段， 向前傳播的沖擊波作用于燃料云團， 界面附近產生Richtmyer-Meshkov 不穩定性(RMI)， 流場誘發出尺寸不一的旋渦， 加速了燃料的混合和擴散.

2.2.1 氣液界面旋渦與液體燃料的相互作用

圖14 給出了顆粒初始粒徑為60 μm 的燃料云團氣液界面流場渦量隨時間的演化. 由圖可見， 由于剪切效應占主導， 燃料界面的內側外緣附近流場產生了大尺寸的旋渦結構. 在40 μs 時刻， 最大渦心強度達到 1.42×105s-1， 極大地促進了顆粒的分散.

此外， 本工作研究結果表明， 當沖擊波驅動液體燃料時， 由于密度和壓力梯度導致斜壓效應， 氣液界面誘導出渦量， 界面產生Richtmyer-Meshkov 不穩定性(RMI)， 并隨時間擾動逐漸發展. 圖15 給出了270 μs 時刻， 粒徑為60 μm 的云團界面附近流場出現的尺寸不一的渦串結構. 在隨后湍流渦的作用下， 擾動進一步發展， 燃料云團與周圍流場發生混合.

2.2.2 誘導渦的發展和演化

下面將針對顆粒初始粒徑為60 μm 的燃料云團， 波誘導出渦對隨時間的演化過程進行分析.

(1) 擾動沿y 方向的發展和演化.

圖16 給出了4 個時刻渦量在界面流場中的分布. 可以看到， 隨時間的發展， 渦的強度逐漸降低， 即 10 和 17 μs 的渦量大約為 2 000 和 1 600 s-1， 大于 82 和 91 μs 時的渦量 1 200 和800 s-1.

圖 14 60 μm 云團流場渦量隨時間的演化Fig.14 Vorticity evolution of flow field over time for 60 μm radius cloud

圖15 270 μs 時刻流場渦量和顆粒的分散狀態Fig.15 Vorticity of flow field and dispersion of particles at 270 μs

圖16 不同時刻云團內側界面流場渦量Fig.16 Flow vorticity in the inner interface of cloud for different time

圖 17 給出了 173 μs 時刻， 分別在 x=0.021 m(side1)和 x=0.024 m(side2)分界線上渦對的分布(分界線位置見圖8)， 可以清晰地觀察到界面附近渦對的結構. 此外， 在渦對的卷吸作用下， 液體燃料的蒸發動力學效應增強. 在y =0 m 的直線上， x=0.020 8 ～0.021 6 m 之間的區域燃料蒸汽體積的體積分數最大， 約為0.526.

圖17 173 μs 時刻云團內側界面流場渦量Fig.17 Flow vorticity in the inner interface of cloud at 173 μs

(2)擾動沿主流方向的發展和演化.

擾動在燃料介質中傳播， 渦的強度逐漸降低. 在云團的外側界面， 由于存在壓力梯度和密度梯度， 界面同樣出現類似的不穩定現象， 渦對結構不斷演化. 圖18 給出了y =0 m 軸線上流場渦沿x 軸的演化. 順時針和逆時針旋轉的渦分布在氣體和液體燃料內側、外側界面.

2.2.3 誘導渦對液體燃料的影響

圖19 給出了顆粒初始粒徑分別為60 和80 μm 的云團內側界面流場渦量的演化情況. 可以看到， 對于60 μm 顆粒組成的云團界面附近流場誘導的旋渦， 渦運動的方向為順時針， 渦量為負值; 而對于80 μm 顆粒組成的云團界面附近流場誘導的旋渦， 渦運動的方向為逆時針， 流場性質略有不同. 顆粒分布如圖20 所示.

旋渦運動沿著液體燃料界面伸展， 隨著不穩定性的發展， 氣液混合增長速率隨小尺度渦結構的發展得到增強. 因此， 擾動產生的大尺度旋渦和小尺度旋渦使界面不斷變形， 加劇了顆粒的運動和分散.

圖18 y =0 m 軸線上渦量的分布Fig.18 Vorticity distribution along the axis of y =0 m

圖19 y =0 m 軸線上流場渦量隨時間的演化Fig.19 Flow vorticity over time along the axis of y =0 m

圖 20 400 μs 時刻顆粒初始粒徑分別為60， 80 μm 的燃料云團狀態Fig.20 Fuel cloud with initial particle radius 60， 80 μm at 400 μs

3 結 論

爆炸驅動燃料分散過程存在沖擊波、湍流、顆粒動力學等復雜物理場的相互作用， 本工作開展了爆炸作用沖擊波與燃料界面相互作用過程的數值模擬， 重點分析最初期波、誘導湍流渦對燃料顆粒分散的作用機制.

(1) 爆炸驅動液體燃料分散是云團在空間展開的、逐漸形成燃料和空氣混合的可爆性云團過程. 在爆炸作用力下， 云團尺寸發展較快; 當云團內的壓力接近外界環境壓力時， 云團運動變慢， 云團尺寸變化趨緩. 計算所得的云團尺寸隨時間按照指數律發展， 滿足云霧膨脹的相似率: 半徑的時間演化， R = 5.026 46-5.137 63e-0.079t; 高度的時間演化，H =1.774 59-1.853 92e-0.084t(t:0 ～60 ms). 云團半徑增長系數(0.079)和云團高度增長系數(0.084)反映了中心裝藥釋放的總能量與爆炸驅動力密切相關.

(2) 在爆炸驅動最初期， 波在兩相介質中傳播時發生衰減， 氣相和燃料介質中超壓峰值的變化滿足冪數律的關系. 此外， 由于受到空氣和燃料介質的阻抗作用， 云團中波傳播速度較慢，因此相對于空氣介質， 燃料介質中波陣面發生彎曲. 同時， 波后云團狀態發生變化， 燃料的運動速度提高. 如在112 μs 時間內， 初始粒徑為60 μm 的燃料云團在x 方向上的運動速度增加至 22.8 m/s， 最大加速度為 5.0×105m/s2.

(3) 波繞過燃料云團的外緣產生大尺度的旋渦. 對于顆粒初始粒徑為60 μm 的燃料云團，大尺度渦的強度最大達到1.42×105s-1， 驅動著顆粒向外分散. 同時， 氣液界面附近流場密度和壓力梯度誘導出渦量， 界面產生不穩定性現象， 并隨時間逐漸發展成沿界面排列的多個渦對結構， 流場湍流得到發展.

因此， 在爆炸驅動最初階段， 流場湍流度的增加疊加繞射波的雙重作用， 加劇了氣液的混合， 也為燃料顆粒的湍流分散提供了較大的初始動能.