聲子晶體板復能帶計算方法

陳圣兵　張浩　宋玉寶

摘要： 聲子晶體具有彈性波帶隙，可以用于結構振動與噪聲控制。聲子晶體傳統能帶算法一般給定波矢k在不可約布里淵區邊界取值，然后求解特征頻率ω，得到ω-k曲線。因而，傳統方法中波矢k只取實數，只能求解實能帶。為了求解復能帶，一般需要給定頻率ω，求解特征波矢k，從而得到k-ω曲線。提出了一種參數變換方法，解決了特征波矢求解中復雜的非線性特征值問題，實現了復能帶的快速求解。最后，采用兩個算例對文中算法進行了驗證，包括布拉格聲子晶體板和局域共振聲子晶體板，研究了帶隙內衰減常數隨波傳播方向的變化和阻尼對帶隙的影響。

關鍵詞： 聲子晶體; 超材料; 復能帶; 參數變換

中圖分類號： O735; TB535? 文獻標志碼： A? 文章編號： 1004-4523（2019）03-0415-06

DOI：10.16385/j.cnki.issn.1004-4523.2019.03.005

引 言

聲子晶體是由兩種或兩種以上介質組成的具有彈性波帶隙的周期性復合材料或結構。當彈性波在聲子晶體中傳播時，某些頻率范圍內的彈性波傳播將被抑制，相應的頻率范圍稱為帶隙。由于周期性結構廣泛存在于實際工程結構中，人們對周期性結構的研究有著悠久的歷史[1-3]。1992年，Sigalas和Economou研究了球形散射體埋入某一基體材料中形成的三維周期性復合介質中彈性波的傳播特性，首次從理論上證實了三維周期點陣結構中存在彈性波帶隙[4]。1993年，Kushwaha等在研究鎳/鋁二維固體周期復合介質時第一次提出了聲子晶體的概念，類比光子晶體分析了聲子完全帶隙在理論研究中的意義[5]。1995年，Martínez-Sala等對西班牙馬德里的一座具有兩百多年歷史的雕塑進行了聲學特性測試，該雕塑是由直徑為2.9 cm的中空不銹鋼圓柱周期性排布在一個4 m直徑的圓形平臺上，形成的晶格常數為10 cm，通過測試他們第一次從實驗角度證實了彈性波帶隙的存在[6]。2000年，劉正猷等將包覆軟硅橡膠材料的鉛球按立方晶格嵌入環氧樹脂基體中形成了一種三維三組元聲子晶體，理論和實驗結果均表明該聲子晶體帶隙所對應的波長遠大于晶格常數，由此提出了局域共振聲子晶體的概念[7]。近年來，局域共振聲子晶體超常物理特性的揭示引起了學者的關注，并將這種具有超常物理特性的復合材料統稱為超材料，其特性包括帶隙、負折射、負模量和聲學斗篷等[8-13]。

聲子晶體傳統能帶算法一般將Bloch定理作為周期邊界條件施加到求解域，先給定波矢k在不可約布里淵區邊界取值，然后求解特征頻率ω，從而得到能帶曲線（ω-k）[14-16]。由于傳統方法中波矢k只能取實數，所以這種求解方法存在一定的局限性，即只能得到實能帶。但是為了研究帶隙衰減能力，復能帶的計算常常更有價值，不僅可以得到實能帶，還可以得到帶隙內衰減常數。尤其是引入阻尼后，通帶內受到阻尼的耗散作用也會呈現一定的衰減，只有復能帶能夠很好地描述阻尼對帶隙的影響。而反過來，如果給定頻率ω，然后求解波矢k，就可以求解出復能帶（k-ω）。但是，特征波矢k將出現在邊界條件中，此時原問題變為非線性特征值問題，求解困難。這也是傳統方法一般都是給定波矢，然后求解特征頻率的原因。

為了能夠計算聲子晶體復能帶，本文提出了一種參數變換方法，可以將周期邊界條件中的未知波矢變換到偏微分方程的系數中，從而將原來復雜的非線性特征問題變換成線性特征值問題，能夠實現復能帶的快速求解。

3 算例與討論

本文將采用兩個算例來對文中算法進行驗證，其中算例一為布拉格聲子晶體板，算例二為局域共振聲子晶體板，算例中所用到的材料參數如表1所示。

兩個算例中，本文提出算法的理論計算都采用4節點Kirchhoff板單元，為了保證有限元計算結果收斂，分別對不同密度網格得到的結果進行了對比，結果表明采用12×12的網格已經能夠得到收斂的計算結果。為了驗證本文算法的正確性，還采用COMSOL軟件計算了兩個算例的實波矢，并與理論計算結果進行了對比。COMSOL軟件中元胞建模為三維模型，網格劃分采用二階拉格朗日單元，同樣對網格密度進行了分析，使得計算結果收斂。

3.1 布拉格聲子晶體板

為了形成布拉格聲子晶體板，基體材料選用環氧樹脂，散射體為結構鋼，材料參數如表1所示。晶格常數a=40 mm，散射體邊長l=20 mm，板厚為h=3 mm。采用本文算法得到復能帶如圖3所示，并且采用COMSOL軟件計算了該聲子晶體的實波矢。從圖中可以看出本文計算結果與COMSOL仿真結果符合的很好，但本文算法不僅可以計算實波矢，還可以計算波矢虛部從而得到帶隙內衰減常數。在0-3000 Hz范圍內，該聲子晶體在ΓX方向形成了一個布拉格帶隙（如圖3（a）所示），而在ΓM方向沒有產生帶隙（如圖3（b）所示），因此該帶隙具有方向性。在ΓX方向，帶隙內最大衰減常數約為δ=0.14，那么該頻率的彎曲波經過一個周期傳播后幅值衰減為e-0.14≈0.87。

為了研究帶隙隨方向的變化，計算了衰減常數隨方向角θ的變化， 如圖4所示。 從圖中可以看出帶隙寬度和帶隙位置都隨著方向而變化，帶隙寬度隨著角度θ的增大先減小后增大，帶隙頻率升高，而帶隙內衰減也是先減小后增大。特別是在2200 Hz附近，方向帶隙消失。

以1900 Hz為例，衰減常數隨方向的變化如圖5所示，從圖中可以看出衰減常數在0°，±90°和180°方向最大（約為0.13），而在±45°和±135°方向衰減常數為零。因此，該頻率的波在聲子晶體中傳播時具有明顯的方向特性。

由于本文算法可以計算出復能帶，所以可以用來分析阻尼對聲子晶體帶隙影響，這是傳統實數能帶無法描述的。在基體材料中引入結構阻尼系數η=0.01，計算得到復能帶如圖6所示。對比圖3和6可以看出，阻尼對布拉格帶隙的影響較小，尤其是實波矢基本沒有變化。但通帶內受到阻尼的影響，衰減不再為零，特別是帶邊附近出現了較強的阻尼衰減。

3.2 局域共振聲子晶體板

為了形成二組元局域共振聲子晶體板，基體材料選用結構鋼，散射體選用較為柔軟的橡膠，材料參數如表1所示。晶格常數為a=40 mm，散射體邊長為l=30 mm，板厚為h=3 mm。計算得到復能帶如圖7所示，同樣利用COMSOL軟件計算了該聲子晶體的實波矢。從圖中可以看出本文計算結果與COMSOL仿真結果也符合的很好，但本文方法同樣能夠給出局域共振帶隙內衰減常數。在0-200 Hz范圍內，該局域共振聲子晶體產生了一個局域共振帶隙。對比圖7（a）和（b）可以看出，該帶隙頻率范圍在ΓX方向和ΓM方向基本一致，證明該帶隙為一個完全帶隙，但從衰減常數可以看出，帶隙內衰減在這個兩個方向并不相同。可以看出，ΓX方向的帶隙內衰減要大于ΓM方向的帶隙內衰減，但衰減常數都是在接近下帶邊的頻率附近達到最大。

同樣，為了研究帶隙隨方向的變化，計算了衰減常數隨方向角θ的變化，如圖8所示。從圖中可以看出帶隙寬度和帶隙位置不隨方向變化，但帶隙內衰減隨著角度θ的增大而減小，在45°方向達到最小。

以150 Hz為例，衰減常數隨方向的變化如圖9所示，從圖中可以看出衰減常數在0°，±90°和180°方向最大（約為0.48），而在±45°和±135°最小（約為0.23）。那么150 Hz的波傳播一個周期后，在0°，±90°和180°方向幅值衰減為e-0.48≈0.62，而在±45°和±135°方向幅值衰減為e-0.23≈0.79。

為了分析阻尼對聲子晶體局域共振帶隙影響，在散射體中引入結構阻尼系數η=0.05，計算得到復能帶如圖10所示。對比圖7和10可以看出，增加阻尼后帶隙內波矢不再為純虛數，而帶隙外較大頻率范圍內波矢也不再為純實數，即帶隙內外較大范圍都將是復數，只有復能帶能很好地進行描述。阻尼使得帶隙內衰減減小，但帶邊附近卻出現了較強的阻尼衰減，證明引入適當的阻尼可以增大帶隙的衰減帶寬。與布拉格帶隙不同，阻尼對實波矢的影響很大，抑制了實波矢在帶隙附近的劇烈變化。

可以看出，阻尼對局域共振帶隙的影響遠大于布拉格帶隙，這是由兩種帶隙的形成機理決定的。布拉格聲子晶體板的基體和散射體材料參數相差較小，因而帶隙內振動能量在整個元胞內分布，基體阻尼只對部分機械能形成耗散;局域共振聲子晶體板的基體和散射體材料參數相差較大，因而帶隙內振動能量基本都被局域在散射體內，散射體阻尼將對幾乎全部機械能形成耗散，其帶隙內阻尼效果要遠比布拉格聲子晶體強。

4 結 論

本文提出了一種參數變換方法，解決了特征波矢求解中復雜的非線性特征問題，實現了復能帶的快速求解。文中通過兩個算例對算法進行了驗證，分別為布拉格聲子晶體板和局域共振聲子晶體板。布拉格聲子晶體板在0-3000 Hz內產生了一個布拉格帶隙，該帶隙具有方向性，在ΓX方向存在較大衰減，但在ΓM無衰減。基體阻尼對布拉格帶隙影響較小，但在帶邊處產生了顯著的阻尼衰減，擴大了帶隙的衰減范圍。局域共振聲子晶體在0-200 Hz內產生了一個局域共振帶隙，該帶隙為完全帶隙，但帶隙內衰減在不同方向稍有不同。散射體阻尼對局域共振帶隙影響顯著，帶隙內衰減減小，但帶邊處阻尼衰減較大。本文參數變換方法不限于聲子晶體板，也可推廣于其他二維和三維聲子晶體復能帶的計算。

參考文獻：

[1] Mead D J. Free wave propagation in periodic-supported，infinite beams[J].Journal of Sound and Vibration， 1970， 11（2）： 181-197.

[2] Mead D J， Markus S. Coupled flexural-longitudinal wave motion in a periodic beam[J]. Journal of Sound Vibration， 1983， 90（1）： 1-24.

[3] Mead D J. A new method of analyzing wave propagation in periodic structures： Applications to periodic Timoshenko beams and stiffened plates[J]. Journal of Sound Vibration， 1986， 104（1）： 9-27.

[4] Sigalas M M， Economou E N. Elastic and acoustic wave band structure[J]. Journal of Sound Vibration， 1992， 158（2）： 377-382.

[5] Kushwaha M S， Halevi P， Dobrzynski L， et al. Acoustic band structure of periodic elastic composites[J]. Physics Review Letter，1993，71（13）： 2022-2025.

[6] Martinez-Salar R， Sancho J， Sanchez J V， et al. Sound attenuation by sculpture[J]. Nature， 1995，378（6554）： 241.

[7] Liu Z Y， Zhang X， Mao Y， et al. Locally resonant sonic materials[J]. Science， 2000，289（5485）： 1734-1736.

[8] Li J， Chan C T. Double-negative acoustic metamaterial[J]. Physics Review E， 2004， 70（5）： 055602

[9] Fang N， Xi D， Xu J， et al. Ultrasonic metamaterials with negative modulus[J]. Nature Materials， 2006， 5（6）： 452-456.

[10] Milton G W. New metamaterials with macroscopic behavior outside that of continuum elastodynamics[J]. New Journal of Physics， 2007， 9（39）： 46-50.

[11] Li Y F， Meng F， Li S， et al. Designing broad phononic band gaps for in-plane modes[J]. Physics Letters A， 2018， 382（10）： 679-684.

[12] Srivastava A， Lu Y. Variational methods for phononic calculations[J]. Wave Motion， 2016， 60： 46-61.

[13] Coffy E， Euphrasie S， Addouche M， et al. Evidence of a broadband gap in a phononic crystal strip[J]. Ultrasonics， 2017， 78： 51-56.

[14] Chen S B， Wen J H， Wang G， et al. Tunable band gaps in acoustic metamaterials with periodic arrays of resonant shunted piezos[J]. Chinese Physics B， 2013， 22（7）： 074301.

[15] Yu D L， Wang G， Liu Y Z， et al. Flexural vibration band gaps in thin plates with two-dimensional binary locally resonant structures[J]. Chinese Physics， 2006， 15（2）： 266-271.

[16] Xiao Y， Wen J， Wen X. Flexural wave band gaps in locally resonant thin plates with periodically attached spring-mass resonators[J]. Journal of Physics D： Applied Physics， 2012， 45（19）： 195401.

Abstract： Phononic crystals possess elastic wave band-gaps， which can be used to control the vibration and noise of structures. To obtain the band structures of the phononic crystals， the conventional procedures are as follows： given the wave vector k， whose value sweep the boundary of Brillouin zone， the eigenfrequency ω can be evaluated， resulting in the ω-k curve. However， this method can only yield the real band structure. In order to get the complex band structure， the frequency ω usually is given and then the eigenvalue of wave vector k is calculated， resulting in the k-ω curve. This work proposes a parameter transformation method， which can resolve the complicated non-linear eigenvalue problem and achieve the rapid solution of complex band structure. Finally， two examples， i.e.， a Bragg phononic plate and a locally resonant phononic plate， are adopted to validate the proposed method. The variation of the attenuation constant along with the wave direction in the band gap and the influence of damping on the band gap are investigated in detail.

Key words： phononic crystals; metamaterials; complex band structure; parameter transformation

作者簡介： 陳圣兵（1984-），男，助理研究員。電話： （0816）2465610;E-mail：nudt_chen@163.com