Non Uniform Rational B-Splines and Lagrange approximations for time-harmonic acoustic scattering: accuracy and absorbing boundary conditions

ABSTRACT The paper aims to evaluate the performance of the Lagrange-based finite element method and the non-uniform rational B-splines isogeometric analysis of time-harmonic acoustic exterior scattering problems using high-order local absorbing boundary conditions, in particular based on the Karp’s and Wilcox’s far-field expansions. The analysis of accuracy and convergence of both methods is achieved by observing the effect of the order of the approximating polynomial, the number of degrees of freedom, the wave number, and the absorbing boundary conditions tuning parameters. It is concluded that, regardless of the polynomial order, IGA provides a higher accuracy per degree of freedom compared to the traditional Lagrange-based finite element method.


Introduction
Various numerical methods have been designed over the years for wave propagation arising in solid mechanics, geophysics, meteorology, acoustics and electromagnetic applications.Broadly speaking, the numerical techniques can be classified as domainand boundary-based approaches.Before focusing on some specific numerical methods related to the paper, we refer to the references [1][2][3][4][5][6][7][8], where the interested reader can find some overviews on general numerical methods for solving acoustic scattering problems.
Spectral finite element method is one of the most commonly used methods for wave propagation [9,10].The distribution of the domain nodes is such that oscillations that occur because of the Runge phenomenon can be reduced [11].This method requires fewerdegrees of freedom (DOFs) per wavelength compared to the conventional Finite Element Method (FEM) [12].Using vector and parallel computing algorithms [13], it easily reduces the storage requirements and the computational complexity of this method.The Chebyshev polynomials [14] are typically employed to minimize the dispersion error [15].Higher accuracy is obtained compared to the low-order p-FEM (quadratic) [16,17].The major drawback of SEM is its difficulties to deal with complex geometries.The finite difference method (FDM) based on the Taylor series expansion is a classical way to solve the wave equation.The terms are truncated to arbitrary number and the dominant power of the truncated terms determines the accuracy.The Cartesian grids are usually necessary to obtain the solution unless generalized finite differences (GFD) are used.The GFD is a meshless technique used on the domain, eliminating the mesh generation and the numerical quadrature.This results in requiring a high number of degrees of freedom to obtain higher accuracy.In the case of curved complex geometries which is a major characteristic and difficulty of scattering problems, this method cannot produce accurate results due to the well-known staircase effect [18].The semi analytical finite element (SAFE) [19,20] is another way of combining the advantages of numerical and analytical methods.This technique which uses Fourier Transforms to recover time-domain analysis was used to study Lamb wave propagation, whose behaviour is independent of the direction of propagation [21].For a given accuracy, the computational cost for SAFE is less than that of FEM.The disadvantage here is that the wave propagation over complex geometrical features cannot be handled.When the geometries are complex, it is needed to make some of the approximations and assumptions which affect the accuracy.In these cases, the FEM becomes more useful.
Boundary Element Methods (BEM) are another way of solving the infinite domain problem.The method employs discretization of the boundary, thus reducing the ddimensional problem to the ðd À 1Þ-dimension.The advances in the computing power has catalysed the use of this method.The fundamental solution represented by the Green's function is applied at the boundary and implicitly satisfies the Sommerfeld's radiation condition at infinity.The method can be used when the domain has a complex geometry but cannot be applied usually in heterogeneous media where the Green's function is not known.The resulting system of equations from BEM results in a dense and large linear system to resolve since the operator is non-local.Thus, this drastically increases the computational effort and memory storage required [8].In addition, going to very high-order BEM schemes, most particularly for high-frequencies, still remains unclear, in particular because it is usually combined with fast iterative Krylov subspace solvers [22], preconditioners [23,24] and matrix compression algorithms (see, e.g.[25,26]).
The domain-based approaches, in particular the FEM, rely on discretizing the bounded d-dimensional domain with non-overlapping regions and employing polynomial basis functions to approximate the unknown fields.The advantage of the FEM is that it can handle complex geometries, anisotropic properties, and requires no fundamental solution to represent the unknown fields.Depending on the differential operator, the matrices obtained are sparse, thus requiring less storage space while being very large and indefinite positive.In practice, in the difficult high-frequency regime, the FEM requires the use of very fine meshes and/or to increase the polynomial order to capture the fast oscillations of the unknown [27], increasing then the computational cost.In addition, the linear system solution is non-trivial since the system is highly indefinite positive and makes the Krylov subspace solvers diverging or slowly converging [28].Specific Domain Decomposition Methods (DDM) such as the optimized Schwarz's DDM [29] are then necessary for solving large scale engineering and industrial problems.Finally, let us remark that another way to use both the FEM and BEM is to combine these two methods through FEM-BEM coupling strategies.Another method that shares the advantages of both the FEM and the BEM is the scaled boundary finite element method (SBFEM) [30,31].Like the FEM, it does not require any fundamental solution and like the BEM, only the boundary is discretized.Further, the radiation boundary condition can exactly be satisfied.
Using Delaunay triangulation, in one minute 3 billion tetrahedral meshes can be generated using one machine.However, it usually does not provide quality meshes for analysis purposes [32].The meshing of the computer-aided design (CAD) generated models is time consuming, especially optimization of the mesh to obtain the quality mesh.The original geometry has to be consulted for the mesh adaptation or remeshing procedure in case of refinement of the solution on the given domain.Thus, the mesh generation process becomes tedious [33].Furthermore, mesh generation in FEM takes 80% of the total analysis time, due to the lack of a direct connection with the models generated by CAD platforms.To attempt to alleviate these difficulties, isogeometric analysis (IGA) was introduced in 2005 by Hughes et al. [34].The idea was to construct one geometrical representation for all levels of mesh refinement.In IGA, the NURBS basis used for the CAD models is also used for the solution space.The mesh directly interacts with the geometry and the remeshing process can be fully automated.This property combined with unique refinement possibilities makes IGA an attractive platform for shape optimization of devices relying on wave propagation phenomena [35][36][37].One of the difficulties is that the CAD model provides only surface mesh and an extension of the surface mesh to a volumetric mesh is non-trivial.The mesh requires a tensor product description and also, inhomogeneous Dirichlet boundary conditions need to be applied weakly.Adaptive meshing is another issue but can be addressed by using the other variants of NURBS, viz., T-splines and PHT-splines.Within this framework, NURBS are used as basis functions but are truncated.The performance of B-splines FEM and adaptive PHT-spline IGA for solving exterior time-harmonic acoustic problems were studied in [38,39].In particular, it was shown that the pollution error is well controlled for a fixed discretization density and order of the basis functions p � 3.
The k-refinement strategy offers robust results and optimal convergence rate which is also achievable using conventional FEM at higher computational cost, but the convergence rates remain the same as that of higher-order FEM.The improved accuracy can be better achieved in k-refinement than with p-refinement for vibrations and wave propagation problems, where the solutions are smooth [40].Therefore, in this study, we employ k-refinement IGA.
When using the IGA or the FEM for exterior problems, the unbounded domain must be truncated with an artificial/absorbing boundary condition to approximate the Sommerfeld's radiation condition and to ensure that the scattered wave is outgoing to the computational domain boundary while minimizing the spurious unphysical reflection.The two common choices for satisfying this on the truncated boundary are by employing (a) a perfectly matched layer (PML) [41][42][43][44][45][46] or (b) an absorbing boundary condition (ABC) [47][48][49][50][51][52][53][54].In PML, an additional absorbing layer is added that surrounds the computational domain.However, obtaining accurate results using PML usually demands to tune some specific physical and non-physical parameters which is far from being trivial [55].It is especially challenging to tune these parameters when no exact solution is available for the problem under consideration.Apart from these two truncation techniques, a single infinite element [56][57][58] is used instead of non-reflecting boundary conditions.The basis functions of this element are in terms of radial functions of the outgoing wave.The test or the weighting function is conjugate or unconjugate to the trial function chosen.Some of the examples of these functions are Bettess-Burnett unconjugated type, Burnett conjugated and Ashtey-Leis conjugated [59][60][61].For unconjugated type, the Gauss integration can be used to integrate the radial functions.For conjugated type, there is no oscillatory term and the integration can be performed analytically.The unconjugated types are less accurate at the far-field and offer higher precision at the near-field.In case of conjugated type, it is other way round.The accuracy of both methods deteriorates when the wave number increases.
The domain truncation introduces an error even at the continuous level which can be reduced by considering high-order ABCs.However, they are mostly tedious to implement and require the evaluation of higher order derivatives [62,63].Also, non-local highorder ABCs will result in fully populated linear system greatly reducing the sparsity of the FEM.Feng's higher order boundary condition [64] was studied within the framework of the FEM.It was inferred that the framework yielded accurate results for a larger domain of computation, which directly increases the computational cost.Another approach is to use enrichment techniques [65,66], but from the numerical study it is shown that the truncation error had higher influence than the discretization error.Khajah et al. [67] employed the recently developed higher order local ABC based on the Karp's far-field expansion (KFE) and Wilcox's far-field expansion (WFE) within IGA for low to high frequencies for 2D and 3D problems, respectively.Unlike other ABCs, the salient feature of these higher order boundary conditions is that they do not require higher order derivative evaluations and hence impose no additional constraint on the regularity of the basis functions.This implies that the same basis functions can be used to estimate the field and the auxiliary ABC functions [67].This is accomplished by introducing two families of auxiliary variables on the fictitious boundary.More recently, the KFE-ABC was used successfully to eliminate the domain truncation error in IGA collocation context [68].To the best of the author's knowledge, the conventional Lagrange-based FEM with KFE and WFE has not been studied in the literature.
The purpose of this paper is to investigate the performance and to compare the accuracy of the higher order FEM with the NURBS-based IGA with and without domain truncation error.This is accomplished in three steps: (a) first, we artificially eliminate the domain truncation error in the numerical results by imposing Robin and BGT-2 absorbing boundary conditions on the artificial boundaries in duct and circular cylinder problems and, (b) later, we apply BGT-2 and KFE-ABC on the artificial boundary, and (c) compare numerical errors considering the full exact solution.In order to have a fair comparison between FEM and IGA, we compare them on the basis of the number of degrees of freedom for a given basis order.We also use tensor product meshes with equidistant nodes, Gaussian quadrature integration rule with identical number of Gaussian points and equal number of degrees of freedom.In the third part, we apply Wilcox's far-field expansion to 3D problems.
The paper is organized as follows.Section 2 presents the governing equations, the corresponding weak forms and a brief overview of the BGT-2, KFE and WFE ABCs.Section 3 presents an analysis by eliminating the domain truncation error, followed in Section 4 by a systematic parametric study on the numerical solution for the circular cylinder from low to high frequencies, that demonstrates the efficiency of the KFE-ABC.In Section 5, the higher order WFE-ABC is applied to 3D problems for illustrating the method.Major conclusions are presented in the last section.

Governing equations-ABCs
We define Ω À as a d-dimensional scattering bounded domain of R d , with boundary Γ :¼ @Ω À .We introduce the corresponding exterior unbounded domain of propagation Ω þ :¼ R d nΩ À .Then, we assume that an inhomogeneous Neumann boundary condition is prescribed with a function g on the boundary Γ and solve the scattered wave field u in the following Boundary-Value Problem (BVP): find u such that where Δ is the Laplacian operator, Ñ is the gradient operator and n is the outwardly directed unit normal vector to Ω À .The wave number k is related to the wavelength λ by the relation: λ :¼ 2π=k.Denoting by a � b the Hermitian inner-product of two complexvalued vector fields a and b, then the last equation of system ( 1) is known as the Sommerfeld's radiation condition [5,7,69,70] at infinity, ensuring the uniqueness of the solution to the BVP.We also define the discretization density n λ as the number of degrees of freedom per wavelength along each spatial direction.
In domain-based computational methods, the Sommerfeld's radiation condition at infinity is approximately satisfied by artificial truncation of the computational domain at a fictitious boundary � and by replacing the Sommerfeld's condition with an ABC.The resulting bounded computational domain with boundary Γ and � is denoted by Ω.In what follows, we consider three types of ABCs to truncate the computational domain.We use the Dirichlet-to-Neumann map to describe the general form of the ABC: Now, we write the weak form of Equation 1 and find u 2 H 1 ðΩÞ such that aðu; vÞ ¼ ,ðvÞ "v 2 V; (3) ,ðvÞ ¼ ð Γ gv dΓ: (5) In the simplest form, the low-order Robin boundary condition Bu ¼ À iku can be used to approximate the Sommerfeld's radiation condition.In the present paper, we analyse the BGT-2 and Karp's far-field expansion ABCs in 2D, the Wilcox's far-field expansion ABC in 3D, and compare the performance of FEM with IGA, both with and without eliminating domain truncation error.

2D BGT-2 ABC
The symmetrical second-order Bayliss-Gunzburger Turkel ABC (BGT-2) can be defined by using Bu :¼ @ s ðα@ s uÞ À βu, and where κ is the curvature over the surface � and the curvilinear derivative is @ s .The BGT-2 boundary condition can be written in the weak form ð Hence, the bilinear form aðu; vÞ is When the domain is truncated with a circular fictitious boundary � of radius R 1 , the curvature is then κ ¼ 1=R 1 .
We define the following space |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } 2L times

;
to rewrite the form a to solve (3) as: find ðu; where We denote the number of terms on the artificial boundary with m.We add 2m columns to include the Karp's expansion to the linear system by using ð Next, we write 2m additional rows to satisfy ð where and Finally, we add 2mðL À 1Þ rows to incorporate the recurrence formulas, for l This procedure increases the size of the stiffness matrix by 2mL to simultaneously solve the unknowns introduced through the far-field expansion.We note that KFE-ABC requires the evaluation of the first derivative only, similarly to BGT-2.The matrix structure is discussed in detail in Ref. [67].
The far-field pattern can be obtained from the numerical scattered field at the fictitious boundary �.Indeed, the far-field is known from the angular integral function f 0 involved in the dominant term of the asymptotic expansion of u when the radial variable r tends to infinity as uðr; θÞ � e ikr ffi ffi r p f 0 ðθÞ: The efficient and accurate numerical evaluation of the far-field is then as follows.First, we compute the volume solution u on the computational domain Ω, e.g. with IGA.This means that the mesh onto � is a priori not uniform.Then, we sample the values of the trace u � at m uniformly distant points on � to get u j� , 1 � j � m.Using this, the far-field can be evaluated efficiently by using a Fourier series expansion.More concretely, the points on the boundary should be equidistant in order to compute the far-field pattern by using b q ðÀ iÞ q e iqθ ; (19) where b q ¼ c q =H ð1Þ q ðkRÞ, and H ð1Þ q is the Hankel's function of first-kind and order q.We set ðc q Þ q¼À m=2;���m=2À 1 , as the discrete Fourier transform vector of the numerical solution on the artificial boundary, u � being interpolated at the points ðθ j� ; u j� Þ on �, and m being the (even) number of points on the fictitious boundary �.The exact solution at the far-field is obtained from where R 0 is the radius of a circular scatterer centred at origin and Jn 0 designates the derivative of the Bessel's functions of the first-kind.The parameter 2 n is the Jacobi symbol: 2 0 ¼ 0, and 2 n ¼ 2 for n ¼ 1; 2; 3; . . .

3D Wilcox's far-field expansion ABC
The higher order absorbing boundary condition for the three-dimensional scattering problem can be given by the Wilcox's expansion.The sphere is of radius R 1 with surface � :¼ S. Let us introduce The Wilcox's expansion has only one recurrence formula.Thus, only the radial derivative continuity equation is sufficient to define the boundary conditions.In the above expressions, Δ S is the Laplace-Beltrami operator in spherical coordinates.We define  3) is: find ðu; F 0 ; . . .
For m degrees of freedom on the artificial boundary, we add m equations to the linear system of equations based on ð and fix mðL À 1Þ rows for the recurrence formula, for l The implementation is similar to the KFE-ABC in 2D, but the Wilcox's expansion has only one recurrence formula.

Comparison between FEM and IGA with no domain truncation error
In this [p = [htd[[n[nn[tion, we study the performance of high-order FEM and compare it to IGA by solving two benchmark problems: (a) the propagation in an infinite waveguide along the x-axis and (b) the scattering by a sound-hard circular cylinder subject to an incident plane wave.By adopting modified exact solutions, the domain truncation error is removed.Therefore, the comparison performed in this section is an indicator of the performance of IGA and FEM in solving time-harmonic acoustic problem regardless of the ABC used to truncate the domain.For discussing the results, the relative L 2 -error is employed where f h is the numerical solution and f ex is the exact solution.Such comparison sheds light on the levels of the pollution error of each method and clarify the minimum error one can expect for a given refinement and frequency without domain truncation error.
Before that, we provide a few elements about the Lagrange and NURBS interpolation basis.

Basics on NURBS and Lagrange interpolation
In this study, the domain is discretized with Lagrange finite element elements and nonuniform rational B-splines functions.In case of Lagrange finite element method, the finite element basis construction requires the definition of the nodes and the nonoverlapping regions are defined by its nodal coordinates.The elements are tied together by C 0 continuity.Figure 1 shows the Lagrange polynomials of degree p ¼ 5 for the onedimensional domain: 0 � � � 1.The domain is divided into two equal parts.From Figure 1, it can be seen that the Lagrange basis functions are interpolatory and the basis can take both positive and negative values.
On the other hand, the isogeometric analysis uses non-uniform rational B-splines (NURBS) as the basis functions.The construction of NURBS basis functions requires a non-decreasing sequence of parameter values called the knot vector, � i � � iþ1 ; i ¼ 0; 1; � � � ; n À 1), the degree of the curve p, the control points P i , and the weighting of each of the control points w.For IGA, there are two notions of mesh, the knot vector denotes the physical mesh and the control points designate the control mesh.The ith B-spline basis function of degree p, denoted by N i;p is defined as A NURBS curve of pth degree is defined as follows: The control mesh does not conform to the geometry, but is instead like a scaffold of the geometry.Unlike the Lagrange basis functions, the NURBS basis functions are non interpolatory in nature.The entire geometry can be a single or multiple elements or patches.In case the entire geometry is divided into subdomains, they are joined together with C 0 continuity, like the FEM.In order to illustrate the difference, the B-Splines functions are shown in Figure 2a -2b, for order p ¼ 5 with C 4 and C 0 , respectively.The one-dimensional line is divided into two patches, on Figure 2a, with the knot vector, Ξ ¼ f0; 0; 0; 0; 0; 0; 0:5; 1; 1; 1; 1; 1; 1g for the domain divided into two elements.From Figure 2a, it is observed that the basis are discontinuous at the boundary and C 4 continuous at the patch boundary, i.e. at � ¼ 0:5. Figure 2b shows the NURBS p ¼ 5 basis with a knot vector Ξ ¼ f0; 0; 0; 0; 0; 0; 0:5; 0:5; 0:5; 0:5; 0:5; 1; 1; 1; 1; 1; 1g.The knot 0.5 is repeated p times to obtain C 0 continuity at that location.Thus, the patches are discontinuous at the location � ¼ 0:5.We see that the NURBS functions are pointwise positive.Thus, the geometric representation using NURBS can be exact at all levels of discretization, whereas the FEM uses piecewise continuous basis in order to approximate the geometry.

Waveguide problem
In the first example, we analyse the case of an infinite waveguide along the x-axis.To compare the performance between FEM and IGA, we consider the truncated domain Ω ¼ ½0; 2� � ½0; 1� (see Figure 3).We assume rigid lower and upper walls and denote the outward unit normal vector by n.The domain is discretized with Lagrange and NURBS basis functions with order p, for a discretization density n λ .We solve the Helmholtz problem to find the acoustic pressure u in the duct with the following boundary conditions in Ω; @ n u ¼ cosðcπyÞ; on x ¼ 0; @ n u À Bu ¼ 0; ðBu ¼ À ikuÞ; on x ¼ 2; @ n u ¼ 0; on y ¼ 0; 1; where c 2 N is the mode number.The inlet is subject to inhomogeneous Neumann boundary conditions and absorbing (and transparent for c ¼ 0) boundary condition on the outlet boundary (x ¼ 2).Since the boundaries at y ¼ 0; 1 are assumed to be perfectly rigid, the normal derivative of the acoustic pressure vanishes on these boundaries.The exact solution of problem (31) with ABC is given as [71] Figure 2. NURBS basis functions of polynomial order p ¼ 5 in one dimension discretised by using two patches, defined by the knot vector: (a) Ξ ¼ f0; 0; 0; 0; 0; 0; 0:5; 1; 1; 1; 1; 1; 1g.The basis are C 4 continuous at 0.5 and (b) Ξ ¼ f0; 0; 0; 0; 0; 0; 0:5; 0:5; 0:5; 0:5; 0:5; 1; 1; 1; 1; 1; 1g.The basis are C 0 continuous at 0.5.
where k x ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi k 2 À ðcπÞ 2 q and the coefficients A 1 and A 2 are obtained from

� � :
We get the expression of the modified exact solution by solving this 2 � 2 linear system.Since this is a modified exact solution which excludes the ABC error on the duct outlet, this provides the possibility to study the finite-dimensional approximation regardless of the ABC applied.The duct cut-off frequency is c cutÀ off ¼ k=π.The wave in the duct propagates when c � c cutÀ off , and it represents evanescent modes when c > c cutÀ off .The real parts of estimated FEM and IGA solutions for k ¼ 40, c ¼ 2, p ¼ 3, n λ ¼ 10 are presented in Figure 4a and 4b respectively.The corresponding absolute errors are compared in Figure 5a and 5b, respectively.The geometry is exactly represented in both cases, but the error obtained by using the IGA is about two orders lower than that of FEM.This could possibly be attributed to the smooth and higher order continuous NURBS basis functions.
The evolution of the relative L 2 -error with discretization density n λ for the wave number k ¼ 40 and according to the wave number k for n λ ¼ 10 are shown in Figure 6a -6b, respectively.In both cases, the mode number c is taken as 2. The results are shown for different orders of basis functions, p ¼ 1; � � � ; 5. The slope of the graphs are noted with s in the legend.It is seen that the error of the FEM and the IGA are identical when p ¼ 1, as expected.Indeed, the IGA basis functions are then not different from linear Lagrangian basis functions.It is evident that IGA yields more accurate results per DOF for a given basis order.For p > 1, the relative error in IGA decreases sharply, most particularly for p � 3. The increase in the error with respect to k is observed for both the Lagrange and IGA basis functions, when p � 3.This is an indicator of the pollution error.Interestingly, the pollution in IGA seems to be (almost) fully under control for orders p � 3 and a fixed discretization density n λ .This could be attributed to the higher order and globally continuous basis functions employed in the IGA.

Scattering by a sound-hard circular cylinder
We consider now the acoustic scattering of a circular cylinder with radius R 0 ¼ 1 centred at the origin.The BGT-2 type boundary condition is applied on the fictitious circular boundary at radius R 1 ¼ 2 centred at ð0; 0Þ shown in Figure 7.We consider the Neumann boundary conditions at the scatterer boundary Γ.The BVP is then given by: find u satisfying on ΓðR 0 ¼ 1Þ; @ n u À Bu ¼ 0; on �ðR 1 ¼ 2Þ: (33) We apply the BGT-2 boundary condition on the fictitious boundary � following the weak form developed in Section 2.1 and consider an incident plane wave u inc ðxÞ ¼ e ikd�x ,  where d is the incidence direction d ¼ ðcosðθ inc Þ; sinðθ inc ÞÞ T and θ inc is the scattering angle.Since this problem is symmetric, we can fix the incidence direction to d¼ ð1; 0Þ T .The modified exact solution can be obtained as a summation of modal solutions [38] where Here, H where We plot the real part of the numerical solutions obtained with FEM and IGA for k ¼ 40, p ¼ 3, n λ ¼ 10 in Figure 8a -8b, respectively.The corresponding absolute errors are also plotted in Figure 9a for FEM, and Figure 9b for IGA.The average absolute error obtained for FEM and IGA are in the range of 5 � 10 À 3 and 1:5 � 10 À 4 , respectively.Thus, the wave field is properly estimated over the domain considering its frequency and selected discretization density.Equal number of degrees of freedom (397 � 64) were used for both FEM and IGA.The error of IGA is considerably lower which is consistent with the results obtained in Section 3.2.This is an indicator that IGA shape functions can conform better to the high oscillatory nature of the solution.The evolution of the relative L 2 -error with discretization density n λ in FEM and IGA is shown in Figure 10 for k ¼ 40, p ¼ 1; � � � ; 5.Both FEM and IGA meshes are structured and an identical number of radial and circumferential degrees of freedom were used in both methods.It is clear that IGA yields much higher accuracy per DOF for a given basis order and discretization density.The error levels shown in Figure 11 are the minimum we can expect for the given wave number, basis order and discretization density regardless of the ABC applied to truncate the computational domain.This yields an improved computational cost when using IGA since an equivalent accuracy can be achieved with lower refinement.The evolution of the relative L 2 -error with wave number k is shown in Figure 11, for p ¼ 1; � � � ; 5. Again, the BGT-2 type ABC is applied on the fictitious boundary at R 1 ¼ 2 and the contribution of the domain truncation error is eliminated from the error calculation by adopting the modified exact solution given in [35].Here, we can clearly see the advantage of IGA over FEM in solving high-frequency problems.The pollution error of IGA is really lower than that of FEM for the same order p and same number of degrees of freedom, providing excellent accuracy in mid-to high-frequencies.The low pollution error of IGA was also observed in [38] for very high frequencies but was not compared to the pollution error in FEM before.

2D acoustic scattering with KFE-ABC: from low to high frequencies
In this section, we solve the acoustic scattering from a circular cylinder problem again after applying now the Karp's far-field expansion boundary condition on the fictitious boundary �.The boundary value problem reads: We adopt the weak form for KFE developed in Section 2.2, where B is a compact notation for representing the KPE ABC.The parameter NT denotes the number L of terms of the Karp's expansion.We compare the performance of FEM and IGA in the very low-, midand high-frequency regimes.First, we consider the wave number k ¼ 5 and compute the far-field pattern from the numerical solution of the scattered field by using (19).We calculate the relative L 2 -error in the far-field with the exact solution (21).The evolution of the relative L 2 -error in the computational domain and the far-field with discretization density n λ and the number of Karp's expansion terms ðNTÞ are compared in Figure 12a -12b, respectively, for both the FEM and IGA.The comparison is done a wave number k ¼ 5 and for the basis order p ¼ 5. Unlike the analysis performed with the modified exact solution, the relative error of the FEM and the IGA are not significantly different for the selected range of the discretization density and the number of Karp's expansion terms.This is a clear indicator that for k ¼ 5 and basis order p ¼ 5, the ABC error is dominating the numerical solution for NT � 7.Moreover, with IGA, higher accuracy can be achieved by adding more terms in the KFE-ABC such that the accuracy matches that of the interior domain.
Next, we compare the relative error in L 2 -norm of FEM and IGA in the interior domain for k ¼ 10 and p ¼ 5, with varying discretization density n λ and number of Karp's expansion terms in Figure 13 .It is evident that the error obtained using IGA is smaller compared to FEM.As expected, the optimal number of Karp's expansion terms NT required depends on the accuracy of the numerical method in the interior.Therefore, it is a function of the wave number k, basis order p and discretization density n λ .This can be confirmed by observing the optimal number of terms required in the FEM and IGA in Figure 13 .For the discretization density n λ ¼ 5, we see that the accuracy is not improved even after adding more terms.This is true for both methods and can be attributed to the coarse discretization in the interior.
However, increasing the discretization density n λ improves the accuracy of the FEM and IGA in the interior which should be matched by adding more ABC terms.With higher discretization densities, such as n λ ¼ 25, we remark a similar convergence for both methods up to a certain number of Karp's expansion terms (NT ¼ 10).Beyond, adding more terms only increases the accuracy in IGA but not in the FEM.The accuracy on the ABC is matching the accuracy of the FEM but not the IGA at NT ¼ 10 and n λ ¼ 25.As a result, the IGA error is still dominated by the ABC.Hence, increasing the number of Karp's expansion terms to NT ¼ 20 consistently reduces the error in the IGA.The continuity of the NURBS basis function used in this study is C pÀ 1 , where p is the order of the basis function.If we keep the number of elements fixed and reduce the continuity to C 0 , then the total number of degrees of freedom will increase to reduce the continuity which is sub-optimal.The discretization density n λ in this study is defined based on the number of degrees of freedom as a fair basis for comparison.Let us consider a mesh with 37 � 10 degrees of freedom of order p ¼ 5 in IGA which is C 4 continuous corresponding to a mesh made of 120 (20 � 6) meshes as shown in Figure 14.
The number of meshes should be reduced to 24 (8 � 3) to maintain the number of degrees of freedom and reduce the continuity to C 0 in IGA.This reduction in number of meshes to maintain 120 degrees of freedom is shown in Figure 15 .
In order to better understand the effect of reducing the continuity in IGA, we plot in Figure 16 the relative L 2 -norm error calculated on the fictitious boundary for k ¼ 10, p ¼ 5, maintaining fixed the number of degrees of freedom for both the C 4 and C 0 analyses.
The evolution of the relative L 2 -error in the interior domain with wavenumber k in the FEM is depicted in Figure 17a when BGT-2 ABC is imposed on the artificial boundary.This is compared with only one Karp's expansion term NT ¼ 1 in Figure 17a, where n λ ¼ 10 and p ¼ 1; � � � ; 5.For the low-to mid-frequency range, the BGT-2 and Karp's expansion with one term NT ¼ 1 can be considered as equivalent.However, we note that BGT-2 or NT ¼ 1 not accurate enough for high-order FEM analysis.In order to find the minimum possible error achievable by increasing the number of Karp's expansion terms, we plot the evolution of the relative L 2 -error in the interior domain by using the modified exact solution in Figure 17b .It is clear that NT ¼ 1 provides an adequate  accuracy for p ¼ 1; 2; 3 since the error calculated with the modified exact solution is not lower for p � 3.However, when using a basis order p � 4, the accuracy in the interior is increased beyond BGT-2 and NT ¼ 1 accuracy and only one term of Karp's expansion is no longer enough.Hence, the error is dominated by the ABC error and we see that the error stagnates for p � 4 in Figure 17b .
The proposed boundary condition is now applied to the boundary R 1 ¼ 1:05.The radius of the scatterer is still R 0 ¼ 1.This makes the domain of interest in which the computations are required very small, thus decreasing the computational effort mously.Figure 18 shows the error obtained for FEM and IGA for p ¼ 5, k ¼ 5, NT ¼ 25 and R 1 ¼ 1:05.For the same discretization density and number of terms, the accuracy of IGA is higher than for FEM.
Since IGA is more precise than FEM, we now only retain IGA for the simulations.The advantage of Karp's ABC over BGT-2 to truncate the computational domain close to the boundary of the scatterer is demonstrated in Figure 19, where R 1 ¼ 1:2, enclosing the unit circular cylinder with R 0 ¼ 1.The increased truncation radius from R 1 ¼ 1:05 in Figure 18 to R 1 ¼ 1:2 in Figure 19 allows a large enough computational domain for comparing the absolute error obtained by using BGT-2 with the one obtained with Karp's expansion ABC.An incident plane wave for k ¼ 10 was considered and the numerical results were found by using the basis order p ¼ 5 and n λ ¼ 12.A considerable improvement is observed by switching to Karp's ABC.We note that the fictitious boundary R 1 can be considered very close to the scattering surface, leading to a highly reduced computational cost and an accurate solution.
Let us apply the BGT-2 and KFE-ABC to the circular boundary at R 1 ¼ 2, enclosing a 2D submarine-like shape.We consider the wave number k ¼ 50 and an incident plane wave.We truncate the domain at R 1 ¼ 2 and R 1 ¼ 1:5, without changing the discretization density.The estimated total field is shown in Figure 20 (R 1 ¼ 2) and Figure (R 1 ¼ 1:5) for the IGA basis order p ¼ 5, and five terms in the Karp's expansion, i.e.NT ¼ 5. Reducing the size of the computational domain yields less degrees of freedom.More precisely, the number of DOFs is reduced by 30%, passing from 25,600 (320 � 80) to 17,920 (320 � 56).We note that this reduction can be more significant for scatterers of smaller slenderness ratio since the truncation boundary can be placed very close to the  boundary of the scatterer.Getting Karp-like ABCs for more adapted shapes like for an elliptical boundary would be probably extremely efficient, most particularly for an IGA implementation.

3D acoustic scattering: numerical examples with IGA and WFE-ABC
To end, we now report some preliminary examples of simulations for 3D sound-soft obstacles.The frequency range studied here is not too high since the problems lead to very large size linear systems to resolve, needing some specifically designed DDM solvers.In addition, we only report results related to IGA which has been proved to be more accurate than FEM in 2D as well as for the WFE-ABC which is the extension to 3D of the accurate KFE-ABC.To start, we plot the real part of the scattered field around the sound-soft sphere in Figure 22 (R 0 ¼ 1), for p ¼ 2 and n λ ¼ 6 (16,530 degrees of freedom), where we choose the mid-frequency k ¼ 10.We retain NT ¼ 10 terms in the Wilcox's expansion set on a spherical boundary of radius R 1 ¼ 2.
Next, we demonstrate in Figure 23 the effect of choosing the basis order and the need to increase the accuracy of the ABC by adding more terms in the Wilcox's expansion, for the wavenumber k ¼ 2π.We observe that, for p ¼ 1, the accuracy of IGA is rather low for the discretization densities n λ ¼ 4 � � � 10.Therefore, implementing a more accurate ABC by increasing the number of terms in the Wilcox's expansion does not reduce the error.On the other hand, when the basis order p ¼ 4 is used, IGA is much more accurate in the  computational domain for the same discretization density.This improved accuracy is limited by the error related to the ABC employed.The error is not considerably reduced by increasing the basis order when n λ ¼ 10, and NT ¼ 1.This shows that the error of the high order analysis is bounded by that of the ABC.However, by adding more terms in the Wilcox's expansion, the error reduces to the minimal IGA accuracy for the selected basis order and discretization density, thus effectively avoiding the domain truncation error.
The real part of the scattered field by a sound-soft ellipsoid with semi-axes 1:5, 0:5 and 0:5, along the x-, y-and z-directions, respectively, is shown in Figure 24 .The incident plane wave has a wavenumber k ¼ 10.The analysis was performed by IGA with basis order p ¼ 5, for 16,530 degrees of freedom.The number of terms of the WFE-ABC is equal to NT ¼ 10.This shows that the scattering for a complex geometry can also be studied by using the high-order WFE-ABC with IGA.However, more investigations still remain to be performed to confirm this general behaviour.In addition, it would be extremely interesting, similarly to the 2D submarine case, to derive WPE-ABCs for a spheroidal fictitious boundary since this would lead to the possibility of diminishing the number of DOFs and then to better prospect high-frequency wave scattering.Finally, the full 3D methodology still needs to be investigated with more details concerning many aspects such as, convergence rate for complex problems, conditioning and spectral distribution related to the matrix defining the linear system.However, this is out of the scope of the present paper.

Conclusions
The performance of high-order FEM and IGA for solving exterior acoustic scattering problems is presented.Initially, the comparison of both methods was rendered by applying it to 2D problems, and then applied to numerically solve 3D cases.First, the   absorbing boundary condition on the fictitious boundary.The IGA results, compared with the FEM, show that IGA yields higher accuracy per degree of freedom and suffers less from pollution for a given basis order and discretization density.Therefore, if the domain truncation error is not dominant, a lower numerical error can be obtained using IGA.Next, KFE-and WFE-ABCs were implemented in both IGA and FEM, for 2D and 3D problems, respectively.The evolution of the error was studied with h-and prefinements.Additionally, the number of ABC terms is also analysed.Again, higher order IGA is shown to achieve a superior accuracy per degree of freedom compared to conventional FEM.It is opined that the error introduced by ABC can be conveniently reduced to match the accuracy of the higher-order method employed in the computational domain.All the findings in this paper indicate an improved accuracy in the NURBS-based IGA results, regardless of the domain truncation method used, and compared with the Lagrange finite elements.

Disclosure statement
No potential conflict of interest was reported by the authors.

ORCID
ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } L times ; the bilinear form aðu; vÞ involved in (

Figure 1 .
Figure 1.Finite element Lagrange basis functions of polynomial order p ¼ 5 in one dimension discretised by using two elements, with equidistant nodes.

Figure 5 .
Figure5.Absolute error ju ex À u h j for k ¼ 40, p ¼ 3, n λ ¼ 10.The numerical error of IGA is considerably lower for the identical basis order, number of DOFs, and integration points even when there is no geometrical errors.

ð1Þ m and H ð2Þ m
are the m th order first-and second-kind Hankel functions, respectively.The Neumann boundary condition is applied at the scatterer boundary Γ at R 0 and the ABC is applied on the fictitious boundary � at R 1 .This results in the linear system of equations with two unknowns a m and b m .

Figure 7 .
Figure 7. Circular cylinder of radius R 0 and the truncated cylinder around the domain of radius R 1 , the problem is solved in 2D and the computational domain considered is Ω.

Figure 8 .
Figure 8. Real part of the numerical solution for k ¼ 40, p ¼ 3, n λ ¼ 10.The numerical solution in both FEM and IGA captures high oscillatory wave field.

Figure 9 .
Figure 9. Absolute error ju ex À u h j for k ¼ 40, p ¼ 3, n λ ¼ 10.The numerical error with IGA is considerably lower for the identical basis order, number of degrees of freedom and the integration points.The contribution of the domain truncation error was eliminated by considering BGT-2 type ABC and modified exact solution.

Figure 10 .
Figure 10.Evolution of relative L 2 -error with n λ for k ¼ 40 for p ¼ 1; � � � ; 5. The BGT-2 type ABC is applied on the fictitious boundary and a modified exact solution is used to evaluate the error.

Figure 12 .
Figure12.Evolution of the relative L 2 -error in the domain with discretization density n λ and number of terms NT in Karp's expansion for k ¼ 5, and basis order ¼ 5.By increasing NT, it is possible to match the higher accuracy of the IGA in the interior.

Figure 14 .
Figure 14.(a) the mesh and (b) location of degrees of freedom of an analysis performed using p ¼ 5, C 4

Figure 15 .
Figure 15.(a) the mesh and (b) of degrees of freedom of an analysis performed using p ¼ 5, C 0 basis.

Figure 16 .
Figure 16.Relative error in L 2 -norm for p ¼ 5 and k ¼ 10 for IGA C 0 and C 4 NURBS basis.

Figure 17 .
Figure 17.Comparison of BGT-2 and KFE ABCs with a single term in mid-to high-frequency regime and the minimum relative L 2 -error in the interior domain for p ¼ 1; � � � ; 5, k ¼ 10; � � � ; 100 and n λ ¼ 10.More ABC terms should be added to match the accuracy of the interior solution when p � 4.

Figure 18 .
Figure 18.Relative L 2 -error for the artificial boundary truncated at 5% increase in radius compared to the scatterer radius.

Figure 20 .
Figure 20.Total wave field scattered by a 2D submarine with a KFE-ABC set on a circular fictitious boundary with R 1 ¼ 2.

Figure 21 .
Figure 21.Total wave field scattered by a 2D submarine with a KFE-ABC set on a circular fictitious boundary with R 1 ¼ 1:5.

Figure 22 .
Figure 22.Real part of the scattered field produced by a sound-soft spherical scatterer of radius R 1 ¼ 1, for p ¼ 5 and 16,530 degrees of freedom.The number of terms for the WPE-ABC is set to NT ¼ 10.The wavenumber is k ¼ 10.

Figure 23 .
Figure 23.Relative L 2 -norm error according to the discretization density n λ and number of terms NT in the Wilcox's expansion for the scattering problem by the unit sphere.

Figure 24 .
Figure 24.Real part of the scattered field produced by an ellipsoid using p ¼ 5 and n λ ¼ 6 for IGA.The number of terms for the WPE-ABC is set to NT ¼ 10.The wavenumber is k ¼ 10.
List of symbols used in this paper.