Electric circuit element boundary conditions in the finite element method for full-wave passive electromagnetic devices

A natural coupling of a circuit with an electromagnetic device is possible if special boundary conditions, called Electric Circuit Element (ECE), are used for the electromagnetic field formulation. This contribution shows how these ECE boundary conditions can be implemented into the 3D-finite element method for solving coupled full-wave electromagnetic (EM) field-circuit problems in the frequency domain. The frequency response allows the extraction of a reduced order model of the analyzed device, accounting for all the EM field effects. The implementation is based on a weak formulation that uses the electric field strength E strictly inside the domain and a scalar potential V defined solely at the boundary. Edge elements for E are used inside the three-dimensional domain and nodal elements for V are used on its two-dimensional boundary. The weak formulation is described and implemented in the free environment Open Numerical Engineering LABoratory (onelab). The validation is carried out on 3D examples.

Embedded field/circuit systems. Left: Coupling of electric circuits and EM device models are naturally ensured by means of terminals. Right: To ensure the coupling, "node voltages" (potentials) and electric currents of non-isolated circuits must have a correspondent in the EM device model The power transferred to it is if i m is expressed according to Kirchhoff current law for a cutset and the terminal m is connected to ground. This power expression shows that the state of a m-terminal circuit is characterized by 2(m -1) independent quantities: m -1 currents and m -1 voltages. The assumption v m = 0 is not a restriction for the purpose of this paper, which is stated at the end of Sect. 2. A natural coupling of this sub-circuit with an EM device is possible if some connecting surfaces are defined on the device boundary, for which currents and potentials are defined, in order to satisfy Kirchhoff relationships and provide the same transmitted power formula (1) as subcircuits do. The conditions that satisfy these requirements are the ones proposed in [20], used in [12,17] and referred to Electric Circuit Element (ECE) boundary conditions. The ECE boundary conditions, combined with current excited terminals, are the "realistic boundary conditions" used in [2] to solve eddy current problems with the finite element method (FEM) using a formulation in H and an ungauged Tϕ, ϕ one in [1]. Similar conditions, although with a different definition for the terminal voltages are proposed in [15] and used for A, V eddy current formulations [13]. In [8] the coupling of MQS (formulation in H and formulation in A, V ) FEM models with circuits is done by using electric global quantities that are naturally coupled with local quantities. The use of ECE in MQS problems for inductance extraction with an A, V formulation is discussed in [18].
Our aim is to use ECE boundary conditions to solve full-wave (FW) problems with FEM. We have successfully used ECE to model passive on-chip components such as resistors, inductors, capacitors, interconnects or RF-MEMS switches in FW [5], with the Finite Integration Technique as numerical method. To the best of our knowledge, the ECE conditions are not available in FEM codes which implement the formulation of microwave ports for FW. Setting such ports in FW requires complicate modeling e.g. as explained in [6]. Theoretical studies exists, e.g. in [12], based on an E, V formulation for the whole domain. In [4] we used E strictly inside the domain and V solely on the boundary. During the reviewing process of [4], Hiptmair and Ostrowski released a relevant report published in [14] proving the interest for this subject. In [19] the same authors propose a FW FEM formulation coupled with circuits, with a strong emphasis on the stability of the field solution inside the domain at low frequencies. An EQS formulation is herein used to gauge the proposed FW formulation in terms of potentials A and ϕ. Our approach in [4] states the Electric terminals are disjoint surfaces on the domain boundary. The non-grounded terminals can be either voltage excited (its potential is given) or current excited (its total current is given) problem in terms of the fields and does not need gauging. We focused only on coupling field devices with circuits at high frequencies.
This paper is an extended version of [4], where the validation was carried out on 2D examples. The proposed formulation is implemented in the Open Numerical Engineering LABoratory (onelab) environment, consisting of the free mesh generator Gmsh [9,10] and the finite element solver GetDP [7]. Here the validation is carried out on 3D examples. The implementation and the described test problems will be added to the collection available at https://gitlab.onelab.info/doc/models.
According to Faraday's law, (ECE3) implies (ECE1) for the terminals, the inclusion of the terminals in (ECE1) is only kept to emphasize the physical meaning. By definition, the currents and potentials of any terminal are: where, in order to ensure conservation, each terminal current is the total current (conductive and displacement) and the potential is properly defined as the voltage between this terminal and the reference one, along a path C k included in the domain boundary. Due to (ECE1) the voltage between two points placed on the boundary surface is independent of Figure 3 Each non-grounded terminal of the EM device with ECE boundary conditions can be either current excited or voltage excited. Its hybrid transfer matrix is obtained after computing voltages of the current excited terminals and currents of the voltage excited terminals in linear problems the path of the integration line connecting these points, provided that this path is included in the surface. Thus, the potential on the surface is well defined, although this is not the case for a general time-varying EM field. Under these conditions, (currents and voltages satisfy Kirchhoff equations and (1) holds for the EM device, where i k and v k are given by (2), and thus the ECE boundary conditions are compatible with the power transferred through its terminals by a multipolar circuit [17,20]. If we assume that the potentials at the terminals are known, then it can be proved that the problem of EM field analysis described by the Maxwell equations in a passive linear domain with ECE boundary conditions has a unique solution in the frequency domain. Consequently, the terminal currents are output signals and are obtained by solving the field problem [20]. As the domain is linear, so are the equations, hence the device with ECE conditions is a linear system, defining a multiple input multiple output (MIMO) type dynamic system with m -1 inputs and m -1 outputs (Fig. 3). In the frequency domain, the input-output relationship is: The problem to be solved is: "Find H(f ), where f is the working frequency in a given frequency range of interest, defined by its minimum and maximum values f min and f max f ∈ [f min , f max ], from the EM field solution. " If this hybrid matrix H(f ) is known, then the "field" element can be realized with common circuit elements and included in any circuit simulator. This final step can be carried out with vector fitting (VF), a very efficient model order reduction procedure which computes a rational approximation of a given order for the frequency response [11]. Moreover, VF can be included in a loop that increases the order of the reduced model while choosing new appropriate frequency points for evaluation in an adaptive frequency sampling strategy (AFS). Thus, for an imposed accuracy, the minimal reduced order model which describes high frequency field effects is obtained [3].

ECE in FEM (a) Strong formulation for E with classical boundary conditions
It is useful to recall the formulation in E with classical boundary conditions, since the newly proposed formulation inherits part of it. We will assume a frequency domain for-mulation, and therefore complex representation of the vector field quantities will be used and denoted by 1 E(r) = C(E(r, t)).
The well-known FW Maxwell equations in the frequency domain, for passive linear media and no internal field sources are: where permittivity ε, permeability μ and conductivity σ are positive, space dependent material parameters. The reluctivity ν = 1/μ might be used instead of μ. The solution of these equations is unique if in any point of ∂ , either exclusively E t or H t are known (given). The subscript t indicates the tangential component of the vector on the surface. It is useful to denote a disjoint partition of the boundary: In what follows we will name them classical boundary conditions. The uniqueness of the field solution can be proven on the basis of the complex form of the Poynting theorem that gives the expression of the transmitted power (assuming a linear field domain, with no moving parts): The proof assumes that there exist two such fields that satisfy the same boundary conditions. This means that the Poynting theorem in complex form is valid for the difference field, which satisfies the Maxwell equations (due to linearity) and zero boundary conditions. This implies that the real part is zero which results in zero difference electric field (conductivity of the domain is assumed non-zero everywhere) and the imaginary part is zero which amounts to zero difference magnetic field. The second order equation satisfied by E, obtained from Maxwell equations by elimination of H is a curl-curl type PDE: This equation has an unique solution when the following boundary conditions are given:

formulation for E with classical boundary conditions
In general, solving (5) implies a numerical approach, e.g. FEM, which is based on weak formulations. The needed functionals are obtained by projecting (5) onto a set of appropriate test functions E , then integrating by parts and applying Gauss-Ostrogradski formula: Replacing the expression of the magnetic field strength in the right hand side we get 1 Notation: An underlined symbol represents a complex quantity in order to avoid confusions with the root mean square (rms) value of the quantity it represents. For instance, the time representation of the potential is V : × (0, T) → R and its complex representation is V : → C, where the instantaneous value is V(r, t) = √ 2V(r) sin(ωt + θ), and its complex representation is V(r) = V(r) exp(jθ ), with V(r) the rms value, j the imaginary unit and ω = 2π f .

With classical boundary conditions, the integral in the right hand side is
E t are essential boundary conditions that is why the test functions are chosen so that E t is zero on S E . Thus, the weak equation for the trial functions E is: The boundary conditions H t = n × (H × n) = -n × (∇ × E × n)/(jωμ) are natural, they appear in the functional equation.
Finally, the weak formulation in E with classical boundary conditions is: where are the curl-conform Sobolev spaces.

(c) Weak formulation for E with ECE boundary conditions
If we use ECE boundary conditions, the unknowns are the electric field inside the domain and an electric scalar potential solely defined on ∂ . That is why the formulation is still named E, V , but is different from other formulations, such as the E, V in [12] where V is defined also inside the domain. An E, V interpretation of the ECE boundary conditions (ECE 1, 2, 3) is: • (ECE1b) E · dl = 0, ∀ ∈ ∂ , is a closed curve; • (ECE2b) n · E(r) = 0, ∀r ∈ ∂ -m k=1 S k ; • (ECE3b) E t (r) = 0∀r ∈ S k , k = 1, . . . , m. From (ECE1b) an electric scalar potential V can be defined on the boundary ∂ , such that E t = -∇ 2 V . Condition (ECE3b) requires that the electric terminals are equipotential. For uniqueness reasons, at least one terminal has to be fixed to a value. Without lack of generality we can assume it is grounded in what follows. For the other terminals the uniqueness implies that, exclusively, either their voltages or currents are known.
Using (5) we get the weak equation for E: where I c is the set of indices of current excited terminals. Similarly, we will denote by I v is the set of indices of voltage excited terminals.
In [4] a separate equation for the electric potential on the boundary was proposed. However, the coupling between the unknowns strictly inside the domain and those on the boundary can be done at the level of the basis functions, as follows.
Find E ∈ H E and V ∈ H V , so that and (d) Discrete formulation in E In [4] we used a simplicial mesh (tetrahedrons in 3D, triangles in 2D), numerical test functions N k that correspond to edge elements of order (0,1), and degrees of freedom that represent the complex representations of voltages U k along the edges. In the case of using classical boundary conditions, the numerical trial function is approximated as where Ne is the total number of edges in the domain, including its boundary.
In the case of using ECE boundary conditions, the function space where the trial function is searched for is curl free on the domain boundary, where nodal unknowns V k and test functions ϕ k are needed. The connection between the approximations inside and on the boundary can be done at the level at test functions. For instance, since for one element it follows that the numerical trial function when using ECE boundary conditions is approximated as where NeInt is the total number of edges that are strictly inside the domain and NnBnd is the total number of nodes on the boundary. Some of the nodes that are on the boundary belong to the same terminal, which must be equipotential. The corresponding terms in (20) have to be grouped together, and the final expression of numerical solution with respect to the trial functions is: where m is the total number of terminals, and NnTermK are the number of nodes that are covered by terminal k. The discrete formulation (15)- (17) was implemented in getdp [7]. The most relevant one for the formulations above are: the discrete function spaces (18) for classical boundary conditions and (21) for ECE; the weak equation (9) or (15), with the same bilinear functional (10) or (16), but with the linear functional (11) for classical boundary conditions (the tangential component of the magnetic field strength is a natural boundary condition) and (17) for ECE (currents of current excited terminals are natural boundary conditions); constraints -essential boundary conditions -tangential component of the electric field in the case of classical boundary conditions, from which appropriate voltages along edges on the boundary are computed, and potentials of voltage excited terminals in the case of ECE boundary conditions.

Numerical results
The results in [4] were obtained with an in-house code for two simple though relevant academic 2D test problems: a homogeneous rectangular domain corresponding to a single input single output system (m = 2) and a non-homogeneous rectangular domain corresponding to a multiple input multiple output system (m = 3). Two 3D test cases are considered herein. The new implementation has been done in onelab, using first order tetrahedral elements.
The first test is a cylindrical domain with radius a and length l, having linear and homogeneous material properties 2 (ε, μ, σ ). Its ends are two terminals, one grounded and Figure 4 Cylinder benchmark: Quantitative validation of the implementation for a 3D case with analytic solution, extracted resistance R (left) and inductance L (right) versus the frequency. Data used: radius a = 2.5μm, length l = 10μm, conductivity σ = 6.6 · 10 7 S/m, μ r = 1, ε r = 1. These results illustrate also the importance of the mesh adaptation to the field solution the other excited either in current or voltage. This configuration has the advantage that a formulation with classical boundary conditions 3 is equivalent to a formulation with ECE boundary conditions. 4 The classical boundary conditions formulation admits an analytic solution in terms of Bessel functions for the current excitation case. This is used to validate the numerical solution of FEM, in 3D-FW regime with ECE boundary conditions. Figure 4 shows a quantitative validation in terms of the associated frequency-dependent resistance (R) and inductance (L), the relative errors (vs. the analytic exact solution) are less than 2% for the whole frequency range when using a mesh with 2 elements per skin depth. The figure includes also the results obtained with the formulation proposed in [12] that uses V defined in the whole domain. The two numerical solutions overlap (up to machine precision), but the computational effort for a formulation with V inside is greater, since the number of DoFs increases with the number of inner nodes. For instance at 60 GHz, if 2 elements per skin depth are imposed, then a mesh with 128,219 elements is generated, leading to 12,659 and 23,441 unknowns for V when it is defined only on the boundary, or in the whole domain, respectively.
The second benchmark is the LC 3D field device proposed in [19]. Not only this test was solved with FEM, 5 in FW regime with ECE boundary conditions, but the solving was embedded into a model order reduction procedure based on Vector Fitting (VF) [11] and Adaptive Frequency Sampling (AFS) [3]. Figure 5 shows a part of the surface mesh and the color map of the real part of the potential on the boundary, as well as a color map on a cut of the rms value of the electric field strength for f = 1 kHz. The mesh used has 161,485 first order tetrahedral elements leading to 187,560 DoFs. Note that in [19] a mesh with 1.12 million elements was used to solve this problem. Figure 6 shows the results of the model extraction using FEM for field evaluation combined with model reduction with Vector Fitting (VF) and Adaptive Frequency Sampling (AFS). The reference result corresponds to a series connection of three lumped elements:  [19], especially if instead or the DC value of resistance R ref = 24 m is used, which considers the eddy-current effects effect at the resonance frequency. 6 We have set the AFS local relative error to 0.01. A reduced model of order 4 was obtained, needing only 7 frequencies evaluations in FEM to achieve this error. The evaluation of the transfer function (TF) in the frequency range shows that the resonance frequency is at 8889 kHz, which is 2.5% different from the reference resonance frequency. Moreover, a new strictly proper rational approximation of order 2 was obtained, starting from a set of 100 points evaluated for the reduced model of order 4. From the two poles and residuals, the values of a RLC series circuit were extracted: R extr = 26.6 m ; L extr = 9.18 mH; C extr = 33.8 mF, which correspond to relative errors versus the reference circuit of 10.8 %, 9.1% and 1.1%, respectively. As in the previous cylinder benchmark, better numerical results are expected for a refined mesh. The circuit simulation is also shown in Fig. 6. Figure 7 shows the results obtained when connecting the extracted reduced models to the outside circuit consisting of a resistance R ext = 450 m from [19] and an ideal voltage source having an imposed voltage of 1 V (rms value). The results are extremely good, Figure 7 LC benchmark, models of the 3D field device connected to the external circuit: admittance (modulus and phase) obtained by the reference lumped model from [19], AFS points (displayed here only to recall their position, no field solution is needed here), reduced model of order 4 (transfer function of VF), extracted RLC model especially considering that the mesh we used is about seven times coarser than the mesh in the reference paper. No instabilities were noticed at terminal values for frequencies lower than 1 kHz, down to 50 Hz, even if the conductivity in the dielectric domain is zero.
The external resistance which is relatively high makes that the accuracy of the extracted resistance of the LC device have no relevance.

Conclusions
The main advantage of ECE BC for the Maxwell equations is that the ports are well defined, without ambiguity, and compatible with the circuit terminals, even for RF devices. There is no restriction on the field regime (DC to full wave, even including nonlinear media). For MIMO systems, the hybrid excitation is obtained in a natural way. It is important to be aware that ECE BC for parameter extraction can be applied only to a simply connected subdomain, obtained after partitioning the domain corresponding to a whole system in parts that do not overlap or do not have holes.
This paper proposed a FEM formulation for ECE, which E strictly inside the domain and V on the boundary. The degrees of freedom are the electric voltages on the inner edges and the potentials of the floating nodes on the boundary (nodes outside terminals and current excited terminals). The paper presents its theoretical background (weak formulation) and its implementation details in the onelab environment. This new approach provides the same results but is more data-efficient than using internal node-potentials as described in [12]. The simulated frequency response allows extracting of a reduced order model of the analyzed device, accounting for all the EM field effects [3]. The formulation was validated for two 3D test problems, the frequency responses being saved in Touchstone file format. The numerical experiments showed that besides fundamental physical and mathematical aspects related to the weak formulation and functional framework, the choice of the mesh and its adaptation to the solution are key to ensure a highly accurate of the numerical model. Besides the appropriate choice of the field regime and boundary conditions which are apriori model order reduction (MOR) techniques, and compact "equivalent" circuit extraction which is an aposteriori MOR technique, the choice of an optimal mesh is an "on the fly" MOR technique. Effective methodologies for solving real-life field problems should use MOR at every modeling step [16].