Hybrid Solvers for the Maxwell Equations in Time-Domain

(1)

FREDRIK EDELVIK

Hybrid Solvers for the Maxwell Equations in

Time-Domain

(2)

presented at Uppsala University in 2002 Abstract

Edelvik, F. 2002. Hybrid Solvers for the Maxwell Equations in Time-Domain. Acta Univ. Ups. Uppsala Dissertations from the Faculty of Science and Technology 40. 164 pp. Uppsala. ISBN 91-554-5354-6.

The most commonly used method for the time-domain Maxwell equations is the Finite- Difference Time-Domain method (FDTD). This is an explicit, second-order accurate method, which is used on a staggered Cartesian grid. The main drawback with the FDTD method is its inability to accurately model curved objects and small geometrical features. This is due to the Cartesian grid, which leads to a staircase approximation of the geometry and small details are not resolved at all.

This thesis presents different ways to circumvent this drawback, but still take advantage of the benefits of the FDTD method. An approach to avoid staircasing errors but still retain the efficiency of the FDTD method is to use a hybrid grid. A few layers of unstructured cells are used close to curved objects and a Cartesian grid is used for the rest of the domain. For the choice of solver on the unstructured grid two different alternatives are compared: an explicit Finite-Volume Time-Domain (FVTD) solver and an implicit Finite-Element Time-Domain (FETD) solver.

The hybrid solvers calculate the scattering from complex objects much more efficiently compared to using FDTD on highly resolved Cartesian grids. For the same accuracy in the solution roughly a factor of 10 in memory requirements and a factor of 20 in execution time are gained.

The ability to model features that are small relative to the cell size is often important in electromagnetic simulations. In this thesis a technique to generalize a well- known subcell model for thin wires, in order to take arbitrarily oriented wires in FETD and FDTD into account, is proposed. The method gives considerable modeling flexibility compared to earlier methods and is proven stable. The results show excellent consistency and very good accuracy on different antenna configurations.

The recursive convolution method is often used to model frequency dispersive materials in FDTD. This method is used to enable modeling of such materials in the unstructured FVTD and FETD solvers. The stability of both solvers is analyzed and their accuracy is demonstrated by computing the radar cross section for homogeneous as well as layered spheres with frequency dependent permittivity.

Keywords: Maxwell’s equations, time-domain, finite volume methods, finite element methods, hybrid solver, dispersive materials, thin wires.

Fredrik Edelvik, Department of Information Technology, Scientific Computing, Uppsala University, Box 337, SE-751 05 Uppsala, Sweden

c

° Fredrik Edelvik 2002 ISSN 1104-2516

ISBN 91-554-5354-6

Printed in Sweden by Elanders Gotab, Stockholm 2002

(3)

Acknowledgments

First of all I would like to thank Dr. Gunnar Ledfelt for all our endless discussions on the material in this thesis. I have really enjoyed working together with you!

I thank my advisor Professor Per L¨otstedt for sharing his deep knowledge in numerical analysis. I am grateful to Dr. Ulf Andersson for having the main responsibility of the FDTD code which is the basis for the time-domain solvers in GEMS and for answering my questions regarding the computer system at PDC.

I thank Dr. Bo Strand for leading the GEMS project and for fruitful discussions on frequency dispersive materials. I am grateful to Anders ˚Alund for his hard work on the GEMS solvers and for always finding a compiler where my code does not work.

Furthermore, I would like to express my sincere gratitude to:

Dr. Douglas Riley for promptly answering my questions on computational electromagnetics and for our joint work on thin wires. Thanks also for pushing me into the world of finite elements.

All the people involved in the GEMS and GEMS2 projects. It has been inspiring to work together with you and to understand that there is an industrial need for the results in this thesis. In particular I would like to thank Jonas Gustafsson for making most of the hybrid grids used in the results sections.

Lars Eriksson for implementing parts of the FETD solver.

Martin Nilsson for his work on the fast multipole method. To be able to validate my results against yours has been very valuable.

All the people involved in tdboll activities.

My wife Anna for love and support. I love you!

This work was financed by the Parallel and Scientific Computing Institute (PSCI), the Swedish Agency for Innovation Systems (VINNOVA) and the National Aero- nautical Research Program (NFFP). I have been able to participate in a number of conferences thanks to: PSCI, Bernt J¨armarks stiftelse f¨or vetenskaplig forskn- ing, the European Space Agency, Ericsson’s Research Foundation and C F Lil- jewalchs resestipendiefond.

(4)

(5)

Introduction

1.1 Computational Electromagnetics

The rapid increases in computing capabilities and the development of new numerical methods have facilitated the solutions of larger, more complex electromagnetic scattering and radiation problems. Numerical methods are today at a level where they can be used to accurately simulate the electromagnetic re- sponse from very complicated objects such as aircraft, satellites, cars and various electronic devices. A numerical simulation can often replace expensive and very time consuming measurements especially in the design and construction phase.

Computational Electromagnetics (CEM) can for instance be used to deter- mine the Radar Cross Section (RCS) for objects. The RCS is a far-field parameter, which characterizes the scattering properties of radar targets. For a target, there are monostatic and bistatic RCS. For monostatic RCS the transmitter and the receiver are at the same location, which is not the case for bistatic RCS.

When designing a low-observable target, the RCS is the parameter you attempt to minimize.

Electronic systems are increasing rapidly in numbers and complexity as well as their sensitivity to electromagnetic disturbances. Such a disturbance could for example be direct or indirect effects of lightning. The high surface currents that are induced can damage safety critical electronic systems in the aircraft. Test- ing of Electromagnetic Compatibility (EMC) properties for large and complex objects is an area where CEM is used successfully.

Antenna design is another important field of application for CEM. For wire- less communication systems the antenna is one of the most critical components and a good design of the antenna is crucial for the overall system performance.

Software that can be used to simulate an entire antenna including transmission line accurately is an important tool for engineers.

1

(10)

1.2 Numerical methods

Electromagnetic fields are governed by the Maxwell equations, see Chapter 2.

They can be solved either in the time-domain or the frequency-domain. Fur- thermore, the numerical method can be applied either to the Partial Differential Equation (PDE) formulation of the Maxwell equations or to an integral formulation.

Since, our focus is on time-domain methods we will concentrate on them and only briefly comment upon frequency-domain methods. The advantage of time-domain methods is that we can solve for a broad range of frequencies at one single simulation.

1.2.1 Frequency-domain methods, integral formulation

The most well-known frequency-domain method is the Method of Moments (MoM) [68]. It is based on integral equations describing interactions between surface currents. The integral equations transform the Maxwell equations into the problem of finding currents on the surface of the object. This implies that the dimen- sion of the problem is reduced by one. Another advantage is that several angles of incidence are easily treated through the use of multiple right-hand sides in the resulting system of equations. A drawback is that problems for which the volume to surface ratio is low, which is often the case for EMC problems, and problems involving varying material properties, are not solved efficiently.

MoM results in a dense linear system of equations. Solving this system with Gaussian elimination requires O(N³) arithmetic operations, where N is the number of unknowns. If we keep the number of elements per wavelength constant, N is proportional to f², where f is the frequency. The work to solve the MoM system directly is therefore O(f⁶). One way to diminish this workload is to solve the system using iterative methods. Iterative methods are usually based on matrix-vector multiplication, which has a complexity of O(N²). Hence, the work to solve the MoM system iteratively is O(f⁴) if the iterative method converges rapidly. An even better complexity can be achieved if the communication between the unknown current sources is based on multipole expansions. In this case the linear system of equations can be solved with O(N log(N )) arith- metic operations if a multilevel method is used, which results in a workload of O(f²log(f²)) [12, 13, 49, 61].

1.2.2 Frequency-domain methods, PDE formulation

The most common PDE formulation of the Maxwell equations in the frequency- domain is the vector Helmholtz equation. For the electric field it is given by

∇ × (1

µ∇ × E) = ω²²E , (1.1)

(11)

1.2. NUMERICAL METHODS 3

where ω is the angular frequency, and ² and µ are space-dependent material properties. Helmholtz equation is usually solved with Finite Element Meth- ods (FEM) [66] due to their geometrical flexibility and ability to handle varying material properties. FEM in frequency-domain is efficiently used for calculating the RCS of deep open cavities [43] and it is also commonly used within the microwave technology industry. The widespread commercial code HFSS [32] uses FEM.

1.2.3 Time-domain methods, integral formulation

Time-domain methods for the integral formulation of the Maxwell equations have so far not been widely used. However, in the last few years there has been increased efforts in this area. Most methods are so-called Marching-On- in-Time (MOT) methods. The complexity of traditional MOT methods is O(NtN_s²), where Nt is the number of time steps and Ns is the number of surface unknowns. This complexity can be improved by using so-called Plane-Wave Time-Domain (PWTD) methods [12, 24]. PWTD methods are an extension of the frequency-domain fast multipole method to the time-domain. The complex- ity of the multilevel PWTD is O(NtNslog(Ns)).

Integral equation methods in time-domain have some advantages compared to PDE methods in the time-domain: they do not suffer from dispersion errors, they only discretize a surface and no absorbing boundary condition is needed. In fact, recently time-domain integral equation methods have been used to create absorbing boundary conditions for FDTD and FETD [25, 36].

The drawback with MOT methods is that they have been shown prone to late time instabilities. This issue has been studied in detail by Walker’s group and they demonstrate in [9, 17] that MOT schemes for solving magnetic field integral equations can be stabilized for “all practical purposes” by relying on accurate spatial integration rules and implicit timestepping.

1.2.4 Time-domain methods, PDE formulation

The most commonly used method for the time-domain Maxwell equations is the Finite-Difference Time-Domain method (FDTD) [72]. It has several advantages compared to other time-domain methods. The simple data structures make it straightforward to implement and facilitate vectorization and parallelization.

The main drawback with the FDTD method is its inability to accurately model curved objects and small geometrical features. This is due to the Cartesian grid, which leads to a staircase approximation of the geometry and small details are not resolved at all. In many applications, e.g. RCS predictions for complex objects, the staircased approximations of the objects are not accurate enough.

There are several suggested remedies to circumvent the effects of the errors introduced by the staircase approximations. These include Cartesian subgridding [50, 64], conformal modeling [15], staircase free Cartesian grid meth-

(12)

ods [16, 18], overlapping grids [73, 74, 75, 76], multiblock body fitted structured grids [2, 60] and unstructured grids [31, 41].

Cartesian subgridding suffers from stability problems and the interpolations required at grid interfaces violate the divergence relations. The idea of the conformal modeling approach is to change the updating stencils in all cells that are intersected by the boundary. The main drawback with this method is its extreme complexity in 3D. Furthermore, the question of stability is not sufficiently inves- tigated. The staircase free Cartesian grid methods are also based on a change of the FDTD updating stencil close to metallic or dielectric boundaries. However, stability in 3D remains to be solved.

The overlapping grids approach has serious stability problems due to the interpolations between the grids [74, 76]. Another disadvantage is the grid generation which often requires a high user intervention. The methods based on multiblock body fitted structured grids originate from CFD. They are designed to be able to capture shocks and are in general very dissipative, which is not a good property for a solver of Maxwell’s equations. Moreover, the generation of the grid is a major hassle for complex objects and could take several man months even for a skilled engineer. Using an unstructured grid for the whole domain makes the grid generation considerably easier for complex objects. But an unstructured solver generally implies a more complex data structure and more arithmetic operations per cell compared to a structured solver. However, the higher-order discontinuous Galerkin method presented in [31] has shown promising results.

Another approach to avoid staircasing that has gained popularity lately is to use unstructured grids near curved objects and around small geometrical details, but revert to Cartesian grids as quickly as possible for the rest of the computational domain [20, 21, 55, 57]. This combines the efficiency of Cartesian grids with the flexibility of unstructured grids. Another advantage of using a hybrid grid is that it gives access to accurate absorbing boundary conditions developed for Cartesian grids [6, 26]. However, recently a Perfectly Matched Layer (PML) absorbing boundary condition has been developed also for the finite element method in time-domain [35].

1.3 The GEMS and GEMS2 projects

All the research in this thesis has been performed in close cooperation with Swedish industry in the Parallel and Scientific Computing Institute (PSCI) projects GEMS and GEMS2. GEMS is an abbreviation for General Electro- Magnetic Solvers and was a Swedish three-year code development project that was supported by an extensive research program [62]. A substantial part of the funding was supplied by the National Aeronautical Research Program (NFFP).

The GEMS project has now been superceded by the GEMS2 project.

The main objective of the GEMS project was to develop a software suite

(13)

1.4. OUTLINE AND MAIN RESULTS 5

for solving the Maxwell equations. This software suite aims to be state-of-the- art and to form a platform for future development by Swedish industry and academia. The codes are currently used in an industrial environment.

The core of the software suite is two hybrid codes, one for the time-domain and one for the frequency-domain. The time-domain code is a hybrid between FDTD, an explicit Finite-Volume Time-Domain (FVTD) solver and an implicit Finite-Element Time-Domain (FETD) solver. The frequency-domain code is a hybrid between MoM, Physical Optics (PO) and Geometrical Theory of Diffrac- tion (GTD) [7, 22, 23].

The industrial partners in GEMS were Ericsson Microwave Systems, Saab Ericsson Space and Saab Avionics. Code developers were PSCI, the Swedish Institute of Applied Mathematics (ITM) and the Swedish Defence Research Agency (FOI). The GEMS2 project involves PSCI, Ericsson Microwave Systems, Saab Avionics and ITM.

1.4 Outline and main results

The main part of this thesis is devoted to the unstructured FVTD and FETD solvers. The description and analysis of the FVTD solvers were the basis of my Licentiate thesis presented in October 2000 [19]. In addition to the material in the Licentiate thesis this thesis includes the FETD solver and subcell modeling of thin wires. The coupling of the solvers has been improved and the chapter on frequency dispersive materials now includes FETD and a stability analysis.

Chapters 2, Governing equations

Electromagnetic fields are governed by the Maxwell equations which are stated in Chapter 2.

Chapters 3 and 4, The FVTD solvers

In Chapters 3 and 4 the FVTD solvers in 2D and 3D are presented and analyzed. The 2D FVTD solver is in space very similar to the solver used by Riley and Turner [55], whereas the solver in 3D differs in several ways. The main improvements of the 3D solver compared to [55] include the generation of the dual grid, the update of the node values and the treatment of material interfaces. The time integration method we use is a third-order staggered Adams–Bashforth scheme (ABS3) proposed by Ghrist et al. [27] for scalar wave equations. Stability analysis of this scheme shows that it works very well in our case. Furthermore, using ABS3 makes the time-coupling between FDTD and the FVTD solver straightforward. A spatial filter of Laplace type is introduced in 3D to improve the stability characteristics of the solver.

(14)

Chapter 5, The FETD solver

Chapter 5 is devoted to the FETD solver. The solver is based on a discretization of the second-order wave equation for the electric field using edge-elements as basis functions [37, 47]. A stability proof is given which shows that the FETD solver is unconditionally stable when the Newmark-Beta scheme is used for timestepping.

Chapter 6, Hybridization

In Chapter 6 the coupling between the unstructured solvers and FDTD is presented. Coupling of two stable methods does not necessarily result in a stable hybrid method and in fact, the hybridization technique described in Section 6.2 is not stable. The source of the instabilities is the diagonal component that is not present in FDTD and therefore is calculated by interpolation. This breaks the symmetry of the hybrid method and makes it unstable.

To develop a symmetric interface for the coupling of FDTD and FETD we follow [57] and use a two-cell thick layer of hexahedral elements, such that cells neighboring the tetrahedral region are subdivided into pyramidal and tetrahedral elements. The stability of the coupling follows from the fact that the resulting hybrid solver can be viewed as an Explicit-Implicit algorithm. Such algorithms are proven stable in [34].

Chapter 7, Numerical experiments in 2D

In Chapter 7 we show that the staircasing errors in FDTD destroy the second- order accuracy on metallic, dielectric and diamagnetic cylinders in 2D. Further- more, we show that second-order accuracy is retained by our hybrid solvers with a correct averaging of the permeability and permittivity parameters. Finally, we perform an experiment on a PMC wall, which is a classical case demonstrating the staircasing errors in FDTD [10].

Chapter 8, A comparison of the hybrid solvers

In this chapter we compare the two hybrid solvers in terms of accuracy, efficiency and stability. Results from several scattering cases that show the superiority of the hybrid solvers compared to stand-alone FDTD are presented. The conclusion from the comparison is that the hybrid solver combining FETD and FDTD is the preferred choice. The major reason is that it is stable as long as the time step satisfies the CFL condition in the FDTD region, whereas the solver combining FVTD and FDTD may exhibit late time instabilities.

(15)

1.4. OUTLINE AND MAIN RESULTS 7

Chapter 9, Frequency dispersive materials

The ability to model frequency dispersive materials is an important feature for a general purpose electromagnetic code. They are for instance used as radar absorbing materials to minimize the RCS for an aircraft. One approach to model these materials in FDTD is to use the Recursive Convolution (RC) method [44].

In Chapter 9 the piecewise constant centered RC method presented in [59] is generalized to our FVTD and FETD solvers. A RC method for FETD has been proposed in [77]. However, our method includes fundamental improvements because it preserves the definiteness properties of the matrices as well as reducing to FDTD on hexahedral elements. Both these issues are important for the stability of the FDTD-FETD solver.

Chapters 10 and 11, Modeling thin wires in FETD and FDTD

The ability to model features that are small relative to the cell size is often important in electromagnetic simulations. Thus, the development of accurate models that characterize the physics of the features without the need for a highly resolved grid is essential.

Thin wires are often important parts of electromagnetic compatibility and antenna problems. A subcell model for thin wires in the FDTD method using modified telegraphers equations has been developed [33]. A generalization of this model to arbitrarily oriented wires has been proposed by Ledfelt [40]. Riley has developed a method that incorporates thin wires into the FETD method by discretizing a second-order wave equation for the current [54].

In the FETD chapter we present a new method, which is also based on a second-order equation for the current derived from modified telegraphers equations. A discretization of this equation with linear nodal basis functions and the use of a radial weighting function result in a symmetric spatial coupling between field and wire. We prove using the energy method that the fully discrete field-wire system is unconditionally stable. Furthermore, our method allows for arbitrarily located and oriented wires with respect to the tetrahedral grid, which give considerable modeling flexibility. Results on different configurations of dipole and loop antennas show good accuracy and excellent consistency. Note that by consistency in this context we mean that the results are independent of the location and orientation of the antenna within the grid.

In the FDTD chapter we present a similar method for modeling arbitrarily oriented wires in FDTD. The method is shown to be stable under a CFL limit.

The consistency of the results are also in this case excellent, and the accuracy is very similar to the FETD results.

(16)

Chapter 12, Complex applications in 3D

In this chapter results from a few very challenging scattering cases are presented. The results obtained by the FDTD-FETD hybrid solver are compared to measurements, frequency-domain solutions and highly resolved FDTD solutions. Two EMC examples are also included namely a visualization of a shielding enclosure and a simulation of a lightning strike to a Saab 2000 aircraft.

1.5 List of papers

This thesis in partly based on material from the following papers:

I F. Edelvik, Analysis of a Finite Volume Solver for Maxwell’s Equations, In R. Vilsemeier (editor), Finite Volumes for Complex Applications II, IVG, University of Duisburg, Germany, pp. 141-148, July 1999.^∗

II E. Abenius, U. Andersson, F. Edelvik, L. Eriksson, G. Ledfelt, Hybrid Time Domain Solvers for the Maxwell Equations in 2D, International Jour- nal for Numerical Methods in Engineering, 53(9):2185-2199, March 2002.^∗∗

III F. Edelvik, G. Ledfelt, Explicit Hybrid Time Domain Solver for the Maxwell Equations in 3D, Journal of Scientific Computing, 15(1):61-78, March 2000.^{∗ ∗ ∗}

IV F. Edelvik, B. Strand, Frequency Dispersive Materials for 3D Hybrid Solvers in Time Domain, to appear in IEEE Transaction on Antennas and Prop- agation, 2002.^{∗ ∗ ∗∗}

V F. Edelvik, G. Ledfelt, A comparison of time-domain hybrid solvers for complex scattering problems, International Journal of Numerical Modeling, 15(5), Sept/Oct 2002.^{∗ ∗ ∗}

VI F. Edelvik, G. Ledfelt, P. L¨otstedt, D. J. Riley, An Unconditionally Stable Subcell Model for Arbitrarily Oriented Thin Wires in the FETD Method, Technical Report 2002-007, Department of Information Technology, Up- psala University, Sweden. Submitted to IEEE Transaction on Antennas and Propagation.^{∗ ∗ ∗∗}

VII F. Edelvik, A New Technique for Accurate and Stable Modeling of Arbitrar- ily Oriented Thin Wires in the FDTD Method, Technical Report 2002-016, Department of Information Technology, Uppsala University, Sweden. Sub- mitted to IEEE Transaction on Electromagnetic Compatibility.^{∗ ∗ ∗∗}

∗Copyright c° 1999 Hermes, Paris, France

∗∗Copyright c° 2000 Plenum Publishing Corporation

∗ ∗ ∗Copyright c° 2002 John Wiley & Sons, Ltd

∗ ∗ ∗∗Copyright c° 2002 IEEE

(17)

Chapter 2

Governing equations

2.1 The Maxwell equations

Electromagnetic fields are governed by the Maxwell equations. In general form they are given by

∇ · D = ρf, (2.1)

∂B

∂t + ∇ × E = 0 , (2.2)

∇ · B = 0 , (2.3)

∂D

∂t − ∇ × H + Jf = 0 , (2.4) where D is the electric flux density vector, E is the electric field, B is the magnetic flux density vector, H is the magnetic field, ρf is the free charge, and Jf is the total free current density. Equation (2.1) is Gauss’ law for the electric field and it summarizes Coulomb’s law of force between point charges plus the electrical effects of matter, while (2.2) represents Faraday’s law of induction. The third equation reflects that free magnetic charges have not yet been observed.

Finally, (2.4) includes Amp`ere’s law of force between currents plus the magnetic effects of matter and conservation of free charge.

For linear, isotropic and non-dispersive media we have

D = ²E , B = µH , Jf = σE + J , (2.5)

where ² is the electric permittivity, µ is the magnetic permeability, σ is the electric conductivity and J is the free current density arising from sources other

9

(18)

than conductivity. Substituting these relations into (2.1–2.4) implies

∇ · ² E = ρf, (2.6)

µ∂H

∂t + ∇ × E = 0 , (2.7)

∇ · µ H = 0 , (2.8)

²∂E

∂t − ∇ × H + σ E = −J . (2.9) Hence, the Maxwell equations in E and H are a system of eight scalar equa- tions in six variables. At a first glance it may seem that the system is overdeter- mined, but it is readily shown that if the two divergence relations are satisfied initially they are always satisfied. Together with boundary conditions and ini- tial values the Maxwell equations constitute a well posed system of hyperbolic partial differential equations.

If we take the time derivative of (2.9) and use (2.7) we obtain the vector wave equation

²∂²E

∂t² + ∇ × µ1

µ∇ × E

¶

+ σ∂E

∂t = −∂J

∂t . (2.10)

For lossless homogeneous materials without sources this reduces to

∂²E

∂t² − c²∆E + σ

²

∂E

∂t = 0 , (2.11)

where c = 1/√

µ² is the speed of propagation for the electromagnetic wave. In a similar manner, we can show that

∂²H

∂t² − c²∆H +σ

²

∂H

∂t = 0 , (2.12)

for lossless homogeneous materials without sources.

2.2 Integral formulation

Using Gauss’ divergence theorem and Stokes’ law the integral formulation of Maxwell’s equations is obtained

I

A

² E · n dA = ρf, (2.13)

∂

∂t Z

A

µ H · n dA + I

Γ

E · dl = 0 , (2.14) I

A

µ H · n dA = 0 , (2.15)

∂

∂t Z

A

² E · n dA − I

Γ

H · dl + Z

A

σ E · n dA = − Z

A

J · n dA , (2.16)

(19)

2.3. REDUCTION TO TWO DIMENSIONS 11

where A is an arbitrary area with unit normal n and Γ is the contour that encloses A.

2.3 Reduction to two dimensions

If we have no variation in the z-direction all derivatives with respect to z equal zero. In that case the Maxwell curl equations in Cartesian coordinates reduce to

µ∂Hx

∂t = −∂Ez

∂y , (2.17)

µ∂Hy

∂t = ∂Ez

∂x , (2.18)

²∂Ez

∂t = ∂Hy

∂x −∂Hx

∂y − σ Ez− Jz, (2.19)

²∂Ex

∂t = ∂Hz

∂y − σEx− Jx, (2.20)

²∂Ey

∂t = −∂Hz

∂x − σEy− Jy, (2.21)

µ∂Hz

∂t = ∂Ex

∂y −∂Ey

∂x , (2.22)

where the first three equations constitute the 2D Transverse Magnetic (TM) Maxwell’s equations and the last three constitute the 2D Transverse Electric (TE) Maxwell’s equations. The TM and TE modes are decoupled, i.e. they contain no common field component. The two modes are completely independent for isotropic materials and they can exist simultaneously with no mutual interac- tion.

2.4 Material properties

In the Maxwell equations we have three parameters. For vacuum they are µ0= 4π · 10⁻⁷, ²0≈ 8.8541878 · 10⁻¹² and σ = 0. For other materials the permeability and permittivity are defined relative to those of vacuum, i.e. we have ² = ²r²0

and µ = µrµ0. For most materials, ²r and µr are frequency dependent, see Chapter 9. Materials for which we assume that ²r and µr are independent of frequency are referred to as simple materials. Materials for which µ_r 6= 1 are metals with high conductivity and we treat them as perfect electric conductors.

(20)

2.5 Boundary conditions

At the interface between two lossless media we have [11]

n · (D1− D2) = 0 , n × (E1− E2) = 0 , n · (B1− B2) = 0 , n × (H1− H2) = 0 ,

(2.23)

where the subscripts indicate which region the field belongs to, and n is the interface normal. Using the relations in (2.5), we get

n · (²1E1− ²2E2) = 0 , n × (E1− E2) = 0 , n · (µ1H1− µ2H2) = 0 , n × (H1− H2) = 0 .

(2.24)

For Perfect Electric Conductors (PEC) we have [11]

n · ²E = ρs, n × E = 0 , n · H = 0 , n × H = Js,

(2.25)

where ρsis the surface charge density and Jsis the surface current density. The normal n is pointing out of the PEC region.

The tangential electric field is zero at the surface of a PEC. This is a conse- quence of the term perfect conductor. A nonzero tangential electric field would drive an infinite surface current which is clearly unphysical. However, this does not imply that the surface current must be zero. In fact, an external field always implies nonzero surface current since the magnetic field only has tangential components at the PEC surface and the surface current is related to the tangential magnetic field through the fourth condition in (2.25).

(21)

Chapter 3

Finite Volume method in 2D

The Finite-Volume Time-Domain (FVTD) method is based on the following integral formulation of Faraday’s and Amp`ere’s laws:

∂

∂t Z

A

µ H dA = − I

Γ

n × E dl , (3.1)

∂

∂t Z

A

² E dA = I

Γ

n × H dl − Z

A

σE dA , (3.2)

where A is an arbitrary area with unit normal n and Γ is the path that encloses A. We have chosen to focus on the TM Maxwell’s equations (2.17–2.19). Due to the duality of the Maxwell equations a solver for the TE equations (2.20–2.22) could be obtained through a change of variables.

3.1 Spatial discretization

The integral formulations (3.1) and (3.2) are discretized on a staggered grid by introducing a dual grid to the primary triangular grid. The magnetic components are situated at the nodes of the primary grid and the electric components, in the 2D TM case only Ez, are situated at the nodes of the dual grid. The dual grid is created at the preprocessing stage by defining dual nodes at the barycenters of the primary cells; see [19] for a detailed description.

In the 2D TM case integrating the magnetic and electric fields over each dual

13

(22)

and primary cell gives the following integral form:

∂

∂t Z

A^d_j

˜

µ^d_jH dA = −X

k

Z

Γ^d_j,k

n^d_j,k× Ezdl , (3.3)

∂

∂t Z

A^p_i

²^p_iEzdA = X

m

Z

Γ^p_i,m

n^p_i,m× H dl − Z

A^p_i

σ^p_i EzdA , (3.4)

where A^d_j is the area of dual cell j, Γ^d_j is the path that encloses A^d_j and n^d_j,kare the unit edge normals for the dual edges k in dual cell j. The variables belonging to the primary cell i are defined similarly.

All materials are defined relative to primary grid cells. For dielectric ma- terials that is a natural definition. However, the magnetic permeability, µ, is associated with the dual grid cells. Therefore, µ requires averaging. The aver- age permeability of dual cell j is computed as

˜ µ^d_j =

P

q

µ^p_qA^d_j,q

A^d_j , (3.5)

where A^d_j,q is the part of the area A^d_j that is inside primary cell q. Performing the material averaging in this manner preserves the second-order accuracy of the solver for inhomogeneous materials, see Section 7.2.

The area integrals in (3.3) and (3.4) are evaluated by taking the average values of the fields multiplied by the areas of the respective cells. Simplifying the two integrands in the line integrals implies

˜ µ^d_j ∂

∂tHj= 1 A^d_j

X

k

Z

Γ^d_j,k

E˜z|j,kt^d_j,kdl , (3.6)

²^p_i ∂

∂tEz|i+ σ_i^pEz|i= 1 A^p_i

X

m

Z

Γ^p_i,m

H · t˜ ^p_i,mdl , (3.7)

where t^d_j,k are unit vectors for the dual edges k in dual cell j and t^p_i,m are unit vectors for the edges m in primary cell i. The line integral in (3.6) is evaluated by assuming that the electric field is piecewise linear along the dual edges. Hence, E˜z|j,k is computed by taking the arithmetic mean value of the electric field at the two nodes defining the dual edge, t^d_j,k. However, we cannot use the same approach when calculating the integral in (3.7) since that does not guarantee that the divergence is preserved on a local cell level. This has been found to be very critical if spurious modes in the numerical solution are to be suppressed.

The divergence is preserved if we incorporate an “FDTD”-correction along the edges in the primary grid (see Section 3.3 for a proof). Following Riley et al. [55],

(23)

3.1. SPATIAL DISCRETIZATION 15

the magnetic field projected along the primary edge t^p_i,m is evaluated as H · t˜ ^p_i,m= H · n^d_j,k¡

n^d_j,k· t^p_i,m¢ +1

2

£(Hj+ Hr) −¡

(Hj+ Hr) · n^d_j,k¢ n^d_j,k¤

· t^p_i,m, (3.8) where Hj and Hr are the magnetic field at the two nodes defining the primary edge and H · n^d_j,k is the “FDTD”-correction in the direction n^d_j,k orthogonal to the dual edge t^d_j,k, which crosses the primary edge t^p_i,m, see Figure 3.1. The

“FDTD”-correction is updated according to

∂ H · n^d_j,k

∂t =E_z|_q− E_z|_i

∆^d_j,kµ˜^d_j,k , (3.9) where ∆^d_j,k is the length and ˜µ^d_j,k is the average permeability of the dual edge t^d_j,k. The average permeability, ˜µ^d_j,k, is approximated as

˜

µ^d_j,k= µ^p_q∆q+ µ^p_i

³

∆^d_j,k− ∆q

´

∆^d_j,k , (3.10)

where i and q are the two primary cells sharing dual edge t^d_j,k, ∆q is the part of the length of the dual edge that is inside primary cell q, see Figure 3.1.

n_j,k^d t_i,m^p

H_j H_r

E_{z i}

Ez q

∆^q t_j,k^d

Figure 3.1: The magnetic field along t^p_i,m is approximated by H_j, H_r and an

“FDTD”-correction in direction n^d_j,k.

Taking a closer look at (3.8) we see that if the primary and dual edges are orthogonal, the vectors n^d_j,k and t^p_i,m are parallel and the second part of (3.8)

(24)

vanishes. Hence, the name “FDTD”-correction is somewhat misleading since that term is actually the important one. The magnetic node values are only used to give a better approximation of the edge-projected field in the case when n^d_j,k and t^p_i,m do not align.

The boundary condition for a Perfect Magnetic Conductor (PMC) gives us that the tangential component of the magnetic field at the object is zero. A complication occurs whenever the computation of ˜H · t^p_i,m does not reduce to H · n^d_j,k, where t^p_i,m denotes a primary edge with one node on the boundary.

When that is not the case the following alternative is used (see Figure 3.2):

H · t˜ ^p_i,m = H · n^d_j,k¡

n^d_j,k· t^p_i,m¢ +1

2

£Hj+ (Hj· ñr) ñr− (Hj+ (Hj· ñr) ñr) · n^d_j,k¤

· t^p_i,m, (3.11)

where ˜nr denotes an average normal at the boundary node. This normal is defined by taking the average direction of the two boundary edges including the boundary node and then taking the cross product with that direction and ˆz.

n_j,k^d n~_r

p

H t

r i,m

Hj

Ez i

PMC boundary

Figure 3.2: Primary and dual cells at a PMC boundary.

For a Perfect Electric Conductor (PEC) the tangential electric components, in the TM case only Ez, should equal zero at the boundary. However, the Ezcomponents are not situated on the boundary. The boundary condition is enforced by setting the electric field below the conductor equal to the value of the electric field directly above the conductor with a change of sign, see Figure 3.3.

Hybrid Solvers for the Maxwell Equations in Time-Domain

FREDRIK EDELVIK