923阅读
0回复

一个关于FDTD 2D的问题

没有

UID ：58796

注册：2010-05-05
登录：2010-06-05
发帖：24
等级：仿真新人

0楼发表于: 2010-05-19 22:44:04

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

下面是我根据人家改写的关于FDTD 2D正弦波传播的问题，下面是程序：
clear all ;
clc;
IE=100;
JE=100;                                        % 网格数量
L=200;                                         % 迭代次数
Lambda=1.0e-2;                                   % 波长
eps=1.0;                                       % 散射体相对介电常数
% eps_r2=3.5;                                  % 散射体相对介电常数
sigma=3.72*1.0e+7;                                % 散射体电导率
% sigma_e2=0.0;                                % 散射体电导率
WL=40;                                         % 每个波长等分数，剖分越细，
mu_0 = 4.0*pi*1.0e-7;                          % 真空磁导率
eps_0 = 8.854*1.0e-12;                         % 真空介电常数
c_0 = 1/sqrt(mu_0*eps_0);                      % 自由空间光速
Z_0=sqrt(mu_0/eps_0);                          % 自由空间波阻抗
Freq=c_0/Lambda;                               % 频率
Dx=Lambda/WL;                                  % 空间步长
Dt=Dx/(2.0*c_0);                               % 时间步长
ic=3;                                          %平面波源位置

ia=7;
ib=IE-ia+1;
ja=7;
jb=JE-ja+1;

eaf =Dt*sigma/(2.0*eps_0*eps);                 % 电场分量系数
CA=(1.0-eaf)/(1.0+eaf);                        %系数CA
CB=Dt/(eps_0*eps)/(Dx*(1.0+eaf));              %系数CB
haf =Dt*sigma/(2.0*mu_0);                      % 磁场分量系数
CP=(1.0-haf)/(1.0+haf);                        %系数CP
CQ=Dt/mu_0/(Dx*(1.0+haf));                     %系数CQ

for i=1:IE
    for j=1:JE
        ga(i,j)=1/(eps+sigma*Dt/eps_0);         % 求和参量
        gb(i,j)=sigma*Dt/eps_0;                 % 求和参量
    end
end

ez_inc=zeros(1,IE);
hx_inc=zeros(1,JE);

ez_low1=0;
ez_low2=0;
ez_low3=0;
ez_low4=0;

dz=zeros(IE,JE);                     % z方向电荷密度
ez=zeros(IE,JE);                     % z方向电场
iz=zeros(IE,JE);                     % z方向电场求和参量
hx=zeros(IE,JE);                     % x方向磁场
hy=zeros(IE,JE);                     % y方向磁场
ihx=zeros(IE,JE);                    % x方向磁场参量
ihy=zeros(IE,JE);                    % y方向磁场参量

%%%%%%%%%%%%%%%%%%%%%%%%%迭带求解电磁场%%%%%%%%%%%%%%%%%
for T=1:L
    for j=2:JE
        ez_inc(j)=CA*ez_inc(j)+CB*(hx_inc(j-1)-hx_inc(j));
    end

    ez_inc(1)=ez_low2;
    ez_low2=ez_low1;
    ez_low1=ez_inc(2);
    ez_inc(JE)=ez_low3;
    ez_low3=ez_low4;
    ez_low4=ez_inc(JE-1);
%%%%%%%%%%%%%% 电荷密度dz%%%%%%%%%%%%%%%%%%%%%%%%%
    for i=2:IE
        for j=2:JE
            dz(i,j)=CA*dz(i,j)+CB*( hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1));
        end
    end

    source=sin(2*pi*Freq*Dt*T);
    ez_inc(ic)=source;

    for i=ia:ib
        dz(i,ja)=CA*dz(i,ja)+CB*hx_inc(ja-1);
        dz(i,jb)=CA*dz(i,jb)-CB*hx_inc(jb);
    end
%%%%%%%%%%%%%% 电场ez%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:IE
        for j=1:JE
            ez(i,j)=ga(i,j)*( dz(i,j)-iz(i,j) );
            iz(i,j)=iz(i,j)+gb(i,j)*ez(i,j) ;
        end
    end
%%%%%%%%%%%%%%%%%%%%边界电场%%%%%%%%%%%%%%%%%%%%%%%
    for j=1:JE
        ez(1,j)=0;
        ez(IE,j)=0;
    end


    for i=1:JE
        ez(i,1)=0;
        ez(i,JE)=0;
    end

    for j=1:JE-1
        hx_inc(j)=CP*hx_inc(j)+CQ*(ez_inc(j)-ez_inc(j+1));
    end
%%%%%%%%%%%%%%%%%%%%%%%%磁场Hx%%%%%%%%%%%%%%%%%%%%%%
    for i=1:IE
        for j=1:JE-1
            hx(i,j)=CP*hx(i,j)+CQ*(ez(i,j)-ez(i,j+1));
        end
    end

    for i=ia:ib;
        hx(i,ja-1)=CP*hx(i,ja-1)+CQ*ez_inc(ja);
        hx(i,jb)=CP*hx(i,jb)-CQ*ez_inc(jb);
    end ..

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

图片:fig1.jpg
图片:fig2.jpg

没有

发帖回复