流域網(wǎng)格中嵌入式復雜聲源的時域傳播仿真

申肖雪，盧浩，肖友剛

?

流域網(wǎng)格中嵌入式復雜聲源的時域傳播仿真

申肖雪，盧浩，肖友剛

(中南大學交通運輸工程學院，湖南長沙 410075)

針對聲學商業(yè)軟件較難模擬任意形狀、時變復雜聲源的聲輻射問題，使用有限體積法在時域內(nèi)求解無聲源項亥姆霍茲(Helmholtz)方程，將復雜聲源嵌入到有限體積單元節(jié)點，推導了由給定聲源表面聲壓或振動位移得到速度勢公式，提高了聲源處理的靈活性和計算效率。該方法允許對初始場問題及復雜時變聲源聲輻射進行仿真。對常見聲源及二階圓柱體聲源聲輻射進行了數(shù)值模擬，結(jié)果與解析解及商業(yè)軟件結(jié)果進行了對比，誤差均小于15%。程序具有良好的封裝性及通用性，可以靈活地對不同聲源進行組合，得出任意復雜聲源時域的傳播規(guī)律，為復雜聲源聲輻射等線性聲學問題的研究提供了一個可靠的平臺。

復雜聲源；時域；有限體積法；聲輻射

0 引言

要研究復雜聲源在復雜密閉空間內(nèi)的輻射聲波特性，如何容易、方便地設置聲源便成為必須要面對的問題。目前聲學商業(yè)軟件中聲源類型大多為單/偶/四極子聲源、平面波聲源、混響聲源等，而對復雜形狀、任意變化的時變聲源則較難進行仿真，在研究論文中也比較少見。對于聲源的仿真及聲場計算，MILLAN[1]使用了有限元素法(Finite Element Method，F(xiàn)EM)和邊界元法(Boundary Element Method，BEM)對時域內(nèi)聲波傳播進行仿真，對單極子聲源的處理是在給定球形網(wǎng)格并在球形表面設定位移。Takuya Oshima[2]開發(fā)了一款時域內(nèi)聲波傳輸?shù)耐ㄓ们蠼馄鱬otentialWaveFoam。Schmalz J[3]在OpenFOAM[4-5]下實現(xiàn)了Lighthill和Curle聲類比，聲源作用分別體現(xiàn)在Lighthill方程及FW-H方程。Nilsson[6]在OpenFOAM和CALFEM中實現(xiàn)聲類比，對聲源的處理是在求解器每一個時間步計算完成后根據(jù)公式計算聲源項。Ji H[7]使用有限體積法計算聲傳播，對聲源的處理是將聲源項加入Helmholtz方程，然后對整個方程進行離散求解。

以上文獻中無聲源項的Helmholtz方程適用于初始場問題的求解，而對于時變聲源，則是將聲源項加入到聲波方程，然后對其進行離散、矩陣求解等處理，大大占用了計算資源。另外，大部分聲學商業(yè)軟件可以對簡易聲源進行處理，但對復雜形狀或復雜時變聲源的處理技術(shù)較為缺乏。本文將聲源以邊界條件的形式巧妙地嵌入到數(shù)值計算的網(wǎng)格中，將求解時域內(nèi)含聲源項的Helmholtz方程的求解問題轉(zhuǎn)化為求解無聲源項Helmholtz方程，將聲源以邊界條件的形式加入到計算過程。在計算初期對流域劃分網(wǎng)格時即考慮聲源作用及位置，然后根據(jù)聲源類型的不同在現(xiàn)有計算域網(wǎng)格中“掏去”不同形狀的網(wǎng)格，使之形成邊界，然后再在此邊界上設定速度勢，以模擬聲源。在提高對復雜聲源處理的靈活性的同時，提高了計算效率。

1 聲波控制方程

理想氣體的Helmholtz方程[8]為

將式(2)～(4)代入式(5)，可得

由于本文側(cè)重于聲源的實現(xiàn)方式，不對其方程的離散方法及修正算法進行細述。采用有限體積法計算聲波的傳播，對方程(1)進行空間離散，對式(6)進行時間離散，離散格式及修正算法參考文獻[2]，算法迭代過程、誤差分析見文獻[9]。

2 嵌入式網(wǎng)格表面設置速度勢作聲源

空氣噪聲聲源可以看作是單極子、偶極子、四極子聲源的疊加。此節(jié)對時域內(nèi)單極子、偶極子和四極子聲源的聲輻射進行計算，將結(jié)果進行FFT變換，并與商業(yè)軟件的頻域計算結(jié)果進行對比。然后對復雜時變二階圓柱聲源聲輻射進行了時域仿真，結(jié)果與理論值吻合。

2.1 單極子聲源

2.1.1 單極子聲源

如圖1所示，在邊長為1 m的立方體聲腔中心存在一個直徑為0.01 m的脈動球形聲源。脈動球源是進行均勻漲縮振動的球面聲源，也就是在球源表面上各點沿著徑向做同振幅、同相位的振動。

圖1 在立方體聲腔中心的半徑為0.01m的球形聲源

其中：為聲源的強度，定義=0.2 N/m，=0.01 m，則球形聲源的表面上的聲壓為20 Pa。

2.1.2 仿真模型

在立方體中心有一球形聲源，由于整個模型關于立方體中心點對稱，所以取立方體的1/4和球形聲源的1/4進行仿真即可，如圖2所示。原立方體的3個表面都設置為吸聲邊界，吸聲系數(shù)為0.01，根據(jù)吸聲系數(shù)的定義有[9]

圖2 計算區(qū)域網(wǎng)格，原尺寸的1/4

Fig.2 The meshes in computation region, 1/4 of origin geometry

使用OpenFOAM自帶的groovyBC邊界條件設置工具，設置隨時間變化的邊界條件，在聲腔左側(cè)采用Dirichlet邊界條件，根據(jù)式(2)的離散形式：

其中：為聲波振幅，為聲波頻率。這樣就能保證在球形表面上的聲壓以正弦方式變化。

2.1.3 計算結(jié)果

在初始時刻，聲場中各點的聲壓值為0，聲源聲壓在初始時刻開始以正弦方式變化，模擬時間為1.5 s。P點聲壓變化如圖3所示。從圖3可以看出，在1.0 s時間后，聲壓趨于穩(wěn)定。

計算時間步長為10-6s，每20個時間步長取一次數(shù)據(jù)，總共有7 500個數(shù)據(jù)。取其中最后的2 048個時域數(shù)據(jù)進行FFT變換，得到其頻域譜如圖4所示。從圖4可以看到，該點的聲壓在200 Hz時為0.17 Pa，商業(yè)軟件中有限元法得到的該點聲壓為0.18 Pa，誤差為5.6%。考慮到信號的采樣頻率及采樣長度，頻譜中幅值最大值所在的頻率在200 Hz左右存在偏移。由于截斷信號并非整周期，頻譜在200 Hz以外發(fā)生頻率泄漏。根據(jù)單極子聲源的聲壓變化規(guī)律，設置邊界上速度勢，模擬單極子聲源的聲輻射。

圖3 1.5 s內(nèi)P點的聲壓變化

圖4 時域數(shù)據(jù)的FFT變換

2.2 偶極子聲源、四極子聲源

本文對偶極子和四極子聲源的仿真同樣采用在流域聲腔網(wǎng)格內(nèi)嵌入球形網(wǎng)格的方式進行處理。在邊長為1 m的立方體聲腔中心，偶極子聲源由兩個聲壓幅值相同、相位相反的關于中心對稱的脈動球源組成，兩個聲源的球心位置為(0.5，0.5，0.487 5)和(0.5，0.5，0.512 5)，半徑為0.006 25 m，球形表面設置聲壓幅值均為20 Pa，頻率為200 Hz，但是相位相反。聲腔表面設為無反射邊界，在時域內(nèi)對聲場進行計算。在聲壓穩(wěn)定后，兩聲源的中垂面上由于兩聲源幅值相等而相位相反，聲壓相互抵消為0。在兩球源中心連線上的聲壓疊加增強，聲壓最大。四極子聲源由一對極性相反的偶極子組成。同單極子聲源一樣，使用snappyHexMesh命令指定四個球源的球心為(0.5，0.5，0.487 5)、(0.5，0.5，0.512 5)和(0.5，0.487 5，0.5)、(0.5，0.512 5，0.5)，半徑為0.006 25 m及相關劃分網(wǎng)格參數(shù)，劃分好后施加邊界條件即可將聲源以邊界條件的方式容易地加入到數(shù)值計算。同之前單極子聲源的處理相似，將偶極子聲源和四極子聲源輻射聲壓時域內(nèi)的結(jié)果經(jīng)FFT變換后，同解析解及商業(yè)軟件頻域結(jié)果進行對比，結(jié)果吻合。使用嵌入式聲源方法大大簡化了計算過程。

對比MATLAB中k-wave工具箱中對偶極子、四極子聲源的處理，將聲源項加入二階Helmholtz方程[10]，得

2.3 復雜時變聲源

對簡易聲源的聲輻射理論研究有不少[11]，但對復雜時變聲源的數(shù)值計算研究并不多。在數(shù)值計算中，大多數(shù)將聲源項加入聲波方程然后將整個方程進行格式離散。對于單、偶極子等簡單聲源可以將其在網(wǎng)格點上設置為點源，將數(shù)據(jù)導入到計算過程。但是對于復雜聲源(尤其近場輻射問題)，不能在方程及求解中設定幾個網(wǎng)格點以代表整個聲源，其表面每個點的變化規(guī)律也可能不同，因此需要盡量細致地描述聲源。另外，對于復雜聲源的聲輻射仿真大多使用邊界元、無限元或格林函數(shù)方法，要求仿真者對聲輻射理論公式等有一定的了解，相關的數(shù)值計算方法不多。本文將任意形式的聲源輻射仿真問題簡化為求解無聲源項的Helmholtz方程問題，將求解過程中的聲源處理問題轉(zhuǎn)化到前處理過程，使用者只需關注前處理過程，這樣使得求解器具有良好的通用性。二階圓柱體聲源變化規(guī)律較單極子、偶極子等聲源復雜，本節(jié)推導了極坐標系下圓柱體表面上速度勢與位移關系式，對給定位移變化的二階振動圓柱體聲輻射案例進行仿真，并對理論值和數(shù)值解進行了對比。

2.3.1 問題描述與解析解

本小節(jié)對圓柱體振動輻射聲場進行了仿真，對數(shù)值計算結(jié)果和理論值進行了對比。假設有一無限長圓柱體[9]，在柱狀坐標系下，半徑隨角度的變化規(guī)律為

徑向速度的表達式可寫為

其中：為振動圓柱體的平均半徑，為圓柱體徑向振動幅值。根據(jù)文獻[12]和文獻[13]，階圓柱體聲源輻射聲場的表達式為

2.3.2 圓柱體聲源輻射仿真

已知圓柱體表面的位移變化，推導聲源表面速度勢變化。當圓柱體振動時，圓柱體表面速度與空氣粒子速度相同。根據(jù)式(3)，圓柱體面上各點徑向振動速度等于速度勢的負梯度，在極坐標中可以表示為

在OpenFOAM中使用blockMesh工具對邊長為1 m的立方體劃分為20*20*20=8 000個六面體單元。然后使用cellSet和refineMesh工具對圓柱體所在的中心區(qū)域進行4層網(wǎng)格細化，最后使用snappyHexMesh工具劃出半徑為0.01 m、高度為1 m的圓柱體網(wǎng)格，并將其移去從而形成圓柱體邊界。劃分好后，將圓柱面的速度勢使用groovyBC設置成為時變量，程序代碼對應為

outCylinder_region0

{ type groovyBC;

valueExpression "-(cos(2*pi*200*time())-1)*2*pi*200*

(pow((pos().x-0.5),2)-pow((pos().y-0.5),2))";

//后半段代表cos(2)

value uniform 0;

gradientExpression "0";

}

其他邊界設置為無反射邊界，以模擬圓柱體聲源在無限空間的聲輻射。

2.3.3 數(shù)值計算結(jié)果與解析值對比

(a)= 0.01 m

(b)= 0.32 m

圖5 距離圓柱體軸心= 0.01 m和= 0.32 m處的聲壓變化

Fig.5 The sound pressures at distance= 0.01 m and= 0.32 m in time domain

由圖5可見，在= 0.01 m處，聲壓幅值為186 Pa。在= 0.32 m的地方，聲壓幅值為0.21 Pa，理論值為0.18 Pa，相差14.3%，吻合度較高。符合聲壓與距離平方成反比的規(guī)律。

3 結(jié)論

(1)本文提出了通過在邊界上設置速度勢的方式將聲源引入數(shù)值計算的方法，推導了由給定的聲源表面聲壓或振動位移得到速度勢的公式。該方法在建立網(wǎng)格初期即考慮到聲源的作用，聲源以邊界條件的方式出現(xiàn)，求解器只需要求解Helmholtz方程，不需要將聲源項加入到Helmholtz方程，使得聲場仿真變得簡單、直觀。

(2)在時域內(nèi)對單/偶/四極子聲源等簡易聲源及復雜時變二階圓柱體聲源的聲輻射進行了數(shù)值計算，單極子聲源聲輻射結(jié)果與有限元結(jié)果相差5.6%，圓柱體聲源的聲輻射結(jié)果與解析解相差14.3%，誤差較小。為更復雜聲源的聲輻射及聲場仿真提供了思路。

(3)該方法將前處理、求解、后處理三個過程全部分開，程序具有良好的封裝性及通用性，可用于研究任意形狀、任意變化的物體的聲輻射及傳播等線性聲學問題。

(4) 該方法對于遠場仿真具有較好效果，而對近場仿真尚有一定的不準確性。方程(7)和方程(15)是自由聲場下的聲源聲輻射公式。假設空間中存在反射波，即非自由聲場下，以算例1為例，=0.01 m的球形表面聲壓為20 Pa，而商業(yè)軟件得到的=0.01 m的結(jié)果并非20 Pa，而是有一定的誤差。說明存在壁面反射時，近場聲場仿真并不精確。

[1] MILLAN BEV. Acoustic time-domain simulation with BEM and FEM[D]. Stuttgart: University of Stuttgart, 2011: 36-39.

[2] OSHIMA T, IMANO M. A full finite-volume time-domain approach towards general-purpose code development for sound propagation prediction with unstructured mesh[C]//Proceedings of Inter-Noise 2008, Shanghai, 2008, 1511-1525.

[3] SCHMALZ J, KOWALCZYK W. Implementation of acoustic analogies in openFOAM for computation of sound fields[J]. Open Journal of Acoustics, 2015, 5(2): 29-44.

[4] JASAK H. OpenFOAM: open source CFD in research and industry[J]. International Journal of Naval Architecture and Ocean Engineering, 2015, 1(2): 89-94.

[5] WELLER G, TABOR G, JASAK H. A tensorial approach to computational continuum mechanics using object-oriented techniques[J]. Computers in Physics, 1998, 12(6): 620-631.

[6] NILSSON, J. Implementation of acoustical analogies in openFOAM and CALFEM[D]. Lund: Lund University, 2010: 13-14.

[7] JI H, XU X, GUO X. Direct FVM simulation for sound propagation in an ideal wedge[J]. Shock & Vibration, 2016, 2016(9): 1-9.

[8] 杜功煥, 朱哲民, 龔秀芬. 聲學基礎[M]. 3版, 南京: 南京大學出版社, 2012. DU Gonghuan, ZHU Zhemin, GONG Xiufen. Fundamentals of Acoustics[M]. 3rdEdition, Nanjing: Nanjing University Press, 2012, 147-152.

[9] MECHEL P. Formulas of Acoustics[M]. Springer Berlin Heidelberg, 2008, 147-152.

[10] TREEBY B E, JAROS J, ROHRBACH D. Modelling Elastic Wave Propagation Using the K-Wave MATLAB Toolbox[C]// 2014 IEEE International Ultrasonics Symposium,Chicago, 2014, 146-149.

[11] 曹景記, 唐曉明, 蘇遠大. 多極聲源在裸眼井外的縱橫波輻射特性[J]. 聲學技術(shù), 2016, 35(6): 487-492. CAO Jingji, TANG Xiaoming, SU Yuanda. Radiation characteristics of a multipole acoustic source in an open borehole[J]. Technical Acoustics, 2016, 35(6): 487-492.

[12] RUSSELL A. On the sound field radiated by a tuning fork[J]. American Journal of Physics, 2000, 68(12): 1139-1145.

[13] JACOBSEN F, JUHL P. Radiation of Sound[R]. Lyngby: Technical University of Denmark, 2011, 25-27.

Time domain propagation characteristics of complex embedded acoustic sources in fluid meshes

SHEN Xiao-xue, LU Hao, XIAO You-gang

(School of Traffic and Transportation Engineering, Central South University, Changsha 410075, Hu’nan, China)

Commercial acousticsoftware is generally hard to simulate sound radiation of complex acoustic sources, which are of arbitrary shape or time-varying. In order to solve this problem, the Helmholtz equation (without acoustic source term) is solved in time domain via Finite Volume Method with embedding acoustic sources to Finite-volume nodes, and the relationship between velocity potential and the pressure/displacement of acoustic source is acquired. This method makes it more convenient and efficient to deal with the acoustic source in numerical calculation. In addition, it allows simulating initial value problem and time-varying source problem. The sound radiation of common acoustic sources and the second order cylindrical acoustic source has been simulated, and the results agree well with analytic solutions and results of commercial code, with numerical errors smaller than 15%. This code has good encapsulation and versatility, via which different types of acoustic sources can be combined to form an arbitrary complex acoustic source and the propagation characteristics can be obtained. It is a good platform of carrying out research on linear acoustics like sound radiation of complex acoustic sources.

complex acoustic source; time domain; Finite Volume Method(FVM); sound radiation

TB533

A

1000-3630(2017)-05-0405-05

10.16300/j.cnki.1000-3630.2017.05.002

2017-01-18;

2017-02-28

國家自然科學基金資助項目(50975289、51275531)

申肖雪(1993－), 男, 河北景縣人, 碩士研究生，研究方向為鐵路噪聲預測及控制。

肖友剛, E-mail: csuxyg@163.com