# Méthode de réduction de modèles pour les vibrations de plaques sandwich viscoélastiques

Faiza Boumédiene, Jean-Marc Cadou, Laëtitia Duigou, El Mostapha Daya

► To cite this version:

Faiza Boumédiene, Jean-Marc Cadou, Laëtitia Duigou, El Mostapha Daya. Méthode de réduction de modèles pour les vibrations de plaques sandwich viscoélastiques. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717095

HAL Id: hal-01717095

https://hal.science/hal-01717095v1

Submitted on 25 Feb 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Méthode de réduction de modèles pour les vibrations de plaques sandwich viscoélastiques

Faiza Boumédiene<fnref n="1" />, Jean-Marc Cadou<fnref n="2" />*, Laëtitia Duigou<fnref n="2" />, El Mostapha Daya<fnref n="3" />

**Résumé** — Ce travail concerne le calcul des propriétés amortissantes de structures sandwich viscoélastiques. La résolution de ce problème aux valeurs propres non linéaires nécessite des techniques numériques performantes, telles que les méthodes de Newton d'ordre élevé. Cette dernière consiste en l'association d'une technique d'homotopie et d'une méthode de perturbation. La méthode de Newton d'ordre élevé est très efficace [8] mais demande des temps de calculs relativement importants. Une technique de réduction de modèles est alors proposée pour diminuer ces temps de calculs tout en conservant les propriétés de convergence de la méthode de Newton d'ordre élevé. Les applications numériques concernent des plaques sandwich viscoélastiques.

**Mots clés** — Structure sandwich, viscoélasticité, réduction de modèles, méthodes numériques

## 1. Introduction

Pour réduire les vibrations et le bruit, les industriels proposent des structures sandwich amortissantes, dans lesquelles une couche de matériau viscoélastique est comprise entre deux couches de matériau métallique, tel que l'aluminium. L'analyse des vibrations libres de ces structures amène à résoudre un problème aux valeurs propres complexes et non-linéaires. La non-linéarité est due au matériau viscoélastique de la couche centrale qui dépend de la fréquence. De nombreux modèles et méthodes ont été développés pour résoudre ce problème non-linéaire. Par exemple, Daya et Potier-Ferry [6] proposèrent une méthode de continuation. A partir des modes non-amortis et en utilisant la technique de perturbation, ils déterminent les fréquences et amortissements. Duigou et al. [8] développèrent deux méthodes itératives basées sur la technique d'homotopie, la méthode asymptotique numérique et les approximants de Padé. L'une de ces techniques est une sorte de méthode de Newton d'ordre élevé. D'autres auteurs, comme Damanpack et al. [5] ou Banerjee et al. [1] utilisèrent la méthode de rigidité dynamique. Cependant, quelque soit la méthode utilisée, la taille des matrices à manipuler dans un problème de vibrations non-linéaires peut devenir rapidement très grande et donc amener à des temps de calculs importants. Pour réduire ce coût, différentes techniques existent comme la méthode de condensation proposée par Park et al. [10] ou Lima et al. [9], l'algorithme de Lanczos utilisée par Chen [4].

Dans ce travail, nous considérons la méthode itérative d'ordre élevé proposée dans Duigou et al [8]. Cette méthode est très efficace mais elle demande un temps de calcul important principalement du à la taille des matrices à inverser. C'est pourquoi nous proposons d'utiliser une technique de réduction de modèles pour diminuer ces temps de calculs.

<footnote n="1">1 Laboratoire de Mécanique Avancée, Faculté de Génie Mécanique & Génie des Procédés, USTHB, BP 32, El Alia, 16111
Bab Ezzouar, Alger, Algeria, boumedienefaiza@yahoo fr</footnote>

<footnote n="2">2 Laboratoire d'Ingénierie des Matériaux de Bretagne, Université Européenne de Bretagne, Université de Bretagne Sud, Rue de Saint Maudé, BP 92116, 56321 Lorient Cedex, France, {jean-marc.cadou;Laetitia.duigou}@univ-ubs.fr</footnote>

<footnote n="3">3 LEM3, Ile du Saulcy 57045 Metz - cedex 01, daya@lpmm.sciences.univ-metz.fr</footnote>

<footnote>* Auteur correspondant</footnote>
## 2. Problème de vibration

On considère une plaque sandwich viscoélastique dont la couche centrale est composée d'un matériau viscoélastique. Ce dernier induit un module d'Young complexe. L'équation à résoudre pour déterminer les valeurs propres de ce type de structure sandwich est alors écrite sous forme discrète:

$$ R_0 U = [K_0 + E(p) K_v] U = 0 \quad (1) $$

Où **M** est la matrice de masse, **K**(0) est la matrice de rigidité élastique, *E*(*ω*) est le module d'Young dépendant de la pulsation *ω* et **U** représente le vecteur des déplacements nodaux généralisés. La matrice **K**<sub>v</sub> ne dépend pas de *ω*. Duigou [7] a proposé une méthode numérique pour résoudre le problème aux valeurs propres non linéaires (1). Cette méthode est basée sur une technique d'homotopie et une méthode de perturbation. Le principe est le suivant :

1- Un point de départ (**U**₀, *p*₀) avec ( *p* = *iω*) est choisi et une solution incrémentale est introduite de la façon suivante :

$$ \begin{cases} \Delta \mathbf{U} \\ \Delta p \end{cases} = \begin{cases} \mathbf{U} - \mathbf{U}_0 \\ p - p_0 \end{cases} \qquad (2) $$

2- Après insertion de (2) dans (1) et quelques manipulations, la technique d'homotopie est appliquée. Celle-ci consiste à introduire un paramètre *ε* tel que pour *ε* = 0, la correction (Δ**U**, Δ*p*) est égale à zéro et pour *ε* = 1, le problème non-linéaire exact (1) est résolu.

Ce problème modifié peut par exemple être défini par :

$$ \mathbf{L} \cdot \Delta \mathbf{U} + \mathbf{G} \cdot \mathbf{U}_0 \Delta p + 2 p_0 \Delta p \mathbf{M} \Delta \mathbf{U} + \varepsilon (\Delta E(p) - E'(p_0) \Delta p) \mathbf{K}_v \cdot \mathbf{U}_0 + \Delta E(p) \mathbf{K}_v \cdot \Delta \mathbf{U} + (\Delta p)^2 \mathbf{M} \cdot \mathbf{U}_0 + (\Delta p)^2 \mathbf{M} \cdot \Delta \mathbf{U} + \varepsilon \mathbf{R}_0 = 0 \quad (3) $$

Avec $\mathbf{L} = \mathbf{K}(0) + E(p_0)\mathbf{K}_v + p_0^2\mathbf{M}$ opérateur linéaire et $\mathbf{R}_0 = (\mathbf{K}(0) + E(p_0)\mathbf{K}_v + p_0^2\mathbf{M})\cdot\mathbf{U}_0$, résidu initial.

3- La technique de perturbation est enfin appliquée. C'est à dire que l'incrément (Δ**U**, Δ*p*) est cherché sous la forme de séries tronquées en fonction du paramètre de perturbation *ε* :

$$ \Delta \mathbf{U} = \sum_{i=1}^{N} \varepsilon^i \mathbf{U}_i \quad \text{et} \quad \Delta p = \sum_{i=1}^{N} \varepsilon^i p_i \qquad (4) $$

4- Après insertion de ces séries dans le problème modifié (3), une série de problèmes linéaires à chaque ordre est obtenue et résolue.

Comme le nombre d'inconnues est plus grand que le nombre d'équations, une équation pour normaliser le mode est introduite :

$$ \Delta \mathbf{U}^t \mathbf{G} \mathbf{U}_0 = 0 \qquad (5) $$

Alors, les problèmes linéaires peuvent être écrits au ième ordre de troncature, sous la forme matricielle suivante :

$$ \begin{bmatrix} \mathbf{L}_t \\ \mathbf{U}_i \\ p_i \end{bmatrix} = \begin{bmatrix} \mathbf{L} & (\mathbf{G} \mathbf{U}_0) \\ (\mathbf{G} \mathbf{U}_0)^T & 0 \end{bmatrix} \begin{bmatrix} \mathbf{U}_i \\ p_i \end{bmatrix} = \begin{bmatrix} \mathbf{F}_i^{nl} \\ 0 \end{bmatrix} \quad (6) $$

$$ \text{Avec } \mathbf{L} = \mathbf{K}(0) + E(p_0)\mathbf{K}_v - p_0^2\mathbf{M} $$
Où $\mathbf{L}_t$ désigne l'opérateur tangent complexe, $\mathbf{F}_i^{nl}$ est un vecteur second membre qui diffère à chaque ordre 'i' mais qui ne dépend que des variables calculées aux ordres précédents et $(\cdot)^T$ indique l'opérateur transposé.

Les paramètres caractéristiques de cette méthode sont donc l'ordre de troncature des séries (N) et un petit paramètre de tolérance qui permet de vérifier si la méthode a convergé vers la solution. Cette technique de résolution, appelée méthode de Newton d'ordre élevé donne de bons résultats comparés aux résultats trouvés par exemple expérimentalement. L'inconvénient de cette méthode est le temps de calcul, puisque les problèmes à résoudre demandent une triangulation de matrice de grande taille. A noter cependant que la méthode de Newton d'ordre élevé ne demande qu'une seule triangulation de matrice et la résolution de 'N' systèmes linéaires (6) par itération pour une solution initiale (U0, p0). Afin de réduire les temps de calcul, les auteurs proposent dans la référence [8] d'utiliser des matrices réelles au lieu des matrices consistantes complexes, correspondant à la discrétisation de l'opérateur Lt. Cette méthode est plus rapide mais diverge dans quelques cas. L'objectif de ce travail est d'appliquer un modèle de réduction à l'algorithme de Newton d'ordre élevé afin de réduire les temps de calcul tout en conservant les propriétés de convergence de l'algorithme.

### 3. Réduction de modèle

Le modèle de réduction est appliqué à l'équation (1) : le vecteur de déplacement est projeté sur une base de petite dimension:

$$ \mathbf{U} = \mathfrak{R} \mathbf{u} \qquad (7) $$

Où $\mathbf{u}$ est le vecteur réduit et $\mathfrak{R}$, la matrice de projection. Cette matrice est construite en utilisant les vecteurs propres du problème de vibrations linéaires $\Phi$. Elle est, également, enrichie par des vecteurs obtenus en considérant la partie viscoélastique du problème (1), soit d'après la référence [9]:

$$ \mathfrak{R} = [\mathbf{\Phi}, \mathbf{\Psi}] \qquad (8) $$

avec les vecteurs $\mathbf{\Psi}$ solutions des problèmes linéaires suivants :

$$ \mathbf{\Psi} = \mathbf{K}(0)^{-1} \mathbf{K}_v \mathbf{\Phi} \qquad (9) $$

En introduisant la relation (7) dans l'équation (1) et en multipliant à gauche par $\mathfrak{R}^T$, on peut écrire l'équation (1) sous sa forme réduite:

$$ \mathbf{r}(\omega, \mathbf{u}) = [\mathbf{k}(0) + E(\omega)\mathbf{k}_v - \omega^2 \mathbf{m}] [\mathbf{u}] = \mathbf{0} \qquad (10) $$

$$ \text{Où } \mathbf{k}(0) = \mathbf{\Phi}^T \mathbf{K}(0) \mathbf{\Phi}, \quad \mathbf{k}_v = \mathbf{\Phi}^T \mathbf{K}_v \mathbf{\Phi} \text{ et } \mathbf{m} = \mathbf{\Phi}^T \mathbf{M} \mathbf{\Phi} $$

Cette équation réduite est alors résolue en utilisant comme pour le problème de taille réelle, une méthode de Newton d'ordre élevé. Un paramètre de perturbation est alors introduit dans l'équation (10) et les inconnues du problème réduit sont cherchées sous la forme de séries entières tronquées à l'ordre N par exemple.
## 4. Exemples numériques

Pour comparer les algorithmes d'ordre élevé avec et sans méthode de réduction, la méthode des éléments finis est appliquée. L'élément fini utilisé est l'élément coque sandwich à huit nœuds et huit degrés de liberté par nœud développé par Duigou [7] et basé sur l'élément de Büchter et al. [1]. Dans un premier temps, le module d'Young de la poutre est considéré complexe constant et dans un deuxième temps, il est considéré dépendant de la température. Pour s'assurer de la convergence des calculs, un critère basé sur le calcul du résidu est réalisé à chaque pas de la façon suivante :

$$RES \le \delta$$

tel que $\delta$ est la tolérance qui doit être fixée pour chaque méthode et $RES$ est la norme du vecteur résidu calculé pour la méthode de Newton d'ordre élevé à l'aide de l'expression suivante:

$$1. \text{ sans réduction: } RES = \frac{\|\mathbf{R}(\mathbf{U}, p)\|}{\|\mathbf{K}(0)\mathbf{U}\|}$$

$$2. \text{ avec réduction: } RES = \frac{\|\mathbf{r}(\mathbf{u}, p)\|}{\|\mathbf{k}(0)\mathbf{u}\|}$$

### 4.1. Poutre 1 ou exemple avec un module d'Young complexe constant

L'exemple considéré est une poutre sandwich encastrée libre, composée d'une fine couche centrale d'un matériau viscoélastique dont le module d'Young est constant et défini par :

$$E = E_0(1 + i\eta_c)$$

Dans cette étude, le paramètre $\eta_c$ est choisi égal à 1,5 ce qui représente un amortissement relativement élevé. La poutre a les propriétés mécaniques et les caractéristiques géométriques données dans le Tableau 1. Celle-ci est discrétisée avec 424 d.d.l. comme dans la référence [1]. Les paramètres utilisés pour la méthode de Newton d'ordre élevé sont : un ordre de troncature de N=20, une tolérance égale à $1.10^{-7}$. Pour appliquer la réduction, 20 vecteurs sont utilisés. La base (8) est alors composée de 10 vecteurs $\Phi$ (issus d'un calcul de vibrations linéaires) et de 10 vecteurs $\Psi$ solutions du problème linéaire (9). Les résultats obtenus par la méthode de Newton d'ordre élevé avec (désignée par R-HONA) et sans réduction (HONA) sont donnés dans le tableau 2 pour les 2 premiers modes de flexion de la poutre définie dans le Tableau 1.

Le tableau 2 montre que la technique de réduction proposée permet de retrouver les valeurs obtenues à l'aide de la méthode de Newton d'ordre élevé.

<table><caption>Tableau 1. Propriétés mécaniques et géométriques de la poutre 1.</caption><thead><tr><th>Propriétés mécaniques</th><th>Couches élastiques – Aluminium</th><th>Couche viscoélastique</th></tr></thead><tbody><tr><td>Module d'Young</td><td>Ea = 6,9 10<sup>10</sup> N/m<sup>2</sup></td><td>E0 = 1794 10<sup>3</sup> N/m<sup>2</sup></td></tr><tr><td>Coefficient de Poisson</td><td>υa = 0,3</td><td>υp = 0,3</td></tr><tr><td>Masse volumique</td><td>ρa = 2766 kg/m<sup>3</sup></td><td>ρp = 968,1 kg/m<sup>3</sup></td></tr><tr><td>Epaisseur</td><td>ha = 1,524 mm</td><td>hp = 0,172 mm</td></tr></tbody></table>
<table><caption>Tableau 2. Les 2 premiers modes de flexion de la poutre 1. $E=E_0(1+i\eta_c)$= cst, $\eta_c=1.5$ et 20 vecteurs dans la base de réduction.</caption><thead><tr><th>Pulsation non amortie</th><th></th><th>HONA [8]</th><th>R-HONA</th></tr></thead><tbody><tr><td rowspan="4">64.32</td><td>Itération</td><td>2</td><td>3</td></tr><tr><td>Pulsation amortie</td><td>70.19647</td><td>70.19646</td></tr><tr><td>Coeff. d'amortissement</td><td>0.228979</td><td>0.228978</td></tr><tr><td>η<sub>m</sub> / η<sub>c</sub></td><td>0.152653</td><td>0.152652</td></tr><tr><td rowspan="4">297.85</td><td>Itération</td><td>2</td><td>3</td></tr><tr><td>Pulsation amortie</td><td>310.4359</td><td>310.436</td></tr><tr><td>Coeff. d'amortissement</td><td>0.295438</td><td>0.295435</td></tr><tr><td>η<sub>m</sub> / η<sub>c</sub></td><td>0.196959</td><td>0.196957</td></tr></tbody></table>

## 4.2. Poutre 2 ou poutre avec un module d'Young dépendant de la fréquence

Une deuxième poutre sandwich encastrée-libre de dimension (178*10 mm2) est considérée. Mais cette fois, le module d'Young de la couche centrale viscoélastique est variable. Il est représenté par un modèle de Maxwell généralisé :

$$E(\omega) = k_0 + \eta_0 i\omega + \sum_{j=1}^{N_{\max}} \frac{i\omega}{\left(\frac{i\omega}{k_j} + \frac{1}{\eta_j}\right)} \quad (4)$$

Le nombre de Maxwell $N_{\max}$ est égal à 129. Les autres caractéristiques de la poutre sont notées dans le tableau 3.

Pour la simulation élément fini, la température utilisée ici est de 20°C et le nombre de nœuds pour la discrétisation est de 53 soit 424 d.d.l.. Les paramètres utilisés par la méthode de Newton d'ordre élevé avec et sans réduction sont les suivants : un ordre de troncature de N=20, une tolérance de $1.10^{-7}$ et 20 vecteurs dans la base de projection. Les trois premiers modes de flexion sont étudiés.

Dans le tableau 4, on peut noter que le calcul réalisé par la méthode de Newton d'ordre élevé avec ou sans la réduction diverge lorsque l'on prend un ordre de troncature égal à N=20 pour les séries. C'est pourquoi dans le tableau 5, le calcul de fréquences et des amortissements est réalisé en utilisant un ordre de troncature de 5. Dans ce cas, les deux algorithmes donnent de bons résultats.
<table><caption>Tableau 3. Propriétés de la poutre 2.</caption><thead><tr><th>Propriétés mécaniques</th><th>Couche élastique - Acier</th><th>Couche visco-élastique</th></tr></thead><tbody><tr><td>Module d'Young</td><td>E<sub>a</sub> = 2,1 10<sup>11</sup> N/m<sup>2</sup></td><td>E<sub>p</sub> = 27,216 10<sup>6</sup> N/m<sup>2</sup></td></tr><tr><td>Coefficient de Poisson</td><td>υ<sub>a</sub> = 0,3</td><td>υ<sub>p</sub> = 0,44</td></tr><tr><td>Masse volumique</td><td>ρ<sub>a</sub> = 7800 kg/m<sup>3</sup></td><td>ρ<sub>p</sub> = 1200 kg/m<sup>3</sup></td></tr><tr><td>Epaisseur</td><td>h<sub>a</sub> = 0,6 mm</td><td>h<sub>p</sub> = 0,045 mm</td></tr></tbody></table>

<table><caption>Tableau 4. Premiers modes de flexion de la poutre 2. N=20.</caption><thead><tr><th>Pulsation<br>non amortie</th><th></th><th>HONA</th><th>R-HONA</th></tr></thead><tbody><tr><td rowspan="5">29.96</td><td>Iteration</td><td>2</td><td>5</td></tr><tr><td>Pulsation amortie</td><td>32.7725859</td><td>32.7725747</td></tr><tr><td>Coefficient<br>d'amortissement</td><td>0.033747927</td><td>0.033751093</td></tr><tr><td>η<sub>m</sub> / η<sub>c</sub></td><td>43.2200279</td><td>43.2240829</td></tr><tr><td>RES</td><td>3.72228416E-11</td><td>1.14015107E-08</td></tr><tr><td rowspan="5">143.49</td><td>Iteration</td><td>3</td><td>4</td></tr><tr><td>Pulsation amortie</td><td>197.870522</td><td>197.864831</td></tr><tr><td>Damping coefficient</td><td>0.0981405456</td><td>0.0995990691</td></tr><tr><td>η<sub>m</sub> / η<sub>c</sub></td><td>125.68586</td><td>127.553749</td></tr><tr><td>RES</td><td>2.53576098E-12</td><td>3.25600381E-09</td></tr><tr><td>349.78</td><td>Iteration</td><td>Diverge</td><td>Diverge</td></tr></tbody></table>
<table><caption>Tableau 5. Premiers modes de flexion pour la poutre 2. N=5.</caption><thead><tr><th>Mode</th><th></th><th>HONA</th><th>R-HONA</th></tr></thead><tbody><tr><td rowspan="4">29.99</td><td>Iteration</td><td>2</td><td>4</td></tr><tr><td>Pulsation amortie</td><td>32.7842474</td><td>32.785868</td></tr><tr><td>Coefficient d'amortissement</td><td>0.0335581189</td><td>0.0335753581</td></tr><tr><td>η<sub>m</sub> / η<sub>c</sub></td><td>42.9769478</td><td>42.9990273</td></tr><tr><td rowspan="4">143.75</td><td>Iteration</td><td>2</td><td>3</td></tr><tr><td>Pulsation amortie</td><td>198.384908</td><td>198.635266</td></tr><tr><td>Coefficient d'amortissement</td><td>0.0928710797</td><td>0.0931472569</td></tr><tr><td>η<sub>m</sub> / η<sub>c</sub></td><td>118.937401</td><td>119.291092</td></tr><tr><td rowspan="4">349.78</td><td>Iteration</td><td>3</td><td>4</td></tr><tr><td>Pulsation amortie</td><td>548.090895</td><td>548.325729</td></tr><tr><td>Coefficient d'amortissement</td><td>0.0922039246</td><td>0.117323568</td></tr><tr><td>η<sub>m</sub> / η<sub>c</sub></td><td>118.082993</td><td>150.253021</td></tr></tbody></table>

### 3. Conclusion

Ce travail porte sur l'étude des vibrations de structures sandwich viscoélastiques. La couche centrale entraîne la résolution d'un problème aux valeurs propres non-linéaires. Pour cela, deux techniques sont comparées ici : la méthode de Newton d'ordre élevé avec et sans réduction. Deux exemples de poutres sandwich viscoélastiques encastrée-libre sont étudiés. La première a un module d'Young complexe constant ($E= E_0(1+i\eta_c)$) et la seconde un module d'Young qui dépend de la fréquence (modèle de Maxwell pour une température de 20°C). Dans les deux cas, les méthodes proposées donnent de bons résultats. Les tests numériques de ce travail montrent que l'utilisation d'un ordre de troncature proche de 5 conduit à une meilleure convergence que si l'on choisit un ordre élevé. Des tests numériques sur des exemples géométriques plus complexes nécessitant des discrétisations spatiales plus fines, sont actuellement en cours afin de quantifier le gain obtenu avec la méthode de réduction en terme de temps de calculs.
Références

[1] J.R. Banerjee, C.W. Cheung, R. Morishima, M. Perera, J. Njuguna, 2007, Free vibration of a three-layered sandwich beam using the dynamic stiffness method and experiment. International Journal of Solids and Structures 44, 7543–7563, 2007.

[2] N. Büchter, E. Ramm and D. Roehl. Three dimensional extension of non-linear shell formulation based on the enhanced assumed strain concept. International Journal of Numerical Methods in Engineering, 37:2551-2568, 1994.

[3] Yu-Ren Chen, Lien-Wen Chen, 2007, Vibration and stability of rotating polar orthotropic sandwich annular plates with a viscoelastic core layer. Composite Structures 78, 45–57, 2007.

[4] X. Chen, H. L. Chen and X.L. Hu. Damping predication of sandwich structures by order-reduction-iteration approach, Journal of Sound and Vibration 222(5), 803-812, 1999.

[5] A.R. Damanpack, S.M.R. Khalili. High-order free vibration analysis of sandwich beams with a flexible core using dynamic stiffness method. Composite Structures 94, 1503–1514, 2012.

[6] E.M. Daya, M. Potier-Ferry. A numerical method for nonlinear eigenvalue problems application to vibrations of viscoelastic structures. Computers and Structures 79, 533-541, 2001.

[7] L. Duigou. Modélisation numérique de l'amortissement passif des tôles sandwich comportant des couches viscoélastiques ou piézoélectriques. Thèse de Doctorat. Metz. 2002.

[8] L. Duigou, E.M. Daya, M. Potier-Ferry. Iterative algorithms for non-linear eigenvalue problems. Application to vibrations of viscoelastic shells. Comput. Methods Appl. Mech. Engrg. 192, 1323–1335, 2003.

[9] A.M.G. de Lima et al. Component mode synthesis combining robust enriched Ritz approach for viscoelastically damped structures. Engineering Structures 32, 1479–1488, 2010.

[10] C. H. Park, D. J. Inman and M. J. Lam. Model reduction of viscoelastic finite element models, Journal of Sound and Vibration 219(4), 619-637, 1999.