許團(tuán)委，田維平，王建儒

(1西北工業(yè)大學(xué)航天學(xué)院，西安 710072;2中國(guó)航天科技集團(tuán)公司第四研究院第41所，西安 710025)

0 引言

固體發(fā)動(dòng)機(jī)不穩(wěn)定燃燒現(xiàn)象的主要特征是燃燒室壓強(qiáng)、燃速以及推力等參數(shù)做周期或近似周期性的變化。諸多因素可誘發(fā)不穩(wěn)定燃燒，如推進(jìn)劑壓強(qiáng)耦合響應(yīng)、速度耦合響應(yīng)、分布燃燒等[1]。近年來國(guó)外針對(duì)多個(gè)發(fā)動(dòng)機(jī)開展了不穩(wěn)定性的影響分析研究[2－3]，以法國(guó)為主的歐洲啟動(dòng)兩大計(jì)劃研究壓強(qiáng)振蕩問題，分別是:分段式固體發(fā)動(dòng)機(jī)氣體動(dòng)力學(xué)(ASSM)和壓強(qiáng)振蕩計(jì)劃(POP)。ASSM計(jì)劃的主要目標(biāo)是對(duì)渦脫落深入理解和建模，促進(jìn)數(shù)值模擬技術(shù)的發(fā)展。POP計(jì)劃是利用P230的縮比模型發(fā)動(dòng)機(jī)開展實(shí)驗(yàn)研究，獲得實(shí)驗(yàn)和數(shù)值數(shù)據(jù)庫。美國(guó)在此也投入了大量的人力和經(jīng)費(fèi)，開展了多學(xué)科大學(xué)研究倡議(MURI)，試圖從基礎(chǔ)化學(xué)、燃燒和流體動(dòng)力學(xué)的角度深入研究燃燒不穩(wěn)定課題[4]。但大多停留在已有結(jié)構(gòu)的驗(yàn)證分析方面，對(duì)于具體結(jié)構(gòu)設(shè)計(jì)參數(shù)對(duì)流動(dòng)穩(wěn)定性的影響未見報(bào)道。隨著我國(guó)針對(duì)分段式固體發(fā)動(dòng)機(jī)技術(shù)研究的不斷深入，開展單因素對(duì)壓強(qiáng)振蕩的影響研究至關(guān)重要。文中運(yùn)用大渦模擬方法，進(jìn)行了絕熱環(huán)及潛入式噴管空腔背壁區(qū)對(duì)壓強(qiáng)振蕩的影響分析，獲得了關(guān)鍵設(shè)計(jì)參數(shù)對(duì)分段式發(fā)動(dòng)機(jī)壓強(qiáng)振蕩特性的一般影響規(guī)律，對(duì)于工程應(yīng)用及后續(xù)的多因素耦合分析有一定參考意義。

1 計(jì)算模型與數(shù)值算法

1)控制方程

針對(duì)氣動(dòng)聲學(xué)、燃燒室中湍流動(dòng)量輸運(yùn)及分離流或渦流區(qū)流動(dòng)等問題，LES以其耗散小、精度高等特點(diǎn)比RANS具有更大的優(yōu)勢(shì)。因此，可利用LES來直接計(jì)算壓強(qiáng)振蕩引起的聲波運(yùn)動(dòng)。文中的研究涉及的正是分離流、旋渦脫落、壓強(qiáng)振蕩及聲學(xué)。LES采用過濾器對(duì)N-S方程進(jìn)行過濾，將小尺度脈動(dòng)過濾掉，得到亞格子應(yīng)力項(xiàng)。然后，直接計(jì)算大尺度脈動(dòng)。文中利用空間濾波器G將流場(chǎng)變量分為可解尺度脈動(dòng)量與亞格子尺度脈動(dòng)量:

考慮到氣體的可壓縮性，利用Favre平均對(duì)控制方程進(jìn)行簡(jiǎn)化:

式中:“～”表示Favre平均，“—”表示Reynold平均。

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

式中:μ為動(dòng)力粘性系數(shù)，Cp為定壓比熱。將氣體工質(zhì)做理想氣體處理，物性參數(shù)如表1所示。

表1 物性參數(shù)表

2)亞格子模型選擇及控制方程離散

選擇WALE模型對(duì)亞格子應(yīng)力張量τsgsij、亞格子通量張量Hsgsi以及亞格子尺度粘性力變形功Θsgsi進(jìn)行封閉。

對(duì)于連續(xù)方程與動(dòng)量方程，采用BCD(bounded central differencing)格式進(jìn)行離散，該格式可以發(fā)現(xiàn)求解區(qū)域中波長(zhǎng)小于2△X的波動(dòng)，并加以抑制;能量方程則采用Power Law格式以加速收斂。

2 研究對(duì)象的幾何模型及邊界條件

鑒于國(guó)外對(duì)于VKI冷流模型，從實(shí)驗(yàn)和數(shù)值模擬方面已經(jīng)做了大量工作，并公開了模型的幾何尺寸、實(shí)驗(yàn)和計(jì)算結(jié)果。文中亦選用VKI冷流模型發(fā)動(dòng)機(jī)進(jìn)行數(shù)值計(jì)算，計(jì)算區(qū)域如圖1所示，含有絕熱環(huán)及潛入式噴管。

圖1 模型發(fā)動(dòng)機(jī)計(jì)算區(qū)域網(wǎng)格劃分

旋渦具有強(qiáng)三維特性，雖然軸對(duì)稱模型不能很好的模擬旋渦拉伸現(xiàn)象，但三維模擬計(jì)算量太大。文獻(xiàn)[5]使用二維大渦模擬成功的預(yù)估了 ETM-03與RSRM的振蕩特性，因此文中仍采用二維軸對(duì)稱模型。為了更好的識(shí)別邊界層區(qū)域，對(duì)網(wǎng)格在徑向進(jìn)行加密，其余地方確保網(wǎng)格尺度在1mm之內(nèi)，對(duì)于各種工況，計(jì)算網(wǎng)格數(shù)量約為11萬。

根據(jù)文獻(xiàn)[6]，發(fā)動(dòng)機(jī)采用軸向進(jìn)氣，將冷氣簡(jiǎn)化為理想氣體。入口靜壓Ps=0.178MPa(與試驗(yàn)一致)，出口設(shè)定為壓強(qiáng)出口，出口截面參數(shù)可通過內(nèi)部外推得到，在壁面上選取無滑移邊界條件。

3 結(jié)構(gòu)特性對(duì)壓強(qiáng)振蕩的影響分析

四種設(shè)計(jì)條件下的發(fā)動(dòng)機(jī)結(jié)構(gòu)二維軸對(duì)稱后半部幾何結(jié)構(gòu)如圖2所示，其網(wǎng)格劃分均在熱隔板處及轉(zhuǎn)角處進(jìn)行網(wǎng)格加密。其中Mode3為基準(zhǔn)模型，在此基礎(chǔ)上分別引入絕熱環(huán)(Mode4)及潛入噴管(Mode2)或者二者都引入，如表2所示。為了便于比較不同結(jié)構(gòu)對(duì)燃燒室壓強(qiáng)振蕩的影響，對(duì)燃燒室內(nèi)的不同位置進(jìn)行了監(jiān)測(cè)點(diǎn)的布置，以獲得壓強(qiáng)－時(shí)間曲線，并對(duì)其進(jìn)行FFT分析。

表2 四種不同發(fā)動(dòng)機(jī)的結(jié)構(gòu)特征

3.1 穩(wěn)態(tài)流場(chǎng)尾部流動(dòng)分析

在進(jìn)行模擬計(jì)算中，首先開展的是穩(wěn)態(tài)流場(chǎng)計(jì)算。經(jīng)計(jì)算，發(fā)動(dòng)機(jī)穩(wěn)態(tài)流場(chǎng)尾部流線圖如圖2所示。1)絕熱環(huán)的存在對(duì)流場(chǎng)的影響，Mode1和Mode4在絕熱環(huán)后端形成了回流旋渦;2)潛入式噴管對(duì)流場(chǎng)的影響，Mode1與Mode2的潛入空腔內(nèi)也形成了回流旋渦，而Mode3穩(wěn)態(tài)流場(chǎng)中沒有形成回流旋渦。初步表明，絕熱環(huán)是引起障礙渦脫落的主要因素，絕熱環(huán)后端及潛入空腔是旋渦運(yùn)動(dòng)的主要區(qū)域。

圖2 四種不同結(jié)構(gòu)發(fā)動(dòng)機(jī)穩(wěn)態(tài)流場(chǎng)尾部流線圖

3.2 壓強(qiáng)振蕩的頻率與振幅特性分析

對(duì)四種模型各自頭部監(jiān)測(cè)點(diǎn)所獲得的壓強(qiáng)隨時(shí)間變化曲線進(jìn)行FFT分析，以便開展壓強(qiáng)振蕩的頻率與振幅特性分析，Mode2～Mode4的計(jì)算結(jié)果如圖3～圖5所示。

圖3 Mode2發(fā)動(dòng)機(jī)頭部壓強(qiáng)隨時(shí)間變化曲線及FFT分析結(jié)果

圖4 Mode3發(fā)動(dòng)機(jī)頭部壓強(qiáng)隨時(shí)間變化曲線及FFT分析結(jié)果

圖5 Mode4發(fā)動(dòng)機(jī)頭部壓強(qiáng)隨時(shí)間變化曲線及FFT分析結(jié)果

1)頻率特性分析

從圖6可以看出，對(duì)于四種幾何結(jié)構(gòu)局部不同的發(fā)動(dòng)機(jī)前兩階壓強(qiáng)振蕩頻率整體上有逐漸增大的趨勢(shì)。對(duì) Mode1和Mode2進(jìn)行對(duì)比分析，認(rèn)為絕熱環(huán)的存在，使得前兩階頻率有所降低，主要是各自產(chǎn)生旋渦脫落的部位不同，前者絕熱環(huán)引起的障礙旋渦脫落是主因，后者潛入空腔入口引起的轉(zhuǎn)角渦脫落是主因;對(duì)Mode1和Mode4進(jìn)行對(duì)比分析，認(rèn)為潛入式空腔的增加相當(dāng)于延長(zhǎng)了聲腔軸向長(zhǎng)度，依據(jù)軸向固有頻率預(yù)估公式前兩階頻率有所降低。

2)振幅特性分析

從旋渦強(qiáng)度分布云圖7以及振幅比分析圖8、圖9，可以看出Mode1～Mode4中均出現(xiàn)不同強(qiáng)度的渦脫落和不同程度的壓強(qiáng)振蕩。下面分兩類對(duì)結(jié)果進(jìn)行分析:

圖6 四種不同結(jié)構(gòu)發(fā)動(dòng)機(jī)前兩階壓強(qiáng)振蕩頻率對(duì)比

圖7 四種不同結(jié)構(gòu)發(fā)動(dòng)機(jī)某時(shí)刻旋渦強(qiáng)度分布云圖

圖8 絕熱環(huán)對(duì)振幅比的影響

①絕熱環(huán)的影響

從圖8可以看出，振幅比大小順序?yàn)镸ode1＞Mode2＞Mode3。 分 析認(rèn)為主要有兩方面的原因:1)旋渦的產(chǎn)生源不同，Mode1中旋渦脫落來自突出在流場(chǎng)中的絕熱環(huán)引起的障礙渦脫落和由于內(nèi)流不穩(wěn)定而產(chǎn)生在裝藥表面的脫落渦，而Mode2中的旋渦主要是后者——表面渦脫落;2)Mode1中絕熱環(huán)形成的障礙渦脫落頻率可能與一階聲模態(tài)非常接近，從而誘發(fā)了較大程度的共振，而Mode2誘發(fā)的共振程度可能較小。因此初步認(rèn)為，障礙渦脫落對(duì)流場(chǎng)產(chǎn)生的不穩(wěn)定程度高于表面渦脫落，與國(guó)外學(xué)者所提出的一致。

②噴管形式的影響

從圖9可以看出，振幅比大小順序?yàn)镸ode1＞Mode4＞Mode3。分析認(rèn)為Mode1中引入了潛入噴管，所形成的潛入空腔類似于共振器，當(dāng)脫落的旋渦穿越空腔入口處的聲速線時(shí)將會(huì)增大聲功率，進(jìn)而增大壓強(qiáng)振幅，由于上游具有穩(wěn)定的障礙渦脫落源，從而空腔可以一直維持聲渦耦合誘發(fā)的聲能增益。

其中，Mode3中壓強(qiáng)振幅比非常小(最大為0.096%)，對(duì)發(fā)動(dòng)機(jī)工作穩(wěn)定性的影響可忽略不計(jì)。

圖9 噴管形式對(duì)振幅比的影響

4 結(jié)論

1)在分段式固體發(fā)動(dòng)機(jī)中，絕熱環(huán)是引起燃燒室旋渦脫落的主要因素;

2)潛入式噴管空腔背壁區(qū)維持并加強(qiáng)了燃燒室聲渦耦合引起的聲能增益，對(duì)壓強(qiáng)振蕩有較大的影響;

3)多種渦脫落共同存在的燃燒室中，由障礙渦脫落產(chǎn)生的流動(dòng)不穩(wěn)定程度大于表面渦脫落。

[1]陳曉龍，何國(guó)強(qiáng)，劉佩進(jìn)，等.固體火箭發(fā)動(dòng)機(jī)燃燒不穩(wěn)定的影響因素分析和最新研究進(jìn)展[J].固體火箭技術(shù)，2009，32(6):600 －605.

[2]Fred S Blomshield.Lessons learned in solid rocket combustion instability[C]//43rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference ＆ Exhibit，AIAA 2007 －5803，July 2007.

[3]Stany GALLIER，Michel PREVOST，Jouke HIJLKEMA，Michaёl ROUMY .Effects of cavity on thrust oscillations in subscale solid rocket motors[C]//45st AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit，AIAA－2009－5253.

[4]Fred S Blomshied.Summary of multi-disciplinary university research initiative in solid propellant combustion instability，AIAA 2000 －3172[R].2000.

[5]Mason D R，Morstadt R A，Cannon S M.Pressure oscillation and structural vibrations in space shuttle RSRM and ETM-3 motors，AIAA 2004 －3898[R].2004.

[6]Stella F，Paglia F.Pressure oscillations in solid rocket motors:Numerical study[J].Aerospace Science and Technology，2011，15(1):53－59.