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

HAL
open science

# Energy Localization through Locally Resonant Materials

Marco Moscatelli, Claudia Comi, Jean-Jacques Marigo

► To cite this version:

Marco Moscatelli, Claudia Comi, Jean-Jacques Marigo. Energy Localization through Locally Resonant Materials. Materials, 2020, 13 (13), pp.3016. 10.3390/ma13133016. [https://doi.org/10.3390/ma13133016]

HAL Id: hal-03591208

https://hal.science/hal-03591208v1

Submitted on 11 Jan 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Energy Localization through Locally Resonant Materials

Marco Moscatelli<sup>1,2 ID</sup>, Claudia Comi<sup>1,* ID</sup> and Jean-Jacques Marigo<sup>2 ID</sup>

Received: 14 May 2020; Accepted: 1 July 2020; Published: 6 July 2020

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

check for
updates

**Abstract:** Among the attractive properties of metamaterials, the capability of focusing and localizing waves has recently attracted research interest to establish novel energy harvester configurations. In the same frame, in this work, we develop and optimize a system for concentrating mechanical energy carried by elastic anti-plane waves. The system, resembling a Fabry-Pérot interferometer, has two barriers composed of Locally Resonant Materials (LRMs) and separated by a homogeneous internal cavity. The attenuation properties of the LRMs allow for the localization of waves propagating at particular frequencies. With proper assumptions on the specific ternary LRMs, the separation of scales (between the considered wave lengths and the characteristic dimension of the employed unit cells) enables the use of a two-scale asymptotic technique for computing the effective behavior of the employed LRMs. This leads to a complete analytic description of the motion of the system. Here we report the results achieved by optimizing the geometry of the system for obtaining a maximum focusing of the incoming mechanical energy. The analytic results are then validated through numerical simulations.

**Keywords:** locally resonant material; energy harvesting; homogenization

## 1. Introduction

The conversion of vibrations into low-power electricity has been extensively explored since decades in view of the wide range of applications of the harvesters in distributed wireless sensors for structural health monitoring and for security systems, in embedded and implanted medical sensors and in recharging of batteries in various systems, see, e.g., the reviews [1,2]. To obtain an efficient energy harvester, one has to focus the mechanical energy produced by external vibrations in a given region of the system, where it can be then converted in electricity by means, e.g., of piezoelectric materials.

The use of structured materials for Energy Harvesters (EH) has more recently attracted the attention of the researchers and a number of possible configurations have been proposed [3–7]. Localized defects in phononic crystals [8] and in micro-structured plates [9] have been exploited to focus vibration energy in a small region where a piezoelectric material converts the mechanical energy into the electric one. However, good performance is obtained only for frequencies close to the defect eigen-frequency which is related with the typical size of the unit cell of the lattice. Even though proper methods can be used to optimize the geometry of the cell [10], for small devices, this frequency is very high when compared with the frequency of ambient vibrations in common applications. Localized defects in tensegrity materials have also been studied in [11].

Locally resonant metastructures, which are structures that comprise Locally Resonant Material (LRM) components, enable band gap formation at wavelengths much longer than the lattice size

<footnote n="1">¹ Department of Civil and Environmental Engineering, Politecnico di Milano, Piazza Leonardo da Vinci 32, 20133 Milan, Italy; marco.moscatelli@polimi.it</footnote>

<footnote n="2">² Laboratoire de Mécanique des Solides, École Polytechnique, Route de Saclay, 91120 Palaiseau, France; jean-jacques.marigo@polytechnique.edu</footnote>

<footnote>* Correspondence: claudia.comi@polimi.it</footnote>
(see, e.g., [12–16]), therefore they can represent good candidates for efficient, small EH. In [17] the localized deformation pattern achieved in the frequency neighborhood of the band gaps was employed in order to harvest a certain amount of the kinetic energy available in the oscillating members of the lattice. Recently [18] proposed a prototype of locally resonant energy harvesting-metastructure composed of a primary beam with several small secondary cantilever beams with tip masses acting as mechanical resonators. The design of innovative metasurfaces, with graded resonators, that trap waves was proposed and developed in [5] for enhanced piezoelectric energy harvesting.

In the present work we explore the possibility to couple the advantage of the LRM mechanism to create band gaps at low frequencies with the energy localization mechanism in local defects of regular lattices to design and optimize a Resonant Energy Harvester (REH). An initial attempt in this direction was presented in [19], considering binary LRMs. The ternary LRM here considered, endowed with two-dimensional periodicity, is constituted by a stiff matrix with almost rigid cylindrical inclusions coated by a compliant material. The two-scale homogenization approach, proposed in [20] for high-contrast binary composite materials in the long wavelength regime and developed by many authors ([21–23]), provides a powerful tool to define equivalent material properties. The authors recently studied through homogenization the spectral properties and the band gaps of binary and ternary LRM in [16,24], respectively. The REH here proposed is constituted by two barriers made of the above LRM and a cavity, made of the matrix alone, which represents a line defect in the two dimensional regular lattice. The use of barriers made of LMR and of phononic crystals was also recently proposed to obtain an acoustic diode [25]. To optimize the REH, the behavior of the LRM is here characterized by its homogenized properties and the dynamic response is analytically obtained. The conditions to obtain the maximal mechanical energy inside the cavity are explicitly given and the role of the different configurations of the unit cell of the metamaterial are highlighted.

Beside extending the analytic results obtained in [19] for a REH with binary LRM to the case of a REH with three components LRM, the paper provides a new numerical validation, introduces the concept of indexes of concentration and develops maps of these indexes at varying frequency and geometry of the REH. The maps can be of practical use in the design of a real REH that, with a fixed geometry, should work in a wide range of frequencies.

The paper is organized as follows. The proposed configuration of the REH and the dynamic problem of shear wave propagation is set in Section 2. In Section 3 the energy in each part of the system is computed and the optimal width of the cavity to maximize the energy trapped in the cavity is obtained. Section 4 provides a comparison of the analytic results with the numerical results obtained for a special choice of the LRM, thus validating the analytic results. These latter allows to perform a parametric study of the system. Some conclusions are given in Section 5. The results of the asymptotic homogenization, obtained in [24], which are used in the paper, are recalled in the Appendix A.

Notation, vectors and tensors are indicated by bold face letters. Complex numbers are denoted with sans serif letters, while italic is utilized for the Real numbers. The complex conjugate of a complex number is indicated by a superposed bar. The angular frequency $\omega$ is often called just frequency. The notation $\langle \cdot \rangle$ denotes the average over the period.

## 2. The Resonant Cavity in Metamaterials

### 2.1. *Problem Formulation*

We consider the system represented in Figure 1. It is composed of five domains: parts $\Omega_1$, $\Omega_3$ and $\Omega_5$ are constituted by an isotropic homogeneous medium (defined as $\Omega_h$ in the following), whereas parts $\Omega_2$ and $\Omega_4$, which will be denoted as “barriers” ($\Omega_b$), are made of the same LRM. This system, which is similar to a Fabry-Pérot interferometer in optics [26], will be referred to as Resonant Energy Harvester (REH). A similar configuration has been already studied by the
authors in [19], although now the two barriers are made of a ternary LRM instead of the binary heterogeneous material previously employed.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Top-view of the studied system. The wave shown in the figure represents an out-of-plane wave propagating throughout all the five domains. The dots denote the fact that the system must be extended along the $x_2$ direction.</figcaption></figure>

The LRM is composed of a two dimensional periodic repetition of the square unit cell represented in Figure 2, with a stiff matrix $Y_m$ containing an almost rigid circular fiber $Y_f$, coated by a very compliant layer $Y_c$. The developed system has a thickness in the out-of-plane direction and along the axis $x_2$ which is much larger than the characteristic size $L$ of the unit cell; parts $\Omega_1$ and $\Omega_5$ are infinitely extended towards $-\infty$ and $+\infty$ along the axis $x_1$.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Lattice of the LRM and zoom over a unit cell composing it. $R_f$ and $R_c$ are respectively the fiber and coating external radii.</figcaption></figure>

We analyze the propagation of anti-plane waves along the $x_1$ direction, allowing for the decoupling from the in-plane wave propagation problem. In the harmonic regime, this results in the following two dimensional Helmholtz scalar equation:

$$ \begin{cases} \mathrm{div}\sigma + \rho\omega^2 u = 0, & \sigma = \mu \mathrm{grad} u \\ \text{with } u(x,t) \text{ and } \sigma \cdot n \text{ continuous at each interface} \end{cases} \qquad (1) $$

where the elastic displacement $u(x,t)$ has a time dependence $e^{i\omega t}$, which from now on will be omitted and $\sigma$ collects the non vanishing stress components $\sigma_{31}$ and $\sigma_{32}$. The mass density $\rho$ and the shear elastic modulus $\mu$ are spatially varying parameters, depending on the position $x$ in the system, as specified in the following:

$$ \mu = \begin{cases} \mu_m & \text{in } Y_m \\ \mu_f & \text{in } Y_f, \\ \mu_c & \text{in } Y_c \end{cases}, \qquad \rho = \begin{cases} \rho_m & \text{in } Y_m \\ \rho_f & \text{in } Y_f \\ \rho_c & \text{in } Y_c \end{cases} \qquad (2) $$

Denoting by $k_m = \omega\sqrt{\rho_m/\mu_m}$ the wave number in the matrix, we limit our study by considering a low frequency regime such that the dimensionless parameter $\epsilon = k_mL$ be very small ($\epsilon \ll 1$). This assumption can be equivalently set by requiring that the characteristic size $L$ of the cell be much smaller than the wave length, in the matrix, of the considered wave. This hypothesis enables the application of a two-scale asymptotic homogenization technique for the description of the motion
of the barriers. By assuming $ρ_m/ρ_c = O(1)$ and the ratio $μ_c/μ_m$ of order $O(ε^2)$ the effective mass density $ρ_{eff}$ and the effective shear modulus $μ_{eff}$ of the ternary LRM were computed in [24] through homogenization. In the Appendix A we summarize the main results derived in [24], which are used for treating the problem here analyzed. The reader is thus addressed to that paper for a more complete derivation.

## 2.2. Solution of the Homogenized Problem

By virtue of the homogenization, the REH can be treated as an equivalent system composed of five homogeneous and isotropic parts. Since we are interested in the propagation of waves whose wave front are perpendicular to the $x_1$ direction, there is no dependence on the $x_2$ direction, $σ_{32} = 0$ and the problem becomes one-dimensional, obtaining the simplified system sketched in Figure 3 (from now on thick lines will represent a LRM).

<figure><img src="image_6.png" /><figcaption>**Figure 3.** One-dimensional counterpart of the two-dimensional homogenized problem. The system is simplified by fixing a position along the $x_2$ axis and by considering $x_1 = x$. The dashed ends are used to indicate that the domain is infinitely extended toward $-∞$ and $+∞$.</figcaption></figure>

To summarize, we report here below the homogenized problem governing the motion inside each of the five parts composing our system:

$$ \left\{ \begin{array}{ll} \text{div}\sigma(x) + \rho_m \omega^2 U(x) = 0, & \sigma = \mu_m \partial U(x)/\partial x \quad \text{in } \Omega_h \\ \text{div}\sigma(x) + \rho_{eff}(\omega) \omega^2 U(x) = 0, & \sigma = \mu_{eff} \partial U(x)/\partial x \quad \text{in } \Omega_b \\ \text{with } U(x) \text{ and } \sigma(x) \text{ continuous at each interface} \end{array} \right. \quad (3) $$

where $U(x)$ is the first term in the expansion of the displacement field in the matrix (denoted by $U^0(x)$ in the Appendix A) and $σ = σ_{31}$ denotes the only non-zero stress component.

As specified in the Appendix A, the sign of the effective mass density $ρ_{eff}$ depends whether the frequency $ω$ is inside or outside a band gap; this means that within the barriers $Ω_2$ and $Ω_4$ the second of the Equation (3) has a general solution whose form varies depending on the frequency, as follows:

$$ \left\{ \begin{array}{ll} U_j(x) = A_j e^{-isx} + B_j e^{isx} & \text{when } \rho_{eff} \ge 0 \\ U_j(x) = A_j \cosh sx + B_j \sinh sx & \text{when } \rho_{eff} < 0 \end{array} \right. \quad (4) $$

where $s = ω \sqrt{|ρ_{eff}|/μ_{eff}}$, $A_j$ and $B_j$ (with $j = 2, 4$) are complex integration constants.

The general solution of the motion problem inside the regions $Ω_h$ is instead always given by:

$$ U_j(x) = A_j e^{-ikx} + B_j e^{ikx} \quad (5) $$

with $j = 1,3$ or $5$ and $k = k_m = ω \sqrt{ρ_m/μ_m}$.

Considering a propagating input wave, the integration constants are obtained imposing the continuity of the displacement and of the stress between the various regions and the conditions at infinite. In particular, when considering an incoming wave which travels from the left toward the right (see Figure 3) of amplitude 1, the displacement in the different parts $Ω_i$ reads (see [19] for details):
$$
\left\{
\begin{aligned}
U_1(x) &= e^{-ik(x+d+l)} + R e^{ik(x+d+l)} \\
U_2(x) &= (1+R) \cosh s(x+d+l) + i a (R-1) \sinh s(x+d+l) \\
U_3(x) &= \frac{T}{2} \beta e^{ik(x-d)} + \frac{T}{2} \alpha e^{-ik(x-d)} \\
U_3(x) &= T \cosh s(x-d-l) - i a T \sinh s(x-d-l) \\
U_5(x) &= T e^{-ik(x-d-l)}
\end{aligned}
\right.
\qquad (6)
$$

where R and T are respectively the amplitudes of the reflected and transmitted wave:

$$
T = \frac{4}{\alpha^2 e^{2ikd} - \beta^2 e^{-2ikd}} \qquad (7)
$$

$$
R = \frac{\alpha\beta e^{2ikd} - \bar{\alpha}\bar{\beta} e^{-2ikd}}{\alpha^2 e^{2ikd} - \beta^2 e^{-2ikd}} \quad (8)
$$

with

$$
a = \frac{\mu_m k}{\mu_{eff} s'}, \quad \alpha = 2 \cosh sl + i \left(a - \frac{1}{a}\right) \sinh sl, \quad \beta = i \left(a + \frac{1}{a}\right) \sinh sl \quad (9)
$$

3. Energy Localization in the Cavity

3.1. Energy in the Homogeneous Parts

We are now interested in finding the mechanical energy density, averaged over a period, along the
whole REH ⟨e(x, t)⟩ and the total energy in parts Ω₂, Ω₃ and Ω₄ again averaged over a time period,
⟨E(t)⟩.

The mechanical energy in regions Ωⱼ with j = 1, 3 or 5 of homogeneous material, given by the
sum of the potential energy p and the kinetic energy c, reads:

$$
e_j(x,t) = p_j + c_j = \frac{1}{2} \mu_m u_j^2 + \frac{1}{2} \rho_m \dot{u}_j^2 \qquad (10)
$$

where the dependency over time is again accounted for, a superimposed dot denotes the time
derivative and a prime denotes the space derivative; the displacement field $u_j$ is given by $u_j(x,t) =$
$Re\{U_j(x) e^{i\omega t}\}$.

Averaging Equation (10) over time gives:

$$
\langle e_j \rangle = \frac{\mu_m k^2}{2} [|A_j|^2 + |B_j|^2] \quad (11)
$$

Notice that ⟨e_j⟩ is independent from the position x. Finally, the total mechanical energy ⟨E(t)⟩ inside part Ω_3 reads:

$$
\langle E_3 \rangle = 2d \langle e_3 \rangle \tag{12}
$$

While the energy inside the regions Ωh can be easily obtained, a particular treatment must be
devoted to the calculation of the energy inside the two barriers.

3.2. Energy Inside the LRM: Homogenization Approach

Each cell is composed of a matrix, a coating layer and a fiber considered as rigid, hence the
mechanical energy density of one cell reads:

$$
e_j(x,t) = p_j^m + c_j^m + p_j^c + c_j^c + c_j^f \qquad (13)
$$
The homogenization technique enables to find the displacement field and the stress over a cell composing the metamaterial. It should be noted that, while at the macro-scale the problem is one-dimensional, so that only the variable $x_1 = x$ must be retained, at the micro-scale (inside each cell) the problem remains 2D and the coordinates $\mathbf{y} = y_1, y_2$ or the polar coordinates $r, \theta$ must be retained. Similarly, inside the cell the stress $\sigma_{32}$ is non-zero and a stress vector should be considered. At the leading order, one has:

$$ u(x,t) \approx u^0 = \begin{cases} U(x) e^{i\omega t} & \text{in } Y_m \\ U(x)\eta(r) e^{i\omega t} & \text{in } Y_c \\ U(x)\eta(R_f) e^{i\omega t} & \text{in } Y_f \end{cases} \qquad (14) $$

and

$$ \sigma(x,t) \approx \sigma^0 = \begin{cases} \mu_m (\text{grad}_y \chi + \mathbf{e}_1) \frac{\partial U(x)}{\partial x} e^{i\omega t} & \text{in } Y_m \\ \mathbf{0} & \text{in } Y_c \end{cases} \qquad (15) $$

with $\eta(r)$ and $\chi$ defined in the Appendix A.

Considering a frequency inside the band gap of the metamaterial, in regions $\Omega_j, j = 2, 4$, one has $U(x) = U_j(x), j = 2, 4$, with $U_j$ given by Equation (6). Denoting the real part of the displacement in the matrix by $u_j(x,t) = \text{Re}\{U_j(x) e^{i\omega t}\}$ the contributions to the energy density in Equation (13) read

$$ p_j^m = \frac{1}{2} \mu_{eff} u_j^2 \qquad (16) $$

$$ c_j^m = \frac{|Y_m| \rho_m}{|Y| 2} \dot{u}_j^2 \qquad (17) $$

$$ p_j^c = \frac{\mu_c}{2} \int_{Y_c} \text{grad}_y \eta \cdot \text{grad}_y \eta \, dy \frac{u_j^2}{|Y|} \qquad (18) $$

$$ c_j^c = \frac{\rho_c}{2} \int_{Y_c} \eta^2 dy \frac{\dot{u}_j^2}{|Y|} \qquad (19) $$

$$ c_j^f = \frac{\rho_f}{2} \eta(R_f)^2 \frac{|Y_f|}{|Y|} \dot{u}_j^2 \qquad (20) $$

From Equations (16) to (20), after some manipulation, the energy $\langle e_j \rangle$, with $j = 2$ or 4, can be written as:

$$ \langle e_2 \rangle = \frac{1}{4} \left\{ \mu_{eff} s^2 \left| (1+R) \sinh sl + a i (R-1) \cosh sl \right|^2 + \omega^2 \gamma \left| (1+R) \cosh sl + a i (R-1) \sinh sl \right|^2 \right\} \qquad (21) $$

$$ \langle e_4 \rangle = \frac{1}{4} \left\{ \mu_{eff} s^2 \left| T \sinh sl - a i T \cosh sl \right|^2 + \omega^2 \gamma \left| T \cosh sl - a i T \sinh sl \right|^2 \right\} \qquad (22) $$

where $\gamma$ is given by:

$$ \gamma = \rho_{eff} + \frac{2}{|Y|} \left[ \rho_c \int_{Y_c} \eta^2 - \eta dy + \rho_f |Y_f| \eta(R_f) (\eta(R_f) - 1) \right] \qquad (23) $$
The averaged mechanical energy density for the barriers is thus dependent on the position x along them. Integrating over their thickness l, one can derive the averaged total mechanical energies $\langle E_2 \rangle$ and $\langle E_4 \rangle$ as:

$$ \langle E_2 \rangle = \frac{1}{4} \left\{ \left[ \frac{l}{2}(a^2 - 1)(\mu_{eff}s^2 - \omega^2\gamma) + \frac{\sinh 2sl}{4s}(a^2 + 1)(\mu_{eff}s^2 + \omega^2\gamma) \right] + \right. \\ \left. - 2 Re\{R\} \left[ \frac{l}{2}(a^2 + 1)(\mu_{eff}s^2 - \omega^2\gamma) + \frac{\sinh 2sl}{4s}(a^2 - 1)(\mu_{eff}s^2 + \omega^2\gamma) \right] + \right. \\ \left. - a Im\{R\}(1 - \cosh 2sl)(\mu_{eff}s^2 + \omega^2\gamma) \right\} \quad (24) $$

$$ \langle E_4 \rangle = \frac{1}{4} |T|^2 \left[ \frac{l}{2}(a^2 - 1)(\mu_{eff}s^2 - \omega^2\gamma) + \frac{\sinh 2sl}{4s}(a^2 + 1)(\mu_{eff}s^2 + \omega^2\gamma) \right] \quad (25) $$

### 3.3. Optimal Cavity Width

Our objective is now to exploit the attenuation capabilities of the employed LRMs for obtaining a concentration of mechanical energy inside the cavity. Let us focus our attention on frequencies inside a band gap for the analyzed LRM.

Fixing the materials used for the system, the width l of the barriers and the frequency of the propagating wave, the modulus of the coefficient T given in Equation (7) is maximized and becomes equal to 1 for a discrete set of optimal cavity widths $\tilde{d}_n, n \in \mathbb{N}^*$ and, at the same time, $|R|$ is null. The optimal cavity widths $\tilde{d}_n$ can be found from the following relation (see [19] for the details):

$$ \tan 4kd = \frac{2 [2 \cosh sl] [(a - \frac{1}{a}) \sinh sl]}{[(a - \frac{1}{a}) \sinh sl]^2 - [2 \cosh sl]^2} \quad (26) $$

One should notice that Equation (26) gives also the set of d that minimize T.

Using the expression of the displacement $U_3$ into Equation (11), the total mechanical energy Equation (12) inside part $\Omega_3$ reads:

$$ \langle E_3 \rangle = d \frac{\mu_m k^2}{4} |T| [|\alpha|^2 + |\beta|^2] \quad (27) $$

From Equation (27), as both $\alpha$ and $\beta$ are independent from d, the mechanical energy inside part $\Omega_3$ is maximized whenever T = 1, i.e., when T is maximum.

When $d = \tilde{d}_n$, the total energy of the regions $\Omega_2, \Omega_3$ and $\Omega_4$ reads:

$$ \begin{cases} \langle E_2 \rangle = \langle E_4 \rangle = \frac{1}{4} \left[ \frac{l}{2}(a^2 - 1)(\mu_{eff}s^2 - \omega^2\gamma) + \frac{\sinh 2sl}{4s}(a^2 + 1)(\mu_{eff}s^2 + \omega^2\gamma) \right] \\ \\ \langle E_3 \rangle = d \frac{\mu_m k^2}{4} [|\alpha|^2 + |\beta|^2] \end{cases} \quad (28) $$

### 3.4. Limit Case of the Well

By extending the barriers toward $-\infty$ and $+\infty$ the homogeneous parts $\Omega_1$ and $\Omega_5$ disappear and the above system tends to a well, as sketched in Figure 4.

<figure><img src="image_7.png" /><figcaption>Figure 4. Sketch of the studied limit case. The barriers are infinitely extended toward $-\infty$ and $+\infty$.</figcaption></figure>
The difference with respect to the REH case, expressed by Equation (3), is that now the two conditions at $x \to -\infty$ and $x \to +\infty$ are imposed to the metamaterial. Since no energy can be accumulated at $+\infty$ and $-\infty$, the motion inside parts $\Omega_2$ and $\Omega_4$ of the well must respect the following statements:

$$ U_2(x \to -\infty) = 0 \quad U_4(x \to +\infty) = 0 \qquad (29) $$

By applying the above conditions to the general solution of Equation (3), two of the six unknown integration constants can be given as functions of the remaining four. Furthermore, the continuity of displacements and stresses at each interface leads to a homogeneous system of four equations. This system admits a non-trivial solution only for a set of $kd$ values, which can be found from the following relation:

$$ \tan 2kd = \frac{2a}{a^2 - 1} \qquad (30) $$

By fixing for instance the amplitude $W$ of the left-wards traveling wave inside the well, the remaining three amplitudes are found and the motion inside each part is finally defined.

The only solutions are now modes for which the frequency is a resonant frequency given by Equation (30). Note that condition Equation (30) can also be obtained directly from Equation (26) by taking the limit for $l \to \infty$. The total mechanical energy stored by each part can be expressed as a function of $|W|$ and reads:

$$ \begin{aligned} \langle E_3 \rangle &= 2\theta |W|^2 \\ \langle E_2 \rangle &= \langle E_4 \rangle = \delta |W|^2 \end{aligned} \qquad (31) $$

with

$$ \theta = \mu_m k^2 d, \quad \delta = \frac{a^2}{2s(1+a^2)} (\mu_{eff}s^2 + \omega^2\gamma) \qquad (32) $$

Fixing the amount of total energy of the system, one can find the coefficient $|W|$ and, thus, the motion inside each part.

## 4. Results

The results of the previous sections are valid provided several hypotheses on the materials composing the REH are fulfilled. For clarity, we summarize them here below:

* Parts $\Omega_1$, $\Omega_3$ and $\Omega_5$ are constituted by the same material utilized for the matrix composing the LRM;

* The coating layer of the fibers must be very compliant with respect to the matrix;

* The fibers must be very stiff so that they can be treated as rigid in the homogenization procedure.

Since we are now interested in considering a real possible application, the material properties shall be fixed by respecting the three conditions just specified. We choose to use the following combination: an epoxy material is employed for the matrix, with lead inclusions coated by a layer of rubber. The material properties are specified in Table 1.

<table><caption>**Table 1.** Material properties.</caption><thead><tr><th>Material</th><th>E [MPa]</th><th>ν</th><th>ρ [Kg/m³]</th></tr></thead><tbody><tr><td>Epoxy</td><td>3600</td><td>0.370</td><td>1180</td></tr><tr><td>Rubber</td><td>0.118</td><td>0.469</td><td>1300</td></tr><tr><td>Lead</td><td>14000</td><td>0.420</td><td>11340</td></tr></tbody></table>
We consider a lattice with square unit cells with side L = 1 mm, the radius of the coated circular inclusion is $R_c = 0.415L$, while the radius of the internal part is $R_f = 0.7R_c$. A scaling of the cell would simply scale the effective mass density when plotted with respect to the frequency, without changing the band gap structure; considering different filling fractions ($\pi R_c^2/L^2$) or thicknesses of the coating, ($R_c - R_f$) would instead modify the dynamic behavior and the band structure (see [24] for other examples).

## 4.1. Energy Localization: Analytic and Numerical Results

To check the validity of the analytic expressions derived in the previous sections, we have carried out some numerical analyses by using the commercial software COMSOL Multiphysics 5.4.

First of all, for the chosen material constants, from the homogenization technique we expect the presence of band gaps, as one can see from Figure 5a, where we have plotted the $\rho_{eff}$ for a range of frequencies $\omega$ from 0 to 25 kHz. The effective mass density becomes indeed negative between 2.4 and 5.7 kHz and between 22.5 and 24.1 kHz (within the range of frequencies shown by the plot). These intervals correspond to the band gaps obtained numerically by applying Bloch-Floquet’s periodic conditions at the unit cell boundaries of the LRM, as shown in Figure 5b. For the dispersion analysis, the cell is meshed with 2D triangular quadratic Lagrange finite elements and the formal analogy with the acoustic problem is exploited. The real stiffness of lead, see Table 1, is used, nonetheless the agreement with the analytic prediction, with the rigid inclusion assumption, is very good. One should remark that in the second band gap there are two flat modes which correspond to local resonances of the coated inclusion endowed with almost zero displacements in the matrix as already shown in [16].

<figure><img src="image_8.png" /><figcaption>Figure 5. (a) Effective mass density vs. frequency, shaded areas correspond to negative effective mass, (b) dispersion plot with shaded band gaps obtained from the numerical Bloch-Floquet’s analysis; Irreducible Brillouin Zone and path followed for the numerical analysis is reported in the inset.</figcaption></figure>

Fixing to 40 the number of cells employed for each of the two barriers (hence fixing the thickness $l$) and optimizing the REH to work at a mid-band gap frequency $\tilde{\omega} = 4.07$ kHz, Equation (26) gives a set of $\tilde{d}_n$, each of them giving a value of $\langle e_3 \rangle$ 23 times bigger than the energy $\langle e_{in} \rangle$ carried by an incoming wave of unit amplitude. The energy of the incoming wave is obtained from Equation (11) and for a REH with the optimal cavity reads:

$$ \langle e_{in} \rangle = \frac{\mu_m k^2}{2} \qquad (33) $$

The dynamic behaviour of this REH and the energy densities in the different domains were also numerically computed. In the numerical analyses, the actual two-dimensional system sketched in Figure 1 is considered. In the $x_2$ direction, only one row of the metamaterial is discretized and symmetry boundary conditions are imposed to simulate the ideal case. In the $x_1$ direction perfectly matched layers (PML) are added to consider the extension of the matrix towards infinity. 2D triangular quadratic elements are used to mesh the barriers, while 2D rectangular quadratic Lagrange finite
elements are employed for the remaining homogeneous parts. Use is made again of the formal analogy with the acoustic problem, hence a background pressure field, analogous to the incoming out of plane wave, is imposed inside part $\Omega_1$ and the energy density is integrated in $x_2$ and in time over a period to be compared with the analytic results.

The plot in Figure 6 shows $\langle e_j \rangle / \langle e_{in} \rangle$ (with $j$ from 1 to 5) along the whole system $\cup_{j=1}^5 \Omega_j$ for $\tilde{d}_1 = 19.9$ mm. The orange lines represent the energy density computed by the analytic Equations (11) and (21), while the results of the numerical analysis are shown in blue. The oscillations of the numerical response are due to the intrinsic heterogeneity of the LRM composing the barriers, the analytic results, based on an homogenized material, give a mean value in good agreement.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Ratio of the averaged mechanical energy density along the whole system, with respect to the incoming one. Orange: analytic results and blue: numerical results. Each part composing the REH is separated by the vertical dashed lines.</figcaption></figure>

In Figure 7a, we report the transmission coefficient modulus $|T|$ vs. frequency, from a transmission analysis of our REH. The results coming from the asymptotic technique agrees very well with those of the real case numerically studied. The expected peak of transmission at the frequency $\tilde{\omega}$ inside the band gap is well captured both by the analytic and numerical results

For showing how the presence of a cavity modifies this behavior, we have plotted in Figure 7b the results coming from a transmission analysis in absence of the cavity, i.e., for a simple layer of LRM with a thickness equal to $2l$ (obtained by attaching together the two barriers and getting rid of the cavity in the middle). Comparing these results with those in Figure 7a, it is clear how the peak inside the band gap disappears when no cavity is present.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Transmission coefficient modulus $|T|$ as a function of the frequency $\omega$ in presence of the cavity (a) and without any cavity (b). The vertical red lines define the band gap. Orange: analytic results and blue: numerical results of analyses in the frequency domain, considering a sweep of frequencies from 0 to 20 KHz.</figcaption></figure>
## 4.2. Towards the Optimization of the Harvester: Parametric Study

The results of the previous subsection refer to a particular case, with all the parameters involved in the problem fixed. In this subsection we consider the effect of different possible configurations of the system and, to compare them, we introduce two indexes as a “measure” of the harvesting capabilities of the REH.

By keeping for the unit cell the same geometric dimensions employed in the previous subsection, there are three parameters which can be left free to vary, namely the barrier width $l$, the optimal cavity width $\tilde{d}_n$ and the frequency $\omega$. To compare different configurations we introduce the following two indexes:

$$IC = \frac{\langle E_3 \rangle}{\langle E_2 \rangle + \langle E_3 \rangle + \langle E_4 \rangle} \qquad (34)$$

$$AIC = \frac{\langle E_3 \rangle / (\tilde{d}_n)}{\langle (E_2 + \langle E_3 \rangle + \langle E_4 \rangle) \rangle / (\tilde{d}_n + l)}$$

where IC stands for “index of concentration” and AIC for “averaged index of concentration”. The former represents a measure of the concentration level of energy reached by the system, whereas the latter gives a ratio between values of energy averaged over the dimensions of the corresponding REH components.

As stated previously, Equation (26) gives a set of optimal cavity widths $2\tilde{d}_n, n \in \mathbb{N}^*$, that maximize the focused energy for a given frequency inside the band gap. In the contours of Figure 8, we show the variation of the two smallest elements of this set $\tilde{d}_1$ and $\tilde{d}_2$ with respect to the frequency $\omega$ inside the band gap and the width $l$ of the barriers. One can observe that the optimal cavity width is almost independent from the width of the barriers $l$ while it strongly depends on the frequency. This dependence is slightly attenuated when considering thin barriers. This is important for the real design of an harvester, with a fixed cavity width, which should work in a wide range of frequencies.

Keeping the frequency inside the band gap and the barrier width as variables, the two indexes IC and AIC are plotted in Figure 9 using the first two optimal cavity widths $2\tilde{d}_1$ and $2\tilde{d}_2$. The index IC (panels Figure 9a,c) which measures the quantity of energy localized in the cavity with respect to the global energy in the metamaterial with the cavity is higher when considering the larger cavity ($2\tilde{d}_2$). The dependence on the width of the barriers $l$ is very limited when considering a frequency close to the band gap opening, while IC decreases with $l$ when considering a higher frequency, close to the band gap closing. A different pattern is exhibited by the index AIC (Figure 9b,d) which measures the density of energy localized in the cavity with respect to mean energy density in the metamaterial with the cavity. Values of AIC greater than one correspond to systems effectively concentrating energy in the cavity. High concentration is obtained, for both cavity dimensions, considering high frequency and large barriers; the smaller cavity gives a higher concentration (see Figure 9b,d).

<figure><img src="image_11.png" /><figcaption>**Figure 8.** First (a) and second (b) optimal $\tilde{d}_n$ plotted with respect to the frequency $\omega$ inside the first band gap and the width $l$ of the barriers.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** The contours (a,b) show respectively the indexes IC and the AIC for the first optimal cavity widths $2\tilde{d}_1$, while (c,d) for the second optimal cavity widths $2\tilde{d}_2$. The plots are obtained using analytic expressions.</figcaption></figure>

As stated before, a change of the geometry of the unit cell, in terms of both filling fraction and coating thickness, would change the expression of the effective mass density (hence the frequencies of the band gap). Moreover, a change in the filling fraction (or equivalently in the ratio between the radius $R_c$ and the cell size $L$) would also modify the value of the effective stiffness of the barriers. To explore the effect of these variations on the energy localization of the system, we study, in particular, the four cases specified in Table 2.

<table><caption>**Table 2.** Chosen geometry for the unit cell.</caption><thead><tr><td>Case</td><td>R<sub>c</sub>/L</td><td>R<sub>f</sub>/R<sub>c</sub></td></tr></thead><tbody><tr><td>a</td><td>0.475</td><td>0.9</td></tr><tr><td>b</td><td>0.475</td><td>0.4</td></tr><tr><td>c</td><td>0.355</td><td>0.9</td></tr><tr><td>d</td><td>0.355</td><td>0.4</td></tr></tbody></table>

Figure 10 shows the effect of the geometric modifications of the cell on the indexes IC and AIC; the cavity width is fixed to $2\tilde{d}_1$ in all cases. Note that the change of the cell also affects the band gap frequencies. In Figure 10 the represented frequency range always corresponds to the first band gap of the LRM. By comparing Figure 10a,b with Figure 10c,d, one can observe that the concentration of energy is improved by increasing the filling fraction of the LRM. At equal filling fraction, a smaller thickness of the coating (Figure 10a,c) also leads to a higher concentration with respect to the case with a large thickness of the coating (Figure 10b,d, respectively).
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Each panel (a–d) shows the index IC on the first row and AIC on the second row, for the first optimal cavity widths $2\tilde{d}_1$. The following geometrical dimensions are used: (a) $R_c = 0.475L, R_f = 0.9R_c$, (b) $R_c = 0.475L, R_f = 0.4R_c$, (c) $R_c = 0.355L, R_f = 0.9R_c$ and (d) $R_c = 0.355L, R_f = 0.4R_c$.</figcaption></figure>

**5. Conclusions**

This paper investigates the possibility to localize the vibration mechanical energy in a cavity between two barriers constituted by ternary LRMs. The concentration of energy is possible for waves having frequency inside the bandgap of the LRM. The proposed system can be reduced to a one-dimensional problem of wave propagation, thus allowing for a complete analytic solution. We consider the long-wave length regime, and we apply the asymptotic homogenization to obtain the dynamic effective properties of the LRM. The optimal dimension of the cavity for energy localization is given in close form. Results are applied to a specific LRM and validated by comparison with numerical results. The analytic solution allows to highlight the influence of the different material and geometric parameters of the metamaterial on the energy concentration. In particular it is shown that, high filling fraction and small coating thickness of the LRM, besides resulting in a larger band gap at low frequency as shown in [24], improve the energy concentration.

The results obtained in this paper could represent a first step to design a REH with optimal performance in terms of localization of the mechanical energy carried by waves within a specific frequency interval.

The analysis was here restricted to the out-of-plane wave propagation, the case of in-plane propagation, even though leading to more complex expressions, could be treated in a similar way.

The treatment of a real REH, with finite out-of-plane dimension, deserves further studies, as the decoupling between the two cases does not hold anymore. The experimental validation should then be carried out.

**Author Contributions:** Conceptualization, methodology and formal analysis, M.M., C.C. and J.-J.M.; software and validation, M.M.; writing and editing, M.M., C.C. and J.-J.M.; supervision, C.C. and J.-J.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially funded by the Italian MIUR (PRIN project nr.2015LYYXA8 on Multi-scale mechanical models for the design and optimization of micro-structured smart materials and metamaterials).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
## Appendix A. Dynamic Effective Properties of the LRM

In the sub-wavelength regime, the motion problem of the ternary LRM is characterized by a macroscopic scale x and by a microscopic scale y = x/ε, see [21,24,27]. This means that both the displacement field u and the stress field σ can be expanded in the following way:

$$u(\mathbf{x}) = u^0(\mathbf{x}, \mathbf{x}/\epsilon) + \epsilon u^1(\mathbf{x}, \mathbf{x}/\epsilon) + \epsilon^2 u^2(\mathbf{x}, \mathbf{x}/\epsilon) + \dots \quad (A1)$$

$$\sigma(\mathbf{x}) = \sigma^0(\mathbf{x}, \mathbf{x}/\epsilon) + \epsilon\sigma^1(\mathbf{x}, \mathbf{x}/\epsilon) + \epsilon^2\sigma^2(\mathbf{x}, \mathbf{x}/\epsilon) + \dots \quad (A2)$$

By substituting Equations (A1) and (A2) and their derivatives in Equation (1), an effective description of the LRM is achieved. The terms of the same order in ε are then considered. At order −1 one obtains that the displacement field u_m^0(\mathbf{x}, \mathbf{x}/\epsilon) of the matrix is independent from the microscopic variable and will be denoted as U^0(\mathbf{x}). At the leading order, and considering the symmetry of the lattice, the homogenized problem is thus governed by the following equation:

$$\mu_{eff} \operatorname{div} \operatorname{grad} U^0(\mathbf{x}) + \rho_{eff}(\omega) \omega^2 U^0(\mathbf{x}) = 0 \quad (A3)$$

where μ eff is the effective shear modulus. This latter is defined as

$$\mu_{eff} = \frac{1}{|Y|} \int_{Y_m} \mu_m (\operatorname{grad}_y \chi + \mathbf{e}_1) \cdot (\operatorname{grad}_y \chi + \mathbf{e}_1) \, dy \quad (A4)$$

where e_1 is the unit vector of axis x_1 and χ is the solution of the cell problem, associated to the matrix only

$$\left\{
\begin{array}{ll}
\operatorname{div} \operatorname{grad}_y \chi = 0 & \text{in } Y_m \\
\mu_m (\operatorname{grad}_y \chi + \mathbf{e}_1) \cdot \mathbf{n} = 0 & \text{on } \partial Y_c \\
\chi \text{ periodic, } \mu_m \operatorname{grad}_y \chi \cdot \mathbf{n} \text{ anti-periodic on } \partial Y
\end{array}
\right.
\quad (A5)$$

Note that for a general case of a 2D lattice, without the symmetries of the present case, an effective 2 × 2 stiffness matrix intervenes in Equation (A3) whose entries are defined similarly to Equation (A4) by introducing the solution χ² of a second cell problem (with e_2 instead of e_1 in Equation (A5))

In Equation (A3), ρ eff (ω) is the effective mass density, given by:

$$\rho_{eff}(\omega) = \rho_m \frac{|Y_m|}{|Y|} + \rho_f \frac{|Y_f|}{|Y|} \eta(R_f) + \frac{\rho_c}{|Y|} 2\pi \int_{R_f}^{R_c} \eta(r) dr \quad (A6)$$

where the function η(r) reads:

$$\eta(r) = \frac{k_c R_f \rho_f (J_0(k_c r) Y_0(k_c R_f) - Y_0(k_c r) J_0(k_c R_f)) + 2 \rho_c (Y_0(k_c r) J_1(k_c R_f) - J_0(k_c r) Y_1(k_c R_f))}{k_c R_f \rho_f (J_0(k_c R_c) Y_0(k_c R_f) - Y_0(k_c R_c) J_0(k_c R_f)) + 2 \rho_c (Y_0(k_c R_c) J_1(k_c R_f) - J_0(k_c R_c) Y_1(k_c R_f))} \quad (A7)$$

where J_0 and J_1 are Bessel functions of the first kind, Y_0 and Y_1 are Bessel functions of the second kind and k_c = ω√ρ_c/μ_c.

Denoting by ω_n the set of frequencies which make the denominator of Equation (A7) equal to zero, the Equation (A6) tends to −∞ and to +∞ when ω tends to ω_n from above and from below respectively. The intervals where the effective mass is negative represent band gaps, i.e., intervals of frequencies where waves traveling across the metamaterial are attenuated.
References

1. Priya, S.; Inman, D. *Energy Harvesting Technologies*; Springer: Boston, MA, USA, 2009.

2. Ahmed, R.; Mir, F.; Banerjee, S. A review on energy harvesting approaches for renewable energies from ambient vibrations and acoustic waves using piezoelectricity. *Smart Mater. Struct.* **2017**, *26*, 085031.
[CrossRef]

3. Tol, S.; Degertekin, F.; Erturk, A. 3D-printed phononic crystal lens for elastic wave focusing and energy harvesting. *Addit. Manufacturing* **2019**, *29*, 100780. [CrossRef]

4. Deng, T.; Zhang, S.; Gao, Y. Tunability of band gaps and energy harvesting based on the point defect in a magneto-elastic acoustic metamaterial plate. *Appl. Phys. Express* **2020**, *13*, 015503. [CrossRef]

5. Ponti, J.M.D.; Colombi, A.; Ardito, R.; Braghin, F.; Corigliano, A.; Craster, R.V. Graded elastic metasurface for enhanced energy harvesting. *New J. Phys.* **2020**, *22*, 013013. [CrossRef]

6. Carrara, M.; Cacan, M.; Toussaint, J.; Leamy, M.; Ruzzene, M.; Erturk, A. Metamaterial-inspired structures and concepts for elastoacoustic wave energy harvesting. *Smart Mater. Struct.* **2013**, *22*, 065004. [CrossRef]

7. Chen, Z.; Guo, B.; Yang, Y.; Cheng, C. Metamaterials-based enhanced energy harvesting: A review. *Phys. B* *Phys. Condens. Matter* **2014**, *438*, 1–8. [CrossRef]

8. Lv, H.; Tian, X.; Wang, M.Y.; Li, D. Vibration energy harvesting using a phononic crystal with point defect states. *Appl. Phys. Lett.* **2013**, *102*, 2013–2016. [CrossRef]

9. Qi, S.; Oudich, M.; Li, Y.; Assouar, B. Acoustic energy harvesting based on a planar acoustic metamaterial. *Appl. Phys. Express* **2016**, *10*, 263501. [CrossRef]

10. D'Alessandro, L.; Bahr, B.; Daniel, L.; Weinstein, D.; Ardito, R. Shape optimization of solid-air porous phononic crystal slabs with widest full 3D bandgap for in-plane acoustic waves. *J. Comput. Phys.* **2017**, *344*, 465–484. [CrossRef]

11. Yin, X.; Zhang, S.; Xu, G.K.; Zhang, L.Y.; Gao, Z.Y. Bandgap characteristics of a tensegrity metamaterial chain with defects. *Extrem. Mech. Lett.* **2020**, *36*, 100668. [CrossRef]

12. Liu, Z. Locally Resonant Sonic Materials. *Science* **2000**, *289*, 1734–1736. [CrossRef] [PubMed]

13. Krushynska, A.O.; Kouznetsova, V.G.; Geers, M.G. Towards optimal design of locally resonant acoustic metamaterials. *J. Mech. Phys. Solids* **2014**, *71*, 179–196. [CrossRef]

14. Krushynska, A.O.; Miniaci, M.; Kouznetsova, V.G.; Geers, M.G. Multilayered inclusions in locally resonant metamaterials: Two-dimensional versus three-dimensional modeling. *J. Vib. Acoust. Trans. ASME* **2017**, *139*, 3–6. [CrossRef]

15. Moscatelli, M.; Ardito, R.; Driemeier, L.; Comi, C. Band-gap structure in two- and three-dimensional cellular locally resonant materials. *J. Sound Vib.* **2019**, *454*, 73–84. [CrossRef]

16. Comi, C.; Marigo, J.J. Homogenization Approach and Bloch–Floquet Theory for Band-Gap Prediction in 2D Locally Resonant Metamaterials. *J. Elast.* **2020**, *139*, 61–90. [CrossRef]

17. Gonella, S.; To, A.C.; Kam, W. Interplay between phononic bandgaps and piezoelectric microstructures for energy harvesting. *J. Mech. Phys. Solids* **2009**, *57*, 621–633. [CrossRef]

18. Sugino, C.; Erturk, A. Analysis of multifunctional piezoelectric metastructures for low-frequency bandgap formation and energy harvesting. *J. Phys. D Appl. Phys.* **2018**, *51*, 215103. [CrossRef]

19. Moscatelli, M.; Comi, C.; Marigo, J.J. Locally Resonant Materials for Energy Harvesting at Small Scale. In Proceedings of the XXIV AIMETA Conference 2019, Rome, Italy, 15–19 September 2019; Carcaterra, A., Paolone, A., Graziani, G., Eds.; Springer Nature: Cham, Switzerland, 2020; pp. 606–626.

20. Auriault, J.L.; Bonnet, G. Dynamique des composites élastiques périodiques. *Arch. Mech.* **1985**, *37*, 269–284.

21. Sanchez-Hubert, J.; Sanchez-Palencia, E. *Introduction aux Méthodes Asymptotiques et à l'homogénéisation*; Dunod: Masson, Paris, 1992.

22. Auriault, J.L.; Boutin, C. Long wavelength inner-resonance cut-off frequencies in elastic composite materials. *Int. J. Solids Struct.* **2012**, *49*, 3269–3281. [CrossRef]

23. Pham, K.; Maurel, A.; Marigo, J.J. Two scale homogenization of a row of locally resonant inclusions—The case of anti-plane shear waves. *J. Mech. Phys. Solids* **2017**, *106*, 80–94. [CrossRef]

24. Comi, C.; Moscatelli, M.; Marigo, J.J. Two scale homogenization in ternary locally resonant metamaterials. *Mater. Phys. Mech.* **2020**, *44*, 8–18. [CrossRef]
25. Gliozzi, A.S.; Miniaci, M.; Krushynska, A.O.; Morvan, B.; Scalerandi, M.; Pugno, N.M.; Bosia, F. Proof of concept of a frequency-preserving and time-invariant metamaterial-based nonlinear acoustic diode. *Sci. Rep.* **2019**, *9*, 1–9. [CrossRef] [PubMed]

26. Fabry, C.; Pérot, A. Théorie et applications d'une nouvelle méthode de spectroscopie interférentielle. *Ann. Chim. Phys. 1899*, **16**, 115-144.

27. Auriault, J.L.; Boutin, C.; Geindreau, C. (Eds.) *Homogenization of Coupled Phenomena in Heterogenous Media*; ISTE: Arlington, VA, USA, 2009. [CrossRef]

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

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).