A Full-Wave Discontinuous Galerkin Time-Domain Finite Element Method for Electromagnetic Field Mode Analysis

The accurate analysis and the characterization of guided-wave structures are important steps in the design of microwave and optical systems. A numerical methodology based on a three-dimensional discontinuous Galerkin time-domain (DGTD) finite element method is used to estimate the resonant frequencies of various structures in a single run, with several polychromatic sources (dipole, plane wave, etc). The proposed approach leads to simulations that are fast and accurate, this being the result of the local discretization strategy of the full-wave time-domain Maxwell’s equations on unstructured tetrahedral meshes that is easy to parallelize. The stability of the overall numerical scheme is obtained using an upwind numerical flux and an implicit second-order accurate time integration scheme. The numerical methodology is verified based on two problems with closed-form solutions. First, homogeneous and non-homogeneous waveguides are studied using the proposed methodology. The WR-284 waveguide is used, and cut-off frequencies of the first 7 modes are computed and compared with analytical solutions. Propagation characteristics of mode-selective transmission lines (MSTL) are then investigated. These MSTL simulation results are verified using measurement results of fabricated MSTL, microstrip, and RWG prototypes up to 110 GHz. Simulation results show the flexibility, the robustness, and the high accuracy of the proposed numerical strategy.


I. INTRODUCTION
The analysis of electromagnetic fields has been the subject of numerous studies over the past decades. The main two conductor structures that can support verious TE, TM, and TEM mode wave propagations are waveguides and transmission lines. In this paper, a numerical method is proposed for analyzing waveguides and mode-selective transmission line (MSTL).
Several authors have analyzed waveguides problem in the literature. Their shape and cut-off frequency influence the The associate editor coordinating the review of this manuscript and approving it for publication was Guido Lombardi . design of microwave, millimeter-wave, and optic components, such as filters [1], [2], [3], [4], [5], tuners [6], [7], phase shifters [8], [9], attenuators [10], and circulator elements [11], [12]. Waveguides can be filled with a homogeneous or non-homogeneous dielectric medium. Propagation characteristics of MSTL are then investigated. MSTL is studied theoretically and experimentally by some researchers [13], [14]. This type of transmission line has frequency-dependent mode-switching behavior and is among low-loss and lowdispersion transmission line [13] The MSTL is selected to show the accuracy of the proposed numerical method to analyze the losses in materials and caused by leakage and radiation in the interface between different layers in multiplayer VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ transmission lines. For this purpose, the time-domain fullwave Maxwell' equations are discretized on bounded 3D domains.
Maxwell's equations can give multiple solutions for a given problem. With homogeneous waveguides, two different modes exist: the transverse electric mode (TE) and the transverse magnetic mode (TM). For homogeneous rectangular or circular waveguides, the partial differential equations can be solved analytically. But for waveguides with complex shapes, or filled with non-homogeneous materials, we cannot find an analytic solution. For multilayered printed transmission line, computing the losses due to radiation and leakage are also challenging. Numerical methods, such as the finite difference method [15], [16], [17], the boundary element method [18], [19], [20], the method of moments [21], [22], [23], or the finite element method [24], [25], [26] must therefore be considered.
Various software packages, such as HFSS, ADS, CST, and FEKO, are used to perform the numerical analysis for Maxwell's equations. They are in the frequency domain (FD), or in the time domain (TD). FD methods are often more computationally costly for computing wideband frequency responses since, for each frequency, the simulation must be run completely. Some of commercial software packages are designed based on modified version of Maxwell's equations (green's functions, Helmholtz equation, H-formulation, etc. ), therefore they disregarded important physics phenomena that can only be modeled using full-wave Maxwell's equations. Commercial simulation software is usually designed to study either low frequency, or high frequency problems. Using the full-wave Maxwell's equations will make the simulation software usable for various frequencies and media.
Some software packages are based on the finite-difference time-domain (FDTD) numerical method which is not very accurate for complex geometries. The finite element method (FEM) is the most popular numerical method that can be used to discretize Maxwell's equations on unstructured meshes. The FEM is a useful simulation tool for studying problems with complex geometries. But waveguides with anisotropic materials cannot be handled easily. The numerical simulation of these problems can lead to unstable approximations [27], [28]. The development of stable high-order finite element discretizations for Maxwell's equations is still an active research topic.
In order to help control numerical instabilities, some studies propose to introduce modifications to Maxwell's equations, which lead to nonphysical discrete approximations. These techniques are known for not completely succeeding in eliminating numerically induced oscillations [29], [30], [31]. Edge elements, also known as Whitney elements of degree 1, provide a cure to many problems associated with nodal-based finite element approximations [32], [33], [34]. First, they make it easier to manage discontinuous permittivities and permeabilities, since edge elements only enforce the continuity of the tangential component of the dependant variables, allowing the normal component of the vector fields to be discontinuous. They also provide more accurate approximations of vector fields near geometric singularities, such as in regions in the vicinity of sharp corners. Finally, since the numerical approximations are taken in the proper functional space, spurious modes associated with nonzero eigenvalues are reduced or eliminated. These properties come at a price, since edge elements necessitate an increased number of degrees of freedom when compared to nodal discretization techniques, which translates into simulations that necessitate more time and memory.
A fast and accurate numerical methodology based on a three-dimensional discontinuous Galerkin time-domain (DGTD) finite element method with vector basis functions (Whitney elements of degree 1) is proposed for analyzing waveguides and mode-selective transmission line (MSTL). The DGTD method does not necessitate the assembly of the contributions of the elements of the mesh in a large matrix. Instead, small local problems are discretized on each element of the mesh, which makes the method easier to parallelize [35], [36]. The numerical approximation of the discrete problem is therefore discontinuous, which adds to the flexibility provided by the Whitney element. Both homogeneous and non-homogeneous waveguides are studied to verify the proposed numerical methodology, where the cut-off frequency is computed for various modes in a single run. Propagation characteristics of MSTL are also investigated in this work. These analysis results are verified by the measured results of fabricated prototypes of MSTL, microstrip, and RWG up to 110 GHz.

A. MAXWELL's EQUATIONS FOR ELECTROMAGNETIC FIELDS MODELING
The differntial form of Maxwell's equations in a limited region of interest are expressed as where H is the magnetic field, E is the electric field, ε and µ represent the permittivity and permeability of the material, ρ is the charge density, and J is the current density. Dirichlet and Neumann boundary conditions must be defined as boundary conditions to have complete definition of the problem. Dirichlet condition imposed the value and Neumann condition imposed the normal derivative of the electromagnetic fields on the surface of the region.

B. INITIAL AND BOUNDARY CONDITIONS
Since the problems under study are time dependent, initial conditions E(x, t 0 ) = E 0 (x) and H(x, t 0 ) = H 0 (x), for all x ∈ , allow us to compute the discrete variables E h (x, t) and The appropriate boundary values must be imposed on the surface area or perpendicular to the surface depending on the problem under study. They can be defined for the electromagnetic fields in media 1 and 2 as where n is the unit normal vector to and g E , h E , g H and h H are known functions defined on × I (or with a subset of ). A well known example is given by the perfect electric conductor (PEC) boundary condition: which can be used, as an example, with conductors at microwave frequencies. Another example is the perfect magnetic conductor (PMC) boundary condition: which is often used as a symmetry boundary condition.

C. NONDIMENSIONALIZATION
Dimensionless forms are usually used to help identify the dominant physics of a given system. But from a numerical modeling point of view, it can also act as a natural preconditionner. The order of magnitude of each term of system (1) can vary significantly, which can lead to an ill-conditioned discrete problem. Given the reference quantities t 0 , 0 , µ 0 , ε 0 ,, ρ 0 , E 0 and H 0 , we define the nondimensional variables The dimensionless form of system (1), after replacing J by σ E, then becomes: H 0 t 0 , for conducting materials 2 ≈ 0, and 3 = E 0 0 σ 0 H 0 . In this study, for a unit reference magnetic field strength of H 0 = 1 A/m, we define the reference electric field strength as E 0 = H 0 Z 0 for an impedance of Z 0 = 120π , the reference permeability µ 0 is chosen as 4π × 10 −7 H/m, and the reference permittivity ε 0 is given by 8.85 × 10 −12 F/m. Finally, t 0 and 0 are problem dependent. The nondimensional quantities will be written without ''∼'' in the remainder of the paper.

D. NUMERICAL METHODOLOGY
Industrial problems often involve complex geometries. This leads us to base this numerical study on the finite element method (FEM). Many authors discretize a modified version of system (1), which sometimes leads to ill-conditioned discrete problems. This explains why we choose to perform a direct finite element discretization of the third and fourth equations of system (1) using the discontinuous Galerkin FEM, also known as the discontinuous Galerkin time-domain (DGTD) method. In addition to the fact that it can provide a stable discrete problem, it's local nature makes it a numerical method that is easy to parallelize, which results in fast solvers for the discretization of Maxwell's equations. After discretizing the domain of definition of the problem using tetrahedral elements, we obtain the finite element mesh h = The base functional space associated with the weak forms of Maxwell's equations is expressed as The classical approach consists in using the piecewise N-th order polynomial Lagrange approximation of the form where the E j and H j are scalar degrees of freedom in each node which are four for linear and ten for quadratic polynomial function. The nodal elements have the disadvantage that imposes continuous functions of the spatial variable, and they also cause difficulties near sharp corners. To overocome these problems the vectorial approximations or vectorial degrees of freedom E j and H j are used instead of nodal degree of freedoms. Linear edge element, also called Whitney element, is applied which is the most popular and simlest basis function were proposed by Nédélec [40].
To obtain the weak form of Faraday's law and Ampere's law, the test functions v and w must be multiplies in the equations of system (4) taken in an appropriate functional space, and we integrate over each element. The details about the process of obtainng the weak form of Maxwell's equation can be found in [41]. The weak form associated with the third and fourth equations of system (1) are where ∂K is the boundary of element K. The right-hand side of the equations in system (6) The H * , and E * are expressed in the form where {{u}} represents the average and [u] is used for defining the jumps along a normal vector between two elements. These notations are defined as here ''+'' and ''−'' refers to the left and right side of the interface which is shown in Fig. 1. The right-hand side of the equations in system (6) can be written as The M, S, and F matrices can be defined as where all these matrices are local and based on integrals of the basis functions.
The second-order accurate backward differentiation formula (BDF) is used to discretize the transient derivative of Maxwell's equations. The general formula for a second order BDF can be written as where t denotes the step size. The BDF is a populare implicit methods based on Lagrange interpolation polynomial for the numerical integration of ordinary differential equations, see [37]. More details of the DG-based strategy can be found in [38] and [39].

III. NUMERICAL AND EXPERIMENTAL RESULTS
The proposed numerical methodology is assessed based on a two-part analysis. The developed DGTD solver is first used to study rectangular homogeneous and non-homogeneous waveguides with known cut-off frequencies, at two different TE and TM modes. The computed numerical results are compared with closed-form analytical solutions. This approach is then used to study the behavior of MSTLs, which meets the stringent requirements for high-quality guided-wave signal propagation. The proposed MSTL is implemented and characterized for an experimental verification.

A. WAVEGUIDE
The rectangular waveguide illustrated in Fig. 2 is first studied. It is assumed that it is homogeneous with perfectly electric conducting (PEC) walls (cf. eqn (2)). A basic sinusoidal plane wave is often used as a source. Although plane wave source can be practical in some situations, it is not worth running a complete time-domain simulation to only get a system response for a single frequency. Frequency-domain methods are usually better suited in these situations. In order to obtain a wider frequency spectrum with good control over the Gaussian function, modulation using a sinusoidal function is considered, where the duration σ , centered at t 0 , is chosen with a frequency range of 0 -10 GHz, ω = π(f max + f min ) = 10π × 10 9 , and t 0 = 2 × σ = 0.2 × 10 −9 . The boundary conditions at the input (for z = 0 mm) and the output (for z = 150 mm) of the waveguide are therefore given by Fourier transforms of the time response is used to compute the frequency response of rectangular waveguides. The Gaussian function and its Fourier transforms are illustrated in Fig. 3. A Gauss pulse is applied at a corner of the geometry, and the components of the electrical and magnetic fields are computed at the center of the waveguide. Numerical and analytical time and frequency responses of the TE mode for three components of the electromagnetic field (H z , E x , and E y ) are given and illustrated in Fig. 4. Numerical and analytical values of the TM mode for three components of the electromagnetic field (E z , H x and H y ) are given and illustrated in Fig. 5. Each peak shown in the 0 -10 GHz range of the frequency response represents the mode cutoff frequency. All quantities of Figs 4 and 5 are normalized. The computed results agree with analytical solutions. As expected, the mesh grid size and the time-step size influence the accuracy of the computed approximations. In our case, the mesh grid size used is h = 2 mm, for a total of 0 tetrahedra, and the time-step size is t = 5 ps, which result in a relative nodal error of less than 1%. The second verification problem consists of a rectangular waveguide with dimensions 8 × 4 × 32 cm, which is excited for the lowest mode TE10 at a frequency of f = 3 GHz. In this case, the TE10 electromagnetic field is such that Two cases are considered. The first waveguide is filled with air (ε r = 1 and µ r = 1). The relative permittivity of the second waveguide varies linearly from 1 to 5 over the length of the waveguide, i.e. ε r (x) = 1 8 z + 1. The computational domain h is discretized into 40, 000 tetrahedral elements, and the time-step size used is t = 10 −3 . The reference quantities for the normalization are: H 0 = 1, E 0 = 120π, 0 = 0.1, f 0 = t −1 0 = 3 GHz, and 1 = 2 = 1. Fig. 6 shows the three different components of electromagnetic field (H x , H z , E y ) for air field waveguides and non-homogeneous waveguides. The phase constant can be calculated from the magnitude of H x components.
All the values in figures are normalized amount. As a source, the basic sinusoidal plane wave is chosen. Dielectric filled metallic waveguide has the same property as the airfilled one with the only exception that the cut-off frequency for the TE modes is varies in different part of waveguide due to change od dielectric properties.
The phase constant calculations are done in nine different slices where the magnetic field has their maximum values which are presented in Fig. 7a. The electromagnetic field component (H x and E y ) along z-axis are shown in Figs. 7a, 7b. The results in Fig. 7 are in good agreement with the analytical solutions for TE10 mode for different values of relative dielectric permittivities.

B. MSTL
In this section, the numerical results of MSTL are shown to experimentally verify the capability of the proposed DGTD methods. To show the accuracy of the method, the cross-sectional electromagnetic (EM) fields in MSTL are examined, and its propagation constant is calculated and then compared with the measured and HFSS simulated results. Fig. 8 shows the structure of the proposed MSTL. The geometry dimensions and dielectric substrate properties are presented in Table 1. It consists of a homogeneous dielectric substrate of thickness h and width d with ε r . The top and bottom planes are metals with the thickness of t and two parallel slots of size s on the top. The metal strip at the center of MSTL has the width of w. The dimensions are chosen in such a way that MSTL supports a dominant quasi-TEM mode of microstrip at lower frequency range covering DC. At the cut-off frequency, a mode conversion from the propagating quasi-TEM mode to a quasi-TE10 mode takes place. In fact, it behaves similarly to the microstrip line cases over the low frequency region and the same way as a rectangular waveguide over the higher frequency region. For experimental verification, a standard printed circuit board (PCB) process technology was adopted to implement the MSTL. The MSTL prototype is fabricated on RT/duroid 6010LM laminate from Rogers Corporation. Measurements are carried out from 10 MHz-to-110 GHz using a vector network analyzer.  More details about the fabricated MSTL can be find in [14] and [42] For numerical verification, the solution domain is discretized into 229, 282 tetrahedral elements. The problem geometry is decomposed into eight subdomains and Message Passing Interface (MPI) and domain decomposition techniques are applied to make the simulation distribute on various processors of a parallel computer [43]. The Whitney basis function or edge elements are used as basis functions because of stability and consistency properties specially for eliminating spurious solutions at the interface of two layers in MSTL systems. In oreder to have a good compromise between accuracy and computational cost the time step of 0.5 × 10 −3 is chosen.  The different components of the magnetic field (H x , H y , H z ) and electric field (E x , E y , E z ) at two different frequencies of 100 GHz, and 10 GHz along x-axis are depicted in Figs. 10, and 11. As expected, all the components of electric and magnetic fields have a noticeable change in their distribution shape and magnitude in various time steps specially around slots. A sketch of the electric and magnetic field lines for the dominant quasi-TE10 waveguide mode of   MSTL at 100 GHz are shown in Fig. 12. It can be observed that MSTL should have the guided wave performance similar to that of a dielectric-filled RWG when frequency is higher than cut of frequency.
Phase and attenuation constant curves calculated using different method are plotted in Figs. 13a, 13b. As it shows, excellent agreement between numerical, HFSS, and experimental results is observed. It is confirmed that the low-dispersion and low-loss behavior of MSTL makes it an outstanding integrated waveguide in support of high-performance superbroadband signal transmission and/or ultra-fast pulse propagation in a fully-integrated platform. In order to evaluate    dispersion error reduce notably by decreasing the mesh size (h) for second order nodal and Whitney basis functions. However, using Whitney edge element requires a large number of DoFs which resulting increase in computation cost (memory and time). To compensate this issue, the geometry is divided in not very fine mesh and Whitney basis function is applied. The presented method has no limitation for applications, and the processing time is all depend to the processor we used. For presented problem, we use 4 GHz Intel Core i7, memory: 32 GB 1867 MHz DDR3. For more complicated problem we need more power full processor and more memorry. This paper, for the first time, analyze the full-wave timedomain Maxwell's equations in a 3D structure for modeling electromagnetic field propagation in three layers MSTL system. The electromagnetic field components at the interface between layers are computed accurately. Several authors have studied the DGTD method for Maxwell's equations, and their results prove that the DGTD exhibits better consistency, stability and convergence properties when compared to the standard FEM. These advantages are because of the element-based computation nature of DGTD and having a local elements strategy [44], [45]. In this section, it is proven that the acceptable results can be obtained (in comparison to measurement results) with a relatively larger time step and mesh size (h = 0.01, and t = 0.5 × 10 −3 ) which demonstrate the high stability of DGTD in compared to FEM. Using MPI parallel computer library reduces the numerical computational time significantly.

IV. CONCLUSION
In this paper, a 3D time domain DGTD has been introduced to analyze the electromagnetic field propagation problem. As a first example, the electromagnetic field inside homogeneous and non-homogeneous waveguides are approximated and the resonant frequencies are estimated. The developed code has capability to give us all required information in a single run for different polychromatic sources. The accuracy and efficiency of present algorithms are assessed in comparison with the analytical and measurement results. In the second verification problem, the developed simulator is used for the first time to study the behavior of MSTL. The phase constants and attenuation are measured for specific MSTL geometry and compared with measurement and HFSS software results. The DGTD method can analyze accurately the losses in matrials and caused by leakage and radiation in the interface between different layers in multiplayer transmission lines. The DGTD is an implicit variant of FEM that equations are discretized using an element-by-element strateg. This property of DGTD leads to better convergence rates and lower computational cost specially for complex multilayered geometries. In addition, the proposed method can be used in various frequency ranges, complex geometries, non-homogeneous materials, different excitation modes, etc. MPI paralle computers technique makes the proposed implementation of the DGTD solver significantly faster. The combination of the DGTD method, Whitney basis function, and applying MPI techniques are the components of the proposed numerical methodology that makes the developed solver faster and more stable than commercial codes. The developed code was validated based on both theoretical and measurement results. He is also with the School of Information Science and Engineering, Ningbo University, Ningbo, China, on leave from his home institution, leading a special 5G and future wireless research program. He has authored or coauthored over 1100 refereed articles and a number of books/book chapters. He has filed more than 50 patents. His current research interests include substrateintegrated circuits, antenna arrays, field theory and joint field/circuit modeling, ultrafast interconnects, wireless power transmission and harvesting, megahertz-through-terahertz transceivers and sensors for wireless systems and biomedical applications, and modeling and design of microwave and terahertz photonic circuits and systems.
Dr He is teaching applied mathematics to engineering students. His current research interests include the development of new finite-element based numerical methodologies for modeling fluid mechanics problems found in free surface and turbulent flows, for modeling electromagnetism phenomena found in supraconductors and for the optimization of antennas, and for the study of magnetohydrodynamics. VOLUME 10, 2022