??Hyx?t??mxHyx??(Ezx?Ezy)?x (2-22j)
??(Eyz?Eyx)?Hzx??mxHzx?? (2-17k) ?t?x??Hzy?t??myHzy??(Exy?Exz)?y (2-22l)
式中?x，?y，?z，?mx，?my，?mz为电导率和磁导率，描述了PML介质的各向异性。当?x??y??z??且?mx??my??mz??m时， PML介质退化为普通有耗介质;当??0且?m?0时，退化为自由空间中的Maxwell方程。同时，要满足PML介质的重要基本条件即阻抗匹配条件，如下式所示：
??m? (2-23) ?0?0式中?0和?0分别为真空的介电常数和真空的磁导率。
对式(2-18)做差分处理，可得到PML介质中的FDTD差分方程表达式：
Exy(i?12,j,k)?CAy(j)Exy(i?12,j,k)?CBy(j) ?HZn?12n?1n(i?12,j?12,k)?HZ?yn?12(i?12,j?12,k)(2-23a)
n?1nExz(i?12,j,k)?CAz(k)Eyz(i,j?12,k)?CBz(k)n?12Hx(i,j?12,k?12)?Hxn?12(i,j?12,k?12)(2-23b)
??zn?1nEyz(i,j?12,k)?CAz(k)Exy(i,j?12,k)?CBz(k)n?12Hxn?12(i,j?12,k?12)?HZ(i,j?12,k?12)(2-23c)
??zn?1nEyx(i,j?12,k)?CAx(i)Eyx(i,j?12,k)?CBx(i)n?12n?12Hz(i?12,j?12,k)?HZ(i?12,j?12,k)(2-23d)
??xn?1nEzx(i,j,k?12)?CAx(i)Ezx(i,j,k?12)?CBx(i)n?12Hy(i?12,j,k?12)?Hyn?12(i?12,j,k?12)(2-23e)
??xn?1nEzy(i,j,k?12)?CAy(j)Ezy(i,j,k?12)?CBx(j)n?12Hx(i,j?12,k?12)?Hxn?12(i,j?12,k?12) (2-23f)
??yn?12n?12Hxy(i,j?12,k?12)?DAy(j?12)Hxy(i,j?12,k?12)Ezn(i,j?1,k?12)?Ezn(i,j,k?12) (2-23g)
?DBy(j?12)??yn?12n?12Hxz(i,j?12,k?12)?DAz(k?12)Hxz(i,j?12,k?12)nEy(i,j?12,k?1)?Eyn(i,j?12,k) (2-23h)
?DBz(k?12)??zn?12n?12Hyz(i?12,j,k?12)?DAz(k?12)Hyz(i?12,j,k?12)Exn(i?12,j,k?1)?Exn(i?12,j,k) (2-23i)
?DBz(k?12)??zn?12n?12Hyx(i?12,j,k?12)?DAx(i?12)Hyx(i?12,j,k?12)nEz(i?1,j,k?12)?Ezn(i,j,k?12) (2-23j)
?DBx(i?12)??xHzxn?12(i?12,j?12,k)?DAx(i?12)Hzxnn?12(i?12,j?12,k)n ?DBx(i?12)?Ey(i?1,j?12,k)?Ey(i,j?12,k) (2-23k)
?xn?12n?12Hzy(i?12,j?12,k)?DAy(j?12)Hzy(i?12,j?12,k)nEx(i?12,j?1,k)?Exn(i?12,j,k) (2-23l)
?DBy(j?12)??y式中
CAq(?)?2???q(?)?t2???q(?)?tt，CBq(?)?2?t (2-24) 2???q(?)?ttDAq(??12)?2???m(??12)?t，2?t (2-25) DBq(??12)?2???m(??12)?t2???m(??12)?t其中，q?x,y,z；??i,j,k
在实际计算中，在FDTD计算中PML的设置不可能延伸到半无限空间，只能是有限厚度，因此PML的外侧边界需要特殊处理，通常情况下采用理想导体边界截断。这样，透入PML中的外向波到达理想导体边界处会反射回来，重新进入计算区域，PML的反射系数不再等于零。PML介质内沿q方向的电导率分布通常采取以下函数形式
?r? ?(r)???qqma?x???m （2-26）
式中，?为PML层的厚度，r为PML层靠近FDTD分界面的距离，m是电导率分布阶数，表示PML中电导率变化的程度，取整数。式（2-26）说明了实际计算中PML介质层厚度必须取若干个空间步长，使电导率从FDTD-PML分界面的0渐变到PML最外面的?qmax，避免电导率跃变太大，尽量消除数值反射。这样，如果相对于FDTD-PML分界面定义的外向波入射角为?，则PML内侧表面反射系数为
?2??qmaxcos?R????exp??m?1??0?? （2-27） ?0?使用PML吸收边界条件时，首先要选定三个参数：PML层数N，电导率分布阶数m，PML表面反射系数R(?)。大量的数值实验表明：（1）当层数N固定时，减小R(?)，即增加PML的衰减，可以使局部及总体误差都单调地减小。然而，当R(?)减小到一定程度后，这种现象不再出现，原因是存在由空间网格引起的固有误差。（2）增加PML层数N，可以使局部及总体误差都单调地减小，但是N过大又会使计算量剧增，需折衷考虑吸收效果和计算量。（3）电导率分布阶数m的改变不会影响计算量，却会影响PML的吸收效果，一般情况下，m越大，吸收效果越好。因此，在实际计算中参数的选择随所处理问题的不同而不同，需综合考虑，并由数值实验寻找最佳值。
2.5 FDTD计算所需时间步的估计
为了使计算达到稳定,通常计算所需要时间步按照电磁波往返穿越FDTD计算
区对角线3～5次来估计。若FDTD计算区总元胞数为N，则对角线上元胞约为
3N13 （三维)。按照Courant稳定条件，设计算中心区?t??(2c)，即穿
越对角线一次需要时间步为213总计算时间步约需23??6~10?N步。3N13。
12对于二维情况则约为22??6~10?N。或者说，计算时间步大约等于FDTD
计算区对角线上元胞数目的12～20倍。实际上，计算所需时间步还与散射体具体形状、结构有关。
图2-9给出了应用FDTD分析电磁场问题时的程序流程图
变量定义变量赋值 初始化波导建模加入激励源计算磁场HzStep=step+1边界和角点的处理计算电场Ex,EyStep=Max step?YES输出结果NO
图2-9 FDTD 程序流程图
第三章 MATLAB的仿真的程序及模拟
3.1 MATLAB程序及相应说明
% 二维FDTD TE波仿真 clear all; % 定义常数 pi=3.1415;
c=3.0e10; %高斯制下光速 f=1.0e15; %频率 lambda=c/f; %波长 nmax=400; %时间步数
del_s=lambda/20; %每最小波长20个采样点 del_t=0.5*del_s/c; %迭代时间步长 n=182; %真空区域网格数 np=9; %pml层数 N1=n+2*np; %总网格数 N=N1+1; %采样点数 M=4; %导电率渐变指数 sigma_max=(M+1)/1.50/pi/del_s; %最大导电率
% TE波的分量初始化 tic; figure(1);
axis([0 N 0 N -0.5 0.5]);
Ex=zeros(N1,N); %x方向为横向，采样点为网格的横向边，故行数+1 Ey=zeros(N,N1); %y方向为纵向，采样点为网格的纵向边，故列数+1 Bz=zeros(N1,N1); %矩阵行为纵向网格数，矩阵列为横向网格数,循环中用j表示行数，i表示列数 Bzx=zeros(N1,N1);