李映坤，韓珺禮，2，陳 雄，鞠玉濤

(1.南京理工大學(xué) 機(jī)械工程學(xué)院，南京 210094;2.北京機(jī)電研究所，北京 100012)

0 引言

固體火箭發(fā)動(dòng)機(jī)具有結(jié)構(gòu)簡單、使用方便、可靠性好等優(yōu)點(diǎn)，在國防領(lǐng)域受到廣泛重視及應(yīng)用。與傳統(tǒng)的單室單推、單室雙推發(fā)動(dòng)機(jī)相比，雙脈沖固體火箭發(fā)動(dòng)機(jī)具有改善固體火箭發(fā)動(dòng)機(jī)能量可控性的突出優(yōu)勢，可多次點(diǎn)火，并提供不連續(xù)推力，有利于提高導(dǎo)彈武器的作戰(zhàn)能力［1］。

級間隔離裝置是雙脈沖固體火箭發(fā)動(dòng)機(jī)的關(guān)鍵部件之一。根據(jù)隔離裝置承力情況，隔離裝置可分為隔板式(硬隔離)和隔層式(軟隔離)。其中，隔板式隔離裝置又包括了隔塞式隔板、陶瓷艙蓋隔板、金屬膜片式隔板等，它們都具有結(jié)構(gòu)簡單、易于加工、承壓以及密閉可靠性好等優(yōu)點(diǎn)，國內(nèi)外對類似結(jié)構(gòu)的雙脈沖發(fā)動(dòng)機(jī)已進(jìn)行了試驗(yàn)研究［2－3］，但試驗(yàn)測得的推力-時(shí)間曲線與計(jì)算結(jié)果有一定的偏差，難以重現(xiàn)內(nèi)彈道的性能。這是因?yàn)樵诙壝}沖燃燒室工作期間，一級脈沖燃燒室中存在多個(gè)隔塞(或者陶瓷碎片)，隔塞從噴管排出的過程中對噴管的內(nèi)流場產(chǎn)生了影響，造成發(fā)動(dòng)機(jī)瞬時(shí)推力的波動(dòng)，而國內(nèi)外未見對隔塞在噴管中運(yùn)動(dòng)過程仿真的報(bào)道。

本文基于AUSM-PW矢通量分裂格式，求解圓柱坐標(biāo)系下的二維軸對稱非定?？蓧嚎s雷諾平均Navier-Stokes方程，采用基于格心的有限體積法及雙時(shí)間步LU-SGS方法，結(jié)合運(yùn)動(dòng)嵌套網(wǎng)格技術(shù)，研究了隔塞沿著噴管軸線運(yùn)動(dòng)過程中噴管內(nèi)流場的演變過程，并分析了隔塞在噴管內(nèi)運(yùn)動(dòng)對發(fā)動(dòng)機(jī)推力的影響，為相關(guān)設(shè)計(jì)提供了參考。

1 控制方程及數(shù)值求解方法

1.1 流動(dòng)控制方程

采用ALE(Arbitrary Lagrangian Eulerian)有限體積方法描述的可壓縮非定常Navier-Stokes方程為

式中 U為守恒變量的一般形式;Fc為無粘通量;Fv為粘性通量;H為軸對稱幾何源項(xiàng);Fc=Fi+Gj;Fv=Fvi+Gvj;(i，j)表示直角坐標(biāo)系2個(gè)坐標(biāo)方向(x，y)的單位矢量。

式中 ρ為氣體密度;u、v為氣體運(yùn)動(dòng)速度矢量的2個(gè)分量;xt、yt分別為網(wǎng)格運(yùn)動(dòng)速度矢量的2個(gè)分量;T為氣體的溫度;E為單位體積氣體的總能量;τ表示應(yīng)力張量，其具體形式參考文獻(xiàn)［4］。

對于湍流問題，本文采用Menter F R提出的k-ε剪切應(yīng)力輸運(yùn)(shear-stress-transport)模式［5］，該模型通過混合函數(shù)F1將k-ε模型和k-ε模型結(jié)合起來，這樣充分發(fā)揮了k-ε模型對自由流和k-ω模型對壁面受限流動(dòng)的處理優(yōu)勢。具體描述如下:

式中 k為湍動(dòng)能;ω為比耗散率;μt為湍流粘性系數(shù)，其他它參數(shù)的具體形式見參考文獻(xiàn)［5］。

1.2 數(shù)值求解方法

本文采用基于格心的有限體積法，對方程進(jìn)行離散。其中，無粘項(xiàng)采用AUSM-PW失通量分裂格式。以i方向上的通量F為例，i+1/2邊界上的通量可寫為［6］

式中 c為單元界面聲速;Φ為守恒通量;p為壓力項(xiàng)。

粘性項(xiàng)采用Jameson中心差分法離散，為了提高非定常流動(dòng)的時(shí)間計(jì)算精度，同時(shí)又要具有較高的計(jì)算效率，本文采用Jameson提出的一種雙時(shí)間步［7］的計(jì)算方法，在凍結(jié)的真實(shí)時(shí)刻點(diǎn)上，引入類似牛頓迭代的虛擬時(shí)間迭代過程，通過這種內(nèi)迭代過程來提高LU-SGS隱式算法處理過程中損失的時(shí)間精度。

2 嵌套網(wǎng)格技術(shù)

隔塞在噴管內(nèi)運(yùn)動(dòng)的過程為典型的移動(dòng)邊界問題。目前，處理移動(dòng)邊界問題較成熟的方法包括非結(jié)構(gòu)動(dòng)網(wǎng)格技術(shù)及動(dòng)態(tài)嵌套網(wǎng)格技術(shù)等，但非結(jié)構(gòu)網(wǎng)格在處理具有較大相對位移問題時(shí)，網(wǎng)格變形與重構(gòu)的表現(xiàn)能力較差。因此，本文選擇動(dòng)態(tài)嵌套網(wǎng)格方法，其優(yōu)點(diǎn)是對于復(fù)雜幾何體，不要求流場各個(gè)計(jì)算域共享邊界，減輕了網(wǎng)格生成難度，適用于運(yùn)動(dòng)幅度較大的相對運(yùn)動(dòng)。

2.1 洞邊界的確定

在嵌套網(wǎng)格系統(tǒng)中，建立人工插值內(nèi)邊界的過程稱之為挖洞過程。為了提高確定洞邊界的效率，本文采用 Chiu I T和 Meakin R L提出的“Hole-Map”方法［8］，其核心思想是對于給定的嵌套網(wǎng)格體系，若已知其拓?fù)浣Y(jié)構(gòu)，就能用均勻的笛卡爾網(wǎng)格單元去近似每個(gè)曲面邊界，從而得到該曲面邊界的笛卡爾網(wǎng)格近似。在此基礎(chǔ)上，網(wǎng)格點(diǎn)P與曲面邊界的位置關(guān)系可近似轉(zhuǎn)換為網(wǎng)格點(diǎn)P與笛卡爾網(wǎng)格的關(guān)系，而這種關(guān)系十分容易確定。

2.2 貢獻(xiàn)單元的搜尋

在嵌套網(wǎng)格系統(tǒng)中，通過貢獻(xiàn)單元將流場解的信息插值到插值邊界面網(wǎng)格上。因此，尋找貢獻(xiàn)單元技術(shù)的優(yōu)劣，對整個(gè)嵌套網(wǎng)格方法有很大影響，有時(shí)甚至是嵌套網(wǎng)格方法成敗的關(guān)鍵。本文結(jié)合Wang Z J的矢量判別法［9］，提出了一種簡單易行的尋找插值點(diǎn)貢獻(xiàn)單元的方法，具體步驟如下:

(1)對于任意給定網(wǎng)格單元，尋找距離其最近的網(wǎng)格單元;

(2)采用Wang Z J的矢量判別法，尋找包圍該點(diǎn)的貢獻(xiàn)單元。對于任意網(wǎng)格單元，如果其所有邊界均滿足rfc·n≥0，則點(diǎn)C位于網(wǎng)格單元內(nèi)，從而該網(wǎng)格單元就為點(diǎn)C的貢獻(xiàn)單元，其中f為網(wǎng)格單元邊界的中心，n為邊界外法線矢量;

(3)使用雙線性插值方法進(jìn)行插值，完成該點(diǎn)物理信息的傳遞。

3 算例驗(yàn)證

3.1 靜態(tài)NACA0012翼型繞流

為了驗(yàn)證本文靜態(tài)嵌套網(wǎng)格計(jì)算程序，根據(jù)文獻(xiàn)［10］，計(jì)算了攻角 α0=2.05°，來流馬赫數(shù) Ma∞=0.755，雷諾數(shù)為 9.9×106工況下 NACA0012 翼型繞流。為了便于對比分析，本文分別采用了單塊網(wǎng)格和嵌套網(wǎng)格。圖1所示的是采用單塊網(wǎng)格和嵌套網(wǎng)格計(jì)算的翼型表面壓力系數(shù)分布，并與實(shí)驗(yàn)結(jié)果進(jìn)行了對比。從圖1中可看出，本文的計(jì)算與實(shí)驗(yàn)值吻合較好，而且嵌套網(wǎng)格計(jì)算的結(jié)果與單塊網(wǎng)格計(jì)算的結(jié)果相差不大。因此，證明本文所編制的計(jì)算程序可信、可靠。

圖1 翼型表面壓力系數(shù)對比Fig.1 Comparison of steady pressuredistributions for the airfoil

3.2 俯仰振動(dòng)的NACA0012翼型繞流

以NACA0012翼型繞1/4弦點(diǎn)作簡諧俯仰振動(dòng)為算例，驗(yàn)證本文動(dòng)態(tài)嵌套網(wǎng)格計(jì)算程序的準(zhǔn)確性，來流馬赫數(shù)為0.755，翼型攻角的變化規(guī)律是

式中 初始攻角(平均攻角)α=0.016°;振幅 αm=2.51°;無量綱角頻率 k=wc/u∞;c表示弦長;ω 是角頻率。

圖2給出了翼型的升力系數(shù)隨攻角的變化規(guī)律，并與Batina J T［11］的結(jié)果以及實(shí)驗(yàn)結(jié)果進(jìn)行了比較。從圖2中可看出，本文的計(jì)算結(jié)果與文獻(xiàn)［11］的計(jì)算結(jié)果吻合很好，但都與實(shí)驗(yàn)值稍有差異。Batina J T認(rèn)為導(dǎo)致這一差異的原因是實(shí)驗(yàn)數(shù)據(jù)可能在稍大的平均攻角下獲得的。

圖2 升力系數(shù)隨攻角的變化Fig.2 Lift coefficient vs angle of attack

圖3 計(jì)算區(qū)域網(wǎng)格(上)和邊界條件(下)Fig.3 Computational domain(upper)and prescribed boundary conditions(lower)

4 隔塞運(yùn)動(dòng)對噴管內(nèi)流場的影響及結(jié)果分析

4.1 物理模型和邊界條件

不考慮發(fā)動(dòng)機(jī)工作過程中噴管的形變，計(jì)算區(qū)域取噴管和一段圓柱形外流場區(qū)域。噴管入口直徑為220 mm，噴管喉部直徑為72 mm，擴(kuò)張比為2.5，收斂半角和擴(kuò)張半角分別為30°和15°，一般隔塞式脈沖發(fā)動(dòng)機(jī)多使用多個(gè)小隔塞，為了重點(diǎn)研究小隔塞在噴管內(nèi)運(yùn)動(dòng)對發(fā)動(dòng)機(jī)性能影響的規(guī)律，并簡化仿真模型，本文研究了一個(gè)小隔塞。在噴管內(nèi)運(yùn)動(dòng)的過程，隔塞的直徑為30 mm，長度為20 mm。外流場區(qū)域徑向取5倍的噴管出口直徑，軸向取10倍的噴管出口直徑。噴管入口的質(zhì)量流率為18 kg/s，總溫為3 300 K，總壓為10 MPa。由于噴管入口為亞音速，因此在計(jì)算中采用局部準(zhǔn)一維特征分析確定入口邊界條件，壁面采用無滑移邊界條件，絕熱壁假設(shè)，計(jì)算域出口邊界條件根據(jù)馬赫數(shù)判定。當(dāng)出口為超聲速時(shí)，此時(shí)所有物理量外推;當(dāng)出口為亞聲速時(shí)，給定環(huán)境反壓，其它參數(shù)由內(nèi)向外插值。

首先，計(jì)算定常狀態(tài)下噴管的流場;然后，讓隔塞從初始位置(此時(shí)記時(shí)間t=0)沿著軸線運(yùn)動(dòng)。根據(jù)文獻(xiàn)［12］的實(shí)驗(yàn)結(jié)果，隔塞從隔板飛出后，在一級燃燒室中的運(yùn)動(dòng)方向與燃燒室的軸線基本保持平行，隔塞到達(dá)噴管后，靠近中心線的一小部分隔塞直接排出噴管［12］，且實(shí)驗(yàn)得到隔塞到達(dá)噴管之前的飛行速度約為30 m/s。本文重點(diǎn)研究隔塞在噴管內(nèi)運(yùn)動(dòng)過程中噴管內(nèi)流場結(jié)構(gòu)的演化，暫不考慮隔塞在噴管中的運(yùn)動(dòng)規(guī)律對流場結(jié)構(gòu)的影響。因此，根據(jù)文獻(xiàn)［12］的實(shí)驗(yàn)，假設(shè)隔塞以30m/s的速度在噴管中運(yùn)動(dòng)。計(jì)算區(qū)域網(wǎng)格如圖3所示，包圍隔塞的貼體網(wǎng)格整體在噴管的網(wǎng)格上運(yùn)動(dòng)，總網(wǎng)格數(shù)目為49 654。

圖4 不同時(shí)刻的馬赫數(shù)云圖Fig.4 Mach number contour at different times

圖 5 馬赫數(shù)云圖(局部)t=3.1 ms(上)t=6.0 ms(下)Fig.5 Mach number contour(local)at 3.1 ms(upper)and 6.0 ms(lower)

圖6 t=3.1 ms時(shí)刻軸線上有無隔塞時(shí)馬赫數(shù)分布Fig.6 Axial mach number distributions with plug and without plug at 3.1 ms

4.2 噴管內(nèi)流場結(jié)構(gòu)演化分析

以馬赫數(shù)云圖來揭示隔塞運(yùn)動(dòng)過程中噴管內(nèi)流場的演變過程，并與不含隔塞的云圖進(jìn)行了對比。圖4所示的是隔塞在噴管內(nèi)運(yùn)動(dòng)的過程中，6個(gè)時(shí)刻的馬赫數(shù)云圖。圖4中，上半部分為噴管中含隔塞時(shí)不同時(shí)刻的馬赫數(shù)云圖;下半部分為噴管中不含隔塞的云圖。由圖4可見，圖4(a)、圖4(b)中，隔塞位于噴管的收斂段內(nèi)，距離噴管的喉部還有一段距離，其流場結(jié)構(gòu)與不含隔塞時(shí)的結(jié)構(gòu)相差不大。

在t=3.1 ms時(shí)刻(圖4(c)所示)，隔塞正好位于噴管的喉部。此時(shí)，隔塞前端的燃?xì)饬鳛閬喴羲?，后端的燃?xì)饬鳛槌羲?，附著在隔塞上的燃?xì)饬麟x開物面，流動(dòng)出現(xiàn)分離，并在分離區(qū)上游引起分離激波;同時(shí)，燃?xì)饬髟诜蛛x點(diǎn)后形成渦流，由于此時(shí)隔塞附近的流場為軸對稱的，因此隔塞后端分離形成的渦流交匯形成尾跡，很大程度上影響了燃?xì)饬髟趪姽軘U(kuò)張段的流場結(jié)構(gòu)，如圖5上半部分所示。軸線上有無隔塞時(shí)，燃?xì)饬黢R赫數(shù)對比如圖6所示。由圖6可見，由于隔塞的影響，噴管軸線上的馬赫數(shù)均比無隔塞的時(shí)候要小，尤其是在噴管的擴(kuò)張段，平均要低0.5;同時(shí)，隔塞后端的馬赫數(shù)先增大后減小，有一段很小的波動(dòng)，這是因?yàn)槿細(xì)庠诟羧暮蠖诵纬傻臏u流所引起的。

圖4(d)、圖4(e)、圖4(f)中，隔塞已經(jīng)完全通過了噴管喉部，位于噴管的擴(kuò)張段中，此時(shí)噴管擴(kuò)張段里的燃?xì)饬鳛槌羲?。因此，在隔塞的前端形成了弓形脫體激波，如圖5下半部分所示。隨著隔塞的后移，隔塞后端形成的尾跡與噴管出口的馬赫錐逐漸相互作用，形成了另一道激波，如圖4(e)、4(f)所示。同時(shí)，當(dāng)隔塞位于噴管擴(kuò)張段上游時(shí)，形成的弓形脫體激波還作用在噴管壁面上，造成此處壓力的變化，如圖7所示。由圖7可見，與初始時(shí)刻相比，隔塞在噴管內(nèi)運(yùn)動(dòng)的過程中，作用在噴管擴(kuò)張段壁面上的壓力劇烈波動(dòng)。

圖7 不同時(shí)刻噴管壁面上的壓力Fig.7 Pressure distributions in nozzle wall at different time

4.3 隔塞運(yùn)動(dòng)對發(fā)動(dòng)機(jī)推力的影響

火箭發(fā)動(dòng)機(jī)的推力定義為靜止條件下發(fā)動(dòng)機(jī)所受內(nèi)、外壁面壓強(qiáng)的合力，發(fā)動(dòng)機(jī)內(nèi)壁面壓強(qiáng)的合力由兩部分組成，即

式中 FM為作用在燃燒室殼體上的力，F(xiàn)N為作用在噴管上的力。

每種力都可由作用在壁面上的壓強(qiáng)和粘性力積分得到，但粘性力僅占到總推力的1%以下。因此，計(jì)算中忽略粘性力，同時(shí)還假設(shè)作用在燃燒室殼體上的力FM和發(fā)動(dòng)機(jī)外壁面受到的力保持不變。所以，在隔塞運(yùn)動(dòng)的過程中，只有FM變化。

圖8中所示的為發(fā)動(dòng)機(jī)無量綱推力隨時(shí)間的變化曲線(取火箭運(yùn)動(dòng)的方向?yàn)檎较?。從圖8中可看出，在隔塞沿著噴管軸線運(yùn)動(dòng)的過程中，發(fā)動(dòng)機(jī)的瞬時(shí)推力先逐漸減小，在t=3.1 ms時(shí)刻降到最小，并在接下來的一段時(shí)間內(nèi)，瞬時(shí)推力的大小變化不顯著，t=3.5 ms之后，瞬時(shí)推力開始急劇升高，并在t=4.5 ms時(shí)刻推力達(dá)到最大值，隨后緩慢下降，最終趨于初始時(shí)刻的值。隔塞的運(yùn)動(dòng)是造成發(fā)動(dòng)機(jī)瞬時(shí)推力變化的主要原因，隔塞通過噴管喉部后，由于噴管擴(kuò)張段中燃?xì)饬鳛槌羲?，因此在隔塞前端形成了弓形脫體激波，并與噴管壁面相互作用，造成噴管壁面上壓力的變化，如圖7所示。t=3.1 ms時(shí)刻，作用在噴管擴(kuò)張段內(nèi)的壓力急劇減小，而此時(shí)收斂段內(nèi)的壓力變化不大，積分后的合力減小。所以，此時(shí)瞬時(shí)推力降至最小。t=4.5 ms時(shí)刻，作用在噴管擴(kuò)張段內(nèi)的壓強(qiáng)顯著增大，積分后的瞬時(shí)推力也達(dá)到最大。

圖8 無量綱推力隨時(shí)間的變化曲線Fig.8 Normalized thrust vs time

同時(shí)，圖8中還顯示隔塞在噴管收斂段運(yùn)動(dòng)時(shí)，作用在噴管壁面上的力變化不大，發(fā)動(dòng)機(jī)的瞬時(shí)推力幾乎不變;隔塞位于噴管喉部時(shí)，發(fā)動(dòng)機(jī)瞬時(shí)推力降至最小，相比于初始時(shí)刻減小了27.7%;隔塞完全通過噴管喉部后，并位于噴管擴(kuò)張段上游時(shí)，發(fā)動(dòng)機(jī)瞬時(shí)推力達(dá)到最大值，比初始時(shí)刻增大了8.9%。一般隔塞式雙脈沖發(fā)動(dòng)機(jī)二脈沖工作時(shí)，一級燃燒室中的隔塞不止一個(gè)。因此，當(dāng)這些隔塞依次排出噴管時(shí)，會(huì)使發(fā)動(dòng)機(jī)瞬時(shí)推力發(fā)生巨大改變，造成推力-時(shí)間曲線的波動(dòng)，使發(fā)動(dòng)機(jī)的內(nèi)彈道性能難以重現(xiàn)，嚴(yán)重影響發(fā)動(dòng)機(jī)的正常工作。因此，有必要針對隔塞的形狀及運(yùn)動(dòng)規(guī)律進(jìn)行大量數(shù)值計(jì)算，對其特性進(jìn)行詳細(xì)研究，以總結(jié)出規(guī)律。

5 結(jié)論

(1)本文所采用的動(dòng)態(tài)結(jié)構(gòu)嵌套網(wǎng)格方法能處理運(yùn)動(dòng)幅度較大的相對運(yùn)動(dòng)問題，具有較高的工程應(yīng)用價(jià)值，為數(shù)值模擬隔塞在一級燃燒室中的運(yùn)動(dòng)軌跡奠定基礎(chǔ)。

(2)隔塞通過噴管喉部后，燃?xì)饬髟诟羧岸诵纬晒蚊擉w激波，在隔塞后端形成尾跡，作用在噴管收斂段上的壓強(qiáng)變化顯著。

(3)第二級燃燒室中的多個(gè)隔塞依次從噴管排出的過程中，發(fā)動(dòng)機(jī)的瞬時(shí)推力逐漸增大，到達(dá)最大值后保持一段很小的時(shí)間，然后開始下降，最終緩慢增加至初始時(shí)刻的值。

(4)發(fā)動(dòng)機(jī)正常工作時(shí)，當(dāng)直徑為30 mm的隔塞運(yùn)動(dòng)到直徑為72 mm噴管喉部時(shí)，發(fā)動(dòng)機(jī)瞬時(shí)推力降至最小，比初始時(shí)刻減小了27.7%;隔塞完全通過噴管喉部后，并位于噴管擴(kuò)張段上游時(shí)，發(fā)動(dòng)機(jī)瞬時(shí)推力達(dá)到最大值，比初始時(shí)刻增大了8.9%。

［1］Naumann K W，Stadler L.Double-pulse solid rocket motor technology-applications and technical solutions［R］.AIAA 2010-6754.

［2］Nishii S，F(xiàn)ukuda K，Kubota N.Combustion tests of two-stage pulse rocket［R］.AIAA 89-2426.

［3］Wang Chang-hui，Liu Yu ，Liu Ya-bing.Design and experimental studies on ceramic port cover for dual pulse motor［J］.Acta.Astronautica，2011，68(11):1881-1890.

［4］張涵信，陳堅(jiān)強(qiáng)，高樹椿.H2/O2燃燒的超聲速非平衡流動(dòng)的數(shù)值模擬［J］.宇航學(xué)報(bào)，1994，15(2):14-23.

［5］Menter F R.Two-equation eddy-viscosity turbulence models for engineering applications［J］.AIAA Journal，1994，32(08):1598-1605.

［6］Liou Meng-sing.Methods for the accurate computations of hypersonic flows［J］.Journal of Computational Physics，2001，174(1):38-80.

［7］Jameson Antony.Time dependent calculations using multigrid，with applications to unsteady flows past airfoils and wings［R］.AIAA 91-1596.

［8］Chiu Ing-Tsau ，Robert L Meakin.On automating domain connectivity for overset grids［R］.AIAA 95-0854.

［9］Wang Z J，Parthasarathy V.A fully automated Chimera methodology for multiple moving body problems［J］.International Journal for Numerical Methods in Fluids，2000，33(7):919-938.

［10］John B McDevitt，Arthur F Okuno.Static and dynamic pressure measurements on a NACA0012 airfoil in the ames high Reynolds number facility［R］.NACA-TP-2485.

［11］John T Batina.Unsteady euler airfoil solutions using unstructured dynamic meshes［J］.AIAA Journal ，1990，28(9):1381-1388.

［12］李江，肖育民，何國強(qiáng)，等.雙脈沖固體火箭發(fā)動(dòng)機(jī)二次點(diǎn)火內(nèi)視研究［J］.推進(jìn)技術(shù)，1998，19(3):61-64.