Optimal High-Order Method of Moments combined with NURBS for the scattering by a 2 D cylinder

This paper deals with the High-Order Method of Moments (HO-MoM) combined with Non-Uniform Rational Basis Splines (NURBS) segments to evaluate the scattering by a 2D cylinder. The authors mainly focus upon the influence of the different parameters (polynomial basis, order, mesh length, curvature, polarization,...) and try to determine if a optimal choice exists for the convergence speed.


Introduction
The Boundary Element Method (BEM) usually called the Method of Moments (MoM) has definitely become a well established numerical approach to model the electromagnetic waves scattering.The reduction of the problem dimension by one clearly appears as a major reason for the success of MoM [1].Basically, this method induces a discretization of field integral equations in the vicinity of boundaries, and usually approximates the surface currents using constant or linear subsectional basis functions (pulse or triangle functions for 2D and RWG or rooftop functions for 3D).It is generally admitted that MoM provides a reliable estimation where the size of the elements in the mesh is lower than ten percent of the wavelength (λ/10).
For the sake of optimal performance, High-Order MoM (HO-MoM), involving high order polynomial functions, has been intensively explored for more than a decade [2,3,4,5].In many configurations, the relative acceleration with respect to the basic MoM has been evaluated [6].It is noteworthy that, in those publications, the high order models only concern the electromagnetic components (currents and scattered fields) of the scattering issue.The geometry of the scatterer always remains a mesh of piecewise linear functions.More precisely, the standard meshing technique decomposes the scatterer's boundaries to sub-domains of straight lines (2D problem) or planar triangles/quadrilaterals (3D problem).For plane and low curvature objects, this approach does not raise particular problems.
For a high curvature scatterer, the criterion λ/10 for the mesh elementary size often imposed by MoM is far from being enough to obtain reliable estimation of the scattered fields.In this situation, MoM and standard HO-MoM require very small sized elements and the number of unknowns can severely increase.A solution consists in modeling the geometry of the scatterer's boundaries with Non-Uniform Rational Basis Splines (NURBS) elements.NURBS elements approximate the geometry using a high order polynomial basis and easily model high curvature objects.Then, the combination of the HO-MOM with NURBS meshing constitutes a really powerful approach whose advantages are quite recently studied [7,8].
Nevertheless, NURBS meshing combined with HO-MoM constitutes a quite complex method.The number of unknowns depends on the number of the elements in the mesh (C) and the order of polynomial basis function (N ).More, this approach can be applied with continuous or discontinuous polynomial functions.Finally, the determination of the best compromise that optimize the accuracy of this numerical estimation for a given scatterer obviously remains a difficult task.As a matter of fact, this optimal compromise notably depends upon the characteristics of the scatterers, the electromagnetic characteristics (polarization for instance) or the position of the observation (far field versus near field).For the sake of clarity, we consider that it is important to evaluate the influence of the different parameters for the HO-MOM approach combined with NURBS in the case of canonical scatterers.
With this purpose in mind, this paper deals with the HO-MoM+NURBS modeling of the scattering by a very basic 2-dimensional scatterer.In fact, HO-MoM+NURBS is used to simulate the scattering in TM and TE polarization for a PEC circle.The accuracy of the method is investigated by comparing the surface currents and the Radar Cross Sections (RCS) with the analytical solution obtained by modal decomposition (2D Mie theory).The Radar Cross Section appears as a reliable description of the far field scattering.On the opposite, the surface current is far more related to the very near field scattering.
In this paper, we specifically focus on the convergence to the numerical simulation at different levels, for different polarizations and for two components (surface currents and RCS).To a lesser extent, the computational time as a function of segment number (C) and the order of basis function (N ) is also presented.
In the following section, we point out the notion of curvature and the relation with NURBS curve expression adapted for our computational electromagnetic method.In section 3, the formulation of HO-MoM on NURBS curve is introduced and the polynomial expressions are explicitly presented.The choice of polynomial basis function is evaluated in section 4. Finally, the numerical results is given and commented in the last section.

Scatterer with high curvature
Standard MoM and HO-MoM are commonly used in many applications and the so-obtained numerical results can be the most of the time considered as a reliable model to simulate the scattered field.However, where the geometry of the scatters is complex with high curvature component and with no particular symmetry [9], the MoM approaches often fail.given by Cp = 1 r (where r is radius of the osculating circle), must be evaluated with respect to the wavelength λ of the electromagnetic wave.For high curvature scatterers (Cp 1 λ ), the scattering induces complex creeping wave phenomenons [10,12,13].In these situations, the numerical methods such as standard MoM and HO-MoM constitute a very slowly convergent algorithms, and an accurate simulation requires a mesh with a very high density of nodes in the curved area.A fully detailed analysis of the error induced by the curvature for different MoM approaches in 2 dimension is studied by Davis et al [11].
To overcome the convergence issue for curved scatterers, a more recent way is to explicitly take into account the curvature and the global geometry of the scatterer into the electromagnetic model, see [14] for instance.In this paper, we show that NURBS forms a relevant tool to model curved scatterers.

NURBS curve
Non Uniform Rational Basis Splines (NURBS) is a mathematical model for representing curves and surfaces especially for the complex geometry [15].Each point of the NURBS curve is computed by taking a weighted average of given control points.The control points determine the shape of the curve.Figure (2) illustrates the algorithm to construct a NURBS curve from a set of control points.
By the lack of numerical stability to determine its derivation, NURBS curve is needed to be transformed to Bézier format [7].The mathematical expression for a where P k are the control points, w k are the weights and The NURBS curve is here expressed for v ∈ ]0 1[.The basis functions used in this work are the modified Legendre polynomials, defined in u ∈ ]−1 1[ (section 4).Thus we need to convert the equation (1) to this interval by the relation u = 2v − 1 [7]: From (3), we can easily find the covariant unitary vector, defined as the derivation of the NURBS curve in u.
NURBS preserves the surface normal and radius of curvature at every point by adjusting the control points.Figure (3) shows the comparison of the meshing by the standard technique (straight line) and NURBS curve.Despite a larger mesh length, we notice that the NURBS model provides a far better accuracy.
Development of NURBS began in the 1950s, but the integration of the geometrical tools given by NURBS for the electromagnetic scattering models was first introduced in 1994 by Valle et al [16].

High-Order Method of Moments on NURBS Curve
The first step in MoM application involves the development of the integral equation of the system [17].Using Huygens' principle, the electromagnetic system is represented in the form of a Electric Field Integral Equation (EFIE) and a Magnetic Field Integral Equation (MFIE).In 2D problem, the vector field can be simplified to the scalar field according to Transverse Magnetic (TM) and Transverse Electric (TE) polarization.So, EFIE or MFIE can be used to model these two polarization problems.Better in numerical condition, we suggest to solve the TM polarization with EFIE and TE polarization with MFIE.
First and foremost, comparing some articles and books on propagation theory, let us note the inconsistent use of TM and TE polarization [18].So, it is necessary to mention that the convention applied in this article is the same as used by Peterson [19] or Gibson [1] in contrast to Tsang [20] or Naqvi [18].If the propagation vector k lies in x − ẑ plane, the TM/TE polarization refers to the electric/magnetic field being in ŷ direction as shown in figure (4).The position in the x − ẑ plan is given by the vector ρ(xx, z ẑ).
x y In the next section, we will briefly present the EFIE for TM polarization and MFIE for TE polarization in case of PEC object.The formulation for a dielectric can be obtained by combining of these two form.

EFIE for TM Polarization
The EFIE in TM polarization for a PEC object is presented as [1]: where µ is the permeability of free space, ω is the angular frequency of the wave, G(ρ, ρ ) is the 2D Green function, ρ and ρ are the source and observation points.The incident plane wave is given by E inc y (ρ) = e −ik(x sin θ i −z cos θ i ) .For the HO-MoM on NURBS curves, the current is expanded in arbitrary u coordinate as: where |a(u )| is jacobian of coordinate transformation, P n (u ) are polynomial basis function of order n and b n are the unknown coefficients.
Using the Galerkin procedure, the EFIE is transformed to the matrix equation m, n are the order of testing and basis function respectively and p, q are observation and source curve respectively.The component for each matrix are:

MFIE for TE polarization
The MFIE for TE polarization for a PEC object is presented as [1]: The incident magnetic plane wave is given by H inc y (ρ) = e −ik(x sin θ i −z cos θ i ) .The current is expanded in coordinate u as: The Galerkin procedure is used to transform the MFIE to the matrix equation: with the matrix V : When p = q and u = u , the matrix Z components are: [Z m,n ] p,q = −0.5 Cpq P m (u)P m (u)du (14) and for the other conditions: The matrices V and Z for EFIE and MFIE are constructed using Gaussian-Quadrature integration method [21].According to [22], the polynomials of order 2N − 1 can be integrated exactly by N points Gaussian quadrature rules.
The fundamental difference between equations ( 10) and ( 5) is the fact that the MFIE equation involves the derivative of the Green function.We could think that this point have a significant importance for the choice between a continuous and a discontinuous functional basis.

High-Order Basis function
Basis functions for MoM can be categorized to entiredomain and subsectional [1].Entire-domain basis functions are not practical and rarely used.Subsectional basis functions involves the discretization of the domain in to some subsections.The polynomials are introduced to approach the surface currents in each subsection.In the basic MoM, the polynomial are limited to lower order N = 0 (pulse) if the continuity is not required and N = 1 (triangle) if the continuity need to be imposed.High-Order MoM use higher-order polynomial N ≥ 2 as the basis functions.
To ensure the computational convergence, the basis function must to be in the same function space as the surface currents.Since the currents in EFIE equation (5) and MFIE (10) are in the integral form, the basis functions have to be integrable, but do they have to be continuous?The advantage of the EFIE for TM polarization and the MFIE for TE polarization is the extinction of the currents divergence from the integral equations so that the basis functions can be discontinuous (L 2 space).This is not the case in EFIE for TE polarization which requires the continuity of the basis.
Since the basis functions in our case are in L 2 space, a priory, all type of polynomials can be chosen.Hamilton et all [23] used the product of Jacoby and Legendre polynomials.Although they did not impose the continuity, the results are corrects for EFIE in TM polarization.
However, the continuity is still important tools to reduce the number of unknowns [2] and it is interesting to see the influence of the imposition of the continuity.Djordevic and Narsos [24] have defined three polynomial for this propose.We can also include the Bernstein polynomial in to this category.However, those polynomials are not orthogonal and leads to severely ill-conditioned matrices [2].The modified Legendre polynomials developed by Jorgensen seems to be the best solution.They ensure the orthogonality and can be used to impose the continuity.N -order modified Legendre polynomials are given by: where n = 1, 2 • • • , N and u ∈ ]−1 1[ and L n are the Legendre polynomials that can be expressed in Rodrigues' formulas: Figure (5) shows the 5 th order modified Legendre polynomials introduced in one subsection.The continuity is imposed by adjusting two first order of polynomials with neighboring subsection as shown in figure (6).The higher order term (n ≥ 2) are zeros in the extremities and do not contribute to continuity of the currents.In this article, we will study the two case: discontinuous and continuous basis function.The pulse basis function are added as lowest order (N = 0) for the discontinuous basis function.For the continuous basis, the lowest order is the triangle (N = 1) function.

Numerical results
By definition, the notion of curvature is linked to the radius of the circle and can be analyzed as a function of the wavelength.It is justifiable to investigate the performance of HO-MoM+NURBS curve for this canonical geometry.The analytical solution of electromagnetic scattering for a circle is given by modal theory (Mie theory [25]) and is expressed by the sum of Bessel and Hankel functions [1,26,27,28,29] (the series will be truncated at 30 terms).The surface currents and RCS computed by this analytical method, called in the sequel "exact solution", will be used as the reference in our analysis.
The HO-MoM in other side is a numerical method whose precision depends on many factors: mesh length (L), basis function orders (N ), basis function types (discontinuous or continuous), wavelength (λ) etc.To investigate the accuracy of HO-MoM from these factors, we apply the Normalized Root Mean Square (NRMS) error criterion: where N S is the number of sampling points (360 points in this work), X are the surface currents or RCS computed by HO-MoM and X ref are the surface currents or RCS computed from analytical solution.

Surface Current
To illustrate how HO-MoM+NURBS approximates the surface currents, we take the example of an incident wave (λ = r) which radiates a PEC circle as shown figure (7).The surface is represented in 6 NURBS curves with the constant mesh length L.
Incident wave   Zooming the figure (8) for the order N = 2, we see clearly the effect of the continuity imposition (figure (9)).In facts, it reduces the oscillation in extremities of the sections and gives the best approach of the surface currents in this zone.
The similar results are found for TE polarization as shown in figure (10).Obviously, we see the effect of the imposition of the continuity to reduce the oscillation in extremities zones.
To study more precisely the influence of the basis func- First, it sounds logical to say that the accuracy of surface currents is better for the higher order basis function and smaller mesh length.In a global way, the figure (11) seems to confirm this intuition.Nevertheless, in the TM case and when L = 1 dL = λ/10, the best result is obtained for an order N = 2. Thus in this situation, an optimal order exists.
Not surprisingly, the discontinuous basis functions give more accurate results than the continuous basis especially for the lower orders.In fact, by imposing the continuity, we loss one degree of freedom in the computation of the coefficients of basis function.It reduces the oscillations in the vicinity of the segment extremities, but paradoxically, it decreases the accuracy in global term.
In the figure (12), we focus on the case L = 1 dL = λ/10, and we evaluate the influence of the ratio wavelength/radius.For TM polarization, in the same ways as the figure (11), the error curves reach a minimum points.The position of this points depends on the wavelength (N = 1 for λ = r/2 and N = 2 for λ = 3r).For TE polarization, higher order seems to be better and no minimum appears.
In fact, the figure (12) highlights the effect of the curvature.In both cases (TM and TE), a higher curvature (lower wavelength to radius ratio) induces a loss of precision using the standard MoM approach (order 0) for the same wavelength to mesh length ratio (dL = λ/10).For the TE polarization, the increase of the HO-MoM+NURBS order clearly compensates the curvature effect.For the TM polarization, the curve-crossing let us think that he HO-MoM+NURBS order partly compensates the influence of the curvature, and a optimal order must be determine as a

Radar Cross Section
Once the surface currents have been estimated, the Radar Cross Section (RCS) of the scattered can be evaluated by a very well known integration process, see for instance [19].RCS in TM polarization for discontinuous and continuous basis functions is given in figure (13).It is important to notice a significant difference with regard to the surface current estimation: the lower orders do not generate oscillations.Indeed, the RCS computation is obtained by a  integration that considerably attenuates the oscillation phenomena.Comparing the order N = 2 for these both basis, we note that discontinuous basis function give better results than the continuous basis.For N = 5, the RCS for two basis give the good precision.The same behaviors of RCS are obtained for TE polarization (figure (14)).
In figure (15), we show the accuracy of RCS as a function of the order (N ), mesh lengths (L) and the basis function types.For TM case, we can underline the fact that RCS error curves reach a minimum accuracy point contrary to the TE case.In a global way, the computation of the RCS involves a far better precision than one obtained for the surface current.The speed of convergence in the RCS case is significantly greater, and the order of the HO-MoM+NURBS approach does not need to be so high, especially with discontinuous function basis.

Computational Time
The numerical simulations in this work were done using a standard personal computer.To highlight the influence of order (N ) and the mesh length (L) upon the CPU time consumption, the figure (17) presents the CPU time for TM case in λ = r.To obtain a standard evaluation, the computation times are normalized by the time of the standard MoM (N = 0) with 1dL.
In a first approximation, figure (17) suggests that the computational time is given by the following heuristic formula:  where A is the constant for the CPU.In fact, the integration point number is one of mainly factor that determine the computational time.Gaussian Quadrature method allows us to take a few point number and decrease significantly the time even for the higher order basis.
Nevertheless, this paper mainly deals with the accuracy of our HO-MoM+NURBS approach.To avoid biases induced by numerical acceleration process, no such acceleration was applied.In addition, it is noteworthy that the integrations for HO-MOM can be easily realized in an independent way.So the parallel computing significantly modify the estimation the time consumption in the present case.This is the reason why our analysis of the CPU time consumption probably needs to be completed by further stud-

Conclusion
Our paper presents the HO-MoM combined with NURBS segments for the scattering by a circle (2D problem).In this canonical configuration, we show that the optimal parameters for the HO-MoM with NURBS depend on what the observer want to evaluate.The high order basis is not always the best in the convergence speed.In many situation, the optimal order exists and do not exceed several units.Even  for a basic canonical scatterer, the optimal order (N ) is significantly influenced by the polarization, the mesh length (L), the wavelength (λ), the curvature and the distance of the observation (far versus very near field).
We have also shown that the influence of the polynomial basis (discontinuous/continuous) is only relevant where the order is null or very weak.
Finally, to establish the best compromise in practice, we have introduced an evaluation of the computational time and we have pointed the fact that the HO-MoM approach

Figure 4 :
Figure 4: Convention of TM and TE polarization

Figure 11 :
Figure 11: Log plot of surface current error versus the order and L for λ = r

Figure 15 :
Figure 15: Log plot of RCS error versus the order and dL for λ = r

Figure 16 :
Figure 16: Log plot of RCS error versus the order and λ for L = 1 dL = λ/10