平面波在时域有限差分法中的引入方法研究

柴焱杰 孙继银 孙东阳 胡 寅

(1.第二炮兵工程学院，陕西 西安 710025； 2.西北核技术研究所，陕西 西安 710024； 3.第二炮兵指挥学院，湖北 武汉 430012)

1.引 言

时域有限差分法(FDTD，The finite difference time domain method)是研究电磁问题的一种迅速发展的仿真计算方法，能够在时域直接计算得到宽带结果，适合于分析复杂电磁系统。引入平面波激励源是各类电磁仿真研究的重要内容，是进行场效应计算的基础前提[1-5]。FDTD中，仿真平面波激励源时使用总场-散射场(TF-SF)连接边界条件[6]，如图1所示。

图1 FDTD总场-散射场计算模型

根据等效原理，在连接边界上设置平面波的等效面电磁流，并设平面外的场为零，就可将入射波只引入到总场区。连接边界上的等效电磁流为

(1)

式中：en为面的外法向；Ei、Hi分别为入射电场(V/m)和入射磁场(A/m)。通常使用球坐标系描述平面波激励源的传播方向和极化状态，使用三维直角坐标系进行FDTD迭代运算，因此，需要解决平面波激励源的投影问题，如图2所示[7]。

图2 HEMP平面波使用的球坐标系

图2中α是极化角，坐标系中的原点O位于连接边界贴近三维直角坐标系中FDTD计算空间原点的那个角点。由于平面波传播方向k由角度(θ，φ)标定，而一般情况下FDTD仿真计算模型位于第一卦限中，因此，当入射角超越这个范围，将会引起平面波不能“照射”(引入)到计算区域内的情况。为方便研究，使计算模型更为灵活、符合实际情况，研究两种任意角度引入平面波激励源的方法， 并对其仿真效果进行评价。

2.理论分析

2.1 基于角点延迟引入任意角度平面波[6-9]

为使任意角度平面波能够入射到仿真空间，将三维FDTD连接边界上的八个角点编号O1～O8，如图3所示。角点延迟的思想是，当平面波以任意角度(θ，φ)入射该区域时，无论该角度是多少，总可以认为该平面波是以其中某一角点出发引入进来的。

图3 连接边界上的角点编号

由于平面波传播方向k在直角坐标系统中的投影为

er=(kx,ky,kz)

=(sinθcosφ,sinθsinφ,cosθ)

(2)

因此，三个分量kx、ky、kz反映了从原点O出发的k在哪个卦限的信息，从而可选择某一角点作为基准点Os，用以计算一维平面波相对一维源点的延迟距离。设点变量Os(INCX,INCY,INCZ)是平面波的出发点，依据kx、ky、kz的符号可按表1的选择条件确定基准点Os.

表1 基准点Os的确定

图3中同时示出了以O2点作为基准点Os的情况。当入射平面波的方向矢量er满足条件：kx<0,ky≥0,kz≥0，此时平面波将入射到第二卦限。选择O2点作为基准点Os，即平移原入射方向，就能够保证平面波进入到FDTD总场空间中。

(3)

根据延迟距离r′计算出某点的场量后，即可转换为连接边界上该点的切向场分量：

(4)

(5)

由式(4)、(5)计算出的电、磁切向场分量可实现任意角度平面波在三维FDTD中的传播。根据延迟距离r′计算场量的方法主要有两种——一维平面波推进法和解析法。需要注意的是，无论是一维平面波推进法还是解析法，由于电场E和磁场H在空间上相隔半个空间步，而在八个基准点Os上均没有场量，因此需要仔细推算连接边界上的场点与基准点之间的实际距离r′，这是实现正确投影和坐标系统转换的前提。

2.2 一维平面波推进法

一维平面波推进法将入射平面波用一维FDTD随时间逐步推进的方式提供连接边界上的切向场分量。一维FDTD平面波由式(6)表述。

(6)

一维FDTD场点的网格其时间步长、空间步长与三维FDTD计算空间相同。一维FDTD区域的边界应留有若干网格以便加入吸收边界条件。一维FDTD总长度与总场空间(连接边界内的区域)有关，应满足条件

(7)

式中：ZNX、ZNY、ZNZ是三维FDTD中总场空间三个方向上的长度(网格数)；UNZ是一维FDTD吸收边界的厚度(网格数)。

r′标明了一维FDTD平面波上与三维FDTD连接边界中的P点所对应的坐标位置。r′有时可能不是整数，设r′=(p+w)△x，0

(8)

2.3 解析法

(9)

3.实验结果与分析

为验证算法的正确性，选择高空核爆电磁脉冲(HEMP，High-altitude electromagnetic pulse)作为研究对象。HEMP是指发生在30 km以上的高空核爆炸产生的电磁脉冲，是电磁兼容等领域研究的重要内容。当前应用比较广泛的描述HEMP波形的标准有1976年出版物标准、Bell实验室标准和国际电工委员会(IEC)制定的HEMP标准等[10]。这些对HEMP的描述中，将HEMP辐射波形拟合为双指数函数表达式，其中IEC制定的HEMP标准波形为

E(t)=1.3×5×104×(e-4×107×t-e-6×108×t)

(10)

IEC定义的HEMP时域波形和归一化频谱如图4所示。从图4可见，HEMP峰值达到50,000 V/m；上升时间和衰落时间分别是2.5 ns和55 ns；波形能流分布的主要频段范围(能流比例在98%)是100 kHz～100 MHz[11]。因此，HEMP具有高峰值场强、快上升沿、宽频带等特点，会对电子设备构成严重威胁。

图4 IEC标准定义的HEMP波形及其归一化频谱

当使用一维平面波推进法引入HEMP平面波时，一维FDTD源点的HEMP脉冲波形离散化为

e-6×108×(n+1)×Δt)

(11)

式中，k0是源点的位置，某位置的一维平面波场量由式(6)递推得出。为使用解析法引入HEMP平面波，延迟距离为r′的某点其HEMP的电场值离散化为

(12)

3.1 计算实例

为满足Courant稳定性条件以及控制数值色散，可设置网格大小Δs=c/(20fm)=0.15 m，时间步长Δt=Δs/(2c)=0.25 ns.为便于观察，将总网格数目设置为90×90×90，其中总场区的网格数目为60×60×60。图5(见413页)是θ=90°，α=180°，φ=0°时使用一维平面波推进法计算场值Ez随时间步推进的两个网格图片段。

为比较两种引入HEMP平面波方法的效果，设置θ=90°，α=180°，φ=0°、135°，此时电场E只有z分量。取z=NZ/2(位于计算空间一半高度的平面)为观察面。两种方法计算HEMP平面波传播的等高线图如图6(见413页)所示。最后，设置几种网格尺寸，对几种情况下入射电场在散射场区泄露的最大值进行统计和比较，结果如表2、表3所示。

表3 φ=135°时，入射电场在散射场区泄露的最大值 /(V/m)

3.2 仿真结果分析

由图5～6，表2～3可以看出：

1)φ=0°时，基准点为O1，HEMP平面波的k应沿x轴方向；φ=135°时，基准点为O2，HEMP平面波的k应沿-x轴与y轴中间45°方向。从仿真结果可以看出，基于两种方法的场波形均沿正确的方向传播。

2) 一维平面波推进法中，由于一维FDTD存在的数值色散能够很大程度上抵消三维FDTD存在的固有数值色散，因而总场区域内的结果比较准确，总场区域外则较为纯净。解析法虽然计算出了连接边界上各点的准确值，却因不能抵消这种数值色散，而在波前附近出现了幅值较大的电场泄露；总场区域内也出现了不均匀的现象：等高线图中出现了较多的噪点，总场区中亦出现较多的曲线，这有悖于平面波的均匀特性。

3) 在一维平面波推进法的计算结果中可以发现，当φ=0°时，一维FDTD平面波与三维FDTD空间的网格大小吻合，没有实际的插值现象，入射电场在散射场区的泄露非常小，接近于计算机的基本计算误差，结果很理想；当φ=135°时，由于在r′为非整数的位置使用了插值的方法(如式(8))，此时亦出现了少许波的溢出现象，然而，入射电场在散射场区的泄露会随着网格尺寸的减小而减小，即减小网格尺寸能够明显降低HEMP的泄露。从泄露的最大值上看，一维平面波推进法要明显优于解析法。

4.结 论

在论述基于角点延迟思想的基础上，从原理上推导了将HEMP平面波引入到三维FDTD中的两种方法——一维平面波推进法和解析法，解决了任意入射角度HEMP平面波引入三维FDTD计算空间的问题。对HEMP平面波在自由空间的传播过程进行了仿真验证和比较。结果表明:两种方法都能以任意角度将平面波引入到总场区域，计算代价小；一维平面波推进法很大程度上抵消了FDTD法的数值色散误差，溢出波较少，而解析法出现了较多的溢出波。工程应用中如果使用一维平面波推进法，可通过适当减小网格尺寸进一步改善平面波的仿真效果。

[1] 刘顺坤, 聂 鑫, 陈向跃. 电磁脉冲对电缆耦合问题的理论研究[J]. 电波科学学报, 2010, 25(2): 348-352.

LIU Shunkun, NIE Xin, CHEN Xiangyue. Numerical study on cable coupling effects excited by electromagnetic pulse [J]. Chinese Journal of Radio Science, 2010, 25(2): 348-352. (in Chinese)

[2] TESCHE F M, IANOZ M V, KARLSSON T. EMC Analysis Methods and Computational Models[M]. New York: John Wiley & Sons, Inc., 1997.

[3] SULLIVAN D M. Electromagnetic Simulation Using the FDTD Method[M]. Idaho: IEEE Press, 2000.

[4] 高本庆, 薛正辉, 任 武. FDTD计算中关于低频激励源问题的研讨[J]. 电波科学学报, 2009, 24(2): 213-217.

GAO Benqing, XUE Zhenghui, REN Wu. Study on low frequency exciting source in FDTD simulation [J]. Chinese Journal of Radio Science, 2009, 24(2): 213-217. (in Chinese)

[5] 董 慧, 李清亮, 闫玉波, 等. 三维大尺寸介质目标散射问题的总场边界PSTD技术[J]. 电波科学学报, 2005, 20(4): 434-439.

DONG Hui, LI Qingliang, YAN Yubo, et al. Application of TB-PSTD algorithm to the scattering of 3-D large dielectric objects[J]. Chinese Journal of Radio Science, 2005, 20(4): 434-439. (in Chinese)

[6] TAFLOVE A, HAGNESS S C. Computational Electrodynamics:The Finite-Difference Time-Domain Method[M]. 2nd ed. Boston: Artech House, 2000.

[7] 葛德彪, 闫玉波. 电磁波时域有限差分方法[M]. 2版. 西安电子科技大学出版社, 2006.

GE Debiao, YAN Yubo. The finite differrence time domain method of electromagnetic[M]. 2nd ed. Xi’an Electronic Science and Technology University Press, 2006. (in Chinese)

[8] 王长清. 现代计算电磁学基础[M].北京： 北京大学出版社, 2005.

WANG Changqing. Foudation of Computational Electromagnetism[M].Beijing: Beijing University Press, 2005. (in Chinese)

[9] 庄钊文, 袁乃昌, 莫锦军, 等. 军用目标雷达散射截面预估与测量[M]. 北京: 科学出版社, 2007: 163-164.

ZHUANG Zhaowen, YUAN Naichang, MO Jinjun, et al. Estimate and measurement of martial object’s RCS[M]. Beijing: Science Press, 2007: 163-164. (in Chinese)

[10] 谢彦召, 孙培云, 周 辉, 等. 地面附近高空核爆电磁脉冲环境[J]. 强激光与粒子束, 2003, 15(7): 680-684.

XIE Yanzhao, SUN Beiyun, ZHOU Hui, et al. High-altitude electromagnetic pulse environment over the lossy ground[J]. High Power Laser and Particle Beams, 2003, 15(7): 680-684. (in Chinese)

[11] 谢彦召, 王赞基, 王群书, 等. 高空核爆电磁脉冲波形标准及特征分析[J]. 强激光与粒子束, 2003, 15(8): 781-787 .

XIE Yanzhao, WANG Zanji, WANG Qunshu, et al. High altitude nuclear electromagnetic pulse waveform standards: a review [J]. High Power Laser and Particle Beams, 2003, 15(8): 781-787. (in Chinese)