© The Author(s), Published by EDP Sciences, 2024

1 Introduction

The concept of acoustic black hole (ABH) was introduced in 2002 by Mironov and Pislyakov in their seminal paper [1]. The authors proposed a retarding structure consisting of a piece of tapered waveguide with a locally compliant wall, which serves as an anechoic termination for a fluid-filled duct. According to their analysis, in which the waveguide radius was assumed to decrease linearly towards the end of the structure, the tapering together with the wall compliance cause that the impinging acoustic waves slow down, theoretically to the zero velocity at the ABH end, and thus, they cannot reach the end of the structure and reflect back in a finite time. Even if the presented theory has obvious limits, for example, the reduction of the wave speed to the zero value requires the waveguide with null radius at its end, the structure can still significantly reduce the wave reflection. As the impinging wave slows down, its energy increases, and it is partly absorbed thanks to weak, but omnipresent thermoviscous processes.

An implementation of this ABH, proposed in paper [1], is depicted in Figure 1. It consists of a set of N rigid circular ribs, they are attached to a waveguide with inner radius R, and have central holes, the radii of which vary smoothly along the structure whose total length is L. The end of the structure is terminated with a rigid wall. The aim of the structure is to absorb the acoustic plane wave $p_i^{\text{'}}$ impinging from the right side from an infinite waveguide, as it is depicted in Figure 1. The tapered compliant wall of the ABH is formed by the cavities between the ribs; the authors derived a 1D wave equation describing the wave propagation within the ABH.

Figure 1 Model of an acoustic black hole.

In subsequent theoretical works [2, 3], the behavior of ABHs with quadratic radius profile has been modeled employing the transfer matrix method (TMM), the application of fractional calculus has been also explored [4]. Acoustic black holes have also been studied experimentally [5, 6]. The authors of these works constructed ABHs with linear and quadratic radius profiles and demonstrated a reduction of the reflection coefficient when the ABH was used as an absorbing duct termination. In the paper [7], two types of retarding structures for which, in an ideal case, wave speed decreases to zero have been proposed. Some early experimental results have also been provided. Mi et al. [8] studied the performance of an ABH with a quadratic profile and closed or open end. They performed experimental measurements on a rather fine (compared to the previous works [5–7]) 3D-printed structure. Various modifications and applications of ABHs have also been explored, such as open mufflers inspired by ABHs [9], mufflers based on ABHs combined with micro-perforated plates [10, 11], a duct with embedded periodic ABHs [12], or thin ABH-inspired metamaterials [13].

In 2022 (that is more than 20 years after the original paper [1] by Mironov and Pislyakov), three papers [14–16] have been published, where the following mechanisms important for the understanding of ABHs were pointed out. First, acoustic energy dissipation in the acoustic boundary layer adjacent to the solid-fluid interfaces, especially in the cavities between the ribs, must be properly accounted for. Second, “local” resonance taking place in the cavities between the ribs plays an important role in the mechanism of the acoustic energy absorption (the mechanism is similar to the rainbow trapping [17]), significantly improving the ABHs performance especially at higher frequencies. Moreover, the role of acoustic energy radiation from the walls of an ABH was also discussed [14]. In paper [16], the authors also studied a wider range of shapes of ABHs, and the role of vibration of the ribs on the ABHs’ absorbing properties. Note, that properties of ABHs with rectangular cross-section have also been studied [18–21].

In most of the previously referred works, ABHs with simple profiles (mainly linear and quadratic) have been studied. As far as we are aware of, there has been only a few papers published attempting to optimize an ABH profile to enhance its performance.

Bravo and Maury [22] studied a fully-opened ABH-type silencer made up of side-branch annular cylindrical cavities (resonators) with axially increasing depths. As the inner radius of this retarding structure is constant, it only half-slows down the incident sound wave. Employing an optimization, they determined the structure profile, wall porosity, and the density of the resonators to maximize the acoustic energy absorption.

Serra et al. [23] used an optimization algorithm to determine a suitable profile of an ABH. Following work [1], they employed a simple ad-hoc mathematical model of losses in the ABH and also neglected the effect of local resonances, which may strongly influence the ABH behavior at higher frequencies. On the other hand, they studied the effect of added absorbing material and its distribution inside the slits.

In a very recent paper, Mousavi et al. [24] proposed a new type of ABH, which impedance-matches the external waveguide to a small perfect absorber in a wide frequency rage, which prevents the reflection of the incident wave from the structure. To achieve this goal, the authors employed topology optimization to design the ABH internal structure, which very much differs from the ribbed geometry proposed by Mironov and Pislyakov [1].

Within this work, we employ evolutionary optimization to find the optimum ABH shape maximizing its absorption properties without using any other means, such as added absorbing materials or perforated plates. In our model, we only take into account, and take advantage of, the intrinsic losses in thermoviscous acoustic boundary layer. We also identify other parameters influencing the ABH performance as a broadband acoustic energy absorber.

The rest of this paper is organized as follows. In Section 2 we discuss three mathematical models for the calculation of absorbing properties of ABHs. In Section 3 we introduce an objective function for the ABH performance optimization, its shape parametrization, and the numerical procedure. In Section 4 we validate the computationally-efficient mathematical model used for the optimization and conduct a parametric study on the influence of the ABH shape and other parameters on its absorbing properties. Some concluding remarks are drawn in Section 5.

2 Mathematical model

In this section, we describe two mathematical models. The first is a detailed and computationally intensive 2D model that resolves the inner structure of the ABHs. The second is a simplified, computationally efficient 1D model used for optimization purposes. The 2D model is used to validate the numerical results obtained from the 1D model.

Throughout this work, all the acoustic-field quantities are assumed to be harmonic in time, their complex amplitudes denoted with hat are introduced as $q(\mathit{r}, t) = \mathfrak{R}\left[\widehat{q}\right(r) \exp (\mathrm{i}\omega t\left)\right]$, r is the position vector, t is the time, ω is the angular frequency, and $\mathrm{i} =\sqrt{-1}$ is the imaginary unit.

2.1 2D model

The studied structures, as shown in Figure 1, possess an axial symmetry, the acoustic field is also assumed to have this symmetry, and thus a 2D mathematical model formulated in r–z axi-symmetric cylindrical coordinates is adopted here. The model outlined within this section is the same as the one developed in our previous work [15] so that only a brief description is given here.

Detail of the 2D model of an ABH is schematically depicted in Figure 2. The individual circular ribs are modeled as rigid, they are placed with the spatial period H = L/N. The rigid wall terminating the structure is positioned at z = 0. The width of the slits between the ribs is set as h = ξH, where ξ ∈ (0, 1). From here, it follows that the width of the ribs h_rib = H – h = (1 − ξ)H. The centers of the individual ribs are positioned at

$$z_j=[j+(\xi -1)/2]H,$$(1)

where j = 1, 2, ..., N, and the inner radii r_j of the ribs, see Figure 2, are given as r_j = r(z_j), where r(z) is a function defining the ABH profile. The outer radius of the ribs R is the same as the waveguide inner radius, the waveguide wall is also assumed to be rigid.

Figure 2 Detail of an acoustic black hole.

As it has been demonstrated in recent works [14–16], thermoviscous losses in the boundary layer adjacent to the fluid-solid boundaries, especially in the slits between the ribs play an important role in the mechanism of the acoustic-energy absorption in ABHs.

In order to capture properly these losses, the acoustic field is modeled by the linearized Navier-Stokes equations (LNS), which in the frequency domain, see e.g. [25, 26], have the form

$$\mu \left[{\nabla }^2\widehat{\mathit{v}}+\frac{1}{3}\nabla (\nabla \cdot \widehat{\mathit{v}})\right]-\nabla \widehat{p}=\mathrm{i}\omega {\rho }_0\widehat{\mathit{v}},$$(2a)

$$-{\rho }_0\nabla \cdot \widehat{\mathit{v}}=\mathrm{i}\omega \widehat{\rho },$$(2b)

$$\kappa {\nabla }^2\widehat{T}=\mathrm{i}\omega {\rho }_0{c}_p\widehat{T}-\mathrm{i}\omega \widehat{p},$$(2c)

where $\widehat{\mathit{v}, } \widehat{T}, \widehat{\rho }$ are the acoustic velocity vector, acoustic temperature, and acoustic density complex amplitudes, respectively, ρ₀ is the equilibrium fluid density, μ is the coefficient of shear viscosity, κ is the coefficient of thermal conduction, and c_p is the specific heat at constant pressure. Equations (2) are supplemented with the linearized state equation for an ideal gas $\widehat{p}/p_0= \widehat{\rho }/{\rho }_0+ \widehat{T}/T_0$, where T₀, p₀ are the equilibrium values of the gas temperature and pressure. Equations (2) are solved with the isothermal ($\tilde{T} = 0$) and no-slip ($\tilde{\mathit{v}} = \mathbf{0}$) boundary conditions on the rigid walls.

Equations (2) were solved numerically employing commercial FEM software COMSOL Multiphysics (Acoustic Module, Thermoviscous Acoustics Interface, frequency domain, 2D axi-symmetric geometry). Further details are given in our previous work [15].

2.2 1D models

The 2D model accurately captures the thermoviscous losses, as it resolves the boundary layer, as well as other effects, such as evanescent coupling of the adjacent slits, or the “end effects” for the slits and short waveguide segments formed by the central holes in the ribs. On the other hand, this model is rather computationally demanding and thus it is not suitable for the purpose of ABH profile optimization. That is why a simplified, computationally efficient 1D mathematical model of an ABH is needed. Within the 1D model, the acoustic field quantities are assumed to be only the functions of the spatial coordinate z, measuring the position along the symmetry axis. Therefore, the model is only valid up to the cut-on frequency

$$f_{\mathrm{cut}-\mathrm{on}} =\frac{1.84}{2\pi } \frac{c_0}{R},$$(3)

where c₀ is the sound speed.

2.2.1 Wave equation

A one-dimensional wave equation describing the acoustic wave propagation in an ABH was first introduced in the work [1] by Mironov and Pislyakov. In their paper, the ABH was modeled as a piece of tapered waveguide with a locally compliant wall characterized by its complex acoustic admittance ${\tilde{Y}}_w$, see Figure 3.

Figure 3 Model of an acoustic black hole.

The wave equation adopted within this work is similar to the one proposed by Umnova et al. [16]. It is represented by a generalized Webster equation for a waveguide with a compliant wall and filled with an inhomogeneous fluid. The inhomogeneous equivalent fluid is introduced to capture the thermoviscous losses in the main duct of the ABH, caused by isothermal and no-slip boundary conditions at the ends of the ribs, forming a boundary layer. The effective density ${\tilde{\rho }}_{0\mathrm{g}}$ and sound speed ${\tilde{c}}_{0\mathrm{g}}$ in this equivalent fluid are introduced employing the Stinson’s model [27] for a circular pore with radius r. As in this case this radius is represented by the ABH profile r(z), the equivalent fluid properties are also position-dependent, namely, ${\tilde{\rho }}_{0\mathrm{g}}={\tilde{\rho }}_{0\mathrm{g}}\left(z\right),$ ${\tilde{c}}_{0\mathrm{g}}={\tilde{c}}_{0\mathrm{g}}\left(z\right)$. The wave equation reads

$$\frac{{\tilde{\rho }}_{0\mathrm{g}}}{r^2}\frac{\mathrm{d}}{\mathrm{d}z}\left(\frac{r^2}{{\tilde{\rho }}_{0\mathrm{g}}}\frac{\mathrm{d}\widehat{p}}{\mathrm{d}z}\right)+{\tilde{k}}^2\widehat{p}=0,$$(4)

where the complex wavenumber $\tilde{k} = \tilde{k}\left(z\right)$ reads

$$\tilde{k} =\sqrt{{\tilde{k}}_{\mathrm{g}}^2-2\mathrm{i}{\tilde{k}}_{\mathrm{g}}{\tilde{Z}}_{0\mathrm{g}}{\tilde{Y}}_{\mathrm{w}}/r},$$(5)

where ${\tilde{k}}_{\mathrm{g}}=\omega /{\tilde{c}}_{0\mathrm{g}}$, and ${\tilde{Z}}_{0\mathrm{g}}={\tilde{\rho }}_{0\mathrm{g}}{\tilde{c}}_{0\mathrm{g}}$.

It has been recently demonstrated [14–16] that acoustic resonance in slits between the ribs forming the ABH structure plays an important role in the mechanism of the acoustic energy absorption. That is why more accurate formula for the slits’ input acoustic admittance than the one used by Mironov and Pislyakov [1] was employed here. The wall acoustic admittance ${\tilde{Y}}_{\mathrm{w}}\left(z\right)$ is set as ${\tilde{Y}}_{\mathrm{w}}\left(z\right) = \xi {\tilde{Y}}_{\mathrm{s}}\left(z\right)$, where the slits’ input acoustic admittance, see e.g. [28], reads

$${\tilde{Y}}_{\mathrm{s}}\left(z\right)= -\mathrm{i}{\tilde{Y}}_{0\mathrm{s}}\frac{{\mathrm{J}}_1\left({\tilde{k}}_{\mathrm{s}}r\right){\mathrm{Y}}_1\left({\tilde{k}}_{\mathrm{s}}R\right)– {\mathrm{Y}}_1\left({\tilde{k}}_{\mathrm{s}}r\right){\mathrm{J}}_1\left({\tilde{k}}_{\mathrm{s}}R\right)}{{\mathrm{J}}_0\left({\tilde{k}}_{\mathrm{s}}r\right){\mathrm{Y}}_1\left({\tilde{k}}_{\mathrm{s}}R\right)– {\mathrm{Y}}_0\left({\tilde{k}}_{\mathrm{s}}r\right){\mathrm{J}}_1\left({\tilde{k}}_{\mathrm{s}}R\right)},$$(6)

where ${\tilde{k}}_{\mathrm{s}}=\omega /{\tilde{c}}_{0\mathrm{s}}$, and ${\tilde{Y}}_{0\mathrm{s}}=1/{\tilde{\rho }}_{0\mathrm{s}}{\tilde{c}}_{0\mathrm{s}}$. The quantities ${\tilde{\rho }}_{0\mathrm{s}}$ and ${\tilde{c}}_{0\mathrm{s}}$ are the density and sound speed in an equivalent fluid modeling the thermoviscous losses in boundary layer adjacent to the boundaries of the slits. For this purpose, Stinson’s model [27] for a slit with width h was used. The coefficient ξ in the wall admittance formula represents the fact that the filling fraction of the slits in the “waveguide wall” in the 1D model is ξ and the ribs are considered rigid.

2.2.2 Riccati equation

If the acoustic field properties in an ABH are not needed to be known, the solution of the two-point boundary value problem for the wave equation (4) can be avoided. Absorption properties of an ABH can be inferred from the solution of the Riccati equation [16], which is a first-order differential equation and thus its numerical solution is easier. The Riccati equation for acoustic admittance inside an ABH, $\tilde{Y}=\widehat{v}/\widehat{p}$, where $\widehat{v}$ is the complex amplitude of acoustic velocity, can be written as

$$\frac{\mathrm{d}\tilde{Y}}{\mathrm{d}z}=\mathrm{i}{\tilde{k}}_{\mathrm{g}}{\tilde{Z}}_{0\mathrm{g}}{\tilde{Y}}^2-\frac{\mathrm{i}{\tilde{k}}_{\mathrm{g}}}{{\tilde{Z}}_{0\mathrm{g}}}-\frac{2}{r}\frac{\mathrm{d}r}{\mathrm{d}z}\tilde{Y}-\frac{2{\tilde{Y}}_{\mathrm{w}}}{r}.$$(7)

The Riccati equation (7) is numerically solved (integrated) with the initial condition $\tilde{Y}(z=0)=0$, representing the rigid wall terminating the ABH, in the interval z ∈ 〈0, L〉. This way, the ABH input acoustic admittance ${\tilde{Y}}_{\mathrm{ABH}}=\tilde{Y}\left(L\right)$. As it is assumed that the characteristic admittance of the infinite waveguide with radius R in front of the ABH is Y₀ = 1/ρ₀c₀, the reflection coefficient from the ABH is

$$\mathcal{R}=\frac{R^2+{\rho }_0{c}_0{r}_{\mathrm{ABH}}^2{\tilde{Y}}_{\mathrm{ABH}}}{R^2-{\rho }_0{c}_0{r}_{\mathrm{ABH}}^2{\tilde{Y}}_{\mathrm{ABH}}}{\mathrm{e}}^{ 2\mathrm{i}\omega L/c_0},$$(8)

where r_ABH = r(L) represents the input radius of the ABH (at z = L), which may differ from (be smaller than) the waveguide radius R.

3 ABH shape optimization

3.1 ABH performance evaluation

In order to be able to optimize the ABH performance, it is necessary to choose a suitable quantity describing the ABH absorbing performance. A plausible choice here is the reflection coefficient modulus which can be calculated employing equation (8). As it is a frequency-dependent quantity, its mean value in a given frequency range was used in this work as the ABH performance quantifier

$$\langle \left|\mathcal{R}\right|\rangle =\frac{1}{f_{\max }-f_{\min }}{\int }_{f_{\min }}^{f_{\max }}\mathrm{}\left|\mathcal{R}\right(f\left)\right| \mathrm{d}f,$$(9)

where the frequency limits can be set as f_min = 0 Hz, f_max = f_cut-on (unless explicitly stated otherwise, $\langle \left|\mathcal{R}\right|\rangle$ will always be calculated in this frequency range). In theory, the value $\langle \left|\mathcal{R}\right|\rangle$, depending on the ABH parameters, can vary from zero (perfect absorption) to unity (perfect reflection).

As all the calculations are performed numerically at discrete frequencies, the integral (9) is evaluated numerically employing the trapezoidal rule.

Note that there are other approaches how to quantify the performance of a wide-band absorber, causal-based optimization [29] can serve as an example.

3.2 ABH shape parametrization

It has been demonstrated in the previous works [1, 2, 15, 16] that the absorbing performance of an ABH is strongly dependent on its profile shape and the minimum inner radius. Within this work, the ABH shape is introduced as

$$r\left(z\right)=\left(R-R_{\min }\right)F\left(z/L\right)+R_{\min },$$(10)

where the ABH shape function 0 ≤ F(Z) ≤ 1 for Z ∈ 〈0, 1〉.

Here, two types of shape functions are employed. First, it is a single-parameter power function

$$F\left(Z\right)=Z^{\alpha },$$(11)

where the values of α = 1, 2 correspond to the previously studied cases [1, 2] of linear and quadratic profiles.

Second, in order to be able to define the ABH profile more freely, the shape function was defined employing the cubic-spline interpolation of a set of M control points. The cubic splines are piecewise-third-order-polynomial smooth functions, whose shape is defined using these control points and the derivatives at the end-points. The higher number of the control points results in the greater variability of the shape function, but, this leads to a more-dimensional parameter space in which it is harder to find the parameters optimizing the ABH performance.

In all the studied cases, the left-most control point has the coordinates [Z₁, F₁] = [0, 0], the right-most one then [Z_M, F_M] = [1, F_M], where the value of F_M is either set to one, or its value is determined by the optimization algorithm. Also, derivatives at the end-points are used as control parameters, namely,

$$A={\frac{\mathrm{d}F}{\mathrm{d}Z}|}_{Z=0}, \hspace{1em}\mathrm{and}\hspace{1em}B={\frac{\mathrm{d}F}{\mathrm{d}Z}|}_{Z=1}.$$

For example, the state vectors X defining the shape function in a 2–7 dimensional parameter space and corresponding control parameters C for the spline functions are given as

$$\begin{array}{lll}{\mathit{X}}_2=[A,B]& \to & {\mathit{C}}_2=\left[A,\left[\mathrm{0,0}\right],\left[\mathrm{1,1}\right],B\right],\\ {\mathit{X}}_3=[A,F_2,B]& \to & {\mathit{C}}_3=\left[A,\left[\mathrm{0,0}\right],[1,F_2],B\right],\\ {\mathit{X}}_4=[A,Z_2,F_2,B]& \to & {\mathit{C}}_4=\left[A,\left[\mathrm{0,0}\right],[Z_2,F_2],\left[\mathrm{1,1}\right],B\right],\\ {\mathit{X}}_5=[A,Z_2,F_2,F_3,B]& \to & {\mathit{C}}_5=\left[A,\left[\mathrm{0,0}\right],[Z_2,F_2],[1,F_3],B\right],\\ {\mathit{X}}_6=[A,Z_2,F_2,Z_3,F_3,B]& \to & {\mathit{C}}_6=\left[A,\left[\mathrm{0,0}\right],[Z_2,F_2],[Z_3,F_3],\left[\mathrm{1,1}\right],B\right],\\ {\mathit{X}}_7=[A,Z_2,F_2,Z_3,F_3,F_4,B]& \to & {\mathit{C}}_7=\left[A,\left[\mathrm{0,0}\right],[Z_2,F_2],[Z_3,F_3],[1,F_4],B\right].\\ & & \end{array}$$(12)

An example of a shape function defined by three control points and five parameters is shown in Figure 4. The values of the control parameters X are limited by the condition 0 ≤ F(Z; X) ≤ 1 for Z ∈ 〈0, 1〉.

Figure 4 ABH shape function defined using cubic splines with 3 control points and 5 parameters.

It should be noted that although the ABH shape described by a 7-dimensional parameter vector is quite versatile, especially in comparison to the power-law shape, cubic splines are limited to creating smooth shapes that are not overly complex.

3.3 Optimization algorithm

As it has been described in Section 3.1, the ABH performance is quantified by means of the mean value of the reflection coefficient modulus $\langle \left|\mathcal{R}\right|\rangle$ introduced using equation (9). Within this work, unless explicitly stated otherwise, the optimization is performed as a search for the parameters of a shape function F(Z), which minimize the objective function represented by the mean value of the reflection coefficient modulus $\langle \left|\mathcal{R}\right|\rangle$. Within the optimization process, the other parameters (L, R_min, N, ξ) are considered constant, and their influence is examined separately.

The numerical value of the objective function results from a rather complex numerical calculation. The ABH reflection coefficient at a given frequency is calculated employing the numerical integration of the Riccati equation (7), and then its mean value is determined by numerical evaluation of a definite integral. As a result, the objective function is “noisy” and numerical calculation of its gradient (with respect to the state vector) is rather inaccurate. Moreover, the parameter space is not continuous, as it is necessary to reject the values of X, for which F(Z; X) < 0 and F(Z; X) > 1 for Z ∈ 〈0, 1〉. These facts, together with the one that the objective function is multimodal (there are multiple local minima), proven an evolutionary heuristic approach to be the most efficient one in searching for the global minimum of the objective function. For this purpose, genetic algorithm implemented in the MATLAB’s Global Optimization Toolbox [30] was used.

4 Numerical results

4.1 Parameters of the numerical simulations

In all the numerical simulations, unless explicitly stated otherwise, the length of the ABH was set to L = 25 cm, and the radius of the waveguide R = 3 cm. The fluid filling the waveguide and ABH is assumed to be air at normal conditions with the following material parameters: the sound speed c₀ = 343.2 m s⁻¹, the ambient density ρ₀ = 1.204 kg m⁻³, the adiabatic exponent γ = 1.4, the dynamic viscosity μ = 1.83 × 10⁻⁵ Pa s, the specific heat capacity at constant pressure c_p = 1004 J kg⁻¹ K⁻¹, and the thermal conductivity κ = 2.59 W m⁻¹ K⁻¹. For these conditions, the cut-on frequency f_cut-on = 3350 Hz.

4.2 Validation of the numerical model

Before the simplified mathematical model based on the Riccati equation is used for the optimization of an ABH performance, it is necessary to validate its predictions by comparison with the 2D mathematical model based on the linearized Navier-Stokes equations and FEM introduced in Section 2.1. In all these comparisons, the minimum ABH inner radius was set to R_min = 0.25 mm.

The frequency characteristics of the reflection coefficient modulus for ABHs with linear profile for N = 250 and 500, and different values of ξ are shown in Figure 5.

Figure 5 Reflection coefficient spectra of linear ABHs calculated employing the 1D model based on the Riccati equation (blue lines), and 2D finite element model (red lines), for different parameters: (a) −N = 500, ξ = 0.5, (b) −N = 250, ξ = 0.25, (c) −N = 250, ξ = 0.5, (d) −N = 250, ξ = 0.75.

Panel (a) of Figure 5 shows the case of N = 500 (H = 0.5 mm), and ξ = 0.5. As the ABH internal geometry is very fine, the 1D model provides essentially the same results as the 2D model. The other panels of Figure 5 show the case of N = 250 (H = 1 mm), and different values of ξ. It can be observed that for N = 250, there is a small discrepancy between the 2D FEM and 1D Riccati results, which increases with the value of ξ.

The difference between the predictions of the 2D model and the 1D model is mainly caused by the fact that the 1D model only accounts for plane-wavemodes and does not capture effects such as evanescent coupling of adjacent slits, end-effects, etc. Moreover, the 1D model employing the Riccati equation entirely ignores the discrete character of the structure. The influences of all these effects decreases with increasing N and decreasing ξ.

Table 1 shows a quantitative comparison of the above-discussed numerical results, namely, the comparison of the mean value of the reflection coefficient modulus $\langle \left|\mathcal{R}\right|\rangle$. It can be seen that in all the cases the difference between the predictions of the 2D and 1D model are less than 1%. It can be also observed that the best (and very similar) performance (the lowest value of $\langle \left|\mathcal{R}\right|\rangle$) is provided by ABHs with (N = 500, ξ = 0.5), and (N = 250, ξ = 0.75).

Taking into account the facts mentioned above, all the results presented below, unless explicitly stated otherwise, are obtained using the 1D model that employs the Riccati equation. The evaluation of the reflection coefficient at a single frequency using the 1D model achieves a typical speed-up factor of 10,000 compared to the 2D model.

4.3 Performance of ABHs with power-law shape function

Within this section, the dependence of the absorption properties of ABHs with shape function given by equation (11) is studied.

Figure 6 shows the dependence of the mean value of the reflection coefficient modulus as a function of the power-law exponent α for different values of R_min; N = 500, ξ = 0.5. It can be seen that the minimum value of $\langle \left|\mathcal{R}\right|\rangle$ and the corresponding exponent α depend on the value of R_min; the lower the R_min, the lower the minimum of $\langle \left|\mathcal{R}\right|\rangle$ and the corresponding value of α. For example, for R_min = 0.05 mm, min$(\langle |\mathcal{R}|\rangle )=0.1130$ (for α = 0.89), for R_min = 0.25 mm, min$(\langle |\mathcal{R}|\rangle )=0.1181$ (for α = 0.96), and for R_min = 0.45 mm, min$(\langle |\mathcal{R}|\rangle )=0.1235$ (for α = 1.03).

Figure 6 Mean value of the reflection coefficient modulus as a function of the power-law exponent for different values of Rmin; N = 500, ξ = 0.5.

The numerical results show that the optimum value of α also depends on the ABH’s internal structure, for example, for N = 250 and ξ = 0.75 and R_min = 0.25 mm, min$(\langle |\mathcal{R}|\rangle )=0.1079$ for α = 1.07, which is very close to the linear ABH profile.

The performance of an ABH also strongly depends on its total length L, which is demonstrated in Figure 7. In all the cases H = 0.5 mm (the number of the ribs N depends on the ABH length), ξ = 0.5, and R_min = 0.25 mm.

Figure 7 Panel (a) – mean value of the reflection coefficient modulus as a function of the power-law exponent α for different values of the ABH length L; Panel (b) – reflection coefficient spectra for different values of L and α; H = 0.5 mm, ξ = 0.5.

It can be observed in Panel (a) of Figure 7 that the minimum value of $\langle \left|\mathcal{R}\right|\rangle$ decreases with increasing ABH’s length, whereas the value of the corresponding α depends on L only weakly. For example, for L = 21 cm, min$(\langle |\mathcal{R}|\rangle )=0.1354$ (for α = 0.95), for L = 29 cm, min$(\langle |\mathcal{R}|\rangle )=0.1049$ (for α = 0.97), It can be observed in Panel (b) of Figure 7 that the reason for lower values of $\langle \left|\mathcal{R}\right|\rangle$ in ABHs with larger length is caused by the fact that the larger L, the lower frequency at which the ABH starts to effectively absorb the acoustic energy. This can be explained by the fact that the low-frequency absorption in these ABHs is connected with global resonances in the main cavity, see the discussion in [16], and the larger the ABH length, the lower these resonant frequencies are. Numerical results show that ABHs with different internal-structure parameters (R_min, H, ξ) behave similarly.

4.4 Performance of ABHs with shape function defined using cubic splines

Parametrization of the ABH shape function employing cubic splines offers greater variability in ABH shape definition, which may result in improved absorbing properties.

Top part of Table 2 shows the values of the $\langle \left|\mathcal{R}\right|\rangle$ for optimized ABHs with cubic-spline shape functions with different number of control parameters (see Eq. (12)). The bottom part, for the sake of comparison, shows the values for linear (non-optimized) ABH, and optimized ABH with power-law shape (α = 0.96). In all the presented cases, N = 500, ξ = 0.5, and R_min = 0.25 mm.

There can be observed a general trend that the higher number of parameters results in lower $\langle \left|\mathcal{R}\right|\rangle$ and the fact that odd number of parameters provides better results. As it can be seen in equation (12), the shape function with odd number of parameters allows F(1) < 1, which introduces a small cross-sectional jump at the ABH opening. This slightly improves the ABH performance at the lower frequencies, probably by increasing the quality factor of the global resonances, and it slightly decreases the ABH performance at the higher frequencies by introducing the impedance mismatch. The overall effect on the value of $\langle \left|\mathcal{R}\right|\rangle$ is positive though.

For example, the mean of the reflection coefficient modulus for the case of 5 control parameters is only 77.3% of the value for the non-optimized linear ABH.

Panel (a) of Figure 8 shows the shape functions for the case of optimization with 5 and 6 parameters, and for the case of a non-optimized linear ABH. It can be seen, that in the case of 5 parameters, F(1) < 1 [F(1) = 0.9881], and for both the optimized cases, A ≈ 0.

Figure 8 Panel (a) – cubic-spline shape functions of optimized ABHs, Panel (b) – reflection coefficient modulus for optimized ABH (cubic-splines, 5 pars.), and linear ABH for comparison; N = 500, ξ = 0.5, Rmin = 0.25 mm.

Panel (b) of Figure 8 shows the reflection coefficient modulus for the optimized ABH with 5 parameters, and the non-optimized linear one for comparison. Figure shows that the shape optimization mainly reduces the reflection coefficient at higher frequencies.

If the mean value of the reflection coefficient modulus is minimized in another frequency range than 0 Hz–f_cut-on (=3350 Hz), the shape function changes correspondingly. Figure 9 shows the shape function and the reflection coefficient modulus for an ABH optimized in a frequency range 0 Hz–1000 Hz (cubic-spline shape function with 5 parameters, N = 500, ξ = 0.5, R_min = 0.25 mm), and for the sake of comparison, the same optimized in the full frequency range (0 Hz–f_cut-on). The low-frequency optimization significantly reduced the reflection coefficient at most frequencies up to 1 kHz, however, the reflection coefficient increased at the higher frequencies. Panel (a) of Figure 9 shows that the reduction of the reflection coefficient at the lower frequencies has been achieved by the reduction of the shape function value at Z = 1: F(1) = 0.876. This cross-sectional jump introduces an impedance mismatch which increases the wave reflection from the ABH at higher frequencies.

Figure 9 Panel (a) – shape functions, and Panel (b) – reflection coefficients of ABHs optimized in full frequency range (0 Hz–fcut-on) and low frequency range (0 Hz–1000 Hz); N = 500, ξ = 0.5, Rmin = 0.25 mm.

As the benefit of employing more control parameters then 5 for the cubic-spline shape function is negligible, 5 control parameters are used within the rest of this work.

Similarly as it has been shown for the optimized ABHs with power-law shape function, the absorbing performance of optimized ABHs with cubic-spline shape function weakly depends on the minimum radius R_min; the smaller R_min, the smaller the value of $\langle \left|\mathcal{R}\right|\rangle$.

For example, for N = 500 and ξ = 0.5, it holds: for R_min = 0.05 mm, $\langle \left|\mathcal{R}\right|\rangle =0.0908$, for R_min = 0.25 mm, $\langle \left|\mathcal{R}\right|\rangle =0.0919$, and for R_min = 0.45 mm, $\langle \left|\mathcal{R}\right|\rangle =0.0943$.

Figure 10 shows in Panel (a) the shape functions, and in Panel (b) the reflection coefficient spectra of optimized ABHs with different lengths. In all the cases, H = 0.5 mm, ξ = 0.5, and R_min = 0.25 mm. It can be seen that the optimal shape function does not change very much with the length L, except for the fact that the value of F(L) decreases with the decreasing length L. Similarly as in the case of ABHs with power-law shape functions, see Figure 7, the absorbing performance of ABHs with cubic-spline shape functions rapidly increases with the increasing length L. For example, increasing the ABH length from 15 cm to 30 cm results in decreasing the value of $\langle \left|\mathcal{R}\right|\rangle$ from 0.1463 to 0.0775.

Figure 10 Panel (a) – shape functions, and Panel (b) – reflection coefficient spectra of optimized ABHs with different lengths; H = 0.5 mm, ξ = 0.5, Rmin = 0.25 mm.

Absorbing performance of optimized ABHs is also influenced by their inner structure. Figure 11 shows the dependence of the mean value of the reflection coefficient modulus of optimized ABHs as a function of the number of the ribs, varying from N = 500 to N = 200. In all the cases, the ABH length is fixed at L = 25 cm, and R_min = 0.25 mm. Three cases are studied: (a) the width of the ribs is kept constant–h_rib = 0.25 mm (the width h of the slits increases with decreasing N), (b) the coefficient ξ is set as ξ = 0.5 (the width of the ribs and slits is the same, both increase with decreasing N), and (c) the width of the slits is kept constant–h = 0.25 mm (the width h_rib of the ribs increases with decreasing N). Note that for N = 500, it holds for all three cases that h = h_rib = 0.25 mm. It can be observed in Figure 11 that in the case of fixed h_rib, the absorption performance increases with the decreasing number of ribs N, which seems to be practical from the manufacturing point of view. In the other two cases, the mean value of the reflection coefficient modulus increases with the decreasing number of the ribs.

Figure 11 Mean value of the reflection coefficient modulus of ABHs optimized for different number of ribs and different internal-structure parameters; L = 25 cm, Rmin = 0.25 mm.

Figure 12 shows the frequency characteristics of the reflection coefficient modulus of ABHs optimized for N = 500, 400, and 300 ribs, respectively; in all the cases, L = 25 cm, h_rib = 0.25 mm, and R_min = 0.25 mm. The solid lines denote the predictions by the 1D model, the dashed lines represent the results obtained employing the 2D FEM model (with the shape function resulting from the optimization based on the 1D model).

Figure 12 Reflection coefficient modulus as a function of frequency for optimized ABHs with different number of ribs N; solid lines–prediction employing the Riccati equation, dashed lines–validation employing 2D FEM model. In all the cases, L = 25 cm, hrib = 0.25 mm, Rmin = 0.25 mm. Inset figure shows the corresponding shape functions.

It can be observed that the reflection coefficient decreases faster at lower frequencies in the case of ABHs with the lower number of ribs, in the middle and high frequency range the behavior is similar. Similarly as it has been shown in Figure 5, the higher the value N the better the correspondence between the prediction by the 1D and 2D model. What is interesting is the fact that even if in the case of N = 300, the predictions of the models are slightly frequency-shifted in the middle frequency range, the corresponding values of the mean value of the reflection coefficient modulus are the same, see the legend in Figure 12. The inset figure in Figure 12 shows the optimum shapes for the individual cases, it can be seen that they slightly differ, namely, with decreasing N, the length of the most contracted part of the ABH increases.

It is also interesting to compare the performance of the optimized ABHs with the performance of optimal absorbing structures, which is limited by the causality constraint [31]. If we calculate the minimum width of absorbing structures with the reflection coefficient spectrum $\left|\mathcal{R}\right(f\left)\right|$ as

$$L_{\min }=-\frac{c_0}{2{\pi }^2}{\int }_0^{f_{\mathrm{cut}-\mathrm{on}}}\mathrm{}\ln \left|\mathcal{R}\right(f\left)\right|\frac{\mathrm{d}f}{f^2},$$(13)

we get for the ABHs with N = 500, 400, 300, whose spectra are shown in Figure 12, L_min = 19.6 cm, 21.4 cm, and 23.3 cm, respectively, which, especially for N = 300, is close to the ABHs’ length L = 25 cm.

As it has been mentioned earlier, in all the previous numerical examples involving the optimization, only the ABH shapes were optimized, and the other parameters (L, R_min, H, h_rib) were fixed. Nevertheless, it might be useful to determine the optimal values of these additional parameters as well.

As an example, the ABH shape, described by the 5-dimensional parameter vector (see Eq. (12)), as well as the other parameters, were searched for within the following limits, following from the practical considerations, and the validity of the mathematical model, see Section 4.2: L ∈ 〈15 cm, 30 cm〉, R_min ∈ 〈0.25 mm, 1.00 mm〉, H ∈ 〈0.5 mm, 1.0 mm〉, and h_rib ∈ 〈0.25 mm, H − 0.25 mm〉.

The optimization algorithm converged to the following geometrical parameters, which are given to three significant figures: L = 30.0 cm, R_min = 0.250 mm, H = 0.999 mm, h_r = 0.253 mm. Consistently with the previous discussion, the total ABH length L and the spatial period of ribs H were maximized, and the minimum ABH radius R_min, and the thickness of the ribs h_rib were minimized. The resulting mean value of the reflection coefficient modulus $\langle \left|\mathcal{R}\right|\rangle =0.0694$ is slightly lower than the one obtained for the 30 cm-long ABH from Figure 10b, as in this case, the parameters H and ξ have higher values.

Figure 13 shows the reflection coefficient modulus spectrum calculated employing the 1D model (red line), and employing the 2D FEM model for the validation (the shape function, see the figure inset, results from the optimization based on the 1D model). Similarly as it has been shown in Figures 5d and 12, there is some discrepancy between the results at the middle-frequency range. Nevertheless, the relative difference between the values of $\langle \left|\mathcal{R}\right|\rangle$ is less then 0.15%. As predicted by equation (13), the minimum length of an absorber with this reflection coefficient spectrum is L_min = 29.2 cm.

Figure 13 Reflection coefficient modulus spectrum of an ABH, for which its shape as well as the other parameters (L, Rmin, H, hrib) were optimized; red line–prediction employing the Riccati equation, blue line–validation employing the 2D FEM model. Inset figure shows the corresponding optimal shape function.

5 Conclusions

This study provides a comprehensive numerical investigation into the properties of acoustic black holes, which are used as an anechoic termination in circular waveguides. Our research was focused on how the absorbing performance of ABHs is influenced by their profile and other parameters. We compared numerical predictions by a 1D model that uses the Riccati equation with a detailed 2D model based on the linearised Navier-Stokes equations and FEM. Our results demonstrate good agreement among the numerical models, especially for ABHs with a fine internal structure that corresponds to the original structure proposed by Mironov and Pislyakov [1]. The key components of the 1D model that enables this good agreement are accounting for the resonance in slits, boundary-layer losses in the slits, and boundary-layer losses in the ABH’s main cavity.

We quantified the absorbing performance of ABHs by calculating the mean value of the reflection coefficient modulus. Our results show that if the ABH profile is described by a power-law function, a profile close to linear provides the best results. Additionally, we found that the absorbing performance slightly improves with decreasing the minimum radius of the ABH and improves more rapidly with increasing the total length. The optimum profile for all these studied cases was found to be close to a linear one.

We also used a set of control points interconnected with cubic splines to allow for greater variability in ABH shape-function definition. Employing an evolutionary algorithm, we demonstrated that this increased variability can result in the reduction of the mean value of the reflection coefficient modulus by tens of percent compared with the optimum power-law shape-function, but with the same parameters. The number of ribs forming the ABH and the ratio of the widths of ribs and slits also play an important role, and better results are achieved with narrower ribs and wider slits. The numerical results indicate that the optimum shape functions are relatively simple, have derivatives close to zero at their extremities, and only weakly depend on the other parameters. The shape optimization significantly reduces the reflection coefficient at higher frequencies, which is caused by perfect impedance matching of the optimized ABH to the external waveguide.

Acoustic black holes provide an interesting and unconventional way of implementing passive anechoic termination in ductwork, as they do not require fibrous materials and can be used in harsh environments. Structures similar to those described in this work may soon be achievable with additive manufacturing techniques. Accurate and computationally efficient mathematical models, in conjunction with optimization methods, can help to determine the performance limits of ABHs and guide their further development. We hope that this work will contribute to this pursuit.

Funding

The authors are grateful for financial support from GACR grant GA22-33896S.

Conflicts of interest

The authors declare no conflict of interest in relation to this article.

Data availability statement

Data are available on request from the authors.

References

All Tables

All Figures

	Figure 1 Model of an acoustic black hole.
In the text

	Figure 2 Detail of an acoustic black hole.
In the text

	Figure 3 Model of an acoustic black hole.
In the text

	Figure 4 ABH shape function defined using cubic splines with 3 control points and 5 parameters.
In the text

	Figure 5 Reflection coefficient spectra of linear ABHs calculated employing the 1D model based on the Riccati equation (blue lines), and 2D finite element model (red lines), for different parameters: (a) −N = 500, ξ = 0.5, (b) −N = 250, ξ = 0.25, (c) −N = 250, ξ = 0.5, (d) −N = 250, ξ = 0.75.
In the text

	Figure 6 Mean value of the reflection coefficient modulus as a function of the power-law exponent for different values of Rmin; N = 500, ξ = 0.5.
In the text

	Figure 7 Panel (a) – mean value of the reflection coefficient modulus as a function of the power-law exponent α for different values of the ABH length L; Panel (b) – reflection coefficient spectra for different values of L and α; H = 0.5 mm, ξ = 0.5.
In the text

	Figure 8 Panel (a) – cubic-spline shape functions of optimized ABHs, Panel (b) – reflection coefficient modulus for optimized ABH (cubic-splines, 5 pars.), and linear ABH for comparison; N = 500, ξ = 0.5, Rmin = 0.25 mm.
In the text

	Figure 9 Panel (a) – shape functions, and Panel (b) – reflection coefficients of ABHs optimized in full frequency range (0 Hz–fcut-on) and low frequency range (0 Hz–1000 Hz); N = 500, ξ = 0.5, Rmin = 0.25 mm.
In the text

	Figure 10 Panel (a) – shape functions, and Panel (b) – reflection coefficient spectra of optimized ABHs with different lengths; H = 0.5 mm, ξ = 0.5, Rmin = 0.25 mm.
In the text

	Figure 11 Mean value of the reflection coefficient modulus of ABHs optimized for different number of ribs and different internal-structure parameters; L = 25 cm, Rmin = 0.25 mm.
In the text

	Figure 12 Reflection coefficient modulus as a function of frequency for optimized ABHs with different number of ribs N; solid lines–prediction employing the Riccati equation, dashed lines–validation employing 2D FEM model. In all the cases, L = 25 cm, hrib = 0.25 mm, Rmin = 0.25 mm. Inset figure shows the corresponding shape functions.
In the text

	Figure 13 Reflection coefficient modulus spectrum of an ABH, for which its shape as well as the other parameters (L, Rmin, H, hrib) were optimized; red line–prediction employing the Riccati equation, blue line–validation employing the 2D FEM model. Inset figure shows the corresponding optimal shape function.
In the text