某固體火箭發(fā)動機過載環(huán)境下內(nèi)流場研究①

張衛(wèi)平，楊　丹，熊永亮

(1.中國航天科工集團公司六院四十一所，呼和浩特　010010；2.華中科技大學(xué) 船舶與海洋工程學(xué)院，武漢　430074；3.華中科技大學(xué) 力學(xué)系，武漢　430074)

?

某固體火箭發(fā)動機過載環(huán)境下內(nèi)流場研究①

張衛(wèi)平1，楊丹2，熊永亮3

(1.中國航天科工集團公司六院四十一所，呼和浩特010010；2.華中科技大學(xué) 船舶與海洋工程學(xué)院，武漢430074；3.華中科技大學(xué) 力學(xué)系，武漢430074)

文中對某固體火箭發(fā)動機工作末期出現(xiàn)的壓力振蕩現(xiàn)象進行數(shù)值研究，采用大渦模擬對裝藥燃燒后期時刻進行數(shù)值模擬，得到了燃燒室內(nèi)裝藥轉(zhuǎn)角產(chǎn)生的旋渦脫落現(xiàn)象。通過對燃燒室內(nèi)不同位置的流場參數(shù)進行監(jiān)測，得到了由于旋渦脫落產(chǎn)生的壓力振蕩頻率。同時，對發(fā)動機軸向過載工作條件進行了模擬，得到了過載條件對發(fā)動機內(nèi)壓力振蕩的影響。結(jié)果表明，過載瞬間對發(fā)動機內(nèi)部產(chǎn)生一個大幅壓力脈沖，且發(fā)動機內(nèi)不同位置處的影響是不同的，過載使得發(fā)動機尾部旋渦脫落的主頻發(fā)生小幅度增大。最后對發(fā)動機頭部裝藥結(jié)合進行優(yōu)化，改進后的內(nèi)流場旋渦強度得到大幅減小。

固體火箭發(fā)動機；不穩(wěn)定燃燒；壓力振蕩；過載

0　引言

固體火箭發(fā)動機不穩(wěn)定燃燒，又稱振蕩燃燒，是研制固體火箭發(fā)動機經(jīng)常會遇到的一個具有挑戰(zhàn)性的技術(shù)難題。其典型的特征是發(fā)動機內(nèi)部工作壓力-時間曲線不規(guī)則，進而使得發(fā)動機的推力振蕩與設(shè)計要求不符，結(jié)構(gòu)發(fā)生振蕩等多種危害[1]。經(jīng)過半個世紀的深入研究，人們對固體火箭發(fā)動機的不穩(wěn)定燃燒取得了一些理論與經(jīng)驗性的共識。然而，隨著發(fā)動機的多樣化需求以及工作環(huán)境的惡化，目前又有越來越多的火箭發(fā)動機遇到不穩(wěn)定燃燒的問題[2]。

通過對法國Ariane 5號發(fā)生不穩(wěn)定燃燒的P230固體火箭助推器的大量研究，Dotson[3]等提出旋渦脫落是引起固體火箭發(fā)動機壓強振蕩的主要原因，旋渦脫落的頻率與發(fā)動機的固有聲振頻率耦合，激發(fā)了發(fā)動機內(nèi)流動的不穩(wěn)定性。最初有關(guān)聲渦耦合機理的描述是由Flandro和Jacobs[4]一起提出的，他們在1973年首次指出旋渦脫落有可能激發(fā)大型分段式固體火箭發(fā)動機燃燒室的聲模態(tài)產(chǎn)生共振，從而引起壓強振蕩。流體中的旋渦既是一種擾動，也是釋放聲能的流動聲源，流場中大量存在的旋渦與發(fā)動機內(nèi)部聲振之間的復(fù)雜非線性耦合作用，也是近年來發(fā)動機不穩(wěn)定燃燒研究的熱點問題。Space Shuttle RSRM,Titan IV SRMU 與Ariane 5號助推器P230發(fā)動機等一些大型分段式固體火箭發(fā)動機，由于燃燒室凸出的障礙物或者裝藥的突擴結(jié)構(gòu)，在其下游誘導(dǎo)出強剪切流作用的旋渦脫落。發(fā)動機內(nèi)的旋渦脫落主要可歸結(jié)為轉(zhuǎn)角旋渦脫落[5]、障礙物旋渦脫落[5]以及表面旋渦脫落[6-7]三類。Dotson[3]針對發(fā)生推力振蕩的Tian IV SRMU發(fā)動機，通過試驗測量建立了旋渦從形成、脫落、碰撞、聲反饋的4步模型。利用數(shù)值計算[8-11]，人們也捕捉到了發(fā)動機內(nèi)的旋渦脫落，結(jié)合流場與聲場的能量交換機理來研究渦/聲耦合現(xiàn)象，并取得了一些結(jié)論，但目前對渦/聲耦合機理的理解仍存在分歧。

特別是戰(zhàn)術(shù)火箭通常燃燒室更短，不應(yīng)輕易地認為發(fā)生了渦/聲耦合。即使很長的戰(zhàn)術(shù)火箭，基頻也在150 Hz以上。渦的作用和聲振的作用不容易使發(fā)動機達到極大振幅，渦聲耦合是極難出現(xiàn)的。引起發(fā)動機不穩(wěn)定燃燒的因素是多方面的，旋渦與外界干擾或者過載產(chǎn)生的發(fā)動機簡單殼體振動響應(yīng)及流動響應(yīng)，也是有可能激發(fā)發(fā)動機不穩(wěn)定燃燒的。因此，對發(fā)動機內(nèi)旋渦流動進行研究，對深入理解發(fā)動機壓力振蕩具有重要意義，本文的主要研究內(nèi)容為發(fā)動機內(nèi)旋渦脫落現(xiàn)象。

另外，發(fā)動機在工作過程中，由于加速和機動性能改變等因素，不僅會影響總體內(nèi)彈道性能以及產(chǎn)生局部燒蝕[12]，同時也會對內(nèi)部壓力振蕩產(chǎn)生影響。本文針對發(fā)動機內(nèi)旋渦流動現(xiàn)象，對某固體火箭發(fā)動機內(nèi)流場進行數(shù)值模擬，并對該發(fā)動機加載前后瞬間進行對比，研究發(fā)動機發(fā)生縱向過載條件下的流動穩(wěn)定性，目的是為了得到發(fā)動機過載條件對其內(nèi)流場不穩(wěn)定流動的影響。最后，通過改變裝藥結(jié)構(gòu)，對比發(fā)動機尾部的旋渦強度。

1　流場數(shù)學(xué)模型與計算方法

通過大渦模擬方法求解可壓縮Navier-Stokes方程組，對發(fā)動機內(nèi)流場進行數(shù)值模擬，流場控制方程可通過對原始N-S方程進行空間濾波后得到。對于可壓縮流動，引入密度加權(quán)平均(Favre)，可表示為

(1)

(2)

式中D為流場區(qū)域；G為濾波函數(shù)。

濾波函數(shù)定義如下：

(3)

式中V為計算單元的體積。

對N-S方程進行濾波后，得到控制方程如下：

(4)

(5)

(6)

式中“—”表示雷諾(Reynolds)平均；“～”表示密度加權(quán)平均。

分子粘性應(yīng)力張量定義為

(7)

式中μ為動力粘度。

亞格子應(yīng)力相表示為

(8)

采用壁掛式本地渦粘性(Wall-adapting local eddy-iscosity,WALE)亞格子模型，在WALE模型中，渦粘度由式(9)計算：

(9)

Ls=min(kd,C?V1/3)

本文選取發(fā)動機工作末期發(fā)生壓力振蕩的某一時刻進行數(shù)值計算，計算區(qū)域及邊界條件設(shè)置如圖1所示。圖1中，同時給出了改進后的裝藥結(jié)構(gòu)。發(fā)動機噴管網(wǎng)格見圖2。計算中，分別選取發(fā)動機燃燒室內(nèi)部沿軸線以及壁面位置進行數(shù)值監(jiān)測，并對時域信號進行分析，監(jiān)測點如圖3所示。

圖1　某固體火箭發(fā)動機工作末期計算區(qū)域及相應(yīng)邊界條件(對稱軸下側(cè)為改進后裝藥)

圖2　發(fā)動機噴管網(wǎng)格

圖3　發(fā)動機內(nèi)各監(jiān)測點位置

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

2.1發(fā)動機無過載情況

圖4所示為發(fā)動機內(nèi)部的壓力、馬赫數(shù)與溫度的瞬態(tài)結(jié)果，由于這些參數(shù)的空間脈動值相比總的壓力、馬赫數(shù)與溫度是小量，因而在圖4中觀察不出小的空間脈動分布?？傮w上看，這些參數(shù)的結(jié)果分布與穩(wěn)態(tài)計算的差別較小，圖中所示的結(jié)果符合發(fā)動機內(nèi)部流動的物理行為，燃氣在燃燒室內(nèi)積聚進行壓力填充，通過噴管流出時燃氣膨脹做功，壓力在噴管中下降、馬赫數(shù)升高、溫度下降。

圖5為監(jiān)測點的壓力隨時間變化情況。

圖4　發(fā)動機內(nèi)部燃氣的壓力云圖、馬赫數(shù)云圖與溫度云圖

(a)軸線上各監(jiān)測點

(b)壁面附近各監(jiān)測點

由于不同位置各點的脈動幅度不相同，為了顯示清楚脈動的變化規(guī)律，各圖縱坐標(biāo)的幅度是不同的。通過脈動幅度來看，各點沿著發(fā)動機軸線方向，脈動幅度逐漸增加。在相同的軸線位置，靠近發(fā)動機頭部部分，軸線與壁面處的脈動幅度相當(dāng)，但在靠近噴管的發(fā)動機尾部，壁面處的脈動幅度要大于軸線處的脈動幅度。從總的壓力脈動信號來看，脈動在幾千帕的范圍以內(nèi)，且并不規(guī)則。

為了觀測壓力脈動的頻率，圖6給出了壓力脈動的頻譜。傅立葉變換的結(jié)果顯示了壓力脈動具有一系列的特征頻譜，其中各點在60 Hz附近存在一個較為明顯的特征頻率。然而，該頻率卻并不固定，隨各點的軸線位置不同，頻率也有所變化，總的趨勢是在發(fā)動機軸向上往噴管方向頻率逐漸降低。

從各點的壓力譜同樣可觀察到，靠近發(fā)動機頭部位置，壓力振蕩幅值較小，軸線以及壁面位置處，其主頻位于111.82 Hz，見圖6中點P1～P4，而這一頻率在到達發(fā)動機中部后消失。另外，發(fā)動機頭部另一高頻位置位于62.99 Hz處，與主頻近似倍頻關(guān)系，見圖6中點P1～P6。隨著監(jiān)測點向噴管方向移動，壓力主要振蕩頻率變?yōu)?4.20 Hz，相對于這一頻率，此時62.99 Hz的振蕩幅值較小，見圖6中點P7～P12，且此時倍頻關(guān)系消失。除此之外，在低于20 Hz范圍內(nèi)，也有一系列低頻存在，對于尾部各點，由于60 Hz左右的頻率較強，這一頻率的譜密度要弱許多。各圖的縱坐標(biāo)的變化范圍也都不相同，總的趨勢是發(fā)動機越往后，其幅值越大。

(a)軸線上各監(jiān)測點　　　　(b)壁面附近各監(jiān)測點

圖7所示為各監(jiān)測點的速度頻域曲線。從圖7中可看出，除了發(fā)動機頭部的小幅度速度振蕩頻率與壓力振蕩頻率不同以外，速度譜在發(fā)動機后半段同樣捕捉到了主頻為54.199 22 Hz的速度振蕩，且速度譜在發(fā)動機壁面位置處，相對于軸線位置，其頻譜較復(fù)雜，說明壁面處的擾動強于軸線處。從振蕩幅值來看，靠近發(fā)動機噴管的壁面位置處，其振蕩幅值高于軸線上的振蕩幅值。

(a)軸線上各監(jiān)測點　　　　(b)壁面附近各監(jiān)測點

由于發(fā)動機裝藥結(jié)構(gòu)存在轉(zhuǎn)角，發(fā)動機內(nèi)出現(xiàn)了轉(zhuǎn)角旋渦脫落現(xiàn)象。圖8為發(fā)動機內(nèi)渦量在時間間隔為0.004 5 s的瞬態(tài)計算結(jié)果。由于噴管處流體的速度非常高，與噴管無滑移壁面之間形成極強的剪切流動，因而發(fā)動機內(nèi)部最大渦量值出現(xiàn)在噴管壁面附近，但這些極強的渦被壓縮在很小范圍內(nèi)，且隨高速流動流出發(fā)動機，對燃燒室聲腔影響最大的旋渦，還是在燃燒室內(nèi)部。

圖8　發(fā)動機內(nèi)瞬時渦量云圖(時間間隔為0.004 5 s)

圖8中，顯示了渦量從0～600(1/s)的范圍?？梢姡谘b藥尾部的突擴幾何結(jié)構(gòu)附近，產(chǎn)生了較復(fù)雜的流動結(jié)構(gòu)，裝藥的轉(zhuǎn)角結(jié)構(gòu)使得此處流動發(fā)生分離，產(chǎn)生較強較大的旋渦，這里產(chǎn)生的旋渦在主流的作用下，向下游移動發(fā)生脫落，但圖中脫落的渦在強度上受流體粘性影響，逐漸減弱，但其尺寸并沒有縮小，甚至有所增長，這是由軸對稱計算導(dǎo)致的湍流能量反向串級所造成的。結(jié)合圖5～圖7的壓力時間曲線以及壓力、速度的頻譜可發(fā)現(xiàn)，發(fā)動機內(nèi)轉(zhuǎn)角位置脫落的旋渦并非均勻一致的，但其脫落的頻率約為54.2 Hz。

另外，在突擴處下游靠近壁面處，也能觀察到有數(shù)量較多、尺度較小的漩渦脫落，這些小渦會與幾何擴張引起的大渦在軸對稱計算情況發(fā)生合并。需要特別指出的是，在固體燃面附近剪切流動的強度是隨空間變化的，越往下游剪切作用越強。由于當(dāng)前計算方法的限制，燃面的無滑移壁面本質(zhì)無法體現(xiàn)，實際的流量入口條件相當(dāng)于采用了一個有滑移壁面，這很大程度地降低了此處渦的強度和數(shù)量。實際的燃面附近應(yīng)比現(xiàn)在的計算結(jié)果產(chǎn)生更多更強的旋渦。

2.2軸向過載情況

為考察過載工作條件對發(fā)動機內(nèi)流場的影響，在發(fā)動機內(nèi)流場計算達到統(tǒng)計平衡態(tài)后，對流場施加23g的軸向過載條件，并對加載前后的流動情況進行對比。圖9所示為發(fā)動機內(nèi)各監(jiān)測點的壓力在加載條件下的變化情況。由圖9可知，加載條件使得靠近發(fā)動機頭部位置處的壓力發(fā)生階躍性下降，其壓力下降幅度可達到約0.1 bar(見圖9中P1～P4)。對于下降幅度較大的各點(如P1、P2)，其壓力經(jīng)過約0.2 s的振蕩后，逐漸趨于穩(wěn)定。隨著監(jiān)測點的位置向噴管喉部靠近，壓力的下降幅度逐漸減小，而此時過載條件對壓力的影響表現(xiàn)為壓力類似脈沖信號的振蕩，經(jīng)過同樣約0.2 s后，壓力的平均值趨于穩(wěn)定(見圖9中P7～P10，P12)。監(jiān)測點P11位于噴管收斂段，由圖9中發(fā)現(xiàn)，過載對P11處壓力的影響并不明顯。

(a)軸線上各監(jiān)測點　　　　(b)壁面附近各監(jiān)測點

圖10為加載后，發(fā)動機內(nèi)流場趨于穩(wěn)定后的壓力頻域曲線，與加載前的壓力頻譜相比，此時發(fā)動機頭部主頻位置是112.79 Hz，其次為61.04 Hz，見圖10中點P1～P4。同樣，隨著監(jiān)測點的位置逐漸向噴管移動，其壓力振蕩主頻下降為57.13 Hz，見圖10中點P5～P12。通過對比可發(fā)現(xiàn)，發(fā)動機軸向過載對靠近頭部位置的壓力小幅度振蕩的影響較小，而對于靠近發(fā)動機尾部受旋渦脫落影響較強的各點，過載條件使得低于60 Hz的振蕩頻寬增加，且使得靠近噴喉位置的振蕩幅度有較大地下降。

圖11對比了軸向過載分別為11.5g、23g、34.5g時，發(fā)動機內(nèi)不同位置處的壓力時間曲線。從圖11中可發(fā)現(xiàn)，加載瞬間壓力振蕩隨加載強度的增大而增大。

(a)軸線上各監(jiān)測點　　　　(b)壁面附近各監(jiān)測點

圖11　不同加載情況下監(jiān)測點P1、P11、P12位置處的壓力-時間曲線(圖中黑色實線為11.5 g，紅色點劃線為

2.3改變裝藥計算結(jié)果

圖12和圖13所示為改變裝藥結(jié)構(gòu)發(fā)動機工作末期燃燒室內(nèi)各監(jiān)測點壓力譜以及燃燒室渦量場，新裝藥方式的頭部采用快燃藥，則在其頭部出現(xiàn)了一個更大的空腔，燃面則相應(yīng)的減小，因而改變裝藥后的燃燒室壓力要更加小，馬赫數(shù)與溫度場同改變前差別較小。圖12中的壓力脈動頻譜顯示各點的主頻均為62 Hz，對比圖6中裝藥修改之前的頻譜曲線，此時燃燒室內(nèi)各點的主頻單一，相對于主頻，其他頻率并不明顯。

圖13顯示，改變裝藥的渦量被大大地降低，而且渦的尺寸明顯減小，呈現(xiàn)了較弱的流動不穩(wěn)定特征。這可能是因為此時加質(zhì)與改進前相比要更加少，剪切作用更弱。當(dāng)然，也有全尺寸實驗總結(jié)得出的結(jié)論指出，分段發(fā)動機中頭部空腔促進發(fā)動機不穩(wěn)定，因而對頭部空腔的影響，也值得更進一步細致地研究。

(a)軸線上各監(jiān)測點　　　　(b)壁面附近各監(jiān)測點

圖13　改變裝藥結(jié)構(gòu)發(fā)動機內(nèi)瞬時渦量云圖(時間間隔為0.008 s)

3　結(jié)論

(1)針對某固體火箭發(fā)動機工作過程中內(nèi)流場出現(xiàn)的旋渦流動現(xiàn)象，本文通過對發(fā)動機工作末期進行大渦模擬，得到了發(fā)動機后半段由于裝藥構(gòu)型產(chǎn)生的旋渦脫落現(xiàn)象，并得到了由于旋渦脫落產(chǎn)生的發(fā)動機內(nèi)流場壓力振蕩頻率。從壓力與速度的脈動幅度來看，沿發(fā)動機軸線方向脈動幅度逐漸增加。在相同的軸線位置，靠近發(fā)動機頭部部分，軸線與壁面處的脈動幅度相當(dāng)，但在靠近噴管的發(fā)動機尾部，壁面處的脈動幅度要大于軸線處的脈動幅度。

(2)在軸向過載條件下，過載條件使得發(fā)動機頭部壓力發(fā)生階躍性下降，而對于受旋渦脫落影響較大的發(fā)動機尾部附近，過載條件使得壓力產(chǎn)生較強脈沖式振蕩，振蕩的強度隨加載強度的增加而增大，同時過載條件使得發(fā)動機尾部壓力振蕩頻率有小幅度增長。

(3)發(fā)動機軸向過載時，能引起非常明顯的壓力脈動，本文的計算中，沒有考慮燃燒、過載時結(jié)構(gòu)振動等耦合因素，在過載一定時間后，過載所引起的壓力振蕩、速度振蕩影響逐漸消失。此時，處于過載情況的流動與無過載情況差別較小。最后，通過改變發(fā)動機頭部裝藥結(jié)構(gòu)，使得發(fā)動機尾部靠近壁面處的旋渦限度大幅減小。

[1]Blomshield F S.Historical perspective of combustion instability in motors-case studies[C]//37th Joint Propulsion Conference and Exhibit，AIAA 2001-3875，Salt Lake City，USA，2001.

[2]Ribéreau D，Pevergne T，Le Helley Ph，et al.Thrust oscillations passive control on large solid rocket motors[C]//41st Joint Propulsion Conference and Exhibit，AIAA 2005-4168，Tucson，Arizona，USA，2005.

[3]Dotson K W，Koshingoe 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)：197-206.

[4]Flandro G A，Jacobs H.Vortex generated sound in cavities[R].AIAA 73-1014.

[5]Kourta A.Computation of vortex shedding in solid rocket motors using time-dependent turbulence model[J].Journal of Propulsion and Power，1999，15(3)：390-400.

[6]Luppoglazoff N，Vuillot F.Parietal vortex shedding as a cause of instability for long solid propellant motors-numerical simulation and comparison with firing tests[R].AIAA 96-0761.

[7]Avalon G.Etude experimentale des instabilités naturelles d'un écoulement avec injection a la paroi[R].Campagne d'Essais VECLA ASSM7，Rapport Technique No.4/6163 DSNA/Y/DEFA，March 1999.

[8]Kourta A.Accoustic-mean flow interaction and vortex shedding in solid rocket motors[J].International Journal for Numerical Methods in Fluids，1996，22(6)：449-465.

[9]Lupoglazoff N，Vuillot F.Numerical simulation of vortex shedding phenomenon in 2D test case solid rocket motors[R].AIAA 92-0766.

[10]張翔宇，何國強，劉佩進.固體火箭發(fā)動機轉(zhuǎn)角渦脫落分類[J].航空動力學(xué)報，2014，29(8)：2003-2011.

[11]張嶠，李軍偉，王偉臣，等.固體火箭發(fā)動機渦聲耦合特征數(shù)值研究[J].推進技術(shù)，2011，32(3)：348-354.

[12]許團委，田維平，王建儒，等.固體發(fā)動機地面過載試驗?zāi)嗔Ｗ蛹铀俣确治鯷J].固體火箭技術(shù)，2014，37(4)：486-489.

(編輯：崔賢彬)

Internal flow in a solid rocket motor at the end of burning with axial acceleration

ZHANG Wei-ping1, YANG Dan2, XIONG Yong-liang3

(1.The 41st Institute of the Sixth Academy of China Aerospace Science & Industry Corporation,Huhhot010010,China;2.School of Naval Architecture and Ocean Engineering,Huazhong University of Science & Technology,Wuhan430074,China;3.School of Civil Engineering and Mechanics,Huazhong University of Science & Technology,Wuhan430074,China)

The pressure oscillation observed at the end of burning of solid rocket motor was studied numberically.The large eddy simulation (LES) was adopted to calculate the internal flow of solid rocket motor before it burns out.The angle vortex shedding was obtained due to the angular configuration of the propellant.By monitoring the flow parameters at different positons inside the motor,the frequency of the pressure oscillation due to the vortex shedding was captured.An acceleration condition at the axis direction was introduced to investigate its effect on the pressure oscillation.The results indicate that a pressure pulse was excited and the pressure changes are distinct at different position inside the motor due to the acceleration.The vortex shedding frequency increases slightly due to the overload condition.An improved propellant structure was adopted at the front of the motor and the results show that the vorticity in the end part of the chamber reduces significantly.

solid rocket motor;combustion instability;pressure oscillation;acceleration

2015-05-30；

2015-08-28。

國家自然科學(xué)基金(11372079，11502086，11502087)；中央高?；究蒲袠I(yè)務(wù)費專項資金(2015QN018， 2015QN141)。

張衛(wèi)平(1973—)，男，高級工程師，研究方向為發(fā)動機工作過程仿真。E-mail：yjy0372@sina.com

V435

A

1006-2793(2016)03-0306-06

10.7673/j.issn.1006-2793.2016.03.002