基于URANS/TPP 模型的風扇轉靜干涉單音噪聲預測

蔡堉楠，紀良，李旦望，夏燁，倪臻

中國航發商用航空發動機有限責任公司 先進技術研究部，上海 200241

航空發動機噪聲是飛機最主要的噪聲源，對其進行降噪設計是“安靜型”飛機研制的關鍵。對于當代大涵道比構型發動機，風扇噪聲占據主導［1］，其轉子葉片的非定常周期性尾跡與出口導流葉片（Outlet Guide Vane，OGV）干涉會產生強烈的單音噪聲。對于亞聲速工況，風扇轉靜干涉噪聲通常決定了噪聲窄帶頻譜中聲壓級的最高值，即使是在超聲速情況下，轉靜干涉噪聲也依然是風扇離散噪聲的重要組成部分［2］。隨著對噪聲抑制機理認知的深入和設計方法的進步，通過精細化葉型設計實現聲源降噪成為進一步提升風扇低噪聲性能的有效途徑［1，3-4］，代表性的手段包括風扇轉子低噪聲設計［5-7］、彎掠靜子設計［8-12］、仿生葉片設計［13-15］等。為了實現對設計的支撐，要求預測模型能夠反映詳細設計參數對聲源強度的影響。

風扇噪聲的產生伴隨著流場和聲場的耦合效應，給其聲源的計算帶來了極大挑戰。近幾十年來，該領域學者采用了計算流體力學（Computational Fluid Dynamics，CFD）方法和計算氣動聲學（Computational Aeroacoustics，CAA）方法對三維流場/聲場進行真實建模［16-19］。由于由轉靜干涉產生的離散單音噪聲擾動量級很高，而且只在特定的頻率上出現，因此為采用非定常雷諾平均數值模擬（Unsteady Reynolds Average Numerical Simulation，URANS）模型的CFD 方法求解這類問題提供了可能。相比于直接模擬或者大渦模擬方法，URANS 模型計算量大幅降低，同時該模型能夠刻畫葉型設計參數對噪聲的影響［20］。因此，URANS 模型成為目前工程中葉輪機單音噪聲源預測的有效手段［21-24］。另外，風扇流動具有時-空周期性特征，非線性諧波（Non-Linear Harmonic，NLH）方法利用截斷的傅里葉級數對流場進行逼近，可通過單通道計算實現非定常流場預測，具有比URANS 模型更高的計算效率［25-26］。但是當考慮周向非均布靜子構型、進氣攻角效應或者安裝效應等存在非周期情況時，需要開展對所有風扇葉片數值建模的全環URANS 計算。

由于聲波具有小尺度特征，其對數值方法的離散誤差敏感，因此如何在計算量可接受的前提下保證聲場預測精度，是進行風扇噪聲預測的關鍵［27］。同時，風扇噪聲在管道內聲場以模態波形式傳播。為了獲取聲源特性，需要從聲源的脈動量中提取聲模態信息［28-30］。

本文主要工作是針對大涵道比風扇構型，在對數值方法進行誤差分析的基礎上，基于URANS 模型進行風扇三維非定常流場模擬，結合三平面壓力（Triple Plane Pressure）模態匹配方法（簡稱TPP方法）從流場中提取聲源模態信息，并采用部件級聲學試驗數據驗證了預測方法的可靠性。此外，通過本文數值模型，對比驗證了彎掠靜子葉型設計的降噪效果，分析了轉靜干涉噪聲抑制基本規律。

1 預測模型

本文的風扇單音噪聲預測模型以Navier-Stokes 方程和管道聲模態傳播方程為基礎，結合URANS 模型和聲源提取方法，實現了轉靜干涉單音噪聲的預測。

1.1 基于URANS 模型的非定常流場計算

為了獲得聲源特性，需要通過CFD 方法計算獲得風扇流場中的脈動壓力場。描述非定常黏性流動的控制方程為可壓縮的Navier-Stokes 方程組。在URANS 模型中，將Navier-Stokes 方程組進行時間積分平均，得到一組以時均物理量和脈動量乘積的時均值作為未知量的非封閉方程。本文采用k-ε模型，基于渦黏性假設，建立雷諾應力項與湍動能和湍流耗散率等物理量的關系，使方程組封閉。基于數值求解方法，能夠獲得流場壓力脈動信息，從而為聲源提取提供輸入。

1.2 TPP 方法基本方程

對于聲源信息的提取，Wilson［31］最早提出波分離（Wavesplitting）方法，通過在一個交界面上對壓力和速度進行模態分解，實現對傳播波和反射波的分離。基于該思路，Ovenden 和Rienstra［28］提出TPP 方法，實現了聲模態的有效分解，同時可分辨傳播波和反射波。近些年來，相關學者對TPP 方法做了進一步擴展，Wohlbrandt 等［30］提出了XTPP 模型，能夠分離對流壓力擾動對聲源信號的影響；Guérin［32］提出改進TPP 模型，進一步考慮了旋流和徑向流動的影響。

本文 采用Ovenden 和Rienstra 的TPP 方法［28］，將發動機的進氣道和外涵道簡化成圓形或者環形管道，基于波動方程，某個頻率的聲波在管道內的傳播可采用以下形式：

式中：z、θ、r分別表示軸向、周向和徑向坐標；m和n分別為周向和徑向模態階數；為聲模態幅值；為軸向波數；Ψ(κmnr)為徑向分布函數；f為聲波頻率；c為聲速；κmn為管道特征值；A和B為常系數；Jm和Ym分別為m階的第1 類和第2 類貝塞爾函數。其中，上標“+”表示順流傳播，例如風扇外涵道出口的后傳噪聲；上標“-”表示逆流傳播，例如風扇進氣道的前傳噪聲。

常系數A和B的計算公式為

式中：上標“′”表示求導；rD表示管道半徑。

TPP 方法通過設置3 個監測面，能夠在非定常流場中提取周向、徑向模態階數及對應的幅值。對于某一周向模態m，式（1）可寫為

式中：Pn為第n個聲模態幅值；Ψn表示徑向分布函數；Φn為軸向分布函數；i=1，2，3 表示TPP監測面序號。通過在式（6）等號左右兩邊分別乘以第j個模態的徑向分布函數的共軛，并對徑向取積分可得

式（7）可改寫成以下形式：

定義代價函數：

其中：

為了使代價函數值最小，須滿足以下條件：

式中：下標l表示目標聲模態的序號。由于Pl為復數形式，因此有

根據式（12）和式（13）可得

根據式（14），得到最終形式：

式中：χln表示系數矩陣。式（16）共有N++N-個未知數。如圖1 所示，P1、P2和P3分別為3 個監測平面，其中上標“+”和“-”分別表示沿軸向和逆軸向傳播，且n≠0。此時，未知數個數與方程個數相同，通過求解式（16），可得到各階模態對應的幅值，從而完成聲模態提取。

圖1 TPP 方法示意圖Fig.1 Sketch of TPP method

對于（m，n）階模態對應的聲功率，可由式（17）得到。

式中：ρ0為流體密度；c0為流體聲速；R和r分別為環形管道的外徑和內徑；k0和kmn分別為總波數和軸向波數。

2 聲場計算誤差分析

在管道中傳播的聲場相對于流場擾動值是小量，對數值離散產生的耗散誤差極為敏感。而離散誤差與空間和時間步長相關。過大的網格尺寸會導致計算失真，而網格過密會使計算量太大。因此，在保證計算可靠前提下盡量降低網格數量是本文模型進行工程級噪聲預測的基礎。

采用模態波在環形管道中傳播的算例，討論網格尺寸與時間步長對聲場誤差的影響。根據模態波傳播特征，定義單位波長點數NPPW：

式中：Nz和Nm分別為單個波長軸向和周向投影距離所包含的網格點數。采用單頻（3，1）模態聲波入射，考慮Ma分別為0.5、0.3、0.1 和-0.3 時的工況，選取不同NPPW計算傳播單位波長后的聲壓級誤差ΔSPL，結果如圖2 所示。可以看到，隨著網格的加密，聲壓級誤差逐漸減小，其與單位波長點數的擬合公式為

圖2 網格尺寸與誤差擬合曲線Fig.2 Fit line of mesh size and error

下面討論時間步長的影響。定義無量綱單位周期時間步長數NT為關注單音聲波的周期T與CFD 物理時間步長dt的比值，即

數值計算的耗散誤差與NT存在相關性。圖3為NT與聲壓級誤差ΔSPL 的關系。可以看到，對于不同的單位周期時間步長數，在NT大于60 時，單位周期時間步長數對聲壓級誤差的影響很小，并且繼續增加單位周期時間步長數對精度提升有限，綜合考慮計算經濟性和準確性，在以下的計算中均采用NT=60。由于風扇單音噪聲通常包含多個頻率單音，越高頻的成分要求的物理時間步長越小。因此在實際計算中，式（22）中T取關注的最高頻率單音周期，以保證各頻率單音分量的準確計算。

圖3 時間步長與誤差擬合曲線Fig.3 Fit line of time step and error

通過以上分析，驗證了URANS 模型具備聲場計算可行性，同時提出空間和時間步長選取方法，以指導風扇噪聲計算。

3 風扇轉靜干涉單音噪聲預測驗證

采用國內設計的某型大涵道比風扇縮尺試驗件構型，風扇的氣動和聲學試驗在德國 Wildau 的AneCom 航空試驗（AneCom Aero Test，ACAT）噪聲試驗臺（如圖4 所示）上完成。

圖4 AneCom 航空試驗噪聲試驗臺Fig.4 AneCom AeroTest noise test facility

試驗件的風扇轉子葉片數為18，OGV 數為55。在進氣道和外涵出口的管道外壁面分別布置1 圈100 個麥克風，采用壁面齊平安裝方式并且周向非均布，用以進行聲源周向模態分解。同時，在外涵出口布置旋轉桶麥克風，旋轉桶中周向等距設置有3 列麥克風，每列60 個且軸向等距分布，用以進行聲源周向和徑向模態分解。通過試驗測得給定工況下的麥克風聲壓脈動信號和風扇軸的轉速脈沖信號，結合轉速脈沖信號對聲壓脈動信號開展重采樣，并按葉片通過頻率（Blade Passing Frequency，BPF）的整周期倍截斷，以準確獲取聲壓測點在BPF 位置的聲壓幅值和相位特性。在給定的頻率下，管道內各階模態在測點位置的聲壓特性是1 組線性無關基向量，基向量的線性組合即為測得的聲壓幅相特性。因此，利用最小二乘法求解測得的聲壓在給定聲模態構成的線性空間下的坐標，即可得到各階聲模態在給定頻率下的幅值，從而可獲取聲模態，這為噪聲數值評估方法提供了試驗數據對比驗證。

風扇試驗件直徑約為0.45 m。管道聲模態結果與風扇/增壓級的轉子、靜子數相關，若改變葉片數進行約化計算，預測的模態階數將與實際情況不符，因此本文采用全環計算方案。風扇轉靜干涉噪聲主要由風扇轉子與外涵出口導葉干涉產生，且為了盡量減少網格點數，對于內涵僅考慮第1 級靜子導葉S0，風扇轉子葉片數為18，OGV 數為55，S0 葉片數為67。風扇進口設置總溫總壓，外涵和內涵出口設置為靜壓出口。采用Numeca AUTOGRID 完成風扇部件網格劃分，流場計算及后處理采用ANSYS CFX 完成。

計算網格設置如圖5 所示，在風扇轉子前進口和OGV 后出口通過拉伸網格的方法設置耗散段，以實現無反射邊界。對于計算網格，根據式（21）可知，保證預測的最小目標波長的NPPW應大于25，同時每個徑向波長內網格點數應大于15，最小目標聲波周期的時間步長數NT為60。在計算網格前端和后端分別設置TPP 監測面，為前傳和后傳聲源提取提供輸入。采用URANS方法進行非定常流場計算，湍流模擬選用k-ε模型，在風扇轉子和外涵靜子表面y+接近于30，總網格量約為5×107，其中y+為無量綱數，表示網格節點到壁面的距離與流動特征尺度之比。

圖5 計算網格設置Fig.5 Computing mesh settings

計算了適航規定的邊線（OP1，103.7%）、飛越（OP2，93.9%）和進場（OP3，60.0%）轉速工況。首先進行定常計算，將其結果作為非定常模擬的初場以保證計算穩定性。非定常計算總計時間步數約為4 000 步，在每個時間步內進行10 次內循環。采用了448 核CPU 并行運算，總計算時間為341 h。

圖6 為103.7%轉速下非定常計算迭代過程某一軸向位置不同徑向高度處脈動壓力隨時間變化曲線。在迭代大于2 800 步時，脈動壓力趨近于穩定周期解，表明非定常流在計算過程中已趨近于收斂。

圖6 非定常壓力隨時間變化曲線Fig.6 Time history of fluctuating pressure

觀察外涵出口處流量和壓比值并與試驗數據對比。如圖7 所示，3 個適航工況流量誤差分別為0.58%、-0.21% 和0.41%，壓比誤差為2.77%、0.12% 和-0.10%。通過試驗對比分析，表明流場計算的正確性。

圖7 流量和壓比預測值與試驗數據對比Fig.7 Comparison of prediction and experiment results of mass flow and pressure ratio

圖8 為80%葉高處轉子和靜子間的靜熵云圖，圖8 中采用了相同的云圖標尺。可以觀察到轉子的尾跡進入外涵后與下游靜子產生相互作用，它是非定常壓力擾動的重要來源之一。轉子尾跡被靜子前緣分割，并與葉片表面邊界層相互作用，在葉片表面產生高熵增區域，造成邊界層的損失增加。在靜子下游，轉子尾跡與外涵通道中的主流繼續發生摻混，這種影響一直發展到外涵出口。

圖8 不同轉速工況瞬時靜熵云圖Fig.8 Instantaneous contours of static entropy of different operation conditions

觀察風扇靜子后方外涵道內BPF 的瞬時脈動壓力實部。從圖9 中可以看到，由于轉靜干涉產生的脈動壓力在周向和徑向上呈明顯周期分布，符合管道聲模態波特征。需要指出的是，對于OP1 和OP2，其1 階和2 階BPF 均導通，但對于低轉速工況OP3，1BPF 截止，因此觀察了OP3的2BPF 和3BPF 下的脈動壓力。

圖9 OGV 出口脈動壓力實部云圖Fig.9 Contours of real part of fluctuating pressure at OGV outlet

下面對邊線（OP1，103.7%）、飛越（OP2，93.9%）和進場（OP3，60.0%）3 個工況的后傳噪聲和進場（OP3，60.0%）工況前傳噪聲進行聲模態提取。在對監測面時域壓力脈動數據進行頻域轉換后采用TPP 方法進行聲模態分解和聲功率計算。最終，采用聲模態階數和幅值信息的形式，給出風扇轉靜干涉單音聲源。

根據Tyler 和Sofrin［33］的轉靜干涉噪聲理論，聲波周向模態階數m符合以下特征：

式中：j為諧波階數；k為整數；F為轉子葉片數；S為OGV 數或者S0 葉片數，分別代表轉子與外涵OGV 或者第1 級內涵靜子導葉的干涉噪聲模態。

由于URANS/TPP 模型的限制，本文僅對BPF 下的主要占優模態進行討論。圖10 和圖11為不同轉速下風扇后傳聲模態與AneCom 試驗數據的對比，干涉噪聲周向模態階數和誤差如表1 所示。為了保證與試驗數據一致，OP1 和OP2 工況聲模態幅值采用聲功率表示，OP3 采用聲壓級表示。可以觀察到，由轉靜干涉產生的主要導通單音模態被準確捕捉到，聲模態階數試驗結果一致。其中，在OP1、OP2 和OP3 后傳工況，由轉子與OGV 干涉產生的模態占主導；在OP3前傳工況3BPF 下，除了轉子/OGV 干涉產生的m=-19 模態，預測模型還捕捉了由轉子/S0 干涉產生的m=-13 模態。計算出各工況下幅值絕對值平均誤差為4.7 dB，驗證了本文方法對占優單音噪聲預測的有效性。

表1 各工況占優聲模態階數和預測誤差Table 1 Order and prediction error of dominant acoustic modal under different operation conditions

圖10 不同工況下后傳聲模態與試驗數據對比Fig.10 Comparison of prediction data with experiment data under different operation conditions for rear acoustic mode

圖11 OP3 前傳聲模態與試驗數據對比Fig.11 Comparison of prediction data with experiment data at OP3 for front acoustic mode

4 低噪聲彎掠靜子設計的數值分析

風扇靜子OGV 彎掠造型是抑制轉靜干涉噪聲的有效手段。彎掠葉片主要增加了上洗速度沿徑向的相位變化，從而使不同徑向位置的非定常載荷相互抵消以降低聲源強度。圖12 為彎掠OGV 示意圖，從葉根到葉尖，當葉片前緣隨徑向位置增大而向下游變化時，葉片后掠，掠角α取為正；反之，葉片前掠，掠角為負。當葉片隨徑向位置增大而沿周向相位角發生變化時，稱之為彎葉片，彎角β與風扇轉子旋轉方向相同時為正。

圖12 彎掠OGV 示意圖Fig.12 Sketch of leaned and swept OGV

基于三維升力面理論對風扇彎掠靜子對噪聲的影響進行快速解析計算［11］，結果如圖13 所示。考慮結構強度的限制，選取側彎角5°和后掠角23°作為設計方案，進行三維數值分析。需要指出的是基準構型靜子存在15°的后掠角，但無側彎角。

圖13 OGV 彎掠角對聲功率影響分析Fig.13 Impact of leaned and swept angle of OGV on acoustical power

采用與第3 節中一樣的URANS/TPP 耦合數值計算方法，計算了103.7%轉速工況的風扇轉靜干涉噪聲。首先觀察彎掠設計對氣動性能的影響。如表2 所示，采用彎掠設計后，外涵流量和壓比變化分別為0.09%和0.26%，未對氣動性能造成不良影響。

表2 彎掠靜子設計對氣動性能影響Table 2 Impact of leaned and swept design on aerodynamic performance

風扇干涉噪聲的強度取決于葉片對上洗速度的響應而產生的表面非定常載荷分布。圖14 為靜子葉片表面非定常壓力幅值分布，可以觀察到1BPF 幅值大于2BPF。同時，吸力面幅值大于壓力面。吸力面幅值熱區集中在靜子前緣和葉身靠近前緣處，壓力面幅值主要集中在靜子前緣。

圖14 彎掠靜子表面脈動壓力幅值分布Fig.14 Fluctuating pressure amplitude contours on surfaces of leaned and swept OGV

提取彎掠設計和基準構型的靜子葉片吸力面非定常力實部，從圖15 可以觀察到，非定常力在葉片的徑向和軸向呈周期分布，在不同單音頻率下分布呈現不同形式。根據Envia 和Nallasamy［8］的理論分析，聲場強度由葉片表面壓力差決定，可以通過對非定常力在葉片表面積分獲得。徑向相位變化越大，對葉片不同區域積分時的相位相消越明顯，從而可降低管道聲場強度。通過數值結果分析可以觀察到，采用彎掠設計后非定常力在徑向上的周期性相比于基準構型更加明顯。例如在2BPF 下，在葉片中下段前緣，彎掠構型的力分布周期數顯著增大。

圖15 基準構型和彎掠構型靜子表面脈動壓力實部分布Fig.15 Real part of fluctuating pressure contours on surfaces of leaned/swept OGV and baseline OGV

為了驗證降噪效果，對靜子后傳噪聲進行TPP 聲源提取。圖16 給出聲模態分布，可以觀察到彎掠設計對主要單音模態階數未產生影響，但是降低了聲模態幅值。對1BPF 下m=18 以及2BPF 下m=-19 和m=36 主要模態有明顯的抑制效果。通過功率計算，根據表3 預測結果，本文數值方法和三維升力面模型存在一定差異，主要原因可能在于建模方法上的區別，但是在1BPF 和2BPF 下，不同方法的結果均表明聲功率幅值得到抑制。上述結果表明，通過三維數值分析，實現了對低噪聲彎掠靜子設計的降噪基本規律的分析和降噪效果驗證。

表3 降噪量預測結果Table 3 Prediction results of noise reduction

圖16 基準構型和彎掠構型靜子后傳聲模態Fig.16 Rear acoustic mode of baseline OGV and leaned/swept OGV

5 結論

本文基于URANS/TPP 模型進行了風扇轉靜干涉單音噪聲預測研究，主要結論如下：

1）URANS 模型能夠對三維管道模態聲波傳播進行計算，通過合理地選擇網格尺寸和時間步長，可以實現對型號中實際風扇構型的可靠預測。

2）基于URANS 模型對風扇非定常流場進行了數值計算，結合TPP 模型可得到轉靜干涉噪聲模態信息。該耦合模型能夠反映流場和噪聲的細節變化，通過部件級試驗數據對比，驗證了該方法的可靠性。

3）通過本文數值模型計算，反映了彎掠靜子設計對葉片表面非定常力分布的影響。通過增加載荷力的分布周期，達到相位相消的效果，從而實現轉靜干涉聲源的抑制，在不影響氣動性能的前提下實現風扇的降噪設計。

4）需要指出的是，聲源預測結果中出現了本底噪聲模態干擾，這可能是在聲源提取面處由旋流、徑向流動和對流壓力波動產生的，此外聲源提取處的非聲壓擾動也可能會產生額外誤差。因此本文僅對主占優模態進行討論，采用XTPP 等改進聲源提取方法［30，32］能夠提升預測結果精確度。