BEM with the Real-Valued Fundamental

Solutions for the Helmholtz Equation

© Mykola Yas'ko

Dnipropetrovsk, 49050, Ukraine

e-mail: nyasko@hotmail.com

Abstract

The numerical experiment for the 2D/3D problems shows that both the complexity, the memory storage and the required CPU time of the method are approximately same as for the corresponding Laplace equation.

INTRODUCTION

The Helmholtz equation arises in many physical applications, in particular in acoustics. The boundary element method (BEM) is the powerful tool for the solution of many linear partial differential equations and, in particular, for the Helmholtz equation.

The traditional BEM is derived from a boundary integral equation (BIE) by dividing the domain boundary into boundary elements and applying collocation method to obtain the solution. The modern variants of BIE for the Helmholtz equation are based on the complex-valued fundamental solutions (Green's functions) in the two-dimensional (2D) and three-dimensional (3D) cases. This approach made the computational problem for the 2D/3D Helmholtz equation more complicated than, for example, for the corresponding Laplace equation. The main difficulties are connected with the computations for the high wavenumbers, when we have to use a large boundary element mesh. The memory storage is very large and CPU time is very long for solving of the large linear system with complex-valued coefficients. The modern state-of-the-art for the BEM-based solutions is shown in many reviews, in particular Brebbia [6], Banerjee and Butterfild [3], Amini and Kirkup [2].

In this work a new formulation of the boundary element method is suggested for the fast computing of the Helmholtz equation. It is a variant of the direct integral equation formulation based on the real-valued fundamental solutions and the application of the conjugate gradient methods for a resulting linear system.

The mathematical formulation of the new real-valued fundamental solutions is shown in section 2. Basics of the direct boundary element method are presented in section 3. Benefits of conjugate gradient methods and count of iterations are discussed in section 4, where the obtained results are presented for interior and exterior 2D/3D boundary-valued problems. Solution times are compared for two variants of the conjugate gradient method and Gaussian elimination method. The precision of obtained numerical results for the some 2D/3D problems is discussed in section 5, where comparison are given for the exact solutions.

MATHEMATICAL FORMULATION

The Helmholtz equations for the complex-valued potential u may be written in the 2D/3D space as

Lkn [u] = Dn u + k2u = 0

(1)

where k is a real-valued wavenumber; n is space dimension (2 or 3). The equation is valid in the domain D with boundary S. We consider the linear boundary-valued problems with the linear boundary conditions of the irst, second or third kind. The boundary condition on the boundary S takes the following form:

au + bu,n = f,

(2)

where a, b are real-valued coefficients; n is the unit outward normal at the point P on the boundary.

Today for solving of the Helmholtz equation by the boundary element methods the fundamental solutions (Green's functions for the infinite medium) in the complex variables [1,2,3] are used as a rule. In two-dimensional case it is

FkMP = i 4 H(1)0 (krMP),

(3)

where H⁽¹⁾₀ is the Hankel function of the first kind of order zero; r_MP is a distance between the points M and P; k is a wavenumber.

In three-dimensional case the fundamental solution is

FkMP = 1 4p e-ikrMP rMP

(4)

This approach takes a lot of memory storage and CPU time in comparison with analogous Laplace equation. Utilization of the real-valued fundamental solutions will give a reduction of memory storage in half and simplification and speeding-up program codes for computation of solution of the boundary-valued problems for the Helmholtz equation.

In the present work for two-dimensional case a new real-valued fundamental solution was used

FkMP = 1 4 Y0 (krMP),

(5)

where Y₀ is the Bessel function of the second kind of order zero (called as Neumann function of order zero also) [1]. This function is a real part of the complex-valued Hankel function of a real argument. The function has a logarithmic singularity for r® 0 and the suitable conditions at ¥.

In a tree-dimensional case the fundamental solution was

FkMP = 1 4p coskrMP rMP

(6)

i.e., it was a real part of the corresponding complex-valued fundamental solution (4) also.

Let us j is a real or imaginary part of complex-valued potential u. Then an integral representation for j in domain D may be written as

jM = ó õ S  (jP,nFMP - jPFMP,n) dSP,

(7)

where point P Î S and M Î D. Below we'll name point M as viewpoint.

BOUNDARY ELEMENT METHOD

The boundary element method is derived from the integral representation (7) by discretizing the boundary and through the application of collocation method. This approach is well known and described in details in references [1,2,5]. The discretization of (7), which leads to the classic 'boundary element method' technique is described shortly below. In the constant boundary element method, the above integral equation is solved numerically by dividing the boundary S into N boundary elements, in each of which j and j_,n are approximated by constants.

In a 2D case the boundary was approximated by straight line elements. In order that normal to the boundary outward the two nodes that define each element must be numbered in the anti-clockwise direction when it viewed from just inside the domain.

In a 3D case the boundary was represented as a set of triangular panels. In order that normal to the boundary outward the tree nodes that define each element must be numbered in the anti-clockwise direction when it viewed from just outside the domain.

We denote these values by j_i and j_i,n, i = 1,¼N; and apply equation (7) at one nodal point M_i in the center of each boundary element to obtain

1 2 ji = N å j = 1  (jj,n ó õ Sj  F iPdS P-jj ó õ Sj  F iP,ndS P ) = N å j = 1  (Aijjj,n+Bijjj),

(8)

where S_j denotes integration over the jth boundary element.

Coefficients A_ij and B_ij of the linear system of equations (8) are integrated numerically over boundary elements using the following approach.

When the viewpoint M lies far from the boundary element, the integrand is continuously differentiable and hence standard quadrature rules are satisfactory. The quadrature order decreased with the increasing distance from viewpoint to the center of the boundary element.

However, when M lies near the element or on the element, the integrands are became near-singular ones and standard quadrature rules are not applicable. If the viewpoint is situated near from boundary element then the special techniques applied involve `subtracting out' the singularity and evaluating the singular part analytically and remaining regular part numerically by the same quadrature rules.

For example, for the 3D case it was

Aij = ó õ Sj  FkiP dSP = ó õ Sj  (FkiP- 1 4priP ) dSP + 1 4p ó õ Sj  1 riP dSP

(9)

Bij = ó õ Sj  FkiP,n dSP = ó õ Sj  [FkiP,n- 1 4p ( k2 2riP + riP nP r3iP )] dSP + 1 4p ó õ Sj  ( k2 2riP + riP nP r3iP ] dSP

(10)

where, in (9) and (10) the first integral is non-singular. Evaluation in this way requires the computation of the regular integral (amenable to standard quadrature) and the determination of the subtracted out part.

The regular integrals in (9, 10) that arise are approximated by a quadrature rule defined on a triangle. Paper of Laursen and Gellert [7] contains a selection of Gauss-Legendre quadrature rules for the triangle. The remaining parts in (9) and (10) are integrated analytically over boundary elements.

The described procedure was obtained very precision and very fast method for coefficient computation.

Eliminating the j_i or j_i,n from each element on the boundary by applying the boundary-valued condition in each nodal point, we thus obtain from (8) a system of the N linear algebraic equations with N unknowns. The system of N linear algebraic equations with the full matrix was solved by the different methods that described in the next section.

SOLUTION OF A LINEAR SYSTEM

Traditionally, the direct methods have been employed to solve the resulting linear system such as the Gaussian elimination method. However, in recent years, interest has grown in iterative solvers, and in particular, the use of the conjugate gradient methods has been investigated by a number of researches ([5], [8], [10]). The field of iterative methods for solving systems of linear equations is in constant flux, with new methods and approaches continually being created, modified, tuned, and some eventually discarded. At the present work the three methods were used for the numerical solving of the linear equation system.

• Conjugate gradient method (CGM) modified for the unsymmetrical systems. The description of this method was took from Peyret and Taylor [9]. The convergence of this method may be slow, since the spectrum of the modified matrix will be less favourable than the spectrum of the initial matrix.
• Biconjugate gradient method (BiCGJ) with Jacobi preconditioner. Algorithm of this method is described in book [4].
• Gaussian elimination method (GaussElim). This direct method based on Gaussian elimination with partial pivoting.

BiCG method, like CG, uses limited storage. It is suitable when the matrix is asymmetric and nonsingular; however, convergence may be irregular [4], and there is a possibility that the method will break down [10]. CGM requires the four matrix-vector multiplication at each iteration. BiCG requires the three matrix-vector multiplication at each iteration. The iterative solvers were used with Jakobi preconditioning; that is, in general, the simplest form of preconditioner.

The comparison of the different methods efficiency for the linear system was made for the 3D problems. Domain boundary was a sphere with radius r = 1. Boundary-valued conditions of the first and the second kinds were

u = sinkx,                                                             u,n = k nx coskx,

i.e., the real part just presented in a solution.

In general for all computations the numerical solutions were obtained by all methods (we did not compute the cases when k coincided with eigenfrequencies). Comparison of CPU time and number of iterations are shown in table 1 for the boundary-valued problems of the first and second kind for the 3D Helmholtz equation. For the iterative methods the iteration is stopped when difference between two approximate solutions was less than 10^-9. All computations were performed on the IBM PC with Pentium-166 processor. Numerical examples are presented for some interior and exterior boundary value problems of the first and the second kind for some wavenumbers. There was no goal to investigate the influence of wavenumbers on the iterative methods. The results are shown in table 1 demonstrate the considerable benefits of conjugate gradient methods in comparison with Gaussian elimination. This benefits increase with the increasing the number of equations.

The following trends are evident in this table:

• Iteration count are relatively insensitive to mesh size. Indeed, the counts often decreased when number of the boundary elements increased.
• Iteration count are especially small for the Neumann boundary problem.
• The results are clear demonstrated the benefits from conjugate gradient methods in comparison with Gaussian elimination method.

Many of these trends are the same as for finite elements (see [8]), but the sensitive to wavenumber is considerably less pronounced.

NUMERICAL EXAMPLES

In this section, we present the results of numerical experiments to solve the 2D/3D Helmholtz equation. The objective for all computations was the test of the wavenumber influence for accuracy and reliability of the suggested method. The influence of mesh size on performance of the computed solution was not investigate.

For the test of the computation coefficients (9) and (10) the two boundary problems of the first and second kinds were computed in the circular domain with analytical solution

u = eikx.

(11)

The mesh with 500 boundary elements was used in both cases. The maximum of the absolute error was approximately 10^-3 for the wavenumbers in range [1¸10].

For the exterior problem a test problem for the scattering wave was computed for the circular cylinder. There is a case of scattering problem, when incoming harmonic wave interacted with the obstacle (circle cylinder). The incoming harmonic wave u_i defined by (11) with the Neumann boundary condition for the scattering wave

¶us ¶n = - ¶ui ¶n

The problem was computed for the scattering wave. Cylinder was divided in 1000 boundary elements also. The real part of the computed solution for the exterior problem is shown in figure 1. The figure are drown as greyscale picture, where the white color are corresponding the maximum of function and black color is corresponding the minimum.

Figure 1: The real part of the complex-valued acoustic potential u for the case of the scattering by cylinder (k = 12).

For the 3D case the test problem was for the interior of sphere with the boundary represented by 512 triangular elements. The analytical potential defined by (11) also. The absolute error maximums

Dumax = max |unum-uanal|

are presented in table 2 for some wavenumbers. We can see for the long waves (k < 2) the method gives a good precision. But with the increasing wavenumber to 4 the maximum error increases dramatically. It is connected with the properties of the solutions of the Helmholtz equation and it signify that size mesh is not enough.

wavenumber	maximum error
k	1st kind problem	2nd kind problem
0.5	0.008	0.003
1	0.018	0.005
2	0.038	0.033
3	0.088	0.060
4	4.301	1.050

This problem exists for all boundary elements methods but the described approach will be more economical than other methods.

SUMMARY

Boundary element methods are becoming increasingly popular as methods for the numerical solution of linear elliptic partial differential equations such as the Helmholtz equation. The advantages of the proposed variant of the BEM for the 2D/3D Helmholtz equation arises from the fact that we can separately solve boundary-valued problems for the real and imaginary parts of the complex-valued potential. In this case the BEM technique for the Helmholtz equation is not more complex than for the corresponding Laplace equation. The computational time and memory storage for this approach are comparative as for the Laplace equation too.

Performance of conjugate gradient method is relatively insensitive to number of boundary elements and wavenumber, and the corresponding algorithms are highly parallelizable. Comparisons with analytical data suggest that the solutions obtained by the present numerical method are quite accurate. The proposed method will be especially effective one for the high wavenumbers when size of the boundary element mesh has to be very large.

REFERENCES