1. 引言

随着微电子技术的发展，微电子设备已在各种军用设备和系统中获得广泛的应用。瞬息万变的信息化战场时刻离不开微电子设备，而微电子设备对电磁脉冲又极为敏感。高功率电磁脉冲对信息设备与系统构成严重威胁。电磁波在坑道中的传播具有特定的规律，当波的频率高于坑道截止频率时才可以传播，坑道截止频率一般在几十兆赫兹左右[1] 。对于非核电磁脉冲，研究表明HPM的载波频率一般在500 MHz以上，UWB的覆盖频段一般在100 MHz以上，均高于坑道的截止频率。因此，在坑道内部传播非常容易，从而对内部的电子系统及设备构成了威胁 [2] 。因此，研究高功率电磁脉冲电磁效应及工程防护技术具有重要意义。

自1966年K. S. Yee首次提出时域有限差分(Finite Difference Time Domain，FDTD)方法，经过近四十多年的发展，已经成为一种非常成熟的电磁学计算方法，具备了强大的解决多种复杂问题的能力。作为时域方法，FDTD法只需一次计算即可获得宽频带的丰富信息，可以方便的得到电磁场传播的过程，这种宽频带信息与可视化结果使其得到广泛的应用。

本文针对电磁脉冲在坑道中传播的问题，提出了一种基于FDTD方法的解决模型。仿真中采用卷积完全匹配层截断计算区域；对坑道的拱顶，采用共形的方案进行模拟；激励源采用超宽带脉冲；为克服仿真资源的占用对计算机硬件的要求，采用了并行的方案。最后我们建立的模型仿真了电磁脉冲在坑道中的传播，得到了一些规律，验证了本模型的有效性。

2. 数值方法

把带时间变量的Maxwell旋度方程转化为差分方程，差分格式中每个网格点上的电场(或磁场)分量仅与它相邻的磁场(或电场)分量以及上一时间步该点的场值有关。在每个时间步上分别计算网格空间各点的电场和磁场分量，随着时间步的推进，就能模拟出电磁脉冲与机箱的相互作用过程 [3] 。

电场以分量为例，其差分方程为

(1)

其中，，。为真空中的介电常数，为相对介电常数，为电导率。

磁场以分量为例，其差分方程为

(2)

其中，，。为真空中的磁导率，为相对磁导率。其余四个电磁场分量的差分方程可类似得到。

3. 吸收边界条件

除了一些全封闭导体的问题，一般用FDTD法求解电磁场问题时的计算域总是开放的。所取的问题域越大，要求存储量越大。由于计算机容量的限制，计算只能在有限区域进行。为了能模拟开区域电磁过程，在计算区域的截断边界处必须设置吸收边界条件，使电磁波在边界处无明显的反射以免影响计算域内场的分布。

这里我们采用一种基于卷积的吸收边界(Convolution Perfectly Matched Layers, CPML) [4] 。CPML的FDTD实现形式完全独立于所截断媒质类型，可以不加修改的应用于非均匀、有耗、各向异性、色散和非线性等媒质。在CPML中，每个场分量只需两个辅助变量，可显著节省内存。该边界条件不仅可以吸收传播模，还可以吸收低频凋落模，能够克服PML长时间计算带来的晚时反射，因此CPML可贴近散射体设置，减小计算区域。

在CPML中，电场分量迭代公式为

(3)

其中，

(4)

(5)

(6)

(7)

类似地，可以得到其他电磁场分量的迭代公式。

4. 计算模型及曲面共形处理

坑道工程通常沿轴向延伸，如图1所示，区别主要是截面的形状，此处以具有代表性的平底拱顶形截面为计算模型，如图2(a)所示。坑道底边长及侧壁高均为4.0 m，顶部拱形用半径为2.0 m的半圆近似，岩土介质层电磁参数取值：电导率，相对介电常数，相对磁导率。坑道外侧岩土介质及两端的自由空间用完全匹配层(CPML)截断，对于底边及侧壁的电磁参数取平均值，对于坑道顶部曲面部分采用共形技术处理，如图2(b)所示。

5. 激励源及其设置

波导系统中常采用强迫激励的方法引入激励源，由于坑道周围岩土为非理想导体，且截面为非典型波导横截面(如圆形、矩形)，无法获取其闭式解的场强分布表达式，可以将上述截面自由空间区域近似为矩形，采用近似的模式分布来激励。在此近似模式强迫激励下，随着波的传播将不断修正，经过一段距离的传播，坑道中的场强分布将接近其在此激励模式下固有的、满足坑道边界条件的分布模式。考虑高次模在坑道中衰减较快，这里主要以矩形波导中的两个典型模式TE10、TM11模为例，设长边大小为a，短边为b，下面作简要介绍。

1) TM11模有五个场分量，纵向磁场分量为零，场方程为：

(8)

2) TE10模有三个场分量，纵向的电场分量，横向磁场分量、横向电场分量均为零，场方程为：

(9)

设坑道纵向沿z轴方向，激励源位于处，截面场分布函数为如式(8)和(9)所示，激励源时间函数为。于是，强迫激励源可设置如下：

(10)

式中下标“tan”表示横向场分量。为激励源的横向场分量，计算过程中激励源非核电磁脉冲的时间函数采用超宽带电磁脉冲(UWB) [5] - [7] 。

UWB入射波电场近似表示为：

(11)

式中，为电场峰值，为脉冲宽度，为时间常数，时域波形如图3所示。为了计算结果具有可比性，两种脉冲电场峰值取。

6. 并行技术在坑道仿真中的运用

时域有限差分法并行的基本思路是，将总区域分割成若干个子区域，每个子区域对应一个进程，通过各个进程间的数据交换来完成整个区域的计算[8] 。对于坑道这种一维方向计算区域较大的问题，在区域分割时采用沿z轴方向一维线性分割，如图4示。

将FDTD区域划分成多个子区域后，连接边界、吸收边界、散射体都被划分到某些子区域中，由于这些边界需要特别处理，因此增加了编程的复杂性。编程实现的办法是将各种边界的位置传输给所有的子域，各个子域根据自身的网格范围来判断本子域是否包含某边界。而各个进程分界面处的场量则必须通过相互通信来传递，其计算流程与场量通信如图5所示。

进程分界面处的电场、在相邻两个进程中同时进行计算，在左进程中计算时需要用到右进程中的磁场、；在右进程中计算时也要用到左进程中的磁场、。故在计算电场之前需要先在两个相邻子域之间传输边界的切向磁场。通信的原理可以按照图5来解释：对每个时间步来说，编号为N的进程，需要从进程N − 1接收、的场值，同时给进程N + 1发送、的场值；还要从进程N + 1接收、的场值，同时给进程N − 1发送、的场值。

由于采用的是一维分割，通信量相对较少，故采用阻塞通信方式。编程中采用MPI_SENDRECV命令，使程序简化，更重要的是在通信过程中不需要考虑哪个进程先发送还是先接收，从而可以避免消息传递过程中的死锁。通信顺序与程序如图6示。

在编程中，采用主从模式(Master-slave)：有一个主进程，其它为从进程。主进程负责整个并行程序的数据控制，向从进程分配任务并输入输出。从进程负责对数据的处理和计算，主进程也可参与对数据的处理和计算。

当激励源从左侧处引入后，由于采用的是强迫激励的方法，将向两侧传播，向左传播的波将被左侧的匹配层吸收，向右传播的波沿坑道继续向前传播，直到右匹配层完全吸收，整个程序计算结束。

7. 计算结果及分析

计算中截面场强分布函数采用TE10模，其截面电场分量分布如式(9)所示，式中，时间函数为UWB，如式 (11)所示，脉冲宽度，，其他参数选取不变，根据稳定性条件及数值色散对空间步长与时间步长的要求，，，为脉冲在坑道中的传播速度，考虑坑道截面大小及激励源的高频特性，以后均以光速近似代替。

在上述激励源及参数选取下，模拟的坑道长度为40 m，图7(a)所示为某一时刻记录的分量空间分布，图7(b)为相对源位置不同距离观察点的输出波形。从脉冲波空间分布分析，可知传播过程中传播能量集中在一个较窄的区域内；从观察点输出波形发现随传播的增加，在时域其波形逐渐变宽，对应的频域变窄。因为UWB是一个宽带信号，频率丰富，而坑道中不同频率成分的波随距离的传播衰减率不同，因此在传播过程中UWB产生了严重的色散现象，从而导致脉冲波时域变宽。在频域随着脉冲波传播，低频分量较高频衰减快而逐渐消失，能量主要靠截止频率以上的高频部分传播，另一方面脉冲波在坑道

中传播，其能量将不断向周围半有耗介质扩散，对应的幅值不断减小。这些现象与实际的物理过程是相一致的。

本文针对电磁脉冲在坑道中传播的问题，提出了一种基于FDTD方法的解决模型。仿真中采用卷积完全匹配层截断计算区域；对坑道的拱顶，采用共形的方案进行模拟；激励源采用超宽带脉冲；为克服仿真资源的占用对计算机硬件的要求，采用了并行的方案。最后我们建立的模型仿真了电磁脉冲在坑道中的传播，得到了一些规律，验证了本模型的有效性。

基金项目

受中国人民解放军理工大学电磁环境效应与光电工程重点实验室资助，资助号：FD201506。

参考文献