[资料共享]发一个我自己在matlab下编写的高斯制fdtd+pml程序代码

        前段时间一直在论坛里面向大家求助有关于高斯制下完全匹配层的差分公式的问题，在大家的热心帮助下，现在我基本上在高斯制下实现了同国际制下同等的吸收效果，这里把我所用的控制方程和程序代码传上来回馈论坛中各位热心坛友的帮助！由于我的编程思路不是很规范，所以在吸收效果上与理论上的完全匹配层吸收效果还有差距，但是已经基本达到了我的预期目的，我在这里也附上几张效果图。目前存在的问题是角点的反射比较严重，也恳请有熟悉这方面的编程高手帮我看看，该怎么改动才能取掉角点的反射，实现更好的吸收。请多多指教！！ O! (85rp/  
    程序代码如下： L!vWRwZwC  
% 二维FDTD  TE波仿真 K0 QH?F  
clear all; r0uJ$/!  
% 定义常数 |0]YA  
pi=3.1415; 6>u
Qt:e  
c=3.0e10;                     %高斯制下光速 D6:DrA:  
f=1.0e15;                     %频率 Q@Dkl F  
lambda=c/f;                   %波长 eI?HwP{m  
nmax=400;                     %时间步数 `
u=oeM:  
del_s=lambda/20;              %每最小波长20个采样点           %ET # z!  
del_t=0.5*del_s/c;            %迭代时间步长 Gw`/.0
 
n=182;                         %真空区域网格数 tWX7dspx/  
np=9;                         %pml层数 Gym#b{#":  
N1=n+2*np;                    %总网格数 9tAE#A  
N=N1+1;                       %采样点数 s9 &)Fv-#V  
M=4;                          %导电率渐变指数 Z9f/-|r5  
sigma_max=(M+1)/1.50/pi/del_s;%最大导电率 b=G4MZQ  
h[y*CzG  
% TE波的分量初始化 <(?' s9  
tic; xD^wTtT  
figure(1); %2YN,a4  
axis([0 N 0 N -0.5 0.5]); m{O Dz:  
Ex=zeros(N1,N);               %x方向为横向，采样点为网格的横向边，故行数+1  +Lhe,  
Ey=zeros(N,N1);               %y方向为纵向，采样点为网格的纵向边，故列数+1 -)>(8f  
Bz=zeros(N1,N1);              %矩阵行为纵向网格数，矩阵列为横向网格数,循环中用j表示行数，i表示列数 W1`Dx(g  
Bzx=zeros(N1,N1); J@gm@ jLc  
Bzy=zeros(N1,N1); UQnBqkE  
Bzxx=zeros(nmax,2); $>![wZ3  
xyvND  
%进入电磁场迭代计算 T+(M8qb  
for tt=1:nmax AHWh}~Yi  
  for i=1:N1 R.O  
        if i>=np+1&&i<=N1-np I}_;A<U  
           di=0; Cy6%S).c
 
        elseif i<=np $$\V2%v  
           di=np-i+0.5; bw020@O*  
        elseif i>=N1-np+1 G ~A$jStm
 
           di=np+i-N1-0.5; []=_<]{   
        end                                 %di是采样点横向距PML内边界的距离 ka8$dfC  
         sigma_mx=sigma_max*(di/np)^M; WV@X@]U  
     for j=1:N1   i)[kubM  
         if j>=np+1&&j<=N1-np nfa_8  
                dj=0; !YY6o V  
         elseif j<=np W7$s5G,  
                dj=np-j+0.5; X*)?LxTj  
         elseif j>=N1-np+1 gY%OhYtF2  
                dj=np+j-N1-0.5; |zh +  
         end                                 %dj是采样点纵向距PML内边界的距离 y]7%$* <  
            sigma_my=sigma_max*(dj/np)^M;                                     X'xUwT|_+  
         T:6K?$y?  
              if sigma_mx+sigma_my==0        %真空区 #)FDl70S8  
                  if j==100&&i==100 \*Ts)EW  
                      t=30;                %可选择高斯脉冲 6UO$z-e  
                      term=(tt-t); %fBP:5%K  
                      pulse=exp(-4*pi*term^2/20^2) ; Enu!u~1]F  
%                       pulse=sin(2*pi*tt/40); %可选正弦时谐源 R qnWtE  
                      c_miu=c*del_t/del_s; !*5_pGe  
                      Eterm1=c_miu*(Ex(i,j+1)-Ex(i,j)); 3Qk/ Ll  
                      Eterm2=c_miu*(Ey(i+1,j)-Ey(i,j)); m^'~&!ba  
                      Bz(i,j)=Bz(i,j)+Eterm1-Eterm2+pulse;%加入脉冲源 o:H'r7N  
                  else z>7=k`x`:  
                      c_miu=c*del_t/del_s; "% SX@  
                      Eterm1=c_miu*(Ex(i,j+1)-Ex(i,j)); gU+ss  
                      Eterm2=c_miu*(Ey(i+1,j)-Ey(i,j)); JBvk)
ogM  
                      Bz(i,j)=Bz(i,j)+Eterm1-Eterm2; 9@Q&B+!  
                  end WX ,p`>n  
              else                            %PML区
B#QL
M^  
                    cpm=(1-2*c*sigma_mx*del_t)/(1+2*c*sigma_mx*del_t); tg3zXJ4k_  
                    cqm=c*del_t/(1+2*c*sigma_mx*del_t)/del_s; jmSt?M0.xV  
                    Bzx(i,j)=cpm*Bzx(i,j)-cqm*(Ey(i+1,j)-Ey(i,j));               3Fgz)*Gu]  
                    cpm=(1-2*c*sigma_my*del_t)/(1+2*c*sigma_my*del_t); ~Po\ En  
                    cqm=c*del_t/(1+2*c*sigma_my*del_t)/del_s; eVrnVPkM  
                    Bzy(i,j)=cpm*Bzy(i,j)+cqm*(Ex(i,j+1)-Ex(i,j)); <r_3obRC  
                    Bz(i,j)=Bzx(i,j)+Bzy(i,j); [A|(A$jl  
              end ?THa5%8f  
     end K[*h+Y
O  
  end `t%|.=R  
  for i=2:N1 k{_1r;  
      if i>=np+1&&i<=N1-np IeH^Wm&^  
         di=0; 40R"^*  
      elseif i<=np ^u3*hl}YKy  
         di=np-i+1; <>VIDE  
      elseif i>=N1-np+1 )dqR<)  
         di=np+i-N1-1; k&yQ98H$K"  
      end                                   %di是采样点横向距PML内边界的距离 >CH  
      sigma_ex=sigma_max*(di/np)^M; p?}f|mQS)  
     for j=1:N1 >B`Cch/'U  
         cam=(1-2*c*sigma_ex*del_t)/(1+2*c*sigma_ex*del_t); d3og?{i<}&  
         cbm=c*del_t/(1+2*c*sigma_ex*del_t)/del_s; (VfwLo>#  
         Ey(i,j)=cam*Ey(i,j)-cbm*(Bz(i,j)-Bz(i-1,j)); |y klT  
     end j{)fC]8H  
  end aXX,Zu^  
  for i=1:N1 7,f:Qi@g  
      for j=2:N1 hH#lTye  
          if j>=np+1&&j<=N1-np CcBQo8!G  
             dj=0; kp'b>&9r  
          elseif j<=np ]F !'M  
             dj=np-j+1; $y8mK|3.3u  
          elseif j>=N1-np+1 J`4Z<b53  
             dj=np+j-N1-1; 3\,MsoAl  
          end                               %dj是采样点纵向距PML内边界的距离 ,tau9>!  
          sigma_ey=sigma_max*(dj/np)^M; ji
qi!*  
          cam=(1-2*c*sigma_ey*del_t)/(1+2*c*sigma_ey*del_t); ZTr:xX{R6  
          cbm=c*del_t/(1+2*c*sigma_ey*del_t)/del_s;< .. dr.**fGYde  
K!9y+%01  

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

角点反射比较严重，这个情况在国际单位制和高斯制中都存在吗？

       嗯，我编的程序国际制下角点的反射也比较严重，可能是我的分区方法不对的缘故，我看看换个思路再尝试一下吧。 \K_!d]I {  
   前段时间我采用的差分格式和现在差不多，但是一直实现不了较好的吸收效果，原因就是出在sigma_max的设置上面，在高斯制下长度单位变为cm，所以同样的空间步长在高斯制下数量级就扩大了100倍，而sigma_max的计算式里面分母含有空间步长，所以如果不加以修正，sigma_max在高斯制下就变成了国际制下的1%，太大的导电率和太小的导电率都不能实现较好的吸收效果，现在我把它修正后就差不多了

呵呵，和我最开始预想的一样，长度单位还是有影响的。 -m x3^  
角点的反射率大概有多少？1%?

那个我还没有进行量级的分析，看图就可以看得出来，怕是不只1%吧

很好，多谢分享！研究一下！

留个脚印，编程的时候参考

图形很好看啊，我的MRTD什么时候能达到这个程度就好了！

坚持+努力 ，加油！ C2T,1=  

引用第7楼ghost32于2008-12-11 21:42发表的  : FeJ5^Gh.  
图形很好看啊，我的MRTD什么时候能达到这个程度就好了！

角点反射比较严重，这个情况在国际单位制和高斯制中都存在吗？ l;fH5z  
您正在看的文章来自RFEDA微波社区--微波仿真论坛 原文地址:http://bbs.rfeda.cn/read.php?tid=17082