欽 宇，周楓林,，王煒佳，張玉良，袁小涵

（1.湖南工業(yè)大學(xué) 機(jī)械工程學(xué)院，湖南 株洲 412007；2.衢州學(xué)院 機(jī)械工程學(xué)院，浙江 衢州 324000）

1 背景知識(shí)

離心風(fēng)機(jī)出口氣動(dòng)噪聲、機(jī)械設(shè)備運(yùn)轉(zhuǎn)時(shí)摩擦噪聲等聲源的主要頻率成分在頻譜上具有一定寬度，而多頻聲源的輻射問(wèn)題一直是工程師關(guān)注的重要問(wèn)題。適用于聲輻射分析的工程數(shù)值方法以邊界元法（boundary element method，BEM）[1-3]、有限元法（finite element analysis，F(xiàn)EM）[4]和矩量法（method of moments，MoM）[5]為主要代表，其中BEM 和MoM 都以邊界積分方程為理論基礎(chǔ)，特別適用于無(wú)限域的聲輻射分析。

邊界元法是一種基于邊界積分方程為數(shù)學(xué)基礎(chǔ)的數(shù)值分析方法，其邊界積分方程能夠自動(dòng)滿足無(wú)窮遠(yuǎn)場(chǎng)的輻射條件，并且僅需要對(duì)邊界進(jìn)行離散。較之其他的數(shù)值計(jì)算方法，邊界元法引入了基本解，使其具有了解析和數(shù)值相結(jié)合的特點(diǎn)，從而使其計(jì)算精度相對(duì)較高，是一種比較適用于分析聲學(xué)問(wèn)題的計(jì)算方法[6]，目前被廣泛應(yīng)用于求解裂紋問(wèn)題[7]、擴(kuò)散波問(wèn)題[8]、熱傳導(dǎo)問(wèn)題[9]等。

目前對(duì)聲學(xué)問(wèn)題的分析主要采用頻域法[10]和時(shí)域法[11]，時(shí)域類方法具有直接且通用性好的優(yōu)點(diǎn)，然而在選擇時(shí)間步長(zhǎng)時(shí)會(huì)遇到數(shù)值不穩(wěn)定問(wèn)題。頻域法通過(guò)Helmholtz 控制方程和基本解，推導(dǎo)出邊界積分方程，再將邊界積分方程進(jìn)行線性離散，建立矩陣，求解線性方程組來(lái)獲得相應(yīng)的解[12]。頻域法在頻率成分預(yù)判較準(zhǔn)確的基礎(chǔ)上能夠獲得較高的計(jì)算精度，通過(guò)反變換獲得時(shí)域解。

目前，數(shù)值反變換的方法主要有 Laplace 數(shù)值反變換法[13-14]，Stehfest 算法[15]，離散傅里葉反變換法[16-18]，工程中常利用 Laplace 數(shù)值變換結(jié)合邊界元法計(jì)算結(jié)構(gòu)動(dòng)力響應(yīng)問(wèn)題，首先需要得到變換域中的一組所對(duì)應(yīng)的位移和應(yīng)力，然后通過(guò) Laplace 數(shù)值反變換獲得時(shí)域中相應(yīng)量的解[19-21]。然而對(duì)于多頻噪聲的分析，需要使用多個(gè)采樣頻率并借助數(shù)值反變換技術(shù)才能獲得滿足要求的時(shí)域解答。

本研究中先應(yīng)用頻域邊界元法分析二維多頻聲場(chǎng)輻射問(wèn)題，將時(shí)域描述的控制方程通過(guò) Fourier變換轉(zhuǎn)換為頻域描述的控制方程（Helmholtz 方程），再求解等間隔采樣頻率下的 Helmholtz 方程， 然后采用離散傅里葉反變換（inverse discrete Fouriertransform，IDFT）方法獲得時(shí)域解。最后以兩個(gè)不同區(qū)域的內(nèi)聲場(chǎng)問(wèn)題算例，通過(guò)與參考聲壓進(jìn)行比較，論證了所提方法的正確性。

2 二維聲學(xué)頻域邊界元法

在均勻理想介質(zhì)E中，二維聲學(xué)問(wèn)題在空間內(nèi)傳播的波動(dòng)方程表示形式為

x為配置點(diǎn)空間坐標(biāo)；

c為介質(zhì)中聲音傳播的速度。

對(duì)式（1）進(jìn)行Fourier 變換，假設(shè)時(shí)間項(xiàng)選取為e-jωt，可將瞬態(tài)問(wèn)題從時(shí)間域轉(zhuǎn)化到頻率域進(jìn)行簡(jiǎn)單穩(wěn)態(tài)分析，則頻域聲壓可以表示為

j 為虛數(shù)單位，j2=-1；ω為聲源簡(jiǎn)諧振動(dòng)角頻率。

將式（2）代入（1），消去時(shí)間項(xiàng)，即可得如下Helmholtz 控制方程。

式中k為波數(shù)，且k=ω/c。

控制方程的邊界條件一般滿足以下3 類。

Dirichlet 邊界條件，即聲壓已知：

式中符號(hào)ˉ表示已知值。

Neumann 邊界條件，即法向速度已知：

式中：q(x)為聲壓在x處的法向?qū)?shù)；n為單位外法向，指向背離聲域；ρ為聲學(xué)介質(zhì)的質(zhì)量密度；vn為法向速度。

Robin 邊界條件，即阻抗邊界條件已知：

式中Z為聲阻抗。

將方程（3）轉(zhuǎn)化為等效積分形式，并將試函數(shù)取為二維聲學(xué)問(wèn)題的基本解：

式中：c(x)為自由項(xiàng)；

y為場(chǎng)點(diǎn)空間坐標(biāo)；

q(y)為聲壓在y處的法向?qū)?shù)；

Sy為結(jié)構(gòu)邊界S的子邊界；

Gk(x,y)和為二維和三維聲學(xué)問(wèn)題的基本解，且

r為源點(diǎn)與場(chǎng)點(diǎn)之間的距離，且r=|x-y|；

x為配置點(diǎn)（即計(jì)算點(diǎn)或源點(diǎn)）；

y為場(chǎng)點(diǎn)。

為了得到式（7）的數(shù)值解，需要根據(jù)不同的要求與情況，將邊界劃分為一些小單元。本文主要采用線性離散方法，將邊界S劃分為Ne個(gè)線性單元，則邊界積分方程可以離散為

式中：u為任一邊界的節(jié)點(diǎn)，u=1, 2, …,Ne；

Sv為邊界上的單元v；

qv為在單元v上的法向?qū)?shù)；

將源點(diǎn)置于所有場(chǎng)點(diǎn)上可得到：

式中Huv和Guv均為影響系數(shù)矩陣，為遍及所有節(jié)點(diǎn)單元上的積分總和。

也可將式（1）寫(xiě)為

可以將上述方程組表示為如下矩陣形式：

此時(shí)，利用邊界條件，將已知量移到方程右邊，未知量移至左邊可得：

式中：A為未知聲壓值和位置法向?qū)?shù)組成的Ne×Ne系數(shù)矩陣；

b為已知量組成的Ne階已知向量；

h為由包含未知聲壓值和未知法向?qū)?shù)組成的Ne階向量。

3 離散傅里葉逆變換

為實(shí)現(xiàn)離散傅里葉反變換運(yùn)算，需將連續(xù)信號(hào)進(jìn)行截?cái)唷⒉蓸印?/p>

在頻域邊界元法中，首先取N個(gè)等間隔不同的采樣頻率，分別求解線性方程組（14），得到一組解x(ω)，再通過(guò)離散傅里葉逆變換（IDFT）得到時(shí)域解x(t)。設(shè)x(t)是周期為T(mén)的離散函數(shù)，每個(gè)周期內(nèi)有N個(gè)等間隔頻率離散點(diǎn)。其中Δω=ω0=2π/T。為實(shí)現(xiàn)IDFT 運(yùn)算，必須通過(guò)頻率采樣使頻域函數(shù)為有限離散值，此處采用δ函數(shù)：

式中：Δt為時(shí)間步長(zhǎng)，其中Δt=T/N；

ω0為初始頻率點(diǎn)。

為避免產(chǎn)生頻率混疊現(xiàn)象，信號(hào)采樣頻率

同時(shí)為避免頻混現(xiàn)象，可通過(guò)選擇較小的采樣間隔T（即較高的采樣頻率）來(lái)減小這種誤差。

為保證離散后的信號(hào)能唯一確定原連續(xù)信號(hào)，采樣一般等間隔取值，其頻率采樣信號(hào)為

其Fourier 變換：

根據(jù)δ函數(shù)的篩選性，化簡(jiǎn)求得其Fourier 系數(shù)Ck為

將式（20）代入式（19），可求得：

一個(gè)周期內(nèi)N個(gè)采樣點(diǎn)的復(fù)數(shù)值為

式中N為采樣頻率點(diǎn)的數(shù)量。

由T=2π/ω0=NΔt可得：

傅里葉級(jí)數(shù)可以寫(xiě)成復(fù)指數(shù)形式，根據(jù)歐拉公式

可得到：

4 數(shù)值算例

為驗(yàn)證二維頻域邊界元法結(jié)合離散傅里葉反變換求解多頻聲輻射問(wèn)題的準(zhǔn)確性，設(shè)計(jì)了兩個(gè)驗(yàn)證算例，將計(jì)算結(jié)果與解析解進(jìn)行對(duì)比以驗(yàn)證該方法的正確性。

4.1 旋轉(zhuǎn)L 形邊界區(qū)域內(nèi)部聲壓變化

在一個(gè)旋轉(zhuǎn)L 形邊界區(qū)域條件下進(jìn)行內(nèi)點(diǎn)聲壓變化情況分析，其邊界節(jié)點(diǎn)之間間距設(shè)為0.1 cm，x、y分別為內(nèi)部點(diǎn)的橫坐標(biāo)與縱坐標(biāo)（算例中都將如此使用），計(jì)算求解過(guò)程中一共采用了64 個(gè)邊界節(jié)點(diǎn)，旋轉(zhuǎn)L 形邊界區(qū)域節(jié)點(diǎn)分布情況如圖1所示。

圖1 旋轉(zhuǎn)L 形邊界節(jié)點(diǎn)分布圖Fig.1 Distribution diagram of rotating L-shaped boundary nodes

圖1所示邊界上的聲壓分布按照如下公式給出：

在y=0 的邊界上，

在x=0.8 的邊界上，

在y=0.8 的邊界上，

在x=1.6 的邊界上，

在y=1.6 的邊界上，

在x=0 的邊界上，

域內(nèi)聲壓精確解為

計(jì)算采用在頻率區(qū)間1~4 內(nèi)等間隔取7 個(gè)頻率點(diǎn)，計(jì)算時(shí)間區(qū)間為0.5~3.0 s，在域內(nèi)選取兩個(gè)觀察點(diǎn)，分別為Q1(0.4, 1.2)與Q2(0.5, 1.2)，計(jì)算兩個(gè)內(nèi)部觀察點(diǎn)聲壓值隨時(shí)間變化的過(guò)程，Q1、Q2聲壓值隨時(shí)間變化的情況，如圖2所示。

圖2 Q1、Q2 聲壓值隨時(shí)間變化情況Fig.2 Variation of Q1 and Q2 sound pressure value with time change

從圖2中可以看到，在該方法下兩個(gè)取樣點(diǎn)上的聲壓計(jì)算值與精確值吻合較好。為更準(zhǔn)確直觀地展示計(jì)算結(jié)果和精確解的差異，于表1~2 中列出了隨機(jī)選取的部分時(shí)間點(diǎn)上，對(duì)于Q1和Q2兩個(gè)內(nèi)部觀察點(diǎn)的聲壓以及通量計(jì)算結(jié)果與精確解數(shù)值。

從圖2以及表1與表2中可以得出，在該分析方法下旋轉(zhuǎn)L 形邊界區(qū)域域內(nèi)點(diǎn)的計(jì)算結(jié)果和精確解不僅在變化趨勢(shì)上基本保持一致，而且其具體計(jì)算數(shù)值與精確解下的精確值較為吻合，誤差較小，充分說(shuō)明該方法的準(zhǔn)確度較高。

表1 內(nèi)部取樣點(diǎn)Q1 和Q2 上聲壓計(jì)算值與精確值Table 1 Calculated and exact values of sound pressure at internal sampling points Q1 and Q2

表2 內(nèi)部取樣點(diǎn)Q1 和Q2 上通量計(jì)算值與精確值Table 2 Calculated and exact values of flux at internal sampling points Q1 and Q2

4.2 凸形邊界區(qū)域內(nèi)部聲壓變化

在一個(gè)凸形邊界區(qū)域條件下進(jìn)行內(nèi)點(diǎn)聲壓變化情況的分析，其邊界節(jié)點(diǎn)之間間距同樣設(shè)為0.1 cm，區(qū)域聲速設(shè)為1 個(gè)標(biāo)準(zhǔn)單位，計(jì)算求解過(guò)程中一共采用了80 個(gè)邊界節(jié)點(diǎn)，凸形邊界區(qū)域的節(jié)點(diǎn)分布情況如圖3所示。

圖3 凸形邊界節(jié)點(diǎn)分布圖Fig.3 Distribution diagram of convex boundary nodes

圖3所示凸形邊界上聲壓分布按如下公式給出：

在y=0 的邊界上，

在x=2.4 的邊界上，

在y=0.8 的邊界上，

在x=1.6 的邊界上，

在y=1.6 的邊界上，

在x=0.8 的邊界上，

在x=0 的邊界上，

此域內(nèi)聲壓精確解為

計(jì)算時(shí)間區(qū)域?yàn)?0.5~4.0 s，取兩個(gè)域內(nèi)點(diǎn)Q3(1.5,0.4)、Q4(1.6, 0.4) 作為該邊界條件下的觀察點(diǎn)，在頻率區(qū)間 1~4 內(nèi)等間隔取 7 個(gè)頻率點(diǎn)計(jì)算。分析Q3及Q4聲壓值隨時(shí)間變化的過(guò)程，比較結(jié)果如圖 4所示。表3中列出Q3和Q4的部分計(jì)算結(jié)果與精確解的對(duì)比。

由圖4和表3可知，在凸形邊界條件下選取的兩個(gè)域內(nèi)取樣點(diǎn)上聲壓計(jì)算值與精確值吻合情況也是較為理想的。

表3 內(nèi)部取樣點(diǎn)Q3 和Q4 上聲壓計(jì)算值與精確值Table 3 Calculated and exact values of sound pressure at internal sampling points Q3 and Q4

圖4 內(nèi)部觀察點(diǎn)Q3、Q4 聲壓值隨時(shí)間變化情況Fig.4 Variation of sound pressure value at internal observation point Q3 and Q4 with time change

比較旋轉(zhuǎn)L 形與凸形算例的結(jié)果，發(fā)現(xiàn)在不同邊界條件下其域內(nèi)點(diǎn)的計(jì)算結(jié)果都能夠和精確解基本保持一致，而且在合理的采樣區(qū)間內(nèi)能達(dá)到較高的計(jì)算精度，這充分說(shuō)明了所提方法的可行性以及其具有較高的時(shí)域解計(jì)算準(zhǔn)確度。

5 結(jié)語(yǔ)

針對(duì)含有復(fù)雜頻率成分聲源的輻射問(wèn)題，采用離散傅里葉反變換和頻域邊界元法相結(jié)合，研究了一種二維多頻聲輻射問(wèn)題的頻域邊界元分析方法。采用傅里葉變換將時(shí)域描述的波動(dòng)方程轉(zhuǎn)化為Helmholtz 方程，并選取多個(gè)等間隔頻率點(diǎn)作為采樣頻率，應(yīng)用邊界單元法求解各特征的Helmholtz 方程，獲得不同位置在各采樣頻率下的聲壓，采用離散傅里葉反變換將頻域聲壓幅值和相位轉(zhuǎn)化成時(shí)域聲壓。兩個(gè)數(shù)值算例分析下的域內(nèi)聲壓計(jì)算值和精確解吻合較好，誤差較小，證明該方法具有可行性，并具有較高的時(shí)域解計(jì)算準(zhǔn)確度。