應用于風扇后傳噪聲的三維降噪聲襯優(yōu)化設計技術

邱昇

(中國航發(fā)商用航空發(fā)動機有限責任公司， 上海 200240)

隨著航空運輸量的日益增加、大型客機功率和體積的不斷增大等，飛機噪聲污染問題越來越嚴重，給機場附近居民帶來越來越多的煩惱。人們提出嚴格適航要求來控制起飛和降落時飛機噪聲。“美國聯(lián)邦適航條例(Federal Airworthiness Regulations,FAR)第36部[1]”“歐洲聯(lián)合適航標準第36部”以及“ 國際民航公約(International Civil Aviation Organization,ICAO)附件16”等都對商用飛機做出了非常嚴格的適航噪聲規(guī)定。如今，為了追求高效率，民用渦扇發(fā)動機的涵道比越來越大，風扇噪聲變成了現(xiàn)代高涵道比發(fā)動機的主要的噪聲來源。

目前，大型民機無一例外均采用降噪聲襯來滿足噪聲控制要求。傳統(tǒng)的穿孔板-蜂窩聲襯被廣泛地應用到渦扇發(fā)動機中。但傳統(tǒng)聲襯的吸聲頻帶較低，降噪能力有限。近些年來，許多研究致力于吸聲頻譜較寬的寬頻聲襯。最近幾年，中外進行了一系列多自由度聲襯降噪設計研究[2-7]。同時，Jones等[8]在美國國家航空航天局(National Aeronautics and Space Administration， NASA)蘭利聲襯實驗室中驗證了多自由度聲襯的降噪效果。新型聲襯包括在共振腔中加薄膜的聲襯[9]、變腔深的非均勻[10]聲襯、泡沫金屬聲襯、具有多孔腔壁和柔性阻尼腔壁的混合聲襯、可調(diào)節(jié)空腔體積的變幾何自適應聲襯、偏流聲襯以及應用超材料的聲襯[11]。這些寬頻聲襯有效地擴寬了聲襯的吸聲帶寬，但其制造、維護和安裝等方面面臨新的問題。

研究者一方面追求擴寬聲襯的降噪頻帶，發(fā)展多自由度降噪聲襯，另一方面解決聲襯在進氣道和外涵道中的應用問題。進氣道聲襯主要降低向上游傳播的風扇葉片通過頻率單音噪聲，而外涵道聲襯則用于衰減燃燒室、渦輪有關的噪聲和向下游傳播的風扇后傳噪聲(也稱“外涵道噪聲”)[12]。與緩慢變化的環(huán)形進氣道相比，三維外涵道管道結(jié)構復雜，流道細長、狹窄且曲率變化大，呈非環(huán)形的C形橫截面。管內(nèi)周向存在非均勻流動、多個支桿和分墻，且聲襯在周向和徑向存在分段和不連續(xù)現(xiàn)象。這些三維復雜特征均制約了聲襯在外涵道中的應用和增加了三維外涵道聲襯優(yōu)化設計的難度[13]。

已有的航空發(fā)動機降噪工作主要是針對風扇前傳噪聲，進行降噪聲襯優(yōu)化設計工作[14-16]。研究者已進行大量風扇前傳噪聲聲傳播預測工作。與進氣道問題相比，通過發(fā)動機外涵道向后傳播的管道噪聲的預測，對于常規(guī)數(shù)值格式來說是個挑戰(zhàn)。其主要復雜性在于外涵道流體和外部流之間出現(xiàn)了一個剪切層。它會與噪聲計算相互干擾。經(jīng)常會污染噪聲解且使邊界處理變得復雜。這使得，可計算風扇前傳噪聲的方法，如頻域有限元和無限元等，無法直接應用于風扇后傳噪聲。早期，Gabard等[17]發(fā)展了半無限長發(fā)動機管道輻射噪聲的解析解，用于預測風扇后傳噪聲的遠場指向性。在數(shù)值求解方面， Bailly等[18]通過在控制方程中加入非線性項來抑制不穩(wěn)定波。Bogey等[19]和Zhang等[20]通過人為刪除一些剪切項來抑制風扇后傳噪聲的剪切層失穩(wěn)問題。然而，這些數(shù)值處理都是非物理和不合理的。Redonnet等[21]通過求解三維全Euler方程來求解真實三維發(fā)動機外涵道噪聲。商業(yè)軟件ACTRAN DGM模塊通過采用離散加遼金方法，在時域求解線性歐拉方程來獲得風扇后傳噪聲解[22]。但是，這些外涵道噪聲預測方法的顯著缺點是計算開銷太大。

綜上可知，風扇噪聲產(chǎn)生后，既可通過發(fā)動機進氣道向前傳播，也可通過外涵道向后傳播。為了降低風扇噪聲，故作為其主要組成部分之一的風扇后傳噪聲，值得研究和進行降噪設計。

民用航空發(fā)動機的真實外涵道結(jié)構復雜，存在非均勻流動、邊界反射等現(xiàn)象。基于解析預測模型的方法在計算效率和預測精度上可以滿足工程需要，但是只能開展二維聲傳播分析。基于數(shù)值方法的有限元方法和CAA(computational aeroacoustics)方法，雖然計算精度滿足工程要求，也能考慮三維聲傳播問題，但是計算效率低，采用該方法進行三維聲襯設計，其計算開銷讓設計者無法承受。風扇后傳噪聲預測問題的復雜性，進一步增加了其降噪聲襯設計的難度。

目前，針對三維外涵道管道聲傳播問題，主要的研究是外涵道聲傳播預測技術。由于其計算復雜且計算量大，故而很少有人研究三維外涵道聲襯優(yōu)化設計問題，外涵道聲襯設計問題尚未解決。首先發(fā)展精度滿足要求的高效外涵道管道聲傳播方法，然后基于此方法，發(fā)展三維外涵道聲襯設計方法。在聲傳播方法和聲襯設計方法方面有如下特點。

(1)抓住外涵道內(nèi)聲波主要沿著負方向傳播的特性，合理地將雙向聲傳播的耦合問題近似成單向聲傳播問題。建立聲源邊界條件和聲襯阻抗邊界條件，推導得到新的簡化噪聲控制方程。然后，兼顧精度和計算效率，發(fā)展高精度求解拋物線近似方程技術，即采用高階色散關系保持(dispersion relation preserving， DRP)格式進行空間離散和采用低耗散、低色散龍格庫塔方法(low dispersion Runge Kutta,LDDRK)進行時間推進計算，實現(xiàn)了高效求解風扇后傳噪聲問題。

(2)提出風扇后傳噪聲三維降噪聲襯設計方法和流程。聯(lián)合Kriging代理模型技術，構建代理模型來減少單次風扇后傳噪聲聲傳播計算時間，以進一步降低優(yōu)化計算開銷。在滿足計算精度前提下，高效率實現(xiàn)了三維外涵道聲襯設計。

1 風扇后傳噪聲預測方法

1.1 聲源建模

管道中聲場的分布可以表示成一系列周向、徑向模態(tài)的疊加，表達式為

p(x,r,θ)=∑∑[AJm(αmnr)+BYm(αmnr)]×

(1)

1.2 拋物線近似技術推導新方程

頻域?qū)α鞑▌臃匠虨?/p>

(2)

式(2)中：p為聲壓；x1表示軸向坐標；k0表示聲波波數(shù)；β是與馬赫數(shù)相關的系數(shù)；γ表示橫截面；M表示流場馬赫數(shù)。

為了求解全三維風扇后傳噪聲問題，變換坐標系，在正交曲線坐標系ξ={ξ1,ξ2,ξ3}中求解。在曲線坐標系中，拉普拉斯算子可寫為

(3)

式(3)中：ξ1、ξ2、ξ3表示曲線坐標系；h1、h2、h3分別表示笛卡爾坐標向曲線坐標的投影系數(shù)。

將這些坐標變換關系代入原控制方程[式(2)]中，可得

(4)

引入三個算子σ、R、W，即

(5)

拋物線控制方程可寫為

(6)

此方程為p的簡單二次方程，其解可寫為

(7)

(8)

對于式(8)，通過因式分解，可得

(9)

對于風扇后傳噪聲問題，聲波主要沿著負方向傳播，故可以把雙向聲傳播的耦合問題近似成單向聲傳播問題，合適的控制方程為

(10)

(11)

采用泰勒技術展開方法，有

(12)

代入風扇后傳噪聲控制方程，得到新方程為

(13)

式(13)中：an為泰勒展開的系數(shù)。

寫成顯示形式為

(14)

該方程為一階偏微分方程，可采用數(shù)值方法進行求解。

1.3 空間離散格式

本文中采用Tam等[23]發(fā)展的4階色散關系保持(DRP)格式來進行空間離散。一階導數(shù)有限差分的一般形式為

(15)

式(15)中：x為坐標；f為函數(shù)；Δx為網(wǎng)格間距；aj為格斯系數(shù)；M和N為網(wǎng)格節(jié)點數(shù)目。

通過對x取傅里葉變換，方程變?yōu)?/p>

(16)

通過比較式(16)等號兩邊，很明顯有效波數(shù)為

(17)

為了確保微分格式的色散關系在感興趣波數(shù)范圍內(nèi)與偏微分方程的色散關系相近，必須選擇aj使下述積分誤差最小，積分函數(shù)式為

(18)

對E求偏導數(shù)，給出一個線性代數(shù)方程，可求得系數(shù)aj的值。

1.4 時間推進格式

對于氣動聲學問題，僅僅考慮穩(wěn)定性問題是不夠的，還需考慮聲波的耗散和色散誤差。本文中采用Hu等[24]提出的4階4/6級低耗散、低色散龍格庫塔方法(LDDRK)解決這個問題。以下列控制方程為例，即

(19)

式(19)中：U和F分別表示速度和壓力。

第一步為

(20)

式(20)中：h(1)、h(2)、h(3)、h(4)、h(5)、h(6)分別表示每一步的推進步長；Un表示第n個時間步的物理量；Δt表示時間步長。

第二步為

(21)

1.5 聲襯邊界條件

對于含聲襯的情況，邊界條件考慮聲襯的位移，故需包括動力學條件和運動學條件。前者與聲襯處法向速度和壓力有關，后者表示流動是連續(xù)的，無流動分離。動力學條件即阻抗邊界條件可表示為

(22)

式(22)中：p為聲壓；ψ為聲襯表面位移；un為聲襯表面法向方向的聲速度；Z表示聲襯阻抗;w表示頻率；負號表示應用在外表面，正號表示應用于內(nèi)壁面。

通過減去表面位移，式(22)可合并為

(23)

動量方程為

(24)

式(24)中：ρ0是來流密度值；xn是軸向坐標。

于是，可得到聲襯阻抗邊界條件為

(25)

式(25)中:x1表示x坐標;ρ0表示密度;xn表示法向坐標。

2 風扇后傳噪聲三維降噪聲襯優(yōu)化設計

2.1 優(yōu)化代理模型建模

Kriging模型的典型表達式為

y(x)=β+Z(x)(26)

式(26)中：β是常數(shù)；Z(x)是期望為0、方差為σ2,且服從高斯分布的隨機數(shù)。Z(x)的協(xié)方差為

cov[Z(xi),Z(xj)]=σ2R(xi,xj)=

(27)

式(27)中:σ2是隨機過程的方差；R(xi,xj)是兩個樣本點xi和xj之間的高斯相關函數(shù);n代表設計參數(shù)的個數(shù)，θk>0及1

(28)

當θk給定后，σ2和β的計算公式為

(29)

超參數(shù)的最大似然估計通過數(shù)值優(yōu)化得到。當通過基礎樣本點求的超參數(shù)后，在設計空間的一個新設計點的函數(shù)值計算公式為

(30)

式(30)中：x*表示新設計點；1表示單位矢量；r為相關向量。

方差s2(x*)的估算公式為

(31)

式(31)中：r(x)=R(x,x1),R(x,x2)…,R(x,xn)是新設計點與原基礎樣本點之間的相關向量；s2為方差。

2.2 設計變量和目標函數(shù)

本文中優(yōu)化變量為聲襯阻值R和聲襯抗值X。其上下限范圍為

(32)

式(32)中：R為聲阻；X為聲抗。

優(yōu)化目標函數(shù)定為降噪聲襯的傳遞損失(transmisson loss,TL)值，即風扇后傳噪聲通過降噪聲襯后的能量損失，其計算公式為

(33)

式(33)中：Win為進口邊界上的入射聲能量;Wra為出口邊界的輻射總能量。

2.3 三維降噪聲襯設計流程

本文發(fā)展了基于拋物線近似技術和Kriging代理模型的三維聲襯優(yōu)化設計技術，以進一步降低優(yōu)化計算開銷。提出詳細設計流程如下(圖1)。

圖1 風扇后傳噪聲三維降噪聲襯設計流程Fig.1 The three-dimensional bypass liner optimization framework based on parabolic approximation model and Kriging model for aft fan noise

(1)確定聲襯阻值和抗值的上下限制，通過實驗設計方法得到初始樣本點。

(2)基于拋物線近似方法計算所有樣本點的風扇后傳噪聲聲傳播，得到傳遞損失值。

(3)基于初始樣本點的設計變量和傳遞損失值，建立響應面函數(shù)。

(4)采用多島遺傳算法，尋找聲襯阻值和抗值的優(yōu)化值，該遺傳算法是針對傳統(tǒng)遺傳算法早熟、易陷入局部最優(yōu)的缺點，在傳統(tǒng)遺傳算法基礎上增加許多“島”，產(chǎn)生更優(yōu)秀個體基因，幫助算法跳出局部最優(yōu)點，達到全局最優(yōu)[25]。本文算例中，島的數(shù)目設為30，個體數(shù)量設為50，進化代數(shù)設為50、交叉概率為0.5，變異概率為0.6，遷徙率為0.5。

(5)優(yōu)化判斷準則定為前后兩次優(yōu)化算例的傳遞損失值的差值是否小于0.001。判斷得到的優(yōu)化值是否滿足最優(yōu)判斷準則，若滿足，結(jié)束優(yōu)化計算，否則將得到的優(yōu)化值加入已有的初始樣本點中，重復步驟(3)～步驟(5)，直到滿足最優(yōu)判斷準則。

3 三維降噪聲襯設計結(jié)果

3.1 驗證算例一

首先采用NASA流管試驗數(shù)據(jù)[26]驗證本文噪聲預測方法。測試段的截面為寬50.8 mm，高50.8 mm，聲襯樣本長度為406.4 mm。其他試驗參數(shù)詳見表1。管內(nèi)噪聲聲傳播計算結(jié)果如圖2所示。從中可知，管內(nèi)敷設聲襯區(qū)域，噪聲被有效吸收。沿軸向分布的聲壓級曲線如圖3所示，編寫的程序預測結(jié)果與試驗結(jié)果吻合，驗證了本文方法精度滿足要求。

圖2 噪聲頻率2 000 Hz下的聲傳播云圖Fig.2 The sound propagation results at 2 000 Hz

圖3 沿軸向分布的聲壓級對比曲線Fig.3 Comparison of SPL results between the predictions and the experimental data

表1 NASA流管試驗參數(shù)Table 1 The parameters of NASA tube experiment

此外，常規(guī)計算聲學方法，如CAA方法，其聲壓項偏微分方程[27]為

(34)

采用64位操作系統(tǒng)，Intel Xeon 2.8 GHz處理器，16 GB內(nèi)存的小型工作站計算該三維管道聲傳播算例，至少需要10 min，而本文方法僅需不到1 min時間，在滿足計算精度前提下，節(jié)省了90%的計算時間。

3.2 驗證算例二

其次，本文中采用已有的項目試驗數(shù)據(jù)驗證本文外涵道聲傳播計算方法。外涵道聲傳播試驗在中國飛機強度研究所完成。整個試驗噴管由以下部分組成：流管收縮段、聲源部分及聲源模態(tài)發(fā)生器、出口導向葉片(outlet guide vane，OGV)支板、聲傳播段及模態(tài)測量分析裝置、聲襯段、尾椎和遠場麥克風陣列等。整個試驗裝置如圖4所示。

圖4 試驗噴管裝置及其幾何Fig.4 The test nozzle device and its geometry

流場馬赫數(shù)為0.3，噪聲頻率為3 000 Hz，管道聲模態(tài)(8,0)，本文方法通過高精度求解拋物線近似方程，獲得了外涵道聲傳播結(jié)果，如圖5中實線所示，每一小格表示20 dB。試驗測試了觀測角30°～120°內(nèi)的50倍噴管直徑處遠場聲壓，其數(shù)據(jù)如圖5含正方形的實線所示。從兩者的對比可知，本文方法計算精度高，與試驗數(shù)據(jù)結(jié)果吻合。

圖5 遠場聲壓級結(jié)果對比Fig.5 Comparison of far-field sound pressure level

3.3 典型外涵道三維聲襯設計

采用Iob等[28]中的典型航空發(fā)動機外涵道構型(圖6)，背景流動特征如表2所示。噪聲頻率為8 282 Hz，可傳播聲模態(tài)為(9,1)，基于聲源公式，計算得到其聲源分布如圖7所示。模態(tài)采用商業(yè)軟件ANSYS FLUENT計算外涵道流場，風扇處設置為壓力出口邊界條件，外部為壓力遠場邊界條件。計算得到的平均背景流場結(jié)果如圖8所示。

圖6 航空發(fā)動機外涵道示意圖Fig.6 The sketch of aeroengine bypass duct

圖7 聲源分布云圖Fig.7 The distribution of sound source

圖8 管內(nèi)馬赫數(shù)分布云圖Fig.8 The mach number distribution of bypass duct at approach condition

在三維聲襯聲阻抗優(yōu)化計算過程中，首先，建立一些聲襯聲阻抗初始樣本點，進行風扇后傳噪聲傳播計算，得到所有初始樣本點設計變量值和傳遞損失值之后，建立響應面模型，然后，采用多島遺傳算法在響應面上尋優(yōu)，多次迭代后最終得到最優(yōu)阻值和抗值。基于上述理論，發(fā)展了三維外涵道聲襯優(yōu)化設計程序，并開展了典型外涵道三維聲襯設計。

表2 背景流動條件Table 2 The mean flow properties at bypass duct

由3.1節(jié)可知，其他常規(guī)方法，如CAA方法，其聲壓項偏微分方程含9個偏導數(shù)，其三個速度項的偏微分方程中均含7個偏導數(shù)，而本文求解的控制方程僅含1項偏導數(shù)，顯著降低了計算開銷。在已開展的二維軸對稱發(fā)動機管道聲傳播算例中[29]，噪聲頻率為f=866.2 Hz，聲模態(tài)為(13,1)，網(wǎng)格節(jié)點數(shù)約20萬，采用64位操作系統(tǒng)，Intel Xeon 2.8 GHz處理器，16 GB內(nèi)存的小型工作站進行計算，單個聲傳播算例耗時約5 min。故可知，計算與二維外涵道網(wǎng)格同等網(wǎng)格密度下的單個三維聲傳播算例，耗時約1 h。以30個初始三維樣本算例為例，至少需要1～2 d才能完成一輪優(yōu)化計算。而本文方法僅需計算1項偏導數(shù)，單個三維聲傳播計算時間大約1 min，初始樣本點為120個，完成整個三維聲襯優(yōu)化計算過程，僅耗時兩個多小時，體現(xiàn)了本文發(fā)展的三維聲襯優(yōu)化方法的高效率性。

優(yōu)化歷史如圖9所示，從中可知，前100步優(yōu)化迭代，目標函數(shù)TL值已經(jīng)接近最大值，后面迭代運算中目標函數(shù)的改進較少，最后50步左右收斂到最大TL值。圖10給出了聲襯最優(yōu)化結(jié)果，最優(yōu)聲襯阻值和抗值為Z=1.08-0.91i，在該頻率和傳播模態(tài)下最大可取得23.76 dB的降噪效果。硬壁面構型和含最優(yōu)聲襯構型的管內(nèi)聲傳播計算結(jié)果如圖11所示，可知，風扇后傳噪聲傳播經(jīng)過最優(yōu)聲襯時，聲能量被有效吸收，證實了本文中三維聲襯最優(yōu)阻抗值結(jié)果的有效性。

圖9 傳遞損失值隨優(yōu)化迭代次數(shù)的變化Fig.9 The optimization history of TL objective functions

圖10 聲襯最優(yōu)阻值和抗值云圖Fig.10 The optimized impedance values of acoustic liner

圖11 管內(nèi)聲傳播結(jié)果對比Fig.11 Comparison of sound pressure in bypass duct

4 結(jié)論

由于航空發(fā)動機外涵道流體和外部流之間出現(xiàn)了剪切層，導致風扇后傳噪聲模擬變得復雜。此外，可抑制風扇后傳噪聲的三維降噪聲襯設計計算開銷巨大，進一步增加了其降噪聲襯設計的難度。通過本文研究工作，得到以下結(jié)論。

(1)基于線化歐拉方程，通過數(shù)學變換，將線化歐拉方程簡化成新的拋物線近似方程，建立了聲源邊界條件和聲襯阻抗邊界條件，合理地把雙向聲傳播的耦合問題，近似成單向聲傳播問題。最后，采用高階DRP空間離散和LDDRK時間推進格式數(shù)值求解拋物線方程，兼顧精度和計算效率問題，發(fā)展得到高效的風扇后傳噪聲聲傳播預測方法。

(2)基于拋物線近似技術和Kriging代理模型技術，發(fā)展了可抑制風扇后傳噪聲的三維降噪聲襯設計方法和流程。建立了代理模型進一步降低優(yōu)化尋優(yōu)開銷，實現(xiàn)了三維外涵道聲襯快速優(yōu)化設計。

(3)將該發(fā)展的方法應用于設計典型外涵道降噪聲襯。最優(yōu)聲襯阻值和抗值為Z=1.08-0.91i，在該頻率和傳播模態(tài)下，最大可取得23.76 dB 的降噪效果，設計出了一個全局最優(yōu)的三維降噪聲襯，驗證了本文中三維外涵道降噪聲襯的高效優(yōu)化設計框架的有效性和高效率。