1. 稳定SBP-SAT FDTD子网格方法解析
在电磁场数值模拟领域，有限差分时域（FDTD）方法因其直观的物理意义和广泛的适用性，已成为解决复杂电磁问题的标准工具。然而，当面对包含精细几何结构或复杂材料分布的电磁问题时，传统FDTD方法面临两个主要挑战：一是阶梯误差导致的几何失真，二是全局网格细化带来的计算资源剧增。
1.1 传统FDTD方法的局限性
标准Yee网格FDTD方法在处理曲面边界时，会因阶梯近似引入几何误差。以圆形结构为例，当网格尺寸与曲率半径相当时，离散化后的结构会呈现明显的锯齿状轮廓。这种几何失真不仅影响场分布的视觉呈现，更会导致谐振频率偏移、散射特性改变等实质性误差。
全局网格细化虽能缓解阶梯误差，但代价高昂。根据CFL稳定性条件，时间步长与空间步长成正比，网格加密意味着时间步长必须相应减小。对于三维问题，计算复杂度与网格数的四次方成正比，这使得全局细化在实际工程中往往不可行。
1.2 子网格技术的优势与挑战
子网格技术通过在关键区域局部加密网格，实现了精度与效率的平衡。其核心思想可类比于摄影中的"焦点堆叠"技术——只在需要清晰呈现的区域使用高分辨率，其余区域保持基础分辨率。这种方法特别适用于以下场景：
- 天线设计中辐射单元的精细结构
- 生物电磁学中不同组织的交界面
- 超材料器件的亚波长结构
然而，子网格技术的实现面临严峻的稳定性挑战。在粗细网格交界处，电磁场分量具有不同的空间分辨率，需要精心设计的插值方案。传统方法如对称算子法、惠更斯面法等，虽然能保证短时精度，但在长时间仿真中会因误差累积导致数值发散。
2. SBP-SAT框架的数学基础
2.1 求和部分(SBP)算子原理
SBP算子的核心在于离散微分运算的能量守恒特性。在连续域中，微分运算满足分部积分公式： ∫(du/dx)v dx = uv|边界 - ∫u(dv/dx)dx
SBP算子通过在离散域精确重现这一性质，确保数值解的能量行为与物理系统一致。具体实现时，微分矩阵D与范数矩阵P需满足： D^T P + PD = B 其中B为边界矩阵，仅包含边界点贡献。这种结构将系统能量变化严格限制在物理边界，避免内部数值耗散。
2.2 同步近似项(SAT)边界处理
SAT技术通过弱施加边界条件来调控能量流动。其本质是在控制方程中引入惩罚项，驱动数值解向物理边界条件收敛。对于子网格问题，SAT项实现粗细网格场的耦合，形式为： SAT = σP⁻¹L^T (u_fine - Iu_coarse) 其中σ为惩罚参数，L为提取算子，I为插值矩阵。适定的σ选择能确保能量在接口处既不增长也不衰减。
2.3 投影SBP算子的创新设计
传统SBP-SAT方法需要添加辅助边界节点来满足SBP性质，这会增加未知量数目并限制拓扑灵活性。本文提出的投影SBP算子通过两个关键创新解决了这一问题：
- 在标准Yee网格上构造满足SBP性质的差分算子，无需修改网格结构
- 通过范数兼容插值矩阵支持任意网格比例
投影算子的数学表达为： D_proj = M⁻¹D_std 其中M为根据局部材料特性调整的加权矩阵。这种设计保持了原网格的拓扑简洁性，同时严格满足离散能量守恒。
3. 非分裂拓扑的实现细节
3.1 整体计算域分解策略
传统子网格方法需要将粗网格域分割为多个对齐块来容纳细化区域，如图1(a)所示的15块分割。即使改进的T型连接配置[图1(b)]仍需7个块。本文方法[图1(c)]直接将细化区域嵌入单一连通粗网格域，仅需4个SAT接口。
这种非分裂拓扑的优势体现在：
- 计算复杂度降低：接口数从12(T型)或8(对齐)减少到4
- 内存需求下降：消除冗余边界节点
- 编程实现简化：无需处理复杂块连接逻辑
3.2 场分量在交错网格的分布
如图1(d)所示，细化网格直接嵌入粗网格的拓扑空隙中。电场分量(Ez)位于整数网格节点，磁场分量(Hx,Hy)位于半整数节点。这种自然嵌入避免了传统方法中的场分量重定义问题。
关键实现步骤包括：
- 构建粗网格的全局差分和范数矩阵
- 定义细化区域的局部算子
- 通过范数兼容插值矩阵连接两者
3.3 范数兼容插值矩阵构造
对于含Nh个粗网格点和Nf个细网格点的西向接口，关键步骤为：
定义粗网格范数矩阵： P'_y = diag[1,1,...,1] (Nh×Nh)
细网格采用标准SBP范数： P̂_y = diag[1/2,1,...,1,1/2] (Nf×Nf)
通过转换矩阵Bc桥接范数差异： Bc = diag[1/2,1,...,1,1/2] (Nh×Nh)
最终插值矩阵： TW = Bc Tf2c T̂W = Tc2f
这种构造严格满足稳定性条件T_W^T P'_y = P̂_y T̂W，即使对于分数网格比(如2:3)也能保证长期稳定。
4. 数值验证与性能分析
4.1 长期稳定性验证
在6m×6m的PEC腔体中设置中心细化区域(2-4m)，测试两种网格比：
- 整数比1:5：粗网格Δx=5cm，细网格1cm
- 分数比2:3：粗网格5cm，细网格7.5cm
施加150MHz高斯脉冲激励，仿真10^6时间步。图3显示两种情况下电场幅度严格有界，验证了理论分析的稳定性。
4.2 数值反射系数测量
通过波导模型量化接口反射。设置3.6m×0.27m计算域，中心嵌入精细区域，测试网格比2:5至1:10。图5显示所有情况反射系数低于-60dB，1:10比2:5仅高约2dB，证明接口处理的优越性。
4.3 双C-SRR阵列仿真比较
模拟50mm腔体中两个2×2开口环谐振器阵列，比较不同方法在1:10网格比的性能：
| 方法 | 网格点数 | 相对误差 | 加速比 |
|---|---|---|---|
| 全局细网格 | 1,002,001 | - | 1× |
| 本文方法 | 191,403 | 0.35% | 5.29× |
| 文献[31]方法 | 197,725 | 0.34% | 4.57× |
| 文献[19]方法 | 191,403 | 0.35% | 5.14× |
本文方法在保持精度的同时，因减少接口数获得更高效率。图9的空间误差分布显示，相比传统方法，本文方案在谐振器边缘区域的误差显著降低。
4.4 人头模型电磁传播仿真
1m×1m腔体中放置生物电磁模型，比较五种离散策略：
| 方法 | 网格点数 | 计算时间(s) | 相对误差 |
|---|---|---|---|
| 全局细网格 | 1,002,001 | 512.77 | - |
| 本文1:2 | 18,001 | 9.85 | 5.37% |
| 对齐块1:2 | 18,209 | 11.95 | 5.83% |
| T型连接1:2 | 18,105 | 10.90 | 6.13% |
非分裂拓扑在保持精度的同时，因减少冗余接口获得52倍加速。图12展示1:10网格比能准确分辨不同生物组织的界面。
5. 工程应用中的实施建议
5.1 网格划分策略
- 几何特征分析：识别需要高分辨率的区域（如曲面、薄层结构）
- 网格比例选择：建议整数比(1:2,1:3)优先，分数比需验证插值矩阵
- 过渡区宽度：至少3个粗网格单元，避免场突变
5.2 参数调优经验
- 惩罚参数：初始值设为-0.5，可微调±10%优化收敛
- 时间步长：取0.95×细网格CFL限制，兼顾稳定与效率
- 材料界面处理：在ε/μ突变处设置网格对齐节点
5.3 常见问题排查
- 后期发散：检查范数兼容条件是否满足，重新生成插值矩阵
- 异常谐振：验证几何模型与实际结构的一致性
- 性能下降：检查稀疏矩阵存储格式(建议CSR)
关键提示：在生物电磁仿真中，组织界面处的网格线应尽量与解剖结构对齐，可减少高达30%的误差。
6. 方法优势与适用范围
相比传统子网格技术，本方法具有三大创新点：
- 拓扑简化：单连通粗网格域直接嵌入细化区域，避免复杂分块
- 稳定性保证：严格的离散能量分析，适合长时间瞬态仿真
- 计算高效：4接口设计减少30%以上内存消耗
典型应用场景包括：
- 5G天线阵列的多尺度仿真
- 植入式医疗设备的电磁安全评估
- 超表面器件的宽带响应分析
未来工作将扩展至三维非共形网格和色散媒质，进一步提升工程适用性。
