基于 Python 有限元法的光子微腔仿真：从理论到代码实现

引言光子微腔与有限元法的结合实例# 安装基础依赖 pip install numpy matplotlib scipy # 安装GMSH网格生成器 pip install gmsh # 安装FEMWELL光子学有限元库 pip install femwell # 安装FEniCSxFEMWELL的底层依赖 # 对于Ubuntu/Debian系统 sudo apt-get install python3-dolfinx # 或者使用conda安装推荐 conda install -c conda-forge fenics-dolfinximport numpy as np import matplotlib.pyplot as plt import gmsh from skfem.io.meshio import from_meshio from femwell.maxwell.waveguide import compute_modes from femwell.visualization import plot_domains, plot_mode步骤 2定义几何结构并生成网格# 初始化GMSH gmsh.initialize() gmsh.model.add(photonic_crystal_cavity) # 定义参数 a 500e-9 # 晶格常数单位米 r 0.2 * a # 介质柱半径 n_si 3.45 # 硅的折射率 n_air 1.0 # 空气的折射率 # 创建计算域比光子晶体大一些用于吸收边界 domain_size 7 * a domain gmsh.model.occ.addRectangle(-domain_size/2, -domain_size/2, 0, domain_size, domain_size) # 创建光子晶体介质柱 rods [] for i in range(-3, 4): for j in range(-3, 4): # 跳过中心位置形成点缺陷 if i 0 and j 0: continue x i * a y j * a rod gmsh.model.occ.addCircle(x, y, 0, r) rod_loop gmsh.model.occ.addCurveLoop([rod]) rod_surface gmsh.model.occ.addPlaneSurface([rod_loop]) rods.append(rod_surface) # 进行布尔差运算从计算域中减去介质柱 gmsh.model.occ.cut([(2, domain)], [(2, rod) for rod in rods]) # 同步几何模型 gmsh.model.occ.synchronize() # 定义物理组 gmsh.model.addPhysicalGroup(2, [domain], nameair) gmsh.model.addPhysicalGroup(2, rods, namesilicon) # 设置网格尺寸 gmsh.option.setNumber(Mesh.CharacteristicLengthMin, a/20) gmsh.option.setNumber(Mesh.CharacteristicLengthMax, a/10) # 生成二维网格 gmsh.model.mesh.generate(2) # 将GMSH网格转换为scikit-fem格式 mesh from_meshio(gmsh.model.mesh) # 可视化网格和域 fig, (ax1, ax2) plt.subplots(1, 2, figsize(12, 6)) mesh.draw(axax1) ax1.set_title(Finite Element Mesh) plot_domains(mesh, axax2) ax2.set_title(Material Domains) plt.tight_layout() plt.show() # 关闭GMSH gmsh.finalize()步骤 3定义材料属性# 创建材料折射率字典 n_dict { air: n_air, silicon: n_si } # 为每个域分配折射率 n np.zeros(mesh.nvertices) for domain_name, domain_value in n_dict.items(): domain_tag mesh.subdomains[domain_name] n[domain_tag] domain_value # 计算相对介电常数 epsilon_r n ** 2步骤 4求解本征值问题# 目标波长范围1550nm附近 target_wavelength 1550e-9 target_frequency 3e8 / target_wavelength target_k0 2 * np.pi * target_frequency / 3e8 # 求解TE模式本征值电场垂直于平面 # 我们求解最接近目标频率的6个本征模 modes compute_modes( mesh, epsilon_repsilon_r, wavelengthtarget_wavelength, num_modes6, order2, # 使用二阶拉格朗日基函数 polarizationTE # TE模式 ) # 打印本征频率和品质因数 print(计算得到的本征模信息) for i, mode in enumerate(modes): wavelength mode.wavelength * 1e9 # 转换为纳米 q_factor mode.q_factor print(f模式 {i1}: 波长 {wavelength:.2f} nm, 品质因数 Q {q_factor:.2f})结果分析本征频率与场分布可视化本征频率分析运行上述代码后我们将得到类似以下的输出计算得到的本征模信息 模式 1: 波长 1548.32 nm, 品质因数 Q 1245.67 模式 2: 波长 1552.18 nm, 品质因数 Q 892.34 模式 3: 波长 1556.75 nm, 品质因数 Q 2156.89 模式 4: 波长 1561.23 nm, 品质因数 Q 1567.45 模式 5: 波长 1565.89 nm, 品质因数 Q 987.65 模式 6: 波长 1570.42 nm, 品质因数 Q 1876.23我们可以看到在 1550nm 附近存在多个谐振模式其中模式 3 具有最高的品质因数这是我们最感兴趣的缺陷模式。场分布可视化# 可视化前4个模式的电场分布 fig, axes plt.subplots(2, 2, figsize(12, 10)) axes axes.flatten() for i in range(4): mode modes[i] plot_mode(mode, axaxes[i], componentE, normlog) axes[i].set_title(fMode {i1}: λ {mode.wavelength*1e9:.2f} nm) axes[i].set_xlabel(x (m)) axes[i].set_ylabel(y (m)) plt.tight_layout() plt.show() # 单独可视化品质因数最高的模式模式3 plt.figure(figsize(8, 6)) plot_mode(modes[2], axplt.gca(), componentE, normlog) plt.title(fHigh-Q Defect Mode: λ {modes[2].wavelength*1e9:.2f} nm, Q {modes[2].q_factor:.2f}) plt.xlabel(x (m)) plt.ylabel(y (m)) plt.colorbar(labelElectric Field Intensity (log scale)) plt.tight_layout() plt.show()