Numerical Stability of Dual Full-Wave Formulations With Electric Circuit Element Boundary Conditions

We have successfully implemented full-wave (FW) frequency-domain <inline-formula> <tex-math notation="LaTeX">$E$ </tex-math></inline-formula>- and <inline-formula> <tex-math notation="LaTeX">$H$ </tex-math></inline-formula>-based finite-element (FE) formulations with electric circuit element (ECE) boundary conditions (BCs) for high frequency (HF) applications. These formulations are stable and accurate at HF, with an unknown field strictly inside the domain and a scalar potential with support only on the boundary. Aiming at having all-purpose Maxwell’s solvers at hand, we investigate their numerical stability at very low frequencies (LFs). We show that these FW formulations can be used at LF provided that techniques to ensure numerical stability are applied. If only global (energy-based) stability is sought, then just using the excitation which is essential for the finite element method (FEM) might be enough. If a local (field spatial distribution) stability is needed, we show that scaled hybrid versions of these formulations ensure LF stability. Models of increasing complexity are considered: a homogeneous cylinder (2-D axi/3-D), a coaxial cable (2-D axi/3-D), a coplanar waveguide (3-D).


I. INTRODUCTION
T HE modeling of coupled field-circuit problems is still a topic of interest, see e.g., [1] for quasi-statics (QSs).A natural coupling, appropriate not only for QS but also for the full wave (FW) electromagnetic (EM) regime is obtained when the field part of the model uses a formulation with special boundary conditions (BCs) called electric circuit element (ECE).For FW, a theory for weak E-/H-based formulations for linear materials and ECE BC can be found in [2].Our approach differs from other relevant works e.g., [3], where circuit models for the EM domain are coupled with external circuits in circuit simulators.The use of ECE allows an easy field-circuit coupling either in the field simulator, thus not losing the spatial information and providing direct access for postprocessing fields, or in a circuit simulator after model reduction and circuit synthesis, as in [4].
This article is motivated by the need for robust formulations for large frequency ranges [5], [6], [7].Several approaches are proposed in the literature to improve the numerical stability of the finite element method (FEM) solutions for the FW Maxwell equations in the frequency domain (FD).One possibility is to use so-called hybrid field formulations involving two or more fields as main unknowns as in [8] with (E, B) or in [9] with (D, H).An E, V mixed approach is proposed in [5] where a dummy additional scalar field is used to stabilize the numerical solution, which eventually vanishes in the computational domain.In [6] a monolithic formulation is proposed, without auxiliary variables, based on function spaces splitting.Scaling equations is another option, as in [10] for electro-QS and models with perfect insulators.The conditioning number of the assembled system constitutes an experimental indicator for stability [5], [6], [10].
Previously we have successfully implemented ECE BC in FEM for dual FW FD: a formulation in E (FWE) with a scalar potential V on the boundary in [4] and a H-based one (FWH) with a scalar potential ϕ on the boundary in [11].Here we investigate their numerical behavior at low frequeny (LF).

II. E -AND H -BASED FORMULATIONS WITH ECE BC
Let be a simply connected domain, linear and isotropic, 1with a Lipschitz boundary ∂ that includes m disjoint terminals S k , k = 1, . . ., m, either voltage (k ∈ I v ) or current excited (k ∈ I c ).

A. Weak E Formulation With V on the Boundary
The functional spaces are ∈ H V,0 .We implement (1) with these spaces. 2he coupling of unknowns happens at discrete level where Neint-number of edges inside the domain; NE-number of nodes on the boundary but not on the terminals; NnTermK-number of boundary nodes placed on the terminal k; ϕ j -first order nodal elements on the boundary; N j -first order3 edge elements inside the domain (which are curl-conform).Symmetry conditions E t = 0 is essential, H t = 0 is natural.The same FWE formulation holds in 2-D and 3-D.

B. Weak H Formulation With ϕ on the Boundary
The functional spaces are: The coupling between the formulation inside and at the boundary is done at the discrete function space level where Ne are the edges inside the domain and on the surface of terminals; NH are the nodes on the boundary but not on the terminals; NeCutTermK are the edges that belong to the cut that corresponds to the terminal k.The cuts are automatically generated with the cohomology solver in gmsh.Symmetry condition H t = 0 is essential, E t = 0 is natural.Different formulations 4 are needed in 2-D versus 3-D.In 2-D FWH does not need for ϕ, because the boundary nodes either belong to a cut, or to a classical surface.The inner edges are perpendicular to the domain and thus the unknowns are associated with nodes.

III. HYBRIDIZATION OF FW ECE BC FORMULATIONS
Note that, strictly speaking, the above formulations are also hybrid, as we have potential scalar unknowns (V /ϕ) with support limited to the boundary and fields inside the domain.In order to ensure stability at low frequency, further hybridization of formulations (1) and ( 3) is needed.
A first approach consists of extending the support of the existing scalar fields V /ϕ (number of unknowns NE/NH) to the inside of the domain, together with the fields E/H.
A second hybridization possibility is using a dummy scalar field inspired by [5].The uniqueness theorem for the FW formulation in FD requires that the material properties are strictly positive and frequency is non-zero.Therefore, the hybridization of an FWE-based formulation proposed in [5] assumes non-zero conductivity everywhere.However, even for a homogeneous domain with non-zero conductivity the FW formulation is unstable at LF due to the non-uniqueness of the FW field solution at a zero frequency, caused by the non-trivial kernel of the curl operator.In the case of domains with high parameter contrast, this LF instability worsens.Inspired by [5], we have modified the initial formulation, by introducing a dummy scalar field.A system of two weak equations is solved, with unknowns the field inside the domain, the scalar potential on the boundary, and a new dummy scalar potential both inside the domain and on the boundary.The discrete field is still given by ( 2)/(4).The discrete expression of the dummy field V d is obtained only with nodal elements.In the hybrid EV the curl-curl ( 5) is enhanced with a new term involving the dummy field as unknown.Equation ( 6) is obtained by explicitly imposing the continuity law, and projecting it on test functions from the dummy field function space.The hybrid Hϕ is analogously obtained, with magnetic Gauss's law explicitly imposed.Scaling with a factor p is easily implemented, and does not affect the field IV. BENCHMARKS AND HF VALIDATION Five test cases of increasing complexity are considered.The first two are 2-D axi/3-D models of the conducting homogenous cylinder described in [4], single input single output (SISO) models, with two terminals (one grounded).It reveals the inherent LF instabilities FWE and FWH.Then, we consider 2-D axi/3-D models of an air-filled coaxial cable as in [8].The coaxial is a multiple input multiple output (MIMO) model, with three terminals (one grounded), two homogeneous subdomains, one highly conductive and the other perfect insulator.Higher instabilities are expected due to the zero/low conductivity of the air.Transmission lines (TLs) theory provides an exact analytical solution−the reference.The validation5 of this test case is shown in Fig. 1.
The final test is a coplanar waveguide (CPW) with an analytical but approximate solution based on TL theory (for per unit length parameters) taken as reference.The frequency range spans from LF to high frequency (HF), including capacitive and inductive effects even if propagation is not apparent.The

V. LF TO HF NUMERICAL STABILITY
The LF results are investigated both locally (field distributions in the domain) as well as globally (terminal quantities: admittances and impedances in case of voltage and current excitation, respectively).The terminal quantities are linked to the power and stored energy since the complex power transferred through the system terminal is Note that hereafter we adopt MUMPS solver with default settings in PETSc, i.e., with reordering and scaling.

A. Quick Hybridization With the Existing Scalar Field
This approach has a stabilization effect at LF.This is illustrated in Figs. 3 and 4 for the coaxial cable and the CPW, respectively, showing FWE results at 1 Hz with voltage excitation.The field E is unstable with V only at the boundary, but becomes stable with V inside.

B. Chose an Appropriate Excitation Type
The global solution at the terminals can sometimes be stabilized by choosing an appropriate excitation, even when the   scalar potential is only defined on the boundary.In the MIMO examples, we considered both terminals excited in the same way, voltage or current. 6Fig. 5 shows results for the 3-D coax test.The worst results are obtained for FWH voltage excitation (bottom left) and FWE current excitation (Fig. 5 top right), the terminal solution being stable only over 100 kHz.Note that the result at 1 Hz, for FWE voltage excitation is very accurate also when using V only on the boundary, despite the instability noticed in the field (Fig. 3 top).For FWE voltage excitation (Fig. 5 top left), instabilities occur for very LF (less than 1e−5 Hz, much lower than a practical limit).No instabilities are seen at the terminals for FWH with current excitation (Fig. 5

bottom right).
For the CPW test, instabilities are not as dramatic as in the coaxial case, but the behavior is similar (Fig. 6).Using the excitation that corresponds to essential BC is better.Recall that, in FEM, essential BC are imposed by constraining the values of the corresponding DoFs, whereas natural BC are only satisfied on average by the final numerical solution.

C. Use a Non-Zero Conductivity in Perfect Insulators
As suggested by the hypotheses of the uniqueness theorem for a well-posed FW-EM problem, the numerical instability   can be solved by means of a strictly positive value for the conductivity in the air (Fig. 7 left).This value can be correlated with the value of the dielectric constant for the maximum frequency for which the instability happens (here 10 kHz).This is illustrated by the regime map (represented as explained in [12]) of a point in air (Fig. 7 right).To stabilize the solution, σ air must be chosen so that the maximum frequency of the instability zone enters the electrostatic regime (ES): τ | max( f unstable ) = 1/(max( f unstable )) ≈ 10τ e = 10ε/σ air .

D. Hybridization With a Dummy Scalar Field
When using solvers from well-known libraries, one must pay attention to the options used.Here, the default setting means some inner scaling of the solved system.However, a better conditioned assembled system is expected to produce a more stable solution [5], [10].Figs. 8 and 9 show the conditioning numbers of the assembled systems for the implemented formulations.The legends also indicate the size of the assembled system.No difference is observed at the level of conditioning number for different excitations.

VI. CONCLUSION
At HF (QS to FW), the FW FEM formulations with ECE BC and scalar potential only on the boundary are stable and robust even if air has zero conductivity.In 3-D the computational effort of FWE and FWH are similar, but in 2-D FWH is much cheaper.In 3-D FWE leads to better conditioned systems, but FWH results are also accurate.In 2-D FWH is better conditioned.At LF, for 2-D models FWH is the best option, and to stabilize it is enough to choose a small non-zero conductivity for perfect insulators.As to the 3-D models, the best option is the scaled hybrid version with a dummy field of FWE, which is able to ensure almost constant conditioning numbers.This formulation must be applied only after setting a small non-zero conductivity for perfect insulators.The optimum scaling can be obtained by 1-D optimization for a selected LF.
To use FWE/H toward LF, we recommend: 1) start with boundary scalar potential and use essential excitation; you may use σ air = 0; 2) use a quick hybridization, check fields; and 3) set σ air ̸ = 0 and use scaled hybrid formulations.However, an FW solver used in LF to extract static behavior can never outperform a dedicated solver.We advocate the use of ECE BC especially for HF applications, including radiation [13].

Fig. 3 .
Fig. 3. Coaxial cable 3-D−LF results with FWE at 1 Hz: real (E) and real (J) with V on the boundary (two figs.on the left: cross and longitudinal sections), and with V inside (two figs.on the right: cross and longitudinal sections).explicit setting of BC (left) and validation (right) are depicted in Fig. 2.

Fig. 6 .
Fig. 6.CPW−LF to HF characteristics.There are no instabilities at LF if V is inside the domain, and excitation with essential BC is used.

Fig. 7 .
Fig. 7. Coaxial cable−left: Y 11 phase versus frequency obtained with FWH and voltage excitation.The small conductivity in air (fictitious dielectric losses) stabilizes the solution at LF. Right: regime map for a point in the air with frequency range [1 Hz(rhombus), 600 MHz(square)].The star, triangle, and circle correspond to 10 kHz, 50 MHz, and 300 MHz, respectively.Red and black markers are for σ air = 10 −6 and 10 −12 S/m, respectively.

Fig. 8 .
Fig. 8. Cylinder−2-D axi (left), 3-D (right).Conditioning number of assembled systems.In 2-D, FWH does not need any stabilization, and the model is at least three times cheaper than any FWE model.In 3-D, FWE is the best option, at HF with boundary V , at LF after hybridization with a dummy field.

Fig. 9 .
Fig. 9. Coaxial cable−2-D axi (left), 3-D (right).Conditioning number of assembled systems (missing points correspond to ∞).Best results in 2-D are obtained with the H formulation.Best results in 3-D are given by the hybridization with the dummy scalar field with a suitable scaling factor p.