非緊致邊界氣動噪聲數(shù)值預測方法

宋亞輝，秦 浩，劉秋洪

（1.中國飛行試驗研究院，西安 710089； 2.西北工業(yè)大學 航空學院，西安 710072）

非緊致邊界氣動噪聲數(shù)值預測方法

宋亞輝1，秦 浩1，劉秋洪2

（1.中國飛行試驗研究院，西安 710089； 2.西北工業(yè)大學 航空學院，西安 710072）

通過尋找滿足非緊致邊界的精確格林函數(shù)，發(fā)展了一種基于聲模擬理論的非緊致邊界氣動噪聲數(shù)值預測方法。該方法首先采用邊界積分方法計算滿足相應邊界條件的精確格林函數(shù)，然后采用精確格林函數(shù)求解FW-H方程，進行遠場氣動噪聲預測?？紤]了非緊致邊界對聲波傳播的散射作用，適用于尺寸較大且?guī)缀瓮庑螐碗s的邊界與非定常流動相互作用誘發(fā)的氣動噪聲的數(shù)值預測。以雷諾數(shù)90 000、馬赫數(shù)0.2的圓柱繞流誘發(fā)的氣動噪聲為例，數(shù)值預測了非緊致邊界條件下的氣動噪聲，并與采用自由空間格林函數(shù)求解FW-H的計算結(jié)果進行了對比，對該方法進行了驗證。

聲學；氣動噪聲；FW-H方程；非緊致邊界；數(shù)值預測

隨著計算機技術(shù)的進步和氣動聲學研究的深入，計算氣動聲學逐漸發(fā)展成為一種探索氣動噪聲機理、尋找噪聲源位置、預測噪聲大小的重要工具。計算流體力學(CFD)和聲模擬理論[1](一般采用FWH方程[2]和自由空間格林函數(shù))相結(jié)合的方法是迄今為止工程應用中最受歡迎的氣動噪聲預測方法[3]。其基本思想是將聲場的計算分解為兩步：首先，近場部分采用CFD技術(shù)進行數(shù)值計算，得到聲源信息數(shù)據(jù)；然后，采用聲模擬理論來求解聲波從近場到遠場部分的傳播。聲模擬理論由Lighthill[1]提出，經(jīng)過研究者[2，4]不斷的發(fā)展，目前被廣泛用于氣動噪聲預測問題研究。Lighthill的聲模擬方程為

式中c0為聲速，p′是聲壓，Tij表示Lighthill湍流應力張量，定義為

式(1)解的通用表達式為[5]

式(3)即是時域的FW-H方程。其中，x和y分別表示觀察點和源點位置，t和τ分別表示觀察點和源點時刻，Ω是產(chǎn)生聲源的流動區(qū)域，S是物體邊界，n是物體邊界上指向流動區(qū)域的法線，ph為單位面積物體邊界作用在流體上的力，ρ0是密度，vn表示物體邊界運動的法向速度，G(x,t|y,τ)是格林函數(shù)，一般采用自由空間格林函數(shù)是Heaviside函數(shù)。

FW-H方程實質(zhì)上是Naiver-Stokes方程的恒等變形，求解氣動噪聲時，理論上要求聲源必須采用滿足低色散、低耗散的高精度方法計算流動[6]。而目前非定常流動多采用低精度的有限容積法求解，由于非定常流動激發(fā)的氣動噪聲在向空間傳播時會在物體邊界上發(fā)生反射和散射，采用自由空間格林函數(shù)直接求解FW-H方程不能完全捕捉邊界對聲波傳播的影響，因而在流場求解精度不足夠高的情況下，要求物體邊界需要滿足聲學緊致條件(即聲波波長遠大于邊界特征尺寸)。對于外形尺寸較大的物體或關(guān)注的聲波頻率較高時，邊界不再滿足聲學緊致條件，直接應用FW-H方程會產(chǎn)生較大誤差[7]。

針對非緊致邊界問題，大量學者進行了研究，主要有以下幾種方法：

(1)基于Lighthill方程變分形式的有限元法[8]，能模擬聲波從近場到遠場的傳播過程，需要對整個求解域進行離散求解，當進行遠場噪聲計算時，計算量巨大，限制了在工程中的應用；

(2)基于自由空間格林函數(shù)的邊界積分法[9]，該方法先求物體邊界處的散射聲壓，然后統(tǒng)一計算遠場的噪聲，計算量小而適用于工程應用，但當前的研究一般適用于不可壓流動誘發(fā)的氣動噪聲預測；

(3)尋找滿足邊界散射條件的精確格林函數(shù)(Exact/Tailored Green’s Function)法[10-12]，該方法不需要對流動進行高精度計算，簡明且易于實現(xiàn)數(shù)值計算。規(guī)則空間的精確格林函數(shù)可以通過理論解析法[11，13]求得，而任意空間的精確格林函數(shù)則需要采用數(shù)值方法進行求解。Takaishi等[12]基于邊界元思想提出了一種求解精確格林函數(shù)的邊界積分方法，適用于任意物體邊界形狀，但Takaishi等只計算了精確格林函數(shù)的1階偏導數(shù)，且只將該方法應用于渦聲理論。本文在文獻[13]中對Takaishi等的數(shù)值計算方法做了進一步研究，求解了精確格林函數(shù)的2階偏導數(shù)，并通過理論解進行了驗證。

在文獻[13]的基礎(chǔ)上，采用數(shù)值方法得到的精確格林函數(shù)來求解FW-H方程，提出了非緊致邊界氣動噪聲數(shù)值預測方法。隨后，數(shù)值預測了圓柱繞流誘發(fā)的氣動噪聲，對提出的方法進行了驗證。

1 數(shù)值計算方法

1.1 非緊致邊界氣動噪聲數(shù)值計算方法

實質(zhì)上，求解聲學波動方程的格林函數(shù)滿足Helmholtz方程，而頻域的自由空間格林函數(shù)是相應的Helmholtz方程的基本解。精確格林函數(shù)的表達式可以應用格林定理演繹相應的Helmholtz方程得到。在頻域下，方程(1)的解也可以寫成通用表達式

其中i為虛數(shù)單位，ω為圓頻率表示頻域下的值，C(f)是立體角函數(shù)，是頻域的精確格林函數(shù)，滿足

對于非緊致邊界，G(x，y，ω)在邊界上滿足

式中n是邊界S的單位法向量，正方向由邊界指向聲源區(qū)域。將式(6)代入式(4)，即可得到非緊致邊界氣動噪聲數(shù)值預測方法的基本方程

假設點z為邊界S上的觀察點，對于精確格林函數(shù)和自由空間格林函數(shù)有

式 (9)乘 以G(x，z，ω)減 去 式 (8)乘 以G0(z，y，ω)，再將等式兩邊同時積分得

根據(jù)格林定理，式(10)等號左端可化為

將式(11)代入式(10)，并利用δ函數(shù)的性質(zhì)，可得精確格林函數(shù)G(x,y,ω)的最終表達式為

式中

式(12)即是精確格林函數(shù)的表達式，等式右邊第一項是自由入射項，意義為當邊界不存在時，觀察點x處的格林函數(shù)值；等式右邊第二項是散射項，意義為當邊界存在時，觀察點x處由邊界的散射產(chǎn)生的格林函數(shù)值。

當物體邊界為靜止剛性邊界時，將式(12)代入式(7)得非緊致邊界氣動噪聲數(shù)值預測方法的基本方程的最終表達式為

式(14)等式右邊第一項為四極子源噪聲，為體積聲源直接輻射到觀察點的噪聲，稱為直接入射噪聲，等式右邊第二項是由于邊界存在產(chǎn)生的噪聲，這部分噪聲包含有邊界與流動相互作用產(chǎn)生的噪聲和聲波傳播到邊界上后散射的噪聲，稱為邊界輻射噪聲。

1.2 精確格林函數(shù)數(shù)值計算方法

采用式(14)計算氣動噪聲關(guān)鍵的一步是求出式(12)所示的精確格林函數(shù)的解。將源點y無限趨近于邊界上的點zˉ，對式(12)有

式中Φ(zˉ)是位于邊界上的點zˉ處的內(nèi)固體角，S′是邊界S除去積分奇點后的部分。將邊界S離散為M個單元，式(15)離散為線性方程組如下

式中E是單位對角矩陣，H的子項滿足下式

其中m=1,2,3,…M，n=1,2,3,…M。

通過數(shù)值方法求解線性方程組(16)，將結(jié)果代入式(12)的等號右端第二項，即可實現(xiàn)對精確格林函數(shù)的數(shù)值計算。

式(14)用到了精確格林函數(shù)的2階偏導數(shù)，下面說明其求解方式。根據(jù)公式(12)，可得精確格林函數(shù)2階偏導數(shù)的表達式為

對于自由空間格林函數(shù)G0(z,y,ω)，其偏導數(shù)存在理論解析解，式(18)右端自由空間格林函數(shù)各偏導數(shù)項均可由理論解析解求出，同時邊界S上各點的G(x,z,ω)值也已經(jīng)由式(15)求出，將這些已知結(jié)果代入式(18)的右端，即可求得精確格林函數(shù)的二階偏導數(shù)值。

2 圓柱非緊致邊界氣動噪聲計算算例

2.1 流動計算結(jié)果及分析

選擇流動計算參數(shù)如下：圓柱直徑D=0.1 m，雷諾數(shù)ReD=90 000，空氣來流馬赫數(shù)Ma=0.2。流場區(qū)域采用分塊結(jié)構(gòu)化網(wǎng)格進行劃分，圓柱壁面第一層網(wǎng)格的厚度為2×10-4D，保證近壁面網(wǎng)格的無量綱距離y+在1的數(shù)量級，尾渦區(qū)域及圓柱附近區(qū)域做了加密處理。流動采用二維非定常不可壓縮流體的雷諾時均Naiver-Stokes方程（RANS）求解，湍流模型采用k-ωSST，時間步長取Δt=5×10-5，保證一個脫落渦周期內(nèi)有不少于100個時間步。

圖1為圓柱壁面升力系數(shù)Cl和阻力系數(shù)Cd隨時間的變化歷程?？梢钥闯觯S圓柱繞流的Cl和Cd呈典型的周期性波動。圖2為對Cl和Cd進行的功率譜分析，Cl的基頻為f0=157.5 Hz，相應的無量綱量斯特勞哈數(shù)St=0.231，Cl明顯的諧波頻率為mf0(m=3，5，7…)，Cd的基頻為2f0，Cd明顯的諧波頻率為nf0(n=2，4，6……)。這一結(jié)果與Takaishi等[12]和Reinaldo等[14]所得結(jié)果一致。

圖3為得到的圓柱壁面時均壓力系數(shù)分布與Reinaldo等的計算結(jié)果及Cantwell實驗結(jié)果[15]的對比?？梢钥闯?，本文結(jié)果與Reinaldo等的結(jié)果基本一致。但與實驗值相比，二者在轉(zhuǎn)捩點附近的壓力系數(shù)都偏小，這是由于采用二維計算模型，忽略了圓柱繞流的三維展向效應所致。本文的重點在于提出和驗證氣動噪聲計算方法，這里不再對流動做進一步的計算和驗證。

圖1 升力系數(shù)和阻力系數(shù)隨時間變化情況

圖2 升力系數(shù)和阻力系數(shù)的功率譜密度

圖3 圓柱壁面平均壓力系數(shù)分布

2.2 氣動噪聲計算結(jié)果及分析

選取26個周期的流場數(shù)據(jù)，總采樣步數(shù)為3 303步，根據(jù)采樣定理，聲場計算頻率范圍選取為50 Hz～2 000 Hz，頻率分辨率為6.1 Hz。選取的聲壓指向性觀察點在以圓柱圓心為原點、半徑為128D的圓上，聲場分布云圖觀察點選取在以圓柱圓心為原點、半徑為20D的圓外。

圓柱空間存在精確格林函數(shù)的理論解析解[13]，可以用來驗證精確格林函數(shù)的數(shù)值計算方法。圖4是采用理論解析方法和本文提出的數(shù)值方法計算精確格林函數(shù)后代入式(14)得到的聲壓指向性圖的對比，可以看出，兩種計算方法的聲壓一致，說明精確格林函數(shù)的數(shù)值計算方法是可靠的。

圖4 數(shù)值結(jié)果和理論解析結(jié)果的聲壓指向性圖

圖5是采用本文數(shù)值計算的精確格林函數(shù)求解FW-H方程得到的聲壓和采用自由空間格林函數(shù)求解FW-H方程得到的聲壓的指向性圖對比。低頻情況下，當頻率為f=f0和f=2f0時(對應的聲波波長分別為21.6D和10.8D)，兩種方法得到的結(jié)果基本一致。隨著頻率增高，當頻率為f=3f0和f=4f0時(對應的聲波波長λ分別為7.2D和5.4D)，兩種方法得到的結(jié)果出現(xiàn)了較為明顯的差異。在低頻情況下，聲波波長遠大于圓柱邊界特征尺寸，圓柱邊界近似可看作滿足聲學緊致條件，噪聲預測方法得到的結(jié)果與采用自由空間格林函數(shù)求解FW-H方程得到的結(jié)果基本一致。隨著頻率的增高，聲波波長和圓柱邊界特征尺寸越來越接近，圓柱邊界已是聲學非緊致邊界，邊界對聲波傳播的影響將變得明顯且不能忽略，本文噪聲預測方法得到的結(jié)果與采用自由空間格林函數(shù)求解FW-H方程得到的結(jié)果出現(xiàn)了明顯的差異，表明噪聲預測方法能夠適用于當邊界為聲學非緊致邊界時的氣動噪聲預測，考慮了邊界對聲波的散射效應。

圖5 不同計算方法得到的聲壓指向性圖

圖6和圖7分別給出了噪聲預測方法得到的聲壓分布云圖和各成分指向性圖?？梢钥闯?，圓柱繞流誘發(fā)的氣動噪聲總聲壓呈明顯的偶極子特征，總聲壓和邊界輻射聲壓相當，直接入射聲壓對總聲壓的貢獻小，說明低頻時的噪聲主要由邊界輻射噪聲組成，這與Takaishi等[12]的結(jié)果一致。

圖6 總噪聲聲壓分布云圖

通過以上分析可知，邊界存在時，非定常繞流誘發(fā)的氣動噪聲的邊界輻射噪聲從本質(zhì)上看有兩部分：一是通常認為的等效偶極子聲源產(chǎn)生的噪聲，這是邊界與流動相互作用產(chǎn)生的；二是邊界對聲波傳播的散射效應，這部分噪聲來源于聲波傳播與邊界的相互作用。從本質(zhì)上看，偶極子源噪聲和邊界對聲波的散射都是邊界對聲波傳播的散射效應的外在表現(xiàn)，所不同的是，偶極子源噪聲是流場計算所能夠捕捉到的聲學脈動量，噪聲預測方法和采用自由空間格林函數(shù)求解FW-H方程都能預測，而非緊致邊界對聲波傳播的散射是噪聲預測方法中波動方程所能捕捉到的聲學脈動量，本文噪聲預測方法的意義就在于能準確預測這部分噪聲，而采用自由空間格林函數(shù)求解FW-H方程則不能捕捉到這部分聲學脈動量。

圖7 噪聲各成分指向性的對比圖

3 結(jié)語

通過數(shù)值方法尋找滿足非緊致邊界的精確格林函數(shù)，提出了一種基于FW-H方程的非緊致邊界氣動噪聲數(shù)值預測方法，適用于幾何外形復雜的邊界與非定常流動相互作用誘發(fā)的氣動噪聲的數(shù)值預測。通過對圓柱繞流誘發(fā)的氣動噪聲的數(shù)值預測，對該方法進行了驗證，結(jié)果表明：

(1)通過采用圓柱邊界的精確格林函數(shù)的解析解來求解提出的氣動噪聲數(shù)值預測基本方程，與本文數(shù)值計算方法的結(jié)果進行對比，驗證了精確格林函數(shù)數(shù)值計算方法的正確性；

(2)與采用自由空間格林函數(shù)求解FW-H方程的傳統(tǒng)方法相比，提出的氣動噪聲數(shù)值預測方法能夠考慮非緊致邊界對聲波的散射效應，能更為準確的預測流動誘發(fā)的氣動噪聲。

[1]Lighthill M J.On sound generated aerodynamically[M].I. General Theory.Proc.Roy.London Soc.1952,211A, 1107:564-587.

[2]Ffowcs Williams J E,Hawkings D L.Sound generation by turbulence and surface in arbritary motion[J].Proceedings of the Royal Society of London,Series A,Mathematics and Physical Sciences,1969,264:321-342.

[3]Farassat F,Jay H C.Towards an airframe noise prediction methodology:Survey of current approaches[R].AIAA Paper,2006-0210.

[4]Curle N.The influence of solid boundaries upon aerodynamic sound[J].Proceedings of the Royal Society of London,Series A,Mathematics and Physical Sciences, 1955,231,1187:505-514.

[5]Goldstein M E.Aeroacoustics[M].New York:McGraw-Hill International Book Co.,1976.

[6]Hu F Q,Hussaini Y H,Manthey J L.Low-dissipation and low-dispersion Runge-Kutta schemes for computational acoustics[J].Journal of Computational Physics,1996, 124:177-191.

[7]Ffowcs Williams J E,Hall L H.Aerodynamic sound generation by turbulent flow in the vicinity of a scattering half plane[J].Journal of Fluid Mechanics.1970,40(4):657-670.

[8]Ali I,Escobar M,Kaltenbacher M,et al.Time domain computation of flow induced sound[J].Computer&Fluids,2008,37:349-359.

[9]Mao Y,Qi D,Gu Y.Prediction of airfoil noise induced by low Mach number flow[J].Proc.ImechE,Part G,Journal ofAerospace Engineering,2012,226(5):561-573.

[10]Hu F Q,Guo Y P,and Jones A D.On the computation and application of exact greens function in acoustic analogy[C].In 11th CEAS/AIAA Aeroacoustics Conference. Monterey,California,AIAApaper 2005-2986,2005.

[11]Golerfelt X,Perot F,Bailly C,et al.Flow-induced cylinder noise formulated as diffraction problem for low Mach numbers[J].Journal of Sound and Vibration, 2005,287:129-151.

[12]Takaishi T,Miyazawa M,Kato C.Computational method of evaluating noncompact sound based on vortex sound theory[J].The Journal of Acoustical Society of America,2007,121:1353-1361.

[13]宋亞輝，劉秋洪，蔡晉生.氣動噪聲預測中精確格林函數(shù)的數(shù)值計算方法[J].噪聲與振動控制，2013，33(2)：5-10.

[14]Reinaldo M O,Julio R M,Fabio S.Two and three-dimensional simulation of sound generated by flow around a circular cylinder[J].AIAAPaper,2009,3270.

[15]Cantwell B,Coles D.An experimental study of entrainment and transport in the turbulent near wake of circular cylinder[J].Fluid Mech.,1983,136:321-374.

Computational Method ofAerodynamic Noise Induced by Airflow over Non-compact Boundary

SONG Ya-hui1,QIN Hao1,LIU Qiu-hong2
(1.Chinese Flight Test Institute,Xi’an 710089,China; 2.School ofAeronautics,Northwestern Polytechnical University,Xi’an 710072,China)

A computational method was proposed to evaluate the aerodynamic noise generated by unsteady flow over non-compact boundary.Firstly,the boundary element method was applied to compute the exact Green’s Function that satisfies the specific boundary conditions.Then,the numerical results of the exact Green’s Function were used to solve the FWH equation to predict the far field noise.This numerical method can predict the radiation noise from the aerodynamic sources as well as the scattering noise from the non-compact surfaces.Especially,this approach is suitable for the boundaries with large sizes and complex geometric shapes.The present method was applied to predict aerodynamic noise generated by turbulent flows at Reynolds number of 90 000 and Mach’s number of 0.2 over a circular cylinder with non-compact boundary conditions.Results were compared with those of FW-H equations based on the free-space Green’s function to verify the effectiveness of this method.

acoustics;aerodynamic noise;FW-H equation;non-compact boundary;numerical method

O42；O355

：A

10.3969/j.issn.1006-1335.2015.02.011

1006-1355(2015)01-0052-06

2014－09－12

國家自然科學基金(11002116)；西北工業(yè)大學基礎(chǔ)研究基金(GCKY1006)

宋亞輝（1985－），男，河南省鹿邑縣人，碩士，工程師，主要研究方向：飛機結(jié)構(gòu)動力學。E-mail:songyahuilym@163.com

劉秋洪，男，副教授，主要研究方向：計算流體力學和氣動聲學計算。E-mail:qhliu@nwpu.edu.cn