Direct and Inverse Methods for Waveguides and Scattering Problems in the Time Domain

Full text

(1)

(2)

(3)

(4)

(5) .

(6)

(7)

(8)

(9)

(10) !" #$ %&.

(11) . % '()'* $% +)((*)','*)+ -- -- )',.

(12)

(13)

(14)

(15)

(16)

(17) !!" #!$#" % & ' % % (& &) *&

(18)

(19) +

(20)

(21) ,+&) -

(22) .) !!")

(23)

(24) & % / '

(25) ,

(26) ' (

(27) & *

(28) ) -

(29)

(30) )

(31)

(32)

(33)

(34) #!) 0 ) ) ,1 2#3""3432)

(35)

(36)

(37)

(38)

(39)

(40)

(41) ' & '

(42) %

(43) & +

(44) + & &

(45)

(46)

(47) ) %%

(48) % 3

(49) %

(50) 3%%

(51) 5*6.

(52) %

(53) 3

(54) 5.*6 & % +7 8

(55)

(56) & &) , %

(57) % & %%

(58)

(59) ' % 9 &

(60) & ') -

(61) % &

(62) &

(63) ' &

(64) ) *&. & & &

(65) ' % :

(66) '

(67)

(68)

(69)

(70)

(71) &

(72)

(73)

(74) % & %

(75) :

(76)

(77)

(78)

(79)

(80) 516) ,

(81) & . % & &

(82) &) / ' % %

(83)

(84)

(85)

(86)

(87) +

(88) %

(89)

(90)

(91)

(92) %) *& : %

(93)

(94)

(95)

(96) % + ' ) %8

(97) &% % % &

(98) % & 5(;6

(99) '

(100)

(101)

(102) 5-16

(103) .* ) (

(104) . ' '

(105)

(106)

(107)

(108) &

(109) & %

(110) ) <

(111) +& %

(112) + 30!1

(113) % %%

(114) % + '

(115)

(116) '

(117) & '

(118) ) .

(119)

(120)

(121) &

(122) ' &

(123) ' + ' '

(124) %8

(125) 3

(126) % & =& > 8

(127) ) ='

(128) %

(129) & *

(130) .* %

(131)

(132) + ' )

(133) &

(134)

(135)

(136) & % +% ) -

(137)

(138)

(139)

(140)

(141) 3

(142) % ) 3

(143)

(144) >

(145) &

(146) +& '

(147)

(148) ' & 9

(149) ) *& & '

(150) % % +7 8

(151)

(152)

(153) ' .

(154) (;) ,%

(155)

(156)

(157) &

(158) '

(159) &

(160)

(161)

(162) . ) (

(163) % '

(164) ' %

(165)

(166)

(167) ' &

(168)

(169) ) ?

(170) %

(171) %

(172)

(173)

(174) %

(175) &

(176) ' ) +7 8

(177) 3

(178) * %

(179) 3

(180) & &

(181) & + ' (;

(182)

(183) ' 9

(184) '

(185) 3

(186) >

(187) ! "#$

(188) %

(189) $ & ' (()$ $ *)+,-+ $ @ .: -

(190) !!" ,, #"#3# ,1 2#3""3432

(191) $

(192)

(193) $$$ 3!#4 5& $AA

(194) ):)A B

(195) C

(196) $

(197)

(198) $$$ 3!#46.

(199) To Viktor and Vera.

(200)

(201) List of Papers. This thesis is based on the following papers, which are referred to in the text by their Roman numerals. I. II. III. IV. V. ∗ Reprints. Abenius, E., Andersson, U., Edelvik, F., Eriksson, L., Ledfelt, G. (2002) Hybrid Time Domain Solvers for the Maxwell Equations in 2D. International Journal for Numerical Methods in Engineering, 53:2185-2199.∗ Abenius, E., Edelvik, F., Johansson, C. (2005) Waveguide Truncation Using UPML in the Finite-Element Time-Domain Method. Technical Report, 2005-026, Department of Information Technology, Uppsala University. Abenius, E., Strand, B. (2005) Solving Inverse Electromagnetic Problems Using FDTD and Gradient-based Minimization. Submitted to International Journal for Numerical Methods in Engineering. Abenius, E., Edelvik, F. (2005) Thin Sheet Modeling Using Shell Elements in the Finite-Element Time-Domain Method. Accepted for publication in IEEE Transactions on Antennas and Propagation. Abenius, E., Johansson, C. (2005) Modeling of Inhomogeneous Waveguides using Hybrid Methods. Submitted to International Journal of Numerical Modelling: Electronic Networks, Devices and Fields. were made with permission from the publishers.. v.

(202)

(203) Contents. 1 2. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Time-domain hybrid solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 The finite-difference time-domain method . . . . . . . . . . . . . . . . 2.2.1 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.3 Plane wave excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.4 Absorbing boundary condition (ABC) . . . . . . . . . . . . . . . 2.2.5 Dispersive media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.6 Subcell models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.7 Parallelization by domain decomposition . . . . . . . . . . . . . 2.3 The finite-element time-domain method . . . . . . . . . . . . . . . . . . 2.3.1 Absorbing boundary conditions . . . . . . . . . . . . . . . . . . . . 2.3.2 Subcell models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Hybridization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Waveguides . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Inverse scattering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Difficulties in the inverse problem . . . . . . . . . . . . . . . . . . . . . . 4.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Minimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Computing gradients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Summary in Swedish: Direkta och inversa metoder för vågledare och spridningsproblem i tidsdomänen . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 1 5 6 8 8 10 10 11 11 12 12 14 17 18 19 21 25 26 27 28 29 31 35 37 41. vii.

(204)

(205) 1. Introduction. In the modern society electromagnetic fields are used in many different devices and for many different purposes. Common examples include wireless communication in antenna systems for mobile phones and satellites, heating in microwave ovens and medical X-ray imaging equipment and tomography. As a result it has become increasingly important to have a good understanding of the electromagnetic field and how it interacts with the environment. The propagation of the field is described by a system of partial differential equations (PDEs) called Maxwell’s equations, [48]. Even though the equations are linear for most applications, the added complexity of boundary conditions and materials precludes the exact solution of most realistic problem. Since experiments are usually very costly and time-consuming, numerical simulation often represents the most efficient tool for analysis and design. Computational electromagnetics (CEM) has experienced a tremendous growth in the last decades due to the availability of powerful computers. Some of the applications that receive most attention are • • • • •. prediction of radar cross-section (RCS), antenna coupling and analysis, microwave devices and guided structures, electromagnetic compatibility (EMC), bioelectromagnetics, e.g. tomography or determination of the specific absorption rate (SAR) for radiation from mobile phones.. Due to the computational power of modern computers the area of design optimization and inverse problems has received an increased interest. In these problems the physical properties are not fully characterized but contain unknown parameters. Examples include antenna design and optimization, and the design of radar absorbing materials to minimize the RCS. Also, reflection measurements are used to determine the characteristics of unknown media. Similar problems arise in medical applications where X-rays are used to examine the human body, for instance to locate tumors. The work presented in this thesis has been carried out within the Parallel and Scientific Computing Institute (PSCI), mainly in the projects GEMS (General Electromagnetic Solver) and IMPACT (Inverse Methods for Wave Propagation Applications in Time-Domain). The PSCI projects are joint ef1.

(206) forts between industry and academia, supported by the Swedish Agency for Innovation Systems (Vinnova). The objectives are motivated by industrial demands and a high level of academic research. In the GEMS projects (including the succeeding GEMS2, GEMS3 and SMART) a general software was developed, consisting of hybrid solvers for Maxwell’s equations. The frequency-domain solver combines a state-of-theart multipole boundary integral method with different asymptotic techniques for high-frequency problems. The time-domain solver consists of a coupled structured-unstructured method with an extensive library of modules including subcell objects, such as thin wires and sheets, dispersive materials, absorbing boundary conditions, waveguide modeling, etc.. Both codes are parallelized and have been extensively tested on challenging industrial problems on high performance computers. Several theses [5, 45, 19, 38, 55] including the current one are direct results of the research performed within GEMS. The solvers are in production use in the industry and a company is being founded for the maintenance and further development of the codes. The industrial partners that were involved in the GEMS projects are Ericsson Microwave Systems, Saab Ericsson Space, AerotechTelub (former Saab Avionics), the Swedish Defence Research Agency (FOI), the Fraunhofer-Chalmers Research Centre for Industrial Mathematics (FCC) and the academic partners were KTH and Uppsala University. In the IMPACT project, an inverse scattering method for wave propagation problems was developed, with applications in electromagnetics, ultrasonics and elastodynamics. The problem is solved in the time-domain using a deterministic, gradient-based minimization approach where the gradient is computed using the adjoint solution. This approach has the advantage of covering a range of frequencies in a single computation, and of being generic and applicable to existing solvers for the direct problem. Applications in electromagnetics include determining the characteristics of an unknown material using reflection measurements and model reduction. The electromagnetic approach has been described in [1]. The industrial partners were Aerospatiale (now part of EADS), AerotechTelub, Esaote (Italy) and QSW (UK) and the code providers were KTH and CRS4 (Italy). The work presented in [1] was conducted at KTH. In the following chapters a background to the topics addressed in papers I– V is presented. The three-dimensional FEM-FDTD hybrid solver is described in Chapter 2. The solver is a further development of the two-dimensional hybrid solver discussed in Paper I and has been carefully described in [19]. The primary contributions to the hybrid solver in this thesis consist of the development of a UPML absorbing boundary condition in Paper II and a novel subcell model for thin material sheets in Paper IV. In Chapter 3 time-domain waveguide modeling is addressed. This is the subject of Paper V where a 2.

(207) general approach for time-domain analysis of inhomogeneous waveguides is developed. Inverse scattering is a large subject in electromagnetics, with many different research areas. The background given in Chapter 4 is intended as a brief introduction to the minimization approach used for reconstruction problems in Paper III. In this paper a general form of Maxwell’s equations covering both PML and dispersive media is formulated and the derivation of exact adjoint gradients is described. Further, model reduction using subcell models is discussed and gradients are derived for a combined thin sheet and slot geometry.. 3.

(208)

(209) 2. Time-domain hybrid solver. The most commonly used time domain method for the Maxwell equations is the finite-difference time-domain method (FDTD) [77]. It has several advantages compared to other time domain methods. The simple data structures make it straightforward to implement and facilitates vectorization and parallelization. With respect to time and memory consumption, the FDTD method is an almost optimal second-order method for homogeneous materials. The main drawback of the FDTD method is its inability to accurately model curved objects and small geometrical features. The Cartesian grid leads to a staircase approximation of the geometry and small details are not resolved at all. An efficient approach to avoid staircasing is to use unstructured grids near curved objects and around small geometrical details but revert to structured grids as quickly as possible for the rest of the computational domain. This approach combines the efficiency of structured grids with the flexibility of unstructured grids. An example of such a hybrid grid is shown in Figure 2.1.. Figure 2.1: Hybrid grid for a Clavin antenna. A structured FDTD grid is used for the main part of the antenna and an unstructured tetrahedral grid resolves the slot and monopoles.. Wu and Itoh [76] were first to present a combination of the FDTD method and an implicit finite element (FE) method. They have been followed by several others [51, 43, 67, 78]. A combination of explicit finite volumes (FV) and FDTD have been proposed by Riley and Turner [63]. 5.

(210) The coupling of two stable methods does not necessarily lead to a stable hybrid and many of the early efforts suffered from severe stability problems. An important contribution in this area is [64] where a stable coupling between FEM and FDTD is presented. In this chapter we will give a brief description of the two main components of the time-domain hybrid solver, the FETD and FDTD methods, as well as explain the stable coupling between them. For more detailed descriptions we refer to [70] for the FDTD method, [37] for the finite-element method and [19] for the hybrid FEM-FDTD method.. 2.1 Governing equations The propagation of electromagnetic fields is described by the Maxwell equations, ∂D =∇×H −J (Ampère’s law) (2.1) ∂t ∂B =∇×E (Faraday’s law) (2.2) − ∂t ∇·D =ρ (Gauss’ law) (2.3) ∇·B =0. (Gauss’ law). (2.4). where B is the magnetic flux density in V s/m2 , D is the electric flux density in As/m2 , H is the magnetic field in A/m, E is the electric field in V /m, J is the electric current density in A/m2 and ρ is the electric charge density in As/m3 . The divergence equations (2.3) and (2.4) may be seen as resulting from the first two equations in the following sense. Taking the divergence of (2.1) and (2.2) and changing the order of differentiation we obtain ∂(∇ · D − ρ) =0 ∂t ∂(∇ · B) =0 ∂t. (2.5) (2.6). where we have used the continuity equation for the charge ∇·J +. ∂ρ = 0. ∂t. (2.7). The equations (2.5) and (2.6) show that if the divergence equations are satisfied initially, they will continue to be so for all later times. The fields in Maxwell’s equations are related by so-called constitutive relations. In linear, isotropic and non-dispersive (i.e., the material properties are 6.

(211) independent of the frequency of the field) media they are B = µH. (2.8). D = εE. (2.9). where ε is the permittivity, µ is the permeability. The corresponding properties of vacuum are denoted by ε0 and µ0 . Lossy media are included through the current density by J = σE. (2.10). where σ is the conductivity. Using (2.8) and (2.9) the equations (2.1) and (2.2) may be written as ∂H ∂t ∂E ε ∂t. −µ. = ∇×E. (2.11). = ∇ × H − σE − J s .. (2.12). where J s denotes an applied current source. The solution for a given problem is determined by the above equations together with initial and boundary conditions. A common type of boundary is the perfect electric conductor (PEC) used for metals. In this case the proper boundary condition is given by n×E =0. (PEC).. (2.13). In two dimensions the system (2.11)–(2.12) decouples into two independent systems. The transverse magnetic (TM) mode: ∂Ez ∂y. µ. ∂Hx ∂t. = −. µ. ∂Hy ∂t. =. ∂Ez ∂x. (2.15). ∂Ez ∂t. =. ∂Hy ∂Hx − − σEz − Jzs ∂x ∂y. (2.16). ε. (2.14). 7.

(212) and the transverse electric (TE) mode: ∂Hz − σEx − Jxs ∂y. ε. ∂Ex ∂t. =. ε. ∂Ey ∂t. = −. µ. ∂Hz ∂t. =. ∂Hz − σEy − Jys ∂x. ∂Ex ∂Ey − ∂y ∂x. (2.17) (2.18) (2.19). where the field subscripts denote the respective components in a Cartesian coordinate system.. 2.2 The finite-difference time-domain method The minimization approach in Paper III uses an FDTD scheme for the Maxwell’s equations. In this section an introduction is given to the FDTD method, including modeling of dispersive media in Section 2.2.5 and subcell models in Section 2.2.6. These subjects are further discussed in the context of inverse problems in Paper III. The accurate modeling of electrically large problems may result in computationally very large problems, even in 2D. Therefore, the full minimization algorithm described in Paper III was parallelized using MPI. The implementation is briefly described in Section 2.2.7. Absorbing boundary conditions (ABC) in Section 2.2.4 and plane wave excitation in Section 2.2.3 are important parts of the waveguide modeling discussed in Papers II and V.. 2.2.1. Discretization. FDTD is based on the discretization of the Maxwell curl equations (2.11) and (2.12) using centered finite differences on a uniform Cartesian grid. The resulting scheme is explicit and second-order accurate. The discrete E - and H -fields are staggered in both time and space. In Figure 2.2 the unit cell in the three-dimensional grid is illustrated to the left and the two-dimensional grid for the TM equations (2.14)–(2.16) to the right. 8.

(213) x j=Ny y j=Ny−1 j=2. z. Ez. Hx. Hy. y x. j=1. Hz. Ey. Hx. y. z. Ex x. j=0. Ez i=0. Hy i=1. i=2. i=3 i=Nx−1. i=Nx. Figure 2.2: The FDTD unit cell in 3D (left) and the FDTD grid for the 2D TM case (right).. In an explicit scheme the solution is propagated by time-stepping. We illustrate the FDTD scheme for the TM case. Let the spatial grid index (i, j) correspond to the physical coordinates (i∆x, j∆y) and the temporal index n denote the time t = n∆t. Half indices are used for the components that are centered between the grid points. Using this notation, the time-stepping expressions for the TM case in a source free region become Ez |ni,j+1 − Ez |ni,j n+ 21 n− 12 , (2.20) Hx |i,j+ 1 = Hx |i,j+ 1 − Ci,j+ 12 2 2 ∆y Ez |ni+1,j − Ez |ni,j n+ 21 n− 21 , (2.21) Hy |i+ 1 ,j = Hy |i+ 1 ,j + Ci+ 12 ,j 2 2 ∆x n Ez |n+1 i,j = Da |i,j Ez |i,j   n+ 1 n+ 1 n+ 1 n+ 1 Hx |i,j+2 1 − Hx |i,j−2 1 Hy |i+ 12,j − Hy |i− 12,j 2 2 2 2  − (2.22) + Db |i,j  ∆x ∆y. where Ci,j = ∆t/µi,j , 1 − σi,j ∆t/2εi,j , Da |i,j = 1 + σi,j ∆t/2εi,j ∆t/εi,j . Db |i,j = 1 + σi,j ∆t/2εi,j. (2.23) (2.24) (2.25). The material properties (ε, µ, σ ) are assumed to be constant within each cell. In order to preserve the second order accuracy in inhomogeneous media, interpolation has to used at the interfaces between different media, as discussed in [5] 9.

(214) The electric field in the σEz -term in (2.16) is approximated using the centered time average n+ 21. Ez |i,j. 2.2.2. ≈. n Ez |n+1 i,j + Ez |i,j. 2. .. (2.26). Stability. The use of an explicit time-stepping scheme for a hyperbolic system implies that the time step has to be chosen sufficiently small for the solution to remain bounded. The stability condition in the three-dimensional case (see for example [69]) is: 1 ∆t < (2.27) 1 1 1 cmax (∆x)2 + (∆y) 2 + (∆z)2 where cmax is the maximum value of the phase speed c = √1εµ in the computational domain. Note that the given stability condition does not take into account boundary conditions, dispersive material models, subcell models,etc.. Such modifications of the basic scheme may result in further restrictions on the time step.. 2.2.3. Plane wave excitation. Plane waves are generated by impressing electric and magnetic current sources on a virtual surface, called Huygens’ surface, surrounding the excited region. This surface divides the computational domain into a scattered-field and a total-field region. To see how the excitation is implemented in FDTD, consider the Ez -equation (2.22) in the 2D-TM case for a Huygens’ surface with normal in the x-direction. The electric current source is added to the Ez -field according to n+1 Ez |n+1 i,j = Ez |i,j −. where n+ 21. Jzs |i,j. ˆ· =z. ∆t s n+ 21 J | ε0 z i,j. 1 inc n+ 21 1 n+ 1 ˆ × H inc |i,j 2 = x H | 1 . ∆x ∆x y i+ 2 ,j. (2.28). (2.29). The incident magnetic field H inc is assumed to be a known function. The scaling 1/(∆x) converts the surface current density into an equivalent volume current density. Similarly, magnetic current sources are added to the Hy components on the Huygens’ surface. Due to the staggered grid, the excitation surfaces for the electric and mag10.

(215) netic currents are separated by a half cell. This has to be taken into account when evaluating the current source to avoid perturbation of the incident wave, see [16].. 2.2.4. Absorbing boundary condition (ABC). In order to truncate the computational domain, boundary conditions are needed. For open problems non-reflecting, artificial boundaries have to be introduced to avoid outgoing waves from reentering the problem. Among the most popular ABCs that have been used in FDTD is the Mur ABC [52] which implements the Engquist-Majda boundary condition [24]. It is computationally cheap but in many cases the reflection errors are not sufficiently low. The perfectly matched layer (PML) technique introduced in [8] is a different approach where an outer layer of highly absorbing media is used to attenuate outgoing waves. This approach is computationally expensive but allows very low levels of reflection independent of the angle of incidence of the incident wave. In [66] an anisotropic PML absorber was described, called uniaxial PML (UPML). In this case the PML consists of an anisotropic medium where the permittivity and permeability tensors are of the form     a 0 0 a 0 0     ε = ε0 εr  0 a 0  and µ = µ0 µr  0 a 0  (2.30) 0 0 a−1. 0 0 a−1. where εr and µr are the relative permittivity and permeability and a is chosen σ as a(ω) = 1 + jωε . In order to avoid numerical reflections the conductivity 0 is gradually increased in the PML. Usually, the layer is terminated by an outer PEC wall. σ In [44] the modified attenuation function a(ω) = κ + α+jωε was in0 troduced. This complex frequency shifted PML (CFS-PML) has later been shown to remedy problems of linear growth [7] and further the attenuation of evanescent waves is improved [9].. 2.2.5. Dispersive media. Most real materials are dispersive, i.e., the constitutive relations (2.8)–(2.10) depend on the frequency of the wave. In case of dispersive permittivity the constitutive relation for the electric field may be written as ˆ ˆ ˆ D(ω) = εˆ(ω)E(ω) = ε0 (ε∞ + χ(ω))E(ω). (2.31) 11.

(216) where εˆ is the complex dielectric function, ε∞ is the asymptotic permittivity when ω → ∞ and χ(ω) is called the susceptibility function. In the timedomain the corresponding relation becomes

(217) t D(t) = ε0 ε∞ E(t) + ε0 E(t − τ )χ(τ ) dτ (2.32) 0. where the second term is a convolution. For the FDTD modeling of dispersive media two different approaches have been suggested in the literature. In the auxiliary differential equation (ADE) method [39], equation (2.31) is turned into an ordinary differential equation by an inverse Fourier transform. The resulting equation is discretized using an FDTD approximation and solved for the D -field. The electric field is then obtained from (2.1) which is solved in place of (2.12). In the recursive convolution method [41], the integral in (2.32) is discretized and computed numerically. The approach uses the fact that for most common dispersive material models it is possible to compute the convolution recursively, thereby avoiding a summation over the entire time-history of the field.. 2.2.6. Subcell models. A uniform grid is poor when the grid size is limited by resolution of small geometrical details and not the wavelength. In this case most of the grid may be unnecessarily fine. An alternative to grid refinement is local modifications of the update equations. These subcell models are often derived from the integral form of Maxwell’s curl equations

(218)

(219) ∂E ε + σE · dS = H · dl (2.33) ∂t S C

(220)

(221) ∂H · dS = − µ E · dl. (2.34) ∂t S C The above equations are applied for the discrete fields in the cells containing the subcell geometry. Integration over the cell results in a time-stepping expression where the material coefficients depend on the subcell geometry. In Figure 2.3 we see an example of a thin sheet with a narrow gap, located in a TE FDTD grid.. 2.2.7. Parallelization by domain decomposition. Due to the local property of the FDTD scheme, it is well suited for parallelization using domain decomposition where the computational domain is 12.

(222) (ε, σ, µ) g. d (εs , σs , µ). Ey. Hz Ex. Figure 2.3: Subcell thin sheet with a narrow gap in the FDTD TE grid.. divided along the Cartesian axes, as illustrated for a two-dimensional problem to the left in Figure 2.4. Note that the processor grid does not have to be uniform. To illustrate the parallelization strategy, consider the 2D-TM equations (2.20)–(2.22). The Ez -field is advanced in time using the surrounding Hx - and Hy -fields and the H -fields are correspondingly updated using the surrounding Ez -fields. In Figure 2.4 (to the right) the communication over a processor boundary in the y -direction is illustrated. The following operations take place during one time step ∆t, 1. Processor Pi+1,j sends the Hy -fields to the right of the boundary to Pi,j , 2. Pi,j updates the Ez -fields, t = t + ∆t/2, 3. Pi,j sends the Ez -fields to the left of the boundary to Pi+1,j , 4. Pi+1,j updates the Hy -fields, t = t + ∆t/2.. y processor boundary P(i,j). Pi,j. P(i+1,j). Pi+1,j communication. PML H−field E−field. x. ghost cell. Figure 2.4: Domain decomposition using a Cartesian processor grid (left) and communication over a processor boundary (right).. The efficiency of a parallel FDTD implementation is often very good due to the small amount of data in the communication. As an example, see Figure 2.5 where the parallel 2D TM solver is evaluated. Two different measures were 13.

(223) Figure 2.5: Speed-up and scale-up for the direct solver (left) and speed-up for the different parts of a full minimization (right).. used, scale-up, where the local problem size on each processor is constant and speed-up where the global problem size is constant. In the first case the problem size for each processor was 2500 × 2500 cells and in the second case the total problem size was 2500 × 2500 cells. The results are illustrated to the left in Figure 2.5. The results to the right in Figure 2.5 show the parallel efficiency for the inverse method discussed in Paper III.. 2.3 The finite-element time-domain method The unstructured solver is based on a time-domain finite-element method for the vector wave equation. Taking the time derivative of (2.12) and substituting the magnetic field using (2.11) results in the vector wave equation for the electric field. 1 ∂J s ∂2E ∂E +∇× ∇×E =− , (2.35) ε 2 +σ ∂t ∂t µ ∂t The initial conditions are taken to be homogeneous E=. ∂E = 0 at t = 0 ∂t. (2.36). Depending on the problem, different boundary conditions are applied. We distinguish between Dirichlet boundary conditions where the tangential electric 14.

(224) field is prescribed n × E = g d (r, t),. r ∈ Γd. (2.37). and Neumann boundary conditions where n × (∇ × E) = g n (r, t),. r ∈ Γn .. (2.38). Examples of Dirichlet boundary conditions are perfect electric conductors (PEC) where g d (r, t) = 0 and hybrid boundaries where the FEM domain connects to FDTD. The perfect magnetic conductor is an example of a Neumann boundary condition where g n (r, t) = 0. Impedance boundary conditions and absorbing boundary conditions are also of this type. Defining the space W = H(curl, Ω) = {w : w ∈ L2 (Ω), ∇×w ∈ L2 (Ω)} the weak form of (2.35) is stated as: Find E ∈ W such that for all w ∈ W ,

(225) 2 ∂ E 1 ∂E 2 +σ · w + (∇ × E) · (∇ × w) dΩ = ∂t ∂t µ Ω

(226)

(227) 1 ∂J s − (n × (∇ × E)) · w ds − · wdΩ. (2.39) Γ µ Ω ∂t The corresponding finite element problem is to find E ∈ W h such that (2.39) is satisfied for all w ∈ W h where W h ⊂ W is a suitable finite dimensional subspace. In order to avoid the appearance of spurious solutions, it is important that the finite basis preserves the divergence free nature of the electric field. This property is inherent in the edge or Whitney basis functions which were introduced in [75] and extended to tetrahedra in [53]. Compared to nodal elements, edge elements also have the advantage that tangential continuity over element boundaries is enforced and singularities at edges and corners in the geometry are avoided since the unknowns are located on the edges, see Figure 2.6. Edge basis functions have been developed for most of the common element shapes used in unstructured grids. The edge elements in the current FETD implementation are described in [37] (tetrahedra and hexahedra), [29] (prisms) and [15, 28, 79] (pyramids).. Figure 2.6: Examples of linear edge elements where the unknowns are located at the midpoint of the edges, tetrahedra (left), right angular prism element (middle) and brick element (right).. 15.

(228) Let us denote by N i the i:th edge basis function and expand the electric field in (2.39) according to E=. N . Ei N i. (2.40). i=1. where N is the number of basis functions. Following the Galerkin approach the test functions w are chosen as the basis functions. The equation (2.39) may then be written as the linear system of ordinary differential equations Ne 2 ed E e dE e M +S E =f +K (2.41) dt2 dt e=1 where Ne is the number of elements, M e , K e and S e are the elemental mass and stiffness matrices and E is a vector containing the unknowns Ei (t), i = 1, . . . , N . The local mass and stiffness matrices are given by:

(229) e (M )ij = ε (N i · N j ) dV , (2.42) Ve.

(230). e. (K )ij =. σ (N i · N j ) dV ,. (2.43). 1 (∇ × N i ) · (∇ × N j ) dV µ. (2.44). Ve.

(231). (S e )ij = Ve. where Ve is the volume of the element. Assembling the local matrices gives the system M. d2 E dE + SE = f +K 2 dt dt. (2.45). The vector f = f d +f n +f s contains contributions from boundary conditions and sources. Dirichlet conditions (2.37) give contributions of the form fid = −. j∈Γd. Mij. d2 gjd dgjd + Sij gjd + K ij dt2 dt. (2.46). where gjd denotes the function g d (r, t) evaluated on edge j . Neumann conditions (2.38) are included by

(232) 1 n g · N i ds. fin = − (2.47) µ Γn. 16.

(233) The source terms enter the right hand side by

(234) ∂ s fi = N i · J s dV. ∂t. (2.48). V. The system of ODEs in (2.45) is discretized using the second-order NewmarkBeta method [34] resulting in. ∆t ∆t2 ∆t2 n+1 M+ =2 M− K+ S E S En− 2 4 4. ∆t ∆t2 M− K+ S E n−1 − ∆t2 f |t=n∆t (2.49) 2 4 where ∆t is the time step. The numerical scheme is unconditionally stable due to the symmetry and positive definiteness of the matrices [26].. 2.3.1. Absorbing boundary conditions. In the hybrid FEM-FDTD solver the external boundary is typically located in the structured FDTD region and thus the need for an ABC in the unstructured region is limited. In some cases however, it is necessary to truncate the unstructured domain without causing reflections. For example consider a coaxial waveguide feeding an antenna. In this case the waveguide needs to be resolved by an unstructured FEM grid while the exterior problem is most efficiently modeled by FDTD. To avoid that backscattered fields in the waveguide perturb the solution in the exterior domain an ABC has to be used to truncate the waveguide. A simple form of ABC is given by. ∂E 1 =0 (2.50) n × (∇ × E) + n × n × c ∂t √ where c = 1/ µε is the speed of propagation of the electromagnetic wave. Such a boundary condition enters the FEM problem through contributions of the form (2.47). In case the backscattered field is a plane wave with orthogonal angle of incidence, the above ABC works rather well. However, in many common types of waveguides, such as hollow rectangular and circular waveguides, the propagating field consists of plane waves that are almost tangential to an orthogonal surface. For such cases (2.50) is not sufficiently accurate. Recently, the UPML method has been applied to open problems in the finite-element time-domain method [36], [65]. These implementations are based on tetrahedral discretization of the UPML. However, there are problems of late-time instability. In [65] they use higher order edge elements in 17.

(235) the UPML to reduce the effect of the instability. In Paper II we describe a UPML in a complex frequency-shifted form using prism elements. Although the instability is not completely avoided, our results indicate that the problems are less severe for the prism element PML compared to a tetrahedral PML. The strategy employed in Paper II is to create a truncating UPML by extending the triangulated waveguide port into layers of prism elements with constant height, as illustrated in Figure 2.7. An advantage with this approach is that the grid for the UPML does not have to be created by the grid generator but can easily be generated in a preprocessing step in the solver. In the paper it is shown that reflection errors as low as -70dB can be obtained. Waveguide (tetrahedral elements). UPML (prism elements). inner conductors. different media. 000 111 111 000 000 111 000 111 000 111 000 111 000 111. Figure 2.7: Waveguide truncation using UPML and prism elements.. 2.3.2. Subcell models. In principle, a tetrahedral grid could be used to resolve fine system details. However, in practice, the number of unknowns can be prohibitive and may lead to worse conditioning of the system. Thus, the development of accurate models that characterize the physics of small features without the need for a highly resolved grid is essential. In recent publications subcell models for the FETD method have been developed for thin wires and thin slots [61, 21, 22]. In Paper IV a method for subcell modeling of thin sheets and coatings is described. The approach is based on a representation of the sheet using degenerated prism elements, or shell elements. The advantage of this approach compared to impedance boundary condition (IBC) approximations [62] is the possibility to take the discontinuous normal electric field component into account by introducing additional degrees of freedom. This is similar to what is achieved by the Maloney-Smith thin sheet model for FDTD [47]. Shell elements are frequently used in solid mechanics, see e.g. [35] and 18.

(236) references therein, and have also been applied to eddy current problems [60]. Figure 2.8 illustrates the construction of an edge shell element following the approach in [60]. 3’. z y. n. x. 3 1’. 2’. 1. 2. 3. δ 2. 1. Figure 2.8: A right prism edge element (left) and the corresponding shell element (right).. As an example of the results that can be achieved by the shell element modeling, consider a dielectric sphere with permittivity ε = 4ε0 and radius a covered by a thin (a/30) dielectric coating with ε = 10ε0 . In Figure 2.9 the bistatic RCS using the shell elements are compared with the exact solution for the case ka = 1 where k = 2π/λ. 5. 0. 2. σ [dB/π a ]. −5. −10. −15. −20 Mie series FEM shell −25 0. 20. 40. 60. 80 100 Azimuth (deg). 120. 140. 160. 180. Figure 2.9: Bistatic RCS for the coated dielectric sphere for ka = 1.. 2.4 Hybridization The coupling between the FEM and FDTD domains follows the approach in [64]. In order for the coupling to be stable it is necessary that interpolation is avoided. This may be achieved by the use of pyramids and hexahedra in the 19.

(237) transition from the unstructured tetrahedral grid to the structured FDTD grid. The tetrahedral grid is thus phased into a layer of hexahedral cells identical with the FDTD cells via pyramidal elements, as illustrated in Figure 2.10. Dirichlet boundary for FEM. Dirichlet boundary for FDTD. Tet Hexa. Hexa. Pyra. Tet Tet. Tet. Tet. Hexa. FDTD. FEM+FDTD. FEM. Figure 2.10: The hybrid interface between the FEM and FDTD domains. The shaded region show the overlap between the two domains. Each time step, the fields on the respective boundaries are exchanged between the domains and treated as time-dependent Dirichlet boundary conditions. By using trapezoidal integration in FEM for the hexahedral elements, the FEM and FDTD fields in the overlapped region becomes identical. The resulting interface is stable subject to the FDTD stability condition. In [20] results are presented where the FEM-FDTD hybrid solver is applied to a number of complex problems.. 20.

(238) 3. Waveguides. Waveguides are used to transfer electromagnetic energy, typically at microwave frequencies. Applications include all forms of guided signal transmission, for example the feeding of antennas. We will assume that the cross-section of the waveguide is uniform in the direction of the energy transport. The most common shapes are rectangular and circular cross-sections. The field in the waveguide can only propagate in certain patterns, called waveguide modes, determined by the geometry and media in the cross-section. Assume a waveguide oriented in the z -direction. Let us divide the electric ˆ, field in the waveguide into a tangential and a normal part E = E t + Ez z as illustrated in Figure3.1. The tangential field can be written as a linear Et Ez. Figure 3.1: Decomposition of the field in a waveguide oriented in the zdirection.. combination of forward and backward propagating modes. Vn+ en e−jβn z + Vn− en ejβn z Et =. (3.1). n. where en = en (ω, x, y) denotes the electric field pattern for mode n and βn = βn (ω) is the associated propagation constant. Similarly, the tangential magnetic field may be decomposed into modes, hn . The modes are orthogonal

(239) (em × hn ) · dS = 0, βn = ±βm (3.2) S. and are often normalized to unity, i.e.

(240) (en × h∗n ) · dS = 1.. (3.3). S. 21.

(241) where S denotes the cross-sectional surface. The energy transport in system of waveguides is often characterized by the scattering parameters. The energy is assumed to enter and exit the system through ports, cross-sectional surfaces with a field decomposition according to (3.1). For a thorough description of this concept see [14]. We simply note that S11 denotes the reflection coefficient for mode 1 in port 1 and S12 denotes the transmission of mode 1 to port 2, S11 =. V1− , V1+. S12 =. V2− . V1+. (3.4). The amplitude of mode n, propagating backwards through a port in a lossless waveguide is given by

(242) − Vn = (E t × hn ) · dS (3.5) S. where we have used (3.3). For some simple cases, including homogeneous rectangular and circular waveguides as well as coaxial waveguides, the mode fields are given by analytical expressions. However, for more complex geometries and inhomogeneous waveguides they have to be computed numerically. The frequency domain mode field is a solution of the vector Helmholtz equation for the waveguide cross-section. 1 ∇× ∇ × E − k02 εr E = 0, (3.6) µr where k02 = ω 2 µ0 ε0 . Using a FEM discretization a generalized eigenvalue problem is obtained which may be solved iteratively, for example using the Arnoldi method [6]. In Paper V we describe how waveguide modes computed by a frequencydomain mode solver can be used for excitation and registration in the 3D hybrid FEM-FDTD code. The work was also presented in [2]. The suggested approach relies on Huygens’ surfaces for the excitation and the time-domain excitation is generated by inverse Fourier transform of the mode fields. The following problem was analyzed in Paper II. The geometry consists of a rectangular waveguide of dimensions 0.2m×0.1m where the bottom half of the waveguide is filled with a dielectric with εr = 4 and an interior conductor of width 0.1m is suspended on top of the dielectric. The waveguide is enclosed by a PEC boundary. This type of waveguide is referred to as a shielded microstrip. Figure 3.2 illustrates the geometry and the computed mode field for the fundamental quasi-TEM mode. A short section of the waveguide, of length 0.1m, was analyzed in the 3D time-domain FEM-FDTD solver. The 22.

(243) 0.1. 0 0. 0.2. Figure 3.2: The mode field at f=2GHz for the shielded microstrip.. waveguide was terminated in both ends by 8 cells thick UPMLs. The S11 and S12 parameters were computed and are found in Figure 3.3. In this case we should have S11 = 0 and S12 = 1 (note that the results in the figure are given in dB). The RMS error in the frequency interval 2–3GHz was −33dB for the S11 parameter and −57dB for S12 . The main part of the error in the S11 is due to the excitation while the error in S12 is mainly from the UPML, as discussed further in Paper II. 0. −25. −27.5. −0.005. |S | [dB]. 12. −32.5. 11. |S | [dB]. −30. −0.01. −35. −0.015 −37.5. −40 2. 2.2. 2.4. 2.6 f [GHz]. 2.8. 3. −0.02 2. 2.2. 2.4. 2.6. 2.8. 3. f [GHz]. Figure 3.3: The scattering parameters for the basic quasi-TEM mode of the shielded microstrip, S11 (left) and S12 (right). 23.

(244)

(245) 4. Inverse scattering. incident field. receiver. scattering geometry. scattered field. Figure 4.1: The inverse scattering problem.. For any electromagnetic problem we may distinguish between the fundamental part, consisting of the governing equations, and the data i.e. boundary and initial conditions, sources and physical properties, such as geometry and media. Finding the fields from the governing equations with given data is the direct or forward problem. In the inverse problem the field is partially known, for example from measurements or design objectives, and we wish to determine the data from which this observed field originates. In particular, the investigation of an object by measuring the scattered field outside the object is called inverse scattering. Figure 4.1 illustrates an inverse scattering experiment where an unknown object is illuminated by an incident plane wave and the scattered field is observed in a number of receivers. We distinguish between two common types of inverse scattering problems. The inverse media problem where we wish to determine the material parameters of the scatterer and the inverse obstacle problem where the geometrical shape of a perfectly reflecting body is sought. Here, we will only consider the inverse media problem where the unknown parameters appear in the coefficients of the governing equations. Figure 4.2 illustrates an experimental set-up for reflection measurements at Saab Bofors Dynamics, Linköping. 25.

(246) Figure 4.2: Anechoic chamber used for reflection experiments at Saab Bofors Dynamics, Linköping. (Courtesy of Erik Söderström, AerotechTelub AB.). 4.1 Difficulties in the inverse problem The inverse problem is more difficult to solve than the direct problem mainly due to ill-posedness and non-linearity. A problem is called well-posed if the solution exists, is unique and continuous with respect to the data. If any of these properties is not fulfilled the problem is called ill-posed. Inverse problems are usually ill-posed due to a number of reasons. The inverse scattering problem, where the field is only observed outside the scatterer, has inherently non-unique solution since evanescent fields resulting from lossy media or small geometrical features are not observed. Other common difficulties are: • Lack of data. The observed data are not complete, for example due to limited spatial or spectral information in measurement data. • Noisy data. Experimental data are contaminated with random perturbations. • Unreachable observation data. The desired response in an optimization problem may be non-physical or otherwise unattainable by the direct model. • Inexact method. Numerical approximations may lead to instability and limited resolution of discrete models may result in non-physical solutions. There is a number of ways to stabilize the inverse problem. Usually some kind of regularization is used, taking into account a priori information on the parameters. For example, the solution may be required to stay close to an a 26.

(247) priori known value or by penalizing local variations. The observation data often contains unnecessary or unwanted information and some filtering has to be done to extract the pertinent part, for instance by time gating time-domain data or filtering out high frequency noise components. The frequency of the incident field also affects the properties of the inverse problem. The non-linearity is a serious practical difficulty since it precludes direct solution methods and leads to non-convex minimization. However, in some cases it is possible to linearize the inverse problem. The theoretical analysis of inverse problems has primarily been concerned with the issue of instability, for which a unified treatment using regularization is possible. Most of the available theory is limited to linear inverse problems. In [23] a thorough description of the regularization method is given including convergence results.. 4.2 Methods Analytical solutions of inverse electromagnetic scattering problems are possible only in a few special cases, such as scattering from layered media in one dimension. More interesting problems require approximative techniques. The early methods typically involved approximations of the physics in order to obtain a simplified model for the direct problem. In the Born and Rytov approximations [10] linear relations between the material properties of the scatterer and the scattered field are obtained by assuming that the contrast between the scatterer and the background media is low, so-called weak scattering. The error of the approximation depends on the frequency of the incident wave, see [40] for a discussion of the validity of the approximations. In both cases the scattered field is related to the object through a Fourier integral. The Fourier space is covered by varying the position of transmitter and receiver as well as the frequency and the object is recovered using an inverse transform. This inverse algorithm is called diffraction tomography [17]–[18]. An iterative variant of the Born inverse method, called the Born-iterative method [74], takes into account multiple scattering effects which are omitted using a linear model. In X-ray tomography the WKB approximation is used to linearize the inverse problem, see [12], and the object is reconstructed using the Radon transform [58]. In case the object is perfectly conducting, the physical optics approximation (PO) may be applied to linearize the inverse problem. The resulting problem may be solved in a diffraction tomography manner, see [46] and [68]. However, the PO approximation is feasible only for convex objects and high frequencies. In addition to the severe physical approximations, all of the above methods require great flexibility in the measurement procedure since the inverse 27.

(248) transform typically relies on observation data for many different angles of incidence and many different frequencies. Further, multiple scattering effects, which are important for non-convex geometries with high contrast between background media and scatterer, are lost in the linearization. A more general approach is minimization where the difference between simulated and measured field is minimized by varying the unknown parameters in an iterative optimization algorithm. This method is costly, in particular for higher dimensions, since the direct problem needs to be solved in each iteration. The minimization approach was introduced for acoustic problems in [72]. In most references, deterministic algorithms such as the conjugate gradient and quasi-Newton methods are used but global optimization has also been applied, for applications in electromagnetics, see [25], [49] and [56]. The possibility to cover a large range of frequencies using pulse excitation makes time-domain approaches attractive. For time-domain minimization approaches using FDTD, see for example [31], [4], [59], [13], [54] and [71]. In [57] FDTD is used in a layer-stripping reconstruction algorithm. We also mention inverse scattering using wave-splitting where the field is decomposed into waves traveling in opposite directions. The inverse problem is solved in a layer-stripping manner using invariant embedding as described in [32] and [30]. In Paper III, [1] and [3] deterministic gradient-based minimization is used where the cost function is the difference between estimated and observed fields. The gradient is derived using a Lagrange multiplier formulation for the discrete system. By considering the completely discretized problem we ensure that the resulting gradient is algebraically exact. The numerical scheme for the adjoint equations as well as the gradient expressions are thus determined by the direct scheme and not chosen independently. We show that it is possible to state the Maxwell equations in a general form which is valid in both dispersive media (of Debye and Lorentz type) as well as in the perfectly matched layer (PML) and thereby avoid such interface problems. Subcell models are used in FDTD to avoid resolving fine geometrical details by refining the grid. We show how the gradient for parameters in such discrete models may be derived.. 4.3 Minimization In the minimization approach for inverse scattering the cost function is a measure of the difference between the observed field and the computed field for some trial parameters and an optimizer is applied to improve the solution. Since each evaluation of the cost function requires the solution of the direct problem the computational cost is high, particularly for three-dimensional 28.

(249) problems. The resulting optimization problem is non-linear and may contain many local minima. Regularization using a priori information is therefore often employed to obtain a smoother cost function and limit the parameter space. In non-linear optimization there are two basic choices, global optimization where one attempts to find (approximately) the global minimum using stochastic methods, and local, deterministic optimization where gradient information is used to find a local minima. Global techniques, such as simulated annealing [42], [50] and genetic algorithms [33], [27], often require a large number of cost function evaluations and are therefore considered more costly than local methods. In Paper III we use deterministic, gradient-based optimization where the gradients are obtained by an adjoint formulation. The resulting gradients require a lot of disk space but the computational cost is independent of the number of parameters. The adjoint approach for the parameter identification problem was introduced in [11] and is also described in for example [73].. 4.4 Computing gradients The efficiency of deterministic algorithms relies on accurate gradient information. In our case, the cost function depends implicitly on the parameters through Maxwell’s equations and the gradient may not be computed directly. However, using an adjoint formulation the Fréchet derivative of the cost function may be computed using an additional simulation to solve the adjoint problem. An alternative approach is to approximate the gradient using finite differences. Assume a direct problem of the form L(p)u = f , x ∈ Ω, t ∈ [0, T ],. (4.1). u = 0 at t = 0. where L(p) is a linear differential operator in time and space including boundary conditions and depending on some parameters p. An inverse problem can be stated as: given an observed field uobs , determine p such that the solution to (4.1) coincides with the observed field. Let us define a cost function according to

(250) T

(251) j= 0. u − uobs δ obs dΩdt. (4.2). Ω. where δ obs denotes the Dirac delta function which has the value 1 in the observation points and is 0 everywhere else. The minimization approach to the 29.

(252) Initial guess p = p 0. Direct solution [0,T] u (p) Observed data. u obs Compute residual u (p) u obs Minimize j( p) p = p +dp Adjoint solution [T,0] u *(p). Compute gradient j(p). no. j =0 yes. p* = p. Figure 4.3: The minimization procedure using adjoint gradients.. inverse problem is: minimize j(p) subject to the governing equations (4.1). The gradient of the cost function can be derived using a Lagrangian formulation. The resulting gradient becomes ∂j = ∂pi.

(253) T

(254) < λ, 0. Ω. ∂L u > dΩdt. ∂pi. (4.3). where < ·, · > is the scalar product between two vectors and λ is the solution of the adjoint problem. L∗ λ = − u − uobs δ obs . (4.4) The adjoint operator L∗ satisfies the relation

(255) T

(256) 0. 30. Ω. < L∗ λ, u > dΩdt =.

(257) T

(258) < λ, Lu > dΩdt. 0. Ω. (4.5).

(259) 7 iterations. 15 iterations. 20. 20. 15. 15. 10. 10. 5. 5. 5. 10. 15. 20. 200 iterations. 5. 10. 15. 20. cost function and rel. error. 2. 10 20. rel. error cost function 0. 10 15. −2. 10. 10. −4. 10. 5. −6. 5. 10. 15. 20. 10. 0. 50. 100 150 iteration. 200. 250. Figure 4.4: Reconstruction of the “z”-profile. Note that the integrand in (4.3) is non-zero only in the parts of the domain where the direct operator depends on the parameter, pi . To be able to compute the gradient we need to solve both the direct and adjoint problems. Figure 4.3 illustrates the resulting algorithm. 4.5 Numerical examples In this example the inverse algorithm is applied to the reconstruction of a conductivity profile i.e., εr = εr (x, y). Synthetic observation data is used and the geometry to be reconstructed consists of a “z”-shaped object with permittivity εr = 5.0. The field is sampled in twelve points surrounding the unknown domain. The size of the FDTD grid is 40 by 40 cells the cell size is ∆ = 0.001 m and the permittivity is varied within a 21 by 21 cells square resulting in 441 unknown parameters. The dimensions of the unknown region corresponds to one by one wavelength for the cutoff frequency of the incident Gaussian pulse. The permittivity is initially set to free space. 31.

(260) low frequency. high frequency. 3. 0.06. 2. 0.04. 1. 0.02. 0 0. 10. 0 0. 10. 5. 5 10. ε1. ε3. 0. 8. 6. 6. ε3. 8. ε3. 10. 4. 4. 2. 2 4. ε1. 6. 8. 10. ε1. 10. 2. 5. 5. 10. 2. 4. ε1. ε3. 0. 6. 8. 10. Figure 4.5: The cost function using a high-frequency incident field (left) and low-frequency (right).. In Figure 4.4 the permittivity profile in different minimization iterations is illustrated. We note that the “z”-shape is apparent after just a few iterations. The convergence of each parameter is however quite slow and the accuracy is about three digits after 800 iterations. The frequency spectrum of the incident wave strongly affects the characteristics of the minimization problem. In general a short wavelength gives sharper minima, in the sense that a small deviation from the minima results in a large increase of the cost function. The drawback is that short wavelengths also tend to generate many local minima. Longer wavelengths on the other hand give a smoother cost function with less distinct minima. To illustrate this behavior we studied the cost function for a multilayer problem where the parameters were the relative permittivities in two of the layers, ε1 and ε3 . Two different excitations were used, both Gaussian plane waves with cutoff frequencies of about 1 GHz in the low-frequency case and 30 GHz in the high-frequency case. Figure 4.5 illustrates the cost function for the two excitations. The minimum is located at ε1 = 5,ε3 = 3. The figure indicates that for long wavelengths the convergence will be initially fast and then slow down as the iterations approach the minima. This is preferable if little or nothing is known about the media such that the initial guess may be far away from the correct solution. In case the quality of 32.

(261) the a priori data is better, a short wavelength will give faster convergence. One could also think of using a combination of two excitations, a longer wavelength to get a rough estimate of the solution and then a short wavelength to improve the accuracy. In Paper III successful examples of material reconstruction are presented, where both synthetic data obtained by simulation and measurement data from practical experiments are used. Furthermore, a synthetic reconstruction example is given, where parameters in a thin sheet and slot subcell model are recovered using the inverse method.. 33.

(262)

(263) 5. Acknowledgments. There are a number of people who deserve my gratitude for contributing to this thesis in different ways. First of all I would like to thank my co-authors, Ulf Andersson, Fredrik Edelvik, Lars Eriksson, Christer Johansson, Gunnar Ledfelt and Bo Strand for their contributions to the presented papers. Most of the work has been conducted in close collaboration with the industry within the GEMS, GEMS2, GEMS3 and IMPACT projects. This experience has been very rewarding and in particular I would like to thank Anders Höök and Bo Wästberg, Ericsson Microwave Systems, Stephane Alestra, EADS, Bo Strand and Erik Söderström, AerotechTelub for generously sharing their experience and providing results for various test cases. My advisor Prof. Per Lötstedt and co-advisor Dr. Fredrik Edelvik both have a large part in the successful finish of my Ph.D. studies. Thank you for your guidance and support and not the least for the time spent on proofreading! Part of the work presented in this thesis was carried out at Nada, KTH and I am very grateful for the support of my advisors during that period, Prof. Björn Engquist and Prof. Jesper Oppelstrup. I also thank colleagues and friends at KTH and Uppsala University for support and sharing all kinds of extra-curricular activities. Financial support has been provided by the Parallel and Scientific Computing Institute (PSCI), a Vinnova Competence Center, KTH, Uppsala University and the Swedish Foundation for Strategic Research (SSF) and is gratefully acknowledged. Finally, to Nina for your love and support and to Viktor and Vera for making my home a safe haven, completely void of electromagnetic fields (figuratively speaking, at least).. 35.

(264)

(265) 6. Summary in Swedish: Direkta och inversa metoder för vågledare och spridningsproblem i tidsdomänen. I det moderna samhället används elektromagnetiska vågor inom en rad olika områden. De flesta av oss har ju till exempel en mikrovågsugn, i vilken elektromagnetiska vågor sätter vattenmolekyler i maten i rörelse, och en uppvärmning sker. I mobiltelefoner används elektromagnetiska vågor för trådlös kommunikation, ett brett användningsområde som även innefattar radio, TV, trådlösa telefoner, GPS (Global Positioning System), nätverk för datorer samt en mängd andra tillämpningar. I den medicinska världen används högfrekventa elektromagnetiska vågor både i traditionell röntgenutrustning och för tomografiska undersökningar. Att radar används för identifiering inom det militära området känner många till men det finns även flera civila användningsområden inom området fjärranalys, till exempel vid geologiska och meteorologiska undersökningar. Figur 6.1 visar hur olika typer av elektromagnetiska fält kategoriseras efter frekvens. AM− radio. Kraftledningar. 0. 10. DC 10Hz. 2. 10 1kHz. Lågfrekvens. 4. 10. 6. 1MHz. Radiovågor. Radar. TV. 10. 8. 10. Infrarött. 10. 1GHz. 10. 12. 10 14. 1THz. Mikrovågor. Ultraviolett. 10 16 1PHz. Gammastrålning Röntgen 10. 18. 10. 20. 1EHz. 10 1ZHz. 22. frekvens (Hertz). Joniserande strålning. Figure 6.1: Det elektromagnetiska spektrumet.. Den omfattande användningen av elektromagnetiska fält har lett till ett ökat behov av kunskap om hur de beter sig. Det system av partiella differentialekvationer som matematiskt beskriver vågornas utbredning formulerades redan 1873, av James Clerk Maxwell. Att lösa ekvationerna analytiskt med papper och penna är möjligt bara för mycket enkla problem. Ett annat alternativ är att genom experiment och mätningar ta fram den information som eftersöks. I många fall är detta dock alldeles för kostsamt, både vad gäller tidsåtgång och i pengar. Ett bättre alternativ är numerisk lösning av ekvationerna med hjälp av datorer. Detta innebär oftast att den studerade geometrin delas upp i 37.

No results found