A Pollution Eﬀect in the Boundary Element Method for Acoustic Problems

The pollution eﬀect is a well-known and well-investigated phenomenon of the ﬁnite element method for wave problems in general and for acoustic problems in particular. It is understood as the problem that a local mesh reﬁnement cannot compensate the numerical error which is generated and accumulated in other regions of the model. This is the case for the phase error of the ﬁnite element method which leads to dispersion resulting in very large numerical errors for domains with many waves in them and is of particular importance for low order elements. Former investigations have shown that a pollution eﬀect resulting from dispersion is unlikely for the boundary element method. However, numerical damping in the boundary element method can account for a pollution eﬀect. A further investigation of numerical damping reveals that it has similar consequences as the phase error of the ﬁnite element method. One of these consequences is that the number of waves within the domain may be controlling the discretization error in addition to the size and the order of the boundary elements. This will be demonstrated in computational examples discussing traveling waves in rectangular ducts. Diﬀerent lengths, element types and mesh sizes are tested for the boundary element collocation method. In addition to the amplitude error which is due to numerical damping, a rather small phase error is observed. This may indicate numerical dispersion.


Introduction
It is common practice to use a certain fixed number of boundary elements per wavelength in engineering simulations. The author himself has run many test cases on this and published a number of papers 19,22,28 discussing this topic. In one of them, 28 there has been the statement This is an Open Access article published by World Scientific Publishing Company. It is distributed under the terms of the Creative Commons Attribution 4.0 (CC-BY) License. Further distribution of this work is permitted, provided the original work is properly cited. that in a duct example, the authors "could not find any indication of a pollution effect in the frequency range of 0 < kl ≤ 160π" which has been equivalent to 80 waves in the duct. Actually, there was no indication that convergence had been influenced by frequency. Large wavenumbers kl did neither influence convergence rates nor the linearity of the (logarithmic) error as a function of the (logarithmic) frequency. A similar statement is found in the recent mathematical literature. 8 However, in these papers investigating the influence of the normalized wavenumber kl, the emphasis was usually put on the dependence on k. In particular, the author's papers have not investigated a dependence on the duct length l while keeping k (and the element size) constant.
The dependence of the numerical error on the length of the computational domain was first reported in a paper by Bayliss et al. 4 Investigation of this problem became very popular during the decade of the 1990ies, when many papers on the pollution effect as an accumulated phase error were published, see for example. 2,3,5,7,15,16,35 Reviews of this, further findings and conclusions were discussed thereafter. 1,9,10,14 However, the most comprehensive and detailed review of this topic was published in the very interesting book by Ihlenburg. 13 Many methods have been proposed to overcome the unwanted dispersion error. In recent years, the smoothed finite element method 11,37 has received much attention in literature. In a recent study on the accuracy of eigenfrequencies for structural components using commercial software (Abaqus), this effect could not be found. 17 The phase error as a result of a dispersive solution is also an important problem in the finite difference method. A very recent publication on this topic is by Wang and Wong 36 who directly correspond to the finite element community by choosing the title of their work very similar to the paper by Babuška and Sauter. 3 As mentioned above, there has not been any indication for a significantly dispersive solution in the boundary element method (BEM). Apparently, this is the reason why the author was unable to find any indication of a pollution effect in the boundary element literature. Rather recently, the author published a brief paper on the problem of numerical damping in the conventional boundary element method. 24 This paper describes numerical damping observations of both, findings in literature 39 and the author's own simulations. 18 Furthermore, it concludes that a real fundamental solution (Green's function) may solve the problem of numerical damping. Unfortunately, the real fundamental solution method is subject to numerical instabilities, 30 see also the discussion in. 24 Independently of the author's activities, Fahnline published a paper on the effect of numerical damping in the boundary element method some years before. 6 While the effect of numerical damping has been described several times before, the conclusion that it may be understood as a pollution effect has not been drawn so far. The computational test cases in this paper, i.e. problems of traveling waves in a duct, are well suited to demonstrate the effect and have not been presented and discussed in a similar way in the references mentioned above.
In the undamped case of a duct with excitation at one end, 24 infinite resonance peaks are expected at the eigenfrequencies. In reality, these peaks are finite and loose height with frequency, i.e. numerical damping is increasing with frequency. It was shown in a recent paper by Zheng et al. 39 that the eigenvalues of the boundary element formulation are not purely real but exhibit a small imaginary part. This small imaginary part is not physical but results in a numerical damping causing the amplitude peaks to be limited in height. Therefore, it is obvious that the numerical error very close to the resonance peaks will always be quite high and thus, controlling this error by a fixed number of elements will be rather difficult. However, it was further shown in the author's paper 24 that even a traveling wave in a duct is damped. This effect is the starting point for the current paper which approaches the numerical damping and, with it, the pollution effect from the numerical example of traveling waves in a long duct. The author has extensively studied this problem before 18,19,22,[24][25][26]28,33 and it was even proposed as a suitable benchmark problem 12 for the benchmark platform of the European Acoustics Association. The main motivation of studying this problem in the context of a pollution effect results from the thought that, according to the damping, the numerical error in the duct should become greater the longer the duct is. Consequently, the numerical error would depend not alone on the element size but also on the length of the domain or the number of waves within the domain. This is exactly what happens in the finite element method and what was described in many papers, see for example 9,10,14 and references therein.
This paper is organized as follows: It starts with a formulation of the boundary value problem and the boundary element formulation with discretization by collocation and a selection of boundary elements. The main part reviews and investigates the problem of traveling waves in a long duct and behaviour of different types of elements again. Finally, it roughly estimates the phase error compared to the finite element method. The paper will be completed by a brief conclusion.

Boundary Element Formulation
The derivation of the boundary element formulation will be presented in an abridged way and is reduced to a few points which are relevant for understanding of this paper. A more detailed presentation of the author's approach to BEM for acoustic problems is found in the concept chapter 27 and in the lecture notes. 25

Helmholtz equation, boundary conditions and integral equation
The time-harmonic wave equation, i.e. the Helmholtz equation, for the sound pressure p in a domain Ω with wavenumber k = ω/c is considered. The angular frequency and the speed of sound are denoted as ω and c, respectively. A harmonic time dependence of e −iωt is assumed. Admittance boundary conditions are applied on the boundary Γ. They are equivalent to Robin boundary conditions which may degenerate to Neumann boundary conditions if the admittance is zero. This condition is written as Y represents the boundary admittance which relates the sound pressure to the difference between the normal components of the fluid particle velocity v f and the underlying structural particle velocity v s . The normal component of the fluid particle velocity is related to the normal derivative of the sound pressure p by means of the linearized Euler equation in frequency domain as where 0 represents the ambient density of the fluid. A residual or weak formulation, partial integration and substituting the Green's function for the test function leads to the Kirchhoff-Helmholtz integral equation with admittance boundary conditions as The Green's function (or fundamental solution) is given as with r being the Euclidean distance between field point x and source point y as r( x, y) = | x − y|. The coefficient c in Eq. (4) is c = 1 for y ∈ Ω and 0 < c < 1 for y ∈ Γ. It is c = 0.5 for a point y on a smooth part of the surface Γ. Since G is a particular solution of the inhomogeneous Helmholtz equation, it is not unique. For example, the real part of G can be used as a stand-alone Green's function too, see the examples and the discussion in. 24,30

Discretization by collocation and choice of elements
For discretization, Eq. (4) is tested with delta functions. Furthermore, all physical quantities are approximated by using Lagrangian interpolation polynomials. As a result, the system of linear algebraic equations is yielded as 25,27 (H − GY )p = Gv s .
The system matrices G and H are complex, fully populated and non-hermitian. Y , p and v s are the (close to) diagonal matrix of the nodal admittance values and the column matrices of the nodal values of sound pressure and structural particle velocity, respectively. For this study, continuous and discontinuous Lagrangian boundary elements are used. Among the continuous elements, linear and quadratic elements are denoted as P 1c and P 2c , respectively. Furthermore, discontinuous elements will be applied as constant elements, i.e. P 0 , and linear elements with nodes at the zeros of the Legendre polynomials, i.e. P 1d . These elements have been described and discussed extensively in former work of the author. 22,25,28 The system of equations in (6) is solved iteratively using a GMRes algorithm without preconditioning, cf. Saad and Schultz 34 and also the author's contribution. 29 A very low residual of R n ≤ 10 −10 , see Eq. (4) of Marburg and Schneider, 29 is demanded to avoid any influence of the iterative solver on the accuracy of the numerical solution.

Description of boundary element models
Air-filled ducts of lengths l = n · l Ref with l Ref = 3.4 m and a square cross section of w 2 = 0.2 × 0.2 m 2 account for the examples in this section. They are well suited test cases since the three-dimensional numerical solution can be compared to the one-dimensional analytical solution, at least up to a frequency where no modes perpendicular to the length of the duct occur. Assuming a speed of sound of 340 m/s, an ambient density of 1.3 kg/m 3 and the width of w = 0.2 m, the lowest perpendicular modes are observed at a frequency of 1700 Hz. a A solution with traveling plane waves is found if a particle velocity is applied to one end, e.g. at x = 0 it is v s (x = 0) = v 0 , and a fully absorbing boundary condition is applied to the other end, e.g. at x = n · l Ref . Full absorption is achieved by an admittance of 0 cY (x = n · l Ref ) = 1. All other walls are considered acoustically rigid and at rest, i.e. Y = 0 and v s = 0. A vivid description of this example is given in Fig. 1. This configuration leads to a solution of constant sound pressure magnitude everywhere in the duct and at all frequencies. Only the phase angle varies. The solution of the one-dimensional problemp is given asp Independence of the sound pressure magnitude from position and frequency makes this example an ideally suited one to compare different configurations. This has been the motivation to propose this example as a benchmark case for computational acoustics. 12 Setting the particle velocity v 0 = 1 m/s over the entire frequency range is not a realistic assumption but just, for convenience, used for the computational test case of this section.
Travelling wave (plane wave) in duct a Theoretically, a half wave perpendicular mode occurs at 850 Hz. However, in all these and former studies of this duct problem, this mode has hardly ever been observed. The author assumes that the invisibility of this mode is due to symmetry conditions in this example because the entire configuration shows a double symmetry over the cross section and the half wave is asymmetric.

S. Marburg
The boundary element mesh consists of square elements only. For a specific model, all of them have the same size which is controlled by the parameter m s . This parameter defines the number of elements over the width of the duct. This results in element edge lengths of 0.2 m, 0.1 m and 0.0667 m for m s = 1, m s = 2 and m s = 3, respectively. Because of the square shape of the elements, it is easy to evaluate the number of elements over the length. It is 17 for m s = 1 and n = 1 and goes up to 816 for m s = 3 and n = 16, see Table 1. Figure 2 shows a selection of different meshes for different configurations.
As mentioned above, the ducts may be of different length while the cross section is always the same. The basic configuration is the same duct which has been investigated in many papers of the author 18,19,22,[24][25][26]28,33 with n = 1, i.e. with a length of l Ref = 3.4 m. In addition, the cases of n = 2, 4, 8 and 16 are considered. For preparation of this manuscript, the author investigated other (integer) values of n < 16 without qualitatively new findings. The example has been created such that at 100 Hz and n = 1, one full wave is found along the duct. This is equivalent to the normalized wave number of kl = 2π. Increasing the frequency increases the number of waves with the same rate as by increasing the length. Consequently, in both cases, for 200 Hz with n = 1 and for 100 Hz with n = 2, two waves are found along the length, which is equivalent to kl = 4π. One case will be presented for 1500 Hz and n = 16. This results in 15 · 16 = 240 waves over the length of the duct, meaning that kl = 480π.
A selection of numerical solutions is shown in Fig. 3. For n = 1 (see top two subfigures), the duct is filled with 1.5 and with 15 waves at 150 Hz and 1500 Hz, respectively. 60 waves   are found for the four times longer duct at 1500 Hz, while the longest duct, i.e. n = 16, shows 24 waves at the low frequency of 150 Hz already. The results stem from computations with a mesh of P 2d elements and m s = 3. The color bars indicate minimum and maximum values of the real part of the sound pressure. Analytically, these values would be given as ±442 Pa.

Local numerical error over the duct length
In this subsection, the sound pressure and its numerical error at field points inside the duct are investigated. Recently published results on numerical damping in BEM indicated that numerical damping is not just visible in the undamped case but is also observed as a decaying sound pressure of a traveling wave, 24 see also Fig. 7 in Marburg. 19 This is the motivation of closer looking at the wave's amplitude in the duct and how it behaves with respect to the length of the duct. Figure 4 shows the sound pressure magnitude in terms of the location in the duct for different types of elements and different lengths. The frequencies are chosen such that in all four cases the number of degrees of freedom over the wavelength (number of nodes) remains the same. Note that P 1d and P 2c elements have twice as many nodal supports in x-direction as P 0 and P 1c elements. Furthermore, according to Table 1, the overall degree of freedom is twice as large for P 1d and P 2c elements. Thus, 4.25 and 8.5 elements per wavelength are counted. The normalized wave numbers are between kl = 6π and 96π for P 0 and P 1c elements and between kl = 12π and 192π for P 1d and P 2c elements. With respect to the element size, it is kh = 0.74π for P 0 and P 1c elements and kh = 1.48π for P 1d and P 2c elements, where h the element size. It is clearly visible that the sound pressure magnitude  It is quite obvious that the decay of the magnitude is much stronger for the P 0 and the P 1c elements. While for the latter, the solution at the left end seems to be close to the actual value of 0 cv 0 , it exhibits an initial offset for the P 0 elements. At the right end of the duct, the amplitude has decayed for the longest duct by almost 50% for P 0 elements and clearly even more for the P 1c elements. The decay is much smaller for the other two element types. It is even such that the linear discontinuous elements show a smaller numerical damping than the continuous quadratic elements. In general, the qualitative behaviour of the decay is very similar for the P 1c and the P 2c elements, also for the constant, i.e. P 0 elements and somewhat different for the P 1d elements. This will be discussed again at the end of this subsection. Figures 5 and 6 present the numerical error of the (complex) sound pressure solution in terms of the location in the duct for different types of elements and different lengths. This error function e(x) is evaluated as wherep(x) denotes the analytical solution according to Eq. (7) and p( x) represents the numerical solution at a field point x = (x, w/2, w/2) with w being half the width of the square cross section of the duct, cf. Fig. 1. The values of p are determined by a field point evaluation using the representation formula which is Eq. (4) for y ∈ Ω and, hence, c( y) = 1.  For investigating the numerical error, again, the frequencies are chosen such that, in all four cases, the number of degrees of freedom over the wavelength (number of nodes) remains the same. Two cases are investigated: Fig. 5 shows the local numerical error for the case that 34 and 17 elements P 0 , P 1c and P 1d , P 2c per wavelength are used, respectively. These numbers are much larger than the usual recommendations. Figure 6 shows the local numerical error for the case that 6.8 and 3.4 elements P 0 , P 1c and P 1d , P 2c per wavelength are used, respectively. This accounts for the common recommendations of discretization using boundary elements.
Analyzing the results shown in Fig. 5 reveals a surprise. While the numerical error of a few percent for the P 0 and P 1c elements and the n = 1 case could of acceptable order, when 34 elements per wavelength are used, numerical errors of up to 30% for P 0 elements and 45% for P 1c elements are really unexpected. According to what is known from the literature, that many elements per wavelength should be save to remain within an accuracy which is acceptable for technical problems, i.e. mostly between 1% and 10%. Such an accuracy is achieved when using P 1d and P 2c elements which are clearly providing accuracies of less than 1% even for the longest duct. Figure 6 presents the results of the case where 6.8 and 3.4 elements P 0 , P 1c and P 1d , P 2c per wavelength are used, respectively. This is what is usually done in practice. The common rule demands six boundary elements per wavelength. However, the author has successfully tested even fewer low order elements per wavelength in the past. 19,22,28 But staying with this number, it becomes clear that P 0 and P 1c elements may provide solutions with acceptable error for the shorter version of the duct. However, in case of the longest duct, P 0 elements end up with 100% error after traveling 80% of the length of the duct. P 1c elements even  manage to provide a numerical error of less than 100% up to only 40% of the length of the duct. The other two element types, i.e. P 1d and P 2c , show a qualitatively similar behaviour but at a somehow acceptable level. However, the local numerical error of the sound pressure at the right end is of orders of magnitude higher than at the left end. It is also substantially larger at the right end for the longer ducts and again, the local numerical error at the right end is most likely not acceptable for technical applications. It has been indicated above that P 1d elements are not always behaving exactly the same as the other three element types which are investigated. To further investigate this anomaly, the sound pressure and the local error are plotted for three different frequencies in Fig. 7. It is obvious that at 600 Hz, numerical damping cannot be observed. Actually, the sound pressure amplitude is increasing from the left to the right and this effect is even emphasized for the longer ducts at this frequency. At the frequency of 750 Hz, the sound pressure and its numerical error remain almost constant over the length of the duct and even for the ducts of different size. The behaviour at 900 Hz is very similar to what was observed for the other element types. A diminishing magnitude of the traveling wave and an increasing error are found in particular for the longest duct. This anomaly can be explained with the variation of the optimal location of the collocation points. The zeros of the Legendre polynomials account for the optimal location of the collocation points only for the Neumann problem,

1850018-11
i.e. for the case of the entire surface being acoustically rigid and, hence, Y = 0. The current example is a mixed problem since most of the surface is rigid but the right end is not. Therefore, the zeros of the Legendre polynomials are not necessarily to optimal location for the collocation points but they are close to the optimum. It was shown in former work of the author 22,28 that the optimum is moving with frequency. Apparently, the optimal location of the collocation points coincides with the zeros of the Legendre polynomials in the vicinity of the frequency of 750 Hz. Since this additional effect superposes the numerical damping effect, the otherwise very popular P 1d elements 21,23,29,31,32 will not be applied in what follows in this paper. When comparing the performance of the different element types, it becomes quite obvious that quadratic elements P 2c and linear discontinuous elements P 1d perform much better than the low order constant P 0 and continuous linear elements P 1d . This is due to the exponential convergence when increasing the polynomial degree. This means that doubling the number of degrees of freedom per wavelength is not just balancing the numerical error. The numerical error is much smaller when using higher order elements, at least up to order two. This has been described in former papers of the author 19,22,25,28 and it is obvious here.
So far, it can be concluded that there is no doubt that the length of the duct is a substantial parameter controlling the numerical error of the example. This is a contradiction to the widely accepted discretization rule of using a fixed number of boundary elements to control the discretization error and shows similarities to the pollution effect known from the finite element method.

Effect of local mesh refinement
In the previous subsection, only uniform meshes were investigated. However, when recalling the so-called pollution effect, it is necessary to check whether a local mesh refinement is able to compensate for the numerical error accumulated in another region of the model. For this, the model is divided into two regions of different mesh sizes. The model of l = 8l Ref = 27.2 m is chosen. One half of it is meshed rather coarsely with m s = 1 and the other half rather fine with m s = 6. Doing this, the entire mesh consists of 10101 boundary elements where in the coarse half, only 273 are found. The use of quadratic elements P 2c leads to 40426 nodes. Two configurations are investigated: In the first configuration, the coarse mesh is on the left side of the duct and the fine mesh on the right side whereas in the other configuration, it is the other way around. The highest frequency of 500 Hz is equivalent to 3.4 elements per wavelength for m s = 1 and 20.4 for m s = 6. This is equivalent to kh = 1.85 in the region of the coarse mesh and kh = 0.31 for the fine mesh. The overall wave number is kl = 80π.
The results are shown in Fig. 8. It can be clearly recognized that the wave is rather strongly damped in the region of the coarse mesh. In the case, where the wave is first traveling through the coarse mesh, it looks as if part of the error which is generated by the coarse grid is compensated by the finer mesh. However, the compensation seems to converge at a certain level, i.e. in the upper right subfigure, the error at 500 Hz in the middle of the duct is approximately 21% and it goes down to 15% at the outlet of the wave. At the same time the amplitude (upper left subfigure) increases from approximately 408 Pa to 419 Pa. In the case of the wave traveling through the fine mesh region first and then entering the region of the coarse mesh, the error in the middle of the duct is quite low and then increases with the same rate as shown in the left part of the upper subfigure. Even the sound pressure magnitude at 500 Hz at the right end is approximately 408 Pa and hence, almost the same as in the middle of the duct for the other case. Similarly, the numerical error is approximately 21% which is almost the same as in the other case. The conclusion for this example confirms the findings of the last subsection. It adds that a mesh refinement in the downstream region can somewhat compensate the accumulated error but, apparently, only to a certain level. Although being very similar meshes, a discretization error determined in the L 2 norm (as it is common for BEM) will be much higher if the wave travels through the coarse grid region first.

Phase error
It has been mentioned in the introduction that the phase error is the main reason for the pollution effect in the finite element method. In the previous subsections, the suspect candidate for the pollution effect in the boundary element method was numerical damping. Therefore, it will be examined in what follows how much influence on this results from dispersion. For this the longest duct is investigated at 1500 Hz, i.e. the highest frequency considered in this paper. This is equivalent to a normalized wavenumber of kl = 480π. As mentioned before, 240 waves are counted in the duct in the analytical solution at this particular frequency. The meshes considered are different for the three different element types. The low order polynomials, P 0 and P 1c are applied with the mesh parameter m s = 6. Quadratic elements P 2c use the mesh with m s = 3. This gives 6.8 elements per wavelength (kh ≈ 0.92) for the former and 3.4 for the latter case (kh ≈ 1.85). All of the models have approximately 39170 degrees of freedom. Coming back to the normalized wavenumbers, it is kl = 480π. Figure 9 gives an impression how the numerical solutions behave in comparison to the analytical solutions. All these subfigures show the real part of the traveling wave solutions for the left end and for the right end. Similar to the previous results, the solutions coincide very well at the left end and they deviate at the right end. It can be noted that both, numerical damping and phase errors are observed. Taking the results of the constant elements P 0 , the amplitude reaches 272 Pa (where 442 Pa are expected from the analytical solution) which is 61.5% of the original wave amplitude. The wave is shifted at the right end by about 16% of the wavelength. The results are somewhat worse for the linear continuous elements P 1c , for which the amplitude has already diminished to 175 Pa, which is only 39.6% of the original one. The wave is shifted by approximately 25% of the wavelength. Both, numerical damping and phase error are significantly higher for the linear continuous elements. Quadratic elements P 2c are much more accurate than the other two. While the amplitude of 395 Pa is still more than 89% of the original one, the wave is shifted by approximately 4.7% with respect to the analytical solution.
While the numerical damping, as it occurs for BEM, is not known for the finite element method applied to time-harmonic problems in acoustics, the phase error may be easily compared. For this, the one-dimensional problem is discretized by a simple and straightforward Galerkin method. The matrices can be found in the author's paper. 20 The meshes assume the same resolution as for the three-dimensional BEM mesh, i.e. 6.8 linear elements and 3.4 quadratic elements per wavelength. In these cases, the solution is much more dispersive. In the case of linear finite elements, 232.2 waves are counted instead of 240 and in the case of quadratic finite elements, 238.4 waves are counted. This results in a mismatch of approximately 7.8 and 1.6 waves, respectively. In comparison with these finite element results, the BEM phase error seems to be rather low and, at least, one order of magnitude lower than for similar cases in the finite element method.
Although small, the phase error in the boundary element solution exists and can be measured. It can be assumed that this error is due to a numerical dispersion. However, this assumption cannot be proved without an analytical investigation of the dispersion effect for boundary element solutions.

Conclusion
This paper has discussed the problem of the pollution effect for the boundary element method applied to acoustic problems. Other than known for the finite element method, the pollution effect for BEM is not primarily caused by a phase error but by numerical damping. The problem is very clearly exhibited for the long and slender duct problem. The author has observed indications that the pollution effect is less evident if the dimensions of the model are of the same order. However, this effect requires further investigation. A local mesh refinement for the long duct problem can only partly compensate the error which has been accumulated by numerical damping in another region of the boundary element model.
In addition to the amplitude error which is due to numerical damping, a small phase error has been observed. This may be caused by numerical dispersion. However, further investigations are required to confirm this.
It must remain unsolved at this point how relevant the pollution effect is for the practical use of the boundary element method. However, since it could be clearly shown that for the duct problems, the numerical error depends on the length and/or the number of waves within the problem, at least, for a long and slender computational domain. Therefore, users of the boundary element method should be very careful with the rule of thumb in which they are using a fixed number of elements per wavelength.
Finally, it was shown again that continuous linear boundary elements should not be applied for collocation. Whenever tested in this study, their performance was always worse than the performance of constant elements which were also not reaching the efficiency of discontinuous linear and continuous quadratic elements. This confirms earlier results of the author. 19,22,28 Although not tested, it is likely that, similar to the finite element method, 35 higher order elements perform even better than the elements tested here.