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

ELSEVIER

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

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

# Vibration damping of fiberglass honeycomb sandwich panels based on the acoustic black hole effect

Alex Besse<sup>1,*</sup>, Patrick O'Donoughue, Omar Aklouche, Frédéric Ablitzer, Adrien Pelat<sup>1</sup>, François Gautier<sup>1</sup>

Laboratoire d'Acoustique de l'Université du Mans (LAUM), UMR CNRS 6613, Institut d'Acoustique - Graduate School (IA-GS), Avenue Olivier Messiaen, 72085 Le Mans, France

## ARTICLE INFO

**Keywords:**

Honeycomb sandwich panel

Acoustic Black Hole effect

Bending waves

Vibration damping

Analytical model

## ABSTRACT

An Acoustic Black Hole (ABH) is a scatterer, embedded in a panel, allowing passive vibration control without adding mass. In practice, it is achieved by means of a local reduction in thickness (axisymmetric pit with a parabolic profile) and the addition of a thin viscoelastic coating in a central region of uniform thickness (plateau). The vibration absorption induced by the ABH, enables the design of stiff, light and non-resonant panels. In this paper, a study of the ABH effect for a sandwich panel made of a honeycomb core and fiberglass skins, which gives rise to both bending and shear effects is conducted. In particular, it is shown that the equations of motion of the thick, symmetrical sandwich panel with variable characteristics is obtained within the framework of the Nilsson's theory. These equations lead to a sixth-order analytical model that is used to determine the dispersion curves of the sandwich and also to establish an analytical model of the ABH inserted in a sandwich panel. Experimental tests postproceed by two complementary inverse methods are used to characterize the sandwich's mechanical properties. In particular, the model powered by experimentally determined parameters is used to analyze the effect of shear on the ABH effect, and to demonstrate that shear-induced softening is beneficial to ABH-induced absorption.

## 1. Introduction

In transport engineering (aeronautical, ground or rail transport), reducing vehicle mass is a major concern, as it helps to limit energy consumption. However, lightening structures generally makes them more susceptible to higher vibration and noise levels. Vibration control under mass reduction constraints is therefore an important research topic. Among the possible strategies for increasing the damping of structures without increasing their mass, the Acoustic Black Hole (ABH) effect is a good candidate. An ABH in a beam or in a plate is achieved by means of a local reduction in thickness and the addition of a thin viscoelastic coating. Integrating such a device enables the design of stiff, lightweight and non-resonant panels. The mechanism of bending wave attenuation by ABH is effective over a wide frequency band and results from the combination of a wave trapping effect in the zone where the thickness of the structure is reduced, and a localized damping effect induced by the viscoelastic film located in the thin zone. This effect has been the subject of a recent and very active literature [1], with the first demonstration in beams by Mironov [2], based on geometrical acoustics and subsequent developments by Krylov [3]. This pioneering work has its limitations in ignoring the effect of evanescent waves, prompting the development of more refined models by rewriting the equation for bending

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: Alex.Besse@univ-lemans.fr (A. Besse).</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. (a) Honeycomb fiberglass panel (b) ABH consisting of a pit of power law profile embedded in a uniform panel.</figcaption></figure>

wave propagation in a beam as a Riccati equation [4], leading to the exact calculation of the impedance and reflection matrices of the ABH termination. The latter wave-based model has been extended to rectangular ABH plates considered as waveguides (1D configuration), in order to take into account local imperfections and tears in the ABH edge, which favor wave scattering within the region of varying thickness and, consequently, wave attenuation [5]. A modal approach has also been proposed to investigate the role of nonlinearities in ABH beams due to large-amplitude vibrations [6] or impacts with a rigid support [7], and the potential advantage that can be gained by using them as a vector for transferring energy from low frequencies, where ABH is known to be generally inefficient, to high frequencies. The case of circular ABH embedded in plate (2D configuration) can be treated by wave approaches, enabling the modeling of ABH behavior in response to a point force [8] or in response to an incident plane wave [9]. The effect of the viscoelastic layer is linked to shear deformations within the thickness of the viscoelastic material as studied in [10]. In addition to these models describing the behavior of an ABH, several practical implementations have been studied. The effectiveness of the ABH effect as a function of the exact distribution of damping properties was studied in [11], where a temperature gradient was applied to the damping layer covering the ABH termination, showing that a controlled temperature gradient is beneficial to the ABH effect. Studies on the fine control of the bending stiffness of the ABH termination have also been carried out: an ABH 3D-printed with different materials was studied in [12], while an ABH placed in a fibrous structure with controlled fiber orientation was also examined in [13]. Sandwich structures are also becoming increasingly popular in a variety of industrial applications, particularly in the transport sector, as they offer the advantages of high rigidity and low weight. However, they introduce new challenges from a vibro-acoustic perspective due to their very low mass, which is detrimental to acoustic insulation, and a more complex analysis of vibration wave propagation in comparison with panels made of homogeneous material. This complexity arises from the existence of both bending waves and shear waves, as well as their possible coupling. These specificity have been extensively studied [14–18]. These studies have identified two behaviors in a sandwich panel. At low frequencies, pure bending motion is dominant and the shear effect can be ignored. At higher frequencies, the movements of the sandwich skins become out of phase and induce shear waves that propagate. At these higher frequencies, the coupling between bending and shear waves leads to a reduction in apparent bending stiffness as frequency increases. Studies on the insertion of ABHs into sandwich structures have already been conducted [19,20], showing that the induced attenuation is interesting. However, these studies are mainly focused on experimental exploration. In this paper, an investigation is performed in order to analyze how the ABH embedded in a sandwich as described in Fig. 1 can be enhanced by the softening induced by the shear effect. To this end, a hybrid model combining wave decomposition and finite difference discretization is developed. This model uses mechanical parameters estimated from a preliminary experimental test. This experimental characterization is derived from experimental data processed by two inverse methods, whose results are fitted to the Nilsson model [14,15]. The model is used to compare ABH-induced attenuation in panels with sandwich architecture or made of homogeneous material. First, Section 2 presents a procedure for sandwich characterization, i.e. for identifying the equivalent mechanical parameters. A wave-based model of the ABH panel is detailed in Section 3. Finally, the benefit of the softening sandwich structure due to shear effects on ABH performance is discussed in Section 4.

## 2. Mechanical properties of a honeycomb sandwich panel

### 2.1. Model of the bending stiffness

Considering a lightweight honeycomb panel used in aeronautical applications, composed of three superimposed layers comprising a Nomex honeycomb core and two fiberglass facings, as shown in Fig. 2(a). The panel comprises a core of thickness $H_c$, modeled as a uniform material where $E_c$ is the Young's modulus, $G_c$ the shear modulus, $\rho_c$ the density and $v_c$ the Poisson's ratio. The two fiberglass skins of thickness $h_s$ are assumed to be identical with Young's modulus $E_s$, density $\rho_s$ and Poisson's ratio $v_s$. Such a composite panel can be described as an equivalent single layer panel as proposed by Nilsson [14,15]. Nilsson's model concerns the vibration of a thick panel, taking into account both shear and bending movements [14,15,21] (see Fig. 2(a)).

In this model, the transverse motion of the two thin skins is described using Kirchhoff's plate theory [22]. The honeycomb core is modeled as a thick equivalent plate and the energy stored in this thick core is described by Mindlin's theory [22], for which the bending stiffness $D_s$ of one skin is defined as

$$D_s = E_s \frac{h_s^3}{12(1 - v_s^2)}, \qquad (1)$$
<figure><img src="image_5.png" /><figcaption>Fig. 2. (a) Side view of a three-layer composite panel. (b) Bending stiffness of an equivalent single leaf panel. The low-frequency limit of bending stiffness is $D_b$ (static bending stiffness). The high-frequency asymptotic value of the bending stiffness is $2D_s$ (two times the bending stiffness of a skin).</figcaption></figure>

the total static stiffness $D_b$ of the sandwich,

$$D_b = E_c \frac{H_c^3}{12(1 - v_c^2)} + \frac{E_s}{(1 - v_s^2)} \left( \frac{H_c^2 h_s}{2} + H_c h_s^2 + \frac{2}{3} h_s^3 \right), \quad (2)$$

the mass moment of inertia $I_ρ$,

$$I_ρ = \frac{\rho_c H_c^3}{12} + \rho_s \left( \frac{H_c^2 h_s}{2} + H_c h_s^2 + \frac{2h_s^3}{3} \right), \quad (3)$$

the mass per unit area $\sigma$ of the panel,

$$\sigma = 2h_s \rho_s + H_c \rho_c, \quad (4)$$

the slope in the x-direction is described by pure bending $\beta_x$ and by shear in the core $\gamma_{xz}$ giving,

$$\frac{\partial w}{\partial x} = \beta_x + \gamma_{xz}, \quad (5)$$

and the slope in the y-direction can be written in a similar manner,

$$\frac{\partial w}{\partial y} = \beta_y + \gamma_{yz}, \quad (6)$$

The equation governing the transverse vibration $w(x, t)$ of the sandwich structure is obtained by an energetic approach using Hamilton's principle. This principle uses kinetic, $\mathcal{T}$, and potential energy, $\mathcal{V}$, of the sandwich plate. The kinetic energy given by the lateral motion of the plate and the kinetic energies caused by rotation around the x- and y-axes are given by,

$$\mathcal{T} = \frac{\sigma}{2} \left(\frac{\partial w}{\partial t}\right)^2 + \frac{I_\rho}{2} \left[ \left(\frac{\partial \beta_x}{\partial t}\right)^2 + \left(\frac{\partial \beta_y}{\partial t}\right)^2 \right]. \quad (7)$$

The potential energy can be split into two categories, $\mathcal{V} = \mathcal{V}_c + \mathcal{V}_s$. With $\mathcal{V}_c$ corresponding to the bending and shear of the core and $\mathcal{V}_s$ represents the bending of both skins caused by the shear deformation of the core. The potential energy induced by pure bending of the structure is given by,

$$\mathcal{V}_{c1} = \frac{D_b}{2} \left[ \left( \frac{\partial \beta_x}{\partial x} \right)^2 + \left( \frac{\partial \beta_y}{\partial y} \right)^2 + 2\nu_c \left( \frac{\partial \beta_x}{\partial x} \right) \left( \frac{\partial \beta_y}{\partial y} \right) \right], \quad (8)$$

$$\mathcal{V}_{c2} = \frac{D_b}{4} (1 - \nu_c) \left( \frac{\partial \beta_x}{\partial y} + \frac{\partial \beta_y}{\partial x} \right)^2, \quad (9)$$

and the potential energy caused by shear in the core,

$$\mathcal{V}_{c3} = \frac{G_c H_c}{2} \left[ \left( \frac{\partial w}{\partial x} - \beta_x \right)^2 + \left( \frac{\partial w}{\partial y} - \beta_y \right)^2 \right]. \quad (10)$$

Finally,

$$\mathcal{V}_c = \mathcal{V}_{c1} + \mathcal{V}_{c2} + \mathcal{V}_{c3}, \quad (11)$$

As for $\mathcal{V}_s$, it is written as,

$$\mathcal{V}_s = \mathcal{V}_{s1} + \mathcal{V}_{s2} \quad (12)$$

with

$$\mathcal{V}_{s1} = D_s \left[ \left( \frac{\partial \gamma_{xz}}{\partial x} \right)^2 + \left( \frac{\partial \gamma_{yz}}{\partial y} \right)^2 + 2\nu_l \left( \frac{\partial \gamma_{xz}}{\partial x} \frac{\partial \gamma_{yz}}{\partial y} \right) \right], \quad (13)$$
$$V_{s2} = \frac{D_s(1 - v_s)}{2} \left( \frac{\partial \gamma_{xz}}{\partial y} + \frac{\partial \gamma_{yz}}{\partial x} \right)^2 . \qquad (14)$$

Using these formulations of energy and Hamilton's principle it is possible to obtain the following differential equation of the transverse motion of the sandwich panel $w(x,t)$ as done in [15].

$$-2D_b D_s \nabla^2[\nabla^2(\nabla^2 w)] + 2D_s I_\rho \nabla^2(\nabla^2 \frac{\partial^2 w}{\partial t^2}) + G_c H_c D_b \nabla^2(\nabla^2 w) - [(D_b + 2D_s) \sigma + G_c H_c I_\rho] \nabla^2 \frac{\partial^2 w}{\partial t^2} \\ + G_c H_c \sigma \frac{\partial^2 w}{\partial t^2} + \sigma I_\rho \frac{\partial^4 w}{\partial t^4} = 0. \qquad (15)$$

Applying Fourier transforms in space and time to equation (15), the following dispersion equation is obtained:

$$2D_b D_s k_f^6 - (2D_s I_\rho \omega^2 - G_c H_c D_b) k_f^4 - [\sigma (D_b + 2D_s) + I_\rho G_c H_c] k_f^2 \omega^2 - \sigma \omega^2 (G_c H_c - I_\rho \omega^2) = 0. \qquad (16)$$

In Kirchhoff's theory, the bending wave number is defined by

$$k_f^4 = \frac{\sigma \omega^2}{D}, \qquad (17)$$

with D the bending stiffness of the panel. Inserting (17) into (16) yields

$$2D_b D_s - (2D_s I_\rho \omega^2 - G_c H_c D_b) \frac{D(\omega)^{1/2}}{\omega \sqrt{\sigma}} - [\sigma (D_b + 2D_s) + I_\rho G_c H_c] \frac{D(\omega)}{\sigma} - (G_c H_c - I_\rho \omega^2) \frac{D(\omega)^{3/2}}{\omega \sqrt{\sigma}} = 0, \qquad (18)$$

where $D(\omega)$ is the equivalent bending stiffness of the single layer panel [23], i.e. the stiffness of a Kirchhoff plate that would have at any frequency the same wave number as that provided by Nilsson's model. The bending stiffness $D(\omega)$ calculated from Eq. (18) is represented in Fig. 2(b). For the simulation, material parameters are those given in Table 1 and geometrical parameters are those given in Section 2.2.4. The equivalent bending stiffness decreases with frequency, revealing two asymptotic limits. The low-frequency asymptotic value corresponds to $D_b$, the static bending stiffness of the entire sandwich panel. The high-frequency asymptotic value is $2D_s$, where $D_s$ is the bending stiffness of one isolated skin.

## 2.2. Experimental characterization

### 2.2.1. Setup and typical results

In order to build a semi-analytical model, the mechanical properties of the sandwich panel are identified using a sample panel. The panel studied experimentally is a 0.94 m × 0.84 m honeycomb sandwich panel made of two fiberglass skins with thickness $h_s = 0.28$ mm and a honeycomb core made of NOMEX-75 with thickness $H_c = 7.5$ mm. The cells of the core are hexagons with sides of $s_c = 5$ mm. This panel is characterized by vibrometric response measurements performed using a scanning laser vibrometer (Polytec PSV500 Xtra), an impedance head (PCB 288D01) and an electrodynamic shaker (Modal Shop 2075E). Fig. 3(a–b) shows

<figure><img src="image_6.png" /><figcaption>Fig. 3. (a) Experimental setup for vibration field measurements on a uniform sandwich panel using a shaker, an impedance head and a laser vibrometer. (b) Close-up view of shaker connection to plate. (c) Mobility at the center of the panel. (d), (e) and (f) are deflection shapes at 3 frequencies of interest.</figcaption></figure>
the experimental setup. The panel is excited by a periodic chirp in the range [25–10000] Hz. Its response is measured on a regular meshgrid of (81 × 73) = 5913 points, covering the entire surface of the panel. The mobility curve at the center of the panel is shown in Fig. 3(c). Three typical operating deflection shapes for low, mid and high frequencies are shown in Fig. 3(d–f). These results show different behaviors of the panel depending on the frequency of interest. In high frequency, only the propagating field can be observed due to the panel having a high intrinsic damping. The effect of the damping in the plate is also visible in the mid frequency range, showing a hybrid behavior of the vibration map: the modal field is starting to fade to leave only the propagating field.

To identify the mechanical properties of the experimental panel, Nilsson's model is fitted to the value of the identified bending stiffness. Two approaches are considered for the identification of bending stiffness as a function of frequency. The first one is an adaptation of the Force Analysis Technique. The second one is a wavenumber identification technique, Inhomogeneous Wave Correlation. Both of these techniques are applied considering a non-restrictive orthotropic plate model.

### 2.2.2. Force Analysis Technique (FAT)

The Force Analysis Technique (FAT) is an inverse method that was originally developed to identify external force distributions on simple structures such as beams and plates [24,25]. This method is based on the assumption that the measured displacement field obeys the local equation of motion of the structure, in which the spatial derivatives are approximated using a finite difference scheme. With the inverse method framework, it is also possible to solve for material parameters if the displacement field is measured in an area where no external excitation applies. In [26], this approach has been used to identify the frequency-dependent bending stiffness and principal orthotropy angle of composite panels. The same approach is applied here to the vibrometric measurements on the uniform sandwich panel.

The vibration field $\hat{w}(X, Y)$ measured under harmonic excitation is supposed to verify a local equilibrium equation of the form:

$$D_x \frac{\partial^4 \hat{w}}{\partial x^4} + 2 \sqrt{D_x D_y} \frac{\partial^4 \hat{w}}{\partial x^2 \partial y^2} + D_y \frac{\partial^4 \hat{w}}{\partial y^4} - \rho h \omega^2 \hat{w}(x, y) = 0, \quad (19)$$

where the coordinate system $(x, y)$ corresponds to the orthotropy directions of the material. Since the coordinate system $(X, Y)$ of the measurement mesh is possibly different from that of the material, a procedure to simultaneously find the orthotropy angle $\theta$ and the frequency-dependent bending stiffnesses $D_x$ and $D_y$ is implemented, following the procedure described in [26].

### 2.2.3. Inhomogeneous Wave Correlation (IWC)

Inhomogeneous Wave Correlation (IWC) is a wavenumber identification technique based on finding a damped (inhomogeneous) plane wave that is maximally correlated with the vibration field at a given frequency. This technique has been applied in the past to identify material parameters of a composite sandwich structure [23]. It assumes a harmonic displacement field $\hat{w}(x, y)$. An inhomogeneous wave $\tilde{w}_{k_C, \theta}$ is defined by a complex-valued wavenumber $k_C$ and direction $\theta$ in the coordinate system of the measurement mesh:

$$\tilde{w}_{k_C, \theta}(x, y) = e^{-j k_C (x \cos(\theta) + y \sin(\theta))}. \quad (20)$$

The use of a complex wavenumber $k_C = k(1 + j\gamma)$ enables this wave to encode both the bending wavenumber $k$ and the spatial damping $\gamma$. A correlation coefficient between the measured field and the inhomogeneous wave is then calculated using the following expression:

$$\mathrm{IWC}(k_C, \theta) = \frac{\left| \iint_S \hat{w} \tilde{w}_{k_C, \theta}^* dxdy \right|}{\sqrt{\iint_S \hat{w} \hat{w}^* dxdy \times \iint_S \tilde{w}_{k_C, \theta} \tilde{w}_{k_C, \theta}^* dxdy}}. \quad (21)$$

At each frequency $\omega$, material properties are identified through a two-step process. First, the values of the wavenumber $k_{IWC}(\theta, \omega)$ and the attenuation factor $\gamma_{IWC}(\theta, \omega)$ that maximize the correlation coefficient are found for each direction $\theta$ in the interval $[0, 2\pi]$. Then, the identified wavenumber $k_{IWC}(\theta, \omega)$ is inserted into the dispersion relation of a plate with elliptic orthotropy,

$$k_{IWC}^4(\theta, \omega) \left( \sqrt{\frac{D_x}{\sigma}} \cos^2 \theta + \sqrt{\frac{D_y}{\sigma}} \sin^2 \theta \right)^2 = \omega^2, \quad (22)$$

in which $\sigma$ is the mass per unit area. Finally, the two bending stiffnesses $D_x$ and $D_y$ are identified using a least-squares procedure [23].

### 2.2.4. Identification of mechanical parameters

To identify the mechanical parameters of the panel, the two inverse methods previously described are applied to the experimental data. As shown in Fig. 4(a), (b) and (c), a broadband identification of the bending stiffnesses is achieved by combining the information from the two methods, which have their own frequency ranges of validity. The IWC method is known to achieve accurate results in the mid- and high-frequency range but performs poorly at lower frequencies due to the method requiring several wavelength in order to accurately evaluate the correlation coefficient. In contrast, the validity of the FAT method is focused towards low and mid frequencies because the finite difference scheme used in the identification of the local behavior needs to stay small in regard of the wavelength as stated in [24–26]. Using both methods and an initial estimation of the mechanical properties it is possible to apply a least squares method of optimization onto the initial values to find the best fitted curve based on (18). Once the analytical bending stiffness fit the experimental data, the values of the mechanical properties are found. The data used for the fit are
<table><caption>**Table 1**</caption><caption>Identified material parameters of a honeycomb sandwich panel of uniform thickness from FAT and IWC methods.</caption><thead><tr><th></th><th>ρ</th><th>ν</th><th>E<sub>x</sub></th><th>G<sub>x</sub></th><th>E<sub>y</sub></th><th>G<sub>y</sub></th></tr></thead><tbody><tr><td>Aluminum</td><td>2700 kg/m<sup>3</sup></td><td>0.3</td><td>69 GPa</td><td>26.5 GPa</td><td>-</td><td>-</td></tr><tr><td>Fiberglass Skin</td><td>1541 kg/m<sup>3</sup></td><td>0.35</td><td>19.3 GPa</td><td>7.1 GPa</td><td>16.4 GPa</td><td>6 GPa</td></tr><tr><td>Homogenized honeycomb Core</td><td>37 kg/m<sup>3</sup></td><td>0.3</td><td>55 MPa</td><td>21 MPa</td><td>45 MPa</td><td>17 MPa</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** (a) Stiffness along the first orthotropy axis $D_x$, (b) Stiffness along the second orthotropy axis $D_y$, (c) Orthotropy angle estimation from IWC at one frequency, and (d) Combined orthotropy angle estimations from IWC and FAT.</figcaption></figure>

provided by the FAT method below 1 kHz, $D_{x/y}$ FAT, and completed by the IWC method, $D_{x/y}$ IWC, for the rest of the frequency range. The bending stiffness of the fit is $D_{x/y}$ Nilsson. The plate being orthotropic, bending stiffness in each directions, x and y, are experimentally found. In order to estimate the mechanical properties of the panel in both directions, two least square fit are performed with a tolerance of $\epsilon = 1 \cdot 10^{-13}$. The bending stiffness in both directions are presented in Fig. 4(a) and (b).

Both methods enable the identification of the orthotropic angle. In the IWC method, the orthotropic angle at each frequency is obtained from the wavenumber identified for each direction of propagation, as illustrated in Fig. 4(c). In the FAT method, the orthotropic angle is determined by systematically testing different angle values until the equation of motion (19) most closely fits with the experimental vibration field. Fig. 4(d) shows the angle identified by both methods at each frequency. The orthotropic angle of the plate $\theta = 6^\circ$ is obtained by averaging the identified values over the whole frequency range. The mechanical parameters of the equivalent homogeneous core and the skins, are estimated from Eqs. (19) and (1)-(4) and are synthesized in Table 1.

## 2.3. Equivalent rigidity of ABH profiles based on experimental characterization

Considering the unknown nature of the panel at the start of the study, the panel was assumed to be orthotropic in order to take a non-restrictive approach. The results given by Fig. 4 show the weak orthotropy of the panel, and thus, the orthotropy is ignored in the rest of the article. This choice is made in order to simplify the study of the sandwich panel as the orthotropic effect is not crucial to the rest of the study. Based on the previous experimental characterization giving the x-directions mechanical properties in Table 1 and the equivalent stiffness model in Section 2.1, this section presents the equivalent radial stiffness profile of a circular acoustic black hole placed in a honeycomb sandwich panel. In such a panel, ABH consists in locally thinning the honeycomb core according to a power law (quadratic law). The top skin of constant thickness is then formed to follow the shape of the bowl created by the
<figure><img src="image_8.png" /><figcaption>Fig. 5. (a) Bending stiffness versus frequency and radius for a solid ABH plate. (b) Bending stiffness versus frequency and radius for a sandwich ABH plate.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Side view of circular sandwich panel (at $\theta = 0$) with ABH placed at the center of the panel. The plate is divided into 5 zones in which the mechanical properties are different.</figcaption></figure>

honeycomb core. In Fig. 5, the frequency-dependent bending stiffness profile of the ABH sandwich panel (Fig. 5(b)) is compared with that of an ABH panel having the same geometry but a constant bending stiffness, equal to the static bending stiffness of the sandwich panel, as given by Eq. (2). In the following, this panel is called solid panel (Fig. 5(a)). As expected, the stiffness of the solid panel decreases locally by a factor of $10^4$ between the thick outer zone and the very thin core of the ABH. In the case of a sandwich panel, this radial reduction is identical and is complemented by a second, frequency-dependent reduction that is significant over the entire audible spectrum, due to the appearance of a shear movement above 1 kHz for the considered materials and dimensions. For the sandwich panel the wavelength at $f = 5$ kHz is $\lambda = 7.76$ cm. In regards to the size of a honeycomb cell $s_c$, it is more than 10 times larger meaning that the cells are not individually seen by the flexural wave. The effect of this additional stiffness reduction on the ABH effect is further investigated in Section 4.

### 3. Wave-based model of an ABH embedded in a honeycomb panel

This section focuses on establishing a semi-analytical model, based on the work of Aklouche [9], that considers a circular sandwich plate of radius $r_{ext}$ equipped with a circular acoustic black hole and loaded by a point force of amplitude $F$ at a location marked in the polar plane by the coordinates $(r_m, \theta_m)$ (see Fig. 6). Inside the ABH ($b < r < c$), the thickness profile varies radially according to the $m$th order power law given by the following equation:

$$h(r) = (h_0 - h_t) \left( \frac{r-b}{c-b} \right)^m + h_t, \qquad (23)$$

where $h_0$ is the total thickness of the plate outside the ABH, including the two skins and the honeycomb core, and $h_t$ is the truncation thickness at $r = b$. In practice a quadratic profile is used ($m = 2$)
## 3.1. Governing equations

As shown in Sections 2.1 and 2.3, the apparent bending stiffness of the plate is a function of space and frequency for the ABH zone (zone III) and only a function of frequency in the zones of constant thickness (zones I, II, IV and V). In zone III, where the thickness varies axisymmetrically, the displacement field W resulting from a point force loading satisfies the following bending equation of motion [22]:

$$ \nabla^2 [D_{III}(r, f) \nabla^2 W(r, \theta)] - (1 - v_{III}) \diamond^4 \{D_{III}(r, f), W(r, \theta)\} - \rho_{III} h(r) \omega^2 W(r, \theta) = F \frac{\delta(r - r_m)}{2\pi|\vec{r} - \vec{r}_m|} \delta(\theta - \theta_m), \quad (24) $$

where $D(r, f)$, $\rho_{III}$ and $v_{III}$ represent respectively the apparent bending stiffness, the density and the Poisson's ratio in zone III, and $\diamond^4$ is a partial derivative operator given as

$$ \diamond^4 \{D, W\} = \frac{\partial^2 D}{\partial r^2} \left[ \frac{1}{r} \frac{\partial W}{\partial r} + \frac{1}{r^2} \frac{\partial^2 W}{\partial \theta^2} \right] + \left[ \frac{1}{r} \frac{\partial D}{\partial r} \right] \frac{\partial^2 W}{\partial r^2}. \quad (25) $$

In the zones of constant thickness (zones I, II, IV and V), the bending stiffness is only frequency-dependent. In this case, the operator $\diamond^4$ becomes zero and Eq. (24) is modified to give the classical Love-Kirchhoff equation as:

$$ D_l(f) \Delta^2 W(r, \theta) + \rho_l h_l \frac{\partial^2 W(r, \theta)}{\partial t^2} = F \frac{\delta(r - r_m)}{2\pi|\vec{r} - \vec{r}_m|} \delta(\theta - \theta_m), \quad (26) $$

where $D_l(f)$, $\rho_l$ and $h_l$ ($i$ = I, II, IV, V) are the apparent bending stiffness, the density and the constant thickness in each zone.

### 3.1.1. Continuity and boundary conditions

To accurately represent the vibration field in the plate, continuity relations and boundary conditions are formulated using a state vector $\Phi_i$, defined for each zone as

$$ \Phi_i = \begin{Bmatrix} W^i(r, \theta) \\ \Theta^i(r, \theta) \\ M^i(r, \theta) \\ V^i(r, \theta) \end{Bmatrix}, \quad (27) $$

where $W^i(r, \theta)$ represents the bending displacement of zone $i$ (with $i$ = I, II, III, IV, V) and $\Theta^i(r, \theta)$, $M^i(r, \theta)$ and $V^i(r, \theta)$ are the slope, bending moment and shear force, given respectively by

$$ \Theta^i(r, \theta) = \frac{\partial W^i(r, \theta)}{\partial r}, \quad (28) $$

$$ M^i(r, \theta) = -D_l \left[ \frac{\partial^2 W^i(r, \theta)}{\partial r^2} + v_l \left( \frac{1}{r} \frac{\partial W^i(r, \theta)}{\partial r} + \frac{1}{r^2} \frac{\partial^2 W^i(r, \theta)}{\partial \theta^2} \right) \right], \quad (29) $$

$$ V^i(r, \theta) = -D_l \frac{\partial^3 W^i(r, \theta)}{\partial r^3} - \left( \frac{\partial D_l}{\partial r} + \frac{D_l}{r} \right) \frac{\partial^2 W^i(r, \theta)}{\partial r^2} + \left( \frac{D_l}{r^2} - \frac{v_l}{r} \frac{\partial D_l}{\partial r} \right) \frac{\partial W^i(r, \theta)}{\partial r} \\ - \left( \frac{D_l(2 - v_l)}{r^2} \right) \frac{\partial^3 W^i(r, \theta)}{\partial r \partial \theta^2} + \left( \frac{D_l(3 - v_l)}{r^3} - \frac{v_l}{r^2} \frac{\partial D_l}{\partial r} \right) \frac{\partial^2 W^i(r, \theta)}{\partial \theta^2}. \quad (30) $$

The state vector components $\Phi_i$ are continuous at the interfaces between the different zones. At the interface between zones IV and V (at $r = r_m$), the shear force $V^i(r, \theta)$ is discontinuous and must satisfy the following equation:

$$ V^{IV}(r_m, \theta) + F = V^V(r_m, \theta). \quad (31) $$

As the studied plate is free all around ($r = r_{ext}$), the bending moment and shear force are consequently equal to zero at the free edge:

$$ M^V(r_{ext}, \theta) = 0, \quad (32) $$

$$ V^V(r_{ext}, \theta) = 0. \quad (33) $$

### 3.1.2. Wave representation of the fields

In zones of constant thickness, the bending movement of the plate is described by Eq. (26). This equation admits homogeneous solutions that can be written in the form of a Fourier series such that [27]:

$$ W^i(r, \theta) = \sum_{n=0}^{\infty} W_n^i(r) \epsilon_n \cos(n\theta), \quad \begin{cases} \text{if } n = 0, & \epsilon_n = 1 \\ \text{else}, & \epsilon_n = 2 \end{cases}, \quad (34) $$

where

$$ W_n^i(r) = [u_n J_n(k_i r) + v_n I_n(k_i r)], \quad (35) $$
$$W_n^{II}(r) = [m_n J_n(k_{II}, r) + n_n Y_n(k_{II}, r) + s_n I_n(k_{II}, r) + t_n K_n(k_{II}, r)], \quad (36)$$

$$W_n^{IV}(r) = [e_n J_n(k_{IV}, r) + f_n Y_n(k_{IV}, r) + g_n I_n(k_{IV}, r) + l_n K_n(k_{IV}, r)], \quad (37)$$

and

$$W_n^V(r) = [a_n J_n(k_V, r) + b_n Y_n(k_V, r) + c_n I_n(k_V, r) + d_n K_n(k_V, r)], \quad (38)$$

represent the radial displacement field in each zone. The functions $J_n$ and $Y_n$ correspond to ordinary Bessel functions [28,29] of the first and second kind of order $n$ that can be interpreted as propagating waves. The functions $I_n$ and $K_n$ are the modified Bessel functions of the first and second kind, that can be interpreted as evanescent waves. Finally, $H_n^{(1)}$ is the first kind Hankel function of order $n$. These functions involve the flexural wave number $k_i$ (with $i = I, II, IV, V$) corresponding to each zone.

To obtain the displacement field in the zone III, Eq. (24) is numerically solved using the finite difference method. Same as in the constant thickness zones, the displacement field can be written as a Fourier series as follows

$$W^{III}(r, \theta) = \sum_{n=0}^{\infty} W_{n,j}^{III}(r) \epsilon_n \cos(n\theta) \text{ with } \begin{cases} \text{if } n=0, & \epsilon_n = 1 \\ \text{else}, & \epsilon_n = 2 \end{cases}, \quad (39)$$

in which, $W_{n,j}^{III}$ is the discrete displacement of a point $j$ ($j = 1, 2, ..., J$) inside the ABH. Inserting Eq. (39) in Eq. (24) gives

$$\lambda_4(r) \frac{\partial^4 W_{n,j}^{III}}{\partial r^4} + \lambda_3(r) \frac{\partial^3 W_{n,j}^{III}}{\partial r^3} + \lambda_2(r) \frac{\partial^2 W_{n,j}^{III}}{\partial r^2} + \lambda_1(r) \frac{\partial W_{n,j}^{III}}{\partial r} + \lambda_0(r) W_{n,j}^{III} = 0, \quad (40)$$

where

$$\lambda_0(r) = - \left[ \frac{n^2 v_{III}}{r^2 D_{III}} \frac{\partial^2 D_{III}}{\partial r^2} - \frac{3n^2}{r^3 D_{III}} \frac{\partial D_{III}}{\partial r} + \frac{n^2(4-n^2)}{r^4} + \frac{\rho_{III} h(r) \omega^2}{D_{III}} \right], \quad (41)$$

$$\lambda_1(r) = \left[ \frac{v_{III}}{r D_{III}} \frac{\partial^2 D_{III}}{\partial r^2} - \frac{2n^2 + 1}{r^2 D_{III}} \frac{\partial D_{III}}{\partial r} + \frac{2n^2 + 1}{r^3} \right], \quad (42)$$

$$\lambda_2(r) = \left[ \frac{1}{D_{III}} \frac{\partial^2 D_{III}}{\partial r^2} + \frac{2 + v_{III}}{r D_{III}} \frac{\partial D_{III}}{\partial r} - \frac{2n^2 + 1}{r^2} \right], \quad (43)$$

$$\lambda_3(r) = \left[ \frac{2}{D_{III}} \frac{\partial D_{III}}{\partial r} + \frac{2}{r} \right], \quad (44)$$

$$\lambda_4(r) = 1. \quad (45)$$

Finally, using a centered finite-difference scheme, with discretization step $\delta r$, a discrete form of Eq. (40) is obtained as

$$\Lambda_1(r_j) W_{j-2} + \Lambda_2(r_j) W_{j-1} + \Lambda_3(r_j) W_j + \Lambda_4(r_j) W_{j+1} + \Lambda_5(r_j) W_{j+2} = 0, \quad (46)$$

where

$$\Lambda_1(r_j) = \left[ \frac{\lambda_4(r_j)}{\delta_r^4} - \frac{\lambda_3(r_j)}{2\delta_r^3} \right], \quad (47)$$

$$\Lambda_2(r_j) = \left[ -4 \frac{\lambda_4(r_j)}{\delta_r^4} + \frac{\lambda_3(r_j)}{\delta_r^3} + \frac{\lambda_2(r_j)}{\delta_r^2} - \frac{\lambda_1(r_j)}{2\delta_r} \right], \quad (48)$$

$$\Lambda_3(r_j) = \left[ 6 \frac{\lambda_4(r_j)}{\delta_r^4} - 2 \frac{\lambda_2(r_j)}{\delta_r^2} + \lambda_0(r_j) \right], \quad (49)$$

$$\Lambda_4(r_j) = \left[ -4 \frac{\lambda_4(r_j)}{\delta_r^4} - \frac{\lambda_3(r_j)}{\delta_r^3} + \frac{\lambda_2(r_j)}{\delta_r^2} + \frac{\lambda_1(r_j)}{2\delta_r^3} \right], \quad (50)$$

$$\Lambda_5(r_j) = \left[ \frac{\lambda_4(r_j)}{\delta_r^4} + \frac{\lambda_3(r_j)}{2\delta_r^3} \right]. \quad (51)$$

Applying the continuity and boundary conditions described in Section 3.1.1 to the interfaces $r = b$ and $r = c$ implies the use of right and left finite-difference schemes respectively (instead of a centered scheme). All of the continuity and boundary conditions can be formulated as follows

$$M \cdot X = S, \quad (52)$$

where $M$ is a square matrix, $X$ is a vector grouping all the unknowns and $S$ a source vector. These quantities are described in detail in Appendix. Inverting Eq. (52) provide $X$ and so the unknown amplitudes $a_n, b_n, c_n, d_n, e_n, f_n, g_n, l_n, m_n, n_n, s_n, t_n, u_n, v_n$ and $W_{n,j}^{III}$.
# 4. Improvement of the ABH effect induced by sandwich shear effects

## 4.1. Validation of the wave-based model

This section is devoted to the validation of the wave-based model (WbM), by means of a comparison with results obtained using the commercial software Comsol Multiphysics®. The validation is carried out on a circular panel of radius $r = 0.5$ m, uniform thickness $h = 8.06$ mm, in which zones I to V are defined as shown in Fig. 6. In particular, zone III, which is uniform for this validation check, is nevertheless described using the finite-difference scheme given in Appendix. In this model, a circumferential order of $n = 20$ is taken into account, and 100 points are used for the finite difference scheme. The two plate configurations are used with the same parameters. Configuration 1 corresponds to an aluminum plate, with a frequency-independent loss factor $\eta = 0.001$. Configuration 2 corresponds to a circular sandwich plate, with the same total thickness and radius as configuration 1, consisting of a core $H_c = 7.5$ mm thick and two skins $h_s = 0.28$ mm thick. The loss factor $\eta = 0.04$, constant with frequency, is assumed to be identical in the core and skins. To describe this medium, an equivalent material whose Young's modulus depends on frequency is used. Recall that to obtain the equivalent Young's modulus, the parameters identified experimentally (see Table 1) for the core and skins are introduced into Eq. (18), whose resolution provides the equivalent stiffness $D(f)$, which is used to obtain the equivalent Young's modulus $E(f)$ by means of the relationship,

$$E_{eq}(f) = D(f) \frac{12(1 - v_{eq}^2)}{h_{eq}^3}. \qquad (53)$$

In this way, the weak orthotropy of the material, characterized experimentally in Section 2.2.4, is ignored and only the parameters relating to the x-direction in Table 1 are taken into account. The Poisson's ratio of the equivalent material is assumed to be $v_{eq} = 0.33$. The mass density of the material ($\rho_{eq} = 141$ kg/m³) is obtained using the mixing law Eq. (4).

For configuration 1 (aluminum plate), the response to a harmonic point force located at a distance $r = 0.25$ m from the center is shown in Fig. 7(a). The result obtained with a standard FEM simulation (FEM) using 2D shell elements meshed by 1600 structured quadrilateral, is in good agreement with the wave-based model considering Kirchhoff's plate assumptions (WbM), at low frequencies ($f < 2500$ Hz); discrepancies between the two models are clearly observed at higher frequencies ($f > 2500$ Hz). However, these discrepancies are linked to the limits of Kirchhoff's thin plate theory, on which the wave-based model is based. The plate can therefore be described as a thick plate, and the effects of shear and rotational inertia, not described by Kirchhoff's theory but taken into account by the FEM simulation, explain the discrepancies observed. Validation of the wave-based model using the FEM method can therefore not be carried out directly. To enable this validation, a similitude principle is proposed.

## 4.2. Principle of similitude for flexural waves

Due to limitations of the commercial software Comsol Multiphysics®, a direct validation of the WbM with FEM results cannot be performed, as the WbM was developed based on Kirchhoff's thin plate theory while the FEM uses Mindlin's plate theory. In order to fully validate the WbM, an alternative solution needs to be used. For this purpose, a principle of similitude can be established within the framework of Kirchhoff's theory: consider two plates of thickness $h$ and $ah$, with $\alpha$ an arbitrary scaling factor. These two plates, having the same arbitrary shape and the same arbitrary boundary conditions are excited by a harmonic point force $F$, applied at the same location. The Kirchhoff equation for the first plate can be written as

$$- \rho h \omega^2 w + \frac{E h^3}{12(1 - v^2)} \nabla^4 w = F, \qquad (54)$$

and for the second plate as:

$$- \rho(\alpha h)(\alpha \omega)^2 \frac{w}{\alpha^3} + \frac{E(\alpha h)^3}{12(1 - v^2)} \nabla^4 \frac{w}{\alpha^3} = F. \qquad (55)$$

Comparison of Eqs. (54) and (55) shows that the displacement of the plate of thickness $h$ at frequency $f$ is equal to the displacement of the plate of thickness $ah$ at frequency $\alpha f$ multiplied by the factor $\alpha^3$. Consequently, if $Y_h$ designates the input mobility of the plate of thickness $h$ and $Y_{ah}$ designates the input mobility of the plate of thickness $ah$, then these two quantities are related as follows:

$$Y_h(f) = \alpha^3 Y_{ah}(\alpha f). \qquad (56)$$

This principle of similarity between plates is illustrated in Fig. 7(d), which shows that the mobilities $Y_h$ and $Y_{ah}$ verify property (56). It is possible to take advantage of the similitude law in FEM simulation to guarantee that the thin-plate assumptions of Kirchhoff theory are verified even if the plate is thick. This special type of simulation is called FEM-K. It consists of simulating a thin plate of thickness $\alpha h$ ($\alpha < 1$), in order to have valid Kirchhoff kinematic assumptions. Using the similitude principle, Eq. (56), it is possible to obtain the mobility of the original thick plate as it would be described by Kirchhoff's plate theory.

In Fig. 7(a), the simulation using the similitude principle is used to obtain the curve FEM-K. This figure shows good agreement between the three models wave-based model (WbM), the Finite Element Model (FEM), and FEM-K at low frequencies where the Kirchhoff assumptions are valid for a plate of thickness $h = 8.06$ mm (see Fig. 7(b)). Going above the frequency of $f = 2500$ Hz, the model FEM shows disagreement due to the shear effect and rotational inertia effect due to large thickness (see Fig. 7(c)). This
<figure><img src="image_10.png" /><figcaption>Fig. 7. (a) Comparison of the driving point mobility of an aluminum panel of uniform thickness between the wave-based model and two FEM models. The solid black curve corresponds to the wave-based model, the dashed blue curve corresponds to the FEM model for a plate of thickness of h and the orange dotted line corresponds to the FEM model for a plate of thickness ah. (b) Zoom on low frequency peaks showing good agreement between the three models. (c) Zoom on high frequency peaks showing the shift of the standard FEM model and the good agreement between the WbM and the FEM-K model. (d) Illustration of the similitude principle: the mobility of a plate of thickness h is equal to the FRF of a plate of thickness ah that is shifted in amplitude by a factor of a³ and in frequency by a factor of a. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

disagreement is no more present if the similitude principle is used. While the real behavior of the panel is described by the FEM results, the good agreement between FEM-K and WbM allows for a validation of the wave-based model. The panel is accurately described as a Kirchhoff plate by the wave based model. The frequency error of the Kirchhoff model does not affect the rest of the study, since it is a comparative study using the same model and it does not impact the behavior of the ABH.

## 4.3. Sandwich model validation

The validation of the wave-based model, in the case of a thick sandwich plate, corresponds to the study of configuration 2. The sandwich is described by an equivalent material whose Young's modulus, $E(f)$, depends on frequency. Therefore, following the similitude principle, a Young's modulus $E(\alpha f)$ must be attributed to the plate of thickness $ah$ to carry out the simulation using the FEM and FEM-K method, for both of these models, the plate is represented using 2D shell elements meshed with 6241 structured quadrilaterals. A convergence study of the WbM is presented in Fig. 8. In this figure, the error, $\epsilon_n$, is plotted for each circumferential order $n$ at a few selected frequencies. The error $\epsilon_n$ is given by,

$$ \epsilon_n = \sqrt{\frac{(\int_S (W_{WbM} - W_{FEM-K}) dS)^2}{(\int_S W_{FEM-K} dS)^2}}, \quad (57) $$

with $W_{WbM}$ the displacement field given by the WbM, $W_{FEM-K}$ the displacement field given by FEM-K and S the surface of the plate. This result show that to accurately represents the displacement field using the WbM, at least 15 circumferential orders are needed. This convergence study is only shown for the sandwich panel as it is the more complex problem. In Figs. 7 and 9, 20 circumferential order are taken into account.
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Value of the error, $ε_n$, between the displacement field obtained by the WbM and the FEM-K for several frequencies at each circumferential order. The displacement field for $f = 1500$ Hz at $n = 3$ and $n = 17$ are shown. The results show that at $n = 20$ the model is fully converged for the whole frequency band.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Comparison of the driving point mobility obtained for WbM, FEM and FEM-K for a sandwich panel.</figcaption></figure>

Fig. 9 shows the input mobility of the panel excited at a distance $r_{ext} = 0.17$ m from the center by a harmonic force, calculated with the WbM model, the FEM model and the FEM-K model. Looking at the FEM results, the curve is in good agreement with the FEM-K and WbM in the low frequency. As expected, the frequency shifts appear sooner than for the aluminum plate due to the sandwich panel having a larger wavelength. Looking at FEM-K and WbM, the curves are almost identical, showing a good representation of the sandwich, in the framework of Kirchhoff's plate theory, by the WbM while reducing the computation time. The FEM-K takes 55 min for 1080 frequencies, while it only takes 40 min for the WbM model. Note that the computation time of the WbM could be further reduced by improving the software architecture of the model.

### 4.4. Effect of the shear movements on the ABH-induced absorption

The variable thickness plate model (see Section 4.1), in conjunction with the use of an equivalent material, allows to examine the benefits of integrating the ABH into a sandwich panel. The shear effect within the sandwich material tends to reduce the stiffness of the panel, which is also described by an equivalent Young's modulus that decreases with increasing frequency. This softening effect is expected to be beneficial to the ABH effect. To quantify the interest of this effect, a parametric study is carried out, described in Fig. 10. The system consists of the circular plate already studied for validation (thickness $h = 8.06$ mm, radius $r_{ext} = 0.5$ m), into which is inserted a black hole with the following characteristics: radius of zone I (central plate) $a = 0.0475$ m, radius of zone II (skin overlap) $b = 0.05$ m, radius of zone III (zone ABH) $c = 0.15$ m, radius of zone IV (excitation position) $r_m = 0.17$ m, radius of zone V (disc radius) $r_{ext} = 0.5$ m, thickness of the central plate $h_t = 0.28$ mm. In the center of the ABH (zones I, II, III), a viscoelastic coating is added, which locally increases damping. The loss factor of the coated sandwich is described by the RKU model [11,30], which describes a mixing law between the sandwich and the coating. The viscoelastic coating has a uniform thickness $h_{visco} = 0.05$ mm
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Mean quadratic mobility, as a function of frequency, and also calculated in one-third octave bands for solid panel and sandwich panel for two material damping values, $\eta = 5 \cdot 10^{-3}$ for (a) and (b) and $\eta = 0.04$ for (c) and (d). Figures (a) and (c) correspond to the mean quadratic mobility of a solid panel (with bending stiffness constant in frequency). Figures (b) and (d) correspond to the mean quadratic mobility of a sandwich panel (with bending stiffness dependant of frequency). The reduction of the mean quadratic mobility provided by the ABH is greater in configurations (b) and (d) than in configurations (a) and (c). This translate to better efficiency of the ABH when the bending stiffness is frequency dependant. Figures (e) and (f) are the displacement field, at $f = 2885$, Hz of the highly damped solid panel with and without ABH. Figures (g) and (h) are the displacement field, at $f = 2885$, Hz of the highly damped sandwich panel with and without ABH.</figcaption></figure>

and is made of a material with density $\rho_{visco} = 1000 \text{ kg/m}^3$, Young's modulus $E_{visco} = 7\text{e9 Pa}$, and loss factor $\eta_{visco} = 1.6$. The mean square mobility averaged over the whole plate ($\langle Y^2 \rangle$) is used as an indicator of the panel's vibratory level, considering that it is proportional to the kinetic energy of the plate when excited by a unit force.

Fig. 10 shows the variations of $\langle Y^2 \rangle$ as a function of frequency and the values averaged over one-third octave bands. The calculation is carried out for a circular panel of uniform thickness (reference case) and a circular panel incorporating an ABH at its center. In the first column of Fig. 10(a) and (c), the softening effect of shear movements is ignored. The stiffness therefore depends only on the thickness and not on the softening effect in frequency. In this case, it is equal to its static limit value, as described in Fig. 5(a). In the second column (Fig. 10(b), (d)), the softening effect is taken into account as described by Nilsson's model (equation (18)). The stiffness variations resulting from both thickness variations and shear effect are those described in Fig. 5(b).
In the first line of Fig. 10(a) and (b)), core and skin damping is low ($\eta = 5 \cdot 10^{-3}$). In the second line (Fig. 10(c) and (d)), core and skin damping is high ($\eta = 0.04$). In all cases, the coating damping, described by the RKU mixing law, is added.

All the simulations show that thickness variations induce a smoothing of the mobility peaks, which is characteristic of the ABH effect. This smoothing is primarily due to the substantial damping provided by the coating. A comparison of Figs. 10(a) and (b) shows that the absorption effect of the ABH is higher when embedded in a sandwich panel. This result also holds for panels with higher damping as shown by Figs. 10(c) and (d). By way of illustration, in the one-third octave band centered at 3500 Hz, the reduction in mobility averaged by third octave band, brought about by the ABH is 4 dB for the highly damped sandwich panel (Fig. 10(d)) whereas it is only 1 dB when this softening effect is ignored (Fig. 10(c)). On the other hand, the dynamic range of mobility, which reflects the resonant nature of the panel, is reduced by the black hole effect. By way of illustration, the dynamic range is of the order of 20 dB for the reference sandwich panel and only 12 dB for the ABH sandwich panel (see (Fig. 10(b)) whereas it is of the order of 17 dB for the reference solid panel and 11 dB for the ABH solid panel (see (Fig. 10(a)). This reduction in dynamic range is all the more noticeable when the panel is poorly damped as shown by the comparison of Fig. 10(b) and (d).

It should also be noted that softening induces an increase in the modal density of the panel, which has the effect of increasing the smoothing effect and therefore the ABH effect. The displacement fields, simulated at a frequency of interest (2885 Hz), show the focusing effect induced by the ABH as well as the local complexity of the vibration field resulting from the added damping. This focusing and local complexity are greater in the softer central zone, this is shown by the Figs. 10(e)–(h), in which the displacement field with and without ABH is plotted for a solid plate (see (e) and (f)) and for a sandwich plate (see (g) and (h)).

## 5. Conclusion

Sandwich structures are characterized by a frequency-dependent equivalent bending stiffness reduction that comes from the shear movement between the skins. The main objective of this paper was to study in what manner this frequency-dependent stiffness reduction impacts the Acoustic Black Hole (ABH) effect induced by a spatial stiffness gradient. In order to do so, the equivalent bending stiffness of a uniform sandwich panel has been identified by two inverse methods applied to vibration measurements. Using the identified mechanical properties, a wave-based model, considering Kirchhoff's thin plate theory, of an isotropic circular panel equipped with an ABH has been developed to explore the implications of the frequency-dependent stiffness reduction on the ABH effect. This model uses Nilsson's assumption to describe the behavior of a sandwich panel and is thus limited to thin plates as it is based on Kirchhoff's bending stiffness. It has also been developed for isotropic axisymmetric configurations and does not allow for comparison with the sample used in characterization. This limitation could be addressed by redeveloping the model with other geometrical considerations. Finally, this model has been used to compare the sandwich panel with a solid panel, both having the same static bending stiffness and being equipped with the same ABH.

The main contribution of this paper is the demonstration that the softening effect resulting from a sandwich structure is beneficial for ABH, as demonstrated using the wave-based model. This model has been validated using a FEM model on which a similitude principle was applied to ensure conformity with the assumptions of Kirchhoff's plate theory. This validation shows that the wave-based model correctly represents the stiffness reduction effect described by Nilsson's theory. Through the mean quadratic mobility of a panel with an ABH, it is shown that the attenuation effect of the ABH is higher in a sandwich panel than in a solid panel. The reduction provided by the ABH (compared to a reference panel of uniform thickness) in a solid panel is in the order to 4 dB, whereas in a sandwich panel the reduction goes up to 7 dB. This increased attenuation shows that the inclusion of an ABH in sandwich structures can be beneficial in terms of noise reduction, in particular in the context of aeronautics, where sandwich structures are widely used due to their low mass. The reduction in rigidity also highlights the phenomenon of a reduction in wavelength both outside and inside the ABH. Finally, as the wave-based model is relatively inexpensive in terms of computation time, it enables optimization processes to be envisaged in order to adapt the ABH design to specific applications.

### CRediT authorship contribution statement

**Alex Besse:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation. **Patrick O'Donoghue:** Writing – review & editing, Resources, Methodology, Investigation, Data curation. **Omar Aklouche:** Writing – review & editing, Resources, Methodology, Formal analysis, Conceptualization. **Frédéric Ablitzer:** Writing – review & editing, Supervision, Resources, Methodology, Investigation, Data curation. **Adrien Pelat:** Writing – review & editing, Supervision, Resources, Project administration, Methodology, Conceptualization. **François Gautier:** Writing – review & editing, Supervision, Resources, Project administration, Methodology, Funding acquisition, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

This work is part of the MAMBO project, financially supported by the Direction Générale de l'Aviation Civile (DGAC), 2021–2025. The authors sincerely thank F. Besluau and M. Roisin for preparing and supplying the sandwich panel, and Julien Nicolas and Mathieu Sécail-Géraud for their fruitful contributions to the experiments.
Appendix. Matrix form of continuity and boundary conditions

The problem solved in Section 3 consists in solving equation (24) which describes the bending motion of a finite Love–Kirchhoff plate with an acoustic black hole. The plate is divided into 5 zones as shown in Fig. 6. Eq. (24) admits analytical solutions for zones I, II, IV and V and are given by Eqs. (34) to (38). The vibratory field in zone III is given by Eq. (39) and is obtained numerically using the finite-difference method. The continuity of the state vector at the interfaces of each zone and boundary conditions (see Section 3.1.1) enables the formulation of a linear system in Eq. (52). The amplitudes of the fields are determined by inversion of this system. The components of this system are given below:

$$
\mathbf{M} = \begin{bmatrix} \mathbf{M}_1 & \mathbf{M}_2 & \mathbf{M}_3 \\ \mathbf{M}_4 & \mathbf{M}_5 & \mathbf{M}_6 \\ \mathbf{M}_7 & \mathbf{M}_8 & \mathbf{M}_9 \end{bmatrix}, \qquad (A.1)
$$

where M is of size ((J + 14) × (J + 14)), with

$$
\mathbf{M}_1 = \begin{bmatrix}
S_{J_n}(\alpha_e) & S_{Y_n}(\alpha_e) & S_{I_n}(\alpha_e) & S_{K_n}(\alpha_e) & 0 & 0 & 0 & 0 \\
T_{J_n}(\alpha_e) & T_{Y_n}(\alpha_e) & T_{I_n}(\alpha_e) & T_{K_n}(\alpha_e) & 0 & 0 & 0 & 0 \\
J_n(\alpha_m) & Y_n(\alpha_m) & I_n(\alpha_m) & K_n(\alpha_m) & -J_n(\alpha_m) & -Y_n(\alpha_m) & -I_n(\alpha_m) & -K_n(\alpha_m) \\
J'_n(\alpha_m) & Y'_n(\alpha_m) & I'_n(\alpha_m) & K'_n(\alpha_m) & -J'_n(\alpha_m) & -Y'_n(\alpha_m) & -I'_n(\alpha_m) & -K'_n(\alpha_m) \\
S_{J_n}(\alpha_m) & S_{Y_n}(\alpha_m) & S_{I_n}(\alpha_m) & S_{K_n}(\alpha_m) & -S_{J_n}(\alpha_m) & -S_{Y_n}(\alpha_m) & -S_{I_n}(\alpha_m) & -S_{K_n}(\alpha_m) \\
T_{J_n}(\alpha_m) & T_{Y_n}(\alpha_m) & T_{I_n}(\alpha_m) & T_{K_n}(\alpha_m) & -T_{J_n}(\alpha_m) & -T_{Y_n}(\alpha_m) & -T_{I_n}(\alpha_m) & -T_{K_n}(\alpha_m) \\
0 & 0 & 0 & 0 & J_n(\alpha) & Y_n(\alpha) & I_n(\alpha) & K_n(\alpha) \\
0 & 0 & 0 & 0 & k J'_n(\alpha) & k Y'_n(\alpha) & k I'_n(\alpha) & k K'_n(\alpha) \\
0 & 0 & 0 & 0 & S_{J_n}(\alpha) & S_{Y_n}(\alpha) & S_{I_n}(\alpha) & S_{K_n}(\alpha) \\
0 & 0 & 0 & 0 & T_{J_n}(\alpha) & T_{Y_n}(\alpha) & T_{I_n}(\alpha) & T_{K_n}(\alpha)
\end{bmatrix}, \qquad (A.2)
$$

$$
\mathbf{M}_2 = \begin{bmatrix}
0 & \cdots & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & \cdots & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & \cdots & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & \cdots & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & \cdots & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & \cdots & 0 & 0 & 0 & -1 & 0 & 0 \\
0 & \cdots & 0 & 0 & -p_1 & -p_2 & 0 & 0 \\
0 & \cdots & 0 & 0 & -c^2 p_3 & -c^2 p_4 & -c^2 p_5 & 0 \\
0 & \cdots & 0 & 0 & -c^3 p_6 & -c^3 p_7 & -c^3 p_8 & -c^3 p_9
\end{bmatrix}, \qquad (A.3)
$$

$$
\mathbf{M}_5 = \begin{bmatrix}
\Lambda_1(\hat{r}_1) & \Lambda_2(\hat{r}_1) & \Lambda_3(\hat{r}_1) & \Lambda_4(\hat{r}_1) & \Lambda_5(\hat{r}_1) & 0 & 0 \\
0 & \cdots & \cdots & \cdots & \cdots & 0 & 0 \\
0 & 0 & \Lambda_1(\hat{r}_J) & \Lambda_2(\hat{r}_J) & \Lambda_3(\hat{r}_J) & \Lambda_4(\hat{r}_J) & \Lambda_5(\hat{r}_J)
\end{bmatrix}, \qquad (A.4)
$$

$$
\mathbf{M}_8 = \begin{bmatrix}
1 & 0 & 0 & 0 & 0 & \cdots & 0 \\
q_1 & 0 & 0 & 0 & 0 & \cdots & 0 \\
b^2 q_3 & b^2 q_4 & b^2 q_5 & 0 & 0 & \cdots & 0 \\
b^3 q_6 & b^3 q_7 & b^3 q_8 & b^3 q_9 & 0 & \cdots & 0 \\
0 & 0 & 0 & 0 & 0 & \cdots & 0 \\
0 & 0 & 0 & 0 & 0 & \cdots & 0 \\
0 & 0 & 0 & 0 & 0 & \cdots & 0
\end{bmatrix}, \qquad (A.5)
$$

and

$$
\mathbf{M}_9 = \begin{bmatrix}
-J_n(\beta) & -Y_n(\beta) & -I_n(\beta) & -K_n(\beta) & 0 & 0 \\
-k_{II} J'_n(\beta) & -k_{II} Y'_n(\beta) & -k_{II} T'_n(\beta) & -k_{II} K'_n(\beta) & 0 & 0 \\
-\frac{D_{II}}{D_{III}} |_b S_{J_n}(\beta) & -\frac{D_{II}}{D_{III}} |_b S_{Y_n}(\beta) & -\frac{D_{II}}{D_{III}} |_b S_{I_n}(\beta) & -\frac{D_{II}}{D_{III}} |_b S_{K_n}(\beta) & 0 & 0 \\
-\frac{D_{II}}{D_{III}} |_b T_{J_n}(\beta) & -\frac{D_{II}}{D_{III}} |_b T_{Y_n}(\beta) & -\frac{D_{II}}{D_{III}} |_b T_{I_n}(\beta) & -\frac{D_{II}}{D_{III}} |_b T_{K_n}(\beta) & 0 & 0 \\
J'_n(\gamma) & Y'_n(\gamma) & I'_n(\gamma) & K'_n(\gamma) & -J_n(\zeta) & -I_n(\zeta) \\
J'_n(\gamma) & Y'_n(\gamma) & I'_n(\gamma) & K'_n(\gamma) & -\frac{k_I}{k_{II}} J'_n(\zeta) & -\frac{k_I}{k_{II}} I'_n(\zeta) \\
S_{J_n}(\gamma) & S_{Y_n}(\gamma) & S_{I_n}(\gamma) & S_{K_n}(\gamma) & -\frac{D_I}{D_{II}} |_a S_{J_n}(\zeta) & -\frac{D_I}{D_{II}} |_a S_{I_n}(\zeta) \\
T_{J_n}(\gamma) & T_{Y_n}(\gamma) & T_{I_n}(\gamma) & T_{K_n}(\gamma) & -\frac{D_I}{D_{II}} |_a T_{J_n}(\zeta) & -\frac{D_I}{D_{II}} |_a T_{I_n}(\zeta)
\end{bmatrix}. \qquad (A.6)
$$
The matrices M₃, M₄, M₆, and M₇ are matrices of zeros with sizes (10 × 6), ((J - 4) × 8), ((J - 4) × 6) and (8 × 8) respectively. The matrix M₅, resulting from the discretization of the equation of motion, is of size ((J - 4) × J). The unknown and source vectors are respectively written as

$$ \mathbf{X} = \begin{bmatrix} a_n & b_n & c_n & d_n & e_n & f_n & g_n & l_n & W_{n,1} & W_{n,2} & \dots & W_{n,J-1} & W_{n,J} & m_n & n_n & s_n & t_n & u_n & v_n \end{bmatrix}^T, \quad (A.7) $$

$$ \mathbf{S} = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & F & 0 & \dots & 0 \end{bmatrix}^T. \quad (A.8) $$

The coefficients S and T in the matrices are derived from the bending moment (Eq. (29)) and shear force (Eq. (30)) operators, given by:

$$ S_J(\chi) = [n^2(1-\nu) - \chi^2] J_n(\chi) - (1-\nu)\chi J'_n(\chi), \quad (A.9) $$

$$ S_Y(\chi) = [n^2(1-\nu) - \chi^2] Y_n(\chi) - (1-\nu)\chi Y'_n(\chi), \quad (A.10) $$

$$ S_I_n(\chi) = [n^2(1-\nu) + \chi^2] I_n(\chi) - (1-\nu)\chi I'_n(\chi), \quad (A.11) $$

$$ S_{K_n}(\chi) = [n^2(1-\nu) + \chi^2] K_n(\chi) - (1-\nu)\chi K'_n(\chi), \quad (A.12) $$

$$ T_J(\chi) = [n^2(1-\nu)] J_n(\chi) - [\chi(n^2(1-\nu) + \chi^2)] J'_n(\chi), \quad (A.13) $$

$$ T_Y(\chi) = [n^2(1-\nu)] Y_n(\chi) - [\chi(n^2(1-\nu) + \chi^2)] Y'_n(\chi), \quad (A.14) $$

$$ T_I(\chi) = [n^2(1-\nu)] I_n(\chi) - [\chi(n^2(1-\nu) - \chi^2)] I'_n(\chi), \quad (A.15) $$

$$ T_{K_n}(\chi) = [n^2(1-\nu)] K_n(\chi) - [\chi(n^2(1-\nu) - \chi^2)] K'_n(\chi), \quad (A.16) $$

where χ takes the value of the dimensionless wavenumber corresponding to each interface. These values are $k = k_{IV} = k_V$, $\alpha_e = kr_{ext}$, $\alpha_m = kr_m$, $\alpha = kc$, $\beta = k_{II}b$, $\gamma = k_{II}a$, et $\zeta = k_Ia$. The Poisson's ratio $\nu$ depends on the considered zone: $\nu = \nu_{IV} = \nu_V$ for $r > c$ and $\nu = \nu_I = \nu_{II}$ for $r < b$.

The coefficients $q_i$ and $p_i$ ($i = 1...9$) are obtained by a right and left finite-difference schemes as:

$$ q_1 = -\frac{1}{\delta_r}, \quad (A.17) $$

$$ q_2 = +\frac{1}{\delta_r}, \quad (A.18) $$

$$ q_3 = +\frac{1}{\delta_r^2} - \frac{\nu_{III}}{b} \frac{1}{\delta_r} - \frac{\nu_{III} n^2}{b^2}, \quad (A.19) $$

$$ q_4 = -\frac{2}{\delta_r^2} + \frac{\nu_{III}}{b} \frac{1}{\delta_r}, \quad (A.20) $$

$$ q_5 = +\frac{1}{\delta_r^2}, \quad (A.21) $$

$$ q_6 = -\frac{1}{\delta_r^3} + \frac{1}{\delta_r^2} \left( \frac{D'_{III}}{D_{III}} |_b + \frac{1}{b} \right) - \frac{1}{\delta_r} \left( \frac{\nu_{III}}{b} \frac{D'_{III}}{D_{III}} |_b - \frac{(2 - \nu_{III})n^2 + 1}{b^2} \right) \\ - \frac{\nu_{III} n^2 D'_{III}}{b^2 D_{III}} |_b + \frac{(3 - \nu_{III})n^2}{b^3}, \quad (A.22) $$

$$ q_7 = +\frac{3}{\delta_r^3} - \frac{2}{\delta_r^2} \left( \frac{D'_{III}}{D_{III}} |_b + \frac{1}{b} \right) + \frac{1}{\delta_r} \left( \frac{\nu_{III}}{b} \frac{D'_{III}}{D_{III}} |_b - \frac{(2 - \nu_{III})n^2 + 1}{b^2} \right), \quad (A.23) $$

$$ q_8 = -\frac{3}{\delta_r^3} + \frac{1}{\delta_r^2} \left( \frac{D'_{III}}{D_{III}} |_b + \frac{1}{b} \right), \quad (A.24) $$

$$ q_9 = +\frac{1}{\delta_r^3}, \quad (A.25) $$

$$ p_1 = -\frac{1}{\delta_r}, \quad (A.26) $$

$$ p_2 = +\frac{1}{\delta_r}, \quad (A.27) $$

$$ p_3 = +\frac{1}{\delta_r^2}, \quad (A.28) $$
$$p_4 = -\frac{2}{\delta_r^2} - \frac{v_{III}}{c} \frac{1}{\delta_r}, \qquad (A.29)$$

$$p_5 = +\frac{1}{\delta_r^2} + \frac{v_{III}}{c} \frac{1}{\delta_r} - \frac{v_{III} n^2}{c^2}, \qquad (A.30)$$

$$p_6 = -\frac{1}{\delta_r^3}, \qquad (A.31)$$

$$p_7 = +\frac{3}{\delta_r^3} + \frac{1}{\delta_r^2} \left( \frac{D'_{III}}{D_{III}} |_c + \frac{1}{c} \right), \qquad (A.32)$$

$$p_8 = -\frac{3}{\delta_r^3} - \frac{2}{\delta_r^2} \left( \frac{D'_{III}}{D_{III}} |_c + \frac{1}{c} \right) - \frac{1}{\delta_r} \left( \frac{v_{III}}{c} \frac{D'_{III}}{D_{III}} |_c - \frac{(2 - v_{III}) n^2 + 1}{c^2} \right), \qquad (A.33)$$

$$p_9 = +\frac{1}{\delta_r^3} + \frac{1}{\delta_r^2} \left( \frac{D'_{III}}{D_{III}} |_c + \frac{1}{c} \right) + \frac{1}{\delta_r} \left( \frac{v_{III}}{c} \frac{D'_{III}}{D_{III}} |_c - \frac{(2 - v_{III}) n^2 + 1}{c^2} \right) \\ - \frac{v_{III} n^2 D'_{III}}{c^2 D_{III}} |_c + \frac{(3 - v_{III}) n^2}{c^3}. \qquad (A.34)$$

## Data availability

Data will be made available on request.

## References

[1] A. Pelat, F. Gautier, S.C. Conlon, F. Semperlotti, The acoustic black hole: A review of theory and applications, J. Sound Vib. 476 (2020) 115316.

[2] M. Mironov, Propagation of a flexural wave in a plate whose thickness decreases smoothly to zero in a finite interval, Sov. Phys. Acoust. 34 (1988) 318–319.

[3] V.V. Krylov, Localized acoustic modes of a quadratic solid wedge, Mosc. Univ. Phys. Bull. 45 (1990) 65–69.

[4] 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, J. Sound Vib. 330 (11) (2011) 2497–2508.

[5] V. Denis, A. Pelat, F. Gautier, Scattering effects induced by imperfections on an acoustic black hole placed at a structural waveguide termination, J. Sound Vib. 362 (2016) 56–71.

[6] V. Denis, A. Pelat, C. Touzé, F. Gautier, Improvement of the acoustic black hole effect by using energy transfer due to geometric nonlinearity, Int. J. Non-Linear Mech. 94 (2017) 134–145.

[7] H. Li, P. O'Donoughue, F. Masson, A. Pelat, F. Gautier, C. Touzé, Broadband shock vibration absorber based on vibro-impacts and acoustic black hole effect, Int. J. Non-Linear Mech. 159 (2024) 104620.

[8] D.J. O'Boy, V.V. Krylov, Damping of flexural vibrations in circular plates with tapered central holes, J. Sound Vib. 330 (10) (2011) 2220–2236, Dynamics of Vibro-Impact Systems.

[9] O. Aklouche, A. Pelat, S. Maugeais, F. Gautier, Scattering of flexural waves by a pit of quadratic profile inserted in an infinite thin plate, J. Sound Vib. 375 (2016) 38–52.

[10] J. Deng, Z. Ling, P. Zeng, Y. Zuo, O. Guasch, Passive constrained viscoelastic layers to improve the efficiency of truncated acoustic black holes in beams, Mech. Syst. Signal Process. 118 (2018) 461–476.

[11] G. Raybaud, M. Quisse, J. Leng, A. Pelat, J.P. Groby, V. Romero-García, R. Picó, F. Gautier, Control of bending wave reflection at beam terminations by thermally tunable subwavelength resonators, J. Sound Vib. 530 (2022) 116918.

[12] W. Huang, H. Zhang, D. Inman, C. Cesnik, H. Ji, Low reflection effect by 3D printed functionally graded acoustic black holes, J. Sound Vib. 450 (2019).

[13] K. Mizukami, M. Shiratori, K. Ogi, Fiber-steered acoustic black hole beam with low cut-on frequency and high stiffness, J. Sound Vib. (2024) 118396.

[14] A.C. Nilsson, B. Liu, Prediction of some vibro-acoustic properties of sandwich plates with honeycomb and foam cores, J. Acoust. Soc. Am. 144 (3) (2018) 1600–1614.

[15] A. Nilsson, B. Liu, Vibro-Acoustics, Vol. 2, Springer Berlin Heidelberg, Berlin, Heidelberg, 2016.

[16] E. Nilsson, A.C. Nilsson, Prediction and measurement of some dynamic properties of sandwich structures with honeycomb and foam cores, J. Sound Vib. 251 (3) (2002) 409–430.

[17] K. Ege, N. Roozen, Q. Leclère, R. Rinaldi, Assessment of the apparent bending stiffness and damping of multilayer plates; modelling and experiment, J. Sound Vib. 426 (2018) 129–149.

[18] N. Peake, S.V. Sorokin, On the behavior of fluid-loaded sandwich panels with mean flow, J. Sound Vib. 242 (4) (2001) 597–617.

[19] E.P. Bowyer, V.V. Krylov, Experimental investigation of damping flexural vibrations in glass fibre composite plates containing one- and two-dimensional acoustic black holes, Compos. Struct. 107 (2014) 406–415.

[20] T. Zhou, L. Tang, H. Ji, J. Qiu, L. Cheng, Dynamic and static properties of double-layered compound acoustic black hole structures, Int. J. Appl. Mech. 9 (5).

[21] M. Geradin, D.J. Rixen, Mechanical Vibrations: Theory and Application to Structural Dynamics, John Wiley & Sons, 2014.

[22] A.W. Leissa, Vibration of plates, 1969, NASA Headquarters Number: NASA-SP-160.

[23] M.N. Ichchou, O. Bareille, J. Berthaut, Identification of effective sandwich structural properties via an inverse wave approach, Eng. Struct. 30 (10).

[24] C. Pézerat, J.-L. Guyader, Two inverse methods for localization of external sources exciting a beam, Acta Acust. 3 (1995) 1–10.

[25] C. Pézerat, J.-L. Guyader, Force Analysis Technique: Reconstruction of Force Distributions on Plates, Acta Acust. United Ac. 86 (2) (2000) 322–332.

[26] F. Ablitzer, C. Pézerat, B. Lascoup, J. Brocail, Identification of the flexural stiffness parameters of an orthotropic plate from the local dynamic equilibrium without a priori knowledge of the principal directions, J. Sound Vib. 404 (2017) 31–46.

[27] A.N. Norris, C. Vemula, Scattering of flexural waves on thin plates, J. Sound Vib. 181 (1) (1995) 115–125.

[28] K.F. Graff, Wave Motion in Elastic Solids, Dover Publications, New York, 1991, OCLC 829198880.

[29] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions With Formulas, Graphs, and Mathematical Tables, vol. 55, 1972.

[30] D. Ross, E. Ungar, E.K. Kerwin Jr., Damping of plate flexural vibrations by means of viscoelastic laminae, Struct. Damping (1959) 49–97.