3225阅读
3回复

[讨论]为什么我的PML吸收边界不起作用呀！求高手帮忙

UID ：3862

注册：2007-07-15
登录：2008-06-06
发帖：8
等级：旁观者

0楼发表于: 2007-11-01 20:17:42

哪位高手看看我的二维FDTD-PML程序的吸收边界为什么不起作用呀

程序如下：
clear;
clc;
pi=3.14159265;
eps0=8.854e-12;
mu0=pi*4.0e-7;
c0=3.0e8;
f=2.0*10^9;
bochang=c0/f; %波长
wl=10;
dx=bochang/wl; %网格大小
dy=dx;
dt=dx/(2*c0);
nn=50;
npml=8;
ic=npml+nn/2; %波源位置
jc=npml+nn/2;
Re=0.0001;
p=3;
sigmamax=-log(Re)*(p+1)*sqrt(mu0/eps0)/2/dx/npml;
imax=nn+2*npml;
jmax=nn+2*npml;
hx=zeros(imax,jmax);
hy=zeros(imax,jmax);
ez=zeros(imax,jmax);
ezx=zeros(imax,jmax);
ezy=zeros(imax,jmax);
ezz=zeros(nn,nn);
nt=200; %时间步数
for k=0:nt
for j=2:jmax
for i=2:imax
if i>npml&j>npml&i<=nn+npml&j<=nn+npml
if i==ic&j==jc %或用平面波：if j==jc
ez(i,j)=sin(2*pi*f*k*dt);
else
ez(i,j)=ez(i,j)+(dt/eps0)*((hy(i,j)-hy(i-1,j))/dx-(hx(i,j)-hx(i,j-1))/dy);
end
end
if i<=npml %左边PML
if j<=npml|j>nn+npml %左边两角
di=npml-i+1;
if j<=npml
dj=npml-j+1;
else
dj=j-nn-npml;
end
sigmaX=sigmamax*(di/npml)^p;
sigmaY=sigmamax*(dj/npml)^p;
ezx(i,j)=ezx(i,j)*exp(-sigmaX*dt/eps0)+(1-exp(-sigmaX*dt/eps0))*(hy(i,j)-hy(i-1,j))/sigmaX/dx;
ezy(i,j)=ezy(i,j)*exp(-sigmaY*dt/eps0)-(1-exp(-sigmaY*dt/eps0))*(hy(i,j)-hy(i,j-1))/sigmaY/dy;
ez(i,j)=ezx(i,j)+ezy(i,j);
else %左边非角点
di=npml-i+1;
sigmaX=sigmamax*(di/npml)^p;
sigmaY=0.;
ezx(i,j)=ezx(i,j)*exp(-sigmaX*dt/eps0)+(1-exp(-sigmaX*dt/eps0))*(hy(i,j)-hy(i-1,j))/sigmaX/dx;
ezy(i,j)=ezy(i,j)-dt*(hy(i,j)-hy(i,j-1))/dy/eps0;
ez(i,j)=ezx(i,j)+ezy(i,j);
end
elseif i>nn+npml %右边PML
if j<=npml|j>nn+npml %右边两角
di=i-nn-npml;
if j<=npml
dj=npml-j+1;
else
dj=j-nn-npml;
end
sigmaX=sigmamax*(di/npml)^p;
sigmaY=sigmamax*(dj/npml)^p;
ezx(i,j)=ezx(i,j)*exp(-sigmaX*dt/eps0)+(1-exp(-sigmaX*dt/eps0))*(hy(i,j)-hy(i-1,j))/sigmaX/dx;
ezy(i,j)=ezy(i,j)*exp(-sigmaY*dt/eps0)-(1-exp(-sigmaY*dt/eps0))*(hy(i,j)-hy(i,j-1))/sigmaY/dy;
%ez(i,j)=ezx(i,j)+ezy(i,j);
else %右边非角点
di=i-nn-npml;
sigmaX=sigmamax*(di/npml)^4;
sigmaY=0.;
ezx(i,j)=ezx(i,j)*exp(-sigmaX*dt/eps0)+(1-exp(-sigmaX*dt/eps0))*(hy(i,j)-hy(i-1,j))/sigmaX/dx;
ezy(i,j)=ezy(i,j)-dt*(hy(i,j)-hy(i,j-1))/dy/eps0;
%ez(i,j)=ezx(i,j)+ezy(i,j);
end
end
if j<=npml %上边PML
if i>npml&i<=nn+npml %上边非角点
dj=npml-j+1;
sigmaY=sigmamax*(dj/npml)^p;
sigmaX=0.;
ezx(i,j)=ezx(i,j)+dt*(hy(i,j)-hy(i-1,j))/dx/eps0;
ezy(i,j)=ezy(i,j)*exp(-sigmaY*dt/eps0)-(1-exp(-sigmaY*dt/eps0))*(hy(i,j)-hy(i,j-1))/sigmaY/dy;
%ez(i,j)=ezx(i,j)+ezy(i,j);
end
elseif j>nn+npml %下边PML
if i>=npml&i<nn+npml %下边非角点
dj=j-nn-npml;
sigmaY=sigmamax*(dj/npml)^p;
sigmaX=0.;
ezx(i,j)=ezx(i,j)+dt*(hy(i,j)-hy(i-1,j))/dx/eps0;
ezy(i,j)=ezy(i,j)*exp(-sigmaY*dt/eps0)-(1-exp(-sigmaY*dt/eps0))*(hy(i,j)-hy(i,j-1))/sigmaY/dy;
%ez(i,j)=ezx(i,j)+ezy(i,j);
end
end
end
end
%以下为磁场循环
for j=1:jmax-1
for i=1:imax-1
if i>npml&j>npml&i<=nn+npml&j<=nn+npml
hx(i,j)=hx(i,j)-(dt/mu0)*(ez(i,j+1)-ez(i,j))/dy;
hy(i,j)=hy(i,j)+(dt/mu0)*(ez(i+1,j)-ez(i,j))/dx;
end
if i<=npml %左边PML
if j<=npml|j>nn+npml %左边两角
di=npml-i+0.5;
if j<=npml
dj=npml-j+0.5;
else
dj=j-nn-npml+1.5;
end
sigmaX=sigmamax*(di/npml)^p;
sigmaMx=sigmaX*mu0/eps0;
sigmaY=sigmamax*(dj/npml)^p;
sigmaMy=sigmaY*mu0/eps0;
hx(i,j)=hx(i,j)*exp(-sigmaMy*dt/mu0)-(1-exp(-sigmaMy*dt/mu0))*((ezx(i,j+1)+ezy(i,j+1))-(ezx(i,j)+ezy(i,j)))/sigmaMy/dy;
hy(i,j)=hy(i,j)*exp(-sigmaMx*dt/mu0)+(1-exp(-sigmaMx*dt/mu0))*((ezx(i+1,j)+ezy(i+1,j))-(ezx(i,j)+ezy(i,j)))/sigmaMx/dx;
else %左边非角点
di=npml-i+0.5;
sigmaX=sigmamax*(di/npml)^p;
sigmaMx=sigmaX*mu0/eps0;
sigmaMy=0.;
hx(i,j)=hx(i,j)-dt*(ez(i,j+1)-ez(i,j))/dy/mu0;
hy(i,j)=hy(i,j)*exp(-sigmaMx*dt/mu0)+(1-exp(-sigmaMx*dt/mu0))*((ezx(i+1,j)+ezy(i+1,j))-(ezx(i,j)+ezy(i,j)))/sigmaMx/dx;
end
elseif i>nn+npml %右边PML
if j<=npml|j>nn+npml %右边两角
di=i-nn-npml+1.5;
if j<=npml
dj=npml-j+0.5;
else
dj=j-nn-npml+1.5;
end
sigmaX=sigmamax*(di/npml)^p;
sigmaMx=sigmaX*mu0/eps0;
sigmaY=sigmamax*(dj/npml)^p;
sigmaMy=sigmaY*mu0/eps0;
hx(i,j)=hx(i,j)*exp(-sigmaMy*dt/mu0)-(1-exp(-sigmaMy*dt/mu0))*((ezx(i,j+1)+ezy(i,j+1))-(ezx(i,j)+ezy(i,j)))/sigmaMy/dy;
hy(i,j)=hy(i,j)*exp(-sigmaMx*dt/mu0)+(1-exp(-sigmaMx*dt/mu0))*((ezx(i+1,j)+ezy(i+1,j))-(ezx(i,j)+ezy(i,j)))/sigmaMx/dx;
else %右边非角点
di=i-nn-npml+1.5;
sigmaX=sigmamax*(di/npml)^p;
sigmaMx=sigmaX*mu0/eps0;
sigmaMy=0.;
hx(i,j)=hx(i,j)-dt*(ez(i,j+1)-ez(i,j))/dy/mu0;
hy(i,j)=hy(i,j)*exp(-sigmaMx*dt/mu0)+(1-exp(-sigmaMx*dt/mu0))*((ezx(i+1,j)+ezy(i+1,j))-(ezx(i,j)+ezy(i,j)))/sigmaMx/dx;
end
end
if j<=npml %上边PML
if i>npml&i<=nn+npml %上边非角点
dj=npml-j+0.5;
sigmaY=sigmamax*(dj/npml)^p;
sigmaMy=sigmaY*mu0/eps0;
sigmaMx=0.;
hx(i,j)=hx(i,j)*exp(-sigmaMy*dt/mu0)-(1-exp(-sigmaMy*dt/mu0))*((ezx(i,j+1)+ezy(i,j+1))-(ezx(i, ..

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