High Order Finite Difference Approximations of Electromagnetic Wave Propagation Close to Material Discontinuities

High Order Finite Difference Approximations of

Electromagnetic Wave Propagation Close to

Material Discontinuities

Jan Nordström and Rikard Gustafsson

The self-archived postprint version of this journal article is available at Linköping

University Institutional Repository (DiVA):

http://urn.kb.se/resolve?urn=urn:nbn:se:liu:diva-68595

N.B.: When citing this work, cite the original publication.

The original publication is available at www.springerlink.com:

Nordström, J., Gustafsson, R., (2003), High Order Finite Difference Approximations

of Electromagnetic Wave Propagation Close to Material Discontinuities, Journal of

Scientific Computing, 18, 215-234. https://doi.org/10.1023/A:1021149523112

Original publication available at:

https://doi.org/10.1023/A:1021149523112

Copyright: Springer (part of Springer Nature) (Springer Open Choice Hybrid

Journals)

High Order Finite Difference Approximations of

Electromagnetic Wave Propagation Close to

Material Discontinuities

Jan Nordström1 , 2_{and Rikard Gustafsson}2

1_{The Swedish Defense Research Agency (FOI), Aerodynamics Division (FFA),}

Computa-tional Aerodynamics Department, SE-172 90 Stockholm, Sweden. E-mail: Jan.Nordstrom @foi.se

2_{The Department of Scientific Computing, Information Technology, Uppsala University, Box}

337, SE-751 05 Uppsala, Sweden.

Electromagnetic wave propagation close to a material discontinuity is simulated by using summation by part operators of second, fourth and sixth order accuracy. The interface conditions at the discontinuity are imposed by the simultaneous approximation term procedure. Stability is shown and the order of accuracy is verified numerically.

KEY WORDS: Maxwell’s equations; computational electromagnetics; high

accuracy; finite difference methods; interface conditions.

1. INTRODUCTION

One approach to minimize radar reflections is to introduce one or several layers of absorbing materials at the surface of the scattering object. In this paper we simulate electromagnetic wave propagation close to such material discontinuities. High order finite difference methods, see [1–11], that guarantee bounded error growth in time for realistic meshes will be used.

More precisely we will apply the technique developed in [12] for boundary conditions and block interfaces to the problem with interface conditions at material discontinuities. We focus on geometrically relatively simple problems and aim for high accuracy. Related work for less accurate methods with focus on complex geometries are presented in [13–17].

The Laplace transform technique for stability (see [18]) is to complex for high order methods and multiple-dimensions. Our strategy is to combine the energy-method, finite difference operator of summation-by-parts (SBP) form with a penalty formulation for the interface conditions and gain control over the boundary terms at the interface. That will lead to long time stability and accurate results.

2. EQUATIONS AND INTERFACE CONDITIONS

The physics of electromagnetic fields are described by the Maxwell’s equations m“H “t=−N × E, E “E “t=N × H − J, (1) N · EE=r, N · mH=0, (2) and the equation of continuity

“r

“t+N · J=0. (3)

Here E is the electric field, H the magnetic field, J is the electric current density and r is the charge density. E and m are permittivity and

permea-bility respectively.

The interface conditions for E at a surface of discontinuity are

E2t− E1t=0, E2E2n− E1E1n=sf, (4)

where the subscripts 1, 2 referes to the different sides of the discontinuity. The subscripts t, n referes to the tangential and normal components respectively. In (4),sfis the free surface charge density.

The interface conditions for H are

H2t− H1t=K, m2H2n− m1H1n=0, (5)

where K is the surface current density at the interface.

3. THE CONTINUOUS PROBLEM

With J=0 we can formulate the 2D version of (1) in TE (transverse electric) mode as

Slut+Aux+Buy=0 Srvt+Avx+Bvy=0,

north south east west x=--1 x=1 y=1 u v

Fig. 1. The computational domain and the grid.

where the subscripts l, r refer to the left and right part of the computational domain (see Fig. 1). In (6), u, v are the solutions in the left and right part of the computational domain, and

A=

R

S

R

S

R

S

R

S

R

S

R

S

Also, the conditions at the material discontinuity with sf=0, K=0 (see

Eqs. (4) and (5)) are

u(x=0)=Gv(x=0), where G=

R

S

3.1. The Interface Conditions

Three physical interface conditions are given in (7). To determine how many and which of these conditions to impose at the interface x=0, we apply the standard Laplace (and Fourier) transform technique (see [18]). The general solution is obtained from the ansatzkexp(ox) inserted into the

Laplace and Fourier transformed version of (6). The decaying part (away from x=0) of the solutions can be written,

u˜ˆ=sl

R

S

R

S

where s,w are the dual variables to t, y respectively, R(s) > 0 and o+

l=+`Elmls2+w2, o−r=−`Ermrs2+w2.

The form (8) of the solution has two undetermined coefficients which leads to the first result: two of the three interface conditions in (7) should be imposed at x=0.

The two interface conditions lead to an equation of the form Es=0 where E(s) is a 2 × 2 matrix and s=(sl, sr)T. The first and third condition,

the first and second condition and the second and third condition in (7) leads to,

det(E)13=s(El`Ermrs2+w2+Er`Elmls2+w2),

det(E)12=0,

det(E)23=iw(El`Ermrs2+w2+Er`Elmls2+w2),

(9)

respectively. A well-posed bounded solution is obtained only if det(E(s)) ]0 for R(s) > 0. This requirement (see (9)) is fulfilled only for the

combi-nation of the first and third condition in Eq. (7). We conclude that the interface conditions that lead to a unique and bounded solution are,

Cu(x=0)=Cv(x=0), where C=

R

1 0 0 0 0 0 0 0 1

S

Remark. Another way to arrive at (10) is to apply the energy method

on (6) with periodic outer boundary conditions, the result

(||u||2 S_l+||v|| 2 S_r)t=F 1 0 (vT_{Av − u}T_Au) x=0 dy,

The Maxwell’s equations without material discontinuities (as well as any system of equations with constant coefficients) are conservative and it is interesting to see whether the discontinuity and the interface conditions (10) preserve that property. Multiply Eq. (6) with a smooth continuous test functionF(x, y, t) that vanishes on the outer boundaries and at t=0, ..

Integration in time and space, the fact that F is continuous over the

material interface and vanishes at the outer boundaries and at t=0, ., leads to F . 0 F1 0 F0 −1 (FT tSlu+F T xAu+F T yBu) dx dy dt +F . 0 F1 0 F1 0 (FT tSrv+FTxAv+F T yBv) dx dy dt +F . 0 F1 0 FT_{A(v − u)} x=0dy dt=0. (11)

The interface condition (10) removes the interface terms at x=0.

3.2. The Energy Estimate

By introducing W=(u, v)T_{, A¯ =diag(A, A) and B¯=diag(B, B), S=}

diag(Sl, Sr), we can write Eq. (6) as

SWt+A¯ Wx+B¯ Wy=0. (12)

Multiplication of Eq. (12) with WT_{, integration in space and the fact that}

the matrices A and B are symmetric, leads to

d dt||W|| 2 S=− 2 F 1 0 HzEy|x=1dy+2 F 1 0 HzEy|x=−1dy

z

z

z

z

where ||W||2

S=WTSW. The problem (12) augmented with initial condition

is well-posed if we specify the Ey component at the west and east boundary

and the Excomponent at the south and north boundary.

Remark. By specifying the ingoing characteristic variable we can

obtain a strongly well-posed method, see [18].

3.3. The Divergence of E and H

Consider the continuous problem, away from the interface. The divergence of (1) and Eq. (3) leads to

(N · EE − r)(t)=(N · EE − r)(0), (N · mH)(t)=(N · mH)(0). (14) This means that ifN · EE − r=0, N · mH=0 initially, Eq. (2) will be satisfied for all t. Equation (14) show that Eq. (2) follows from Eqs. (1) and (3).

4. THE DISCRETE PROBLEM

The spatial discretization is made by introducing a grid, see Fig. 1, with node points (xi, yj), where i=0, 1,..., N and j=0, 1,..., M. The

unknowns are organized in a vector u with the following enumeration

u=

R

S

R

S

R

S

R

S

The vector ui, jcorresponds to the grid point (xi, yj). At the interface x=0, uNand v0correspond to the same grid points.

Let U, DU be the numerical approximation of the scalar quantities

u, ux respectively. The approximation of the first derivative DU is

intro-duced as

DU=P−1_{QU, (15)}

where P and Q are matrices. If a spatial operator is of the form (15) and the conditions (i) and (ii) below are full-filled the operator is referred to as a Summation By Parts operator (SBP) (see [1, 2]).

(i) The matrix P is symmetric, positive definite and bounded,

DxpI [ P [ DxqI, where p > 0 and q are bounded independent of N.

(ii) The matrix Q is almost skew-symmetric. Q+QT₌

diag(−1, 0,..., 0, 1).

The boundary and interface conditions will be imposed by the SAT (Simultaneous Approximation Term) procedure, see [9].

The semi-discrete equation for the two domain problem becomes

Slut+(P −1 xlQxlé IMé A) u+(INlé P −1 y Qyé B) u =BC_l+(P−1

xlENlNlé IMé Sx=0l)(eNlNlé ((IMé C)(uNl− v0))) Srvt+(P −1 xrQxré IMé A) v+(INré P −1 y Qyé B) v =BC_r+(P−1 xrE0Nré IMé Sx=0r)(e0Nré ((IMé C−1)(v0− uNl))). (16)

The symbol é denotes the Kronecker product, see Appendix A. In (16),

IL are identity matrices of size (L+1) × (L+1). The structure of E0L and ELL with size (L+1) × (L+1), where L can be Nl, Nr or M are,

E0L=

R

S

R

S

e0Land eLL are vectors of length L+1.

e0L=(1 0 · · · 0)T, eLL=(0 0 · · · 1)T. (18)

BCl and BCr are the penalty terms including the outer boundary

conditions, BCl=(P

z

z

z

BCr=(P

z

z

z

Here Eij and eij are defined as in Eq. (17) and (18). The boundary data are

stored in gi.

To visualize how the two regions are connected to each other we can write the system (16) in the form

Wt=M(W+G), W=(u, v)T, (19)

where G includes the outer boundary data. In the 6th order case, with

M=20, Nl=20, Nr=20, the matrix M is schematically depicted in Fig. 2.

The square in the center of the matrix is the 36M × 36M matrix which couples the two regions.

4.1. The Interface Conditions

Multiply Eq. (16) with FT_P

l and FTPr respectively. We have used

the notation Pl=(Pxlé Pyé I3) and Pr=(Pxré Pyé I3) and introduced

Fig. 2. The matrix in the 6th order case with M=20, Nl=20, Nr=20. The 36M × 36M square matrix in the center couples the two regions.

F(x, y, t) which is the restriction to the grid of a smooth continuous test

function that vanishes at the outer boundaries and at t=0, .. The sum-mation by parts properties of Ql and Qr, integration in time, and the fact

thatF is continuous over the interface, leads to

F . 0 (FT tPlSlu+uT(Qxlé Pyé A) F+uˆT(Qyé Pxlé B) Fˆ ) dt +F . 0 (FT tPrSrv+vT(Qxré Pyé A) F+vˆT(Qyé Pxré B) Fˆ ) dt +F . 0 F T_(x=0)(P yé (Sx=0lC − Sx=0rC − A)) uNldt +F . 0 F T_(x=0)(P yé (Sx=0 rC − Sx=0lC+A)) v0dt=0. (20)

Note the close similarity between Eqs. (11) and (20). The relation

Sx=0 r=Sx=0l− A (21)

eliminates the interface terms. In (21) we have introduced the scaling

S_x=0r=S_x=0rC, S_{x=0 l}=S_x=0lC.

4.2. The Energy Estimate

With the use of the so called Q-formulation, see [12], Eq. (16) can be written in the following way

PW_t+(Q_x+S_int) W+Q_yW=BC, (22)

where W=(u, v)T_and

P=diag((Pxlé Pyé I3) Sl, (Pxré Pyé I3) Sr)

Qx=diag((Qxlé Pyé A), (Qxré Pyé A))

Qy=diag((Pxlé Qyé B), (Pxré Qyé B)).

The penalty terms BC including boundary conditions imposed by the SAT procedure are,

BC=(Elé E0Nlé Pyé Sx=−1)(Eu₀W − elé e0Nlé g1)

+(Elé E0Mé Pxlé Sy=0l)(Euˆ₀W − elé g2lé e0M)

+(Elé EMMé Pxlé Sy=1l)(Euˆ_MW − elé g3lé eMM)

+(Eré E0Mé Pxré Sy=0r)(Evˆ₀W − eré g2ré e0M)

+(Eré EMMé Pxré Sy=1r)(Evˆ_MW − eré g3ré eMM)

+(E_ré E_NrNré P_yé Sx=1)(EvNrW − eré eNrNré g4). (23) Eij in Eq. (23) is defined as in Eq. (17), Eui is defined such that the only

non-zero elements in Eu_iW is the ones corresponding to ui, and

El=

r

s

r

s

r

s

r

s

The penalty matrices at the outer boundaries are Si. Sint is the penalty

matrix at the interface between the two regions.

Now let a ‘‘tilde’’ sign indicate the (6M × 6M) ( for the sixth order case) block that couples the solutions in the left and right domains, i.e.,

Sint=

r

s

r

s

r

s

Now we can split Q4x+S4int into a symmetric and a skew-symmetric part as Q4x+S4int= (Q4_x+S4_int) − (Q4_x+S4_int)T 2 + (Q4_x+S4_int)+(Q4_x+S4_int)T 2

z z

r

s

r

s

By introducing the interface condition (21), requiring Sx=0l and Sx=0 r to be symmetric, i.e., Sx=0l=

r

s

r

s

x and D4. The final form of the operators becomes Q4sk x= 1 2

r

s

r

s

To get an energy estimate we multiply Eq. (22) with WT_{and add its}

trans-pose. We obtain d dt||W|| 2 P=WTPWt+W T tPTW, − WT_Q yW − (WTQyW)T+WTBCn, s+(WTBCn, s)T =uˆT

0(Pxlé (B+Sy=0+STy=0)) uˆ0

+uˆT

M(Pxlé ( − B+Sy=1+STy=1)) uˆM

+vˆT

0(Pxré (B+Sy=0+STy=0)) vˆ0

+vˆT

M(Pxré ( − B+Sy=1+STy=1)) vˆM

− 2uˆT

0(Pxlé Sy=0) g2l− 2uˆTM(Pxlé Sy=1) g3l

− 2vˆT 0(Pxré Sy=0) g2r− 2vˆMT(Pxré Sy=1) g3r, (26) − WT_(Q x+Sint) W − (WT(Qx+Sint) W)T+WTBCe, w+(WTBCe, w)T =uT 0(Pyé (A+Sx=−1+STx=−1)) u0 +vT Nr(Pyé ( − A+Sx=1+S T x=1)) vNr +2uT 0(Pyé Sx=−1) g1− 2vTM(Pyé Sx=1) g4− WTint2D4Wint, (27)

where Wint=(uNl, v0)T. We now choose our penalty parameters, at the

outer boundary, such that the quadratic terms in (26) and (27) vanishes, i.e.,

Sx=−1=

r

s

r

s

r

s

r

s

The final version of the energy estimate has the form

d dt||W|| 2 P=BT+IT, (29) where BT=− 2uT 0(Pyé Sx=−1) g1− 2vTnr(Pyé Sx=1) g4

z

z

z

z

Note the similarity between Eqs. (13) and (30). To get a stable method we must make sure that D4 is positive semidefinite. We achieve this by choosing

l1, l2[0 in Eq. (24).

Remark. Note that only the two colocated vectors uNl, v0 are

involved in the stability condition (31). For other types of methods, when the operators do not posses the SBP property, or even worse when using staggered grids, the stability condition involves vectors at a number of locations which makes it virtually impossible to prove stability.

Remark. A strictly and strongly stable method can be obtained by

specifying the ingoing characteristic variable at the outer boundaries, see [18].

4.3. The Spectrum

To analyze the numerical stability we consider the matrix A − 2Sx=0l

in Eq. (25) and the spectrum of the matrix M in Eq. (19). Note that a posi-tive semi definite D4 or

A − 2Sx=0l=

r

s

is required for stability. The spectrum depends on the penalty parameters in Sx=0l, Sx=0r and is depicted in Figs. 3 and 4. The parameter values

l1=0, l2=0 leads to D4=0 which means that we have zero dissipation.

The spectrum is shown in Fig. 3. The parameter values l1=−0.1,

l2=−0.1 introduces dissipation as can be seen in Eq. (32) and Fig. 4.

5. NUMERICAL EXPERIMENTS

To evaluate the method we consider a plane wave which propagates over the discontinuity. The outer boundary data are given by the analytic solution. For the time integration we use the classical fourth order Runge– Kutta method.

--1 -- 0.5 0 0.5 1

-- 50 0 50

-- 8 -- 6 -- 4 -- 2 0 2 -- 50

0 50

Fig. 4. Spectrum for matrix M in Eq. (19) forl1=−0.1, l2=−0.1.

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4

2nd order, whole domain

t

||

⋅ ||

|| ∇ ⋅ E || || error ||

5.1. The Divergence

Consider Eq. (1) in the semi discrete case. The fact that our SBP operators satisfy the property DxDy=DyDxleads to

(Nd· EE − r)(t)=(Nd· EE − r)(0),

i.e., the same result as in the continuous case. In the computations the ana-lytic solution is used as initial data and r=0. Away from the interface

where E=const., we have Nd· E(0)=N · E3(0)+O(Dx)p. Here E3 is the exact

solution projected on the grid. This means that Nd· E=O(Dx)p initially,

and thusNd· E=O(Dx)pfor all t.

Figures 5–7 show the development of the divergence and the L2 error

in the whole computational domain for long times. The divergence stays constant and the growth of the L2error is limited. As mentioned above, the

analytic solution at t=0 is divergence free. Hence the divergence shown in the figure is introduced by the numerical computation of N · E as was discussed above. 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05

4th order, whole domain

t

||

⋅ ||

|| ∇ ⋅ E ||

|| error ||

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5x 10

-- 3 _{6th order, whole domain}

t

||

⋅ ||

|| ∇ ⋅ E || || error ||

Fig. 7. ||N · E||2and ||error||2as a function of time, 6th order accuracy. 5.2. The Accuracy

The order of accuracy in space, will be computed by mesh refinement in space for a given time. The solutions are computed with a very small time step. When the wave hits the discontinuity with an angle of incidence that prohibits a reflecting wave, that special angle is known as the Brewster

angle or angle of polarization. Gp is the result of Brewster’s law [19]

tanGp=n2/n1. n1 and n2 are the refraction index of respective medium.

The wave is depicted in Fig. 8. The accuracy in space is given by Fig. 9. The angleGi=p6, lead to a significant reflection at the interface between the

two medium, see Fig. 10. The accuracy in space is shown Fig. 11.

6. CONCLUSIONS AND DISCUSSION

Wave propagation close to a material discontinuity have been simu-lated. SBP operators of second, fourth and sixth order accuracy have been used. The boundary and interface conditions have been imposed by the

Fig. 8. (a) the wave propagating to the left, (b) the wave propagating to the right, (c) the total wave,Gi=Gp. 101 102 10-- 8 10-- 6 10-- 4 10-- 2 100 Θi = Brewster angle ∆ x --1_, ∆ y --1 error acc 2 acc 4 acc 6 ref. slope 2 ref. slope 4 ref. slope 6

Fig. 10. (a) the wave propagating to the left, (b) the wave propagating to the right, (c) the total wave,Gi=p6. 101 102 10--8 10-- 6 10-- 4 10-- 2 100 Θi = π/6 ∆ x--1_{,∆ y}-- 1 error acc 2 acc 4 acc 6 ref. slope 2 ref. slope 4 ref. slope 6

SAT procedure. We have shown that we can control the stability by choosing penalty parameters in the SAT procedure.

The simulations were made for two angles of incidence to the material interface between the two media. For each case we have verified the order of accuracy. The divergence has been shown to be constant in time for uniform meshes. A material discontinuity will cause no problem when simulating electromagnetic wave propagation using SBP operators and a weak form of interface conditions.

Extension to three dimensions poses no theoretical problems, stability and accuracy follows by using the technique discussed above. However, added complexity will be introduced by geometries that require transfor-mation of the piecewise constant coefficient problems into variable coeffi-cient problems with dicontinuities via the metric coefficoeffi-cients. That problem is of a different kind and is considered in ongoing work.

APPENDIX A. THE KRONECKER PRODUCT

Definition 1. Let A be a p × q matrix and let B be an m × n matrix,

then A é B=

R

S

The p × q block matrix A é B is called a Kronecker product.

There are a number of rules for Kronecker products, see [20], we will present some of them. Let A, B, C and D be matrices of arbitrary sizes, such that the specified operations are defined, then

(A é B)(C é D) = (AC) é (BD) (A+B) é C = A é C+B é C (A é B)T_{= A}T_{é B}T (A é B)−1_{= A}−1_{é B}−1 A > 0, B > 0 S (A é B) > 0. REFERENCES

1. Scherer, G. (1977). On Energy Estimates for Difference Approximations to Hyperbolic

Partial Differential Equations, Ph.D. thesis, Department of Scientific Computing, Uppsala

2. Kreiss, H. O., and Scherer, G. (1974). Finite Element and Finite Difference Methods for

Hyperbolic Partial Differential Equations, Mathematical Aspects of Finite Elements in

Partial Differential Equations, Academic Press, New York.

3. Olsson, P. (1992). High-Order Difference Methods and Data-Parallel Implementation, Ph.D. thesis, Uppsala University, Department of Scientific Computing.

4. Strand, B. (1996). High-Order Difference Approximations for Hyperbolic Initial Boundary

Value Problems, Ph.D. thesis, Uppsala University, Department of Scientific Computing.

5. Strand, B. (1994). Summation by parts for finite difference approximations for d/dx.

J. Comput. Phys. 110, 47–67.

6. Olsson, P. (1995). Summation by parts, projections, and stability I. Math. Comp. 64, 1035–1065.

7. Olsson, P. (1995). Summation by parts, projections, and stability II. Math. Comp. 64, 1473–1493.

8. Carpenter, M. H., Gottlieb, D., and Abarbanel, S. (1994). The stability of numerical boundary treatments for compact high-order finite-difference schemes. J. Comput. Phys.

108, 272–295.

9. Carpenter, M. H., Gottlieb, D., and Abarbanel, S. (1994). Time-stable boundary condi-tions for finite-difference schemes solving hyperbolic systems: Methodology and applica-tion to high-order compact schemes. J. Comput. Phys. 111, 220–236.

10. Carpenter, M. H., Nordström, J., and Gottlieb, D. (1999). A stable and conservative interface treatment of arbitrary spatial accuracy. J. Comput. Phys. 148, 341–365.

11. Nordström, J., and Carpenter, M. H. (1999). Boundary and interface conditions for high order finite difference methods applied to the Euler and Navier–Stokes equations.

J. Comput. Phys. 148, 621–645.

12. Nordström, J., and Carpenter, M. H. (2001). High order finite difference methods, multi-dimensional linear problems and curvilinear coordinates. J. Comput. Phys. 173, 149–174. 13. Dridi, K. H., Hesthaven, J. S., and Ditkowski, A. (2001). Staircase-free finite-difference

time-domain formulation for general materials in complex domains. IEEE Trans.

Anten-nas Propagat. 49, No. 5.

14. Ditkowski, A., Dridi, K. H., and Hesthaven, J. S. (2001). Convergent cartesian grid methods for maxwell’s equations in complex geometries. J. Comput. Phys. 170, 39–80. 15. Yefet, A., and Petropoulus, P. G. (2001). A staggered fourth-order accurate explicit finite

difference scheme for the time-domain Maxwell’s equations. J. Comput. Phys. 168, 286–315.

16. Yefet, A., and Turkel, E. (2001). Fourth-order accurate compact implicit method for the Maxwell’s equations with discontinuous coefficients. Appl. Numer. Math. 33, 125–134. 17. Turkel, E., and Yefet, A. (2001). On the construction of a high order difference scheme

for complex domains in a cartesian grid. Appl. Numer. Math. 33, 113–124.

18. Gustafsson, B., Kreiss, H.-O., and Oliger, J. (1995). Time Dependent Problems and

Dif-ference Methods, Wiley.

19. Wangsness, R. K. (1986). Electromagnetic Fields, Wiley.