基于COMSOL的聲子晶體帶結構計算新方法

付子義，王晨旭，王立國，長谷川弘治

?

基于COMSOL的聲子晶體帶結構計算新方法

付子義1*，王晨旭1，王立國1，長谷川弘治2

（1. 河南理工大學電氣工程與自動化學院，河南 焦作 454150；2. 室蘭工業大學情報電子工學系，日本 室蘭 0500071）

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

聲子晶體；能帶結構；有限元；PDE；本征方程

0 引言

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

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

本文主要的研究對象是聲子晶體的二維結構，不僅因為它們比三維結構更簡單和容易地構造，而且因為它們具有實際有用的性質。首先總結并比較了聲波和電磁波的基本特性，構建出光子晶體與聲子晶體之間物理量的聯系。其次，在眾多方法中，有限元方法憑借其良好的運算速度、精度、和收斂性，在各領域獲得了廣泛的應用，并已開發了許多商業化的有限元軟件，因此基于有限元仿真軟件COMSOL的系數偏微分模塊，由聲波方程出發，結合布洛赫態，推導出偏微分形式的標量方程，開發出快速，準確的聲子晶體能帶結構計算方法，簡化了數學建模過程，并且為進一步實現聲子晶體的優化提供一個簡單易行的辦法。

1 聲波與電磁波的聯系

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

本文總結并比較了聲波和電磁波的基本特性。首先，我們回顧一下描述聲波和電磁波行為的運動方程的歸一化版本。

聲波歸一化后的運動方程：

用于分析聲子晶體的聲波運動方程，可以由等式（2）導出以下形式：

電磁波歸一化后的運動方程如下：

用于分析光子晶體的電磁波運動方程，由等式（4）可以導出以下形式：

表1 二維聲波和電磁波之間的對應關系

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

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

2 聲子晶體能帶結構的數學推導

光子晶體能帶結構的求解過程是直接由空間部分的亥姆霍茲方程出發，分析電磁波的傳播特性，將矢量形式的布洛赫態轉化為標量形式的布洛赫態，推導出偏微分形式的本征方程，然后基于有限元仿真軟件COMSOL的系數偏微分模塊求解本征值。類比光子晶體能帶的求解，以下對聲子晶體建立數學模型。

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

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

上一節已經推導出用于分析聲子晶體的聲波運動方程，假設聲壓處于時諧模式：

則等式（3）可轉化為以下形式：

則等式（7）可轉化為：

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

將上式布洛赫態代入聲波運動方程（9）中，可得將等式轉化為如下標量形式：

3 二維聲子晶體能帶結構仿真

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

圖2 二維聲子晶體能帶結構，圖中實線和離散點分別為平面波法和限元法所得結果。

圖3 兩種方法對能帶結構中A點收斂性比較，圖中實線和離散點分別為有限元法和平面波法計算所得結果。

4 結論

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

[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] 齊共金, 楊盛良, 白書欣, 等. 基于平面波算法的二維聲子晶體帶結構的研究[J]. 物理學報, 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]. 固體力學學報(英文版), 2008, 21(6): 517-528.

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

[7] 劉啟能. 固-流結構聲子晶體中彈性波能帶的色散研究[J]. 人工晶體學報, 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

國家自然科學基金61501175

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

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

付子義，王晨旭，王立國，等. 基于COMSOL的聲子晶體帶結構計算新方法[J]. 軟件，2018，39（12）：06-09