基于RKADE-SS-FDTD 方法的无条件稳定高阶CFS-PML 算法

李建雄，魏之栋

（1.天津工业大学电子与信息工程学院，天津 300387；2.天津工业大学 天津市光电检测技术与系统重点实验室，天津 300387）

随着计算机技术与电磁基础理论的结合日益紧密，时域有限差分法（finite-difference time-domain，FDTD）在电磁学的各个领域得到了广泛的应用。然而，传统FDTD 方法的计算速度受到柯朗-弗里德里希斯-列维（Courant-Friedrichs-Lewy，CFL）条件的限制。因此，学者们提出了一些无条件稳定的FDTD 方法，如交替方向隐式FDTD（alternatingdirectionimplicitFDTD，ADI-FDTD）[1-3]。但是，当时间步长较大时，ADI-FDTD方法存在较大的离散误差和色散误差[4]。分裂步FDTD（split-step FDTD，SS-FDTD）方法[5-7]也是一种无条件稳定的算法。SS-FDTD 方法根据“Strang”分裂格式将1 个时间步分裂为3 个子步。与ADI-FDTD 方法相比，该方法具有更好的精度[8]。

等离子体作为一种色散介质，在天线、等离子体光子晶体器件和隐身技术中发挥着重要作用[9-12]。而二维场景下的等离子体算法也在很多领域中有广泛的应用，如确定反射计系统的各种参数以及应对航天器返回时的通信中断问题[13-14]等。到目前为止，学者们已经提出了几种解决色散介质的方法。在这些方法中，辅助微分方程（auxiliary differential equation，ADE）方法[15-17]避免了复杂的转换和卷积，是一种简单有效的色散介质仿真方法。近年来，Song 等[18-19]学者提出了一种龙格-库塔辅助微分方程（Runge-Kutta ADE，RKADE）的算法。该算法采用龙格-库塔技术代替线性差分技术对非磁化等离子体的微分方程进行离散。与传统的ADE 方法相比，RKADE 具有更高的精度，并且不会增加存储和计算负担。但是，Song 等仅将该算法和ADI 算法进行了结合，依然存在ADI 算法在大时间步下计算精度较低的问题。

在模拟开放或半开放的电磁辐射和散射问题时，为了用有限的网格空间模拟无限空间，必须在计算区域的截断边界处给出吸收边界条件（absorbing boundary condition，ABC）。在实际计算中，完美匹配层（perfectly matched layer，PML）是目前最常用也是吸收效果最好的一种ABC 算法。其中，复频率偏移PML（complex frequency-shifted PML，CFS-PML）具有实现简单、对倏逝波吸收效果好等优点，深受学者们的欢迎。此外，由于高阶PML（higher-order PML，HO-PML）的吸收效果优于一阶PML（first-order PML，FO-PML），近年来CFS-PML，特别是HO-CFS-PML 受到越来越多的关注以及广泛的研究[20-21]。但是目前缺少对结合RKADE算法的PML 算法研究，这限制了RKADE 算法在开放空间下的应用。

本文提出了一种截断等离子体介质的二维无条件稳定HO-RKADE-SS-CFS-PML 算法，将RKADE算法与SS-FDTD 算法相结合，消除了CFL 稳定性条件的限制；并且通过采用ADE 方法在计算域边界处实现HO-CFS-PML。通过数值算例验证该算法的有效性，并与FO-RKADE-SS-CFS-PML 算法和HORKADE-ADI-CFS-PML 算法相比较。

1 HO-RKADE-SS-CFS-PML 算法公式推导

截断非磁化等离子体介质的HO-CFS-PML 在二维TMz 波下的麦克斯韦微分方程组可以写为：

式中：Ez为电场分量；Hx、Hy为磁场分量；ε0为真空下的介电常数；μ0为真空下的磁导系数；Jz为极化电流密度分量；v 为电子碰撞频率；ωp为等离子体角频率；Sη（η=x，y）为HO-CFS-PML 的拉伸坐标变量，定义为：

为了得到Jz分量的值，式（4）采用二阶龙格库塔格式进行离散，以第1 个子时间步为例，可以得到：

为了保证局部截断误差尽可能小，各参数需要满足如下条件：

将式（12）代入式（6）中，可以得到第1 个子时间步中的电磁场分量公式为：

式中：Fzxn和Gyxn（n=1，2）均为辅助变量。下面以Fzxn为例对辅助变量的更新公式进行推导：

使用克兰克-尼科尔森（Crank-Nicolson）格式对式（16）进行离散，则其表达式为：

2 数值算例分析

为了展示HO-RKADE-SS-CFS-PML 算法的有效性，采用了一个二维空间下的数值算例。在该算例中，激励源设为调制高斯脉冲，最大频率为70 GHz，中心频率为35 GHz。每个波长下的采样点为100。计算域为100Δx×100Δy 的二维网格空间，空间步长取Δx=Δy=Δ=4.3×10-5m。计算域边界设置10 层的PML 以截断非磁化等离子体介质。传统FDTD 算法的最大时间步长为=10-13s，而无条件稳定的算法时间步长为Δt=CFLN·，CFLN 为正数。等离子体参数为ωp=2π×28.7 Grad/s，v=20 GHz。观察点设置在计算域左下角，距离两边PML 仅1 个空间单位，算例的模型图如图1 所示。

图1 数值算例草图Fig.1 Sketch of numerical example

相对反射误差定义为：

图2 各CFS-PML 算法的相对反射误差Fig.2 Relative reflection errors of different CFS-PML algorithms

由图2（a）可知，在CFLN=1 时，HO-RKADE-SSCFS-PML 算法拥有最小的相对反射误差，为-116 dB，但相比其他算法优势不够明显。由图2（b）可知，随着CFLN 的增大，RKADE-ADI-CFS-PML 算法的相对反射误差会明显增大，而HO-RKADE-SS-CFS-PML 算法则基本保持不变，其BRRE 值在不同CFLN 的时候均在-116 dB 附近。这说明HO-RKADE-SS-CFS-PML算法不仅拥有更好的电磁波吸收效果，而且不会随着CFLN 的增大而减弱，因此，该算法在CFLN 较大时会有更明显的优势。

表1 给出了各个PML 算法所需的CPU 时间、相对传统CFS-PML 算法的时间减少率以及BRRE 值。

表1 各CFS-PML 算法的CPU 时间、时间减少率以及BRRETab.1 CPU time，time reduction and BRRE of different CFS-PML algorithms

由表1 可知，RKADE-SS-CFS-PML 和RKADEADI-CFS-PML 算法都可以通过增大时间步长来减少计算时间，当CFLN=8 时，各算法消耗的时间相差不大，均能减少80%左右的时间。但是此时HO-RKADESS-CFS-PML 算法的BRRE 为-115 dB，远远小于其他的算法。由此说明，HO-RKADE-SS-CFS-PML 算法在大时间步长下拥有更好的电磁波吸收效率。

算法的稳定性可以通过波形的长时间表现[21]来证明。图3 所示为在CFLN=8 时，HO-RKADE-SSCFS-PML 算法的Ez、Hx和Hy分量。该仿真的时间扩大为上述算例的1 000 倍。由图3 可知，各分量的波形均长期收敛于零点，证明该算法具有无条件稳定性，并且PML 具有有效的吸收效果。

图3 CFLN=8 时电场和磁场的波形Fig.3 Waveforms of electric field and magnetic field when CFLN=8

3 结论

本文提出了一种无条件稳定HO-RKADE-SSCFS-PML 算法，用于截断非磁化等离子体。数值算例结果表明：

（1）该算法克服了CFL 稳定条件的限制，在CFLN=8 时可以节省79.2%的计算时间，与传统RKADECFS-PML 算法相比，提高了计算效率。

（2）在CFLN 相同时，本算法的BRRE 值小于FO-RKADE-SS-CFS-PML 算法和2 种RKADE-ADICFS-PML 算法，证明本文所提出的算法拥有更强的电磁波吸收能力。

（3）与2 种RKADE-ADI-CFS-PML 方法不同的是，HO-RKADE-SS-CFS-PML 的吸收性能不会随着CFLN 的增加而产生明显的下降。在CFLN=8 时，本算法的BRRE 值可以达到-115 dB，远远大于HORKADE-ADI-CFS-PML 算法的-88 dB，证明其在大时间步长下拥有更大的优势。