[已解决]二阶mur边界吸收条件为什么发散了

我是fdtd新手 ，编了一个电磁波在矩形波导中传播的fdtd代码。开始的时候截断边界也是PEC，虽然有点发散，但在1000个时间步里还可以忍受。现在我学习二阶mur边界吸收条件，于是在原程序的基础上在截断处加了二阶mur，但是当电磁波一到截断边界就发散相当厉害，值变得很大。如图： >jH%n(TcC  
第165步，电场最大值为1.365 ~f5g\n;  
E Zh.*u@^r  
第264步，电场最大值为3.650 =/zb$d cz  
U2TR>0l  
第363步，电场最大值为285.3 {M&Vh]  
iB0#Z
_  
LzE$z,  
以下是部分二阶mur边界条件代码，没有包含棱边 /ij)[WK@  
do i=x2,x2,1 4te QG  
     do j=1,(ny-1),1 .{ +Obi  
        do k=1,(nz-1),1 66 @#V  
           ey(i,j,k) = -ey_scratch(i-1,j,k) + (LIGHT_SPEED*dt-dx)*(ey(i-1,j,k)+ey_scratch(i,j,k))/(LIGHT_SPEED*dt+dx)+2*dx*(ey_current(i,j,k)+ey_current(i-1,j,k))/(LIGHT_SPEED*dt+dx)+dx*LIGHT_SPEED*dt*LIGHT_SPEED*dt*(ey_current(i,j+1,k)-2*ey_current(i,j,k)+ey_current(i,j-1,k)+ey_current(i-1,j+1,k)-2*ey_current(i-1,j,k)+ey_current(i-1,j-1,k-1))/(2*dy*dy*(LIGHT_SPEED*dt+dx))+dx*LIGHT_SPEED*dt*LIGHT_SPEED*dt*(ey_current(i,j,k+1)-2*ey_current(i,j,k)+ey_current(i,j,k-1)-ey_current(i-1,j,k+1)-2*ey_current(i-1,j,k)+ey_current(i-1,j,k-1))/(2*dz*dz*(LIGHT_SPEED*dt+dx)) 2T >K!jS  
               ez(i,j,k) = -ez_scratch(i-1,j,k) + (LIGHT_SPEED*dt-dx)*(ez(i-1,j,k)+ez_scratch(i,j,k))/(LIGHT_SPEED*dt+dx)+2*dx*(ez_current(i,j,k)+ez_current(i-1,j,k))/(LIGHT_SPEED*dt+dx)+dx*LIGHT_SPEED*dt*LIGHT_SPEED*dt*(ez_current(i,j+1,k)-2*ez_current(i,j,k)+ez_current(i,j-1,k)+ez_current(i-1,j+1,k)-2*ez_current(i-1,j,k)+ez_current(i-1,j-1,k-1))/(2*dy*dy*(LIGHT_SPEED*dt+dx))+dx*LIGHT_SPEED*dt*LIGHT_SPEED*dt*(ez_current(i,j,k+1)-2*ez_current(i,j,k)+ez_current(i,j,k-1)-ez_current(i-1,j,k+1)-2*ez_current(i-1,j,k)+ez_current(i-1,j,k-1))/(2*dz*dz .. G,fh/E+  
\#O}K  

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

你用啥编的,Matlab吗

引用第1楼chenxl0426于2008-10-08 12:07发表的  : 8AQ__&nT  
你用啥编的,Matlab吗

/Os6i&;  
我用fortran编的

Ey中的第二个系数怎么会有+dx，应该是 -dy 吧？

顺便问一句，你是用VisIt绘制的图吗？电场值在什么范围算正常？你的激励源是怎么加入的？

引用第3楼air81666于2008-10-11 12:31发表的  : @?>5~  
Ey中的第二个系数怎么会有+dx，应该是 -dy 吧？

eX1_=?$1P  
jx{ fel  
哪儿？葛德彪老师的书上就是这样的啊。

引用第4楼air81666于2008-10-11 19:49发表的  : |lk:(~DM  
顺便问一句，你是用VisIt绘制的图吗？电场值在什么范围算正常？你的激励源是怎么加入的？

p$ETAvD  
Jw>na _FJ  
不是哈。 ;)SWwhQ  
电场值的正常范围在几千个时间步里面，应该在激励源的电场幅值左右。因为mur边界吸收条件也不能完全吸收波，所以会有小的反射。 gyPwNE  
我的激励源很简单，就是在x=0的界面令ez等于一个正弦函数，ey和ex都为0。

问题解决了！ CJwzjH  
我总结了两点原因： rx@i.+  
1. 我所参考的书上误把一个加号写成了减号。我开始编写的时候的没有注意到这个问题，所以mur边界吸收条件出现错误。 ZG{#CC
=  
2. 之前编写程序的时候，边界棱边没有处理好，进一步造成部分发散。 94]i|2qj*  
谢谢大家的帮助！

有QQ吗？在QQ上聊吧，我的：86693550

引用第7楼shawny于2008-10-12 12:12发表的  : g[~J107%A  
问题解决了！ PlT_]p  
我总结了两点原因： \OWxf[  
1. 我所参考的书上误把一个加号写成了减号。我开始编写的时候的没有注意到这个问题，所以mur边界吸收条件出现错误。 'z)cieFKP  
2. 之前编写程序的时候，边界棱边没有处理好，进一步造成部分发散。 eb6y-TwY  
谢谢大家的帮助！

%K06owV(S)  
gS`Z>+V5!c  
是哪里的符号错了？我参考的是葛德彪老师的那本书