二维光子晶体带隙结构计算与分析系统
项目介绍
本项目是一个基于 MATLAB 开发的专业仿真工具,旨在通过数值计算揭示二维光子晶体的能带结构特性。光子晶体是由具有不同介电常数的材料在空间呈周期性排列构成的人工微结构。本项目利用平面波展开法(Plane Wave Expansion Method, PWE),将麦克斯韦方程组转化为矩阵本征值问题,从而精确获取光子在特定晶格结构中的色散关系。
该系统核心价值在于能够计算并识别光子禁带(Photonic Bandgap),即某些特定频率范围的光波无法在该结构中传播。这种特性对于设计光子晶体光纤、低阈值激光器、光波导以及各类集成光子器件具有极其重要的科研和工程意义。
功能特性
- 多晶格支类型持:内置了正方形晶格和三角(六角)晶格的数学模型,适应不同的物理场景需求。
- 灵活的散射体建模:支持圆形和正方形两种主要的散射体形状,并可调节其在晶格中的填充比例。
- 高敛散性算法:采用了 Ho-Chan 方法(H-method),通过计算介电常数矩阵的逆矩阵提高计算能带时的收敛速度和精度。
- 全参数化设计:用户可以自由定义背景介质和散射体的介电常数对比度、平面波展开阶数、求解能带数量以及布里渊区采样密度。
- 自动化路径规划:系统能自动构建第一布里渊区内高对称点之间的扫频路径,包括 Gamma, X, M(正方形晶格)以及 Gamma, M, K(三角晶格)。
使用方法
- 设置物理参数:在程序开头配置晶格类型、散射体形状、介电常数(如硅为12.0,空气为1.0)以及填充率参数(半径或边长与晶格常数的比值)。
- 定义计算精度:通过调整平面波展开阶数 N_order 来平衡计算精度与运行时间。
- 设置能带数:指定需要求解的本征频率数量。
- 运行仿真:执行脚本后,系统将依次完成空间基矢计算、倒格矢生成、介电常数矩阵构建及能带路径采样。
- 结果分析:程序将处理本征值方程,并最终生成直观的能带结构图,展示扫频路径上的光子频率分布。
系统要求
- 运行环境:MATLAB R2016b 或更高版本。
- 硬件要求:由于涉及大规模矩阵运算(尤其是平面波阶数较高时),建议配备 8GB 以上内存。
功能实现逻辑
程序遵循严谨的电磁场数值计算流程,具体实现逻辑如下:
- 晶格空间定义:根据选定的晶格类型(正方形或三角),定义实空间基矢 a1, a2 和倒空间基本矢量 b1, b2。同时根据几何关系计算单胞面积。
- 倒格矢网络构建:在倒空间内基于展开阶数生成网格点。对于阶数为 N_order 的计算,将生成总量为 (2*N_order+1) 的平方个倒格矢点。
- 介电常数傅里叶变换:
- 程序根据散射体在单胞内的几何形状及其占据的面积比(填充率 f)进行傅里叶分量计算。
- 对于圆形散射体,利用一阶 Bessel 函数计算形状因子(Form Factor)。
- 对于正方形散射体,利用 sinc 函数形式计算空间结构因子。
- 构建介电常数矩阵 epsilon(G-G'),并执行逆矩阵运算生成 eta 矩阵,这是提高 PWE 稳定性的关键步骤。
- 定义高对称点坐标。对于正方形晶格为 [0,0] 到 [pi/a, 0] 再到 [pi/a, pi/a] 并回到原点。
- 在这些点之间进行线性插值,生成密集的采样点路径,确保能带图的平滑度。
关键算法与实现细节
这是本系统的算法核心。它假设电磁场可以在平面波基组下展开,将偏微分方程转化为代数特征值问题。
- 在计算非零倒格矢差(G - G')处的傅里叶系数时,代码区分了两种形状。
- 圆形逻辑:通过
besselj(1, x) / x 实现,这代表了圆盘结构的傅里叶变换。
- 正方形逻辑:采用
sin(x)/x(即 sinc 函数)的乘积形式,准确描述了矩形势垒在频域的分布。为了防止分母为零,代码中引入了极小值偏移 eps。
在处理 G - G' = 0 的特殊位置时,代码直接取用平均介电常数(由填充率加权计算),确保了零频分量的准确性。
代码通过预分配矩阵并使用双重循环遍历倒格矢集合来构建矩阵。利用 MATLAB 强大的线性代数库进行逆矩阵运算,为后续求解本征方程做准备。