基于嵌套網格的脈沖發動機噴管內流場數值模擬①

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

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

0 引言

固體火箭發動機具有結構簡單、使用方便、可靠性好等優點，在國防領域受到廣泛重視及應用。與傳統的單室單推、單室雙推發動機相比，雙脈沖固體火箭發動機具有改善固體火箭發動機能量可控性的突出優勢，可多次點火，并提供不連續推力，有利于提高導彈武器的作戰能力［1］。

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

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

1 控制方程及數值求解方法

1.1 流動控制方程

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

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

式中 ρ為氣體密度;u、v為氣體運動速度矢量的2個分量;xt、yt分別為網格運動速度矢量的2個分量;T為氣體的溫度;E為單位體積氣體的總能量;τ表示應力張量，其具體形式參考文獻［4］。

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

式中 k為湍動能;ω為比耗散率;μt為湍流粘性系數，其他它參數的具體形式見參考文獻［5］。

1.2 數值求解方法

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

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

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

2 嵌套網格技術

隔塞在噴管內運動的過程為典型的移動邊界問題。目前，處理移動邊界問題較成熟的方法包括非結構動網格技術及動態嵌套網格技術等，但非結構網格在處理具有較大相對位移問題時，網格變形與重構的表現能力較差。因此，本文選擇動態嵌套網格方法，其優點是對于復雜幾何體，不要求流場各個計算域共享邊界，減輕了網格生成難度，適用于運動幅度較大的相對運動。

2.1 洞邊界的確定

在嵌套網格系統中，建立人工插值內邊界的過程稱之為挖洞過程。為了提高確定洞邊界的效率，本文采用 Chiu I T和 Meakin R L提出的“Hole-Map”方法［8］，其核心思想是對于給定的嵌套網格體系，若已知其拓撲結構，就能用均勻的笛卡爾網格單元去近似每個曲面邊界，從而得到該曲面邊界的笛卡爾網格近似。在此基礎上，網格點P與曲面邊界的位置關系可近似轉換為網格點P與笛卡爾網格的關系，而這種關系十分容易確定。

2.2 貢獻單元的搜尋

在嵌套網格系統中，通過貢獻單元將流場解的信息插值到插值邊界面網格上。因此，尋找貢獻單元技術的優劣，對整個嵌套網格方法有很大影響，有時甚至是嵌套網格方法成敗的關鍵。本文結合Wang Z J的矢量判別法［9］，提出了一種簡單易行的尋找插值點貢獻單元的方法，具體步驟如下:

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

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

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

3 算例驗證

3.1 靜態NACA0012翼型繞流

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

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

3.2 俯仰振動的NACA0012翼型繞流

以NACA0012翼型繞1/4弦點作簡諧俯仰振動為算例，驗證本文動態嵌套網格計算程序的準確性，來流馬赫數為0.755，翼型攻角的變化規律是

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

圖2給出了翼型的升力系數隨攻角的變化規律，并與Batina J T［11］的結果以及實驗結果進行了比較。從圖2中可看出，本文的計算結果與文獻［11］的計算結果吻合很好，但都與實驗值稍有差異。Batina J T認為導致這一差異的原因是實驗數據可能在稍大的平均攻角下獲得的。

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

圖3 計算區域網格(上)和邊界條件(下)Fig.3 Computational domain(upper)and prescribed boundary conditions(lower)

4 隔塞運動對噴管內流場的影響及結果分析

4.1 物理模型和邊界條件

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

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

圖4 不同時刻的馬赫數云圖Fig.4 Mach number contour at different times

圖 5 馬赫數云圖(局部)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時刻軸線上有無隔塞時馬赫數分布Fig.6 Axial mach number distributions with plug and without plug at 3.1 ms

4.2 噴管內流場結構演化分析

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

在t=3.1 ms時刻(圖4(c)所示)，隔塞正好位于噴管的喉部。此時，隔塞前端的燃氣流為亞音速，后端的燃氣流為超音速，附著在隔塞上的燃氣流離開物面，流動出現分離，并在分離區上游引起分離激波;同時，燃氣流在分離點后形成渦流，由于此時隔塞附近的流場為軸對稱的，因此隔塞后端分離形成的渦流交匯形成尾跡，很大程度上影響了燃氣流在噴管擴張段的流場結構，如圖5上半部分所示。軸線上有無隔塞時，燃氣流馬赫數對比如圖6所示。由圖6可見，由于隔塞的影響，噴管軸線上的馬赫數均比無隔塞的時候要小，尤其是在噴管的擴張段，平均要低0.5;同時，隔塞后端的馬赫數先增大后減小，有一段很小的波動，這是因為燃氣在隔塞的后端形成的渦流所引起的。

圖4(d)、圖4(e)、圖4(f)中，隔塞已經完全通過了噴管喉部，位于噴管的擴張段中，此時噴管擴張段里的燃氣流為超音速。因此，在隔塞的前端形成了弓形脫體激波，如圖5下半部分所示。隨著隔塞的后移，隔塞后端形成的尾跡與噴管出口的馬赫錐逐漸相互作用，形成了另一道激波，如圖4(e)、4(f)所示。同時，當隔塞位于噴管擴張段上游時，形成的弓形脫體激波還作用在噴管壁面上，造成此處壓力的變化，如圖7所示。由圖7可見，與初始時刻相比，隔塞在噴管內運動的過程中，作用在噴管擴張段壁面上的壓力劇烈波動。

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

4.3 隔塞運動對發動機推力的影響

火箭發動機的推力定義為靜止條件下發動機所受內、外壁面壓強的合力，發動機內壁面壓強的合力由兩部分組成，即

式中 FM為作用在燃燒室殼體上的力，FN為作用在噴管上的力。

每種力都可由作用在壁面上的壓強和粘性力積分得到，但粘性力僅占到總推力的1%以下。因此，計算中忽略粘性力，同時還假設作用在燃燒室殼體上的力FM和發動機外壁面受到的力保持不變。所以，在隔塞運動的過程中，只有FM變化。

圖8中所示的為發動機無量綱推力隨時間的變化曲線(取火箭運動的方向為正方向)。從圖8中可看出，在隔塞沿著噴管軸線運動的過程中，發動機的瞬時推力先逐漸減小，在t=3.1 ms時刻降到最小，并在接下來的一段時間內，瞬時推力的大小變化不顯著，t=3.5 ms之后，瞬時推力開始急劇升高，并在t=4.5 ms時刻推力達到最大值，隨后緩慢下降，最終趨于初始時刻的值。隔塞的運動是造成發動機瞬時推力變化的主要原因，隔塞通過噴管喉部后，由于噴管擴張段中燃氣流為超音速，因此在隔塞前端形成了弓形脫體激波，并與噴管壁面相互作用，造成噴管壁面上壓力的變化，如圖7所示。t=3.1 ms時刻，作用在噴管擴張段內的壓力急劇減小，而此時收斂段內的壓力變化不大，積分后的合力減小。所以，此時瞬時推力降至最小。t=4.5 ms時刻，作用在噴管擴張段內的壓強顯著增大，積分后的瞬時推力也達到最大。

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

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

5 結論

(1)本文所采用的動態結構嵌套網格方法能處理運動幅度較大的相對運動問題，具有較高的工程應用價值，為數值模擬隔塞在一級燃燒室中的運動軌跡奠定基礎。

(2)隔塞通過噴管喉部后，燃氣流在隔塞前端形成弓形脫體激波，在隔塞后端形成尾跡，作用在噴管收斂段上的壓強變化顯著。

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

(4)發動機正常工作時，當直徑為30 mm的隔塞運動到直徑為72 mm噴管喉部時，發動機瞬時推力降至最小，比初始時刻減小了27.7%;隔塞完全通過噴管喉部后，并位于噴管擴張段上游時，發動機瞬時推力達到最大值，比初始時刻增大了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，Fukuda 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］張涵信，陳堅強，高樹椿.H2/O2燃燒的超聲速非平衡流動的數值模擬［J］.宇航學報，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］李江，肖育民，何國強，等.雙脈沖固體火箭發動機二次點火內視研究［J］.推進技術，1998，19(3):61-64.