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

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 aluminium deposited periodically on the aluminium 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 does not 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.


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 porpagations. 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 behaviour was developed in the 1960s by physicist Viktor Veselago [1]. It was not until the 2000s that the first experimental realisation [2,3] was achieved. Mechanical metamaterials have been recently hailed as a new class of structural concept able to bring novel multifunctionalities [4] by changes of compliance, shape, or by embedding oscillators or smart material inserts. Some examples (to name a few) are multiscale architecturally structured topologies [5], zig-zag folded sheets [6], pentamodal lattices [7], systems with distributed resonators [8,9], smart/magnetic materials [10], tunable connectivity [11], phononic stubbed plates [12], flat lenses (super lenses) [13], 3D tunable phononic crystals [14], tunable metamaterial beam with shape memory alloy resonators [15,16], auxetic periodic structure (with negative Poisson ratio) [17,18] and nonlinear auxetic dampers [19].
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 [20].
The control of the elastic waves can be performed by combining Bragg's bandgap [21] (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 realisation 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 [22]. 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 [17,23,24]. The methods currently used are most of the time based on those derived from wave propagation in crystals [25], 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 behaviour 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 supercellbased 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 twocell with adequate boundary conditions. The methodology is described and validated using the full finite model and experimental tests on a 2D metamaterial structure. The reference structure used in this paper is presented in [26]. The system consists of an infinite 3D periodic bidirectional

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 PDE r w is the strain tensor, ρ is the material density and ω is the circular frequency.
The periodicity is defined on the borders of the domain using the Floquet-Bloch boundary conditions, that write where u R (resp. v R ) is the displacement on the right border and u L (resp. v L ) is the displacement on the left border in x (resp. y) axis, k x and k y are respectively the wavenumbers in the x and y directions [27,28] as shown in figure 1(b). A parametric eigenvalue analysis is performed using the Pardiso solver [29], two parameters (wavenumbers) are considered, namely k x =[0 π/r 1 ] and k y =[0 π/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 2(a)) and on the contours ( figure 2(b)). The frequencies defining the bandgaps can be found by considering only the contour of the irreducible Brillouin zone for nonsingular 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 2(b) 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 [33], where the damping provided by a generalised Maxwell model is included in the stiffness matrix, and [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.

Shifted-cell operator method
The shifted-cell operator, alternative formulation proposed by [35][36][37], 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 [37] 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: where l = jk i i is the ith eigenvalue, f i r is the right eigenvector associated to l M , i and K are respectively the standard symmetric definite mass and symmetric semi-definite stiffness matrices, -L L T is a skew-symmetric matrix and H is a symmetric semi-definite positive matrix for details). In this formulation, all matrices can depend on ω. A parametric eigenvalue analysis is then performed where the pulsation ω and the wave propagation angle θ (such that q = k k cos x and q = k k sin y ) are fixed as real parameter, allowing introduction of damping effects. The wavenumbers l = jk i i and the associated right eigenvectors f 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 Conversely, a left-eigenvector for the same eigenvalue satisfies Equation (8) can be developed to find the link between y i r and y i l . After a few steps the following relationship is obtained : Hence, by solving the right eigenvalue problem, the ith mode ( So, by solving the right eigenvalue problem, the adjoint solution is found too.

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 [38,39].
Equation (4) is differentiated and multiplied by the left eigenvector y i lT following the procedure proposed by [38] such that y y y y y y y y y y According to equation (8) After several steps, isolating the l w ¶ ¶ i expression, the following formulation is obtained: In this case, the group slowness (see appendix for details) is expressed as The complex group velocity is the inverse of the group slowness 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.

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 wave- • 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 = • Ratio between real and imaginary part of the group velocity Only the waves corresponding to C 1 >τ 1 ; C 2 >τ 2 et C 3 >τ 3 are considered propagative. In practice, the thresholds τ 1 , τ 2 et τ 3 are chosen such as τ 1 =τ 2 =1; τ 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.  1(a)) is used. Figure 3 presents dispersion curves along the Γ−X (f=0°) 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 to 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 shiftedcell operator is first performed on a conservative structure to validate the implementation of the shifted-cell operator technique. Figure 4(a) presents the dispersion curves obtained along the Γ−X direction when f is equal to zero. Shiftedcell 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 section 2.3.5.  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 [37]. 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 [41]. . This point is defined as next point and step by step the branch is identified. The sorting algorithm has been validated with damping. 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 colour 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.

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 aluminium 6063-T83.
As it will be shown in the next section, the polymer thickness has been chosen in order to open a resonant bandgap around 40 kHz, which means below the Bragg bandgap frequency. The polymer has a density r =    In the harmonic regime, the complex modulus of the polymer writes 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  (table 1).
The translation factor a T comes from the Williams-Landel-Ferry law [45], the temperature dependency of a T is described by ). The relaxation time τ is linked to the translation factor : The storage modulus and the loss factor are plotted in figures 9(a) and (b) for two temperatures of interest, namely 25°C and 90°C. At ambient temperature, the polymer is stiff and its loss factor is almost equal to zero. At 90°C, the loss factor is over 1.5 on the frequency range (0-100 kHz), with a maximum value of 2.5 and the stiffness is very low.

Dynamical properties
3.2.1. Dispersion analysis. The shifted-cell operator technique is used to obtain the dispersion curves along the Γ−X direction. Damping is included in the analysis using fractional derivative Zener model presented above. Figures 10(a)-(c) show dispersion curves along the Γ−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 Γ−X direction at ambient temperature 25°C and 90°C respectively, using branch tracking shown in the section 2.3, is plotted figures 11(a) and (b).
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.

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 2 , the total size of the structure is 21×7 cm 2 .
A finite element model of the structure is build (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.
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 (Δ 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 As expected, energy can propagate through the lattice at 34.4 kHz ( figure 15(a)) and 50.0 kHz ( figure 15(d)), these two frequencies being outside of the predicted bandgap. At      15(c)), 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 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 aluminium plate. After drying and loading processes, at ambient temperature during 24 h, aluminium cylinders are glued on the polymer cylinders. A piezoelectric patch (sintered ceramic PI, d31 effect, 8×7×0.2 mm 3 ) provides the harmonic excitation up to 50 kHz. Figure 17(a) shows the experimental facility composed by the metamaterial with its mounting bracket, a thermal chamber, a 3D vibrometer, a piezoelectric amplifier and a thermocouple.

Experimental validation.
The metamaterial is suspended to reproduce free boundary condition as shown in figure 17(b). The mounting bracket is introduced in a thermal chamber which presents a glass wall allowing 3D vibrometer measurements.
The white noise generator provides random input voltage (5 V) between 500 Hz and 50 kHz, 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 [46].  Figure 18(a) 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 U 2 2 for the input plate (IN). Purple and red curves correspond, respectively, to experimental and numerical spatially averaged squared velocity/voltage amplitude | | Vz U 2 2 for the output plate (OUT). In the bandgap region, the attenuation predicted by numerical simulation is around 100 dB, 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 18(a), transmissibility functions IN/OUT are presented in figure 18(b). This transfer is defined, for both experimental and numerical data, as follows: 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. Figure 19 show experimental deformed shapes at 33.4, 38.2, 42.2 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 19(c) 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.  The effect related to the periodicity is no longer observed at this temperature. The aluminium 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.

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 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 22(a) 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 10(b)). The dispersion for the 7-cell system is plotted in yellow, two branches appear in the bandgap and the

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 aluminium 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 behaviour 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.

Appendix. Group velocity
In this section, the details for the expression of the group velocity is given. Equation The matrices w ( ) A 1 and w ( )