雙縫激波誘導推力矢量噴管的非定常流動特性研究

(沈陽航空航天大學航空發動機學院遼寧省航空推進系統先進測試技術重點實驗室，沈陽 110136)

1 引言

未來飛行器對機動性和可操縱性的要求越來越高，而推力矢量技術可以有效提升發動機的機動性，近年來矢量噴管已成為學者們的研究熱點[1]。相比于傳統的機械矢量噴管，流體矢量噴管質量減少24%～80%，發動機推重比提高7%～12%，成本和維護費用降低37%～53%[2]。根據控制原理的不同，流體推力矢量技術有激波誘導法[3-4]、喉道偏斜法[5]、逆流控制法[6]、雙喉道法[7-8]等多種控制方法。相比于其他方法，激波誘導法動態響應更快且結構簡單[9]。

在激波誘導法的相關研究中，李麗等[10]采用試驗設計方法進行激波矢量控制噴管氣動性能數值分析與優化設計，綜合研究了二次流幾何參數對二元收擴噴管氣動性能的影響。屈高敏等[11]根據CFD計算結果，分析了NPR(落壓比)、主次流能量比以及次流速度對流體矢量裝置主流矢量角的影響。張相毅等[12]通過計算得到激波誘導矢量噴管的雙縫存在一個最佳間距和開縫位置，可以使噴管的矢量效率達到最大值。王曉明等[13]模擬了不同種類次流氣體噴射時激波誘導矢量噴管的推力矢量性能，得到噴管的推力矢量效率與次流氣體的質量分數成反比。王占學等[14]試驗研究了主次流氣動參數及幾何參數與矢量噴管流場結構和性能影響的關聯。張曉博等[15]研究得出，在滿足二次流壓力需求的前提下，宜采用風扇引氣作為二次流氣源。張少麗等[16]研究了主流落壓比、擴張段二次流落壓比、擴張段二次流角度和引射對固定幾何氣動矢量噴管軸向推力系數和矢量角的影響。Waithe等[17]針對二元收擴矢量噴管內部性能的研究表明，當NPR<4 時，在高次流壓力比條件下，多縫射流的推力矢量效果明顯好于單縫射流；當NPR>4時，多縫射流的推力矢量效果與單縫射流的相比并無過大差別。Giulinao 等[18-19]通過在噴管擴張段的上下壁面和側壁面進行射流控制來實現俯仰矢量控制和偏航矢量控制，其研究表明，噴管的俯仰矢量角可達到15°，偏航矢量角達到10°，但推力性能相比傳統噴管有所損失。

現如今對矢量噴管的研究主要集中于其推力矢量性能，對噴管的非定常流動特性研究較少，有待進一步研究分析。研究非定常流動需要先進的技術手段和方法，而傅里葉模態分解(FMD)技術是一種新的非定常流場分析方法[20]，可以滿足此方面的研究要求。本文數值模擬雙縫激波誘導推力矢量噴管的非定常流動特性，然后利用FMD技術對尾噴管的內外流場進行模態分解，得出非定常流場的頻率特性，獲得特定頻率對應模態的振蕩規律和潛在的動力學信息。

2 計算模型及方法

2.1 計算模型及網格

雙縫非軸對稱激波誘導矢量噴管內流場的模型和尺寸參考文獻[17]，外流場的模型和尺寸參考文獻[21]。考慮到外流場對計算結果可靠性的影響，外流場的高度和寬度分別為噴管出口寬度和高度的6倍和3倍，長度為噴管全長的10倍。噴管設計壓比為8.78，喉部面積為2.785×10-3m2，寬度為0.101 m，兩個二次流口尺寸均為0.101 m×0.001 m。計算域網格如圖1所示，圖中位于噴管入口中心處+位置為坐標原點。網格劃分參考文獻[21]。結構化網格應用于整個流場，總數約為290 萬。對兩個二次流口進行局部加密，外流場網格沿x 軸逐漸稀疏。近壁面第一層網格的y+<1。

圖1 計算域和網格Fig.1 Computational domain and grids

2.2 計算參數及設置

運用商業軟件FLUENT進行模擬計算，尾噴管落壓比為4.6，二次流壓比為0.7，出口靜壓為101 325 Pa，出口和入口總溫均為300 K，噴管側壁、上下壁面以及二次流口側壁皆為靜止無滑移壁面。采用分離渦模擬(DES)湍流模型求解流場[22-23]。時間步長為5×10-5s，共計算30個流動周期。

2.3 傅里葉模態分解技術

基于離散傅里葉變換(DFT)，FMD技術運用快速傅里葉變換(FFT)來快速、準確地提取復雜流場中的流動信息。與動力學模態分解(DMD)[24-25]和本征模態分解(POD)相比，FMD能大大節省計算時間，降低計算難度。FMD方法的核心內容見文獻[20]。

3 計算結果與分析

3.1 計算結果的準確性驗證

選擇激波誘導推力矢量噴管的上、下壁面中心線的無量綱壓力與文獻[17]中的試驗結果進行對比，如圖2所示。圖中，橫坐標是壁面監測點的軸向坐標x與噴嘴喉道坐標xt之比，縱坐標是監測點的壁壓p與噴管主流入口總壓pt,m的比值?？梢?，計算曲線與試驗曲線吻合良好，這說明本文后續的計算結果準確可靠。

圖2 計算結果與試驗結果的無量綱壓力對比Fig.2 Comparison between experimental results and calculation results of pressure distribution

3.2 噴管的流動結構

圖3示出了雙縫激波誘導推力矢量噴管對稱面的馬赫數等值線。欠膨脹狀態的二次流從噴管擴張段上壁面射入超聲速的主流，二次流的膨脹作用會使其在主流中形成一定的入射深度。二次流與主流的相互作用在噴口上游形成誘導激波(弓形激波)，而誘導激波的存在使得其下游形成強逆壓梯度，且逆壓梯度沿亞聲速邊界層向上游傳遞，進而導致誘導激波前流動分離，形成回流區?；亓鲄^的存在導致此位置的超聲速主流被壓縮，形成一道斜激波(分離激波)。誘導激波與分離激波相交，形成了λ型激波。

圖3 對稱面馬赫數分布Fig.3 Mach number contour at z=0

本文的矢量噴管擁有兩個二次流口，兩個噴口對應兩個λ型激波(圖3)。上游λ激波(即λ激波1)恰好延伸到下壁面唇口處并反射出一道激波，而反射激波延伸到噴管外流場的上剪切層的自由邊界上并反射為一道膨脹波。同時，反射激波與下游λ激波(即λ激波2)相交，使得此處激波分布呈現X型。λ激波2到達尾流下層剪切層的自由邊界上會反射出一道膨脹波，而膨脹波會繼續在自由邊界上反射為壓縮波，并繼續向下游交替變化運動。由于氣流的黏性，尾流與周圍環境氣體不斷進行動量、質量交換，使得這種交替壓縮膨脹現象逐漸消失。由圖2(a)可知，噴管上壁面存在兩個壓力突升處，對應于兩道分離激波。主流經過兩道激波后實現了兩次有效的矢量偏轉。上游二次流在主流的沖擊下附著于噴管上壁面，表現為閉式分離；下游二次流在噴口后并未附著于壁面，表現為開式分離[26]。

對圖4中6個y-z截面的馬赫數進行分析，結果如圖5所示。圖中，只截取了y-z截面中半個流場的數據，速度矢量為z和y方向速度矢量的合成。圖5(a)位于分離激波1 上游，流動受側壁影響產生的三維特性較弱。圖5(b)～圖5(f)中皆可以觀測到上壁面處的流動分離及激波。分離激波1延伸到下壁面唇口處并未造成流動分離，如圖5(f)所示。圖5(d)和圖5(e)中可同時觀測到兩道分離激波。激波在靠近側壁面位置顯示為偏轉狀態，是因為激波與側壁相交時側壁的邊界層發生了分離，進而生成了分離激波；兩道激波相交則表現為圖中的偏轉激波[26]。圖5(f)右上角可觀測到逆時針的流向渦，此流向渦的存在會強化尾流對周圍環境大氣的卷吸作用。

圖4 y-z截面的位置示意圖Fig.4 Specific position of the y-z plane

圖6 是尾噴管沿z軸方向不同截面的馬赫數分布?？梢?，不同位置截面中，二次流口前后分離區的尺度有所變化，這體現了流動的三維性，也印證了圖5(d)～圖5(f)中馬赫數的波浪狀分布。同時，不同截面的流動結構在靠近對稱面附近并沒有明顯變化，邊界效應引起的三維特性較弱，主要的流動結構發生在x-y截面上。因此，可選取噴管對稱面z=0來對流場主要流動結構做進一步分析。

圖5 噴管1/2 y-z面的馬赫數分布Fig.5 Mach number distribution of the nozzle 1/2 y-z plane

圖6 沿z軸不同截面的馬赫數分布Fig.6 Mach number distribution on different z-sections

圖7是矢量噴管下游二次流口處流線結構的局部放大圖。由圖可知，二次流口處的流動結構十分復雜，二次流口前的兩個渦分別是由誘導激波產生的逆時針分離渦，以及由于流動的引射作用生成在分離渦與二次流之間的順時針誘導渦。二次流口后是開式分離，可觀測到上壁面唇口下部的順時針繞流渦，以及由于繞流渦的特殊位置而在二次流與繞流渦之間生成的兩個誘導渦。大量研究[27-28]表明，噴管唇口附近渦的存在使得流場壓力波動變大，放大了流場的振蕩能力。

圖7 噴管對稱面下游二次流口處的流線分布Fig.7 Streamline distribution near downstream injection port at z=0

下面選取一個流動周期(T≈0.017 s)內的流線和壓力系數分布來詳細說明流場的流動結構變化。此處的流線是指由z=0 截面上x和y方向速度分量組成的矢量線。圖8示出了矢量噴管對稱面一系列瞬時流場的流線分布，從圖中可直觀看出剪切層把尾流和外界流體分成了速度相差極大的兩個區域。上剪切層的擺動程度強于下剪切層，下剪切層擺動較弱。這是因為上剪切層受到了上壁面唇口處誘導渦和繞流渦振蕩以及脫落渦的擾動影響，下壁面唇口處結構簡單，受到的擾動較小。T/6時刻，噴管上壁面唇口處出現渦脫落，生成的Vortex 1 運動至x=0.120 m 處；2T/6～3T/6時刻，Vortex 1 沿剪切層從x=0.160 m 位置運動至x=0.190 m 處，空間尺度增大；4T/6時刻，Vortex 1運動至x=0.220 m處，上壁面唇口處誘導渦產生的新的脫落渦Vortex 2 運動至x=0.120 m處；5T/6時刻，Vortex 1運動至x=0.250 m處且空間尺度增大，Vortex 2運動至x=0.160 m；T時刻，運動至x=0.270 m 處的Vortex 1 空間尺度繼續增大，Vortex 2 運動至x=0.170 m，此時唇口處產生新的脫落渦Vortex 3。T時刻流場分布與T/6時刻的類似，流場完成一次循環。上壁面唇口處的渦在向下游脫落過程中其空間尺度不斷變化。

圖8 噴管對稱面流線分布Fig.8 Streamline distribution at z=0

圖9 是z=0 截面上的壓力系數分布。所選用的壓力系數(其定義見式(1))與上文的流線分布對應的時刻相同?？梢姡瑖姽苁諗慷魏秃淼朗冀K保持高的壓力系數分布。分離激波前后壓差顯著，分離激波1 的位置和強度并無明顯變化。分離激波2 發生了輕微振蕩，這是由于下游二次流口前后渦系的振蕩和脫落導致分離區尺度改變，進而改變了分離激波2 的強度。分離區尺度和分離激波2 共同作用使尾流產生波動，尾流波動影響外流場的低壓區和高壓區，使得壓力區的壓力產生微小變化。剪切層由脫落的旋渦構成，因此可用脫落渦來表示剪切層處的壓力變化。

圖9 噴管對稱面壓力系數分布Fig.9 Pressure coefficient distribution at z=0

式中：pin為矢量噴管入口壓力，ρin、vin分別為噴管入口流體的密度和速度。

3.3 傅里葉模態分解功率譜分析

對計算數據的連續3 000個壓力系數流場進行FMD 處理，可得圖10 所示的功率譜圖。圖中，St為各階模態對應的斯特勞哈爾數，St=fD/Uref，Uref為入口流速且Uref=80 m/s，D為水力直徑且D=0.028 m；‖A‖為FMD處理后的壓力系數幅值。圖中存在多階模態，其中St=0.021對應的尖峰尤為突出，選取此模態為一階模態。二階模態選取為St=0.644，此模態對應的尖峰也較為突出。選擇St=0.644 為二階模態，是因為隨著St的增大，一階模態和二階模態之間的突出尖峰處的幅值是遞減的，而直到St=0.644，突出尖峰處的幅值增加。一階模態和二階模態分別對應各自的FMD流動結構振蕩模態，一階模態的幅值明顯大于二階模態的幅值，說明一階模態流動結構蘊含的能量更高、振蕩更強。

圖10 功率譜圖Fig.10 Power spectrum

圖11示出了以DEM(數字高程模型)形式顯示的噴管對稱面的一階模態和二階模態的幅值分布變化。對比兩階模態的幅值分布可知：兩階模態的幅值在λ激波2位置都較為突出，且沿著流動方向振蕩幅值都是越向下游越小。不同的是一階模態激波下游尾流處的振蕩幅值比激波處的小，二階模態激波下游尾流處的振蕩幅值在局部會大于其λ激波2 處的振蕩幅值，同時尾流振蕩的空間分布范圍很廣，λ激波2 處脈動能量只在局部有影響，較尾流擾動的空間范圍和整體擾動能量弱很多。綜上，一階模態的脈動能量主要源于下游二次流口的激波2和上剪切層的振蕩，二階模態的脈動能量主要源于尾流上下剪切層渦的脫落。

圖11 兩個模態DEM形式下的幅值云圖Fig.11 Amplitude contours from the DEM for the two modes

3.4 平均流模態

平均流模態代表原始流場的時均處理結果，是原始流場的基本流動結構[20]。圖12 示出了平均流模態的壓力系數分布，從圖中可以觀測到兩道分離激波和反射激波的位置，以及尾流處由剪切層反射的壓縮波和膨脹波分割而成的交替的高壓區和低壓區。

圖12 平均流模態的壓力系數分布Fig.12 Cp distribution of the mean flow mode

3.5 一階傅里葉模態分解模態

FMD模態代表不同頻率對應的流場壓力脈動，流場的一階FMD 模態在時域的發展如圖13 所示。上游二次流口激波1 處并未觀測到壓力脈動，只在上游二次流口后有小尺度的壓力振蕩；下游二次流口激波2處存在強烈的壓力振蕩，直觀印證了3.2節的分析。尾流上剪切層處的壓力振蕩也比較明顯，且各處的壓力脈動狀態并不相同。上剪切層中圓圈圈住部位(x=0.18 m)的壓力脈動正負與其余部位恰好相反，此處正好是反射激波與上剪切層相交之處，該現象應該是反射激波與上剪切層共同作用的結果。圖13(a)中上剪切層x=0.18 m 處的壓力為負脈動，其上游壓力為正脈動。圖13(b)中x=0.18 m 處的壓力仍為負脈動，其上下游剪切層壓力皆為正脈動。圖13(c)中x=0.18 m處的壓力即將由負脈動過渡到正脈動，趨近于平衡；其上游剪切層壓力恰好由正脈動過渡到負脈動，正處于平衡狀態(振幅趨近于0)，下游剪切層壓力的正脈動增強。圖13(d)中x=0.18 m 處的壓力過渡為正脈動，其上游剪切層壓力已過渡為負脈動，下游壓力振幅趨近于0。圖13(e)中，x=0.18 m 上游的剪切層壓力恰好由負脈動過渡到正脈動，處于平衡狀態，振幅趨近于0，下游剪切層壓力為負脈動。圖13(f)與圖13(a)的流場分布類似，流動完成1次循環。綜上，一階模態主要表現為下游二次流口激波處和尾流上剪切層處的壓力振蕩。

圖13 一階FMD模態的壓力系數分布Fig.13 Pressure coefficient contour of the first FMD mode

3.6 二階傅里葉模態分解模態

不同相位壓力系數分布的二階FMD 模態如圖14所示?？梢?，二階模態和一階模態之間存在很大差異。二階模態激波振蕩很弱，主要變化出現在尾流的上下剪切層中。上下剪切層都有渦不斷向下游脫落，且脫落過程中渦的空間尺度不斷發生變化，相鄰渦之間的壓力波動恰好相反，上剪切層中渦的空間尺度明顯比下剪切層渦的大。圖14(a)中，Vortex 1的渦核位于x=0.230 m處；圖14(b)～圖14(e)中，Vortex 1 的運動軌跡是x=0.250 m→x=0.280 m→x=0.300 m→x=0.320 m，在此過程中渦的空間尺度不斷增大；圖14(f)中，Vortex 1 脫落到x=0.350 m 處，Vortex 2運動至圖14(a)中Vortex 1的位置，圖14(f)與圖14(a)的流場結構類似，流動完成一次循環。二階模態主要表現為尾流上下剪切層處渦的脫落。

圖14 二階FMD模態的壓力系數分布Fig.14 Pressure coefficient contour of the second FMD mode

4 結論

采用DES方法數值模擬了雙縫激波誘導推力矢量噴管的非定常流場，研究了其原始流場的流動結構；利用FMD 方法對噴管對稱面流場進行了分解，分析了各階模態對應壓力系數分布。主要得到如下結論：

(1)噴管上壁面唇口處不斷有渦向下游脫落，尾流上剪切層的振蕩幅度強于下剪切層的振蕩幅度；流場在λ激波1處未發生振蕩，在λ激波2和尾流處發生了壓力脈動。

(2)噴管流場中存在兩個主FMD模態，一階模態主要表現為下游二次流口激波處和尾流上剪切層處的壓力振蕩，二階模態主要表現為上下剪切層處渦的脫落。