考慮徑向加質(zhì)的固體火箭發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)整體穩(wěn)定性分析①

楊尚榮，劉佩進(jìn)，魏少娟，魏祥庚

(西北工業(yè)大學(xué)燃燒、熱結(jié)構(gòu)與內(nèi)流場(chǎng)重點(diǎn)實(shí)驗(yàn)室，西安 710072)

0 引言

大型分段固體火箭發(fā)動(dòng)機(jī)中的壓強(qiáng)振蕩主要是渦脫落現(xiàn)象引起的［1］。理論上，CFD(例如LES)可完整地模擬渦脫落與聲場(chǎng)的耦合過(guò)程，預(yù)示燃燒室內(nèi)壓力振蕩的頻率和幅值［2］。但實(shí)際預(yù)估時(shí)，由于發(fā)動(dòng)機(jī)幾何尺寸較大，工作時(shí)間又很長(zhǎng)(大于100 s)，對(duì)每一時(shí)刻的工作狀態(tài)進(jìn)行數(shù)值模擬計(jì)算是不現(xiàn)實(shí)的，因此有必要從穩(wěn)定性理論角度出發(fā)，對(duì)渦脫落的機(jī)理進(jìn)行研究。

縮比發(fā)動(dòng)機(jī)［3］和冷流實(shí)驗(yàn)［4－5］都證明單純表面渦脫落也可引起發(fā)動(dòng)機(jī)內(nèi)的壓強(qiáng)振蕩，并且其振蕩幅值強(qiáng)于由障礙渦脫落導(dǎo)致的振蕩幅值［6］。與由剪切不穩(wěn)定產(chǎn)生的障礙渦脫落相比，表面渦脫落被認(rèn)為是由固體推進(jìn)劑燃燒生成的徑向加質(zhì)流的“本質(zhì)不穩(wěn)定”決定的［7］。相似的是，它們都可以用流動(dòng)不穩(wěn)定性理論來(lái)研究并獲得不穩(wěn)定特征以及演化規(guī)律。

Casalis等［7］利用局部非平行方法討論了平面徑向加質(zhì)流的不穩(wěn)定模式，發(fā)現(xiàn)平均流的橫向分量對(duì)穩(wěn)定性結(jié)果有很大的影響。計(jì)算結(jié)果表明分別用主變量和勢(shì)函數(shù)表示的擾動(dòng)方程得到的結(jié)果會(huì)有差異［8－9］，說(shuō)明局部非平行方法是不一致的。文獻(xiàn)［10］同樣利用局部方法研究了圓柱構(gòu)型的徑向加質(zhì)流不穩(wěn)定特征，但其與實(shí)驗(yàn)結(jié)果的差距較大，認(rèn)為該方法不太適合于圓柱構(gòu)型的研究［11］。Chedevergne［11－13］采用了整體穩(wěn)定性方法重新討論了該問(wèn)題，解決了局部方法帶來(lái)的不一致問(wèn)題，得到的特征值虛部即時(shí)間增長(zhǎng)率皆為負(fù)值，即它們是時(shí)間穩(wěn)定的。因此，其認(rèn)為實(shí)驗(yàn)中出現(xiàn)的表面不穩(wěn)定現(xiàn)象應(yīng)該是由某些外在的噪聲激發(fā)的。

本文使用更合理的邊界條件，利用整體流動(dòng)穩(wěn)定性方法研究了圓柱構(gòu)型的徑向加質(zhì)流整體不穩(wěn)定特征，來(lái)解釋上述理論與實(shí)驗(yàn)的差異，以加深對(duì)固體火箭發(fā)動(dòng)機(jī)內(nèi)由于表面渦脫落引起的不穩(wěn)定現(xiàn)象的理解，并為分段發(fā)動(dòng)機(jī)穩(wěn)定性分析提供一種快速、有效的手段。

1 控制方程

采用不可壓縮粘性流體N－S方程，在軸對(duì)稱坐標(biāo)系下(圖 1)，無(wú)量綱速度向量 U=(Ur，Uθ，Uz)和壓強(qiáng) p滿足［11］:

其中，無(wú)量綱參考量為半徑R(長(zhǎng)度)、徑向加質(zhì)速度Vinj(速度分量)、R/Vinj(時(shí)間)、ρV2inj(壓強(qiáng))，密度 ρ假定為常數(shù)，雷諾數(shù)Re=VinjR/ν依賴于徑向加質(zhì)速度和半徑，ν為運(yùn)動(dòng)粘性系數(shù)。

把瞬態(tài)解分解為穩(wěn)態(tài)和擾動(dòng)部分的和，U=Ub+u，p=pb+p，帶入方程(1)，并消去 Ub，pb滿足的穩(wěn)態(tài)方程，可以得到擾動(dòng)方程:

流動(dòng)穩(wěn)定性分析便是在一定的邊界條件下去求解上述方程(2)，得到擾動(dòng)量u和p隨時(shí)間的變化(特征值)及其在空間的分布(特征函數(shù))。方程(2)中用到了穩(wěn)態(tài)流場(chǎng)解Ub，需要提前計(jì)算。在無(wú)粘且前端無(wú)滑移條件下，Taylor和Culick分別獨(dú)立地得到了方程(1)的穩(wěn)態(tài)精確解［10］，數(shù)值計(jì)算表明，在 Re＞100時(shí)，精確解與數(shù)值解相當(dāng)接近［7］。因此，本文采用該精確解進(jìn)行穩(wěn)定性計(jì)算，如下式所示:

2 整體穩(wěn)定性分析

局部方法假定擾動(dòng)量 q=(p，ur，uθ，uz)為如下形式:

整體穩(wěn)定性分析則把擾動(dòng)量 q=(p，ur，uθ，uz)寫成如下的形式:

式中 m為整數(shù)，代表角波數(shù);ω為復(fù)數(shù)，實(shí)部ωr代表無(wú)量綱角頻率，ωi代表時(shí)間增長(zhǎng)率。

本文只討論軸對(duì)稱模態(tài)，此時(shí)m=0，uθ=0。把式(3)代入微分方程(2)，便得到了如下的廣義特征值問(wèn)題:

其中，A與B為3×3矩陣，各分量在柱坐標(biāo)下的表達(dá)式為

求解特征值(4)需要合適的邊界條件，速度擾動(dòng)在頭部和壁面應(yīng)為0，r=0處設(shè)為軸對(duì)稱條件。文獻(xiàn)［11］在嘗試了多種出口邊界條件之后，認(rèn)為出口處取插值條件較為理想，但后來(lái)發(fā)現(xiàn)其結(jié)果與 DNS［14－15］結(jié)果相差較大。在研究后向臺(tái)階的瞬態(tài)增長(zhǎng)時(shí)，Blackburn［16］認(rèn)為出流邊界處取法向應(yīng)力為0更為合理，本文出口處便采用了該種邊值條件。計(jì)算所用邊界條件總結(jié)如下:

3 結(jié)果分析

利用譜配置方法［17－18］和邊值條件(5)求解特征值問(wèn)題(4)，徑向×軸向網(wǎng)格量為20×50。數(shù)值計(jì)算了Re=2 100時(shí)的特征模態(tài)，如圖2所示。可看到長(zhǎng)徑比L/R分別為10和8時(shí)，兩者特征值的實(shí)部，即無(wú)量綱角頻率ωr相差不大，但前者的虛部，即能量時(shí)間增長(zhǎng)率ωi明顯大于后者，并且有些越過(guò)了臨界值0，由穩(wěn)定模態(tài)變成了不穩(wěn)定模態(tài)。這證明了文獻(xiàn)［15］中的猜測(cè)，即當(dāng)長(zhǎng)徑比大到一定程度后，ωi將為正值，時(shí)間特征由衰減變?yōu)樵鲩L(zhǎng)。已有的冷流實(shí)驗(yàn)［19－21］中長(zhǎng)徑比都大于10，且都出現(xiàn)了表面不穩(wěn)定現(xiàn)象，這與上述的理論結(jié)果一致，證明了外部噪聲不是產(chǎn)生表面不穩(wěn)定的必要條件。

不同長(zhǎng)徑比及特征值的擾動(dòng)速度的分布見圖3。從圖3(a)看出，無(wú)論軸向擾動(dòng)速度(上圖)，還是徑向擾動(dòng)速度(下圖)，其沿軸向都是逐漸增大的，且向內(nèi)部擴(kuò)張。同一徑向位置，接近壁面的地方擾動(dòng)達(dá)到最大，這與文獻(xiàn)［19］中實(shí)驗(yàn)結(jié)果一致。徑向擾動(dòng)速度在z=2的地方就已存在，其幅值比軸向擾動(dòng)速度幅值小了一個(gè)數(shù)量級(jí)。圖3(b)為L(zhǎng)/R=8時(shí)擾動(dòng)速度的分布，該特征值與圖3(a)中的特征值有著相近的實(shí)部，但其虛部要小一些。兩圖中的擾動(dòng)速度分布很相似，幅值也相差不大，只不過(guò)L/R=8時(shí)擾動(dòng)存在的位置更靠前。與圖3(a)相比，圖3(c)不穩(wěn)定波的波長(zhǎng)較短，徑向擾動(dòng)速度出現(xiàn)位置明顯靠后。

4 實(shí)驗(yàn)驗(yàn)證

為了與實(shí)驗(yàn)結(jié)果比較，需把理論求解得到的一系列特征值ωM(其實(shí)部為角頻率)，量綱化為

文獻(xiàn)［19］利用徑向加質(zhì)軸對(duì)稱冷流實(shí)驗(yàn)器對(duì)固體發(fā)動(dòng)機(jī)燃燒室內(nèi)的氣動(dòng)聲學(xué)做了深入的研究。實(shí)驗(yàn)器半徑 R=30 mm，空氣徑向加質(zhì)速度 Vinj可在 0.8 ～2.0 m/s之間變化。在加質(zhì)邊界附近測(cè)得的不穩(wěn)定頻率隨徑向加質(zhì)速度的變化如圖4所示，圖4(a)為徑向加質(zhì)速度遞增時(shí)的擾動(dòng)速度頻率，圖4(b)為徑向加質(zhì)速度遞減時(shí)的擾動(dòng)速度頻率。從圖4看到，理論和實(shí)驗(yàn)吻合得很好。實(shí)驗(yàn)中的不穩(wěn)定頻率有數(shù)條軌跡，但它們都在理論解的范圍內(nèi)。在徑向加質(zhì)速度大于1.7 m/s時(shí)，實(shí)驗(yàn)結(jié)果開始偏離理論值，應(yīng)歸因于實(shí)驗(yàn)中使用的供氣系統(tǒng)，壓力過(guò)高時(shí)，其與徑向加質(zhì)速度的線性關(guān)系可能已不再適合。

Dunlap［20］設(shè)計(jì)的冷流實(shí)驗(yàn)器半徑R=50 mm，氮?dú)鈴较蚣淤|(zhì)馬赫數(shù)見表1，實(shí)驗(yàn)中氣體溫度在257～286 K之間，本文計(jì)算中選取了溫度的兩邊界值，取上界時(shí)，偏差在2%以內(nèi)，取下界時(shí)，偏差在6%內(nèi)。

表1 理論結(jié)果和文獻(xiàn)［20］中實(shí)驗(yàn)不穩(wěn)定頻率的比較Table 1 Comparison between theoretical and experimental results［20］ on frequencies of instability

VKI流體力學(xué)研究中心Anthoine［21］設(shè)計(jì)了Ariane 5 MPS的1/30縮比冷流實(shí)驗(yàn)器來(lái)研究表面渦脫落引起的壓力振蕩。實(shí)驗(yàn)器 L=0.63 m，D=0.076 m，噴管入口處馬赫數(shù)Ma0=0.09，按下式計(jì)算得到徑向加質(zhì)速度［21］:

實(shí)驗(yàn)中在L/R=10處測(cè)得的不穩(wěn)定頻率為270 Hz和350 Hz，本文理論計(jì)算結(jié)果為273 Hz和337 Hz，偏差在4%以內(nèi)。

5 結(jié)論

(1)圓柱形徑向加質(zhì)流的整體不穩(wěn)定分析發(fā)現(xiàn)特征譜是離散的，且隨著長(zhǎng)徑比的增加，特征值的虛部明顯變大，而其實(shí)部變化不大，在L/R=10時(shí)虛部從時(shí)間穩(wěn)定變?yōu)椴环€(wěn)定。

(2)特征函數(shù)表明，軸向和徑向擾動(dòng)速度沿軸向逐漸增大，且向內(nèi)部擴(kuò)張，同一徑向位置，接近壁面的地方擾動(dòng)達(dá)到最大。

(3)理論結(jié)果準(zhǔn)確地估計(jì)了實(shí)驗(yàn)中出現(xiàn)的不穩(wěn)定頻率，證明流動(dòng)穩(wěn)定性方法可以解釋固體火箭發(fā)動(dòng)機(jī)內(nèi)的表面不穩(wěn)定現(xiàn)象，也可以對(duì)發(fā)動(dòng)機(jī)內(nèi)的表面渦脫落現(xiàn)象和流動(dòng)穩(wěn)定性進(jìn)行預(yù)估。

［1］Dotson K W，Koshigoe S，Pace K K.Vortex shedding in a large solid rocket motor without inhibitors at the segment interfaces［J］.Journal of Propulsion and Power，1997，13(2).

［2］Zhang Q，Wei Z J，Su W X，et al，Theoretical modeling and numerical study for thrust oscillation characteristics in solid rocket motors［J］.Journal of Propulsion and Power，2012，28(2).

［3］Lupoglazoff N，Vuillot F.Parietal vortex shedding as a cause of instability for long solid propellant motors:numerical simulations and comparisons with firing tests［R］.AIAA 96－0761.

［4］Avalon G，Casalis G，Griffond J.Flow instabilities and acoustic resonance of channels with wall injection［R］.AIAA 98－3218.

［5］Avalon G，Ugurtas B，Grisch F，et al，Numerical computations and visualization tests of the flow inside a cold gas simulation with characterization of a parietal vortex shedding［R］.AIAA 2000－3387.

［6］Anthoine J，Lema M R.Passive control of pressure oscillations in solid rocket motors:cold－flow experiments［J］.Journal of Propulsion and Power，2009，25(3).

［7］Casalis G，Avalon G，Pineau J P.Spatial instability of planar channel flow with fluid injection through porous walls［J］.Physics of Fluids，1998，10(10).

［8］Griffond J，Casalis G.On the dependence on the formulation of some nonparallel stability approaches applied to the taylor flow［J］.Physics of Fluids，2000，12(2).

［9］Griffond J，Casalis G.On the nonparallel stability of the injection induced two－dimensional taylor flow［J］.Physics of Fluids，2001，13(6).

［10］Griffond J，Casalis G，Pineau J P.Spatial instability of flow in a semi－infinite cylinder with fluid injection through its porous walls［J］.Eur.J.Mech.B－Fluids，2000，19.

［11］Chedevergne F，Casalis G，F(xiàn)eraille T.Biglobal linear stability analysis of the flow induced by wall injection［J］.Physics of Fluids，2006，18.

［12］Chedevergne F，Casalis G.Thrust oscillations in reduced scale solid rocket motors，part Ⅱ:a new theoretical approach［R］.AIAA 2005－4000.

［13］Chedevergne F，Casalis G.Detailed analysis of the thrust oscillations in reduced scale solid rocket motors［R］.AIAA 2006－4424.

［14］Casalis G，Boyer G，Radenac E.Some recent advances in the instabilities occurring in long solid rocket motors［R］.AIAA 2011－5642.

［15］Chedevergne F，Casalis G，Majdalani J.DNS investigation of the taylor－culick flow stability［R］.AIAA 2007－5796.

［16］Blackburn H M，Barkley D，Sherwin S J.Convective instability and transient growth in flow over a backward－facing step［J］.Journal of Fluid Mechanics，2008:603.

［17］Weideman J A C，Reddy S C.A matlab differentiation matrix suite［J］.ACM Trans.Math.Softw.，2000，26(4).

［18］Trefethen L N.Spectral methods in matlab［J］.Society for Industrial and Applied Mathematics，Philadelphia，2000.

［19］Cerqueira S，Avalon G，F(xiàn)eyel F.An experimental investigation of fluid－structure interaction inside solid propellant rocket motors［R］.AIAA 2009－5427.

［20］Dunlap R，Blackner A，Waugh R，et al.Internal flow field studies in a simulated cylindrical port rocket chamber［J］.Journal of Propulsion and Power，1990，6(6).

［21］Anthoine J.Experimental and numerical study of aeroacoustic phenomena in large solid propellant motors［D］.PhD thesis，Universite libre，Bruxelles，October，2000.