二维FDTD有限元仿真

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=1:IE+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=2:IE; % 
for j=2:JE; 
ga(i,j)=1; 
end; 
end; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%PML参数的设置 
for i=1:IE; 
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=2:IE; % 为了使每个电场周围都有磁场进行数组下标处理 
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=1:IE; % 为了使每个电场周围都有磁场进行数组下标处理 
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; % 电荷密度循环结束 
%figure(1);                                %可视化处理
      %clf;   
      %mesh(ez);                                 %电场的幅值
      %axis([0 IE 0 JE -0.5 0.5]);
      %xlabel('i')
      %ylabel('j')
  %drawnow;                %完成114_ez
for j=1:JE; 
ez(1,j)=0; 
ez(IE,j)=0; 
end 
for i=1:IE; 
ez(i,1)=0; 
ez(i,JE)=0; 
end; 
for i=1:IE; % 为了使每个磁场周围都有电场进行数组下标处理 
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循环结束 
      figure(2);                                %可视化处理
      clf;   
      mesh(hx);                                 %磁场的幅值
      axis([0 IE 0 JE -0.5 0.5]);
      xlabel('i')
      ylabel('j')
      drawnow;                %完成114_hx
for i=1:IE-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循环结束 
figure(3);                                %可视化处理
          %clf;   
          %mesh(hy);                                 %磁场的幅值
          %axis([0 IE 0 IE -0.5 0.5]);
          %xlabel('i')
          %ylabel('j')
          %drawnow;              %完成114_hy
end 
 

 

 D139