企业官网建设流程全解析

简介：这是一份面向光子晶体研究者与Comsol/MATLAB用户的陈数计算专题资料，内容聚焦光子晶体能带结构模拟与拓扑陈数计算，解决从Comsol建模到MATLAB数据后处理的核心问题。压缩包共5个文件，内含两篇PDF说明文档、一个mph模型源文件、一段m计算脚本以及txt配套文本，总体约20MB；PDF侧重原理解读与文件说明，mph与m代码则提供可运行的三维建模和能带提取示例。已有1043人学习下载。资料兼顾概念与实操：既讲解光子晶体周期性结构、光子禁带、布里渊区划分等基础，也演示利用频域仿真求解麦克斯韦方程获取电磁场分布，并通过LiveLink将结果导入MATLAB绘制能带图，最终用于陈数分析。对设计光子晶体光纤、光子晶体激光器等器件，或从事拓扑光子学研究的读者具有直接参考价值。

1. 陈数计算分享：这个拓扑整数到底怎么从Comsol里“写”出来

陈数（Chern number）在拓扑能带理论里算是最重要的一个整数：二维周期结构每一条能带都带一个 C∈Z，C 不为 0 意味着边界上必然出现单向输运的边界态。Comsol 用户把能带图扫得很漂亮，卡就卡在“陈数的数怎么写”——翻遍后处理菜单也找不到一个叫 Berry 曲率的变量。这篇文章就把两件事讲透：用 Comsol 把布洛赫波函数导出来，再在后处理脚本里用格点规范场把陈数算出来；同时解释为什么能带打开、带隙清晰，陈数照样可能算成 0.5。适合已经在 Comsol 里跑过周期结构本征模、现在想算拓扑数和边界态的工程师、研究生。下面按我实际跑通的一条路径来写，所有参数都给出可以直接改的值。

2. 把陈数拆成能上机的公式：先从 Berry 曲率写到格点相位求和

2.1 陈数的标准写法与为什么不能直接在 Comsol 后处理里积分

布洛赫态写成：

ψ_{n,k}(r) = e^{ik·r} u_{n,k}(r)

其中u是周期因子，注意它跟 Comsol 直接解出的场ψ差了一个平面波因子。Berry 曲率的标准写法是：

Ω_n(k) = i ( ⟨∂_kx u | ∂_ky u⟩ − ⟨∂_ky u | ∂_kx u⟩ )

陈数就是 Berry 曲率在布里渊区上的积分：

C_n = (1/2π) ∬_{BZ} Ω_n(k) dk_x dk_y

这个式子看着简洁，直接拿到 Comsol 里却写不出来，原因有三个。第一，Comsol 解出的布洛赫波函数相位是不固定的，每次重新求解都会乘一个任意的e^{iχ(k)}；直接用数值微分算∂_k u，微商会被人工相位噪声淹没。第二，你需要在相邻k点拿到同一物理能带上的波函数，但特征值求解器默认按频率排序，简并点附近很容易串带。第三，有限元节点值的内积并不是简单的逐点相乘，后续会专门讲。

所以业内做陈数几乎都用“格点规范场”写法，也就是 Fukui–Hatsugai–Suzuki 方法：把布里渊区切成N_x × N_y个小格子，每个格子边上定义：

U1(i,j) = ⟨u(i,j)|u(i+1,j)⟩ / |⟨u(i,j)|u(i+1,j)⟩|

U2(i,j) = ⟨u(i,j)|u(i,j+1)⟩ / |⟨u(i,j)|u(i,j+1)⟩|

然后陈数写成：

C = (1/2πi) Σ_ij ln[ U1(i,j) · U2(i+1,j) · conj(U1(i,j+1)) · conj(U2(i,j)) ]

这才是“陈数的数”在后处理代码里的标准写法：四个 U 围成一个小闭合回路，取整体辐角，把所有小格子的辐角累加再除以2π。直接用这段公式，避开了 Berry 曲率里所有的规范依赖问题。

2.2 Comsol 四个必须对齐的配置：Floquet 边界、复值解、扫描顺序和导出列

要把上面的公式跑起来，Comsol 模型里有四个配置不能错，错了后面代码再漂亮也算不出整数。

扫描网格我一般取Nx = Ny = 24起步。kx = -π/a + (0:Nx-1) * 2π/(a*Nx)，意思是最后一个点不取π/a，让首尾天然闭合。如果你用linspace(-π/a,π/a,Nx)那样的写法，布里渊区边界会被重复采样，后面mod周期索引会出现错位。

特征值求解时要多求两条能带。只求目标能带很容易搅进非物理的零模或数值模式；多算两条带，做模式连接时才有余量。对于光子晶体、声子晶体，物理场接口本身自带复数求解；如果用 PDE 模式自己写方程，记得确认“允许复数本征值”没有被关掉。这些配置对齐后，求解器会给你一组带kx/ky/mode/freq/u/im(u)的记录，下一步就是导数据。

3. 把 Comsol 波函数倒进 MATLAB：从导出列到陈数积分的最小脚本

3.1 导出文件的行列布局：一列复数场分量怎么摆

在 Comsol 的“导出”节点里，把数据集指到特征值求解的解。导出的文本通常是这种结构：

# kx ky mode freq re_u im_u 0.0000 0.0000 1 2.1345e11 1.0243e-1 0.0000e0 0.0000 0.0000 2 3.2011e11 3.2210e-2 4.1103e-1 ... 0.2618 0.0000 1 2.1410e11 5.1110e-2 8.0030e-2

同一组(kx, ky)会有mode=1,2,3,...多行，每一行是该波矢下的一个本征模。这里最容易犯的错是把mode列当成能带编号：它只是求解器输出顺序，不是物理能带顺序。物理上同一条能带要靠后面讲的模式连接来重新排。

读取时按kx、ky重排成数组。以下代码假设数据列依次是 kx、ky、mode、freq、实部、虚部：

% 读取Comsol导出的文本，跳过1行表头 raw = readmatrix('chern_data.txt', 'NumHeaderLines', 1); kx = raw(:,1); ky = raw(:,2); mode = raw(:,3); u = raw(:,5) + 1i*raw(:,6); % 拼成复数波函数 band = 1; % 先取第1条物理能带 sel = mode == band; % 这一句只在模式已经对齐时可用 kxb = kx(sel); kyb = ky(sel); ub = u(sel); kx_list = unique(kxb); ky_list = unique(kyb); Nx = numel(kx_list); Ny = numel(ky_list); Psi = zeros(Nx, Ny, 1); % 先存单分量波函数 for q = 1:numel(kxb) ix = find(kx_list == kxb(q), 1); iy = find(ky_list == kyb(q), 1); Psi(ix, iy, 1) = ub(q); end

如果你是弹性波、压电或多分量 PDE，就把每个分量的实部/虚部都导出来，在循环里拼成Psi(ix,iy,:)的一个复数向量再传入后面的函数。unique要求kx/ky网格是规则的，如果扫描里混入了失败解，数组重排前先检查numel(kx_list)*numel(ky_list)是否等于实际行数。

3.2 Fukui 格子规范：把陈数的数用四行相位凑出来

重排好波函数后，主菜就是一页纸的陈数函数。这个函数不计算 Berry 曲率，只计算相邻k点之间的 U(1) 相位链接：

function C = chern_from_comsol(Psi) % Psi: 复数波函数数组 size=(Nx, Ny, Ndof) % 输出: 整数陈数 Nx = size(Psi, 1); Ny = size(Psi, 2); ix_next = [(2:Nx), 1]; % kx 方向周期相邻 iy_next = [(2:Ny), 1]; % ky 方向周期相邻 U1 = zeros(Nx, Ny); U2 = zeros(Nx, Ny); for i = 1:Nx for j = 1:Ny v0 = Psi(i, j, :); v0 = v0(:); v1 = Psi(ix_next(i), j, :); v1 = v1(:); v2 = Psi(i, iy_next(j), :); v2 = v2(:); O1 = v0' * v1; % 复内积 O2 = v0' * v2; U1(i, j) = O1 / abs(O1); % 只保留相位，丢掉模长 U2(i, j) = O2 / abs(O2); end end F = U1 .* U2(ix_next, :) .* conj(U1(:, iy_next)) .* conj(U2); C = round( sum(angle(F), 'all') / (2*pi) ); end

这段代码的逻辑就三句话。第一，U1/U2是相邻波矢上波函数的内积再归一化，它把求解器随机乘的e^{iχ}消掉了；这就是为什么不需要手动固定规范。第二，F把四个方向链接合成一个闭合回路，angle(F)是每格 Berry 相位。第三，全部小格相位累加后除以2π，round只负责清掉1e-13量级的数值尾巴。

调用时直接说C = chern_from_comsol(Psi)。如果算出来是0.9998，取整后1，正常；如果算出来是0.63或1.35，别急着 round，先去看是不是模式串了、网格不够细、或者内积近似出了问题。

4. 陈数算成 0.5 的翻车现场：五个让结果变成废数的坑

4.1 能带在交叉处互换，积分给出 0.63

现象：同一个模型跑了三遍，陈数分别是 0.63、1.02、0.5，怎么改参数都不稳定。

原因：Comsol 特征值求解器按频率大小返回本征模。能带交叉或反交叉时，mode序号在某个k点突然跳变，你取的第 1、2、3 条“带”已经不是同一个物理能带。

解决：在算法里做模式连接。具体做法是用上一个k点的波函数跟当前点所有候选模式做内积，取模最大的候选作为同一条带的延续。模式连接做完，陈数通常立刻回到整数。这也是我不建议按mode==band直接筛数据的原因，第 5 章会给可操作步骤。

4.2 导出只带实部，陈数恒等于 0

现象：所有U1/U2都等于 1，angle(F)全为 0，C 永远为 0。

原因：Comsol 默认后处理视图只显示实数部分，导出节点如果没有显式加imag(u)，你拿到的是一组纯实波函数。相位信息全丢，陈数自然写不出来。

解决：导出表达式中单独加一列imag(u)。用 PDE 多分量时，记得是imag(u1)、imag(u2)，不是imag(u)一个量。检查导出文件里波形是否带明显的复数虚部。

4.3 布里渊区两端点重合，陈数对 N 不稳定

现象：Nx=20时 C=0.87，Nx=30时 C=0.92，算出来的数离整数越来越近但就是不落在整点上。

原因：扫描范围写成了[-π/a, π/a]闭区间。边界点-π/a和π/a在布里渊区里是同一个点，等于是重复采样。重复点破坏周期相邻关系，整个闭合回路被多算或少算了一个格点。

解决：扫描区间用半开区间。比如kx = -pi/a + (0:Nx-1) * 2*pi/(a*Nx)。代码里的ix_next=[(2:Nx) 1]已经默认了周期闭合，扫的时候别再把同一个点扫两次。

4.4 高对称点简并，相位因子乱跳

现象：C 在 0 和 1 之间反复横跳，angle(F)在某个局部区域出现类似棋盘格的乱值。

原因：在 Γ、K、M 这类高对称点上，两条能带简并。简并子空间里的本征矢可以任意旋转，单条能带的 Berry 相位在简并点附近根本不唯一，Fukui 公式在简并点上给不出稳定结果。

解决：把网格整体偏移半个格点，避开高对称点。具体做法是kx = -pi/a + (i-0.5) * 2*pi/(a*Nx)，让采样点落在高对称点之间的位置。如果简并不是孤立的而是整条能带简并，单能带陈数本身没有定义，需要做多能带的 Wilson loop，那不是几行脚本能解决的问题。

4.5 节点值直乘当内积：网格一改符号就翻

现象：用一样的网格算某个模型 C=+1，加密网格后符号变 −1，甚至出现 0.5 附近的数。

原因：有限元的自由度是节点值，严格的两个波函数内积应该是u1ᵀ M u2，其中 M 是质量矩阵。直接用节点值逐点相乘等于把 M 当成单位阵。网格很细时逼近得还可以，但粗网格或者波函数剧烈振荡时会出问题。

解决：有三个常见做法。第一，用 LiveLink for MATLAB 从 Comsol 里拿质量矩阵做完整内积；第二，把解插值到规则采样网格，再算逐点内积，能减小质量矩阵带来的偏差；第三，先做收敛性检查，在 Nx=24、32、40 三种网格上各跑一次陈数，如果对 N 不敏感就可以接受。陈数本身是拓扑整数，抗噪声能力很强，但如果出现基础积分权重错误，收敛性检查是能揪出来的。

5. 把结果钉死：用模式连接和边界态验证陈数，而不是只看整数

5.1 模式连接：让能带编号在扫描里不换脸

拿到整数 C 之前，我建议先跑一遍模式连接，至少能确认你的“第 1 条带”从头到尾是同一条物理能带。做法很简单：对每个(i,j)格点，先拿上一个k点已经确认的波函数v_old，跟当前格点所有候选模式算内积：

G(m) = | ⟨v_old | v_m(k)⟩ |

模最大的m就是这条能带的延续。把所有格点都这样扫一遍，重排后的Psi再丢进chern_from_comsol。如果重排前后 C 不一致，说明之前mode列表确实串带了。

这个步骤的关键参数是内积阈值：G_max如果小于 0.5，说明当前格点已经不在同一物理带附近，要么扫描步长太大，要么存在简并。这时候加密kx/ky网格比怀疑公式更有效。高对称点附近建议用 0.5 格偏移的网格来处理。

5.2 边界态穿越体带隙的次数是一台“验钞机”

陈数是整数量子数，单看一个整数也可能算错符号。我一般的验证习惯是把同一几何改成一维带状结构：x 方向仍用 Floquet 条件扫kx，y 方向取 20 个周期并截断成有限宽，边界保留开放。算完特征频率后画f(kx)能带图。

这时你会看到体带围出的带隙里出现一条或几条边界态色散。边界态从下体带穿越带隙、连到上体带的次数，就是整条介质边界处的拓扑陈数。比如 C=1，边界态必须跨带隙一次；C=2，要跨两次。方向对不上，说明符号写反了，通常回去检查U1和U2的定义顺序即可。我还习惯画出每个小格子的angle(F)热图，陈数非 0 时热图上曲率会集中在几个动量点附近，而不是均匀铺开；看到热图就基本知道这个整数是从哪里挣来的。

这套流程我从最初算光子晶体、声子晶体一路用到现在，最大的教训就是：陈数不能靠肉眼看出能带边缘的绕向，也不能只看后处理里的某个“趋势”，必须把布洛赫波函数的相位老老实实导出来，用格点规范场凑成闭合回路。把模式连接和边界态验证两个习惯一起用上，即便不依赖别人现成的代码包，也能在 Comsol 里稳定复现别人文章里的陈数。希望帮到你。

本文还有配套的精品资源，点击获取