This article provides insight into the modeling of waveguides using the discontinuous Galerkin time-domain (DGTD) method. The spatial domain decomposition in the DGTD method is controlled by changing the resolution of the finite elements (FEs) which in turn have effect on the stable time marching step. Different resolutions of FEs support different resonant modes. This paper investigates the effects of the resolution of FEs and the effect of Final Time on appearance of low and high–frequency resonant modes in the waveguides. The experimentation is performed by considering the both transverse magnetic (TM) and transverse electric (TE) polarizations. The unstructured triangular mesh type is considered for the analysis. Integration in the temporal domain is achieved by using five-stage 4th order low-storage Runge-Kutta (LSERK) method and the stability of the numerical method is ensured by stable time marching step which is calculated by using Courant-Friedrichs-Levy (CFL) condition. The operation is performed on WR90 waveguide and the analytical and numeric values for various resonant modes are compared to validate the study.

Keywords: DGTD, modeling of waveguide, CFL, effect of resolution of finite elements, effect of final time, resonant modes of waveguide, WR90 waveguide

1D, one-dimensional; 2D, two-dimensional; 3D, three-dimensional; DGTD, discontinuous galerkin time-domain; Fes, finite elements; TM, transverse magnetic; TE, transverse electric; ODE, ordinary differential equations; LSERK, low-storage explicit runge-kutta; CFL, courant-friedrichs-levy; TD, time-domain; PEC, perfect electric conductor; FFT, fast fourier transformation

Waveguides have ability to handle very high power which makes them suitable for applications in RADAR and satellite communication. Unlike free space, rectangular waveguides support transverse magnetic (TM) and transverse electric (TE) modes – if air-filled or fully loaded with some dielectric material.¹ While modeling waveguides by using numerical methods, certain factors like digitization of the computational domain, accuracy and execution time, need considerations. The variation of these parameters can lead to unwanted results and accurate modeling cannot be performed if the significance of these parameters being used in the numerical method remains unknown. The numerical methods have ability to deal with complex electromagnetic (EM) problems which cannot be molded using traditional analytical methods. These methods require due attention for their development otherwise the variation in parameters can generate unwanted results. Discontinuous Galerkin time-domain (DGTD) method is one of the recent methods in computational electromagnetic (CEM). It is being widely studied^2–9 since its first appearance to solve the neutron transport equation using unstructured triangular mesh.¹⁰ The DGTD method is the combination of finite element method (FEM) and finite volume method (FVM)^11,12 but offers better performance as compared to these two methods.¹³ Since the space discretization of the DGTD has effect on the bandwidth coverage,^14,15 so it is vital to study the effect of resolution of the finite elements (FEs) on the appearance of the resonant modes of the waveguides. Also, the value of Final Time, the time required to terminate the simulation sequence, has its own significance regarding the resonant modes. To the knowledge of authors, these kinds of investigations have never been performed before and require due attention in order to properly model and analyze waveguides for their resonant frequencies. The organization of this article is as follows: Section ‎5 provides the formulation required for the study. It includes the formulation of the DGTD method for the both polarization i.e., TM and TE, upwind numerical fluxes are used to derive the expression, and truncation of the boundary by using perfect electric conductor (PEC) condition is provided. The formula for resonant mode calculation using exact method is also given. Section 6 provides the numerical analysis and discusses the effect of resolution of FEs and Final Time on the resonant modes, followed by comparison of the values of these modes by using the exact and the numerical methods. Finally, the conclusions are drawn in section ‎7.

The DGTD Method Consider Maxwell’s equation in three-dimensions (3D), as given in Eq.(1):

$\epsilon \frac{\partial {\rm E}}{\partial t}-\nabla \times {\rm H}=0,\mu \frac{\partial {\rm H}}{\partial t}+\nabla \times {\rm E}=0$ (1)

where

1. E=(Ex, Ey, Ez) and H=(Hx, Hy, Hz) are the electric and magnetic fields;
2. respectively, ε is the permittivity and μ is the permeability.

To solve this numerically, the standard procedure is to treat space and time domain separately. In the standard DGTD method, space is discretized into K non overlapping finite elements (FEs). The local solutions are estimated for each element and the global solution is obtained by integrating all the local solutions. The communication between neighboring elements is ensured by using numerical fluxes. From the theory of Riemann solvers,^16,17 the penalty terms inflowing into the normal $\widehat{n}$ are given below:

$\widehat{n}\cdot (F_E-F_E^{*})=[\widehat{n}\times E-{(\widehat{n}\times E)}^{*}]=-\frac{1}{2\left\{\left\{Z\right\}\right\}}\widehat{n}\times \left[Z^{+}(H^{-}-H^{+})-\alpha \widehat{n}\times (E^{-}-E^{+})\right]$ (2)

$\widehat{n}\cdot (F_H-F_H^{*})=-[\widehat{n}\times H-{(\widehat{n}\times H)}^{*}]=\frac{1}{2\left\{\left\{Y\right\}\right\}}\widehat{n}\times \left[Y^{+}(E^{-}-E^{+})+\alpha \widehat{n}\times (H^{-}-H^{+})\right]$ (3)

where

1. $F_{E,H}^{*}$ is the numerical flux;
2. The ‘+’ and ‘–’ signs in the superscript indicate the exterior and interior of an interface, respectively;
3. The $Z^{\pm }=\sqrt{{\mu }^{\pm }/{\epsilon }^{\pm }}={(Y^{\pm })}^{-1}$ is the impedance;
4. ‘α’ controls the dissipation (i.e., α = 1 and 0 are for classic upwind and non-dissipative central fluxes respectively).

The function {{∙}} corresponds to the average value at an interface due to exterior and interior values e.g., $\left\{\left\{Z\right\}\right\}=\left(Z^{+}+Z^{-}\right)/2$ . After some mathematical calculations and evaluating the vector identities, Eq. (2) and Eq. (3) generate the following numerical fluxes for the electric and magnetic fields, respectively:

$\begin{array}{l}\widehat{n}\cdot (E-E^{*})=-\frac{1}{2\left\{\left\{Z\right\}\right\}}\left[\widehat{n}\times Z^{+}dH-\left(\widehat{n}\cdot dE\right)\alpha \widehat{n}+\left(\alpha \right)dE\right],\because dU=U^{-}-U^{+};forU=[E,H]\\ \end{array}$ (4)

$\widehat{n}\cdot (H-H^{*})=\frac{1}{2\left\{\left\{Y\right\}\right\}}\left[\widehat{n}\times Y^{+}dE+\alpha \left(\left(\widehat{n}\cdot dH\right)\widehat{n}-dH\right)\right],\because dU=U^{-}-U^{+};forU=[E,H]$ (5)

The above two equations are for 3D and very useful to obtain the numerical flux for any component along their respective axis.

1) TMz case in 2D: Maxwell’s equations given in Eq. (1) takes the following form:

$\begin{array}{l}\mu \frac{\partial H^x}{\partial t}=-\frac{\partial E^z}{\partial y},\\ \mu \frac{\partial H^y}{\partial t}=\frac{\partial E^z}{\partial x},\\ \epsilon \frac{\partial E^z}{\partial t}=\frac{\partial H^y}{\partial x}-\frac{\partial H^x}{\partial y}.\end{array}$ (6)

In 2D, nx and ny exist while nz = 0 and for the TMz mode we have only Hx, Hy and Ez components. The numerical fluxes for each component are obtained by using Eq. (4) and Eq. (5)and are given as follows:

$\widehat{n}\cdot (H_x-H_x^{*})=\frac{1}{2\left\{\left\{Y\right\}\right\}}\left[n_y{Y}^{+}d{E}_z+\alpha \left(ndotdH.n_x-d{H}_x\right)\right],$ (7)

$\widehat{n}\cdot (H_y-H_y^{*})=\frac{1}{2\left\{\left\{Y\right\}\right\}}\left[-n_x{Y}^{+}d{E}_z+\alpha \left(ndotdH.n_y-d{H}_y\right)\right],$ (8)

$\widehat{n}\cdot (E_z-E_z^{*})=\frac{1}{2\left\{\left\{Z\right\}\right\}}\left[n_y{Z}^{+}d{H}_x-n_x{Z}^{+}d{H}_y-\alpha d{E}_z\right].$ (9)

where

$ndotdH=\left[n_{x}d{H}_x+n_{y}d{H}_y\right]$ is the dot product between dH and $\widehat{n}$ .

$\begin{array}{l}\frac{d{H}_h^x}{dt}=\frac{1}{\mu }\left(-D_y{E}_h^z+{\left(JM\right)}^{-1}{\displaystyle \underset{\partial D^k}{\int }\left(\frac{{\widehat{n}}_y{Y}^{+}\left[E_h^z\right]+\alpha \left({\widehat{n}}_x〚H_h〛-\left[H_h^x\right]\right)}{2\left\{\left\{Y\right\}\right\}}\right)}l(r)dr\right),\\ \frac{d{H}_h^y}{dt}=\frac{1}{\mu }\left(D_x{E}_h^z+{\left(JM\right)}^{-1}{\displaystyle \underset{\partial D^k}{\int }\left(\frac{-{\widehat{n}}_x{Y}^{+}\left[E_h^z\right]+\alpha \left({\widehat{n}}_x〚H_h〛-\left[H_h^x\right]\right)}{2\left\{\left\{Y\right\}\right\}}\right)}l(r)dr\right),\\ \frac{d{E}_h^z}{dt}=\frac{1}{\epsilon }\left(D_x{H}_h^y-D_y{H}_h^x+{\left(JM\right)}^{-1}{\displaystyle \underset{\partial D^k}{\int }\left(\frac{{\widehat{n}}_y{Z}^{+}\left[H_h^x\right]-{\widehat{n}}_x{Z}^{+}\left[H_h^y\right]-\alpha \left[E_h^z\right]}{2\left\{\left\{Z\right\}\right\}}\right)}l(r)dr\right).\end{array}$ (10)

where the subscript

1. h indicates to the approximate solution;
2. l(r) is two-dimensional Lagrange polynomialL
3. D, J and M are differential, Mass, and Jacobian matrices, respectively.

The details on these parameters can be found in.¹⁸ The functions shown in the brackets $〚\cdot 〛$ , $[\cdot ]$ and $\left\{\{\cdot \}\right\}$ are characterized for any parameter q as follows:

$\left[q\right]=q^{-}-q^{+},〚q〛=\widehat{n}·q,\left\{\left\{q\right\}\right\}=(q^{-}+q^{+})/2$

2) TEz case in 2D: The set of Maxwell’s equation, numerical fluxes of individual component and their corresponding semidiscrete DGTD formulation for the TEz mode are summarized in Eq. (11) to Eq. (15).

$\begin{array}{l}\epsilon \frac{\partial E^x}{\partial t}=\frac{\partial H^z}{\partial y},\\ \epsilon \frac{\partial E^y}{\partial t}=-\frac{\partial H^z}{\partial x},\\ \mu \frac{\partial H^z}{\partial t}=\frac{\partial E^x}{\partial y}-\frac{\partial E^y}{\partial x}.\end{array}$ (11)

$\widehat{n}\cdot (E_x-E_x^{*})=\frac{1}{2\left\{\left\{Z\right\}\right\}}\left[-n_y{Z}^{+}d{H}_z+\alpha \left(ndotdE.n_x-d{E}_x\right)\right],$ (12)

$\widehat{n}\cdot (E_y-E_y^{*})=\frac{1}{2\left\{\left\{Z\right\}\right\}}\left[n_x{Z}^{+}d{H}_z+\alpha \left(ndotdE.n_y-d{E}_y\right)\right],$ (13)

$\widehat{n}\cdot (H_z-H_z^{*})=\frac{1}{2\left\{\left\{Y\right\}\right\}}\left[n_x{Y}^{+}d{E}_y-n_y{Y}^{+}d{E}_x-\alpha d{H}_z\right].$ (14)

$\begin{array}{l}\frac{d{E}_h^x}{dt}=\frac{1}{\epsilon }\left(D_y{H}_h^z+{\left(JM\right)}^{-1}{\displaystyle \underset{\partial D^k}{\int }\left(\frac{-{\widehat{n}}_y{Z}^{+}\left[H_h^z\right]+\alpha \left({\widehat{n}}_x〚E_h〛-\left[E_h^x\right]\right)}{2\left\{\left\{Z\right\}\right\}}\right)}l(r)dr\right),\\ \frac{d{E}_h^y}{dt}=\frac{1}{\epsilon }\left(-D_x{H}_h^z+{\left(JM\right)}^{-1}{\displaystyle \underset{\partial D^k}{\int }\left(\frac{{\widehat{n}}_x{Z}^{+}\left[H_h^z\right]+\alpha \left({\widehat{n}}_y〚E_h〛-\left[E_h^y\right]\right)}{2\left\{\left\{Z\right\}\right\}}\right)}l(r)dr\right),\\ \frac{d{H}_h^z}{dt}=\frac{1}{\mu }\left(D_y{E}_h^x-D_x{E}_h^y+{\left(JM\right)}^{-1}{\displaystyle \underset{\partial D^k}{\int }\left(\frac{{\widehat{n}}_x{Y}^{+}\left[E_h^y\right]-{\widehat{n}}_y{Y}^{+}\left[E_h^x\right]-\alpha \left[H_h^z\right]}{2\left\{\left\{Y\right\}\right\}}\right)}l(r)dr\right).\end{array}$ (15)

The semidiscrete expression are ordinary differential equations (ODE) in the time-domain (TD) which can be solved by using five-stage 4th order low-storage explicit Runge-Kutta (LSERK) method, given in.¹⁹ To ensure the stable time steps while integrating the solution in TD, Courant-Friedrichs-Levy (CFL) given in²⁰ is used. This study uses unstructured triangular mesh elements. The marching step (Δt) as per CFL condition for this type of mesh is given in Eq. (16):

$\Delta t\le CFL(\frac{2}{3}\min \Delta r_i){\min }_{\Omega }(r_D)$ (16)

where rD is the normalized radius of the inscribed circle in a given triangular element and Δri is the grid spacing for one-dimensional (1D) standard interval of.11 For 4th order LSERK method, choose CFL=1. Applying PEC boundary: In the DGTD method, the perfect electric conductor (PEC) boundary is very easy to implement. The tangential components are set to zero i.e., $\widehat{n}\times E=0$ only by changing the values at the boundary interfaces. The values of all nodes appearing at the boundary are changed by using mirror principle given in,¹⁸ i.e., ${(E^z)}^{*}=-{(E^z)}^{-}$ such that ${(E^z)}^{-}+{(E^z)}^{*}=0.$

This is implemented as given in Eq. (17).

$\left[E^z\right]=\widehat{n}\cdot 〚E^z〛=2{(E^z)}^{-},$ (17)

Exact method

The cutoff frequency fc of any waveguide having dimensions a × b is given in Eq. (18). The permittivity ε and permeability μ are the material parameters used for filling the waveguide either be air or vacuum. This expression holds true for either TM^z or TE^z mode.²¹

${(f_c)}_{mn}=\frac{1}{2\pi \sqrt{\mu \epsilon }}\sqrt{{\left(\frac{m\pi }{a}\right)}^2+{\left(\frac{n\pi }{b}\right)}^2}$ (18)

Simulation setup

The waveguide WR90 having dimensions a=0.9 in. (2.286 cm), b=0.4 in. (1.016 cm) is selected for the analysis of modes in the numerical method. The coarse triangular mesh is obtained from commercial mesh-generating software. The initial conditions used for TMz case are given in Eq. (19). For the TEz mode, the conditions are similar. So we do not repeat them.

$\begin{array}{l}{E}_z=\sin (m\pi fx).\sin (n\pi fy);\\ Hx=0;Hy=0;\end{array}$ (19)

where

1. m and n are modal values;
2. f is frequency;

The x and y are N×K matrices of the nodes along x- and y-axis, respectively, in 2D.

Here, N is the order of polynomial interpolation used in each Kth element. All the units are made non dimensional by using set of expressions given in Eq. (20).

$\frac{x}{L}\to \widehat{x},\frac{y}{L}\to \widehat{y},\frac{ct}{L}\to \widehat{t},Z_{0}H\to \widehat{H},E\to \widehat{E}$ (20)

where

1. c is speed of light, t is time;
2. $Z_0=\sqrt{{\mu }_0/{\epsilon }_0}$ is the free-space impedance;
3. L is reference length for a given problem.

In this study, we choose L=1.0160 cm which is the smallest dimension of the waveguide. Hence, all the parameters with hat (^) are non dimensional. The time-domain data for electric field is recorded at some test points in the given computational domain i.e., [–1.1430 1.1430]×[–0.5080 0.5080] and their frequency response is obtained by employing fast Fourier transformation (FFT) on the recorded data. Since the unstructured triangular mesh suffers from accurate positioning of the test point in the Cartesian coordinate system, this issue is resolved by using adaptive mesh method of node displacement mentioned in.¹⁵

Effect of resolution of finite elements

This section discusses the effect of resolution of number of non overlapping elements (K) used to discretize the computational domain. The simulation is performed and data is recorded at test point located at (x₀, y₀)=(0, 0). The effect of different modes appearing with respect to K is summarized in (Figure 1).

Figure 1 Effect of K on resonating modes with N=4: (a) K=30, (b) K=68, (c) K=134, and (d) K=546.

Figure 1 shows that as the resolution of the FEs used to discretized computational domain is increased from K=30 to K=546, the resonant modes are shifted from lower frequency modes to higher frequency modes. The reason is: as the resolution increases, the time step required for stable marching also decreases. The decreased value of Δt supported higher frequency components when FFT is used. So, if we are interested in lower frequency modes of a waveguide structure, we must use lower resolution of FEs in the computational domain. In the next section, we discuss the effect of Final Time on these models.

Effect of final time

This section explains the significance of using different final time for estimating the resonant frequency of the waveguide. The findings are shown in (Figure 2) where data for two values of Final Time=[200,500] is recorded. (Figure 2(a-b)) shows the time-domain data of electric field recorded at test point (x0, y0)=(0, 0) for the two final times and (Figure 2(c-d)) presents the spectrum of the recorded data which is obtained by taking FFT of the time-domain data. From (Figure 2), it is obvious that we need enough time to obtain different resonant modes. But, the longer time only captures more energy at the already existing modes, as can be seen from the variation in the peaks of normalized amplitudes in (Figure 2(c-d)). By increasing the Final Time we cannot obtain more modes. The appearance of the modes depended on the resolution of FEs used, as has been discussed in the previous section.

Figure 2 Effect of final time on resonating modes for K=68, N=4 and Δt=0.01427: (a) K=30, (b) K=68, N=4, (c) Spectrum View at normalized Final Time=200, and (d) Spectrum View at normalized Final Time=500. (Dimensional Final Time=Final Time×L×Δt /c).

Comparison of analytical and numerical results This section summarizes the results of both TMz mode and TEz mode. The exact results are calculated by using the analytical formulation and the numerical results are obtained by using the DGTD method with appropriate settings K. The modes calculated by using TM polarization are presented in (Table 1). The effect of variation with K is also given. Similarly, the values obtained by using TE polarization are given in (Table 2). The lower frequency modes in the TE polarization are obtained by using K=30.

In this paper, modeling of waveguide is performed by using the DGTD method in 2D by incorporating the both TM and TE polarization. The space discretization of the computational domain by changing resolution of finite elements (FEs) i.e., K has effect on stable time marching step (Δt). This resulted in the finding that higher the value of K, hence leading to smaller time step, supports higher frequency modes. In order to obtain the lower frequency modes of waveguides, we need to reduce the resolution of FEs. Secondly, increasing the Final Time too much does not provide all the frequency modes rather energy in the existing modes keeps increasing. In order to get all the low- and high-frequency modes, we needed to perform multiple simulations with different low- and high-resolution of the FEs in the computational domain.

This work has been partially funded by Postgraduate Technology Innovation Project of Jiangsu Province (No: KYLX15_0286), Natural Science foundation of Jiangsu Province (No: BK20140810), and Natural Science Foundation of China (No: 61401199).

The author declares no conflict of interest.

1. Collin RE. Field Theory of Guided Waves. Electromagn Waves. 2nd ed. USA: Wiley Press; 1990. p. 1–864.
2. Yan S, Greenwood AD, Jin JM. Modeling of Plasma Formation During High–Power Microwave Breakdown in Air Using the Discontinuous Galerkin Time–Domain Method. IEEE J Multiscale Multiphysics Comput Tech. 2016;1:2–13.
3. Li P, Jiang LJ, Shi Y, et al. A discontinuous galerkin time domain–boundary integral method for analyzing transient electromagnetic scattering. IEEE Antennas Propag Soc Int Symp. 2014:1893–1894.
4. Gedney SD, Roden JA, Crawford RD, et al. A discontinuous galerkin finite element time domain method with PML. IEEE Antennas Propag Soc Int Symp. 2008:1–4.
5. Fahs H. High–order discontinuous Galerkin method for time–domain electromagnetics on non–conforming hybrid meshes. Math Comput Simul. 2015;107:134–156.
6. Gao S, Cao Q. Discontinuous Galerkin time domain method for SOI thin–ridge waveguide problem. Trans Nanjing Univ Aeronaut Astronaut. 2013;30:162–168.
7. Li P, Shi Y, Jiang LJ, et al. DGTD analysis of electromagnetic scattering from penetrable conductive objects with IBC. IEEE Trans Antennas Propag. 2015;63:5686–5697.
8. Cfdrc VK. A Discontinuous Galerkin Time–Domain Method with Dynamically Adaptive Cartesian Meshes for Computational Electromagnetics. IEEE Trans Antennas Propag. 2016;65(6):1–18.
9. Chen J, Liu QH. Discontinuous Galerkin Time–Domain Methods for Multiscale Electromagnetic Simulations: A Review. Proc IEEE. 2013;101:242–254.
10. Reed WH, Hill TR. Triangular Mesh Methods for the Neutron Transport Equation. Los Alamos Sci Lab. 1973;LA–UR–73–4:1–23.
11. Cockburn B, Shu C–W. Runge–Kutta Discontinuous Galerkin Methods for Convection–Dominated Problems. J Sci Comput. 2001;16(3):173–261.
12. LeVeque RJ. Finite volume methods for hyperbolic problems. India: Cambridge University Press; 2002.
13. Hesthaven JS, Warburton T. Nodal Discontinuous Galerkin Methods, Algorithms, Analysis and Applications, Springer, UK; 2008.
14. Hussain I, Li H, Wang Y, et al. Modeling of excitation source for time–domain EM solvers. In: Prog Electromagn Res Symp St Petersburg, Russia; 2017.
15. Hussain I, Li H, Cao Q. Multiscale Structure Simulation Using Adaptive Mesh using DGTD Method. IEEE J Multiscale Multiphysics Comput Tech. 2017;2:115–123.
16. Le Veque RJ. Finite Volume Methods for Hyperbolic Problems. India: Cambridge Univ Press; 2002. 54:258.
17. Toro EF. Riemann solvers and numerical methods for fluid dynamics, Springer, UK; 2009.
18. Hesthaven JS, Tim Warburton, Warburton T. Nodal Discontinuous Galerkin Methods. Algorithms, Analysis, and Applications. Editors. Springer, New York; 2008.
19. Niegemann J, Busch K, Singh MR, et al. Time–Stepping and Convergence Characteristics of the Discontinuous Galerkin Time–Domain Approach for the Maxwell Equations. AIP Conf Proc AIP. 2009;1147(1):22–29.
20. Toulorge T, Desmet W. CFL Conditions for Runge–Kutta discontinuous Galerkin methods on triangular grids. J Comput Phys. 2011;230(12):4657–4678.
21. Balanis CA. Advanced Engineering Electromagnetics. 2nd ed. USA: John Wiley & Sons; 2012.