王 革，張 瑩，李冬冬，孫 娜

(1.哈爾濱工程大學 航天與建筑工程學院，哈爾濱 150001；2.上海航天動力技術研究所，上海 201109)

0 引言

由于工藝限制及發(fā)動機運輸問題，大型固體火箭發(fā)動機多采用分段式裝藥，每段裝藥用絕熱環(huán)來支撐。在發(fā)動機工作過程中，隨著燃面推移，絕熱環(huán)便成為主流中的障礙物，當主流流過絕熱環(huán)時便會形成障礙渦脫落。渦脫落作為發(fā)動機內一種常見的現(xiàn)象，受到越來越多的國內外專家學者的研究。

Flandro和Jacobs[1]于1975年提出發(fā)動機內的渦脫落現(xiàn)象可能是燃燒不穩(wěn)定[2]的另一源頭，渦脫落會導致燃燒室內壓力振蕩。早期的實驗和數(shù)值分析論證了潛入式噴管中流動和聲場的耦合，并且證明聲壓等級隨噴管空腔體積線性增加。Anthoine J等[3-4]用CPS模擬了流動和聲學耦合現(xiàn)象，結果表明空腔的存在會導致更高的壓力振蕩。Stella F等[5]對SRM中產生的壓力振蕩展開了數(shù)值研究，證明了空腔對流動渦影響的重要性。認識聲場結構異常重要，Jiang X[6]通過對亞聲速軸對稱射流內的聲場算法進行了研究，發(fā)現(xiàn)軸對稱Lilley方程的預測與聲場的DNS結果之間存在很好的一致性。劉佩進等[7]對分段SRM中渦聲耦合的實驗研究現(xiàn)狀進行了總結，并且通過研究發(fā)現(xiàn)頭部空腔、段間空腔以及潛入式噴管處空腔的存在都會對壓強振蕩產生較大影響。王健儒等[8]通過LES模型對不同燃面下的準穩(wěn)態(tài)內流場進行了模擬，發(fā)現(xiàn)在點火初期表面渦脫落占主導地位，隨著燃面逐漸推移，障礙渦脫落占主導。楊尚榮等[9]從流動穩(wěn)定性的角度進行數(shù)值模擬，并根據(jù)參考文獻[10-12]的理論方法進行理論求解，二者對比結果表明加質表面會產生表面渦脫落；此外，蘇萬興[13]對大長徑比SRM中不穩(wěn)定燃燒進行了線型預估，結果表明，隨著燃面的退移，壓力耦合響應增益與噴管阻尼均隨之減小。張嶠等[14]以VKI實驗發(fā)動機作基礎，通過LES方法揭示了SRM內的渦聲耦合機理，結果表明二維軸對稱大渦模擬方法基本滿足發(fā)現(xiàn)渦聲耦合規(guī)律的精度要求。李鵬飛[15]采用LES方法，通過改變障礙物的高度和位置、空腔大小和來流速度大小對固體火箭發(fā)動機中渦脫落現(xiàn)象展開研究，結果表明：障礙物越接近噴管，渦脫落頻率越大；當來流速度不斷增加時，渦脫落頻率也越大；渦脫落頻率與空腔體積成反比。Raheem[16]通過將燃燒室流動簡化成二維瞬態(tài)不可壓流動，借助Fluent軟件對Shanbhogue[17]的實驗進行模擬，李強等[18]又用大渦模擬方法針對Shanbhogue的部分實驗工況展開數(shù)值模擬，通過和實驗結果以及Raheem的計算結果進行對比發(fā)現(xiàn)：其計算結果與實驗基本相符，究其原因是渦結構與聲波有較強的三維性，通過二維模型很難模擬出其特點，并且LES能更準確地模擬渦脫落現(xiàn)象。張翔宇等[19]采用三維LES方法對固體火箭發(fā)動機C1x模型展開了數(shù)值模擬。楊羽卓[20]對帶潛入式噴管的含一個障礙的Ariane5縮比模型發(fā)動機展開了研究，結果表明，壓力場更能反映聲場特點，速度場與渦量場更多表現(xiàn)出的是流動特點。其還對渦脫落撞擊冷流縮比模型的噴管可能會導致的壓力振蕩展開了數(shù)值模擬[21]。

顯然，國內外專家學者對發(fā)動機內的渦脫落現(xiàn)象是非常關注的，但是這方面的研究局限于冷流實驗模型和準穩(wěn)態(tài)內流場模擬。本文就大長徑比固體火箭發(fā)動機內的障礙渦脫落現(xiàn)象展開了研究，重點分析了渦脫落過程中的壓力響應規(guī)律，結合開發(fā)的含包覆層界面識別方法，模擬了兩端包覆裝藥大長徑比發(fā)動機燃面動態(tài)推移過程，并且從準穩(wěn)態(tài)和燃面動態(tài)推移兩個角度展開了對比。

1 物理模型和計算方法

1.1 物理模型

聲學不穩(wěn)定性會引起內彈道壓力和推力振蕩，導致火箭發(fā)動機性能降低，以Ariane5助推器為例，Ariane5固體火箭發(fā)動機中分段燃燒室由三段推進劑組成，并采用絕熱環(huán)隔開，確保熱防護。參考文獻[3]，在本文中只采用了一個絕熱環(huán)，分析絕熱環(huán)障礙對聲場的影響。

參考冷流模型，基于含包覆層界面識別法，結合文獻[3]中信息設計出裝藥，以絕熱環(huán)暴露在主流中9 mm作為初始燃面，藥柱兩端被包覆，只有側面燃燒。本文設計的含裝藥大長徑比發(fā)動機幾何模型如圖1所示，陰影部分為藥柱，其余為流場區(qū)域，圖中尺寸單位為mm。其中，噴管喉部直徑為26 mm，此發(fā)動機的膨脹比為6.05。在燃面推移的過程中燃面附近不斷加質加熱，從而模擬燃面推移過程中燃氣的產生。

圖1 燃面推移幾何模型示意圖

一般研究大長徑比發(fā)動機內的渦脫落現(xiàn)象都采用準穩(wěn)態(tài)計算，從而得出不同型面下的內流場變化。雖然本文可以采用改進的LSPM方法實現(xiàn)燃面動態(tài)推移下的發(fā)動機內工作過程數(shù)值模擬，但是進行準穩(wěn)態(tài)內流場研究仍很有必要。對準穩(wěn)態(tài)內流場的分析可以對燃面動態(tài)推移結果起到一定參考作用。為了解準穩(wěn)態(tài)計算過程中的壓力響應規(guī)律，本文對圖1中裝藥推移到兩個時刻燃面下的流場進行了準穩(wěn)態(tài)計算，分別是初始燃面(型面A)和裝藥被燃去一半處燃面(型面B)，對應的計算模型如圖2所示。其中，型面A條件下絕熱環(huán)高度為9 mm，型面B條件下絕熱環(huán)高度為15 mm。

圖1、圖2中物理模型均采用結構網格，網格質量較高。為能夠更加準確地捕捉流場，在近壁面和絕熱環(huán)附近均進行加密，加密后的網格如圖3所示。為能將準穩(wěn)態(tài)計算結果與瞬態(tài)結果展開對比，均采用如表1所示的推進劑和燃氣參數(shù)。

1.2 計算方法

本文采用大渦模擬方法進行數(shù)值模擬，此方法本質是大渦直接求解，小渦用模型。小渦對大渦的影響用近似的模型體現(xiàn)，這稱為亞格子雷諾應力，文中采用WALE亞格子模型。大渦模擬的兩個重要步驟分別為濾波和建模。濾波方程如式(1)所示。

文中基于多孔介質模型[22]，結合Level-set方法實現(xiàn)燃面推移，通過添加質量、動量和能量源項模擬燃氣的生成。其中Level-set方程見式(2)，質量、動量和能量源項詳見式(3)～式(5)。

方法的具體細節(jié)參考文獻[22]。

后文中的計算涉及到裝藥兩端包覆的情況，為此基于參考文獻[22]的算法，開發(fā)了含包覆層復雜界面識別方法。

(a)型面A

(b)型面B

圖3 近壁面及障礙處網格加密

性能參數(shù)數(shù)值推進劑密度ρp/(kg/m3)1800燃速系數(shù)a0.006燃速指數(shù)n0.2比定壓熱容cp/[J/(kg·K)]1965.55相對分子量M30.303燃氣總溫T/K3532

注：對應的壓強單位為MPa。

(1)

(2)

(3)

(4)

(5)

圖4為復雜含包覆層的計算模型，只有側面推移，即只有一條界面隨時間推移，包覆層不沿著其中垂線向固體域內部推移，而是隨著界面的推移而消失。圖中固體域兩側為包覆層，在實際建模過程中包覆層厚度不體現(xiàn)，只是象征性的內部邊界，將流體域分成如圖所示無填充的4塊方便說明。將圖4中的1和2區(qū)域設為虛擬固體域，3是為了方便識別界面而劃分的流體區(qū)域，這三個區(qū)域即為改進LSPM的計算域，其大小可以隨意設定，但都是越小越好，能夠極大的縮短初始化和界面更新時間。(虛擬)固體域和流體域相交的部分即識別成初始界面。很明顯，界面不包括含包覆層的部分，但是卻多出了兩段界面，如圖4中紅色虛線所示。因此，燃面面積的計算方式需進行處理，處理的方式為只計算固體域和流體域3中識別到的界面。這種處理方式要求在網格劃分時將計算域分成更多的子區(qū)域，但提高了算法的適用范圍，并縮短了計算時間。

為避免中心差分格式引起的數(shù)值振蕩，動量方程采用BCD格式進行離散；連續(xù)性方程和能量方程采用二階迎風格式；為使速度與壓力能更好地滿足動量和連續(xù)性方程，采用PISO算法進行壓力速度耦合求解。

圖4 含包覆層界面識別模型示意圖

2 數(shù)值校驗

國內外研究渦脫落現(xiàn)象通常以冷流實驗為主，而本文旨在模擬實際發(fā)動機的工作過程中的渦脫落現(xiàn)象，發(fā)動機工作在高溫燃氣下，與冷流實驗的285 K不同，因此對冷熱流注入方式下的計算結果展開對比顯得尤為重要。本節(jié)數(shù)值校驗旨在研究流動溫度可能帶來的影響，為之后研究高溫燃氣下的流動過程作鋪墊。

2.1 發(fā)動機固有頻率

發(fā)動機能夠近似看做兩端封閉圓柱形腔體，聲在其中的傳播以一維縱向駐波為主，燃燒室內產生的渦在下游與噴管碰撞產生聲波，波經過傳播又在燃燒室頭部被反射回來，此反射波與原先的瞬時波干涉產生駐波[21]。因此，發(fā)動機可近似看做一個共振器，其縱向固有頻率計算公式如下：

f縱=nc/2L

(6)

式中f縱為縱向固有頻率；n為模態(tài)階數(shù)；c為當?shù)芈曀?；L為燃燒室長度。

當渦脫落頻率與固有頻率相近時，發(fā)動機內可能會發(fā)生渦聲耦合現(xiàn)象，渦脫落頻率會在各階聲頻間跳躍，并且渦和噴管收斂段的碰撞又會將聲壓提高到共振條件，并將能量以一階聲模態(tài)反饋至邊界層，把渦脫落頻率調整到對應聲頻，導致壓力振蕩被放大。

2.2 冷熱流注入結果分析

數(shù)值校驗的計算模型與圖2(a)一致，為冷流模型，參照文獻[3]采用頭部加質，以0.3 kg/s空氣頭部注入，而非圖中的側面加質。冷熱流數(shù)值模擬除了空氣溫度不同，其他計算設置均一致。其中冷流溫度為285 K，熱流溫度為3300 K。

對以冷、熱流頭部注入方式得出的內流場結果分別展開分析，計算過程中對頭部壓力進行監(jiān)測，圖5為冷熱流注入方式下頭部壓力振蕩曲線的FFT分析。從圖5可看出，冷流注入方式下，各階峰值及其對應的響應頻率較熱流注入條件下低很多。

(a)冷流頭部注入

(b)熱流頭部注入

從冷熱流注入條件下的數(shù)值分析結果可看出：熱流對應的壓力振蕩響應頻率集中在高頻，而冷流注入對應的壓力振蕩響應頻率集中在中頻。

冷熱流對比旨在研究音速對壓力波動響應頻率的影響，式(6)中已經交代了音速對固有頻率的影響，針對本節(jié)模型，分別計算出了冷熱流固有一階和二階頻率。表2為冷熱流數(shù)值計算中各階壓力峰值對應的響應頻率、冷流實驗值及發(fā)動機固有頻率的對比。從表2數(shù)據(jù)可看出：

(1)冷流數(shù)值模擬下壓力振蕩一階和二階峰值對應的頻率分別為406、840 Hz，與冷流實驗值基本吻合，認為數(shù)值計算結果可信。

(2)熱流數(shù)值模擬與冷流數(shù)值模擬之間唯一的區(qū)別就在于溫度，因此計算結果的不一致說明了溫度會對壓力響應頻率產生較大影響。從熱流數(shù)值模擬下壓力振蕩一階和二階響應頻率分別為1447、2935 Hz，可以推測出壓力響應頻率與溫度可能是開方倍關系，即壓力響應頻率與音速成正比。

(3)根據(jù)式(6)計算出的冷熱流條件下的發(fā)動機固有頻率，其大小與數(shù)值模擬結果比較接近，很容易產生渦聲耦合。但是未發(fā)現(xiàn)渦聲耦合現(xiàn)象，說明響應頻率與發(fā)動機固有頻率接近并不是產生渦聲耦合的必要條件。

表2 頭部壓力振蕩響應頻率對比

3 大長徑比SRM準穩(wěn)態(tài)內流場研究

針對1.1節(jié)中設計的型面A和型面B兩種物理模型，對大長徑比固體火箭發(fā)動機中準穩(wěn)態(tài)內流場展開了研究，方便與燃面動態(tài)推移過程中的瞬態(tài)內流場形成對比。

3.1 型面A條件下準穩(wěn)態(tài)內流場模擬擬

采用圖2(a)中計算模型，與燃面動態(tài)推移過程中加質方式一致，根據(jù)該型面下的燃面位置和大小進行側面加質，實現(xiàn)了型面A條件下的準穩(wěn)態(tài)內流場的數(shù)值模擬，通過對發(fā)動機頭部及凹腔位置壓力的監(jiān)測，對該型面下的壓力響應規(guī)律展開了分析。

分別對燃燒室頭部和凹腔監(jiān)測點壓力作FFT分析，其結果如圖6所示。從圖6可看出，發(fā)動機頭部監(jiān)測點平均壓力振蕩幅值在1000 Pa，而凹腔監(jiān)測點平均壓力振蕩幅值能達到2000 Pa。

從圖6統(tǒng)計出監(jiān)測點的前四階峰值對應的頻率見表3。從表3可見，雖然監(jiān)測點所處發(fā)動機中的位置不一致，但是前四階峰值對應的響應頻率相同，在型面A條件下的前四階壓力響應頻率分別為1100、2631、3970、5262 Hz。結合圖6中峰值大小發(fā)現(xiàn)，雖然不同位置處的各階響應頻率一致，但其對應的峰值大小不一。凹腔處頻譜分析的幅值更大。

3.2 型面B條件下準穩(wěn)態(tài)內流場模擬

采用圖2(b)中計算模型，與燃面動態(tài)推移過程中加質方式一致，根據(jù)該型面下的燃面位置和大小進行側面加質，實現(xiàn)了型面B條件下的準穩(wěn)態(tài)內流場的數(shù)值模擬，通過對發(fā)動機頭部及凹腔位置壓力的監(jiān)測，對該型面下的壓力響應規(guī)律進行了分析。

(a)發(fā)動機頭部

(b)凹腔處

階數(shù)頭部凹腔處階數(shù)頭部凹腔處一階11001100三階39703970二階26312631四階52625262

分別對監(jiān)測點壓力做FFT分析，其結果見圖7，表4為圖中各階峰值對應的響應頻率統(tǒng)計結果。結合圖7、表4中的信息可發(fā)現(xiàn)，發(fā)動機頭部和凹腔處的壓力頻譜分析曲線規(guī)律一致。表4中信息也表明兩個監(jiān)測點前四階響應頻率也相同，其中前四階峰值對應的各階響應頻率分別為1321、2541、3965、5337 Hz。與型面A條件下的準穩(wěn)態(tài)計算結果相同，凹腔處壓力振蕩平均幅值大約為2000 Pa，而發(fā)動機頭部壓力平均振蕩幅值只有1000 Pa，可見凹腔處的擾動較大，渦的運動也比較復雜。

(a)發(fā)動機頭部

(b)凹腔處

階數(shù)頭部凹腔處階數(shù)頭部凹腔處一階13211321三階39653965二階25412541四階53375337

4 大長徑比SRM燃面動態(tài)推移工作過程模擬

對準穩(wěn)態(tài)計算過程中渦脫落過程中壓力響應的分析是為研究含包覆層裝藥動態(tài)推移過程中壓力響應規(guī)律作鋪墊。研究燃面推移過程中瞬態(tài)內流場，是為探索燃面動態(tài)推移過程會對大長徑比發(fā)動機內渦脫落現(xiàn)象產生的影響，進而判定展開燃面動態(tài)推移過程中瞬態(tài)內流場計算的重要性。

與準穩(wěn)態(tài)求解模型不同，采用圖1計算模型對大長徑比固體火箭發(fā)動機工作過程中渦脫落現(xiàn)象開展研究，根據(jù)不同時刻燃面面積的大小實現(xiàn)燃面處的加質加熱，從而模擬燃氣的產生，以此來實現(xiàn)燃面推移過程中發(fā)動機內流場的計算。由于藥柱被包覆，且包覆面為凸曲面，運用自主開發(fā)的含包覆層界面識別方法可以很好地實現(xiàn)燃面推移。通過對發(fā)動機頭部及凹腔位置壓力的監(jiān)測，對燃面動態(tài)推移過程中的壓力響應規(guī)律展開了分析。

由于裝藥含包覆層，燃面隨著推移時間的增加而不斷變大，而計算時間步長較小，燃速又較低，所以整個裝藥燃燒完需要的計算量較大，即監(jiān)測點保存的數(shù)據(jù)也較多。為能夠更加直觀地看出頻譜分析的規(guī)律，對圖像進行了濾波處理，圖8即為濾波后監(jiān)測點的頻譜分析結果，表5給出了各階峰值對應的響應頻率統(tǒng)計結果。

(a)發(fā)動機頭部

(b)凹腔處

結合圖8、表5分析可發(fā)現(xiàn)，與傳統(tǒng)的準穩(wěn)態(tài)計算不同，瞬態(tài)計算伴隨著燃面的不斷推移，其中不僅有發(fā)動機型面及加質的不斷改變，還有時間疊加的因素包含在其中，因此整個燃面推移過程中的壓力響應與固定型面下的準穩(wěn)態(tài)計算有本質的區(qū)別。從圖8還可看出，發(fā)動機頭部和凹腔處壓力振蕩平均幅值較準穩(wěn)態(tài)計算結果均有大幅度下降。但是兩個監(jiān)測點的各階峰值對應的響應頻率相同，前四階響應頻率分別為1385、2620、4215、5527 Hz。凹腔處的壓力振蕩幅值比燃燒室頭部大，這與準穩(wěn)態(tài)計算得出的規(guī)律一致。

表5 燃面動態(tài)推移過程中不同監(jiān)測點的壓力振蕩響應頻率對比

5 瞬態(tài)與準穩(wěn)態(tài)計算結果對比

準穩(wěn)態(tài)計算已經成為國內外研究人員通用的一種計算發(fā)動機內流場的方式，其流場的準確性也已得到驗證。但裝藥燃面動態(tài)推移可能會對發(fā)動機內流場產生一定的影響，因此本節(jié)將會對燃面動態(tài)推移結果以及準穩(wěn)態(tài)結果展開對比，分析二者的相同點與區(qū)別。

燃面動態(tài)推移過程中關于壓力的頻譜分析是對發(fā)動機整個工作過程而言的，因此選取瞬態(tài)計算過程中與型面A和型面B相對應時刻的燃面展開分析，并且與準穩(wěn)態(tài)計算結果進行對比。

圖9為型面A、B條件下不同監(jiān)測點的瞬態(tài)和準穩(wěn)態(tài)壓力計算方式下的頻譜分析圖。表6統(tǒng)計出了兩種不同型面下的監(jiān)測點的瞬態(tài)和準穩(wěn)態(tài)一階響應頻率。結合圖9、表6可看出：

(1)燃面動態(tài)推移過程的模擬不會改變監(jiān)測點壓力頻譜分析曲線的變化趨勢，但是會令壓力峰值及其對應的響應頻率產生數(shù)值上的差異。此規(guī)律從側面驗證了瞬態(tài)計算結果的合理性，且具有一定的參考價值。

(2)瞬態(tài)計算對凹腔壓力振蕩產生的影響較大，型面A條件下瞬態(tài)計算凹腔處的一階峰值比準穩(wěn)態(tài)要高出2000 Pa，型面B條件下瞬態(tài)計算凹腔處的一階峰值比準穩(wěn)態(tài)要高出7000 Pa。

(3)瞬態(tài)和準穩(wěn)態(tài)壓力頻譜分析結果都表明，凹腔處的壓力振蕩幅值比發(fā)動機頭部要大，其中一階峰值在數(shù)值上尤為突出。

(4)初始燃面下的準穩(wěn)態(tài)計算中監(jiān)測點壓力一階峰值對應的頻率均為1100 Hz，該燃面下的瞬態(tài)計算結果中壓力一階峰值對應的頻率均為1150 Hz；型面B條件下的準穩(wěn)態(tài)和瞬態(tài)計算中監(jiān)測點一階峰值對應的頻率分別為1321 Hz和1200 Hz。結果表明，燃面推移帶來的一階響應頻率的改變是非線型的。

(a)型面A，頭部

(b)型面A，凹腔處

(c)型面B，頭部

(d)型面B，凹腔處

型面監(jiān)測點準穩(wěn)態(tài)瞬態(tài)A頭部11001150凹腔11001150B頭部13211200凹腔13211200

6 結論

(1)溫度(音速)會對發(fā)動機內的壓力響應頻率產生較大影響，其響應頻率的大小與溫度呈開方倍關系。

(2)針對發(fā)動機內同一位置監(jiān)測點，與準穩(wěn)態(tài)計算結果相比，燃面動態(tài)推移過程中的壓力響應幅值較小，大長徑比發(fā)動機工作過程中不僅有由于燃面推移導致的型面不斷改變，且還有較長的時間累積過程，因此壓力頻譜分析中幅值有大幅度下降。

(3)通過對不同型面下的瞬態(tài)和準穩(wěn)態(tài)FFT結果展開對比，其結果表明，凹腔處的壓力振蕩幅值比發(fā)動機頭部要大。瞬態(tài)計算不會改變監(jiān)測點壓力頻譜分析曲線的變化趨勢，但是會令壓力峰值及其對應的響應頻率產生數(shù)值上的差異，且燃面推移帶來的一階響應頻率的改變是非線型的。