彎曲波寬頻分波超柵拓撲優化設計和表征1)

李林 張雪彬 劉濤 章俊

(重慶大學航空航天學院,重慶 400044)

引言

自由地控制波的傳播一直是波動調控領域中的熱門話題,如利用透鏡[1]、光纖[2]、消聲器[3]等對光波/聲波的調控.作為一種亞波長準二維超材料,超表面[4]的出現使得波的各種反常、新穎調控變得可能,最近十幾年得到了眾多科研工作者的關注.超表面研究始于電磁波領域[5],隨后被引入到聲波[6-7]和固體彈性波領域[8-9].作為一種經典的彈性波,彎曲波在工程中十分常見,是梁板類工程結構在橫向動態載荷激勵下產生的一種波動形式.利用超表面/超材料對彎曲波進行調控在梁板類結構隔振減振[10-11]、基于梁板類結構開發振動能量收集裝置[12-15]和各種傳感器以及結構健康檢測等領域有重要應用.截至目前,科研工作者利用超表面/超材料已經實現了彎曲波的聚焦[16-21]、反常折射和反射[22-25]、隱身[26-35]、完美吸收[36]、分波[37-38]、非對稱傳輸[39-40]和拓撲保護[41-42]等調控.

通過設計結構或材料某種參數可改變的微結構來構成超表面的子功能單元是目前常用的手段,利用這一參數的改變來調控彈性波在子功能單元中的等效相速度,從而實現超表面中每個子功能單元的透射/反射波相位滿足一定的空間分布,達到波前調控目的.截止目前,基于上述思想設計的彎曲波超表面子功能單元主要有以下幾類: 組合梁型[23]、附著臺柱直梁型[24,43]、變厚度梁型[25]和鋸齒形聲子晶體梁型[21,44].組合梁型子功能單元由兩部分不同材料的直梁軸向拼接而成,通過改變這兩部分直梁的長度比來調節彎曲波在其中的傳播速度.附著臺柱直梁型子功能單元由直梁上附著不同幾何尺寸的柱狀體構成,附著柱狀體幾何尺寸的不同會改變彎曲波有效相速度.變厚度梁型子功能基于梁中彎曲波相速度與梁厚度相關這一事實而提出.鋸齒形聲子晶體梁型子功能單元為波紋型梁,通過改變波紋狀輪廓曲線的高度來調節梁的等效彎曲剛度,從而調節彎曲波有效相速度.有關超表面的更多細節,讀者可以參考北京交通大學陳阿麗教授課題組最近發表的綜述文章[4].

隨著對超表面研究的深入,研究人員發現當入射角大于某一臨界角時,用來預測折射型超表面行為的廣義斯涅耳定理(GSL)不再適用.因此,表達式如下的高階衍射理論[45-47]被提了出來,即

其中kt,r和ki分別為透射/反射波和入射波波數,G=2π/Γ,Γ為超表面中單個超胞在其排布方向的寬度.一般情況下,單個超胞由多個子功能單元(subunit)組成,且這些子功能單元的透射/反射波相位覆蓋[0,2π]區間且滿足一定的空間分布規律.n表示衍射波階次,為第n階衍射波角度.在本文中入射波和衍射波具有相同的波數,將ki=2π/λi代入式(1)可以得到如下關系

其中 η=λi/Γ,λi為入射波波長.當n=1 時,高階衍射定理即退化為GSL.與GSL 相比,高階衍射定理除了可以預測第1 階衍射波,還可以預測其他階次的衍射波,因此高階衍射定理可被用來設計分波超表面.基于高階衍射定理設計的超表面又稱為超柵(metagratings),因其具有類似光柵的周期結構.根據式(2)中 η值,基于高階衍射定理設計的超柵可以分為三類: η >2.0,1.0<η<2.0和η <1.0.由式(2) 可知,第一類超柵中只存在0 階衍射模式,而且與入射角無關,因此這一類超柵可以實現全角度反射[8,10].第二類超柵中存在0 階和±1 階模式,但由于0 階在子功能單元中傳播次數最多[36],很容易被阻尼耗散掉,因此0 階可以忽略,+1 階對應透射波模式,?1 階模式為透射波還是反射波取決于超胞中子功能單元的數目,如果子功能單元的數目為奇數,為透射波,反之為反射波.因此,這類超柵可用來實現非對稱傳輸調控[48].第三類超柵中的衍射模式比前兩類超柵都多,因此這類超柵更容易實現分波功能[38,49].

不論是超表面還是超柵,目前的設計大多基于理論直覺、實踐經驗和兩者的結合,這樣設計出的超表面和超柵調控性能往往不是最優.此外,此類方法往往只針對單一頻率進行設計,導致設計出來的超表面和超柵工作頻率范圍非常窄,嚴重限制它們在實際工程中的應用.為了解決這個問題,目前提出了兩種辦法: 一種是引入主動控制技術[50];另一種為實施拓撲優化設計.主動調控技術通常需要引入額外的結構或裝置,這將增大超表面的復雜性.作為一種被動寬頻調控技術,拓撲優化設計[51]將數學優化算法引入到超表面的設計中,通過對子功能單元微結構拓撲優化設計來滿足透射/反射波相位和幅值要求.例如,Rong 等[52-53]建立了基于多目標遺傳算法NSGA-II 的彈性波超表面拓撲優化設計框架,通過在目標函數中設置多個目標相移,優化得到了不同頻率下具有不同相位的子功能單元,提出了一種在不同頻率有不同調控功能的彈性波超表面.Li等[38]基于遺傳算法設計了一類超胞由兩個子功能單元組成的結構簡單超表面,實現了單一頻率下彎曲波的分波和聚焦調控.目前,基于拓撲優化設計超表面得到了越來越多的關注,但寬頻超表面的拓撲優化設計工作還較少.拓撲優化所得到的結構形狀一般比較復雜,加工是個難題,慶幸的是隨著3D打印技術的發展,這類結構的加工難度在逐漸降低.因此,結合拓撲優化和3D 打印技術開展彈性波寬頻超表面研究變得可行.

本文基于拓撲優化方法提出一種彎曲波寬頻分波超柵設計方法.具體地,利用多目標算法NSGAII 搭建一套二維平面應力穩態動力學問題拓撲優化框架,然后利用這套優化平臺優化超柵子功能單元,最后對優化得到的超柵進行數值表征、利用3D 打印技術制備試件開展實驗驗證,并與其他設計方法進行比較,進一步驗證本文提出的設計方法的優勢.

1 基于高階衍射定理的彎曲波寬頻分波超柵拓撲優化設計和表征

在前面的引言中已經提到,當 η ≥2 時,n≡0,此時衍射波場中只有0 階通道存在.隨著η的減小,衍射階次會增多.圖1 所示為根據式(2)計算得到的兩個不同 η值情況下衍射角與入射角之間的關系.從圖中可知,當 η <1時,存在多個衍射階次,且不同階次的衍射波有不同的傳播方向.因此,可以設計η值小于1.0 的超柵來實現分波功能.本文考慮設計一種彎曲波寬頻分波超柵.

圖1 不同衍射階次的衍射角與入射角的關系Fig.1 Relationship between diffraction angles of different diffraction orders and incident angles

為了制備簡單起見,本文設計的分波超柵每個超胞由2 個子功能單元組成,且第一個子功能單元選為與主板等厚度的直梁,其尺寸為長度L=50 mm,沿y方向的寬度W=34.36 mm,厚度h=2 mm.第二個子功能單元由優化得到,優化目標為在設定頻率范圍內由這兩個子功能單元所組成的分波超柵功能保持不變.本文設定頻率范圍為: 中心頻率為 ωc=2.4 kHz,上下限在中心頻率基礎上各偏移10%,即[2160,2640] Hz.圖2 所示為第二個子功能單元的二維平面應力優化模型,其中左右兩邊灰色區域代表與超柵相連接的主板,中間黃色區域為第二個子功能單元的優化區域.此優化區域由5 個胞元周期排列構成,優化區域的總長度為L,總高度為H.考慮到拓撲優化可能得到的結構復雜性,灰色區域和黃色區域采用同一種材料.本文具體采用的是3D 打印材料VeroPureWhite,這樣即使優化得到的結構較復雜,也可通過現有的3D 打印技術進行制備.VeroPureWhite 材料的實測楊氏彈性模量E=3.2394 GPa,密度ρ=1185 kg/m3,泊松比 ν=0.4185[21].在有限元模擬中,沿著z方向的單位位移加載在圖中的黑色虛線上.為了減少反射波,模型的左右兩端設定為完美匹配層(PML).本文考慮了三種不同的優化區域總高度H,目的是想研究不同H情況下能否得到工作機制不同的胞元結構.子功能單元確定之后,按照圖3 排列即可形成超柵整體結構.超胞中相鄰子功能單元間留有間隙,本文中間隙設為1 mm.

圖2 超柵子功能單元二維拓撲優化模型示意圖Fig.2 Schematic of topology optimization 2D models for the subunits in a supercell

圖3 薄板彎曲波超柵整體結構俯視圖,本文中子功能單元數為2Fig.3 Top view of metagratings for flexural waves in thin plates.In this paper,the number of subunits is 2

1.1 定義適應度函數

優化胞元區域被離散成一定密度的像素塊集合,每個像素塊的質量密度都是一個獨立的設計變量,并通過適應度函數最小化來確定像素塊的密度分布.在本文中,適應度函數定義如下

其中,S為頻率區間[2160,2640] Hz 內的采樣頻率點數.考慮到數值計算成本,本文中S=3,對應的頻率分別為頻率區間的下限、中心頻率和頻率區間的上限.其中和是第j個采樣頻率下的目標相移和目標透射率;?j和Tj是優化過程中第j個采樣頻率下的可行相移和可行透射率.可行相移 ?j=?s·k0,?s為經過第一和第二兩個子功能單元的透射波相應波峰或波谷的位置差,Tj等于透射波幅值之比.由于本文研究的是透射分波超柵,所以目標透射率=1(j=1,2,···,S).不同頻率下,目標相移是不一樣的,具體由以下公式確定

式(3)中的vios為罰函數.在本文中,它由兩部分組成.第一部分為

其中nd,Ad分別為非連通區域的數量和非承載區域的像素數量.f1,f2為相應的罰參數.在本文中,f1=0.5,f2=0.1.罰參數滿足f1>f2以便有效地促進結構的連通性.這一部分罰函數的作用是對不可行的拓撲結構進行懲罰,比如斷開的拓撲形狀、不能提供承載作用的像素[54-55].第二部分為

其中ft=0.1 是罰參數,nt為薄弱連接區域的個數.在本文中,通過判斷沿x和y方向單列連續的拓撲結構的長度是否小于預設的最小加工尺寸來確定薄弱連接.這一部分的作用是為了提高拓撲結構的可制造性和可靠性.除了這些罰函數外,在計算適應度函數之前,還應對拓撲形狀中的鉸接情況進行處理以提高優化所得結構的承載能力和可加工性.

本文采用了“粗?細”策略的遺傳算法[56],在第一階段,優化區域x方向被劃分為20 個網格,在第二階段則被劃分為 40 個網格.y方向的網格密度與x方向始終保持一致.在這兩個階段,染色體種群的個體數均為Np=100,染色體交叉概率為Pc=0.9,變異概率為Pm=0.02,總的迭代次數設定為1000 次.

具體的拓撲優化流程與已有文獻報道中的流程基本一致[57],但是本文在子代的選擇方面做了一些改變,并沒有完全采用精英化策略[54],而是在輪盤賭和精英化策略之間采用概率選擇.具體地,在本文中,輪盤賭和精英化選擇的概率各為50%.這樣處理的好處是可以避免遺傳算法過早地收斂.

1.2 胞元優化結果

圖4 所示為在中心頻率2.4 kHz 下,優化區域高度為H=3h情況下的胞元拓撲優化結構進化歷程,從此圖可以看到,優化從隨機生成的連接性較差的結構開始,經過一百代左右優化之后,拓撲結構的連通性明顯改善.之后,遺傳算法搜索到有效的拓撲特征并慢慢演化結構,結構連通性在整個演化過程中逐漸提高.

圖4 優化區域高度H=3h,在2.4 kHz 中心頻率激勵下的胞元拓撲結構的進化歷程Fig.4 Evolutionary history of topological structure of a unit cell as the optimization region height H=3h and under excitation of a 2.4 kHz loading

圖5 所示為三種不同優化區域高度下胞元的最終拓撲優化結構.圖6 所示為以圖5 中三種拓撲結構為胞元分別形成的三種一維聲子晶體梁的能帶圖,其中最低的黑色擴散曲線對應彎曲波模式.擴散曲線結果顯示這三種聲子晶體梁的工作機理相同,均為通過改變等效彎曲剛度從而實現有效相速度的改變.圖7 所示為由圖5 中三種優化結構分別組成的子功能單元在2.4 kHz 中心頻率激勵下的透射模擬結果.圖8 所示為由這三種拓撲結構所構成的子功能單元在三個采樣點頻率下的透射率和相移與目標值的比較.在中心頻率下,每個胞元準確地實現了所需的相移;在其他兩個頻率點,實際相移與目標相移存在小的偏差,這樣的結果符合預期;透射率均大于0.9,表明所設計的透射超柵將具有較高的透射率.總的來說,H=2h情況下優化得到的胞元(見圖5(a))性能最優.

圖5 三種不同優化區域高度下所得到的胞元最終拓撲優化結構Fig.5 The optimized topological structures of unit cells under three different optimization region heights

圖6 以圖5 中三種拓撲結構為胞元分別形成的三種一維聲子晶體梁的能帶圖Fig.6 Band structures of the three 1D phononic crystals constituted by the three units presented in Fig.5

圖6 以圖5 中三種拓撲結構為胞元分別形成的三種一維聲子晶體梁的能帶圖(續)Fig.6 Band structures of the three 1D phononic crystals constituted by the three units presented in Fig.5 (continued)

圖7 三種拓撲優化結構分別組成的子功能單元在中心頻率激勵下透射波行為數值模擬表征結果Fig.7 Numerical characterization of wave transmission behavior through three subunits composed by each of the three optimized structures under the central frequency excitation

圖8 分別由三種拓撲優化結構構成的三種子功能單元在三個采樣點頻率下的(a)相移和(b)透射率Fig.8 (a) Phase shifts and (b) transmittances of subunits composed by the three optimized structures at three different frequencies

1.3 超柵整體結構數值模擬表征

基于圖5 中優化得到的胞元,形成第二個子功能單元(見圖7),然后將超胞(第一和第二個子功能單元)沿y方向周期排列即可形成超柵(見圖3).考慮到數值模擬和實際加工成本以及加工技術限制,本文中的超柵由四個超胞組成,胞元沿y方向的寬度為34.36 mm.因此,中心頻率下超柵的 η=0.707.圖9 所示為由圖5 中三種胞元分別形成的三種超柵整體結構.為了后面表述方便,這三種超柵分別被稱為超柵1、超柵2 和超柵3.為了驗證這三種分波超柵的性能,建立如圖10 所示有限元數值模型,彎曲波由線源激發,單位載荷沿z方向加載.圖11 展示了這三種分波超柵在三個采樣點頻率激勵下的響應結果.從圖中可以清晰地看出,正入射的彎曲波在經過超柵之后被分成了兩束,傳播角度與高階衍射理論預期角度(圖中標注值)非常吻合,且在三個頻率下的角度基本一致,說明超柵在本文關注的頻率范圍內功能保持不變.

圖9 由圖5 中三種胞元分別形成的三種超柵整體結構Fig.9 Three metagratings constituted by each of the three unit cells presented in Fig.5

圖10 超柵結構有限元數值模擬模型Fig.10 Numerical simulation models of metagratings with the finite element method

圖11 三種超柵在三個采樣點頻率激勵下響應的數值模擬結果Fig.11 Numerical simulation results for the three metagratings under excitation of three different frequencies

1.4 實驗驗證

為了進一步驗證本文基于優化設計所得到的寬頻超柵性能,本文利用3D 打印技術制備了兩種試件.一種是與圖9(b)所對應的超柵整體結構,另外一種是基于優化所得到的單個子功能單元,如圖12 所示.實驗測試平臺如圖13 所示,主要設備包含一套激光掃描測振儀系統(Polytec PSV-600)、一個功率放大器、一個相位轉換裝置(自制)和一臺信號發生器(ATA-308,Agitek).為了產生彎曲波,在試件正反面相同位置各粘一排壓電片(無錫海鷹).由信號發生器產生的電壓信號經過功率放大器放大,然后經過相位轉換裝置出來電壓相反的兩路信號.隨后,這兩路信號分別輸入到正反面的壓電片上,在相反的兩個電壓作用下,上下表面的壓電片將產生反向的面內變形,相當于對試件施加一諧變彎矩,從而產生彎曲波.離面質點振動速度由激光測振儀測量.為了減少邊界的反射波,試件的邊界被粘上一層阻尼材料(Blu-Tack).整個試件使用彈性繩懸掛在框架上.

圖12 實驗試件Fig.12 Experiment specimens

圖12 實驗試件(續)Fig.12 Experiment specimens (continued)

圖13 本研究實驗測試平臺Fig.13 Experimental setup in this study

圖14～ 圖16 所示為超柵整體結構在三個采樣點頻率激勵下的實驗結果.結果表明,在三個頻率下,入射波經過超柵之后均被分成兩束,這兩束透射波相對入射波的折射角度為 ±45°,與理論值和前面的數值模擬結果均一致,表明本文基于拓撲優化所設計的超柵確實在設定的寬頻范圍內具有所設計的功能.

圖14 激勵頻率為2160 Hz 時,測得的三個時刻的超柵離面速度場Fig.14 Snapshots of measured out-of-plane velocity fields under the excitation of 2160 Hz at three instants

圖15 激勵頻率為2400 Hz 時,測得的三個時刻的超柵離面速度場Fig.15 Snapshots of measured out-of-plane velocity fields under the excitation of 2400 Hz at three instants

圖16 激勵頻率為2640 Hz 時,測得的三個時刻的超柵離面速度場Fig.16 Snapshots of measured out-of-plane velocity fields under the excitation of 2640 Hz at three instants

除了上面的超柵整體結構,本文還測試了單個子功能單元的結果.圖17 所示為第一個子功能單元和由圖5(b)胞元所構成的第二個子功能單元在三個不同頻率激勵下的實驗結果.從圖中可以清楚地看到經過它們之后的透射波相移與理論預期非常吻合,進一步驗證了本文基于拓撲優化方法設計寬頻超柵的可行性.

圖17 三個采樣點頻率激勵下,拓撲優化梁與直梁實測相位比較Fig.17 Comparison of phases of transmitted waves between optimized beams and straight beams at three different frequencies

2 與其他方法設計的超柵的寬頻優勢比較

為了進一步展示本文拓撲優化設計超柵的寬頻特性,把本文的方法與已有的兩種彎曲波超柵設計方法進行了對比.第一種方法為前面引言中提及的鋸齒形聲子晶體梁方法;另一種方法為變厚度梁方法[10,25].基于這兩種方法本文設計了相應的超柵結構,分別稱為鋸齒型超柵和直梁型超柵.然后對超柵整體結構進行數值表征,并與本文的超柵進行對比.兩種超柵的超胞均由兩個不同高度的子功能單元構成,根據相位要求,兩直梁高度分別為2 mm 和6.4 mm,兩個鋸齒形聲子晶體梁的高度分別選為1h和3.9h,其他參數與本文的超柵一致.圖18 和圖19 分別展示了直梁型和鋸齒型分波超柵在[2160,2640] Hz 范圍內的數值模擬結果.可以看出,直梁型分波超柵的性能較差,鋸齒型分波超柵在前兩個頻率激勵下的性能也不太理想,兩者均無法在本文所考慮的寬頻范圍內實現穩定的分波功能.本節內容表明本文提出的基于拓撲優化的寬頻超柵設計方法非常有必要,為寬頻超柵設計提供了一種可行路徑.

圖18 直梁型分波超柵數值模擬結果Fig.18 Numerical simulation results of straight-beam-type wave-splitting metagratings

圖19 鋸齒型分波超柵數值模擬結果Fig.19 Numerical simulation results of zigzag-type wave-splitting metagratings

3 結論

本文提出了一種基于拓撲優化的彎曲波寬頻超柵設計方法.利用遺傳優化算法設計一系列滿足高階衍射理論的亞波長胞元結構,這些胞元組成的超胞產生的相位能夠覆蓋[0,2π]范圍,并在目標寬頻范圍內具有高透射率.然后將這些超胞組合以形成具有特定相位梯度的超柵,并通過數值模擬和實驗驗證了超柵的分波功能和寬頻特性.此外,通過與直梁型分波超柵和鋸齒型分波超柵進行對比,發現本文的超柵寬頻性能明顯優于這兩種超柵,顯示出本文提出的基于拓撲優化的寬頻超柵設計方法的必要性和有效性.本文提出的方法同樣也適用于其他寬頻超柵的設計.