二維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