1 FDTD discretization
In finite-difference time-domain (FDTD), the discretization of the
spatial domain is driven by the mesh definition. The mesh or grid is
comprised of grid cells, with the size of these cells establishing the
simulation(‘)s spatial
resolution. The conventional grid cell dimensions in FDTD varies from
(lambda/10) to (lambda/20), where (lambda) is the smallest wavelength in the
simulation source. It is possible for a coarse mesh using the Yee grid
to introduce a nonphysical numerical dispersion, which results in the
simulated wave travelling with a slower phase velocity. While it is
important to note that structures often have fine features or curved
surfaces that may require a finer mesh, this article will focus on the
simulation resolution relative to the source wavelength.
1.2 Numerical dispersion
Solving Maxwell(‘)s equations
using the Yee grid [1] results in a discrepancy between the analytical
and simulated solutions. In this regard, the phase velocity of the
numerical wave deviates slightly from its true value depending on the
frequency, grid resolution, and the direction of the propagating wave.
To compare the analytic and simulated phase velocities, one must first
compute the wavevector (k)), which
is related to the phase velocity as [begin{equation}
v_{p}=frac{omega}{k}
end{equation}] where (omega) is the angular frequency.
1.2.1 Analytical Dispersion
Equation
For the comparison of the analytical and simulation results this
article will consider a propagating monochromatic planewave in
freespace. In the 2D domain, this planewave ((E_{y}=E_{y_ {0}}e^{jleft(omega
t-k_{x}x-k_{z}zright)})) along with Maxwell(‘)s equations yields the analytic
dispersion equation[1-2] [begin{equation}
label{analytic}k_{x}^{2}+k_{z}^{2}=left(frac{omega}{c}right)^{2}
end{equation}] where (c)
represent the speed of light and the analytical wavevector components
along x and z are defined as (k_{x})
and (k_{z}), respectively.
1.2.2 Numerical Dispersion
Equation
For the simulated case, Maxwell(‘)s equations need to be expressed in
their discretized form. The source free expressions for the TE
polarization (i.e. Hx, Hz, Ey) are written as [1-2] [begin{eqnarray}
label{continuous1}frac{partial E_{y}}{partial
t}&=&frac{1}{varepsilon_{0}}left
(frac{partial H_{x}}{partial z}-frac{partial H_{z}}{partial
x}right)\
label{continuous2}frac{partial H_{x}}{partial
t}&=&frac{1}{mu_{0}}left(frac
{partial E_{y}}{partial z}right)\
label{continuous3}frac{partial H_{z}}{partial
t}&=&frac{1}{mu_{0}}left(frac
{-partial E_{y}}{partial x}right)
end{eqnarray}] Discretizing the partial derivatives in the
above equations is achieved through the use of the Yee cell [1-2] and
yields [begin{eqnarray}
label{discretized1} E_{y}^{n+1}(i,k) & = &
E_{y}^{n}(i,k)+frac{Delta t}{varepsilon_{0},Delta
z}left[H_{x}^{n+0.5}(i,k)-H_{x}^{n+0.5}(i,k-1)right]nonumber \
& & -frac{Delta t}{varepsilon_{0},Delta
x}left[H_{z}^{n+0.5}(i,k)-H_{z}^
{n+0.5}(i-1,k)right]\
label{discretized2}H_{x}^{n+0.5}(i,k)&=&H_{x}^{n-0.5}(i,k)+frac{Delta
t}{mu_{0},Delta z}left[E_
{y}^{n}(i,k+1)-E_{y}^{n}(i,k)right]\
label{discretized3}H_{z}^{n+0.5}(i,k)&=&H_{z}^{n-0.5}(i,k)-frac{Delta
t}{mu_{0},Delta x}left[E_
{y}^{n}(i+1,k)-E_{y}^{n}(i,k)right]
end{eqnarray}] where the discretized field components are
written as [2] [begin{eqnarray}
label{component1}E_{y}mid_{I,K}^{n}&=&E_{y_{0}}e^{jleft(omega
nDelta t-tilde{k}_{x}
IDelta x-tilde{k}_{z}KDelta zright)}\
label{component2}H_{x}mid_{I,K}^{n}&=&H_{x_{0}}e^{jleft(omega
nDelta t-tilde{k}_{x}
IDelta x-tilde{k}_{z}KDelta zright)}\
label{component3}H_{z}mid_{I,K}^{n}&=&H_{z_{0}}e^{jleft(omega
nDelta t-tilde{k}_
{x}IDelta x-tilde{k}_{z}KDelta zright)}
end{eqnarray}] where (tilde{k}_{x}) and (tilde{k}_{z}) are numerical wavevectors
in x and z direction. The indices I and K denote the spatial indices in
the x and z directions, respectively. Inserting the discretized fields
into the discretized Maxwell(‘)s
equations leads to [begin{equation}
left[left(frac{1}{cDelta t}right)sinleft(frac{omegaDelta
t}{2}right)right]^
{2}=left[frac{1}{Delta x}sinleft(frac{tilde{k_{x}}Delta
x}{2}right)right]^{2}
+left[frac{1}{Delta z}sinleft(frac{tilde{k_{z}}Delta
z}{2}right)right]^{2}
end{equation}] It can be shown that the discretized dispersion
equation (12) simplifies to the analytical expression (2) for a square
grid ((Delta x = Delta z = Delta))
under the limit that (Delta->0)
and using the Courant stability criteria ((Delta t = frac{Delta}{sqrt{2}c})).
1.2.3 Numerical Dispersion Relative
to Frequency and Resolution
To explore the deviation between analytical and discretized
dispersion equations relative to frequency and resolution it is easier
to reduce the problem to 1D ((k_z=0))
. This action modifies the analytical (1) and discretized (12) equations
to be [begin{eqnarray}
frac{omegaDelta x}{c} & = & |k_xDelta
x|qquadtext{analytical}\
frac{omegaDelta x}{c} &=&
left|frac{2}{S}sin^{-1}left(Scdot sinleft(frac
{k_xDelta x}{2}right)right)right|qquadtext{discretized}
end{eqnarray}] where (S=frac{cDelta t}{Delta x}). Note that
the analytical expression has been normalized by (Delta x) for the purpose of comparison.
Plotting both of these expressions demonstrates a number of trends in
how the discretized dispersion equation deviates from the analytical
one, see figure 1.
The observations that can be made are that the deviation:
- proportional to the frequency.
- proportional to the size of the grid cells (inversely proportional
to the number of points per wavelength). - inversely proportional to the ratio of (Delta t) to (Delta_x), i.e. (S). Note that while S = 1 would match the
analytical and discretized results this is at the edge of numerical stability and as such is not
advisable.
1.2.4 Numerical Dispersion Relative
to Propagation Direction
To explore numerical dispersion relative to propagation direction the
analytical and discretized expressions can be rewritten relative to a
general wavevector (k) such that
[begin{eqnarray}
k_x &=& kcosphi\
k_z &=& ksinphi
end{eqnarray}] With this identity the analytical expression (2)
becomes [begin{eqnarray}
left(frac{omegaDelta_x}{c}right)^2 & = & k^2(cos^2phi +
sin^2phi) \
left(frac{omegaDelta_x}{c}right)^2 & = & k^2
end{eqnarray}] and similarly the discretized expression (12),
on a square grid, becomes [begin{eqnarray}
left(frac{Delta}{cDelta t}right)^2sin^2left(frac{omegaDelta
t}{2}right) & = &
sin^2left(frac{tilde{k}
Delta cosphi}{2}right)+sin^2left(frac{tilde{k}
Delta sinphi}{2}right)
end{eqnarray}] Solving (19) with respect to (tilde{k}) requires the use of Newton’s
method on the following expression [begin{equation}
label{dispersion}fleft(tilde{k}right)=sin^{2}left(frac{tilde{k}Deltacosvarphi}{2}right)+sin^{2}left(frac{tilde{k}Deltasinvarphi}{2}right)-left(frac{Delta}{cDelta
t}right)^{2}sin^{2}left(frac{pi cDelta t}{lambda_{0}}right)=0
end{equation}] The algorithm [2] is given as [begin{equation}
label{Newton}tilde{k}_{n+1}=tilde{k}_{n}-frac{fleft(tilde{k}right)}{f^{‘}left(tilde{k}right)}=tilde{k}_{n}-frac{sin^{2}left(C_{1}tilde{k}_{n}right)+sin^{2}left(C_{2}tilde{k}_{n}right)-C_{3}}{C_{1}sinleft(2C_{1}tilde{k}_{n}right)+C_{2}sinleft(2C_{2}tilde{k}_{n}right)}
end{equation}] where (f^{‘}) represents the first derivative
of the function (f) and
[begin{eqnarray}
C_{1}&=&frac{Deltacosvarphi}{2}\
C_{2}&=&frac{Deltasinvarphi}{2}\
C_{3}&=&left(frac{Delta}{cDelta
t}right)^{2}sin^{2}left(frac{pi cDelta t}
{lambda_{0}}right)
end{eqnarray}] The numerical phase velocity ({tilde{v}_{p}}=frac{omega}{tilde{k}})
can now be calculated from (1). Figure 2 shows the normalized numerical
phase velocity versus the angle of propagation, (varphi), for three cases i.e. when the
number of grid cells per wavelength are 4, 8, and 16. In all plots the
time step is set as (S=frac{cDelta
t}{Delta}=frac{1}{sqrt{2}}). As can be seen from the results,
the numerical phase velocity is lower than its true value, c (shown by
the dashed line), and that the deviation between the analytical and
discretized results can be minimized by increasing the number of grid
cells.
