麦克斯韦方程组的差分格式推导
在电磁场数值分析中，直接求解连续形式的麦克斯韦方程组（Maxwell's Equations）是不可能的，因为计算机只能处理有限数量的离散数据。因此，将连续的偏微分方程（PDEs）转化为离散的代数方程组是数值计算的核心步骤。本文将重点介绍最经典的差分格式——时域有限差分法（FDTD）的核心思想，即基于 Yee 网格（Yee Cell）的麦克斯韦方程组差分格式推导过程。
在各向同性、线性、无损介质中，麦克斯韦方程组的微分形式可以表示为：
- 法拉第电磁感应定律：
$$\nabla \times \mathbf{E} = -\mu \frac{\partial \mathbf{H}}{\partial t}$$ - 安培环路定律：
$$\nabla \times \mathbf{H} = \epsilon \frac{\partial \mathbf{E}}{\partial t} + \mathbf{J}$$ - 高斯电场定律：
$$\nabla \cdot \mathbf{E} = \frac{\rho}{\epsilon}$$ - 高斯磁场定律：
$$\nabla \cdot \mathbf{H} = 0$$
在进行时域数值模拟时，我们主要关注前两个涉及时间导数和旋度的方程，因为它们描述了电磁波在空间和时间上的演化过程。
离散化策略：Yee 网格与交错时空
直接在规则的笛卡尔网格上对 $E$ 和 $H$ 进行中心差分会导致数值不稳定或无法准确表达旋度关系。为了解决这个问题，Yee 提出了一种“交错网格”（Staggered Grid）的概念。
空间交错
在 Yee 网格中，电场分量 $E_x, E_y, E_z$ 和磁场分量 $H_x, H_y, H_z$ 并不位于同一个空间点上。
- $E_x$ 分量位于网格边线上，其位置坐标为 $(i+1/2, j, k)$。
- $H_y$ 分量位于网格边线上，其位置坐标为 $(i, j+1/2, k)$。
这种空间上的错位使得在计算某个分量的旋度时，所需要的场分量恰好位于该点周围，从而实现了空间上的自然匹配。
时间交错
为了实现时间上的二阶精度，我们采用“蛙跳法”（Leapfrog method）：
- 电场 $\mathbf{E}$ 在半整数时刻 $n+1/2, n+3/2, \dots$ 进行更新。
- 磁场 $\mathbf{H}$ 在整数时刻 $n, n+1, n+2, \dots$ 进行更新。
这种时间上的交错使得在计算时间导数时，可以使用中心差分，从而获得二阶精度。
差分格式的具体推导
我们以三维空间中的法拉第定律为例进行推导：$\nabla \times \mathbf{E} = -\mu \frac{\partial \mathbf{H}}{\partial t}$。
1. 空间算子的离散化
考虑 $\nabla \times \mathbf{E}$ 的 $y$ 分量，其表达式为：
$$(\nabla \times \mathbf{E})_y = \frac{\partial E_x}{\partial z} - \frac{\partial E_z}{\partial x}$$
根据泰勒级数展开，对于空间导数 $\frac{\partial E_x}{\partial z}$，在点 $(i, j, k+1/2)$ 处的中心差分近似为：
$$\left. \frac{\partial E_x}{\partial z} \right|_{i,j,k+1/2} \approx \frac{E_x(i, j, k+1) - E_x(i, j, k)}{\Delta z}$$
同理，对于 $\frac{\partial E_z}{\partial x}$，在点 $(i, j, k+1/2)$ 处的中心差分为：
$$\left. \frac{\partial E_z}{\partial x} \right|_{i,j,k+1/2} \approx \frac{E_z(i+1, j, k) - E_z(i, j, k)}{\Delta x}$$
2. 时间算子的离散化
对于时间导数 $\frac{\partial \mathbf{H}}{\partial t}$，由于 $\mathbf{H}$ 在整数时刻 $n$ 定义，而 $\nabla \times \mathbf{E}$ 在半整数时刻 $n+1/2$ 定义，我们使用中心差分：
$$\left. \frac{\partial H_y}{\partial t} \right|_{n+1/2} \approx \frac{H_y^{n+1} - H_y^n}{\Delta t}$$
3. 组合得到更新方程
将上述离散化结果代入法拉第定律，得到 $H_y$ 的更新公式：
$$\frac{H_y^{n+1} - H_y^n}{\Delta t} = -\frac{1}{\mu} \left( \frac{E_x(i, j, k+1) - E_x(i, j, k)}{\Delta z} - \frac{E_z(i+1, j, k) - E_z(i, j, k)}{\Delta x} \right)$$
整理后得到 $H_y$ 的显式更新格式：
$$H_y^{n+1} = H_y^n - \frac{\Delta t}{\mu} \left( \frac{E_x^{n+1/2}(i, j, k+1) - E_x^{n+1/2}(i, j, k)}{\Delta z} - \frac{E_z^{n+1/2}(i+1, j, k) - E_z^{n+1/2}(i, j, k)}{\Delta x} \right)$$
同理，利用安培定律可以推导出 $E$ 场的更新格式：
$$E_x^{n+1/2} = E_x^{n-1/2} + \frac{\Delta t}{\epsilon} \left( \frac{H_z^{n}(i, j+1/2, k) - H_z^{n}(i, j-1/2, k)}{\Delta y} - \frac{H_y^{n}(i, j, k+1/2) - H_y^{n}(i, j, k-1/2)}{\Delta z} \right)$$
算法执行流程
在实际编程实现时，算法遵循以下迭代步骤：
- 初始化：设置初始时刻的 $E$ 和 $H$ 场分布（通常为零），设置边界条件。
- 更新磁场：利用当前时刻的 $E$ 场，根据上述推导的 $H$ 场更新方程，计算 $n+1$ 时刻的 $H$ 场。
- 更新电场：利用刚刚更新得到的 $n+1$ 时刻的 $H$ 场，计算 $n+1/2$ 时刻的 $E$ 场。
- 边界处理：应用吸收边界条件（如 PML）或反射边界条件。
- 循环：重复步骤 2-4，直到达到预设的仿真时间。
数值稳定性与 CFL 条件
由于该差分格式是显式格式，其稳定性并非无条件的。为了保证数值解不发散，时间步长 $\Delta t$ 必须满足 Courant-Friedrichs-Lewy (CFL) 条件。
在三维笛卡尔网格中，CFL 条件要求：
$$\Delta t \le \frac{1}{c \sqrt{\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} + \frac{1}{\Delta z^2}}}$$
其中 $c$ 是介质中的波速。如果 $\Delta t$ 超过这个阈值，数值误差会随时间指数级增长，导致计算崩溃。这反映了一个物理直觉：数值波速必须大于或等于物理波速，信息才能在网格间正确传递。
总结
麦克斯韦方程组的差分格式推导，本质上是将连续的微分算子通过泰勒展开转化为离散的差分算子。通过 Yee 网格的交错设计，我们不仅实现了空间和时间上的二阶精度，还巧妙地满足了麦克斯韦方程组内部的旋度耦合关系。这种方法简单、高效，是现代电磁场数值仿真软件（如 CST, HFSS 的时域模块）的理论基石。