基于COMSOL的声子晶体能带结构仿真方法
声子晶体作为一种周期性人工复合结构，能够调控弹性波传播。COMSOL软件为计算其能带结构（即色散曲面）提供了完整解决方案。本文详细介绍从几何建模到结果分析的全流程。
模型建立基础流程
1. 几何结构定义
以二维正方晶格为例，包含散射体（如圆柱形）和基体材料。晶格常数设为a，散射体半径为r。在COMSOL中通过几何模块直接绘制：先创建矩形单元区域（尺寸a×a），再在中心添加圆形域（半径r）。利用阵列功能复制单元，形成周期结构。
下面是用MATLAB生成类似几何参数的示意代码（仅展示逻辑，非COMSOL直接代码）：
a = 1e-3; % 晶格常数，单位m r = 0.35e-3; % 散射体半径 % 计算填充比 f = pi*r^2 / a^2; % 生成单个单元的几何坐标 x = linspace(-a/2, a/2, 100); y = linspace(-a/2, a/2, 100); [X, Y] = meshgrid(x, y); mat = (X.^2 + Y.^2) <= r^2; % 散射体区域
在COMSOL实际建模时，可直接使用"几何"工作台中的"圆"和"矩形"工具，并应用"布尔运算"实现差值或并集。
2. 材料参数输入
假设基体为硅（Si），散射体为空气。需在"材料"节点中定义密度、弹性模量和泊松比。等效参数如下：
- 硅：密度ρ=2330 kg/m³，杨氏模量E=1.69×10¹¹ Pa，泊松比ν=0.28
- 空气：密度ρ=1.2 kg/m³，声速c=343 m/s（使用流体模型时）
对于固-气系统，通常将空气等效为极低阻抗材料，或直接使用"压力声学"物理场。
3. 物理场与方程配置
选择"固体力学"接口，弹性波传播满足Navier方程：
ρ ∂²u/∂t² = ∇·σ + f
其中u为位移场，σ为应力张量，f为体力。COMSOL自动离散方程，用户只需添加"周期性条件"于晶胞边界。
4. 边界条件和求解设置
在晶胞的四个边界上应用Floquet周期性边界条件：
// 示意命令（非COMSOL脚本）： // 左侧边界 u_left = u_right * exp(i*kx*a) // 底部边界 u_bottom = u_top * exp(i*ky*a)
在"研究"节点中，选择"频域"或"特征频率"求解。若求解能带，需设置波矢k沿不可约布里渊区路径扫描（如Γ-X-M-Γ）。每个k点求解特征频率，即可获得色散曲面。
求解器建议采用"MUMPS"直接求解器，对于三维或大模型可换用迭代求解器（如GMRES）。
结果解读与应用
计算后，通过"一维绘图组"可绘制频率vs波矢的色散曲线。图中清晰的频率禁带（带隙）表明弹性波在该频段无法传播。带隙位置和宽度直接影响声子晶体的滤波或隔声性能。
例如，若在3000-5000 Hz频率范围内无色散曲线，则说明存在约2 kHz的带隙。通过调整几何参数（如填充比、散射体形状）可调控带隙。
