Two-Dimensional FDTD Simulation

Once again, we will start with the normalized Maxwell’s equations, assume that the structure being modeled extends to infinity in the z-direction. If the incident wave is also uniform in the z-direction, then all partial derivatives of the fields with respect to z must equal zero. Under these conditions, the full set of Maxwell’s curl equations for TMz mode.

As in one-dimensional simulation, it is important that there is a systematic interleaving of the fields to be calculated. This is illustrated in Fig.1. Putting Eqs. (1)~(3) into the finite differencing scheme results in the following difference equations.

Fig.1 Interleaving of the E and H fields for the two-dimensional TM formulation

The program fd2d_001 implements the above equations. It has a simple Gaussian pulse source that is generated in the middle of the problem space. Fig.2 demonstrates a simulation for the different time steps
!* FD2D_001 2D FDTD simulation in free space *!
Module common_data
implicit none
!Calculate region
integer,save:: Imin,Imax,Jmin,Jmax
!constant
real,parameter:: Eps0=8.854e-12,Mu0=1.25663706e-6
real,parameter:: C0=3.e8,Z0=377.,Pi=3.1415926536d0
!source parameter
real,save ::SourceT0,SourceT
!MAXIMUM NUMBER OF TIME STEPS
integer,parameter:: timestop=250
!grid size,timestep
real,save:: dx,dy,dt
!variables for Field
real,dimension(:,:),allocatable,save:: Ez,Hx,Hy
endmodule common_data

program main
use common_data
implicit none
integer ntime
call pre_processing
!MAIN LOOP FOR FIELD COMPUTATIONS AND DATA SAVING
do Ntime=1,timestop
call Ezfd2d(Ntime)
call Hxfd2d
call Hyfd2d
enddo
call savedata
end

subroutine pre_processing
use common_data
implicit none
Imin=-150;Imax=150
Jmin=-150;Jmax=150
allocate (Ez(Imin:Imax,Jmin:Jmax))
allocate (Hx(Imin:Imax,Jmin:Jmax))
allocate (Hy(Imin:Imax,Jmin:Jmax))
Ez=0.;Hx=0.;Hy=0.
dx=1.e-2;
dy=dx
dt=dx/2./c0
SourceT0=120;SourceT=50
return
end

subroutine Ezfd2d(n)
use common_data
implicit none
integer i,j,n
do i=Imin,Imax-1
do j=Jmin,Jmax-1
Ez(i,j)=Ez(i,j)+(Hy(i+1,j)-Hy(i,j))*dt/dx/Eps0
* -(Hx(i,j+1)-Hx(i,j))*dt/dy/Eps0
enddo
enddo
Ez(0,0)=exp(-((n-1)*1.-SourceT0)**2/SourceT**2)
return
end

subroutine Hxfd2d
use common_data
implicit none
integer i,j
do i=Imin,Imax-1
do j=Jmin+1,Jmax-1
Hx(i,j)=Hx(i,j)-(Ez(i,j)-Ez(i,j-1))*dt/dy/Mu0
enddo
enddo
return
end

subroutine Hyfd2d
use common_data
implicit none
integer i,j
do i=Imin+1,Imax-1
do j=Jmin,Jmax-1
Hy(i,j)=Hy(i,j)+(Ez(i,j)-Ez(i-1,j))*dt/dx/Mu0
enddo
enddo
return
end

subroutine savedata
use common_data
implicit none
integer j
open(1,file='ez.m')
open(2,file='xgrid')
open(3,file='ygrid')
write(1,*) 'z=['
do j=Jmin,Jmax-1
write(1,100) Ez(Imin:Imax-1,j)
enddo
write(1,*) ']'
write(2,*) (j,j=Imin,Imax-1)
write(3,*) (j,j=Jmin,Jmax-1)
close(1)
close(2)
close(3)
100 format(99999(1x,E12.4))
return
end

Fig.2 Results of a simulation using the program fd2d_001.