噴管對水下爆轟燃?xì)馍淞鬟^程影響的數(shù)值模擬

邵新科,康 楊,胡錦橋,李 寧,黃孝龍,邊哲民,高榮維,翁春生

(南京理工大學(xué) 1.瞬態(tài)物理國家重點實驗室;2.能源與動力工程學(xué)院,江蘇 南京 210094)

爆轟是一種超聲速燃燒現(xiàn)象,爆轟燃燒近似等容燃燒,與爆燃相比,其具有熱效率高的優(yōu)勢[1]。更高的熱循環(huán)效率使得脈沖爆轟發(fā)動機[2-3]、旋轉(zhuǎn)爆轟發(fā)動機[4-5]和斜爆轟發(fā)動機[6-7]等基于爆轟循環(huán)的發(fā)動機成為空天推進領(lǐng)域的研究熱點之一[8]。脈沖爆轟水沖壓發(fā)動機(pulse detonation hydro-ramjet,PDH)是一種基于爆轟燃燒,依靠周期性爆轟波和高溫高壓爆轟燃?xì)馀c水介質(zhì)相互作用產(chǎn)生推力的新型水下推進裝置,具有結(jié)構(gòu)簡單、比沖大、燃燒效率高等諸多優(yōu)勢[9],在水下推進領(lǐng)域有著廣闊的應(yīng)用前景。PDH將爆轟發(fā)動機的應(yīng)用范疇由空天推進領(lǐng)域拓展至水下推進領(lǐng)域,逐漸成為水下新型動力系統(tǒng)研究的前沿與熱點。

針對水下爆轟推進的研究最早由俄羅斯科學(xué)院的FROLOV教授在2013年展開[10]。文獻[11-12]基于兩相可壓縮有黏流方程,開展了激波從氣體向液相水或含氣水中傳播過程的二維數(shù)值模擬和實驗研究,并指出激波沖擊可壓含氣水的動量轉(zhuǎn)換率要高于激波沖擊不可壓純水的轉(zhuǎn)換率,驗證了以含氣水為工質(zhì)的PDH的可行性。FROLOV等[9,13]搭建了PDH的實驗樣機,以實驗方式得到了采用汽油燃料的PDH的比沖：工作頻率為20 Hz時,比沖為370 s;工作頻率為1 Hz時,比沖為1 200 s。并指出通過實驗測量的比沖超過了液體火箭發(fā)動機。

國內(nèi)有關(guān)水下爆轟的研究較少,起步也較晚。由于水下爆轟燃?xì)馍淞鞯难芯渴撬卤Z發(fā)動機研發(fā)的基礎(chǔ),因此國內(nèi)的相關(guān)研究更側(cè)重水下爆轟過程的基礎(chǔ)流場機理。其中劉威等[14-16]利用高速攝影技術(shù)和PIV技術(shù),以實驗手段研究了水下爆轟燃?xì)馀菪螒B(tài)變化規(guī)律,應(yīng)用小波分析方法獲取了水下爆轟壓力場的時頻特性。同時也對水下爆轟燃?xì)馍淞鬟M行了數(shù)值模擬研究,揭示了水中激波和燃?xì)馍淞饕鸬膲毫_動的傳播特點和爆轟燃?xì)馀莸男螒B(tài)成長規(guī)律等流場特征。侯子偉[17-18]針對水下爆轟流場的激波傳播過程和爆轟燃?xì)馀菪螒B(tài)等流場特征進行了數(shù)值模擬和實驗研究,得到了水下爆轟燃?xì)馀莸男螒B(tài)成長規(guī)律和水中激波的傳播規(guī)律,并指出使激波在爆轟管零度軸線方向保持較高強度可以提升發(fā)動機的性能,同時還利用歐拉雙流體模型和時空守恒元與求解元方法對水下雙管爆轟的近場復(fù)雜波系發(fā)展過程進行了三維數(shù)值模擬,研究了點火延遲、管間距和填充條件等因素對壓力場特征和激波傳播過程的影響。周幃等[19]基于Mixture模型,開展了水下爆轟外流場的二維軸對稱數(shù)值模擬,得到了出口爆轟波與外界水介質(zhì)作用下誘發(fā)的壓力場、速度場和溫度場的變化規(guī)律。

綜上所述,現(xiàn)有研究初步探究了水下爆轟過程的相關(guān)典型物理現(xiàn)象的機理,并進行了簡單樣機的實驗,驗證了將爆轟發(fā)動機用于水下推進的可行性。尾噴管是發(fā)動機研究的重要組成部分,尾噴管的加裝對發(fā)動機的出口流場特性和推進性能均有重要影響,但目前有關(guān)加裝噴管的水下爆轟燃?xì)馍淞鞯南嚓P(guān)研究較少。本文采用計算流體力學(xué)手段,基于VOF多相流模型建立了水下爆轟過程的二維軸對稱數(shù)值計算模型,并對無噴管和采用擴張噴管、收斂噴管的爆轟管水下爆轟過程的內(nèi)外流場進行數(shù)值模擬研究,重點分析水下爆轟過程中噴管構(gòu)型對透射激波、反射激波的傳播過程和爆轟燃?xì)馀菪螒B(tài)尺度發(fā)展過程的影響。

1 數(shù)值計算模型

1.1 基本假設(shè)與數(shù)值方法

水下爆轟流場是包含激波在氣液兩相中傳播的反應(yīng)流、多相流,流場結(jié)構(gòu)十分復(fù)雜。為了簡化模型,作如下假設(shè)：

①對于氣相,不考慮分子間作用力,采用理想氣體狀態(tài)方程;對于液相,采用可壓縮液體的Tait狀態(tài)方程。

②不考慮表面張力,重力;不考慮兩相交界面的相變和質(zhì)量交換。

③爆轟管管壁與內(nèi)外流場均無熱量交換,熱通量處處為0。

根據(jù)以上假設(shè),基于VOF多相流模型的控制方程如下。

連續(xù)性方程：

(1)

式中：αq,ρq,vq分別為多相流的第q相的體積分?jǐn)?shù)、密度和速度矢量。

動量方程：

(2)

式中：ρ為體積分?jǐn)?shù)平均密度,μ為體積平均黏性。

(3)

(4)

能量方程：

(5)

(6)

(7)

式中：hq為第q相的焓。

本文采用VOF多相流模型進行爆轟管內(nèi)外氣-液兩相流動的耦合計算,利用RNGk-ε兩方程湍流模型來封閉RANS方程的雷諾應(yīng)力項,并采用標(biāo)準(zhǔn)壁面函數(shù)。湍流模型的輸運方程詳見文獻[20]。壓力與速度耦合采用分離式PISO(pressure-implicit with splitting of operators,PISO)算法,采用PRESTO!(pressure staggering option)格式對壓力項進行空間離散,其余各項均采用QUICK(quadratic upstream interpolation for convective kinematics,QUICK)格式進行空間離散,時間項采用一階隱式格式進行離散。

1.2 計算模型與邊界條件

圖1給出了采用擴張噴管的爆轟管的二維軸對稱計算域和邊界條件。計算域以爆轟管的管壁為界,分為內(nèi)流場和外流場。內(nèi)流場長度AB=1 000 mm,半徑AD=15 mm。外流場寬度KM=600 mm,長度JK=1 200 mm,CM=600 mm。計算中所用的收斂、擴張噴管的收斂或擴張角度均為30°,長度均為15 mm。

圖1 計算域與邊界條件的示意圖Fig.1 Schematic diagram of the computational domain and boundary conditions

初始條件：初始時刻爆轟管外流場充滿了靜止的水介質(zhì),壓力為101 325 Pa,溫度為300 K。內(nèi)流場充滿了當(dāng)量比為1的甲烷和氧氣預(yù)混氣,填充壓力為0.3 MPa,溫度為300 K。計算域左側(cè)AEFD區(qū)域為點火區(qū),寬度為5 mm,溫度為3 000 K,壓力為3 MPa。

內(nèi)流場的氣相介質(zhì)為多相流的第一相,采用有限反應(yīng)速率模型,反應(yīng)速率由Arrhenius 公式控制,即：

(8)

式中：A為指前因子;Ea為活化能;R為理想氣體常數(shù);T為絕對溫度;b為溫度指數(shù),b=0。采用甲烷-氧氣兩步反應(yīng),相關(guān)參數(shù)如表1所示。外流場的液相介質(zhì)為多項流的第二相,不考慮化學(xué)反應(yīng)。

表1 甲烷-氧氣兩步反應(yīng)Table 1 Methane/oxygen two-step reaction

邊界條件：圖1中爆轟管壁處,采用無滑移絕熱固壁邊界條件;在爆轟管軸線處,采用軸對稱邊界條件。為避免透射激波傳播至壓力出口邊界時反射形成假波,在外流場出口處采用非反射壓力出口邊界,分兩種情況處理：①當(dāng)出口為亞聲速時,壓力出口邊界處的靜壓恒為101 325 Pa;②當(dāng)壓力出口邊界局部出現(xiàn)超聲速流動時,則該處的靜壓由流場內(nèi)部外推得到。壓力出口邊界處的所有其他守恒變量均由流場內(nèi)部外推得到。

爆轟管口處的氣-液交界面BC的處理方式如下：由于爆轟波是超聲速燃燒波,波后擾動無法影響波前狀態(tài),因此在爆轟波傳出爆轟管之前,僅對內(nèi)流場進行計算,這一階段將BC作為固壁處理。當(dāng)爆轟波傳播至管口后,BC不再采用固壁條件,而是成為流場內(nèi)部的一部分,隨即開始內(nèi)-外流場的耦合計算,以爆轟波到達管口BC處作為內(nèi)-外流場耦合計算開始的0時刻。

為了捕捉外流場的壓力場特征,以圖1中爆轟管口中心C點為極坐標(biāo)原點,以外軸線lCM為零度方向建立局部極坐標(biāo)系。分別在極角θ=0°,15°,30°,45°,60°,75°,90°,極徑r=100 mm,200 mm,300 mm,400 mm,500 mm處預(yù)先布置35個壓力監(jiān)測點。

1.3 模型驗證

點火后爆轟波傳播過程的軸線壓力分布如圖2所示,其中縱坐標(biāo)為爆轟管軸線壓力幅值p與該工況下的理論C-J爆壓pC-J的比值。由圖2可見,在點火區(qū)域的作用下,爆轟管內(nèi)形成了穩(wěn)定傳播的爆轟波。表2將內(nèi)流場爆轟參數(shù)計算結(jié)果與NASA CEA軟件在相同填充條件下的理論C-J爆轟參數(shù)進行比較,由表2可知,計算結(jié)果與理論C-J爆轟參數(shù)相對誤差較小,驗證了內(nèi)流場模型的可靠性和計算的準(zhǔn)確性。

圖2 點火后不同時刻爆轟管軸線的壓力分布Fig.2 Pressure distribution along the detonation tube axis

表2 數(shù)值模擬計算結(jié)果與理論值的對比Table 2 Comparison of computational results and C-J theoretical values

圖3 氣泡無量綱面積的實驗與計算結(jié)果對比Fig.3 Comparison of bubble dimensionless area between experimental and computational results

以無噴管的工況為例,對外流場計算域進行網(wǎng)格無關(guān)性驗證和時間步長無關(guān)性驗證。圖4(a)給出了采用尺寸為0.8 mm,1.0 mm和2.0 mm的四邊形結(jié)構(gòu)化網(wǎng)格進行計算時,中心軸線上距離管口10 mm處監(jiān)測點的壓力曲線?？梢?種網(wǎng)格均能有效捕捉透射激波,且激波壓力幅值幾乎不受網(wǎng)格大小影響,進一步增加網(wǎng)格數(shù)量的意義不大。圖4(b)對比了采用3種時間步長(Δt)進行計算時,監(jiān)測點處的壓力曲線?？梢?種時間步長均可有效捕捉激波的傳播過程,監(jiān)測點壓力曲線幾乎不受時間步長的影響,進一步縮小時間步長收效甚微。綜合考慮計算精度和計算成本,最終選擇1.0 mm的四邊形結(jié)構(gòu)化網(wǎng)格劃分外流場計算域,時間步長取Δt=1×10-7s。計算中內(nèi)流場網(wǎng)格總數(shù)約為13.5萬,外流場網(wǎng)格總數(shù)約為54.7萬。

圖4 監(jiān)測點處的壓力變化曲線Fig.4 Pressure evolutions of measurement points

2 計算結(jié)果及其分析

2.1 噴管構(gòu)型對透射與反射激波的影響

在水下爆轟過程中,爆轟波傳播至爆轟管出口時失去化學(xué)反應(yīng)的能量支持,解耦為不能自持的激波,并與環(huán)境水介質(zhì)發(fā)生碰撞。圖5給出了爆轟波碰撞氣-液交界面后t=0.05 ms時刻,加裝擴張噴管、收斂噴管和無噴管的爆轟管流場壓力云圖。

圖5 各式爆轟管流場壓力云圖Fig.5 Pressure contour of detonation tube with different nozzles at t=0.05 ms

由于爆轟管外環(huán)境水介質(zhì)的聲阻抗遠(yuǎn)大于管內(nèi)高溫高壓的爆轟燃?xì)?流場中均形成了向下游傳播的透射激波與向上游傳播的反射激波。隨后高溫高壓爆轟燃?xì)獬掷m(xù)向水中噴射,形成水下高速氣體射流。向下游傳播的透射激波的強度和指向性會直接影響沖擊波與水介質(zhì)的動量交換效率;向上游傳播的反射激波會提升爆轟管的頭部壓力,并最終作用在爆轟管的推力壁上,增大平均推力。

采用擴張噴管的爆轟管不同時刻軸線上的壓力分布如圖6所示,圖中以x=1.0 m為界,左側(cè)為管內(nèi),右側(cè)為管外軸線壓力分布。和無噴管爆轟管水下單次爆轟過程[17]相似：激波碰撞氣-液交界面后初期,管內(nèi)靠近管口附近的爆轟燃?xì)馐芊瓷浼げǘ螇嚎s,出現(xiàn)了短暫高于C-J爆壓的高壓區(qū)。透射激波和反射激波的壓力幅值均隨時間逐漸衰減,且同一時刻管內(nèi)反射激波壓力幅值大于管外透射激波。t=0.1 ms時反射激波壓力幅值為10.81 MPa,略大于C-J爆轟壓力pC-J=9.06 MPa,而透射激波的幅值僅為1.89 MPa,遠(yuǎn)小于前者。由于環(huán)境水介質(zhì)的慣性和密度極大,因此燃?xì)馀莩跗谂蛎浭茏?氣-液交界面的發(fā)展過程遠(yuǎn)遠(yuǎn)滯后于透射激波引起的壓力擾動在水介質(zhì)中的傳播過程。以t=0.2 ms為例,此時透射激波已傳播至x=1.298 m處,而氣-液交界面僅發(fā)展至x=1.009 m處,遠(yuǎn)遠(yuǎn)滯后于前者。

圖6 采用擴張噴管爆轟管的內(nèi)外軸線壓力分布圖Fig.6 Pressure distribution along the axis of detonation tube with diverging nozzle at various time

透射激波的壓力幅值隨傳播距離的變化如圖7所示。圖中,d為到管口的距離,pTS為透射激波的壓力幅值。

圖7 透射激波壓力幅值變化曲線Fig.7 Amplitudes of transmitted shock waves at different distance on axis with different nozzles

由圖7可見,各型爆轟管的透射激波在離開管口后均迅速衰減,但采用收斂噴管時透射激波強度明顯更小。t=0.1 ms時,無噴管和采用擴張噴管時透射激波壓力幅值分別為2.510 MPa和2.784 MPa,而采用收斂噴管時透射激波壓力幅值僅為1.278 MPa。這是因為透射激波在收斂噴管內(nèi)傳播時會不斷碰撞收斂噴管內(nèi)壁并向管內(nèi)反射激波,這些二次反射激波消耗了部分透射激波的能量,造成了透射激波強度的減小。

圖8給出了透射激波的指向性分布。圖中不同角度的激波壓力幅值均以該工況下的零度軸線方向進行歸一化。由圖8可知,透射激波在傳播初期均具有明顯的軸線方向的指向性。當(dāng)加裝擴張噴管時,透射激波沿軸線方向的指向性更強,擴張噴管對透射激波有一定的“聚焦”作用。而收斂噴管的出口尺寸小,出口處更相似于“點”擾動源,采用收斂噴管時透射激波的壓力幅值被削弱,透射激波的能量分布更加分散,透射激波的指向性不顯著。對比圖8(a)和圖8(b)可知：隨著透射激波在水中傳播距離的增加,噴管對透射激波的指向性的影響逐漸減弱。所有透射激波均隨著傳播距離的增加而逐漸失去了沿零度軸線的指向性,并逐漸發(fā)展為近乎均勻傳播的球面波。

圖8 管口不同距離處的透射激波指向性分布Fig.8 Directivity of transmitted shock wave of detonation tube with different nozzles at different distance

反射激波的壓力幅值隨傳播時間的變化曲線如圖9所示。由圖9可見所有反射激波均隨著時間逐漸衰減。這種現(xiàn)象在激波碰撞氣-液界面后的前0.2 ms內(nèi)尤為明顯,由圖中曲線斜率可知前0.2 ms內(nèi)反射激波強度衰減速度更快,這是因為向上游傳播的反射激波與向下游傳播的泰勒膨脹波迎面相碰(圖6中A點),兩波相遇處軸向速度梯度巨大,因而造成了更大的壓力損失。對比圖9中3條曲線可知,反射激波的強度受噴管影響有限,但采用收斂噴管時,反射激波強度略大于其他兩種,這是因為收斂噴管內(nèi)壁的二次反射激波追趕并疊加增強了反射激波的強度。t=0.1 ms時,采用收斂噴管的反射激波壓力幅值為11.343 MPa,而無噴管與采用擴張噴管的反射激波壓力幅值分別為10.063 MPa和10.534 MPa,后兩者略小于前者。

圖9 反射激波壓力幅值變化曲線Fig.9 Amplitudes evolutions of reflected shock waves on detonation tube with different nozzles

2.2 噴管構(gòu)型對燃?xì)馍淞餍螒B(tài)的影響

在水中傳播的透射激波于t=1 ms時已傳出計算域,但爆轟波后的高溫高壓爆轟燃?xì)鈺^續(xù)向管外的水介質(zhì)中噴射,形成水下高溫高速燃?xì)馍淞?。隨著爆轟燃?xì)庵饾u進入水中,氣-液界面逐漸成長。氣-液界面的發(fā)展和爆轟燃?xì)馀莸男螒B(tài)變化也會影響下一循環(huán)爆轟波的出口邊界條件,因此燃?xì)馀莸男螒B(tài)特性是描述水下爆轟流場特性的重要參數(shù),具有重要的研究意義。

采用不同噴管時,爆轟燃?xì)馀莸男螒B(tài)演化過程如圖10～圖12所示,圖中黑色輪廓線為燃?xì)馀輾庀囿w積分?jǐn)?shù)φg=0.5的等值線。圖中矢量場為速度矢量場,箭頭方向為當(dāng)?shù)厮俣仁噶糠较?。由于出口爆轟燃?xì)獾妮S向速度分量遠(yuǎn)大于徑向速度分量,因此,圖中可見3種燃?xì)馀菥暂S向發(fā)展為主。爆轟燃?xì)庋刂淞髦魍ǖ老蛳掠伟l(fā)展,造成氣液兩相之間的裹挾與摻混,所以射流下游也形成了大尺度的渦。但相比之下,采用收斂噴管的爆轟管管口半徑更小,造成了管口處更大的雷諾數(shù),由于氣液兩相之間存在巨大的軸向速度梯度,爆轟燃?xì)馐軅?cè)邊水介質(zhì)的強剪切作用,迅速由均勻流轉(zhuǎn)變?yōu)榧羟辛?在其管口附近形成了更強的渦系。t=1 ms時,水介質(zhì)的剪切作用尤為顯著,燃?xì)庀騼蓚?cè)翻轉(zhuǎn)明顯,其氣泡呈“W”型,而后隨著氣-液交界面的逐漸成長,兩相之間的軸向速度梯度減小,兩側(cè)水介質(zhì)的剪切作用減弱,氣泡膨脹速度逐漸放慢,并逐漸轉(zhuǎn)變?yōu)闄E球形。采用擴張噴管和無噴管[14-17]爆轟管的氣泡形態(tài)變化較為相似,初期由于外流場水介質(zhì)與內(nèi)流場爆轟燃?xì)饩薮蟮膲毫Σ钆c密度差,氣泡因膨脹受阻而呈扁平狀,隨著管內(nèi)爆轟燃?xì)獬掷m(xù)噴出,逐漸發(fā)展為近大遠(yuǎn)小的橢球形,相較于采用收斂噴管的爆轟管,其管口出現(xiàn)更明顯的回流現(xiàn)象,燃?xì)馀菹蛏嫌畏D(zhuǎn)的距離更遠(yuǎn)。

圖10 加裝收斂噴管的爆轟管燃?xì)馀葑兓^程Fig.10 Bubble shape evolutions of detonation tube with converging nozzle

圖11 無噴管爆轟管燃?xì)馀葑兓^程Fig.11 Bubble shape evolutions of detonation tube with no nozzle

圖12 采用擴張噴管的爆轟管燃?xì)馀葑兓^程Fig.12 Bubble shape evolutions of detonation tube with diverging nozzle

為了定量描述氣泡的形態(tài)尺度特征,并考察噴管構(gòu)型對燃?xì)馀莅l(fā)展過程的影響,定義了兩個氣泡特征尺度：軸向發(fā)展尺度Ra和徑向發(fā)展尺度Rr,如圖13所示。其中軸向尺度Ra為爆轟管出口截面至爆轟燃?xì)馀葸吘壍妮S向距離,徑向尺度Rr為爆轟管中心軸線到燃?xì)馀葸吘壍膹较蚓嚯x。

圖14為燃?xì)馀萏卣鞒叨萊a和Rr隨時間的變化曲線。可見在t=6 ms之前,噴管構(gòu)型對燃?xì)馀莸妮S向特征尺度Ra的影響較小,但燃?xì)馀莸膹较虺叨萊r受噴管構(gòu)型影響較大,采用收斂噴管時徑向尺度Rr明顯更小,噴管對出口燃?xì)馍淞鞯膹较蚣s束作用遠(yuǎn)大于軸向。由圖14中曲線斜率的變化可知,3種燃?xì)馀莸腞a和Rr增長速度均隨著時間逐漸放慢。在t=6 ms之前,由于出口爆轟燃?xì)馍淞鞯妮S向速度分量遠(yuǎn)大于徑向分量,造成了氣泡的軸向和徑向成長速度規(guī)律不同：3種燃?xì)馀莸妮S向尺度成長速度均遠(yuǎn)大于徑向尺度的成長速度。以t=5 ms為例,此時無噴管和采用收斂噴管、采用擴張噴管的氣泡軸向尺度平均生長速度分別為13.4 m/s,10.1 m/s,14.6 m/s;徑向尺度平均生長速度分別為4.1 m/s,6.0 m/s,5.6 m/s。

圖14 加裝不同噴管的爆轟管的燃?xì)馀葺S向發(fā)展尺度和徑向發(fā)展尺度Fig.14 Ra and Rr of bubbles of detonation tube with different nozzles

3 結(jié)束語

本文基于VOF多相流模型,建立了水下爆轟過程的流場計算模型,并分別對采用擴張噴管、收斂噴管和無噴管爆轟管的水下爆轟過程進行了數(shù)值模擬,重點分析了噴管構(gòu)型對透射激波、反射激波的傳播過程和爆轟燃?xì)馍淞餍螒B(tài)演化過程的影響,得到以下結(jié)論：

①對比數(shù)值計算結(jié)果和實驗結(jié)果,文中所建立的水下爆轟數(shù)值計算模型可以相對準(zhǔn)確地描述水下爆轟初期的流動特征和流場特性。

②爆轟波解耦形成的激波沖擊氣-液交界面會形成向下游傳播的透射激波和向上游傳播的反射激波。反射激波與向下游傳播的泰勒膨脹波相撞,導(dǎo)致反射激波在前0.2 ms衰減最為劇烈。擴張噴管可以加強透射激波沿軸線方向的指向性,而收斂噴管會減弱透射激波的強度,增強向上游傳播的反射激波強度。

③爆轟管出口爆轟燃?xì)獾妮S向速度分量遠(yuǎn)大于徑向速度分量,水下爆轟初期燃?xì)馍淞鞯陌l(fā)展以軸向為主。爆轟燃?xì)庋刂淞髦魍ǖ老蛳掠伟l(fā)展,造成氣、液兩相之間的裹挾與摻混,在射流下游形成了大尺度的渦。燃?xì)馀莩跗诘妮S向和徑向發(fā)展速度均隨著時間逐漸減小,噴管對燃?xì)馀莸妮S向尺度影響較小,但收斂噴管能夠顯著抑制燃?xì)馀莸膹较虺叨取?/p>