Hx(i,j,k)?CP(i,j,k)Hx(i,j,k)?CQ(i,j,k)?EZ(i,j?1,k)?EZ(i,j,k)Ey(i,j,k?1)?Ey(i,j,k)?(2-8d)
?????y?z??Hy(i,j,k)?CP(i,j,k)Hx(i,j,k)?CQ(i,j,k)?E(i,j,k?1)?Ex(i,j,k)EZ(i?1,j,k)?EZ(i,j,k)?(2-8e)
??x???z?x??HZ(i,j,k)?CP(i,j,k)HZ(i,j,k)?CQ(i,j,k)?Ey(i?1,j,k)?Ey(i,j,k)Ex(i,j?1,k)?Ex(i,j,k)?(2-8f)
?????x?y??根据上述FDTD差分方程组可得出计算电磁场的时域推进计算方法，如图2-2所示。
图2-2 FDTD在时域的交叉半步逐步推进计算
已知t1=t0=n?t 0时刻空间各处的电磁场初始值 计算t2=t1+?t2 时刻空间各处的磁场值 循环n次 计算t1=t2+?t2 时刻空间各处的电场值 式(2-8a)~(2-8c)的等号左边的电场值是第n次循环的电场值，等号右边的电场值是第n-1次循环存储在内存中的电场值，磁场值是本次循环计算得到的磁场值；式(2-8d)~(2-8f) 等号左边的磁场值是第n次循环的磁场值，等号右边的磁场值和电场值都是第n-1次循环存储在内存中的场值。这样，就解决了半个时间步在程序中无法表示的问题，而且也没有破坏电磁场在时间上逐步推进的逻辑关系。
2.3 时域有限差分法相关技术 2.3.1 数值稳定性问题
上述FDTD方程是一种显式差分方程，在执行时，存在一个重要的问题：即算法的稳定性问题。这种不稳定性表现为在解显式方程时，随着时间步数的继续
增加，计算结果也将无限制地增加。Taflove等[4]于1975年对Yee氏差分格式的稳定性进行了讨论，并导出了对时间步长的限制条件。数值解是否稳定主要取决于时间步长?t与空间步长?x、?y、?z的关系。对于非均匀媒质构成的计算空间选用如下的稳定性条件：
?t?1 (2-9)
111v()2?()2?()2?x?y?z(2-9)式是空间和时间离散之间应当满足的关系，又称为Courant稳定性条件。若采用均匀立方体网格：?x??y??z??s
?s?t? (2-10)
v3其中，v为计算空间中的电磁波的最大速度。
2.3.2 数值色散
FDTD方程组是对Maxwell旋度方程进行差分近似，在进行数值计算时，将
会在计算网格中引起数字波模的色散，即在FDTD网格中，电磁波的相速与频率有关，电磁波的相速度随波长、传播方向及变量离散化的情况不同而改变。这种关系由非物理因素引起，且色散将导致非物理因素引起的脉冲波形畸变、人为的各向异性和虚假折射等现象。显然，色散与空间、时间的离散间隔有关，如下式所示：
k?y?111112???t?2?kx?x?2?y2?kz?z?2?kz?z? (2-11) ??sin?sin?sin?sin?sin????????2???z?2222?c?t?2?2???x?2?2???y?2????z??????式(2-11)是三维情况下在FDTD方法中的单色平面波数值色散关系的一般形式，它表明FDTD计算中波的传播速度与传播方向有关。式中kx、ky、kz分别是波矢量沿x、y、z方向的分量，?是角频率，c是被模拟的均匀介质中的光速。与数值色散关系相对应，在无耗介质中的单色平面波，色散解析关系是：
??c?2?kx2?ky2?kz2 (2-12)
由式(2-11)可知，当式(2-11)中的?t、?x、?y、?z均趋于零时，它就趋于式(2-12)。也就是说数值色散是由于用近似差分替代连续微分而引起的，而且在理论上可以减小到任意程度，只要此时时间步长和空间步长都足够小，但这将大大增加所需的计算机存储空间和计算时间，并使累积误差增加。因此，在实际计算中要根据问题的性质和计算机的软硬件条件来选择合适的时间步长和空间步长。为获得理想的色散关系，问题空间分割应按照小于正常网格的原则进行。一般选取的最大空间步长为?max??min20，?min为所研究范围内电磁波的最小波长。由上分析说明，数值色散在用FDTD法分析电磁场传播中的影响是不可能避免的，但我们可以尽可能的减小数值色散的影响。
2.3.3 离散网格的确定
无论是简单目标还是复杂目标，在进行FDTD离散时网格尺寸的确定，除了
受计算资源的限制不可能取得很小外，还需要考虑以下几个因素：
1．目标离散精确度的要求。网格应当足够小以便能精确模拟目标几何形状和电磁参数。
2．FDTD方法本身的要求。主要是考虑色散误差的影响。设网格为立方体
?x??y??z??，所关心频段的频率上限为fmax，对应波长为?min，则考虑
FDTD的数值色散要求
???minN （2-13）
通常N?10。上式是根据已知所关心频率上限情况下来确定FDTD网格尺寸?的；反之，若给定?，则FDTD计算结果可用的上限频率也随之确定。
3．入射波的要求。入射波的上限截止频率fc应包含所关心频率范围，即
fc?fmax。
2.4 吸收边界条件
由时域有限差分法的基本原理可知，在利用时域有限差分法研究电磁场时，
需在全部问题空间建立Yee氏网格空间，并存储每个单元网格上任一时间步的六个场分量用于下一时间步的计算。而在对于辐射、散射这类开放系统的实际研究
中，不可能有无限大的存储空间。因此，必须在某处将网格空间截断，且在截断边界网格点处运用特殊的场分量计算方法，使得向边界面行进的波在边界处保持“外向行进”特性、无明显的反射现象，并且不会使内部空间的场产生畸变，从而用有限网格空间模拟电磁波在无界空间中传播的情况。具有这种功能的边界条件称之为吸收边界条件，或辐射边界条件，或网格截断条件[29~31]，如图2-3所示。
截断边界散射体入射波散射波
图 2-3 附加截断边界使计算区域变为有限域
从FDTD的基本差分方程组可以看出，在截断边界面上切向场分量的计算需要利用计算空间以外的电磁场分量，因此FDTD基本差分方程对这些截断边界面上的场分量失效。如何处理截断边界上的场分量，使之与需要考虑的无限空间有尽量小的差异，是FDTD中必须很好解决的一个重要问题。实际上，这是要求在误差可容忍的范围内，计算空间中的外向波能够顺利通过截断边界面而不引起波的明显反射，使有限计算空间的数值模拟与实际情况趋于一致，对外向波而言，就像在无限大空间中传播一样。所以，需要一种截断边界网格处的特殊计算方法，它不仅要保证边界场计算的必要精度，而且还要大大消除非物理因素引起的波反射，使得用有限的网格空间就能模拟电磁波在无限空间中的传播。但是如果处理不当，截断边界面可能造成较大反射，构成数值模拟误差的一部分，甚至可能造成算法不稳定。
加于截断边界场分量符合上述要求的算法就称为吸收边界条件（Absorbing Boundary Conditions）。
2.4.1 一阶和二阶近似吸收边界条件
在截断边界附近通常没有激励源。考虑齐次波动方程
??2?2?21?2????x2??y2??z2?c2?t2??f?0 (2-14) ??式中，f表示直角坐标系下任意电磁场分量。
B.Engquist和A.Majda[15]利用偏微分算子对式(2-13)作因式分解，并分别取其Taylor级数展开式中的第一项和前两项近似，导出了适合直角坐标系下FDTD吸收边界条件的单向波动方程，这就是Engquist-Majda吸收边界条件。设三维长方体FDTD区域0
一阶近似 二阶近似 x?0 ??1?????f??xc?t?x?0?0 ?0 ?1?21?21??2?2???22???2???f2??2??y?z???c?x?tc?t?1?21?21??2?2???22???2???f2??2??y?z???c?x?tc?t?1?21?21??2?2???c?y?t?c2?t2?2??z2??x2??f????x?0?0 ?0 ?0 x?a ??1?????f??xc?t?x?ax?ay?0 ??1?????y?c?t??f??y?b ??1?????y?c?t??f??z?0 ??1?????f??zc?t?z?d ??1?????f??xc?t?
y?0?0 ?0 ?0 ?0 y?0y?b?1?21?21??2?2???c?y?t?c2?t2?2??z2??x2??f?????1?21?21??2?2???22???2???f2??2??y?x???c?z?tc?ty?b?0 ?0 z?0z?0z?d?1?21?21??2?2???22???2???f2??2??y?x???c?z?tc?tz?d?0