On the formulation of a BEM in the Bézier–Bernstein space for the solution of Helmholtz equation

On the formulation of a BEM in the B´ezier-Bernstein space for the solution of Helmholtz equation A. Romeroa,∗, P. Galv´ına, J.C. C´amara-Molinaa, A. Tadeub,c aEscuela T´ecnica Superior de Ingenier´ıa, Universidad de Sevilla, Camino de los Descubrimientos, 41092 Sevilla, Spain bITeCons - Institute for Research and Technological Development in Construction, Energy, Environment and Sustainability, Rua Pedro Hispano s/ n., 3030-289 Coimbra, Portugal cADAI - LAETA, Department of Civil Engineering, University of Coimbra, P´olo II, Rua Lu´ıs Reis Santos, 3030-788 Coimbra, Portugal Abstract This paper proposes a novel boundary element approach formulated on the B´ezier-Bernstein basis to yield a geometry-independent field approximation. The proposed method is geometrically based on both computer aid design (CAD) and isogeometric analysis (IGA), but field variables are independently approximated from the geometry. This approach allows the appropriate approximation functions for the geometry and variable field to be chosen. We use the B´ezier-Bernstein form of a polynomial as an approximation basis to represent both geometry and field variables. The solution of the element interpolation problem in the B´ezier-Bernstein space defines generalised Lagrange interpolation functions that are used as element shape functions. The resulting Bernstein-Vandermonde matrix related to the B´ezier-Bernstein interpolation problem is inverted using the Newton-Bernstein algorithm. The applicability of the proposed method is demonstrated solving the Helmholtz equation over an unbounded region in a two-and-a-half dimensional (2.5D) domain. Keywords: subparametric method, B´ezier-Bernstein curve, Newton-Bernstein algorithm, computer-aided design, isogeometric analysis 1. Introduction We present a geometry-independent field approximation of the boundary element method (BEM). The proposed formulation shares some characteristics with the standard collocation methods and isogeometric analysis (IGA), exploting the advantages of each method. A similar concept was recently introduced by Atroshchenko and Bordas [1]. The two strategies are briefly introduced to put this approach into context. The standard collocation methods are commonly used in engineering to solve boundary value problems [2]. Typically, the objective is to find a function which satisfies a given differential equation with prescribed boundary conditions. The procedure for reducing the continuous problem to its discrete formulation is the ∗Corresponding author Email address: [email protected] (A. Romero) Preprint submitted to Applied Mathematical Modelling May 10, 2019 discretisation of both domain geometry and field variables. The geometry and problem solution are then both approximated as a linear combination of polynomial basis functions defined on a given mesh. The solution accuracy is conditioned by the domain discretisation. Mesh refinement and polynomial based order increase, hp−refinement, are strategies used to improve the problem solution. Mesh refinement is a key process that starts with a CAD model of the domain. Each successive mesh is refined in the CAD model. The analysis procedure could fail when the domain geometry needs to be accurately represented. Nowadays, the number of papers dedicated to the development of the BEM and other mesh reduction methods covers a wide field of applications [3–5]. These advances include the development boundary elements for solving thin-walled structures [6] or meshless based methods [7, 8], among others. In the past decade isogeometric analysis (IGA) has been proposed to deal with geometry approximation. IGA has proved to be very attractive in solving a wide range of problems because it can represent geometries exactly [9–13]. IGA was first introduced by Hughes et al. [14] “to be geometrically exact no matter how coarse the discretisation”. Also, it simplifies successive mesh refinement since CAD communication is not necessary once a coarse mesh has been defined. It is thus a flexible method that dispenses the polynomial geometry approximation. The major findings reported in [14] are summarised as follows: i) IGA formulation is based on non-uniform rational B-splines basis (NURBS); ii) a mesh for a NURBS patch is defined by the knot vectors and the span between knots divides the representation domain into elements; iii) the geometry is defined as a linear combination of B-spline basis with control points or de Boor points; iv) field variables are represented by the same function as the geometry by control variables; v) mesh refinement is based on knot insertion and order elevation techniques without CAD model interaction; vi) element matrices are integrated for each NURBS patch and they are assembled into a global matrix. One of the main drawbacks of IGA is that the geometry is represented in terms of control points that do not necessarily belong to the boundary, so related control variables do not have a clear physical meaning. Moreover, defining boundary conditions is not straightforward because control variables do not lie on the boundary. The application of inhomogeneous Dirichlet boundary conditions might affect the solution accuracy [14, 15]. Some authors have proposed alternatives based on constraining equations, Lagrange multipliers, penalty methods, or some other method. This fact complicates the isogeometric formulation of the BEM, however, since nodal collocation points do not belong to the boundary [16–19]. Alternatively, the Greville abscissae are used to define the collocation points [20]. The proposed strategy can be understood as a formulation that lies between the standard collocation method and the isogeometric analysis. The main thrust of this formulation is the use of an independent basis to represent the geometry and the field variables [1]. This approach was successfully used before the isogeometric development. The works of Zieniuk [21], Zieniuk and Boltuc [22], and Zieniuk and Szersze´n [23], among others, proposed the concept of parametric integral equation system (PIES) for accurate representation of boundary geometry. In these works, the boundary geometry was described as a parametric 2 function that was included in the PIES kernels, rather than in the boundary integral as is done in standard BEM. Then, the solution took into account the boundary geometry defined by a parametric function as it is usually done in CAD. Moreover, the authors overcame the difficulties of the original IGA related to the isoparametric concept, computing the solution of PIES by an independent polynomial basis for the field variable approximation [22]. In this work, we started with the idea that a B-spline can be converted to B´ezier curve using the Boehm algorithm [24]: a B-spline of order nis transformed to a composite B´ezier curve by inserting knots until all interior knots have a multiplicity equal to n−1. Therefore, taking this identity into account, the proposed method will be derived in the B´ezier-Bernstein space hereinafter. The geometry is represented by B´ezier curves defined as a linear combination of Bernstein basis with control points. The Bernstein basis is widely used in CAD because of its good properties for building and modifying geometric shapes [25]. On the other hand, the use of orthogonal basis allows simple and efficient field approximation. To define an orthogonal basis for the field interpolation, therefore, a least-square problem in terms of Bernstein-Vandermonde matrix is generated by the interpolation constraints whose solution is a Lagrange interpolant relative to the Bernstein basis [26, 27]. This methodology is an efficient and robust approach that combines the Bernstein basis for geometry representation and the Lagrange interpolating polynomial relative to the Bernstein basis for the field approximation. The proposed method is powerful because it includes the main characteristics of standard BEM formulation and isogeometric analysis: i) the geometry is exactly defined as B´ezier curves; ii) field variables are represented at nodal points instead of control variables, which enables easy definition of the boundary conditions; iii)h−refinement is performed directly on the geometry without CAD interaction, as it is performed by IGA; and iv)p−refinement only affects the field approximation, rather than the geometry discretisation. Table 1 summarises the characteristics of the BEM approach, isogeometric methods and the proposed formulation of the BEM. In this paper, the capability of the proposed method is demonstrated for solving the Helmholtz equation in longitudinally invariant problems. This situation is referred to as a two-and-a-half dimensional (2.5D) problem [28], in which the three-dimensional (3D) solution is computed by means of a Fourier transform. The structure of the rest of the paper is as follows. First, the B´ezier-Bernstein space is briefly described and the Lagrange interpolant relative to the Bernstein basis is derived. Next, the boundary element formulation in the B´ezier-Bernstein space is presented, after which, the proposed methodology is verified by means of two benchmark problems and a numerical example is presented. 3 Table 1: Comparison of the boundary element method, based on standard collocation method, isogeometric analysis and geometry-independent field approximation analysis. Collocation Isogeometric Proposed method Classification Sub-/iso-/superparametric Isoparametric Sub-/iso-/superparametric Geometry Approximate Exact Exact Nodal points Control points Control points Polynomial basis B-splines B´ezier curves Field Nodal variables Control variables Nodal variables Interpolation basis Approximation basis Interpolation basis Boundary condition/ Nodes Control points/ Nodes Collocation point Greville abscissae 2. Numerical model This approach is based on the Helmholtz integral equation for solving 3D interior, exterior and scatterering problems in the frequency-wavenumber domain. The pressure field for the acoustic waveguide at point x=x(x, y, z) is computed as the superposition of 2D solutions with different wavenumbers κz[29]: p(x, ω) = Z+∞ −∞ ep(e x, κz, ω)e−ικzzdκz(1) where ep(e x, κz, ω) is the frequency-wavenumber representation of the sound pressure, e x=x(x, y) and ι= √−1. The boundary integral representation is derived from the Helmholtz equation. The fundamental solution e Ψ(e x, κz, ω;e xi) for sound pressure at receiver position e xdue to a source acting at e xiis the solution to: ∇2e Ψ(e x, κz, ω;e xi) + κ2e Ψ(e x, κz, ω;e xi) = −δ(r) (2) where κ=ω/cfis the problem wavenumber, cfis the sound propagation velocity, δis the Dirac delta function and r=||e xi−e x|| is the distance from the source to the receiver. The solution to this equation defines the fundamental solution for an unbounded region in the frequency-wavenumber: e Ψ(e x, κz, ω;e xi) = −ι 4H(2) 0(καr) (3) where κα=pκ2−κ2 zand H(2) 0is the Hankel function. The following sections deal with: i) the B´ezier-Bernstein space for geometry representation, ii) the Lagrange interpolant relative to Bernstein basis for variable field approximation, and iii) the boundary element formulation in the B´ezier-Bernstein space. 4 2.1. The B´ezier-Bernstein space The Bernstein polynomial basis was introduced to prove the Weierstrass theorem: given any continuous function f(x) on an interval [a, b] and a tolerance  > 0, there is a polynomial Pn(x) of sufficiently degree n, such that [25]: |f(x)−Pn(x)| ≤ , x ∈[a, b] (4) The Bernstein basis of degree nis defined over the interval t∈[0,1] as: Bn k(t) = n ktk(1 −t)n−k, k = 0, . . . , n (5) The change of variable to map x∈[a, b] to t∈[0,1] allows the definition of the Bernstein polynomial related to a continuous function f(t) as follows: Pn(t) = n X k=0 fk nBn k(t) (6) Then, Equation (6) converges to f(t) as ngoes to infinity due to the Bernstein basis being non-negative and forming a partition unity, although the polynomial form does not interpolate the sample values f(k/n). The basis defined by Equation (5) is motivated by its capacity to represent a polynomial in Bernstein form [26] through coefficients bk: Pn(t) = n X k=0 bkBn k(t) (7) The Bernstein basis has partition unity property, is defined as non-negative and symmetric about the interval mid-point (Bn n−k(1−t) = Bk n(t)). Also, the Bernstein form satisfies the following recursive definition: Bn k(t) = (1 −t)Bn−1 k(t) + tBn−1 k−1(t) (8) The Bernstein basis derivatives are defined by: d dtBn k(t) = nBn−1 k−1(t)−Bn−1 k(t)(9) where Bn−1 −1(t) = 0 and Bn n−1(t) = 0. Thus, the derivative of a polynomial in Bernstein form is obtained by combining Equations (7) and (9): d dtPn(t) = n−1 X k=0 n[bk+1 −bk]Bn−1 k(t) (10) 5 The application of polynomials in Bernstein form growth with the development of B´ezier curves rn(t) in computer-aided design: rn(t) = n X k=0 bkBn k(t) (11) where bkare the control points used to approximate the geometry. A control polygon is obtained by connecting the control points and this is used to modify the curve geometry. The de Casteljau algorithm is often used for evaluating and splitting a B´ezier curve rn(t) at a given point t[25]. The recursive property of the Bernstein basis defines successive control points bl k(t) as a recursive linear interpolation: bl k(t) = (1 −t)bk(t)l−1+tbk+1(t)l−1, l = 1, . . . , n , k = 0,1, . . . , n −l(12) Starting with the control polygon given by control points b0 k=bk, de Casteljau algorithm defines a new polygon whose vertex is smaller than the original one. The last point in the recursive procedure corresponds to the B´ezier representation r(t) = b0 n(t). Moreover, two new B´ezier curves can be defined by control points {bl k(t)}l=0,...,n k=0 and {bl k(t)}l=n,...,0 k=n−l. Although the de Casteljau algorithm allows an easy evaluation of a B´ezier curve, it is computationally expensive. An efficient curve computation is achieved using the polar form (or blossom) of a B´ezier curve rn(t) [30], which defines a multi-affine transformation satisfying: bk=R(0,...,0 | {z } n−k ,1,...,1 | {z } k ) (13) where R(t1, . . . , tn) is computed as: R(t1, . . . , tn) = X I∩J=∅ I∪J={1,2,...,n}Y i∈I (1 −ti)Y j∈J tjb|J|(14) Thus, a polynomial in Bernstein form can be formulated in polar form substituting Equation (13) into Equation (11) as follows: rn(t) = n X k=0 R(0,...,0 | {z } n−k ,1,...,1 | {z } k )Bn k(t) = R(t, . . . , t) (15) The polar form of a B´ezier curve rn(t) of degree ndefines a multi-affine polynomial R(t1, . . . , tn) with n variables that satisfy rn(t) = R(t, . . . , t). Moreover, the following expression can be derived from Equations 6 (12) and (13): R(0,...,0 | {z } n−(l+k) , t, . . . , t | {z } l ,1,...,1 | {z } k ) = (1 −t)R( 0,...,0 | {z } n−(l+k)+1 , t, . . . , t | {z } l−1 ,1,...,1 | {z } k ) +tR(0,...,0 | {z } n−(l+k) , t, . . . , t | {z } l−1 ,1,...,1 | {z } k+1 ), l = 1, . . . , n , k = 0,1, . . . , n −l (16) which it is in agreement with Equation (14) when l=n. The new B´ezier curves defined by control points {bl k(t∗)}l=0,...,n k=0 and {bl k(t∗)}l=n,...,0 k=n−lat a given point t∗ according to the de Casteljau algorithm are represented by: r∗ n(t) = n X k=0 R(0,...,0 | {z } n−k , t∗, . . . , t∗ | {z } k )Bn k(t) (17) and, r∗ n(t) = n X k=0 R(t∗, . . . , t∗ | {z } n−k ,1,...,1 | {z } k )Bn k(t) (18) respectively. The BEM formulated within the B´ezier-Bernstein space needs the derivative of rn(t) to compute tangent and normal vectors and the Jacobian of the transformation. The first derivative is evaluated according to Equations (10) and (11) as: d dtrn(t) = n−1 X k=0 n[bk+1 −bk]Bn−1 k(t) (19) This expression can be rewritten using the polar form as follows [31]: d dtrn(t) = d dtR(t, . . . , t) = nR(t, . . . , t | {z } n−1 ,1) −R(0, t, . . . , t | {z } n−1 )(20) Once the first derivative of rn(t) is computed, the tangent and normal vectors and the Jacobian are easily obtained. Cubic B´ezier approximations are widely used to represent smooth patches in CAD. More complex shapes would be defined by either increasing the curve order or using a piecewise B´ezier curve. In this case, control points must satisfy C1and C2end-point interpolation conditions [32, 33]. Figure 1 shows a third-order 7 B´ezier patch and its related control points represented in polar form as: R(t1, t2, t3) = [(1 −t1)(1 −t2)(1 −t3)] b0 + [(1 −t2)(1 −t3)t1+ (1 −t1)(1 −t3)t2+ (1 −t1)(1 −t2)t3]b1 + [(1 −t3)t1t2+ (1 −t1)t2t3+ (1 −t2)t3t1]b2 + [t1t2t3]b3 (21) The resulting new curves are defined through control points {R(0,0,0),R(0,0, t),R(0, t, t),R(t, t, t)}and {R(t, t, t),R(t, t, 1),R(t, 1,1),R(1,1,1)}according to the de Casteljau algorithm. R(0,0,0) R(0,0,1) R(0,1,1) R(1,1,1) R(0,0, t) R(0, t, 1) R(t, 1,1) R(0, t, t) R(t, t, 1) R(t, t, t) Figure 1: Third-order B´ezier patch. The geometry approximation presented in this section is used to formulate the BEM in the B´ezier- Bernstein space. The boundary geometry is defined by patches with an arbitrary approximation degree, that depends on the geometry shape complexity. Then, each patch is subdivided into elements according to h−refinement. Section 3.2 analyses the accuracy of the proposed method studying the scattered wavefield by a circular inclusion using C1third and C2fourth order B´ezier curves. 2.2. Field approximation The B´ezier-Bernstein space allows the approximation of a curve geometry through control points. Moreover, the element shape function can be derived from the Lagrange interpolant relative to the Bernstein basis presented in the previous section. The field approximation consists of variable interpolation from known nodal values through element shape functions. The approximation of a field variable is given by finding a polynomial Pn(x) that interpolates a function f(x) at discrete points xj: Pn(xi) = f(xi), i = 0, . . . , n (22) where xiare n+ 1 nodes and f(xi) are given data. The problem has a solution that only depends on the interpolation data. The solution can be defined in Lagrange form [34], but each evaluation of Pn(x) 8 results in an O(n2) complexity algorithm and its computation is numerically unstable [35]. Berrut and Trefethen [34] presented a barycentric Lagrange formula which has O(n) complexity. Moreover, they stated that setting equidistant nodes for high-order interpolation polynomials causes the appearence of the socalled Runge phenomenon, which results in a highly ill-conditioned basis. The polynomial basis becomes a well-conditioned process for polynomial interpolation when elements with nodes mainly concentrated at the endpoints with an asymptotic density proportional to (1 −x)−1/2are used. Legendre-Gauss-Lobatto (LGL) integration points are used in spectral formulations for this purpose [36]. Alternatively, we propose simpler node distributions such as the families of Chebyshev points to develop the BEM in the B´ezier-Bernstein space. Chebyshev points of the first kind are given by: xi= cos (2i+ 1)π 2n+ 2 , i = 0, . . . , n (23) and, of the second kind by: xi= cos iπ n, i = 0, . . . , n (24) In these cases, if f(x) is a smooth function defined at the interval [0,1] interpolated by polynomials in Chebyshev points, the rate of convergence is remarkably faster. The interpolant Pnsatisfies an error estimation, thus: max x∈[0,1] |f(x)−Pn(x)| ≤ CK−n(25) for constants Cand K > 1 [37]. The Lagrange interpolant relative to the Bernstein basis has become useful for high-order element development [38]. Although a polynomial in Bernstein form satisfies Pn(0) = f(0) and Pn(1) = f(1), these properties do not remain valid at internal element points as has been stated. Therefore, it is necessary to properly modify the control points bkto ensure the following condition: Pn(xi) = n X k=0 bkBn k(xi) = f(xi) (26) This condition is commonly expressed as a linear system of equations through the Bernstein-Vandermonde matrix Aij =Bj i(xi), for i, j = 0, . . . , n. The control points b= [b0, . . . , bn]Tare given by the solution of Ab =f, where the vector f= [f(x0), . . . , f(xn)]Tcollects the interpolation data. Some algorithms have been presented for solving Bernstein-Vandermonde linear systems [35, 39–41]. Marco and Mart´ınez [39] studied the condition number of linear systems defined by the Bernstein-Vandermonde matrix and concluded that they are exponentially ill-conditioned. In this work, we have used the Newton- Bernstein algorithm proposed by Ainsworth and S´anchez [35] to compute the control points. 9 normals at this location. However, the compatibility of normal velocities is not achieved when assembling matrix e G(e x, κz, ω;e xi). The assembled matrix has one extra column for each discontinuity point at the boundary and the particle normal velocity vector also has an extra component. Another possibility commonly used for the treatment of normal discontinuities in potential problems consists of shifting the endpoint nodes inside the element. Thus, the computed solution is discontinuous between elements. Parreira [47] analysed the approximation error for different nodal distribution in discontinuous boundary elements and concluded that the optimal position for nodes inside the elements coincided with Gauss-Legendre integration points. In the next section we show that Chebyshev points of the first kind can be used to define discontinuous elements which simplify the treatment of particle normal velocity at boundary discontinuities. Furthermore, discontinuous elements allow an easy evaluation of the free term beacause the element geometry is smooth at the collocation points, to give ci(e xi)=1/2. Therefore, boundary elements with nodes located at Chebyshev point of the first kind provide a reliable formulation without loss of accuracy as will be seen in next section. Finally, to conclude this section the proposed hp−refinement is described. The solution refinement is quite similar to a standard BEM approach. Boundary elements are subsequently divided using the de Casteljau algorithm and preserving the exact geometry in each refinement without CAD interaction. Moreover, the element order is increased up without comprising the B´ezier-Bernstein basis used to describe the boundary geometry. Table 2 summarises the proposed method. Table 2: Summary of the boundary element method formulated in the B´ezier-Bernstein space AThe B´ezier-Bernstein space settings: A.1 Import geometry from CAD model. A.2 Path subdivision for h−refinement through de Casteljau algorithm. BField approximation: B.1 Set element nodes distribution and element order p(p−refinement). B.2 Compute element shape function from the Lagrangian interpolant relative to the Bernstein basis. CBoundary element integration: C.1 Set the critical radius r∗for element integration. C.2 For each collocation point, integration elements further than r∗are regularly integrated. C.3 If any integration element is closer than r∗, then a smoothing transformation is used for singular integration. DSolution: D.1 Boundary solution is obtained by imposing appropriate boundary conditions. D.2 Domain solution is computed following the same boundary element integration scheme. 16 3. Numerical verification In this section, we analyse the performance of the proposed method with a benchmark problem concerning the wave propagation in an unbounded domain in two different configurations: i) open boundary, and ii) rigid boundary. An incident pressure field given by epI(e x, ω) = exp(ικd·e x) was considered in both cases, where the polarised direction was d= [1,0]. This problem has an analytical solution for the two-dimensional Helmholtz equation that is equivalent to setting the longitudinal wave number at κz= 0. The remaining properties were fluid density ρ= 1.225 kg/m3and sound propagation velocity cf= 340 m/s. Numerical results were compared with reference solutions and the L2scaled error, 2, was used to assess the accuracy [48]: 2=kepex(e x, κz, ω)−ep(e x, κz, ω)k kepex(e x, κz, ω)k(51) where epex denotes the analytical solution and epis the numerical result computed by the proposed methodology. Different hp−strategies were investigated to get the optimal discretisation with the lowest computational effort. The nodal density per wavelength is used to describe the mesh density [49]: dλ=2πp κh (52) where κis the characteristic wavenumber, his the element size and pis the order of the shape functions. A convergence analysis was carried for several element sizes hwith successive p−enrichment. Boundary elements were subsequently divided to perform a h−refinement and the element order was increased until convergence is reached. For this purpose, we considered that the problem solution was properly approximated if errors satisfied log(2(h, p −1)/2(h, p + 1)) ≤1. 3.1. Case 1: open boundary In this example, the Helmholtz equation was solved over a square boundary with sides of 2 m at a wavenumber κ= 7.5 rad/m (see Figure 5). Four linear B´ezier patches were used to define the boundary geometry. The open boundary was set with a Dirichlet condition at the left and right sides and a Neumann condition at the other sides. The problem solution was computed over a grid of receivers inside the boundary allowing a density of twelve points per wavelength. The accuracy of the computed results was assessed using different families of shape functions defined at: i) equidistant nodes, ii) LGL integration points, iii) Chebyshev points of the first kind, and iv) Chebyshev points of the second kind. Boundary patches were discretised for different element sizes defined by κh = 7.5, κh = 3 and κh = 1. Figure 6 shows the errors in computed results plotted versus the element order and the nodal density per wavelength. The lowest error was O(10−10) for interpolation functions defined at Chebyshev points of the second kind for the finest discretisation. However, the solution computed by equidistant nodes had the highest error for a given element size and element order due to deterioration of the numerical results 17 -2 -1 0 1 x[m] -2 -1.5 -1 -0.5 0 0.5 1 1.5 y[m] -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 Figure 5: Wave propagation over an unbounded domain. The open boundary is represented in red line. induced by the Runge phenomenon. The computed results were quite similar for the other interpolant families. This analysis showed that the boundary element formulation based on discontinuous elements defined at Chebyshev points of the first kind is suitable for the easy treatment of normal discontinuities without loss of accuracy. The convergence was faster for finer discretisation than coarser discretisation. Numerical results converged to the analytical solution when the nodal densities per wavelength were higher than dλ= 3 for the coarsest mesh, dλ= 4 for the medium element size, and dλ= 12 for the smallest element size. Moreover, the solution error for a fixed nodal density was lower for a coarse mesh than a fine mesh, as can be seen in Figure 6.(b). Therefore, a coarse mesh with an adequate nodal density gave better results than a fine mesh with the same density, or equally, the use of high-order elements in coarse meshes (p−refinement) provided better results than low-order elements in a fine discretisation (h−refinement) for a fixed nodal density. 3.2. Case 2: rigid boundary Next, the accuracy in representing boundary geometries in the B´ezier-Bernstein space was evaluated by studying the wave scatterering by a fixed cylindrical cavity with radius r= 1 m. The accuracy of the computed results is constrained by the level of precision used to approximate the geometry [50]. In this case, the circular geometry was approximated either by a C1cubic or by a C2quartic B´ezier curves [33]. Control points for accurate representation of circular arcs can be found in [33]. The performance of the proposed method was compared with the isogeometric approach presented by Peake et al. [51]. The analytical solution of this problem can also be found in [33]. First, the problem solution was computed for a wavenumber κ= 40 rad/m, κh = 3 and nodal density dλ= 10. The geometry was defined by four cubic patches subdivided into 21 boundary elements of order p= 5. The element nodes were located at the Chebyshev points of the second kind. The problem solution 18 0 2 4 6 8 10 12 14 16 18 20 Element order p 10−10 10−8 10−6 10−4 10−2 100 Scaled L2error ǫ2 Equidistant nodes LGL integration points Chebyschev 1st kind Chebyschev 2nd kind (a) 100101102 Nodal density dλ 10−10 10−8 10−6 10−4 10−2 100 Scaled L2error ǫ2 Equidistant nodes LGL integration points Chebyschev 1st kind Chebyschev 2nd kind (b) Figure 6: Convergence of sound pressure using different families of shape functions for κh = 7.5 (dotted lines), κh = 3 (dashed lines) and κh = 1 (solid lines) plotted versus (a) element order pand (b) nodal density per wavelength dλ. shows a shadow region behind the cavity, as can be seen in Figure 7.(a). The pressure amplitude over the cavity is represented in Figure 7.(b). Maximum amplitude was reached when the incident wave struck the cavity at point θ=πrad, whereas null pressure was found at the opposite side (θ= 0 rad). The computed results were in very good agreement with the analytical solution [51]. (a) 0 30 60 90 120 150 180 210 240 270 300 330 0.5 1 1.5 Analytical solution Numerical results (b) Figure 7: (a) Real part of wavefield scatterered by a fixed circular cavity at κ= 40 rad/m. (b) Pressure amplitude at the boundary. Thereafter, the solution error over the wavenumber range κ∈(1,400) rad/m was compared with those obtained by Peake et al. [51] using an isogeometric boundary formulation (IGABEM) and a conventional BEM. The procedure followed by these authors consisted of adding new degrees of freedom to the problem discretisation as the wavenumber increased, to ensure a fixed nodal density dλ= 10. We also followed this 19 scheme, starting with three patches at the initial wavenumber κ= 1 rad/m. Each patch was subsequently subdivided into elements according to κh = 3 as the wavenumber increases. The element order was fixed to achieve a nodal density dλ= 10. As an example, Figure 8 represents the related boundary discretisation for κ= 1 rad/m and κ= 2.7 rad/m using C1cubic and C2quartic B´ezier curves. -1.5 -1 -0.5 0 0.5 1 1.5 x[m] -1.5 -1 -0.5 0 0.5 1 1.5 y[m] κ= 1 rad/m (a) -1.5 -1 -0.5 0 0.5 1 1.5 x[m] -1.5 -1 -0.5 0 0.5 1 1.5 y[m] κ= 2.7 rad/m (b) -1.5 -1 -0.5 0 0.5 1 1.5 x[m] -1.5 -1 -0.5 0 0.5 1 1.5 y[m] κ= 1 rad/m (c) -1.5 -1 -0.5 0 0.5 1 1.5 x[m] -1.5 -1 -0.5 0 0.5 1 1.5 y[m] κ= 2.7 rad/m (d) Figure 8: Circular approximation (black line) with (a,b) C1cubic B´ezier curves and (c,d) C2quartic B´ezier curves. Control points and their related control polygons are represented by grey circles and grey lines. Figure 6 shows the scaled L2error 2of the pressure solution at the boundary. The analysis of the solution error obtained from cubic B´ezier curves (n= 3) showed a decreasing tendency at low wavenumber until the element size was small enough to accurately represent the circular geometry. Afterwards, the error became uniform of O(10−5). The previously mentioned initial effect was avoided when C2quartic B´ezier curves (n= 4) were used because they can better represent circular arcs. The problem solution deteriorated slightly for higher wavenumbers in both cases, although good results were nonetheless achieved. In general, the proposed method gave more accurate results than the isogeometric formulation and the conventional 20 BEM over the studied wavenumber range. However, the solution errors obtained by these methods at high wavenumbers were more uniform than those obtained with the proposed method. Figure 9: Scaled L2error for different wavenumber compared to IGABEM and BEM [51] 3.3. A remark on the condition number Finally, in this section, we study the condition number for the system of equations in the previous examples. The condition number of BEM matrix is assumed to be moderate for practical applications. It is well known that it increases with the wavenumber kin the solution of the Helmholtz equation due to the pollution effect. Nevertheless, the BEM does not seem to be susceptible to pollution as referred in Reference [52]. L¨ohndorf and Melenk [52] proved that the hp-BEM accuracy is independent of the wavenumber if the condition kh/p is sufficiently small and the element order is p≥ O(log k). Although this condition makes it possible to obtain well-posed BEM matrices, in practical applications can produce large systems as the wavenumber increases. Figures 10 and 11 show the condition number in 2-norm for cases 1 and 2, respectively. The analysis in the first case showed that the condition number remained constant regardless of the element order for a fixed wavenumber. This is true for all families of elements studied except for equidistant nodes. The condition number rapidly increased for element order higher than six due to the Runge phenomenon if the equidistant family of shape function was used. In the second case, the wavenumber varied in the range 1 −400 rad/m. The analysis was carried out with a constant nodal density per wavelength dλ= 10 regardless the wavenumber. This condition leads to a relation kh/p = 0.6 with a fixed element order p= 5 (kh = 3). Figure 11 shows that the condition wavenumber increased with kand it did not depend on the basis order used for the geometry. This analysis allows us to conclude that the proposed method is at least as well-posed as other formulations to accurately solve this problem, even though the condition given by L¨ohndorf and Melenk [52] was not achieved due to 21 limitations in the system size. This analysis reveals that the use of independent polynomial basis produces well conditioned systems independently of the approximation developed for geometry and field variable. Figure 10: Condition number using different families of shape functions for κh = 7.5 (dotted lines), κh = 3 (dashed lines) and κh = 1 (solid lines) plotted versus element order p. Figure 11: Condition number for different wavenumber. 4. Numerical example In this section, the capability of the proposed method is shown with a numerical example. The problem solved herein studies the three-dimensional wave propagation in a fluid channel with a fixed scatterer. The fluid channel is bounded by two reflection planes that were assumed to be rigid with null particle normal velocity. The total channel height was H= 6 m and fluid air properties were taken to be those described in the previous section. 22 A modified Green’s function for the fluid channel was used to limit the discretisation to the scatterer boundary. The method of images was used to construct the modified Green’s function [44]. The scatterer and the fluid channel were centred at the origin. The fluid channel and the scatterer were excited by a point source acting at location x= (−3,0,0) m with κ= 20 rad/m. The source was defined by the modified Green’s function. The problem solution was computed in the frequency-wavenumber domain for a longitudinal wavenumber varying over a range from zero to κz= 2κin 128 equidistant steps. The pressure field was evaluated over a grid of receivers with a density of six points per wavelength. Two geometries were considered: i) a cylindrical scatterer as shown in Figure 8 and ii) the complex shaped scatterer represented in Figure 12.(a) directly imported from a CAD model. The geometry of the second example was defined by twelve C1cubic patches. Problem discretisation was given by κh = 6 and dλ= 6, which leads to an element order p= 6. The element nodes were located at the Chebyshev points of the second kind. Once the boundary discretisation had been defined, the geometry patches were subdivided as shown in Figure 12.(b). The cylindrical scatterer was discretised following a similar procedure. -1.5 -1 -0.5 0 0.5 1 1.5 x[m] -1.5 -1 -0.5 0 0.5 1 1.5 y[m] (a) -1.5 -1 -0.5 0 0.5 1 1.5 x[m] -1.5 -1 -0.5 0 0.5 1 1.5 y[m] (b) Figure 12: (a) Complex geometry imported from CAD model and (b) related discretisation for κh = 6. Control points and their related control polygons are represented by grey circles and grey lines. Figure 13 shows the incident wavefield and the scatterered wave. The representation is limited to a half model, according to the problem symmetry. The incident pressure field shows a spherical distribution that was gradually lost as the distance from the source point increased due to channel wave refraction (Figure 13.(a)). The amplitude decreased with the distance to the excitation point, according to the radiation condition. The pattern of the incident wave changed under the effect of the scatterers. A complex pattern was found, caused by waves reflected from both the channel borders and the scatterer boundary. A shadow region was found behind the scatterers where pressure amplitude was considerably lower. This example has shown that the proposed method is, therefore, a useful tool for studying boundary 23 problems with exact CAD geometry. An independent geometry and variable field approximation allowed efficient modelling that took both the complexity of the boundary shape and the studied wavenumber into account. In this case, it was enough to use cubic curves for the geometry approximation, while it was necessary to use sixth order elements for the variable field. 24 (a) (b) (c) Figure 13: Real part of (a) incident wavefield, and the scatterered wavefield by (b) a cylindrical cavity, and (c) a complex shape boundary. 25