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

HAL
open science

# A Floquet-Bloch decomposition of the elastodynamical equations for the computation of wave's dispersion in damped mechanical systems

Manuel Collet, Morvan Ouisse, Massimo Ruzzene, Mohamed Ichchou

► To cite this version:

Manuel Collet, Morvan Ouisse, Massimo Ruzzene, Mohamed Ichchou. A Floquet-Bloch decomposition of the elastodynamical equations for the computation of wave's dispersion in damped mechanical systems. ISMA2010 International Conference and Vibration and Noise Engineering, Sep 2010, Louvain, Belgium. pp.2837 - 2848. hal-01512527

HAL Id: hal-01512527

https://hal.science/hal-01512527

Submitted on 25 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# A Floquet-Bloch decomposition of the elastodynamical equations for the computation of wave's dispersion in damped mechanical systems

M. Collet¹, M. Ouisse¹, M. Ruzzene², M. Ichchou³

¹ FEMTO-ST, Applied Mechanics Dept, UMR-CNRS 6174,
24 chemin de l'épitaphe
25000 Besançon, France

e-mail: manuel.collet@univ-fcomte.fr

² Georgia Institute of Technology
Aerospace School of Engineering
Atlanta, GA 30332-0405, USA

³ LTDS, UMR5513, Ecole Centrale de Lyon
36 avenue Guy de Collongue
69134 Ecully, France

## Abstract

In the context of mid-frequency elastodynamical analysis of periodic structures, the Floquet-Bloch theorem has been recently applied. The latter allows the use of limited-size modeling of a representative cell to characterize waves dispersion properties of the assembled structure. The theorem provides a rigorous and well posed spectral problem representing waves dispersion in viscoelastic media. In particular, the wave-FEM (WFE) method and associated techniques are based on the Floquet-Bloch approach. Most of the current applications are limited either to undamped or to slightly damped systems. Through the mesh of the k-space (k being the wave number vector) inside or at the boundary of the first Brillouin zone is performed. For each wave headings a linear spectral problem is solved in order to obtain the corresponding eigenfrequency and the corresponding Floquet vector. Such process is more easily implemented for undamped structures but face some drawbacks when damping is considered. In this paper, the Floquet-Bloch theorem is used to set up an alternative technique in order to estimate the dispersion characteristics of a periodic structure. The approach is based on the same basic assumption, but the generalized eigenvalue problem which has to be solved differs from the one usually considered in WFE. Instead of discretising the k-space, the harmonic frequency and the wave heading are used to scan the k-space. Then the formulated eigenvalue problem is solved and the dispersion characteristics are obtained, including spatial attenuation terms. Some fundamentals properties of the eigensolutions are discussed, and the methodology is finally applied on a 2D waveguide application which can be found in the literature for an undamped case. The same example is considered with various damping levels, in order to illustrate the performances and specificities the efficiency of proposed approach. The proposed approach finds application in the analysis of wave propagation in the presence of damping materials or shunted piezoelectric patches, as well as in actively controlled systems, where equivalent damping terms are associated with the considered control scheme.
# 1 Introduction

Tailoring the dynamical behavior of wave-guide structures can provide an efficient and physically elegant means to optimize mechanical components 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 direction. 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 then MF optimization. For example, piezoelectric materials and other adaptive and smart systems are employed to improve the vibroacoustic quality of structural components, especially in the LF range (Preumont 1997, Nelson & S.J. 1992, Banks, Smith & Wang 1996). 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 (Collet, Cunefare & Ichchou 2009, Thorp, Ruzzene & Baz 2001). However, there is still a lack of studies in the literature for MF optimization of structural vibration. To that end, the focus of this study is to provide a suitable numerical tool for computing wave dispersion in 2D periodic systems incorporating controlling electronics devices. The main final aim is to allow their optimization in order to optimize vibroacoustic diffusion in 2D wave's guides. Two numerical approaches can be distinguished for computing that dispersion: 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 (FE). 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, a large amount of FE has been developed since 1975 to compute dispersion curves of rails (Mead 1996), laminated composite plates (Mace & Manconi 2008, Gonella & Ruzzene 2008) and viscoelastic laminated composite plates (Akrout 2005). 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 (Mead 1996), an eigenvalue problem can be formulated from the stiffness and mass matrix of the FE model to find wave numbers of all the propagating waves. Contrary to SAFE method, the displacement field is now 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 (Mace & 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 (Manconi 2008, Berthaut, Collet & Ichchou 2008). One of the main problem all these approaches is the difficulty to compute the damped wave numbers in the whole Brillouin domain necessary for optimizing vibroacoustic behavior of such periodic structures. After recalling in the two first parts of this paper the Floquet then the Bloch theorems, we introduce a new numerical formulation for computing the multi-modal damped wave numbers in the whole first Brillouin domain of a periodical structures. Then a 2D numerical application is presented allowing us to validate the method and to use it for estimating the 2D band gaps as well as a suitable evanescence's indicator that could be used for control optimization.

# 2 Elasto-Dynamical application of the Floquet-Bloch Theorem

In this section the application of the celebrated Floquet-Bloch theorem is presented in the context of elasto-dynamic. First we recall the main results in a one-dimensional setting by Floquet (Floquet 1883) and later rediscovered by Bloch (Bloch 1928) in multidimensional problems. These results are recalled here since these references may be difficult to find. Application to bi-dimensional elastodynamical problem is proposed leading to very general numerical implementation for computing waves dispersion for periodically
distributed damped mechanical systems.

## 2.1 The Floquet Theorem

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

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

where is the unknown vector ($w(x) : \mathbb{R} \to \mathbb{C}^n$) and $A(x)$ is a given matrix of continuous periodic functions with period $r_1$, i.e. $A(x + r_1) = A(x)$.

The Floquet Theorem indicates that any solution of this linear system can be expressed as a linear combination of functions $V(x)e^{kx}$, where $V(x)$ is a $r_1$-periodic function and $k \in \mathbb{C}$ is a scalar complex value. The theory provides a way to evaluate $V$ and $k$ from the solution of an eigenvalue problem.

Among the many mathematical aspects of the theory, some points should be mentioned to the reader for proper understanding:

* Starting from any given basis $W(x) \in \mathbb{C}^n$ of fundamental solutions of (1), a new basis of solutions normalized so that $P_0(0) = I_n$ can be defined, $I_n$ being the $n \times n$ identity matrix:

$$ P_0(x) = W(x)W^{-1}(0). \qquad (2) $$

It is then possible to search for $W(x + r_1)$ on this basis. The calculations lead to:

$$ W(x + r_1) = P_0(x)W(r_1) = W(x)W^{-1}(0)W(r_1). \qquad (3) $$

$P_0$ is then called Floquet propagator since it allows one to evaluate $W(x + r_1)$ from the knowledge of $W(r_1)$. The next step in the theory is to provide a technique to obtain $P_0$.

* First of all, a diagonalization of the Floquet propagator is performed in $x = r_1$:

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

where $\Lambda$ and $Z$ are solutions of the following eigenproblem:

$$ \lambda_j Z_j = P_0(r_1) Z_j. \qquad (5) $$

$\Lambda$ is then a diagonal matrix with $\lambda_j$ terms, and $Z$ is the matrix grouping eigenvectors $Z_j$. The eigenvalues can also be written as:

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

where $k$ is a diagonal matrix whose generic term is $k_j$ such that $\lambda_i = e^{k_j r_1}$, $k_j$ being called the Floquet or characteristic exponents, while $\lambda_j$ are the Floquet or characteristic multipliers.

* The computation is not performed directly on $P_0$: a more convenient way to obtain the solution of the set of ordinary differential equations is to try to identify $Y$, which is the Floquet propagation of basis vectors $Z$ such that $Y(x) = P_0(x)Z$. It can be shown that:

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

$Y(x)$ are solutions of the initial periodic problem (1) restricted to elementary cell $[0, r_1]$, with fixed boundary conditions at $x = 0$ and $x = r_1$. So the eigenvectors $Z_j$ and eigenvalues $\lambda_j = e^{k_j r_1}$ are solutions of the following generalized eigenvalue problem:

$$ \left\{ \begin{array}{l} \displaystyle \frac{dY}{dx}(x) = A(x)Y(x) \quad \forall x \in [0, r_1], \\ Y(0) = Z, \\ Y(r_1) = Z\Lambda. \end{array} \right. \qquad (8) $$
The Floquet propagators are then obtained using the backward relationship:

$$P_0(x) = Y(x)Z^{-1}, \qquad (9)$$

and a basis of solutions is given by:

$$W(x) = Y(x)W(0). \qquad (10)$$

* An alternative way to compute the Floquet propagators is to introduce the undamped Floquet Vectors $V_i$ defined as:

$$V(x) = Y(x)e^{-kx}, \qquad (11)$$

where $V$ is the matrix regrouping the vectors $V_j$. It can be shown that these function are $r_1$ periodic, and that they are solutions of the following problem:

$$\left\{ \begin{array}{l} \frac{dV}{dx}(x) = A(x)V(x) - V(x)k \quad \forall x \in [0, r_1], \\ V(0) = Z, \\ V(r_1) = Z. \end{array} \right. \qquad (12)$$

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

$$P_0(x) = V(x)e^{kx}Z^{-1}, \qquad (13)$$

while a basis of solutions is given by:

$$W(x) = V(x)e^{kx}W(0). \qquad (14)$$

Equation (14) is called the Floquet Normal form of the fundamental basis $W(x)$. We underline on equation (14) that the characteristic multipliers are also the eigenvalues of the linear Poincaré maps defined as the function $w(x) \to w(x + r_1)$, $w(x)$ being solution of (1).

Some interesting properties can be derived from the theory:

* The Floquet exponents are not unique since $e^{(k+i\frac{2m\pi}{r_1})r_1} = e^{kr_1}$ if $m$ is an integer.

* The Floquet vectors are periodic, so they are bounded on all $\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.

## 2.2 The Bloch Theorem

The Bloch theorem gives the form of homogeneous states of Schrdinger equation with periodic potential. This theorem can be considered as a multidimensionnal application of the Floquet theorem (Joannopoulos, Meade & Winn 1995). The periodic medium (or potential) properties satisfy $M(x + R.m) = M(x)$, $m \in \mathbb{Z}$ where $R = [r_1, r_2, r_3] \in \mathbb{R}^{3 \times 3}$ the matrix grouping the three lattice's basis vectors (in 3D). We can also define the primitive cell as a convex polyhedron of $\mathbb{R}^3$ called $\Omega_x$. The reciprocal unit cell is denoted by $\Omega_k$ limited by the reciprocal lattice vector defined by the three vectors $g_j$ so that: $r_i.g_j = 2\pi\delta_{ij}$ ($\delta_{ij}$ the Kronecker index). We note $G = [g_1, g_2, g_3]$ 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. One can see (Kittel 1986) for details.
<figure><img src="image_3.png" /><figcaption>Figure 1: Generic 3D periodic cells</figcaption></figure>

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

$$ u(x) = \int_{\Omega_k} e^{ikx} \tilde{u}(x, k) dk \quad (15) $$

where the Bloch amplitude $\tilde{u}(x, k)$ is $\Omega_x$-periodic in $x$ and has the representations

$$ \tilde{u}(x, k) = \sum_{n \in \mathbb{Z}^3} \hat{u}(k + Gn)e^{iGn \cdot x}, \quad u(x) = \frac{|\Omega_x|}{(2\pi)^3} \sum_{n \in \mathbb{Z}^3} u(x + Rn)e^{ik(x+Rn)} \quad (16) $$

where $\hat{u}(k)$ stands for the Fourier transform of $u(x)$. One can also demonstrate that the mean value of the Bloch amplitude is the Fourier amplitude of $u(x)$ for the corresponding wave vector: $\langle \tilde{u}(., k) \rangle_{\Omega_x} = \hat{u}(k)$. Using the Bloch theorem to represent the solutions of periodical partial derivative equations implies that all derivatives are shifted by $k$ in the sense given by the used spatial operator.

Based on that theorem one can define the expansion functions $v_m(x, k)$, called the Bloch eigenmodes, such they can be used to represent the Bloch amplitudes of any solution of the corresponding partial derivative equation as

$$ \tilde{u}(x, k) = \sum_{m} u_m(k) v_m(x, k) \quad (17) $$

and at the same time diagonalize the partial derivative equations. One notes that the expansion coefficients $u_m(k)$ depend on the applied disturbance and also on the induced wave vector.

## 2.3 Application to Elastodynamic

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

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

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

$$ w(x) = w_{n,k}(x, k)e^{ik \cdot x} \quad (19) $$
<figure><img src="image_4.png" /><figcaption>Figure 2: (a) Rectangular parallelepiped primitive Lattice (b) Corresponding rectangular parallelepiped reciprocal Lattice</figcaption></figure>

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

$$ \rho\omega_n(\mathbf{k})2w_{n,k}(\mathbf{x}) + \nabla C \nabla_{\text{sym}}(w_{n,k}(\mathbf{x})) -iC\nabla_{\text{sym}}(w_{n,k}(\mathbf{x})).\mathbf{k} -i\nabla C\frac{1}{2}(w_{n,k}(\mathbf{x}).\mathbf{k}^T + \mathbf{k}.w_{n,k}^T(\mathbf{x})) +C\frac{1}{2}(w_{n,k}(\mathbf{x}).\mathbf{k}^T + \mathbf{k}.w_{n,k}^T(\mathbf{x})).\mathbf{k} = 0 \quad \forall \mathbf{x} \in \Omega_R \quad (20) $$

$$ w_{n,k}(\mathbf{x} - \mathbf{R}.\mathbf{n}) = w_{n,k}(\mathbf{x}) \quad \forall \mathbf{x} \in \Gamma_R \quad (21) $$

The first equation is simply obtained by introducing equation (19) into elastodynamic equation (18). The second equation represents the symmetrical boundary conditions expressed on boundary faces of the lattice polyhedron as described on figure 2(a).a for a rectangular parallelepiped cell. In this equation $n$ stands for the outpointing unitary normal vector. It corresponds to a complex Quadratic Eigenvalue Problem (QEO) that can be solved neither 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 also based on the computation of the Floquet vectors (see equations (12) for 1D application or (20) in 3D), instead of computing the Floquet propagators (as in equation (8)) commonly used for elastodynamic applications. Our approach allows to obtain the full 3D waves dispersions functions and to clearly introduce damping part in the elastodynamic operator. The adopted methodology allows the computation of the complete complex map of the dispersion curves incorporating computation of evanescent waves and allowing the introduction of damping operator if any.

## 3 Weak Formulation and computation of waves dispersion functions in periodical lattice

Let us consider the partial derivative equations (20) on a unit cell $\Omega_R$. It stands for a generalized eigenvalue problem leading to compute the dispersion functions $\omega_n(\mathbf{k})$ and the corresponding Floquet eigenvectors $w_{n,k}(\mathbf{x})$.
For 3D applications, the wave vectors supposed to be complex if damping terms are added into equation (20), can be written as $k = k \begin{bmatrix} \sin(\theta)\cos(\phi) \\ \sin(\theta)\sin(\phi) \\ \cos(\theta) \end{bmatrix}$ where $\theta, \phi$ represent the direction angles into the reciprocal lattice domain as shown in figure 2(b). This decomposition assumes that real and imaginary parts of vector $k$ are co-linear. In the following, we note $\Phi = \begin{bmatrix} \sin(\theta)\cos(\phi) \\ \sin(\theta)\sin(\phi) \\ \cos(\theta) \end{bmatrix}$ that direction vector.

## 3.1 Weak Formulation

If $w_{n,k}(\mathbf{x})$ is a solution of equations (20), also $\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\}$ we have:

$$ \int_{\Omega_R} \rho \omega 2n(\mathbf{k}) \tilde{w}_{n,k}(\mathbf{x}) w_{n,k}(\mathbf{x}) - \tilde{\varepsilon}_{n,k}(\mathbf{x}) C \varepsilon_{n,k}(\mathbf{x}) \\ +ik\tilde{\kappa}_{n,k}(\mathbf{x}) C \varepsilon_{n,k}(\mathbf{x}) - ik\tilde{\varepsilon}_{n,k}(\mathbf{x}) C \kappa_{n,k}(\mathbf{x}) + k2\tilde{\kappa}_{n,k}(\mathbf{x}) C \kappa_{n,k}(\mathbf{x}) d\Omega \\ + \int_{\Gamma_R} \tilde{w}_{n,k}(\mathbf{x}) (C(\varepsilon_{n,k}(\mathbf{x}) + ik\kappa_{n,k}(\mathbf{x}))). \mathbf{n} d\Gamma = 0 \quad (22) $$

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

This weak formulation is simply obtained by integrating equation (20) projected onto any test function $\tilde{w}_{n,k}(\mathbf{x})$. The boundary integral vanishes as the test functions are chosen so that $\tilde{w}_{n,k}(\mathbf{x} - R\mathbf{n}) = \tilde{w}_{n,k}(\mathbf{x})$ on $\Gamma_R$ that implies $C(\varepsilon_{n,k}(\mathbf{x}+R)+ik\kappa_{n,k}(\mathbf{x}+R)).\mathbf{n}(\mathbf{x}+R) = -C(\varepsilon_{n,k}(\mathbf{x})+ik\kappa_{n,k}(\mathbf{x})).\mathbf{n}(\mathbf{x})$. That corresponds to the exact compensation of the boundary applied generalized constrains $C(\varepsilon_{n,k}(\mathbf{x})+ik\kappa_{n,k}(\mathbf{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 called $w_{n,k}(\mathbf{x} - R\mathbf{n}) = w_{n,k}(\mathbf{x})$ explicitly link these associated surfaces. As the corresponding normal vector $\mathbf{n}$ are opposite, $\kappa_{n,k}(\mathbf{x} + R) = \kappa_{n,k}(\mathbf{x})$ and the stress condition can be restricted to $C(\varepsilon_{n,k}(\mathbf{x} + R)).\mathbf{n}(\mathbf{x} + R) = -C(\varepsilon_{n,k}(\mathbf{x})).\mathbf{n}(\mathbf{x})$ on the two opposite surfaces. Thus, all boundary integrations vanish and the used weak formulation is:

$$ \begin{array}{l} \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\} \\ \int_{\Omega_R} \rho \omega 2n(\mathbf{k}) \tilde{w}_{n,k}(\mathbf{x}) w_{n,k}(\mathbf{x}) - \tilde{\varepsilon}_{n,k}(\mathbf{x}) C \varepsilon_{n,k}(\mathbf{x}) \\ +ik\tilde{\kappa}_{n,k}(\mathbf{x}) C \varepsilon_{n,k}(\mathbf{x}) - ik\tilde{\varepsilon}_{n,k}(\mathbf{x}) C \kappa_{n,k}(\mathbf{x}) + k2\tilde{\kappa}_{n,k}(\mathbf{x}) C \kappa_{n,k}(\mathbf{x}) d\Omega = 0 \quad (23) \end{array} $$

## 3.2 Numerical Computation

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

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

where $\lambda = ik$, $M$ and $K$ are respectively the standard symmetric definite mass and symmetric semi-definite stiffness matrices coming from $\int_{\Omega_R} \rho \omega_n^2(\mathbf{k}) \tilde{w}_{n,k}(\mathbf{x}) w_{n,k}(\mathbf{x}) d\Omega$ and $\int_{\Omega_R} \tilde{\varepsilon}_{n,k}(\mathbf{x}) C \varepsilon_{n,k}(\mathbf{x}) d\Omega$ terms
in the weak formulations. L is an skew-symmetric matrix associated with term $\int_{\Omega_R} -\tilde{\kappa}_{n,k}(\boldsymbol{x})C\varepsilon_{n,k}(\boldsymbol{x}) + \tilde{\varepsilon}_{n,k}(\boldsymbol{x})C\kappa_{n,k}(\boldsymbol{x})d\Omega$ and H is a symmetric semi-definite positive matrix linked to $\int_{\Omega_R} \tilde{\kappa}_{n,k}(\boldsymbol{x})C\kappa_{n,k}(\boldsymbol{x})d\Omega$.

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

This approach has been widely used for developing homogenization techniques and spectral asymptotic analysis (Allaire & Congas 1998). It can also be applied for computing wave's dispersion even if Floquet propagators is preferred for 1D or quasi 1D computation (Ichchou, Akrout & Mencik 2007, Houillon, Ichchou & Jezequel 2005, Mencik & Ichchou 2005). Nevertheless these approaches has been only developed for undamped mechanical system that is to say representing 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 λ and Φ) inside the first Brillouin zone for obtaining the corresponding frequency dispersion and the associated Floquet vectors. For undamped system only propagative or evanescent waves exist corresponding to a family of eigen solutions purely real or imaginary. Discrimination between each class of waves is easy. If damped system is considered, that is to say if matrices K, L, H are complex, evanescent part of propagating waves appears as the imaginary part of $\omega_n^2(\lambda, \Phi)$ and vice versa. It also becomes very difficult to distinguish the two family of wave but also to compute the corresponding physical wave's movements by applying spacial deconvolution.

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

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

In this problem, the pulsation $\omega$ is a real parameter corresponding to the harmonic frequency. Wave's numbers and Floquet vector are also computed. 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 (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 pairs 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.

As previously mentioned, the real part of $k = k\Phi$ vector is restricted to stand inside the first Brillouin zone (cf figure 2(b)). In the quadratic eigen value problem (25) nothing restricts computation to only find eigen values satisfying this condition. For direction vector $\Phi$ orthogonal to the lattice facelets (i.e for $\Phi_{p_1} = [1, 0]^T$ and $\Phi_{p_2} = [0, 1]^T$ in 2D rectangular cell) we have the same periodical conditions as expressed for one dimensional wave guide in equation (11): if $\lambda_i(\Phi_p)$ is an eigen value associated to $w_{i,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_{i,k}(\Phi_p)e^{-i \cdot \Phi_p^T(G \cdot m)x}$. Thus, for undamped system, all obtained eigenvalues are periodically distributed in the k-space along its principal directions.

## 4 Applications for computing 2D waves dispersion

We only present here 2D wave guide application. Thus, we can easily find on literature comparative work to validate this new computational methodology. Two different systems are considered in this section. The first one corresponds to the undamped system used in (Wu, Wu & Hsu 2009) to validate our computation and the second one correspond to the damped version of the same system.
<figure><img src="image_5.png" /><figcaption>Figure 3: (a) Schematic of the considered 2D wave guide (b) Description of the unit cell with periodic stubbed surface</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: 2D physical cell and the corresponding first Brillouin zone. Shaded area is its irreducible part</figcaption></figure>

**4.1 Undamped wave dispersion and band-gap computation in thin plate with periodic stubbed surface**

The considered system is exactly those treated in the work of Wu (Wu et al. 2009). It consists of an infinite periodic 2D wave guide as described in 3(a). The system is made of a 1 mm thick thin aluminum plate with periodic cylindrical stubs on one of its face as described on 3(b). The used material for the whole system is an isotropic Aluminum 6063-T83 ($\nu = 0.33$, $E = 69e9[Pa]$, $\rho = 2700[kg/m3]$).

By using symmetry of the unit cell, the corresponding first Brillouin zone is described on figure 4 where the irreducible zone is the shaded area. Our computation method allows us to compute eigen frequencies corresponding to any *k* vector described in cylindrical coordinates system by its radius *k* and its angle *φ*.

**4.1.1 Numerical implementation**

The numerical implementation is based on the 3D weak formulation (23), using a 2D 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 $\Gamma_{r_1^+} \Gamma_{r_1^-}$ and $\Gamma_{r_2^+} \Gamma_{r_2^-}$ as depicted in figure 3(b). The implementation is made with COMSOL Multiphysics® platform and parametric computation allowing to obtain $k(\omega, \phi)$ is carried out with Matlab® routines.
<figure><img src="image_7.png" /><figcaption>Figure 5: Unrefined (a) and refined (b) Meshes case</figcaption></figure>

For each parameter $\omega$ and $\phi$, one reformulates the quadratic eigen value problem as a linear 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 (Lehoucq, Sorensen & Yang 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 they accomplish by solving the linear system $Ax = Bv$ using here a Pardiso solver (Schenk & Gärtner 2002).

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

The used mesh is shown in figure 4.1.1. 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.

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

The first computation has been made to compare our algorithm with the previously obtained results presented in (Wu et al. 2009). Let us apply our method for computing the wave's dispersion functions of the undamped system along the $\Gamma - X$ direction (i.e for $\phi = 0$). We present in figure 6 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 (Aberg & Gudmundson 1997, Mace & Manconi 2008) based on equation (20). The second and third dispersion curves (in dotted and crossed lines in the figure 6) show, respectively, the obtained results with the unrefined and refined meshes cases. The results show a really good agreement between a standard computation method (Aberg & Gudmundson 1997, Mace & Manconi 2008) used to obtained the reference results as proposed by (Wu et al. 2009) and our method with the refined mesh. We point 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 show up the convergence of the refined model compare to the unrefined one. These computations validates the numerical implementation of the proposed method.

For evaluating the band-gap of the periodic system, we propose to use an indicator of the minimal evanescence ratio of all the computed waves for each considered frequencies, defined as :

$$Ind(\omega, \phi) = \min_{n} \left| \frac{Real(\lambda_{n})}{|(\lambda_{n})|} \right| \qquad (26)$$

We present in figure 7 the plot of this indicator for both mesh cases. We can also figure out location of the first two stop bands of the system : the first one is from 40 to 50 MHz and the second from 156 to 176 MHz. Precision of these results depends on the frequency discretization rate. The obtained band-gap is
<figure><img src="image_8.png" /><figcaption>Figure 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_9.png" /><figcaption>Figure 7: Evanescence ratio $Ind(\omega)$. plain lines : refined mesh, dashed line : unrefined mesh</figcaption></figure>

totally comparable with those computed by Wu (Wu et al. 2009). Figure 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 figure 4. The obtained wavenumbers are symmetrical according to the boundary conditions of the corresponding polyhedral surface. This properties is observed on figure 8 where the whole set of obtained wavenumbers (i.e the imaginary parts of $\lambda_n(\omega)$) is plotted. We observe that they are symmetrical with respect to the vertical axes on $\pm\frac{\pi}{dx} = 100\pi$ when $\phi = 0$.

## 4.2 Dispersion of the damped system in the whole 2D K-space

The proposed computational method allows us to compute multi modal wave's propagations in the complete 2D K-space in the first Brillouin zone. The proposed methodology is based on the direct computation of complex wave numbers as a function of frequency. The Bloch theorem is expanded in the case of damped system 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 equation (5).
<figure><img src="image_10.png" /><figcaption>Figure 8: Whole propagative wave numbers ($imag(\lambda_n(\omega))$) when $\phi = 0$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9: Propagative wave numbers of damped system ($imag(\lambda_n(\omega))$) when $\phi = 0$ (left) and $\phi = 18^\circ$</figcaption></figure>

The damping behavior is introduced by assuming a complex Hooke elasticity tensor. The same methodology could have been realized by introducing any kind of linear viscoelastic modeling such as viscous behavior or any other complex frequency-dependent terms.

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 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. When new waves appear or some of them loose correlation, the frequency step is adapted in order to enhance the ability to follow the waves.

The figure 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 figure 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 side 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 base. This point is currently under investigation.

Figure 10 illustrates the stop bands directivity of the damped system using the evanescence indicator sat-
<figure><img src="image_12.png" /><figcaption>Figure 10: Directivity of damped system using evanescence indicator saturated at unit value</figcaption></figure>

urated at unit value for a sake of visualization. The full red areas correspond to stop band 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.

# 5 Concluding remarks

We propose here a validated numerical procedure able to compute the damped wave's dispersion functions in the whole first Brillouin domain of elastodynamical wave's guides. The method was applied for determining the 2D band-gaps of the well known 2D periodic structures studied by (Wu et al. 2009) when damping is considered. Based on this approach, we also propose a suitable criterion indicating the evanescence ratio of computed waves. It can be also used for optimizing electronics circuits and transducers for controlling vibroacoustic behavior of such system. The introduced damping operator can be frequency dependent as a viscous one but can also be much more complicated. It can compass particular dissipation phenomenoms such as those induced by distributed shunted piezoelectric patches as proposed in (Beck, Cunefare & Ruzzene 2008, Casadei, Beck & Ruzzene 2009). The proposed method furnishes an efficient tool for future optimization of distributed shunted piezoelectric cells as proposed in the case of 1D wave guide in (Collet et al. 2009).

# 6 Acknowledgement

This work was carried out with a grant of French agency ANR number NT09-617542. We gratefully acknowledge Georgia Tech and the French ANR and CNRS for supporting such international collaborations.
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 dterministe et large bande des structures guides*, PhD thesis, cole Centrale Lyon.

Allaire, G. & Congas, C. (1998). Bloch waves homogenization and spectral asymptotic analysis, *Journal de Mathmatiques Pures et Appliques* **77**: 153–208.

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

Beck, B., Cunefare, K. A. & Ruzzene, M. (2008). Broadband vibration suppression assessment of negative impedance shunts, in *ASME (ed.), Proceedings of SMASIS08*.

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

Bloch, F. (1928). ber die Quantenmechanik der Electron in Kristallgittern, *Zeitschrift fr Physik* **52**: 550–600.

Casadei, F., Beck, B. & Ruzzene, M. (2009). Vibration control of plates featuring periodic arrays of hybrid shunted piezoelectric patches, in *SPIE (ed.), Proceedings of SPIE - Smart Structures and Materials*.

Collet, M., Cunefare, K. & Ichchou, N. (2009). Wave Motion Optimization in Periodically Distributed Shunted Piezocomposite Beam Structures, *Journal of Int Mat Syst and Struct* **20**(7): 787–808.

Floquet, G. (1883). Sur les équations diffrentielles linéaires á coefficients périodiques, *Annales de l'Ecole Normale Supérieure* **12**: 47–88.

Gonella, S. & Ruzzene, M. (2008). Analysis of in-plane wave propagation in hexagonal and re-entrant lattices, *J. Sound Vibrat.* **312**: 125–139.

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.

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*, SIAM Philadelphia, PA.

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*, PhD thesis, University of Parma.

Mead, D. (1996). A general theory of harmonic wave propagation in linear periodic systems with multiple coupling, *J. Sound Vibrat.* **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. & S.J., E. (1992). *Active Control of Sound*, Pub. Academic Press, London, San Diego.

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

Schenk, O. & Gärtner, K. (2002). Solving unsymmetric sparse systems of linear equations with PARDISO, *Computational ScienceICCS 2002* pp. 355–363.

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.

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.