一种在二维声子晶体的三维单胞能带中区分面内与面外模态的方法

1. 前言

        在形式多样的声子晶体研究中，能带结构是判断声子晶体基本性质的一个极其重要的参考指标。不过能带曲线通常只能反映出方向带隙、全向带隙的位置以及带隙上下边界所产生的位置，如果我们想了解更多，比如哪些能带是表示波动位移场集中在面内的模式，而哪些又是集中在面外的，普通的能带绘制方法就不能满足需求了。此处参考2020年发表在Extreme Mechanics Letters上的名为Flexural wave energy harvesting by multi-mode elastic metamaterial
cavities的文章中所提到的公式，提供一种可以用COMSO和MATLAB绘制区分面内和面外模态的能带曲线的方法，为能带曲线中模态极化方向的直观判断提供一种便利手段。

（注：标题中所说的二维声子晶体是指声子晶体单胞在二维平面内周期性排列，三维单胞能带是指考虑了x,y,z三个方向位移场的能带。）

2. COMSOL中的设置

FIG.1 单胞的几何构型

        单胞的几何构型如FIG.1所示，为一三叶草型孔隙结构，所用材料为铝。施加周期性边界条件之后进行特征频率的计算。计算结束之后，在结果中选取“派生值”(FIG.2)，右键插入体积分，在其中把数据集设置为“研究1/参数化解1”，计算域选择单胞整体。

FIG.2 体积分的计算

        根据前述文章所提供的的公式(FIG.3)，计算每个模态下位移场的面外极化率β。ux,uy,uz分别代表微元体在x,y,z方向的位移分量，ux*,uy*,uz*分别代表ux,uy,uz的共轭复数。此处之所以要使用共轭复数，是因为COMSOL在进行有限元计算时，会因数值计算产生虚部。为了抵消虚部的影响，故给每个分量乘以自己的虚部。如果数值计算不产生虚部的话，这里可以把乘以的共轭复数直接换为各个分量本身，形成一个平方项，再进行后续积分。根据FIG.3所提供的公式，在COMSOL中输入积分项，其中u2,v2,w2代表x,y,z方向的位移，conj()表示计算某个数共轭复数的算符。注意，此处的u2,v2,w2中的“2”表示这个能带计算在我的COMSOL模型中属于第二个组件，各位在计算的时候如果只有一个组件，这三个量应该是u1,v1,w1。

FIG.3 面外位移场极化率的计算公式

 

FIG.4 COMSOL中位移场积分的表达式

        设置后表达式之后，点击计算，很快就会在COMSOL界面的右下方生成一个表格(FIG.5)。点击红框中的导出键，将表格输出为.txt格式，以备MATLAB绘图之用。

FIG.5 计算之后的表格

3. MATLAB代码

        打开MATLAB，输入如下代码并点击运行：

n = 1 %要绘制的能带数量
A = cell(n,1) %整体能带数据的提取
for i = 1:n
    A{i,1} = readmatrix(['E:\My Blog\位移场积分信息',num2str(i),'.txt'])
    A{i,1}(:,5) = A{i,1}(:,4)./A{i,1}(:,3) %求出面外极化因子，并作为每一原胞数组的第五列
end

%区分面内及面外模态的能带曲线绘制
for i = 1:n
    x = A{i,1}(:,1); %数据第一列，表示波矢
    y = A{i,1}(:,2); %数据第二列，表示频率
    z = A{i,1}(:,5); %数据第五列，表示面外极化因子
    figure
    scatter(x,y,20,z,'filled')
    colorbar
    title(['Colored bands of unit cell ',num2str(i)])
    xlabel('K Vector','Fontname', 'Times New Roman','FontSize',18)
    ylabel('Frequency (Hz)','Fontname', 'Times New Roman','FontSize',18)
    set(0,'defaultfigurecolor','w')
    set(gca,'XTick',0:1:3)
    set(gca,'XTicklabel',{'M','\Gamma','X','M'})
    set(gca,'Fontname','Times New Roman','FontSize',15)
end

        运行之后，显示如FIG.6所示的彩色能带图，曲线中每一点的颜色表示该波动模态下面外位移场的极化率，该点颜色越黄，表示此模态下面外位移越占主要成分；相反，如果该点颜色越蓝，表示该模态下面内位移越占主要成分。FIG.7展示了普通绘制方法下相同能带信息的能带曲线图，对比之后可以明显看出本文中彩色能带图绘制方法的优势：即可以清晰地看出能带曲线中哪些点是面外或面内模态为主。希望本文所提供的方法能为各位研究声子晶体的朋友以及常用COMSOL和MATLAB的朋友带来一丝技术上的便利和启发：）

FIG.6 相应单胞的彩色能带图

 

FIG.7 相应单胞的普通能带图