3、FDTD方法入门：时域有限差分法原理、Yee网格、Courant稳定性条件、边界条件(PML/Periodic)
各位同学，今天我们来聊聊FDTD。全称是时域有限差分法，英文叫Finite-Difference Time-Domain。说实话，我刚入行那会儿，觉得这名字挺唬人的。但用久了你会发现，它其实就是把麦克斯韦方程组在时间和空间上拆成小块，然后一步一步往前推。
我个人习惯把FDTD比作「拍电影」。你把电磁波传播的过程切成无数帧，每一帧算一算电场和磁场长什么样，然后下一帧接着算。只要帧率够高，你就能看到完整的物理过程。嗯，就是这么回事。
3.1 FDTD的核心思想
FDTD不搞什么频域变换，它直接在时域里干活。你给一个光源，比如一个脉冲，它就在计算区域里跑。跑着跑着，碰到你的MicroLED结构，就会反射、透射、散射。FDTD把这些过程全部记录下来。
为什么要这么做？因为MicroLED的光学行为很复杂。你想想看，一个几微米大小的LED，里面有量子阱、布拉格反射镜、各种介质层。光在里面来回反射、干涉、衍射。用解析方法算？头都大了。FDTD直接暴力求解，反而简单。
核心要点：FDTD把麦克斯韦旋度方程中的偏微分，用中心差分近似替代。时间上一步步推进，空间上一个个网格计算。
3.2 Yee网格——这个网格有点意思
1966年，有个叫Yee的美国人提出了一个网格方案。这个方案有多经典？到现在快60年了，几乎所有FDTD软件还在用。
Yee网格的精髓在于：电场和磁场不在同一个点上。它们错开半个网格步长。为什么这么干？因为麦克斯韦方程里，电场的旋度产生磁场，磁场的旋度产生电场。你让它们错开，正好可以用中心差分来近似旋度。
我举个例子你就明白了。假设你有一个三维空间，每个网格点我们放一个电场分量。那磁场分量放在哪里？放在网格的棱边上。具体来说：
- 电场E的分量放在网格棱边的中心
- 磁场H的分量放在网格面心的中心
- 时间上，E和H也错开半个时间步
这样做的好处是什么？自动满足法拉第定律和安培定律的积分形式。说白了，就是物理上自洽。我在项目中遇到过有人自己改网格方案，结果算出来的东西不守恒，能量越算越多。嗯，这就是没理解Yee网格的精髓。
我的经验：初学者最容易犯的错误是搞混E和H的位置。画个图贴在显示器旁边，很有帮助。
3.3 Courant稳定性条件——别让计算炸掉
FDTD是显式算法，什么意思？就是每一步的计算只依赖上一步的结果。这带来一个问题：如果时间步长太大，数值误差会指数级增长，最后结果全是NaN。
Courant稳定性条件就是告诉你：时间步长不能超过某个上限。公式长这样：
Δt ≤ 1 / (c * sqrt(1/Δx² + 1/Δy² + 1/Δz²))其中c是光速，Δx、Δy、Δz是空间步长。对于三维均匀网格，可以简化为：
Δt ≤ Δx / (c * √3)这个条件叫Courant-Friedrichs-Lewy条件，简称CFL条件。我建议你取安全系数0.9，也就是算出来的上限再打九折。为什么？因为我在项目中吃过亏，取0.99的时候，跑了十万步才炸，浪费了整整两天计算时间。
警告：如果你用非均匀网格，Courant条件会更严格。网格最密的地方决定了时间步长。我曾经见过有人只在局部加密网格，忘了调整时间步长，结果算到一半直接崩了。
3.4 边界条件——让波「消失」的艺术
FDTD的计算区域是有限的，但电磁波传播是无限的。你总不能让波碰到边界就反射回来吧？那算出来的结果全是驻波，跟实际情况差远了。
3.4.1 PML边界——完美匹配层
PML，全称Perfectly Matched Layer，完美匹配层。这是1994年Berenger提出来的，堪称FDTD领域最重要的突破之一。
PML的原理是什么？在计算区域的外围加一层「吸收材料」。这层材料经过特殊设计，波进去之后会被指数衰减，几乎不产生反射。我打个比方，就像你走进一个消声室，声音进去就出不来了。
实际使用中，PML的厚度一般是10到20个网格。太薄了吸收不干净，太厚了浪费计算资源。我个人习惯用16层，效果不错。
避坑指南：我曾经遇到过PML参数设置不对，导致低频波反射严重。后来发现是PML的导电率分布函数没选对。建议用多项式渐变，阶数取3或4。
3.4.2 Periodic边界——模拟无限周期结构
MicroLED经常做成阵列，比如一个4×4的微透镜阵列。你不可能把整个阵列都仿真，计算量太大了。这时候就用Periodic边界条件。
Periodic边界的意思是：左边出去的波，从右边进来。这样你只需要仿真一个周期单元，就能得到整个阵列的行为。但要注意，Periodic边界只适用于垂直入射或者特定角度入射的情况。斜入射时要用Bloch边界条件，那是Periodic的推广。
我做过一个项目，用Periodic边界仿真MicroLED的出光效率。一个周期单元只有5微米×5微米，算出来的结果跟实测误差在3%以内。这就是Periodic边界的威力。
3.5 知识体系总览
下面这张图是我自己画的，把FDTD的核心要素串起来了。你仔细看看，应该能对整体有个把握。
3.6 一个简单的FDTD代码框架
说了这么多理论，我们来点实际的。下面是一个二维FDTD的核心循环代码，用Python写的。你不需要完全理解每一行，但看看整体结构，能帮你建立直觉。
# 二维TM波FDTD核心循环
# Ez, Hx, Hy 分别代表电场和磁场分量
# 初始化
Ez = np.zeros((Nx, Ny))
Hx = np.zeros((Nx, Ny))
Hy = np.zeros((Nx, Ny))
# 时间步进
for n in range(max_steps):
# 更新磁场Hx
Hx[:-1, :] += dt / (mu0 * dy) * (Ez[:-1, :] - Ez[1:, :])
# 更新磁场Hy
Hy[:, :-1] += dt / (mu0 * dx) * (Ez[:, 1:] - Ez[:, :-1])
# 更新电场Ez
Ez[1:-1, 1:-1] += dt / (eps0 * dx) * (Hy[1:-1, 1:-1] - Hy[1:-1, :-2]) \
- dt / (eps0 * dy) * (Hx[1:-1, 1:-1] - Hx[:-2, 1:-1])
# 添加光源
Ez[source_x, source_y] += source_pulse(n * dt)
# 应用PML边界（这里省略具体实现）
apply_pml(Ez, Hx, Hy)
# 记录结果
if n % record_step == 0:
save_field(Ez, n)你看，核心循环就这么几行。但就是这几行代码，能模拟从微波到光波的几乎所有电磁问题。我当年第一次跑通这个循环，看到电磁波在屏幕上传播的时候，说实话，挺激动的。
重要提醒：上面的代码只是教学演示。实际工程中，你需要考虑材料色散、各向异性、非线性效应等。但这些都是在核心循环基础上加「插件」。
3.7 我的几点建议
最后，给刚开始学FDTD的同学几点建议：
- 先跑通一个简单例子。比如真空中的点光源传播，看看波前是不是圆的。如果连这个都跑不对，后面别想了。
- 网格尺寸要合适。一般取波长的1/10到1/20。太粗了精度不够，太细了算不动。MicroLED的可见光波长400-700nm，网格取20-40nm比较合适。
- PML边界一定要用。别想着计算区域够大就不用PML。我试过，反射回来的波会污染你的结果，尤其是远场计算。
- 先做收敛性测试。把网格加密一倍，看看结果变不变。如果变化很大，说明你的网格还不够密。
好了，FDTD入门就讲到这里。内容不少，但都是干货。你回去之后，建议动手写一个最简单的二维FDTD程序，哪怕只算一个点光源也行。亲手跑起来的感觉，跟看书完全不一样。