考慮材料性能不確定性的聲振強耦合系統穩健拓撲優化

鄭文治, 陳海波, 操小龍

(1. 中國科學技術大學 近代力學系 中國科學院材料力學行為與設計重點實驗室,合肥 230027;2. 北京機電工程研究所,北京 100074)

拓撲優化由于其可以相當靈活地改變結構的拓撲構型,在振動與噪聲控制領域受到了廣泛關注。Du等[1]以輻射聲功率最小為目標,對平板的材料分布進行了拓撲優化設計;之后Du等[2]將此優化框架擴展到微結構的拓撲優化設計;Zhang等[3]以目標點的聲壓最小為目標,基于固體各向同性材料懲罰(solid isotropic material with penalization,SIMP)模型對振動結構的阻尼材料分布進行了優化設計;吳振云等[4]使用非負聲強識別結構表面對遠場聲輻射起主要貢獻的區域,對約束阻尼板阻尼材料的最優分布進行研究。目前,大多數聲振優化的研究僅考慮弱耦合,即僅考慮結構對聲場的作用而忽略了聲場對結構的反作用。對于像空氣這種流體密度遠低于結構的聲學介質來說,這種簡化通常是合理的,但對于像水這種重流體,則必須考慮結構和聲學介質之間的相互作用,即強耦合作用。已有學者針對強耦合情況開展了研究工作。Shu等[5]基于水平集方法對聲振耦合系統進行了拓撲優化設計;Chen等[6]以目標點的聲壓級最小為目標,對聲振耦合系統中的復合材料板的微結構進行了優化。上述兩項研究都基于有限元法,針對有限域的內聲場問題開展的。對于無限域外聲場問題,有限元法需要對較大的聲場域進行離散,計算效率大大降低,且無限遠處的無反射條件難以精確滿足,需要進行麻煩的特殊處理[7-8]。邊界元法只需在聲場邊界進行離散,且自動滿足無限遠處聲場的無反射條件,對于分析外聲場問題有很大的優勢,因此基于有限元-邊界元耦合的分析方法可以有效地用于解決外聲場與結構的耦合問題[9-10]。近些年,針對外聲場強耦合問題的拓撲優化的研究工作也已開展起來[11-12]。Zhao等采用有限元-邊界元耦合方法,以輻射聲功率級最小為目標,基于材料屬性的有理近似方法(rational approximation of material properties, RAMP)模型和移動漸近優化算法(method of moving asymptotes,MMA)對聲振耦合系統的材料分布進行了拓撲優化設計。

上述的拓撲優化都是在確定性的條件下進行的。然而在實際工程中,由于制造、裝配、外部載荷的不可完全預測等因素,導致材料屬性、幾何尺寸、外部載荷等存在不確定性。傳統的確定性拓撲優化忽略了系統的不確定性,可能導致優化的結果是次優的甚至完全不符合實際。因此,開展考慮不確定性的拓撲優化研究工作是很有必要的。

考慮不確定性因素的拓撲優化通常分為兩種:可靠性拓撲優化[13-14]和穩健拓撲優化。可靠性拓撲優化關注的是在極端條件下系統的工作狀況,而穩健拓撲優化側重于在一定程度上減小目標函數對不確定性參數的敏感性。趙清海等[15]考慮載荷的不確定性,針對多材料結構進行了穩健拓撲優化;李冉等[16]基于隨機有限元方法進行不確定性分析,考慮增材制造工藝中引起的表面層厚度不確定性進行了穩健拓撲優化;Yin等[17]基于隨機攝動方法進行不確定性分析,將聲介質的物理參數、外部載荷和材料特性視為隨機參數,對聲振耦合系統的微結構進行了穩健拓撲優化設熟計。Keshavarzzadeh等[18]基于混沌多項式展開(polynomial chaos expansion,PCE)方法量化不確定性的傳播和響應的靈敏度,建立可靠性拓撲優化設計方法和穩健拓撲優化設計方法;Zhang等[19]將PCE方法拓展到聲子晶體的穩健拓撲優化中;Torii等[20]使用該方法處理了考慮不確定量的發熱、不確定位置的發熱和不確定位置損傷的熱傳導穩健拓撲優化;Kang等[21]考慮材料梯度界面梯度的不確定性,基于PCE方法和水平集方法對多材料結構進行了穩健拓撲優化設計。以上研究均未涉及外聲場強耦合聲振系統問題,我們的文獻調研也未見相關報道。

本文針對外聲場與結構強耦合問題,建立了考慮材料性能空間不確定性的穩健拓撲優化模型。使用RAMP模型來描述材料分布;假設材料的彈性模量服從高斯隨機場分布,通過級數最優線性估值(expansion optimal linear estimation,EOLE)方法將彈性模量隨機場離散成有限個隨機變量;基于PCE方法進行隨機響應預測和隨機響應的靈敏度分析;采用MMA算法進行設計更新,最終得到最優的材料布局。最后通過數值算例說明了考慮材料性能不確定性的必要性和本文所建立方法的有效性。

1 有限元-邊界元分析方法

1.1 結構有限元分析

在簡諧激勵作用下,有限元離散后的結構振動方程可表示為

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

(1)

式中:K,C和M分別為結構剛度矩陣、阻尼矩陣和質量矩陣;u為節點位移;f為等效節點載荷;ω為激勵圓頻率。

考慮聲場對結構的反作用時,外載荷表示為

f=fs+fp

(2)

式中:fs為直接加載到結構的載荷;fp為聲場反作用到結構上的載荷。將式(2)代入式(1)可得

Kdu=fs+fp

(3)

采用瑞利阻尼模型,結構的阻尼矩陣可以表示為剛度矩陣和質量矩陣的線性組合

C=α0M+β0K

(4)

式中,α0和β0為瑞利阻尼模型的相關系數。

1.2 聲學邊界元分析

在均勻理想流體介質中,聲場的控制方程為Helmholtz方程

?2p(x)+k2p(x)=0, ?x∈Ωf

(5)

式中:p為聲壓;k=ω/cf為波數,ω為圓頻率,cf為波速。應用格林第二公式,并將源點趨近于邊界,可得常規邊界積分方程(conventional boundary integral equation, CBIE)

(6)

式中:x為源點;y為場點;q(y)=?p(y)/?n(y)為通量,n(y)為邊界點y處的外法線方向向量;系數c(x)的大小由點x處的幾何性質決定,當點x處為光滑邊界時,c(x)=0.5;對于三維聲學問題,其基本解為

(7)

式中,r=|x-y|為源點x和場點y的距離。

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

(8)

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

CBIE+αHBIE=0

(9)

式中,α為耦合系數,本文取為-i/k[23]。對聲學邊界進行離散,得到式(9)的離散格式

Hp=Gq

(10)

式中:H和G為邊界元系數矩陣;p和q分別為邊界節點聲壓及其法向通量。

1.3 有限元與邊界元耦合

結構與聲學耦合面上需滿足位移連續與力的平衡條件

q=iωρfvf=ω2ρfS-1Cfsu

(11)

fp=Csfp

(12)

式中:ρf為流體介質密度;vf為聲場法向速度,Cfs和Csf為結構與流場之間的耦合矩陣;S為邊界質量矩陣。

(13)

(14)

(15)

將式(11)、式(12)分別代入式(10)、式(3)中,得到聲振耦合系統控制方程

(16)

通過求解上述方程,就可以完成對聲振耦合系統的分析。整個系統對外輻射做功的大小可以通過輻射聲功率表示

(17)

式中:W為輻射聲功率; 上標()H為共軛轉置; Re()為取實部。輻射聲功率級可以表示為

(18)

式中,W0為參考聲功率,本文取W0=2×10-12W。

2 隨機場模型及隨機響應分析

2.1 隨機場離散

本文考慮采用高斯隨機場t(x)描述材料彈性模量的不確定性,其均值為μ(x),其協方差函數為

(19)

式中: cov(x1,x2)為隨機場的協方差函數;σ為隨機場的標準差;l為隨機場的相關長度; |x1-x2|為點x1和點x2的距離。

首先,將隨機場域離散為一系列節點x1,…,xn,則隨機場的協方差矩陣為

(20)

基于EOLE方法[24],隨機場可以表示為

(21)

式中:Cv=[cov(x1,x),cov(x2,x),…,cov(xn,x)]T;λk,ψk分別為協方差矩陣Cm的特征值和與之對應的特征向量;ξ1,ξ2,…,ξk,…,ξn為互不相關的隨機變量,其均值和方差滿足

E[ξi]=0, E[ξiξj]=δij

(22)

式中: E[]為期望;δij為Kronecker delta符號。

為了減小問題的維數,將特征值λk降序排列,對式(21)取前s(s≤n)項進行截斷

(23)

隨機場的截斷精度誤差為

(24)

根據式(23)就可以實現對隨機場的有限離散。

2.2 基于PCE方法的隨機響應分析

對隨機場離散后,得到一系列互不相關的隨機變量ξ=[ξ1,ξ2,…,ξs]T, 聲振耦合系統的隨機響應通過PCE方法進行預估

(25)

式中:Y(ξ)為系統響應;yi為多項式的系數;φi(ξ)為多項式的正交基函數,正交于隨機變量ξ的概率密度函數。對于不同的概率分布,多項式的形式會有不同,例如當隨機變量ξ服從均勻分布時,φi(ξ)為Legendre多項式; 當隨機變量ξ服從高斯分布時,φi(ξ)為Hermite多項式。正交基函數φi(ξ)滿足

(26)

對于多維隨機變量ξ=[ξ1,ξ2,…,ξs]T,多項式函數φi為

(27)

式中:καj為隨機變量ξj的αj階多項式;φi的階數為

(28)

為了在保證計算精度的前提下,有效降低計算量,一般將式(25)截斷為有限項

(29)

式中,多項式系數的項數N+1取決于隨機變量的維數s和PCE的階數m

(30)

求解混沌多項式展開的關鍵在于求解其中的系數yi,利用多項式φi的正交特性可得

(31)

由式(31)可得

(32)

(33)

式中:ρ(ξ)為隨機變量ξ的概率密度函數;ξ(q),w(q)和nq分別為積分點、權值和積分點數目。由于本文采用的是高斯隨機變量,所以采用Gauss-Hermite數值積分計算,由此可求解系數。

計算得到式(29)的系數yi后,響應Y(ξ)的均值E[Y(ξ)]和標準差σ[Y(ξ)]可以通過下式獲得

(34)

(35)

3 穩健拓撲優化設計

3.1 確定性拓撲優化數學模型

拓撲優化是通過優化材料的空間分布,在滿足約束的前提下最大程度的改善結構性能。聲振耦合系統一般選擇位移幅值、聲壓級、聲功率級作為目標函數,其數學模型表達為

(36)

式中: 目標函數為結構位移u和聲壓p的實函數;ρ=[ρ1,ρ2,…,ρNe]T為設計變量;Ne為設計變量總數;ve為第e個單元的體積;fv為體積比約束。對于本文研究的聲振耦合問題,除了滿足上述的材料體積約束之外,還要滿足式(16)的控制方程。

3.2 穩健拓撲優化數學模型

與確定性拓撲優化相比,穩健拓撲優化考慮了隨機變量對結構響應的影響,其目標是尋找對隨機變量不敏感的優化設計,目標函數的形式通常是隨機響應的均值和其標準差的加權和,屬于多目標優化設計。其數學模型表述為

(37)

式中,k為隨機響應標準差的權值因子。

3.3 材料插值模型

本文通過優化雙材料的分布,在滿足約束的條件下達到減振降噪的效果。本文采用RAMP插值方法來描述雙材料的分布并建立設計變量ρ和材料分布的關系,可以將單元矩陣表示為

(38)

(39)

式中:Ke,Me分別為第e個單元的剛度矩陣、質量矩陣;Ke1和Ke2,Me1和Me2是與分別給定的兩個材料相關的剛度矩陣和質量矩陣;ρe=1為該單元完全由材料1(強材料)組成,ρe=0為該單元完全由材料2(軟材料)組成;η為讓中間密度趨近于0或1的懲罰因子,通常取5。

3.4 確定性條件下的靈敏度分析

確定性拓撲優化的目標函數對于設計變量ρe的導數可以表示為

(40)

(41)

式中,λ1,λ2為伴隨乘子。則目標函數對于設計變量的靈敏度可以表示為

(42)

由于G,H只與頻率、幾何形狀有關;Csf,Cfs,S只與幾何現狀有關;均與設計變量無關。本文考慮的載荷fs也與設計變量無關。因此,式(42)可以簡化為

(43)

(44)

(45)

則,式(43)進一步簡化為

(46)

3.5 隨機響應的靈敏度分析

隨機響應相對于設計變量的靈敏度也采用PCE方法進行表示

(47)

式中,系數zi可根據與式(32)類似的方法獲得

(48)

(49)

(50)

由此可得,穩健拓撲優化問題的目標函數的靈敏度為

(51)

3.6 穩健拓撲優化設計流程

基于以上隨機響應的靈敏度信息,采用MMA算法進行設計變量更新。收斂條件為當兩個相鄰迭代步的目標函數的相對比值小于一個特定的常數τ

(52)

式中,Ji為第i步迭代對應的目標函數。本文中取τ=5.0×10-5。為了有效消除中間密度單元,使單元的相對密度趨近于0～1分布,本文采用保體積密度過濾方法[25],且將過濾半徑設為0.05 m。該過濾方法中的β參數會有效影響過濾的程度,β越大,過濾越明顯,但是β取值過大可能導致最終結果收斂至局部最優解。本文中,在前20個迭代步,β=0;在迭代20步之后,迭代每增加一步,β增加0.5,最大不超過100。詳細的優化過程如圖1所示。

圖1 聲振耦合系統穩健拓撲優化流程圖Fig.1 Flow chart of robust topology optimization of coupled structural-acoustic systems

4 數值算例

本章中,考慮水下立方殼在點激勵作用下的拓撲優化問題。立方殼邊長為1 m,上表面為厚度0.02 m的四邊固支板,其余各面是剛性的。在上表面中心處有幅值為1 N的簡諧力作用,力的方向垂直于上表面豎直向下,如圖2所示。在殼的上表面劃分80×80的四節點板單元進行結構分析,對于聲場部分,在殼的每一個面采用40×40的四節點常量元進行聲場分析。立方殼的外部介質為水,其聲速為cf=1 482 m·s-1,密度為ρf=1 000 kg·m-3,不考慮立方殼內部的耦合影響。

圖2 立方殼結構Fig.2 The cubic shell

4.1 克服網格依賴性驗證及體積約束影響測試

Sigmund等[26]指出采用獨立于網格的過濾技術可以克服網格依賴性問題。本文采用了保體積過濾技術進行密度過濾,其過濾半徑是獨立于網格尺寸的。這里,我們先通過一個算例測試來說明本文的優化方法是沒有網格依賴性問題的。以輻射聲功率級為目標函數,材料1的體積分數約束fv=0.5;設計變量的初始值全部設置為1,即設計域全由材料1構成。分別以80×80和60×60的結構網格對在300 Hz激勵下的立方殼進行拓撲優化,最終得到的拓撲優化結果如圖3所示。可以看出,采用不同網格密度進行優化,最終得到的拓撲結構是基本相同的。80×80和60×60網格的最終優化結果對應的輻射聲功率級分別為60.328 dB和60.340 dB, 相差很小。因此在本文中網格密度的差異不會對拓撲優化的結果產生影響。本文后續算例將采用80×80的結構網格進行計算。

圖3 不同網格密度下的拓撲優化結果Fig.3 Topology optimization results under different grid densities

為了說明體積約束分數fv對優化結果的影響,這里我們以輻射聲功率級為目標函數,分別在材料1的體積分數約束fv取0.5,0.6,0.7時對300 Hz激勵下的立方殼進行拓撲優化。表1給出了不同體積分數約束下的優化結果對應的輻射聲功率級以及優化迭代步數。可以看出隨著fv的取值變大,輻射聲功率級是降低的,這是因為強材料的剛度更大,強材料的應用有利于降低振動幅值,進而降低輻射聲功率級。

表1 不同體積分數約束fv下拓撲優化結果對應的輻射聲功率級和優化迭代步數

4.2 隨機響應分析的精度驗證

在介紹穩健拓撲優化之前,我們先通過一個數值算例驗證本文所用的PCE方法進行隨機響應分析的正確性。

我們選擇4階的PCE方法進行隨機響應分析,其中計算PCE系數過程中需要計算的積分,即式(33),我們通過5階Gauss-Hermite積分的張量積格式進行計算。表2給出了當設計域全為材料1時,即全為強材料時,采用PCE方法分別在頻率為50 Hz,100 Hz,150 Hz,200 Hz,250 Hz,300 Hz,400 Hz激勵下的輻射聲功率級的均值和標準差的計算結果,并將其與取10 000個樣本點下的蒙特卡洛方法(Monte Carlo method, MCM)得到的結果進行比較。從表2結果可以看出,在以上所述的各個頻率激勵下,PCE方法得到的隨機響應的均值和標準差和MCM得到的結果非常相近。由此可知,采用4階PCE方法可以獲得相當精確的隨機響應結果,因此在后續的穩健拓撲優化設計中使用PCE方法進行隨機響應分析是可靠的。

表2 PCE和MCM隨機響應結果對比

4.3 穩健拓撲優化設計結果

首先對立方殼在300 Hz激勵下的情況進行穩健拓撲優化設計,其目標函數取輻射功率級的均值和標準差的加權和,如式(37)所示,其中權重因子k=3;材料1的體積分數約束fv=0.5;設計變量的初始值全部設置為1,即設計域全由材料1構成。在優化44步迭代后目標函數和體積約束函數收斂。目標函數和體積約束以及輻射聲功率級的均值和標準差的迭代歷史如圖4所示。在最初的幾次迭代中,強材料的體積分數從1降至0.5,并保持在0.5附近。由于我們設計變量的初始值全部設為1,目標函數在開始時有很大的增加,但是在接下來的迭代中,目標函數、均值、標準差基本在穩定的減小,從而實現對材料分布不確定性的不敏感設計。穩健拓撲優化設計的輻射聲功率級的均值為60.336 dB,標準差為0.218 0 dB。為了體現穩健優化設計的優越性,我們以輻射聲功率級作為目標函數,同樣基于RAMP材料插值模型和MMA算法進行了確定性拓撲優化設計。考慮材料性能的不確定性,對獲得的確定性拓撲優化結果進行隨機響應分析,其輻射聲功率級的均值為60.351 dB,標準差為0.236 7 dB。穩健拓撲優化設計的輻射聲功率級的均值和標準差均小于確定性拓撲優化結果,穩健拓撲優化結果比確定性拓撲優化具有更好的隨機性能,并且對材料性能的不確定性不太敏感。穩健拓撲優化和確定性拓撲優化的結果,如圖5所示。圖5中,黑色表示材料1,淺灰色表示材料2,兩者的結果僅有在一些細節上有所差別。

圖4 300 Hz激勵下穩健拓撲優化的迭代歷史Fig.4 Iterative history of robust topology optimization under 300 Hz excitation

圖5 300 Hz激勵下的穩健優化和確定性優化結果Fig.5 Robust and deterministic optimization results under 300 Hz excitation

為了進一步研究穩健拓撲優化設計的性能,我們在與上述相同的條件和參數下又對立方殼分別在200 Hz,250 Hz,400 Hz激勵下進行了穩健拓撲優化設計和確定性拓撲優化設計。圖6給出了相應的確定性拓撲優化和穩健拓撲優化的結果。其中200 Hz激勵下的穩健拓撲優化結果和確定性拓撲優化僅在一些細節上存在差別,但是對于250 Hz和400 Hz,穩健優化拓撲結果和確定性拓撲優化存在明顯差異,這說明考慮材料彈性模量的不確定性對聲振系統的拓撲優化結果有較大的影響,也說明考慮材料彈性模量的不確定性很有必要。為了更詳細地比較確定性優化和穩健優化的結果,表3給出了在各個頻率激勵下穩健拓撲優化設計和確定性拓撲優化設計對應的輻射聲功率級的均值和標準差。可以看出,與確定性拓撲優化相比,穩健拓撲優化結果的標準差都較小,這表明穩健拓撲優化設計對材料彈性模量的不確定性的敏感性更低,不僅如此,穩健拓撲優化設計的均值也都比確定性拓撲優化的均值更低。

表3 確定性設計和穩健設計的輻射聲功率級的均值和標準差

圖6 不同激勵頻率下的穩健和確定性優化結果Fig.6 Robust and deterministic optimization results at different excitation frequencies

為了研究穩健拓撲優化的目標函數中的標準差的權值因子k對優化結果的影響,以5種不同的權重因子對在400 Hz激勵下的水下立方殼進行穩健拓撲優化設計,其他參數設置均相同。表4中給出了采用不同權值穩健拓撲優化結果對應的輻射聲功率級的均值和標準差。結果顯示優化結果的輻射聲功率級的均值隨著權值因子k的逐漸增大而增大,標準差逐漸減小,這和雙準則優化問題帕累托解的基本特征是吻合的。同時也表明隨著權值因子k的增大,目標函數中標準差的權重變得更大,有利于降低優化結果對材料彈性模量的不確定性的敏感性。圖7給出了k取1,2,4,5時的穩健拓撲優化結果,k取3時的結果已在圖6(e)中給出。從這些穩健優化得到的結果上可以看出,不同的權值因子k得到的優化結果存在明顯的差異。綜上所述,權值因子k對優化結果有較大的影響,選擇合適的權值因子k對于穩健優化設計至關重要。

表4 400 Hz激勵下不同權值因子k下的穩健設計對應的輻射聲功率級的均值和標準差

圖7 400 Hz激勵下不同權值因子k的穩健優化結果Fig.7 Robust optimization results of different weight factors k under 400 Hz excitation

為了研究兩材料彈性模量的變異系數的影響,在其他參數保持不變的情況下,以4種不同的變異系數γ對在400 Hz激勵下的立方殼進行了穩健拓撲優化設計,目標函數中的權重因子k都取3。圖8給出了變異系數為0.06,0.16,0.25時對應的穩健拓撲優化結果,γ=0.10對應的優化結果已在圖6(e)給出。當不確定性較小時(γ=0.06),穩健拓撲優化結果和確定性拓撲優化結果較為相近,隨著變異系數的增大,穩健拓撲優化結果和確定性拓撲優化的差異也越來越大。表5中給出了變異系數γ取不同值時,穩健拓撲優化設計對應的輻射聲功率級對應的均值和標準差。

表5 400 Hz激勵下不同變異系數γ下的穩健設計對應的輻射聲功率級的均值和標準差

圖8 400 Hz激勵下不同變異系數γ下的穩健優化結果Fig.8 Robust optimization results with different coefficients of variation γ under 400 Hz excitation

5 結 論

本文研究了考慮材料性能的不確定性情況下,外聲場與結構強耦合情況下的穩健拓撲優化問題。通過隨機場模型描述兩材料的彈性模量的不確定性,采用EOLE方法將隨機場離散為有限個不相關的隨機變量;采用PCE方法進行隨機響應預估和隨機響應的靈敏度分析;采用RAMP模型進行密度插值,建立單元相對密度和彈性模量、材料密度之間的關系;穩健拓撲優化的目標函數為輻射聲功率級的均值和標準差的加權和;采用MMA算法對優化問題進行更新迭代求解。數值結果表明:

(1) 由于采用了獨立于網格尺寸的過濾方法,提出的拓撲優化方法沒有網格依賴性問題。

(2) 與確定性拓撲優化設計相比,穩健拓撲優化設計對應的輻射聲功率級的標準差更低,實現了對材料性能不確定性更加不敏感的設計,證實了本文建立的考慮材料性能不確定性的穩健拓撲優化模型的有效性。

(3) 隨著權值因子k的增大,穩健拓撲優化設計的輻射聲功率級的均值逐漸增大,標準差逐漸減小,且權值因子k對穩健拓撲優化結果的影響也較大。選擇合適的權值因子k很重要。