Accurate Nonstandard Path Integral Models for Arbitrary Dielectric Boundaries in 2-D NS-FDTD Domains

An efficient path integral (PI) model for the accurate analysis of curved dielectric structures on coarse grids via the two-dimensional nonstandard finite-difference time-domain (NS-FDTD) technique is introduced in this paper. In contrast to previous PI implementations of the perfectly electric conductor case, which accommodates orthogonal cells in the vicinity of curved surfaces, the novel PI model employs the occupation ratio of dielectrics in the necessary cells, providing thus a straightforward and instructive means to treat an assortment of practical applications. For its verification, the reflection from a flat plate and the scattering from a cylinder using the PI model are investigated. Results indicate that the featured methodology can enable the reliable and precise modeling of arbitrarily shaped dielectrics in the NS-FDTD algorithm on coarse grids.


Introduction
The nonstandard finite-difference time-domain (NS-FDTD) method is a high-accuracy analysis tool of the finite-difference (FD) type, proposed by Cole [1][2][3], to improve the performance of the FDTD algorithm [4,5].The scheme, based on an FD Laplacian featured via the nonstandard (NS) concept, introduced by Mickens [6], can reduce the overall propagation error of a typical FDTD implementation by a factor of 10 −4 on a coarse grid at a desired frequency [1,2].Through these excellent characteristics, the NS-FDTD method is particularly suitable for the accurate computation of the radar cross section (RCS) of electrically large objects.Amid the most prominent applications of the NS-FDTD method, one may discern the full-scale analysis of aircrafts-whose electrical size is very large, namely, more than 500 λ (λ is the wavelength)-with partially non-metal parts, including the radome, the canopy, the windows, and several radar-absorbing materials.
Nevertheless, the NS-FDTD method uses discrete space points to simulate electromagnetic wave propagation on orthogonal grids as in the FDTD technique [1][2][3][4][5].Therefore, the discrete space treatment is unsuitable for realistic problems with arbitrarily shaped objects and fine details, not aligned to the grid axes [4,5,[7][8][9][10][11], owing to the use of the insufficient staircase approximation on orthogonal grids in an effort to model the realistic object under study.Such structures can be, frequently, encountered in various applications, ranging from electromagnetic compatibility configurations [12][13][14] and microwave devices [15][16][17] to antennas [18][19][20], optical arrangements [21][22][23][24][25], and designs of low observability, including RCS scenarios.To circumvent such a drawback, a path integral (PI) model, based on the path integral form of Ampere's and Faraday's laws, has been previously presented [26][27][28].The specific path integral scheme can effectively handle real-world objects on orthogonal meshes, as the integral path may be intuitively placed on the object's surface.Thus, and owing to the PI model, the applicability of the NS-FDTD technique can be drastically expanded.However, the prior technique has been formulated only for the case of perfectly electric conductors (PECs); a fact that prohibits the analysis of various dielectric components found in many contemporary avionic and microwave telecommunication systems.
In this paper, a robust PI model for the rigorous treatment of dielectric objects, whose shape does not fit to any orthogonal grid, via the 2-D NS-FDTD method is developed.The proposed PI concept launches a simple scheme that considers the occupation ratio of all dielectrics at a specific cell.In this way, the manipulation of realistic dielectric structures becomes feasible on a coarse lattice.Firstly, to certify the key competence of the proposed PI model, we perform a reflection analysis of a flat plate and compare it with the exact solution.Next, an RCS problem is comprehensively explored via a real-world application.It is emphasized that for the NS-FDTD method, a practical grid width, ∆, is ∆ ≤ λ/8.On the other hand, if dielectrics are involved in the computational domain, the modified grid width should be ∆ ′ ≤ λ/(8 √ ε r ) (ε r is the relative dielectric permittivity).Furthermore, the use of a much finer lattice, i.e., ∆ ′ ≪ ∆, in order to satisfactorily treat the shape of a complicated structure, poses prohibitive issues in the implementation of the NS-FDTD method due to the unduly increase of the necessary system memory requirements and the extremely prolonged simulation times.These considerable drawbacks can be drastically circumvented by means of the proposed PI method, without the need for any additional conventions or non-physical assumptions, as occurs in existing schemes.Finally, all numerical outcomes reveal that the featured technique exhibits superior accuracy and convergence, as opposed to other computational approaches (with much finer spatial resolutions), for the trustworthy study of dielectric devices with intricate shape and elaborate geometrical features.

Basics and Formulation
First, let us review the concept of the PI model for the NS-FDTD method [26][27][28].The model uses both basic and complementary paths, as depicted in Figure 1, for the evaluation of the H z magnetic-field component.Herein, we intend to establish a complementary relation between these two paths.The idea is based on the wave propagation characteristics of the Yee algorithm [1,2,5], where maximum and minimum errors occur at an angle of 45 • on a square lattice.This indicates that we can mutually cancel these errors by means of two meshes for one node calculation.In an effort to attain the optimal cancellation of the inevitable discretization errors, the complementary (C) path is rotated 45 • with respect to the basic (B) one.Specifically, the two integral paths are expressed as for the basic path, with S = ∆ 2 (Figure 1a), and for the complementary path (Figure 1b).Note that in (1) and (2), we have replaced ∆ with the nonstandard correction function s k (∆) = 2 sin(k∆/2)/k, where k is the physical wavenumber [1,2].Based on this substitution, the propagation accuracy, (k n − k)/k, of the basic and complementary path model is shown in Figure 2 for two lattice resolutions and k n , representing the numerical wavenumber of the proposed PI model.The results are obtained via the method presented in [3,27], and as can be promptly observed, both models exhibit opposite propagation errors for the desired mutual cancellation.In this manner, the new PI model is derived by means of where, for homogeneous media, β 0 is given by [26] On the dielectric surface, we select the first term of (4), i.e., β 0 = 2/3, which is independent of the k, as the latter varies across the dielectric boundary.Additionally, in (3), we consider ∂H z /∂t = (H n+1 z − H n z )/s ω (∆t), with s ω (∆t) = 2 sin(ω∆/2)/ω, from the nonstandard formulation of [1,2], ω is the angular frequency, ∆t is the time step, and n = t/∆t.On the other hand, the E C y term in ( 2) is obtained from the already known E x,y values on the basic path as where e x,y are the corresponding unit vectors.Bear in mind that an analogous expression for the E C x term can be similarly obtained.

Treatment of the Magnetic-Field Components
The new PI model of ( 1) and ( 2) is, now, extended to the treatment of the definitely more demanding dielectric boundaries.Therefore, based on the geometric depiction of the proposed PI for such a scenario, the µ (∂H/∂t) • dS = − E • dl along the basic path in Figure 3 is given by where l 0,1 ≡ s k 0,1 (∆) = 2 sin(k 0,1 ∆/2)/k 0,1 with k 0,1 as the wavenumber of the ε 0,1 dielectric medium, S B0,B1 are the areas occupied by the ε 0,1 dielectric medium inside the basic path, and δ B1,B2 are path lengths along the basic path in the ε 1 medium.Also, the E x,y (x, y) ε 0,1 notation refers to the E x,y component in the ε 0,1 dielectric medium.Note that if E x,y is not available at a specific node on the grid, we use the E x,y extrapolation values at the nearest-neighbor node; for example, Similarly, the ∂H C z /∂t term using the complementary path in Figure 3 is given by It should be stressed that for the practical application of ( 7), we assume that the E C x,y (x, y) ε 0,1 terms exist on the integral path that straddles the dielectric boundary.In this manner, we can project the already known E x,y (x, y) values to the corresponding e x,y directional unit vector on the integral path.For illustration, the E C x (x − ∆/2, y + ∆/2) ε 0 term in ( 7) is obtained from the E x (x, y + ∆/2)e x • (e x + e y )/ √ 2 projection.Then, and similarly to (5), the E  6) and (7), respectively, the required ∂H z /∂t derivative is evaluated by means of (3).In this scheme, the path length is measured by correction function of ( )

Treatment of the magnetic field components
Here, we extend the basic PI model described by the formulas of (1a) and (1b) for the treatment of a dielectric boundary.The proposed PI model is shown in Figure 3.Then, the with ∆ , k0,1 the wave number, SB0,B1 and B0,B1 δ the occupied areas and path ratios in the 0,1 ε dielectric at the basic path region.Moreover, the 0 ,1 ( , ) x E x y ε notation re- fers to the Ex component in the 0,1 ε dielectric.If Ex is not available at a specific node, we use the projection by the inner product with the unit vector , x y e or the extrapolation values of the electric field at the nearest-neighbor node.Similarly, the term using the complementary path in Figure 2, is given by Finally, the t z H ∂ term is obtained from (2) via ( 5) and (6).

Treatment of the electric field components
Two typical cases at a dielectric boundary are shown in Figure 4. Specifically, in Figure 4(a), via with On the other hand, for the scenario in Figure 4(b), we derive S B0 S B1 Figure 3.The proposed PI model for the calculation of the H z component at a ∆ × ∆ lattice in the case of a dielectric boundary between an ε 0 and an ε 1 medium.The path length is measured via the s k (∆) nonstandard correction function [1,2].Moreover, S B0,B1 are the areas occupied by the ε 0,1 dielectric medium inside the basic path, δ B1,B2 are path lengths along the basic path in the ε 1 medium, and δ C1,C2 are path lengths along the complementary path in the ε 0 medium.

Treatment of the Electric-Field Components
The most frequently encountered cases for the computation of electric-field quantities at a dielectric boundary are those of a normal E x and a parallel E y component with regard to the boundary, as described in Figure 4.In particular, concentrating on Figure 4a and in terms of ε (∂E/∂t) for Moreover, for the scenario of Figure 4b, one acquires for ch p = s ω (∆t) On the other hand, when the case of an obliquely aligned, with respect to the lattice, dielectric boundary is taken into account, we introduce the cell division scheme of Figure 5, where A δ,ε is the occupied area ratio.Explicitly, Figure 5a is considered as the parallel case of Figure 4b, while the entire model in Figure 5a,b is treated as the normal case of Figure 4a.Hence, the corresponding ∂E x /∂t can be expressed as for ch a = s ω (∆t) with a similar formula holding for the ∂E y /∂t temporal derivative.In addition, the E C x,y components on the complementary path can be extracted through the inner products of these E x,y components with the respective unit vector u x,y = (e y ± e x )/ √ 2. Notice that ( 8)-( 13) are the best models, derived theoretically and numerically from the previous analysis.Finally, as a supplement, a special case of the right-angle edge can be treated by (i) selecting the ∆ width, so that the edge fits on the grid, (ii) moving the edge to the grid line, so that it fits on, and (iii) applying the modified PI models of Figures 4 and 5.   Furthermore, for an obliquely-placed boundary, we adopt the cell division scheme of Figure 5, where , A δ ε is the occupied area ratio.Explicitly, the in Figure 5(a) is considered as the parallel case of Figure 4(b), while the entire model in Figure 4(a) and 4(b) is treated as the perpendicular case of Figure 4(a).Hence, the corresponding t x E ∂ term

Reflectivity Analysis of a Flat Dielectric Plate
In order to substantiate the merits of the featured PI model, we focus on the reflectivity analysis of a flat dielectric plate as a basic example with an analytic solution.The computational domain is shown in Figure 6, whereas the PI models of ( 5)-( 13) are employed only on the dielectric surface.Moreover, the plate is infinite along the ±y directions, and an incident plane wave impinges from the right side in Figure 6.Therefore, the sine-cosine method [29] is used for our oblique incidence calculations.For our (TM = (H x , H y , E z ) and TE = (E x , E y , H z )) analysis, we examine two dielectric media, namely, an acrylic resin with ε 1 = 3ε 0 and a glass epoxy with ε 1 = 4ε 0 .The remaining implementation aspects are: λ = 1 m, ∆ = λ/20, and ∆t = T/30, with T as the corresponding wave period, while open-space truncation is conducted in terms of a 15∆-thick perfectly matched layer (PML) [30].To facilitate our comparisons, we use, as our reference, the analytical solution of the specific problem described in [31].In this context, Figure 7 illustrates the reflectivity of the flat dielectric plate for the aforementioned materials and two ∆w configurations, i.e., ∆w = 0 (the right-hand side of the plate surface is aligned to the grid) and ∆w = ∆/4 (the plate surface is not aligned to the grid).As observed, all the outcomes derived via the enhanced PI models are in promising agreement with the reference solution, thus proving the competence of the new scheme to effectively treat arbitrary dielectric boundaries.In this case, the flat plate model can be treated accurately for a geometrical integral path length and occupied ratio in the cell since the dielectric boundary is straight and parallel to the grid line (y-axis).Hence, we can detect that the PI model closely follows the analytic solution if we utilize an exact path length and occupied ratio in the PI cell.This fact proves the validity of the proposed PI model in the analysis of arbitrary dielectric boundaries.

RCS Analysis of a Dielectric Cylinder
Probing further, we proceed to a more realistic application and compute the RCS of an infinite dielectric cylinder with a radius of 5∆ and a relative permittivity of ε r = 3.For the curved boundary, we utilize the oblique dielectric models of Figure 5, formulated in ( 6)- (13).The grid layout is given in Figure 8, where the novel PI model is applied only on the cylinder surface and the rest of the domain is a regular NS-FDTD region.In particular, Figure 8a depicts the PI model of the H z component, expressed via (3), (6), and (7), while Figure 8b presents the PI model for the E x component, described by ( 8)- (13).This configuration is attributed to the occupation ratio of the dielectrics and the direction of the electric field in the shaded (cyan) cells.Note that the treatment of the E y component is the same as the E x one, owing to the rotational symmetry of the problem.Furthermore, for the excitation of the structure through an incident plane wave (impinging with an angle of θ = 45 • ) and the evaluation of the scattered waves from the dielectric cylinder, the total-field/scattered-field formulation [5] is employed.Finally, the scattered waves, so obtained, are converted to the appropriate RCS data by means of the near-to-far-field transformation technique [5,31].
In this framework, Figure 9 illustrates the comparative results between our model (PI model combined with the NS-FDTD technique: PI+NS-FDTD) and other approaches (typical FDTD and NS-FDTD method without the PI model) for λ = 1 m, ∆ = λ/14, and ∆t = T/20.Also, the outcomes of a much finer FDTD implementation (∆ = λ/84, ∆t = T/120) serve as our reference solution.It becomes apparent that the PI+NS-FDTD method agrees very well with the reference data within the range of 0.5 dBm, which is, actually, a negligible difference in most real-world RCS configurations.Herein, the reduction effect is (14/84) 2 for the overall computer memory and the same analysis space size since the reference FDTD solution uses a ∆ = λ/84 and the NS-FDTD simulation with the PI model employs a ∆ = λ/14.Additionally, our method achieves a reduction effect of (14/84) 2 × (20/120) for the total CPU time, compared with the FDTD algorithm, because the reference FDTD solution uses a ∆t = T/120, while the NS-FDTD simulation with the PI model employs a ∆t = T/20.In Figure 9, a small fluctuation of 0.5 dBm is observed in the PI results for the observation angle θ.This small discrepancy can be attributed to the geometrical layout of the PI cell in Figure 8, as the PI models periodically vary regarding the cylinder surface position.To mitigate this RCS fluctuation, one may use finer cells, as deduced from the FDTD results with ∆ = λ/14 and ∆ = λ/84, since the specific PI model robustly replaces the complicated scattering phenomena on the cylinder surface (smaller than ∆ × ∆) with one cell calculation on the coarse grid.However, it should be emphasized that the RCS fluctuation versus the observation angle is much smaller than that of the NS-FDTD simulation without the PI model.Lastly, regarding the bistatic RCS for θ = 45 • , our scheme offers, again, a very satisfactory performance, as shown in Figure 10, for the same set of implementation parameters as in the prior example.It should be stated that only a small RCS fluctuation is detected in the PI outcomes.Bear in mind that if a sufficiently fine (i.e., ∆ = λ/20) cell is selected in Figure 8, such fluctuations are not at all observed.This implies that the fluctuations in Figure 10 are attributed to the use of the ∆ = λ/(8 √ ε r )| ε r =3 ∼ λ/14 maximum grid size [1,2] and the geometrical layout of the PI cell, as already mentioned in the monostatic case.Consequently, the novel formulation can reliably and accurately treat the curved shape of dielectric structures, even when their surface is not aligned to the axes of orthogonal NS-FDTD grids.8)- (13).The proposed scheme is used solely on the cylinder surface, while the remaining domain is a typical NS-FDTD region.Note that due to rotational symmetry, the E y case is similar to the E x one.Finally, in the case of the bistatic RCS, our scheme offers, again, a very satisfactory per-232 formance, as shown in Figure 10.Consequently, the featured formulation can reliably 233 and accurately treat the curved shape of dielectric structures, even when their surface is 234 not aligned to the axes of orthogonal NS-FDTD grids.

Conclusions
Conventional grid-based computational methods exhibit significant precision issues since they are not able to efficiently model curved media boundaries.To alleviate these shortcomings and accurately handle dielectric objects of an arbitrary shape, not aligned to the grid axes, a consistent and efficient PI form for the 2-D NS-FDTD technique has been presented in this paper.The key idea of the proposed PI scheme is to utilize the occupation ratio of dielectrics in the necessary cells, surrounded by basic and complementary paths, thus attaining a notable reduction of the overall system burden.Conducting the numerical RCS analysis of a dielectric plate and cylinder, it has been validated that the new concept improves the accuracy of the regular NS-FDTD algorithm in the case of dielectric boundaries, considerably curved interfaces, and fine geometric characteristics.Thus, the universality of the NS-FDTD technique can be seriously expanded, even on coarse grids, without having to resort to any additional complicated approaches or artificial assumptions.Specifically, for the demanding RCS study of electrically large configurations (partially occupied by various dielectrics), which opt for excessive memory requirements and elongated simulations, the new PI model is proven to offer remarkable computational savings.Nonetheless, it should be stated that generating the appropriate geometrical data for the PI model is a laborious task from a calculation point of view.Therefore, for a beneficial real-world use of the featured scheme, the development of a CAD tool is deemed necessary.

Figure 1 .Figure 1 .Figure 1 .Figure 2 .
Figure 1.The two paths involved in the proposed integral NS-FDTD form on a square grid for the 49 Hz magnetic field component.The two paths are rotated 45 o , one with respect the other, in order 50 to mutually cancel their errors.5152

Figure 3 .
Figure 3.The PI model at a ∆× ∆ lattice in the presence of a dielectric boundary.

Figure 4 .
Figure 4. Typical dielectric boundary cases for Ex and Ey components.

Figure 5 .
Figure 5. Division model of an Ex cell in the presence of an oblique dielectric boundary.

Figure 4 .Figure 4 .
Figure 4. Treatment of a (a) normal E x component and (b) parallel E y component with respect to the dielectric boundary.Here, δ is the occupation ratio equal to the area of the ε 1 medium over the entire area of a cell.

Figure 5 .Figure 5 .
Figure 5. Division model of an Ex cell in the presence of an oblique dielectric boundary.

Figure 6 .Figure 7 .
Figure 6.Reflectivity evaluation of a flat dielectric plate by means of the proposed PI model, where, except for the PI area, the rest of the domain is a differential type NS-FDTD region.Note that ∆w is the distance between the grid line of the basic path and the plate edge, with ∆w = 0 at the left side.The dielectric plate is infinite along the ±y directions with a thickness of λ/4 + ∆w.

Figure 7 .
Figure 7. Reflectivity (TM =(H x , H y , E z ) and TE = (E x , E y , H z )) analysis versus the incidence angle,θ, for the flat dielectric plate of Figure6, performed in terms of the proposed PI model and compared to the analytical (reference) solution of[31].The analysis involves two dielectric media (acrylic resin and glass epoxy) in conjunction with two ∆w arrangements (∆w = 0: the right-hand side of the plate surface is aligned to the grid and ∆w = ∆/4: the plate surface is not aligned to the grid).(a) Acrylic resin with ε 1 = 3ε 0 and ∆w = 0, (b) glass epoxy with ε 1 = 4ε 0 and ∆w = 0, (c) acrylic resin with ε 1 = 3ε 0 and ∆w = ∆/4, and (d) glass epoxy with ε 1 = 4ε 0 and ∆w = ∆/4.

Figure 8 .
Figure 8.The PI model of an infinite dielectric cylinder (first quadrant) for the TE mode with the positions of the involved (a) H z components, computed through (3),(6), and(7), and (b) E x components, calculated via (8)-(13).The proposed scheme is used solely on the cylinder surface, while the remaining domain is a typical NS-FDTD region.Note that due to rotational symmetry, the E y case is similar to the E x one.

Figure 9 . 225 Figure 9 .
Figure 9. Monostatic RCS results derived from the PI model and compared to other 221 schemes.In a practical RCS design, a difference of 0.5 dB is negligible, namely, the 222 PI+FDTD outcomes practically agree with the reference data.223224 225

Figure 9 .Figure 10 .
Figure 9. Monostatic RCS results derived from the PI model and compared to other 221 schemes.In a practical RCS design, a difference of 0.5 dB is negligible, namely, the 222 PI+FDTD outcomes practically agree with the reference data.223224 225

Figure 10 .
Figure 10.Bistatic RCS results derived from the proposed (PI+NS-FDTD) technique and the NS-FDTD algorithm without the PI model.The reference solution is acquired from an FDTD realization with a very fine (∆ = λ/84, ∆t = T/120) lattice resolution.The incidence angle of the impinging plane wave is 45 • , while the geometric PI cell data of the cylinder are depicted in Figure 8.