陳俊屹, 王 兵, 田小濤

(1.清華大學(xué)， 北京 10084； 2.西安現(xiàn)代控制技術(shù)研究所，西安 710065)

【航天工程】

某型固體火箭發(fā)動(dòng)機(jī)燃燒不穩(wěn)定性數(shù)值模擬研究

陳俊屹1, 王 兵1, 田小濤2

(1.清華大學(xué)， 北京 10084； 2.西安現(xiàn)代控制技術(shù)研究所，西安 710065)

針對(duì)某工程單位某改性雙基裝藥長(zhǎng)尾管固體火箭發(fā)動(dòng)機(jī)試車過程，建立了發(fā)動(dòng)機(jī)一維模型，在一維準(zhǔn)靜態(tài)燃燒模型的基礎(chǔ)上，對(duì)發(fā)動(dòng)機(jī)燃燒與聲流場(chǎng)耦合的過程進(jìn)行了數(shù)值模擬。計(jì)算結(jié)果表明，低溫下壓力耦合的響應(yīng)函數(shù)峰值頻率恰好處于燃燒室固有頻率附近，誘發(fā)了不穩(wěn)定燃燒。通過提高燃速，使得壓力耦合的響應(yīng)函數(shù)頻率峰值偏離燃燒室固有頻率，可以有效地抑制不穩(wěn)定燃燒的出現(xiàn)。

固體火箭發(fā)動(dòng)機(jī)；燃燒不穩(wěn)定性；數(shù)值模擬；改性雙基裝藥

隨著導(dǎo)彈總體對(duì)固體火箭發(fā)動(dòng)機(jī)性能要求的提高，在工程上更多地采用高能推進(jìn)劑。隨著推進(jìn)劑能量的提高，近年來發(fā)動(dòng)機(jī)頻繁出現(xiàn)燃燒不穩(wěn)定現(xiàn)象。固體火箭發(fā)動(dòng)機(jī)燃燒穩(wěn)定性，在工程中大多采用后驗(yàn)的方法，實(shí)驗(yàn)修正的周期長(zhǎng)、成本高。對(duì)于發(fā)動(dòng)機(jī)燃燒穩(wěn)定性的預(yù)估，持續(xù)受到固體發(fā)動(dòng)機(jī)研究工作者的關(guān)注。發(fā)動(dòng)機(jī)燃燒穩(wěn)定性預(yù)估方面已經(jīng)有許多已經(jīng)完成的工作[1-3]。在工程實(shí)踐中處理發(fā)動(dòng)機(jī)燃燒不穩(wěn)定性的方法，多為線性理論預(yù)估與實(shí)驗(yàn)相結(jié)合?？梢园l(fā)現(xiàn)，對(duì)于發(fā)動(dòng)機(jī)整機(jī)燃燒不穩(wěn)定性的數(shù)值模擬研究較少。因此有必要發(fā)展一種發(fā)動(dòng)機(jī)整機(jī)數(shù)值模擬方法，為更深入的穩(wěn)定性預(yù)估數(shù)值研究方法提供基礎(chǔ)。

本文發(fā)展了一種固體火箭發(fā)動(dòng)機(jī)理論預(yù)估的數(shù)學(xué)模型與數(shù)值方法。針對(duì)某型固體火箭發(fā)動(dòng)機(jī)試車過程中出現(xiàn)的燃燒不穩(wěn)定性現(xiàn)象進(jìn)行分析，提出了不穩(wěn)定燃燒的抑制方法。

1 數(shù)學(xué)模型與數(shù)值方法

1.1 數(shù)學(xué)模型

1.1.1 發(fā)動(dòng)機(jī)幾何模型

針對(duì)某工程單位某型發(fā)動(dòng)機(jī)出現(xiàn)的縱向聲振燃燒不穩(wěn)定，將發(fā)動(dòng)機(jī)化簡(jiǎn)為一維。其燃燒室內(nèi)流場(chǎng)幾何構(gòu)型如圖1所示。

圖1 發(fā)動(dòng)機(jī)幾何構(gòu)型示意圖

圖1中所示的整個(gè)左端面為推進(jìn)劑燃面，右端面為噴管出口。發(fā)動(dòng)機(jī)總長(zhǎng)為406 mm。燃燒室長(zhǎng)度約為167 mm。

1.1.2 氣相控制方程

忽略發(fā)動(dòng)機(jī)中燃?xì)獾恼承耘c凝聚相，認(rèn)為燃?xì)獾母黜?xiàng)熱力學(xué)常數(shù)不隨溫度變化。推進(jìn)劑的燃燒過程通過源項(xiàng)的形式添加到方程中。發(fā)動(dòng)機(jī)氣相控制方程如下：

(1)

燃?xì)獾臒崃W(xué)參數(shù)關(guān)系如下：

(2)

方程的左邊界為固壁邊界條件，右邊界為無反射邊界條件。方程的初始條件為流場(chǎng)穩(wěn)態(tài)解，并在左端部施加一個(gè)1%室壓的壓力脈沖。

1.1.3 固相導(dǎo)熱控制方程

(3)

其初始條件由穩(wěn)態(tài)解決定，為：

(4)

非穩(wěn)態(tài)邊界條件為：

(5)

1.1.4 一維準(zhǔn)靜態(tài)燃燒模型

在本文采用的一維準(zhǔn)靜態(tài)燃燒模型如圖2所示，模型中將燃燒的推進(jìn)劑分為固相、氣相和反應(yīng)薄層三個(gè)部分[4,5]。其中反應(yīng)薄層包含氣相反應(yīng)區(qū)與固相反應(yīng)區(qū)。對(duì)于氣相反應(yīng)區(qū)，認(rèn)為反應(yīng)的特征時(shí)間遠(yuǎn)小于固相導(dǎo)熱與流動(dòng)的特征時(shí)間，也就是認(rèn)為氣相反應(yīng)是準(zhǔn)穩(wěn)態(tài)的，反應(yīng)放出的熱量為Q。假設(shè)氣相反應(yīng)的劉易斯數(shù)為1，則高溫燃?xì)鈪^(qū)的溫度梯度為零，且無化學(xué)反應(yīng)。固相反應(yīng)區(qū)對(duì)應(yīng)于穩(wěn)態(tài)燃燒模型中的嘶嘶區(qū)與亞表面反應(yīng)區(qū)，包含了推進(jìn)劑的融化、氣化、熱解、固相表面反應(yīng)、氣泡中的氣相反應(yīng)等復(fù)雜的物化過程，整個(gè)過程放出的熱量為Qs。

圖2 一維準(zhǔn)靜態(tài)燃燒模型

推進(jìn)劑燃速由阿倫尼烏斯公式?jīng)Q定：

(6)

對(duì)氣相反應(yīng)區(qū)建立能量守恒方程，有：

(7)

對(duì)固相反應(yīng)區(qū)建立能量守恒方程，有：

(8)

氣相對(duì)固相的導(dǎo)熱量[6-7]為：

(9)

聯(lián)立式(7)、式(8)、式(9)，消去Ts與Tf，得：

(10)

1.2 數(shù)值方法

1.2.1 氣相數(shù)值方法

氣相歐拉方程采用有限體積方法進(jìn)行計(jì)算。對(duì)于方程(1)的半離散格式可以寫為：

(11)

空間重構(gòu)使用一階MUSCL格式，數(shù)值通量的計(jì)算采用Roe格式。網(wǎng)格數(shù)經(jīng)過網(wǎng)格無關(guān)性分析，取1 600。

1.2.2 固相數(shù)值方法

固相導(dǎo)熱控制方程采用有限差分的方法進(jìn)行離散。網(wǎng)格數(shù)取100。

(12)

其邊界條件為：

(13)

方程(12)的半離散形式為：

(14)

1.2.3 氣固耦合與時(shí)間推進(jìn)方法

固相導(dǎo)熱的邊界溫度Ts由式(6)與(10)確定。邊界上，單側(cè)的有限差分格式為：

(15)

(16)

時(shí)間推進(jìn)使用二階Runge-Kutta方法。

1.3 理論預(yù)估方法

根據(jù)發(fā)動(dòng)機(jī)穩(wěn)定性線性理論[3]，發(fā)動(dòng)機(jī)聲振基頻模態(tài)聲壓振幅的線性增長(zhǎng)率為：

α=αb+αN

(17)

其中αb為燃面增益系數(shù)，αN為噴管阻尼系數(shù)。其具體的表達(dá)式為[8]：

(18)

燃面聲導(dǎo)納函數(shù)[3]為：

(19)

壓力耦合的燃燒響應(yīng)函數(shù)[4,8]可以寫為：

(20)

其中：

對(duì)于短噴管，噴管聲導(dǎo)納函數(shù)為：

(21)

2 計(jì)算結(jié)果與分析

2.1 發(fā)動(dòng)機(jī)試車工況

某工程單位改性雙基裝藥小推力固體火箭發(fā)動(dòng)機(jī)試車中出現(xiàn)壓力振蕩現(xiàn)象。發(fā)動(dòng)機(jī)裝藥與內(nèi)部構(gòu)型如圖3所示。

圖3 發(fā)動(dòng)機(jī)裝藥與內(nèi)部結(jié)構(gòu)示意圖

發(fā)動(dòng)機(jī)在低溫點(diǎn)火時(shí)，實(shí)驗(yàn)結(jié)果顯示，在一、二級(jí)過渡段位置出現(xiàn)推力振蕩，與此同時(shí)，壓強(qiáng)曲線也存在類似的振蕩。分析認(rèn)為，低溫點(diǎn)火時(shí)工作壓強(qiáng)過低，誘發(fā)了不穩(wěn)定燃燒現(xiàn)象。具體的曲線如圖4、圖5與圖6、圖7所示。根據(jù)推力頻譜曲線(如圖8)，可知振蕩頻率主要分布于1 324 Hz、2 500～2 800 Hz、4 000Hz等三個(gè)位置；根據(jù)壓強(qiáng)頻譜曲線(如圖9)，可知振蕩主要位于2 500～2 800 Hz范圍，在1 324 Hz和4 000 Hz附近的振蕩幅度較小。

圖4 低溫點(diǎn)火時(shí)的推力曲線

圖5 低溫點(diǎn)火時(shí)振蕩段的推力曲線

圖6 低溫點(diǎn)火時(shí)的壓強(qiáng)曲線

圖7 低溫點(diǎn)火時(shí)振蕩段的壓力曲線

圖8 推力振蕩段的頻譜分析

圖9 壓強(qiáng)振蕩段的頻譜分析

根據(jù)發(fā)動(dòng)機(jī)幾何與振動(dòng)頻率，進(jìn)行簡(jiǎn)單的模態(tài)分析后認(rèn)為，發(fā)動(dòng)機(jī)出現(xiàn)的燃燒燃燒不穩(wěn)定為縱向聲振不穩(wěn)定燃燒。

2.2 發(fā)動(dòng)機(jī)典型工況計(jì)算結(jié)果

取發(fā)動(dòng)機(jī)常溫工況為例，計(jì)算參數(shù)如表1所示。

表1 發(fā)動(dòng)機(jī)常溫計(jì)算參數(shù)

先求解穩(wěn)態(tài)解，不考慮固相導(dǎo)熱的影響。將穩(wěn)態(tài)解作為非穩(wěn)態(tài)求解的初始條件，并在初始時(shí)刻對(duì)燃面端部施加一個(gè)1%室壓的壓力脈沖，進(jìn)行非穩(wěn)態(tài)計(jì)算。

2.2.1 振型與頻譜分析

發(fā)動(dòng)機(jī)聲腔對(duì)于不同頻率的聲振蕩有調(diào)制作用。在計(jì)算29 ms后，發(fā)動(dòng)機(jī)聲振振幅的空間分布如圖10所示。

取發(fā)動(dòng)機(jī)燃面、距燃面168 mm與287.1 mm三處壓力為監(jiān)測(cè)點(diǎn)。其中距離燃面168 mm位置在后封頭處。287.1 mm位置為長(zhǎng)尾管的中部，為聲振駐波波節(jié)。三處聲壓隨時(shí)間變化如圖11、圖12和圖13所示。

圖10 聲壓振幅與速度脈動(dòng)空間分布(29 ms)

圖11 燃面聲壓隨時(shí)間變化

圖12 后封頭聲壓隨時(shí)間變化

圖13 長(zhǎng)尾管波節(jié)處聲壓隨時(shí)間變化

對(duì)上述三個(gè)位置20～30 ms時(shí)間段的聲壓做頻譜分析，結(jié)果如圖14所示。

圖14 聲壓振蕩頻譜分析

三處聲壓峰值頻率為2 988 Hz。根據(jù)發(fā)動(dòng)機(jī)穩(wěn)態(tài)解，知道發(fā)動(dòng)機(jī)燃燒室平均音速為a=983.4 m/s。燃燒室長(zhǎng)度為L(zhǎng)=167 mm。根據(jù)f=a/2L，可以估算發(fā)動(dòng)機(jī)縱向聲振頻率為f=2 944.2 Hz。與數(shù)值計(jì)算值相當(dāng)。根據(jù)聲學(xué)原理知道，速度梯度可以有效的反射聲能。對(duì)于長(zhǎng)尾管發(fā)動(dòng)機(jī)，存在兩處速度梯度較大的區(qū)域，一處是后封頭收斂段，另一處是噴管收斂段。聲能在流場(chǎng)中傳播，在后封頭處大部分聲能被反射，在燃燒室內(nèi)形成一個(gè)2 900 Hz左右的駐波。計(jì)算所得壓力振蕩頻率與實(shí)驗(yàn)試車中壓力振蕩頻率以及理論估計(jì)頻率都符合較好。

2.2.2 發(fā)動(dòng)機(jī)聲振線性衰減率

取20～30 ms燃面處的聲壓振幅，從而得到聲壓振蕩振幅隨時(shí)間變化關(guān)系，在對(duì)數(shù)坐標(biāo)系下做出下圖(見圖15)：

圖15 燃面處聲壓振幅隨時(shí)間變化關(guān)系

根據(jù)式(17)～(21)計(jì)算所得聲壓振幅理論衰減率為α=-23.1。理論計(jì)算與數(shù)值計(jì)算所得的衰減率符合較好。

3.3 發(fā)動(dòng)機(jī)低溫不穩(wěn)定原因分析

考慮發(fā)動(dòng)機(jī)初始溫度為高(323 K)、低(233 K)、常溫(293 K)的3種工況。除初始溫度外，其他參數(shù)由表1給出。發(fā)動(dòng)機(jī)燃面處聲壓隨時(shí)間變化關(guān)系如圖16和圖17。

圖16 高溫工況燃面聲壓隨時(shí)間變化關(guān)系

圖17 低溫工況燃面聲壓隨時(shí)間變化關(guān)系

不同溫度下，峰值頻率都為2 941.2 Hz。區(qū)別只在于聲壓振幅變化率不同。聲壓振幅的變化率如表2。

表2 不同初溫聲壓變化率

發(fā)動(dòng)機(jī)在低溫下是不穩(wěn)定的，在常溫與高溫下是穩(wěn)定的。根據(jù)式(20)做出三種工況下的燃燒響應(yīng)函數(shù)如圖18所示。

圖18 燃燒響應(yīng)函數(shù)實(shí)部隨頻率變化關(guān)系

圖18中豎直虛線為發(fā)動(dòng)機(jī)聲振頻率?？梢钥闯?，隨著初溫的提高，燃燒相應(yīng)函數(shù)實(shí)部峰值頻率右移，逐漸偏離了發(fā)動(dòng)機(jī)聲腔的固有頻率。低溫工況下壓力耦合的響應(yīng)函數(shù)峰值頻率正好處在發(fā)動(dòng)機(jī)燃燒室的固有頻率下，所以發(fā)動(dòng)機(jī)在低溫下出現(xiàn)了不穩(wěn)定燃燒的現(xiàn)象。

2.4 發(fā)動(dòng)機(jī)低溫不穩(wěn)定抑制方法

對(duì)燃速系數(shù)進(jìn)行調(diào)節(jié)，從而使得穩(wěn)態(tài)燃速不同。采用穩(wěn)態(tài)燃速公式，折算到1.4 MPa壓力、293 K的工況條件下，分別令折算燃速為10、11、12 mm/s。在低溫(233 K)工況下，進(jìn)行穩(wěn)態(tài)計(jì)算，再采用穩(wěn)態(tài)計(jì)算結(jié)果作為初始條件，進(jìn)行非穩(wěn)態(tài)計(jì)算。其他參數(shù)與低溫工況相同。表3為不同燃速下低溫工況聲壓振幅變化率：

表3 不同燃速下低溫工況聲壓變化率

從表3可以看出隨燃速增加，低溫下發(fā)動(dòng)機(jī)工作逐漸穩(wěn)定。這是因?yàn)樘岣呷妓?，燃燒響?yīng)函數(shù)的峰值頻率右移，從而偏離了發(fā)動(dòng)機(jī)聲腔的固有頻率。

實(shí)際工程中，通過增加發(fā)動(dòng)機(jī)二級(jí)燃速、提高了燃燒穩(wěn)定劑的比例、調(diào)整燃燒穩(wěn)定劑的粒徑，成功抑制了燃燒不穩(wěn)定現(xiàn)象的發(fā)生。

3 結(jié)論

1) 基于發(fā)動(dòng)機(jī)一維模型與一維準(zhǔn)靜態(tài)燃燒模型，發(fā)展了一種數(shù)學(xué)模型和數(shù)值方法用于發(fā)動(dòng)機(jī)線性穩(wěn)定性預(yù)估。計(jì)算結(jié)果與線性理論符合較好。

2) 針對(duì)某工程單位某型改性雙基裝藥長(zhǎng)尾管固體火箭發(fā)動(dòng)機(jī)燃燒不穩(wěn)定現(xiàn)象開展了數(shù)值模擬研究。認(rèn)為低溫工況下，發(fā)動(dòng)機(jī)燃燒響應(yīng)函數(shù)實(shí)部在發(fā)動(dòng)機(jī)固有頻率下較大，容易誘發(fā)不穩(wěn)定燃燒。計(jì)算分析表明，提高燃速可以有效抑制不穩(wěn)定燃燒。

3) 本文中的數(shù)學(xué)模型沒有考慮多維效應(yīng)、兩相流動(dòng)、與氣體粘性。在將來的研究工作中，可以將本文中模型進(jìn)行擴(kuò)展，考慮兩相流動(dòng)、粘性等的影響?？梢宰龅礁鼮槿娴貙?duì)發(fā)動(dòng)機(jī)整機(jī)的穩(wěn)定性進(jìn)行預(yù)估。

[1] CULICK F E C.Combustion Instability in Solid Rocket Motors.Volume II.A Guide for Motor Designers[J].Combustion Instability in Solid Rocket Motors, 1980.

[2] 劉佩進(jìn), 何國(guó)強(qiáng).固體火箭發(fā)動(dòng)機(jī)燃燒不穩(wěn)定及控制技術(shù)[M].西安：西北工業(yè)大學(xué)出版社, 2015.

[3] 孫維申.固體火箭發(fā)動(dòng)機(jī)不穩(wěn)定燃燒[M].北京:北京工業(yè)學(xué)院出版社, 1988.

[4] SUMMERFIELD M, PRICE E W, LUCA L D.Nonsteady Burning and Combustion Stability of Solid Propellants[M].1971：601-641.

[5] ZEL’DOVICH Y B.Theory of propellant combustion in a gas flow[J].Combustion Explosion & Shock Waves, 1971, 7(4):399-408.

[6] HANZAWA M.A theoretical study on depressurization induced extinction of solid propellant[C]// Propulsion Conference.2006.

[7] SHIMADA T, HANZAWA M, MORITA T, et al.Stability Analysis of Solid Rocket Motor Combustion by Computational Fluid Dynamics[J].Aiaa Journal, 2015, 46(46):947-957.

[8] CULICK F E C.A review of calculations for unsteady burning of a solid propellant.[J].Aiaa Journal, 2012, 6(12):2241-2255.

StabilityAnalysisofSolidRochetMotorCombustionwithDoubleBasePropellantbyComputationalFluidDynamics

CHEN Junyi1，WANG Bing1，TIAN Xiaotao2

(1.Tsinghua University, Beijing 100084, China;2.Xi’an Modern Control Technology Research Institute, Xi’an 710065, China)

The acoustic combustion instability of a double base propellant rocket motor is investigated by computational fluid dynamics.The quasi one dimensional Euler equation and quasi steady flame model was used to simulate the unsteady flow.Calculation results show that, at low temperature, the peak frequency of pressure coupled response function is close to the natural frequency of combustion chamber, which induced unsteady combustion.By increasing the burning rate, unstable combustion can be effectively suppressed.The numerical results agree well with the linear theory.

solid rocket motor; instability of combustion; numerical simulation; modified double base propellant.

2017-10-05；

2017-10-29

陳俊屹(1993—)，男，碩士研究生，主要從事固發(fā)不穩(wěn)定燃燒研究。

王兵(1977—)，男，博士，副教授，博士生導(dǎo)師，主要從事航空宇航領(lǐng)域的前沿基礎(chǔ)科學(xué)問題、燃燒不穩(wěn)定性、先進(jìn)數(shù)值模擬方法及工程應(yīng)用研究。

10.11809/scbgxb2017.12.044

本文引用格式：陳俊屹,王兵,田小濤.某型固體火箭發(fā)動(dòng)機(jī)燃燒不穩(wěn)定性數(shù)值模擬研究[J].兵器裝備工程學(xué)報(bào)，2017(12):195-200,252.

formatCHEN Junyi，WANG Bing，TIAN Xiaotao.Stability Analysis of Solid Rochet Motor Combustion with Double Base Propellant by Computational Fluid Dynamics[J].Journal of Ordnance Equipment Engineering,2017(12):195-200,252.

TJ7

A

2096-2304(2017)12-0195-06

(責(zé)任編輯楊繼森)