【求助】看看我编的二维FDTD程序
以下是我编写的二维TM波的FDTD 的FORTRAN程序是自由空间中的电磁波的传播模拟，没有其他的限制条件，只是想实现模拟电磁波在二维真空中的传播情况，但这个程序的结果却是不稳定的，就是电场随时间无限增大，我看了下这个程序是符合时间稳定条件的，现在怎么也找不到毛病出在哪里了，请高手帮忙看看问题出在哪里：
program fdtd2D
real,dimension(:,:),allocatable::hx,hy,ez
real eps,mu,dt,dx,dy,c0,pi,f
pi=3.14159265
c0=3.0e8 !光速
f=3.e12!激励源频率
dx=0.05
dy=dx
dt=dx/(2*c0)
allocate (hx(1:nn,1:nn))
allocate (hy(1:nn,1:nn))
allocate (ez(1:nn,1:nn))
hx=0.
hy=0.
ez=0.
eps=8.854e-12 !真空
mu=4.0e-7 !真空
nt=100 !时间步
nn=100 !网格数
open(1,file='D:\MSDEV\ex.txt')
!open(2,file='D:\MSDEV\cell.txt')
do n=1,nt!外围循环是时间
do j=2,nn
do i=2,nn
if (i.eq.2.and.j.eq.2) then
ez(i,j)=sin(2*pi*f*n*dt) !激励源 !
else
ez(i,j)=ez(i,j)+(dt/eps)*((hy(i,j)-hy(i-1,j))/dx-(hx(i,j)-hx(i,j-1))/dy)
end if
end do
end do
do j=1,nn-1
do i=1,nn-1
hx(i,j)=hx(i,j)-(dt/mu)*(ez(i,j+1)-ez(i,j))/dy
hy(i,j)=hy(i,j)+(dt/mu)*(ez(i+1,j)-ez(i,j))/dx
end do
end do
end do
write(1,'(f20.10)') (ez(i,15),i=1,nn) !输出电场
!write(2,'(i5)') (i,i=1,nn)
end
[ 本帖最后由 lanzhiyu520 于 2007-10-18 20:47 编辑 ] 我把数组的初值重新定义了一下，就是赋初值0
do i=1,nn
do j=1,nn
hx(i,j)=0.
end do
end do
do i=1,nn
do j=1,nn
hy(i,j)=0.
end do
end do
do i=1,nn
do j=1,nn
ez(i,j)=0.
end do
end do
替换掉下面的代码：
hx=0.
hy=0.
ez=0.
结果还是会增大，只是增大的速度慢了 hx(i,j)=hx(i,j)-(dt/mu)*(ez(i,j+1)-ez(i,j))/dy
好象应为
hx(i,j)=hx(i,j)+(dt/mu)*(ez(i,j+1)-ez(i,j))/dy
-------------------------------------------------------
do j=1,nn-1
do i=1,nn-1
hx(i,j)=hx(i,j)-(dt/mu)*(ez(i,j+1)-ez(i,j))/dy
hy(i,j)=hy(i,j)+(dt/mu)*(ez(i+1,j)-ez(i,j))/dx
end do
end do
end do 回3楼的好心人：书上的是“-”号，不过我按你的建议试了下，结果还是会增大
公式应该没有错呀 现在我把源改在中心了，网格也改成300nm了，并且把电场和磁场放在纵坐标内部循环了，结果还是变大，就是不收敛！！
program fdtd2D
real*8,dimension(:,:),allocatable::hx,hy,ez
real eps0,mu0,dt,dx,dy,c0,pi,f
integer pmin,pmax
pi=3.14159265
c0=3.0e8!光速
f=3.e12 !源频率
dx=3.e-7
dy=dx
dt=dx/(2*c0)
eps0=8.854e-12
mu0=4.0e-7
nn=100 !网格数
pmin=0
pmax=nn
open(1,file='D:\MSDEV\ex.txt')
allocate (hx(pmin:pmax,pmin:pmax))
allocate (hy(pmin:pmax,pmin:pmax))
allocate (ez(pmin:pmax,pmin:pmax))
do i=pmin,pmax
do j=pmin,pmax
hx(i,j)=0.
end do
end do
do i=pmin,pmax
do j=pmin,pmax
hy(i,j)=0.
end do
end do
do i=pmin,pmax
do j=pmin,pmax
ez(i,j)=0.
end do
end do
nt=20 !时间步数
do n=1,nt
do j=pmin+1,pmax
if(j>=pmin+1) then
do i=pmin+1,pmax
if (i.eq.50.and.j.eq.50) then
ez(i,j)=sin(2*pi*f*n*dt) !激励源 !
else
ez(i,j)=ez(i,j)+(dt/eps0)*((hy(i,j)-hy(i-1,j))/dx-(hx(i,j)-hx(i,j-1))/dy)
end if
end do
end if
if(j<=pmax-1) then
do i=pmin,pmax-1
hx(i,j)=hx(i,j)-(dt/mu0)*(ez(i,j+1)-ez(i,j))/dy
hy(i,j)=hy(i,j)+(dt/mu0)*(ez(i+1,j)-ez(i,j))/dx
end do
end if
end do
end do
write(1,'(f20.10)') (hx(i,52),i=pmin,pmax)
end
[ 本帖最后由 lanzhiyu520 于 2007-10-24 09:37 编辑 ]
Reply
你不用花时间编二维程序，二维程序葛老师的那本书（西安电子科技大学出版社）后面已经编好了，你可以直接用了，sannxi@163.com
页:
[1]