不同二維聲學(xué)邊界單元在敏感度分析上的數(shù)值表現(xiàn)

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

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

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

最后本文使用MMA方法進行典型聲屏障結(jié)構(gòu)優(yōu)化分析，以驗證本文發(fā)展算法的有效性。

1 二維聲學(xué)邊界元

1.1 無奇異邊界積分方程

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

(1)

(2)

式中：

(3)

(4)

(5)

(6)

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

(7)

(8)

(9)

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

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

(1)常量單元

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

(2)線性單元

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

(10)

式中ζ表示單元內(nèi)部局部坐標(biāo)，β表示非連續(xù)單元內(nèi)部插值節(jié)點的局部坐標(biāo)。當(dāng)β=1時，式(10)為連續(xù)線性單元的插值函數(shù)表達式。

(3)二次單元

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

(11)

當(dāng)β=1時，上式為連續(xù)二次單元的插值形函數(shù)。

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

(1) 常單元

(12)

式中，系數(shù)C1和C2分別表達如下：

(13)

(14)

(2)線性單元：

(15)

式中， xl表示邊界單元中不同于節(jié)點xj的另外一節(jié)點。α表示節(jié)點xj的局部坐標(biāo)，當(dāng)β≠1 時，β=|α|，上式中的系數(shù)可表達如下：

(16)

(17)

(18)

(19)

式中：

(20)

(21)

(22)

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

(23)

(24)

(25)

(26)

(3)3二次單元：

(27)

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

(28)

(29)

式中，

(31)

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

(32)

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

(33)

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

聲學(xué)設(shè)計敏感度分析能夠提供結(jié)構(gòu)形狀變化和結(jié)構(gòu)聲學(xué)表現(xiàn)之間的關(guān)系，因此敏感度分析時結(jié)構(gòu)聲學(xué)優(yōu)化設(shè)計中極其重要的一個環(huán)節(jié)。任意目標(biāo)函數(shù)W可以定義為：

(34)

(35)

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

(36)

(37)

(38)

(39)

(40)

伴隨方程可以定義為：

(42)

根據(jù)下面邊界條件：

(43)

(44)

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

同樣的，

(46)

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

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

3 算例分析

3.1 單元類型和誤差定義

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

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

3.2 無限長圓柱殼散射算例

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

(48)

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

(49)

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

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

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

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

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

圖4 DBE22單元面誤差隨節(jié)點位置參數(shù)變化關(guān)系Fig.4 Surface error in terms of nodal position for DBE22 element

圖5 DBE33單元面誤差隨節(jié)點位置參數(shù)變化關(guān)系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給出了相似自由度下不同類型單元數(shù)值精度對比。觀察該圖可以發(fā)現(xiàn)，同等階次的非連續(xù)單元比連續(xù)單元表現(xiàn)的更有效，比如DBE33比CBE33精度高，DBE22比CBE22精度高。另外，通過對比DBE33與DBE22可以發(fā)現(xiàn)，二次單元比線性單元表現(xiàn)更好。因此，同等自由度下二次非連續(xù)單元表現(xiàn)最有效。

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

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

圖7 DBE22單元面誤差隨節(jié)點位置參數(shù)變化關(guān)系Fig.7 Error in terms of nodal position for DBE22 element

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

(50)

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

圖8 半徑r=2m圓環(huán)上計算點聲壓幅值敏感度分布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單元用于數(shù)值計算，波數(shù)為k=4。圖10中的圓柱殼分布在1m×1m的方格子里，每個圓環(huán)被離散成1000個單元，總共形成4×105個單元，快速多極算法被用于加速邊界元系數(shù)矩陣與向量相乘的計算和系統(tǒng)方程的求解。這兩幅圖顯示出在進行大規(guī)模多域問題的分析時，本文發(fā)展的二維快速敏感度算法具有較好的應(yīng)用潛力。

圖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型聲屏障結(jié)構(gòu)優(yōu)化分析

隨著我國汽車保有量的迅速增加，引發(fā)的交通噪聲污染問題也日益嚴(yán)重。因此，降低降低噪聲是一項越來越引起重視的研究課題。在噪聲源和居民區(qū)之間放置聲屏障可以有效的降低保護區(qū)內(nèi)的聲壓級[13]。雖然聲屏障的高度對其降噪效果影響很大，然而考慮到車內(nèi)人員的視覺感受和建筑成本，聲屏障的高度是受到限制的。在高度一定情況下，聲屏障頂端結(jié)構(gòu)對其降噪表現(xiàn)也有很大的影響，常用的頂端結(jié)構(gòu)是Y型。本文對典型的Y型聲屏障進行二維優(yōu)化分析，以得到降噪效果更好的頂端結(jié)構(gòu)外形。如圖11所示，聲源距離地面高1.0m，離聲屏障10.4m，Y型聲屏障頂端左右兩邊角度分別為θ1和θ2，長度分別為l2和l1，低端高為3.0m，分布在聲隱區(qū)的6個計算點坐標(biāo)分別為(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型聲屏障進行優(yōu)化分析，目標(biāo)函數(shù)為6個計算點處的平均聲壓級，設(shè)計變量為頂端角度和長度。約束條件如下：

(51)

設(shè)定設(shè)計變量初始值為

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

圖12 設(shè)計變量及目標(biāo)函數(shù)隨迭代次數(shù)變化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 總結(jié)

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

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