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

申肖雪，盧浩，肖友剛

?

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

申肖雪，盧浩，肖友剛

(中南大學(xué)交通運(yùn)輸工程學(xué)院，湖南長(zhǎng)沙 410075)

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

復(fù)雜聲源；時(shí)域；有限體積法；聲輻射

0 引言

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

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

1 聲波控制方程

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

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

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

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

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

2.1 單極子聲源

2.1.1 單極子聲源

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

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

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

2.1.2 仿真模型

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

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

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

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

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

2.1.3 計(jì)算結(jié)果

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

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

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

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

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

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

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

2.3 復(fù)雜時(shí)變聲源

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

2.3.1 問題描述與解析解

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

徑向速度的表達(dá)式可寫為

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

2.3.2 圓柱體聲源輻射仿真

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

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

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";

}

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

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

(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è)置速度勢(shì)的方式將聲源引入數(shù)值計(jì)算的方法，推導(dǎo)了由給定的聲源表面聲壓或振動(dòng)位移得到速度勢(shì)的公式。該方法在建立網(wǎng)格初期即考慮到聲源的作用，聲源以邊界條件的方式出現(xiàn)，求解器只需要求解Helmholtz方程，不需要將聲源項(xiàng)加入到Helmholtz方程，使得聲場(chǎng)仿真變得簡(jiǎn)單、直觀。

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

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

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

[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] 杜功煥, 朱哲民, 龔秀芬. 聲學(xué)基礎(chǔ)[M]. 3版, 南京: 南京大學(xué)出版社, 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] 曹景記, 唐曉明, 蘇遠(yuǎn)大. 多極聲源在裸眼井外的縱橫波輻射特性[J]. 聲學(xué)技術(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

國(guó)家自然科學(xué)基金資助項(xiàng)目(50975289、51275531)

申肖雪(1993－), 男, 河北景縣人, 碩士研究生，研究方向?yàn)殍F路噪聲預(yù)測(cè)及控制。

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