A Two-Dimensional Acoustic Tangential Derivative Boundary Element Method Including Viscous and Thermal Losses

In recent years, the boundary element method has shown to be an interesting alternative to the ﬁnite element method for modeling of viscous and thermal acoustic losses. Current implementations rely on ﬁnite-diﬀerence tangential pressure derivatives for the coupling of the fundamental equations, which can be a shortcoming of the method. This ﬁnite-diﬀerence coupling method is removed here and replaced by an extra set of tangential derivative boundary element equations. Increased stability and error reduction is demonstrated by numerical experiments.


Introduction
The isentropic acoustic wave equation is appropriate for modeling a vast variety of applications, but fails to give accurate solutions, especially if the acoustic domain is small and contains narrow gaps. In such situations, the effects of viscous and thermal dissipations must be considered. Two numerical approaches can be used for modeling of viscous and thermal acoustic dissipations, namely the finite element method (FEM) and the boundary P. R. Andersen et al. element method (BEM). FEM implementations for arbitrary geometries rely on a direct evaluation of the full linearized Navier-Stokes equations. This was initially proposed by Bossart et al., 1 realized by Malinen et al. 2 and later presented in different variations by Cheng et al., 3 Joly et al. 4 and Kampinga et al. 5 FEM can be computationally demanding, since the loss mechanisms require proper meshing of the associated boundary layers. Further, temperature and velocity are added as extra degrees of freedom. A numerically less costly method was recently proposed by Kampinga named the sequential linearized Navier-Stokes. 6 It is realized by removing some contributions, deemed negligible after an order of magnitude analysis and thus making it possible to separately solve uncoupled scalar viscous and thermal fields. While being more cost-efficient, it is an approximation and also requires careful meshing of the viscous and thermal boundary layers.
As opposed to modeling viscous and thermal losses with FEM, BEM is an interesting alternative, avoiding cumbersome meshing of boundary layers. Early BEM implementations relied on the work by Bruneau et al.,7 turned into BEM by Dokumaci 8,9 and Karra and Ben Tahar. 10 Their contributions were either having some restrictions or completely neglecting viscosity. Later, Cutanda Henríquez extended the ideas of Karra and Ben Tahar, 10 relying on the Kirchhoff's decomposition of the Navier-Stokes equations to include the effects of viscosity. 11 This approach was further developed into more general axisymmetric and threedimensional formulations. [12][13][14] While recent publications have shown the capabilities of this method, tackling large complicated problems, 15,16 the method might still have some shortcomings. Especially, the coupling of the Kirchhoff's decomposition equations is troublesome. In the current formulation of Cutanda Henríquez the coupling is handled through first-and second-order tangential derivatives evaluated with the use of finite differences. While this approach works for a large range of problems, it may be problematic for interior problems at low frequency, where the acoustic pressure is nearly uniform, potentially making finite differences inaccurate. The low-frequency implications were discussed in a recent conference paper by the authors, where a combined FEM and BEM approach was presented. 17 In the following, a new two-dimensional BEM approach is developed avoiding the use of first-and second-order finite-difference schemes in the coupling of the fundamental equations. This is shown possible by using an extra set of BEM tangential derivative equations for a more natural coupling of the equations. Finally, the new formulation will be evaluated through two simple academic test cases comparing it to the original finite-difference formulation.

Two-Dimensional Dissipative Boundary Element Formulation
Previous BEM implementations including the effects of viscous and thermal dissipations were based on the Kirchhoff's decomposition of the Navier-Stokes equations to achieve a form suitable for BEM discretization. The same approach is adopted in the following derivations. A large part of the present work relies on the formulation found in Ref. 12 and this work uses the same notation. In Refs. 12 and 7, the parameters (τ a , τ h , φ a , k a , k h and k v ) of the Kirchhoff's decomposition are discussed in more detail.
A 2D Acoustic Tangential Derivative BEM Including Viscous and Thermal Losses

Kirchhoff 's decomposition and boundary conditions
The final form of the decomposed Navier-Stokes equations is given in Eqs. (1)-(3), resulting in three equations of the Helmholtz form, where p a and p h are the so-called acoustic and thermal pressures, respectively. 18,19 Their sum p = p a + p h represents the total pressure. The vector v n is the rotational part of the velocity, also known as viscous velocity. In this paper, magnitudes written with an arrow vector on top are vectorial fields, while bold letters indicate matrices. While p a resembles an actual acoustic pressure wave, the equations containing p h and v v can be considered as heavily damped wave equations. k a , k h and k v denote the acoustic, thermal and viscous wavenumbers, respectively: Equations (1)-(3) are coupled at the boundary through boundary conditions. It is reasonable to assume isothermal boundary conditions. The heat capacity is very much higher at the boundary than in the fluid that temperature variations at the boundary can be neglected. An isothermal boundary condition can be expressed in terms of the acoustic and thermal pressures and two frequency-dependent terms τ a and τ h relating the acoustic and thermal pressures to the temperature fluctuations T . This results in a frequencydependent thermal boundary layer forming at the boundary, with a thickness ranging from micrometers to millimeters in the audible frequency range. Isolating p h in the isothermal boundary condition, Eq. (4), is discretized as, which will be used as the first condition to couple the equations of Kirchhoff's decomposition. A second boundary condition states that the fluid will tend to stick to the boundaries (noslip condition). A viscous boundary layer will form, having similar thickness as the thermal boundary layer. This boundary condition is expressed as where v b is the boundary velocity and φ a and φ h are constants depending on physical parameters and frequency. For the further development, it is convenient to describe the no-slip conditions in a local boundary coordinate system: with the subscripts n and t denoting the normal and tangential components to the boundary. P. R. Andersen et al.

Boundary element discretization
The boundary element implementation for regular isentropic acoustic problems starts with the integral form of the harmonic Helmholtz equation: where C(P ) is the integral-free term. The first term on the right-hand side is the doublelayer potential and the second term is the single-layer potential. P is a calculation point, Q is an integration point on the generator, R = |Q − P | is the distance between calculation and integration points and G(R) is the fundamental solution in free space. An e iωt time convention is assumed. It should be noted that inclusion of a source term is omitted for simplicity, but could be added as an extra term. Meshing and collocation of the Helmholtz integral equation in matrix form is given by where C is a diagonal matrix containing the integral-free term, H and B are the discretized double-and single-layer potential matrices, respectively, and p and ∂p ∂n are vectors. A more compact notation with A = H − C will be used in the following development. Since the Kirchhoff's decomposition produces equations formally equivalent to the Helmholtz wave equation, it is straightforward to apply BEM to Eqs. (1)-(3) using collocation and replacing the isentropic wavenumber k with k a , k h and k v : where Eq. (3) is split into its Cartesian components forming Eqs. (12) and (13). Note that we here and in the following limit the formulation to two dimensions (x, y) for simplicity. In two dimensions, the fundamental solution is given by is the modified Bessel function of the second kind of order j. The argument of the fundamental solution is the imaginary unit i and the wavenumber k can be either acoustic, thermal or viscous.

Coupling of equations
The previous boundary element implementations with losses make use of finite differences to couple the fundamental equation set. The new approach presented here utilizes the boundary element itself to estimate tangential derivatives. Taking the tangential derivative thus allowing for the evaluation of the tangential derivative of the pressure on the collocation points from boundary pressures p(Q) and their normal derivatives ∂p(Q) ∂n . Equation (14) contains second derivatives of the fundamental solution. A similar situation arises when dealing with irregular frequencies for exterior problems using the Burton-Miller formulation, which contains double normal derivative kernels that are said to be hypersingular and require special treatment. 20 However, the integrals arising from the tangential derivative equations are Cauchy principal value (CPV) integrals. The kernel description can for example be found in the work by Gallego and Martínez-Castro. 21 The viscothermal implementation presented here relies on an adaptive integration scheme for the evaluation of near-singular integrals arising in the case of narrow gaps, 22 but also for the evaluation of the CPV integrals. Initial convergence studies have shown that this approach is feasible for the evaluation of the tangential kernels.
A second set of discretized equations can now be formed containing the tangential derivative kernels: It is possible to describe the no-slip condition in terms of the discretized equations by isolating ∂p a ∂n , ∂p h ∂n , ∂p a ∂t and ∂p h ∂t in Eqs. (10), (11), (15) and (16), respectively. The boundary condition equations, i.e. Eqs. (6) and (7), can then be described by Isolating p h in the isothermal boundary condition (4), the thermal pressure can be eliminated from Eqs. (19) and (20). Finally, the boundary conditions of the discretized equations can be stated as where

Null-divergence of the viscous velocity
The viscous velocity forms a rotational vector field, meaning that ∇ · v v = 0 must be fulfilled. This is ensured through the null-divergence in a normal and tangential coordinate system (note that the tangential direction only has a single component in a 2D problem): which will require a relation between local and global quantities on the boundary, achieved in an all-geometry approach for the individual discrete nodes by where • is the element-wise Hadamard product (see Appendix A) and n x , n y , t x and t y are the Cartesian components of the normal and tangential vectors at each node. The transform of the viscous velocity back to the Cartesian coordinates is The goal is now to establish an expression for each of the terms in Eq. (26) using the discretized equations, containing only local normal and tangential components of the viscous velocity. This requires taking the derivative of Eqs. (27) and (28) with respect to the normal and tangential components, respectively. Doing so yields By isolating ∂vv,x ∂n and ∂vv,y ∂n in Eqs. (12) and (13) and substituting into Eq. (31), we obtain Using Eqs. (29) and (30) to describe the viscous velocity in terms of the normal and tangential components, Eq. (33) becomes Applying the properties of the Hadamard product for diagonal terms, the terms in Eq. (34) can be rearranged so that where The result of Eq. (35) will be used later to ensure that the null-divergence condition is fulfilled. Proceeding with the second term in Eq. (26) (17) and (18) and substituting the result into Eq. (28), an expression for the tangential derivative of the tangential viscous velocity is found: The Cartesian viscous quantities can be transformed into the local form by combining Eqs. (29), (30), (12) and (13) into Rearranging the terms and using the Hadamard product properties, the following equation is obtained: where The result of the normal and tangential derivatives of the local viscous velocity components, Eqs. (35) and (40), is substituted into the local form of the null-divergence, Eq. (26), forming where V and T are defined as The result of Eq. (43) will be used in the following to ensure that the viscous velocity has null-divergence.

System of equations and boundary quantities
It is now possible to establish a system of equations containing boundary normal and tangential velocities and the acoustic pressure. Isolating v v,n and v v,t in Eqs. (43) and (23), respectively, and substituting them into Eq. (22) yields the final system of equations: which allows for applying boundary velocity conditions and solving for the nodal acoustic pressures p a . The total pressure can be described as the sum of the acoustic and thermal pressures, by assuming isothermal boundary conditions: The normal and tangential viscous velocities can be found by rearranging Eqs. (22) and (23), so that The evaluation of field points can be done by following the approach in Ref. 12. The expression using the finite-difference approach originally developed by Cutanda Henríquez, is reproduced here for comparison: where DT 1 and DT 2 denote the first and second tangential surface finite-difference matrices, respectively. Equation (50) is developed by taking the tangential derivative of the boundary condition, Eq. (7), whereas the new approach presented in this paper achieves this goal through a set of tangential boundary element matrices, Eq. (38). As a consequence, the use of finite difference but also the evaluation of second tangential pressure derivatives are avoided.

Test Cases
Only few analytical solutions are available for testing of acoustic viscothermal implementations. In the first example, the implementation will be compared with an analytical solution for an oscillating cylinder. 23 This example considers only viscosity, but this is not an important limitation since it is the implementation of viscosity that poses a challenge in BEM with losses. Thermal losses are easily removed from the formulation by neglecting the terms containing thermal quantities in E n and E t . The second example is a narrow duct. This example highlights a problematic behavior that might be present in the finite-difference implementation for some cases.

Discretization
One requirement when dealing with the tangential derivative integration kernels is the need for C 1 continuity at the collocation points. This can be achieved by using discontinuous elements. Discontinuous elements are also known to perform well for boundary element implementations. 24 In order to evaluate the effect of the new tangential derivative formulation on equal footing, the finite-difference formulation and the new formulation will use the same boundary element matrices created using discontinuous quadratic elements. In the second example, however, the finite-difference implementation uses the full approach found in Ref. 12, with regular continuous quadratic Lagrangian elements.

Oscillating infinite cylinder
The infinite cylinder geometry is shown in Fig. 1. The radius R is set to 1 m and the amplitude of oscillation is 1 m/s. A relative error measure will be defined as, with the index j denoting the individual nodal solutions, N n is the total number of nodes and p j,ref is the corresponding reference solution. Figure 2 shows the convergence of the pressure  Higher frequency solutions could be obtained by using for example the combined Helmholtz integral equation formulation, 25 but this is considered beyond the scope of the paper. The normal viscous velocity is usually much smaller than the corresponding tangential component and more difficult to compute accurately. The error in pressure is similar for the two implementations above approximately 100 degrees of freedom, but higher mesh densities seem to favor TD-BEM. The convergence situation is slightly different in the case of the normal viscous velocity. In this case, TD-BEM experiences a larger error for low degrees of freedom, as compared to FDD-BEM, but for higher mesh densities TD-BEM is again a better choice. It is not fully understood why TD-BEM experiences a higher error at low mesh densities. A plausible explanation could be related to the use of BEM to estimate the tangential derivatives of the viscous velocity, which might require more elements to be captured appropriately.

Traveling wave in narrow duct
As a second test case, a traveling wave in a narrow duct is considered. Approximate analytical solutions for this case exist, that are applicable in the extreme case of very narrow or wide tubes. To highlight a problematic behavior of FDD-BEM, the case of slightly overlapping boundary layer is studied, therefore an FEM solution with a very fine mesh is chosen as a reference. The geometry is shown in Fig. 3. The left end of the duct is excited with a normal boundary velocity with an amplitude of 1/(ρc), where ρ is the density and c is the speed of sound for air. The other end of the duct is fitted with an impedance of value ρc.
For isentropic acoustic problems this would resemble a fully absorbing boundary, but this is not necessarily true when the duct is narrow and the boundary layers fill up most of the domain. The reference FEM solution is implemented in COMSOL Multiphysics with the same boundary conditions and approximately 300,000 degrees of freedom, using an equally spaced structured quadrilateral mesh. The boundary layers are covered with about 50 elements in the direction perpendicular to the boundary. To ensure similar conditions in the Kirchhoff's decomposition and the COMSOL implementation, simulations were carried out assuming an ideal gas and no flow. The length of the duct, L, will be five times the viscous boundary layer thickness, δ v , which for air is 2.21/ √ f (mm), 18 with f being the frequency of the oscillation. The pressure magnitudes are plotted in Fig. 4 along the boundary in the x-direction for three different duct heights, h. The heights are half, equal and two times the viscous boundary layer thickness. Simulations were carried out at 2 kHz, meaning that the wavelength is much larger than the computational domain, resulting in a relatively uniform pressure. It is seen how FDD-BEM fails to give accurate solutions, especially when the height is equal to the viscous boundary layer thickness. This result is invoked by using slightly smaller elements on the top and bottom parts of the duct compared to the sides (30 elements on the top and bottom parts and five elements on the sides). To further investigate this behavior, the worst case with h = δ v is tested with different numbers of elements. The sides of the duct are fixed to five elements on each boundary, and the top and bottom boundaries keep the same element count, but with varying size. The error measure of   Eq. (51) is used, but with the FEM solution as a reference. The results are shown in Fig. 5, where the total number of elements in the geometry is ranging from 30 to 90. Surprisingly, FDD-BEM shows an initial low error but a rapid error increase followed by a low error at 60 elements, recovering the high error at higher element counts. In the 60-element geometry, all elements, including those on the lid, are of equal size. This means that FDD-FEM is more likely to behave with low errors, for a controlled environment with elements of equal size. This is also the case for the oscillating infinite cylinder. It appears that FDD-BEM might result in instability problems if the element sizes are dissimilar and the pressure is nearly uniform, as seen in the example. On the other hand, the new TD-BEM shows a much more stable error response to the change in element size.

Conclusions
A new approach to the numerical implementation of the coupling of the Kirchhoff's decomposition describing the propagation of sound waves with viscous and thermal losses is shown possible through development of an extra set of tangential boundary element equations. The new formulation completely removes finite differences as a coupling method. This leads to an improvement over the previous BEM implementation with losses. The method is compared to an earlier implementation through a simple convergence study of an infinite cylinder, showing that for higher element counts the new method can be expected to give slightly lower errors. Coarser meshes might lead to a larger error in the viscous velocity computation.
In the second example, a narrow duct, the original finite-difference implementation presents some instability issues, leading to an unwanted behavior when element sizes are different.
On the other hand, the presented new approach shows a much more stable error behavior.
There exists a hope that the new coupling concept can be extended to three dimensions, where the use of tangential finite-difference derivatives is more cumbersome.