劉 超，劉秋洪，蔡晉生

（西北工業(yè)大學 翼型葉柵空氣動力學國防科技重點實驗室，西安710072）

用快速多極方法預測圓柱繞流的氣動噪聲

劉 超，劉秋洪，蔡晉生

（西北工業(yè)大學 翼型葉柵空氣動力學國防科技重點實驗室，西安710072）

采用聲模擬理論預測氣動噪聲時需要大量的計算時間，快速多極方法將傳統(tǒng)點對點計算轉(zhuǎn)變?yōu)辄c集之間的相互作用，可以有效加速計算。基于二維自由空間格林函數(shù)的分波展開方式，推導了FW-H方程應用快速多極變換后的積分核函數(shù)與計算公式。計算了低馬赫數(shù)圓柱繞流的非定常流場；并由此預測了氣動聲源。隨后，分別采用傳統(tǒng)方法和快速多極方法計算其聲場分布。結(jié)果表明，基于分波展開方式的快速多極方法能準確計算圓柱繞流氣動噪聲，在頻率較低時能大幅減少聲場計算時間，且觀測點數(shù)越多，加速效果越明顯。

聲學；氣動噪聲；數(shù)值預測；快速多極；聲模擬

快速多極方法[4]對積分核函數(shù)引入多極擴展，將源點和場點分離，用節(jié)點集之間相互作用取代節(jié)點之間一對一計算，使計算量和存儲量大幅減少，目前該方法已廣泛應用于力學、電磁學等諸多領(lǐng)域。在聲學方面，Rokhlin[5]最早使用快速多極方法求解聲散射問題，隨后國內(nèi)外很多學者將其與邊界元結(jié)合用于振動聲學計算[6]。然而，在氣動聲學方面，關(guān)于將快速多極方法應用于聲模擬理論的文獻較少，目前已知Wolf[7—9]利用寬頻快速多極方法與聲模擬理論相結(jié)合，數(shù)值計算了圓柱和翼型等流場的氣動噪聲，但在其文章[7]中僅給出了多極擴展系數(shù)的轉(zhuǎn)移公式，而未給出相應的FW-H方程的變換公式。在國內(nèi)方面，目前為止還沒有快速多極方法應用于氣動聲學計算的公開報道。

本文借鑒Wolf[7]的思路，將快速多極方法應用于聲模擬理論，詳細推導了FW-H方程應用快速多極變換后的積分核函數(shù)與計算公式，并編寫相應的計算程序，數(shù)值預測了二維圓柱繞流的氣動噪聲，與傳統(tǒng)方法對比驗證了快速多極方法的有效性和快速性。

1 氣動聲學計算的快速多極方法

1.1 傳統(tǒng)氣動聲學計算方法

聲模擬理論在忽略聲場對流場影響的假設(shè)下，將近場非線性流動區(qū)域的聲源和遠場線性區(qū)域的聲輻射分離求解，其計算流場中固體物繞流誘發(fā)氣動噪聲的FW-H方程[3]為

其中H(f)是Heaviside函數(shù)，c0是介質(zhì)中聲速，Tij為Lighthill應力張量，G0為自由空間格林函數(shù)，p為壁面壓力，vn為壁面運動的法向速度，Ω是聲源區(qū)域，S是固體邊界。當邊界靜止時，式(1)的頻域表達式為

1.2 快速多極方法應用于聲學控制方程

快速多極方法通過構(gòu)造自適應四叉樹將所有離散的聲源點劃分成多個點集，然后將積分核函數(shù)在樹結(jié)構(gòu)的葉子節(jié)點處多極展開并計算多級擴展系數(shù)，之后通過多極擴展系數(shù)轉(zhuǎn)移公式計算每個葉子節(jié)點處的局部擴展系數(shù)，最后由所有葉子節(jié)點處的局部擴展系數(shù)計算不同觀測點處聲壓，具體流程如下圖1所示：

對于二維問題，式(2)中自由空間格林函數(shù)為

設(shè)聲場點x和聲源點y在極坐標系中的表達式為x=(ρx,θx)和y=(ρy,θy)。根據(jù)Graf加法定理[10]得零階第1類漢克爾函數(shù)的級數(shù)表達式為

圖1 快速多極方法示意圖

頻率ω固定時，將G0(x,y)在聲源點y周圍的結(jié)構(gòu)中心yc多極展開為

將上述二維自由空間格林函數(shù)的多極展開式(5)代入到原頻域聲學控制方程(2)，得到兩個需要重新計算的積分核函數(shù)如下

將式(7)代入控制方程(2)，并令多極擴展系數(shù)

則式(2)變換為

1.2.1 多極擴展系數(shù)轉(zhuǎn)移

當多極展開中心由子正方形中心yc轉(zhuǎn)移到父正方形中心yc'時，多極擴展系數(shù)轉(zhuǎn)移按下式計算

可以看出新的多極擴展系數(shù)Mn(yc')可以從老的多極擴展系數(shù)Mm(yc)計算得到，從而避免了重復計算。通過多極擴展系數(shù)的轉(zhuǎn)移可以把同層的子正方形的多極擴展系數(shù)總和轉(zhuǎn)移到包含更大邊界的父正方形的多極擴展系數(shù)中去。此時的頻域聲學控制方程(9)變換為

1.2.2 局部擴展

當多極擴展系數(shù)向局部擴展系數(shù)轉(zhuǎn)移，把原配置點y的父正方形的多極擴展系數(shù)（中心在yc'）轉(zhuǎn)移到了包含x的父正方形的局部擴展系數(shù)（中心在xc）時，局部擴展系數(shù)按下式計算

其中Ln(xc)表示中心在xc的局部擴展系數(shù)。此時的頻域聲學控制方程(11)變換為

1.2.3 局部擴展系數(shù)轉(zhuǎn)移

當局部擴展中心由父正方形的中心xc轉(zhuǎn)移到子正方形的中心xc'時，局部擴展系數(shù)轉(zhuǎn)移按下式計算

通過局部擴展系數(shù)轉(zhuǎn)移，將父正方形的局部擴展系數(shù)轉(zhuǎn)移到了所包含的子正方形局部擴展系數(shù)。此時的頻域聲學控制方程(13)變換為

在實際計算中，以上公式的無窮級數(shù)必須截斷，由此引入截斷誤差并影響聲場計算精度。截斷項數(shù)越多，截斷誤差越小，聲場計算精度越高，但計算時間變長，計算效率下降；若截斷項數(shù)不足，聲場計算則會出現(xiàn)失真。因此需要合理的選擇截斷項數(shù)以同時兼顧計算精度和計算效率。

1.3 快速多極方法的計算過程

1.3.1 建立分層樹結(jié)構(gòu)

首先建立根節(jié)點使其包含所有流場和聲場網(wǎng)格，之后根據(jù)每個節(jié)點內(nèi)包含流場和聲場網(wǎng)格的數(shù)量遞歸劃分四叉樹，最后建立鄰居節(jié)點、相互作用節(jié)點和遠場節(jié)點的索引。

1.3.2 上行遍歷計算多極擴展系數(shù)

上行遍歷指由樹結(jié)構(gòu)的葉子節(jié)點開始，逐層上行計算所有節(jié)點的多極擴展系數(shù)，直至第二層。對于葉子節(jié)點，直接使用公式(8)計算其多極擴展系數(shù)，此時需要計算貝塞爾函數(shù)的偏導數(shù)及其積分，所需計算時間較長；對于非葉子節(jié)點，直接使用公式(10)將子節(jié)點的多極擴展系數(shù)轉(zhuǎn)移到父節(jié)點的多極擴展系數(shù)。

1.3.3 下行遍歷計算局部擴展系數(shù)

下行遍歷指由樹結(jié)構(gòu)的第二層開始，逐層下行計算出所有節(jié)點的局部擴展系數(shù)，直至葉子節(jié)點。所有節(jié)點的局部擴展系數(shù)來源于相互作用節(jié)點和遠場節(jié)點兩部分貢獻量的和。前者的貢獻通過公式(12)計算，后者的貢獻通過公式(14)計算。

1.3.4 計算總積分

樹結(jié)構(gòu)的多極擴展系數(shù)及局部擴展系數(shù)計算完成后，針對不同的聲場觀測點，計算兩個積分核函數(shù)在固體壁面和流體流動區(qū)域積分的總和。其中，對于遠離聲場觀測點的流場聲源點，使用快速多極算法計算積分；對于臨近觀測點的聲源點則直接計算其積分。

2 算例驗證與分析

將上述方法應用于二維圓柱繞流的氣動聲學計算，使用Fortran 90編寫相應的計算程序。

2.1 流場計算模型與結(jié)果

流場計算選取直徑D=0.019 m的圓柱，均勻來流馬赫數(shù)為Ma=0.2，雷諾數(shù)為Re=9×104，遠場采用無反射邊界條件，時間推進步長取Δt=2×10-5s，使用k-w SST兩方程模型求解湍流雷諾應力，圓柱局部網(wǎng)格如圖2所示。

圖2 圓柱局部網(wǎng)格

圖3給出了圓柱升力系數(shù)Cl和阻力系數(shù)Cd隨時間的變化?？梢钥闯?，二維模型所得到的升力系數(shù)和阻力系數(shù)有規(guī)整的周期性波動，其中升力系數(shù)的脈動幅值遠大于阻力系數(shù)，因此卡門渦街中渦脫落形成升力激發(fā)的噪聲要遠大于壓差阻力及摩擦阻力引起的噪聲。

圖4為升阻力系數(shù)變化的功率譜密度(PSD)，其中升力系數(shù)的基頻也即尾渦脫落頻率為f0= 897.53 Hz，相應的無量綱斯特勞哈數(shù)為St=0.251，與Orselli[11]的計算結(jié)果(St=0.235)較符合。此外升力系數(shù)較明顯的諧波分量為基頻的奇數(shù)倍，而阻力系數(shù)較明顯的諧波分量為基頻的偶數(shù)倍，與Orselli[11]和Takaishi[12]得到的結(jié)論一致。

圖4 升阻力系數(shù)隨時間變化的功率譜密度

2.2 聲場計算結(jié)果與分析

選取36個周期的流場數(shù)據(jù)輸出作為氣動聲源，采樣步長與流場計算的物理時間推進步長一致，總采樣步數(shù)2 048步，頻率分辨率約為24.2 Hz。聲場云圖計算使用3 000個觀測點，布置在以圓柱中心點為圓心，半徑10D到128D的環(huán)形區(qū)域內(nèi)。自適應四叉樹劃分時每個葉子節(jié)點內(nèi)最多包含50個聲場觀測點，得到的自適應樹結(jié)構(gòu)共有5層，149個節(jié)點，128個葉子節(jié)點。

圖5為使用傳統(tǒng)方法和使用快速多極方法得到的渦脫落頻率f=f0下聲場觀測點位于128D處的聲場指向性對比圖，其中多極展開的截斷項數(shù)取p= 20。圖中TM-dip和TM-quad表示傳統(tǒng)方法得到的偶極子和四極子噪聲，F(xiàn)M-dip和FM-quad則表示使用快速多極方法得到的偶極子和四極子噪聲。從圖中可以看出，偶極子噪聲在空間呈正8字形，聲壓級幅值占優(yōu)，而四極子噪聲則呈現(xiàn)光滑規(guī)則的四花瓣狀，聲壓級較小。此外，兩種方法得到聲場結(jié)果一致，驗證了本文快速多極方法的準確性。

圖5 渦脫落1倍頻時聲場指向性對比圖

圖6和圖7給出了升力系數(shù)變化的基頻即渦脫落頻率f=f0下聲場總聲壓云圖和聲壓實部，從中可以看出升力脈動引起聲場的聲壓級在空間呈上下對稱并沿豎直方向傳播，此外，多極展開的截斷項數(shù)取p=20可以準確計算一倍頻下聲場聲壓云圖，下面將討論2倍頻時截斷項數(shù)取值過小導致聲場失真的情況。

圖6 1倍頻截斷項數(shù)p=20時聲場云圖

圖7 一倍頻截斷項數(shù)p=20時聲壓實部

圖8給出了截斷項數(shù)繼續(xù)取p=20時渦脫落2倍頻f=2 f0下聲場總聲壓云圖，可以看出此時部分聲場失真，這是由于頻率提高，聲壓量級減小，多極擴展系數(shù)所需計算精度提高，截斷項數(shù)需要相應的增加才能準確計算聲場分布。

圖8 二倍頻截斷項數(shù)p=20時聲場云圖

圖9和圖10給出了p=35時二倍頻率下聲場總聲壓云圖和聲壓實部，此時聲場呈倒8字，主要由阻力脈動引起，與圖4功率譜密度中計算得到阻力系數(shù)較明顯的諧波分量為基頻的偶數(shù)倍的結(jié)論相符，進一步驗證了本文快速多極方法的準確性。

圖9 二倍頻截斷項數(shù)p=35時聲場云圖

圖10 二倍頻截斷項數(shù)p=35時聲場實部

圖11給出了渦脫落三倍頻f=3 f0時聲場總聲壓云圖，此時聲場恢復正8字形狀，對應升力系數(shù)脈動的2次諧波分量，但聲壓量級遠遠小于一倍頻，使用快速多極方法取截斷項數(shù)p=50計算得到聲場結(jié)果仍有輕微的失真。

圖11 三倍頻截斷項數(shù)p=50時聲場云圖

表1列出了在圓柱繞流渦脫落的前3個整數(shù)倍頻下使用傳統(tǒng)方法和使用快速多極方法求解FW-H方程所需要的計算時間對比。從表中可以看出，使用傳統(tǒng)方法計算3個頻率下的聲場所需時間差別不大，而使用快速多極方法在頻率較低時能大幅縮減聲場計算時間，有很高的實用價值和工程應用前景。

表1 快速多極方法加速效率

在1.3.2節(jié)提到，對于葉子節(jié)點，計算其包含的所有聲源點在該葉子節(jié)點處的多極擴展系數(shù)時，需要計算貝塞爾函數(shù)的偏導數(shù)及其積分，所需計算時間較長，而由葉子節(jié)點的局部擴展系數(shù)計算觀測點處聲壓時需要的時間較短，所以觀測點數(shù)量增多對整體計算時間影響較小。下表列出了渦脫落一倍頻時，觀測點數(shù)量變化導致快速多極方法加速計算效率變化的比較，從中可以看出，基于分波展開方式的快速多極方法特別適合于廣域空間內(nèi)大量聲場點的聲壓計算。

表2 聲場觀測點數(shù)量對加速效率的影響

3 結(jié)語

(1)快速多極方法通過格林函數(shù)的多極展開將聲源點和觀測點分離，通過構(gòu)造自適應四叉樹使傳統(tǒng)FW-H方程中格林函數(shù)點對點積分轉(zhuǎn)化為點集之間的相互作用，從而使聲場計算時間大幅縮短；(2)本文使用的分波展開方式由流場信息計算葉子節(jié)點的多極擴展系數(shù)時所需時間較長，由葉子節(jié)點的局部擴展系數(shù)計算觀測點處聲壓時需要的時間較短，因此聲場觀測點數(shù)越多，加速效果越明顯；

(3)隨著頻率提高，聲壓量級降低，需要增大多極截斷項數(shù)來準確計算聲場分布，導致高頻時加速效率下降，但工程中往往更關(guān)注聲壓級較高的低頻噪聲問題，因此本文快速多極加速計算方法具有實際應用價值。

[1]Lighthill M J.On sound generated aeroynam ically.I.general theory[J].Proceedings of the Royal Society of London,SeriesA,1952,A211:564-587.

[2]Lighthill M J.On sound generated aerodynam ically.II.Turbulence as a source of sound[J].Proceedings of the Royal Society of London,Series A,1954,222:1-32.

[3]Curle N.The influence of solid boundaries on aerodynam ic sound[J].Proceedings of the Royal Society of London,SeriesA,1952,213,1187:505-514.

[4]Greengard L,Rokhlin V.A fast algorithm for particle simulations[J].Journal of Computational Physics.1987, 73:325-348.

[5]Rokhlin V.Rapid solution of integral equations of scattering theory in two dimensions[J].Journal of Computational Physics,1990,86:414-439.

[6]崔曉兵，季振林，武 耀.快速多極子邊界元法預測船舶艙室噪聲[J].噪聲與振動控制，2012，32(6)：179-183.

[7]Wolf W R,Lele S K.Acoustic analogy formulations accelerated by fast multipole method for two-dimensional aeroacoustic problems[J].AIAA Journal,2010,48(10): 2274-2285.

[8]Wolf W R,Lele S K.Aeroacoustic integrals accelerated by fast multipole method[J].AIAA Journal,2011,49(7): 1466-1477.

[9]Wolf W R,Lele S K,Jothiprasad G,Cheung L.Investigation of noise generated by a DU 96 airfoil[C].18 th AIAA/CEAS Aeroacoustics Conference,Colorado Springs,CO,2012,2012-2055.

[10]Abramow ita M,Stegun I A.Handbook of mathematical function[M].New York:Dover,1965.

[11]Orselli R M,Meneghini J R,Saltara F.Two and threedimensional simulation of sound generated by flow around a circular cylinder[J].AIAA paper,2009,3270.

[12]Takaishi T,M iyazawa M,Kato C.Computational method of evaluating noncompact sound based on vortex sound theory[J].Journal of Acoustical Society of America, 1998,233-253.

Fast Multipole Method Applied in Prediction ofAeroacoustics Induced by a Circular Cylinder

LIU Chao,LIU Qiu-hong,CAI Jin-sheng

(National Key Laboratory of Aerodynam ic Design and Research, Northwestern Polytechnical University,Xi’an 710072,China)

Prediction of aeroacoustics using acoustic analogy is a time-consum ing process,while the fast multipole method which changes traditional way of node-to-node computing into set-to-set interaction can accelerate the process effectively.In this paper,the FW-H equation and its integral kernel function are derived w ith the fast multipole method based on the partial-wave expansion formulation of free-space Green’s function.The unsteady flow field w ith low Mach’s number near a two-dimensional circular cylinder is computed and exported as the sound source,and the sound field is obtained via traditional method and the fast multipole method.The results show that the fast multipole method based on partial-wave expansion can calculate the aerodynam ic noise accurately,and reduce the computing time greatly for relatively low frequencies.And the acceleration effect is more obvious w ith larger number of observers.

acoustics;aerodynamic noise;numerical prediction;fast multipole method;acoustic analogy

1006-1355(2014)04-0123-05+133

O42 < class="emphasis_bold">文獻標識碼：A DOI編碼：

10.3969/j.issn.1006-1335.2014.04.027

隨著國內(nèi)外商用飛機的迅速發(fā)展，飛機噪聲預測與控制受到人們的普遍關(guān)注，基于計算流體力學的計算氣動聲學逐漸成為研究熱點。然而，傳統(tǒng)的氣動聲學混合計算方法[1—3]需要計算所有流場聲源點和遠場觀測點之間格林函數(shù)的偏導數(shù)及其積分，這個過程需要大量的計算時間，而應用于飛機噪聲預測等大規(guī)模聲學計算問題的計算耗時更是難以接受，因此發(fā)展快速和有效的氣動噪聲計算方法越來越受到人們的重視，同時也具有十分重要的工程應用價值。

2013-10-30

國家自然科學基金(基金編號：11002116)

劉超（1990-），男，河北滄州人，碩士研究生，主要研究方向為計算流體力學和氣動聲學計算。

E-mail:lc1990@mail.nwpu.edu.cn