火箭燃?xì)馍淞髁鲌龅臒o網(wǎng)格方法模擬①

卓長飛，封 鋒，武曉松

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

0 引言

研究火箭燃?xì)馍淞髟诠こ碳夹g(shù)領(lǐng)域有著重要意義，如機(jī)載火箭導(dǎo)彈發(fā)射時燃?xì)馍淞魉纬傻膲毫_動往往使得進(jìn)氣道內(nèi)增壓，它產(chǎn)生的溫度與壓力畸變都是引起壓氣機(jī)失速或發(fā)動機(jī)停車的重要因素。燃?xì)馍淞鞯奈肱c藥柱燃燒后產(chǎn)物的再燃燒也對飛機(jī)發(fā)動機(jī)正常工作有影響。此外，發(fā)動機(jī)燃?xì)馍淞鲗︼w行火箭外部流分離和底壓也有影響，可以改變它們的空氣動力特性［1］。

近年來，火箭燃?xì)馍淞鲃恿W(xué)問題的數(shù)值模擬水平不斷提高，使得火箭燃?xì)馍淞鲃恿W(xué)獲得迅速發(fā)展。然而，在數(shù)值研究火箭燃?xì)馍淞髁鲌龇矫?，均是建立在網(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ǔ)上，進(jìn)一步拋棄了網(wǎng)格化的思維，徹底地打破了結(jié)構(gòu)化和網(wǎng)格化的思想約束，放棄了網(wǎng)格而不需要像網(wǎng)格方法那樣形成網(wǎng)格單元。無網(wǎng)格算法的區(qū)域離散只涉及離散點，在點云結(jié)構(gòu)上離散控制方程，從而不存在網(wǎng)格質(zhì)量和拓?fù)浣Y(jié)構(gòu)限制等問題。因此，無網(wǎng)格算法比網(wǎng)格算法具有更大的幾何靈活性，在處理復(fù)雜外形和復(fù)雜流動過程方面具有獨(dú)特的優(yōu)勢。在過去的20多年中，國外計算流體力學(xué)研究者對無網(wǎng)格算法的格式精度和計算效率進(jìn)行了大量研究［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］。

火箭燃?xì)馍淞骱懈蝗細(xì)怏w，能與外流中的空氣發(fā)生化學(xué)反應(yīng)，整個流場是典型的化學(xué)非平衡流問題。為了研究采用無網(wǎng)格算法進(jìn)行火箭燃?xì)馍淞髁鲌龅目尚行?，本文詳?xì)給出在無網(wǎng)格條件下求解帶化學(xué)反應(yīng)的多組分Euler方程具體過程，并驗證了采用無網(wǎng)格算法進(jìn)行化學(xué)非平衡流數(shù)值模擬的可行性。最后，對不同工況火箭燃?xì)馍淞髁鲌鲞M(jìn)行了數(shù)值研究。結(jié)果表明，無網(wǎng)格方法模擬火箭燃?xì)馍淞髁鲌鍪强尚械摹?/p>

1 數(shù)值方法

1．1 控制方程

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

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

這里僅列出主要的變量:

式中 ρ為密度;u、v分別為軸向、徑向速度;p為壓力;E 為單位體積的總能;ci(i=1，．．．，N－1)為 i組分質(zhì)量分?jǐn)?shù);ωi為i組分生成速率;N為總組分?jǐ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ù)逼近原理

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

式中 下標(biāo)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的守恒變量進(jìn)行重構(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ī)的雙基固體推進(jìn)劑產(chǎn)生的燃?xì)庵饕怯?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)表達(dá)式和系數(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)疊加到物理變量過渡值，最終得到下一時刻體現(xiàn)整體效應(yīng)的物理量值。其中，在計算化學(xué)反應(yīng)對流場的貢獻(xiàn)時，需把求解流動偏微分方程時采用的時間步長進(jìn)一步細(xì)分，作為求解化學(xué)反應(yīng)剛性常微分方程的步長。具體做法是先凍結(jié)化學(xué)反應(yīng)，求解得到流場參數(shù);然后，將化學(xué)反應(yīng)看做等容吸熱或放熱過程，保持內(nèi)能、速度參數(shù)不變，計算各組分的質(zhì)量變化率;最后，迭代求解溫度。

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

本文取燃?xì)馍淞鲌鲚S向計算長度為彈體直徑的40倍，徑向計算長度彈體直徑的10倍。火箭燃?xì)馍淞髁鲌鲇嬎銋^(qū)域中噴管出口附近空間布點如圖3所示，布點總數(shù)為28 049個。在本文計算中，不考慮火箭發(fā)動機(jī)燃燒室和噴管內(nèi)流動，直接在噴管出口給定燃?xì)獬曀偕淞鳁l件，并保持噴管出口燃?xì)馍淞鞯鸟R赫數(shù)Maj=2．5、靜溫Tj=1 500 K、燃?xì)饨M分與質(zhì)量分?jǐn)?shù)(如表2所示)不變。

表2 燃?xì)庵饕M分與質(zhì)量分?jǐn)?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翼型的跨聲速流場進(jìn)行計算。圖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所示。高速可燃?xì)怏w通過尖劈會發(fā)生偏轉(zhuǎn)，形成斜激波，波后混合氣體受斜激波壓縮的作用，壓力和溫度有所提高，但不能直接點燃混合氣體。壓縮后的混合氣體經(jīng)過一定的誘導(dǎo)期后形成斜爆轟波。整個流場是由斜激波、橫波、斜爆轟波、滑移線組成，這些波系相交于三波點。本文采用無網(wǎng)格算法的數(shù)值結(jié)果與文獻(xiàn)［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，火箭噴管出口的燃?xì)馍淞鲏毫Ψ謩e為 1．0、2．5、5．0 atm。當(dāng)噴管出口的燃?xì)馍淞鲏毫Υ笥谕饬鲏毫r為欠膨脹狀態(tài)，且隨著燃?xì)馍淞鲏毫Φ脑龃?，欠膨脹程度增加。不同噴管出口燃?xì)馍淞鲏毫ο氯細(xì)馍淞髁鲌龅鸟R赫數(shù)云圖如圖8所示，軸線馬赫數(shù)、壓力、溫度如圖9所示。

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

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

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

由圖8、圖9可看出，相交射流激波在軸線上的交點，隨著欠膨脹程度的增加而遠(yuǎn)離噴管出口。這是由于欠膨脹程度的增加，導(dǎo)致燃?xì)馍淞鲗ν饬鞯淖饔迷鰪?qiáng)，即相交射流激波強(qiáng)度增強(qiáng)，射流激波張角增大，從而交點更靠后。由流動參數(shù)軸線分布圖還可看出，隨著欠膨脹程度的增加，燃?xì)饬鞒鰢姽芎笈蛎洺潭纫泊蟠筇岣?，在相交射流激波前的馬赫數(shù)更高，壓力、溫度更低。

射流場軸線上離開噴管一定距離后，不同噴管出口燃?xì)馍淞鲏毫?yīng)下的軸線壓力也基本一致，且均約為外流的靜壓。這主要是由于燃?xì)馍淞髋c外流之間混合層的相互作用，通過剪切層的傳輸，導(dǎo)致燃?xì)馍淞鬏S線上的靜壓與外流的靜壓保持一致，而馬赫數(shù)、溫度仍保持較高水平?？梢?，火箭燃?xì)馍淞鲗ξ擦鲌龅挠绊懛秶浅４蟆?/p>

圖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個大氣壓，火箭噴管出口的燃?xì)馍淞鲏毫Ψ謩e為0．8、0．6、0．4 atm。當(dāng)噴管出口的燃?xì)馍淞鲏毫π∮谕饬鲏毫r屬于過膨脹狀態(tài)，且隨著燃?xì)馍淞鲏毫Φ慕档?，過膨脹程度將增加。不同噴管出口燃?xì)馍淞鲏毫ο律淞髁鲌龅鸟R赫數(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)下相交射流激波呈曲線。

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

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

計算工況為來流馬赫數(shù)分別為 1．5、2．0、2．5 Ma，來流靜壓為0．4 atm，火箭噴管出口的燃?xì)馍淞鲏毫? 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)格計算方法能有效地進(jìn)行化學(xué)非平衡流的數(shù)值模擬，能準(zhǔn)確反映流場的流動特性，具有布點靈活、求解效率高、精確模擬化學(xué)非平衡流中如激波等復(fù)雜流動現(xiàn)象的特點。采用無網(wǎng)格方法模擬火箭燃?xì)馍淞髁鲌?，為?shù)值研究燃?xì)馍淞髁鲌鲩_辟了一種新途徑。

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

(3)在欠膨脹或過膨脹狀態(tài)下，射流流場軸線上離開噴管一定距離后，不同噴管出口燃?xì)馍淞鲏毫?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］ 張福祥．火箭燃?xì)馍淞鲃恿W(xué)［M］．哈爾濱:哈爾濱工業(yè)大學(xué)出版社，2004．

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

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

［4］ 常見虎，周長省，李軍．高低空環(huán)境下火箭發(fā)動機(jī)射流場的數(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］ 周星，許厚謙．計算含運(yùn)動邊界非定常流動的無網(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．