[RFEDA原创]发一个电磁波在梯度等离子体中传播的程序回馈一下论坛

前些天一直在赶计算电磁学的期中大作业，由于本人将来的研究方向可能是等离子，所以就选了电磁波在等离子体中传播这个题目。由于等离子物理这门课这学期才开始上，原理还不是太理解，随便找了些公式来拼凑，但发现自己做的仿真老出问题。后来参考了gwzhao前辈的一维程序，并经其指点发觉CFL条件在等离子中是不稳定的，顿时豁然开朗，马上把程序写好了。在此感谢一下gwzhao前辈，并把程序贴出来回馈一下论坛。 T)QZ9a  
oSGx7dj+  
clear all; /{|<3CEe  
close all; kVu8/*Q  
l5?fF6#j  
N_x=301; 2?=R_&0Q  
N_y=101; ^;on  
timesteps=1500; 3I|&}+Z6  
er=ones(N_y,N_x); y]1:IJL2;  
un1=zeros(N_y,N_x); }d<x
bL!#  
un0=zeros(N_y,N_x); Jxyeh1zqB  
un_1=zeros(N_y,N_x); E:EXp7  
p8yn? ~]^  
c=3e8;%velocity of light in vacuum aD5jy  
f=8e10;%frequency of the incident wave :wlX`YW+e  
']'H8Y-M  
dx=.5e-3; iF.e
BL%  
dy=1.5e-3; l# -4}95  
dt=1/(1/dx^2+1/dy^2)^.5/c/2;%one have of the CFL coefficent RL0#WBR  
"rnZ<A}
  
x=-dx/2:dx:(N_x-1)*dx-dx/2; 3Zy$NsY3  
y=0:dy:(N_y-1)*dy; G=rgL'{  
x_inc=0.03;%incident plane g[Tl#X7F  
Q@M>DA!d^V  
%O9kq  
e=1.6e-19;%charge of an electron ,mpvGvAI  
me=9.1e-31; V1aWVLltj  
epsilon=1/(4*pi*9*10^9); % permittivity of vacuum `jl 1Q,~2r  
for(i=1:N_x-x_inc/dx-1)%set the gradient of plasma density \tL9`RKpg  
n0(i)=1e19*(1+.08*i);%particle density in the plasma o;.6Y `-fJ  
wpe(i)=(n0(i)*e^2/epsilon/me)^.5;%calculate the plasma frequency from the plasma density z^tws*u],5  
end r3OTU$t?  
' KX'{Gy  
W=0.03e-9;%temperal width of the gaussian pulse < 0S+[7S"  
T0=0.1e-9;%center time of the gaussian pulse :g1C,M~  
S1_X@[t  
a1=c^2*(dt/dx)^2; RPb/U8  
a2=c^2*(dt/dy)^2; q]1HCWde  
gY)NPi}!`  
w=2*pi*f; ql Z()  
k=2*pi*f/c; #Y)Gos  
YX||\  
startnodeleft=round(x_inc/dx);%computate the incident nodes near the incident plane /^#8z(@B  
startnoderight=round(x_inc/dx)+1; XaFu(Xu7  
a(!_3i@  
pstart=0.05;%calculate the shape of the plasma. DA>_9o/l  
%it is a circle center at (N_x,(N_y+1)/2),the radius is iw`,\V&  
% the length of the computational region minus 0.05cm {Ac5(li_  
for(i=1:N_y) -!X,MDO  
    for (j=5/dx+N_x) f7y a0%N  
    Npstart(i)=round(N_x-(((N_x-pstart/dx-1)*dx)^2-((i-51)*1.5e-3)^2)^.5/dx); [&IJy  
    end i5oV,fiZo  
end C~C}b  
Oa[G #  
for (n=1:timesteps) #bFJ6;g=V  
    for(i=2:N_y-1)%calculate the incident wave using TF/SF method tS$^k)ZXip  
        un1(i,startnoderight)=2*(1-(a1+a2))*un0(i,startnoderight)+(a1+a2)*... j%
|#8oV  
            (un0(i,startnoderight-1)+sin(w*(n-1)*dt-k*(x_inc-dx/2))*exp(-... T#L/HD  
            (((n-1)*dt-T0)/W*2)^2)+un0(i,startnoderight+1))-un_1(i,startnoderight); 6\8 lx|w  
        un1(i,startnodeleft)=2*(1-(a1+a2))*un0(i,startnodeleft)+(a1+a2)*... g>])O  
            (un0(i,startnodeleft-1)-sin(w*(n-1)*dt-k*(x_inc+dx/2))*exp(-... /FzO9'kj  
            (((n-1)*dt-T0)/W*2)^2)+un0(i,startnodeleft+1))-un_1(i,startnodeleft); #LP38wE  
    end z(-j%?  
    for(i=2:N_y-1) ROWb:tX}  
        for(j=2:startnodeleft-1)%calculate the node on the left of the plasma nEtG(^N  
        un1(i,j)=1/er(i,j)*(a1*(un0(i,j-1)+un0(i,j+1))+a2*(un0(i-1,j)+... {J`]6ba  
            un0(i+1,j))-2*(a1+a2)*un0(i,j))+2*un0(i,j)-un_1(i,j); wUnz D)  
        end =(!&8U9  
        for(j=startnoderight+1:Npstart(i)-1)%calculate the nodes between the incident wave and the plasma RR*eq.;  
        un1(i,j)=1/er(i,j)*(a1*(un0(i,j-1)+un0(i,j+1))+a2*(un0(i-1,j)+... wHGiN9A+  
            un0(i+1,j))-2*(a1+a2)*un0(i,j))+2*un0(i,j)-un_1(i,j); I-@A{vvPK  
        end fgmu*\x<  
         UIhU[
f]  
        for(j=Npstart(i):N_x-1)%caluclate the nodes in the plasma |AY`OVgcKD  
        un1(i,j)=1/er(i,j)*(a1*(un0(i,j-1)+un0(i,j+1))+a2*(un0(i-1,j)+... 3['aK|qk.  
            un0(i+1,j))-2*(a1+a2)*un0(i,j))+2*un0(i,j)-dt^2*wpe(j-Npstart(i)+1)^2*un0(i,j)-un_1(i,j); bn
so+cA  
        end <"Ox)XG3]W  
    end !/1a
ot^(  
     Eb9 eEa<W  
for (i=2:N_y-1)%cauculate the nodes on the left and right boundary using Mur's second order ABC O<A$,<67  
    un1(i,N_x)=-un_1(i,N_x-1)+(c*dt-dx)/(c*dt+dx)*(un_1(i,N_x)+un1(i,N_x-1))+... :>@6\  
        2*dx/(c*dt+dx)*(un0(i,N_x)+un0(i,N_x-1))+(c*dt)^2*dx/2/(c*dt+dx)/dy^2*... 3<5E254N  
        (un0(i+1,N_x)-2*un0(i,N_x)+un0(i-1,N_x)+un0(i+1,N_x-1)-2*un0(i,N_x-1)+un0(i-1,N_x-1)); 2WKYf0t  
    un1(i,1)=-un_1(i,2)+(c*dt-dx)/(c*dt+dx)*(un_1(i,1)+un1(i,2))+2*dx/(c*dt+dx)*... ccLTA  
        (un0(i,1)+un0(i,2))+(c*dt)^2*dx/2/(c*dt+dx)/dy^2*(un0(i+1,1)-2*un0(i,1)+... 7=$@bHEF#*  
        un0(i-1,1)+un0(i+1,2)-2*un0(i,2)+un0(i-1,2)); ~Ibq,9i  
end )qgcz<p?W  
OxGS{zs  
for (j=2:N_x-1)%cauculate the nodes on the top and bottom boundary using Mur's second order ABC %'}zr>tx:  
    un1(N_y,j)=-un_1(N_y-1,j)+(c*dt-dy)/(c*dt+dy)*(un_1(N_y,j)+un1(N_y-1,j))+... B@v\tpR  
        2*dy/(c*dt+dy)*(un0(N_y,j)+un0(N_y-1,j))+(c*dt)^2*dy/2/(c*dt+dy)/dx^2*... \KBE+yj  
        (un0(N_y,j+1)-2*un0(N_y,j)+un0(N_y,j-1)+un0(N_y-1,j+1)-2*un0(N_y-1,j)+un0(N_y-1,j-1)); `WjRb  
    un1(1,j)=-un_1(2,j)+(c*dt-dy)/(c*dt+dy)*(un_1(1,j)+un1(2,j))+2*dy/(c*dt+dy)*... O8&=qZ6T  
        (un0(1,j)+un0(2,j))+(c*dt)^2*dy/2/(c*dt+dy)/dx^2*(un0(1,j+1)-2*un0(1,j)+...  ,\HZIl[8  
       & .. wi7a_^{  
^A4bsoW  

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

程序使用的是电磁波在等离子体中的波动公式，但是没有将等离子碰撞频率加进去，因为查文献发现碰撞频率是用于计算等离子电介常数的，但公式的推导中并没有出现电介常数项，不知道如何代进去。个人估计这个介电常数的作用大概和等离子振荡频率的作用差不多，可能可以把碰撞项去掉，然后用介电常数来代替它的作用。不知道对不对，如果gwzhao前辈有空的话还望指点一下。 G
i9s*v,s  
!U?Z<z
h  
另外边界条件只是用了Mur's二阶吸收边界，PML不懂怎么用，如果版内高人有空的话指教一下怎么改进好吗

看到我的资料对别人有用，我还是很开心的。呵呵。    HI` q!LPv  
但是好像不对啊。 Xx=jN1=,  
1）“程序使用的是电磁波在等离子体中的波动公式” ？？ 使用了什么公式？？ U]aH4N  
2）如果电磁场的数值模拟的话，应该看到电场和磁场或电流变量啊，但你程序里面只有un0,un1,un_1，这个分别表示什么意思呢？ s$V'|Pt  
3）你加入的电磁波是什么波啊，我怎么有点看不明白呢？

这里没有公式编辑器，大概描述一下公式吧 5kRP Sfh  
梯度²E=时间求导²E+c²Wpe²E -W#-m'Lvu  
在这个公式的推导过程中好像没看到碰撞频率啊和等离子介电常数啊

这里只是计算了电场的波动公式，所以只有电场出现了 `X=2Ff  
un0,un1,un_1 就是迭代的数组，因为用了二阶的分解，所以计算新的时间数组是需要涉及到之前的两个时间数据

加入的电波只是一个高斯脉冲，指向等离子方向，所以另外一个方向就用总场/散射场方法滤掉了

晕，不能这么搞的吧。 ~dIb>[7wy  
电磁场数值模拟一般用来解决数值计算无法解决的问题。如果你直接利用这个公式的话，我觉得就没什么意义了啊。在使用这个公式前，你得了解这个公式是怎么来的，每一步都做了什么样的近似？ kXj%thDx  
电磁场数值模拟应该直接从麦克斯韦方程着手，推导电磁场的递推公式。 M!=WBw8Y]a  
,tEvz  
数值模拟是建立最基本的方程递推公式，设置好边界条件，激励源，剩下的事情就是观察现象，得到结果了。 U`j[Ni}"  
数值计算一般是最基本方程，做大量的近似，在某些特定的条件下，得到解析解。不能把数值计算中的方程随便拿到数值模拟中去使用的，这样就失去了模拟的意义了。 sFEkxZi<  

引用第3楼itnice于2008-10-22 12:00发表的  : -jZP&8dPH  
这里没有公式编辑器，大概描述一下公式吧 ^ExA  
梯度²E=时间求导²E+c²Wpe²E ;x[F4d  
在这个公式的推导过程中好像没看到碰撞频率啊和等离子介电常数啊

不过我这个公式就是从麦克斯韦方程直接推导出来的啊，所有附带参量最后都可以用等离子频率来表示了，中间都没有作过近似. ed4:r/Dpo  
因为大作业的要求是仿真一个可以用波动方程解的问题，所以就把方程组整合成一个公式了，参考的文献是高教出的等离子体物理学，以及F.F. Chen的Introduction to Plasma Physics and Controlled Fusion 。仿真的时候可以模拟出w与wpe的关系对传播造成的影响啊

你把公式推导过程上传上来么？我看在等离子体和方程那边做了什么样的近似

图片:DSCF0435.jpg

推导过程来了，设备比较简陋，希望看得清楚吧