許 聰， 劉琪麟， 賴煥新

（華東理工大學(xué)機械與動力工程學(xué)院，上海 200237）

作為氣動聲學(xué)的基本問題，噴射流動及其噪聲的準(zhǔn)確預(yù)測一直是國內(nèi)外該領(lǐng)域的研究熱點。對于實際的航空發(fā)動機，尾氣噴流從渦輪到噴管出口的距離較短，而且在噴口處即為完全湍流狀態(tài)，但相關(guān)的邊界層厚度或湍流強度等湍流參數(shù)難以用實驗測量得知[1]。數(shù)值計算方面，由于計算資源的限制，目前還難以一體化地計算上游葉輪機械和下游噴流的流-聲場；噴射流動及其噪聲問題仍只能單獨考慮。而對于噴流模擬，上述湍流實驗參數(shù)的缺乏使得數(shù)值合成湍流來流缺少依據(jù)，Bogey 等[2]的研究表明，噴口初始狀態(tài)相關(guān)來流條件的準(zhǔn)確性嚴重影響遠場噪聲的預(yù)測結(jié)果。實際噴管出口往往是湍流狀態(tài)，而過去關(guān)于噴流計算的初始條件多為層流狀態(tài)。文獻[3-4]的結(jié)果顯示層流的初始剪切層可導(dǎo)致遠場總聲壓級的預(yù)測值與實驗值有5～10 dB 的誤差。層流的初始狀態(tài)導(dǎo)致噴流發(fā)展中存在層流向湍流轉(zhuǎn)捩，該現(xiàn)象會產(chǎn)生額外的低頻噪聲[5]，但實際噴流下游并不存在這種流動過程。因此，數(shù)值合成湍流的噴管出口，從而準(zhǔn)確預(yù)測湍流的混合層流動特性至今是噴流噪聲研究的重點[6]。文獻[7-8]將噴管包含在計算域內(nèi)，試圖讓其自由發(fā)展得出湍流來流，但商用飛機發(fā)動機的噴流和全尺寸實驗中的噴流中并沒有這樣的發(fā)展段[2]。因此，在噴流數(shù)值計算中如何利用人工擾動觸發(fā)轉(zhuǎn)捩，使其在盡可能短的范圍內(nèi)接近自由湍流條件，是當(dāng)前噴流模擬普遍采用的策略，也是實現(xiàn)噴流大渦模擬準(zhǔn)確預(yù)測的關(guān)鍵。

目前文獻中對于噴管流動常用的方法有兩種。一種是在流動中加入大尺度渦結(jié)構(gòu)相關(guān)特征量的合成湍流[9-10]，這一類方法可以使流動快速進入自相似區(qū)域進而改善預(yù)測結(jié)果，但是擾動作用隨著大尺度湍流結(jié)構(gòu)破碎而消失，不能準(zhǔn)確預(yù)測下游剪切層流動特性。另一種較普遍采用的是基于不穩(wěn)定本征模態(tài)的擾動方法[11-12]，該方法借助于噴流中不穩(wěn)定模態(tài)來觸發(fā)湍流。相對于前者，后者合成的來流更符合實際噴管的流動條件，使噴流下游的特性也與實驗一致，但是這種方法會增強近場湍流水平，進而導(dǎo)致遠場噪聲預(yù)測過高。Bogey 等[13-14]提出的多種周向模態(tài)組合的渦環(huán)擾動克服了這一問題，該方法對于主導(dǎo)流動的軸對稱模態(tài)不敏感[15]，并且可以降低模擬中預(yù)測過高的脈動和聲壓水平。但是受到擾動方式和施加位置的限制，只能在來流邊界的下游一段距離開始起作用，且在與實驗對比時，需要人為地減去渦環(huán)流向位置。

為了改善從噴口位置開始的全場預(yù)測結(jié)果，更早地觸發(fā)向湍流過渡，本文將Hu 等[11]基于線性不穩(wěn)定性分析得出的平行射流入口擾動擴展應(yīng)用到圓口噴流，同時結(jié)合周向模態(tài)快速觸發(fā)的優(yōu)勢，得出一種圓口噴流的多模態(tài)線性不穩(wěn)定性來流擾動方法，以克服渦環(huán)擾動的上述不足。在此基礎(chǔ)上，本文采用大渦模擬（LES）方法研究出口聲學(xué)馬赫數(shù)Ma=0.75、雷諾數(shù)ReD=8.7×105的噴流，將該擾動與Bogey 等[13]的渦環(huán)擾動以及來流無擾動的3 種工況進行對比分析，并采用已有的實驗結(jié)果進行驗證，研究比較3 種流動條件對該噴流流動和噪聲預(yù)測準(zhǔn)確性和有效性的影響。

1 計算設(shè)置

1.1 物理模型及網(wǎng)格劃分

本文計算的圓口噴管出口直徑Dj=50 mm[16]，基于噴口速度（Uj）及遠場聲速（α∞）的Ma=Uj/a=0.75 ，基于噴口直徑的ReD= ρUjDj/μ=8.7× 1 05，計算域示于圖1 中。其中， ρ 為密度， μ 為動力黏度。坐標(biāo)原點位于噴管出口的圓心處，計算域總長度為97.6Dj，其中噴管長度7.6Dj，漸擴段長度50Dj，緩沖區(qū)長度40Dj，計算域入口寬度為20Dj，出口寬度為40Dj。采用結(jié)構(gòu)化網(wǎng)格，對剪切層區(qū)域進行了網(wǎng)格加密，網(wǎng)格節(jié)點數(shù)為 3 .5×106。本文前期已經(jīng)通過外推法評估網(wǎng)格質(zhì)量，結(jié)果表明該套網(wǎng)格可以分辨出唇線上80%以上的湍流脈動動能，滿足LES 的質(zhì)量指標(biāo)要求[17]。

1.2 計算方法

大渦模擬計算采用商用軟件Fluent，計算時選用Density-Based 求解器，擴散項和對流項分別采用中心差分和高階迎風(fēng)格式，時間采用隱式推進。初場由雷諾時均模擬給出，后續(xù)LES計算采用Smagorinsky-Lilly 亞網(wǎng)格模型，量綱為一時間推進步長為0.005Dj/Uj，對應(yīng)全場CFL（Courant-Friedrichs-Lewy）≤1。經(jīng)過時長為300Dj/Uj的發(fā)展，認為漸擴段已排除初始流場的影響，之后開始進行湍流量的采樣統(tǒng)計，統(tǒng)計采樣時間持續(xù)150Dj/Uj。采用FW-H 方程計算遠聲場，Mendez 等[18]提出該方法要求聲源面需要包含絕大部分含能結(jié)構(gòu)，但又不能過大而影響截止頻率?；谕牧鹘鼒龅臏u量軸截面分布，取聲源面流向范圍25Dj，入口區(qū)域半徑2Dj，出口區(qū)域半徑4Dj。監(jiān)測點位于r=30Dj和r=50Dj的兩個同心圓面，圓心為噴管出口中心，流向監(jiān)測角θ取30°～150°，測點和聲源面位置如圖1 所示，內(nèi)部虛線為聲源面。

圖1 計算域及測點示意圖Fig. 1 Computational domain and measuring point positions

1.3 邊界與擾動條件

本文計算了3 種不同來流條件下的圓口噴流，邊界條件均為：入口邊界給定總溫320.4 K 和總壓147 059.9 Pa，出口為靜壓邊界，噴管壁面采用絕熱無滑移條件，其他邊界均給定遠場靜壓進口條件，環(huán)境壓力p∞=101 300 Pa，環(huán)境溫度T∞=288 K，對應(yīng)a∞=340 m/s。其中工況Nopert 的流場中無擾動施加是本文用作對照的基礎(chǔ)工況，更多計算細節(jié)前期工作[17]已給出，工況Pert1 中加入來自Bogey 等[13]的渦環(huán)擾動，工況Pert2 中加入改進的多模態(tài)線性不穩(wěn)定擾動。

工況Pert1 擾動計算式：

工況Pert2 的來流擾動即本文提出的以平行噴流擾動[11]為基礎(chǔ)改進多模態(tài)線性不穩(wěn)定擾動方案。Hu 等[11]基于線性不穩(wěn)定性分析，得出了一種平行射流入口的單一模態(tài)擾動方法，為使其可應(yīng)用于圓口噴流，本文考慮原擾動參數(shù)已經(jīng)優(yōu)化，因此首先對該擾動的流向分量保留，橫向分量則改造為圓口噴流的徑向。考慮多模態(tài)快速觸發(fā)湍流的優(yōu)勢，本文將原擾動方法的單一模態(tài)修正為周向多模態(tài)組合，得到一種改進的多模態(tài)線性不穩(wěn)定入口擾動，軸向和徑向擾動速度的計算式如下：

在式(3)和(4)中，a1=0.33,a2=0.4,b=6, ω =2Uj/Dj，t為時間，擾動施加位置為噴口處，其余參數(shù)均采用渦環(huán)擾動中的優(yōu)化數(shù)值[13]。采用軟件Fluent UDF 將兩種速度擾動施加在計算域中相應(yīng)位置。

2 流場特性

2.1 平均流動特性

為方便表述，本文做如下規(guī)定：u和u0分別表示軸向速度及其脈動，v和v0分別表示徑向速度及其脈動； h ·i 表示平均值，下標(biāo)t和θ分別表示時間和周向平均。圖2 和圖3 分別示出了3 種工況下軸向速度的流向和徑向分布，各計算工況與實驗中的勢流核長度如表1所示。由表1 可得，3 種工況的軸向速度預(yù)測結(jié)果與文獻[16]中實驗結(jié)果吻合良好，Pert1 的勢流核長度Lc/Dj為5.52，最接近實驗值6.5，Pert2 次之，Nopert 差距最大。此外，徑向分布中，兩種擾動工況比無擾動結(jié)果也有明顯改善。以上結(jié)果表明本文所采用的計算方法及邊界條件合理有效。

表1 勢流核長度Table 1 Lengths of potential core

圖2 中心線上軸向速度分布Fig. 2 Axial velocity distribution along centerline

圖3 x/Dj 為1、2.5、5 處軸向速度徑向分布Fig. 3 Radial distribution of axial velocity at position x/Dj=1, 2.5, 5

剪切層速度半值寬r0.5是衡量剪切層發(fā)展的一個重要指標(biāo)，其斜率反映剪切層的增長速率。從圖4剪切層半值寬的軸向分布中可以看出，Nopert 中存在斜率明顯不同的兩個階段，這種轉(zhuǎn)變是層流剪切層中渦配對引起的轉(zhuǎn)捩[19]。加入擾動后，兩擾動工況的剪切層增長速率在近噴口區(qū)域升高，而在勢流核下游降低，使剪切層半值寬整體上呈現(xiàn)無突變的緩慢增長，表明轉(zhuǎn)捩在噴口附近短距離內(nèi)已經(jīng)完成，符合自由湍流的發(fā)展機制，形成了接近真實噴流的湍流來流。

圖4 剪切層半值寬的軸向分布Fig. 4 Axial distribution of shear layer half-width

為了驗證該現(xiàn)象，圖5 示出了350Dj/Uj時刻的瞬態(tài)渦強的軸截面分布。從圖5 渦量的軸截面分布中可以看出，Nopert 的渦配對在流向3Dj～4Dj發(fā)生，圖中用虛線框標(biāo)識出了這一區(qū)域，該位置對應(yīng)剪切層半值寬斜率的轉(zhuǎn)變，標(biāo)志著層流向湍流過渡的開始。圖5(b)和圖5(c)中的Pert1 和Pert2 工況并沒有配對現(xiàn)象，在x/Dj＜3 之前，噴流剪切層已劇烈混合，處于完全發(fā)展的湍流狀態(tài)，與均衡發(fā)展的剪切層一致。圖6 示出了350Dj/Uj時刻x/Dj=2 截面瞬態(tài)渦強分布。從圖6 可以看出，渦量在橫截面x/Dj=2 的分布中，也可觀察到Nopert 的噴流仍處于勢流狀態(tài)，渦量呈現(xiàn)規(guī)則的環(huán)狀分布；而兩擾動工況的旋渦周向不穩(wěn)定性增強，這是由于加入的速度擾動破壞了流動穩(wěn)定性，使流動與周圍流體之間的摻混加劇。這些現(xiàn)象表明加入的擾動有效觸發(fā)了流動快速轉(zhuǎn)捩，形成了符合實際的湍流噴流。

圖5 t=350 Dj/Uj 時刻瞬態(tài)渦強的軸截面分布Fig. 5 Instantaneous vorticity magnitude distribution in axial section at t=350 Dj/Uj

圖6 t=350 Dj/Uj 時刻x/Dj=2 截面瞬態(tài)渦強分布Fig. 6 Instantaneous vorticity magnitude distribution in the plane x/Dj=2 at t=350 Dj/Uj

2.2 噴流渦結(jié)構(gòu)

噴射流動中，相干結(jié)構(gòu)主導(dǎo)動力學(xué)以及聲學(xué)特性，為了研究擾動對渦結(jié)構(gòu)演化過程的影響，采用Q判據(jù)等值面進行可視化。圖7 為Nopert、Pert1、Pert2 這3 種工況下Q判據(jù)等值面以及近噴口流域的局部放大圖，顏色由藍到紅表征為從0 到Uj進行著色。從圖7(a)示出的Nopert 工況中可以觀察到環(huán)狀渦形成、渦配對到失穩(wěn)破碎的完整級串過程，在流向x/Dj為1～4 范圍內(nèi)有軸對稱的渦環(huán)。而對于圖7(b)和7(c)中兩擾動工況，渦環(huán)的對稱特征已不再明顯，分別從Pert2 出口和Pert1 出口下游x≈Dj破碎成小尺度結(jié)構(gòu)。施加擾動的射流均識別出了更細致的渦結(jié)構(gòu)，物理現(xiàn)象更符合湍流來流的實際，再次證明了施加的兩種擾動的有效性。

圖7 t=350 Dj/Uj 時刻Q 判據(jù)等值面Fig. 7 Iso-surface of Q-criterion at t=350 Dj/Uj

2.3 湍流特性比較

2.3.1 空間分布 圖8(a)和8(b)分別給出了軸向速度脈動和徑向速度脈動沿中心線分布，其中rms 為均方根。3 種計算工況的速度脈動水平均高于實驗結(jié)果，且峰值位置向上游移動，這與中心線上快速的速度衰減而導(dǎo)致的黏性耗散有關(guān)。在x/Dj＜3 流域，Pert1 和Pert2 的速度脈動水平高于Nopert，這是由于Pert1、Pert2 和Nopert的剪切層在該區(qū)域處于不同的流動狀態(tài)，這與2.2 節(jié)的現(xiàn)象一致，并且Pert2 噴流向湍流過渡早于Pert1，這是由于擾動渦分別從不同的流向位置開始破碎。當(dāng)流動向下游發(fā)展，Pert1 和Pert2 的預(yù)測結(jié)果有不同程度的改善，均低于Nopert 預(yù)測的脈動水平，其中Pert1 的峰值水平最接近實驗值。從圖9 示出的脈動速度的徑向分布也可以看出，雖然在噴口附近本文計算的脈動峰值低于實驗測量值，但隨著流動向下游進行，兩種擾動工況的預(yù)測結(jié)果均有所改善，與測量的脈動水平相當(dāng)。相對于參考文獻[16]的LES 結(jié)果，本文兩種擾動的加入，使噴流湍流特性更接近實驗值。

圖8 中心線上脈動速度的均方根分布Fig. 8 Root mean square of fluctuating velocity distribution along centerline

圖9 脈動速度徑向分布Fig. 9 Radial profiles of fluctuating velocity

軸向速度脈動和徑向速度脈動沿唇線分布如圖10(a)和10(b)所示。對比3 種工況，在唇線分布上同樣可以觀察到，其中Pert1 工況峰值水平最低。在靠近噴口區(qū)域x/Dj＜2 ，Pert1 和Pert2 的速度脈動水平快速升高，且Pert2 脈動水平增長早于Pert1，而Nopert 在該區(qū)域發(fā)展緩慢，隨后在x/Dj＞2 的流域開始快速增長，Zhu 等[19]提出這是層流剪切層產(chǎn)生的自然現(xiàn)象。勢流核末端及其下游區(qū)域，擾動工況的徑向速度脈動低于無擾動工況，而軸向速度脈動區(qū)別較小。兩條流向線上的預(yù)測結(jié)果表明，施加的擾動能夠觸發(fā)出逼近實際噴嘴出口的湍流水平，并且降低下游的脈動水平，改善流動發(fā)展的計算準(zhǔn)確性。

圖10 唇線上脈動速度均方根分布Fig. 10 Root mean square of fluctuating velocity distribution along lip line

2.3.2 譜密度(PSD)分布 高頻聲源位于中心線速度開始衰減的位置，低頻聲源位于勢流核末端下游[20]，因此選取了唇線上x/Dj為2.5、5、10 這3 個流向位置，分析速度脈動譜密度PSD。圖11 和12 分別是唇線上軸向脈動和徑向脈動的PSD 圖。在流向位置x/Dj=2.5 處，Nopert 工況中軸向和徑向速度脈動的峰值量綱為一的頻率f D/Uj為0.6 左右，接近線性不穩(wěn)定性理論中f Dj/Uj=0.68 ，該頻率下的入口剪切層處于最不穩(wěn)定軸對稱模態(tài)，這代表了不穩(wěn)定波的發(fā)展[21]，造成與理論有差距的原因可能是非線性效應(yīng)。隨著流動向下游進行，脈動峰值逐漸向低頻移動，這與噴射流動中兩種不穩(wěn)定模態(tài)有關(guān)，剪切層模態(tài)主導(dǎo)噴口區(qū)域，由剪切層渦卷曲的高頻振蕩引起。下游則受噴流優(yōu)先模態(tài)控制，與低頻的大尺度湍流結(jié)構(gòu)的發(fā)展有關(guān)[22]。加入擾動后，在x/Dj=2.5 位置，Pert1 和Pert2 的PSD幅值在頻率f Dj/Uj≈0.6 處變化不大，但其他在頻段上均升高，因此與來流無擾動對比，擾動工況下的頻譜沒有明顯的峰值，全頻段的脈動水平都有所增強，說明Nopert 中的峰值可能與層流流動中的渦配對有關(guān)，而兩種擾動加入以后，快速轉(zhuǎn)捩得到的流動接近自然湍流，因此呈現(xiàn)出明顯的多尺度特征，頻譜表現(xiàn)為寬頻特性。在勢流核末端x/Dj=5 處，3 種工況下的PSD 趨于一致；而下游x/Dj=10的流向位置，Pert1 和Pert2 的PSD 在整個頻段上都低于Nopert，與2.3.1 節(jié)中脈動沿軸向分布結(jié)果一致，快速轉(zhuǎn)捩使得擾動工況的速度脈動峰值以及下游的脈動水平下降，頻譜分析也進一步表明擾動能促使噴流在近噴口區(qū)域快速轉(zhuǎn)捩以及在下游湍流發(fā)展。

圖11 唇線上軸向脈動的譜密度Fig. 11 PSD of axial velocity fluctuation along lip line

圖12 唇線上徑向脈動的譜密度Fig. 12 PSD of radial velocity fluctuation along lip line

3 遠聲場

3.1 總聲壓級（OASPL）

圖13 和圖14 分別示出了r=30Dj和r=50Dj兩圓面的OASPL 沿測角θ的分布。對比計算與實驗結(jié)果，可以看出Nopert 的OASPL 高于測量結(jié)果5～10 dB，加入擾動后，與實驗吻合程度明顯改善。表明無擾動工況渦配對導(dǎo)致轉(zhuǎn)捩的附加噪聲引起很大誤差，而采用擾動促成快速轉(zhuǎn)捩，可以有效地消除這種誤差。在r=30Dj的圓面上，Pert1 在θ為30°～110°的預(yù)測優(yōu)于Pert2，而在θ為110°～150°時，Pert2 更接近實驗值。而在更遠位置r=50Dj，兩種擾動工況預(yù)測相當(dāng)。因此兩種擾動雖有細微差別，但都能夠大幅度地改善無擾動對聲場計算造成的誤差。

圖13 遠場r=30Dj 測點總聲壓級Fig. 13 OASPL at far-field r=30Dj

圖14 遠場r=50Dj 測點總聲壓級Fig. 14 OASPL at far-field r=50Dj

3.2 頻譜特性

為了進一步研究擾動對聲場預(yù)測的影響，測角選取了60°、100°和140°這3 個角度進行頻譜特性分析。圖15 和圖16 分別示出了r=30Dj和r=50Dj兩圓面上的聲壓級（SPL）頻譜圖。與實驗結(jié)果對比發(fā)現(xiàn)，Nopert 的SPL 在低頻段（f Dj/Uj＜1 ）高于實驗結(jié)果，而高頻段（f Dj/Uj＞1 ）低于實驗結(jié)果，由于低頻大尺度湍流噪聲是總聲場主要貢獻者，其中就包括了渦配對帶來的附加噪聲誤差，因此在前面總聲壓級分布中，Nopert 的預(yù)測高于實驗值。相對于Nopert，隨著擾動的加入，Pert1 和Pert2 的高頻噪聲有所升高，但低頻噪聲的降低有效地使整體噪聲減弱，也證明了亞音速噴流中低頻噪聲占優(yōu)勢的地位。此外，Pert1和Pert2 的高頻SPL 與實驗仍有較大的視覺誤差，可能是計算網(wǎng)格對最小尺度湍流結(jié)構(gòu)的分辨仍然不是很充分而引起的。

圖15 遠場r=30Dj 聲壓級頻譜圖Fig. 15 Spectra of SPL at far-field r=30Dj

圖16 遠場r=50Dj 聲壓級頻譜圖Fig. 16 Spectra of SPL at far-field r=50Dj

4 結(jié) 論

(1)觸發(fā)噴口流動的快速轉(zhuǎn)捩，使之接近實際的湍流來流，是噴流大渦模擬的一個關(guān)鍵。本文改進的多模態(tài)線性不穩(wěn)定性擾動和渦環(huán)擾動都成功達到了這一效果。

(2)來流快速轉(zhuǎn)捩的實現(xiàn)，有效地消除了無擾動時渦配對造成的聲場誤差，改善了湍流噴流遠聲場的預(yù)測結(jié)果。

(3)本文改進的多模態(tài)線性不穩(wěn)定擾動方法，克服了渦環(huán)擾動湍流起始位置須人為移動的缺點。多模態(tài)的加入，更符合湍流來流的多尺度特性，體現(xiàn)了湍流頻譜的寬頻特征。