The Dilemma of Resolving the Low-Frequency Breakdown Problem in Microwave Components via Traditional and Improved Finite-Element Time-Domain Techniques

This paper developed an improved numerical methods for calculating electromagnetic field at low frequencies, and compared the performance with the traditional method. Traditional finite-element time-domain (FETD) is known to suffer from low frequency breakdown (LFB) problem, where system matrix becomes ill-conditioned, and leads to instability of small size electrical structures. In this paper, we proposed a mixed FETD (mFETD) method in resolving the LFB anomaly in the traditional FEM, with a particular reference to transmission line, inductor, and coaxial cable. We considered wave equation of E-field with incorporation of divergence constraint equation (DCE) as a function of Lagrange multiplier using current continuity equation (CCE) and Gauss’s law. For spatial discretization, both nodal basis functions and Curl conforming vector basis functions were selected, and we employed implicit Newmark beta algorithm (NBA) for integration of time. We describe how the components constructed via DCE in system matrix mitigate the singularity impact of the stiffness matrix, which consequently led to significant improvement of the system matrix. Numerical experiment and results demonstrate how the mFETD obtains stable numerical solution and attains faster rate of convergence while using an iterative solver, as such, improves the computational efficiency. Therefore, the mFETD method handles the transient problems in transmission line, which causes LFB anomaly in the traditional FETD, but not efficient in terms of computational time and iteration at solving the coaxial cable problem.


I. INTRODUCTION
Transmission line segments finds wide application in the design of circuits, RF, and microwave components, such as power combiners/dividers, phase shifters, delay lines, and filters. Characteristic impedance and electrical length are the major parameters that determine the behavior of transmission line (TL) segment [1]- [5]. At microwave frequencies, the electrical length needed for the implementation of the above listed components can be achieved by TL of few centimeters length at most. At lower frequencies, the design of circuits and components based on TL segments is not achievable in practice except for a compaction procedure, which usually arrive at a practical structure with no reduction in electrical length. Meandering is the popular approach to The associate editor coordinating the review of this manuscript and approving it for publication was Wen-Sheng Zhao . achieve such goal in planar technologies. In such situation, a straight, long TL is transformed into a set of short TL connected by meanders. There is reduction in the physical length of but with greater component width for the same electrical length. On the other hands, all meander along the path of signal causes a discontinuity that require consideration during the modeling of TL segment behavior [6]. In the simulation of integrated circuits, radiation antenna, signal integrity, etc. numerical methods, such as, finite-difference time-domain (FDTD), method of moments, spectral-element method (SEM), and finite-element method (FEM) have been widely used. However, all the aforementioned traditional fullwave solvers suffer from ''low frequency breakdown (LFB) problem'' [7].
The FETD method has been demonstrated as a powerful and versatile tool for analysis of various electromagnetic devices. However, while solving low frequency transmission 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/ line problem in which the size of the transmission line in finite-element discretization is very small compared to the wavelength, there is a problem of instability in FETD, which is referred to as ''(LFB)'' in finite-element method or ''non-physical dc modes'' phenomenon [7]- [11]. LFB in electromagnetics is related to spurious modes for eigenvalue problems. It leads to slow convergence in an iteration or unstable numerical solution of a direct solver for radiation phenomena in both time and frequency domains. The reason for this LFB is the irrotational vector space spanned by the edge basis functions leading to non-zero divergence of electric flux density. That is, the divergence of E-field is not constrained by Gauss's law in the traditional methods.
Overcoming the LFB phenomenon, Gauss's law has to be considered as discussed in [8]- [13]. These works rely on orthogonality solution of gradient matrix (null space), where free DCE is introduced using the tree cotree partitioning or projection operator of the mesh of FE to eliminate spurious modes for frequency domain eigenvalue problems. An analytical method is developed in [14] via decomposition of the solution into its complementary and null space term with no computational complexity to solve the LFB problem in FEM. Kikuchi [15] proposed mixed FEM (mFEM) method enforcing the constraint of free divergence for eigenvalue issues, where free divergence condition is considered the equation of constraint and added into the system equation as a function of Lagrange multiplier. Mixed SEM is another method [15]- [17] that is appropriate for large scale problems with higher convergence order. Furthermore, in [18], a mixed FEM/SEM is used to solve the LFB phenomenon in sub-surface sensing problems with low-frequency electromagnetic fields. mFEM has been formulated in frequency domain in order to solve the LFB anomaly. This paper extends the mFEM in time domain called mixed finite-element time-domain (mFETD) to enhance matrix conditioning and numerical stability when solving very small sizes transmission line problems. Of cause, this may be extended to other microwave components. We focus on transmission line because it is a fundamental component in antenna or field-circuit coupled analysis. Specifically, we considered wave equation of E-field with incorporation of divergence constraint equation (DCE) as a function of Lagrange multiplier using current continuity equation (CCE) and Gauss's law. For spatial discretization, both nodal basis functions and Curl conforming vector basis functions were selected, and we employed implicit Newmark beta algorithm (NBA) in order to integrate time. We describe how the components constructed via DCE in system matrix mitigate the singularity impact of the stiffness matrix, which consequently led to significant improvement of the system matrix. In the end, the proposed mFETD method is able to handle the transient problems in transmission line, which causes LFB anomaly in the traditional FETD. Finally, it is important to clearly state that Cheng et al. [19], first looked into and developed mFETD method. The novelty of our work as against/compared with [19], [20] is that this paper demonstrates why the mFETD is not efficient at solving the problem of coaxial cable, under the circumstance of higher computation time and iteration as traditional FETD method. This is the major gap filled in our work.
The rest of this paper is organized as follows. In Section II, we formulated the problem, showing the LFB in the traditional FETD method. The proposed solution of mFETD method is presented in Section III. Section IV presents numerical experiment, which verifies how the proposed mFETD method is more advantageous than the traditional FETD. Conclusion is drawn in Section V.

II. PROBLEM FORMULATION
Consider In traditional FETD, the wave equation (second order) for E-field in domain is given as follows: which subjects to the outer boundary conditions with any combination of perfect magnetic conductor (PMC) boundary, PMC , a perfect electric conductor (PEC) boundary PEC , or/and the first order absorbing boundary condition ABC , where µ and σ e represent the permeability, permittivity, and electric conductivity of the material, respectively, B i and J i denote the magnetic and electric current densities, respectively, whilen denotes the unit normal of the outward boundary [8].
The variation formula for the initial boundary value problem can be gotten by finding E ∈ H (curl; ) in such a way where H (curl; ) = w ∈ L 2 ( ) 3 : ∇× w ∈ L 2 ( ) 3 w (r) represents the vector basis functions.
Combining the Galerkin technique with mixed order of the curl-conforming vector basis functions within H (curl; ) for the expansion of E-field, such that where M, D and S depict mass, damping matrices, and stiffness, respectively. It is important to know that the ABC is incorporated into the damping matrix D. e represents the column vector holding the unknown parameters of the E-field in domain. m and j denote column vector associated with enforced magnetic and electric current density, respectively, following the basis functions test.  In traditional FETD method, adoption of NBA for which β = 1 4 towards the integration of time gives rise to a resultant time updating system matrix G that combines M, S, and D t is the increment in time step. In Eqn. (7), S is singular while M is well conditioned. The generated weight by stiffness matrix in Eqn. (7) affects the G matrix, which consequently decides the accuracy and convergence of the solution.
In this paper, the spatial sampling density (SD x ) and temporal sampling density (SD t ) are defined as a function of PPP (points per period) and PPW (points per wavelength) at maximum frequency selected for simulations. They were employed in the analysis demonstrating the existence of the LFB in the traditional FETD. Furthermore, for any given transient relationship, SD x and SD t can be calculated using where x and λ min are the typical element size and minimum wavelength, respectively. f max represents the maximum frequency of choice. According to literature [16], S increases based on (1 x) and M grows with x. The ratio of S to M approaches c 2 x 2 , where c is the speed of light. Considering the relationship in Eqn. (8), it can be seen that t 2 S and (SD 2 x SD 2 t )M share equal order of magnitude. Then, as SD 2 x SD 2 t get bigger, the system matrix gets more ill-conditioned as a result of the increase in weight of the fitness matrix. Applying the stable unconditional NBA time integration technique, t is arbitrarily big with no impact on the stability; as such, SD t is set to a small value (e.g. , in order to ensure the requirement for accuracy. Thereby, SD x majorly determines the behavior of system matrix. The traditional FETD method can provide appropriate results as long as SD x is not too big. That is, the effect of M is not ignored by the computer. On the other hands, when the simulated transmission lines or frequency of operation is low enough, SD x arrives at a reasonable value. It then makes G nearly singular, and breaks down the solution. Therefore, SD x is considered as the condition for the LFB phenomenon in transmission lines in an inherent manner.

III. MIXED VARIATIONAL FORMULATION FOR THE PROPOSED mFETD
In order to overcome the LFB phenomenon, we combine Gauss's law with the CCE to formulate the constraint model as Thus, Eqn. (9) is substituted into Eqn. (5) by Lagrange multiplier τ (r, t) as given in [18] within the framework of finiteelement. The variational model for mFETD is to compute E ∈ H (curl; ) and τ ∈ H 1 ( ) in such a way that = − w (r) , ∀w ∈ H(curl; )    As against the frequency domain, here Lagrange multiplier τ (r, t) depends on space and time. With separation of spatial and temporal dependencies, τ (r, t) is expressed in form: where M, D, S, j, m, and e are as defined in Eqn. (6). G 1 and G 2 denote the integrations based on Lagrange multiplier τ ,  while T represents the transpose operator of matrix. y is a column vector that represents the unknowns of τ at the entire nodes, apart from the PEC nodes of the loop. j 1 is a column vector initiated by Eqn. (11), and it is a term that denotes the electric current source. The expression of each parameter is as presented in Table 1.
Employing NBA of β = 1 4 for integration of time, the last matrix G for mFETD of such time updating equation is The system matrix G in Eqn. (14) ensures the sparsity and symmetry of the FEM. Making e and y in Eqn. (13) of the same order of magnitude, it becomes important to multiply G by a scaling factor of constant η. η is selected as a ratio of the largest element in the (M + 1 2 tD + 1 4 t 2 S) matrix to maximum element in (G 1 + 1 2 tG 2 ) matrix. Here, G is now Moreover, a simple preconditioner called diagonal scaling concept is incorporated into the mFETD for better performance. Simple comparison between Eqn. (15) and Eqn. (7), we observe that the added terms in Eqn. (15) which is imposed by the divergence constraint Eqn. (11) remove the impact of stiffness matrix S singularity nature on our system matrix G. Hence, mFETD basically enhances the system matrix and eliminates the LFB phenomenon in transmission lines, which demand large spatial SD x

IV. NUMERICAL SIMULATIONS AND RESULTS
In order to show and compare the performance of the proposed mFETD method against the traditional FETD method while overcoming LFB, a numerical simulation example is presentedin this section in proof-of-concept. The simulation is conducted with Matlab R2018a installed on personal HP PAVILION laptop computer (1-TB memory, 16G RAM, and Intel core i7-8565U CPU @1.80GHz 1.99 GHz).

A. COAXIAL CABLE
For illustration, we consider a coaxial cable filled with air and of 1m length as shown in Figure 1. The coaxial cable is made up of inner and outer perfect electric conductors of radius 5 and 10 mm, respectively. We excited the cable using a lumped port loaded with 10 −3 inductor (H) at near end, and a current source, parallel to a resistor, 7-, at the far end. Furthermore, the lump port has its first and the second side positioned in the inner and outer PEC, respectively. Gaussian pulse is used to define the time function of the source, which VOLUME 10, 2022 is expressed as (16) with γ = 0.484/f max , where f max denotes highest frequency of 200 kHz. For this scenario, it becomes easy to formulate transmission line model (TLM), which is used to validate mFETD. In order to emulate the TLM so well, we enforce a PMC boundary at the two ends in 3D FE model to generate an open circuited plane [18]- [23].
Firstly, an LU decomposition from a direct solver generated from Matlab is adapted for the evaluation of the resulting matrix system in the mFETD and traditional FETD. The input impedance and the transient voltage at the excited end are plotted for clear comparison between TLM, mFETD, and the traditional FETD in Figure 2 and Figure 3, respectively. It is clearly observed that the three methods agree quite well. The computed relative error in the transient voltage executed by TLM and mFETD is nearly 0.51%.
To evaluate and justify mFETD performance at low frequencies, we adopt a biconjugate gradients approach and the insufficient LU decomposition preconditioner, which is an iterative solver generated by Matlab in the traditional FETD and mFETD. The condition number and convergence of the resultant system matrix with f max changing from 3 GHz to 30kHz are examined. It is seen and observed in Table 2 how SD x proportionally increases as f max decreases. Also, when SD x is increased to 1000 PPW, the model under consideration can be considered electrically small. SD t is set at 20 all the time due to the unconditionally stable time stepping framework. When SD x is less than 1000 PPW, the system matrix of mFETD and that of the traditional FETD exhibit appreciable properties, showing small condition number and small average number of iterations for each time-step to achieve 10 −6 convergence tolerance. In other words, the condition number of our proposed mFETD is by far smaller than the traditional FETD counterpart, for such problem. The traditional FETD require big average number of iterations for each time step before convergence or non-convergence of the iterative solver. However, in mFETD method, few iterations are needed for convergence to exact numerical solution even when SD x approaches 10 6 at 30 kHz f max . Note that the constraint equation adds some level of nodal degree of freedom (DoF), and the DoF of mFETD is biggerthan the traditional FETD method. The condition numbers and average number iterations for each time-step against f max is as intuitively depicted in Figure 4 and Figure 5, respectively.
Finally, it is important to report the computational resource of the methods. The CPU time of the traditional FETD method is 21 minutes while that of mFETD method is 27 minutes. mFETD requires 23 GHz memory space while the traditional FETD requires the 21 GHz memory space. These facts show the bottleneck of the proposed method. mFETD methodexhibits higher CPU time and memory space than the traditional FETD because it has larger number of parameters.

B. AN INDUCTOR MODEL
We further show the application of the proposed mFETD method by designing a spiral circular inductor with 2.9 × 10 −7 H inductance value as depicted in Figure 6. The inductance value shares a relationship with the overall width of the coil [24]. Assuming the number of turns of the spiral circular inductor equals 16 and we vary w to get the wanted inductance value. The port 1 is the origin or source that is made up of 100 MHz highest frequency Blackman Harris window (BHW) function, while port 2 serves as the destination or receiver. To achieve impedance matching, the two ports have internal resistances R s = 50 . At maximum frequency, there is sampling density of 5500 points per wavelength (PPW). In general, whenever the sampling density is higher than 10 3 PPW, it is considered electrically small model with low frequency. The FETD system matrix property get worse and easy to breakdown as sampling density becomes larger; on the other hand, the specific sampling density as a function of PPW that makes the conventional FETD technique to breakdown depends on certain condition. We applied an ABC in this model and air box with 25mm×25 mm×2 mm dimension. The time-domain simulation results obtained are presented in Figure 7. In time domain, there is a good agreement between the traditional FETD in COMSOL and the proposed mFETD methods. The relative voltages errors at port 1 and port 2 are 3.38 % and 4.05 %, respectively. In order to estimate the S-parameters, we can use Fourier transform technique to transform the transient results into the frequency domain. Hence, the value of inductance L can be computed using when the value of w is changed, it is noted that when w is 50µm, the value of the inductance is 2.9×10 −7 H . Therefore, w is set at 50µm.
For a reduced maximum frequency of Blackman Harris window pulse of about 1 kHz, the sampling density becomes 165000 PPW. This means the representation of electrically small model at low frequency. R s = 0.0005 is chosen to ensure impedance matching. Figure 8 depicts the results of the traditional FETD and mFETD techniques. In this category, investigating inductance value of the inductor between 100 Hz and 100 MHz is the goal; however, there is breakdown at low frequencies in the traditional FETD method (i.e. COMSOL). Hence, at higher frequencies, the computed inductance value via the traditional FETD is considered as a reference and compared against mFETD technique, and as shown in Figure 9 (a), the relative error is 2.33%. At less than 1 kHz frequency, the COMSOL breaks down, consequently unable to generate good results. Conversely, the relative error that exist between L value computed by mFETD at high and low frequency results fall in the confines of 3.55 %, based on the result presented in Figure 9 (b). Therefore, there is a constant inductance value.  The computational resources of mFETD and COMSOL at 100 MHz maximum frequency of BHW is investigated. The COMSOL requires 1177034 DoFs, CPU time of 67 min. and 37 GB memory, while mFETD requires 1478313, CPU time of 54 min. and 49 GB memory. We can observe that mFETD technique requires more degree of freedoms, even so, it requires reduced time compared to the COMSOL. This is because, at low frequencies, the mFETD matrix properties is good. At very low frequencies, the iteration from the COM-SOL solution does not converge due to poor system matrix condition.
Finally, we compared the computational parameters for the proposed mFETD and FETD-based method in [17]. Running the two algorithms under the same conditions of the Inductor model Section IV B. The result is as presented in Table 3. As can be seen, more DoF in MFETD [17] method causes a longer computation time than the proposed mFETD method. Also, it takes lesser time to compute the gradient matrix than the total time in MFETD [17] method, which implies the proposed mFETD method, obtains well-conditioned system matrices with little additional computation in overcoming the low-frequency breakdown.

V. CONCLUSION
In conclusion, we presented in this paper, a new mFETD method capable of resolving the LFB anomaly that cannot be handled by the traditional FEM method, while solving transient problems in transmission line. We incorporated the DCE into wave equation for E-field as a function of Lagrange multiplier using CCE and Gauss's law. The nodal basis functions and Curl conforming vector basis functions were selected in order to achieve spatial discretization, also, we employed implicit NBA technique for integration of time. The numerical experiment show that the proposed mFETD is capable of overcoming the LFB anomaly in transmission lines, which is difficult for the traditional FETD. With improvement in the system matrix properties, mFETD show high level of stability from dc level to high frequency bands, which evidently implies that mFETD is more adequate and effective alternative approach in modeling transmission lines. The proposed mFETD method can be extended to microwave circuits and electromagnetic scattering problems.
Finally, it is worth mentioning that the proposed mFETD is not effective in the coaxial cable example, according to the computational time and iteration as traditional FETD method.