姬金祖

(北京航空航天大学 航空科学与工程学院，北京100191)

黄大庆

(中国航空工业集团公司北京航空材料研究院，北京100095)

黄沛霖 鲁振毅

(北京航空航天大学 航空科学与工程学院，北京100191)

飞行器隐身技术是现代战争中的一项重要军事技术，对于提高飞行器的战场生存力意义重大［1－2］.时域有限差分法(FDTD，Finite－Difference Time－Domain)是计算金属和复杂介质目标电磁散射、传播特性的一种重要方法，在隐身飞机概念设计、吸波材料设计等方面具有重要作用［3－5］.

用FDTD计算目标电磁散射的过程中，一般要把计算区域分成总场区和散射场区两个区域［6］.总场区的场量包括入射场和散射场，散射场区的场量只有散射场，没有入射场，因此由散射场区进行远场外推即可得到目标的雷达散射截面RCS(Radar Cross Section)［7］.总场区和散射场区的边界为连接边界，入射场通过连接边界引入到总场区.

散射体位于总场区内，当总场区不存在散射目标时，理想情况下散射场区的场量为0，总场区只有入射场.但在实际计算过程中，即使总场区没有散射目标，散射场区的场一般也不会为0，这是因为在连接边界引入入射场时的数值计算产生了误差，导致电磁泄漏.

文献［8］采用分裂平面波时域有限差分法(SP－FDTD，Splitting Plane Wave FDTD)引入入射场，在一维网格中将场量分裂，建立一维辅助网格和三维网格的色散关系.文献［9］在电磁散射计算中使用交错布置的方法引入入射波，该方法基于分裂场麦克斯韦方程组，扩展了多点辅助入射场的作用.文献［10］采用一维辅助多点模型，通过CN(Crank－Nicholson)差分方法引入入射波.该方法将二维色散方程投影到一维，使电磁泄漏在时间步长取CFL(Courant Friedrichs Levy)时达到很小.

实际计算表明，时间步长以及电磁波入射方向都对电磁泄漏有一定影响，本文研究这些参数对电磁泄漏的影响的大小.

1 FDTD计算区域划分

计算目标电磁散射过程中，FDTD计算区域一般划分为总场区和散射场区，在散射场区还设置辐射边界，用于将近场外推至远场，进一步计算RCS.FDTD计算区域划分如图1所示.

计算过程中，在连接边界处引入入射场，在此边界上的时间迭代公式要做一些修改，以保证该点周围的环量均为散射场或均为总场.时间迭代结束后，散射场区只有散射场.在辐射边界根据等效原理，将电磁场外推至远场，得到目标的RCS.

图1 FDTD计算区域划分

当总场区没有目标时，散射场区仍然会有电磁场，这是数值计算误差引入的电磁泄漏.本文对一维和二维情形采取不同的研究方法.一维情形下，用散射场区电场的平方和来衡量电磁泄漏的大小.在二维情形下，将散射场区的电磁场通过远场外推得到“无目标”下的远场RCS，用这种方法得到的RCS类似于微波暗室的背景电平.研究不同时间步长、入射角度下的远场RCS来研究电磁泄漏的程度.

2 一维情形电磁泄漏

2.1 差分方程

一维情形下电磁场微分方程具有最简单的形式，设电磁波传播方向为z，电场方向为x，磁场方向为y，则电磁波在真空中传播的微分方程如式(1)所示.

式中，Ex为电场强度;Hy为磁场强度;μ0为真空磁导率;ε0为真空介电常数.设c为真空中光速，令，结合真空中光速公式，将上式改写成较简单的形式，如式(2)所示.

式(2)显得更加简单、对称，本文涉及的电场、磁场也均采用上式简单、对称的形式.

用一阶Mur吸收边界条件在两个端点进行截断.采用一阶差分格式，对微分方程进行离散.将计算区域划分成100个网格，在距边界20个网格处设置连接边界条件.将ex在网格离散点处采样，hy在网格中点处采样，则得到101个ex离散点和100个 hy离散点，设离散点为 ex(i，n)和hy(i，n)，其中，i表示网格点编号，n表示时间步.设入射波为正谐波，波长为λ.

设网格长度为Δz，时间步长为Δt，定义空间网格因子为 fs=z/λ，时间网格因子为 ft=Δz/(cΔt).根据色散条件和稳定性条件的要求，fs应大于15，一维情形下ft应大于1.fs主要与色散误差相关，而本文主要研究时间网格因子ft的影响，故fs一律取为20.

用fs、ft表示真空中电磁波传播的差分格式，如式(3)所示.

该差分格式将介电常数、磁导率、网格尺寸、时间步长等参数用空间网格因子和时间网格因子进行归一化，形式更加简洁，便于编写程序和进行相应的计算.

2.2 时间网格因子影响

设电磁波传播方向为从左向右，设置不同时间网格因子，对比计算结果.运行1000个时间步后，两种情形下电场分布对比如图2所示.

图2 不同时间网格因子电场分布

由图2可见，ft的取值对电磁泄漏影响严重.ft=1时，电磁泄漏最小，散射场区电场强度基本为0.而ft=10时，电磁泄漏导致散射场区电场强度有较大起伏.

为定量研究ft对电磁泄漏的影响，以散射场区的电场平方和来表征泄漏的能量.取ft的值从1到10，间隔为1，分别计算其电磁泄漏能量.电磁泄漏能量随时间网格因子的变化如图3所示.

由计算结果可得，当 ft=1时，泄漏能量为 1.5179×10－29，非常接近 0，可以认为此时泄漏能量很小，可以忽略不计.由图3可得，泄漏能量基本随ft的增加而增加.当ft＜4时，泄漏能量增加很快，但当ft＞4时，能量泄漏增加很慢，基本保持在0.06 到0.07 之间.

图3 电磁泄漏能量随时间网格因子的变化

下面研究ft在1到2之间变化时，电磁泄漏的变化情况.ft步长取为0.1，计算结果如图4所示.

图4 电磁泄漏随时间网格因子的变化

由图4可见，时间网格因子在1到2之间变化时，电磁泄漏能量基本上线性增加.

可见，在一维FDTD中，时间网格因子在满足稳定性条件的基础上，取值应尽量小，即接近于1.

3 二维情形电磁泄漏

3.1 二维情形电磁泄漏比较

下面研究二维情形下的电磁泄漏.为保证稳定性，二维情形ft的下限是计算区域网格数为160×160，连接边界到吸收边界距离为20个网格，外推边界到吸收边界为10个网格.fs与一维情形相同，取为20.电磁波在真空中传播，无散射体时，不涉及到边界条件，故横磁波 TM(Transverse Magnetic)和横电波TE(Transverse Electric)没有本质区别.本文以TM波为例来进行分析研究.

入射角度定义为入射波来波方向与水平线夹角，逆时针为正，如图5所示.

ft取1.5 和10 两个值，φ 取0°和45°，计算不存在散射目标时的远场双站RCS，得到4条曲线.双站角为0°～359°，步长1°.将4条曲线在一张图上进行对比，如图6所示.二维RCS单位为m，在图中换算为对数值，单位为dBm

图5 入射角度示意图

图6 无目标下远场双站RCS

上面计算的RCS可以理解为类似微波暗室的“背景电平”.由上图可见，斜向45°照射时双站RCS明显小于正入射时双站RCS，一定程度上反映了斜向入射情况下电磁泄漏要小于正入射.将360个双站RCS值进行算术平均，RCS均值表如表1所示.

表1 双站RCS全向均值统计 dBsm

由表1可知，入射角45°下RCS远低于入射角0°，ft=1.5时 RCS远低于 ft=10.可见，ft与入射角对电磁泄漏影响很大，且由表1可知其差别甚至接近30 dB.

3.2 时间网格因子的影响

下面研究ft对电磁泄漏的影响.计算过程中取入射角45°，ft取值从 1.5 到10，步长 1.5，其他参数同前.RCS均值随ft的变化曲线如图7所示.

图7 RCS均值随ft变化

由图7可知，ft从1.5 变化到2.5，RCS均值迅速增加，而ft超过2.5之后，RCS均值的变化区域平缓，逐渐趋于一个固定的值，最后维持在约－18 dB左右.ft取值接近稳定性条件，可使远场RCS降低约30 dB.

3.3 入射角度影响

二维情形下电磁泄漏与入射角度也密切相关.下面研究随入射角度变化的影响，网格划分同上.由于对称性，入射角取为0°～45°就能够代表所有的情形.入射角步长取为5°.

用远场外推方法计算双站RCS，双站RCS均值随入射角的变化如图8所示.

图8 双站RCS均值随入射角度的变化

由图8可见，双站RCS均值随入射角度增大而减小.在入射角较小时，双站RCS均值随入射角改变比较缓慢，入射角大于15°后，双站RCS均值改变逐渐加快.当入射角为45°时，达到最小值，约－37 dB.不同的入射角双站RCS均值差别可达30 dB.

4 结论

本文用不同的衡量方法对一维、二维情形下连接边界导致的电磁泄漏进行了研究，主要结论如下:

1)为减小电磁泄漏，时间步长在满足稳定性要求的前提下，应越小越好.在二维情形下，时间步长过大引起的电磁泄漏可比时间步长接近稳定性条件时大30 dB.电磁泄漏随时间步长变化的趋势先慢后快，时间网格因子小于2时，基本呈线性关系.时间网格因子大于3后，电磁泄漏随时间步长的增加缓慢变化，趋于稳定

2)入射角对电磁泄漏也有很大影响.当入射波斜入射到连接边界时，不同的入射角对电磁泄漏也有很大影响.当入射角与计算边界成斜角时，电磁泄漏较小.对于正方形区域，入射角度45°时电磁泄漏最小.垂直入射和斜入射时，电磁泄漏随入射角度变化较为缓和，其他角度变化较剧烈.

References)

［1］郑奎松，葛德彪，魏兵.导弹目标的FDTD建模与RCS计算［J］.系统工程与电子技术，2004，26(7):896－899 Zheng Kuisong，Ge Debiao，Wei Bing.FDTD modeling of missile target and RCS computation［J］.Systems Engineering and Elec－tronics，2004，26(7):896－899(in Chinese)

［2］李军，武振波，武哲.稳定因子对FDTD数值计算的影响［J］.北京航空航天大学学报，2004，30(1):70－73 Li Jun，Wu Zhenbo，Wu Zhe.Effects of stability factor on FDTD numerical calculation［J］.Journal of Beijing University of Aeronautics and Astronautics，2004，30(1):70－73(in Chinese)

［3］黄志祥，沙威，吴先良，等.辛FDTD 算法［J］.系统工程与电子技术，2009，31(2):456－458 Huang Zhixiang，Sha Wei，Wu Xianliang，et al.Scheme of symplectic FDTD［J］.Systems Engineering and Electronics，2009，31(2):456－458(in Chinese)

［4］胡晓娟，卢兆林，葛德彪，等.基于三角面元的涂层目标FDTD共形网格生成技术［J］.系统工程与电子技术，2010，32(9):1884－1888 Hu Xiaojuan，Lu Zhaolin，Ge Debiao，et al.Conformal FDTD mesh－generating scheme for coated targets based on trianglepatch［J］.Systems Engineering and Electronics，2010，32(9):1884－1888(in Chinese)

［5］葛德彪，闫玉波.电磁波时域有限差分方法［M］.2版.西安:西安电子科技大学出版社，2005 Ge Debiao，Yan Yubo.Finite－difference time－domain method for electromagnetic waves［M］.2nd ed.Xi′an:Xidian University Press，2005(in Chinese)

［6］王长清.现代计算电磁学基础［M］.北京:北京大学出版社，2004 Wang Changqing.Basic of modern computational electromagnetics［M］.Beijing:Peking University Press，2004(in Chinese)

［7］姜彦南，葛德彪，魏兵.时域有限差分并行算法中的吸收边界研究［J］.系统工程与电子技术，2008，30(9):1636－1640 Jiang Yannan，Ge Debiao，Wei Bing.Study on absorbing boundary condition in parallel FDTD algorithm［J］.Systems Engineering and Electronics，2008，30(9):1636－1640(in Chinese)

［8］Wang Hui，Huang Zhixiang，Wu Xianliang，et al.Perfect plane wave injection into 3D FDTD(2，4)scheme［C］//Cross Strait Quad－Regional Radio Science and Wireless Technology Conference(CSQRWC)Volume 1.Harbin:IEEE，2011:40－43

［9］Hadi M F.A versatile split－field 1－D propagator for perfect FDTD plane wave injection［J］.IEEE Transactions on Antennas and Propagation，2009，57(9):2691－2697

［10］Huang Z，Pan G，Pan H K.Perfect plane wave injection for Crank－Nicholson time－domain method［J］.Microwaves，Antennas ＆ Propagation，IET，2010，4(11):1855－1862