Combined FVTD / PSTD Schemes with Enhanced Spectral Accuracy for the Design of Large-Scale EMC Applications

A generalized conformal time-domain method with adjustable spectral accuracy is introduced in this paper for the consistent analysis of large-scale electromagnetic compatibility problems. The novel 3-D hybrid schemes blend a stenciloptimized finite-volume time-domain and a multimodal Fourier-Chebyshev pseudo-spectral time-domain algorithm that split the overall space into smaller and flexible areas. A key asset is that both techniques are updated independently and interconnected by robust boundary conditions. Also, combining a family of spatial derivative approximators with controllable precision in general curvilinear coordinates, the proposed method launches a conformal field flux formulation to derive electromagnetic quantities in regions with fine details. For advanced grid reliability at dissimilar media interfaces, dispersion-reduced adaptive operators, which assign the proper weights to each spatial increment, are developed. So, the resulting discretization yields highly rigorous and computationally affordable simulations, devoid of lattice errors. Numerical results, addressing detailed comparisons of various realistic applications with reference or measurement data verify our methodology and reveal its significant applicability.


Introduction
The development of powerful time-domain solvers for contemporary electromagnetic compatibility (EMC) applications is solidly related to the correct representation of their size by the suitable spatial sampling rates [1][2][3][4][5][6][7][8][9][10].Considering that their overall size should be frequently reconfigured to comply with modern standards and several prototypes must be devised prior to the selection of the final design parameters, the use of rigorous discrete models seems a very attractive means to restrict high fabrication costs.Such a task, though, is often rather difficult, as most devices comprise electricallylarge homogeneous areas and fine details that require excessive overheads.To alleviate this hindrance, the finite-volume time-domain (FVTD) [11][12][13][14][15] and pseudospectral time-domain (PSTD) [16][17][18] techniques can be proven quite beneficial.However, the broadband operation of the preceding structures along with the usual frequency-dependence of their media properties, entail very large resolutions which are still expensive as well as potentially error prone.This is chiefly regular when variable-cell complexes are utilized near different mate-rial interfaces, as occurs in many EMC problems.Consequently, grid reflection arises, even if these interfaces are located in free-space regions.The strength of such undesired phenomena depends on the magnitude and abruptness of the grid wave mode velocity distribution and may be crucial in diverse large-scale applications.Pursuing the constant improvement of these important drawbacks, an assortment of proficient algorithms has been hitherto presented [19][20][21][22][23][24][25][26].
In this paper, a hybrid FVTD/PSTD method in general curvilinear coordinates with enhanced spectral precision for the consistent design and modeling of large-scale 3-D EMC structures is introduced.Established on a robust decomposition concept, the novel multimodal formulation launches a family of finite-volume operators for abrupt field variation and a Fourier-Chebyshev interpolation for periodic details.To this aim, spatial derivatives are evaluated via advanced operators of adjustable accuracy order to establish a conformal flux framework for all electromagnetic quantities under study.Both time-domain processes are separately updated, while robust continuity conditions guarantee their smooth boundary interconnection.Moreover, a stencil management technique sets the optimal lattice density and diminishes dispersion errors, while to retain consistency at media interfaces, adaptive nodal patterns are developed.In particular, by inserting auxiliary nodes at the interface, spatial increments can be computed through the convex combination of several possible candidates.Therefore, and unlike the unduly needs of existing approaches, discretization artifacts receive a notable reduction and total accuracy overwhelms typical thresholds.The proposed schemes, except their theoretical verification, are successfully applied to complex arrangements like large-scale planar microwave devices, antennas, waveguide structures, and anechoic/reverberation test facilities.

The hybrid 3-D time-domain method
A critical issue in large-scale EMC problems are the fine, regarding the minimum wavelength, geometric details with irregularly curved cross-sections [27].Despite their limited appearance compared to homogeneous areas, they create nonseparable wavefronts which demand high mesh densities.

Explicit multimodal formulation
To alleviate this serious defect in subwavelength regions, our hybrid algorithm launches an efficient multimodal concept with controllable resolution.Let us consider the discontinuity of Fig. 1, spanning over an angle θ T , with inner and outer mean radii r A (θ, φ) and r B (θ, φ).The basic point is the efficient separation of all propagating modes at κ preselected transverse planes, which fulfills the pertinent physical conditions.Thus, each component G can be expressed via infinite series as ( , , ) ( , , ) ( , , ) G r r g r ( , ) ( , , ) cos ( , ) r M (θ, φ) = r B (θ, φ) -r A (θ, φ), ξ κ the respective eigenfunctions and g κ amplitude coefficients.An important asset of (1)-( 3) is their conformal modeling capabilities, which enable the accurate description of the above details through the selection of r Having attained the required decoupling of the complicated modes, the parametric method projects Maxwell's laws on the prefixed κ planes.Hence, a system of differential equations is derived, whose solutions act as transverse intermediate excitation surfaces in the discontinuity.Such an idea is proven promising and the most significant; it retains lattice duality even near the geometric oddity, in contrast to typical algorithms, which cannot provide adequate treatment.Using the matrix notation, Ampere's and Faraday's laws become where E and H are the electric and magnetic field intensities, defined at the general coordinates (u,v,w).The elements of Λ i (i = A,…,E) impedance matrices provide the main details of the structure, under investigation, and are obtained by means of all field continuity boundary constraints, as

Enhanced spectrally-accurate FVTD/PSTD schemes
The principal part of the proposed methodology lies on the hybridization of generalized FVTD and PSTD schemes.Starting from the former, we establish a curvilinear dual-cell grid, where a new class of derivative approximators of adjustable spectral accuracy is developed.Concerning temporal integration, the process uses a stable integration as the one of the finite-difference time-domain (FDTD) approach [28].To increase its applicability, the 3-D algorithm involves cells with the four corners of their faces not necessarily lying on the same plane.Thus, the suitable unit vector direction cosines are employed to transform Maxwell's equations in the form of which may be regarded analogous to the main conservationlaw formulation.In (11), J s is the source vector that gathers all possible external excitations imposed to our problem.Moreover, the parametric operator [.] is denoted as where G = [E, H] T with the second term defining the novel Kth-order approximation forms that correct existing stencils around curved interfaces or small geometric details.Real parameters a l (l = 1,2,3) assure the stability of the schemes and subdue late-time exponential oscillations during the handling of discontinuities [30], by satisfying Kth-order constraint Actually, the above constraints constitute an alternative expression of the usual continuity conditions at the different media interfaces of the problem.Therefore, they can annihilate vector parasites that corrupt traveling waves in the domain.Note that for the sake of clarity, their values will be provided for every numerical simulation in the paper.Also, spatial derivatives in (12) are acquired via enhanced operators, whose spectral resolution is of order K. So, we define with η(u,v,w), Δη(Δu,Δv,Δw) is the spatial increment and M η,P [.] allowing the fine regulation of spectral precision by diverse stencils and smoothness summation limit P. Particularly, considering the system metrics, , , , , , , 1 1 Structural parameters T s and R s,p increase algorithmic consistency in the case of complex discontinuities, frequently encountered in realistic EMC devices.
Bearing in mind the above analysis of spatial derivative calculation, Maxwell's equations can, now, be expressed as μΙ}, and I the identity matrix.If our computational domain Ω is divided into a prefixed number of cells Ω i of boundary ∂Ω i , then ( 16) can be integrated over each Ω i and V i the volume of Ω i , and n 0 the outward unit vector perpendicular to surface ∂Ω i .Observe that Ξ[n 0 ]G represents the flux through ∂Ω i .To exploit all fluxes through the cell faces, the concept of the cell average -i.e. the average of G placed at the barycenter of every cell [30] -is described by 1 ( ) where r is the position vector.For the required flux splitting, based on our FVTD formulation, matrix which can be decomposed to a sum of matrices with positive and negative eigenvalues, is also used.So, recalling that the eigenvalues of F[n 0 ] are given by diag{0,0,υ,υ,-υ,-υ} for υ = (με) -1/2 and = [n 0 ], one arrives at The surface integral on the left-hand side of ( 17) can be evaluated if flux F[n 0 ]G is already known on all faces of ∂Ω i .So, two subvectors of G are considered; one on the internal and the other on the external side of each face.Both terms are acquired through an interpolation process centered on the respective face side.Replacing ( 17)-( 19) into (16), we derive for Ζ a bloc-diagonal matrix with all media features and Y containing the source contributions of every grid cell.It is mentioned that Z is inverted only once at the beginning of our computations [5].Regarding the explicit temporal integration of the proposed FVTD procedure, time derivatives in (20) are approximated by the simple leapfrog-like predictorcorrector scheme.In its general form, this scheme is given by 1 whose important property is the incorporation of ( 12)-( 19) in the temporal integration process.Herein, to achieve higherorder accuracy, a 3-stage scheme is devised.Denoting as K[.] the bracketed term on the right-hand side of (21), we get (1) ( , ) ( , ) ( , ), ( , ) 4 ( , ) 3 ( , ) ( , ) ( , ) ) that can be expressed in the more compact form of 1/ 2 0.5 where A is the typical update matrix for vector G.This stable 3-stage scheme enhances the precision of the FVTD method, since it can track abrupt field variations in a smoother and more physical way than the usual integration techniques.Focusing on the second part of the hybridization, a curvilinear Fourier-Chebyshev PSTD interpolation is developed for sections of homogeneous media and periodical features.Since they are the majority in Ω, their analysis via a pseudo-spectral approach is anticipated to accelerate the total simulation.This is attained by evaluating spatial derivatives in terms of with FFT and FFT -1 the Fourier and inverse Fourier transform, L η the node number towards η-axis, and k η the transformed wavevector component.Additionally, Γ describes all η coordinates along a straight-line cut through the transverse plane of the other two coordinates.In fact, (24) implies differentiation by the trigonometric Chebyshev-type polynomial with j 2 = -1 and the barred G(m) being the Fourier series of G.Note that field values are located at specific non uniform mesh points ψ i = πi/L η , i = 0,1,…,2L η -1.This permits space partition into collocated elements that conform to the boundaries of media distributions.Next, and in terms of a curvilinear coordinate transformation, each element is mapped into a cube, while all field quantities adjacent to internal interfaces are corrected to resolve the physical jump conditions.Regarding temporal update, all temporal derivatives appearing in the PSTD forms are computed through a similar leapfrog notion, as in (23) for the FVTD forms, given by thus facilitating the fully stable, non-oscillatory, and uniform time-advancement for the entire FVTD/PSTD method.A closer inspection of ( 11)-( 20) and ( 24), shows their ability to treat frequency-dependent setups owing to the direct tensor-oriented classification of constitutive parameters.So, dispersion errors are greatly minimized and broadband studies are easy to conduct.Finally, it should be emphasized that, except some limited storage needs, the proposed method does not increase the total CPU and memory burden.

Theoretical analysis
In the case of arbitrarily-curved interfaces formed by materials with considerably different constitutive parameters and losses, special attention must be drawn in order to circumvent severe instabilities and lattice discretization errors.Returning to ( 6), (7) or ( 16) and re-expressing cell average in (18) as with respect to an (iu, jv, kw) mesh, where V i = uvw is the volume of the p i,j,k cell and m c are functions of ,  constitutive parameters.The integral of ( 27) is calculated by a Gaussian quadrature along the edges and over the faces of p i,j,k .Hence, for a general cell extending in the region of the cell average at the (i, j+1/2, k) face is obtained via for s, q referring to the N extra integration points u s , w q and structural coefficients V s , V q .So, instead of ( 12) or ( 26), all spatial derivatives in ( 16) are effectively approximated by and the piecewise polynomials P η (r,t).Subsequently, we introduce a weighted convex combination S r of all stencils related to the p i,j,k cell to reconstruct P η (r,t).As a matter of fact, the novel (2M -1)th-order schemes assign an optimal weight to each stencil in the vicinity of the interface, thus offering increased accuracy and reduced anisotropy and dissipation errors.Taking into account that v is the direction of interest, the proposed stencil classification is denoted as for r = 0,1,…,M -1.Thus, a set of M interpolating polynomials, Y M,v (r,t), is recursively associated with these stencils, i.e.
, , 1 /2 , 0 , , , , at time-step nΔt, while a similar process holds for directions u and v. Coefficients b M , τ ξ,v are computed by 1 ( , , ) ( , / 2, ) 2 In (31), spatial increments are formed by the initial G quantities in contrast to (32) where cell averages are employed.For example, for M = 3, the first member of Y ξ,v (r,t) is Therefore, P v (r,t) are extracted by the summation of for a κ e below 10 -5 to evade inconsistencies.Parameters ρ l typify every material in the structure, while where Figure 2: Integration process for an electric field component at the interface of two dissimilar lossy materials by means of the complex-interface modeling method (grey area).
In this way, the hybrid method achieves fast and very accurate outcomes with prominent savings prior to actual implementation, since it can decrease large-scale EMC problems to much smaller interface sections unlike traditional schemes.

Hybrid realization and boundary treatment
To complete the interface modeling approach, a word on its application to the adjacent regions with the hybrid method should be devoted.Focusing on every complex interface or object, our algorithm, first, inserts the suitable nodes around it, in order to construct its mesh.Next, the convex combination of all stencils (31), as dictated by these nodes, is generated and the process is ready to be employed.It is important to stress that the scheme is implemented independently near any discontinuity, without the need to use the FVTD/PSTD hybrid technique [21,30].Thus, instead of applying (12) or (26) throughout such critical details, the piecewise polynomials of (30) along with the optimized stencils are utilized.This discretization leads to more robust and simple meshes, which do not generate any spurious modes at their interface with the FVTD/PSTD lattices, as elaborately explained in Section 4.
For illustration, let us inspect the two lossy dielectric regions of Fig. 2, where across the interface, both the tangential electric and magnetic field components are deemed continuous.In fact, the policy via which the PSTD method sets its variables [16] facilitates the simulation of the boundary condition and concurrently offers the sufficient means to handle the hybrid FVTD/PSTD area in a numerical sense.Specifically, region A is characterized by ε Α , μ Α , σ Α , ρ Α , while region B by ε Β , μ Β , σ Β , ρ Β , with σ and ρ denoting electric and magnetic losses, respectively.Also, the unit surface normal vectors pointing into the respective areas are depicted as n A and n B .
Setting n A = -n B = n, one obtains the n×E A = n×E B and n×H A = n×H B boundary conditions, that lead to If t indicates the tangential direction on the interface of the two media, (38) become ADVANCED ELECTROMAGNETICS, VOL.X, NO.Y, MONTH 2012 (a) (b) Figure 3: (a) A transverse cut of the hybrid interface for an orthogonal discretization between the complex-interface (prism region), the FVTD (light grey region), and the PSTD (white region) lattices.Dotted (black) lines notify evaluations via the complexinterface scheme, dashed (blue) lines indicate computations through the FVTD method, and straight (black) lines represent calculations via the PSTD technique.(b) A non-orthogonal dual-grid cell complex with the primary (blue) and secondary (red) mesh.
where S A = S A t and S B = S B t are specific area vectors defining contour C in Fig. 2. So, (39) can supply the required field quantities at almost any interface, even at the demanding ones in large-scale shielding/immunity or printed circuit board setups.

Conformal dual-mesh curvilinear modeling
The most convenient way to built our conformal lattice in a general coordinate system, properly selected to guarantee cells of optimal quality without face skewing or deformations, is the use of a non-overlapping dual-grid discretization [28][29][30][31].Actually, the precision of the proposed modeling algorithm relies on the systematic construction of the mesh in order to obtain the correct curves which enclose the object, under study, and the proper surfaces that surround its volume.
To this aim, we start with an adaptive primary lattice, whose refinement is performed according to the geometry of the involved structures and the attainment of the minimum dispersion error at all neighboring cell interfaces.This generalized concept offers a finite number of uniquely defined nodal patterns, which, if appropriately connected with their adjacent vertices, enable the formation of an equivalently well-posed secondary grid.The lattice, so developed, may comprise a variety of cells, like hexahedra or prisms, as indicatively shown in Figs 3a and 3b for an orthogonal and a general nonorthogonal discretization.This mesh is designated as primary and its vertices as primary vertices, where E and D vectors are typically positioned.Similarly, each face of a primary cell is called a primary face, while the center of the cell represents a secondary vertex.Two secondary vertices are connected if their two respective primary cells share a common face, so forming a set of secondary cells, namely the secondary lattice, where H and B vectors are located.For a non-planar face, the idea of the barycenter (i.e. point whose coordinates are the average of those of the face vertices), is incorporated [27,32].
Based on the above aspects, the leapfrog-like scheme of ( 21)-( 23) updates the unknown field variables of the problem as follows: (a) The PSTD update of H and E is conducted at time-steps n + 1/2 and n + 1 by their past PSTD E and H constituents at n and n -1/2, respectively.(b) The evaluation of H components at n + 1/2 near the interface requires the FVTD E values at time-step n.These values are supplied from the red-lined cell edges of Fig. 3a.(c) The rest of the components situated in the prism (dotted) lattice are explicitly advanced.Particularly, all E quantities at time-step n + 1 are obtained via past FVTD E components at n and n -1 along with those already derived at n + 1 via the FVTD method and are located on the blue-lined (dashed) cell edges of Fig. 3a.
In this context, the resultant hybrid technique employs a composite grid which consists of a rectangular mesh intersecting a conformal one.Essentially, for arbitrary objects or complex media interfaces, the primary lattice includes a small number of prisms perpendicularly erected along the objects' surfaces (Fig. 3a).This is exactly where the stencil-optimized procedure of Section 3 is applied, thus efficiently avoiding the use of the hybrid FVTD/PSTD method in this hard-to-model region.The height of the prisms is approximately the same as the PSTD spatial increment.In the internal primary or secondary vertices, the FVTD algorithm is used for time integration, with the interior boundary of the rectangular mesh set sufficiently close to fine geometric details.On the other hand, the outer-boundary E (or H) components on the FVTD primary lattice are acquired via interpolation of their computed PSTD counterparts, whereas E (or H) quantities on the PSTD mesh at the inner boundary of the rectangular grid are evaluated by interpolating the already calculated FVTD values.It is stated that this process, due to its local application, does not increase the memory or CPU burden to a noteworthy extent.

Performance and implementation aspects
A major aspect in the analysis of EMC structures with symmetries or periodicities is the separation of the domain into smaller regions, where a numerical scheme can be economically applied.Herein, to evade unnecessary discretizations of identical geometric areas, a domain decomposition algorithm is utilized.Our computational space is divided into N nonoverlapping subdomains.Specifically, for planar symmetries or periodicities, these regions include one unit cell.For instance, when two subdomains (i = 1,2) are considered, at their interface E 1 and H 2 (dually E 2 and H 1 ) vectors are related by Also, for the artificial boundary oscillations, an equivalent surface current is employed, i.e.
The proposed method along with ( 27)-( 37) is proven stable for all types of numerical simulations.Indeed von Neumann's analysis shows that ) sin (0.74) with the summations denoting a consecutive cyclic permutation of (u,v,w) and g h (u,v,w) the system metrics.Note the less strict character of (42) than the FDTD Courant condition [28], due to the dependence on orders K and M.This implies that time-steps can be chosen quite larger, without arising any oscillatory exponential modes.Alternatively, stability can be certified via the energy inequality method [5], which gives where a is a growth factor.Inequality (43) guarantees that G(t) does not surpass a certain limit namely the system's energy is conserved.Also the dispersion relation of the hybrid method is greatly improved due to the parameterized adjustment of its spectral resolution, as compared to the one, W FDTD , of the usual FDTD approach.Therefore, we derive To substantiate these merits, we present the maximum values of the L 2 error norm versus time in Fig. 4a and versus spatial increment in Fig. 4b for different K and M orders.
Obviously, the hybrid schemes exhibit a promising performance, as compared to two popular FDTD implementations.

Numerical results
The new method is verified by means of several real-world large-scale EMC applications, whose unbounded domain is truncated by a suitably adjusted 8-cell perfectly matched layer (PML) [33,22] with complex frequency shifting attributes.
Their selection has been mainly based on the dissimilarity of the involved media, fine details, electrical size, and modeling difficulty.The results, so extracted, are compared with those acquired by means of efficiently programmed (with regard to memory management) and very fine FDTD simulations as well as other popular techniques (used as reference solutions) or available measurements retrieved from the literature.Also, the precision of our FDTD realizations is confirmed via comparisons with two commercial packages, i.e. the Remcom XFdtd ® and the CST MWS ® for the same setup of all problems.Finally, all simulations are conducted on a 3.47 GHz Intel TM Xeon dual-processor X5690 computer with 96 GB RAM.It should be stressed, herein, that despite its large cell resolution and dense grids the FDTD method cannot adequately handle the aforesaid structural issues.Differently speaking, although one would anticipate that the FDTD scheme would accomplish acceptably precise outcomes for such ample computational setups, realistic problems reveal some important inaccuracies.These defects are mainly attributed to the high level of geometric details, the multitude of dissimilar media, the arbitrarily varying interfaces, and the large electrical size (in the order of 30λ to 100λ) of the EMC applications.Thus, the different spatial and temporal increments employed to treat these heavy details along with the significantly different constitutive parameters of the involved materials lead to detrimental dispersion and anisotropy errors as well as to serious late-time spurious artifacts.In this manner, the regular FDTD algorithm is not able to provide sufficient results, as would have occurred with other simpler configurations and similar grids.This is, in fact, the key motive for the development of the proposed hybrid methodology, namely to provide a general and flexible tool for the efficient alleviation of the preceding FDTD deficiencies in large-scale EMC arrangements.
The first structure is a 251 patch array, an indicative 21 element of which is given in Fig. 5, with its cross-polarized traits attained by two perturbed segments (0.65%) on every circular patch.Owing to its frequent use in many systems, the  optimal design of this large-scale array needs a trustworthy simulation means.Dimensions are: L x = 396 cm, L y = 14 cm, r = 18 mm, L = 26 mm, h = 0.19 mm, w = 3.1 mm, d = 2.8 mm and s = 6.2 mm.The grid, for the modeling regions of Fig. 5, comprises 2229462 ≈ 1.29 million cells (64.5 MB RAM for K = 3 and M = 3) in contrast to the much (around 88%) denser 670328260 ≈ 57.14 million-cell FDTD mesh and 2.05 GB RAM.Based on the first constraint of (13), we also get a 1 = 0.139874, a 2 = 0.256871, and a 3 = 0.103255 for (12).To certify the enhanced accuracy and economical profile of the new algorithm, Fig. 6 gives the return loss of the array.As detected, the accuracy of the FDTD simulations (confirmed by the Remcom XFdtd), even for the prior fine lattice is not satisfactory, unlike the results of the hybrid schemes which are in very good agreement with the reference solution (semianalytical integral-equation technique) of [34].This deduction can be justified by the maximum dispersion error (L 2 error norm) of the two techniques (FDTD: 2.8436×10 -1 ; proposed: 3.7109×10 -10 ), i.e. practically 9 orders of magnitude.
Next, we proceed to the coaxial waveguide of Fig. 7, which involves a set of microelectromechanical (MEMS) actuators for selective mode propagation.Its inner slab is based on a dielectric medium with ε r = 5.2 and the basic dimensions are: L = 410 mm, a = 32.3mm, and b = 58.4mm.The selection of this problem is primarily attributed to the increasing use of complicated MEMS in several high-end arrangements.To this aim, the variation of the S 21 -parameter magnitude versus the number of actuators is illustrated in Fig. 8.The specific setup with the modeling domain areas of Fig. 7 -especially for many actuators -leads to a fairly large-scale problem; e.g. for 24 actuators our grid has 134186232 ≈ 5.78 million cells (0.29 GB of RAM), while the corresponding FDTD one reaches the size of 342584648 ≈ 129.4 mil-   lion cells and almost 4.65 GB RAM.For operator (12), we select the values of a 1 = 0.260617, a 2 = 0.893701, and a 3 = -0.654318.Note that comparisons are conducted via the Sparameter theoretical formulation of [35] combined with the electromechanical model of [7], which is used as our reference solution.Moreover, Table 1 summarizes the resonance frequencies for four cases (with K = 4, M = 3 as our approximation orders) and several implementation details.It appears that, the proposed technique is very precise and cost-effective, attaining seriously diminished dispersion errors (up to 10 orders of magnitude) and sufficiently reduced simulation times.It is mentioned that the contribution of the complex-interface concept in the treatment of the intricate geometric details between the actuators and the slab, is critical.Remaining in the area of electrically-large microwave devices, our analysis moves to the coaxial four-port microwave splitter with a 3.23.10.7 cm dielectric (ε r = 25.2) region, depicted in Fig. 9.The outer conductor's cross-section is 5.22.4cm and the inner's 2.11.3 cm.Should we have selected the FDTD method, its lattice for a λ/80 resolution would consist of 720718456 ≈ 235.7 million cells (nearly 8.48 GB RAM).Nevertheless, our hybrid technique along with the complex-interface approach -used for the switches as in Fig. 9 -needs a total mesh of only 208214152 ≈ 6.76 million cells and almost 0.34 GB RAM.All simulation results are obtained for a 1 = 0.617981, a 2 = -0.902416,and a 3 = 0.215565.Figure 10 displays the behavior of S 13 and S 24 parameters with respect to the number of switches.Our reference solution is obtained from the modal analysis method of [34] along with the excitation concept of [28].As observed, the FVTD/PSTD algorithm overwhelms the FDTD method, since it successfully copes with the fine switches and the rather dissimilar dielectric region.This deduction can be, alternatively, confirmed by the maximum dispersion error of Figure 12: Geometry of a reconfigurable MEMS-controlled antenna on a double-negative metamaterial substrate.Modeling regions: complex-interface schemes for the radiating elements and the routers, optimized FVTD technique for the open region, and PSTD algorithm for the substrate.the two realizations, namely 9.6052×10 -2 for the FDTD and 1.3287 ×10 -10 for the proposed technique.Also, the required CPU time versus the number of switches is, again, proven to be very limited, as described in Fig. 11.
Let us, now, investigate a 120×116 monolithic antenna on a metamaterial substrate shown in its general form in Fig. 12.The directivity and main lobe of the structure are controlled by a set of MEMS power routers and phase shifters to obtain the necessary gain and radiation characteristics at a preselected frequency spectrum.For its discretization, the PSTD method is used for the substrate, whereas the FVTD algorithm is applied in the open region, thus leading to an overall mesh of only 246240122 ≈ 7.2 million cells (0.36 GB RAM).Special attention is drawn on the interfaces between the MEMS and the substrate, where the complex-interface process ( 27)-( 37) is applied.Notice that for (12), the values of a 1 = 0.679012, a 2 = 0.912748, and a 3 = -1.09176have been chosen.In this context, Fig. 13a gives the return loss of the antenna and Fig. 13b its radiation pattern for K = M = 4.For comparison, we derive the reference solution from the approximation models of [8] merged with the fast multipole schemes of [34,18].Again, the FDTD method (certified by the Remcom XFdtd) -requiring the rather large lattice of 832826450 ≈ 309.2 million cells and 11.2 GB RAM -does not manage to provide quite accurate results.In contrast, the hybrid algorithm is very precise for a 91% coarser grid.Similar observations can be drawn by inspecting the maximum dispersion error of the two approaches.Specifically, one re-Figure 14: A hemi-ellipsoidal H-slot dielectric-loaded antenna with two patches.Modeling regions: complex-interface schemes for the patches, the slot and the micro-strip line, optimized FVTD technique for the cover and the open space, and PSTD algorithm for the substrate.ceives 2.6701×10 -1 for the FDTD and 6.8934×10 -11 for the proposed method respectively.
Another large-scale structure is the hemi-ellipsoidal H-slot antenna with a 5.64 mm thick dielectric (ε r = 4.4) substrate and two patches for extra axial ratio (AR) and gain control (Fig. 14).The relative permittivity of the cover is ε rc = 6.2, while L = 16.5 cm, l 2 = 25.5 mm, and d 1 = 36.2mm.Applying the hybrid method with a 1 = 0.921742, a 2 = 0.706519, a 3 = -1.128261,we build a mesh of 196192186 ≈ 6.9 million cells (0.34 GB RAM). Figure 15 gives the 3-dB bandwidth versus AR frequency (f AR ) for two l 1 , d 2 sets and the radiation pattern.As discerned, our technique is proven more accurate, even with a far smaller (almost 89%) grid, than the FDTD one (716700596 ≈ 298.7 million cells and 10.8 GB RAM), whose precision is validated by means of the CST MWS.
Apart from the above applications, our investigation will also focus on several EMC test facilities, whose construction is very costly and thus any reliable simulation model is indeed necessary.The first is the inclined-wall dual reflector compact range anechoic chamber of Fig. 16 To model this large-scale facility, the proposed technique (for a 1 = 0.068923, a 2 = 0.175408, a 3 = 0.255669) requires a grid of 196198214 ≈ 8.3 million cells (0.42 GB RAM) for the increased orders of K = 7, M = 6, due to its very fine details. Figure 17a shows the deviation from open area test site (OATS) and Fig. 17b the normalized site-attenuation (NSA) of the facility's semi-anechoic version.Results indicate that notwithstanding its very large lattice (788792820 ≈ 511.7 million cells and 18.42 GB RAM), the FDTD approach (accuracy confirmed by the Remcom XFdtd) lacks to ensure measurement suitability in contrast to the novel algorithm which attains a significant 93% memory reduction.In this manner, we are able to correctly validate the chamber -namely full competence of performing reliable EMC or EMI measurements -prior to any actual construction attempt.
To improve the absorption performance of the compact range chamber, we introduce the generalized absorber linings of Fig. 18, which as observed consists of many dissimilar material interfaces.Hence, the use of the proposed methodology seems very appropriate.For our analysis, we select Κ = 7 and M = 7, with a 1 = -1.36827, a 2 = 0.455267, and a 3 = 1.413003.The grid has 202208216 ≈ 9.1 million cells (0.45 GB RAM), almost 85% coarser than the corresponding FDTD mesh with 804810846 ≈ 550.9 million cells and 19.83 GB RAM, for the minimum required resolution of λ/60. Figure 19a shows the reflectivity of several absorbers and Fig. 19b the electric field behavior in the quite zone for diverse chamber design parameters.The superior behavior of our absorber along with the accuracy of the hybrid technique is clearly visible.Furthermore, the maximum dispersion errors of the two techniques were found to differ around 10 orders of magnitude (FDTD: 4.89302×10 -1 ; proposed: 5.60214×10 -11 ).18b).The size of the pyramids behind the quiet zone follows a curved Chebyshev regime.For its precise discretization via our technique and a choice of a 1 = 0.183053, a 2 = 0.249421, a 3 = 0.067526, this complicated facility needs a 194198206 ≈ 7.9 million-cell grid and about 0.39 GB RAM.So, to characterize the chamber, field uniformity is depicted in Fig. 21a, which reveals the capability of the structure to perform measurements, unlike the FDTD method (precision validated by the CST MWS) with its 772784812 ≈ 491.46 million-cell mesh and practically 17.69 GB RAM.Similar deductions are drawn from Fig. 21b, where the effects of several chamber dimensions are thoroughly investigated.
Closing the category of test facilities, the hard-to-model reverberation chamber of Fig. 22 equipped with quadratic diffusers, i.e. gratings that diffract incident waves, is studied.).Herein, its movements are modeled with an interval of 22 o and a 60 temporal sampling rate.For this highly-detailed structure and K = 7, M = 8 (with a 1 = -0.870521, a 2 = 0.618902, a 3 = 0.751619) the proposed scheme leads to a mesh of 204216186 ≈ 8.19 million cells (0.41 GB RAM), while the FDTD method yields a 798812784 ≈ 508.1 million cells (92% finer and 18.29 GB RAM). Figure 23a illustrates the normalized E z variation, while Fig. 23b provides the shielding effectiveness of a 1.5 mm thick fiber-glass-fiber plate.It stressed that the reference solution has been acquired via the method of moments, described in [35].Moreover, the measurement procedure of [3] has been followed for the reverberation chamber of our laboratory to extract the necessary data for our comparisons.As detected, the FVTD/PSTD approach is superior to all FDTD realizations -verified in accuracy by the pertinent Remcom XFdtd simulations -in a fairly wideband sense.
For the case of a large-scale open-region problem, the validation procedure explores the thin shielded cable of Fig. 24, which connects two cylindrical cavities.Due to its low frequency nature, the discretization of this structure normally leads to very fine grids for the existing time-domain methods.Regarding its dimensions, a real-world scenario is:  The next application focuses on the analysis of a vehicleto-vehicle communication link coverage at the frequency of 5.9 GHz (Fig. 26a).Essentially, the main concept of this system Normalized Ez (dB) Frequency (GHz) Refence solution [35] Proposed lays in the exchange of information between moving vehicles or between vehicles and infrastructure in order to achieve high levels of efficiency and electromagnetic coupling (Fig. 26b).The antennas of the vehicles are vertically-polarized monopoles that transmit according to the IEEE 802.11p standard, i.e. with a power of 1.0 W EIRP (equivalent isotropic radiated power).For this very large-scale (56λ×30λ×100λ) problem, the overall mesh of our method comprises 192168446 ≈ 14.4 million cells (nearly 0.72 GB), while that of the FDTD technique requires 936×582×1528 ≈ 832.4 million cells and around 30 GB.Notice that regarding (12), a 1 = 0.221067, a 2 = 0.038512, and a 3 = 0.240421.Figure 27a presents the electromagnetic coupling between two vehicles versus their distance and compares the results with the reference solution acquired through the formulation of [36] and Fig. 27b shows the electric field distribution at the surface of the vehicle when the antenna is mounted on its roof.Again, the proposed approach offers a significantly more reliable performance than the FDTD one (verified by the CST MWS), despite its finer (almost 93%) discretization.
The final application is concerned with the wave propagation study of an indoor environment, where a wireless local area network (WLAN) operation has been installed.Our site is a typical 9.1×8.5×2.8 m office area with three independent rooms, as depicted in Fig. 28, whose walls are 20 cm thick and have a σ = 30 mS/m conductivity.Several potential positions (blue numbers in Fig. 29) of the transmitter -modeled as an omnidirectional radiating source -are considered, while both the 2.44 GHz and the 5.8 GHz case are explored.To implement the proposed method, we select a 1 = 0.031365, a 2 = 0.950086, a 3 = -0.418721,and construct a grid with 388342 192 ≈ 25.5 million cells (about 1.27 GB).The corresponding FDTD discretization for a resolution of λ/25 requires a mesh of 1768×1624×652 ≈ 1.87 billion cells and around 67.32 GB. Figure 29a illustrates the received (averaged) power levels along a diagonal path sketched in Fig. 28.Evidently, the hybrid schemes attain a very satisfactory accuracy, compared to the reference solution via the geometrical theory of diffraction [34] and the outcomes of the FDTD approach (validated by the Remcom XFdtd).Moreover, Fig. 29b and 29c show the Reference solution [36] power-coverage maps at 2.44 GHz and 5.8 GHz at source locations 2 and 4, respectively, and a dynamic range level from -40 to -15 dB.Observe the different penetration power levels through the walls of the site for the two frequencies and the smoothness of the plots achieved by our algorithm.Generally, interpreting the results from the applications studied in this section, it is derived that the conformal hybrid methodology can successfully cope with complex electricallylarge structures and realistic arrangements.It is emphasized that all problems can be safely deemed large-scale, since the corresponding FDTD meshes for their adequate discretization range from 55 million to 1.9 billion cells.As already deduced, the proposed technique leads to remarkably coarser (almost 90-93%) lattices and minimizes the undesired dispersion errors up to 10 orders of magnitude, without any exponential instabilities even for simulations lasting for several million time-steps.Consequently, high levels of accuracy are accomplished in contrast to the typical time-domain realizations, which require excessive numbers of cells without being able to provide reliable outcomes.Such assets are proven more significant as the application size becomes larger due to the noteworthy computational savings, thus enabling the analysis of demanding problems by usual limited-memory systems.

Conclusions
A hybrid time-domain methodology, based on advanced FVTD and PSTD schemes, has been proposed in this paper for the precise design and modeling of large-scale 3-D EMC structures.Established through a low-complexity generalized rationale, the parametric technique introduces new operators with optimal nodal density and media-sensitive interface conditions.In essence, the FVTD approach uses a leapfrog integration scheme -conditionally stable and nondissipative -that conserves the total charge in a global and local sense.Thus, canonical staggered-meshed formalisms (as in the FDTD algorithm) are efficiently circumvented with the computational overhead limited in reasonable levels.On the other hand, the PSTD concept uses Chebyshev polynomials, through a fast Fourier transform (FFT), instead of the customary finite-difference approximators, to evaluate spatial derivatives in areas of abrupt field variation.Due to its enhanced spectral accuracy, the PSTD technique entails only two cells per wavelength according to the Nyquist sampling theorem and, therefore, can deal with more challenging scenarios.Hence, and as realistic numerical evidence determines, artificial reflection errors are greatly suppressed and very rigorous as well as inexpensive solutions are acquired.

Figure 1 :
Figure 1: Transverse-plane view of a general discontinuity with constantly changing geometric features.

Figure 4 :
Figure 4: Variation of the maximum L 2 error norm (a) versus time in sec and (b) spatial increment in mm.

Figure 5 :
Figure 5: An indicative 2×1 element of the electrically-large 25×1 circularly-polarized wideband patch array.Modeling regions: complex-interface schemes for the corners and the segments, optimized FVTD technique for the patches and the PML, and PSTD algorithm for the rest of the domain.

Figure 7 :Figure 8 :
Figure 7: A two-port coaxial waveguide with a dielectric slab and several pairs of MEMS actuators.Modeling regions: complex-interface schemes for the actuators, optimized FVTD technique for the slab, and PSTD algorithm for the interior of the waveguide.Magnitude of S21 (dB)

Figure 9 :
Figure9: A four-port microwave splitter with a dielectric core and a number of switches.Modeling regions: complexinterface schemes for the switches, optimized FVTD technique for the dielectric and the inclined parts of the structure, and PSTD algorithm for the rest of its interior.

Figure 10 :
Figure 10: Magnitude of different S-parameters as a function of the switches number.

Figure 11 :
Figure 11: CPU-time required for the simulation of several splitter configuration as a function of the switches number.

Figure 13 :
(a) Return loss and (b) radiation pattern of the reconfigurable 120×116 antenna.

Figure 15 :
(a) 3-dB bandwidth versus axial ratio frequency and (b) radiation pattern for the hemi-ellipsoidal antenna.

Figure 16 :Figure 17 :
Figure 16: Geometry of an inclined-wall dual reflector compact range anechoic chamber.Modeling regions: complexinterface schemes for the antennas and the absorbers, optimized FVTD technique for the reflectors and the walls, and PSTD algorithm for the interior of the facility.

Figure 18 :
Figure 18: General structure of two enhanced absorber setups.(a) A curvilinear pyramidal multilayer lining and (b) a set of adjustable alternating wedges.

Figure 20 :
Figure 20: A double-horn semi-anechoic chamber.Modeling regions: complex-interface schemes for the absorbers, optimized FVTD technique for the walls, and PSTD algorithm for the interior of the facility.

Figure 22 :
Figure 22: A general mode-stirred reverberation chamber.Modeling regions: complex-interface schemes for the stirrer and the diffusers, optimized FVTD technique for the chamber's equipment and the walls, and PSTD algorithm for the interior of the facility.

LFieldLFigure 23 :
A /L D = 1.3, L B /L C = 0.72 L A /L D = 1.5, L B /L C = 0.58 L A /L D = 1.7,L B /L C = 0.63 L A /L D = 1.9, L B /L C = 0.44 Deviation from OATS (dB) , d E = 0.7 m) Twisted pyramids (FDTD, d C = 0.88 m) Prop.absorber (Hybrid, d C = 0.7 m) Prop.absorber (Hybrid, d C = 0.88 m) Prop.absorber (Hybrid, d C = 1.26 m) F /L B = 1.9, l F = 3.6 m L F /L B = 2.1, l F = 2.4 m L F /L B = 2.4, l F = 3.6 m L F /L = F = 2(a) Amplitude of the normalized E z component versus frequency in the interior of the chamber and (b) shielding effectiveness for a 1.5 mm thick fiber-glass-fiber plate.A set of dimensions is: l A = 6.22 m, l B = 12.84 m, l C = 6 m, l D = 3.2 m, l E = 5.8 m and l F = 0.8 m.A point where existing schemes fail is the stirrer's treatment (θ T = 34 o a = 36 m, b = 8 m, c = 12 m, d = 20 cm, d 1 = 1.2 m, d 2 = 0.45 cm, and d 3 = 0.7 cm.Also, a typical FDTD lattice for this configuration

Figure 24 :Figure 25 :
Figure24:A thin shielded cable.Modeling regions: complex-interface schemes for the cylinders and the cable, optimized FVTD technique for the cylinders' interior and the PML, and PSTD algorithm for the rest of the domain.Voltage (V)

Figure 26 :
Figure 26: (a) A vehicle-to-vehicle communication link scenario and (b) the 3-D vehicle geometry.Modeling regions: complex-interface schemes for the outer surfaces, optimized FVTD technique for the external space and the PML, and PSTD algorithm for the interior of the vehicle.

Figure 27 :
(a) Coupling between two vehicles versus their distance for two communication link cases (A: antennas mounted on the roofs and B: antennas mounted at the back of the vehicles).(b) Electric field distribution (in dB) at the outer surface of the vehicle for a roof-mounted antenna.

Figure 28 :Figure 29 :
Figure 28: Perspective view of the investigated office area.Modeling regions: complex-interface schemes for the inner walls, the windows, and the wooden structures, optimized FVTD technique for the external walls, and PSTD algorithm for the interior of the domain.

Table 1 :
Resonant frequencies of the coaxial waveguide