1487阅读
4回复

[求助]三维pml matlab 程序，有问题，大侠帮忙啊，急！

UID ：18655

注册：2008-10-07
登录：2009-10-21
发帖：9
等级：旁观者

0楼发表于: 2008-10-21 12:27:17

经过进一步检验，这是修改之后的程序，不过仍存在问题：1、激励源的加法，我的激励源修改后分别加到两个分量上，但距离源点几个网格后仍产生毛刺，说明源的加入还是有问题；2、我的程序是把自由空间和pml层一起做了运算和迭代，其实边角问题已经不存在，不需要考虑了（二维已经验证了），这样编程的好处是缩短程序，全域统一编程。
现在的程序仍有毛刺，请高手帮我改正，万分感激！！！谢谢！！
clear
c=3e8;
muz=4*pi*1e-7;
epsz=1/(c*c*muz);
f=1e+9;
lambda=c/f;
omega=2*pi*f;
ds=lambda/12;
dt=ds/(2*c);
muz=4*pi*1e-7;
R0=1e-5;

%************************************************
ie=50;                                                                    %x向空间步
je=50;                                                                    %y向空间步
ke=50;                                                                    %z向空间步
ib=ie+1;
jb=je+1;
kb=ke+1;
is=25;                                                                     %激励源x向位置
js=25;                                                                     %激励源y向位置
ks=25;                                                                     %激励源z向位置
iebc=8;                                                                    %左、右PML层层数
jebc=8;                                                                    %上、下PML层层数
kebc=8;                                                                    %前、后PML层层数
ibbc=iebc+1;
jbbc=jebc+1;
kbbc=kebc+1;
nmax=100;                                                                  %时间步
dlta=6;                                                                    %观测点距激励源步数

%************************************************
Exy=zeros(ie,jb,kb);
Exz=zeros(ie,jb,kb);
Eyx=zeros(ib,je,kb);
Eyz=zeros(ib,je,kb);
Ezx=zeros(ib,jb,ke);
Ezy=zeros(ib,jb,ke);
Hxy=zeros(ib,je,ke);
Hxz=zeros(ib,je,ke);
Hyx=zeros(ie,jb,ke);
Hyz=zeros(ie,jb,ke);
Hzx=zeros(ie,je,kb);
Hzy=zeros(ie,je,kb);

fid1=fopen('m1.txt','wt');
fid2=fopen('m2.txt','wt');
%************************************************
sigxmax=-log(R0)/(2*iebc*ds)*3*epsz*c;                                     %x向最大电导率
sigymax=-log(R0)/(2*jebc*ds)*3*epsz*c;                                     %y向最大电导率
sigzmax=-log(R0)/(2*kebc*ds)*3*epsz*c;

for i=1:ib                                                                 %x向网格节点处电导率分布
                if      i<=iebc
                    sigx(i)=sigxmax*((iebc+1-i)/iebc)^2;                      %左侧PML层
                elseif i>=ib-iebc+1                                           %右侧PML层
                    sigx(i)=sigxmax*((i-ib+iebc)/iebc)^2;
                else
                    sigx(i)=0;                                              %中间自由空间
                end
      cax(i)=(1-(dt*sigx(i))/(2*epsz))/(1+(dt*sigx(i))/(2*epsz));
      cbx(i)=(dt/epsz/ds)/(1+(dt*sigx(i))/(2*epsz));

end
for i=1:ie                                                                 %x方向网格边中点处电导率分布
                if      i<=iebc
                    sigxh(i)=sigxmax*((iebc+1-i-0.5)/iebc)^2;                 %左侧PML层
                elseif i>=ie-iebc+1                                           %右侧PML层
                    sigxh(i)=sigxmax*((i-ie+iebc-0.5)/iebc)^2;
                else
                    sigxh(i)=0;                                              %中间自由空间
                end
      roux(i)=muz*sigxh(i)/epsz;                                              %磁导率
      cpx(i)=(1-(dt*roux(i))/(2*muz))/(1+(dt*roux(i))/(2*muz));
      cqx(i)=(dt/muz/ds)/(1+(dt*roux(i))/(2*muz));
end

for j=1:jb                                                                 %y方向网格节点处电导率分布
                if      j<=jebc
                    sigy(j)=sigymax*((jebc+1-j)/jebc)^2;                      %下侧PML层
                elseif j>=jb-jebc+1                                           %上侧PML层
                    sigy(j)=sigymax*((j-jb+jebc)/jebc)^2;
                else
                    sigy(j)=0;                                              %中间自由空间
                end
      cay(j)=(1-(dt*sigy(j))/(2*epsz))/(1+(dt*sigy(j))/(2*epsz));
      cby(j)=(dt/epsz/ds)/(1+(dt*sigy(j))/(2*epsz));
end
for j=1:je                                                                 %y方向网格边中点处电导率分布
                if      j<=jebc
                    sigyh(j)=sigymax*((jebc+1-j-0.5)/jebc)^2;                 %下侧PML层
                elseif j>=je-jebc+1                                           %上侧PML层
                    sigyh(j)=sigymax*((j-je+jebc-0.5)/jebc)^2;
                else
                    sigyh(j)=0;                                              %中间自由空间
                end
      rouy(j)=muz*sigyh(j)/epsz;                                              %磁导率
      cpy(j)=(1-(dt*rouy(j))/(2*muz))/(1+(dt*rouy(j))/(2*muz));
      cqy(j)=(dt/muz/ds)/(1+(dt*rouy(j))/(2*muz));
end

for k=1:kb                                                                 %z方向网格节点处电导率分布
                if      k<=kebc
                    sigz(k)=sigzmax*((kebc+1-k)/kebc)^2;                      %前侧PML层
                elseif k>=kb-kebc+1                                           %后侧PML层
                    sigz(k)=sigzmax*((k-kb+kebc)/kebc)^2;
                else
                    sigz(k)=0;                                              %中间自由空间
                end
      caz(k)=(1-(dt*sigz(k))/(2*epsz))/(1+(dt*sigz(k))/(2*epsz));
      cbz(k)=(dt/epsz/ds)/(1+(dt*sigz(k))/(2*epsz));
end
for k=1:ke                                                                 %z方向网格边中点处电导率分布
                if      k<=kebc
                    sigzh(k)=sigzmax*((kebc+1-k-0.5)/kebc)^2;                 %前侧PML层
                elseif k>=ke-kebc+1                                           %后侧PML层
                    sigzh(k)=sigzmax*((k-ke+kebc-0.5)/kebc)^2;
                else
                    sigzh(k)=0;                                              %中间自由空间
                end
      rouz(k)=muz*sigzh(k)/epsz;                                              %磁导率
      cpz(k)=(1-(dt*rouz(k))/(2*muz))/(1+(dt*rouz(k))/(2*muz));
      cqz(k)=(dt/muz/ds)/(1+(dt*rouz(k))/(2*muz));
end

% Ex=Exy+Exz;
% Ey=Eyx+Eyz;
% Ez=Ezy+Ezx;

%************************************************
tic

schedule_box=waitbar(0,'Please wait...'); % display the schedule box

for n=1:nmax

           for i=1:ib
               for j=1:je
                 for k=1:ke
                    if(i==is&&j==js&&k==ks)
                        Hxy(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励
                    else
                        Hxy(i,j,k)=cpy(j).*Hxy(i,j,k)+cqy(j).*(Ezx(i,j,k)+Ezy(i,j,k)-Ezx(i,j+1,k)-Ezy(i,j+1,k));
                    end

                 end
               end
           end


            for i=1:ib
               for j=1:je
                 for k=1:ke
                   if(i==is&&j==js&&k==ks)
                       Hxz(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励
                   else
                       Hxz(i,j,k)=cpz(k).*Hxz(i,j,k)+cqz(k).*(Eyx(i,j,k+1)+Eyz(i,j,k+1)-Eyx(i,j,k)-Eyz(i,j,k));
                   end
%                    Hxz(i,j,k)=cpz(k).*Hxz(i,j,k)+cqz(k).*(Eyx(i,j,k+1)+Eyz(i,j,k+1)-Eyx(i,j,k)-Eyz(i,j,k));
               end
               end
          end
%           Hx=Hxy+Hxz;
          for i=1:ie
              for j=1:jb
                 for k=1:ke
                   Hyz(i,j,k)=cpz(k).*Hyz(i,j,k)+cqz(k).*(Exy(i,j,k)+Exz(i,j,k)-Exy(i,j,k+1)-Exz(i,j,k+1));
                 end
               end
           end
          for i=1:ie
               for j=1:jb
                 for k=1:ke
                   Hyx(i,j,k)=cpx(i).*Hyx(i,j,k)+cqx(i).*(Ezx(i+1,j,k)+Ezy(i+1,j,k)-Ezx(i,j,k)-Ezy(i,j,k));
               end
               end
          end
%           Hy=Hyx+Hyz;
          for i=1:ie
               for j=1:je
                 for k=1:kb
                   Hzx(i,j,k)=cpx(i).*Hzx(i,j,k)+cqx(i).*(Eyx(i,j,k)+Eyz(i,j,k)-Eyx(i+1,j,k)-Eyz(i+1,j,k));
               end
               end
           end
          for i=1:ie
               for j=1:je
                 for k=1:kb
                   Hzy(i,j,k)=cpy(j).*Hzy(i,j,k)+cqy(j).*(Exy(i,j+1,k)+Exz(i,j+1,k)-Exy(i,j,k)-Exz(i,j,k));
               end
               end
          end
% Hz=Hzx+Hzy;

          for i=1:ie
               for j=2:je  % PEC at j=1 & j=jb
                 for k=2:ke  % PEC at k=1 & k=kb  后面相同
                   Exy(i,j,k)=cay(j).*Exy(i,j,k)+cby(j).*(Hzx(i,j,k)+Hzy(i,j,k)-Hzx(i,j-1,k)-Hzy(i,j-1,k));
                 end
               end
          end
          for i=1:ie
               for j=2:je
                 for k=2:ke
                   Exz(i,j,k)=caz(k).*Exz(i,j,k)+cbz(k).*(Hyx(i,j,k-1)-Hyx(i,j,k)+Hyz(i,j,k-1)-Hyz(i,j,k));
                 end
               end
          end

            for i=2:ie
               for j=1:je
                 for k=2:ke
                   Eyz(i,j,k)=caz(k).*Eyz(i,j,k)+cbz(k).*(Hxy(i,j,k)-Hxy(i,j,k-1)+Hxz(i,j,k)-Hxz(i,j,k-1));
                 end
               end
            end

            for i=2:ie                &nb ..

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

离线zks800

UID ：16797

注册：2008-08-15
登录：2017-02-21
发帖：188
等级：积极交流六级

1楼发表于: 2008-10-21 20:59:46

什么问题也不说清楚，只列code，汗。。。。。。。。。

离线gwzhao

方恨少

UID ：17098

注册：2008-08-24
登录：2019-01-09
发帖：1374
等级：荣誉管理员

2楼发表于: 2008-10-21 22:03:04

哈，这个问题我问过了，楼主说，加了源，结果没反应，只有在is,js,ks那一点有电磁场，其它地方没有。

逆流而上

离线cem-uestc

UID ：9061

注册：2008-03-07
登录：2019-01-05
发帖：2575
等级：荣誉管理员

3楼发表于: 2008-10-21 22:11:02

有很多问题，最好按照公式重新写一写

欢迎光临
http://www.mwtee.com/home.php?mod=space&uid=13535

离线zks800

UID ：16797

注册：2008-08-15
登录：2017-02-21
发帖：188
等级：积极交流六级

4楼发表于: 2008-10-22 13:33:43

也不一定是错误的。如果没有记错的话，3D点源的辐射场是按照r的平方衰减的。并且点源的辐射场有解析式的，可以比较一下就知道对错了。

引用第2楼gwzhao于2008-10-21 22:03发表的 : p/B&R@%
哈，这个问题我问过了，楼主说，加了源，结果没反应，只有在is,js,ks那一点有电磁场，其它地方没有。

发帖回复