(信陽(yáng)師范學(xué)院 建筑與土木工程學(xué)院,河南 信陽(yáng) 464000)
優(yōu)化方法可分為直接方法和近似方法[1],直接方法的代表有遺傳算法和退火算法[2-3],這類(lèi)方法是通過(guò)合適的參數(shù)變化尋找全局最優(yōu)值,但該類(lèi)方法往往需要大量的函數(shù)計(jì)算。作為對(duì)比,近似方法只需少量的函數(shù)計(jì)算就可以達(dá)到收斂,但該方法往往只收斂到局部最小值。大多數(shù)局部近似方法都是基于目標(biāo)函數(shù)在某一設(shè)計(jì)點(diǎn)附件的泰勒展開(kāi)進(jìn)行尋優(yōu)計(jì)算的,因此在每個(gè)迭代步都需要進(jìn)行敏感度計(jì)算,這一步往往是十分耗時(shí)的。有限差分法,因其很容易的被實(shí)施,已被廣泛應(yīng)用于敏感度計(jì)算,然而該方法計(jì)算精度不高,并對(duì)每一個(gè)設(shè)計(jì)變量都需要進(jìn)行差分計(jì)算,計(jì)算量大。解析與半解析方法可以有效用于敏感度計(jì)算,例如直接微分法[4-6]。該方法將目標(biāo)函數(shù)對(duì)設(shè)計(jì)變量直接進(jìn)行微分計(jì)算,計(jì)算精度高,然而對(duì)于多設(shè)計(jì)變量的優(yōu)化分析,在每一個(gè)迭代步對(duì)每一個(gè)設(shè)計(jì)變量都要進(jìn)行微分計(jì)算,這會(huì)大大降低計(jì)算效率。伴隨變量法[7],通過(guò)引入與設(shè)計(jì)變量無(wú)關(guān)的伴隨方程,在進(jìn)行敏感度分析時(shí),該伴隨方程只需求解一次,然后間接計(jì)算目標(biāo)函數(shù)的敏感度值,因此有效提高了計(jì)算效率。本文推導(dǎo)出基于伴隨變量法的聲學(xué)邊界元敏感度表達(dá)式,然后采用不同類(lèi)型邊界單元進(jìn)行聲學(xué)敏感度計(jì)算,并對(duì)比其計(jì)算效率。
實(shí)際上在早期的積分方程應(yīng)用中,包括最近的許多研究[8-10],常量單元得到了廣泛的應(yīng)用,因其實(shí)施起來(lái)簡(jiǎn)單。但其計(jì)算精度差,在實(shí)際工程計(jì)算中,為了達(dá)到一定的精度要求,往往要加大離散,增加了計(jì)算量。很多學(xué)者使用連續(xù)線性和二次單元進(jìn)行數(shù)值計(jì)算,以提高計(jì)算精度。然而連續(xù)單元對(duì)網(wǎng)格的適應(yīng)性不高,在角點(diǎn)處需要特殊處理,難以進(jìn)行復(fù)雜邊界的積分計(jì)算,然而高階非連續(xù)單元可以很好的克服這些問(wèn)題。文獻(xiàn)[8]和文獻(xiàn)[9]分別給出了二維和三維非連續(xù)高階單元的計(jì)算過(guò)程。非連續(xù)單元的插值節(jié)點(diǎn)位于單元內(nèi)部,插值型函數(shù)的表達(dá)依賴(lài)節(jié)點(diǎn)位置。文獻(xiàn)[10]給出了不同類(lèi)型連續(xù)與非連續(xù)邊界單元的計(jì)算精度對(duì)比,并研究了單元誤差與頻率、網(wǎng)格尺寸和插值節(jié)點(diǎn)位置的關(guān)系,并得到在相似自由度下非連續(xù)元比連續(xù)元表現(xiàn)更有效的結(jié)論。上述文獻(xiàn)都是基于傳統(tǒng)邊界元法進(jìn)行不同類(lèi)型單元的計(jì)算精度比較的,使用傳統(tǒng)邊界元法進(jìn)行外聲場(chǎng)分析時(shí),在某些虛假本征頻率處可能得不到精確結(jié)果。另外,這些文獻(xiàn)都沒(méi)有進(jìn)行敏感度分析。本文使用Burtion-Miller[11]法克服虛假本征頻率問(wèn)題,首先推導(dǎo)出基于不同類(lèi)型邊界單元離散下的無(wú)奇異Burtion-Miller邊界積分表達(dá)式,然后使用帶有解析解的簡(jiǎn)單算例對(duì)比不同類(lèi)型邊界單元在進(jìn)行聲學(xué)敏感度分析時(shí)的計(jì)算表現(xiàn),并研究在進(jìn)行敏感度分析時(shí)的數(shù)值誤差隨非連續(xù)單元內(nèi)部插值節(jié)點(diǎn)位置變化關(guān)系,以考察得到優(yōu)化節(jié)點(diǎn)位置參數(shù)值。
最后本文使用MMA方法進(jìn)行典型聲屏障結(jié)構(gòu)優(yōu)化分析,以驗(yàn)證本文發(fā)展算法的有效性。
流體中聲場(chǎng)分布滿足Helmholtz控制方程,通過(guò)格林公式變換后,可得到如下的邊界積分方程:
(1)
(2)
式中:
(3)
(4)
(5)
(6)
如果在點(diǎn)處邊界光滑,系數(shù)為1/2。單獨(dú)使用式(5)或(6)求解外聲場(chǎng)問(wèn)題,在虛擬頻率處會(huì)產(chǎn)生非唯一解問(wèn)題,通過(guò)組合這兩個(gè)方程構(gòu)成Burton-Miller表達(dá)式可以有效克服這一問(wèn)題[11-12]。然而式(5)和(6)含有奇異積分項(xiàng),為了得到精確結(jié)果必須對(duì)此做特殊處理,本文使用Cauchy主值積分和Hadamard有限積分法直接精確的計(jì)算奇異積分。接下來(lái)給出不同類(lèi)型單元離散下的,無(wú)奇異積分表達(dá)式:
(7)
(8)
(9)
式中m表示每個(gè)單元的插值節(jié)點(diǎn)數(shù)目,Φ表示插值函數(shù)。函數(shù)f(x)可以為聲壓 φ(x)或聲壓通量q(x)等。
接下來(lái)首先給出不同類(lèi)型邊界單元的插值形函數(shù)表達(dá)式:
(1)常量單元
對(duì)于常量單元來(lái)說(shuō),滿足m=1,Φ1=1。
(2)線性單元
對(duì)于線性單元來(lái)說(shuō),插值形函數(shù)Φ表達(dá)為:
(10)
式中ζ表示單元內(nèi)部局部坐標(biāo),β表示非連續(xù)單元內(nèi)部插值節(jié)點(diǎn)的局部坐標(biāo)。當(dāng)β=1時(shí),式(10)為連續(xù)線性單元的插值函數(shù)表達(dá)式。
(3)二次單元
對(duì)于二次單元來(lái)說(shuō)m=3,插值函數(shù)Φ如下表達(dá):
(11)
當(dāng)β=1時(shí),上式為連續(xù)二次單元的插值形函數(shù)。
使用不同單元離散邊界時(shí),可得到方程(7)與(8)中的不同b(xj)和D(xj)的表達(dá)式,分別如下:
(1) 常單元
(12)
式中,系數(shù)C1和C2分別表達(dá)如下:
(13)
(14)
(2)線性單元:
(15)
式中, xl表示邊界單元中不同于節(jié)點(diǎn)xj的另外一節(jié)點(diǎn)。α表示節(jié)點(diǎn)xj的局部坐標(biāo),當(dāng)β≠1 時(shí),β=|α|,上式中的系數(shù)可表達(dá)如下:
(16)
(17)
(18)
(19)
式中:
(20)
(21)
(22)
當(dāng)β=1時(shí),非連續(xù)線性單元退化為連續(xù)線性單元,式(15)中的系數(shù)表達(dá)式變化為:
(23)
(24)
(25)
(26)
(3)3二次單元:
(27)
式中,xl和xm表示邊界單元中除了源點(diǎn)xj外另兩個(gè)插值節(jié)點(diǎn)。當(dāng)β≠1 時(shí),式(27)中系數(shù)表達(dá)式為:
(28)
(29)
式中,
(31)
在此,我們推出了基于不同類(lèi)型邊界單元離散下的無(wú)奇異邊界積分表達(dá)式,離散并組合(7)式與(8)式后,可得到如下的線性系統(tǒng)方程:
(32)
將未知量轉(zhuǎn)移到方程左邊,已知量轉(zhuǎn)移到方程右邊,式(32)可重新表達(dá)如下:
(33)
聲學(xué)設(shè)計(jì)敏感度分析能夠提供結(jié)構(gòu)形狀變化和結(jié)構(gòu)聲學(xué)表現(xiàn)之間的關(guān)系,因此敏感度分析時(shí)結(jié)構(gòu)聲學(xué)優(yōu)化設(shè)計(jì)中極其重要的一個(gè)環(huán)節(jié)。任意目標(biāo)函數(shù)W可以定義為:
(34)
(35)
同時(shí)邊界和域的微分表達(dá)式可寫(xiě)為:
(36)
(37)
(38)
(39)
(40)
伴隨方程可以定義為:
(42)
根據(jù)下面邊界條件:
(43)
(44)
將式(42)、(43)和(44)代入式(41)中,可以得到下面表達(dá)式:
同樣的,
(46)
將式(46)代入式(45),可以得到:
可以看出式(47)不含有邊界上未知量的敏感度。邊界聲壓的梯度和伴隨函數(shù) 被引入到目標(biāo)函數(shù)的敏感度方程中。
非連續(xù)單元的插值節(jié)點(diǎn)位于單元內(nèi)部,插值節(jié)點(diǎn)位置參數(shù)值決定了插值形函數(shù)的表達(dá),并決定了數(shù)值計(jì)算精度,因此研究節(jié)點(diǎn)位置對(duì)計(jì)算精度的影響以得到優(yōu)化節(jié)點(diǎn)位置參數(shù)值是十分重要的。在圖1中'DBEmn'和'CBEmn'分別表示帶有'm'個(gè)幾何插值點(diǎn)和'n'個(gè)物理插值點(diǎn)的非連續(xù)和連續(xù)單元類(lèi)型。例如,圖中'DBE21'表示線性幾何近似的常量邊界單元,'DBE22'表示非連續(xù)線性邊界單元,'CBE22'表示連續(xù)線性邊界單元,而'CBE33'表示連續(xù)二次邊界單元,'DBE33'表示非連續(xù)二次邊界單元。對(duì)于三維不同連續(xù)和非連續(xù)單元類(lèi)型的比較可以參照文獻(xiàn)[10]。常單元的物理插值節(jié)點(diǎn)位于單元中心處,而非連續(xù)單元的物理插值節(jié)點(diǎn)位置由參數(shù) (0< <1)決定,見(jiàn)圖1。另外,本文使用基于復(fù)數(shù)形式的聲壓面誤差來(lái)表征數(shù)值計(jì)算精度,對(duì)于面誤差的詳細(xì)定義,請(qǐng)參考文獻(xiàn)[10]。
圖1 不同類(lèi)型二維邊界單元的幾何與物理插值點(diǎn)分布Fig.1 Distribution of geometrical nodes and interpolation nodes
本小節(jié)進(jìn)行無(wú)限長(zhǎng)圓柱殼在平面波入射時(shí)的散射聲場(chǎng)分析,平面波入射方向垂直于圓柱殼軸線,該模型可以簡(jiǎn)化成二維問(wèn)題,平面波幅值為1Pa,圓柱殼半徑 m。文獻(xiàn)給出了該模型的聲場(chǎng)解析表達(dá),如下:
(48)
式中 代表紐曼符號(hào), 代表對(duì) 進(jìn)行求導(dǎo)。選取圓柱半徑 為設(shè)計(jì)變量,將上式對(duì)設(shè)計(jì)變量進(jìn)行微分,可得到如下聲壓敏感度表達(dá)式:
(49)
圖2給出相似自由度下不同類(lèi)型邊界單元計(jì)算精度的對(duì)比,用于計(jì)算面誤差的曲線選取為距圓心 處的圓環(huán)上,計(jì)算頻率為5000Hz。觀察該圖可以發(fā)現(xiàn),面誤差隨著計(jì)算節(jié)點(diǎn)數(shù)增加而降低,連續(xù)線性單元CBE22精度最差,非連續(xù)二次單元DBE33精度最高,非連續(xù)單元比連續(xù)單元表現(xiàn)更好。
圖2 不同類(lèi)型單元計(jì)算精度對(duì)比Fig.2 Comparison of numerical performance for different types of boundary elements
圖3給出了分別使用單Helmholtz邊界積分方程(5)式(亦稱(chēng)作CBIE方程)和Burton-Miller方程進(jìn)行數(shù)值計(jì)算得到的結(jié)果與解析解的對(duì)比,DBE33單元用于數(shù)值離散,離散單元數(shù)為10000.觀察該圖可以發(fā)現(xiàn)使用CBIE方法在某些虛假頻率處計(jì)算結(jié)果偏離解析解,然后使用Burton-Miller方法得到的結(jié)果與解析解符合一致。
對(duì)于非連續(xù)單元,選取合適的節(jié)點(diǎn)位置參數(shù)值以提高單元在計(jì)算精度方面的表現(xiàn)是十分關(guān)鍵的。圖4呈現(xiàn)非連續(xù)線性邊界單元DBE22插值節(jié)點(diǎn)位置參數(shù)值對(duì)計(jì)算精度的影響,依次選取4個(gè)計(jì)算頻率進(jìn)行數(shù)值誤差分析,離散單元數(shù)為10000。觀察該圖可以發(fā)現(xiàn),隨著計(jì)算頻率的增加面誤差相應(yīng)增加;對(duì)于每條頻率誤差曲線面誤差都是先隨著節(jié)點(diǎn)參數(shù)值增加緩慢降低,然后迅速降低,在達(dá)到最低值時(shí)迅速增加,最后又開(kāi)始緩慢增加。優(yōu)化節(jié)點(diǎn)位置參數(shù)值一致逼近0.57,該值非常接近勒讓德多項(xiàng)式的零值。圖5給出DBE33單元插值節(jié)點(diǎn)位置參數(shù)值對(duì)計(jì)算精度的影響,優(yōu)化節(jié)點(diǎn)位置參數(shù)值一致逼近0.77,該值非常接近勒讓德多項(xiàng)式的零值。
圖3 計(jì)算點(diǎn) 處聲壓敏感度實(shí)部與虛部分別隨頻率變化關(guān)系Fig.3 Acoustic pressure sensitivity at point with different frequencies
圖4 DBE22單元面誤差隨節(jié)點(diǎn)位置參數(shù)變化關(guān)系Fig.4 Surface error in terms of nodal position for DBE22 element
圖5 DBE33單元面誤差隨節(jié)點(diǎn)位置參數(shù)變化關(guān)系Fig.5 Error in terms of nodal position for DBE33 element
圖6 不同單元聲壓敏感度誤差對(duì)比Fig.6 Surface error of acoustic pressure sensitivity for different types of elements
圖6給出了相似自由度下不同類(lèi)型單元數(shù)值精度對(duì)比。觀察該圖可以發(fā)現(xiàn),同等階次的非連續(xù)單元比連續(xù)單元表現(xiàn)的更有效,比如DBE33比CBE33精度高,DBE22比CBE22精度高。另外,通過(guò)對(duì)比DBE33與DBE22可以發(fā)現(xiàn),二次單元比線性單元表現(xiàn)更好。因此,同等自由度下二次非連續(xù)單元表現(xiàn)最有效。
在這一小節(jié)里,我們進(jìn)行多級(jí)散射問(wèn)題分析。首先考慮四個(gè)無(wú)限長(zhǎng)圓柱體在平面波作用下的聲散射問(wèn)題,圓柱半徑為0.2m,圓心位置和坐標(biāo)原點(diǎn)位置如圖7所示。
圖7 DBE22單元面誤差隨節(jié)點(diǎn)位置參數(shù)變化關(guān)系Fig.7 Error in terms of nodal position for DBE22 element
360個(gè)計(jì)算點(diǎn)均勻分布在以原點(diǎn)為圓心,半徑為2m的圓環(huán)上。聲壓幅值敏感度值如下表達(dá):
(50)
圖8給出在這些計(jì)算點(diǎn)處的聲壓幅值敏感度值分布,設(shè)計(jì)變量為圓柱殼的半徑,DBE33單元被用于數(shù)值計(jì)算,波數(shù)k=4。觀察該圖可以發(fā)現(xiàn),使用快速多極邊界元法得到的計(jì)算結(jié)果與傳統(tǒng)邊界元法得到的結(jié)果符合一致,表明快速算法的使用保持了傳統(tǒng)邊界元法的高計(jì)算精度特點(diǎn)。
圖8 半徑r=2m圓環(huán)上計(jì)算點(diǎ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個(gè)和400個(gè)圓柱殼在平面波入射時(shí)的散射聲壓敏感度分布 ,DBE33單元用于數(shù)值計(jì)算,波數(shù)為k=4。圖10中的圓柱殼分布在1m×1m的方格子里,每個(gè)圓環(huán)被離散成1000個(gè)單元,總共形成4×105個(gè)單元,快速多極算法被用于加速邊界元系數(shù)矩陣與向量相乘的計(jì)算和系統(tǒng)方程的求解。這兩幅圖顯示出在進(jìn)行大規(guī)模多域問(wèn)題的分析時(shí),本文發(fā)展的二維快速敏感度算法具有較好的應(yīng)用潛力。
圖9 4圓柱散射時(shí)的聲壓幅值敏感度分布Fig.9 Sensitivity field for the amplitude of sound pressure from 4 rigid cylinders scattering
圖10 400圓柱散射時(shí)的聲壓幅值敏感度分布Fig.10 Sensitivity field for the amplitude of sound pressure from 400 cylinders scattering
隨著我國(guó)汽車(chē)保有量的迅速增加,引發(fā)的交通噪聲污染問(wèn)題也日益嚴(yán)重。因此,降低降低噪聲是一項(xiàng)越來(lái)越引起重視的研究課題。在噪聲源和居民區(qū)之間放置聲屏障可以有效的降低保護(hù)區(qū)內(nèi)的聲壓級(jí)[13]。雖然聲屏障的高度對(duì)其降噪效果影響很大,然而考慮到車(chē)內(nèi)人員的視覺(jué)感受和建筑成本,聲屏障的高度是受到限制的。在高度一定情況下,聲屏障頂端結(jié)構(gòu)對(duì)其降噪表現(xiàn)也有很大的影響,常用的頂端結(jié)構(gòu)是Y型。本文對(duì)典型的Y型聲屏障進(jìn)行二維優(yōu)化分析,以得到降噪效果更好的頂端結(jié)構(gòu)外形。如圖11所示,聲源距離地面高1.0m,離聲屏障10.4m,Y型聲屏障頂端左右兩邊角度分別為θ1和θ2,長(zhǎng)度分別為l2和l1,低端高為3.0m,分布在聲隱區(qū)的6個(gè)計(jì)算點(diǎn)坐標(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]對(duì)Y型聲屏障進(jìn)行優(yōu)化分析,目標(biāo)函數(shù)為6個(gè)計(jì)算點(diǎn)處的平均聲壓級(jí),設(shè)計(jì)變量為頂端角度和長(zhǎng)度。約束條件如下:
(51)
設(shè)定設(shè)計(jì)變量初始值為
圖12依次給出設(shè)計(jì)變量和目標(biāo)函數(shù)隨迭代次數(shù)變化關(guān)系,可以看出在第500次迭代后收斂的比較好。優(yōu)化后的角度值為θ1=27°和θ1=¥2°,目標(biāo)函數(shù)優(yōu)化值為60.42dB。在優(yōu)化分析過(guò)程中,為了保證頂端結(jié)構(gòu)可被用于邊界元離散,強(qiáng)迫設(shè)置邊長(zhǎng)大于0.1m,觀察圖12中的第二幅圖可以發(fā)現(xiàn),頂端邊長(zhǎng)的優(yōu)化值趨向于滿足l1=0和l2=3,該結(jié)果表明在滿足一定材料成本限制時(shí),Y型聲屏障的優(yōu)化后頂端結(jié)構(gòu)外形是半Y型,實(shí)際上半Y型聲屏已被廣泛應(yīng)用于工程應(yīng)用中。圖13給出半Y型聲屏障周?chē)穆晧杭?jí)分布,頻率為200Hz。使用基于MMA的優(yōu)化分析時(shí),在每一個(gè)迭代步中都需要進(jìn)行目標(biāo)函數(shù)的敏感度計(jì)算,實(shí)際上這一步是十分耗時(shí)的,本文采用伴隨變量法進(jìn)行敏感度計(jì)算并使用FMM算法加速優(yōu)化分析,使得本文發(fā)展算法對(duì)大規(guī)模實(shí)際問(wèn)題有較高的應(yīng)用潛力。
圖12 設(shè)計(jì)變量及目標(biāo)函數(shù)隨迭代次數(shù)變化Fig.12 Values of design variables and objective function with number of iteration
圖13 半Y型聲屏障周?chē)穆晧杭?jí)分布Fig.13 SPL contour plot of half-Y noise barrier
本文首先給出基于Burton-Miller方法的不同類(lèi)型單元無(wú)奇異計(jì)算表達(dá)式,然后給出基于伴隨變量法的聲學(xué)敏感度方程。通過(guò)帶有解析解的數(shù)值算例對(duì)比不同類(lèi)型邊界單元的計(jì)算精度,發(fā)現(xiàn)非連續(xù)單元比連續(xù)單元表現(xiàn)的更有效,同時(shí)非連續(xù)單元的優(yōu)化節(jié)點(diǎn)位置參數(shù)值非常接近勒讓德多項(xiàng)式的零值。最后通過(guò)多極算射例子和聲屏障算例驗(yàn)證本文發(fā)展算法的有效性,尤其在結(jié)構(gòu)優(yōu)化分析上的高應(yīng)用潛力。
致謝:國(guó)家自然科學(xué)基金項(xiàng)目(11172291);高等學(xué)校博士學(xué)科點(diǎn)專(zhuān)項(xiàng)科研基金項(xiàng)目(20133402110036);中國(guó)科學(xué)技術(shù)大學(xué)校青年創(chuàng)新基金項(xiàng)目(WK2090000007)。