基于混合數(shù)值算法的鋸齒形噴管氣動噪聲仿真

閆國華，馮叔陽

(中國民航大學(xué)航空工程學(xué)院，天津 300300)

1 引言

隨著大涵道比渦扇發(fā)動機(jī)的廣泛應(yīng)用，作為現(xiàn)代渦扇發(fā)動機(jī)主要噪聲源之一的噴流噪聲已經(jīng)降低到與風(fēng)扇噪聲相當(dāng)?shù)牧考塠1]。但隨著國際民航組織(International Civil Aviation Organization，ICAO) 開始采用更加嚴(yán)厲的第五階段噪聲標(biāo)準(zhǔn)[2]，我國國產(chǎn)大飛機(jī)也將面臨更嚴(yán)峻的適航審定。因此在建立自己的擁有噪聲預(yù)測評估能力的民用航空器噪聲審定技術(shù)研究體系方面，在整機(jī)噪聲中依舊占主導(dǎo)地位的噴流寬頻噪聲依舊是研究重點，

對噴流噪聲的產(chǎn)生機(jī)理、噪聲預(yù)測以及降噪手段等方面的探索有很多的理論和實驗研究。TAM、VISWANATHAN等人[3][4]揭示了亞音速噴流噪聲的產(chǎn)生主要是由于湍流的摻混導(dǎo)致的。Kawai S[5]等分析了亞音速到超音速下自由剪切層的不穩(wěn)定性對噴管整體基礎(chǔ)流場的影響。Samimy[6]等人對Tab齒形的降噪原理進(jìn)行分析，揭示了鋸齒產(chǎn)生的對旋渦對噴管流場剪切層內(nèi)的摻混的增加作用。NASA開展多項實驗[7][9]，在消聲室對一系列不同鋸齒結(jié)構(gòu)小型噴嘴(Small Metal Chevron Nozzles)在不同測試條件下進(jìn)行噪聲測量，這一系列試驗既取得了噴管流場的時均值和湍流脈動量的信息，也得到了噴管遠(yuǎn)場檢測點的噪聲總聲壓級和頻譜信息，其實驗數(shù)據(jù)常被用來檢驗噴流噪聲預(yù)測方法的可靠性和精度。

通常使用兩種計算方法來計算預(yù)測噴流噪聲，一種基于半經(jīng)驗公式的預(yù)測模型，如SAE算法[10]、Pao算法[11]和Stone算法[12]等。二是基于聲類比的數(shù)值計算方法，如基于FW-H方程的Ligthill聲類比方法[13]以及考慮流體可壓縮性的Mohring聲類比[14]。盡管半經(jīng)驗?zāi)Ｐ驮诳焖僭u估噪聲級方面具有先天的優(yōu)勢，但是其無法精確的解析噪聲頻譜結(jié)構(gòu)，不能對噪聲的產(chǎn)生機(jī)理進(jìn)行分析，且有嚴(yán)格的使用范圍。如今，伴隨著計算機(jī)技術(shù)的發(fā)展，數(shù)值計算已經(jīng)漸漸成為了計算噴流噪聲的主流方法。

采用計算流體力學(xué)(Computational Fluid Dynamics，CFD) 和計算航空聲學(xué)( Computational Aeronautical Acoustics，CAA) 混合方法對噴管進(jìn)行噪聲數(shù)值計算[15]。首先對噴管的噴流流場進(jìn)行瞬態(tài)數(shù)值計算，在瞬態(tài)流場的計算結(jié)果的基礎(chǔ)上提取噴流噪聲體聲源，對其進(jìn)行聲傳播計算，得到近場和遠(yuǎn)場的聲學(xué)求解，并采用 NASA 的實驗數(shù)據(jù)[7]對流場以及噪聲計算結(jié)果進(jìn)行對比驗證。

2 噴管幾何模型

研究對象是NASA實驗用的SMC002號4鋸齒小型金屬噴管[7]，根據(jù) NASA 提供的噴管數(shù)據(jù)完成建模，其三維模型如圖1所示。噴管為單流通通道的收縮噴管，尺寸如圖2所示，噴管進(jìn)口直徑65.85mm，出口直徑53.6mm，出口安裝4個Chevron型鋸齒，齒長32mm，齒彎角5°。

圖1 噴管模型

圖2 噴管尺寸

3 噴流流場數(shù)值計算

3.1 流場計算方法

精確的流場解析是聲學(xué)仿真的基礎(chǔ)，為了保證流場計算量適中以及流場計算的精度要求。瞬態(tài)噴流流場的計算選擇介于直接求解N-S方程和雷諾時均模擬之間大渦模擬(LES)，來提取噴流噪聲聲源。LES首先將N-S 方程進(jìn)行濾波處理，其中對平均流動影響較大的大尺度量渦流通過N-S方程直接求解，而主要起耗散作用的小尺度渦流因其各項同性的特點，則通過亞格子模型求解，從而保證流場解析的精度。

(1)

其中由濾波操作而得到的大尺度量定義為

(2)

式中積分遍布整個流體域，G(r，x)是空間濾波函數(shù)，噴流流場的小尺度運(yùn)動的尺寸和結(jié)構(gòu)由其決定。

將過濾函數(shù)作用于N-S方程，得到紊流的控制方程

(3)

(4)

3.2 流體計算域及網(wǎng)格

為保證噴流流場充分發(fā)展，以半徑為500 mm、長2070 mm的圓柱體作為噴流計算域，將SMC002噴管包裹。采用結(jié)構(gòu)網(wǎng)格對計算域進(jìn)行六面體網(wǎng)格劃分，噴流流通區(qū)域進(jìn)行O型拓?fù)洹榱藵M足大渦模擬計算要求，設(shè)置壁面附面層第一層網(wǎng)格距離y+值為1，使得第一個網(wǎng)格位于湍流粘性底層區(qū)域。為了更好捕捉噴管出口剪切層內(nèi)的渦流，在噴管出口區(qū)域進(jìn)行局部加密，生成初始網(wǎng)格，總節(jié)點數(shù)184萬，總網(wǎng)格數(shù)190萬。計算域網(wǎng)格分布細(xì)節(jié)如圖3所示。

圖3 噴嘴周圍的網(wǎng)格

為了確定噴流流場瞬態(tài)計算中網(wǎng)格數(shù)量與計算結(jié)果之間的無關(guān)性，再準(zhǔn)備一種90萬數(shù)量網(wǎng)格和150萬數(shù)量網(wǎng)格，分別對三種不同密度網(wǎng)格進(jìn)行試算，監(jiān)測出口中心線馬赫數(shù)并進(jìn)行比較，結(jié)果如圖4所示。稀疏網(wǎng)格計算結(jié)果與初始網(wǎng)格和加密網(wǎng)格結(jié)果有微小差別，而初始網(wǎng)格與加密網(wǎng)格計算結(jié)果幾乎相同。因此，初始網(wǎng)格可以滿足LES計算需求。

圖4 網(wǎng)格無關(guān)性驗證

3.3 邊界條件及求解設(shè)置

為了加快噴流流場的收斂，首先使用RNGk-e模型進(jìn)行穩(wěn)態(tài)噴流流場計算，將此穩(wěn)態(tài)流場作為LES瞬態(tài)求解的物理初場。選取NASA實驗中的冷噴流實驗來確定邊界條件[7]，具體邊界條件如下：

1) 進(jìn)口邊界條件：噴嘴入口采用壓力入口，進(jìn)口總壓178200 Pa，總溫288.15 K；

2) 出口邊界條件：噴嘴出口采用壓力出口，總壓97700 Pa，總溫280.2 K；

3) 遠(yuǎn)場邊界條件：遠(yuǎn)場采用壓力遠(yuǎn)場，總溫總壓與出口邊界條件相同。

3.4 流場計算結(jié)果

圖5為SMC002噴流速度時均分布云圖，其中噴口流速達(dá)到297.773 m/s，接近 NASA[7]實驗要求的出口流速0.9馬赫的工況。由此可見，該實驗為亞音速無限空間內(nèi)淹沒射流。整個噴流流場可大致劃分為三段，混合區(qū)包括了保持與噴管出口速度相當(dāng)?shù)膭萘骱诵囊约盎旌蠈樱渲袆萘骱诵膮^(qū)長度0.25 m；過渡區(qū)域內(nèi)即為時均速度衰減區(qū)，這一區(qū)域內(nèi)噴流速度急速衰減到100 m/s左右，區(qū)域長度與混合區(qū)大致相同；之后即為流體自維持流動的充分發(fā)展區(qū)，流體保持低速流動。

圖5 SMC002噴流流場結(jié)構(gòu)

圖6為SMC002噴流湍流動能的時均分布云圖，可以明顯看出勢力核心區(qū)域湍流動能極小，可視為層流流動區(qū)域；而整個混合區(qū)域內(nèi)圍繞在勢流核心周圍的混合層為強(qiáng)湍流區(qū)域，中間區(qū)域湍流動能最大，并一直向后延神，并超過勢流核心區(qū)延伸到過渡區(qū)內(nèi)，強(qiáng)湍流區(qū)域圍繞勢流核心上下對稱，沿軸線方向面積增大，在過渡區(qū)內(nèi)開始衰減，在最后的充分發(fā)展區(qū)維持低水平狀態(tài)。

圖6 SMC002噴流湍流動能結(jié)構(gòu)

圖7是噴口出口截面時均速度云圖，噴流高速區(qū)域在兩鋸齒之間向外部低速區(qū)域擴(kuò)散，相比于鋸齒頂尖由噴流最大速度到外界靜止空氣間巨大的速度梯度，高速噴流從兩鋸齒之間逐漸地向外滲出，兩鋸齒之間的流體速度梯度減小。Chevron型鋸齒有效的增大了噴流面積，從而起到了更好的摻混作用。

圖7 噴管出口速度云圖

圖7為噴管出口截面的湍流動能云圖，鋸齒頂尖湍流動能極小，而強(qiáng)湍流區(qū)域出現(xiàn)鋸齒兩邊，沿鋸齒頂尖的軸向?qū)ΨQ，并在兩個鋸齒之間開始向外部衰減擴(kuò)散，只要的湍流分布區(qū)域向兩鋸齒之間集中。

圖8 噴管出口截面湍流動能云圖

在定常計算基礎(chǔ)上進(jìn)行瞬態(tài)計算，圖9噴管出口附近渦量云圖，噴管出口之后，高速氣流與外部低速氣流形成了脫落渦，渦流主要分布在混合區(qū)的混合層內(nèi)和過渡層和充分發(fā)展層內(nèi)，及渦流集中強(qiáng)湍流區(qū)域以及湍流擴(kuò)散區(qū)域內(nèi)，勢流核心區(qū)域內(nèi)渦流分布不大，整體的渦流伴隨氣流向后高速流動；圖10為噴管出口截面渦量云圖，出口處，高速氣流與外部低速氣流之間形成了強(qiáng)烈的剪切環(huán)，渦流集中大聲在兩鋸齒之間，向四周擴(kuò)散。

圖9 噴管出口附近的渦量云圖

圖10 噴管出口截面渦量云圖

由此可見，高速射流噴入靜止空間，強(qiáng)烈沖擊靜止空氣，與其進(jìn)行了急劇的摻混，高速的噴流氣體與周圍靜止空氣形成了環(huán)形的混合層，混合層內(nèi)的渦流帶來強(qiáng)烈湍流，渦流伴隨高速流體向后流動衰減，湍流也在過渡發(fā)展區(qū)之后開始衰減；而鋸齒的加入，增大了兩鋸齒間的摻混，也使得兩鋸齒間的渦流增大，湍流增大。

圖11為瞬態(tài)0.025 LES計算的噴管出口之后中心線速度分布曲線。其中縱坐標(biāo)u/Uj為噴流速度u與噴流最大速度Uj的比值，橫坐標(biāo)X/Dj為噴口出口中心向方向距離X與噴口出口直徑Dj的比值。與NASA實驗[7]SMC002號噴管結(jié)果進(jìn)行對比，可以看出，LES方法的計算結(jié)果在靠近噴口出口的位置與實驗值符合得更好，即核心區(qū)和實驗結(jié)果吻合更好。

圖11 SMC002噴管中心線速度分布

4 聲場計算

4.1 CFD/CAA 混合求解流程

CFD計算得到產(chǎn)生聲源的流場脈動量，以此為基礎(chǔ)進(jìn)行CAA聲傳播計算，完成聲場求解。基于瞬態(tài)CFD和CAA聯(lián)合求解氣動噪聲已經(jīng)成為成熟且廣泛使用的計算方法。圖12為混合數(shù)值算法求解噴流噪聲的基本流程圖。

圖12 聲場求解流程圖

SMC002號噴管聲學(xué)求解過程如下：

1) 使用LES方法進(jìn)行噴流瞬態(tài)流場計算，獲得噴流流場密度、速度矢量；

2) 利用瞬態(tài)流場密度、速度矢量進(jìn)行時域的噴流噪聲體生源提取，在Actran中icfd模塊完成；

3) 利用離散傅里葉變換( Discrete Fourier Transform，DFT)，將時域的聲源信號變換到頻域信號；

4) 將噴流噪聲的聲源點集映射到聲學(xué)網(wǎng)格，基于Ligthill聲類比理論進(jìn)行噴流噪聲聲場求解。

4.2 Ligthill聲類比方法

考慮到流動和聲場的耦合以及方程的非線性給求解帶來困難，因此將聲學(xué)現(xiàn)象假設(shè)為線性問題，把聲場分為近場聲源區(qū)域和遠(yuǎn)場聲傳播區(qū)域，使流動與聲場相分離。由此從N-S方程得到Ligthill方程

(5)

對(5)進(jìn)行頻域轉(zhuǎn)[16]換得

(6)

式(6)中：ω為頻率；Tij為Ligthill應(yīng)力張量，可表示為：

(7)

式(7)中：v為流體流速；c0為靜止介質(zhì)的聲速；p為流體壓力；xij為空間固定坐標(biāo)。對(7)式進(jìn)行有限元變分離散化即可得到以下聲學(xué)解

(8)

4.3 離散傅里葉變換

由于聲傳播的計算是以頻域形式進(jìn)行的，因此，必須將時域的噴流體聲源進(jìn)行離散傅里葉變換，公式如下

(9)

式(9)中，N為非定常流場采樣的時間步數(shù)；Qn為頻域信號，qk為時域信號。

4.4 聲學(xué)網(wǎng)格和聲學(xué)計算域

近場聲壓采用有限元計算噴流噪聲聲源以及近場聲傳播，遠(yuǎn)場場點采用無限元插值求解，從而精確求解整個噴流噪聲聲場。

4.4.1 網(wǎng)格尺度

4.4.2 聲學(xué)計算域

在有限元計算區(qū)域，以流體計算域作為聲源區(qū)；并向外增加10層網(wǎng)格區(qū)域，以此區(qū)域為聲傳播區(qū)域，將聲傳播區(qū)域邊界面全部定義為無限元邊界，噪聲沒有能量損耗以及反射，全透射傳播出去。聲學(xué)近場有限元計算域如圖13所式。

圖13 聲學(xué)計算域

聲學(xué)遠(yuǎn)場區(qū)域部分，按照NASA實驗要求建立遠(yuǎn)場場點[7]。噴流為軸向X軸，徑向為YZ方向。在XY平面，以噴口出口中心為圓點，2.54 m為半徑，在180°范圍內(nèi)每隔10°建立19個遠(yuǎn)場聲壓接收點，如圖14所示。

圖14 遠(yuǎn)場場點示意圖

4.5 聲場計算結(jié)果

近場噪聲聲壓級如圖15所示。與流場結(jié)果對比可以看出，聲壓級比較大地方集中在噴流混合區(qū)域內(nèi)，特別是勢流核心周圍的強(qiáng)湍流區(qū)域的混合層內(nèi)，相比于核心勢流的聲壓，混合層的聲壓級要高其10dB以上。噴管出口的混合層內(nèi)，大量渦量間的相互作用產(chǎn)生各種脈動量，因此強(qiáng)烈脈動湍流造成了噴流噪聲的產(chǎn)生。整體噴流噪聲沿著噴流方向發(fā)展，聲壓集中分布區(qū)域與湍流動能較大區(qū)域比較耦合，噴流噪聲實質(zhì)上反應(yīng)了噴流湍流混合過程的渦旋變化。但是伴隨頻率升高，整體聲壓級逐漸降低，聲壓級的分布也更均勻。

圖15 近場聲壓級

噴管出口近場聲壓級如圖16所示，可以明顯看出，聲壓級較高的地方集中在兩鋸齒之間周圍，圍繞齒缺處形成環(huán)狀，并向外擴(kuò)散。在1500.94 Hz下，鋸齒之間周圍聲壓級達(dá)到最大143 dB，而齒頂尖區(qū)域只有119 dB；2501.56 Hz下，鋸齒之間周圍聲壓級達(dá)到最大144d B，而齒頂尖區(qū)域只有121 dB；3302.6 Hz下，鋸齒之間周圍聲壓級達(dá)到最大130 dB，而齒頂尖區(qū)域只有110 dB；4552.85 Hz下，鋸齒之間周圍聲壓級達(dá)到最大124 dB，而齒頂尖區(qū)域只有109 dB；由此可見兩鋸齒相比于鋸齒頂尖，鋸齒周圍聲壓級要高20 dB左右。隨著頻率升高，噴管出口整體聲壓級減小，而鋸齒與齒頂尖的聲壓差距也有所減小。對比噴管出口的流場結(jié)果以看出，強(qiáng)聲壓區(qū)域同樣于強(qiáng)湍流區(qū)域耦合；而鋸齒的安裝則大大增加了流場的摻混，使得流場的大尺度渦量破碎成小尺度渦量，從而使得齒頂尖噪聲減小。

圖17為遠(yuǎn)場場點噪聲總聲壓級(OSPL)指向性圖。噪聲總聲壓級越靠近噴流流動方向越大，正對噴口方向的噴流噪聲可達(dá)113 dB，從30°到120°范圍內(nèi)聲壓級變化平穩(wěn)，在105 dB上下浮動。

圖16 噴口出口近場聲壓級

圖17 噪聲指向性特征圖

取場點中0°、60°、90°、120°、150°為主要測量點，其頻譜圖如圖18所示。其中，正對噴流方向的0°場點各個頻率下的聲壓級都是最大。可以看出，整個噴流噪聲聲壓級較大的部分都集中在2000 Hz以下的中低頻段；2000Hz到8000 Hz段噪聲變化平穩(wěn)；而從8000Hz以上開始，聲壓級急劇衰減，高頻段噪聲因其在空氣中容易被損耗有明顯降低。

圖18 主要測量點的噪聲頻譜曲線

為驗證結(jié)果準(zhǔn)確型，取4個主要場點與文獻(xiàn)[7]的實驗結(jié)果進(jìn)行對比。使用實驗中不同極坐標(biāo)角度下的場點總聲壓級(OSPL)進(jìn)行驗證，數(shù)值計算結(jié)果與實驗值的結(jié)果比較如表1所示。

表1 數(shù)據(jù)對比與計算誤差

由表1可見，數(shù)值計算結(jié)構(gòu)大多小于實驗值，但是最大誤差不超過3.5%，且大部分結(jié)果都偏小，出現(xiàn)此誤差的原因可能是計算量的限制，使得在噴流流場的瞬態(tài)仿真計算上，網(wǎng)格不夠密，從而無法捕獲更破碎的渦，脈動渦量的計算誤差從而使得聲學(xué)仿真計算出現(xiàn)誤差。但在航空器噪聲適航噪聲評估的工程計算上，這個計算精度可以接受。

5 結(jié)論

1)文章綜合非定常噴流流場數(shù)值模擬方法、Ligthill聲類比方法、有限元結(jié)合無限元方法實現(xiàn)并完成了對單通道的Chevron型鋸齒噴管噴流噪聲聲輻射數(shù)值模擬，相比于適航上常用的半經(jīng)驗公式以及傳統(tǒng)基于FW-H方程的聲輻射數(shù)值模擬，能夠得到近場噪聲結(jié)果，以及詳細(xì)聲輻射特性，并且與流場信息得到更好的耦合，可以為發(fā)動機(jī)噪聲結(jié)構(gòu)設(shè)計提供參考。

2)噴流產(chǎn)生噪聲的重要原因是高速氣流沖擊周圍靜止空氣，形成的湍流帶來了強(qiáng)烈脈動量形成噪聲，噴流噪聲主要在噴流流場中的混合層內(nèi)產(chǎn)生，渦量大的地方湍流動能大，同時噴流噪聲大。

3)Chevron型鋸齒的安裝使得噴管的噴流出口與外界靜止空氣的接觸面積增大，摻混增大，兩鋸齒之間湍流動能相比于齒頂尖增大，脈動量也增大，在中地頻段，齒頂尖聲壓級比鋸齒之間聲壓級小20 dB左右。

4)噴流噪聲主要是低頻噪聲，并且具有明顯指向型，指向噴流流動方向，在遠(yuǎn)場噪聲計算方面，噴流噪聲數(shù)值計算結(jié)果與實驗結(jié)果相比，誤差小于3.5%，結(jié)果可靠。可用來進(jìn)行飛機(jī)發(fā)動機(jī)噴流噪聲預(yù)測，為噴管設(shè)計以及航空器噪聲適航審定提供參考。