# Sur le calcul numérique des modes non linéaires

François Blanc, Kerem Ege, Cyril Touzé, Jean-François Mercier, Anne-Sophie

Bonnet-Ben Dhia

► To cite this version:

François Blanc, Kerem Ege, Cyril Touzé, Jean-François Mercier, Anne-Sophie Bonnet-Ben Dhia. Sur le calcul numérique des modes non linéaires. CFM 2011 - 20e Congrès Français de Mécanique, Aug 2011, Besançon, France. pp.1-6. hal-00627519

HAL Id: hal-00627519

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

Submitted on 28 Sep 2011

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
# Sur le calcul numérique des modes non-linéaires

F. Blanc<sup>a</sup>, K. Ege<sup>b</sup>, C. Touzé<sup>a</sup>, J.-F. Mercier<sup>c</sup>, A.-S. Bonnet Ben Dhia<sup>c</sup>

a. UME-ENSTA-ParisTech, Chemin de la Hunière, F-91761 PALAISEAU CEDEX, FRANCE

b. LVA-INSA de Lyon, 25 bis avenue Jean Capelle, F-69621 VILLEURBANNE CEDEX, FRANCE

c. UMA-ENSTA-ParisTech, 32 Boulevard Victor, F-75739 PARIS CEDEX 15, FRANCE

## Résumé :

Nos travaux portent sur le calcul numérique de modes non linéaires. L'approche adoptée consiste à résoudre par différences finies l'équation aux dérivées partielles (EDP) décrivant la géométrie d'un mode non linéaire dans l'espace des phases. Cette EDP est vue comme un problème de transport dont on recherche les conditions initiales donnant des solutions périodiques. Les algorithmes de résolution et d'optimisation sont testés sur un système mécanique à deux degrés de liberté et à non linéarité cubique. Cet exemple nous permet de discuter de la convergence des algorithmes et des problèmes de mise en œuvre. Les résultats sont également comparés à des calculs par continuation.

## Abstract :

This paper deals with non linear normal modes (NNM) computation. The partial differential equations (PDE) governing the geometry of the NNM in phase space is solved via a finite difference scheme. The PDE is integrated as a transport equation, the initial conditions of which are searched so as to result in periodic solutions. The solving and optimisation algorithms are tested on a 2 dof system with cubic non-linearities. This example let us discuss the algorithms' convergence and implementation. The results are also compared to computations with more standard continuation methods.

## 1 Introduction

Les modes non linéaires (MNL) ont été introduits en théorie des vibrations par les travaux de Rosenberg [1], les définissant comme une famille d'orbites périodiques, approche reprise plus tard par Vakakis et al. [2]. En utilisant le théorème de la variété centrale, Shaw et Pierre ont défini un MNL comme une variété invariante de l'espace des phases, tangente à l'origine au sous-espace propre linéaire [3], étendant ainsi le champ d'applications du concept à une classe de systèmes plus large. Depuis, les MNL sont largement utilisés dans le contexte du pompage énergétique [2] et de la réduction de modèles [4]. Cependant, la majorité des méthodes de calcul des MNL repose sur des développements asymptotiques, limitant le champ d'application aux amplitudes modérées, si bien que des techniques de calcul numériques sont aujourd'hui à développer.

Pour ce faire, différentes stratégies ont déjà été proposées : méthodes de tir (Slater [5]), couplées à des méthodes de continuation pour un suivi automatique de branches [6] ; résolution numérique de l'équation de la variété invariante par méthode de Galerkin [7] ; ou encore différences finies sur une formulation du MNL en amplitude-phase [8].

Dans la méthode présentée ici, nous utilisons la définition d'un MNL comme variété invariante. L'équation aux dérivées partielles (EDP) décrivant la géométrie de la variété dans l'espace des phases est résolue par différences finies. Le MNL s'exprime alors comme un problème de transport à une dimension, la variable angulaire jouant le rôle du temps.

La section 2 présente le principe et les équations de la méthode ; un exemple de calcul de MNL d'un
système à deux degrés de liberté est donné en section 3, dont les résultats sont comparés avec des résultats de calculs par continuation (section 4).

## 2 Principe de la méthode

Le principe général de la méthode consiste à calculer une surface dans l'espace des phases ; cette surface est paramétrée par les coordonnées d'un mode maître : le mode linéaire auquel est tangent le MNL au point d'équilibre. En exprimant la dynamique du système en fonction de ce mode maître, on aboutit à deux équations différentielles à résoudre pour chaque degré de liberté du système.

### 2.1 Équations de la variété invariante

Le point de départ des calculs est un ensemble d'oscillateurs non linéaires dans le cas conservatif et sans forçage extérieur. En supposant que le terme linéaire est diagonalisé, les équations représentant la dynamique d'une structure en non linéaire géométrique s'écrivent, $\forall k = 1, ..., N, \ddot{\eta}_k + \omega_k^2 \eta_k = f_k(\eta)$, où $f_k(\eta) = \sum_{i=1}^{N} \sum_{j \ge i}^{N} g_{ij}^k \eta_i \eta_j + \sum_{i=1}^{N} \sum_{j \ge i}^{N} \sum_{\kappa \ge j}^{N} h_{ijk}^k \eta_i \eta_j \eta_\kappa$ regroupe les non-linéarités quadratiques et cubiques.

Les variables $\eta_k(t)$ sont homogènes à un déplacement. Comme l'on a supposé le problème linéaire diagonalisé, $\eta_k$ représente l'amplitude modale du mode propre numéro $k$.

Pescheck et coll. [7] établissent les EDP définissant un MNL (une variété invariante) en coordonnées polaires, choix que nous retenons ici. Soit $p$ le mode maître dont on souhaite calculer le MNL, on pose $\eta_p = a \cos \phi$ et $\dot{\eta}_p = -a\omega_p \sin \phi$ pour passer de la représentation cartésienne usuelle déplacement-vitesse à un problème où l'amplitude et la phase sont nos inconnues. Calculer le $p^{ème}$ mode non linéaire revient à exprimer les $N-1$ autres paires de coordonnées, $\eta_k = P_k(a, \phi)$ et $\dot{\eta}_k = Q_k(a, \phi)$, en fonction du $p^{ème}$ degré de liberté.

Par ailleurs la dynamique du mode maître $p$ fournit les relations $\dot{a} = -\frac{f_p \sin \phi}{\omega_p}$ et $\dot{\phi} = \omega_p - \frac{f_p \cos \phi}{a\omega_p}$.

Au final, on obtient une EDP décrivant la projection du $p^{ème}$ mode non linéaire sur le $k^{ème}$ degré de liberté, $\forall k = 1...N, k \neq p$:

$$ \left\{ \begin{array}{l} \left(1 - \frac{f_p}{\omega_p^2} \frac{\cos \phi}{a}\right) \partial_\phi P_k = \frac{Q_k}{\omega_p} + \frac{f_p}{\omega_p^2} \sin \phi \partial_a P_k \\ \\ \left(1 - \frac{f_p}{\omega_p^2} \frac{\cos \phi}{a}\right) \partial_\phi Q_k = -\frac{\omega_k^2 P_k}{\omega_p} + \frac{f_p}{\omega_p^2} \sin \phi \partial_a Q_k + \frac{f_k}{\omega_p} \end{array} \right. \Leftrightarrow \left\{ \begin{array}{l} \vec{V} \cdot \vec{\nabla} P_k - Q_k = 0 \\ \\ \vec{V} \cdot \vec{\nabla} Q_k + \omega_k^2 P_k = f_k \end{array} \right. \quad (1) $$

### 2.2 Mise en œuvre

Les équations (1) ont la forme d'équations de transport avec termes sources non linéaires, dont la vitesse de transport $\vec{V} = (\dot{a}, \dot{\phi})$ est fonction de $\phi$ et des termes non linéaires. On se propose d'en faire la résolution par différences finies.

#### 2.2.1 Schéma numérique

Pour la résolution par différences finies d'une équation de transport, l'obtention de solutions stables se fait usuellement en décentrant le schéma dans le sens du flot. Ici, le flot est donné par $\dot{a}$ et $\dot{\phi}$. En pratique, $\dot{\phi}$ est positif et un schéma explicite en $\phi$ est utilisé. En revanche, $\dot{a}$ change de signe selon $\phi$; le décentrement du schéma doit suivre ces changements de signe. Il est important de noter que le flot suit la dynamique du système : la figure 1 présente le flot calculé ainsi que différentes trajectoires d'un système à deux degrés de liberté projetées dans le plan du premier degré de liberté, calculées par continuation (voir les sections 3 et 4).
<figure><img src="image_2.png" /><figcaption>Figure 1: Flot de l'équation de transport (flèches bleues) et trajectoires calculées par continuation (---).</figcaption></figure>

### 2.2.2 Conditions aux limites

Le schéma selon $\phi$ étant explicite, des conditions initiales pour $P_k$ et $Q_k$ en $\phi = 0$ suffisent. La résolution selon $a$ est plus délicate car le sens du flot change quand $\phi$ croît. La condition naturelle est en $a = 0$: le système est à l'équilibre, on a donc $P_k(0, \phi) = 0$ et $Q_k(0, \phi) = 0 \ \forall k$. Cette condition est adaptée sur les lignes $\phi = \text{cste}$ pour lesquelles $\dot{a} > 0$ (flot entrant en $a = 0$). La difficulté est lorsque $\dot{a} < 0$ car il faut alors faire entrer dans le domaine de calcul $a \le A$ de l'information non connue venant de $a > A$ (flot entrant en $a = A$). On choisit ici d'extrapoler linéairement les valeurs de $P_k$ et $Q_k$.

### 2.2.3 Optimisation : problème adjoint

Dans la suite de l'article, on se restreint au cas de deux oscillateurs, le mode maître sera désigné par 1, et le mode esclave par 2; on note de plus $P = P_2$ et $Q = Q_2$.

Pour trouver le MNL numériquement, il suffit de déterminer les conditions initiales en $\phi = 0$, notées $P_0 = P(a, 0)$ et $Q_0 = Q(a, 0)$ telles que la solution trouvée soit $2\pi$- périodique, ce qui assure la continuité de la surface calculée. Ainsi, voir ce calcul comme un problème de transport présente des analogies avec les méthodes de tir, à ceci près que nous calculons directement une fonction pour $a \in [0, A]$ au lieu de procéder point par point.

L'approche adoptée est itérative : en partant d'une condition initiale donnée par le calcul asymptotique, conduisant en général à une solution non périodique, il s'agit de l'optimiser afin de réduire l'erreur de périodicité en $\phi$ de la solution calculée. Dans ce but, on se donne le critère $J$, quantifiant l'erreur de périodicité :

$$J(P_0, Q_0) = \frac{1}{2} \int_0^A \left( [P(a, 2\pi) - P(a, 0)]^2 + [Q(a, 2\pi) - Q(a, 0)]^2 \right) da, \quad (2)$$

Pour minimiser $J$, la méthode choisie est une *descente de gradient*, où ce dernier est calculé *via* un problème adjoint. Le gradient est défini à partir des variations de $J$ :

$$J(P_0 + P_0^h, Q_0 + Q_0^h) \approx J(P_0, Q_0) + \int_0^A \nabla_P J(P_0, Q_0) P_0^h da + \int_0^A \nabla_Q J(P_0, Q_0) Q_0^h da.$$

Pour déterminer $\nabla J(P_0, Q_0)$, on utilise la méthode de l'état adjoint : on résout le problème non linéaire direct (équation (1)), ce qui nous donne $P(a, 2\pi)$ et $Q(a, 2\pi)$, puis on résout dans le sens rétrograde le problème adjoint d'inconnues $\tilde{P}$ et $\tilde{Q}$ associé au système (1) linéarisé autour de $(P_0, Q_0)$, avec pour conditions initiales $\tilde{P}(a, 2\pi) = P(a, 2\pi) - P_0(a)$ et $\tilde{Q}(a, 2\pi) = Q(a, 2\pi) - Q_0(a)$. Finalement, le gradient est défini par :

$$\vec{\nabla} J(P_0, Q_0) = \begin{bmatrix} \tilde{P}(a, 0) - \tilde{P}(a, 2\pi) \\ \tilde{Q}(a, 0) - \tilde{Q}(a, 2\pi) \end{bmatrix}.$$
Il faut ainsi résoudre numériquement le problème adjoint rétrograde (en $\phi$), dont les conditions initiales sont données par la résolution du problème direct.

Ce gradient correspond à la correction à apporter aux conditions initiales pour tendre vers des solutions périodiques; en pratique, il convient d'être prudent en appliquant à chaque itération une correction modérée, de la forme $\theta \nabla J(P_0, Q_0)$, où $\theta$ est suffisamment petit.

## 3 Résultats

Cette section présente les résultats de calculs de modes non linéaires pour le système à deux degrés de liberté présenté par Pescheck et al., afin de comparer nos résultats à une solution de référence. Ce système contient les termes non linéaires suivants (voir section 2.1) :

$$f_1 = -0.405\eta_1^3 - 1.34\eta_1^2\eta_2 - 1.51\eta_1\eta_2^2 - 0.349\eta_2^3 \quad (3)$$

$$f_2 = -0.448\eta_1^3 - 1.51\eta_1^2\eta_2 - 1.05\eta_1\eta_2^2 - 4.58\eta_2^3$$

La figure 2 présente les résultats donnés par notre algorithme pour un rayon $a \le A = 2.2$. À cette amplitude, la différence entre les solutions asymptotique et numérique est déjà notable. Sur ces trois représentations (surface $P(a, \phi)$ et coupes en $a = A$ et $\phi = 0$), on observe que la solution numérique s'éloigne de la solution asymptotique et tend vers la solution de référence proposée par Pescheck et al.

<figure><img src="image_3.png" /><figcaption>Figure 2: (a) Nappe représentant le déplacement du second degré de liberté en fonction du premier degré de liberté; (b) coupe $P(A, \phi)$, $A = 2.24$ à différentes itérations; (c) coupe $P(a, 0)$ aux mêmes étapes du calcul.</figcaption></figure>

La figure 3 présente la convergence du calcul. Les normes des gradients selon $P$ et $Q$ décroissent en présentant des oscillations amorties. Il est à noter que les oscillations selon $P$ et $Q$ sont déphasées, ce qui empêche le calcul de converger en bloc. Dans ce cas précis, et en particulier avec le facteur $\theta = 2 \times 10^{-5}$ appliqué au gradient, le calcul converge en quelques 300 itérations.

Le paramètre $\theta$ est déterminant sur la convergence des calculs: il y a une valeur optimale en deçà de laquelle la convergence est lente, et au-delà de laquelle le calcul diverge. Plus $A$ est petit (plus le comportement du système se rapproche d'un comportement linéaire) et plus la valeur de $\theta$ optimale se rapproche de 1, ce qui accélère grandement la convergence du calcul.

Cependant, en considérant $Q(a, 0)$, il apparaît que ce dernier n'est pas identiquement nul (figure 4). Or sait que $Q(a, \phi = 0) = 0$, tous les termes de son développement asymptotique s'annulant. Ce phénomène est lié aux oscillations en $P$ et $Q$ de la correction à apporter aux conditions initiales, et semble se résorber au cours des itérations; cependant, la convergence est d'autant plus lente que $a$ tend vers $A$.

## 4 Comparaison avec des calculs par continuation

Les prédictions données par la méthode de transport sont comparées aux solutions fournies par continuation d'orbites périodiques via la Méthode Asymptotique Numérique (MAN) implémentée dans le logiciel MANLAB [9].
<figure><img src="image_4.png" /><figcaption>Figure 3: (a) Composantes sur P et Q du gradient minimisant le critère d'erreur de périodicité ; (b) erreur de périodicité à différentes étapes de calcul pour P (—) et Q (---).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Malgré des conditions initiales exactes, $Q(a, 0)$ s'éloigne de la solution (ici après 5000 itérations).</figcaption></figure>

La figure 5 (a) présente la coupe dans le plan $(a, P)$, correspondant à l'espace des configurations, des variétés calculées par continuation, transport et calcul asymptotique. On constate immédiatement que la méthode de continuation des orbites périodiques permet de suivre ces dernières au delà des points de rebroussement, comme c'est le cas ici en $a = A_{max} = 2.251$. Par nature, la méthode de la variété invariante utilisée ici comme point de départ pour le calcul par transport ne sera pas en mesure d'aller au delà de ce point. Cependant, même à l'approche de $a = A_{max}$, la convergence devient très lente à mesure que la pente augmente, et on observe que le calcul diverge si l'on souhaite calculer sur un domaine tel que $A > A_{max}$.

<figure><img src="image_6.png" /><figcaption>Figure 5: (a) Coupe dans le plan $\phi = 0$ des variétés prédites par continuation, par transport et par le calcul asymptotique ; les croix indiquent les zones où la variété est instable ; (b) les calculs par continuation montrent que toute la variété ne peut pas être décrite en effectuant le calcul sur une grille régulière.</figcaption></figure>

Enfin, à l'aide des orbites calculées par continuation, on peut reconstruire la variété invariante. La figure 5 (b) montre la continuation des orbites dans le plan $(a, \phi)$. Cette représentation permet de comprendre une différence cruciale entre les approches *continuation* et *variété invariante*, découlant de la définition utilisée pour le calcul des MNL. À cause de la nature de la dynamique, qui rend les orbites fortement elliptiques, tout calcul utilisant une grille régulière dans le plan $(a, \phi)$ ne pourra pas
suivre les déformations imposées par le flot.

# 5 Conclusion et perspectives

Une méthode de calcul numérique des modes non linéaires a été proposée. Elle utilise la définition d'un MNL comme variété invariante de l'espace des phases. Voyant l'EDP gouvernant sa géométrie dans l'espace des phases comme un problème de transport, le problème se ramène à la recherche d'une condition initiale donnant une solution périodique, assurant la continuité de la variété. Le calcul utilise donc un critère de périodicité que l'on optimise *via* un calcul adjoint.

La comparaison avec les méthodes de calcul par continuation d'orbites périodiques montre que ces dernières semblent les plus adaptées au traitement des cas conservatifs, puisqu'elles permettent un meilleur suivi des branches, ainsi que la détection de la stabilité. Cependant, les méthodes de continuation ne s'étendent pas au cas dissipatif, alors que la définition en termes de variété invariante permet de traiter ces cas, plus proches de la réalité.

Les prochaines étapes de ce travail consisteront donc à adapter notre algorithme aux cas dissipatifs et à utiliser les MNL comme base réduite permettant de simplifier les dynamiques forcées de structures minces de type coques.

# References

[1] R. M. Rosenberg. On non-linear vibrations of systems with many degrees of freedom. *Advances in Applied Mechanics*, 9:155-242, 1966.

[2] A.F. Vakakis, O.V. Gendelman, L.A. Bergman, D.M. McFarland, G. Kerschen, and Y.S. Lee. *Nonlinear targeted energy transfer in mechanical and structural systems*. Springer, New-York, 2008.

[3] S. W. Shaw and C. Pierre. Non-linear normal modes and invariant manifolds. *Journal of Sound and Vibration*, 150(1):170-173, 1991.

[4] C. Touzé, M. Amabili, and O. Thomas. Reduced-order models for large-amplitude vibrations of shells including in-plane inertia. *Computer Methods in Applied Mechanics and Engineering*, 197(21-24):2030-2045, 2008.

[5] J.C. Slater. A numerical method for determining nonlinear normal modes. *Nonlinear Dynamics*, 10(1):19-30, 1996.

[6] M. Peeters, R. Viguié, G. Sérandour, G. Kerschen, and J.C. Golinval. Non-linear normal modes, part II: toward a practical computation using numerical continuation techniques. *Mechanical Systems and Signal Processing*, 23(1):195-216, 2009.

[7] E. Pescheck, C. Pierre, and S. Shaw. A new Galerkin-based approach for accurate non-linear normal modes through invariant manifolds. *Journal of Sound and Vibration*, 249(5):971-993, 2002.

[8] D. Noreland, S. Bellizzi, C. Vergez, and R. Bouc. Nonlinear modes of clarinet-like musical instruments. *Journal of Sound and Vibration*, 324(3-5):983-1002, 2009.

[9] B. Cochelin, B. and C. Vergez. A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions. *Journal of sound and vibration*, 324:243-262, 2009.