一个二维的FDTD程序

一个二维的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方向磁场

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

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)-h

y(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;

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;

在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]实质上代表的是),21(j i H y -，即对于y H 而言i 实际上代表2

1-i 。比较（7）和（4）亦可看出，hx[i][j+1]实际上代表的是)2

1,(+j i H x ，即对于x H 而言j ＋1实际上代表2

1+j 。这样用此处的hx[i][j+1],hx[i][j],hy[i+1][j]和hy[i][j]一起代入（6）中计算下一时间步中的电场值ez[i][j]。就是说虽然程序中没有出现网格的中间点，但这些点上的场值实际上被计算了出来。现在来考察时间关系，还是从磁场分量的计算开始，假设现在是第k 步，则在（6）－（8）

中hy[k][i][j]实际上表示的为),(2

1

j i H k y +，hx[k][i][j]实际上表示的为),(21j i H k x +，式中对于y H ，x H 而言，实际上k 表示2

1+k ,首先由（7）和（8）推算出hx[k+1][i][j]，hy[k+1][i+1][j]，然后把这两式代入下一时间步中，由式（6）便可计算出ez[k+1][i][j]的值。这就完成了各场量的时间递推关系。