for the Computation of
Electromagnetic Fields
for the Computation of
Electromagnetic Fields
in Strongly Inhomogeneous Media
Proefschrift
ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,
op gezag van de Rector Magnificus prof.dr.ir. J.T. Fokkema, voorzitter van het College van Promoties,
in het openbaar te verdedigen op maandag 7 februari 2005 om 13:00 uur
door
Pieter Jorna
wiskundig ingenieur geboren te Doetinchem.
Prof.dr.ir. P.M. van den Berg
Samenstelling promotiecommissie:
Rector Magnificus voorzitter
Prof.dr.ir. P.M. van den Berg Technische Universiteit Delft, promotor Prof.dr.ir. P. Wesseling Technische Universiteit Delft
Prof.dr.ir. A.T. de Hoop Technische Universiteit Delft Prof.dr.ir. P. Zwamborn Technische Universiteit Eindhoven Prof.dr.ir. F. Olyslager Universiteit Gent
Dr.ir. E.C. Slob Technische Universiteit Delft Dr.ir. I.E. Lager Technische Universiteit Delft
ISBN 90-9019073-2 NUR 950
Cover design by: Dennie Rensink & Pieter Jorna
“A mathematician’s 3D perspective on Pink Floyd’s ‘Dark Side of the Moon’ ” Copyright c 2005 by P. Jorna
All rights reserved. No part of this thesis may be reproduced, stored in a retrieval system, or transmitted, in any form or by any means, electronic, mechanical, photocopying, recording, or otherwise, without the prior written permission from the copyright owner.
1 Introduction 1
1.1 Finite-Element Modelling . . . 2
1.2 Domain-Integrated Field Equations Approach . . . 5
1.3 Outline of the Thesis . . . 7
2 The Electromagnetic Field Equations 11 2.1 Time-Domain Maxwell’s Equations in Differential Form . . . 12
2.2 Time-Domain Constitutive Relations . . . 13
2.3 Frequency-Domain Representations . . . 15
2.4 Frequency-Domain Maxwell’s Equations . . . 16
2.5 Frequency-Domain Constitutive Relations . . . 17
2.6 Surface-Integrated Field Equations . . . 18
2.7 Volume-Integrated Field Equations . . . 19
2.8 Interface Conditions . . . 20
2.9 Boundary Conditions . . . 23
3 Discretization of the 2D Domain of Computation 29 3.1 The 2D Electromagnetic Field . . . 29
3.2 Geometrical Discretization . . . 31 3.3 Local Discretization of the Electromagnetic Field Quantities . 36
3.4 Global Discretization of the Electromagnetic Field Quantities 46
3.5 Discretization of the Material Related Quantities . . . 59
3.6 Discretization of the Source Densities . . . 61
4 Discretization of the 2D Electromagnetic Field Equations 63 4.1 Geometrical Elements of the Prism . . . 63
4.2 Volume-Integrated Field Equations Method . . . 66
4.3 Surface-Integrated Field Equations Method . . . 73
4.4 Constitutive Relations . . . 77
4.5 Conclusions . . . 80
5 Implementation and Numerical Results of the 2D Methods 83 5.1 System of Equations . . . 83
5.2 Implementation . . . 88
5.3 Geometrical Discretization . . . 89
5.4 Building the System’s Matrix . . . 90
5.5 Solving the System of Linear Equations . . . 93
5.6 The Four-Media Test Problem . . . 100
6 Discretization of the 3D Domain of Computation 111 6.1 3D Geometrical Discretization . . . 111
6.2 Local Discretization of the Electromagnetic Field Quantities . 116 6.3 Global Discretization of the Electromagnetic Field Quantities 124 6.4 Discretization of the Material Related Quantities . . . 135
6.5 Discretization of the Source Densities . . . 137
7 Discretization of the 3D Electromagnetic Field Equations 139 7.1 Volume-Integrated Field Equations Method . . . 139
7.2 Surface-Integrated Field Equations Method . . . 144
7.3 Constitutive Relations . . . 147
7.4 Conclusions . . . 150
8 Implementation and Numerical Results of the 3D Methods 151 8.1 System of Equations . . . 151
8.2 Implementation . . . 155
8.3 Geometrical Discretization . . . 156
8.4 The Four-Media Test Problem . . . 158
8.5 Testing the PMLs . . . 168
8.6 Scattering Computation Using PMLs . . . 176
9 Discussion and Conclusions 183
Bibliography 189
Summary 195
Samenvatting 199
Curriculum Vitae 203
Introduction
Electromagnetic theory is an indispensable part of science and engineering and is inextricably bound up with present-day technology. Some examples of technologies that rely heavily on electromagnetic theory and of which the importance is growing each day are radar, remote sensing, geoelectro-magnetics, bioelectrogeoelectro-magnetics, antennas, wireless communication, optics, high-frequency circuits, and so on. In many of these applications it is very important to accurately predict the electromagnetic field. This widely ap-plicable and elaborate electromagnetic theory is in essence fully described by a set of partial differential equations, the Maxwell equations. Hence, calcu-lating the electromagnetic field is actually the problem of solving Maxwell’s equations subject to given boundary conditions. The range of problems for which analytical solutions can be obtained is unfortunately very small and to obtain solutions for problems of practical interest we usually turn to nu-merical methods.
In this thesis we concentrate on time-harmonic electromagnetic field problems. The Finite-Element (FE) method is a numerical method that is often considered to be very suitable for solving this type of problem in complex configurations consisting of inhomogeneous media. In the following section we shortly describe some of the difficulties encountered when applying the FE method and discuss developments in FE research that try to eliminate these obstacles. The domain-integrated field equations method developed by De Hoop and Lager [19], [20], to solve static and stationary magnetic field
problems for the class of strongly heterogeneous media aims to be as close to the physics as possible and attempts to overcome the problems related to application of the FE method. In this thesis we present an extension of the method to the time-harmonic problem where the electric and magnetic field strengths and the electric and magnetic flux densities are ingeniously coupled via Maxwell’s equations and the constitutive relations.
1.1
Finite-Element Modelling
The Finite-Element (FE) method is a local discretization method with re-spect to both the geometry and physical quantities that serves to generate numerical values of the field quantities characterizing some desired physical state. Historically, the application of FE techniques to structural mechanics has preceded the use of the method in electrical engineering. Thus some of the standard works focus on structural problems, such as Norrie and De Vries [43] and Zienkiewicz [57]. Through the years, the latter book has developed into a three-volume standard text on finite elements: Zienkiewicz and Tay-lor [58], [59] and [60]. Excellent texts on the FE method in computational electromagnetics are presented by Jin [24] and Silvester and Ferrari [47].
The FE computation of electromagnetic fields is often carried out in terms of (vector) potentials. In most applications, however, we are interested in the electric and/or magnetic field strength which follow from the (vector) potential by means of a numerical differentation. This differentiation causes a loss of accuracy of the order of the mesh size, and hence a FE approach based on direct computation of the field strength is advantageous. Methods that simultaneously compute both the electric and magnetic field strength are discussed by Mur [34], [39], but most of the time-harmonic field prob-lems are solved using a FE formulation in terms of either the electric field strength or the magnetic field strength. This implies that we need numerical differentiations to obtain the magnetic field strength in case calculations are performed in terms of the electric field strength and vice versa. This is a serious drawback if we are interested in an accurate solution of both field strengths because, as mentioned before, numerical differentiations cause a loss of accuracy of the order of the mesh size.
When applying the FE method it is essential to use a set of expansion functions that can accurately and consistently approximate the quantity to be computed. While the tangential components of the electric and magnetic field
strength are always continuous, the normal components may show a jump discontinuity across interfaces between different media. Cartesian elements, also referred to as nodal elements, cannot accommodate this behaviour and hence they yield unacceptably large representation errors as is demonstrated by Mur [33]. By assigning degrees of freedom to the edges rather than to the nodes, one creates elements that automatically account for the continuity of the tangential field component across interfaces and that allow for a jump in the normal component. These elements are referred to as edge elements. This type of elements was already described in 1957 by Whitney [55], but their use and importance in electromagnetics were realized only after the fundamental paper by N´ed´elec [42] and the application of the elements by Bossavit and Verit´e [7], [8].
Local expansions using the standard edge elements introduced by N´ed´elec, also referred to as Whitney elements, carry only six unknowns, one on each edge, and as a consequence the field component in the direction of these edges is constant along each of the edges. The low degree of approxima-tion yields large local approximaapproxima-tion errors; Bandelier and Rioux-Damidau [2], Mur and De Hoop [32], and Trabelsi et al. [49] give experimental nu-merical verification of the fact that correspondingly large errors are found in global solutions. Edge elements that do not suffer from large local approxi-mation errors are consistently linear edge elements. These edge elements are linear vector functions of the position in all spatial directions and were first introduced by Mur and De Hoop [32].
Expansions over a tetrahedron that use consistently linear vectorial elements carry twelve unknowns: three at each of its four vertices. Compu-tationally, however, edge elements are less efficient than nodal elements due to the fact that the span of global edge expansion functions is smaller than that of global nodal expansion functions, which causes the edge elements to generate more unknowns for the same mesh. In regions where the media are homogeneous or where the differences in electromagnetic properties are small, the application of edge elements is computationally inefficient, and what is more, discontinuities in the normal component of the field strength appear where this component should be continuous. Methods that alternately em-ploy edge elements and Cartesian elements, depending on the difference in material properties of adjacent subdomains, are described by Mur [33], [34], [35], [36], by Lager [27], and by Lager and Mur [28].
A third type of finite element that is encountered in electromagnetic field computations is the face (or facet) element. By assigning degrees of
freedom to the faces, elements are created that automatically account for the continuity of the normal field component across interfaces and that allow for a jump in the tangential component. Bearing in mind these properties, it is clear that face elements are especially suited to model electric and magnetic fluxes. The span of global face expansion functions is restricted to at most two adjacent elementary subdomains and this explains why consistently lin-ear face elements are computationally more expensive than the consistently linear edge elements previously described. The relative unpopularity of face elements is to a large extent due to their inefficiency.
There is no doubt that FE direct field computation has significantly improved with the introduction of these new types of elements. For a short while it was believed that these new elements were the answer to the feared and unwanted spurious solutions, but Mur showed in [37] and [40] that these elements by themselves do not rule out the existence of spurious solutions. In addition he pointed out several disadvantages of applying edge elements. In another article by Mur [38] it is made clear that to eliminate spurious solutions it is imperative that numerical solutions not only satisfy the dis-cretized Maxwell equations but also all existing compatibility relations. Edge elements provide the proper continuity conditions with respect to the tangen-tial component of the field strength, but other compatibility relations, such as continuity of the corresponding normal flux, are not guaranteed. Simi-larly, face elements provide the proper continuity conditions with respect to the normal component of the flux but other compatibility relations, such as tangential continuity of the corresponding field strength, are not guaranteed. To obtain reliable results, one need consider alternatives to the naive use of either edge or face elements.
An alternative, recognized by Mur in his paper [40], is presented by so-called mixed formulations. In these formulations, also referred to as dual or complementary formulations, complementary vectorial field quantities are chosen that together allow a consistent representation of the electromagnetic field inside the domain of computation. Although they require a larger num-ber of degrees of freedom, dual approaches have the advantage of satisfacto-rily modelling both the field equations and the relevant compatibility rela-tions. The domain-integrated field equations method developed by De Hoop and Lager [19], [20] for static and stationary magnetic field computation in strongly heterogeneous configurations is an example of a dual approach with the unusual property that it does not rely on a variational formulation. In the next section we describe in more detail the class of strongly
heteroge-neous media and present the domain-integrated field equations approach in its original formulation.
1.2
Domain-Integrated Field Equations Approach
The class of strongly heterogeneous media and the domain-integrated field equations method, developed for static and stationary magnetic field compu-tation in configurations containing these media, are best described by quoting De Hoop and Lager from their paper [20], where they give a physical justifi-cation of the method:
“The domain-integrated field equations method can be viewed upon as to be the computational method that models, as closely as possible, the un-derlying physics in strongly heterogeneous (and anisotropic) media. Strongly heterogeneous media are characterized by the property that on any scale on which the macroscopic laws of nature maintain to hold, the granular struc-ture persists. Now, in any practical situation one is only interested in the values of the field quantities in such a medium down to the scale where actual measurements are still feasible. Let us denote this scale as the mesoscopic scale (to distinguish it from the microscopic or atomic scale where for the field description the quantum laws of physics are required). In view of the fact that on the mesoscopic scale heterogeneity still persists, with its as-sociated discontinuous behavior of (at least some components of) the field quantities involved, the field quantities as such do not belong to the class of spatially differentiable functions and, hence, corresponding field equations in differential form cannot serve as a point of departure. On the mesoscopic scale, the domain-integrated field equations are the alternatives to serve as the basis for field computation. In the latter, only the field quantities them-selves occur, whose values should on physical grounds remain bounded and hence integrable over any bounding surface of a computational cell of the discretized geometry. A remarkable feature of the domain-integrated field equations is further that in them only the field components occur that are continuous across an interface of discontinuity in material properties. (For the quasi-static magnetic field these are the tangential components of the magnetic field strength and the normal component of the magnetic flux den-sity.) Computationally, we now choose the mesoscopic level such that spatial polynomial expansions of degree one suffice to approximate these continu-ous components on the boundary surface of any elementary cell. Algebraic topology learns that consistently linear edge expansions (for field quantities
whose tangential components are to be continuous) and consistently linear face expansions (for field quantities whose normal components are to be con-tinuous) on a simplicial mesh form a coherent framework in this respect.
So far, the physics of the field behavior on a mesoscopic scale dictates the computational formulation of the problem. The pertaining constitutive relations present, however, a difficulty of a different nature. Since, by assump-tion, the medium is on the mesoscopic level still heterogeneous, we cannot assign to the medium (a set of scalar or tensorial) constitutive coefficients that are to relate the magnetic flux density to the magnetic field strength. To model the magnetic constitutive behavior of the medium under consid-eration, we consistently linearly extrapolate the field quantities from their (edge or face) representations on the surface to the interior of each (simpli-cial) cell and assume that the extrapolated values allow for the introduction of ‘smoothed’ constitutive coefficients that follow from a minimization pro-cedure applied to the norm in the discrepancy in a constitutive relation that would apply to a locally homogeneous medium. In our case, we have chosen for this the L2-norm, but this choice is, of course, open to debate. Since sim-ilar arguments hold for prescribed, field independent, source magnetizations, we have followed the same procedure for this quantity.
The procedure indicated above leads to an overdetermined system of linear, algebraic equations in the field expansion coefficients. Such a system has, as we know, no solution. This, in itself, we consider as indicative that we are, again, truly modeling the physics. To substantiate this argument, let us for a moment consider the field in a homogeneous piece of matter. Here, we know that the field components are analytic functions of the spatial co-ordinates, which functions require an infinity of polynomial terms to expand them. Hence, a local representation in the form of a polynomial of degree one can, even if sufficiently accurate, only be approximate and never exact. Therefore, the local expansions used cannot (exactly) solve the problem, which property can be rephrased as that the field problem as discretized has no solution. The best we can do is to construct the ‘best’ values of the field expansion functions in accordance with the physical properties accounted for, i.e., minimizing, subject to a certain error criterion, the discrepancy between the values that are prescribed and the combinations of field expansion coef-ficients they should be equal to. Here, too, we have chosen for an L2-norm. This choice, as well, is subject to debate.
It is clear that nowhere in our analysis the differentiability of the field quantities plays a role. This sets the domain-integrated field equa-tions method apart from other methods such as the finite-difference method
(where spatial derivatives of the field vectors are somehow approximated) and the finite-element method (where a reformulation of the field problem to a weighted residual one is used, where the procedure of reformulation leans on the use of spatial derivatives of the field vectors).
Finally, it is to be remarked that if on the mesoscopic scale the medium should be (piecewise) continuous, the field quantities in any subdomain of continuity are differentiable vector function of position and the formulation in terms of field equations in differential form, with a pointwise constitutive relation, holds. Our formulation meets the criterion that under this condition the equivalence of the domain-integrated field equations formulation and the formulation in terms of field differential equations and a local constitutive re-lation is evident. Given the relevant conditions, straightforward applications of Gauss’ theorem and a continuity argument applied to the constitutive law link the two formulations.”
In this thesis we present an extension of the domain-integrated field equations method to the time-harmonic problem where both the magnetic and electric properties of the media may be strongly heterogeneous. This extended version of the domain-integrated field equations method is referred to as the Volume-Integrated Field Equations (VIFE) method. This name is chosen to distinguish the method from a second method that we discuss in this thesis, a method that uses Amp`ere’s circuital law and Faraday’s law instead of the volume-integrated field equations. The latter method is appropriately referred to as the Surface-Integrated Field Equations (SIFE) method. In the following section we describe in more detail the content and structure of the thesis.
1.3
Outline of the Thesis
We start the thesis by recapitulating the basic equations that govern clas-sical electromagnetic field theory. As the point of departure we take the time-domain Maxwell equations in differential form and the time-domain constitutive relations. From these equations the equivalent counterparts in the frequency domain are derived. Since the basic electromagnetic field equa-tions in the frequency domain are available, all other important equaequa-tions can then be derived.
The integral formulations of the basic equations with their direct phys-ical significance preceded the differential formulation which was derived from
the integral formulations assuming certain smoothness conditions. This is exactly the reason why we revert to integral formulations for the problem at hand. If the electromagnetic field is sufficiently smooth, one can establish equivalence between integral formulations and the differential formulation of the basic equations by integrating the differential equations over a surface or volume.
The equations and relations described in chapter 2 include the con-ditions imposed on the field quantities at material interfaces. According to Mur, see [38], these conditions are compatibility relations, but mostly we refer to them as interface conditions and reserve the term compatibility relations for the divergence conditions on the field quantities with respect to a volume element. To obtain a unique solution for the field problem, we must provide the outer boundary of a configuration under consideration with boundary conditions, and therefore chapter 2 is concluded by a discussion on several of these boundary conditions. If the configuration is unbounded, a situation that is frequently encountered in electromagnetics, this configuration is trun-cated using Absorbing Boundary Conditions (ABCs) or Perfectly Matched Layers (PMLs) so that computational methods can be applied. In this thesis we only indicate briefly the application of ABCs but extensively address the powerful concept of PMLs.
The first part of the remainder of the thesis, consisting of chapters 3, 4 and 5, is devoted to two-dimensional (2D) electromagnetic field problems, while the second part, chapters 6, 7 and 8, deals with very complex three-dimensional (3D) field problems. In 2D configurations, the electromagnetic field can be decoupled into a parallelly polarized field and a perpendicularly polarized field, which reduces the overall complexity considerably.
In chapter 3 we describe the geometrical discretization of the 2D com-putational domain, and we give a detailed explanation of the local and global discretization of the field quantities for both polarization cases. To satisfy the requirements imposed by the interface conditions, a combination of 2D Cartesian, 2D edge and 2D face elements is used for the discretization of the field quantities. To gain a clear understanding of the properties of the differ-ent expansions we have added to the text the construction of the expansions. In chapter 4 we describe the Volume-Integrated Field Equations (VIFE) method and the Surface-Integrated Field Equations (SIFE) method for 2D electromagnetic field problems, each proposing criteria to be satisfied by the approximated field quantities, constructed in the previous chapter. The
methods differ from each other in the way they apply Maxwell’s equations, but they both employ the constitutive relations in the same unconventional manner. At the end of the chapter we have at our disposal, for an arbi-trary elemental volume in the mesh, different sets of constraints on the local expansion coefficients for the different methods.
In chapter 5 we explain the way to obtain a global system of equa-tions for the unknown global expansion coefficients from the local sets of constraints. Subsequently, we discuss some implementation aspects and de-scribe techniques to solve the large system of equations. Direct solution methods are used when the number of unknowns is not too large. Iterative solution methods, in combination with preconditioning, are used when the direct methods are computationally too expensive. The chapter is concluded by showing how the 2D numerical methods perform on a test configuration consisting of four subdomains with different medium properties.
In chapter 6 we describe the geometrical discretization of the 3D com-putational domain and we give a detailed explanation of the local and global discretization of the field quantities. To satisfy the requirements imposed by the interface conditions, we use a combination of 3D edge and 3D face elements for the discretization of the field quantities. To gain a clear under-standing of the properties of the expansions, we have added their construction to the text.
In chapter 7 we describe the Volume-Integrated Field Equations (VIFE) method and the Surface-Integrated Field Equations (SIFE) method for 3D electromagnetic field problems, each proposing criteria to be satisfied by the approximated field quantities, constructed in the previous chapter. At the end of the chapter we have derived for an arbitrary elemental volume in the mesh different sets of constraints on the local expansion coefficients for the different methods.
In chapter 8 we follow the same reasoning as in the 2D case, to explain how to obtain a global system of equations in the unknown global expansion coefficients from the local sets of constraints. The increase in complexity of implementation is considerable when going from 2D to 3D, but the structure of the methods and hence of the Fortran program is retained and does not need further explanation. Solution methods for large systems of linear equa-tions were already described in chapter 4. We note that the direct solution method is too costly for 3D field problems and only iterative techniques are feasible. The chapter is concluded by showing how the numerical methods
perform on a number of test configurations. The four-media problem is revis-ited, Perfectly Matched Layers (PMLs) are tested, and a scattering problem is investigated, where the domain is truncated with PMLs.
In the final chapter of the thesis we evaluate the accomplishments of this thesis and offer suggestions for further research.
The Electromagnetic Field Equations
In this chapter we recapitulate the basic equations that govern classical elec-tromagnetic field theory. The summary presented is based on the theory discussed by Cheng [9], De Hoop [18], Lager [27], Mur [39] and Verweij et al. [50]. As the point of departure we take the time-domain Maxwell equations in differential form and the time-domain constitutive relations. From these equations the equivalent counterparts in the frequency domain are derived. In this thesis we will restrict ourselves to a class of problems described using the frequency-domain analysis. A theoretical discussion using a time-domain analysis is given by De Hoop et al. [21]. Since the basic electromagnetic field equations in the frequency domain are available, all other equations can then be derived. The compatibility relations follow directly from Maxwell’s equa-tions and hence a field solution must satisfy them. Integral formulaequa-tions of the basic equations are given enabling a proper treatment of discontinuities in the field quantities caused by discontinuities in the material parameters. The conditions imposed on the field quantities at the material interfaces are derived using these integral formulations and are referred to as the interface conditions. To obtain a unique solution for the electromagnetic field problem we must provide the outer boundary of the configuration under consideration with boundary conditions. If the configuration is unbounded, a situation frequently encountered in electromagnetics, this configuration is truncated using Absorbing Boundary Conditions (ABCs) or Perfectly Matched Layers (PMLs) so that computational methods can be applied.
2.1
Time-Domain Maxwell’s Equations
in Differential Form
Let D ⊂ R3 be a bounded domain with piecewise smooth boundary surface
∂D. The unit vector along the outward normal to ∂D is denoted by n. Let T ⊂ R be a time interval T = {t ∈ R | t0 < t < t1} with boundary points t0
and t1, where t1=∞ is allowed. In any (x, t) ∈ D × T ⊂ R3× R where the
electromagnetic field quantities are continuously differentiable with respect to both space and time coordinates, Maxwell’s equations
∇ × H = Jimp+ Jcon+ ∂
tD , (2.1)
∇ × E = −Kimp− Kcon− ∂
tB , (2.2)
are satisfied. In this set of first-order partial differential equations E(x, t) = electric field strength (V/m) ,
H(x, t) = magnetic field strength (A/m) , D(x, t) = electric flux density (C/m2) , B(x, t) = magnetic flux density (T) ,
Jimp(x, t) = volume density of impressed electric current (A/m2) , Jcon(x, t) = volume density of electric conduction current (A/m2) , Kimp(x, t) = volume density of impressed magnetic current (V/m2) , Kcon(x, t) = volume density of magnetic conduction current (V/m2) . Here Jimp and Kimp describe the action of sources that generate the field, their detailed physical behaviour is either irrelevant to or outside the scope of the analysis at hand. Note that we have introduced magnetic current densities to obtain symmetry and ease of modelling. Now we introduce
J (x, t) = volume density of electric current (A/m2) , K(x, t) = volume density of magnetic current (V/m2) , these quantities are defined as
J = Jimp+ Jcon+ ∂tD , (2.3)
Included in the volume density of electric current are not only the volume density of impressed electric current and the volume density of electric con-duction current, also included is the so-called volume density of electric dis-placement current. Similarly, included in the volume density of magnetic current is the so-called volume density of magnetic displacement current. Using these new defined quantities, we rewrite Maxwell’s equations (2.1) and (2.2) to get a very succinct and convenient form
∇ × H = J , (2.5)
∇ × E = −K . (2.6)
2.2
Time-Domain Constitutive Relations
We assume the media present in the configurations are linear, time-invariant and locally reacting in their electromagnetic behaviour. Then the constitutive relations are Jcon(x, t) = ∞ t=0 ¯ σe(x, t)· E(x, t − t) dt, (2.7) Kcon(x, t) = ∞ t=0 ¯ σm(x, t)· H(x, t − t) dt, (2.8) D(x, t) = ∞ t=0 ¯ ε(x, t)· E(x, t − t) dt, (2.9) B(x, t) = ∞ t=0 ¯ µ(x, t)· H(x, t − t) dt, (2.10) where ¯
σe(x, t) = electric conduction relaxation tensor (S/ms) , ¯
σm(x, t) = magnetic conduction relaxation tensor (Ω/ms) , ¯
ε(x, t) = permittivity relaxation tensor (F/ms) , ¯
µ(x, t) = permeability relaxation tensor (H/ms) .
Note that the causality of the medium’s response has been enforced by restricting the interval of integration of the elapse time t to the interval {t ∈ R | t > 0}, which implies that only the values of {E, H} prior to t
contribute to the values of {Jcon, Kcon, D, B} at the instant t. For the lim-iting case of an instantaneously reacting material, the relaxation functions
approach an impulse (delta-distribution) time behaviour. That is ¯ σe(x, t) = ¯σe(x)δ(t) , (2.11) ¯ σm(x, t) = ¯σm(x)δ(t) , (2.12) ¯ ε(x, t) = ¯ε(x)δ(t) , (2.13) ¯ µ(x, t) = ¯µ(x)δ(t) , (2.14)
where δ(t) is the Dirac impulse distribution operative at t = 0 and ¯σe(x),
¯
σm(x), ¯ε(x) and ¯µ(x) are the electric conductivity, magnetic conductivity, permittivity and permeability tensor, respectively. For instantaneously re-acting media the constitutive relations (2.7), (2.8), (2.9) and (2.10) become
Jcon(x, t) = ¯σe(x)· E(x, t) , (2.15)
Kcon(x, t) = ¯σm(x)· H(x, t) , (2.16)
D(x, t) = ¯ε(x)· E(x, t) , (2.17)
B(x, t) = ¯µ(x)· H(x, t) . (2.18)
Although not necessary, in this thesis we assume that only the diagonal ele-ments of the material tensors are nonzero. If the three diagonal eleele-ments are different, the material is biaxially anisotropic. If two of the three elements have the same value, the material is uniaxially anisotropic. In case all diago-nal elements are the same, the material is isotropic. For an isotropic medium with relaxation we have
¯ σe(x, t) = σe(x, t) ¯I , (2.19) ¯ σm(x, t) = σm(x, t) ¯I , (2.20) ¯ ε(x, t) = ε(x, t) ¯I , (2.21) ¯ µ(x, t) = µ(x, t) ¯I , (2.22)
and for an instantaneously reacting isotropic medium ¯ σe(x) = σe(x) ¯I , (2.23) ¯ σm(x) = σm(x) ¯I , (2.24) ¯ ε(x) = ε(x) ¯I , (2.25) ¯ µ(x) = µ(x) ¯I , (2.26)
2.3
Frequency-Domain Representations
An electromagnetic phenomenon is transient in nature. Essential feature is the physical condition of causality that requires that an effect does not occur before its cause, this implies that a field that is causally related to the action of a source that is switched on at the instant t = t0 necessarily vanishes prior
to this instant, that is for t < t0. In linear time-invariant configurations the
behaviour of such transient fields can advantageously be described using a frequency-domain representation. The dependence on t is hereby replaced by a dependence on a transform parameter.
Laplace Transformation
Let f = f (x, t) denote an electromagnetic field quantity that is causally related to the action of a source that has been switched on at the instant t = t0, then its Laplace transform (with respect to time) is defined as
ˆ
f (x, s) = ∞
t=t0
exp(−st)f(x, t) dt , for Re(s) > s0, (2.27)
where s is the complex Laplace transform parameter and ˆf is analytic in the right half of the complex s-plane Re(s) > s0. An important property of the
Laplace transformation is that the Laplace transform of ∂tf (x, t) is s ˆf (x, s)
if f (x, t) is causal and f (x, t0) = 0.
Once ˆf = ˆf (x, s) has been evaluated, f = f (x, t) can be recovered by the application of the Bromwich inversion integral
f (x, t) = 1 2πj
s0+j∞
s=s0−j∞
exp(st) ˆf (x, s) ds , for all t, (2.28) that automatically yields the value zero for t < t0 and for which numerical
routines are available. In equation (2.28), j is the imaginary unit (j2 =−1). Fourier Transformation
In case f = f (x, t) decays sufficiently fast in magnitude as t → ∞, we are allowed to consider ˆf = ˆf (x, s) for imaginary values s = jω, where ω is the (real) angular frequency. Then equation (2.27) can be written as
ˆ
f (x, jω) = ∞
t=t0
exp(−jωt)f(x, t) dt , for all real ω. (2.29) Important property of the Fourier transformation is that the Fourier trans-form of ∂tf (x, t) is jω ˆf (x, jω) if f (x, t0) = 0.
Once ˆf = ˆf (x, jω) has been evaluated, f = f (x, t) can be recovered by the application of the inverse Fourier transform
f (x, t) = 1 2π
∞
ω=−∞
exp(jωt) ˆf (x, jω) dω , for all t. (2.30) Steady State Analysis
In a steady state analysis all electromagnetic field quantities are taken to depend sinusoidally on time with a common angular frequency ω. In such an analysis a real quantity f (x, t) can be represented by a phasor ˆf (x, jω) and time factor exp(jωt). The original quantity is obtained from its complex counterpart by f (x, t) = Re ˆ f (x, jω) exp(jωt) . (2.31)
It is easily seen that ∂tf (x, t) is represented by the phasor jω ˆf (x, jω).
2.4
Frequency-Domain Maxwell’s Equations
We leave the spatial domain unchanged, that isD ⊂ R3 is the same bounded domain with piecewise smooth boundary surface ∂D. The unit vector along the outward normal is again denoted by n. Upon applying the Laplace transformation (2.27) to equations (2.1) and (2.2), and the property ˆ∂t= s,
we obtain for all x ∈ D where the electromagnetic field quantities are con-tinuously differentiable with respect to the space coordinates, the following relations
∇ × ˆH = Jˆimp+ ˆJcon+ s ˆD , (2.32) ∇ × ˆE =− ˆKimp− ˆKcon− s ˆB . (2.33) The frequency-domain counterparts of J and K are obtained by transforming equations (2.3) and (2.4) giving
ˆ
J = ˆJimp+ ˆJcon+ s ˆD , (2.34) ˆ
K = ˆKimp+ ˆKcon+ s ˆB . (2.35) Upon using these expressions, Maxwell’s equations (2.32) and (2.33) can be given a concise and orderly form, namely
∇ × ˆH = J ,ˆ (2.36)
Compatibility Relations
Let us assume the field quantities are sufficiently smooth. Upon applying the divergence operator ∇· to Maxwell’s equations (2.36) and (2.37) it follows, using ∇ · (∇ × . . .) = 0, that
∇ · ˆJ = 0 , (2.38)
∇ · ˆK = 0 . (2.39)
Equations (2.38) and (2.39) are known as the electromagnetic compatibility relations. We conclude that a field solution that satisfies Maxwell’s equa-tions (2.36) and (2.37) must satisfy the compatibility relaequa-tions. Logically equivalent with the preceding conclusion is the observation that if a possible field solution does not satisfy the compatibility relations, it cannot satisfy Maxwell’s equations. In the remainder of this chapter we choose to apply the Laplace transformation. Note that applying the Fourier transformation or a steady state analysis provides us with similar equations, we only need to replace s by jω.
2.5
Frequency-Domain Constitutive Relations
The frequency-domain constitutive relations for linear, time-invariant and locally reacting media are obtained by transforming equations (2.7) through (2.10). Application of the Laplace transform to a convolution of two functions results in a product of the transforms of these functions, hence
ˆ Jcon(x, s) = ˆσ¯e(x, s)· ˆE(x, s) , (2.40) ˆ Kcon(x, s) = ˆσ¯m(x, s)· ˆH(x, s) , (2.41) ˆ D(x, s) = ˆε(x, s)¯ · ˆE(x, s) , (2.42) ˆ B(x, s) = ˆµ(x, s)¯ · ˆH(x, s) . (2.43) For instantaneously reacting media we transform constitutive relations (2.15) through (2.18) to arrive at ˆ Jcon(x, s) = ¯σe(x)· ˆE(x, s) , (2.44) ˆ Kcon(x, s) = ¯σm(x)· ˆH(x, s) , (2.45) ˆ D(x, s) = ¯ε(x)· ˆE(x, s) , (2.46) ˆ B(x, s) = ¯µ(x)· ˆH(x, s) . (2.47)
For an isotropic medium with relaxation we have ˆ ¯ σe(x, s) = ˆσe(x, s) ¯I , (2.48) ˆ ¯ σm(x, s) = ˆσm(x, s) ¯I , (2.49) ˆ ¯ ε(x, s) = ˆε(x, s) ¯I , (2.50) ˆ ¯ µ(x, s) = ˆµ(x, s) ¯I , (2.51)
and for an instantaneously reacting isotropic medium relations (2.23) through (2.26) are still valid. By substituting the constitutive relations (2.40) through (2.43) in the expressions (2.34) and (2.35) for ˆJ and ˆK we obtain
ˆ
J = ˆJimp+ (ˆσ¯e+ sˆε)¯ · ˆE , (2.52) ˆ
K = ˆKimp+ (ˆσ¯m+ s ˆµ)¯ · ˆH . (2.53) Now the electric conduction current, magnetic conduction current, electric flux density and magnetic flux density are implicitly present in our modelling. Let us introduce ˆ ¯ η = ˆσ¯e+ sˆε ,¯ (2.54) ˆ ¯ ζ = ˆσ¯m+ s ˆµ ,¯ (2.55)
where ˆη and ˆ¯ ζ are the (tensorial) transverse admittance per unit length of¯ the medium (S/m) and the longitudinal impedance per unit length of the medium (Ω/m), respectively. This way, equations (2.52) and (2.53) become
ˆ
J = ˆJimp+ ˆη¯· ˆE , (2.56)
ˆ
K = ˆKimp+ ˆζ¯· ˆH . (2.57)
Usually (2.56) and (2.57) are substituted in (2.36) and (2.37) reducing the number of variables in the model. Often, a further reduction of the number of variables in the model is obtained by turning the set of first-order Partial Differential Equations (PDEs) into a second-order PDE in either ˆE or ˆH. Essential in our modelling will be the proper treatment of the interface con-ditions described in section 2.8, therefore we will use the more elaborate set of equations (2.36), (2.37), (2.56) and (2.57).
2.6
Surface-Integrated Field Equations
In explaining electromagnetic phenomena in a physical environment we have to deal with finite objects of specified shapes and boundaries. It is convenient
to convert the differential forms into their integral equivalents. Let us con-sider an arbitrary open surfaceA with surface-bounding contour ∂A. Taking the surface integral (of the component normal to the surface) of both sides of the curl equations (2.36) and (2.37) overA and applying Stokes’s theorem gives ∂A ˆ H· dl = A ˆ J· ν dA , (2.58) ∂A ˆ E· dl = − A ˆ K· ν dA , (2.59)
where the unit normal ν and the direction of circulation are related by the right-hand rule. Equation (2.58) is also known as Amp`ere’s circuital law whereas (2.59) is well-known as Faraday’s law.
Compatibility Relations
It is easily derived using equations (2.58) and (2.59) that for an arbitrary closed surface A with outwardly oriented unit normal ν
A ˆ J · ν dA = 0 , (2.60) A ˆ K· ν dA = 0 . (2.61)
These equations are in fact the compatibility relations in integral form.
2.7
Volume-Integrated Field Equations
Integrating Maxwell’s equations in differential form (2.36) and (2.37) over an arbitrary volume V with piecewise smooth boundary surface ∂V, and applying the Curl theorem gives
∂V ν× ˆH dA = V ˆ J dV , (2.62) ∂V ν× ˆE dA =− V ˆ K dV , (2.63)
where the unit vector along the outward normal to ∂V is denoted by ν. Compatibility Relations
over an arbitrary volume V with piecewise smooth boundary surface ∂V and applying Gauss’s theorem gives
∂V ˆ J · ν dA = 0 , (2.64) ∂V ˆ K· ν dA = 0 , (2.65)
where the unit vector along the outward normal to ∂V is denoted by ν. Note that in section 2.6 these relations were obtained by appropriately applying Amp`ere’s circuital law and Faraday’s law.
2.8
Interface Conditions
In almost all relevant configurations the media are inhomogeneous. The pres-ence of materials with different electromagnetic properties close to each other is often best described by introducing an abrupt finite jump in the material parameters. This is modelled by allowing for discontinuities in the constitu-tive parameters. We will assume the discontinuity interfaces are piecewise smooth. At a point x on such an interface, Maxwell’s equations in differen-tial form (when classically interpreted) cease to hold. Maxwell’s equations in integral form and the corresponding compatibility relations are still well defined. Due to the discontinuity of the material properties some components of the electromagnetic field quantities may show finite jumps at the interface, but appropriately applying Maxwell’s equations in integral form will link the fields on both sides.
At an interface between media 1 and 2, we denote by ν the unit vector normal to the interface that points from medium 1 to 2 . Applying Amp`ere’s circuital law (2.58) at the interface gives
ν× ˆH2
1 = ˆJS, (2.66)
where the left-hand side represents the jump in the tangential component of ˆH across the interface between media 1 and 2. Here ˆJS is the surface
density of electric current (A/m) at the interface; ˆJS satisfies the condition
ˆ
JS· ν = 0. In figure 2.1(a) we have depicted the situation described above.
The plane S1 is located in medium 1 and is parallel to the material interface
at a distance ∆h/2 from the material interface. The plane S2 is located in
∆h 2 ∆h2 ˆ Ht1 ˆ Ht2 ˆ JS S1 S S2 (a) ∆h 2 ∆h2 ˆ Ht1 Hˆ t 2 S1 S S2 (b)
Figure 2.1. The tangential magnetic field at a material interface with a surface electric current (a) and without a surface electric current (b).
∆h 2 ∆h2 ˆ Jn1 Jˆ n 2 S1 S S2 (a) ∆h 2 ∆h2 ˆ Kn1 Kˆ n 2 S1 S S2 (b)
Figure 2.2. The normal component of the volume density of electric current (a) and the normal component of the volume density of magnetic current (b) at a material interface. ∆h 2 ∆h2 ˆ Et1 ˆ Et2 ˆ KS S1 S S2 (a) ∆h 2 ∆h2 ˆ Et1 Eˆt2 S1 S S2 (b)
Figure 2.3. The tangential electric field at a material interface with a surface mag-netic current (a) and without a surface magmag-netic current (b).
the material interface and hence at a distance ∆h from S1. In figure 2.1(a)
we also show the jump in the tangential magnetic field strength going from plane S1 to S2 due to the surface density of electric current present at the
material interface. In the limiting case, where ∆h→ 0, equation (2.66) holds. Assuming the interface is free of surface electric currents, equation (2.66) is rewritten as
ν× ˆH2
1 = 0 ⇐⇒ ν × ˆH = continuous . (2.67)
In figure 2.1(b) we have depicted the situation without surface electric cur-rents. Shown is the continuity of the tangential magnetic field strength going from plane S1 to S2. In the limiting case where, ∆h → 0, equation (2.67)
holds. Applying compatibility relation (2.60) at an interface between media 1 and 2 gives
ˆ J· ν2
1 = 0 ⇐⇒ ˆJ · ν = continuous . (2.68)
In figure 2.2(a) we have shown the continuity of the normal component of the volume density of electric current going from plane S1 toS2. In the limiting
case, where ∆h→ 0, equation (2.68) holds.
Analogously applying Faraday’s law (2.59) at an interface between me-dia 1 and 2 gives
ν× ˆE2
1=− ˆKS, (2.69)
where the left-hand side represents the jump in the tangential component of ˆE across the interface between media 1 and 2. Here ˆKS is the surface
density of magnetic current (V/m) at the interface; ˆKSsatisfies the condition
ˆ
KS· ν = 0. In figure 2.3(a) we have depicted the behaviour of the tangential
electric field strength at a material interface with a surface magnetic current. Assuming the interface is free of surface magnetic currents, equation (2.69) is rewritten as
ν× ˆE2
1= 0 ⇐⇒ ν × ˆE = continuous . (2.70)
The latter case is depicted in figure 2.3(b). Applying compatibility relation (2.61) at an interface between media 1 and 2 gives
ˆ K· ν2
In figure 2.2(b) we have depicted the behaviour of the normal component of the volume density of magnetic current at a material interface.
2.9
Boundary Conditions
Since numerical computations can only be carried out in some bounded sub-domain of space, the question arises as to what boundary conditions on the boundary ∂D of the domain D of computation should supplement the elec-tromagnetic field equations that are solved inD such that the physics of the radiation of the field into the embedding D of D is adequately modelled. Mathematically such a modelling must in any case lead to the existence of a unique solution, moreover it is desired that the solution in the bounded do-main D accurately approximates the solution of the original unbounded field problem. Note that uniqueness of the electromagnetic wavefield in the orig-inal unbounded field problem is guaranteed, even for dispersive anisotropic media, via the time-domain uniqueness theorem by De Hoop [23].
In this thesis prescribed boundary conditions are considered, while the topic of Absorbing Boundary Conditions (ABCs) is shortly addressed. A uniqueness proof in the case of prescribed boundary conditions requires that the Laplace transforms of the electric and magnetic relaxation tensors of rank two are positive definite for all real, positive values of s. ABCs are succinctly described since they are extensively used in electromagnetics even though proof of the uniqueness of the solution to the corresponding electromagnetic wavefield problem is often not rigorously established. We conclude the section by introducing the so-called Perfectly Matched Layers (PMLs) that offer an alternative means of truncating a computational domain, and that satisfy the uniqueness properties as is demonstrated by De Hoop et al. [22].
Prescribed Tangential Electric Field
The tangential component of the electric field strength ˆE may be given a prescribed value ˆEextt at the outer boundary ∂D of the configuration D, where
n· ˆEextt = 0 , (2.72)
and n the unit vector along the normal to the outer boundary. Thus we have for the electric field strength ˆE at the outer boundary
In order to ensure compatibility with the prescribed tangential electric field strength ˆEextt at the outer boundary ∂D, the normal component of ˆK has to satisfy
ˆ
K· n = −(∇S× ˆEextt )· n , (2.74)
where∇S is the nabla operator restricted to the plane tangential to the outer
boundary in the point under consideration.
In an electrically impenetrable object the electric constitutive para-meters have infinitely large values so that the electric field strength has a negligibly small value. Consequently, if D is bounded by such an object, the tangential electric field strength ˆEextt at the outer boundary ∂D is prescribed to be zero, thus we have for the electric field strength ˆE at the boundary
n× ˆE = 0 . (2.75)
For this special case equation (2.74) becomes ˆ
K· n = 0 . (2.76)
The latter boundary conditions are referred to as Perfect Electric Conducter (PEC) boundary conditions.
Prescribed Tangential Magnetic Field
The tangential component of the magnetic field strength ˆH may be given a prescribed value ˆHextt at the outer boundary ∂D of the configuration D, where
n· ˆHextt = 0 , (2.77)
and n the unit vector along the normal to the outer boundary. Thus we have for the magnetic field strength ˆH at the outer boundary
n× ˆH = n× ˆHextt . (2.78)
In order to ensure compatibility with the prescribed tangential magnetic field strength ˆHextt at the outer boundary ∂D, the normal component of ˆJ has to satisfy
ˆ
J· n = (∇S× ˆHextt )· n , (2.79)
where∇S is the nabla operator restricted to the plane tangential to the outer
In a magnetically impenetrable object the magnetic constitutive parameters have infinitely large values so that the magnetic field strength has a negligibly small value. Consequently, ifD is bounded by such an object, the tangential magnetic field strength ˆHextt at the outer boundary ∂D is prescribed to be zero, thus we have for the magnetic field strength ˆH at the boundary
n× ˆH = 0 . (2.80)
For this special case equation (2.79) becomes ˆ
J· n = 0 . (2.81)
The latter boundary conditions are referred to as Perfect Magnetic Conducter (PMC) boundary conditions.
Prescribed Tangential Electric and Magnetic Field
A combination of previously described boundary conditions is allowed, that is on some part ∂De of ∂D the tangential electric field has a prescribed value and on some other part ∂Dm of ∂D the tangential magnetic field has a prescribed value, provided that ∂De∪ ∂Dm= ∂D and ∂De∩ ∂Dm=∅. It is obvious that equations (2.73) through (2.76) apply at ∂De and equations (2.78) through (2.81) apply at ∂Dm.
Absorbing Boundary Conditions
In electromagnetic field problems many geometries of interest are defined in open regions where the spatial domain of the field is unbounded in one or more coordinate directions. Application of numerical methods on a bounded computational domainD is only possible if we truncate the unbounded con-figuration with an artificial boundary. The computational domainD must be large enough to enclose the structure of interest and the boundary condition on the outer perimeter of the domain must adequately simulate the extension to infinity. Such a condition should make the boundary appear as transpar-ent as possible, that is, it should minimize the nonphysical reflections from that boundary. A class of boundary conditions designed for this purpose is called Absorbing Boundary Conditions (ABCs).
ABCs provide us with a relation between the field quantities at the boundary ∂D such that radiation of the electromagnetic field in D into the surrounding D of D is simulated. Unlike boundary integral equations or eigenfunction expansions, ABCs lead to localized relations between the boundary fields, and the highly sparse and banded nature of the system matrices of local methods like the Finite-Element (FE) method and the
region of interest background medium χx= χy= 1 χx= 1 χy= n−y − jm−y χx= 1 χy= n+ y − jm+y χx= n−x− jm−x χy= 1 χx= n+ x − jm+x χy= 1 χx= n−x − jm−x χy= n−y − jm−y χx= n+ x − jm+x χy= n−y− jm−y χx= n−x − jm−x χy= n+ y− jm+y χx= n+ x − jm+x χy= n+ y− jm+y x−PEC x−PML x+PML x+PEC y−PEC yPML− y+PEC y+ PML
Figure 2.4. Use of PEC-backed PMLs for truncating the computational domain.
Finite-Difference (FD) method is retained. The performance of these ar-tificial boundaries depends, among other things, on the order of the ABC and the distance of the boundary from the region of interest.
Perfectly Matched Layers
Nowadays so-called Perfectly Matched Layers (PMLs), introduced by B´erenger [5], offer an alternative way of truncating a computational domain. In the literature we find different interpretations of the PML. In particular Chew and Weedon [10] have interpreted the PML as a coordinate stretching in the frequency domain. Here we adopt the interpretation of the PML as an anisotropic absorber. The anisotropic absorber model of the PML was first derived by Sacks et al. [46] and is especially suited for Finite-Element (FE) or Finite-Difference (FD) implementation. We use the formulation derived by Gedney [14]. A short introduction on PMLs, the coordinate stretching interpretation, the anisotropic absorber formulation derived by Gedney and references to relevant literature on the topic are given in the book by Jin [24].
A perfectly matched interface is an interface between two half-spaces, that does not reflect an incident field for all frequencies and all angles of incidence and polarizations. Each homogeneous half-space with arbitrary medium parameters allows for the construction of a neighbouring lossy half-space, consisting of a so-called PML-medium, such that the interface is perfectly matched. The field in the PML-medium is attenuated in the direction nor-mal to the interface. As such, the PML-medium can be used to truncate the computational domain. An artificial boundary is placed in the half-space with the PML-medium at a distance from the perfectly matched interface where the field is expected to have been sufficiently attenuated. At this outer boundary often simple PEC/PMC boundary conditions will suffice but additional absorption is acquired with an ABC boundary condition.
In 3D field problems where the spatial domain is unbounded in all coordinate directions, we need six perfectly matched layers to truncate the computational domain. The basic scheme, showing a 2D intersection of a 3D domain, is sketched in figure 2.4. The region of interest, embedded in a homogeneous background medium, is surrounded by PMLs which are termi-nated by a metal (PEC). In the article by Gedney [14] and the book by Jin [24] the material parameters of the anisotropic PML-medium are chosen as follows ˆ ¯ ηPML = χyχz χx ηˆx,x 0 0 0 χxχz χy ηˆy,y 0 0 0 χxχy χz ηˆz,z , (2.82) ˆ ¯ ζPML = χyχz χx ζˆx,x 0 0 0 χxχz χy ˆ ζy,y 0 0 0 χxχy χz ζˆz,z . (2.83)
Here ˆηx,x, ˆηy,y, ˆηz,z and ˆζx,x, ˆζy,y, ˆζz,z are the material parameters of the
background medium. Note that the PML-medium is anisotropic even if the background medium is isotropic. The choice of the PML-parameters depends
on the so-called stretch factors χx(x), χy(y), χz(z) that depend on position χx(x) = n−x(x)− jm−x(x) x−PEC< x < x−PML, 1 x−PML< x < x+PML, n+x(x)− jm+x(x) x+PML< x < x+PEC, (2.84) χy(y) =
n−y(y)− jm−y(y) y−PEC< y < y−PML, 1 y−PML< y < y+PML, n+y(y)− jm+y(y) y+PML< y < y+PEC,
(2.85) χz(z) = n−z(z)− jm−z(z) z−PEC < z < z−PML, 1 z−PML < z < z+PML, n+z(z)− jm+z(z) z+PML < z < z+PEC, (2.86)
this is also illustrated in figure 2.4. Here n−x, n+x, n−y, n+y, n−z, n+z ≥ 1 are real numbers that cause faster decay of evanescent waves whereas m−x, m+x, m−y, m+y, m−z, m+z ≥ 0 are real numbers that determine the atten-uation of the propagating waves.
Discretization of the 2D Domain of
Computation
In this chapter we describe the discretization of the computational domain and of the field quantities for 2D electromagnetic field problems. First we dis-cuss the decoupling of the 2D electromagnetic field in the parallelly polarized field and the perpendicularly polarized field. Next we describe the geometri-cal discretization of the domain. The remainder of the chapter is devoted to the discretization of the field quantities where a strong emphasis is laid on the requirements imposed by the interface conditions. A detailed explanation of the local and global discretization is provided for both polarization cases.
3.1
The 2D Electromagnetic Field
In chapters 3 through 5 we describe electromagnetic field problems in con-figurations that have no variation in the z-direction. In these concon-figurations the permittivity, the permeability, the electric conductivity, and the magnetic conductivity are independent of the z-coordinate. Hence
ˆ ¯ ε = ˆε(x, y, s) ,¯ σˆ¯e = ˆσ¯e(x, y, s) , (3.1) ˆ ¯ µ = ˆµ(x, y, s) ,¯ σˆ¯m= ˆσ¯m(x, y, s) . (3.2)
Then ˆη and ˆ¯ ζ defined in terms of these material tensors in (2.54) and (2.55)¯ are independent of the z-coordinate
ˆ ¯
η = ˆη(x, y, s) ,¯ ζ = ˆˆ¯ ζ(x, y, s) .¯ (3.3) We further assume that the sources are independent of z, that is
ˆ
Jimp= ˆJimp(x, y, s) , Kˆimp = ˆKimp(x, y, s) . (3.4) If the boundary conditions at the outer boundary are also independent of z, then the generated electromagnetic field will be independent of z, thus
ˆ
E = ˆE(x, y, s) , H = ˆˆ H(x, y, s) , (3.5) ˆ
J = ˆJ (x, y, s) , K = ˆˆ K(x, y, s) . (3.6) Therefore ∂z = 0 in (2.36), (2.37). With this information and the assumption
that the off-diagonal elements of the material tensors in (2.56) and (2.57) are zero, we separate the electromagnetic field problem in two independent problems reducing the overall complexity considerably.
3.1.1 Parallel Polarization
In the first set of equations the components ˆEx, ˆEy, ˆJx, ˆJy, ˆHz, ˆKzare coupled.
These equations are
∂yHˆz= ˆJx, (3.7) −∂xHˆz= ˆJy, (3.8) ∂xEˆy− ∂yEˆx=− ˆKz, (3.9) ˆ Jx= ˆJximp+ ˆηx,xEˆx, (3.10) ˆ Jy = ˆJyimp+ ˆηy,yEˆy, (3.11) ˆ Kz= ˆKzimp+ ˆζz,zHˆz. (3.12)
Since the electric field vector is parallel to the (x, y)-plane, the plane of observation, this electromagnetic field is called to be parallelly polarized. Since the magnetic field vector is pointing in the z-direction, perpendicular to the (x, y)-plane, this electromagnetic field is also referred to as the H-polarized field.
3.1.2 Perpendicular Polarization
In the second set of equations the components ˆHx, ˆHy, ˆKx, ˆKy, ˆEz, ˆJz are
coupled. These equations are
−∂yEˆz= ˆKx, (3.13) ∂xEˆz= ˆKy, (3.14) ∂xHˆy− ∂yHˆx= ˆJz, (3.15) ˆ Kx= ˆKximp+ ˆζx,xHˆx, (3.16) ˆ Ky = ˆKyimp+ ˆζy,yHˆy, (3.17) ˆ Jz= ˆJzimp+ ˆηz,zEˆz. (3.18)
Since the electric field vector is pointing in the z-direction, perpendicular to the (x, y)-plane, the plane of observation, this electromagnetic field is called to be perpendicularly polarized or is referred to as the E-polarized field.
3.2
Geometrical Discretization
In the 2D electromagnetic field problem the configuration is infinitely ex-tended in both the positive and negative z-direction and is invariant with respect to the z-direction. Using the invariance of the field in the z-direction, the domain can be reduced to the 2D-domain D obtained by projecting the original domain on the horizontal (x, y)-plane. The geometrical discretization of the domain of computationD is achieved by subdividing it into small sub-domains called finite elements. The union of all finite elements represents the mesh. There are many possible choices for the shape of the finite elements. For topological reasons, see Naber [41], simplexes are chosen and we speak of a simplicial decomposition of the domain of computation. The 2D simplex is the triangle and correspondingly the mesh is referred to as triangulation. In figure 3.1 we give an example showing a finite-element discretization of the domain D. Note that a triangle in the mesh can be interpreted as the projection on the (x, y)-plane of a prism with arbitrary length lz in the
in-variant direction. The notation S is used to refer to an unspecified triangle. Just as in Naber [41] we takeS to be open in D. The closure of S is denoted by ¯S. The boundary ∂S = ¯S\S of the triangle S consists of the edges and vertices that delimit the triangle. The surface area within the edges, that is the area ofS, is denoted by A. The total number of triangles in the mesh is
!["7599 dni Drugiej Rzeczypospolitej : antologia reportażu międzywojennego", wybór i oprac. E. Sabelanka, K. Koźniewska, Warszawa 1983 : [recenzja]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==)