基于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