Marching On-In-Time Unstructured PEEC Method for Electrically Large Structures with Conductive, Dielectric, and Magnetic Media

: The Marching On-In-Time (MOT) unstructured Partial Element Equivalent Circuit (PEEC) method for time domain electromagnetic problems is presented. The method allows the transient analysis of electrically large electromagnetic devices consisting of conductive, dielectric, and magnetic media coupled with external lumped circuits. By re-formulating PEEC following the Coulombian interpretation of magnetization phenomena and by using electric and magnetic vector potentials, the proposed approach allows for a completely equivalent treatment of electric and magnetic media and inhomogeneous and anisotropic materials are accounted for as well. With respect to the recently proposed Marching On-In-Time PEEC approach, based on the standard (structured) discretization of PEEC, the method presented in this paper uses a different space and time MOT discretization, which allows for a reduction in the number of the unknowns. Analytical and industrial test cases consisting in electrically large devices are considered (e.g., the model of a Neutral Beam Injector adopted in thermonuclear fusion applications). Results obtained from the simulations show that the proposed method is accurate and yields good performances. Moreover, when rich harmonic content transient phenomena are considered, the unstructured MOT–PEEC method allows for a signiﬁcant reduction of the memory and computation time when compared to techniques based on Inverse Discrete Fourier Transform applied to the frequency domain unstructured PEEC approach.


Introduction
Integral equation (IE) methods are particularly suited for the study of electromagnetic (EM) devices in unbounded domains since they do not require the introduction of artificial absorbing boundary conditions for truncating the computational domain and they avoid the discretization of regions with the characteristic of vacuum [1][2][3]. In particular, the Partial Element Equivalent Circuit (PEEC) method, introduced by Ruehli in 1972 [4], has been shown to be well suited for the analysis of several classes of EM devices, such as interconnecting [5,6] antennas [7][8][9], and other complex electric and electronic devices [10][11][12][13]. IE methods such as the PEEC approach are widely used for time-harmonic EM problems and they are also particularly suited for transient analysis, usually performed by means of traditional time-stepping algorithms, e.g., forward and backward Euler, Crank-Nicolson, and Runge-Kutta methods.
However, in the context of transient analyses with electrically large devices, where time delay effects on the propagation of the EM fields need to be considered, traditional time-stepping algorithms cannot be easily used. Among all the different possibilities to develop Time Domain Integral Equation (TDIE) methods, the Marching On-In-Time (MOT) scheme [14] allows for naturally considering the effects of the time delay on the EM fields propagation and such scheme has been been widely adopted over the years [15].
The possibility of combining the MOT procedure with the structured version of PEEC has been shortly discussed in [16,17], whereas, in the literature, time domain PEEC methods under the assumption of electrically small structures have been mostly proposed (i.e., without time delay effects). Recently, a MOT-PEEC method based on the structured PEEC approach has been presented in [18] for the study of conductive and dielectric media only (i.e., non-magnetic media). Another PEEC based ad hoc approach for the study of a wire structure excited by a distant lightning channel has been proposed in [19].
The aim of this work is to provide a general and efficient numerical tool by combining the MOT scheme with the unstructured version of the PEEC approach [20]. Indeed, differently from existing MOT-PEEC approaches, the formulation proposed herein is based on the unstructured version of PEEC [10], and the MOT procedure is differently applied to the PEEC formulation with respect to [16]. As a further innovative feature, unlike [18], magnetic media are also treated by the proposed MOT-PEEC method by following the recently proposed harmonic PEEC formulation based on the Coulombian interpretation of magnetization phenomena [21]. This also allows for a reduction in the number of unknowns. Furthermore, the proposed formulation provides the same well-known circuit interpretation of the harmonic PEEC method, allowing for a very easy coupling of discretized devices and lumped circuits' components. Moreover, inhomogeneous and anisotropic media can be considered as well. Therefore, the typical target problems for the MOT-PEEC approach are electrically large devices consisting of a wide class of materials. For instance, the method is well suited for the study of scattering problems. Moreover, thanks to the features of integral equation methods which avoid the air meshing, the time domain response of large and complex devices to rich harmonic content excitations in unbounded domains can been efficiently evaluated. Indeed, as shown in the numerical results section, with respect to approaches based on standard Inverse Discrete Fourier Transform applied to frequency domain methods, the proposed MOT-PEEC technique allows for a significant reduction in the computation time and memory requirements.
The rest of the paper is organized as follows. First, in Section 2, the delayed time domain PEEC formulation based on the frequency domain PEEC formulation proposed in [21] is presented. Then, in Section 3, the spatial and temporal MOT discretization is applied to the unstructured-PEEC formulation. In the last part of Section 3, some numerical aspects concerning the proposed unstructured MOT-PEEC scheme are discussed and possible extensions and modifications of the numerical approach based on the vast literature concerning MOT and TDIE are discussed. In Section 4, to show the usefulness of the MOT-PEEC method for the study of different problems which include electrically large structures, several numerical test cases simulated with the proposed unstructured MOT-PEEC approach are presented and the numerical results show good performances in terms of both computation time and memory. First, the MOT-PEEC is validated by studying the analytical case of a dielectric shell excited by an EM plane wave. Then, the case of a Neutral Beam Injector used for Thermonuclear Fusion applications is considered [22,23]. Finally, the numerical results obtained from a MOT-PEEC method for equivalent surface models are presented and the same stability analysis proposed in [24] is applied to one of these cases.

MOT-PEEC Formulation
In this section, the formulation of the Unstructured MOT-PEEC method is presented. In this work, magnetic media are considered by following the approach proposed in [21], based on magnetic currents and charges (i.e., the Coulombian interpretation of the magnetization phenomena [21,25,26]).
Let Ω c , Ω d , and Ω m be the conductive (σ > 0), the dielectric (ε r = 1), and the magnetic (µ r = 1) domains, respectively, where σ is the electric conductivity, ε r is the relative electric permittivity, and µ r is the relative magnetic permeability. The electric domain is also defined as Ω e = Ω c ∪ Ω d and the following boundaries are defined: Γ c = ∂Ω c , Γ d = ∂Ω d , Γ m = ∂Ω m , and Γ e = ∂Ω e . The computational domain is Ω = Ω e ∪ Ω m . As in [21], where a time harmonic PEEC method is proposed, the electric, E, and magnetic, H, fields can be written in terms of electric and magnetic potentials, i.e., where A e and A m are electric and magnetic vector potentials, ϕ e and ϕ m are electric and magnetic scalar potentials, and E 0 and H 0 are external electric and magnetic fields produced by the source domain Ω 0 . Electromagnetic quantities in (1) and (2) are evaluated in the field point r and at the present time instant t. Moreover, as in [21], the conduction, J c , the polarization, J d , and the magnetic, J m , currents are introduced together with the following time domain constitutive relations: where ρ c = σ −1 is the electric resistivity, is the equivalent electric resistivity of dielectric media and is the equivalent magnetic resistivity of magnetic media. ε 0 and µ 0 are the vacuum permittivity and permeability, respectively.
In (1) and (2), A e and its time derivative are given by where J e = J c + J d is the electric current density, r is the integration point, and t = t − |r−r | c 0 is the retarded time; c 0 is the speed of light in vacuum.
Similarly, A m and ∂ ∂t A m are given by: The scalar electric potential is given by: where e and ς e are the volume and surface electric charges [21], respectively. e and ς e are related to J e by means of charge conservation relations: where n is the outward unit normal vector of Γ e . Equivalently, ϕ m is given by where m and ς m are the volume and surface magnetic charges, respectively. Like (13) and (14), m and ς m are related to J m by means of charge conservation relations: Finally, as usual with PEEC discretization scheme, (1) and (2) are combined with (3)-(17) resulting in final integral equations where the current densities and the scalar potentials are considered as the problem unknowns.

Spatial and Time Discretization
The computational domain Ω is now discretized into a tetrahedral (or 8-node hexahedral elements with planar quadrilateral faces) unstructured grid G, as proposed in [21]. Moreover, a temporal discretization T of the time interval of interest is also introduced consisting of N T time steps, equally spaced by ∆ T .
Then, J c is expanded in terms of degrees of freedoms (DoFs) by using face functions w k for the spatial discretization and temporal basis functions T s for the time discretization [3]: where j (s) c k , with k = 1, · · · , N f c , is the flux of J c through the kth face of the mesh (i.e., the kth branch of the equivalent circuit) at the sth time step, N f c and N t are the number of spatial (i.e., faces of the mesh) and temporal basis functions, respectively.
At each sth time instant (i.e., t = s∆ T ), with s = 0, 1, · · · , N T , the DoFs are stored in the array j For simplicity reasons, from now on, T is considered as the well-known hat shape function with support [−∆ T , ∆ T ], where ∆ T is the chosen time step. This leads to a causal MOT scheme. However, as discussed in the following, different choices of the basis function can be made, leading to non-causal MOT systems [3].
Following [14], a Galerkin approach for the spatial discretization and a collocation method for the temporal discretization are applied to the equations obtained by combining (1) and (2) combined with (3)- (17). At each sth time step, this operation leads to the following matrix equations: for the conductive media, for the dielectric media, and for the magnetic media.
In (19)-(21), e The coefficients of the resistance (R) and inductance (L) matrices in (19)-(21) are given by where t u = u∆ T − |r − r |/c 0 and T(0) = 1 when temporal hat shape functions are used. T is the temporal (hat) basis function centered in the origin (i.e., T = T 0 ). Matrices R α and L αβ u have dimension N f α × N f β , where, for the sake of conciseness, α = c, d, m and β = c, d, m depending on the domain.
The coefficients of the mutual interaction matrices between electric and magnetic domains (K) are instead given by: Matrices L and K in (19)- (21), with u = 0, · · · , H T , represent the electromagnetic interactions between the unknowns at the time instant s and unknowns at the previous time instant s − u. When u = 0, these matrices represent instantaneous interactions between unknowns at the same time instant. Instead, when u > 0, these matrices represent retarded interactions between unknowns that are separated in time by a quantity equal to u∆ T . H T = 1 + D max /(∆ T c 0 ) indicates how many previous time steps actually interact with the present solution, where D max is the maximum distance between two mesh elements of Ω and • is the ceiling operator.
These instantaneous and retarded interactions are exemplified in Figure 1: unknowns which are close to each other interact instantaneously (orange cross vs. red unknowns). Instead, unknowns which are far from each other interact with a retardation which increases with the spatial distance (orange cross vs. green, blue, and black unknowns). Following a leap-frog scheme [17] and according to the classical Tonti diagram, ϕ e and ϕ m are expanded onto the dual spatial and temporal gridsG andT obtained by the barycentric subdivision of G and T , respectively. Figure 2 shows the primal and dual temporal grids: T andT , respectively. As in [21], φ where, as in [21], the potential matrices (P) are subdivided into volume, surface, and volume-surface matrices, i.e., P ss αβ u,kh P vs αβ u,kh where, for the sake of conciseness, α = c, d and β = c, d, and | · | indicates the measure (i.e., volume or area) of either volume or surface elements.
where subscript l in (46) Finally, combining (34)-(36) with (46)-(51), one obtains: Equations (19)- (21) and (52)-(54) (scaled by ∆ T to improve the conditioning) are finally combined, leading to the final system of linear equations to be solved at each sth and lth time step where (for each sth and lth time step) the unknowns are j The left-hand side of (55) consists of resistance, incidence, and instantaneous matrices only (i.e., those ones with u = 0) and is the same for each time step. The right-hand side is instead updated at each time step by multiplying the marching matrices (i.e., ones with u = 1, · · · , H T ) with the previous H T solutions, as shown in (56)-(61). It is worth noting that instantaneous and marching matrices are in general sparse with a sparsity ratio which depends on the choice of ∆ T (i.e., the smaller ∆ T , the lower the sparsity ratio and the greater H T ). As a distinctive feature of PEEC, connections with lumped circuit elements are easily considered by adding the circuit currents and node potentials of the lumped components to the unknowns and solving Kirchhoff's Voltage and Currents Laws together with (55).
With the aim of increasing the stability and the accuracy of the MOT scheme applied to integral equation methods, several authors have proposed different temporal shape functions (e.g., higher order with a more extended [27] and/or not-causal support [3]) and different choices of the unknowns (e.g., electric and magnetic polarizations P e and M or densities fields D and B). Indeed, these different choices (which also allow for reducing the number of unknowns by neglecting the scalar potentials from the independent unknowns [3]) seem to generally improve the accuracy and the stability of the MOT method. Indeed, the use of a smoothly varying temporal shape functions alleviates the instability issues caused by the numerical integration of integral matrix coefficients subject to the encroachment of the propagating electromagnetic wave [28,29] (see Figure 3). However, the use of temporal shape functions with large support deteriorates the sparsity ratio of the matrices. Moreover, when temporal basis functions with non-causal support are chosen, several future solutions must be predicted as proposed in [30].
In the context of the proposed MOT-PEEC method, in accordance with the analysis presented in [18], hat shape functions have been chosen since they provide a good trade-off between stability, accuracy, and numerical performances. However, the discretization process proposed herein can be easily modified by considering different temporal shape functions and/or predictor-corrector schemes [30]. Indeed, the proposed MOT-PEEC method (which is intended as a basic and general formulation) can be however modified by following more sophisticated approaches [31][32][33][34] and parallelization techniques, as proposed in [15].
Lastly, the choice of considering the scalar potentials as independent unknowns makes it possible to separate the inductive (L) and capacitive (P) effects, avoiding instabilities which can be attributed to round-off errors equivalent to those due to the breakdown in frequency issue in harmonic problems.

Numerical Results
The MOT-PEEC code has been developed with MATLAB R and parallel MEX-FORTRAN functions based on OpenMP libraries. Thin wire, triangular, and hexahedral mesh elements are considered. Simulations were run on a Linux machine equipped with dual 6-core/12-thread processors (Xeon E5-2643 v4 @3.40 GHz) and 512 GB RAM.
The results obtained from the MOT-PEEC code are compared with those obtained from a frequency-domain PEEC (FD-PEEC) code [21]. For case 1 only, the results are also compared with the MIE solution. Discrete Fourier Transform (DFT) and Inverse-DFT (IDFT) are adopted to extract frequency domain results from time simulations and vice versa. Figure 4 shows the scattered magnetic field magnitude (in time domain) evaluated in the target point [5 m, 5 m, 5 m] for case 1). For cases 1, 2, and 3, Figure 5 shows the frequency spectra of the three components of the scattered magnetic field. The same hexahedral mesh with 3072 elements is adopted for both MOT-PEEC and FD-PEEC.
Results show good agreement and some small discrepancy can be observed only for the smallest components of H.

Neutral Beam Injector
Neutral Beam Injectors (NBIs) are electrically large devices adopted in Thermonuclear Fusion for plasma heating. In these devices, the protection against the grid breakdown and other possible failures is a very sensitive issue [22,23]. NBIs are complex devices fed by very long transmission lines (e.g., ∼ 70 m for ITER), which, with good approximation, can be modeled as coaxial cables. With the aim of protecting NBIs against possible breakdowns (i.e., damping dangerous voltage oscillations), magnetic core snubbers (combined also with other technologies) are usually adopted. Here, the simplified case of NBI shown in Figure 6 is considered.  Figure 7 shows the results obtained by truncating the sweep after 10, 20, 30, 70, and 150 frequency simulations. Since retardation effects must be considered, the dense integral matrices of the FD-PEEC must be updated at each frequency simulation. Thus, to avoid prohibitive computation time of FD-PEEC, the exponential term of the (harmonic) dynamic Green's function is considered to be constant within every single cell (with a good approximation) in order to avoid the re-evaluation of the integrals at each frequency step. Thus, the inductance and potential PEEC matrices are updated at each frequency value by multiplying each matrix entry with the related (retardation) coefficient.
The Peak Memory Usage (PMU) of the MOT-PEEC method is 9.7 GB, the time for the generation of the matrices is T g = 324 s, and the time for a single time solution is T s = 1.4 s. For the sake of comparison, the computational details of the IDFT FD-PEEC are 15 GB of PMU (due to the need of storing the non-retarded matrices, the retardation coefficients, and the system to be solved), T g = 254 s, and T s = 49 s (i.e., time for the updating of the matrices and the assembling and solution of the system). Thus, the total simulation time of MOT-PEEC is T t = 954 s, whereas the one of FD-PEEC is T t = 7604 s. Indeed, due to the rich harmonic content of the problem, the MOT-PEEC code allows for significantly reducing the computational effort with respect to the FD-PEEC.

Equivalent Surface Models
The MOT-PEEC formulation can be easily modified for the study of equivalent surface models. In this paragraph, some of the numerical results obtained from a MOT-PEEC formulation for equivalent surface models (i.e., thin structures) are briefly presented. Figure 8 shows a snapshot of the electric potential wave propagating on an airplane excited by a GMP (62) ( f 0 = f bw = 100 MHz, ς = 6 2π f bw , t 0 = 4ς). The model consists of 33,290 triangular elements, resulting in 83,240 unknowns (i.e., the edges of the mesh). The resistivity of the plane is ρ c = 7 · 10 −3 Ωm, with a thickness of 1 mm. A time step of ∆ T = 0.5 ns has been used and the simulation consisted of N T = 800 time steps. The computation time required by the simulation was T g = 927 s, for the generation of the MOT matrices, and T s = 6.7 s for each time step solution. The PMU reached by the simulation was 78 GB.  Finally, the MOT-PEEC formulation for the case of equivalent surface conductive models is validated on the case of a 1 × 1 m 2 square plate with normal [−0.34u x , 0.81u y , 0.34u z ] excited by the GMP (62) with f 0 = 10 8 Hz, ς = 6 2π f bw , t 0 = 4ς, and f bw = 10 8 Hz. The thickness of the plate is 1 mm and the resistivity is ρ = 4 · 10 −5 Ωm. The chosen time step is ∆ T = 10 −9 s. Figure 10 shows the temporal evolution of the current density (x-component) in one of the triangular mesh elements obtained from the MOT-PEEC code for equivalent surface models. The results are compared with the ones obtained from the IDFT applied to a FD-PEEC code [20] and the well-known commercial software FEKO R based on IE methods.

Stability Analysis
Finally, when the dimension of the problem is relatively small, the stability analysis proposed in [24] and based on the eigenvalues distribution of the MOT transmission matrix can be applied (i.e., the method is unconditionally stable if the magnitude of all eigenvalues is smaller than 1). Figure 11 shows the eigenvalue distribution of the MOT-PEEC transmission matrix for the case of the square plate considered above.

Conclusions
The unstructured Marching On-In-Time Partial Element Equivalent Circuit method for the study of conductive, dielectric, and magnetic media has been presented with several application examples including inhomogeneous and anisotropic media. When fast transients need to be analyzed, the proposed method shows itself to be competitive with frequency domain methods combined with Inverse Discrete Fourier Transform. The proposed MOT-PEEC method is also suitable for the application of more sophisticated temporal discretization approaches [3] and parallelization techniques, as proposed in [15].

Conflicts of Interest:
The authors declare no conflict of interest.