官方淘宝店 易迪拓培训 旧站入口
首页 > 微波射频 > 射频工程师交流 > FDTD的PML边界问题救助

FDTD的PML边界问题救助

12-09
编了一个三维的程序,使用PML吸收边界,没有到达边界的时候运行
情况良好,进入PML边界也是有很好的吸收效果。但是一到PML外边界,
就开始反射,这是什么原因造成的?
(我参照《electromagnetic simulation using the fdtd method》
这本书的,书上设置了PML外边界上的E场都为0)
万分感谢!

没有人回答我的问题么?:(

pml的外边界确实会反射的,如果设置成pec
但是pml内部,衰减是指数形式,如果你的pml边界厚度足够大的,反射波应该可以忽略?

大牛啊。。。

杀了我吧

杀了你,那个希腊女人非砍了我不可,呵呵

多谢指导!
我取了八个网格的,40×40×40的,运行50步以后就会反射,而且还
很厉害。但是有人用C语言编的程序可以有很好的吸收效果,我用的是
Matlab。
但是类似的二维的PML边界的吸收效果就相当好,
why?

当然很好的吸收效果啊,不然怎么叫pml

呵呵,你的pml怎么吸收的啊?

o 那就先检查检查你的程序有没有错误咯
某一时刻的场分布画出来看看
在pml边界上是不是完全匹配?
在pml内部是不是成指数衰减的

不瞒您说,我已经耗了两天了,自己检查不出什么东西了。
要不我把程序贴出来,各位大虾帮我看看?
万分感谢!
%fdtd42 3D fdtd  radiation from a dipole antenna  have the PML
clc
clear
IE=40;
JE=40;
KE=40;
ic=IE/2;
jc=JE/2;
kc=KE/2;
ia=7;
ja=7;
ka=7;
ib=IE-ia+1;
jb=JE-ia+1;
kb=KE-ka+1;
ddx=0.01;
dt=ddx/6e8;
epsz=8.8e-12;
muz=4*pi*1e7;
pi=3.14159;
%initial value
dx=zeros(IE,JE,KE);
dy=zeros(IE,JE,KE);
dz=zeros(IE,JE,KE);
ex=zeros(IE,JE,KE);
ey=zeros(IE,JE,KE);
ez=zeros(IE,JE,KE);
hx=zeros(IE,JE,KE);
hy=zeros(IE,JE,KE);
hz=zeros(IE,JE,KE);
ix=zeros(IE,JE,KE);
iy=zeros(IE,JE,KE);
iz=zeros(IE,JE,KE);
gax=ones(IE,JE,KE);
gay=ones(IE,JE,KE);
gaz=ones(IE,JE,KE);
gbx=zeros(IE,JE,KE);
gby=zeros(IE,JE,KE);
gbz=zeros(IE,JE,KE);
%specify the dipole
gaz(ic,jc,11:30)=0;
gaz(ic,jc,kc)=1;
curl_hx=zeros(IE,JE,KE);
curl_hy=zeros(IE,JE,KE);
curl_hz=zeros(IE,JE,KE);
curl_ex=zeros(IE,JE,KE);
curl_ey=zeros(IE,JE,KE);
curl_ez=zeros(IE,JE,KE);
idxl=zeros(IE,JE,KE);
idxh=zeros(IE,JE,KE);
ihxl=zeros(IE,JE,KE);
ihxh=zeros(IE,JE,KE);
idyl=zeros(IE,JE,KE);
idyh=zeros(IE,JE,KE);
ihyl=zeros(IE,JE,KE);
ihyh=zeros(IE,JE,KE);
idzl=zeros(IE,JE,KE);
idzh=zeros(IE,JE,KE);
ihzl=zeros(IE,JE,KE);
ihzh=zeros(IE,JE,KE);
%boundary conditions
gi1=zeros(IE,JE,KE);
fi1=zeros(IE,JE,KE);
gi2=ones(IE,JE,KE);
fi2=ones(IE,JE,KE);
gi3=ones(IE,JE,KE);
fi3=ones(IE,JE,KE);
gj1=zeros(IE,JE,KE);
fj1=zeros(IE,JE,KE);
gj2=ones(IE,JE,KE);
fj2=ones(IE,JE,KE);
gj3=ones(IE,JE,KE);
fj3=ones(IE,JE,KE);
gk1=zeros(IE,JE,KE);
fk1=zeros(IE,JE,KE);
gk2=ones(IE,JE,KE);
fk2=ones(IE,JE,KE);
gk3=ones(IE,JE,KE);
fk3=ones(IE,JE,KE);
npml=7;
for i=1:npml
    xxn1=(npml+1-i)/npml;
    xn1=0.33*xxn1^3;
    fi1(i,1,1)=xn1;
    fi1(IE-i+1,1,1)=xn1;
    gi2(i,1,1)=1/(1+xn1);
    gi2(IE-i+1,1,1)=1/(1+xn1);
    gi3(i,1,1)=(1-xn1)/(1+xn1);
    gi3(IE-i+1,1,1)=(1-xn1)/(1+xn1);
    xxn2=(npml-i+0.5)/npml;
    xn2=0.33*xxn2^3;
    gi1(i,1,1)=xn2;
    gi1(IE-i,1,1)=xn2;
    fi2(i,1,1)=1/(1+xn2);
    fi2(IE-i,1,1)=1/(1+xn2);
    fi3(i,1,1)=(1-xn2)/(1+xn2);
    fi3(IE-i,1,1)=(1-xn2)/(1+xn2);
end
for j=1:npml
    xxn3=(npml+1-j)/npml;
    xn3=0.33*xxn3^3;
    fj1(1,j,1)=xn3;
    fj1(1,JE-j+1,1)=xn3;
    gj2(1,j,1)=1/(1+xn3);
    gj2(1,JE-j+1,1)=1/(1+xn3);
    gj3(1,j,1)=(1-xn3)/(1+xn3);
    gj3(1,JE-j+1,1)=(1-xn3)/(1+xn3);
    xxn4=(npml-j+0.5)/npml;
    xn4=0.33*xxn4^3;
    gj1(1,j,1)=xn4;
    gj1(1,JE-j,1)=xn4;
    fj2(1,j,1)=1/(1+xn4);
    fj2(1,JE-j,1)=1/(1+xn4);
    fj3(1,j,1)=(1-xn4)/(1+xn4);
    fj3(1,JE-j,1)=(1-xn4)/(1+xn4);
end
for k=1:npml
    xxn5=(npml+1-k)/npml;
    xn5=0.33*xxn5^3;
    fk1(1,1,k)=xn5;
    fk1(1,1,KE-k+1)=xn5;
    gk2(1,1,k)=1/(1+xn5);
    gk2(1,1,KE-k+1)=1/(1+xn5);
    gk3(1,1,k)=(1-xn5)/(1+xn5);
    gk3(1,1,KE-k+1)=(1-xn5)/(1+xn5);
    xxn6=(npml-k+0.5)/npml;
    xn6=0.33*xxn6^3;
    gk1(1,1,k)=xn6;
    gk1(1,1,KE-k)=xn6;
    fk2(1,1,k)=1/(1+xn6);
    fk2(1,1,KE-k)=1/(1+xn6);
    fk3(1,1,k)=(1-xn6)/(1+xn6);
    fk3(1,1,KE-k)=(1-xn6)/(1+xn6);
end
t0=20;
spread=6;
for n=1:53
%dx
curl_hx(2:ia,2:JE,2:KE)=hz(2:ia,2:JE,2:KE)-hz(2:ia,1:JE-1,2:KE)-hy(2:ia,2:JE
,2:KE)+hy(2:ia,2:JE,1:KE-1);
idxl(2:ia,2:JE,2:KE)=idxl(2:ia,2:JE,2:KE)+curl_hx(2:ia,2:JE,2:KE);
dx(2:ia,2:JE,2:KE)=gj3(2:ia,2:JE,2:KE)gk3(2:ia,2:JE,2:KE)dx(2:ia,2:JE,2:
KE)+gj2(2:ia,2:JE,2:KE)gk2(2:ia,2:JE,2:KE)*0.5(curl_hx(2:ia,2:JE,2:KE)+g
i1(2:ia,2:JE,2:KE)idxl(2:ia,2:JE,2:KE));
curl_hx(ia+1:ib,2:JE,2:KE)=hz(ia+1:ib,2:JE,2:KE)-hz(ia+1:ib,1:JE-1,2:KE)-hy(
ia+1:ib,2:JE,2:KE)+hy(ia+1:ib,2:JE,1:KE-1);
dx(ia+1:ib,2:JE,2:KE)=gj3(ia+1:ib,2:JE,2:KE)gk3(ia+1:ib,2:JE,2:KE)dx(ia+
1:ib,2:JE,2:KE)+gj2(ia+1:ib,2:JE,2:KE)gk2(ia+1:ib,2:JE,2:KE)*0.5curl_hx(
ia+1:ib,2:JE,2:KE);
for i=ib+1:IE
    ixd=i-ib;
    curl_hx(i,2:JE,2:KE)=hz(i,2:JE,2:KE)-hz(i,1:JE-1,2:KE)-hy(i,2:JE,2:KE)+h
y(i,2:JE,1:KE-1);
    idxh(ixd,2:JE,2:KE)=idxh(ixd,2:JE,2:KE)+curl_hx(i,2:JE,2:KE);
dx(i,2:JE,2:KE)=gj3(i,2:JE,2:KE)gk3(i,2:JE,2:KE)dx(i,2:JE,2:KE)+gj2(i,2:
JE,2:KE)gk2(i,2:JE,2:KE)*0.5(curl_hx(i,2:JE,2:KE)+gi1(i,2:JE,2:KE)idxh
(ixd,2:JE,2:KE));
end
%dy
curl_hy(2:IE,2:ja,2:KE)=hx(2:IE,2:ja,2:KE)-hx(2:IE,2:ja,1:KE-1)-hz(2:IE,2:ja
,2:KE)+hz(1:IE-1,2:ja,2:KE);
idyl(2:IE,2:ja,2:KE)=idyl(2:IE,2:ja,2:KE)+curl_hy(2:IE,2:ja,2:KE);
dy(2:IE,2:ja,2:KE)=gi3(2:IE,2:ja,2:KE)gk3(2:IE,2:ja,2:KE)dy(2:IE,2:ja,2:
KE)+gi2(2:IE,2:ja,2:KE)gk2(2:IE,2:ja,2:KE)*0.5(curl_hy(2:IE,2:ja,2:KE)+g
j1(2:IE,2:ja,2:KE)idyl(2:IE,2:ja,2:KE));
curl_hy(2:IE,ja+1:jb,2:KE)=hx(2:IE,ja+1:jb,2:KE)-hx(2:IE,ja+1:jb,1:KE-1)-hz(
2:IE,ja+1:jb,2:KE)+hz(1:IE-1,ja+1:jb,2:KE);
dy(2:IE,ja+1:jb,2:KE)=gi3(2:IE,ja+1:jb,2:KE)gk3(2:IE,ja+1:jb,2:KE)dy(2:I
E,ja+1:jb,2:KE)+gi2(2:IE,ja+1:jb,2:KE)gk2(2:IE,ja+1:jb,2:KE)*0.5curl_hy(
2:IE,ja+1:jb,2:KE);
for j=jb+1:JE
    jyd=j-jb;
    curl_hy(2:IE,j,2:KE)=hx(2:IE,j,2:KE)-hx(2:IE,j,1:KE-1)-hz(2:IE,j,2:KE)+h
z(1:IE-1,j,2:KE);
    idyh(2:IE,jyd,2:KE)=idyh(2:IE,jyd,2:KE)+curl_hy(2:IE,j,2:KE);
dy(2:IE,j,2:KE)=gi3(2:IE,j,2:KE)gk3(2:IE,j,2:KE)dy(2:IE,j,2:KE)+gi2(2:IE
,j,2:KE)gk2(2:IE,j,2:KE)*0.5(curl_hy(2:IE,j,2:KE)+gj1(2:IE,j,2:KE)idyh
(2:IE,jyd,2:KE));
end
%dz
curl_hz(2:IE,2:JE,2:ka)=hy(2:IE,2:JE,2:ka)-hy(1:IE-1,2:JE,2:ka)-hx(2:IE,2:JE
,2:ka)+hx(2:IE,1:JE-1,2:ka);
idzl(2:IE,2:JE,2:ka)=idzl(2:IE,2:JE,2:ka)+curl_hz(2:IE,2:JE,2:ka);
dz(2:IE,2:JE,2:ka)=gi3(2:IE,2:JE,2:ka)gj3(2:IE,2:JE,2:ka)dz(2:IE,2:JE,2:
ka)+gi2(2:IE,2:JE,2:ka)gj2(2:IE,2:JE,2:ka)*0.5(curl_hz(2:IE,2:JE,2:ka)+g
k1(2:IE,2:JE,2:ka)idzl(2:IE,2:JE,2:ka));
curl_hz(2:IE,2:JE,ka+1:kb)=hy(2:IE,2:JE,ka+1:kb)-hy(1:IE-1,2:JE,ka+1:kb)-hx(
2:IE,2:JE,ka+1:kb)+hx(2:IE,1:JE-1,ka+1:kb);
dz(2:IE,2:JE,ka+1:kb)=gi3(2:IE,2:JE,ka+1:kb)gj3(2:IE,2:JE,ka+1:kb)dz(2:I
E,2:JE,ka+1:kb)+gi2(2:IE,2:JE,ka+1:kb)gj2(2:IE,2:JE,ka+1:kb)*0.5curl_hz(
2:IE,2:JE,ka+1:kb);
for k=kb+1:KE
    kzd=k-kb;
    curl_hz(2:IE,2:JE,k)=hy(2:IE,2:JE,k)-hy(1:IE-1,2:JE,k)-hx(2:IE,2:JE,k)+h
x(2:IE,1:JE-1,k);
    idzh(2:IE,2:JE,kzd)=idzh(2:IE,2:JE,kzd)+curl_hz(2:IE,2:JE,k);
dz(2:IE,2:JE,k)=gi3(2:IE,2:JE,k)gj3(2:IE,2:JE,k)dz(2:IE,2:JE,k)+gi2(2:IE
,2:JE,k)gj2(2:IE,2:JE,k)*0.5(curl_hz(2:IE,2:JE,k)+gk1(2:IE,2:JE,k)idzh
(2:IE,2:JE,kzd));
end
%source
pulse(n)=exp(-0.5*((t0-n)/spread).^2);
dz(ic,jc,kc)=pulse(n);
%caculate the E from D field
ex(2:IE-1,2:JE-1,2:KE-1)=gax(2:IE-1,2:JE-1,2:KE-1)(dx(2:IE-1,2:JE-1,2:KE-1
)-ix(2:IE-1,2:JE-1,2:KE-1));
ix(2:IE-1,2:JE-1,2:KE-1)=ix(2:IE-1,2:JE-1,2:KE-1)+gbx(2:IE-1,2:JE-1,2:KE-1).
*ex(2:IE-1,2:JE-1,2:KE-1);
ey(2:IE-1,2:JE-1,2:KE-1)=gay(2:IE-1,2:JE-1,2:KE-1)(dy(2:IE-1,2:JE-1,2:KE-1
)-iy(2:IE-1,2:JE-1,2:KE-1));
iy(2:IE-1,2:JE-1,2:KE-1)=iy(2:IE-1,2:JE-1,2:KE-1)+gby(2:IE-1,2:JE-1,2:KE-1).
*ey(2:IE-1,2:JE-1,2:KE-1);
ez(2:IE-1,2:JE-1,2:KE-1)=gaz(2:IE-1,2:JE-1,2:KE-1)(dz(2:IE-1,2:JE-1,2:KE-1
)-iz(2:IE-1,2:JE-1,2:KE-1));
iz(2:IE-1,2:JE-1,2:KE-1)=iz(2:IE-1,2:JE-1,2:KE-1)+gbz(2:IE-1,2:JE-1,2:KE-1).
*ez(2:IE-1,2:JE-1,2:KE-1);
%part of the PML is E=0 at the edges
ex(1:IE,1,1:KE)=0;
ex(1:IE,JE,1:KE)=0;
ex(1:IE,1:JE,1)=0;
ex(1:IE,1:JE,KE)=0;
ey(1,1:JE,1:KE)=0;
ey(IE,1:JE,1:KE)=0;
ey(1:IE,1:JE,1)=0;
ey(1:IE,1:JE,KE)=0;
ez(1,1:JE,1:KE)=0;
ez(IE,1:JE,1:KE)=0;
ez(1:IE,1,1:KE)=0;
ez(1:IE,JE,1:KE)=0;
%caculate the hx field
curl_ex(1:ia,1:JE-1,1:KE-1)=ey(1:ia,1:JE-1,2:KE)-ey(1:ia,1:JE-1,1:KE-1)-ez(1
:ia,2:JE,1:KE-1)+ez(1:ia,1:JE-1,1:KE-1);
ihxl(1:ia,1:JE-1,1:KE-1)=ihxl(1:ia,1:JE-1,1:KE-1)+curl_ex(1:ia,1:JE-1,1:KE-1
);
hx(1:ia,1:JE-1,1:KE-1)=fj3(1:ia,1:JE-1,1:KE-1)fk3(1:ia,1:JE-1,1:KE-1)hx(
1:ia,1:JE-1,1:KE-1)+fj2(1:ia,1:JE-1,1:KE-1)fk2(1:ia,1:JE-1,1:KE-1)*0.5(c
url_ex(1:ia,1:JE-1,1:KE-1)+fi1(1:ia,1:JE-1,1:KE-1)ihxl(1:ia,1:JE-1,1:KE-1)
);
curl_ex(ia+1:ib,1:JE-1,1:KE-1)=ey(ia+1:ib,1:JE-1,2:KE)-ey(ia+1:ib,1:JE-1,1:K
E-1)-ez(ia+1:ib,2:JE,1:KE-1)+ez(ia+1:ib,1:JE-1,1:KE-1);
hx(ia+1:ib,1:JE-1,1:KE-1)=fj3(ia+1:ib,1:JE-1,1:KE-1)fk3(ia+1:ib,1:JE-1,1:K
E-1)hx(ia+1:ib,1:JE-1,1:KE-1)+fj2(ia+1:ib,1:JE-1,1:KE-1)fk2(ia+1:ib,1:JE
-1,1:KE-1)*0.5curl_ex(ia+1:ib,1:JE-1,1:KE-1);
for i=ib+1:IE
    ixh=i-ib;
    curl_ex(i,1:JE-1,1:KE-1)=ey(i,1:JE-1,2:KE)-ey(i,1:JE-1,1:KE-1)-ez(i,2:JE
,1:KE-1)+ez(i,1:JE-1,1:KE-1);
    ihxh(ixh,1:JE-1,1:KE-1)=ihxh(ixh,1:JE-1,1:KE-1)+curl_ex(i,1:JE-1,1:KE-1)
;
    hx(i,1:JE-1,1:KE-1)=fj3(i,1:JE-1,1:KE-1)fk3(i,1:JE-1,1:KE-1)hx(i,1:J
E-1,1:KE-1)+fj2(i,1:JE-1,1:KE-1)fk2(i,1:JE-1,1:KE-1)*0.5(curl_ex(i,1:JE-
1,1:KE-1)+fi1(i,1:JE-1,1:KE-1)ihxh(ixh,1:JE-1,1:KE-1));
end
%calculate the hy field
curl_ey(1:IE-1,1:ja,1:KE-1)=ez(2:IE,1:ja,1:KE-1)-ez(1:IE-1,1:ja,1:KE-1)-ex(1
:IE-1,1:ja,2:KE)+ex(1:IE-1,1:ja,1:KE-1);
ihyl(1:IE-1,1:ja,1:KE-1)=ihyl(1:IE-1,1:ja,1:KE-1)+curl_ey(1:IE-1,1:ja,1:KE-1
);
hy(1:IE-1,1:ja,1:KE-1)=fi3(1:IE-1,1:ja,1:KE-1)fk3(1:IE-1,1:ja,1:KE-1)hy(
1:IE-1,1:ja,1:KE-1)+fi2(1:IE-1,1:ja,1:KE-1)fk2(1:IE-1,1:ja,1:KE-1)*0.5(c
url_ey(1:IE-1,1:ja,1:KE-1)+fj1(1:IE-1,1:ja,1:KE-1)ihyl(1:IE-1,1:ja,1:KE-1)
);
curl_ey(1:IE-1,ja+1:jb,1:KE-1)=ez(2:IE,ja+1:jb,1:KE-1)-ez(1:IE-1,ja+1:jb,1:K
E-1)-ex(1:IE-1,ja+1:jb,2:KE)+ex(1:IE-1,ja+1:jb,1:KE-1);
hy(1:IE-1,ja+1:jb,1:KE-1)=fi3(1:IE-1,ja+1:jb,1:KE-1)fk3(1:IE-1,ja+1:jb,1:K
E-1)hy(1:IE-1,ja+1:jb,1:KE-1)+fi2(1:IE-1,ja+1:jb,1:KE-1)fk2(1:IE-1,ja+1:
jb,1:KE-1)*0.5curl_ey(1:IE-1,ja+1:jb,1:KE-1);
for j=jb+1:JE
    jyh=j-jb;
    curl_ey(1:IE-1,j,1:KE-1)=ez(2:IE,j,1:KE-1)-ez(1:IE-1,j,1:KE-1)-ex(1:IE-1
,j,2:KE)+ex(1:IE-1,j,1:KE-1);
    ihyh(1:IE-1,jyh,1:KE-1)=ihyh(1:IE-1,jyh,1:KE-1)+curl_ey(1:IE-1,j,1:KE-1)
;
    hy(1:IE-1,j,1:KE-1)=fi3(1:IE-1,j,1:KE-1)fk3(1:IE-1,j,1:KE-1)hy(1:IE-
1,j,1:KE-1)+fi2(1:IE-1,j,1:KE-1)fk2(1:IE-1,j,1:KE-1)*0.5(curl_ey(1:IE-1,
j,1:KE-1)+fj1(1:IE-1,j,1:KE-1)ihyh(1:IE-1,jyh,1:KE-1));
end
%caculate the hz field
curl_ez(1:IE-1,1:JE-1,1:ka)=ex(1:IE-1,2:JE,1:ka)-ex(1:IE-1,1:JE-1,1:ka)-ey(2
:IE,1:JE-1,1:ka)+ey(1:IE-1,1:JE-1,1:ka);
ihzl(1:IE-1,1:JE-1,1:ka)=ihzl(1:IE-1,1:JE-1,1:ka)+curl_ez(1:IE-1,1:JE-1,1:ka
);
hz(1:IE-1,1:JE-1,1:ka)=fi3(1:IE-1,1:JE-1,1:ka)fj3(1:IE-1,1:JE-1,1:ka)hz(
1:IE-1,1:JE-1,1:ka)+fi2(1:IE-1,1:JE-1,1:ka)fj2(1:IE-1,1:JE-1,1:ka)*0.5(c
url_ez(1:IE-1,1:JE-1,1:ka)+fk1(1:IE-1,1:JE-1,1:ka)ihzl(1:IE-1,1:JE-1,1:ka)
);
curl_ez(1:IE-1,1:JE-1,ka+1:kb)=ex(1:IE-1,2:JE,ka+1:kb)-ex(1:IE-1,1:JE-1,ka+1
:kb)-ey(2:IE,1:JE-1,ka+1:kb)+ey(1:IE-1,1:JE-1,ka+1:kb);
hz(1:IE-1,1:JE-1,ka+1:kb)=fi3(1:IE-1,1:JE-1,ka+1:kb)fj3(1:IE-1,1:JE-1,ka+1
:kb)hz(1:IE-1,1:JE-1,ka+1:kb)+fi2(1:IE-1,1:JE-1,ka+1:kb)fj2(1:IE-1,1:JE-
1,ka+1:kb)*0.5curl_ez(1:IE-1,1:JE-1,ka+1:kb);
for k=kb+1:KE
    kzh=k-kb;
    curl_ez(1:IE-1,1:JE-1,k)=ex(1:IE-1,2:JE,k)-ex(1:IE-1,1:JE-1,k)-ey(2:IE,1
:JE-1,k)+ey(1:IE-1,1:JE-1,k);
    ihzh(1:IE-1,1:JE-1,kzh)=ihzh(1:IE-1,1:JE-1,kzh)+curl_ez(1:IE-1,1:JE-1,k)
;
    hz(1:IE-1,1:JE-1,k)=fi3(1:IE-1,1:JE-1,k)fj3(1:IE-1,1:JE-1,k)hz(1:IE-
1,1:JE-1,k)+fi2(1:IE-1,1:JE-1,k)fj2(1:IE-1,1:JE-1,k)*0.5(curl_ez(1:IE-1,
1:JE-1,k)+fk1(1:IE-1,1:JE-1,k)ihzh(1:IE-1,1:JE-1,kzh));
end
end
E=(ez(1:IE,1:JE,kc));
contour(E)
hold on

才两天就让别人给自己调程序
太f了吧

Top