# Optimal profile design for acoustic black holes using Timoshenko beam theory

Kasper S. Sørensen<EMAIL_ADDRESS>Department of Mathematical Sciences,
Aalborg University
Skjernvej 4A, 9220 Aalborg Ø, Denmark.  ,  Horia D. Cornean
<EMAIL_ADDRESS>Department of Mathematical Sciences, Aalborg University
Skjernvej 4A, 9220 Aalborg Ø, Denmark.  and  Sergey Sorokin<EMAIL_ADDRESS>Department of Materials and Production, Aalborg University
Fibigerstræde 16, 9220 Aalborg Ø, Denmark.

###### Abstract.

We revisit the problem of constructing one-dimensional acoustic black holes.
Instead of considering the Euler-Bernoulli beam theory, we use Timoshenko’s
approach instead, which is known to be more realistic at higher frequencies.
Our goal is to minimize the reflection coefficient under a constraint imposed
on the normalized wave number variation. We use the calculus of variations in
order to derive the corresponding Euler-Lagrange equation analytically and
then use numerical methods to solve this equation in order to find the optimal
height profile for different frequencies. We then compare these profiles to
the corresponding ones previously found using the Euler-Bernoulli beam theory
and see that in the lower range of the dimensionless frequency $\Omega$
(defined using the largest height of the plate), the optimal profiles almost
coincide, as expected. For higher such frequencies, even for values where
Euler-Bernoulli theory should still be marginally valid, the profiles
predicted using Euler-Bernoulli differ substantially from the correct ones
predicted by Timoshenko theory. One explanation for this phenomenon is that
unlike in the constant height case, in our setting the wave numbers also
depend on the ratio between the smallest and the largest heights.

##  1\. Introduction

The study of acoustic black holes originates from the seminal work of Mironov
[10]. He showed, using analytic methods, that if one can construct a plate
with an ideal wedge (i.e. the thickness of the wedge goes smoothly to zero in
a finite interval (see Fig. 1)), then the velocity of the flexural wave moving
towards the tip would go to zero and thus never reach the end of the rod. This
is equivalent to saying that no reflection of the flexural wave would occur
and hence an acoustic black hole has been created.

Unfortunately, it is not possible to create an ideal wedge in practice, thus
the current research concerning acoustic black holes is focused around
minimizing the reflection of the wave when it hits the end of the wedge.
Acoustic black holes for other geometries and higher dimensions have also been
considered, see eg. [8] and [4] for different geometries and [6] for the
2-dimensional case. Several different directions have been studied in
connection to this minimization problem, e.g. profile optimization by adding a
thin dampening layer on the wedge. To the best of the authors knowledge, the
study of acoustic black holes has, up until now, only been considered using
the Euler-Bernoulli beam theory. We refer the reader to the recent survey
paper [11] and the references therein for a thorough overview of the state of-
the-art in this field.

In our paper, inspired by some recent work [16], we apply the mathematical
theory of calculus of variations in order to determine an optimal wedge
profile which minimizes the reflection of flexural waves. For the first time
regarding acoustic black holes, we consider waves using the Timoshenko beam
theory. Timoshenko’s approach can be considered as an extension of the one of
Euler-Bernoulli, in the sense that the their predictions almost coincide for
dimensionless frequencies (defined in (1.2) where $h_{1}$ is the largest
possible height of the plate) up to $\Omega\approx 0.3$, see e.g. [9] and
[14].

For dimensionless frequencies higher than this value, it is known that the
Euler-Bernoulli beam theory starts deviating from the true one, while the
Timoshenko approach is known (from the comparison with the exact solution of
the Rayleigh-Lamb problem [9],[1]) to hold true up to $\Omega\approx 3.5$ (see
eg. [15]). This is in part due to the fact the the Timoshenko beam theory
takes into account shear deformation along with rotational bending effects,
which are neglected in the Euler-Bernoulli approximation.

We stress that the above values for $\Omega$ are derived for a constant height
of the plate. Because we let the height of the plate to vary, we will see that
the wave numbers also depend on the ratio between the smallest and the largest
heights. This means that in our setting, $\Omega$ is no longer the only
parameter which decides the transition between the two regimes.

Due to this much higher validity range of the Timoshenko approach, we find it
natural to investigate the possibility of constructing acoustic black holes
within this more complete theory. Moreover, it is also quite relevant to
compare our results with the earlier ones based on Euler-Bernoulli, and to
find the frequency threshold above which the optimal profiles predicted by two
approaches start to substantially differ from each other.

The Timoshenko beam theory also predicts the appearance of a second wave. This
second wave is considered analytically in our paper, but we have not yet been
able to identify real materials where the second wave appears in an acoustic
frequency range. Thus the numerical investigations in this paper are only with
regard to the first wave. A numerical investigation of the second wave, along
with trying to find some material (or meta-material) where the second wave
appears in the acoustic frequency range, is very interesting but it will be
postponed for a future work.

###  1.1. Problem statement

The primary aim of this paper is to determine a height profile for a truncated
plate (see Fig. 2), which minimizes the reflection of a flexural wave, when
the wave propagation is considered within Timoshenko’s beam theory. The system
of partial differential equations describing the vibration of a Timoshenko
beam is given by (see [9] for a derivation)

$\displaystyle-\rho A\frac{\partial^{2}w}{\partial\tau^{2}}(x,\tau)+\kappa
GA\Big{(}\frac{\partial^{2}w}{\partial
x^{2}}(x,\tau)-\frac{\partial\psi}{\partial x}(x,\tau)\Big{)}+q(x,\tau)$
$\displaystyle=0,$ (1.1) $\displaystyle-\rho
I\frac{\partial^{2}\psi}{\partial\tau^{2}}(x,\tau)+EI\frac{\partial^{2}\psi}{\partial
x^{2}}(x,\tau)+\kappa GA\Big{(}\frac{\partial w}{\partial
x}(x,\tau)-\psi(x,\tau)\Big{)}$ $\displaystyle=0,$

where $\rho$ is the density of the material, $A=Bh_{d}(x)$ is the cross
section area where $B$ is the constant width, while $h_{d}(x)$ is the height,
$E=E_{0}(1-i\eta)$ is the complex elastic modulus with loss $\eta>0$,
$G=\tfrac{E}{2(1+\nu)}$ is the shear modulus with $\nu$ its Poisson ratio,
$I=\tfrac{B(h_{d}(x))^{3}}{12}$ is the second moment area, $\kappa$ is the
Timoshenko shear coefficient and $q$ is the distributed load.

0$x_{1}$0$h_{1}$$-h_{1}$
Figure 1. Non-truncated wedge $x_{0}$$x_{1}$$h_{0}$$h_{1}$$h_{d}(x)$
Figure 2. Truncated wedge

We will see in the next section that when looking for propagating solutions
with an angular frequency $\omega$ to the ”locally homogeneous” system (1.1)
(i.e. with $q=0$), the ”local” wave number $k_{d}$ (”$d$” from dimensional)
must depend on the profile through the height $h$ of the profile.

We also denote with $c=\sqrt{\tfrac{E}{\rho}}$ the P-wave propagation speed,
and with $\tilde{W}$ the amplitude of the lateral displacement $w(x,t)$. We
define some dimensionless variables as follows (see Fig. 3):

$\displaystyle h_{\ell}$ $\displaystyle=\frac{h_{0}}{h_{1}},\quad
t=\frac{x_{1}-x_{0}}{h_{1}},\quad 0\leq\xi\leq t,\quad\Omega=\frac{\omega
h_{1}}{c},\quad W=\frac{\tilde{W}}{h_{1}},$ (1.2) $\displaystyle
h(\xi)=\frac{h_{d}(x_{0}+\xi h_{1})}{h_{1}},\quad k(\xi)=k_{d}(x_{0}+\xi
h_{1})h_{1}.$

0$t$$h_{\ell}$$1$$h(\xi)$

Figure 3. Dimensionless, truncated wedge

In order to quantify the performance of an acoustic black hole, we use the
reflection coefficient (see eg. [5] or [10])

$\displaystyle
R=\exp\Big{(}-2\int_{0}^{t}\mathrm{Im}(k(\xi))\;\mathrm{d}\xi\Big{)},$ (1.3)

a formula whose validity is based on a slow variation assumption [7]:

$\displaystyle\Big{|}\frac{1}{k^{2}}\frac{\;\mathrm{d}k}{\;\mathrm{d}\xi}\Big{|}\ll
1.$ (1.4)

We refer to the left-hand side of (1.4) as the normalized wave number
variation. In practice, a normalized wave number variation less than $0.3$ is
often seen as reasonable [2].

To summarize, given $k(\xi)$ as a function of the profile (to be found in the
next section), our aim is to find a height function
$h\colon[0,t]\to[h_{\ell},1]$ (see Fig 3) which satisfies:

  1. (i)

$h\in C^{1}$, i.e. differentiable with a continuous derivative,

  2. (ii)

$h(0)=h_{\ell}$, $h(t)=1$, and $h$ is non-decreasing,

  3. (iii)

$h$ is a minimizer of (1.3) under the constraint (1.4); this statement will be
made mathematically more precise later.

Note, that $h$ being a minimizer of (1.3) is equivalent to $h$ being a
maximizer of

$\displaystyle\int_{0}^{t}\mathrm{Im}(k(\xi))\;\mathrm{d}\xi.$ (1.5)

###  1.2. Structure of the paper

The rest of this paper is as follows. In Section 2 we develop an analytic
framework in order to maximize (1.5) under some further conditions. This
framework is based on the calculus of variations and we formulate an Euler-
Lagrange equation which must be satisfied by a maximizer.

Section 3 is dedicated to the numerical implementation of the analytic results
from the first part. At the end, we summarize our results and formulate a
number of open questions which we believe to be worth considering in the near
future.

##  2\. The analytic results

Our first task is to derive an expression for the ”local” wave numbers,
reasoning like in the WKB approximation under a slow variation assumption.
Namely, we look for a solution of (1.1) when $q=0$ of the form

$\displaystyle w(x,t)=W(x)\;h_{1}\;e^{ik_{d}(x)x-i\omega t}\qquad\textup{ and
}\qquad\psi(x,t)=\Psi(x)e^{ik_{d}(x)x-i\omega t}.$

The zeroth-order component in the WKB-like expansion leads to the following
system of linear equations in $W$ and $\Psi$:

$\displaystyle\Big{(}-k^{2}+\frac{2(1+\nu)}{\kappa}\Omega^{2}\Big{)}W-ik\Psi=0,$
$\displaystyle\frac{6ik\kappa}{(1+\nu)(h(\xi))^{2}}W+\Big{(}-k^{2}+\Omega^{2}-\frac{6\kappa}{(1+\nu)(h(\xi))^{2}}\Big{)}\Psi=0.$

In order to get a non-trivial solution for the pair $W$ and $\Psi$, the
coefficient matrix of the above linear system must be singular, hence its
determinant must vanish. This gives the Timoshenko dispersion relation

$\displaystyle
k^{4}-\Big{(}1+\frac{2(1+\nu)}{\kappa}\Big{)}k^{2}\Omega^{2}-\frac{12}{(h(\xi))^{2}}\Omega^{2}+\frac{2(1+\nu)}{\kappa}\Omega^{4}=0.$
(2.1)

Next we want to find the solutions to this equation. To do so, we observe that
the left-hand side is an ordinary second order polynomial with $k^{2}$ as
variable. In our paper the complex square root is defined as follows:
$z=re^{i\phi}$ with $r>0$ and $0\leq\phi<2\pi$, then
$\sqrt{z}=\sqrt{r}\;e^{i\phi/2}$.

Even though there are four algebraic solutions, we are only interested in
those which remain bounded for $x\geq 0$. This demands that the imaginary part
of the physical solutions is non-negative. Thus the solution(s) are among the
expressions below:

$\displaystyle
k(\xi)=\pm\frac{\sqrt{(\kappa+2(1+\nu))\Omega^{2}\pm\sqrt{(\kappa-2(1+\nu))^{2}\Omega^{4}+\frac{48\kappa^{2}}{h(\xi)^{2}}\Omega^{2}}}}{\sqrt{2\kappa}},$
(2.2)

with the extra condition ${\rm Im}(k)\geq 0$. On the other hand, we are
neither interested in evanescent waves where the imaginary part is large. That
is why it is important to analyze how the imaginary part of these solutions
behaves as a function of $\Omega$. We can already see from the above formula
that, by our choice of scaling, the term
$\frac{48\kappa^{2}}{h(\xi)^{2}}\Omega^{2}$ can become very large when
$h(\xi)$ gets very small, hence the transition between Euler-Bernoulli and
Timoshenko is no longer uniquely determined by $\Omega$, but also by the ratio
$h_{\ell}$ between the smallest and largest values of the height.

###  2.1. Frequency ranges

We now make a short detour in order to study the behaviour of the wave numbers
in different frequency ranges. The dimensionless frequency $\Omega$ has -by
construction- an infinitesimal positive imaginary part, which can be traced
back to the negative imaginary part of the elastic modulus $E$; the same holds
true for $\Omega^{2}$ and $\Omega^{4}$. It is very important to warn the
reader that in our setting, the transition between Euler-Bernoulli and
Timoshenko theories is not solely determined by $\Omega$, but also by
$h_{\ell}$.

We start our analysis in the regime $|\Omega|\ll 1$, where Euler-Bernoulli’s
and Timoshenko’s beam theories should agree (see [9] and [14]). Let us confirm
this. In this frequency regime, we have $|\Omega|^{4}\ll|\Omega|^{2}$, hence
the inner square root in (2.2) ”almost” equals
$\frac{\sqrt{48}\kappa}{h(\xi)}\Omega$. This linear term is also much larger
in absolute value than $|\Omega|^{2}$, hence:

$\displaystyle
k(\xi)\approx\pm\frac{\sqrt[4]{12}}{\sqrt{h(\xi)}}\sqrt{\pm\Omega}.$

The solution with $\sqrt{-\Omega}$ is evanescent, hence the only solution with
a non-negative small imaginary part is the one with two $+$ signs, which
coincides with the Euler-Bernoulli solution.

In other words, we just showed that the solution

$\displaystyle
k_{+}(\xi)=\frac{\sqrt{(\kappa+2(1+\nu))\Omega^{2}+\sqrt{(\kappa-2(1+\nu))^{2}\Omega^{4}+\frac{48\kappa^{2}}{(h(\xi))^{2}}\Omega^{2}}}}{\sqrt{2\kappa}}$
(2.3)

is the one which recovers the Euler-Bernoulli wavenumber for small
frequencies. When the frequency grows, the term with $\Omega^{4}$ becomes
dominant and there exists a critical frequency $\Omega_{\mathrm{crit}}$, such
that when $|\Omega|>\Omega_{\mathrm{crit}}$, the wave behaves like a Rayleigh
wave, i.e.

$k_{+}\approx\sqrt{\frac{|\kappa-2(1+\nu)|}{2\kappa}}\;\Omega.$

This means that for very high frequencies, the shape profile will have no
impact on this wave. Note, that $\Omega_{\mathrm{crit}}$ can be chosen to be
the solution to the equation
$(\kappa-2(1+\nu))^{2}\Omega_{\mathrm{crit}}^{2}=\frac{48\kappa^{2}}{h_{\ell}^{2}}$.

On the other hand, when $|\Omega|$ grows, there exists a second wave number
which obeys our two conditions, namely:

$\displaystyle
k_{-}(\xi)=\frac{\sqrt{(\kappa+2(1+\nu))\Omega^{2}-\sqrt{(\kappa-2(1+\nu))^{2}\Omega^{4}+\frac{48\kappa^{2}}{(h(\xi))^{2}}\Omega^{2}}}}{\sqrt{2\kappa}}.$

which appears when $|\Omega|>\Omega_{\mathrm{cut-on}}$, where
$\Omega_{\mathrm{cut-on}}$ solves
$8\kappa(1+\nu)\Omega^{2}=48\kappa^{2}h_{\ell}^{-2}$. In this frequency range
the behaviour of the second wave number is given as

$\displaystyle
k_{-}(\xi)\approx|\kappa-2(1+\nu)|^{-1/2}\Big{(}8\kappa(1+\nu)\Omega^{2}-\frac{48\kappa^{2}}{(h(\xi))^{2}}\Big{)}^{1/2}.$
(2.4)

It should be noted that in general, it is expected that the validity of the
Timoshenko beam theory for the second wave only holds “close” to
$\Omega_{\mathrm{cut-on}}$.

###  2.2. The imaginary part of the wave number

We now return to solving the minimization problem. To get an expression for
the imaginary part of the wave numbers, we make a first order approximation of
(2.2) with respect to $\eta$. This leads to

$\displaystyle
k(\xi)\approx\pm\Big{(}\frac{\omega^{2}h_{1}^{2}}{2c_{0}^{2}\kappa}\Big{)}^{1/2}\sqrt{\kappa^{+}(1+i\eta)\pm\sqrt{(\kappa^{-})^{2}+\tilde{h}_{\xi}+i\eta\Big{(}2(\kappa^{-})^{2}+\tilde{h}_{\xi}\Big{)}}},$
(2.5)

where $c_{0}^{2}=\tfrac{E_{0}}{\rho}$, $\kappa^{\pm}=\kappa\pm 2(1+\nu)$ and
$\tilde{h}_{\xi}=\tfrac{48\kappa^{2}c_{0}^{2}}{(h(\xi))^{2}\omega^{2}h_{1}^{2}}$.
It has the form $\sqrt{a+bi\pm\sqrt{c+di}}$, where $|b|,|d|\ll 1$ and $c>0$.
By a first order Taylor approximation we can write this as

$\displaystyle\sqrt{a+bi\pm\sqrt{c+di}}\approx\sqrt{a\pm\sqrt{c}}+i\frac{b\pm\frac{d}{2\sqrt{c}}}{2\sqrt{a\pm\sqrt{c}}},$

when $\kappa^{+}\pm\sqrt{(\kappa^{-})^{2}+\tilde{h}_{\xi}}>0$. Thus under this
condition, we can approximate the non-negative imaginary part of the two
solutions by

$\displaystyle\mathrm{Im}(k_{\pm})\approx\eta\Big{(}\frac{\omega^{2}h_{1}^{2}}{8c_{0}^{2}\kappa}\Big{)}^{1/2}\frac{\kappa^{+}\pm\frac{2(\kappa^{-})^{2}+\tilde{h}_{\xi}}{2\sqrt{(\kappa^{-})^{2}+\tilde{h}_{\xi}}}}{\sqrt{\kappa^{+}\pm\sqrt{(\kappa^{-})^{2}+\tilde{h}_{\xi}}}}.$

###  2.3. Dampening layer

Before we continue with the method of maximizing this integral, we make a
comment about attaching an absorbing layer to the wedge surfaces. Here we
follow the ideas from [5] and [7].

By attaching a visco-elastic layer of constant thickness $\delta$ to one of
the sides of a plate of constant thickness $h$, we get the following
additional loss factor

$\displaystyle\Xi=\frac{\tilde{\eta}}{1+\big{(}\alpha_{2}\beta_{2}(\alpha_{2}^{2}+12\alpha_{21}^{2})\big{)}^{-1}},$

where $\tilde{\eta}$ is the loss factor of the visco-elastic layer,
$\alpha_{2}=\delta/h$, $\beta_{2}=E_{1}/E_{0}$, where $E_{0}$ is Young’s
modulus for the plate and $E_{1}$ is Young’s modulus for the visco-elastic
layer and $\alpha_{21}=(1+\alpha_{2})/2$. This formula has been derived under
the assumption $(\delta/h)(E_{1}/E_{0})\ll 1$.

If we instead cover both sides of the wedge with a a visco-elastic layer,
which is thin compared to the thickness of the plate (which here is allowed to
be non-constant), i.e. $\delta/h(\xi)\ll h_{1}$, one can derive a formula for
the wave numbers, where $\eta$ is simply replaced by

$\displaystyle\eta+\frac{3}{2}\frac{\delta}{h(\xi)h_{1}}\frac{E_{1}}{E_{0}}\tilde{\eta}.$

Following the same first order expansion as above, the imaginary part of the
physical wave numbers are:

$\displaystyle\mathrm{Im}(k_{\pm})\approx\Big{(}\frac{\omega^{2}h_{1}^{2}}{8c_{0}^{2}\kappa}\Big{)}^{1/2}\frac{\kappa^{+}\pm\frac{2(\kappa^{-})^{2}+\tilde{h}_{\xi}}{2\sqrt{(\kappa^{-})^{2}+\tilde{h}_{\xi}}}}{\sqrt{\kappa^{+}\pm\sqrt{(\kappa^{-})^{2}+\tilde{h}_{\xi}}}}\Big{[}\eta+\frac{3}{2}\frac{\delta}{h(\xi)h_{1}}\frac{E_{1}}{E_{0}}\tilde{\eta}\Big{]},\quad\kappa^{+}\pm\sqrt{(\kappa^{-})^{2}+\tilde{h}_{\xi}}>0.$
(2.6)

An interesting open question here would be to consider how one could get a
better performing acoustic black hole by including the thickness $\delta$ of
the visco-elastic layer in the optimization problem by allowing it to vary as
a function of $\xi$, but still under the condition that $\delta/h(\xi)\ll
h_{1}$. We will not address this issue here.

###  2.4. Lagrange optimization inspired method

The problem of maximizing the integral (1.5) under a pointwise constraint like
in (1.4) is quite complicated, for which no general recipe is available.
Instead of a pointwise constraint, we will impose a weaker integral condition,
implemented by a penalty term, which guarantees at least that (1.4) cannot be
violated on a too large interval in $\xi$. Our method is related to the
Lagrange multiplier method from calculus of variation (see [3] and [12]).

More precisely, given some $n\geq 1$ and $\beta>0$, we will maximize the
functional

$\displaystyle
I(h)=\int_{0}^{t}\mathrm{Im}(k(\xi))\;\mathrm{d}\xi-\int_{0}^{t}\Big{(}\beta^{-1}\Big{|}\frac{1}{k^{2}}\frac{\;\mathrm{d}k}{\;\mathrm{d}\xi}\Big{|}\Big{)}^{2n}\;\mathrm{d}\xi,$
(2.7)

over $C^{1}$ dimensionless, non-negative and non-decreasing height functions
$h$ which obey $h(0)=h_{\ell}$ and $h(1)=1$. The maximizer of such a
functional should obey

$\displaystyle\Big{|}\frac{1}{k^{2}}\frac{\;\mathrm{d}k}{\;\mathrm{d}\xi}\Big{|}\leq\beta$

on a large piece of the interval $[0,t]$, larger as $n$ grows, otherwise the
penalty term would decrease the value of the functional. This penalty is
similar to the one introduced in topology optimization [13].

###  2.5. Euler-Lagrange equation

A great advantage coming from working with (2.7) is that we can rewrite the
integrand in terms of $h(\xi)$ and $h^{\prime}(\xi)$ and identify an Euler-
Lagrange equation for $h$. The integrand does not explicitly depend on $\xi$,
which makes that the ”energy” is conserved and we can reduce the problem to a
first order equation. Using the expression from (2.5) of $k(\xi)$ where we
neglect the imaginary part after differentiation, by elementary calculations
it follows that

$\displaystyle\Bigg{|}\frac{1}{k^{2}}\frac{\mathrm{d}k}{\mathrm{d}\xi}\Bigg{|}^{2n}$
$\displaystyle\approx
h^{\prime}(\xi)^{2n}\frac{1}{4^{2n}}\Big{(}\frac{2c_{o}^{2}\kappa}{\omega^{2}h_{1}^{2}}\Big{)}^{n}\Bigg{|}\kappa^{+}\pm\sqrt{(\kappa^{-})^{2}+\tilde{h}_{\xi}}\Bigg{|}^{-3n}\Big{(}(\kappa^{-})^{2}+\tilde{h}_{\xi}\Big{)}^{-n}\Bigg{(}\frac{96\kappa^{2}c_{0}^{2}}{(h(\xi))^{3}\omega^{2}h_{1}^{2}}\Bigg{)}^{2n}.$

To ease notation in the following let
$\tilde{h}(x)\coloneqq\frac{48\kappa^{2}c_{0}^{2}}{x^{2}\omega^{2}h_{1}^{2}}$.
We define the function
$\tilde{F}_{\pm}\colon(0,\infty)\times\mathbb{R}\to\mathbb{R}$ by

$\displaystyle\tilde{F}_{\pm}(x,y)$
$\displaystyle\coloneqq\Big{(}\frac{\omega^{2}h_{1}^{2}}{8c_{0}^{2}\kappa}\Big{)}^{1/2}\frac{\kappa^{+}\pm\frac{2(\kappa^{-})^{2}+\tilde{h}(x)}{2\sqrt{(\kappa^{-})^{2}+\tilde{h}(x)}}}{\sqrt{\kappa^{+}\pm\sqrt{(\kappa^{-})^{2}+\tilde{h}(x)}}}\Big{[}\eta+\frac{3}{2}\frac{\delta}{xh_{1}}\frac{E_{1}}{E_{0}}\tilde{\eta}\Big{]}$
$\displaystyle-y^{2n}\beta^{-2n}\frac{1}{4^{2n}}\Big{(}\frac{2c_{0}^{2}\kappa}{\omega^{2}h_{1}^{2}}\Big{)}^{n}\Bigg{|}\kappa^{+}\pm\sqrt{(\kappa^{-})^{2}+\tilde{h}(x)}\Bigg{|}^{-3n}\Big{(}(\kappa^{-})^{2}+\tilde{h}(x)\Big{)}^{-n}\Bigg{(}\frac{96\kappa^{2}c_{0}^{2}}{x^{3}\omega^{2}h_{1}^{2}}\Bigg{)}^{2n}.$

Our original functional can now be rewritten as

$\int_{0}^{t}\tilde{F}_{\pm}(h(\xi),h^{\prime}(\xi))\;\mathrm{d}\xi.$

Since the maximizer is invariant when the functional is multiplied by a
positive constant, we choose to divide our functional with all the constants
on the second term of $\tilde{F}$ and consider the function
$F_{\pm}\colon(0,\infty)\times\mathbb{R}\to\mathbb{R}$, given by

$\displaystyle F_{\pm}(x,y)$
$\displaystyle=b\frac{\kappa^{+}\pm\frac{2(\kappa^{-})^{2}+\tilde{h}(x)}{2\sqrt{(\kappa^{-})^{2}+\tilde{h}(x)}}}{\sqrt{\kappa^{+}\pm\sqrt{(\kappa^{-})^{2}+\tilde{h}(x)}}}\Big{[}\eta+\frac{3}{2}\frac{\delta}{xh_{1}}\frac{E_{1}}{E_{0}}\tilde{\eta}\Big{]}$
$\displaystyle\qquad\qquad-y^{2n}\Bigg{|}\kappa^{+}\pm\sqrt{(\kappa^{-})^{2}+\tilde{h}(x)}\Bigg{|}^{-3n}\Big{(}(\kappa^{-})^{2}+\tilde{h}(x)\Big{)}^{-n}x^{-6n},$

where
$b=2^{-(7n+3/2)}3^{-2n}\beta^{2n}\omega^{6n+1}h_{1}^{6n+1}(c_{0}^{2})^{-(3n+1/2)}\kappa^{-(5n+1/2)}$.
From now on we will separately investigate the existence of a maximizer for
the functionals

$\displaystyle I_{\pm}(h)=\int_{0}^{t}F_{\pm}(h(\xi),h^{\prime}(\xi))d\xi.$

Since $F_{\pm}$ do not depend explicitly on $\xi$, it follows by a standard
result from Calculus of Variation (see eg. [3]), that if a maximizer exists,
then it must obey the first order initial value problem given by:

$\displaystyle
F_{\pm}(h(\xi),h^{\prime}(\xi))-h^{\prime}(\xi)(\partial_{y}F_{\pm})(h(\xi),h^{\prime}(\xi))=a,\quad
h(t)=1,$

where $a$ is some constant which has to be chosen so that a $C^{1}$ non-
decreasing solution exists on $(0,t)$ and which obeys $h(0)=h_{\ell}$. We
have:

$\displaystyle\frac{\partial}{\partial
y}F_{\pm}(x,y)=-2ny^{2n-1}\Bigg{(}\kappa^{+}\pm\sqrt{(\kappa^{-})^{2}+\tilde{h}(x)}\Bigg{)}^{-3n}\Big{(}(\kappa^{-})^{2}+\tilde{h}(x)\Big{)}^{-n}x^{-6n}.$

Isolating $h^{\prime}(\xi)$ and asking that $h^{\prime}\geq 0$, we obtain a
separable differential equation:

$\displaystyle h^{\prime}(\xi)$
$\displaystyle=\Bigg{\\{}\Bigg{[}a-b\frac{\kappa^{+}\pm\frac{2(\kappa^{-})^{2}+\tilde{h}_{\xi}}{2\sqrt{(\kappa^{-})^{2}+\tilde{h}_{\xi}}}}{\sqrt{\kappa^{+}\pm\sqrt{(\kappa^{-})^{2}+\tilde{h}_{\xi}}}}\Big{[}\eta+\frac{3}{2}\frac{\delta}{h(\xi)h_{1}}\frac{E_{1}}{E_{0}}\tilde{\eta}\Big{]}\Bigg{]}$
$\displaystyle\quad\times\Bigg{[}(2n-1)\Bigg{|}\kappa^{+}\pm\sqrt{(\kappa^{-})^{2}+\tilde{h}_{\xi}}\Bigg{|}^{-3n}\Big{(}(\kappa^{-})^{2}+\tilde{h}_{\xi}\Big{)}^{-n}(h(\xi))^{-6n}\Bigg{]}^{-1}\Bigg{\\}}^{1/2n}.$
(2.8)

We have not been able to explicitly solve this differential equation and we
had to appeal to ”semi-analytical” methods coupled with numerical tools,
described in what follows. First, remember that
$\tilde{h}_{\xi}=\frac{48\kappa^{2}c_{0}^{2}}{h(\xi)^{2}\omega^{2}h_{1}^{2}}$.
Second, we notice that the right-hand side of (2.5) can be seen as a function
$g_{a}:(0,\infty)\mapsto(0,\infty)$ of $h(\xi)$, which also depends
parametrically on $a$. Hence (2.5) can be rewritten as

$h^{\prime}(\xi)=g_{a}(h(\xi)),\quad h(t)=1.$

Any local solution near $\xi=t$ must obey the functional equation

$t-\xi=\int_{h(\xi)}^{1}\frac{1}{g_{a}(y)}\;\mathrm{d}y.$ (2.9)

If we start by choosing a very large $a$, then $1/g_{a}$ stays bounded on the
whole interval $y\in[h_{\ell},1]$ and we may compute the map

$a\mapsto\Phi(a):=\int_{h_{\ell}}^{1}\frac{1}{g_{a}(y)}\;\mathrm{d}y$

where the right-hand side goes to zero when $a\to\infty$. This means that if
$a$ is large enough we have $t>\Phi(a)$, thus no such $a$ can insure the other
boundary condition $h(0)=h_{\ell}$.

Third, by reducing the value of $a$, $\Phi(a)$ increases continuously and
there can be at most one value $a_{*}>0$ for which $\frac{1}{g_{a_{*}}(y)}$ is
still integrable on $[h_{\ell},1]$ and also $\Phi(a_{*})=t$. With this value
of $a=a_{*}$ we go back to (2.9) and find the unique value of
$x\in[h_{\ell},1]$ for which $\int_{x}^{1}\frac{1}{g_{a_{*}}(y)}\;\mathrm{d}y$
equals $t-\xi$ for every $\xi\in[0,t]$.

For some high frequency values, another phenomenon can appear. It might happen
that by decreasing $a$ we arrive to a critical value $\tilde{a}_{*}$ such that
$\frac{1}{g_{a}(y)}$ is no longer well-defined/integrable for
$a<\tilde{a}_{*}$, while $\Phi(\tilde{a}_{*})$ is still less than $t$. In this
case, the optimal profile hits the value $h_{\ell}$ earlier than at $\xi=0$,
namely at the point $t-\Phi(\tilde{a}_{*})$, and then stays constant equal to
$h_{\ell}$ between this point and zero.

###  2.6. The second Timoshenko wave in the large frequency regime

If $|\Omega|$ is large enough, we saw that $k_{+}$ becomes independent of the
profile, while $k_{-}$ is given by (2.4). Hence the relevant question here is
to optimize the profile in order to minimize the reflection coefficient of the
second wave.

Denote by

$X(u):=\frac{8\kappa(1+\nu)\omega^{2}h_{1}^{2}\rho}{E_{0}}-\frac{48\kappa^{2}}{u^{2}},\quad
Y(u):=\frac{8\kappa(1+\nu)\omega^{2}h_{1}^{2}\rho}{E_{0}}\Big{(}\eta+\frac{3}{2}\frac{\delta}{uh_{1}}\frac{E_{1}}{E_{0}}\tilde{\eta}\Big{)}.$
(2.10)

Then using (2.4) we have

${\rm
Im}(k_{-}(\xi))\approx|\kappa-2(1+\nu)|^{-1/2}\frac{Y(h(\xi))}{2\sqrt{X(h(\xi))}}.$
(2.11)

Also, by neglecting the imaginary part we have:

$|k_{-}^{-2}(\xi)k_{-}(\xi)^{\prime}|^{2}\approx|\kappa-2(1+\nu)|\;\frac{48^{2}\kappa^{4}|h^{\prime}(\xi)|^{2}}{X(h(\xi))^{3}h(\xi)^{6}}$
(2.12)

By reasoning like in the previous case, and working with $n=1$ in the
functional, there must exist a constant $a$ such that

$|\kappa-2(1+\nu)|^{-1/2}\frac{Y}{2\sqrt{X}}+\beta^{-2}|\kappa-2(1+\nu)|\;\frac{48^{2}\kappa^{4}|h^{\prime}|^{2}}{X^{3}h^{6}}=a,$

or put differently, any $C^{1}$ increasing solution must obey

$h^{\prime}(\xi)=g_{a}(h(\xi)),\quad h(t)=1,$

where

$g_{a}(u):=\Big{(}a-|\kappa-2(1+\nu)|^{-1/2}\frac{Y(u)}{2\sqrt{X(u)}}\Big{)}^{1/2}\frac{\beta
X(u)^{3/2}u^{3}}{48\kappa^{2}|\kappa-2(1+\nu)|^{1/2}},\quad h_{\ell}\leq u\leq
1.$

When $n>1$, the formula for $g_{a}$ becomes:

$g_{a}(u):=(2n-1)^{1/(2n)}\Big{(}a-|\kappa-2(1+\nu)|^{-1/2}\frac{Y(u)}{2\sqrt{X(u)}}\Big{)}^{1/(2n)}\frac{\beta
X(u)^{3/2}u^{3}}{48\kappa^{2}|\kappa-2(1+\nu)|^{1/2}}.$

##  3\. Numerics

###  3.1. Comparison between Euler-Bernoulli and the first Timoshenko wave

The dimensionless parameters of the wedge (see Fig. 3) are as follows

$\displaystyle t=14.4,\quad h_{\ell}=0.0016,\quad n=1,$

and the material parameters are those of steel. The loss factor of the wedge
is taken as $\eta=0.001$. The loss factor of the visco-elastic layer
$\tilde{\eta}=0.25$, while $\delta/h_{1}=0.0012$ and
$\beta_{2}=E_{1}/E_{0}=0.3$. Lastly we take $\beta=\sqrt{2}/5$ as our penalty
parameter.

For clarity, we here present the results in terms of both dimensional and
dimensionless frequencies. As we have already explained in the introduction,
working with these parameters makes that the second Timoshenko wave (the one
with $k_{-}$) only appears at very high ultrasound frequencies, when the first
Timoshenko wave number (the one with $k_{+}$) becomes independent of the
profile. That is why the predictions of the Euler-Bernoulli theory, valid at
low frequencies, can only be compared with those of the first Timoshenko wave.

Figure 4. Profiles for $\Omega$ equal to $0.0025$ (dashed), $0.0124$ (dotted),
$0.0248$ (dash-dotted), $0.0495$ (line), where both the length and height are
dimensionless.

Figure 5. Zoomed in profiles for $\Omega$ equal to $0.0025$ (dashed), $0.0124$
(dotted), $0.0248$ (dash-dotted), $0.0495$ (line), where both the length and
height are dimensionless.

In Fig. 5 and Fig. 5, we can see that the optimal profile does not change much
for $\Omega$ between $0.0025$ and $0.0495$ (1000Hz and 20000Hz) although the
profile does get a little steeper as the frequency increases. The critical
behaviour seen in [16] for the Euler-Bernoulli beam theory, where the optimal
profile has a constant part, also shows op in Timoshenko’s framework, but for
much larger frequencies. This is caused by the appearance of an
$\tilde{a}_{*}$ as explained before. With the material parameters we have used
here, one gets a constant part in the optimal profile for $\Omega\approx 0.06$
(24546Hz). For frequencies which lies in the range where the height profile
found using Timoshenko’s beam theory is smooth, the profile can be well
approximated by a second order polynomial. More precisely, the optimal profile
for steel and $\Omega\approx 0.495$ (20000 HZ) is well approximated by the
polynomial

$\displaystyle p(x)=0.004673x^{2}+0.001673x+0.003816.$

In the rest of this section we will refer many times to the Normalized Wave
number Variation (NWV) given by the left-hand side of (1.4). In Fig. 7\- 8 we
have plotted the NWV for the optimal profile for the following values of
$\Omega$, $0.0025$, $0.0124$ and $0.0495$, respectively (1000 Hz, 5000Hz and
20000Hz). The white graph is set at $NWV=0.3$, i.e. the upper bound for which
the constraint for the NWV is satisfied. Similar for all three profiles is
that the condition for the NWV is satisfied over the whole profile for
frequencies above $5500$Hz, but one can also note that the constraint is
satisfied for much lower frequencies in the low part of the profile, as seen i
Fig. 8.

It is interesting to note, that the optimal profile we find for $\Omega$ below
$0.0136$ (5500Hz) does not satisfy the pointwise NWV constraint all the way.
This is due to the fact that we have only worked with $n=1$, which for low
frequencies is not enough to transfer the integral penalty into a pointwise
one. We cannot implement a high $n$ with our current numerical approach, since
already at $n=5$, the constant $b=9.65\cdot 10^{-100}$, which numerically
cancels out the imaginary part in our functional and creates serious
instabilities. It would therefore be interesting to develop other numerical
methods to find the optimal profile for lower frequencies and see how the
optimal profile changes for these frequencies, if the penalty for violating he
NWV constraint is harsher.

In Fig. 10\- 11 the optimal profile using both the Euler-Bernoulli and the
Timoshenko beam theory is plotted for the following values of $\Omega$,
$0.0248$, $0.0371$ and $0.0495$, respectively (10000Hz, 15000Hz and 20000Hz).
Interestingly, the two profile coincides for frequencies up to $\Omega\approx
0.03$ (as seen eg. in Fig. 10) and then they start deviating from each other
as seen in Fig. 10 and 11, while the two theories agrees up to $\Omega\approx
0.3$

Figure 6. The NWV for optimal profile for $\Omega\approx 0.0025$.

Figure 7. The NWV for optimal profile for $\Omega\approx 0.0124$.

Figure 8. The NWV for optimal profile for $\Omega\approx 0.0495$

Figure 9. Wedgeprofile Euler-Bernoulli beam theory (blue) and Timoshenko beam
theory (red) at $\Omega\approx 0.0248$, where both the length and height are
dimensionless.

Figure 10. Wedge profile Euler-Bernoulli beam theory (blue) and Timoshenko
beam theory (red) at $\Omega\approx 0.0371$, where both the length and height
are dimensionless.

Figure 11. Wedge profile Euler-Bernoulli beam theory (blue) and Timoshenko
beam theory (red) at $\Omega\approx 0.0495$, where both the length and height
are dimensionless.

###  3.2. Considerations on the second Timoshenko wave

We have shown in Section 2.6 that the second wave stops being evanescent and
becomes relevant at frequencies for which the first wave becomes a Rayleigh
wave. In other words, while the first wave can no longer be influenced by a
clever choice of the profile, we can do something for the second wave.
Mathematically, the optimal profiles for the second wave will have the same
general features as for the first wave, which means that starting from some
critically high frequency, the profiles will become piecewise constant near
$\xi=0$, and also more and more steep. We decided not to add more plots
because no qualitatively new phenomenon appears in this case.

##  4\. Conclusions and open questions

We have considered the problem of determining the optimal profile for a wedge
in a beam which minimizes the reflection coefficient of the first Timoshenko
wave. This has been done within the mathematical framework of calculus of
variations. We analytically derived an associated Euler-Lagrange equation
which takes into account an integral penalty constraint on the normalized wave
number variation (NWV), and then numerically solved this Euler-Lagrange
equation in order to get the optimal profile. In the numerical investigation
where we used the parameters of steel and a height ratio $h_{\ell}=0.0016$ we
also show, that in the dimensionless frequency range $\Omega/h_{\ell}\ll 1$
where the Euler-Bernoulli and the Timoshenko beam theories certainly agree,
the optimal profiles also agree. The two profiles start to differ at
$\Omega>0.0371$ (15000Hz). One interesting remark is that this change happens
at $\Omega\approx 0.0371$ which is less than $0.3$, which is the threshold for
which the two theories coincides for a plate with constant height. This shows
that the transition between the two theories is more subtle in the variable
height case.

Finally, we identified two open questions, which we believe it would be
interesting to investigate in the near future:

  * •

First, to try to find some material (or meta-material) where the second
Timoshenko wave appears in the acoustic frequency range, and to numerically
investigate whether the optimal profile for this wave compares to the one for
the first wave.

  * •

Second, to find numerical methods for solving the Euler-Lagrange equation
which better take into account the pointwise NWV constraint, like for example
allowing a large $n$ in the integral penalty term.

## References

  * [1] J.D. Achenbach.  Wave Propagation in Elastic Solids.  North-Holland, Amsterdam, 1973.
  * [2] P. Feurtado and S. Conlon.  Investigation of boundary-taper reflection for acoustic black hole design.  J. Noise Cont. Eng., 63(5):460–466, 2015.
  * [3] I.M. Gelfand and S.V. Fomin.  Calculus of Variations.  Dover publications, 2000.
  * [4] O. Guasch, M. Arnela, and P Sanchez-Martin.  Transfer matrices to characterize linear and quadratic acoustic black holes in duct terminations.  Journal of Sound and Vibration, 395:65–79, 2017.
  * [5] V.V Krylov.  New type of vibration dampers utilising the effect of acoustic ’black holes’.  Acta Acust. United, 90:830–837, 2004.
  * [6] V.V. Krylov.  Propagation of plate bending waves in the vicinity of one- and two-dimensional acoustic black holes.  ECCOMAS Thematic Conference on Computational Methods in Structural Dynamics and Earthquake, pages 460–466, 2007.
  * [7] V.V. Krylov and F.J.B.S. Tilman.  Acoustic ’black holes’ for flexural waves as effective vibration dampers.  Journal of Sound and Vibration, 274:605–619, 2004.
  * [8] J.Y. Lee and W. Jeon.  Vibration damping using a spiral acoustic black hole.  The Journal of the Acoustical Society of America, 141(3):1437–1445, 2017.
  * [9] J. Miklowitz.  The Theory of Elastic Waves and Waveguides.  North-Holland Publishing Company, 1980.
  * [10] M.A. Mironov.  Propagation of a flexural wave in a plate whose thickness decreases smoothly to zero in a finite interval.  Sov. Phys. Acoust., 34:318–319, 1988.
  * [11] A. Pelat, F. Gautier, S.C. Conlon, and F. Semperlotti.  The acoustic black hole: A review of theory and applications.  Journal of Sound and Vibration, 476:1–24, 2020.
  * [12] I.H. Shames and C.L. Dym.  Energy and Finite Element Methods in Structural Mechanics.  Taylor & Francis, New York, 2003.
  * [13] O. Sigmund and K. Maute.  Topology optimization approaches.  Struct. Multidic. Optim., 48:1031–1055, 2013.
  * [14] S. Sorokin and C.J. Chapman.  A hierarchy of rational timoshenko dispersion relations.  Journal of Sound and Vibration, 330:5460–5473, 2011.
  * [15] N.G. Stephen and S. Puchegger.  On the valid frequency range of timoshenko beam theory.  Journal of Sound and Vibration, 297:1082–1087, 2006.
  * [16] B.B. Støttrup, S. Sorokin, and H.D. Cornean.  A rigorous approach to optimal profile design for acoustic black holes.  The Journal of the Acoustical Society of America, 149:447–456, 2021\.