A new fast multipole boundary element method for two dimensional acoustic problems

Comput. Methods Appl. Mech. Engrg. 200 (2011) 1333–1340

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

A new fast multipole boundary element method for two dimensional acoustic problems Shande Li, Qibai Huang ⇑ State Key Laboratory of Digital Manufacturing Equipment and Technology, Huazhong University of Science and Technology, Wuhan 430074, PR China

a r t i c l e

i n f o

Article history: Received 28 June 2010 Received in revised form 2 November 2010 Accepted 3 November 2010 Available online 12 November 2010 Keywords: Fast multipole method Boundary element method Burton–Miller formulation Helmholtz equation Approximate inverse preconditioner Acoustic problems

a b s t r a c t A new fast multipole boundary element method (BEM) is presented in this paper for solving large-scale two dimensional (2D) acoustic problems based on the improved Burton–Miller formulation. This algorithm has several important improvements. The fast multipole BEM employs the improved Burton–Miller formulation, and successfully overcomes the non-uniqueness difficulty associated with the conventional BEM for exterior acoustic problems. The improved Burton–Miller formulation contains only weakly singular integrals, and avoids the numerical difficulties associated to the evaluation of the hypersingular integral, it leads to the numerical implementations more efficient and straightforward. Furthermore, the fast multipole method (FMM) and the approximate inverse preconditioned generalized minimum residual method (GMRES) iterative solver are adopted to greatly improve the overall computational efficiency. The numerical examples with Neumann boundary conditions are presented that clearly demonstrate the accuracy and efficiency of the developed fast multipole BEM for solving large-scale 2D acoustic problems in a wide range of frequencies. Ó 2010 Elsevier B.V. All rights reserved.

1. Introduction Boundary element method (BEM) has been widely used to solve acoustic problems, since it offers an excellent accuracy and easy mesh generation. For exterior acoustic problems, the radiation condition is automatically satisfied at infinity. However, the conventional BEM has some drawbacks. One drawback is that for exterior acoustic problems the conventional integral equation suffers from the problem of interior resonances [1]. These frequencies are due only to the numerical method, which has no unique solution at some characteristic frequencies of the corresponding interior problem. Two major methods [2,3] have been proposed to remove the non-uniqueness difficulty. In the first method, a combined Helmholtz integral equation formulation (CHIEF) [2] successfully removed the non-uniqueness difficulty by adding some additional Helmholtz integral relations in the interior domain. This resulted in an over-determined system of equations, which was then solved using a least-squares technique. However, selecting the optimum number and suitable positions of the interior points may become difficult for the high-frequency problems or complicated structures [4]. Another well-known formulation to overcome ⇑ Corresponding author. Address: School of Mechanical Science and Engineering, Huazhong University of Science and Technology, Wuhan 430074, PR China. Tel.: +86 27 87557664; fax: +86 27 87544175. E-mail addresses: [email protected], [email protected] (S. Li), qbhuang@ mail.hust.edu.cn (Q. Huang). 0045-7825/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2010.11.005

the non-uniqueness difficulty at characteristic frequencies was proposed by Burton and Miller [3]. This method consisted of a linear combination of the Helmholtz equation and its normal derivative equation. The major difficulty of the Burton–Miller formulation is that it includes a hypersingular integral. The hypersingular integral evaluation is an expensive task, a weakly singular form of the hypersingular integral is necessary. Various methods may effectively implement hypersingular problem, including the concept of Hadamard finite-part integral [5], the high-order Galerkin method [6], singularity subtraction technique and properties from associated Laplace equations [7,8], and so on. In contrast to other formulations, the Burton–Miller method works very reliable and offers efficiency and robustness to the exterior acoustic problems of arbitrary shape. The second drawback of the conventional BEM is that discretization of the integral equation leads to a dense, non-symmetric system matrix which is expensive to store and solve. For an acoustic problem with N unknowns, a direct solution of such linear system requires O(N3) operations and O(N2) memory storage. As a result, the conventional BEM is prohibitively expensive for largescale engineering problems. Improving the overall computational efficiency is the main task in the implementation of the conventional BEM. One possibility is to solve the matrix equation with iterative methods, such as the generalized minimum residual method (GMRES) [9] and the conjugate gradient method (CGM) [10]. Iterative methods can reduce the cost to O(NiterN2) operations, where Niter is the number of iterations required, and the O(N2) each

1334

S. Li, Q. Huang / Comput. Methods Appl. Mech. Engrg. 200 (2011) 1333–1340

iteration cost arises from the dense matrix–vector product. Still the conventional BEM was not used for solving large-scale problems, until the recent development of the fast multipole method (FMM) [11,12], drastically changed the circumstances around the BEM. The FMM was initially introduced by Rokhlin and Greengard [11,12] as a fast solution method for boundary integral equations, the operations and memory requirement can be drastically reduced from O(N2) to O(N). The key idea of the FMM is to combine the effect of sources far away from a collocation point in a far-field term using the multipole expansion, that is, cell-to-cell interactions of far-field translations instead of element-to-element interactions. The introduction of the FMM in the conventional BEM has generated enormous interests, and the fast multipole BEM (FMBEM) has been intensively applied to really solve large-scale problems of various scientific fields, including those arising from the Helmholtz, Maxwell and elasticity equations [13–19]. A comprehensive review of the fast multipole BEM can be found in Ref. [20]. Besides the FMM, another fast method is hierarchical matrices [21] is used to improve the efficiency of the BEM, this method is of pure algebraic nature and represents partitions of the original system matrix by low-rank approximations in outer-product form. A comparison of the fast multipole with hierarchical matrices for the Helmholtz integral equation is presented in Ref. [22]. The development of the fast multipole BEM for solving largescale acoustic problems is the important advantage compared with other methods. The FMM has been extended to the field of acoustics for quite some time. Rokhlin [16] first applied the FMM for solving the Helmholtz equation, and proposed diagonal form of the translation matrices for high-frequency problems [17,18]. A very detailed description of the FMM for the Helmholtz equation was presented by Epton and Dembart [19]. Darve [23] provided more precise error estimates introduced by the FMM. Gumerov and Duraiswami [24,25] presented the FMM accelerated iterative solution of the BEM based on Burton–Miller formulation for the 3D Helmholtz equation, and proposed recurrence relations to develop a general recursive method for obtaining the translation matrices, these relations provided faster translations. Chen and Chen [26] solved the 2D exterior acoustic problems and expanded the four kernels in the dual BEM into degenerate kernels that separated the field and source points by the FMM. An adaptive fast multipole BEM for 3D acoustic problems was investigated by Shen and Liu [27] where the Burton–Miller formulation was applied to overcome the non-uniqueness difficulty, the large-scale practical acoustic model with 200,000 unknowns was solved on a laptop PC. More recently, several research works have been published for extending the applicability of the fast multipole BEM for a wide range of frequencies [28,29]. New engineering applications of the fast multipole BEM for solving the half-space acoustic problems [30] and simulating the noise distribution around the 2D acoustic barrier [31] have been developed. In this paper, a new fast multipole BEM formulation is presented for 2D acoustic problems based on the improved Burton–Miller formulation. The developed algorithm has several important advantages. First, the fast multipole BEM adopts the improved Burton–Miller formulation, and successfully overcomes the non-unique difficulty associated with the conventional BEM for exterior acoustic problems. The improved form of the hypersingular integral is regularized by using the weakly singular integrals, which are obtained analytically and do not introduce any approximations, thus making the numerical implementations more efficient and straightforward. Second, the computational efficiency of the developed algorithm is further improved by adopting the FMM and the preconditioned GMRES iterative solver to solve system matrix equation. The approximate inverse preconditioner is very efficient and results in a large reduction of required iteration

steps. Thus, the overall solution efficiency of the developed algorithm is further improved. Several numerical examples of Neumann type problems are presented to study the accuracy and efficiency of the developed fast multipole BEM. These results clearly demonstrate the advantages of the developed fast multipole BEM for solving large-scale 2D acoustic problems in a wide range of frequencies. This paper is organized as follows: the basic BEM formulation is reviewed in Section 2. The formulation and algorithm of the fast multipole BEM for 2D acoustic problem are presented in Section 3. Numerical examples demonstrate the accuracy and efficiency of the fast multipole BEM for large-scale acoustic problems in Section 4. Section 5 concludes the paper with further discussions. 2. The improved Burton–Miller formulation for acoustic problems The problem under consideration in this paper is the solution of the Helmholtz equation in the domain E (which can be either finite or infinite). To be precise, we consider propagation of timeharmonic acoustic waves in a homogeneous isotropic acoustic medium is described by the Helmholtz equation

r2 /ðxÞ þ k2 /ðxÞ ¼ 0;

x 2 E:

ð1Þ

For the exterior acoustic problems it is necessary to introduce a condition at infinity. This ensures the physical requirement that all scattered and radiated waves are outgoing. This is termed the Sommerfeld radiation condition 1

lim r 2

r!1



 @/  ik/ ¼ 0; @r

ð2Þ

where / is the total acoustic wave (velocity potential or acoustic pressure), k = x/c the wave number, x the angular frequency, and c the wave speed in the acoustic medium E, r the distance from a pffiffiffiffiffiffiffi fixed origin to a general field point and i ¼ 1. The integral representation of the solution to the Helmholtz equation is

CðxÞ/ðxÞ ¼

 Z  @Gk ðx; yÞ  /ðyÞ þ Gk ðx; yÞqðyÞ dSy þ /in ðxÞ; @ny S

ð3Þ

where x is the collocation point, y is the source point, q(y) is the normal derivative of / at the point y, Sy denotes a subset of surface S, /in(x) is a prescribed incident wave (for scattering problems only), the coefficient C(x) = 1, 1/2 or 0 when the collocation point x is in the exterior region E (acoustic medium), on the boundary S (if it is smooth) or in the interior region B(a body or scatterer), respectively, ny is the outward normal at y, the fundamental solution Gk for 2D acoustic problems is given by

Gk ðx; yÞ ¼

i ð1Þ h ðkrÞ; 4 0 ð1Þ

with r ¼ jx  yj

ð4Þ

in which h0 ðÞ denotes the Hankel function of the first kind [32]. Eq. (3) with x 2 S is the commonly used form of the conventional BEM for acoustic wave problems. If the normal derivative is given on the boundary S, the acoustic pressure can be computed at any point in E by using Eq. (3). However, the solution of Eq. (3) does not possess a unique solution at certain characteristic frequencies of the corresponding interior problems. A number of different methods [2–8] have been proposed for overcoming this non-unique problem. One of the most effective and robust method is the Burton–Miller method [3] and for this reason, this is the method we shall use here. Burton–Miller formulation is a linear combination of the Helmholtz integral Eq. (3) and its normal derivative with an appropriate coupling coefficient b, it is can be written as

1335

S. Li, Q. Huang / Comput. Methods Appl. Mech. Engrg. 200 (2011) 1333–1340

Z S

Z

@Gk ðx; yÞ /ðyÞdSy  Gk ðx; yÞqðyÞdSy @ny S Z Z @G2k ðx; yÞ @Gk ðx; yÞ þb /ðyÞdSy  b qðyÞdSy @nx S @ny @nx S 1 b @/ ¼  /ðxÞ  qðxÞ þ /in þ b in x 2 S; 2 2 @nx

resulting values are sometimes called influence coefficient. They will be called Hij and Gij, i.e.

Hij ¼ b ð5Þ

Eq. (5) has a unique solution for all real and positive k if the imaginary part of the coupling coefficient b is non-zero. We choose the parameter b [33] as

 b¼

i;

k 6 1;

ð6Þ

i=k; k > 1:

S

@G2k ðx; yÞ /ðyÞdSy ¼ @ny @nx

" # @ 2 Gk ðx; yÞ @ 2 G0 ðx; yÞ  dSy @nx @ny @nx @ny S Z Z @G0 ðx; yÞ þ r/ðxÞ  ny dSy þ ½/ðyÞ @nx S S Z

/ðyÞ

 /ðxÞ  r/ðxÞ  ðy  xÞ

Z S

@Gk ðx; yÞ qðyÞdSy ¼ @nx

 @Gk ðx; yÞ @G0 ðx; yÞ qðyÞdSy þ @nx @ny S   Z @G0 ðx; yÞ @/ðxÞ dSy :  qðyÞ  @ny @nX S

ð7Þ

Z 

ð8Þ

Z

Sj

Sj

@Gk ðx; yÞ 1 dSy þ dij ; @ny 2

ð10Þ

b Gk ðx; yÞdSy  dij ; 2

ð11Þ

Hij /j ¼

j¼1

N X

Gij qj þ fbi g;

ð12Þ

j¼1

with i = 1, 2, . . ., N. In the conventional BEM, a standard linear system of equations is formed by applying the boundary conditions at each node and moving the unknown terms to the left-hand side and the known terms to the right-hand side (e.g., for Neumann type problems, the pressure / is unknown values and the normal derivative q is specified) in Eq. (12), leads to the following matrix equation

ð13Þ

where A is a square system matrix, x the unknown vector, b the known vector. Obviously, the construction of matrix A requires O(N2) operations, as well as O(N2) memory storage since A is a dense and non-symmetric matrix. That is why the conventional BEM is in general slow and inefficient for large-scale problems. However, using the fast multipole method will reduce the operations from O(N2) to O(N), and memory requirement from O(N2) to O(N). In the following, the basic fast multipole BEM formulations for the 2D Helmholtz equation are reviewed. 3. The fast multipole BEM 3.1. The fast multipole method formulation

where G0 is the free-space Green’s function for the Laplace equation. A 2D form of G0 can be expressed as G0 = log(1/r)/(2p), r/(x) denotes the domain gradient of / at the point x. Substituting Eq. (7) and (8) into (5), we obtain the improved weakly singular boundary integral equations for acoustic problems, in which has not only a unique solution (like Burton–Miller method) but also no hypersingularity (unlike Burton–Miller method). The improved Burton–Miller formulation contains only weakly singular integrals and is evaluated analytically. The analytical integration has the benefit of accuracy and efficiency and is well suited for integration with the fast multipole BEM. The boundary S is assumed to be divided into N elements (e.g. using linear element in this study). Eq. (5) can be discretized for a given point i before applying any boundary conditions and rearranged as follows:

" Z # Z @G2k ðx; yÞ @Gk ðx; yÞ 1 b /ðyÞdSy þ /ðyÞdSy þ /i @n @n 2 @n y x y S S j j j¼1 " Z # Z N X @Gk ðx; yÞ ¼ b qðyÞdSy þ Gk ðx; yÞqðyÞdSy @nx Sj Sj j¼1

N X

b  qi þ fbi g; 2

Sj

@Gk ðx; yÞ dSy þ @nx

Z

Ax ¼ b;

@ 2 G0 ðx; yÞ dSy @nx @ny

1  r/ðxÞ  nx ; 2

Z

Sj

@G2k ðx; yÞ dSy þ @ny @nx

where d is Kronecker’s delta. The coefficients Hij and Gij are sums of the integrals on the elements around node j, respectively. The shape functions for different types of elements are given in Ref. [34]. Substituting Eqs. (10) and (11) into (9) yields N X

However, the major difficulty in Eq. (5) is that the normal derivative of the Helmholtz integral equation involves a hypersingular integral operator. The hypersingular integral evaluation is an expensive task, a weakly singular form of the hypersingular integral is necessary. Using singularity subtraction technique and properties from associated Laplace equations developed by authors recently [8], the hypersingular and singular integrals in Eq. (5) can be regularized as a weakly form

Z

Gij ¼ b

Z

In this section, the FMM is employed to solve the improved Burton–Miller formulation. Iterative solver is used in the FMM, where matrix–vector products are calculated using fast multipole expansions. Using the FMM for the BEM, the computational cost is reduced significantly to O(N) operations. The memory requirement can also be reduced to O(N) since iterative solver does not need to store the entire matrix in the memory. The appropriate theories of the expansions and translations used in the FMM for 2D Helmholtz equation are scattered in the Refs. [16,20]. 3.1.1. Multipole expansion First, the fundamental solution can be obtained from the boundary integral equations. It clearly shows the relationship between field and source points. Therefore, the fundamental solution can be separated the field points and source points into two terms. In 2D acoustic problems, letting x = (qx, hx) and y = (qy, hy) in polar coordinates, Graf’s addition theorem gives [32] ð1Þ

h0 ðkkx  ykÞ ¼

1 X

ð1Þ

hn ðkqx Þjn ðkqy Þeinðhx hy Þ ;

ð14Þ

n¼1

ð9Þ

where Sj is element j of the boundary. {bi} comes from the incident wave only for scattering problems. The point i is one of the boundary nodes. There are now two type of integrals to be carried out over the elements. These integrals relate the node i where the fundamental solution is acting on any other node j. Because of this their

where qx > qy. Thus, the fundamental solution can be expressed as a multipole expansion around the center yc near y (Fig. 1)

8 1 P ! ! i > > Sn ðyc xÞRn ðyc yÞ < Gk ðx; yÞ ¼ 4 n¼1

! 1 P ! > > k ðx;yÞ : @G@n ¼ 4i Sn ðyc xÞ @Rn@nðyyc yÞ y

n¼1

;

!

!

yc x > yc y;

ð15Þ

1336

S. Li, Q. Huang / Comput. Methods Appl. Mech. Engrg. 200 (2011) 1333–1340

Fig. 1. The fast multipole method: M2M, M2L and L2L translations.

where the functions !

! Sn ðyc xÞ

and

! Rn ðyc yÞ

!

n ð1Þ

Sn ðyc xÞ ¼ i hn ðkj yc x jÞeinh ; !

n

!

!

are defined as

ð16Þ inh

Rn ðyc yÞ ¼ ð1Þ jn ðkj yc y jÞe ; !

ð17Þ ð1Þ

!

in these formulae jn ðkj yc y jÞ and hn ðkj yc x jÞ are the spherical Bessel function and spherical Hankel function of the first kind, respectively. h is the polar angle of yc respect to x and y, respectively. In Eq. (5), there are a total of four integrals that need to be evaluated, we use the first and second integrals for the formulation, the third and fourth integrals have a similar procedure as used for the first and second integrals. For integrals on a subset Sy which is far ! ! away from a collocation point xi ðj yc x j > j yc y jÞ, the first and second integrals in Eq. (5) can be written as the multipole expansions around yc as follows:

 Z  @Gk ðxi ; yÞ /ðyÞ  Gk ðxi ; yÞqðyÞ dSy ; @ny Sy X ! i 1 Sn ðyc xi ÞMn ðyc Þ; fornode i ¼ 1; 2; . . . ; N; ¼ 4 n¼1

ð18Þ

where Mn(yc) is the multipole moment centered at yc defined as

M n ðyc Þ ¼

! ! ! @Rn ðyc yÞ /ðyÞ  Rn ðyc yÞqðyÞ dSy : @ny

Z Sy

ð19Þ

!

!

!

where j xc y j > j xc x j and j yc x j > j yc y j. This translation is called multipole-to-local (M2L) translation. 3.1.4. Local moment to local moment translation The coefficients of the local expansion are translated according to the following formula as the center of the local expansion is shift from xc to xc0 using local-to-local (L2L) translation

Ln ðxc0 Þ ¼

1 X

 ! Rnm xc xc0 Lm ðxc Þ;

ð23Þ

m¼1 !

!

where j xc0 y j > j xc0 x j. Another the third and fourth integrals in Eq. (5) can also be shown that the moments, M2M, M2L, and L2L translations have a similar procedure as used for the first and second integrals. M2M, M2L, L2L translations are showed in Fig. 1. It is noted that the infinite sums in Eq. (15), (18), (21) and (23) must be truncated with suitable terms for practical computation [16,17]. Finally, evaluation of the integrals in Eq. (5) is the sum of two contributions. For a group of source points which are far away from the collocation point xi, the contributions can be evaluated from local expansion coefficients obtained from Eq. (23). For source points which are adjacent to the collocation point xi, the contributions are evaluated directly from Eq. (5) as in the conventional BEM. 3.2. Numerical implementation of the fast multipole BEM

3.1.2. Multipole moment to multipole moment translation The multipole moments are translated according to the following formula:

M n ðyc0 Þ ¼

1 X

!

Rnm ðyc0 yc ÞMm ðyc Þ;

ð20Þ

m¼1

Eq. (20) is used to translate the multipole moment as the center of the multipole expansion is shifted from yc to yc0 , and this translation is called multipole-to-multipole (M2M) translation. 3.1.3. Local expansion In the local expansion, the first and second integrals in Eq. (5) are expanded with respect to the collocation point x around a selected point xc (Fig. 1)

 Z  @Gk ðxi ; yÞ /ðyÞ  Gk ðxi ; yÞqðyÞ dSy @ny Sy 1 ! i X ¼ ð1Þn Rn ðxc xi ÞLn ðxc Þ; 4 n¼1

ð21Þ

where Ln(xc) is called the local expansion coefficients centered at xc, can be obtained from the multipole moments centered at yc as

Ln ðxc Þ ¼

1 X m¼1

!

Snm ðyc xc ÞM m ðyc Þ;

ð22Þ

In the numerical implementation of the fast multipole BEM, firstly the boundaries of model are discretized. Secondly a hierarchical tree structure is constructed. Then an iterative process is executed to solve the equation system. A very detailed description of the algorithm for the fast multipole BEM was given in Refs. [13,20,35]. 3.2.1. Discretization The fast multipole BEM has the same discretizations as the conventional BEM. The boundary S is discretized using boundary elements. 3.2.2. Tree construction A quad-tree structure is constructed hierarchically to store both the multipole and local moments. The root of the tree represents a square that contains the boundary S, we call this square the cell of level 0 (Fig. 2). This parent cell is divided into four equal child cells at level 1 whose the edge length is half of the parent cell. Each child cell is continuously divided in the same way until every cell contains at most some predefined number of elements. A childless cell is called a leaf. An example is shown in Fig. 2, where each leaf contains at most one element at the finest level 4. Furthermore, the concepts of adjacent cells, interaction list and far cells should be defined. We say that two cells at the same level l

S. Li, Q. Huang / Comput. Methods Appl. Mech. Engrg. 200 (2011) 1333–1340

Cell C

Cells in interaction list of C (M2L) Far cells to C (L2L) Adjacent cells to C (direct BEM)

Fig. 2. A hierarchical cell structure covering all the boundary elements.

are adjacent cells if they share at least one common vertex (Fig. 2). The interaction list of cell C at level l is those cells whose parent cells are adjacent to C’s parent at level l  1, but themselves are not adjacent to C at level l (Fig. 2). Far cells of C at level l are those cells whose parent cells are not adjacent to the parent cell of C (Fig. 2). 3.2.3. Iterative process Step 1: Upward stage. Computation of the multipole moments. For a leaf, the multipole moments are obtained from all the elements in it using Eq. (19). For a non-leaf cell, the multipole moments are shifted from its child cells using M2M translation (Eq. (20)). This procedure is repeated upward to level 2. Step 2: Downward stage. Computation of the local expansion coefficients. After multipole moments are obtained for all nodes in the quad-tree structure, a down stage is performed in order to obtain local expansion coefficients for each node. For each cell C, the local expansion coefficients are the sum of two parts. One part is obtained from the multipole moments of C’s interaction list using M2L translation (Eq. (22)). The other part is obtained from the local moments of C’s parent cell using L2L translation (Eq. (23)). Starting from level 2, this procedure is repeated downward to the finest level. Step 3: Evaluation of the integrals in Eq. (5). For a collocation point xi is in leaf C, evaluation of the integrals in Eq. (5) are the sum of two parts. The contributions from elements in leaf C and its adjacent cells are evaluated directly in the way as conventional BEM (Fig. 2). The contributions from all other elements are obtained from the local moments of C using Eq. (21). Step 4. Update the iterative vector. The iterative vector in the system Ax = b corresponding to Eq. (5) is updated, and next iteration is started from Step 1 until the desired accuracy is achieved.

to the inverse of system matrix A for which kAP1  Ik is small in some norm. Then P1 is used as a preconditioner. Let P1 = [m1, . . ., mn] and I = [e1, . . ., en]. Using the usual Frobenius norm, P minP1 kAP1  Ik2F ¼ nj¼1 minmj kAmj  ej k22 , where we let 1 define the sparsity pattern of P1. The pattern 1 is constructed by adjacent elements of a collocation point. For each column j, minmj kAmj  ej k22 is a least squares problem for an n  g system where g is the number of columns of A involved via pattern 1. We can solve the approximate problem minmj kAj mj  ej k22 where Aj is the g  g principal submatrix of A induced by the columns of A specified by 1. For approximate inverse preconditioner of GMRES, the system Eq. (13) is simply replaced by

P1 Ax ¼ P1 b:

ð24Þ

For a problem with N equations, the computational complexity of the approximate inverse preconditioner isO(N). 4. Numerical examples We present several examples with Neumann boundary conditions to demonstrate the accuracy and efficiency of the fast multipole BEM for solving 2D acoustic problems. The developed algorithm has been implemented in Fortran 90, and has been tested on a laptop PC with an Intel 2.2 GHz CPU and 2 GB RAM. In all the examples, we use linear elements to discretize boundaries, and truncate the infinite sums in Eq. (15), (18), (21) and (23) with 10 terms. The maximum number of boundary elements in a leaf is set to 20. The conventional BEM code uses the iterative solver for solving the linear system. 4.1. Verification and performance The first numerical example is used to validate the accuracy of the developed algorithm for acoustic problems. The rigid cylinder where the normal derivative @//@n = 0 on the surface S is impinged upon by an incident wave of unit amplitude along the negative xaxis as shown in Fig. 3, the radius a of the cylinder is 1 m. The available analytical solution [37] of the pressure / can be used to check the fast multipole BEM results. The whole cylinder is considered for modeling the problem with 3144 elements. The iterative solver (GMRES) will terminate the iterations when the tolerance for convergence of the solution is less than 105. As shown in Fig. 4, the acoustic pressure j/j at computational point (0, 2) is plotted versus the nondimensionalized wave numbers ka = 0–10 with

3.3. Approximate inverse preconditioner Preconditioners for the fast multipole BEM are crucial for its convergence and computing efficiency. We remark that a potentially useful application of approximate inverse preconditioner is to combine it with the FMM where the available near-field elements provide sufficient information for preconditioner construction. The application of the sparse approximate inverse preconditioned GMRES to the BEM system for solving acoustic problems is presented by Chen and Harris [36]. The basic idea of the approximate inverse method is operator splitting. Briefly the objective is to construct approximations,P1,

1337

Fig. 3. Plane wave scattering from a rigid cylinder with radius a.

1338

S. Li, Q. Huang / Comput. Methods Appl. Mech. Engrg. 200 (2011) 1333–1340 100

Analytical solution Fast multipole BEM Conventional BEM

1.6

Approximate inverse preconditioned GMRES Unpreconditioned GMRES

10-2 1.4

|φ |

Residual

1.2

1

10-4

10-6 0.8

10-8 0.6

0.4

0

1

2

3

4

5

6

7

8

9

10

10-10 0

5

10

15

ka

20

25

30

35

Iteration

Fig. 4. Pressure from an acoustic scattering cylinder at the computational point.

Fig. 6. Convergence of the approximate inverse preconditioned GMRES.

100 frequency steps. The fictitious frequencies for the conventional BEM can be identified clearly at ka = 2.41, 5.53, and 8.67, near which the conventional BEM results deviate substantially from the analytical solution (a clear indication of the non-uniqueness of the conventional BEM solution). However, the fast multipole BEM results provide very good agreement with the analytical solutions with a relative error of less than 0.6% throughout all wave numbers. Next we study the efficiency of the fast multipole BEM by the same cylinder. The total CPU time used for solving scattering sphere as the number of elements increase from 632 to 12,568 are plotted in Fig. 5, which shows significant advantage of the fast multipole BEM compared with the conventional BEM. Two curves for the fast multipole BEM and conventional BEM show the order O(N) and O(N2) computational efficiencies, respectively. This numerical example clearly demonstrates the developed algorithm is very efficient for solving large-scale acoustic problems. Then we further explore the efficiency of the approximate inverse preconditioner for the fast multipole BEM. The whole cylinder is considered for modeling the problem with 3144 elements. The ka is taken to be 1. Fig. 6 shows the convergence of the approximate inverse preconditioned GMRES and unpreconditioned GMRES. As shown in Fig. 6, the number of iterations reduces dra-

matically with the use of the approximate inverse preconditioner, fast convergence provides finally 4 times faster time (computational cost) to solve the same problem. This numerical example shows fast convergence rate of the developed algorithm with the efficient preconditioner. Finally we test the convergence rate of the approximate inverse preconditioner for a wide range of frequencies. As shown in Fig. 7, the convergences of the approximate inverse preconditioned GMRES and unpreconditioned GMRES are plotted versus the frequencies from 50 to 500 Hz. The approximate inverse preconditioned GMRES shows little influence as frequency increases. This numerical example further demonstrates the developed algorithm is very efficient to solve large-scale high-frequency acoustic problems.

4500

4.2. Large-scale acoustic problem Many problems in acoustics require computations for largescale complex models, which include multiple scatterers, acoustic barrier, room acoustics, sound propagation in dispersed media, etc. We test the developed fast multipole BEM for solving more realistic example, and further explore the efficiency and applicability of the developed algorithm for large-scale acoustic problems.

Fast multipole BEM Conventional BEM

4000

90

Approximate inverse preconditioned GMRES Unpreconditioned GMRES

80 70

3000 60

2500

Iteration

Total CPU time (s)

3500

2000

50 40

1500 30

1000 20

500 10

0

0

2000

4000

6000

8000

10000

12000

14000

Number of elements

0

0

50

100

150

200

250

300

350

400

450

500

Frequency (Hz) Fig. 5. Comparison of the CPU time used by the fast multipole BEM and conventional BEM.

Fig. 7. Frequency dependency of the approximate inverse preconditioned GMRES.

S. Li, Q. Huang / Comput. Methods Appl. Mech. Engrg. 200 (2011) 1333–1340

Fig. 8. Exterior multiple scattering for ten cylinders involving 123,000 elements when k = 1.

1339

the weakly singular integrals, which are obtained analytically and do not introduce any approximations, thus making the numerical implementations more efficient and straightforward. Furthermore, the computational efficiency of the developed algorithm is further improved by adopting the FMM and the approximate inverse preconditioned GMRES iterative solver to solve system matrix equation. The approximate inverse preconditioner proves more efficient and results in a large reduction of required iteration steps. The numerical examples including a model with 123,000 elements are solved successfully on a laptop PC. These examples clearly demonstrate the accuracy and efficiency of the developed fast multipole BEM for solving large-scale 2D acoustic problems in a wide frequency range, and show the great potential of engineering application. Further works need to further improve the developed fast multipole BEM. Firstly, one problem is valuably discussed. All the integrals in the improved Burton–Miller formulation are at worst weakly singular and so can be readily evaluated using higher-order elements. Thus the improved Burton–Miller formulation opens up new possibility of using higher-order elements to replace constant or linear elements for computational models in the developed algorithm. There will be a considerable improvement in the computational accuracy. Research works are in progress and focus on whether using higher-order elements may evidently reduce the efficiency of the developed algorithm for large-scale problems. Furthermore, using a new adaptive tree structures proposed by Bapat and Liu [38] may further improve and accelerate the present algorithm. These further works are in progress. The developed fast multipole BEM can be readily extended to solve 3D acoustic problems, and can also be extended to other applications, including those arising from the Maxwell and elasticity equations. Acknowledgments Financial supports from the National Nature Science Foundation of China (No. 50075029) and Research Fund for the Doctoral Program of Higher Education of China (No. 20070487403) are gratefully acknowledged. References

Fig. 9. Exterior multiple scattering for ten cylinders involving 123,000 elements when k = 10.

We study acoustic scattering from multiple scatterers. Ten randomly distributed rigid cylinders where the normal derivative @/ / @n = 0 on the surface S are modeled with 123,000 elements, the radius a of the cylinder is 5 m. The incoming unit plane wave travels along the negative x-axis. Figs. 8 and 9 show the pressure field when wave number k = 1 and k = 10, respectively. All the fast multipole BEM results converge in about 40 iterations. The total CPU time used is less than 3000 s. This example clearly demonstrates the potentially useful application of the developed algorithm for solving large-scale engineering acoustic problems. 5. Conclusion A new fast multipole BEM for solving 2D acoustic problems is presented in this paper based on the improved Burton–Miller formulation. Some improvements can be obtained. The fast multipole BEM employs the improved Burton–Miller formulation, and successfully overcomes the non-uniqueness difficulty associated with the conventional BEM for exterior acoustic problems. The improved form of the hypersingular integral is regularized by using

[1] D. Colton, R. Kress, Integral Equation Methods in Scattering Theory, Wiley, New York, 1983. [2] H.A. Schenck, Improved integral formulation for acoustic radiation problems, J. Acoust. Soc. Am. 44 (1968) 41–58. [3] A.J. Burton, G.F. Miller, The application of integral equation methods to the numerical solution of some exterior boundary value problems, Proc. R. Soc. Lond. A 323 (1971) 201–210. [4] I.L. Chen, J.T. Chen, M.T. Liang, Analytical study and numerical experiments for radiation and scattering problems using the CHIEF method, J. Sound Vib. 248 (2001) 809–828. [5] C.C. Chien, H. Rajiyah, S.N. Atluri, An effective method for solving the hypersingular integral equations in 3D acoustics, J. Acoust. Soc. Am. 88 (1990) 918–937. [6] P.J. Harris, K. Chen, On efficient preconditioners for iterative solution of a Galerkin boundary element equations for the three dimensional exterior Helmholtz problem, J. Comput. Appl. Math. 156 (2003) 303–318. [7] Y.J. Liu, F.J. Rizzo, A weakly singular form of the hypersingular boundary integral equation applied to 3D acoustic wave problems, Comput. Method Appl. Mech. Engrg. 96 (1992) 271–287. [8] S. Li, Q. Huang, An improved form of the hypersingular boundary integral equation for exterior acoustic problems, Engrg. Anal. Bound. Elem. 34 (2010) 189–195. [9] Y. Saad, M.H. Schultz, GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Statist. Comput. 7 (1986) 856–869. [10] R.W. Freund, Conjugate gradient-type methods for linear system with complex symmetric matrices, SIAM J. Sci. Statist. Comput. 13 (1992) 425–448. [11] V. Rokhlin, Rapid solution of integral equations of classical potential theory, J. Comput. Phys. 60 (1985) 187–207. [12] L. Greengard, V. Rokhlin, A fast algorithm for particle simulations, J. Comput. Phys. 73 (1987) 325–348.

1340

S. Li, Q. Huang / Comput. Methods Appl. Mech. Engrg. 200 (2011) 1333–1340

[13] K. Yoshida, N. Nishimura, S. Kobayashi, Application of new fast multipole boundary integral equation method to crack problems in 3D, Engrg. Anal. Bound. Elem. 25 (2001) 239–247. [14] W.C. Chew, H.Y. Chao, T.J. Cui, et al., Fast integral equation solvers in computational electromagnetics of complex structures, Engrg. Anal. Bound. Elem. 27 (2003) 803–823. [15] Y.J. Liu, A new fast multipole boundary element method for solving large-scale two-dimensional elastostatic problems, Int. J. Numer. Methods Engrg. 65 (2006) 863–881. [16] V. Rokhlin, Rapid solution of integral equations of scattering theory in two dimensions, J. Comput. Phys. 86 (1990) 414–439. [17] R. Coifman, V. Rokhlin, S. Wandzura, The fast multipole method for the wave equation: a pedestrian prescription, IEEE Ant. Propag. Mag. 35 (1993) 7–12. [18] V. Rokhlin, Diagonal forms of translation operators for the Helmholtz equation in three dimensions, Appl. Comput. Harmon. Anal. 1 (1993) 82–93. [19] M. Epton, B. Dembart, Multipole translation theory for the three dimensional Laplace and Helmholtz equations, SIAM J. Sci. Comput. 16 (1995) 865–897. [20] N. Nishimura, Fast multipole accelerated boundary integral equation methods, Appl. Mech. Rev. 55 (2002) 299–324. [21] M. Bebendorf, Hierarchical Matrices: A Means to Efficiently Solve Elliptic Boundary Value Problems, Springer, Berlin, Heidelberg, 2008. [22] D. Brunner, M. Junge, P. Rapp, M. Bebendorf, L. Gaul, Comparison of the fast multipole method with hierarchical matrices for the Helmholtz-BEM, CMES 58 (2010) 131–160. [23] E. Darve, The fast multipole method I: error analysis and asymptotic complexity, SIAM J. Numer. Anal. 38 (2000) 98–128. [24] N.A. Gumerov, R. Duraiswami, Recursions for the computation of multipole translation and rotation coefficients for the 3-D Helmholtz equation, SIAM J. Sci. Comput. 25 (2003) 1344–1381. [25] N.A. Gumerov, R. Duraiswami, Fast Multipole Methods for the Helmholtz Equation in Three Dimensions, Elsevier, Amsterdam, 2004.

[26] J.T. Chen, K.H. Chen, Applications of the dual integral formulation in conjunction with fast multipole method in large-scale problems for 2D exterior acoustics, Engrg. Anal. Bound. Elem. 28 (2004) 685–709. [27] L. Shen, Y.J. Liu, An adaptive fast multipole boundary element method for three-dimensional acoustic wave problems based on the Burton–Miller formulation, Comput. Mech. 40 (2007) 461–472. [28] H. Cheng, W.Y. Crutchfield, Z. Gimbutas, et al., A wideband fast multipole method for the Helmholtz equation in three dimensions, J. Comput. Phys. 216 (2006) 300–325. [29] N.A. Gumerov, R. Duraiswami, A broadband fast multipole accelerated boundary element method for the three dimensional Helmholtz equation, J. Acoust. Soc. Am. 125 (2009) 191–205. [30] M.S. Bapat, L. Shen, Y.J. Liu, Adaptive fast multipole boundary element method for three-dimensional half-space acoustic wave problems, Engrg. Anal. Bound. Elem. 33 (2009) 1113–1123. [31] H. Wu C, C.N. Wang, A study of fast multipole method on the analysis of 2D barrier, J. Mech. 25 (2009) 233–240. [32] M. Abramowita, I.A. Stegun, Handbook of Mathematical Function, Dover, New York, 1965. [33] S. Amini, On the choice of coupling parameter in boundary integral formulations of the acoustic problem, Appl. Anal. 35 (1990) 75–92. [34] T.W. Wu, Boundary Element Acoustics: Fundamentals and Computer Codes, WIT, Southampton, 2000. [35] Y.J. Liu, N. Nishimura, The fast multipole boundary element method for potential problems: a tutorial, Engrg. Anal. Bound. Elem. 30 (2006) 371–381. [36] K. Chen, P.J. Harris, Efficient preconditioners for iterative solution of the boundary element equations for the three-dimensional Helmholtz equation, Appl. Numer. Math. 36 (2001) 475–489. [37] S.A. Yang, A boundary integral equation method for two-dimensional acoustic scattering problems, J. Acoust. Soc. Am. 105 (1999) 93–105. [38] M.S. Bapat, Y.J. Liu, A new adaptive algorithm for the fast multipole boundary element method, CMES 58 (2010) 161–184.