[原创]等离子体数值模拟

引用第9楼gwzhao于2008-08-31 14:38发表的  : kj`h{
Wc[)  
本想上传论文中等离子体FDTD的递推公式的推导过程的资料的，word文档，提示附件类型不匹配，是不是只能上传图片？

'L2[^iF9  
可以压缩成rar文件上载。如果方便的话，是否可以发些资料到我信箱？zqtong@hotmail.com。先谢谢了。

冷等离子体FDTD的Z变换

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% "QCViR  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Bhw|!Y&%  
%%%%%%%%%%%%%%%%%%                1D                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% '@eH)wh@m)  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  FK|q*  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% bGa"r  
%%%%%%%%%初始化 "SN+ ^`  
clear; &p:GB_  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% $nW^Gqwj]1  
%%%%%%%%%%%%%%%%系统参数 GYNLyd)  
g+M& _n  
TimeT=200;%迭代次数 XLsOn(U\&  
KE=2000;%网格树木  iU a `<  
kc=450;%源的位置 ;n=A245W\  
kpstart=500;%等离子体开始位置 {Qj7?}xW  
kpstop=1000;%等离子体终止位置 LkIbvJCV  
DispE=zeros(1,TimeT); pt- 1>Ui  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% P};GcV-  
%%%%%%%%%%%%%物理参数 5 _ a-nWQ  
c0=3e8;%真空中波速 VNWa3`w  
f=150e6; dSIZsapH  
f1=75e6; {-)*.
l=  
lamda=c0/f; M]M(E) *5  
WL=50; RaLc}F)9  
OMIGA=pi/WL; t'L#8MJ  
zdelta=lamda/WL;%网格大小 QLe<).S1B2  
dt=zdelta/(2*c0);%时间间隔 UvM_~qo  
J9yB'yE8  
qn|~z@"  
u0=1e7;%碰撞频率 2q*aq%  
fpe=1e7;%等离子体频率
${)s ~[  
wpe=2*pi*fpe;%等离子体圆频率 z+6%Ya&
ls  
epsz=1/(4*pi*9*10^9); % 真空介电常数 TeWpdUCO  
mu=1/(c0^2*epsz);%磁常数 +~1~f'4J  
ex_low_m1=0; txXt<]N  
ex_low_m2=0; bdkxCt  
ex_high_m1=0; o
Y0`igH  
ex_high_m2=0; APT/z0X>  
gVI2{\a  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% uQtwh08i  
%%%%%%%%%%%%%初始化电磁场 ,]|*~dd>G  
Ex=zeros(1,KE); 'K|tgsvgme  
Ex_Pre=zeros(1,KE); X"3Za[9j  
Hy=zeros(1,KE); n0CS=  
Hy_Pre=zeros(1,KE); >PMLjXK  
Dx=zeros(1,KE); 088"7 s  
Dx_Pre=zeros(1,KE); R3Ka^l8R|  
Sx1=zeros(1,KE); D!CuE7}  
Sx2=zeros(1,KE); TkSeDP  
Sx3=zeros(1,KE); [u}2xsSx  
Sx=zeros(1,KE); 6b9&V`  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% zV"oB9\9O  
%%%%%%%%%%%%%%%%%%开始计算 !`EhVV8u-_  
for T=1:TimeT A J"/T+g_  
    %%%保存前一时间的电磁场 = ~^  
    Ex_Pre=Ex; N"rZK/@}  
    Hy_Pre=Hy; om}jQJ]KH  
    %%%%中间差分计算Dx ~7Jj\@68  
    for i=2:KE [p[C45d=<  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); hz~jyH.h_  
    end {kp^@  
    %%%%%%%%加入源 !E|k#c9
 
    %Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); ,7nA:0P  
       3Y+ bIz!  
    Dx(kc)=sin(OMIGA*T)+sin(0.5*OMIGA*T); f6zS_y9gn  
    if T<WL rYw
UD7ip  
         Dx(kc)=Dx(kc)*0.5*(1-cos(OMIGA*T));%开关函数，升余弦函数 B!zqvShF  
    end H{nYZOf/  
    %%%计算电场Ex zaoC  
    for i=1:kpstart-1 m5rJY/  
        Ex(i)=Dx(i)/epsz; N"q C-h  
    end @%sr#YqY  
    for i=kpstop+1:KE p7kH"j{xD  
        Ex(i)=Dx(i)/epsz; bem-T`>'  
    end *QT|J6ng  
    if T>50 T&PLvyBL  
        gw=0; hX
x
.  
    end J=C63YB  
    i=kpstart; !wZ9P  
         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); A/#Xr  
         Ex(i)=Dx(i)/epsz-Sx1(i); "ZLujpZcG  
         Ex(i)=Ex(i); Ps_q\R  
     e
<[0H 8  
    for i=kpstart+1:kpstop A T'P=)F@  
         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); y:3d`E4Xw  
         Ex(i)=Dx(i)/epsz-Sx1(i); 2>!?EIE7  
         %Ex(i)=Ex(i); *.ZU" 5e  
    end [#-!&>  
    Sx3=Sx2; iu(+ N~  
    Sx2=Sx1; 37QXML  
    Ex(1)=ex_low_m2; .a *^6TC.  
    ex_low_m2=ex_low_m1; .&Ok53]b  
    ex_low_m1=Ex(2); &9F(uk=X  
zr5(nAl  
    Ex(KE)=ex_high_m2; j1
{\nP/  
    ex_high_m2=ex_high_m1; {#'M3z=  
    ex_high_m1=Ex(KE-1); l`A4)8Y@  
    %%%%%%%%%%%%%%%%%%计算磁场 'fkaeFzOl  
    for i=1:KE-1 pRLs*/Bw  
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i)); X ?lF,p  
    end @*rED6zH  
    DispE(T)=Ex(kpstart+50);  {!9i8T  
    %if round(T/100)*100==T 5};$>47m  
      plot(Ex); oK3aW6  
      grid on; r9d dVD  
      pause(0.00001); J\V.J/  
      T 3Ta<7tEM  
      % end qTl/bFD  
     t[-0/-4  
end

对这个一维模型可以做的工作有: DJ9;{,gm  
当模型为 4P>tGO&*x  
真空 == 等离子体 == 真空 时，可以研究等离子体对电磁波的反射，折射和衰减。 4#qjRmt  
     同时研究不同等离子体密度分布对电磁波的影响。 v^lR]9;  
28j=q-9Z  
模型为： QxGQF|  
  真空 == 等离子体 == 金属 时，可以研究等离子体隐身，主要想得到加和不加等离子体，对反射电磁波的吸收，以及何种密度分布对电磁波的吸收最大。 'ZHu=UT7_  
] _/d  
同时，等离子体的FDTD算法有很多种，每一种都对应着不同的递推公式，在一维下，可以每一种都尝试一下，看看精度什么的，变化大不大。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% lq@Vb{Z  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% co3H=#2a  
%%%%%%%%%%%%%%%%%%%%        1D              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4(4JQ(5  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% HHZ!mYr  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  2H<?  
%%%%%%%%%初始化  ZW2#'$b  
clear; Snt=Hil`  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2LYd # !i  
%%%%%%%%%%%%%%%%系统参数 ,l<-*yMD  
TimeT=3000;%迭代次数 7/vr!tbL`p  
KE=2000;%网格树木 p _e-u-  
kc=450;%源的位置 E|9LUPcb  
kpstart=500;%等离子体开始位置 LM'` U-/e$  
kpstop=1000;%等离子体终止位置 G 7)D+],{Y  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% {NK>9phoB  
%%%%%%%%%%%%%物理参数 Ut-6!kAm  
c0=3e8;%真空中波速 ) WIl
j  
zdelta=1e-9;%网格大小 DuvP3(K  
dt=zdelta/(2*c0);%时间间隔 0iXqAa  
f=900e12;%Gause脉冲的载频 =] 5;=>(  
d=3e-15%脉冲底座宽度 ~~1~_0?e  
t0=2.25/f;%脉冲中心时间 ]@J}f}Mjo  
u0=57e12%碰撞频率 ~rCn
ST  
fpe=2000e12;%等离子体频率 tL+8nTL  
wpe=2*pi*fpe;%等离子体圆频率 eq/5$b(  
epsz=1/(4*pi*9*10^9); % 真空介电常数 U>;itHW/  
mu=1/(c0^2*epsz);%磁常数 +>qBK}`  
ex_low_m1=0; _3p:q.  
ex_low_m2=0; z K8#gif@  
ex_high_m1=0; -R'p^cMA  
ex_high_m2=0; ):78GVp  
a0=2*u0/dt+(2/dt)^2; Re1@2a>  
a1=-8/(dt)^2; &SN$D5U'  
a2=-2*u0/dt+(2/dt)^2; dr6 dK
  
b0=wpe^2+2*u0/dt+(2/dt)^2; CHpDzG>]4  
b1=2*wpe^2-8/(dt)^2; (&MSP  
b2=wpe^2-2*u0/dt+(2/dt)^2; |V9%@ Y?  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% &up/`8  
%%%%%%%%%%%%%初始化电磁场 8 (^2  
Ex=zeros(1,KE); k7{fkl9|#  
Ex_Pre=zeros(1,KE); r-YQsu&  
Hy=zeros(1,KE); wI}'wALhA  
Hy_Pre=zeros(1,KE); *A\NjXJl~  
Dx=zeros(1,KE); Dlj=$25  
Dx_Pre=zeros(1,KE); J-PzIFWd  
Sx1=zeros(1,KE); !P Gow  
Sx2=zeros(1,KE); _R(5?rG,  
Sx3=zeros(1,KE); MF\n@lX  
Sx=zeros(1,KE); 'v|2}T*  
Dx=Ex; ($(1KE  
Dx1=Ex; =w A< F  
Dx2=Ex; 3/rvSR!  
Dx3=Ex; IVNNiNN*5  
Ex1=Ex; )Elr8XLw  
Ex2=Ex; el|t6ZT*  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% =cC]8Pz?  
%%%%%%%%%%%%%%%%%%开始计算 >!)VkDAG  
for T=1:TimeT , D1[}Lr=K  
    %%%保存前一时间的电磁场 KR^peWR  
    Ex_Pre=Ex; 1yB;"q&Xd  
    Hy_Pre=Hy; `4EOy:a  
    %%%%中间差分计算Dx T4!]^_t^  
    for i=2:KE x#8=drh.:C  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); LVEVCpp@  
    end MZjiJZaO:L  
    %%%%%%%%加入源 C-Fp)Zs{0  
    Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); hTG d Uw]  
     Ic=V:  
    %%%计算电场Ex _SP u`=~K  
    for i=1:kpstart-1 Xi{(1o4%  
        Ex(i)=Dx(i)/epsz; UiZ61lw  
    end uBE,z>/,;  
    for i=kpstop+1:KE +~fu-%,k  
        Ex(i)=Dx(i)/epsz; e4)gF*  
    end 6$6NVq  
    Dx3=Dx2; .sG,TLE[<  
Dx2=Dx1; ,;d9uG2  
Dx1=Dx; 93WYZ
NpX  
    for i=kpstart:kpstop <dYk|5AdLF  
                   Ex(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex2(i)); wO!hVm,Ta  
    end >Au<y,Tw  
    Ex2=Ex1; NUnP'X=J,  
Ex1=Ex; (x^|  
    Sx3=Sx2; E+1j3Q;  
    Sx2=Sx1; G>RYQ{O  
    Ex(1)=ex_low_m2; CLkVe  
    ex_low_m2=ex_low_m1; \7j)^  
    ex_low_m1=Ex(2); Zj7XmkL  
FoG<$9  
    Ex(KE)=ex_high_m2; PRz oLzr  
    ex_high_m2=ex_high_m1; \~)573'  
    ex_high_m1=Ex(KE-1); =!CuCV7$1O  
    %%%%%%%%%%%%%%%%%%计算磁场 Ez$5wY^J  
    for i=1:KE-1 BkZ%0rw%  
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i)); H0(zE*c~  
    end QNJG}Upl  
    plot(Ex); (&*F`\  
    grid on;
D1Sl+NOV  
    pause(0.01); tJu<#hX  
end

上面则是另一种方法得到的递推公式。 CP#79
=1  
等离子体的众多方法我就不一一例举了。

太牛了,能否上传一些源码,让我们学习啊

在FDTD那边有一些源代码的。但我感觉还是讨论问题收获大一些。

关于代码我基本都可以看懂，但是我有一点疑惑，为什么在等离子体区域，波形的传播没有了？波形只在等立体之前的区域传播。这应该不是吸收的问题吧 ]^iFqQe  

经过几天的调试，发现赵博士的代码中计算等离子体区域电磁场时的代码 k%{ l4  
Ex(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex2(i)); 9uuta4&u
I  
漏写了真空介电常数，应该改为： 5gO/-Zj  
Ex(i)=(1/b0)*(a0*(Dx1(i)/epsz)+a1*(Dx2(i)/epsz)+a2*(Dx3(i)/epsz)-b1*Ex1(i)-b2*Ex2(i)); }BA9Ka#%  
如果没有epsz这一项，Dx1,Dx2,Dx3都很小，以致Ex都很小，近似为0，故得到的结果是，波不能再等离子体区域穿过去。修改后，可以看到，波形穿过等离子区域，波形会被吸收而衰减。 6{5q@9F  
作为一名新手，搞得好辛苦啊……