聲振耦合聲場分析與聲輻射結構優化

王旭東，朱擁勇，王德石

(1.海軍工程大學 兵器工程學院，武漢 430033；2.中國人民解放軍92064部隊，廣東 東莞 523900)

同時具有內部聲場和外部自由聲場的聲振耦合系統廣泛存在于汽車、飛機和潛艇等運載工具和武器平臺之中。結構與聲介質因振動相互作用而產生的中低頻噪聲會對人員及設備產生負面影響，輻射噪聲影響人員舒適性和造成環境污染，在軍事上，聲輻射過大會影響武器系統的隱身性和自身的聲吶探測距離[1]。因此在準確計算聲輻射的基礎上，通過適當設計材料的結構分布以使耦合系統在所關心的激勵頻段下具有較低的聲輻射值的研究十分必要。

作為新興的結構動力優化方法，不同于尺寸優化[2]和形狀優化[3]僅能改變已知設計部件的尺寸，拓撲優化的目標是尋求最優的結構材料分布以滿足性能和結構工藝等多方面要求[4]。國、內外諸多研究者在減振降噪問題上采用拓撲優化法開展了研究。在內聲場優化問題中，Akl等[5]采用有限元法建立內聲場聲輻射模型，采用MMA法優化柔性板厚度，降低了聲腔內的聲輻射，數值計算結果與實驗吻合較好；Dühring等[6]結合混合有限元法和MMA優化法實現在無法明確耦合邊界的情況下解決聲-結構耦合優化問題；隨后Kook[7]用Comsol軟件實現將混合有限元法擴展至BESO優化方法中，優化結果與采用傳統有限元法的文獻[8]結果一致。在外聲場優化問題的研究中，陳爐云等[9]結合ESO優化法和單向耦合的邊界元法研究了簡支板在低頻率下外場聲輻射優化問題；張軍[10]推導了空氣為聲介質的有限元-邊界元耦合方法的聲壓級和輻射聲功率靈敏度公式，并以板厚為設計變量對薄板聲輻射進行優化；針對兩相材料板結構的聲輻射拓撲優化問題，文獻[11-12]通過優化布局兩種彈性材料的方式降低聲輻射；此后，黃其柏在上述基礎上研究了板結構聲輻射計算方法和輻射特性，提出約束阻尼復合板的拓撲優化方法；Zhang等[13]采用狀態空間中的復模態疊加法求解非比例全局阻尼的復合板聲響應，伴隨變量法進行聲壓靈敏度分析，通過合理分配阻尼層材料達到最小化指定點聲壓的目標。

目前，聲振耦合方法的研究大多針對單一的有限內聲場[5-8]或無限大/半無限大自由聲場問題[9-13]，且在處理自由聲場靈敏度分析時往往忽略了聲介質對結構的反作用，無法適用于內聲場和外部自由聲場均存在的聲輻射問題和強耦合問題。基于上述問題，本文結合聲振耦合有限元方程和邊界元法，推導了內外聲場聲壓的耦合計算公式，并提出一種新的強耦合聲功率的靈敏度求解方法，通過不同設計條件下的聲輻射優化設計結果驗證了聲輻射預報的準確性和靈敏度方法在優化過程中的適用性。

1 聲輻射問題

1.1 聲學基本方程

聲振耦合示意圖如圖1所示，假定流體是可壓縮的，無黏且無旋，忽略自由表面重力效應的均勻理想聲學介質，聲壓滿足如下的波動方程：

(1)

?2pf+(ω/c0)2pf=0

(2)

在流固交界面上，應滿足如下的運動學和力學條件：

(3)

usns+ufnf=0

(4)

σs|ns=-pf

(5)

式中：ρf為流體密度，nf為耦合表面從聲場指向結構的法向量，ns為耦合表面從結構指向聲場的法向量，σs為結構應力，us為結構位移向量。

圖1 聲振耦合示意

1.2 內聲場聲振耦合的有限元法

結構域的平衡狀態方程為

(6)

式中：ρs為結構質量密度，fs為結構所受外激勵。

聯立式(2)、(3)～式(6)，并采用加權余量的伽遼金法可得到聲振耦合的有限元方程[14]：

(7)

式中：Ks、Kf分別為結構域和流體域的系統剛度矩陣，Ms、Mf分別為結構域和流體域系統質量矩陣，C為系統耦合矩陣。

為方便后文推導計算，式(7)可簡寫為

(Kfe-ω2Mfe)Xfe=Ffe

(8)

1.3 自由場聲振耦合的邊界元法

穩態諧激勵下，為求解結構表面上受到介質的壓力，需要建立介質中任一點的聲壓積分方程。對聲壓和格林函數應用格林第二公式，由Sommerfeld輻射定律，外域自由聲場通過Helmholtz積分方程求解：

(9)

式中：px為源點聲壓；py為場點聲壓;c為空間角系數，光滑面上取0.5，聲場內取1;G(x,y)為自由空間格林函數;Γbe為自由場聲輻射邊界。

以二維空間為例，對于二維聲場問題，格林函數G(x,y)和其對場點處的法向導數的基本階函數如下：

(10)

(11)

聯立式(3)，式(9)～式(11)，并通過Kronecker函數合并聲壓項得到邊界積分方程：

Hpbe-Gvbe=0

(12)

式中：H、G分別為采用邊界積分方程得到的頻率相關的矩陣，pbe、vbe分別為邊界上的聲壓和法向速度向量。

1.4 聲振耦合方程及輻射聲功率

當聲場由有限大內聲場和無限大自由聲場組成時，式(8)的平衡方程須引入自由場聲介質對結構表面的反作用力項為

(Kfe-ω2Mfe)Xfe=Ffe+Lpbe

(13)

式中L為結構與自由聲場的邊界壓強與邊界法向作用力的耦合矩陣。

(14)

由穩態諧波激勵下的界面位移與界面法向速度的界面連續性條件可知：

(15)

式中：Nu為結構單元的位移形函數，Γe為自由聲場與結構的單元耦合邊界，m為耦合邊界單元總數。上式求積分并求逆移項后，將結構位移從邊界單元擴充至全局響應向量Xfe得到：

(16)

式中：Mbe為邊界質量矩陣，其表達式即為式(15)左端積分項；聯立式(12)、(13)和式(16)，得到修正后含外部自由聲場的聲振耦合方程

(17)

直接求解式(17)即可得到內外聲場聲壓值以及結構邊界位移。同式(14)，采用形函數Np對自由場邊界聲速和聲壓進行采用統一的插值離散并代入上式，得到聲功率W的數值計算公式為[15]

(18)

式中Re為取計算值的實部。

式(18)在聲功率靈敏度求解中不便于解耦變量，為方便后文求解單元靈敏度，需要將聲功率公式修改為用結構位移響應表達的方式。將邊界積分方程式(12)代入上式并由復共軛矩陣乘法規則得到采用邊界振速表示的結構聲功率的表達式：

(19)

聲功率靈敏度是求聲功率對結構單元狀態參數xi求微分，為方便計算，應使式(19)盡量簡單。由實部操作符得到式(19)的簡單表達式：

(20)

其中

將式(16)代入式(20)，得到輻射聲功率關于Xfe的關系式：

(21)

其中

2 單元靈敏度分析與優化算法

拓撲優化問題需要對目標函數相對設計變量的靈敏度進行分析，由輻射聲功率的表達式可知，在外載荷和輻射邊界不發生變化的情況下，輻射聲功率僅與結構表面法向速度vbe有關，即是結構響應Xfe的函數。Kang等[16-17]利用伴隨變量法(AVM)推導出廣義結構響應函數g(Xfe)的靈敏度分析方法，該方法明確指出靈敏度依賴于結構位移響應的幅值，隨后他們團隊采用這種伴隨變量法計算僅考慮單向耦合的聲壓靈敏度[13]，取得了一定的成果。本文將該方法推廣到強耦合的聲功率靈敏度分析中。

2.1 聲功率單元靈敏度的AVM方法

在本文中，首先將AVM方法擴展到重介質的靈敏度分析中。對于依賴于位移響應Xfe的結構響應g(Xfe)，其目標函數可以通過引入兩個拉格朗日乘子的和改寫為

(22)

式中，Kd=Kfe-ω2Mfe。上式中存在與Xfe不統一的變量pbe，聯立式(12)和式(16)可得

Lpbe=BXfe

(23)

(24)

將式(24)兩邊對設計變量xi求導，可得

(25)

令伴隨向量滿足以下方程：

(26)

式中μ1、μ2滿足μ1=(μ2)H。式(25)可以化簡為

(27)

由聲振耦合方程(17)求解的Xfe是復向量，將其分解為實部與虛部之和的形式代入式(21)，得到輻射聲功率的表達式：

(28)

式中( )R、( )I分別為取實部和虛部。式(28)對Xfe求導得到：

(29)

(30)

將式(29)、(30)代入式(26)可得到伴隨向量的表達式：

(31)

求解式(31)即可得到拉格朗日乘子μ1，進而輻射聲功率靈敏度可由式(27)得出

(32)

由上述伴隨變量法的推導過程可以發現，該方法推導過程較為繁瑣，中間變量較多，容易出錯。因此，本文提出更為簡潔高效的直接求導法求解輻射聲功率靈敏度。

2.2 聲功率單元靈敏度的直接求導法

式(21)中的輻射聲功率W直接對單元變量xi求導，并由實部操作符得到聲功率的簡單形式：

(33)

式(13)和式(23)分別對xi求導得到：

(34)

由式(34)得到：

(35)

將式(35)代入式(33)得到輻射聲功率的單元靈敏度公式：

(36)

對比本文的直接求導法和AVM法可以發現，本文方法的推導過程較為簡潔高效，中間變量少。最終的靈敏度公式中的各矩陣都在AVM法的最終公式中出現，不需要計算額外的變量，說明靈敏度計算量也并未增加。

2.3 目標函數與優化算法

體積約束條件下，以輻射聲功率C最小化為目標函數，則結構拓撲優化數學模型為：

(37)

式中：V*為結構有效體積，xi為單元狀態變量(取xmin或1)。

在涉及到特征值問題的優化中，低密度區域容易產生局部模態。在使用BESO優化方法時，Huang等[18]對質量和剛度采用相同的罰函數冪率，該模型對部分優化問題是有效的。但是在計算中發現，采用這種方法后有些優化問題在失效單元區域仍然會出現局部模態，這導致優化過程中結構拓撲突變，迭代不穩定難以收斂。Pedersen[19]提出線性化單元剛度的方法，該方法可以使單元在極低密度時仍然能保持較高的剛度質量比，可有效避免偽本征模，使優化結果更加合理，罰函數模型如下：

(38)

式中：ρ0為結構初始密度，E0為結構初始彈性模量，p為罰因子。

結合BESO方法中單元狀態變量xi僅可取xmin或1的特點，同時為保證罰函數模型的連續性，本文提出以下改進的罰函數模型：

(39)

拓撲優化過程中需要采用基于半徑濾波法和單元靈敏度歷史平均方法等來避免“孤立網格”和加速收斂。具體的耦合優化的計算流程如下：

1)初始化結構參數。定義初始耦合域的物理參數、定義BESO算法的各初始參數；

2)生成數值網格、生成單元矩陣和耦合矩陣并組裝生成有限元聲振耦合方程，生成結構外表面邊界元的各矩陣；

3)根據步驟2)計算各矩陣，求解系統響應。根據式(38)計算單元靈敏度。進行靈敏度過濾和靈敏度歷史平均[20]；

4)判斷失效單元是否與聲學單元相鄰，如果相鄰，則將失效單元轉變為聲學單元，并更新結構域和流體域；

5)根據更新后的各域，識別結構域和流體域的耦合單元和耦合邊界；

6)進行收斂性判斷，如不滿足收斂條件，則重復步驟2)～步驟5)，直到滿足收斂條件。

3 數值算例

3.1 聲輻射數值計算結果驗證

為驗證本文文建立的聲振耦合方程的正確性，設置驗證算例。無限大水中放置一個內徑為1.0 m外徑為1.2 m的鋁質圓環，結構密度為2 700 kg/m3，彈性模量為69 GPa，泊松比為0.3，內聲場和外聲場介質為水，其密度為1 000 kg/m3，聲速為1 450 m/s。圓環內側頂點處施加頻率為10 Hz，幅值為1 N，方向為y軸方向的激勵力。

采用Matlab編程計算圓環內外的耦合聲場，并與無反射邊界法的計算計算結果進行對比。在無反射邊界中，取圓環外側半徑5 m的同心圓為外側計算域，外邊界設置為無反射邊界，材料屬性和其他邊界條件與聲振耦合法相同，兩種方式的聲壓計算結果如圖2所示。對比可知，本文的聲振耦合法與無反射邊界法計算的內、外場聲壓分布基本一致。

圖2 圓環內外聲場聲壓分布

圖3給出圓環外邊界聲壓實部的分布曲線，為方便比較，設置起始點為圓環左側中點(與圖3中弧度為0處相對應)，沿順時針方向。兩種數值方法計算的整體匹配度較好，但在谷值處誤差相對較大，原因在于有限元的無反射邊界對正入射波吸收性較好，但對傾斜入射波有一定的反射作用，因此與實際值有一定誤差。此外可明顯看出無反射邊界法的曲線在對稱節點上的聲壓值略微不對稱，而邊界條件和外載荷均是對稱，其邊界聲壓實部也應該是對稱的。出現這種現象的原因是劃分有限元網格無法做到完全對稱，造成一定的誤差。因此，可以認為采用聲振耦合方程可獲得良好的聲壓預報結果。

圖3 結構表面聲壓實部分布

3.2 算例1

圖4為算例1的設計模型，設計域為一外部半徑3.0 m，內部半徑2.0 m的半圓環，外圈厚度0.1 m部分為非設計域，初始設計域為全設計域的80%(即圖中半徑2.2 m部分)，半圓環兩側的底部固定約束，內部流體的底部邊界為零聲壓邊界。內外聲場的聲介質均為水，聲介質屬性和結構材料屬性同聲輻射數值計算結果。內聲場中點Pin點為一幅值大小為100 Pa的聲壓激勵，激勵頻率為125 Hz。本算例以結構外表面輻射聲功率最小化為目標函數，設置目標體積為全結構域的80%不變。計算域用最大0.03 m尺寸的四邊形四節點等參數單元離散，過濾半徑rmin=0.1 m，罰函數因子p=3，xmin=0.000 1。

圖4 拱形梁設計模型

為驗證本文所提靈敏度方法的設計效果，同條件下采用AVM靈敏度法進行對比。迭代過程的聲功率級和靈敏度計算消耗時間如圖5所示，其中基準聲功率為1×10-12W。兩種方法的目標函數在迭代過程中除了在開始有小幅波動之外，均很快達到平穩狀態，說明兩種靈敏度方法的魯棒性較好，且最終的優化目標函數一致。對比兩種方法在迭代過程中靈敏度部分消耗的時間可以發現，兩種靈敏度求解方法在每一個迭代步中的分析時間沒有明顯差異。從式(32)和式(36)看出，兩個靈敏度計算公式都包含相似的變量，因此計算時間沒有太大差異。

圖5 目標函數和靈敏度計算時間的迭代曲線

結構對外聲輻射情況也體現在外部聲場的聲壓級分布情況上。圖6為優化前、后結構拓撲圖與內部聲場和部分外部聲場的聲壓級響應分布，外部聲場截取10 m×6 m范圍。對比內部聲場，優化前、后僅分布范圍改變，強度未發生明顯變化，但從外部聲場的聲壓級分布可以看出，初始設計的外聲場存在大范圍的高聲壓級區域。優化設計模型中，該區域聲壓級明顯低于初始設計，這也印證了優化后結構的聲輻射得到有效的降低。

圖6 優化前、后結構拓撲和聲壓級分布

3.3 算例2

算例1為鋁在水介質中的聲輻射優化，為研究不同聲介質和結構材料的優化效果，算例2的設計模型同算例1，改變聲介質和結構材料。初始設計域為全設計域，半圓環兩側的底部固定約束，流體下邊界為剛性邊界。圓環內部和外部輻射聲場的聲介質均為空氣，密度為1.21 kg/m3，聲速為343 m/s。結構的材料密度為800 kg/m3，彈性模量為2 GPa，泊松比為0.3。聲壓激勵源位于圓心處Pin點，幅值大小為100 Pa，激勵頻率為205 Hz。以結構外表面輻射聲功率最小化為目標函數，設置目標體積為全結構域的65%。單次迭代的體積縮減率為2%，過濾半徑rmin=0.1 m，罰函數因子p=3，xmin=0.000 1。

為方便對比外聲場情況，截取外部半徑為6 m的圓環區域。圖7為優化前、后同色度范圍的聲壓級分布圖和結構拓撲形態。結構拓撲的邊界清晰，且較為規整，內部聲場和外部聲場的聲壓級得到大幅的降低。這與圖8的迭代曲線相對應，結構拓撲漸變穩定，迭代曲線未出現大幅波動，聲功率級從初始的104.8 dB降低至62.1 dB，降低幅度達到40.7%。從優化前、后內部聲場的分布情況看，由激勵源聲波、邊界反射聲波以及結構振動輻射聲波相互疊加形成環形的低聲壓帶，由于結構拓撲的改變，低聲壓帶分布位置和個數發生改變。為方便觀察內部聲壓帶與外部輻射聲場的關系，圖7還給出聲壓級高對比圖，從圖中圈出部分可以看出，圓環內邊界上的聲壓帶個數與外聲場高(或低)聲壓是對應的。

圖7 優化前、后結構拓撲和聲壓級分布

圖8 目標函數和約束體積的優化迭代曲線

3.4 算例3

圖9為算例3的設計模型，結構域為1.0 m×0.2 m的固支梁，梁頂部0.02 m部分為非設計域，初始設計域即為全設計域，結構密度為2 700 kg/m3，彈性模量為69 GPa，泊松比為0.3，設置目標體積為全結構域的85%。結構域下側聲介質域為1.0 m×0.2 m的矩形域，聲介質為水，梁上側聲介質為空氣，各介質屬性與上文相同。聲壓激勵源為矩形流體域的下邊界，幅值大小為100 Pa，激勵頻率為1 940 Hz。同樣以聲功率最小化為目標函數。設計域離散為200×40個四邊形單元，體積縮減率設置為1%，過濾半徑rmin=0.015 m，罰函數因子p=3，xmin=0.000 1。

圖9 固支梁設計模型

截取固支梁上方2.0 m×1.0 m區域的聲壓級分布來比較優化結果，如圖10所示。可以看出，優化后的外聲場聲壓得到明顯降低。圖中的等值線上數值表明，初始設計的輻射聲場最大聲壓級達到100 dB，且占據固支梁上方的大部分區域，優化構型的最大輻射聲壓級降低至65 dB，降低幅度約35%。圖11顯示了兩種模型在不同激勵頻率下的聲功率級曲線，激勵頻率位于初始設計模型的第2個響應峰值處，優化后的結構輻射聲功率峰值頻率降低，避開激勵頻率所在位置，使其位于響應谷值附近，聲功率級從初始的70.6 dB降低至54.7 dB，降低幅度達22.5%。圖12為優化迭代曲線，與前文一致，迭代曲線平緩，結構拓撲形式未發生突變，說明本文算法的穩定性較好。

圖10 優化前、后結構拓撲和聲壓級分布

圖11 聲功率頻率響應曲線

圖12 目標函數和約束體積的優化迭代曲線

4 結 論

1)建立的含內部聲場的強耦合有限元/邊界元模型計算精度較高，可用于求解耦合系統的聲輻射問題。

2)提出的直接求導法與經典的伴隨變量法得到的優化結果一致，且魯棒性較好。

3)本文的線性化剛度模型和靈敏度分析法適用于不同結構材料與聲介質問題，因此使用本文的優化算法可以得到最優拓撲。