李 陽，劉佩進(jìn)，金秉寧

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

0 引言

國內(nèi)外大型分段固體火箭發(fā)動機(jī)經(jīng)常出現(xiàn)低頻率的壓強(qiáng)振蕩現(xiàn)象[1-4]。研究發(fā)現(xiàn)，大部分的壓強(qiáng)振蕩都是由于周期性渦脫落所引起的[5]。由于對渦脫落自身形成和發(fā)展規(guī)律的認(rèn)識不足，且對渦脫落與發(fā)動機(jī)聲場之間相互作用的研究不足，目前工程上還沒有實用方法預(yù)示渦脫落和壓強(qiáng)振蕩。渦脫落與聲場的耦合可利用高精度CFD模擬[6]，但由于發(fā)動機(jī)幾何尺寸大，工作時間長，對每一時刻的工作狀態(tài)進(jìn)行CFD計算并不現(xiàn)實。從流動穩(wěn)定性理論出發(fā)，對渦脫落與聲場耦合的機(jī)理進(jìn)行研究，獲得流動穩(wěn)定特征以及演化規(guī)律，是當(dāng)前學(xué)術(shù)界的研究熱點[7]。

Flandro利用流動穩(wěn)定性理論和能量平衡方法，建立了固體火箭發(fā)動機(jī)內(nèi)的渦聲耦合模型，在動量方程中引入局部體積力來描述渦的反饋效果，計算了由渦脫落產(chǎn)生的能量增長指數(shù)[8]。而Vuillot將Flandro的方法用到實驗發(fā)動機(jī)C1x的設(shè)計過程中時，發(fā)現(xiàn)局部體積力的引入具有較大的任意性，計算結(jié)果敏感地依賴渦從下游撞擊點到剪切層原點的距離(此距離并非確定值，且有一定估計偏差)。因此，結(jié)果的可靠性不足[9]。Griffond采用局部非平行方法，計算了特定工況下Taylor-Culick流的穩(wěn)定性，并與Dunlap[10]的冷流實驗做了比較，得到了兩個特定工況下的不穩(wěn)定模態(tài)。但中性曲線的分析僅描述了特定工況下方程解的分布，并不適合描述物理過程[11]。Akiki研究可壓縮Taylor平面流的穩(wěn)定性時，數(shù)值方法與解析方法分別計算得到的幅值出現(xiàn)偏差，他認(rèn)為是解析方法把有旋量和無旋量分開的求解過程造成的[12]。但Chedevergne發(fā)現(xiàn)解析方法得到的流動穩(wěn)定性模態(tài)可疊加，并能準(zhǔn)確地重現(xiàn)DNS結(jié)果。因此，Akiki的解釋并不令人信服[13]。楊尚榮[6]應(yīng)用局部非平行理論，分析了由主變量公式和流函數(shù)公式導(dǎo)出的擾動方程的差異，發(fā)現(xiàn)在不發(fā)生共振的情況下，理論方法可預(yù)估發(fā)動機(jī)內(nèi)的流動不穩(wěn)定頻率，但并未討論理論方法在火箭發(fā)動機(jī)流動穩(wěn)定性預(yù)估計算中的實際應(yīng)用。

本文利用譜配置方法，對不同結(jié)構(gòu)及工況下的Taylor-Culick流進(jìn)行分析，研究結(jié)構(gòu)參數(shù)、流動擾動參數(shù)變化對穩(wěn)定性的影響。本文討論了對固體火箭發(fā)動機(jī)進(jìn)行穩(wěn)定性預(yù)估時，側(cè)向注入雷諾數(shù)、長徑比、頻率、環(huán)向波數(shù)與空間穩(wěn)定性的相互關(guān)系，得到了不穩(wěn)定模態(tài)的變化規(guī)律。

1 控制方程與求解方法

本文研究的徑向加質(zhì)Taylor-Culick流模型，是固體火箭發(fā)動機(jī)燃燒室前部內(nèi)流場的一種簡化模型，模型中不包含噴管附近的高馬赫數(shù)流動部分，幾何構(gòu)型與坐標(biāo)如圖1所示。固體火箭發(fā)動機(jī)前部的燃燒室，流動馬赫數(shù)小于0.3，適合采用不可壓流體控制方程。

利用歸一化參考量(半徑R、徑向加質(zhì)速度Vinj、密度ρ*、運(yùn)動粘性系數(shù)υ)對不可壓粘性流體N-S方程進(jìn)行無量綱化。將歸一化N-S方程中的瞬時變量分解為平均量和擾動量之和，得到關(guān)于擾動量的方程：

▽·u′=0

(1)

平均量作為基本流可以求出解析解：

(2)

利用分離變量法，假設(shè)擾動量為簡正模態(tài)形式：

=(ur,uθ,uz,p)(r)ei(mθ+αz-ωt)

(3)

其中，ω為實數(shù)，代表無量綱擾動頻率；α=αr-iαi為復(fù)數(shù)，實部αr為軸向擾動波數(shù)，負(fù)虛部-αi為軸向擾動增長率。m為正整數(shù)，代表環(huán)向波數(shù)。進(jìn)行空間穩(wěn)定性分析時，給定頻率ω和環(huán)向波數(shù)m，利用譜配置方法可求得任意軸向位置z處的增長率-αi和波數(shù)αr。

將簡正形式的擾動量代入擾動方程，可化為如下關(guān)于α的多項式的特征值問題：

(4)

Taylor-Culick流的邊界條件為壁面上和對稱軸上的速度邊界條件：頭部壁面速度為零、側(cè)壁上注入速度為常數(shù)、側(cè)壁上軸向速度無滑移以及對稱軸上的徑向速度對稱條件。利用擾動量表示物理邊界條件，得到擾動方程的齊次邊界條件。利用擾動方程，還可得到壓力邊界條件。

(ur，uθ，uz)(±1)=0,Dp(±1)=0

(5)

利用譜配置方法，選擇切比雪夫多項式Tk(x)=cos(karccos(x))進(jìn)行空間離散；配置點選擇Gauss-Lobatto積分點：

(6)

權(quán)函數(shù)α(x)=1。通過MatLab編程計算得到相應(yīng)結(jié)果[14-15]。

2 計算結(jié)果與討論

2.1 計算結(jié)果驗證

作為驗證，計算了長徑比z=10，側(cè)向注入雷諾數(shù)Re=4500，頻率ω=80，環(huán)向波數(shù)m=0下的軸向波數(shù)、軸向增長率與不穩(wěn)定模態(tài)，得到的特征譜如圖2所示。圖2中，兩個不穩(wěn)定點即為兩個不穩(wěn)定模態(tài)，對應(yīng)的橫坐標(biāo)為軸向波數(shù)、縱坐標(biāo)為軸向增長率。Griffond通過設(shè)定初值，并利用泰勒展開的牛頓-拉夫森迭代法計算了相應(yīng)結(jié)果，其定性分析與dunlop實驗結(jié)果[10]沒有矛盾。本文結(jié)果與Griffond的結(jié)果[11]比較，如表1所示。

本文所得結(jié)果與文獻(xiàn)結(jié)果基本吻合，可驗證本文所用方法可行。

利用譜方法計算的結(jié)果，通過側(cè)向注入雷諾數(shù)Re、長徑比z、頻率ω及環(huán)向波數(shù)m的變化與軸向波數(shù)的關(guān)系，來討論流動不穩(wěn)定現(xiàn)象。計算結(jié)果基于歸一化擾動方程，量綱為1，對于不同的物理構(gòu)型及物理量，可換算成真實數(shù)值。定性分析利用無量綱參變量代替物理數(shù)值。

模態(tài)Griffond的結(jié)果[11]本文結(jié)果誤差/‰16．09529456566．09529597-1．0787998140-1．078801377<123．32642853663．326428826-0．1095525589-0．1095526851<1

2.2 側(cè)向注入雷諾數(shù)的影響

側(cè)向注入雷諾數(shù)定義為Re=RVinj/υ，其中R為圓管半徑，不可壓縮流體運(yùn)動粘度系數(shù)為常數(shù)。不可壓縮流體的密度為常數(shù)，動力粘度確定后，運(yùn)動粘度系數(shù)也隨之確定。由定義可知，當(dāng)圓管半徑一定時，加質(zhì)速度越大，側(cè)向注入雷諾數(shù)越大。因此，側(cè)向注入雷諾數(shù)的變化反映了徑向加質(zhì)速度變化對流動不穩(wěn)定的影響。

在不同側(cè)向注入雷諾數(shù)下，流動會產(chǎn)生不同的狀態(tài)。低側(cè)向注入雷諾數(shù)不易出現(xiàn)不穩(wěn)定情況，隨著側(cè)向注入雷諾數(shù)增大產(chǎn)生的不穩(wěn)定模態(tài)增多(對于圖3中情況，對應(yīng)不穩(wěn)定模態(tài)增加的側(cè)向注入雷諾數(shù)臨界值約為100、750)。不同的不穩(wěn)定點隨側(cè)向注入雷諾數(shù)變化的趨勢并不相同。不穩(wěn)定點1是最早出現(xiàn)的對應(yīng)于小側(cè)向注入雷諾數(shù)的波動，隨著側(cè)向注入雷諾數(shù)增大，其軸向波數(shù)及波動增長率都是不斷增大的，側(cè)向注入雷諾數(shù)在1500以內(nèi)時，頻率和增長率增幅較大，而隨后增幅趨緩，并逐漸趨近于有限值。這種現(xiàn)象符合物理規(guī)律，不穩(wěn)定波動并不是無限增大的。不穩(wěn)定點2點對應(yīng)于中側(cè)向注入雷諾數(shù)波動，此點在側(cè)向注入雷諾數(shù)達(dá)到一定程度時才會出現(xiàn)，軸向增長率逐漸增大，并趨近于有限值，而軸向波數(shù)逐漸減小逐漸由高頻變?yōu)榈皖l。所以，隨著側(cè)向注入雷諾數(shù)的增加，高低頻不穩(wěn)定模態(tài)會發(fā)生轉(zhuǎn)換。這意味著某些工作狀態(tài)下(如側(cè)向加質(zhì)雷諾數(shù)的變化過程中對應(yīng)的波動頻率正好跨過某一階聲腔頻率)，更易引發(fā)流動不穩(wěn)定。

固體火箭發(fā)動機(jī)內(nèi)是大雷諾數(shù)流動，一般在103量級，上文討論范圍選取0～7000，僅為全面反映TC流模型側(cè)向加質(zhì)雷諾數(shù)的影響，固體火箭發(fā)動機(jī)的討論限定在圖3中大側(cè)向加質(zhì)雷諾數(shù)的部分。

2.3 長徑比的影響

對于本文所述模型，長徑比是固體火箭發(fā)動機(jī)內(nèi)燃?xì)饬鲃涌涨坏奈锢肀碚?。固體火箭發(fā)動機(jī)長徑比超過10稱為大長徑比發(fā)動機(jī)，但一般不會超過15。為直觀地反映長徑比變化引起不穩(wěn)定性變化的趨勢、描述流動不穩(wěn)定的物理現(xiàn)象，本文選取長徑比范圍0～15。圖4為長徑比的變化對特征值的影響，反映了物理構(gòu)型與流動穩(wěn)定性的關(guān)系。

由圖4可看出，當(dāng)長徑比增加，特征值有兩段變化趨勢。當(dāng)長徑比小于3.5時，隨著長徑比的增加，軸向波數(shù)急劇下降，而軸向增長率急劇增長；當(dāng)長徑比大于3.5時，軸向波數(shù)及增長率都隨著長徑比的增加而緩慢降低。這與物理現(xiàn)象是內(nèi)在一致的，當(dāng)長徑比小于3.5時，均勻流場流線曲率大，加質(zhì)注入的速度激勵及速度方向的改變對不穩(wěn)定的增長作用很大；當(dāng)長徑比大于3.5時，均勻流場流線近似平行于軸向，表現(xiàn)為均勻流場對軸向不穩(wěn)定的阻尼增強(qiáng)，軸向增長率以指數(shù)形式變化，并趨于穩(wěn)定值。

在發(fā)動機(jī)設(shè)計中，除了要避開聲共振的固有頻率外，還要注意構(gòu)型對流動不穩(wěn)定的影響，特別要注意長徑比為3.5的構(gòu)型。此時，擾動的放大效果最強(qiáng)，較小的擾動便可產(chǎn)生較大幅度的振動，不利于發(fā)動機(jī)正常工作。長徑比對不穩(wěn)定性的影響體現(xiàn)在對擾動的放大效果上，并不意味著長徑比為3.5時，流動不穩(wěn)定性最易發(fā)生。

2.4 頻率的影響

ω表示時間方向上的波動頻率，其量綱化公式為

頻率在一定范圍內(nèi)對于空間穩(wěn)定性產(chǎn)生影響，計算得到頻率ω影響范圍為無量綱區(qū)間[50, 190]，見圖5。

不穩(wěn)定模態(tài)數(shù)量(即不穩(wěn)定點數(shù))會隨著頻率的增長而變化，見圖5。頻率達(dá)到66～67時，不穩(wěn)定點由出現(xiàn)一個變?yōu)閮蓚€(1→2)；頻率在之后達(dá)到116～117、156～157、182～183以及188～189時，不穩(wěn)定點數(shù)分別出現(xiàn)由(2→3)、(3→2)、(2→1)以及(1→0)的變化。但不管出現(xiàn)多少不穩(wěn)定點，其變化規(guī)律是相同的，軸向不穩(wěn)定波數(shù)不斷增長，新的低頻軸向不穩(wěn)定模態(tài)不斷出現(xiàn)；而當(dāng)頻率達(dá)到一定的值后，流動不穩(wěn)定模態(tài)消失，流動變?yōu)榭臻g穩(wěn)定；不穩(wěn)定的軸向增長率都是先增加后達(dá)到極值，然后快速減小?？梢?，時間波動與空間波動并不是同步的，隨著頻率的增強(qiáng)，空間波動先增強(qiáng)、再減小，然后變?yōu)榉€(wěn)定流動。在工程中發(fā)展對流動不穩(wěn)定的抑制方法時，可利用此特性。

2.5 環(huán)向波數(shù)的影響

環(huán)向波數(shù)m在擾動方程的參量中表現(xiàn)為eimθ，環(huán)向波數(shù)增大，則環(huán)向波動頻率增大。對于小擾動，理論上環(huán)向與軸向波數(shù)應(yīng)在同一量級，數(shù)值計算得到同樣的結(jié)果。

圖6為環(huán)向波數(shù)變化對軸向波數(shù)與軸向增長率的影響。隨著環(huán)向波數(shù)的變化，不穩(wěn)定模態(tài)的個數(shù)也發(fā)生變化。從特征譜(圖7)中可了解到，隨著環(huán)向波數(shù)的增加，不穩(wěn)定點(圖中右下角三點沿軸向依次為Ⅰ、Ⅱ、Ⅲ點)逐漸由不穩(wěn)定模態(tài)變?yōu)榉€(wěn)定模態(tài)，不同的不穩(wěn)定點其變化規(guī)律是相同的，其中不穩(wěn)定點Ⅱ、Ⅲ最具有代表性。隨著環(huán)向波數(shù)的增加，軸向波數(shù)在一定范圍內(nèi)增加(見圖6(a))，這說明各坐標(biāo)上的空間不穩(wěn)定是協(xié)同的，表現(xiàn)為頻率的共同增長；但隨著環(huán)向波數(shù)繼續(xù)增大，軸向波數(shù)卻隨之迅速減小(見圖6(b))，表現(xiàn)為頻率的遷移。軸向增長率卻是逐漸減小，最后變?yōu)榉€(wěn)定點。綜合來看，流動不穩(wěn)定效果在軸向、徑向以及切向上表現(xiàn)為頻率共同增長，波數(shù)此消彼長。

3 結(jié)論

(1)側(cè)向注入雷諾數(shù)決定基頻不穩(wěn)定模態(tài)，隨著側(cè)向注入雷諾數(shù)增大，振蕩模態(tài)可能增加，振蕩頻率可能發(fā)生由高頻變?yōu)榈皖l的轉(zhuǎn)換；側(cè)向注入雷諾數(shù)影響振蕩幅值的放大比率，側(cè)向注入雷諾數(shù)越大，軸向不穩(wěn)定振幅放大率越大。

(2)構(gòu)型長徑比為3.5的固體火箭發(fā)動機(jī)在發(fā)生壓強(qiáng)不穩(wěn)定振蕩時的軸向增長率最大，較小的擾動便有可能產(chǎn)生較大振動。

(3)頻率范圍影響空間不穩(wěn)定出現(xiàn)的可能性，可利用此特性發(fā)展對固體火箭發(fā)動機(jī)發(fā)生流動不穩(wěn)定主動抑制的方法；環(huán)向波數(shù)影響不穩(wěn)定模態(tài)數(shù)目，環(huán)向波數(shù)與軸向波數(shù)具有相關(guān)性，幅度總和趨于有限值。

[1] 劉佩進(jìn),魏少娟,楊尚榮.固體火箭發(fā)動機(jī)非線性燃燒穩(wěn)定分析方法[J].固體火箭技術(shù),2013,36(3):302-309.

LIU Peijin,WEI Shaojuan,YANG Shangrong.Nonlinear combustion instability analysis method of solid rocket motor[J].Journal of Solid Rocket Technology,2013,36(3):302-309.

[2] 吳亞可,何國強(qiáng),劉佩進(jìn),等．潛入式噴管背區(qū)空腔對壓強(qiáng)振蕩及旋渦運(yùn)動影響的實驗研究[J].固體火箭技術(shù),2010,33(6):621-625.

WU Yake,HE Guoqiang,LIU Peijin,et al.Experimental study of influence of cavity formed by submerged nozzle on pressure oscillation and vortices motion[J].Journal of Solid Rocket Technology,2010,33(6):621-625.

[3] 王建儒,何國強(qiáng),許團(tuán)委,等.燃燒室對接狹縫設(shè)計參數(shù)對壓強(qiáng)振蕩的影響研究[J].固體火箭技術(shù),2012,35(4):474-478+484.

WANG Jianru,HE Guoqiang,XU Tuanwei,et al.Effect of internal joint design parameters on pressure oscillation in combustion chamber[J].Journal of Solid Rocket Technology,2012,35(4):474-478+484.

[4] 張翔宇,何國強(qiáng),劉佩進(jìn).軸對稱后向臺階不穩(wěn)定流動及壓強(qiáng)振蕩數(shù)值研究[J].固體火箭技術(shù),2012,35(2):171-176+182.

ZhANG Xiangyu,HE Guoqiang,LIU Peijin.Numerical study on unsteady backward-facing step flow and pressure oscillations[J].Journal of Solid Rocket Technology,2012,35(2):171-176+182.

[5] 劉佩進(jìn),楊尚榮.分段固體火箭發(fā)動機(jī)中聲渦耦合現(xiàn)象的實驗研究現(xiàn)狀[J].固體火箭技術(shù),2012,35(6):726-731.

LIU Peijin,YANG Shangrong.An overview on experimental study of vortex-acoustic coupling in segmented solid rocket motors[J].Journal of Solid Rocket Technology,2012,35(6):726-731.

[6] 李強(qiáng),何國強(qiáng),劉佩進(jìn),等.模型固體火箭發(fā)動機(jī)內(nèi)渦脫落過程的大渦模擬[J].推進(jìn)技術(shù),2005,26(6):481-484.

LI Qiang,HE Guoqiang,LIU Peijin,et al.Large eddy simulation of vortex shedding for model SRM[J].Journal of Propulsion Technology,2005,26(6):481-484.

[7] 楊尚榮,劉佩進(jìn),魏祥庚,等.固體火箭發(fā)動機(jī)徑向加質(zhì)流局部穩(wěn)定性分析方法研究[J].推進(jìn)技術(shù),2013,34(7):925-931.

YANG Shangrong,LIU Peijin,WEI Xianggeng,et al.A study of analysis method for hydrodynamic stabilities of solid rocket motors with side-induced flow[J].Journal of Propulsion Technology,2013,34(7):925-931.

[8] Flandro G A.Vortex driving mechanisms in oscillatory rocket flows[J].AIAA Journal,1986,2(3):206-214.

[9] Vuillot F.Vortex-shedding phenomena in solid rocket motors[J].Journal of Propulsion and Power,1995,11(4):626-639.

[10] Dunlap R,Blackner A M,Waugh R C,et al.Internal flow field studies in a simulated cylindrical port rocket chamber[J].Journal of Propulsion and Power,1990,6(6):690-704.

[11] Griffond J,Casalis G and Pineau J P.Spatial instability of flow in a semi infinite cylinder with fluid injection through its porous walls[J].European Journal of Mechanics.B/Fluids,2000,19(1):69-87.

[12] Akiki M,Batterson J W and Majdalani J.Biglobal stability of compressible flow fields.part 2:application to solid rocket motors[R].AIAA 2013-3866.

[13] Chedevergne F,Casalis G and Majdalani J.Direct numerical simulation and biglobal stability investigations of the gaseous motion in solid rocket motors[J].Journal of Fluid Mechanics,2012,706(3):190-218.

[14] 何松年.譜方法及其應(yīng)用[J].中國民航學(xué)院學(xué)報(綜合版),2000,18(6):30-33.

HE Songnian.Spectral method and its application[J].Journal of Civil Aviation University of China,2000,18(6):30-33.

[15] 孫德軍,尹協(xié)遠(yuǎn),童秉綱.關(guān)于Orr-Sommerfeld方程的Chebyshev譜方法的討論[J].上海力學(xué),1998(1):1-8.SUN Dejun,YIN Xieyuan,TONG Binggang.On Chebyshev spectral methods for the Orr-Sommerfeld equation[J].Chinese Quarterly of Mechanics,1998(1):1-8.