时域有限差分法(FDTD算法)的基本原理及仿真

1、时域有限差分法（FDTD算法）时域有限差分法是1966年K.S.Yee发表在AP上的一篇论文建立起来的，后 被称为Yee网格空间离散方式。这种方法通过将Maxwell旋度方程转化为有限差 分式而直接在时域求解，通过建立时间离散的递进序列，在相互交织的网格空 间中交替计算电场和磁场。FDTD算法的基本思想是把带时间变量的 Maxwell旋度方程转化为差分形 式，模拟出电子脉冲和理想导体作用的时域响应。需要考虑的三点是差分格式、 解的稳定性、吸收边界条件。有限差分通常采用的步骤是：采用一定的网格划分 方式离散化场域；对场内的偏微分方程及各种边界条件进行差分离散化处理，建立差分格式，得到差分方程组；

2、结合选定的代数方程组的解法，编制程序，求边 值问题的数值解。1. FDTD的基本原理FDTD方法由Maxwell旋度方程的微分形式出发，利用二阶精度的中心差分 近似，直接将微分运算转换为差分运算，这样达到了在一定体积内和一段时间上 对连续电磁场数据的抽样压缩。Maxwell方程的旋度方程组为：:h(1)g E =卩CT Hm1 1t在直角坐标系中，（1）式可化为如下六个标量方程:汨zHyJa：Ex：:yz；:tHx：Hz ：E y:zt:x；：tHy：Hx ：Ez:xy；：tE2Ey，:y-Ex.zEy:x.:E：Ez.x.:E；:H_=t；=H-:t(2)-:t上面的六个偏微分方程是FDTD

3、算法的基础。 Yee首先在空间上建立矩形差分网格，在时刻n t时刻，F（x,y,z）可以写成F(x,y,z,t)二 F(i x, j y,k z,n=t Fn(i, j,k)(3)用中心差分取二阶精度: 对空间离散：F(x,y,z,t):xFn(i 12,j,k)-Fn(i -12,j,k) oLJlx：F(x,y, z,t)Fn(i,j 1 2,k)-Fn(i,j -1 2,k)：F(x,y,z,t) zFn(i,j,k12)-Fn(i,j,k-12) . o，：z25对时间离散:F(x, y, z,t):tFn12(ijk)-Fn2(ijk) fl(4)Yee把空间任一网格上的E和H的六个

4、分量，如下图放置图1 Yee氏网格及其电磁场分量分布在FDTD中，空间上连续分布的电磁场物理量离散的空间排布如图所示。 由图可见，电场和磁场分量在空间交叉放置，各分量的空间相对位置也适合于 Maxwell方程的差分计算，能够恰当地描述电磁场的传播特性。同时，电场和磁 场在时间上交替抽样，抽样时间间隔相差半个时间步，使Maxwell旋度方程离散 以后构成显式差分方程，从而可以在时间上迭代求解，而不需要进行矩阵求逆运 算。因此，由给定相应电磁问题的初始条件，FDTD就可以逐步推进地求得以后 各个时刻空间电磁场的分布。根据这一原则可以写出六个差分方程：i+1/2,j,k2t讥 FXTOjt 尼(iT

5、/2jk)2 (i 1/2, j,k)+1( 5)；(i 1/2 jk) .1/2 j,k) t2叩+1/2, j,k)卢丁/20+1/2*)-Hj/2(i+1/2,j -1/2,k)丄 H2(i+1/2j,k-1/2)-出乜勺 + “冇飞牛 yz其余的也如法可以写出，每个网格点上的个场分两的新值依赖于该点在前一 时间步长时刻的值机该点周围的临近点上另一场量在早半个时间步长时的值。因此任一时刻可一次算出一个点，并行算法可计算出多个点。通过这些运算可以交 替算出电场磁场在各个时间步的值。根据上述FDTD差分方程组可得出计算电磁场的时域推进计算方法， 如图2已知t1 =tn ：t 0时刻空间各处的

6、电磁场初始所示。图2 FDTD在时域的交叉半步逐步推进计算2.数值稳定性条件时间步长氏，空间步长.=x，：y，z必须满足一定的关系，否则就使得数值表现不稳定，表现为：随着计算步数的增加，计算场量的数值会无限的增大， 这种增大不是由于误差积累造成的，而是由于电磁波的传播关系被破坏造成的。所以.述，x，細，.:z必须满足一定的关系以保证稳定性。Taflove等在1975年对Yee氏差分格式的稳定性进行了讨论，并导出了对时间步长的限制条件。数 取决于时间步长 厶t与空间步长 x、诃 z的关系。对于非:1 : 1定主匕要均匀媒质构成的计算空间选用如下的稳定性条件:=(:)2(6)匚s若米用均匀立方体网

7、格：泳=诃二厶z二厶s， t -:W3(7)而一般取：.芸，c为光速。当x， y， z不相等时，2c(8)3.数值色散FDTD网格中，会导致数字波模在网格中发生改变，这种改变是由于计算网 格本身引起的，而非物理因素，所以必须考虑。即在 FDTD网格中，电磁波的 相速与频率有关，电磁波的相速度随波长、传播方向及变量离散化的情况不同而 改变。色散将导致非物理因素引起的脉冲波形畸变、人为的各向异性和虚假折射 等现象。显然，色散与空间、时间的离散间隔有关，如下式所示：1. 2偽攻1sin I二丄 sin2匚 sin2厘+亠卅込 (9)徴广气2八。气2丿与数值色散关系相对应，在无耗介质中的单色平面波，色

8、散解析关系是:由式（9）可知，当式（9）中的述、x、厶y、厶z均趋于零时，它就趋于 式（10）。也就是说数值色散是由于用近似差分替代连续微分而引起的，而且在 理论上可以减小到任意程度，只要此时时间步长和空间步长都足够小。为获得理 想的色散关系，问题空间分割应按照小于正常网格的原则进行。一般选取的最大 空间步长为厶maxm20，韦n为所研究范围内电磁波的最小波长。由上分析 说明，数值色散在用 FDTD法分析电磁场传播中的影响是不可能避免的，但我 们可以尽可能的减小数值色散的影响。现在适当选取时间和空间步长，传播方向，可以得到理想情况，如下所示：3-D方形网格：（数值稳定的极限状态，可得理想色散关

9、系）取波沿对角线传播 kx 二 ky 二 kz = k/ 3， _x = y = z = t （ 11）L曲2-D方形网格：也是沿对角线传播kkkk/ . 2，t=二（1-D网格: 4.吸收边界条件在电磁场的辐射和散射问题中，边界总是开放的，电磁场占据无限大空间， 而计算机内存是有限的，所以只能模拟有限空间。即:时域有限差分网格将在某 处被截断。这要求在网格截断处不能引起波的明显反射，因而对向外传播的波而言，就像在无限大的空间传播一样，一种行之有效的方法是在截断处设置一种吸 收边界条件。使传播到截断出的波被边界吸收而不产生反射。下面只给出Engquist-Majda吸收边界条件，采用 Mur差

10、分格式，其总体虚 假反射在1%5%之间。一维一阶近似情形，x=0 边界：un 1（0） = uW） u0 1（1） - u（0） （ 14）二维二阶近似情形,x=0边界:wn1(0，j)m wsmWn(0, j -12Wn(0, j) Wn(0, j -1) Wn(1,j 1) -2Wn(1, j) Wn(1, j -1)三维二阶近似情形，x=0边界：Wn#(O,j,k) = W七 j, k) +CXWn%, j,k) W(0, j,k)cit x242Wn(0，j，k) Wn(1,j,k)(：)+咲(c ：t ： -x2( ：y)2(c ：t ： -x)(c ：t)2 x2( z) 2(c t x)Wn(0, j +1,k)如(0, j,k)