1739阅读
13回复

色散性的FDTD仿真求助

方恨少

UID ：17098

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

10楼发表于: 2008-09-02 18:52:03

一维等离子体FDTD的Matlab源代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%                1D                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%初始化
clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%系统参数

TimeT=200;%迭代次数
KE=2000;%网格树木
kc=450;%源的位置
kpstart=500;%等离子体开始位置
kpstop=1000;%等离子体终止位置
DispE=zeros(1,TimeT);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%物理参数
c0=3e8;%真空中波速
f=150e6;
f1=75e6;
lamda=c0/f;
WL=50;
OMIGA=pi/WL;
zdelta=lamda/WL;%网格大小
dt=zdelta/(2*c0);%时间间隔

u0=1e7;%碰撞频率
fpe=1e7;%等离子体频率
wpe=2*pi*fpe;%等离子体圆频率
epsz=1/(4*pi*9*10^9); % 真空介电常数
mu=1/(c0^2*epsz);%磁常数
ex_low_m1=0;
ex_low_m2=0;
ex_high_m1=0;
ex_high_m2=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%初始化电磁场
Ex=zeros(1,KE);
Ex_Pre=zeros(1,KE);
Hy=zeros(1,KE);
Hy_Pre=zeros(1,KE);
Dx=zeros(1,KE);
Dx_Pre=zeros(1,KE);
Sx1=zeros(1,KE);
Sx2=zeros(1,KE);
Sx3=zeros(1,KE);
Sx=zeros(1,KE);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%开始计算
for T=1:TimeT
    %%%保存前一时间的电磁场
    Ex_Pre=Ex;
    Hy_Pre=Hy;
    %%%%中间差分计算Dx
    for i=2:KE
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1));
    end
    %%%%%%%%加入源
    %Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2);

    Dx(kc)=sin(OMIGA*T)+sin(0.5*OMIGA*T);
    if T<WL
         Dx(kc)=Dx(kc)*0.5*(1-cos(OMIGA*T));%开关函数，升余弦函数
    end
    %%%计算电场Ex
    for i=1:kpstart-1
        Ex(i)=Dx(i)/epsz;
    end
    for i=kpstop+1:KE
        Ex(i)=Dx(i)/epsz;
    end
    if T>50
        gw=0;
    end
    i=kpstart;
         Sx1(i)=(1+exp(-1*u0*dt))*Sx2(i)-exp(-1*u0*dt)*Sx3(i)+(wpe^2*dt/(2*u0))*(1-exp(-1*u0*dt))*Ex(i);
         Ex(i)=Dx(i)/epsz-Sx1(i);
         Ex(i)=Ex(i);

    for i=kpstart+1:kpstop
         Sx1(i)=(1+exp(-1*u0*dt))*Sx2(i)-exp(-1*u0*dt)*Sx3(i)+(wpe^2*dt/(u0))*(1-exp(-1*u0*dt))*Ex(i);
         Ex(i)=Dx(i)/epsz-Sx1(i);
         %Ex(i)=Ex(i);
    end
    Sx3=Sx2;
    Sx2=Sx1;
    Ex(1)=ex_low_m2;
    ex_low_m2=ex_low_m1;
    ex_low_m1=Ex(2);

    Ex(KE)=ex_high_m2;
    ex_high_m2=ex_high_m1;
    ex_high_m1=Ex(KE-1);
    %%%%%%%%%%%%%%%%%%计算磁场
    for i=1:KE-1
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i));
    end
    DispE(T)=Ex(kpstart+50);
    %if round(T/100)*100==T
      plot(Ex);
      grid on;
      pause(0.00001);
      T
      % end

end

逆流而上

离线gwzhao

方恨少

UID ：17098

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

11楼发表于: 2008-09-02 18:56:43

这是一维等离子体FDTD源代码，Matlab程序，其中的介电系数为epsonr= epson0 * ( 1 - wpe^2 / (w * ( w - i * up))) 。所以你会发现在求解的时候用了exp( - u0 * dt) ，而不是1 - u0 * dt ，这就是我前面说的，其实1 - u0 * dt 是exp(-u0*dt)的一阶近似。

逆流而上

离线gwzhao

方恨少

UID ：17098

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

12楼发表于: 2008-09-02 19:15:25

一维等离子体FDTD的Matlab源代码（方法二）

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%        1D              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%初始化
clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%系统参数
TimeT=3000;%迭代次数
KE=2000;%网格树木
kc=450;%源的位置
kpstart=500;%等离子体开始位置
kpstop=1000;%等离子体终止位置
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%物理参数
c0=3e8;%真空中波速
zdelta=1e-9;%网格大小
dt=zdelta/(2*c0);%时间间隔
f=900e12;%Gause脉冲的载频
d=3e-15%脉冲底座宽度
t0=2.25/f;%脉冲中心时间
u0=57e12%碰撞频率
fpe=2000e12;%等离子体频率
wpe=2*pi*fpe;%等离子体圆频率
epsz=1/(4*pi*9*10^9); % 真空介电常数
mu=1/(c0^2*epsz);%磁常数
ex_low_m1=0;
ex_low_m2=0;
ex_high_m1=0;
ex_high_m2=0;
a0=2*u0/dt+(2/dt)^2;
a1=-8/(dt)^2;
a2=-2*u0/dt+(2/dt)^2;
b0=wpe^2+2*u0/dt+(2/dt)^2;
b1=2*wpe^2-8/(dt)^2;
b2=wpe^2-2*u0/dt+(2/dt)^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%初始化电磁场
Ex=zeros(1,KE);
Ex_Pre=zeros(1,KE);
Hy=zeros(1,KE);
Hy_Pre=zeros(1,KE);
Dx=zeros(1,KE);
Dx_Pre=zeros(1,KE);
Sx1=zeros(1,KE);
Sx2=zeros(1,KE);
Sx3=zeros(1,KE);
Sx=zeros(1,KE);
Dx=Ex;
Dx1=Ex;
Dx2=Ex;
Dx3=Ex;
Ex1=Ex;
Ex2=Ex;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%开始计算
for T=1:TimeT
    %%%保存前一时间的电磁场
    Ex_Pre=Ex;
    Hy_Pre=Hy;
    %%%%中间差分计算Dx
    for i=2:KE
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1));
    end
    %%%%%%%%加入源
    Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2);

    %%%计算电场Ex
    for i=1:kpstart-1
        Ex(i)=Dx(i)/epsz;
    end
    for i=kpstop+1:KE
        Ex(i)=Dx(i)/epsz;
    end
    Dx3=Dx2;
Dx2=Dx1;
Dx1=Dx;
    for i=kpstart:kpstop
                  Ex(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex2(i));
    end
    Ex2=Ex1;
Ex1=Ex;
    Sx3=Sx2;
    Sx2=Sx1;
    Ex(1)=ex_low_m2;
    ex_low_m2=ex_low_m1;
    ex_low_m1=Ex(2);

    Ex(KE)=ex_high_m2;
    ex_high_m2=ex_high_m1;
    ex_high_m1=Ex(KE-1);
    %%%%%%%%%%%%%%%%%%计算磁场
    for i=1:KE-1
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i));
    end
    plot(Ex);
    grid on;
    pause(0.01);
end

逆流而上

离线随风のkenzo

Blowing In The Wind

UID ：8986

注册：2008-03-06
登录：2018-10-23
发帖：314
等级：荣誉管理员

13楼发表于: 2008-09-02 20:55:16

超过48小时，移动至相关板块继续讨论
此帖锁定
鼠标轻点这里即可参与讨论

发帖回复