r5.06.01 Réduction de modèle en dynamique linéaire et non-linéaire : Méthode de RITZ

Résumé :

Ce document présente le principe de réduction de modèle par projection sur base réduite (méthode de Ritz).

La base le plus couramment utilisée est la base modale.

Les problèmes de troncature dus à l’utilisation d’une base réduite sont évoqués. Des corrections de troncature sont proposées.

La description et les propriétés des algorithmes de résolution du système d’équations différentielles du second ordre obtenue en analyse transitoire sont présentées dans le document [R5.06.04].

Méthodes de réduction de Ritz en linéaire

Description générale

Formulation continue

La méthode de Ritz consiste à projeter le déplacement sur une base restreintede fonctions vérifiant les conditions cinématiques du problème:

\(\tilde{u}(M,t)=\sum_{i=1}^{n}{\eta}_{i}(t){\psi}_{i}(M)\) éq 2.1.1-1

Le déplacement est décrit par une série de formes indépendantes \(\left\lbrace {\psi}_{i}(M);i=1\dots n\right\rbrace\) multipliées par des amplitudes fonctions du temps \(\left\lbrace {\eta}_{i}(t);i=1\dots n\right\rbrace\) .

La difficulté consiste à définir cette famille de forme \(\left\lbrace {\psi}_{i}(M);i=1\dots n\right\rbrace\) qui contrairement aux fonctions de forme de la méthode des éléments finis sont non nulles sur la plus grande partie de la structure.

La qualité de l’approximation est liée au fait que les déplacements obtenus ont une bonne approximation dans le sous-espace engendré par \(\text{Vect}\left\lbrace {\psi}_{i}(M),i=1\dots n\right\rbrace\) .

Projection sur base modale

On sait que les modes propres \(\left\lbrace {\varphi}_{i}(M);i=1\dots \infty \right\rbrace\) engendrent l’espace des champs cinématiquement admissibles. Le déplacement se décompose selon :

\(u(M,t)=\sum_{i=1}^{\infty}{\eta}_{i}(t){\varphi}_{i}(M)\) éq 2.1.1-2

L’option la plus couramment utilisée pour la méthode de Ritz consiste alors à prendre comme base de projection les \(n\) premiers modes:

\(\tilde{u}(M,t)=\sum_{i=1}^{n}{\tilde{\eta}}_{i}(t){\varphi}_{i}(M)\) éq 2.1.1-3

Le déplacement obtenu est une approximation du déplacement réel.

Il peut être intéressant d’ajouter au \(n\) premiers modes, d’autres formes (voir [§2.6.2]).

Approximation éléments finis

Dans le cas d’une approximation du déplacement par éléments finis le déplacement est déjà approché dans l’espace des fonctions de formes:

\({u}^{h}(M,t)=\sum_{i=1}^{\text{Nh}}{q}_{i}(t){N}_{i}(M)\) éq 2.1.2-1

On note \(U\) le vecteur des degrés de liberté du déplacement: \(U(t)=\left[{q}_{1}(t),{q}_{2}(t),\dots {q}_{\text{Nh}}(t)\right]\) ;

Méthode de Ritz en dimension finie

Si \(n<\text{Nh}\) , la méthode de Ritz appliquée au champ \(u(M,t)\) vient alors comme une seconde approximation:

\(\tilde{U}(t)=\sum_{i=1}^{n}{\eta}_{i}(t){\psi}_{i}\) éq 2.1.2-2

avec \(\left\lbrace {\psi}_{i,i=1\dots n}\right\rbrace\) la base des \(n\) vecteurs indépendants et cinématiquement admissibles. On pose

\(\Psi =\left[{\psi}_{1},{\psi}_{2},{\psi}_{3,\dots },{\psi}_{n}\right]\) . D’où l’écriture matricielle: \(U=\Psi \eta\) éq 2.1.2-3

Projection sur base réduite

Considérons le système différentiel suivant obtenu par une méthode d’éléments finis :

\(\mathrm{M}\ddot{U}+\mathrm{C}\dot{U}+\mathrm{K}U=\mathrm{F}\text{}U\in {R}^{\text{Nh}}\) éq 2.2-1

La solution recherchée sous la forme [éq 2.1.2-3]. En considérant la même forme pour le déplacement virtuel, il vient:

\({\Psi}^{T}M\psi \ddot{\eta}+{\Psi}^{T}C\Psi \dot{\eta}+{\Psi}^{T}K\Psi \eta ={\Psi}^{T}F\text{}\eta \in {R}^{n}\) éq 2.2-2

où : \(\eta\) est le vecteur des déplacements généralisés, \(\stackrel{ˉ}{K}={\Psi}^{T}K\Psi \text{et}\stackrel{ˉ}{M}={\Psi}^{T}M\Psi\) sont appelées respectivement matrices de raideur et de masse généralisées.

Le système [éq 2.2-2] est généralement un système différentiel couplé, les matrices généralisées qui le composent sont dans le cas général pleines même si au départ les matrices \(M\) et \(K\) étaient creuses. On perd donc la structure particulière au profit d’une taille de problème beaucoup plus réduite \(n\times n\) .

Dans le cas général, le système [éq 2.2-2] ne fournit qu’une solution approchée du système [éq 2.2-1]. L’erreur que l’on commet est appelée erreur de troncature.

On n’a aucune information sur la valeur de cette erreur. Elle peut être très grande si le sous-espace de projection est mal choisi. On sait seulement que cette erreur diminue quand la taille de la base de projection augmente.

Si l’on dispose d’information a priori sur la forme de la solution, on peut choisir de façon efficace la base de projection de façon à minimiser cette erreur. Par exemple, si l’on sait que la solution n’est constituée que de mouvements de corps solide, on peut restreindre à 6 la dimension de l’espace.

Par la suite, on choisit la base des modes propres comme base de projection.

Projection sur base modale

Modes propres

Les modes sont définis comme les couples \(\left\lbrace {({\omega}^{{h}_{i}},{\Phi}^{{h}_{i}})}_{i=1\dots \text{Nh}}\right\rbrace\) solutions de l’équation:

\((K-{\omega}^{2}M)\Phi =0\) éq 2.3-1

Remarque:

Il convient de vérifier queles modes calculés par approximation éléments finis sont suffisamment représentatifs : \(({\omega}^{{h}_{i}},{\Phi}^{{h}_{i}})\approx ({\omega}_{i},{\phi }_{i})\) *.* On peut considérer que l’approximation éléments finis est correcte lorsque les déformées modales présentent une longueur d‘onde supérieure à la taille des mailles du maillage (la notion de longueur d’onde est une généralisation de la notion définie sur l’équation des ondes, on peut la définir comme deux fois la longueur entre deux nœuds de la déformée modale).

Par la suite on omet, l’indice \(h\) correspondant à l’approximation éléments finis.

Quotient de Raleigh : interprétation énergétique

Les pulsations et formes propres peuvent être définies comme les solutions du problème de minimisation suivant :

\(\forall i\in \left[1,\text{Nh}\right]:{\Phi}_{i}\) minimise dans le sous espace \({R}^{\text{Nh}}-\text{Vect}\left\lbrace {\Phi}_{j},j\in \left[0,i-1\right]\right\rbrace\) la fonctionnelle:

\(R(X)=\frac{{X}^{t}\text{KX}}{{X}^{t}\text{MX}}\text{on}\text{pose}:\) \({\omega}_{i}^{2}=\frac{{\Phi}_{{i}^{t}}K{\Phi}_{i}}{{\Phi}_{{i}^{t}}M{\Phi}_{i}}=R({\Phi}_{i})\) éq 2.3-2

Méthode de réduction

Une méthode de réduction très largement employée pour les problèmes linéaires est la méthode de recombinaison modale. Elle consiste à choisir comme base de projection les \(n\) premiers modes propres de la structure \(\left\lbrace {\Phi}_{i,i=1\dots n}\right\rbrace\) .

\(\tilde{U}(t)=\sum_{i=1}^{n}{\eta}_{i}(t){\Phi}_{i}\) éq 2.3-3

Considérons toujours le système différentiel suivant:

\(M\ddot{U}+C\dot{U}+KU=F\text{}U\in {R}^{\text{Nh}}\) éq 2.3-4

Les modes propres \(\left\lbrace {\Phi}_{i,i=1\dots \text{Nh}}\right\rbrace\) ont la propriété d’être \(M\) et \(K\) orthogonaux, c’est-à-dire qu’on a les relations suivantes:

\(\begin{array}{}{\Phi}_{i}^{T}M{\Phi}_{j}={m}_{i}{\delta}_{ij}\\ {\Phi}_{i}^{T}K{\Phi}_{j}={k}_{i}{\delta}_{ij}\end{array}\)

\(\delta\) est le symbole de KRONECKER

\({m}_{i}\) est appelée masse modale ou masse généralisée du mode \(i\)

\({k}_{i}\) est appelée rigidité modale ou rigidité généralisée du mode \(i\)

Les matrices projetées de \(M\) et \(K\) sur la base des modes propres sont donc diagonales ; c’est un des avantages qui a motivé l’emploi de la base modale comme base de projection. Le système [éq2.3‑4] projeté sur la base des premiers modes propres du système s’écrit:

\((\begin{array}{ccc}\setminus & 0& 0\\ 0& {m}_{i}& 0\\ 0& 0& \setminus \end{array})\ddot{\eta}+{\Phi}^{T}C\Phi \dot{\eta}+(\begin{array}{ccc}\setminus & 0& 0\\ 0& {k}_{i}& 0\\ 0& 0& \setminus \end{array})\eta ={\Phi}^{T}{F}_{\text{ext}}\) éq 2.3-5

La projection de la matrice \(C\) n’a aucune raison en toute généralité d’être également diagonale. Si le système est fortement amorti (présence d’amortisseurs sur la structure), cette matrice ne sera pas diagonale.

Remarque :

Contrairement à ce que font beaucoup de logiciels, Code_Aster permet dans ce cas d’intégrer le système d’équations modales couplées sans diagonalisation de la matrice d’amortissement généralisé. La méthode d’intégration est dans ce cas une méthode implicite de NEWMARK ou explicite EULER.

Par contre, si le seul amortissement entrant en jeu est un amortissement structurel (dissipation interne du matériau pour une structure homogène) il est alors licite de faire l’hypothèse d’un amortissement proportionnel, encore appelé hypothèse de BASILE, dans ce cas \(C\) s’exprime comme combinaison linéaire de \(M\) et \(K\) (amortissement de RAYLEIGH), et sa projection sur les modes propres est diagonale(cf doc [R5.05.04] sur la modélisation de l’amortissement).

Dans ce cas, le système [éq 2.3-4] se scinde en p équations différentielles linéaires du second ordre découplées. La réponse du système est alors la recombinaison de la réponse de p oscillateurs simples associés aux modes propres, d’où l’expression de « superposition modale » utilisée couramment.

Chaque équation différentielle s’écrit \({m}_{i}\) :

\({m}_{i}{\ddot{\eta}}_{i}+{c}_{i}{\dot{\eta}}_{i}+{k}_{i}{\eta}_{i}={f}_{i}\) éq 2.3-6

ou encore en divisant par la masse modale:

\({\ddot{\eta}}_{i}+2{\xi}_{i}{\omega}_{i}{\dot{\eta}}_{i}+{\omega}_{i}^{2}{\eta}_{i}=\frac{{f}_{i}}{{m}_{i}}\) éq 2.3-7

avec:

\({\xi}_{i}\text{amortissement modal réduit}=\frac{{c}_{i}}{{c}_{\text{critique}}}=\frac{{c}_{i}}{2{m}_{i}.{\omega}_{i}}\)

Cette équation peut être résolue très simplement dans le domaine fréquentiel:

\({\stackrel{ˆ}{\eta}}_{i}=\frac{{\stackrel{ˆ}{\mathrm{f}}}_{i}(\omega )}{{m}_{i}.(-{\omega}^{2}+2.{\xi}_{i}{\omega}_{i}\omega +{\omega}_{i}^{2})}\) éq 2.3-8

où \(\stackrel{ˆ}{}\) représente la transformée de FOURIER et \(\omega\) la fréquence d’excitation.

Des méthodes numériques particulières telles l’intégrale de DUHAMEL permettent de passer cette expression du domaine fréquentiel au domaine temporel. (voir par exemple doc [R5.05.01] sur une méthode d’intégration temporelle).

Erreur de troncature modale

Dans le cas de la recombinaison modale avec amortissement proportionnel, on peut mettre en évidence l’erreur de troncature que l’on commet en projetant sur la base des premiers modes propres du système. En effet, si l’on considère la base complète des n modes propres du problème discrétisé, il y a équivalence entre le problème initial et le problème projeté. Donc la solution exacte du problème discrétisé par éléments finis s’écrit:

\(U=\sum_{i}^{\text{Nh}}{\eta}_{i}{\Phi}_{i}\)

où les coordonnées généralisées sont solution de:

\({\ddot{\eta}}_{i}+2{\xi}_{i}{\omega}_{i}{\dot{\eta}}_{i}+{\omega}_{i}^{2}{\eta}_{i}=\frac{{f}_{i}}{{m}_{i}}\)

la sommation s’étendant sur tous les modes propres du système (de taille finie).

En résolvant le problème avec un nombre réduit de modes propres, \(n<\text{Nh}\) . La solution obtenue est la suivante:

\(\tilde{U}=\sum_{i=1}^{n}{\eta}_{i}{\Phi}_{i}\)

L’erreur commise en tronquant la base de représentation de la solution est donc:

\(E=U-\tilde{U}=\sum_{i=n+1}^{\text{Nh}}{\eta}_{i}{\Phi}_{i}\) éq 2.4-1

Dans le domaine fréquentiel l’expression de l’erreur est:

\({\stackrel{ˆ}{E}}_{(\omega )}=\stackrel{ˆ}{U}-\stackrel{ˆ}{\tilde{U}}=\sum_{i=n+1}^{\text{Nh}}\frac{{\Phi}_{i}^{t}{\stackrel{ˆ}{F}}_{(\omega )}}{{m}_{i}}.\frac{1}{{\omega}_{i}^{2}-{\omega}^{2}+\mathrm{2j}{\xi}_{i}\omega {\omega}_{i}}.{\Phi}_{i}\) éq 2.4-2

la sommation s’effectue sur tous les modes négligés du système.

Étudions la réponse relative \(\eta /{\eta}_{\text{statique}}\) d’un oscillateur à une excitation purement sinusoïdale de fréquence variable (schéma ci-dessous), avec \({\eta}_{\text{statique}}\) les coefficients de la réponse statique correspondant à une force statique. On peut distinguer trois intervalles dans le spectre où l’oscillateur a un comportement différent. En basse fréquence \((\omega \text{<<}{\omega}_{0})\) l’oscillateur a une réponse statique. Autour de \({\omega}_{0}\) l’oscillateur a une réponse dynamique (amplification du mode), et à haute fréquence l’oscillateur répond de façon inertielle (\(\frac{1}{{\omega}^{2}}\) terme prépondérant).

Supposons que l’excitation du système, définie par le vecteur \(F(\omega )\) , soit à bande étroite, notamment qu’elle soit nulle pour des fréquences supérieures à \({\omega}_{\max}\) donné.

Dans ce cas, pour représenter correctement la réponse du système linéaire, il faut assurément prendre en compte tous les modes ayant une pulsation inférieure à \({\omega}_{\max}\) , car ces derniers vont répondre de façon dynamique à l’excitation.

Par contre, les modes tels que \({\omega}_{i}\text{>>}{\omega}_{\max}\) ont quand même une contribution statique à la réponse du système. Ce sont souvent ces modes que l’on ne prend pas en compte.

En faisant un développement limité en \(\omega\) au voisinage de \(0\) . On obtient la partie principale de l’erreur qui est:

\({\stackrel{ˆ}{E}}_{(\omega )}=\stackrel{ˆ}{U}-\stackrel{ˆ}{\tilde{U}}=\sum_{i=p+1}^{n}\frac{{\Phi}_{i}^{T}.{\stackrel{ˆ}{F}}_{(\omega )}}{{k}_{i}}.(1-2j{\xi}_{i}\frac{\omega}{{\omega}_{i}}+0(\frac{\omega}{{\omega}_{i}})).{\Phi}_{i}\) éq 2.4-3

L’erreur est d’autant plus petite que les rigidités généralisées des modes négligés sont grandes. En principe donc, il faudra prendre tous les modes les plus souples jusqu’à ce que la souplesse résiduelle d’un mode supplémentaire soit en valeur relative négligeable par rapport à la somme des souplesses déjà prises en compte:

\(\frac{1}{{k}_{n+1}}\text{<<}\sum_{i=1}^{n}\frac{1}{{k}_{i}}\)

Cependant, on observe qu’en négligeant les modes de fréquence élevée on commet une erreur systématique sur la réponse du système (même en basse fréquence). Il existe différentes possibilités que nous allons détailler maintenant pour corriger la réponse dans la plage \(\left[0,{\omega}_{\max}\right]\) où l’on a choisi les modes.

Corrections de la troncature modale

Pour pallier au problème de troncature dû aux modes négligés, il faut essayer d’estimer leur effet dans le domaine de fréquence \(\left[0,{\omega}_{\max}\right]\) qui nous intéresse. Nous avons vu que les modes négligés ayant une pulsation propre telle que \({\omega}_{i}\text{>>}{\omega}_{\max}\) ont une contribution dite statique à la réponse du système dans le domaine \(\left[0,{\omega}_{\max}\right]\) . Les techniques de correction consistent à calculer cette contribution statique.

Correction statique a posteriori

L’erreur de troncature, en ne considérant que la réponse statique des modes négligés (transformée inverse de la partie principale de l’erreur) est:

\({E}_{(t)}=U-\tilde{U}\approx \sum_{i=n+1}^{\text{Nh}}\frac{{\Phi}_{i}^{T}.{F}_{(t)}}{{k}_{i}}.{\Phi}_{i}\) éq 2.5.1-1

Mais a priori les modes négligés ainsi que leurs rigidités généralisées sont inconnus. Par contre, on sait déterminer la réponse statique complète du système à un chargement \(F\) , cette dernière vaut:

\(U={K}^{-1}.{F}_{(t)}\approx \sum_{i=1}^{\text{Nh}}\frac{{\Phi}_{i}^{T}.{F}_{(t)}}{{k}_{i}}.{\Phi}_{i}\)

La correction qu’il faut apporter est donc:

\(\sum_{i=n+1}^{\text{Nh}}\frac{{\Phi}_{i}^{T}.{F}_{(t)}}{{k}_{i}}.{\Phi}_{i}\approx {K}^{-1}.{F}_{(t)}-\sum_{i=1}^{n}\frac{{\Phi}_{i}^{T}.{F}_{(t)}}{{k}_{i}}.{\Phi}_{i}\)

La solution corrigée de la réponse du système vaut donc:

\(U=\tilde{U}+E\approx \tilde{U}+{K}^{-1}.{F}_{(t)}-\sum_{i=1}^{n}\frac{{\Phi}_{i}^{T}.{F}_{(t)}}{{k}_{i}}.{\Phi}_{i}\) éq 2.5.1-2

Cette correction est appelée a posteriori, car elle n’intervient pas dans la résolution dynamique du système linéaire et peut n’être calculée que bien ultérieurement. Si \({F}_{(t)}\) se décompose en \(k\) produits de fonctions du temps par des fonctions des coordonnées d’espace, cette correction nécessite une factorisation de \(K\) et \(k\) résolutions.

Cette méthode a l’avantage de ne pas augmenter le nombre de vecteurs pris en compte dans la base. Cette méthode est applicable dans le cas d’une excitation à bande étroite, ou tout au moins ayant une fréquence de coupure connue. La correction est exacte dans le domaine basse fréquence mais peut fausser la réponse du système en haute fréquence [§Annexe1].

Adjonction de modes statiques à la base

Supposons que le chargement \(F(t)\) s’écrive:

\(\mathrm{F}(t)=\sum_{j}{\alpha}_{j}(t).{\mathrm{F}}_{j}\)

La seconde façon de corriger l’erreur de troncature consiste à ajouter à la base des modes propres initiaux des modes statiques \({\Psi}_{j}\) définis comme la déformée à chaque effort \({F}_{j}\) donné:

\({\Psi}_{j}={K}^{-1}.{F}_{j}\) éq 2.5.2-1

La nouvelle base de projection à considérer est la suivante:

\(\stackrel{ˆ}{\Phi}=[{\Phi}_{1},{\Phi}_{2},...,{\Phi}_{p},{\Psi}_{1},{\Psi}_{2},...,{\Psi}_{m}]=[\Phi \Psi ]\) éq 2.5.2-2

Les composantes généralisées à utiliser sont les suivantes:

\(\stackrel{ˆ}{\eta}=[{\eta}_{1},{\eta}_{2},...,{\eta}_{p},{\mu}_{1},{\mu}_{2},...,{\mu}_{m}]=[\eta ,\mu ]\) éq 2.5.2-3

Le problème projeté sur la base complétée est:

\(\left[\begin{array}{cc}(\begin{array}{ccc}\setminus & 0& 0\\ 0& {m}_{i}& 0\\ 0& 0& \setminus \end{array})& {\Phi}^{T}.M.\Psi \\ {\Psi}^{T}.M.\Phi & {\Psi}^{T}.M.\Psi \end{array}\right]\left\lbrace \begin{array}{c}\ddot{\eta}\\ \ddot{\mu}\end{array}\right\rbrace +\left[\begin{array}{cc}(\begin{array}{ccc}\setminus & 0& 0\\ 0& {k}_{i}& 0\\ 0& 0& \setminus \end{array})& {\Phi}^{T}.K.\Psi \\ {\Psi}^{T}.K.\Phi & {\Psi}^{T}.K.\Psi \end{array}\right]\left\lbrace \begin{array}{c}\eta \\ \mu \end{array}\right\rbrace =\left\lbrace \begin{array}{c}{\Phi}^{T}.F.\alpha \\ {\Psi}^{T}.F.\alpha \end{array}\right\rbrace\) éq 2.5.2-4

On constate que l’on a perdu le caractère diagonal des matrices généralisées, mais l’avantage obtenu est que la base complétée avec des modes statiques permet de représenter correctement le comportement à basse fréquence du système initial.

Par exemple il est simple de montrer qu’à fréquence nulle la solution de ce système est:

\(\eta =0\) et \(\mu =\alpha\) qui est la solution exacte du problème statique initial.

On présente en annexe 1, la comparaison sur un système discret à 3 degrés de liberté entre la solution exacte, la solution projetée sur 1 mode, celle projetée sur un mode avec une correction statique et le solution constituée par 1 mode propre et 1 mode statique.

On s’aperçoit que l’adjonction de modes statiques permet d’étendre au delà de l’intervalle \(\left[0,{\omega}_{\max}=\max({\omega}_{j})\right]\) la bonne représentation dynamique du système. Cette technique semble donc très intéressante, elle a le mérite de procéder immédiatement à la correction ce qui sera intéressant pour les méthodes non linéaires où l’on a besoin de la connaissance des déplacements physiques à chaque pas de temps.

La mise en œuvre pratique de cette technique est précisée dans [U2.06.04].

Extension des méthodes de réduction de Ritz en non‑linéaire

Problème général

Le problème de dynamique non-linéaire discrétisé sans amortissement peut le plus souvent se mettre sous la forme suivante:

\(M.\ddot{X}+G(X)=F(t)\text{}X\in {R}^{n}\) éq 3.1-1

\(G(X)\) est une fonction non-linéaire de \(X\) qui représente les forces internes du système ainsi que toutes les autres forces qui sont dépendantes du déplacement, \(F\) le vecteur des forces externes et \(M\) la matrice de masse du système.

La matrice de rigidité tangente du système est par définition:

\({K}_{(x)}^{\text{tg}}=\frac{\partial G}{\partial X}(x)\) éq 3.1-2

Elle permet de définir à chaque instant une base modalepar :

\((-{\omega}_{i(X)}^{{\text{tg}}^{2}}M+{K}^{\text{tg}}).{\Phi}_{i(X)}^{\text{tg}}=0\) éq 3.1-3

Les modes ainsi définis dépendent de \(X\) , donc de l’instant \(t\) .

Sachant que le calcul de la base modale est très coûteux en temps de calcul, l’idée de vouloir projeter à chaque pas de temps le modèle sur une base modale, puis de résoudre, est sans intérêt par rapport à une résolution directe.

La méthode la plus couramment utilisée consiste à définir une base de projection en ajoutant aux modes calculés sur une configuration initiale des formes permettant de projeter la non linéarité.

Exemple: dans le cas où la non linéarité provient d’un choc ponctuel, on propose d’enrichir la base modale avec des modes statiques permettant de projeter l’effort subi par la structure durant le choc [R5.06.04].

La méthode de Ritz reste toujours pertinente dans les calculs non linéaires, si la base choisie permet de projeter correctement les déplacements et les efforts.

Le problème non linéaire projeté sur une base \(\Psi\) quelconque s’écrit:

\({\Psi}^{T}M\Psi .\ddot{\eta}+{\Psi}^{T}G(X)={\Psi}^{T}F\text{}\eta \in {R}^{n}\) éq 3.1-4

Deux possibilités sont alors possibles:

• les non linéarités sont localisées et l’on peut évaluer la non linéarité sur la base de projection : \(G(X)=G(\Psi \eta )\) .

• Le problème à résoudre est un système différentiel non linéaire en \(\eta\) de taille plus petite. Différentes stratégies sont possibles pour résoudre ce problème, dépendant essentiellement de la technique d’intégration que l’on souhaite utiliser.

• les non linéarités sont globales, et il faut repasser dans l’espace des degrés de liberté physiques pour calculer les forces internes : \(G(X)\) .

Cette seconde méthode est plus coûteuse elle est beaucoup moins courante.

On présente en annexe 2, la comparaison sur un système à 3 degrés de liberté avec une non linéarité en \({x}^{3}\) entre la solution exacte et la solution obtenue par la méthode ci-dessus avec 1 puis 2 modes. On s’aperçoit qu’il est nécessaire de prendre plus de modes en compte que pour le problème linéaire. Par contre, sur cet exemple 2 modes suffisent très bien à décrire le système.

Indication de l’erreur de projection

Pour les problèmes non linéaires le sens physique du nombre de modes à prendre en compte est complètement perdu, et si les méthodes de réduction donnent toujours une solution il faut savoir le degré de confiance que l’on peut leur accorder. Une façon de procéder, qui est un peu coûteuse mais indispensable est de calculer le résidu du système initial à chaque pas de temps. Il se définit par:

\(R=M.\ddot{X}+G(X)-F(t)\)

Ce vecteur résidu n’est malheureusement pas nul, c’est uniquement sa projection sur la base utilisée qui l’est.

Une norme peut alors être calculée pour ce résidu ; plus la norme du résidu sera petite plus on pourra accorder de confiance à la solution.

Pour utiliser une valeur relative, on a intérêt à calculer la fraction suivante:

\(r=\frac{\parallel R\parallel }{\max(\parallel F(t)\parallel ,\parallel G(X)\parallel ,\parallel M.\ddot{X}\parallel )}\) éq 3.2-1

Remarque :

Cet indicateur n’est pas implanté actuellement dans Code_Aster.

Utilisation dans Code_Aster

Dans Code_Aster, les méthodes de Ritz sont utilisables en transitoire essentiellement par l’opérateur DYNA_TRAN_MODAL [U4.53.21].

Une phase de projection des matrices de rigidité et de masse sur une base de vecteurs est réalisée par les opérateurs PROJ_MATR_BASE [U4.63.12] et PROJ_VECT_BASE [U4.63.13].

Le problème dynamique généralisé est ensuite résolu dans l’opérateur DYNA_TRAN_MODAL par un schéma d’intégration explicite (EULER ou DEVOGELEARE) ou implicite (NEWMARK). Les caractéristiques et propriétés des schémas d’intégration sont présentés dans la note [R5.06.04]. Pour les structures pour lesquelles l’hypothèse de BASILE ne s’applique pas (amortissement non proportionnel) on projettera aussi la matrice d’amortissement qui ne devient pas diagonale. L’intégration du système couplé se fait alors obligatoirement avec le schéma implicite (NEWMARK) ou explicite (EULER).

Les non-linéarités localisées sont spécifiées directement dans l’opérateur DYNA_TRAN_MODAL. On peut introduire des non-linéarités localisées du type choc et frottement (voir [R5.06.03] Modélisation des chocs et frottements), des forces modales fonction du déplacement ou de la vitesse.

Les corrections statiques de troncature a posteriori sont disponibles dans le cas d’une excitation unique (voir [R4.05.01] Réponse sismique).

L’adjonction de modes statiques peut se faire en utilisant au préalable les opérateurs MODE_STATIQUE [U4.52.14] et DEFI_BASE_MODALE [U4.64.02]. Quand le problème comporte des non-linéarités seuls les schémas explicites peuvent être utilisés.

Pour les non-linéarités globales [éq 3.1-4], il est possible d’utiliser la commande DYNA_NON_LINE [U4.53.01] avec l’opérande PROJ_MODAL, qui permet de calculer à chaque pas de temps les forces internes en fonction des degrés de liberté physiques, puis de projeter le problème sur une base modale.

Une opération de retour à la base physique est ensuite nécessaire pour obtenir les grandeurs physiques tels que déplacement, vitesse ou accélération sur la structure. Cette opération est réalisée par l’opérateur REST_GENE_PHYS [U4.63.31] pour l’ensemble de la structure ou RECU_FONCTION, (mot clef facteur RESU_GENE) [U4.32.01] pour le suivi d’une grandeur en un nœud. Dans le cas des non-linéarités globales (calcul avec DYNA_NON_LINE) on utilise l’opérateur REST_COND_TRAN [U4.63.33].

Plus généralement, l’approche de Ritz peut être utilisée en calcul harmonique par la commande DYNA_LINE_HARM [U4.53.22] et en densité spectrale de puissance par la commande DYNA_ALEA_MODAL [U4.53.23].

Enfin la sous-structuration dynamique peut être considérée comme une méthode de Ritz spécifique [R4.06.02].

Bibliographie

1. BATHE - WILSON: Finite Element procedures in Engineering Analysis

2. 13. GERADIN, D. RIXEN: Théorie des vibrations: application à la dynamique des structures – Masson 1993

3. J.F.IMBERT: Analyses des structures par éléments finis - cepadues éditions 1979

4. BELYTSCHKO - LIU - PARK: Innovative methods for non linear problems

5. EWINS D. J. Modal Testing : Theorie and practice Reserch Studies Press LTD

6. R.E. RICKELL: Non-linear dynamics by mode superposition Computer Methods in Applied Mechanics and Engineering (1976) vol 7

7. 16. LUKKUNAPRASIT: Dynamic response of an elastic, viscoplastic system in modal coordinates. Earthquake Engineering & Structural Dynamics (1980)

8. 7. JACQUART: Méthodes de RITZ en dynamique non-linéaire - Application à des systèmes avec choc et frottement localisé - Rapport EDF-DER HP-61/91.105

Description des versions du document

Version Aster	Auteur(s) Organisme(s)	Description des modifications
6	G. JACQUART (EDF R&D/AMV)	Texte initial
7.4	L.RATIER, G. JACQUART (EDF R&D/AMA, EDF/CNPE de Tricastin)

Considérons le système discret à trois masses suivant:

Les matrices de rigidité et de masse sont:

\(M=(\begin{array}{ccc}m& 0& 0\\ 0& m& 0\\ 0& 0& m\end{array})\text{}K=(\begin{array}{ccc}k& -k& 0\\ -k& \mathrm{2k}& -k\\ 0& -k& \mathrm{2k}\end{array})\)

Soit: \({\omega}_{0}^{2}=\frac{k}{m}\)

Les modes propres et leur pulsation valent:

\(\begin{array}{c}\begin{array}{}{\omega}_{1}^{2}=0,198{\omega}_{0}^{2},{m}_{1}=1,841\text{}m,{\Phi}_{1}=(\begin{array}{c}1\\ 0,802\\ 0,445\end{array})\\ \end{array}\\ \begin{array}{}{\omega}_{2}^{2}=1,555{\omega}_{0}^{2},{m}_{2}=2,863\text{}m,{\Phi}_{2}=(\begin{array}{c}1\\ -0,555\\ -1,247\end{array})\\ \end{array}\\ {\omega}_{3}^{2}=3,247{\omega}_{0}^{2},{m}_{3}=9,296\text{}m,{\Phi}_{3}=(\begin{array}{c}1\\ -2,247\\ 1,802\end{array})\end{array}\)

Comparons les réponses du système modélisé par un seul mode propre avec ou sans correction statique:

On constate que la correction statique permet de corriger la réponse à basse fréquence, le modèle avec 1 mode plus correction colle parfaitement à la solution exacte en basse fréquence. Par contre, en haute fréquence (au delà du premier mode), cette correction conduit à surestimer énormément la réponse. L’utilisation de la correction statique devra donc être utilisée avec prudence et dans le cadre d’une excitation à bande étroite.

Regardons ce que donne la méthode d’adjonction d’un mode statique.

Si l’on applique une force unitaire au point 1, la déformée statique vaut:

\({\Psi}_{s}=\frac{1}{k}(\begin{array}{c}3\\ 2\\ 1\end{array})\)

Les matrices projetées de masse et rigidité que l’on obtient sont les suivantes:

\(\stackrel{ˆ}{M}=(\begin{array}{c}1,841m\frac{5,049}{{\omega}_{0}^{2}}\\ \frac{5,049}{{\omega}_{0}^{2}}14\text{}\frac{m}{{k}^{2}}\end{array}\text{})\text{}\text{et}\text{}\stackrel{ˆ}{K}=(\begin{array}{c}0,365k\text{}1\\ \text{}1\text{}\frac{3}{k}\text{}\end{array})\)

ce système a pour fréquences propres:

\({\omega}_{1}^{2}=0,198{\omega}_{0}^{2}\) et \({\omega}_{1}^{2}=1,667{\omega}_{0}^{2}\)

La réponse du système modélisé avec un mode propre et un mode statique est la suivante:

On s’aperçoit que l’on corrige très bien en basse fréquence, (effet de correction statique), on modélise la dynamique du système bien au delà du premier mode pris en compte. Par contre l’effet du second mode est mal représenté (shift sur la fréquence).

Considérons le système discret à trois masses suivant:

Les matrices de raideur et masse sont:

\(M=(\begin{array}{ccc}m& 0& 0\\ 0& m& 0\\ 0& 0& m\end{array})\text{}K=(\begin{array}{ccc}k& -k& 0\\ -k& \mathrm{2k}& -k\\ 0& -k& \mathrm{2k}\end{array})\)

Rendons ce système non-linéaire en ajoutant un terme de force interne entre \(\mathrm{x1}\) et \(\mathrm{x2}\) cubique:

\(F=k.{(\mathrm{x1}-\mathrm{x2})}^{3}\)

Cherchons à évaluer la réponse de ce système à une excitation forcée de fréquence proche de la première fréquence propre du système linéaire (on a choisi \(\omega =0,18{\omega}_{0}\) ), avec une amplitude importante \({F}_{m}=3.k\) .

Dans cette configuration, la réponse du système ne peut être évaluée par la réponse du système linéaire (le terme cubique est bien trop important), il faut mettre en œuvre un calcul non-linéaire avec pseudo-forces comme on l’a montré en [§3.2].

On peut voir dans le rapport [bib8] les courbes des résultats transitoires de cette méthode, en prenant en compte un ou 2 modes propres du système linéaire initial.

Avec un seul mode propre, on s’aperçoit que l’on commet une erreur relativement importante (atteignant parfois 50%), par contre il est satisfaisant de constater que les extrema des vibrations sont assez bien prévus. On aurait pu espérer qu’en excitant en deçà de la première fréquence propre il aurait suffit d’un seul mode propre pour modéliser la réponse du système, on voit ici que ce n’est pas le cas. Comme on le constate souvent pour des systèmes non-linéaires, le système répond également avec des surharmoniques de la fréquence d’excitation.

Par contre, en prenant 2 modes propres pour modéliser la réponse de cette structure à 3 ddl, on obtient un résultat très satisfaisant (quelques% d’erreur sur l’amplitude), à l’œil on a du mal à distinguer la différence. Ceci montre qu’en choisissant une base de projection suffisamment riche on peut grâce à une méthode de pseudo‑forces très bien modéliser un système dynamique complexe avec des non-linéarités.