氣體與液體兩相連續(xù)旋轉(zhuǎn)爆轟波二維數(shù)值模擬研究

王放， 翁春生， 武郁文， 白橋棟， 鄭權(quán)

(南京理工大學(xué) 瞬態(tài)物理國(guó)家重點(diǎn)實(shí)驗(yàn)室， 江蘇 南京 210094)

0 引言

爆轟是一個(gè)由強(qiáng)反應(yīng)激波誘導(dǎo)的燃燒過程，它與傳統(tǒng)的等壓燃燒相比具有更高的熱效率[1]。基于爆轟燃燒的發(fā)動(dòng)機(jī)還具有結(jié)構(gòu)緊湊、推重比大的優(yōu)點(diǎn)。連續(xù)旋轉(zhuǎn)爆轟發(fā)動(dòng)機(jī)(RDE)是一種利用一道或多道爆轟波在環(huán)形燃燒室內(nèi)連續(xù)傳播從而產(chǎn)生推力的發(fā)動(dòng)機(jī)，近年來已成為國(guó)內(nèi)外學(xué)者研究的熱點(diǎn)。燃料和氧化劑通常采用分開噴注方式從頭部的環(huán)縫或者孔進(jìn)入燃燒室，連續(xù)旋轉(zhuǎn)爆轟波(RDW)在燃燒室頭部周向傳播，傳播方向與燃料噴注方向垂直，高壓爆轟產(chǎn)物側(cè)向膨脹，通過噴管排放到大氣中產(chǎn)生推力。

實(shí)驗(yàn)研究方面，俄羅斯Bykovskii等[2-3]使用不同燃料(包括氫氣、甲烷、丙烷、煤油、汽油等)在多種尺寸RDE上進(jìn)行了大量實(shí)驗(yàn)，獲得兩相RDW結(jié)構(gòu)如圖1所示。魏萬里等[4]通過實(shí)驗(yàn)研究了氧化劑噴注面積對(duì)氣態(tài)RDE的影響。文獻(xiàn)[5-6]成功起爆了汽油與富氧空氣RDE，獲得的RDW波速為1 022.2～1 171.8 m/s.

圖1 實(shí)驗(yàn)所得不同燃料下RDW結(jié)構(gòu)[2]Fig.1 Experimental results of RDW structure with different mixtures[2]

數(shù)值模擬研究方面，國(guó)內(nèi)外學(xué)者對(duì)氣態(tài)燃料RDE進(jìn)行了廣泛深入的研究，取得了豐富的成果。Liu等[7]使用5階加權(quán)本質(zhì)非振蕩(WENO)方法和9組分化學(xué)反應(yīng)模型對(duì)H2/Air RDE進(jìn)行了數(shù)值模擬，發(fā)現(xiàn)了低頻不穩(wěn)定性和高頻不穩(wěn)定性同時(shí)存在于單波和雙波模態(tài)RDW的現(xiàn)象。Fujii等[8]通過non-MUSCL-type 2階迎風(fēng)方法研究了C2H4和O2噴射條件對(duì)RDW波速的影響，發(fā)現(xiàn)在非勻混進(jìn)氣條件下已燃?xì)怏w對(duì)RDW波速有著重要的影響。Sun等[9]使用ANSYS/Fluent軟件研究了非預(yù)混進(jìn)氣條件的H2/Air RDE，發(fā)現(xiàn)空氣噴孔寬度和總質(zhì)量流量會(huì)影響RDW的傳播模態(tài)，RDW的自持傳播需要較高的總質(zhì)量流量。以上研究揭示了氣態(tài)燃料RDW流場(chǎng)結(jié)構(gòu)演變規(guī)律，分析了氣態(tài)燃料旋轉(zhuǎn)爆轟的特性。李寶星等[10-11]用結(jié)構(gòu)化網(wǎng)格守恒元與求解元(CE/SE)方法對(duì)汽油與富氧空氣兩相爆轟進(jìn)行了數(shù)值模擬，獲得了單波模態(tài)RDW的爆轟特性。目前國(guó)內(nèi)外學(xué)者關(guān)于液態(tài)燃料RDW的數(shù)值研究較少，缺少對(duì)于氣體與液體(簡(jiǎn)稱氣液)兩相連續(xù)旋轉(zhuǎn)爆轟流場(chǎng)結(jié)構(gòu)和爆轟特性的認(rèn)識(shí)。

非結(jié)構(gòu)化CE/SE方法最早由Chang[12]提出，Shen等[13-14]拓展了非結(jié)構(gòu)化CE/SE方法，獲得了很好的計(jì)算結(jié)果。Nisar等[15]提出了新的守恒元和求解元分裂方法。這些研究帶來了非結(jié)構(gòu)化網(wǎng)格CE/SE方法研究的熱潮。非結(jié)構(gòu)化CE/SE方法除了具有計(jì)算格式簡(jiǎn)單、精度高、捕捉間斷能力強(qiáng)等優(yōu)點(diǎn)外，還可以靈活劃分計(jì)算域內(nèi)的網(wǎng)格，在三維復(fù)雜流場(chǎng)計(jì)算上更具潛質(zhì)。目前還沒有發(fā)現(xiàn)該方法用于計(jì)算氣液兩相連續(xù)旋轉(zhuǎn)爆轟的文獻(xiàn)報(bào)道。

為了解決氣液兩相旋轉(zhuǎn)爆轟數(shù)值計(jì)算問題，探究氣液兩相旋轉(zhuǎn)爆轟的流場(chǎng)結(jié)構(gòu)和爆轟波參數(shù)變化規(guī)律，本文利用非結(jié)構(gòu)化三角形網(wǎng)格CE/SE方法計(jì)算汽油與空氣兩相旋轉(zhuǎn)爆轟過程，開展并行化程序設(shè)計(jì)，計(jì)算在不同當(dāng)量比和不同進(jìn)氣總壓下的RDW傳播特性，獲得兩相RDW不同模態(tài)的流場(chǎng)結(jié)構(gòu)，分析不同模態(tài)RDW的工作特性。

1 物理模型

1.1 計(jì)算模型

由于RDE燃燒室的厚度通常遠(yuǎn)小于其半徑，為方便計(jì)算，徑向氣流變量可以忽略，RDE的數(shù)值計(jì)算可以簡(jiǎn)化為二維問題[16-18]。 沿著環(huán)形燃燒室的1條母線剪開，可以得到二維矩形計(jì)算區(qū)域。

如圖2所示，計(jì)算區(qū)域?yàn)殚L(zhǎng)0.3 m、寬0.1 m的矩形計(jì)算域。左右兩邊是周期邊界，參數(shù)可以在二者之間互換。上邊界是無反射出口邊界，當(dāng)出口速度為超聲速時(shí)，出口邊界參數(shù)由2階外推得到；當(dāng)出口速度為亞聲速時(shí)，出口邊界輔助網(wǎng)格的壓力和溫度分別為0.1 MPa和298 K，其余物理量包括密度、周向速度和軸向速度均與邊界網(wǎng)格相同。下邊界為入口邊界， Schwer等[19]以及Frolov等[20]已經(jīng)表明燃燒產(chǎn)物很少穿透至上游，因此此入口條件也忽略了回火的影響。入口邊界各參數(shù)包括壓力、溫度和噴射速度由等熵流動(dòng)計(jì)算所得，即將每個(gè)網(wǎng)格點(diǎn)視為等截面流管，入口參數(shù)根據(jù)當(dāng)?shù)貕毫腿肟诳倝河傻褥亓鲃?dòng)關(guān)系計(jì)算。設(shè)燃燒室內(nèi)入口頭部的壓力為ph，燃料噴射總壓和總溫分別為p0和T0，壅塞條件下臨界壓力為

(1)

式中：γ為比熱比。根據(jù)ph和p0的大小關(guān)系，入口邊界的壓力p、溫度T和軸向速度v的賦值可以分為以下3種情況：

圖2 計(jì)算區(qū)域示意圖Fig.2 Computational domain

1)若ph≥p0，則新鮮燃料停止噴射，入口邊界設(shè)為固壁邊界；

2)若pcr

(2)

式中：R為氣體常數(shù)。

3)若ph

(3)

為了模擬真實(shí)點(diǎn)火過程，在燃料噴射條件下，通過一股熱射流切向噴入燃燒室完成點(diǎn)火過程。如圖2所示，初始時(shí)計(jì)算區(qū)域內(nèi)充滿環(huán)境大氣，計(jì)算開始后，先通過入口邊界向燃燒室內(nèi)填充燃料和氧化劑的混合物，當(dāng)形成可燃預(yù)混層后點(diǎn)火程序啟動(dòng)。點(diǎn)火位置位于左邊界靠近入口處，高度為0.01 m. 在一個(gè)極短的時(shí)間內(nèi)，一股p=20p0、T=20T0和周相速度u=1 200 m/s的熱射流從點(diǎn)火位置噴射進(jìn)入燃燒室，完成點(diǎn)火。

1.2 控制方程

由于兩相爆轟包含復(fù)雜的兩相流動(dòng)和燃燒過程，液滴微粒的描述采用兩流體模型[21-22]。為了簡(jiǎn)化問題，本文做出以下假設(shè)：

1)爆轟過程是二維、無黏的，忽略與壁面的熱傳導(dǎo)、熱輻射和摩擦力的影響；

2)液滴保持球形，溫度分布均勻，液滴間互不影響；

3)液滴通過蒸發(fā)和剝離變?yōu)闅鈶B(tài)，并且與其他氣體瞬間混合均勻。

根據(jù)以上假設(shè)，得到二維氣液兩相連續(xù)旋轉(zhuǎn)爆轟的控制方程：

(4)

(5)

(6)

(7)

(8)

(9)

(10)

式中：φ、ρ、u、v、E分別為體積分?jǐn)?shù)、密度、周向速度、軸向速度和總能，p為氣相壓力，下標(biāo)g代表氣相參數(shù)，下標(biāo)l代表液相參數(shù)；Yk(k=1，2，3，4，5)分別為汽油蒸汽、氧氣、二氧化碳、水蒸氣和氮?dú)獾馁|(zhì)量分?jǐn)?shù)；Cv為汽油液滴的定容比熱；Id為單位體積的液滴質(zhì)量變化率[23-24]，

(11)

|vg-vl|=[(ug-ul)2+(vg-vl)2]1/2，

(12)

(11)式中第1項(xiàng)是液滴的剝離項(xiàng)，第2項(xiàng)是液滴的蒸發(fā)項(xiàng)，r為液滴半徑，(11)式已在多篇文獻(xiàn)中得到驗(yàn)證[21-22]，本文所用液滴初始半徑為40 μm，N為單位體積內(nèi)的液滴數(shù)，μg和μl分別為氣相和液相的黏性系數(shù)，vg和vl分別為氣相和液相的速度矢量，L為燃料液滴的蒸發(fā)潛熱，λ為氣體的熱傳導(dǎo)系數(shù)，Nu為Nusselt數(shù)，

Nu=2+0.6Re1/2Pr1/3，

(13)

Re為雷諾數(shù)，

(14)

Pr為普朗特?cái)?shù)；ωkp為每個(gè)組分的生成速率；ωkc為每個(gè)組分的消耗速率；Qc為汽油蒸氣燃燒的釋放熱；Qd為單位體積內(nèi)氣液兩相之間的交換熱，

Qd=4πr2NλNu(Tg-Tl)/2r；

(15)

Fx、Fy分別為氣體作用于單位體積混合物中液滴的拖拽力[21-22]，

Fx=0.5Nπr2CDρg|vg-vl|(ug-ul),

(16)

Fy=0.5Nπr2CDρg|vg-vl|(ug-ul),

(17)

CD為拖曳力系數(shù)，

(18)

本文采用辛烷的一步化學(xué)反應(yīng)代替汽油燃燒過程[10-11]，這種方法可以較為準(zhǔn)確地反映燃燒室內(nèi)組分的變化過程，并且相比于詳細(xì)的基元化學(xué)反應(yīng)可以節(jié)約大量時(shí)間。辛烷液滴經(jīng)蒸發(fā)和剝離后變?yōu)闅鈶B(tài)辛烷參加化學(xué)反應(yīng)，氣態(tài)辛烷的一步化學(xué)反應(yīng)式為

C8H18+12.5O2+N2→8CO2+9H2O+N2，

(19)

(19)式中的5個(gè)氣相組分分別是氣態(tài)辛烷、氧氣、氮?dú)?、二氧化碳、水蒸氣。根?jù)Arrhenius定律，汽油蒸汽的消耗速率為

(20)

式中：A為化學(xué)反應(yīng)指前因子；n1和n2分別為反應(yīng)級(jí)數(shù)；Ea為活化能；R0為通用氣體常數(shù)。其余組分的質(zhì)量變化速率可根據(jù)該反應(yīng)公式的當(dāng)量系數(shù)計(jì)算得到。

2 計(jì)算方法

2.1 CE/SE方法

CE/SE方法是一種用于求解雙曲型偏微分方程的數(shù)值計(jì)算方法，尤其適用于求解包含強(qiáng)間斷的問題。它將時(shí)間和空間統(tǒng)一進(jìn)行處理，設(shè)立SE和CE，保證了計(jì)算格式在計(jì)算域內(nèi)滿足物理守恒。

根據(jù)Chang[12]對(duì)時(shí)間和空間統(tǒng)一處理的思想，在CE的各個(gè)中心上對(duì)SE進(jìn)行泰勒展開，分別計(jì)算CE各個(gè)邊界面上時(shí)間- 空間密度矢量的積分通量，通過對(duì)(1)式運(yùn)用高斯定理, 可以得到非結(jié)構(gòu)化三角形網(wǎng)格CE/SE方法的計(jì)算格式：

(21)

(22)

(23)

(24)

(25)

(26)

(27)

兩相爆轟包含了復(fù)雜的源項(xiàng)，本文對(duì)于源項(xiàng)的處理方法是，將源項(xiàng)和流動(dòng)控制方程解耦，先凍結(jié)源項(xiàng)求解流場(chǎng)參數(shù)，然后不考慮流場(chǎng)的影響，采用4階龍格- 庫(kù)塔法求解源項(xiàng)，4階龍格- 庫(kù)塔法的時(shí)間步長(zhǎng)取為

(28)

式中：NR-K取為10.

在計(jì)算效率方面，對(duì)非結(jié)構(gòu)化三角形網(wǎng)格CE/SE方法進(jìn)行并行化處理。基于直接控制共享內(nèi)存式并行編程的應(yīng)用程序接口OpenMP并行化方法設(shè)計(jì)并行優(yōu)化方案，將計(jì)算程序?qū)崿F(xiàn)并行化，節(jié)省了大量時(shí)間。

2.2 方法驗(yàn)證

平面激波反射算例是檢驗(yàn)二維算法的一個(gè)經(jīng)典算例。該算例初始時(shí)在平面的左上角給定一個(gè)與上邊界夾角為29°的懸掛斜激波，激波左側(cè)來流參數(shù)為ρ=1.0，u=2.9，v=0，p=1/1.4，激波右側(cè)氣體參數(shù)ρ=1.7，u=2.619 3，v=-0.528 2，p=1.528 2.ρ、u、v、p在此算例中分別是無量綱化的密度，x軸、y軸方向的速度以及壓力。左邊界、上邊界采用來流邊界，下邊界為固壁邊界，右邊界為自由出口邊界。

本文用于計(jì)算激波反射算例的計(jì)算區(qū)域如圖3所示，邊界上劃分三角形網(wǎng)格邊長(zhǎng)為2 mm. 該區(qū)域共劃分為22 722個(gè)三角形，計(jì)算時(shí)間為無量綱時(shí)間2.5. 圖4給出了使用非結(jié)構(gòu)化三角形網(wǎng)格CE/SE方法計(jì)算所得的激波反射云圖和等值線圖，圖5為y=0.5處的壓力系數(shù)Cp分布圖和精確解的對(duì)比，Cp=2(p/p∞-1)/(γMa2)，p∞為大氣壓力，此處設(shè)置為0.1 MPa，Ma為馬赫數(shù)。從圖5可以看出，穩(wěn)定后的入射和反射激波具有較高的分辨率，壓力系數(shù)基本與精確解吻合，證明了該方法在使用2 mm三角形網(wǎng)格時(shí)具有準(zhǔn)確捕捉激波的能力。

圖3 非結(jié)構(gòu)化三角形網(wǎng)格Fig.3 Unstructured triangular meshes

圖4 壓力系數(shù)Cp云圖和線圖Fig.4 Contour and line graph of Cp

圖5 中線y=0.1 m上壓力系數(shù)分布Fig.5 Distribution of Cp for y=0.1 m

3 計(jì)算結(jié)果與討論

采用非結(jié)構(gòu)化CE/SE方法對(duì)氣液兩相RDE燃燒室內(nèi)流場(chǎng)進(jìn)行了數(shù)值計(jì)算，研究進(jìn)氣總壓和當(dāng)量比對(duì)兩相RDW的影響，分析不同模態(tài)RDE的工作特性。

3.1 進(jìn)氣總壓對(duì)兩相RDW的影響

進(jìn)氣總壓對(duì)RDW的傳播過程具有較大的影響。在恰當(dāng)量比下，分別計(jì)算了進(jìn)氣總壓為0.5 MPa、0.6 MPa、0.7 MPa、0.8 MPa時(shí)兩相RDW的傳播特性，如表1所示。

在計(jì)算RDE的總推力時(shí)，采用火箭式發(fā)動(dòng)機(jī)公式，靜推力即等于總推力，燃燒室推力為

(29)

式中：ρ、v、po分別為出口處的燃?xì)饷芏?、速度和壓力。此時(shí)基于燃料的比沖Isp為

(30)

表1 RDE各參量隨進(jìn)氣總壓變化統(tǒng)計(jì)表Tab.1 Parameters of RDE with the increase in total inlet pressure

(31)

為了計(jì)算推力等參數(shù)，假設(shè)計(jì)算模型有一個(gè)0.01 m的厚度，即本文所求推力、比沖等參量，均是基于燃燒室厚度為0.01 m、頭部周長(zhǎng)為0.3 m的RDE模型。在改變進(jìn)氣總壓條件下，各算例獲得穩(wěn)定傳播的爆轟波時(shí)間各不相同，但最終均形成了穩(wěn)定的RDW. 在各算例燃燒室內(nèi)流場(chǎng)進(jìn)入穩(wěn)定傳播階段后，選取離散的若干時(shí)刻來計(jì)算流場(chǎng)的推進(jìn)性能物理量，最后以各時(shí)刻的平均值作為最后的結(jié)果。

計(jì)算所得RDE傳播特性如表1所示，爆轟壓力選取20周內(nèi)RDW的波動(dòng)范圍。在當(dāng)量比為1的條件下，僅在進(jìn)氣總壓為0.5 MPa時(shí)燃燒室內(nèi)形成了穩(wěn)定傳播的同向雙波RDW，如圖6所示。其余各算例均形成了同向三波RDW. 在此僅給出進(jìn)氣壓力為0.6 MPa時(shí)同向三波模態(tài)RBW，如圖7所示。圖7中的三波模態(tài)RBW不穩(wěn)定性較高，每個(gè)波頭的強(qiáng)度和速度不同，因此也導(dǎo)致每2個(gè)波頭之間的間隔不一致，但3個(gè)波頭仍然可以保持動(dòng)態(tài)平衡和自持傳播。對(duì)比同向雙波和同向三波，三波模態(tài)的進(jìn)氣壓力要高于雙波模態(tài)，但是雙波模態(tài)的波頭高度仍然明顯高于三波模態(tài)。

圖6 同向雙波RDWFig.6 Double-wave RDW

圖7 同向三波RDWFig.7 Triple-wave RDW

(32)

從表1中波速標(biāo)準(zhǔn)差可以看出，傳播速度的不穩(wěn)定性與爆轟波壓力波動(dòng)的不穩(wěn)定性一致，即爆轟波的壓力波動(dòng)越大，其傳播速度的波動(dòng)也越大。

鄭權(quán)等[6]通過試驗(yàn)測(cè)量所得汽油與富氧空氣RDW傳播波速為1 022.2～1 171.8 m/s，通過數(shù)值計(jì)算所得氣兩相RDW傳播速度為1 204.01～1 265.82 m/s，與試驗(yàn)結(jié)果較為接近，爆轟波壓力和溫度峰值耦合，流場(chǎng)結(jié)構(gòu)與Bykovskii通過數(shù)據(jù)補(bǔ)償法監(jiān)測(cè)到的兩相旋轉(zhuǎn)爆轟流場(chǎng)結(jié)構(gòu)[2]定性一致(見圖1)，計(jì)算結(jié)果與試驗(yàn)研究結(jié)果吻合較好。

3.2 當(dāng)量比對(duì)兩相RDW的影響

為了研究當(dāng)量比對(duì)RDW的影響，在進(jìn)氣總壓為0.5 MPa條件下，分別計(jì)算當(dāng)量比為0.8、1.0、1.2、1.4共4種情況下RDW的傳播特性。計(jì)算所得RDE性能參數(shù)隨當(dāng)量比變化如表2所示。

表2 RDE各參量隨當(dāng)量比變化統(tǒng)計(jì)表Tab.2 Parameters of RDE with increase in equivalent ratio

隨著當(dāng)量比的增加，燃燒室內(nèi)形成了不同模態(tài)的RDW，在當(dāng)量比為0.8時(shí)獲得了同向單波模態(tài)的RDW，逐步提高當(dāng)量比，分別在當(dāng)量比為1.0時(shí)形成了同向雙波模態(tài)RDW，在當(dāng)量比為1.2和1.4時(shí)均形成了同向三波模態(tài)RDW. 由此可見，當(dāng)量比對(duì)RDW的波頭數(shù)影響較大，當(dāng)量比越低越容易形成單波模態(tài)，當(dāng)量比越高、越容易產(chǎn)生多個(gè)波頭。隨著當(dāng)量比的提升，多波頭模態(tài)的壓力峰值波動(dòng)范圍也逐漸擴(kuò)大，表明多波頭模態(tài)每個(gè)波頭的強(qiáng)度并不均勻，其強(qiáng)度仍處于動(dòng)態(tài)變化中。對(duì)比單波、雙波和三波的波速可以發(fā)現(xiàn)，隨著當(dāng)量比的提高，爆轟波的波速略有提升，但是，單波模態(tài)所形成的爆轟波速度略高于雙波和三波。

圖8為當(dāng)量比為0.8時(shí)所獲得的同向單波模態(tài)監(jiān)測(cè)點(diǎn)(0.001,0.001)處壓力和溫度曲線圖。由圖8可見，監(jiān)測(cè)點(diǎn)處溫度和壓力規(guī)律波動(dòng)，表明發(fā)動(dòng)機(jī)處于穩(wěn)定工作中。壓力峰值約為5.4 MPa，溫度峰值約為2 500 K，二者的上升沿非常陡峭并且緊密耦合在一起，溫度峰值到達(dá)時(shí)間略延遲于壓力峰值到達(dá)時(shí)間，二者的差距在10%以內(nèi)，符合液態(tài)燃料爆轟的典型特征。

圖8 單波模態(tài)監(jiān)測(cè)點(diǎn)處壓力和溫度Fig.8 Pressure and temperature of single wave at the monitoring point

圖9為當(dāng)量比為0.8時(shí)所獲得的單波模態(tài)RDW. 由于氣液兩相RDW傳播速度較慢，進(jìn)氣速度對(duì)燃燒室流場(chǎng)結(jié)構(gòu)的影響也更為顯著。由圖9(a)中可見斜激波角度較小，幾乎與爆轟波平行。在斜激波作用下，經(jīng)過斜激波的燃?xì)庵芟蛩俣确较虬l(fā)生了突躍變化，在斜激波下游和爆轟波上端產(chǎn)生一個(gè)渦，將爆轟波上端部分燃?xì)饩砣肫湎掠巍Ｍ瑫r(shí)由于氣相密度遠(yuǎn)低于液相，卷入斜激波下游的燃?xì)獯蟛糠譃槲慈伎諝?見圖9(c))，由于缺少燃料熱能的釋放，導(dǎo)致此處溫度較低，即為圖9(b)中的爆轟波上端溫度較低區(qū)域。

圖9 單波模態(tài)RDWFig.9 Single wave RDW

圖10為沿直線y=0.01 m處壓力、溫度和燃料液滴半徑的分布圖，由圖10可見，壓力峰值和溫度峰值之間存在2.4 cm的距離，燃料液滴半徑在溫度峰值處仍有21 μm. 圖11為L(zhǎng)iu等[7]使用二維5階WENO方法獲得的H2/Air單波模態(tài)RDW燃燒室頭部周向壓力和溫度分布曲線。由圖11可見，氣態(tài)燃料爆轟波的壓力峰值和溫度峰值之間的距離極小，二者完全耦合在一起，燃料燃燒釋放的熱量可直接支持旋轉(zhuǎn)爆轟波的傳播。對(duì)比圖10和圖11可以發(fā)現(xiàn)，液滴的蒸發(fā)和剝離過程延緩了液滴燃料的燃燒和熱量的釋放，導(dǎo)致了兩相爆轟波溫度峰值滯后壓力峰值2.4 cm，大量燃料燃燒釋放的熱量需要穿透2.4 cm的液滴層才得以支撐前導(dǎo)激波的傳播，因此溫度峰值和壓力峰值的不完全耦合也降低了兩相爆轟波的傳播速度。

圖10 壓力、溫度和液滴半徑分布圖(y=0.01 m)Fig.10 Distribution graph of pressure, temperature and droplet radius(y=0.01 m)

圖11 氣態(tài)燃料RDW壓力和溫度分布圖[7]Fig.11 Distribution graph of pressure and temperature of gaseous RDW[7]

3.3 不同模態(tài)RDE工作特性分析

為了分析不同模態(tài)RDW的壓力波動(dòng)特性，將表1和表2中所示壓力峰值波動(dòng)范圍繪制成圖12. 從圖12(a)中可見，隨著進(jìn)氣總壓的升高，同向三波模態(tài)RDW壓力波動(dòng)范圍總體呈現(xiàn)出增大的趨勢(shì)，但是同向雙波模態(tài)RDW壓力峰值波動(dòng)范圍明顯低于同向三波模態(tài)。圖12(b)反映出RDW壓力波動(dòng)范圍隨當(dāng)量比增大而增大的趨勢(shì)，同時(shí)也反映出單波模態(tài)的壓力峰值波動(dòng)最小，雙波居中，三波模態(tài)壓力波動(dòng)最大的規(guī)律，表明多波頭模態(tài)雖然可以維持穩(wěn)定傳播，但是每個(gè)波頭強(qiáng)度并不恒定，爆轟波傳播的不穩(wěn)定性較大。

圖12 不同模態(tài)RDW壓力波動(dòng)范圍圖Fig.12 Range of detonation pressure in different modalities

圖13給出了進(jìn)氣總壓為0.6 MPa時(shí)三波模態(tài)監(jiān)測(cè)點(diǎn)處的壓力曲線。由圖13可見，3個(gè)爆轟波頭的強(qiáng)度并不固定，呈現(xiàn)出規(guī)律變化的趨勢(shì)，每個(gè)波頭的壓力都在高、中、低不斷變換。以3.0 ms波頭1為例，由于此時(shí)波頭強(qiáng)度較低，速度落后于波頭2，波頭前可燃預(yù)混層較高，爆轟波得到更多能量的支持，其強(qiáng)度不斷增加，逐漸達(dá)到頂峰。但在此過程中波頭1逐漸接近波頭2，其前端可燃預(yù)混層變小，爆轟波逐漸減弱，爆轟壓力開始降低。3個(gè)波頭共同維持一個(gè)動(dòng)態(tài)平衡，使得三波模態(tài)RDW能夠持續(xù)穩(wěn)定傳播。

圖13 三波模態(tài)監(jiān)測(cè)點(diǎn)處壓力曲線Fig.13 Pressure graph of triple waves at monitoring point

由于RDE通常是一個(gè)環(huán)形燃燒室，二維數(shù)值模擬研究雖然可以反映RDW的傳播特性，但是想要深入探究曲率效應(yīng)和燃燒室寬度對(duì)兩相旋轉(zhuǎn)爆轟波流場(chǎng)結(jié)構(gòu)和傳播不穩(wěn)定性的影響，后續(xù)研究仍然需要進(jìn)行三維數(shù)值仿真以及與試驗(yàn)研究驗(yàn)證的對(duì)照。

圖14為推力和燃料比沖隨當(dāng)量比變化曲線圖。由圖14可見：發(fā)動(dòng)機(jī)的推力隨著當(dāng)量比提升而提升，但提升幅度較小；發(fā)動(dòng)機(jī)比沖隨當(dāng)量比的增加呈現(xiàn)出下降的趨勢(shì)。

圖14 推力和燃料比沖隨當(dāng)量比變化曲線Fig.14 Thrust and specific impulse versus equivalent ratio

圖15為按照表2順序繪制的不同模態(tài)RDE推力曲線圖。由圖15可見，點(diǎn)火之后發(fā)動(dòng)機(jī)推力均出現(xiàn)一個(gè)突躍的上升，而后迅速降低，經(jīng)歷數(shù)次振蕩之后逐漸穩(wěn)定。其中：?jiǎn)尾B(tài)進(jìn)入穩(wěn)定階段最快，僅用時(shí)1.81 ms；雙波模態(tài)次之，用時(shí)2.50 ms；三波模態(tài)最慢，用時(shí)2.82 ms和3.16 ms. 進(jìn)入穩(wěn)定階段之后，單波模態(tài)發(fā)動(dòng)機(jī)推力波動(dòng)最大，為0.981 kN，是其平均推力的42.84%，雙波模態(tài)推力波動(dòng)為0.387 kN，占其平均推力的16.33%，三波模態(tài)發(fā)動(dòng)機(jī)推力波動(dòng)最小，為0.179 kN和0.231 kN，是其平均推力的7.28%和9.31%. 因此，多波頭模態(tài)有利于降低RDE的推力波動(dòng)，增加RDE的工作穩(wěn)定。

圖15 不同模態(tài)推力曲線Fig.15 Thrust graphs of different modalities

4 結(jié)論

本文采用非結(jié)構(gòu)化三角形網(wǎng)格CE/SE方法對(duì)氣液兩相RDW進(jìn)行了數(shù)值模擬，建立了氣液兩相連續(xù)旋轉(zhuǎn)爆轟計(jì)算模型，獲得了單波模態(tài)、同向雙波模態(tài)以及同向三波模態(tài)的氣液兩相RDW，計(jì)算了不同進(jìn)氣總壓和不同當(dāng)量比下的RDW傳播過程，分析了不同模態(tài)下RDW的壓力和推力特性。得出以下主要結(jié)論：

1)非結(jié)構(gòu)化CE/SE方法可以準(zhǔn)確捕捉激波、爆轟波等強(qiáng)間斷，具有良好的應(yīng)用前景；數(shù)值計(jì)算所得兩相RDW傳播波速與試驗(yàn)結(jié)果吻合較好，流場(chǎng)結(jié)構(gòu)與試驗(yàn)結(jié)果定性一致，驗(yàn)證了數(shù)值模擬研究的準(zhǔn)確性。

2)液滴的蒸發(fā)和剝離過程，延緩了液滴燃料的燃燒，導(dǎo)致兩相爆轟波壓力峰值和溫度峰值的不完全耦合，降低了兩相爆轟波的傳播速度。

3)當(dāng)量比對(duì)RDW的傳播模態(tài)影響較大，當(dāng)量比越低越容易形成單波模態(tài)，當(dāng)量比越高越容易產(chǎn)生多個(gè)波頭。

4)單波模態(tài)的爆轟波壓力和傳播速度波動(dòng)最小，雙波模態(tài)居中，三波模態(tài)壓力和傳播速度波動(dòng)最大；其推力波動(dòng)規(guī)律與壓力波動(dòng)規(guī)律相反，單波模態(tài)推力波動(dòng)最大，雙波模態(tài)次之，三波模態(tài)推力波動(dòng)最小。