付子義，王晨旭，王立國(guó)，長(zhǎng)谷川弘治

?

基于COMSOL的聲子晶體帶結(jié)構(gòu)計(jì)算新方法

付子義1*，王晨旭1，王立國(guó)1，長(zhǎng)谷川弘治2

（1. 河南理工大學(xué)電氣工程與自動(dòng)化學(xué)院，河南 焦作 454150；2. 室蘭工業(yè)大學(xué)情報(bào)電子工學(xué)系，日本 室蘭 0500071）

結(jié)合聲波與電磁波的物理特性，建立光子晶體與聲子晶體物理量之間的聯(lián)系，通過(guò)類比光子晶體能帶結(jié)構(gòu)的求解過(guò)程，對(duì)聲子晶體能帶結(jié)構(gòu)進(jìn)行仿真計(jì)算。以二維聲子晶體為研究對(duì)象，通過(guò)重新建立數(shù)學(xué)模型，利用Born-von-karman周期性邊界條件和布洛赫態(tài)，推導(dǎo)出聲子晶體偏微分形式的本征方程，基于有限元仿真軟件COMSOL Multiphysics系數(shù)偏微分模塊（coefficient form PDE），求解出相應(yīng)的本征頻率從而求解出能帶。與已有計(jì)算方法相比，該方法在適用性，難易程度等方面具有明顯的優(yōu)越性。

聲子晶體；能帶結(jié)構(gòu)；有限元；PDE；本征方程

0 引言

聲子晶體是由嵌入均勻主體材料中的聲波散射體組成的有限尺寸周期陣列的人造晶體。聲波與周期性排列的散射體之間的相互作用導(dǎo)致形成聲波不能通過(guò)該結(jié)構(gòu)的頻帶。這種頻帶被稱為帶隙?；谶@一特性使得聲子晶體具有廣闊的應(yīng)用前景。

聲子晶體的概念是由光子晶體的概念演繹而來(lái)的，因此與光子晶體相類似的是，能帶結(jié)構(gòu)也是聲子晶體最重要的研究?jī)?nèi)容之一。文獻(xiàn)[1]基于TMM方法，對(duì)聲子晶體能帶結(jié)構(gòu)進(jìn)行求解計(jì)算，但是TMM方法主要用于一維聲子晶體的能帶計(jì)算。雖然結(jié)合PWE方法可用于二維情況，但計(jì)算量較大。文獻(xiàn)[2]基于PWE方法，但是PWE方法對(duì)于組分材料聲學(xué)參數(shù)相差較大的聲子晶體需要更多的計(jì)算時(shí)間和內(nèi)存，而且不易給出正確結(jié)果。文獻(xiàn)[3]采用MST方法，此方法理論推導(dǎo)復(fù)雜，并且只限于處理圓柱及球形散射體單元結(jié)構(gòu)的聲子晶體，應(yīng)用上存在較大的局限性。文獻(xiàn)[4]采用的FDTD方法可以模擬各種復(fù)雜的聲子晶體結(jié)構(gòu)但當(dāng)組分材料的聲學(xué)參數(shù)差異較大時(shí)，需要更細(xì)的網(wǎng)格來(lái)達(dá)到更高的精度，計(jì)算量較大。

本文主要的研究對(duì)象是聲子晶體的二維結(jié)構(gòu)，不僅因?yàn)樗鼈儽热S結(jié)構(gòu)更簡(jiǎn)單和容易地構(gòu)造，而且因?yàn)樗鼈兙哂袑?shí)際有用的性質(zhì)。首先總結(jié)并比較了聲波和電磁波的基本特性，構(gòu)建出光子晶體與聲子晶體之間物理量的聯(lián)系。其次，在眾多方法中，有限元方法憑借其良好的運(yùn)算速度、精度、和收斂性，在各領(lǐng)域獲得了廣泛的應(yīng)用，并已開發(fā)了許多商業(yè)化的有限元軟件，因此基于有限元仿真軟件COMSOL的系數(shù)偏微分模塊，由聲波方程出發(fā)，結(jié)合布洛赫態(tài)，推導(dǎo)出偏微分形式的標(biāo)量方程，開發(fā)出快速，準(zhǔn)確的聲子晶體能帶結(jié)構(gòu)計(jì)算方法，簡(jiǎn)化了數(shù)學(xué)建模過(guò)程，并且為進(jìn)一步實(shí)現(xiàn)聲子晶體的優(yōu)化提供一個(gè)簡(jiǎn)單易行的辦法。

1 聲波與電磁波的聯(lián)系

聲子晶體實(shí)際上是光子晶體的聲波版本[5-7]，其主體材料可以是固體，在固體基質(zhì)材料中，同時(shí)存在縱波和橫向剪切波并且彼此耦合。相反，當(dāng)主體材料為流體時(shí)，聲波在流體中傳播是只存在縱波，橫波不必考慮[8]。作為本文理論模型的驗(yàn)證，我們采用文獻(xiàn)[2]所用參數(shù)，以便對(duì)比，因此本文采用水銀/和四氯化碳組成的二維雙組分液相體系，在此只需考慮縱波即可。

本文總結(jié)并比較了聲波和電磁波的基本特性。首先，我們回顧一下描述聲波和電磁波行為的運(yùn)動(dòng)方程的歸一化版本。

聲波歸一化后的運(yùn)動(dòng)方程：

用于分析聲子晶體的聲波運(yùn)動(dòng)方程，可以由等式（2）導(dǎo)出以下形式：

電磁波歸一化后的運(yùn)動(dòng)方程如下：

用于分析光子晶體的電磁波運(yùn)動(dòng)方程，由等式（4）可以導(dǎo)出以下形式：

表1 二維聲波和電磁波之間的對(duì)應(yīng)關(guān)系

Tab.1 A complete correspondence between two-dimensional sound waves and electromagnetic waves

聲波的運(yùn)動(dòng)方程，等式（2）僅包括梯度或散度的操作，而電磁波的運(yùn)動(dòng)方程，等式（3）通過(guò)旋轉(zhuǎn)的矢量操作來(lái)描述。這些方程本質(zhì)上看起來(lái)完全不同，沒有對(duì)應(yīng)關(guān)系。然而，在二維情況下，在笛卡爾坐標(biāo)中發(fā)現(xiàn)它們之間完全對(duì)應(yīng)，如表1所示。TE波的磁場(chǎng)H以及TM波的電場(chǎng)E表現(xiàn)得像標(biāo)量的聲壓P。TE波的橫向分量E以及TM波的H垂直于傳播方向，而聲波的粒子速度與傳播方向平行。

2 聲子晶體能帶結(jié)構(gòu)的數(shù)學(xué)推導(dǎo)

光子晶體能帶結(jié)構(gòu)的求解過(guò)程是直接由空間部分的亥姆霍茲方程出發(fā)，分析電磁波的傳播特性，將矢量形式的布洛赫態(tài)轉(zhuǎn)化為標(biāo)量形式的布洛赫態(tài)，推導(dǎo)出偏微分形式的本征方程，然后基于有限元仿真軟件COMSOL的系數(shù)偏微分模塊求解本征值。類比光子晶體能帶的求解，以下對(duì)聲子晶體建立數(shù)學(xué)模型。

本文研究對(duì)象為水銀/四氯化碳體系的簡(jiǎn)單立方晶格結(jié)構(gòu)的聲子晶體，每一個(gè)聲子晶體原胞只含有一個(gè)散射體。所采用的文獻(xiàn)[2]的材料特性參數(shù)_ = 0.1,ˉ = 20。假設(shè)聲子晶體是由圓柱體A以正方點(diǎn)陣形式排列與基體B中組成的二維雙組分復(fù)合材料體系。晶格常數(shù)為a，圓柱體橫截面半徑為r=0.3a。圓柱體軸線與z軸平行，二維周期性結(jié)構(gòu)分布在x-y坐標(biāo)平面內(nèi)，如圖1所示。

圖1 （a）二維聲子晶體的橫截面示意圖，（b）波矢沿第一布里淵區(qū)域移動(dòng)

上一節(jié)已經(jīng)推導(dǎo)出用于分析聲子晶體的聲波運(yùn)動(dòng)方程，假設(shè)聲壓處于時(shí)諧模式：

則等式（3）可轉(zhuǎn)化為以下形式：

則等式（7）可轉(zhuǎn)化為：

周期結(jié)構(gòu)中的波傳播可以通過(guò)使用布洛赫定理限制為單個(gè)單元，在這種周期性結(jié)構(gòu)中傳播的波由以下形式所描述：

將上式布洛赫態(tài)代入聲波運(yùn)動(dòng)方程（9）中，可得將等式轉(zhuǎn)化為如下標(biāo)量形式：

3 二維聲子晶體能帶結(jié)構(gòu)仿真

利用上一部分建立的數(shù)學(xué)模型，結(jié)合COMSOL系數(shù)偏微分（coefficient form PDE）模塊，可以開展對(duì)二維聲子晶體能帶結(jié)構(gòu)的仿真計(jì)算。使用者可以根據(jù)實(shí)際問(wèn)題在數(shù)學(xué)方程上進(jìn)行修改，具有極大的靈活性，因此某些特殊的不便于用現(xiàn)成專題模塊描述的問(wèn)題可通過(guò)數(shù)學(xué)模塊獲得合理的解決。

圖2 二維聲子晶體能帶結(jié)構(gòu)，圖中實(shí)線和離散點(diǎn)分別為平面波法和限元法所得結(jié)果。

圖3 兩種方法對(duì)能帶結(jié)構(gòu)中A點(diǎn)收斂性比較，圖中實(shí)線和離散點(diǎn)分別為有限元法和平面波法計(jì)算所得結(jié)果。

4 結(jié)論

通過(guò)總結(jié)并比較聲波與電磁波之間的物理特性和聯(lián)系，從光子晶體能帶結(jié)構(gòu)求解出發(fā)，建立求解聲子晶體能帶結(jié)構(gòu)的數(shù)學(xué)模型，結(jié)合有限元軟件COMSOL Multiphysics系數(shù)偏微分模塊（coefficient form PDE）對(duì)此數(shù)學(xué)模型進(jìn)行仿真計(jì)算，數(shù)值結(jié)果與傳統(tǒng)的平面波展開法吻合較好，能大大提高收斂性和計(jì)算速度。不同于以往需要大量復(fù)雜的數(shù)學(xué)推導(dǎo)，本文從聲波的運(yùn)動(dòng)方程出發(fā)，結(jié)合布洛赫態(tài)，推導(dǎo)出偏微分形式的本征方程，通過(guò)COMSOL系數(shù)偏微分模塊內(nèi)置有限元算法直接對(duì)歸一化后的本征方程進(jìn)行數(shù)值計(jì)算，具有極大的靈活性。此方法對(duì)于周期性人造晶體結(jié)構(gòu)能帶求解具有較好的借鑒意義，為今后研究更復(fù)雜聲子晶體模型提供了一種更為便捷的數(shù)值仿真手段。

[1] Hou Z, Kuang W, Liu Y. Transmission property analysis of a two-dimensional phononic crystal[J]. Physics Letters A, 2004, 333(1): 172-180.

[2] 齊共金, 楊盛良, 白書欣, 等. 基于平面波算法的二維聲子晶體帶結(jié)構(gòu)的研究[J]. 物理學(xué)報(bào), 2003, 52(3): 000668-671.

[3] Mei J, Liu Z, Shi J, et al. Theory for elastic wave scattering by a two-dimensional periodical array of cylinders: An ideal approach for band-structure calculations[J]. Physical Review B, 2003, 67(24): 841-845.

[4] Tanaka Y, Tamura S I. Band structures of acoustic waves in phononic lattices[J]. Physica B Physics of Condensed Matter, 2002, 316(1): 237-239.

[5] Chen A, Wang Y, Yu G, et al. Elastic wave localization in two-dimensional phononic crystals with one-dimensional quasi-periodicity and random disorder[J]. 固體力學(xué)學(xué)報(bào)(英文版), 2008, 21(6): 517-528.

[6] 溫激鴻, 王剛, 劉耀宗, 等. 基于集中質(zhì)量法的一維聲子晶體彈性波帶隙計(jì)算[J]. 物理學(xué)報(bào), 2004, 53(10): 3384-3388.

[7] 劉啟能. 固-流結(jié)構(gòu)聲子晶體中彈性波能帶的色散研究[J]. 人工晶體學(xué)報(bào), 2009, 38(1): 112-117.

[8] Miyashita T. Sonic crystals and sonic wave-guides[J]. Measurement Science & Technology, 2005, 16(5): 47-63.

A New Method for Computation of Phononic Crystals Band Structure by COMSOL

FU Zi-yi1*, WANG Chen-xu1, WANG Li-guo1, HASEGAWA Koji2

(1. Henan Polytechnic University, Henan Jiaozuo 454150; 2. Muroran Institute of Technology, Japen Shilan 0500071)

Combining the physical properties of acoustic waves and electromagnetic waves, the relationship between the photonic crystal and the physical quantity of phononic crystal is established. The phononic crystal band structure is calculated by analogy photonic crystal band structure. Taking the two-dimensional phononic crystal as the research object, the eigenequation of the partial differential form of the phononic crystal is derived by reconstructing the mathematical model and using the Born-von-karman periodic boundary condition and the Bloch state. Based on the finite element simulation software COMSOL The Multiphysics coefficient form (PDE) solves the corresponding eigenfrequency to solve the energy band. Compared with the existing calculation methods, this method has obvious advantages in terms of applicability and difficulty.

Phononic crystals; Band structure; FEM; PDE; Eigenfunction

0735

A

10.3969/j.issn.1003-6970.2018.12.002

國(guó)家自然科學(xué)基金61501175

王晨旭(1994-)，女，碩士研究生，研究方向：聲子晶體與光子晶體。

付子義(1958-)，男，教授，博士生導(dǎo)師，研究方向：智能信號(hào)處理。

付子義，王晨旭，王立國(guó)，等. 基于COMSOL的聲子晶體帶結(jié)構(gòu)計(jì)算新方法[J]. 軟件，2018，39（12）：06-09