# Réduction de modèles dynamiques non-linéaires en grands déplacements

Youssef Gerges, Emeline Sadoulet-Reboul, Mohamed Guedri, Morvan Ouisse, Noureddine Bouhaddi

► To cite this version:

Youssef Gerges, Emeline Sadoulet-Reboul, Mohamed Guedri, Morvan Ouisse, Noureddine Bouhaddi. Réduction de modèles dynamiques non-linéaires en grands déplacements. 20e Congrès Français de Mécanique, Aug 2011, besançon, France. hal-01512524

HAL Id: hal-01512524

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

Submitted on 23 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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Réduction de modèles dynamiques non-linéaires en grands déplacements - Application aux plaques minces

Y. GERGESª, E. SADOULET-REBOULª, M. GUEDRIᵇ, M. OUISSEª, N. BOUHADDIª

a. Institut FEMTO ST, Départ. Méc'Appli, 24 rue de l'Epitaphe, 25030 BESANCON

b. Institut Préparatoire aux Etudes d'Ingénieurs de Nabeul, 8000 M'rezgua, NABEUL, TUNISIE

## Résumé :

Cet article présente une méthode de réduction adaptée aux problèmes de dynamique non-linéaire. La non-linéarité étudiée est de type géométrique en grands déplacements. La réduction est basée sur la méthode des approximations combinées qui consiste à construire une base de série binomiale actualisée. Cette base est issue du comportement non-linéaire considéré comme une modification du modèle linéaire. La simulation numérique de la réponse temporelle d'une plaque mince modélisée en éléments finis (MEF) montre l'intérêt de la méthode proposée.

## Abstract :

This paper presents a reduced order method adapted to a nonlinear dynamic problem. Geometrical non-linearity with large displacements behaviour is considered. The reduction tool is based on the combined approximation method that consists in expanding a binomial series. This basis is issued from the nonlinear behaviour considered as a modification of the linear model. The numerical simulation of the time response of a thin plate modeled by the finite element method (FEM) shows the applicability of the proposed method.

**Mots clefs :** vibration non-linéaire ; réduction de modèle ; intégration temporelle

## 1 Introduction

La formulation et la résolution numérique d'un problème de dynamique non-linéaire est un problème complexe. Il dépend du type de non-linéarité étudié ; dans notre cas, on s'intéresse aux non-linéarités géométriques de type grands déplacements. Un des problèmes majeurs est le coût de calcul qui, comparé à un problème linéaire, peut s'avérer prohibitif. La réduction de modèle est un outil envisageable pour tenter de réduire ce temps de calcul. Des outils bien spécifiques, [3, 7] basés sur des méthodes dédiées à la dynamique non-linéaire comme les modes normaux non-linéaire (NNM) ou la décomposition propre orthogonale (POD) sont proposés dans la littérature. Réduire le modèle avec une base de Ritz est une solution concevable à condition de prendre en compte le comportement non-linéaire du modèle qui nécessite une actualisation de la base de réduction suivant l'évolution de la non-linéarité. Dans ce papier le comportement non-linéaire est considéré comme une modification du comportement linéaire et selon cette modification, une base de réduction est construite. Après un rappel de la formulation du problème non-linéaire en grands déplacements, la méthode de réduction est présentée. La simulation de la réponse temporelle d'une plaque encastrée illustre les avantages et les inconvénients de la méthode.

## 2 Formulation du problème non-linéaire - Application aux plaques minces

Dans ce paragraphe, on rappelle la formulation d'un problème de dynamique non-linéaire dans le cas de structures minces. On peut trouver la formulation détaillée dans plusieurs ouvrages notamment [6, 9].
On considère une plaque de surface moyenne $\Gamma_m$ et d'épaisseur $h$ dans le repère $\mathbb{R}(x, y, z)$. On se place dans l'hypothèse de contraintes planes et on utilise le modèle de Love-Kirchoff dédié aux plaques minces. Les déplacements sont ainsi déterminés par la connaissance du mouvement de la surface moyenne dont les composantes sont notées $(u, v, w)$.

Le champ de déplacements dans le repère global s'écrit sous la forme :

$$ \underline{u}(x, y, z, t) = \begin{cases} u(x, y, t) \\ v(x, y, t) \\ w(x, y, t) \end{cases} - z \begin{cases} w_{,x}(x, y, t) \\ w_{,y}(x, y, t) \\ 0 \end{cases}. \qquad (1) $$

Le tenseur du 2<sup>nd</sup> ordre des déformations de Green-Lagrange est lié au champ de déplacements par la relation suivante :

$$ E(\underline{u}) = \frac{1}{2}(\nabla \underline{u} + \nabla \underline{u}^T) + \frac{1}{2}\nabla \underline{u}^T \nabla \underline{u}. \qquad (2) $$

L'hypothèse des petites déformations et des rotations modérées permet d'écrire le tenseur des déformations de Green-Lagrange sous forme vectorielle :

$$ E = \begin{Bmatrix} E_{xx} \\ E_{yy} \\ 2E_{xy} \end{Bmatrix} = \underbrace{\begin{Bmatrix} u_{,x} \\ v_{,y} \\ u_{,y} + v_{,x} \end{Bmatrix} + \begin{Bmatrix} \frac{1}{2}(w_{,x})^2 \\ \frac{1}{2}(w_{,y})^2 \\ w_{,x} w_{,y} \end{Bmatrix}}_{E^m} - Z \underbrace{\begin{Bmatrix} w_{,xx} \\ w_{,yy} \\ 2w_{,xy} \end{Bmatrix}}_{K^f}, \qquad (3) $$

où $E^m$ représente les déformations dues à l'effet membrane et $E^f = ZK^f$ les déformations dues à l'effet flexion.

Soient $U$ l'énergie potentielle totale égale à la somme de l'énergie de déformation du corps et de l'énergie potentielle des forces exercées et $T$ l'énergie cinétique, d'après le principe des travaux virtuels on a :

$$ \delta \int_{t_1}^{t_2} (U + T) dt = 0. \qquad (4) $$

Le développement de l'équation (4) en introduisant $S$, 2<sup>nd</sup> tenseur des contraintes de Piola-Kirchoff, ainsi que $D^m$ et $D^f$, matrices de raideur associées respectivement aux effets de membrane et de flexion s'écrit :

$$ \int_{\Gamma_m} (\delta E^m)^T D^m E^m dS + \int_{\Gamma_m} (\delta K^f)^T D^f K^f dS + \delta \Pi_{\text{inertie}} - \delta \Pi_{\text{ext}} = 0, \qquad (5) $$

où

$$ \left\{ \begin{array}{l} D^m = \frac{Eh}{1 - \nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{bmatrix} \\ D^f = \frac{h^2}{12} D^m \end{array} \right. \qquad (6) $$

avec

$$ \left\{ \begin{array}{l} \delta\Pi_{\text{inertie}} = \int_{\Omega} \rho \ddot{u} \delta u d\Omega \quad \text{puissance virtuelle associée aux effets inertiels,} \\ \delta\Pi_{\text{ext}} = \int_{\Omega} F \delta u d\Omega \quad \text{puissance virtuelle associée aux efforts extérieurs.} \end{array} \right. \qquad (7) $$

**Discrétisation par la méthode des éléments finis** : le problème est discrétisé par éléments finis en introduisant le développement suivant :

$$ \left\{ \begin{matrix} u \\ v \\ w \end{matrix} \right\} = \begin{bmatrix} N & 0 & 0 & 0 & 0 \\ 0 & N & 0 & 0 & 0 \\ 0 & 0 & N^w & N^{\theta_x} & N^{\theta_y} \end{bmatrix} \underbrace{\left\{ \begin{matrix} \tilde{u} \\ \tilde{v} \\ \tilde{w} \\ \tilde{\theta}_x \\ \tilde{\theta}_y \end{matrix} \right\}}_{\tilde{a}}, \qquad (8) $$
où $N, N^w, N^{\theta_x}$ et $N^{\theta_y}$ désignent les fonctions d'interpolation.

La discrétisation du tenseur des déformations de Green-Lagrange conduit à :

$$ E^m = \underbrace{\begin{bmatrix} N, x & 0 \\ 0 & N, y \\ N, y & N, x \end{bmatrix}}_{B^m} \underbrace{\left\{ \begin{matrix} \tilde{u} \\ \tilde{v} \end{matrix} \right\}}_{\mathbf{u}} + \frac{1}{2} \underbrace{\left( \begin{bmatrix} w, x & 0 \\ 0 & w, y \\ w, y & w, x \end{bmatrix} \right)}_{B^{nl}(\mathbf{w})} \underbrace{\left\{ \begin{bmatrix} N, w & N, \theta_x & N, \theta_y \\ N, y & N, \theta_x & N, \theta_y \\ \tilde{w} & \tilde{\theta}_x & \tilde{\theta}_y \end{bmatrix} \right\}}_{\mathbf{w}} \quad (9) $$

$$ K^f = \underbrace{\begin{bmatrix} N_{,xx} & N_{,\theta_x} & N_{,\theta_y} \\ N_{,yy} & N_{,\theta_x} & N_{,\theta_y} \\ N_{,xy} & N_{,\theta_x} & N_{,\theta_y} \end{bmatrix}}_{B^f} \left\{ \begin{matrix} \tilde{w} \\ \tilde{\theta}_x \\ \tilde{\theta}_y \end{matrix} \right\}. \quad (10) $$

Le développement de la puissance virtuelle de déformation conduit à :

$$ \delta\Pi_{def} = \delta\tilde{a}^T K \tilde{a} = \delta\tilde{a}^T \int_{\Gamma_m} \left[ B^{mT} D^m B^m \quad \frac{1}{2} B^{mT} D^m B^{nl}(\mathbf{w}) \quad B^{nlT}(\mathbf{w}) D^m B^m \quad \frac{1}{2} B^{nlT}(\mathbf{w}) D^m B^{nl}(\mathbf{w}) + B^f T D^f B^f \right] \tilde{a} dS. \quad (11) $$

Finalement, la discrétisation de la formulation variationnelle conduit à l'équation générale suivante :

$$ M \ddot{a} + K(\tilde{a}) \tilde{a} - F(t) = 0. \quad (12) $$

L'équation (12) est l'équation d'équilibre dynamique du système conservatif avec $M$ matrice de masse et $K$ matrice de rigidité non-linéaire. La méthode de résolution d'un tel système dépend de la nature de $F(t)$. S'il s'agit d'une excitation harmonique, la résolution peut alors être effectuée par la méthode de la balance harmonique. S'il s'agit d'une excitation quelconque, la résolution de (12) se fait par intégration temporelle.

Dans les deux cas, cette résolution s'avère lourde et coûteuse, en particulier pour les MEF de grande taille. La réduction de modèle adaptée aux problèmes non-linéaires permet de réduire le temps de calcul. Ce modèle réduit peut par la suite être utilisé pour des procédures itératives de type optimisation, analyse de robustesse et fiabilité.

## 3 Réduction de modèle par une variante de la méthode des approximations combinées

La méthode des approximations combinées a été développée à l'origine par Kirsch [4] pour des stratégies de réanalyse nécessaires en études d'optimisation et de robustesse. Elle consiste à exprimer chaque vecteur propre du problème initial sous la forme d'une série binomiale fonction de la modification apportée au problème. Les sous-bases développées sont ensuite réunies afin de construire une base de réduction globale. Pour les analyses non-linéaires, le comportement non-linéaire est considéré comme une modification du comportement linéaire, ce qui permet de construire la base de réduction [2, 5]. D'autre part, le comportement non-linéaire étant évolutif, la base de réduction doit l'être aussi. La méthode des approximations combinées est présentée dans la suite.

Soit une structure quelconque modélisée par les matrices masse et raideur, notées respectivement $M_0$ et $K_0$. Le problème aux valeurs propres linéaire de cette structure vérifie l'équation suivante :

$$ K_0 r_0^\nu = \lambda_0^\nu M_0 r_0^\nu, \quad (13) $$

où $r_0^\nu$ et $\lambda_0^\nu$ sont respectivement le vecteur propre et la valeur propre du $\nu$-ième mode.
On apporte une modification à cette structure donnant lieu à des variations de masse $\Delta M$ et de raideur $\Delta K$. Le nouveau problème aux valeurs propres s'écrit sous la forme :

$$ (K_0 + \Delta K) r^{\nu} = \lambda^{\nu} (M_0 + \Delta M) r^{\nu}. \quad (14) $$

$r^{\nu}$ peut être approximé par un développement en série binomiale qui dépend de la structure initiale et de la modification apportée. Soit $n$ le nombre de modes de la base modale tronquée, soit $s$ l'ordre du développement souhaité dans la série binomiale; pour un mode $\nu$ donné, la construction de la sous-base de réduction se présente comme suit :

$$ \begin{cases} r_1^{\nu} = K_0^{-1} (M_0 + \Delta M) r_0^{\nu} \\ r_i^{\nu} = -Br_{i-1}^{\nu} \quad \text{avec } (i = 2, 3, \dots, s) \end{cases} \quad (15) $$

avec

$$ B = K_0^{-1} \Delta K. \quad (16) $$

Pour les analyses non-linéaires, $\Delta K$ est la variation de la raideur non-linéaire par rapport à la raideur linéaire. On peut l'écrire sous la forme :

$$ \Delta K = K(\tilde{a}) - K_0. \quad (17) $$

On construit ainsi une sous-base $r_B^{\nu}$ à partir de chaque vecteur $r_0^{\nu}$ tel que :

$$ r_B^{\nu} = [r_1^{\nu} \ r_2^{\nu} \ \dots \ r_i^{\nu} \ \dots \ r_s^{\nu}]. \quad (18) $$

On assemble ensuite toutes les sous-bases pour former une base globale :

$$ r_B = [r_B^1 \ r_B^2 \ \dots \ r_B^{\nu} \ \dots \ r_B^n]. \quad (19) $$

La propriété d'orthogonalité entre les vecteurs n'est pas prise en compte dans la construction de la base, ceci peut engendrer un mauvais conditionnement des matrices réduites qui peut fausser le modèle prédictif. Pour assurer le bon conditionnement du modèle, on effectue une normalisation de ces vecteurs par rapport à la matrice de masse; une décomposition en valeurs singulières (SVD) de la base de réduction est aussi nécessaire si le rang de cette base est inférieur à sa taille.

Le problème réduit s'écrit alors sous la forme :

$$ \bar{M}\ddot{q} + \bar{f}_{nl}(\tilde{a}) = \bar{F}(t) \quad (20) $$

avec :

$$ \left\{ \begin{array}{l} \bar{M} = r_B^T M r_B \\ \bar{f}_{nl}(\tilde{a}) = r_B^T f_{nl}(\tilde{a}) \\ \bar{F}(t) = r_B^T F(t). \end{array} \right. \quad (21) $$

# 4 Application

On s'intéresse à la réponse temporelle d'une plaque encastrée sur les bords et excitée en son centre par une force harmonique. La résolution s'effectue par intégration implicite en utilisant le schéma de prédiction-correction de Newmark [1]. Un pas de temps adaptatif est utilisé pour optimiser le temps de calcul [8].

L'exemple académique proposé permet d'illustrer l'intérêt de la méthode. Il s'agit d'une plaque rectangulaire, de dimensions $0.654 * 0.527 * 0.003m^3$. La taille de ce modèle est de 3000 ddls environ valable dans la bande fréquentielle $[0-300Hz]$. La structure est excitée au centre par une force harmonique de pulsation $\omega = 135rad/s$, d'amplitude $100N$. La dissipation est introduite sous forme d'amortissement visqueux proportionnel à la raideur linéaire $K_0$ avec un coefficient $\alpha = 0.01$ et $B = \alpha K$. La méthode de réduction proposée avec un choix de 10 vecteurs propres et un développement à l'ordre 5 est comparée
au résultat donné par le modèle complet. La taille du modèle réduit varie selon la SVD itérative mais reste aux alentours de 30 à 40 vecteurs ; ce qui représente une réduction de taille de 99%. Les résultats sont également comparés à un modèle prédictif utilisant une base modale actualisée à chaque pas de temps comme méthode de réduction. Les figures 1 et 2 montrent respectivement la variation du ratio entre le déplacement transversal et l'épaisseur de la plaque en fonction du temps ainsi que la vitesse en fonction du temps au point d'excitation.

<figure><img src="image_2.png" /><figcaption>FIG. 1 – *w/h* en fonction du temps : Modèle complet, Approximations Combinées (CA), Base Modale Actualisée (BMA)</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 2 – Vitesse en fonction du temps : Modèle complet, Approximations Combinées (CA), Base Modale Actualisée (BMA)</figcaption></figure>
En visualisant les deux figures, on constate que les méthodes de réduction prédisent correctement la réponse du modèle complet. L'avantage de la méthode des approximations combinées par rapport à la base modale actualisée est essentiellement au niveau du coût, puisque la résolution du problème aux valeurs propres à chaque itération est nécessaire.

Concernant le temps de calcul CPU de chaque modélisation, le modèle complet est le plus rapide en terme de convergence. En effet, ce modèle présente une meilleure prédiction dans le schéma d'intégration. Le nombre d'itérations de correction étant inférieur aux autres, le pas de temps adaptatif est plus grand. En effet, pour une période de calcul T donnée, le pas de temps est discrétisé en n itérations pour le modèle complet alors qu'avec la CA on a 2.5 * n itérations environ. La procédure la plus coûteuse étant l'évaluation de la force non-linéaire à chaque instant, elle occupe 99% du temps de calcul pour les différents modèles. Cette contrainte empêche la comparaison en terme de coût des différents modèles. Des travaux sont en cours pour corriger ce problème.

# 5 Conclusion

Une méthode de réduction de modèle a été présentée. Elle est basée sur la méthode des approximations combinées où il s'agit de développer les vecteurs propres du système linéaire sous forme d'une série binomiale en fonction de la non-linéarité.

Cette méthode conduit à des résultats comparables à ceux obtenus avec le modèle complet, ce qui valide sa mise en oeuvre. Le temps CPU de cette méthode reste à optimiser en comparaison avec le calcul complet. Toutefois, cette méthode avec son bon niveau de prédicтивité demeure incontournable pour mener des études d'optimisation ou de robustesse sur des modèles non linéaires couplés ou non de grande taille où le calcul complet peut devenir impraticable ou tout simplement incompatible avec des délais de conception raisonnables. Elle a été appliquée dans le domaine temporel avec une méthode d'intégration implicite. Les travaux à venir visent à adapter cette méthode au domaine fréquentiel en utilisant la balance harmonique (BH) comme outil de résolution. L'application à un problème vibro-acoustique non-linéaire est aussi envisagée; les phénomènes de couplage devront être pris en compte dans la base de réduction.

# Références

[1] GERADIN, M. 1993 Théorie des Vibrations.

[2] GUEDRI, M., WEISSER, T., BOUHADDI, N. 2010 Reanalysis of Nonlinear Structures by a Reduction Method of Combined Approximations *Proceedings of The Tenthe International Conference on Computational Structures Technology*.

[3] KERSCHEN, G., GOLINVAL, J.C., VAKAKIS, A.F., BERGMAN, L.A. 2005 The Method of Proper Orthogonal Decomposition for Dynamical Characterization and Order Reduction of Mechanical Systems : An Overview *Nonlinear Dynamics* **41** 147-169.

[4] KIRSCH, U. 2003 A Unified Reanalysis Approach for Structural Analysis, Design and Optimisation *Structural and Multidisciplinary Optimization* **25** 67-85.

[5] KIRSCH, U., BOGOMOLNI, M., SHEINMAN, I. 2006 Nonlinear Dynamic Reanlaysis of Structures by Combined Approximations *Comput. Methods Appl. Mech. Engrg* **195** 4420-4432.

[6] PERIGNON, F. 2004 Vibrations Forcées De Structures Minces, Elastiques, Non Linéaires *Thèses de l'Université Aix-Marseille II*.

[7] PESHECK, E., PIERRE, C. 2002 A New Galerkin-Based Approach for Accurate Non-Linear Normal Modes Through Invariant Manifolds *Journal of Sound and Vibration* **249**(5) 971-993.

[8] PONTHOT, J.P. 1995 Mécanique de Milieux Continus Solides en Grandes Transformations et Traitement Unifié par la Méthode des Eléments Finis *Thèses de l'Université de Liège*.

[9] ZIENKIEWICZ, O.C., TAYLOR, R.L. 2000 The Finite Element Method, Fifth Edition, Volume 2 : Solid Mechanics.