# Simulation numérique de la dynamique non linéaire de micro-structures minces axisymétriques stratifiées piézoélectriques

Arthur Givois, Jean-François Deü, Olivier Thomas

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

## To cite this version:

Arthur Givois, Jean-François Deü, Olivier Thomas. Simulation numérique de la dynamique non linéaire de micro-structures minces axisymétriques stratifiées piézoélectriques. CFM 2017 - 23ème Congrès Français de Mécanique, Aug 2017, Lille, France. hal-03465736

HAL Id: hal-03465736

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

Submitted on 3 Dec 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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 CC BY 4.0 - Attribution - International License
# Étude des vibrations non-linéaires de nano-structures électromécaniques par une démarche de réduction de modèles

A. GIVOIS<sup>a,b</sup>, O THOMAS<sup>a</sup>, J. F. DEÜ<sup>b</sup>

a. Arts et Métiers ParisTech, LSIS UMR CNRS 7296, 8 bd. Louis XIV 59046 Lille, France

b. Laboratoire de Mécanique des structures et des systèmes couplés, Conservatoire National des Arts et Métiers, 2 Rue Conté, 75003 Paris, France

arthur.givois@ensam.eu

## Résumé

On propose dans cette étude une démarche de réduction de modèles pour la prise en compte de non-linéarités géométriques dans les vibrations de structures minces (poutres, plaques, coques) équipées de patchs piézoélectriques. L'application visée est la simulation du comportement vibratoire d'un nano-système électromécanique. Un modèle par éléments finis est employé, et la démarche s'appuie sur une projection des déplacements sur une base de modes linéaires et l'estimation des coefficients non-linéaires par une méthode non-intrusive. Les problèmes statiques et dynamiques non-linéaires sont résolus par une méthode de continuation. Cette démarche est dans un premier temps appliquée à des structures homogènes isotropes et validée par des résultats issus de modèles analytiques, avant d'être étendue à des structures équipées de patchs piézoélectriques.

## Abstract

We propose in this paper a reduced order model procedure to investigate the geometrical nonlinear vibratory response of thin structures (beams, plates and shells) with piezoelectric patches. The main application of this work is to simulate vibrational behavior of an electromechanical nanosystem. This work uses a finite-element formulation, the reduced-order model is based on a projection of displacements on a linear modal basis. Non-linear modal stiffness coefficients are estimated thanks to a non-intrusive method. Static and dynamic problems are solved with continuation method. This procedure is first validated to homogeneous and isotropic structures by comparing with results based on analytical models, and then extended to structures equipped with piezoelectric patches.

**Mots clefs** Vibrations non-linéaires, MEMS, Réduction de modèles, éléments finis, méthodes de continuation.

## Introduction

Le développement des microsystèmes électromécaniques (MEMS) a entraîné l'émergence de nouvelles technologies, parmi lesquelles la mesure de masse de biomolécules à de très hautes résolutions par
l'étude du comportement dynamique de ces capteurs [1]. Si le rapport signal à bruit peut être amélioré par l'augmentation de l'amplitude des vibrations, des phénomènes non-linéaires géométriques deviennent alors non négligeables. Ces phénomènes sont encore mal connus pour de tels systèmes [2] et justifient la nécessité d'établir de nouveaux modèles. Les capteurs MEMS présentent des géométries complexes et variées avec divers types de conditions aux limites (plaques circulaires non-axisymétriques [3], poutres encastrées-libres [4] ou encastrées-encastrées [5]). Dans ce cadre d'étude, les modèles par éléments finis sont privilégiés par rapport à des modèles analytiques, et les méthodes de réduction de modèles par projection des vibrations dans une base de modes linéaires sont bien adaptées [6]. Par ailleurs, les MEMS présentent des facteurs de qualité élevées : la durée élevée des phases transitoires rend ainsi délicate une résolution du système non-linéaire par des méthodes d'intégration temporelle car celles-ci posent régulièrement des problèmes de stabilité et de convergence. Ainsi, une méthode de continuation asymptotique numérique [7] sera préférée pour la résolution des problèmes. Les coefficients non-linéaires du système sont déterminés de coefficients à l'aide de méthodes non-intrusives [8, 9]. L'ensemble de cette démarche a été appliquée sur des MEMS modélisés par des poutres et discrétisés par éléments finis [10], et cette étude propose de l'étendre à d'autres géométries : en particulier, les plaques circulaires sont préférées aux poutre en raison de leur plus grand facteur de qualité [11].

Le formalisme éléments finis avec des non-linéarités géométriques est présenté dans un premier temps. La procédure d'estimation des coefficients non-linéaires fait ensuite l'objet d'une seconde section. Enfin, une démarche de validation systématique sur des structures simplifiées de la méthode est proposée, les coefficients non-linéaires sont pouvant être calculés à partir d'une formulation analytique [12].

Formalisme

Le problème des vibrations non amorties et non-linéaires géométriques d'une structure discrétisée par la méthode des éléments finis (EF) s'écrit :

$$M\ddot{u} + f_{int}(u) = f_e, \quad (1)$$

avec $f_{int}(u) = Ku + f_{nl}$, et où $u, K, M, f_{int}$ et $f_e$ désignent respectivement le vecteur des déplacements nodaux, les matrices de raideur et de masse, les vecteurs des efforts extérieurs et internes et $f_{nl}$ le vecteur des efforts internes dus aux non-linéarités géométriques. Le problème est projeté sur une base modale linéaire : après l'estimation des pulsations propres $\omega_i$ et vecteurs propres $\Phi_i$ par résolution du problème aux valeurs propres

$$(K - \omega_i^2 M)\Phi_i = 0, \quad (2)$$

la matrice modale $\Phi$ est obtenue par l'assemblage des vecteurs propres en colonnes $\Phi = (\Phi_1, \Phi_2, \dots, \Phi_N)$ et les déplacements modaux $q$ sont définis par $u = \Phi q$. Les $\Phi_i$ sont normés de telle sorte que $\Phi_i^T M \Phi_i = 1$ et $\Phi_i^T K \Phi_i = \omega_i^2 \quad \forall i = 1, \dots, N$. Ainsi, l'équation 1 multipliée par $\Phi^T$ s'écrit dans le cadre d'une approximation des champs de déplacement par des fonctions polynomiales ou pour une cinématique de Kirchhoff-Love sous la forme d'un système d'oscillateurs couplés :

$$\ddot{q}_i + \omega_i^2 q_i + \sum_{j=1}^{N} \sum_{k \ge j}^{N} \beta_{jk}^i q_j q_k + \sum_{j=1}^{N} \sum_{k \ge j}^{N} \sum_{l \ge k}^{N} \Gamma_{jkl}^i q_j q_k q_l = f_i, \quad (3)$$

avec $\beta_{jk}^i$, $\Gamma_{jkl}^i$ les coefficients non-linéaires quadratiques et cubiques, tandis que $i$ et $f_i$ désignent res-
pectivement la ième composante modale et l'effort généralisé. Les non-linéarités géométriques sont notamment dus au couplage entre les vibrations transverses et de membrane, ainsi les N modes retenus contiennent à la fois des modes linéaires transverses et des modes linéaires de membranes bien que ces derniers présentent des fréquences plus élevées.

## Estimation des coefficients non-linéaires et résolution

La procédure d'estimation des coefficients non-linéaires consiste à résoudre différents problèmes statiques par un code de calcul par éléments finis, où un déplacement combinaison linéaire des modes propres du système est imposé et depuis lequel les forces nodales sont extraites. Ce champ est choisi de façon à ce qu'un système linéaire d'équations algébriques permette d'obtenir les coefficients quadra-tiques et cubiques. A titre d'illustration, deux champs de déplacement $\mathbf{u} = +\mathbf{\Phi}_1 a$ et $\mathbf{u} = -\mathbf{\Phi}_1 a$ où $a$ est un scalaire permettent d'estimer les coefficients non-linéaires $\beta_{11}^1$ et $\Gamma_{111}^1$. En effet, après avoir noté $\mathbf{f}_e^+ \mathbf{f}_e^-$ les efforts aux noeuds résultants du code de calcul éléments finis, l'équation 3 donne :

$$ 
\begin{align} 
a\omega_1^2 + \beta_{11}^1 a^2 + \Gamma_{111}^1 a^3 &= \mathbf{\Phi}_1^T \mathbf{f}_e^+ \tag{4a} \\ 
-a\omega_1^2 + \beta_{11}^1 a^2 - \Gamma_{111}^1 a^3 &= \mathbf{\Phi}_1^T \mathbf{f}_e^- \tag{4b} 
\end{align} 
$$

qui constituent un système algébrique linéaire. Ici $a$ est choisi volontairement grand, de sorte que le maximum $\mathbf{u} = \pm a \mathbf{\Phi}_1$ soit de l'ordre de 100 fois l'épaisseur de la structure, afin que les termes non-linéaires soient non-négligeables. Les coefficients $\Gamma_{iii}^i$ et $\beta_{ii}^i$ peuvent également être obtenus en imposant deux déplacements proportionnels au mode $\mathbf{\Phi}_i$. Des combinaisons impliquant plusieurs modes doivent être appliquées pour obtenir les autres coefficients.

Une fois que les coefficients non-linéaires sont estimés, les problèmes statique et dynamique sont alors résolus par continuation : le problème est formulé sous la forme

$$ R(U, \lambda) = 0 \tag{5} $$

où U est un vecteur contenant les variables modales $q_i$, $i = 1, ..., N$ et R le système d'équation non-linéaire algébrique dans le cas statique. $\lambda$ désigne le paramètre de chargement extérieur au système ($f_i$ dans le cas de l'équation 3. Le logiciel MANLAB est utilisé pour la résolution d'un tel système d'équations. Dans le cas dynamique, des solutions périodiques du système sont recherchées et MANLAB met en œuvre une procédure de continuation de solutions périodiques.

## Procédure de validation

La méthode est comparée avec des résultats issus de modèles analytiques pour des poutres et des plaques circulaires. Cette comparaison permettra de sélectionner les coefficients non-linéaires non négligeables parmi ceux obtenus après application de la méthode décrite à la section précédente. Les modèles analytiques employés reposent sur les hypothèses suivantes :

- Les mouvements de la poutre et de la plaque circulaire sont respectivement régis par une cinématique d'Euler-Bernoulli et de Kirchhoff-Love.
— Une relation non-linéaire de Von Kàrmàn relie les déplacements aux déformations.

— Les rotations peuvent être linéarisées.

— Les structures ne sont pas précontraintes.

— Les vibrations sont non amorties.

— L'inertie des vibrations de membranes sont négligeables.

— les sections et épaisseurs sont constantes respectivement dans le cas de la poutre et de la plaque, et celles-ci sont constituées d'un matériau homogène isotrope de module de Young E, de coefficient de Poisson ν et de masse volumique ρ.

Ces modèles font intervenir les vibrations de membrane sous la forme d'un effort : il s'agit d'un effort normal N dans le cas d'une poutre et une fonction de force ou d'Airy F pour une plaque. Les formulations liant les déplacements transverses w et les efforts de membrane s'écrivent pour une poutre de direction principale x, de moment quadratique I, de section S et de longueur l :

$$ \rho S \ddot{w} + EI \frac{\partial^4 w}{\partial x^4} - \frac{\partial}{\partial x} \left( N \frac{\partial w}{\partial x} \right) = p(x, t) \quad (6a) $$

$$ N = \frac{ES}{2l} \int_{0}^{l} \left( \frac{\partial w}{\partial x} \right)^2 dx, \quad (6b) $$

et dans le cas d'une plaque d'épaisseur h :

$$ \rho h \ddot{w} + \frac{Eh^3}{12(1 - \nu^2)} \Delta\Delta w - L(w, F) = p(x, t) \quad (7a) $$

$$ \Delta\Delta F = -\frac{Eh}{2} L(w, w), \quad (7b) $$

où Δ désigne l'opérateur laplacien, et L un opérateur bilinéaire défini dans un repère cylindrique par $L(f, g) = f,rr(\frac{g_{rr}}{r} + \frac{g_{r\theta\theta}}{r^2}) + g,rr(\frac{f_{rr}}{r} + \frac{f_{r\theta\theta}}{r^2}) - 2(\frac{f_{rr\theta}}{r} - \frac{f_{r\theta\theta}}{r^2})(\frac{g_{r\theta\theta}}{r} - \frac{g_{\theta\theta}}{r^2})$. En suivant les étapes décrites en [12], les équations précédentes sont adimensionnées, et les déplacements transverses et la fonction d'Airy sont projetés sur une base modale :

$$ w(x, t) = \sum_{k=1}^{K} \Phi_k(x) q_k(t) \quad \text{pour une poutre et} \quad (8a) $$

$$ w(\vec{r}, t) = \sum_{k=1}^{N_w} \Phi_k(\vec{r}) q_k(t) \quad \text{et} \quad F(\vec{r}, t) = \sum_{j=1}^{N_F} \Psi_j(\vec{r}) \eta_j(t) \quad \text{pour une plaque,} \quad (8b) $$

où $\vec{r}$ est la position d'un point du plan médian de la plaque. Les pulsations propres associées sont notées $w_k$ pour les vibrations transverses et $\zeta_j^2$ pour les efforts de membrane. Pour la poutre, les expressions 6 deviennent
$$ \ddot{q}_k + \omega_k^2 q_k + \varepsilon N \sum_{i=1}^{K} E_i^k q_i = Q_k(t) \quad (9a) $$

$$ N = \frac{1}{2} \sum_{i=1}^{K} \sum_{j=1}^{K} H_{ij} q_i q_j, \quad (9b) $$

$$ \text{avec } E_i^k = - \frac{\int_0^1 \Phi_{i,xx} \Phi_k dx}{\int_0^1 \Phi_k^2 dx}, \quad Q_k(t) = \frac{\int_0^1 \Phi_k p(x,t) dx}{\int_0^1 \Phi_k^2 dx}, \quad H_{ij} = \int_0^1 \Phi_{i,x} \Phi_{j,x} dx, $$

et $\varepsilon$ dépendant du choix de l'adimensionnement du déplacement transverse ($\varepsilon = 12$ si $w$ est de l'ordre de l'épaisseur de la poutre). Dans le cas de la plaque, les équations 7 deviennent

$$ \ddot{q}_k + \omega_k^2 q_k - \varepsilon \sum_{p=1}^{N_w} \sum_{j=1}^{N_F} G_{pi}^k q_p \eta_i = K_k(t) \quad (10a) $$

$$ \zeta_j^4 \eta_j = \sum_{p=1}^{N_w} \sum_{q=1}^{K} J_{pq}^j q_p q_q, \quad (10b) $$

$$ \text{avec } G_{pi}^k = \frac{\iint_S \Phi_k L(\Phi_p, \Psi_i) dS}{\int_0^1 \Phi_k^2 dS}, \quad K_k(t) = \frac{\iint_S \Phi_k p(x,t) dx}{\iint_S \Phi_k^2 dx}, \quad J_{pq}^j = \frac{\iint_S \Psi_j L(\Phi_p, \Phi_q) dS}{\iint_S \Psi_j^2 dS}, $$

et $\varepsilon = 12(1 - \nu^2)$ si $w \sim h$. Les coefficients non-linéaires du modèle analytique sont alors calculés à partir des expressions des modes d'une poutre à section constante et d'une plaque circulaire. Les problèmes statiques et dynamique peuvent alors également être résolus avec les méthodes asymptotique numérique, et les résultats sont alors directement comparés avec ceux issus du modèle réduit.

## Conclusion et perspectives

L'article propose de généraliser une méthode de réduction de modèle par éléments finis dans le cadre de l'étude du comportement mécanique des MEMS. Une procédure de validation est en cours pour comparer directement les résultats du modèle réduit avec ceux issus de modèles analytiques. L'intérêt de cette démarche est d'établir un critère pour déterminer quels sont les termes négligeables parmi les coefficients non-linéaires. La modélisation par éléments finis est complexifiée progressivement pour prendre en compte la structure stratifiée puis les patchs piézoélectriques.

## Références

[1] K. Jensen, K. Kim, A. Zettl, An atomic-resolution nanomechanical mass sensor, Nature nanotechnology, 3(9), 2008, 533-537.
[2] O. Thomas, L. Nicu, C. Ayela, C. Touzé, Flambage et vibrations non-linéaires d'une plaque stratifiée piézoélectrique. Application à un bio-capteur MEMS, In 8e Colloque national en calcul de structure, Giens, 2007.

[3] A. K. Ismail, J. S. Burdess, A. J. Harris, C. J. McNeil, J. Hedley, S. C. Chang, G. Suarez, The principle of a MEMS circular diaphragm mass sensor, Journal of micromechanics and microengi- neering, 16(8), 2006, p. 1487.

[4] D. Dezest, O. Thomas, F. Mathieu, L. Mazenq, C. Soyer, J. Costecalde, D. Remiens, J. F. Deü, L. Nicu, Wafer-scale fabrication of self-actuated piezoelectric nanoelectromechanical resonators based on lead zirconate titanate (PZT), Journal of Micromechanics and Microengineering, 25(3), 2015, pp. 035002.

[5] H. Li, S. Preidikman, B. Balachandran C. D. Mote Jr, Nonlinear free and forced oscillations of piezoelectric microresonators, Journal of Micromechanics and Microengineering, 16(2), 2006, p. 356.

[6] A. Nayfeh, M. I. Younis, E. M. Abdel-Rahman, Reduced-order models for MEMS applications, Nonlinear dynamics, 41(1), 2005 pp. 211–236.

[7] B. Cochelin, N. Damil, M. Potier-Ferry, Méthode asymptotique numérique, Hermes Lavoisier, 2007.

[8] A. A. Muravyov, S. A. Rizzi, Determination of nonlinear stiffness with application to random vibra- tion of geometrically nonlinear structures. Computers & Structures, 81 (15), 2003, pp. 1513–1523.

[9] M. P. Mignolet, A. Przekop, S.A. Rizzi, S.M. Spottswood, A review of indirect/non-intrusive re- duced order modeling of nonlinear geometric structures, Journal of Sound and Vibration, 332(10), 2013, pp. 2437–2460.

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

[11] C. Ayela, L. Nicu, Micromachined piezoelectric membranes with high nominal quality factors in newtonian liquid media : A Lamb's model validation at the microscale, Sensors and Actuators B : chemical, 123(2), 2007, pp. 860–868.

[12] C. Touzé, O. Thomas, A. Chaigne, Asymmetric non-linear forced vibrations of free-edge circular plates. Part 1 : Theory, Journal of Sound and Vibration, 258(4), 2002, pp. 649–676.