亞音速軸流風扇靜轉子相互作用線譜噪聲預報*

武星宇，魏應三，靳栓寶，王 東，祝 昊，胡鵬飛，孫方旭

(海軍工程大學 艦船綜合電力技術國防科技重點試驗室， 湖北 武漢 430033)

軸流風扇被廣泛應用于民用和軍用航空領域，例如渦扇發動機，其低噪聲設計一直以來都受到國內外許多學者的關注和研究[1]。從頻譜上看，風扇氣動噪聲包含線譜噪聲和寬帶噪聲。其中，寬帶噪聲源于風扇靜轉子葉片與隨機湍流脈動的相互作用，線譜噪聲主要源于周期性轉子尾流與下游靜子的相互作用[2-3]。對于大涵道比航空發動機，線譜噪聲為其主要噪聲源，因此深入研究轉子尾流與下游靜子相互作用產生的線譜噪聲對抑制風扇噪聲具有重要意義。

隨著計算機技術的快速發展，數值分析法開始廣泛用于軸流機械非定常氣動特性和噪聲機理的研究，Verdon等[4]通過數值法求解三維線性歐拉方程，得到軸流渦輪機械靜子葉柵葉頻和倍葉頻下葉片表面非定常脈動壓力的分布。Elhadidi等[5]采用線性的歐拉方程和片條理論，計算分析了葉片彎掠設計對風扇線譜噪聲的影響。為考慮流體黏性和靜子對轉子尾流的影響，Rumsey等[6]通過數值法求解非定常雷諾平均N-S方程，計算分析了風扇靜轉子相互作用噪聲。Grace等[7]通過三維非定常雷諾平均N-S方程模擬計算靜子表面的非定常壓力，并據此預測葉柵出口線譜輻射聲功率。通過數值法可較準確地預報軸流風扇的氣動噪聲[8-9]，但數值法建模時間和計算時間較長，在風扇初期選型與參數設計時不適宜采用數值法。相對于數值計算方法，解析法在計算時間上具有很大的優勢，適用于風扇初期的低噪聲設計。

通過升力面理論計算流場與風扇葉片相互作用產生的表面壓力脈動，結合聲類比方程可進一步計算得到風扇的輻射噪聲。張偉光等[10]利用三維升力面理論計算轉子尾流與靜子葉片相互作用產生的線譜噪聲，并分析葉片彎掠設計對風扇線譜噪聲的影響，但是升力面理論不能考慮風扇靜子的安裝角[11]。為此，Carazo等[12]采用數值法模擬轉子尾流，結合Amiet[13]葉片非定常響應理論模型，計算得到對轉螺旋槳線譜聲壓值，計算結果與試驗結果相近，但是該模型未考慮葉片之間的相互作用。Posson等[14]在Glegg[15]葉柵響應函數的基礎上，推導得到葉片表面壓力脈動分布解析表達式。在此基礎上，de Laborderie等[16-17]推導出風扇線譜噪聲計算公式，結合試驗得到的轉子葉片尾流參數，得到的風扇輻射聲功率級計算值與試驗值在一倍葉頻處相近。此外，Rozenberg[18]、Sanjosé[19]、Magne[20]等也采用解析的葉柵非定常響應模型計算風扇線譜噪聲，但上述研究需要先計算葉片的非定常響應，然后通過聲類比方程計算風扇輻射噪聲，且需要數值法或試驗法提供的轉子尾流場參數。Philbrick等[21]試驗測量了葉片尾流參數，擬合曲線得到葉片經驗尾流模型，并通過該模型預報風扇線譜噪聲，發現采用試驗尾流數據和采用經驗尾流模型得到的風扇線譜噪聲結果相近。

本文將針對風扇靜轉子葉柵相互作用線譜輻射噪聲進行研究。在文獻[15]中的簡諧湍流波作用下葉柵散射場速度勢解析表達式的基礎上，考慮風扇轉子尾流的作用，通過傅里葉分解將轉子尾流分解為若干簡諧湍流波的疊加[21]，從而推導得到轉子尾流作用下靜子葉柵散射場的速度勢計算表達式，通過聲強與速度勢以及輻射聲功率與聲強的關系進一步推導得到了靜子葉柵線譜噪聲計算公式，并通過與NASA風扇試驗模型和后掠、傾斜靜子升力面計算模型對比驗證了該公式的適用范圍。并在此基礎上分析靜轉子間距、靜子葉片參數等對風扇線譜噪聲的影響。

1 葉柵線譜噪聲公式推導

1.1 轉子葉柵尾流模型

忽略葉片厚度，在半徑r處截取風扇靜轉子葉片并展開，可得到如圖1所示的靜轉子葉柵。圖1中：Ω為轉子的轉速，Sr為轉子葉柵柵距，αr為轉子葉柵安裝角，L為轉子尾流半尾跡寬度，DS為轉子后緣流向距離，(x,y)為固定在轉子葉柵上的坐標系；αs為靜子葉柵的安裝角，(xs,ys)為固定在靜子葉柵上的坐標系。

圖1 靜轉子葉柵示意Fig.1 Schematic diagram of stator and rotor cascades

大量試驗表明，單個葉片尾流可用高斯分布描述，在圖1中的o-xy坐標系下其表達式[21-23]為

(1)

式中，w0(x,y)為葉片尾流場速度梯度分布，wr為轉子中心尾跡區軸向速度虧損。

對于轉子葉柵，轉子尾流場為單個葉片尾流場的疊加，即

(2)

式中，wR(x,y)為轉子尾流場速度梯度分布，表示的是轉子尾流場中某位置的流速與平均流速的差值，Z為整數。

以T=Srcosαr為周期將式(2)傅里葉級數展開，得

(3)

式中，xr·nr=ycosαr+xsinαr，F為xr·nr的函數。

(4)

(5)

由于式(4)中的半尾跡寬度L與DS有關，為此F可改寫為F′(DS)。 結合式(3)得到

Philbrick等[21]給出轉子中心尾跡區軸向速度虧損wr和半尾跡寬L經驗計算公式。

(6)

(7)

其中：W0為轉子尾流平均速度；cr為轉子葉片弦長；Cd為轉子葉片阻尼系數，根據經驗[21]有

(8)

式中，κ為葉片稠度，α1為進口相對氣流角，α2為出口相對氣流角，λ為壓力損失系數。

考慮轉子與靜子坐標變換，在y方向上有

ys=y+Ωrt

(9)

得到周期性轉子尾流垂直于靜子葉柵的速度分量為

(10)

式中，B為轉子葉片數。

1.2 葉柵線譜輻射聲功率

Glegg[15]針對圖2所示葉柵模型提出葉柵響應函數理論，該理論通過Wiener-Holf方法求解葉柵散射場速度勢函數積分表達式，最后得到簡諧湍流脈動作用下葉柵散射場的速度勢函數解析表達式，推導得到的表達式可由式(11)表示。

圖2 平面葉柵結構示意Fig.2 Schematic diagram of plane cascade structure

ikzzd-iωt]

(11)

通過片條理論，在半徑r處截取環形風扇中Δr段，假設在分段內的流體入流參數一致[16]，考慮葉柵入流為轉子尾流，葉柵的速度勢函數解析表達式可改寫為

i(σ-2πq)yd/h+ikzzd-imBΩt]

(12)

聲強與速度勢的關系[15]為

(13)

式中，“*”表示共軛，I為聲強。

代入速度勢函數，可得到

(14)

輻射聲功率Π與聲強的關系[15]為

(15)

對于圖1所示的葉柵結構，其積分曲面定義在曲線x-yd/h=const和r

(16)

式中，Re表示取實部，這是由于輻射聲功率為實數時才能輻射出來。 通過式(16)可得到靜子第m階輻射聲功率計算表達式，即

(17)

相鄰葉片間湍流相位差角σ=kxdd+kydh，其中kxd為Glegg坐標系下的軸向波數，kyd為側向波數。不考慮葉片后掠和傾斜，波數在Glegg坐標系od-xdyd與坐標系os-xsys之間的轉換關系[16]為

(18)

1.3 葉片后掠和傾斜的影響

圖3為后掠和傾斜葉片示意圖。圖3(a)中φ為后掠角，可表示為沿葉片徑向的函數φ(z)。圖3(b)中ψ為后掠角，可表示為沿葉片徑向的函數ψ(z)。此時固定在靜子葉片前緣的坐標系os-xsyszs與Glegg坐標系od-xdydzd之間的轉換關系如式(19)所示。

(a) 后掠(a) Swept (b) 傾斜(b) Lean圖3 后掠與傾斜葉片示意Fig.3 Schematic diagram of swept and lean blade

(xd,yd,zd)T=Q×(xs,ys,zs)T

(19)

式中，Q為轉換矩陣，其定義[16]為

(20)

通過Q可得到坐標系od-xdydzd下的葉柵入流參數。

(21)

其中，[U0,Uyd,Uzd]T為od-xdydzd坐標下的入流速度，[Uxs,Uys,Uzs]T為os-xsyszs坐標下的入流速度，kzd為od-xdydzd坐標下z向入流波數，kzs為os-xsyszs坐標下z向入流波數。靜子葉柵來流波數K=[kxs,kys,kzs]T與流速U=[Uxs,Uys,Uzs]T滿足如式(22)所示關系[16]。

K·U=mBΩ

(22)

考慮葉片在zs方向上波數的變化，在坐標系os-xsyszs下的波數為

(23)

其中，p為整數，Rt為葉梢半徑，Rh為輪轂半徑。通過式(21)可得到Glegg定義坐標系下的葉柵入流參數。此外，由于傾斜和后掠角的影響，靜子在Glegg坐標系下的弦長c、垂向間距h和軸向間距d也發生了變化。

(24)

其中，cs為os-xsyszs坐標系下的弦長，Ss為靜子葉柵柵距。

2 結果和討論

2.1 不考慮葉片后掠與傾斜

NASA[22]在半消聲室中測量了不同轉速和靜子葉片數下風扇的線譜噪聲，其試驗模型如圖4所示，主要參數如表1所示。試驗測量了轉子轉速分別為1 700 r/min，1 750 r/min，1 800 r/min，1 850 r/min和1 887 r/min時的風扇線譜噪聲，在轉速為1 800 r/min時，轉子葉片葉梢馬赫數為0.34，軸向馬赫數為0.15。

圖4 NASA試驗風扇結構[22]Fig. 4 NASA′s fan test structure[22]

表1 NASA風扇試驗參數

Sutliff等[24]給出了轉子流場參數試驗測量結果，通過轉子入流流速，圖5給出了轉速為1 850 r/min時的轉子尾流場馬赫數梯度經驗公式計算云圖。通過經驗公式計算得到，在葉片中部轉子尾流中心線速度虧損約為18 m/s，Sutliff試驗測量的結果約為20 m/s，經驗結果與試驗結果相近。

圖5 周期性轉子經驗尾流公式計算云圖Fig.5 Cloud diagram of periodic rotor wake by empirical

假設導管壁面為硬壁面，即不考慮導管的吸聲作用，通過式(17)即可得到風扇靜子的線譜輻射聲功率。圖6、圖7分別給出13葉靜子和14葉靜子的線譜輻射聲功率級(參考聲壓級為10-12W)試驗值與計算值對比曲線圖，計算值包括本文公式計算值與Laborderie[16]模型計算值，式(18)中的m=1時可得到葉頻處的輻射聲功率，m=2時可得2倍葉頻處的輻射聲功率。圖6(a)與圖6(b)分別為葉頻處風扇上游和下游線譜輻射聲功率級，圖6(c)與圖6(d)分別為風扇上游和下游2倍葉頻處線譜輻射聲功率級。

(a) 風扇上游輻射聲功率級(葉頻)(a) Sound power level of the fan′s upstream in BPF

通過圖6可得，當靜子葉片數為13、轉子轉速為1 700～1 887 r/min范圍時，在葉頻處本文公式計算得到的風扇線譜輻射聲功率級與試驗值相差2～4 dB，特別是在葉頻下游工況，本文公式計算結果與試驗值相近。

通過圖7對比14葉靜子的線譜輻射聲功率級可得，在葉頻下游本文公式計算得到的輻射聲功率級與試驗值相差3 dB以內；但是在2倍葉頻處，計算值與試驗值相差較大，接近10 dB?？傮w而言，在上述計算工況，相比于Laborderie模型，本文公式在葉頻下游處的預報結果與試驗值更加接近，但是在2倍的葉頻處，本文的計算結果偏大，這是由于使用經驗模型預報轉子尾流參數時時未考慮轉子尾流沿靜子弦長的變化，導致傅里葉分解后2倍葉頻處的速度梯度偏大；而在葉頻處，本身速度梯度基數較大，弦長對轉子尾流速度梯度的影響較小。

(a) 風扇上游輻射聲功率級(葉頻)(a) Sound power level of the fan′s upstream in BPF

2.2 考慮葉片后掠

Schulten[25]采用升力面法分析了后掠葉片對風扇線譜噪聲的影響，模型參數如表2所示。

表2 Schulten后掠靜子風扇模型主要參數

表2中，歸一化葉頻為轉子葉梢馬赫數與轉子葉片數的乘積。此外，轉子阻尼系數為0.01，靜轉子間距在葉中位置保持2.6倍的靜子弦長。在采用經驗公式計算轉子尾流參數時，考慮后掠角導致的靜轉子間距的變化。圖8給出靜子的線譜輻射聲功率級本文公式計算值與Schulten升力面理論計算值對比曲線圖，由于風扇葉頻噪聲被截止，因此圖中計算頻率為2倍的葉片通過頻率。

由圖8可得，在當后掠角小于8°時，本文公式計算得到的輻射聲功率與Schulten升力面計算模型計算得到的輻射聲功率相近，但是當后掠角大于8°時，本文公式計算得到的輻射聲功率快速下降，與Schulten計算模型得到的結果相差較大。這是由于本文公式基于平面葉柵推導而來，在大后掠角下，軸向波數變小，葉柵部分模態的輻射噪聲被“截止”，但對于環形葉柵，在平面葉柵中被“截止”的輻射模態不一定會在環形葉柵中被“截止”，這導致本文計算得到的結果與Schulten計算得到的結果相差較大。

(a) 靜子上游輻射聲功率(a) Sound power level of the stator′s upstream

2.3 考慮葉片傾斜

張偉光等[10]采用升力面法分析了傾斜角對風扇線譜噪聲的影響，計算參數如表3所示，其中弦長與靜轉子間距參數通過導管直徑做歸一化處理，導管直徑為0.4 m。

表3 Schulten傾斜靜子風扇模型主要參數

圖9給出靜子線譜輻射聲功率級本文公式計算值與文獻[10]升力面理論計算值對比曲線圖，計算頻率為1倍的葉片通過頻率，其中縱坐標為計算結果與上游0°傾斜角的輻射聲級計算結果相對值。

(a) 靜子上游輻射聲功率(a) Sound power level of the stator′s upstream

由圖9可得，當傾斜角在0°～10°范圍內變化時，本文公式計算值與升力面計算值相差5 dB以內，且兩者發展趨勢一致。對于上述模型，葉片的傾斜設計反而使靜子的線譜噪聲變大。

3 參數影響分析

在2.1節NASA 13葉靜子風扇試驗模型的基礎上，圖10給出了靜轉子間距、靜子弦長、靜子傾斜角和靜子葉片數等參數變化對靜子輻射噪聲的影響。

(a) 靜轉子間距的影響(a) Influence of the stator-rotor space

考慮到下游葉頻噪聲為主要噪聲貢獻量，圖10計算的均是下游葉頻噪聲，圖10(a)為靜轉子間距變化對靜子線譜噪聲的影響曲線圖，由圖可得，靜轉子間距對靜子線譜噪聲影響較大，隨著靜轉子間距增大，風扇線譜噪聲降低，但同時降低的幅值越來越小。

圖10(b)為弦長變化對靜子線譜噪聲的影響曲線圖，由圖可得，弦長變化對NASA試驗模型的線譜噪聲影響較小。

圖10(c)為傾斜角變化對靜子線譜噪聲的影響曲線圖，由圖可得，當傾斜角小于4°時，風扇線譜噪聲變化較小，當傾斜角大于4°時，傾斜角越大，風扇線譜噪聲越小。

圖10(d)為葉片數變化對靜子線譜噪聲的影響圖，由圖可得，15葉靜子的線譜噪聲最小，其次是17靜子，20葉靜子的線譜噪聲最大。此外，在圖10(d)中16葉靜子由于與轉子葉片數相等，導致計算中出現奇點；當葉片數大于21時，葉頻噪聲被“截止”。

4 結論

本文在葉柵響應函數的基礎上推導了風扇靜轉子相互作用線譜輻射聲功率預報公式，通過轉子經驗尾流模型為該公式提供入流參數，并考慮靜子葉片后掠和傾斜設計對線譜噪聲的影響。通過與NASA試驗模型和升力面計算模型對比，驗證了本文公式適用性和適用范圍，得到如下結論：

1)本文公式能夠較準確地預報風扇線譜噪聲，但是本文公式不適用于后掠葉片線譜噪聲的計算，這是因為平面葉柵中被“截止”的聲輻射模態在環形葉柵中不一定會被“截止”。

2)葉片的傾斜設計并非都能抑制線譜噪聲，對于部分模型葉片的傾斜設計反而會增大線譜噪聲。

3)針對NASA風扇試驗模型，增大靜轉子間距和傾斜設計能夠抑制線譜噪聲，但是靜子葉片弦長對線譜噪聲影響很??；此外，葉片數對風扇的線譜噪聲影響較大，特別是可通過增加葉片數的方式來“截止”靜子葉頻噪聲。因此，合理的靜轉子葉片數匹配設計對風扇噪聲被動抑制具有重要意義。

綜上所述，本文解析公式可預報風扇靜轉子相互作用線譜噪聲，相對于常用的數值法，在計算時間上具有一定的優勢，適用于風扇葉片的初期選型和低噪聲設計。與升力面理論相比，本文公式直接通過速度勢函數計算得到風扇線譜噪聲，無須先求解葉片表面的非定常壓力脈動，其物理意義更加清晰，計算過程更加簡單。此外，升力面法無法考慮安裝角的影響。本文公式則可考慮安裝角和傾斜角對風扇線譜噪聲的影響，但是也存在一定的缺陷，對于后掠角過大的模型其預報結果不是很準確，還需要進一步優化本文公式。

致謝

清華大學黃振衛博士和海軍潛艇學院段嘉希博士在數據處理提供了幫助和指導，謹致謝意！