王建明，齊曉航，欒思琦，王成軍，陳寶嶺

(1.沈陽航空航天大學(xué)航空發(fā)動機學(xué)院,遼寧省航空推進(jìn)系統(tǒng)先進(jìn)測試技術(shù)重點實驗室,沈陽 110136;2.沈陽航天新樂有限責(zé)任公司,沈陽 110034)

超燃沖壓發(fā)動機作為高超聲速飛行器的動力裝置，相比于傳統(tǒng)的渦噴和渦扇發(fā)動機，它在高馬赫數(shù)下能夠擁有較大的推重比[1]。它的結(jié)構(gòu)簡單，由進(jìn)氣道、隔離段、燃燒室和尾噴管四部分構(gòu)成，構(gòu)型上采用飛行器機體一體化設(shè)計，將飛行器的下表面作為進(jìn)氣道和尾噴管的上壁面，以減少迎風(fēng)面積和重量。超燃沖壓發(fā)動機尾噴管是提供發(fā)動機推力的主要部件，它產(chǎn)生推力的原理是：燃燒室高溫高壓氣體在噴管入口處膨脹加速，利用尾噴管上壁面作為外膨脹面以增大燃?xì)獾呐蛎浢娣e[2]，讓出口處氣流和進(jìn)口處產(chǎn)生動量差，從而使發(fā)動機獲得凈推力。

影響超燃沖壓發(fā)動機尾噴管性能的因素眾多，對尾噴管性能影響規(guī)律的研究是近年來該領(lǐng)域的重要研究方向。Huang等[3]做了氫燃料超燃沖壓發(fā)動機噴管的化學(xué)非平衡流動分析，方法是利用CHEMKIN軟件對Hyshot型超燃沖壓發(fā)動機不同工況下噴管內(nèi)的化學(xué)非平衡現(xiàn)象進(jìn)行了數(shù)值分析。通過對凍結(jié)流、平衡流和非平衡流三種化學(xué)狀態(tài)的測試和比較，論證了化學(xué)反應(yīng)對噴嘴流場的影響。結(jié)果表明：真實非平衡流模擬計算結(jié)果介于凍結(jié)流和平衡流兩種極限情況之間，更接近于凍結(jié)流。在噴管入口燃?xì)馔耆紵那闆r下，自由基復(fù)合反應(yīng)會釋放出大量熱量，這些熱量主要用于提高氣體溫度，對推力作用很小。在噴管入口不完全燃燒情況下，噴管內(nèi)既有燃燒反應(yīng)，又有自由基復(fù)合，反應(yīng)放熱對推力的影響取決于噴管進(jìn)口處燃燒的完全程度，推力相比于凍結(jié)流則大幅度提高。Mo等[4]做了非均勻進(jìn)氣超燃沖壓發(fā)動機噴管設(shè)計與冷態(tài)流動的試驗。提出了一種考慮非均勻入流影響的非對稱超燃沖壓發(fā)動機噴管設(shè)計計算程序。給出了給定非均勻馬赫數(shù)分布和相應(yīng)的質(zhì)量加權(quán)平均均勻馬赫數(shù)分布的典型設(shè)計實例。然后，通過計算流體力學(xué)分析和冷態(tài)流動實驗測量，對非均勻設(shè)計的噴管性能改善進(jìn)行了量化。計算和實驗的結(jié)果表明，非均勻馬赫數(shù)噴管設(shè)計法設(shè)計出的噴管較均勻馬赫數(shù)方法下整體性能更優(yōu)。Lü等[5]做了考慮橫向展開和幾何約束的三維超燃沖壓發(fā)動機噴管設(shè)計的研究，提出了一種基于準(zhǔn)二維超音速流動和最大推力理論的考慮橫向展開和幾何約束的三維噴管設(shè)計方法。為了生成三維噴管的外形，采用特征線法計算無粘流場以及參考溫度法修正附面層厚度，并用計算流體力學(xué)方法計算了噴管的氣動性能。結(jié)果表明：初始圓弧半徑對軸向推力系數(shù)影響較小，而橫向輪廓、下整流罩長度和初始膨脹角的變化對軸向推力系數(shù)有顯著影響。并和用流線追蹤法設(shè)計的三維噴管進(jìn)行了對比研究，驗證了新方法的優(yōu)越性。黃偉等[6]也用數(shù)值模擬的方法論述了尾噴管構(gòu)型對性能的影響，研究表明對于俯仰力矩，上壁面切線角越小越好，外罩長度越長越好，而對于升力則剛好相反，所以在設(shè)計尾噴管構(gòu)型時要充分平衡升力和俯仰力矩的需求。Zhang等[7]做了非均勻入口馬赫數(shù)的超燃沖壓發(fā)動機噴管三維流場仿真實驗，分析某種構(gòu)型的噴管內(nèi)的馬赫數(shù)和壓力等值線，結(jié)果表明相比于均勻馬赫數(shù)入口條件，非均勻入口條件下的推力和陀螺力矩都會降低。

目前關(guān)于超燃沖壓發(fā)動機尾噴管流場的數(shù)值模擬大多都是基于二維流場，討論的內(nèi)容也是以入口參數(shù)對發(fā)動機推力、升力等性能參數(shù)的影響為主，很少有對噴管流場三維結(jié)構(gòu)的研究，而噴管流場的結(jié)構(gòu)特性也是影響噴管性能的重要因素?，F(xiàn)著重于研究噴管三維流場的結(jié)構(gòu)，通過觀察各個截面上密度、馬赫數(shù)等值線、渦量及速度矢量線圖來分析三維流場中的波系和渦系的分布特性。

1 計算模型與結(jié)果驗證

1.1 計算模型和網(wǎng)格

噴管計算模型選自文獻(xiàn)[8]，三維模型及結(jié)構(gòu)組成如圖1(a)所示，二維結(jié)構(gòu)示意圖如圖1(b)所示，主要設(shè)計參數(shù)及尺寸如表1所示，x是主流方向。

計算網(wǎng)格如圖2所示，圖2(a)是整個計算域的網(wǎng)格，圖2(b)是噴管處的網(wǎng)格。在噴管近壁面處的網(wǎng)格均進(jìn)行了加密處理，網(wǎng)格原點取在噴管入口的左下角。網(wǎng)格大部分由O型網(wǎng)格組成，局部位置采用非結(jié)構(gòu)網(wǎng)格，近壁面第一層網(wǎng)格高度設(shè)置為0.005 mm，y+＜1，網(wǎng)格總數(shù)約230萬。

圖1 噴管模型圖Fig.1 Model of nozzle

表1 噴管主要設(shè)計參數(shù)Table 1 Main design parameters of nozzle

1.2 數(shù)值方法

流動可視為理想氣體可壓縮定常流動。在三維坐標(biāo)系下微分形式的N-S方程如下：

(1)

(2)

(3)

式(1)～式(3)分別是質(zhì)量、動量和能量的守恒形式方程[9-10]。式中：t為時間；i、j代表方向；ρ為流體密度；p為靜壓；u為流體速度；E為流體微團總能量；T為溫度；gi和Fi分別為i方向上控制體所受的重力體積力和外部體積力；τij為應(yīng)力張量；k為有效熱傳導(dǎo)系數(shù)；Sh為其他熱源項。

圖2 計算域與網(wǎng)格Fig.2 Computing domain and grid

現(xiàn)選用基于SST(shear-stress transport)k-ω湍流模型。其中湍動能k的輸運方程為

(4)

特定湍流耗散率ω的輸運方程為

(5)

式中：Γk和Γω分別為k和ω的有效耗散系數(shù)；Gk和Gω分別為k和ω的產(chǎn)生項；Yk和Yω分別為k和ω由于湍流引起的耗散項；Xω是交叉擴散項。湍動能的耗散項Yk可表示為

Yk=ρβkωFDES

(6)

式(6)中：FDES為湍流強度項。

數(shù)值模擬用FLUENT軟件，求解器選擇壓力基，計算方法采用Coupled算法，選取Green-Gause節(jié)點格式，密度、壓力等參數(shù)均選取二階迎風(fēng)格式做定常計算，共迭代14 000步。設(shè)置流體為理想氣體，流體黏性采用Sutherland公式。所有壁面皆為絕熱、無滑移壁面。根據(jù)參考文獻(xiàn)設(shè)定邊界條件為噴管入口靜壓35 574 Pa，靜溫298 K，馬赫數(shù)2.5，自由流靜壓2 940 Pa，總溫673 K，馬赫數(shù)7.1。

1.3 網(wǎng)格穩(wěn)定性和計算結(jié)果驗證

為了保證計算精度的基礎(chǔ)上加快計算速度，對網(wǎng)格的獨立性進(jìn)行驗證，分別計算了180萬、200萬、230萬和270萬的網(wǎng)格，并對計算結(jié)果進(jìn)行無量綱化，計算結(jié)果如圖3所示，橫坐標(biāo)為噴管上壁面中心線位置坐標(biāo)，縱坐標(biāo)是噴管上壁面中心線上壓力無量綱化值，根據(jù)文獻(xiàn)[8],P是上壁面中心線的壓力值，Pc是試驗臺模擬的燃燒室出口氣體的壓力值，按照文獻(xiàn)[8]中無量綱化的方法取(P/Pc)以10為底的對數(shù)值。如圖3可知當(dāng)網(wǎng)格在230萬時結(jié)果已經(jīng)不再隨網(wǎng)格量的增加而改變，因此選取230萬的網(wǎng)格進(jìn)行計算。

仿真結(jié)果與文獻(xiàn)[8]的實驗結(jié)果的對比圖如圖4所示，圖4中的橫坐標(biāo)為噴管上壁面中心線坐標(biāo)，縱坐標(biāo)是噴管上壁面中心線處壓力的無量綱化的數(shù)值，圖4表明仿真的曲線結(jié)果與實驗測得的九個測量點的結(jié)果基本吻合，進(jìn)一步說明了仿真結(jié)果的準(zhǔn)確性和可信性。

圖3 網(wǎng)格獨立性驗證Fig.3 Grid independence verification

圖4 數(shù)值結(jié)果與實驗值對比Fig.4 Comparison of numerical and experimental values

2 三維流場結(jié)構(gòu)的結(jié)果與分析

2.1 三維流場中膨脹波分布

為了觀察三維流場結(jié)構(gòu)，除了截取噴管對稱面(z=0.025 m)以外，另取y=0.005 7 m、x=0.02 m、x=0.13 m三處截面，截面位置示意圖如圖5所示。

流場各個截面的密度分布圖如圖6所示。根據(jù)流場中密度分布情況可以確定膨脹波的位置。圖6(a)中，發(fā)現(xiàn)噴管內(nèi)存在沿流向密度下降明顯的邊界，流場中密度下降、壓力必然下降，因此該處即為管內(nèi)膨脹波的位置。圖6(a)中可以看到這兩道膨脹波的起始位置分別在下壁面噴管入口處和上壁面內(nèi)噴管出口處，這些位置產(chǎn)生膨脹波的原因是這兩處都是超聲速氣流發(fā)生外折的位置。根據(jù)圖6(a)中所示兩道膨脹波會在內(nèi)噴管出口處附近交匯，管內(nèi)氣流在經(jīng)過這兩道膨脹波時膨脹加速，進(jìn)而產(chǎn)生推力和升力。圖6(b)中可以看到圖6(a)中噴管內(nèi)膨脹波在y方向上的位置，同時可以觀察到產(chǎn)生于兩側(cè)壁面出口處的膨脹波，側(cè)壁面產(chǎn)生的兩道膨脹波均存在內(nèi)折角。根據(jù)圖6(c),可以得知膨脹波之所以發(fā)生內(nèi)折，是因為側(cè)壁尾端的射流邊界與側(cè)壁形成了外鈍角,而且射流邊界處的密度小于噴管出口的密度，即射流邊界處壓強小于入口處壓強，因此會在兩個側(cè)壁尾端產(chǎn)生內(nèi)折的膨脹波。

圖5 截面位置示意圖Fig.5 Schematic diagram of section position

2.2 三維流場激波和剪切層分布

圖7是流場中各個截面的馬赫數(shù)等值線，由于激波處馬赫數(shù)等值線會發(fā)生匯聚，因此根據(jù)馬赫數(shù)等值線匯聚的位置可以確定激波的大致位置。圖7(a)中可以觀察到上、下壁面出口處的羽流激波，羽流激波的存在是因為存在高空羽流效應(yīng)。高空羽流效應(yīng)的產(chǎn)生是因為超聲速噴流離開噴管時會在外部的超聲速介質(zhì)中流動，二者的氣流參數(shù)不同，此時噴管處的超聲速氣流在模擬狀態(tài)下的超聲速介質(zhì)中擴散，就會產(chǎn)生類似于發(fā)動機在高空飛行環(huán)境下所產(chǎn)生的高空羽流效應(yīng)，高空羽流效應(yīng)會導(dǎo)致飛行器產(chǎn)生附體斜激波，這類激波被稱為羽流激波[11-12]。在上壁面尾部，由于上方自由流穿過羽流激波，導(dǎo)致上壁面尾部附近的壓力明顯增高，進(jìn)而導(dǎo)致此處附近的自由流對主流產(chǎn)生壓縮，進(jìn)而形成一道由壓力決定的管內(nèi)斜激波。圖7(b)中，可以看到噴管側(cè)壁面同樣存在羽流激波。因此可以推斷，外噴管流場周圍的羽流激波是呈環(huán)狀分布在噴管流場周圍。并且根據(jù)圖中馬赫數(shù)等值線的分布情況可以看出隨著流動的發(fā)展，激波處馬赫數(shù)等值線變得越來越稀疏，這表明激波會隨著流動的發(fā)展強度逐漸減弱。

觀察馬赫數(shù)等值線圖，發(fā)現(xiàn)除了激波位置外，仍存在其他的等值線發(fā)生匯聚的位置。實際上流場中除了激波處，剪切層處的馬赫數(shù)等值線也會發(fā)生匯聚，因此需要判定這些位置上是否存在剪切層。剪切層最明顯的特征是它的內(nèi)部渦量很高，大多數(shù)研究者根據(jù)這個特征來判定流場中是否存在剪切層。渦量云圖如圖8所示，根據(jù)圖8(a)，發(fā)現(xiàn)了噴管上、下壁面尾端的剪切層包含z方向的渦。圖8(b)中可以觀察到噴管主流兩側(cè)剪切層則包含y方向的渦。

圖6 流場密度云圖Fig.6 Cloud map of density

圖7 流場馬赫數(shù)等值線圖Fig.7 Contour map of Mach number

2.3 三維流場渦流分布

流向的渦流比較難捕捉，為了確定流場中流向渦的具體位置，做出了各個截面上的速度矢量線圖，通過觀察各個截面的速度矢量線圖來確定渦流的位置，矢量圖中速度的大小值取基于該截面下的y方向和z方向的合成速度的值。根據(jù)圖9(a)，發(fā)現(xiàn)流場中有渦流的存在，由于該截面方向垂直于主流，因此這些渦流為流向渦，根據(jù)矢量線的分布，發(fā)現(xiàn)渦產(chǎn)生于內(nèi)噴管出口附近。渦流的形成原因是由于內(nèi)噴管出口處產(chǎn)生的高空羽流效應(yīng)，使得主流的超聲速氣流受到羽流的影響，氣流會沿橫向分離并發(fā)生卷折，進(jìn)而形成流向渦。圖9(b)中，上壁面出口處也觀察到了流向渦的存在。它的形成原因是由于上壁面的邊界層存在橫向分離的趨勢，這一趨勢導(dǎo)致氣流離開上壁面后，進(jìn)而在尾流中生成流向渦[13-14]。

圖8 流場渦量云圖Fig.8 Vorticity cloud map of flow field

2.4 三維流場主要結(jié)構(gòu)示意圖

圖10是綜合前面結(jié)果得到了三維流場結(jié)構(gòu)示意圖。圖10結(jié)果與文獻(xiàn)[8]的實驗結(jié)果相比，文獻(xiàn)[8]所得到的流場結(jié)構(gòu)示意圖并沒有討論流場中的膨脹波，討論的流場結(jié)構(gòu)也都是基于二維截面下的。本文揭示的流動結(jié)構(gòu)則是基于三維流場，比文獻(xiàn)[8]中的更全面，并且對所有物理現(xiàn)象的產(chǎn)生原因加以論述。

3 結(jié)論

通過數(shù)值模擬的方法對超燃沖壓發(fā)動機三維定常流場結(jié)構(gòu)進(jìn)行分析，結(jié)論如下。

圖9 速度矢量線圖Fig.9 Velocity vector line diagram

(1)噴管內(nèi)存在膨脹波，管內(nèi)的膨脹波會發(fā)生交匯形成膨脹波面。同時側(cè)壁也產(chǎn)生內(nèi)折的膨脹波，該處膨脹波會隨著氣流膨脹，強度逐漸變?nèi)酢?/p>

(2)噴管外流場周圍存在羽流激波和剪切層，羽流激波會形成一個激波面分布在外流場周圍，剪切層則分布在羽流激波和主流之間。在上壁面出口處存在羽流激波和因邊界層分離產(chǎn)生的斜激波，并且在這兩道激波之間也存在剪切層。

(3)內(nèi)噴管出口處氣流會在噴管主流兩側(cè)形成流向渦。外噴管出口處氣流也會在上壁面出口兩側(cè)形成渦流。內(nèi)噴管出口附近的流向渦主要是受到羽流效應(yīng)的影響，尾流中的流向渦則是由于上壁面氣流沿橫向分離造成的。