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

前些天一直在赶计算电磁学的期中大作业，由于本人将来的研究方向可能是等离子，所以就选了电磁波在等离子体中传播这个题目。由于等离子物理这门课这学期才开始上，原理还不是太理解，随便找了些公式来拼凑，但发现自己做的仿真老出问题。后来参考了gwzhao前辈的一维程序，并经其指点发觉CFL条件在等离子中是不稳定的，顿时豁然开朗，马上把程序写好了。在此感谢一下gwzhao前辈，并把程序贴出来回馈一下论坛。 
ps_q3Cyp  
CB^.N>'  
clear all; 2 9#jKh  
close all; P09f  
x@3" SiC  
N_x=301; Q<>b3X>O  
N_y=101; 5tl($j  
timesteps=1500; Sjp
]TWj  
er=ones(N_y,N_x); F%IvgXt5  
un1=zeros(N_y,N_x); :nS$cC0x*  
un0=zeros(N_y,N_x); n7/>+V+  
un_1=zeros(N_y,N_x); W1_.wN$,5  
#QNa| f#=  
c=3e8;%velocity of light in vacuum $X,dQ]M  
f=8e10;%frequency of the incident wave "1Y'VpKm(~  
mtmTlGp6Lc  
dx=.5e-3; "\+.S]~  
dy=1.5e-3; SS6K7  
dt=1/(1/dx^2+1/dy^2)^.5/c/2;%one have of the CFL coefficent %JL P=(  
ha?M[Vyw4Q  
x=-dx/2:dx:(N_x-1)*dx-dx/2; C`=YGyj=TL  
y=0:dy:(N_y-1)*dy; +_1sFH`  
x_inc=0.03;%incident plane ,lcSJ^yr  
g.pR4Mf=Z  
5x"eM=  
e=1.6e-19;%charge of an electron g#:P cl  
me=9.1e-31; ckYT69U  
epsilon=1/(4*pi*9*10^9); % permittivity of vacuum JbAmud,  
for(i=1:N_x-x_inc/dx-1)%set the gradient of plasma density *WfQi8  
n0(i)=1e19*(1+.08*i);%particle density in the plasma mXs.@u/  
wpe(i)=(n0(i)*e^2/epsilon/me)^.5;%calculate the plasma frequency from the plasma density Rga *68s|&  
end ^[6el_mj  
FfNUFx2N  
W=0.03e-9;%temperal width of the gaussian pulse Y94S!TbB  
T0=0.1e-9;%center time of the gaussian pulse 8j)*T9  
<Hr~|oG  
a1=c^2*(dt/dx)^2; cH6++r  
a2=c^2*(dt/dy)^2; 0X:$ASocU  
4L_)@n}  
w=2*pi*f; &grqRt  
k=2*pi*f/c; %-woaj   
o Vpq*"  
startnodeleft=round(x_inc/dx);%computate the incident nodes near the incident plane h [@}}6  
startnoderight=round(x_inc/dx)+1; r+k g$+%b  
xqb*;TBh*  
pstart=0.05;%calculate the shape of the plasma. saTS8p z  
%it is a circle center at (N_x,(N_y+1)/2),the radius is AsI\#wL)  
% the length of the computational region minus 0.05cm '+\t,>nRkl  
for(i=1:N_y) C55Av%-=  
    for (j=5/dx+N_x) ;0l
Y_ii  
    Npstart(i)=round(N_x-(((N_x-pstart/dx-1)*dx)^2-((i-51)*1.5e-3)^2)^.5/dx); 1Qw_P('}  
    end i1ScXKO  
end wC?$P  
d ehK#8  
for (n=1:timesteps) uUb`Fy9  
    for(i=2:N_y-1)%calculate the incident wave using TF/SF method H?rCIS0  
        un1(i,startnoderight)=2*(1-(a1+a2))*un0(i,startnoderight)+(a1+a2)*... waU2C2!w
 
            (un0(i,startnoderight-1)+sin(w*(n-1)*dt-k*(x_inc-dx/2))*exp(-... Y5c[9\'\  
            (((n-1)*dt-T0)/W*2)^2)+un0(i,startnoderight+1))-un_1(i,startnoderight); g`C\pdX"B  
        un1(i,startnodeleft)=2*(1-(a1+a2))*un0(i,startnodeleft)+(a1+a2)*... k [LV^oEg  
            (un0(i,startnodeleft-1)-sin(w*(n-1)*dt-k*(x_inc+dx/2))*exp(-... @N]]Cf>x  
            (((n-1)*dt-T0)/W*2)^2)+un0(i,startnodeleft+1))-un_1(i,startnodeleft); 6AdC  
    end K#Zv>x!to  
    for(i=2:N_y-1) _2Mpzv  
        for(j=2:startnodeleft-1)%calculate the node on the left of the plasma d(yTz&u)  
        un1(i,j)=1/er(i,j)*(a1*(un0(i,j-1)+un0(i,j+1))+a2*(un0(i-1,j)+... Sc]G7_  
            un0(i+1,j))-2*(a1+a2)*un0(i,j))+2*un0(i,j)-un_1(i,j); Z~].v._YV)  
        end \CX6~  
        for(j=startnoderight+1:Npstart(i)-1)%calculate the nodes between the incident wave and the plasma Y,Lx6kU  
        un1(i,j)=1/er(i,j)*(a1*(un0(i,j-1)+un0(i,j+1))+a2*(un0(i-1,j)+... "W~vSbn7  
            un0(i+1,j))-2*(a1+a2)*un0(i,j))+2*un0(i,j)-un_1(i,j); 'Y.6sB  
        end GVhy }0|  
         zdm2`D;~p  
        for(j=Npstart(i):N_x-1)%caluclate the nodes in the plasma fejC,H4I  
        un1(i,j)=1/er(i,j)*(a1*(un0(i,j-1)+un0(i,j+1))+a2*(un0(i-1,j)+... ]T6pH7~  
            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); ,m8mh)K?0>  
        end E\ls- (,  
    end 7"F29\  
     1m5*MY  
for (i=2:N_y-1)%cauculate the nodes on the left and right boundary using Mur's second order ABC Y?3f Fg  
    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))+... 2;(+]Ad<  
        2*dx/(c*dt+dx)*(un0(i,N_x)+un0(i,N_x-1))+(c*dt)^2*dx/2/(c*dt+dx)/dy^2*... '50}QY_R.  
        (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)); b:cy(6G(  
    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)*... v-BQ>-&
s  
        (un0(i,1)+un0(i,2))+(c*dt)^2*dx/2/(c*dt+dx)/dy^2*(un0(i+1,1)-2*un0(i,1)+... %>$Puy\U  
        un0(i-1,1)+un0(i+1,2)-2*un0(i,2)+un0(i-1,2)); 74 &q2g{  
end `FEa(Q+s  
9F[_xe@  
for (j=2:N_x-1)%cauculate the nodes on the top and bottom boundary using Mur's second order ABC 7]9,J(:Ed  
    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))+... Sy']fGvx  
        2*dy/(c*dt+dy)*(un0(N_y,j)+un0(N_y-1,j))+(c*dt)^2*dy/2/(c*dt+dy)/dx^2*... d8iq9AP\o  
        (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)); 6bPl(.(3  
    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)*... je5[.VTM  
        (un0(1,j)+un0(2,j))+(c*dt)^2*dy/2/(c*dt+dy)/dx^2*(un0(1,j+1)-2*un0(1,j)+... |Sm/s;&c6  
       & .. >7PNl\=gG  
4,>9N9.?9  

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

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

看到我的资料对别人有用，我还是很开心的。呵呵。    dS <*DP  
但是好像不对啊。 <HN{.p{  
1）“程序使用的是电磁波在等离子体中的波动公式” ？？ 使用了什么公式？？ 9iM%kY#)W  
2）如果电磁场的数值模拟的话，应该看到电场和磁场或电流变量啊，但你程序里面只有un0,un1,un_1，这个分别表示什么意思呢？ h~CLJoK<  
3）你加入的电磁波是什么波啊，我怎么有点看不明白呢？

这里没有公式编辑器，大概描述一下公式吧 S5y.H  
梯度²E=时间求导²E+c²Wpe²E Px \cT  
在这个公式的推导过程中好像没看到碰撞频率啊和等离子介电常数啊

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

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

晕，不能这么搞的吧。 E_HB[9  
电磁场数值模拟一般用来解决数值计算无法解决的问题。如果你直接利用这个公式的话，我觉得就没什么意义了啊。在使用这个公式前，你得了解这个公式是怎么来的，每一步都做了什么样的近似？ CI|lJ  
电磁场数值模拟应该直接从麦克斯韦方程着手，推导电磁场的递推公式。 kmuksT\)a  
!,WGd|oJ  
数值模拟是建立最基本的方程递推公式，设置好边界条件，激励源，剩下的事情就是观察现象，得到结果了。 "(koR Q  
数值计算一般是最基本方程，做大量的近似，在某些特定的条件下，得到解析解。不能把数值计算中的方程随便拿到数值模拟中去使用的，这样就失去了模拟的意义了。 fn#8=TIDf  

引用第3楼itnice于2008-10-22 12:00发表的  : BxY t*b%  
这里没有公式编辑器，大概描述一下公式吧 ,w }Po  
梯度²E=时间求导²E+c²Wpe²E TQ Vk;&A  
在这个公式的推导过程中好像没看到碰撞频率啊和等离子介电常数啊

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

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

图片:DSCF0435.jpg

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