# Design and experimental validation of a temperature-driven adaptive phononic crystal slab

Kévin Billon, Morvan Ouisse, Emeline Sadoulet, Manuel Collet, Pauline

Butaud, Gael Chevallier, Abdelkrim Khelif

► To cite this version:

Kévin Billon, Morvan Ouisse, Emeline Sadoulet, Manuel Collet, Pauline Butaud, et al.. Design and experimental validation of a temperature-driven adaptive phononic crystal slab. Smart Materials and Structures, 2019, 28 (3), 10.1088/1361-665X/aaf670. [hal-02160947](https://hal.science/hal-02160947)

HAL Id: hal-02160947

https://hal.science/hal-02160947

Submitted on 20 Jun 2019

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
# Design and experimental validation of a temperature-driven adaptive phononic crystal slab

K. Billon<sup>a</sup>, M. Ouisse<sup>a</sup>, E. Sadoulet-Reboul<sup>a</sup>, M. Collet<sup>c</sup>, P. Butaud<sup>a</sup>, G. Chevallier<sup>a</sup>, and A. Khelif<sup>b</sup>

## ABSTRACT

In this paper, an adaptive phononic crystal slab based on the combination of metallic parts and highly dissipative polymeric interfaces is designed. Cylindrical pillars are composed of shape memory polymer and aluminum deposited periodically on the aluminum slab. The mechanical properties of the polymer depend on both temperature and frequency and can radically change from glassy to rubbery state, with various combinations of high/low stiffness and high/low dissipation. A 3D finite element model of the cell is developed for the design of the metamaterial. The shifted-cell operator technique is used to correctly handle damping effects in the dispersion analysis. In order to validate the design and the adaptive character of the metamaterial, results issued from a full 3D model of a finite structure embedding an interface composed by a distributed set of the unit cells are presented. Various driving temperatures are used to change the behaviour of the system, and numerical results obtained on the adaptive structure are compared to experimental ones. Two states are obtained by changing the temperature of the polymeric interface: at 25°C a bandgap is visible around a selected resonance frequency, and it doesn't exist anymore above the glass transition temperature, where the phononic crystal slab tends to behave as an homogeneous plate. Numerical and experimental results show energy propagation along the borders of the slab in the bandgap.

**Keywords:** Periodic structures, Metamaterial, Dispersion, Dissipation, Vibroacoustics.

## 1. INTRODUCTION

A metamaterial is an artificial multiscale architECTured material designed to control and manipulate waves in gases, liquids or solids. This paper focuses on elastic wave propagations. In the case of elastic waves, metamaterials can exhibit particular elastodynamic properties that are not found in a natural material. The term appeared in 1999. However, physics governing their behavior was developed in the 1960s by physicist Viktor Veselago.¹ It was not until the 2000s that the first experimental realization<sup>2,3</sup> was achieved. Mechanical metamaterials have been recently hailed as a new class of structural concept able to bring novel multifunctionalities<sup>4</sup> by changes of compliance, shape, or by embedding oscillators or smart material inserts. Some examples (to name a few) are multiscale architecturally structured topologies,<sup>5</sup> zig-zag folded sheets,<sup>6</sup> pentamodal lattices,<sup>7</sup> systems with distributed resonators,<sup>8,9</sup> smart/magnetic materials,<sup>10</sup> tunable connectivity,<sup>11</sup> phononic stubbed plates,<sup>12</sup> flat lenses (super lenses),<sup>13</sup> 3D tunable phononic crystals,<sup>14</sup> tunable metamaterial beam with shape memory alloy resonators,<sup>15,16</sup> auxetic periodic structure (with negative Poisson ratio)<sup>17,18</sup> and nonlinear auxetic dampers.<sup>19</sup>

<footnote><sup>a</sup>Univ. Bourgogne Franche-Comté FEMTO-ST Institute CNRS/UFC/ENSMM/UTBM, Department of Applied Mechanics, 24 rue de l'Epitaphe, 25000 Besançon, France</footnote>

<footnote><sup>b</sup>Univ. Bourgogne Franche-Comté FEMTO-ST Institute CNRS/UFC/ENSMM/UTBM, Department of Micro Nano Sciences & Systems, 15B Avenue des Montboucons, 25000 Besançon, France</footnote>

<footnote><sup>c</sup>Laboratory of Tribology and Systems Dynamics (LTDS), École Centrale de Lyon, 36 avenue Guy de Collongue, 69134 Ecully, France</footnote>

<footnote>Further author information: (Send correspondence to Morvan Ouisse)
E-mail: morvan.ouisse@femto-st.fr, Telephone: +33 3 81 66 60 00</footnote>
A phononic crystal is a metamaterial with a periodic structure, that exhibits spatial periodicity. Historically, they have been introduced in order to extend to the domain of elastic waves phenomena highlighted during the propagation of electromagnetic waves in photonic crystals. When photonics and phononics effects cross and there is a change in the propagation of elastic and optical waves, the term phoXonic crystal is used.²⁰

The control of the elastic waves can be performed by combining Bragg's bandgap²¹ (wave interferences), resonant's bandgap (resonance of a component embedded in the unit cell), damping and/or active control. The first occurs when the wavelength is near with the characteristic length of the periodic network. The local resonance phenomenon occurs when the frequency of the wave corresponds to the resonance frequency of the resonator. The periodic structures can stop the wave propagation, but can also lead to other properties such as confining waves or guiding waves in a particular direction or guiding along a chosen path of propagation. These phenomena occur at a longer wavelength than the crystal period. The energy can then be reflected, transmitted, damped, focused or confined in a specific zone of the structure. However, the practical realization of real-life 2D or 3D finite systems may lead to some situations where energy transfers are not in accordance with those predicted by the infinite system considered in the design, because of reflections on the boundaries conditions of the finite structure.

On the other hand, the study of periodic structures has a long history in the field of vibrations and acoustics.²² This topic has interested researchers over the years, and a growing activity on this field is observed on the last decade, with the objective of designing structures exhibiting properties that conventional ones cannot possess.¹⁷, ²³, ²⁴ The methods currently used are most of the time based on those derived from wave propagation in crystals,²⁵ where almost no dissipation occurs. Reaching the upper scale for structural dynamics implies that damping effects have to be included in the analysis.

In this paper, first, some numerical tools for dispersion analysis of periodic structures are presented. The classical Floquet-Bloch approach is presented, as a reference. This technique uses proper boundary conditions on the unit cell, but dealing with damping is not easy for 2D or 3D cases. Secondly, the shifted-cell operator technique is described. It consists of a reformulation of the partial differential equation (PDE) problem by shifting -in terms of wavenumber- the space derivatives appearing in the mechanical behavior operator inside the cell, while imposing continuity boundary conditions on the borders of the domain. Damping effects can be introduced in the system. This strategy makes it possible to solve the problem with an arbitrary frequency dependency of the physical properties of the cell. A focus is proposed on tools for the post-processing of dispersion diagrams in damped configurations, e.g., group velocity. Third, an adaptive metamaterial based on the combination of metallic parts with a highly dissipative polymeric interface is designed. In order to validate the design and the adaptive character of the metamaterial, results issued from a full 3D model of a finite structure embedding an interface made of a distributed set of the unit cells are presented. After this step, a comparison of the results obtained using the tunable structure simulation and the experimental results is presented. Finally, a supercell-based approach is proposed to handle finite system boundary conditions in order to be able to identify situations in which energy transfer may arise because of reflections on the border of the elastic domain. Computations are performed on a two-cell with adequate boundary conditions. The methodology is described and validated using the full finite model and experimental tests on a 2D metamaterial structure.

## 2. NUMERICAL TOOLS FOR EFFICIENT SIMULATIONS OF WAVE PROPAGATION IN DAMPED PERIODIC STRUCTURES

### 2.1 Reference structure

The reference structure used in this paper is presented in reference.²⁶ The system consists of an infinite 3D periodic bidirectional waveguide shown in figure 1a. It is a 1 mm thick plate with periodic cylindrical pillars made of isotropic aluminium 6063-T83 ($\nu = 0.33$, $E = 69e9[Pa]$ and $\rho = 2700[kg/m^3]$). The unit cell and the corresponding first Brillouin zone are shown in figure 1b.

### 2.2 Classical method based on Floquet-Bloch theorem

The Floquet-Bloch approach is a method commonly used for the analysis of periodic structures. The material is supposed to be linear, elastic and isotropic. For a 2D periodic structure, the harmonic homogeneous dynamical equilibrium of the system is characterised by the following partial differential equation
<figure><img src="image_2.png" /><figcaption>Figure 1. a) Infinite reference wave guide. b) Real and reciprocal lattices.</figcaption></figure>

$$ \rho(\mathbf{x})\omega^2\mathbf{u}(\mathbf{x}) + \nabla[\mathbf{C}(\mathbf{x})\cdot\nabla_{\text{sym}}(\mathbf{u}(\mathbf{x}))] = 0, \quad (1) $$

where $\mathbf{u}(\mathbf{x}) \in \mathbb{R}^3$ is the displacement vector, $\mathbf{C}(\mathbf{x})$ is the elastic tensor, $\varepsilon(\mathbf{x}) = \nabla_{\text{sym}}(\mathbf{u}(\mathbf{x})) = \frac{1}{2}(\nabla\mathbf{u}^T(\mathbf{x}) + \mathbf{u}(\mathbf{x})\nabla^T)$ is the strain tensor, $\rho$ is the material density and $\omega$ is the circular frequency.

The periodicity is defined on the borders of the domain using the Floquet-Bloch boundary conditions, that write $\mathbf{u}_R = e^{-j\mathbf{k}_x r_1}\mathbf{u}_L$ and $\mathbf{v}_R = e^{-j\mathbf{k}_y r_2}\mathbf{v}_L$ where $\mathbf{u}_R$ (resp. $\mathbf{v}_R$) is the displacement on the right border and $\mathbf{u}_L$ (resp. $\mathbf{v}_L$) is the displacement on the left border in $x$ (resp. $y$) axis, $\mathbf{k}_x$ and $\mathbf{k}_y$ are respectively the wavenumbers in the $x$ and $y$ directions$^{27,28}$ as shown in figure 1b. A parametric eigenvalue analysis is performed using the Pardiso solver,$^{29}$ two parameters (wavenumbers) are considered, namely $k_x = [0 \pi/r_1]$ and $k_y = [0 \pi/r_2]$. The eigenfrequencies are obtained by solving the problem. The wave's dispersion curves of the undamped system are plotted on the whole first Brillouin zone (figure 2a) and on the contours (figure 2b). The frequencies defining the bandgaps can be found by considering only the contour of the irreducible Brillouin zone for non-singular system.$^{25}$ It is worth noticing that this property is extensively used in open literature, although no formal proof of its validity is available, and therefore the results obtained need to be checked carefully.$^{30}$ So, figure 2b is sufficient to observe partial bandgaps along $x$ direction and complete bandgap in all the directions. It appears that the structure exhibits partial bandgaps for frequencies around 20 and 100 kHz in some specific directions and a complete bandgap between 114 and 143 kHz. In this formulation, the wavenumbers are parameters of the eigenvalue problem, the solutions being the frequencies. Therefore, handling frequency-dependency of the mechanical properties is not easy, since it requires the resolution of a nonlinear and non polynomial eigenvalue problem, or condensation procedures providing ill-conditioned matrices.$^{31,32}$ Alternative procedures can be found in the literature, such as ref.,$^{33}$ where the damping provided by a generalized Maxwell model is included in the stiffness matrix, and ref.,$^{34}$ where an EBSM (Extended Bloch Mode Synthesis) with modal reduction is applied for fast calculation. The next section recalls a suitable method for handling frequency-dependency in dispersion analysis.
<figure><img src="image_3.png" /><figcaption>Figure 2. a) Dispersion surfaces on the Brillouin zone. b) Dispersion curves on the contour of the Brillouin zone. Partial bandgaps and complete bandgap respectively in blue and green.</figcaption></figure>

## 2.3 Shifted-cell operator method

The shifted-cell operator, alternative formulation proposed by,³⁵⁻³⁷ is recalled in this subsection with details on its numerical implementation and the link between left and right eigenvalue problems.

### 2.3.1 Numerical implementation

The shifted-cell operator³⁷ technique consists of a reformulation of the PDE problem by shifting -in terms of wavenumber- the space derivatives appearing in the operators inside the cell, while imposing continuity boundary conditions on the borders of the domain. The formulation leads to the following eigenvalue problem :

$$ [(\mathbf{K} - \omega^2 \mathbf{M}) + \lambda_i(\mathbf{L} - \mathbf{L}^T) - \lambda_i^2 \mathbf{H}] \phi_i^r = 0, \quad (2) $$

where $\lambda_i = jk_i$ is the $i-th$ eigenvalue, $\phi_i^r$ is the right eigenvector associated to $\lambda_i$, $\mathbf{M}$ and $\mathbf{K}$ are respectively the standard symmetric definite mass and symmetric semi-definite stiffness matrices, $\mathbf{L} - \mathbf{L}^T$ is a skew-symmetric matrix and $\mathbf{H}$ is a symmetric semi-definite positive matrix for details). In this formulation, all matrices can depend on $\omega$. A parametric eigenvalue analysis is then performed where the pulsation $\omega$ and the wave propagation angle $\theta$ (such that $k_x = k \cos\theta$ and $k_y = k \sin\theta$) are fixed as real parameter, allowing introduction of damping effects. The wavenumbers $\lambda_i = jk_i$ and the associated right eigenvectors $\phi_i^r$ are computed by solving the quadratic eigenvalue problem (equation 2).

For the computation of the group velocity that will be described in the next section, both left and right eigenvectors of the non-symmetric eigenvalue problem are required. We propose to rewrite (equation 2), using the state-space, as

$$ \left( \begin{array}{cc} 0 & \mathbf{I}_d \\ (\mathbf{K} - \omega^2 \mathbf{M}) & \mathbf{L} \end{array} \right) \left( \begin{array}{c} \phi_i^r \\ \lambda_i \phi_i^r \end{array} \right) = \lambda_i \left( \begin{array}{cc} \mathbf{I}_d & 0 \\ \mathbf{L}^T & \mathbf{H} \end{array} \right) \left( \begin{array}{c} \phi_i^r \\ \lambda_i \phi_i^r \end{array} \right), \quad (3) $$

which is denoted

$$ \mathbf{A}_1(\omega)\psi_i^r = \lambda_i \mathbf{A}_2(\omega)\psi_i^r, \quad (4) $$

where
$$ \mathbf{A}_1(\omega) = \begin{pmatrix} 0 & \mathbf{I}_d \\ (\mathbf{K} - \omega^2\mathbf{M}) & \mathbf{L} \end{pmatrix}, \qquad (5) $$

$$ \mathbf{A}_2(\omega) = \begin{pmatrix} \mathbf{I}_d & 0 \\ \mathbf{L}^T & \mathbf{H} \end{pmatrix}, \qquad (6) $$

and

$$ \psi_i^r = \begin{pmatrix} \phi_i^r \\ \lambda_i \phi_i^r \end{pmatrix}. \qquad (7) $$

Conversely, a left-eigenvector for the same eigenvalue satisfies

$$ \psi_i^{lT} \mathbf{A}_1(\omega) = \lambda_i \psi_i^{lT} \mathbf{A}_2(\omega). \qquad (8) $$

Equation 8 can be developed to find the link between $\psi_i^r$ and $\psi_i^l$. After a few steps the following relationship is obtained:

$$ \psi_i^l = \begin{pmatrix} \lambda_i(\mathbf{H}^T \phi_i^l) - \mathbf{L}^T \phi_i^l \\ \phi_i^l \end{pmatrix} = \begin{pmatrix} \lambda_i(\mathbf{H}^T \phi_{-i}^r) - \mathbf{L}^T \phi_{-i}^r \\ \phi_{-i}^r \end{pmatrix}. \qquad (9) $$

Hence, by solving the right eigenvalue problem, the i-th mode ($i \in \mathbb{N}^*$) is obtained, $\lambda_i \in \mathbb{C}$ being the eigenvalue associated to the right eigenvector $\phi_i^r$. For each mode $i$, a mode $-i$ is associated with $\lambda_{-i}$ such that $\lambda_{-i} = -\lambda_i$ and $\phi_{-i}^r$ such that $\phi_{-i}^r = \phi_i^l$. So, by solving the right eigenvalue problem, the adjoint solution is found too.

### 2.3.2 Group velocity

In this section, the expression of the group velocity is derived. It will be used in the next section, on the one hand, as a sorting criterion for distinguishing a propagating wave from an evanescent wave and, on the other hand, as an indicator for branch tracking.

For frequency-dependent systems, the estimation of the group velocity is not trivial.³⁸,³⁹

Equation 4 is differentiated and multiplied by the left eigenvector $\psi_i^{lT}$ following the procedure proposed by³⁸ such that

$$ \psi_i^{lT} \frac{\partial \mathbf{A}_1(\omega)}{\partial \omega} \psi_i^r + \psi_i^{lT} \mathbf{A}_1(\omega) \frac{\partial \psi_i^r}{\partial \omega} = \psi_i^{lT} \frac{\partial \lambda_i}{\partial \omega} \mathbf{A}_2(\omega) \psi_i^r + \psi_i^{lT} \lambda_i \frac{\partial \mathbf{A}_2(\omega)}{\partial \omega} \psi_i^r + \psi_i^{lT} \lambda_i \mathbf{A}_2(\omega) \frac{\partial \psi_i^r}{\partial \omega}. \qquad (10) $$

According to equation 8,

$$ \psi_i^{lT} \mathbf{A}_1(\omega) \frac{\partial \psi_i^r}{\partial \omega} = \psi_i^{lT} \lambda_i \mathbf{A}_2(\omega) \frac{\partial \psi_i^r}{\partial \omega}, \qquad (11) $$

hence

$$ \psi_i^{lT} \frac{\partial \mathbf{A}_1(\omega)}{\partial \omega} \psi_i^r = \psi_i^{lT} \frac{\partial \lambda_i}{\partial \omega} \mathbf{A}_2(\omega) \psi_i^r + \psi_i^{lT} \lambda_i \frac{\partial \mathbf{A}_2(\omega)}{\partial \omega} \psi_i^r. \qquad (12) $$

After several steps, isolating the $\frac{\partial \lambda_i}{\partial \omega}$ expression, the following formulation is obtained:

$$ \frac{\partial \lambda_i}{\partial \omega} = \frac{\psi_i^{lT} \left[ \frac{\partial \mathbf{A}_1(\omega)}{\partial \omega} - \lambda_i \frac{\partial \mathbf{A}_2(\omega)}{\partial \omega} \right] \psi_i^r}{\psi_i^{lT} \mathbf{A}_2(\omega) \psi_i^r}. \qquad (13) $$
In this case, the group slowness (see section 6 for details) is expressed as

$$S_g = \frac{\partial k_i}{\partial \omega} = -j \frac{\phi_i^{LT} [-2\omega M + \frac{\partial K}{\partial \omega} - \lambda_i \frac{\partial L^T}{\partial \omega} + \lambda_i \frac{\partial L}{\partial \omega} - \lambda_i^2 \frac{\partial H}{\partial \omega}] \phi_i^r}{\phi_i^{LT} [-L + L^T + 2\lambda_i H] \phi_i^r}. \quad (14)$$

The complex group velocity is the inverse of the group slowness

$$C_g = \frac{1}{\frac{\partial k_i}{\partial \omega}} = \frac{1}{-j \frac{\phi_i^{LT} [-2\omega M + \frac{\partial K}{\partial \omega} - \lambda_i \frac{\partial L^T}{\partial \omega} + \lambda_i \frac{\partial L}{\partial \omega} - \lambda_i^2 \frac{\partial H}{\partial \omega}] \phi_i^r}{\phi_i^{LT} [-L + L^T + 2\lambda_i H] \phi_i^r}}. \quad (15)$$

This expression is consistent with the estimation of the energy velocity $v = I/E^{40}$ where I is the flow of energy and E is the total mechanical energy density.

### 2.3.3 Sorting criteria for distinguishing a propagating wave from an evanescent wave

The problem being damped, all the wavenumbers are complex, and the distinction between "propagative" and "evanescent" waves becomes difficult. This is why sorting criteria are proposed in the following.

First of all, all the waves are shifted to the first irreducible Brillouin zone. Then, criteria are defined as:

* ratio between real and imaginary part of each wavenumbers

$$C_1 = real(k)/imag(k); \quad (16)$$

* ratio between real and imaginary part of the velocity of energy transport$^{40}$ $v = I/E$ where I is the flow of energy and E the total mechanical energy density, approximated from the kinetic energy $E = 2E_c$

$$C_2 = real(v)/imag(v); \quad (17)$$

* ratio between real and imaginary part of the group velocity$^{38}$

$$C_3 = real(C_g)/imag(C_g). \quad (18)$$

Only the waves corresponding to $C_1 > \tau_1$; $C_2 > \tau_2$ et $C_3 > \tau_3$ are considered propagative. In practice, the thresholds $\tau_1$, $\tau_2$ et $\tau_3$ are chosen such as $\tau_1 = \tau_2 = 1$; $\tau_3 = 2$. This is an arbitrary choice that provides good results for the cases investigated by the authors. However, for others cases, alternative values of the thresholds may be required.

### 2.3.4 Shifted-cell operator implementation validation for conservative case

The reference structure (figure 1a) is used. Figure 3 presents dispersion curves along the $\Gamma - X$ ($\phi = 0^\circ$) direction obtained with the shifted-cell operator method for conservative case with raw results in black (unsorted results) and sorted results in blue, the real part of wavenumbers is represented on the line from 0 to 1 and the imaginary part of wavenumbers from 0 -1. Sorted results are obtained after the application of sorting criteria describe above. No damping being used in the simulations, consequently, the propagating branches possess purely real wavenumbers hence the blue points along the ordinate axis ($imag(k) = 0$).

A comparison between the results obtained using the method based on the Floquet-Bloch theorem and the shifted-cell operator is first performed on a conservative structure to validate the implementation of the shifted-cell operator technique. Figure 4a presents the dispersion curves obtained along the $\Gamma - X$ direction when $\phi$ is equal to zero. Shifted-cell operator and Floquet-Bloch results are respectively in blue dots and dotted lines. Both methods lead to similar results and the shifted-cell operator method is thus verified.

It can be noted that around 70 kHz on the red dotted lines, the branches seem to get closer without crossing, there is indeed a crossing from which a method of branch tracking is developed in the next subsection 2.3.5.
<figure><img src="image_4.png" /><figcaption>Figure 3. Dispersion curves along the $\Gamma - X$ direction obtained with the shifted-cell operator method for conservative case, real part of wavenumbers is represented on the line from 0 to 1 and the imaginary part of wavenumbers from 0 -1, unsorted and sorted results respectively in black and blue. Grey shapes represent bandgaps.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Dispersion curves obtained with the method based on the Floquet-Bloch theorem and the shifted-cell operator technique for conservative case along the $\Gamma - X$ direction. Grey shapes represent bandgaps.</figcaption></figure>
### 2.3.5 Branch tracking

The dispersion diagrams being obtained point by point, the tracking of the branches is required to identify automatically the curves corresponding to each physical mode. The group velocity constitutes a pertinent indicator for the branch tracking from one computational point to another. In order to do that, an algorithm allows to compare the computed group velocity ($\frac{\partial\omega}{\partial k_i}$) computed by equation 34 and the post-processed group velocity ($\frac{\Delta\omega}{\Delta k_i}$) which is evaluated using finite difference.

Solutions for associating points on the basis of the nature of the displacement field to identify the branches are proposed in the literature, MAC criteria or eigenvectors orthogonality.³⁷ These methods need to store the eigenvectors at each step, which may correspond to a lot of data. An alternative is proposed here on the basis of the group velocity value. This indicator ensures the correct plotting of the dispersion diagram, particularly in the case of crossing branches, veering or bifurcation.⁴¹

The indicator operates from left to right (figure 5) : starting from a $(\omega, k)$ point with known group velocity ($\frac{\partial\omega}{\partial k_i}$), the routine compares this velocity with all the post-processed group velocities by finite difference of the right-hand points. A minimization is done to identify the point with the post-processed group velocity closest to ($\frac{\partial\omega}{\partial k_i}$). This point is defined as next point and step by step the branch is identified. The sorting algorithm has been validated with damping.

<figure><img src="image_6.png" /><figcaption>Figure 5. Scheme of the branch tracking routine</figcaption></figure>

For example, in figure 6, around 70 kHz, there is a cross between the light blue branch (n°6) and the green branch (n°5). The branches on this diagram are determined from the figure 7, showing the evolution of the group velocity with a color code corresponding to the dispersion diagram of the figure 6. It can be seen that the group velocity goes to zero near the upper or the lower frequency of the band gaps, meaning that the energy distribution becomes stationary.
<figure><img src="image_7.png" /><figcaption>Figure 6. Dispersion curves along the Γ − X direction obtained with the shifted-cell operator method for conservative case with branch tracking.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Group velocity associated to dispersion curves (figure 6)</figcaption></figure>

3. DESIGN AND VALIDATION OF A TEMPERATURE-DRIVEN ADAPTIVE PHONONIC CRYSTAL SLAB

3.1 Description of the concept

The structure presented in section 2.1 is modified in order to obtain an adaptive structure. As shown in figure 8, the pillars are now made of combination between a highly dissipative polymer tBA/PEGDMA 42 and aluminum 6063-T83.

As it will be shown in the next section, the polymer thickness has been chosen in order to open a reso- nant bandgap around 40 kHz, which means below the Bragg bandgap frequency. The polymer has a density ρ$_{poly}$ = 1004[kg/m$^{3}$] and a Poisson ratio ν$_{poly}$ = 0,37. The elastic properties are highly dependent on the temperature and the frequency.

In the harmonic regime, the complex modulus of the polymer writes

$$
E_{\text{poly}}^{*}(\omega, T) = E' + jE'' = E'(1 + j\tan(\delta)), \qquad (19)
$$

where tan( δ ) = E ′′ /E ′ is the loss factor, E ′ is the storage modulus and E ′′ is the loss modulus. In this work, a fractional derivative Zener model is used. The expression of the elastic complex modulus is
<figure><img src="image_9.png" /><figcaption>Figure 8. a) Reference structure, 10 mm height aluminum cylinder. b) Metamaterial with 3 mm high polymeric interface between base and 7 mm high aluminum cylinder.</figcaption></figure>

<table><caption>Table 1. tBA/PEGDMA Zener parameters</caption><thead><tr><th>E<sub>0,poly</sub> [MPa]</th><th>E<sub>∞,poly</sub> [Mpa]</th><th>α</th><th>τ<sub>0</sub></th></tr></thead><tbody><tr><td>0.67</td><td>2211</td><td>0.78</td><td>0.91</td></tr></tbody></table>

$$E_{poly}^{*}(\omega, T) = \frac{E_{0,poly} + E_{\infty,poly}(j\omega\tau(T))^{\alpha}}{1 + (j\omega\tau(T))^{\alpha}}, \qquad (20)$$

where $E_{0,poly}$ et $E_{\infty,poly}$ are respectively the static elastic modulus and the asymptotic value of the dynamical modulus, $\tau$ is the relaxation time and $\alpha$ is the order of fractional derivative. The values of the four parameters $E_{0,poly}$, $E_{\infty,poly}$, $\alpha$ and $\tau$ is obtained by the fitting of experimental measurements$^{43,44}$ (tableau 1).

The translation factor $a_T$ comes from the Williams-Landel-Ferry law (WLF),$^{45}$ the temperature dependency of $a_T$ is described by

$$\log_{10}(a_T(T)) = \frac{-C_1^0(T - T_0)}{C_2^0 + (T - T_0)}, \qquad (21)$$

with $C_1^0$ et $C_2^0$ the law constants for reference temperature $T_0 = 40^{\circ}C$ ($C_1^0 = 10,87$ and $C_2^0 = 32,57 K$). The relaxation time $\tau$ is linked to the translation factor :

$$\tau(T) = a_T(T)\tau_0. \qquad (22)$$

The storage modulus and the loss factor are plotted in figure 9a and figure 9b for two temperatures of interest, namely $25^{\circ}C$ and $90^{\circ}C$. At ambient temperature, the polymer is stiff and its loss factor is almost equal to zero. At $90^{\circ}C$, the loss factor is over 1.5 on the frequency range (0 to 100kHz), with a maximum value of 2.5 and the stiffness is very low.
<figure><img src="image_10.png" /><figcaption>Figure 9. Mechanical properties of the tBA/PEGDMA at 25°C and 90°C. a) Storage modulus. b) Loss factor.</figcaption></figure>

## 3.2 Dynamical properties

### 3.2.1 Dispersion analysis

The shifted-cell operator technique is used to obtain the dispersion curves along the $\Gamma - X$ direction. Damping is included in the analysis using fractional derivative Zener model presented above.

Figures 10a to 10c show dispersion curves along the $\Gamma - X$ direction obtained with the shifted-cell operator method. A comparison between the results obtained using the reference structure and the adaptive structure is presented. The two states are obtained by changing the temperature of the polymeric interface: at 25°C, the resonant bandgap is visible around the selected frequency (40 kHz), as defined during the design of the resonator. Above the glass transition of the polymer, the phononic crystal slab tends to behave as an homogeneous damped plate.

The group velocity associated to dispersion curves along the $\Gamma - X$ direction at ambient temperature 25°C and 90°C respectively, using branch tracking shown in the subsection 2.3, is plotted figure 11a and 11b.

At this step, we have used dispersion analysis to design the temperature-driven adaptive phononic crystal slab, which is expected to follow the physical behaviour described by the diagrams. The purpose of the next section is to check the validity of the design when the crystal is embedded in a finite structure.
<figure><img src="image_11.png" /><figcaption>Figure 10. Dispersion curves along the Γ – X direction obtained with the shifted-cell operator method. a) Reference structure, all in aluminum. b) Metamaterial at ambient temperature 25°C. c) Metamaterial at 90°C and aluminum plate. Grey shapes represent bandgaps.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Group velocity associated to dispersion curves a) at 25°C, b) at 90°C. Grey shapes represent bandgaps.</figcaption></figure>

3.2.2 Embedding of the crystal in a finite structure : numerical modelling

The crystal (7 × 7 cells) is used as interface between two plates, as shown in figure 12. Cell size is 1 × 1 cm², the
total size of the structure is 21 × 7 cm².

<figure><img src="image_13.png" /><figcaption>Figure 12. Finite structure with an interface composed by 7 × 7 unit cells distributed. Piezoelectric patch (PZT) for excitation in the bottom left corner.</figcaption></figure>

A finite element model of the structure is built (with quadratic Lagrange elements). The structure is free.
A piezoelectric patch (transversal excitation) with harmonic voltage (|U| = 100 V) is included in the model
in order to be close to the experimental set up and covers all the frequency range from 0 to 50 kHz. Figure
13 shows numerical frequency responses obtained with the finite elements model for the two temperatures of interest. Squared velocity amplitudes $|Vz|^2$ are averaged on the domain of the input plate (IN) and the output plate (OUT). Blue and red curves are numerical frequency response at ambient temperature $25°C$ for the input plate (IN) and the output plate (OUT), respectively. Yellow and purple curves correspond to the same quantities at $90°C$. The grey shape represents the bandgap predicted by the dispersion analysis at ambient temperature.

<figure><img src="image_14.png" /><figcaption>Figure 13. Numerical frequency responses. Average squared velocity amplitude $|Vz|^2$ on the domain of the input plate (IN) and the output plate (OUT) at ambient temperature $25°C$ respectively in blue and red. At $90°C$, respectively in yellow and purple. The grey shape represents the bandgap predicted by the dispersion analysis at ambient temperature.</figcaption></figure>

At ambient temperature $25°C$, an output attenuation around 40 kHz is observed. Its width (about 3.1 kHz) is smaller than the frequency range predicted by the bandgap ($\Delta f = 8.9$ kHz). At $90°C$, after the glass transition, the output attenuation is not visible anymore.

In order to understand this unexpected bandgap width reduction, the vibration field is investigated at 4 frequencies of interest, illustrated in figure 14: a first point below the bandgap (34.4 kHz), two points in the expected bandgap (39.2 kHz, 43.3 kHz) and the last one above the bandgap (at 50.0 kHz). The corresponding operational deflections are shown in figures 15a, 15b, 15c and 15d.

As expected, energy can propagate through the lattice at 34.4 kHz (figure 15a) and 50.0 kHz (figure 15d), these two frequencies being outside of the predicted bandgap. At 39.2 kHz (figure 15b), inside the bandgap, the attenuation is clearly visible, all the energy is confined in the input plate (IN): the stop band effect is obtained.

However, at 43.3 kHz (figure 15c), the lattice is not efficient despite the fact that the frequency of interest is located inside the predicted bandgap. The operational deflection field clearly shows that in the central part of the crystal the vibration level is close to zero while the energy propagates along the edges of the periodic lattice: this phenomenon can not be predicted by the dispersion analysis. This phenomena is experimentally investigated in the next section.
<figure><img src="image_15.png" /><figcaption>Figure 14. Zoom on numerical frequency responses. Average squared velocity amplitude $|Vz|^2$ on the domain of the input plate (IN) and the output plate (OUT) at ambient temperature 25°C respectively in blue and red. The grey shape represents the bandgap predicted by the dispersion analysis at ambient temperature.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15. Several points are selected with corresponding deformed shapes a) at 34.4 kHz, b) 39.2 kHz, c) 43.3 kHz and d) 50.0 kHz.</figcaption></figure>

### 3.2.3 Experimental validation

Figure 16 shows the structure which has been tested. Polymer parts are realised by laser cutting and metallic parts by classical mechanical processing. The metamaterial is assembled in two steps by bonding parts together. The polymer cylinders are glued on the support aluminum plate. After drying and loading processes, at ambient temperature during 24 hours, aluminum cylinders are glued on the polymer cylinders. A piezoelectric patch (sintered ceramic PI, d31 effect, 8 × 7 × 0.2 mm³) provides the harmonic excitation up to 50 kHz.

Figure 17a shows the experimental facility composed by the metamaterial with its mounting bracket, a thermal chamber, a 3D vibrometer, a piezoelectric amplifier and a thermocouple.

The metamaterial is suspended to reproduce free boundary condition as shown in figure 17b. The mounting bracket is introduced in a thermal chamber which presents a glass wall allowing 3D vibrometer measurements.
<figure><img src="image_17.png" /><figcaption>Figure 16. Metamaterial after several manufacturing steps.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17. a) Experimental facility; b) focus on the test structure with free boundary conditions.</figcaption></figure>

The white noise generator provides random input voltage (5 V) between 500Hz and 50kHz, this signal is amplified 20 times by the piezoelectric amplifier. Temperature in the thermal chamber is controlled by a thermocouple. A 3D vibrometer is used to measure the velocity, and frequency responses are calculated with $H1$ estimator.⁴⁶

## Experimental results at 25°C

Figure 18a presents numerical and experimental frequency responses at ambient temperature 25°C. Curves need to be compared by pairs (yellow-blue and purple-red). Yellow and blue curves correspond, respectively, to experimental and numerical spatially averaged squared velocity/voltage amplitude (after ×20 amplification) $|Vz|^2/U^2$ for the input plate (IN). Purple and red curves correspond, respectively, to experimental and numerical spatially averaged squared velocity/voltage amplitude $|Vz|^2/U^2$ for the output plate (OUT). In the bandgap region, the attenuation predicted by numerical simulation is around 100dB, which corresponds to a ratio between output and input squared velocity amplitude of the order of $10^5$. Numerical amplitudes in this area are considerably lower than experimental amplitudes taking into account the experimental environment limitations. Reaching such a dynamic range in the measurement is not possible with the used experimental setup. However, the noise level estimation is performed and the green curve shows the background noise generated by the surrounding equipments (acquisition chain, room fan, etc.). It is hence impossible to measure vibration levels below this limit. The general trends (yellow-blue and purple-red) are similar except in the areas where the measurement is close to ambient noise.

Additionally to the velocity/voltage transfer functions presented in figure 18a, transmissibility functions
<figure><img src="image_19.png" /><figcaption>Figure 18. a) Numerical and experimental frequency responses at ambient temperature 25°C. Experimental spatially averaged squared velocity/voltage amplitude $|Vz|^2/U^2$ on the domain of the input plate (IN) and the output plate (OUT) respectively in blue and red. Numerical ones respectively in yellow and purple. The green curve shows the noise measurement. The grey shape represents the bandgap predicted by the dispersion analysis at ambient temperature. b) Numerical and experimental transmissibility functions IN/OUT at ambient temperature respectively in blue and red. Red areas show when velocity are in the same order of magnitude as noise.</figcaption></figure>

IN/OUT are presented in figure 18b. This transfer is defined, for both experimental and numerical data, as follows:

$$H = 20 \times \log_{10} \left( \frac{\Sigma^{out}(|Vz|^2/U^2)}{\Sigma^{in}(|Vz|^2/U^2)} \right), \qquad (23)$$

where $\Sigma^{out}(|Vz|^2/U^2)$ and $\Sigma^{in}(|Vz|^2/U^2)$ are, respectively, the average squared velocity amplitude $|Vz|^2$ on the domain of the input plate (IN) and the output plate (OUT).

Red areas show when velocity are in the same order of magnitude as noise. These results confirm that the transmissibility is well estimated by the model, however the bandgap is not clearly observed in the experiments.

Next, the experimental analysis is focused on the 4 frequency points pointed out in the numerical part. Figures 19 show experimental deformed shapes at 33.4 kHz, 38.2 kHz, 42.2 kHz and 50.0 kHz corresponding to numerical deformed shapes (Figure 14). The experimental mesh is not as fine as the finite element one, however the patterns can clearly be identified. Shapes are similar between numerical and experimental results. In particular, it can clearly be observed in figure 19c that the energy is transported along the borders at 42.2 kHz. This observation being experimentally validated, the section 4 aims will propose a new methodology to predict such configurations using cell-based computation.

## Experimental results at 90°C

In order to switch to the second state of the crystal, the structure is heated at 90°C. Figure 20a presents the corresponding numerical and experimental frequency responses. Curves need to be compared by pairs (yellow-blue and purple-red). Yellow and blue curves correspond, respectively, to experimental and numerical spatially averaged squared velocity/voltage amplitude $|Vz|^2/U^2$ on the domain of the input plate (IN). Purple and red curves correspond, respectively, to experimental and numerical spatially averaged squared velocity/voltage amplitude $|Vz|^2/U^2$ on the domain of the output plate (OUT). A good coherence is observed between the numerical and experimental results. Additionally to the velocity/voltage transfer functions presented in figure 20a, transmissibility functions IN/OUT are presented in figure 20b.
<figure><img src="image_20.png" /><figcaption>Figure 19. Experimental deformed shapes a) at 33.4 kHz, b) 38.2 kHz, c) 42.2 kHz and d) 50.0 kHz corresponding to numerical deformed shapes (Figures 14).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 20. a) Numerical and experimental frequency responses at ambient temperature 90°C. Experimental spatially averaged squared velocity/voltage amplitude $|Vz|^2/U^2$ on the domain of the input plate (IN) and the output plate (OUT) respectively in blue and red. Numerical ones respectively in yellow and purple. b) Numerical and experimental transmissibility functions IN/OUT at 90°C respectively in blue and red.</figcaption></figure>

The effect related to the periodicity is no longer observed at this temperature. The aluminum pillars are completely decoupled from the plate and the polymer still plays a damping role. Frequency response functions are smoothed at this temperature.

The prediction of the resonances at low frequencies is excellent. The discrepancies increase with frequency, the numerical model being stiffer than the real structure (no model updating has been performed). Resonance levels are quite well estimated. A poor damping model probably should explains the difference. Anyway, the physics is well captured by the model and the structure behaves as expected.

## 4. SUPERCELL ANALYSIS TO HANDLE FINITE BOUNDARY CONDITION

In this section, a supercell analysis is performed in order to predict the propagation of waves along the borders of the finite structure, without requiring the modelling of the full finite structure, which could induce high calculation costs. The principle consists in considering a set of seven connected cells in the analysis, as shown in figure 21. Boundary conditions are as follow:
* in the *x* direction, classical Floquet-Bloch conditions are used, namely $u_R = e^{-jk_{x}r}u_L$ where $u_R$ is the displacement on the right border and $u_L$ is the displacement on the left border in *y* axis,

* in the *y* direction, the left and right borders use the effective boundary conditions (free end in our case).

<figure><img src="image_22.png" /><figcaption>Figure 21. Boundary conditions of the 7-cell system.</figcaption></figure>

Figure 22a shows dispersion curves along the Γ – X path of the Brillouin zone obtained with the shifted-cell operator method at ambient temperature (25°C) for the 1-cell (classical Bloch analysis) and 7-cell systems. The blue curve corresponds to the 1 cell dispersion curve (figure 10b). The dispersion for the 7-cell system is plotted in yellow, two branches appear in the bandgap and the corresponding deformed shapes are visible in the figures 22b and 22c. Energy can propagate in the *x* direction along the free border. The interference effect leads to a non vibrating bulk and vibrating edges, which propagate vibration energy.

The proposed methodology hence allows identification of propagating waves along the borders inside the bandgap, cancelling the expected stop band effect. This method can highlight that the attenuation frequency range in finite structures is smaller than predicted by the bandgap (figure 13).
<figure><img src="image_23.png" /><figcaption>Figure 22. a) Dispersion curves along the Γ − X path of the Brillouin zone obtained with the shifted-cell operator method for the 1 cell and 7-cell system. b) Deformed shape at 40.1 kHz, c) Deformed shape at 40.2 kHz.</figcaption></figure>

5. CONCLUSION

In this article, numerical tools for the dispersion analysis in periodic damped structures are used, with a focus
on the ability of the shifted-cell operator method to deal with the dissipative behaviour of the system. An
adaptive metamaterial based on aluminum pillars glued on a dissipative polymer interface is designed. The
mechanical properties of the polymer depend on the frequency and the temperature, they change radically at
the glass transition. The viscoelastic behavior of tBA/PEGDMA is approximated by a fractional Zener model.
Two states are obtained by changing the temperature of the polymer interface to ambient temperature, the
bandgap predicted by the dispersion diagram is visible around the selected frequency. The attenuation in the
finite structure is not as important as that predicted by the bandgap in the infinite structure. Energy propagates
along each edge of the periodic lattice. This point is confirmed through experimental measurements, and the
dispersion analysis can not predict this effect. Finally, a supercell-based approach is used to handle finite system
boundary conditions in order to identify situations where energy transfer may arise because of reflections at the
border of the domain.
ACKNOWLEDGMENTS

This work was financed by The French National Research Agency under grant number ANR-12-JS09-008-COVIA.
It has been performed in cooperation with the Labex ACTION program (ANR-11-LABX-0001-01) and the hybrid
microfabrication platform MIFHySTO.

REFERENCES

[1] V. Veselago, “The electrodynamics of substances with simultaneously negative values of epsilon and mu,” *Soviet physics uspekhi* **10**(4), p. 509, 1968.

[2] J. Pendry, “Negative refraction makes a perfect lens,” *Physical review letters* **85**(18), p. 3966, 2000.

[3] V. Veselago and E. Narimanov, “The left hand of brightness: past, present and future of negative index materials,” *Nature materials* **5**(10), pp. 759–762, 2006.

[4] M. Kadic, T. Bückmann, R. Schittny, and M. Wegener, “Metamaterials beyond electromagnetism,” *Reports on Progress in Physics* **76**(12), p. 126501, 2013.

[5] J. Lee, J. Singer, and E. Thomas, “Micro-/nanostructured mechanical metamaterials,” *Advanced materials* **24**(36), pp. 4782–4810, 2012.

[6] M. Eidini and G. Paulino, “Unraveling metamaterial properties in zigzag-base folded sheets,” *Science advances* **1**(8), p. e1500224, 2015.

[7] M. Kadic, T. Bückmann, R. Schittny, and M. Wegener, “On anisotropic versions of three-dimensional pentamode metamaterials,” *New Journal of Physics* **15**(2), p. 023029, 2013.

[8] H. Huang and C. Sun, “Wave attenuation mechanism in an acoustic metamaterial with negative effective mass density,” *New Journal of Physics* **11**(1), p. 013003, 2009.

[9] M. Nouh, O. Aldraihem, and A. Baz, “Wave propagation in metamaterial plates with periodic local resonances,” *Journal of Sound and Vibration* **341**, pp. 53–73, 2015.

[10] J. Grima, R. Caruana-Gauci, M. Dudek, K. Wojciechowski, and R. Gatt, “Smart metamaterials with tunable auxetic and other properties,” *Smart Materials and Structures* **22**(8), p. 084016, 2013.

[11] P. Wang, F. Casadei, S. Kang, and K. Bertoldi, “Locally resonant band gaps in periodic beam lattices by tuning connectivity,” *Physical Review B* **91**(2), p. 020103, 2015.

[12] M. Oudich, M. Senesi, M. Assouar, M. Ruzenne, J. Sun, B. Vincent, Z. Hou, and T. Wu, “Experimental evidence of locally resonant sonic band gap in two-dimensional phononic stubbed plates,” *Physical Review B* **84**(16), p. 165136, 2011.

[13] M. Al-Lethawe, M. Addouche, A. Khelif, and S. Guenneau, “All-angle negative refraction for surface acoustic waves in pillar-based two-dimensional phononic structures,” *New Journal of Physics* **14**(12), p. 123030, 2012.

[14] L. Wang and K. Bertoldi, “Mechanically tunable phononic band gaps in three-dimensional periodic elastic-tomeric structures,” *International Journal of Solids and Structures* **49**(19-20), pp. 2881–2885, 2012.

[15] V. Candido de Sousa, C. Sugino, C. De Marqui Junior, and A. Erturk, “Adaptive locally resonant metamaterials leveraging shape memory alloys,” *Journal of Applied Physics* **124**(6), p. 064505, 2018.

[16] V. C. de Sousa, D. Tan, C. De Marqui Jr, and A. Erturk, “Tunable metamaterial beam with shape memory alloy resonators: Theory and experiment,” *Applied Physics Letters* **113**(14), p. 143502, 2018.

[17] K. Billon, I. Zampetakis, F. Scarpa, M. Ouisse, E. Sadoulet-Reboul, M. Collet, A. Perriman, and A. Hetherington, “Mechanics and band gaps in hierarchical auxetic rectangular perforated composite metamaterials,” *Composite Structures* **160**, pp. 1042–1050, 2017.

[18] L. DAlessandro, V. Zega, R. Ardito, and A. Corigliano, “3D auxetic single material periodic structure with ultra-wide tunable bandgap,” *Scientific reports* **8**(1), p. 2262, 2018.

[19] Y. Ma, F. Scarpa, D. Zhang, B. Zhu, L. Chen, and J. Hong, “A nonlinear auxetic structural vibration damper with metal rubber particles,” *Smart Materials and Structures* **22**(8), p. 084012, 2013.

[20] V. Laude, J.-C. Beugnot, S. Benchabane, Y. Pennec, B. Djafari-Rouhani, N. Papanikolaou, J. Escalante, and A. Martinez, “Simultaneous guidance of slow photons and slow acoustic phonons in silicon phoxonic crystal slabs,” *Optics express* **19**(10), pp. 9690–9698, 2011.

[21] W. Bragg, “Bakerian lecture: X-rays and crystal structure,” *Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character* **215**, pp. 253–274, 1915.
[22] M. Hussein, M. Leamy, and M. Ruzzene, “Dynamics of phononic materials and structures: Historical origins, recent progress, and future outlook,” *Applied Mechanics Reviews* **66**(4), p. 040802, 2014.

[23] S. Shan, S. Kang, Z. Zhao, L. Fang, and K. Bertoldi, “Design of planar isotropic negative Poissons ratio structures,” *Extreme Mechanics Letters* **4**, pp. 96–102, 2015.

[24] A. Slann, W. White, F. Scarpa, K. Boba, and I. Farrow, “Cellular plates with auxetic rectangular perforations,” *physica status solidi* (b) **252**(7), pp. 1533–1539, 2015.

[25] J. Joannopoulos, S. Johnson, J. Winn, and R. Meade, *Photonic crystals: molding the flow of light*, Princeton university press, 2011.

[26] T. Wu, T. Wu, and J. Hsu, “Waveguiding and frequency selection of Lamb waves in a plate with a periodic stubbed surface,” *Physical Review B* **79**(10), p. 104306, 2009.

[27] G. Floquet, “Sur les équations différentielles linéaires à coefficients périodiques,” in *Annales scientifiques de l'École normale supérieure, 12*, pp. 47–88, 1883.

[28] F. Bloch, “Über die quantenmechanik der elektronen in kristallgittern,” *Zeitschrift für physik* **52**(7-8), pp. 555–600, 1929.

[29] O. Schenk and K. Gärtner, “Solving unsymmetric sparse systems of linear equations with Pardiso,” *Future Generation Computer Systems* **20**(3), pp. 475–487, 2004.

[30] R. Craster, T. Antonakakis, M. Makwana, and S. Guenneau, “Dangers of using the edges of the Brillouin zone,” *Physical Review B* **86**(11), p. 115130, 2012.

[31] E. Manconi and B. Mace, “Estimation of the loss factor of viscoelastic laminated panels from finite element analysis,” *Journal of Sound and Vibration* **329**(19), pp. 3928–3939, 2010.

[32] B. Mace and E. Manconi, “Modelling wave propagation in two-dimensional structures using finite element analysis,” *Journal of Sound and Vibration* **318**(4), pp. 884–902, 2008.

[33] A. Krushynska, V. Kouznetsova, and M. Geers, “Visco-elastic effects on wave dispersion in three-phase acoustic metamaterials,” *Journal of the Mechanics and Physics of Solids* **96**, pp. 29–47, 2016.

[34] A. Palermo and A. Marzani, “Extended bloch mode synthesis: Ultrafast method for the computation of complex band structures in phononic media,” *International Journal of Solids and Structures* **100**, pp. 29–40, 2016.

[35] A. Bensoussan, J. Lion, and G. Papanicolaou, “Asymptotic analysis for periodic structure,” *North Holland*, 1978.

[36] C. Wilcox, “Theory of Bloch waves,” *Journal d'analyse mathématique* **33**(1), pp. 146–167, 1978.

[37] M. Collet, M. Ouisse, M. Ruzzene, and M. Ichchou, “Floquet–Bloch decomposition for the computation of dispersion of two-dimensional periodic, damped mechanical systems,” *International Journal of Solids and Structures* **48**(20), pp. 2837–2848, 2011.

[38] R. Moiseyenko and V. Laude, “Material loss influence on the complex band structure and group velocity in phononic crystals,” *Phys. Rev. B* **83**(6), p. 064301, 2011.

[39] V. Laude, J. Escalante, and A. Martínez, “Effect of loss on the dispersion relation of photonic and phononic crystals,” *Physical Review B* **88**(22), p. 224302, 2013.

[40] G. Vladimir and S. Marek, “Complex group velocity and energy transport in absorbing media,” *Phys. Rev. E* **81**(5), p. 056602, 2010.

[41] H. J. Morand and R. Ohayon, *Interactions fluides-structures*, Issy-les-Moulineaux (Hauts-de-Seine): Masson, 1992.

[42] P. Butaud, *Contribution à l'utilisation des polymères à mémoire de forme pour les structures à amortissement contrôlé*. PhD thesis, 2015.

[43] P. Butaud, M. Ouisse, V. Placet, and E. Foltête, “Experimental investigations on viscoelastic properties of a shape memory polymer,” in *ASME 2014 Conference on Smart Materials, Adaptive Structures and Intelligent Systems*, pp. V001T01A029–V001T01A029, American Society of Mechanical Engineers, 2014.

[44] P. Butaud, E. Foltête, and M. Ouisse, “Sandwich structures with tunable damping properties: On the use of shape memory polymer as viscoelastic core,” *Composite Structures* **153**, pp. 401–408, 2016.

[45] M. Williams, R. Landel, and J. Ferry, “The temperature dependence of relaxation mechanisms in amorphous polymers and other glass-forming liquids,” *Journal of the American Chemical society* **77**(14), pp. 3701–3707, 1955.
[46] G. E. Box, G. M. Jenkins, G. C. Reinsel, and G. M. Ljung, Time series analysis: forecasting and control, John Wiley & Sons, 2015.

## 6. ANNEX: GROUP VELOCITY

In this section, the details for the expression of the group velocity is given. Equation 2.3.2 gives the expression of the group slowness using $\lambda_i = jk_i$

$$ \frac{\partial k_i}{\partial \omega} = -j \frac{\psi_i^{lT} \left[ \frac{\partial \mathbf{A}_1(\omega)}{\partial \omega} - \lambda_i \frac{\partial \mathbf{A}_2(\omega)}{\partial \omega} \right] \psi_i^r}{\psi_i^{lT} \mathbf{A}_2(\omega) \psi_i^r}. \quad (24) $$

The matrices $\mathbf{A}_1(\omega)$ and $\mathbf{A}_2(\omega)$ are used

$$ \frac{\partial k_i}{\partial \omega} = -j \frac{\psi_i^{lT} \left( \begin{array}{ccc} 0 & & 0 \\ \frac{\partial(\mathbf{K}-\omega^2 \mathbf{M})}{\partial \omega} & -\lambda_i \frac{\partial \mathbf{L}^T}{\partial \omega} & \frac{\partial \mathbf{L}}{\partial \omega} - \lambda_i \frac{\partial \mathbf{H}}{\partial \omega} \end{array} \right) \psi_i^r}{\psi_i^{lT} \begin{pmatrix} \mathbf{I}_d & 0 \\ \mathbf{L}^T & \mathbf{H} \end{pmatrix} \psi_i^r}. \quad (25) $$

The expression left-eigenvector ($\psi_i^l$) and right-eigenvector ($\psi_i^r$) are replaced in the following expression

$$ \frac{\partial k_i}{\partial \omega} = -j \frac{\left( \lambda_i (\phi_i^{lT} \mathbf{H}) - \phi_i^{lT} \mathbf{L} \quad \phi_i^{lT} \right) \left( \frac{0}{\partial \omega} - \lambda_i \frac{\partial \mathbf{L}^T}{\partial \omega} \quad \frac{0}{\partial \omega} - \lambda_i \frac{\partial \mathbf{H}}{\partial \omega} \right) \begin{pmatrix} \phi_i^r \\ \lambda_i \phi_i^r \end{pmatrix}}{\left( \lambda_i (\phi_i^{lT} \mathbf{H}) - \phi_i^{lT} \mathbf{L} \quad \phi_i^{lT} \right) \begin{pmatrix} \mathbf{I}_d & 0 \\ \mathbf{L}^T & \mathbf{H} \end{pmatrix} \begin{pmatrix} \phi_i^r \\ \lambda_i \phi_i^r \end{pmatrix}}, \quad (26) $$

after the multiplication of the matrices, the expression of the slowness is found

$$ S_g = \frac{\partial k_i}{\partial \omega} = -j \frac{\phi_i^{lT} [-2\omega \mathbf{M} + \frac{\partial \mathbf{K}}{\partial \omega} - \lambda_i \frac{\partial \mathbf{L}^T}{\partial \omega} + \lambda_i \frac{\partial \mathbf{L}}{\partial \omega} - \lambda_i^2 \frac{\partial \mathbf{H}}{\partial \omega}] \phi_i^r}{\phi_i^{lT} [-\mathbf{L} + \mathbf{L}^T + 2\lambda_i \mathbf{H}] \phi_i^r}. \quad (27) $$

The group velocity is the inverse of the group slowness

$$ C_g = \frac{1}{\frac{\partial k_i}{\partial \omega}} = \frac{1}{-j \frac{\phi_i^{lT} [-2\omega \mathbf{M} + \frac{\partial \mathbf{K}}{\partial \omega} - \lambda_i \frac{\partial \mathbf{L}^T}{\partial \omega} + \lambda_i \frac{\partial \mathbf{L}}{\partial \omega} - \lambda_i^2 \frac{\partial \mathbf{H}}{\partial \omega}] \phi_i^r}{\phi_i^{lT} [-\mathbf{L} + \mathbf{L}^T + 2\lambda_i \mathbf{H}] \phi_i^r}}. \quad (28) $$

For frequency-independent cases, a simplification is possible

$$ C_g = \frac{1}{\frac{-2j\omega\phi_i^{lT}\mathbf{M}\phi_i^r}{[\phi_i^{lT}(\mathbf{L}-\mathbf{L}^T)\phi_i^r - 2\lambda_i\phi_i^{lT}\mathbf{H}\phi_i^r]}}. \quad (29) $$

In the particular case of on homogeneous cases where the frequency dependency is characterized by a Young's modulus such that $E = f(\omega)E_0$ and a constant Poisson's ratio. Hence $\mathbf{K} = f(\omega)\mathbf{K}_0$, $\mathbf{H} = f(\omega)\mathbf{H}_0$ and $\mathbf{L} = f(\omega)\mathbf{L}_0$. The equivalent quadratic eigenvalue problem is

$$ [f(\omega)(\mathbf{K}_0 + \lambda_i(\mathbf{L}_0 - \mathbf{L}_0^T) - \lambda_i^2 \mathbf{H}_0) - \omega^2 \mathbf{M}] \phi_i^r = 0 \quad (30) $$

In this case, the slowness is expressed as
$$S_g = \frac{\partial k_i}{\partial \omega} = -j \frac{\phi_i^{lT} [-2\omega \mathbf{M} + \frac{\partial f}{\partial \omega} (\mathbf{K}_0 + \lambda_i (\mathbf{L}_0 - \mathbf{L}_0^T) - \lambda_i^2 \mathbf{H}_0)] \phi_i^r}{\phi_i^{lT} [f(\omega)(-\mathbf{L}_0 + \mathbf{L}_0^T + 2\lambda_i \mathbf{H}_0)] \phi_i^r}. \quad (31)$$

Using the expression of the quadratic eigenvalue problem

$$\left[ \frac{\partial f}{\partial \omega} (\mathbf{K}_0 + \lambda_i (\mathbf{L}_0 - \mathbf{L}_0^T) - \lambda_i^2 \mathbf{H}_0) \right] \phi_i^r = \left[ \frac{\partial f}{\partial \omega} \frac{\omega^2 \mathbf{M}}{f(\omega)} \right] \phi_i^r, \quad (32)$$

so, the slowness can be written as

$$S_g = \frac{\partial k_i}{\partial \omega} = \frac{j \phi_i^{lT} [\omega^2 \frac{\partial f}{\partial \omega} - 2\omega] \mathbf{M} \phi_i^r}{\phi_i^{lT} [f(\omega)(\mathbf{L}_0 - \mathbf{L}_0^T - 2\lambda_i \mathbf{H}_0)] \phi_i^r}, \quad (33)$$

Therefore, the group velocity is expressed as

$$C_g = \frac{1}{\frac{\partial k_i}{\partial \omega}} = \frac{1}{\frac{j \phi_i^{lT} [\omega^2 \frac{\partial f}{\partial \omega} - 2\omega] \mathbf{M} \phi_i^r}{\phi_i^{lT} [f(\omega)(\mathbf{L}_0 - \mathbf{L}_0^T - 2\lambda_i \mathbf{H}_0)] \phi_i^r}} \quad (34)$$