Comparison between Specialized Quadrature Rules for Method of Moments with NURBS Modelling Applied to Periodic Multilayer Structures

Abstract

A comparison between Ma-Rokhlin-Wandzura (MRW) and double exponential (DE) quadrature rules for numerical integration of method of moments (MoM) matrix entries with singular behavior is presented for multilayer periodic structures. Non Uniform Rational B-Splines (NURBS) modelling of the layout surfaces is implemented to provide high-order description of the geometry. The comparison is carried out in order to show that quadrature rule is more suitable for MoM matrix computation in terms of sampling, accuracy of computation of MoM matrix, and CPU time consumption. The comparison of CPU time consumption shows that the numerical integration with MRW samples is roughly 15 times faster than that numerical integration using DE samples for results with similar accuracies. These promising results encourage to carry out a comparison with results obtained in previous works where a specialized approach for the specific analysis of split rings geometries was carried out. This previous approach uses spectral MoM version with specific entire domain basis function with edge singularities defined on split ring geometry. Thus, the previous approach provides accurate results with low CPU time consumption to be compared. The comparison shows that CPU time consumption obtained by MRW samples is similar to the CPU time consumption required by the previous work of specific analysis of split rings geometries. The fact that similar CPU time consumptions are obtained by MRW quadrature rules for modelling of general planar geometries and by the specialized approach for split ring geometry provides an assessment for the usage of the MRW quadrature rules and NURBS modelling. This fact provides an efficient tool for analysis of reflectarray elements with general planar layout geometries, which is suitable for reflectarray designs under local periodicity assumption where a huge number of periodic multilayer structures have to be analyzed.

1. Introduction

Many electrical devices built using planar multilayer structures such as frequency selective surfaces (FSS) [1], reflectarray/transmitarray [2,3] antennas, phased array antennas, and metasurface (MTS) leakywave antennas [4,5,6] are analyzed many times assuming a periodic environment. This assumption is known as local periodicity assumption and it was experimentally validated in the literature, and was enabled for accurate design of many antennas [3]. In fact, when the local periodicity assumption is used in the design of a reflectarray, the scattering problem of a plane wave obliquely incident on a periodic multilayered structure has to be analyzed a huge number of times. Therefore, very efficient numerical tools are required for the solution of this scattering problem.

Since most of these devices work with resonant elements, an accurate modelling of the geometry of the resonant elements is required. NURBS surfaces show high-order description of the geometry of layout for complex geometries [7]. Mature Computer Aided Geometric Design (CAGD) based on NURBS meshing is available. In this work, the NURBS meshing technique used is based on a redesigned and optimized Paving Algorithm [8]. Although a mesh of NURBS surfaces of an arbitrary layout is accurate for constructing and storing the surface of layout, more suitable surface meshing for numerical computation of parameters associated with the surface is preferred (curvature, derivatives, integrations, etc). Bézier patches provide these desired properties with parametric surfaces defined by Bernstein polynomials [9]. Moreover, NURBS surfaces can be efficiently treated using Bézier patches [9,10] by the Cox-de Boor transformation algorithm [11]. In this context, generalized subsectional rooftops functions are usually defined between pairs of adjacent Bézier patches to approximate the surface current densities induced on the layout [7].

The most used method for determination of the induced current densities in periodic multilayered structures is the Method of Moments (MoM) in the spectral domain [12,13]. Unfortunately, when the MoM in the spectral domain is applied, the MoM matrix entries are a double infinite series with slow convergence, and larger CPU time consumption is required for an accurate computation of these double series. In order to speed up this convergence, several approaches have been proposed. Two-dimensional fast Fourier transforms were implemented by Chan and Mittra to compute the series when subsectional rectangular or triangular basis functions (BF) are used to model the current density induced on the layout [14,15]. However, the edges of the layout have to exactly coincide with the cells of a uniform mesh (rectangular or triangular) defined in the unit cell of the periodic structure. In [16], the computation is accelerated by means of the combined use of Kummer’s transformation, Poisson’s formula, and Chebyshev polynomial interpolation. However, a uniform mesh is required in all these proposals and, unlike NURBS surfaces, these types of meshes provide bad modelling of complex resonant geometries, particularly curved geometries. In a recent work, a sophisticated and efficient treatment for MoM in spectral domain was developed for specific entire domain basis functions which accounts for the singularities of the surface current densities at edges of split ring geometries [17]. In that treatment, the Fourier transform of basis functions, which appear in the double infinite series, are expressed in terms of Hankel transforms for this specific geometry of split rings. These Hankel transforms can be efficiently interpolated in terms of Chebyshev polynomials. Moreover, asymptotic term of these Fourier transforms is easily identifiable. This fact provides a controlled truncation of the infinite series in terms of desired accuracy. However, this development is only suitable for the analysis of the split-rings geometry.

An interesting alternative to the MoM in the spectral domain for the analysis of multilayered periodic structures is the Mixed Potential Integral Equation (MPIE) formulation of the MoM in the spatial domain [18,19,20]. In this latter approach, one has to face the computation of multilayered periodic Green’s functions where a double infinite series with slow convergence has to be computed. Moreover, in this approach, integrals with singular integrands appear in the computation of MoM matrix entries. Fortunately, there are efficient and accurate techniques for fast computation of periodic Green’s functions for multilayer medium. In [21], the periodic Green’s functions of the periodic multilayer structure are judiciously interpolated in the spatial domain in terms of 2-D Chebyshev polynomials after extracting the singular behaviour of the Green’s functions around the source points (which includes the source singularities plus the images through the closest layers). The evaluation of singular integrals is a problem that has been addressed in detail in [22,23,24] for subsectional triangular BF and free space Green’s functions. Moreover, sophisticated and specific techniques have been proposed for periodic multilayer structures using entire domain basis functions, which account for the singularities of the surface current densities at edges of patches or coplanar dipoles [21,25]. However, these approaches are not suitable for NURBS surfaces in MoM context.

When the layout geometries are modelled by NURBS, these singular integrals have to be computed by quadrature rules. The singular integrals can be accurately computed by means of specialized quadrature rules for integration of singular behavior located in the limit of the integration domain. Ma-Rokhlin-Wandzura (MRW) quadrature rule [26] and Double Exponential (DE) quadrature rule [27] are specialized quadrature rules for singular integrals. The MRW quadrature rule is designed for integrals of functions with logarithmic singularities in the lower limit of the integration interval, while the DE quadrature rule is suitable for integrals of functions with a wide variety of singular behaviors located at endpoints of the integration interval. Both MRW [21,25,28] and DE [27,29,30] quadrature rules have been successfully used for the determination of integrals with singularities that arise in the application of the MoM for multilayer problems. In more recent works results obtained by MRW [31,32,33] or by DE [34,35], quadrature rules are used as comparison patterns. However, the performance of these two types of quadrature rules have never been compared in an MoM context. In fact, the MRW quadrature has only been used for integration of the logarithmic singular behavior of integrands that arise in MoM. In this work, we show a comparison between both quadrature rules for integration of MoM matrix entries with singular behavior for NURBS modelling of general planar geometry of layout embedded in multilayer periodic structures. The comparison is carried out in terms of sampling of quadrature rules, accuracy, and CPU time consumption for computation of MoM matrix entries.

This paper is organized as follows. Section 2 shows a description of the solution of the MPIE using generalized rooftop defined on a pair of adjacent Bézier patches. In this section, the formulation of numerical integration involved in the MoM matrix entries is shown by MRW and DE quadrature rules. Section 3 shows convergent patterns of results for computation of MoM matrix entries using both types of quadrature rules. In this section, validations and comparisons of CPU time consumptions are also carried out for analysis of a multilayer periodic structure based on dual concentric split-rings geometry. Finally, the conclusions are provided in Section 4.

2. Description of the Problem

Let us consider a periodic multilayer structure of p_x × p_y cell size with NC dielectric layers of thickness h_p (p = 1,…, NC) and complex permittivity ε_p= ε₀ε_r,p(1-jtanδ_p). The multilayer medium is upper limited by free space and lower limited by either free space or ground plane. This multilayer medium hosts an interface with periodic planar layout with negligible thickness and arbitrary geometry (see unit cell in Figure 1). For simplicity in this work, an unique interface with periodic layout is considered but more interfaces with periodic layout can be considered. These surfaces of planar layouts will be considered as PEC throughout. On these surfaces, surface current densities J(x,y) and surface charge densities σ(x,y) are induced by either linearly or circular polarized incident plane wave with incidence direction given by the angular spherical coordinate θ_inc and φ_inc. Both, current and charge densities are related to each other by a known continuity equation. Since the periodic surfaces are planar, the induced current densities J(x,y) have no z-component. These induced currents and charge densities satisfy a MPIE from boundary conditions on PEC [19,20]. In order to solve the MPIE, expansion in term of known BFs weighted by unknown coefficients, and method of weighted residual with weighting functions (WFs) are applied.

In this work, the planar surface of the layout hosted of the unit cell of the periodic structure is modelled by NURBS surfaces. The NURBS are efficiently discretized as piecewise Bezier patches [7,9,10] using the Cox-de Boor transformation algorithm [11]. Figure 1 shows in grey color the planar surface with generic geometry, which is discretized in terms of Bézier patches. Each Bézier patch is a parametric surface defined by rectangular parametric coordinates u and v. Both rectangular parametric coordinates, u and v, vary from 0 to 1 along the Bézier patch. Thus, known relationships between Cartesian coordinates (x, y) and parametric coordinates (u, v) are defined in each Bézier patch by a pre-processing procedure [7]:

$$\begin{array}{l}x(u,v)={\displaystyle \sum _{p=0}^2{\displaystyle \sum _{q=0}^2{x}_{pq}{B}_p^2(u)B_q^2(v)}}\\ y(u,v)={\displaystyle \sum _{p=0}^2{\displaystyle \sum _{q=0}^2{y}_{pq}{B}_p^2(u)B_q^2(v)}}\end{array}\}$$

(1)

where B_p²(·) is the pth-Bernstein polynomial [9] of degree 2 and values x_pq and y_pq (p,q = 0, 1, 2) are the control points of Bézier patches, which are obtained in the pre-processing procedure. Once these relationships given in (1) are available for each Bézier patch, subsectional rooftops J_j(u,v) (j = 1,…, N_b) functions and razor blade functions can be defined between a pair of adjacent Bézier patches from their parametric coordinates u and v, as proposed in [7]. The razor-blade joins the centers of each adjacent Bezier patches (i.e., values of x and y given by (1) when u = v = 0.5 in each adjacent Bezier patch) by means the isoparametric curved lines C_i (i = 1,…, N_b).

Therefore, the rooftops functions are used as BFs (j = 1,…, N_b) and razor-blade are used as WFs [7]. Taking into account the relationships given in (1), the expansion of the surface current densities J(x,y) in term of the known BFs J_j(u,v) (j = 1,…, N_b) are given by:

$$J[x(u,v),y(u,v)]={\displaystyle \sum _{j=1}^{{\mathrm{N}}_{\mathrm{b}}}{c}_{j,}{J}_j(u,v)}$$

(2)

where c_j (j = 1,…, N_b) are unknown coefficients.

When the expansion given in (2) is substituted in the MPIE and razor-blade are used as WFs, the resultant system of equations for the unknown coefficients c_j (j = 1,…, N_b) is obtained:

$${\displaystyle \sum _{j=1}^{{\mathrm{N}}_{\mathrm{b}}}\left(Z_{ij}^{\mathrm{ind}}+Z_{ij}^{\mathrm{cap}}\right)c_j=e_i}(i=1,\dots ,{\mathrm{N}}_{\mathrm{b}})$$

(3)

where the coefficients e_i of the system of the linear equations can be computed by the next line integral:

$$e_i={\displaystyle \underset{C_i}{\int }{E}_t^{\mathrm{exc}}(x,y,z=z_k)\cdot dr}$$

(4)

With E_t^exc(x,y,z = z_k), the tangential electric field generated in the observation point (x,y,z = z_k) by multilayer substrate in the absence of the patches when a plane wave impinges on the multilayer medium. The coefficients Z_ij^ind and Z_ij^cap are obtained as follows:

$$Z_{ij}^{\mathrm{ind}}=\mathrm{j}\omega {\displaystyle \underset{C_i}{\int }{f}_j^{\mathrm{ind}}(x,y,z=z_k)\cdot dr}$$

(5)

$$Z_{ij}^{\mathrm{cap}}=[f_j^{\mathrm{cap}}(x_i^{+},y_i^{+},z=z_k)-f_j^{\mathrm{cap}}(x_i^{-},y_i^{-},z=z_k)]$$

(6)

where:

$$f_j^{\mathrm{ind}}(x,y,z=z_k)={\displaystyle \underset{0}{\overset{1}{\int }}{\displaystyle \underset{0}{\overset{1}{\int }}{G}_{xx}^{\mathrm{A}}[x-x^{\prime }(u,v),y-y^{\prime }(u,v),z=z^{\prime }=z_k]J_j(u,v)|J_{uv}(u,v)|dudv}}$$

(7)

$$f_j^{\mathrm{cap}}(x,y,z=z_k)={\displaystyle \underset{0}{\overset{1}{\int }}{\displaystyle \underset{0}{\overset{1}{\int }}{G}^{\Phi }[x-x^{\prime }(u,v),y-y^{\prime }(u,v),z=z^{\prime }=z_k]\nabla J_j(u,v)|J_{uv}(u,v)|dudv}}$$

(8)

With |J_uv(u,v)| the determinant of Jacobian matrix of the relationships given in (1). G^A_xx is the periodic Green’s function for the x-component of the vector potential and G^Φ is the periodic Green’s function for the scalar potential, respectively, both for the periodic multilayer structure of Figure 1. In this work, the periodic Green’s functions of the periodic multilayer structure are efficiently obtained by the pre-processing procedure of interpolation described in [21]. In this procedure, the periodic Green’s functions of the periodic multilayer structure are judiciously interpolated in the spatial domain in terms of 2-D Chebyshev polynomials after extracting the singular behaviour of the Green’s functions around the source points (which includes the source singularities plus the images through the closest layers). The singular behaviour of the periodic Green’s functions is proportional to 1/(4πρ) where ρ = [(x − x^’)² + (y − y^’)²]^1/2 [21]. Since the periodic Green’s functions provide singular behaviour when the Cartesian coordinates x and y of the observation point and the Cartesian coordinates x^’ and y^’ of the source point are close, the integrands of the bi-dimensional integrals shown in (7) and (8) show a similar singular behaviour. This singular behaviour of these integrands is produced in parametric domain of Bézier patch around values u_obs and v_obs, which satisfy the following system of equation:

$$\begin{array}{l}x=x^{\prime }(u_{\mathrm{obs}},v_{\mathrm{obs}})\\ y=y^{\prime }(u_{\mathrm{obs}},v_{\mathrm{obs}})\end{array}\}$$

(9)

where the relationships between the Cartesian variables x^’ and y^’ and the parametric variables u and v are given by (1) in each Bézier patch. Once the values u_obs and v_obs are available for each observation point in each Bézier patch, the integrals (7) and (8) can be handled in order to take into account the singular behavior of the integrands around the values u_obs and v_obs.

Let us consider G[x-x^’(u,v),y-y^’(u,v)] stands either Green’s function G^A_xx or G^Φ for fixed values of Cartesian coordinates x and y of the observation point in the kth-interface (i.e., when z = z^’ = z_k). Let us also consider that K(u,v) stands either x or y component of a fixed BF or the divergence of fixed BF. Then the integrals (7) and (8) can be expressed as:

$$\begin{array}{l}{\displaystyle \underset{0}{\overset{1}{\int }}{\displaystyle \underset{0}{\overset{1}{\int }}F[u,v]dudv}}=\\ \underset{\mathrm{I}}{\underbrace{{\displaystyle \underset{0}{\overset{v_{\mathrm{obs}}}{\int }}{\displaystyle \underset{0}{\overset{u_{\mathrm{obs}}}{\int }}F[u,v]dudv}}}}+\underset{\mathrm{II}}{\underbrace{{\displaystyle \underset{v_{\mathrm{obs}}}{\overset{1}{\int }}{\displaystyle \underset{0}{\overset{u_{\mathrm{obs}}}{\int }}F[u,v]dudv}}}}+\underset{\mathrm{III}}{\underbrace{{\displaystyle \underset{v_{\mathrm{obs}}}{\overset{1}{\int }}{\displaystyle \underset{u_{\mathrm{obs}}}{\overset{1}{\int }}F[u,v]dudv}}}}+\underset{\mathrm{IV}}{\underbrace{{\displaystyle \underset{0}{\overset{v_{\mathrm{obs}}}{\int }}{\displaystyle \underset{u_{\mathrm{obs}}}{\overset{1}{\int }}F[u,v]dudv}}}}\end{array}$$

(10)

where,

$$F[u,v]=G[x-x^{\prime }(u,v),y-y^{\prime }(u,v)]K(u,v)|J_{uv}(u,v)|$$

(11)

Since the Green’s functions provide singular behaviour (i.e., G[x − x^’(u,v),y − y^’(u,v)] for fixed values of x and y is proportional to 1/(4πρ)) when (9) is satisfied, the square domain 0 < u,v < 1 of the bi-dimensional integral of F(u,v) is divided in four rectangular regions where the values of the values u_obs and v_obs are located in the corner of each region, as shown in Figure 2:

Each integral of the right side of (10) can be accurately computed by means of specialized quadrature rules for integration of singular behavior located in the limit of the integration domain using an appropriate trivial variable change. In this work, we compare the performance of two specialized quadrature rules for the numerical integration of (7) and (8) by the division of regions given in (10): MRW quadrature rules [26] and DE quadrature rules [27]. The MRW quadrature rules are specifically designed for integrals of functions with logarithmic singularities in the lower limit of the integration interval (0,1), while the DE quadrature rules are suitable for integrals of functions with a wide variety of singular behaviors located at endpoints of the integration interval (−1,1). Let us consider t_m, w_m abscissas and weights, respectively, for either MRW or DE quadrature rules. In order to show an example, let us consider the term II of the right side of (10). If we define the variable changes u = (1 − t_u)u_obs, v = (1 − v_obs)t_v + v_obs, then integration domain of the new variables t_u and t_v is (0,1) × (0,1). The singular behavior is produced in the integration domain around of the values t_u = 0 and t_v = 0 (i.e., singular behavior is produced in the lower limit of integration interval of new variables t_u and t_v). In this condition, the integral can be accurately computed from t_m, w_m abscissas and weights obtained by MRW quadrature rules, as shown below:

$${\displaystyle \underset{v_{\mathrm{obs}}}{\overset{1}{\int }}{\displaystyle \underset{0}{\overset{u_{\mathrm{obs}}}{\int }}F[u,v]dudv}}=u_{\mathrm{obs}}(1-v_{obs}){\displaystyle \sum _{m=1}^{M_{\mathrm{MRW}}}{\displaystyle \sum _{n=1}^{N_{\mathrm{MRW}}}{w}_m{w}_{n}F[u(t_m),v(t_n)]}}$$

(12)

where M_MRW × N_MRW are the number of samples used in the MRW quadrature rule. Similar treatment can be made for the remaining integrals to the right side of (10).

On the other hand, if we define the variable changes u = (t_u + 1)u_obs/2, v = [(1 − v_obs)t_v + (1 + v_obs)]/2, then integration domain of the new variables t_u and t_v is (−1,1) × (−1,1). In this case, singular behavior is produced in the integration domain around of the values t_u = 1 and t_v = −1 (i.e., singular behavior is produced at endpoints of the integration interval of new variables t_u and t_v). In this condition, the integral can be accurately computed from t_m, w_m abscissas and weights obtained by DE quadrature rules, as shown below:

$${\displaystyle \underset{v_{\mathrm{obs}}}{\overset{1}{\int }}{\displaystyle \underset{0}{\overset{u_{\mathrm{obs}}}{\int }}F[u,v]dudv}}=\frac{u_{\mathrm{obs}}(1-v_{obs})}{4}{\displaystyle \sum _{m=1}^{M_{\mathrm{DE}}}{\displaystyle \sum _{n=1}^{N_{\mathrm{DE}}}{w}_m{w}_{n}F[u(t_m),v(t_n)]}}$$

(13)

where M_DE × N_DE are the number of samples used in the DE quadrature rule. Similar treatment can be done for the remaining integrals to the right side of (10). Both MRW [21,25,28] and DE [27,29,30] quadrature rules have been successfully used for the determination of integrals with singularities that arise in the application of the MoM for multilayer problems. In more recent works, results obtained by MRW [31,32,33] or by DE [34,35] quadrature rules are used as comparison patterns. However, the performance of these two types of quadrature rules have never been compared to each other in an MoM context. In fact, the MRW quadrature has been used in the literature only for integration of logarithmic singular behavior of integrands that arise in the MoM context and not for integration of singular behavior proportional to 1/(4πρ). In this work, we show a comparison of both types of quadrature rules in terms of sampling, accuracy, and CPU time consumption for integration of (7) and (8) with singular behavior proportional to 1/(4πρ) when the observation and source point are close. These comparisons are obtained by sequential programming of home-made FORTRAN code (parallel programming is not implemented).

3. Numerical Results

Let us consider the reflectarray element based on dual concentric split-rings, as shown in Figure 3a:

Figure 3a shows a periodic cell of p_x × p_y = 5 × 5 mm² size, which consists of dual concentric split rings centered in the center. The split-rings of fixed w width are printed on 0.787 mm thick Rogers Duroid 5880 (ε_r,1 = 2.2 and tanδ₁ = 0.0009) [36]. The outer split-ring consists of symmetrical arcs with subtended angle ψ₂ and inner radious ρ₂. The inner split-ring also shows symmetrical arcs but, with subtended angle ψ₁ and inner radious ρ₁. Fixed simmetrical axis are considered for outer split-ring, while the simmetrical axis of the inner split-ring are able to rotate an α angle, as shown in Figure 3a. This reflectarray element has provided satisfactory electrical reflexion properties under local perodicity assumption (i.e., the reflectarray element is in periodic enviroment) in circular polarized reflectarray designs [17,37] by means of a variable rotation technique [38,39]. The outer split-rings work at 19.95 GHz, while the inner split-rings work at 29.75 GHz. Let us consider the following values of the geometrical parameters ψ₁ = 162.8 deg, ψ₂ = 150.4 deg, α = 0 deg, ρ₂ = 1.85 mm, ρ₁ = 1.20 mm, w = 0.2 mm for the following results that will be shown at 29.75 GHz under normal incident (i.e., θ_inc = φ_inc = 0). These values of geometrical parameters are not arbytrary since these values ensure that the phase of the complex number of copolar reflexion coefficients for ‘x’ and ‘y’ direction are shifted 180 deg to each other at 29.75 GHz under normal incidence [17]. This condition of 180 deg in the shifting of the phases is a necessary condition to apply VRT [38,39] in the reflectarray design. Figure 3b shows the discretization of the dual concentric split-rings in terms of Bézier patches. Since the inner split-ring works at 29.75 GHz and the outer split-rings work at lower frequency of 11.95 GHz, the inner split-ring is discretized in terms of more Bézier patches than the outer split rings to provide accure results at 29.75 GHz.

Figure 4a,b show the magnitude and phase of the complex number of the integrand of (8) in the integration domain given for region I shown in Figure 2 for a fixed Bézier patch under study. Similar results are obtained for the integrand of (7) (not shown). We can see that there is singular behaviour of the magnitude of the integrand around the values u_obs = v_obs = 0.485 (i.e., the values of the coordinates ‘x’ and ‘y’ of the observation point have been selected in order to satisfy the condition (9) when u_obs = v_obs = 0.485). The distribution of the 7 × 7 DE samples and 5 × 5 MRW samples are also shown in the integration domain. We would like to point out that the 7 × 7 DE samples are the minimum number of samples provided by DE quadrature rules when the level of the quadrature is null [27]. In Figure 4a, we can see that the singular behaviour is captured pretty well by the 7 × 7 DE samples. However, this fact is obtained by the cost of oversampling the rest of corners of the rectangular integration domain. Thus, the rest of the integration domain is badly sampled by DE quadrature rules. This is because the DE quadrature rules are suitable for integrals of functions with singular behaviors located at endpoints of the integration domain. Note that, due to this, the rest of the integration domain is badly sampled by DE quadrature. Moreover, this fact could provide a bad capture of the soft behavior of the phase of the integrand in the integration domain, as shown in Figure 4b. On the contrary, the MRW sampling is specifically designed for integrals of functions with singularities in the lower limit of the integration interval. In this way, the 5 × 5 MRW samples keep an optimal sampling between the sampling of the singular behavior of the magnitude of the integrand and the sampling of the magnitude and phase of integrand in the rest of the integration domain.

Table 1. Convergence pattern with respect to the number of MRW samples per region of the inductive, Z₁₁^ind, and capacitive, Z₁₁^cap, contributions of the MoM coefficient matrix.

Number of Samples MMRW × NMRW per Region	Z11ind	Z11cap
3 × 3	(0.2167181 + 14.19454j)	(−0.0061680079 − 712.6991j)
5 × 5	(0.2167180 + 14.29379j)	(−0.0056559443 − 699.2062j)
10 × 10	(0.2167187 + 14.32665j)	(−0.0058996081 − 695.4257j)
15 × 15	(0.2167188 + 14.32226j)	(−0.0057702065 − 695.0881j)
20 × 20	(0.2167189 + 14.31678j)	(−0.0057617426 − 695.0215j)
25 × 25	(0.2167189 + 14.31555j)	(−0.0057619214 − 695.0031j)

shows the convergence pattern of the inductive, Z₁₁^ind, and capacitive, Z₁₁^cap, contributions of the MoM coefficient matrix of the system of linear equations given in (3). This convergence pattern is shown with respect to the number of MRW samples per region used in the integration of (7) and (8). Since the functions f_j^ind (j = 1,…,N_b) have no singular behavior, the line integral given in (5) can be numerically computed by conventional Gauss-Legendre (GL) quadrature rules [40]. In this work, we use 4 GL samples for integration of the line integral shown in (5). A similar convergence pattern is shown in

Table 2. Convergence pattern with respect to the number of DE samples per region of the inductive, Z₁₁^ind, and capacitive, Z₁₁^cap, contributions of the MoM coefficient matrix.

Number of Samples MDE × NDE per Region	Z11ind	Z11cap
7 × 7 (level 0 of DE quadrature)	(0.2233992 + 14.89468j)	(−0.0060498714 − 718.6270j)
13 × 13 (level 1 of DE quadrature)	(0.2167199 + 14.45327j)	(−0.0056549907 − 697.3384j)
25 × 25 (level 2 of DE quadrature)	(0.2167186 + 14.35121j)	(−0.0057088733 − 695.0247j)
51 × 51 (level 3 of DE quadrature)	(0.2167187 + 14.32173j)	(−0.0057591200 − 694.9873j)
101 × 101 (level 4 of DE quadrature)	(0.2167180 + 14.31637j)	(−0.0057511926 − 694.9833j)
203 × 203 (level 5 of DE quadrature)	(0.2167156 + 14.31532j)	(−0.0057668090 − 694.9734j)

but, in this case with respect to the number of DE samples used in the integration of (7) and (8).

If we compare the results shown in Table 1 and Table 2, we can see that a lesser number of samples are required by MRW quadrature than that required by DE quadrature to reach the same number of significant digits. For example, in order to reach three significant digits of imaginary part of Z₁₁^ind and Z₁₁^cap (note that the imaginary part is significantly higher than the real part), MRW quadrature requires 10 × 10 samples per each region shown in Figure 2, while the DE quadrature requires 51 × 51 samples per region. Note that this fact is produced in spite of MRW quadrature rules specifically being designed for integrals of functions with logarithmic singularities and not with singularities proportional to 1/(4πρ). These results anticipate a CPU time savings when MRW quadrature rules are used instead of DE quadrature rules.

In order to compare the CPU time consumption by MRW and DE quadrature rules, the phase curves shown at 29.75 GHz in [17] have been reproduced by MRW and DE quadrature rules in the numerical integration of (7) and (8). The geometrical variables ψ₁ and α used to reproduce results given in [17] are shown in Figure 5 (i.e., these values of the geometrical parameters ψ₁ and α ensure that the required condition of 180 deg of phase difference is satisfied at 29.75 GHz to apply VRT).

In our analysis, 284 BFs have been used (i.e., N_b = 284 has been used). Figure 6 shows the reproduced results using the 3 × 3 MRW samples per region, 7 × 7 DE samples per region and 13 × 13 DE samples per region. The results obtained in [17] by the proposed method in [17] and by CST simulations carried out in [17] are also shown. We would like to point out that the proposed method in [17] consists of a sophisticated and efficient treatment in the spectral domain for specific entire domain basis functions, which account for the singularities of the surface current densities at the edges of the arcs of the split rings. In that treatment, the MoM matrix entries are expressed as infinites series where Fourier transforms of the entire domain basis functions with edge singularities are involved. In order to compute these Fourier transforms, the part of the basis functions, which depends on an angular polar coordinate, is expanded as Fourier series. This expansion leads Fourier transforms of the basis function in terms of Hankel transforms. These Hankel transforms can be efficiently interpolated in terms of Chebyshev polynomials. Moreover, due to Fourier expansion, the asymptotic term of these Fourier transforms is easily identifiable. This fact provides a controlled truncation of the infinite series in terms of desired accuracy. Thus, important CPU time savings is expected when the proposed method in [17] is applied since it is specifically designed and speeded up for analysis of the split-rings geometry, while the implementation of specialized quadrature rules for integration of (7) and (8) consider the general planar surface of complex geometries.

In Figure 6, we can see that the results obtained by 3 × 3 MRW samples and 13 × 13 DE samples per region have good agreements with the results provided by the proposed method in [17] and CST simulations shown in [17]. However, the results provided by the 7 × 7 DE samples per region show distorted phase curves (i.e., 7 × 7 DE samples per region are insufficient to achieve convergence).

Table 3. CPU time consumption required per point of the phase curves of ∠R_RHCP,RHCP and ∠R_LHCP,LHCP for results shown in Figure 6.

Type of Analysis Method	CPU Time (s) per Point	Relative Percent
3 × 3 MRW samples per region	6.42	6.60%
13 × 13 DE samples per region	97.28	100%
Proposed method in [17]	4.93	5.07%
CST simulations in [17]	80.19	82.4%

shows the CPU time consumption required per point of each pair of phase curves (∠R_RHCP,RHCP and ∠R_LHCP,LHCP) shown in Figure 6 by each type of analysis method. Relative percent of CPU time consumption is also shown with respect to the maximum CPU time required in order to ease comparisons. The CPU time consumption required by 7 × 7 DE samples per region is not shown since these samples are insufficient to achieve convergence. These results were obtained on a laptop computer with processor Intel Core i7-6700HQ, 2.6 GHz of clock frequency with 32 GB of RAM. We would like to point out that the laptop computer used for our results is the same laptop computer used in [17]. In this way, an honest comparison of CPU time consumption can be carried out with that required by the shown simulations in [17].

In Table 3, we can see that the analysis carried out by the numerical integration of (7) and (8) using 3 × 3 MRW samples per region is roughly 15 times faster than that numerical integration using 13 × 13 DE samples per region (i.e., the relative percent of CPU time consumption obtained with 3 × 3 MRW samples with respect to that obtained with 13 × 13 DE samples is roughly 6.6%). Moreover, we can see that the CPU time consumption obtained by 3 × 3 MRW samples per region is comparable with the CPU time consumption required by the proposed method in [17] (i.e., the relative percent of CPU time consumption obtained with proposed method in [17] with respect to that obtained with 13 × 13 DE samples is roughly 5.07%). On the other hand, the CPU time consumption obtained by 13 × 13 DE samples per region is comparable with the CPU time consumption required by CST simulations in [17] (i.e., the relative percent of CPU time consumption obtained with CST simulations in [17] with respect to the obtained 13 × 13 DE samples is roughly 82.4%). We would like to remember that the implementation of MRW quadrature rules for integration of (7) and (8) consider a general planar surface of complex geometries, while the proposed method in [17] is designed for specific analysis of the split-rings geometry. In this way, the fact that similar CPU time consumptions are obtained by both general purpose software with MRW quadrature rules and proposed method in [17] provides favourable results for the general purpose software. Thus, this fact justifies the effort to implement MRW quadrature in MoM context with NURBS modelling of geometries.

4. Conclusions

A comparison of the accuracy in terms of samples used by MRW and DE quadrature rules for integration of MoM matrix entries with singular behavior have been carried out for multilayer periodic structures which host a general planar layout. A decomposition of the integration domain in four regions are carried out for efficient implementation of both types of quadrature rules. In order to provide high-order description of the geometry of layout for complex geometries, NURBS modelling of layout surface has been implemented and generalized rooftops defined on pair of adjacent Bézier patches have been used in the approximation of induced surface current densities on the layout.

Samplings of MRW and DE quadrature have been compared to capture the behavior of the integrands of the integrals that arise in the computation of MoM matrix entries. The singular behaviour of the integrand is captured pretty well by the DE samples, but at the cost to oversample the rest of corners of the rectangular integration domain. Thus, the rest of the integration domain is badly sampled by DE quadrature rules. On the contrary, MRW samples keep an optimal sampling between the sampling of the singular behavior and the sampling in the rest of the integration domain. Convergence patterns of MoM matrix entries with respect to the number of samples are compared for both types of quadrature rules. A smaller number of samples are required by MRW quadrature than those required by DE quadrature to reach the same number of significant digits.

Finally, the obtained promising results encourage a comparison with results obtained in previous works, where a specialized approach for the specific analysis of split rings geometries was carried out. The results shown in this previous work for the analysis of split rings geometries have been reproduced by MRW and DE quadrature rules and NURBS modelling of layout geometry. Good agreements with results provided by previous works are obtained when 3 × 3 MRW samples and 13 × 13 DE samples per region are used. Comparison of CPU time consumption shows that the numerical integration with 3 × 3 MRW samples per region is roughly 15 times faster than that numerical integration using 13 × 13 DE samples per region. Moreover, the CPU time consumption obtained by 3 × 3 MRW samples per region is similar to the CPU time consumption required by previous work designed for specific analysis of the split-rings geometry.

Therefore, this fact provides an efficient tool for the analysis of reflectarray elements with general planar layout geometries, which is suitable for reflectarray designs under local periodicity assumption where a huge number of periodic multilayer structures have to be analyzed.