Quantiﬁcation of Numerical Damping in the Acoustic Boundary Element Method for Two-Dimensional Duct Problems

Spurious numerical damping in the collocation boundary element method is considered for plane sound waves in two-dimensional ducts subjected to rigid and absorbing boundary conditions. Its extent is quantiﬁed in both conditions based on a damping model with exponential decay, and meshes of linear and quadratic continuous elements are studied. An exponential increase of numerical damping with respect to frequency is found and the results suggest an upper bound for given element-to-wavelength ratios. The quantiﬁcation of numerical damping is required for evaluation of meshes covering a large number of waves, and real damping phenomena can be prevented from overestimation.


Introduction
The boundary element method (BEM) is a popular numerical tool for solving linear timeharmonic acoustic problems arising in many technical and scientific applications. [1][2][3] Using the BEM, modeling is limited to the boundary of the wave-carrying fluid domain. Hence, boundary element discretizations involve fewer degrees of freedom than corresponding discretizations obtained by the finite element method (FEM). Moreover, the implicit satisfaction of the Sommerfeld radiation condition is particularly advantageous for exterior problems compared to the FEM. 4 Concerning the linear time-harmonic FEM, it is widely known that it suffers from the so-called pollution effect 5 -a cumulative numerical error growing with the number of waves in the domain. 6 Different from the common discretization error, the local error caused by the pollution effect can only be reduced by a global mesh refinement. The pollution and the related dispersion effect, i.e. the numerical error introduced by the wavenumber in the discrete setting, have been rigorously discussed in literature. 5,[7][8][9] Moreover, error estimates 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.
have been presented [10][11][12] and modified finite element formulations have been suggested to address this issue. [13][14][15][16] However, much less work has been carried out in the past regarding a similar effect that occurs when using the BEM. This effect will be denoted as numerical damping in what follows. It has been studied in the last years for acoustic interior problems, [17][18][19] and first approaches towards its quantification exist. 20 The present paper is built upon these works, and it is organized as follows.
Section 2 serves as a short introduction to numerical damping in the BEM and demonstrates its implications to the reader. In Secs. 3 and 4, two methods are presented for the quantification of numerical damping arising for plane sound waves in two-dimensional ducts subjected to rigid and fully absorbing boundary conditions, respectively. The first method is based on the full width at half maximum (FWHM) of resonance peaks in the course of frequency response curves. In the second approach, the sound pressure amplitude decay of a traveling wave is related to an analytical solution subjected to losses. Both methods assume a qualitative equivalence between numerical damping and linear viscous damping. They are applied in Sec. 5, and the extent of numerical damping is studied over frequency for different meshes of linear and quadratic continuous boundary elements. The results are discussed in Sec. 6, and the paper concludes by pointing out the benefits, applicability and the limitations of the presented approaches for quantifying numerical damping.

Brief Review of Numerical Damping in the Boundary Element Method
In the following, the phenomenon of numerical damping is introduced using the example of an acoustic interior problem. A plane sound wave in a closed two-dimensional duct of length l = 3.4 m and width w = 0.2 m is considered, as shown in Fig. 1. It is filled entirely with air with a density of ρ = 1.3 kg/m 3 , and the speed of sound is c = 340 m/s. A sound wave is generated due to a harmonic excitation by a particle velocity of v 0 = 1 mm/s at x = 0 and is fully reflected at the end x = l of the rigid duct. The system is free of dissipation, and, hence, resonances are expected at the integer multiples of 50 Hz. At those frequencies, standing waves with infinite amplitudes are formed, and the phase jumps equal 180 • . Comprehensive studies on error estimates of this problem are available in the literature. 21,22 The threedimensional equivalent of this model is defined as a benchmark problem of computational acoustics. 23  The frequency response of the system is studied using the BEM in conjunction with a collocation discretization. The mesh consists of 32 one-dimensional boundary elements with quadratic Lagrangian pressure approximation, of which 15 elements are employed for both edges parallel to the direction of wave propagation, respectively. Regarding the numerical integration, an adaptive scheme is used, where the domain of singular integrals is subdivided into smaller intervals. 24 The number of Gaussian points is chosen according to a relative distance function. 25 The frequency response is evaluated up to 520 Hz. At 500 Hz, a standing wave composed of 10 half-waves emerges, corresponding to a ratio of three elements per wavelength. A numerical error of less than 5 % in the L ∞ -vectornorm can be expected. 21 In Fig. 2, the sound pressure level p L and the phase angle φ at the center of the duct, i.e. x = 1.7 m, y = 0.1 m, are shown. Resonances only occur every 100 Hz since the chosen point is associated with a node of the modes at the uneven integer multiples of 50 Hz. Figure 3 provides close-ups of the frequency response curves at around 100 and 500 Hz, respectively. Obviously, the peak at 500 Hz is wider and less pronounced than the peak at 100 Hz. Moreover, the phase jumps at the low-frequency resonances are sharp, while at 500 Hz, they are rounded.
Those phenomena clearly indicate damping that is introduced numerically, since the underlying physical system is free of dissipation. Furthermore, the extent of numerical damping seems to grow with frequency. However, numerical damping does not only occur for the case of a fully reflected sound wave as presented here, but it is also encountered for traveling waves. 19 The upcoming sections present methods to quantify the numerical damping in both scenarios.

Full Width at Half Maximum for Interior Acoustic Problems
The occurrence of artificial numerical damping can lead to an overestimation of real damping phenomena. Moreover, in the case of large computational domains accommodating a large number of waves, excessive numerical damping leads to unacceptable errors. Hence, methods for the quantification of its extent are desirable. For this purpose, it is assumed that numerical damping behaves similarly to viscous fluid damping, i.e. it is proportional to the time derivative of the sound pressure. The method presented herein is based on the FWHM of the sound pressure peak.
The FWHM method is widely applied in experimental damping estimation by means of frequency response analysis and extensively discussed for mechanical systems with a single degree of freedom. 26,27 The starting point for its derivation in the context of acoustic interior problems is the scalar, damped, homogeneous wave equation in one dimension, i.e.
Here, µ denotes the unknown damping coefficient. Separating the variables, the sound pressure is decomposed into solely time-and position-dependent functions, respectively, and it is written as p(x, t) = g(x)p(t). For the purpose of damping estimation from a frequency response, only the frequency ranges around resonances are of interest. At these frequencies, standing waves form, that can be described by a sinusoidal position-dependence under the assumption of light damping. Hence, in which the wave number k = 2πf /c is introduced. Double differentiation, inserting p(x, t) into (3.1), and imposing a time-harmonic excitation on the right-hand side yields 3) describes the sound pressure in the time domain at a discrete position in the fluid and is similar to the differential equation of a mechanical system with a single degree of freedom. In (3.3), ω 0 is the angular eigenfrequency of the undamped system. A e and ω e denote the amplitude and angular frequency of the excitation, respectively. In steadystate, the time-dependent pressure function readsp(t) =p cos(ω e t + φ) with the pressure amplitudep and the phase shift φ with respect to the excitation. Insertingp(t) into (3.3) and rearranging giveŝ for the sound pressure amplitude and its time derivativeṗ. In the case of light damping, the resonance frequency accords approximately with the eigenfrequency (ω e ≈ ω 0 ), and the time derivative can be written asṗ Furthermore, two angular frequencies ω I and ω II are defined at which the amplitude drops to √ 2/2ṗ r , which corresponds to a decrease of 3 dB. Relating that expression to (3.4) and solving the quadratic equation leads to If ω I and ω II are given from the frequency response curve, the unknown damping coefficient µ can be determined by The application of the FWHM method to acoustic interior problems is schematically depicted in Fig. 4. However, awareness of the underlying assumptions is important when using the method: (a) The damping behaves qualitatively as viscous damping, i.e. it is proportional to the time derivative of the sound pressure. (b) The extent of damping is small. Hence, (3.2) is valid. Moreover, the resonance frequencies accord with the eigenfrequencies of the undamped system. (c) Damping between ω I and ω II is independent of frequency.

Quantification of Numerical Damping for a Traveling Sound Wave
The phenomenon of numerical damping is not only observable in the case of resonances but also for a traveling sound wave. The recent paper 18 demonstrates the amplitude decay of an one-dimensional sound wave in a duct with an absorbing boundary condition at the outlet. Since the underlying physical system is free of dissipation, the decay clearly indicates numerical damping in the BEM. In the following, a method for the quantification of numerical damping occurring for a traveling wave will be presented. Assuming again that it qualitatively behaves similarly to viscous damping, the boundary element solution will be related to an analytical solution to obtain the unknown damping coefficient. The derivation of the analytical solution again starts with the one-dimensional wave equation (3.1). Introducing the harmonic timedependence e −iωt leads to the damped Helmholtz equation, which is given as where the complex wavenumberk has been introduced. Note that in contrast to Sec. 3, p(x) denotes the solely position-dependent complex sound pressure. Equation (4.1) is accompanied by Robin boundary conditions relating the structural particle velocity v s to the sound pressure on the boundary, i.e.

iωρ
where Y (x) is the boundary admittance. A traveling wave in a duct of length l is excited by a particle velocity v 0 at the rigid left end and absorbed at the outlet on the right-hand side, i.e. v s (0) = v 0 , Y (0) = 0, v s (l) = 0, Y (l) = 1/ρc. The corresponding boundary value problem is solved by integrating (4.1), and the analytical solution emerges as The extent of numerical damping will be quantified by relating the amplitude decay of (4.2) to the decay of the boundary element solution. However, the boundary element solution is additionally subjected to a discretization error that needs to be excluded in the following study. Considering that discretization induces a local numerical error independent of the number of waves in the computational domain, the numerical solution of a wave traveling along one dimension in a uniformly meshed setting will yield a spatially constant discretization error. In contrast, the numerical damping subjects the amplitude of the boundary element solution to a spurious decline with respect to the number of waves. an arrayp be . It clearly exhibits a decay along the whole length l = 100 m of the duct while the rate of decay decreases towards the propagation direction of the wave. The amplitude of the analytical solutionp a (µ e ) with an exemplary damping coefficient µ e = 0 is also given in Fig. 5. Note that viscous damping, as modeled here for the analytical solution, results in an exponential decay e −αx with the attenuation coefficient However, the exponential behavior is hardly recognizable in Fig. 5.
With the analytical solution at hand, the unknown numerical damping coefficient µ corresponding to a particular boundary element model is estimated based on the following idea: The constant offset, i.e. the mean difference betweenp be andp a (µ e ), is identified as the discretization error. This offset, and hence the discretization error, is excluded in the process of finding µ. Thereby, the spurious oscillations of the boundary element solution are also ruled out. The final coefficient µ is found such that the absolute difference between the boundary element solution and the adjusted analytical solution is minimal. The algorithm for finding the numerical damping coefficient for a traveling wave and the emerging minimization problem are given below.
(1) Calculate the boundary element solution for a given mesh and evaluatep be ∈ R n along n interior points along the duct.

S. K. Baydoun & S. Marburg
However, in view of the curves in Fig. 5, it has to be noted that the decay characteristics of the analytical and the boundary element solutions differ significantly. While the slope of the numerical solution clearly flattens, the slope of the analytical sound pressure is almost constant in the considered case. As a consequence, the resulting damping coefficient of the above algorithm is dependent on the traveled distance of the wave and can be interpreted as a spatially averaged damping coefficient. Shorter captured distances will yield higher damping coefficients and vice versa.

A Benchmark Study: Numerical Damping in Boundary Element Methods for Two-Dimensional Acoustic Duct Problems
Having introduced two approaches for quantifying numerical damping, its extent is studied in what follows for acoustic waves in two-dimensional ducts. Linear and quadratic continuous boundary elements and different mesh sizes are investigated. For the sake of comparability, the mesh size h refers to the distance between adjacent nodes in the direction of wave propagation. One quadratic and two linear elements are employed across the width of the duct, respectively, for all meshes. Hence, for a given h, meshes of linear and quadratic elements yield the same number of degrees of freedom. Numerical damping of fully reflected sound waves is evaluated using the same fluid properties, geometrical setup and conditions as in Sec. 2. Figure 6(a) shows the damping coefficients at 10 resonance peaks between 50 Hz and 500 Hz. They increase with increasing   mesh size and numerical damping is generally higher for linear elements of same mesh size in the considered frequency range. Moreover, linear and quadratic elements exhibit different slopes with respect to frequency. However, the slopes seem to be independent of the mesh size. The linearity of the graphs in Fig. 6 indicates an exponential relation µ ∝ f β between the numerical damping coefficient and the frequency with an exponent β, which depends on the element type. The curves in Fig. 6(a) correspond to β 1 = 1.5 for linear and β 2 = 3.5 for quadratic elements, respectively.
Similar analyses are performed for the case of a traveling wave and presented in Fig. 6(b). The damping coefficients are obtained as described in Sec. 4 by relating the sound pressure decay over the course of l = 100 m to the analytical solution. At 500 Hz, this corresponds to 150 waves in the duct. Note that it is hardly possible to obtain converged damping coefficients for finer meshes of quadratic elements in the lower frequency range due to weakly pronounced decays. These meshes would require longer traveled distances than l = 100 m to be captured.
A comparison of both evaluation methods is shown in Fig. 7. As already discussed in Sec. 4, the results obtained for the traveling wave are dependent on the choice of the captured decay length l, and they are to be interpreted as spatially averaged damping coefficients. Therefore, only a qualitative comparison of the two methods regarding the relation of numerical damping to frequency, element type and mesh size is possible. In this regard, the results in Fig. 7 show a good agreement.
The increase of numerical damping with respect to frequency in Figs. 6 and 7 is not surprising, as the approximation quality generally deteriorates with an increasing number    of waves for a given mesh. Therefore, it is common to study numerical errors for constant ratios of elements per wavelength in the context of wave simulations. This has been justified in a both descriptive 21 and normative 28 manner. Figure 8 shows the numerical damping for eight, 12 and 16 linear as well as four, six and eight quadratic elements, respectively, evaluated with both methods. Note that eight linear elements per wavelength correspond to the same number of degrees of freedom as four quadratic elements and so forth. Clearly, the quadratic elements perform significantly better in the considered frequency range. Moreover, it can be seen that for the constant element-to-wavelength ratio, the numerical damping is not constant but decreases with increasing frequency. Consequently, an upper bound for numerical damping in the whole frequency range could be obtained by evaluating the numerical damping at the lowest frequency of interest. For example, the widely applied rule of six quadratic boundary elements per wavelength 21 is subjected to a numerical damping of µ < 1.22 · 10 −6 s/m 2 for frequencies above 50 Hz and the given duct geometry. According to (4.3), this value corresponds to an attenuation coefficient ofα < 2.1 · 10 −4 1/m, and the values for the other ratios of elements per wavelength are given in Table 1. The conversion of attenuation coefficients to other units such as dB/m and Np/m is described in a monograph by Knudsen, 29 though their notation is not consistent in the literature. 30 Typical values in air 31 and sea water 32 in standard conditions at low frequencies are given in Table 2. It can be seen that the attenuation induced by numerical damping in the BEM is of the same order of magnitude as attenuation in air in the considered frequency range. The damping ratio ζ is another widely used dimensionless measure for damping, in particular, for mechanical systems with a single degree of freedom. For acoustic interior problems  at resonance, it can be determined from ζ = µc 2 /2ω 0 . The above example of a damping coefficient µ = 1.22 · 10 −6 s/m 2 at 50 Hz corresponds to a damping ratio of ζ = 2.4 · 10 −4 .

Conclusion
The phenomenon of spurious numerical damping in the acoustic BEM has been shown phenomenologically for plane waves in two-dimensional ducts subjected to rigid and fully absorbing boundary conditions. Under the assumption that numerical damping can be characterized similarly to linear viscous damping, simple methods for quantifying its extent in both conditions have been presented. In the case of acoustically rigid conditions, the sound pressure wave is fully reflected, and standing waves form at the resonance frequencies. The FWHM method is then applied in order to deduce the numerical damping coefficients at the sound pressure peaks. Secondly, in the case of full absorption at the outlet, numerical damping is quantified based on comparisons to analytically determined decay characteristics of damped traveling waves. Using these methods, numerical damping has been evaluated for different meshes of linear and quadratic continuous boundary elements. In the considered frequency range, the unknown numerical damping coefficient increases exponentially with respect to the frequency. Its extent has also been studied for constant ratios of elements per wavelength, and the results suggest that an upper bound of numerical damping may be found for given ratios. Thereby, real damping phenomena, such as the structural damping of fluid-filled cavities due to sound radiation, can be prevented from being overestimated in the future. The quantification of numerical damping can also be essential in the case of large computational domains accommodating a large number of waves. However, while the presented methods are suitable for evaluating the appropriateness of a given mesh in regard to artificial damping, there are clearly some limitations. First and foremost, the presented results revealed that the rate of decay due to numerical damping differs qualitatively from an exponential decay. Hence, the mathematical description of numerical damping similar to viscous damping provides a spatially averaged measure. Numerical damping as defined here can neither be estimated uniquely in terms of element sizes nor in terms of element-to-wavelength ratios. Instead, it is a property of each individual boundary element model as a whole involving its geometry. Therefore, the results in this paper can only provide qualitative insights to the relation of numerical damping to mesh sizes and element types. The question, to what extent their actual values are valid for other geometries, is beyond the scope of this paper and requires further research. However, preliminary analyses of ducts of different lengths and widths have shown similar slopes of the damping coefficient with respect to frequency as given in Figs. 6 and 7 for both, linear and quadratic boundary elements.
Three-dimensional ducts have been qualitatively investigated in recent papers 18,19 of one of the authors. Though no actual values were determined there, it was also shown that a refinement of the mesh decreases the effect of numerical damping, and that its extent increases with frequency. Both quantification methods presented here could be applied straight-forwardly to three-dimensional duct problems in order to obtain the damping coefficients.
Since the quantification method for traveling waves relies on a comparison to an analytical solution, it is not generally applicable. On the other hand, the FWHM method provides a general tool for acoustic interior problems subjected to resonances.
Although exterior problems account for the main field of application of the BEM, they have not been considered in this paper. So far, the authors did not identify the occurrence of numerical damping for exterior problems. However, further research is required to confirm this.
Moreover, the studies in this paper were limited to the boundary element method with collocation discretizations. Preliminary investigations of a duct problem using the noncommercial software BEM++, 33 which is based on the Galerkin method, also showed indications of numerical damping, though the extent is smaller than for the collocation method. The investigations in the recent paper 18 of one of the authors indicate that the origin of the numerical damping stems from the complex-valued fundamental solution of the Helmholtz equation. It was shown that no numerical damping is introduced when using real-valued fundamental solutions. However, boundary element formulations associated with purely real fundamental solutions can lead to numerical instabilities. 18 Finally, it has to be noted that truly normative estimations of numerical damping in the BEM can only be carried out mathematically, as previously done for the pollution effect in the FEM. 5 Meanwhile, the presented results provide the boundary element community with useful qualitative and quantitative insights into this phenomenon.