一个二维的FDTD程序(总6页)

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*108; % 波速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*109); % 真空介电常数epsilon=1; %

2、相对介电常数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方向磁场 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=

3、1E; gi2(i)=1;gi3(i)=1;fi1(i)=0;fi2(i)=1.0;fi3(i)=1.0;endfor j=1:JE;gj2(j)=1;gj3(j)=1;fj1(j)=0;fj2(j)=1;fj3(j)=1;endfor 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;f

4、i1(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);endfor 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

5、/(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)+pu

6、lse; % 软源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;endfor i=1E;ez(i,1)=0;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

7、(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; 在Maxwell旋度方程的差分表示中，按照Yee氏的空间

8、网格设置，将出现半空间步长，通过前一时刻的磁、电场值得到时刻的电、磁场值，并在每一时刻上，将此过程算遍整个空间中随时间变化的电、磁场值的解，但在编程计算中不使用1/2的空间表示，而要通过一定的相互关系把它表达出来，在自制的C程序中采用数组来表示上面的表达式中各场值及系数，其表达式在程序中表示如下所示：ezk+1ij=CAij*ezkij+CBij*CD*(hyki+1j-hykij （6）+hxkij-hxkij+1)hxk+1ij+1=hxkij+1+CD*(ezkij-ezkij+1) （7）hyk+1i+1j=hyki+1j+CD*(ezki+1j-ezkij) （8） 由于各场值起始均赋零，时间步数从零开始，每一时间步均按上面的顺序在整个模拟区计算一遍，这样场的实际那关系和空间关系就完全被体现出来。首先我们看空间关系，比较（6）（8）和（3）（5）很容易看处，在（6）中的hyij实质上代表的是，即对于而言i实际上代表。比较（7）和（4）亦可看出，hxij+1实际上代表的是，即对于而言j1实际上代表。这样用此处的hxij+1,hxij,hyi+1j和hyij一起代入（6）中计算下一时间步中的电场值ezij。就是说虽然程序中没有出现网格的中间点，但这些点上的场值实际上被计算了出来。现在来考察时间关系，还是从磁场分量的计算开始，假