二維多頻聲輻射問題的頻域邊界元分析方法

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

（1.湖南工業大學 機械工程學院，湖南 株洲 412007；2.衢州學院 機械工程學院，浙江 衢州 324000）

1 背景知識

離心風機出口氣動噪聲、機械設備運轉時摩擦噪聲等聲源的主要頻率成分在頻譜上具有一定寬度，而多頻聲源的輻射問題一直是工程師關注的重要問題。適用于聲輻射分析的工程數值方法以邊界元法（boundary element method，BEM）[1-3]、有限元法（finite element analysis，FEM）[4]和矩量法（method of moments，MoM）[5]為主要代表，其中BEM 和MoM 都以邊界積分方程為理論基礎，特別適用于無限域的聲輻射分析。

邊界元法是一種基于邊界積分方程為數學基礎的數值分析方法，其邊界積分方程能夠自動滿足無窮遠場的輻射條件，并且僅需要對邊界進行離散。較之其他的數值計算方法，邊界元法引入了基本解，使其具有了解析和數值相結合的特點，從而使其計算精度相對較高，是一種比較適用于分析聲學問題的計算方法[6]，目前被廣泛應用于求解裂紋問題[7]、擴散波問題[8]、熱傳導問題[9]等。

目前對聲學問題的分析主要采用頻域法[10]和時域法[11]，時域類方法具有直接且通用性好的優點，然而在選擇時間步長時會遇到數值不穩定問題。頻域法通過Helmholtz 控制方程和基本解，推導出邊界積分方程，再將邊界積分方程進行線性離散，建立矩陣，求解線性方程組來獲得相應的解[12]。頻域法在頻率成分預判較準確的基礎上能夠獲得較高的計算精度，通過反變換獲得時域解。

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

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

2 二維聲學頻域邊界元法

在均勻理想介質E中，二維聲學問題在空間內傳播的波動方程表示形式為

x為配置點空間坐標；

c為介質中聲音傳播的速度。

對式（1）進行Fourier 變換，假設時間項選取為e-jωt，可將瞬態問題從時間域轉化到頻率域進行簡單穩態分析，則頻域聲壓可以表示為

j 為虛數單位，j2=-1；ω為聲源簡諧振動角頻率。

將式（2）代入（1），消去時間項，即可得如下Helmholtz 控制方程。

式中k為波數，且k=ω/c。

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

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

式中符號ˉ表示已知值。

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

式中：q(x)為聲壓在x處的法向導數；n為單位外法向，指向背離聲域；ρ為聲學介質的質量密度；vn為法向速度。

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

式中Z為聲阻抗。

將方程（3）轉化為等效積分形式，并將試函數取為二維聲學問題的基本解：

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

y為場點空間坐標；

q(y)為聲壓在y處的法向導數；

Sy為結構邊界S的子邊界；

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

r為源點與場點之間的距離，且r=|x-y|；

x為配置點（即計算點或源點）；

y為場點。

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

式中：u為任一邊界的節點，u=1, 2, …,Ne；

Sv為邊界上的單元v；

qv為在單元v上的法向導數；

將源點置于所有場點上可得到：

式中Huv和Guv均為影響系數矩陣，為遍及所有節點單元上的積分總和。

也可將式（1）寫為

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

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

式中：A為未知聲壓值和位置法向導數組成的Ne×Ne系數矩陣；

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

h為由包含未知聲壓值和未知法向導數組成的Ne階向量。

3 離散傅里葉逆變換

為實現離散傅里葉反變換運算，需將連續信號進行截斷、采樣。

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

式中：Δt為時間步長，其中Δt=T/N；

ω0為初始頻率點。

為避免產生頻率混疊現象，信號采樣頻率

同時為避免頻混現象，可通過選擇較小的采樣間隔T（即較高的采樣頻率）來減小這種誤差。

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

其Fourier 變換：

根據δ函數的篩選性，化簡求得其Fourier 系數Ck為

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

一個周期內N個采樣點的復數值為

式中N為采樣頻率點的數量。

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

傅里葉級數可以寫成復指數形式，根據歐拉公式

可得到：

4 數值算例

為驗證二維頻域邊界元法結合離散傅里葉反變換求解多頻聲輻射問題的準確性，設計了兩個驗證算例，將計算結果與解析解進行對比以驗證該方法的正確性。

4.1 旋轉L 形邊界區域內部聲壓變化

在一個旋轉L 形邊界區域條件下進行內點聲壓變化情況分析，其邊界節點之間間距設為0.1 cm，x、y分別為內部點的橫坐標與縱坐標（算例中都將如此使用），計算求解過程中一共采用了64 個邊界節點，旋轉L 形邊界區域節點分布情況如圖1所示。

圖1 旋轉L 形邊界節點分布圖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 的邊界上，

域內聲壓精確解為

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

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

從圖2中可以看到，在該方法下兩個取樣點上的聲壓計算值與精確值吻合較好。為更準確直觀地展示計算結果和精確解的差異，于表1~2 中列出了隨機選取的部分時間點上，對于Q1和Q2兩個內部觀察點的聲壓以及通量計算結果與精確解數值。

從圖2以及表1與表2中可以得出，在該分析方法下旋轉L 形邊界區域域內點的計算結果和精確解不僅在變化趨勢上基本保持一致，而且其具體計算數值與精確解下的精確值較為吻合，誤差較小，充分說明該方法的準確度較高。

表1 內部取樣點Q1 和Q2 上聲壓計算值與精確值Table 1 Calculated and exact values of sound pressure at internal sampling points Q1 and Q2

表2 內部取樣點Q1 和Q2 上通量計算值與精確值Table 2 Calculated and exact values of flux at internal sampling points Q1 and Q2

4.2 凸形邊界區域內部聲壓變化

在一個凸形邊界區域條件下進行內點聲壓變化情況的分析，其邊界節點之間間距同樣設為0.1 cm，區域聲速設為1 個標準單位，計算求解過程中一共采用了80 個邊界節點，凸形邊界區域的節點分布情況如圖3所示。

圖3 凸形邊界節點分布圖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 的邊界上，

此域內聲壓精確解為

計算時間區域為 0.5~4.0 s，取兩個域內點Q3(1.5,0.4)、Q4(1.6, 0.4) 作為該邊界條件下的觀察點，在頻率區間 1~4 內等間隔取 7 個頻率點計算。分析Q3及Q4聲壓值隨時間變化的過程，比較結果如圖 4所示。表3中列出Q3和Q4的部分計算結果與精確解的對比。

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

表3 內部取樣點Q3 和Q4 上聲壓計算值與精確值Table 3 Calculated and exact values of sound pressure at internal sampling points Q3 and Q4

圖4 內部觀察點Q3、Q4 聲壓值隨時間變化情況Fig.4 Variation of sound pressure value at internal observation point Q3 and Q4 with time change

比較旋轉L 形與凸形算例的結果，發現在不同邊界條件下其域內點的計算結果都能夠和精確解基本保持一致，而且在合理的采樣區間內能達到較高的計算精度，這充分說明了所提方法的可行性以及其具有較高的時域解計算準確度。

5 結語

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