arXiv:1111.4289v2 [physics.comp-ph] 21 Nov 2011

Interfacial Numerical Dispersion and New
Conformal FDTD Method
Axman Fisher
yjliu@fudan.edu.cn

Abstract
This article shows the interfacial relation in electrodynamics shall
be corrected in discrete grid form which can be seen as certain numerical dispersion beyond the usual bulk type. Furthermore we construct
a lossy conductor model to illustrate how to simulate more general
materials other than traditional PEC or simple dielectrics, by a new
conformal FDTD method which main considers the effects of penetrative depth and the distribution of free bulk electric charge and current.

Contents
1 Introduction

2

2 Grid representation of interfacial relation

3

3 Lossy conductor model and conformal FDTD

4

4 Discussion

7

1

1

Introduction

Bulk dynamics and boundary conditions are two essential ingredients
for solving a concrete physical problem. In computational electromagnetic, it is well known finite difference form of wave equation in
bulk leads to numerical dispersion which is measured by the factor
sin x
x , where x is proportional to the grid size. If the grid length is
λ
smaller than 12
, then the dispersive factor approximates to 1 which
means the error due to this numerical dispersion can always be ignored in problem. However, beyond this traditional bulk numerical
dispersion (BND), we find discrete process in simulation as in FiniteDifference Time-Domain (FDTD) approach also generates certain unphysical anti-dispersion effects on the boundary which can be defined
as the interfacial numerical dispersion (IND). If we ignore this effect,
we shall get the staircase approximation grid representation of the
curved boundary, which is just one of the shortcomings of traditional
FDTD method. Instead, we find interfacial relation in usual electrodynamics must be corrected in grid variable form due to IND, and some
treatments in conformal FDTD method before as in perfect electric
conductor (PEC) or dielectric materials can be seen as special cases of
the new relation.
Traditional conformal FDTD (CFDTD) often deal with PEC [110] and simple dielectrics [11-15]. In former case, the electric field in
material can be ignored due to the penetrative depth is far less than
σ
 1 and in later
one grid length due to the lossy parameter p = ω
σ
case simple denote near lossless, i.e. p = ω  1.
For general dispersive and lossy media like Lorentz model, the related Conformal FDTD approach is difficult due to the permittivities
in the interfacial relation involve complex dependence of the frequency
and in corresponding time-domain algorithm this shall turns to be involve with high order time derivatives.
However, for a large class of lossy dielectrics, whose penetrative
depth δ is between PEC and perfect electric insulator (PEI), more exσ
actly in which the lossy parameter p = ω
∼ 1 or the δ is about several
grids’ length, if we can ignore the influence of the resonance frequency
in the problem or the involved electromagnetical frequency is far away
from the material resonance frequencies, we can use certain conductivity model as the dispersive model and get a (theoretical) more effective
CFDTD method called lossy conductor CFDTD (LCCFDTD), which
consider the effects of IND, surface/volume free electric current/charge,
penetrative depth and conformal angle (defined as the angle between
the cross-section of interface and the grid line) in more physical way.
Especially our method is more suitable for the FDTD simulation of
the problems like the design of microwave dark chamber with strong
absorber materials.

2

2 Grid representation of interfacial relation
Consider the 2-dimensional TE case below without loss of generality.
First we define the angle θ between the tangent line at the interfacial point and the grid line as the conformal angle. The grid line at
the boundary is divided into two parts on which two electric fields is
uniform distributed as E in air and E 0 in materia. General case of
arbitrary two medias is similar.
For PEC model, the inside field can be assumed to be vanished. But
for general materials or even as in nondispersive, lossless dielectrics, to
obtain a conformal algorithm instead of naive staircase approximation,
we have to discuss the interfacial relation between E and E 0 . In continuous electrodynamics theory we have following boundary relations:
Eq = Eq0
Hq =
D⊥ =
B⊥ =

Hq0 + α
0
D⊥
+ρ
0
B⊥

(1)
(2)
(3)
(4)

where α is line density of free surface electrical current and ρ is surface
density of free electrical charge. Here we only consider the simplest case
in which α, ρ vanish and permittivity is a real const. In next section
we may show how to deal with nonzero α, ρ and complex permittivity
for a more general model of materials.
As we know, for certain frequency only two of above interfacial
relations are independent and here we use (1) (3), i.e.
Eq = Eq0 ,

0
E⊥ = εr E⊥
, εr = ε/ε0 .

This interfacial relation indicate E, E 0 can’t be along same direction
through the interface, i.e. the refraction phenomena exist or for general
case if the permittivity as function of the frequency, electromagnetical
waves of different frequency shall refract to different directions, i.e.
dispersive phenomena appears. However, in discrete grid simulation,
E, E 0 have to be along same grid line and this in fact introduces
extra error. Like the discrete treatment for wave equation in bulk
will produces numerical dispersion, the effect of this elimination of
physical refraction or dispersion on interface due ro simulation method
as FDTD can be seen as certain interfacial numerical dispersion (IND).
Hence we should take one of two fields on the same boundary grid line
as projection of the refraction field along the grid line, take E for
instance and the relation between E, E 0 give the grid representation
of physical interfacial relation (1), (3).

3

By a simple geometrical analysis indicated by figure 1 above, we
find this grid interfacial relation,
E = (cos2 θ + εr sin2 θ)E 0 = εr (θ)E 0 .

(5)

It’s direct to see several special correct results. First take εr = 1 i.e. for
grid in same materia, E = E 0 as our basic assumption in simulation.
Second, consider θ = 0, i.e. the coincident case of grid choice, E = E 0 ,
as same as traditional FDTD. Third, if θ = π2 , as the translation
case for coincident grid choice, we get E = εr E 0 or equivalently as
0
Eq = Eq0 = 0, D⊥ = D⊥
.

3 Lossy conductor model and conformal
FDTD
To deal with more general materials as needed in most applications,
we construct a lossy conductor model and for simplicity, we assume
the boundary of materia is flat, i.e the tangent line at the interfacial
point is coincident with the boundary. New model main considers
a reasonable generalization of the surface distribution of free electric
charge and current in well conductor model (WEC) via a special discrete treatment.
First we introduce a kind of lossy conductor conformal grid (LCCG),
and take a 2-dimensional cross-section of a rotated planar interface for
example as figure 2 shows. Except usual field nodes, this grid have
four material points a, b, c, d and a new grid parameter B called grid
penetrative number. For the grid outside of the material B = 0 and
on the geometrical boundary of media B = 1. For the grids in the
material, B = min([Di /ds]) + 1, where Di is the distance between the
grid’s center node and the tangent line of interface, [x] is the maximal
integer less than x, ds denote the size of grid. In fact we can take its
relative coordinate of the nearest interfacial grid as its B value.
Only two material points in a grid are independent which are determined by the dividing character of the boundary grids, while other
two points if needed can be obtained via translation of the line which
connects the two independent points as figure 2 shows. In fact like

4

traditional subgrid method which introduces four field nodes in a grid
to construct the relative small grids, our new conformal method also
consider new sub nodes which describe the geometry of material boundary and this process don’t introduces new field nodes to increase the
complexity of computation.
The start point of conformal method is the integrated form of
Maxwell equation below:
Z
I
∂
µH · dS = −
E · dl.
(6)
∂t S
∂S
For the case of figure 2, by discretion of above equation we have iterated
equation for magnetic field H,
n+ 12

Hz

0
4t
1
1
[Exn (i, j + )lx (i, j) + Exn (i, j + )lx0 (i, j)
µS(i, j)
2
2
0
1
1 0
n
n
− Ex (i, j − )lx (i, j − 1) − Ex (i, j − )lx (i, j − 1)
2
2
0
1
1
n
n
+ Ey (i − , j)ly (i − 1, j) + Ey (i − , j)ly0 (i − 1, j)
2
2
0
1
1
n
n
− Ey (i + , j)ly (i, j) − Ey (i + , j)ly0 (i, j)]
2
2

n− 21

(i, j) =Hz

(i, j) +

Note we can just take some terms in above formula to be zero for
other cases of division. In the situation of figure 2, we also have some
constraints on the length variables due to two independent materia
points in a grid,
ly0 (i, j)
,
ly (i − 1, j)
ly (i, j)
lT (i, j)
= 0
.
lx + lx0 = ∆x, ly + ly0 = ∆y; 0
lT (i, j)
ly (i, j)

ly0 (i, j) = ly0 (i − 1, j), lx0 (i, j − 1) = lx (i, j) ·

It’s easy to see that the usual PEC conformal algorithm is just as the
special case for which S(i, j) denote the area of the outside cell and
1
1
1
1
Ex (i, j + ) = Ex0 (i, j + ), Ey (i − , j) = Ey0 (i − , j);
2
2
2
2
1
1
0
0
Ex (i, j − ) = 0, Ey (i + , j) = 0.
2
2
By the symmetry of electric grids and magnetic grids in FDTD
method, it’s not necessary to do conformal treatment for electric grids,
i.e. the iteration of electrical field is same as in traditional FDTD
which is just by discretion of another Maxwell integrated equation:
I
Z
Z
∂
H · dl =
D · ds +
J · ds.
∂t S
∂S
S
Note in our lossy conductor model, the conductivity σ is not zero which
require the term of free electric current.
To deal with the effects of the free electric charge and penetrative
depth, we take a special discrete treatment for the distribution of free

5

bulk electric current in our grid system, which assumes the free electric
charge only distribute uniformly on the transverse grid line lT and the
electric current just conduct along the perpendicular direction of the
grid plane. Thus no electric charge and current on the grid line in
our model and especially the transverse line lT in one grid just denote
the unit transverse line which appears in the definition of free surface
electric current α in (2).
By a simple deduction like for the lossless model, we can get the
generalization of formula (5) in lossy conductor model,
E = (cos2 θ + (εr −

|α|
σr
2
0
+
0 ) sin θ)E
jω ε0 vE⊥

= (cos2 θ + ε̄rc sin2 θ)E 0 = ε̄rc (θ)E 0 .
where |α| = σE 0 lT as the electric current through the unit transverse
line lT and v is the velocity of the electromagnetical wave in the mac
terial and clearly v = <(n
= √εrcµr L , where the refraction lossy
r)
q √
1+ 1+p2
, p = εσrrω , σr = σ/ε0 . Finally we obtain
parameter L =
2
the grid interfacial relation of lossy conductor model:
√
lT σr εr µr L
σr
+
) sin2 θ)E 0
(7)
E = (cos2 θ + (εr −
jω
c sin θ
To apply above formula in conformal FDTD scheme, we have to
make some approximations and furthermore turn it to the iteration
relation in time domain.
• Half order and IND corrected approximation (HIND): lT = 0
(ignore the distribution of free electrical current)
σr
) sin2 θ)E 0
(8)
E = (cos2 θ + (εr −
jω
∂E
∂E 0
= (cos2 θ + εr sin2 θ)
− σr sin2 θ · E 0 .
(9)
∂t
∂t
• First order and IND corrected approximation (FIND):
lT 6= 0, L ' 1
√
lT σr εr µr
σr
2
E = (cos θ + (εr −
+
) sin2 θ)E 0
jω
c sin θ
√
lT σr εr µr
∂E
∂E 0
= (cos2 θ + (εr +
) sin2 θ)
− σr sin2 θ · E 0 .
∂t
c sin θ
∂t

(10)
(11)

• Second order and IND corrected approximation (SIND):
2
lT 6= 0, L ' 1 + p8
√
lT σr εr µr
σr
σ2
2
E = (cos θ + (εr −
+
(1 + 2r 2 )) sin2 θ)E 0 (12)
jω
c sin θ
8εr ω
√
2
2
lT σr εr µr
∂ E
∂ E0
∂E 0
2
2
2
)
sin
θ)
=
(cos
θ
+
(ε
+
−
σ
sin
θ
·
r
r
∂t2
c sin θ
∂t2
∂t
√
lT σr3 εr µr
−
sin θE 0 .
(13)
8cε2r

6

• Geometric boundary conformal FDTD:
Only consider above conformal treatment for the geometrical boundary grid in which B = 1. Thus in the iteration equation of magnetic
field we have
1
1
1
1
Ex (i, j − ) = Ex0 (i, j − ), Ey (i + , j) = Ey0 (i + , j)
2
2
2
2
and Ex (i, j + 12 ), Ey (i− 12 , j) can be determined by Ex0 (i, j + 12 ), Ey0 (i−
1
2 , j) via the discrete form of the equation (9), (11) and (13). For B > 1,
we set E = E 0 which means we have ignored the effect of penetrative
depth in this scheme.
• Penetrative boundary conformal FDTD:
Consider above conformal treatment for the penetrative boundary
δ
], ∆l = min(∆x, ∆y). For such grids we
grid in which 1 ≤ B ≤ [ ∆l
0
get all E from E via discrete form of (9), (11) and (13). The grid
δ
with B = [ ∆l
] + 1 will be treated as geometric boundary grid while for
δ
δ
B > [ ∆l ] + 1 we just set E = E 0 . Especially for PEC case, [ ∆l
]=0
0
and E = 0 as in traditional way.
Although the penetrative depth δ is an important physical quantity
which can be determined by the electromagnetic parameters of the
material exactly, we can take it as an adjustable parameter called lossy
depth δl for practical use. It can be set by the precision and cost of
computation in a problem and penetrative depth δ could just be its
default value.

4

Discussion

We know electrodynamics manifest Lorentz symmetry or Maxwell equations in continuous spacetime is covariant under the translation and
rotation, while in discrete grid spacetime such Poincare symmetry is
broken. For flat boundary material, we find formula (5) considers the
exact effect of conformal angle in the interfacial relation such that the
simulation physical result is independent of grid system in the sense
of global translation and rotation, or in other words for flat boundary
case the shape of material in simulation is invariant of grid translation
and rotation which is more physical unlike in traditional FDTD. In this
manner we could say our IND corrected conformal FDTD have restored
certain 3-dimensional space symmetry and especially it manifests the
consistence between the conformal and nonconformal algorithm.
There are several directions to generalize above section’s result.
First, if the boundary is curved, the tangent line at the interfacial
point shall not be along lT . One simple way to deal with this case is
to define the angle between the tangent line and transversal line lT as
relative conformal angle θr and the conformal angle still to be the one
between lT and the grid line. If θr ≤ θ2 , such boundary’s treatment
can be approximated by transversal line boundary as done in previous,
while θr > θ2 , we’d better approximate the boundary in this direction
by the complete grid line, i.e. E = E 0 .

7

Second, our LCCFDTD method only considers the electrical loss,
it seems the magnetic loss and some special anisotropy or nonlinear or
dispersive medias also can be treated in similar way [16]. It may also
have applications in other numerical computational areas.
Finally, we may consider the stability problem as in traditional conformal FDTD method [17-18] and an improved absorb boundary may
be obtained by the combination of a large lossy conductor material and
the ideas of PML. Of course all above discussions are main in theoretical manner and it can be verified by some numerical experiments. The
results may be similar with [14], i.e. the simulation of rotated problem
should has good agreement with the coincidence problem as reference
solution.

References
[1] Thomas G. Jurgens, Allen Taflove et al., Finite-Difference TimeDomain Modeling of Curved Surfaces, IEEE Trans. Antennas
Propagat., 40(1992) 357.
[2] T. Jurgens and A. Taflove, ”Three-dimensional contour FDTD
modeling of scattering from single and multiple bodies,” IEEE
Trans. Antennas and Prop., vol. 41, pp. 1803-1708, Dec. 1993.
[3] S.Dey, R.Mittra, S.Chebolu, A technique for implementing the
FDTD algorithm on a nonorthogonal grid, IEEE Trans. Antennas
Propagat., 40(1992)357.
[4] S. Dey and R. Mittra, ”A locally conformal FDTD algorithm for
modeling 3D PEC objects,” IEEE Microwave and Guided Wave
Letters, vol. 7, pp. 273-275, Sept. 1997.
[5] C. J. Railton and J. Schneider, ”An analytical and numerical analysis of several locally conformal FDTD schemes,” IEEE Trans. On
Microwave Theory and Techniques, vol. 47, pp. 56-66, Jan. 1999.
[6] Mark W.Steeds, et al., A Comparison of Two Conformal Methods
for FDTD Modeling, IEEE Trans. Electr. Compat., 38(1996)181.
[7] Y. Hao and C.J. Railton, Analyzing electromagnetic structures
with curved boundaries on Cartesian FDTD meshes, IEEE Trans
Microwave Theory Tech 46(1998), 82-88.
[8] I.A.Zagorodnov, R.Schuhmann, A uniformaly stable conformal FDTD-method in Cartesian grids, Int. J. Numer. Model.,
16(2003)127.
[9] I.A.Zagorodnov, R.Schuhmann, Conformal FDTD-methods to
avoid time step reduction with and without cell enlargement, J.
Comput. Phys., 225(2007)1493.
[10] Tian Xiao and Qing Huo Liu, A 3-D Enlarged Cell Technique
(ECT) for the Conformal FDTD Method, IEEE Trans. Antenn.
Prop., 225(2007)1493.

8

[11] R. D. Meade, et al., ”Accurate theoretical analysis of photonic
band-ap materials,” Phys. Rev. B, vol. 48, pp. 8434-8437, 1993.
[12] N. Kaneda, B. Houshmand, and T. Itoh, FDTD analysis of dielectric resonators with curved surfaces, IEEE Trans Microwave
The- ory Technol 45 (1997), 1645-1649.
[13] S. Dey and R. Mittra, A locally conformal finite difference time
domain technique for modeling arbitrary shaped objects, IEEE
Antennas Propagat Soc Int Symp, June 1998, vol. 1, pp. 584-587.
[14] J.-Y. Lee and N.-H. Myung, ”Locally tensor conformal FDTD
method for arbitrary dielectric surfaces,” Microwave and Optical
Technology Letters, vol. 23, No. 4, pp. 245-249, Nov. 20, 1999.
[15] Nadobny J, Sullivan D, Wlodarczyk W, et al.,A 3-D tensor FDTD
formulation for treatment of sloped interfaces in electrically inhomogeneous media, IEEE T. Antennas and Propagation, vol. 51
(8), pp. 1760-1770, August 2003.
[16] Luebbers R., Lossy dielectrics in FDTD, IEEE transactions on
antennas and propagation, 1993.
[17] Pereda J.A., Garcia O., Vegas A., Prieto A., Numerical dispersion
and stability analysis of the FDTD technique in lossy dielectrics,
Microwave and Guided Wave Letters, IEEE, Jul. 1998.
[18] A.C. Cangellaris and D.B. Wright, Analysis of the numerical error caused by the stair-stepped approximation of a conduction
boundary in FDTD simulations of electromagnetic phenomena,
IEEE Trans Antennas Propagat 39(1991), 1518-1525.

9

