An improved boundary point interpolation method for exterior acoustic radiation problem

Engineering Analysis with Boundary Elements 103 (2019) 11–21

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

An improved boundary point interpolation method for exterior acoustic radiation problem Linchong Chen a, Xiaolin Li a,b,∗ a b

School of Mathematical Sciences, Chongqing Normal University, Chongqing 400047, PR China Key Laboratory for Optimization and Control Ministry of Education, Chongqing Normal University, Chongqing 400047, PR China

a r t i c l e

i n f o

Keywords: Meshless Boundary point interpolation method Acoustic radiation Exterior problem Hypersingular integral Nearly singular integral

a b s t r a c t In this paper, an improved boundary point interpolation method (BPIM) is developed for solving exterior acoustic radiation problems. With the help of Lagrange basis polynomials, a Lagrange point interpolation method (PIM) is presented to avoid the inversion of the moment matrix in the PIM. Then combining the acoustic double layer potential with the Lagrange PIM, the improved BPIM for exterior acoustic radiation problems is presented. Detailed and efficient computational formulas are established to compute acoustic hypersingular and nearly hypersingular integrals. Numerical results show that the current BPIM can produce efficient numerical solutions to exterior acoustic radiation problems.

1. Introduction The exterior acoustic radiation problem has many applications in science and engineering to describe radiation patterns and wave propagations of vibrating surfaces. A number of computational methods such as the finite element method (FEM) [1] and the boundary element method (BEM) [2–4] have been widely used to simulate acoustic problems. To alleviate meshing-related drawbacks in the FEM and the BEM, the radial point interpolation method [5], the meshless local Petrov–Galerkin (MLPG) method [6], the smoothed FEM [7,8] and the element-free Galerkin (EFG) method [9] have been developed to simulate acoustic problems. For exterior acoustic problems, the infinite domain and the Sommerfeld radiation condition need special attention and treatment in these meshless methods. Because boundary integral equations (BIEs) require only boundary discretization and satisfy the Sommerfeld radiation condition naturally, some BIEs-based meshless methods, such as the singular meshless method [10], the boundary knot method [11,12], the method of fundamental solutions [13–15], the singular boundary method [16–20] and the boundary element-free method (BEFM) [21,22] have been developed to simulate acoustic problems. These meshless methods require only nodes on the boundary [23,24]. In many meshless methods, such as the MLPG, the EFG and the BEFM, the moving least squares (MLS) approximation [25] is used to form meshless shape functions. The lack of interpolation properties of the MLS shape functions causes inconveniences of implementing boundary conditions in these methods [26–30]. The point interpolation

∗

method (PIM) [26,31] is another efficient and frequently used technique to form meshless shape functions. Unlike the MLS approximation, the PIM shape functions fulfill the desired interpolation properties. Many meshless methods have been developed by using the PIM [26]. Based on radial basis functions, the radial PIM has also been widely used to form meshless shape functions [32–35]. The boundary point interpolation method (BPIM) [26,36] is a meshless method that combines BIEs and the PIM. The BPIM has been used to solve elasticity problem [26,36], biharmonic-type problem [37], Signorini problem [38], seismic soil-tunnel interaction problems [39], etc. Excellent numerical results are obtained for these problems. In the BPIM, the PIM is used to interpolate unknown boundary functions in BIEs. In the PIM, basis functions are monomials, and an algebra equations system must be solved for every computation point. Then, the computation of shape functions involves the inversion of the moment matrix, which increases CPU time. In this paper, by using Lagrange basis polynomials, a Lagrange PIM is presented to ameliorate this drawback. Shape functions in the Lagrange PIM are explicitly expressed and can be directly computed with no inverse matrix. Thus, shape functions in the Lagrange PIM are much simpler and require less CPU time than those in the polynomial PIM. With the aid of the Lagrange PIM, an improved BPIM is developed in this paper for boundary-only analysis of two-dimensional exterior acoustic radiation problems. In the improved BPIM, the exterior acoustic radiation problem is reduced to a hypersingular BIE by using the acoustic double layer potential, while the Lagrange PIM is used to interpolate the unknown function in the BIE. Due to the complexity of the

Corresponding author at: School of Mathematical Sciences, Chongqing Normal University, Chongqing 400047, PR China. E-mail address: [email protected] (X. Li).

https://doi.org/10.1016/j.enganabound.2019.02.002 Received 9 December 2018; Received in revised form 24 January 2019; Accepted 15 February 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.

L. Chen and X. Li

Engineering Analysis with Boundary Elements 103 (2019) 11–21

fundamental solution of two-dimensional acoustic problem, boundary integrals in the discretized hypersingular BIE are analyzed in detail. It can be found that there have logarithmic weakly singular integrals, hypersingular integrals and nearly hypersingular integrals. The success of BIEs-based methods depends mainly on the efficient computation of these singular integrals. In the BEM, a number of techniques have been developed to compute singular integrals [40–45]. In this paper, detailed and efficient computational formulas are established to compute acoustic hypersingular and nearly hypersingular integrals in the current BPIM. Especially, a power series expansion method is developed to compute hypersingular integrals, while a cell subdivision method is developed to compute the nearly hypersingular integrals. Numerical results show that the current BPIM can produce efficient numerical solutions to exterior acoustic radiation problems.

2.2. Lagrange PIM

2. Point interpolation method (PIM)

𝓁𝑗 (𝑥(𝑠)) =

In the polynomial PIM, Eq. (4) shows that the shape function involves the inversion of the moment matrix Pm . Clearly, the computation of the inverse matrix will consume some computing times. To ameliorate this drawback, a Lagrange PIM is presented in this section. Given m distinct data points (s1 , 𝜇 1 ), (s2 , 𝜇 2 ), … , (sm , 𝜇 m ) surrounding the point x(s) ∈ Γ, according to the classical Lagrange interpolation [47], the interpolation polynomial in the Lagrange form for the function 𝜇(x) is 𝜇(𝑥(𝑠)) =

𝑗=1

𝑗=1

𝑎𝑗 𝑝𝑗 (𝑠) = 𝐩T (𝑠)𝐚,

𝑥 (𝑠 ) ∈ Γ

(1)

𝜆 = max 𝑥∈Γ

Φ𝑗 (𝑥(𝑠))𝜇𝑗 ,

𝑥 (𝑠 ) ∈ Γ

where Φj (x(s)) is the PIM shape function, [ ] Φ𝑗 (𝑥(𝑠)) = 𝐩T (𝑠)𝐏−1 𝑗 = 1, 2, … , 𝑚 𝑚 𝑗,

𝑗 = 1, 2, … , 𝑚

(6)

𝑗 = 1, 2, … , 𝑚

𝑚 ∑ | | |𝓁𝑗 (𝑥)| | | 𝑗=1

𝑢(𝑥) = (sin 𝑥 + 𝑒𝑥 )|𝑥|3.5 ,

𝑥 ∈ [−1, 1]

In computation, uniformly distributed nodes are used, and the support domain of each evaluation point is chosen such that the support domain contains five nodes. Fig. 1 gives the 2-norm condition number of the matrix Pm and the Lebesgue constant 𝜆 in terms of the nodal spacing h. The results are computed at a sample point 𝑥 = 0. Fig. 1(a) shows that the condition number in the polynomial PIM increases quickly as h decreases, while Fig. 1(b) shows that the Lebesgue constant in the Lagrange PIM is about 2.2 in all cases.

Inserting Eq. (2) into Eq. (1) produces

𝑗=1

𝑘=1,𝑘≠𝑗

𝑠 − 𝑠𝑘 , 𝑠𝑗 − 𝑠𝑘

On the other hand, the condition number of the moment matrix Pm affects the stability and performance of the polynomial PIM. To show the stability and performance, the two PIMs are considered to fit the function

(2)

𝑚 ∑

(5)

Since no inverse matrix appears in Eq. (6), the Lagrange PIM in Eq. (5) is much simpler and requires less CPU time than the polynomial PIM in Eq. (3). A useful measure of the stability and performance of the Lagrange interpolation is the Lebesgue constant [46],

where 𝜇 = [𝜇1 , 𝜇2 , … , 𝜇𝑚 ]T with 𝜇𝑗 = 𝜇(𝑠𝑗 ), and [ ( ) ( ) ( )]T 𝐏𝑚 = 𝐩 𝑠1 , 𝐩 𝑠2 , … , 𝐩 𝑠𝑚

𝜇(𝑥(𝑠)) =

𝑚 ∏

Φ𝑗 (𝑥) = 𝓁𝑗 (𝑥),

where 𝐩T (𝑠) = [1, 𝑠, 𝑠2 , … , 𝑠𝑚−1 ] and 𝐚 = [𝑎1 , 𝑎2 , … , 𝑎𝑚 ]T . To determine the coefficient vector a, a support domain is formed for the point x(s) ∈ Γ, and m nodes y1 , y2 ,…, ym are included in the support domain [26,31–35]. Then, the vector a can be obtained by collocation, i.e., 𝐩T (𝑠𝑗 )𝐚 = 𝜇(𝑠𝑗 ) for 𝑗 = 1, 2, … , 𝑚. Thus 𝐚 = 𝐏−1 𝑚 𝜇

𝑥 (𝑠 ) ∈ Γ

Because a polynomial of degree 𝑚 − 1 interpolating m distinct data points is unique, the interpolation polynomials in Eqs. (3) and (5) are the same. Then,

The point interpolation in the BPIM is constructed on the boundary Γ of a two-dimensional domain Ω. Let {𝑦𝑗 }𝑁 be N nodes on Γ, and let 𝑗=1 h be the nodal spacing. The curvilinear coordinates of x ∈ Γ and yj on Γ are denoted as s and sj , respectively. The point interpolation for a function 𝜇(x) is [26] 𝑚 ∑

𝓁𝑗 (𝑥(𝑠))𝜇𝑗 ,

where the shape function 𝓁 j (x(s)) is the Lagrange basis polynomial,

2.1. Polynomial PIM

𝜇(𝑥(𝑠)) =

𝑚 ∑

(3)

(4)

Fig. 1. Results of (a) condition number of the moment matrix Pm in the polynomial PIM and (b) Lebesgue constant 𝜆 in the Lagrange PIM. 12

L. Chen and X. Li

Engineering Analysis with Boundary Elements 103 (2019) 11–21

prescribed boundary condition, ∫Γ

𝜇(𝑥)

𝜕 2 𝑢∗ (𝑥, 𝑦) dΓ𝑥 = 𝑔 (𝑦), 𝜕 𝑛𝑦 𝜕 𝑛𝑥

𝑦∈Γ

(9)

To solve Eq. (9) numerically, let {𝑦𝑗 }𝑁 be N nodes on Γ. If 𝑦(𝑠𝐼1 ), 𝑗=1 𝑦(𝑠𝐼2 ), … , 𝑦(𝑠𝐼𝑚 ) are the m nodes with support domains that include the point x(s) ∈ Γ, then from Eq. (5) the Lagrange PIM for 𝜇(x) can be represented as 𝜇(𝑥(𝑠)) =

𝐼𝑚 ∑ 𝑡=𝐼1

Φ𝑡 (𝑥(𝑠))𝜇𝑡 =

𝑁 ∑ 𝑡=1

Φ𝑡 (𝑥(𝑠))𝜇𝑡 ,

𝑥 (𝑠 ) ∈ Γ

(10)

where 𝑚 ⎧ ∏ 𝑠 − 𝑠𝐼𝑙 , ⎪ Φ𝑡 (𝑥(𝑠)) = ⎨ 𝑙=1,𝑙≠𝑗 𝑠𝐼𝑗 − 𝑠𝐼𝑙 ⎪ ⎩ 0,

{ } 𝑡 ∉ 𝐼1 , 𝐼2 , … , 𝐼𝑚 ,

𝑡 = 1, 2, … , 𝑁

Substituting Eq. (10) into Eq. (9) leads to ( ) 𝑁 ∑ 𝜕 2 𝑢∗ 𝑥, 𝑦𝑗 ( ) 𝜇 𝑡 Φ 𝑡 (𝑥 ) dΓ𝑥 = 𝑔 𝑦𝑗 , 𝑗 = 1, 2, … , 𝑁 ∫Γ 𝜕 𝑛 𝜕 𝑛 𝑦 𝑥 𝑗 𝑡=1

Fig. 2. Convergence and stability of the polynomial PIM and the Lagrange PIM.

Let 𝑒 (𝑢 ) =

{ } 𝑡 = 𝐼𝑗 ∈ 𝐼1 , 𝐼2 , … , 𝐼𝑚 ,

( ) ( ) 𝑁𝑡 1 ∑ || 𝑢 𝑥𝑗 − 𝑢ℎ 𝑥𝑗 || | | ( ) | 𝑁𝑡 𝑗=1 || 𝑢 𝑥𝑗 |

(11)

Dividing Γ into NC non-overlapping integration cells Γ𝓁 [48,49] yields 𝑁 ∑

where xj , 𝑗 = 1, 2, … , 𝑁𝑡 , are test points, and u(xj ) and uh (xj ) are analytical and numerical results, respectively. To compare the performance of the polynomial PIM and the Lagrange PIM, Fig. 2 gives the errors of the two PIMs in terms of the nodal spacing h. The error results are computed using 10001 test points 𝑥𝑗 = 0.0002(𝑗 − 1) − 1, 𝑗 = 1, 2, … , 10, 001. It can be found that the error of the Lagrange PIM decreases monotonously as h decreases. However, when h is small enough, i.e., a large number of nodes are used, the error of the polynomial PIM does not decrease along with the decrease of h. This may be due to the large condition number as shown in Fig. 1(a). Besides, to compute shape functions, the polynomial PIM consumes 2.960 s, while the Lagrange PIM consumes 0.578 s. Clearly, the Lagrange PIM consumes much less CPU time than the polynomial PIM. Therefore, the Lagrange PIM outperforms the polynomial PIM in terms of convergence, stability and computational efficiency.

𝑡=1

𝜇𝑡

𝑁𝐶 ∑ 𝓁=1

∫Γ𝓁

Φ 𝑡 (𝑥 )

( ) 𝜕 2 𝑢∗ 𝑥, 𝑦𝑗 𝜕 𝑛𝑦𝑗 𝜕 𝑛𝑥

( ) dΓ𝑥 = 𝑔 𝑦𝑗 ,

𝑗 = 1, 2, … , 𝑁

(12)

which contains a set of N linear algebraic equations. As discussed in Section 4, Eq. (12) contains regular integrals, logarithmic weakly singular integrals, hypersingular integrals and nearly hypersingular integrals. Regular integrals can be computed using Gaussian quadrature, while logarithmic weakly singular integrals can be computed using the integration technique in Ref. [40]. In Sections 5 and 6, a power series expansion method and a cell subdivision method are developed to compute hypersingular and nearly hypersingular integrals, respectively. Solving Eq. (12) obtains all unknowns 𝜇 t . Finally, the acoustical potential u in Ω can be approximated using Eq. (8) as 𝑢ℎ (𝑦) =

𝑁 ∑ 𝑡=1

𝜇𝑡

𝑁𝐶 ∑ 𝓁=1

∫Γ𝓁

Φ 𝑡 (𝑥 )

𝜕𝑢∗ (𝑥, 𝑦) dΓ𝑥 , 𝜕𝑛𝑥

𝑦∈Ω

(13)

3. An improved BPIM for exterior acoustic radiation problem 4. Analysis of boundary integrals

This paper considers an exterior acoustic radiation problem in a twodimensional infinite acoustic fluid domain Ω radiated by a closed radiator with boundary Γ. This problem can be represented as ⎧ Δ𝑢(𝑥) + 𝑘2 𝑢(𝑥) = 0, ⎪ ⎪ 𝜕𝑢(𝑥) = 𝑔 (𝑥), ⎨ 𝜕𝑛𝑥 ⎪ ) √ ( ⎪ lim 𝑟 𝜕𝑢 − i𝑘𝑢 = 0, ⎩ 𝑟→∞ 𝜕𝑟

In Eq. (12), we have to compute the following boundary integrals ( ) 𝜕 2 𝑢∗ 𝑥, 𝑦𝑗 𝐼= Φ 𝑡 (𝑥 ) dΓ𝑥 (14) ∫Γ𝓁 𝜕 𝑛𝑦 𝜕 𝑛𝑥

( )T 𝑥 = 𝑥1 , 𝑥2 ∈ Ω

𝑗

𝑥 ∈ Γ = 𝜕Ω

(7)

where 𝑥 = (𝑥1 , 𝑥2 )T , 𝑦𝑗 = (𝑦𝑗1 , 𝑦𝑗2 )T , 𝑡 = 1, 2, … , 𝑁, 𝑗 = 1, 2, … , 𝑁, and 𝓁 = 1, 2, … , 𝑁𝑐 . If yj ∉ Γ𝓁 , the integral in Eq. (14) is regular and can be computed using Gaussian quadrature. However, if yj ∈ Γ𝓁 , the integral in Eq. (14) may be singular. Let 𝑟 = |𝑥 − 𝑦𝑗 |. Then, ( ) 𝜕𝑢∗ 𝑥, 𝑦𝑗 ( ) i i𝑘 𝜕𝑟 𝑢∗ 𝑥, 𝑦𝑗 = 𝐻0(1) (𝑘𝑟), = − 𝐻1(1) (𝑘𝑟) 4 𝜕𝑛𝑥 4 𝜕𝑛𝑥

𝑟 = |𝑥|

where u(x) is the unknown acoustic pressure, k is the acoustic wavenumber, nx is the unit outward normal to Γ at the point x, g(x) is a given √ boundary function on Γ, and i = −1. Using the acoustic double layer potential, the solution of the exterior acoustic radiation problem (7) can be expressed as [2,3] 𝑢 (𝑦 ) =

∫Γ

𝜇(𝑥)

𝜕𝑢∗ (𝑥, 𝑦) dΓ𝑥 , 𝜕𝑛𝑥

𝑦∈Ω

where 𝐻0(1) and 𝐻1(1) are the Hankel functions of the first kind of orders zero and one, respectively. Thus, [ ] ( ) (1 ) 𝜕 2 𝑢∗ 𝑥, 𝑦𝑗 i𝑘 𝜕𝐻1 (𝑘𝑟) 𝜕𝑟 𝜕𝑟 𝜕2 𝑟 =− + 𝐻1(1) (𝑘𝑟) 𝜕 𝑛𝑦𝑗 𝜕 𝑛𝑥 4 𝜕𝑟 𝜕𝑛𝑦𝑗 𝜕𝑛𝑥 𝜕 𝑛𝑦𝑗 𝜕 𝑛𝑥

(8)

where 𝜇 denotes the jump through Γ of u, and u∗ (x, y) is the fundamental solution of two-dimensional Helmholtz equation. According to the continuousness of the normal derivative of the double layer potential, the unknown 𝜇 can be solved from the following BIE by invoking the 13

L. Chen and X. Li

Engineering Analysis with Boundary Elements 103 (2019) 11–21

where

Let ( ) 𝑥 − 𝑦𝑗 ⋅ 𝑛𝑦𝑗

( ) 𝑥 − 𝑦𝑗 ⋅ 𝑛𝑥 𝜕𝑟 𝜕𝑟 = , =− 𝜕𝑛𝑥 𝑟 𝜕𝑛𝑦𝑗 𝑟 ] [( ) ][( ) 𝑥 − 𝑦𝑗 ⋅ 𝑛𝑥 𝑥 − 𝑦𝑗 ⋅ 𝑛𝑦𝑗 + 𝑟3

,

𝑛𝑥 ⋅ 𝑛𝑦𝑗 𝜕2 𝑟 =− 𝜕 𝑛𝑦𝑗 𝜕 𝑛𝑥 𝑟

𝐼3 = 𝐼11 + 𝐼21 = −

1

𝜕𝑟

1 = 𝑘𝐻0(1) (𝑘𝑟) − 𝐻1(1) (𝑘𝑟), 𝑟

When yj ∈ Γ𝓁 , the integrals in I12 and I23 are regular, the integral in I3 is weakly singular, while the integral in I22 is hypersingular. Regular integrals can be computed using Gaussian quadrature, but singular integrals need to be tackled efficiently. To compute integrals on Γ𝓁 numerically, let 𝜉 ∈ [−1, 1] be the intrinsic curvilinear coordinate of 𝑥 = (𝑥1 (𝜉), 𝑥2 (𝜉))T ∈ Γ𝓁 . The value of 𝜉 at yj is expressed as 𝜉 j . Then, the weakly singular integral I3 in Eq. (20) reads

𝑛𝑥 ⋅ 𝑛𝑦𝑗 𝜕2 𝑟 1 𝜕𝑟 𝜕𝑟 =− − 𝜕 𝑛𝑦𝑗 𝜕 𝑛𝑥 𝑟 𝑟 𝜕𝑛𝑥 𝜕𝑛𝑦𝑗

we have ( ) 𝜕 2 𝑢∗ 𝑥, 𝑦𝑗 𝜕 𝑛𝑦𝑗 𝜕 𝑛𝑥

=−

[ ( )] i𝑘 𝜕𝑟 𝜕𝑟 1 1 𝜕𝑟 𝜕𝑟 1 𝑘𝐻0( ) (𝑘𝑟) − 𝐻1( ) (𝑘𝑟) 2 + 𝑛𝑥 ⋅ 𝑛𝑦𝑗 4 𝜕𝑛𝑥 𝜕𝑛𝑦𝑗 𝑟 𝜕𝑛𝑥 𝜕𝑛𝑦𝑗 (15)

𝐼3 =

Let zmax be a positive number and let 𝐶 = 0.5772156649 … be the Euler’s constant. When yj ∈ Γ𝓁 , the integration cell Γ𝓁 can be chosen such that |kr| ≤ zmax . The choice 𝑧max = 20 can yield acceptable results [21]. Then [21] 𝐻0(1) (𝑘𝑟) =

𝐼12

𝐼21 𝐼22

√ [𝑥1 (𝜉) − 𝑦𝑗1 ]2 + [𝑥2 (𝜉) − 𝑦𝑗2 ]2, and q(𝜉) is the non-singular

1

[ ( )] ( ) ln [𝑟(𝜉)] 𝑔 (𝜉) − 𝑔 𝜉𝑗 d𝜉 + 𝑔 𝜉𝑗

1

and then computed using Gaussian quadrature, √( )2 ( )2 | d𝑥1 (𝜉) + d𝑥d2𝜉(𝜉) || is the Jacobean at yj . d𝜉 |𝜉=𝜉𝑗

where

𝐽𝑗 =

5. Treatment of hypersingular integrals (17)

When 𝑦𝑗 = (𝑦𝑗1 , 𝑦𝑗2 )T ∈ Γ𝓁 , the integral I22 in Eq. (19) is hypersingular. With the aid of the power series expansion technique [43], a power series expansion method (PSEM) is developed in this section to compute hypersingular integrals in the current BPIM. In Ref. [21], the PSEM is developed to compute strongly singular integrals in the BEFM, and numerical results verify the efficiency of the PSEM. We assume that the curvilinear coordinates of 𝑥 = (𝑥1 , 𝑥2 )T ∈ Γ𝓁 and yj are 𝜉 ∈ [−1, 1] and 𝜉 j , respectively. Then, the hypersingular integral I22 in Eq. (19) reads

where ( ) 𝑘𝑟 𝜕𝑟 𝜕𝑟 2i𝑘 Φ𝑡 (𝑥) ln dΓ 𝜋 ∫Γ𝓁 2 𝜕𝑛𝑥 𝜕𝑛𝑦𝑗 𝑥

( ) 2i𝐶 ∑ (−1)𝑚 𝑘𝑟 2𝑚 =𝑘 Φ (𝑥 ) 1 + + 2 ∫Γ𝓁 𝑡 𝜋 2 𝑚=1 (𝑚!) [ ( )]} 𝑚 ( ) ∑1 2i 𝜕𝑟 𝜕𝑟 𝑘𝑟 × 1+ ln dΓ𝑥 +𝐶 − 𝜋 2 𝑗 𝜕𝑛 𝑥 𝜕𝑛𝑦𝑗 𝑗=1 ) ( ( ) i𝑘 𝜕𝑟 𝜕𝑟 𝑘𝑟 =− Φ𝑡 (𝑥) ln 2 + 𝑛𝑥 ⋅ 𝑛𝑦𝑗 dΓ𝑥 𝜋 ∫Γ𝓁 2 𝜕𝑛𝑥 𝜕𝑛𝑦𝑗 ( ) 2i 1 𝜕𝑟 𝜕𝑟 = Φ 𝑡 (𝑥 ) 2 + 𝑛𝑥 ⋅ 𝑛𝑦𝑗 dΓ𝑥 𝑘 ∫Γ𝓁 𝜕𝑛𝑥 𝜕𝑛𝑦𝑗 𝜋𝑟2

(21)

⎤ ⎡ 𝑟(𝜉) ⎥ ln ⎢ d𝜉 ∫−1 ∫−1 ⎢ 𝐽 ||𝜉 − 𝜉 || ⎥ 𝑗| ⎦ ⎣ 𝑗| ( ){ ( ) ( )[ ( ) ] ( )[ ( ) ]} + 𝑔 𝜉𝑗 2 ln 𝐽𝑗 + 1 + 𝜉𝑗 ln 1 + 𝜉𝑗 −1 + 1 − 𝜉𝑗 ln 1 − 𝜉𝑗 − 1 (22)

𝐼3 =

Eqs. (16) and (17) indicate that the Hankel functions in Eq. (15) are singular when 𝑟 = 0, i.e., 𝑥 = 𝑦𝑗 . Inserting Eq. (15) into Eq. (14) yields ( ) 𝜕 2 𝑢∗ 𝑥, 𝑦𝑗 ) i𝑘 ( 𝐼 + 𝐼12 + 𝐼21 + 𝐼22 + 𝐼23 (18) Φ 𝑡 (𝑥 ) dΓ𝑥 = − ∫Γ𝓁 𝜕 𝑛𝑦𝑗 𝜕 𝑛𝑥 4 11

{

ln [𝑟(𝜉)]𝑞 (𝜉)d𝜉

part of the integrand. Applying the integration technique in Ref. [40], the weakly singular integral can be regularized as

(16)

𝐼11 =

1

∫−1

where 𝑟(𝜉) =

( ) 𝑘𝑟 2i ln +1 𝜋 2 [ )] ( ∞ 𝑚 ( ) ( ) ∑ 2i𝐶 ∑ (−1)𝑚 𝑘𝑟 2𝑚 2i 1 𝑘𝑟 + 1+ ln +𝐶 − + 2 𝜋 2 𝜋 2 𝑗 𝑚=1 (𝑚!) 𝑗=1

( ) 𝑘𝑟 1 (1) i𝑘 2i 𝐻1 (𝑘𝑟) = − ln 𝑟 𝜋 2 𝑘𝜋𝑟2 { ∞ i(2𝐶 − 1) ∑ (−1)𝑚 ( 𝑘𝑟 )2𝑚 𝑘 + 1+ + 2 𝜋 𝑚!(𝑚 + 1)! 2 𝑚=1 [ ( )]} 𝑚 −1 ( ) ∑ 2i 1 1 𝑘𝑟 × 1+ ln − +𝐶 − 𝜋 2 𝑗 + 1 2𝑚 + 2 𝑗=0

(20)

Then Eq. (18) can be rewritten as ( ) 𝜕 2 𝑢∗ 𝑥, 𝑦𝑗 ) ] i𝑘 [( Φ (𝑥) dΓ𝑥 = − 𝐼12 + 𝐼23 + 𝐼3 + 𝐼22 ∫Γ𝓁 𝑡 𝜕 𝑛𝑦𝑗 𝜕 𝑛𝑥 4

Then, according to 𝜕𝐻1( ) (𝑘𝑟)

( ) i𝑘 𝑘𝑟 Φ𝑡 (𝑥) ln 𝑛 ⋅ 𝑛 dΓ ∫ 𝜋 Γ𝓁 2 𝑥 𝑦𝑗 𝑥

𝐼22 =

1 𝑓 (𝜉) d𝜉 ∫−1 𝑟2 (𝜉)

(23)

where f(𝜉) is the non-singular part of the integrand, and √ [ ]2 [ ]2 𝑟(𝜉) = 𝑥1 (𝜉) − 𝑦𝑗1 + 𝑥2 (𝜉) − 𝑦𝑗2

∞

(24)

is the global distance between x and yj . Let the local distance between x and yj be | | 𝜌(𝜉) = |𝜉 − 𝜉𝑗 |, | |

𝜉 ∈ [−1, 1]

(25)

We assume that r can be expanded in 𝜌 as [43] (19)

𝑟=𝜌

∞ ∑ 𝑚=0

{

∞ i(2𝐶 − 1) ∑ (−1)𝑚 (𝑘𝑟)2𝑚 [ 2i ( (𝑘𝑟) 𝑘 Φ 𝑡 (𝑥 ) 1 + 1+ + ln ∫ 2 Γ𝓁 𝜋 𝑚!(𝑚 + 1)! 2 𝜋 2 𝑚=1 )]}( ) 𝑚 −1 ∑ 1 1 𝜕𝑟 𝜕𝑟 − 2 +𝑛𝑥 ⋅ 𝑛𝑦𝑗 dΓ𝑥 + 𝐶− 𝑗 + 1 2 𝑚 + 2 𝜕𝑛 𝑥 𝜕𝑛𝑦𝑗 𝑗=0

𝐼23 = −

𝑐 𝑚 𝜌𝑚

(26)

To determine cm in Eq. (26), quadratic integration cells are adopted. Then, for any 𝑥 = (𝑥1 , 𝑥2 )T ∈ Γ𝓁 , ( ) 𝑥𝑡 (𝜉) = 0.5𝜉(𝜉 − 1)𝑥2𝑡 𝓁−1 + 1 − 𝜉 2 𝑥2𝑡 𝓁 + 0.5𝜉(𝜉 + 1)𝑥2𝑡 𝓁+1 ,

𝑡 = 1, 2 (27)

14

L. Chen and X. Li

Engineering Analysis with Boundary Elements 103 (2019) 11–21

where (𝑥21𝓁−1 , 𝑥22𝓁−1 )T and (𝑥21𝓁+1 , 𝑥22𝓁+1 )T are the ends of Γ𝓁 and (𝑥21𝓁 , 𝑥22𝓁 )T is its mid-point. Eq. (27) can be rearranged as 𝑥𝑡 (𝜉) = 0.5𝑎𝑡 𝜉 2 − 𝑏𝑡 𝜉 + 𝑥2𝑡 𝓁 , 𝑥2𝑡 𝓁−1

2𝑥2𝑡 𝓁

where 𝐵𝑡0 = 𝑓𝑡 (0) = 𝑓 (𝜉𝑗 )∕𝑐02 , and Bt1 , Bt2 and Bt3 are obtained by solving

𝑥2𝑡 𝓁+1

0.5(𝑥2𝑡 𝓁−1

𝑥2𝑡 𝓁+1 )

where 𝑎𝑡 = − + and 𝑏𝑡 = − for 𝑡 = 1 and 2. Since the value of 𝜉 at 𝑦𝑗 = (𝑦𝑗1 , 𝑦𝑗2 )T ∈ Γ𝓁 is 𝜉 j , we have 𝑦𝑗𝑡 = 0.5𝑎𝑡 𝜉𝑗2 − 𝑏𝑡 𝜉𝑗 + 𝑥2𝑡 𝓁 ,

𝑡 = 1, 2

𝑡 𝐼22 = lim

𝜌𝑡

𝜀𝑡 →0 ∫𝜀

(

= 𝐵𝑡0

√ √ √ 2 √ 2 √∑ [ ∑ {( ]2 √ )[ ( ) ]}2 √ 𝑟= 𝑥𝑡 (𝜉) − 𝑦𝑗𝑡 = √ 𝜉 − 𝜉𝑗 0.5𝑎𝑡 𝜉 + 𝜉𝑗 − 𝑏𝑡

𝜌𝑡

lim

Let 𝜌,𝜉 = d𝜌∕d𝜉. Then, using Eq. (25) yields

𝜀𝑡 →0 ∫𝜀

𝐵 𝑡 0 + 𝐵 𝑡 1 𝜌 + 𝐵 𝑡 2 𝜌2 + 𝐵 𝑡 3 𝜌3

𝑡

𝑐1 1 − 𝑐 0 𝜌𝑡

)

𝜌2

d𝜌

𝐵𝑡3 𝜌2𝑡 ( ) + 𝐵𝑡1 ln 𝑐0 𝜌𝑡 + 𝐵𝑡2 𝜌𝑡 + 2

𝜉 > 𝜉𝑗 𝜉 < 𝜉𝑗

𝑡

𝑐 1 1 d𝜌 = 1 − , 𝑐 0 𝜌𝑡 𝜌2

𝜌𝑡

lim

𝜀𝑡 →0 ∫𝜀

𝑡

( ) 1 d𝜌 = ln 𝑐0 𝜌𝑡 𝜌

( )𝑐 ( ) ( ) 1 1 𝐼22 = 𝐵10 + 𝐵20 1 − 𝐵10 − 𝐵20 + 𝐵12 1 + 𝜉𝑗 + 𝐵22 1 − 𝜉𝑗 𝑐0 1 + 𝜉𝑗 1 − 𝜉𝑗 ( )2 ( )2 𝐵23 1 − 𝜉𝑗 [ ( )] [ ( )] 𝐵13 1 + 𝜉𝑗 +𝐵11 ln 𝑐0 1 + 𝜉𝑗 +𝐵21 ln 𝑐0 1 − 𝜉𝑗 + + 2 2

√ √ 2 √∑ { [ ) ( ]}2 𝑟=√ 𝜌,𝜉 𝜌 0.5𝑎𝑡 2𝜉𝑗 + 𝜌,𝜉 𝜌 − 𝑏𝑡 𝑡=1

(36)

√ √ 2 [ ] √∑ ( ) = 𝜌√ 𝑎2𝑡 𝜉𝑗2 − 2𝑎𝑡 𝑏𝑡 𝜉𝑗 + 𝑏2𝑡 + 𝑎2𝑡 𝜉𝑗 − 𝑎𝑡 𝑏𝑡 𝜌,𝜉 𝜌 + 0.25𝑎2𝑡 𝜌2

Remark 5.1. In the case of 𝜉𝑗 = 0, i.e., yj is the center of the cell Γ𝓁 , Eq. (36) reduces to ) ( ( ) 𝑐1 ( ) ) 1( 𝐼22 = 𝐵10 + 𝐵20 − 1 + 𝐵12 + 𝐵22 + 𝐵11 +𝐵21 ln 𝑐0 + 𝐵13 + 𝐵23 𝑐0 2

𝑡=1

√ 𝑑 0 + 𝑑 1 𝜌 + 𝑑 2 𝜌2

(28)

Remark 5.2. According to Eqs. (19) and (23), the computation of 𝐵𝑡0 = 𝜕𝑟 𝑓 (𝜉𝑗 )∕𝑐02 involves lim 𝜕𝑛 and lim 𝜕𝑛𝜕𝑟 . From Eq. (27), we have

where 𝑑0 = (𝑎21 + 𝑎22 )𝜉𝑗2 − 2(𝑎1 𝑏1 + 𝑎2 𝑏2 )𝜉𝑗 + 𝑏21 + 𝑏22 , 𝑑1 = [(𝑎21 + 𝑎22 )𝜉𝑗 − (𝑎1 𝑏1 + 𝑎2 𝑏2 )]𝜌,𝜉 , and 𝑑2 = 0.25(𝑎21 + 𝑎22 ). Eqs. (26) and (28) imply √ ∑∞ 𝑚 𝑑0 + 𝑑1 𝜌 + 𝑑2 𝜌2 . Therefore, the coefficients cm can be de𝑚=0 𝑐𝑚 𝜌 = termined approximately as 𝑐0 =

√

𝑑0 ,

𝑐1 =

𝑑1 , 2𝑐0

𝑐2 =

) 1 ( 𝑑 − 𝑐12 , 2𝑐0 2

𝑐𝑚 = 0,

𝜌→0

𝑥𝑡 − 𝑦𝑗𝑡 𝜕𝑟 = 𝜕𝑥𝑡 𝑟

𝑚 = 3, 4, …

(30)

𝑥𝑡 (𝜉) = 0.5(1 − 𝜉)𝑥𝓁𝑡 + 0.5(1 + 𝜉)𝑥𝓁+1 , 𝑡

(31)

where

𝜉𝑗 𝑓 (𝜉) 1 d𝜉, ∫−1 𝜌2 (𝑐 + 𝑐 𝜌 + 𝑐 𝜌2 )2 0 1 2

2 𝐼22 =

1 𝑓 (𝜉) 1 d𝜉 ∫𝜉𝑗 𝜌2 (𝑐 + 𝑐 𝜌 + 𝑐 𝜌2 )2 0 1 2

Eq. (25) shows 𝜌 = 𝜉𝑗 − 𝜉 for 𝜉 < 𝜉 j and 𝜌 = 𝜉 − 𝜉𝑗 for 𝜉 > 𝜉 j . Then, using 1 and 𝜉 = 𝜉 + 𝜌 in 𝐼 2 , we obtain 𝜉 = 𝜉𝑗 − 𝜌 in 𝐼22 𝑗 22 𝑡 𝐼22 = lim

𝜀𝑡 →0 ∫𝜀

𝜌𝑡 𝑡

1 𝑓𝑡 (𝜌)d𝜌, 𝜌2

𝑡 = 1, 2

𝑡 = 1, 2

(𝑥𝓁1 , 𝑥𝓁2 )T

𝑡 = 1, 2

(37)

(𝑥𝓁+1 , 𝑥𝓁+1 )T 1 2

and are the ends of Γ𝓁 . Substituting √ Eq. (37) into Eq. (24) yields 𝑟 = 𝜌 𝑑0 , where 𝑑0 = 𝑏21 + 𝑏22 with 𝑏𝑡 = √ ∑ 𝑚 0.5(𝑥𝓁𝑡 − 𝑥𝓁+1 ) for 𝑡 = 1 and 2. Then, Eq. (26) implies ∞ 𝑑0 𝑡 𝑚=0 𝑐𝑚 𝜌 = √ and thus 𝑐0 = 𝑑0 and 𝑐𝑚 = 0 for 𝑚 = 1, 2, …. Finally, the PSEM formula for the hypersingular integral in Eq. (23) in the case of linear integration cells is

where 1 𝐼22 =

𝜕𝑟 𝜕𝑟 =− , 𝜕𝑦𝑗𝑡 𝜕𝑥𝑡

Remark 5.3. If linear integration cells are adopted, then Eq. (27) becomes

1

𝑓 (𝜉) 𝑓 (𝜉) 1 1 2 d𝜉 = d𝜉 = 𝐼22 + 𝐼22 ∫−1 𝑟2 (𝜉) ∫−1 𝜌2 (𝑐 + 𝑐 𝜌 + 𝑐 𝜌2 )2 0 1 2

𝑦𝑗

( ) 2 2 ∑ ∑ 𝑎𝑡 𝜉𝑗 − 𝑏𝑡 𝜌,𝜉 𝜕𝑟 𝜕𝑟 = lim 𝑛𝑥𝑡 = 𝑛𝑥𝑡 𝜌→0 𝜕𝑛𝑥 𝜌→0 𝜕𝑥𝑡 𝑐0 𝑡=1 𝑡=1 ( ) 2 2 ∑ ∑ 𝑎𝑡 𝜉𝑗 − 𝑏𝑡 𝜌,𝜉 𝜕𝑟 𝜕𝑟 lim = lim 𝑛 =− 𝑛𝑦𝑗𝑡 𝜌→0 𝜕𝑛𝑦 𝜌→0 𝜕𝑦𝑗𝑡 𝑦𝑗𝑡 𝑐0 𝑗 𝑡=1 𝑡=1

Inserting Eq. (30) into Eq. (23) yields 𝐼22 =

𝜌→0

( ) 0.5𝑎𝑡 𝜌 + 𝑎𝑡 𝜉𝑗 − 𝑏𝑡 𝜌,𝜉 = √ , 𝑑 0 + 𝑑 1 𝜌 + 𝑑 2 𝜌2

lim

Then, r can be expanded approximately in 𝜌 as 𝑟 = 𝑐 0 𝜌 + 𝑐 1 𝜌2 + 𝑐 2 𝜌3

𝑥

Hence,

(29)

1

(35)

Finally, substituting Eq. (35) into Eq. (31), the PSEM formula for the hypersingular integral in Eq. (23) in the case of quadratic integration cells is

from which we have 𝜉 = 𝜉𝑗 + 𝜌,𝜉 𝜌 and (𝜌,𝜉 )2 = 1. Hence,

=𝜌

(34)

where we have used the following results [43],

𝑡=1

{ 𝜉 − 𝜉𝑗 1, 𝜌,𝜉 = = −1, 𝜌

( )3 ( ) 0.25𝜌𝑡 ⎤ ⎡𝐵𝑡1 ⎤ ⎡𝑓𝑡 0.25𝜌𝑡 − 𝐵𝑡0 ⎤ ) ( )3 ⎥ ⎢ ⎥ ⎢ ( 0.50𝜌𝑡 ⎥ 𝐵𝑡2 = 𝑓𝑡 (0.50𝜌𝑡 ) − 𝐵𝑡0 ⎥ ⎢ ⎥ ⎢ ⎥ ( )3 0.75𝜌𝑡 ⎥⎦ ⎣𝐵𝑡3 ⎦ ⎣𝑓𝑡 0.75𝜌𝑡 − 𝐵𝑡0 ⎦

Inserting Eq. (33) into Eq. (32) yields

Consequently, Eq. (24) becomes

𝑡=1

( )2 0.25𝜌𝑡 ( )2 0.50𝜌𝑡 ( )2 0.75𝜌𝑡

⎡0.25𝜌 𝑡 ⎢ ⎢0.50𝜌𝑡 ⎢0.75𝜌 ⎣ 𝑡

𝑡 = 1, 2

( ) ( ) 1 1 − 𝐵20 + 𝐵12 1 + 𝜉𝑗 + 𝐵22 1 − 𝜉𝑗 1 + 𝜉𝑗 1 − 𝜉𝑗 ( )2 ( )2 𝐵23 1 − 𝜉𝑗 [ ( )] [ ( )] 𝐵13 1 + 𝜉𝑗 + +𝐵11 ln 𝑐0 1 + 𝜉𝑗 +𝐵21 ln 𝑐0 1 − 𝜉𝑗 + 2 2

𝐼22 = −𝐵10

(32)

where 𝜌𝑡 = 1 + (−1)𝑡+1 𝜉𝑗 , 𝜀t is a local infinitesimal distance [43], 𝑓1 (𝜌) = 𝑓 (𝜉𝑗 − 𝜌)∕(𝑐0 + 𝑐1 𝜌 + 𝑐2 𝜌2 )2 , and 𝑓2 (𝜌) = 𝑓 (𝜉𝑗 + 𝜌)∕(𝑐0 + 𝑐1 𝜌 + 𝑐2 𝜌2 )2 . As in Ref. [43], we can expand ft (𝜌) as 𝑓𝑡 (𝜌) = 𝐵𝑡0 + 𝐵𝑡1 𝜌 + 𝐵𝑡2 𝜌2 + 𝐵𝑡3 𝜌3 ,

𝑡 = 1, 2

(33) Fig. 3. Dividing the integration cell Γ𝓁 corresponding to the boundary node yj . 15

L. Chen and X. Li

Engineering Analysis with Boundary Elements 103 (2019) 11–21

Fig. 4. Results of u(r, 𝜋/4) obtained using (a) 𝑘 = 1 and (b) 𝑘 = 10 for the first example.

6. Treatment of nearly hypersingular integrals According to Eqs. (18) and (19), when the boundary node yj is close to the integration cell Γ𝓁 but not on Γ𝓁 , nearly hyper-singularity appears in Eq. (14). Such integral cannot be computed efficiently by using standard Gaussian quadrature, since the integrand varies very sharply within the integration cell. Inspired by the method in Refs. [41,42], a cell subdivision method is developed in this section to compute the nearly hypersingular boundary integral in Eq. (14). Let l be the length of the integration cell Γ𝓁 , and let d be the distance between the boundary node yj and the center of Γ𝓁 . Eq. (19) indicates that the nearly hyper-singularity is caused basically by the term 1/(𝜋r2 ). Thus, if 𝜋d2 < l, the cell Γ𝓁 is divided into four equal sub-cells (see Fig. 3). Otherwise, Γ𝓁 is considered to be a regular integration cell. Finally, Gaussian quadrature is used for all cells and sub-cells. Numerical results in Section 7 show that this integration method can compute the nearly hypersingular integral in Eq. (14) efficiently.

Fig. 5. Error |𝑢 − 𝑢ℎ | at (𝑟, 𝜃) = (4.0, 𝜋∕4) produced by the BPIM and the BEM for the first example.

Fig. 6. Convergence of the BPIM with and without treatment of nearly hypersingular integrals obtained using (a) 𝑘 = 1 and (b) 𝑘 = 10 for the first example. 16

L. Chen and X. Li

Engineering Analysis with Boundary Elements 103 (2019) 11–21

Fig. 7. Convergence of the BPIM using linear and quadratic cells and (a) 𝑘 = 1 and (b) 𝑘 = 10 for the first example.

Fig. 8. Result of u(4, 𝜃) obtained using (a) 𝑘 = 5 and (b) 𝑘 = 5 − i for the second example.

Fig. 9. Result of u(r, 𝜋/4) obtained using (a) 𝑘 = 5 and (b) 𝑘 = 5 − i for the second example.

17

L. Chen and X. Li

Engineering Analysis with Boundary Elements 103 (2019) 11–21

Table 1 Error |𝑢 − 𝑢ℎ | at (𝑟, 𝜃) = (4.0, 𝜋∕4) produced by the BPIM and the BEM for the first example. k

0.500

𝑁 = 64 𝑁 = 128

1.000 −3

BEM BPIM BEM BPIM

2.42×10 1.20×10−5 1.22×10−3 3.27×10−6

1.205 −3

7.74×10 1.80×10−5 4.07×10−3 3.13×10−6

1.930 −2

1.11×10 2.29×10−5 5.90×10−3 3.79×10−6

2.400 −3

2.25×10 1.42×10−3 6.08×10−4 2.18×10−4

2.760 −2

3.25×10 6.14×10−5 4.52×10−4 9.30×10−6

3.100 −2

2.69×10 5.87×10−5 9.62×10−3 8.81×10−6

3.540 −2

4.327 −3

2.70×10 9.22×10−5 7.66×10−3 1.38×10−5

4.47×10 1.22×10−3 1.25×10−3 1.81×10−4

4.880 −2

2.25×10 9.87×10−5 1.17×10−2 1.45×10−5

2.51×10−2 2.87×10−4 9.92×10−3 4.31×10−5

7. Numerical examples The algorithm of the BPIM for the exterior acoustic radiation problem (7) is as follows. 1. Choose N nodes on the boundary Γ of the infinite domain Ω, and divide Γ into NC non-overlapping integration cells Γ𝓁 . 2. Compute shape functions using Eq. (6). 3. Compute integrals in Eq. (12). Regular integrals are computed using Gaussian quadrature, logarithmic weakly singular integrals are computed using Eq. (22), hypersingular integrals are computed using Eq. (36) for quadratic integration cells and Eq. (37) for linear integration cells, and nearly hypersingular integrals are computed using the cell subdivision method in Section 6. 4. Form the linear system given in Eq. (12) and solve the system to obtain all unknowns 𝜇 t . 5. Compute the numerical solution of the acoustic problem from Eq. (13). To illustrate the efficiency of the current improved BPIM, numerical results of some exterior acoustic radiation problems are given in this section. In computations, unless otherwise specified, quadratic cells are adopted for numerical integrations, hypersingular integrals are treated using the PSEM in Section 5, and nearly hypersingular integrals are treated using the cell subdivision method in Section 6. To show the convergence, let

Fig. 10. Error |𝑢 − 𝑢ℎ | at (𝑟, 𝜃) = (4.0, 𝜋∕4) produced by the BPIM and the BEM for the second example.

can be observed that both precision and convergence rate of quadratic cells are much higher than those of linear cells.

( ) ( ) 𝑁𝑡 1 ∑ || 𝑢 𝑥𝑗 − 𝑢ℎ 𝑥𝑗 || 𝑒 (𝑢 ) = | | ( ) | 𝑁𝑡 𝑗=1 || 𝑢 𝑥𝑗 |

Example 7.2. The second example considers a radiation problem of an infinite cylinder with boundary condition 𝜕 𝑢∕𝜕 𝑛 = −𝑒i𝜃 . The analytical solution is [4,21],

where xj , 𝑗 = 1, 2, … , 𝑁𝑡 , are test points in Ω, and u(xj ) and uh (xj ) are analytical and numerical results, respectively.

𝑢(𝑟, 𝜃) =

1

Example 7.1. The first example considers an exterior acoustic radiation problem with boundary condition 𝜕 𝑢∕𝜕 𝑛 = i𝑘. The analytical solution in the polar coordinate system (r, 𝜃) is [4], 𝐻0( ) (𝑘𝑟) 𝐻1(1) (2𝑘)

,

𝑟 ≥ 2,

1

1

𝐻1( ) (2𝑘) − 2𝑘𝐻0( ) (2𝑘)

,

𝑟 ≥ 2,

0 ≤ θ < 2𝜋

In Figs. 8 and 9, the results of u(r, 𝜃) are displayed along the lines 𝑟 = 4 and 𝜃 = 𝜋∕4, respectively. The numerical results are obtained by the current BPIM with 64 boundary nodes. Obviously, numerical results match analytical results. In Fig. 10, errors of u(r, 𝜃) at (4.0, 𝜋/4) produced by the current BPIM and the BEM [4] are shown with various k. The error results are also given in Table 2. Again, the BPIM has less error than the BEM. To show the convergence of the current BPIM, errors obtained by the BPIM and the BEFM [21] are displayed in Fig. 11. The error results are also given in Table 3. The two methods produce convergent results, but the BPIM produces less error than the BEFM.

1

𝑢(𝑟, 𝜃) = −i

2𝑒iθ 𝐻1( ) (𝑘𝑟)

0 ≤ θ < 2𝜋

In Fig. 4, the results of u(r, 𝜃) are displayed along the line 𝜃 = 𝜋∕4. The numerical results are obtained by the current BPIM with 64 boundary nodes. It is true that numerical results match analytical results. In Fig. 5, the error |𝑢 − 𝑢ℎ | of u(4.0, 𝜋/4) produced by the current BPIM with 64 and 128 boundary nodes is shown with various wave number k. The error produced by the BEM [4] is also shown for comparison. The error results are also given in Table 1. It can be observed that the BPIM has much less error than the BEM. To show the importance of the treatment of nearly hypersingular integrals, Fig. 6 displays the convergence of the current BPIM with and without the treatment. The results are obtained by the BPIM using quadratic cells with wave number 𝑘 = 1 and 10. From Fig. 6 it can be observed that this treatment yields much more stable and convergent results. Thus, the treatment of nearly hypersingular integrals is absolutely essential for the efficient application of the current BPIM to exterior acoustic radiation problems. To show the influence of the integration cells, errors obtained by the current BPIM using linear and quadratic cells are displayed in Fig. 7. It

Example 7.3. The third example considers a radiation problem of an infinite cylinder with boundary condition 𝜕 𝑢∕𝜕 𝑛 = 𝑘 cos(4𝜃). The analytical solution is [16], 𝑢(𝑟, 𝜃) = −

𝑘𝐻4(1) (𝑘𝑟) 𝑘𝑎𝐻3(1) (𝑘) − 4𝐻4(1) (𝑘)

cos 4𝜃,

𝑟 ≥ 1,

0 ≤ 𝜃 < 2𝜋

where 𝐻3(1) and 𝐻4(1) are the Hankel functions of the first kind of orders three and four, respectively. In Figs. 12 and 13, the results of u(r, 𝜃) are displayed along the lines 𝑟 = 4 and 𝜃 = 𝜋∕4, respectively. The numerical results are obtained by the current BPIM with 64 boundary nodes. Again, numerical results match analytical results. 18

L. Chen and X. Li

Engineering Analysis with Boundary Elements 103 (2019) 11–21

Fig. 11. Convergence of the BPIM and the BEFM obtained using (a) 𝑘 = 5 and (b) 𝑘 = 5 − i for the second example.

Fig. 12. Result of u(4, 𝜃) obtained using (a) 𝑘 = 2 and (b) 𝑘 = 10 for the third example.

Fig. 13. Result of u(r, 𝜋/4) obtained using (a) 𝑘 = 2 and (b) 𝑘 = 10 for the third example.

19

L. Chen and X. Li

Engineering Analysis with Boundary Elements 103 (2019) 11–21

Table 2 Error |𝑢 − 𝑢ℎ | at (𝑟, 𝜃) = (4.0, 𝜋∕4) produced by the BPIM and the BEM for the second example. k 𝑁 = 64 𝑁 = 128

1.000 BEM BPIM BEM BPIM

1.205 −4

9.22×10 2.56×10−4 5.10×10−4 4.38×10−5

1.500 −3

3.11×10 8.41×10−5 1.70×10−3 1.38×10−5

1.930 −3

7.57×10 5.53×10−5 2.99×10−2 8.73×10−6

2.400 −3

7.89×10 5.92×10−5 2.98×10−3 9.09×10−6

2.690 −3

6.94×10 1.50×10−4 5.05×10−3 2.26×10−5

3.500 −3

1.12×10 1.77×10−3 3.93×10−4 2.65×10−4

3.800 −3

4.305 −3

1.00×10 1.11×10−4 5.31×10−3 1.64×10−5

9.28×10 1.53×10−4 4.34×10−3 2.25×10−5

5.100 −2

1.18×10 1.94×10−3 1.33×10−3 2.84×10−4

9.90×10−3 1.73×10−4 4.12×10−3 2.48×10−5

Table 3 Error e(u) of the BPIM and the BEFM for the second example. N

8

16

32

64

128

256

𝑘=5

3.14×10−1

3.11×10−2

3.90×10−3

5.05×10−4

7.89×10−5

4.50×10−2 1.71×10−1 1.75×10−3

7.32×10−3 8.88×10−3 2.74×10−4

1.22×10−3 8.65×10−4 4.52×10−5

1.69×10−4 9.84×10−5 6.34×10−6

2.34×10−5 1.27×10−5 9.37×10−7

2.46×10−5 2.83×10−6 2.37×10−6 3.47×10−7

𝑘=5−i

BEFM BPIM BEFM BPIM

Fig. 14. Distributions of (a) Re(u), (b) Im(u) and (c) |𝑢 − 𝑢ℎ | obtained by the BPIM with 𝑘 = 2 for the third example.

Fig. 15. Distributions of (a) Re(u), (b) Im(u) and (c) |𝑢 − 𝑢ℎ | obtained by the BPIM with 𝑘 = 10 for the third example.

Figs. 14 and 15 display the distributions of the acoustic pressure u for 𝑘 = 2 and 10, respectively. The associated errors |𝑢 − 𝑢ℎ | are also displayed. When the cylinder boundary is discretized with 64 nodes, the errors for 𝑘 = 2 and 10 are less than 9 × 10−5 and 3.5 × 10−4 , respectively. Hence, the current BPIM produces very good numerical results, and the distributions show the correct acoustic wave propagation.

radiation problems are presented. After analyzing the boundary integrals in the discretized hypersingular BIE, a PSEM and a cell subdivision method are presented to compute acoustic hypersingular and nearly hypersingular integrals, respectively. Both precision and convergence rate of quadratic integration cells are much higher than those of linear cells. Numerical examples demonstrate the accuracy, efficiency and convergence of the current BPIM. The method has great capabilities for exterior acoustic radiation problems. The proposed methods possess some features. Firstly, the improved BPIM requires only boundary nodes and satisfies the Sommerfeld radiation condition naturally. Numerical results reveal that the current BPIM has less error than the BEM and the BEFM. Secondly, shape functions in the Lagrange PIM are explicitly expressed, and thus the Lagrange PIM is much simpler and requires less CPU time than the polynomial

8. Conclusion An improved BPIM is developed in this paper to numerical simulate exterior acoustic radiation problems. To avoid the inversion of the moment matrix in the PIM, a Lagrange PIM is presented by using Lagrange basis polynomials. Based on the acoustic double layer potential and the Lagrange PIM, the formulas of the improved BPIM for exterior acoustic 20

L. Chen and X. Li

Engineering Analysis with Boundary Elements 103 (2019) 11–21

PIM. Numerical results reveal that the Lagrange PIM outperforms the polynomial PIM in terms of convergence, stability and computational efficiency. Thirdly, the PSEM provides a direct and efficient quadrature formula to compute hypersingular integrals. The final computational formulation does not involve any integration and additive term such as function gradient or tangential derivative. As in many BIEs-based meshless methods such as the BEFM, curvilinear coordinates are required in the BPIM to generate shape functions. For complicated boundary problems, it is bothersome to gain curvilinear coordinates. On the other hand, the PSEM for the computation of hypersingular integrals relies on integration cells and some expansions given in Eqs. (26), (30) and (33). Thus, there may have some limitations on the solution domain. The geometrical limitations of the BPIM need to be analyzed theoretically and thoroughly. In addition, the proposed methods encounter some issues to be further studied. Firstly, in the PSEM, the coefficients in Eq. (36) need to be computed by solving the linear system (34). Secondly, although numerical results show that the current BPIM possesses good convergence characteristic, a theoretical convergence analysis needs to be addressed.

[17] Qu WZ, Chen W, Zheng CJ. Diagonal form fast multipole singular boundary method applied to the solution of high-frequency acoustic radiation and scattering. Int J Numer Methods Eng 2017;111:803–15. [18] Fu ZJ, Chen W, Wen PH, Zhang CZ. Singular boundary method for wave propagation analysis in periodic structures. J Sound Vib 2018b;425:170–88. [19] Li JP, Chen W, Qin QH, Fu ZJ. A modified multilevel algorithm for large-scale scientific and engineering computing. Comput Math Appl 2019b. doi:10.1016/j.camwa.2018.12.012. [20] Lin J, Zhang CZ, Sun LL, Lu J. Simulation of seismic wave scattering by embedded cavities in an elastic half-plane using the novel singular boundary method. Adv Appl Math Mech 2018;10:322–42. [21] Chen LC, Liu X, Li XL. The boundary element-free method for 2d interior and exterior Helmholtz problems. Comput Math Appl 2019;77:846–64. [22] Chen LC, Li XL. Boundary element-free methods for exterior acoustic problems with arbitrary and high wavenumbers. Appl Math Model 2019. in press [23] Tang ZC, Fu ZJ, Zheng DJ, Huang JD. Singular boundary method to simulate scattering of SH wave by the canyon topography. Adv Appl Math Mech 2018;10:912–24. [24] Fu ZJ, Yang LW, Zhu HQ, et al. A semi-analytical collocation Trefftz scheme for solving multi-term time fractional diffusion-wave equations. Eng Anal Bound Elem 2019;98:137–46. [25] Li XL. Error estimates for the moving least-square approximation and the element-free Galerkin method in n-dimensional spaces. Appl Numer Math 2016;99:77–97. [26] Liu GR. Meshfree methods: moving beyond the finite element method. second ed. Boca Raton: CRC Press; 2009. [27] Cheng YM. Meshless methods. Beijing: Science press; 2015. [28] Li XL, Dong HY. The element-free Galerkin method for the nonlinear p-Laplacian equation. Comput Math Appl 2018;757:2549–60. [29] Zhang T, Li XL. A generalized element-free Galerkin method for stokes problem. Comput Math Appl 2018;75:3127–38. [30] Li XL. Three-dimensional complex variable element-free Galerkin method. Appl Math Model 2018;63:148–71. [31] Liu GR, Gu YT. A point interpolation method for two-dimensional solids. Int J Numer Methods Eng 2001;50:937–51. [32] Dehghan M, Abbaszadeh M, Mohebbi A. The use of element free Galerkin method based on moving Kriging and radial point interpolation techniques for solving some types of turing models. Eng Anal Bound Elem 2016;62:93–111. [33] Dehghan M, Abbaszadeh M. A reduced proper orthogonal decomposition (POD) element free Galerkin (POD-EFG) method to simulate two-dimensional solute transport problems and error estimate. Appl Numer Math 2018;126:92–112. [34] Dehghan M, Ghesmati A. Numerical simulation of two-dimensional sine-gordon solitons via a local weak meshless technique based on the radial point interpolation method (RPIM). Comput Phys Commun 2010;181:772–86. [35] Dehghan M, Haghjoo-Saniji M. The local radial point interpolation meshless method for solving maxwell equations. Eng Comput 2017;33:897–918. [36] Gu YT, Liu GR. A boundary point interpolation method for stress analysis of solids. Comput Mech 2002;28:47–54. [37] Li XL, Zhu JL, Zhang SG. A meshless method based on boundary integral equations and radial basis functions for biharmonic-type problems. Appl Math Model 2011;35:737–51. [38] Ren YL, Li XL. A meshfree method for signorini problems using boundary integral equations. Math Probl Eng 2014;2014. Article ID 490127 [39] Hassanzadeh M, Tohidvand HR, Hajialilue-Bonab M, Javadi AA. Scaled boundary point interpolation method for seismic soil-tunnel interaction analysis. Comput Geotech 2018;101:208–16. [40] Telukunta S, Mukherjee S. An extended boundary node method for modeling normal derivative discontinuities in potential theory across edges and corners. Eng Anal Bound Elem 2004;28:1099–110. [41] Li XL, Zhu JL, Zhang SG. A hybrid radial boundary node method based on radial basis point interpolation. Eng Anal Bound Elem 2009;33:1273–83. [42] Zhang JM, Qin XY, Han X, Li GY. A boundary face method for potential problems in three dimensions. Int J Numer Meth Eng 2009;80:320–37. [43] Gao XW. An effective method for numerical evaluation of 2d and 3d high order singular boundary integrals. Comput Methods Appl Mech Eng 2010;199:2856–64. [44] Shen Q, Chen HR, Li S. Adaptive integration method based on sub-division technique for nearly singular integrals in near-field acoustics boundary element analysis. J Low Freq Noise Vib Act Contr 2014;33:27–45. [45] Li JP, Chen W, Fu ZJ, Qin QH. A regularized approach evaluating the near-boundary and boundary solutions for three-dimensional Helmholtz equation with wideband wavenumbers. Appl Math Lett 2019c;91:55–60. [46] Szabados J, Vértesi P. Interpolation of functions. Singapore: World Scientific; 1990. [47] Burden RL, Faires JD. Numerical analysis. Ninth ed. Boston: Brooks/Cole Cengage Learning; 2010. [48] Li XL, Dong HY. Analysis of the element-free Galerkin method for signorini problems. Appl Math Comput 2019;346:41–56. [49] Zhang T, Li XL. Meshless analysis of Darcy flow with a variational multiscale interpolating element-free Galerkin method. Eng Anal Bound Elem 2019;100:237– 245.

Acknowledgments This work was supported by the Science and Technology Research Program of Chongqing Municipal Education Commission, China (Grant No. KJZD-M201800501), the Chongqing Research Program of Basic Research and Frontier Technology, China (Grant Nos. cstc2015jcyjBX0083, cstc2018jcyjAX0266), and the National Natural Science Foundation of China (Grant No. 11471063). References [1] Ihlenburg F. Finite element analysis of acoustic scattering. New York: Springer; 1998. [2] Zhu JL, Yuan ZQ. Boundary element analysis. Beijing: Science press; 2009. [3] Kress R. Boundary integral equations in time-harmonic acoustic scattering. Math Comput Model 1991;15:229–43. [4] Ma JJ, Zhu JL, Li MJ. The Galerkin boundary element method for exterior problems of 2-d Helmholtz equation with arbitrary wavenumber. Eng Anal Bound Elem 2010;34:1058–63. [5] Wenterodt C, von Estorff O. Dispersion analysis of the meshfree radial point interpolation method for the Helmholtz equation. Int J Numer Methods Eng 2009;77:1670–89. [6] Tadeu A, Stanak P, Sladek J, Sladek V. Coupled BEM-MLPG acoustic analysis for non-homogeneous media. Eng Anal Bound Elem 2014;44:161–9. [7] Chai YB, Li W, Gong ZX, Li TY. Hybrid smoothed finite element method for two-dimensional underwater acoustic scattering problems. Ocean Eng 2016;116:129–41. [8] Chai YB, Gong ZX, Li W, Li TY, Zhang QF, Zou ZH, Sun YB. Application of smoothed finite element method to two-dimensional exterior problems of acoustic radiation. Int J Comput Methods 2018;15:1850029. [9] Wu SW, Yang X. A coupled interpolating meshfree method for computing sound radiation in infinite domain. Int J Numer Methods Eng 2018;113:1466–87. [10] Young DL, Chen KH, Lee CW. Singular meshless method using double layer potentials for exterior acoustics. J Acoust Soc Am 2006;119:96–107. [11] Jiang XR, Chen W, Chen CS. Fast multipole accelerated boundary knot method for inhomogeneous Helmholtz problems. Eng Anal Bound Elem 2013;37:1239–43. [12] Fu ZJ, Xi Q, Chen W, Cheng A. A boundary-type meshless solver for transient heat conduction analysis of slender functionally graded materials with exponential variations. Comput Math Appl 2018a;76:760–73. [13] Lin J, Chen CS, Liu CS, Lu J. Fast simulation of multi-dimensional wave problems by the sparse scheme of the method of fundamental solutions. Comput Math Appl 2016;72:555–67. [14] Wang FJ, Liu CS, Qu W. Optimal sources in the MFS by minimizing a new merit function: energy gap functional. Appl Math Lett 2018;86:229–35. [15] Li JP, Fu ZJ, Chen W, Liu XT. A dual-level method of fundamental solutions in conjunction with kernel-independent fast multipole method for large-scale isotropic heat conduction problems. Adv Appl Math Mech 2019a;11:501–17. [16] Fu ZJ, Chen W, Gu Y. Burton-miller-type singular boundary method for acoustic radiation and scattering. J Sound Vib 2014;333:3776–93.

21