亞音速圓口噴流的大渦模擬

林 健, 賴煥新

(華東理工大學(xué)承壓系統(tǒng)與安全教育部重點實驗室,上海 200237)

?

亞音速圓口噴流的大渦模擬

林 健, 賴煥新

(華東理工大學(xué)承壓系統(tǒng)與安全教育部重點實驗室,上海 200237)

采用大渦模擬方法研究了馬赫數(shù)M=0.75,雷諾數(shù)Re=50 000的圓口可壓縮噴流問題。在與實驗數(shù)據(jù)比對驗證的基礎(chǔ)上,分析了噴流剪切層不穩(wěn)定性發(fā)展的兩個主要階段,著重考察了剪切層中軸向速度脈動u′、徑向速度脈動v′、壓力脈動p′在空間任意兩點上的時空相關(guān)系數(shù)隨時間和空間發(fā)展分布的特點,并由此分析剪切層中旋渦脈動的遷移特性,為進(jìn)一步研究與控制噴流噪聲提供基礎(chǔ)。

亞音速; 噴管噴流; 大渦模擬; 剪切層; 相關(guān)性分析

亞音速噴流噪聲是現(xiàn)代民航客機(jī)中主要的噪聲來源之一。隨著世界各國對民航客機(jī)噪聲適航標(biāo)準(zhǔn)與要求的日益苛刻,進(jìn)一步深入研究噴流流場中的渦運動與發(fā)聲過程,從而為研究和改善大型民航客機(jī)噪聲水平提供可能,依然是當(dāng)今民用航空噪聲控制的必由之路。

自從Lighthill在1952年創(chuàng)立氣動聲學(xué)[1-2]這一學(xué)科以來,噴流噪聲一直是氣動聲學(xué)的基本問題,聲源機(jī)理的解析是問題的關(guān)鍵。Lilley[3]總結(jié)了自Lighthill至上世紀(jì)末共50多年中噴流噪聲預(yù)測的Lighthill方程及各種版本的聲比擬理論方程,指出這些方程中的源項本身都是未知量,必須結(jié)合聲源流場的準(zhǔn)確預(yù)測。Tam[4]指出具有物理意義的聲源是流場中的旋渦,而聲比擬中的所謂四極子、偶極子等僅是用于數(shù)學(xué)上重構(gòu)某一給定聲場的等價聲源,它們并不具有物理意義。Tam[4-5]在總結(jié)分析了大量實驗數(shù)據(jù)之后,進(jìn)一步提出噴流中存在兩類噪聲源,前者是大尺度渦結(jié)構(gòu),后者是湍流的小尺度脈動,它們各自可用一種相似頻譜來描述。但由于實驗測量的困難,噴流流場中的渦運動迄今并未被完全清楚地認(rèn)識。例如對于噴管出口下游處極薄的剪切層尚無有效的測量方法,而該區(qū)域的流動狀態(tài)對于噴流噪聲又非常關(guān)鍵。

本文選取Jordan[6]和Andersson[7]研究的亞音速噴管噴流為對象,采用大渦模擬(LES)方法進(jìn)行計算分析,在方法驗證的基礎(chǔ)上,著重考察分析剪切層中不穩(wěn)定性發(fā)展過程。在此基礎(chǔ)上,分析了兩點兩時刻流場中脈動速度和脈動壓力的相關(guān)系數(shù),討論它們隨時間和空間發(fā)展分布的特點,為下一步分析亞音速噴流的噪聲場提供基礎(chǔ)。

1 亞音速噴管噴流計算

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

本文計算的噴管幾何尺寸與流動參數(shù)來源于文獻(xiàn)[6-7]。其中噴管出口直徑Dj為50 mm,出口馬赫數(shù)為0.75。計算域如圖1所示,噴管外長度取為90Dj,上游寬度為20Dj,下游寬度為40Dj。整個計算區(qū)域均采用結(jié)構(gòu)化網(wǎng)格,節(jié)點總數(shù)為1 640 854。對剪切層區(qū)域進(jìn)行加密處理,噴管內(nèi)部采用O型網(wǎng)格劃分,如圖1所示。

1.2 控制方程與邊界條件

為對噴管中的可壓縮流動進(jìn)行大渦模擬,采用Favre[8]過濾得到的LES控制方程組:

連續(xù)性方程:

(1)

(2)

(3)

圖1 計算域與網(wǎng)格劃分Fig.1 Computational domain and mesh

(4)

(5)

(6)

為使方程組封閉補(bǔ)充狀態(tài)方程:

(7)

式中R=287.1 J/(kg·K)。

本文選用動態(tài)Smagorinsky-Lilly模型,計算亞格子尺度應(yīng)力采用式(8):

(8)

(9)

(10)

本文先采用RNGk-ε湍流模型進(jìn)行穩(wěn)態(tài)計算,當(dāng)整體質(zhì)量流量誤差小于1%時認(rèn)為達(dá)到收斂,將所得結(jié)果作為LES的初始流場。計算時選用Density-Based求解器,流體介質(zhì)設(shè)置為理想氣體,黏度采用Sutherland模型。瞬態(tài)時間步長設(shè)置為5×10-7s,確保計算域內(nèi)90%區(qū)域迭代步長CFL數(shù)小于0.5。本文從初場開始,計算時長約130(Dj/Uj),此時非穩(wěn)態(tài)流動沿時間已充分發(fā)展,可排除初場的影響,在此基礎(chǔ)上再存取數(shù)據(jù)進(jìn)行統(tǒng)計分析。

上游噴管入口設(shè)置:給定總壓147 152Pa,總溫320.4K,同時給出參考靜壓101 325Pa; 噴管外計算域入口以及四周環(huán)境界面設(shè)置為無反射的壓力邊界條件,給定總壓101 325Pa,總溫288K和參考靜壓101 325Pa;下游出口設(shè)置為無發(fā)射的壓力出口邊界條件,給定靜壓101 325Pa,總溫288K; 噴管固壁采用無滑移、絕熱壁面邊界條件。

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

2.1 數(shù)值計算與實驗結(jié)果的比較

為方便表述,對本文的符號作出如下說明:〈…〉t表示時間平均; 〈…〉θ表示軸向平均; Uj為噴管出口速度; u和u′分別表示軸向速度和軸向速度脈動,v和v′分別表示徑向速度和徑向速度脈動; uc為中心線軸向速度; 坐標(biāo)原點為噴嘴出口圓心。

圖2示出了本文計算所得噴管中心線時均軸向速度,圖中給出了計算結(jié)果與文獻(xiàn)[6]實驗值、文獻(xiàn)[7]計算值的分布對比。中心線速度〈uc〉t≥0.95Uj時對應(yīng)的距離為勢流核長度。從中可以看出,本文結(jié)果與文獻(xiàn)[7]計算結(jié)果基本一致,預(yù)測勢流核長度均比實驗值偏短,這主要與亞格子模型的選取和網(wǎng)格質(zhì)量有關(guān)。

圖3示出了在噴口下游x/Dj為1、2.5和5這3個位置處軸向速度沿徑向分布對比,可以看出在x/Dj為1和2.5這兩處,本文計算結(jié)果更接近實驗值; 但在x/Dj=5處文獻(xiàn)[7]計算值稍好,可能是由于本文采用的網(wǎng)格在下游位置比文獻(xiàn)[7]的網(wǎng)格較少導(dǎo)致的。

圖2 噴流中心線上軸向速度的分布Fig.2 Axial velocity along the centerline

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

圖4給出了不同軸向位置上湍流脈動的徑向分布,可以看出當(dāng)前結(jié)果和實驗值基本吻合,整體上與文獻(xiàn)[7]的計算結(jié)果一致,表明本文的計算方法在模擬亞音速噴流時是可行的。

2.2 斷面流速分布的相似性

圖5示出了在噴管自由噴流軸截面上,各個不同斷面〈u〉tθ/Uj的徑向分布。從圖中可以看出,不同斷面軸向流速的分布顯示出一定程度的相似性,即軸線上流速最大,往下游發(fā)展流速逐漸減小,距離噴管出口越遠(yuǎn)流速衰減越慢。為進(jìn)一步分析斷面流速分布的相似特性,定義剪切層厚度δ[10]為

(11)

式(11)中r0.10、r0.95以及下文中的r0.5分別表示不同斷面上軸向流速為10%〈uc〉t、95%〈uc〉t和50%〈uc〉t(即中心線軸向流速的10%、95%和50%)處對應(yīng)的半徑值。

圖6中離散點為本文剪切層厚度的計算結(jié)果。由于數(shù)據(jù)量太大,本文僅選取了7個位置的數(shù)據(jù)結(jié)果。使用最小二乘法對這些離散點進(jìn)行線性擬合,延長該擬合直線使之與x軸相交于x0,即為噴流虛源。

將流速〈u〉tθ和斷面上徑向坐標(biāo)r分別以量綱為一坐標(biāo)〈u〉tθ/〈uc〉t和(r-r0.5)/(x-x0)來表示,則所有斷面上量綱為—流速均落在一定的區(qū)域之內(nèi),如圖7所示?？梢钥闯霰疚挠嬎憬Y(jié)果基本上已經(jīng)體現(xiàn)了自由噴流斷面流速分布相似性的這個重要特性。但可以看出圖中量綱為—流速一致性還有待改善,這主要是受到數(shù)據(jù)樣本容量的限制以及虛源定位方法準(zhǔn)確性限制[11]所致。

圖4 不同軸向位置上湍流脈動的徑向分布Fig.4 Radial profiles of turbulent fluctuations at axial positions x/Dj=1,2.5 and 5

圖5 不同軸向斷面上〈u〉tθ/Uj的徑向分布Fig.5 Radial profiles of 〈u〉tθ/Uj at different axial position

圖6 剪切層厚度δ/Dj沿軸向坐標(biāo)的分布Fig.6 Shear layer thickness δ/Dj with axial coordinate

圖7 〈u〉tθ/〈uc〉t隨(r-r0.5)/(x-x0)分布Fig.7 Variation of 〈u〉tθ/〈uc〉twith (r-r0.5)/(x-x0)

2.3 剪切層不穩(wěn)定性的發(fā)展過程

剪切層不穩(wěn)定性的發(fā)展過程直接影響著噴流流場的湍流特性,對于亞音速噴流噪聲機(jī)理的研究具有重要的意義。許常悅等[12]通過分析渦量厚度δΩ沿剪切層的分布來對該過程進(jìn)行了分析。Wan等[13]用噴管半值速度時的半徑r0.5和渦量厚度δΩ對自由流噴管和旋流式噴管進(jìn)行了對比研究。兩者的物理意義基本相似,均表示剪切層增長率。

圖8示出了噴管r0.5沿軸向坐標(biāo)的分布(rj為噴管出口半徑),圖中曲線主要呈現(xiàn)出兩個斜率:第1階段(515)斜率快速增長。為進(jìn)一步分析該現(xiàn)象,圖9示出了兩個時刻(t1=223Dj/Uj和t2=224Dj/Uj)渦強(qiáng)Ω0的瞬時分布,Ω0=‖×V‖。

圖8 噴管半值速度寬r0.5沿軸向坐標(biāo)的分布Fig.8 Half-width of jet r0.5 along x

從圖9可以看出當(dāng)515時,剪切層開始混合、擾動渦合并,導(dǎo)致了曲線第2個階段斜率快速增長的變化。

2.4 時空關(guān)聯(lián)分析

遠(yuǎn)場噪聲的強(qiáng)度及能譜與聲源中任意兩點兩個時刻的時空相關(guān)系數(shù)有密切關(guān)系[14]。如圖10所示,對于剪切層中軸向速度的脈動值,點x與其下游點x+ξ處的兩點時空相關(guān)系數(shù)Ruu定義如下:

(12)

式中:x為測點在流場中的軸向坐標(biāo);ξ為兩個測點間的距離;τ為延遲時間。時空相關(guān)系數(shù)越大,表明兩個現(xiàn)象相關(guān)的程度越高。當(dāng)兩個測點在湍流的同一個旋渦內(nèi)時,湍流脈動速度的相關(guān)程度就高,對應(yīng)相關(guān)系數(shù)就越接近于1。反之,兩個脈動值的相關(guān)程度就低,時空相關(guān)系數(shù)會趨于零[15]。

圖11示出了剪切層中r=0.5Dj,x=2.5Dj處的點與其下游x+ξ處的點兩點間的時空相關(guān)系數(shù)分布。當(dāng)ξ=0時,對應(yīng)圖中相關(guān)系數(shù)從最大值1起始的一條曲線,該曲線表示x=2.5Dj位置上t時刻與t+τ時刻的兩個脈動的相關(guān)系數(shù)隨τ逐步變化的關(guān)系。圖中其他每一條連續(xù)曲線均代表固定點x=2.5Dj處與x+ξ(ξ以0.1Dj遞增)處的時空相關(guān)系數(shù)?？梢钥闯鰰r空相關(guān)系數(shù)的峰值隨著ξ的增加而下降,表明兩點之間脈動的相關(guān)程度隨著空間距離的增大而降低。當(dāng)空間距離增大時,旋渦在向下游發(fā)展的過程中不斷地卷入周圍的流體介質(zhì),互相之間發(fā)生質(zhì)量、動量和能量的交換,并伴隨著能量的耗散,因此兩個脈動現(xiàn)象之間的相關(guān)關(guān)系逐步被解除。

圖9 渦強(qiáng)的瞬時分布Fig.9 Instantaneous vortical shown by Ω0

圖10 剪切層中脈動相關(guān)系數(shù)定義的對應(yīng)位置x和測點間距離ξ的示意圖

Fig.10 Two-point space-time correlations in the shear layer are obtained for locationxusing a spatial separationξ

圖12示出了噴流剪切層中r=0.5Dj,x=2.5Dj處的兩點時空相關(guān)系數(shù)Ruu隨ξ和τ的二維分布圖。圖中等值線表示Ruu,點劃線是Ruu極大值的連線,其斜率表征了旋渦的遷移速度UA。

圖11 剪切層中r=0.5Dj,x=2.5Dj處的軸向 速度脈動時空相關(guān)系數(shù)Fig.11 Two-point space-time correlation of axial velocity fluctuation at a location x in the shear layer r=0.5Dj,x=2.5Dj

圖12 剪切層中r=0.5Dj,x=2.5Dj處Ruu云圖Fig.12 Contours of Ruu obtained in the shear layer ofr=0.5Dj,x=2.5Dj

由于旋渦在遷移過程中發(fā)生質(zhì)量、動量和能量的變化,因此理論上只有在原點處的斜率才等于當(dāng)?shù)氐男郎u遷移速度UA[16]:

(13)

式中τRmax表示兩點時空相關(guān)系數(shù)取極值時的τ值。本文在圖11中處理時空相關(guān)系數(shù)時,ξ的增量取值為0.1Dj,因此近似采用Δξ=0.1Dj估算噴管唇線上x=2.5Dj、5Dj及10Dj3個點上的旋渦遷移速度,結(jié)果列于表1。

由表1可看出本文計算結(jié)果與文獻(xiàn)[7]計算結(jié)果基本一致。距離噴管出口x/Dj<5內(nèi),UA趨于恒定,為0.65Uj。隨著流動向下游發(fā)展和剪切層的增長變厚,使得旋渦的遷移速度降低,在x/Dj=10時UA減小為0.53Uj。計算結(jié)果略低于文獻(xiàn)[7]的0.56Uj,可能是本文的計算網(wǎng)格較少造成的。

表1 剪切層噴嘴唇線(r=0.5Dj)上各點的遷移速度UATable 1 Convection velocity UA at different positions in the shear layer r=0.5Dj

在可壓縮流動中,聲波的傳遞與壓力的脈動值p′密切相關(guān),為此本文進(jìn)一步對壓力脈動p′、徑向流速的脈動v′、以及軸向速度脈動u′之間的時空相關(guān)進(jìn)行了分析。參照Ruu的定義式(12),定義以下時空相關(guān)系數(shù)Rvv、Rpp、Rup、Rvp分別如下:

圖14示出了剪切層中r=0.5Dj,x/Dj為2.5、3和5這3處Rup與Rvp的對比?？梢钥闯鯮up與Rvp基本同相位,但各峰值點處Rvp數(shù)值均高于Rup,顯示了p′與v′之間存在較強(qiáng)的相關(guān)關(guān)系。

由于Lighthill聲源的主要貢獻(xiàn)來自軸向速度,壓力項對總聲源貢獻(xiàn)較小[14]。所以v′增大雖然會伴隨出現(xiàn)p′增大,但并不會對總聲源造成明顯的增強(qiáng),因此在考慮降噪措施時,如果某一種方法能適當(dāng)?shù)卦龃髒′而減小u′,則可減弱主要噪聲源的貢獻(xiàn),該方法將有可能是一種可行的降噪措施。

圖13 剪切層中r=0.5Dj,x=2.5Dj處Rvv、Rpp云圖Fig.13 Contours of Rvv、Rpp obtained in the shear layer (r=0.5D,x=2.5Dj)

圖14 剪切層中不同位置處壓強(qiáng)脈動與速度脈動相關(guān)系數(shù)Fig.14 Two-point space-time correlation of pressure and velocity fluctuation obtained in the shear layer

3 結(jié) 論

本文采用大渦模擬的方法對馬赫數(shù)為0.75的亞音速圓口噴管噴流進(jìn)行研究,得到如下結(jié)論:

(1) 亞音速噴流剪切層不穩(wěn)定性的發(fā)展過程主要分為兩個階段:第1階段擾動渦開始失穩(wěn)并進(jìn)行配對,剪切層厚度緩慢增加; 第2階段剪切層開始混合、擾動渦合并,剪切層厚度快速增長。

(2) 剪切層中軸向速度脈動u′、徑向速度脈動v′和壓力脈動p′的時空相關(guān)分布規(guī)律類似,而p′的遷移速度比u′和v′的遷移速度要快。

(3) 剪切層中Rup與Rvp基本同相位,p′與v′之間存在較強(qiáng)的相關(guān)關(guān)系。如果某一種方法能適當(dāng)?shù)卦龃髒′而減小u′,則該方法將有可能是一種可行的降噪措施。

[1] LIGHTHILL M J.On sound generated aerodyne-mically:I.General theory[J].Proceedings of the Royal Society,A:Mathematical,Physical and En-gineering Sciences,1952,211(1107):564-587.

[2] LIGHTHILL M J.On sound generated aerodyne-mically:II.Turbulence as a source of sound[J].Proceedings of the Royal Society：A.Mathematical,Physical and Engineering Sciences,1954,222(222):1-32.

[3] LILLEY G M.Jet noise and classical theory and experiments[C]//Hubbard H H.Aeroacoustics of Flight Vehicles:Theory and Practice.Hampton,Virginia:NASA,1991:211-289.

[4] TAM C K W,VISWANATHAN K,AHUJA K K,etal.The sources of jet noise:Experimental evidence[J].Journal of Fluid Mechanics,2008,615:253-292.[5] TAM C K W.Jet noise:Since 1952[J].Theoretical and Computational Fluid Dynamics,1998,10(14):393-405.

[6] JORDAN P,GERVAIS Y,VALIERE J C,etal.Final results from single point measurements[R].[s.l.]:[s.n.],2002.

[7] ANDERSSON N,ERIKSSON L E,DAVIDSON L.Largeeddy simulation of subsonic turbulent jets and their radiated sound[J].AIAA Journal,2005,43(9):1899-1912.

[8] VREMAN B.Direct and largeeddy simulation of the compressible turbulent mixing layer[D].The Netherlands:University of Twente,1995.

[9] LARCHEVEQUE L,SAGAUT P,MARY I,etal.Large-eddy simulation of a compressible flow past a deep cavity[J].Physics of Fluids,2003,15(1):193-210.

[10] ARAKERI V H,KROTHAPALLI A,SIDDAVARAM V,etal.On the use of microjets to suppress turbulence in a mach 0.9 axi-symmetric jet[J].Journal of Fluid Mechanics,2003,490:75-98.

[11] Tide P S,Babu V.Numerical predictions of noise due to subsonic jets from nozzles with and without chevrons[J].Applied Acoustics,2009,70(2):321-332.

[12] 許常悅,吳丹,孫建紅.可壓縮湍流圓孔射流的大渦模擬[J].南京航空航天大學(xué)學(xué)報,2010,42(5):583-587.

[13] WAN Z H,ZHOU L,YANG H H,etal.Large eddy simulation of flow development and noise generation of free and swirling jets[J].Physics of Fluids,2013,25(12):1-27.

[14] 李棟.可壓縮湍流時空關(guān)聯(lián)與混合層噪聲的數(shù)值研究[D].北京:中國科學(xué)院大學(xué),2013.

[15] 張強(qiáng).氣動聲學(xué)基礎(chǔ)[M].北京:國防工業(yè)出版社,2012.

[16] 張兆順,崔桂香,許春曉.湍流理論與模擬[M].北京:清華大學(xué)出版社,2005.

Large Eddy Simulation of Subsonic Round Jet

LIN Jian, LAI Huan-xin

(Key Laboratory of Pressurized Systems and Safety,Ministry of Education,East China University of Science and Technology,Shanghai 200237,China)

A compressible round jet is investigated by large eddy simulation.The jet Mach number is 0.75 and the Reynolds number is 50 000.Numerical results are validated by available experimental data.Two major stages for the shear layer instability development are observed and analyzed.After this,two-point space-time correlations of axial velocity fluctuationu′,radial velocity fluctuationv′ and pressure fluctuationp′ in the shear layer are investigated.Their distribution and evolution in time and space are presented and discussed.The migration characteristics of vortices are analyzed.The study provides a base for further investigation and controlling of jet noise.

subsonic; jet flow; large eddy simulation; shear layer; correlation analysis

1006-3080(2016)05-0715-07

10.14135/j.cnki.1006-3080.2016.05.020

2016-02-23

國家自然科學(xué)基金(51576067)

林 健 (1990-),男,江蘇江陰人,碩士生,主要從事流體機(jī)械氣動聲學(xué)的研究。E-mail:lj307585298@163.com

賴煥新,E-mail:hlai@ecust.edu.cn

V211

A