陳磊磊， 申曉偉， 劉 程， 徐延明

(1. 信陽師范學(xué)院 建筑與土木工程學(xué)院，河南 信陽 464000； 2. 中國科學(xué)技術(shù)大學(xué) 近代力學(xué)系， 合肥 230026)

聲屏障作為一種有效、經(jīng)濟(jì)的降噪工具，在交通噪聲治理中已被廣泛采用，如何有效地利用好這一降噪措施，使其發(fā)揮出更大的經(jīng)濟(jì)技術(shù)效果，對(duì)于改善人們的生活質(zhì)量具有重要的意義。聲屏障的降噪效果與聲屏障的形狀、尺寸和材料屬性有關(guān)，進(jìn)行聲屏障結(jié)構(gòu)優(yōu)化設(shè)計(jì)可以有效地提高其降噪效果。目前大多研究都集中在采用遺傳算法等啟發(fā)式算法進(jìn)行聲屏障的形狀優(yōu)化例如，Baulac等[1]采用遺傳算法對(duì)T型聲屏障的頂端結(jié)構(gòu)進(jìn)行了優(yōu)化設(shè)計(jì)；Toledo等[2]采用了進(jìn)化算法對(duì)聲屏障頂端形狀進(jìn)行了優(yōu)化；Mun等[3]采用了模擬退火方法對(duì)聲屏障的幾何尺寸進(jìn)行了優(yōu)化。這些研究大多采用傳統(tǒng)幾何插值方法描述結(jié)構(gòu)形狀，在優(yōu)化過程中很難靈活地控制形狀變化，并需進(jìn)行網(wǎng)格再劃分。由于這些缺點(diǎn)，現(xiàn)有研究大多對(duì)幾何形狀較為簡單的頂端結(jié)構(gòu)進(jìn)行優(yōu)化[4]。等幾何方法(Isogeometric Analysis，IGA)可以克服這些缺點(diǎn)，該方法采用CAD(Computer Aided Design)領(lǐng)域中廣泛使用的非均勻有理B樣條(Non Uniform Rational B Sample，NURBS)進(jìn)行幾何建模，并在CAE(Computer Aided Engineering)分析中采用相同的樣條基函數(shù)進(jìn)行物理場插值計(jì)算，實(shí)現(xiàn)了幾何模型與分析模型的同一表達(dá)。與傳統(tǒng)算法對(duì)比，具有以下優(yōu)點(diǎn)：①幾何模型是精確的，避免了幾何模型誤差，提高計(jì)算精度；②幾何模型是由控制點(diǎn)決定的，避免了網(wǎng)格再劃分的工作，提高了計(jì)算效率。此外，在用于模型形狀優(yōu)化分析的敏感度求解方面也具有一些優(yōu)勢：①對(duì)具有高階效應(yīng)的復(fù)雜幾何模型的敏感度分析更精確。NURBS基函數(shù)具有高階連續(xù)特性，比拉格朗日形函數(shù)包含更多的信息，更有利于進(jìn)行結(jié)構(gòu)響應(yīng)分析和敏感度計(jì)算；②在基于敏感度分析的優(yōu)化過程中，大大簡化了修改模型與原CAD模型幾何數(shù)據(jù)的交互，由于NURBS基函數(shù)的使用，僅需簡單地移動(dòng)等幾何節(jié)點(diǎn)位置，即可實(shí)現(xiàn)分析模型的形狀優(yōu)化[5]。因此，IGA方法的建立，使復(fù)雜實(shí)際結(jié)構(gòu)的形狀拓?fù)鋬?yōu)化得以有效開展。

邊界元法已廣泛應(yīng)用于流體力學(xué)、聲學(xué)、電磁學(xué)、斷裂力學(xué)等工程領(lǐng)域。與其他數(shù)值方法相比，邊界元法(Boundary Element Method，BEM)在求解無限域或半無限域問題上，具有獨(dú)特的優(yōu)勢。近些年，一些學(xué)者將IGA與BEM結(jié)合起來，發(fā)展了等幾何邊界元法(IGABEM)，并有效應(yīng)用于位勢問題和彈性力學(xué)領(lǐng)域。例如，Gu等[6]將IGA與邊界面法相結(jié)合，進(jìn)行三維位勢問題分析；Gong等[7]將該方法應(yīng)用到三維位勢問題領(lǐng)域；Simpson等[8]詳細(xì)給出了該方法的計(jì)算流程，并應(yīng)用到二維彈性力學(xué)領(lǐng)域。近幾年，IGABEM開始被應(yīng)用到聲學(xué)領(lǐng)域。例如，Peake等[9]采用該方法進(jìn)行二維Helmholtz問題分析，隨后，Coox等[10-11]對(duì)該方法進(jìn)行深入研究，從二維到三維，從NURBS樣條到T樣條。另外，Kostas等[12]采用IGABEM進(jìn)行彈性力學(xué)問題的敏感度計(jì)算和形狀優(yōu)化分析，這些研究都顯示出IGA在形狀優(yōu)化中具有巨大的優(yōu)勢。然而，在聲學(xué)領(lǐng)域的敏感度分析及其結(jié)構(gòu)聲學(xué)優(yōu)化分析卻研究較少。對(duì)于本文所關(guān)心的外聲場問題，采用IGABEM可自動(dòng)滿足無限域處的索莫菲邊界條件，只需對(duì)結(jié)構(gòu)外表面進(jìn)行離散和數(shù)值積分計(jì)算；單元上的積分面和變量都采用NURBS基函數(shù)插值得到，但是NURBS基函數(shù)是通過遞推公式計(jì)算得到的，在邊界元系數(shù)矩陣元素計(jì)算過程中，需要增加一次循環(huán)計(jì)算，這會(huì)降低IGABEM的計(jì)算效率，導(dǎo)致相同自由度數(shù)的IGABEM比傳統(tǒng)BEM消耗更多的計(jì)算時(shí)間，但由于結(jié)構(gòu)模型的精確建立也產(chǎn)生更高精度的計(jì)算結(jié)果。另外對(duì)于帶有很多間斷不光滑結(jié)構(gòu)面的描述，基于NURBS的IGA方法還是很難實(shí)現(xiàn)的，因此限制了該方法對(duì)于實(shí)際復(fù)雜問題的應(yīng)用。但一些方法可以用來改善這一缺點(diǎn)，例如T樣條插值、細(xì)分曲面等。雖然IGA有這些難以克服的缺點(diǎn)，但由于其在結(jié)構(gòu)形狀表征和形狀優(yōu)化上的獨(dú)特優(yōu)勢，國內(nèi)外很多學(xué)者致力于該方法的開發(fā)和應(yīng)用研究上。

本文采用IGABEM進(jìn)行聲學(xué)敏感度分析，使用Burton-Miller方法克服外聲場的非唯一解的問題[13-14]，并推導(dǎo)出無奇異敏感度邊界積分方程。然后以等幾何控制點(diǎn)坐標(biāo)作為設(shè)計(jì)變量，以聲屏障降噪?yún)^(qū)域的聲壓作為目標(biāo)函數(shù)，以結(jié)構(gòu)的面積大小作為約束，建立二維聲屏障結(jié)構(gòu)聲學(xué)優(yōu)化數(shù)學(xué)模型，并采用MMA(Moving Approximation Algorithm)方法進(jìn)行尋優(yōu)計(jì)算。

1 基于邊界元法的二維聲學(xué)狀態(tài)分析

1.1 基于傳統(tǒng)邊界元法的二維聲學(xué)狀態(tài)分析

對(duì)于半空間二維聲學(xué)問題，傳統(tǒng)邊界積分方程(Conventional Boundary Integral Equation，CBIE)和法向?qū)?shù)邊界積分方程(Normal Derivative Boundary Integral Equation，NDBIE)分別表示為

(1)

(2)

式中：r0為聲源點(diǎn)；β為聲導(dǎo)納。格林函數(shù)為

(3)

將CBIE方程和NDBIE方程線性組合后形成的Burton-Miller方程可以有效的克服數(shù)值解的非唯一性問題。將邊界S離散后，可得到如下的線性代數(shù)方程組

[H-ikGB]P=Pi

(4)

(5)

式中：H和G為邊界元系數(shù)矩陣；P為節(jié)點(diǎn)聲壓向量；Pi為入射聲壓向量；βi為結(jié)構(gòu)表面第i個(gè)邊界單元上的吸聲材料介質(zhì)導(dǎo)納參數(shù)。

實(shí)際上，式(1)和式(2)的積分計(jì)算需要克服奇異性問題。Sx代表包含源點(diǎn)r的邊界單元，在該邊界單元上的積分具有奇異性，需要做出特別處理。在這里，柯西主值積分和Hadamard有限部分積分法被用于克服這一問題。 在Sx上的兩項(xiàng)奇異積分分別為

(6)

(7)

式中：ξ為每個(gè)單元的局部坐標(biāo)；J為雅可比。

(8)

(9)

式中：m為邊界元內(nèi)插值節(jié)點(diǎn)的個(gè)數(shù)； 當(dāng)拉格朗日函數(shù)插值用于數(shù)值解時(shí)，pi為單元節(jié)點(diǎn)聲壓向量。而對(duì)于NURBS基函數(shù)，則表示在對(duì)應(yīng)的控制點(diǎn)上的聲壓值。將式(8)代入式(9)，我們可以得到

(10)

使用不同類型的拉格朗日函數(shù)進(jìn)行邊界元積分計(jì)算，可得到Bi和Di的不同表達(dá)。當(dāng)使用NURBS基函數(shù)進(jìn)行邊界元插值計(jì)算時(shí)，需作出不同處理。首先，NURBS基函數(shù)為

(11)

(12)

(13)

(14)

然而， 對(duì)于p=1,2,3。

(15)

將式(11)代入式(10)， 可得到基于NURBS插值計(jì)算的Bi和Di的表達(dá)式為

(16)

(17)

在精確得到奇異積分的計(jì)算表達(dá)式后，亦是給出系數(shù)Bi和Di的表達(dá)式，進(jìn)而可以得到無奇異的邊界積分方程式(6)和式(7)。最后使用式(4)可得到節(jié)點(diǎn)處未知變量值。

1.2 基于傳統(tǒng)邊界元法的二維聲學(xué)形狀靈敏度分析

將式(1)和式(2)分別對(duì)設(shè)計(jì)變量進(jìn)行求導(dǎo)，可得

(18)

(19)

(20)

(21)

式(20)和式(21)等號(hào)右邊第二項(xiàng)是非奇異的，可以用高斯積分法求解。然而，右邊的第一項(xiàng)是奇異性的。通過使用式(9)，我們可以得到以下表達(dá)式

(22)

(23)

當(dāng)使用NURBS基函數(shù)進(jìn)行數(shù)值計(jì)算時(shí)，邊界曲線可表示為

(24)

式中：X(ξ)為邊界點(diǎn)的坐標(biāo)；n為控制點(diǎn)數(shù)量；Pi為控制點(diǎn)的坐標(biāo)。將式(24)對(duì)設(shè)計(jì)變量進(jìn)行求導(dǎo)，可得如下邊界曲線敏感度表達(dá)式

(25)

同樣地，可得到坐標(biāo)導(dǎo)數(shù)的敏感度表達(dá)式邊界積分解

(26)

結(jié)合式(18)和式(19)，消除奇異邊界積分后，最終可以得到如下敏感度邊界積分方程

(27)

利用上述方程，我們可以得到用于物理場近似的特殊控制點(diǎn)處的聲壓敏感度值。然后利用式(18)得到域內(nèi)任意點(diǎn)處的聲壓敏感度pf。

2 聲學(xué)結(jié)構(gòu)形狀優(yōu)化模型的建立

對(duì)于結(jié)構(gòu)聲學(xué)形狀優(yōu)化問題，等幾何控制點(diǎn)決定結(jié)構(gòu)形狀及幾何尺寸，將其選做優(yōu)化設(shè)計(jì)變量，能靈活的體現(xiàn)形狀的變化。選取參考區(qū)內(nèi)若干點(diǎn)上的聲壓作為目標(biāo)函數(shù)。在二維聲學(xué)結(jié)構(gòu)形狀優(yōu)化中，結(jié)構(gòu)面積是重要的約束條件，以聲影區(qū)參考點(diǎn)聲壓幅值在一定頻帶上的均值為目標(biāo)函數(shù)，滿足多約束條件下的目標(biāo)函數(shù)最小為設(shè)計(jì)目標(biāo)，優(yōu)化過程中采用MMA法更新設(shè)計(jì)變量[15-17]。綜上所述，該優(yōu)化問題的數(shù)學(xué)模型可表示為

(28)

目標(biāo)函數(shù)∏(X)對(duì)形狀設(shè)計(jì)變量x的靈敏度表示為

(29)

式中： R為取復(fù)數(shù)的實(shí)部值。 當(dāng)采用基于梯度分析的MMA方法進(jìn)行優(yōu)化分析時(shí)， 約束函數(shù)A0(封閉結(jié)構(gòu)的面積)的求解是非常重要的，因此結(jié)構(gòu)的面積可由邊界積分計(jì)算出，表達(dá)式為

(30)

另外結(jié)構(gòu)面積對(duì)設(shè)計(jì)變量的敏感度為

(31)

式中：Ne為NURBS單元數(shù)； [ξe,ξe+1]為第e個(gè)NURBS單元在參數(shù)空間內(nèi)的位置。

3 直立型聲屏障的形狀優(yōu)化

3.1 直立型聲屏障模型描述

圖1 直立型聲屏障初始模型Fig.1 Initial model of vertical shape sound barrier

利用NURBS建模的初始控制點(diǎn)為圖中6個(gè)點(diǎn)，坐標(biāo)分別為p0=(5.2,0)，p1=(5.2,3)，p2=(5.2,3)，p3=(5,3)，p4=(5,3)，p5=(5,0)。 對(duì)應(yīng)的節(jié)點(diǎn)向量為Ξ=[0.0,0.0,0.0,0.45,0.45,0.55,0.55,1.0,1.0,1.0]， 階數(shù)pg=2。

3.2 單頻激勵(lì)下優(yōu)化結(jié)果

本節(jié)研究了聲屏障在單頻激勵(lì)下的形狀優(yōu)化。需要注意的是，選擇適當(dāng)?shù)目刂泣c(diǎn)數(shù)量對(duì)于基于IGA的形狀優(yōu)化分析是十分必要的。更多的控制點(diǎn)數(shù)可以生成更靈活更復(fù)雜的的結(jié)構(gòu)外形，自然可得到更低目標(biāo)函數(shù)值，但設(shè)計(jì)變量的增加往往大大提高CPU(Central Processing Unit)計(jì)算時(shí)間，因此找到合適的優(yōu)化控制點(diǎn)數(shù)是很有必要的。在圖2中，縱坐標(biāo)表示目標(biāo)函數(shù)值，橫坐標(biāo)表示MMA迭代次數(shù)，4個(gè)不同控制點(diǎn)數(shù)被用于優(yōu)化分析。起初隨著MMA迭代次數(shù)的增加，不同設(shè)計(jì)變量數(shù)目的曲線都會(huì)快速收斂到一個(gè)固定值，同時(shí)隨著控制點(diǎn)數(shù)N由5～20的增加，目標(biāo)函數(shù)最終收斂值也在逐一下降表明隨著控制點(diǎn)數(shù)的增加，優(yōu)化自由度在不斷提高，可得到更好的降噪效果。與此同時(shí)，需要更多的CPU時(shí)間，并且此時(shí)目標(biāo)函數(shù)隨N增加時(shí)下降速率逐漸變緩，表明隨著N的增加，每增加一個(gè)控制點(diǎn)目標(biāo)函數(shù)值降低的幅度是越來越小的。這說明當(dāng)優(yōu)化自由度增加到一定程度后，再增加設(shè)計(jì)變量數(shù)目降噪效果是不會(huì)有太大提高的。圖3顯示了不同優(yōu)化控制點(diǎn)數(shù)目下面積約束收斂的效果，設(shè)計(jì)變量數(shù)目越多，最終面積的利用率就越大。圖4表示計(jì)算頻率400 Hz時(shí)優(yōu)化控制點(diǎn)數(shù)目越多，優(yōu)化后聲屏障的形狀就越復(fù)雜，降噪效果也越顯著。從圖5中可看出，當(dāng)控制點(diǎn)數(shù)目相同計(jì)算頻率不同時(shí)優(yōu)化結(jié)果并不相同，計(jì)算頻率越大優(yōu)化后的結(jié)構(gòu)形狀越復(fù)雜，說明單頻優(yōu)化結(jié)果具有很強(qiáng)的頻率依賴性，因此需要進(jìn)一步在整個(gè)頻帶上進(jìn)行形狀優(yōu)化研究。在圖6中，橫坐標(biāo)表示目標(biāo)區(qū)域內(nèi)參考點(diǎn)的編號(hào)，縱坐標(biāo)表示聲壓級(jí)值。Initial value表示聲屏障初始形狀時(shí)目標(biāo)區(qū)域內(nèi)參考點(diǎn)的聲壓級(jí)值，Optimal value表示聲屏障形狀優(yōu)化后參考點(diǎn)的聲壓級(jí)值。聲壓級(jí)值與聲壓值的關(guān)系為

圖2 400 Hz時(shí)不同優(yōu)化控制點(diǎn)數(shù)目下目標(biāo)函數(shù)隨迭代次數(shù)變化Fig.2 Objective function in terms of iteration step for different numbers of design variables

圖3 400 Hz時(shí)不同優(yōu)化控制點(diǎn)數(shù)目下面積約束隨迭代次數(shù)變化Fig.3 Area constraint in terms of iteration step for different numbers of design variables

(32)

式中：φ0為參考聲壓，空氣中參考聲壓取2×10-5Pa。從圖6中可以看出聲屏障形狀在優(yōu)化后聲壓級(jí)值明顯降低，表明降噪效果有了很大的改善。與此同時(shí)優(yōu)化結(jié)果對(duì)頻率的依賴性依然較強(qiáng)。通過對(duì)不同頻率下聲屏障形狀的優(yōu)化，聲壓級(jí)值降低的幅度也是不同的。如頻率400 Hz時(shí)降噪改善在5 dB，頻率700 Hz時(shí)降噪改善達(dá)到8 dB等。從聲源波長的角度分析，聲源頻率越高，波長就越短，進(jìn)而波長與聲屏障的零碎尺寸比值就越小，從而越不容易發(fā)生衍射，表明聲屏障對(duì)高頻聲源的降噪效果更具敏感性。

圖4 400 Hz 時(shí)不同優(yōu)化控制點(diǎn)數(shù)目下聲屏障最終優(yōu)化形狀Fig.4 The optimized shape of the sound barrier with different numbers of optimization control points

圖5 不同頻率下聲屏障最終優(yōu)化形狀Fig.5 The final shape of the sound barrier with 15 control points at different frequencies

圖6 優(yōu)化前后不同頻率下參考區(qū)內(nèi)計(jì)算點(diǎn)處聲壓級(jí)變化Fig.6 Improvement of noise reduction at different frequencies

總之，基于IGABEM方法的形狀優(yōu)化分析，選取一定數(shù)量的形狀控制點(diǎn)為設(shè)計(jì)變量，數(shù)量越多，優(yōu)化后的形狀就越復(fù)雜，雖然可得到更優(yōu)的計(jì)算結(jié)果，但由于設(shè)計(jì)變量數(shù)量的增加，也消耗了更多的計(jì)算時(shí)間，因此，選取一個(gè)合適的設(shè)計(jì)變量數(shù)是重要的。經(jīng)驗(yàn)上，在計(jì)算量可以接受的情況下，盡可能在邊界線或面上布置更多的優(yōu)化控制點(diǎn)。下一步，我們會(huì)進(jìn)一步考察優(yōu)化控制點(diǎn)數(shù)與計(jì)算量之間的對(duì)應(yīng)關(guān)系，給出經(jīng)驗(yàn)表達(dá)式。另外，優(yōu)化控制點(diǎn)數(shù)的增加，還導(dǎo)致優(yōu)化過程中的收斂迭代步數(shù)大大增加，大幅降低優(yōu)化算法的計(jì)算效率，下一步，通過發(fā)展有效的預(yù)處理方法和迭代求解算法，提高優(yōu)化算法的計(jì)算效率。

3.3 頻帶激勵(lì)下優(yōu)化結(jié)果

由于現(xiàn)實(shí)生活中聲源的頻率是復(fù)雜的，并不是以單頻存在，通常分布在某一個(gè)特定的頻率帶內(nèi)。為此需要在一個(gè)頻率范圍內(nèi)建立一個(gè)平均目標(biāo)函數(shù)，頻帶內(nèi)的目標(biāo)函數(shù)為

(33)

式中：ω1為下限頻率；ω2為上限頻率； ∏為單頻下的目標(biāo)函數(shù)，可參見式(28)。頻帶內(nèi)目標(biāo)函數(shù)敏感度表示為

(34)

圖7 200～500 Hz頻帶內(nèi)聲屏障最終優(yōu)化形狀Fig.7 The final shape of the sound barrier in certain frequency band

另外，得到的優(yōu)化聲屏障結(jié)構(gòu)外形由于表面凸凹不定，復(fù)雜多變，是難以直接應(yīng)用于實(shí)際工程問題，實(shí)際上可通過增加約束條件達(dá)到符合實(shí)際規(guī)范的優(yōu)化結(jié)構(gòu)外形，這也是我們接下來的重點(diǎn)工作。本文的重點(diǎn)是進(jìn)行基于等幾何的結(jié)構(gòu)聲學(xué)優(yōu)化設(shè)計(jì)算法研究，使用聲屏障結(jié)構(gòu)主要是為了驗(yàn)證開發(fā)算法的有效性。

圖8 不同頻帶下聲屏障優(yōu)化后降噪效果Fig.8 Noise reduction effect after noise barrier optimization in terms of frequencies

4 結(jié) 論

本文主要介紹了二維聲場等幾何邊界元的聲學(xué)狀態(tài)分析、聲學(xué)形狀敏感度分析及建立結(jié)構(gòu)形狀優(yōu)化模型，導(dǎo)出了帶有設(shè)計(jì)變量的目標(biāo)函數(shù)靈敏度分析表達(dá)式，并通過實(shí)例驗(yàn)證了基于等幾何邊界元法的聲學(xué)結(jié)構(gòu)形狀優(yōu)化算法，結(jié)果表明優(yōu)化后的結(jié)構(gòu)形狀對(duì)降噪效果確有改善，同時(shí)聲學(xué)結(jié)構(gòu)的形狀優(yōu)化與頻率是密切相關(guān)的，在不同的計(jì)算頻率下最優(yōu)解不同，即降噪效果不同。如與實(shí)際相結(jié)合，通?？疾煲欢l帶內(nèi)的優(yōu)化結(jié)果往往更加有效。該方法的突出特點(diǎn)是將設(shè)計(jì)變量作為等幾何建模的控制點(diǎn)，采用MMA方法更新優(yōu)化過程中的設(shè)計(jì)變量，靈活準(zhǔn)確地控制形狀變化，并通過敏感度分析找到形狀變化與聲場聲壓分布的對(duì)應(yīng)關(guān)系，從而解決最優(yōu)解問題。