# Improvement of the acoustic black hole effect by using energy transfer due to geometric nonlinearity

Vivien Denis, Adrien Pelat, Cyril Touzé, François Gautier

► To cite this version:

Vivien Denis, Adrien Pelat, Cyril Touzé, François Gautier. Improvement of the acoustic black hole effect by using energy transfer due to geometric nonlinearity. International Journal of Non-Linear Mechanics, 2017, A Conspectus of Nonlinear Mechanics: A Tribute to the Oeuvres of Professors G. Rega and F. Vestroni, 94, pp.134-145. 10.1016/j.ijnonlinmec.2016.11.012. hal-01442428

HAL Id: hal-01442428

https://ensta-paris.hal.science/hal-01442428

Submitted on 20 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Improvement of the acoustic black hole effect by using energy transfer due to geometric nonlinearity

V. Denis *¹, A. Pelat¹, C. Touzé² and F. Gautier¹

¹Laboratoire d'Acoustique de l'Université du Maine, UMR CNRS 6613, Avenue Olivier Messian, 72085 Le Mans, Cedex 09, France

²IMSIA, ENSTA ParisTech, CNRS, CEA, EDF, Université Paris-Saclay, 828 Boulevard des Maréchaux, 91762 Palaiseau Cedex, France

December 20, 2016

## Abstract

Acoustic Black Hole effect (ABH) is a passive vibration damping technique without added mass based on flexural waves properties in thin structures with variable thickness. A common implementation is a plate edge where the thickness is locally reduced with a power law profile and covered with a viscoelastic layer. The plate displacement in the small thickness region is large and easily exceeds the plate thickness. This is the origin of geometric nonlinearity which can generate couplings between linear eigenmodes of the structure and induce energy transfer between low and high frequency regimes. This phenomenon may be used to increase the efficiency of the ABH treatment in the low frequency regime where it is usually inefficient. An experimental investigation evidenced that usual ABH implementation gives rise to measurable geometric nonlinearity and typical nonlinear phenomena. In particular, strongly nonlinear regime and wave turbulence are reported. The nonlinear ABH beam is then modeled as a von Kármán plate with variable thickness. The model is solved numerically by using a modal method combined with an energy-conserving time integration scheme. The effects of both the thickness profile and the damping layer are then investigated in order to improve the damping properties of an ABH beam. It is found that a compromise between the two effects can lead to an important gain of efficiency in the low frequency range.

## 1 Introduction

Controlling vibration is a major concern in many industrial applications today. Classical methods for passive vibration damping include the use of heavy viscoelastic layer [1] or tuned mass dampers [2]. These techniques are efficient and largely studied but their implementation is often limited by an important added mass or a narrow frequency range. The Acoustic Black Hole (ABH) effect is a particular surface damping method adapted to mitigate vibrations on a wide frequency range without added mass. The name ABH in this context has been first introduced by Mironov [3] and Krylov [4] and is now used in an increasing number of references in the literature [5, 6, 7]. Although the name may not be completely adapted — as referring to an astrophysical phenomenon — it is kept in this paper in order to be consistent with the current developments around this particular vibration absorber.

As described by Mironov [8], the ABH effect takes advantage of the flexural wave properties in a plate with variable thickness. It can be shown that waves propagating in a beam or plate with thickness $h(x)$

$$h(x) = \epsilon x^m, \quad m \ge 2, \qquad (1)$$

slow down without being reflected and that their travel time in the profile tends to infinity if thickness tends to zero. Due to the small thicknesses involved, manufacturing this "Acoustic Black Hole" profile is difficult and the residual thickness at the edge may be too large to obtain the desired effect. The addition of a thin viscoelastic layer [4] can compensate the effect of the profile truncation and the reflection coefficient of the extremity is then

<footnote>* vivien.denis@univ-lemans.fr</footnote>
considerably reduced. Further theoretical works propose a more refined model for the reflection coefficient [5] of an ABH termination and demonstrate the resulting increase of the modal overlap [9]. Experimental evidences of the effect are numerous [10, 11, 12]. While it is very efficient at high frequencies, several works point out the inefficiency of the ABH technique in the low frequency regime [9, 13], which may be a drawback in its potential applications.

The minimum thickness of the ABH extremity needed to obtain an important damping is very small (about hundred microns) and the displacement amplitudes at the extremity become large with respect to the thickness. Structures presenting this behavior are usually studied in the framework of nonlinear vibrations [14]. To the authors' knowledge, all numerical and theoretical models of ABH are based on linear assumptions such that the effects of the large amplitudes are neglected. Nonlinear effects include couplings between linear eigenmodes of the structure [14, 15]. Such couplings could be used for transferring energy from low frequency regime, where the ABH is inactive, to high frequency regime where vibrations can be efficiently dissipated, thus improving the ABH performance.

This paper presents an investigation of the nonlinear behavior of ABH beams and evaluates the interest of these nonlinear effects to improve the ABH damping efficiency at low frequency. Sec. 2 reports an experimental highlight of nonlinear response of an ABH beam. Sec. 3 develops a von Kármán plate model with variable thickness for the ABH beam and its numerical resolution from a modal projection and an energy-conserving time-marching scheme. The results of the model are detailed in Sec. 4. Conclusions and perspectives are given in Sec. 5.

## 2 Experiment

This section is devoted to an experimental investigation of the nonlinear behavior of a typical ABH beam. Dedicated forced vibration experiments are conducted in order to reveal the nonlinearities which can induce a variety of typical nonlinear phenomena including wave turbulence, characterized by a broadband Fourier spectrum and an energy transfer from low to high frequency range.

### 2.1 Experimental setup

The experimental study concerns two aluminium beams : the first is uniform with dimensions 1.3 m × 20 mm × 5 mm. Given its thickness, this beam is awaited to behave linearly. The second selected beam, with dimensions 1.5 m × 20 mm × 5 mm, has an Acoustic Black Hole termination of length 0.35 m; the thickness of this extremity decreases quadratically down to 50 µm approximately (see the geometry of the sample in Fig. 3). A thin viscoelastic tape (for example an electric tape) is applied on the ABH extremity. It is approximately 50 µm thick and its length is 0.2 m.

The experimental setup consists in a shaker (LDS V201), an amplifier (LDS PA25E), a force sensor (PCB 208C03), an accelerometer (PCB 352C29) and an acquisition card (NI-USB4431). It is driven by a Labview interface and the results are exported and processed with MATLAB. The beam to be tested is vertically suspended and excited at its midpoint with the shaker. The accelerometer is placed at an arbitrary location on the beam in order to record the vibrations at a given point.

### 2.2 Nonlinear regimes

The reference and ABH beams are excited at 83 Hz with a forcing amplitude increasing from 0 to 17 N in 40 s. By increasing gradually the forcing amplitude, the transition from periodic to chaotic regime can be observed [16], and compared to the transition scenario given in [15]. Figs. 1(a-b) show spectrograms of output acceleration for the two configurations. The spectrum of the response for the uniform beam shows that most of the energy is concentrated in the excitation frequency at 83 Hz, which is an indication that the beam behaves linearly (see Fig. 1(a)). Additional components in the spectrum may be observed. Their amplitudes are small compared to the peak at the excitation frequency. They can be attributed to other elements in the experimental setup, or residual nonlinearities in the shaker. The behavior is significantly different for the ABH beam, see Fig. 1(b), the spectrum of which rapidly enriches with a very large number of harmonics present in the response. The response of the beam is periodic but nonlinear, with a spread of energy in a large harmonic series which, in this reported case, goes up to 3 kHz. Similar results have been obtained with numerous different forcing
frequencies. Assuming that the beam's material constitutive behavior is linearly elastic, the nonlinearity present in the system comes from large amplitude oscillations of the beam, with respect to the thickness, and is denoted as geometric nonlinearity in the literature [14]. It significantly influences the system's response.

Fig. 1(c) has been obtained for a different forcing frequency, namely 113 Hz. It reveals the case of a complete transition to the turbulent, chaotic behavior. A direct transition from the periodic regime, composed of numerous harmonics, to the turbulent regime, characterized by a broadband Fourier spectrum and a cascade of energy to the high frequencies [17, 18, 19], is observed after 30 seconds of increasing the forcing amplitude. This observation is in the line of the transition to turbulence scenario for thin vibrating plates, already commented and analyzed in details with experiments and numerical simulations in [15, 16]. A more peculiar feature is shown in Fig. 1(d), where the ABH beam excited at 102 Hz first transits to the turbulent regime, clearly observed between 10 and 20 seconds, and then, while the forcing amplitude still increasing, recovers a periodic motion, characterized by a Fourier spectrum composed of discrete peaks. Note that the regime is still highly nonlinear as the number of harmonics is very large with energy up to 8 kHz. This peculiar regime has been observed a significant number of times during all the experiments realized, so that it appeared meaningful to report this behavior even if it is not analyzed specifically. During these large amplitude periodic nonlinear oscillations, it has been observed that the ABH beam vibrates around a deflected position, different from the position of the beam at rest. Hence these oscillations can be interpreted as strongly nonlinear normal modes (NNM), i.e. periodic orbits, existing at large energy levels, and far from the rest (straight) position. They can be brought closer to earlier experimental observations shown on vibrating beams in [20]. In these strongly nonlinear regimes, the system may have a number of large-amplitude NNMs that are isolas in the phase space, and for some particular conditions, nonlinear oscillations can be trapped in one of these regimes, as observed in Fig. 1(d).

<figure><img src="image_2.png" /><figcaption>Figure 1: Spectrograms of output acceleration for the reference beam and the ABH beam excited at various frequencies. (a) reference beam, 83 Hz, (b-d) : ABH beam, excited respectively at 83 Hz (b), 113 Hz (c), and 102 Hz (d). The forcing amplitude is raised from 0 to 15 N in 40 s. Color scale in dB (ref. 1 m.s⁻²)</figcaption></figure>
## 2.3 Energy transfer between low and high frequencies

In order to better understand the nonlinearity and the benefit that can be gained from it, the beams are now excited with a filtered white noise in the 100-500 Hz range, and for several increasing amplitudes. This frequency band has been selected as being the band in which the ABH is generally not efficient enough, and where the nonlinearity can bring an improvement by transferring the energy to higher frequencies. The parameter $\sigma$ denotes the standard deviation of the exciting force for each case. For a white noise, $\sigma$ is proportional to the square root of the power. The Power Spectral Density (PSD) of the acceleration $a(t)$, normalized by $\sigma^2$, is defined as:

$$S_{aa}(f) = \frac{1}{\sigma^2} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} a^*(t)a(t - \tau)e^{-j2\pi f\tau} dt d\tau. \quad (2)$$

It is displayed in Fig. 2 for uniform and ABH beams. Fig. 2(a) shows that the uniform beam exhibits a linear behavior since the power is concentrated in the excitation range, except for few artefacts that may come from the setup and happen however at a much lower scale. Most importantly, the spectrum does not depend on the amplitude. On the contrary, in the ABH case (see Fig. 2(b)), the power largely leaks outside the excitation range for high levels of excitation. On the zoom of Fig. 2(c) a small reduction of the vibratory levels with the increase of amplitude can actually be observed for the ABH beam, more specifically on the two peaks around 360 and 450 Hz. This is interpreted as an energy transfer from low to high frequency domains, the latter being inherently more damped. This is particularly interesting as it could be used in order to increase the efficiency of the ABH treatment at low frequencies. Note that here the transfer is not important enough and the gain in the damping of vibrations is very limited. The aim of the present study is thus to develop a model of the ABH including geometric nonlinearity, in order to analyze how this transfer can be optimized for a better efficiency of the ABH in the low frequency range.

# 3 Model of nonlinear acoustic black hole

## 3.1 A von Kármán plate model with variable thickness

The ABH beam is modelled as a nonlinear von Kármán plate [14, 21] with variable thickness. It is emphasized that the plate profile is assumed to vary along the $x$ direction only in this study, hence plate characteristics such as mass density $\rho$ and thickness $h$ only depend on $x$. The equation of motion for the transverse displacement $w(x, y, t)$ can then be written, following the developments proposed in [22, 23], as:

$$\rho(x)h(x)\ddot{w} + \Delta(D(x)\Delta w) - (1 - \nu)L(D(x), w) = p + L(w, F), \quad (3)$$

where $D(x)$, $\rho(x)$, $h(x)$ and $\nu$ are respectively bending stiffness, mass density, thickness and Poisson ratio of the plate. The term $p(x, y, t)$ represents the exciting force per unit surface. The Airy stress function $F(x, y, t)$ respects the following compatibility equation:

$$\Delta\left(\frac{1}{E(x)h(x)}\Delta F\right) - (1 + \nu)L\left(\frac{1}{E(x)h(x)}, F\right) = -\frac{1}{2}L(w, w), \quad (4)$$

where $E(x)$ is the Young modulus of the plate. Finally $L(f, g)$ is the Monge-Ampère bilinear operator which can be written, in cartesian coordinates, as:

$$L(f, g) = f_{xx}g_{yy} + f_{yy}g_{xx} - 2f_{xy}g_{xy}, \quad (5)$$

where $f_{xx}$ denotes the second derivative of $f$ with respect to $x$.

The geometry of the plate is shown in Fig. 3. It is composed of three different regions. The first one on the right is a uniform region (primary structure) with constant thickness $h_0$. The second region is the ABH where the thickness decreases from $h_0$ to $h_t$, which represents the minimum thickness of the plate. An extension of the small thickness region is considered in the left part of the beam, which is parameterized by its length $l_{add}$ that can be tuned in order to increase the nonlinear behavior of the structure. The thickness $h_p(x)$ of the plate can thus be written as:

$$h_p(x) = \begin{cases} h_t & \forall x \le l_{add}, \\ h_t + (h_0 - h_t) \frac{(x - l_{add})^m}{l_{ABH}^m} & \forall x \in [l_{add}, l_{add} + l_{ABH}], \\ h_0 & \forall x \ge l_{add} + l_{ABH}, \end{cases} \quad (6)$$
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

Figure 2: Power spectral density of acceleration response to a 100-500 Hz white noise of (a) Uniform beam for $\sigma=0.39$ N (yellow), 0.62 N (red), 0.84 N (brown) and 1.4 N (black) and (b) ABH beam for $\sigma=0.26$ N (yellow), 0.50 N (red), 0.68 N (brown) and 1.1 N (black), normalized by the variance $\sigma^2$ of the excitation, in dB (ref. 1 m$^2$.s$^{-4}$.Hz$^{-1}$.N$^{-2}$). Vertical lines indicate the excitation range. (c) Zoom of (b) on the 100-500 Hz range.
where m is the exponent of the ABH profile and $l_{ABH}$ is the length of the region of variable thickness. The plate has a width b and a total length L. In all the cases tested in this study, the exponent m is selected as m=2.

The effect of the damping layer is represented with equivalent mechanical properties that modify the bending stiffness, the mass density and the thickness, following the model of Ross-Ungar-Kerwin [1]. As shown in Fig. 3, the added viscoelastic layer is assumed to be of thickness $h_l$, and located in the interval $[x_{l1}, x_{l2}]$. Its effect on the flexural rigidity is taken into account more easily in the frequency domain, where the damping is represented as an imaginary part of the complex bending stiffness $D^*(x)$, which can be written as:

$$D^*(x) = \begin{cases} \frac{E_p h_p(x)^3}{12(1-\nu^2)} (1 + j\eta_p) = D_p(x)(1 + j\eta_p), & \forall x \notin [x_{l1}, x_{l2}], \\ D_p(x) \left[ (1 + j\eta_p) + \frac{E_l}{E_p} \left( \frac{h_l}{h_p(x)} \right)^3 (1 + j\eta_l) \right. \\ \left. + \frac{3 \left( 1 + \frac{h_l}{h_p(x)} \right)^2 \frac{E_l h_l}{E_p h_p(x)} (1 - \eta_p \eta_l + j(\eta_p + \eta_l))}{1 + \frac{E_l h_l}{E_p h_p(x)} (1 + j\eta_l)} \right], & \forall x \in [x_{l1}, x_{l2}], \end{cases} \quad (7)$$

where j is the imaginary unit, $D_p(x)$, $E_p$, $\eta_p$ and $h_p(x)$ are respectively the bending stiffness, the Young modulus, the loss factor and the thickness of the plate without the layer, while $E_l$, $\eta_l$ and $h_l$ are the Young modulus, the loss factor and the thickness of the damping layer, respectively. The real bending stiffness $D(x)$ is defined as the real part of $D^*(x)$:

$$D(x) = \mathrm{Re}(D^*(x)), \quad (8)$$

which is finally inserted into the equations of motion (3)-(4) in order to take into account all these effects in the final model. The effect of the mass added by the presence of the viscoelastic layer is taken into account by introducing the following equivalent mass density $\rho(x)$:

$$\rho(x) = \begin{cases} \rho_p, & \forall x \notin [x_{l1}, x_{l2}], \\ \frac{\rho_p h_p(x) + \rho_l h_l}{h_p(x) + h_l}, & \forall x \in [x_{l1}, x_{l2}], \end{cases} \quad (9)$$

where $\rho_p$ and $\rho_l$ are the mass densities of the structure and the layer respectively. The total thickness $h(x)$ of the sandwich plate can finally be written as:

$$h(x) = \begin{cases} h_p(x), & \forall x \notin [x_{l1}, x_{l2}], \\ h_p(x) + h_l, & \forall x \in [x_{l1}, x_{l2}]. \end{cases} \quad (10)$$

The equivalent Young modulus of the sandwich plate that have to be inserted in (4) for the compatibility equation on the Airy stress function, is deduced from the real bending stiffness $D(x)$ and the total thickness as:

$$E(x) = \frac{12D(x)(1 - \nu^2)}{h(x)^3}. \quad (11)$$

<figure><img src="image_6.png" /><figcaption>Figure 3: Schematic view of the ABH plate thickness profile. $h_0$ is the thickness of the plate thickness in the uniform region, $h_t$ is the minimum thickness of the plate, m is the exponent of the ABH profile, $l_{ABH}$ is the length of the region of variable thickness and $l_{add}$ is the length of the extension of small constant thickness. The damping layer stands between abscissas $x_{l1}$ and $x_{l2}$.</figcaption></figure>

More concisely, the equations of motion can be written as:

$$\rho h \ddot{w} + \Box(D, w) = p + L(w, F), \quad (12a)$$

$$\diamond(A, F) = -L(w, w)/2, \quad (12b)$$
with $A(x) = 1/E(x)h(x)$, and □ and ◊ the spatial operators applied to w and F:

$$ \Box(D, w) = \Delta(D(x)\Delta w) - (1 - \nu)L(D(x), w), \quad (13) $$

$$ \diamond(D, w) = \Delta\left(\frac{1}{E(x)h(x)}\Delta F\right) - (1 + \nu)L\left(\frac{1}{E(x)h(x)}, F\right). \quad (14) $$

In the case of a plate with free edges, the boundary conditions write, on an edge with normal line $n$ and tangent one $\tau$ [9, 23]:

$$ -\frac{\partial}{\partial n} \left( D \left( \frac{\partial^2 w}{\partial n^2} + \nu \frac{\partial^2 w}{\partial \tau^2} \right) \right) - 2(1 - \nu) \frac{\partial}{\partial \tau} \left( D \frac{\partial^2 w}{\partial n \partial \tau} \right) = 0, \quad (15a) $$

$$ \frac{\partial^2 w}{\partial \tau^2} + \nu \frac{\partial^2 w}{\partial n^2} = 0, \quad (15b) $$

$$ F = \frac{\partial F}{\partial n} = 0. \quad (15c) $$

One can note that Eq. (15c) are not the classical free conditions but are sufficient conditions to impose an edge free of loads in the in-plane direction [23].

## 3.2 Modal approach

A modal approach is used for the spatial discretization of the problem, following the general procedure described in [24]. Here the difficulty comes from the plate's variable thickness which, in turn, gives rise to a more complex formulation for the linear operators given in Eqs. (12). The eigenmodes are computed for the unforced, conservative problem. For that purpose, the real bending stiffness $D$, defined in (8) is used. The conservative and linear associated problems to be solved finally read:

$$ - \rho h \omega^2 w + \Box(D, w) = 0, \quad (16) $$

and

$$ \diamond(A, F) = \zeta^4 F, \quad (17) $$

to which are added the boundary conditions (15). These are eigenvalue problems with solutions $(\omega_k, \Phi_k(x, y))$ and $(\zeta_n, \Psi_n(x, y))$, respectively. The displacement and the Airy stress function are then expressed as:

$$ w(\mathbf{x}, t) = \sum_{k=1}^{N_{\Phi}} \Phi_k(x, y) q_k(t), \quad (18) $$

and

$$ F(\mathbf{x}, t) = \sum_{n=1}^{N_{\Psi}} \Psi_n(x, y) \beta_n(t), \quad (19) $$

where $q_k(t)$ and $\beta_n(t)$ are the modal coordinates for the transverse displacement and the Airy stress function, respectively. The integers $N_{\Phi}$ and $N_{\Psi}$ represent the number of transverse and in-plane modes kept in the modal truncation.

It is emphasized that the modes are orthogonal with respect to the mass:

$$ \int_S \rho(x)h(x)\Phi_k(x, y)\Phi_p(x, y)dS = 0, \quad p \neq q, \quad (20) $$

and normed such that:

$$ \int_S \Phi_k(x, y)^2 dS = 1. \quad (21) $$

The projection on the linear modes gives two equivalent formulations of the problem. The formulation (22) is quadratic in $(q, \beta)$ and can be written as:

$$ \beta_p(t) = -\frac{1}{2\zeta_p^4} \sum_{i=1}^{N_{\Phi}} \sum_{j=1}^{N_{\Phi}} H_{i,j}^p q_i(t) q_j(t), \quad (22a) $$

$$ \ddot{q}_s(t) + \omega_s^2 q_s(t) = F_s(t) + \frac{1}{M_s} \sum_{k=1}^{N_{\Phi}} \sum_{n=1}^{N_{\Psi}} E_{k,n}^s q_k(t) \beta_n(t), \quad (22b) $$
with $M_s = \int_S \rho(x)h(x)\Phi_s(x, y)^2 dS$ the modal masses, $F_s = \frac{1}{M_s} \int_S p\Phi_s dS$ the modal force, and the coupling coefficients:

$$H_{i,j}^p = \int_S \Psi_p L(\Phi_i, \Phi_j) dS, \qquad (23)$$

and

$$E_{k,l}^s = \int_S L(\Phi_k, \Psi_l) \Phi_s dS. \qquad (24)$$

The formulation (25), where variables $\beta$ have been eliminated to keep only the modal amplitudes $q$ of the transverse displacement $w$, shows a cubic nonlinearity in $q$:

$$\ddot{q}_s(t) + \omega_s^2 q_s(t) = F_s(t) - \frac{1}{M_s} \sum_{k=1}^{N_{\Phi}} \sum_{m=1}^{N_{\Phi}} \sum_{n=1}^{N_{\Phi}} \Gamma_{k,m,n}^s q_k(t) q_m(t) q_n(t). \qquad (25)$$

The $\Gamma_{k,m,n}^s$ coefficients are here calculated from $H_{m,n}^l$ and $E_{k,l}^s$:

$$\Gamma_{k,m,n}^s = \sum_{l=1}^{N_{\Psi}} \frac{H_{m,n}^l E_{k,l}^s}{2\zeta_l^4}. \qquad (26)$$

Although completely equivalent, the two different formulations can be used for different purposes. As shown in section 3.5, an energy-conserving scheme can be adopted for a robust and accurate time integration of the system. As shown in [24], this scheme particularly suits the quadratic formulation. On the other hand, the cubic formulation (25) shows that in-plane variables $\beta$ can be eliminated from the procedure. As two different integers have been introduced in the modal truncation given in Eqs. (18)-(19), one can use the cubic formulation in order to verify the convergence of the solution as function of the number of in-plane modes $N_{\Psi}$. Indeed, monitoring the convergence of the terms of the tensor $\Gamma_{k,m,n}^s$ versus $N_{\Psi}$ gives the sufficient number of in-plane modes needed for convergence. This will be shown in more details in the result section, see 4.2.

One can remark that the properties $H_{i,j}^p = H_{j,i}^p$ and $H_{i,j}^p = E_{p,j}^i$ already described in [23, 24] for the coupling tensors are also valid in this case and allow one to reduce the number of computed coefficients and thus the computational burden associated with this preliminary step of offline calculations.

## 3.3 Computation of the linear modes

Due to the variable thickness, the solution to the eigenvalue problems (16-17) is particularly difficult and needs a dedicated numerical approach. The solutions are obtained using a finite difference method thoroughly described in [9]. The use of a uniform grid for a structure where the wavelength varies spatially can result in a high computational cost. This problem is solved using a coordinate change and a transformed problem, which is equivalent to using a grid adapted to the flexural wavelength variation. Since the wavelength $\lambda = 2\pi\sqrt[4]{Eh^2/12\rho(1-\nu^2)\omega^2}$ varies as the square root of the thickness according to classical plate theory [25], the coordinate change is selected as follows:

$$d\tilde{x} = \frac{dx}{\tilde{x}_{av}\sqrt{h(x)}}, \quad \text{with} \quad \tilde{x}_{av} = \int_{0}^{L} \frac{1}{\sqrt{h(x)}} dx, \qquad (27)$$

where $\tilde{x}$ defines the transformed coordinate.

First the displacement $w(x, y)$, the bending stiffness $D(x)$, the mass density $\rho(x)$ and the thickness $h(x)$ are approximated with grid functions $w_{n,q}$, $D_n$, $\rho_n$ and $h_n$, respectively, where $n$ and $q$ are the spatial indexes such that $x = n\Delta_{\tilde{x}}$ and $y = q\Delta_y$, with $\Delta_{\tilde{x}}$ and $\Delta_y$ the spatial steps along the $\tilde{x}$ and $y$ directions. Then the eigenvalue problem in transverse displacement (16) is approximated with a transformed discrete problem, which can be written as[26, 9]:

$$\frac{h_n^{-1/2}}{\tilde{x}_{av}^2} \delta_{\tilde{x}+} \left[ (\mu_{\tilde{x}-}h_n^{-1/2}) \delta_{\tilde{x}-} \left[ D_n \left( \frac{h_n^{-1/2}}{\tilde{x}_{av}^2} \delta_{\tilde{x}+} \left( (\mu_{\tilde{x}-}h_n^{-1/2}) \delta_{\tilde{x}-} w_{n,q} \right) + \delta_{yy} w_{n,q} \right) \right] \right. \\ + D_n \delta_{yy} \left[ \frac{h_n^{-1/2}}{\tilde{x}_{av}^2} \delta_{\tilde{x}+} \left( (\mu_{\tilde{x}-}h_n^{-1/2}) \delta_{\tilde{x}-} w_{n,q} \right) + \delta_{yy} w_{n,q} \right] \\ \left. - (1-\nu) \left( \frac{h_n^{-1/2}}{\tilde{x}_{av}^2} \delta_{\tilde{x}+} \left( (\mu_{\tilde{x}-}h_n^{-1/2}) \delta_{\tilde{x}-} D_n \right) \delta_{yy} w_{n,q} \right) - \rho_n h_n \omega^2 w_{n,q} = 0. \right] \qquad (28)$$
185 where $\delta_{\tilde{x}-}$, $\delta_{\tilde{x}+}$, $\mu_{\tilde{x}-}$, $\mu_{\tilde{x}+}$ and $\delta_{yy}$ are the finite difference operators defined in A.

The in-plane problem is treated in the same manner. The convergence of the method has been extensively studied in [9]. The results presented in this article have been obtained using a 2000 × 100 grid. It is reminded that finite difference methods (especially those using second-order approximations), although converged, still presents an inherent numerical dispersion [26]. Other numerical methods, such as Finite Element Analysis (FEA), are indeed also suitable to find the desired solutions, as it has for instance been recently done in the current context in [13] for circular ABH.

The resolution of the discrete eigenvalue problems yields the out-of-plane and in-plane modes and the associated eigenvalues. The coefficients $H_{i,j}^p$, $E_{k,l}^s$ and $\Gamma_{k,m,n}^s$ are computed with the expressions (23, 24, 26) from the numerical mode shapes, following the general procedure also described in [27] for finite element shell models.

195

## 3.4 Losses

In accordance with the general modal framework used to solve the problem, the losses are taken into account by using a modal damping ratio for each oscillator equation. Assuming they remain small, the modal damping ratios $\xi_s$ are obtained by solving the eigenvalue problem in displacement with complex bending stiffness:

$$ -\rho h \omega^{*2} w + \Box(D^*, w) = 0. \qquad (29) $$

200

According to the literature [28], the eigenvalues $\omega_k^*$ of the problem defined in Eq. (29) can be written as:

$$ \omega_s^* = -\bar{\omega}_s \left( \xi_s - j\sqrt{1 - \xi_s^2} \right), \qquad (30) $$

where $\bar{\omega}_s$ are the natural frequencies and $\xi_s$ are the modal damping ratios. The modal damping ratios are then directly used in the system (22a-22b), which can finally be written as:

$$ \beta_p(t) = -\frac{1}{2\zeta_p^4} \sum_{i=1}^{N_{\Phi}} \sum_{j=1}^{N_{\Phi}} H_{i,j}^p q_i(t) q_j(t), \qquad (31) $$

205

$$ \ddot{q}_s(t) + 2\xi_s\omega_s\dot{q}_s(t) + \omega_s^2 q_s(t) = F_s(t) + \frac{1}{M_s} \sum_{k=1}^{N_{\Phi}} \sum_{n=1}^{N_{\Psi}} E_{k,n}^s q_k(t)\beta_n(t). \qquad (32) $$

This method is advantageous since it yields a dissipative time evolution problem that can be solved in the same manner as for the conservative one. Note that modal damping have been introduced only on Eq. (32) for the transverse motions $w$. This can be justified by considering the cubic formulation given in Eq. (25). Indeed the in-plane motions can be condensed in order to use a closed formulation for the modal displacements $q_s$ only where damping has to be introduced. Secondly, inertia effects are not taken into account for the in-plane motions, which gives a second justifications for neglecting possible losses in Eq. (31).

## 3.5 Time integration scheme

The system of Eqs. (31-32) is solved using an energy-conserving scheme described in [29, 24], that is conditionally stable. The first step is to discretize the modal coefficients $q_s(t) \to q_s^n$ (where $t = n\Delta t$ and $\Delta t$ is the time step) and the derivation operators with the proper finite difference operators. The conservative scheme can then be written [24]:

$$ \delta_{tt} q_s^n + 2\xi_s \omega_s \delta_t q_s^n + \omega_s^2 q_s^n = \frac{1}{M_s} \sum_{k=1}^{N_{\Phi}} \sum_{l=1}^{N_{\Psi}} E_{k,l}^s q_k^n [\mu_t \cdot \beta_l^n] + F_s^n, \qquad (33) $$

$$ \mu_{t-} \beta_p^n = \frac{1}{2\zeta_p^4} \sum_{i=1}^{N_{\Phi}} \sum_{j=1}^{N_{\Phi}} H_{i,j}^p q_i^n [e_{t-} q_j^n], \qquad (34) $$

where the operators $\delta_{tt}$, $\delta_t$, $\mu_{t-}$, $\mu_t$, and $e_{t-}$ are defined in A.

The stability condition of such a scheme is [24]:

$$ f_e > \pi f_{\max}, \qquad (35) $$

210

where $f_e = 1/\Delta t$ is the sampling rate and $f_{\max}$ is the largest frequency retained in the modal truncation for the transverse motion $w$, i.e. $f_{\max} = f_{N_{\Phi}}$.
# 4 Numerical results

## 4.1 Parameters of the simulated beams

The model is used to analyze and compare the damping performance of thirteen different beams, denoted U, E0, E0D0, E1, E2, E3, E4, E4D0, E4D1, E4D2, E4D3 and E4D4, and represented schematically in Fig. 4. A special emphasis is put on their ability to improve the efficiency of the ABH in the low frequency regime. All the tested beams have common parameters that are gathered in Tab. 1, while the distinctive features of each of these configurations are given in Tab. 2. Configuration U is a uniform beam and will be used as a reference case for comparing the different effects brought respectively by the ABH, the nonlinear additional length and the damping layer. Beams E0 to E4 are ABH beams without viscoelastic layer. Their additional length $l_{add}$ is gradually increased in order to extend the interval of the beam where the thickness is minimal. Their behavior is expected to be more and more nonlinear. Configuration E0D0 has the same geometry than E0, the only difference being the presence of the viscoelastic layer. Finally configurations E4D0 to E4D4 are all based on E4, with the presence of the damping layer of increasing length. As a memory aid, E stands for extension, D stands for damping and the number immediately after the letter is an indication for the concerned length.

<table><caption>Table 1: Values of the parameters shared by all the configurations simulated with the numerical model.</caption><thead><tr><th colspan="7">Parameters of the beam</th></tr><tr><th>E<sub>p</sub> [GPa]</th><th>ρ<sub>p</sub> [kg m<sup>-3</sup>]</th><th>ν</th><th>η<sub>p</sub> [%]</th><th>h<sub>0</sub> [mm]</th><th>L [m]</th><th>b [mm]</th></tr></thead><tbody><tr><td>70</td><td>2700</td><td>0.3</td><td>0.2</td><td>5</td><td>1.5</td><td>20</td></tr><tr><th colspan="7">Parameters of the ABH</th></tr><tr><th>l<sub>ABH</sub> [mm]</th><th>h<sub>t</sub> [μm]</th><th>m</th><th></th><th></th><th></th><th></th></tr><tr><td>200</td><td>100</td><td>2</td><td></td><td></td><td></td><td></td></tr><tr><th colspan="7">Parameters of the damping layer</th></tr><tr><th>E<sub>l</sub> [GPa]</th><th>ρ<sub>l</sub> [kg m<sup>-3</sup>]</th><th>η<sub>l</sub> [%]</th><th>h<sub>l</sub> [μm]</th><th></th><th></th><th></th></tr><tr><td>7</td><td>1000</td><td>40</td><td>100</td><td></td><td></td><td></td></tr></tbody></table>

<table><caption>Table 2: Values of the parameters specific to the configurations simulated with the numerical model.</caption><thead><tr><th></th><th>U</th><th>E0</th><th>E0D0</th><th>E1</th><th>E2</th><th>E3</th><th>E4</th><th>E4D0</th><th>E4D1</th><th>E4D2</th><th>E4D3</th><th>E4D4</th></tr></thead><tbody><tr><th>l<sub>add</sub> [mm]</th><td>/</td><td>0</td><td>0</td><td>50</td><td>100</td><td>150</td><td>200</td><td>200</td><td>200</td><td>200</td><td>200</td><td>200</td></tr><tr><th>x<sub>l1</sub> [mm]</th><td>/</td><td>/</td><td>0</td><td>/</td><td>/</td><td>/</td><td>/</td><td>200</td><td>150</td><td>100</td><td>50</td><td>0</td></tr><tr><th>x<sub>l2</sub> [mm]</th><td>/</td><td>/</td><td>200</td><td>/</td><td>/</td><td>/</td><td>/</td><td>400</td><td>400</td><td>400</td><td>400</td><td>400</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Figure 4: Schematic representations of the ABH extremities of the simulated beams. As a memory aid, E stands for extension, D stands for damping and the number immediately after the letter is an indication for the concerned length.</figcaption></figure>

## 4.2 Convergence study

The convergence of the numerical scheme is first investigated. The main free parameters that have to be set are the number of modes kept in the truncation, both for transverse and in-plane modes, $N_Φ$ and $N_Ψ$, and the time step $\Delta t$, or, equivalently, the sampling rate $f_e$. In order to ensure a converged representation of the in-plane motions, the convergence of coefficient $\Gamma_{p,p,p}^p$ with respect to $N_Ψ$ is investigated. An example is shown in Fig. 6, for $p=29$ (which correspond to the transverse mode (23,0) with 23 nodal lines along $x$ and no nodal lines along $y$, the modes being sorted by increasing eigenfrequencies), for two different configurations: the reference beam U and the ABH E0 with a small length of decreasing thickness and no damping layer. The
<figure><img src="image_8.png" /></figure>

(a)

<figure><img src="image_9.png" /></figure>

(b)

<figure><img src="image_10.png" /></figure>

(c)

<caption>Figure 5: (a) transverse modal shape for mode (23,0) of the E0 beam. (c)-(d) in-plane modal shapes for E0, (b) mode $\Psi_{101,2}$ belonging to family 1 and (d) mode $\Psi_{61,3}$ belonging to family 2.</caption>
corresponding modal shape for E0 is plotted in Fig. 5(a).

In all cases, the coupling term $\Gamma_{p,p,p}^p$ for the simple configuration U converges very quickly: approximately 100 in-plane modes are sufficient in order to ensure a converged value and thus a fine representation of the in-plane motions. On the other hand, the convergence for the beam E0 is much slower, as evidenced in Fig. 6(a). Note also that the value of the coupling coefficient is also much larger in case E0 than U, which was awaited since the small thickness region is designed in order to increase the nonlinearity in the system. The first consequence is thus the values of the nonlinear coupling coefficients, which are substantially increased. The convergence of $\Gamma_{p,p,p}^p$ for mode $p=29$ is slow and behaves incrementally; which suggests that some in-plane modes contributes strongly to the final value of the coefficients while others do not. In fact, the convergence can be made much faster by sorting the in-plane modes by families instead of sorting them by increasing eigenfrequencies. Let us define in-plane family $n$ the family containing the modes with $n+1$ nodal lines across the $y$-dimension (including the boundaries). Two examples are shown in Fig. 5(c-d): mode $\Psi_{101,2}$ in Fig. 5(c) belongs to family 1 while mode $\Psi_{61,3}$ in Fig. 5(d) belongs to family 2. Fig. 6(b) shows the result of the computation when the in-plane modes are sorted by families. Remarkably enough, only modes of family 1 actually contribute to the value of the coupling term for this specific transverse mode $p=29$. The same remarks can be made for other value of $p$ and it is verified that these result hold for $N_{\Phi} = 200$ out-of-plane modes in the modal truncation. The sorting procedure allows to eliminate non-contributing families from the computation, which greatly reduces the computational cost.

<figure><img src="image_11.png" /><figcaption>Figure 6: Coupling term $\Gamma_{p,p,p}^p$ for configuration U (black dashed line) and E0 (red full line). (a) Convergence of $\Gamma_{p,p,p}^p$ for $p=29$ (which correspond to transverse mode (23,0)), as function of the in-plane modes, sorted by eigenfrequencies of increasing order. (b) same convergence, now as function of the in-plane modes sorted by families. Vertical line on (b) separates the first family of modes $\Psi$ of E0 from the others.</figcaption></figure>

The number of transverse modes $N_{\Phi}$ depends on the frequency band one would like to cover in the dynamical simulations. It depends on the nonlinear regime at hand and the largest frequency present in the spectrum of the solution. For all the presented simulations, $N_{\Phi}=200$ has been selected, which gives an accurate representation of the solution up to the frequency 7 kHz for the most nonlinear configuration E4. This number has been found
to be sufficient for the tested regimes. Finally, once the number of in-plane modes chosen, the sampling rate $f_e$ can be simply selected thanks to the stability rule given by Eq. (35). However a larger value has been retained with $f_e=100$ kHz in all the computations, in order to have a more refined time representation.

## 4.3 Nonlinear regimes

In this section, the behavior of the model is investigated with respect to a harmonic forcing:

$$p(x, y, t) = \delta(x - x_0)\delta(y - y_0)G(t)\sin(\Omega t), \quad (36)$$

localized at point $(x_0, y_0)$, with fixed frequency $\Omega$ and increasing amplitude $G(t)$, in order to show the ability of the model in retrieving the different nonlinear phenomena experimentally observed. Two configurations are selected for that purpose, U and E0, with a forcing frequency of 115 Hz. This forcing frequency has been arbitrarily chosen to underline how the nonlinear effects appear, it is nonetheless close to 120 Hz which corresponds to an eigenfrequency of configuration E0; it is also close to 109 Hz which corresponds to an eigenfrequency of configuration U. The forcing amplitude $G(t)$ is linearly increasing from 0 to 20 N in 20 s. The spectra of output displacement at the excitation point $(x_0=0.6 \text{ m}, y_0=0.01 \text{ m})$ are plotted in Fig. 7. Fig. 7(a) shows that the uniform beam behaves linearly as expected, since the spectrum only contains the excitation frequency. It can be observed in Fig. 7(b) that the spectrum of E0 contains two distinct regimes: there is first a periodic regime showing an enrichment of the spectrum with odd harmonics, which qualitatively reproduces the results of the experiment of Sec. 2.2. Note that, as compared to the experiments, numerical spectra display only odd harmonics while experimental spectra contain both even and odd harmonics. The presence of only odd harmonics in the simulation is related to the fact that the model contains only cubic nonlinearities. Experimentally, geometric imperfections are unavoidable in real plates, especially in the ABH case where the thickness at the termination is very small. The presence of these imperfections creates in turn quadratic nonlinearity, which gives rise to even harmonics [30, 31]. The second regime observed in Fig. 7(b) is a chaotic one, associated with the Wave Turbulence (WT) phenomenon, as discussed in Sec. 2.2. Note that the von Kármán model is a moderate perturbation of the motions around the equilibrium position, with linearized angles. Hence the model is not able to produce periodic orbits around other positions than the configuration at rest. Consequently in all the simulated cases, the return to periodic orbits sometimes observed experimentally when still increasing the forcing amplitude, has not been retrieved. The choice of not considering imperfections in the model has been made both for the sake of simplicity and also mainly because in turbulent vibratory regimes with large amplitude, the cubic nonlinearities dominate the quadratic ones, as shown theoretically and numerically in [32].

<figure><img src="image_12.png" /><figcaption>Figure 7: Spectrograms of output velocity at point $(x_0=0.6 \text{ m}, y_0=0.01 \text{ m})$ for a forcing at 115 Hz from 0 to 20 N in 20 s, for beams (a) U and (b) E0. Color scale in dB (ref. $1 \text{ m.s}^{-1}$)</figcaption></figure>
## 4.4 Low frequency ABH performance

In order to observe the potential for an increased efficiency of the ABH in low frequency, the beams U, E0 and E0D0 are excited with a filtered white noise on the 100-500 Hz range during 10 s, for several amplitudes, i.e. several standard deviation $\sigma$, as it is done experimentally in Sec. 2.3. The Power Spectral Density of the output velocity of beams E0 and E0D0 at the excitation point ($x_0=0.6$ m, $y_0=0.01$ m) is displayed in Fig. 8 for different amplitudes of excitation. For the covered ABH beam E0D0, experimental results are qualitatively reproduced since energy leaks outside the excitation range, and the transfer becomes more important with the amplitude level, without however affecting the resonance peaks in the excitation range. This energy transfer towards higher frequency (energy cascade due to the Wave Turbulence [15, 33]) is much more noticeable and important in the E0 case, since the level outside excitation range easily reaches the level within, and a larger frequency range is concerned (see Fig. 8(c)). At large amplitude, the level of the resonance peaks within the excitation range is slightly reduced because this energy is transferred to higher frequencies where it is dissipated. It is important to note that while very little damping is present in the beam, nonlinearities allow for an interesting reduction of vibration in this case.

In order to evaluate the performance gain of the ABH due to the aforementioned nonlinear effects, an indicator $I^{c,\sigma}$ is built based on the responses to the white noise excitation of Fig. 8. It is defined as:

$$I^{c,\sigma} = \int_{100}^{500} \frac{S_{aa}^{c,\sigma}(f)}{S_{pp}^{c,\sigma}(f)} df, \qquad (37)$$

where c refers to the current configuration, $S_{aa}$ is the PSD of acceleration output $a(t)$ and $S_{pp}$ is the PSD of the forcing $p(t)$. The comparison is drawn out with a reference case denoted as $I^{\text{ref}}$. This reference case is selected as the uniform beam, configuration U, excited with the minimum amplitude ($\sigma = 0.11N$), hence giving rise to the indicator denoted as $E^{c,\sigma}$:

$$E^{c,\sigma} = 10 \log_{10} \left( \frac{I^{c,\sigma}}{I^{\text{ref}}} \right). \qquad (38)$$

Indicator $E^{c,\sigma}$ is plotted in Fig. 9 for the three configurations U, E0 and E0D0. A quasi-constant value is obtained for the uniform beam, as it is expected. Although the PSD is averaged over 10 s, the variability of the white noise excitation induces that the same exact value cannot be obtained for a linear case at different excitation levels, as it is shown by the small variations of case U. The value of the indicator for beam E0 is constant for small $\sigma$ and then decreases, indicating that less energy is present in the 100-500 Hz range, as it was shown by the spectrum of Fig. 8(c). Finally, the indicator has a constant and smaller value for beam E0D0, reflecting a better performance due to the damping layer but also smaller nonlinear effects and no significant energy transfer towards the high frequency domain, as it is shown in Fig. 8(b). This first result shows that a compromise has to be found between two competing effects: nonlinearity and damping. A too large value of damping prevents the nonlinear effect to settle down and thus creates an overdamped situation with no dependence on vibration amplitude. The next sections are devoted to a quantitative study of these two effects in order to optimize the low-frequency efficiency of the ABH.

## 4.5 Effect of the thickness profile

The importance of the nonlinear effects is first investigated by simulating the responses of configurations E1 to E4, which have increasing additional length of constant and minimum thickness. Increasing this length should enhance the nonlinear behavior. Moreover, the modal density will also increase in the excited frequency band, hence favoring the turbulent behavior and the energy transfer.

The procedure described in Sec. 4.4 is thus applied to beams E1 to E4 with gradually increasing length of extensions, according to Fig. 4. The results for the indicator $E^{c,\sigma}$ are plotted in Fig. 10. Configuration E1 shows that adding a small extension positively affects the nonlinear behavior since the structure becomes more sensitive to the excitation level: Indicator $E^{c,\sigma}$ decreases with respect to the amplitude of excitation and reaches values close to -9 dB. Configurations E2 to E4 are even more interesting from this point of view since the indicator reaches values around -12 dB. However, the similar values obtained for these configurations seem to indicate a saturation effect, i.e. that a longer extension does not result in a better energy transfer or a better dissipation due to nonlinearities.

As known from wave turbulence theory, finite-size effects alters the energy cascade due to the fact that a continuum of wavevectors is not at hand anymore, so that the resonance relationships shall then be verified
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

Figure 8: Power spectral density of velocity response to a 100-500 Hz white noise of (a) U for σ=0.11 N (yellow), 0.98 N (red) and 2.2 N (black), (b) E0D0 for σ=0.11 N (yellow), 1.1 N (red) and 2.4 N (black), and (c) E0 for σ=0.13 N (yellow), 0.99 N (red) and 2.7 N (black), normalized by the variance σ² of the excitation, in dB (ref. 1 m².s⁻².Hz⁻¹.N⁻²).
<figure><img src="image_16.png" /><figcaption>Figure 9: Indicator of performance $E^{c,\sigma}$ for configurations U (black dashdotted line ○), E0 (black dashed line ∇), E0D0 (red full line □).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 10: Indicator of performance $E^{c,\sigma}$ for configurations with different additional lengths: case U (black dashdotted line ○), E0 (black dashed line ∇), E1 (brown full line ▷), E2 (red full line △), E3 (yellow full line ◁) and E4 (black full line ◁).</figcaption></figure>
on a discrete set composed of the modes of the system [19, 34, 35]. In our case where energy transfer has to be optimized, the best situation would then be to have a nonlinear termination as large as possible in order to ensure the best modal density in the low-frequency range, and thus optimize the energy flux toward high frequencies. This situation is out of range for the selected beams where the number of eigenmodes in the tested frequency band 100-500 Hz increases from 8 (E0 configuration) to 17 (configuration E4). Fig. 11 shows the energy spectra for beams E0 to E4, and for the largest excitation amplitude tested, $\sigma=2.6$ N, highlighting the fact that the energy is redistributed over the prominent peaks defined by the eigenfrequencies, when extending the additional length, which explains the observed saturation on the indicator. However, better performance could be awaited for two-dimensional ABH as described in [5, 36], where the extension of the minimum thickness region can be significantly larger.

<figure><img src="image_18.png" /><figcaption>Figure 11: Power spectral density of velocity response to a 100-500 Hz white noise of E0 (blue dashed line), E1 (brown dashdotted line), E2 (red full line), E3 (yellow full line), E4 (black thick full line) for $\sigma=2.6$ N, normalized by the variance $\sigma^2$ of the excitation, in dB (ref. $1 \text{ m}^2 \cdot \text{s}^{-2} \cdot \text{Hz}^{-1} \cdot \text{N}^{-2}$).</figcaption></figure>

## 4.6 Balance of damping and nonlinearity

As shown in the previous section, the additional length $l_{add}$ has to be chosen as large as possible in order to enhance the nonlinear effects. The length of the damping layer is now investigated, in order to find the best compromise between nonlinearity and damping, and to ensure an increasing performance of the ABH thanks to the nonlinear behavior.

In order to observe the competition between damping and nonlinearities in the ABH beam, different lengths of damping layer are added to configuration E4, resulting in configurations E4D0 to E4D4 (see Fig. 4). The behavior of the indicator $E^{c,\sigma}$ for the associated simulations with a filtered white noise excitation are plotted in Fig. 12. Beam E4D0 reveals an interesting combination since the indicator $E^{c,\sigma}$ is lower than for E4 for small excitation levels, due to the more important damping in the structure. When increasing the vibration amplitude, indicator $E^{c,\sigma}$ for E4D0 shows then a decreasing behavior indicating that the nonlinearity is activated and a better efficiency of the ABH is obtained. The effect can clearly be seen in Fig. 13, which reveals reductions of the PSD down to 14 dB in the 100-500 Hz range (see Fig. 13(b)), with the $\sigma^2$ normalization. More specifically one can observe that the resonance frequencies tend to larger values when increasing the nonlinearity, reflecting the hardening behavior of the plate. Secondly, all resonance peaks widen due to the appearance of the strongly nonlinear behavior associated with the chaotic dynamics, and globally decrease as energy has flown to higher frequencies.

When the damping layer length is increased, gradually covering the nonlinear extension, nonlinear effects tend to disappear, as shown in Fig. 12 for E4D2 to E4D4 configurations, where the indicator is seen to stay at a constant value whatever the excitation level. A longer layer, hence a more important damping, however results in a better performance of the ABH. This study shows that both damping and nonlinear effect can be decoupled as the additional length (or the length of the ABH region) and the length (or the material characteristics) of the
damping layer can be independently adjusted in order to find an optimal behavior, depending of the context.

<figure><img src="image_19.png" /><figcaption>Figure 12: Indicator of performance $E^{c,\sigma}$ for structures with different damping layer length : case U (black dashdotted line ◦), E4 (black full line♦) E4D0 (red dashdotted line ◊), E4D1 (brown full line ∇), E4D2 (brown dashed line ◁), E4D3 (yellow full line △) and E4D4 (yellow dashed line ▷).</figcaption></figure>

# 5 Conclusions and perspectives

This study is devoted to the nonlinear effects in ABH beams, and the potential benefit that can be obtained from using them as a vector for transferring energy from the low frequencies, where ABH is known to be usually inefficient, to the high frequencies; this transfer is provided by the wave turbulence regime and its associated energy cascade. Experimental results on an ABH beam have been first reported, evidencing that nonlinearities are activated, giving rise to typical nonlinear phenomena. The thinness of the ABH termination are the support for the geometric nonlinearity, given the amplitude of motions involved.

A von Kármán plate model with variable thickness has then been derived in order to study numerically the main effects involved, with the aim of optimizing the energy transfer and improving the efficiency of the ABH in the low frequency range. A full numerical model has been developed, including the geometry variations, the damping due to the added viscoelastic layer, and the geometric nonlinearity due to large-amplitude vibrations. A modal approach with an energy-conserving scheme has been used for the resolution, this method allowing one to ensure a very good accuracy and convergence of the model for reasonable computational burden, together with a fine description of the frequency-dependent losses present in the structure.

A parametric study on a number of selected geometries has then been led, in order to investigate how the two competing effects represented by both the nonlinearity and the damping, can be tuned for obtaining an important gain in efficiency for the ABH. A good compromise has been found by using a long additional termination such that the nonlinearity can develop into a wave turbulence regime, together with a moderate length of the damping layer. Indeed, using a too large value of the damping overdamps all vibrations and counteracts the nonlinear behavior such that no improvement in efficiency is obtained with increased amplitudes. On the other hand, letting the chaotic dynamics settle down allows for the energy transfer in high frequencies where vibrations are damped, and thus improves the efficiency as amplitude increases. More generally, the results presented in this study show that previous linear models for the ABH may not apply anymore to some particular configurations, or need to be used with great care as they neglect an important nonlinear effect present in the dynamics. Further work dealing with plates embedding two-dimensional circular ABH with nonlinear effects could also be very interesting and closer to possible applications in mechanical engineering. As noted in the course of this study, two-dimensional ABH can present larger extensions with minimal thickness, ensuring a larger modal density in
<figure><img src="image_20.png" /><figcaption>Figure 13: (a) Power spectral density of velocity response to a 100-500 Hz white noise for configuration E4D0 for $\sigma=0.12$ N (yellow), 0.81 N (red) and 2.6 N (black), normalized by the variance $\sigma^2$ of the excitation, in dB (ref. $1 \text{ m}^2 \cdot \text{s}^{-2} \cdot \text{Hz}^{-1} \cdot \text{N}^{-2}$). (b) Zoom of (a) on the 100-500 Hz range.</figcaption></figure>
the low-frequency range. This fact should optimize the energy cascade, and could be more advantageous for evacuating the energy toward high frequencies.

References

395

[1] D. Ross, E. E. Ungar, E. M. Kerwin, Damping of plate flexural vibrations by means of viscoelastic laminae, in: Structural damping, Pergamon Press Edition, J.E. Ruzicka, Oxford, 1960, pp. 49—57.

400

[2] D. I. G. Jones, A. D. Nashif, R.L. Adkins, The effect of tuned dampers on vibrations of simple structures, American Institute of Aeronautics and Astronautics Journal 5 (2) (1967) 310–315.

[3] M. A. Mironov, V. V. Pislyakov, One-dimensional acoustical "black hole", in: VI session of the Acoustical Society of Russia, 1997, pp. 339–342.

405

[4] V. V. Krylov, F. J. B. S. Tilman, Acoustic 'black holes' for flexural waves as effective vibration dampers, Journal of Sound and Vibration 274 (3-5) (2004) 605–619.

[5] V. B. Georgiev, J. Cuenca, F. Gautier, L. Simon, V.V. Krylov, Damping of structural vibrations in beams and elliptical plates using the acoustic black hole effect, Journal of Sound and Vibration 330 (11) (2011) 2497–2508.

410

[6] L. Zhao, Stephen C. Conlon, F. Semperlotti, Experimental verification of energy harvesting performance in plate-like structures with embedded acoustic black holes, in: Proc. Internoise 2015, San Francisco, CA, United States, 2015, pp. 4043–4050.

[7] W. Huang, H. Ji, J. Qiu, L. Cheng, Flexural wave focalization in plates with imperfect two-dimensional acoustic black hole, in: Proc. of Internoise 2016, Hamburg, Germany, 2016, pp. 2392–2402.

[8] M. A. Mironov, Propagation of a flexural wave in a plate whose thickness decreases smoothly to zero in a finite interval, Soviet Physics: Acoustics 34 (3) (1988) 318–319.

[9] V. Denis, A. Pelat, F. Gautier, B. Elie, Modal overlap factor of a beam with an acoustic black hole termination, Journal of Sound and Vibration 333 (2014) 2475–2488.

415

[10] D. J. O'Boy, V. V. Krylov, V. Kralovic, Damping of flexural vibrations in rectangular plates using the acoustic black hole effect, Journal of Sound and Vibration 329 (22) (2010) 4672–4688.

[11] E. P. Bowyer, Daniel J. O'Boy, V. V. Krylov, J.L. Horner, Effect of geometrical and material imperfections on damping flexural vibrations in plates with attached wedges of power law profile, Applied Acoustics 73 (5) (2012) 514–523.

420

[12] V. Denis, F. Gautier, A. Pelat, J. Poittevin, Measurement and modelling of the reflection coefficient of an acoustic black hole termination, Journal of Sound and Vibration 349 (2015) 67–79.

[13] S. C. Conlon, J. B. Fahnline, Numerical analysis of the vibroacoustic properties of plates with embedded grids of acoustic black holes, The Journal of the Acoustical Society of America 137 (1) (2015) 447–457.

[14] A. H. Nayfeh, D. T. Mook, Nonlinear oscillations, Wiley New York, 1995.

425

[15] C. Touzé, S. Bilbao, O. Cadot, Transition scenario to turbulence in thin vibrating plates, Journal of Sound and Vibration 331 (2012) 412–433.

[16] C. Touzé, O. Thomas, M. Amabili, Transition to chaotic vibrations for harmonically forced perfect and imperfect circular plates, International Journal of Non-Linear Mechanics 46 (2011) 234–246.

430

[17] G. During, C. Josserand, S. Rica, Weak turbulence for a vibrating plate: can one hear a Kolmogorov spectrum?, Physical Review Letters 97 (2006) 025503.

[18] A. Boudaoud, O. Cadot, B. Odille, C. Touzé, Observation of wave turbulence in vibrating plates, Physical Review Letters 100 (2008) 234504.

435

[19] O. Cadot, M. Ducceschi, T. Humbert, B. Miquel, N. Mordant, C. Josserand, C. Touzé, Wave turbulence in vibrating plates, in: Handbook of Applications of Chaos theory, Chapman and Hall/CRC Edition, Skiadas, C.H. and Skiadas, C., 2016, pp. 425–448.
[20] J. P. Cusumano, F. C. Moon, Chaotic non-planar vibrations of the thin elastica, part I: experimental observation of planar instability, Journal of Sound and Vibration 179 (2) (1995) 185-208.

[21] T. von Kármán, Festigkeitsprobleme im maschinenbau, Enzyclopadie der Mathematischen Wissenschaften 4 (4) (1910) 311-385.

[22] G. J. Efstathiades, A new approach to the large-deflection vibrations of imperfect circular disks using Galerkin's procedure, Journal of Sound and Vibration 16 (2) (1971) 231-253.

[23] O. Thomas, S. Bilbao, Geometrically nonlinear flexural vibrations of plates: In-plane boundary conditions and some symmetry properties, Journal of Sound and Vibration 315 (3) (2008) 569-590.

[24] M. Duceschi, C. Touzé, Modal approach for nonlinear vibrations of damped impact plates: application to sound synthesis of gongs and cymbals, Journal of Sound and Vibration 344 (2015) 313-331.

[25] A. Leissa, Vibration of plates, Acoustical society of America, London, 1993.

[26] S. Bilbao, Numerical sound synthesis : finite difference schemes and simulation in musical acoustics, John Wiley & Sons, Hoboken, 2009.

[27] C. Touzé, M. Vidrascu, D. Chapelle, Direct finite element computation of non-linear modal coupling coefficients for reduced-order shell models, Computational Mechanics 54 (2014) 567-580.

[28] M. Géradin, D. Rixen, Mechanical Vibrations: Theory and Applications to Structural Dynamics, John Wiley & Sons, 1997.

[29] S. Bilbao, A family of conservative finite difference schemes for the dynamical von Kármán plate equations, Numerical Methods for Partial Differential Equations 24 (1) (2008) 193-216.

[30] C. Camier, C. Touzé, O. Thomas, Non-linear vibrations of imperfect free-edge circular plates and shells, European Journal of Mechanics A/Solids 28 (2009) 500-515.

[31] G. L. Ostiguy, S. Sassi, Effects of initial geometric imperfections on dynamic behaviour of rectangular plates, Nonlinear Dynamics 3 (1992) 165-181.

[32] M. Duceschi, O. Cadot, C. Touzé, S. Bilbao, Dynamics of the wave turbulence spectrum in vibrating plates: A numerical investigation using a conservative finite difference scheme, Physica D 280-281 (2014) 73-85.

[33] M. I. Auliel, B. Miquel, N. Mordant, Wave turbulence buildup in a vibrating plate, European Physical Journal B 88 (2015) 276.

[34] E. Kartashova, Discrete wave turbulence, EPL (Europhysics Letters) 87 (4) (2009) 44001.

[35] B. Miquel, N. Mordant, Nonlinear dynamics of flexural wave turbulence, Phys. Rev. E 84 (2011) 066607.

[36] D. J. O'Boy, V. V. Krylov, Damping of flexural vibrations in circular plates with tapered central holes, Journal of Sound and Vibration 330 (10) (2011) 2220-2236.

A Finite difference operators

The spatial finite difference operators $\delta_{\tilde{x}-}$, $\delta_{\tilde{x}+}$, $\mu_{\tilde{x}-}$, $\mu_{\tilde{x}+}$ and $\delta_{yy}$ appearing in Eq. (28) are defined as [26]:

$$ \delta_{\tilde{x}-}w_{n,q} = \frac{w_{n,q} - w_{n-1,q}}{\Delta_{\tilde{x}}}, \quad (39) $$

$$ \delta_{\tilde{x}+}w_{n,q} = \frac{w_{n+1,q} - w_{n,q}}{\Delta_{\tilde{x}}}, \quad (40) $$

$$ \mu_{\tilde{x}-} = \frac{w_{n,q} - w_{n-1,q}}{2}, \quad (41) $$

$$ \mu_{\tilde{x}+} = \frac{w_{n+1,q} - w_{n,q}}{2}, \quad (42) $$

$$ \delta_{yy} = \frac{w_{n,q+1} - 2w_{n,q} + w_{n,q-1}}{\Delta_y}. \quad (43) $$
In Eq. (34), the finite difference operators in time are defined as:

$$ \delta_{tt}q^n = \frac{q^{n+1} - 2q^n + q^{n-1}}{\Delta t^2}, \qquad (44) $$

$$ \delta_t q^n = \frac{q^{n+1} - q^{n-1}}{2\Delta t}, \qquad (45) $$

$$ \mu_{t-}q^n = \frac{q^n + q^{n-1}}{2}, \qquad (46) $$

$$ \mu_t q^n = \frac{q^{n+1} + q^{n-1}}{2}, \qquad (47) $$

$$ e_{t-}q^n = q^{n-1}. \qquad (48) $$