王建明， 欒思琦， 齊曉航， 毛曉東， 王成軍

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

推力矢量噴管分為機械推力矢量噴管[1-2]和流體推力矢量噴管。流體推力矢量噴管是指噴管在非機械偏轉(zhuǎn)狀態(tài)下，通過注入二次流來達到控制主流方向進而實現(xiàn)氣流的矢量偏轉(zhuǎn)[3]。根據(jù)原理不同，可分為多種控制方式：激波矢量控制[4]、雙喉道控制[5]、逆向流法[6]、同向流法[7]等。激波矢量控制技術(shù)由于結(jié)構(gòu)簡單、矢量效率高已成為航空屆的重點研究方向[8-9]。激波誘導矢量噴管的工作機理是在噴管擴張段引入二次流，二次流和主流相互作用生成誘導激波，誘導激波的存在導致了激波后的強逆壓梯度的生成，逆壓梯度沿亞聲速邊界層傳遞至激波前，導致誘導激波前回流區(qū)的生成，進而生成了分離激波，兩道激波相交稱為“λ”激波，主流經(jīng)過激波后實現(xiàn)了矢量偏轉(zhuǎn)。

李志杰等[10]對激波控制的二元收擴流體推力矢量噴管進行了數(shù)值模擬，研究自由流馬赫數(shù)(Ma)，噴管落壓比(nozzle pressure ratio, NPR)，主次流總壓比(secondary pressure ratio, SPR)相同時，不同二次流噴射位置對噴管矢量角及其他噴管性能參數(shù)的影響，研究指出二次流噴射位置對推力矢量角有較大影響，但置對推力系數(shù)的影響不明顯，落壓比是影響推力系數(shù)變化的主要因素。徐學邈等[11]研究了9組不同射流角下的二元矢量噴管流場，得出射流垂直壁面注入時產(chǎn)生的推力矢量角最大，逆流注入氣流對噴管流場的改變要明顯大于順流注入氣流的影響, 但產(chǎn)生的損失也較大。吳曉軍等[12]采用數(shù)值模擬手段對流體二次引射推力轉(zhuǎn)向參數(shù)影響規(guī)律進行了研究，得出不同參數(shù)對干擾流場結(jié)構(gòu)和推力轉(zhuǎn)向偏角的影響，最終得出流體二次引射實現(xiàn)推力轉(zhuǎn)向的基本設(shè)計原則及較優(yōu)的參數(shù)組合方案。王猛杰等[13]通過研究噴管落壓比影響矢量角機理及基于多縫腔體和多縫輔助注氣方法的分離區(qū)控制，得知矢量性能隨落壓比增大而下降主要是由次流下游分離區(qū)的結(jié)構(gòu)變化導致的，并最終得到環(huán)境壓力的多縫輔助注氣可以在不增加系統(tǒng)額外引氣量的條件下使得大落壓比時的矢量性能得到提高。宋亞飛等[14]對基于激波控制的二維射流推力矢量噴管非定常流場進行研究，分析二次流總壓連續(xù)變化對噴管流場和推力性能的影響，結(jié)果表明噴管擴張段注入二次流的總壓增加到一定程度才可以實現(xiàn)推力矢量控制。史經(jīng)緯等[15]考慮到變比熱合及溫度對傳熱系數(shù)、黏性系數(shù)的影響，對次噴噴口幾何構(gòu)型在不同噴管NPR、SPR工況下的二維流動特性及推力矢量特性進行了研究，指出在大的二次流壓比下，誘導激波會與噴射壁面相交，導致推力矢量效率的惡化。Song等[16]提出了一種新型的混合SVC噴管，在傳統(tǒng)SVC噴管的基礎(chǔ)上將超聲速氣流通過噴管的上下壁面橫射入超聲速主流，并對噴管的基本流動機理、氣動參數(shù)和幾何參數(shù)對其性能的影響進行了計算研究。Neely等[17]采用小尺度超聲速噴管研究射流推力矢量控制技術(shù)，測量了一定NPR(10～30)、SPR(0.012～0.11)范圍內(nèi)噴管的紋影可視化和雙組分推力。在二次流的噴射壓力足以在二次流噴口對側(cè)壁的唇口處造成誘導激波的情況下，獲得約5°的最大矢量角。He等[18]對圓-矩形噴管進行三維數(shù)值模擬，研究了不同SPR下二次噴射角和槽長對噴管流動特性和推力性能的影響。Forghany等[19]研究了不同自由流馬赫數(shù)(Ma=0.05～1.1)對收擴噴管推力矢量控制的影響。研究指出增加自由流馬赫數(shù)會降低推力矢量控制效果，優(yōu)化射流角度可以減少外部自由流對推力矢量控制性能的負面影響，在不同自由流工況下增大SPR和減小NPR，都可以縮短啟動射流推力矢量過程的響應時間。

目前的推力矢量噴管研究大多針對噴管的定常性能參數(shù)，而關(guān)于研究噴管的非定常流動結(jié)構(gòu)的報道較少。為此，采用數(shù)值模擬方法求解2DCD矢量噴管的非定常流場，并通過對噴管對稱面上的流線、渦量系數(shù)、密度梯度和熵的分布等的具體分析，得出噴管流場的非定常流動特性。

1 計算模型和數(shù)值方法

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

激波誘導矢量噴管模型選自文獻[20]。噴管主要設(shè)計參數(shù)如表1所示。

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

噴管計算域和全局網(wǎng)格如圖1(a)所示。二次流噴口和噴管壁面處進行了加密處理，如圖1(b)所示。全場采用六面體結(jié)構(gòu)化網(wǎng)格，總數(shù)約250×104。近壁面第一層網(wǎng)格高度為0.01 mm，y+<1。

圖1 計算域和網(wǎng)格Fig.1 Computational domain and grids

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

數(shù)值模擬利用FLUENT軟件，采用分離渦(detached eddy simulation, DES)[21]進行湍流模擬。DES湍流模型是RANS(reynolds averaged Navier-Stokes)和LES(large eddy simulation)的混合方法，是模擬復雜工程問題中分離流動的一種主要方法[22]。選用基于SST(shear-stress transport)k-ω的DDES(delayed detached eddy simulation)湍流模型。其中，湍動能k的輸運方程：

(1)

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

(2)

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

Yk=ρβ*kωFDES

(3)

引入DDES概念，通過因子FDES對湍動能消耗項加以修改，F(xiàn)DES表示如式(4)所示：

(4)

式(4)中：CDES為模型常數(shù)；Δ為單元網(wǎng)格的最大尺度；Lt為湍流長度尺度；FSST為混合函數(shù)。

速度和壓力的耦合采用Coupled算法，物理時間步長 Δt=5×10-5s，共計算150個流動周期。計算所用流體為理想氣體，流體黏性采用Sutherland公式。所有壁面皆為絕熱、無滑移壁面。噴管入口總壓為466 095 Pa，二次流入口總壓為326 266.5 Pa。

2 計算結(jié)果和分析

2.1 計算結(jié)果驗證

網(wǎng)格數(shù)量直接影響著數(shù)值模擬的計算精度和收斂速度。為了在保證計算精度的基礎(chǔ)上盡可能地加快計算速度，對網(wǎng)格進行獨立性驗證。如圖2所示，分別計算了150×104、200×104、250×104和300×104的網(wǎng)格，橫軸為噴管的擴張段坐標，縱軸為擴張段上壁面的壓力分布。從圖2中可知，當網(wǎng)格數(shù)達到250×104時，壁面壓力基本不再變化，就此可認定網(wǎng)格數(shù)量≥250×104后，計算結(jié)果不再發(fā)生變化。以此同時計算速度也較為合理。因此選擇250×104網(wǎng)格進行計算。

圖2 網(wǎng)格獨立性驗證Fig.2 Verification of grids independence

圖3為計算所得噴管的上、下壁面中心線無量綱壓力分布與文獻[20]實驗結(jié)果的對比。其中，Xt為噴管喉道的軸向坐標，Pt,m為噴管主流的入口總壓。圖3說明計算曲線和實驗所得曲線的結(jié)果基本一致，這進一步證明計算結(jié)果具有一定準確性和可信性。

圖3 數(shù)值計算所得壁面無量綱壓力與實驗[20]對比Fig.3 Comparison of dimensionless wall pressure between numerical and experimental[20]results

2.2 噴管的流場結(jié)構(gòu)

圖4為噴管流場X方向和Z方向系列截面的壓力系數(shù)分布。壓力系數(shù)定義為

(5)

式(5)中：Pin為噴管入口壓力；ρin為主流入口流體密度，vin為主流入口的流速。

從圖4(a)可看到：噴管喉道和收斂段處壓力系數(shù)顯示為紅色，為高壓分布。超聲速主流進入擴張段發(fā)生膨脹加速，壓力逐漸減小。由于二次流的作用主流被壓縮，壓力增大。圖4(a)可觀測到分離激波的位置。圖4(b)中，距離二次流較遠的上游流場并未受到激波的影響，流場壓力系數(shù)分布保持上下對稱性；受到激波影響區(qū)域以X=0.108 m的截面為例，此截面貫穿了二次流噴口前的分離區(qū)，截面的壓力系數(shù)從上到下越來越小，可分為三個區(qū)域，這三個區(qū)域的生成原因分別為：由二次流引起的流動分離生成的高壓區(qū)、流經(jīng)分離激波后增壓減速的流動以及未流經(jīng)分離激波的主流。未流經(jīng)分離激波的主流相較其他區(qū)域仍然保持低壓高速狀態(tài)。沿流向，未流經(jīng)分離激波的主流區(qū)域越來越小。從圖4可知噴管的對稱面處結(jié)構(gòu)是接近對稱的，因此可以用對稱面的流動結(jié)構(gòu)揭示三維流場，從而更直觀地分析流場的結(jié)構(gòu)變化。

運用Q準則[23-24]來清晰識別噴管內(nèi)流場的旋渦結(jié)構(gòu)。Q準則是速遞梯度張量的二次不變量[25]，速度梯度可表示流體運動中的角變形率，Q可定義為

(6)

式(6)中：Ωij為速度梯度的反對稱張量；Sij為速度梯度的對稱張量，分別定義為

(7)

(8)

圖4 噴管內(nèi)流場壓力系數(shù)分布Fig.4 CP distribution of flow field in the nozzle

圖5是噴管流場的Q準則等值面圖。在圖5中可觀測到二次流噴口前后分離區(qū)中渦系的存在。圖6為二次流入口附近流線分布。需要說明的是，三維流線一般不會處于同一個二維平面上，但如前所述由于流場結(jié)構(gòu)的對稱性，因此用噴管對稱面上的X-Y速度矢量線來近似代替流線。對比圖5、圖6可知，二次流噴口前渦的生成機理是誘導激波的存在生成的分離渦，以及由于流動引射生成的誘導渦，二次流噴口后的主流在高壓二次流的作用下不再附著于上壁面，因此外流場氣流繞過噴管上壁面唇口進入內(nèi)流場，進而生成了繞流渦，以及在繞流渦和主流間生成的誘導渦。

圖5 噴管內(nèi)流場旋渦結(jié)構(gòu)Fig.5 Vortex structure of flow field in the nozzle

圖6 二次流入口附近流線分布Fig.6 Streamlines distribution near the secondary flow inlet

圖7 FMD幅頻曲線Fig.7 FMD amplitude-frequency curve

接下來分別從一個周期中噴管的流線、渦量系數(shù)、密度梯度和熵的變化分析噴管的流動結(jié)構(gòu)。周期取決于對噴管流場做傅里葉模態(tài)分析(Fourier mode decomposition, FMD)得出的一階頻率，為512.8 Hz(圖7)；FMD的具體內(nèi)容見文獻[26]。流線、渦量系數(shù)、密度梯度和熵增所選用的圖片是由相同時刻的數(shù)據(jù)后處理得到的。

圖8為一個周期內(nèi)噴管對稱面上的流線分布。流線不同顏色代表了X、Y方向合成速度的大小。由于尾流和外部氣流的速度差和流體黏性的影響，尾流上、下側(cè)會生成剪切層。T/6時刻，從上壁面唇口處脫落的上層渦1運動至X=0.12 m處，從下壁面唇口處脫落的下層渦1位于X=0.155 m處；2T/6 時刻，上層渦1繼續(xù)向下游脫落，已運動至X=0.125 m處，下層渦1運動至X=0.17 m處；(3/6- 4/6)T時刻，上層渦1從X=0.14 m運動至X=0.155 m處，下層渦1從X=0.175 m運動至X=0.19 m處，同時可觀測到下壁面唇口處新生成的下層渦2已脫落至X=0.13 m處；5T/6時刻，上層渦1運動至X=0.16 m處，下層渦1運動至X=0.195 m處，下層渦2運動至X=0.135 m處；T時刻，上壁面唇口處新生成的上層渦2即將脫落至T/6 時刻上層渦1的位置，下層渦2即將脫落至T/6 時刻下層渦1的位置。T/6時刻和T時刻流場結(jié)構(gòu)類似。通過流線分布可知上、下壁面的唇口處不斷有渦沿著剪切層外側(cè)向下游脫落，脫落的過程中伴隨著旋渦空間尺度的變化。

渦量系數(shù)定義為

(9)

式(9)中：Vvor為渦量；L為噴管特征長度。

圖9為一個周期內(nèi)噴管對稱面渦量系數(shù)的變化。二次流噴口前側(cè)渦量的波動是由于分離區(qū)的分離渦和誘導渦的振蕩所致。尾流的上剪切層處渦量系數(shù)顯示為紅色， 說明上剪切層內(nèi)的旋渦為逆時針旋轉(zhuǎn)的渦，下剪切層處渦量系數(shù)為負，下剪切層內(nèi)的旋渦為順時針旋轉(zhuǎn)的渦。由于上壁面唇口處繞流渦和誘導渦的振蕩等因素，上剪切層的擺動角度明顯強于下剪切層。上剪切層的擺動是尾流振蕩的主要原因之一。

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

圖9 噴管對稱面渦量系數(shù)分布Fig.9 Vorticity coefficient distribution at Z=0

圖10、圖11分別為噴管對稱面一個周期的X方向和Y方向密度梯度分布。由于下壁面唇口處回流區(qū)的存在，使得此處剪切層向內(nèi)凹陷，主流在λ波前被壓縮，生成了一道反射激波，反射激波與λ波在靠近下壁面處相交。λ波與尾流的自由邊界相交，反射出膨脹波系。在圖10、圖11中皆可觀測到λ波、反射激波和膨脹波的具體位置。激波隨著流動的發(fā)展是不斷發(fā)生振蕩的λ波的振蕩是由于二次流噴口前的分離區(qū)中渦的不斷振蕩，改變了分離區(qū)尺度引起的，反射激波的振蕩是因為下壁面唇口處渦系的振蕩和脫落引起的。圖10的喉道處的密度梯度為負，清晰地觀測到主流在經(jīng)過喉道處時發(fā)生的膨脹。λ波作為強激波在圖10的X方向的密度梯度分布和圖11的Y方向的密度梯度分布皆為紅色，顯示為密度增大，流體在此處被壓縮。反射激波在圖10的X方向密度梯度分布為紅色，顯示為壓縮狀態(tài)；在圖11的Y方向密度梯度分布為藍色，在Y方向發(fā)生膨脹；反射激波作為一道壓縮波，X方向密度梯度占主導位置。膨脹波與反射激波恰好相反，在圖10的X方向密度梯度分布為藍色，在圖11的Y方向密度梯度分布為紅色，X方向密度梯度占主導位置。圖11的Y方向的密度梯度可觀測到尾流處剪切層的密度梯度，上剪切層為負，下剪切層為正，說明上剪切層處于膨脹狀態(tài)，下剪切層被壓縮；圖10的X方向尾流剪切層的密度梯度趨近于0。

圖10 噴管對稱面X方向密度梯度分布Fig.10 The X-direction density gradient distribution

圖11 噴管對稱面Y方向密度梯度分布Fig.11 The Y-direction density gradient distribution

圖12為噴管對稱面一個周期相對于主流入口的熵增分布，可觀測到λ波和反射波的存在。λ波的上游流場熵分布均勻，主流經(jīng)過λ波后實現(xiàn)了熵增，同時熵分布不再均勻，忽略剪切層部分，主流熵增為50～110 J/(mol·K)。經(jīng)過反射激波后主流并無明顯熵變化。噴管內(nèi)流場中分離區(qū)的熵增最大，因為分離區(qū)存在渦的振蕩導致流場混亂度增高。隨著流動時間的發(fā)展，尾流在發(fā)生能量耗散的同時伴隨著熵增，如圖12中圈示部分。外流場中，噴管尾流熵的最大值位于剪切層處，剪切層由內(nèi)到外熵越來越大，熵值跨度約為210 J/(mol·K)。剪切層熵增是由于剪切層的擺動增大了其湍流強度，流場的混亂增強；同時，剪切層中的渦不斷進行著脫落和振蕩，也是熵增的原因之一。

圖12 噴管對稱面熵增分布Fig.12 Entropy increase distribution

3 結(jié)論

運用數(shù)值模擬方法研究流體推力矢量噴管的非定常氣動特性。通過對矢量噴管的流線、渦量系數(shù)、密度梯度和熵的分布的分析，得到如下結(jié)論。

(1)隨著流動的發(fā)展，上下壁面唇口處生成的渦沿著剪切層外側(cè)不斷向下游脫落。

(2)二次流噴口前分離區(qū)渦系的振蕩改變了分離區(qū)的尺度，導致了λ波隨著流動的發(fā)展不斷發(fā)生振蕩；下壁面分離區(qū)渦系的振蕩和脫落導致了反射激波的振蕩。

(3)λ波是強激波。尾流上剪切層為膨脹狀態(tài)，下剪切層處于壓縮狀態(tài)；上剪切層的擺動角度和熵增皆大于下剪切層。