# Dynamic analysis of periodic structures and metamaterials via wave approaches and finite element procedures

Jean-Mathieu Mencik, Denis Duhamel

► To cite this version:

Jean-Mathieu Mencik, Denis Duhamel. Dynamic analysis of periodic structures and metamaterials via wave approaches and finite element procedures. 8th International Conference on Computational Methods in Structural Dynamics and Earthquake Engineering, Jun 2021, Athens, Greece. hal-03368952

HAL Id: hal-03368952

https://hal.science/hal-03368952

Submitted on 7 Oct 2021

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
DYNAMIC ANALYSIS OF PERIODIC STRUCTURES AND
METAMATERIALS VIA WAVE APPROACHES AND FINITE ELEMENT
PROCEDURES

J.-M. Mencik¹ and D. Duhamel²

¹INSA Centre Val de Loire,
Université d'Orléans, Université de Tours,
Laboratoire de Mécanique Gabriel Lamé, EA 7494,
Rue de la Chocolaterie, 41000 Blois, France
e-mail: jean-mathieu.mencik@insa-cvl.fr

²Ecole des Ponts, laboratoire Navier,
6 et 8 Avenue Blaise Pascal, Cité Descartes, Champs sur Marne,
77455 Marne la Vallée, cedex 2, France
e-mail: denis.duhamel@enpc.fr

**Keywords:** Periodic Structures, Metamaterials, Wave Approaches, Finite Elements, Model Reduction, Nearly Periodic Structures.

**Abstract.** Two kinds of numerical approaches able to predict the dynamic response of periodic structures and metamaterials are presented. Both of them are model reduction techniques which can be used to obtain the response functions of the structures at a low computational cost. The first kind of approaches employs the wave finite element (WFE) method for modeling 1D periodic structures made up of complex substructures, e.g., 2D cells of arbitrary shapes. Resonant metamaterials are particular cases of periodic structures where the substructures possess local resonances — e.g., layered substructures involving soft layers and heavy layers/core —, which in turn lead to band gap effects and low vibration levels. Although well-suited for modeling periodic structures in the frequency domain, the WFE method also works for periodic structures subject to time-dependent excitations. The second kind of approaches employs finite element (FE) procedures for modeling 2D periodic structures, and 2D nearly periodic structures involving substructures with slightly varying (random) geometrical properties. For nearly periodic structures like plates with disordered resonant 2D cells, results show that the vibrational energy is localized around the excitation sources.
# 1 INTRODUCTION

Periodic structures are frequently encountered in engineering applications, e.g., in the aero-nautic (fuselages, turbines) and railway industries. Also, they appear interesting to passively control the vibration levels of mechanical systems. Metamaterials are periodic structures which possess such interesting features, i.e., the fact that they do not convey vibrational energies on some frequency bands (band gaps) which as such yield low vibration levels. Well-known exam-ples of metamaterials are periodic structures made up of resonant cells where band gaps occur in the vicinity of the resonance frequencies of the cells. There exist plenty of works in the liter-ature about the analysis of the wave propagation in 1D or 2D periodic structures and the related band gap effects [1, 9, 8, 29, 25, 4, 28]. However, the forced response analysis of periodic structures is much less reported. This especially means developing numerical models able to predict the dynamic response of periodic structures of finite dimensions with various kinds of boundary conditions, or assemblies of finite dimensions involving several periodic structures and other non-periodic components which are systems of practical interest in engineering ap-plications. In this case, the analysis of the band gap effects does not only rely on the wave propagation properties, but also on the boundary and coupling conditions which induce energy conversion between waves.

This paper aims at investigating the modeling of 1D and 2D periodic structures, of finite dimensions, subject to various kinds of excitations and boundary conditions. Periodic structures made up of complex cells — i.e., substructures which can be of various shapes and which are modeled with a moderately high number of degrees of freedom (DOFs) — are especially dealt with. Numerical results are brought which concern the dynamic response of periodic structures and resonant metametarials.

The paper starts with an overview of the wave finite element (WFE) method and its recent applications. The method is well suited for modeling 1D periodic structures, i.e., structures made up of identical substructures which are connected to each other along a straight or cir-cumferential direction. The WFE method works by considering the FE model of a substructure, and computing the eigenvalues and eigenvectors of the related transfer matrix to determine the waves traveling along a periodic structure. To predict the forced response of periodic structures, the WFE method involves expanding the displacement and force vectors on the basis of the eigenvectors (wave shapes) of the transfer matrix of a substructure. The methodology is de-tailed in [16, 18] (see also [30, 6, 32, 15]). Among the recent extensions of the method are: the development of model reduction strategies based on reduced wave bases to target substructures with large-sized FE models [22]; the modeling of periodic structures with cyclic symmetry, such as turbines [19]; the time response analysis of periodic structures via wave-based absorb-ing boundary conditions [7].

Aside from the study of 1D periodic structures, 2D periodic structures with resonant sub-structures are well-known examples of metamaterials to passively control the vibration levels of mechanical systems. As mentioned earlier, the analysis of the wave propagation (e.g., Bloch waves) in 2D periodic structures is well documented in the literature. However, contrary to the study of 1D periodic structures, the forced response analysis of 2D periodic structures of finite dimensions, based on wave expansion (displacement and force vectors), is prone to tough theoretical and numerical modeling issues which make it difficult even impossible to carry out. The issue is linked to the multi-directional behavior of the propagating waves, and the descrip-tion of the wave reflection at boundaries that arise from these multi-directional waves. Here, two FE-based approaches are proposed as alternative model reduction techniques. The first
one is standard and makes use of the Craig-Bampton (CB) method and the related fixed interface modes of the substructures for modeling purely periodic structures. The second approach targets the modeling of 2D nearly periodic structures made up of resonant substructures with slightly varying geometrical properties [20]. Plates with disordered resonant 2D cells, made up of soft layers and heavy layers/core with random geometrical properties, are examples of 2D nearly periodic structures. Such structures possess interesting features, e.g., the fact that the vibrational energy remains localized around the excitation points. The proposed strategy involves computing the reduced matrices of the substructures for some particular distorted FE meshes (a few number), and interpolating these matrices between these “interpolation points” for modeling substructures with random FE meshes.

## 2 WFE METHOD FOR 1D PERIODIC STRUCTURES

### 2.1 Wave propagation

The WFE method originates from the transfer matrix theory and the Bloch’s theorem. Let us first consider an infinite periodic structure under harmonic disturbance $e^{i\omega t}$ made up of identical substructures with similar FE meshes as shown in Fig. 1. Also, let us denote by M, C and K the conventional mass, damping and stiffness matrices of the substructures (identical for all the substructures). From the practical point of view, these matrices can be obtained via usual FE softwares, Matlab FE mesh generators (e.g., DistMesh [27]) or Matlab in-house FE codes. The related dynamic equilibrium equation, for a given substructure k, is given by:

$$Dq^k = F^k, \qquad (1)$$

where $q^k$ and $F^k$ refer to the displacement vector and the force vector (respectively) of the substructure k; D is the dynamic stiffness matrix expressed by $D = -\omega^2 M + i\omega C + K$.

<figure><img src="image_2.png" /><figcaption>Figure 1: Periodic structure of infinite length, and FE mesh of a substructure.</figcaption></figure>

The FE mesh of the substructure involves left (L) and right (R) boundaries (see Fig. 1) which are here assumed to be described with a same number n of DOFs. From Eq. (1), the following transfer matrix relation between the right and left boundaries can be proposed:

$$\begin{bmatrix} q_R^k \\ F_R^k \end{bmatrix} = S \begin{bmatrix} q_L^k \\ -F_L^k \end{bmatrix}, \qquad (2)$$
where S is the so-called transfer matrix (size of $2n \times 2n$), expressed by:

$$S = \left[ \frac{-(D_{LR}^*)^{-1}D_{LL}^*}{D_{RL}^* - D_{RR}^*(D_{LR}^*)^{-1}D_{LL}^*} \middle| \frac{-(D_{LR}^*)^{-1}}{-D_{RR}^*(D_{LR}^*)^{-1}} \right]. \quad (3)$$

Here, $D^*$ represents the condensed dynamic stiffness matrix, i.e., on the left and right boundaries of the substructure [18]. Let us consider two consecutive substructures $k$ and $k + 1$ which are rigidly connected at a common interface. In this case, the coupling conditions write $q_R^k = q_L^{k+1}$ and $F_R^k = -F_L^{k+1}$, and Eq. (2) gives:

$$\left[ \begin{matrix} q_L^{k+1} \\ -F_L^{k+1} \end{matrix} \right] = S \left[ \begin{matrix} q_L^k \\ -F_L^k \end{matrix} \right] \quad \text{and also} \quad \left[ \begin{matrix} q_R^{k+1} \\ F_R^{k+1} \end{matrix} \right] = S \left[ \begin{matrix} q_R^k \\ F_R^k \end{matrix} \right], \quad (4)$$

which relates the displacement and force vectors between two consecutive substructures $k$ and $k + 1$.

It can be proven that the transfer matrix $S$ is symplectic [23], which means that its eigenvalues and eigenvectors come in pairs $(\mu_j, \phi_j)$ and $(\mu_j^* = 1/\mu_j, \phi_j^*)$ with $|\mu_j| < 1$. Numerical issues about the computation of the eigensolutions of $S$ are discussed in [22, 32, 36]. According to the Bloch's theorem [34], the eigenvalues of $S$ have the meaning of wave parameters which are given by $\mu_j = e^{-i\beta_j d}$ for the waves traveling to the right direction ($\beta_j$: wave number for a wave $j$; $d$: substructure length) and $\mu_j^* = e^{i\beta_j d}$ for the waves traveling to the left direction. Also, the eigenvectors $\phi_j$ and $\phi_j^*$ have the meaning of wave shapes for the waves traveling to the right and left directions, respectively. These are vectors of size $2n \times 1$, expressed by:

$$\phi_j = \begin{bmatrix} \phi_{qj} \\ \phi_{Fj} \end{bmatrix}, \quad \phi_j^* = \begin{bmatrix} \phi_{qj}^* \\ \phi_{Fj}^* \end{bmatrix}, \quad (5)$$

where $\phi_{qj}$ and $\phi_{qj}^*$ (resp. $\phi_{Fj}$ and $\phi_{Fj}^*$) are $n \times 1$ displacement (resp. force) component vectors for the right-going and left-going waves. In matrix form, these vectors give:

$$\mathbf{\Phi}_q = [\phi_{q1} \cdots \phi_{qn}], \quad \mathbf{\Phi}_q^* = [\phi_{q1}^* \cdots \phi_{qn}^*], \quad \mathbf{\Phi}_F = [\phi_{F1} \cdots \phi_{Fn}], \quad \mathbf{\Phi}_F^* = [\phi_{F1}^* \cdots \phi_{Fn}^*]. \quad (6)$$

## 2.2 Forced response

The analysis of the harmonic response of periodic structures of finite length — i.e., with a finite number $N$ of substructures — can be easily performed with the WFE method. The advantage, against the conventional FE method, is in the modeling of a whole periodic structure via a small matrix system of size $2n \times 2n$ that is linked to the number of DOFs on the left and right boundaries of a substructure. In comparison, the FE method would involve assembling the condensed dynamic stiffness matrices of all the substructures, yielding a bigger matrix system of size $(N+1)n \times (N+1)n$ whose computation is likely to take much more time. The main idea behind the WFE modeling is to expand the displacement and force vectors, on a particular interface ($k$) connecting two substructures $k-1$ and $k$, on the wave shapes $\phi_{qj}$ and $\phi_{qj}^*$, and $\phi_{Fj}$ and $\phi_{Fj}^*$, of a periodic structure [16, 18]:

$$q_L^{(k)} = q_R^{(k)} = \mathbf{\Phi}_q \mu^{k-1} \mathbf{Q} + \mathbf{\Phi}_q^* \mu^{N-k+1} \mathbf{Q}^* \quad k = 1, \dots, N+1, \quad (7)$$

$$-\mathbf{F}_L^{(k)} = \mathbf{F}_R^{(k)} = \mathbf{\Phi}_F \mu^{k-1} \mathbf{Q} + \mathbf{\Phi}_F^* \mu^{N-k+1} \mathbf{Q}^* \quad k = 1, \dots, N+1. \quad (8)$$
Here, $q_L^{(k)}$ and $F_L^{(k)}$ represent the displacement and force vectors at the left boundary of the substructure $k$ for $k = 1, \dots, N$; $q_R^{(k)}$ and $F_R^{(k)}$ represent the displacement and force vectors at the right boundary of the substructure $k-1$ for $k = 2, \dots, N+1$; $Q$ is the wave amplitude vector for the right-going waves at the substructure interface (1) (left end of the periodic structure); $Q^*$ is the wave amplitude vector for the left-going waves at the substructure interface ($N+1$) (right end of the periodic structure). Finally, $\mu = \text{diag}\{\mu_j\}_{j=1}^n$ is the diagonal matrix of the wave parameters $\mu_j$ for the right-going waves, with the property that $\|\mu\|_2 < 1$ [16].

Consider for the sake of clarity a finite periodic structure as shown in Fig. 2 whose left end is excited by some forces (vector $F_0$), and whose right end is constrained by prescribed displacement vector $q_0^*$. In this case, by considering Eqs. (7) and (8), the boundary conditions write $-F_0 = \Phi_F Q + \Phi_F^* \mu^N Q^*$ and $q_0^* = \Phi_q \mu^N Q + \Phi_q^* Q^*$ which, after some simple manipulations, yield the following wave-based matrix equation (size $2n \times 2n$):

$$ \left[ \begin{array}{cc} I_n & (\Phi_F)^{-1} \Phi_F^* \mu^N \\ (\Phi_q^*)^{-1} \Phi_q \mu^N & I_n \end{array} \right] \left[ \begin{array}{c} Q \\ Q^* \end{array} \right] = \left[ \begin{array}{c} -(\Phi_F)^{-1} F_0 \\ (\Phi_q^*)^{-1} q_0^* \end{array} \right]. \quad (9) $$

Solving the wave-based matrix equation (9) yields the wave amplitude vectors $Q$ and $Q^*$. The determination of the displacement and force vectors, at any substructure interface ($k$), follows from Eqs. (7) and (8). If needed, the displacements at the internal nodes (I) of the substructures can be retrieved by considering the full dynamic stiffness matrix, i.e., $q_I = -(D_{II})^{-1}(D_{II}q_L^{(k)} + D_{IR}q_R^{(k+1)})$ for a substructure enclosed between two interfaces ($k$) and ($k+1$).

<figure><img src="image_3.png" /><figcaption>Figure 2: Finite periodic structure with prescribed forces $F_0$ (left end) and prescribed displacements $q_0^*$ (right end).</figcaption></figure>

The main numerical tasks involved in the WFE method may be summed up as follows:

(i) Computation of the condensed dynamic stiffness matrix $D^*$ of a substructure;

(ii) Computation of the eigensolutions of the transfer matrix $S$;

(iii) Computation of the matrix system (9).

Task (iii) is not computationally cumbersome given that the matrix system (9) has a small size, as explained earlier. As for the computation of the condensed dynamic stiffness matrix $D^*$ (task (i)), this can be easily performed via the CB method [18]. For periodic structures involving substructures whose number of boundary DOFs is not excessive (e.g., smaller than 300), task (ii) is not cumbersome too. Things are more complicated for substructures with large-sized FE models, and in this case, it is advised to only compute a few low-order waves (rather than all the waves) using appropriate eigensolvers like the Lanczos algorithm [22].

The WFE method can be advantageously applied to structural assemblies made up of different periodic structures and other non-periodic structural components, i.e., mechanical systems which are of practical interest in engineering. The reader is referred to [23, 24, 17, 30]. The
WFE modeling of those structural assemblies involves considering coupling conditions between periodic structures (two or more), across a rigid interface c or a non-periodic structural component c (coupling element), which in wave-based form are written as:

$$Q^{c*} = \mathbb{C}^c Q^c + \mathbb{F}^c, \qquad (10)$$

where $Q^c$ and $Q^{c*}$ represent the wave amplitudes vectors for the outgoing waves and incoming waves, respectively, for the periodic structures connected at interface c or coupling element c; $\mathbb{C}^c$ is the so-called scattering matrix (reflection and transmission coefficients); $\mathbb{F}^c$ is a vector representing the excitation sources imposed at the interface or inside the coupling element. Hence, a global wave-based matrix equation which takes into account coupling conditions between several periodic structures, and classical boundary conditions (prescribed forces or displacements), can be proposed whose size is not big, i.e., linked to the numbers of boundary DOFs of the substructures used for modeling the periodic structures.

## 2.3 Some recent and ongoing works

1D periodic structures are often considered as straight structures, made up of simple or more complicated substructures, vibrating under harmonic disturbance. The purpose here is to show that the WFE method can handle more sophisticated engineering problems. Two examples of recent WFE strategies are proposed. The purpose is to give the main ideas and outlines of these strategies rather than fully detailing all the theoretical developments (see [19, 7]).

### 2.3.1 Periodic structures with cyclic symmetry

The first kind of problems concerns the analysis of multi-stage bladed disks, which is well known in the aeronautic industry. The proposed approach involves computing the dynamic flexibility modes of a given periodic structure (cyclic symmetry) via the WFE method, and modeling assemblies made up of several periodic structures via dynamic substructuring [19].

<figure><img src="image_4.png" /><figcaption>Figure 3: (left) Periodic structure with one of its substructures subject to one unit force on the boundary. (right) Two connected periodic structures with arbitrary forces applied to the coupling interface.</figcaption></figure>

A dynamic flexibility mode of a periodic structure (N substructures) represents the displacement vector that results from applying one unit force to a given DOF of the boundary, the others being free (see Fig. 3). For instance, assume that the unit force (force vector $\mathbb{F}_B$) is applied to a boundary DOF i of a substructure enclosed between two substructure interfaces (N) and (1).
The resulting dynamic flexibility mode is given by:

$$ \chi_i^{(N)} = \begin{bmatrix} (\mathbf{I}_{(n+n_{\mathrm{IB}})/3} \otimes \mathcal{R})\mathbf{q}^{(1)} \\ (\mathbf{I}_{(n+n_{\mathrm{IB}})/3} \otimes \mathcal{R}^2)\mathbf{q}^{(2)} \\ \vdots \\ (\mathbf{I}_{(n+n_{\mathrm{IB}})/3} \otimes \mathcal{R}^N)\mathbf{q}^{(N)} \end{bmatrix}, \quad (11) $$

where $\mathbf{q}^{(k)}$ ($k = 1, \dots, N-1$) represents the displacement vector of a substructure enclosed between two substructure interfaces ($k$) and ($k+1$) and expressed in the local cartesian coordinate system of the substructure; $\mathcal{R}$ is the $3 \times 3$ rotation matrix for the angle $-2\pi/N$, with $\mathcal{R}^N = \mathbf{I}_3$; $n$ is the number of interface DOFs (between the substructures), and $n_{\mathrm{IB}}$ is the number of internal DOFs which belong to the boundary. Recall that the unit force is applied to the boundary of a substructure enclosed between interfaces ($N$) and (1). The displacement vector of this substructure is denoted by $\mathbf{q}^{(N)}$ and is given by:

$$ \mathbf{q}^{(N)} = \begin{bmatrix} \mathbf{I}_n & \mathbf{0} \\ \tilde{\mathbf{X}}_{\mathrm{el}|\mathrm{B}} & \tilde{\mathbf{D}}_{\mathrm{el}-\mathrm{el}}^{-1} \tilde{\mathbf{D}}_{\mathrm{st}-\mathrm{el}}^{T} \\ \mathbf{X}_{\mathrm{st}|\mathrm{B}} & -\tilde{\mathbf{X}}_{\mathrm{el}|\mathrm{B}} \tilde{\mathbf{D}}_{\mathrm{el}-\mathrm{el}}^{-1} \tilde{\mathbf{D}}_{\mathrm{st}-\mathrm{el}}^{T} \end{bmatrix} \mathcal{H}^T \begin{bmatrix} \mathbf{\Phi}_q \boldsymbol{\mu}^{N-1} & \mathbf{\Phi}_q^* \\ \mathbf{\Phi}_q & \mathbf{\Phi}_q^* \boldsymbol{\mu}^{N-1} \end{bmatrix} \begin{bmatrix} \mathbf{Q} \\ \mathbf{Q}^* \end{bmatrix} + \begin{bmatrix} \mathbf{0} \\ \tilde{\mathbf{X}}_{\mathrm{el}|\mathrm{B}} \tilde{\mathbf{D}}_{\mathrm{el}-\mathrm{el}}^{-1} \tilde{\mathbf{X}}_{\mathrm{el}}^T \mathcal{L}_{\mathrm{IB}} \end{bmatrix} \mathbf{F}_{\mathrm{B}}. \quad (12) $$

Here, $\mathbf{X}_{\mathrm{st}}$ and $\tilde{\mathbf{X}}_{\mathrm{el}}$ are, respectively, the matrices of static modes and fixed interface modes of the substructure (CB method), $\tilde{\mathbf{D}}_{\mathrm{el}-\mathrm{el}}$ represents the dynamic stiffness matrix of the substructure “projected” on $\tilde{\mathbf{X}}_{\mathrm{el}}$, and $\mathcal{H}$ is a transformation matrix which relates the unit vectors in cylindrical coordinates (i.e., $\tilde{\mathbf{e}}_r$, $\tilde{\mathbf{e}}_\theta$ and $\tilde{\mathbf{e}}_z$) to those in cartesian coordinates (i.e., $\mathbf{e}_x$, $\mathbf{e}_y$ and $\mathbf{e}_z$) at the interface nodes of the substructure. Also, for a non-excited substructure, one has:

$$ \mathbf{q}^{(k)} = \begin{bmatrix} \mathbf{I}_n & \mathbf{0} \\ \tilde{\mathbf{X}}_{\mathrm{el}|\mathrm{B}} & \tilde{\mathbf{D}}_{\mathrm{el}-\mathrm{el}}^{-1} \tilde{\mathbf{D}}_{\mathrm{st}-\mathrm{el}}^{T} \end{bmatrix} \mathcal{H}^T \begin{bmatrix} \mathbf{\Phi}_q \boldsymbol{\mu}^{k-1} & \mathbf{\Phi}_q^* \boldsymbol{\mu}^{N-k} \\ \mathbf{\Phi}_q^* \boldsymbol{\mu}^k & \mathbf{\Phi}_q^* \boldsymbol{\mu}^{N-k-1} \end{bmatrix} \begin{bmatrix} \mathbf{Q} \\ \mathbf{Q}^* \end{bmatrix} \quad \text{for } k \neq N. \quad (13) $$

In Eqs. (12) and (13), $\mathbf{Q} = \mathbf{Q}^{(1)}$ and $\mathbf{Q}^* = \mathbf{Q}^{*(N)}$ are wave amplitude vectors, expressed at substructure boundaries (1) and (N) (see Fig. 3) and obtained via WFE procedure. The matrix of flexibility modes which results from applying successive unit forces to the boundary DOFs of the substructure enclosed between interfaces (N) and (1) is given by $\mathbf{X}^{(N)} = [\chi_1^{(N)} \chi_2^{(N)} \cdots \chi_{n_B}^{(N)}]$ where $n_B$ is the number of boundary DOFs of a substructure. The whole matrix of flexibility modes which results from applying successive unit forces to the boundary DOFs of the whole periodic structure follows as $\mathbf{X} = [\mathbf{X}^{(1)} \mathbf{X}^{(2)} \cdots \mathbf{X}^{(N)}]$, where:

$$ \mathbf{X}^{(k)} = (\mathbf{I}_{n_T/3} \otimes \mathcal{R}^k) \mathcal{P}^{(k)} \mathbf{X}^{(N)} (\mathbf{I}_{n_B/3} \otimes \mathcal{R}^k)^T, \quad (14) $$

where $\mathcal{P}^{(k)}$ is a permutation matrix, and $n_T = N(n + n_{\mathrm{IB}})$. The modeling of several periodic structures with cyclic symmetry follows from domain decomposition procedure. For instance, for two structures (displacement vectors $\mathbf{q}_1$ and $\mathbf{q}_2$) subject to force vectors $\mathbf{F}_{\mathrm{B}1}$ and $\mathbf{F}_{\mathrm{B}2}$ on their boundaries, this yields:

$$ \begin{bmatrix} \mathbf{I} & \mathbf{0} & -\mathbf{X}_1 \mathcal{L}_{\mathrm{B}1}^T \\ \mathbf{0} & \mathbf{I} & \mathbf{X}_2 \mathcal{L}_{\mathrm{B}2}^T \\ \mathcal{L}_1 & -\mathcal{L}_2 & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{q}_1 \\ \mathbf{q}_2 \\ \mathbf{\lambda} \end{bmatrix} = \begin{bmatrix} \mathbf{X}_1 \mathbf{F}_{\mathrm{B}1} \\ \mathbf{X}_2 \mathbf{F}_{\mathrm{B}2} \\ \mathbf{0} \end{bmatrix}, \quad (15) $$
where $X_1$ and $X_2$ are the matrices of flexibility modes of the structures, and $\mathcal{L}_1, \mathcal{L}_2, \mathcal{L}_{B1}, \mathcal{L}_{B2}$ are localization matrices. Also, $\lambda$ represents a vector of Lagrange multipliers to enforce the displacement continuity conditions (i.e., $\mathcal{L}_1 q_1 - \mathcal{L}_2 q_2 = 0$) for some nodes of the boundaries of the structures. Solving Eq. (15) yields:

$$ \lambda = -(\mathcal{L}_1 X_1 \mathcal{L}_{B1}^T + \mathcal{L}_2 X_2 \mathcal{L}_{B2}^T)^{-1} (\mathcal{L}_1 X_1 F_{B1} - \mathcal{L}_2 X_2 F_{B2}). \quad (16) $$

The displacement vectors of the structures follow as $q_1 = X_1(F_{B1} + \mathcal{L}_{B1}^T \lambda)$ and $q_2 = X_2(F_{B2} - \mathcal{L}_{B2}^T \lambda)$.

Some results are proposed hereafter which concern two coupled bladed disks 1 and 2 (two stages) with $N_1 = 36$ and $N_2 = 60$ substructures (respectively), excited by non-symmetric forces. Here, the harmonic response function (transverse displacement) at some node at the tip of one blade of structure 1 is assessed up to 10,000 Hz (see Fig. 4), which is high enough to include local resonance effects for each blade. Those results show that the proposed approach perfectly agrees with the conventional FE method.

<figure><img src="image_5.png" /><figcaption>Figure 4: (left) Two periodic structures and FE meshes of the related substructures ($\Gamma_-$ and $\Gamma_+$ are coupling interfaces with previous and subsequent substructures). (right) Response functions, conventional FE (black solid line) and WFE (pink dotted line).</figcaption></figure>

### 2.3.2 Time response of periodic structures via wave-based absorbing boundary conditions

The second kind of problems concerns the modeling of infinite periodic structures (or periodic structures which are supposed to be long enough) subject to localized time-dependent loadings. Potential applications could be, for instance, structures like pipelines or tanks subject to shocks. A periodic structure with a finite number $N$ of substructures and subject to time dependent forces (vector $F(t)$) is shown in Fig. 5. The semi-infinite behavior at the left and right ends — i.e., the fact that no wave is coming from infinity — is described by means of absorbing boundary conditions expressed in the time domain [7]. Although easy to understand in the frequency domain, the determination of absorbing conditions for time response analysis is not a straightforward task and, for this reason, this study hasn't been conducted extensively so far.

The proposed approach can be summarized as follows. Let us start with the impedance matrices (frequency domain) at the left (L) and right (R) ends of a periodic structure ($N$ substructures), which can be simply obtained via the WFE method as follows:

$$ Z_L = -\Phi_F^*(\Phi_q^*)^{-1}, \quad Z_R = \Phi_F(\Phi_q)^{-1}. \quad (17) $$
<figure><img src="image_6.png" /><figcaption>Figure 5: Schematic of a periodic structure with N substructures subject to time-dependent forces and semi-infinite boundary conditions.</figcaption></figure>

Then, decompose the impedance matrices $Z_L$ and $Z_R$ via rational approximations:

$$ 
\begin{aligned}
Z_L &= \sum_{k=1}^{Q} 2 \frac{i \omega \Re\{\mathbf{R}_{L(2k)}\} - \Re\{\overline{p_{L(2k)}} \mathbf{R}_{L(2k)}\}}{-\omega^2 - 2i \omega \Re\{p_{L(2k)}\} + |p_{L(2k)}|^2} + \sum_{k=2Q+1}^{P} \frac{\mathbf{R}_{Lk}}{i\omega - p_{Lk}} + \mathbf{K}_L, \\
Z_R &= \sum_{k=1}^{Q} 2 \frac{i \omega \Re\{\mathbf{R}_{R(2k)}\} - \Re\{\overline{p_{R(2k)}} \mathbf{R}_{R(2k)}\}}{-\omega^2 - 2i \omega \Re\{p_{R(2k)}\} + |p_{R(2k)}|^2} + \sum_{k=2Q+1}^{P} \frac{\mathbf{R}_{Rk}}{i\omega - p_{Rk}} + \mathbf{K}_R,
\end{aligned} 
\quad (18) $$

where $(p_{Lk}, p_{Rk})$ and $(\mathbf{R}_{Lk}, \mathbf{R}_{Rk})$ denote poles and matrices of residues ($k = 1, \dots, P$), respectively. Usually, some of these terms appear in conjugate pairs, say, for $k = 1, \dots, Q$. In the frequency domain, the absorbing boundary conditions are expressed by $\mathbf{F}_L = Z_L q_L$ and $\mathbf{F}_R = Z_R q_R$ where $q_L$ (resp. $q_R$) and $\mathbf{F}_L$ (resp. $\mathbf{F}_R$) are the displacement vector and force vector at the left (resp. right) end of the structure, respectively. By considering vectors of supplementary variables $\mathbf{X}_{Lk}$ and $\mathbf{X}_{Rk}$ together with Eq. (18), it can be shown that the absorbing conditions can be rewritten as follows:

$$ 
\begin{aligned}
\mathbf{F}_L &= \sum_{k=1}^{Q} 2 \left( i\omega \Re\{\mathbf{R}_{L(2k)}\} - \Re\{\overline{p_{L(2k)}} \mathbf{R}_{L(2k)}\} \right) \mathbf{X}_{Lk} + \sum_{k=2Q+1}^{P} \mathbf{R}_{Lk} (i\omega) \mathbf{X}_{L(k-Q)} + \mathbf{K}_L \mathbf{q}_L, \\
\mathbf{F}_R &= \sum_{k=1}^{Q} 2 \left( i\omega \Re\{\mathbf{R}_{R(2k)}\} - \Re\{\overline{p_{R(2k)}} \mathbf{R}_{R(2k)}\} \right) \mathbf{X}_{Rk} + \sum_{k=2Q+1}^{P} \mathbf{R}_{Rk} (i\omega) \mathbf{X}_{R(k-Q)} + \mathbf{K}_R \mathbf{q}_R,
\end{aligned} 
\quad (19) $$

where:

$$ 
\begin{aligned}
(-\omega^2 - 2i\omega\Re\{p_{L(2k)}\} + |p_{L(2k)}|^2) \mathbf{X}_{Lk} &= \mathbf{q}_L \quad \text{for } k = 1, \dots, Q, \\
(-\omega^2 - 2i\omega\Re\{p_{R(2k)}\} + |p_{R(2k)}|^2) \mathbf{X}_{Rk} &= \mathbf{q}_R \quad \text{for } k = 1, \dots, Q, \\
(-\omega^2 - i\omega p_{Lk}) \mathbf{X}_{L(k-Q)} &= \mathbf{q}_L \quad \text{for } k = (2Q+1), \dots, P, \\
(-\omega^2 - i\omega p_{Rk}) \mathbf{X}_{R(k-Q)} &= \mathbf{q}_R \quad \text{for } k = (2Q+1), \dots, P.
\end{aligned} 
\quad (20) $$

The interesting feature of Eq. (19) is that the force vectors are described in terms of polynomials of $i\omega$, of order 1. In comparison, Eq. (18) is more complicated with rational terms of $i\omega$. As for Eq. (20)—i.e., the relations between the displacement vectors and the vectors of supplementary variables — it also involves polynomials of $i\omega$ (up to order 2). By separating the terms of identical powers of $i\omega$ in Eqs. (19) and (20), and by invoking the classical time-frequency transforms $\mathbf{q}(\omega) \to \mathbf{q}(t)$, $i\omega\mathbf{q} \to \dot{\mathbf{q}}$, $-\omega^2\mathbf{q} \to \ddot{\mathbf{q}}$ and $\mathbf{X}(\omega) \to \mathbf{X}(t)$, $i\omega\mathbf{X} \to \dot{\mathbf{X}}$, $-\omega^2\mathbf{X} \to \ddot{\mathbf{X}}$
(where dot and double-dot notations mean single and double time derivatives, respectively), two classical second-order differential equations can be obtained:

$$ M_L \begin{bmatrix} \ddot{q}_L \\ \ddot{X}_L \end{bmatrix} + C_L \begin{bmatrix} \dot{q}_L \\ \dot{X}_L \end{bmatrix} + K_L \begin{bmatrix} q_L \\ X_L \end{bmatrix} = \begin{bmatrix} F_L \\ 0 \end{bmatrix}, \quad M_R \begin{bmatrix} \ddot{q}_R \\ \ddot{X}_R \end{bmatrix} + C_R \begin{bmatrix} \dot{q}_R \\ \dot{X}_R \end{bmatrix} + K_R \begin{bmatrix} q_R \\ X_R \end{bmatrix} = \begin{bmatrix} F_R \\ 0 \end{bmatrix}, \quad (21) $$

where $X_L$ and $X_R$ are vectors built from all the vectors of supplementary variables $X_{Lk}$ and $X_{Rk}$. Eq. (21) can be simply integrated into a FE model of a periodic structure (N substructures). Thus, the FE model of a periodic structure with absorbing boundary conditions can be written in the time domain as follows:

$$ M_{tot} \ddot{y} + C_{tot} \dot{y} + K_{tot} y = F_{tot}, \quad (22) $$

where

$$ y = y(t) = \begin{bmatrix} q_I(t) \\ q_L(t) \\ q_R(t) \\ X_L(t) \\ X_R(t) \end{bmatrix}, \quad F_{tot} = F_{tot}(t) = \begin{bmatrix} F_I(t) \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}, \quad (23) $$

where $q_I$ and $F_I$ represent the displacement vector and force vector (time loadings) for the internal DOFs (between the left and right ends). Eq. (22) is a classical second-order differential equation which can be solved in a standard way via time integration numerical schemes (e.g., Newmark scheme). For instance, the time response of an infinite periodic structure, with 2D substructures, subject at its center to a shock (transverse Ricker wavelet force pulse) is analyzed. Here, a periodic structure consisting of only 20 substructures (squares with holes, meshed with 2D plane stress triangular elements) is considered. Absorbing boundary conditions are invoked at the left and right ends of the periodic structure, as discussed earlier. For comparison purposes, a FE model of a much longer periodic structure with free ends and 200 substructures, aiming at modeling an infinite structure, is simulated. Results are shown in Fig. 6 for the time history and the displacement field at $t = 0.01$ s, which fully validate the proposed approach.

<figure><img src="image_7.png" /><figcaption>Figure 6: Time response of the periodic structure: (left) time history at the right end (reference FE: red line; proposed approach: blue crosses); (right) deformed shape.</figcaption></figure>

## 2.4 Some results about locally resonant metamaterials

The analysis and design of locally resonant metamaterials is a research topic of growing interest, which is mostly due to their ability to create band gaps and strongly attenuate the vibration levels of structures [31, 1, 35, 8, 33]. Band gaps refer to frequency bands in which
one or several waves do not propagate, and which are here induced by the local resonances of the substructures (cells).

Among the resonant metamaterials are those made up of multi-layered continuous substruc- tures combining soft and heavy materials (e.g., soft and stiff/heavy layers/core) which behave like a single or multi-DOF mass-spring system. For instance, a periodic structure with periodic inclusions of viscoelastic layers (spring and damping effects) and heavy parts (mass effect) is shown in Fig. 7. The goal of the present study is to analyze the band gap effect in such periodic structures. Also, it aims at investigating whether periodic structures with viscoelastic proper- ties can be used to create band gaps with larger bandwidths (compared to when purely elatic properties are considered), to make them useful for a wider range of engineering applications.

<figure><img src="image_8.png" /><figcaption>Figure 7: Periodic structure with N = 10 squared substructures with two viscoelastic layers (rubber), a heavy layer and a heavy core (tungsten).</figcaption></figure>

Following the WFE procedure, the wave properties of the periodic structure and the related forced response can be obtained straightforwardly. Attention must be paid, however, to the FE modeling of the viscoelastic layers and the fact that their stiffness matrices depend on the frequency (see [21] for further details):

$$
\mathbf{K}_{v}^{\text{vel}}(\omega) = \mathbf{K}_{v\infty} + \sum_{i=1}^{n_v} \mathbf{K}_{vi} \left( \frac{\omega^2 \tau_i^2}{1 + \omega^2 \tau_i^2} + i \frac{\omega \tau_i}{1 + \omega^2 \tau_i^2} \right), \quad (24)
$$

where $\tau_i$ are relaxation times with $\tau_i = 10^{-i}$ s for $i = 1, 2, \dots, n_v$. The stiffness matrix of a whole substructure involves assembling several FE viscoelastic layers ($v$) together with FE elastic layers ($e$), i.e.:

$$
\mathbf{K}(\omega) = \mathbf{K}_1 + \mathbf{K}_2(\omega), \tag{25}
$$

where

$$
\mathbf{K}_1 = \sum_{e=1}^{N^{\text{el}}} (1 + i\eta_e^{\text{el}})(\mathcal{L}_e^{\text{el}})^T \mathbf{K}_e^{\text{el}} \mathcal{L}_e^{\text{el}} + \sum_{v=1}^{N^{\text{vel}}} (\mathcal{L}_v^{\text{vel}})^T \mathbf{K}_{v\infty} \mathcal{L}_v^{\text{vel}}, \quad (26)
$$

and

$$
\mathbf{K}_2(\omega) = \sum_{v=1}^{N^{\text{vel}}} \sum_{i=1}^{n_v} (\mathcal{L}_v^{\text{vel}})^T \mathbf{K}_{vi} \mathcal{L}_v^{\text{vel}} \left( \frac{\omega^2 \tau_i^2}{1 + \omega^2 \tau_i^2} + i \frac{\omega \tau_i}{1 + \omega^2 \tau_i^2} \right), \quad (27)
$$

where $\mathcal{L}_e^{\text{el}}$ and $\mathcal{L}_v^{\text{vel}}$ are localization (Boolean) matrices, and $\eta_e^{\text{el}}$ denotes the loss factors of the elastic parts. As it turns out, the dynamic stiffness matrix of the substructure is expressed by $\mathbf{D} = -\omega^2\mathbf{M} + \mathbf{K}_1 + \mathbf{K}_2(\omega)$. Numerical issues about the condensation of the dynamic stiffness matrix on the left and right boundaries of the substructure (see Sec. 2.1) are discussed in [21]. Results about the dispersion curves of the bending wave are shown in Fig. 8, e.g., when $n_v = 5$
(number of terms in the series in Eq. (24)) and $K_{vi} = 0.05 \times K_{v\infty}$ ($i = 1, \dots, 5$). The structure of interest is shown in Fig. 7. Here, squared substructures with two viscoelastic layers (rubber, $N^{\text{vel}} = 2$) and elastic parts including an external layer in aluminum, a central layer and a core in tungsten ($N^{\text{el}} = 3$), are considered. Here, each substructure is meshed using 20 x 20 eight-node plane stress quad elements. For comparison purposes, a purely equivalent elastic structure, when $K_v^{\text{vel}} = K_{v\infty}(1 + 5 \times 0.05)$, is also considered.

Band gaps occur when the real parts of the wavenumbers are equal to 0 or $\pi/d$ ($d$ being the substructure length). The dispersion curves displayed in Fig. 8 represent plots of the functions $\omega \mapsto \Re\{\beta d\}$ and $\omega \mapsto \Im\{\beta d\}$ where $\beta$ is the wavenumber for the bending wave. Associated to the band gap behavior — when $\Re\{\beta d\} = 0$ or $\Im\{\beta d\} = \pi$ — are local increases of the magnitude of the imaginary part of the wavenumber, which means that the wave is becoming evanescent. The band gap effect in resonant metamaterials mostly relies on the presence of substructure local modes. For instance, the fixed interface mode shapes of a purely equivalent elastic substructure are shown in Fig. 8, which could help understand the physical behavior of the structure in the vicinity of band gaps. As expected, band gap phenomena yield strong reduction of the vibration levels of the structure. Here, a structure with $N = 10$ substructures with prescribed transverse harmonic displacement at the right end is considered. For instance, two locally resonant band gaps appear at 742 Hz and 984 Hz which are close to each other. In the purely elastic case, the vibration levels show two local minima around these frequencies which appear to be well separated to each other. The interesting feature in the viscoelastic case is that these frequency bands with minimum vibration levels tend to merge to produce a wide frequency band (i.e., between 700 Hz and 1100 Hz) where, overall, the vibration levels remain low. This opens interesting prospects for the passive vibration control of structures.

<figure><img src="image_9.png" /><figcaption>Figure 8: Periodic structure with $N = 10$ substructures. (top) Dispersion curves (elastic: blue line, viscoelastic: green line); (middle) Quadratic velocity at the left end; (bottom) fixed interface modes of the substructures.</figcaption></figure>

Additional results about 1D periodic structures and resonant or non-resonant (Bragg) band gaps can be found in [30] (aircraft structures) and [2] (rotating phononic crystal ring).
# 3 2D PERIODIC STRUCTURES

## 3.1 Wave approaches and the modeling of structures of finite dimensions

The dynamic analysis of 2D periodic structures like metamaterial plates with resonant sub-structures constitutes a research topic of growing interest which mostly rely, with the era of 3D printing, on the ability to design and manufacture such structures at a low cost. A schematic of a 2D periodic structure made up of resonant layered substructures (soft layer and heavy core) is shown in Fig. 9. From the FE point of view, there exist plenty of works analyzing the propagation of Bloch waves in infinite 2D periodic structures and the related band gap effects. These approaches can be classified in two ways:

1. Those where a Bloch solution $u(x, \beta, t) = \tilde{u}(x, \beta)e^{i(\omega t - \beta^T x)}$ ($\beta = [\beta_x \ \beta_y]^T$: wavenumber vector with $\beta_x = \beta \cos(\theta)$ and $\beta_y = \beta \sin(\theta)$; $x = [x \ y \ z]^T$: position vector; $\tilde{u}(x, \beta)$: periodic function) is considered in the dynamic equilibrium equation of a cell (see for instance [9, 4]). The resulting dynamic equation of the cell together with the periodicity conditions for $\tilde{u}(x, \beta)$ are then converted into a FE model (Galerkin method), yielding a generalized eigenproblem of the form $A(\beta_x, \beta_y)\tilde{q} = \omega^2 B(\beta_x, \beta_y)\tilde{q}$ for the eigenvalue $\omega^2$ and the eigenvector $\tilde{q}$. The computation of the so-called band diagrams then involves plotting the pulsation $\omega$ for some pre-determined values of $\beta_x$ and $\beta_y$ following the first irreducible Brillouin zone.

2. Those where a FE model of a substructure is considered and where pseudo-periodic conditions are applied to its boundary. For substructures with rectangular boundaries (e.g., of lengths $L_x$ and $L_y$), the pseudo-periodic conditions write $q_R = e^{-i\beta_x L_x} q_L$ between the left (L) and right (R) edges, $q_T = e^{-i\beta_y L_y} q_B$ between the bottom (B) and top (T) edges, $q_{RB} = e^{-i\beta_x L_x} q_{LB}$ between the left-bottom and right-bottom corners, $q_{LT} = e^{-i\beta_y L_y} q_{LB}$ between the left-bottom and left-top corners, and $q_{RT} = e^{-i(\beta_x L_x + \beta_y L_y)} q_{LB}$ between the left-bottom and right-top corners (see for instance [14, 25]). By considering the FE model of a substructure together with these pseudo-periodic conditions, a generalized eigenproblem of the form $A(\beta_x, \beta_y)\tilde{q} = \omega^2 B(\beta_x, \beta_y)\tilde{q}$ can be proposed (again) for computing the band diagrams of the structure.

<figure><img src="image_10.png" /><figcaption>Figure 9: 2D periodic structure and FE mesh of a substructure (red spots highlight boundary nodes).</figcaption></figure>

Although well detailed in the literature, the dynamic analysis of 2D periodic structures that uses wave approaches appears quite restricted to infinite systems where waves are supposed to
propagate freely, or infinite systems subject to point forces (see for instance [5]). In contrast, the analysis of 2D periodic structures of finite dimensions — i.e., systems which are frequently encountered in engineering applications — is not well reported. The thing is, the modeling of 2D periodic structures of finite dimensions via wave approaches appears to be a tough task (as opposed to the 1D periodic case) mostly due to the multi-directional behavior of the propagating waves, and the description of the boundary conditions (or coupling conditions) using these multi-directional waves.

Fortunately, there remain several FE-based reduced models which can be used to simulate the dynamic behavior of finite periodic structures, at a low computational cost. Indeed, by considering the periodicity hypothesis into account, the FE models of the substructures are identical and, therefore, they can be efficiently reduced via component mode synthesis like the Craig Bampton (CB) method. In fact, since the substructures are identical, their component modes are identical too which means that they only have to be computed once, i.e., regardless of the number of substructures considered. Hence, the following transformation matrix, for all the substructures, can be proposed [20]:

$$ \tilde{T} = \begin{bmatrix} \mathbf{I}_{N_B} & \mathbf{0} \\ -(\mathbf{K}_{II})^{-1}\mathbf{K}_{IB} & \tilde{\mathbf{X}} \end{bmatrix}, \qquad (28) $$

where $-(\mathbf{K}_{II})^{-1}\mathbf{K}_{IB}$ and $\tilde{\mathbf{X}}$ are the matrix of static modes and the matrix of fixed interface modes (a reduced set) of the substructures, respectively; $N_B$ is the number of boundary DOFs of the substructures. The reduced dynamic equation of a substructure $s$ follows as:

$$ \left[ (-\omega^2 + i\omega a)\tilde{\mathbf{M}} + (i\omega b + 1)\tilde{\mathbf{K}} \right] \begin{bmatrix} \tilde{\mathbf{u}}_B^s \\ \tilde{\alpha}^s \end{bmatrix} = \begin{bmatrix} \mathbf{F}_B^s - \mathbf{K}_{BI}(\mathbf{K}_{II})^{-1}\mathbf{F}_I^s \\ \tilde{\mathbf{X}}^T \mathbf{F}_I^s \end{bmatrix}, \qquad (29) $$

where notation “u” is used to denote displacement vectors (instead of “q”); $\tilde{\mathbf{u}}_B^s$ is the approximated displacement vector for the boundary nodes; $\tilde{\alpha}^s$ is a vector of generalized coordinates; $\mathbf{F}_B^s$ and $\mathbf{F}_I^s$ are force vectors for the boundary and internal nodes (respectively); $a$ and $b$ are damping coefficients (Rayleigh damping). Also, $\tilde{\mathbf{M}}$ and $\tilde{\mathbf{K}}$ are the reduced mass and stiffness matrices of the substructure (similar for all the substructures), expressed by:

$$ \tilde{\mathbf{M}} = \tilde{\mathbf{T}}^T \mathbf{M} \tilde{\mathbf{T}}, \quad \tilde{\mathbf{K}} = \tilde{\mathbf{T}}^T \mathbf{K} \tilde{\mathbf{T}}. \qquad (30) $$

The reduced global mass and stiffness matrices of a whole 2D periodic structure — namely, $\tilde{\mathbf{M}}_{\text{per}}$ and $\tilde{\mathbf{K}}_{\text{per}}$ — follow from conventional FE assembly procedures by considering the displacement continuity conditions at the boundary DOFs. The global dynamic equation of the periodic structure follows as:

$$ \left[ (-\omega^2 + i\omega a)\tilde{\mathbf{M}}_{\text{per}} + (i\omega b + 1)\tilde{\mathbf{K}}_{\text{per}} \right] \begin{bmatrix} (\tilde{\mathbf{u}}_B)_{\text{per}} \\ \tilde{\alpha}_{\text{per}} \end{bmatrix} = \tilde{\mathbf{F}}_{\text{per}}, \qquad (31) $$

where $(\tilde{\mathbf{u}}_B)_{\text{per}}$ is the displacement vector of the boundary/interface nodes of the substructures. In addition, the boundary conditions of the periodic structure may be written in a general way as:

$$ \mathbf{Y}_B(\tilde{\mathbf{u}}_B)_{\text{per}} + \mathbf{Z}_B(\tilde{\mathbf{F}}_B)_{\text{per}} = \mathbf{Y}_B^0(\tilde{\mathbf{u}}_B)_{\text{per}}^0 + \mathbf{Z}_B^0(\tilde{\mathbf{F}}_B)_{\text{per}}^0, \qquad (32) $$

where $(\tilde{\mathbf{u}}_B)_{\text{per}}^0$ and $(\tilde{\mathbf{F}}_B)_{\text{per}}^0$ are vectors of prescribed displacements and prescribed forces, respectively, and where $\mathbf{Y}_B$, $\mathbf{Z}_B$, $\mathbf{Y}_B^0$ and $\mathbf{Z}_B^0$ are matrices whose expressions depend on the kind of applications considered. A second reduced mode expansion can be proposed to speed up the computation of the matrix system (31) by expressing the displacement vector $(\tilde{\mathbf{u}}_B)_{\text{per}}$ on some “boundary modes” (eigenvectors) for the matrix pencil $((\tilde{\mathbf{K}}_{BB})_{\text{per}}, (\tilde{\mathbf{M}}_{BB})_{\text{per}})$.
## 3.2 2D nearly periodic structures

It is well known that nearly periodic structures — e.g., periodic structures with disordered resonant substructures with slightly varying geometrical properties — possess interesting properties, e.g., the fact that the vibrational energy is localized around the excitation points [3, 13]. For instance, a schematic of a nearly periodic plate with layered substructures (soft layer and heavy core), embedded in a floor panel, is shown in Fig. 10.

<figure><img src="image_11.png" /><figcaption>Figure 10: 2D nearly periodic plate embedded in a floor panel, and related transverse displacement field resulting from harmonic excitations.</figcaption></figure>

A strategy for modeling such nearly periodic structures has been recently proposed in [20] and is briefly summarized here. The key idea for describing disordered substructures is to make use of distorted FE meshes, with node positions which randomly vary around a baseline (undistorted) mesh:

$$x_j^{se} = x_j^e + \epsilon_x^s f_x(x_j^e, y_j^e), \quad y_j^{se} = y_j^e + \epsilon_y^s f_y(x_j^e, y_j^e), \qquad (33)$$

where $(x_j^{se}, y_j^{se})$ and $(x_j^e, y_j^e)$ are the node coordinates of the distorted and undistorted meshes, respectively; $\epsilon_x^s$ and $\epsilon_y^s$ are uniform random variables with support $[-\delta, \delta]$ where $\delta$ is a dispersion parameter; $f_x(x, y)$ and $f_y(x, y)$ are two arbitrary deterministic functions of $(x, y)$, identical for all the substructures, which are supposed to cancel out on the boundary. Examples of distorted substructures are highlighted in Fig. 10.

When it comes to nearly periodic structures, the modeling of substructures raises two main issues compared to the purely periodic case. First, since the substructures are not identical, their mass and stiffness matrices $M^s$ and $K^s$ have to be computed several times (as many times as the number of substructures considered); second, the matrices of static modes and fixed interface modes $-(K_{II}^s)^{-1}K_{IB}^s$ and $\tilde{X}^s$, and also the reduced mass and stiffness matrices $\tilde{M}^s$ and $\tilde{K}^s$ (matrix multiplications with the transformation matrix, see Eq. (30)), have to be computed several times. For structures with many substructures, this makes the CB method ineffective. To address these issues — i.e., to avoid computing the reduced mass and stiffness matrices many times — it is proposed to compute these matrices at some points $\epsilon_x^s = (\epsilon_x)_p$ and $\epsilon_y^s = (\epsilon_y)_p$ (a small number) for some particular distorted FE meshes of the substructures, and approximate these matrices between these points via matrix interpolation for any distortion $\epsilon_x^s$ and $\epsilon_y^s$. Although easy to understand, the procedure for interpolating matrices is not straightforward, however, and requires the reduced matrices to be expressed using coordinate systems which are
compatible to each other [26]. For this task, the following alternative reduced matrices of fixed interface modes, at the interpolation points, can be considered [20]:

$$ \hat{\mathbf{X}}_p = \tilde{\mathbf{X}}_p (\mathbf{\Psi}^T \tilde{\mathbf{X}}_p)^{-1}, \qquad (34) $$

where

$$ \mathbf{\Psi} = \left( (\mathbf{M}_{\mathrm{II}}^{0})^{\frac{1}{2}} \right)^{T} \tilde{\mathbf{X}}^{0}. \qquad (35) $$

Here, notations $\mathbf{M}_{\mathrm{II}}^{0}$ and $\tilde{\mathbf{X}}^{0}$ mean that the matrices $\mathbf{M}_{\mathrm{II}}^{s}$ and $\tilde{\mathbf{X}}^{s}$ are expressed at $\epsilon_{x}^{s} = 0$ and $\epsilon_{y}^{s} = 0$. As a result, the following transformation matrix, reduced mass matrix and reduced stiffness matrix can be proposed:

$$ \hat{\mathbf{T}}_p = \begin{bmatrix} \mathbf{I}_{N_B} & \mathbf{0} \\ -(\mathbf{K}_{\mathrm{II}})_p^{-1} (\mathbf{K}_{\mathrm{IB}})_p & \hat{\mathbf{X}}_p \end{bmatrix}, \qquad (36) $$

$$ \hat{\mathbf{M}}_p = \hat{\mathbf{T}}_p^T \mathbf{M}_p \hat{\mathbf{T}}_p, \quad \hat{\mathbf{K}}_p = \hat{\mathbf{T}}_p^T \mathbf{K}_p \hat{\mathbf{T}}_p. \qquad (37) $$

The determination of the interpolated reduced mass and stiffness matrices — namely, $\hat{\mathbf{M}}^{s}$ and $\hat{\mathbf{K}}^{s}$ — for any substructure $s$ with mesh parameters $(\xi^{s}, \eta^{s})$ follows from classic interpolation, e.g., using an interpolation scheme based on eight interpolation points $(\xi_{p}, \eta_{p})$ and eight Serendipity interpolation functions $N_{p}(\xi^{s}, \eta^{s})$, as shown in Fig. 11. Thus:

$$ \hat{\mathbf{M}}^{s} = \sum_{p=1}^{8} N_{p}(\xi^{s}, \eta^{s}) \hat{\mathbf{M}}_{p}, \quad \hat{\mathbf{K}}^{s} = \sum_{p=1}^{8} N_{p}(\xi^{s}, \eta^{s}) \hat{\mathbf{K}}_{p}, \qquad (38) $$

where:

$$ \xi^{s} = \sqrt{2} \frac{\epsilon_{x}^{s}}{\delta}, \quad \eta^{s} = \sqrt{2} \frac{\epsilon_{y}^{s}}{\delta}. \qquad (39) $$

The modeling of a whole 2D nearly periodic structure follows from classic FE assembly procedure, as reported earlier.

Numerical experiments can be carried out to validate the proposed strategy. Let us consider, for instance, a nearly periodic plate with 8 × 4 substructures embedded in a floor panel and subject to random forces, as shown in Fig. 12 where simply supported boundary conditions are considered. Typical distorted substructures are shown in Fig. 10 and are meshed using 20 × 20 isoparametric quadratic Mindlin elements. Concerning the distorted meshes of the substructures, $f_{x}(x_{j}^{e}, y_{j}^{e})$ and $f_{y}(x_{j}^{e}, y_{j}^{e})$ are chosen as trigonometric functions, while $\delta = 0.01$ m (mesh dispersion).

The frequency responses of the structure (quadratic velocities at the location of the excitation points) over a frequency band of [0, 150] Hz, obtained from the interpolation strategy and the conventional FE modeling, are shown in Fig. 13 where good matches are observed. The eigenfrequencies of the undistorted substructures (fixed interfaces) are shown with dotted vertical red lines. For the purely periodic case, these frequencies correspond to band gap effects where vibration levels are low.

Also, the transverse displacement fields of the plate/panel at 125 Hz are shown in Fig. 14 for the nearly and purely periodic cases. For the nearly periodic case, the displacement field shows localized peaks of high magnitude in the plate, and small displacement levels in the floor. In comparison, the periodic case involves uniform (diffuse) field in the plate, and higher displacement levels in the floor. In this sense, the nearly periodic plate appears to be a suitable means to passively control the vibration levels of the floor. Finally, in terms of computational times, the interpolation strategy provides an overall time reduction of 82% compared to the FE method.
<figure><img src="image_12.png" /><figcaption>Figure 11: Eight point interpolation scheme based on serendipity interpolation functions $N_p(\xi^s, \eta^s)$ ($p = 1, \dots, 8$).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Schematic of a nearly periodic plate with 8 × 4 substructures with a floor panel subject to random forces.</figcaption></figure>

4 CONCLUSION

An overview of some wave approaches and FE-based procedures, for modeling periodic structures and metamaterials, has been proposed. For 1D periodic structures, the WFE method has been investigated and some recent advances in this field have been presented. For 2D periodic structures, a reduced FE modeling of the substructures based on the CB method has been proposed. The procedure has been extended to the dynamic analysis of 2D nearly periodic structures (with substructures having randomly varying geometrical properties) via the use of distorted FE meshes and the interpolation of the reduced matrices of the substructures between “interpolation points” (some particular distorted meshes). Those wave approaches and FE procedures have been proven relevant for predicting the dynamic behavior of periodic structures including metamaterials with resonant cells that could be of various shapes. The dynamic analysis of periodic structures is a research topic in constant progress. This concerns the wave
<figure><img src="image_14.png" /><figcaption>Figure 13: FRF (mean value of the quadratic velocity) of the plate with 8 × 4 substructures with the floor panel. (black solid line) FE method, nearly periodic structure; (red dotted line) interpolation strategy, nearly periodic structure; (blue dotted line) FE method, purely periodic structure.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14: Transverse displacement field (real part) of the plate with 8 × 4 substructures with the floor panel at 125 Hz, obtained with the interpolation strategy. (left) Periodic case; (right) nearly periodic case.</figcaption></figure>

propagation analysis in infinite systems, but also, the study of systems of finite dimensions (one or several, possibly coupled together). Follow-on works could include the analysis of periodic structures with localized nonlinear effects, e.g., metamaterials with local nonlinear resonant devices, or structures locally subject to fast loadings of high magnitudes (shocks, blast).

## REFERENCES

[1] D. Beli, A.T. Fabro, M. Ruzzene, and J.R.F. Arruda. Wave attenuation and trapping in 3D printed cantilever-in-mass metamaterials with spatially correlated variability. *Scientific Reports*, 9:5617, 2019.

[2] D. Beli, J.-M. Mencik, P.B. Silva, and J.R.F. Arruda. A projection-based model reduction strategy for the wave and vibration analysis of rotating periodic structures. *Computational Mechanics*, 62:1511–1528, 2018.
[3] M.P. Castanier and C. Pierre. Individual and interactive mechanisms for localization and dissipation in a mono-coupled nearly-periodic structure. *Journal of Sound and Vibration*, 168(3):479-505, 1993.

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

[5] D. Duhamel. Finite element computation of Green's functions. *Engineering Analysis with Boundary Elements*, 31(11):919-930, 2007.

[6] D. Duhamel, B.R. Mace, and M. J. Brennan. Finite element analysis of the vibrations of waveguides and periodic structures. *Journal of Sound and Vibration*, 294(1-2):205-220, 2006.

[7] D. Duhamel and J.-M. Mencik. Time response analysis of periodic structures via wave-based absorbing boundary conditions. *European Journal of Mechanics - A/Solids (submitted)*, 2021.

[8] G.L. Huang and C.T. Sun. Band Gaps in a Multiresonator Acoustic Metamaterial. *Journal of Vibration and Acoustics*, 132:031003, 2010.

[9] M.I. Hussein. Reduced Bloch mode expansion for periodic media band structure calculations. *Proceedings of the Royal Society A*, 465:2825-2848, 2009.

[10] M.I. Hussein, G.M. Hulbert, and R.A. Scott. Dispersive elastodynamics of 1D banded materials and structures: Design. *Journal of Sound and Vibration*, 307:865-893, 2007.

[11] J.S. Jensen. Phononic band gaps and vibrations in one- and two-dimensional mass-spring structures. *Journal of Sound and Vibration*, 266:1053-1078, 2003.

[12] D. Krattiger and M.I. Hussein. Generalized Bloch mode synthesis for accelerated calculation of elastic band structures. *Journal of Computational Physics*, 357:183-205, 2018.

[13] F.-M. Li and Y.-S. Wang. Study on wave localization in disordered periodic layered piezoelectric composite structures. *International Journal of Solids and Structures*, 42:6457-6474, 2005.

[14] B.R. Mace and E. Manconi. Modelling wave propagation in two-dimensional structures using finite element analysis. *Journal of Sound and Vibration*, 318(4-5):884-902, 2008.

[15] D.J. Mead. The forced vibration of one-dimensional multi-coupled periodic structures: An application to finite element analysis. *Journal of Sound and Vibration*, 319:282-304, 2009.

[16] J.-M. Mencik. On the low- and mid-frequency forced response of elastic systems using wave finite elements with one-dimensional propagation. *Computers and Structures*, 88(11-12):674-689, 2010.

[17] J.-M. Mencik. A wave finite element-based formulation for computing the forced response of structures involving rectangular flat shells. *International Journal for Numerical Methods in Engineering*, 95(2):91-120, 2013.
[18] J.-M. Mencik. New advances in the forced response computation of periodic structures using the wave finite element (WFE) method. *Computational Mechanics*, 54(3):789-801, 2014.

[19] J.-M. Mencik. A wave finite element approach for the analysis of periodic structures with cyclic symmetry in dynamic substructuring. *Journal of Sound and Vibration*, 431:441-457, 2018.

[20] J.-M. Mencik. Model reduction based on matrix interpolation and distorted finite element meshes for dynamic analysis of 2D nearly periodic structures. *Finite Elements in Analysis and Design*, 188:103518, 2021.

[21] J.-M. Mencik, V. Denis, M. Caliez, and S. Méo. Application of the wave finite element method to the analysis of locally resonant band gaps in periodic structures with viscoelastic properties. *Proceedings of the 28th International Conference on Noise and Vibration engineering (ISMA2018)*, Leuven, Belgium, pages 2469-2481, 2018.

[22] J.-M. Mencik and D. Duhamel. A wave-based model reduction technique for the description of the dynamic behavior of periodic structures involving arbitrary-shaped substructures and large-sized finite element models. *Finite Elements in Analysis and Design*, 101:1-14, 2015.

[23] J.-M. Mencik and M. N. Ichchou. Multi-mode propagation and diffusion in structures through finite elements. *European Journal of Mechanics - A/Solids*, 24(5):877-898, 2005.

[24] J.-M. Mencik, M. N. Ichchou, and L. Jézéquel. Propagation multimodale dans les systèmes périodiques couplés (in French). *Revue Européenne de Mécanique Numérique*, 15(1-3):293-306, 2006.

[25] 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-101:29-40, 2016.

[26] H. Panzer, J. Mohring, R. Eid, and B. Lohmann. Parametric model order reduction by matrix interpolation. *at-Automatisierungstechnik*, 58(8):475-484, 2010.

[27] P.-O. Persson and G. Strang. A Simple Mesh Generator in MATLAB. *SIAM Rev.*, 46(2):329-345, 2004.

[28] A.S. Phani, J. Woodhouse, and N. Fleck. Wave propagation in two-dimensional periodic lattices. *The Journal of the Acoustical Society of America*, 119(4):1995-2005, 2006.

[29] V.F. Dal Poggeto and A.L. Serpa. Flexural wave band gaps in a ternary periodic metamaterial plate using the plane wave expansion method. *Journal of Sound and Vibration*, 495:115909, 2021.

[30] P.B. Silva, J.-M. Mencik, and J.R.F. Arruda. Wave finite element-based superelements for forced response analysis of coupled systems via dynamic substructuring. *International Journal for Numerical Methods in Engineering*, 107(6):453-476, 2016.
[31] C. Sugino, Y. Xia, S. Leadenham, M. Ruzzene, and A. Erturk. A general theory for bandgap estimation in locally resonant metastructures. *Journal of Sound and Vibration*, 406:104–123, 2017.

[32] Y. Waki, B.R. Mace, and M.J. Brennan. Numerical issues concerning the wave and finite element method for free and forced vibrations of waveguides. *Journal of Sound and Vibration*, 327(1-2):92–108, 2009.

[33] Z. Wang, P. Zhang, and Y. Zhang. Locally Resonant Band Gaps in Flexural Vibrations of a Timoshenko Beam with Periodically Attached Multiscillators. *Mathematical Problems in Engineering*, Article ID 146975:104110, 2013.

[34] K. Willcox and J. Peraire. Balanced model reduction via the proper orthogonal decomposition. *AIAA Journal*, 40(11):2323–2330, 2002.

[35] Y. Xiao, B.R. Mace, J. Wen, and X. Wen. Formation and coupling of band gaps in a locally resonant elastic system comprising a string with attached resonators. *Physics Letters A*, 375:1485–1491, 2011.

[36] W. X. Zhong and F. W. Williams. On the direct solution of wave propagation for repetitive structures. *Journal of Sound and Vibration*, 181(3):485–501, 1995.