# Floquet-Bloch decomposition for the computation of dispersion of two-dimensional periodic, damped mechanical systems

Manuel Collet, Morvan Ouisse, Massimo Ruzzene, Mohamed Ichchou

► To cite this version:

Manuel Collet, Morvan Ouisse, Massimo Ruzzene, Mohamed Ichchou. Floquet-Bloch decomposition for the computation of dispersion of two-dimensional periodic, damped mechanical systems. International Journal of Solids and Structures, 2011, 48, pp.2837-2848. 10.1016/j.ijsolstr.2011.06.002. hal-01516397

HAL Id: hal-01516397

https://hal.science/hal-01516397

Submitted on 30 Apr 2017

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.

Public Domain
# Floquet-Bloch decomposition for the computation of dispersion of two-dimensional periodic, damped mechanical systems

M. Collet<sup>a</sup>, M. Ouisse<sup>a</sup>, M. Ruzzene<sup>b</sup>, M.N. Ichchou<sup>c</sup>

<sup>a</sup> FEMTO-ST, Applied Mechanics Department, UMR-CNRS 6174, 24 chemin de l'épitaphe, 25000 Besançon, France

<sup>b</sup> D. Guggenheim School of Aerospace Engineering, Georgia Institute of Technology, Atlanta, GA 30332, USA

<sup>c</sup> LTDS UMR-CNRS 5513, École Centrale de Lyon, 36 avenue Guy de Collongue, 69134 Écully, France

Floquet-Bloch theorem is widely applied for computing the dispersion properties of periodic structures, and for estimating their wave modes and group velocities. The theorem allows reducing computational costs through modeling of a representative cell, while providing a rigorous and well-posed spectral problem representing wave dispersion in undamped media. Most studies employ the Floquet-Bloch approach for the analysis of undamped systems, or for systems with simple damping models such as viscous or proportional damping. In this paper, an alternative formulation is proposed whereby wave heading and frequency are used to scan the *k*-space and estimate the dispersion properties. The considered approach lends itself to the analysis of periodic structures with complex damping configurations, resulting for example from active control schemes, the presence of damping materials, or the use of shunted piezoelectric patches. Examples on waveguides with various levels of damping illustrate the performance and the characteristics of the proposed approach, and provide insights into the properties of the obtained eigensolutions.

Keywords:

Waves

Bloch theorem

Damping

Mid-frequency

## 1. Introduction

Tailoring the dynamical behavior of one or two-dimensional waveguides can provide efficient and physically elegant means to optimize mechanical structures with regards to vibration and acoustic criteria, among others. However, achieving this objective may lead to different outcomes depending on the context of the optimization. In the preliminary stages of a product's development, one mainly needs optimization tools capable of rapidly providing global design directions. Such optimization will also depend on the frequency range of interest. One usually discriminates between the low frequency (LF) range and the medium frequency (MF) range, especially if vibration and noise are considered. However, it should be noted that LF optimization of vibration is more common in the literature than MF optimization. For example, piezoelectric materials and other adaptive and smart systems are employed to improve the vibroacoustic quality of structural com-

ponents, especially in the LF range (see Preumont, 1997; Nelson, 1992 or Banks and RCSmith, 1996 among many others).

Recently, much effort has been spent on developing new multi-functional structures integrating electro-mechanical systems in order to optimize their vibroacoustic behavior over a larger frequency band of interest, among which Thorp et al. (2001) or Collet et al. (2009). However, there is still a lack of studies in the literature for MF optimization of structural vibration. To that end, the aim of this study is to provide a suitable numerical tool for computing wave dispersion in two-dimensional periodic systems incorporating damping and/or active devices (visco-, poro-elastic materials, controlling electronics devices, etc). The final aim is to allow their optimization in terms of vibroacoustic diffusion in two-dimensional waveguides.

The two most popular numerical approaches that can be distinguished for computing dispersion are the Semi-Analytical Finite Element method (SAFE) and the wave finite element (WFE) method. In the former approach, the displacement field is modeled exactly in the direction of wave propagation by using a harmonic function and approximately in the directions perpendicular by using finite elements. An eigenvalue problem is then formulated by introducing the displacement field into the governing equations. Solving the eigenvalue problem for a given frequency gives the wave numbers of all the propagating modes. The main
disadvantage of the SAFE method is that FE used are not standard so they must be specifically defined for each application. Nevertheless, many specific finite elements have been developed since 1975, among which those dedicated to computation of different types of structures. (Gavrić, 1995) used this technique to calculate the dispersion relationship in a free rail by using triangular and quadrilateral elements, with those elements obtained from Hamilton's principle. (Hayashi et al., 2003) derived the SAFE formulation for the waveguides of complex cross-sectional shape through virtual work principles. By using Lagrange's equations (Damljanovic and Weaver, 2004) developed the linear triangular elements for SAFE method to investigate the elastic waves in waveguides of arbitrary cross-section. The wave mode characteristics in damped waveguides were examined by Bartoli et al. (2006), with the elements derived also from Hamilton's principle. The SAFE technique has also been adopted to investigate the wave propagation characteristics for thin-walled structures in Finnveden (2004). In order to avoid development of specific FE, the WFE method considers the structures as periodic in order to model, with standard FE, a period of the structure. By using the periodic structure theory (PST) introduced by Mead (1996), an eigenvalue problem can be formulated from the stiffness and mass matrices of the FE model to find wave numbers of all the propagating waves. Contrary to SAFE method, the displacement field is approximated in the direction of propagation. Thus, some numerical issues can arise when the size of FE are too coarse. As recommended by Mace and Manconi (2008), a minimum of six elements per wavelength is a good rule of thumb to ensure a reliable analysis. The WFE method has been successfully used to deal with wave propagation in two dimensional structures by Manconi (2008), Ichchou et al. (2008) or Akrout (2005).

One of the main problems of these two approaches is the difficulty to deal with complex mechanical wave propagation specifically of multi-modal nature. Indeed, the existence at each frequency step of a number of wavemodes that potentially exchange energy make the computation and characterization of wave attenuation a delicate task. This task is of fundamental interest in the optimization process of energy dissipation features. This paper will then concentrates in the computation of the damped wave numbers in the whole Brillouin domain defined in Brillouin (1953) that are necessary for vibroacoustic optimization behavior of such periodic structures. By using WFE technique, one can introduce structural or viscous damping, as indicated in Duhamel (2007) or Mace and Manconi (2008). Such computations involve resolution of complex polynomial or transcendental eigenvalue problem as underlined in Mace and Manconi (2008). The treatment of the obtained damped wave numbers also necessitates specific tools for defining and estimating the wave loss factors as in Manconi and Mace (2010).

In this paper, after recalling the Floquet and the Bloch theorems, a new numerical formulation is introduced for computing the multi-modal damped wave numbers in the whole first Brillouin domain of periodical structures with non homogeneous and generic frequency dependent damping terms. Then a bi-dimensional numerical application is presented in order to validate the method and to use it for estimating the bi-dimensional band gaps as well as a suitable evanescence's indicator in the context of strongly damped systems. The validated methodology can also be used for optimizing damping layers or active/semi-active elements to control vibroacoustic power flow into mechanical systems. The paper is precisely structured as follows. Section 2 reminds The Floquet-Bloch theorem for elasto-dynamical system and offers the mathematical and physical context of its application. Section 3 deals with a two-dimensional application of Section 2 main finding. The considered example corresponds to a periodic distributed passive system. Section 4 concludes the paper.

## 2. Floquet and Bloch theorems for elasto-dynamic dispersion analysis

This section summarizes Floquet and Bloch theorems and their application to elasto-dynamics. The well-known formulations of Floquet (1883) and Bloch (1928) respectively for one dimensional (1D) and two dimensional (2D) systems governed by differential equations with periodic coefficients are here specifically revisited in light of their application to the analysis of damped periodic mechanical systems.

### 2.1. The Floquet theorem

The Floquet theory is a methodology to solve ordinary differential equations of the form:

$$ \frac{d\mathbf{w}(x)}{dx} = \mathbf{A}(x)\mathbf{w}(x), \quad \forall x \in \mathbb{R}, \qquad (1) $$

where $\mathbf{w}(x): \mathbb{R} \to \mathbb{C}^n$ is the unknown function, and $\mathbf{A}(x)$ is a given matrix of continuous periodic functions with period $r_1$, i.e. $\mathbf{A}(x+r_1) = \mathbf{A}(x)$. Floquet Theorem dictates that any solution of this system of equations can be expressed as a linear combination of functions $\mathbf{v}(x)e^{ikx}$, where $\mathbf{v}(x)$ is $r_1$-periodic, while $k \in \mathbb{C}$ is a scalar complex quantity. The theory provides a way to evaluate $\mathbf{v}$ and $k$ from the solution of an eigenvalue problem.

Among the many mathematical aspects of the theory, some points should be mentioned for proper understanding. First, for any given basis $\mathbf{W}(x) \in \mathbb{C}^{n \times n}$ of fundamental solutions of (1), a new basis $\mathbf{P}_0(x)$ of solutions normalized so that $\mathbf{P}_0(0) = \mathbf{I}_n$ can be defined:

$$ \mathbf{P}_0(x) = \mathbf{W}(x)\mathbf{W}^{-1}(0), \qquad (2) $$

where $\mathbf{I}_n$ denotes the $n \times n$ identity matrix. It is possible to search for $\mathbf{W}(x+r_1)$ from the expression:

$$ \mathbf{W}(x+r_1) = \mathbf{P}_0(x)\mathbf{W}(r_1) = \mathbf{W}(x)\mathbf{W}^{-1}(0)\mathbf{W}(r_1), \qquad (3) $$

where $\mathbf{P}_0$ is the Floquet propagator which allows the evaluation of $\mathbf{W}(x+r_1)$ from knowledge of $\mathbf{W}(r_1)$. The estimation of $\mathbf{P}_0$ is based on its diagonalization performed for $x=r_1$:

$$ \mathbf{P}_0(r_1) = \mathbf{Z}\mathbf{\Lambda}\mathbf{Z}^{-1}, \qquad (4) $$

where $\mathbf{\Lambda}$ and $\mathbf{Z}$ contain the solutions of the following eigenvalue problem:

$$ \mathbf{P}_0(r_1)\mathbf{z}_j = \lambda_j\mathbf{z}_j, \qquad (5) $$

so that $\mathbf{\Lambda}$ is a diagonal matrix with $\lambda_j$ terms, and $\mathbf{Z}$ is the matrix grouping eigenvectors $\mathbf{z}_j$ as columns. The eigenvalues can also be written as:

$$ \mathbf{\Lambda} = e^{\mathbf{K}r_1}, \qquad (6) $$

where $\mathbf{K}$ is a diagonal matrix whose generic element is $k_j$ such that $\lambda_j = e^{k_j r_1}$. The parameter $k_j$ is the jth Floquet (characteristic) exponent, while $\lambda_j$ is the corresponding Floquet multiplier. The computation of the eigenvalues is not performed directly on $\mathbf{P}_0$, since a more convenient approach identifies $\mathbf{Y}$ as the Floquet propagation of basis $\mathbf{Z}$ such that:

$$ \mathbf{Y}(x) = \mathbf{P}_0(x)\mathbf{Z}. $$

Also, It may be shown that:

$$ \mathbf{Y}(x+r_1) = \mathbf{Y}(x)e^{\mathbf{K}r_1}. \qquad (7) $$

The vectors included as columns in $\mathbf{Y}(x)$ are solutions of the initial periodic problem (1) restricted to the elementary cell $[0, r_1]$, with fixed boundary conditions at $x=0$ and $x=r_1$. Accordingly, the
eigenvectors $z_j$ and eigenvalues $\lambda_j = e^{k_j r_1}$ are solutions of the generalized eigenvalue problem:

$$ egin{cases} \frac{d\mathbf{Y}(x)}{dx} = \mathbf{A}(x)\mathbf{Y}(x) & \forall x \in [0, r_1], \\ \mathbf{Y}(0) = \mathbf{Z}, \\ \mathbf{Y}(r_1) = \mathbf{Z}\mathbf{\Lambda}. \end{cases} \quad (8) $$

The Floquet propagators are then obtained from:

$$ \mathbf{P}_0(x) = \mathbf{Y}(x)\mathbf{Z}^{-1} \quad (9) $$

and a basis of solutions is given by:

$$ \mathbf{W}(x) = \mathbf{Y}(x)\mathbf{W}(0). \quad (10) $$

An alternative way computes the Floquet propagators by considering the undamped Floquet vectors $\mathbf{v}_j$ defined from the following relation:

$$ \mathbf{V}(x) = \mathbf{Y}(x)e^{-\mathbf{K}x}, \quad (11) $$

where $\mathbf{V}$ is the matrix containing the vectors $\mathbf{v}_j$ as columns. It can be shown that these function are $r_1$-periodic, and they are solutions of the following problem:

$$ egin{cases} \frac{d\mathbf{V}(x)}{dx} = \mathbf{A}(x)\mathbf{V}(x) - \mathbf{V}(x)\mathbf{K} & \forall x \in [0, r_1], \\ \mathbf{V}(0) = \mathbf{Z}, \\ \mathbf{V}(r_1) = \mathbf{Z}. \end{cases} \quad (12) $$

The generalized eigenvalue problem (12) is equivalent to (8) and gives eigenvectors $z_j$ and eigenvalues $k_j$. The solution in terms of Floquet propagator can then be expressed as:

$$ \mathbf{P}_0(x) = \mathbf{V}(x)e^{\mathbf{K}x}\mathbf{Z}^{-1}, \quad (13) $$

while a basis of solutions is given by:

$$ \mathbf{W}(x) = \mathbf{V}(x)e^{\mathbf{K}x}\mathbf{W}(0). \quad (14) $$

Eq. (14) is the Floquet normal form of the fundamental basis $\mathbf{W}(x)$. The characteristic multipliers in Eq. (14) are also the eigenvalues of the linear Poincaré maps defined as the function $\mathbf{w}(x) \to \mathbf{w}(x + r_1)$, where $\mathbf{w}(x)$ is a solution of (1).

Based on the above discussion, it is clear that two different approaches for the calculation of wave solution are possible: the Floquet propagators can be obtained by solving the eigenvalue problem (8) or by estimating the Floquet vectors from (12). The two resulting eigensolutions are related to each other and permit the computation of wave solutions for (1). However, they are obtained from two different eigenvalue problems. The first approach computing the Floquet propagators (8) leads to compute non standard eigen solutions of a problem where the eigenvalues (i.e. the Floquet multipliers) appear in the expression of the applied boundary conditions, while the second approach considers only standard Dirichlet periodic boundary conditions and computes the eigenvalues (i.e. the Floquet exponents) by solving a standard eigenvalue problem by introducing another expression of the used operator inside the cell domain (1). The distinction between these two approaches is essential to the understanding of the numerical implementation presented in what follows. Of note and relevant to the upcoming discussion is that Floquet exponents are not unique since $e^{(k_j + i\frac{2\pi m}{r_1})r_1} = e^{k_j r_1}$ where $m$ is an integer. Also, Floquet vectors are periodic, and therefore bounded on $\mathbb{R}$. The stability of homogeneous solutions of (1) are also given by the value of the Lyapunov exponents, which are the real parts of the Floquet exponents: the solutions are asymptotically stable if all Lyapunov exponents are negative, Lyapunov stable if the Lyapunov exponents are nonpositive and unstable otherwise. These properties remain valid when multi-dimensional problems are considered.

## 2.2. The Bloch theorem

Bloch theorem was originally introduced to represent the form of homogeneous states of Schrödinger equation with periodic potential. This theorem can be considered as a multidimensional application of the Floquet theorem, as indicated by Joannopoulos et al. (1995).

For illustration purposes, we consider a medium whose generic property $M$ satisfies the periodicity condition:

$$ M(\mathbf{x} + \mathbf{Rm}) = M(\mathbf{x}), $$

where $\mathbf{m} \in \mathbb{Z}^3$, and $\mathbf{R} = [\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3] \in \mathbb{R}^{3 \times 3}$ is a matrix containing the three lattice vectors $\mathbf{r}_j, j = 1, \dots, 3$, as illustrated in Fig. 1. The primitive cell is defined as a convex polyhedron of $\mathbb{R}^3$ called $\Omega_R$. The reciprocal unit cell, denoted by $\Omega_G$ is defined by the reciprocal lattice vector basis $\mathbf{g}_k$ for which the following holds:

$$ \mathbf{r}_j \cdot \mathbf{g}_k = 2\pi\delta_{jk}, $$

where $\delta_{jk}$ the Kronecker delta. Also, $\mathbf{G} = [\mathbf{g}_1, \mathbf{g}_2, \mathbf{g}_3]$ is the reciprocal lattice matrix in the later. If $\Omega_R$ is the irreducible primitive cell, $\Omega_G$ corresponds to the first Brillouin zone of the lattice (see (Kittel, 1986) for details).

The Bloch theorem stipulates that any functions $\mathbf{u}(\mathbf{x}) \in L^2(\mathbb{R}^3, \mathbb{C}^n)$ can be expressed as

$$ \mathbf{u}(\mathbf{x}) = \int_{\Omega_G} \tilde{\mathbf{u}}(\mathbf{x}, \mathbf{k})e^{i\mathbf{k}\cdot\mathbf{x}}d\mathbf{k}, \quad (15) $$

where the Bloch amplitude $\tilde{\mathbf{u}}(\mathbf{x}, \mathbf{k})$ is $\Omega_G$-periodic in $\mathbf{k}$ and can be represented as:

$$ \begin{aligned} \tilde{\mathbf{u}}(\mathbf{x}, \mathbf{k}) &= \sum_{\mathbf{n} \in \mathbb{Z}^3} \tilde{\mathbf{u}}(\mathbf{k} + \mathbf{G}\mathbf{n})e^{i\mathbf{G}\mathbf{n}\cdot\mathbf{x}}, \\ &= \frac{|\Omega_R|}{(2\pi)^3} \sum_{\mathbf{n} \in \mathbb{Z}^3} \mathbf{u}(\mathbf{x} + \mathbf{R}\mathbf{n})e^{-i\mathbf{k}\cdot(\mathbf{x} + \mathbf{R}\mathbf{n})}, \end{aligned} \quad (16) $$

where $\tilde{\mathbf{u}}(\mathbf{k})$ stands for the Fourier transform of $\mathbf{u}(\mathbf{x})$. It can also be demonstrated that the mean value of the Bloch amplitude is the Fourier amplitude of $\mathbf{u}(\mathbf{x})$ for the corresponding wave vector:

$$ \langle \tilde{\mathbf{u}}(\mathbf{x}, \mathbf{k}) \rangle_{\Omega_x} = \hat{\mathbf{u}}(\mathbf{k}), $$

where $\langle \cdot \rangle_{\Omega_x}$ denotes for the spatial mean value computed on domain $\Omega_x$.

The application of Bloch theorem for the representation of solutions of partial differential equations with periodic coefficients

<figure><img src="image_1.png" /><figcaption>Fig. 1. Generic 3D periodic cells.</figcaption></figure>
allows for all derivatives to be shifted by **k** in the sense given by the considered spatial operator. Consider for example the second order elliptic operator:

$$A(\mathbf{x}) = - \sum_{p,q=1}^{N} \frac{\partial}{\partial x_p} \left( \mathbf{A}_{pq}(\mathbf{x}) \frac{\partial}{\partial x_q} \right) + \mathbf{W}(\mathbf{x}) \quad (17)$$

defined on smooth functions of $\mathbb{R}^N$ (i.e. $C^\infty(\mathbb{R}^N)$) where $A_{pq}(\mathbf{x})$ is a symmetric, smooth and uniformly positive matrix with $\Omega_R$-periodic coefficients. Indices $p$ and $q$ denote each basis vector of the considered $\mathbb{R}^N$ domain. Also, $\mathbf{W}(\mathbf{x}) \ge 0$ is a real smooth function $\Omega_R$-periodic. A spectral resolution of the closure of this operator can be found in $L^2(\mathbb{R}^N)$, and can be expressed in terms of Bloch waves associated with $A(\mathbf{x})$. Indeed, let $\mathbf{k} \in \Omega_G$, and define:

$$A(\mathbf{x}, \mathbf{k}) = - \sum_{p,q=1}^{N} \left( \frac{\partial}{\partial x_p} + ik_p \right) \left( a_{pq}(\mathbf{x}) \frac{\partial}{\partial x_q} + ik_q \right) + W(\mathbf{x}). \quad (18)$$

The shifted cell eigenvalue problem is then considered:

$$A(\mathbf{x}, \mathbf{k}) \mathbf{v}_j(\mathbf{x}, \mathbf{k}) = \omega_j^2 \mathbf{v}_j(\mathbf{x}, \mathbf{k}), \quad (19)$$

for $\mathbf{k} \in \Omega_k$ with $\mathbf{v}(\mathbf{x}, \mathbf{k})$ in $H^1(\Omega_R)$ and $\Omega_R$-periodic.

The essentially self-adjoint operator $A(\mathbf{x}, \mathbf{k})$ is non negative and the eigenvalue problem (19) has a discrete sequence of eigenvalues $\omega_1^2$, with $j \in \mathbb{N}$ whose corresponding eigenfunctions are the Bloch waves $\mathbf{v}_j(\mathbf{x}, \mathbf{k})$. These eigenvalues are smooth functions of $\mathbf{x}$ and are complete in $L^2(\Omega_R)$. More details in mathematical properties of this eigenvalue problem can be found in Bensoussan et al. (1978) and Wilcox (1978). Based on these results, the Bloch expansion of any function $\mathbf{u}(\mathbf{x}) \in L^2(\mathbb{R}^3, \mathbb{C}^n)$ can be expressed as:

$$\mathbf{u}(\mathbf{x}) = \int_{\Omega_G} \sum_{j=1}^{+\infty} u_j(\mathbf{k}) e^{ikx} \mathbf{v}_j(\mathbf{x}, \mathbf{k}) d\mathbf{k} \quad (20)$$

and

$$u_j(\mathbf{k}) = \int_{\mathbb{R}^N} \mathbf{u}(\mathbf{x}) e^{-ikx} \mathbf{v}_j(\mathbf{x}, \mathbf{k}) d\mathbf{x}, \quad (21)$$

where $\tilde{\mathbf{v}}_j$ is the complex conjugate of $\mathbf{v}_j$. Moreover, Parseval's identity holds

$$\int_{\mathbb{R}^N} |\mathbf{u}(\mathbf{x})|^2 d\mathbf{x} = \int_{\Omega_G} \sum_{j=1}^{+\infty} |u_j(\mathbf{k})|^2 d\mathbf{k}. \quad (22)$$

The spectral resolution of operator $A$ can also be expressed as:

$$A(\mathbf{x})\mathbf{u}(\mathbf{x}) = \int_{\Omega_G} \sum_{j=1}^{+\infty} u_j(\mathbf{k}) e^{ikx} \omega_j^2(\mathbf{k}) \mathbf{v}_j(\mathbf{x}, \mathbf{k}) d\mathbf{k}. \quad (23)$$

Examples of applications of these results can be found in Bensoussan et al. (1978).

## 2.3. Application to elastodynamic

Let us consider an infinite periodic elastodynamic problem as presented in Fig. 1. The harmonic homogeneous dynamical equilibrium of system is driven by the following partial derivative equation:

$$\rho(\mathbf{x})\omega^2\mathbf{w}(\mathbf{x}) + \nabla\mathbf{C}(\mathbf{x})\nabla_{sym}(\mathbf{w}(\mathbf{x})) = 0 \quad \forall \mathbf{x} \in \mathbb{R}^3, \quad (24)$$

where $\mathbf{w}(\mathbf{x}) \in \mathbb{R}^3$ is the displacement vector, $\mathbf{C}(\mathbf{x})$ stands for the Hook elasticity tensor and $\varepsilon(\mathbf{x}) = \nabla_{sym}(\mathbf{w}(\mathbf{x})) = \frac{1}{2}(\nabla\mathbf{w}^T(\mathbf{x}) + \mathbf{w}(\mathbf{x})\nabla^T)$ is the strain tensor. By considering a primitive cell of the periodic problem $\Omega_R$ and by using the Bloch theorem, the associated Bloch eigenmodes (19) and the dispersion functions can be found by searching the eigen solutions of the homogeneous problem (24) as:

$$\mathbf{w}(\mathbf{x}) = \mathbf{w}_{n,k}(\mathbf{x}, \mathbf{k})e^{ikx}, \quad (25)$$

where $\mathbf{w}_{n,k}(\mathbf{x}, \mathbf{k})$ are $\Omega_R$-periodic functions. In that case $\mathbf{w}_{n,k}(\mathbf{x}, \mathbf{k})$ and $\omega_n(\mathbf{k})$ are the solutions of the following generalized eigenvalues problem:

$$\rho(\mathbf{x})\omega_n(\mathbf{k})^2 \mathbf{w}_{n,k}(\mathbf{x}) + \nabla\mathbf{C}(\mathbf{x})\nabla_{sym}(\mathbf{w}_{n,k}(\mathbf{x})) - i\mathbf{C}(\mathbf{x})\nabla_{sym}(\mathbf{w}_{n,k}(\mathbf{x})) \cdot \mathbf{k} \\ - i\nabla\mathbf{C}(\mathbf{x}) \frac{1}{2}(\mathbf{w}_{n,k}(\mathbf{x}) \cdot \mathbf{k}^T + \mathbf{k} \cdot \mathbf{w}_{n,k}^T(\mathbf{x})) + \mathbf{C}(\mathbf{x}) \frac{1}{2}(\mathbf{w}_{n,k}(\mathbf{x}) \cdot \mathbf{k}^T + \mathbf{k} \cdot \mathbf{w}_{n,k}^T(\mathbf{x})) \cdot \mathbf{k} = 0 \quad \forall \mathbf{x} \in \Omega_R, \quad (26)$$

$$\mathbf{w}_{n,k}(\mathbf{x} - \mathbf{R}\mathbf{n}) - \mathbf{w}_{n,k}(\mathbf{x}) = 0 \quad \forall \mathbf{x} \in \Gamma_R. \quad (27)$$

The first equation is simply obtained by introducing Eq. (25) into elastodynamic Eq. (24). The second equation represents the symmetrical boundary conditions expressed on boundary faces of the lattice polyhedron as described in Fig. 2(a) for a rectangular parallelepiped cell. In this equation $\mathbf{n}$ stands for the outpointing unitary normal vector. It corresponds to a complex Quadratic Eigenvalue Problem (QEP) that can be solved by fixing two of the constants $\omega$, $|\mathbf{k}|$ (the complex amplitude) or cosine directions of $\mathbf{k}$ and compute the last one.

The proposed formulation is based on the computation of the Floquet vectors from Eq. (26), instead of computing the Floquet propagators commonly used for elastodynamic applications. The methodology allows the computation of the full complex map of

<figure><img src="image_2.png" /><figcaption>Fig. 2. (a) Rectangular parallelepiped primitive lattice (b) Corresponding rectangular parallelepiped reciprocal lattice.</figcaption></figure>
the dispersion curves incorporating computation of evanescent waves and allowing the introduction of damping operator if any.

### 2.3.1. Weak formulation and computation of waves dispersion curves in periodical lattice

Let us consider the partial derivative Eq. (26) on a unit cell $\Omega_R$. It stands for a generalized eigenvalue problem leading to computation of the dispersion curves $\omega_n(k)$ and corresponding Floquet eigenvectors $w_{n,k}(x)$.

For 3D applications, the wave vectors are supposed to be complex if damping terms are added into Eq. (26), they can be written

as $k = k \begin{bmatrix} \sin(\theta) \cos(\phi) \\ \sin(\theta) \sin(\phi) \\ \cos(\theta) \end{bmatrix}$ where $\theta$ and $\phi$ represent the direction angles into the reciprocal lattice domain as shown in Fig. 2(b). This decomposition assumes that real and imaginary parts of vector $k$

are co-linear. In the following, $\Phi = \begin{bmatrix} \sin(\theta) \cos(\phi) \\ \sin(\theta) \sin(\phi) \\ \cos(\theta) \end{bmatrix}$ indicates that direction vector.

### 2.3.2. Weak formulation

If $w_{n,k}(x)$ is a solution of Eq. (26), then:

$$ \begin{aligned} \forall \tilde{w}_{n,k}(\mathbf{x}) \in \{H_1(\Omega_R, \mathbb{C}^3) / \tilde{w}_{n,k}(\mathbf{x} - R\mathbf{n}) = \tilde{w}_{n,k}(\mathbf{x}) \forall \mathbf{x} \in \Gamma_R\}, \\ \times \int_{\Omega_R} \rho(\mathbf{x}) \omega_n^2(\mathbf{k}) \tilde{w}_{n,k}(\mathbf{x}) w_{n,k}(\mathbf{x}) - \tilde{\varepsilon}_{n,k}(\mathbf{x}) C(\mathbf{x}) \varepsilon_{n,k}(\mathbf{x}) \\ + ik \tilde{\kappa}_{n,k}(\mathbf{x}) C(\mathbf{x}) \varepsilon_{n,k}(\mathbf{x}) - ik \tilde{\varepsilon}_{n,k}(\mathbf{x}) C(\mathbf{x}) \kappa_{n,k}(\mathbf{x}) \\ + k^2 \tilde{\kappa}_{n,k}(\mathbf{x}) C(\mathbf{x}) \kappa_{n,k}(\mathbf{x}) d\Omega + \int_{\Gamma_R} \tilde{w}_{n,k}(\mathbf{x}) (C(\mathbf{x}) (\varepsilon_{n,k}(\mathbf{x}) \\ + ik \kappa_{n,k}(\mathbf{x}))) \cdot n d\Gamma = 0, \end{aligned} \quad (28) $$

where $\varepsilon_{n,k}(\mathbf{x}) = \nabla_{\text{sym}}(w_{n,k}(\mathbf{x}))$ is the strain tensor, $\kappa_{n,k}(\mathbf{x}) = \frac{1}{2}(w_{n,k}(\mathbf{x})\Phi^T + \Phi w_{n,k}^T(\mathbf{x}))$ is the symmetric dyadic tensor or the dyadic product of the displacement $w_{n,k}(\mathbf{x})$ and the direction vector $\Phi$. ~means that the specified operator is applied to the test functions and $n$ is the unit outpointing normal vector on the considered boundary.

This weak formulation is simply obtained by integrating Eq. (26) projected onto any test function $\tilde{w}_{n,k}(x)$. The boundary integral vanishes as the test functions are chosen so that $\tilde{w}_{n,k}(x - Rn) = w_{n,k}(x)$ on $\Gamma_R$ that implies $C(x)(\varepsilon_{n,k}(x - Rn) + ik\kappa_{n,k}(x - Rn)) \cdot n(x - Rn) = -C(x)(\varepsilon_{n,k}(x) + ik\kappa_{n,k}(x)) \cdot n(x)$. That corresponds to the exact compensation of the boundary applied generalized constraints $C(x)(\varepsilon_{n,k}(x) + ik\kappa_{n,k}(x))$. For a polyhedron cell, each boundary is a polyhedral plane sub-domain that can be associated with a parallel opposite one. The symmetry conditions $w_{n,k}(x - Rn) = w_{n,k}(x)$ explicitly link these associated surfaces. As the corresponding normal vector $n$ are opposite, $\kappa_{n,k}(x - Rn) = \kappa_{n,k}(x)$ and the stress condition can be restricted to $C(x)(\varepsilon_{n,k}(x - Rn)) \cdot n(x - Rn) = -C(x)(\varepsilon_{n,k}(x)) \cdot n(x)$ on the two opposite surfaces. Thus, all boundary integrations vanish and the weak formulation can be written as:

$$ \begin{aligned} \forall \tilde{w}_{n,k}(\mathbf{x}) \in \{H_1(\Omega_R, \mathbb{C}^3) / \tilde{w}_{n,k}(\mathbf{x} - R\mathbf{n}) = \tilde{w}_{n,k}(\mathbf{x}) \forall \mathbf{x} \in \Gamma_R\}, \\ \times \int_{\Omega_R} \rho(\mathbf{x}) \omega_n^2(\mathbf{k}) \tilde{w}_{n,k}(\mathbf{x}) w_{n,k}(\mathbf{x}) - \tilde{\varepsilon}_{n,k}(\mathbf{x}) C(\mathbf{x}) \varepsilon_{n,k}(\mathbf{x}) \\ + ik \tilde{\kappa}_{n,k}(\mathbf{x}) C(\mathbf{x}) \varepsilon_{n,k}(\mathbf{x}) - ik \tilde{\varepsilon}_{n,k}(\mathbf{x}) C(\mathbf{x}) \kappa_{n,k}(\mathbf{x}) \\ + k^2 \tilde{\kappa}_{n,k}(\mathbf{x}) C(\mathbf{x}) \kappa_{n,k}(\mathbf{x}) d\Omega = 0. \end{aligned} \quad (29) $$

### 2.3.3. Numerical computation

The numerical implementation is obtained by using a standard finite elements method to discretize the weak formulation (29). The assembled matrix equation is given by:

$$ (K + \lambda L(\Phi) - \lambda^2 H(\Phi) - \omega_n^2(\lambda, \Phi) M) w_{n,k}(\Phi) = 0, \quad (30) $$

where $\lambda = ik$, $M$ and $K$ are respectively the standard symmetric definite mass and symmetric semi-definite stiffness matrices, $L$ is a skew-symmetric matrix and $H$ is a symmetric semi-definite positive matrix:

$$ \begin{aligned} M &\rightarrow \int_{\Omega_R} \rho(\mathbf{x}) \omega_n^2(\mathbf{k}) \tilde{w}_{n,k}(\mathbf{x}) w_{n,k}(\mathbf{x}) d\Omega, \\ K &\rightarrow \int_{\Omega_R} \tilde{\varepsilon}_{n,k}(\mathbf{x}) C(\mathbf{x}) \varepsilon_{n,k}(\mathbf{x}) d\Omega, \\ L &\rightarrow \int_{\Omega_R} -\tilde{\kappa}_{n,k}(\mathbf{x}) C(\mathbf{x}) \varepsilon_{n,k}(\mathbf{x}) + \tilde{\varepsilon}_{n,k}(\mathbf{x}) C(\mathbf{x}) \kappa_{n,k}(\mathbf{x}) d\Omega, \\ H &\rightarrow \int_{\Omega_R} \tilde{\kappa}_{n,k}(\mathbf{x}) C(\mathbf{x}) \kappa_{n,k}(\mathbf{x}) d\Omega. \end{aligned} \quad (31) $$

When $k$ and $\Phi$ are fixed, the system (30) is a linear eigen value problem allowing us to compute the dispersion functions $\omega_n^2(\mathbf{k}, \Phi)$ and the associated Bloch eigenvector $w_{n,k}(\Phi)$.

This approach has been widely used for developing homogenization techniques and spectral asymptotic analyses like in the work of Allaire and Congas (1998). It can also be applied for computing wave's dispersion even if Floquet propagators are preferred for 1D or quasi 1D computation, as indicated by Ichchou et al. (2007), Houillon et al. (2005) or Mencik and Ichchou (2005). Nevertheless these approaches have been only developed for undamped mechanical systems that is to say represented by a set of real matrices. In this case, most of the previously published works present techniques based on the mesh of a real $k$-space (i.e $k$ or $\lambda$ and $\Phi$) inside the first Brillouin zone for obtaining the corresponding frequency dispersion diagrams and the associated Floquet vectors. For undamped systems, only propagative or evanescent waves exist, corresponding to families of eigen solutions purely real or imaginary. Discrimination between each class of waves is easy. If a damped system is considered, that is to say if matrices $K, L, H$ are complex, evanescent part of propagating waves appear as the imaginary part of $\omega_n^2(\lambda, \Phi)$ and vice versa. It then becomes very difficult to distinguish the two families of waves but also to compute the corresponding physical wave's movements by applying spatial deconvolution.

Another possibility much more suitable for computing damped system, dedicated for time/space deconvolution and for computation of diffusion properties as defined by Collet et al. (2009) or Mencik and Ichchou (2005), is to consider the following generalized eigen value problem:

$$ (K - \omega^2 M) + \lambda_n(\omega, \Phi) L(\Phi) - \lambda_n^2(\omega, \Phi) H(\Phi) w_{n,k}(\Phi) = 0. \quad (32) $$

In this problem, the pulsation $\omega$ and the propagative angle $\Phi$ are fixed real parameters. Wave's numbers $\lambda_n = ik_n$ and associated Floquet vectors $w_{n,k}$ are then computed by solving the quadratic eigen problem. This approach allows introduction of frequency dependent matrices corresponding to generalized damping terms (viscoelasticity), multiphysic coupling (especially electromechanical with electronic ordinary differential equation), foam (Biot-Allard model) or open domain boundary conditions (Sommerfeld condition).

Based on this approach, an inverse Fourier transformation in the $k$-space domain can lead us to evaluate the physical wave's displacements and energy diffusion operator when the periodic distribution is connected to another system, like in the work by Collet et al. (2009). Another temporal inverse Fourier transformation can furnish a way to access spatio-temporal response for non-homogeneous initial conditions. As $L$ is skew-symmetric, the obtained eigen values are quadruple $(\lambda, \bar{\lambda}, -\lambda, -\bar{\lambda})$ collapsing into real or imaginary pairs (or a single zero) when all matrices are real (i.e. for an undamped system). In this case a real pair of eigen values correspond to evanescent modes oriented in two opposite directions on the $k$-space and imaginary values to two traveling waves propagating in opposite direction. The obtained eigen solutions are similar, in 1D for homogeneous material (non periodic), to those given by SAFE method and additional important properties can be extrapolated from Gavrić (1994).
<figure><img src="image_3.png" /><figcaption>Fig. 3. (a) Schematic of the considered bi-dimensional waveguide (b) Description of the unit cell with periodic stubbed surface.</figcaption></figure>

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

<figure><img src="image_5.png" /><figcaption>Fig. 4. Bi-dimensional physical cell and corresponding first Brillouin zone. Shaded area is its irreducible part.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Unrefined (a) and refined (b) Mesh cases.</figcaption></figure>

As previously mentioned, the real part of $k = k\Phi$ vector is re-
stricted to stand inside the first Brillouin zone (see Fig. 2(b)). In
the quadratic eigen value problem (2.3.3) nothing restricts compu-
tation to only find eigen values satisfying this condition. For direc-
tion vector $\Phi$ orthogonal to the lattice facelets (i.e. for $\Phi_{p1} = [1,0]^T$ and $\Phi_{p2} = [0,1]^T$ in bi-dimensional rectangular cell), the periodical
conditions expressed for one dimensional waveguide are still valid:
if $\lambda_j(\Phi_p)$ is an eigen value associated to $w_{j,k}(\Phi_p)$ then $\forall m \in \mathbb{Z}^3$,
$\lambda + i \cdot \Phi_p^T(G \cdot m)$ is also an eigen value associated to $w_{j,k}(\Phi_p)$
$e^{-i \cdot \Phi_p^T(G \cdot m)x}$. Thus, for undamped systems, all obtained eigenvalues
are periodically distributed in the k-space along its principal
directions.

**3. Applications for computing bi-dimensional waves dispersion**

Illustrations in this paper are limited to bi-dimensional wave-
guide applications. Thus, it can easily be found in literature
comparative works to validate this new computational methodol-
ogy. Two different systems are considered in this section. The first

<figure><img src="image_7.png" /><figcaption>Fig. 6. Dispersion curves for undamped system (imaginary part of $\lambda_n(\omega)$). Plain lines: standard method, dot: unrefined modeling of the proposed procedure, cross: refined modeling of the proposed procedure.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Evanescence ratio $Ind(\omega)$. Plain lines: refined mesh, dashed line: unrefined mesh.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 8. Whole propagative wave numbers ($imag(\lambda_n(\omega))$) when $\phi = 0$.</figcaption></figure>

one corresponds to the undamped system used by Wu et al. (2009) to validate our computation and the second one corresponds to the damped version of the same system.

## 3.1. Undamped wave dispersion and band-gap computation in thin plate with periodic stubbed surfaces

The system is presented in the work of Wu et al. (2009). It consists of an infinite periodic bi-dimensional waveguide shown in Fig. 3(a). The system is made of a 1 mm thick aluminum plate with periodic cylindrical stubs on one of its faces as shown in Fig. 3(b). The whole system is made of isotropic Aluminum 6063-T83 ($v = 0.33, E = 69e9[Pa], \rho = 2700[kg/m^3]$).

By using symmetry of the unit cell, the corresponding first Brillouin zone is described in Fig. 4 where the irreducible zone is the shaded area. The method allows us to compute eigen frequencies corresponding to any $k$ vector described in cylindric coordinates system by its radius $k$ and its angle $\phi$.

### 3.1.1. Numerical implementation

The numerical implementation is based on the 3D weak formulation (29), using a bi-dimensional orientation in the $k$-space by imposing $\Phi = \begin{bmatrix} \cos(\phi) \\ \sin(\phi) \\ 0 \end{bmatrix}$. The applied boundary conditions are equalities of all 3D displacements on the two pairs of lateral faces

<figure><img src="image_10.png" /><figcaption>Fig. 9. Propagative wave numbers of damped system ($imag(\lambda_n(\omega))$) when $\phi = 0$ (left) and $\phi = 18^\circ$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 10. Directivity of damped system using evanescence indicator (33) saturated at unit value in the case of 1% of damping ratio.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Directivity of damped system using evanescence indicator (33) saturated at unit value in the case of 10% of damping ratio.</figcaption></figure>
$\Gamma_{r_1^+}$, $\Gamma_{r_1}$ and $\Gamma_{r_2^+}$, $\Gamma_{r_2}$ as depicted in Fig. 3(b). To impose such Dirichlet boundary conditions an extrusion coupling variable maps of displacements from the source face is exported to the destination face corresponding to the opposite one (i.e. from $\Gamma_{r_1^+}$ to $\Gamma_{r_1}$). As the domains are of the same space dimension, we typically use a point-wise mapping. The exported mapping is also coupled to the destination displacement by using dedicated Lagrange multipliers. The implementation is made with COMSOL Multiphysics® platform and parametric computation to obtain $k(\omega, \phi)$ is carried out with Matlab® routines.

For each parameters $\omega$ and $\phi$, the quadratic eigenvalue problem can be reformulated as a first order one by doubling the state dimension. After constraint handling, it is possible to write the system in the form $Ax = \lambda Bx$. The algorithm computes the largest eigenvalues of the matrix $C = A^{-1}B$. To do this, the solver uses the ARPACK FORTRAN routines for large-scale eigenvalue problems which is described by Lehoucq et al. (1998). This code is based on a variant of the Arnoldi algorithm: the implicitly restarted Arnoldi method (IRAM). The ARPACK routines must perform several matrix-vector multiplications $Cv$, which are accomplished here by solving the linear system $Ax = Bv$ using the PARDISO solver developed by Schenk and Gärtner (2004). This procedure uses double precision floating point numbers and is implemented using out of core memory management in order to avoid any memory problem even when dense (and converged) mesh is considered as shown in the following.

For all presented examples, computations have been carried out with $\omega = 2\pi \cdot [1000:1000:200000]$ (frequency between 1 and 200 kHz) and $\phi = [0: \frac{\pi}{20} : \frac{\pi}{2}]$.

The mesh of the cell is shown in Fig. 5. The first mesh case consists of 296 tetrahedral Lagrange quadratic elements for 1947 degrees of freedom and the refined one of 1550 tetrahedral Lagrange quadratic elements for 23913 degrees of freedom.

### 3.1.2. Dispersion along $\Gamma - X$ direction of the undamped system

A first computation has been made to compare our results with those presented by Wu et al. (2009). The proposed method is applied for computing the wave's dispersion curves of the undamped system along the $\Gamma - X$ direction (i.e for $\phi = 0$, see Fig. 4). Fig. 6 shows three different computations of the same dispersion curves. The first one (plain red line) corresponds to the direct simulation of the undamped system by fixing $k$ along the $\Gamma - X$ segment in the Brillouin zone and computing the corresponding eigenfrequencies $\omega$ by using a standard numerical method based on Eq. (26) (Aberg and Gudmundson, 1997; Mace and Manconi, 2008). The second and third dispersion curves (in dotted and crossed lines in Fig. 6) show, respectively, the results obtained with the unrefined and refined meshes cases. The results show a really good agreement between the standard computation method used to obtain the reference results as proposed by Wu et al. (2009) and the proposed method with the refined mesh. It can be pointed out that the evanescent modes are included into the computation and are represented by crosses points with a null imaginary parts located along the frequency axis. It also shows the convergence of the refined model compared to the unrefined one. The running time for the refined mesh case is 13.897 s and 0.621 s for the unrefined one by using an Intel Core i7 CPU running at 2.67 GHz with a RAM of 8 Go. The convergence of the proposed method is also obtained

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Directivity of damped system using wave numbers in the first Brillouin zone in the case of 1% of damping ratio: each mark indicates a computed solution, the size of the mark is a measure of the propagative nature of the wave ("evanescent" waves correspond to small radius, "propagative" waves correspond to large radius).</figcaption></figure>
by using 12 time more degrees of freedom for 22 time more computing time which 3 more ARPACK iterations. These computations validate the numerical implementation of the proposed method.

For evaluating the band-gap of the periodic system, an indicator of minimal evanescence ratio of all the computed waves for each considered frequency can be used, defined as:

$$Ind(\omega, \phi) = \min_n \frac{|\text{Real}(\lambda_n)|}{|(\lambda_n)|}. \qquad (33)$$

Fig. 7 shows the plot of this indicator for both mesh cases. The location of the first two stop bands of the system can be observed: the first one is from 40 to 50 kHz and the second from 156 to 176 kHz. Precision of these results depends on the frequency discretization rate. The obtained band-gap is totally comparable with those computed by Wu et al. (2009). Fig. 7 allows us to observe convergence of the obtained results when refined mesh density is improved.

As previously mentioned, the Bloch theorem only allows computation of waves dispersion into the first irreducible Brillouin zone, here for $k \cos(\phi)$ and $k \sin(\phi)$ inside the shaded area in Fig. 4. The obtained wavenumbers are symmetrical according to the boundary conditions of the corresponding polyhedral surface. This property is observed in Fig. 8 where the whole set of obtained wavenumbers (i.e the imaginary parts of $\lambda_n(\omega)$) is plotted. It can be observed that they are symmetrical with respect to the vertical axes on $\pm \frac{\pi}{r_1} = 100\pi$ when $\phi = 0$.

## 3.2. Dispersion of the damped system in the whole bi-dimensional k-space

The proposed computational method allows us to compute multi modal wave's propagations in the complete bi-dimensional

k-space in the first Brillouin zone. The proposed methodology is based on the computation of complex wave numbers as a function of frequency. The Bloch theorem is expended in the case of damped systems and the obtained values become complex integrating phase velocity and evanescent part for each computed wave number associated to the real and imaginary parts of the obtained eigenvalues of (2.3.3). The damping behavior is introduced by assuming a complex Hook elasticity tensor. The same methodology could have been realized by introducing any kind of linear visco-elastic modeling such as viscous behavior or any other complex frequency-dependent terms.

The first calculations have been done considering a 1% damping ratio on the same structure as the one presented before. A specific procedure has been developed to track the waves from one frequency to another in the dispersion curve, in order to follow the characteristic propagating waves: starting from a set of waves which are considered as propagative (typically such that the ratio of the real part of eigenvalue to its modulus is lower than 5%), a MAC-based correlation criteria is computed to associate the waves from one frequency step to another:

$$MAC(\mathbf{u}, \mathbf{v}) = \frac{|\mathbf{u}^H \mathbf{v}|^2}{|\mathbf{u}|^2 |\mathbf{v}|^2}, \qquad (34)$$

where $\mathbf{u}^H$ is the hermitian of $\mathbf{u}$. Even if this criteria does not constitute a scalar product for the considered basis, it nevertheless gives a good estimation of the closest shape to a given reference vector. The use of a correlation indicator to track the waves between two calculation steps gives the opportunity to plot confident dispersion diagrams, in particular when curves are crossing together that is to say when veering or bifurcation phenomena occur in dispersion curves. Moreover, in particular situations like apparition of a new

<figure><img src="image_14.png" /><figcaption>Fig. 13. Directivity of damped system using wave numbers in the first Brillouin zone in the case of 10% of damping ratio: each mark indicates a computed solution, the size of the mark is a measure of the propagative nature of the wave ("evanescent" waves correspond to small radius, "propagative" waves correspond to large radius).</figcaption></figure>
wave, loosing of correlation, or even for vanishing of group velocity, the frequency step is adapted in order to enhance the ability to follow the waves.

Fig. 9 illustrates the typical results of the analysis. Propagative wave numbers of the damped system are shown for $\phi = 0$ and $\phi = 18^\circ$. It can be observed that if $\phi = 0$, the symmetry illustrated in Fig. 8 still exists, while as soon as other directions are considered, the symmetry in the dispersion diagram does not exist anymore. This can be explained by the fact that the periodicity of the initial pattern is lost when the orientation is not parallel to one of the sides of the initial cell. Concerning the correlation, some surprising results can be observed: in some cases the correlation indicator fails to follow a given mode, even for small frequency steps. It is not yet clear if this is a numerical artefact or if this can be explained physically. One should emphasized that the MAC-based correlation is not mathematically justified since it does not constitute a scalar product for the considered basis. This point is currently under investigation.

Fig. 10 illustrates the stop bands directivity of the damped system using the evanescence indicator saturated at unit value for a sake of visualization. The full dark areas correspond to stop bands in which only evanescent waves can exist. The stop bands can exist even in the case of lightly damped system. These bands can be angle-dependent and exist only for particular directions.

A second computation has been made using a damping ratio of 10%. It can be observed in Fig. 11 that the evanescence ratio is modified and that larger band gaps can be observed. In that case, the existence of specific frequencies for which wave propagation is only possible along the main system axis (Ox) and (Oy) is observed. This particular behavior is reinforced by the damping effect smoothing the stop wave domains as observed in Fig. 10. The high frequency behavior is almost always 'evanescent' and cancels all wave transmissibility inside the system.

These results can also be interpreted visually by using directivity diagrams in which the solutions are plotted in terms of wave numbers, for a given set of frequency points. This graphical representation is shown in Fig. 12 for a damping ratio of 1%. One can observe in particular that for frequencies in pass-bands several propagative solutions can be found, while for frequencies in stop-bands, all solutions correspond to evanescent waves. This can clearly be observed at 49 kHz for example. For particular frequencies like 179 kHz, the directivity is strong: waves traveling along x and y axes are propagative, while when $\phi$ belongs to $[30^\circ,60^\circ]$, the corresponding waves are becoming evanescent.

When the damping ratio becomes larger, the real part of the solutions are increasing in amplitude and the waves tend to be attenuated, as it can be observed in Fig. 13. The directive effect observed at 179 kHz for a lower value of damping does no longer appear for that particular frequency since the former "propagative" waves have been changed in evanescent ones.

All these aspects constitute interesting results that will be confronted to experimental results in the next part of this work.

## 4. Concluding remarks

This paper presents a validated numerical procedure able to compute the damped wave's dispersion functions in the whole first Brillouin domain of multi dimensionnal elastodynamical wave's guides. The method has been applied for determining the bidimensional band-gaps of the well known periodic structures studied by Wu et al. (2009) when damping effect is considered. Based on this approach, a suitable criterion indicating the evanescence ratio of computed waves is proposed. It can be used for optimizing structured damping layers or electronics circuits and transducers for controlling vibroacoustic behavior of the systems. The damping

operator introduced in the formulation can be frequency dependent as viscous one but can also be much more complicated. It can compass specific dissipation phenomenon such as those induced by distributed shunted piezoelectric patches as proposed by Beck(2008), Casadei et al. (2009), or even foams or complex polymers behaviors. The proposed method furnishes an efficient tool for future optimization of distributed smart cells as proposed in the case of 1D waveguides by Collet et al. (2009).

## Acknowledgement

This work was carried out with a grant of French agency ANR number NT09-617542. The authors gratefully acknowledge the French ANR and CNRS for supporting this program.

## References

* Aberg, M., Gudmundson, P., 1997. The usage of standard finite element codes for computation of dispersion relations in materials with periodic microstructure. Journal of the Acoustical Society of America 102 (4), 2007-2013.

* Akrout, S., 2005. Comportement dynamique deterministe et large bande des structures guidees. Ph.D. Thesis, +cole Centrale Lyon.

* Allaire, G., Congas, C., 1998. Bloch waves homogenization and spectral asymptotic analysis. Journal de Mathématiques Pures et Appliquées 77, 153-208.

* Banks, H., Smith, R.C., Wang, Y., 1996. Smart material structures Modeling Estimation and Control. Masson and Wiley.

* Bartoli, I., Marzani, A., di Scalea, F.L., Viola, E., 2006. Modeling wave propagation in damped waveguides of arbitrary cross-section. Journal of Sound and Vibration 295, 685-707.

* Beck, B., Cunefare, K., Ruzzene, M., 2008. Broadband vibration suppression assessment of negative impedance shunts. In: Proceedings of SMASIS08. vol. 6928. ASME.

* Bensoussan, A., Lions, J., Pananicolaou, G., 1978. Asymptotic Analysis for Periodic Structures. North Holland.

* Bloch, F., 1928. Über die Quantenmechanik der Electron in Kristallgittern. Zeitschrift für Physik 52, 550-600.

* Brillouin, L., 1953. Wave Propagation in Periodic Structures. Dover Publications.

* Casadei, F., Ruzzene, M., Beck, B., Cunefare, K., 2009. Vibration control of plates featuring periodic arrays of hybrid shunted piezoelectric patches. In: Proceedings of SPIE - Smart Structures and Materials. vol. 7288. SPIE.

* Collet, M., Cunefare, K., Ichchou, N., 2009. Wave motion optimization in periodically distributed shunted piezocomposite beam structures. Journal of International Mat System and Structure 20 (7), 787-808.

* Damljanovic, V., Weaver, R.L., 2004. Propagating and evanescent elastic waves in cylindrical waveguides of arbitrary cross section. Journal of the Acoustical Society of America 115, 1572-1581.

* Duhamel, E., 2007. Finite element computation of green's functions. Engineering Analysis with Boundary Elements 31, 919-930.

* Finnveden, S., 2004. Evaluation of modal density and group velocity by a finite element method. Journal of Sound and Vibration 273, 51-75.

* Floquet, G., 1883. Sur les équations diffTrentielles linéaires a coefficients périodiques. Annales de l'Ecole Normale Supérieure 12, 47-88.

* Gavrić, L., 1994. Finite element computation of dispersion properties of thin-walled waveguides. Journal of Sound and Vibration 173 (1), 113-124.

* Gavrić, L., 1995. Computation of propagative waves in free rail using a finite element technique. Journal of Sound and Vibration 185, 531-543.

* Hayashi, T., Song, W.J., Rose, J.L., 2003. Guided wave dispersion curves for a bar with an arbitrary cross-section, a rod and rail example. Ultrasonics 41, 175-183.

* Houillon, L., Ichchou, M., Jezequel, L., 2005. Wave motion in thin-walled structures. Journal of Sound and Vibration 281 (3-5), 483-507.

* Ichchou, M.N., Akrout, S., Mencik, J., 2007. Guided waves group and energy velocities via finite elements. Journal of Sound and Vibration 305 (4-5), 931-944.

* Ichchou, M., Berthaut, J., Collet, M., 2008. Multi-mode wave propagation in ribbed plates : Part I k-space characteristics. International Journal of Solids and Structures 45 (5), 1179-1195.

* Joannopoulos, J., Meade, R., Winn, J., 1995. Photonic Crystals: Molding the Flow of Light. Princeton University Press.

* Kittel, C., 1986. Introduction to Solid State Physics. John Wiley and Sons, New York.

* Lehoucq, R., Sorensen, D., Yang, C., 1998. ARPACK users' guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods. Siam.

* Mace, B., Manconi, E., 2008. Modelling wave propagation in two-dimensional structures using finite element analysis. Journal of Sound and Vibrations 318, 884-902.

* Manconi, E., 2008. The wave finite element method for 2-dimensional structures. Ph.D. Thesis, University of Parma.

* Manconi, E., Mace, B., 2010. Estimation of the loss factor of viscoelastic laminated panels from finite elements analysis. Journal of Sound and Vibration 329, 3928-3039.

* Mead, D., 1996. A general theory of harmonic wave propagation in linear periodic systems with multiple coupling. Journal of Sound Vibration 27 (2), 429-438.
Mencik, J., Ichchou, M., 2005. Multi-mode propagation and diffusion in structures through finite elements. *European Journal of Mechanics A-Solids* 24(5), 877-898.

Nelson, P.A., Elliott, S.J., 1992. Active Control of Sound. Academic Press, London, San Diego.

Preumont, A., 1997. Vibration control of structures : An introduction. Kluwer.

Schenk, O., Gärtner, K., 2004. Solving unsymmetric sparse systems of linear equations with PARDISO. *Future Generation Computer Systems* 20 (3), 475-487.

Thorp, O., Ruzzene, M., Baz, A., 2001. Attenuation and localization of wave propagation in rods with periodic shunted piezoelectric patches. *Proceedings of SPIE – The International Society for Optical Engineering Smart Structures and Materials* 4331, 218-238.

Wilcox, C., 1978. Theory of bloch waves. *Journal d'Analyse Mathématique* 33, 146-167.

Wu, T., Wu, T., Hsu, J., 2009. Waveguiding and frequency selection of Lamb waves in plate with periodic stubbed surface. *Physical Review B* 79, 104306.