903阅读
3回复

这是一个关于FDTD的对左手材料的仿真，哪位高手能对其中的场值迭代给小弟说明一下

UID ：58843

注册：2010-05-06
登录：2013-09-08
发帖：10
等级：仿真新人

0楼发表于: 2010-05-11 13:53:23

— 本帖被 destroyer 从【互助速问速答有求速应】移动到本区(2010-05-28) —

   %1D FDTD 双负介质
   % 定义计算参数
    nz=200;    %总计算格数
    TIME=1000;    %总计算时间
    hyi=zeros(1,nz);exi=zeros(1,nz+1);    %网格剖分
    kyi=zeros(1,nz);jxi=zeros(1,nz);     % 等效电磁流
     d=10;       %PML总层数
     m=2;        %PML中参数衰减速率阶数
     C=3.0e8;       %光速
     R0=1.0e-5;     %设定的波吸收率
     pi=3.1415926;
     ee=8.851e-12;
     mu=4*3.1415926e-7;
     f=30*10^9;
     dx=0.0005;dy=0.0005;dz=0.0005;dt=1*1e-12;      %网格剖分与时间剖分
     dez=zeros(1,nz+1);dhz=zeros(1,nz);       %网格中的介质参数


    c1=(m+1)*C*ee/2.0/(d*dz)*log(1/R0);             %最大煤质电参数
    m1=mu*c1/ee;     %最大煤质磁参数

    %DNG频率设定
    wpe=5*1e11;wpm=5*1e11;
    w0=2*pi*f;
    gaie=1e8;gaim=1e8;


    %计算各网格点的介质参数
            for i=1:nz+1
               if  i<=d
                   dez(i)=c1*(d+1-i)*(d+1-i)/d/d;
               elseif i>=nz-d+2
                   dez(i)=c1*(i+d-nz-1)*(i+d-nz-1)/d/d;
               else
                   dez(i)=0;
               end
           end

           for i=1:nz
               if i<=d
                   dhz(i)=m1*(d-i+1.0/2.0)*(d-i+1.0/2.0)/d/d;
               elseif i>=nz-d+1
                   dhz(i)=m1*(i+d-nz-1.0/2.0)*(i+d-nz-1.0/2.0)/d/d;
               else
                   dhz(i)=0;
               end
           end

    %场值迭代
  %   for (t=0;t<=TIME;t++)
    for t=0:TIME

%       for(i=0;i<=nx-1;i++)
           for i=1:nz                                     %更新Hy
                   if i==80
                hyi(i)=(1-dhz(i)*dt/(2*mu))/(1+dhz(i)*dt/(2*mu))*hyi(i)+dt/(mu*(1+dhz(i)*dt/(2*mu)))*(exi(i)-exi(i+1))/dz+sin(2*pi*f*t*dt);
            elseif i>nz/2&i<nz*3/4
                hyi(i)=(1-dhz(i)*dt/(2*mu))/(1+dhz(i)*dt/(2*mu))*hyi(i)+dt/(mu*(1+dhz(i)*dt/(2*mu)))*(exi(i)-exi(i+1)-kyi(i)*dz)/dz;
            else
                hyi(i)=(1-dhz(i)*dt/(2*mu))/(1+dhz(i)*dt/(2*mu))*hyi(i)+dt/(mu*(1+dhz(i)*dt/(2*mu)))*(exi(i)-exi(i+1))/dz;
               end
           end

           for i=nz/2+1:nz*3/4-1                                    %更新ky
                   kyi(i)=(1-0.5*gaim*dt)/(1+0.5*gaim*dt)*kyi(i)+mu*wpm^2*dt/(1+0.5*gaim*dt)*hyi(i);
            end

           for i=2:nz                                 %更新ex
             if i>nz/2&i<nz*3/4
                exi(i)=(1-dez(i)*dt/2/ee)/(1+dez(i)*dt/2/ee)*exi(i)+dt/ee/(1+dez(i)*dt/2/ee)*(hyi(i-1)-hyi(i)-0.5*(jxi(i)+jxi(i-1))*dz)/dz;
            else
                exi(i)=(1-dez(i)*dt/2/ee)/(1+dez(i)*dt/2/ee)*exi(i)+dt/ee/(1+dez(i)*dt/2/ee)*(hyi(i-1)-hyi(i))/dz;
            end
           end
        &nbs ..

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