COMSOL 光栅衍射仿真工具

English | 中文

基于 COMSOL Multiphysics Wave Optics 模块的 2D 矩形光栅衍射仿真 Python 工具。采用全场 + 周期端口 + 衍射阶 + Floquet + TE方案，可精确计算总反射率/透射率及各级次衍射效率。

🏆 验证状态：经过 66 个案例全面验证（N_teeth=5..26, θ ∈ {±66.36°, 0°}），能量守恒 R+T = 1.00000，精度优于 1×10⁻⁵。

✨ 功能特性

功能	说明
总反射/透射率	R 和 T，精度达 5-6 位有效数字
各级次衍射效率	提取每个传播 Bragg/超胞级次的 R_m / T_m
多层光栅支持	基底、匹配层、DBR 多层结构
Littrow 光栅设计验证	判断是否命中目标角度
场分布可视化	生成可在 COMSOL GUI 中查看的已求解模型

🚫 不适用场景

3D 光栅（本工具仅支持 2D）
TM 偏振（存在已知未解决的端口 S 参数问题，详见 BUGS_AND_FIXES.md Bug-9）
时域或本征频率分析
有损材料（本工具假设无损 k=0 材料；可传入复数折射率 n，但需自行验证 R+T+A=1）

📁 项目结构

comsol-skill/
├── README.md                       ← 英文说明
├── readme_zh.md                    ← 中文说明（本文件）
├── SKILL.md                        ← Agent 使用规范
├── BUGS_AND_FIXES.md               ← 9 个已知 API 陷阱及修复方案
├── COMSOL_MATLAB_Pitfalls_Guide.md ← 遗留笔记
└── grating-skill/
    ├── comsol_grating_template.py  ← 主脚本模板
    ├── example_config.json         ← 最小示例配置
    ├── example_run.bat             ← Windows 双击运行
    └── requirements.txt            ← Python 依赖

🔧 环境要求

软件依赖

COMSOL Multiphysics 6.x（需要 Wave Optics 模块许可证）
Python 3.8+

Python 包

mph>=1.2
jpype1>=1.4
numpy
scipy

一次性安装：

pip install -r grating-skill/requirements.txt

🚀 快速开始

1. 配置仿真参数

编辑 grating-skill/example_config.json，或修改模板脚本顶部的 CFG 字典：

{
  "lambda_nm": 1053.0,           // 真空波长 (nm)
  "theta_deg": 0.0,              // 入射角 (度，相对表面法线)

  "N_teeth": 1,                  // 光栅周期数 (1 = 标准无限周期 RCWA 单元)
  "period_nm": 500.0,            // 光栅周期 (nm)
  "duty_cycle": 0.5,             // 占空比 (齿宽/周期)
  "d_groove_nm": 200.0,          // 齿高 (nm)
  "n_tooth": 2.0,                // 齿折射率
  "n_groove": 1.0,               // 槽折射率 (通常为空气)

  "layer_stack": [               // 层栈：从底部（基底）向上
    [1000.0, 1.46, "substrate"]  // [厚度(nm), 折射率, 标签]
  ],

  "d_air_nm": 2106.0,            // 上方空气层厚度 (≥ 2λ 推荐)
  "out_prefix": "grating",       // 输出文件名前缀
  "save_mph": true               // 是否保存 .mph 模型文件
}

2. 运行仿真

cd grating-skill
python comsol_grating_template.py --config example_config.json

Windows 用户可双击 example_run.bat。

单次仿真耗时约 30 秒 ~ 5 分钟，取决于 N_teeth 和层数。

📤 输出文件

以 out_prefix="grating", N_teeth=1, theta_deg=0 为例：

文件	大小	用途
`grating_N1_theta0.00.mph`	~10-25 MB	已求解的 COMSOL 模型（可打开查看场分布）
`grating_N1_theta0.00.mat`	~50 KB	仿效率 + 各级次谱（MATLAB/Python 可读）
`grating_N1_theta0.00_orders.csv`	~1 KB	各级次 R/T 纯文本表格

控制台输出示例

[INFO] === Grating: N_teeth=1, theta=0 deg, lambda=1053 nm, W=500 nm ===
[INFO] propagating orders: top air=3 (-1..1), bot substrate n=1.46: 5 (-2..2)
[INFO] total ports: 8
[INFO] mesh elements: 12384
[INFO] solving ...
[INFO] solved.
[RESULT] R=0.21345  T=0.78655  R+T=1.00000
[INFO] top-3 reflected orders:
        m=  0  eta=0.21345  sin_t=+0.0000
        m= +1  eta=...
[CHECK] Poynting net flux: top=-0.21344 W/m, bot=-0.78656 W/m

⚠️ 重要：对于无损材料，R+T 必须为 1.0（精度约 5 位）。如果偏差超过 0.005，请先诊断问题（检查网格或几何），不要直接使用结果。

📊 结果读取

MATLAB

data = load('grating_N1_theta0.00.mat');
sim = data.grating_simulation;

R = sim.efficiencies_port_S.R;
T = sim.efficiencies_port_S.T;

% 读取反射级次数据
m_top   = sim.spectrum_top_portS.m;
R_m     = sim.spectrum_top_portS.eta;
sin_out = sim.spectrum_top_portS.sin_theta_out;

% 按效率排序显示
[R_sorted, idx] = sort(R_m, 'descend');
disp(table(m_top(idx), sin_out(idx), R_sorted, ...
    'VariableNames', {'m','sin_theta_out','R_m'}));

Python

import scipy.io as sio
import numpy as np

d = sio.loadmat("grating_N1_theta0.00.mat",
                squeeze_me=True, struct_as_record=False)
sim = d["grating_simulation"]

print("R =", sim.efficiencies_port_S.R)
print("T =", sim.efficiencies_port_S.T)

top = sim.spectrum_top_portS
order = np.argsort(-top.eta)
for i in order[:5]:
    print(f"m={top.m[i]:+d}: R={top.eta[i]:.5f}, sin_t={top.sin_theta_out[i]:+.4f}")

CSV 文件

直接用 Excel 打开 *_orders.csv，即可查看各级次衍射效率：

side, m, sin_theta_out, eta
top,  0,  +0.000000,    2.134527e-01
top, +1,  +0.707107,    3.456789e-02
top, -1,  -0.707107,    3.456789e-02
...

🔍 验证检查清单

仿真完成后，请检查以下三项：

能量守恒：[RESULT] R+T 应为 1.0000 ± 0.001（无损材料）。若偏离过大，网格太粗或几何有误。
Poynting 交叉验证：[CHECK] Poynting net flux 的 top 值应为负数（能量向下传播），bot 值应在 ~1% 内匹配。
几何可视化：在 COMSOL GUI 中打开 .mph 文件，确认几何结构正确。

📐 .mat 文件结构

grating_simulation
├── method                    "Periodic_Port_full_field_TE"
├── lambda_nm                 标量 [nm]
├── theta_i_deg               标量 [deg]
├── N_teeth                   整数
├── period_nm                 标量 [nm]
├── W_nm                      仿真单元宽度 [nm] = N_teeth × period
├── n_air                     标量
├── n_substrate               标量（底层折射率）
├── efficiencies_port_S       struct {R, T, total}      ← 可信赖数据
├── spectrum_top_portS        struct {m[], eta[], sin_theta_out[]}
├── spectrum_bot_portS        struct {m[], eta[], sin_theta_out[]}
├── poynting_check            struct {top_y_nm, top_Poavy, bot_y_nm, bot_Poavy}
└── field_Ez                  struct {x_nm[nx], y_nm[ny], Ez[ny×nx] 复数,
                                       abs_Ez[ny×nx], note}   （关闭导出时为空 struct）

其中 spectrum_*_portS.eta 是可信赖的各级次效率（来自 Port-S 数据，而非 FFT）。

对于每个传播级次 m：

m[i] = 超胞衍射级次索引（整数）
eta[i] = |S_{port_m, port_excited}|² = R_m 或 T_m
sin_theta_out[i] = sin(θ_i) + m·λ/W（介质中的传播常数比 kₓ_m/k₀）

Σeta 应等于 efficiencies_port_S.{R 或 T}。

Ez 场（复数）

field_Ez 是在整个仿真单元上规则网格采样的复数面外场 Ez(x, y)（频域相量，含实部/虚部，另附 abs_Ez 方便取模）：

x_nm[j] = 采样点 x 坐标 [nm]，范围 0 .. W
y_nm[i] = 采样点 y 坐标 [nm]，范围 0 .. (Σ各层厚 + d_groove + d_air)
Ez[i, j] = (x_nm[j], y_nm[i]) 处的复数 Ez [V/m]，行索引 y、列索引 x
落在几何之外的网格点为 NaN

网格分辨率由配置项 Ez_nx（x 采样数，默认 201）和 Ez_ny（y 采样数；0 = 自动，约方形像素）控制。设 export_Ez: false 可跳过采样，此时 field_Ez 保存为空 struct。

Ez = sim.field_Ez.Ez;   % 复数 (ny, nx)
imagesc(sim.field_Ez.x_nm, sim.field_Ez.y_nm, abs(Ez)); axis xy equal tight; colorbar

🧩 配置示例

单周期 RCWA 风格光栅（无多层）

{
  "lambda_nm": 1053,
  "theta_deg": 0,
  "N_teeth": 1,
  "period_nm": 500,
  "duty_cycle": 0.5,
  "d_groove_nm": 200,
  "n_tooth": 2.0,
  "n_groove": 1.0,
  "layer_stack": [[1000.0, 1.5, "substrate"]],
  "d_air_nm": 2106
}

玻璃基底减反射光栅（含匹配层）

{
  "lambda_nm": 1550,
  "theta_deg": 30,
  "N_teeth": 1,
  "period_nm": 700,
  "duty_cycle": 0.4,
  "d_groove_nm": 350,
  "n_tooth": 2.05,
  "n_groove": 1.0,
  "layer_stack": [
    [1000.0, 1.45, "substrate"],
    [120.0,  1.62, "matching"]
  ],
  "d_air_nm": 3100
}

多齿超胞（有限宽度边缘效应研究）

{
  "lambda_nm": 1053,
  "theta_deg": 66.36,
  "N_teeth": 13,
  "period_nm": 574,
  "duty_cycle": 0.30,
  "d_groove_nm": 380,
  "n_tooth": 1.46,
  "n_groove": 1.0,
  "layer_stack": [
    [1000.0, 1.46, "substrate"],
    [464.38, 1.46, "matching"]
  ],
  "d_air_nm": 2106
}

对于真实的 DBR 多层膜，只需在 layer_stack 中依次列出每层 HfO₂/SiO₂ 即可。层标签必须为合法的 Python 标识符（无空格、无特殊字符）。

⚠️ 关键 API 陷阱

本模板防御了 9 个已知的 COMSOL mph API 陷阱。最常见且容易被"修复"的问题：

1. 不要使用 `property()` 方法

# ❌ 错误：PhysicsClient 没有 .property() 方法，静默失败
ewfd.property("Components", "OutofPlane")

# ✅ 正确：使用 .prop()
ewfd.prop("components").set("components", "outofplane")  # 注意全小写！

2. 不要遗漏 Port 的 E0 设置

# ❌ 错误：默认 E0="t1x" 会导致 TM 偏振！
p.set("LinearPol", "S")  # 这个设置无法覆盖 E0

# ✅ 正确：显式设置 E 沿 z 方向（TE）
p.set("E0", ["0", "0", "1"])

3. 不要将 lengthUnit 改为 "nm"

# ❌ 错误：会导致 kx×x 表达式的单位不匹配（偏差 10⁹ 倍！）
geom.lengthUnit("nm")

# ✅ 正确：使用 SI 单位，长度值加 [nm] 后缀
geom.lengthUnit("m")
P.set("period", "574[nm]")

4. 不要省略 FreeTri 网格特征

# ❌ 错误：省略后网格元素数为 0，求解静默返回 NaN

# ✅ 正确：显式添加 FreeTri
ft = mesh.feature().create("ftri1", "FreeTri")

完整列表请参阅 BUGS_AND_FIXES.md。修改模板前请先阅读。

🤝 贡献

欢迎提交 Issue 和 Pull Request！

📄 许可证

MIT License

📚 参考资料

BUGS_AND_FIXES.md — 完整的 API 陷阱列表及修复方案
SKILL.md — Agent 使用规范（输入/输出格式详细说明）
COMSOL_MATLAB_Pitfalls_Guide.md — 历史笔记
COMSOL Wave Optics Module 文档

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

COMSOL 光栅衍射仿真工具

✨ 功能特性

🚫 不适用场景

📁 项目结构

🔧 环境要求

软件依赖

Python 包

🚀 快速开始

1. 配置仿真参数

2. 运行仿真

📤 输出文件

控制台输出示例

📊 结果读取

MATLAB

Python

CSV 文件

🔍 验证检查清单

📐 .mat 文件结构

Ez 场（复数）

🧩 配置示例

单周期 RCWA 风格光栅（无多层）

玻璃基底减反射光栅（含匹配层）

多齿超胞（有限宽度边缘效应研究）

⚠️ 关键 API 陷阱

1. 不要使用 `property()` 方法

2. 不要遗漏 Port 的 E0 设置

3. 不要将 lengthUnit 改为 "nm"

4. 不要省略 FreeTri 网格特征

🤝 贡献

📄 许可证

📚 参考资料

FilesExpand file tree

readme_zh.md

Latest commit

History

readme_zh.md

File metadata and controls

COMSOL 光栅衍射仿真工具

✨ 功能特性

🚫 不适用场景

📁 项目结构

🔧 环境要求

软件依赖

Python 包

🚀 快速开始

1. 配置仿真参数

2. 运行仿真

📤 输出文件

控制台输出示例

📊 结果读取

MATLAB

Python

CSV 文件

🔍 验证检查清单

📐 .mat 文件结构

Ez 场（复数）

🧩 配置示例

单周期 RCWA 风格光栅（无多层）

玻璃基底减反射光栅（含匹配层）

多齿超胞（有限宽度边缘效应研究）

⚠️ 关键 API 陷阱

1. 不要使用 property() 方法

2. 不要遗漏 Port 的 E0 设置

3. 不要将 lengthUnit 改为 "nm"

4. 不要省略 FreeTri 网格特征

🤝 贡献

📄 许可证

📚 参考资料

1. 不要使用 `property()` 方法