# Generalized analytic model for rotational and anisotropic metasolids

Elie Favier, Navid Nemati, Camille Perrot, Qi-Chang He

► To cite this version:

Elie Favier, Navid Nemati, Camille Perrot, Qi-Chang He. Generalized analytic model for rotational and anisotropic metasolids. Journal of Physics Communications, 2018, 2 (3), pp.35035 - 35035. 10.1088/2399-6528/aab5a5. http://hal-01743930

HAL Id: hal-01743930

https://hal.science/hal-01743930

Submitted on 26 Mar 2018

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

Generalized analytic model for rotational and anisotropic metasolids

To cite this article: Elie Favier et al 2018 J. Phys. Commun. 2 035035

View the [article online](#) for updates and enhancements.
<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

4 January 2018

ACCEPTED FOR PUBLICATION

9 March 2018

PUBLISHED

21 March 2018

Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.

Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

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

PAPER

# Generalized analytic model for rotational and anisotropic metasolids

Elie Favier, Navid Nemati¹, Camille Perrot and Qi-Chang He¹

Université Paris-Est, Laboratoire Modélisation et Simulation Multi Echelle, MSME UMR 8208 CNRS, 5 bd Descartes, 77454 Marne-la-Vallée, France

¹ The authors to whom any correspondence should be addressed.

E-mail: navid.nemati@u-pem.fr, nnemati@mit.edu and qi-chang.he@u-pem.fr

**Keywords:** elastic metamaterials, acoustic metamaterials, local resonances, rotational motions, effective density of moment of inertia, homogenization, effective-medium parameters

## Abstract

An analytical approach is presented to model a metasolid accounting for anisotropic effects and rotational mode. The metasolid is made of either cylindrical or spherical hard inclusions embedded in a stiff matrix via soft claddings. It is shown that such a metasolid exhibits negative mass densities near the translational-mode resonances, and negative density of moment of inertia near the rotational resonances. As such, the effective density of moment of inertia is introduced to characterize the homogenized material with respect to its rotational mode. The results obtained by this analytical continuum approach are compared with those from discrete mass-spring model, and the validity of the latter is discussed. Based on derived analytical expressions, we study how different resonance frequencies associated with different modes vary and are placed with respect to each other, in function of the mechanical properties of the coating layer. We demonstrate that the resonances associated with additional modes that are taken into account, i.e., axial translation for cylinders, and rotations for both cylindrical and spherical systems, can occur at lower frequencies compared with the previously studied plane-translational modes.

## 1. Introduction

Elastic or acoustic metamaterials are structures with subwavelength units that exhibit unusual macroscopic parameters within certain frequencies, such as negative mass density [1–3], negative elastic bulk modulus [4, 5], simultaneous negative mass density and elastic bulk modulus [6–11], or negative density and shear modulus [12, 13], and negative index of refraction [14]. Metamaterials with unconventional constitutive parameters have broadly extended the ability of manipulation and control of mechanical wave propagation. Wave control can arise from the formation of band gaps, that are the frequency bands where the propagation becomes forbidden. Differently from the Bragg band-gaps in phononic crystals [15, 16] that are produced based on the collective effects of the periodically-arranged scatterers in the medium, the spectral gaps manifested in metamaterials originates in localized resonances of the material building-blocks. Thus, in metamaterials the production of stop bands depends on the internal structure of the building units. Another mechanism emerged to achieve extreme material-parameters is based on space coiling [17, 18].

Over the last decade, research on anisotropic metamaterials has enhanced the ability to control sound and elastic wave propagation, leading to the proposal and fabrication of new devices [19–26]. The first acoustic metamaterial that was realized is a block of stiff material (epoxy) with uniformly (or periodically) dispersed locally resonant structural units that consist of a high-density lead sphere coated with a soft material. This material exhibits resonance-based spectral-gaps in low-frequency regime where the sonic wavelength is of two orders of magnitude larger than the lattice constants. Thus, this designed material was a major progress as sound attenuation is challenging in low frequencies. An analytical model has been proposed to capture the essence of physics dealing with this type of materials [2]. Since in practice the lead sphere and the host are much stiffer compared with the coating material, in this model the host and coated sphere have been approximated to be rigid within the long-wavelength limit. Considering the two-dimensional (2D) system with coated cylindrical
inclusions, and three-dimensional (3D) with coated spheres, effective dynamics of the metasolid was derived in resonance-induced band-gap regimes related to the translational motions. Analytical expressions were given for the effective mass density which becomes negative near the resonance.

In this paper, generalizing the simple model proposed in [2], we describe the dynamics of the 3D three-component composite taking into account the rotational motions of the lead cylinder or sphere, and the matrix. It is shown that when the inclusions are cylindrical, because of their geometrical anisotropy, wave propagation in the direction perpendicular to the axis of cylinders gives rise to different effective density in comparison with the effective density describing the propagation perpendicular to that axis. Components of the anisotropic effective density are given through analytical expressions. Furthermore, we introduce and analytically calculate the effective density of the moment of inertia, as a parameter that describes the effective rotational mode of the material. It has been previously demonstrated that, when the cladding is very soft, the stop-band frequency can be very low [1]. This makes such a composite with a practical size interesting for various applications including sound and vibration isolations. Here, allowing for additional degrees of freedom (DOFs) related to rotational motions and translation along cylinder axis, we show that the model predicts spectral gaps characteristic for each mode of propagation. In particular, with the material components chosen in [2], we demonstrate that in the case of spherical inclusions, the band gaps related to rotational mode occur in lower frequencies compared with those corresponding to translational modes, where the respective effective parameters become negative. Also, for the case of cylindrical inclusions, comparison between the stop bands for translational modes along and perpendicular to the cylinders shows that those that occur along the cylinder axis are formed at lower frequencies. In general, for materials made of spherical or cylindrical inclusions, we analyze the occurrence of local resonances related to each mode on the frequency axis, as function of coating-material properties. It turns out that accounting for additional DOFs that results in additional propagation modes improve notably the tunability of the material in terms of size and relevant frequencies. This is significant for applications such as sound and vibration proofing, which requires small-size materials enabling the low-frequency wave attenuation.

Rotational modes in periodic solid composites have been previously studied [10, 19, 27–31]. Rotary resonances have been estimated by a simple model in 2D square lattice of glass cylinders in epoxy [27]. Numerical analysis was performed to obtain double negative properties, for density and bulk modulus, by utilizing simultaneously the local translational and rotational resonances of a chiral microstructure with a unit cell of three-component solid media composed of a chirally soft-coated heavy cylinder core embedded in a elastic matrix [10]. Periodic system of rotating resonators coated by an anisotropic heterogeneous elastic material embedded in a matrix, which was modeled by an asymptotic approach, has been reported to produce negative refraction [29]. Thus, by modelling the anisotropic material that is placed between the core and matrix as ligaments, and by coupling the translational motions with rotations, it has been numerically demonstrated that this structured medium could be employed for lensing and mode localization. Due to simultaneous translational and rotational local-resonances, a metasolid made of a chiral microstructure with a single phase material served as a system to experimentally achieve negative refraction [30]. Accounting for rotational modes, a mass-spring model was used to reproduce band gaps in a 2D phononic crystal made of rigid cylinders in epoxy [28]. This model was enhanced to describe dispersion relation for rotational modes of a 2D square array of rubber-coated steel cylinders in epoxy [31]. However, no analytical approach for rotational metasolids based on full elastodynamics has been reported so far.

Here, we study rotational modes and analyze their associated homogenized properties by modeling analytically the classical three-component subwavelength structures through taking general assumptions and using direct first-principle elastodynamics without any fitting parameters. This offers us not only essential physical insights for the effective dynamics of the metasolid but also a quick and efficient guide to design materials with local rotational-resonances. Additionally and in parallel, for each mode we have systematically calculated the effective parameters of the metasolid based on the discrete mass-spring modeling, and compared the results with those arising from the continuum elastodynamic model. This comparison clarifies the limits of the mass-spring based model, in particular in terms of describing the effective-material dynamics related to the second local-resonance phenomenon that manifests itself for all translational and rotational modes.

In the following, we introduce the elastodynamics equations at microscale in section 2 for an arbitrary shape of the inclusions. We then study the case of cylindrical inclusions in section 3, followed by the analysis for spherical inclusions in section 4. In section 5, the effective parameters for translational and rotational modes are calculated for the homogenized media including microscopically either cylindrical or spherical resonators. In the concluding section 6, the main results of the paper are briefly summarized.
<figure><img src="image_4.png" /><figcaption>Figure 1. Illustration of the cross-section of a 3D unit cell with coated inclusion of arbitrary shape: rigid inclusion $\Omega_a$ surrounded by the elastic cladding $\Omega$ and embedded in the rigid matrix $\Omega_b$. The boundary regions of the cladding are denoted by $\Gamma_a$ and $\Gamma_b$.</figcaption></figure>

## 2. Microscale equations

The unit cell under consideration consists of a hard inclusion $\Omega_a$ embedded in a hard matrix $\Omega_b$ via a soft cladding $\Omega$ (figure 1). The interfaces $\Gamma_a$ and $\Gamma_b$ located between $\Omega_a$ and $\Omega$, and between $\Omega_b$ and $\Omega$, respectively, are assumed to be perfect. The materials constituting $\Omega_b$ and $\Omega_a$ are very stiff with respect to the material forming $\Omega$. For this reason, and within the long-wavelength limit, we make the assumption that the inclusion $\Omega_a$ and matrix $\Omega_b$ are rigid while the cladding $\Omega$ is linearly elastic [2, 32].

The elastodynamic formulation of the unit-cell system is based on the assumptions that only the cladding is deformable and that the displacements of the inclusion and matrix are small and the strains inside the soft layer are infinitesimal. Consequently, the kinematics of the inclusion and the matrix can be formulated by

$$ \hat{u}_{\alpha}(r, t) = \hat{V}_{\alpha}(t) + \hat{\Theta}_{\alpha}(t) \times r, \quad \text{on } \Gamma_{\alpha}, \qquad (1) $$

where $\alpha = a, b, r$ is the position vector relative to the origin $O$ of a 3D space, and $\hat{u}$ represents the displacement field generated by a translational displacement vector $\hat{V}_\alpha$ and an infinitesimal rotation specified by the 3D vector $\hat{\Theta}_\alpha$. Considering the time-harmonic motion of the foregoing unit system, we write $\hat{u}(r, t) = u(r, \omega)e^{i\omega t}$, $\hat{V}_\alpha(t) = V_\alpha(\omega)e^{i\omega t}$, and $\hat{\Theta}_\alpha(t) = \Theta_\alpha(\omega)e^{i\omega t}$, where $i = \sqrt{-1}$ and $\omega$ is the angular frequency. The motion of the cladding made of a linearly elastic isotropic homogeneous material is described by the following Navier equation together with the displacement continuity across the interfaces $\Gamma_b$ and $\Gamma_a$:

$$ (\lambda + 2\mu)\nabla(\nabla \cdot u) - \mu \nabla \times (\nabla \times u) = -\rho\omega^2 u, \quad \text{in } \Omega \qquad (2a) $$

$$ u(r, \omega) = V_{\alpha}(\omega) + \Theta_{\alpha}(\omega) \times r, \quad \text{on } \Gamma_{\alpha}, \quad \alpha = a, b \qquad (2b) $$

where $\lambda$ and $\mu$ are the Lamé constants and $\rho$ is the mass density of the coating material.

For later use, it is convenient to introduce the Helmholtz decomposition:

$$ u = \nabla\Phi + \nabla \times \Psi, \quad \nabla \cdot \Psi = 0, \qquad (3) $$

where $\Phi$ is a scalar field and $\Psi$ a vector field. Using equation (3), the equation (2a) can be recast into the following ones

$$ \nabla^2 \Phi + h^2 \Phi = 0, \qquad (4a) $$

$$ \nabla^2 \Psi + \kappa^2 \Psi = 0 \qquad (4b) $$

where $h = \omega/c_l$ and $\kappa = \omega/c_b$, with $c_l = \sqrt{(\lambda + 2\mu)/\rho}$ and $c_t = \sqrt{\mu/\rho}$ being the celerities of the longitudinal and transverse elastic waves, respectively, in an infinite linearly elastic isotropic medium. In the following, we shall analytically solve the Helmholtz equations (4) with the boundary conditions (2b) for two special cases of major interest.

## 3. Media with cylindrical inclusions

The first special case concerns a material with structural units consisting of three coaxial cylinders of height $L$, uniformly or periodically distributed. The inner cylinder $\Omega_a$ is a long cylinder of revolution one of radius $a \ll L$, while the outer cylinder $\Omega_b$ is bounded by an inner cylinder of circular cross-section of radius $b > a$ and an outer cylindrical surface with square cross-section (figure 2(a)).
<figure><img src="image_5.png" /><figcaption>Figure 2. Illustration of a unit cell with cylindrical or spherical inclusion in xy plane (a), and coordinate system in 3D space (b).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 3. Illustration of boundary conditions at the interfaces of the cladding, $r = a$ and $r = b$, for plane translation (a), rotation (b), and axial translation (c) motions.</figcaption></figure>

Regarding the geometry of the unit cell, it is convenient to use the cylindrical coordinates with associated orthonormal vectors $e_r$, $e_\theta$, and $e_z$ (figure 2(b)). As it was assumed in section 2, the inner and outer cylinders are rigid, whereas the cylindrical cladding is deformable. In addition, since the dimensions of the unit cell are such that $L \gg b > a$, the six DOFs of the motion of $\Omega_a$ relative to $\Omega_b$ can be reduced to four: two translations in the transverse plane, one translation along the cylinder axis and one rotation about the latter. Thus, with no loss of generality, the boundary conditions (2b) for the cladding $\Omega$ (figure 3) can be written as ($\alpha = a, b$)

$$ \mathbf{u}|_{r=\alpha} = U_{\alpha} \{ \cos \theta_{\alpha} \mathbf{e}_x + \sin \theta_{\alpha} \mathbf{e}_y \} + \alpha \Theta_{\alpha} \mathbf{e}_z \times \mathbf{e}_r + T_{\alpha} \mathbf{e}_z, \quad (5) $$

where $U_\alpha$ represents the amplitude of the displacement in the plane direction $e_\alpha$ resulting from the rotation of $e_x$ by the angle $\theta_\alpha$ around $e_z$. The last term in the above expression stands for the axial translation, with $T_\alpha$ the amplitude of the displacement along $e_z$. The last DOF, i.e. the rotation around the common axis $e_z$, is proportional to $\Theta_\alpha e_z \times e_r$, where $\Theta_\alpha$ is the infinitesimal amplitude of the rotation. To link the general form of the boundary conditions (2b) to equation (5) associated with the case of cylindrical inclusions, it is obvious that $V_\alpha \equiv U_\alpha e_\alpha + T_\alpha e_z$ and $\Theta_\alpha \equiv \Theta_\alpha e_z$ (figure 3).

Considering the cylindrical nature of both the geometry and boundary conditions, the displacement field $\mathbf{u}$ is clearly independent of the axial coordinate $z$, thereby the potentials $\Phi$ and $\Psi$ are also independent of $z$. Moreover, following Love [33], $\Psi$ can be chosen to take the form $\Psi = \psi_\theta(r, \theta)e_\theta + \Psi_z(r, \theta)e_z$. The condition $\nabla \cdot \Psi = 0$ in equation (3) is now satisfied if and only if $\Psi_\theta$ is independent of $\theta$. Thus, the potential $\Psi$ can be expressed as $\Psi = \Psi_\theta(r)e_\theta + \Psi_z(r, \theta)e_z$.

The general form of the solutions to the Helmholtz equations (4) can be written as

$$ \Phi(r, \theta) = \sum_{n=0}^{+\infty} \{A_n J_n(hr) + B_n Y_n(hr)\} \times \{C_n \cos(n\theta) + D_n \sin(n\theta)\}, \quad (6a) $$

$$ \Psi_z(r, \theta) = \sum_{n=0}^{+\infty} \{E_n J_n(\kappa r) + F_n Y_n(\kappa r)\} \times \{G_n \cos(n\theta) + H_n \sin(n\theta)\}, \quad (6b) $$
$$ \Psi_{\theta}(r) = A J_{1}(\kappa r) + B Y_{1}(\kappa r), \qquad (6c) $$

where $J_n$ (resp. $Y_n$) stands for the n-th-order Bessel function of the first (resp. second) kind. In these solutions, $A_n$, $B_n$, $C_n$, $D_n$, $E_n$, $F_n$, $G_n$, $H_n$, $A$, and $B$ are unknown constants to be determined.

Boundary conditions (5) and Helmholtz decomposition (3) expressed in cylindrical coordinates lead to

$$ \left[ \frac{\partial \Phi}{\partial r} + \frac{1}{r} \frac{\partial \Psi_z}{\partial \theta} \right]_{r=\alpha, \theta} = U_{\alpha} \cos(\theta_{\alpha} - \theta), \qquad (7a) $$

$$ \left[ \frac{\partial \Phi}{\partial \theta} - \frac{\partial \Psi_z}{\partial r} \right]_{r=\alpha, \theta} = U_{\alpha} \sin(\theta_{\alpha} - \theta) + \alpha \Theta_{\alpha}, \qquad (7b) $$

$$ \left[ \frac{1}{r} \frac{d(r\Psi_{\theta})}{dr} \right]_{r=\alpha, \theta} = T_{\alpha}. \qquad (7c) $$

Now, we are able to determine all unknown constants involved in equations (6).

Owing to the linearity of the problem, general solution for the displacement field inside the cladding $\Omega$ can be written as $\mathbf{u} = \mathbf{u}^{(p)} + \mathbf{u}^{(a)} + \mathbf{u}^{(r)}$, where $\mathbf{u}^{(p)}$ refers to the translation in xy plane, $\mathbf{u}^{(a)}$ to the axial displacement in $e_z$ direction, and $\mathbf{u}^{(r)}$ to the rotation around the axis (O, $e_z$). The first displacement field $\mathbf{u}^{(p)}$ can be expressed as

$$ \begin{aligned} \mathbf{u}^{(p)} = & A_1 C_1 \left\{ \frac{J_1(hr)}{r} e_x - h J_2(hr) \cos \theta e_r \right\} + B_1 C_1 \left\{ \frac{Y_1(hr)}{r} e_x - h Y_2(hr) \cos \theta e_r \right\} \\ & + E_1 H_1 \left\{ \frac{J_1(\kappa r)}{r} e_x + \kappa J_2(\kappa r) \sin \theta e_\theta \right\} + F_1 H_1 \left\{ \frac{Y_1(\kappa r)}{r} e_x + \kappa Y_2(\kappa r) \sin \theta e_\theta \right\} \\ & + A_1 D_1 \left\{ \frac{J_1(hr)}{r} e_y - h J_2(hr) \sin \theta e_r \right\} + B_1 D_1 \left\{ \frac{Y_1(hr)}{r} e_y - h Y_2(hr) \sin \theta e_r \right\} \\ & + E_1 G_1 \left\{ -\frac{J_1(\kappa r)}{r} e_y + \kappa J_2(\kappa r) \cos \theta e_\theta \right\} + F_1 G_1 \left\{ -\frac{Y_1(\kappa r)}{r} e_y + \kappa Y_2(\kappa r) \cos \theta e_\theta \right\} \end{aligned} \qquad (8) $$

where the first four terms, and the last four terms describe the displacements generated by the translations along x-axis and y-axis, respectively. The unknown coefficients in the above expression are the solutions to the linear systems (A1) and (A2) given in appendix A. As such, we find in different terms, the solution of the plane translation given in [2] corresponding to the two DOFs related to translational motions.

The displacement field due to the translation along the cylinders' axis, i.e. $e_z$, is expressed as

$$ \mathbf{u}^{(a)} = \kappa \{ A J_0(\kappa r) + B Y_0(\kappa r) \} e_z. \qquad (9) $$

Finally the displacement field related to the rotational motions is provided by

$$ \mathbf{u}^{(r)} = \kappa \{ E_0 G_0 J_1(\kappa r) + F_0 G_0 Y_1(\kappa r) \} e_{\theta}. \qquad (10) $$

After applying the boundary conditions (7), we obtain explicitly<fnref n="2" />

$$ \begin{aligned} \mathbf{u}^{(a)} &= \left( J_0(\kappa r) \quad Y_0(\kappa r) \right) \begin{pmatrix} J_0(\kappa a) & Y_0(\kappa a) \\ J_0(\kappa b) & Y_0(\kappa b) \end{pmatrix}^{-1} \begin{pmatrix} T_a \\ T_b \end{pmatrix} e_z \\ &= \begin{pmatrix} J_0(\kappa r) Y_0(\kappa b) - J_0(\kappa b) Y_0(\kappa r) \\ J_0(\kappa a) Y_0(\kappa b) - J_0(\kappa b) Y_0(\kappa a) \end{pmatrix} \begin{pmatrix} J_0(\kappa a) Y_0(\kappa r) - J_0(\kappa r) Y_0(\kappa a) \\ J_0(\kappa b) Y_0(\kappa b) - J_0(\kappa b) Y_0(\kappa a) \end{pmatrix} \begin{pmatrix} T_a \\ T_b \end{pmatrix} e_z, \end{aligned} \qquad (11) $$

$$ \begin{aligned} \mathbf{u}^{(r)} &= \left( J_1(\kappa r) \quad Y_1(\kappa r) \right) \begin{pmatrix} J_1(\kappa a) & Y_1(\kappa a) \\ J_1(\kappa b) & Y_1(\kappa b) \end{pmatrix}^{-1} \begin{pmatrix} a\Theta_a \\ b\Theta_b \end{pmatrix} e_\theta \\ &= \begin{pmatrix} a \frac{J_1(\kappa r) Y_1(\kappa b) - J_1(\kappa b) Y_1(\kappa r)}{J_1(\kappa a) Y_1(\kappa b) - J_1(\kappa b) Y_1(\kappa a)} & b \frac{J_1(\kappa a) Y_1(\kappa r) - J_1(\kappa r) Y_1(\kappa a)}{J_1(\kappa a) Y_1(\kappa b) - J_1(\kappa b) Y_1(\kappa a)} \\ \frac{J_1(\kappa r) Y_1(\kappa b) - J_1(\kappa b) Y_1(\kappa r)}{J_1(\kappa a) Y_1(\kappa b) - J_1(\kappa b) Y_1(\kappa a)} & \frac{\Theta_a}{\Theta_b} \end{pmatrix} e_\theta \end{aligned} \qquad (12) $$

It is important to note that, these axial and rotational contributions of the motion involves only $\kappa$ (and not $h$). This is due to the state of shear dictating such kinematics. This constitutes a key point in order to achieve an anisotropic effective density, and obtain different characteristics for different directions.

We now use the equations of motion and the stresses that are known in the cladding, to deduce the relation between the displacements of the core cylinder and those of the embedding matrix. The equations of motion are written as

$$ -M_a \omega^2 V_a = \int_{\Gamma_a} \sigma \cdot n \, dS \qquad (13a) $$

<footnote n="2">² In the limit $a \to 0$, i.e., making disappear the rigid core, the results presented in [34] can be easily obtained as a particular case of these general expressions.</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 4. Displacement of the cylindrical core normalized to and exited by that of the matrix, based on the continuum model and discrete mass-spring model, and related to three modes: plane translation, axial translation, and rotation. The displacement field patterns inside the elastic coating layer (between the core and the matrix) for the frequency points (a), (b), ..., (h) will be shown in figure 5.</figcaption></figure>

$$-I_a \omega^2 \Theta_a = \int_{\Gamma_a} r \times (\sigma \cdot n) dS \qquad (13b)$$

where $M_a = \rho_a(\pi a^2 L)$ is the mass of the core cylinder with the density $\rho_a$, and $I_a = M_a a^2/2$ is its moment of inertia. The stress tensor $\sigma = \mu(\nabla u + \nabla u^T) + \lambda(\nabla \cdot u)\mathcal{I}$, with $u^T$ the transpose of $u$ and $\mathcal{I}$ the identity matrix, is known since the displacement $u$ in the cladding is fully determined. After some simple calculations, we obtain the relationship between the displacements of the core and those of the matrix, as follows

$$\left\{
\begin{aligned}
U_a &= \frac{-\gamma g^{(p)}(\omega)}{R^{(p)}(\omega) - \rho_a / \rho} & U_b &:= H^{(p)}(\omega) U_b \\
T_a &= \frac{-g^{(a)}(\omega)}{R^{(a)}(\omega) - \rho_a / \rho} & T_b &:= H^{(a)}(\omega) T_b \\
\Theta_a &= \frac{-\gamma g^{(r)}(\omega)}{R^{(r)}(\omega) - \rho_a / \rho} & \Theta_b &:= H^{(r)}(\omega) \Theta_b
\end{aligned}
\right.
\qquad (14)$$

where $\gamma := b/a$. The expressions of $R^{(p)}(\omega)$ and $g^{(p)}(\omega)$, involved in plane translational motion are given in the appendix A (equation (A4)), which can be found also in [2]. The other components related to axial translation and rotation are expressed as

$$\left\{
\begin{aligned}
g^{(a)}(\omega) &= \frac{2}{\kappa a} \frac{J_0(\kappa a) Y_1(\kappa a) - J_1(\kappa a) Y_0(\kappa a)}{J_0(\kappa a) Y_0(\kappa b) - J_0(\kappa b) Y_0(\kappa a)} \\
R^{(a)}(\omega) &= \frac{2}{\kappa a} \frac{J_1(\kappa a) Y_0(\kappa b) - J_0(\kappa b) Y_1(\kappa a)}{J_0(\kappa a) Y_0(\kappa b) - J_0(\kappa b) Y_0(\kappa a)}
\end{aligned}
\right.
\qquad (15a)$$

$$\left\{
\begin{aligned}
g^{(r)}(\omega) &= \frac{4}{\kappa a} \frac{J_1(\kappa a) Y_2(\kappa a) - J_2(\kappa a) Y_1(\kappa a)}{J_1(\kappa a) Y_1(\kappa b) - J_1(\kappa b) Y_1(\kappa a)} \\
R^{(r)}(\omega) &= \frac{4}{\kappa a} \frac{J_2(\kappa a) Y_1(\kappa b) - J_1(\kappa b) Y_2(\kappa a)}{J_1(\kappa a) Y_1(\kappa b) - J_1(\kappa b) Y_1(\kappa a)}
\end{aligned}
\right.
\qquad (15b)$$

Figure 4 shows the evolution of the displacement of the core cylinder for the three modes in function of frequency, represented by $H^{(p)}(\omega)$, $H^{(a)}(\omega)$ and $H^{(r)}(\omega)$, using the same material properties as in [1] (shown in table 1) and geometrical parameters with the filling fraction for coated cylinders being 40%, 5.0 mm for the radius of the lead cylinder, and the coating thickness is 2.5 mm. For these modes of plane translation, axial translation, and rotation, the displacement field patterns inside the cladding close to the first and second resonance frequencies are depicted in figure 5. Through figure 5, we observe that the field amplitudes are maximum on the boundary of the rigid cylinder at the first resonance frequency, while it is maximum inside the cladding at the second resonance frequency. The field patterns in figures 5(a) to 5(h) correspond to the frequency points (a) to (h) in figure 4. Behaviors of $H^{(a)}(\omega)$ and $H^{(r)}(\omega)$ are similar to the behavior of $H^{(p)}(\omega)$, although obviously the resonances relating to different types of motions, when $R(\omega) \equiv \rho_a/\rho$, are located at different frequencies. Indeed, the first resonances for axial, rotational and plane modes are obtained, respectively, at $f_0^{(a)} \approx 128$ Hz, $f_0^{(r)} \approx 217$ Hz and $f_0^{(p)} \approx 355$ Hz. These resonances concern rigid-body mode insofar as the
<table><caption>**Table 1.** Material properties of the three-component unit-cell.</caption><thead><tr><th>Material</th><th>Epoxy (matrix)</th><th>Silicone (cladding)</th><th>Lead (core)</th></tr></thead><tbody><tr><td>λ[Pa]</td><td>4.43 × 10<sup>9</sup></td><td>6 × 10<sup>5</sup></td><td>4.23 × 10<sup>10</sup></td></tr><tr><td>μ [Pa]</td><td>1.59 × 10<sup>9</sup></td><td>4 × 10<sup>4</sup></td><td>1.49 × 10<sup>10</sup></td></tr><tr><td>ρ [kg.m<sup>-3</sup>]</td><td>1.18 × 10<sup>3</sup></td><td>1.3 × 10<sup>3</sup></td><td>11.6 × 10<sup>3</sup></td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Displacement field distributions inside the coating layer close to the first and second resonance frequencies for the modes related to plane translation, axial translation, and rotation. The displacements are normalized such that they have fixed unit value on the outer boundary of the coating, i.e., at $r = b$ in figure 3. The field patterns in (a), (b), ..., (h) are related to respective frequency points shown in figure 4.</figcaption></figure>

core vibrates as a whole within the elastic cladding. In this case, the dynamics of the structural unit suggests that the unit can be reasonably assimilated to a simple mass-spring system for the translational [35] as well as rotational mode [32], by taking the hard cylinder as the mass, and the cladding as the spring; the latter being attached to the mass on one side and to the rigid matrix on its other side. The mass-spring system is exited by the motion of the matrix.

Based on this mass-spring model the first resonance frequencies $f_0^{(p)}$, $f_0^{(a)}$, and $f_0^{(r)}$ for the motions related to plane translation, axial translation, and rotation, respectively, can be obtained as follows

$$f_0^{(p)} = f_0 \left\{ \frac{8(1 - \nu)(3 - 4\nu)\gamma^2}{(3 - 4\nu)^2 \ln \gamma - \frac{\gamma^2 - 1}{\gamma^2 + 1}} \right\}^{1/2}, \quad f_0^{(a)} = f_0 \left( \frac{2\gamma^2}{\ln \gamma} \right)^{1/2}, \quad f_0^{(r)} = f_0 \left( \frac{8\gamma^4}{\gamma^2 - 1} \right)^{1/2}, \quad (16)$$

where $f_0 = (1/2\pi)\sqrt{\mu/\rho_a b^2}$. With our present configuration, the above formulas give $f_0^{(a)} \approx 131$ Hz, $f_0^{(r)} \approx 224$ Hz, and $f_0^{(p)} \approx 360$ Hz. The good agreement of the latter frequencies with the resonance frequencies from the continuum model allows a very fast and simple estimation of the first resonance frequencies, especially in terms of geometrical and material parameters associated with each constituent. To obtain the displacement based on the mass-spring model, we need to calculate the effective spring constants for each of the translational and rotational modes. The corresponding expressions for the effective spring constants in terms of microscale material parameters, relating to the three modes, are found in appendix B.

The study of relations (16) makes it possible to demonstrate that $f_0^{(a)}$ always remains smaller than the other two resonance frequencies, regardless of the geometrical parameters and material properties of the coating layer. This can be observed in figure 6. The relations (16) also show that unlike the Poisson's ratio $\nu = \lambda/\{2(\lambda + \mu)\}$, the Lam coefficient $\mu$ does not change the positions of the frequencies $f_0^{(p)}$, $f_0^{(a)}$, and $f_0^{(r)}$ relative to each other. More precisely, for Poisson's ratio laying in the interval $0 \le \nu \le 1/4$, the rotational resonance always manifests itself at higher frequency compared with the plane translation resonance ($f_0^{(p)} < f_0^{(r)}$), whereas for Poisson's ratio $1/4 < \nu \le 1/2$ the relative position of these two frequencies evolves with the form factor $\gamma = b/a$. This implies that for small inclusions such that $\gamma_c(\nu) < \gamma$, we have $f_0^{(p)} < f_0^{(r)}$, and vice versa otherwise; where $\gamma_c(\nu)$ is the unique solution to the equation $f_0^{(p)}(\gamma) = f_0^{(r)}(\gamma)$, i.e.,
$$ \ln \gamma = \frac{1}{(3 - 4\nu)^2} \left\{ (1 - \nu)(3 - 4\nu) \left(1 - \frac{1}{\gamma^2}\right) + 1 - \frac{2}{\gamma^2 + 1} \right\}. $$ We note that $\gamma_c(\nu)$ is an increasing function, varying from 1 for $\nu = 1/4$, to about 3.80 when $\nu = 1/2$. Figure 6 illustrates these remarks; in particular it shows the evolution of $f_0^{(p)}$ with $\gamma$ for three representative values of $\nu$: $1/8$ ($0 \le \nu \le 1/4$), $3/8$ ($1/4 < \nu \le 1/2$),
<figure><img src="image_9.png" /><figcaption>Figure 6. Evolution of the first resonance frequencies normalized by $f_0$ and estimated by the discrete mass-spring model (equation (16)), with the form factor $\gamma = b/a$. The dependence of plane-translation mode on the Poisson's ratio $\nu$ of the coating is also illustrated. The Poisson's ratio takes two representative values 1/8 and 3/8 in the intervals [0, 1/4] and [1/4, 1/2], respectively. The value of $\nu = 0.47$ corresponds to the material parameters studied as an example throughout this paper (see table 1). The axial-translation and rotation modes are independent of $\nu$.</figcaption></figure>

and 0.47 that is the value for the structure studied in particular in this paper (table 1), laying also in the interval $1/4 < \nu \le 1/2$.

Although, as seen in figure 4, the mass-spring model describes well the dynamics of the material at lower frequencies including the first resonances, it fails to predict the second resonance behavior for all of the three modes. This is due to the fact that the displacement field at the second resonance is mainly produced by the standing waves inside the coating layer, while the lead cylinder follows the motion of the matrix and remains practically at rest (see [2] for plane translation mode). We know that the spring, modeling the coating layer, has been assumed to be massless, and consequently this discrete model by its construction cannot describe the second resonance. That being said, if we consider a more complex discrete model by taking into account a mass for the mode-related springs a second local resonance in the elastic medium will appear. However, the outcome of latter model does not correspond to that of continuum model; particularly because in the mass-spring model the springs function in one dimension, while at the second local resonance, the displacement fields in the coating layer can be distributed in three dimensions taking, e.g. the shear effects into account. The distribution of the displacement field in the xy plane for each of the three modes is shown in figure 5, close to the first and to the second resonance frequencies. In this figure, comparison of the field pattern at the first resonance for either of the modes with those arising from the same mode but at the second resonance, shows that at the second resonance the displacement field is distributed in space to a significantly greater extent than that related to the first resonance. In addition, in the mass-spring model, the evaluation of the spring constants is performed in a quasi-static regime, which is different from the frequency-dependent nature of the second resonance occurring inside the elastic layer at higher frequencies. The resonance frequencies in figures 5(a) to 5(h), correspond to the points (a), (b), ..., and (h) in figure 4. The second resonance of the structure almost coincides with the first cladding resonance. This first cladding frequency can be estimated by searching for the first root of the determinant, which is a classical technique for modal analysis. In fact, we can proceed to find the solutions, by imposing a zero-displacement for the lead cylinder within the continuum model, assuming that the resonance occurs only inside the elastic medium. Thus, we are led to determine the first root of $\Delta^{(p)}$, $\Delta^{(a)}$, and $\Delta^{(r)}$ in order to find the second resonance frequency of the plane translation, axial translation, and rotation mode, respectively, through the following equations

$$ 
\begin{aligned}
\Delta^{(p)}(\omega) &:= \det \mathcal{D} = \det \begin{pmatrix} [E(ha)] & [F(\kappa a)] \\ [E(hb)] & [F(\kappa b)] \end{pmatrix} = 0, \\
\Delta^{(a)}(\omega) &:= J_0(\kappa a) Y_0(\kappa b) - J_0(\kappa b) Y_0(\kappa a) = 0, \\
\Delta^{(r)}(\omega) &:= J_1(\kappa a) Y_1(\kappa b) - J_1(\kappa b) Y_1(\kappa a) = 0.
\end{aligned} 
\quad (17) $$

We refer to the equation (A1) in appendix A for the derivation of the first relation and also the definition of the associated block matrix $\mathcal{D}$. The two last relations arise from the equations (11) and (12). The first and second local-resonance frequencies related to different modes, and calculated by different methods, are collected in table 2, in order to easily compare these values associated with their respective method of calculation. In particular, it can be noticed that the relations (17) provide good estimations of the second resonance frequencies.
<table><caption>**Table 2.** Fist resonance (R1) and second resonance (R2) frequencies obtained by the continuum model, discrete model (equation (16)), and estimation (equation (17)), for plane translation, axial translation, and rotation modes.</caption><thead><tr><th>Frequency</th><th>Plane (Hz)</th><th>Axial (Hz)</th><th>Rotation (Hz)</th></tr></thead><tbody><tr><td>Continuum model: R1</td><td>355</td><td>128</td><td>217</td></tr><tr><td>Discrete model: R1</td><td>360</td><td>131</td><td>224</td></tr><tr><td>Continuum model: R2</td><td>1236</td><td>1120</td><td>1140</td></tr><tr><td>Estimation: R2</td><td>1236</td><td>1107</td><td>1116</td></tr></tbody></table>

### 4. Media with spherical inclusions

Here, we study the same kind of problem as in the previous section with the only difference that the geometry of the structural unit is changed: the cylinder of radius $a$ and height $L \gg a$ is replaced by a hard ball of radius $a$. An embedding matrix has now spherical cavities of radius $b$, that are filled with the hard spheres and a concentric elastic cladding such that $a < R < b$ (figure 2(a)). All of these structural components are made of the same materials as in the case of cylindrical system (see table 1), so that the same simplifications can be applied to the corresponding kinematics. By virtue of the analysis in the previous section with coaxial cylinders, it is evident to claim, without any calculations, that for an arbitrary motion of the embedding matrix as a rigid body, the same type of motion is induced for the embedded hard ball. The superposition principle, verified for the previous case, suggests us that translations and rotations are totally uncoupled. A translation in a fixed direction or a rotation around a fixed direction, will lead, respectively, to the same translation or rotation for the hard ball, but with possibly different amplitude and/or phase.

Moreover, assuming a total isotropy of the problem, it is now convenient to consider only one translation in a certain direction and, independently, one rotation around a fixed axis to describe the whole kinematics of this problem. The former problem regarding the translational mode based on the continuum model has been solved in [2]. In fact, regarding the translational mode, we only need to replace the Bessel functions in the cylindrical system by spherical Bessel functions and use the formalism explained in appendix A. The aim of this section is to describe the kinematics of the rotational mode according to the continuum model, and compare with its counterpart built upon the discrete mass-spring model. Related equations of motion to be solved, for the cladding zone and its boundaries are stated as follows

$$ (\lambda + 2\mu)\nabla(\nabla \cdot \mathbf{u}) - \mu \nabla \times (\nabla \times \mathbf{u}) = -\rho\omega^2\mathbf{u}, \quad (18a) $$

$$ \mathbf{u}|_{R=\alpha} = \Theta_{\alpha} \times \mathbf{r}|_{R=\alpha}, \quad \alpha = a, b. \quad (18b) $$

Without loss of generality, we can choose $\Theta_{\alpha} = \Theta_{\alpha} e_z$ (see figure 2(b) for the coordinate system), so that the solution is of the form $\mathbf{u} = u_{\theta}(R, \varphi)e_{\theta} \equiv f(R)\sin\varphi e_{\theta}$. The equation (18a) becomes $f''(R) + (2/R)f'(R) + (\kappa^2 - 2/R^2)f(R) = 0$, with the boundary conditions $f(R = \alpha) = \alpha\Theta_{\alpha}$, following (18b). We can find the solution explicitly in terms of the first order spherical Bessel functions of the first $-j_1$—and second $-y_1$—kind:

$$ \begin{aligned} \mathbf{u} &= \begin{pmatrix} j_1(\kappa R) & y_1(\kappa R) \\ j_1(\kappa b) & y_1(\kappa b) \end{pmatrix} \begin{pmatrix} j_1(\kappa a) & y_1(\kappa a) \\ j_1(\kappa b) & y_1(\kappa b) \end{pmatrix}^{-1} \begin{pmatrix} a\Theta_a \\ b\Theta_b \end{pmatrix} \sin \varphi e_{\theta} \\ &= \begin{pmatrix} a \frac{j_1(\kappa R)y_1(\kappa b) - j_1(\kappa b)y_1(\kappa R)}{j_1(\kappa a)y_1(\kappa b) - j_1(\kappa b)y_1(\kappa a)} & b \frac{j_1(\kappa a)y_1(\kappa R) - j_1(\kappa R)y_1(\kappa a)}{j_1(\kappa a)y_1(\kappa b) - j_1(\kappa b)y_1(\kappa a)} \\ \frac{j_1(\kappa R)y_1(\kappa b) - j_1(\kappa b)y_1(\kappa R)}{j_1(\kappa a)y_1(\kappa b) - j_1(\kappa b)y_1(\kappa a)} & \frac{j_1(\kappa a)y_1(\kappa b) - j_1(\kappa b)y_1(\kappa a)}{j_1(\kappa a)y_1(\kappa b) - j_1(\kappa b)y_1(\kappa a)} \end{pmatrix} \begin{pmatrix} \Theta_a \\ \Theta_b \end{pmatrix} \sin \varphi e_{\theta} \end{aligned} \quad (19) $$

As in section 3, we derive the relation between the displacement of the lead sphere and that of the embedding matrix by applying the equation of motion for the coated sphere, and using the fact that the stresses in the cladding can be known through the displacements (equation (19)):

$$ -I_a \omega^2 \Theta_a = \int_{\Gamma_a} \mathbf{r} \times (\boldsymbol{\sigma} \cdot \mathbf{n}) \, dS \quad (20) $$

where the moment of inertia of the hard sphere $I_a = 2M_a a^2/5$ around $e_z$. We finally obtain

$$ \Theta_a = \frac{-\gamma g_*^{(r)}(\omega)}{R_*^{(r)}(\omega) - \rho_a / \rho} \Theta_b := H_*^{(r)}(\omega) \Theta_b, \quad (21) $$
<figure><img src="image_10.png" /><figcaption>Figure 7. Displacement of the spherical core normalized to and exited by that of the matrix, based on the continuum model and discrete mass-spring model, and related to two modes: plane translation and rotation.</figcaption></figure>

provided that

$$ \left\{ \begin{aligned} g_{*}^{(r)}(\omega) &= \frac{5}{\kappa a} \frac{j_1(\kappa a)y_2(\kappa a) - j_2(\kappa a)y_1(\kappa a)}{j_1(\kappa a)y_1(\kappa b) - j_1(\kappa b)y_1(\kappa a)} \\ R_{*}^{(r)}(\omega) &= \frac{5}{\kappa a} \frac{j_2(\kappa a)y_1(\kappa b) - j_1(\kappa b)y_2(\kappa a)}{j_1(\kappa a)y_1(\kappa b) - j_1(\kappa b)y_1(\kappa a)} \end{aligned} \right. \qquad (22) $$

Evolution of the rotational displacement $H_{*}^{(r)}$, and translational displacement of the lead sphere with frequency are shown in figure 7, for a medium with the filling fraction 47%, the lead sphere radius $a = 0.5$ cm, and the thickness of the coating layer is 0.25 cm. We note that, here also, with our material and geometrical parameters, the first two local resonances of the rotational mode occur in lower frequencies compared to those of translational mode, as in the case of cylindrical inclusions. The first resonance frequencies for the translation $f_{0}^{(p)} \approx 380$ Hz and for rotational $f_{0}^{(r)} \approx 261$ Hz. Similarly to the case of cylindrical system, the frequency of these resonances that correspond mainly to the vibration of the core sphere can be estimated through mass-spring models, leading to the following expressions

$$ f_{0}^{(p)} = f_{0} \left\{ 18 \frac{(1 - \nu)(2 - 3\nu)\gamma^2}{(5 - 6\nu)(2 - 3\nu)\left(\frac{\gamma - 1}{\gamma}\right) - \frac{5(\gamma^2 - 1)^2}{4\gamma^5 - 1}} \right\}^{1/2}, \quad f_{0}^{(r)} = f_{0} \left(\frac{15\gamma^5}{\gamma^3 - 1}\right)^{1/2}. \qquad (23) $$

With our specific material and geometrical parameters, the above formulas give $f_{0}^{(p)} \approx 381$ Hz and $f_{0}^{(r)} \approx 272$ Hz, which are quite close to those obtained by the continuum model, in particular for the translation mode. Based on the mass-spring model, the displacement of the hard sphere can be obtained following the calculation of spring constants equivalent of the elastic layer (see equation (B2) in appendix B). Because of the physics of the second resonances in the cladding, we see here also from figure 7, that the mass-spring model ignores these resonances.

The study of the equations (23) leads to the similar conclusions as in section 3. Indeed, as long as $0 \le \nu \le 1/4$, the rotational resonance frequency $f_{0}^{(r)}$ is higher than the translational resonance frequency $f_{0}^{(p)}$. For $1/4 < \nu \le 1/2$, the relative positions of these resonances evolves with the parameter $\gamma = b/a$: for the small size inclusions $\gamma_c(\nu) < \gamma$, we have $f_{0}^{(p)} < f_{0}^{(r)}$, and vice versa otherwise (figure 8). The quantity $\gamma_c(\nu)$ is the unique solution to the equation $f_{0}^{(p)}(\gamma) = f_{0}^{(r)}(\gamma)$, which is an increasing function between 1 for $\nu = 1/4$, and 2.71 when $\nu = 1/2$. Figure 8 shows the evolution of $f_{0}^{(p)}$ with $\gamma$ for three values of $\nu$: 1/8 ($0 \le \nu \le 1/4$), 3/8 ($1/4 < \nu \le 1/2$), and 0.47 that is the value corresponding to the displacement behavior in figure 7.

## 5. Effective-medium properties: metasolids

The knowledge of the kinematics of the cladding is of great interest to go forward in the dynamical homogenization of the whole structure, with the hard coated cylinder or sphere. Indeed, our purpose is now to describe the homogenized behavior of the hard inclusion, the cladding and the matrix; and to obtain the key effective-parameters of the medium.
<figure><img src="image_11.png" /><figcaption>Figure 8. Evolution of the first resonance frequencies normalized by $f_0$ and estimated by the discrete mass-spring model (equation (23)), with the form factor $\gamma = b/a$. The dependence of plane-translation mode on the Poisson's ratio $\nu$ of the coating is also illustrated. The Poisson's ratio takes two representative values 1/8 and 3/8 in the intervals [0, 1/4] and [1/4, 1/2], respectively. The value of $\nu = 0.47$ corresponds to the material parameters studied as an example throughout this paper (see table 1). The rotational mode are independent of $\nu$.</figcaption></figure>

## 5.1. Homogenization for translational modes

For a given translation of the embedding matrix $V_b$, we aim at finding the equivalent density for the embedded solid (sphere or cylinder) and the cladding, subject to the induced translation $V_a$. The knowledge of the stress distribution in the cladding permits us to evaluate the total force applied on the cladding by the matrix, and subsequently to derive the homogenized density. Indeed we can write

$$ -\rho_e \mathcal{V}_e \omega^2 V_b = \int_{\Gamma_b} \sigma \cdot n \, dS \quad (24) $$

where $\mathcal{V}_e$ is the volume of the region occupied by the hard inclusion $\Omega_a$ and the cladding $\Omega$, i.e. the region $\Omega_a \cup \Omega$, that is bounded by the surface $\Gamma_b$. The frequency-dependent effective density $\rho_e$ is defined for each point inside $\Omega_a \cup \Omega$ and must coincide with the static effective value given by $\tilde{\phi}_a \rho_a + \tilde{\phi} \rho$, with $\tilde{\phi}_a$ and $\tilde{\phi}$ being the filling fraction of the coated solid and the cladding, respectively, such that $\tilde{\phi}_a + \tilde{\phi} = 1$.

In order to obtain the effective density for the whole unit cell, it remains to include the density of the matrix $\rho_b$, so that we can finally set $\rho_{eff} = (\phi + \phi_a)\rho_e + \phi_b\rho_b$, where $\phi_a, \phi$ and $\phi_b$ stand for the filling fraction of the coated solid, the cladding and the embedded matrix, respectively, such that $\phi_a + \phi + \phi_b = 1$.

The case of spherical inclusions has been studied in [2], through a result with an analytical expression<fnref n="3" />. Therefore, here we only deal with the cylindrical inclusions, where in particular an anisotropic aspect appears due to the fact that we take into account not only the plane-translation mode, but also the axial mode. Inserting $V_b \equiv U_b + T_b$ in equation (24), the decoupling between the plane and axial translations leads to two different expressions for the effective density $\rho_e$, depending on whether we are dealing with the axial or plane translation. The anisotropy due to the geometrical configuration of the unit cell is obviously the cause for the anisotropy of the effective mass density. We must consider the effective density $\rho_{eff}$ as a second order tensor, and rewrite equation (24) as $-\mathcal{V}_e \omega^2 \rho_{eff} \cdot V_b = \int_{\Gamma_b} \sigma \cdot n \, dS$. Thus, the final expression will be

$$ \rho_{eff}(\omega) = \rho_{eff}^{(p)}(\omega)(e_x \otimes e_x + e_y \otimes e_y) + \rho_{eff}^{(a)}(\omega)e_z \otimes e_z $$

with $(\beta = p, a)$

$$ \rho_{eff}^{(\beta)} = (\phi + \phi_a)\rho_e^{(\beta)} + \phi_b\rho_b, \quad (25) $$

$$ \rho_e^{(\beta)}(\omega) = \rho\{g_1^{(\beta)}(\omega) + g_2^{(\beta)}(\omega)H^{(\beta)}(\omega)\}, \quad (26) $$

$$ \left\{ \begin{aligned} g_1^{(p)}(\omega) &= \begin{pmatrix} [G(hb)] \\ [G(\kappa b)] \end{pmatrix}^T \begin{pmatrix} [E(ha)] & [F(\kappa a)] \\ [E(hb)] & [F(\kappa b)] \end{pmatrix}^{-1} \begin{pmatrix} [0] \\ [1] \end{pmatrix} \\ g_2^{(p)}(\omega) &= \frac{1}{\gamma} \begin{pmatrix} [G(hb)] \\ [G(\kappa b)] \end{pmatrix}^T \begin{pmatrix} [E(ha)] & [F(\kappa a)] \\ [E(hb)] & [F(\kappa b)] \end{pmatrix}^{-1} \begin{pmatrix} [1] \\ [0] \end{pmatrix} \end{aligned} \right. $$

<footnote n="3">³ Note that there is a mistake in the given expressions of $g(\omega)$, $R(\omega)$, $g_1(\omega)$ and $g_2(\omega)$ (equations (24), (25) (32) and (33) in this paper). The terms $(T_{11} + T_{12})j_1(0) + (T_{21} + T_{22})n_1(0)$ involved in $g(\omega)$ and $g_1(\omega)$ and $(T_{13} + T_{14})j_1(0) + (T_{23} + T_{24})n_1(0)$ involved in $R(\omega)$ and $g_2(\omega)$ should be multiplied by $(1 + (\alpha/\beta)^2)$. However, the graphics in figure 2 seem to be correct.</footnote>
<figure><img src="image_12.png" /><figcaption>Figure 9. Effective mass densities based on continuum model and also discrete mass-spring model, associated with plane translation, axial translation, and rotation modes for media with cylindrical inclusions; and with plane translation and rotation modes for media with spherical inclusions.</figcaption></figure>

$$ \left\{ \begin{aligned} g_1^{(a)}(\omega) &= \frac{2}{\kappa b} \frac{J_0(\kappa a) Y_1(\kappa b) - J_1(\kappa b) Y_0(\kappa a)}{J_0(\kappa a) Y_0(\kappa b) - J_0(\kappa b) Y_0(\kappa a)} \\ g_2^{(a)}(\omega) &= \frac{2}{\kappa b} \frac{J_1(\kappa b) Y_0(\kappa b) - J_0(\kappa b) Y_1(\kappa b)}{J_0(\kappa a) Y_0(\kappa b) - J_0(\kappa b) Y_0(\kappa a)} \end{aligned} \right. $$

where the elements of the block matrices are defined in appendix A. Figure 9 shows the evolution of $\rho_{eff}^{(\beta)}(\omega) / \rho_{eff}^{(\beta)}(0)$ with frequency, where we assume that $\rho_{eff}^{(\beta)}(0)$ matches with the static value $\rho_0 \equiv \phi_a \rho_a + \phi \rho + \phi_b \rho_b$. The evolutions for the axial or plane translations are qualitatively similar, but quantitatively different in that $\rho_{eff}^{(a)}$ becomes negative at a lower frequency compared with $\rho_{eff}^{(p)}$ in accordance with the fact that the first resonance occurs for the axial mode (figure 4). We explained in section 3 that the (first) resonance of the axial mode should emerge at lower frequency, independently of the micro-structural parameters. This can be of particular interest, given that it is more challenging to make attenuate low-frequency waves. We note, however, that the resonance band where the mass density becomes negative, may also be less wide for the axial mode. With the present material and geometrical parameters, the band width for the axial mode $\Delta f^{(a)} \approx 118$ Hz, while that of plane-translational mode $\Delta f^{(p)} \approx 280$ Hz.

The effective masses related to the plane translations for media with cylindrical or spherical inclusions, and axial translation modes for the cylindrical system, can also be calculated through the discrete mass-spring model. Once the effective spring constants $k_{eff}^{(p)}$ (plane translation) and $k_{eff}^{(a)}$ (axial translation) are obtained (Appendix B), and the displacements $U_a/U_b$ (plane translation), and $T_a/T_b$ (axial translation) calculated (sections 3 and 4), the equations of motion for the respective mass-spring systems lead easily to the frequency-dependent effective-masses for plane and axial translations. Figure 9 shows these effective masses divided by the volume of the structural unit, in function of frequency. We see that the discrete model predicts well the dynamic effective densities in low frequencies including the first negative region for both plane and axial translation modes. Expectedly, from the analysis on the resonance behaviors in sections 3 and 4, the second negative regions of effective mass densities that correspond to the second resonances (figures 4 and 7), are ignored by the discrete model for both modes and medium micro-geometries (figure 9).

## 5.2. Homogenization for rotational modes

In order to obtain the effective parameter related to the rotations, we apply first the angular equation of motion to the sub-structure $\Omega a \cup \Omega$ formed by the hard inclusion and the cladding. It reads

$$ -I_e \omega^2 \Theta_b = \int_{\Gamma_b} r \times (\sigma \cdot n) dS \qquad (27) $$

where $I_e$ denotes the resulting moment of inertia around $e_z$ for the embedded rigid body and its surrounded cladding that move in phase. Following a similar procedure regarding translation modes, we can establish the expression of the density of the moment of inertia $i_e(\omega)$ such that $dI_e = i_e dV_e$. Since this unusual quantity is obviously intensive in thermodynamic sense, its homogenization in the medium should be possible. Tensorial behavior is involved to ensure the anisotropic aspects of the problem, hence we write
$$ : -\omega^2 V_e i_e \cdot \Theta_b = \int_{\Gamma_b} r \times (\sigma \cdot n) dS. $$
<figure><img src="image_13.png" /><figcaption>Figure 10. Density of moment of inertia based on continuum model and discrete mass-spring model, for media with cylindrical inclusions and spherical inclusions.</figcaption></figure>

Similar to the case of translations, the formal expression of the effective density of the moment of inertia can be written as $i_{eff}(\omega) = (\phi_a + \phi)i_e(\omega) + \phi_b i_b I$.

The spherical system gives rise to an effective density of moment of inertia $i_{eff} \equiv i_{eff} I$, with

$$i_{eff}(\omega) = (\phi_a + \phi)i_e(\omega) + \phi_b i_b, \qquad (28)$$

provided that

$$i_e(\omega) = \frac{2}{5}\rho b^2 \{h_1^*(\omega) + h_2^*(\omega)H_*^{(r)}(\omega)\} \qquad (29)$$

where,

$$\left\{ \begin{aligned} h_1^*(\omega) &= \frac{5}{\kappa b} \frac{j_1(\kappa a) y_2(\kappa b) - j_2(\kappa b) y_1(\kappa a)}{j_1(\kappa a) y_1(\kappa b) - j_1(\kappa b) y_1(\kappa a)} \\ h_2^*(\omega) &= \frac{5}{\gamma \kappa b} \frac{j_2(\kappa b) y_1(\kappa b) - j_1(\kappa b) y_2(\kappa b)}{j_1(\kappa a) y_1(\kappa b) - j_1(\kappa b) y_1(\kappa a)} \end{aligned} \right. $$

The isotropy of $i_{eff}$ is clear and it is evident that $i_{eff}(\omega)$ coincides with the static value $i_0 := \phi_a i_a + \phi_i + \phi_b i_b$ in the quasistatic limit $\omega \to 0$; with $i_a = \frac{2}{5}\rho_a a^2$, $i = \frac{2}{5}\rho \frac{b^5 - a^5}{b^3 - a^3}$, and $i_b = \frac{1}{6}\rho_b \frac{L^5 - 16\pi b^5}{L^3 - 4\pi b^3}$.

Employing the same method that was used in the case of spherical system, for media with cylindrical inclusions, we obtain

$$i_e(\omega) = \frac{1}{2}\rho b^2 \{h_1(\omega) + h_2(\omega)H^{(r)}(\omega)\} \qquad (30)$$

with

$$\left\{ \begin{aligned} h_1(\omega) &= \frac{4}{\kappa b} \frac{j_1(\kappa a) y_2(\kappa b) - j_2(\kappa b) y_1(\kappa a)}{j_1(\kappa a) y_1(\kappa b) - j_1(\kappa b) y_1(\kappa a)} \\ h_2(\omega) &= \frac{4}{\gamma \kappa b} \frac{j_2(\kappa b) y_1(\kappa b) - j_1(\kappa b) y_2(\kappa b)}{j_1(\kappa a) y_1(\kappa b) - j_1(\kappa b) y_1(\kappa a)} \end{aligned} \right. $$

where $i_a = \frac{1}{2}\rho_a a^2$, $i = \frac{1}{2}\rho(b^2 + a^2)$, and $i_b = \frac{1}{6}\rho_b \frac{L^4 - 3\pi b^4}{L^2 - \pi b^2}$.

Figure 10 presents the behavior of $i_{eff}$ in function of frequency, for media with cylindrical and spherical inclusions. This behavior with the negative values near rotational stop-band frequencies is similar to that we observed earlier for the effective mass densities. Here also, after calculating the spring constants for the analogous mass-spring system (appendix B), and then obtaining the displacement $\Theta_a/\Theta_b$ (figures 4 and 7), the effective mass densities can be expressed according to the discrete model. Again, from figure 10 we observe that negative mass densities related to the lower resonance frequency band are appropriately described by the discrete model, while those associated with the higher frequency band, that is largely due to the rotational resonance inside the coating, are disregarded.
## 6. Summary and conclusions

A metasolid composed of a distribution of either cylindrical or spherical hard inclusions coated by a soft material has been studied. Generalizing an existing continuum model [2] that describes resonance-induced band gaps related to plain-translational modes, we provide analytical analysis to investigate anisotropic properties in the case of cylindrical inclusions, and also rotational modes for both cylindrical and spherical systems.

We showed that the translational mode along the cylinders' axis and the rotational mode in both cylindrical and spherical systems could generate band gaps at low frequencies. The band gaps that originate in local resonance phenomena at long-wavelength regime, are tunable in terms of frequency by changing the coating's elastic properties and geometrical parameters of the structural unit. Furthermore, for all translational and rotational modes we found that the discrete mass-spring model could provide appropriate description of the material's kinematics including first resonances, but neglects the second resonances that occur at higher frequencies and mainly arise from resonances inside the coating. The mass-spring model also provides simple analytical expressions allowing us to study easily the kinematics of first resonances in terms of the material parameters. We demonstrated that, while the macroscopic material dynamics is dictated by the effective mass densities for translational modes, it is characterized by its effective moment of inertia regarding the rotational modes. Through analytical expressions, it was shown that components of frequency-dependent anisotropic effective densities for translational modes, and effective density of moment of inertia for rotational mode, become negative near corresponding resonance frequencies. These effective parameters have been also calculated based on the mass-spring model, and their limits have been clarified by comparison to the continuum model.

These results establish a step forward for metasolid characterization via introducing the density of moment of inertia as an effective-medium parameter. They give insight to improving the tunability of the studied metasolid for various applications, such as sound and vibration isolation, and can serve as a physical-based quick guide for optimal design of the material to exhibit desired properties.

## Acknowledgments

The authors acknowledge support from the French National Research Agency (grant no. ANR-13-RMNP-0003-03 and ANR-11-LABX-022-01).

## Appendix A. Expressions of the solutions to Helmholtz equations

Here, we will give the explicit form of the solutions to the Helmholtz equations (4), combined with the boundary conditions for the case of cylindrical inclusions in section 3. The boundary conditions (7) and the solutions (6) lead to the following first conclusions

$$ 
\left\{ 
\begin{array}{l}
C_n = 0 = D_n = G_n = H_n, \forall n \ge 2; \\
A_0 C_0 = 0 = B_0 C_0 \\
E_0 G_0 \text{ and } F_0 G_0 \text{ to be determined;} \\
A_1 C_1, B_1 C_1, E_1 H_1, \text{ and } F_1 H_1 \text{ to be determined;} \\
A_1 D_1, B_1 D_1, E_1 G_1, \text{ and } F_1 G_1 \text{ to be determined.}
\end{array}
\right.
$$

The expressions of the unknowns $E_0 G_0$ and $F_0 G_0$ are established in the article. The coefficients $A_1 C_1, B_1 C_1, E_1 H_1$ and $F_1 H_1$ involved in (8) are solutions to the following linear system

$$ \underbrace{\begin{pmatrix} [E(ha)] & [F(\kappa a)] \\ [E(hb)] & [F(\kappa b)] \end{pmatrix}}_{\mathcal{D}} \begin{pmatrix} A_1 C_1 \\ B_1 C_1 \\ E_1 H_1 \\ F_1 H_1 \end{pmatrix} = \begin{pmatrix} [C_a] \\ [C_b] \end{pmatrix}, \qquad (A1) $$

where the matrices used in the bloc-matrix formalism are

$$ 
\begin{gathered}
[E(x)] = \begin{pmatrix} xJ_1'(x) & xY_1'(x) \\ J_1(x) & Y_1(x) \end{pmatrix}, \quad [F(x)] = \begin{pmatrix} J_1(x) & Y_1(x) \\ xJ_1'(x) & xY_1'(x) \end{pmatrix} \\
\text{and } [C_x] = \begin{pmatrix} xU_x \cos \theta_x \\ xU_x \cos \theta_x \end{pmatrix}.
\end{gathered}
$$
Similarly, the coefficients $A_1 D_1$, $B_1 D_1$, $E_1 G_1$ and $F_1 G_1$ in (8) are determined by solving the linear system:

$$ \begin{pmatrix} [E(ha)] & -[F(\kappa a)] \\ [E(hb)] & -[F(\kappa b)] \end{pmatrix} \begin{pmatrix} A_1 D_1 \\ B_1 D_1 \\ E_1 G_1 \\ F_1 G_1 \end{pmatrix} = \begin{pmatrix} [S_a] \\ [S_b] \end{pmatrix} \qquad (A2) $$

with the same notations for the matrices $[E(x)]$ and $[F(x)]$, whereas $[S_x] = \begin{pmatrix} xU_x \sin \theta_x \\ xU_x \sin \theta_x \end{pmatrix}$. The solutions of the above linear systems can be obtained with bloc-matrix formalism, which leads to the final expression for the displacement field involved in (8).

It is easy now to compute the ratio between the relative plane translations $U_a$ and $U_b$. Using (13), after some formal calculations, we obtain

$$ U_a = H^{(p)}(\omega) U_b = \frac{-\gamma g^{(p)}(\omega)}{R^{(p)}(\omega) - \rho_a / \rho} U_b \qquad (A3) $$

with

$$ \left\{ \begin{aligned} g^{(p)}(\omega) &= \begin{pmatrix} [G(ha)] \\ [G(\kappa a)] \end{pmatrix}^T \begin{pmatrix} [E(ha)] & [F(\kappa a)] \\ [E(hb)] & [F(\kappa b)] \end{pmatrix}^{-1} \begin{pmatrix} [0] \\ [1] \end{pmatrix} \\ R^{(p)}(\omega) &= \begin{pmatrix} [G(ha)] \\ [G(\kappa a)] \end{pmatrix}^T \begin{pmatrix} [E(ha)] & [F(\kappa a)] \\ [E(hb)] & [F(\kappa b)] \end{pmatrix}^{-1} \begin{pmatrix} [1] \\ [0] \end{pmatrix} \end{aligned} \right. \qquad (A4) $$

where $[G(x)] = \begin{pmatrix} J_1(x) \\ Y_1(x) \end{pmatrix}$, $[0] = \begin{pmatrix} 0 \\ 0 \end{pmatrix}$, and $[1] = \begin{pmatrix} 1 \\ 1 \end{pmatrix}$.

We note that the relation (A3) is not only obtained for the plane translation in $e_x$ direction (given after solving (A1)), that is $U_a \cos \theta_a = H^{(p)}(\omega) U_b \cos \theta_b$, but also for the plane translation in $e_y$ direction (given after solving (A2)), i.e., $U_a \sin \theta_a = H^{(p)}(\omega) U_b \sin \theta_b$. Although this seems clear in the physical sense, the mathematical treatment is relatively tedious.

## Appendix B. Expressions for effective stiffnesses according to mass-spring model

The effective stiffnesses based on the discrete mass-spring model are obtained by using the static regime of general equations (2a) for the cladding with $\omega = 0$, and by imposing zero displacement on the inner boundary of the cladding $\Gamma_a$, with $V_a = 0$ and $\Theta_a = 0$, regarding equation (2b). The resulting force gives the spring constants for translations (plane and axial), and the resulting torque leads to the spring constant for rotational motion. The corresponding equations are

$$ \int_{\Gamma_a} \boldsymbol{\sigma}.\boldsymbol{n} \, dS = k_{\text{eff}} V_b \quad \text{and} \quad \int_{\Gamma_a} \boldsymbol{r} \times (\boldsymbol{\sigma}.\boldsymbol{n}) \, dS = C_{\text{eff}} \boldsymbol{\Theta}_b, $$

for translational and rotational motions, respectively. In the above equations $k_{\text{eff}}$ is the spring constant, equivalent of the coating-layer that is represented by a spring in 1D tension/compression, and related to the (plane and axial) translational motion of the layer. Likewise, $C_{\text{eff}}$ is the spring constant describing the torsional stiffness of the coating layer that is taken to be analogous to a spring in 1D tension/compression, related to the rotational motion of the layer. For cylindrical inclusions, the expressions of these constants are obtained as

$$ k_{\text{eff}}^{(p)} = \mu \frac{8\pi L(1 - \nu)(3 - 4\nu)}{(3 - 4\nu)^2 \ln \gamma - \frac{\gamma^2 - 1}{\gamma^2 + 1}}, \quad k_{\text{eff}}^{(a)} = \mu \frac{2\pi L}{\ln \gamma}, \quad C_{\text{eff}}^{(r)} = \mu \frac{4\pi b^2 L}{\gamma^2 - 1}, \qquad (B1) $$

where $k_{\text{eff}}^{(p)}$, $k_{\text{eff}}^{(a)}$, and $C_{\text{eff}}^{(r)}$ refer to plane translation, axial translation, and rotation, respectively. For the spherical inclusions, we obtain

$$ k_{\text{eff}}^{(p)} = \mu \frac{24\pi a(1 - \nu)(2 - 3\nu)}{(5 - 6\nu)(2 - 3\nu)\left(\frac{\gamma - 1}{\gamma}\right) - \frac{5(\gamma^2 - 1)^2}{4\gamma^5 - 1}}, \quad C_{\text{eff}}^{(r)} = \mu \frac{8\pi b^3}{\gamma^3 - 1}. \qquad (B2) $$

The expressions for spring constants $k_{\text{eff}}^{(p)}$ and $C_{\text{eff}}^{(r)}$ for media with cylindrical and spherical inclusions can also be found in [32].

## ORCID iDs

Navid Nemati https://orcid.org/0000-0002-1370-0677
References

[1] Liu Z, Zhang X, Mao Y, Zhu Y Y, Yang Z, Chan C T and Sheng P 2000 Locally resonant sonic materials *Science* **289** 1734

[2] Liu Z, Chan C T and Sheng P 2005 Analytic model of photonic crystals with local resonances *Phys. Rev. B* **71** 014103

[3] Yang Z, Mei J, Yang M, Chan N H and Sheng P 2008 Membrane-type acoustic metamaterial with negative dynamic mass *Phys. Rev. Lett.* **101** 204301

[4] Fang N, Xi D, Xu J, Ambati M, Srituravanich W, Sun C and Zhang X 2006 Ultrasonic metamaterials with negative modulus *Nat. Mater.* **5** 452

[5] Nemati N, Kumar A, Lafarge D and Fang N X 2015 Nonlocal description of sound propagation through an array of helmholtz resonators *C. R. Mecanique* **343** 656

[6] Li J and Chan C T 2004 Double-negative acoustic metamaterial *Phys. Rev. E* **70** 055602(R)

[7] Ding Y, Liu Z, Qiu C and Shi J 2007 Metamaterial with simultaneously negative bulk modulus and mass density *Phys. Rev. Lett.* **99** 093904

[8] Cheng Y, Xu J Y and Liu X J 2008 One-dimensional structured ultrasonic metamaterials with simultaneously negative dynamic density and modulus *Phys. Rev. B* **77** 045134

[9] Lee S H, Park C M, Seo Y M, Wang Z G and Kim C K 2010 Composite acoustic medium with simultaneously negative density and modulus *Phys. Rev. Lett.* **104** 054301

[10] Liu X N, Hu G K, Huang G L and Sun C T 2011 An elastic metamaterial with simultaneously negative mass density and bulk modulus *Appl. Phys. Lett.* **98** 251907

[11] Yang M, Ma G, Yang Z and Sheng P 2013 Coupled membranes with doubly negative mass density and bulk modulus *Phys. Rev. Lett.* **110** 134301

[12] Wu Y, Lai Y and Zhang Z Q 2011 Elastic metamaterials with simultaneously negative effective shear modulus and mass density *Phys. Rev. Lett.* **107** 105506

[13] Lai Y, Wu Y, Sheng P and Zhang Z Q 2011 Hybrid elastic solids *Nat. Mater.* **10** 620

[14] Brunet T, Merlin A, Mascaro B, Zimny K, Leng J, Poncelet O, Aristégui C and Mondain-Monval O 2015 Soft 3D acoustic metamaterial with negative index *Nat. Mater.* **14** 384

[15] Sigalas M and Economou E 1993 Band structure of elastic waves in two-dimensional systems *Solid State Commun.* **86** 141

[16] Kushwaha M S, Halevi P, Dobrzynski L and Djafari-Rouhani B 1993 Acoustic band structure of periodic elastic composites *Phys. Rev. Lett.* **71** 2022

[17] Liang Z and Li J 2012 Extreme acoustic metamaterial by coiling up space *Phys. Rev. Lett.* **108** 114301

[18] Xie Y, Popa B I, Zigoneanu L and Cummer S A 2013 Measurement of a broadband negative index with space-coiling acoustic metamaterials *Phys. Rev. Lett.* **110** 175501

[19] Milton G W and Willis J R 2007 On modifications of Newton's second law and linear continuum elastodynamics *Proc. R. Soc. A* **463** 855

[20] Cummer S A and Schurig D 2007 One path to acoustic cloaking *New J. Phys.* **9** 45

[21] Ao X and Chan C T 2008 Far-field image magnification for acoustic waves using anisotropic acoustic metamaterials *Phys. Rev. E* **77** 025601(R)

[22] Li J, Fok L, Yin X, Bartal G and Zhang X 2009 Experimental demonstration of an acoustic magnifying hyperlens *Nat. Mater.* **8** 931

[23] Zhu J, Christensen J, Jung J, Martin-Moreno L, Yin X, Fok L, Zhang X and Garcia-Vidal F G 2010 A holey-structured metamaterial for acoustic deep-subwavelength imaging *Nat. Phys.* **7** 52

[24] Zhu X, Liang B, Kan W, Zou X and Cheng J 2011 Acoustic cloaking by a superlens with single-negative materials *Phys. Rev. Lett.* **106** 014301

[25] Chen Y, Liu H, Reilly M, Bae H and Yu M 2014 Enhanced acoustic sensing through wave compression and pressure amplification in anisotropic metamaterials *Nat. Comm.* **5** 5247

[26] Torrent D, Pennec Y and Djafari-Rouhani B 2015 Resonant and nonlocal properties of phononic metasolids *Phys. Rev. B* **92** 174110

[27] Zhao H, Liu Y, Wang G, Wen J, Yu D, Han X and Wen X 2005 Resonance modes and gap formation in a two-dimensional solid phononic crystal *Phys. Rev. B* **72** 012301

[28] Peng P, Mei J and Wu Y 2012 Lumped model for rotational modes in phononic crystals *Phys. Rev. B* **86** 134304

[29] Bigoni D, Guenneau S, Movchan A B and Brun M 2013 Elastic metamaterials with inertial locally resonant structures: application to lensing and localization *Phys. Rev. B* **87** 174303

[30] Zhu R, Liu X N, Hu G K, Sun C T and Huang G L 2014 Negative refraction of elastic waves at the deep-subwavelength scale in a single-phase metamaterial *Nat. Commun.* **5** 5510

[31] Peng P, Asiri S, Zhang X, Li Y and Wu Y 2013 A lumped model for rotational modes in periodic solid composites *EPL* **104** 26001

[32] Bonnet G and Monchiet V 2015 Low frequency locally resonant metamaterials containing composite inclusions *J. Acoust. Soc. Am.* **137** 3263

[33] Love A E H 1944 *A Treatise on the Mathematical Theory of Elasticity* (New York: Dover)

[34] Bonnet G and Monchiet V 2017 Dynamic mass density of resonant metamaterials with homogeneous inclusions *J. Acoust. Soc. Am.* **142** 890

[35] Huang H H and Sun C T 2009 Wave attenuation mechanism in an acoustic metamaterial with negative effective mass density *New J. Phy.* **11** 013003