聲子晶體板復(fù)能帶計(jì)算方法

陳圣兵　張浩　宋玉寶

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

關(guān)鍵詞： 聲子晶體; 超材料; 復(fù)能帶; 參數(shù)變換

中圖分類號(hào)： O735; TB535? 文獻(xiàn)標(biāo)志碼： A? 文章編號(hào)： 1004-4523（2019）03-0415-06

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

引 言

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

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

為了能夠計(jì)算聲子晶體復(fù)能帶，本文提出了一種參數(shù)變換方法，可以將周期邊界條件中的未知波矢變換到偏微分方程的系數(shù)中，從而將原來復(fù)雜的非線性特征問題變換成線性特征值問題，能夠?qū)崿F(xiàn)復(fù)能帶的快速求解。

3 算例與討論

本文將采用兩個(gè)算例來對(duì)文中算法進(jìn)行驗(yàn)證，其中算例一為布拉格聲子晶體板，算例二為局域共振聲子晶體板，算例中所用到的材料參數(shù)如表1所示。

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

3.1 布拉格聲子晶體板

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

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

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

由于本文算法可以計(jì)算出復(fù)能帶，所以可以用來分析阻尼對(duì)聲子晶體帶隙影響，這是傳統(tǒng)實(shí)數(shù)能帶無法描述的。在基體材料中引入結(jié)構(gòu)阻尼系數(shù)η=0.01，計(jì)算得到復(fù)能帶如圖6所示。對(duì)比圖3和6可以看出，阻尼對(duì)布拉格帶隙的影響較小，尤其是實(shí)波矢基本沒有變化。但通帶內(nèi)受到阻尼的影響，衰減不再為零，特別是帶邊附近出現(xiàn)了較強(qiáng)的阻尼衰減。

3.2 局域共振聲子晶體板

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

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

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

為了分析阻尼對(duì)聲子晶體局域共振帶隙影響，在散射體中引入結(jié)構(gòu)阻尼系數(shù)η=0.05，計(jì)算得到復(fù)能帶如圖10所示。對(duì)比圖7和10可以看出，增加阻尼后帶隙內(nèi)波矢不再為純虛數(shù)，而帶隙外較大頻率范圍內(nèi)波矢也不再為純實(shí)數(shù)，即帶隙內(nèi)外較大范圍都將是復(fù)數(shù)，只有復(fù)能帶能很好地進(jìn)行描述。阻尼使得帶隙內(nèi)衰減減小，但帶邊附近卻出現(xiàn)了較強(qiáng)的阻尼衰減，證明引入適當(dāng)?shù)淖枘峥梢栽龃髱兜乃p帶寬。與布拉格帶隙不同，阻尼對(duì)實(shí)波矢的影響很大，抑制了實(shí)波矢在帶隙附近的劇烈變化。

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

4 結(jié) 論

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

參考文獻(xiàn)：

[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