# Calcul de vibrations non linéaires de micro/nano-structures piézoélectriques stratifiées par modèles réduits avec correction quasi-statique

Olivier Thomas, Saeid Nezamabadi, Jean-François Deü

► To cite this version:

Olivier Thomas, Saeid Nezamabadi, Jean-François Deü. Calcul de vibrations non linéaires de micro/nano-structures piézoélectriques stratifiées par modèles réduits avec correction quasi-statique. 11e Colloque National en Calcul des Structures (CSMA 2013), May 2013, Giens, France. hal-01399615

HAL Id: hal-01399615

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

Submitted on 19 Nov 2016

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 CC BY 4.0 - Attribution - International License
# Calcul de vibrations non linéaires de micro/nano structures piézoélectriques stratifiées par modèles réduits avec correction quasi-statique

Olivier THOMAS<fnref n="1" />, Saeid NEZAMABADI<fnref n="2" />, Jean-François DEÜ<fnref n="3" />

**Résumé** — On s'intéresse dans ce travail au calcul de vibrations de structures avec non linéarités géométriques par la méthode des éléments finis. On cherche à déterminer les régimes vibratoires périodiques du système pour obtenir une cartographie de son comportement en fonction d'un paramètre de contrôle, comme par exemple la fréquence d'excitation, ou pour calculer ses modes non linéaires. Pour réduire les temps de calcul, la taille du modèle est réduite par un approche modale. Dans cet article, on cherche à tester les techniques classiques d'enrichissement de base par ajout de modes quasi-statiques, et notamment à les étendre au cadre non linéaire. Des exemples sur des micro/nano poutres stratifiées piézoélectriques sont proposés.

**Mots clés** — vibrations, modes non linéaire, piézoélectrique, stratification, périodique.

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Exemple de micro-poutre étudiée dans ce travail : photographie prise au microscope électronique à balayage, vue en coupe de l'une des micro-poutres et réponse non linéaire.</figcaption></figure>

L'objectif de ce travail est de fournir un outil de simulation efficace du comportement vibratoire non linéaire de structure minces stratifiées piézoélectriques. On s'intéresse en particulier à des micro/nano-systèmes électromécaniques (M/NEMS en anglais) résonants, potentiellement utilisables comme filtre électronique ultra haute fréquence [1], comme composants logiques élémentaires [2] ou comme capteur de masse moléculaire. Dans ce dernier cas, il s'agit de mesurer avec précision la fréquence de résonance d'un mode donné et de détecter une éventuelle variation de celle-ci lorsque la structure est alourdie par la molécule dont on veut détecter la présence ou mesurer la masse. Certains auteurs sont capables de descendre à des résolutions massiques de l'ordre du zeptogramme ($10^{-21}$ grammes) [3]. Dans chacun des cas, on tire parti de deux propriétés : (i) la petite taille, qui garantie la résolution massique et qui produit des fréquences de résonance allant jusqu'au megahertz et (ii) le très faible amortissement des micro structures mécaniques, qui produit des résonances extrêmement aiguës. Pour coupler les oscillations mécaniques de ces systèmes au circuit électronique environnant (mise en vibration et détection), on utilise des couches de matériaux piézoélectriques.

Le cœur de cette étude est de simuler le comportement vibratoire non linéaire géométrique de ces structures. En effet, d'une part les conditions d'utilisation de ces systèmes imposent des grandes amplitudes de vibrations, suffisantes pour créer des courbes de résonance incurvées (Fig. 1 et [4]). Il s'agit alors de prévoir avec précision le seuil d'apparition des non-linéarités pour prévoir la plage de fonctionnement du système en régime linéaire. D'autre part, il est aussi possible de tirer partie des non-linéarités, notamment à travers une technique d'amplification paramétrique de la structure, qui permet d'augmenter

<footnote n="1">¹ LSIS / Arts et Métiers ParisTech, Lille, olivier.thomas@ensam.eu</footnote>

<footnote n="2">² LMGC / Université de Montpellier 2, saeid.nezamabadi@univ-montp2.fr</footnote>

<footnote n="3">³ LMSSC / Conservatoire National des Arts et Métiers, Paris, jean-francois.deu@cnam.fr</footnote>
artificiellement le facteur de qualité apparent du système et d'améliorer ainsi ses performances [5, 4]. En plus de cela, le modèle proposé prend en compte la géométrie stratifiée des structures, incluant des couches piézoélectriques, ainsi que le couplage électromécanique avec le circuit électronique environnant, pour obtenir des résultats les plus prédictifs possible.

# 1 Formulation non linéaire du problème électro-mécanique

On considère une poutre stratifiée piézoélectrique (Fig. 1) modélisée en utilisant la théorie classique des structures stratifiés, équivalente aux hypothèses d'Euler-Bernouilli [6]. Les non linéarités géométriques sont prisent en compte à partir d'un modèle de von Kármán. L'effet piézoélectrique d'une couche donnée, supposée isotrope transverse, est pris en compte à travers une loi de comportement piézoélectrique linéaire classique. L'état électrostatique des couches piézoélectriques et quantifié par des variables globales : la différence de potentiel entre ses électrodes et la charge électrique, notées respectivement $V^{(k)}$ et $Q^{(k)}$ pour la $k$e couche. L'utilisation d'une procédure standard éléments finis unidimensionnels permet de discrétiser la partie mécanique des équations du mouvement. Elle conduit à la formulation non linéaire suivante, décrivant les mouvements de la poutre piézoélectrique stratifiée [7, 8] :

$$ \left\{ \begin{array}{l} M\ddot{u} + D\dot{u} + Ku + f_{nl}(u) + \sum_{p=1}^{K} \left[ f_c^{(k)} + P^{(k)}u \right] V^{(k)} = f_m, \end{array} \right. \quad (1) $$

$$ C^{(k)}V^{(k)} - Q^{(k)} - \left[ f_c^{(k)} + P^{(k)}u \right]^T u = 0 \quad \forall k = 1, \dots, K. \quad (2) $$

où $u$ est le vecteur de taille $N$ contenant le déplacement axial $u$, le déplacement transverse $w$ et $\theta = \partial w/\partial w$ aux noeuds du maillage. D'autre part, $M, D$ et $K$ sont les matrices de masse, d'amortissement et de rigidité de la partie mécanique, de taille $N \times N$. La partie non linéaire de l'énergie interne résultant des non linéarités géométriques est représentée par un vecteur colonne $f_{nl}(u)$, de dimension $N$. $f_m$ est le vecteur colonne des forces mécaniques extérieures, de taille $N$. Lorsque la structure contient $K$ couches piézoélectriques, $f_c^{(k)}$ et $P^{(k)}u$ sont les vecteurs de couplage électromécanique, de taille $N$, qui couplent les degrés de liberté mécaniques à la différence de potentiel $V^{(k)}$ du $k$e patch. Celui-ci est lié à la charge électrique $Q^{(k)}$ via la capacité électrique $C^{(k)}$ de la $k$e couche piézoélectrique. On peut noter que le terme $P^{(k)}u$ correspond à un forçage paramétrique de la structure. Les définitions de ces matrices et vecteurs sont données dans [8]. Les deux équations précédentes correspondent aux deux effets piézoélectriques : l'effet inverse pour (1) (une tension électrique provoque des déformations électriques) et l'effet direct pour (2) (une déformation provoque un effet électrique).

# 2 Problèmes à résoudre

Dans les applications N/MEMS, on s'intéresse en général à la réponse du système au voisinage d'une de ses résonances $\omega_i$ (avec $\omega_i$ la $i$e pulsation propre du système). Plusieurs problèmes pratiques sont alors à considérer.

- Le cas le plus courant est celui du calcul de la réponse de la structure, en régime permanent, sous l'action d'un forçage sinusoïdal à la fréquence $\Omega$. Celui-ci est soit créé par un forçage mécanique ($f_m = f_0 \cos \Omega t$), soit par une tension électrique imposée sur un patch piézoélectrique ($V^{(k)} = V_0 \cos \Omega t$). On s'intéresse alors au cas $\Omega \simeq \omega_i$.

- Pour certaines applications, le calcul de la réponse résonante de la structure sous excitation paramétrique est nécessaire [2]. Dans ce cas, il faut exciter la structure au double de sa fréquence de résonance : $V^{(k)} = V_p \cos 2\Omega t$, avec $\Omega \simeq \omega_i$.

- Lorsqu'il s'agit d'améliorer les performances du système, la technique d'amplification paramétrique consiste à diminuer artificiellement la largeur de bande des résonances (augmentation artificielle du facteur de qualité) à partir d'une excitation paramétrique de la structure superposée à l'excitation directe : $V^{(k)} = V_0(\cos \Omega t + \varphi) + V_p \cos 2\Omega t$, avec toujours $\Omega \simeq \omega_i$ [4].

- Enfin, un moyen de caractérisation de la réponse non linéaire géométrique d'une structure vibrante donnée est le calcul de ses modes non linéaires, qui permettent d'évaluer l'évolution des
résonances en fonction de l'amplitude des vibrations et ainsi de fournir à moindre coût le squelette de la réponse vibratoire. Une méthode consiste alors à calculer des familles de solutions périodiques du système en vibrations libres, c'est-à-dire de trouver $u(t)$, T-périodique, avec $T = 2\pi/\Omega$ inconnu, solution de l'équation (1) avec $V^{(k)} = 0$ et $f_m = 0$ [9, 10, 11, 12].

- Le calcul de la charge électrique $Q^{(k)}$, pour par exemple se servir de certaine couches piézoélectriques comme capteur, peut se faire en post-traitement à partir de (2), lorsque $u$ est calculé.

Le point commun de cette série de problèmes est la recherche de solutions périodiques de (1). Pour réaliser cela, une méthode consiste à chercher les inconnues du problème sous la forme de séries de Fourier (méthode d'équilibrage harmonique, HBM) et à résoudre ensuite le système algébrique non linéaire obtenu par une méthode de continuation. C'est cette stratégie qui est utilisée ici, avec la Méthode Asymptotique Numérique (MAN) comme méthode de continuation [13].

# 3 Réduction modale et correction quasi-statique

## 3.1 Projection modale

Pour obtenir une formulation réduite du problème discrétisé obtenu à l'équation (1), on décompose le vecteur des inconnues de déplacement mécanique sur la base des modes propres du système en court-circuit, qui sont solutions de l'équation (1) avec $V^{(k)} = 0$ et $f_m = f_{nl}(u) = 0$. Les fréquences propres $\omega_i$ et les déformées modales $\Phi_i$ sont donc les N solutions propres du problème suivant :

$$ (K - \omega^2 M) \Phi = 0, \quad (3) $$

qui ne dépend que des propriétés élastiques du système et est donc facilement calculable [14]. On projette alors le vecteur des déplacements sur les M premier modes propres :

$$ u(t) = \sum_{i=1}^{M} \Phi_i q_i(t), \quad (4) $$

ce qui permet de remplacer l'équation (1) par un système d'oscillateurs couplés avec non linéarités quadratiques et cubiques [8] : $\forall i \in \{1, \dots, M\}$,

$$ \ddot{q}_i + 2\xi_i \dot{q}_i + \omega_i^2 q_i + \sum_{j=1}^{M} \sum_{k \ge j} \beta_{jk}^i q_j q_k + \sum_{j=1}^{M} \sum_{k \ge j} \sum_{l \ge j} \gamma_{jkl}^i q_j q_k q_l + \sum_{k=1}^{K} \chi_i^{(k)} V^{(k)} + \sum_{k=1}^{K} \sum_{j=1}^{M} \Theta_{ij}^{(k)} q_j V^{(k)} = f_i, \quad (5) $$

où

$$ \chi_i^{(k)} = \Phi_i^T f_c^{(k)}, \quad f_i = \Phi_i^T f_m \quad \Theta_{ij}^{(k)} = \Phi_i^T P^{(k)} \Phi_j, \quad \forall i, j \in \{1 \dots M\}, \quad \forall k \in \{1 \dots K\}. \quad (6) $$

et les $\xi_i$ sont les coefficients d'amortissement modaux. Les coefficients non linéaires $\beta_{jk}^i$ et $\gamma_{jkl}^i$ sont calculés avec la méthode expliquée dans [15]. De plus, les déformées modales ont été normées de telle sorte que $\Phi_i^T M \Phi_i = 1$.

Comparé au modèle complet (1), ce système est de taille réduite, puisqu'au lieu des N degrés de liberté éléments-finis, il ne comporte que M oscillateurs non linéaires couplés, avec $M \ll N$. On peut alors résoudre ce système réduit en combinant la méthode d'équilibrage harmonique et la méthode asymptotique numérique, à l'aide du logiciel Manlab (voir [13, 8, 16]). En revanche, du fait de la structure stratifiée des systèmes à l'étude, les déformées modales, qui servent de fonction de projection, sont couplées entre les mouvements axial et transverse, et la convergence est lente [8]. Les sections suivantes explore l'idée d'étendre au cas non linéaire la technique de correction quasi-statique bien connue pour les problèmes de calcul de structures linéaires.

## 3.2 Correction statique du modèle linéaire

On s'intéresse ici au problème du calcul de la réponse dynamique d'une structure linéaire :

$$ M \ddot{u} + D \dot{u} + Ku = f_e, \quad (7) $$
Dans le cas présent des poutres N/MEMS, il s'agit de l'équation (1) linéarisée, où le vecteur de forçage extérieur $f_e$ regroupe le forçage piézoélectrique et le forçage mécanique. L'idée de la correction quasi-statique est de remplacer les $N - M$ modes de haute fréquence négligés dans le développement (4) par leur contribution quasi-statique. Précisément, il s'agit d'écrire :

$$ \mathbf{u}(t) = \sum_{i=1}^{N} \mathbf{\Phi}_i q_i(t) \approx \sum_{i=1}^{M} \mathbf{\Phi}_i q_i(t) + \underbrace{\sum_{i=M+1}^{N} \mathbf{\Phi}_i q_i^{qs}(t)}_{\mathbf{u}^{\text{cor}}(t)}, \quad (8) $$

avec $q_i$ et $q_i^{qs}$ sont définis par :

$$ \left\{ \begin{array}{ll} \ddot{q}_i + \xi_i \omega_i \dot{q}_i + \omega_i^2 q_i = f_i & \forall i \in \{1, \dots, M\}, \\ \omega_i^2 q_i^{qs} = f_i & \forall i \in \{M+1, \dots, N\} \end{array} \right. \quad (9) \quad (10) $$

On notera que l'approximation vient de du fait que la partie dynamique de (10) a été négligée.

En pratique – c'est d'ailleurs ce point qui justifie l'intérêt de cette méthode –, le second terme $\mathbf{u}^{\text{cor}}(t)$ de l'équation (8) peut être calculé en une seule opération, en considérant préalablement $\mathbf{u}^{qs}(t)$ solution du problème linéaire :

$$ K \mathbf{u}^{qs}(t) = f_e(t), \quad (11) $$

En effet, en appliquant le développement modal (8) à $\mathbf{u}^{qs}$, on a :

$$ \mathbf{u}^{qs}(t) = \sum_{i=1}^{N} \mathbf{\Phi}_i q_i^{qs}(t) \quad \Rightarrow \quad \mathbf{u}^{\text{cor}}(t) = \mathbf{u}^{qs}(t) - \sum_{i=1}^{M} \mathbf{\Phi}_i q_i^{qs}(t). \quad (12) $$

Ainsi, en insérant (12) dans (8) et en tenant compte de (10), on montre que la présente méthode de correction quasi-statique consiste à écrire le déplacement $\mathbf{u}(t)$ sous la forme :

$$ \mathbf{u}(t) = \sum_{i=1}^{M} \mathbf{\Phi}_i [q_i(t) - q_i^{qs}(t)] + \mathbf{u}^{qs}(t), \quad \text{avec} \quad q_i^{qs}(t) = f_i(t)/\omega_i^2 = \mathbf{\Phi}_i^{\text{T}} \mathbf{f}_e(t)/\omega_i^2. \quad (13) $$

En pratique, le calcul de $\mathbf{u}^{qs}(t)$ est immédiat dans le cas où le forçage extérieur s'écrit

$$ f_e(t) = f_0 r(t) $$

où $r(t)$ est une fonction scalaire temporelle donnée et le vecteur $f_0$ représente la répartition spatiale du forçage. On a alors :

$$ \mathbf{u}^{qs}(t) = \mathbf{\Psi} r(t), \quad \text{avec} \quad K \mathbf{\Psi} = f_0. \quad (14) $$

Finalement, en définissant le vecteur $\tilde{\mathbf{\Psi}}$ de la façon suivante :

$$ \mathbf{u}^{\text{cor}} = \tilde{\mathbf{\Psi}} r(t), \quad (15) $$

(écriture qui se justifie avec (12), (13) et (14)), on obtient finalement :

$$ \mathbf{u}(t) = \sum_{i=1}^{M} \mathbf{\Phi}_i q_i(t) + \tilde{\mathbf{\Psi}} r(t), \quad \text{avec} \quad \tilde{\mathbf{\Psi}} = \mathbf{\Psi} - \sum_{i=1}^{M} \mathbf{\Phi}_i \mathbf{\Phi}_i^{\text{T}} \mathbf{f}_0 / \omega_i^2. \quad (16) $$

Le modèle réduit ci-dessus est beaucoup plus riche que celui de l'équation (4), il comporte le même nombre de degrés de libertés (les $M$ premières coordonnées modales $q_i(t)$ solutions de (9)), et cela à un prix réduit : celui du calcul de $\mathbf{\Psi}$ par le système linéaire (14) en plus des $M$ premiers modes $\mathbf{\Psi}_i$.

Cette méthode de correction statique est connue depuis les années 1945 [17] et largement utilisée depuis, notamment sous le nom de méthode des accélérations modales [18, 19, 20, 21]. Une idée identique est par ailleurs utilisée dans le cadre d'analyse modales expérimentales, sous la forme de résidus qui sont identifiés en parallèle des paramètres modaux [22, §5.4].
Il est aussi possible d'étendre cette méthode à un cadre plus général, où le vecteur $\tilde{\Psi}$ est calculé a priori avec l'équation (16) et est utilisé comme vecteur de base supplémentaire, sous la forme :

$$ u(t) = \sum_{i=1}^{N} \Phi_{i} q_{i}(t) + \tilde{\Psi} \zeta(t), \qquad (17) $$

où cette fois, $\zeta(t)$ est une fonction inconnue quelconque qui remplace $r(t)$. En constatant que $\tilde{\Psi}$ est par définition orthogonal aux $\Phi_i$ de basse fréquence ($i = 1, \dots, M$) (voir Eq. (12) où $u^{cor}$ se développe sur les $N-M$ modes de haute fréquence), l'équation définissant $\zeta(t)$ est :

$$ \ddot{\zeta} + 2\mu_{cor}\dot{\zeta} + \omega_{cor}^2\zeta = \frac{\tilde{\Psi}^{\mathrm{T}}f_e}{m_{cor}}, \qquad (18) $$

avec $m_{cor} = \tilde{\Psi}^{\mathrm{T}}M\tilde{\Psi}$, $\mu_{cor} = \tilde{\Psi}^{\mathrm{T}}D\tilde{\Psi}/2m_{cor}$ et $\omega_{cor}^2 = \tilde{\Psi}^{\mathrm{T}}K\tilde{\Psi}/m_{cor}$. Cette nouvelle stratégie, appelée méthode d'augmentation de troncature modale, est très légèrement meilleur que la précédente [21, 20], et n'ajoute qu'un seul degré de liberté supplémentaire, $\zeta(t)$, qui se calcule à partir de (18). Bien évidemment, il est aisé de vérifier que si la partie dynamique de (18) est négligée ($\ddot{\zeta} = \dot{\zeta} = 0$), $\zeta(t) = r(t)$ et on retombe sur la méthode précédente.

## 3.3 Correction statique du modèle non linéaire

On s'intéresse ici à étendre les résultats du paragraphe précédent au cas d'une structure avec non linéarités géométriques, dont la dynamique vérifie :

$$ M\ddot{u} + D\dot{u} + Ku + f_{nl}(u) = f_e, \qquad (19) $$

qui correspond au modèle (1) pour une poutre M/NEMS sans forçage paramétrique ($P = 0$). Appliquer la correction quasi-statique introduite au paragraphe précédent consiste là encore à séparer le développement modale entre une partie basse fréquence sur $M$ modes et une partie haute fréquence (les $N-M$ modes restants), à partir de l'équation (8) qui est toujours valable. En revanche, à partir de (5), on montre que les coordonnées modales vérifient :

$$ \left\{ \ddot{q}_i + \xi_i \omega_i \dot{q}_i + \omega_i^2 q_i + \varphi_i(q_j, q_j^{qs}) = f_i \quad \forall i \in \{1, \dots, M\}, \right. \qquad (20) $$

$$ \omega_i^2 q_i^{qs} + \varphi_i(q_j, q_j^{qs}) = f_i \quad \forall i \in \{M + 1, \dots, N\} \qquad (21) $$

où $\varphi_i(q_j, q_j^{qs})$ est une fonction non linéaire quadratique et cubique des ($q_j, q_j^{qs}$), pour tout $j = 1, \dots, N$. Les $\varphi_i$ correspondent au termes facteurs des $\beta_{jk}^i$ et $\gamma_{jkl}^i$ dans (5). La première conséquence des non linéarités est que le système (20) vérifié par les $\tilde{M}$ coordonnées modales $q_j$ de basse fréquence n'est plus découplé de (21), vérifié par les $N-M$ coordonnées modales $q_j^{qs}$ de haute fréquence. Il est ainsi impossible de reprendre telle quelle la méthode des accélérations modales décrite par les équations (13) et (16), car le calcul des $q_i^{qs}(t)$ n'est plus indépendant de celui des $q_i(t)$ dans (21). C'est même précisément cet effet de couplage entre les coordonnées de basse fréquence et celles de hautes fréquence qu'on cherche à capturer.

Néanmoins, il est possible de tenter des méthodes analogue à celle de l'équation (17). On a ici plusieurs choix pour le vecteur $\tilde{\Psi}$ :

- Utiliser le vecteur $\Psi$ solution du problème statique linéaire :

$$ K\Psi = f_0 \qquad (22) $$

- Utiliser un vecteur $\Psi_{nl}$ solution du problème statique non linéaire :

$$ K\Psi_{nl} + f_{nl}(\Psi_{nl}) = f_0. \qquad (23) $$

- Utiliser un vecteur $u_{nl}^{qs}(t)$ solution du problème fréquentiel non linéaire :

$$ K u_{nl}^{qs} + f_{nl}(u_{nl}^{qs}) = f_0 \cos \Omega t, \qquad (24) $$

Comme il n'y a pas de partie dynamique dans l'équation, si $u_{nl}^{qs}(t)$ est décomposé en série de Fourier, la valeur de chacune des composantes de Fourier ne dépend pas de $\Omega$ et un seul calcul non linéaire est nécessaire, quel que soit $\Omega$. Ce sont ces composantes de Fourier qui peuvent être utilisées comme vecteurs d'enrichissement de la base.
Dans tous les cas, il peut être judicieux d'orthogonaliser $\tilde{\Psi}$ avant la projection modale :

$$ \Psi = \sum_{i=1}^{N} \Phi_i q_i^{qs} \quad \Rightarrow \quad \Phi_i^T K \Psi = \omega_i^2 q_i^{qs} \qquad (25) $$

si bien que $\tilde{\Psi}$ est obtenu en enlevant la contribution des $M$ modes de basse fréquence à $\Psi^1$ :

$$ \tilde{\Psi} = \Psi - \sum_{i=1}^{M} \Phi_i \Phi_i^T K \Psi / \omega_i^2. \qquad (26) $$

Ainsi, en introduisant l'équation (17) dans (19) et en orthogonalisant le résultat par la base $\{\Phi_1, \dots, \Phi_M \tilde{\Psi}\}$, on obtient le problème réduit non linéaire vérifié par les coordonnées modales $\{q_1, \dots, q_M \zeta\}$ :

$$ \left\{ \ddot{q}_i + \xi_i \omega_i \dot{q}_i + \omega_i^2 q_i + \varphi_i(q_i, \zeta) = f_i \quad \forall i \in \{1, \dots, M\}, \right. \qquad (27) $$

$$ \left. \ddot{\zeta} + 2\mu_{cor}\dot{\zeta} + \omega_{cor}^2\zeta + \varphi_{cor}(q_i, \zeta) = f_{cor} \right. \qquad (28) $$

avec $\varphi_i = \Phi_i^T f_{nl}(u)$, $\varphi_{cor} = \tilde{\Psi}^T f_{nl}(u)/m_{cor}$, $f_{cor} = \tilde{\Psi}^T f_e/m_{cor}$, $m_{cor} = \tilde{\Psi}^T M \tilde{\Psi}$, $\mu_{cor} = \tilde{\Psi}^T D \tilde{\Psi}/2m_{cor}$ et $\omega_{cor}^2 = \tilde{\Psi}^T K \tilde{\Psi}/m_{cor}$. Là encore, ce problème réduit est de taille voisine de son équivalent sans correction statique ($M + 1$ degrés de liberté au lieu de $M$).

# 4 Calcul de la charge électrique

Dans ce paragraphe, on s'intéresse au calcul de la charge électrique $Q(t)$ générée par une couche piézoélectrique durant les vibrations de la poutre. On considère uniquement des vibrations linéaires et on cherche à tester le bénéfice d'une correction statique du modèle à $M < N$ modes (Eq. (16)). Pour les calculs, on considère une poutre encastrée/libre similaire à celles représentées sur la figure 1 (de longueur $l_b = 700 \ \mu m$), avec un patch piézoélectrique du coté encastré, de longueur $l_p = l_b/4 = 175 \ \mu m$. La couche piézoélectrique est utilisée à la fois comme actionneur, en imposant la différence de potentiel $V(t) = V_0 \cos \Omega t$, et comme capteur, en mesurant la charge générée $Q(t)$.

On utilise les équations (1) et (2), avec $K = 1$, $f_{nl} = f_m = 0$ et $P^{(k)} = 0$, qui prend la forme de l'équation (7) avec $f_e = -f_c V(t)$ et $f_c = f_e^{(1)}$. Le déplacement $u(t)$ est alors solution du problème suivant :

$$ u(t) = \sum_{i=1}^{M} \Phi_i q_i(t) + \tilde{\Psi} V(t), \qquad (29) $$

avec

$$ \ddot{q}_i + 2\xi_i \omega_i \dot{q}_i + \omega_i^2 q_i = -\chi_i V, \quad \tilde{\Psi} = \Psi + \sum_{i=1}^{M} \Phi_i \Phi_i^T f_c / \omega_i^2, \quad K \Psi = -f_c \qquad (30) $$

Pour le calcul de la charge, on utilise en post-traitement :

$$ Q(t) = CV(t) - f_c^T u(t). \qquad (31) $$

En pratique, comme on cherche à évaluer les déplacements de la poutre $u(t)$ à partir de la mesure de charge $Q(t)$, le terme $CV(t)$ de l'équation précédente est annulé par un dispositif électronique [23]. Dans la suite, on considère simplement $Q(t) = -f_c^T u(t)$.

La figure 2 montre que la correction statique n'a d'influence sur le déplacement $u(t)$ qu'en dehors des résonances. Elle est notamment extrêmement efficace pour prédire les anti-résonances. À l'opposé, si l'on ne s'intéresse qu'aux résonances, un modèle à un seul degré de liberté sans correction statique est suffisant. En revanche, dans le cas du calcul de la charge électrique $Q(t)$, on montre qu'un grand nombre de modes est nécessaire pour obtenir une bonne précision, même aux résonances, et qu'une correction statique offre à moindre coût une excellente approximation. Ceci peut s'expliquer en constatant que l'action piézoélectrique imposé par la tension $V(t)$ est de moduler la longueur de la couche piézoélectrique.

<footnote n="1">1. On peut noter que l'équation (26) montre aussi que dans le cas linéaire, la procédure qui à conduit à utiliser $\tilde{\Psi}$ au lieu de $\Psi$ dans (16) est en fait une orthogonalisation.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 2 – Convergence modale, en fonction du nombre *M* de modes retenus dans le développement modal, pour le déplacement transverse à l'extrémité du patch piézoélectrique (gauche) et pour la charge électrique (milieu, droite). ('—') : modèle avec correction statique ; ('- -') modèle sans correction statique. Poutre encastrée-libre de longueur *l*<sub>p</sub> = 700 µm, oscillant sous une tension d'amplitude *V*<sub>0</sub> = 1 V, avec un amortissement ξ<sub>1</sub> = 0.01, autour de la première résonance. Calcul avec 100 éléments-finis.</figcaption></figure>

Cela crée par conséquent un mouvement de flexion de la poutre, mais aussi des déformations longitudinales. La charge électrique $Q(t)$, de manière identique, résulte à la fois du mouvement de flexion et du mouvement longitudinal de la poutre. L'effet de ce dernier, non négligeable sur la charge, n'est pas capturé par les premiers modes $\Psi_i$ du développement modal, car ce sont des modes de flexion. Il est donc nécessaire d'inclure dans la base de projection des modes axiaux, qui arrivent dans le spectre à des fréquences bien plus hautes, et dont l'effet est bien capturé par la correction quasi-statique. Cela se vérifie sur le diagramme d'erreur de la Fig. 2 (droite) : ce sont deux modes axiaux (de rang 18 et 27) qui sont responsable des sauts qui améliore la convergence du modèle.

L'expression analytique suivante, tirée de [24] :

$$Q = -\Xi u(x_+) + \Theta w_{,x}(x_+) + CV \quad (32)$$

confirme l'analyse : la charge électrique est proportionnelle au déplacement axial $u$ à l'extrémité du patch piézoélectrique en $x = x_+$ et au déplacement transverse $w_{,x}$ au même point.

## 5 Conclusion

En linéaire, on a mis en évidence l'effet sur la charge électrique du couplage piézoélectrique avec le mouvement axial de la poutre, qui est avantageusement calculé avec une correction statique. En non linéaire, seules des simulations numériques du problème (27-28), en court, permettront d'évaluer la performance d'une correction statique.

## Références

[1] V. Gouttenoire, T. Barois, S. Perisanu, J.-L. Leclercq, S. T. Purcell, P. Vincent, and A. Ayari. Digital and fm demodulation of a doubly clamped single-walled carbon-nanotube oscillator : Towards a nanotube cell phone. Small, 6(9):1060-1065, 2010.

[2] I. Mahboob and H. Yamaguchi. Bit storage and bit flip operations in an electromechanical oscillator. Nature nanotechnology, 3:275-279, 2008.

[3] Y. T. Yang, C. Callegari, X. L. Feng, K. L. Ekinci, and M. L. Roukes. Zeptogram-scale nanomechanical mass sensing. Nano Letters, 6(4):583-586, 2006.

[4] O. Thomas, F. Mathieu, W. Mansfield, C. Huang, S. Trolier-McKinstry, and L. Nicu. Efficient parametric amplification in mems with integrated piezoelectric actuation and sensing capabilities. Applied Physics Letters, 2013. Submitted.

[5] O. Thomas, S. Nezamabadi, F. Mathieu, L. Nicu, and J.-F. Deü. Finite elements reduced order models of piezoelectric micro-resonators : electric charge computation and nonlinear parametric amplification. In Actes du XVIIIe. congrès "Vibrations, Chocs & Bruit", Clamart, France, July 2012.
[6] J. N. Reddy. *Mechanics of laminated composite plates*. CRC press, 1997.

[7] O. Thomas, J.-F. Deü, and J. Ducarne. Dynamics of an elastic structure with piezoelectric patches: finite-element formulation and electromechanical coupling coefficients. *International Journal of Numerical Methods in Engineering*, 80(2):235-268, 2009.

[8] A. Lazarus, O. Thomas, and J.-F. Deü. Finite elements reduced order models for nonlinear vibrations of piezoelectric layered beams with applications to NEMS. *Finite Elements in Analysis and Design*, 49(1):35-51, 2012.

[9] C. Touzé, O. Thomas, and A. Chaigne. Hardening/softening behaviour in non-linear oscillations of structural systems using non-linear normal modes. *Journal of Sound Vibration*, 273(1-2):77-101, 2004.

[10] R. Arquier, S. Bellizzi, R. Bouc, and B. Cochelin. Two methods for the computation of nonlinear modes of vibrating systems at large amplitude. *Computers & Structures*, 84:1565-1576, 2006.

[11] G. Kerschen, M. Peeters, J. C. Golinval, and A. F. Vakakis. Nonlinear normal modes, part I: A useful framework for the structural dynamicist. *Mechanical Systems and Signal Processing*, 23(1):170-194, 2009.

[12] O. Thomas. *Dynamique linéaire et non linéaire de structures élastiques et piézoélectriques*. Mémoire d'Habilitation à diriger des recherches, École Normale Supérieure de Cachan, November 2011.

[13] B. Cochelin and C. Vergez. A high order purely frequential harmonic balance formulation. *Journal of Sound and Vibration*, 324(1-2):243-262, 2009.

[14] A. Sénéchal, O. Thomas, and J.-F. Deü. Optimization of shunted piezoelectric patches for vibration reduction of complex structures - application to a turbojet fan blade. In *Proceedings of the ASME 2010 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference, IDETC/CIE 2010*, Montreal, Canada, August 2010.

[15] A. A. Muravyov and S. A. Rizzi. Determination of nonlinear stiffness with application to random vibration of geometrically nonlinear structures. *Computers and Structures*, 81(15):1513-1523, 2003.

[16] R. Arquier, S. Karkar, A. Lazarus, O. Thomas, C. Vergez, and B. Cochelin. Manlab 2.0: an interactive pathfollowing and bifurcation analysis software. Technical report, Laboratoire de Mécanique et d'Acoustique, CNRS, http://manlab.lma.cnrs-mrs.fr, 2005-2011.

[17] D. E. Williams. Dynamic loads in aeroolanes under given impulsive loads with particular reference to landing and gust loads on a large flying boat. *Great Britain Royal Aircraft Establishment, SME Reports* 3309 and 3316, 1945.

[18] É. Balmès. Methods for vibration design and validation - méthodes de conception et de validation en vibration. http://www.sdtools.com, 2010.

[19] M. Géradin and D. Rixen. *Théorie des vibrations*. Masson, 2nd. edition, 1996.

[20] D. J. Rixen. Generalized mode acceleration methods and modal truncation augmentation. In *42nd AIAA/ASME/ASCE/AHS/ASC structures, structural dynamics, and materials conference and exhibit proc.*, Seattle, USA, April 2001.

[21] J. M. Dickens, J.M. Nakagawa, and M.J. Wittbrodt. A critique of mode acceleration and modal truncation augmentation methods for modal response analysis. *Computers & Structures*, 62(6):985-998, 1997.

[22] D. J. Ewins. *Modal testing: theory and practice*. Research studies press ltd., 1984.

[23] F. Mathieu, F. Larramendy, D. Dezest, C. Huang, G. Lavallee, S. Miller, C. M. Eichfeld, W. Mansfield, S. Trolier-McKinstry, and L. Nicu. *J. Micromech. Microeng*, 2012, submitted.

[24] J. Ducarne, O. Thomas, and J.-F. Deü. Placement and dimension optimization of shunted piezoelectric patches for vibration reduction. *Journal of Sound and Vibration*, 331(14):3286-3303, 2012.