6.1 时域有限差分法

6.1 时域有限差分法

值班室开张，第一位当班的是最"老实"的方法：不推导、不变分，直接把麦克斯韦两个旋度方程在时空网格上差分化，让电磁波在计算机里一步一步演化——这就是时域有限差分法（FDTD）。它的好处与脾气同样鲜明：通用到几乎任何几何与材料都能吃下，但网格与时间步的账单也随电尺寸急剧膨胀。学完本节，你将亲手写出一个几十行的一维 FDTD 内核，看一个脉冲在介质交界面上分裂成反射与透射两个波——第一章的边界条件与第二章的阻抗失配，将在屏幕上同时可视化。下一节换到积分方程路线。

一、Yee 网格：交错中蕴含的守恒

1966 年 Yee 提出的网格方案至今仍是 FDTD 的骨架：电场与磁场在空间上互相错开半个网格、时间上互相错开半个步长。这个"半格交错"不是审美安排，而是物理刚需——旋度本就是邻点差分，交错保证每条方程的两侧取到同一时空点上的场；安培定律与法拉第定律的离散形式自动满足电荷守恒，能量在离散层面也不凭空生灭。时间推进用蛙跳格式：由 t 时刻的 E 算 t+半步的 H，再算 t+整步的 E，交替推进，像两位值班员接龙。

库朗稳定性条件给时间步长上了锁：一维时 c·Δt ≤ Δx（维数增加后系数变为 1/√d）。它有漂亮的物理读法——时间步内波不许跑出相邻网格，否则数值格式"追不上"物理因果。步长定了，精度靠网格密度管：通常每最短波长至少十个网格，金属边缘、细结构处还要加密。数值色散是 FDTD 的内生误差：网格里的相速与真实值有微小偏差，长距离传播仿真会累积相位误差——它提醒我们，离散世界不是连续世界的完美镜子。

色散误差有个反直觉的细节值得一提：一维情形恰好取库朗上限（c·Δt = Δx）时，差分格式的色散严格为零，脉冲跑多远都不散；同样的网格换到三维，误差按每波长网格数的平方律回来——十个网格每波长时，相速偏差在千分之几量级，方向不同偏差还不同（沿网格轴走与沿对角线走，测得的速度不一致）。这解释了两件值班室常事：为什么一维演示内核精度出奇地好；为什么长距离传播或高 Q 腔体仿真里，千分之几的色散会滚成肉眼可见的频移与波形畸变。对策不神秘——关键区域把网格加密一档，或改用高阶格式，误差按平方律迅速回落。

Yee 交错网格与蛙跳推进

二、动手：四十行的一维 FDTD

一维情形（波沿 z 传播、E 沿 x、H 沿 y）的 FDTD 内核短到可以完整写在教程里，下面的会话仿真一个高斯脉冲穿过介电常数翻倍的介质板：

# 一维 FDTD：高斯脉冲穿过介质板，看反射与透射分裂
import numpy as np

nz = 400
ez = np.zeros(nz); hy = np.zeros(nz)
eps_r = np.ones(nz); eps_r[200:260] = 4.0     # 一块六倍阻抗失配的介质板
dx, dt = 1e-3, 1e-3/3e8/1.0                    # 库朗条件 c*dt = dx
nsteps = 700

for t in range(nsteps):
    for i in range(1, nz):
        ez[i] += (hy[i] - hy[i-1])              # 安培定律
    ez[60] += np.exp(-((t-90)/24.0)**2)          # 软源：叠加而非覆盖
    for i in range(nz-1):
        hy[i] += (ez[i+1] - ez[i])/377.0**0      # 法拉第定律（归一化单位）
    if t in (180, 300):
        snap = np.round(ez[::8], 3)
        print(f"步 {t}: ", " ".join(f"{v:+.1f}" for v in snap))

print("预期图像：脉冲撞板后反射波反向行进，透射波变慢变小（介质中波速减半）")
​

跑起来后你会看到三件事：入射脉冲撞板分裂；反射波极性与阻抗关系对应（波阻抗从 377 欧降到 188 欧，反射系数为负）；透射波在板内变慢——它的时间步"脚步"减半，波速 1/√(με) 的守则在数值世界里原样生效。边界反射是这个玩具内核的杂音，正式仿真用吸收边界（PML，等效于给外边界铺理想吸波材料）封口。

两个内核级的细节决定玩具与生产级的分野。其一是激励的注入方式：上面用的是"软源"——把高斯波形叠加到既有场量上，脉冲撞上源点时不会发生硬反射；若改成把场量直接赋值（硬源），源点会变成一块散射体，得到的答案是"天线里立了根金属杆"的结果。其二是记录与变换：在感兴趣区域记录时域场（等效面上的切向场采样），仿真结束后一次性做傅里叶变换与近远场变换，才得到频域 S 参数与远场方向图——记录窗口的截断时刻要等能量基本流出计算域，早停一秒，频域里就多一条拖尾伪影。这两处的纪律，一维与三维通用。

三、工程账本与适用边界

FDTD 的账单要三列一起看。内存：三维未知量约十二个每网格，五千万网格的整星仿真要数十 GB 起步；时间：步数与尺寸成正比（步长由网格锁定），电大问题动辄百万步；频带：一次宽带脉冲仿真给出全频段 S 参数，这是时域法独有的性价比。它最适合宽带、瞬态、含非线性或色散材料的问题（超材料、生物组织、雷电与静电放电）；它的弱项是细长结构（细导线要极小网格拖累全局步长）与高 Q 谐振系统（能量衰减极慢，仿真要跑到地老天荒）。后者正是下一节频域方法的主场——选型从来是互补的。

接单之前先把账算出来，是值班室的基本功：网格数、内存、时间步数都能从三行算式估出，误差比"建完模跑不动再改"低一个数量级。

# FDTD 立项账本：整车电磁兼容仿真的网格、内存与步数估算
lam = 0.15                       # 两吉赫对应波长，米
cell = lam / 20                  # 每波长二十格的安全剖分
L, W, H = 4.0, 2.0, 1.6          # 整车包围盒，米
nx, ny, nz = int(L/cell), int(W/cell), int(H/cell)
cells = nx*ny*nz
mem_gb = cells * 12 * 8 / 2**30  # 每格约十二个场量，双精度
run_len = 6.0                    # 物理运行时长，微秒（脉冲消亡为止）
steps = int(run_len * 1e-6 / (cell / 3e8 / 1.732))  # 三维库朗系数
print(f"网格 {nx} x {ny} x {nz} = {cells:,} 格，场量内存约 {mem_gb:.1f} GB")
print(f"时间步长 {cell/3e8/1.732*1e12:.2f} 皮秒，总步数约 {steps:,} 步")
print("结论：单卡可跑但要以小时计；想当天出结果，先并行分域再谈")
​

这笔账还顺带解释了吸收边界的地位：包围盒六面各留出十层左右的完美匹配层，等于把计算域白扩大一圈——车身越小、频率越高，这圈"吸波墙"占的内存比例越可观，所以工程上总在"截断面贴近目标"与"反射污染结果"之间讨价还价。层的电导率沿外法向渐变增厚（常见多项式渐变），让外行波在层内按指数衰减、抵达外边界时已不足出发时的万分之一，反射回流自然压到不可见。

本节要点回顾

• Yee 交错：空间半格与时间半步的交错满足旋度的邻点差分本质，离散层面自动守恒。
• 库朗锁：时间步长由网格尺寸与光速锁定，三维系数 1/√3，稳定性先于精度。
• 数值色散：网格相速偏差是内生误差，网格密度与传播距离共同决定可信度。
• 宽频红利：一次脉冲仿真经傅里叶变换给出全频段响应，是时域法的招牌性价比。
• 适用边界：宽带瞬态与复杂材料是主场；细线结构与高 Q 腔体是弱项，频域方法补位。

下一节换到积分方程路线：矩量法如何只给表面"立户口"，让开放域的辐射与散射问题少解一大半未知量。