On Consistent Boundary Conditions for the Yee Scheme in 3D

for the Yee Scheme in 3D

B. Engquist

, J. H¨ aggblad

, O. Runborg

, and A.-K. Tornberg

1 Department of Mathematics and Institute for Computational Engineering and Sciences, The University of Texas at Austin, 1 University Station C1200, Austin TX

78712, U.S.A.

2 Department of Numerical Analysis, CSC, KTH Royal Institute of Technology, SE-10044 Stockholm, Sweden.

3 Swedish e-Science Research Centre, Sweden

4 jonhagg@kth.se

Abstract. The standard staircase approximation of curved boundaries in the Yee scheme is inconsistent. Consistency can however be achieved by modifying the algorithm close to the boundary. We consider a technique to consistently model curved boundaries where the coefficients of the update stencil is modified, thus preserving the Yee structure. The method has pre- viously been successfully applied to acoustics in two and three dimension, as well as electromagnetics in two dimensions. In this paper we general- ize to electromagnetics in three dimensions. Unlike in previous cases there is a non-zero divergence growth along the boundary that needs to be pro- jected away. We study the convergence and provide numerical examples that demonstrates the improved accuracy.

1 Introduction

One of the most popular computational methods for solving Maxwell’s equations is the Yee scheme [1] from 1966, also sometimes called the Finite-Difference Time- Domain method (FDTD) [2]. While more accurate methods have since been de- veloped, the Yee scheme is still the method of choice for many electromagnetic simulations due to its simplicity and efficiency.

The method is based on a straight-forward discretization using compact explicit centred differences on a structured staggered grid. This provides relatively low dispersion—a property important for any simulation involving wave propagation—

and a very low memory footprint.

While it is easy to generate a grid for the method, a downside is that the

structured grid is inefficient in approximating boundaries that do not coincide with

the grid alignment, e.g. curved boundaries. Not only is it inefficient, it is inconsis-

tent [3]. This aspect of the Yee scheme is well-studied [4–9], and many methods

to improve this have been devised, including contour path FDTD [10, 11], locally

conformal FDTD [12–17], and hybrid methods using an unstructured grid close to the boundary [18–20].

In this paper we generalize a method developed in [3], where a simple method was proposed that removes the inconsistencies by modifying the weights on the stencils close to the staircase approximation. This can be done without restricting the CFL-condition and while preserving the standard Yee structure. This is an ap- pealing property, making it simple to incorporate the method into already existing codes.

This paper is organized as follows. After describing the necessary background material in Section 2, we discuss the boundary approximation in Section 3, as well as the generalization from two to three dimensions. Section 4 deals with the issue of divergence preservation, and Section 5 include systematic numerical tests.

2 Maxwell’s equations and the Yee scheme

Maxwell’s equations governing all electromagnetic phenomena, are given by

∂

B + ∇ × E = 0, (1)

∂

D − ∇ × H = −J

, (2)

∇ · D = ρ

, (3)

∇ · B = 0. (4)

where E, D, H, B, J

denote the electric field, electric displacement field, magnetic field, magnetic flux density, and free current density. These fields are for linear isotropic and non-dispersive materials related through

D = E, B = µM,

where  and µ are the permittivity and permeability of the medium, respectively.

We will assume the absence of charges and currents.

The traditional Yee scheme is an explicit discretization of Maxwell’s time- dependent curl-equations (1)–(2), together with alternating leapfrog time-updates,

H

= H

− ∆t

µ ∇ × E

, (5)

E

= E

+ ∆t

 ∇ × H

(6)

where t

= n∆t and E

, H

approximates E and H at t

and t

, respec-

tively. The Yee cell for the spatial discretization is shown in Fig. 1. See e.g. [21] for

a detailed description.

E

E

E

H

H

H

x

y z

Fig. 1. The Yee unit cube, illustrating the spatial discretization of the fields.

3 Consistent boundary modelling

A structured grid leads to efficient simulation of wave propagation in free space, but means that boundaries are modelled by staircase approximations. This is not only a poor approximation, it is inconsistent. Errors of order O(1) are generated at the boundary, destroying convergence in L

and reducing it to O( √

h) in L

. The work in [3] for the two-dimensional case proposes a method where the difference stencils adjacent to a staircase approximated boundary are modified by introducing extra coefficients. Note, this does not change the Yee structure. These are then chosen such that the largest error term in the error expansion disappears, making the boundary approximation consistent.

To illustrate this approach, consider the two-dimensional discretization in Fig. 2.

The relevant boundary condition for the fields on a perfect electric conductor (PEC) surface is

ˆ

n × E = 0, (7)

which says that the tangential electric field on the boundary is zero. Thus the boundary in Fig. 2 is defined by setting the tangential electric field E

and E

to zero. Hence the update stencils for the H

field borders the staircase approximation.

These stencils are in the standard Yee scheme (projected to the xy-plane) given by H

n+¹₂

i+¹₂,j+¹₂

= H

n−¹₂

i+¹₂,j+¹₂

+ ∆t µ

 D

E

n

i+¹₂,j

− D

E

n i,j+¹₂

 . (8)

This reduces to H

n+¹₂

i+¹₂,j+¹₂

= H

n−¹₂

i+¹₂,j+¹₂

+ ∆t µh

 E

n

i+¹₂,j+1

+ E

n i,j+¹₂

 ,

a

a

H

= E

= E

= (Inside)

(Outside)

0 1/2 1 x/h

0 1/2 1 y/h

Fig. 2. Typical staircase discretization of a PEC boundary. The tangential electric field along the dashed line is set to zero.

since the tangential electric field is zero. Focusing on the discretization of the spatial derivative ∂

E

− ∂

E

in (8), one can show that the truncation error for the cell shown in Fig. 2 is

τ

= 1

µh (E

(x

, y

) + E

(x

, y

)) + O(1), (9) where, unfortunately, the 1/h term is not zero. The fields are evaluated at a point (x

, y

) on the boundary. However, if we introduce the generalized stencil

∂

E

− ∂

E

≈ 1 h

 a

E

n

i+¹₂,j+1

− a

E

n

i+¹₂,j

− a

E

n

i+1,j+¹₂

+ a

E

n i,j+¹₂

 , then the truncation error for the same cell takes the form

τ

= 1

µh a

E

(x

, y

) + a

E

(x

, y

)
+ O(1).

If we now consider the rotated coordinate system (ξ, η), where ξ is along the tangent line of the boundary at (x

, y

), we obtain

τ

= 1

µh a

sin α + a

cos α
E ˜

(ξ

, η

) + O(1), (10)

since ˜ E

(ξ

, η

) = 0 by (7), and where α is the angle of the boundary relative the

x-axis. Thus we see that we can set the leading order term in the trunction error

(a) NW (b) NE (c) SW (d) SE (e) W (f) E (g) N (h) S

Fig. 3. The eight different cases for how the boundary can intersect the update stencil for p.

a^Ecos α − a^Ssin α = 0 (NW) a^Wcos α + a^Ssin α = 0 (NE) a^Ecos α + a^Nsin α = 0 (SW) a^Wcos α − a^Nsin α = 0 (SE) a^Ecos α + a^N− a^S
sin α = 0 (W) a^Wcos α − a^N− a^S
sin α = 0 (E) a^E− a^W
cos α − a^Ssin α = 0 (N) a^E− a^W
cos α + a^Nsin α = 0 (S)

Table 1. Summary of the consistency conditions. The labeling corresponds to the one defined in Fig. 3. The angle α of the tangent of the boundary is with respect to the x-axis, i.e., for the NE-case in Fig. 3 then α ∈ (0, π/2). The compass notation is employed for the coefficients in the different directions.

τ

to zero if we satisfy the condition

a

cos α − a

sin α = 0.

This expression is only valid for the discretization shown in Fig. 2, which we denote the SouthEast (SE) case. There are in total 8 possible cases, which are illustrated in Fig. 3. The different consistency conditions corresponding to these cases are enumerated in Table 1. We refer to [3] for a detailed derivation. The stencils in the inner domain stay the same as standard FDTD.

This illustrates how we can remove the inconsistencies of a staircase approxi- mation in two dimensions by carefully weighting the difference stencils along the boundary. Note that the largest error term only depends on the angle with the boundary, not the distance between the stencil and the analytical boundary. The error due to distance is of higher order.

To generalize the method to three dimensions we simply consider each coordi- nate plane at a time. This is possible since the update equation for the magnetic field (5), gives a two-dimensional expression in each component.

Remark: In the paper [22] it was shown that under some conditions there is a

one order gain in the global error compared to the boundary error. Although our

case is not covered by this theory, in the numerical tests this is what is observed,

i.e., the O(1/h) truncation error of (9) translates to O(1) convergence, and the O(1) error in (10) gives O(h) convergence.

4 On the preservation of divergence

In [3] the method was successfully used for the two-dimensional problem, and in [23]

for acoustics in three dimensions. Here, however, we observe an additional phe- nomenon not present before, which is the accumulation of stationary errors along the inner boundary in the H field. The accumulation speed of the error can be seen from Fig. 11 to be linear.

The mechanism behind this error appears to be the slight non-conservation of divergence in the modified cells, leading to a static field localized there. This error lies in the subspace spanned by the eigenfunctions corresponding to the two zero eigenvalues. That standard FDTD (to machine precision) preserves exactly the discrete divergence is well known [21]. That this is not the case for the generalized update formula follows from

∂

∂t I

Unit cell

D · ndS =  ˆ ∂

∂t

 E

2,k+¹₂

− E

2,k+¹₂

 h

+  ∂

∂t

 E

− E

 h

+  ∂

∂t

 E

2,j+¹₂,k+1

− E

2,j+¹₂,k

 h

. Inserting the stencils gives

∂

∂t I

Unit cell

D · n ˆ dS    ≤ Ch, on the inner boundary.

To understand why this error is stationary we recall some basic results from electromagnetics [24]. First, note that Gauss law for the magnetic field means that we have the standard vector potential representation B = ∇ × A. Thus Faraday’s law gives ∇ × (E + A

) = 0, hence

E + A

= ∇φ, (11)

for some scalar potential φ. To see what happens if Gauss law is not fully satisfied, we recall that in general we can always decompose a vector field into its solenoidal (divergence-free) and irrotational (curl-free) components, e.g., B = ∇ × A + ∇ψ.

For this case then Gauss law becomes ∇ · B = ∆ψ 6= 0. But when we insert this expansion into Faraday’s law we get

∇ × (E + A

) + ∇ × ∇ψ

| {z }

=0

= 0,

i.e., we still have the same equation (11) for the dynamics. The same expansion into Amp`ere’s law gives

∇ × ∇ × A + ∇ × ∇ψ

| {z }

=0

= ∇φ

− A

,

again all influence from the extra irrotational term disappears. Thus we see that the time evolution is not affected (in vacuum).

Since the error is in an orthogonal subspace, we can project it away. This requires the solution of an elliptic problem and, for the case of constant coefficients, can even be done in the post-processing stage. The irrotational part of the field needs to be subtracted away according to

H

= H − ∇ψ,

∇

ψ = ∇ · H, (12)

where we need to solve the elliptic equation for ψ first. A simple way to do this in practice is to timestep the parabolic problem ψ

= ∇

ψ − ∇ · H, which is the same as Jacobi-iterating (12). In Particle-in-Cell methods (PIC) this type of correction is refered to as a Langdon-Marder correction [25, 26]. We implement this by iterating according to

H ˜

= H

+ d∆tD

(∇ · H)

2,j+¹₂,k+¹₂

, H ˜

= H

+ d∆tD

(∇ · H)

2,j+¹₂,k+¹₂

, H ˜

2,j+¹₂,k

= H

2,j+¹₂,k

+ d∆tD

(∇ · H)

2,j+¹₂,k+¹₂

, where

(∇ · H)

2,j+¹₂,k+¹₂

= D

H

2,k+¹₂

+ D

H

2,j,k+¹₂

+ D

H

2,j+¹₂,k

, and d = h

/(6∆t). On the staircase boundary we set the condition ∇ · H = 0.

5 Numerical results

To test the method we perform a number of simulations where we initialize the field

to an exact solution for which we have an analytical expression. Then time stepping

the field followed by a comparison to the analytical solution at the time t = T gives

an accurate characterization of the error in the computation. The field is only

measured in the inner domain, sufficiently far away from the outer boundaries that

the effect from these do not influence the result, while still including the staircase

boundary.

xy z ˆ θ

n y

x ˆ ϕ n

Fig. 4. Illustrating the angles θ and ϕ. The normal ˆn is pointing into the domain.

Hx Ez

Ey

α (a) yz-plane

Hy

Ez

Ex

β

(b) xz-plane Hz

Ex

Ey

γ (c) xy-plane

Fig. 5. The projected angles α, β and γ, as well as how the boundary plane intersects the individual H-stencils.

5.1 Test case 1: constant field

The domain is set to Ω = [0, 2π]

where we let a plane Γ cut through the inner domain. Γ is parametrized by

z = (x − x

) tan β + (y − y

) tan α, which means the normal vector is given by

ˆ

n = 1

p1 + tan

α + tan

β



 tan β tan α

−1



 .

For the cases considered we will formulate the solution in terms of the angles ϕ and θ together with the normal vector

ˆ n =





cos ϕ sin θ sin ϕ sin θ

cos θ



 , 0 ≤ ϕ ≤ 2π, 0 ≤ θ ≤ π,

as well as for simplicity restrict ourselves to θ ∈ (0, π/4), ϕ ∈ (5π/4, 3π/2). To see how these two sets of angles are related we project according to ˆ n

= (ˆ n · y)ˆ ˆ y + (ˆ n · z)ˆ ˆ z = (0, sin ϕ sin θ, cos θ)

, thus giving

α = arccos n ˆ

· (0, 0, 1)

p sin

ϕ sin

θ + cos

θ

!

= arctan(− sin ϕ tan θ).

Similarly we get β = arctan(− cos ϕ tan θ). We note that α is the angle of the line

through the H

update stencil, which is in the (y, z)-plane, and β is the angle of the

line through the H

stencil, which is in the (x, z)-plane. Furthermore, we denote the

angle of the line intersecting the H

stencil by γ, which is in the (x, y)-plane and

is with respect to the x-axis and given in terms of ϕ as γ = ϕ − π/2. See Fig. 4–5.

N kH − Hexactk∞ p kE − Eexactk∞ p

20 3.07(-01) 3.57(-01)

60 1.91(-01) 0.432 5.77(-01) -0.437 180 1.42(-01) 0.269 5.16(-01) 0.101

N kH − Hexactk2 p kE − Eexactk2 p

20 1.02(-01) 1.09(-01)

60 5.94(-02) 0.495 1.01(-01) 0.072 180 3.38(-02) 0.512 5.39(-02) 0.567

Table 2. The resulting error after time stepping the constant solution (5.1–5.1) using the standard Yee method on a N³ grid. Given in the column p is the convergence rate. The parameters for the internal boundary are ϕ = 11π/8, θ = π/7, ¯x = −1/2 and ¯y = −1/2.

N kH − Hexactk∞kE − Eexactk∞

20 2.25(-16) 1.59(-16) 60 7.17(-16) 2.83(-16) 180 2.21(-15) 1.00(-15)

N kH − Hexactk2 kE − Eexactk2

20 6.05(-17) 4.35(-17) 60 1.04(-16) 6.07(-17) 180 1.89(-16) 1.47(-16)

Table 3. Same as Fig. 2 but using modified coefficients. The error is close the machine.

In the first test of the generalization to three dimensions we time step a constant solution. To find such a solution consider the boundary condition for the tangential electric fields (7)





cos ϕ sin θ sin ϕ sin θ

cos θ



 ×



 E

E

E



 =





E

sin ϕ sin θ − E

cos θ E

cos θ − E

cos ϕ sin θ E

cos ϕ sin θ − E

sin ϕ sin θ



 = 0.

Setting E

= 1 gives the solution for both fields as

E = (− tan β, − tan α, 1), H = 0.

The results after time stepping on a grid of resolution N

are shown in Table

2 and 3. We clearly see that we have constant error for the standard Yee scheme,

while the error is down to machine precision when using the modified coefficients.

H

E

k ˆ

n ˆ E

H

k ˆ

xy z

Fig. 6. Illustrating the time-harmonic field vectors used in (5.2–5.2) and (5.2–5.2). The incoming wave has the electric field polarized in the z-direction and the magnetic field polarized in the xy-plane, transverse to the boundary. The incoming wave vector is parallel to the plane of incidence.

5.2 Test case 2: time-harmonic wave field on plane surface

As a second test case we consider the reflection of a plane wave. Thus let the incoming wave be given by

E

= −ik ˆ p

e

, H

= k

µω k ˆ

× E

,

where the polarization and wave vectors are given by ˆ p

= (0, 0, 1)

and ˆ k

=

−ˆ n

= (− cos ϕ, − sin ϕ, 0)

. Thus the reflected wave is given by E

= −ik ˆ p

e

e

, H

= k

µω k ˆ

× E

,

where ˆ k

= (− cos ϕ cos 2θ, − sin ϕ cos 2θ, sin 2θ) and ˆ p

= (ˆ k

× p ˆ

) × ˆ k

. Note that the amplitude for the reflected magnetic field can be simplified to

(−ik)k/(µω) = −iω,

and that the cross products between the polarization and wave vectors have the components

k ˆ

× p ˆ

= (− sin ϕ, cos ϕ, 0)

ˆ

p

= (− cos ϕ sin 2θ, sin ϕ sin 2θ, cos 2θ)

.

These wave and polarization vectors are illustrated in Fig. 6.

In Fig. 7–8 we see some qualitative results of the field after time stepping. Clearly we have a large improvement for modified coefficients. A time-series of the error is shown in Fig. 9. The convergence is shown in Fig. 10, where the E field is evidently O(h) both in L

and L

. As for the H-field, the error is reduced by one order of magnitude, but the convergence seems to flatten out for very small discretization lengths h, indicating that the error is still O(1), albeit much smaller. The unmodified Yee scheme has O( √

h) convergence in L

, and O(1) in L

, as expected. The solutions in Fig. 7–10 are computed without correcting for divergence, as it is still small enough not to influence. To illustrate the growth in the divergence for long simulation times, we plot the L

-norm of the field when initialized with random data in Fig. 11–12. First without, then with a Langdon-Marder correction. We clearly see that the growth disappears in the second case.

5.3 Test case 3: harmonic wave field on a perfectly conduct- ing sphere

As another test case we consider the scattering of plane waves against a metallic (PEC) sphere. This can be solved analytically and is called Mie scattering. See [27]

for an in-depth treatment.

The incoming plane wave and the reflected wave is expanded as E

= E

e

ˆ e

= E

e

∞

X

n=1

i

2n + 1 n(n + 1)



M

(k) − iN

(k) 

H

= H

e

e ˆ

= − kE

µω e

∞

X

n=1

i

2n + 1 n(n + 1)



M

(k) + iN

(k) 

and

E

= E

e

∞

X

n=1

i

2n + 1 n(n + 1)

 a

M

(k) − ib

N

(k)  ,

H

= − kE

µω e

∞

X

n=1

i

2n + 1 n(n + 1)

 b

M

(k) + ia

N

(k)  , where

M

(k) = ± 1

sin θ z

(kr)P

(cos θ) cos

sin ϕˆ e

− z

(kr) ∂P

(cos θ)

∂θ

sin

cos ϕˆ e

,

x/π

y/π

0.5 1 1.5

0.6 0.8 1 1.2 1.4

−1.5

−1

−0.5 0 0.5 1 1.5

(a) Hxstandard Yee.

x/π

y/π

0.5 1 1.5

0.6 0.8 1 1.2 1.4

−1.5

−1

−0.5 0 0.5 1 1.5

(b) Hx coefficients modified.

x/π

y/π

0.6 0.8 1 1.2 1.4 0.6

0.8 1 1.2 1.4

−0.5 0 0.5

(c) Hy standard Yee.

x/π

y/π

0.6 0.8 1 1.2 1.4 0.6

0.8 1 1.2 1.4

−0.6

−0.4

−0.2 0 0.2 0.4 0.6

(d) Hycoefficients modified.

x/π

y/π

0.6 0.8 1 1.2 1.4 0.6

0.8 1 1.2 1.4

−0.1 0 0.1 0.2

(e) Hz standard Yee.

x/π

y/π

0.6 0.8 1 1.2 1.4 0.6

0.8 1 1.2 1.4

−0.02 0 0.02 0.04 0.06 0.08

(f) Hz coefficients modified.

Fig. 7. Comparison of magnetic field for standard Yee and modified coefficients. Again the O(1/h) truncation error generated at the boundary is clearly visible.

x/π

y/π

0.6 0.8 1 1.2 1.4 0.6

0.8 1 1.2 1.4

−0.6

−0.4

−0.2 0 0.2 0.4 0.6

(a) Ex standard Yee.

x/π

y/π

0.6 0.8 1 1.2 1.4 0.6

0.8 1 1.2 1.4

−0.2

−0.1 0 0.1 0.2

(b) Excoefficients modified.

x/π

y/π

0.5 1 1.5

0.6 0.8 1 1.2 1.4

−1

−0.5 0 0.5 1

(c) Eystandard Yee.

x/π

y/π

0.5 1 1.5

0.6 0.8 1 1.2 1.4

−0.6

−0.4

−0.2 0 0.2 0.4 0.6

(d) Eycoefficients modified.

x/π

y/π

0.5 1 1.5

0.6 0.8 1 1.2 1.4

−1.5

−1

−0.5 0 0.5 1 1.5

(e) Ez standard Yee.

x/π

y/π

0.5 1 1.5

0.6 0.8 1 1.2 1.4

−1

−0.5 0 0.5 1

(f) Ez coefficients modified.

Fig. 8. Comparison of electric field for standard Yee and modified coefficients. The qual- itative improvement is noticeable.

0 0.05 0.1 0.15 0.2 0.25 0.3 0

0.05 0.1 0.15 0.2 0.25

t/π

Max norm error

Yee Mod

(a) H.

0 0.05 0.1 0.15 0.2 0.25 0.3

0 0.05 0.1 0.15 0.2

t/π

Max norm error

Yee Mod

(b) E.

0 0.05 0.1 0.15 0.2 0.25 0.3

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08

t/π L2 norm error

Yee Mod

(c) H.

0 0.05 0.1 0.15 0.2 0.25 0.3

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08

t/π L2 norm error

Yee Mod

(d) E.

Fig. 9. Time series of the error in the solution of the reflected time-harmonic wave. Solid line is modified coefficients, dashed line is standard Yee. The resolution is N = 160.

10¹ 10² 10⁻¹

10⁰

PPW

Max norm error

Yee Mod

(a) H.

10¹ 10²

10⁻¹ 10⁰

PPW

Max norm error

Yee Mod

(b) E.

10¹ 10²

10⁻² 10⁻¹

PPW L2 error

Yee Mod

(c) H.

10¹ 10²

10⁻² 10⁻¹

PPW L2 error

Yee Mod

(d) E.

Fig. 10. Convergence of the time-harmonic field reflected on the plane, i.e., the solution given by (5.2–5.2) and (5.2–5.2). The dashed lines indicate slope 1/2 and 1. (◦) is standard Yee and (×) is modified coefficients.

20 40 60 80 100 0.5

1 1.5 2

t L2 norm difference

(a) H.

20 40 60 80 100

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16

t L2 norm difference

(b) E.

Fig. 11. L2 norm of the field when initialized with random data. Shown is the difference compared to the initial norm value. The norm for the standard Yee method is shown as the (in average) constant line, which is included as a comparison.

20 40 60 80 100

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18

t L2 norm difference

(a) H

20 40 60 80 100

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16

t L2 norm difference

(b) E

Fig. 12. Projecting away the divergence using the Langdon-Marder correction. Compare to Fig. 11, where we see that the difference is the removal of the linear growth.

N

1ne

(k) = n(n + 1)

kr z

(kr)P

(cos θ) sin cos ϕˆ e

+ 1 kr

∂

∂r [krz

(kr)] ∂P

(cos θ)

∂θ sin cos ϕˆ e

± 1

kr sin θ

∂

∂(kr) [krz

(kr)]P

(cos θ) cos sin ϕˆ e

. and

a

= − j

(α) h

(α) , b

= − (αj

(α))

 αh

(α) 

, α = ka, with r = a being the radius of the sphere.

These functions M and N are called spherical vector functions. The subscript o,e is used to identify wheter odd (sin) or even (cos) functions are used. The subscript 1n refers to the the index on P

. For o the sign of ± is +, and for e it is −.

The superscript refers to the type of spherical Bessel functions z

used. In (5.3–5.3) spherical Bessel functions j

(kr) of the first kind are used and in (5.3–5.3) spherical Bessel functions of the third kind (Hankel functions of the first kind) h

(kr) are used.

These Bessel functions are defined by z

(x) = r π

2x Z

(x), H

(x) = J

(x) + iY

(x).

Y

(x) = J

(x) cos νπ − J

(x) sin νπ

J

(x) =  x 2



X

k=0

 −



k!Γ (ν + k + 1)

where ν can be noninteger and Γ is the Gamma function. The functions P

are associated Legendre functions of the first kind, obtained from the Legendre poly- nomials by

P

(x) = (1 − x

)

d

dx

P

(x), where P

in turn is defined as

P

(x) = 1 2

n!

 d

dx

(x

− 1)



.

x/π

y/π

−0.5 0 0.5

−0.5 0 0.5

−0.5 0 0.5

(a) Ex

x/π

y/π

−0.5 0 0.5

−0.5 0 0.5

−0.5 0 0.5

(b) Ex, modified coefficients.

x/π

y/π

−0.5 0 0.5

−0.5 0 0.5

−0.15

−0.1

−0.05 0 0.05 0.1 0.15

(c) Ey

x/π

y/π

−0.5 0 0.5

−0.5 0 0.5

−0.15

−0.1

−0.05 0 0.05 0.1 0.15

(d) Ey, modified coefficients.

x/π

y/π

−0.5 0 0.5

−0.5 0 0.5

−0.25

−0.2

−0.15

−0.1

−0.05 0 0.05

(e) Ez

x/π

y/π

−0.5 0 0.5

−0.5 0 0.5

−0.15

−0.1

−0.05 0 0.05

(f) Ez, modified coefficients.

Fig. 13. Comparing the calculated electric field from (5.3–5.3). The grid resultion is N = 140.

x/π

y/π

−0.5 0 0.5

−0.5 0 0.5

−0.3

−0.2

−0.1 0 0.1 0.2 0.3

(a) Hx

x/π

y/π

−0.5 0 0.5

−0.5 0 0.5

−0.15

−0.1

−0.05 0 0.05 0.1 0.15

(b) Hx, modified coefficients.

x/π

y/π

−0.5 0 0.5

−0.5 0 0.5

−1

−0.5 0 0.5 1

(c) Hy

x/π

y/π

−0.5 0 0.5

−0.5 0 0.5

−0.5 0 0.5 1 1.5 2

(d) Hy, modified coefficients.

x/π

y/π

−0.5 0 0.5

−0.5 0 0.5

−0.5 0 0.5

(e) Hz

x/π

y/π

−0.5 0 0.5

−0.5 0 0.5

−1

−0.5 0 0.5 1 1.5

(f) Hz, modified coefficients.

Fig. 14. Comparing the calculated magnetic field from (5.3–5.3). The grid resultion is N = 140.

Note that different normalizations P

occur in the litterature. Here the Con- don–Shortley phase (−1)

is omitted. The derivatives of the associated Legendre functions are obtained by the relation

dP

(cos θ)

dx = 1

2 (n − m + 1)(n + m)P

(cos θ) − P

(cos θ) .

The update stencils are modified in the same way as with a plane inner boundary.

The angle used is the tangent line of the intersection of the sphere and the stencil plane at the grid points for the field components (H

, H

, H

) for which we modify the components.

Qualitative results of the computed fields are shown in Fig. 13–14. We see, similarly to the plane boundary case, that the quality of the solution is noticeably higher when we modify the coefficients.

6 Conclusions

We have shown that the method devised in [3] can successfully be applied to Maxwell’s equations in three dimensions. Despite the added difficulty of divergence preservation, accuracy is increased up to one order of magnitude with minimal modifications of the Yee scheme. For longer simulation times the divergence growth can successfully be projected away. Numerical tests show a clear improvement in quality of the computed solution, where the leading high-frequency errors generated by the staircased boundary are removed. The CFL condition of the Yee scheme can be kept intact, making the method efficient and keeping the dispersion errors low.

An improved convergence rate is also observed, especially in L

.

References

1. Yee, K.S.: Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media. IEEE Trans. Antennas Propag. 14 (1966) 302–307 2. Taflove, A.: Application of the finite-difference time-domain method to sinusoidal

steady-state electromagnetic-penetration problems. IEEE Trans. Electromagn. C.

EMC-22 (1980) 191–202

3. Tornberg, A.K., Engquist, B.: Consistent boundary conditions for the Yee scheme. J.

Comput. Phys. 227(14) (2008) 6922–6943

4. Cangellaris, A., Wright, D.: Analysis of the numerical error caused by the stair-stepped approximation of a conducting boundary in FDTD simulations of electromagnetic phenomena. IEEE Trans. Antennas Propag. 39 (1991) 1518–1525

5. Ditkowski, A., Dridi, K., Hesthaven, J.S.: Convergent cartesian grid methods for Maxwell’s equations in complex geometries. J. Comput. Phys. 170(1) (2001) 39–80 6. Hao, Y., Railton, C.: Analyzing electromagnetic structures with curved boundaries

on Cartesian FDTD meshes. IEEE Trans. Microw. Theory Tech. 46(1) (1998) 82–88

7. Holland, R.: Pitfalls of staircase meshing. IEEE Trans. Electromagn. C. 35(4) (1993) 434–439

8. Jurgens, T., Taflove, A., Umashankar, K., Moore, T.: Finite-difference time-domain modeling of curved surfaces. IEEE Trans. Antennas Propag. 40(4) (1992) 357–366 9. Railton, C., Schneider, J.: An analytical and numerical analysis of several locally

conformal FDTD schemes. IEEE Trans. Microw. Theory Tech. 47(1) (1999) 56–66 10. Jurgens, T., Taflove, A.: Three-dimensional contour FDTD modeling of scattering

from single and multiple bodies. IEEE Trans. Antennas Propag. 41(12) (1993) 1703–

1708

11. Railton, C., Craddock, I.: Stabilised CPFDTD algorithm for the analysis of arbitrary 3D PEC structures. IEE Proc. Microwaves Antennas Propag. 143(5) (1996) 367–372 12. Dey, S., Mittra, R.: A locally conformal finite-difference time-domain (FDTD) al- gorithm for modeling three-dimensional perfectly conducting objects. IEEE Microw.

Guided. W. 7(9) (1997) 273–275

13. Dey, S., Mittra, R.: A modified locally conformal finite-difference time-domain al- gorithm for modeling three-dimensional perfectly conducting objects. Microw. Opt.

Technol. Lett. 17(6) (1998) 349–352

14. Yu, W., Mittra, R.: A conformal FDTD algorithm for modeling perfectly conducting objects with curve-shaped surfaces and edges. Microw. Opt. Technol. Lett. 27(2) (2000) 136–138

15. Tolan, J.G., Schneider, J.B.: Locally conformal method for acoustic finite-difference time-domain modeling of rigid surfaces. J. Acoust. Soc. Am. 114(5) (2003) 2575–2581 16. Zagorodnov, I., Schuhmann, R., Weiland, T.: A uniformly stable conformal FDTD-

method in Cartesian grids. Int. J. Numer. Model. 16(2) (2003) 127–141

17. Zagorodnov, I., Schuhmann, R., Weiland, T.: Conformal FDTD-methods to avoid time step reduction with and without cell enlargement. J. Comput. Phys. 225(2) (2007) 1493–1507

18. Rylander, T., Bondeson, A.: Stable FEM-FDTD hybrid method for Maxwell’s equa- tions. Comput. Phys. Commun. 125(1-3) (2000) 75–82

19. Monorchio, A., Mittra, R.: A hybrid finite-element finite-difference time-domain (FE/FDTD) technique for solving complex electromagnetic problems. IEEE Microw.

Guided. W. 8(2) (1998) 93–95

20. Wu, R.B., Itoh, T.: Hybrid finite-difference time-domain modeling of curved surfaces using tetrahedral edge elements. IEEE Trans. Antennas Propag. 45(8) (1997) 1302–

1309

21. Taflove, A., Hagness, S.C.: Computational Electrodynamics: The Finite-Difference Time-Domain Method. 3 edn. Artech House Publishers (2005)

22. Gustafsson, B.: The convergence rate for difference approximations to mixed initial boundary value problems. Math. Comp. 29(130) (1975) pp. 396–406

23. H¨aggblad, J., Engquist, B.: Consistent modeling of boundaries in acoustic finite- difference time-domain simulations. (submitted) (2012)

24. Jackson, J.D.: Classical Electrodynamics. 3 edn. Wiley (1998)

25. Langdon, A.B.: On enforcing Gauss’ law in electromagnetic particle-in-cell codes.

Comput. Phys. Commun. 70(3) (1992) 447–450

26. Marder, B.: A method for incorporating Gauss’ law into electromagnetic PIC codes.

J. Comput. Phys. 68(1) (1987) 48–55

27. Stratton, J.A.: Electromagnetic Theory. McGraw-Hill Companies (1941)