# Des solveurs implicites d'ordre supérieur pour les problèmes de dynamique non-linéaire des structures

Sanaa Boutmir, Bouazza Braikat, Mohammad Jamal, Noureddine Damil,
Bruno Cochelin, Michel Potier-Ferry

► To cite this version:

Sanaa Boutmir, Bouazza Braikat, Mohammad Jamal, Noureddine Damil, Bruno Cochelin, et al. Des solveurs implicites d'ordre supérieur pour les problèmes de dynamique non-linéaire des structures. Revue Européenne des Éléments Finis, 2004, 13, pp.449-460. 10.3166/reef.13.449-460. hal-00089046

HAL Id: hal-00089046

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

Submitted on 9 Aug 2006

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
# Des solveurs implicites d'ordre supérieur pour les problèmes de dynamique non linéaire des structures

Sanaa Boutmir* — Bouazza Braikat* — Mohammad Jamal*  
Noureddine Damil* — Bruno Cochelin** — Michel Potier-Ferry***

* Laboratoire de Calcul Scientifique en Mécanique
Faculté des Sciences Ben M'Sik, Université Hassan II - Mohammed Media
Avenue Cdt Driss El Harti, BP 9755, Sidi Othman, Casablanca, Maroc
b.braikat@univh2m.ac.ma, m.jamal@univh2m.ac.ma, n.damil@univh2m.ac.ma

**Laboratoire de Mécanique et d'Acoustique, UPR CNRS N° 7051, Ecole Supérieure de Mécanique de Marseille, Technopole de Chateau-Gombert, 13451, Marseille, France
cochelin@imtumn.imt-mrs.fr

*** Laboratoire de Physique et Mécanique des Matériaux UMR CNRS 7554, ISGMP, Université de Metz, Ile du Saulcy, 57045 Metz, France
potier-ferry@lpmm.univ-metz.fr

RÉSUMÉ. On développe dans ce travail quelques solveurs temporels implicites d'ordre élevé pour la résolution des problèmes dynamiques non linéaires des structures élastiques en transformations finies. Ces solveurs se basent sur la méthode de perturbation, la transformation homotopique et sur des techniques de discrétisation temporelle et spatiale. On montre que l'efficacité de ces algorithmes augmente par l'introduction des approximants de Padé. La performance et la comparaison de ces algorithmes avec d'autres solveurs classiques sont testées sur la vibration non linéaire d'une poutre élastique 2D et d'une plaque élastique.

ABSTRACT. We develop in this work, some implicit temporal high order solvers for solving non linear elastic structural dynamic problems involving finite deformations. These solvers are based on the perturbation method, the homotopical transformation and time space discretization techniques. Their accuracy is improved by the introduction of Padé approximants. Numerical calculations, compared with others classical solvers, are illustrated on forced nonlinear vibration problems of a 2D elastic beam and an elastic plate.

MOTS-CLÉS: solveurs implicites d'ordre élevé, dynamique non linéaire des structures, homotopie, perturbation.

KEYWORDS: implicit high order solvers, non linear dynamics of structures, homotopy, perturbation.
# 1. Introduction

Plusieurs activités de recherche ont été consacrées, durant ces deux dernières décades, à la mise au point d’algorithmes efficaces permettant d’accéder numériquement à la réponse d’une structure élastique non linéaire soumise à des sollicitations dynamiques. Il existe dans la littérature, portant sur la thématique, une panoplie d’approches numériques destinées à résoudre cette classe de problèmes non linéaires installationnaires [ARG 91] et [BAT 96]. Certaines d’entre elles, sont actuellement intégrées dans plusieurs codes de calcul. Les plus populaires sont soit de type itérative soit de type prédiction-correction. Ces méthodes sont, généralement, basées sur des intégrations directes des équations qui gouvernent la dynamique non linéaire du problème étudié et font appel aux procédures de discrétisation spatiale classiques MEF, MDF, MVF. L’intégration en temps est effectuée par un schéma en temps du type explicite-implicite. Ces méthodes fonctionnent très bien et permettent d’obtenir la réponse dynamique non linéaire mais elles peuvent exiger des temps de calcul CPU très importants surtout pour les algorithmes implicites.

On propose dans ce travail, deux types d’algorithmes implicites d’ordre supérieur pour la résolution des problèmes de dynamique non linéaires des structures élastiques. L’intégration en temps est faite par le schéma classique implicite de Newmark. Ces algorithmes sont des variantes de la Méthode Asymptotique Numérique (MAN). Ils sont développés en couplant les transformations homotopiques, les procédures de discrétisation temporelle et la discrétisation spatiale et l’approche des développements asymptotiques. L’idée d’associer la technique d’homotopie avec la méthode asymptotique numérique pour la résolution des problèmes installationnaires non linéaires est apparue pour la première fois dans le travail de Cochelin et Compain [COC 00] axé sur l’étude de la vibration non linéaire d’une poutre 2D. Ce couplage a permis d’utiliser la même triangularisation de matrice pour de nombreux pas de temps. Une autre variante a été testée et validée sur des modèles de propagation des ondes non linéaires modélisés par les équations de Korteweg-De-Vries (KdV), de Fisher et sur le problème des oscillations non linéaires d’un pendule élastique bidimensionnel [JAM 02], [BRA 02] et [BRA 04].

Le point important de ces algorithmes est la possibilité de choisir la matrice d’itération et de calculer ainsi plusieurs pas de temps à l’aide d’une seule triangularisation. Les algorithmes proposés diffèrent principalement par le choix de l’ordre d’applications des étapes de base : discrétisation en temps, technique d’homotopie et technique de perturbation.

Nous montrerons que l’introduction des approximants de Padé permet d’augmenter la robustesse des algorithmes proposés. La performance de ces algorithmes est testée sur des exemples de vibrations non linéaires d’une poutre élastique 2D et d’une plaque élastique. Une comparaison des résultats obtenus par ces solveurs avec ceux des méthodes itératives et de Newton Raphson est présentée. Une discussion de l’effet du pas de temps sur la convergence des méthodes utilisées est aussi donnée.
## 2. Formulation du problème en élasticité non linéaire

### 2.1. Les équations dynamiques du système

Considérons une structure élastique 3D, occupant un volume de référence $\Omega$ de frontière $\partial\Omega$. La structure, faite d'un matériau homogène et isotrope de densité volumique $\rho$, est soumise à un déplacement $u_d$ sur une partie $\partial\Omega_u$, à des forces surfaciques $F$ sur la partie complémentaire $\partial\Omega_F$ et reçoit eventuellement des forces de volume $f$. Les quantités $u_d$, $F$ et $f$ sont des données définies respectivement sur $\Omega$ et $\partial\Omega$. La formulation variationnelle traduisant l'équilibre gouvernant le mouvement dynamique de cette structure et la loi de comportement sont :

$$ egin{cases} 
\int_{\Omega} \rho \ddot{u} d\Omega + \int_{\Omega} S : \delta\gamma d\Omega = \int_{\Omega} f \delta\gamma d\Omega + \int_{\partial\Omega_F} F \delta u ds \\ 
S = D : \gamma(u) = D : (\gamma^L(u) + \gamma^{NL}(u, u)) \quad [1] \\ 
u(t=0) = u_0, \quad \dot{u}(t=0) = \dot{u}_0 
\end{cases} 
$$

où $u$ et $\ddot{u}$ représentent respectivement le déplacement et l'accélération, $S$ est le second tenseur des contraintes de Piola-Kirchhoff, $D$ est le tenseur d'élasticité, $\gamma(u)$ est le tenseur non linéaire de déformation de Green-Lagrange, $\gamma^L(u)$ et $\gamma^{NL}(u, u)$ sont respectivement les parties linéaire et non linéaire de $\gamma(u)$ avec $\gamma^L(u) = \frac{1}{2}(\nabla u + {}^t\nabla u)$, $\gamma^{NL}(u, u) = \frac{1}{2} {}^t\nabla u \nabla u$, $u_0$ et $\dot{u}_0$ sont les conditions initiales à l'instant $t=0$.

### 2.2. Discrétisation par éléments finis

La discrétisation spatiale par la méthode des éléments finis des équations (1) conduit à [CRI 91], [BAT 96] :

$$ egin{cases} 
\sum_e \int_{\Omega^e} \rho^t [N] [N] d\Omega^e \{\ddot{q}\} + \sum_e \int_{\Omega^e} {}^t[B(q)] \{S\} d\Omega^e = \{F\} \\ 
\{S\} = [D]\{\gamma\} = [D]([B_L] + \frac{1}{2}[B_{NL}(q)])\{q\} \quad [2] \\ 
\{q(t=0)\} = \{q_0\}, \quad \{\dot{q}(t=0)\} = \{\dot{q}_0\} 
\end{cases} 
$$

où $\{q\}$ et $\{\ddot{q}\}$ sont respectivement les vecteurs de déplacement et d'accélération noyaux, $[N]$ est la matrice des fonctions de forme telle que $\{u\} = [N]\{q\}$, $[B(q)] = [B_L] + [B_{NL}(q)]$ avec $[B_L]$ et $[B_{NL}(q)]$ correspondant respectivement à la discrétisation des parties linéaire et non linéaire de $\gamma$ et $\Omega^e$ est le volume élémentaire. Le vecteur résidu $R$ est donné par l'expression :

$$ \{R\} = \sum_e \int_{\Omega^e} \rho^t [N] [N] d\Omega^e \{\ddot{q}\} + \sum_e \int_{\Omega^e} {}^t[B(q)] \{S\} d\Omega^e - \{F\} \quad [3] $$

Le problème (2) est généralement résolu, après discrétisation en temps, par des algorithmes de prédiction correction de type Newton [BAT 96], [ARG 91], [CRI 97] ou
par des méthodes itératives [BAT 96], [ZIE 87], [GER 93]. L'intégration temporelle est faite par le choix d'un schéma en temps de la famille de Newmark [NEW 59]. Ces méthodes permettent d'accéder à la solution du problème (2) selon une description point par point. Elles sont opérationnelles, quelle que soit la nature des non linéarités impliquées, mais elles peuvent exiger d'importants temps de calcul CPU.

Nous proposons par la suite deux familles d'algorithmes implicites d'ordre élevé pour résoudre le problème dynamique discrétisé (2).

### 3. Algorithmes implicites d'ordre élevé proposés

On propose deux algorithmes qui se différencient par le choix de l'homotopie. Dans le premier, l'homotopie est choisie à partir d'une linéarisation par rapport à l'instant initial et dans le second, elle est choisie à partir d'une linéarisation par rapport à l'instant t.

#### 3.1. Linéarisation par rapport à l'instant initial

Le principe de cet algorithme, pour résoudre le problème (2), repose sur les étapes suivantes : transformation homotopique, discrétisation temporelle et méthode perturbative.

i- Transformation homotopique : on introduit le changement de variable suivant :

$$ egin{cases} \{q(t)\} = \{q_0\} + \{v(t)\} \\ \{S(t)\} = \{S_0\} + \{s(t)\} \end{cases} \quad [4] $$

où $\{v(t)\}$ et $\{s(t)\}$ sont les nouvelles inconnues et $\{S_0\}$ représente la contrainte à l'instant $t = 0$. On injecte (4) dans (2) et on remplace, dans la partie linéaire de l'équation obtenue, la nouvelle expression de la loi de comportement, on obtient :

$$ egin{cases} [M]\{\ddot{v}\} + [K_t]\{v\} + \{FQ_1(v, v)\} + \{FQ_2(v, s)\} = \{g\} \\ \{s\} = [D][B(q_0)]\{v\} + \{fq(v, v)\} \end{cases} \quad [5] $$

où $[M]$ est la matrice masse classique, $[K_t]$ est la matrice tangente évaluée en $\{q_0\}$, condition initiale, $\{FQ_1(v, v)\}$ et $\{FQ_2(v, s)\}$ représentent respectivement les parties quadratiques en (v,v) et en (v,s), $\{g\} = \{F\} - \sum_e \int_{\Omega^e} [B(q_0)] \{S_0\} d\Omega^e$ et $\{fq(v, v)\} = \frac{1}{2} [D] [B_{NL}(v)] \{v\}$.

Le principe de l'homotopie consiste à modifier (5), en introduisant un paramètre $\varepsilon$ dans le terme non linéaire :

$$ [M]\{\ddot{v}\} + [K_t]\{v\} + \varepsilon(\{FQ_1(v, v)\} + \{FQ_2(v, s)\}) = \{g\} \quad [6] $$
Ainsi, la solution {$v(t, \epsilon)$} qui dépend du paramètre $\epsilon$ et du temps $t$ de (6) passe d'une manière continue de la solution {$v(t, \epsilon = 0)$} à la solution {$v(t, \epsilon = 1)$} du problème (5). {$v(t, \epsilon = 0)$} est la solution d'un problème linéaire obtenu à partir de (6) en faisant $\epsilon$ égal à zéro. La solution cherchée est obtenue en écrivant {$q(t)$} = {$q_0$} + {$v(t, \epsilon = 1)$}.

**ii- Discrétisation temporelle :** pour la résolution du problème instationnaire non linéaire (6), on utilise ici le schéma implicite de Newmark, largement utilisé dans la résolution des problèmes dynamiques [NEW 59]. Ce schéma donne l'expression de la vitesse et de l'accélération nodales à l'instant $(n + 1)\Delta t$ [GER 93] :

$$ egin{cases} \{\dot{v}^{n+1}\} = \{\dot{v}^n\} + \Delta t[(1-\beta)\{\ddot{v}^n\} + \beta\{\ddot{v}^{n+1}\}] \\ \{\ddot{v}^{n+1}\} = a_0(\{v^{n+1}\} - \{v^n\}) - a_1\{\dot{v}^n\} - a_2\{\ddot{v}^n\} \end{cases} \quad [7] $$

où $a_0 = \frac{1}{\alpha\Delta t^2}$, $a_1 = \frac{1}{\alpha\Delta t}$ et $a_2 = \frac{1}{2\alpha} - 1$, $\alpha$ et $\beta$ sont les constantes de Newmark. En utilisant (7), le problème (6) discrétisé en temps s'écrit :

$$ [K_T^0]\{v^{n+1}\} + \varepsilon(\{FQ_1(v^{n+1}, v^{n+1})\} + \{FQ_2(v^{n+1}, s^{n+1})\}) = \{g^{n+1}\} + \{H^n\} \quad [8] $$

$$ \text{où } [K_T^0] = a_0[M] + [K_t] \text{ et } \{H^n\} = [M](a_0\{v^n\} + a_1\{\dot{v}^n\} + a_2\{\ddot{v}^n\}). $$

**iii- Méthode de perturbation :** pour résoudre le problème non linéaire (8), on propose ici deux choix de développements asymptotiques :

1er choix : on suppose que les quantités ($\{v^{n+1}\}$, $\{s^{n+1}\}$, $\{H^n\}$) s'écrivent sous la forme d'un développement en séries entières tronquées à l'ordre p par rapport à $\epsilon$ comme dans [COC 00]. Ainsi les termes $\{v_k^{n+1}\}$, $\{s_k^{n+1}\}$ et $\{H_k^n\}$ sont solutions de la séquence de problèmes linéaires suivants :

$$ [K_T^0]\{v_0^{n+1}\} = \{g^{n+1}\} + \{H_0^n\} \quad ; \quad [K_T^0]\{v_k^{n+1}\} = \{H_k^n\} - \{FQ_k^{nl}\} \quad [9] $$

$$ \text{avec } \{FQ_k^{nl}\} = \sum_{r=0}^{k-1} \{FQ_1(v_r^{n+1}, v_{k-r-1}^{n+1})\} + \sum_{r=0}^{k-1} \{FQ_2(v_r^{n+1}, s_{k-r-1}^{n+1})\}. $$

2ème choix : ici, on développe seulement les inconnues $\{v^{n+1}\}$ et $\{s^{n+1}\}$ en série par rapport à $\epsilon$. L'introduction de ces développements dans (8) donne la séquence de problèmes linéaires :

$$ [K_T^0]\{v_0^{n+1}\} = \{g^{n+1}\} + \{H^n\} \quad ; \quad [K_T^0]\{v_k^{n+1}\} = -\{FQ_k^{nl}\} \quad [10] $$

La différence entre les deux choix réside dans le second membre. En effet, le terme $\{H^n\}$ à l'instant précédent apparaît complètement dans le second membre du deuxième choix à l'ordre zéro. Dans la formule similaire du premier choix, il n'y avait que le premier terme $\{H_0^n\}$ du développement en série de $\{H^n\}$. Ainsi, si la série $\{H^n\} = \{H_0^n\} + \varepsilon\{H_1^n\} + \dots + \varepsilon\{H_p^n\}$ converge difficilement, le premier choix récupère cette convergence lente à l'instant suivant.

Les représentations polynomiales ont généralement un domaine de validité fini, pour décrire toute la branche solution, on applique une méthode de continuation pas à pas en définissant un nouveau point de départ pour chaque nouveau pas. Le domaine
de validité de la solution est défini par le critère suivant $\epsilon_{ms} = (\kappa \frac{\|\{v_1\}\|}{\|\{v_p\}\|})^{\frac{1}{p-1}}$ où $\kappa$ est un paramètre qu'on suppose petit.

Même si on détermine un grand nombre de termes de la série, on peut espérer reconstruire la fonction même au-delà du rayon de convergence de la série en introduisant les approximants de Padé [PAD 92]. Dans ce travail, seules les fractions rationnelles ayant le même dénominateur [BRA 97], [NAJ 98] seront considérées. Pour une telle représentation, la fin du pas de continuation est caractérisée par le critère $\frac{\|\{P_p(\epsilon_{mp})\} - \{P_{p-1}(\epsilon_{mp})\}\|}{\|\{P_p(\epsilon_{mp})\} - \{v_0\}\|} = \delta$ où $\delta$ est un paramètre supposé petit, $\{P_p(\epsilon_{mp})\}$ et $\{P_{p-1}(\epsilon_{mp})\}$ correspondent aux Padé des séries respectivement aux ordres successifs $p$ et $p-1$ [ELH 00].

## 3.2. Linéarisation par rapport à l'instant $t = n\Delta t$ et introduction d'un préconditionneur

Cet algorithme est basé sur les mêmes étapes de l'algorithme précédent. Après discrétisation en temps, on utilise une homotopie différente avec introduction d'un préconditionneur $[K^*]$.

### i- Procédure de discrétisation temporelle : en utilisant le schéma implicite de Newmark, on peut écrire le problème (2) sous la forme :

$$ \begin{aligned} & \sum_e \int_{\Omega^e} a_0 \rho^t[N][N] d\Omega^e(\{q^{n+1}\} - \{q^n\}) + \sum_e \int_{\Omega^e} \rho^t[B(q^{n+1})]\{S^{n+1}\} d\Omega^e = \\ & \{F^{n+1}\} + \sum_e \int_{\Omega^e} \rho^t[N][N] d\Omega^e(a_1\{q^n\} + a_2\{q^{n+1}\}) \end{aligned} \quad [11] $$

### ii- Transformation homotopique : on introduit dans (11) et dans la formule de la loi de comportement (2 - b), le changement de variable suivant :

$$ \left\{ \begin{aligned} \{q^{n+1}\} &= \{q^n\} + \{\Delta q\} \\ \{S^{n+1}\} &= \{S^n\} + \{\Delta S\} \end{aligned} \right. \quad [12] $$

où $\{\Delta q\}$ et $\{\Delta S\}$ sont des incrémentes qui représentent les nouvelles inconnues. En remplaçant dans la partie linéaire de l'équation obtenue la nouvelle expression de la contrainte $\{\Delta S\}$, le problème verifié par $\{\Delta q\}$ et $\{\Delta S\}$ s'écrit :

$$ \left\{ \begin{aligned} & [K_T^n]\{\Delta q\} + \{FQ_1(\Delta q, \Delta q)\} + \{FQ_2(\Delta q, \Delta S)\} = \{R^{n+1}\} \\ & \{\Delta S\} = [D]([B(q^n)] + \frac{1}{2}[B_{NL}(\Delta q)])\{\Delta q\} \end{aligned} \right. \quad [13] $$

où $[K_T^n]$ est la même matrice que dans (9) ou (10) mais évaluée à l'instant $n\Delta t$, $\{FQ_1(\Delta q, \Delta q)\}$ et $\{FQ_2(\Delta q, \Delta S)\}$ sont les mêmes que dans le paragraphe précédent et $\{R^{n+1}\} = \{F^{n+1}\} + \sum_e \int_{\Omega^e} \rho^t[N][N] d\Omega^e(a_1\{q^n\} + a_2\{q^{n+1}\}) - \sum_e \int_{\Omega^e} \rho^t[B(q^n)]\{S^n\} d\Omega^e$.

Pour éviter la décomposition de la matrice tangente $[K_T^n]$ à chaque pas de temps, on utilise une transformation homotopique qui consiste à introduire une matrice arbitraire
$[K^*]$ (un préconditionneur) et un paramètre $\epsilon$. L'équation (13) est écrite alors sous la forme :

$$[K^*]\{\Delta q\} + \epsilon([K_T^n] - [K^*])\{\Delta q\} + \{FQ_1(\Delta q, \Delta q)\} + \{FQ_2(\Delta q, \Delta S)\} = \epsilon\{R^{n+1}\} \quad [14]$$

Dans les applications numériques, nous allons choisir le préconditionneur $[K^*]$ égal à la matrice tangente évaluée en temps au début de chaque grand pas de continuation.

### iii- Méthode de perturbation :
on cherche la solution de (14) sous la forme d'une série entière tronquée à un ordre p par rapport à $\epsilon$. L'injection de ce développement dans (14) conduit à la séquence de problèmes linéaires ayant tous la même matrice $[K^*]$ :

$$[K^*]\{\Delta q_1\} = \{R^{n+1}\} \quad ; \quad [K^*]\{\Delta q_k\} = [[K^*] - [K_T^n]]\{\Delta q_{k-1}\} - \{FQ_k^{nl}\} \quad [15]$$

avec $\{FQ_k^{nl}\} = \sum_{r=1}^{k-1}\{FQ_1(\Delta q_r, \Delta q_{k-r})\} - \sum_{r=1}^{k-1}\{FQ_2(\Delta q_r, \Delta S_{k-r})\}$. Finallement la solution du problème de départ est donnée par : $\{q^{n+1}\} = \{q^n\} + \{\Delta q_1\} + \dots + \{\Delta q_p\}$.

## 4. Application numérique

Pour tester la performance des algorithmes proposés, on étudie la vibration non linéaire d'une poutre élastique 2D et d'une plaque carrée élastique. La solution obtenue est comparée avec celle obtenue par la méthode classique utilisant le schéma de Newmark couplé avec l'algorithme de Newton-Raphson (NR), ainsi qu'une méthode itérative semblable à celle utilisée dans divers secteurs de la physique [FOR 78], elle est décrite dans [JAM 02]. Dans ces deux méthodes itératives, la convergence de l'algorithme est évaluée en considérant que l'écart relatif entre deux itérations est inférieur à une tolérance donnée. Dans les deux exemples la matrice $[K^*]$ a la même forme que $[K_T^0]$ définie aux conditions initiales. Dans les calculs présentés, on fixe l'ordre de troncature des séries et donc le nombre de seconds membres à calculer (NSM), ainsi l'efficacité des algorithmes est mesurée d'une part par le nombre de matrices à inverser (NIM) et, d'autre part, par la petitesse de la norme du résidu $||\{R\}||$ (équation 3).

### 4.1. Vibration non linéaire forcée d'une poutre élastique 2D

Les caractéristiques mécaniques et géométriques de la poutre considérée sont : module d'Young $E = 10^5$ MPa, coefficient de Poisson $\nu = 0,3$, masse volumique $\rho = 10^{-4}$ Kg/m$^3$, longueur $L = 200$ mm et largeur $l = 10$ mm. La poutre étudiée est encastrée d'un côté et chargée de l'autre côté d'une force concentrique $\lambda(t)F$ (voir figure 1) où $\lambda(t)$ est un facteur de chargement dépendant du temps $t$ comme indiqué sur la figure 1 et $F = 200$ N est l'intensité de la force appliquée au nœud (100, 0). La
structure est discrétisée en 12 éléments Q8. Les constantes de Newmark et le pas du temps choisis sont : $\alpha = 0,25$, $\beta = 0,5$ et $\Delta t = 2 \cdot 10^{-4}$s. La solution est cherchée dans l'intervalle du temps $[0, 4s]$ pour une condition initiale nulle.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Poutre élastique 2D et facteur de chargement en fonction du temps. (b) Déplacement $U_y$ suivant l'axe oy en fonction du temps</figcaption></figure>

Les résultats numériques obtenus par les algorithmes proposés, LP1 (Linéarisation par rapport à l'instant initial $1^{er}$ choix), LP2 (Linéarisation par rapport à l'instant initial $2^{ème}$ choix) et LPK\* (Linéarisation par rapport à l'instant $t = n\Delta t$ et introduction de $[K^*]$), sont reportés sur le tableau (1) et la figure (1-b). Le tableau (1) illustre une comparaison entre les solutions série et Padé obtenues par les algorithmes proposés tronquées aux différents ordres : 6, 10, 12 et 15 pour $\Delta t = 2 \cdot 10^{-4}$. La petite valeur de $\Delta(t)$ est due au fait que le problème étudié est fortement non linéaire. En effet, l'amplitude du déplacement d'après la figure (1-b) est presque de l'ordre de la largeur de la poutre étudiée.

<table><caption>**Tableau 1.** Comparaison entre les trois algorithmes LP1, LP2 et LPK\* (avec série et Padé) pour les ordres 6, 10, 12 et 15 et pour $\Delta t = 0,0002$s</caption><thead><tr><th rowspan="2">p</th><th colspan="2">LP1</th><th colspan="2">LP2</th><th colspan="2">LPK*</th></tr><tr><th colspan="2">Série</th><th colspan="2">Padé</th><th colspan="2">Série</th></tr></thead><tbody><tr><td>6</td><td>21 10<sup>-2</sup></td><td>21 10<sup>-1</sup></td><td>31 10<sup>-2,5</sup></td><td>31 10<sup>-2,5</sup></td><td>11 10<sup>-1,6</sup></td><td>11 10<sup>-1,6</sup></td></tr><tr><td>10</td><td>11 10<sup>-1,5</sup></td><td>11 10<sup>-2</sup></td><td>11 10<sup>-1,5</sup></td><td>11 10<sup>-2</sup></td><td>1 10<sup>-2,7</sup></td><td>1 10<sup>-4,5</sup></td></tr><tr><td>12</td><td>11 10<sup>-1,5</sup></td><td>11 10<sup>-2</sup></td><td>1 10<sup>-2</sup></td><td>1 10<sup>-3,2</sup></td><td>1 10<sup>-3,7</sup></td><td>1 10<sup>-6,2</sup></td></tr><tr><td>15</td><td>6 10<sup>-1,5</sup></td><td>6 10<sup>-2</sup></td><td>1 10<sup>-3,5</sup></td><td>1 10<sup>-5,7</sup></td><td>1 10<sup>-5,4</sup></td><td>1 10<sup>-7</sup></td></tr></tbody></table>

Ce tableau montre, qu'au-delà de l'ordre 10, l'algorithme LP1 utilise plus de décompositions de matrices que les deux autres. On remarque que l'introduction des approximants de Padé conduit à une bonne qualité de la solution quand l'ordre de troncature augmente surtout pour la solution calculée par LP2 et LPK\*.

L'algorithme couplant le schéma de Newmark et la méthode de Newton-Raphson a eu besoin de 40 000 inversions de matrice pour calculer la solution sur $[0, 4]$ et LP1 a
nécessité 11 inversions de matrices tandis que LP2 et LPK* n'ont eu besoin que d'une seule inversion de matrice pour l'ordre 12.

## 4.2. Vibration non linéaire forcée d'une plaque élastique

On considère ici une plaque mince carrée, de côté $l = 100mm$, d'épaisseur $1mm$, encastrée d'un côté et chargée du côté opposé par une force concentrée $\lambda(t)F$ au nœud (100,0,0) avec $F = 200N$ et $\lambda(t)$ est donné dans la figure (1-a). Elle est caractérisée par $E = 2,0810^{+6}MPa$, $\nu = 0,3$ et $\rho = 2.10^{-3}Kg/m^3$. La plaque considérée est discrétisée en 162 éléments de type DKT18 (voir figure 2 - a). On a choisi pour l'intégration temporelle : $\alpha = 0,25$, $\beta = 0,5$ et $\Delta t = 2.10^{-4}s$. Ici aussi, à cause de la forte non-linéarité, les pas de temps doivent être choisis assez petits.

La solution est cherchée sur l'intervalle temporel $[0, 4s]$. Le tableau (2) illustre une comparaison des solutions série et Padé tronquées aux ordres 7, 9 et 11, obtenues par les algorithmes proposés. D'après ce tableau, on remarque que la solution obtenue par LP2 et LPK*, en introduisant les approximants de Padé, a une qualité meilleure que celle obtenue par la série. La solution obtenue, voir figure (2-b) sur l'intervalle $[0, 4s]$, a nécessité 40 000 inversions pour l'algorithme couplant le schéma de Newmark avec la méthode de Newton-Raphson et 17 seulement pour l'algorithme proposé LP1 tandis que les algorithmes LP2 et LPK* n'ont exigé qu'une seule inversion. On voit aussi qu'on peut améliorer la qualité de la solution (i.e. diminuer la norme du résidu) uniquement en augmentant l'ordre de troncature. Il semble qu'on puisse atteindre une qualité de solution aussi bonne que l'on désire, en particulier avec l'algorithme LPK*.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) Plaque carrée encastrée soumise à un chargement dynamique. (b) Flèche en fonction du temps</figcaption></figure>

Les résultats obtenus par les solveurs proposés, dans le cas de l'exemple de la plaque élastique, sont comparés à ceux calculés par l'algorithme couplant le schéma de Newmark avec la méthode de Newton Raphson et la méthode itérative (tableau (3)). Il est à souligner que la méthode itérative utilisée dans ce travail est une amélioration de la version classique. La différence réside dans l'actualisation de la matrice
<table><caption>**Tableau 2.** Comparaison entre les trois algorithmes LP1, LP2 et LPK* (avec série et Padé) pour les ordres 7, 9 et 11 et pour $\Delta t = 0,0002s$</caption><thead><tr><th rowspan="2">p</th><th colspan="2">LP1</th><th colspan="2">LP2</th><th colspan="2">LPK*</th></tr><tr><th>Serie</th><th>Padé</th><th>Serie</th><th>Padé</th><th>Serie</th><th>Padé</th></tr></thead><tbody><tr><td>7</td><td>40 10<sup>-2,5</sup></td><td>40 10<sup>-2</sup></td><td>2 10<sup>-1,5</sup></td><td>2 10<sup>-1,3</sup></td><td>1 10<sup>-2,6</sup></td><td>1 10<sup>-2,8</sup></td></tr><tr><td>9</td><td>23 10<sup>-2,2</sup></td><td>23 10<sup>-2</sup></td><td>1 10<sup>-2,2</sup></td><td>1 10<sup>-2,8</sup></td><td>1 10<sup>-4</sup></td><td>1 10<sup>-5,2</sup></td></tr><tr><td>11</td><td>17 10<sup>-2</sup></td><td>17 10<sup>-2</sup></td><td>1 10<sup>-3,6</sup></td><td>1 10<sup>-4,8</sup></td><td>1 10<sup>-5,6</sup></td><td>1 10<sup>-7</sup></td></tr></tbody></table>

d'itération dès qu'il y a divergence. Une comparaison de cette version améliorée avec l'algorithme couplant le schéma de Newmark avec la méthode de Newton Raphson est donnée dans le tableau 3 pour un pas de temps $\Delta t = 0,0002s$.

<table><caption>**Tableau 3.** Comparaison de la méthode itérative améliorée et la méthode de Newton Raphson pour $\Delta t = 0,0002s$, NSM : nombre de seconds membres, NIM : nombre d'inversions de matrice et $||{R}||$ : norme euclidienne du vecteur résidu.</caption><thead><tr><th colspan="4">Méthode itérative</th><th colspan="3">Méthode de Newton Raphson</th></tr><tr><th>Tolérance</th><th>NSM</th><th>NIM</th><th>||{R}||</th><th>Tolérance</th><th>NIM</th><th>||{R}||</th></tr></thead><tbody><tr><td>10<sup>-7</sup></td><td>62260</td><td>26</td><td>10<sup>-4</sup></td><td>10<sup>-4</sup></td><td>40000</td><td>10<sup>-6,5</sup></td></tr><tr><td>10<sup>-9</sup></td><td>84469</td><td>19</td><td>10<sup>-5</sup></td><td>10<sup>-5</sup></td><td>40000</td><td>10<sup>-6,5</sup></td></tr><tr><td>10<sup>-10</sup></td><td>113472</td><td>5</td><td>10<sup>-5</sup></td><td>10<sup>-6</sup></td><td>42690</td><td>10<sup>-6,5</sup></td></tr><tr><td>10<sup>-11</sup></td><td>130088</td><td>3</td><td>10<sup>-5,5</sup></td><td>10<sup>-7</sup></td><td>42690</td><td>10<sup>-7</sup></td></tr></tbody></table>

On voit, sur le tableau 3, que le nombre d'inversions de matrice (NIM) est beaucoup plus petit dans le cas de la méthode itérative comparativement à celui de l'algorithme couplant le schéma de Newmark avec la méthode de Newton Raphson. Quant aux qualités des deux méthodes ($||{R}||$), elles sont presque du même ordre de grandeur. La méthode itérative est un bon compromis, peu d'inversions de matrices et un nombre de seconds membres raisonnables (3 ou 4 itérations par pas avec une tolérance de $10^{-7}$ ou $10^{-9}$).

Enfin, pour mettre en lumière l'efficacité des algorithmes itératifs, nous avons repris ce calcul sur un intervalle de temps plus grand $[0, 6s]$ et avec de plus grands pas de temps ($2 \cdot 10^{-2}s$, $4 \cdot 10^{-2}s$, $8 \cdot 10^{-2}s$) (voir tableau (4)). Ces pas de temps sont certes un peu trop grands pour donner une solution très précise, mais cela permet d'avoir une première approximation de la solution dynamique avec assez peu de pas de temps. Les paramètres définissant le calcul (tolérance, ordre, résidu admissible) ont été choisis de telle sorte que le résidu maximal soit le même dans les trois cas (environ $10^{-4}$).
On constate que l'algorithm d'ordre élevé (LP2, ordre 14) converge toujours, même pour le pas de temps le plus grand, alors que l'algorithm couplant le schéma de Newmark avec la méthode de Newton-Raphson diverge dans ce cas. Comme on pouvait s'y attendre, la méthode itérative ne converge que pour le pas de temps le plus petit. Mais dans ce cas, elle demande plus de 10 fois moins d'inversions de matrices avec un nombre de seconds membres raisonnables. En résumé, les algorithmes d'ordre élevé sont des solveurs non linéaires les plus performants, ils sont meilleurs que l'algorithm couplant le schéma de Newmark avec la méthode de Newton-Raphson (NR) qui est à son tour meilleure que la méthode itérative (MI). Toutefois cette dernière, qui est simple et peu coûteuse, méritera d'être davantage utilisée en calcul des structures.

<table><caption>**Tableau 4.** Comparaison des trois algorithmes : LP2, NR et MI pour différents pas de temps sur l'intervalle de temps [0, 6s].</caption><tr><td>Δt(s)</td><td>8.10<sup>-2</sup></td><td>4.10<sup>-2</sup></td><td>2.10<sup>-2</sup></td></tr><tr><td>NR</td><td>diverge</td><td>779 inversions</td><td>1200 inversions</td></tr><tr><td>MI</td><td>diverge</td><td>diverge</td><td>108 inversions<br>4908 SM</td></tr><tr><td>LP2 ordre 14</td><td>9 inversions<br>1125 SM</td><td>9 inversions<br>2250 SM</td><td>2 inversions<br>4500 SM</td></tr></table>

## 5. Conclusion

Nous avons proposé une famille d'algorithmes implicites d'ordre élevé pour la résolution des problèmes de la dynamique non linéaire des structures. Ces algorithmes diffèrent entre eux selon l'ordre de l'application de la transformation homotopique, de la discrétisation en temps et de la méthode de perturbation. Le préconditionneur a introduit beaucoup de souplesse à l'algorithm LPK*. Le résultat le plus intéressant est la capacité d'appliquer un schéma implicite LPK* et LP2 avec peu de matrices à triangulariser. L'introduction des approximants de Padé a permis une amélioration de la qualité de la solution.

## 6. Bibliographie

[ARG 91] ARGYRIS J., MLEJNEK H., *Dynamics of Structures*, Computational Mechanics, 1991.

[BAT 96] BATHE K. J., *Finites elements procedures*, Prentice-Hall, Engwood Cliffs, New Jersey 07632, 1996.

[BRA 97] BRAIKAT B., DAMIL N., POTIER-FERRY M., « Méthodes asymptotiques numériques pour la plasticité », *Revue Européenne des Eléments finis*, vol. 6, n° 3, 1997, p. 337-357.
[BRA 02] BRAIKAT B., JAMAL M., DAMIL N., POTIER-FERRY M., « Algorithmes d'intégration temporelle implicites couplés avec des résoluteurs d'ordre élevé », *Revue Européenne des Eléments finis*, vol. 11, 2002, p. 749-772.

[BRA 04] BRAIKAT B., JAMAL M., DAMIL N., « Utilisation des techniques de la méthode asymptotique numérique pour la résolution des problèmes instationnaires non linéaires », *Revue Européenne des Eléments finis*, vol. 13, 2004, p. 119-139.

[COC 00] COCHELIN B., COMPAIN C., « An asymptotic numerical method for nonlinear transient dynamics », *Revue Européenne des Eléments finis*, vol. 9, 2000, p. 113-128.

[CRI 91] CRISFIELD M., *Non-linear finite element analysis of solids and structures*, vol. 1, John Wiley & Sons, 1991.

[CRI 97] CRISFIELD M. A., *Non-linear finite element analysis of solids and structures*, vol. 2, John Wiley, Chichester, 1997.

[ELH 00] ELHAGE-HUSSEIN A., POTIER-FERRY M., DAMIL N., « A numerical continuation method based on Padé Approximants », *International Journal of Solids and Structures*, vol. 37, 2000, p. 6981-7001.

[FOR 78] FORNBERG B., WHITHAM G. B., « A numerical and theoretical study of certain nonlinear wave phenomena », *Philosophical Transactions of the Royal Society of London, A Mathematical and Physical Sciences*, vol. 289, 1978, p. 373-404.

[GER 93] GERADIN M., RIXEN D., *Théorie des vibrations. Application à la dynamique des structures*, Masson, 1993.

[JAM 02] JAMAL M., BRAIKAT B., BOUTMIR S., DAMIL N., POTIER-FERRY M., « A high order implicit algorithm for solving instationnary non-linear problems », *Computational Mechanics*, vol. 28, 2002, p. 375-380.

[NAJ 98] NAJAH A., COCHELIN B., DAMIL N., POTIER-FERRY M., « A critical review of asymptotic numerical methods », *Archives of Computational Methods in Engineering*, vol. 5, 1998, p. 31-50.

[NEW 59] NEWMARK N. M., « A method of computation for structural dynamics », *The Engineering Mechanics Division, Proceedings of ASCE*, vol. 85, 1959, p. 67-94.

[PAD 92] PADÉ H., « Sur la représentation approchée d'une fonction par des fractions rationnelles », *Ann. de l'Ecole Normale Sup., série 3*, vol. 9, 1892, p. 3-93.

[ZIE 87] ZIENKIEWICZ O. C., TAYLOR R. L., *The finite element method. Solid and Fluid Mechanics and non linearity*, vol. 2, McGraw-Hill, Book Company, 1987.