[FDTD/Fotran] detected reflection near boundary but don't know why

[s_hy]

hi all,

i have wrote codes for 2d fdtd in different permittivity (epsilon). in this code, cell size is 200 x 200, start with eps=1 from center, and different permittivity started at boundary (i1,i2) = (25,75) = (j1,j2), epsilon = 2. the problem is, when the wave propagates and approach (i1,i2) = (25,75) = (j1,j2), i have detected some reflections happened to the wave at timesteps, n after 50, the reflections become obvious...i just wondering, why is that happened. i attached here the animation of the wave and my codes. i really have no idea how to solve this problem.

subroutine test7
implicit none

double precision                                     :: f0,A0
double precision                                     :: delta, deltat
double precision                                     :: eps0,epsr,miu0,miur
double precision                                     :: lamda
double precision                                     :: omega
double precision                                     :: S,k,c
integer                                                     :: i,j
integer                                                     :: ie,je,ic,jc,i1,i2,j1,j2
integer                                                     :: n
double precision,dimension(202,202)   :: Ez,Hy,Hx
double precision,dimension(202,202)   :: Ca,Cb,Da,Db,Dax,Dbx,Day,Dby
!double precision                                     :: Ca,Cb,Da,Db
double precision,dimension(202,202)   :: miu,eps,sigma, sigmat
double precision, parameter                 :: pi = 3.14159265
character(len=20)                                   :: filename

! parameters 
 f0 = 3.0e7
 eps0 = 8.8e-12
 miu0 =  4*pi*1e-7
 c = 3.e8
 lamda = c/f0
 omega = 2*pi*f0
 print *, 'lamda=',lamda
 A0 = 5.0
 k = 2*pi / lamda
 omega = c*k

!courant stability factor 
 !S = 1/sqrt(2.0)
! spatial and time grid steps
 delta = lamda/10
 deltat = delta/(2*c)
 print *, 'deltat=',deltat
 
!model size 
 ie = 100
 je = 100

 !source location
 ic = ie/2
 jc = je/2
 
 !boundary for different permittivity
 i1 = ie/4
 i2 = ie*3/4
 j1 = je/4
 j2 = je*3/4

print*, 'see me?'

!initialize Ez, Hx,Hy to zero at t=0
do i = 1,ie
 do j = 1,je
   Hx(i,j) = 0.0
   Hy(i,j) = 0.0
   Ez(i,j) = 0.0
 end do
end do


!medium's properties
 miur = 1.0
 epsr = 1.0
 
  do i = 1,ie
   do j = 1,je
     !initialize permittivity and permeability
     eps(i,j) = epsr*eps0
     miu(i,j) = miur*miu0
   end do
  end do
 
  do i = i1,i2
    do j = j1,j2
     eps(i,j) = 2*eps0
     miu(i,j) = 1*miu0
    end do
 end do
     
  !medium's properties
! initialize electric conductivity & magnetic conductivity
  do i = 1,ie
   do j = 1,je  
     sigma(i,j) = 0.0
     sigmat(i,j) = 0.0
     Ca(i,j) = (2*eps(i,j)-deltat*sigma(i,j))/(2*eps(i,j)+deltat*sigma(i,j))
     Cb(i,j) = (2*deltat)/(2*eps(i,j)+deltat*sigma(i,j))
     Da(i,j) = (2*miu(i,j)-deltat*sigmat(i,j))/(2*miu(i,j)+deltat*sigmat(i,j))
     Db(i,j) = (2*deltat)/(2*miu(i,j)+deltat*sigmat(i,j))     
   end do
  end do


! beginning of timesteps
 do n = 1,300
   
  write (filename, "('data',I3.3,'.dat')") n
  open (unit=130,file=filename)
 
 
!initiate sinusoidal wavepulse at center
  Ez(ic,jc) = A0*sin(2*pi*f0*n*deltat)
  
! calculate ez-field  
  do i = 2,ie-1
    do j = 2,je-1
          Ez(i,j) = Ca(i,j)*Ez(i,j) + Cb(i,j)*(Hy(i,j) - Hy(i-1,j) - Hx(i,j) + Hx(i,j-1))
          write (130,*) i,j,Ez(i,j)
          if (j == 99) write (130,*) ' '
!          !print *,'i=',i,'j=',j,'Ez(i,j)=',Ez(i,j)
    end do
  end do

  do i = 1,ie-1
    do j = 1,je-1
	  Hx(i,j) = Da(i,j)*Hx(i,j) + Db(i,j)*(Ez(i,j) - Ez(i,j+1))
          !print *,'i=',i,'j=',j,'Hx(i,j)=',Hx(i,j)
     end do
  end do
  
    do i = 1,ie-1
     do j = 1,je-1
          Hy(i,j) = Da(i,j)*Hy(i,j) + Db(i,j)*(Ez(i+1,j) - Ez(i,j))
          !print *,'i=',i,'j=',j,'Hy(i,j)=',Hy(i,j)
     end do
  end do

close (unit=130) 

end do !n

end SUBROUTINE test7

 

[DrClaude]

I'm not familiar with FDTD, but generally speaking grid boundaries act as hard walls, hence waves reaching them will bounce unless some special treatment is used.

Look at the fifth point under http://en.wikipedia.org/wiki/Finite-difference_time-domain_method#Weaknesses_of_FDTD_modeling for references on how to treat that problem.

 

[s_hy]

thank you for the reply. i just found out the problem caused by impedance mismatch. it is happened when i use different epsilon. wave impedance inversely proportional to the square of epsilon. So epsilon =1 wave impedance value is 1.414 times greater than epsilon= 2. If there is any impedance mismatch between the two layers it will create reflection. then, another problem, i found the cause of problem, but have no idea how to correct impedance mismatch in my simulation. do you have any idea?

 

[DrClaude]

I'm sorry, but as I have never used FDTD or worked with that kind of problem, I have no advice to give you.