基于大渦模擬的熱噴流地面效應氣動特性研究

陳鵬飛，史亞林，楊 華，王 衡，劉志友

(1.中國航發四川燃氣渦輪研究院，四川 綿陽 621000；2.西北工業大學 動力與能源學院，西安 710129)

1 引言

近年來民航領域蓬勃發展，但由此帶來的噪聲嚴重影響了機場周圍居民的生活。當前高性能航空發動機工作時燃燒室的溫度可達2 000 K 左右，排氣溫度升高導致的噴流輻射聲強相關問題十分突出[1-4]。上世紀70 年代Hoch 發現，溫度升高使低速噴流的輻射聲強升高和使高速噴流的輻射聲強降低的現象。實際上飛機往往是在起飛和降落階段產生的噪聲最為明顯，在這些階段，噴管工作狀態大，飛機近地飛行，噴管排氣受到地面阻擋，為有限空間排氣，其流場更為復雜。目前，關于高溫下噴流近場的流動特性及近場與遠場的相關性的研究還不夠廣泛[5-8]，研究熱噴流地面效應氣動特性是理解溫度造成噴流輻射聲場特性變化的前提。

國外學者中，Hussian 等[9]研究發現，噴射流動中相干結構的發展由卷吸、配對及破碎等一系列過程組成。Thurow 等[10]觀察到，大尺度渦結構的尺寸相當或大于噴嘴直徑，大尺度渦結構是主要的流體動力學含能結構。國內學者中，李經警等[11-13]為揭示二元圓轉矩噴管尾噴流強化摻混的內在機制，應用大渦模擬(LES)方程對兩種相同進、出口直徑的噴管模型(軸對稱、二元圓轉矩)進行了數值模擬計算，結果表明：與軸對稱噴管相比，圓轉矩噴管射流摻混效應增強，速度衰減快，核心區長度和高溫區域面積減小，同時尾噴流擬序結構(包含渦環、渦辮、發卡渦、螺旋渦)發生變化。因而，從機理上弄清有限空間噴流的發展過程及摻混過程的影響因素極為必要，不僅可以為研究噴流地面效應機理提供依據，也能為建立近場輻射模型提供理論基礎。

本文主要采用數值計算的方法，建立發動機尾噴管排氣模型，對超聲速軸對稱噴管在地面影響下有限空間(Condition A，噴管中心離地面高度一般小于5 倍噴管直徑)和自由空間(condition B，噴管中心離地面高度一般大于5 倍噴管直徑)內噴流流動特性進行數值模擬研究。研究采用可壓縮流大渦模擬控制方程動態亞格子模型進行，對噴管尾噴流核心區長度與卷吸摻混面積變化、渦量和雷諾剪切應力分布，以及射流場擬序結構發展變化等方面進行了比較分析。

2 模擬算法及控制方程

數值模擬采用LES 方程的動態Smagorinsky模型，非穩態條件。利用濾波將瞬時變量φ(x,t)劃分為大尺度量(x,t) 和小尺度量φ′ (x,t)，(x,t) 通過以下加權積分得到：

式中：G(x-x′,Δ) 為濾波函數；Ω 為計算區域；Δ為濾波的寬度，與網格分辨率有關。

式(1)表示的濾波函數處理瞬時狀態下不可壓縮流N-S 方程時，有：

連續性方程

動量方程

能量方程

式中：μSGS為渦黏系數。

3 計算模型及網格

3.1 物理模型

為簡化計算工作量，研究模型為軸對稱拉法爾噴管，如圖1 所示。噴管出口直徑為Dj，噴管中心離地高度為H，坐標原點位于噴管進口圓心處。

圖1 噴管物理模型Fig.1 Physical model of nozzle

3.2 求解方法

計算域如圖2 所示。長度方向為 20Dj，寬度和高度方向均為 10Dj。流動方向為X軸正向，寬度方向為Y軸方向，高度方向為Z軸方向。其中Condition A 狀態H=2Dj，Condition B 狀態H=5Dj。

圖2 不同離地高度狀態下噴管計算域Fig.2 Calculation domain of nozzle at different ground clearance

具體邊界條件為：噴管進口為壓力進口，總溫為870 K，總壓為32 300 Pa。噴管壁面為絕熱無滑移條件。Far1、Far2、Far3、Far4 為壓力遠場。模型采用半模計算，對稱面邊界條件為Symmetry。Condition A 狀態時，Ground 為無滑移條件。Condition B 狀態時，Far5 為壓力遠場。聲學馬赫數基于噴口速度及遠場聲速，噴管出口馬赫數為1.526。

采用可壓縮流動基于密度的方法求解流動控制方程。初場由雷諾時均方法結合SSTk-ω湍流模型得出。經過時長為400 個時間步長的發展，認為流場已排除初始流場的影響，之后開始進行湍流量的采樣統計。計算總時間為0.02 s，取最終0.02 s 的計算結果進行流場分析。

3.3 網格劃分

為更好地捕捉射流與尾噴流相互作用區域的流動細節以及地面效應，在噴管壁面處和地面處對網格加密處理，根據Y+<1，得出第一層附面層厚度為0.01 mm。遠場網格隨著X值增大變大。經網格無關性檢驗和考慮模擬精確程度，最終網格數為600 萬。不同狀態下的網格劃分如圖3 所示。

圖3 不同狀態下網格劃分Fig.3 Grid generation in different states

4 計算結果及分析

4.1 溫度分布

軸向截面溫度分布如圖4 所示。中間溫度高的區域為噴流核心區，由于噴管出口為超聲速，導致其核心區不會與亞聲速噴流核心區一樣光滑。Condition A 的核心區長度大約為 11.8Dj，核心區在X/Dj=10 位置有個間斷面，在此位置地面開始對噴流產生較大影響，導致地面溫度升高。由于地面對噴流的影響，在噴管軸線以上的區域，噴流與環境摻混區域增大，摻混作用增強。Condition B的核心區長度大約為 10Dj，核心區在X/Dj=6 位置時有個間斷面。

圖4 軸向截面溫度分布Fig.4 Temperature distribution of axial section

4.2 速度分布

模型噴流中心線上無量綱速度(U/Uj，Uj為噴管出口截面平均速度)分布如圖5 所示。可見，在噴管內部，2 個狀態速度逐漸增加；在噴管出口初始段，存在1 段速度變化不大的區域(核心區)；在核心區內，隨著激波和膨脹波在自由邊界的交替發展，速度交替增大和減小，但總體上在一定范圍內波動；隨著射流與環境之間能量、動量交換增強，射流動能降低，速度逐漸減小。在X/Dj=9 之前，Condition A 和Condition B 噴流中心線上無量綱速度幾乎相同；在X/Dj=9 之后，由于地面對Condition A 噴流的影響，使得其噴流與空氣的摻混相較于Condition B加劇，速度整體來說降低得更快，如圖6 所示。

圖5 噴流中心線無量綱速度分布Fig.5 Non-dimensional velocity distribution of jet centerline

圖6 Condition A 對稱面馬赫數分布Fig.6 Mach number distribution of symmetry plane of condition A

4.3 摻混特性

采用大渦模型對射流摻混進行模擬，能夠描繪出其射流擬序結構隨時間和空間的變化，噴流剪切層的發展與旋渦運動密切相關。Q準則在渦識別中計算效率較高，效果明顯，是一種常用的渦量識別方法，本文采用Q準則作為旋渦可視化的工具。Q準則定義式為：

根據Q準則顯示出的尾噴流瞬時擬序結構如表1 所示，圖中顏色由速度場著色得到。可見，流場中渦旋結構主要由渦環、渦辮、發卡渦、螺旋渦組成，Condition A 和Condition B 射流場渦核發展過程大致相同。

表1 不同時刻渦結構Q 等值面對比Table 1 Comparison of Q equivalent surfaces of vortex structures at different times

以Condition B為例，當射流從噴管出口流出時，與環境氣流摻混作用較弱，但由于射流剪切層比較穩定，誘導出的渦結構較小，存在一段光滑區。在射流向前發展的過程中，氣流向內螺旋匯聚，破壞速度剪切層，形成渦環結構；其在剪切作用下會逐漸拉伸，而后脫落(脫落頻率受剪切層脈動影響)，隨后向內卷吸的氣流在射流剪切作用下，又形成新的渦環。同時，受渦環外側反向速度拉伸，使連接2 個渦環之間的渦管結構生長——這部分結構被稱為渦辮。該結構是在流動向下游發展的過程中，由于渦環之間的非線性不穩定作用增強，射流脈動較強的條件下逐漸形成的。隨著流動進一步發展，射流柱在剪切作用下脈動特征加強，渦環之間距離縮小，渦環與渦辮之間相互作用增強，最終形成尺度較大的發卡渦，使射流柱表面呈魚鱗狀；隨著射流繼續向下游發展，射流脈動減弱，大尺度發卡渦被耗散成更多小尺度的螺旋渦結構。

從Condition A 可以看出，當X/Dj=9 左右時，地面開始對噴流產生影響。對比t=0.010 2 s 和0.020 0 s 時刻，Condition A 比Condition B 在軸向更早出現發卡渦。Condition A 時，在地面有限空間內的湍流脈動更加強烈，噴流在地面的作用下誘導出二次渦流，渦核分布更加集中，且其中以大尺度發卡渦為主。從Condition B 在t=0.020 0 s 時刻的渦核結構可以完整地看到渦環產生-脫落-產生，渦環反向運動將渦管拉伸為渦辮，隨著渦環距離減小，發卡渦形成，但最終隨著湍流脈動的耗散，發卡渦破碎為螺旋渦。

4.4 雷諾應力分布

雷諾應力能表征湍流脈動強弱，因此通過比較2 個狀態雷諾剪切應力特征，可以反映其射流剪切層內流體脈動狀況，進一步反映地面對噴流摻混特性的區別。水平面無量綱雷諾剪切應力為，垂直面無量綱雷諾剪切應力為。其中，u為X向脈動速度，v為Y向脈動速度，ω為Z向脈動速度。

圖7(a)、圖7(b) 分別給出了Condition A 和Condition B 在各截面水平中心線的雷諾應力分布。可見，兩者雷諾應力分布趨勢類似。在不同截面，雷諾應力從中心軸線隨著Y值的增大先增大后減小，在中心線處值為0，極值只有1 個，分布曲線呈現單峰狀態。

圖7 各工況截面中心線雷諾應力分布Fig.7 Reynolds stress distribution at the center line of section under various working conditions

在Condition B 狀態下，隨著X/Dj增大，雷諾應力的極大值先增大后減小，極大值對應的位置向中心線方向移動。雷諾應力分布趨勢的內在機理為：以X/Dj=5 為例，中心線(Y/Dj=0)上的雷諾應力最小，隨著Y值增大而逐漸增大，在剪切層邊界附近達到最大。在這個區域附近，射流與外界氣流發生劇烈摻混，能量耗散加劇，而后剪切應力逐漸減小。隨著X/Dj增大，射流雷諾應力的極大值先增大后減小。X/Dj=3 時，該截面距離噴口很近，射流處于光滑段，雷諾應力很小，符合前文提到的光滑段擬序結構少的特征。X/Dj=5 時，雷諾應力達到最大，此時射流處于過渡段，射流脈動最強，也佐證了射流在過渡段表面旋渦尺度增大、主要成魚鱗狀的特征。X/Dj=12 以后，雷諾應力隨著大尺度渦被耗散，射流脈動減小，也驗證了此區域以較小尺度的螺旋渦結構為主的擬序特征。在Condition A 狀態下，X/Dj=9 時的雷諾應力極大值大于X/Dj=5 時的極大值。其內在機理為：X/Dj=9 時，地面對噴流的影響較大，射流脈動在有限空間內被擠壓，地面對射流產生了擾動。

圖7(b)、圖7(c)分別為Condition B 各截面水平中心線雷諾應力和Condition A 各截面垂直中心線雷諾應力分布圖。Condition B 為自由射流，各截面水平和垂直中心線雷諾應力保持一致。Condition A 各截面垂直中心線雷諾應力分布呈現雙峰狀態。Z/Dj<0 時，兩者雷諾應力分布趨勢一致，但Condition A 各個截面垂直中心線位置處的雷諾應力比Condition B 的大，其分布曲線更加飽滿。Z/Dj<0 時，Condition A 各個截面垂直中心線位置處的雷諾應力比Condition B 的大，且各截面雷諾應力在Z向范圍更大。計算結果表明，由于地面有限空間對噴流的影響，使得湍流脈動更強，并且使得湍流脈動被擠壓至噴管中心線以上的區域，這也佐證了X/Dj=12 時垂直面速度分布特征。

5 試驗驗證

通過發動機熱噴流場整機試驗，測得發動機離地2Dj高時噴管排氣地面X/Dj=12 處地表的最大溫升約為30 K，其排氣流場瞬時紅外熱像圖如圖8所示。可以看出，試驗和仿真結果一致，充分表明了本文計算模型的有效性。

圖8 熱噴流場瞬時紅外熱像圖Fig.8 Instantaneous infrared thermal image of hot jet field

6 結論

采用大渦模擬計算研究地面熱噴流氣動特性，建立不同離地高度的超聲速噴管模型，對總溫為870 K、馬赫數為1.526 的熱噴流狀態開展了數值模擬，分析了有限空間噴流的發展過程。主要結論如下：

(1) 自由空間噴流的核心區長度大約為10Dj，在X/Dj=6 位置時有個間斷面。發動機有限空間噴流受地面阻擋，排氣能量在地面一側堆積，能量釋放區域變長，核心區長度增大為11.8Dj，核心區間斷面推后到X/Dj=10 位置，在此位置地面開始對噴流產生較大影響，導致地面溫度升高；在X/Dj=9 之前，有限空間噴流與自由空間噴流的無量綱速度幾乎相同，在X/Dj=9 之后，由于地面對噴流的影響，使得噴流與空氣的摻混加劇，其速度整體來說降低得更快。

(2) 有限空間與自由空間噴流場渦核發展過程大致相同。在流場中，渦旋結構主要由渦環、渦辮、發卡渦、螺旋渦組成；在X/Dj=9 左右位置時，地面開始對噴流產生影響，在地面有限空間內的湍流脈動更加強烈；噴流在地面的作用下誘導出二次渦流，其渦核分布更加集中，且其中以大尺度發卡渦為主。

(3) 通過雷諾應力和試驗驗證表明了計算模型的有效性。噴管排氣地面X/Dj=12 處地表最大溫升約為30 K，為研究熱噴流聲場測試提供了參考。