A Nonstandard Path Integral Model for Curved Surface Analysis

: The nonstandard finite-difference time-domain (NS-FDTD) method is implemented in the differential form on orthogonal grids, hence the benefit of opting for very fine resolutions in order to accurately treat curved surfaces in real-world applications, which indisputably increases the overall computational burden. In particular, these issues can hinder the electromagnetic design of structures with electrically-large size, such as aircrafts. To alleviate this shortcoming, a non-standard path integral (PI) model for the NS-FDTD method is proposed in this paper, based on the fact that the PI form of Maxwell’s equations is fairly more suitable to treat objects with smooth surfaces than the differential form. The proposed concept uses a pair of basic and complementary path integrals for H -node calculations. Moreover, to attain the desired accuracy level, compared to the NS-FDTD method on square grids, the two path integrals are combined via a set of optimization parameters, determined from the dispersion equation of the PI formula. Through the latter, numerical simulations verify that the new PI model has almost the same modeling precision as the NS-FDTD technique. The featured methodology is applied to several realistic curved structures, which promptly substantiates that the combined use of the featured PI scheme greatly improves the NS-FDTD competences in the case of arbitrarily-shaped objects, modeled by means of coarse orthogonal grids.


Introduction
The nonstandard finite-difference time-domain (NS-FDTD) method is an accurate simulation tool having excellent isotropic propagation characteristics for the numerical electromagnetic wave analysis at a desired frequency [1][2][3][4][5][6][7].Since its accuracy is 4 10 times higher than that of the FDTD technique [8,9], the NS-FDTD approach is suitable for the precise electromagnetic design of objects with electrically large sizes, such as aircrafts comprising various materials, electronic components, and complex airframe shapes [10][11][12][13].In general, a typical aircraft size is around 500 1800    at the radar frequency.For such dimensions, the phase error of the FDTD method is 7 23   approximation, just as in the typical FDTD method.Nonetheless, it is common knowledge that the optimal design of a realistic arbitrarily-curved structure necessitates the meticulous modeling of every geometric peculiarity, as well as the detailed characterization of any media interface.When these attributes are to be explored by a staircase time-domain method, such as the NS-FDTD one, the requirement of credible discretizations cannot be, easily, fulfilled, because of its inherent restriction to orthogonal lattices.As a consequence, material boundaries, not aligned to the grid axes, are poorly or even erroneously approximated.To evade these difficulties, an assortment of proficient techniques has been proposed (and summarized in [8]), ranging from curvilinear or nonorthogonal FDTD variations to modified Cartesian and conformal algorithms.Nonetheless, the research on this significant topic is constantly escalating, thus leading to various schemes which offer treatment to many contemporary applications from the microwave and optical regime [19][20][21][22][23][24][25][26][27][28][29][30][31] and paving the way to the trustworthy analysis of many future state-of-the-art scientific fields [32][33][34][35][36][37][38][39][40][41][42][43].[6] 2D-3D/staircased Square/cubic cell  2D/3D rectangular cell 4.5 2015 [16] 4D/staircased Cubic cell  4D subgrid model 15 2019 [17] 2D/curved Square grid  2D PI model 100 2020 [18] 3D/curved Cubic cell  3D PI model 216 Therefore, the manipulation of arbitrarily-curved surfaces or material interfaces via a staircase model deteriorates the reliability of any FDTD-based approach and, particularly, of the NS-FDTD scheme, whose stencils must be very carefully selected.Apparently, the use of finer lattices cannot alleviate this serious hindrance due to the excessive system requirements and the undue increase in the overall computational burden.On the other hand, the use of subgridding techniques [8,9] could yield potential solutions, however, such schemes (in their original form) are often prone to long-time instabilities or opt for very small time increments.Such issues, and particularly the latter, may be proven prohibitive for the NS-FDTD method in the case of prolonged simulations.Indeed, the Courant-Friedrich-Levy (CFL) criterion that controls the stability of both the FDTD and the NS-FDTD algorithms has been, in spite of its worth, one of the most limiting factors for the modeling of electrically large structures.Therefore, to model the tiny-scale features of curved-shaped surfaces, such as those of an aircraft, the cell size selected should be adequately small, enforcing an upper threshold on the time-step that, in turn, leads to thousands of time iterations for the steady-state.To overcome this shortcoming, the alternating-direction implicit (ADI) [44][45][46][47][48][49][50][51] and the locally one-dimensional (LOD) [52][53][54][55][56][57] time marching concepts have been introduced in the update mechanism of most FDTD-oriented schemes.The resulting approach is a semi-explicit technique which surpasses the CFL limit and can conduct unconditionally stable simulations, on the condition that the numerical dispersion and dissipation errors are carefully taken into account.
Considering the prior aspects and the theoretical establishment of the NS-FDTD method, a promising candidate for the analysis of curved surfaces could be the contour-path (CP) model, based on Ampere's and Faraday's laws [8].The CP scheme can properly treat curved details, not necessarily aligned to an orthogonal grid, because the path integral can adapt to any geometrical shape.Another advantage is the use of coarse grids for the handling of curvatures, which seriously decrease the computational overhead.In fact, since its initial advent, the CP model has notably increased the applicability of the FDTD method.Nonetheless, in the NS-FDTD technique, the CP scheme has not yet been introduced.Such a combination could be, indeed, proven to be of critical importance, since the development of the optimal path's integral model for the NS-FDTD algorithm is necessary to enhance its accuracy in real-world applications.
In this paper, efficient two-dimensional (2D) and three-dimensional (3D) path integral (PI) forms are introduced for the NS-FDTD method, to facilitate the CP modeling of smooth-surface objects.The new PI model launches two types of integral paths on a square grid for each calculation node [16,17].These path integrals are optimized in terms of certain parameters via the technique's numerical dispersion equation.Extensive numerical simulations verify that the PI model has a superior isotropy and calculation accuracy (almost equivalent to those of NS-FDTD) compared with the conventional FDTD method.Specifically, practical applications involve the 2D scattering analysis of a circular perfectly electric conductor (PEC) cylinder and the study of a 3D antenna via a thin wire approximation, accordingly tailored for the requirements of the PI model.These demanding numerical investigations promptly reveal that the featured PI model, acting as a CP scheme, as in the original FDTD approach, can decisively improve the modeling competences and accuracy of the NS-FDTD method, even by means of very coarse-grid discretizations.

Development of the Path Integral Model
Let us consider a PI model on a square grid, with a cell width  , that offers the highest accuracy for its combined use with the NS-FDTD method.Since, in the conventional FDTD approach, a propagation error arises along the directions of 0  and 45  on the grid [8], we introduce two distinct integral paths (one rotated 45  with respect to the other) that can mutually cancel their errors, as depicted in Figure 1.In this context, for the basic path (Figure 1a) on a square grid, the integral form of Maxwell's equations for the z H magnetic field component is given by while, for the complementary path (Figure 1b),

H S E l
(2) To guarantee the desired isotropy of the wave propagation on the selected grids, the integral forms of Equations ( 1) and ( 2   where 2D 0  is a parameter, derived in the next subsection, that can finely adjust the overall accuracy of the resulting explicit expressions.Essentially, Equations ( 1   with analogous expressions holding for the x E and x E  components can, as well.Re- garding the PI forms of Equations ( 1)-( 5), the cell width  , the time derivative t g  , and the time increment t  are substituted by   ,respectively, according to the NS-FDTD approach [2][3][4].Further, k is the wave number and  the angular fre- quency.That is, for (3), while the PI expressions of Equations ( 4) and (5) Similarly, for the x E and x E  components, we arrive at 3), let us derive the numerical dispersion equation.Firstly, the monochromatic traveling wave of 6)-( 10), where 2  1 j   and r are the position vector, through the analysis of [9].After some calculus, we define the matrix of the coupled system of equations, for T , , , ,

Derivation of Parameter
with Then, the desired numerical dispersion equation is retrieved by imposing the Since, for k    , parameter 2D 0  receives its optimal value for every propagation angle  , we substitute the specific k into Equation ( 15) and solve it with regard to which, after the use of the appropriate Taylor series for the sin terms, provides the compact and convenient formula of

Numerical Accuracy of the 2D PI Model
Taking into account Equation ( 15), we numerically examine the accuracy of our PI model with regard to parameter 2D 0  of Equation (17).The calculated numerical wavenumber, num k , is shown in Figure 2, where the overall error is normalized to the theoret- ical wavenumber, k , in free space, as   . In particular, Figure 2a compares the FDTD method, the NS-FDTD technique, and the new PI model, whereas Figure 2b presents a detail of the previous comparison between the NS-FDTD method and the featured PI model.In all simulations, we have selected a wavelength of 1   .As promptly ob- served, the propagation accuracy of our PI model is almost equivalent to that of the original NS-FDTD technique and almost 4 10 times higher than that of the FDTD approach [8].Consequently, it can be deduced that the PI model is satisfactorily precise to connect with the NS-FDTD scheme.It is stressed that in the next section, we demonstrate the usefulness of the PI scheme, in terms of the prior property, in NS-FDTD regions, so that it can be successfully applied to the treatment of arbitrarily-curved surfaces. .Note that the NS-FDTD method and the PI model have practically the same accuracy, which is almost four orders of magnitude higher (or even more) than that of the FDTD method.

Numerical Verification of the 2D PI Model
In this section, we investigate the efficiency of the proposed 2D PI model.Thus, as an example, let us consider the radar cross section (RCS) of an infinite circular perfectly electric conductor (PEC) cylinder.The PI model, applied on the arbitrary surface, is designated as CP (contour path) in the FDTD algorithm.Therefore, we also use the CP term to refer to the modified PI model, which conforms to curved-surface structures.

CP Formulation through the Proposed 2D PI Model
For the realization of the PI scheme, in order to effectively model arbitrary surfaces, the modified basic and complementary paths are depicted in Figure 3. Specifically, the integral area for the computation of the z H magnetic field component is enclosed by two paths, i.e., the B S and C S paths, apart from the PEC cylinder area.It should be mentioned that the treatment of the prior CP integrals is based on the same principles as those of the FDTD algorithm.Hence, the integral CP form for the basic path in Figure 3a becomes where Actually, this process is the same for the case of the complementary path with an area of C S in Figure 3b.The treatment for the corre- sponding z H formula is also similar to those of Equations ( 3) and ( 6), whereas, at the final step of Equation ( 18 on an integral path which cannot be manipulated via Equations ( 7)- (10).In such a case, the projection or extrapolation values, obtained from the already extracted E field components at the nearest-neighbor nodes, are utilized.Lastly, the remaining E and H components are derived through the PI Equations ( 6)-( 10).4) and ( 5) is based on the same principles as those in the FDTD algorithm.

Scattering Analysis of a Circular PEC Cylinder
Based on the above aspects, the developed CP model is validated through the RCS evaluation of an infinite circular PEC cylinder with a radius of / 2  .To this aim, Figure 4 shows the entire computational domain, where the excitation is a TE plane wave with 1   m and an incident angle of inc 0    .The scattered field is acquired through the total-field/scattered-field (TF/SF) formulation [8] in the NS-FDTD region.A significant asset of the featured PI concept is its ability to reliably employ coarse grids.Therefore, the CP model uses a spatial increment of /12    for the basic path and a /8 2    for the complementary path.Furthermore, the temporal increment is set to . For comparison, the corresponding NS-FDTD simulation is also conducted in the absence of the combined CP and PI model.Figure 5 presents the RCS of the cylinder, where the FDTD results, computed via a staircase model with the much finer resolution of /120 , are also given as a reference.Analyzing the outcomes, it is found that our concept agrees very well with both the reference and the NS-FDTD solution.On the other hand, the accuracy of the NS-FDTD method with the resolution of /12    is much lower than that of the featured model. , concerning the validation of the 2D CP technique applied to the proposed PI concept in the NS-FDTD area terminated by a perfectly matched layer (PML) [58].The structure is excited by a TE plane wave with 1   m and an incident angle of inc 0    , while all scattered-field components are obtained by means of the TF/SF formulation [8] in the NS-FDTD region./ 000  , regarding the time steps, of the respective non-CP case.Consequently, the new PI concept notably enhances the NS-FDTD calculation accuracy, convergence, and overall efficiency for the treatment of curved objects.The small fluctuation of the RCS in Figure 5 depends on the geometric precision of the CP cells at the cylinder surface, as shown in Figure 3, such as the path length and the integral area.This implies that in order to avoid such minor issues, the CP cell model of the curved surfaces must be meticulously designed and implemented.

Development of the Path Integral Model
To derive the 3D extension of our PI concept and cancel the propagation error due to angle dependence, we adopt one basic and one complementary integral path at the z plane and two complementary integral paths at the z   planes (one rotated 45  with respect to the other), as shown in Figure 6.Therefore, for the basic path at the z plane (on a square grid), the integral form for the z H magnetic field component is given by whereas for the complementary paths at the z and z   planes, we obtain the form of Figure 6.The four integral paths for the calculation of the z H magnetic field component on a square grid (similarly for the x H and y H fields). Specifically, we devise one basic and one com- plementary path at the z plane and two complementary paths at the z   planes in order to cancel the propagation error owing to angle dependence.
In order to accomplish high isotropy, we use four integral forms, such as those provided in Equations ( 19) and ( 20), and derive with 0,1  as the parameters for the 3D isotropic wave propagation, whose optimal values are derived in the next subsection.Actually, Equations ( 19)-( 21) constitute the featured 3D PI form.
Probing further in Equation ( 20) and through Figure 6, y E  may be computed in terms of at the z and z   planes.Likewise, x H  , in Equation ( 22), is given by where the cos( /4)  term signifies the projection along the x E  direction and the 0.5 re- fers to the average of the x H or y H magnetic field component at the /2 z   planes, for instance,  .Note that a similar procedure, such as ( 19)-( 23), holds for the evaluation of the x H and y H components. On the other hand, regarding the x E and y E expressions, we employ the integral path concept on a 2   lattice.The final PI formula is acquired by plugging Equations ( 19), ( 20), (22), and ( 23) into (21).

Derivation of Parameters 0
 and 1  For the optimal 0  and 1  values, let us derive the numerical dispersion equation.Hence, using the trial wave solution in [8,9], Equation ( 20) leads to with an analogous expression derived for Equation (19).Moreover, from Equations ( 19)-( 23), we define A as the matrix of the coupled equation set for , , , , , , , , , , , , where superscripts , , xy yz and xz de- fine the plane where E  exists.In this manner, the required numerical dispersion equation [8] of our 3D PI model can be extracted from with coefficients , , , elaborately provided in Appendix A.
Subsequently, we determine 0  and 1  by considering a plane wave propagating along the 90    direction at the xy plane.This case is, essentially, the same as the 2D PI model of Equation (3), i.e., Therefore, comparing Equation (21) with Equation ( 26), we arrive at Substituting Equation ( 27) into Equation ( 25) and setting , /4     , as a typical propagation angle, one gets in which coefficients while, from Equations ( 17), (27), and ( 29), 0

Numerical Accuracy of the 3D PI Model
The propagation accuracy of the novel 3D PI model combined with the NS-FDTD technique is investigated through the numerical dispersion of Equation ( 25), and the results are illustrated in Figure 7, where the analytical 1  values are obtained by means of Equation ( 29) as .Note that  is the grid size, k the physical wavenumber,  the wavelength, and ,   the angles of the spherical coordinate system   , , r   .For comparison, the typical NS-FDTD results [6] are also included.
It should be mentioned that, after the appropriate study, we found that four mantissa digits were enough for the precise representation of 1  .Furthermore, for the sake of comparison, the regular NS-FDTD outcomes are also plotted in Figure 7, with num k extracted through the Newton-Raphson method applied to Equation (25).As can be detected from Figure 7, the proposed PI model offers a high propagation accuracy that is almost the same as the one of the NS-FDTD algorithm toward the axial direction of 0    .In this way, the propagation error of the NS-FDTD method is in the order of 7 10  -8 10  , as presented in Figure 8.Moreover, Figure 9 presents the propagation error of the FDTD method combined with our PI model.Even though the error of the PI model slightly increases at 45    , it is still three orders of magnitude lower than that of the FDTD method for the same  .Regarding the error fluctuation at 45    in the PI results, we can deduce that the error compensation of the PI model was not so sufficient at the z   planes.Thus, it is proven that the new 3D PI model overwhelms the typical NS-FDTD method in terms of isotropic propagation.Note that in this paper, we use the stability condition of /( 3) , since the 3D PI model uses  -and 2 -size cells at the same time, as shown in Figure 6.times lower than that of the FDTD method for the same grid size  .

Numerical Validation of the 3D PI Model
For the certification of the 3D PI model, we analyze the input impedance of a thin wire dipole antenna using the 1/ r approximation approach [8,59].

Study of the Thin Wire Antenna Model
Firstly, we apply the 1/ r thin wire model for the FDTD method to the 3D PI model.The details of the wire antenna given in Figure 10, where the solid black line refers to the basic path and the dashed blue line to the modified complementary path.Moreover, SB is the area enclosed by the basic path and SC the respective area for the complementary path.

In this framework, the integral area for
in Figure 10b; SB = 2 a    in Figure 10c; and SB = 2  and SC = 2 2 7 /4 a   in Figure 10d.In particular, for the complementary path in Figure 10a, where x y

  t = e e
and , x y e e are the unit vectors.Since, in the above expressions, On the other hand, for the complementary path in Figure 10b,d, we get     /2, /2 , /2 cos( /4) where Additionally, for the basic path in Figure 10a, it is derived that with a similar treatment holding for the integral paths with  being the electric conductivity and D J the driving current.

Dipole Antenna Analysis
Bearing in mind the above analysis, we study the impedance of a dipole antenna of length L , using the modified 1/ r models for the featured PI scheme, as in Figure 10.The antenna is placed at the center of the PI region of width 10 around the wire.Outside the PI area, we have a 200 200 300      NS-FDTD region, terminated with a 20 -wide perfectly matched layer (PML) [58].Figure 11 presents the antenna impedance versus the wire length, with 1 , for 41 N  and 81 .For the selection of the time increment, we employ the stability criterion of the FDTD method.In addition, the antenna is driven by a center feed with a gap  along the z axis, as shown in Figure 10d, while the evaluated wire radii are 0.001 a  m and 0.002 m.As our reference in Figure 11, we use the results retrieved via the Numerical Electromagnetics Code (NEC2) computational package [60], which implements the method of moments [61], whereas the corresponding FDTD results [8] are also provided.Note that the antenna impedance, ant , Z R jX   at the feed point is calculated by referring to Figure 10d [59].In Equation ( 44), the denominator gives the contour integration of the magnetic field H around the feed point.It can be discerned from Figure 11 that all the outcomes derived via the PI model are in very good agreement with the NEC2 outcomes.Moreover, the PI concept achieves higher accuracy than the typical FDTD method.As a consequence, the modified thin wire antenna model based on the PI scheme is proven to be very efficient, unlike the original NS-FDTD technique in differential form, offering the opportunity for the precise analysis of demanding thin wire structures.m.Due to the limitations of the 1 / r approximation model [8,59] on the grid, the extracted results are presented for /2 a  , apart from those derived via NEC2 [60].

Conclusions
A novel PI formulation for the rigorous treatment of curved surfaces on square grids via the 2D and 3D NS-FDTD method has been presented in this paper.In essence, the developed model has been extended to a robust CP scheme, which can enhance the efficiency of the main numerical algorithm, even with coarse mesh tessellations.After the pertinent numerical investigation, it has been confirmed that the PI model exhibits a comparable accuracy and isotropy to the NS-FDTD technique.Moreover, all results obtained from 2D RCS calculations and 3D thin-wire antenna studies substantiated that the PI implementation is fully compatible with all the interfaces of the NS-FDTD domain.Therefore, taking into account that the proposed methodology avoids excessively fine lattices, it can be deemed as a precise and reliable means for modeling electrically-large objects with complex curved surfaces or media interfaces, frequently encountered in the design of modern electromagnetic configurations.A14)

Figure 1 .
Figure 1.Integration paths on a square grid for the proposed PI model.(a) Basic path and (b) complementary path.The two paths are rotated 45  , one with respect the other, in order to mutu- ally cancel their errors.Thus, the nodes for the coefficients of the NS-FDTD operators [2-4] are not required.

Figure 2 .
Figure 2. Numerically computed wavenumber, num k , versus propagation angle  for the FDTD method, the NS-FDTD technique, and the proposed 2D PI model.(a) General comparison and (b)detail of the previous plot only between the NS-FDTD method and the PI model.In the analysis, k is the physical wavenumber, the wavelength is chosen to be 1   ,  is the grid size, and ),  is replaced with

Figure 3 .
Figure 3. Modified integral paths for the proposed PI model based on (Figure 1a) and (Figure 1b), respectively.(a) Basic path and (b) complementary path.Dashed lines represent the position of the pair paths, PEC areas are indicated in light blue, and , path length / x y l l   signifies the ratio of the path length over  .Note that the integral path and the area for the calculation of the z H magnetic field component form a part, except for the PEC region, indicated at the right side of the figure.The treatment of the path integrals in Equations (4) and (5) is based on the same principles as those in the FDTD algorithm.

Figure 4 .
Figure 4.The computational domain for the RCS analysis of the infinite circular PEC cylinder with a radius of / 2 , concerning the validation of the 2D CP technique applied to the proposed PI concept in the NS-FDTD area terminated by a perfectly matched layer (PML)[58].The structure is excited by a TE plane wave with 1

Figure 5 . 2 (
Figure 5. RCS results obtained via the FDTD, the NS-FDTD, and the proposed 2D CP model for the computational configuration of Figure 4.The incident angle of the TE plane wave is inc 0 

Figure 7 .
Figure 7. Propagation error of the numerical wavenumber, num k , obtained via the Newton-Raphson method using Equation (25), for the NS-FDTD method combined with the proposed PI model and (a) /12    , (b) /16    , and (c) / 20   .Note that  is the grid size, k the physical wavenumber,  the wavelength, and ,

Figure 8 .
Figure 8. Propagation error of the numerical wavenumber, num k , for the NS-FDTD method combined with the proposed PI model.Note that  is the grid size, k the physical wavenumber,  the wavelength, and ,   the angles of the spherical coordinate system   , , r   .The PI model at- tains a high propagation accuracy, practically equivalent to that of the NS-FDTD method along the axial direction 0    .

Figure 9 .
Figure 9. Propagation error of the numerical wavenumber, num k , for the FDTD method combined with the proposed PI model.Note that ,   are the angles of the spherical coordinate system Figure 10b-d, and 3

Figure 10 .
Figure10.The proposed PI model for a thin wire antenna of radius a .Integral path (a) in the vi- cinity of the wire at a xy cross section, (b) for a wire end on the z axis, (c) in the vicinity of the wire along the z axis (this part uses only the basic path), and (d) for the feed point with a  gap on the z axis.The black solid line signifies the basic path and dashed blue line the modified complementary path, while, for simplicity, y axis, in (b-d), and the z axis, in (a), are omit- ted.Designating as SB the area enclosed by the basic path and SC the area enclosed by the complementary path, the integral area for t S d   H S is: SB = 2 2 / 4 a   

Figure 11 .
Figure 11.Antenna impedance, ant Z R jX   , as in Equation (44), obtained via the proposed PI model, the FDTD method, and NEC2 [60] (reference solution), for 1   m and (a)

Table 1 .
Strategy status of the curved surface analysis in the NS-FDTD method (Efficacy: Decreasing rate of the calculation cell number used for the verification of the model in the paper or digest).
term, the basic integral path is utilized.Note that the feed point with a  gap at k s  .Finally, for the t  E