不同二維聲學邊界單元在敏感度分析上的數值表現

(信陽師范學院 建筑與土木工程學院，河南 信陽 464000)

優化方法可分為直接方法和近似方法[1]，直接方法的代表有遺傳算法和退火算法[2-3]，這類方法是通過合適的參數變化尋找全局最優值，但該類方法往往需要大量的函數計算。作為對比，近似方法只需少量的函數計算就可以達到收斂，但該方法往往只收斂到局部最小值。大多數局部近似方法都是基于目標函數在某一設計點附件的泰勒展開進行尋優計算的，因此在每個迭代步都需要進行敏感度計算，這一步往往是十分耗時的。有限差分法，因其很容易的被實施，已被廣泛應用于敏感度計算，然而該方法計算精度不高，并對每一個設計變量都需要進行差分計算，計算量大。解析與半解析方法可以有效用于敏感度計算，例如直接微分法[4-6]。該方法將目標函數對設計變量直接進行微分計算，計算精度高，然而對于多設計變量的優化分析，在每一個迭代步對每一個設計變量都要進行微分計算，這會大大降低計算效率。伴隨變量法[7]，通過引入與設計變量無關的伴隨方程，在進行敏感度分析時，該伴隨方程只需求解一次，然后間接計算目標函數的敏感度值，因此有效提高了計算效率。本文推導出基于伴隨變量法的聲學邊界元敏感度表達式，然后采用不同類型邊界單元進行聲學敏感度計算，并對比其計算效率。

實際上在早期的積分方程應用中，包括最近的許多研究[8-10]，常量單元得到了廣泛的應用，因其實施起來簡單。但其計算精度差，在實際工程計算中，為了達到一定的精度要求，往往要加大離散，增加了計算量。很多學者使用連續線性和二次單元進行數值計算，以提高計算精度。然而連續單元對網格的適應性不高，在角點處需要特殊處理，難以進行復雜邊界的積分計算，然而高階非連續單元可以很好的克服這些問題。文獻[8]和文獻[9]分別給出了二維和三維非連續高階單元的計算過程。非連續單元的插值節點位于單元內部，插值型函數的表達依賴節點位置。文獻[10]給出了不同類型連續與非連續邊界單元的計算精度對比，并研究了單元誤差與頻率、網格尺寸和插值節點位置的關系，并得到在相似自由度下非連續元比連續元表現更有效的結論。上述文獻都是基于傳統邊界元法進行不同類型單元的計算精度比較的，使用傳統邊界元法進行外聲場分析時，在某些虛假本征頻率處可能得不到精確結果。另外，這些文獻都沒有進行敏感度分析。本文使用Burtion-Miller[11]法克服虛假本征頻率問題，首先推導出基于不同類型邊界單元離散下的無奇異Burtion-Miller邊界積分表達式，然后使用帶有解析解的簡單算例對比不同類型邊界單元在進行聲學敏感度分析時的計算表現，并研究在進行敏感度分析時的數值誤差隨非連續單元內部插值節點位置變化關系，以考察得到優化節點位置參數值。

最后本文使用MMA方法進行典型聲屏障結構優化分析，以驗證本文發展算法的有效性。

1 二維聲學邊界元

1.1 無奇異邊界積分方程

流體中聲場分布滿足Helmholtz控制方程，通過格林公式變換后，可得到如下的邊界積分方程：

(1)

(2)

式中：

(3)

(4)

(5)

(6)

如果在點處邊界光滑，系數為1/2。單獨使用式(5)或(6)求解外聲場問題，在虛擬頻率處會產生非唯一解問題，通過組合這兩個方程構成Burton-Miller表達式可以有效克服這一問題[11-12]。然而式(5)和(6)含有奇異積分項，為了得到精確結果必須對此做特殊處理，本文使用Cauchy主值積分和Hadamard有限積分法直接精確的計算奇異積分。接下來給出不同類型單元離散下的，無奇異積分表達式：

(7)

(8)

(9)

式中m表示每個單元的插值節點數目，Φ表示插值函數。函數f(x)可以為聲壓 φ(x)或聲壓通量q(x)等。

接下來首先給出不同類型邊界單元的插值形函數表達式：

(1)常量單元

對于常量單元來說，滿足m=1,Φ1=1。

(2)線性單元

對于線性單元來說，插值形函數Φ表達為：

(10)

式中ζ表示單元內部局部坐標，β表示非連續單元內部插值節點的局部坐標。當β=1時，式(10)為連續線性單元的插值函數表達式。

(3)二次單元

對于二次單元來說m=3，插值函數Φ如下表達：

(11)

當β=1時，上式為連續二次單元的插值形函數。

使用不同單元離散邊界時，可得到方程(7)與(8)中的不同b(xj)和D(xj)的表達式，分別如下：

(1) 常單元

(12)

式中，系數C1和C2分別表達如下：

(13)

(14)

(2)線性單元：

(15)

式中， xl表示邊界單元中不同于節點xj的另外一節點。α表示節點xj的局部坐標，當β≠1 時，β=|α|，上式中的系數可表達如下：

(16)

(17)

(18)

(19)

式中：

(20)

(21)

(22)

當β=1時，非連續線性單元退化為連續線性單元，式(15)中的系數表達式變化為：

(23)

(24)

(25)

(26)

(3)3二次單元：

(27)

式中，xl和xm表示邊界單元中除了源點xj外另兩個插值節點。當β≠1 時，式(27)中系數表達式為：

(28)

(29)

式中，

(31)

在此，我們推出了基于不同類型邊界單元離散下的無奇異邊界積分表達式，離散并組合(7)式與(8)式后，可得到如下的線性系統方程：

(32)

將未知量轉移到方程左邊，已知量轉移到方程右邊，式(32)可重新表達如下：

(33)

1.2 基于伴隨變量法的聲學敏感度分析

聲學設計敏感度分析能夠提供結構形狀變化和結構聲學表現之間的關系，因此敏感度分析時結構聲學優化設計中極其重要的一個環節。任意目標函數W可以定義為：

(34)

(35)

同時邊界和域的微分表達式可寫為：

(36)

(37)

(38)

(39)

(40)

伴隨方程可以定義為：

(42)

根據下面邊界條件：

(43)

(44)

將式(42)、(43)和(44)代入式(41)中，可以得到下面表達式：

同樣的，

(46)

將式(46)代入式(45)，可以得到：

可以看出式(47)不含有邊界上未知量的敏感度。邊界聲壓的梯度和伴隨函數 被引入到目標函數的敏感度方程中。

3 算例分析

3.1 單元類型和誤差定義

非連續單元的插值節點位于單元內部，插值節點位置參數值決定了插值形函數的表達，并決定了數值計算精度，因此研究節點位置對計算精度的影響以得到優化節點位置參數值是十分重要的。在圖1中'DBEmn'和'CBEmn'分別表示帶有'm'個幾何插值點和'n'個物理插值點的非連續和連續單元類型。例如，圖中'DBE21'表示線性幾何近似的常量邊界單元，'DBE22'表示非連續線性邊界單元，'CBE22'表示連續線性邊界單元，而'CBE33'表示連續二次邊界單元，'DBE33'表示非連續二次邊界單元。對于三維不同連續和非連續單元類型的比較可以參照文獻[10]。常單元的物理插值節點位于單元中心處，而非連續單元的物理插值節點位置由參數 (0< <1)決定，見圖1。另外，本文使用基于復數形式的聲壓面誤差來表征數值計算精度，對于面誤差的詳細定義，請參考文獻[10]。

圖1 不同類型二維邊界單元的幾何與物理插值點分布Fig.1 Distribution of geometrical nodes and interpolation nodes

3.2 無限長圓柱殼散射算例

本小節進行無限長圓柱殼在平面波入射時的散射聲場分析，平面波入射方向垂直于圓柱殼軸線，該模型可以簡化成二維問題，平面波幅值為1Pa，圓柱殼半徑 m。文獻給出了該模型的聲場解析表達，如下：

(48)

式中 代表紐曼符號， 代表對 進行求導。選取圓柱半徑 為設計變量，將上式對設計變量進行微分，可得到如下聲壓敏感度表達式：

(49)

圖2給出相似自由度下不同類型邊界單元計算精度的對比，用于計算面誤差的曲線選取為距圓心 處的圓環上，計算頻率為5000Hz。觀察該圖可以發現，面誤差隨著計算節點數增加而降低，連續線性單元CBE22精度最差，非連續二次單元DBE33精度最高，非連續單元比連續單元表現更好。

圖2 不同類型單元計算精度對比Fig.2 Comparison of numerical performance for different types of boundary elements

圖3給出了分別使用單Helmholtz邊界積分方程(5)式(亦稱作CBIE方程)和Burton-Miller方程進行數值計算得到的結果與解析解的對比，DBE33單元用于數值離散，離散單元數為10000.觀察該圖可以發現使用CBIE方法在某些虛假頻率處計算結果偏離解析解，然后使用Burton-Miller方法得到的結果與解析解符合一致。

對于非連續單元，選取合適的節點位置參數值以提高單元在計算精度方面的表現是十分關鍵的。圖4呈現非連續線性邊界單元DBE22插值節點位置參數值對計算精度的影響，依次選取4個計算頻率進行數值誤差分析，離散單元數為10000。觀察該圖可以發現，隨著計算頻率的增加面誤差相應增加；對于每條頻率誤差曲線面誤差都是先隨著節點參數值增加緩慢降低，然后迅速降低，在達到最低值時迅速增加，最后又開始緩慢增加。優化節點位置參數值一致逼近0.57，該值非常接近勒讓德多項式的零值。圖5給出DBE33單元插值節點位置參數值對計算精度的影響，優化節點位置參數值一致逼近0.77，該值非常接近勒讓德多項式的零值。

圖3 計算點 處聲壓敏感度實部與虛部分別隨頻率變化關系Fig.3 Acoustic pressure sensitivity at point with different frequencies

圖4 DBE22單元面誤差隨節點位置參數變化關系Fig.4 Surface error in terms of nodal position for DBE22 element

圖5 DBE33單元面誤差隨節點位置參數變化關系Fig.5 Error in terms of nodal position for DBE33 element

圖6 不同單元聲壓敏感度誤差對比Fig.6 Surface error of acoustic pressure sensitivity for different types of elements

圖6給出了相似自由度下不同類型單元數值精度對比。觀察該圖可以發現，同等階次的非連續單元比連續單元表現的更有效，比如DBE33比CBE33精度高，DBE22比CBE22精度高。另外，通過對比DBE33與DBE22可以發現，二次單元比線性單元表現更好。因此，同等自由度下二次非連續單元表現最有效。

3.3 多個無限長圓柱殼散射算例

在這一小節里，我們進行多級散射問題分析。首先考慮四個無限長圓柱體在平面波作用下的聲散射問題，圓柱半徑為0.2m，圓心位置和坐標原點位置如圖7所示。

圖7 DBE22單元面誤差隨節點位置參數變化關系Fig.7 Error in terms of nodal position for DBE22 element

360個計算點均勻分布在以原點為圓心，半徑為2m的圓環上。聲壓幅值敏感度值如下表達：

(50)

圖8給出在這些計算點處的聲壓幅值敏感度值分布，設計變量為圓柱殼的半徑，DBE33單元被用于數值計算，波數k=4。觀察該圖可以發現，使用快速多極邊界元法得到的計算結果與傳統邊界元法得到的結果符合一致，表明快速算法的使用保持了傳統邊界元法的高計算精度特點。

圖8 半徑r=2m圓環上計算點聲壓幅值敏感度分布Fig.8 Sensitivity for the amplitude of sound pressure at the sample points distributed on a circle with radius r=2m

圖9和圖10分別給出4個和400個圓柱殼在平面波入射時的散射聲壓敏感度分布 ，DBE33單元用于數值計算，波數為k=4。圖10中的圓柱殼分布在1m×1m的方格子里，每個圓環被離散成1000個單元，總共形成4×105個單元，快速多極算法被用于加速邊界元系數矩陣與向量相乘的計算和系統方程的求解。這兩幅圖顯示出在進行大規模多域問題的分析時，本文發展的二維快速敏感度算法具有較好的應用潛力。

圖9 4圓柱散射時的聲壓幅值敏感度分布Fig.9 Sensitivity field for the amplitude of sound pressure from 4 rigid cylinders scattering

圖10 400圓柱散射時的聲壓幅值敏感度分布Fig.10 Sensitivity field for the amplitude of sound pressure from 400 cylinders scattering

3.4 Y型聲屏障結構優化分析

隨著我國汽車保有量的迅速增加，引發的交通噪聲污染問題也日益嚴重。因此，降低降低噪聲是一項越來越引起重視的研究課題。在噪聲源和居民區之間放置聲屏障可以有效的降低保護區內的聲壓級[13]。雖然聲屏障的高度對其降噪效果影響很大，然而考慮到車內人員的視覺感受和建筑成本，聲屏障的高度是受到限制的。在高度一定情況下，聲屏障頂端結構對其降噪表現也有很大的影響，常用的頂端結構是Y型。本文對典型的Y型聲屏障進行二維優化分析，以得到降噪效果更好的頂端結構外形。如圖11所示，聲源距離地面高1.0m，離聲屏障10.4m，Y型聲屏障頂端左右兩邊角度分別為θ1和θ2，長度分別為l2和l1，低端高為3.0m，分布在聲隱區的6個計算點坐標分別為(15,2)、(30,2)、(45,2)、(15,5)、(30,5)和(45,5)。

圖11 Y型聲屏障截面圖Fig.11 Cross section of the Y-shaped noise barrier

本文采用MMA方法[14]對Y型聲屏障進行優化分析，目標函數為6個計算點處的平均聲壓級，設計變量為頂端角度和長度。約束條件如下：

(51)

設定設計變量初始值為

圖12依次給出設計變量和目標函數隨迭代次數變化關系，可以看出在第500次迭代后收斂的比較好。優化后的角度值為θ1=27°和θ1=￥2°，目標函數優化值為60.42dB。在優化分析過程中，為了保證頂端結構可被用于邊界元離散，強迫設置邊長大于0.1m，觀察圖12中的第二幅圖可以發現，頂端邊長的優化值趨向于滿足l1=0和l2=3，該結果表明在滿足一定材料成本限制時，Y型聲屏障的優化后頂端結構外形是半Y型，實際上半Y型聲屏已被廣泛應用于工程應用中。圖13給出半Y型聲屏障周圍的聲壓級分布，頻率為200Hz。使用基于MMA的優化分析時，在每一個迭代步中都需要進行目標函數的敏感度計算，實際上這一步是十分耗時的，本文采用伴隨變量法進行敏感度計算并使用FMM算法加速優化分析，使得本文發展算法對大規模實際問題有較高的應用潛力。

圖12 設計變量及目標函數隨迭代次數變化Fig.12 Values of design variables and objective function with number of iteration

圖13 半Y型聲屏障周圍的聲壓級分布Fig.13 SPL contour plot of half-Y noise barrier

4 總結

本文首先給出基于Burton-Miller方法的不同類型單元無奇異計算表達式，然后給出基于伴隨變量法的聲學敏感度方程。通過帶有解析解的數值算例對比不同類型邊界單元的計算精度，發現非連續單元比連續單元表現的更有效，同時非連續單元的優化節點位置參數值非常接近勒讓德多項式的零值。最后通過多極算射例子和聲屏障算例驗證本文發展算法的有效性，尤其在結構優化分析上的高應用潛力。

致謝：國家自然科學基金項目(11172291)；高等學校博士學科點專項科研基金項目(20133402110036)；中國科學技術大學校青年創新基金項目(WK2090000007)。