DEVELOPMENT AND APPLICATIONS OF VECTOR GENERALIZED FINITE ELEMENT METHOD IN ELECTROMAGNETICS By Ozgur Tuncer A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Electrical Engineering 2012 ABSTRACT DEVELOPMENT AND APPLICATIONS OF VECTOR GENERALIZED FINITE ELEMENT METHOD IN ELECTROMAGNETICS By Ozgur Tuncer The finite element method (FEM) is widely used in development of today’s electromagnetic technology as it provides accurate and efficient simulation of electromagnetic fields. Recently, techniques such as the vector generalized finite element methods (VGFEM) have been proposed to increase the accuracy and flexibility of the FEM. VGFEM is a general framework that enables inclusion of a large class of basis functions other than interpolatory polynomials in the solution space, and offers several advantages such as: (i) easy implementation of p-adaptivity, (ii) use of mixed orders of polynomial basis functions or use of different basis functions within a simulation domain, and (iii) inclusion of physics in the solution space. Hybridizing VGFEM with FEM can provide a framework for reaping advantages of both methods; resulting in a simulation tool that is highly efficient and accurate. In this dissertation, we study the following topics: (i) A general methodology is developed to adapt the method to polyhedral meshes, and tested for brick and tetrahedral meshes. (ii) Convergence and error characteristics of the method are extensively studied, especially dispersion/phase error characteristics. (iii) Techniques are developed to extend VGFEM to inhomogeneous problems. (iv) Boundary integrals, absorbing boundary conditions, and perfectly matched layers are formulated within the VGFEM framework. (v) VGFEM is hybridized with classical FEM. (vi) Finally, time domain VGFEM is developed for transient electromagnetic problems. To my wife, my daughter, and my parents iii ACKNOWLEDGMENT I would like to express my sincere thanks to several people with whom I have succeeded in completion of my doctoral degree. First, I wish to express my deepest thanks to my advisor, Prof. Shanker Balasubramaniam, for his constant guidance, encouragement, and support during my doctoral study. This dissertation is a result of stimulating discussions that I had with him. More importantly, I would like to thank him for his understanding, enthusiasm, and technical discussions that made me think freely beyond the limits of a dissertation. I would also like to thank each member of my doctoral committee: Prof. Leo Kempel, Prof. Edward Rothwell and Prof. Guowei Wei for their time, constructive comments and advice throughout the Ph.D. examinations. Particularly, I wish to thank Prof. Rothwell for his understanding, supportive attitude to gradate students, and more importantly for being an excellent teacher. I would like to extend my thanks all the teachers who gave me knowledge during my study at Michigan State University. I also wish to gratefully acknowledge the support of the High Performance Computing Center at Michigan State University, the Air Force Office of Scientific Research, and the National Science Foundation. Second, I would like to extend my thanks to all members of Electromagnetic Research Group, especially to my officemates Naveen Nair, Vikram Melapudi, Lu Chuan, Andrew Baczewski, Andrew Pray, Daniel Dault, Collin Meierbachtol, Jose A. Hejase, and Benjamin Crowgey not only for their valuable and fruitful technical discussions but also for their friendship. In particular, I would like to thank Lu Chuan for his invaluable technical support and time that made me experienced in the generalized finite element methods in a very short time. iv Finally, my deepest gratitude goes to my wife and my parents for their endless love and support throughout my life. Especially, I would like to thank my wife, Nevin, for her dedication and patience throughout my doctoral study, and my daughter, Elif Bilge, for brightening my life with her existence in the last moments of my study. v TABLE OF CONTENTS List of Tables . . . . . . . . . . . . . . . List of Figures . . . . . . . . . . . . . . . . . . . viii . . . . . . . . . . . . . . Chapter 1 Introduction . . . . . . . . 1.1 Problem Statement and Motivation 1.2 Background . . . . . . . . . . . . . 1.3 Proposed Approach . . . . . . . . . 1.4 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 2 4 5 Chapter 2 Vector Generalized Finite Element Method . . . . . . . 2.1 Theory of VGFEM . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Statement of the Problem . . . . . . . . . . . . . . . . . . . 2.1.2 Definition of Basis Functions . . . . . . . . . . . . . . . . . . 2.1.3 Partition of Unity Domains and Associated Basis Functions 2.2 Applications of VGFEM . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Waveguide Applications . . . . . . . . . . . . . . . . . . . . 2.2.2 Boundary Integrals . . . . . . . . . . . . . . . . . . . . . . . 2.2.2.1 Differences Between VGFEM-BI and FEM-BI . . . 2.2.2.2 Numerical Results . . . . . . . . . . . . . . . . . . 2.2.3 Absorbing Boundary Conditions . . . . . . . . . . . . . . . . 2.2.4 Perfectly Matched Layers . . . . . . . . . . . . . . . . . . . 2.2.4.1 Additional Basis Functions For Discontinuities . . . 2.2.4.2 Domain Decomposition for Discontinuities . . . . . 2.2.4.3 Numerical Results . . . . . . . . . . . . . . . . . . 2.3 Dispersion Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . 2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 8 8 8 12 18 18 22 25 28 36 39 41 46 48 54 60 67 Chapter 3 Tetrahedral Based VGFEM and its Applications . . . 3.1 Formulation of Tetrahedral Based VGFEM . . . . . . . . . . . . . . 3.1.1 Partition of Unity Domains and Associated Basis Functions 3.1.2 Bilinear Formulation of VGFEM . . . . . . . . . . . . . . . 3.2 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 70 70 75 76 vi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 78 83 90 Chapter 4 Hybrid Vector Generalized Finite Element Methods 4.1 Domain Decomposition for VGFEM . . . . . . . . . . . . . . . . . 4.1.1 Numerical Results . . . . . . . . . . . . . . . . . . . . . . 4.1.1.1 Convergence Analysis . . . . . . . . . . . . . . . 4.1.1.2 Applications . . . . . . . . . . . . . . . . . . . . 4.2 Hybridization of VGFEM with FEM . . . . . . . . . . . . . . . . 4.2.1 Numerical Results . . . . . . . . . . . . . . . . . . . . . . 4.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 94 99 99 104 113 114 129 Chapter 5 VGFEM for Time 5.1 TD-VGFEM Formulation 5.2 Numerical Results . . . . . 5.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 131 134 145 3.3 3.2.1 Convergence Analysis . . . . . . . . . . . . . . . . . 3.2.2 Application of VGFEM to Waveguide Problems . . 3.2.3 Application of VGFEM to Open Domain Problems Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . Domain Electromagnetic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Analysis . . . . . . . . . . . . . . . . . . Chapter 6 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . 146 6.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 6.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 Bibliography . . . . . . . . . . . . . . . vii . . . . . . . . . . . . . . . . . . . 151 LIST OF TABLES Table 2.1 Table 3.1 Eigenvalues of a 1.0 m × 0.4454 m × 0.3 m rectangular cavity filled with the uniaxial anisotropic material . . . . . . . . . . . . . . . . . Eigenvalues for an empty 1.0 m × 0.5 m × 0.75 m rectangular cavity viii 49 77 LIST OF FIGURES Figure 2.1 An illustration of GFEM overlapping PU domains covering the domain Ω. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Figure 2.2 Two different GFEM PU domain schemes. . . . . . . . . . . . . . . 13 Figure 2.3 Definition of the domains of GFEM basis functions. For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this dissertation. . . . . . . . . 15 Geometry of WR-90 waveguide, a = 22.9 mm, b = 10.2 mm, r = 15.24 mm and α = 30◦ . . . . . . . . . . . . . . . . . . . . . . . . . 18 S11 of WR-90 waveguide is simulated using different local approximation functions and compared against MoM data [1] for L = 0mm and L = 25mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 S11 of WR-90 waveguide with L = 25mm is simulated using different PU functions and compared against MoM data [1] . . . . . . . . . . 21 Figure 2.7 Illustration of scattering problem from a cavity-backed aperture. . . 23 Figure 2.8 Illustration of integration steps for singular integrals. . . . . . . . . . 26 Figure 2.9 Backscattered RCS of cavity-backed apertures compared against FEMBI results given in [2] and [3] . . . . . . . . . . . . . . . . . . . . . . 30 Backscattered RCS of the cavity-backed aperture with a = 2.5λ, b = 0.25λ, and c = 0.25λ compared against FEM-BI results [4] . . . . . . 31 Backscattered RCS of an empty cavity with a = 2λ0 , b = 2λ0 , c = 3λ0 forφinc = 0◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Figure 2.4 Figure 2.5 Figure 2.6 Figure 2.10 Figure 2.11 ix Figure 2.12 Figure 2.13 Figure 2.14 Backscattered RCS of the cavity-backed aperture with a = 16.26λ, b = 0.2λ, and c = 0.85λ compared against measurement result at 12 GHz [5]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 Backscattered RCS of a 2.89 in × 2.1 in × 0.057 in cavity is analyzed ˆ at 9.2 GHz for E i = θe−jk.r , and compared against measurement [5] . The cavity is filled with ǫr = 4 . . . . . . . . . . . . . . . . . . . . 34 Backscattered RCS of 1λ0 × 0.25λ0 × 0.25λ0 cavity filled with ǫr = 7 − j1.5 and µr = 1.8 − j0.1 . . . . . . . . . . . . . . . . . . . . . . 35 Backscattered RCS of a cavity with a = 7.32cm, b = 5.2cm, c = 0.158cm over a range of frequency. ǫr = 2.17(1 − j0.001) . . . . . . 36 Backscattered RCS of a partially filled cavity with a = 0.3λ0 , b = 0.1λ0 , c = 0.6λ0 , d = 0.2λ0 , and ǫr = 2 − j2. . . . . . . . . . . . . . 37 Figure 2.17 Scattering from a patch . . . . . . . . . . . . . . . . . . . . . . . . . 39 Figure 2.18 Scattering from a sphere with radius of r = 0.5λ. Spherical ABC surface is placed at rabc = 1.0λ . . . . . . . . . . . . . . . . . . . . . 40 Figure 2.19 Shape functions used for the definition of additional basis functions. 43 Figure 2.20 A vector plot of the additional basis functions . . . . . . . . . . . . 45 Figure 2.21 An illustration of material interfaces and domain decomposition . . 46 Figure 2.22 The normal component of the computed electric field propagating in a 1m × 0.1m × 1m domain filled with ǫr = 4 for z > 0.5m and ǫr = 4 for z < 0.5m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 The normal component of the computed electric field propagating in a 1m×0.1m×1m domain through a slab with a thickness of = 0.25m and dielectric constant of ǫr = 4 . . . . . . . . . . . . . . . . . . . . 52 Figure 2.15 Figure 2.16 Figure 2.23 Figure 2.24 Figure 2.25 Reflection coefficient at the PML interface for the wave propagation in a 1.0 λ0 × 0.5λ0 × 1.0λ0 rectangular waveguide . . . . . . . . . . 53 Reflection coefficient is computed using domain decomposition. . . . 55 x Figure 2.26 Definition of 1D PU function on PU domains separated by nodal distance h for α=1.5. . . . . . . . . . . . . . . . . . . . . . . . . . . 56 Figure 2.27 Overlapping PU domains separated by a nodal distance h in 2D. . . 58 Figure 2.28 h−, and p−convergence of the phase error per wavelength in degrees for scalar GFEM using Legendre polynomials and hat-PU functions. 61 Incident angle dependency of the phase error per wavelength in degrees for scalar GFEM with hat-PU function, h = λ/10 . . . . . . . 62 Figure 2.29 Figure 2.30 h−, and p−convergence and incident angle dependency of the phase error per wavelength in degrees for scalar GFEM with HO-PU function. 64 Figure 2.31 h−, and p−convergence of the phase error per wavelength in degrees for VGFEM with Legendre polynomials and hat-PU function. . . . . 65 h−, and p−convergence and incident angle dependency of the phase error per wavelength in degrees for VGFEM. . . . . . . . . . . . . . 66 An illustration of overlapping PU domains defined on a triangular mesh of the domain Ω. . . . . . . . . . . . . . . . . . . . . . . . . . 70 Illustration of the domains of VGFEM functions, and a plot of the PU function.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 Figure 2.32 Figure 3.1 Figure 3.2 Figure 3.3 h− and p−convergence of VGFEM in 1.0m × 0.5m × 0.75m. . . . . Figure 3.4 S11 of WR-90 S-bend waveguide is simulated using mixed polynomial orders and compared against the MoM data in [1]. . . . . . . . . . . 81 S11 of WR-90 S-bend waveguide is simulated using mixed basis functions and compared against the MoM data in [1]. . . . . . . . . . . . 82 |S11 | of H plane WR-75 waveguide with 180◦ bend is simulated using mixed polynomial orders and compared against the MoM data in [6] 83 Simulated scattering parameters of H-plane T-junction are compared against the reference data in [7]. . . . . . . . . . . . . . . . . . . . . 84 Figure 3.5 Figure 3.6 Figure 3.7 xi 79 Figure 3.8 Bistatic RCS of a 0.3m × 0.3m PEC plate. The plate is illuminated by the electric field E = xejk0 z . . . . . . . . . . . . . . . . . . . . ˆ 85 Bistatic RCS of a PEC sphere with radius of r = 0.5λ. The sphere is illuminated by the electric field E = xejk0 z . . . . . . . . . . . . . . ˆ 86 Backscattered RCS of the cavity-backed aperture with a = 0.7λ, b = 0.1λ, and c = 1.73λ compared with FEM-BI results [2] . . . . . . . . 87 Backscattered RCS of the cavity-backed aperture with a = 2.5λ, b = 0.25λ, and c = 0.25λ compared with FEM-BI results [4]. . . . . . . . 88 Backscattered RCS of the cavity-backed aperture with a = 16.26λ, b = 0.2λ, and c = 0.85λ compared against measurement result at 12 GHz [5]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 Backscattered RCS of the cavity-backed aperture with a = 2λ, b = 2λ, and c = 10λ compared against FEM-BI data given in [3] and HFSS data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 Figure 4.1 An illustration of material interfaces and domain decomposition . . 95 Figure 4.2 Effect of coupling coefficient βtc on convergence. The problems is T E10 mode propagation in a rectangular waveguide with the dimensions of a × 0.5a × 2a. . . . . . . . . . . . . . . . . . . . . . . . . . 100 Figure 4.3 Convergence of DD-VGFEM for different number of partitions. The problems is T E10 mode propagation in a rectangular waveguide with the dimensions of 1m × 0.5m × 2m. . . . . . . . . . . . . . . . . . . 101 Figure 4.4 Effect of DG penalty factor on the wave approximation. 1m × 0.5m × 2m rectangular waveguide is partitioned along z axis . . . . . . . . . 102 Figure 4.5 h− and p− convergence of VGFEM with DG transmission conditions for T E10 mode propagation in a 1.0λ0 × 0.5λ0 × 2.0λ0 rectangular waveguide. Problem domain is partitioned into 4 sub-domains along longitudinal dimension. . . . . . . . . . . . . . . . . . . . . . . . . . 103 Figure 4.6 h− and p− convergence of DG-VGFEM for an inhomogeneous problem. Plane wave with θinc = π/6 and φinc = 0 is propagating in a 1m × 0.1m × 1m domain with ǫr = 4 for z > 0m. . . . . . . . . . . 104 Figure 3.9 Figure 3.10 Figure 3.11 Figure 3.12 Figure 3.13 xii Figure 4.7 VGFEM with DDM transmission conditions applied to the wave propagation in WR-90 S-bend waveguide. MoM result is from [1] . . 105 Figure 4.8 Application of DD-VFEM to scattering from a dielectric object in a rectangular waveguide.The waveguide dimensions are a = 2b, d = 0.8b, c = 3d, and b = 10mm, and relative dielectric constant is ǫr = 6 107 Figure 4.9 VGFEM with DG transmission conditions applied to the wave propagation in WR-90 S-bend waveguide. Legendre polynomials with p = 1 local approximation is used. MoM result is from [1] . . . . . . . . . 108 Figure 4.10 VGFEM with DG transmission conditions applied to the wave scattering in WR-75 waveguide with 180◦ bend. MoM result is from [6] 109 Figure 4.11 Simulated scattering parameters of H-plane T-junction are compared against the reference data in [7]. . . . . . . . . . . . . . . . . . . . . 110 Figure 4.12 Use of mixed basis functions within a simulation in VGFEM with DG transmission conditions. Problem is scattering in a WR-90 Sbend waveguide. MoM result is from [1] . . . . . . . . . . . . . . . . 111 Figure 4.13 VGFEM with DG transmission conditions applied to wave scattering problem in a rectangular waveguide with a dielectric obstacle in it. Relative dielectric constant is ǫr = 6. . . . . . . . . . . . . . . . . . 112 Figure 4.14 Illustrations of domain decomposition for FE-VGFEM . . . . . . . . 113 Figure 4.15 T E10 mode propagation in a 1.0 m × 0.5 m × 0.75 m waveguide is analyzed using FE-VGFEM. FEM is used in z = [0, 0.4m] and VGFEM is used in z = [0.4m, 0.75m] for this study . . . . . . . . . . 115 Figure 4.16 Plane wave propagating in a 1.0 m × 0.5 m × 0.75 m domain partially filled with ǫr = 4 . VGFEM is used in dielectric region z = [0, 0.4m] and FEM is used in z = [0.4m, 0.75m] . . . . . . . . . . . . . . . . . 116 Figure 4.17 Scattering parameters for WR-75 -circular waveguide junction. Measured results are from [8]. . . . . . . . . . . . . . . . . . . . . . . . . 118 Figure 4.18 Reflection coefficient for a rectangular waveguide with a dielectric obstacle in it. Relative dielectric constant is ǫr = 6. . . . . . . . . . 119 xiii Figure 4.19 Bistatic RCS of a dielectric sphere with ǫr = 4 and ka = 1. FEM and VGFEM regions are depicted. The sphere is illuminated by the electric field E = xe−jk0 z . . . . . . . . . . . . . . . . . . . . . . . 120 ˆ Figure 4.20 Bistatic RCS of a 0.3m×0.3m PEC plate. FEM and VGFEM regions are depicted. The plate is illuminated by the electric field E = xejk0 z 122 ˆ Figure 4.21 Bistatic RCS of a PEC patch placed on a PEC backed dielectric substrate with ǫr = 2.17 at f0 = 300MHz. Patch size is 0.27m × 0.19m and substrate size is 0.44m × 0.44m × 0.05m. The geometry is illuminated by the electric field E = xejk0 z . . . . . . . . . . . . 123 ˆ Figure 4.22 Bistatic RCS of a PEC sphere with radius of r = 0.5λ. FEM and VGFEM regions are depicted. The sphere is illuminated by the electric field E = xejk0 z . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 ˆ Figure 4.23 Bistatic RCS of a PEC cube with the edge length of a = 0.755λ. FEM and VGFEM regions are depicted. The sphere is illuminated by the electric field E = y e−jk0 z . . . . . . . . . . . . . . . . . . . . 126 ˆ Figure 4.24 Bistatic RCS of a PEC cube with the edge length of a = 0.755λ. FEM and VGFEM regions are depicted. The sphere is illuminated by the electric field E = y e−jk0 z . . . . . . . . . . . . . . . . . . . . 128 ˆ Figure 5.1 Electric field E = z T (t + y · r/c0 ) is observed at x = 0.5m, y = ˆ ˆ 0.25m, z = 0.375m in 1.0m × 0.5m × 0.75m domain. TD-VGFEM with Neumann BC, p = 0 and h = 0.225m is used. . . . . . . . . . . 135 Figure 5.2 Electric field E = z T (t + y · r/c0 ) is observed at x = 0.5m, y = ˆ ˆ 0.25m, z = 0.375m in 1.0m × 0.5m × 0.75m domain. TD-VGFEM with Dirichlet BC, p = 0 and h = 0.225m is used. . . . . . . . . . . 136 Figure 5.3 Temporal convergence, and h− and p−convergence of TD-VGFEM for Impedance BC. Electric field E = z T (t + y · r/c0 ) propagating in ˆ ˆ 1.0m × 0.5m × 0.75m domain is considered for this study. . . . . . . 137 Figure 5.4 h−, and p− convergence of DG-VGFEM and FE-VGFEM for the ˆ ˆ ˆ ˆ normally incident electric field E = θT (t−k·r/c) with k = r , θinc = 0 and φinc = 0 in 1.0m × 0.5m × 1.0m. . . . . . . . . . . . . . . . . . 139 xiv Figure 5.5 h−, and p− convergence of DG-VGFEM and FE-VGFEM for the ˆ ˆ ˆ obliquely incident electric field E = θT (t − k · r/c) with k = r , ˆ θinc = π/4 and φinc = 0 in 1.0m × 0.5m × 1.0m. . . . . . . . . . . 140 Figure 5.6 h−, and p− convergence of DG-VGFEM and FE-VGFEM for a problem with material variation. ǫr = 4 is used in the domain for z > 0.5m ˆ ˆ of the domain 1.0m×0.5m×1.0m. The electric field E = θT (t−k·r/c) ˆ ˆ with k = r , θinc = π/4 and φinc = 0 is obliquely incident. . . . . . . 141 Figure 5.7 Time signal computed at r = 0.667m, θ = π/4, and φ = 0 for radiation from a dipole. . . . . . . . . . . . . . . . . . . . . . . . . . 143 Figure 5.8 Bistatic RCS of a dielectric sphere with ǫr = 4 and a = 1/(2π). FEM and VGFEM regions are depicted. . . . . . . . . . . . . . . . . . . . 144 Figure 5.9 Scattering from a 0.3m × 0.3m plate at 300MHz for Einc = xT (t + ˆ inc = φinc = φ = 0◦ . . . . . . . . . . . . . . . . . . . . . . 145 z · r/c), θ ˆ xv Chapter 1 Introduction 1.1 Problem Statement and Motivation Computational electromagnetics (CEM) is becoming critical in the developments of today’s technology as the technology shifts to complex microwave and terahertz systems. Accurate and efficient simulations of such systems are very important, especially in the near field, and this depends on the robustness and the accuracy of the computational tools. Although the conventional differential equation based methods such as the finite element method (FEM) and the finite difference method (FDM) have been successfully applied for the solution of moderate electrical size problems, they suffer from some computational limitations such as demand on computational resources and unsatisfactory convergence of iterative solvers. In order to overcome these problems, domain decomposition methods (DDM) [9, 10], and discontinuous Galerkin methods (DGM) [11, 12] have been introduced. DDM enables the solutions of electrically large problems thanks to its divide-and-conquer philosophy. Likewise, DGM enables mesh flexibility and reduces the computational burden in parallelization. 1 This dissertation proposes an alternative approach for the accurate solution to electromagnetic problems as opposed to those provided by higher order basis functions [9–17]; This approach is founded on the use and modification of the Vector Generalized Finite Element Method (VGFEM) in concert with other techniques. VGFEM belongs to the class of the partition of unity methods [18]. The overlapping partition unity framework of the method enables the inclusion of large class of basis functions other than interpolatory polynomials in the solution space. Additionally, it offers other features such as easy implementation of p-adaptivity, use of mixed orders of polynomial basis or use of different basis functions within a simulation domain, and inclusion of physics in the solution space. These features of VGFEM and its hybridization with other differential equation methods can relax the requirement of dense and conformal meshing due to some geometrical complexities, provide better approximation of electromagnetic fields and, perhaps, more efficient EM simulations. This dissertation will demonstrate these features of the method. 1.2 Background The finite element method has been used to solve wide range of scalar and vector electromagnetic problems for decades [13]. FEM is typically based on the space of functions that are defined on an underlying tessellation. These basis functions are typically interpolatory polynomial functions. Both the basis functions and their derivatives satisfy specific conditions at boundary interfaces [13]. For instance, Whitney elements are tangentially continuous across inter-element boundaries. FEM has seen continuous development in order to analyze/solve practical problems both accurately and efficiently. In addition, there has always been an interest in developing a framework to enrich the underlying approximation space. Developing 2 the means to do so can be useful in many different scenarios. They may be tailored to better capture electromagnetic fields in non-Lipschitzian domains; they can possibly permit numerical basis functions; and, they can enable a smaller transition regime in open domain problems and possibly provide a more seamless transition in domain decomposition approaches. Over the years, many methods have been developed including hp-cloud method [15], element-free galerkin method [16], and generalized finite element method (GFEM) [19]. The generalized finite element methods based on the partition of unity (PU) approach has been applied to solve various engineering problems [17, 20, 21]. In electromagnetics, the method was first developed for the solution of scalar Helmholtz equations [17, 19, 21–24], and then extended to the solution of vector electromagnetic problems [25–28]. In [25], a complete and convergent vector basis set was developed along with a methodology for obtaining linearly independent set. The vector generalized finite element was initially proposed as a class of meshless methods. The meshless framework of the method requires the definition of overlapping partition of unity domains on which basis functions are defined. Typically, problem domain is sampled and these overlapping domains are defined around these samples such that the problem domain is completely covered. Although the meshless scheme of VGFEM using canonical PU domains is flexible, one needs to provide information regarding boundaries of domains/inhomogeneities. The need for surface information proves to be a bottleneck in using the method for a larger class of problems. Moreover, in a generic VGFEM procedure, two sets of functions need to be defined; one that forms a partition of unity and another that provides higher order approximation within the PU domain. The conventional vector basis functions developed in [25] are based on functions that are continuous across boundaries. 3 However, it is well known that such function spaces are not appropriate to model electromagnetic fields as material variation appears; here, it is necessary that the function spaces used be normally discontinuous and tangentially continuous across a boundary. In addition, it is required that the normal derivatives of the tangential components be discontinuous as well. This was achieved in [25] by defining an auxiliary set of functions that satisfied all these conditions, where it was assumed the interface was planar and crossing the entire PU domain. Although definition of such additional basis functions is possible for the interfaces that are canonical, it is not a generalizable for arbitrary shaped interfaces. 1.3 Proposed Approach This dissertation addresses several of these issues, and presents solutions through • development of methodologies to use VGFEM with arbitrarily shaped PU domains, • development of a method to understand dispersion characteristics of VGFEM, • development of techniques for complex inhomogenous problems (material discontinuities), • formulations of absorbing boundary conditions and perfectly matched layers in the VGFEM framework for domain truncation, • hybridization of the method with boundary integrals, • hybridization of the method with other differential methods such as FEM, DDM and DGM, 4 • development of time domain VGFEM (TD-VGFEM) for the analysis of transient EM problems, • and applications of all developed techniques to wide range of EM problems. 1.4 Outline The outline of this dissertation is as follows: In Chapter 2, the framework of VGFEM is adapted to polyhedral mesh elements. Techniques are developed for the construction of partition of unity domains and basis functions on these elements. The framework of VGFEM is tested via applications to various open and closed domain problems. A hybrid VGFEMBI technique is developed to study scattering from cavity backed apertures. Bilinear formulations of the method with the perfectly matched layers (PML), absorbing boundary conditions, and boundary integrals are presented in detail. Methodologies to include discontinuities in the VGFEM solution space are developed for the new framework of VGFEM. Finally, dispersion characteristics of the method is studied rigorously by developing a semianalytic technique following the traditional dispersion analysis. All these developments are validated through a large number of simulations of various problems. In Chapter 3, the framework of VGFEM is adapted tetrahedral meshes. Techniques to define the PU domains and basis functions on this tessellation are developed. In addition, VGFEM is applied to wide range of open and closed domain problems using the capabilities of VGFEM, especially use of mixed polynomial orders and mixed basis functions in the same simulation domain. In Chapter 4, VGFEM is hybridized with FEM based methods such as domain decomposition methods and discontinuous Galerkin methods. A technique to define the framework 5 of VGFEM in partitioned domains will be developed. Domain decomposition technique will also be utilized to handle discontinuities at material interfaces. Then, the formulations for the coupling of VGFEM with FEM, DDM or DGM are presented using appropriate transmission conditions at the interfaces. Finally, the proposed hybrid techniques are validated via the analysis of various open and closed domain electromagnetic problems. In Chapter 5, a time domain VGFEM (TD-VGFEM) is developed to utilize all advancements and capabilities for solving time domain EM problems. TD-VGFEM will be formulated using the Newmark time stepping scheme. Temporal and spatial convergences of the method are studied in details. Finally, the contribution of this dissertation is summarized, and future research topics of this method are listed in Chapter 6. 6 Chapter 2 Vector Generalized Finite Element Method In this Chapter, a general framework of VGFEM that is built on polyhedral meshes is presented. Methodologies to construct partition of unity domains and basis functions in polyhedral meshes are developed. These methodologies are formulated for brick elements in Section 2.1. The vector basis functions proposed in [25] is used through this development, and necessary modifications are realized for defining basis functions in polyhedrons. These developments are validated by applying the method to various open and closed domain problems in Section 2.2. Open domains are truncated using the perfectly matched layers (PML), absorbing boundary conditions and boundary integrals, and techniques for such truncations are developed for VGFEM. In addition, the dispersion and convergence characteristics of VGFEM are investigated. Contrary to classical FEM, the presence of overlapping in VGFEM domains precludes simple analysis. Thus, a technique to study the dispersion in VGFEM is introduced in Section 2.3. Finally, the work is summarized in Section 2.4. 7 2.1 2.1.1 Theory of VGFEM Statement of the Problem Consider a linear and homogenous domain Ω whose boundary is denoted by ∂Ω := Γ = i Γi . Interior to the domain, function u(r) satisfies the vector Helmholtz equation ∇ × ∇ × u(r) − k 2 u(r) = f(r) Bi {u(r)} = gi (r), ∀ r ∈ Γi , (2.1) where u(r) is used to denote either the electric field or the magnetic field, f(r) is the impressed source, k is the wavenumber, Bi is a differential operator, and gi (r) is the function imposed on Γi . In the above equation, it is assumed that r ∈ R3 . To solve this problem using VGFEM, we need to develop (i) an appropriate partitioning of the problem domain, and (ii) basis functions that are defined on these domains. These are elucidated next. 2.1.2 Definition of Basis Functions Consider the domain Ω depicted in Fig. 2.1. The domain is covered by sub-domains Ωi defined around N nodes such that they satisfy the point-wise overlap condition ∀ r ∈ Ω card i|r ∈ Ωi ≤ M, (2.2) where M is the maximum number of sub-domains covering the point r. These sub-domains are called partition of unity (PU) domains. On each PU domain Ωi , a PU function ψi (r) is 8 hp i Figure 2.1: An illustration of GFEM overlapping PU domains covering the domain Ω. defined such that N ψi (r) = 1 on Ω. (2.3) i A space of vector basis function can then be defined as [25] ˆ ˆ V = span{ψi (r)∇ × cvi,m (r) , ψi (r)∇ × ∇ × cvi,n (r) }, (2.4) ˆ where vi (r) denotes local approximation space, c is the pilot vector, and m and n are used to identify the function space vi (r) for each basis function appearing in local Helmholtz decomposition. Details of this basis set, its convergence properties and completeness, and the manner in which this can be modified to avoid spurious modes can be found in [25]. In addition, a technique by which one can obtain an independent set is also presented in [25] as well as other candidate function spaces. The PU function can be constructed using any Lipschitz continuous function Wi (r) by 9 ψi (r) = provided that Wi (r) ∀ Ωj ∩ Ωi j Wj (r) (2.5) j Wj (r) > 0 and Wi (r) smoothly vanishes at the PU domain boundaries, where j denotes the index of the domain covering the point r for r ∈ Ωi , [18]. A simple function that can be used to construct PU function in one dimension (1D) is a hat function,    1 + 2(x−xi ) if x ∈ (x − hp , x ]   i 2 i hp    hp 2(x−xi ) Wi (x) = if x ∈ (xi , xi + 2 )  1− hp       0 elsewhere, (2.6) where xi denotes the center of the PU domain i, and hp represents the size of the PU domain. It is apparent that hat-based PU function has discontinuous derivatives. Alternatively, PU functions that have higher order continuous derivatives can be constructed using piecewise polynomials such as    (1 + 2(x−xi ) )2 (1 − 4(x−xi ) ) if x ∈ (x − hp , x ]   i 2 i hp hp    hp 4(x−xi ) 2(x−xi ) 2 Wi (x) = ) (1 + ) if x ∈ (xi , xi + 2 ) (1 −  hp hp       0 elsewhere. (2.7) Three-dimensional (3D) PU function is a product of 1D PU functions. The effects of using different types of PU functions on dispersion will be investigated later. Local approximation functions are defined on each local domain Ωi , and chosen from p a space of functions vi (r) ∈ span vi (r) , where p is the order of local approximation function. Thus, the approximation to the unknown in (2.1) can be written as 10   Pi,1 Pi,2   p p p p ˆ u(r) = ˜ ψi (r)  ai ∇ × cvi,m (r) + bi ∇ × ∇ × ˆvi,n (r)  , c p=1 p=1 i N (2.8) p p where ai and bi are unknown coefficients, and Pi,1 and Pi,2 are the number of basis functions for each type on the PU domain i. Total number of basis functions is Nb = Pi,1 + Pi,2 . i The unknown coefficients are solved from the linear system of equations that is obtained by applying Galerkin’s method to (2.1) subject to boundary conditions ˆ u(r) × n = gd (r) f or r ∈ Γd ˆ n × [∇ × u(r)] = gn1 (r) f or r ∈ Γn ˆ u(r) · n = gn2 (r) f or r ∈ Γn , (2.9) ˆ where n is the unit outward vector normal to the boundary Γ and Γ = Γd ∪ Γn . Specifically, the bilinear form for (2.1) is written as Ω (∇ × w) · (∇ × u) dΩ − k 2 Γd β1 Γd w · u dΩ + β Ω (∇ · w)(∇ · u) dΩ + ˆ ˆ [(w × n) · (∇ × u) + (u × n) · (∇ × w)] dΓ + ˆ ˆ (w × n) · (u × n) dΓ + β2 Γd Γn Ω gd · (∇ × w) dΓ + β1 w · gn1 dΓ + β2 Γn Γn ˆ ˆ (w · n)(u · n) dΓ = Γd ˆ (w × n) · gd dΓ − ˆ (w · n)gn2 dΓ + Ω f · w dΩ, where the constants β, β1 and β2 contribute uniform convergence [25]. 11 (2.10) 2.1.3 Partition of Unity Domains and Associated Basis Functions The definition of basis function, in VGFEM, is intimately tied to the support of the function. As was apparent from the description thus far, classical VGFEM relies on using an open cover of the domain Ω, i.e., i Ωi ⊃ Ω. This is advantageous from a couple of perspectives; (i) it enables a meshless description of the computational domain, and (ii) it permits the use of canonical elements (such as spheres, cubes, etc.) to form the open cover [25, 29]. The latter permits the use of simple functions to construct the PU function. However, the principal drawbacks of this construction are the need to define exterior surface and surfaces that separate piecewise homogeneous domains, and to determine their intersection with canonical PU domains. This is illustrated in Mesh-A in Fig. 2.2(a). This is possible by appropriately identifying these surfaces within each PU domain, expressing it using a set of polynomials, and then fitting this surface using a least square algorithm. This can be a significant bottleneck on two counts. It is difficult to utilize this sequence for geometrically complex problems, and it is difficult to enable analysis of practical problems due to lack of availability of such GFEM meshing schemes. The development here is intended to make this transition easier. In this section, we seek to develop a technique such that existing meshing information (such as that illustrated in Fig. 2.2(b)) can be used to extract the underlying GFEM mesh. Thus, we will develop a framework for construction of PU domains using readily available quadrilateral meshes as illustrated in Fig. 2.2(b). Note, a 2D quadrilateral mesh is shown in the figure to better understand the framework. Extension to forming of PU domains with 3D elements and other types of elements is straightforward. Specifically, the PU domain is 12 (a) Mesh A: Nodes are uniformly distributed. Each PU domain is associated with a node. p=7 p=8 sp=3 p=5 p=4 sp=1 p=1 p=2 patch Sub-patch 1 1 1, 2 2 p=9 3 sp=4 2 1, 3 4 1, 2, 3, 4 5 2, 4 6 sp=2 p=6 3 7 3, 4 8 4 p=3 9 (b) Mesh B: Pathes are constructed from the geometry mesh. Each PU domain is associated with a node, and PU domain domain comprises of the sub-domains connected to this node. Figure 2.2: Two different GFEM PU domain schemes. 13 constructed from the sub-domains that are connected to a node, Ωi = j ¯ Ωj if i ∈ nj (1 : 4), (2.11) ¯ where Ωi denotes the PU domain of node i, Ωj denotes the domain of the sub-domain j, and nj is an array having the indexes of nodes of the sub-domain j. The definition of PU function goes hand-in-hand with the redefinition of the PU domains. The definition of PU function on an arbitrary shape PU domain is difficult due to the condition that the PU function must vanish at the boundaries of the domain. In order to overcome this problem, each PU domain, which is defined around a node and composed of four sub-domains sharing this node, is mapped onto a canonical domain of [−1, 3] × [−1, 3] as shown in Fig. 2.3(a). Each sub-domain belongs to the set T1 = {[−1, 1] × [−1, 1], [−1, 1] × [1, 3], [1, 3] × [−1, 1], [1, 3] × [1, 3]}. (2.12) The PU function then can be defined in this reference domain as depicted in Fig. 2.3(b). The mapping of the PU function to physical domain is performed using the transformation 4 x= n=1 4 y= n=1 1 [1 + (ξn + ξs )(ξ + ξs )][1 + (ηn + ηs )(η + ηs )]xn 4 (2.13) 1 [1 + (ξn + ξs )(ξ + ξs )][1 + (ηn + ηs )(η + ηs )]yn , 4 where (xn , yn ) and (ξn , ηn ) denote node positions in the physical and reference sub-domains respectively, and shifting constants (ξs , ηs ) belong to the set 14 (a) Geometry map 3 1 1 0 1 3 w 1 1 −1 −1 u (b) PU hat function Figure 2.3: Definition of the domains of GFEM basis functions. For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this dissertation. 15 T2 = {(0, 0), (0, −2), (−2, 0), (−2, −2)}, (2.14) which corresponds to the set T1 . If the PU function is continuous and differentiable in this domain, one can readily prove that those properties translate to the physical domain as well. Next, the local approximation functions can be defined in the physical domain as illustrated in Fig. 2.3. As shown in this figure, it is important to realize that what we are interested in is the definition of the basis function ui (r) as a product of the PU and the local approximation functions. Consequently, the domain of the approximation function does not need to be identical to that of the PU domain. In fact, it is necessary for the domain of the PU function Ωi at any node i to be a subset of the approximation function domain Ωa,i at that node, i.e. Ωi ⊂ Ωa,i . This is due to the fact that ψi (r) = 0, ∀ r ∈ {Ωa,i /Ωi }, which ensures the basis function ψi (r) vi (r) = 0, ∀ r ∈ {Ωa,i /Ωi }. For instance, if Legendre polynomials are used, a rectangular domain is needed for the definition of the functions. Thus, we construct the rectangular PU domain by centring it at the PU domain node position (xi , yi) and finding the size of the PU domain such that it encloses the physical PU domain domain. For example, the size of rectangle for linear quadrilateral elements can be easily determined by     lx = 2 max  |xi − xj,n | and ly = 2 max  |yi − yj,n | , j,n j,n (2.15) where (j, n) ∈ {1, 2, 3, 4} are used to denote the nth node of j th sub-domain. Likewise, if exponential functions are used, it makes sense to center them in a circle about (xi , yi ) with 16 a radius of (lx /2)2 + (ly /2)2 , where lx and ly are given in (2.15). Next, it is important to understand the computational complexity of the method presented here. It is assumed that all examples are in 3D, a uniform brick meshing is used, and the number of basis is identical for each PU domain (reduction to 2D is trivial). To this end, assuming that the approximation is of order p, then the local approximation functions in (2.8) are of order (p+1) and (p+2) respectively. It follows that the upper bound on the total number of vector basis functions is (p+3)(p+4)(2p+7)/6. Of course, this number is reduced considerably once the dependent vector basis functions are removed from the approximation space; in our experience, by approximately a factor of two. For instance, if the Legendre polynomials are used for local approximation function, the total number of vector basis functions in each PU domain is Nb = 11 for p = 1 , and Nb = 26 for p = 2. Then, the total number of unknowns is NNb , where N is the number of PU domains. Further, in the scheme described above, each PU domain overlaps M number of PU domains, where the maximum value of M is 27 for a non-boundary PU domain. Consequently, it can be inferred that the number of non-zero matrix entries per row is Nb M, and the maximum number of non-zero matrix entries in a row is 27Nb . As a point of comparison, assuming that a cuboid is discretized uniformly with N nodes forming brick elements, then the number unknowns for higher order 2 FEM is Nf em = 3p2 (p + 1)(Nd − 1)3 + 6p(p + 1)(Nd − 1)2 Nd + 3(p + 1)(Nd − 1)Nd , where Nd = N 1/3 . Likewise, it can be easily proven that the upper bound on the number of non-zero matrix entries per row for FEM is 3(p + 1)[4(p + 2)2 − 1]. 17 2.2 Applications of VGFEM In this subsection, we apply VGFEM to various closed and open domain EM problems, and present additional theoretical developments that are necessary for these applications. To solve open domain problems, we propose methodologies to integrate VGFEM with domain truncation methods such as the perfectly matched layers (PML), absorbing boundary conditions and boundary integrals. In this thesis, we use the anisotropic interpretation of the PML; therefore, we reformulate the VGFEM for anisotropic medium. 2.2.1 Waveguide Applications Here, the new framework of the VGFEM is validated through the analysis of wave propagation in an S-bend rectangular hollow waveguide. Fig. 2.4 shows the geometry of cascaded H-plane bends in a WR-90 waveguide with a = 22.9 mm, b = 10.2 mm, r = 15.24 mm and α = 30◦ . Note, the length of the straight section connecting the cascaded curved Hplane bends is variable, and the geometry of the curved sections is constructed with linear elements, which obviously contributes to some error. While fine meshing is performed to Figure 2.4: Geometry of WR-90 waveguide, a = 22.9 mm, b = 10.2 mm, r = 15.24 mm and α = 30◦ 18 reduce the error from geometry representation at the curved regions, coarse meshing is done for the straight waveguide sections, the lengths of which are chosen to be 3a. The boundary conditions imposed at the input and output ports of the waveguide are respectively given by [13], ˆ z −ˆ × ∇ × u + jkz10 z × ˆ × u = −2jkz10 E10 z (2.16) ˆ ˆ z z × ∇ × u + jkz10 z × ˆ × u = 0, (2.17) and −jkz10 z ˆ where E10 = −j sin( π x) e y and kz10 = a k 2 − ( π )2 . Dirichlet boundary cona dition is imposed at the metallic boundaries using Nitsche’s method. Transmission characteristics of WR-90 waveguide have been well studied with the method of moment and mode matching technique in [1]. Thus, we will compare our results against theirs for L = 0 mm and L = 25 mm cases computing S11 parameter of T E10 mode by S11 = 20 log (|Γ|) = 20 log u S E10 · [˜ − E10 ]dS S E10 · E10 dS , (2.18) where u represents the total electric field computed numerically, and S denotes the surface ˜ of the input port. We investigate the wave propagation in the WR-90 waveguide using two different local approximation functions: plane wave basis and Legendre polynomials. The space of plane wave basis functions is vi (r) ∈ span e−jkl ·(r−ri ) , where ri is the PU domain center and kl is the direction of the plane wave basis functions given by (k, θm , φn ), θm is chosen as the roots of Legendre polynomials of order p, m = 1, 2, .., p, and φn = 2πn/(2p + 1), n = 0, 1, .., 2p, . Thus, the total number of basis functions per PU domain is Nb = p(2p + 1). 19 0 Measurement MoM VGFEM, Plane waves VGFEM, Legendre polynomials −10 −20 −30 −50 S 11 [dB] −40 −60 −70 −80 −90 −100 8 9 10 Freq [GHz ] 11 12 (a) L=0 mm 0 MoM VGFEM, Plane waves VGFEM, Legendre polynomials −10 −20 −40 S 11 [dB] −30 −50 −60 −70 −80 8 8.5 9 9.5 10 Freq [GHz ] 10.5 11 11.5 (b) L=25 mm Figure 2.5: S11 of WR-90 waveguide is simulated using different local approximation functions and compared against MoM data [1] for L = 0mm and L = 25mm. 20 Fig. 2.5 compares S11 obtained using VGFEM against MoM and the measurement data for a range of frequencies. Note, two different VGFEM simulations have been performed by using Legendre polynomials of order p = 2 for one simulation, Nb = 20, and 21 plane waves per PU domain for the other one, Nb = 21. VGFEM with plane waves gives better result than VGFEM with Legendre polynomials as shown Fig. 2.5(a). As the length of intermediate straight section gets longer, the results get much closer to MoM as seen in Fig. 2.5(b). Small deviation between VGFEM and MoM results can be explained with the geometry approximation of the curved sections and boundary conditions that are used to truncate the computational domain. Finally, we examine the contribution of higher order PU function on the accuracy. In Fig. 2.6, it is observed that the HO-PU function better captures the resonance. 0 −10 −20 S11 [dB] −30 −40 −50 MoM VGFEM, Hat−PU VGFEM, HO−PU −60 −70 −80 −90 −100 8 8.5 9 9.5 10 Freq [GHz ] 10.5 11 11.5 Figure 2.6: S11 of WR-90 waveguide with L = 25mm is simulated using different PU functions and compared against MoM data [1] . 21 2.2.2 Boundary Integrals Next, we apply VGFEM to the scattering problems from cavity backed apertures using the boundary integrals. Indeed, solution to such scattering problems has been attempted using mode matching, MoM/modal analysis, high frequency techniques, finite difference-time domain, etc.. However, the stand-alone application of each method has its own drawbacks such as restriction to rectangular cavities, truncation of infinite matrix, slow convergence, computation cost and truncation of the domain [2]. Alternatively, we introduce a hybrid VGFE-BI technique, which can yield an accurate and computationally efficient solution. The hybridization of VGFEM with BI is not as straightforward as that of FEM-BI. Difficulties that need to be overcome include the fact that Dirichlet boundary conditions are not easily imposed as the basis functions are not interpolatory. Also, unlike classical FEM-BI, an auxiliary space of surface basis functions that resides in the aperture can not be defined since these basis functions do not satisfy the Babuska-Brezzi condition. As a result, the formulation has to be done in terms of fields such that the vector basis spaces satisfy the integral boundary conditions. The resolution of these issues will be the subject of this Section. The method is validated by simulating radar crass sections (RCS) of small and large cavity-backed apertures and comparing them against published FEM-BI data and measurements. Consider the scattering from cavity-backed aperture in an infinite ground plane as illustrated in Fig. 2.7. The problem domain is divided into two regions: Region I (z > 0) and Region II (−c < z < 0). In Region I, the electric field satisfies the vector wave equation. Applying equivalence principle and image theory, the fields in both regions are decoupled, and the electric field in Region I is given by [13], 22 Figure 2.7: Illustration of scattering problem from a cavity-backed aperture. uI (r) = ui (r) + ur (r) − 2 Sa ¯ [ˆ × u(r′ )] · [∇′ × G0 (r, r′ )]dS, z (2.19) where ui (r) is the incident electric field radiated by the current source J in the free space, and ur (r) is the reflected electric field by the ground plane without the aperture with surface ¯ area Sa . The free space dyadic Green function G0 in (2.19) is 1 ¯ ¯ G0 = I + ∇∇ G0 (r, r′ ), 2 k0 (2.20) where k0 is the free space wavenumber, and G0 (r, r′ ) is the free space Green’s function, ¯ and I = (ˆx + yy + ˆˆ). The aperture magnetic field in Region I can be obtained by xˆ ˆ ˆ zz taking the curl of (2.19). The aperture fields in both regions are related by the continuity ˆ ˆ ˆ of the tangential electric and magnetic fields, z × uI (r) = z × uII (r) and z × [∇ × uI (r)] = ˆ z ×[∇×uII (r)] respectively. Thus, after taking the curl of (2.19) and applying the continuity equation at z = 0, the electric field in Region II is written as 23 2ˆ ˆ × ∇ × uII (r) − 2k0 z × z Sa ¯ ˆ ˆ z × uII (r) · G0 (r, r′) dS ′ = −2jk0 Z0 z × Hi (r) , (2.21) where Hi (r) is the incident magnetic field, and Z0 is the free space impedance. Equation (2.21) is one of the Neumann boundary conditions imposed at the aperture as the vector Helmholtz equation in Region II is solved using the VGFEM with Galerkin’s method. Additionally, although the normal components of the electric fields are discontinuous at the interface between two different media, they are continuous for this problem at the aperture ˆ due to domain homogeneity. Thus, imposition of the boundary condition ˆ ·uI (r) = z ·uII (r) z ensures smoother convergence of VGFEM-BI, and it can be readily shown as ˆ · uII (r) − 2ˆ· z z Sa ˆ (uII (r) × ˆ) × ∇G0 (r, r′ ) dS ′ = 2Ei (r) · z. z (2.22) As these boundary conditions are used, a bilinear form of VGFEM-BI can be written as V Spec 2 (∇ × w) · (∇ × u) dV − k0 Sa w · u dV + β ˆ ˆ [(w × n) · (∇ × u) + (u × n) · (∇ × w)] dS + β1 2 2k0 β2 V Sa ˆ ˆ (w · n)(u · n) dS − 2β2 ˆ (w · n) 2jk0 Z0 (∇ · w)(∇ · u)dV + (ˆ × w) · (ˆ × u)dS− n n 1 ˆ ¯ (u × n) I + ∇∇ G0 (r, r′ ) dS ′ 2 k0 Sa ˆ w×n· Sa Spec V Sa ˆ n· Sa dS+ ˆ (u × n) × ∇G0 (r, r′) dS ′ dS = ˆ (w × n) · Hi dS + 2β2 Sa ˆ ˆ (w · n)(Ei · n) dS, (2.23) 24 where the vector u(r) denotes the approximation to the electric field in Region II, w(r) represents the testing functions, Ei (r) is the incident electric field, V is the volume of the support domain of the testing function, Spec and Sa are the surfaces of the domain of testing ˆ function on the PEC and aperture respectively, and n is the outward unit normal vector of the support domain. The term with constant β in (2.23) is used to force divergence of the field to vanish. As stated in [25], the constants β, β1 and β2 contribute uniform convergence. Nitsche’s method is used to impose the Dirichlet boundary condition on the PEC walls. 2.2.2.1 Differences Between VGFEM-BI and FEM-BI At this stage, we should point out a key difference between implementation of VGFEM-BI and that of FEM-BI: In FEM-BI, an auxiliary trace space is defined on the boundary to ˆ represent the magnetic field that is related to the electric field via z × u = −M. However, it is well known that the definition of auxiliary space is not possible in VGFEM as the BabuskaBrezzi condition is not satisfied with this auxiliary space. As a result, the formulation has to be done in terms of fields [21]. However, unlike the approach in [21], the definition of patches in this paper helps to avoid the use of auxiliary Kirchoff surfaces. Since the PU domain is constructed from the geometry brick elements, there exists boundary patches that are defined around each boundary node. The vector basis functions defined on these boundary patches are forced to satisfy the boundary conditions. The integral with dyadic term in (2.23) can be rewritten using the vector identities and divergence theorem as 2 −2k0 Sa w×z ˆ 1 ˆ ¯ (u × z) I + ∇∇ G0 (r, r′ )dS ′ dS = 2 k0 Sa 25 2 − 2k0 −2 +2 +2 −2 (w × ˆ) · z Sa C C Sa Sa Sa C Sa ˆ G0 ( u × z) dS ′ dS ′ ˆ ˆ G0 ∇ · (u × ˆ) dS ′ (w × z) · nl dl z ˆ ˆ z ˆ G0 ( u × ˆ) · nl dl′ (w × z) · nl dl Sa C ′ ˆ ˆ G0 ∇ · (u × z) dS ′ ∇ · (w × z) dS ˆ G0 ( u × ˆ) · nl dl′ ∇ · (w × z) dS, z ˆ (2.24) where nl is the outward normal vector of the contour C. Contrary to the FEM-BI forˆ mulation, the line integrals along the boundary C are included in the formulation as the derivatives of the basis functions may not be continuous across the boundaries. Evaluation of the integrals in (2.24) is not straightforward as the PU function is piecewise continuous and the kernel is singular. The integrations are performed on each sub-domain using Gauss-Legendre rule for the prescribed accuracy. To elucidate evaluation of the integrals, consider the sub-domain shown in Fig. 2.8 and let us assume that it is on the aperture of an arbitrary cavity. For non-singular VGFEM terms, the integrals are evaluated on each sub-domain in the reference domain that is shown in Fig. 2.3(a). More specifically, for an integral Figure 2.8: Illustration of integration steps for singular integrals. 26 ¯ Ωj p p [ψi,j (r)vi (r)ˆ ] · [ψi,j (r)vi (r)ˆ ] dS, x x (2.25) ¯ where Ωj denotes the sub-domain of patch i, ψi,j (r) represents the portion of PU function ¯ ψi (r) in Ωj in the physical domain. This integral is performed numerically in the reference domain as N M n=1 m=1 p p wnm ψi,j (rnm )vi (rnm ) · ψi,j (rnm )vi (rnm ) |J|, (2.26) where N and M are the number of quadrature points, wnm are the quadrature weights corresponding to the quadrature points rnm , and |J| is the determinant of the Jacobian matrix. The transformation of the PU function ψi (η, ξ) that is defined in the reference domain, ψi (η, ξ) ⇔ ψi (r), is performed by using (2.13). Evaluation of VGFEM-BI terms in (2.24) is even more challenging due to basis function definitions of VGFEM. The surface-surface integral terms in (2.24) have singularity for the interaction of the patches that share this sub-domain. Specifically, consider an integral form, I= S p G0 (r, r′ )ψi (r′ )vi (r′ ) dS ′ , (2.27) where r and r′ are the source and observation points respectively. Evaluation of this integral is performed carefully in three steps as shown in Fig. 2.8. First, the quadratic sub-domain is divided into two triangles in physical domain, and then the rule presented in [30] is applied by dividing each triangle into three sub-triangles around the observation point r. The integration in each sub-triangle is done numerically using [30] −jkRl ′ )v p (r′ ) e , I ≈J Wl ψi (rl i l 4πRl l 27 (2.28) where J is the Jacobian of transformation and Wl are the weights corresponding to the area coordinate sample points. The reader is referred to [30] for details. In order to evaluate (2.28), we need to know the value of the PU function, ψi (r′ ), at r′ . However, we know the l l PU function is defined in the reference domain. Thus, we find the value of the PU function ¨ at the corresponding point in the reference domain by doing the transforms ψi (r) ⇔ ψi (u, v) ¨ and ψi (u, v) ⇔ ψi (η, ξ) respectively. The former is a simplex coordinate transform u = A1 /A and v = A2 /A, (2.29) where A1 and A2 are the areas of sub-triangles 1 and 2 respectively, A is the area of the original triangle as shown in Fig. 2.8, and (u, v) ∈ T3 = [0, 1] × [0, 1]. The latter transform is done by first scaling the domain T3 , then shifting it to the corresponding reference subdomain that is an element of the set T1 defined in (2.14). For the line-line integrals in (2.24), if an edge is shared by source and observation sub-domains, singularity is removed using Cauchy’s integral; otherwise, the integration is evaluated following the procedure for non-singular integral evaluation that has been mentioned before. 2.2.2.2 Numerical Results Next, we will investigate the capability of the hybrid VGFEM-BI by simulating scattering from cavity backed apertures. The accuracy of the method is validated by comparing the RCS results with FEM-BI data or measurement results. For this purpose, far zone scattered field is first calculated by using [13], Hs = ′ ′ jke−jkr ˆˆ ˆ ˆ z (θ θ + φφ) · (ˆ × E) ejk sin θ[x cos φ+y sin φ] dx′ dy ′, Zo 2πr Sa 28 (2.30) where (r, θ, φ) are the spherical coordinates of the observation point, (x′ , y ′ ) are the source points at the aperture surface Sa . Then, co-polarized and cross-polarized RCS of cavities are computed by using |H s |2 2 φ , σφφ = lim 4πr r→∞ |H i |2 φ and s2 2 |Hθ | . σθφ = lim 4πr r→∞ |H i |2 φ (2.31) Legendre polynomials and hat-PU function are used for VGFEM-BI simulations. First, VGFEM-BI is applied for computing RCS of a deep cavity with a small aperture. ˆ ˆ The dimensions of the cavity are 0.7λ×0.1λ×1.73λ. A θ polarized plane wave, Ei = θ e−jk·r with θinc = 40◦ is incident to the empty cavity as illustrated in Fig. 2.9(a). The vector basis functions of order p = 1 are used for VGFEM-BI. Fig. 2.9(a) compares the co-polarized and cross-polarized RCS of the cavity computed by VGFEM-BI and FEM-BI [2]. As it is seen, the agreement is excellent. Then, RCS of a cavity with larger aperture size, 1.5λ×1.5λ×0.6λ, is simulated. Brick elements with an approximate edge length of h = λ/5 is used. Fig. 2.9(b) shows the -p convergence of VGFEM-BI. The results with higher orders match with FEM-BI data [3] very well. Next, RCS of 2.5λ×0.25λ×0.25λ cavity is simulated for both E-polarized and H-polarized incident fields. The approximate edge length of brick elements and the order of polynomials used are h = λ/8 and p = 1 respectively. VGFEM-BI and FEM-BI results [4] for φinc = 90◦ agree very well as shown in Fig. 2.10(a). Fig. 2.10(b) shows the computed RCS for φinc = 0◦ using the first and second order vector basis functions. The RCS’s of the cavity due to both polarization for different orders of the basis functions agree well with each other. In the next simulation, we analyze a deeper cavity with the dimensions of a = 2λ0 , b = 2λ0 , c = 3λ0 . Brick elements with an approximate edge length of h = λ/8, and local 29 10 5 0 2 σ / λ [dB] −5 −10 −15 Co−pol, FEM−BI X−pol, FEM−BI Co−pol, VGFEM−BI X−pol, VGFEM−BI −20 −25 −30 −35 −40 0 10 20 30 40 50 φ [degrees] 60 70 80 90 (a) small aperture: a = 0.7λ, b = 0.1λ, and c = 1.73λ 20 VGFEM−BI, p=1 VGFEM−BI, p=2 VGFEM−BI, p=3 FEM−BI 10 2 σ / λ [dB] 15 5 0 −5 0 20 40 60 θ [degrees] 80 100 (b) large aperture: a = 1.5λ, b = 1.5λ, and c = 0.6λ Figure 2.9: Backscattered RCS of cavity-backed apertures compared against FEM-BI results given in [2] and [3] 30 20 10 0 2 σ / λ [dB] −10 −20 FEM−BI, H−pol FEM−BI, E−pol VGFEM−BI, E−pol VGFEM−BI, H−pol −30 −40 −50 −60 0 10 20 30 40 50 θ [degrees] 60 70 80 90 (a) φinc = 90◦ 20 VGFEM−BI, H−pol, p=1 VGFEM−BI, H−pol, p=2 VGFEM−BI, E−pol, p=1 VGFEM−BI, E−pol, p=2 10 0 σ / λ2 [dB] −10 −20 −30 −40 −50 −60 0 10 20 30 40 50 θ [degrees] 60 70 80 90 (b) φinc = 0◦ Figure 2.10: Backscattered RCS of the cavity-backed aperture with a = 2.5λ, b = 0.25λ, and c = 0.25λ compared against FEM-BI results [4] 31 30 25 20 10 2 σ / λ [dB] 15 5 VGFEM−BI FEM−BI 0 −5 −10 −15 0 10 20 30 40 50 θ [degrees] 60 70 80 90 Figure 2.11: Backscattered RCS of an empty cavity with a = 2λ0 , b = 2λ0 , c = 3λ0 forφinc = 0◦ . approximation of order p = 2 are used. Fig. 2.11 compares RCS of the cavity computed by VGFEM-BI and FEM-BI given in [31] , and results match. Finally, RCS of a long cavity with the dimensions of 16.26λ × 0.2λ × 0.85λ is simulated at 12 GHz using brick elements with edge length of approximately h = λ/4 and the order of approximation of p = 1. Fig. 2.12 shows the RCS results of VGFEM-BI against the measurement results provided in [5]. As seen in Fig. 2.12, VGFEM-BI results are very close to the measurement results like those of FEM-BI obtained in [5]. Thereafter, we analyze the scattering from cavities filled with materials. First, scattering from a filled cavity with the size of 2.89 in × 2.1 in × 0.057 in is analyzed at 9.2 GHz. The 32 0 −10 σ / m2 [dB] −20 −30 −40 VGFEM−BI, hp=λ/2, p=1 −50 Measurement −60 −70 0 10 20 30 40 50 θ [degrees] 60 70 80 90 Figure 2.12: Backscattered RCS of the cavity-backed aperture with a = 16.26λ, b = 0.2λ, and c = 0.85λ compared against measurement result at 12 GHz [5]. cavity is filled with a material having ǫr = 4. VGFEM-BI results match to the measurement results obtained in [5] as shown in Fig. 2.13. Next, RCS of a 1λ0 ×0.25λ0 ×0.25λ0 cavity filed with a material having ǫr = 7−j1.5 and µr = 1.8 − j0.1 is simulated. Brick elements with edge length of approximately h = λ0 /12 and local approximations order of p = 2 are used for VGFEM-BI. Figure 2.14 shows excellent agreement between VGFEM-BI and FEM-MoM results [32] for both incident angles. In the next example, we simulate the RCS of a cavity with a = 7.32cm, b = 5.2cm, c = 0.158cm over a range of frequency. Cavity is filled with a material ǫr = 2.17(1 − j0.001). The cavity region is meshed using brick elements with average edge length of h = λ0 /4. ˆ Local approximations order of p = 2 is used for VGFEM-BI. A plane wave E i = θe−jk.r , with θinc = 60◦ and φinc = 45◦ is incident. The simulation result match the FEM-BI data 33 0 Measurement VGFEM−BI, p=2 VGFEM−BI, p=3 −10 −30 2 σ/ m [dB] −20 −40 −50 −60 −70 0 10 20 30 40 50 θ [degrees] 60 70 80 90 Figure 2.13: Backscattered RCS of a 2.89 in × 2.1 in × 0.057 in cavity is analyzed at 9.2 ˆ GHz for E i = θe−jk.r , and compared against measurement [5] . The cavity is filled with ǫr = 4 . given in [5] as shown in Figure 2.15. Finally, we simulate the RCS of a partially filled cavity depicted in 2.16(a). Cavity dimensions are a = 0.3λ0 , b = 0.1λ0 , c = 0.6λ0 , d = 0.2λ0 . Cavity is filled with a material of ǫr = 2 − j2. The cavity region is meshed using brick elements with average edge length of h = λ0 /10. Local approximations order of p = 1 is used for VGFEM-BI. A plane wave ˆ E i = θe−jk.r , with θinc = 40◦ is incident. The simulation results match the FEM-BI data given in [33] as shown in Figure 2.16(b). 34 10 0 σ / λ2 [dB] −10 −20 −30 VGFEM−BI, Epol VGFEM−BI, Hpol FEM−MoM, Epol FEM−MoM, Hpol −40 −50 −60 0 10 20 30 40 50 θ [degrees] 60 70 80 90 (a) φinc = 90◦ 10 VGFEM−BI, Epol VGFEM−BI, Hpol FEM−MoM, Epol FEM−MoM, Hpol 0 σ / λ2 [dB] −10 −20 −30 −40 −50 −60 0 10 20 30 40 50 θ [degrees] 60 70 80 90 (b) φinc = 0 Figure 2.14: Backscattered RCS of 1λ0 × 0.25λ0 × 0.25λ0 cavity filled with ǫr = 7 − j1.5 and µr = 1.8 − j0.1 35 −40 −45 σθθ / m2 [dB] −50 −55 −60 −65 −70 2 3 4 5 6 Frequency [GHz] 7 8 Figure 2.15: Backscattered RCS of a cavity with a = 7.32cm, b = 5.2cm, c = 0.158cm over a range of frequency. ǫr = 2.17(1 − j0.001) 2.2.3 Absorbing Boundary Conditions In this subsection, we apply VGFEM to solve some scattering problems of perfect electrically conducting (PEC) scatterers in free space. To truncate open domain, we use absorbing boundary conditions (ABC). Next, we write the bilinear formulation of the VGFEM along with the ABCs. We first start the formulation using the fist order absorbing boundary condition [13] ˆ × ∇ × u(r) + jk0ˆ × ˆ × u(r) = 0, r r r (2.32) where ˆ is the radial unit vector. As scattered field formulation is used, the boundary r condition on PEC scatterer is given by n׈ ×u(r) = −ˆ ׈ ×uinc (r), where n is the outward ˆ n n n ˆ 36 (a) Cavity geometry −10 −15 −20 2 σ [dB / λ ] −25 −30 −35 −40 VGFEM−BI, σEE −45 VGFEM−BI, σHH VGFEM−BI, σEH −50 FEM/MoM −55 −60 0 10 20 30 40 50 φ [Degrees] 60 70 80 90 (b) RCS Figure 2.16: Backscattered RCS of a partially filled cavity with a = 0.3λ0 , b = 0.1λ0 , c = 0.6λ0 , d = 0.2λ0 , and ǫr = 2 − j2. 37 normal vector of the scatterer, u(r) and uinc (r) are the scattered and incident electric fields respectively. With these boundary conditions, a bilinear formulation of VGFEM can then be obtained as Ωi β β1 Γpec Ωi 2 (∇ × w) · (∇ × u) dΩ − k0 (∇ · w)(∇ · u)dΩ + Γpec (w × n) · (u × n)dΓ + jk0 ˆ ˆ = −β1 w · u dΩ+ (w × n) · (∇ × u)dΓ+ ˆ Γabc Γpec Ωi (w × n) · (u × n)dΓ ˆ ˆ (w × n) · (uinc × n)dΓ, ˆ ˆ (2.33) where the constants β, β1 and β2 contribute uniform convergence [25]. Note, we assumed ˆ ∼ n in this formulation [13]. r ˆ Next, the VGFEM with ABC is validated via a number of simulations. First, bistatic RCS of a 0.3m × 0.3m PEC plate is simulated using the first order absorbing boundary condition. The plate is centred at the origin on xy plane, and it is illuminated by the electric field E = xejk0 z at 300MHz. A cubical truncation boundary is placed 0.3m away from ˆ edge of the plate. The problem domain is meshed using 2493 nodes and 12312 tetrahedrons with the average edge length of h = 0.1m. Figure 2.17 shows the RCS plots for p = 1 and p = 2 of local approximations. The RCS result with p = 2 matches the MoM result obtained using our in-house code. Next, bistatic RCS of a sphere with radius of r = 0.5λ, where λ = 1.0m, is simulated. The sphere is centred at the origin, and it is illuminated by the electric field E = xejk0 z . ˆ A spherical truncation boundary is placed 0.5λ away from the sphere, and the first order ABC is imposed on the boundary. The problem domain is meshed using 531 nodes and 38 0 −5 −10 σ [dBsm] −15 MoM VGFEM, p=1 VGFEM, p=2 −20 −25 −30 −35 −40 0 10 20 30 40 50 θ [degrees] 60 70 80 90 Figure 2.17: Scattering from a patch 1762 tetrahedrons with the average edge length of h = λ/5. While p = 1 is used in the interior PU domains associated with the interior nodes, p = 2 is used in the boundary PU domains associated with the boundary nodes. Figure 2.18 shows the RCS results obtained using VGFEM, the commercial software HFSS, and Mie series. Note that HFSS uses 2nd order ABC and, simulation parameters used are h = λ/10 and p = 1. VGFEM with varying p gives accurate enough results despite the coarse tetrahedral mesh. 2.2.4 Perfectly Matched Layers In this dissertation, we use the anisotropic interpretation of the PML; therefore, we reformulate the VGFEM for anisotropic medium. We should remember that the vector basis space in (2.3) is valid only for the continuous field representations. It is apparent that use 39 20 Mie Series HFSS, havg=λ/10, p=1 15 VGFEM, h =λ/5, p=1 and p=2 avg 10 σ [dBsm] 5 0 −5 −10 −15 −20 0 30 60 90 120 θ [degrees] 150 180 Figure 2.18: Scattering from a sphere with radius of r = 0.5λ. Spherical ABC surface is placed at rabc = 1.0λ of PML for domain truncations requires the termination of the problem with complex lossy materials, which means that it is necessary for the basis function spaces used be normally discontinuous and tangentially continuous across material boundaries. In addition, it is required that the normal derivatives of the tangential components be discontinuous as well. Indeed, such basis space has been developed in [25] by defining an auxiliary set of functions that satisfied all these conditions. In [25], it is assumed that the interface is planar and occupied either the entire width or height of the PU domain. However, if this method were to be applied to more practical problems, we should develop other mechanisms that can easily handle the discontinuities. Thus, we develop methodologies to accommodate these inhomogeneities within the VGFEM framework. Consider a linear, anisotropic, and source free region Ω, whose boundary is denoted by 40 ∂Ω := Γ = i Γi . In Ω, the electric field, u(r), satisfies the vector Helmholtz equation 2¯ ∇ × µr −1 · ∇ × u(r) − k0 ǫr · u(r) = 0 ¯ Bi {u(r)} = gi (r), r ∀ ∈ Γi , (2.34) where ǫr and µr are the permittivity and permeability tensors, k0 is the free space wavenum¯ ¯ ber, Bi is a differential operator, and gi (r) is a boundary condition imposed on Γi . To solve the vector wave equation, we need to define appropriate basis functions that satisfy the boundary conditions at material interfaces. Since the vector space defined in the previous Section can be only used to represent continuous fields, we develop two techniques to include the discontinuities in the solution space. The first technique is to augment this continuous basis functions space with additional discontinuous basis functions. The second is to decompose the domain into homogeneous domains, and then couple them via some transmission conditions. Next, we will explain each method. 2.2.4.1 Additional Basis Functions For Discontinuities In Section 2.1.2, we constructed the partition of unity domains as a union of polyhedrons that share a node; i.e., Ωi = L Ω . If any two of the constituents of the domain Ω , ¯ i j j ¯ viz., Ωj have constitutive parameters that differ from each other, then we have to define additional basis functions to account for discontinuities in the normal component of the fields (and the normal derivatives of the tangential components. However, note that one ¯ needs to define basis that accomplishes these goals on the boundary Γj between Ωj and ¯ Ωj+1 that is equipped with a normal n. To define these basis functions, we take advantage ˆ of a salient feature of the partition of unity, viz., one can use any function as long as the 41 local partition of unity goes smoothly to zero at the boundary of the domain. With this, the additional basis function can be defined as p ˜ V = φj1 (r) ∇ × ∇ × (ˆ vj,n (r)) − vn + φj2 (r)vn , n (2.35) vn ˜ vt ˜ where vt , and vn respectively denote the tangential and normal component of the basis ˜ ˜ p functions, n is the vector normal to the junction face, vn = n[ˆ · ∇ × ∇ × (ˆ vj,n (r))], and ˆ ˜ ˆn n the functions φj1 (r) and φj2 (r) are used for forcing v to satisfy the boundary conditions at ˜ the material junction, 0 ˜ ∂n vt |∀r∈Γ j 0 ˜ ∂n vn |∀r∈Γ j =0 k ˜ ∂n vt |∀r∈Γ j =0 =0 k ˜ ∂n vn |∀r∈Γ j = 0, (2.36) where k = 1, 2, .....p − 1. It should be noted that these basis functions are imposed for ¯ r ∈ Ωj and are associated with the boundary Γj . It follows, then that these basis functions should vanish at other faces. For the tangential components, it is easily achieved by defining φj1 (r) as a tensor product of one dimensional functions. Here, we choose raised-cosine filter function,    Tr  if |x| ≤ α     πTr Tr f (x) =  2 [1 + cos( βr [|x| − α])] if α ≤ |x| ≤ β       0 elsewhere, (2.37) 1−β 1+β where α = 2T r and β = 2T r , βr is the roll-off factor, and Tr is the period of the r r impulse response of this filter. Figure 2.19(a) illustrates the function φj1 (r) = f (x)f (y)f (z) at y = 0, which vanishes at the boundaries. A special care is needed for φj2 (r) as it should not vanish on Γj except at the boundaries and this can be done by setting filter parameters 42 φ1 1 0.5 0 1 0.5 0 1 0.5 −0.5 z 0 −1 −0.5 −1 x (a) φ1 φ2 0.4 0.2 0 1 0.5 0 1 0.5 −0.5 0 z −1 −0.5 −1 x (b) φ2 Figure 2.19: Shape functions used for the definition of additional basis functions. 43 such that it abruptly vanishes at the boundaries. Figure 2.19(b) illustrates the function ˜ φj2 (r) = f (x)f (y)f (z) at y = 0 for βr = 0.01 and Tr = 0.505, where Γj is the surface ˜ z = 1, and f (z) = ((z − 1)/2 + 1) In addition, an example of the vector basis functions are demonstrated in Figure 2.20. A material discontinuity at z = 0.5 of the sub-domain (0, 0.5) × (0, 0.1) × (0, 0.5) is considered and the vector functions are plotted on y = 0.05 plane. Note that shape functions are defined on the isoparametric cell and the proposed vector basis function can be used on any linear shaped domains. Next, we formulate the bilinear form of VGFEM using the boundary conditions ˆ u(r) × n = gd (r) , r ∈ Γd ˆ n × (µr −1 · ∇ × u(r)) = gn1 (r), ¯ r ∈ Γn ˆ ǫ n · (¯r · u(r)) = gn2 (r), r ∈ Γn , (2.38) ˆ where n is the outward unit vector normal to the boundary Γ, Γ = Γd ∪ Γn . Then, the bilinear form of the VGFEM for anisotropic medium is written as Ω 2 (∇ × w) · (µr −1 · ∇ × u) dΩ − k0 ¯ β Γd β1 Γd Ω Ω w · (¯r · u) dΩ+ ǫ [∇ · (¯r · w)] [∇ · (¯r · u)] dΩ+ ǫ ǫ ˆ ˆ (w × n) · (µr −1 · ∇ × u) + (u × n) · (µr −1 · ∇ × w) dΓ+ ¯ ¯ ˆ ˆ (w × n) · (u × n) dΓ + β2 Ω Γn [ˆ · (¯r · w)][ˆ · (¯r · u)] dΓ = n ǫ n ǫ f · w dΩ + gd · (µr −1 · ∇ × w) dΓ+ ¯ 44 Γd 0.6 0.5 0.4 z 0.3 0.2 0.1 0 −0.1 −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.5 0.6 x (a) Tangential component 0.6 0.5 0.4 z 0.3 0.2 0.1 0 −0.1 −0.1 0 0.1 0.2 0.3 0.4 x (b) Normal component Figure 2.20: A vector plot of the additional basis functions 45 β1 Γd ˆ (w × n) · gd dΓ + Γn w · gn1 dΓ + β2 Γn [ˆ · (¯r · w)]gn2 dΓ. n ǫ (2.39) In Eq. (2.39), the divergence term is rewritten as 3 ∇ · (¯r · (ϕA)) = ∇ϕ · (¯r · A) + ϕ ǫ ǫ 3 ǫij ∂i Aj , (2.40) i=1 j=1 where ϕA is a general representation of the vector basis functions, ϕ denotes either the PU function ψi (r) or the functions ψj1 (r) and ψj2 (r), the vector A represents the vector parts in those spaces, and ǫrij denotes the entries of the matrix ǫr . ¯ 2.2.4.2 Domain Decomposition for Discontinuities In this approach, we design the PU domains carefully. If there is material variation in a PU domain, then we decompose the PU domain according to the number of materials the PU domain. This separation and redefinition of PU domains in new domains are achieved by duplicating the nodes at the material interfaces. An illustration of node duplication is shown in Fig. 2.21. Once the PU domains are separated, transmission conditions are applied to couple them. Next, we formulate the transmission conditions considering a PU domain with Figure 2.21: An illustration of material interfaces and domain decomposition 46 two different materials (Region I and II ) depicted in Fig. 2.21. Let u1 and u2 denote the electric fields in Region I and II respectively. Then, transmission conditions at Γ12 and Γ21 are respectively given by βtc1 n1 × u1 × n1 + βtc2 n1 (ˆ 1 · ǫr1 u1 ) = βtc1 n2 × u2 × n2 + βtc2 n2 (ˆ 2 · ǫr2 u2 ) ˆ ˆ ˆ n ˆ ˆ ˆ n (2.41) and βtc1 n2 × u2 × n2 + βtc2 n2 (ˆ 2 · ǫr2 u2 ) = βtc1 n1 × u1 × n1 + βtc2 n1 (ˆ 1 · ǫr1 u1 ), ˆ ˆ ˆ n ˆ ˆ ˆ n (2.42) where βtc1 and βtc2 are the penalty factors. As opposed to this conditions, the Robin transmission condition 1 1 ∇ × u1 ) + βrb n1 × u1 × n1 = −ˆ 2 × ( ˆ ˆ n ∇ × u2 ) + βrb n2 × u2 × n2 ˆ ˆ n1 × ( ˆ µr1 µr2 (2.43) and 1 1 n2 × ( ˆ ∇ × u2 ) + βrb n2 × u2 × n2 = −ˆ 1 × ( ˆ ˆ n ∇ × u1 ) + βrb n1 × u1 × n1 ˆ ˆ µr2 µr1 (2.44) √ can also be used, where βrb = jk0 ǫr µr . Then, one can obtain a bilinear formulation using ¯ ¯ domain decomposition for Region I as 47 Ω1 (∇ × w1 ) · ( 1 2 ∇ × u1 ) dΩ − k0 ǫr1 w1 · u1 dΩ + w1 · B1 {u1 } µr1 Ω1 ∂Ω1 \Γ12 βtc1 βtc1 Γ12 w1 × n1 · u1 × n1 dΓ + βtc2 ˆ ˆ Γ12 Γ12 w1 × n2 · u2 × n2 dΓ + βtc2 ˆ ˆ (w1 · n1 )(ǫr1 u1 ) dΓ = ˆ Γ12 (w1 · n2 )(ǫr2 u2 ) dΓ ˆ (2.45) and for Region II as Ω2 (∇ × w2 ) · ( 1 2 ∇ × u2 ) dΩ − k0 ǫr2 w2 · u2 dΩ + w2 · B2 {u2 } µr2 Ω2 ∂Ω2 \Γ21 βtc1 βtc1 Γ21 w2 × n2 · u2 × n2 dΓ + βtc2 ˆ ˆ Γ21 w2 × n1 · u1 × n1 dΓ + βtc2 ˆ ˆ Γ21 (w2 · n2 )(ǫr2 u2 ) dΓ = ˆ Γ21 (w2 · n1 )(ǫr1 u1 ) dΓ, ˆ (2.46) where n1 and n2 are outward normal vectors of Γ12 and Γ21 respectively. ˆ ˆ 2.2.4.3 Numerical Results In this Section, we first validate the new bilinear formulation of the VGFEM for anisotropic medium by simulating cavity modes. Then, we validate the techniques proposed for discontinuity via some wave propagation in inhomogeneous domains. Finally, we investigate the VGFEM with PML by simulating wave propagation and scattering in rectangular waveguides. In all simulations, hat-PU function and Legendre polynomials as local approximation function are used. 48 A 1.0 m × 0.4454 m × 0.3 m rectangular cavity filled with a uniaxial anisotropic material is used to validate the VGFEM formulation. The material parameters are ǫr11 = ǫr22 = 2.1, ǫr33 = 4.6, µr11 = µr33 = 3.2, and µr22 = 4.0. The average size of the brick elements used for meshing is h = 0.15 m. First eight eigenvalues computed by VGFEM are tabulated in Table 2.1 for p = 1 and p = 2. The eigenvalues computed by commercial software are also shown in Table 2.1 for comparison. Note, HFSS uses adaptive solution with higher order basis and the final average element size is h = 0.062 m. Although the coarse meshing is used for VGFEM, it captures the eigenvalues accurately. Table 2.1: Eigenvalues of a 1.0 m × 0.4454 m × 0.3 m rectangular cavity filled with the uniaxial anisotropic material Mode V GF EM, p = 1 V GF EM, p = 2 HF SS CST MW S 1 1.980 1.979 1.978 1.976 2 2.353 2.351 2.349 2.347 3 2.870 2.865 2.863 2.859 4 3.491 3.460 3.456 3.447 5 3.789 3.750 3.746 3.734 6 3.999 3.959 3.955 3.943 7 4.187 4.102 4.094 4.076 8 4.232 4.218 4.214 4.210 Next, we validate the additional basis functions and domain decompositions to approximate the discontinuous fields. A plane with the propagation direction of θinc = π/6 and φinc = 0 is considered for this study . Computation domain of 1m × 0.1m × 1m is filled with ǫr = 4 for z > 0.5m and ǫr = 4 for z < 0.5m. The domain is meshed using brick elements with an average length of h = λ0 /6, where λ0 = 1m. βtc = 103 is for domain decomposition. Figure 2.22 shows the normal component of the computed electric field. Both techniques are captured the field with a relative L2 error of ≈ 0.048. For the next example, a slab with a thickness of = 0.25m and dielectric constant of ǫr = 4 is considered along z axis . The slab is centred in 1m × 0.1m × 1m domain, i.e. at z = 0.5m. Now, the computation domain is partitioned into 3 sub-domains as domain decomposition 49 1 Real {Ez} 0.5 0 −0.5 −1 1 0.8 0.6 0.4 0.2 x (m) 0 0.8 1 0.6 0.4 0.2 0 z (m) (a) Real part of Ez 1 Im {Ez} 0.5 0 −0.5 −1 1 0.5 0 0.2 0.4 0.6 x (m) 0 0.8 1 z (m) (b) Imaginary part of Ez Figure 2.22: The normal component of the computed electric field propagating in a 1m × 0.1m × 1m domain filled with ǫr = 4 for z > 0.5m and ǫr = 4 for z < 0.5m 50 technique is used. The domain is meshed using brick elements with an average length of h = λ0 /6, where λ0 = 1m. βtc = 103 is for domain decomposition. Figure 2.23 shows the normal component of the computed electric field. Both techniques are captured the field with a relative L2 error of ≈ 0.017. Next, the use of one-dimensional PML with VGFEM is investigated for T E10 mode propagation in 1.0 λ0 × 0.5λ0 × 1.0λ0 waveguide. A PEC backed PML layer with the thickness of 0.25λ0 is attached to the end of the waveguide. A uniaxial anisotropic material with µr11 = µr22 = ǫr11 = ǫr22 = s′ − js′′ and µr33 = ǫr33 = 1/(s′ − js′′ ) is used in the PML region, where s′ and s′′ are real variables. For the simulations, Legendre polynomials of order p = 1 is used in all regions; the average element length in the waveguide region is h = 0.25λ0 ; and, the the average element lengths used in the PML region are h = {0.25λ0 , 0.125λ0 , 0.083λ0 } for Nlayer = {1, 2, 3} respectively (Nlayer denotes the number of elements along z direction in the PML region). Fig. 2.24 shows the plots of the reflection coefficient at the PML interface, z = 1.0λ0 , versus s′′ for different layers. As the loss factor s′′ is varied keeping s′ constant, s′ = 1, the reflection is monotonically decreased until s′′ = 7 for Nlayer = 1, and then it starts to increase due to the insufficient field representation in PML region as shown in Fig. 2.24(a). However, as the number elements in PML region is increased, the rate of increase is significantly reduced. In the case of s′ = s′′ , the behaviour is similar to the case s′ = 1 as seen in Fig. 2.24(b); however, the error in PML region is now higher. Finally, we show an application of VGFEM with PML for inhomogeneous problems. Scattering from a dielectric object in a rectangular waveguide depicted in Figure 2.25(a), a = 2b, d = 0.8b, c = 3d, and b = 1mm, is considered for this demonstration. The problem 51 (a) Real part of Ez (b) Imaginary part of Ez Figure 2.23: The normal component of the computed electric field propagating in a 1m × 0.1m × 1m domain through a slab with a thickness of = 0.25m and dielectric constant of ǫr = 4 52 0 N =1 N =2 layer layer −10 Nlayer=3 Γ [dB] −20 −30 −40 −50 −60 0 5 10 s´´ 15 20 15 20 (a) s′ = 1 0 Nlayer=1 Nlayer=2 −10 Nlayer=3 Γ [dB] −20 −30 −40 −50 −60 0 5 10 s´´ ′′ (b) s′ = s Figure 2.24: Reflection coefficient at the PML interface for the wave propagation in a 1.0 λ0 × 0.5λ0 × 1.0λ0 rectangular waveguide 53 domain is partitioned into two homogeneous sub-domains: free space and dielectric regions. Entire domain is meshed using brick elements with an average edge length of h = 2.2mm. Relative dielectric constant of ǫr = 6 and b = 10mm are used. The local approximation order of p = 2 is used. Figure 2.25(b) shows computed reflection coefficients with and without PML. The VGFEM result without PML agrees with the reference data given in [13]. PEC backed PML region with thickness of d/2, and material parameters of s′ = 5 and s′′ = 0 is attached to the end. This PML is not good at high frequency regime as shown in the figure. 2.3 Dispersion Analysis In this Section, dispersion characteristics of GFEM are studied. The motivation herein is that the numerical solution to the Helmholtz equation with a high wave numbers can notably deviate from the exact solution despite a moderate resolution of the geometry. This is largely due to phase error that is endemic in differential equation solvers. Since the phase error is progressive throughout the domain, it especially plays a key role on the accurate solution of electrically large and complex geometries. While dispersion analysis has been performed for classical FEM [34–36], it is a little more difficult for GFEM due to overlapping PU domains. To this end, we develop a semi-analytic technique to analyze phase error in GFEM. This technique presented herein is applicable to the dispersion analysis of both scalar and vector problems. Consider an infinite, linear, homogeneous, isotropic, and source free region. In order to better understand dispersion analysis within the framework of GFEM, we first perform analysis for the wave propagation in 1D. This will then be followed by dispersion analysis 54 (a) Geometry 1 0.9 0.8 VGFEM p=2, h~λ/5 VGFEM p=2, h~λ/5, with PEC−PML Orthogonal expansion HFSS 0.7 Γ 0.6 0.5 0.4 0.3 0.2 0.1 0 8 9 10 11 Frequency [GHz] 12 13 (b) Reflection Coefficient Figure 2.25: Reflection coefficient is computed using domain decomposition. 55 Figure 2.26: Definition of 1D PU function on PU domains separated by nodal distance h for α=1.5. in two-dimensional (2D) scalar and vector problems. An exact solution to the 1D scalar Helmholtz equation, ∇2 u(x) + k 2 u(x) = 0, (2.47) in this region is a plane wave propagating along the x-axis, u(x) = e−jkx , where k is the wavenumber, and unity amplitude is assumed. To solve the problem using GFEM, the computation domain is first covered by the PU domains Ωi surrounding nodes i as shown p in Fig. 2.26. Then, the scalar basis functions are defined by ψi (x)vi (x) on the domains with the size of hp , which is related to the nodal distance h by a constant α, α > 1 and hp = αh. The PU function ψi (x) is constructed using either of the weighting functions given p in (2.6) and (2.7), and the local approximation functions vi (x) are chosen from Legendre polynomials. Assume that there are Nb number of basis functions in a PU domain Ωi ; the same set of basis functions is assumed to exist in all PU domains that extend to infinity to the left and right of this parent PU domain. Since the PU domain Ωi overlaps only with the adjacent PU domains Ωi−1 and Ωi+1 , a sparse and diagonal matrix is constructed applying Galerkin’s method. The system consists of square sub-block matrices [Di−1 ], [Di ] and [Di+1 ] related 56 to the intersection domains, and thus elements are given by Di±1,qp = Ωi ∩Ωi±1 p q −∇ui±1 (x) · ∇wi (x)+ q p k 2 ui±1 (x)wi (x) dx, Di,qp = Ωi p p q q −∇ui (x) · ∇wi (x) + k 2 ui (x)wi (x) dx, (2.48) (2.49) p q where p and q are the indexes of the basis and testing functions (ui , wi ) respectively. Since the system is infinitely periodic, the coefficients are the phase shifted version of each other. N Thus, let [ai ] = [a1 , a2 , ..., ai b ] be the block vector representing the unknown coefficients i i of the basis functions in Ωi , and kh be the numerical wavenumber. Then, the coefficients of the basis functions in the neighbouring PU domains are [ai−1 ] = [ai ]ejkh h and [ai+1 ] = [ai ]e−jkh h . (2.50) Setting the weighted residue integral to 0 for the PU domain Ωi , a block matrix equation is obtained, [Di−1 ]ejkh h + [Di ] + [Di+1 ]e−jkh h [ai ] = 0. (2.51) A A non-trivial solution to (2.51) exists if |A| = 0, which has 2Nb roots. However, there is only one valid root that satisfies the conditions 0 < kh < 2k and kh ∈ R. Then, the phase error per wavelength in degrees is computed using the conventional definition, δp = 360 | k − kh | . k 57 (2.52) Figure 2.27: Overlapping PU domains separated by a nodal distance h in 2D. We should note that the phase error is a function of h and p. Since the dispersion analysis of GFEM is considerably more complicated unlike FEM, it is difficult to get an analytic dispersion expression as a function of k and h by solving the determinant. Generalization of the technique to two-dimension is straightforward but results in more complicated analytical expressions; therefore, we briefly describe the extension of the technique here. Assume a z polarized plane wave with an incident angle θ from the x -axis ˆ propagates over an infinite 2D domain depicted in Fig. 2.27. The PU domain Ω ′ , where ii ′ additional index i is used to denote the discretization along y-axis, now overlaps eight adjacent PU domains. The block coefficient vectors of the adjacent domains are the shifted form of the coefficient vector [a ′ ] of the domain Ω ′ . The amounts of shifts along the x -axis and y-axis are reii ii ±jkh1 h and e±jkh2 h with the corresponding k and k components spectively given by e h1 h2 58 of the numerical wavenumber kh . Then, the block matrix equation for Ω ′ is ii  jk h+jkh2 h + [D jk h [D ′ ]e h1 ′ ]e h2 + i−1,i −1 i,i −1 −jkh1 h+jkh2 h + [D jk h [D ′ ]e ′ ]e h1 + i+1,i −1 i−1,i          [D ′ ] +  ii   −jkh1 h + [D jk h−jkh2 h +  [D ′ ]e ′ ]e h1  i+1,i i−1,i +1   −jkh1 h−jkh2 h [D ′ ]e−jkh2 h + [D ′ ]e i,i +1 i+1,i +1 A          [a ′ ] = 0.  ii       (2.53) Assume that numerical wave propagates along the incident wave. Then, (2.53) is reduced to a one variable block matrix equation by using kh1 = kh cos θ and kh2 = kh sin θ. As in 1D, kh is determined by finding the appropriate root of the equation |A| = 0 for each h and p values. Incident angle-phase error relation is also determined by extracting kh for each incident angle. This technique can also be applied to vector problems. Again, consider a 2D domain that is identical to the one considered earlier. An exact solution for the vector Helmholtz equation in this domain is a plane wave u = [− sin θˆ + cos θˆ]e−jk[cos θx+sin θy] , x y (2.54) where θ is the incident angle, and unity amplitude is assumed. For the solution of the problem with VGFEM, the PU domains are defined as in the scalar case as shown in Fig. 2.27 and p p ˆ vector basis functions are defined by ui (r) = ψi (r)∇ × z vi (r) . Legendre polynomials are used for local approximation function. Following the same procedure in scalar dispersion 59 analysis, the block matrix equation identical to (2.53) is established. However, the elements of the block matrices for non-self and self-PU domain interactions are now respectively given by D qp p q p q = ∇ × u (r) · ∇ × w ′ (r) − k 2 u (r) · w ′ (r) dS, l l l ii ii Ω ′ ∩Ωl ii (2.55) and p q p qp q ∇ × u ′ (r) · ∇ × w ′ (r) − k 2 u ′ (r) · w ′ (r) dS, D ′= ii ii ii ii ii Ωii′ ∩Ωii′ (2.56) where p and q are the indexes of vector basis and testing functions respectively, and l denotes any combination of the indexes of the overlapping domains given in (2.53) providing that l = ii′ . Numerical wavenumber kh is determined following the same procedure in the scalar analysis. Using the valid root, h−, and p−convergence of VGFEM with Legendre polynomials can be analyzed for various incident angles. The proposed semi-analytic technique can be easily extended to 3D scalar and vector dispersion analysis of GFEM. 2.3.1 Numerical Results In this subsection, h−, and p−convergence of the dispersion and its dependence on incidence angle are presented for both scalar GFEM and VGFEM via a set of simulations. First dispersion results for scalar GFEM and then the results for VGFEM are presented. Legendre polynomials and hat-PU function that is given by (2.6) are used unless otherwise is stated. The parameter α relating element size h to PU domain size hp , hp = αh, is set to 1.5 for all simulations. Fig. 2.28(a) shows the 1D dispersion results for scalar GFEM. For the first 60 2 10 0 Phase error [degrees] 10 −2 10 −4 10 −6 10 −8 10 −10 10 5 p=1, FEM p=2, FEM p=3, FEM p=1, GFEM p=2, GFEM p=3, GFEM 7 10 15 20 λ/h 40 80 40 80 (a) FEM vs. GFEM, 1D 0 Phase error [degrees] 10 −5 10 p=1, θi=0 p=2, θi=0 p=3, θ =0 i p=1, θ =π/4 i −10 10 p=2, θi=π/4 p=3, θi=π/4 5 7 10 15 20 λ/h (b) h−, and p−convergence for axial and diagonal incidences, 2D Figure 2.28: h−, and p−convergence of the phase error per wavelength in degrees for scalar GFEM using Legendre polynomials and hat-PU functions. 61 90 8 120 60 6 150 p=1 p=2 p=3 30 4 2 180 0 −4 ×10 210 330 −2 ×10 300 240 270 Figure 2.29: Incident angle dependency of the phase error per wavelength in degrees for scalar GFEM with hat-PU function, h = λ/10 . order polynomials, both scalar FEM and GFEM show approximately the same convergence behaviour, whereas GFEM exhibits considerably better performance for the higher orders. It is well-known that the order of phase error in FEM is O([h/λ]2p ) [36]; therefore, we conclude that the phase error in GFEM with Legendre polynomials is of the same order with an additional p dependent constant coefficient, δp , which shifts the GFEM results down as compared to FEM results. Next, 2D dispersion results of scalar GFEM are shown in Fig. 2.28(b) for axially and diagonally incident wave. The numbers of scalar basis functions (per PU domain) used for the simulations are 3, 6 and 10 for the approximation orders of 1, 2 and 3, respectively. The axial results are identical to those in 1D, as expected. The phase error of diagonally incident wave becomes larger than the phase error of axially incidence wave as the polynomial 62 order increases, and this relation does not change as the mesh size increases as seen in Fig. 2.28(b). Next, incident angle dependency of the phase error is further investigated for various polynomial orders in Fig. 2.29. Note that the results for p = 2 and p = 3 are scaled by the factors of 100 and 1000, respectively, to plot them on the same figure. For p = 1 and h = λ/10 , the phase error for the axially incident wave is about 2◦ more than that for the diagonally incident wave as in FEM with quadrilateral mesh structure [36]. In FEM, the phase error for the diagonally incident wave is always smaller and the shape of the plot is preserved as the polynomial order increases (it is only scaled) [36]. On the other hand, the phase error of GFEM for the diagonally incident wave gets larger as the polynomial order increases and each polynomial order exhibits different incident angle dependency plot as seen in Fig. 2.29. Next, we investigate the dispersion characteristics of scalar GFEM using higher-order (HO) PU function that is obtained using (2.7). Fig. 2.30(a) shows that the h−convergence plots of both HO-PU and hat-PU functions are identical for p = 1, whereas HO-PU exhibits better phase error performance for higher orders of Legendre polynomials. In addition, we investigate the incident angle dependency of the phase error for HO-PU function. As HO-PU results in Fig. 2.30(b) are compared with hat-PU results in Fig. 2.29, it is observed that the phase error for the diagonally incident wave gets larger as the polynomial order increases for both PU functions. However, HO-PU function suppresses the phase error for axially incident wave much more than hat-PU function. Thereafter, the dispersion results for VGFEM is presented and they are compared with those of scalar GFEM. The numbers of vector basis functions per PU domain used for the simulations are 5, 9 and 14 for the polynomials orders of 1, 2 and 3, respectively. h−, and 63 2 10 0 Phase error [degrees] 10 −2 10 −4 10 −6 10 −8 10 −10 10 5 p=1, hat−PU p=2, hat−PU p=3, hat−PU p=1, HO−PU p=2, HO−PU p=3, HO−PU 7 10 15 20 λ/h 40 80 (a) Hat-PU vs. HO-PU, axial incidence (b) Incidence angle dependency, h = λ/10 Figure 2.30: h−, and p−convergence and incident angle dependency of the phase error per wavelength in degrees for scalar GFEM with HO-PU function. 64 2 10 0 Phase error [degrees] 10 −2 10 −4 10 p=1, θ =0 i −6 p=2, θ =0 10 i p=3, θi=0 −8 p=1, θ =π/4 10 i −10 10 5 p=2, θi=π/4 7 10 15 20 λ/h 40 80 Figure 2.31: h−, and p−convergence of the phase error per wavelength in degrees for VGFEM with Legendre polynomials and hat-PU function. p−convergence of the phase error of VGFEM with hat-PU function shown in Fig. 2.31 is exactly the same as that of scalar GFEM for axial incidence. Interestingly, the slope of the convergence is different for diagonal incidence and it is reduced by a factor that is dominant for p = 1. Next, the incident angle dependency is shown in detail in Fig. 2.32(b) for p = 1 and p = 2 with h = λ/10. The shape of dispersion curve of VGFEM for p = 1 is very similar to that of scalar GFEM given in Fig. (2.29), where maximum phase error occurs for axially incidence wave. However, the maximum dispersion for p = 2 appears for diagonal incidence unlike that of scalar GFEM. Finally, we investigate the dispersion in VGFEM using the HO-PU function. Fig. 2.32(a) illustrates the h−, and p−convergence of the VGFEM for both PU functions. h−convergence plots of both HO-PU and hat-PU functions are identical for p = 1, whereas HO-PU con65 2 Phase error [degrees] 10 0 10 −2 10 p=1, hat−PU p=2, hat−PU p=1, HO−PU p=2, HO−PU −4 10 −6 10 5 7 10 15 20 λ/h 40 80 (a) h−, and p−convergence for axial incidence (b) Incidence angle dependency, h = λ/10 Figure 2.32: h−, and p−convergence and incident angle dependency of the phase error per wavelength in degrees for VGFEM. 66 siderably suppresses the phase error for p = 2. Fig. 2.32(b) compares the incident angle dependency of the dispersion for HO-PU function against that for hat-PU function. HO-PU function suppresses the phase error more for p = 1 for near diagonal incidences, but overall the shapes of the curves are similar. For p = 2, the shapes of the curves are, however, completely different. While the maximum error is seen at diagonal incidence for HO-PU function, it is observed at axial incidence for hat-PU function. Overall, HO-PU function considerably reduces dispersion. 2.4 Summary In this Chapter, we have introduced several modifications that lays the foundation for using VGFEM to analyze more complex practical problems. Principally, the mathematics necessary to use arbitrary non-canonical partition of unity domains have been developed together with the manner in which the PU and local approximation functions can be defined on these domains. The new framework of VGFEM on brick elements has been tested via the analysis of various closed and open domain problems. To apply VGFEM to open domain problems, we have developed methodologies to integrate VGFEM with the perfectly matched layers (PML), absorbing boundary conditions, and boundary integrals. We have proposed techniques to include discontinuity in VGFEM, namely additional basis functions and domain decompositions. Both techniques have been tested by simulating the wave propagations in partially filled space and in rectangular waveguide with PML truncation. It has been shown that VGFEM can accurately capture the fields in PML region using only 2 or 3 layers at high attenuation factors. The hybridization of VGFEM with BI is achieved by enforcing the volumetric basis functions to satisfy the boundary conditions. The hybrid VGFEM-BI is 67 applied to solve various scattering problems of cavity backed apertures In addition, a semi-analytic technique has been developed to analyze dispersion characteristics of the scalar and vector GFEM. We have validated the proposed method for a range of practical problems against existing data, and demonstrated excellent agreement. These improvements permit VGFEM to operate in either a “meshless” environment or a “meshed” environment, or a mixture of both. The framework proposed in this Chapter can permit an easy integration of the method with classical hp-adaptive FEM. 68 Chapter 3 Tetrahedral Based VGFEM and its Applications In Chapter 2, the methodologies to adapt the meshless scheme of VGFEM onto brick elements have been presented. In this Chapter, this adaptation is further extended to tetrahedral meshes. This extension is not straightforward since a methodology needs to be developed in order to define the PU domains on tetrahedrons. In addition, the partition of unity functions and local approximation functions need to be redefined on this tessellation for construction of a vector basis set. Section 3.1 will present these developments and definitions along with a bilinear formulation of VGFEM. In addition to these developments, simulations performed will highlight some features of the method such as flexibility in the choice of basis functions, use of different types of basis functions or mixed polynomial orders within a simulation. In Section 3.2, a number of problems are analyzed within the presented framework to validate the developed methodology and demonstrate its efficacy. Finally, the work is summarized in Section 3.3. 69 Figure 3.1: An illustration of overlapping PU domains defined on a triangular mesh of the domain Ω. 3.1 3.1.1 Formulation of Tetrahedral Based VGFEM Partition of Unity Domains and Associated Basis Functions Conventional finite element solution to vector wave equation requires two fundamental steps: discretization of the problem domain and definition of appropriate basis functions on these tessellation. In VGFEM, an additional step is necessary: the construction of partition of unity domains on a given mesh. To this end, we develop (i) an appropriate methodology for this construction, and (ii) basis functions that are defined on these partitioned domains. Partition of unity (PU) domains and associated basis functions are constructed on tetrahedral meshes following the general methodology proposed in Chapter 2 . Consider a tetrahedral mesh with N number of nodes and Ne number tetrahedrons. Each PU domain Ωi is constructed from the tetrahedrons that are connected to the node i such that Nei Ωi = l=1 70 ¯ Ωe(i,l) , (3.1) where Nei is the total number of the tetrahedrons sharing the node i, l is the local index of each of these tetrahedrons, e(i, l) is the array containing the global indexes of the tetra¯ ¯ hedrons (j = e(i, l)), and Ωe(i,l) = Ωj is the domain of the tetrahedron j. For simplicity, an illustration of PU domains is depicted in Fig. 3.1 for a triangular mesh, where the PU domains for node 7 and 12 are explicitly shown. The PU domains Ωi satisfy the point-wise overlapping condition given in [25] such that Ω = N Ω. i=1 i Next, vector basis space given in Chapter 2 is defined on this tessellation by carefully designing partition of unity function ψi (r) and local approximation function vi (r). While defining PU function is straightforward in a regular shaped PU domain such as cuboid and spheroid, it is not so for an arbitrary shaped PU domain since the PU function ψi (r) has to satisfy the following conditions; (i) be Lipschitz continuous in Ωi , (ii) vanish at the domain boundaries ∂Ωi , and (iii) N ψ (r) = 1, r ∀ ∈ Ω. Satisfying all these conditions i=1 i on arbitrary PU domains is very difficult; therefore, we define the PU domains in simplex coordinates for tetrahedral meshes. A methodology to construct ψi (r) on a tetrahedral mesh is described next. Without loss of generality, the PU function can be constructed using any Lipschitz continuous function Wi (r) by ψi (r) = provided Wi (r) , ∀ Ωl ∩ Ωi l Wl (r) (3.2) l Wl (r) > 0 and Wi (r) smoothly vanishes at the PU domain boundaries ∂Ωi [18], where l denotes the index of the PU domain covering the point r for r ∈ Ωi . The function 71 Wi (r) is first defined in the simplex coordinates using Nei Wi (ξ) = l=1 ¯k We(i,l) (ξk ), (3.3) where k ∈ {1, 2, 3, 4} is the local index of the simplex coordinate ξk that corresponds to the ¯ ¯ ¯ (ξ ) is a function of ξk defined in Ωj . And node i in Ωj (remember j = e(i, l)), and W k e(i,l) k then, it is mapped to the global coordinates via the transformation 4 ξ k rk , r= (3.4) k=1 where r is any point in the element, rk is the node position at the k th vertex of the ele4 ¯k k=1 ξk = 1. The scalar function We(i,l) (ξk ) determines the order of the PU ¯ function. For instance, a linear PU function is simply constructed by W k (ξ ) = ξk . Ale(i,l) k ment. Note, ternatively, PU functions that have higher order continuous derivatives can be easily defined ¯ (ξ ). For illustration purposes, a 2D-PU function defined over eight by redefining W k e(i,l) k triangles sharing the node at the center is shown in Fig. 3.2(b). Next, local approximation functions are defined on each local domain Ωi , and chosen p from a space of functions vi (r) ∈ span vi (r) , where p is the local index of local approximation function. The local approximation functions can be defined in the physical domain as illustrated in Fig. 3.2(a). Since the vector basis function is a product of the PU and the local approximation functions, the domain of vector basis is fully determined by the PU lap domain. Consequently, the domain of the local approximation function, denoted by Ωi , lap can be constructed such that Ωi ⊂ Ωi . Local approximation functions based on Legendre polynomials and plane waves will be used in this Chapter. For Legendre polynomials, a 72 (a) Geometry map 0.5 i ψ (r) 1 0 1 0.5 1 y 0.5 0 0 −0.5 −0.5 −1 x −1 (b) PU function Figure 3.2: Illustration of the domains of VGFEM functions, and a plot of the PU function.. 73 cuboid is needed for the definition of the functions. Thus, first the size of the cuboid is lap found such that Ωi ⊂ Ωi , and then center it at the PU node position ri = (xi , yi , zi ). The dimensions of cuboid for each Ωi can be easily determined by Nei 4 li = 2 max( l=1 k=1 |ri − rn (k) |), f or j = ei (l), j (3.5) where nj is an array having the global indexes of the nodes of j th tetrahedron, li = (l1 , l2 , l3 ) is the dimensions of the cuboid, and the operator max applies for each dii i i mension separately. For plane waves, it makes sense to choose the phase center at ri . Since PU functions fully control the support of basis functions, and provide smooth transitions through adjacent PU domains, one has flexibility in the choice of local approximation functions. For instance, one can use Legendre polynomials or plane waves as local approximation functions. Moreover, one can mix these basis functions within a simulation, or one can use different basis functions or polynomial orders in different PU domains. This can be easily achieved via assignment of a basis type or polynomial order to each PU domain. 74 3.1.2 Bilinear Formulation of VGFEM With these preliminaries, a bilinear formulation of VGFEM is presented next along with the boundary conditions u(r) × n = gd (r), f or r ∈ Γd , ˆ n× ˆ 1 ∇ × u(r) = gn (r), f or r ∈ Γn , µr n × ∇ × u(r) + jk0 n × n × u(r) = 0, f or r, ∈ Γabc , ˆ ˆ ˆ (3.6) where n is the unit outward vector normal to the boundary, Γd , Γn , and Γabc denote ˆ the Dirichlet, Neumann, and absorbing boundaries, respectively, and gd (r) and gn (r) are functions corresponding to these boundary conditions. With these boundary conditions, a bilinear formulation of VGFEM can then be written as Ωi β β1 Γd Ωi 2 (∇ × w) · (∇ × u) dΩ − k0 (∇ · w)(∇ · u)dΩ + (w × n) · (u × n)dΓ + ˆ ˆ jk0 β1 Γd (w × n) · gd dΓ + ˆ Γd Γd Ωi w · u dΩ+ (w × n) · (∇ × u)dΓ+ ˆ Γd (∇ × w) · (u × n)dΓ + ˆ Γabc (w × n) · (u × n)dΓ = ˆ ˆ (∇ × w) · gd dΓ − Γn w · gn , dΓ (3.7) where β and β1 are constants used to obtain uniform convergence [25]. In this Chapter, the scattered field formulation is used for the analysis of open domain problems; thus, the 75 boundary condition imposed on a perfect electrically conducting (PEC) scatterer is gd (r) = −uinc (r) × n, where uinc (r) denotes the incident electric field. Note, ˆ = n is assumed in ˆ r ˆ the fist order absorbing boundary condition (ABC) formulation as in [13]. Some of the problems presented in this Chapter include analysis of field scattering from cavity backed apertures. Boundary integrals together with VGFEM are used to analyze these problems. Hybridization of VGFEM with boundary integrals (VGFEM-BI) is achieved by enforcing the volumetric basis functions to satisfy the boundary conditions at the aperture as described in Chapter 2. The bilinear formulation of the VGFEM-BI presented there is also valid for tetrahedrons. For the sake of brevity, we only state the differences that appear in the evaluation of matrix elements here. For VGFEM-BI with brick elements, we have divided each rectangular element into two triangles, and then computed the singular surface integrals using the quadrature rule developed in [30]. Since the surface element in tetrahedral meshes is a triangle, this step is eliminated in the computation of the singular integrals. 3.2 Numerical Results This Section presents a number of results that validate the VGFEM based on tetrahedral elements. The focus of numerical studies will be to demonstrate capabilities of VGFEM: Specifically, p−refinement, and the ability to construct an approximation comprising of different basis functions (polynomials and non-polynomials) without building additional mechanisms to satisfy continuity conditions. These features will be investigated via applications to two classes of problems: (i) wave propagation in various waveguide structures, and (ii) scattering from PEC scatterers and cavity backed apertures. Computational cost of VGFEM was presented Chapter 2 for brick elements. As the scaling is similar, the cost analysis is 76 Table 3.1: Eigenvalues for an empty 1.0 m × 0.5 m × 0.75 m rectangular cavity Mode T E101 T M110 T E011 T E201 T M111 T E111 T M210 T E102 Eigenvalues (h = 0.153) Analytic V GF EM 5.236 5.237 7.025 7.029 7.552 7.556 7.552 7.560 8.179 8.191 8.179 8.194 8.886 8.903 8.947 8.964 −h convergence, Error(%) h = 0.153 0.021 0.059 0.068 0.112 0.146 0.183 0.199 0.182 h = 0.225 0.092 0.170 0.229 0.451 0.323 0.479 7.423 0.165 h = 0.304 0.46 0.49 1.09 1.55 1.27 1.74 4.76 3.66 not repeated here. In this Section, p denotes the order of local approximation in the vector space, and Legendre polynomial orders of p + 1 and p + 2 are used to get pth order approximation. Without loss of generality, the pilot vector is chosen as c = z . In what follows, Nb ˆ ˆ denotes the number of basis functions in each PU domain, and h denotes the average edge length of tetrahedrons. To solve the matrix system, Transpose-Free Quasi-Minimal Residual (TFQMR) iterative solver with incomplete LU (ILU) preconditioner is utilized. The drop tolerance of 10−5 for ILU preconditoner and error tolerance of 10−4 for TFQMR are used in all simulations unless otherwise is stated. 3.2.1 Convergence Analysis We validate the framework of VGFEM defined on tetrahedral elements through cavity mode analysis and convergence study of wave propagations. The hat PU function and Legendre polynomials of various orders are used. First, we study the eigenvalues of a 1.0m × 0.5m × 0.75m rectangular cavity using local approximation of p = 1. First eight eigenvalues computed using average edge length of h = 0.153m are compared against the analytic ones in Table 3.1. Eigenvalues are accurately captured with percentage errors less 77 than 1% for all eight modes as shown in Table 3.1. We also show h− convergence of VGFEM for those eigenvalues in percentage error in the same Table. Next, we investigate h, and p−convergence of the VGFEM for wave propagation problems. A computation domain of 1.0m×0.5m×0.75m, and four different meshes with average edge lengths of h = {0.75m, 0.304m, 0.153m, 0.105m} of the domain are considered for this study. Fig. 3.3(a) shows h−, and p−convergence of VGFEM for T E10 mode propagation in this domain at λ = 1.8m. Fig. 3.3(b) shows the convergence plots for a plane wave propagating with an incident angle of (θinc , φinc ) = (π/4, 0) at λ = 1 m. Neumann BC is imposed at the domain boundaries for the plane wave problem. In Fig. 3.3, we also show the FEM results for linear edge elements. VGFEM with p = 0 exhibits similar convergence behavior as FEM does. VGFEM shows excellent h−, and p− convergence for both problems. 3.2.2 Application of VGFEM to Waveguide Problems In this subsection, we demonstrate capabilities of VGFEM analyzing wave propagations in various hollow waveguides. Boundary conditions are imposed at the input port, output port(s), and on the walls of waveguides as described in Chapter 2, and scattering parameters are computed using the formulation given there. First, we study the use of different orders (p) of polynomial basis functions in different regions of a simulation domain. WR-90 S-bend waveguide with a = 22.9mm, b = 10.2mm, c1 = 2a, L = 25mm, r = 15.24mm, and α = 30◦ as shown in Fig 3.4(a) is considered for this study. The waveguide is meshed using 1067 nodes and 4075 tetrahedrons with the average edge length of h = 4.23mm. Two different simulations are performed using Legendre 78 0 10 −1 Relative L2 error 10 −2 10 FEM −3 10 VGFEM, p=0 VGFEM, p=1 −4 VGFEM, p=2 10 −5 10 0 10 1 2 10 10 3 4 10 Unknowns 10 5 10 (a) T E10 mode in rectangular waveguide. 0 10 −1 Relative L2 error 10 −2 10 FEM VGFEM, p=0 VGFEM, p=1 −3 10 VGFEM, p=2 −4 10 1 10 2 10 3 10 Unknowns 4 10 5 10 (b) Plane wave propagation with Neumann BC. Figure 3.3: h− and p−convergence of VGFEM in 1.0m × 0.5m × 0.75m. 79 polynomials. While p = 1 is used in the entire domain for the first simulation, p = 1 in the straight sections and p = 2 in the bend sections as shown in 3.4(a) are used for the second simulation. Figure 3.4(b) shows the excellent agreement among the VGFEM results and the MoM result given in [1]. Next, we demonstrate another capability of VGFEM that is the flexibility in the choice of basis functions. Plane wave and Legendre polynomial basis functions are chosen for this demonstration. Local approximation order of p = 2 is used for Legendre polynomials such that Nb = 26. The space of plane wave basis functions is chosen to be span{e−jkp ·(r−ri ) }, where ri is the position of the node i that Ωi is defined around, and kp = (k0 , θ, φ) is the direction of the plane wave basis functions. Four particular plane wave directions (θ, φ) ∈ {(π/4, 0), (π/4, π), (π/2, 0), (π/2, π)} are used such that Nb = 8. The basis set includes the projection on the principal mode of the rectangular waveguide, and this will give more accurate results by suppressing the dispersion significantly. The waveguide is meshed using 617 nodes and 1893 tetrahedrons with h = 5.23mm. Two different simulations are performed. For the first simulation, plane wave basis functions are used in the entire region. For the second simulation, while plane wave basis functions are used in the straight sections, Legendre polynomials with the local approximation of p = 2 are used in the bend sections. As shown in Fig. 3.5, both simulation results are in full agreement with MoM data given in [1]. We further demonstrate the use of mixed basis functions and p−refinement within a simulation analyzing the scattering in a H -plane WR-75 waveguide. The waveguide dimensions are a = 19.05mm, b = 9.525mm, c1 = 2a, r = 10.7mm, and the bend angle of α = 180◦ . The top view of the waveguide geometry is shown in Fig. 3.6. The waveguide is meshed using 997 nodes and 3228 tetrahedrons with h = 4.09mm. Legendre polynomials with dif80 (a) Geometry of WR-90 S-bend waveguide. −20 −30 | S11| [dB] −40 −50 −60 VGFEM, p=1 VGFEM, p=1 and p=2 MoM* −70 −80 8 8.5 9 9.5 10 Frequency [GHz ] 10.5 11 11.5 (b) S11 of WR-90 S-bend waveguide Figure 3.4: S11 of WR-90 S-bend waveguide is simulated using mixed polynomial orders and compared against the MoM data in [1]. 81 −20 −30 | S11 | [dB] −40 VGFEM, Plane waves VGFEM, Mixed basis MoM* −50 −60 −70 −80 8 8.5 9 9.5 10 Frequency [GHz ] 10.5 11 11.5 Figure 3.5: S11 of WR-90 S-bend waveguide is simulated using mixed basis functions and compared against the MoM data in [1]. ferent orders and modal basis functions are used. The modal basis function, which is the T E10 mode, is used around the ports, c1 = 1.0a, as depicted in Fig. 3.6, where higher order modes are assumed to vanish. In the rest of the domain, Legendre polynomials with the local approximation order p = 1 (Nb = 11) and p = 2 (Nb = 26) are used. The polynomial order p = 2 is used next to the modal basis for a good coupling of basis functions. Simulated results agree well with the MoM data given in [6] as shown Fig. 3.6 despite the use of a coarse mesh. Finally, use of mixed basis functions and different polynomial orders within VGFEM framework are further investigated simulating the scattering parameters of an H-plane T junction. The waveguide dimensions are a = 15.799mm, b = 7.899mm, and c = 2a as shown in Fig. 3.7. The waveguide is meshed using 1588 nodes and 6477 tetrahedrons with 82 Figure 3.6: |S11 | of H plane WR-75 waveguide with 180◦ bend is simulated using mixed polynomial orders and compared against the MoM data in [6] h = 2.78mm. Two types of basis functions are used: Legendre polynomials and plane waves. The space of plane wave basis functions is the same as the one used in the W R − 90 waveguide problem. The plane wave basis functions are particularly chosen for an accurate representation of the fundamental mode; thus, they are used around the port 1 and port 2, c1 = 1.0a. Legendre polynomials with the local approximation order p = 1 and p = 2 are used in the remaining regions. Figure 3.7 compares absolute values of the scattering parameters at the ports against the reference ones given in [7]. Simulated parameters are in agreement with the reference data. 3.2.3 Application of VGFEM to Open Domain Problems In this subsection, the VGFEM is applied to a set of open domain scattering problems. For domain truncation, either absorbing boundary condition or boundary integrals is used. First, 83 Figure 3.7: Simulated scattering parameters of H-plane T-junction are compared against the reference data in [7]. 84 Figure 3.8: Bistatic RCS of a 0.3m × 0.3m PEC plate. The plate is illuminated by the electric field E = xejk0 z ˆ the bistatic RCS of a 0.3m × 0.3m PEC plate is simulated using the first order ABC and the local approximation order of p = 2. The plate is centred at the origin on xy plane, and it is illuminated by the electric field E = xejk0 z at 300MHz. A cubical truncation boundary ˆ is placed 0.3m away from the edge of the plate. The problem domain is meshed using 2493 nodes and 12312 tetrahedrons with h = 0.1m. Figure 3.8 shows excellent agreement between the simulated RCS and the MoM data. Next, the bistatic RCS of a sphere with radius of r = 0.5λ is simulated for λ = 1.0m to further validate the ABC formulation. The sphere is centred at the origin, and it is illuminated by the electric field E = xejk0 z . A spherical truncation boundary is placed 0.5λ ˆ away from the sphere, and the first order ABC is imposed on the boundary. The problem 85 20 Mie Series HFSS, havg=λ/10, p=1 15 VGFEM, h =λ/5, p=1 and p=2 avg 10 σ [dBsm] 5 0 −5 −10 −15 −20 0 30 60 90 120 θ [degrees] 150 180 Figure 3.9: Bistatic RCS of a PEC sphere with radius of r = 0.5λ. The sphere is illuminated by the electric field E = xejk0 z . ˆ domain is meshed using 531 nodes and 1762 tetrahedrons with the average edge length of h = λ/5. While p = 1 is used in the interior PU domains associated with the interior nodes, p = 2 is used in the boundary PU domains associated with the boundary nodes. Figure 3.9 shows the simulated RCS results obtained using VGFEM and the commercial software HFSS and the analytic RCS obtained using Mie Series. Note that HFSS uses 2nd order ABC and, simulation parameters used are h = λ/10 and p = 1. VGFEM with varying p gives accurate enough results despite the coarse tetrahedral mesh. Thereafter, we present the application of VGFEM with boundary integrals to some scattering from cavity backed apertures. Apertures are assumed be in an infinite ground plane at z = 0, and cavity regions are considered in z < 0. RCS of the apertures are computed using the far field formulation given in [26]. First, the monostatic RCS of the aperture backed 86 10 5 0 σ / λ2 [dB] −5 −10 −15 Co−pol, FEM−BI X−pol, FEM−BI Co−pol, VGFEM−BI X−pol, VGFEM−BI −20 −25 −30 −35 −40 0 10 20 30 40 50 θ [degrees] 60 70 80 90 Figure 3.10: Backscattered RCS of the cavity-backed aperture with a = 0.7λ, b = 0.1λ, and c = 1.73λ compared with FEM-BI results [2] ˆ ˆ by 0.7λ × 0.1λ × 1.73λ cavity is simulated. A θ polarized plane wave, Ei = θ e−jk·r with θi = 40◦ is incident to the empty cavity. The cavity region is meshed using 232 nodes and 644 tetrahedrons with the average edge length of h = λ/8. The order of local approximation is p = 1. Figure 3.10 shows excellent match between RCS results of VGFEM-BI and FEM-BI [2] for both polarizations. Next, the monostatic RCS of the aperture backed by 2.5λ × 0.25λ × 0.25λ cavity is simulated for both E-polarized and H-polarized incident fields. The cavity region is meshed using 249 nodes and 784 tetrahedrons with the average edge size of h = λ/8. The order of local approximation is p = 1. VGFEM-BI and FEM-BI results [4] for φi = 90◦ agree very well as shown in Fig. 3.11(a). Figure 3.11(b) shows the RCS results obtained using tetrahedral elements and brick elements [26] for φi = 0◦ . The results are in excellent agreement for both 87 20 10 0 σ / λ2 [dB] −10 FEM−BI, H−pol FEM−BI, E−pol VGFEM−BI, E−pol VGFEM−BI, H−pol −20 −30 −40 −50 −60 0 10 20 30 40 50 θ [degrees] 60 70 80 90 (a) φinc = 90◦ 20 10 0 2 σ / λ [dB] −10 −20 −30 −40 −50 VGFEM−BI, H−pol, brick VGFEM−BI, E−pol, brick VGFEM−BI, E−pol, tetra VGFEM−BI, H−pol, tetra −60 −70 −80 0 10 20 30 40 50 θ [degrees] 60 70 80 90 (b) φinc = 0◦ Figure 3.11: Backscattered RCS of the cavity-backed aperture with a = 2.5λ, b = 0.25λ, and c = 0.25λ compared with FEM-BI results [4]. 88 0 Measurement VGFEM−BI, Legendre polynomials VGFEM−BI, Plane waves −10 σ / m2 [dB] −20 −30 −40 −50 −60 0 10 20 30 40 50 θ [degrees] 60 70 80 90 Figure 3.12: Backscattered RCS of the cavity-backed aperture with a = 16.26λ, b = 0.2λ, and c = 0.85λ compared against measurement result at 12 GHz [5]. polarizations. Next, RCS of a long cavity with the dimensions of 16.26λ × 0.2λ × 0.85λ is simulated at 12 GHz. The cavity region is meshed using tetrahedrons with an average edge length of h = λ/4. Fig. 3.12 shows the RCS results of VGFEM-BI for two different local approximation functions. The simulation with Legendre polynomials has been performed for p = 1 order of local approximation. For the other simulation, Nb = 10 number of plane wave basis functions in each PU domain is used. As can be seen in Fig. 3.12, both results match and, they are very close to the measurement results like those of FEM-BI obtained in [5]. Finally, the monostatic RCS of a large cavity with the dimensions of 2λ × 2λ × 10λ is simulated to show the suppression of the numerical dispersion in VGFEM. The cavity region is meshed using 7392 nodes and 37342 tetrahedrons with h = λ/5. The incident electric 89 40 30 20 σ / λ2 [dB] 10 0 VGFEM −10 FEM−BI, tetra FEM−BI, mixed −20 HFSS −30 −40 0 10 20 30 40 50 θ [degrees] 60 70 80 90 Figure 3.13: Backscattered RCS of the cavity-backed aperture with a = 2λ, b = 2λ, and c = 10λ compared against FEM-BI data given in [3] and HFSS data. field is φ polarized and propagates with φi = 0◦ . Figure 3.13 shows the monostatic RCS simulations for VGFEM and HFSS. It also shows the FEM-BI reference data given in [3], where 3rd order tetrahedral elements and mixed tetrahedral-prism elements were compared. HFSS simulation has been performed using adaptive option with p = 1 and initial mesh size of h = λ/10 such that the final mesh size is h = λ/13. Although VGFEM uses coarser mesh and p = 1, it suppresses dispersion and approximates the RCS accurately. 3.3 Summary In this Chapter, we have developed a VGFEM framework using tetrahedral elements, and validated it by simulating a number of open and closed domain problems. More importantly, 90 we have sought to highlight some features of VGFEM; principally the ability to use different polynomial orders or mixture of different basis functions as local approximating functions. Realization of these features is as simple as assigning a basis function type or a polynomial order to each PU domain. The framework proposed here can permit an easy hybridization of the method with classical hp−adaptive FEM and discontinuous Galerkin methods. Integration of VGFEM with these methods to create a highly flexible and powerful method is the topic of next Chapter. 91 Chapter 4 Hybrid Vector Generalized Finite Element Methods In this Chapter, we introduce hybridizations of VGFEM with Discontinuous Galerkin (DG) Methods, Domain Decomposition Methods (DDM), and classical FEM for accurate and, perhaps, efficient solutions of EM problems using the advantages of each method. The mathematical framework of VGFEM permits the use of different types of basis functions or mixed polynomial orders within a simulation. Likewise, discontinuous Galerkin (DG) methods enable handling of multi-material problems, flexibility in the mesh, and parallelization [11, 12]. DG methods relax the field continuity requirements across the finite elements and treat each finite element as an individual system. The coupling of these systems is achieved via jump and average conditions across the element boundaries. This makes the method very well suited for multi-material problems and parallelization [12]. The properties of VGFEM and DG are complementary each other. The former uses a space of functions that is continuous across the PU domains, and one has to built functionality to enforce discontinuity for 92 inhomogeneous problems. The latter is discontinuous across interfaces, and one needs to impose conditions to ensure continuity of the fields. We aim to hybridize DG and VGFEM to use the features of both methods for simulations of complex EM problems. For instance, one can partition a large and complex problem domain into simple domains, use VGFEM in each domain, and couple the fields using DG across domain boundaries. Similarly, DD methods enable the solutions of electrically large problems thanks to its divide-and-conquer philosophy [9, 10]. In this approach, the original problem is partitioned into smaller regions and each region is solved separately instead of solving a large and complex problem directly as a whole. This is achieved by using transmission conditions at the interfaces between adjacent sub-domains to enforce the continuity of electromagnetic fields. Furthermore, we can hybridize VGFEM and classical FEM using these domain decomposition methodologies for more accurate and perhaps efficient simulations. In this hybrid approach, a problem domain is decomposed into sub-domains; while FEM is used in some of these sub-domains, VGFEM is utilized in the others, and then the methods are coupled via the a boundary conditions imposed at the sub-domain interfaces. This hybrid technique enables the utilization of advantages of each method within a simulation. For instance, FEM can be used at domain boundaries or dielectric interfaces to impose the boundary conditions easily, and VGFEM can be applied to the regions where more accurate approximations are needed. These hybridizations will be realized by first developing a technique to define the framework of VGFEM in partitioned domains, and then by applying appropriate transmission conditions to couple the methods. Imposition of transmission conditions between sub-domains is not straightforward as vector basis function space in VGFEM is not interpolatory and no auxiliary surface basis is allowed in the solution space. However, they can be imposed 93 by defining PU domains centred at the shared boundary of adjacent sub-domains such that basis functions on these PU domains that are continuous throughout the boundary satisfy both Dirichlet and Neuman transmission conditions. Once we develop the technique for PU definitions in partitioned domains, we can hybridize VGFEM with FEM, DDM and DGM using the related transmission conditions at the interfaces. We will validate the proposed hybrid techniques analyzing wide range of practical EM problems. The Chapter is organized as follows: First, a technique to define the framework of VGFEM in partitioned domains is given in Section 4.1, and then bilinear formulations along with DDM and DGM transmission conditions are presented. The effect of transmission conditions on convergence are studied, and VGFEM with DDM and DGM coupling conditions are validated via simulations of various problems. In Section 4.2, VGFEM is hybridized with classical FEM using DGM and DDM coupling conditions. The hybrid method is applied to various open and closed domain problems. Finally, conclusions and future research are presented in Section 4.3. 4.1 Domain Decomposition for VGFEM Consider a linear, homogeneous, and source free three-dimensional region Ω, whose boundary is denoted by ∂Ω := Γ = i Γi . In Ω, the electric field, u(r), satisfies the vector wave equation ∇× 1 2 ∇ × u(r) − k0 ǫr (r)u(r) = 0 µr (r) Bi {u(r)} = gi (r), r ∀ ∈ Γi , 94 (4.1) Figure 4.1: An illustration of material interfaces and domain decomposition where ko is the free space wavenumber, ǫr (r) and µr (r) are relative constitutive parameters, Bi is a differential operator of boundary conditions, and gi (r) is a function of boundary conditions imposed on Γi . Following the methodologies presented for inhomogeneous problems in Chapter 2, we partition the problem domain into sub-domains, and define the PU domains on these partitions. We will use the same technique to define the PU domains, assuming that the PU domains are defined around each node as a union of the elements sharing this node. The PU functions are then defined on these PU domains as described in Chapter 2 for bricks and in Chapter 3 for tetrahedrons. This definition can be easily extended to any type of elements or mixture of of them. Next, we will present an application of domain decomposition for inhomogeneous problems. Let us consider various configurations of material interfaces shown in Fig. 4.1 within the problem. In order to use VGFEM basis functions to solve this problem, the domain is partitioned into homogeneous ones, and the PU domains and functions are defined as described above. Once the PU domains are defined on partitioned domains, transmission conditions are applied at the domain interfaces. For simplicity and without lost of generality, let us consider a problem domain partitioned into two sub-domains as Ω1 and Ω2 , and let u1 and u2 denote the electric fields in these regions respectively. Next, we formulate the 95 transmission conditions for three different approaches: 1. First approach is to impose the Robin transmission condition at the coupling interfaces Γ12 and Γ21 as 1 1 ∇ × u1 ) + βrb n1 × u1 × n1 = −ˆ 2 × ( ˆ ˆ n ∇ × u2 ) + βrb n2 × u2 × n2 ˆ ˆ n1 × ( ˆ µr1 µr2 (4.2) at the interface Γ12 , and 1 1 n2 × ( ˆ ∇ × u2 ) + βrb n2 × u2 × n2 = −ˆ 1 × ( ˆ ˆ n ∇ × u1 ) + βrb n1 × u1 × n1 ˆ ˆ µr2 µr1 (4.3) √ at the interface Γ21 , where βrb = jk0 ǫr µr . ¯ ¯ 2. Second approach is to impose tangential continuity of the electric field and normal continuity of the flux at the coupling interfaces Γ12 and Γ21 . The transmission conditions are formulated as βtc1 n1 × u1 × n1 + βtc2 n1 (ˆ 1 · ǫr1 u1 ) = βtc1 n2 × u2 × n2 + βtc2 n2 (ˆ 2 · ǫr2 u2 ) ˆ ˆ ˆ n ˆ ˆ ˆ n (4.4) at the interface Γ12 , and βtc1 n2 × u2 × n2 + βtc2 n2 (ˆ 2 · ǫr2 u2 ) = βtc1 n1 × u1 × n1 + βtc2 n1 (ˆ 1 · ǫr1 u1 ) ˆ ˆ ˆ n ˆ ˆ ˆ n (4.5) at the interface Γ21 , where βtc1 and βtc2 are the penalty factors. 3. Third approach is to use the interior penalty discontinuous Galerkin (IPDG) technique 96 for coupling at the domain interface Γdg [12]. The boundary conditions imposed at Γdg in Ω1 are n1 × u1 = (ˆ 1 × u1 + n1 × u2 )/2, ˆ n ˆ 1 1 ∇ × u1 + n1 × ˆ ∇ × u2 )/2, n1 × ∇ × u1 = (ˆ 1 × ˆ n µr1 µr2 (4.6) Likewise, the boundary conditions imposed at Γdg in Ω2 are n2 × u2 = (ˆ 2 × u1 + n2 × u2 )/2, ˆ n ˆ 1 1 n2 × ∇ × u2 = (ˆ 2 × ˆ n ∇ × u1 + n2 × ˆ ∇ × u2 )/2, µr1 µr2 (4.7) By testing (4.6) and (4.7) by using the testing functions n × u and u, and then ˆ summing the resultant equations, one can write the final bilinear formulation for the boundary conditions at Γdg as Γdg 1 βdg [[w]]T · [[u]]T + [[w]]T · { ∇ × u} dΓ = 0 µr (4.8) where βdg is the penalty factor, [[w]]T = w2 ׈ 2 +w1 ׈ 1 with n2 = −ˆ 1 is the jump n n ˆ n 1 operator, and { µr ∇×u} = ( µ1 ∇×u1 + µ1 ∇×u2 )/2 is the average operator. If DG r1 r2 boundaries are dielectric material interfaces, we also impose the boundary condition on the normal components of the electric field βdg Γdg [[ǫr w]]N [[ǫr u]]N = 0 (4.9) where [[ǫr w]]N = ǫr2 w2 · n2 + ǫr1 w1 · n1 . ˆ ˆ Next, bilinear formulations for the second and third approaches are presented along with the 97 coupling boundary conditions formulated above. The formulation for the first approach is extensively studied for DDM [10]; therefore, it is not repeated here. For the second approach, the bilinear formulation is written for Ω1 as Ω1 (∇ × w1 ) · ( 1 2 ∇ × u1 ) dΩ − k0 ǫr1 w1 · u1 dΩ + w1 · B1 {u1 } µr1 Ω1 ∂Ω1 \Γ12 βtc1 βtc1 Γ12 Γ12 w1 × n1 · u1 × n1 dΓ + βtc2 ˆ ˆ w1 × n2 · u2 × n2 dΓ + βtc2 ˆ ˆ Γ12 Γ12 (w1 · n1 )(ǫr1 u1 ) dΓ = ˆ (w1 · n2 )(ǫr2 u2 ) dΓ + ˆ ∂Ω1 \Γ12 w 1 · gi (4.10) and for Ω2 as Ω2 (∇ × w2 ) · ( 1 2 w2 · u2 dΩ + w2 · B2 {u2 } ∇ × u2 ) dΩ − k0 ǫr2 µr2 Ω2 ∂Ω2 \Γ21 βtc1 βtc1 Γ21 Γ21 w2 × n2 · u2 × n2 dΓ + βtc2 ˆ ˆ w2 × n1 · u1 × n1 dΓ + βtc2 ˆ ˆ Γ21 Γ21 (w2 · n2 )(ǫr2 u2 ) dΓ = ˆ (w2 · n1 )(ǫr1 u1 ) dΓ + ˆ ∂Ω1 \Γ21 w 2 · gi (4.11) where n1 and n2 are outward normal vectors of Γ12 and Γ21 respectively. Likewise, for the ˆ ˆ third approach, the bilinear formulation can be written as 1 2 ∇ × u) dΩ − k0 ǫr2 w · u dΩ + w · B{u} µr2 Ω ∂Ω\Γdg Ω 1 [[ǫr w]]N [[ǫr u]]N βdg [[w]]T · [[u]]T + [[w]]T · { ∇ × u} dΓ + βdg µr Γdg (∇ × w) · ( Γdg = 98 ∂Ωl \Γdg w · gi (4.12) Accuracy and convergence characteristics of these coupling conditions will be studied numerically next. 4.1.1 Numerical Results In this Section, we study the convergence characteristics of hybrid DD-VGFEM and DGVGFEM, and test the methods through applications to various problems. The advantageous of VGFEM are used in some simulations such as local p-refinement and use of mixed basis functions, and it is shown that these features are applicable within the hybrid framework. 4.1.1.1 Convergence Analysis First, we study the effect of βtc on the convergence of DD-VGFEM. T E10 mode propagation in a rectangular waveguide with the dimensions of a × 0.5a × 2a is considered for this study. The domain is meshed using brick elements with average element length of h = 0.25a the order of local approximation is p = 1. Figure 4.2 shows the convergence plots for two different a values. For the first simulation with a = 1m and λ0 = 1.8 m, the error converges for βtc 100 for both two and four domain partitions (evenly partitioned along z) as shown in Fig. 4.2(a). As we scale the geometry to a = 0.01m and perform the simulations at λ0 = 0.018 m, we see the convergence for the βtc 1000 as shown inFig. 4.2(b). This simulation clearly shows the h dependency of the βtc . Next, we analyze the convergence of the hybrid method as the number of partitions is increased. The same problem is considered for a = 1m and λ0 = 1.8 m. Simulation parameters are h = 0.25m, βtc = 103 . The domain is evenly partitioned along z. Figure 4.3 shows the convergence plots for two different number of partitions. We can see that the 99 −1 10 2 subdomains 2 subdomains ( βtc= − j βtc) Relative L2 error 4 subdomains −2 10 −3 10 0 10 2 10 4 β 6 10 10 tc (a) a = 1m and λ0 = 1.8 m 0 10 2 subdomains 2 subdomains ( β = − j β ) tc tc Relative L2 error 4 subdomains −1 10 −2 10 −3 10 0 10 2 10 4 β 10 6 10 tc (b) a = 0.01m and λ0 = 0.018 m Figure 4.2: Effect of coupling coefficient βtc on convergence. The problems is T E10 mode propagation in a rectangular waveguide with the dimensions of a × 0.5a × 2a. 100 0 10 0 −2 Relative L error 10 −2 10 2 2 Relative L error 10 −4 10 VGFEM, p=0 VGFEM, p=1 VGFEM, p=2 VGFEM, p=3 −6 10 0 10 10 −6 10 1 λ0 / h VGFEM, p=0 VGFEM, p=1 VGFEM, p=2 DD−VGFEM, p=0 DD−VGFEM, p=1 DD−VGFEM, p=2 −4 10 0 10 λ0 / h (a) 1 subdomain (b) 2 subdomains 0 0 10 Relative L error 10 −2 10 −2 10 2 Relative L2 error 1 10 VGFEM, p=0 VGFEM, p=1 VGFEM, p=2 DD−VGFEM, p=0 DD−VGFEM, p=1 DD−VGFEM, p=2 −4 10 −6 10 0 10 10 −6 10 1 λ0 / h VGFEM, p=0 VGFEM, p=1 VGFEM, p=2 DD−VGFEM, p=0 DD−VGFEM, p=1 DD−VGFEM, p=2 −4 0 10 10 1 λ /h 10 0 (c) 3 subdomains (d) 4 subdomains Figure 4.3: Convergence of DD-VGFEM for different number of partitions. The problems is T E10 mode propagation in a rectangular waveguide with the dimensions of 1m × 0.5m × 2m. 101 −1 Relative L2 error 10 2 subdomains 4 subdomains 8 subdomains 16 subdomains −2 10 −3 10 0 10 2 10 4 β 10 6 10 dg Figure 4.4: Effect of DG penalty factor on the wave approximation. 1m × 0.5m × 2m rectangular waveguide is partitioned along z axis error arising at the domain interfaces increases as the number of partitions rises. Next, we study the effect of βdg on the convergence of the DG-VGFEM. T E10 mode propagation in a rectangular waveguide with the dimensions of 1m×0.5m×2m is considered for this study. The domain is meshed using brick elements with average element length of h = λ0 /6, where λ0 = 1.8m. The order of local approximation is p = 1. Figure 4.4 shows the convergence plots for various domain partitions. We get the convergence for βdg 1000 for all partitions for this problem. Next, we study the convergence characteristics of DG-VGFEM for T E10 mode propagation in a 1.0λ0 × 0.5λ0 × 2.0λ0 rectangular waveguide. Problem domain is partitioned into four identical sub-domains along longitudinal dimension and each of them is meshed using conformal brick elements. Legendre polynomials with different orders, and DG penalty 102 0 −2 10 2 Relative L error 10 −4 10 DG−VGFEM, p=0 DG−VGFEM, p=1 DG−VGFEM, p=2 DG−VGFEM, p=3 −6 10 0 10 1 λ0 / h 10 Figure 4.5: h− and p− convergence of VGFEM with DG transmission conditions for T E10 mode propagation in a 1.0λ0 × 0.5λ0 × 2.0λ0 rectangular waveguide. Problem domain is partitioned into 4 sub-domains along longitudinal dimension. factor of βdg = 103 are used for this study. Note, one needs to use Legendre polynomial orders of p + 1 and p + 2 to get pth order local approximation. Figure 4.5 shows h−, and p−convergence of the method for λ0 = 1.8m. We should point out that these results are identical to those obtained using VGFEM without partitioning. The last simulation for the convergence study is obliquely incident plane wave propagation through an infinite space with a material boundary at x − y plane for θinc = π/6, φinc = 0 and λ0 = 1m. A computation domain of 1m × 0.1m × 1m centred at the origin is considered for this study, half of which is filled with ǫr = 4. Both Neumann and Dirichlet boundary conditions are imposed at the domain boundaries. DG conditions are applied at the dielectric-air interface to represent the discontinuity of the normal component of the electric field. Figure 4.6 shows h− and p−convergence of the method. 103 0 10 −1 Relative L2 error 10 −2 10 DG−VGFEM, p=0 DG−VGFEM, p=1 DG−VGFEM, p=2 DG−VGFEM, p=3 −3 10 −4 10 0 10 1 λ0 / h 10 Figure 4.6: h− and p− convergence of DG-VGFEM for an inhomogeneous problem. Plane wave with θinc = π/6 and φinc = 0 is propagating in a 1m × 0.1m × 1m domain with ǫr = 4 for z > 0m. 4.1.1.2 Applications Next, we will apply the hybrid methods to a number of wave propagation problems. First, we will present the DD-VGFEM results. DD-VGFEM is applied to wave scattering in WR-90 S-bend waveguide. Figure 4.7(a) shows the geometry of the problem, where a = 22.9mm, b = 10.2mm, c1 = 2a, L = 25mm, r = 15.24mm, and α = 30◦ . The geometry is meshed using tetrahedral elements with an average edge length of h = 6.6mm. The waveguide is partitioned into 3 sub-domains as shown in Fig. 4.7(b). Note, meshes at the boundaries are conformal. We first compute S11 parameter of the waveguide via VGFEM over a frequency range. Then, DD-VGFEM simulation is performed for the transmission factor of βtc = 105 . Legendre polynomials with p = 1 order of local approximation is used in both simulations. VGFEM and DD-VGFEM results depicted in Figure 4.7(b) match, and they are in good 104 (a) Geometry of WR-90 S-bend waveguide −10 HFSS VGFEM DD−VGFEM −20 S11 [dB] −30 −40 −50 −60 −70 8 8.5 9 9.5 10 Freq [GHz ] 10.5 11 11.5 (b) S11 of WR-90 S-bend waveguide Figure 4.7: VGFEM with DDM transmission conditions applied to the wave propagation in WR-90 S-bend waveguide. MoM result is from [1] 105 agreement with MoM result given in [1]. Finally, the method is applied for scattering from a dielectric object in a rectangular waveguide depicted in Figure 4.8(a). The waveguide dimensions are a = 2b, d = 0.8b, c = 3d, and b = 10mm. Relative dielectric constant is ǫr = 6. The problem domain is partitioned into two homogeneous sub-domains: free space and dielectric regions. Entire domain is meshed using brick elements with an average edge length of h = 3.3mm. Figure 4.8(b) shows computed reflection coefficients for different p values. It clearly shows the p−convergence of DD-VGFEM for this mesh. VGFEM result with p = 2 agrees with the reference data given in [13]. Hereafter, DG-VGFEM results are presented. First, DG-VGFEM is applied to wave scattering in WR-90 S-bend waveguide. Figure 4.9(a) shows the geometry of the problem, where a = 22.9mm, b = 10.2mm, c1 = 2a, L = 25mm, r = 15.24mm, and α = 30◦ . The geometry is meshed using tetrahedral elements with an average edge length of h = 5.5mm. The waveguide is partitioned into 5 sub-domains with depicted DG transmission boundaries in Fig. 4.9(a). Note, meshes at the boundaries are conformal. We first compute S11 parameter of the waveguide via VGFEM over a frequency range. Then, DG-VGFEM simulation is performed for the DG penalty factor of βdg = 105 . Legendre polynomials with p = 1 order of local approximation is used in both simulations. VGFEM and DG-VGFEM results depicted in Figure 4.9(b) match, and they are in good agreement with MoM result given in [1]. We further demonstrate the application of DG-VGFEM simulating the scattering in a H -plane WR-75 waveguide. The waveguide dimensions are a = 19.05mm, b = 9.525mm, c1 = 2a, r = 10.7mm, and the bend angle of α = 180◦ . The top view of the waveguide 106 (a) Geometry 1 DD−VGFEM, p=1 DD−VGFEM, p=2 Katzier 0.9 0.8 0.7 Γ 0.6 0.5 0.4 0.3 0.2 0.1 0 8 8.5 9 9.5 10 10.5 11 Frequency [GHz] 11.5 12 12.5 (b) Reflection Coefficient Figure 4.8: Application of DD-VFEM to scattering from a dielectric object in a rectangular waveguide.The waveguide dimensions are a = 2b, d = 0.8b, c = 3d, and b = 10mm, and relative dielectric constant is ǫr = 6 107 (a) Geometry of WR-90 S-bend waveguide 0 VGFEM DG−VGFEM MoM* −10 −20 S11 [dB] −30 −40 −50 −60 −70 −80 8 8.5 9 9.5 10 Frequency [GHz ] 10.5 11 11.5 (b) S11 of WR-90 S-bend waveguide Figure 4.9: VGFEM with DG transmission conditions applied to the wave propagation in WR-90 S-bend waveguide. Legendre polynomials with p = 1 local approximation is used. MoM result is from [1] 108 Figure 4.10: VGFEM with DG transmission conditions applied to the wave scattering in WR-75 waveguide with 180◦ bend. MoM result is from [6] geometry is shown in Fig. 4.10. The waveguide is meshed using tetrahedrons with the average edge length of h = 3.18mm. The order of local approximation is p = 2. Figure 4.10 compares the VGFEM result with the HFSS result and MoM data. HFSS and VGFEM results match, and they are in agreement with the MoM data given in [6]. Next, the accuracy of DG-VGFEM is further investigated by simulating scattering parameters of an H-plane T -junction. DG interfaces shown in Fig. 4.11(a) are chosen to be close to the junctions in order to show the capability of DG-VGFEM in capturing the higher order fields at the interfaces. The waveguide with the dimensions of a = 15.799mm, b = 7.899mm, and c1 = 2a as depicted in Fig. 4.11(a) is meshed using 1588 nodes and 6477 tetrahedrons with h = 2.78mm. Local approximation order of p = 1 is used. Figure 4.11(b) compares absolute values of the scattering parameters at the ports against the reference ones 109 (a) Geometry 15 S11 Reference S21 Reference S31 Reference 10 1 / | Si j | S11 DG−VGFEM S21 DG−VGFEM S31 DG−VGFEM 5 0 12 13 14 15 16 Frequency [GHz ] 17 18 (b) Scattering parameters Figure 4.11: Simulated scattering parameters of H-plane T-junction are compared against the reference data in [7]. 110 0 MoM DG−VGFEM, Plane waves DG−VGFEM, Mixed basis −10 −20 S11 [dB] −30 −40 −50 −60 −70 −80 8 8.5 9 9.5 10 Frequency [GHz ] 10.5 11 11.5 Figure 4.12: Use of mixed basis functions within a simulation in VGFEM with DG transmission conditions. Problem is scattering in a WR-90 S-bend waveguide. MoM result is from [1] given in [7]. Simulated parameters are in excellent agreement with the reference data. In the next simulation, we demonstrate use of mixed basis functions within a DG-VGFEM simulation. Again, wave propagation in WR-90 waveguide is considered for this demonstration. Plane waves and Legendre polynomials are used. The space of plane wave basis functions is span{e−jkp ·(r−ri ) }, where ri is the position of the node i, and kp , kp = (k0 , θ, φ), is the direction of the plane wave basis functions. Four particular plane wave directions (θ, φ) ∈ {(π/4, 0), (π/4, π), (π/2, 0), (π/2, π)} are chosen. While plane wave basis functions are used in the straight sections, Legendre polynomials with p = 2 order of local approximation are exploited in the bent sections (see the geometry in Fig. 4.9(a)). Figure 4.12 compares this simulation against the DG-VGFEM simulation with plane wave basis used in entire region. Both results are in agreement with MoM result given in [1]. Note, the 111 numbers of basis functions in each PU domain Ωi are Nbasis = 8 for the plane wave basis, and Nbasis = 26 for the Legendre polynomials with p = 2 order of local approximation. Finally, we show an application of the DG-VGFEM for inhomogeneous problems. Scattering from a dielectric object in a rectangular waveguide depicted in Figure 4.8(a) is considered for this demonstration. The problem domain is partitioned into two homogeneous sub-domains: free space and dielectric regions. Entire domain is meshed using brick elements with an average edge length of h = 2.2mm. Relative dielectric constant of ǫr = 6 and b = 10mm are used. Figure 4.13 shows computed reflection coefficients for different p values. It clearly shows the p−convergence of DG-VGFEM for this mesh. DG-VGFEM result with p = 2 agrees with the reference data given in [13]. 1 DG−VGFEM, p=1 DG−VGFEM, p=2 Reference [Katzier] 0.9 0.8 0.7 Γ 0.6 0.5 0.4 0.3 0.2 0.1 0 8 8.5 9 9.5 10 10.5 11 Frequency [GHz] 11.5 12 12.5 Figure 4.13: VGFEM with DG transmission conditions applied to wave scattering problem in a rectangular waveguide with a dielectric obstacle in it. Relative dielectric constant is ǫr = 6. 112 VGFEM Region Γ FEM Region Γ dg PEC Ω2 Ω1 Figure 4.14: Illustrations of domain decomposition for FE-VGFEM 4.2 Hybridization of VGFEM with FEM In this Section, we hybridize FEM and VGFEM using the decomposition techniques presented above. With this hybridization, it is aimed to use the advantageous of each method within the simulations such that more accurate and perhaps efficient simulations are obtained. The formulation of FE-VGFEM is presented next. For simplicity and without loss of generality, consider a partition of Ω into two sub-domains Ω1 and Ω2 with a common interface Γdg as shown in Fig. 4.14. Let VGFEM be used in Ω1 and FEM be used in Ω2 , and let u1 and u2 denote the electric fields in each sub-domain, respectively. A bilinear formulation for each sub-domain can be obtained by applying Galerkin’s method to (4.1). The bilinear formulation for VGFEM has been presented above for each methodology. The bilinear formulation of FEM is well presented in [13]; therefore we will not repeat the well known FEM formulation here. The linear edge elements, 0th order vector basis, are used in FEM system in this thesis. These bilinear formulations of FEM and VGFEM standalone are not enough to solve the fields in entire domain since the FEM and VGFEM basis functions 113 do not interact with each other. Therefore, an additional boundary condition is needed to enable the coupling of fields at the sub-domain interfaces. Here, we will use either the robin transmission conditions or the interior penalty discontinuous Galerkin conditions for coupling. For the coupling with the IPDG, we will only enforce the tangential continuity of the fields given by Γdg 4.2.1 1 βdg [[w]]T · [[u]]T + [[w]]T · { ∇ × u} dΓ = 0. µr (4.13) Numerical Results In this subsection, the hybridization of classical FEM with VGFEM is validated via the analysis of both open and closed domain problems. The advantages of both methods are utilized for efficient simulations. An FEM code with linear edge elements has been written for this hybridization. For open domain problems, the total field formulation with the first order ABC is used, and the radar cross sections (RCS) are computed by using the simulated field values in the formulations given in [13]. First, we will validate the method analyzing the T E10 mode propagation in a 1.0 m × 0.5 m × 0.75 m hollow waveguide. FEM is used in z = [0, 0.4m] and VGFEM is used in z = [0.4m, 0.75m] for this study. Robin transmission condition is used at the interface. The domain is meshed using tetrahedrons with an average edge length of h = λo /10, where λo = 1.8m. The order of local approximation is p = 0. Figure 4.15 shows the electric field in the waveguide. While VGFEM approximate the field accurately thanks to its continuous framework, FEM introduces discontinuities at the tetrahedron interfaces. Next, we apply the hybrid method to plane wave propagation in inhomogeneous domains. 114 Re { Ey } 0 −0.5 1 −1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 z 0.2 0 x 0 (a) Real part of Ey Im { Ey } 0 −0.5 1 −1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 z 0 x (b) Imaginary part of Ey Figure 4.15: T E10 mode propagation in a 1.0 m × 0.5 m × 0.75 m waveguide is analyzed using FE-VGFEM. FEM is used in z = [0, 0.4m] and VGFEM is used in z = [0.4m, 0.75m] for this study 115 Re { Ez } 1 0 −1 0 0.8 0.6 0.5 0.4 x (m) z (m) 0.2 1 0 (a) Real part of Ez z Im { E } 1 0 0 0.2 −1 0 0.4 0.2 0.6 0.4 z (m) x (m) 0.8 0.6 0.8 1 (b) Imaginary part of Ez Figure 4.16: Plane wave propagating in a 1.0 m × 0.5 m × 0.75 m domain partially filled with ǫr = 4 . VGFEM is used in dielectric region z = [0, 0.4m] and FEM is used in z = [0.4m, 0.75m] 116 The parallel polarized plane wave is propagating with the direction of θinc = π/4 and φinc = 0 in a problem domains of 1.0 m × 0.5 m × 0.75 m. VGFEM is used in dielectric region, z = [0, 0.4m], and FEM is used in free space region,z = [0.4m, 0.75m]. Dielectric constant is ǫr = 4. Robin transmission condition is used at the interface. The domain is meshed using tetrahedrons with an average edge length of h = λo /10, where λo = 1.8m. The order of local approximation is p = 0. Figure 4.16 shows the normal component of the electric field. While VGFEM approximate the field accurately thanks to its framework providing continuity across PU domains, FEM introduces discontinuities at the tetrahedron interfaces. Next, FE-VGFEM is applied to wave scattering at the junction of a rectangular waveguide and a circular cavity. The rectangular waveguide is WR-75 with dimensions a = 19.05 mm, b = 9.525 mm, and c = 2a. The radius and depth of the circular waveguide are r = 13.589 mm and h = 40.945 mm respectively. The waveguide is partitioned into 5 subdomains as shown in Fig. 4.17(a), where FEM and VGFEM regions are depicted. As FEM basis functions can easily handle the boundary conditions, they are used around the junctions. The thickness of FEM regions is t = 3.81 mm. The number of unknowns in FEM and VGFEM regions are Nf em = 16700 and Nvgf em = 29592 respectively, with a total number of unknowns N = 46292. Simulated scattering parameters agree with the measured data given in [8] as shown in Fig. 4.17(b). If only FEM is used in entire region, much finer mesh is required for the same level of accuracy. Figure 4.17(b) also shows the results for such FEM simulation with N = 68735 unknowns. Next, FE-VGFEM is applied to wave scattering in a rectangular waveguide with a dielectric object with ǫr = 6. The problem domain shown in 4.8(a) is partitioned into two 117 n1  dg FEM Region VGFEM Region FEM Region Γ VGFEM Region VGFEM Region n2 ̂ t c (a) Geometry 0 −5 −10 S−parameters −15 S11 measured −20 S21 measured −25 S FEM −30 S FEM −35 S FE−VGFEM S FE−VGFEM 11 21 11 21 −40 −45 −50 8 9 10 11 12 Frequency (GHz) 13 14 15 (b) S paramater Figure 4.17: Scattering parameters for WR-75 -circular waveguide junction. Measured results are from [8]. 118 1 FE−VGFEM Reference [Katzier] Reflection Coefficient 0.8 0.6 0.4 0.2 0 8 8.5 9 9.5 10 10.5 11 Frequency [GHz] 11.5 12 12.5 Figure 4.18: Reflection coefficient for a rectangular waveguide with a dielectric obstacle in it. Relative dielectric constant is ǫr = 6. sub-domains, which are separated by a DG boundary placed inside the dielectric region and 1 mm away from the dielectric-air interface. While VGFEM with p = 0 is used in the inner region, FEM is utilized in the outer region such that Nf em = 11680, Nvgf em = 2169, and N = 13849. Simulated reflection coefficient matches the reference data given in [13] as shown in Fig. 4.18. Next, FE-VGFEM is applied to wave scattering from a dielectric sphere with ǫr = 4 and ka = 1, where a is the radius of the sphere. The sphere is centred at the origin, and the infinite domain is truncated by a spherical boundary with a radius of rabc = 1.5a. The sphere is illuminated by the electric field E = xe−jk0 z . Since VGFEM can approximate ˆ the fields more accurately, it is used in dielectric region while FEM is used in the free space as depicted in Fig. 4.19(a). The average mesh size in both regions is h = λ/12. While p = 1 119 FEM Region Γ dg VGFEM Region ǫr =4 Ω1 Ω2 (a) Geometry 0 Mie series FE−VGFEM −5 −10 0 σ / λ2 [dB] −15 −20 −25 −30 −35 −40 −45 −50 0 20 40 60 80 100 θ [degrees] 120 140 160 180 (b) RCS Figure 4.19: Bistatic RCS of a dielectric sphere with ǫr = 4 and ka = 1. FEM and VGFEM regions are depicted. The sphere is illuminated by the electric field E = xe−jk0 z ˆ 120 is used in the VGFEM region, p = 0 is used in the FEM region. The simulated bistatic RCS agrees well the analytic result as shown in Figure 4.19(b). Next, an advantage of FEM within FE-VGFEM framework is demonstrated by simulating the bistatic RCS of a 0.3m × 0.3m PEC plate for the incident electric field E = xejk0 z at ˆ 300MHz. The plate is centred at the origin on xy plane, and the infinite domain is truncated by a spherical boundary with radius of rabc = 0.45λ. Since imposing the boundary condition on PEC is straightforward for FEM, and results in reduction of number of unknowns, a very thin layer of FEM region with the thickness of t = λ/20 is chosen around the plate as shown in Fig. 4.20(a). For VGFEM, while p = 1 is used in the PU domains associated with inner nodes, p = 2 is utilized in the PU domains associated with the nodes on the ABC surface. The average mesh size in VGFEM region is h = λ/5. Figure 4.20(b) shows excellent agreement between FE-VGFEM and VGFEM, and they match MoM data obtained using our in-house code. Next, another advantage of FEM within FE-VGFEM framework is demonstrated by simulating the bistatic RCS of a PEC patch placed on a PEC backed dielectric substrate with ǫr = 2.17 for the incident electric field E = xejk0 z at 300MHz. Patch size is 0.27m × ˆ 0.19m and substrate size is 0.44m × 0.44m × 0.05m. The infinite domain is truncated by a 0.88m × 0.88m × 0.49m rectangular box as depicted in Fig. 4.21(a). Since imposing the boundary condition on PEC is straightforward for FEM, and FEM basis space satisfies the boundary conditions at the material interface, FEM is used in a 0.64m × 0.64m × 0.25m region that confines the geometry as shown in in Fig. 4.21(a). VGFEM is used in the rest of the domain as shown in the same Figure. While p = 1 and h = 0.1m are used in VGFEM region, p = 0 and h = 0.05m are used in FEM region. Figure 4.21(b) shows the simulated 121 Γ dg FEM Region PEC Ω1 VGFEM Region (a) Geometry 0 MoM VGFEM FE−VGFEM −5 2 σ / λ [dB] −10 −15 −20 −25 −30 −35 −40 0 10 20 30 40 50 θ [degrees] 60 70 80 90 (b) RCS Figure 4.20: Bistatic RCS of a 0.3m × 0.3m PEC plate. FEM and VGFEM regions are depicted. The plate is illuminated by the electric field E = xejk0 z ˆ 122 FEM Region VGFEM Region Γ dg (a) Geometry 10 HFSS, p=1 HFSS, p=0 FE−VGFEM 0 0 σ / λ2 [dB] 5 −5 −10 −15 0 10 20 30 40 50 θ [degrees] 60 70 80 90 (b) RCS Figure 4.21: Bistatic RCS of a PEC patch placed on a PEC backed dielectric substrate with ǫr = 2.17 at f0 = 300MHz. Patch size is 0.27m × 0.19m and substrate size is 0.44m × 0.44m × 0.05m. The geometry is illuminated by the electric field E = xejk0 z ˆ 123 bistatic RCS against the HFSS results. For HFSS simulations, adaptive option is used such that total number of unknowns are N ≈ 106.000 for p = 0 with a final pass number 16, and N ≈ 94.000 for p = 1 with a final pass number 13. On the other hand, total number of unknowns for FE-VGFEM is N = 60343, where the FEM and VGFEM unknowns are Nf em = 36835 and Nvgf em = 23508 respectively. Although FE-VGFEM uses the first order ABC and total field formulation, it provides accurate solutions with fewer unknowns. Next, FE-VGFEM with a spherical DG interface is tested by simulating the bistatic RCS of a sphere with radius of r = 0.5λ. The sphere centred at the origin is illuminated by the electric field E = xejk0 z . A spherical truncation boundary is placed at 0.5λ away from the ˆ sphere, and the first order ABC is imposed at the boundary. A thin spherical shell with the thickness of t = λ/20 around the sphere is considered for FEM region as shown in Fig. 4.22(a). In VGFEM region, while p = 1 is used in the PU domains associated with inner nodes, p = 2 is utilized in the PU domains associated with the nodes on the ABC surface. Figure 4.22(b) shows the comparison of the RCS plots obtained via FE-VGFEM, HFSS, and Mie series. Note that HFSS uses the scattering field formulation with 2nd order ABC, and the simulation is performed using the adaptive option with initial parameters of h = λ/10 and p = 1, which results in a total matrix size of N = 33698. On the other hand, total number of unknowns for FE-VGFEM is N = 18984, where the number of FEM unknowns is only Nf em = 1932. This simulation validates that FE-VGFEM can simulate the problem efficiently and accurately as the advantages of both methods are exploited in the simulation. Next, the capability of FE-VGFEM with a cubical DG interface is tested by simulating the bistatic RCS of a cube with the edge length of a = 0.755λ for λ = 1.0m. The cube centred at the origin is illuminated by the electric field E = y e−jk0 z . A spherical boundary ˆ 124 Γ dg FEM Region PEC Ω1 VGFEM Region (a) Geometry 15 Mie HFSS FE−VGFEM 10 σ [dBsm] 5 0 −5 −10 −15 0 20 40 60 80 100 θ [degrees] 120 140 160 180 (b) RCS Figure 4.22: Bistatic RCS of a PEC sphere with radius of r = 0.5λ. FEM and VGFEM regions are depicted. The sphere is illuminated by the electric field E = xejk0 z . ˆ 125 Γ FEM Region dg PEC Ω1 VGFEM Region (a) Geometry 15 HFSS FE−VGFEM Measured 10 2 σ / λ [dB] 5 0 −5 −10 −15 −20 −25 0 20 40 60 80 100 θ [degrees] 120 140 160 180 (b) RCS Figure 4.23: Bistatic RCS of a PEC cube with the edge length of a = 0.755λ. FEM and VGFEM regions are depicted. The sphere is illuminated by the electric field E = y e−jk0 z . ˆ 126 with the radius of rabc = 1.0λ is used for domain truncation. A thin cubical shell with a thickness of t = λ/20 around the cube is chosen for FEM region as shown in Fig. 4.23(a). In VGFEM region, while p = 1 is used in the PU domains associated with inner nodes, p = 2 is utilized in the PU domains associated with the nodes on the ABC surface. Note that HFSS uses scattering field formulation with the 2nd order ABC, and the simulation is performed using the adaptive option with initial parameters of h = λ/10 and p = 1, which results in a total matrix size of N = 49670. On the other hand, total number of unknowns for FEVGFEM is N = 30383, where the number of FEM unknowns is only Nf em = 3647. Figure 4.23(b) compares the simulated RCS against the HFSS result, and measured data [37], and they are in good agreement. Finally, the features of FE-VGFEM are further demonstrated by simulating radiation from a horn antenna. The horn is fed by the T E10 mode of a rectangular waveguide with dimensions 2.29cm × 1.02cm × 4.58cm as shown in Fig. 4.24(a). The aperture size of the horn is 6.75cm × 4.95cm, and the horn length from the waveguide junction to the aperture is 13.87cm. The infinite domain is truncated by a 10.125cm × 7.425cm × 21.825cm ABC box as shown in Fig. 4.24(a). While VGFEM is used inside the horn and waveguide, FEM is used in the rest of the domain. Two FE-VGFEM simulations have been performed to compute the gain of the antenna at 10GHz. In the first simulation, Legendre polynomials with p = 0 are used in the VGFEM regions. In the second simulation, while Legendre polynomials with p = 0 is used in the horn, T E10 mode is used in the waveguide in order to show accurate simulations with the inclusion of physics in approximation space.Figure 4.24(b) compares the simulation results against the MoM data obtained in [38] and they are in good agreement. The simulation with modal basis is much closer to the MoM data as expected. 127 VGFEM EM VGF Region on Regi FEM Region Γ dg (a) Geometry 25 MoM FE−VGFEM, p=0 FE−VGFEM, p=0 and modal basis 20 15 Gain [dB] 10 5 0 −5 −10 −15 −20 −100 −50 0 θ [degrees] 50 100 (b) RCS Figure 4.24: Bistatic RCS of a PEC cube with the edge length of a = 0.755λ. FEM and VGFEM regions are depicted. The sphere is illuminated by the electric field E = y e−jk0 z . ˆ 128 4.3 Summary In this Chapter, we have hybridized VGFEM with other finite element based methods, namely classical FEM, DDM and DGM. DDM and DGM transmission conditions are used to couple the methods with VGFEM. This development makes the method applicable to inhomogeneous problems with its present continuous framework. Volumetric basis functions are forced to satisfy the transmission conditions. This approach can be extended nonconforming boundaries easily in the VGFEM framework. We have demonstrated the applications of the hybrid methods analyzing some practical problems. h− and p− convergence of the method has been shown using Legendre polynomials. In addition, we have demonstrated some features of VGFEM such as mixed basis functions and mixed polynomial orders in some DG-VGFEM simulation. This hybridization enables us to exploit advantages of each method in solving complex electromagnetic problems accurately, and perhaps, efficiently. Moreover, this decomposition technique has been successfully applied to hybridize FEM with VGFEM. FE-VGFEM has been validated through the analysis of various open and closed domain problems. Further development of these hybrid methods to attain a highly flexible and powerful EM solver is a topic of our current research. 129 Chapter 5 VGFEM for Time Domain Electromagnetic Analysis In this Chapter, a time domain VGFEM (TD-VGFEM) is developed to utilize VGFEM for solving time domain EM problems. Indeed, there was an earlier attempt to develop TDVGFEM [28], where the meshless form of the method was exploited along with different time marching schemes. In this dissertation, TD-VGFEM is developed based on the polyhedral meshes, and the Newmark time stepping scheme is used during this development. Moreover, the method is hybridized with the classical finite element methods and discontinuous Galerkin method (DGM). Through these hybridizations, more accurate and efficient transient simulations are aimed by using the advantageous of each method ins simulations. These hybridizations also make VGFEM applicable to inhomogeneous problems. Accuracy, efficiency, and stability of the hybrid methods are examined in detail. The Chapter is organized as follows: First, the bilinear formulation of TD-VGFEM for different boundary conditions is given in Section 5.1. Formulations for hybrid VGFEM 130 methods are also presented here. In Section 5.2, the application of TD-VGFEM to wave propagations using Dirichlet and Neumann boundary conditions are presented. Temporal convergence, and h−and −pconvergence of the method are shown. Finally, conclusions and future research are presented in Section 5.3. 5.1 TD-VGFEM Formulation This section introduces the framework of VGFEM built on polyhedral meshes briefly, and then presents the formulation of TD-VGFEM for different boundary conditions in details. Consider a linear, homogeneous, and source free three-dimensional region Ω, whose boundary i Γi . In Ω, the time-harmonic electric field, u(r, t), satisfies the is denoted by ∂Ω := Γ = vector wave equation ∇× 1 1 d2 ǫr u(r, t) = 0 ∇ × u(r, t) + µr c2 dt2 0 Bi {u(r, t)} = gi (r, t), r ∀ ∈ Γi , (5.1) where c0 is the speed of light in free space, ǫr and µr are relative constitutive parameters, Bi is a differential operator of boundary conditions, and gi (r, t) is a function of boundary conditions imposed on Γi . Solution to (5.1) using TD-VGFEM is achieved as follow: First, the problem domain is meshed and spatial basis functions are defined on this tessellation. Then, time axis is discretized and appropriate temporal basis functions are defined. Next, a system of equations is set up via Galerkin testing procedure in both time and space. Finally, a time marching scheme with initial conditions is established to solve unknown coefficients at each time step. 131 We will next present this solution procedure in details. A bilinear formulation of TD-VGFEM can be derived using following boundary conditions • u(r) × n = g1 (r), r ∈ Γd , ˆ 1 • n × µr ∇ × u(r) = g2 (r), r ∈ Γn , ˆ 1 • n × µr ∇ × u(r) + c1 n × n × u(r) = g3 (r), r ∈ Γimp ˆ ˆ ˆ 0 • n × ∇ × u(r) + c1 n × n × u(r) = 0, f or r ∈ Γabc , ˆ ˆ ˆ 0 where n is outward unit normal vector, Γd , Γn , Γimp , and Γabc denote Dirichlet, Neumann, ˆ impedance, and absorbing boundaries respectively, and g1 (r), g2 (r), and g3 (r) are the boundary functions for the corresponding boundary conditions. Using the spatial vector space defined in Chapter 2, and applying Galerkin’s method to (5.1), one can obtain a second order differential equation M d d2 a(t) + C a(t) + Ka(t) = f, 2 dt dt (5.2) where a(t) is the time function of unknown coefficients, M, C, and K are sparse matrices, and f is a sparse vector. Entries of these spatial matrices and entry of the vector f are respectively given by 1 ǫr wm (r) · ul (r) dΩm Mml = c2 Ωm 0 132 (5.3) Cml = 1 (wm (r) × n) · (ul (r) × n) dΓm + ˆ ˆ c0 Γimp m 1 (wm (r) × n) · (ul (r) × n) dΓm ˆ ˆ c0 Γabc m Kml = Ωm +β1 Γd m (wm (r) × n) · ( ˆ (∇ × wm (r)) · ( Ωm 1 ∇ × ul (r)) dΩm µr (∇ · ǫr wm (r))(∇ · ǫr ul (r)) dΩm 1 1 ∇ × ul (r))dΓ + ( ∇ × wm (r)) · (ul (r) × n)dΓ+ ˆ µr Γd µ r m +β2 fm = (5.4) Γd m (wm (r) × n) · (ul (r) × n) dΓm ˆ ˆ (5.5) 1 (wm (r) × n) · g1 (r) dΓm − ˆ ∇ × wm (r)) · g1 (r) dΓm + β2 Γd Γd µ r m m ( Γn m wm (r) · g2 (r) dΓm − imp wm (r) · g3 (r) dΓm Γm (5.6) where m is the index of testing function wm (r), l is the index of basis function ul (r), and β1 and β2 are the constants used to get uniform convergence [25]. In the previous formulations, it has been assumed that there is no material variation throughout the domain. If there is a material variation, then the domain decomposition methodology presented in Chapter 4 is applied. As the interior penalty DGM is used, the corresponding boundary terms are added to the matrix K as follow K=K+ Γdg 1 βdg [[w]]T · [[u]]T + [[w]]T · { ∇ × u} dΓ µr 133 (5.7) Likewise, as VGFEM and FEM are hybridized, FEM and VGFEM domains are coupled using the DG boundary condition in (5.7). Note, the bilinear formulation for FEM can be found in [13]. Equation (5.2) can be solved numerically using the Newmark method [14]. To this end, time axis is uniformly discretized using Nt number of time samples and ∆t time step, and quadratic finite element temporal basis functions are defined on this mesh. Then, the Newmark equation is written as [14] 1 ∆t2 − − 1 2M + M+ γ C + βK an+1 = ∆t 1 (1 − 2γ)C + (0.5 − 2β + γ)K an ∆t ∆t2 1 1 M− (1 − γ)C + (0.5 + β − γ)K an−1 − 2 ∆t ∆t + βf n+1 + (0.5 + γ − 2β)f n + (0.5 − γ + β)f n−1 where a is the vector of unknown coefficients, n is the time step index, and β and γ are the parameters that control accuracy and stability of the time scheme. If β ≥ 0.25 and γ = 0.5, then the system is unconditionally stable [14]. 5.2 Numerical Results In this Section, we validate the time domain VGFEM, DG-VGFEM, and FE-VGFEM formulations, and present their convergence characteristics analyzing wave propagation problems. For all simulations, a modulated Gaussian pulse (t−tp )2 − 2σ 2 T (t) = cos(2πf0 t)e 134 (5.8) 1 Analytic Neumann BC, p=0 Neumann BC, p=1 0 z E (V/m) 0.5 −0.5 −1 0 2 4 time (seconds) 6 8 −8 x 10 Figure 5.1: Electric field E = z T (t+ y ·r/c0 ) is observed at x = 0.5m, y = 0.25m, z = 0.375m ˆ ˆ in 1.0m × 0.5m × 0.75m domain. TD-VGFEM with Neumann BC, p = 0 and h = 0.225m is used. with f0 = 200MHz, BW = 100MHz, σ = 3/(2πBW ) and tp = 6σ is used. Unless otherwise stated, time step of ∆t = 1.7 ×10−10 s is chosen in Newmark method with γ = 0.5 and β = 0.4. Tetrahedral meshes and Legendre polynomials are used. Note, one needs to use Legendre polynomial orders of p + 1 and p + 2 to get pth order local approximation. First, we investigate approximations of wave propagation in free space via TD-VGFEM for different boundary conditions. For this study, we consider the electric field E = z T (t + ˆ y · r/c0 ) propagating through 1.0m × 0.5m × 0.75m domain. The domain is discretized using ˆ an average mesh size of h = 0.225m. On the domain boundaries, any of Neumann BC, Dirichlet BC and impedance BC is imposed, and transient behaviour of TD-VGFEM are analyzed for each of them. Figure 5.1 shows analytic and computed electric fields observed at x = 0.5m, y = 0.25m, z = 0.375m. Neumann BC is not stable for p = 0 since the derivative of the local approximation function vanishes and only PU function construct the basis vector. However, Neumann BC approximates the field well enough for p = 1. Figure 135 1 Analytic Dirichlet BC, p=0 0 z E [V/m] 0.5 −0.5 −1 0 2 4 time [seconds] 6 8 −8 x 10 Figure 5.2: Electric field E = z T (t+ y ·r/c0 ) is observed at x = 0.5m, y = 0.25m, z = 0.375m ˆ ˆ in 1.0m × 0.5m × 0.75m domain. TD-VGFEM with Dirichlet BC, p = 0 and h = 0.225m is used. 5.2 shows the electric field computed applying Dirichlet BC for p = 0. TD-VGFEM result matches the analytic result. Next, we study temporal, spatial (h−) and p− convergence for impedance boundary condition. In addition, we investigate the convergences for the use of mixed polynomial basis orders within a domain. For mixed polynomial order simulations, the vector basis functions with p = 1 order approximation is assigned to the PU domains for ri (2) < 0.25m, where ri = (xi , yi, zi ) is the position of the node i. First, we analyze h− and p−convergence for ∆t = 1.7 × 10−10 s. As shown in Fig. 5.3(a), mixed order result is closer to p = 1 result for coarse discretization; however, this behaviour is reverse for fine meshing as expected. We also conclude that reducing the mesh size does not reduce the error after some point since the error in temporal approximation becomes dominant. Hence, we reduce the time step to ∆t = 0.43 × 10−10 and repeat the simulations shown in Fig. 5.3(b). Now, with a lower temporal error threshold, h−and p−convergence is obtained. 136 −1 10 p=0 p=1 p=0 and p=1 −2 −3 10 2 L Relative Error 10 −4 10 −5 10 1 10 2 3 10 10 Number of unknowns (a) ∆t = 1.7 × 10−10 4 10 −1 10 p=0 p=1 p=0 and p=1 −2 −3 10 −4 10 2 L Relative Error 10 −5 10 −6 10 1 10 2 3 10 10 Number of unknowns (b) ∆t = 0.43 × 10−10 4 10 Figure 5.3: Temporal convergence, and h− and p−convergence of TD-VGFEM for Impedance BC. Electric field E = z T (t + y · r/c0 ) propagating in 1.0m × 0.5m × 0.75m ˆ ˆ domain is considered for this study. 137 Hereafter, we study the convergence characteristics of time domains DG-VGFEM and FE-VGFEM. A computation domain of [0, 1] × [0, 0.5] × [0, 1] in meters is divided into two regions with an DG interface at z = 0.5. For FE-VGFEM convergence study, while FEM basis functions are used forz < 0.5, VGFEM basis functions are used for z > 0.5. A normally ˆ ˆ ˆ ˆ incidence time signal E = θT (t− k·r/c) with k = r with θinc = 0 and φinc = 0 is considered. The impedance BC is imposed at the domain boundaries, and time step of ∆t = 4.17×10−11 is used. Figure 5.4 shows h−, and p−convergence plots of DG-VGFEM and FE-VGFEM. The temporal error becomes dominant after certain mesh size as seen in the convergence plots in DG-VGFEM in Fig 5.4(a). Next, the convergence characteristics of both methods are examined for oblique incidence. ˆ ˆ ˆ ˆ The electric fieldE = θT (t − k · r/c) with k = r , θinc = π/4 and φinc = 0 is propagating in the domain presented in the previous example. The domain partitioning is also as before. The impedance BC is imposed at the domain boundaries, and time step of ∆t = 4.17×10−11 is used. Figure 5.5 shows h−, and p−convergence plots of DG-VGFEM and FE-VGFEM. Next example will examine the use of the DG conditions at the material interfaces. The ˆ ˆ ˆ electric fieldE = θT (t − k · r/c) with k = r, θinc = π/4 and φinc = 0 is propagating in a ˆ domain [0, 1] × [0, 0.5] × [0, 1] in meters. A dielectric material with ǫr = 4 is used in the domain for z > 0.5m, while it is free space for z < 0.5m. The Dirichlet BC is imposed at the domain boundaries, and time step of ∆t = 4.17 × 10−11 is used. Figure 5.6 shows h−, and p−convergence plots of DG-VGFEM and FE-VGFEM. Next, the VGFEM formulation with the first order ABC is tested via simulation of radiation from a dipole centred at the origin. Dipole current is J = z d T (t)δ(r), where T (t) ˆ dt is the modulated Gaussian pulse with ∆t = 1.7 × 10−10 . The problem domain is a spherical 138 −1 10 VGFEM, p=0 VGFEM, p=1 VGFEM, p=2 DG−VGFEM, p=0 DG−VGFEM, p=1 DG−VGFEM, p=2 −2 2 L Relative Error 10 −3 10 −4 10 −5 10 −6 10 −7 10 2 3 10 4 10 10 Number of spatial unknowns 5 10 (a) DG-VGFEM −1 L2 Relative Error 10 FEM FE−VGFEM VGFEM −2 10 −3 10 −4 10 2 10 3 10 Number of spatial unknowns 4 10 (b) FE-VGFEM Figure 5.4: h−, and p− convergence of DG-VGFEM and FE-VGFEM for the normally ˆ ˆ ˆ incident electric field E = θT (t − k · r/c) with k = r , θinc = 0 and φinc = 0 in 1.0m × ˆ 0.5m × 1.0m. 139 −1 10 VGFEM, p=0 VGFEM, p=1 VGFEM, p=2 DG−VGFEM, p=0 DG−VGFEM, p=1 DG−VGFEM, p=2 −2 L2 Relative Error 10 −3 10 −4 10 −5 10 −6 10 −7 10 2 3 10 4 5 10 10 Number of spatial unknowns 10 (a) DG-VGFEM −1 2 L Relative Error 10 VGFEM FEM FE−VGFEM −2 10 −3 10 −4 10 2 10 3 10 Number of spatial unknowns 4 10 (b) FE-VGFEM Figure 5.5: h−, and p− convergence of DG-VGFEM and FE-VGFEM for the obliquely ˆ ˆ ˆ incident electric field E = θT (t − k · r/c) with k = r, θinc = π/4 and φinc = 0 in 1.0m × ˆ 0.5m × 1.0m. 140 0 10 DG−VGFEM, p=0 DG−VGFEM, p=1 DG−VGFEM, p=2 −1 L2 Relative Error 10 −2 10 −3 10 −4 10 −5 10 −6 10 2 10 3 4 10 10 Number of spatial unknowns 5 10 (a) DG-VGFEM 0 10 L2 Relative Error FEM FE−VGFEM DG−VGFEM −1 10 −2 10 −3 10 2 10 3 10 Number of spatial unknowns 4 10 (b) FE-VGFEM Figure 5.6: h−, and p− convergence of DG-VGFEM and FE-VGFEM for a problem with material variation. ǫr = 4 is used in the domain for z > 0.5m of the domain 1.0m × 0.5m × ˆ ˆ ˆ ˆ 1.0m. The electric field E = θT (t − k · r/c) with k = r , θinc = π/4 and φinc = 0 is obliquely incident. 141 shell with rin = 0.5m and rout = 0.75m. Dirichlet BC at the inner surface and 1st order ABC at the outer surface are imposed. Legendre polynomials with p = 2 and the average tetrahedral edge length h = 0.3m are utilized. Figure 5.7 compares the computed electric field against the analytic one observed at r = 0.667m, θ = π/4, and φ = 0. Next , FE-VGFEM is applied to the problem of scattering from a dielectric sphere with ǫr = 4 and a = 1/(2π), where a is the radius of the sphere. Incident field of Einc = xT (t + z · r/c) with ∆t = 4.17 × 10−11 is considered for this simulation. The radius of ˆ ˆ the ABC surface is rabc = 1.5/(2π). VGFEM is used in dielectric region and FEM in the rest of the domain. Simulation parameters are p = 0 and h = 0.08λ for FEM, and p = 1 and h = 0.09λ for VGFEM. Figure 5.8 compares the computed RCS against the MoM data obtained in [39], frequency domain (FD) FE-VGFEM data, and time domain FEM data for 300MHz. Finally, the TD-VGFEM is further validated via the simulation of the bistatic RCS of 0.3m × 0.3m plate. Domain is truncated by a cubical boundary located 0.3m away from the edge of the plate. 1st order ABC is imposed on the truncated boundary, and the average mesh size of h = λ0 /5 is used. Incident field of Einc = xT (t+ z ·r/c) with ∆t = 4.17×10−11 ˆ ˆ is considered for this simulation. Figure 5.9 shows computed bistatic RCS at 300MHz for different p values for both time and frequency domain formulations. Since the mesh is not dense enough, VGFEM result with p = 1 is a little off from the MoM results obtained using our in-house MoM code. However, VGFEM result with p = 2 matches the MoM result. 142 −6 4 x 10 Analytic TD−VGFEM 3 2 Ex (V / m) 1 0 −1 −2 −3 −4 0 0.2 0.4 0.6 time (sec) 0.8 1 −7 x 10 (a) Ex −6 4 x 10 Analytic TD−VGFEM 3 2 Ez (V / m) 1 0 −1 −2 −3 −4 0 0.2 0.4 0.6 time (sec) 0.8 1 −7 x 10 (b) Ez Figure 5.7: Time signal computed at r = 0.667m, θ = π/4, and φ = 0 for radiation from a dipole. 143 FEM Region Γ dg VGFEM Region ǫr =4 Ω1 Ω2 (a) Geometry 0 MoM, [Umashankar et al. IEEE 1986 FEM, TD FE−VGFEM, FD FE−VGFEM, TD −5 −10 2 0 σ / λ [dB] −15 −20 −25 −30 −35 −40 −45 −50 0 20 40 60 80 100 θ [degrees] 120 140 160 180 (b) RCS at f0 = 300MHz Figure 5.8: Bistatic RCS of a dielectric sphere with ǫr = 4 and a = 1/(2π). FEM and VGFEM regions are depicted. 144 0 −5 −10 RCS [dBsm] −15 −20 MoM FD−VGFEM, p=1 TD−VGFEM, p=1 FD−VGFEM, p=2 TD−VGFEM, p=2 −25 −30 −35 −40 −45 −50 0 10 20 30 40 50 θ [degrees] 60 70 80 90 Figure 5.9: Scattering from a 0.3m × 0.3m plate at 300MHz for Einc = xT (t + z · r/c), ˆ ˆ inc = φinc = φ = 0◦ . θ 5.3 Summary In this Chapter, time domain VGFEM, DG-VGFEM, and FE-VGFEM have been developed for time domain electromagnetic problems. The Newmark method is used for temporal approximation. The vector wave equation is used for the hybridization of VGFEM with FEM and DGM. The methods have been validated through convergence analysis as well as applications. Temporal, h−convergence, and p−convergence of the methods have been shown for different boundary conditions. The first order ABC has been used for scattering problems. Formulations with higher order ABC and perfectly matched layers, and their applications are the subjects of the future research. Moreover, different time marching schemes, and their convergence and stability characteristics within the VGFEM framework are also our current research topic. 145 Chapter 6 Conclusions and Future Work 6.1 Conclusions This dissertation develops a polyhedral based VGFEM framework for the numerical analysis of electromagnetic problems. This framework is an adaptation of the meshless VGFEM framework onto the finite element meshes. This advancement also provides the framework for integration of VGFEM with other FEM techniques. The methods developed have been validated via simulations of wide range of EM problems in both frequency and time domains such as wave propagation in various waveguide structures, scattering from deep cavities and cavity backed aperture antennas, scattering from dielectric and PEC objects, and radiation from horn antennas. In summary, the contributions of this dissertation are as follows: • The initial VGFEM framework has been modified to utilize VGFEM with advance meshing tools. A general methodology has been developed to define the overlapping partition of unity domains on polyhedral meshes. Without loss of generality, tech146 niques for defining the partition of unity and local approximation functions have been presented for brick and tetrahedral elements. • Techniques to include discontinuities in VGFEM have been developed for the analysis of multi-material problems. Two solutions have been proposed: (i) augmentation of basis functions with discontinuous basis functions, (ii) local domain decompositions with appropriate transmission conditions at material interfaces. Both approaches have been validated via various material interface problems. • Convergence and error characteristics of the method have been rigorously studied, especially dispersion/phase error characteristics. Since the analysis of the dispersion characteristics in overlapping scheme of VGFEM is not as simple as in FEM, a semianalytic technique has been developed by following the traditional dispersion analysis. • Perfectly Matched Layers (PML) and absorbing boundary conditions have been integrated with the VGFEM framework to analyze open domain EM problems. Anisotropic interpretation of the PML has been used in this extension. The capabilities and performance of the VGFEM with ABC and PML have been studied via a number of applications. • VGFEM has been integrated with boundary integrals to analyze scattering from deep cavities and cavity backed antennas. Since the definition of an auxiliary basis function space on the PU domains is not possible in VGFEM, both methods are coupled by forcing the volumetric vector basis functions to satisfy the boundary conditions. • Features of VGFEM such as the use of different basis functions or polynomial orders in different regions of a problem domain have been demonstrated via various applications. 147 Use of mixed basis space or polynomial orders has been easily achieved by assigning a basis type or polynomial type to each PU domain. In that way, more accurate and faster convergent solutions are obtained. • VGFEM has been integrated with other FEM methods such as the classical finite element methods, domain decompositions and discontinuous Galerkin methods. This integration has been achieved by partitioning a problems domain into sub-domains, assigning a method to each sub-domain, and then applying appropriate transmission conditions at the sub-domain interfaces. This technique has also been used to approximate the fields in inhomogeneous domains. • Time domain-VGFEM has been developed to solve transient electromagnetic problems. The Newmark method has been used for time marching scheme. Both spatial and temporal convergences of the method have been shown. 6.2 Future Work Although the method has been advanced through theoretical and computational developments, there are still some problems that need to be solved in order to make VGFEM a powerful and well established technique for solving complex EM problems. These problems are summarized as follows: • Applications of VGFEM to complex and large EM problems using its features in EM simulations: The capabilities of the method have been shown for some practical problems. Its capabilities should be further demonstrated via the analysis of realistic multiscale problems. 148 • Applications of discontinuous Galerkin -VGFEM for complex multi material EM problems: Hybrid DG-VGFEM is a candidate for powerful hybrid methods as advantages of both methods are exploited in EM simulations. We have developed this technique for solving wave equation in frequency domain problems. However, DG methods are well suited for the coupled Maxwell equations that results in efficient parallel solutions. Thus, development of DG-VGFEM using the coupled Maxwell equations would be an important advancement. • Applications of absorbing boundary conditions and perfectly matched layers to solve complex open domain EM problems: ABCs and PMLs have not been studied in details in this thesis; therefore, these domain truncation techniques need detailed analysis within the VGFEM framework, and their capabilities need to be demonstrated through the analysis of complex problems. • Applications and further developments of the time domain VGFEM: Currently, only the Newmark time scheme has been implemented. Different time marching schemes, and their convergence and stability characteristics should be studied. In addition, the method should be tested for large and complex time domain problems. • Development of methodologies for solving ill-conditioned VGFEM system and hybrid VGFEM methods: The overlapping nature of VGFEM and liner dependency are the causes of such ill-conditioned systems. Although techniques to eliminate the linear dependency have been developed [25], the matrix system is still badly conditioned due to the overlapping framework. One can make the system better conditioned by reducing the overlapping among the PU domains. A detailed study for the matrix solutions of VGFEM and hybrid FE-VGFEM is needed. 149 BIBLIOGRAPHY 150 BIBLIOGRAPHY [1] A. Weisshaar, S. M. Goodnick, and V. K. Tripathi, “A rigorous and efficient method of moment solutions for curved waveguides bends,” IEEE Trans. on Ant. and Prop., vol. 40, pp. 2200–2206, 1992. [2] J. M. Jin and L. Volakis, “A finite element- boundary integral fromulation for scattering by three-dimensional-cavity-backed apertures,” IEEE Trans. on Ant. and Prop., vol. 39, pp. 97–104, 1991. [3] J. Liu and J. M. Jin, “A special higher order finite-element method for scattering by deep cavities,” IEEE Trans. on Ant. and Prop., vol. 48, pp. 609–703, 2000. [4] K. Barkeshli and L. Volakis, “Scattering from narrow rectangular filled grooves,” IEEE Trans. on Ant. and Prop., vol. 39, pp. 804–810, 1991. [5] J. M. Jin, J. L. Volakis, and J. D. Collins, “A finite-element-boundary-integral method for scattering and radiation by two- and three-dimensional structures,” IEEE Ant. and Prop. Magazine, vol. 33, pp. 22–32, 1991. [6] P. Cornet, R. Dusseaux, and J. Chandezon, “A rigorous and efficient method of moment solutions for curved waveguides bends,” IEEE Trans. Microwave Theory Tech., vol. 47, pp. 965–972, 1999. [7] J. M. Rebollar, J. Esteban, and J. E. Page, “Fullwave analysis of three and four-port rectangular waveguide junctions,” IEEE Trans. Microwave Theory Tech., vol. 42, pp. 256–263, 1994. [8] R. H. MacPhie and K.-L. Wu, “Scattering at the junction of a rectangular waveguide and a larger circular waveguide,” IEEE Trans. on Microwave Theory and Tech., vol. 43, pp. 2041–2045, 1995. 151 [9] Y. Li and J.-M. Jin, “A vector dual-primal finite element tearing and interconnecting method for solving 3-d large-scale electromagnetic problems,” IEEE Trans. on Ant. and Prop., vol. 54, pp. 3000–3009, 2006. [10] K. Zhao, V. Rawat, S.-C. Lee, and J.-F. Lee, “A domain decomposition method with nonconformal meshes for finite periodic and semi-periodic structures,” IEEE Trans. on Ant. and Prop., vol. 55, pp. 2559–2570, 2007. [11] J. S. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods. New York: Springer, 2007. [12] P. Houston, I. Perugia, and D. Schotzau, “hp-dgfem for maxwell’s equations,” Numerical Mathematics and Advanced Applications ENUMATH 2001, pp. 785–794, 2003. [13] J. M. Jin, The Finite Element Method in Electromagnetics. New York: Wiley, 2002. [14] O. C. Zienkiewicz and K. Morgan, Finite Elements and Approximation. Wiley, 1983. New York: [15] C. A. M. Duarte and J. T. Oden, “Hp clouds - a meshless method to solve boundary value problems,” TICAM Technical Report 95-05, The University of Texas at Austin, Tech. Rep., 1995. [16] T. Belytschko, Y. Y. Lu, and L. Gu, “Element-free galerkin methods,” International Journal For Numerical Methods In Engineering, vol. 37, pp. 229–256, 1994. [17] T. Strouboulis and R. Hidajat, “Partition of unity method for helmholtz equation: qcovergence for plane-wave and wave-band local bases,” Applications of Mathematics, vol. 51, pp. 181–204, 2006. [18] I. Babuska and J. M. Melenk, “The partition of unity method,” International Journal For Numerical Methods In Engineering, vol. 40, pp. 727–758, 1997. [19] T. Strouboulis, I. Babuska, and R. Hidajat, “The generalized finite element method for helmholtz equation: Theory, computation, and open problems,” Compt. Methods Appl. Mech. Engrg, vol. 195, pp. 4711–4731, 2006. [20] J. M. Melenk and I. Babuska, “The partition of unity finite element method: Basic theory and applications,” Comput. Methods Appl. Mech. Eng, vol. 139, pp. 289–314, 1996. 152 [21] C. Lu and B. Shanker, “Hybrid boundary integral-generalized (partition of unity) finiteelement solvers for the scalar helmholtz equation,” IEEE Trans. on Magnetics, vol. 43, pp. 1002–1012, 2007. [22] L. Proekt and I. Tsukerman, “Method of overlapping patches for electromagnetic computation,” IEEE Trans. Magn., vol. 38, p. 741744, 2002. [23] C. Lu and B. Shanker, “Solving boundary value problems using the generalized (partition of unity) finite element method,” IEEE APS 2005, vol. 1B, pp. 125–128, 2005. [24] ——, “Development of hybrid boundary integral-generalized (partition of unity) finite element solvers for scalar problems,” IEEE APS 2005, vol. 1B, pp. 129–132, 2005. [25] ——, “Generalized finite element method for vector electromagnetic problems,” IEEE Trans. on Ant. and Prop., vol. 55, pp. 1369–1381, 2007. [26] O. Tuncer, C. Lu, N. Nair, B. Shanker, and L. C. Kempel, “Further development of vector generalized finite element method and its hybridization with boundary integrals,” IEEE Trans. on Ant. and Prop., vol. 58, pp. 887–899, 2010. [27] I. Tsukerman, “General tangentially continuous vector elements,” IEEE Trans. Magn., vol. 39, p. 12151218, 2003. [28] C. Lu, B. Shanker, and E. Michielssen, “Development of generalized finite element method for vector electromagnetic problems,” IEEE APS 2006, pp. 2813–2816, 2006. [29] M. Griebel and M. A. Schweitzer, “A partical-partition of unity method (part ii: Efficient cover construction and reliable integration),” SIAM J. Sci. Comp., vol. 23, pp. 1655–1682, 2002. [30] M. A. Khayat and D. R. Wilton, “Numerical evaluation of singular and near-singular potential integrals,” IEEE Trans. on Ant. and Prop., vol. 53, pp. 3180–3190, 2005. [31] Y. Xu, C.-F. Wang, Y.-B. Gan, and F.-G. Hu, “A domain decomposition method for analysis of electromagnetic scattering from cavities,” IEEE APS 2005, pp. 1–4, 2005. [32] C. J. Reddy, M. D. Deshpande, C. R. Cockrell, and F. B. Beck, “Electromagnetic scattering analysis of a three-dimensional-cavity-backed aperture in an infinite ground plane using a combined finite element method/method of moments approach,” NASA Technical Paper 3544, Tech. Rep., 1995. 153 [33] T.-M. Wang and H. Ling, “Electromagnetic scattering from three-dimensional cavities via a connection scheme,” IEEE Trans. on Ant. and Prop, vol. 39, pp. 1505–1513, 1991. [34] A. Deraemaeker, I. Babuska, and P. Bouillard, “Dispersion and pollution of the fem solution for the helmholtz equation in one, two and three dimensions,” International Journal For Numerical Methods In Engineering, vol. 46, pp. 471–499, 1999. [35] R. Lee and A. C. Cangellaris, “A study of discritization error in the finite element approximation of wave solutions,” IEEE Trans. on Ant. and Prop., vol. 40, pp. 542– 549, 1992. [36] G. S. Warren and W. R. Scott, “Numerical dispersion of higher order nodal elements in the finite element method,” IEEE Trans. on Ant. and Prop, vol. 44, pp. 317–320, 1996. [37] A. Chatterjee, J. M. Jin, and J. L. Volakis, “Edge-based finite elements and vector abcs applied to 3-d scattering,” IEEE Trans. on Ant. and Prop., vol. 41, pp. 221–226, 1993. [38] Z. Lou and J. M. Jin, “Modeling and simulation of broad-band antennas using the timedomain finite element method,” IEEE Trans. on Ant. and Prop., vol. 53, pp. 4099–4110, 2005. [39] K. Umashankar, A. Taflove, and A. M. Rao, “Electromagnetic scattering by arbitrary shaped three-dimensional homogeneous lossy dielectric objects,” IEEE Trans. on Ant. and Prop., vol. 34, pp. 758–766, 1986. 154