盧蕊, 丁永樂, 余培汛, 譚坤, 潘光

(1．西北工業(yè)大學(xué)航海學(xué)院， 陜西 西安 710072； 2．西安精密機(jī)械研究所， 陜西 西安 710077； 3．西北工業(yè)大學(xué)航空學(xué)院， 陜西 西安 710072)

隨著計(jì)算機(jī)科學(xué)和計(jì)算流體力學(xué)的快速發(fā)展,數(shù)值模擬有效彌補(bǔ)了風(fēng)洞試驗(yàn)的不足,已成為流致噪聲研究不可替代的重要手段。數(shù)值模擬不僅可以精細(xì)化模擬聲源,還可以分析噪聲傳播特性。這有助于對(duì)噪聲的產(chǎn)生、傳播和反射等物理機(jī)制進(jìn)行更加深入的研究。流致噪聲產(chǎn)生的物理本質(zhì)是流動(dòng)引起的壓力脈動(dòng),尾緣流致噪聲包括復(fù)雜的湍流噪聲和流體與固邊界之間相互作用的噪聲。因此,流致噪聲的數(shù)值計(jì)算精度依賴于聲源產(chǎn)生與傳播的求解精度,需采用低耗散、低色散數(shù)值方法進(jìn)行精細(xì)化模擬[1-2]。

在流致噪聲數(shù)值計(jì)算方面,目前主要采用2種數(shù)值方法[3]進(jìn)行模擬分析:直接數(shù)值方法和混合數(shù)值方法。

直接數(shù)值模擬方法的基本思路是將流場(chǎng)與聲場(chǎng)統(tǒng)一處理,通過求解可壓縮納維-斯托克斯(Navier-Stokes,NS)方程來模擬聲源的產(chǎn)生以及傳播過程。雖然該方法在理論上可獲得很高的模擬精度,但其仍然存在一些缺陷。例如:聲學(xué)擾動(dòng)易淹沒在數(shù)值誤差中,整個(gè)計(jì)算域均需采用高精度、低耗散低色散數(shù)值方法和準(zhǔn)確的邊界條件。相比于直接數(shù)值模擬方法,綜合計(jì)算資源與工程應(yīng)用價(jià)值等方面的考慮,混合數(shù)值方法是一種處理復(fù)雜構(gòu)型聲學(xué)問題合理有效的技術(shù)手段。該方法采用一種分區(qū)求解的策略,將整個(gè)物理過程劃分為2個(gè)或3個(gè)計(jì)算區(qū)域[4-5]:聲源區(qū)域,主要是復(fù)雜流動(dòng)變化引起的流場(chǎng)脈動(dòng);聲傳播區(qū)域,主要是描述聲源的傳播過程,可考慮聲波的反射、衍射、折射等現(xiàn)象。

聲源區(qū)域聲源求解通常采用諸如大渦模擬(large eddy simulation,LES)、RANS/LES模擬等[6-7]方法獲得,其難點(diǎn)在于湍流的高精度數(shù)值模擬。此外,在一些寬頻噪聲預(yù)測(cè)問題中,也可通過求解RANS(Reynolds average NS equations)獲得湍流平均信息生成噪聲源項(xiàng),類似于SNGR(stochastic noise generation and radiation)、RPM(random particle mesh)、FRPM(fast random particle mesh)等[7-9]湍流統(tǒng)計(jì)聲源法。湍流聲源統(tǒng)計(jì)方法主要依賴于湍流自身的統(tǒng)計(jì)信息,相比于RANS/LES、LES、DNS等方法,其計(jì)算效率更高,目前已在翼型后緣湍流噪聲、增升裝置縫翼湍流噪聲、發(fā)動(dòng)機(jī)噴流噪聲、燃燒噪聲等問題[10-12]中廣泛應(yīng)用。然而,聲源統(tǒng)計(jì)方法在噪聲產(chǎn)生機(jī)理方面以及周期性聲學(xué)等問題方面有所欠缺。

聲傳播區(qū)域通?？赏ㄟ^求解諸如聲擾動(dòng)方程[11]、線性化歐拉方程[12]、非線性擾動(dòng)方程[13]等控制方程實(shí)現(xiàn)。聲傳播區(qū)域的難點(diǎn)在于如何準(zhǔn)確地反映物理波的運(yùn)動(dòng)(其中包括渦、聲、熱等模態(tài))和傳播特征(包括色散、耗散、方向、相速度、群速度等)[14-15]。其次,為了準(zhǔn)確模擬上述物理波的特性,所使用的離散格式除了具有相容、穩(wěn)定、高精度的特點(diǎn),還應(yīng)該盡可能減小聲場(chǎng)相關(guān)物理量過多的耗散,并保持各向同性。最后,在聲學(xué)遠(yuǎn)場(chǎng)邊界,為了避免擾動(dòng)量產(chǎn)生反射后污染計(jì)算域,需要采用數(shù)值吸收邊界條件(也稱為輻射邊界,無反射邊界等)[16]。這種數(shù)值邊界起2個(gè)作用:①讓計(jì)算區(qū)域內(nèi)的各種波(聲波,渦波,熵波)能無反射地通過邊界;②削弱外界擾動(dòng)污染計(jì)算域。到目前為止,為CFD(computational fluid dynamics)而發(fā)展的離散格式、邊界條件大部分無法直接應(yīng)用于聲傳播方程中。因此,發(fā)展高精度數(shù)值方法越來越重要。

對(duì)于混合流致噪聲計(jì)算方法,另一需要考慮的問題是如何耦合聲源計(jì)算求解與聲傳播求解?，F(xiàn)有的方法主要分為等效聲源法和聲波邊界條件法[17]。等效聲源法是指在傳播域中通過對(duì)聲波傳播方程源項(xiàng)求解直接獲得,該源項(xiàng)常被認(rèn)為是噪聲的主要來源。源項(xiàng)中的湍流速度可通過大渦模擬等高精度數(shù)值方法直接獲得,或隨機(jī)聲源模型構(gòu)造得到。聲波邊界條件法是指將聲源區(qū)域的邊界作為聲傳播方程的邊界條件。相比于等效聲源法,聲波邊界法應(yīng)用范圍更廣,可應(yīng)用于運(yùn)動(dòng)物體的聲學(xué)問題中。

本文闡述了聲擾動(dòng)方程的數(shù)值求解方法(空間離散格式、時(shí)間推進(jìn)格式、邊界條件)、聲源數(shù)據(jù)傳遞方法等?；诼暡ㄟ吔鐥l件思想,改善其使用策略以達(dá)到削弱界面處數(shù)值偽波問題,從而建立基于聲擾動(dòng)方程和大渦模擬的高精度CAA(computational aero-acoustics)混合數(shù)值方法。以NACA0012翼型為研究對(duì)象,開展翼型尾緣聲源的輻射特性分析研究,通過與試驗(yàn)數(shù)據(jù)、LES直接計(jì)算結(jié)果對(duì)比分析,驗(yàn)證混合流致噪聲預(yù)測(cè)方法的可行性及準(zhǔn)確度。

1 聲擾動(dòng)方程數(shù)值方法

1.1 聲擾動(dòng)方程

為了得到直接求解聲學(xué)擾動(dòng)的控制方程,可以將擾動(dòng)形式NS方程轉(zhuǎn)化到頻率-波數(shù)域,對(duì)其本征模態(tài)進(jìn)行分離,忽略擾動(dòng)量的黏性效應(yīng),保留僅會(huì)激發(fā)聲學(xué)模態(tài)響應(yīng)的部分源項(xiàng),最后,通過反變換得到時(shí)域下的聲擾動(dòng)方程[11]。

1.2 空間離散方法

對(duì)于聲波邊界條件方法的應(yīng)用,聲源區(qū)域與傳播區(qū)域交界處往往存在較強(qiáng)的流場(chǎng)間斷問題。為了保證空間離散低耗散、低色散特性,需在交界處增強(qiáng)數(shù)值計(jì)算的穩(wěn)定性。本文空間離散所采用的策略:①內(nèi)部網(wǎng)格通量求解采用四階精度的DRP有限差分格式[15];②交界面處通量求解采用WENO改進(jìn)的有限差分格式,并選用混合形式的Lax-Friedrichs流通矢量分裂方法進(jìn)行通量分裂。其中WENO改進(jìn)格式具體如下:

根據(jù)Tam所提出的DRP思想,對(duì)現(xiàn)有的WENO格式進(jìn)行改善。考慮對(duì)流方程?q/?t+?f(q)/?x=0,等間距的網(wǎng)格分布xi=iΔx(i=0,…,N),其半離散格式可表示為

(3)

其中變量常數(shù)A=?f/?q。

對(duì)于5點(diǎn)插值模板,WENO格式的數(shù)值通量可通過加權(quán)的形式進(jìn)行表示

(4)

將(4)式代入(1)式通量項(xiàng)中,并對(duì)傅里葉變換后的方程右端進(jìn)行泰勒級(jí)數(shù)展開,選取五階精度,整理可得

(5)

以φ0作為自由變量,通過迎風(fēng)形式的積分誤差來進(jìn)行優(yōu)化,積分誤差E定義如下

(6)

求(6)式的誤差函數(shù)的極小值,可轉(zhuǎn)化為?E/?φ0=0的根。表1給出文中所選用的幾組優(yōu)化系數(shù)φ0的具體取值。

表1 優(yōu)化系數(shù)φ0的取值

圖1給出了幾種不同WENO格式及DRP格式的相位誤差特性。從圖中可以看出,DRP對(duì)聲波分辨率最高,所需網(wǎng)格點(diǎn)數(shù)最小。然而,在交界面處極易出現(xiàn)網(wǎng)格及流場(chǎng)間斷現(xiàn)象,為增加該處的穩(wěn)定性,

圖1 不同WENO格式的相位誤差特性

邊界處的離散格式選用了色散特性較差但穩(wěn)定性較高的WENO-opt3格式。

1.3 時(shí)間推進(jìn)方法

在APE方程的時(shí)間離散方面,采用了Vasanth等人所提出的大時(shí)間步長(zhǎng)HALE-RK64格式[18]。以一維偏微分方程為例

(7)

HALE-RK64格式可表示為:

(8)

(8)式中的參數(shù)具體表達(dá)式詳見文獻(xiàn)[18]。

為了說明HALE-RK64時(shí)間離散格式在大時(shí)間步長(zhǎng)下具有較好的穩(wěn)定性,根據(jù)穩(wěn)定性特征公式(9),圖2給出了3種不同時(shí)間離散格式的穩(wěn)定性分析。從圖中可以明顯看出:HALE-RK64格式的穩(wěn)定性區(qū)域最大,而相同階數(shù)的Hu-RK64格式[19]穩(wěn)定性區(qū)域最小。

圖2 不同時(shí)間離散格式的穩(wěn)定性分析

(9)

1.4 無反射邊界條件

無反射邊界條件被人為地劃分為計(jì)算區(qū)域和吸收層,其作用是使得計(jì)算區(qū)域內(nèi)的各種波(聲波、渦波、熵波)能無反射地通過邊界;同時(shí)讓外界的任何擾動(dòng)都不能傳入計(jì)算區(qū)域。已有的無反射邊界條件大致可分為基于特征變量、漸近解、吸收等理念的邊界條件。而理想匹配層邊界(perfectly matched layer,PML)是一種聲反射系數(shù)最小的遠(yuǎn)場(chǎng)邊界條件。

本文采用了無分裂形式的理想匹配層邊界(PML)條件[20],其二維表達(dá)式具體如下

(10)

為了穩(wěn)定地離散PML邊界方程,本文采用了Gao等[14]提出的DRP7544單側(cè)差分格式,其表達(dá)式具體如下

(11)

其穩(wěn)定性的驗(yàn)證分析工作已在文獻(xiàn)[9,15]中給出。

2 LES數(shù)值方法

本文所采用的LES方法計(jì)算的控制方程具體如下

(12)

動(dòng)量方程中亞格子應(yīng)力采用WALE模型進(jìn)行計(jì)算,具體表達(dá)式如下

(13)

式中

3 聲源數(shù)據(jù)傳遞算法

如何快速有效地實(shí)現(xiàn)數(shù)據(jù)從CFD網(wǎng)格到CAA網(wǎng)格的搜索插值,是混合方法預(yù)測(cè)噪聲的一個(gè)重要環(huán)節(jié)。高效網(wǎng)格點(diǎn)搜索和流場(chǎng)插值問題在CAA混合方法中尤為重要,因?yàn)槊窟M(jìn)行一次聲場(chǎng)的時(shí)間推進(jìn)均涉及到耦合層的一次插值。本文使用了一種KD樹數(shù)據(jù)結(jié)構(gòu)進(jìn)行搜索。

這種方法除了滿足全局和高效外,還不要求輸入網(wǎng)格的結(jié)構(gòu)性。具體算法主要包括KD樹的構(gòu)造和搜索。以二維平面上的點(diǎn)為例,介紹基于KD樹結(jié)構(gòu)的搜索方法。二維平面上有6個(gè)點(diǎn),如圖3所示。根據(jù)x,y坐標(biāo)點(diǎn)的方差統(tǒng)計(jì),x方向方差最大,沿x軸分為2個(gè)子區(qū)域。對(duì)左右子空間重復(fù)上述步驟,直到所有子空間中只有一個(gè)坐標(biāo)點(diǎn)形成葉節(jié)點(diǎn)。對(duì)應(yīng)的空間劃分如圖4所示。

圖3 KD樹算法說明

圖4 空間劃分示意圖

雖然該方法在搜索上是有效的,但在節(jié)點(diǎn)的構(gòu)造上卻比較緩慢。在實(shí)際應(yīng)用中,采用場(chǎng)景選擇算法來平衡樹的構(gòu)建和搜索樹的復(fù)雜度。以二維平面中的點(diǎn)(2,4.5)為例,介紹了二叉樹中最近點(diǎn)的查找算法。

通過二叉樹搜索,找到子空間。從二叉樹的根(7,2)開始,因?yàn)?<7,移動(dòng)到左側(cè)子節(jié)點(diǎn)(5,4)。因?yàn)?.5>4,移動(dòng)到右子節(jié)點(diǎn)(4,7)。按此方法移動(dòng),直到當(dāng)前點(diǎn)是葉節(jié)點(diǎn),停止搜索。

回溯以找到最近的點(diǎn)。首先計(jì)算目標(biāo)點(diǎn)與葉節(jié)點(diǎn)(4,7)之間的距離。然后,沿著搜索路徑追溯到父節(jié)點(diǎn)(5,4),并計(jì)算其到目標(biāo)點(diǎn)的距離r。當(dāng)發(fā)現(xiàn)目標(biāo)點(diǎn)靠近父節(jié)點(diǎn)時(shí),以目標(biāo)點(diǎn)為中心,以距離r為半徑,生成圓,如圖5所示。圓與當(dāng)前節(jié)點(diǎn)(4.5,4)另一側(cè)的子空間相交。而另一側(cè)的子空間中可能會(huì)出現(xiàn)最近點(diǎn),因此需要在另一側(cè)的子空間上重復(fù)使用二進(jìn)制搜索和回溯,直到圓不相交于兩邊的子空間。此時(shí),找到的節(jié)點(diǎn)為目標(biāo)點(diǎn)的最近點(diǎn)。

圖5 KD樹搜索示意圖

通過以上分析可以看出,得益于二叉樹的數(shù)據(jù)結(jié)構(gòu),使用KD樹搜索的復(fù)雜度為O(lg(n)),具有較高的搜索效率。在完成CAA網(wǎng)格點(diǎn)的搜索后,將采用基于形函數(shù)的拉格朗日插值方法進(jìn)行CFD信息的插值。

4 基于聲波邊界條件的混合方法

4.1 混合方法的實(shí)施策略

聲波邊界條件法是指將聲源區(qū)域的外邊界作為聲傳播方程的入口邊界條件,具體如圖6所示。在聲源求解和聲傳播的耦合過程中,CFD求解和CAA求解所采用的網(wǎng)格尺度及分布形態(tài)有所區(qū)別,導(dǎo)致兩區(qū)域間存在數(shù)據(jù)插值、時(shí)間步長(zhǎng)匹配、流場(chǎng)間斷等問題。流場(chǎng)間斷引起的數(shù)值偽波問題將在5.3節(jié)進(jìn)行分析及改善。對(duì)于時(shí)間步長(zhǎng)的匹配問題:CFD計(jì)算求解和CAA計(jì)算求解時(shí)間步長(zhǎng)分別為Δt1和Δt2,①當(dāng)Δt1≤Δt2時(shí),CAA聲波邊界區(qū)域的擾動(dòng)值將在一個(gè)CFD計(jì)算時(shí)間間隔中線性插值獲得;②當(dāng)Δt1>Δt2時(shí),CAA聲波邊界區(qū)域的擾動(dòng)值將在3個(gè)CFD計(jì)算時(shí)間間隔中二次曲線插值獲得。

圖6 聲波邊界條件法的耦合說明

4.2 混合方法的數(shù)值驗(yàn)證

1) CFD求解器的數(shù)值驗(yàn)證

下面以“槽道中心阻塞率20%的方柱”為研究對(duì)象,采用LES方法進(jìn)行流場(chǎng)計(jì)算分析,驗(yàn)證本文所使用CFD求解器的可靠性。數(shù)值模擬的計(jì)算域及計(jì)算條件與Nakagawa等人的試驗(yàn)相似,方柱位于槽道的中央。

在本算例中,基于方柱直徑D和來流速度Uinf的雷諾數(shù)Re為3 000,流場(chǎng)域的展向尺寸為3D,在入口處使用了完全一致的來流速度Uinf=0.2。在此流動(dòng)中,柱體下游卡門渦街的發(fā)展將受到壁面的嚴(yán)重影響,而近壁面流動(dòng)結(jié)構(gòu)也會(huì)受到卡門渦街的影響。

圖7給出了Q=0.05的等值面圖,從圖中可以看出方柱后面脫落出復(fù)雜的渦系。

圖7 Q=0.05等值面圖(采用流向速度守恒變量著色)

圖8給出了x/D=1.0,3.5站位處的雷諾應(yīng)力曲線對(duì)比結(jié)構(gòu),從圖中可以看出,相比于DES方法,LES方法與試驗(yàn)數(shù)據(jù)更吻合。

圖8 雷諾應(yīng)力分布圖(x/D=1.0,3.5)

2) CAA求解器的數(shù)值驗(yàn)證

本算例為第四屆CAA基準(zhǔn)算例研討會(huì)[1]中的二維聲波穿過剪切層后的輻射與折射現(xiàn)象,用來驗(yàn)證擾動(dòng)方程的數(shù)值特性與邊界條件的正確性。對(duì)于平行射流剪切層問題,其平均流中包含了黏性和非線性效應(yīng)。

計(jì)算域?yàn)閇-50,150]m×[0,50]m的矩形區(qū)域。在計(jì)算域下方設(shè)置對(duì)稱面邊界,其余3個(gè)方向設(shè)置PML邊界條件。網(wǎng)格內(nèi)點(diǎn)與內(nèi)邊界處空間離散都使用了7點(diǎn)DRP-WENO格式,時(shí)間推進(jìn)使用HALE-RK64格式,無量綱時(shí)間步長(zhǎng)為Δt=0.2。

圖9 射流剪切層問題的計(jì)算條件示意

圖10為t=12T時(shí)刻的擾動(dòng)壓力場(chǎng)模擬結(jié)果,T為聲源周期。為更細(xì)致比較各狀態(tài)間的差異,圖11給出了t=12T時(shí)刻y=15上的擾動(dòng)壓力結(jié)果對(duì)比,其壓力曲線計(jì)算解與解析解基本吻合,驗(yàn)證了聲傳播控制方程數(shù)值離散方法的準(zhǔn)確性與可行性。

圖10 t=12T時(shí)刻的射流問題擾動(dòng)壓力分布等值線云圖

圖11 t=12T時(shí)刻y=15的擾動(dòng)壓力曲線

5 翼型流致噪聲特性

5.1 計(jì)算設(shè)置

為了開展翼型的聲輻射特性研究,選取了BANC(benchmark problems for airframe noise computations)試驗(yàn)中的NACA0012[21]作為研究對(duì)象。NACA0012翼型試驗(yàn)的參考弦長(zhǎng)為0.4 m,單位尺度雷諾數(shù)為3.83×106,轉(zhuǎn)捩位置、來流速度、迎角分別如表2所示,其中SS為翼型吸力面轉(zhuǎn)捩位置,PS為翼型壓力面轉(zhuǎn)捩位置。計(jì)算條件與試驗(yàn)條件保持一致。

表2 BANC試驗(yàn)的測(cè)試條件

CFD網(wǎng)格及CAA網(wǎng)格分布如圖12所示,其中CFD網(wǎng)格量為80萬左右,為了有效模擬附面層的流動(dòng)特征,其網(wǎng)格尺度分布保證了y+<1的數(shù)值計(jì)算要求,此外在翼型尾緣分離流區(qū)域進(jìn)行了局部網(wǎng)格加密處理。而CAA網(wǎng)格在附面層附近則更加注重網(wǎng)格的正交性及長(zhǎng)細(xì)比分布要求,整體網(wǎng)格的分布尺度滿足一個(gè)波長(zhǎng)所需的網(wǎng)格點(diǎn)數(shù)需求,網(wǎng)格量為100萬左右。

圖12 計(jì)算網(wǎng)格分布

5.2 CFD計(jì)算結(jié)果

為了定量對(duì)比數(shù)值計(jì)算結(jié)果的準(zhǔn)確性,圖13給出了翼型表面壓力系數(shù)分布對(duì)比曲線圖。圖中紅色線為本文計(jì)算結(jié)果,黑色實(shí)心圓點(diǎn)為試驗(yàn)測(cè)量結(jié)果[21],藍(lán)色虛線為Xfoil軟件計(jì)算結(jié)果,綠色直線位置為試驗(yàn)固定轉(zhuǎn)捩點(diǎn)位置。從全局圖中可看出:本文所計(jì)算的壓力系數(shù)分布可看出轉(zhuǎn)捩位置附近有明顯的振蕩現(xiàn)象。

圖13 壓力系數(shù)分布對(duì)比

當(dāng)獲得定常狀態(tài)解后,將其作為初始流場(chǎng)進(jìn)行LES非定常計(jì)算,其無量綱時(shí)間步長(zhǎng)為0.005。待計(jì)算收斂后,進(jìn)行非定常源項(xiàng)變量的采樣,采樣步數(shù)30 000。

圖14為非定常計(jì)算求解獲得的流向脈動(dòng)速度分布云圖,圖15為法向脈動(dòng)速度分布云圖。其中,CFD脈動(dòng)速度變量通過瞬時(shí)流場(chǎng)速度變量減去時(shí)均流場(chǎng)速度變量獲得。從這兩幅云圖中可看出:①流向和法向的脈動(dòng)速度主要分布于翼型的尾緣及尾流區(qū);②流向脈動(dòng)速度在尾緣上下端面呈近似對(duì)稱形態(tài),尾流區(qū)呈正負(fù)交叉形態(tài);③法向脈動(dòng)速度在尾緣上下端面呈近似正負(fù)交替的對(duì)稱形態(tài),尾流區(qū)呈正負(fù)交替形態(tài)。④法向脈動(dòng)速度的強(qiáng)度明顯大于流向脈動(dòng)速度,這與升阻力系數(shù)曲線的幅值大小表現(xiàn)一致。

圖14 流向脈動(dòng)速度分布 圖15 法向脈動(dòng)速度分布

5.3 聲傳播區(qū)域計(jì)算結(jié)果分析

基于CFD/CAA混合方法是將整個(gè)計(jì)算域分成聲源產(chǎn)生區(qū)域和聲傳播區(qū)域2個(gè)部分。其中聲源區(qū)計(jì)算采用LES數(shù)值方法進(jìn)行計(jì)算求解,傳播區(qū)采用聲擾動(dòng)方程求解,計(jì)算時(shí)間步長(zhǎng)與LES非定常計(jì)算步長(zhǎng)保持一致,初始的聲源區(qū)域選取如圖16所示。

圖16 初始的聲源區(qū)域選取方式

初始的聲源區(qū)選取范圍主要依靠Lamb矢量分布范圍來決定,本文選取了包含90%Lamb矢量大小的區(qū)域。兩計(jì)算區(qū)域的耦合采用了聲波邊界條件法。

圖17為采用圖16聲源區(qū)域選取方式所獲得的聲壓云圖, 從圖中可以看出存在大量的數(shù)值偽波現(xiàn)象,這主要是由聲波穿過強(qiáng)流場(chǎng)間斷所引起的不穩(wěn)定擴(kuò)散導(dǎo)致[22]。

圖17 CFD/CAA求解的聲壓云圖

為解決上述現(xiàn)象,文中主要采取兩方面的措施:一方面將聲源區(qū)邊界范圍擴(kuò)大至95%Lamb矢量大小區(qū)域,減小聲源區(qū)邊界穿越尾流分離區(qū);另一方面,為了減小聲源區(qū)與傳播區(qū)因聲源變量間斷引起的短波發(fā)散問題,可在初始步賦值時(shí),添加人工過渡層,如圖18所示。其中,灰色區(qū)域?yàn)槁晜鞑^(qū)域,綠色區(qū)域?yàn)槁曉磪^(qū)域,紅色區(qū)域?yàn)檫^渡區(qū)域。過渡區(qū)域采用的函數(shù)是以距離為變量的半高斯函數(shù)。

圖18 改進(jìn)的聲源區(qū)域選取方式

圖19顯示了某時(shí)刻瞬時(shí)聲壓分布云圖,圖19a)為采用項(xiàng)目組自編程序LES直接模擬求解所獲得的聲場(chǎng)分布,圖19b)為基于聲波邊界條件的CFD/CAA混合方法計(jì)算結(jié)果。從聲壓云圖中可以看出,2種方法所獲得聲壓云圖分布形態(tài)及大小基本一致,均呈現(xiàn)沿y=0的對(duì)稱分布形態(tài)。此外,從圖19a)中可看出:翼型的尾緣噪聲主要是由尾緣處上下翼面正負(fù)交替的渦作用產(chǎn)生。

圖19 瞬時(shí)時(shí)刻聲壓云圖

圖20為分別采用LES直接計(jì)算與CFD/CAA方法計(jì)算的噪聲指向性對(duì)比圖。由圖可以看出,采用全場(chǎng)LES直接計(jì)算和CFD/CAA混合方法得到的指向性曲線基本重合,誤差不超過2 dB,聲波主要強(qiáng)度位于1 200到1 500區(qū)域。造成兩者誤差的可能原因:傳播區(qū)域LES方法和CFD/CAA方法的網(wǎng)格分布、數(shù)值離散方法均有所區(qū)別,此外聲源區(qū)的位置選取也可能有一定影響。總的來說,該結(jié)果表明基于聲波邊界條件的流致噪聲混合方法可準(zhǔn)確預(yù)測(cè)尾緣噪聲問題。圖21給出了位于翼型后緣正上方(900)1.5倍翼型參考弦長(zhǎng)處的聲壓級(jí)頻譜曲線對(duì)比結(jié)果。通過對(duì)比可看出:兩者計(jì)算結(jié)果趨勢(shì)和聲壓級(jí)量級(jí)基本相當(dāng),最大聲壓級(jí)相差1 dB左右。

圖20 噪聲指向性對(duì)比

圖21 1/3倍頻聲壓級(jí)頻譜曲線

6 結(jié) 論

本文開展了基于聲擾動(dòng)方程的流致噪聲混合方法的研究分析,并對(duì)混合方法中的各模塊進(jìn)行了驗(yàn)證分析。在此基礎(chǔ)上,以二維NACA0012翼型為研究對(duì)象,分析了其聲學(xué)輻射特性,具體結(jié)論如下:

1) 本文所建立的基于聲擾動(dòng)方程的混合預(yù)測(cè)方法,能夠高精度地預(yù)測(cè)流致噪聲問題中的聲輻射特征,并直觀反映聲波與聲波之間的相互作用現(xiàn)象及空間分布特征。

2) 對(duì)于本文所提及的CFD/CAA混合方法,聲源區(qū)的邊界應(yīng)避免位于尾流的強(qiáng)擾動(dòng)中。此外,為了減小聲源區(qū)與傳播區(qū)間的聲波間斷問題,初始步可通過引入阻尼層達(dá)到順利過渡的作用。

3) 從翼型尾緣的噪聲源地分布特點(diǎn)可看出,聲源的分布位置與旋渦的分布位置基本保持一致。翼型的尾緣噪聲主要是由尾緣處上下翼面正負(fù)交替的渦作用產(chǎn)生。

4) 在聲波的相互作用及與壁面的作用下,翼型噪聲的指向性主要呈現(xiàn)偶極子特征,主要強(qiáng)度位于120°～150°區(qū)域。

5) 基于聲波邊界條件方法的高精度流致噪聲計(jì)算方法由于聲源求解和聲傳播求解區(qū)域部分重合特性,后續(xù)有望將該方法應(yīng)用于類似撲翼、螺旋槳等噪聲問題中。