基于最小化振速重構(gòu)誤差的等效源分布優(yōu)化方法研究

吳紹維，肖程詩，王紅梅，王 艷

（1.重慶交通大學航運與船舶工程學院，重慶 400074；2.船舶動力工程技術(shù)交通運輸行業(yè)重點實驗室，武漢 430063）

0 引 言

邊界元法（Boundary Element Method，簡稱BEM）是預報自由場結(jié)構(gòu)聲輻射最常用的方法，使用BEM 只需離散結(jié)構(gòu)邊界表面，且滿足無窮遠處的Sommerfeld 輻射條件[1]。但在常規(guī)的BEM 中，在與內(nèi)部問題對應(yīng)的特征頻率處將出現(xiàn)非唯一性問題[2]。為克服此問題，Combined Helmholtz Integral Equa?tion Formulation（CHIEF）方法和Burton-Miller 方法是最常用的方法[1,3]。在CHIEF 方法中，當選取的CHIEF 點與結(jié)構(gòu)邊界為界的內(nèi)部Dirichlet 問題的振型節(jié)點重合時，該方法失效[3]。Burton-Miller 法可在全頻段解決此問題，但引入了格林函數(shù)各類奇異積分問題[3]，這些奇異積分需進行復雜的數(shù)學處理[4-5]，導致計算效率顯著下降。使用有限元法（Finite Element Method，F(xiàn)EM）預報自由場結(jié)構(gòu)聲輻射面臨無限域截斷的問題[1,6]，需采用人工邊界將無限域截斷為有限計算域，并在人工邊界處施加無反射邊界條件來代替無窮遠處的Sommerfeld 輻射條件，計算域越大，消除反射聲波的效果越好，但計算耗時增加，因此人工邊界的選取需平衡聲場預報精度和計算效率[6]。

為克服上述方法中的困難，Koopmann 等[7]提出了波疊加法（Wave Superposition Method，WSM），結(jié)構(gòu)輻射的聲場使用一組位于結(jié)構(gòu)內(nèi)部源面上的等效源輻射的聲場線性疊加等效，該方法滿足無窮遠處的Sommerfeld 輻射條件，利用邊界法向振速可求解出等效源的源強，然后通過矩陣運算即可預報任意場點處的聲場量。相比于BEM 和FEM，該方法避免了奇異性問題，且無需對計算域進行截斷，簡化了聲場預報過程。夏雪寶等[8]基于波疊加法研究了聲輻射阻抗計算，分析了結(jié)構(gòu)聲輻射矩陣特性。錢治文等[9-10]基于波疊加法提出了一種預報淺海信道下彈性結(jié)構(gòu)聲輻射的聯(lián)合波疊加法。但WSM 也存在不足，其預報精度受等效源分布及計算頻率的影響。等效源距離結(jié)構(gòu)邊界較遠，系統(tǒng)矩陣呈病態(tài)，使得方程求解不穩(wěn)定而導致預報結(jié)果出現(xiàn)震蕩；過于接近邊界面時，格林函數(shù)產(chǎn)生近奇異性，導致錯誤的聲場預報值。

針對等效源分布問題，國內(nèi)外學者開展了眾多研究。Fahnline[11]基于奇異值分解研究了波疊加法的精度和穩(wěn)定性，通過對聲場進行近似奇異函數(shù)展開，從物理意義上揭示了聲場的本質(zhì)，研究表明：波疊加法中的非唯一性是由于奇異值太小所引起，為確保預報精度，等效源應(yīng)均勻分布。Hwang[12]基于虛擬聲源分布提出了一種聲輻射和聲散射計算的正則積分方程，對一定的網(wǎng)格劃分模式，為確保相應(yīng)的積分收斂，等效源與對應(yīng)的表面單元之間的距離應(yīng)大于此單元的1/4特征長度。這從理論上給出了等效源距離邊界面的最小距離，但僅給出下限值。若等效源距離邊界太遠，此時等效源之間相互距離太近，引起有關(guān)矩陣條件數(shù)增大，導致方程求解過程不穩(wěn)定。Gounot[13-15]提出了一種優(yōu)化確定等效源分布的遺傳算法，此方法使用少量等效源，可較精確地重構(gòu)簡單結(jié)構(gòu)輻射的聲場，但該優(yōu)化方法對多數(shù)量等效源失效。Zellers 及向陽等[16-17]分析了系數(shù)矩陣對角項，使用近似解析表達來替代這些對角項，但近似解析表達僅在低頻率時可代替對角項。Xu等[18]提出了一種改進的波疊加法，該方法采用正則化策略，減小了系數(shù)矩陣的條件數(shù)，在一定程度上提高了計算穩(wěn)定性。

針對等效源位置分布問題，本文從理論上推導了聲壓預報誤差與振速重構(gòu)誤差之間的關(guān)系，揭示了兩者之間的變化規(guī)律，提出了一種優(yōu)化確定等效源分布的方法來提高波疊加法的聲場預報精度；對于一定數(shù)量的等效源，提出了一種頻率閾值準則來確定能夠精確預報聲場的頻率范圍；設(shè)計球型聲源和殼型結(jié)構(gòu)數(shù)值仿真對所提方法的精確性和穩(wěn)定性進行了驗證；最后，通過圓柱殼試驗對所提方法的實用性進行驗證。

1 波疊加法理論背景

波疊加法的原理是結(jié)構(gòu)輻射的聲場由一組位于結(jié)構(gòu)內(nèi)部輔助源面上的點聲源疊加等效表達，這些點聲源稱之為等效源。結(jié)構(gòu)向外部場點r處輻射的聲壓表示為[7]

圖1 波疊加法及等效源分布示意圖Fig.1 Diagram of the wave superposition method and equivalent source array

通常將分布聲源置于結(jié)構(gòu)內(nèi)部的一虛擬球殼面上，經(jīng)離散化處理及代數(shù)運算，式（1）可表示為[7]

式中，?表示關(guān)于邊界面上場點rs的梯度算子，ns為場點rs處的單位外法向，v(rs)為rs處的法向速度。將式（2）代入式（3）并將其運用于邊界上的N個節(jié)點，即可求解出等效源的源強，并將其代入式（2），則結(jié)構(gòu)外部任意場點處的聲壓為

式中，v(rsj)為第j個邊界節(jié)點處的法向速度，Dij為

2 等效源分布優(yōu)化

2.1 振速重構(gòu)誤差

由式（3）可知，邊界處的聲壓與振速并非獨立，根據(jù)式（3）和式（4），邊界上任意點rs處的法向振速可表示為

式中，Nj(rs)為速度插值函數(shù)。當rs=rsh（rsh為第h個邊界節(jié)點坐標位置），Nj(rsh)如下式所示：

質(zhì)量守恒方程的線性形式和狀態(tài)方程為

將式（8）代入式（10）得到：

由上式可知，εp(rs)滿足Helmholtz方程。式（8）可改寫為

從而εp(rs)與εv(rs)滿足自由場聲輻射形式的結(jié)構(gòu)邊界條件。由于p(r)和p?(r)均滿足Sommerfeld 輻射條件，從而在無窮遠處εp(rs)滿足：

式中，r表示柱或球坐標中的極徑，α為空間維數(shù)常數(shù)。式（13）表明，在無窮遠處，εp(r)滿足Sommerfeld形式的輻射條件。

式（11）～（13）與聲壓在無限域自由場滿足的控制方程和邊界條件形式一致，類比結(jié)構(gòu)振動向外輻射聲壓，邊界上的振速重構(gòu)誤差產(chǎn)生聲壓預報誤差。抑制邊界振速可降低向外輻射的聲壓，類比可知，減小εv(rs)能降低εp(rs)，從而εv(rs)可用來反映聲壓預報的精度，通過最小化振速重構(gòu)誤差可優(yōu)化確定等效源位置。在WSM 中，采用單極子型等效源在特征頻率處會產(chǎn)生非唯一性問題，此特征頻率與等效源的位置有關(guān)。當聲壓出現(xiàn)非唯一性時，振速重構(gòu)誤差急劇增加。因此，對于給定的頻率，在確定的等效源優(yōu)化位置處可避免非唯一性問題。

2.2 等效源分布優(yōu)化確定

為精確預報聲壓，偶極矩陣D應(yīng)是對角占優(yōu)的。這就需對等效源和邊界節(jié)點配對，使得每個等效源與所配對的節(jié)點之間的距離小于此等效源與其他節(jié)點的距離。但是，只有球形類結(jié)構(gòu)以及無限長圓柱類結(jié)構(gòu)滿足此條件。為減少確定等效源位置分布的優(yōu)化參數(shù)數(shù)量，同時使等效源源面與結(jié)構(gòu)邊界共形，采用縮比系數(shù)Sc（0

采用兩種離散模式A 和B 離散邊界，由于εv(rs)為矢量，為量化振速重構(gòu)誤差，定義一種無量綱振速重構(gòu)誤差，其解析和數(shù)值形式如下：

式中，Sj為離散模式B 下第j個節(jié)點所占區(qū)域。采用縮比系數(shù)Sc乘離散模式A 下的節(jié)點來獲取等效源的空間坐標，確保等效源源面與結(jié)構(gòu)表面共形，同時使用離散模式A下的節(jié)點法向振速求解等效源源強和計算聲壓。為優(yōu)化確定Sc值，設(shè)計目標函數(shù)為

計算步驟為：對計算頻率k，以離散模式A和B下的節(jié)點振速為輸入，Sc以選定的步長ΔSc遍歷（0，1），對每個Sc值，將離散模式A 下的節(jié)點給定法向振速代入式（6）計算離散模式B 下的節(jié)點重構(gòu)法向振速v?(rsj)，將重構(gòu)法向振速與離散模式B 下的節(jié)點給定法向振速代入式（14）計算無量綱振速重構(gòu)誤差εv；Sc遍歷（0，1）之后，搜索最小εv，其對應(yīng)的Sc值為最優(yōu)Sc值，即可確定等效源的最優(yōu)位置。

2.3 頻率閾值準則

聲波波長隨頻率的升高而減小，對一定數(shù)量的等效源，能精確預報聲場的頻率閾值是有限的，超過此頻率閾值，預報精度急劇下降，此閾值還隨輔助源面的位置改變而變化。因此，在確定等效源分布后需確定頻率閾值，以確保聲場預報精度可接受。

結(jié)構(gòu)向外輻射的聲壓可表示為[7]

式中，σ表示優(yōu)化確定的輔助源面，δτ為厚度常數(shù)。采用四節(jié)點或三節(jié)點單元或兩種單元的任意組合離散σ，從而σ上的積分可表示為這組單元上的積分之和?？臻g坐標和聲場變量可采用單元形函數(shù)插值近似表示，采用合適的高斯積分方法可計算單元上的積分。運用等參單元變換，全局坐標rσ和聲源分布強度q(rσ)可表示為

式中，Nh為已知單元形函數(shù)，rh為單元節(jié)點全局坐標，qh為聲源分布強度在節(jié)點上的值，M為單元節(jié)點數(shù)。從而式（16）可表示為

式中，K為單元數(shù)，Hih(r)可表示為

在局部坐標系s-t中，對四節(jié)點或三節(jié)點單元，式（19）分別等效表示為

式中，J(s,t)為單元映射Jacobian行列式。

對四節(jié)點單元，高斯求積的誤差上限值為[12]

式中，m和n分別為s和v方向的高斯積分點數(shù)，wl和wj分別為與積分點sl和tj對應(yīng)的加權(quán)因子，E1和E2為估計誤差。在式（19）中的被積函數(shù)中，最主要的變化項可歸納為e-ikr/rp這種形式，其中p=1或2。因此，式（19）可近似表示為

由于被積函數(shù)中的Jacobian行列式和場點與單元之間的距離幾乎是不變的，E1和E2可表示為[12]

式中，L1= max(L12,L34)，L2= max(L23,L41)，rmin為場點與單元之間的最小距離。為確保高斯積分收斂，k< min(4/L1,4/L2)。對三節(jié)點單元，同理可得k< min(4/L12,4/L23,4/L31)。令Lmax表示輔助源面上所有單元最大特征長度，為確保聲場預報精度，k應(yīng)滿足

3 數(shù)值仿真

3.1 等效源分布優(yōu)化方法驗證

采用擺動球源（一階球源，半徑a=1，以幅值v=1 m/s 沿z軸橫向振動，坐標原點位于球心）和殼型結(jié)構(gòu)（如圖2 所示）進行數(shù)值仿真，離散模式A 和B 下的單元類型及數(shù)量如表1 所示，其中，空氣的密度ρ=1.21 kg/m3，聲速c=343 m/s。

圖2 殼型結(jié)構(gòu)幾何外形及尺寸Fig.2 Shape and size of the shell structure

表1 離散模式A與B單元類型及數(shù)量Tab.1 Number of elements on the boundary for Modes A and B

其輻射的聲壓解析解為

式中，θ為球面法線方向與z軸的夾角。采用如下范數(shù)作為誤差指標：

式中，ε1m(%)和ε2m(%)分別為在場節(jié)點xm處預報的聲壓實部和虛部相對誤差。

殼型結(jié)構(gòu)無聲壓解析解，但可獲取某些振速分布條件下輻射的聲壓準確值。在場空間中用虛擬結(jié)構(gòu)封閉包裹一點聲源，此點聲源使在虛擬結(jié)構(gòu)邊界處的聲介質(zhì)振動，根據(jù)式（3）可計算出聲介質(zhì)沿邊界外法向方向的振速，若結(jié)構(gòu)按此振速分布向外輻射聲場，則結(jié)構(gòu)輻射的聲場應(yīng)與點聲源的聲場一致。因此，該點聲源輻射的聲場可作為對比標準，此點聲源稱為模擬點聲源。

為研究等效源分布對聲場預報精度的影響規(guī)律，圖3 和圖4 分別給出兩種模型的εv、邊界節(jié)點上的εreal和εimag隨k和Sc的變化，其中殼型結(jié)構(gòu)仿真采用單極子點聲源作為模擬點聲源，其空間坐標為(0.15 m，0 m，0 m)。

圖3 擺動球源邊界節(jié)點處預報誤差隨k與Sc的變化Fig.3 Prediction error on boundary surface of the transversely-oscillating sphere for different k and Sc val?

圖4 殼型結(jié)構(gòu)邊界節(jié)點處預報誤差隨k與Sc的變化Fig.4 Prediction error of εv,εreal and εimag on boundary surface of the shell structure for different k and Sc values

由圖3 和圖4 可知，等效源分布對聲場預報精度影響顯著，聲場預報誤差隨Sc變化而改變。結(jié)果表明：εv、εreal和εimag具有相同的變化規(guī)律，εv很好地反映了聲壓預報誤差的變化趨勢；聲壓預報精度隨εv的減小而提高，在εv取最小值的等效源分布處預報的聲場精度較其他位置處高。對比圖3和圖4可發(fā)現(xiàn)，隨著結(jié)構(gòu)外形趨于復雜，聲場預報精度對等效源位置分布變得更為敏感。對擺動球源模型，為使聲壓預報誤差小于1%且保證方程求解過程穩(wěn)定，Sc需滿足0

圖5 擺動球源在場節(jié)點處輻射聲壓隨k的變化Fig.5 Variation of pressure in the selected field point with k for uniformly-pulsating sphere

圖6 殼型結(jié)構(gòu)在場節(jié)點處輻射聲壓隨k的變化Fig.6 Variation of pressure in the selected field point with k for the shell structure

圖中曲線顯示，等效源位于優(yōu)化分布位置預報的聲壓與準確值吻合良好，聲壓預報精度顯著提高。結(jié)果也再次表明：結(jié)構(gòu)外形越復雜，等效源分布對聲場預報精度的影響越明顯。因此，必須優(yōu)化確定等效源分布來減小聲場預報誤差。根據(jù)數(shù)據(jù)分析，對擺動球源和殼型結(jié)構(gòu)，等效源位于優(yōu)化位置預報的聲壓最大誤差值分別為3.7 × 10-8%和5.6 × 10-3%，優(yōu)化確定的Sc值見表2。由表2 數(shù)據(jù)可知：對不同的結(jié)構(gòu)，優(yōu)化確定的Sc值不同；即使對同一結(jié)構(gòu)，Sc值隨頻率的改變也存在一定的波動。

表2 不同計算頻率下優(yōu)化確定的Sc值Tab.2 Optimal values of Sc determined by the proposed method for different frequencies

為進一步驗證所提等效源分布優(yōu)化方法，對殼型結(jié)構(gòu)分別采用單極子、偶極子和三極子模擬點聲源進行仿真（偶極子與三極子表達式見式（29），其中→n表示偶極方向單位向量）。圖7為等效源位于不同位置預報的y-z平面第一象限、半徑r=0.2 m 場點處的聲壓幅值，其中k=20，圖中pamp表示聲壓幅值。由結(jié)果可知，相比于其他位置，等效源處于優(yōu)化位置預報的數(shù)值解與準確值吻合良好，誤差得以有效抑制。

圖7 殼型結(jié)構(gòu)聲壓幅值對比Fig.7 Comparison of sound pressure amplitudes for the shell structure

3.2 頻率閾值準則驗證

對一定數(shù)量的等效源，等效源位于優(yōu)化分布位置處預報的聲場只在一定頻率范圍內(nèi)具有較好的精度，超過對應(yīng)的頻率閾值kthreshold，誤差將急劇增加。為驗證所提頻率閾值準則，圖8為等效源位于優(yōu)化位置時結(jié)構(gòu)邊界場點處的εreal和εimag隨k的變化圖，對擺動球源和殼型結(jié)構(gòu)兩種模型優(yōu)化確定的Sc值被列于表3 和表4 中，其中殼型結(jié)構(gòu)采用單極子型模擬點聲源。由圖中曲線可知，當計算頻率超過一定值之后，聲場預報精度急劇下降。對擺動球源，kthreshold值以及k

表3 擺動球源在不同計算頻率下優(yōu)化確定的Sc值Tab.3 Optimal values of Sc corresponding to different frequencies for uniformly-pulsating sphere

表4 殼型結(jié)構(gòu)在不同計算頻率下優(yōu)化確定的Sc值Tab.4 Optimal values of Sc corresponding to different frequencies for shell structure

圖8 等效源位于優(yōu)化位置時邊界場點處的εreall和εimag隨k的變化Fig.8 εreal and εimag on boundary field nodes calculated in optimal auxiliary surfaces for different k

4 試驗驗證

為驗證提出的方法，設(shè)計了鋼制圓柱殼自由場聲輻射試驗（見圖9），該試驗是通過測量獲取的結(jié)構(gòu)表面振速來計算結(jié)構(gòu)輻射的聲壓，與試驗測量的聲壓值進行對比。試驗采用揚聲器作為聲激勵源，其在圓柱殼內(nèi)部產(chǎn)生的聲場作用于圓柱殼內(nèi)表面，引起圓柱殼振動向外輻射聲場，結(jié)構(gòu)與揚聲器見圖9（a）和（b），結(jié)構(gòu)尺寸及材料物理屬性見表5，關(guān)于試驗結(jié)構(gòu)的安裝、固定以及密封詳見文獻[19-20]。將圓柱殼外表面劃分為30個均勻分布的四節(jié)點單元，使用單元節(jié)點乘縮比系數(shù)Sc來獲取等效源坐標，節(jié)點上的法向振速用來求解聲源強度和計算聲壓。計算步驟為：對計算頻率f，以單元節(jié)點和單元中心處的實測振速為輸入，Sc以選定的步長ΔSc遍歷（0，1），對每個Sc值，將單元節(jié)點處的測量法向振速代入式（6）計算單元中心處的節(jié)點重構(gòu)法向振速v?(rsj)，將重構(gòu)法向振速與單元節(jié)點處的測量法向振速代入式（14）計算無量綱振速重構(gòu)誤差εv；Sc遍歷（0，1）之后，搜索最小εv并確定其對應(yīng)的Sc值，Sc值對應(yīng)的等效源位置即為最優(yōu)位置。

圖9 鋼制圓柱殼自由場聲輻射試驗圖Fig.9 Test diagram of free field acoustic radiation of steel cylindrical shell

表5 鋼制圓柱殼結(jié)構(gòu)物理參數(shù)Tab.5 Physical data of the steel cylindrical shell

輸入揚聲器的正弦信號的初始電壓和電流為15 V 和1.6 A，頻段為100～700 Hz，步長為Δf=10 Hz，在改變頻率時保持電壓恒定。逐點測量單元節(jié)點和中心處的法向加速度，根據(jù)加速度與速度的轉(zhuǎn)換關(guān)系來獲取這些測量點處的速度，其中使用的加速度傳感器與圓柱殼的質(zhì)量比為0.14%，從而可忽略傳感器質(zhì)量對結(jié)構(gòu)振動的影響。測量如圖10所示的場點聲壓（聲傳感器已采用B&K 4231 聲學校準儀校準），現(xiàn)場聲壓采集如圖9（c）所示，試驗所使用的儀器詳見文獻[19-20]。

圖10 目標場點坐標Fig.10 Coordinates of target field points

圖11給出目標場點處的聲壓計算值和測量值的聲壓級曲線，圖中pamp表示聲壓幅值。采用式（30）計算聲壓級誤差：

圖11 鋼制圓柱殼聲壓幅值測量值與計算值Fig.11 Amplitude of the acoustic pressure of the cylindrical shell at selected field points

式中，pnumericalamp和pmeasuredamp分別表示聲壓幅值計算值和測量值，pref表示參考值。圖12給出等效源位于優(yōu)化位置時聲壓預報誤差隨頻率的變化，表6所示為聲壓預報誤差統(tǒng)計。

表6 等效源位于優(yōu)化分布位置時的聲壓預報誤差Tab.6 Prediction error in optimal auxiliary surface

圖12 等效源位于優(yōu)化分布位置時預報誤差隨頻率的變化Fig.12 Variation of prediction error in optimal auxiliary surface with k

圖11中的曲線表明：在采用所提方法確定的等效源優(yōu)化分布位置處預報的聲壓幅值數(shù)值解與測量值吻合較好，對比Sc=0.75 位置處，具有更高的聲場預報精度。由表6 可知，最大誤差為4.79 dB，平均誤差控制在2.88 dB 范圍內(nèi)。結(jié)果表明：在少數(shù)頻點處預報的聲壓幅值與測量值存在一定的偏離，但平均聲壓預報誤差仍可接受。導致最大誤差偏大的主要原因是：在測量邊界振速時，將所有加速度傳感器底座固定在結(jié)構(gòu)表面各節(jié)點處，改變了結(jié)構(gòu)表面幾何外形，而預報聲場時未考慮這種改變。除此之外，結(jié)構(gòu)周圍的吸聲材料未能完全吸收聲波，地面仍有小部分未能鋪設(shè)吸聲材料，存在少量的反射聲波，這些也導致聲壓預報值與測量值具有一定的偏離。

5 結(jié) 語

本文推導確定了聲場預報誤差與振速重構(gòu)誤差之間的理論關(guān)系，聲場預報誤差可用法向振速重構(gòu)誤差的變化來反映，減小振速重構(gòu)誤差可提高聲場預報精度，通過推導確定了等效源位于優(yōu)化分布位置能夠精確預報聲場的頻率閾值。仿真結(jié)果表明：聲場預報誤差與法向振速重構(gòu)誤差的變化規(guī)律一致，相比其他固定分布位置，等效源位于優(yōu)化分布位置處時，聲場預報精度顯著提高；頻率閾值準則能夠給出精確預報聲場的上限頻率，確保聲場預報精度。試驗結(jié)果表明，提出的等效源分布優(yōu)化方法能有效降低聲場預報誤差，可用于實際聲場計算。