In the acoustics domain, simulating and analyzing the propagation of sound waves through complex structures is essential. Frequency sweep calculations play a pivotal role in understanding sound wave behavior across various frequencies, crucial for designing effective acoustic barriers and noise reduction devices. Traditional methods for these calculations often utilize the finite element method (FEM) or the boundary element method (BEM). It is common knowledge that BEM is frequently utilized to address acoustic issues because of its superior accuracy and simplicity in mesh creation [1–4]. For external acoustic problems, it naturally satisfies the Sommerfeld radiation condition at infinity [5–9]. Conventional sound pressure calculations are typically optimized for specific frequencies, limiting their applicability across a broader frequency spectrum. To address this limitation, broadband analysis is introduced to provide results over a wider range of frequencies [10]. Nevertheless, in broadband analysis, the frequency band within a given range is segmented, and the coefficient matrix along with the boundary elements system equation is recalculated at every distinct frequency point. This procedure results in substantial computational expense. This study presents an efficient approach for rapid frequency sweep calculations in broadband acoustics.

When performing broadband acoustic analysis with the BEM, the frequency dependence of the H n ( 1 ) ( k r ) results in coefficient matrices that change with frequency [11, 12]. Consequently, for large-scale issues, the approach is particularly time-consuming because the boundary element system equations must be recalculated for every discrete frequency point. Researchers came up with a number of techniques to speed up the broadband analysis computational process to address this problem [13–16]. A three-dimensional (3D) axisymmetric multifrequency acoustic analysis technique called the Linear Frequency Interpolation Technique (FIT) was proposed by Vanhille et al. [17]. In this approach, the second-order isoparametric segment of the Helmholtz Integral Equation (HIE) is discretized. Despite its effectiveness, FIT requires substantial storage capacity. Li [18] addressed the issue by separating the frequency terms from the damping function through power series expansions of sine and cosine functions into algebraic polynomials, which significantly reduced computational time for multifrequency problems. Similarly, Zhang et al. [19] tackled multifrequency issues using BEM and incorporated series expansion techniques in their calculations, which offer high accuracy and computational efficiency [20]. The advantages of using series expansion are precision and computational economy [21–23]. In this study, the frequency-related and frequency-unrelated components of the product function in the Boundary Integral Equations (BIE) are separated using a Taylor series expansion.

In the realm of acoustic computations using the BEM, coefficient matrices exhibit characteristics of being asymmetric, full-rank, and dense [24, 25]. These attributes contribute to diminished computational efficiency, particularly evident in scenarios involving large-scale problems [26, 27]. Model Order Reduction (MOR) emerges as a viable solution to this challenge [28–31]. Among the most prominent MOR techniques [32] is Proper Orthogonal Decomposition (POD) [33, 34]. However, the quality of the simplified model is not guaranteed, as it depends on the representativeness of the snapshots (or primary frequencies) selected during the POD process [35]. The second-order Arnoldi (SOAR) algorithm for two-dimensional (2D) linear systems, which was introduced by Bai et al. [36, 37], has garnered interest from numerous scholars. The full-order model (FOM) gets mapped onto the projection space. SOAR facilitates the acquisition of an orthogonal basis within the projection space and the construction of a Reduced-Order Model (ROM) that preserves the characteristics of the initial model. Furthermore, the SOAR method finds extensive application in structural acoustic analysis, second-order dynamical system modeling [19], and quadratic eigenvalue problems [37, 38].

Singular Helmholtz boundary integral equations may fail to yield unique solutions when applied to exterior boundary value problems. To address this issue, two primary approaches have been proposed [39–42]. The Combined Helmholtz Integral Equation Formula (CHIEF), cited in Ref. [43], effectively addresses this issue by introducing additional HIE within the internal domain [44, 45]. The resulting overdetermined system of equations can be solved using the method of least squares. However, determining the optimal number and placement of internal points, particularly for issues of high frequency, remains challenging. The Burton-Miller method [46] presents another useful strategy for handling non-unique solutions, offering a linear formulation of the Classical Boundary Integral Equation (CBIE) and its associated Normal Derivative Boundary Integral Equation (NDBIE). If the boundary is nonsmooth, NDBIE becomes hypersingular, requiring special numerical treatment. This study utilizes Cauchy principal value integrals and Hadamard finite part integrals for handling singular integrals.

To enable the application of the Boundary Element Method (BEM) in 2D acoustic computations across a wide frequency range, this study proposes the following enhancements:

• Frequency-dependence elimination for the 2D acoustic state boundary integral equation using Taylor series expansion.

This article has the following structure. Section 2 details the BEM formulation for acoustic state analysis. Section 3 elaborates on the BEM formulation with Taylor series expansion for acoustic state analysis. MOR built on the adaptive SOAR method is employed in Section 4 to speed up the BEM computation for 2D broadband acoustic situations. The treatment of singular integrals in kernel function border integrals is described in Section 5. Section 6 includes multiple numerical instances that verify the effectiveness of the proposed algorithm. Finally, conclusions and further discussions are drawn in Section 7.

The Helmholtz half-space problem can be represented by the following BIE and normal derivative boundary integral equation (HBIE). C x p x + ∫ S F x , y p y d S y = ∫ S G x , y q y d S y + p inc x (1) and C x q x + ∫ S H x , y p y d S y = ∫ S K x , y q y d S y + ∂ p inc x ∂ n x (2) where y signifies the field point, x denotes the source point, and q is the normal derivative of the sound pressure p : q ( x ) = ∂ p ( x ) / ∂ n ( x ) . When x is located on a border that is smooth S , C ( x ) = 1 / 2 . Acoustic pressure incident at position x is given by p inc . The function of Green G ( x , y ) and its derivatives in Equation 1 and Equation 2 are presented as follows, which can be expressed in Equation 3. G x , y = i 4 H 0 1 k r F x , y = ∂ G x , y ∂ n y = − i k 4 H 1 1 k r ∂ r ∂ n y K x , y = ∂ G x , y ∂ n x = − i k 4 H 1 1 k r ∂ r ∂ n x H x , y = ∂ 2 G x , y ∂ n x ∂ n y = i k 4 r H 1 1 k r n j x n j y − i k 2 4 H 2 1 k r ∂ r ∂ n x ∂ r ∂ n y (3) where the n th order first kind Hankel function is indicated by H n ( 1 ) , k denotes by the wave number, i = − 1 , n j is the Cartesian component of n ( x ) or n ( y ) and the distance between the field and source locations is represented by the formula r = x − y .

For exterior acoustic problems, using either Equation 1 or Equation 2 alone can lead to non-uniqueness of the solution at certain imaginary frequencies. According to the Burton-Miller idea, a linear combination of Equation 1 and Equation 2 can effectively resolve this issue. The Burton-Miller formulation is expressed as follows [46]. C x p x + α q x + ∫ S F x , y + α H x , y p y d S y = ∫ S G x , y + α K x , y q y d S y + p i n c x + α ∂ p i n c x ∂ n x (4) in which α represents the coupling parameter: defined as α = i / k where k > 1 and α = i in other cases.

By discretizing the structural boundary into several elements using constant elements and introducing the coefficient matrix, Equation 4 can be reformulated as follows H p − G q = p f inc (5) where H and G ∈ C N × N ( N indicates the quantity of degrees of freedom.) are the coefficient matrices. They are asymmetric, fully populated, and frequency-related. The column vectors p and q , respectively, represent the sound pressure and the acoustic flux at the collocation locations. p f inc is the vector of the incident wave.

To determine the sound pressure values at the boundary surface nodes, Equation 5 needs to be solved. Subsequently, the sound pressure can be computed at any point within the acoustic domain by using Equation 4 with α = 0 and C ( x ) = 1 .

Its sound pressure, p f , can be written as follows if the computation takes into account the external acoustic field. p f = − H f p − G f q + p f inc (6) where the matrices H f and G f as well as the vector p f are similar to those in Equation 5, except that the source point x is outside the structure domain. And H f and G f ∈ C N ̃ × N ( N ̃ indicates the quantity of degrees of freedom.)

In broadband acoustic analysis using BEM, the frequency-dependent coefficient matrices result in time-consuming repetitive computations and repeated solutions of the system equations, posing challenges for practical engineering applications. This frequency dependence arises because the underlying solution is inherently frequency-related. To address this issue, the frequency dependence of the coefficient matrix is eliminated by applying the Taylor series theorem.

The Taylor series expansion of H n ( 1 ) ( k r ) [47] is given by H n 1 k r = ∑ m = 0 ∞ k r − k 0 r m m ! H n 1 k r k r = k 0 r m (7) where k 0 is a fixed frequency expansion point. As seen in Figure 1, the fixed expansion point in this paper is the midway point of the frequency band.

Considering an incident wave traveling along the ℓ -axis, the Taylor series expansion of the term related to p inc ( x ) in Equation 4 can be expressed as follows p inc x + α ∂ p inc x ∂ n x = e i k x ℓ 1 + α i k ∂ x ℓ ∂ n x (8) where x ℓ denotes the ℓ -axis coordinate of the source point, with ℓ being either x or y .

The terms of incident wave in Equation 8 are Taylor series expanded into e i k r = e i k 0 r ∑ m = 0 ∞ i r m k − k 0 m m ! (9)

The expansion expression for the integral of the kernel functions is obtained by combining Equation 7 and Equation 4. ∫ S G x , y q y d S y = ∑ m = 0 ∞ k − k 0 m m ! I g m ∫ S F x , y p y d S y = ∑ m = 0 ∞ k − k 0 m m ! I f m ∫ S α K x , y q y d S y = ∑ m = 0 ∞ k − k 0 m m ! k I k m ∫ S α H x , y p y d S y = ∑ m = 0 ∞ k − k 0 m m ! k I h 1 m + k 2 I h 2 m (10) where I g m = ∫ S i r m 4 H 0 1 k r k r = k 0 r m q y d S y I f m = − ∫ S i r m − 1 4 k r H 1 1 k r k r = k 0 r m ∂ r ∂ n y p y d S y I k m = − ∫ S α i r m 4 H 1 1 k r k r = k 0 r m ∂ r ∂ n x q y d S y I h 1 m = ∫ S α i r m − 1 4 H 1 1 k r k r = k 0 r m n j x n j y p y d S y I h 2 m = − ∫ S α i r m 4 H 2 1 k r k r = k 0 r m ∂ r ∂ n x ∂ r ∂ n y p y d S y (11) where k r H 1 1 k r m = m H 1 1 k r m − 1 + k r H 1 1 k r m (12) By combining Equation 12 and Equation 9, an expression related to the incident wave in Equation 4 is obtained. p ̃ inc = ∑ m ̄ = 0 ∞ k − k 0 m ̄ m ̄ ! P ̃ inc, 1 m ̄ + k P ̃ inc, 2 m ̄ (13) where the respective components can be expressed as Equation 14. P ̃ inc, 1 m ̄ = e i k 0 x ℓ i x ℓ m ̄ P ̃ inc, 2 m ̄ = e i k 0 x ℓ i x ℓ m ̄ α i ∂ x ℓ ∂ n x (14)

A new formulation of Equation 4 is given below by substituting Equation 10 and Equation 13 into Equation 4. C x p x + α C x q x + ∑ m = 0 ∞ k − k 0 m m ! I f m − I g m + k I h 1 m − I k m + k 2 I h 2 m = ∑ m ̄ = 0 ∞ k − k 0 m ̄ m ̄ ! P ̃ inc, 1 m ̄ + k P ̃ inc, 2 m ̄ (15) By discretizing Equation 15 using constant elements, the following matrix expression is derived. ∑ m = 0 ∞ k − k 0 m m ! I p 1 m + k I p 2 m + k 2 I p 3 m p ̃ − I q 1 m + k I q 2 m q ̃ = ∑ m ̄ = 0 ∞ k − k 0 m ̄ m ̄ ! P ̃ inc, 1 m ̄ + k P ̃ inc, 2 m ̄ (16) Due to the linearization of the BEM system, the solutions p ̃ and q ̃ in Equation 16 with truncation term M ̄ can be expressed as Equation 17 and Equation 18. p ̃ = ∑ m ̄ = 0 M ̄ k − k 0 m ̄ m ̄ ! p ̃ 1 m ̄ + ∑ m ̄ = 0 M k k − k 0 m ̄ m ̄ ! p ̃ 2 m ̄ (17) and q ̃ = ∑ m ̄ = 0 M ̄ k − k 0 m ̄ m ̄ ! q ̃ 1 m ̄ + ∑ m ̄ = 0 M k k − k 0 m ̄ m ̄ ! q ̃ 2 m ̄ (18) where the respective solutions of the following system equations are p ̃ 1 m ̄ , p ̃ 2 m ̄ , q ̃ 1 m ̄ and q ̃ 2 m ̄ . ∑ m = 0 ∞ k − k 0 m m ! I p 1 m + k I p 2 m + k 2 I p 3 m p ̃ j m ̄ − I q 1 m + k I q 2 m q ̃ j m ̄ = P ̃ inc,j m ̄ j = 1,2 (19)

The coefficient matrices I p 1 m , I p 2 m , I p 3 m , I q 1 m and I q 2 m ∈ C N × N in Equation 19 are frequency unrelated. Thus, it only needs to be solved once for wide-frequency acoustic problems, eliminating the need for repetitive calculations of coefficient matrices in BEM systems. However, there is another disadvantage of this method, which is that it is still very difficult to solve the equations directly using GMRES for large problems with multiple frequencies because the coefficient matrices are full-rank and asymmetric and the truncation terms require high storage capacity ( O ( 5 ( M + 1 ) N 2 ) ) . It is evident from observation that Equation 19 is the second-order system equation concerning frequency. An efficient SOAR method is presented in [48] to accelerate the broadband solution of Equation 19.

In this section, an effective SOAR method is proposed to expedite the solution of the coefficient matrix by reducing the dimensionality of the 2D system. This projection technique is predicated on the Krylov subspace of second order. By using the method, a system with an equivalent second-order structure but with a diminished state space dimension is created. In Equation 19, the frequency-independent coefficients I p 1 m , I p 2 m , I p 3 m , I q 1 m and I q 2 m are utilized to construct the frequency-unrelated orthogonal basis iteratively, implementing the SOAR algorithm.

In this work, the scattering of an incident wave by a rigid structural surface is considered. The vectors q in Equation 5 vanish because the particle velocity on the surface of the structure is zero. Thus, Equation 19 can be rewritten as ∑ m = 0 ∞ k − k 0 m m ! I p 1 m + k I p 2 m + k 2 I p 3 m p ̃ j m ̄ = P ̃ inc , j m ̄ j = 1,2 (20)

When m = 0 , the coefficients in Equation 19 are utilized to create frequency-independent orthogonal bases, then Equation 19 is re-expressed as I p 1 0 + k I p 2 0 + k 2 I p 3 0 p ̃ j m ̄ = P ̃ inc , j m ̄ j = 1,2 (21)

It should be noted that Equation 21 is not an approximation of the original system equation; rather, it is used for the construction of the orthogonal basis. Equation 21, approximated around a chosen expansion point k 0 , can be expressed as Equation 22. I p 5 0 + k − k 0 I p 4 0 + k − k 0 2 I p 3 0 p ̃ j m ̄ = P ̃ inc , j m ̄ j = 1,2 (22) where I p 5 0 = k 0 2 I p 3 0 + k 0 I p 2 0 + I p 1 0 and I p 4 0 = 2 k 0 I p 3 0 + I p 2 0 .

Following the SOAR method procedure outlined in Ref. [48], a sequence of frequency-independent orthogonal bases Q j m ̄ was constructed in the second-order Krylov subspace G n ( A , B ; r 0 ) using the coefficients from Equation 19, where n ≪ N as shown in Equation 23. span Q j m ̄ = G n A , B ; r 0 = span r 0 , r 1 , r 2 , … , r n − 1 (23) where the respective quantities are shown in Equation 24. A = − I p 5 0 − 1 I p 4 0 B = − I p 5 0 − 1 I p 3 0 r 0 = I p 5 0 − 1 P ̃ inc , j m ̄ r 1 = A r 0 r ℓ = A r ℓ − 1 + B r ℓ − 2 for ℓ ≥ 2 (24)

To define a reduced system equation of the initial Equation 20, the projection subspace is constructed by spanning a sequence of nonzero columns of Q j m ̄ . Equation 20 is reformulated as Equation 25. p ̃ j m ̄ = ∑ m = 0 M k − k 0 m m ! I p 1 m + k I p 2 m + k 2 I p 3 m − 1 P ̃ inc , j m ̄ j = 1,2 (25) Subsequently, the simplified system equation at the expansion point k 0 can be stated as Equation 26. b j m ̄ = ∑ m = 0 M k − k 0 m m ! I p 1 m , n + k I p 2 m , n + k 2 I p 3 m , n − 1 Q j m ̄ T P ̃ inc , j m ̄ j = 1,2 (26) where the coefficient matrix is shown in Equation 27. I p 1 m , n = Q j m ̄ T I p 1 m Q j m ̄ I p 2 m , n = Q j m ̄ T I p 2 m Q j m ̄ I p 3 m , n = Q j m ̄ T I p 3 m Q j m ̄ (27)

The relation between the solution of FOM and ROM is expressed as p ̃ j m ̄ = Q j m ̄ b j m ̄ (28)

p ̃ j m ̄ is the n th order Padé-type approximation of P ̃ inc , j m ̄ about the fixed expansion point k 0 , as shown in Equation 29. p ̃ j m ̄ = P ̃ inc,j m ̄ + O k − k 0 n (29)

In the ROM, I p 1 m , n , I p 2 m , n and I p 3 m , n are n × n matrices where n ≪ N . This significantly reduced storage requirements and improved computational efficiency. At the collocation points on the structural surface, the sound pressure can be obtained using Equation 28. Afterward, at any point within the acoustic domain, Equation 6 can be used to get the sound pressure value.

The boundary integral of the kernel function in Equation 11 can be expressed as the sum of singular and non-singular terms. ∫ S f y d S y = ∫ S \ S x f y d S y ⏟ non−singular + ∫ S x f y d S y ⏟ singular (30) where S x is the element that contains the source point x . S \ S x denotes the boundary S that does not contain S x . Since the integral does not exhibit singularity on the boundary S \ S x , Gauss quadrature can be employed for its solution. However, the integral over the boundary S x is singular and must be processed to do a numerical calculation.

The steps for handling singular integrals are as follows:

1. Suppose there is a singular function f ^ ( x ) . First, perform a singular order analysis on the integrand to obtain a simple function f ^ 0 ( x ) that shares the same order of singularity.

In this study, singular integrals are handled using the Cauchy principal value and the Hadamard finite part integral methods. With reference to Figure 2, let S ε represent a semicircle of radius ε , and Γ ε represent S x \ S ε . Equation 30 can be used to rewrite the singular integral term in Equation 31. ∫ S x f y d S y = lim ε → 0 ∫ Γ ε f y − D y d S y ⏟ n o n s i n g u l a r + lim ε → 0 ∫ S ε f y d S y ⏟ s i n g u l a r + lim ε → 0 ∫ Γ ε D y d S y ⏟ s i n g u l a r (31)

When y ∈ S x , then ∂ r ∂ n ( x ) = ∂ r ∂ n ( y ) = 0 , hence I f m , I k m and I h 2 m are all zero. Therefore, the singular integrals in Equation 11 are present only in I g m and I h 1 m .

To facilitate the derivation of singular integrals, the singular terms of I g m and I h 1 m are denoted as f g m and f h 1 m respectively. f g m = r m H 0 1 k r k r = k 0 r m f h 1 m = r m − 1 H 1 1 k r k r = k 0 r m (32)

Based on Equation 31 and Equation 32, the integrals I g m and I h 1 m in S x can be restated as ∫ S x f g m d S y = lim ε → 0 ∫ Γ ε f g m − D g m d S y ⏟ n o n s i n g l a r + lim ε → 0 ∫ S ε f g m d S y ⏟ g 0 m + lim ε → 0 ∫ Γ ε D g m d S y ⏟ d 0 m ∫ S x f h 1 m d S y = lim ε → 0 ∫ Γ ε f h 1 m − D h 1 m d S y ⏟ n o n s i n g l a r + lim ε → 0 ∫ S ε f h 1 m d S y ⏟ g 1 m + lim ε → 0 ∫ Γ ε D h 1 m d S y ⏟ d 1 m (33) where n j ( x ) n j ( y ) = 1 and d r = d S ( y ) in S x . The non-singular terms in Equation 33 are computed using Gaussian quadrature. The singular part of H 0 ( 1 ) ( k r ) is 2 i ln ( k r ) / π and that of H 1 ( 1 ) ( k r ) is 2 i / ( − π k r ) + i ( k r ) ln ( k r ) / π Thus, their m-th derivatives are represented as Equations 34, 35. D g m = 2 i π ln k 0 r , m = 0 O r 0 , m ≠ 0 (34) and D h 1 m = − 2 i π k 0 − 1 r − 2 + i k 0 π ln k 0 r , m = 0 2 i π k 0 − 2 r − 2 + i π ln k 0 r , m = 1 − 2 i π − 1 m k 0 − m + 1 r − 2 m ! , m ≥ 2 (35)

By substituting Equation 34 into the first equation of Equation 33, the two singular terms g 0 0 and d 0 0 in Equation 33 can be expressed as Equation 36. g 0 0 = 0 d 0 0 = 2 i L π ln k 0 L / 2 − 1 (36) where L is the length of the element.

Similarly, by substituting Equation 34 into the second equation of Equation 33, the expressions for g 1 m and d 1 m are obtained. g 1 m = lim ε → 0 4 i π − 1 m k 0 − m + 1 m ! 1 ε (37) and d 1 m = 8 i π k 0 L + i k 0 L π ln k 0 L / 2 − 1 − lim ε → 0 4 i π k 0 ε m =  0 − 8 i π k 0 2 L + i L π ln k 0 L / 2 − 1 + lim ε → 0 4 i π k 0 2 ε m =   1 8 i π L − 1 m k 0 − m + 1 m ! − lim ε → 0 4 i π − 1 m k 0 − m + 1 m ! 1 ε m ≥ 2 (38)

By combining Equation 37 and Equation 38, the expression for the singular integral is derived as Equation 39. g 1 m + d 1 m = 8 i π k 0 L + i k 0 L π ln k 0 L / 2 − 1 m =  0 − 8 i π k 0 2 L + i L π ln k 0 L / 2 − 1 m =   1 8 i π L − 1 m k 0 − m + 1 m ! m ≥ 2 (39)

Three computational examples are presented in this section to evaluate the performance of the proposed algorithm. The research employed the Fortran 90 programming language to conduct numerical simulations. 64 GB of RAM and an Intel (R) Core (TM) i9-10900H Central Processing Unit (CPU) were installed on a desktop computer to perform calculations.

In the simulations, some common parameters are as follows: an incident wave with an amplitude of p 0 = 1 propagating in the positive x -axis direction, expressed as p i n c = p 0 e i k r ⁡ cos ⁡ θ . The medium for the acoustic wave is air, with a density of 1.21 k g / m 3 and a speed of 343 m / s .