考慮可制造性約束的聲子晶體多目標拓撲優化1)

曹蕾蕾 武建華 樊 浩 張傳增 孫霖霖

* (長安大學道路施工技術與裝備教育部重點實驗室,西安 710064)

? (德國錫根大學土木工程系,德國錫根D-57068)

引言

聲子晶體是一種新型周期性人工超材料,它可以使得某些頻率范圍內彈性波的傳播衰減或禁止,其特有的帶隙特性使其在負折射、聲聚焦、聲吸收、聲隱身、聲成像領域具有廣闊的應用前景[1-2].帶隙的產生和分布是聲子晶體應用的關鍵問題,因此,利用拓撲優化方法進行聲子晶體帶隙設計與優化引起了國內外學者的廣泛關注.

關于聲子晶體的拓撲優化方法主要有兩類:梯度類算法和非梯度類算法.梯度類算法通過處理目標函數和約束條件的靈敏度信息,使可行解沿問題的梯度方向移動更新,可通過較少迭代次數完成問題的求解.目前,該方面的研究主要采用移動漸近線算法(method of moving asymptotes,MMA)[3-6]和雙向進化結構優化(bi-directional evolutionary structural optimization,BESO)方法[7-8]對聲子晶體進行帶隙優化.雖然梯度類算法求解效率高,但存在容易陷入局部最優解的弊端.非梯度類算法是通過對適應度函數的不斷評估和多次迭代,在設計域內找到滿足要求的近似全局最優解.其中,遺傳算法在聲子晶體拓撲優化中應用最為廣泛,如標準遺傳算法(standard genetic algorithm,SGA)[9-15]、自適應遺傳算法(adaptive genetic algorithm,AGA)[16-18]、多精英遺傳算法(multiple elitist genetic algorithm,MEGA)[19]、快速非支配排序遺傳算法(non-dominated sorting genetic algorithm II,NSGA-II)[20-24]與不同的數值計算方法相結合皆可實現對聲子晶體的帶隙優化.

雖然聲子晶體的拓撲優化近年來已經取得了一系列的重大進展,但優化中存在的“棋盤格效應”使得結構中往往出現某種孤立材料單元、超細而導致剛度較弱的鏈接單元、微結構組元幾何形貌十分復雜、微結構邊界和各組元邊界或界面粗糙等異常情況,這對于實際加工制造十分不利,因而大大限制了聲子晶體的應用.針對優化結果常常出現孤立材料單元這一問題,本文開展考慮可制造性的二維多相聲子晶體拓撲優化研究,通過對優化過程中拓撲構型的連通性分析以限制最小材料連通域的占比,避免出現孤立材料單元,進而建立提高其可制造性的優化模型.在此基礎上,應用NSGA-II 對模型進行求解,得到了多目標帕累托(Pareto)解集.具有代表性的數值算例結果表明本文方法可在滿足帶隙性能等指標的前提下同時兼顧到制造可行性的要求,有利于聲子晶體優化結果的實際加工制造和推廣應用.

1 考慮可制造性的聲子晶體拓撲優化方法

1.1 帶隙計算理論模型

主要考慮二維聲子晶體可視為非均質各向同性彈性材料的情況,忽略體力影響時,其彈性波動方程為

式中,ρ,λ 和 μ 分別是質量密度和拉梅常數,u是位移向量.? 和 ω 分別為梯度算子和角頻率.假定彈性波僅在xy平面內傳播,則有 ?u/?z=0,于是方程(1)可分解為以下兩個解耦的二維波動方程

僅考慮由式(2)描述的面內波動情形而不考慮由式(3)描述的面外或反平面波的傳播.根據周期性結構中的Bloch 定理,位移向量u(r) 可用如下公式表示

式中,uk(r) 為周期性函數,k=k(x,y) 為波矢.利用有限元方法求解,式(2)可轉化為以下特征方程

式中,K和M分別為整體剛度矩陣和整體質量矩陣,U是節點位移向量.給定波矢求解式(5)中的特征頻率可得到能帶結構或頻散曲線.為了方便起見,能帶結構圖的縱坐標為帶隙歸一化頻率 Ω=ωa/(2πct),其中a為晶格常數,ct為基體的剪切波波速.

1.2 優化模型

實施聲子晶體的拓撲優化,主要思路是采用有限元網格對單胞結構進行離散化,通過優化算法在網格中填充合適的材料,得到滿足優化目標要求的拓撲構型.由于遺傳算法采用二進制變量,特別適合建立兩相材料與拓撲構型之間的關系[12,19,24],即采用0 和1 分別代表兩種不同的材料,故聲子晶體的拓撲優化多采用遺傳算法.對于多相聲子晶體,可通過引入“二進制向量”的方式對離散單元的材料進行編碼,即

式中,φi(i=1,2,···,NE) 為任意離散單元的二進制編碼(NE為單元數),bs為二進制編碼值,t為向量長度,NM為材料數.以4 相材料為例,有

對于更多相聲子晶體,可通過增加向量 φi的長度t來進行相應的編碼.由離散單元的“二進制向量”編碼組合可得到整個設計域的材料分布，即

式中，? 為設計域的材料分布.在優化問題求解中,采用N×N的網格對單胞進行離散,由于結構的平移周期性和點群對稱性,通過波矢遍歷不可約布里淵區邊界求解特征頻率可得到能帶結構,則設計域內單元數目減少至 (N/2)(1+N/2)/2,對應遺傳算法染色體長度為t(N/2)(1+N/2)/2 .

為了避免出現孤立材料單元而難以或無法實際制造的問題,需要引入合適的可制造性約束條件.本文試圖通過控制連通單元的占比達到限制孤立單元的目的.引入圖像處理中四連通區域的定義,即從區域內一點出發,可通過上下左右4 個方向的移動組合,在不越出區域的前提下,能達到區域內任一像素的區域定義為四連通區域.將離散后的網格看作像素點,則可將此定義擴展到材料連通單元中,進一步由連通單元組成材料連通域,如圖1 中所示.紅、綠、黑不同標記的網格分別由3 種不同材料填充,其中,紅色、綠色區域分別為獨立的連通域,而黑色區域雖在角點處相連接,但屬于兩個連通域.

圖1 連通域示意圖Fig.1 The connectivity domain

按照編碼方式為元胞網格填充材料后,可進行連通域的識別和其相對于元胞面積的占比計算,即

式中Ajk(?) 和Ijk(?) 分別為第k種材料的第j個連通域的面積和其相對于元胞的面積占比;A為元胞面積,NM為材料數,NQ為每種材料所包含的四連通域數目.Ijk(?) 過小會使得聲子晶體的可制造性下降,故通過限制Ijk(?) 值保障其可制造性,相應的附加約束條件具體如下

式中,I(?) 為Ijk(?) 的最小值,r*為閾值,可根據實際對可制造性的要求和數值試驗確定,r*越大則表示對結構可制造性要求越高.

對于聲子晶體的帶隙優化研究,常以最大化相對帶隙寬度、最大化絕對帶隙寬度、特定帶隙中心頻率等為優化目標.本文主要以輕質結構的減振隔振為應用背景,考慮到實際應用中通常是針對特定場合對特定頻段的減振需求設計具有預期帶隙的聲子晶體,同時,往往還需兼顧結構輕量化設計的需求,故本文以在特定頻率段帶隙最寬和結構質量最小為優化目標.

綜上所述,考慮可制造性約束的二維多相聲子晶體多目標優化優化模型表述如下

式中,? 為設計域材料分布;Δ ωn為第n帶隙的帶隙寬度,為了方便描述,指定第n條能帶和第n+1 能帶之間存在的帶隙為第n帶隙,即當第n+1 條能帶的最小值大于第n條能帶的最大值時,帶隙寬度Δωn＞0,反之帶隙寬度 Δ ωn=0,帶隙不存在;F1(?)為特定頻率范圍段內帶隙寬度之和的占比,其值越大表示帶隙對特定頻率段的覆蓋率越高,ωlow和ωupp分別為特定頻率范圍的上、下界,F2(?) 為優化結構質量,N為網格的離散規模,ρs為單元質量密度,V為單元面積.

1.3 優化算法

對于式(9)的多目標優化問題,可采用快速非支配排序遺傳算法 N SGA-II[25]進行求解,該方法能夠有效生成 P areto 解集,使決策者能夠根據需要從Pareto解集中選擇出最適合的解.算法流程如圖2 所示.

圖2 基于NSGA-II 算法的聲子晶體多目標優化流程圖Fig.2 Flow chart of the multi-objective optimization of phononic crystals based on the NSGA-II algorithm

在計算種群內個體的適應度值時,將種群中不滿足可制造性約束和雖滿足可制造性約束但無帶隙個體的帶隙特性目標函數值賦予極小值0.000 1,可以增加滿足可制造約束且存在帶隙個體在非支配排序時處于更優的非支配排序層的可能性,便于后續選擇操作時對此類優秀個體的保留.

2 數值算例

以二維正方晶格聲子晶體面內模態為例,晶格常數為0.01 m,聲子晶體單胞由4 種材料組成,即金(綠)、環氧樹脂(灰)、硅橡膠(黃)和硫化橡膠(藍),材料名稱、編碼 φi、質量密度 ρ、楊氏模量E和泊松比 ν 如表1 所示.以在給定歸一化頻段0～0.05 帶隙最寬為優化目標,有 ωlow=0,ωupp=0.05 ;考慮可制造性約束要求,取r*=0.005 ;對單胞進行網格劃分,取N=16,則設計域內的單元數目和染色體長度分別為36 和72;遺傳算法的參數選取為:種群規模20,交叉率0.9,變異率0.02,采用二元錦標賽選擇方式,單點交叉方式,設定迭代次數1000 為終止條件.

表1 材料屬性Table 1 Properties of the materials

種群規模Np取值對優化結果不會產生影響,但是會對優化收斂過程產生影響,即Np取值越大,所需得到Pareto 解集的迭代次數越少.對于本文考慮的問題而言,因在每一步迭代中要對種群中所有個體的帶隙進行有限元計算,增大種群規模會導致單次迭代過程中有限元帶隙計算次數大大增加,加大單次迭代的計算成本.因此,綜合考慮以上因素并通過數值試驗,本算例選取Np=20 .為了提高算法的收斂速度,在初始種群中引入具有帶隙的“種子”結構,其余個體隨機產生,“種子”單胞構型如圖3 所示.

圖3 “種子”結構Fig.3 The seed unit-cell

圖4 為進化 1 000 代的Pareto 解集,橫、縱坐標分別為目標函數F1和F2的值,可以看出兩個目標之間的制約關系,即特定頻段的帶隙寬度之和的增大將使得結構質量也在同時變大.

從圖4 中選擇4 個非支配解A,B,C和D進行分析,其拓撲構型、能帶結構和透射譜如圖5 所示.由圖5(a)可以發現,非支配解A,B,C和D的拓撲構型均由金、環氧樹脂和硫化橡膠組成且具有相似的特征,即由金和環氧樹脂組成的夾雜體被硫化橡膠包覆層包裹后嵌入環氧樹脂基體中.圖5(b)與圖5(c)的對比表明能帶結構的帶隙范圍與透射譜頻率衰減范圍吻合.由圖5(b)的能帶結構可知,解A,B,C和D的第3 帶隙充分打開,其歸一化帶寬分別為0.017 1,0.022 0,0.025 2 和0.030 0,計算得到的對應目標函數值如表2 所示.

圖4 附加可制造性約束的多目標Pareto 解集Fig.4 Pareto optimal solutions from the MOOP considering the manufacturing constraint

由表2 可知,結構A,B,C的帶隙特定頻段的覆蓋率比D分別小42%,26% 和15%,結構質量比D分別小52%,44%和35%,因此,決策者可根據實際應用情形的目標權重從Pareto 解集中挑選自己所需要的聲子晶體結構.

表2 附加可制造性約束優化結構目標函數值,S 為單目標優化結果Table 2 The objective function values of the optimized structures from the MOOP and the SOOP considering the manufacturing constraint

為進一步理解優化結果的帶隙形成機理,對圖5(b1～ b4)帶隙的上、下邊界所對應的振動模態進行分析.現以結構D為例進行說明(構型A,B,C的情況與之類似).構型D的第3 個帶隙上、下邊界處對應的振動模態如圖6 所示.

圖5 附加可制造性約束的多目標優化結果Fig.5 Pareto optimal solutions from the MOOP considering the manufacturing constraint

由圖6 可以看出優化結果的帶隙機理為局域共振型.其振動特性可近似的通過簡化的等效“彈簧-質量”系統進行描述,相應的帶隙上、下邊界的頻率可由下式進行估算[26-27]

圖6 結構D 的帶隙邊界所對應的振動模態Fig.6 Vibration modes corresponding to the band-gap edges of the optimized structure D

其中,flow和fupp分別為帶隙的下邊界和上邊界頻率,K為包覆層(彈簧)的動態等效剛度,M1和M2分別為金散射體和環氧樹脂基體的動態等效質量.

由圖5 中的優化結果構型可以看出,優化構型中僅剩下3 種材料,材料初始分布中的硅橡膠在優化過程中被淘汰掉,包覆層全部由彈性模量更大的硫化橡膠構成,使得彈性包覆層動態等效剛度增加.由式(12)可知,當其他參數不變時,K的增大將導致帶隙起始(下邊界)和截止(上邊界)頻率均會出現不同程度的上升.由于截止(上邊界)頻率上升幅度更大,因而帶隙范圍變寬,帶隙對特定頻率段的覆蓋率隨之提高.

另外,通過對比圖5(a),4 種構型發現,金散射體的占比不斷提高,即M1逐步增大,而K與M2基本保持不變.由式(12)可知,M1的增大將導致帶隙起始(下邊界)和截止(上邊界)頻率降低,而起始(下邊界)頻率相比截止(上邊界頻率降幅更大,因而帶隙范圍變寬,帶隙對特定頻率段的覆蓋率隨之提高.即F1與F2同時增大,這種制約關系與圖4 十分吻合.

為了說明多目標和單目標優化結果的差異,本文對只考慮了目標函數F1的單目標帶隙優化問題也進行了求解,圖7 給出了單目標進化曲線、近似最優解S的拓撲構型、能帶結構和透射譜.

從圖7(a)中可明顯看出拓撲構型隨迭代進程的演化過程,即隨著進化代數的增加,包覆層的硅橡膠(黃)分布減小而硫化橡膠(藍)分布增加,金散射體分布變化不大;進化800 代左右后種群的最佳適應度收斂到穩定值,表明種群進化趨于完善,對應的個體即為近似最優解S,其拓撲構型(圖7(b))由硫化橡膠包覆著金嵌入環氧樹脂基體中,能帶結構(見圖7(c))顯示可在第3 和第4 能帶間打開歸一化帶寬為0.030 2 的帶隙,由此計算得到結構S的目標函數值同樣在表2 中給出以作對比.對比D和S的結果發現,兩者對特定頻段的覆蓋率十分接近,雖然S比D的特定頻段覆蓋率高2%,但其質量卻比D大18%,所以從結構整體的角度來講,結構D比S更優,這說明多目標優化不僅可以得到和單目標優化結果接近的非支配解,還可以得到其他的多目標優化非支配解集,因此較單目標優化更具有優越性.

圖7 單目標優化結果Fig.7 The optimized solutions from the SOOP

此外,為說明考慮可制造性優化的有效性和優勢,圖8 給出了未附加可制造性約束的多目標優化Pareto 解集.從圖8 內選擇與考慮可制造性約束的4 個非支配解的目標函數值具有可比擬性的解A′,B′,C′和D′分析,其拓撲構型分別如圖9 所示.

圖8 未附加可制造性約束多目標Pareto 解集Fig.8 Pareto optimal solutions from the MOOP without the manufacturing constraint

從圖9 可以看出,未考慮可制造性約束的結果構型中A′和B′存在多個金(綠)、硅橡膠(黃)和環氧樹脂(灰)孤立單元,C′和D′對角線位置存在多個環氧樹脂(灰)孤立單元.對4 種結構單胞內的孤立材料單元進行統計,結果如表3 所示.

表3 單胞的孤立材料單元數Table 3 The number of the isolated elements in the unit-cells

圖9 未附加可制造性約束多目標優化結果Fig.9 The optimized solutions from the MOOP without the manufacturing constraint

由以上結果分析可見,未考慮可制造性約束的多目標優化結果均存在孤立材料單元,而圖5 中附加可制造性約束的優化結果則有效地避免了這種情況的發生,表明本文所提出的考慮可制造性約束的聲子晶體多目標拓撲優化方法可在滿足帶隙等性能指標的前提下同時兼顧到實際制造加工的可行性,有效提高了優化結果的實用性.

3 結論

本文針對聲子晶體拓撲優化結果常會出現制造性能欠佳的問題,通過引入可制造性附加約束,建立了考慮可制造性約束的二維多相聲子晶體多目標拓撲優化模型,基于NSGA-II 和有限元法對二維四相正方晶格聲子晶體面內模態進行多目標拓撲優化設計,得到了多目標Pareto 解集,數值算例表明了本文模型的有效性,所得主要結論如下.

(1) 通過對比引入可制造性附加約束的多目標優化與相應的單目標拓撲優化結果發現,Pareto 解集內可以找到與單目標優化結果接近的非支配解;但多目標的優化解在一定程度上可優于單目標優化近似最優解,表明多目標優化不僅可以得到和單目標優化結果近似非支配解,還可以得到其他的多目標非支配解集.

(2) 通過對比引入和不考慮可制造性附加約束的多目標優化結果發現,引入可制造性附加約束可有效避免通常優化結果中存在孤立材料單元的情況,其優化結果在滿足帶隙性能指標的同時可兼顧到制造加工的可行性,更具有實用價值.

本文算例針對二維正方晶格,驗證了協同考慮帶隙性能和可制造性約束的多目標優化方法的可行性和有效性,但該方法對其他類型的晶格,比如二維三角晶格和六角蜂窩晶格、甚至三維體心立方晶格等等,也同樣適用.其中,單胞離散及材料編碼方案的選擇、可制造性約束的添加和優化算法的選擇可根據具體問題進行確定.對于三維晶格而言,材料初始分布構型的選取、復雜三維網格的劃分、計算規模大、優化效率低等仍是其拓撲優化問題的主要挑戰[28],同時也是該方面研究工作未來需要加強的一個方向.

最后需要強調指出的是,影響聲子晶體優化結果可制造性的因素很多,比如孤立材料單元問題,超細而導致剛度超弱鏈接單元問題,不同材料連通域的個數問題,微結構奇形怪狀問題和微結構邊界或組元界面非光滑問題等.本文提出的多目標優化模型主要以避免出現孤立材料單元而提高多相聲子晶體優化結構的可制造性能,其他方法比如避免出現超細超弱鏈接單元[20,29-30],力求微結構邊界或組元邊界具有一定的光滑度[31]等也是提高其可制造性的有效措施,值得引起科研和實用工作者的高度重視.