一个二维的FDTD程序

1 一个二一个二维维的的 FDTD 程序程序 本程序实现 2 维 TM 波 FDTD 仿真 此程序用 PML 设置吸收边界条件 FDTD 2D kongqi PML 仅含有 Ez Hx Hy 分量 clear clc 1 初始化 T 200 迭代次数 IE 100 JE 100 npml 8 PML 的网格数量 c0 3 10 8 波速 f 1 5 10 9 频率 lambda c0 f 波长 wl 10 dx lambda wl dy lambda wl pi 3 14159 dt dx 2 c0 时间间隔 epsz 1 4 pi 9 10 9 真空介电常数 epsilon 1 相对介电常数 sigma 0 电导率 spread 6 脉冲宽度 t0 20 脉冲高度 ic IE 2 源的 X 位置 jc JE 2 源的 Y 位置 for i 1E 1 for j 1 JE 1 dz i j 0 z 方向电荷密度 ez i j 0 z 方向电场 hx i j 0 x 方向磁场 hy i j 0 y 方向磁场 2 ihx i j 0 ihy i j 0 iz i j 0 z 方向求和参量 频域卷积转化为时域求和 end end for i 2E for j 2 JE ga i j 1 end end PML 参数的设置 for i 1E gi2 i 1 gi3 i 1 fi1 i 0 fi2 i 1 0 fi3 i 1 0 end for j 1 JE gj2 j 1 gj3 j 1 fj1 j 0 fj2 j 1 fj3 j 1 end for i 1 npml 1 设置 PML 层中的参数 xnum npml 1 i xn 0 33 xnum npml 3 gi2 i 1 0 1 xn gi2 IE 1 i 1 1 xn gi3 i 1 xn 1 xn gi3 IE 1 i 1 xn 1 xn xn 0 25 xnum 0 5 npml 3 fi1 i xn fi1 IE 2 i xn fi2 i 1 0 1 xn fi2 IE 2 i 1 1 xn fi3 i 1 xn 1 xn fi3 IE 2 i 1 xn 1 xn end 3 for i 1 npml 1 xnum npml 1 i xn 0 33 xnum npml 3 gj2 i 1 0 1 xn gj2 JE 1 i 1 1 xn gj3 i 1 xn 1 xn gj3 JE 1 i 1 xn 1 xn xn 0 25 xnum 0 5 npml 3 fj1 i xn fj1 JE 2 i xn fj2 i 1 0 1 xn fj2 JE 2 i 1 1 xn fj3 i 1 xn 1 xn fj3 JE 2 i 1 xn 1 xn end 2 迭代求解电场和磁场 for t 1 T for i 2E 为了使每个电场周围都有磁场进行数组下标处理 for j 2 JE dz i j gi3 i gj3 j dz i j gi2 i gj2 j 0 5 hy i j hy i 1 j hx i j hx i j 1 end end 电场循环结束 pulse sin 2 pi f t dt 正弦波源 dz ic jc dz ic jc pulse 软源 for i 1E 为了使每个电场周围都有磁场进行数组下标处理 for j 1 JE ez i j ga i j dz i j 反映煤质的情况都是放到这里的 iz i j iz i j gb i j ez i j end end 电荷密度循环结束 for j 1 JE ez 1 j 0 ez IE j 0 end for i 1E ez i 1 0 4 ez i JE 0 end for i 1E 为了使每个磁场周围都有电场进行数组下标处理 for j 1 JE 1 curl e ez i j ez i j 1 ihx i j ihx i j fi1 i curl e hx i j fj3 j hx i j fj2 j 0 5 curl e ihx i j end end 磁场 HX 循环结束 for i 1E 1 为了使每个磁场周围都有电场进行数组下标处理 for j 1 JE curl e ez i 1 j ez i j ihy i j ihy i j fj1 j curl e hy i j fi3 i hy i j fi2 i 0 5 curl e ihy i j end end 磁场 HY 循环结束 end end 5 在 Maxwell 旋度方程的差分表示中 按照 Yee 氏的空间网格设置 将出现 半空间步长 通过前一时刻的磁 电场值得到时刻的电 磁场值 并在每一时 刻上 将此过程算遍整个空间中随时间变化的电 磁场值的解 但在编程计算 中不使用 1 2 的空间表示 而要通过一定的相互关系把它表达出来 在自制的 C 程序中采用数组来表示上面的表达式中各场值及系数 其表达式在程序中表 示如下所示 ez k 1 i j CA i j ez k i j CB i j CD hy k i 1 j hy k i j 6 hx k i j hx k i j 1 hx k 1 i j 1 hx k i j 1 CD ez k i j ez k i j 1 7 hy k 1 i 1 j hy k i 1 j CD ez k i 1 j ez k i j 8 由于各场值起始均赋零 时间步数从零开始 每一时间步均按上面的顺序 在整个模拟区计算一遍 这样场的实际那关系和空间关系就完全被体现出来 首先我们看空间关系 比较 6 8 和 3 5 很容易看处 在 6 中的 hy i j 实质上代表的是 即对于而言 i 实际上代表 2 1 jiHy y H 2 1 i 比较 7 和 4 亦可看出 hx i j 1 实际上代表的是 即对于 2 1 jiHx 6 而言 j 1 实际上代表 这样用此处的 hx i j 1 hx i j hy i 1 x H 2 1 j j 和 hy i j 一起代入 6 中计算下一时间步中的电场值 ez i j 就是说 虽然程序中没有出现网格的中间点 但这些点上的场值实际上被计算了出来 现在来考察时间关系 还是从磁场分量的计算开始 假设现在是第 k 步 则在 6 8 中 hy k i j 实际上表示的为 hx k i j 实际上表 2 1 jiH k y 示的为