24小时热门版块排行榜

返回列表

【悬赏金币】回答本帖问题，作者lfm3041将赠送您 8 个金币

lfm3041

金虫 (正式写手)

应助: 5 (幼儿园)
金币: 1202.1
红花: 4
帖子: 591
在线: 55.9小时
虫号: 918027
注册: 2009-12-01
性别: GG
专业: 电磁场与波

[求助] 场域计算中的发散问题

我用FDTD编写了一个往空心金属同轴线填充左手材料的程序，运行中发现数值发散，程序如下，希望高手指点一下！！
clear;clc;
%%%%%%% 设置参数  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c=2.99792458e8;
pi=3.14159267;
epsm_0=8.85e-12;
mu_0=12.56637e-7;
N=80000; %%%%%%%  迭代次数
%%%%%%%%%%  网格剖分 %%%%%%%%%%%%%%%%%%%%%%;
IE=40;
JE=360;
IB=IE+1;
JB=JE+1;
IE1=11;
r3=1.0;%%%%%%%%%%%%%%%%  外圆半径
r2=0.25;%%%%%%%%%%%%%%%%  内圆半径
%%%%% 网格大小 %%%%%%%%
dlr=r3/IE;
dlf=(2*pi)/360;
dlrf=r3*dlf;
%%%%% 源的位置
IS=14;
JS=117;
IS1=14;
JS1=117;
IS2=37;
JS2=217;
%%%%%%%%  采样点的位置
I01=17;
J01=133;
I02=21;
J02=183;
I03=32;
J03=233;
I04=35;
J04=313;
%%%%%%%%% 源的参数设置
s=1/6;
dt=dlr*s/(2*c);
omiga_0=2*pi*28.7e9;
arf=-omiga_0*log(0.01)/(10*pi);
t0=140*dt;
T=20*dt;

%%%% 电磁分量的初始化
hz=zeros(IB,JB);
er=zeros(IE,JB);
ef=zeros(IB,JE);
x=zeros(1,N);%%%%记录傅里叶变换后的值;

%%%%%%%%%%% 电磁参量的设置
Kz=zeros(IB,JB);
Jr=zeros(IE,JB);
Jf=zeros(IB,JE);
omiga=zeros(IB,JB);

%%%%%%  主程序的循环迭代  %%%%%%%%%%%%%%%
for n=1:N;
%%%%%% 加激励源;
  % hz(IS,JS)=hz(IS,JS)+(1-exp(-arf*n*dt))*sin(omiga_0*n*dt);
er(IS1,JS1)=er(IS1,JS1)+exp(-4*pi*((n*dt-t0)/T)^2);
ef(IS1,JS1)=ef(IS1,JS1)+exp(-4*pi*((n*dt-t0)/T)^2);
  % er(IS2,JS2)=er(IS2,JS2)+exp(-4*pi*((n*dt-t0)/T)^2);
  % ef(IS2,JS2)=ef(IS2,JS2)+exp(-4*pi*((n*dt-t0)/T)^2);
  for ii=11:30;
   for jj=1:JE;
      omiga(ii,jj)=sqrt(3)*omiga_0;
      %omiga(ii,jj)=sqrt(2)*2*pi/T;
      Kz(ii,jj)=Kz(ii,jj)+dt*mu_0*omiga(ii,jj)^2*hz(ii,jj);
      Jr(ii,jj)=Jr(ii,jj)+dt*epsm_0*omiga(ii,jj)^2*er(ii,jj);
      Jf(ii,jj)=Jf(ii,jj)+dt*epsm_0*omiga(ii,jj)^2*ef(ii,jj);
   end;
end;


for ii=IE1:IE;
   for jj=1:JE;
      hz(ii,jj)=hz(ii,jj)-dt/((ii-0.5)*mu_0*dlr)*(ii*ef(ii+1,jj)-(ii-1)*ef(ii,jj)-(er(ii,jj+1)-er(ii,jj))/dlf)-dt/mu_0*Kz(ii,jj);%%因为HZ不处在分界面上，所以mu没影响。
   end;
   hz(ii,JB)=hz(ii,1); %%%%%% 周期性边界条件设置
   %hz(15,jj)=hz(16,jj);
   %hz(61,jj)=hz(60,jj);
end;

for ii=IE1:IE;
   for jj=2:JB;
      er(ii,jj)=er(ii,jj)+dt/((ii-0.5)*dlr*epsm_0)*(hz(ii,jj)-hz(ii,jj-1))/dlf-dt/epsm_0*Jr(ii,jj);
   end;
er(ii,1)=er(ii,JB);
end;

for ii=IE1+1:IE;
   for jj=1:JE;
      ef(ii,jj)=ef(ii,jj)-dt/(epsm_0*dlr)*(hz(ii,jj)-hz(ii-1,jj))+dt/epsm_0*Jf(ii,jj);
   end;
end;
HHZ1(n)=hz(I01,J01);
HHZ2(n)=hz(I02,J02);
HHZ3(n)=hz(I03,J03);
HHZ4(n)=hz(I04,J04);
  mesh(hz)
  drawnow;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% 主程序结束!!!!!!!!!

回复此楼