786阅读
1回复

[求助]二维pml差分离散格式的几个问题？

UID ：76973

注册：2011-05-08
登录：2015-11-12
发帖：29
等级：仿真新人

0楼发表于: 2011-06-02 16:05:41

这是论坛里下载的程序，我对其中电场的计算步骤不是很理解，%%%%%%%%%%%%%%%%%%%%%%%%计算电场%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% ex%%%%%%%%%%%%%%%%%%
for ki=1:(ke+2*kp)
for kj=(kp+1):(ke+kp+1)
ex(ki,kj)=ex(ki,kj)+0.5*(hz(ki,kj)-hz(ki,kj-1));
end
end                 计算电场时ex(ki,kj)=ex(ki,kj)+0.5*(hz(ki,kj)-hz(ki,kj-1));左边公式得系数怎么是0.5呢？不应该是dt/(dy*e0)?，还有ki的范围不是从kp+1到kp+ke+1?

完整程序如下
clc,clear
%%%% 二维FDTD TE wave pml吸收边界(指数差分)
ke=80;%元包数
kp=8; %pml 层数
ex(ke+2*kp+1,ke+2*kp+1)=0;%实际电场217
ey(ke+2*kp+1,ke+2*kp+1)=0;%实际电场217
hz(ke+2*kp+1,ke+2*kp+1)=0;%实际磁场216，即两端为电场边界
hzx(ke+2*kp+1,ke+2*kp+1)=0;
hzy(ke+2*kp+1,ke+2*kp+1)=0;
kc=(ke+2*kp)/2;%激励源中心 hz(kc,kc)
e0=8.85*10^-12; %   真空介电常数
mu0=4*3.14159*10^-7; % 真空磁导率
c0=1/sqrt(e0*mu0);% 真空中光速
f=1e9; % 频率
lambda=c0/f; %波长
dx=lambda/10;% dx
dy=lambda/10;% dy
dt=(dx/2)/c0; %dt

sigxm=0.33; % 电导率最大值
sigmxm=sigxm*(mu0/e0);%导磁率最大值
nn=3;    %pml层中电导率.导磁率指数
%%%%%%%% 开始主程序 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=0;
%Nsteps=input('Nsteps=')%循环次数
Nsteps=1000;
for n=1:Nsteps
    T=T+1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hz=hzx+hzy;
hz(kc,kc)=sin(2*pi*f*T*dt);
%%%%%%%%%%%%%%%%%%%%%%%%计算电场%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% ex%%%%%%%%%%%%%%%%%%
    for ki=1:(ke+2*kp)
        for kj=(kp+1):(ke+kp+1)
        ex(ki,kj)=ex(ki,kj)+0.5*(hz(ki,kj)-hz(ki,kj-1));
        end
    end
%%%%%% 下边pml层电场ex%%%%%%%%%%%%%%%%%%%%%

for ki=1:(ke+2*kp)
      for kj=2:kp
        sigx=sigxm*((kp+1-kj)/kp)^nn;
        ex(ki,kj)=ex(ki,kj)*exp(-sigx*dt/e0)+(1-exp(-sigx*dt/e0))/(dy*sigx)*sqrt(e0/mu0)*(hz(ki,kj)-hz(ki,kj-1));
        end
end

%%%%%% 上边pml层电场ex%%%%%%%%%%%%%%%%%%%%%
    for ki=1:(ke+2*kp)
        for kj=(ke+kp+2):(ke+2*kp)
        sigx=sigxm*((kj-ke-kp-1)/kp)^nn;
        ex(ki,kj)=ex(ki,kj)*exp(-sigx*dt/e0)+(1-exp(-sigx*dt/e0))/(dy*sigx)*sqrt(e0/mu0)*(hz(ki,kj)-hz(ki,kj-1));
        end
    end

%%%%%%%%%%%%%% ey%%%%%%%%%%%%%%%%%%
    for ki=(kp+1):(ke+kp+1)
        for kj=1:(ke+2*kp)
            ey(ki,kj)=ey(ki,kj)-0.5*(hz(ki,kj)-hz(ki-1,kj));
        end
    end
%%%%%%%%%%%%%%%左边pml层中ey%%%%%%%%%%%%%%%%%%%
for ki=2:kp
      for kj=1:(ke+2*kp)
        sigx=sigxm*((kp+1-ki)/kp)^nn;
        ey(ki,kj)=ey(ki,kj)*exp(-sigx*dt/e0)-(1-exp(-sigx*dt/e0))/(dx*sigx)*sqrt(e0/mu0)*(hz(ki,kj)-hz(ki-1,kj));
        end
    end
     %%%%%% 右边pml层电场ey%%%%%%%%%%%%%%%%%%%%%
    for ki=(ke+kp+2):(ke+2*kp)
        for kj=1:(ke+2*kp)
        sigx=sigxm*((ki-ke-kp-1)/kp)^nn;
        ey(ki,kj)=ey(ki,kj)*exp(-sigx*dt/e0)-(1-exp(-sigx*dt/e0))/(dx*sigx)*sqrt(e0/mu0)*(hz(ki,kj)-hz(ki-1,kj));
        end
    end

%%%%%%%%%%%%%%%电场循环结束%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%磁场循环
%%%%%%%%%%%%%%下边pml中hzy,含sigmy
for ki=1:(ke+2*kp)
    for kj=1:kp
      sigmx=sigmxm*((kp-kj+0.5)/kp)^nn;
      hzy(ki,kj)=hzy(ki,kj)*exp(-sigmx*dt/mu0)+(1-exp(-sigmx*dt/mu0))/(dx*sigmx)*sqrt(mu0/e0)*(ex(ki,kj+1)-ex(ki,kj));
    end
end
%%%%中间hzy(不含sigmy)
for ki=1:(ke+2*kp)
     for kj=(kp+1):(kp+ke)
         hzy(ki,kj)=hzy(ki,kj)+0.5*(ex(ki,kj+1)-ex(ki,kj));
     end
end
%%%%%%%%%%%%%%%pml磁场%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%上边pml中hzy,含sigmx
for ki=1:(ke+2*kp)
     for kj=(kp+ke+1):(ke+2*kp)
      sigmx=sigmxm*((kj-kp-ke-0.5)/kp)^nn;
      hzy(ki,kj)=hzy(ki,kj)*exp(-sigmx*dt/mu0)+(1-exp(-sigmx*dt/mu0))/(dx*sigmx)*sqrt(mu0/e0)*(ex(ki,kj+1)-ex(ki,kj));
     end
end

%%%%%%%%%%%%%%%左边pml中hzx,含sigmx
for ki=1:kp
    for kj=1:(ke+2*kp)
      sigmx=sigmxm*((kp-ki+0.5)/kp)^nn;
      hzx(ki,kj)=hzx(ki,kj)*exp(-sigmx*dt/mu0)-(1-exp(-sigmx*dt/mu0))/(dx*sigmx)*sqrt(mu0/e0)*(ey(ki+1,kj)-ey(ki,kj));
    end
end
%%%%中间hzx(不含sigmx)
for ki=(kp+1):(kp+ke)
     for kj=1:(ke+2*kp)
         hzx(ki,kj)=hzx(ki,kj)-0.5*(ey(ki+1,kj)-ey(ki,kj));
     end
end
%%%%%%%%%%%%%%右边pml中hzx,含sigmx
for ki=(kp+ke+1):(ke+2*kp)
&nbs ..

未注册仅能浏览部分内容，查看全部内容及附件请先登录或注册

离线zwwj_871215

努力，坚强

UID ：75011

注册：2011-04-04
登录：2012-11-12
发帖：118
等级：仿真二级

1楼发表于: 2011-06-03 12:42:52

对，我也不理解，希望有高手指点

发帖回复