PROCEEDINGS of the
23rd International Congress on Acoustics
9 to 13 September 2019 in Aachen, Germany

A new local boundary integral equation method for meshless
acoustic computations
Hakan DOGAN1 and Martin OCHMANN1;
1

Beuth University of Applied Sciences Berlin, Germany

ABSTRACTThe local boundary integral equation (LBIE) method has been introduced approximately three-decades ago,and has been successfully used in engineering problems such as elasticity, acoustics and electrodynamics.
The method is based on writing the local weak form of the governing equation. If an appropriate form of
the Green’s function can be found considering the governing equation in subdomains, the local integral
equation may have a convenient form, e.g. without domain integrals or the flux term on the subdomainboundaries. Considering the application of the methods to acoustic problems, four different test functions
have been employed so far in the literature. All of the resulting formulations contain domain integrals. In
this paper, we present a new form of Green’s function to be used in LBIE method for acoustic calculations,which yields a simple form of the local integral equation without domain integrals. This feature provides
significant improvement in terms of the computational time (CPU) requirements. The efficiency and theaccuracy of the new method are presented and compared with a Finite Element Method.
Keywords: Meshless methods, local boundary integral

1. INTRODUCTION
The local boundary integral equation (LBIE) method [1, 2, 3, 4, 5] has been introduced about three
decades ago as a robust numerical method for engineering problems; it is one of the version of themeshless local Petrov-Galerkin method [6]. The method is based on writing the local weak form in
spherical subdomains (which may be overlapping) and applying the divergence theorem. The resultant
local boundary and/or domain integrals are numerically discretized. The field variable needs to beapproximated at the numerical integration points using some kind of interpolation technique. The
Moving Least Squares (MLS) approximation has been employed for this purpose in majority of thestudies [1, 2, 3, 4, 5, 6, 7, 8], whereas the Radial Basis Function (RBF) interpolation technique has also
been used in some other works [9, 10, 11, 12, 13, 14].
The MLPG and LBIE methods have been applied to wave propagation problems in Refs. [2, 3, 4, 5, 7,
8, 9, 10, 11, 12], in order to solve for the Helmholtz equation. In Ref. [7] the Gaussian weight function wasused as the test function. Ref. [8] used the Heaviside step function as the test function. Refs. [3, 2, 4, 9, 10]
employed the modified fundamental solution of the Laplace equation as the test function to the Helmholtz
equation. Refs. [5, 11, 12] applied the modified fundamental solution of the Helmholtz equation as the test
function. The final form of the local integrals in all of these studies [2-5, 7-12] contain domain integrals.
In Ref. [15], we have reconsidered the derivation of an appropriate test function to be used in the weakform in meshless methods. Explicitly, we have derived the unique solution of the Dirichlet’s boundary
value problem in the local spherical subdomain with zero boundary condition value. The derived function isa special form of the Green’s function (tailored Green’s function) which satisfies the objective of
eliminating the domain integrals and the flux term from the local integral equation. As a result, a new formof the LBIE method has been established, which has improved features compared to the ones mentioned
above.

1hdogan@beuth-hochschule.de, ochmann@beuth-hochschule.de3229

2. FORMULATIONS
We consider the Helmholtz equation‫׏‬ଶ ‫݌‬ሺ࢞ሻ ൅ ݇ ଶ ‫݌‬ሺ࢞ሻ ൌ Ͳǡ     ࢞ ‫ܸ א‬ǡ

(1)where p is the acoustic pressure and ݇ is the wavenumber. The domain ܸ is enclosed by Ȟ ൌ Ȟ୳ ‫ ׫‬Ȟ୯with the boundary conditions‫݌‬ሺ࢞ሻ ൌ ‫݌‬଴ ǡ          Ȟ୳ ǡ(2)߲‫݌‬ሺ࢞ሻ‫ ݍ ؠ‬ൌ ‫ݍ‬଴ ǡ          Ȟ୯ ǡ   ߲݊where ‫݌‬଴ is the prescribed acoustic pressure on the essential boundary Ȟ୳ , ‫ݍ‬଴ is the prescribed potential
flux on the boundary Ȟ୯ , and ݊ is the outward normal vector on Ȟ.
In the traditional boundary element method (BEM), a weak formulation of the problem is formed in the
global solution domain ܸ using the free space Green’s function. In the LBIE, it is assumed that the
governing equation holds for any local subdomain and a set of ܰ spherical subdomains are distributed in
the domain ܸ and on the boundary Ȟ. The centers of the subdomains are called as source nodes where thevolumes in ࣬ ଷ for each node are denoted as π, such that ‫׊‬π ‫ܸ ؿ‬. The radius of each subdomain isdenoted as ܴπ . The surface (boundary) enclosing each interior subdomain is defined as πୗ ൌ ߲π. Whenthe source node lies on the global boundary Ȟ, the subdomain does not take the shape of a full sphereanymore. Instead, a half sphere or a slice of a sphere is formed depending on the shape of the global
boundary. For such cases, we set πୠ to denote the part of the subdomain boundary which intersects theglobal boundary Ȟ. For instance, for the subdomains which lie on Ȟ, we have ߲π ൌ πୠ ‫ ׫‬πୗ , where πୠ ‫ ؿ‬Ȟ and πୗ ‫( ܸ ؿ‬Fig. 1). The following weak formulation can be written for each local subdomain.
න ‫ ܩ‬ሺ࢞ǡ ࢟ሻ ሺ‫׏‬ଶ ‫ ݌‬ሺ࢞ሻ ൅ ݇ ଶ ‫݌‬ሺ࢞ሻሻ ݀π ൌ ͲǤ(3)π

An appropriate test function G should be used for Eq. (3). As mentioned in the Introduction, the previous
meshless methods have employed different types of test functions, such as the Heaviside step function, theGaussian weight function, the modified fundamental solution of the Laplace equation or the modified
fundamental solution of the Helmholtz equation. In Ref. [15], a new test function has been derived which is
the unique solution of the Dirichlet’s problem for a spherical subdomain with the boundary value zero. The
derived test function is given explicitly as:‫ܩ‬ሺ࢞ǡ ࢟ሻ ൌͳ ݇ሺܴπ െ ȁ࢞ െ ࢟ȁሻǤͶߨȁ࢞ െ ࢟ȁ ܴ݇π(4)

Using the tailored Green’s function (4) in Eq. (3), and applying integration by parts, we obtain the integral݇(5)  ‫݌‬ሺ࢟ሻ െන ‫݌‬ሺ࢞ሻ ߲݀πሺ࢞ሻ ൌ ͲͶߨܴπ  ܴ݇π డπfor local subdomains which lie fully inside the domain.
When the source point y lies on the global boundary, the subdomains are taken as a part of the sphere(e.g. half sphere). In this case, πୱ denotes the boundary which lies within the global volume ܸ, πୠ
indicates the boundary which intersects with Ȟ. The source point (center of the local subdomain) lies onπୠ , i.e. ࢟ ‫ א‬πୠ ‫ ؿ‬Ȟ. The local integral equation for the boundary nodes takes the form:߲݇‫݌‬ሺ࢞ሻߣ ‫݌‬ሺ࢟ሻ െන ‫݌‬ሺ࢞ሻ ߲݀πୗ ൌ න ‫ܩ‬ሺ࢞ǡ ࢟ሻ߲݀πୠ ǡ
(6)Ͷߨܴπ  ܴ݇π 
߲݊
π౏

πౘ

where ߣ is a parameter analogous to the geometric coefficient in the BEM, which is given here asߣሺ࢟ሻ ൌ െ

ͳ
߲‫ݎ‬
ͳ
߲
න
 ൬ ൰ ߲݀πୗ ǤͶߨ ߲݊ሺ࢞ǡ ࢟ሻ ߲‫ݎ ݎ‬πೞ(7)

Eqns. (5) and (6) need to be numerically calculated. This requires, first, the discretization of the boundary
integrals in these equations and then the interpolation of the unknown p at the discretization points. Unlike
the boundary element method, there is no need to use elements to discretize the surface integral in the LBIE

3230

method. The details of the quadrature and interpolation methods in the current meshless method can befound in detail in Ref. [15].

Figure 1 – The sketch of the spherical local subdomains π, and the enclosing surfaces πୗ . For the
subdomains on the global boundary, πୠ defines the surface which intersects with the global boundary Ȟ(the filled blue area).

3. NUMERICAL EXAMPLESIn this section, the numerical examples will be shown in order to demonstrate the efficiency and theaccuracy of the new method. The problem domain is chosen as a rectangular cuboid filled with air. The
density of air is taken as ߩ ൌ ͳǤ͵ kg/m3 and the speed of sound in air is taken as ܿ ൌ ͵ͶͲ m/s.
In the simulations, the nodes are distributed uniformly in the domain with spacing ݄ between two adjacentnodes. The boundary integrals were calculated numerically using ܰୖ୆୊ ൌ ͷͲ neighbor nodes in the radialbasis function interpolation and ܰீ ൌ ͳͲ quadrature points (in each dimension). The computations were
performed with MATLAB software on a desktop workstation computer with 3.20 GHz processor speed and
128 GB of RAM. The resultant system matrices were solved iteratively using the gmres function. The
accuracy of the numerical results is assessed by evaluating the L1 norm error݁௅భ ሺ‫ݔ‬ሻ ൌȁ‫݌‬ሺ࢞ሻ െ ‫݌‬ୣ୶ୟୡ୲ ሺ࢞ሻȁǡȁ‫݌‬ୣ୶ୟୡ୲ ሺ࢞ሻȁ(9)and the L2 norm error݁௅మ ൌ ඩሺ௜ሻሺ௜ሻσே௜ୀଵൣ‫݌‬୬୳୫ୣ୰୧ୡୟ୪ െ ‫݌‬ୣ୶ୟୡ୲ ൧
ଶ

ሺ௜ሻ
σே
௜ୀଵൣ‫݌‬ୣ୶ୟୡ୲ ൧ 

ଶ

Ǥ

(10)

3.1 Standing waves in a duct
We consider a domain with length L=0.4m, width w=0.15m, and height H=0.06m. The inlet and theoutlet of the rectangular duct is located on y-z plane, at x=0 and at x=L, respectively. The remaining
boundaries, which are on the x-y and x-z planes, are rigid walls. In particular, a constant pressure isprescribed at the inlet, i.e.
‫݌‬୧ ȁ௫ୀ଴ ൌ ͳͲͲͲ ܲܽǡ(11)and, a full-reflection condition is given at the outlet:

3231

(12)

߲‫݌‬
ൌ ͲǤ
ฬ
߲݊ ௫ୀ௅

For the conditions given above, the problem reduces to 1D wave propagation in x-direction. One can showthat the exact solution to the problem defined in Eqns. (1), (11) and (12) is given by‫݌‬ୣ୶ୟୡ୲ ሺ‫ݔ‬ሻ ൌ ‫݌‬୧(13)

  ݇ሺ‫ ݔ‬െ ‫ܮ‬ሻǤ ݇‫ܮ‬The computational time and the accuracy of the method are presented in Table 1. The frequency is
set as f=5000 Hz. The second column in Table 1 shows the total number of nodes N for thecorresponding value of h, hence the degrees of freedom (DOF) being solved for. The 3 rd column
indicates the CPU time required for the preparation of the matrix coefficients for the nodes on theboundary and in the domain. The 4 th column presents the total CPU time including the preparation of
the system matrix and the iterative solver. The results are compared to the formulation in Ref. [12].
Note that the method in [12] in the Mach number limit M՜0 becomes the version of the LBIE whichcontains domain integral. It can be seen that the accuracy of the present method and Ref. [12] arecomparable and show only minor difference, which may be because of the different integration
procedures. However, the computational time requirement is significantly improved with the new
method. For a Gaussian quadrature integration with ܰீ weights, a domain integration needs ̱ܰீଷoperations, whereas a surface integration is performed with ̱ܰீଶ operations. Therefore, a time
improvement with factor ͳȀܰீ is achieved. This can be easily noticed by comparing the time required
for the preparation of the system matrix, i.e. the CPU time for the assembly. The CPU time for solving
the matrix is also significantly improved as can be seen in Table 1.
Table 1 – CPU time and L1-norm error for the case f =5000 Hz
Current Method

Ref. [12]

h

Nodes

CPU Time

CPU Time

L1 Error

CPU Time

CPU Time

L1 Error

(mm)

(N)

Assembly (sec)

Total (sec)ࢋࡸ૚ (%)

Assembly (sec)

Total (sec)ࢋࡸ૚ (%)

10

9594

5.4

20.4

1.36 e-349809.89 e-4

8

17952

6.6

29

1.06 e-3591049.44 e-4

5

70875

11

68

5.53 e-4102

2265.36 e-4

2.554434125

2172

8.71 e-5

205

4732

7.49 e-5

3.2 Comparison with COMSOL FEM
The accuracy of the method is further presented in Table 2, and compared with Ref. [12] and with
the COMSOL FEM software. A uniform mesh with h=5 is created for the meshless methods. Similarly,
quadratic cubic elements with 5 mm side length are used in FEM computations, which result in 70875
domain elements (see Fig. 2a). Although equal number of nodes are used both in meshless methods andin FEM, the DOF being solved for is much larger (~5.5 x 10 5 ) in the latter. An important dimensionlessparameter to assess the resolution of the mesh for wave propagation problems is ߣȀ݄, which isreported in the second column. The L 2 -norm errors decrease consistently with increasing DOF for all
formulations, indicating that the methods are convergent. However, the meshless methods aresubstantially more accurate. Although not presented in Table 2, the CPU time requirements of themeshless method is comparable to the FEM solution. The total CPU time spent in the present method
is 68 seconds, whereas the FEM simulation takes 42 seconds. Though, one should note that the meshdiscretization, the mesh connectivity and the system matrix assembly are generally done as
pre-processing operations in the FEM software. As an example, the wave propagation result at 5000
Hz using the FEM is presented in Fig. 2b.

3232

Table 2 – L2-norm error for dx =5 mm.
Current Method

Ref. [12]COMSOL

Freq. (Hz)

λ/h

ࢋ ࡸ૛

ࢋ ࡸ૛

ࢋ ࡸ૛

15000

4.57

0.00460.00470.01694

10000

6.86

0.001

0.00110.1128

5000

13.72

9.67 e-49.50 e-4

0.0438

2000

34.29

5.66 e-4

1.83 e-40.0160

1000

68.60

1.20 e-4

1.09 e-5

0.0073Figure 2 – a) Uniform FEM mesh with quadratic elements with side length 5 mm. b) FEM solution for thereal part of acoustic pressure at f=5000 Hz

4. CONCLUSIONSIn this paper, a new LBIE method for the solution of the Helmholtz equation is presented. Unlike theprevious formulations in this category, the present method does not contain any domain integral or the flux
term. In the final form of the local integral, the acoustic pressure is the only unknown which needs to beapproximated by any kind of mathematical interpolation. The computational time requirements of themethod have shown to be ̱ͳȀܰீ times of that encountered in the predecessors of the method. Moreover,the method yields a real-valued, band-limited system matrix. Provided that more efficient linear solvers
appropriate for these features of the system matrix are used, the method can be used as a better alternative
for the existing numerical FEM and BEM techniques. The comparison of the method with an FEM software(COMSOL v4.3 with quadratic elements) for identical discretization of the problem domain reveals at least
one order of magnitude more accurate results for the meshless method. The proposed new LBIE method
can be applied to many other types of wave propagation problems, provided that circular local subdomains
are used in the implementation.

3233

ACKNOWLEDGEMENTS
This work was funded by the German Ministry of Education and Research (Bundesministerium für
Bildung und Forschung) within the research project High Efficiency Low Noise Heatpump Dryer(HELNOISE).

REFERENCES[1] T. Zhu, J. -D. Zhang, S. N. Atluri. A local boundary integral equation method in computational
mechanics, and a meshless discretization approach. Computational Mechanics, Vol. 21, pp. 223-235,
1998.
[2] T. Zhu, J. Zhang, S. N. Atluri. A meshless numerical method based on the local boundary integral
equation (LBIE) to solve linear and non-linear boundary value problems. Engineering Analysis with
Boundary Elements, Vol. 23, pp. 375-389, 1999.
[3] D. J. Fu, H. B. Chen, P. Q. Zhang. Improved non-singular local boundary integral equation method.
Applied Mathematics and Mechanics (English Edition), Vol. 28 (8), pp. 1093-1099, 2007.
[4] H. B. Chen, D. J. Fu, P. Q. Zhang. An inverstigation of wave propagation with high wavenumbers via
the regularized LBIEM. Computer Modeling in Engineering and Sciences, Vol. 20 (2), pp. 85-98,
2007.
[5] J. Sladek, V. Sladek, S. N. Atluri. Application of the local boundary integral equation method toboundary-value problems. International Applied Mechanics, Vol. 38 (9), pp. 1025-1047, 2002.
[6] S. N. Atluri, T. Zhu. A new meshless local Petrov-Galerkin (MLPG) approach in computational
mechanics. Computational Mechanics, Vol. 22, pp. 117-127, 1998.
[7] D. J. Soares. An iterative time-domain algorithm for acoustic-elastodynamic coupled analysis
considering meshless local Petrov-Galerkin formulations. CMES: Computer Modelling in Engineering& Sciences, Vol. 54 (2), pp. 201-221, 2009.
[8] A. Tadeu, P. Stanak, J. Sladek, V. Sladek. Coupled BEM-MLPG acoustic analysis for
non-homogeneous media. Eng. Anal. Boun. Elem., Vol. 44, pp. 161-169, 2014.
[9] H. Dogan, V. Popov. Numerical simulation of the nonlinear ultrasonic pressure wave propagation in acavitating bubbly liquid inside a sonochemical reactor. Ultrasonics Sonochemistry, Vol. 30, pp. 87-97,
2016.
[10] E. H. Ooi, V. Popov. A simplified approach for imposing the boundary conditions in the local
boundary integral equation method. Computational Mechanics, Vol. 51, pp. 717-729, 2013.
[11] H. Dogan, V. Popov, E. H. Ooi. Dispersion analysis of the meshless local boundary integral equation
and radial basis integral equation methods for the Helmholtz equation. Engineering Analysis with
Boundary Elements, Vol. 50, pp. 360-371, 2015.
[12] H. Dogan, C. Eisenmenger, M. Ochmann, S. Frank. A LBIE-RBF solution to the convected wave
equation for flow acoustics. Engineering Analysis with Boundary Elements, Vol. 92, pp. 196-206,
2018.
[13] E. J. Sellountos, D. Polyzos, S. N. Atluri. A new and simple meshless LBIE-RBF numerical scheme inlinear elasticity. CMES: Computer Modeling in Engineering and Science, Vol. 89 (6), pp. 513-551,
2012.
[14] E. J. Selluntos, A. Sequeira. An advanced meshless LBIE/RBF method for solving two-dimensional
incompressible fluid flows. Computational Mechanics, Vol. 41, pp. 617-631, 2008.
[15] H. Dogan, M. Ochmann. A new test function for acoustic computations with meshless methods.
Journal of Theoretical and Computational Acoustics, Vol. 27 (1), pp. 1940001 (1-18), 2019.

3234