航空工业陕飞公司 设计研究院 陕西汉中 723000
摘要：以正方晶格排布的二维固固布拉格散射型声子晶体为例，介绍了利用有限元方法研究声子晶体能带结构和传输特性的计算过程。基于传输特性计算，研究了声子晶体某一方向有限排布层数对其波动特性的影响，发现有限层声子晶体的传输系数随层数呈现指数线性的形式逐渐衰减。
关键词：声子晶体；有限元方法；能带结构；传输特性
1.引言
声子晶体是继电子晶体、光子晶体之后提出的一种具有弹性波超常特性的新型结构功能材料，其特性之一就是具有弹性波（声波）传播的频率禁带（即某些频率段的弹性波不能在结构中传播）。早期对声子晶体的研究集中在寻求声子晶体带隙计算方法、带隙形成机理、带隙影响因素和寻找宽带隙结构等几个方面。近几年，随着研究的深入，有关声子晶体的带隙调控、表面波、缺陷态、准周期或无序结构等应用基础研究逐渐成为了研究的重点。
目前，有关声子晶体波动特性的计算方法主要有传递矩阵法、平面波展开法、时域有限差分法（FDTD）、多重散射法、小波变换法等。以上几种方法均有各自的使用限制或缺点，如平面波展开法难以应用在大弹性常数差和固流声子晶体带隙的计算上，FDTD方法数值稳定性和收敛性要求高、计算费时等。随着声子晶体应用基础研究的深入，声子晶体的计算模型越来越贴近于实际情况，需要考虑结构的尺寸、边界及外界环境对其波动特性的影响。与之相应的一种普适的、简单的方法，必然受到相关研究人员的青睐。而有限元方法正具有这一特色，它可以计算任意复杂的一维、二维、三维固固和固流声子晶体的波动特性，未来必然在偏于实际应用的声子晶体带隙调控、表面波、缺陷态等方面有广泛的应用。
对于固固声子晶体，有限元方法可以很好的计算无限大声子晶体的能带结构和有限结构声子晶体的传输特性。在以往的研究中，传输特性一般被用来辅助说明结构带隙的存在、或计算能带方法的正确性，却很少被用来严格确定禁带的起始和截止频率。本文给出了利用有限元方法计算固固声子晶体能带和传输特性的理论和过程，计算了二维二组元固固布拉格散射型声子晶体的能带结构和传输特性，说明了可以从传输特性上来确定固固声子晶体禁带起始和截止频率；从传输特性计算方面，研究了某一方向周期排布层数对有限结构声子晶体波动特性的影响，发现有限层声子晶体的传输系数随层数呈指数线性衰减。
2.能带结构计算
对于声子晶体，周期势场下的弹性波波动方程为
.png)
（1）
其中，r=(x,y,z)为质点位置矢量，ρ(r)材料密度，C(r)为材料弹性张量，u(r,t)为弹性波波场位移，R为晶格任意整数倍的平移周期。由Bloch定理可知，式（1）在频域内进行有限元离散求解可得到简约波矢量下的动刚度矩阵方程[1]为：
（K（k）-ω2M（k））×Uk=0 （2）
式（2）为随波矢K取值变化的复数矩阵的广义特征值方程。K沿简约布里渊区的边界取值时，可求得一系列的实特征值ω。所有K取值下同阶特征频率的连线即为色散曲线，曲线之间的空隙即为带隙。
3.传输特性计算
由于一般平面波展开法难于计算固流声子晶体的禁带特性，人们后来将电磁波计算中的FDTD方法拓展来计算声子晶体禁带特性。FDTD方法可以计算声子晶体的能带结构和声子晶体的透射和反射系数，但其计算量大，耗费机时，编程复杂。本文传输特性计算为限元商业软件上的声子晶体频率响应分析,免于手编程序，方便选择频率段进行计算。
一般的无阻尼多自由系统在正弦激励下的频域运动方程为（K-ω2M）x=f ，式中，K和M 为分别为系统刚度矩阵和质量矩阵，f和x分别是输入和响应向量。则多自由度系统的第k, j自由度的频响函数可以写为：
.png)
式中， adj表示伴随矩阵，结构的频率响应是输出点对输入点的响应，不同的输出点可能对同一点输入的响应有所差别。
对于固固声子晶体的频响计算，无限大模型取为有限层进行计算，当考察
.png)
方向的传输特性时，在模型左边界施加沿
.png)
方向的特定幅值的单频正弦加速度激励，上下边界施加位移周期边界条件，以表示模型在
.png)
方向无限大，然后以足够小的频率间隔扫频，以便得到响应幅值的峰值，最后拾取右边界某点的加速度响应幅值，响应幅值与激励幅值之比即为输出点对输入点的频率响应。
4.算例模型
.png)
图1 二维固固布拉格散射型声子晶体
.png)
图 2 二维铝/环氧树脂声子晶体XY模式能带结构
.png)
图3 6层铝/环氧树脂声子晶体x方向传输特性
本文取图1所示的具有相同几何参数的正方晶格排列的二维声子晶体为例用有限元方法计算了相应的固固声子晶体的波动特性。晶格常数a=0.02mm，散射体半径r=0.008mm。当计算能带结构时，其计算单胞如图 1（b）所示，图 1（d）中的灰色区域为晶格倒格矢空间的不可约布里渊区。当计算沿Γx方向的传输特性时，在x方向取有限层，在y方向取单层，如图1（c），在y方向边界施加周期边界条件，以表示y方向为无限大。
对于固固声子晶体，图1（a）中白色部分为基体环氧树脂，灰色部分为散射体铝。具体参数为[2]：ρepo=1180kg/m3，λepo=4.43×109，μepo=1.59×109；铝材料参数为ρalu=2730kg/m3，λalu=6.89×1010，μalu=2.87×1010。ρ和λ，μ分别表示材料的密度和拉梅常数。
5.结果分析
图2为160kHz频率以下的二维铝/环氧树脂声子晶体XY模式能带结构。此图与文献[2]中用平面波展开法得出的XY模式能带结构图相同。图中，灰色区域的部分为XY模式下的第一完全禁带，其带隙频率范围为56.7kHz85.5kHz，三个格子填充的部分为Bloch波沿Γx方向传播的方向禁带，第一Γx方向带隙频率范围为55.7kHz-89.4kHz，略大于第一完全带隙频率。
通过和图 2灰色区域的局部放大图比较可以发现，图 3中响应在100以下的三个区域频率范围和图 2中方向禁带的三个区域频率范围吻合得很好。因此，我们可以从频率响应计算的传输特性图上来确定声子晶体的带隙频率。从第一禁带中心频率附近某频率的声压幅值分布图可以看出，声波的能量主要集中在第一层前半部分，其衰减程度呈指数线性形式。
6.结论
有限元方法既可以计算无限大声子晶体的能带结构，又可以计算有限大声子晶体的传输特性。传输特性图上的频率衰减范围与能带结构图中的带隙频率范围完全相同。传输特性图给出比能带结构图更多的声子晶体的波动信息。声波在有限结构声子晶体中传播衰减呈指数线性形式。
参考文献
[1] Aberg, M. and P. Gudmundson, The usage of standard finite element codes for computation of dispersion relations in materials with periodic microstructure. Journal of the Acoustical Society of America, 1997. 102(4):2007-2013.
[2] 郁殿龙, 基于声子晶体理论的梁板类周期结构振动带隙特性研究. 2006.