这篇笔记在做什么

这篇笔记从一份 Playground 导出会话整理而来。我保留了原始推导里的细节，重点把最简单的三维 PEC 立方体散射问题整理成一个可连续阅读的入门版本。

文章会按下面这条线展开：Maxwell 方程、散射场分解、ABC 截断、弱形式、棱边元、六面体单元、局部到全局编号、矩阵组装。

三维电磁场散射问题的有限元方法——以立方体为例

一、物理问题描述

考虑一个边长为 $a$ 的理想导体(PEC)立方体，置于自由空间中，被一个平面波照射。我们要求解散射场。

        入射平面波 E^inc
             ↓  ↓  ↓
         ___________
        /           /|
       /           / |
      /___________/  |    PEC 立方体
      |           |  |
      |           |  /
      |           | /
      |___________|/

    求：散射场 E^sca

二、理论推导

2.1 从 Maxwell 方程出发

时谐场（$e^{j\omega t}$ 约定）下的 Maxwell 方程：

$$ \nabla \times \mathbf{E} = -j\omega\mu\mathbf{H} $$$$ \nabla \times \mathbf{H} = j\omega\varepsilon\mathbf{E} + \mathbf{J} $$

消去 $\mathbf{H}$，得到矢量波动方程：

$$ \nabla \times \left(\frac{1}{\mu_r} \nabla \times \mathbf{E}\right) - k_0^2 \varepsilon_r \mathbf{E} = -jk_0 Z_0 \mathbf{J} $$

其中：

$k_0 = \omega\sqrt{\mu_0\varepsilon_0}$：自由空间波数
$Z_0 = \sqrt{\mu_0/\varepsilon_0}$：自由空间波阻抗

2.2 散射场分解

总场 = 入射场 + 散射场：

$$ \mathbf{E}^{tot} = \mathbf{E}^{inc} + \mathbf{E}^{sca} $$

散射场满足：

$$ \nabla \times \left(\frac{1}{\mu_r} \nabla \times \mathbf{E}^{sca}\right) - k_0^2 \varepsilon_r \mathbf{E}^{sca} = 0 \quad \text{(散射体外)} $$

2.3 计算域截断

无限空间无法直接计算，需要用人工边界截断：

    +------ ABC/PML 边界 ------+
    |                           |
    |    +---+                  |
    |    |PEC|  散射体          |
    |    +---+                  |
    |                           |
    |    Ω: 有限元求解区域       |
    +---------------------------+

常用方法：

一阶吸收边界条件 (ABC)：$\hat{n} \times (\nabla \times \mathbf{E}) + jk_0 \hat{n} \times (\hat{n} \times \mathbf{E}) = 0$
完美匹配层 (PML)（更精确）

2.4 弱形式（变分形式）

这是有限元的核心步骤。

步骤： 用测试函数 $\mathbf{T}$ 对波动方程做内积，并利用矢量恒等式分部积分。

矢量恒等式：

$$ \mathbf{T} \cdot (\nabla \times \nabla \times \mathbf{E}) = (\nabla \times \mathbf{T}) \cdot (\nabla \times \mathbf{E}) + \nabla \cdot (\mathbf{T} \times \nabla \times \mathbf{E}) $$

对体积积分，并用散度定理：

$$ \boxed{ \int_\Omega \left[ \frac{1}{\mu_r}(\nabla \times \mathbf{T}) \cdot (\nabla \times \mathbf{E}) - k_0^2 \varepsilon_r \mathbf{T} \cdot \mathbf{E} \right] dV - \oint_S \mathbf{T} \cdot \left(\hat{n} \times \frac{1}{\mu_r}\nabla \times \mathbf{E}\right) dS = 0 } $$

其中面积分项通过边界条件处理：

PEC 面：$\hat{n} \times \mathbf{E} = 0$（本质边界条件，直接强加）
ABC 面：代入 ABC 公式，产生附加矩阵项

代入一阶 ABC 后，弱形式变为：

$$ \int_\Omega \frac{1}{\mu_r}(\nabla \times \mathbf{T}) \cdot (\nabla \times \mathbf{E}) , dV

k_0^2 \int_\Omega \varepsilon_r \mathbf{T} \cdot \mathbf{E} , dV
jk_0 \oint_{S_{ABC}} (\hat{n} \times \mathbf{T}) \cdot (\hat{n} \times \mathbf{E}) , dS = \text{右端项} $$

三、为什么用棱边元（Edge Element）而不是节点元？

这是三维电磁 FEM 最关键的概念：

问题	节点元	棱边元（Nédélec元）
自由度位置	节点上	棱边上
物理量	标量	场沿棱边的切向分量
伪解(spurious modes)	❌ 会产生	✅ 不会
切向连续性	难保证	✅ 自动保证
法向不连续性	难处理	✅ 自动允许

结论：三维电磁问题必须用棱边元！

四、立方体网格的离散与编号

4.1 最简单的例子：一个六面体单元

先看一个六面体（就是一个小立方体）：

        7 _____________ 8
         /|            /|
        / |           / |
     5 /__|__________/6 |
       |  |          |  |
       |  3__________|__| 4
       |  /          |  /
       | /           | /
       |/____________|/
      1               2

8 个节点编号：1-8

12 条棱边编号与方向：

棱边编号    起点 → 终点     方向
─────────────────────────────────
  e1        1 → 2          x
  e2        3 → 4          x
  e3        5 → 6          x
  e4        7 → 8          x
  e5        1 → 3          y
  e6        2 → 4          y
  e7        5 → 7          y
  e8        6 → 8          y
  e9        1 → 5          z
  e10       2 → 6          z
  e11       3 → 7          z
  e12       4 → 8          z

每个六面体单元有 12 个自由度（12条棱边，每条1个DOF）。

4.2 棱边基函数

对于线性六面体棱边元，第 $i$ 条棱边的基函数 $\mathbf{N}_i$ 的关键性质：

$$ \int_{e_j} \mathbf{N}_i \cdot d\mathbf{l} = \delta_{ij} $$

即：基函数沿自身棱边的线积分为1，沿其他棱边为0。

在参考单元 $(\xi, \eta, \zeta) \in [0,1]^3$ 中，12个基函数为：

x方向的4条棱（e1-e4）：

$$ \mathbf{N}_1 = (1-\eta)(1-\zeta) \hat{\xi} $$

$$ \mathbf{N}_2 = \eta(1-\zeta) \hat{\xi} $$

$$ \mathbf{N}_3 = (1-\eta)\zeta \hat{\xi} $$

$$ \mathbf{N}_4 = \eta \cdot \zeta \hat{\xi} $$

y方向的4条棱（e5-e8）：

$$ \mathbf{N}_5 = (1-\xi)(1-\zeta) \hat{\eta} $$

$$ \mathbf{N}_6 = \xi(1-\zeta) \hat{\eta} $$

$$ \mathbf{N}_7 = (1-\xi)\zeta \hat{\eta} $$

$$ \mathbf{N}_8 = \xi \cdot \zeta \hat{\eta} $$

z方向的4条棱（e9-e12）：

$$ \mathbf{N}_9 = (1-\xi)(1-\eta) \hat{\zeta} $$

$$ \mathbf{N}_{10} = \xi(1-\eta) \hat{\zeta} $$

$$ \mathbf{N}_{11} = (1-\xi)\eta \hat{\zeta} $$

$$ \mathbf{N}_{12} = \xi \cdot \eta \hat{\zeta} $$

4.3 场的展开

单元内电场用棱边基函数展开：

$$ \mathbf{E} = \sum_{i=1}^{12} E_i \mathbf{N}_i $$

其中 $E_i$ 就是电场沿第 $i$ 条棱边的线积分值（自由度）。

五、单元矩阵的建立

将展开式代入弱形式，得到单元矩阵方程：

$$ \left([S^e] - k_0^2 [T^e]\right) \{E^e\} = \{b^e\} $$

其中 $12 \times 12$ 的单元矩阵为：

$$ S^e_{ij} = \int_{\Omega^e} \frac{1}{\mu_r} (\nabla \times \mathbf{N}_i) \cdot (\nabla \times \mathbf{N}_j) \, dV $$$$ T^e_{ij} = \int_{\Omega^e} \varepsilon_r \, \mathbf{N}_i \cdot \mathbf{N}_j \, dV $$

具体计算 $\nabla \times \mathbf{N}_i$：

以 $\mathbf{N}_1 = (1-\eta)(1-\zeta) \hat{\xi}$ 为例：

$$ \nabla \times \mathbf{N}_1 = \begin{vmatrix} \hat{\xi} & \hat{\eta} & \hat{\zeta} \\ \frac{\partial}{\partial \xi} & \frac{\partial}{\partial \eta} & \frac{\partial}{\partial \zeta} \\ (1-\eta)(1-\zeta) & 0 & 0 \end{vmatrix} $$$$ = (1-\eta) \hat{\eta} + (-(1-\zeta)) \hat{\zeta} $$$$ \nabla \times \mathbf{N}_1 = (1-\eta)\hat{\eta} - (1-\zeta)\hat{\zeta} $$

（注意：实际计算需要通过 Jacobian 变换到物理坐标）

六、多个单元的网格与全局编号

6.1 将立方体剖分为 $2 \times 2 \times 2 = 8$ 个单元

俯视图（z = 常数的一层）:

       7 ----14---- 8 ----15---- 9
       |            |            |
      25    e3     26    e4     27
       |            |            |
       4 ----11---- 5 ----12---- 6
       |            |            |
      22    e1     23    e2     24
       |            |            |
       1 ---- 8---- 2 ---- 9---- 3

  (底层 z=0，节点编号1-9，棱边编号示意)

完整统计（$2\times2\times2$ 网格）：

项目	数量	计算公式
节点	$3 \times 3 \times 3 = 27$	$(n+1)^3$
x方向棱边	$2 \times 3 \times 3 = 18$	$n(n+1)(n+1)$
y方向棱边	$3 \times 2 \times 3 = 18$	$(n+1)n(n+1)$
z方向棱边	$3 \times 3 \times 2 = 18$	$(n+1)(n+1)n$
总棱边数	54	$3n(n+1)^2$
单元数	8	$n^3$

全局自由度数 = 总棱边数 = 54

6.2 局部-全局编号映射（关键！）

这是组装的核心。以两个相邻单元为例：

 单元1          单元2
  4---e2---3  |  4---e2---3
  |        |  |  |        |
  e5  e1  e6  |  e5      e6
  |        |  |  |        |
  1---e1---2  |  1---e1---2

  共享棱边!

映射表示例（单元1和单元2共享一条棱边）：

                局部棱边号    单元1全局编号    单元2全局编号
                ─────────────────────────────────────────
                    e1            1               2
                    e2            4               5
                    e3            8               9
                    e4           11              12
                    e5           22              23  ← 共享!
                    e6           23              24
                    ...          ...             ...

单元2的局部棱边e5 和单元1的局部棱边e6 映射到同一个全局编号23！

6.3 组装过程

全局矩阵 $[A]$ 大小为 $54 \times 54$：

$$ A_{IJ} = \sum_{\text{所有单元 } e} \sum_{i,j=1}^{12} \left(S^e_{ij} - k_0^2 T^e_{ij}\right) \cdot \delta_{I, \text{map}(e,i)} \cdot \delta_{J, \text{map}(e,j)} $$

用伪代码表示：

# 初始化全局矩阵
A = zeros(54, 54)  # 全局棱边数 x 全局棱边数
b = zeros(54)

for elem in range(8):          # 遍历每个单元
    # 计算 12x12 单元矩阵
    Ke = compute_element_matrix(elem)  # S^e - k0^2 * T^e
    fe = compute_element_rhs(elem)
    
    for i in range(12):        # 局部棱边 i
        I = edge_map[elem][i]  # 局部 → 全局映射
        for j in range(12):    # 局部棱边 j
            J = edge_map[elem][j]
            A[I][J] += Ke[i][j]   # 组装!
        b[I] += fe[i]

# 施加 PEC 边界条件
for edge in PEC_edges:
    A[edge, :] = 0
    A[edge, edge] = 1
    b[edge] = 0

# 求解
E = solve(A, b)

七、边界条件的处理

7.1 PEC 边界条件

PEC 表面：$\hat{n} \times \mathbf{E} = 0$

→ PEC面上所有切向棱边的自由度 = 0

    PEC面上的棱边
    ____________
    |          |
    | 这些棱边 |  ← 这些DOF直接设为0
    | 的DOF=0  |
    |__________|

7.2 ABC 边界条件

一阶ABC产生附加面积分矩阵（$B$矩阵）：

$$ B_{ij} = -jk_0 \oint_{S_{ABC}} (\hat{n} \times \mathbf{N}_i) \cdot (\hat{n} \times \mathbf{N}_j) \, dS $$

在ABC面的每个面单元上（$2\times2$ 格子有4条边，4个DOF）：

$$ [A]_{global} = \sum_e ([S^e] - k_0^2 [T^e]) + \sum_{f \in ABC} [B^f] $$

八、完整流程总结

┌─────────────────────────────────┐
│  1. 定义几何：PEC立方体 + 外域  │
└──────────────┬──────────────────┘
               ▼
┌─────────────────────────────────┐
│  2. 网格剖分：六面体/四面体网格  │
│     记录节点坐标、单元连接关系    │
└──────────────┬──────────────────┘
               ▼
┌─────────────────────────────────┐
│  3. 建立棱边表：                 │
│     全局棱边编号                 │
│     局部→全局映射表              │
│     标记边界棱边(PEC/ABC)        │
└──────────────┬──────────────────┘
               ▼
┌─────────────────────────────────┐
│  4. 逐单元计算：                 │
│     [S^e] (旋度-旋度矩阵 12×12) │
│     [T^e] (质量矩阵 12×12)      │
│     用高斯积分计算               │
└──────────────┬──────────────────┘
               ▼
┌─────────────────────────────────┐
│  5. 组装全局矩阵 [A], {b}       │
│     A = Σ(S^e - k₀²T^e) + B    │
└──────────────┬──────────────────┘
               ▼
┌─────────────────────────────────┐
│  6. 施加边界条件：               │
│     PEC: 删除/置零对应行列       │
│     ABC: 加上B矩阵              │
└──────────────┬──────────────────┘
               ▼
┌─────────────────────────────────┐
│  7. 求解线性方程组               │
│     [A]{E} = {b}                │
│     (稀疏、复数、非对称)         │
└──────────────┬──────────────────┘
               ▼
┌─────────────────────────────────┐
│  8. 后处理：                     │
│     近场分布、远场RCS等          │
└─────────────────────────────────┘

九、数值示例的规模感受

参数	值
立方体边长	$1\lambda$（一个波长）
每波长单元数	10
外域厚度	$0.5\lambda$
总网格	$\approx 20 \times 20 \times 20 = 8000$ 单元
总棱边（DOF）	$\approx 25000$
矩阵大小	$25000 \times 25000$（稀疏）

这就是用最简单的立方体几何理解三维电磁散射FEM的完整框架。核心要点是：

矢量波动方程 → 弱形式
棱边元（不是节点元）→ 避免伪解
每条棱边一个自由度 → 场的切向分量
局部-全局映射 → 矩阵组装
ABC/PML → 截断无限域

这篇笔记在做什么#

三维电磁场散射问题的有限元方法——以立方体为例#

一、物理问题描述#

二、理论推导#

2.1 从 Maxwell 方程出发#

2.2 散射场分解#

2.3 计算域截断#

2.4 弱形式（变分形式）#

三、为什么用棱边元（Edge Element）而不是节点元？#

四、立方体网格的离散与编号#

4.1 最简单的例子：一个六面体单元#

4.2 棱边基函数#

4.3 场的展开#

五、单元矩阵的建立#

六、多个单元的网格与全局编号#

6.1 将立方体剖分为 $2 \times 2 \times 2 = 8$ 个单元#

6.2 局部-全局编号映射（关键！）#

6.3 组装过程#

七、边界条件的处理#

7.1 PEC 边界条件#

7.2 ABC 边界条件#

八、完整流程总结#

九、数值示例的规模感受#