Implementation of Isogeometric Fast Multipole Boundary Element Methods for 2D Half-Space Acoustic Scattering Problems with Absorbing Boundary Condition

Isogeometric Analysis (IGA), which tries to bridge the gap between Computer Aided Engineering (CAE) and Computer Aided Design (CAD), has been widely proposed in recent research. According to the concept of IGA, this work develops a boundary element method (BEM) using non-Uniform Rational B-Splines (NURBS) as basis functions for the 2D half-space acoustic problems with absorbing boundary condition. Fast multipole method (FMM) is applied to accelerate the solution of an isogeometric BEM (IGA-BEM). Several examples are tested and it is shown that this advancement on isogeometric fast multipole boundary element method improves the accuracy of simulations.


Introduction
The BEM is a developing computational method for engineering and a commonly used numerical method in different engineering fields, including fluid mechanics, acoustics, electromagnetics and fracture mechanics. 1,2 As an alternative method to the finite element This is an Open Access article published by World Scientific Publishing Company. It is distributed under the terms of the Creative Commons Attribution 4.0 (CC-BY) License. Further distribution of this work is permitted, provided the original work is properly cited.

Isogeometric FMBEM for 2D Half-Space Acoustic Problems
This research extends the IGA-BEM into the 2D half-space acoustic problems with absorbing boundary condition. In order to improve the computing efficiency and decrease the storage, isogeometric fast multipole boundary element method (IGA-FMBEM) obtained by combining the IGA-BEM and FMM is applied for the 2D half-space acoustic problems in this paper. Several examples are tested and it is shown that this advancement on the IGA-FMBEM improves the accuracy of simulations.

Representation of Boundary
It is the key idea of IGA to bring together the fields of design and analysis into a unified framework by using parametric functions that are predominant in CAD. We focus on describing NURBS functions which are used to represent geometries in CAD. B-splines and NURBS can both describe the curves (or surfaces), and they are both parametric.
A single closed curve with B-spline basis functions of degree p is generated to represent the boundary, where the parameter ξ ∈ [0, 1] along the curve to be generated is introduced. A set of knots on the curve are given to generate a knot vector Ξ = [ξ 0 , ξ 1 , . . . , ξ m ], which is a nondecreasing sequence of real numbers given in parameter space, where ξ i ∈ R is the ith knot, m is the length of the knot vector, m = n + p + 1, and n is the number of basic functions. Turning attention toward B-spline and NURBS basis functions, the coordinate at any point on the curve constructed can be expressed by where vector P i denotes the coordinate at the control points, R i,p (ξ) are NURBS basis functions, defined as with where w i denotes a weight associated with the control point P i . B-spline basis functions N i,p (ξ) can be expressed as then, for p = 1, 2, 3, . . . .
Actually, it is hard for B-spline with lower order to model complex objects. For example, it is very difficult to model accurately circles and ellipses using low order B-spline. If we wish L. Chen et al.  to shift the curve closer or farther to the control point without introducing new control points and changing the position of the existing control point and, of course, not losing continuity.
A good method to achieve the aim is to use rational B-splines that are also called NURBS. Rational B-spline means the objects modeled by the spline can be expressed by rational polynomials. The solution is to introduce homogeneous coordinates into the function of spline. It associates a weight with each control point and we can change the curve by changing the weight of the corresponding control point. In Fig. 1, two complex structures, which are called octagonal column and eight blade column respectively, are constructed to show the flexibility and applicability of NURBS curves. Herein, the data about the control points and the associated weights are listed in Table 1.

Conventional BEM
We propose a collocational boundary element method using linear and quadratic Lagrangian functions for the Helmholtz equation in 2D. To make things simple, we focus on the acoustic scattering problem from noise barrier with absorbing material. The boundary integral equation (BIE) and normal derivative equation (NDBIE) based on Helmholtz equation for half-space problem can be expressed as, respectively and c(r)q(r, r 0 ) = ∂p i ∂n(r) + S ∂ 2 G(r s , r) ∂n(r s )n(r) − ikβ(r s ) ∂G(r s , r) ∂n(r) p(r s , r 0 )ds(r s ), Isogeometric FMBEM for 2D Half-Space Acoustic Problems where the coefficient c(r) is determined by the boundary geometry at the source point r, r s denotes the field point on the boundary, r 0 = (x 0 , y 0 ) is the acoustic source point outside the boundary, β is the acoustic admittance, p i is the incident acoustic pressure. The Green function G(r, r 0 ) can be expressed as wherer 0 = (x 0 , −y 0 ) is the mirror point of r 0 . In the conventional BEM, p(r s , r 0 ) in every discretized boundary element can be obtained by using Lagrangian interpolation function.
where m denotes the number of interpolation nodes in every boundary element, Φ denotes the Lagrangian interpolation function, r k denotes the kth nodal point. The Lagrangian interpolation function Φ for constant, linear, and quadratic boundary element are given, respectively 10,43 (1) For constant element m = 1 and Φ = 1 (2) For linear element where ξ means the local coordinate of the field point r s , and α denotes the position of interpolation nodes on the discontinuous element. When α = 1, Eq. (10) denotes expression of the interpolation functions for linear continuous element. (3) For quadratic element.
Similar as linear element, Eq. (11) denotes expression of the interpolation functions for quadratic continuous element when α = 1. geometrical node interpolation node Isogeometric FMBEM for 2D Half-Space Acoustic Problems

FMBEM for Half-Space Acoustic Problem with Absorbing Material
It is well-known that the conventional BEM generates a full coefficient matrix, which takes N 2 calculation operation for a problem with N degree of freedom. However, the use of the FMM avoids the direct solution of the coefficient matrix, and reduces the computing operation from N 2 to N or N log N , and so decreases significantly the required time of numerical solution. Actually, three different procedures of the FMM can be found according to the difference in the form of expansion of the Green's function, see Refs. 45, 46 and 47. For the full-space problem, only the tree structure containing the whole sound barrier structure is produced for the FMM operation. However, for the half-space problem, due to the existence of the mirrored domain, an extra tree structure needs to be produced in order to implement the FMM operation used for the solution of the boundary integral in mirrored domain. The boundary integrals in the two domains are solved separately. This method can be used successfully for the solution of the half-space problem, but it brings extra computing operation and increases the computational complexity. An alternative of this method is to create a mirrored source point, but not a mirrored integral domain. That means the calculation for the original boundary integral in the mirrored domain may be turned into the calculation for the boundary integral with a mirrored source point in the real structure domain.
Implementation process of the FMM is introduced in brief. In a first step, the tree structure is generated. A square containing the boundary S is introduced after discretization of the boundary. This square is divided into four equally sized child squares. Keep dividing a square in this way until the number of elements in that square is less than a specified number. The second step is called upward pass to calculate the multipole moments of each square. For this, the Green function is expanded as follows 7,8 : where r 1 s is very close to r s and |r 1 s r| |r 1 s r s |. The functions O n and I n are expressed as and where J n is the nth Bessel function, and (R, θ) denotes the polar coordinates of a vector. When r is far away from S 0 which is a subset of the boundary S, the boundary integral in Eqs. (6) and (7) can be rewritten as and Substituting Eq. (12) into Eqs. (15) and (16), we obtain the following two equations: where M n denotes the coefficient of the multipole expansion and is expressed as Then, the multipole moment of the higher level square is obtained by using the transfer operation of the multipole moment of the lower level square, that is called the M2M, as follows After obtaining the multipole moments of all divided squares, we need to implement the downward pass operation. The whole boundary integral is divided into two parts. One is called the nearfield integral that is solved by using conventional boundary integral methods.

Isogeometric FMBEM for 2D Half-Space Acoustic Problems
The other one is called the far-field integral that is solved by using the M2L and L2L transfer operation, as follows and where r 2 s is located close to S 0 , r 1 and r 2 close to r. Finally, we can obtain the following new expressions for the boundary integrals where r 3 is close tor. For more details on this, the reader is referred to Refs. 7 and 8.

Isogeometric FMBEM
In the implementation process of the IGA-BEM, the sound pressure at the boundary is interpolated by using the NURBS basis functions. Although a few geometric control points can generate very accurately a simple curve by using NURBS basis functions, large errors will be produced when the same number of control points is used for the approximation of field variables. In order to overcome this problem, h-refinement operations via knot insertion without changing the geometry are implemented. The insertion of a new knotξ ∈ [ξ k , ξ k+1 ] leads to a modification of control points as with whereP stands for the added control points. The existing knot values may be repeated in this algorithm, and the continuity of the basis will be decreased. 16 After that, we obtain a new knot vector Ξ f and a set of particular control points which are used to approximate the field variable, where n f denotes the number of the new control points and is equal to the system's degree of freedom after discretization. By using NURBS basis functions for interpolation, we obtain the sound pressure around the boundary as follows: where p i denotes the sound pressure at the new particular control point which may or may not lie on the boundary, and p f is the order.
In the conventional BEM, the normal practice of collocation at nodal positions is valid. However, this way is no longer valid in the IGA-BEM because the control points may not lie on the boundary. To overcome this difficulty, we use the Greville abscissae definition to define the position of collocation points in parameter space, as follows: Using Eq. (27), we obtain the pressure at the collocation points. By substituting Eq. (27) into Eqs. (6) and (7), we obtain the following boundary integral equation for IGA-BEM and ikβc(r(ξ i )) where ∂p ∂n = ikβp, N e is the number of NURBS element discretizing the boundary curve; [ξ e , ξ e+1 ] denotes a NURBS element and is an interval between two consequent nonrepeating knots in parameter space; J(ξ) stands for the Jacobian. Singular integrals exist in Eqs. (29) and (30) whenξ i ∈ [ξ e , ξ e+1 ], and special treatment needs to be done to eliminate the singular integrals, see Appendix A.
With the FMM applied to accelerate the solution of IGA-BEM, Eq. (19) is rewritten as where the subset S 0 contains N 0 NURBS elements. Similar to the implementation process of the conventional FMBEM, the IGA-FMBEM also needs operation of the M2M, M2L, and L2L. Finally, using Eqs. (23) and (24), we can obtain the solution of the far-field boundary integral based on the NURBS basis functions.

Numerical Examples
In this example, we consider the acoustic scattering of a plane incident wave with a unit amplitude on a rigid cylindrical shell with a radius a = 1 m centered at point (0, 0). The medium in the domain is air with density ρ = 1.2 kg/m 3 . The wave speed is 340 m/sec. The incident wave is traveling along positive x-axis. By comparing the numerical solution with an analytical solution, we can demonstrate the validity and accuracy of the algorithm proposed in this paper. This model can be simplified to be 2D acoustic problem with a circle boundary. The analytical solution of the acoustic pressure at any point located at (r, θ) is given as where ε n denotes the Neumann symbols, i.e. ε 0 = 1; ε n = 2 when n is greater than 0. ( ) stands for the differentiation with respect to ka. Figure 5 shows a comparison between analytical sound pressure and numerical solution based on the IGA-FMBEM, where the test points are located uniformly on a circle with a radius r = 4. Figures 5(a) and 5(b) show the real and the imaginary parts of the sound pressure at test points, respectively. Figures 6 and 7 show the real and the imaginary parts of the sound pressure at a point (4, 0) with frequency, respectively. By observing these figures, we find that the numerical solution obtained using the IGA-FMBEM is in good agreement with the analytical solution, and it confirms the correctness and validity of the algorithm proposed in this paper. Figure 8 shows a comparison of computing accuracy for different interpolation function approximation. The boundary shape is represented by Lagrange linear, quadratic shape function and NURBS curve of 2 degree, respectively. The approximation of physical quantities is implemented by piecewise constant discretization with linear shape approximation (called Lagrange-DBE21), piecewise linear discretization with linear shape approximation (called Lagrange-DBE22), piecewise quadratic discretization with quadratic shape approximation (called Lagrange-DBE33) and NURBS basis function discretization (called NURBS-D2), respectively. Herein, the relative error denotes the surface error based on Euclidean norm, see Ref. 43 for detailed content. A circle with a radius r = 2 m is built. Around 360 computing points distributed on the circle uniformly are used to discretize the circle and compute the surface error. By observing this figure, we can find that: the relative error decreases with increasing DOFs for all different types of the BEM; higher order Lagrange function interpolation produces numerical solution with higher accuracy than lower Lagrange function interpolation; NURBS basis function interpolation performs better than Lagrange function interpolation with same order and DOFs by comparing Lagrange-DBE33 with NURBS-D2. Figure 9 shows a comparison of computing accuracy L. Chen et al.  for different types of BEM with frequency, where the DOFs is set as 5000. The symbol "CBIE" denotes the single boundary integral equation which is used for the solution of the problem, and "BM" denotes the Burton-Miller method which is applied to solve the problem. The result that NURBS interpolation element performs better than Lagrange interpolation element can also be found in Fig. 9. On the other hand, we can find that the peaks originate from the spurious eigen-frequencies. But, Burton-Miller method can be used successfully to overcome this problem by observing the results of NURBS-D2 (CBIE) and NURBS-D2 (BM), where the couple factor is chosen as i/k for k ≥ 1, but i for k < 1.
A key point of this paper consists in applying the FMM into the IGA-BEM for 2D acoustic problem to accelerate the solution of the IGA-BEM. Figure 10 shows that the use of the FMM improves efficiently the computing performance of IGA-BEM, and it demonstrates the high efficiency of the algorithm proposed in this paper. On the other hand, we can find that the use of IGA consumes more computing time by comparing the results of the IGA-FMBEM and the FMBEM. That is because the derived formulas for the NURBS basis functions are recursive and takes much more computing time. Actually, in their current form, the expressions of NURBS basis functions are more expensive than that of the conventional polynomial basis function. The direct solution will reduce the computational efficiency. However, there exist several efficient algorithms for their evaluation, such as the extraction Isogeometric FMBEM for 2D Half-Space Acoustic Problems   operator, 48 and the Cox-de-Boor algorithm. 49 In spite of this, we still develop the applicable of IGA because of its high accuracy and flexibility for shape change. Figure 11 shows the sound pressure contour for scattering field for single and four cylindrical shells problems with f = 1000 Hz, respectively. Figure 12 shows the sound pressure contour for scattering field for octagonal column and eight blade column with f = 1000 Hz, respectively. These figures show the high applicability of the algorithm proposed in this paper for complex problems.   In order to demonstrate the high validity and efficiency of the algorithm proposed in this paper for 2D half-space acoustic problems with absorbing material, an example of noise barrier with curved top structure is tested, see Fig. 13. The data about the coordinate of control points and the associated weight is listed in Table 2. Absorbing material is located on the noise barrier, and the associated normalized surface admittance is set as 1. Figure 14 shows the SPL contour plot for scattering of noise barrier with curved top structure with frequency f = 100 Hz.
Isogeometric FMBEM for 2D Half-Space Acoustic Problems Table 2. Data about control point and the associated weight for noise barrier with curved top structure.

Conclusion
A novel algorithm based on the IGA-FMBEM is presented for the simulation of 2D halfspace acoustic scattering problems with admittance boundary conditions. The Burton and Miller method is used to get correct solutions at all frequencies. A cylinder example for which an analytical solution is available is chosen to demonstrate the correctness and validity of the proposed algorithm, and the performance of different types of boundary elements is presented. The result shows that the proposed method compares very favorably, exhibiting a significantly higher accuracy than the conventional Lagrangian BEM with same degree of freedom. Future work will further extend the proposed algorithm into 2D acoustic optimization analysis and 3D half-space acoustic problems. For 3D acoustic problems, one of the greatest difficulties is to eliminate the singularity of the boundary integrals. we will focus on overcoming this problem.