基于Mohring聲類比發動機噴流噪聲數值計算

閆國華 ，馮叔陽 ，劉 勇 ，李建福

（中國民航大學航空工程學院1，基礎實驗中心2：天津 300300；3.空軍裝備部駐沈陽地區第二軍事代表室，沈陽 110043）

0 引言

商業航空發動機涵道比不斷增大，使得在飛機起飛過程中主要噪聲源之一的噴流噪聲已降低到與風扇噪聲相當的量級[1]。然而，國際民航組織（International Civil Aviation Organization，ICAO）采用的第5 階段噪聲標準[2]對航空發動機的降噪設計以及飛機的噪聲適航審定提出了更加嚴格的要求。為了適應適航更加嚴厲的航空噪聲管制以及完善擁有自主噪聲預測評估能力的民用航空器噪聲審定技術研究體系，噴流噪聲數值計算研究依舊是重點。

在航空器噪聲適航審定時，為了能夠快速地評估噪聲級，常采用基于半經驗公式的噴流噪聲數值計算方法，如：SAE 算法[3]、Pao 算法[4]和Stone 算法[5]等，然而半經驗公式因其固有的局限性，無法精確地解析噪聲頻譜結構、分析輻射特性，且有嚴格的使用范圍。伴隨著計算機技術的發展，航空發動機氣動噪聲的數值計算已經成為計算噴流噪聲的主流方法。使用高精度流體計算的結果，運用各種聲類比方法進行氣動噪聲計算的混合法成為目前最受歡迎的氣動噪聲數值計算方法。作為氣動聲學建立的標志，Lighthill[6-7]采用聲類比方程求解了自由空間下的湍流噪聲；Curl[8]考慮了靜止固體邊界對噪聲的影響，將其等效為在整個固體邊界上分布偶極子源；Williams 等[9]采用廣義函數法擴展到運動固體邊界對聲音的影響，得到著名的福克斯·威廉姆斯-霍金斯方程（Ffowcs Williams-Hawkings equation，FW-H）聲類比方程，成為如今氣動聲學遠場聲學計算最常用的聲類比方法；Powell[10]、Howe 等[11]通過探究渦流與聲源的關系建立和完善了Powell渦聲公式；Mohring[12]對非均質流的氣動噪聲進行研究，提出了適合高馬赫數噪聲的Mohring 聲類比方程。相比于Lighthill 聲類比和FWH 聲類比只考慮聲音在靜止介質中的傳播，Mohring聲類比考慮了聲音在具有流動非均勻介質中的傳播，其應用領域比Lighthill 聲類比更廣泛，更適合高亞聲速氣動噪聲的數值計算；金海波等[13]應用該類比方法對風機噪聲進行數值分析，并驗證了其在旋轉機械噪聲計算中的準確性。

本文采用大渦模擬（Large Eddy Simulation，LES）的方法對大涵道比發動機噴管的簡化縮尺模型進行高亞聲速噴流數值分析，耦合Morhing 聲類比方法提取聲源，對其進行聲傳播計算，得到近場和遠場的聲學求解并進行分析。

1 基本理論

1.1 Mohring聲類比

Mohring 聲類比方程由N-S 方程推導出來，考慮到流動和聲場的耦合以及方程的非線性給求解帶來困難，將聲學現象假設為線性問題，把聲場分為近場聲源區域和遠場聲傳播區域，使流動與聲場相分離。并在Lighthill 方程的基礎上考慮流動特性，在聲學建模中加入了流動效應。忽略粘性張量以及熱傳導的影響，從而得到Mohring方程[12]

式中：Na為與有限元計算中網格相關的量。

式（5）描述了由方程右端的聲源分布產生的聲傳播控制方程，其左邊為聲傳播算子，右邊為聲源項，由面聲源SΓ和體聲源SΩ組成

式中：?為流場中包含參量ρ和v的矢量。

1.2 離散傅里葉變換

由于聲傳播的計算是以頻域形式進行的，而聲源是從時域的流場計算結果中得到的。因此，需要通過離散傅里葉變換處理時域信號的流場結果[14]

式中：N為非定常流場采樣的時間步數；Qn為頻域信號；qk為時域信號。

1.3 大渦模擬（LES）

在瞬態噴流流場的計算中，采用介于直接求解N-S方程和雷諾時均模擬之間大渦模擬（LES）來提取噴流噪聲聲源，該方法的精確性在高亞聲速湍流噪聲的聲源求解[15]和渦聲方程聲源的求解[16]中都得到驗證。LES首先將N-S 方程進行濾波處理，其中對平均流動影響較大的大尺度量渦流通過N-S 方程直接求解。濾波后的控制方程為

2 數值模擬方法的算例驗證

噴流噪聲數值模擬方法的驗證研究采用NASA的SMC002 號4 鋸齒小型金屬噴管的冷噴流試驗[17]完成。SMC002 噴管模型尺寸如圖1 所示。SMC002 噴管為單流通通道的收縮噴管，進、出口直徑分別為65.85、53.6 mm，出口安裝4 個Chevron 型鋸齒，齒長32 mm，齒彎角為5°。

圖1 SMC002噴管模型尺寸

進口總壓為178200 Pa，總溫為288.15 K；出口總壓為97700 Pa，總溫為280.2 K；采用LES 計算流場，保證出口馬赫數可達0.9，該試驗為亞聲速無限空間內淹沒射流。

LES 計算的SMC002 噴管中心線速度分布如圖2所示。其中，縱坐標為噴流速度u與噴流最大速度Uj的比值，橫坐標為噴口出口中心向方向距離X與噴口出口直徑Dj的比值。與NASA試驗SMC002號噴管結果進行對比可知，LES 方法的計算結果在靠近噴口出口的位置與試驗值符合得更好，即在核心區與試驗結果吻合更好。

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

按照NASA 試驗要求，以噴口出口中心為原點，2.54 m 為半徑，在180°范圍內每隔10°建立19 個遠場聲壓接受點，采用Mohring 聲類比提取體聲源并計算聲傳播，同時也采用不考慮非均勻流場的Lighthill 聲類比提取體聲源并計算聲傳播。依次從60°～150°，每隔10°取10個主要場點的數值計算值與文獻[17]中的試驗結果進行對比。使用試驗中不同極坐標角度下的場點總聲壓級（overall sound pressure level，OASPL）進行驗證，數值計算值與試驗值的對比和計算誤差見表1。

表1 數值計算值與試驗值的對比和計算誤差

從表中可見，Morhing 聲類比的數值計算值大多大于試驗值，最大計算誤差為1.97%，在航空器噪聲適航噪聲評估的工程計算中，這一計算精度可以接受。而Lighthill 聲類比數值計算值大多小于試驗值，但是最大誤差達到4.6%。相比之下，Morhing 聲類比計算誤差更小，更接近試驗值。這是由于Morhing 聲類比相比于Lighthill 聲類比加入了非均勻背景流，考慮了流場的速度剪切伴隨的聲傳播折射效應，因此在高馬赫數的流動計算中更接近真實值。

3 噴流流場數值計算

3.1 噴管幾何模型

研究對象為某型大涵道比發動機噴管的簡化縮尺模型，如圖3 所示。將不影響噴流流場的發動機部件進行簡化：去掉進氣道的唇口部分；將表面渦流發生器等部件略去，使表面光滑；去掉外涵道出口支板。以風扇環形出口部分為外涵道進口，進口高度為11 mm，出口高度為6 mm；以渦輪環形出口為內涵道進口，進口高度為4.6 mm，出口高度為3.5 mm。大涵道比發動機噴管尺寸如圖4所示。

圖3 大涵道比發動機噴管的簡化縮尺模型

圖4 大涵道比發動機噴管尺寸

3.2 計算域及網格

為保證噴流流場的充分發展，確定外流域為長600 mm、直徑為140 mm 的圓柱，如圖5 所示。采用結構網格對計算域進行六面體網格劃分，噴流流通區域進行O型拓撲。為了滿足大渦模擬計算要求，設置壁面附面層第1 層網格距離y+=1，使得第1 個網格位于湍流粘性底層區域。將整個噴流核心區域進行網格局部加密，x、y、z方向的網格尺寸為0.2 mm，使其能夠精確地捕捉小尺度渦流。在核心區以外的區域，徑向上以第1層厚度為0.3 mm，增長率為1.15，增大到10 mm為止布置網格。軸向上均勻布置間距為10 mm網格。最終，經網格無關性驗證之后，確定總節點388萬，總網格數為389萬。網格細節分布如圖6所示。

圖5 噴管外流域

圖6 噴管周圍網格

3.3 邊界條件及求解設置

首先使用RNGk-e湍流模型進行時均穩態噴流流場計算，以此穩態流場作為LES 瞬態求解的物理初場，加快噴流流場的收斂。為達到起飛時刻的發動機噴流狀態，確定邊界條件為：（1）外涵入口來流質量流量為0.16023 kg/s，來流溫度347 K；（2）內涵入口來流質量流量為0.03077 kg/s，來流溫度為978 K；（3）遠場邊界為壓力遠場，溫度為300 K，壓力為103125 Pa。

3.4 流場計算結果

噴流速度時均分布如圖7 所示。整個噴流流場可大致劃分為3 個區域：混合區、過渡區和充分發展區。混合區內包含了勢流核心和混合層，勢流核心以內涵噴流為主要部分，其噴流最大速度可達532 m/s，馬赫數為0.9。距離尾錐末端4 mm 的徑向特征截面時均速度分布如圖8 所示。從圖中可見，混合層包含3 部分：外涵環形勢流與周圍靜止空氣形成的外涵混合層；內涵環形勢流與外涵環形勢流形成的內涵混合層；內涵環形勢流圍繞的圓錐形的中心混合層。

圖7 噴流速度時均分布

圖8 距離尾錐末端4 mm的徑向特征截面時均速度分布

中心線速度曲線如圖9 所示。從圖中可見，從尾錐處開始，受粘性的影響，中心混合層內的速度在內涵噴流速度的影響下逐漸增大，混合層面積逐漸減小，在270 mm 處，中心混合層結束，環形的內涵勢力核心融合為一股勢流，速度突然增大并向后發展到過渡區開始衰減，整個混合區長300 mm。過渡區域內即為時均速度衰減區，這一區域內噴流速度迅速衰減到100 m/s 以下，過渡區域長260 mm。之后即為自維持流動的充分發展區，流體保持低速流動。

圖9 發動機噴管中心線速度分布

湍流動能時均分布如圖10 所示。從圖中可見，內外涵的勢流核心區湍流動能幾乎沒有，可視為層流區域。噴流強湍流區域出現在內涵核心層內，伴隨內涵勢流核心向后延伸，在內涵勢流核心區域末端開始集中，在過渡層內集中在呈羽毛狀上下對稱分布的2個區域，主要向后段充分發展區內衰減。

圖10 湍流動能時均分布

瞬態噴管出口渦量分布如圖11 所示。從圖中可見，在外涵勢流核心與周圍靜止空氣的剪切層內，在強迫擾動的作用下，在外涵出口處即發生了渦卷起，從而形成大尺度渦環結構。開始時渦環流動保持層流狀態，渦環在向下對流過程中，剪切層內的周向失穩導致渦環外部形成渦辮結構，擾動的進一步發展使渦環大概在內涵出口處扭曲破碎，渦環結構消失，充分發展為外涵混合層內湍流。隨后大尺度渦破碎為小尺度渦，與內涵剪切層內的小尺度渦在過渡層內融合。

圖11 瞬態噴管出口渦量分布

由于內外涵之間剪切層內速度梯度沒有外涵與空氣之間剪切層內的速度梯度大，其渦卷起形成的大尺度渦環結構一直對流到內涵環形勢流開始匯聚之處，擾動發展明顯增大，渦環外部失穩，小尺度渦辮結構開始形成并增長，加速了渦環破碎為大尺度渦結構。而大尺度渦結構的進一步發展使得上下剪切層內的渦在勢流核心末端融合，渦環結構消失，融合后的大尺度渦結束了穩定的勢力核心的發展。此后，大尺度渦結構螺旋扭曲發展，快速破碎為小尺度渦，在過渡區內形成充分發展的湍流，并向下游對流形成湍流尾跡，噴流則完全融入周圍環境。小尺度渦結構將其從大尺度渦環中獲得機械能在粘性的作用下不斷地耗散，到充分發展區，渦量幾乎消失。

聲功率采用從Lighthill 方程推導而來的Proundman 方程[20]計算，將單位體積內各項同性的湍流產生的聲能表達為

式中：α為常數；u2為均方根速度波動；c2為聲速。

噴流噪聲聲功率級如圖12 所示。從圖中可見，聲功率級最大為163 dB。強聲源區域為圍繞內涵勢流核心環形區域，主要集中在核心勢流末端部分，可達到150 dB 以上，與破碎渦量形成的強湍流區域相當吻合。在整個噴流混合區以及過渡區域內，除層流部分的內外涵勢流核心，聲功率級在渦量分布的地方達到了130 dB 以上。發動機噴流的聲功率整體分布在1 個邊界和中心線夾角11°的錐形區域內，向后方衰減。在靠近噴口出口的邊界上聲功率衰減比較迅速，但伴隨流動，越靠近下游，聲功率的衰減級逐漸緩慢。

圖12 噴流噪聲聲功率級

4 噴流聲場數值計算

4.1 聲學數值計算流程

采用計算流體力學（Computational Fluid Dynamics，CFD）計算得到產生聲源的流場脈動量，以此為基礎根據計算航空聲學（Computational Aeronautical Acoustics，CAA）進行聲傳播計算，完成聲場求解。基于瞬態CFD 和CAA 聯合求解氣動噪聲的方法兼顧了計算資源以及求解精度，成為成熟且廣泛使用的計算方法。混合數值算法求解噴流噪聲的基本流程如圖13所示。

圖13 混合數值算法求解噴流噪聲的基本流程

在采用LES 方法進行瞬態噴流流場求解的基礎上，利用瞬態流場密度、速度矢量的時域結果進行噴流噪聲體聲源的計算，將所得的時域形式的體聲源積分插值到聲學網格上，并對其進行離散性傅里葉變換轉換為頻域形式的體聲源。同時也將穩態流場結果線性插值，以求得Mohring 體聲源。使用Mohring 聲類比理論所得的勢流中的聲傳播算子進行噴流噪聲聲傳播的求解。

4.2 聲學網格及計算域

聲學近場采用有限元計算噴流噪聲聲源以及近場聲傳播，遠場場點聲壓級采用無限元插值求解，插值階次為5，采用Prandtl-Glauert 變換處理后的Helmholtz聲傳播算子進行聲傳播計算[14]，從而精確求解整個噴流噪聲聲場。

在有限元計算區域，以流體計算域作為聲源區；并在此基礎上在x、y、z3 個方向上向外增加10 層網格區域，以此區域為聲傳播區域，將聲傳播區域邊界面全部定義為無限元邊界，噪聲沒有能量損耗以及反射，全透射傳播出去，聲學計算域如圖14 所示。在聲學遠場區域部分，以噴流為軸向X 軸，徑向為YZ 方向。在XY 平面，以內涵噴口出口中心為圓點，1 m 為半徑，在180°范圍內每隔10°建立19 個遠場聲壓接收點，如圖15所示。

圖14 聲學計算域

圖15 遠場場點

4.3 聲場計算結果

選取4 個頻率下的噴流聲場分布，如圖16 所示。聲源區域主要分布在噴流剪切層內以及在噴流勢流核心末端的湍流混合區域內，其中強聲源主要由在勢流核心下游的強湍流脈動形成，最高可達188 dB。隨著頻率的提高，噴流產生的聲輻射有更加均勻廣泛的波輻射陣面，而低頻下的聲傳播則明顯向下游集中。根據Tam 等[20-21]關于亞聲速噴流噪聲的產生機理的研究，亞聲速噴流噪聲主要歸結于湍流的摻混，不同尺寸的渦結構形成了不同頻率的噴流噪聲，由小尺度渦主導形成的高頻噪聲方向性不明顯，因此傳播范圍廣泛；而由大尺度渦主導形成的低頻噪聲則有著相對明顯的方向性，且主導方向為下游方向。

圖16 噴流聲場分布

500 Hz 下噴流流向截面的聲場如圖17 所示。從圖中可見，4 個截面到距離尾椎末端的距離分布為0.1、0.2、0.3、0.4 m。在此方向上，噴流噪聲以噴流的中心軸為圓心呈圓形向外輻射。強聲源區域與渦量的分布呈現很大相似性。在噴管尾椎末端即0.1 m處，聲源主要分布在剛破碎的外涵剪切層內以及內涵渦環結構內。伴隨外涵剪切層大尺度渦對流耗散以及內涵剪切層內渦環破碎；到0.2 m 處聲源分布在中心圓形區域內，并向中心部分集中；到了勢流核心末端0.3 m 處，完全發展為大尺度渦，湍流充分發展，聲源區域均勻分布在圓形區域內；再往下游到0.4 m處，大尺度渦破碎為小尺度渦耗散，聲源區域也明顯衰弱。

圖17 500 Hz下噴流流向截面的聲場

取場點中60°、90°、120°、150°的頻譜特性，如圖18 所示。從圖中可見，發動機噴流噪聲主要峰值集中在低頻范圍內，越靠近噴流流動方向，聲壓級越高。這是由于越靠近下游越能感知到大尺度渦，導致下游方向低頻噪聲明顯，隨著頻率升高，各方向角的聲壓級都在降低，在1000～2500 Hz，從125 dB快速降低到105 dB，之后衰減速度變緩，到100 dB左右趨于穩定。

圖18 噴管各方向角的頻譜特性

遠場場點噪聲總聲壓級（Overall Sound Pressure Level，OASPL）指向性如圖19 所示。從圖中可見，噴流噪聲整體上是沿噴流方向輻射，主要指向噴流流動方向，越靠近噴流方向，聲壓級越大，由于噴流背景流的影響，使其聲傳播呈現這一種狹長單一方向的輻射特性。從發動機所在位置處開始，在發動機壁面上方40°～60°內，由于發動機結構對聲音的影響，使其聲音不太容易輻射出去，聲壓級最小，此范圍內聲壓級均在157.2 dB 左右。在80°～100°內外涵噴流速度還未完全衰減，而在90°方向的內涵噴口又剛開始形成勢流核心，使得此區域內時均速度大，同時外涵剪切先在此區域內形成充分發展的湍流區，使得聲壓級有一個小突升，在100°到達157.5 dB。之后，越靠近噴流方向，湍流發展越充分，流場速度大，聲壓級也越大，在正對噴流處聲壓可達166 dB。同時噪聲的前傳也比較明顯，從40°開始，越靠近上游聲壓級越大，發動機正前端的聲壓級可達161 dB。

圖19 遠場場點噪聲總聲壓級指向性

各場點的總聲功率級曲線如圖20 所示。從圖中可見，總聲功率級的分布與總聲壓級的分布基本一致，聲能量的分布形成了噪聲指向性。在發動機所在位置處，40°～60°內，由于發動機結構的影響，聲功率較小；80°～100°內的場點相較于下游場點離噴流區域較遠，但此處噴流流場時均速度較大，同時從流場的聲功率級云圖中可見，在圍繞外涵環形勢流形成的外涵強湍流區域內形成較大的聲功率，因此使得此處的聲傳播有個小突升，在100°時達到153 dB。此后隨著流場速度的減小和內涵噴流的形成的強聲功率區域還未形成，聲功率級有所降低。但之后場點越靠近下游，越能感知到更快的流動速度帶來的內涵強湍流區域更高的聲能量，聲功率級也在增大，在正對噴流處達到最大157 dB。

圖20 各場點的總聲功率級曲線

5 結論

（1）對于高亞聲速噴流，在強迫擾動下，外涵剪切層內的大尺度渦環結構向下對流過程中，在保持短暫的層流狀態后即破碎為大尺度渦結構，并繼續向下游破碎為小尺度渦耗散；內涵剪切層內渦環結構流動距離相對較長，在擾動作用明顯增大后，破碎為大尺度渦，在勢流核心末端大尺度渦在噴流中心線融合，結束勢流核心，隨后完全破碎為小尺度渦，湍流充分發展。

（2）噴流噪聲的強聲源區域與渦量的分布有很大相似性，集中在由破碎渦形成的強湍流區域內。由大尺度渦主導的低頻噪聲的傳播比較集中，而小尺度渦主導的高頻噪聲傳播更廣泛。噴流噪聲聲壓級峰值主要分布在低頻噪聲范圍內，越靠近噴流下游，低頻下的聲壓級越高。

（3）該數值方法與NASA 試驗結果對比，最大計算誤差為1.97%，結果可靠。相較于Lighthill 聲類比，由于其考慮了聲音在流動介質中的傳播，對于高馬赫數的流動噪聲而言，Morhing 聲類比的計算結果更接近試驗值，能得到更精確的預測結果。同時相比于適航上常用的半經驗公式，能夠得到近場噪聲聲輻射結果，并且能夠更好地耦合流場，可以為發動機噪聲結構設計以及航空器噪聲適航審定提供參考，也可以用于飛機發動機噴流噪聲預測。