固體火箭發(fā)動機聲渦耦合數(shù)值仿真*

韓磊，卞云龍，張衛(wèi)平，孫再庸

(1.內(nèi)蒙動力機械研究所，內(nèi)蒙古 呼和浩特　010010； 2. 中國航天科工集團第二研究院，北京　100038)

?

仿真技術(shù)

固體火箭發(fā)動機聲渦耦合數(shù)值仿真*

韓磊1,2，卞云龍1，張衛(wèi)平1，孫再庸1

(1.內(nèi)蒙動力機械研究所，內(nèi)蒙古 呼和浩特010010； 2. 中國航天科工集團第二研究院，北京100038)

摘要：研究固體火箭發(fā)動機工作中發(fā)生的聲渦耦合導(dǎo)致的不穩(wěn)定燃燒，以及燃燒室空腔變化對壓力振蕩的影響。對VKI縮比發(fā)動機流場采用大渦模擬(LES)方法進行數(shù)值模擬，確定數(shù)值方法可行。采用有限元和大渦模擬結(jié)合的方法，對某固體火箭發(fā)動機聲渦耦合進行仿真分析，對聲場分別采用理論方法和有限元數(shù)值方法對聲模態(tài)和聲學(xué)特性計算，確定有限元方法精度可靠及大渦模擬方法對流場旋渦流動不穩(wěn)定捕捉準(zhǔn)確，計算結(jié)果和試驗中壓力振蕩結(jié)果相符。表明研究的該固體火箭發(fā)動機發(fā)生了聲渦耦合引起的縱向1階聲不穩(wěn)定，且獲得了自由容積對壓力振蕩頻率及幅值的影響規(guī)律。

關(guān)鍵詞：聲渦耦合；不穩(wěn)定燃燒；壓力振蕩；數(shù)值模擬；聲模態(tài)；有限元；大渦模擬

0引言

發(fā)動機燃燒工作過程中，由于某種隨機擾動影響燃燒速率，燃燒室內(nèi)流場和燃燒過程、空腔的反射、聲場等相互作用，引起的不規(guī)則的燃燒室壓力變化，稱為不穩(wěn)定燃燒，也常稱為振蕩燃燒[1]。壓力振蕩會導(dǎo)致較嚴(yán)重的推力振蕩，輕則發(fā)動機振動，中斷燃燒，嚴(yán)重的情況下，和飛行器耦合共振，導(dǎo)致系統(tǒng)可靠性降低、失效乃至災(zāi)難性的爆炸事故發(fā)生[2-3]。

近年來，隨著固體運載要求的提出以及導(dǎo)彈技術(shù)發(fā)展，對固體火箭發(fā)動機提出了高裝填、大推力比、初始大推力、大長細(xì)比等要求。國內(nèi)的發(fā)動機設(shè)計者為了滿足高性能指標(biāo)，提高設(shè)計標(biāo)準(zhǔn)，導(dǎo)致戰(zhàn)術(shù)地空導(dǎo)彈用固體火箭發(fā)動機屢次出現(xiàn)燃燒不穩(wěn)定現(xiàn)象。除此之外，國外研制的大型助推器，如美國的SpaceShuttle和歐洲的Ariane5[4-5]采用的大型分段式固體發(fā)動機相繼出現(xiàn)了較為嚴(yán)重的不穩(wěn)定燃燒，不穩(wěn)定燃燒研究也越來越重要[6]。

渦脫落產(chǎn)生的聲學(xué)振蕩被認(rèn)為是不穩(wěn)定燃燒問題中一個關(guān)鍵問題，最初Flandro和Jacobs[7]提出由于渦脫落激發(fā)的聲模態(tài)是引起不穩(wěn)定燃燒的一個重要因素，早期證據(jù)表明復(fù)雜幾何裝藥的戰(zhàn)術(shù)導(dǎo)彈上面級發(fā)動機的不穩(wěn)定燃燒由旋渦脫落驅(qū)動。后來有研究表明聲渦耦合導(dǎo)致的聲加強剪切層不穩(wěn)定的現(xiàn)象會放大了壓力震蕩，當(dāng)渦脫落頻率與聲振頻率相同時會形成了不穩(wěn)定燃燒[8]。

1問題描述

本文同時采用理論計算和有限元結(jié)合的方法，針對在工作末期發(fā)生了不穩(wěn)定燃燒的某固體火箭發(fā)動機不同燃燒時刻燃燒室聲腔進行分析和數(shù)值計算。并對該發(fā)動機翼柱燒盡、振蕩前、振蕩時、振蕩后等關(guān)鍵工作時刻建模，對不同時刻燃燒室流場采用大渦模擬方法數(shù)值模擬，初步分析認(rèn)為由于流動不穩(wěn)定性產(chǎn)生了渦脫落，導(dǎo)致聲渦耦合，從而產(chǎn)生不穩(wěn)定燃燒，研究該不穩(wěn)定燃燒發(fā)生時刻是否發(fā)生聲渦耦合，對耦合情況下振蕩特性分析，并研究工作過程中自由容積變化對振蕩特性的影響。

2數(shù)值方法與算例校驗

2.1聲學(xué)共振頻率理論方法

燃燒室固有聲學(xué)頻率計算公式[9]:

(1)

(2)

(3)

式中：n為模態(tài)階數(shù)；c為當(dāng)?shù)芈曀佟?/p>

2.2聲學(xué)共振頻率有限元方法

采用基于有限元的方法,直接離散三維波動方程,可得到復(fù)雜裝藥聲腔的聲學(xué)特性。使用迦遼金法對簡化的亥姆霍茲方程進行離散,聲腔模態(tài)方程的單元矩陣[10]形式為

(4)

式中：Kf為聲剛度矩陣；Mf為聲質(zhì)量矩陣；p為聲特征向量；ωa為特征根。模型表面定義0位約束，使用有限元法求得圓周頻率ωa進而求得聲振頻率fa。

2.3大渦模擬方法

通過在傅里葉空間或構(gòu)型空間將隨時間變化的N-S方程進行濾波可得到控制方程。考慮到氣體的可壓縮性，利用Favre平均對控制方程按式(5)簡化。

(5)

本文不考慮化學(xué)反應(yīng)，僅計算單組分工質(zhì)，濾波后連續(xù)方程、動量方程與能量方程分別為

(6)

(7)

(8)

式中：“-”表示Reynolds平均;“～”表示Favre平均。

(9)

(10)

(11)

(12)

式中：常數(shù)Cε為1.1。

2.4大渦模擬算例校驗

為研究障礙物導(dǎo)致渦脫落，VKI實驗室設(shè)計了含有絕熱環(huán)和潛入式噴管空腔縮比冷流實驗發(fā)動機，Anthoine[12]還對該發(fā)動機進行試驗研究。流場計算區(qū)域尺寸結(jié)構(gòu)如圖1所示。VKI發(fā)動機尺寸已知，且經(jīng)過冷流實驗和數(shù)值模擬驗證，可以對本文采用的大渦模擬方法進行校驗，為此選用此發(fā)動機進行方法校驗。

圖1　VKI幾何模型Fig.1　VKI model and domain

旋渦具有三維特性，產(chǎn)生、拉伸過程均為三維空間過程。二維模型計算對旋渦的捕捉和模擬總有欠缺，但三維模擬計算量太大，Mason[13]和張嶠[14]等人采用二維和三維大渦模擬對比后表明二維模擬可以成功的對ETM-03，RSRM，VKI等發(fā)動機進行數(shù)值模擬，且對振蕩頻率等特性很好的預(yù)估，為此本文綜合計算精度和效率采用二維模型。

采用結(jié)構(gòu)化網(wǎng)格，結(jié)構(gòu)化網(wǎng)格質(zhì)量好，易于計算收斂；方便控制局部區(qū)域網(wǎng)格稀疏，易于對流動復(fù)雜區(qū)域和關(guān)注區(qū)域網(wǎng)格加密；減少其他區(qū)域網(wǎng)格，從而減少計算時間。此處對障礙區(qū)域、產(chǎn)生旋渦流動區(qū)域、旋渦和噴管頭部碰撞區(qū)域進行網(wǎng)格加密，局部加密網(wǎng)格如圖2所示。

圖2　局部網(wǎng)格加密Fig.2　Refined mesh zone

在發(fā)動機流場內(nèi)選取關(guān)鍵點監(jiān)測壓力波動，入口為質(zhì)量流率入口，溫度為285K，質(zhì)量流率為66.3kg/m2的理想氣體，采用軸向進氣[15],比熱比為1.4，Pr數(shù)取0.71，出口為壓力出口，壓力為101 325Pa，溫度為285K。流場初始計算先用穩(wěn)態(tài)計算得到一個初始流場，待流場穩(wěn)定后，改湍流模型為大渦模擬，為了避免連續(xù)方程和動量方程中心差分格式產(chǎn)生數(shù)值振蕩，采用BCD格式進行離散，并同時改為時間步長為4e-7的非穩(wěn)態(tài)計算，待流場穩(wěn)定開始監(jiān)測頭部壓力。

圖3為部分時間段內(nèi)發(fā)動機頭部監(jiān)測點壓力曲線，由圖可以看出，發(fā)動機內(nèi)流場存在明顯的壓力振蕩，壓力振蕩區(qū)間為(170 000～210 000Pa)，對壓力曲線進行快速傅里葉變換后頻率分布圖如圖4，在0～2 000BZ頻率內(nèi)，其中振蕩主頻為424BZ，壓力振蕩主要振蕩頻率為424,848,1 272,1 697BZ，呈明顯倍頻，主頻壓力振幅是其他振蕩頻率下振幅的7倍。

圖3　壓力振蕩曲線Fig.3　Pressure-time oscillation curve

圖4　頻譜分析Fig.4　Pressure spectrum analysis

Anthoine等人對此發(fā)動機分別進行過數(shù)值模擬和試驗，本文得到的振蕩頻率與其試驗和數(shù)值模擬結(jié)果對比如圖5和表1所示，四階主頻與Anthoine實驗和數(shù)值結(jié)果吻合，其中一階模態(tài)誤差最大，最大誤差僅為3.4%，表明本文所采用的大渦模擬方法對研究旋渦脫落導(dǎo)致流場不穩(wěn)定可行。

圖5　不同研究方法下頻率對比Fig.5　Comparison of frequencies in differentmethods

模態(tài)階數(shù)LES結(jié)果/BZ實驗結(jié)果/BZ誤差(%)14244103.428488753.13127212800.64169717402.5

3數(shù)值方法與算例校驗

3.1計算模型

本文研究的發(fā)動機為大長徑比固體火箭發(fā)動機，裝藥為前內(nèi)孔后翼，為此選取以下4個關(guān)鍵工作時刻進行研究，發(fā)動機在6s時后翼燒盡，燃燒室空腔變化趨于均勻，在15.5s左右發(fā)生振蕩燃燒，后翼燒盡后很長一段時間均會產(chǎn)生壓力振蕩，壓力振蕩頻率為173BZ。在6～15.5s工作時間段內(nèi)燃燒室內(nèi)自由容積隨著燃面的退移產(chǎn)生顯著的變化，會導(dǎo)致燃通比等參數(shù)明顯變化。為此，本節(jié)主要以后翼燃燒結(jié)束時刻6s，振蕩發(fā)生前的時刻14s，壓力振蕩發(fā)生開始時刻15.5s，振蕩后期時刻16s，4個時刻的發(fā)動機建模，計算模型如圖6。

圖6　關(guān)鍵時刻幾何模型Fig.6　Key models and domain

壓力監(jiān)測點在流場區(qū)域面積最小的16s時刻模型中選取，監(jiān)測點位置分別為：燃燒室頭部，突擴前，突擴后，突擴段中間，突擴段尾部，如圖7所示。

圖7　壓力監(jiān)測點分布Fig.7　Distribution of pressure monitoring points

3.2邊界條件與計算參數(shù)

燃面在不同時刻略有變化，燃燒室工作壓強在不同燃燒時刻也不同，這些因素在本文研究中不是主要因素，為此忽略燃面、壓強的細(xì)微不同，所有時刻燃面入口統(tǒng)一設(shè)置為質(zhì)量流率入口，質(zhì)量流量均為15.5s時刻的預(yù)估質(zhì)量流率，出口為壓力出口101 325Pa，出口溫度300K，燃?xì)鈪?shù)見表2。噴管出口為超聲速，出口界面參數(shù)由外推得到，且由于氣體在噴管加速至超聲速，下游對上游無影響，壁面采用無滑移壁面。

表2　燃?xì)鈪?shù)

3.3燃燒室聲振頻率理論計算

應(yīng)用式(1)計算燃燒室聲振頻率如表3所示。

表3　理論計算縱向頻率

3.4燃燒室聲學(xué)模態(tài)有限元仿真計算

理論計算聲腔聲模態(tài)和振型時要求聲腔幾何構(gòu)型簡單，且邊界條件明確，鑒于許多燃燒室初始聲腔極復(fù)雜，本文中燃燒室相對簡單，但依然存在階梯凸臺等特征，為此對燃燒室三維聲腔建立有限元模型，噴管部分超聲速，相當(dāng)于燃燒室兩端封閉，只選取燃燒室主體部分，流體單元選取fluid30，模型表明定義零位移約束，取聲介質(zhì)密度為4.0kg/m3,平均聲速為1 061m/s。

圖8　各階振型聲壓分布Fig.8　Acoustic pressure distribution of every order axial vibration modes

圖8中為某固體火箭發(fā)動機發(fā)生不穩(wěn)定燃燒時前4階縱向振型聲壓分布云圖，前4階聲振頻率如表4所示，一階縱向頻率165BZ與試驗發(fā)動機壓力振蕩頻率相近，聲壓分布中，發(fā)動機頭尾部為聲壓波腹，表明發(fā)生了與一階縱向聲頻耦合的聲渦耦合。圖9為有限元方法和理論計算方法下聲振頻率對比，對比結(jié)果表明隨著模態(tài)階數(shù)增大，兩者計算誤差也有所增大，最大誤差約為5%。

表4　有限元計算聲振頻率

圖9　不同研究方法下頻率對比Fig.9　Comparison of frequencies with differentmethods

3.5燃燒室各時刻壓力振蕩分析

在發(fā)動機工作6s后翼柱燒盡，此后發(fā)動機內(nèi)腔變化趨于均勻，不同時刻計算入口條件一致的情況下，待監(jiān)測壓力波動穩(wěn)定后，觀測得知從6s到最后16s過程中均會產(chǎn)生一定的壓力振蕩。6s時壓力在波動中下降，14，16s時刻壓力振蕩在波動中上升，分別對以下4個時刻壓力波動周期平均計算后發(fā)現(xiàn)壓力波動周期處于5.86～6.26ms之間，壓力波動頻率均處于159～161BZ之間。若對壓力波動濾去波動量后， 14，15.5，16s時壓力波動類似正弦波動，15.5s時的壓力波動頻率如圖10，周期性明顯。盡管都存在壓力振蕩，可以看出14，15.5s時刻壓力振蕩幅值較大，15.5s壓力振蕩相比其他時刻明顯高一個數(shù)量級，6s時由于后階凸臺明顯，會產(chǎn)生明顯的轉(zhuǎn)角渦脫落，可能會導(dǎo)致頭部壓力波動較大，14s時的壓力波動幅值是前2個時刻的10倍左右。15.5s壓力振蕩幅值遠(yuǎn)大于其他時刻，是前兩時刻壓力波動20倍左右。模擬結(jié)果得到的15.5s的壓力振蕩最大，和發(fā)動機實際工作情況相符。

圖10　15.5 s時刻壓力振蕩Fig.10　Pressure fluctuation at 15.5 s

監(jiān)測點head處的平均壓力隨發(fā)動機工作時間變化如圖11所示，隨著發(fā)動機工作，自由容積的變大，平均壓力先逐漸降低，在振蕩發(fā)生時刻附近上升，15.5s時最大，然后下降，平均壓力的最大波動范圍在5%左右。

對壓力時間曲線，進行開窗快速傅里葉變換后，得到頻率、振幅曲線， 均出現(xiàn)壓力振蕩， 如圖12所示，15.5s時振幅最大，波峰對應(yīng)的前2階頻率分別為167，332BZ，其167BZ與一階縱向聲頻161BZ和165BZ接近。隨著時間變化，燃燒室自由容積增大，1階壓力振蕩頻率幾乎不變，其他頻率峰值逐漸平移，待發(fā)生壓力和聲振耦合共振后，隨著時間推移，產(chǎn)生共振的壓力振蕩頻率消失，到燃燒室工作末期，壓力振蕩幾乎消失。在燃面藥形設(shè)計中避開14～15.5s時刻的變化壓力振蕩明顯時的空腔，即可只發(fā)生低頻、小振幅微弱振蕩，而避免發(fā)生共振。

圖11　不同時刻平均壓力Fig.11　Average pressure at different times

圖12　不同時刻振蕩頻率和振幅Fig.12　Pressure spectrum at different times

4結(jié)束語

大渦模擬和有限元方法結(jié)合的方法可以對某火箭發(fā)動機數(shù)值仿真與分析，表明該發(fā)動機發(fā)生了縱向1階頻率主導(dǎo)的聲渦耦合聲不穩(wěn)定燃燒，壓力振蕩的時刻特性與發(fā)動機試車工作情況相符。隨著發(fā)動機工作，自由容積增大的過程中，平均壓力先增大最然后下降，在振蕩發(fā)生時達到峰值，平均壓力的最大波動范圍在5%左右。表明本文方法可以對聲渦耦合引起的不穩(wěn)定燃燒情況進行預(yù)計，可以為發(fā)動機設(shè)計和抑振提供參考校驗。

參考文獻：

[1]孫維申.固體火箭發(fā)動機不穩(wěn)定燃燒[M].北京：北京工業(yè)學(xué)院出版社，1987.

SUNWei-shen.CombustionInstabilityinSolidRocketMotors[M].Beijing:BeijingInstituteofTechnologyPress, 1987.

[2]劉佩進，金秉寧，李強.戰(zhàn)術(shù)導(dǎo)彈固體發(fā)動機燃燒不穩(wěn)定研究概述[J].固體火箭技術(shù)，2012，35(4)：446-456.

LIUPei-jin,JINBing-ning,LIQiang.ASurveyofCombustionInstabilityinTacticalSRM[J].JournalofChinaOrdance，2012，35(4)：446-456.

[3]王寧飛，張嶠，李軍偉，等.固體火箭發(fā)動機不穩(wěn)定燃燒研究進展[J].航空動力學(xué)報,2011,26(6)：1406-1412.

WANGNing-fei,ZHANGQiao,LIJun-wei,etal.ProgressofInvestigationonCombustionInstabilityofSolidRocketMotors[J].JournalofAerospacePower,2011,26(6):1406-1412.

[4]BLOMSHIELDFS.HistoricalPerspectiveofCombustionInstabilityinMotors:CaseStudies[R].AIAA2001-3875.

[5]胡大寧，何國強，劉佩進, 等. 翼柱型藥柱固體火箭發(fā)動機不穩(wěn)定燃燒研究[J].固體火箭技術(shù)，2010，33(5)：502-506.

HUDa-ning,HEGuo-qiang,LIUPei-jin,etal.StudyonInstableCombustionofSolidRocketMotorwithFinocylGrain[J].JournalofChinaOrdance,2010, 33(5)：502-506.

[6]BLOMSHIELDFS.SummaryofMulti-DisciplinaryUniversityResearchInitiativeinSolidPropellantCombustionInstability[R].AIAA2000-3172.

[7]FLANDROGA,JACOBSHR.Vortex-GeneratedSoundinCavities[R].AIAA1973-1014.

[8]熊永亮，楊丹，郜冶，等.固體火箭發(fā)動機燃燒室內(nèi)聲波振蕩數(shù)值研究[J].推進技術(shù),2006,27(6):484-487.

XIONGYong-liang,YANGDan,GAOYe,etal.NumericalStudyofAcousticOscillationinSolidRocketMotors[J].JournalofPropulsionTechnology,2006,27(6):484-487.

[9]唐金蘭.固體火箭發(fā)動機原理[M].北京：國防科技出版社,2013.

TANGJin-lan.SolidRocketMotorMechanical[M].Beijing:NationalDefenseTechnologyPress,2015.

[10]楊向明,劉佩進,陳曉龍.翼柱型裝藥固體火箭發(fā)動機燃燒室聲場分析[J].宇航學(xué)報，2008，29(5):1593-1597.

YANGXiang-ming,LIUPei-jing,CHENXiao-long.AnalysisofAcousticPropertywithComplicatedGrainShapeinCombustionChamberofSolidRocketMotor[J].JournalofAstronautics, 2008, 29(5): 1593-1597.

[11]NICOUDF,DUCROSF.Subgrid-ScaleStressModelingBasedontheSquareoftheVelocityGradientTensor[J].Flow,TurbulenceandCombustion,1999,62(3):183-200.

[12]ANTHOINEJ,BUNCHLINJM,GUERYJF.EffectofNozzleCavityonResonanceinLargeSRMNumericalSimulation[J].JournalofPropulsionandPower,2003,19(3);374-384.

[13]MASONDR,MORSTADTRA,CANNONSM.PressureOscillationandStructuralVibrationsinSpaceShuttleRSRMandETM-3Motors[R].AIAA2004-3898,2004

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

ZHANGQiao，LIJun-wei，WANGWei-chenetal.NumericalAnalysisonVortexAcousticCouplingCharacteristicsinSolidRocketMotors[J].JournalofPropulsionTechnology,2011,32(3):349-354.

[15]STELLAF,PAGLIAF.PressureOscillationinSolidRocketMotorsNumericalStudy[J].AerospaceScienceandTechnology,2011,15(1):53-59.

Numerical Simulation on Acoustic Vortex Coupling in Solid Rocket Motors

HAN Lei1,2, BIAN Yun-long1, ZHANG Wei-ping1, SUN Zai-yong1

(1. Dynamic Machinery Institute of Inner Mongolia, Inner Mongolia Hohhot 010010, China；2.TheSecondAcademyofChinaAerospaceandIndustryCorporation,Beijing100038,China)

Abstract:In order to study solid rocket motors (SRM) combustion instability induced by acoustic vortex coupling and the effect of combustion chamber free volume on pressure oscillation. Large-eddy-simulation (LES) is used to simulate flow field domain of VKI scale down motor to show that this numerical method is capable of application. A method combining finite-element-method (FEM) and LES is applied to simulate and analyze acoustic vortex coupling for a certain motor. Theoretical method and FEM are used to calculate acoustic modes and acoustic characteristics, and it shows that the accuracy of FEM is reliable. For catching of the flow instability detail, LES is accurate. Pressure oscillation of simulation is consistent with the experiment result. And the result shows that the first-order axial mode acoustic instability induced by acoustic vortex coupling occurs in the SRM. The effect of the free volume on pressure oscillation is also acquired.

Key words:acoustic vortex coupling；unstable combustion；pressure oscillation ；numerical simulation；acoustic mode；finite-element-method(FEM)；Large-eddy-simulation(LES)

*收稿日期：2015-05-20；修回日期：2015-07-06

作者簡介：韓磊(1989-)，男，寧夏銀川人。碩士生，主要從事固體火箭發(fā)動機數(shù)值仿真。

通信地址：010010內(nèi)蒙古呼和浩特市新華東街65號航天科工六院41所E-mail：docarezone@qq.com

doi:10.3969/j.issn.1009-086x.2016.03.030

中圖分類號：V435；TP391.9

文獻標(biāo)志碼：A

文章編號：1009-086X(2016)-03-0194-07