基于分段常數水平集方法的聲振耦合系統拓撲優化

苗曉飛， 趙文暢， 陳海波

(中國科學技術大學 近代力學系 中國科學院材料力學行為與設計重點實驗室，合肥 230026)

近些年來，結構振動和聲學優化問題因其在結構減振降噪方面的重要應用前景而成為國際上的研究熱點之一[1]，減少結構關鍵區域的振動與聲輻射的研究工作也屢見不鮮。拓撲優化作為一種靈活的結構優化手段，可以從本質上改變結構的拓撲形式，在滿足工作要求的前提下改善結構性能和減輕結構質量，經濟且美觀，已成為近年來結構優化領域研究的熱點。

1988年，Bends?e等[2]首先提出將均勻化方法引入到連續體拓撲優化中，揭開了連續體結構拓撲優化的開端。在此之后，變密度法、漸進結構優化方法(evolutionary structural optimization，ESO)和水平集方法(level set method，LSM)等[3-7]越來越多的優化算法不斷地被提出。拓撲優化方法的不斷發展給振動控制提供了新的思路，越來越多的學者開始將各種優化算法運用到振動控制的優化問題中。房占鵬等[8]基于雙向漸進優化算法，以最大化模態損耗因子為目標進行了約束阻尼板的拓撲優化問題研究；Zhang等[9]基于固體各向同性材料懲罰模型 (solid isotropic material with penalization，SIMP)插值，以最小化結構指定位置處的振動水平為目標對自由阻尼板阻尼材料的分布進行了優化；Zheng等[10]利用遺傳算法研究了以減少結構振動幅值為目標的阻尼材料的最優分布問題;Xia等[11]以最大化結構一階頻率為目標，利用水平集算法進行了拓撲優化設計。目前，在結構拓撲優化問題中，大多僅考慮結構自身而忽略了周圍介質對于結構響應的影響。當然，在結構處于類似空氣這類密度相對結構很低的介質環境下，這種忽略是合理且簡便的；而當結構處于水這類相對密度較大的介質中，則必須考慮流體介質和結構之間的相互作用，即強耦合作用。這類流固耦合系統下的拓撲優化工作已開始引起學者們的注意。商林源等[12]利用放松的優化準則法研究了聲結構耦合系統的雙材料拓撲優化問題；Shu等[13]基于水平集方法，考慮了最小化聲振耦合系統內聲場聲響應的優化問題。Akl等[14]以板的厚度為設計變量，考慮了板與聲腔的耦合系統的拓撲優化問題。由于有限元數值模擬方法對于無限大聲場分析具有一定的局限性，前面這些工作考慮的都是有限的聲場問題。邊界元方法由于可以自動滿足無限遠處的輻射邊界條件，對于處理無限大聲場問題具有天生的優勢，近年來基于有限元-邊界元耦合的外聲場問題的拓撲優化工作也得以實現[15-16]。目前上述工作考慮的目標函數都是聲壓值或輻射聲功率這類聲場物理量，忽略了對結構響應的考察；同時也多是基于密度插值和移動漸近線方法展開的，其他優化算法在這類優化問題中的效能還有待考察。

分段常數水平集(piecewise constant level set，PCLS)方法由水平集方法演化而來，最早由Christiansen等[17]提出。由于它可以很好的克服傳統水平集方法存在的依賴初始符號距離函數和無法自動生成新的孔洞等較為明顯的弊端，同時避免了重新初始化的過程，逐漸被一些學者用于結構拓撲優化中[18-19]。但是這類方法一般需要利用拉格朗日乘子法或增廣拉格朗日乘子法來滿足約束條件，而拉格朗日乘子的更新是較為復雜和麻煩的，難以準確表示其變化規律[20]，往往出現隨優化問題的變化而改變的情況，需要通過經驗來判斷，使方法的廣泛使用受到限制。

本文基于分段常數水平集方法，考慮了結構與聲場雙向耦合作用，以有限元-邊界元耦合系統下結構指定位置處振幅的平方為目標函數，進行了雙材料的拓撲優化設計。采用伴隨變量法進行靈敏度分析，引入二次罰函數方法使優化問題變為無約束優化問題，核心處理是基于靈敏度信息對優化參數進行了重新定義，克服了參數的問題依賴性。最后基于最速下降方法進行設計更新，最終得到了雙材料的最優分布。

1 有限元-邊界元聲振耦合分析

考慮結構與無限大聲場的強耦合系統，如圖1所示。結構Ω振動向聲場Ωf輻射聲壓，聲場反作用于結構，影響結構的振動。圖中：Γ=Ωs∩Ωf為結構與聲場的耦合邊界；ns，nf分別為結構與聲場耦合界面處的外法向向量；fs為直接作用在結構上的結構載荷。分別采用有限元和邊界元方法對結構和聲場進行離散近似。

圖1 聲振耦合系統Fig.1 Coupled structural-acoustic systems

1.1 結構有限元分析

對于簡諧激勵作用下的結構，其振動控制方程可以寫為

(-θ2M-iθC+K)u=Kdu=f

(1)

基于瑞利阻尼模型，阻尼陣C可以寫為

C=αM+βK

(2)

式中，α和β為材料的瑞利阻尼系數。

在考慮聲場對結構的反作用時，外載荷可以寫為

f=fs+fp

(3)

式中：fs為結構載荷；fp為聲場載荷。

1.2 聲場邊界元分析

簡諧形式下的聲場Helmhotz方程可以表示為

?2p(x)+k2p(x)=0

(4)

式中：p為聲壓；k=θ/cf為波數，cf為波速。通過格林第二公式并將點x逼近邊界即可得到常規邊界積分方程(conventional boundary integral equation，CBIE)

(5)

將式(5)兩邊同時對源點x所在邊界外法向方向求導，即可得到超奇異積分邊界方程(hypersingular boundary integral equation，HBIE)

(6)

式中：y為場點；q(y)=?p(y)/?n(y)為聲通量；系數c(x)由點x處的幾何性質決定，對于光滑邊界，c(x)=0.5；G(x,y)為格林函數，對于三維聲場問題，可以寫為

(7)

式中，r=|x-y|為源點與場點的歐式距離。

進行外聲場分析時，由于式(5)或式(6)單獨使用會引起解的非唯一性問題，因此Burton等[21]提出了將兩式進行線性組合的Burton-Miller方法

CBIE+αHBIE=0

(8)

式中，α為耦合系數，建議取為-i/k[22]。通過配點法對式(8)進行離散，可得聲場控制方程

Hp=Gq

(9)

式中，H和G為邊界元系數矩陣。

1.3 有限元-邊界元耦合分析

考慮結構與聲場的雙向耦合作用，需要考慮耦合面上位移與力的平衡條件

vf(x)=-iθu(x)·nf(x)

(10)

σ(x)·ns(x)+p(x)·nf(x)=0

(11)

式中：vf為耦合邊界處聲場介質的法向速度；σ(x)為耦合邊界處結構應力。在式(10)、式(11)兩端分別乘以結構和聲場的插值形函數φf，φs并在耦合邊界上積分、離散后得到耦合邊界平衡方程

vf=-iθS-1Cfsu

(12)

fp=Csfp

(13)

式中：Csf和Cfs為結構與聲場耦合轉換矩陣；S為邊界質量矩陣，可以表示為

(14)

(15)

(16)

式中：Ns，Nf為形函數φs，φf的矩陣形式；n為nf的矩陣形式。

又由第二類邊界條件

q=iθρfvf

(17)

式中，ρf為聲場介質密度。將式(12)、式(17)代入式(9)中可得

Hp=θ2ρfGS-1Cfsu

(18)

聯立式(1)、式(18)即可得到基于有限元-邊界元方法的聲振耦合系統控制方程

(19)

2 優化問題描述

如圖2所示，對于雙材料的拓撲優化，優化的目的是要確定兩種材料的最優分布形式。在本文中，我們以材料二作為設計材料，即確定材料二的分布形式，其余部分則全設為材料一。在分段常數水平集方法中，我們將引入PCLS函數φ(x)來表示兩種材料的鋪設區域

圖2 基于PCLS方法的雙材料分布Fig.2 Bi-material distributions based on PCLS method

(20)

即在整個設計區域Ω內，PCLS函數保持為常數，在材料一鋪設區域Ω1為1，在材料二鋪設區域Ω2為2。

基于PCLS函數的單元剛度陣K(e)與質量陣M(e)可分別表示為

K(e)(φ)=(φ-1)K2(e)-(φ-2)K1(e)
M(e)(φ)=(φ-1)M2(e)-(φ-2)M1(e)

(21)

基于瑞利阻尼假設，單元阻尼陣C(e)可以表示為

C(e)(φ)=(φ-1)(α2M2(e)+β2K2(e))-
(φ-2)(α1M1(e)+β1K1(e))

(22)

式中：上標1和2分別為材料一與材料二；α1，β1，α2，β2分別為兩種材料的瑞利阻尼系數。

于是，基于PCLS函數的結構整體剛度陣、質量陣與阻尼陣可以表示為

(23)

式中，Ne為結構單元個數。

而作為設計材料的材料二體積可以表示為

(24)

因此，在考慮材料體積約束時，基于分段常數水平集方法的以結構振幅平方為目標的拓撲優化問題可以表示為

s.t.Kd(φ)u(φ)=fs+Csfp(φ)

Hp(φ)=θ2ρfGS-1Cfsu(φ)

(25)

3 靈敏度分析

基于分段常數水平集方法，目標函數依賴于PCLS函數φ(x)，在計算靈敏度信息時，我們需要對目標函數進行變分。本文中，靈敏度分析采用的是Zhao等使用的伴隨變量法。

首先引入原目標函數的等價形式

(26)

(27)

由于伴隨向量的任意性，有

(28)

(29)

(30)

基于式(1)、式(21)、式(22)可得

(31)

[α2M2(e)+β2K2(e)]-[α1M1(e)+β1K1(e)],

(32)

于是，基于分段常數水平集方法的單元靈敏度可以表示為

(33)

針對本文的目標函數，則有

(34)

聯立式(28)～式(34)，即可得到基于分段常數水平集方法的單元靈敏度。

4 算法流程

對于體積約束條件，我們將引入二次罰函數的方法來滿足體積約束，最終修正的目標函數可以寫為

(35)

式中，μ為體積懲罰系數。修正后目標函數的單元靈敏度可以寫為

(36)

不同于水平集算法的是，在每代優化中，PCLS函數都是定義為1或2，因此需要定義一個算法，使PCLS函數再次滿足約束式(20)，即

(37)

我們通過最速下降方法來更新PCLS函數，即在每個單元中

(38)

式中，τ(k)為時間步長，在第k代迭代步中，設為

(39)

式中：h為單元的最小尺寸;ζ∈(0,1)為一常數。

而在第k代迭代步中，Zhang等(2019年)將體積懲罰系數μ(k)設為

(40)

式中，t和μmax為給定常數，用于保證體積約束得以滿足。然而，μ(k)在這種定義下是與目標函數J(φ(k))相關的，在計算的過程中可以發現，參數t與μmax存在較大的問題依賴性，即隨著優化算例的改變，參數的變化較大，不具有一般性，往往需要通過調試才可獲得。Zhang等研究中不同優化算例選擇不同的優化參數同樣可以看出這一點。因此，我們從靈敏度的角度出發，對參數μ(k)的數值進行了重新定義

(41)

式中，ε∈(0,1)為一常數。由于PCLS函數的更新是直接依賴于靈敏度信息的，式(41)中μ(k)的數值則直接與目標函數靈敏度絕對值的最大值相關，因此更具有一般性。參數ε用于保證優化可以穩定光滑的進行。經過測試，本文將ζ均取為0.05；ε均取為0.2；μmax均取為10。

至此，基于分段常數水平集方法的拓撲優化流程總結如下：

步驟1建立有限元-邊界元耦合模型，初始化PCLS函數，獲得結構有限元與聲場邊界元矩陣M，K,C,H,G,Cfs，Csf確定常系數矩陣Q，給定常數ζ,ε,μmax。

步驟2通過求解式(19)獲得結構位移與聲場聲壓，進而由式(25)計算目標函數J(φ)。

步驟4引入內循環，更新φ(k+1)。

①計算

式中，τ(k)由式(39)計算得到。

步驟5檢查外循環終止條件，如果滿足，則退出循環；否則回到步驟2。

5 算 例

本章中，我們考慮一水下立方殼體在點激勵作用下的優化算例，立方殼外部空間為水，不考慮立方殼內部介質的耦合影響，約束結構底面的平動自由度，如圖3所示。立方殼體邊長為1 m，殼厚度為0.02 m，在其上表面中心有一點簡諧激勵F=F0e-iθt，F0=105N頻率f=120 Hz。結構每個面劃分為40×40個四節點殼單元，共9 600個單元，聲場部分采用匹配的四節點常量元與結構進行耦合。

圖3 立方殼結構及網格劃分Fig.3 The cubic shell structure and the mesh

假設結構由兩種不同的強弱材料組成，弱材料的楊氏模量和密度為強材料的十分之一，具體的材料屬性如表1所示。假設兩種材料的瑞利阻尼系數相同，都設為α=0.1,β=0.001。由于結構上表面對于結構響應的影響最大，因此選擇立方殼的上表面作為設計區域。結構其余部分全部設為弱的基體材料，且其余部分同上表面一樣可以振動向外輻射發聲，以強材料為設計材料。目標函數選為加載點處z向幅值的平方，體積約束V=|Ω|/2。

表1 水下立方殼材料參數

為驗證本文所用伴隨變量法的正確性，我們將PCLS函數φ(x)放松到[1,2]的連續空間，采用差分法計算目標函數的靈敏度同本文的伴隨靈敏度進行比對。圖4給出了120 Hz時初始全為弱材料，攝動值為1×10-4時前200個單元的伴隨靈敏度與差分結果的比對結果。可以看出二者吻合良好，說明了本文伴隨變量方法的正確性。

圖4 伴隨變量法與有限差分法靈敏度結果對比Fig.4 Comparison of sensitivities obtained by the adjoint method and finite difference method

初始PCLS函數設為φ(0)=1，即結構上表面全由弱材料組成。目標函數與體積約束函數的收斂歷史，如圖5所示。為了考察目標函數的收斂性，在已經收斂的情況下繼續往后共計算了90代。從圖5中可以看出，隨著迭代的進行，目標函數迅速下降，經過62次迭代后，目標函數和體積約束函數都已收斂。目標函數從5.019×10-5m2降到了1.408×10-6m2。為進一步說明優化的有效性，計算了結構上表面全為強材料時的目標函數，其值為2.374×10-6m2，高于體積約束為0.5時的優化設計值，也說明了該優化問題的非線性。

圖5 120 Hz時目標函數與體積約束的迭代歷史Fig.5 Iteration histories of objective function and volume constraint at 120 Hz

圖6給出了不同迭代步數下材料的分布情況，圖6中淺灰色部分代表強材料，黑色部分代表弱材料。可以看出強材料的占比在穩步提高，最終優化后的材料分布如圖6(d)所示。優化前后的結構z向振幅云圖，如圖7所示，可以看出優化的效果顯著。

圖6 水下強耦合時不同迭代步數下材料的優化分布Fig.6 The optimal distributions at different iteration steps with strong coupling case when the structure immersed in water

圖7 優化前后結構z向振動幅值Fig.7 z-direction vibration amplitudes before and after the optimization

為了體現PCLS方法的優越性，我們基于SIMP插值和MMA算法進行了同樣的優化設計，優化結果如圖8所示。目標函數降到了1.619×10-6m2，高于本文算法的優化值，這同2019年Zhang等研究中的結論保持一致。

圖8 MMA算法優化的材料分布Fig.8 The optimal distribution based on MMA method

為了說明耦合條件和環境介質對優化結果的影響，我們又考慮了空氣強、弱耦合和水下弱耦合條件。空氣強、弱耦合條件下優化結果完全相同，都為圖9所示的分布。可以看出，對于空氣這類密度相對結構小很多的介質，此時忽略介質對于結構的影響是合理的。弱耦合條件即單向耦合條件，不考慮聲場介質對于結構振動的影響，因此水下弱耦合條件時材料的優化結果與空氣弱耦合時的相同，即如圖9所示。可以看出，水下強耦合條件相對空氣強、弱耦合和水下弱耦合條件時強材料在設計域四周的分布存在明顯的差異。

圖9 空氣強、弱耦合與水下弱耦合條件下的材料優化分布Fig.9 The optimal distribution with strong or weak coupling case when the structure immersed in air, or with weak coupling case when the structure immersed in water

為了進一步說明介質為水時，強、弱耦合對于振動響應的影響，圖10給出了聲場介質為水、初始全為弱材料時，強、弱耦合條件下設計區域z向振幅云圖的對比。可以看出，介質為水時強、弱耦合條件下振動響應明顯不同，這也說明了介質為水這類相對密度較大的介質時，結構振動響應的計算必須考慮強耦合條件，否則會引起較大的誤差甚至計算錯誤。

圖10 120 Hz水下不同耦合條件時設計域z向振動幅值Fig.10 z-direction vibration amplitudes of design domain with different coupling cases when the structure immersed in water at 120 Hz

針對體積約束參數的問題依賴性，我們將參數按式(40)定義方式同樣計算了120 Hz時強耦合的優化結果，設t=1.5,μmax=2×10-7，同本文定義參數的收斂歷史對比如圖11所示，可以看出，本文提出的參數定義方式優化收斂較快。更重要的是，原參數定義方式會出現參數隨算例變化而改變的問題，如在弱耦合時μmax需要設為2×10-5，而本文所有算例的體積約束均在二次罰函數系數取為固定值時得到了很好的滿足，說明改進后的方法明顯改善了參數的問題依賴性，至少在結構振幅優化的問題中有效，有利于算法的進一步推廣。

圖11 強耦合時參數修改前后收斂歷史對比Fig.11 Comparison of convergence histories before and after the parameter modification with strong coupling case

圖12給出了不同激勵頻率下的優化結果。可以明顯看出，優化結果隨頻率變化而變化。在30 Hz時，優化結果較為分散。我們分別對結構為初始設計與優化設計時在30 Hz附近進行了掃頻分析，如圖13所示。可以發現，當激勵頻率為30 Hz時，初始設計下結構響應出現了峰值，這可能由于該頻率與耦合結構的固有頻率較為接近引起的。在優化設計時，結構共振峰頻率出現了右移，由30 Hz提高到了47 Hz左右，遠離了激勵頻率，使加載點位移幅值得到大幅降低。同時可以看出，在其他激勵頻率時優化設計的加載點振幅高于初始設計，說明優化設計只是設計激勵頻率下的最優設計，不適用于其他頻率。

圖12 不同激勵頻率下的材料優化分布Fig.12 The optimal distributions at different excitation frequencies

圖13 30 Hz附近掃頻響應分析Fig.13 Swept-frequency analysis around 30 Hz

當激勵頻率在60～200 Hz時，隨著頻率的提高，強材料有向結構中心聚集的趨勢。

最后，我們考察了瑞利阻尼系數β對于優化結果的影響，如圖14所示。在考察參數修改、強弱耦合條件和介質等因素對于優化結果的影響時，頻率設為120 Hz、α設為0.1，因此這里我們將頻率固定設為120 Hz、α固定設為0.1。同樣可以看出，隨著的增大，強材料有向結構中心聚集的趨勢。

圖14 不同阻尼系數β時的材料優化分布Fig.14 The optimal distributions with different damping coefficents β

6 結 論

本文研究了簡諧激勵作用下，考慮結構與聲場強耦合的水下立方殼體雙材料的優化分布問題。采用四節點殼單元和四節點常量元進行結構和聲場的離散，基于PCLS方法和瑞利阻尼假設，構造出結構整體剛度陣、質量陣與阻尼陣；目標函數選為結構指定位置處振幅的平方，采用伴隨變量方法進行靈敏度計算，降低了計算成本。采用二次罰函數方法使體積約束條件得以滿足，對于算法中的體積約束參數進行了重新定義，降低了體積約束參數對問題的依賴性，有利于算法的推廣。數值結果表明優化設計可以明顯降低結構的振幅，驗證了優化方法的有效性。對頻率以及瑞利阻尼系數β對于優化結果的影響也進行了進一步的討論，發現了在計算頻段內隨頻率和β的提高，設計材料呈現向結構中心聚集的趨勢。

本文針對的是單個頻點時的優化情況，可以進一步考察一個區間內的頻段優化問題。