發動機爆震室氣液兩相爆轟性質

王晨辰,翁培奮,丁 玨,李家驊,陳永杰

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

脈沖爆震發動機是一種基于爆轟燃燒的新概念推進系統,其中爆震室內氣液兩相爆轟是發動機產生推力的重要過程,影響著發動機的動力和熱力性能.

燃料顆粒散布在空氣或氧氣中會形成可燃的混合介質,點燃該系統會發生兩相爆轟.1978年,Gabrijel等[1]研究了在楔形激波管中平均直徑為388μm的煤油液滴與空氣混合物形成的云團的爆轟.此后Dabora[2]建立了初步的理論模型來描述兩相爆轟,模型中采用了較多的經驗公式和假設.隨著計算流體力學的發展,數值模擬逐漸成為研究兩相系統爆轟的一種重要手段,其中Eulerian-Eulerian模擬方法應用較為廣泛.1999年中國工程物理研究院的洪滔等[3]利用一維雙流體模型對正庚烷、鋁粉顆粒在空氣、氧氣中的爆轟現象進行了研究,得出了燃料空氣比、指前因子及液滴半徑等因素對爆轟參數的影響;上海大學的丁玨等[4]分析了在雙流體模型下,兩相爆轟時壓力的最大上升速率以及氣相體積分數對爆轟波壓力的影響;南京理工大學的馬丹花等[5]在一維和二維雙流體模型的基礎上,發展了時空守恒元和解元(space-time conservation element and solution element,CE/SE)方法,較好地對爆轟波進行了研究.由于Eulerian-Eulerian方法對離散相無法較為準確地描述顆粒動力學過程,因此研究者利用Eulerian-Lagrangian模型開展對氣液、氣固兩相爆轟的研究.許厚謙等[6]模擬了激波點燃粉塵顆粒的一維非定常兩相化學反應爆轟模型;Cheatham等[7]對JP-10燃料爆轟進行了數值模擬研究,分析燃料液滴粒徑對爆轟壓力及爆速的影響;蔣弢等[8]對汽油液滴的爆轟情況進行了數值模擬,分析起爆距離與粒徑之間的關系.

兩相爆轟過程涉及液體燃料在前導激波及其氣流作用下變形、破碎、蒸發及化學反應等多個復雜過程.為了揭示爆轟波與顆粒動力學的作用機制,分析爆轟波參數與燃料液滴粒徑的關系,本工作開展了爆震室氣液兩相爆轟過程和顆粒動力學特性的數值研究,為推動爆震推進技術的應用奠定基礎.

1 氣液兩相燃燒爆轟反應模型

1.1 氣相控制方程

爆轟是個非定常的過程.

(1)氣相控制方程

連續方程:

動量方程:

能量方程:

組分輸運方程:

(2)顆粒相控制方程

式中,下標p為液滴;k為第k相組分;為液滴組分k的質量變化率;npk為顆粒數密度,即單位體積內k相顆粒的個數;Fpk為單個液滴所受的力;fpk為場力及其他諸力;qpk為單個液滴從氣相獲得的熱;β為分配因子,表示液滴釋放的熱在兩相間的分配;Q0為單位質量的液滴燃燒釋放的熱;Qs為氣相組分s的反應熱;ωs為反應速率;Df為質量擴散系數;Ys為該組分在氣相中的質量分數;αs為液滴蒸發產物中組分s的質量分數.

此外,氣相遵守理想氣體狀態方程:

式中,Ws為氣相組分s的物質的量.

1.2 化學反應模型

為了提高計算效率,采用一步總包化學反應模型:

假定燃燒速度ωr受指前因子A、氧氣濃度CO2、庚烷氣體濃度CC7H16及活化能Ea的影響,利用Westbrook等[9]的一步反應描述化學反應過程,有

式中,R為普適氣體常數,指前因子A及活化能Ea的值引自Westbrook等[10]的文獻.

燃燒模型采用部分攪拌反應器(partially stirred reactor,PaSR)模型[11-12]來考慮湍流效應的影響.PaSR模型引入了化學反應時間尺度τchem和混合時間尺度τmix,認為在τchem時間內組分基于化學反應速率ωr反應的產生/消耗量,與在(τchem+τmix)時間內以實際化學反應速率反應的產生/消耗量相同.因此,實際的化學反應速率為

式中,K為湍動能,ε為湍動能耗散率,Cmix為常數.各算例中,混合的時間尺度就是湍流時間尺度,因此Cmix取1.

1.3 數值格式

壓力與速度的耦合采用改進版的simple算法——pimple法,壓力場與速度場分別存放在2套不同的網格中,故壓力求解器選用的是代數多重網格求解器(geometric and algebraic multigrid,GAMG).在插值格式方面,時間的離散采用1階隱式Euler形式,梯度格式采用高斯線性插值,擴散相采用2階高斯守恒格式,表面插值采用中心差分格式.關于對流相的選擇,由于需要對爆震波進行捕捉,總變差減小(total variation diminishing,TVD)差分格式具有分辨率高、間斷面處無非物理振蕩、光滑區域精度高等優點,而具有TVD性質的MUSCL(monotonic upwind scheme for conservation law)格式[13]采用一種保持光滑解高階精度的通量限制器,抑制間斷附近的數值解波動,且在解的光滑區域內具有2階精度,因此是一類高分辨率捕捉激波格式.本工作的對流相采用MUSCL差分格式,由于運用拉格朗日法追蹤,故對于燃料液滴采用Euler形式對時間進行積分.

2 氣液兩相爆轟的數值模擬

圖1為計算區域,點火采取高溫高壓區域,該區域長1.00 m,高0.06 m,上、下和左側為固壁,右側為出口.流場內部壓力初始值為一個大氣壓,溫度為300 K.為使燃料氧氣按照等當量比進行化學反應,內部均勻預混分布著9.38×104個正庚烷液滴(化學式為C7H16,密度ρ=675 kg/m3).為便于計算,將計算域劃分為1 000×60個矩形網格.

圖1 計算區域Fig.1 Computed f i eld

2.1 空氣-燃料兩相系統的爆轟性質

為驗證計算模型,將本工作中模擬得到的爆轟波壓力與文獻[8]的實驗結果進行對比(見圖2).圖2中實線為實驗測得的汽油/空氣(汽油的主要成分為正庚烷)爆轟波壓力隨時間演化圖.文中計算的正庚烷與空氣混合物爆轟壓力曲線如圖2的虛線所示.可以看出,MUSCL格式能夠較為準確地捕捉激波,對激波的分辨率是很高的.

相比較而言,數值計算結果與實驗結果波形基本一致,均呈現出在爆轟波到達后流場壓力值先急劇上升,隨后再逐漸下降的現象.計算得到的壓力峰值為3.396 7 MPa,比實驗結果高出約0.5 MPa.偏差的原因主要是本工作中燃料為純庚烷,單位時間、單位體積內釋放的熱量更多.

為了驗證網格的無相關性,將網格數增加至1 500×90和2 000×120,得到的壓強峰值分別為3.408 7 MPa和3.412 5 MPa,相差僅為3.5%和4.7%,因而可以認為當網格數大于6×104后,計算結果受網格數量的影響不大.

(1)爆轟的C-J狀態.

根據ZND(von Zeldorich Neumann-Doering)模型,爆轟波由前導激波和反應區構成,反應區的末端平面稱為C-J面,反應區止于C-J面,該處氣相狀態被稱為C-J狀態.現將由較小粒徑的燃料液滴(25μm)組成的系統進行具體的分析.

圖2 數值模擬與實驗對比Fig.2 Comparison of numerical simulation and experimental data

圖3爆轟過程中的C-J狀態Fig.3 C-J state in detonation process

圖3 (a)為計算不同時刻爆轟場H點的p-v(壓力比體積)關系圖,其中各點表示在爆轟波經過前后,流場H點處氣體的狀態參數.如果將火焰前后氣體參數與釋放熱量Δq的關系式稱為Hugoniot方程,則由Δq值得到的一簇曲線稱為Hugoniot曲線簇.圖中,實線所連的點構成了Δq=0的Hugoniot曲線,亦稱為Hugoniot激波曲線,表示尚未反應的階段;標記點構成的曲線上Δq=1,被稱為Hugoniot產物曲線,表示反應完畢后的階段.以氣相參數初始狀態點O為原點,過該點作Hugoniot爆轟產物曲線的切線,則與Hugoniot曲線的切點(H點)被稱為爆轟的C-J狀態點,表示恰好進入穩定爆轟狀態.切線稱為Rayleigh線,即爆轟波的波速線.

將C-J狀態標記于所對應的壓力p、速度u及密度ρ曲線中(見圖3(b)～(d)).可以看到,C-J狀態H點位于曲線的鞍點處.通過將C-J狀態參數與文獻[14]中C-J理論值的對比(見表1)可以看出,小粒徑的燃料液滴(25μm),計算出的C-J狀態對應的爆轟參數與理論值較接近,驗證了本工作中建立的模型和研究方法是可行的.

表1 正庚烷/空氣爆轟參數Table 1 Detonation parameters of n-heptane/air mixture

(2)爆轟波的傳播過程.

圖4為半徑為25μm的正庚烷液滴與空氣混合發生兩相爆轟時,爆轟波在x方向傳播時的壓力、密度曲線.通過計算可知,穩定后的爆轟波傳播速度約為1 720 m/s.

圖4 軸線上爆轟波壓力、密度隨時間的演化Fig.4 Evolution of detonation pressure and density along axis

2.2 氧氣-燃料兩相系統的爆轟性質

2.2.1 爆轟波后顆粒蒸發域寬度的分析

兩相系統的爆轟波是由前導激波和激波后的化學反應區組成,其中爆轟波后燃料液滴的破碎、蒸發時間和寬度決定了反應區中能量的釋放速率.

圖5為激波到達距離爆震管封閉端0.6 m位置時該處燃料液滴粒徑的分布情況.橫軸坐標為爆震管軸線位置,縱軸坐標代表液滴的直徑.從圖5中可以看出,雖然燃料初始粒徑不同,但是在破碎蒸發等動力學過程中粒徑變化卻相似.將液滴從開始相變到結束的寬度Δx定義為液滴破碎蒸發域寬度.對于初始粒徑為25μm的正庚烷液滴,其爆轟波后顆粒破碎蒸發域寬度Δx=0.010 96 m,此后均相變為氣相的庚烷;而對于初始粒徑為50,75和100μm的液滴,破碎蒸發域寬度則分別為0.022 81,0.032 29和0.041 21 m.

圖5爆轟波后粒徑分布Fig.5 Droplet size distribution behind detonation wave front

圖6 為計算所得的液滴初始半徑與破碎蒸發域寬度.圖中,黑點為計算所得的液滴初始半徑與破碎蒸發域寬度,紅線為擬合曲線.通過擬合分析可知,擬合函數滿足直線變化關系,說明了在本工作討論的范圍內,隨著粒徑的增大,液滴破碎蒸發域寬度增大;液滴破碎蒸發域寬度為粒徑的函數,二者存在線性關系.

圖6 顆粒初始半徑與破碎蒸發域寬度Fig.6 Initial droplet size and width of break and evaporation region

2.2.2 燃料液滴的初始粒徑對爆轟壓力的影響

圖7為當兩相反應的化學當量比為1時,不同粒徑的庚烷液滴在氧氣中爆轟的壓力分布情況.圖中虛線為文獻[3]數值計算值.從圖中可以看出,隨著粒徑的增大波形的寬度不斷增大,而爆轟的壓力峰值逐漸降低.這是因為較大粒徑的庚烷液滴的破碎和蒸發的特征時間較長,使得反應速率降低,導致爆轟波壓力峰值也隨之減小.

圖7 氧氣中爆轟時軸線壓力隨時間變化Fig.7 Time variations of the pressure along the axis while detonation happened in oxygen

當氣液相化學當量比為1時,不同粒徑庚烷液滴在氧氣中爆轟的相關參數如表2所示.

表2 氧氣中兩相爆轟參數與粒徑關系Table 2 Relations between detonation parameters and droplet sizes in oxygen

2.2.3 流場組分的分布及反應區寬度

爆轟波經過x=0.6 m時,液滴在不同初始粒徑下流場化學反應組分體積分數的分布如圖8所示.從圖中爆轟波后化學反應區寬度與液滴半徑的關系可看出,液滴的粒徑越大,反應區寬度也越大.

圖8 爆轟反應氣體組分體積分數的分布Fig.8 Distributions of detonation species volume fraction

3 結論

氣液兩相爆轟是含化學反應的多相流體力學過程.本工作重點研究兩相爆轟波性質,以及液滴初始粒徑對爆轟波結構的影響.

(1)數值分析兩相爆轟的熱力學過程.對于空氣-庚烷燃料兩相系統的爆轟,計算爆轟場H點的p-v圖,得到Hugoniot曲線和Rayleigh曲線.對于初始粒徑為25μm的燃料液滴,計算所得到的在C-J狀態下爆轟壓力、密度、氣流速度、爆轟波傳播速度分別為1.802 2 MPa,2.161 3 kg/m3,721 kg/m3和1 720 m/s,與理論值接近,驗證了本工作中建立的模型和研究方法的正確性.

(2)進行兩相爆轟過程的數值研究.燃料液滴在前導激波及在氣流作用下會發生破碎、蒸發及化學反應的動力學過程.基于此,本工作中數值分析了液滴初始粒徑對其破碎蒸發域寬度、爆轟波反應區寬度的影響.研究結果發現,在均勻預混的爆轟下,燃料液滴初始粒徑越大,液滴蒸發、反應區寬度越大,爆轟波壓力峰值越小,液滴破碎蒸發域的寬度是初始粒徑的函數.

[1]GABRIjEL Z,NICHOLLS J A.Cylindrical heterogeneous detonation waves[J].Acta Astronautica,1978,5:1051-1061.

[2]DABORA E K.A model for spray detonation[J].Acta Astronautica,1978,6:269-280.

[3]洪滔,秦承森.氣體-燃料液滴兩相系統爆轟的數值模擬[J].爆炸與沖擊,1999,19(4):336-342.

[4]丁玨,應夢侃.基于兩相流模型數值研究可燃顆粒燃燒爆轟的特性[J].中國安全生產科學技術,2014,10(3):17-20.

[5]馬丹花,翁春生.二相湍流爆轟流場的數值仿真[J].計算機仿真,2012,29(6):85-88.

[6]許厚謙,于強.激波點燃粉塵的數值模擬研究[J].爆炸與沖擊,1994,14(4):290-297.

[7]CHEATHAM S,KAILASANATH K.Multi-phase detonations in pulse detonation engines[C]//AIAA Aerospace Sciences Meeting and Exhibit.2006.

[8]蔣弢,翁春生.氣液兩相脈沖爆轟發動機的建模和仿真[J].計算機仿真,2012,29(8):77-80.

[9]WESTBROOK C K,DRYER F L.Simplif i ed reaction mechanisms for the oxidation of hydrocarbon[J].Combustion Science and Technology,1981,27:31-43.

[10]WESTBROOK C K,DRYER F L.Simplif i ed reaction mechanisms for the oxidation of hydrocarbon fuels in f l ames[J].Combustion Science and Technology,1981,27:31-43.

[11]GOLOVITCHEV V I,CHOMIAK J.Evaluation of ignition improvers for methane auto-ignition[J].Combustion Science and Technology,1998,135:31-47.

[12]GOLOVITCHEV V I,NORDIN N,CHOMIAK J.Numerical modeling of hydrogen def l agration in multi-room compartments[J].Proceedings of the Computational Technologies for Fluid/Thermal/Chemical Systems with Industrial Applications,1998,377(2):169-176.

[13]VAN LEER B.Towards the ultimate conservative diあerence scheme.V.a second-order sequel to Godunov’s method[J].Journal of Computational Physics,1979,32(1):101-136.

[14]洪滔.兩相爆轟的理論和數值研究[D].綿陽:中國工程物理研究院,2004:38-39.本文彩色版可登陸本刊網站查詢:http://www.journal.shu.edu.cn