基于最小化振速重構誤差的等效源分布優化方法研究

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

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

0 引 言

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

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

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

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

1 波疊加法理論背景

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

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

通常將分布聲源置于結構內部的一虛擬球殼面上，經離散化處理及代數運算，式（1）可表示為[7]

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

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

2 等效源分布優化

2.1 振速重構誤差

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

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

質量守恒方程的線性形式和狀態方程為

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

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

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

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

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

2.2 等效源分布優化確定

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

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

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

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

2.3 頻率閾值準則

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

結構向外輻射的聲壓可表示為[7]

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

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

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

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

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

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

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

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

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

3 數值仿真

3.1 等效源分布優化方法驗證

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

圖2 殼型結構幾何外形及尺寸Fig.2 Shape and size of the shell structure

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

其輻射的聲壓解析解為

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

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

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

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

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

圖4 殼型結構邊界節點處預報誤差隨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變化而改變。結果表明：εv、εreal和εimag具有相同的變化規律，εv很好地反映了聲壓預報誤差的變化趨勢；聲壓預報精度隨εv的減小而提高，在εv取最小值的等效源分布處預報的聲場精度較其他位置處高。對比圖3和圖4可發現，隨著結構外形趨于復雜，聲場預報精度對等效源位置分布變得更為敏感。對擺動球源模型，為使聲壓預報誤差小于1%且保證方程求解過程穩定，Sc需滿足0

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

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

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

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

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

圖7 殼型結構聲壓幅值對比Fig.7 Comparison of sound pressure amplitudes for the shell structure

3.2 頻率閾值準則驗證

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

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

表4 殼型結構在不同計算頻率下優化確定的Sc值Tab.4 Optimal values of Sc corresponding to different frequencies for shell structure

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

4 試驗驗證

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

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

表5 鋼制圓柱殼結構物理參數Tab.5 Physical data of the steel cylindrical shell

輸入揚聲器的正弦信號的初始電壓和電流為15 V 和1.6 A，頻段為100～700 Hz，步長為Δf=10 Hz，在改變頻率時保持電壓恒定。逐點測量單元節點和中心處的法向加速度，根據加速度與速度的轉換關系來獲取這些測量點處的速度，其中使用的加速度傳感器與圓柱殼的質量比為0.14%，從而可忽略傳感器質量對結構振動的影響。測量如圖10所示的場點聲壓（聲傳感器已采用B&K 4231 聲學校準儀校準），現場聲壓采集如圖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給出等效源位于優化位置時聲壓預報誤差隨頻率的變化，表6所示為聲壓預報誤差統計。

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

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

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

5 結 語

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