基于CFD/CAA耦合邊界方法的翼型尾緣噪聲預測

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

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

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

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

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

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

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

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

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

1 聲擾動方程數值方法

1.1 聲擾動方程

為了得到直接求解聲學擾動的控制方程,可以將擾動形式NS方程轉化到頻率-波數域,對其本征模態進行分離,忽略擾動量的黏性效應,保留僅會激發聲學模態響應的部分源項,最后,通過反變換得到時域下的聲擾動方程[11]。

1.2 空間離散方法

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

根據Tam所提出的DRP思想,對現有的WENO格式進行改善?？紤]對流方程?q/?t+?f(q)/?x=0,等間距的網格分布xi=iΔx(i=0,…,N),其半離散格式可表示為

(3)

其中變量常數A=?f/?q。

對于5點插值模板,WENO格式的數值通量可通過加權的形式進行表示

(4)

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

(5)

以φ0作為自由變量,通過迎風形式的積分誤差來進行優化,積分誤差E定義如下

(6)

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

表1 優化系數φ0的取值

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

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

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

1.3 時間推進方法

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

(7)

HALE-RK64格式可表示為:

(8)

(8)式中的參數具體表達式詳見文獻[18]。

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

圖2 不同時間離散格式的穩定性分析

(9)

1.4 無反射邊界條件

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

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

(10)

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

(11)

其穩定性的驗證分析工作已在文獻[9,15]中給出。

2 LES數值方法

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

(12)

動量方程中亞格子應力采用WALE模型進行計算,具體表達式如下

(13)

式中

3 聲源數據傳遞算法

如何快速有效地實現數據從CFD網格到CAA網格的搜索插值,是混合方法預測噪聲的一個重要環節。高效網格點搜索和流場插值問題在CAA混合方法中尤為重要,因為每進行一次聲場的時間推進均涉及到耦合層的一次插值。本文使用了一種KD樹數據結構進行搜索。

這種方法除了滿足全局和高效外,還不要求輸入網格的結構性。具體算法主要包括KD樹的構造和搜索。以二維平面上的點為例,介紹基于KD樹結構的搜索方法。二維平面上有6個點,如圖3所示。根據x,y坐標點的方差統計,x方向方差最大,沿x軸分為2個子區域。對左右子空間重復上述步驟,直到所有子空間中只有一個坐標點形成葉節點。對應的空間劃分如圖4所示。

圖3 KD樹算法說明

圖4 空間劃分示意圖

雖然該方法在搜索上是有效的,但在節點的構造上卻比較緩慢。在實際應用中,采用場景選擇算法來平衡樹的構建和搜索樹的復雜度。以二維平面中的點(2,4.5)為例,介紹了二叉樹中最近點的查找算法。

通過二叉樹搜索,找到子空間。從二叉樹的根(7,2)開始,因為2<7,移動到左側子節點(5,4)。因為4.5>4,移動到右子節點(4,7)。按此方法移動,直到當前點是葉節點,停止搜索。

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

圖5 KD樹搜索示意圖

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

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

4.1 混合方法的實施策略

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

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

4.2 混合方法的數值驗證

1) CFD求解器的數值驗證

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

在本算例中,基于方柱直徑D和來流速度Uinf的雷諾數Re為3 000,流場域的展向尺寸為3D,在入口處使用了完全一致的來流速度Uinf=0.2。在此流動中,柱體下游卡門渦街的發展將受到壁面的嚴重影響,而近壁面流動結構也會受到卡門渦街的影響。

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

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

圖8給出了x/D=1.0,3.5站位處的雷諾應力曲線對比結構,從圖中可以看出,相比于DES方法,LES方法與試驗數據更吻合。

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

2) CAA求解器的數值驗證

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

計算域為[-50,150]m×[0,50]m的矩形區域。在計算域下方設置對稱面邊界,其余3個方向設置PML邊界條件。網格內點與內邊界處空間離散都使用了7點DRP-WENO格式,時間推進使用HALE-RK64格式,無量綱時間步長為Δt=0.2。

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

圖10為t=12T時刻的擾動壓力場模擬結果,T為聲源周期。為更細致比較各狀態間的差異,圖11給出了t=12T時刻y=15上的擾動壓力結果對比,其壓力曲線計算解與解析解基本吻合,驗證了聲傳播控制方程數值離散方法的準確性與可行性。

圖10 t=12T時刻的射流問題擾動壓力分布等值線云圖

圖11 t=12T時刻y=15的擾動壓力曲線

5 翼型流致噪聲特性

5.1 計算設置

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

表2 BANC試驗的測試條件

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

圖12 計算網格分布

5.2 CFD計算結果

為了定量對比數值計算結果的準確性,圖13給出了翼型表面壓力系數分布對比曲線圖。圖中紅色線為本文計算結果,黑色實心圓點為試驗測量結果[21],藍色虛線為Xfoil軟件計算結果,綠色直線位置為試驗固定轉捩點位置。從全局圖中可看出:本文所計算的壓力系數分布可看出轉捩位置附近有明顯的振蕩現象。

圖13 壓力系數分布對比

當獲得定常狀態解后,將其作為初始流場進行LES非定常計算,其無量綱時間步長為0.005。待計算收斂后,進行非定常源項變量的采樣,采樣步數30 000。

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

圖14 流向脈動速度分布 圖15 法向脈動速度分布

5.3 聲傳播區域計算結果分析

基于CFD/CAA混合方法是將整個計算域分成聲源產生區域和聲傳播區域2個部分。其中聲源區計算采用LES數值方法進行計算求解,傳播區采用聲擾動方程求解,計算時間步長與LES非定常計算步長保持一致,初始的聲源區域選取如圖16所示。

圖16 初始的聲源區域選取方式

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

圖17為采用圖16聲源區域選取方式所獲得的聲壓云圖, 從圖中可以看出存在大量的數值偽波現象,這主要是由聲波穿過強流場間斷所引起的不穩定擴散導致[22]。

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

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

圖18 改進的聲源區域選取方式

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

圖19 瞬時時刻聲壓云圖

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

圖20 噪聲指向性對比

圖21 1/3倍頻聲壓級頻譜曲線

6 結 論

本文開展了基于聲擾動方程的流致噪聲混合方法的研究分析,并對混合方法中的各模塊進行了驗證分析。在此基礎上,以二維NACA0012翼型為研究對象,分析了其聲學輻射特性,具體結論如下:

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

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

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

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

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