Topology optimization of a waveguide acoustic black hole for enhanced wave focusing

(Color online) (a) A uniform beam. (b) An ideal ABH where the thickness of the beam goes to zero according to a power-law profile. (c) A practical design when the ideal profile is truncated at a certain thickness and a small amount of damping material (the brown part) is used at the tip of the beam to efficiently attenuate the vibrations. (d) An interpretation of the ABH profile as a particular design of the gray domain $Ω^{D}$ that leads to optimal use of the small amount of damping material at the end of the beam.

FIG. 2.

(Color online) (a) A 3D visualization of the optimization problem where a semi-infinite cylindrical waveguide is considered with an absorbing boundary condition imposed in a small area at the far end of the waveguide. (b) Axisymmetric 2D model of the same semi-infinite cylindrical waveguide.

A. Mathematical model and discretization

Consider the axisymmetric setup illustrated in Fig. 3. The design domain

Ω^{D}

⁠, the region marked with a gray background mesh, may partly be filled with sound-hard material. The aim is to find the distribution of solid material

Ω^{s}

within the design domain that maximizes sound absorption at the lossy boundary

Γ^{L}

marked in brown. The volume outside of

Ω^{D}

is assumed to be filled with air. Consider time-harmonic linear wave propagation in the air-filled region

\hat{Ω}

of the waveguide, and let

P (x, t) = Re {p (x) e^{i ω t}}

denote the acoustic pressure, where ω is the angular frequency and p is the complex amplitude function satisfying Helmholtz equation in cylindrical coordinates

(x = {(r, z)}^{T})

⁠,

- \nabla \cdot (r \nabla p) - k^{2} r p = 0, in and \hat{Ω},

(1)

where

\nabla = (\partial / \partial r, \partial / \partial z), k = ω / c

the wave number, c is the speed of sound in air, and

\hat{Ω}

is the air-filled region of the waveguide. As mentioned earlier, the cylindrical waveguide considered in the optimization problem is infinite to the left. For numerical purposes, we truncate the inlet at artificial boundary

Γ^{in}

as illustrated in Fig. 3.

FIG. 3.

(Color online) The computational domain consists of the truncated waveguide with an arbitrary distribution of sound-hard material in the design domain.

Previous studies on the WAB have shown the significance of the viscous and thermal losses on the wave propagation characteristics of the classical ribbed design (Červenka and Bednařík, 2022; Guasch , 2017; Guasch , 2020; Mi , 2021; Mironov and Pislyakov, 2020; Mousavi , 2022; Umnova , 2023). Thus, we expect that viscothermal losses will play a crucial role in wave propagation also in the design domain for the optimization problem. Here, we specify the optimization formulation considering viscothermal losses in the solid–fluid interface $Γ^{w}$ in the design domain. Later on, in the numerical experiments, we present the results of the optimization problem both for the lossy and lossless cases. A comparison between these cases is provided in the Discussion section.

An accurate model of acoustics in the presence of viscothermal losses is provided by finite-element solutions of the LNS equations [instead of the Helmholtz Eq. (1)]. Unfortunately, this model is ill-suited to design optimization for two reasons. First, the computational cost is high, both in terms of processing time and memory requirements. Second, to limit the computational cost to a reasonable level, it is necessary to use a very aggressive mesh refinement strategy to resolve the typically exceedingly thin viscous and thermal boundary layers close to solid walls. For instance, in the audio regime 20 Hz–20 kHz, the boundary layers are smaller than the wavelength by a factor on the order of $10^{- 5}$ – $10^{- 3}$ ⁠. It would be complicated and computationally expensive to manage and adjust the required extreme mesh refinements in a design optimization context.

To model the effect of viscothermal losses, we therefore instead utilize the boundary condition devised by Berggren (2018). This boundary condition has been shown to provide accurate results—particularly when the geometry is such that the boundary layers do not overlap—to a fraction of the computational cost of solving the LNS equations and without the need for boundary-layer refinements (Andersen , 2023; Berggren , 2018; Billard , 2021). However, to ensure confidence in the computations, we will consistently validate the performance of the final designs using the LNS simulations.

Regarding the rest of the boundaries, the acoustic conditions at the lossy wall

Γ^{L}

will be modeled by the non-dimensional specific admittance

Y

⁠, and a unit-amplitude plane incoming wave will be imposed at

Γ^{in}

⁠. Altogether, we obtain the boundary value problem

\nabla \cdot (r \nabla p) + k^{2} r p = 0 in \hat{Ω},

(2a)

\frac{\partial p}{\partial n} = 0 on Γ^{s},

(2b)

\frac{\partial p}{\partial n} + i k p = 2 i k on Γ^{in},

(2c)

- δ_{V} \frac{i - 1}{2} \nabla_{T} \cdot (r \nabla_{T} p) + δ_{T} k^{2} \frac{(i - 1) (γ - 1)}{2} r p + r \frac{\partial p}{\partial n} = 0 on Γ^{w},

(2d)

\frac{\partial p}{\partial n} + ikYp = 0 on Γ^{L},

(2e)

where

δ_{V} = \sqrt{2 ν / ω}

and

δ_{T} = \sqrt{2 κ / ω ρ_{0} c_{p}}

are the viscous and thermal boundary-layer thicknesses, respectively. Material parameters ν, κ, and c_p are the kinematic viscosity, the thermal conductivity, and the specific heat capacity for air at constant pressure, respectively. Moreover,

\nabla_{T} = \nabla - \partial / \partial n (n)

is the tangential gradient operator, in which n is the outward directed unit normal. Boundary condition models the viscous and thermal losses in the boundary layer at the solid–fluid interface

Γ^{w}

in the design domain. Artificial boundary condition specifies a plane right-going wave with unit amplitude while absorbing the left-going plane waves at

Γ^{in}

⁠. We note that in the numerical experiments, the dimensions of the device and the operating range of frequencies are chosen such that the plane wave is the only propagating mode in the waveguide.

Boundary condition (2e) models a wall with a non-dimensional specific admittance Y at the end of the waveguide on $Γ^{L}$ ⁠. The value of Y depends on the specific application. For instance, $Y = 1$ represents a plane wave absorber; a sensor or a non-ideal absorber can be modeled by a frequency-dependent complex admittance. Alternatively, $- i k Y$ can also denote a Dirichlet-to-Neumann (DtN) operator interfacing $Γ^{L}$ to a region outside the current computational domain, for instance, free space or a semi-infinite damped pipe. The latter case corresponds to a perfect absorbing boundary condition, as defined, for instance, in the article by Wadbro (2014). Such a DtN operator ensures that the outgoing waves are perfectly absorbed by the boundary, mimicking an idealized absorbing material.

The standard variational form of problem (2) reads as follows.

Find p such that

\begin{matrix} \int_{\hat{Ω}} r \nabla q \cdot \nabla p d Ω - k^{2} \int_{\hat{Ω}} r q p d Ω + i k \int_{Γ^{in}} r q p d Γ + i k \int_{Γ^{L}} r q Y p d Γ + δ_{V} \frac{i - 1}{2} \int_{Γ^{w}} r \nabla_{T} q \cdot \nabla_{T} p d Γ + δ_{T} k^{2} \frac{(i - 1) (γ - 1)}{2} \int_{Γ^{w}} r q p d Γ \\ = 2 i k \int_{Γ^{in}} r q d Γ \forall q . \end{matrix}

(3)

For details on the required test and trial space, which is a strict subspace of

H^{1} (\hat{Ω})

⁠, we refer to the publication by Appendix B in Berggren (2018).

To describe the solid material distribution in the design domain $Ω^{D}$ ⁠, we define the material indicator function $α : Ω \to {0, 1}$ such that $α \equiv 0$ in the solid region $Ω^{s}$ and $α \equiv 1$ in the air-filled region $\hat{Ω}$ ⁠. Here, $Ω = \hat{Ω} \cup Ω^{s}$ denotes the computational domain containing both the air and solid regions. Using material indicator function α, we rewrite problem (3) on the computational domain as follows:

Find p such that

\begin{matrix} \int_{Ω} α r \nabla q \cdot \nabla p d Ω - k^{2} \int_{Ω} α r q p d Ω + i k \int_{Γ^{in}} r q p d Γ + i k \int_{Γ^{L}} r q Y p d Γ \\ + δ_{V} \frac{i - 1}{2} \int_{Γ^{w} (α)} r \nabla_{T} q \cdot \nabla_{T} p d Γ + δ_{T} k^{2} \frac{(i - 1) (γ - 1)}{2} \int_{Γ^{w} (α)} r q p d Γ \\ = 2 i k \int_{Γ^{in}} r q d Γ \forall q . \end{matrix}

(4)

Note that for a given design represented by function α, we need to identify the air–solid interface and compute the two boundary integrals on

Γ^{w}

⁠, the fifth and sixth terms in Eq. .

We use the finite element method to discretize and solve problem (4) on a uniform mesh of square elements. Let $V$ be the space of continuous functions that are bi-quadratic on each element. Also, let $φ_{j}, j = 1, 2, \dots, N$ ⁠, be the nodal basis functions of $V$ such that $V = span {φ_{1}, φ_{2}, \dots, φ_{N}}$ ⁠, where N is the number of degrees of freedom in the finite element approximation. By approximating the complex pressure p, test function q and indicator function α by $p_{h} \in V, q_{h} \in V$ ⁠, and the element-wise constant function α_h, respectively, we obtain the discrete form of problem (4) as follows.

Find p_h such that

\begin{matrix} \int_{Ω} α_{h} r \nabla q_{h} \cdot \nabla p_{h} d Ω - k^{2} \int_{Ω} α_{h} r q_{h} p_{h} d Ω + i k \int_{Γ^{in}} r q_{h} p_{h} d Γ + i k \int_{Γ^{L}} r q_{h} Y p_{h} d Γ \\ + δ_{V} \frac{i - 1}{2} \int_{Γ^{w} (α_{h})} r \nabla_{T} q_{h} \cdot \nabla_{T} p_{h} d Γ + δ_{T} k^{2} \frac{(i - 1) (γ - 1)}{2} \int_{Γ^{w} (α_{h})} r q_{h} p_{h} d Γ \\ = 2 i k \int_{Γ^{in}} r q_{h} d Γ \forall q_{h} . \end{matrix}

(5)

To compute the boundary integrals on

Γ^{w}

⁠, we define Θ as the set of all edges in the finite element mesh. At each mesh edge

θ \in Θ

⁠, the squared jump

{[[α_{h} |_{θ}]]}^{2}

of the element-wise constant function α_h is well defined and indicates, when unity, the presence of an air–solid boundary. Using this definition, we extend the integration domain of the two mentioned boundary integrals in Eq. from

Γ^{w}

⁠, the air–solid interface, to Θ, the set of all edges in the mesh. Then, we obtain the extended discrete problem as follows.

Find p_h such that

\begin{matrix} \int_{Ω} α_{h} r \nabla q_{h} \cdot \nabla p_{h} d Ω - k^{2} \int_{Ω} α_{h} r q_{h} p_{h} d Ω + i k \int_{Γ^{in}} r q_{h} p_{h} d Γ + i k \int_{Γ^{L}} r q_{h} Y p_{h} d Γ \\ + δ_{V} \frac{i - 1}{2} \int_{Θ} {[[α_{h}]]}^{2} r \nabla_{T} q_{h} \cdot \nabla_{T} p_{h} d Γ + δ_{T} k^{2} \frac{(i - 1) (γ - 1)}{2} \int_{Θ} {[[α_{h}]]}^{2} r q_{h} p_{h} d Γ \\ = 2 i k \int_{Γ^{in}} r q_{h} d Γ \forall q_{h} . \end{matrix}

(6)

By the definition of the indicator function α_h, for binary values of the material indicator function, the solution p_h to problem is well defined and equal to the solution to problem in

\hat{Ω}

but is undefined in

Ω^{s}

⁠. Following a standard approach in material distribution topology optimization (Bendsøe and Sigmund, 2003; Bokhari , 2021; Kasolis , 2015; Wadbro, 2014; Wadbro and Berggren, 2006), we redefine the lower bound of the material indicator function such that

α_{h} = ϵ

⁠, where ϵ is a small positive number, in the solid region

Ω^{s}

⁠. In this way, the solution to problem is well defined in the whole computational domain

Ω = \hat{Ω} \cup Ω^{s}

⁠.

B. Power balance

Here, we show how different portions of the power of the unit-amplitude incoming wave are either reflected or dissipated in the waveguide. To this end, we choose

q_{h} = {\bar{p}}_{h}

in Eq. , where the bar signifies complex conjugation, and obtain

\begin{matrix} \int_{Ω} α_{h} r \nabla {\bar{p}}_{h} \cdot \nabla p_{h} d Ω - k^{2} \int_{Ω} α_{h} r {\bar{p}}_{h} p_{h} d Ω + i k \int_{Γ^{in}} r {\bar{p}}_{h} p_{h} d Γ + i k \int_{Γ^{L}} r {\bar{p}}_{h} Y p_{h} d Γ \\ + δ_{V} \frac{i - 1}{2} \int_{Θ} {[[α_{h}]]}^{2} r \nabla_{T} {\bar{p}}_{h} \cdot \nabla_{T} p_{h} d Γ + δ_{T} k^{2} \frac{(i - 1) (γ - 1)}{2} \int_{Θ} {[[α_{h}]]}^{2} r {\bar{p}}_{h} p_{h} d Γ \\ = 2 i k \int_{Γ^{in}} r {\bar{p}}_{h} d Γ . \end{matrix}

(7)

Considering only the imaginary part of Eq. , we have

\begin{matrix} k \int_{Γ^{in}} r | p_{h} |^{2} d Γ + k \int_{Γ^{L}} r Re {{\bar{p}}_{h} Y p_{h}} d Γ + \frac{δ_{V}}{2} \int_{Θ} {[[α_{h}]]}^{2} r | \nabla_{T} p_{h} |^{2} d Γ \\ + \frac{(δ_{T} k^{2}) (γ - 1)}{2} \int_{Θ} {[[α_{h}]]}^{2} r | p_{h} |^{2} d Γ = 2 k \int_{Γ^{in}} r Re {p_{h}} d Γ, \end{matrix}

(8)

where

Re {\cdot}

denotes the real part of a complex quantity. By the properties of complex numbers, we rewrite Eq. , as

\begin{matrix} \int_{Γ^{in}} r | p_{h} - 1 |^{2} d Γ - \int_{Γ^{in}} r | 1 |^{2} d Γ + \int_{Γ^{L}} r Re {{\bar{p}}_{h} Y p_{h}} d Γ + \frac{δ_{V}}{2 k} \int_{Θ} {[[α_{h}]]}^{2} r | \nabla_{T} p_{h} |^{2} d Γ \\ + \frac{(δ_{T} k) (γ - 1)}{2} \int_{Θ} {[[α_{h}]]}^{2} r | p_{h} |^{2} d Γ = 0, \end{matrix}

(9)

which can also be written

\begin{matrix} \underset{Input}{\underset{︸}{\int_{Γ^{in}} r | 1 |^{2} d Γ}} = \underset{Reflected}{\underset{︸}{\int_{Γ^{in}} r | p_{h} - 1 |^{2} d Γ}} + \underset{Dissipated at Γ^{L}}{\underset{︸}{\int_{Γ^{L}} r Re {{\bar{p}}_{h} Y p_{h}} d Γ}} \\ + \underset{Dissipated at Γ^{w} as visco - thermal losses}{\underset{︸}{\int_{Θ} {[[α_{h}]]}^{2} (\frac{δ_{V}}{2 k} r | \nabla_{T} p_{h} |^{2} + δ_{T} k \frac{(γ - 1)}{2} r | p_{h} |^{2}) d Γ}} . \end{matrix}

(10)

The balance law shows that the power of the imposed unit-amplitude incoming wave at

Γ^{in}

is split into three portions: the power of the reflected wave on

Γ^{in}

⁠, the viscothermal losses on

Γ^{w}

⁠, and the dissipated power at the lossy boundary

Γ^{L}

⁠. If we divide the power balance with the left-side term, we obtain the normalized power balance

W_{refl}^{norm} + W_{vth}^{norm} + W_{L}^{norm} = 1,

(11)

where

W_{refl}^{norm}

is the power of the reflected wave,

W_{vth}^{norm}

is the viscothermal power loss on

Γ^{w}

⁠, and

W_{L}^{norm}

is the power loss at the lossy boundary

Γ^{L}

⁠, all normalized with the power of the incoming wave.

C. Optimization problem

As mentioned earlier, our aim is to optimize the layout of solid material in the design domain to maximize sound absorption at the lossy boundary

Γ^{L}

⁠. Thus, we aim to maximize

W_{L}^{norm}

⁠, which, due to power balance , is equivalent to minimize the sum of

W_{refl}^{norm}

and

W_{vth}^{norm}

⁠. We define the primary objective function as

\begin{matrix} J_{p} (α_{h}) = \int_{Γ^{in}} r | 1 |^{2} d Γ \sum_{i = 1}^{N_{freq}} (W_{refl}^{norm} + W_{vth}^{norm}) \\ = \sum_{i = 1}^{N_{freq}} (\int_{Γ^{in}} r | p_{h} (k_{i}) - 1 |^{2} d Γ + \int_{Θ} {[[α_{h}]]}^{2} (\frac{δ_{V}}{2 k} r | \nabla_{T} p_{h} (k_{i}) |^{2} + δ_{T} k \frac{(γ - 1)}{2} r | p_{h} (k_{i}) |^{2}) d Γ), \end{matrix}

(12)

where

N_{freq}

is the number of frequencies considered in the optimization, and k_i is the wave number corresponding to the ith frequency subject to optimization.

To be able to use gradient-based algorithms to solve the optimization problem with objective function (12), we allow the indicator function α_h to take values in the range $[ϵ, 1]$ ⁠, following a standard approach in material distribution topology optimization (Bokhari , 2021; Dühring , 2008; Kasolis , 2015; Wadbro, 2014; Wadbro and Berggren, 2006). Nevertheless, we are interested in a final design with a pure solid–air layout, free of intermediate values of α_h. Thus, we use a combination of nonlinear filters and penalty methods to suppress intermediate values of α_h.

Let d be the design variables, a vector of size $M \times 1$ that holds the element values of the material distribution function before filtering, where M is the number of elements in the mesh. We define $M \times 1$ vector $α_{h} = {[α_{1}, α_{2}, \dots, α_{M}]}^{T} : = F_{ρ} (d)$ holding the element values of α_h for a given vector of design variables d, where $F_{ρ}$ is a filtering operator with the filter radius $ρ$ ⁠. In this study, we use the nonlinear filters as described by Bokhari (2021).

Moreover, to suppress intermediate values of α_h, we utilize a standard quadratic penalty term

\frac{γ}{M} \sum_{l = 1}^{M} (α_{l} - ϵ) (1 - α_{l}),

(13)

where γ is the penalty parameter. Thus, the optimization problem in the discrete case reads

\begin{matrix} min_{d \in A} J_{p} (\underset{α_{h}}{\underset{︸}{F_{ρ} (d)}}) + \frac{γ}{M} {(F_{ρ} (d) - ϵ 1_{M})}^{T} (1_{M} - F_{ρ} (d)), \\ A = {d \in ℝ^{M} | ϵ \leq d_{i} \leq 1 \forall i}, \end{matrix}

(14)

where J_p is the primary objective function ,

1_{M}

is

M \times 1

vector with all entries equal to 1, and

A

is the set of admissible designs.

III. NUMERICAL EXPERIMENTS

Here, we present the results of solving design optimization problem (14) for the configuration illustrated in Fig. 3 with $R = 115 mm, L = 255 mm, L^{in} = 155 mm$ ⁠, and $r^{L} = 15 mm$ ⁠, which are the same dimensions used in previous experimental works by El Ouahabi (2015a,b) and in the numerical work by Mousavi (2022). Motivated by these previous studies that adopted the same setup, we aim to maximize the power loss at the small boundary $Γ^{L}$ in the frequency range 400–1000 Hz, considering 13 frequencies spaced 50 Hz apart. That is, $N_{freq} = 13$ and $k_{i} = 2 π (350 + 50 i) / c$ ⁠, for $i = 1, 2, \dots, 13$ ⁠. We note that for the chosen radius of the waveguide, the cut-on frequency of the first non-planar circumferential mode is $f_{c} = 866 Hz$ ⁠, while the cut-on frequency of the first non-planar axisymmetric mode is $1803 Hz$ ⁠. We assume a perfect absorber at $Γ^{L}$ by letting $- i k Y$ be the Dirichlet-to-Neumann (DtN) map for a semi-infinite pipe attached to $Γ^{L}$ ⁠. For more details about DtN type boundary conditions, the reader is referred to the work by Wadbro (2014) or the book by Ihlenburg (1998). This boundary condition models a small perfect absorbing boundary at the end of the waveguide.

We set the lower bound

ϵ = 10^{- 8}

for the design variable, the filter radius to 2 mm, and use a so-called continuation approach for the penalty and non-linearity parameters in the filters. More precisely, we solve problem for a sequence of increasing penalty parameters

γ_{i} = 10^{i}, i = 0, 1, \dots, 5

using the previously computed solution as initial design for the next step. In this way, we gradually move the focus of the optimizer from mainly minimizing the primary objective towards obtaining a pure solid–air design, free of intermediate values of α_h (Bokhari , 2021; Wadbro, 2014; Wadbro and Berggren, 2006). The convergence criterion for the optimization algorithm is based on the residual norm of the Karush–Kuhn–Tucker (KKT) (or the first-order optimality) condition, together with a bound on the number of iterations in each penalty step. At the last step in the continuation approach, a flood-fill algorithm is used to fill up isolated air regions inside the solid parts resulting from the distinct steps in the continuation approach. We also use the normalized power dissipated at the lossy boundary

Γ^{L}

⁠,

W_{L}^{norm} = \frac{\int_{Γ^{L}} r Re {{\bar{p}}_{h} Y p_{h}} d Γ}{\int_{Γ^{in}} r | 1 |^{2} d Γ},

(15)

to compare the performance of different designs.

To demonstrate the effect of viscothermal losses modeled by boundary condition (2d) on the design optimization problem, we will consider a lossless case as well as a lossy case. All computations, involving solutions of the state and adjoint equations as well as the optimization iterations, are carried out using a custom in-house Matlab code. Additionally, all final results are cross-validated using Comsol Multiphysics.

A. Lossless case

In the lossless case, we neglect the viscothermal losses in the device, which corresponds to setting

δ_{V} = δ_{T} = 0

in boundary condition , leading to the vanishing Neumann condition also on the solid–air interface

Γ^{w}

⁠. The terms representing viscothermal losses in primary objective function will also be neglected. Substituting the primary objective function J_p (neglecting viscothermal losses) into problem , we obtain the optimization problem

\begin{matrix} min_{d \in A} \sum_{i = 1}^{N_{freq}} \int_{Γ^{in}} r | p_{h} (k_{i}) - 1 |^{2} d Γ + \frac{γ}{M} {(F_{ρ} (d) - ϵ 1_{M})}^{T} (1_{M} - F_{ρ} (d)), \\ A = {d \in ℝ^{M} | ϵ \leq d_{i} \leq 1 \forall i} \end{matrix}

(16)

in the lossless case.

We employ the least squares formulation of the method of moving asymptotes (MMA), developed by Svanberg (1987) to solve the optimization problem (16). Figure 4 shows the optimized design achieved by solving design optimization problem (16) considering state Eq. (6) with $δ_{V} = δ_{T} = 0$ ⁠. The final design has similarities to the classical ribbed design of WABs. That is, it consists of a few retarding radial cavities getting smaller toward the end of the waveguide; however, some of these cavities are connected through narrow openings. Figure 5 shows the performance of the device presented in Fig. 4. This performance is calculated by importing the final design into Comsol Multiphysics (the “pressure acoustics interface” in the acoustics module) and numerically solving Helmholtz Eq. (3) with $δ_{V} = δ_{T} = 0$ using a classical boundary-fitted mesh. As illustrated in Fig. 5, on average, more than 96% of the power of the incoming wave in the targeted range of frequencies is being dissipated in the lossy boundary $Γ^{L}$ ⁠. Note, however, that this result is achieved by neglecting viscothermal losses in the device.

FIG. 4.

(Color online) (a) A 3D visualization of the optimized design in the lossless case. (b) Cross section of the optimized design in the lossless case.

FIG. 5.

(Color online) The normalized power portions for the optimized design in the lossless case illustrated in Fig. 4. The red dashed line is the normalized power of the reflected wave and the blue line is the normalized power loss at the lossy boundary $Γ^{L}$ ⁠, represented as $W_{refl}^{norm}$ and $W_{L}^{norm}$ in Eq. (11), respectively. Note that we here have neglected the viscothermal losses $W_{vth}^{norm}$ on $Γ^{w}$ ⁠. The targeted range of frequencies from 400 to 1000 Hz is highlighted in gray.

A more realistic model incorporates viscothermal losses modeled by applying boundary condition (2d) at the solid inclusion boundaries, denoted as $Γ^{w}$ ⁠. The simulation results, depicted in Fig. 6, are compared to those obtained from a linearized Navier–Stokes (LNS) simulation. (All LNS computations are carried out in Comsol Multiphysics.) The results reveal that a significant portion of the power of the incoming wave is dissipated as a result of viscothermal losses. Consequently, the performance of the device, defined as having high power loss at $Γ^{L}$ ⁠, exhibits a notable decrease within the targeted frequency range. As a consequence, the lossless optimized design achieves a suboptimal wave-focusing effect. This underscores the crucial role of considering viscothermal losses during the design process for this device.

FIG. 6.

(Color online) The normalized power portions for the optimized design in the lossless case illustrated in Fig. 4. The red dashed line is $W_{refl}^{norm}$ ⁠, the blue line $W_{L}^{norm}$ ⁠, and the gray dash-dotted line $W_{vth}^{norm}$ ⁠, resulting from simulations with the acoustic boundary layer (ABL) model for the viscothermal losses. The green triangles, magenta squares, and gray diamonds represent $W_{refl}^{norm}, W_{L}^{norm}$ ⁠, and $W_{vth}^{norm}$ ⁠, respectively, resulting from linearized Navier–Stokes (LNS) simulations. The targeted range of frequencies from 400 to 1000 Hz is highlighted in gray.

Figure 7 shows the convergence history for the lossless optimization. Note that the increase in the objective function during the late stages of optimization is due to the continuation strategy of increasing the penalty parameter γ to promote integer design variables, leading to an almost pure solid–air design at termination. Furthermore, a significant manufacturing challenge encountered in realizing the final design, as depicted in Fig. 4, relates to the positioning of free-floating solid inclusions within the device. Note that these free-hanging parts in the axially symmetric context, appear as rings in a full three-dimensional (3D) representation as depicted in Fig. 4(a). A feasible approach to address this challenge involves supporting these inclusions using thin strings or bars connected to the outer tube, a strategy suggested by Mousavi (2023). Their work demonstrates that these slender connecting components minimally impact the acoustic performance of the device. Alternatively, an adjustment to the objective function in the optimization problem could be considered to ensure the connectivity of the design inclusions. However, it is crucial to acknowledge that such a modification would significantly limit the available design space for the optimizer in the axially symmetric context, potentially compromising the acoustic performance of the final design. Moreover, employing thin connecting components for structural reinforcement is a common practice in acoustics. For instance, similar approaches are prevalent in providing structural support for loudspeaker frames. Therefore, prioritizing the use of these slender connecting parts is recommended as the primary approach to address the mentioned manufacturing challenge.

FIG. 7.

(Color online) Convergence history for the lossless optimization. The top row shows designs corresponding to different iterations, and the bottom graph shows the value of the primary objective function $J_{p} (α_{h})$ versus the number of iterations.

B. Lossy case

In the previous section, we observed that when neglecting viscothermal boundary losses during the optimization process, we obtain structures with suboptimal wave-focusing effects. To address this limitation, we now consider optimization problem (14) with the primary objective function J_p defined in Eq. (12) with nonzero δ_V and δ_T. This means that our objective is to minimize not only the reflections, as in the lossless case, but also the viscothermal boundary losses as described by Eq. (12). To solve the optimization problem, we use a similar MMA approach as in the lossless case. The gradients necessary for the optimization are obtained through solutions of the corresponding adjoint equations, as outlined in Appendix.

Figure 8 shows the optimized design obtained in the lossy case. Unlike the lossless case, the optimizer appears to limit the perimeter of the solid–air interface and avoid structures with narrow necks in order to minimize the viscothermal boundary losses. The optimized design still incorporates retarding radial cavities, but with a modified structure, where the shape of the cavities is tuned to mitigate local radial resonances and maximize the wave-focusing effect. Figure 9 shows the performance of the optimized design in Fig. 8. Note that in this case, the optimizer succeeded in creating a device that maximizes the power loss in the targeted boundary $Γ^{L}$ by minimizing both the reflections and the viscothermal losses in the targeted range of frequencies. Moreover, LNS simulations validate the optimal wave-focusing effect of the optimized device. Figure 10 illustrates the convergence history of the optimization process for the lossy case. The optimization is initiated with an empty waveguide, and the graph shows, as the design is gradually formed, the progression of the primary objective function (12) value over iterations, indicating the convergence of the optimization algorithm towards an optimal design that minimizes both reflections and viscothermal boundary losses in the device. It is worth noting that the convergence history plot, as shown in Fig. 10, exhibits multiple peaks that correspond to the steps taken in the penalty parameter γ due to the application of a continuation strategy with increasing penalty values. This strategy plays a crucial role in guiding the optimization process towards pure solid–air designs.

FIG. 8.

(Color online) (a) A 3D visualization of the optimized design in lossy case. (b) Cross section of the optimized design in the lossy case.

FIG. 9.

(Color online) The normalized power portions for the optimized design in the lossy case illustrated in Fig. 8. The red dashed line is $W_{refl}^{norm}$ ⁠, the blue line $W_{L}^{norm}$ ⁠, and the gray dash-dotted line $W_{vth}^{norm}$ ⁠, resulting from simulations with the acoustic boundary layer (ABL) modeling for the viscothermal losses. The green triangles, magenta squares, and gray diamonds represent $W_{refl}^{norm}, W_{L}^{norm}$ ⁠, and $W_{vth}^{norm}$ ⁠, respectively, resulting from linearized Navier–Stokes (LNS) simulation. The targeted range of frequencies from 400 to 1000 Hz is highlighted in gray.

FIG. 10.

(Color online) History of convergence for the lossy optimization. The top row shows designs corresponding to different iterations, and the bottom graph shows the values of the primary objective function $J_{p} (α_{h})$ versus the number of iterations.

IV. DISCUSSION

To compare the wave-focusing effect of the optimized devices, we conducted LNS simulations using a boundary-fitted mesh in COMSOL Multiphysics (denoted as “thermoviscous acoustics interface” in the acoustics module of COMSOL Multiphysics). In addition to the simulations of the optimized devices, we also performed simulations for a classical ribbed design of a WAB, as shown in Fig. 11. The ribbed design consists of 42 rings with a thickness of 2 mm, evenly spaced 4 mm apart, and with a linearly retarding inner radius of the rings. These simulations provide a comprehensive comparison between the optimized designs and the classical ribbed design in terms of their wave-focusing performance.

FIG. 11.

(Color online) (a) A 3D visualization of the classical ribbed design of WABs. (b) Axi-symmetrical model of the classical ribbed design of WABs.

Figure 12 presents the performance of the classical ribbed design. The wave-focusing effect of this design is poor, as indicated by the low values of normalized power dissipation at the targeted boundary $Γ^{L}$ ⁠. However, due to significant viscothermal boundary losses, the device exhibits minimal reflections. In fact, this behavior is the reason that in recent literature, this device is often referred to as a sound-absorbing structure rather than a WAB (Červenka and Bednařík, 2022; Mousavi , 2021, 2022). By comparing Figs. 6, 9, and 12, we conclude the following.

The classical ribbed design exhibits low reflections due to significant viscothermal losses in combination with localized radial resonances in the cavities between the rings.
The classical ribbed design demonstrates a minimal wave-focusing effect, indicating its limited ability to achieve an efficient acoustic black hole effect.
The optimized device in the lossless case shows improved wave focusing compared to the classical ribbed design. However, its performance is limited by extensive viscothermal losses due to large areas of solid material and the presence of narrow necks in the structure.
In contrast, the optimized device in the lossy case achieves superior wave-focusing in the targeted frequency range with minimum reflection and viscothermal boundary losses.
The optimized design in the lossy case outperforms, in terms of wave-focusing, both the classical ribbed design and the optimized device in the lossless case, making it a suitable candidate for a true WAB.

FIG. 12.

(Color online) The normalized power portions for the classical ribbed design illustrated in Fig. 11. The red dashed line is $W_{refl}^{norm}$ ⁠, the blue line $W_{L}^{norm}$ ⁠, and the gray dash-dotted line $W_{vth}^{norm}$ ⁠, resulting from simulations with the acoustic boundary layer (ABL) modeling for the viscothermal losses. The green triangles, magenta squares, and gray diamonds represent $W_{refl}^{norm}, W_{L}^{norm}$ ⁠, and $W_{vth}^{norm}$ ⁠, respectively, resulting from linearized Navier–Stokes (LNS) simulation. The targeted range of frequencies from 400 to 1000 Hz is highlighted in gray.

As mentioned previously, two key characteristics of acoustic black holes are the gradual increase in wave amplitude and local wavenumber as the wave propagates toward the end of the black hole profile. To provide a visual representation of this phenomenon in the context of WABs, we plot in Fig. 13 the real part (left column) and absolute value (right column) of the complex pressure for (a) the classical ribbed design, (b) the lossless optimized design, and (c) the lossy optimized design. The graph at the bottom of each subfigure shows the pressure at the center axis for four frequencies and the top image shows the pressure at frequency 1000 Hz. Comparing the results, it is evident that both optimized devices outperform the classical ribbed design in terms of achieving higher centerline pressure amplitudes towards the end of the waveguide. However, as the frequency increases, the optimized device in the lossless case exhibits lower amplitudes towards the end due to the significant increase in viscothermal losses in the narrow necks. On the other hand, consistently for all considered frequencies, the optimized device in the lossy case demonstrates the desired behavior of large pressure amplitudes at the end of the waveguide. This property indicates a successful representation of one of the key characteristics of acoustic black holes. Moreover, as evidenced by the shorter distances between consecutive peaks/deeps towards the end of the waveguide in the plots of the real part of the pressure, the behavior of the centerline pressure in the device optimized in the lossy case indicates a gradual increase in an effective wavenumber in the axial direction. This behavior further highlights the successful realization of the key characteristics of an acoustic black hole in the optimized design in the lossy case.

FIG. 13.

(Color online) On the left is the real part of the complex pressure amplitude $Re {p}$ ⁠, and on the right is the absolute value of the complex pressure amplitude $| p |$ ⁠. On top of each subfigure is a 2D plot at a frequency of 1000 Hz for the corresponding entity; on the bottom of each subfigure is the value of the corresponding entity on the symmetry axis for four different frequencies. (a) The classical ribbed design, (b) the lossless optimized design, (c) the lossy optimized design. The region inside the design domain, the black hole region, is highlighted in gray in the plots.

As discussed in Sec. II A, the choice of the admittance Y of the lossy wall at $Γ^{L}$ depends on the specific application, and it can represent various scenarios, such as a complex frequency-dependent value for a sensor or a porous absorber. Nevertheless, the primary objective of the WAB design is to achieve wave-focusing, regardless of the specific choice of the admittance of the wall at $Γ^{L}$ ⁠. To investigate the robustness of the optimized design with respect to different boundary conditions at $Γ^{L}$ ⁠, we conduct a final numerical experiment. This time, we modify the boundary condition at $Γ^{L}$ from a perfect absorber to yet another extreme case of a perfect hard scatterer by imposing a sound-hard boundary condition at $Γ^{L}$ for all three devices. With this change, we aim to observe the wave-focusing performance of each device under the hard scatterer scenario. This change is analogous to removing the damping material from the end of a beam acoustic black hole, allowing for perfect reflection of waves from the end of the profile. As a result, a standing wave pattern emerges, making it easier to observe changes in the local effective axial wavenumber (the distance between consecutive pressure nodes). Figure 14 presents the absolute value of the complex pressure at the center axis for four frequencies, along with the corresponding two-dimensional (2D) pressure amplitude plots at 1000 Hz displayed at the top for each case. In this scenario, where a sound-hard boundary condition is imposed at $Γ^{L}$ ⁠, the gradual increase in the local effective axial wavenumber is more apparent and can be easily observed as the decrease in distance between consecutive pressure nodes.

FIG. 14.

https://doi.org/10.1016/j.jsv.2022.117455

(Color online) On the top, the absolute value of the complex pressure amplitude $| p |$ at frequency 1000 Hz for the case where sound-hard boundary condition is imposed at $Γ^{L}$ ⁠; on the bottom, $| p |$ on the symmetry axis for four different frequencies in cases (a) the classical ribbed design, (b) the lossless optimized design, (c) the lossy optimized design considering sound-hard boundary condition at $Γ^{L}$ ⁠. The region inside the design domain, the black hole region, is highlighted in gray in the plots.

V. CONCLUSION

In this study, we investigated the feasibility of achieving a wave-focusing acoustic black hole effect in waveguides through computational design optimization. Our findings align with previous studies, indicating that the classical ribbed design of WABs functions more as a sound-absorbing structure rather than a true black hole with wave-focusing capability. In light of this finding, we formulated two separate topology optimization problems—one considering viscothermal losses at the boundaries and the other neglecting these losses. We employed a material distribution approach to solve the optimization problems. The results demonstrate the need to consider boundary losses in the optimization process. The optimized device in the lossy case outperformed the lossless one, exhibiting superior wave-focusing performance. In comparison to the classical ribbed design, both optimized devices showed improvements in achieving two key characteristics of acoustic black holes: increasing wave amplitude and gradual increase in the local effective axial wavenumber towards the end of the profile. Among the optimized devices, the lossy-case design consistently exhibited the best performance. Its ability to minimize both reflections and viscothermal losses in the targeted frequency range resulted in a significant wave-focusing effect. Based on these numerical results, we present the optimized device in the lossy case as the first promising conceptual design for a wave-focusing acoustic black hole in waveguides. Future work can explore further optimization strategies, analytical and numerical methods, and experimental validations to enhance the performance and practical implementation of WABs.

ACKNOWLEDGMENTS

This work was supported by the Swedish strategic research programme eSSENCE and the Swedish Research Council (Grant Nos. 2018-03546 and 2022-03783).

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

DATA AVAILABILITY

The data that support the findings of this study are available within the article.

APPENDIX

The method of MMA requires derivatives of each component of the objective function in optimization problem . Computing derivatives of the quadratic penalty term is straightforward. Also, for the power of the reflected wave, we can follow the adjoint computation of the derivatives as described by Appendix C in Wadbro (2014). The challenge is to compute derivatives for the term that represents viscothermal boundary power losses in primary objective function , that is, derivatives of the loss function

ℓ (α_{h}) = \int_{Θ} {[[α_{h}]]}^{2} (\frac{δ_{V}}{2 k} r | \nabla_{T} p_{h} |^{2} + δ_{T} k \frac{(γ - 1)}{2} r | p_{h} |^{2}) d Γ,

(A1)

where p_h is the solution to the state Eq. associated with the material indicator function

α_{h} \in A_{h}

and

A_{h} : Ω \to [ϵ, 1]

is the set of elementwise constant functions. Consider another, arbitrary material indicator function

{\hat{α}}_{h} \in A_{h}

⁠, define the perturbation field

α_{h}^{'} = {\hat{α}}_{h} - α_{h}

⁠, and consider the family of indicator functions

α_{h}^{t} = α_{h}^{t} = α_{h} + t α_{h}^{'}

parametrized by

t \in [0, 1]

⁠. Moreover, let

p_{h}^{t}

be the solution to state Eq. when α_h is replaced with

α_{h}^{t}

⁠. Computing the directional derivative of

ℓ

at α_h in the direction

α_{h}^{'}

⁠, we obtain

\begin{matrix} d ℓ (α_{h}; α_{h}^{'}) = lim_{t \to 0^{+}} \frac{1}{t} (ℓ (α_{h}^{t}) - ℓ (α_{h})) \\ = 2 \int_{Θ} [[α_{h}]] [[α_{h}^{'}]] (\frac{δ_{V}}{2 k} r | \nabla_{T} p_{h} |^{2} + δ_{T} k \frac{(γ - 1)}{2} r | p_{h} |^{2}) d Γ \\ + Re {\int_{Θ} {[[α_{h}]]}^{2} (\frac{δ_{V}}{k} r \nabla_{T} p_{h}^{'} \cdot \nabla_{T} {\bar{p}}_{h} + δ_{T} k (γ - 1) r p_{h}^{'} {\bar{p}}_{h}) d Γ}, \end{matrix}

(A2)

where

p_{h}^{'} = lim_{t \to 0^{+}} \frac{1}{t} (p_{h}^{t} - p_{h}) .

(A3)

The contribution to the derivative from the first term on the right side of Eq. can directly be computed for a given design α_h and solution p_h by appropriate choices of

α_{h}^{'}

⁠. To compute the contribution from the second term, however, we will employ an adjoint approach as described below.

Differentiating state problem with respect to the perturbation

α_{h}^{'}

of the material indicator function, we obtain that for each test function q_h,

\begin{matrix} \int_{Ω} α_{h}^{'} r \nabla q_{h} \cdot \nabla p_{h} d Ω - k^{2} \int_{Ω} α_{h}^{'} r q_{h} p_{h} d Ω + δ_{V} (i - 1) \int_{Θ} [[α_{h}]] [[α_{h}^{'}]] r \nabla_{T} q_{h} \cdot \nabla_{T} p_{h} d Γ \\ + δ_{T} k^{2} (i - 1) (γ - 1) \int_{Θ} [[α_{h}]] [[α_{h}^{'}]] r q_{h} p_{h} d Γ + \int_{Ω} α_{h} r \nabla q_{h} \cdot \nabla p_{h}^{'} d Ω - k^{2} \int_{Ω} α_{h} r q_{h} p_{h}^{'} d Ω + i k \int_{Γ^{in}} r q_{h} p_{h}^{'} d Γ + i k \int_{Γ^{L}} r q_{h} Y p_{h}^{'} d Γ \\ + δ_{V} \frac{i - 1}{2} \int_{Θ} {[[α_{h}]]}^{2} r \nabla_{T} q_{h} \cdot \nabla_{T} p_{h}^{'} d Γ + δ_{T} k^{2} \frac{(i - 1) (γ - 1)}{2} \int_{Θ} {[[α_{h}]]}^{2} r q_{h} p_{h}^{'} d Γ = 0. \end{matrix}

(A4)

Now let z_h be the solution to the adjoint equation

\begin{matrix} \int_{Ω} α_{h} r \nabla w_{h} \cdot \nabla z_{h} d Ω - k^{2} \int_{Ω} α_{h} r w_{h} z_{h} d Ω + i k \int_{Γ^{in}} r w_{h} z_{h} d Γ + i k \int_{Γ^{L}} r w_{h} Y z_{h} d Γ \\ + δ_{V} \frac{i - 1}{2} \int_{Θ} {[[α_{h}]]}^{2} r \nabla_{T} w_{h} \cdot \nabla_{T} z_{h} d Γ + δ_{T} \frac{i - 1}{2} \int_{Θ} {[[α_{h}]]}^{2} r w_{h} z_{h} d Γ \\ = \int_{Θ} {[[α_{h}]]}^{2} (\frac{δ_{V}}{k} r \nabla_{T} w_{h} \cdot \nabla_{T} {\bar{p}}_{h} + δ_{T} k (γ - 1) r w_{h} {\bar{p}}_{h}) d Γ \forall w_{h} . \end{matrix}

(A5)

Note that the left side of Eq. is the same as for state Eq. , so the matrix factorization for the state equation can be reused when solving adjoint Eq. . Evaluating adjoint Eq. for

w_{h} = p_{h}^{'}

⁠, we obtain

\begin{matrix} \int_{Ω} α_{h} r \nabla p_{h}^{'} \cdot \nabla z_{h} d Ω - k^{2} \int_{Ω} α_{h} r p_{h}^{'} z_{h} d Ω + i k \int_{Γ^{in}} r p_{h}^{'} z_{h} d Γ + i k \int_{Γ^{L}} r p_{h}^{'} Y z_{h} d Γ \\ + δ_{V} \frac{i - 1}{2} \int_{Θ} {[[α_{h}]]}^{2} r \nabla_{T} p_{h}^{'} \cdot \nabla_{T} z_{h} d Γ + δ_{T} \frac{i - 1}{2} \int_{Θ} {[[α_{h}]]}^{2} r p_{h}^{'} z_{h} d Γ \\ = \int_{Θ} {[[α_{h}]]}^{2} (\frac{δ_{V}}{k} r \nabla_{T} p_{h}^{'} \cdot \nabla_{T} {\bar{p}}_{h} + δ_{T} k (γ - 1) r p_{h}^{'} {\bar{p}}_{h}) d Γ . \end{matrix}

(A6)

Choosing q_h = z_h as the test function in Eq. and substituting the terms containing

p_{h}^{'}

using Eq. , we arrive at

\begin{matrix} \int_{Ω} α_{h}^{'} r \nabla z_{h} \cdot \nabla p_{h} d Ω - k^{2} \int_{Ω} α_{h}^{'} r z_{h} p_{h} d Ω + δ_{V} (i - 1) \int_{Θ} [[α_{h}]] [[α_{h}^{'}]] r \nabla_{T} z_{h} \cdot \nabla_{T} p_{h} d Γ \\ + δ_{T} k^{2} (i - 1) (γ - 1) \int_{Θ} [[α_{h}]] [[α_{h}^{'}]] r z_{h} p_{h} d Γ \\ + \int_{Θ} {[[α_{h}]]}^{2} (\frac{δ_{V}}{k} r \nabla_{T} p_{h}^{'} \cdot \nabla_{T} {\bar{p}}_{h} + δ_{T} k (γ - 1) r p_{h}^{'} {\bar{p}}_{h}) d Γ = 0. \end{matrix}

(A7)

Note that the last term on the left side of equality is exactly the second integral of directional derivative . Thus, combining Eqs. and , we find that

\begin{matrix} d ℓ (α_{h}; α_{h}^{'}) = 2 \int_{Θ} [[α_{h}]] [[α_{h}^{'}]] (\frac{δ_{V}}{2 k} r | \nabla_{T} p_{h} |^{2} + δ_{T} k \frac{(γ - 1)}{2} r | p_{h} |^{2}) d Γ \\ - Re {\int_{Ω} α_{h}^{'} r \nabla z_{h} \cdot \nabla p_{h} d Ω - k^{2} \int_{Ω} α_{h}^{'} r z_{h} p_{h} d Ω + δ_{V} (i - 1) \int_{Θ} [[α_{h}]] [[α_{h}^{'}]] r \nabla_{T} z_{h} \cdot \nabla_{T} p_{h} d Γ \\ + δ_{T} k^{2} (i - 1) (γ - 1) \int_{Θ} [[α_{h}]] [[α_{h}^{'}]] r z_{h} p_{h} d Γ} . \end{matrix}

(A8)

Therefore, to compute the derivatives, we follow the subsequent procedure. For a given design α_h, first, we solve the state Eq. to obtain the pressure field p_h. Next, we use this solution to solve the adjoint Eq. and compute the adjoint variable z_h. Once we have the adjoint variable, we can compute the derivatives of the loss function by evaluating Eq. for a sequence of appropriate choices of

α_{h}^{'}

⁠.

REFERENCES

1.

Andersen

,

P. R.

,

Henríquez

,

V. C.

, and

Aage

,

N.

(

2023

). “

On the validity of numerical models for viscothermal losses in structural optimization for micro-acoustics

,”

J. Sound Vib.

547

,

117455

.

2.

Bendsøe

,

M. P.

, and

Sigmund

,

O.

(

2003

).

Topology Optimization. Theory, Methods, and Applications

(

Springer Berlin

,

Heidelberg

).

https://doi.org/10.1016/j.jcp.2018.06.005

3.

Berggren

,

M.

,

Bernland

,

A.

, and

Noreland

,

D.

(

2018

). “

Acoustic boundary layers as boundary conditions

,”

J. Comput. Phys.

371

,

633

–

650

.

https://doi.org/10.1121/10.0018195

4.

Berggren

,

M.

,

Mousavi

,

A.

,

Hägg

,

L.

, and

Wadbro

,

E.

(

2023

). “

Topology optimization of wave-focusing waveguide acoustic black holes

,”

J. Acoust. Soc. Am.

153

,

A70

.

https://doi.org/10.1121/10.0002973

5.

Billard

,

R.

,

Tissot

,

G.

,

Gabard

,

G.

, and

Versaevel

,

M.

(

2021

). “

Numerical simulations of perforated plate liners: Analysis of the visco-thermal dissipation mechanisms

,”

J. Acoust. Soc. Am.

149

(

1

),

16

–

27

.

https://doi.org/10.1007/s00158-020-02832-9

PubMed

6.

Bokhari

,

A. H.

,

Mousavi

,

A.

,

Niu

,

B.

, and

Wadbro

,

E.

(

2021

). “

Topology optimization of an acoustic diode?

,”

Struct. Multidiscipl. Optim.

63

(

6

),

2739

–

2749

.

https://doi.org/10.1016/j.jsv.2023.117825

7.

Bravo

,

T.

, and

Maury

,

C.

(

2023

). “

Broadband sound attenuation and absorption by duct silencers based on the acoustic black hole effect: Simulations and experiments

,”

J. Sound Vib.

561

,

117825

.

https://doi.org/10.1016/j.wavemoti.2022.103039

17.

Červenka

,

M.

, and

Bednařík

,

M.

(

2022

). “

On the role of resonance and thermoviscous losses in an implementation of “acoustic black hole” for sound absorption in air

,”

Wave Motion

114

,

103039

.

https://doi.org/10.1016/j.ijmecsci.2022.107921

8.

Chen

,

X.

,

Zhao

,

J.

,

Deng

,

J.

,

Jing

,

Y.

,

Pu

,

H.

, and

Luo

,

J.

(

2023

). “

Low-frequency enhancement of acoustic black holes via negative stiffness supporting

,”

Int. J. Mech. Sci.

241

,

107921

.

https://doi.org/10.1016/j.compstruct.2022.116434

9.

Chua

,

J. W.

,

Li

,

X.

,

Yu

,

X.

, and

Zhai

,

W.

(

2023

). “

Novel slow-sound lattice absorbers based on the sonic black hole

,”

Compos. Struct.

304

,

116434

.

https://doi.org/10.1016/j.tws.2022.110459

10.

Deng

,

J.

,

Gao

,

N.

, and

Chen

,

X.

(

2023

). “

Ultrawide attenuation bands in gradient metabeams with acoustic black hole pillars

,”

Thin-Walled Struct.

184

,

110459

.

https://doi.org/10.1016/j.jsv.2020.115790

11.

Deng

,

J.

,

Guasch

,

O.

,

Zheng

,

L.

,

Song

,

T.

, and

Cao

,

Y.

(

2021a

). “

Semi-analytical model of an acoustic black hole piezoelectric bimorph cantilever for energy harvesting

,”

J. Sound Vib.

494

,

115790

.

https://doi.org/10.1016/j.apacoust.2020.107744

12.

Deng

,

J.

,

Zheng

,

L.

, and

Guasch

,

O.

(

2021b

). “

Elliptical acoustic black holes for flexural wave lensing in plates

,”

Appl. Acoust.

174

,

107744

.

https://doi.org/10.1016/j.ymssp.2018.08.053

13.

Deng

,

J.

,

Zheng

,

L.

,

Zeng

,

P.

,

Zuo

,

Y.

, and

Guasch

,

O.

(

2019

). “

Passive constrained viscoelastic layers to improve the efficiency of truncated acoustic black holes in beams

,”

Mech. Syst. Signal Process.

118

,

461

–

476

.

https://doi.org/10.1016/j.jsv.2008.03.042

14.

Dühring

,

M. B.

,

Jensen

,

J. S.

, and

Sigmund

,

O.

(

2008

). “

Acoustic design by topology optimization

,”

J. Sound Vib.

317

(

3-5

),

557

–

575

.

15.

El Ouahabi

,

A.

,

Krylov

,

V.

, and

O'Boy

,

D.

(

2015a

). “

Experimental investigation of the acoustic black hole for sound absorption in air

,” in

Proceedings of the 22nd International Congress on Sound and Vibration

,

Florence, Italy

.

16.

El Ouahabi

,

A.

,

Krylov

,

V.

, and

O'Boy

,

D.

(

2015b

). “

Investigation of the acoustic black hole termination for sound waves propagating in cylindrical waveguides

,” in

INTER-NOISE and NOISE-CON Congress and Conference Proceedings, Institute of Noise Control Engineering

,

San Francisco, CA

.

https://doi.org/10.1016/j.jsv.2017.02.007

18.

Guasch

,

O.

,

Arnela

,

M.

, and

Sánchez-Martín

,

P.

(

2017

). “

Transfer matrices to characterize linear and quadratic acoustic black holes in duct terminations

,”

J. Sound Vib.

395

,

65

–

79

.

https://doi.org/10.1016/j.apm.2019.09.039

19.

Guasch

,

O.

,

Sánchez-Martín

,

P.

, and

Ghilardi

,

D.

(

2020

). “

Application of the transfer matrix approximation for wave propagation in a metafluid representing an acoustic black hole duct termination

,”

Appl. Math. Model.

77

,

1881

–

1893

.

https://doi.org/10.1016/j.tws.2023.110821

20.

He

,

C.

,

Zhang

,

F.

,

Lim

,

K. M.

,

Jiang

,

J.

, and

Zhao

,

J.

(

2023

). “

Wave attenuation of a laminated acoustic black hole array in a load-bearing beam structure

,”

Thin-Walled Struct.

188

,

110821

.

21.

Ihlenburg

,

F.

(

1998

). “

Discretization methods for exterior Helmholtz problems

,” in

Finite Element Analysis of Acoustic Scattering

(

Springer

,

New York

), p.

63

.

https://doi.org/10.1016/j.ymssp.2019.06.034

22.

Ji

,

H.

,

Liang

,

Y.

,

Qiu

,

J.

,

Cheng

,

L.

, and

Wu

,

Y.

(

2019

). “

Enhancement of vibration based energy harvesting using compound acoustic black holes

,”

Mech. Syst. Signal Process.

132

,

441

–

456

.

https://doi.org/10.1016/j.jsv.2022.116981

23.

Karlos

,

A.

,

Hook

,

K.

, and

Cheer

,

J.

(

2022

). “

Enhanced absorption with multiple quadratically tapered elastic wedges of different lengths terminating a uniform beam

,”

J. Sound Vib.

531

,

116981

.

https://doi.org/10.1137/140981630

24.

Kasolis

,

F.

,

Wadbro

,

E.

, and

Berggren

,

M.

(

2015

). “

Analysis of fictitious domain approximations of hard scatterers

,”

SIAM J. Numer. Anal.

53

(

5

),

2347

–

2362

.

https://doi.org/10.1016/j.jsv.2003.05.010

25.

Krylov

,

V.

, and

Tilman

,

F.

(

2004

). “

Acoustic ‘black holes’ for flexural waves as effective vibration dampers

,”

J. Sound Vib.

274

(

3-5

),

605

–

619

.

https://doi.org/10.1016/j.ijsolstr.2021.111409

26.

Li

,

H.

,

Doaré

,

O.

,

Touzé

,

C.

,

Pelat

,

A.

, and

Gautier

,

F.

(

2022

). “

Energy harvesting efficiency of unimorph piezoelectric acoustic black hole cantilever shunted by resistive and inductive circuits

,”

Int. J. Solids Struct.

238

,

111409

.

https://doi.org/10.1016/j.jsv.2023.117781

27.

Li

,

S.

,

Xia

,

J.

,

Yu

,

X.

,

Zhang

,

X.

, and

Cheng

,

L.

(

2023

). “

A sonic black hole structure with perforated boundary for slow wave generation

,”

J. Sound Vib.

559

,

117781

.

https://doi.org/10.1016/j.jsv.2022.117047

28.

Liu

,

F.

,

Yang

,

Z.

,

Shi

,

P.

,

Shen

,

Y.

, and

Xu

,

Y.

(

2022

). “

Enhanced vibration absorption of plates with circular metasurface composed of lossy acoustic black hole subunits

,”

J. Sound Vib.

536

,

117047

.

https://doi.org/10.1016/j.ymssp.2022.109991

29.

Liu

,

P.

,

Huang

,

H.

,

Wang

,

X.

,

Tang

,

Q.

,

Qi

,

X.

,

Su

,

S.

,

Xiang

,

Z.

, and

Hu

,

J.

(

2023

). “

Acoustic black hole profiles for high-performance ultrasonic tweezers

,”

Mech. Syst. Signal Process.

188

,

109991

.

https://doi.org/10.1088/1361-665X/ab1f11

30.

Maugan

,

F.

,

Chesne

,

S.

,

Monteil

,

M.

,

Collet

,

M.

, and

Yi

,

K.

(

2019

). “

Enhancement of energy harvesting using acoustical-black-hole-inspired wave traps

,”

Smart Mater. Struct.

28

(

7

),

075015

.

https://doi.org/10.1016/j.jsv.2022.117138

31.

Mi

,

Y.

,

Cheng

,

L.

,

Zhai

,

W.

, and

Yu

,

X.

(

2022

). “

Broadband low-frequency sound attenuation in duct with embedded periodic sonic black holes

,”

J. Sound Vib.

536

,

11738

.

https://doi.org/10.1063/5.0042514

32.

Mi

,

Y.

,

Zhai

,

W.

,

Cheng

,

L.

,

Xi

,

C.

, and

Yu

,

X.

(

2021

). “

Wave trapping by acoustic black hole: Simultaneous reduction of sound reflection and transmission

,”

Appl. Phys. Lett.

118

(

11

),

114101

.

33.

Mironov

,

M.

(

1988

). “

Propagation of a flexural wave in a plate whose thickness decreases smoothly to zero in a finite interval

,”

Sov. Phys. Acoust.

34

,

318

–

319

.

https://doi.org/10.1016/j.jsv.2020.115223

34.

Mironov

,

M.

, and

Pislyakov

,

V.

(

2020

). “

One-dimensional sonic black holes: Exact analytical solution and experiments

,”

J. Sound Vib.

473

,

115223

.

https://doi.org/10.1134/1.1478121

35.

Mironov

,

M. A.

, and

Pislyakov

,

V. V.

(

2002

). “

One-dimensional acoustic waves in retarding structures with propagation velocity tending to zero

,”

Acoust. Phys.

48

(

3

),

347

–

352

.

https://doi.org/10.1121/10.0004664

36.

Mousavi

,

A.

,

Berggren

,

M.

, and

Wadbro

,

E.

(

2021

). “

On the acoustic black-hole effect in waveguides

,”

J. Acoust. Soc. Am.

149

(

4

),

A108

.

https://doi.org/10.1121/10.0011788

37.

Mousavi

,

A.

,

Berggren

,

M.

, and

Wadbro

,

E.

(

2022

). “

How the waveguide acoustic black hole works: A study of possible damping mechanisms

,”

J. Acoust. Soc. Am.

151

(

6

),

4279

–

4290

.

https://doi.org/10.1016/j.matdes.2023.112302

PubMed

38.

Mousavi

,

A.

,

Berggren

,

M.

, and

Wadbro

,

E.

(

2023

). “

Extending material distribution topology optimization to boundary-effect-dominated problems with applications in viscothermal acoustics

,”

Mater. Des.

234

,

112302

.

https://doi.org/10.1016/j.jsv.2019.04.029

39.

Ouisse

,

M.

,

Renault

,

D.

,

Butaud

,

P.

, and

Sadoulet-Reboul

,

E.

(

2019

). “

Damping control for improvement of acoustic black hole effect

,”

J. Sound Vib.

454

,

63

–

72

.

https://doi.org/10.1016/j.jsv.2020.115316

40.

Pelat

,

A.

,

Gautier

,

F.

,

Conlon

,

S. C.

, and

Semperlotti

,

F.

(

2020

). “

The acoustic black hole: A review of theory and applications

,”

J. Sound Vib.

476

,

15316

.

https://doi.org/10.1016/j.jsv.2022.117370

41.

Sheng

,

H.

,

He

,

M.-X.

, and

Ding

,

Q.

(

2023

). “

Vibration suppression by mistuning acoustic black hole dynamic vibration absorbers

,”

J. Sound Vib.

542

,

117370

.

https://doi.org/10.1002/nme.1620240207

42.

Svanberg

,

K.

(

1987

). “

The method of moving asymptotes-a new method for structural optimization

,”

Int. J. Numer. Meth. Eng.

24

(

2

),

359

–

373

.

https://doi.org/10.1016/j.ymssp.2019.106257

43.

Tang

,

L.

, and

Cheng

,

L.

(

2019

). “

Periodic plates with tunneled Acoustic-Black-Holes for directional band gap generation

,”

Mech. Syst. Signal Process.

133

,

106257

.

https://doi.org/10.1016/j.jsv.2022.117377

44.

Umnova

,

O.

,

Brooke

,

D.

,

Leclaire

,

P.

, and

Dupont

,

T.

(

2023

). “

Multiple resonances in lossy acoustic black holes - theory and experiment

,”

J. Sound Vib.

543

,

117377

.

https://doi.org/10.1007/s00158-014-1058-2

45.

Wadbro

,

E.

(

2014

). “

Analysis and design of acoustic transition sections for impedance matching and mode conversion

,”

Struct. Multidiscipl. Optim.

50

,

395

–

408

.

https://doi.org/10.1016/j.cma.2006.05.005

46.

Wadbro

,

E.

, and

Berggren

,

M.

(

2006

). “

Topology optimization of an acoustic horn

,”

Comput. Methods Appl. Mech. Eng.

196

(

1-3

),

420

–

436

.

https://doi.org/10.1016/j.jsv.2021.116401

47.

Xiaoqi

,

Z.

, and

Cheng

,

L.

(

2021

). “

Broadband and low frequency sound absorption by Sonic black holes with Micro-perforated boundaries

,”

J. Sound Vib.

512

,

116401

.

https://doi.org/10.1016/j.apacoust.2022.109078

48.

Zhao

,

J.

,

Huang

,

Y.

,

Yuan

,

W.

,

Zhang

,

J.

,

Song

,

C.

,

Zhang

,

X.

, and

Pan

,

Y.

(

2022

). “

Broadband acoustic black hole for wave focusing and weak signal sensing

,”

Appl. Acoust.

200

,

109078

.

https://doi.org/10.1016/j.ijmecsci.2022.107907

49.

Zhao

,

L.

,

Bi

,

C.

, and

Yu

,

M.

(

2023

). “

Structural lens for broadband triple focusing and three-beam splitting of flexural waves

,”

Int. J. Mech. Sci.

240

,

107907

.

https://doi.org/10.1088/0964-1726/23/6/065021

50.

Zhao

,

L.

,

Conlon

,

S. C.

, and

Semperlotti

,

F.

(

2014

). “

Broadband energy harvesting using acoustic black hole structural tailoring

,”

Smart Mater. Struct.

23

(

6

),

065104

.

https://doi.org/10.1063/1.4998524

51.

Zhu

,

H.

, and

Semperlotti

,

F.

(

2017

). “

Two-dimensional structure-embedded acoustic lenses based on periodic acoustic black holes

,”

J. Appl. Phys.

122

(

6

),

065104

.