卓長飛，封 鋒，武曉松

(南京理工大學(xué)機械工程學(xué)院，南京 210094)

0 引言

研究火箭燃氣射流在工程技術(shù)領(lǐng)域有著重要意義，如機載火箭導(dǎo)彈發(fā)射時燃氣射流所形成的壓力擾動往往使得進氣道內(nèi)增壓，它產(chǎn)生的溫度與壓力畸變都是引起壓氣機失速或發(fā)動機停車的重要因素。燃氣射流的吸入與藥柱燃燒后產(chǎn)物的再燃燒也對飛機發(fā)動機正常工作有影響。此外，發(fā)動機燃氣射流對飛行火箭外部流分離和底壓也有影響，可以改變它們的空氣動力特性［1］。

近年來，火箭燃氣射流動力學(xué)問題的數(shù)值模擬水平不斷提高，使得火箭燃氣射流動力學(xué)獲得迅速發(fā)展。然而，在數(shù)值研究火箭燃氣射流流場方面，均是建立在網(wǎng)格算法基礎(chǔ)上［2－4］。實際上，計算流體力學(xué)的數(shù)值模擬方法可劃分為2種:網(wǎng)格算法和無網(wǎng)格算法。網(wǎng)格算法分為非結(jié)構(gòu)網(wǎng)格算法和結(jié)構(gòu)網(wǎng)格算法，網(wǎng)格算法發(fā)展較早，技術(shù)比較成熟。相比之下，新興的無網(wǎng)格算法在非結(jié)構(gòu)網(wǎng)格算法打破網(wǎng)格結(jié)構(gòu)化思維的基礎(chǔ)上，進一步拋棄了網(wǎng)格化的思維，徹底地打破了結(jié)構(gòu)化和網(wǎng)格化的思想約束，放棄了網(wǎng)格而不需要像網(wǎng)格方法那樣形成網(wǎng)格單元。無網(wǎng)格算法的區(qū)域離散只涉及離散點，在點云結(jié)構(gòu)上離散控制方程，從而不存在網(wǎng)格質(zhì)量和拓撲結(jié)構(gòu)限制等問題。因此，無網(wǎng)格算法比網(wǎng)格算法具有更大的幾何靈活性，在處理復(fù)雜外形和復(fù)雜流動過程方面具有獨特的優(yōu)勢。在過去的20多年中，國外計算流體力學(xué)研究者對無網(wǎng)格算法的格式精度和計算效率進行了大量研究［5－7］。而國內(nèi)對無網(wǎng)格算法研究較少，但無網(wǎng)格算法仍在不斷發(fā)展之中，如將基于網(wǎng)格算法的數(shù)值方法，如Roe格式、AUSM+－UP格式等迎風(fēng)格式推廣到無網(wǎng)格算法［8－9］、無網(wǎng)格的隱式算法研究［8］、無網(wǎng)格算法的精度分析［10、無網(wǎng)格在移動邊界流場中的應(yīng)用研究［11］，并逐漸將無網(wǎng)格算法用于工程問題［12－13］。

火箭燃氣射流含有富燃氣體，能與外流中的空氣發(fā)生化學(xué)反應(yīng)，整個流場是典型的化學(xué)非平衡流問題。為了研究采用無網(wǎng)格算法進行火箭燃氣射流流場的可行性，本文詳細給出在無網(wǎng)格條件下求解帶化學(xué)反應(yīng)的多組分Euler方程具體過程，并驗證了采用無網(wǎng)格算法進行化學(xué)非平衡流數(shù)值模擬的可行性。最后，對不同工況火箭燃氣射流流場進行了數(shù)值研究。結(jié)果表明，無網(wǎng)格方法模擬火箭燃氣射流流場是可行的。

1 數(shù)值方法

1．1 控制方程

在笛卡爾坐標系下，微分守恒形式的二維軸對稱化學(xué)非平衡流Euler方程為

式中 Q為守恒變量;F、G分別為軸向和徑向無粘通量;H為軸對稱源項;S為化學(xué)反應(yīng)源項。

這里僅列出主要的變量:

式中 ρ為密度;u、v分別為軸向、徑向速度;p為壓力;E 為單位體積的總能;ci(i=1，．．．，N－1)為 i組分質(zhì)量分數(shù);ωi為i組分生成速率;N為總組分數(shù)。

1．2 點云的概念

無網(wǎng)格算法只需將計算域離散成一系列點，并按一定規(guī)律組織形成點云，計算域內(nèi)每個點都有自己的點云結(jié)構(gòu)。流動計算過程中，控制方程的離散即建立在點云基礎(chǔ)上。計算域內(nèi)某點的點云結(jié)構(gòu)是由該點和周圍的一些相鄰點所構(gòu)成，該點稱為中心點，而周圍相鄰點稱為衛(wèi)星點。某一點云結(jié)構(gòu)見圖1，點i的點云結(jié)構(gòu)由中心點i和6個衛(wèi)星點組成。

圖1 離散域內(nèi)節(jié)點i的點云結(jié)構(gòu)Fig．1 The structure for clouds of point i in discrete domin

1．3 無網(wǎng)格空間導(dǎo)數(shù)逼近原理

對控制方程進行離散，首先利用曲線擬合方法，確定每個點云結(jié)構(gòu)內(nèi)物理變量的空間導(dǎo)數(shù)。本文采用一階線性函數(shù)進行曲線擬合，函數(shù)為

每個點云的中心點及其衛(wèi)星點都滿足該點云結(jié)構(gòu)的解函數(shù):

式中 k表示i的點云結(jié)構(gòu)中衛(wèi)星點總個數(shù)。

則存在如下矩陣關(guān)系:

解該矛盾方程組，得到該點云結(jié)構(gòu)內(nèi)a0、a1、a2的最小二乘解。則點云結(jié)構(gòu)中i處的空間導(dǎo)數(shù)可表示為

式中 下標j表示i點的點云結(jié)構(gòu)中衛(wèi)星點。

由此可確定該點云結(jié)構(gòu)內(nèi)中心點與每個衛(wèi)星點之間的傳播系數(shù) bij、cij。

多組分化學(xué)非平衡流Euler方程在點云i上半離散形式為

由式(6)可得到多組分化學(xué)非平衡流Euler方程的對流項空間導(dǎo)數(shù)為

其中，Eij表示中心點i與衛(wèi)星點j之間的通量，其值由i和j的守恒量得到，如圖2所示。求解該通量采用高精度高分辨率的AUSMPW+迎風(fēng)格式［14］。在求解通量之前，首先對點i和j的守恒變量進行重構(gòu)到中點的左、右間斷處:

式中 φ表示抑制數(shù)值振蕩的限制器;r表示矢徑。

守恒變量的梯度則由式(6)計算得:

圖2 中心點和衛(wèi)星點通量重構(gòu)Fig．2 The flux reconstruction between center point and satellite point

1．4 化學(xué)反應(yīng)動力學(xué)模型與剛性問題

化學(xué)反應(yīng)動力學(xué)模型是決定成功模擬化學(xué)非平衡流的關(guān)鍵之一。常規(guī)的雙基固體推進劑產(chǎn)生的燃氣主要是由 CO、H2、CO2、H2O、N2、HCl等組成，因此采用CO—H2—O2化學(xué)反應(yīng)系統(tǒng)。本文采用8組分(CO、H2、O2、CO2、H2O、H、OH、O)12 個基元反應(yīng)的 CO—H2—O2系統(tǒng)化學(xué)反應(yīng)模型［15］，基元反應(yīng)表達式和系數(shù)見表1。表1中，M代表第三碰撞體;A為指前因子，其中n為基元反應(yīng)級數(shù);b為溫度指數(shù);E為活化能。

表1 CO—H2—O2基元反應(yīng)模型Table 1 Detailed reaction model for CO—H2—O2

化學(xué)非平衡流控制方程組分為流動部分和化學(xué)反應(yīng)部分，兩者相互耦合，并會產(chǎn)生剛性問題。本文采用時間算子分裂算法處理這種耦合過程，即首先計算流體流動效應(yīng)，得到物理變量的過渡值。然后，繼續(xù)計算，將化學(xué)反應(yīng)的貢獻疊加到物理變量過渡值，最終得到下一時刻體現(xiàn)整體效應(yīng)的物理量值。其中，在計算化學(xué)反應(yīng)對流場的貢獻時，需把求解流動偏微分方程時采用的時間步長進一步細分，作為求解化學(xué)反應(yīng)剛性常微分方程的步長。具體做法是先凍結(jié)化學(xué)反應(yīng)，求解得到流場參數(shù);然后，將化學(xué)反應(yīng)看做等容吸熱或放熱過程，保持內(nèi)能、速度參數(shù)不變，計算各組分的質(zhì)量變化率;最后，迭代求解溫度。

1．5 物理模型與邊界條件

本文取燃氣射流場軸向計算長度為彈體直徑的40倍，徑向計算長度彈體直徑的10倍?；鸺細馍淞髁鲌鲇嬎銋^(qū)域中噴管出口附近空間布點如圖3所示，布點總數(shù)為28 049個。在本文計算中，不考慮火箭發(fā)動機燃燒室和噴管內(nèi)流動，直接在噴管出口給定燃氣超聲速射流條件，并保持噴管出口燃氣射流的馬赫數(shù)Maj=2．5、靜溫Tj=1 500 K、燃氣組分與質(zhì)量分數(shù)(如表2所示)不變。

表2 燃氣主要組分與質(zhì)量分數(shù)Table 2 The main mass fraction of rocket gas

圖3 射流流場空間布點示意圖(噴管出口附近)Fig．3 Clouds of points for rocket gas efflux field

2 數(shù)值驗證

2．1 NACA0012翼型跨聲速繞流

首先，針對NACA0012翼型的跨聲速流場進行計算。圖4給出了NACA0012翼型的空間布點局部放大圖。來流條件為馬赫數(shù)0．8，攻角1．25°，計算區(qū)域布點總數(shù)為6 327個。該算例為二維平面問題，計算時關(guān)閉軸對稱源項，且整個流場氣體為理想完全氣體，不考慮化學(xué)反應(yīng)。計算得到的翼型上下表面壓力系數(shù)圖5所示。由圖5可看出，跨聲速流場存在復(fù)雜的波系，上下表面的壓力系數(shù)與試驗值吻合較好，在激波處沒有出現(xiàn)數(shù)值振蕩，這說明本文發(fā)展的無網(wǎng)格方法數(shù)值模擬精度較高，具有模擬復(fù)雜波系的能力。

圖4 NACA0012翼型空間布點示意圖(表面局部)Fig．4 Clouds of points for NACA0012 airfoils(surface)

圖5 NACA0012翼型上下表面壓力系數(shù)分布Fig．5 Pressure coefficient distributions of NACA0012 airfoils

2．2 楔體誘導(dǎo)斜爆轟流場

計算工況:混合氣體為2H2+O2+3．76N2，來流馬赫數(shù)7．5，靜溫 293 K，靜壓4 000 Pa，尖劈角度為 25°，H2－O2化學(xué)動力學(xué)模型采用7組分8反應(yīng)模型［16］。計算區(qū)域布點如圖6所示，布點總數(shù)為15 074個。

計算得到的溫度場如圖7所示。高速可燃氣體通過尖劈會發(fā)生偏轉(zhuǎn)，形成斜激波，波后混合氣體受斜激波壓縮的作用，壓力和溫度有所提高，但不能直接點燃混合氣體。壓縮后的混合氣體經(jīng)過一定的誘導(dǎo)期后形成斜爆轟波。整個流場是由斜激波、橫波、斜爆轟波、滑移線組成，這些波系相交于三波點。本文采用無網(wǎng)格算法的數(shù)值結(jié)果與文獻［17］數(shù)值計算得到的溫度云圖完全吻合，各種波系均清晰可見，說明無網(wǎng)格算法能適用于具有復(fù)雜化學(xué)非平衡流場的數(shù)值模擬。

圖6 楔體誘導(dǎo)斜爆轟流場的空間布點示意圖Fig．6 Clouds of points for wedge oblique detonation

圖7 斜爆轟流場溫度云圖與流場結(jié)構(gòu)Fig．7 Temperature and fields structure of oblique detonation fields

3 計算結(jié)果與分析

3．1 欠膨脹狀態(tài)下射流流場

計算工況為來流馬赫數(shù)2．0，來流靜壓為0．4 atm，火箭噴管出口的燃氣射流壓力分別為 1．0、2．5、5．0 atm。當噴管出口的燃氣射流壓力大于外流壓力時為欠膨脹狀態(tài)，且隨著燃氣射流壓力的增大，欠膨脹程度增加。不同噴管出口燃氣射流壓力下燃氣射流流場的馬赫數(shù)云圖如圖8所示，軸線馬赫數(shù)、壓力、溫度如圖9所示。

圖8 欠膨脹狀態(tài)下射流流場馬赫數(shù)云圖Fig．8 Mach contours of field flow at under-expanded state

由圖8看出，相交射流激波、反射激波、射流邊界等主要流動特征非常清晰，且由圖9中軸線參數(shù)分布可看出，在軸線上相交射流激波第1次相交形成的強間斷處，各參數(shù)均發(fā)生突變，具有較大的梯度變化率，這完全符合激波物理性質(zhì)。這些均說明本文采用無網(wǎng)格方法模擬得到的激波結(jié)構(gòu)非常清晰且合理，具有較強的捕捉強間斷流場的能力，能較好地模擬火箭燃氣射流流場。

欠膨脹狀態(tài)下，燃氣流出噴管后在噴管唇口附近發(fā)散成一束扇形膨脹波族，形成Prantl－Meyer流，氣體將繼續(xù)膨脹加速，流動參數(shù)中的馬赫數(shù)升高，壓力、溫度降低。同時，燃氣流出噴管后，會受到外流的壓縮作用，而形成相交射流激波，上、下相交射流激波相交于軸線上，并發(fā)生發(fā)射形成反射激波，反射激波又在混合層發(fā)生反射。這樣，在射流軸線上重復(fù)發(fā)生著膨脹－壓縮過程，且強度不斷減弱，這個現(xiàn)象可由軸線上各參數(shù)變化看出，各流動參數(shù)在一定范圍內(nèi)呈振蕩衰減型，最后趨于穩(wěn)定變化。

由圖8、圖9可看出，相交射流激波在軸線上的交點，隨著欠膨脹程度的增加而遠離噴管出口。這是由于欠膨脹程度的增加，導(dǎo)致燃氣射流對外流的作用增強，即相交射流激波強度增強，射流激波張角增大，從而交點更靠后。由流動參數(shù)軸線分布圖還可看出，隨著欠膨脹程度的增加，燃氣流出噴管后膨脹程度也大大提高，在相交射流激波前的馬赫數(shù)更高，壓力、溫度更低。

射流場軸線上離開噴管一定距離后，不同噴管出口燃氣射流壓力對應(yīng)下的軸線壓力也基本一致，且均約為外流的靜壓。這主要是由于燃氣射流與外流之間混合層的相互作用，通過剪切層的傳輸，導(dǎo)致燃氣射流軸線上的靜壓與外流的靜壓保持一致，而馬赫數(shù)、溫度仍保持較高水平。可見，火箭燃氣射流對尾流場的影響范圍非常大。

圖9 欠膨脹狀態(tài)下軸線上主要參數(shù)分布Fig．9 Distribution of main parameter along the axis at under-expanded state

3．2 過膨脹狀態(tài)下射流流場

計算工況為來流馬赫數(shù)2．0，來流靜壓為1．0個大氣壓，火箭噴管出口的燃氣射流壓力分別為0．8、0．6、0．4 atm。當噴管出口的燃氣射流壓力小于外流壓力時屬于過膨脹狀態(tài)，且隨著燃氣射流壓力的降低，過膨脹程度將增加。不同噴管出口燃氣射流壓力下射流流場的馬赫數(shù)云圖如圖10所示，軸線馬赫數(shù)、壓力、溫度如圖11所示。

圖10 不同過膨脹狀態(tài)下射流流場馬赫數(shù)云圖Fig．10 Mach contours of field flow at over-expanded state

圖11 過膨脹狀態(tài)下軸線上主要參數(shù)分布Fig．11 Distribution of main parameter along the axis at over-expanded state

在過膨脹狀態(tài)下，主要流動特征仍然非常清晰可見。從流場結(jié)構(gòu)來看，過膨脹和欠膨脹狀態(tài)相比的重要區(qū)別之一是過膨脹狀態(tài)下相交射流激波比較平直，而欠膨脹狀態(tài)下相交射流激波呈曲線。

隨著過膨脹程度的增大，相交射流激波在軸線上的交點不斷靠近噴管。這是由于隨著燃氣射流壓力降低，外流對燃氣射流的“擠壓”作用更明顯，從而使相交射流激波在軸線上的交點更靠前。隨著過膨脹程度的增大，射流激波會在軸線以及射流與外流的混合層處發(fā)生的反射次數(shù)增多，軸線上的參數(shù)也明顯表現(xiàn)為呈振蕩衰減型，最后趨于穩(wěn)定。值得一提的是在燃氣射流壓力為0．8 atm時，由軸線參數(shù)分布看出，在相交激波交點前，燃氣射流會先進行一定程度的膨脹加速，馬赫數(shù)升高，壓力、溫度降低。這可能是由于噴管出口燃氣射流與外流壓力很接近，處于過膨脹和欠膨脹之間，但不是完全膨脹狀態(tài)，因而產(chǎn)生既含欠膨脹特征，也含過膨脹特征的現(xiàn)象。

3．3 來流馬赫數(shù)對射流流場的影響

計算工況為來流馬赫數(shù)分別為 1．5、2．0、2．5 Ma，來流靜壓為0．4 atm，火箭噴管出口的燃氣射流壓力為1 atm。不同來流馬赫數(shù)下射流流場的馬赫數(shù)云圖如圖12所示，軸線馬赫數(shù)、壓力、溫度如圖13所示。

由計算結(jié)果可看出，超聲速來流條件下，僅改變不同來流馬赫數(shù)，對射流流場總體結(jié)構(gòu)基本無影響，相交射流激波相交點位置基本相同，軸線上參數(shù)分布也基本一致。值得注意的是離開噴管一定距離后，由于射流與外流的充分混合和發(fā)展，軸線上壓力會以不同程度降低，逐漸與外流保持一致。

圖12 不同來流馬赫數(shù)下射流流場馬赫數(shù)云圖Fig．12 Mach number contours of field flow variation with Mach number of incoming flow

4 結(jié)論

(1)發(fā)展的無網(wǎng)格計算方法能有效地進行化學(xué)非平衡流的數(shù)值模擬，能準確反映流場的流動特性，具有布點靈活、求解效率高、精確模擬化學(xué)非平衡流中如激波等復(fù)雜流動現(xiàn)象的特點。采用無網(wǎng)格方法模擬火箭燃氣射流流場，為數(shù)值研究燃氣射流流場開辟了一種新途徑。

(2)在欠膨脹狀態(tài)下，相交射流激波在軸線上的交點隨著欠膨脹程度的增加而遠離噴管出口，且相交射流激波呈曲線;在過膨脹狀態(tài)下，相交射流激波在軸線上的交點，隨著過膨脹程度的增加而不斷靠近噴管出口，且相交射流激波呈直線。

(3)在欠膨脹或過膨脹狀態(tài)下，射流流場軸線上離開噴管一定距離后，不同噴管出口燃氣射流壓力對應(yīng)下軸線壓力也基本一致，且均約為外流的靜壓。

(4)超聲速來流條件下，僅改變不同來流馬赫數(shù)，對射流流場總體結(jié)構(gòu)基本無影響，相交射流激波相交點位置基本相同，軸線上參數(shù)分布也基本一致。

圖13 不同來流馬赫數(shù)下軸線主要參數(shù)分布Fig．13 Distribution of main parameter along the axis variation with Mach number of incoming flow

［1］ 張福祥．火箭燃氣射流動力學(xué)［M］．哈爾濱:哈爾濱工業(yè)大學(xué)出版社，2004．

［2］ 孫振華，徐東來，何國強．飛行參數(shù)對導(dǎo)彈發(fā)動機羽流的影響［J］．固體火箭技術(shù)，2005，28(3):188－191．

［3］ 聶赟，樂貴高，馬大為．矢通量分裂法在噴管超聲速噴流計算中的應(yīng)用［J］．固體火箭技術(shù)，2013，36(1):56－60．

［4］ 常見虎，周長省，李軍．高低空環(huán)境下火箭發(fā)動機射流場的數(shù)值分析［J］．系統(tǒng)仿真學(xué)報，2007，19(16):3672－3675．

［5］ Balasubramanyam S，Raghurama Rao S V．A grid free upwind relaxation scheme for inviscid compressible flows［J］．International Journal for Numerical Methods in Fluids，2006，51(2):159－196．

［6］ Luo H，Baum J D，L?hner R．A hybrid Cartesian grid and gridless method for compressible flows［J］．Journal of Computational Physics，2006，214(2):618－632．

［7］ Jaisankar S，Shivashankar K，Rao S V R．A Grid－free central scheme for inviscid compressible flows［R］．AIAA 2007－3946．

［8］ 陳紅全．隱式無網(wǎng)格算法及其應(yīng)用研究［J］．空氣動力學(xué)報，2002，20(2)．

［9］ 孫迎丹，王剛，葉正寅．AUSM+－UP格式在無網(wǎng)格算法中的推廣［J］．空氣動力學(xué)報．2005，23(4)．

［10］ 胡世祥，李磊，佘春東，等．二維Euler方程無網(wǎng)格算法的精度分析［J］．計算力學(xué)學(xué)報，2005，22(2):232－236．

［11］ 周星，許厚謙．計算含運動邊界非定常流動的無網(wǎng)格算法［J］．力學(xué)與實踐．2010，32(3):16－25．

［12］ 馬新建，譚俊杰，任登鳳．三維無網(wǎng)格及其在超聲速彈丸流場模擬中的應(yīng)用［J］．彈道學(xué)報．2010，22(3):54－57．

［13］ 周星，許厚謙．無網(wǎng)格算法在膛口流場模擬中的應(yīng)用［J］．彈道學(xué)報．2011，23(1):24－26．

［14］ Kim K H．Accurate computation of hypersonic flows using AUSMPW+scheme and shock－induced grid technique［R］．AIAA 98－2442．

［15］ Charles J．Nietubica．Navier－stokes computations for a reacting，M864 base bleed projectile［R］．AIAA 93－0504．

［16］ 劉君，周松柏，徐春光．超聲速流動中燃燒現(xiàn)象的數(shù)值模擬方法及應(yīng)用［M］．長沙:國防科大出版社，2008．

［17］ Viguier C，Silva L F F，Desbordes D，et al．Onset of oblique detonation wave:comparison between experimental and numerical results for hydrogen－air mixtures［C］//Symposium(International)on Combustion．Elsevier，1996，26(2):3023－3031．