On the Free Vibrations of Non-Classically Damped Locally Resonant Metamaterial Plates

In this paper, the focus is on the free vibrations of locally resonant metamaterial plates with viscously damped resonators. Upon formulating a dynamic-stiffness model where the resonators are represented via pertinent reaction forces depending on the deflections of the attachment points, the complex eigenvalues are calculated by a contour-integral algorithm introduced in the literature for general nonlinear eigenvalue problems. The interest in the proposed approach is twofold. The dynamic-stiffness model involves a limited number of generalised coordinates compared to the nodal degrees of freedom of a standard finite-element model, and the contour-integral algorithm proves successful in evaluating all complex eigenvalues, without missing any one, with remarkable computational efficiency. Numerical results are presented for Lévy plates, but are readily extendible to other plate theories. Finally, an ad hoc dynamic-stiffness approach is formulated to calculate the frequency response of the plate under arbitrarily placed loads, which is of particular interest to investigate its elastic wave attenuation properties.


Introduction
There is a considerable body of recent literature on locally resonant metamaterial plates (LRMPs), i.e., plates engineered with a periodic array of small resonators. Indeed, these plates exhibit the inherent attenuation properties of elastic waves that make them ideally suitable for several applications in dynamics. Typically, the resonators may be distributed over the external surface of homogeneous or composite plates [1][2][3][4][5][6][7][8][9][10] or, alternatively, they may be embedded within the core matrix of sandwich plates [11][12][13].
Computational models of LRMPs generally rely on standard finite-element (FE) analysis. The model involves a standard FE for the plate, which is coupled with concentrated mass-spring or mass-spring-damper subsystems representing the resonators [6,8]. In some cases, the resonators are modelled by FEs as well [12,13].
A very appealing approach for modelling bare plates, i.e., plates without resonators, is a dynamic-stiffness approach, as it ensures a very accurate description of the plate dynamics with a limited number of generalised coordinates compared to a standard FE model. For instance, in a comprehensive treatment for plates of various types, Banerjee and coworkers [14][15][16][17][18][19][20][21] developed dynamic-stiffness models where a few coefficients of appropriate Fourier series expansions proved capable of representing, with remarkable accuracy, the plate dynamic response. In this context, natural frequencies of the undamped modes were calculated by the powerful Wittrick-Williams algorithm [17,22], without missing any one and including multiple ones. As for LRMPs, a pertinent dynamic-stiffness model and an ad hoc formulation of the Wittrick-Williams algorithm were recently proposed by Russillo et al. [23].
w(x, y, t) being the transverse deflection, h the thickness, ρ the volumetric mass density, and D the bending rigidity of the plate. A solution that enforces simply supported BCs at y = 0 and y = L is considered: w(x, y, t) = ∞ ∑ m=1 W m (x) sin(α m y)e iωt (2) where α m = mπ L for m = 1, 2, . . . , ∞. Replacing Equation (2) in Equation (1) yields: Introduce the vectors collecting the generalised displacements and forces along the two unconstrained edges of the strip: where Φ ym (x) = −dW m /dx is the rotation in the xz plane and V xm (x) and M xm (x) are the shear force and bending moment per unit length along the edge of the strip. Making use of the Kirchhoff plate equations providing V xm (x) and M xm as functions of W m (x), the following matrix relation is obtained [14]: where D s m is the dynamic-stiffness matrix of the single strip for m = 1, 2, . . . , ∞. The dynamic-stiffness matrix of a plate consisting of n e − 1 strips (n e is number of lines) is assembled in a finite-element fashion, and the following equation of equilibrium holds: Truncating the Fourier series (2) up to N terms, Equation (6) can be equivalently written in matrix form:  It is noticed that the displacement/force q i (y) along the ith line, either at the boundary or between two adjacent strips of the assembled plate, can be expanded into a sine Fourier series: where the coefficients Q im are given as: Now, consider a resonator placed on the ith line at y = y j . Resonators attached along the ith line exert a force per unit length: where δ is Dirac's delta function, n s is the number of resonators attached along the ith line, and k eq (ω) is the frequency-dependent stiffness of the resonator, which can be obtained from its dynamic-stiffness matrix, as demonstrated by the authors in [23,[29][30][31]. By means of Equation (8), the displacement of the attachment point is: Replacing Equation (11) in Equation (10) yields: Applying the transformation in Equation (9) to Equation (12) gives the Fourier series coefficients associated with the reaction forces exerted by the resonators: Equation (13) is readily written in matrix form as: matrix D res being defined as follows: Matrix D res in Equation (15) is the dynamic-stiffness matrix associated with the resonators. Now, consider Equation (15), and define the following diagonal matrix with n e submatrices: Notice that the ith block refers to the ith line of the assembled plate. Using Equation (16) and truncating the Fourier series (13) up to N terms, the following relation is readily obtained, which mirrors Equation (7): Finally, replacing Equation (17) for f in Equation (7) leads to the following nonlinear eigenproblem for the free-vibration response: The matrix D(ω) is the dynamic-stiffness matrix of the Lévy LRMP.

Contour-Integral Algorithm
In order to compute the eigenvalues of the Lévy LRMP, it is necessary to solve the nonlinear eigenproblem given by Equation (18). The presence of viscous dampers within the resonators makes the system non-classically damped and, consequently, the eigenvalues are complex. In this case, the usual algorithm employed to solve the eigenproblems associated with dynamic-stiffness models, that is the Wittrick-Williams algorithm, is no longer applicable. Therefore, the recently introduced contour-integral algorithm formulated by Asakura and coworkers [26] is used here, for the first time, to solve this challenging problem. The application of the algorithm follows these steps:
Compute two complex random source matrices U and V having dimensions n 0 × L 0 , n 0 being the size of the dynamic-stiffness matrix D(ω) in Equation (18) and L 0 the number of source vectors collected in U and V; 3.
Compute the shifted and scaled moments M k using the N 0 -point trapezoidal rule: with K the maximum moment degree considered for the moment and U H the Hermitian transpose of U; 4.
Omit small singular-value components σ i < · max i σ i ; set the numberm of remaining singular value components (m < KL 0 ); constructĤm andĤ < m extracting the principal submatrix with maximum indexm fromĤ KL 0 andĤ < KL 0 , i.e., Compute the eigenvalues ζ j of the linear pencil: Calculate the eigenvalues: The algorithm allows calculating all the eigenvalues ω j of the nonlinear eigenproblem (18) that fall within the selected circle Γ, including multiple roots [25][26][27].

Frequency Response
The frequency response of a Lévy LRMP can be computed by considering the equilibrium equation: A distributed force p i (y) and moment m i (y) acting along the ith line, either at the boundary or between two adjacent strips of the assembled plate, can be expanded into a sine Fourier series by Equation (8), whose coefficients V im and M im are: The coefficients in Equation (20) are collected in the subvector the force vector f m in Equation (6) (represented in Figure 2), which in turn is collected in the global force vector f in Equation (19). Once the vector f is built, the Fourier coefficients of the frequency response are given as: where u = u T The displacement u i (y) and rotation φ i (y) along the ith line are simply given as: where W im and Φ im are respectively the coefficients of the Fourier series expansion of the displacement and rotation along the ith line. Notice that concentrated loads can readily be handled within the framework above: for this, two strips separated by a line passing through the load application point can be considered, and the load can be modelled as a standard 1D Dirac's delta.

Numerical Applications
To validate the proposed method, consider a simply supported locally resonant steel plate coupled with 2-DOF viscously damped resonators and having dimensions 0.20 m × 0.20 m × 0.01 m. The plate and resonator parameters are: Young's modulus E = 200 GPa, Poisson's ratio ν = 0.3, volumetric mass density ρ = 7750 kg/m 3 The proposed dynamic-stiffness approach is implemented by applying the contourintegral algorithm to solve the eigenvalue problem involving the dynamic-stiffness matrix in Equation (18), where the frequency-dependent stiffness of the resonators is calculated as in [30]. The parameters of the contour-integral algorithm are N 0 = 36, L 0 = 60, K = 15.
The first 150 eigenvalues are computed using the proposed method and, in order to assess the accuracy, modelling the plate in Figure 3 with the FE code ABAQUS using a mesh of 200 × 200 S4R5 elements. Some eigenvalues are reported in Table 1, while the complete list of the first 150 is given in Tables A1-A3 in Appendix A, for conciseness. The calculated eigenvalues are in excellent agreement and the maximum relative error computed is equal to = 3.65% for the real part and = 0.78% for the imaginary part. Furthermore, it is seen that the contour-integral algorithm is capable of capturing eigenvalues very close to each other, which may differ even by a few digits. An excellent agreement is found in terms of mode shapes, as shown in Figures 4-8. Notice that the mode shapes are complex, as they are associated with complex eigenvalues; see Tables 1 and A1-A3 in Appendix A; in particular, Figures 4-8 show the dimensionless real parts, as obtained upon dividing all components by the one with the maximum absolute value.
As for the size of the two models, the proposed dynamic-stiffness model involves a 144 × 144 dynamic-stiffness matrix (corresponding to eight strips with nine generalised coordinates each), while the FE model involves 11,8791 × 11,8791 stiffness and mass matrices (corresponding to 40,401 nodes). It is seen that the eigenvalues and the maximum relative error do not change appreciably by increasing the sizes of the two models. Moreover, no significant differences are found in the eigenvalues if the parameters of the contourintegral algorithm are changed. As for the computational effort, it is noteworthy that the contour-integral algorithm is implemented in an in-house MATLAB code. For the parameters assumed, i.e., N 0 = 36, L 0 = 60, K = 15, the computation time required to calculate the eigenvalues in Tables 1 and A1-A3 of Appendix A is 44.4 s. On the other hand, the corresponding eigenvalues of the finite-element model in ABAQUS are calculated by the default Lanczos eigensolver, the computation time of which is 111.2 s. Therefore, it can be concluded that the dynamic-stiffness model in conjunction with the contour-integral algorithm provides very accurate results compared to a standard FE approach, requiring a very limited size of the model and computational effort. Table 1. Complex eigenvalues (from 1 to 150) of the LRMP in Figure 3 with 2-DOF resonators, computed by means of the proposed dynamic-stiffness approach (DS) and the FE method (FE).

Mode
Eigenvalue (   Finally, the dynamic-stiffness approach proposed in Section 4 is applied to calculate the transmittance of the plate in Figure 3, considering two sets of BCs: (i) all edges are simply supported; (ii) two edges are simply supported and two are free. Specifically, a unit harmonic concentrated load is applied at (x 0 , y 0 ) = (0.025, 0.10), and the transmittance is evaluated as the ratio of the deflection at (x 1 , y 1 ) = (0.175, 0.10) to the deflection at the load application point (x 0 , y 0 ). For completeness, Figure 9 includes the band gaps of the corresponding infinite LRMP, as computed by the standard FE approach [32]. The transmittance varies significantly with the BCs and, within the band gaps of the corresponding infinite plate, better wave attenuation properties of the finite plate are found when two edges of the plate are free and two simply supported, compared to the case where all the edges are simply supported.

Discussion
This paper has addressed the free vibrations of non-classically damped locally resonant metamaterial plates, focusing on Lévy plates. The main novelties of the proposed approach can be summarized as: (1) the formulation of a reduced-order dynamic-stiffness model where multi-degree-of-freedom resonators are represented via frequency-dependent reaction forces, involving the deflection of the attachment points only and no resonator DOFs; (2) the application of the contour-integral algorithm to calculate the complex eigenvalues; (3) a dynamic-stiffness approach to the calculation of the frequency response under arbitrarily placed concentrated loads, which is of interest to investigate the elastic wave attenuation properties of the plate. Considering a simply supported plate coupled with 2-DOF viscously damped resonators, it has been demonstrated that the proposed approach provides very accurate eigenvalues and mode shapes, with a very limited size of the model and computational costs compared to a standard FE approach implemented in ABAQUS. It is noteworthy that the proposed approach requires only that the dynamic-stiffness matrix of the bare plate is available. As such, it is readily extendible to locally resonant metamaterial Kirchhoff plates with arbitrary BCs, using the appropriate dynamic-stiffness matrix formulated in [18,19,21].
It is important to remark that the LRMPs under study are of interest not only as structural/mechanical components at the macroscale, but also at much smaller scales. For instance, the concept of periodic arrays of resonators coupled with a primary host system was proposed at the microscale, e.g.,: Reference [33] investigated numerically and experimentally the formation of locally resonant band gaps in a 2D surface phononic crystal with inverted conical pillars; specifically, the inverted conical pillars were deposited on a semi-infinite lithium-niobate substrate and arranged in a honeycomb lattice array for applications in low-frequency guiding, acoustic wave isolation, acoustic absorbers, and acoustic filters. Moreover, numerical and experimental investigations demonstrated the existence of band gaps at a multi-GHz frequency range in a pillar-based hypersonic 2D phononic crystal with nanoscale dimensions [34]; in this case, the fabricated phononic crystal consisted of a periodic array of nanopillars arranged according to the triangular lattice structure of the crystal. An interesting discussion on the potential applications of locally resonant phononic nanostructures with other pertinent references may be found in the work by Guo et al. [35]. These recent studies demonstrate the interest in LRMPs at various scales, for which the proposed approach may represent a valuable analysis tool.

Data Availability Statement:
The raw/processed data required to reproduce the above findings cannot be shared at this time as the data also forms part of an ongoing study.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

BC
Boundary condition DOF Degree of freedom FE Finite element LRMP Locally resonant metamaterial plate Appendix A. Eigenvalues Table A1. Complex eigenvalues (from 1 to 50) of the LRMP in Figure 3 with 2-DOF resonators, computed by means of the proposed dynamic-stiffness approach (DS) and the FE method (FE).

Mode
Eigenvalue ( Table A3. Complex eigenvalues (from 101 to 150) of the LRMP in Figure 3 with 2-DOF resonators, computed by means of the proposed dynamic-stiffness approach (DS) and the FE method (FE).