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

HAL
open science

# Analyse multi-échelle du flambage de grandes structures

Philippe Cresta, Christian Rey, Olivier Allix, Stéphane Guinard

► To cite this version:

Philippe Cresta, Christian Rey, Olivier Allix, Stéphane Guinard. Analyse multi-échelle du flambage de grandes structures. 17e Congrès Français de Mécanique, Aug 2005, Troyes, France. hal-04089442

HAL Id: hal-04089442

https://hal.science/hal-04089442

Submitted on 4 May 2023

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Analyse multi-échelle du flambage de grandes structures

Philippe Cresta¹, Christian Rey², Olivier Allix² & Stéphane Guinard¹

¹EADS Centre Commun de Recherche (CCR), Département Ingénierie des Structures (DCR/SE)
12 rue Pasteur, BP 76
92152 SURESNES Cedex, France

²Laboratoire de Mécanique et Technologie (LMT), ENS Cachan
61 avenue du Président Wilson
94235 CACHAN Cedex, France

## Résumé :

Un challenge industriel d'aujourd'hui est d'être à même d'effectuer des calculs de post-flambage de grandes structures de type fuselage à une échelle fine (~10⁷-10⁸ d.d.l.), le but étant à terme d'éviter des essais à l'échelle 1 jusqu'à rupture. Un calcul direct étant hors de portée, ce travail explore les possibilités liées aux approches multi-échelles et parallèles dans le contexte des non-linéarités géométriques. La formulation non-linéaire retenue est tout d'abord introduite, et les méthodes classiques de résolution rappelées. Dans un deuxième temps, les méthodes spécifiques que nous avons développées sont présentées. Le potentiel des différentes stratégies numériques est ensuite comparé sur un cas test et les résultats sont analysés.

## Abstract :

Performing non-linear analysis of large structures at a fine scale (~10⁷-10⁸ d.o.f.) is one of the industrial challenges of our times, the goal herein being to avoid testing large structures up to failure. Since computing that kind of problem by a classical technique remains inaccessible at the moment, we explore the possibilities related to multiscale and parallel approaches in the context of geometrically non-linear problems. The retained non-linear formulation is first introduced as well as the classical computational methods used to solve such problems. In a second part, the specific methods we developed are presented. The potential of the various numerical strategies is then compared on a simplified study case and the results are analyzed.

## Mots clefs :

non-linéaire géométrique ; post-flambage ; décomposition de domaine ; calcul parallèle.

## 1 Problématique industrielle

Les structures aéronautiques sont constituées d'éléments élancés susceptibles de flamber au-delà d'une charge critique. Lors des essais structuraux de certification, réalisés sur des fuselages complets d'aéronefs, on observe localement l'apparition de zones de flambage de la peau, formant des cloques entre les raidisseurs. Avec l'augmentation du chargement, ces zones non-linéaires peuvent s'étendre et provoquer des redistributions de contraintes dans la structure. Pour les charges de service couramment rencontrées, ces phénomènes sont réversibles, le matériau restant dans le domaine élastique. Cependant, ils peuvent provoquer des concentrations de contrainte au niveau des bases des raidisseurs et être à l'origine de décollements locaux menant à la rupture globale. Le calcul de telles structures (~10⁷-10⁸ d.d.l.) par des techniques classiques est aujourd'hui hors de portée du fait des limitations en puissance de calcul et mémoire des stations de travail traditionnelles. Les ingénieurs des bureaux d'études réalisent cependant des calculs approximatifs au moyen d'analyses descendantes, à savoir des calculs globaux grossiers et linéaires suivis d'analyses fines non-linéaire sur des zones d'intérêt local. Ces méthodes n'offrent toutefois qu'un dialogue limité entre les échelles introduites et ne peuvent traiter que les non-linéarités localisées n'ayant pas
d'influence sur la réponse globale.

Les progrès de l'informatique, avec l'arrivée des machines parallèles et des clusters de PC, permettent
aujourd'hui de partager les coûts de calcul et les besoins en mémoire entre plusieurs processeurs/machines.
Des problèmes de plus en plus importants deviennent accessibles à conditions de mettre en place des
méthodes parallèles adaptées. Aussi cette dernière décennie a vu le développement de solveurs linéaires
basés sur les méthodes de décomposition de domaines et des solveurs itératifs de Krylov (BDD, FETI,...)
[1] [2] [3]. Un premier pas vers le non-linéaire est de les exploiter dans le cadre d'un schéma de Newton [2].
Les résultats obtenus montrent toutefois des limitations en termes de convergence pour les problèmes où
les effets non-linéaires sont non-équitablement répartis [4].

2 Formulation et méthodes classiques de résolution

2.1 Formulation non-linéaire

Cette étude concerne les non-linéarités d'origine géométrique, dues aux grands déplacements et grandes
rotations de structures minces constituées par des assemblages de plaques ou coques. Le comportement du
matériau est supposé rester linéaire, les déformations restant petites. Une formulation de type
corotationnelle a alors été retenue car bien adaptée aux hypothèses de grands déplacements et de petites
déformations des éléments constitutifs [6]. Dans cette phase de développements numériques, des
assemblages de poutres sont considérés afin de tester les méthodes numériques sur des problèmes
représentatifs de petite taille. Les outils retenus seront ultérieurement déployés sur le cas simplifié d'un essai
statique de tronçon Airbus.

<figure><img src="image_3.png" /><figcaption>FIG. 1 –Définition des différentes configurations pour un élément de poutre.</figcaption></figure>

Les déplacements $\tilde{v}_d$ et efforts $\tilde{f}$ sont exprimés dans un repère tournant attaché aux éléments (figure 1). Ils sont reliés par l'opérateur de Hooke élémentaire classique de l'élasticité linéaire, $\tilde{K}_e$, exprimé dans le même référentiel. Les efforts sont ensuite équilibrés par un projecteur $P$ puis ramenés dans la configuration initiale par la rotation $R$ :

$$
f = R^T P^T \tilde{K}_e \tilde{v}_d
$$

Par dérivation de cette expression, on obtient l'opérateur tangent cohérent $K_T$ qui peut s'exprimer sous la
forme suivante [6]:

$$
K_T = R^T (P^T \tilde{K}_e HP - F_{mn}G - G^T F_n^T P)R
$$

où $H, F_{mn}, G, F_n$ sont des termes explicités dans [6] représentant les raideurs géométriques de la structure
précontrainte. Cette formulation peut être mise en œuvre pour tous les types d'éléments de manière assez
simple et indépendante. Elle permet en outre de s'affranchir des problèmes rencontrés par les formulations
lagrangiennes classiques dans le cas de grandes rotations.

2.2 Algorithmes de résolution parallèle classiques

2.2.1 Linéarisation par méthode de Newton-Raphson

Le problème non-linéaire obtenu est résolu en général par une méthode de Newton-Raphson ou une de ses variantes (méthode sécante, etc.). Dans tous les cas, le système est linéarisé autour d'une position connue et on se ramène à la résolution itérative de systèmes linéaires du type :

$$
K_T \Delta u = r
$$

où $K_T$ est l'opérateur tangent et $r$ le résidu. La convergence de ces algorithmes est quadratique mais dépend
de la taille du pas de chargement. On met donc en œuvre cette méthode dans le cadre d'un processus incrémental. Pour une méthode de Newton-Raphson classique avec pilotage du chargement en effort, on adapte la taille du pas en fonction de la convergence du pas précédent. Dans notre cas, en cas de non-convergence, le pas est divisé par 2 ; après plusieurs convergences successives, sa taille est augmentée.

### 2.2.2 Méthodes de décomposition de domaines

Afin de pouvoir traiter des problèmes à grands nombres de ddl, les problèmes linéaires obtenus par linéarisation du problème non-linéaire initial peuvent être résolus par des méthodes parallèles. La plupart de ces méthodes reposent sur une décomposition de la structure en sous-domaines qui peuvent alors être répartis sur différents processeurs [3]. Les compléments de Schur sur les interfaces des opérateurs tangents forment un problème linéaire global qui est en général résolu par un solveur itératif de Krylov, ne nécessitant pas l'assemblage du problème. Une implémentation avec la méthode FETI est proposée dans [2], les inconnues aux interfaces étant les flux d'effort.

Une méthode classique de décomposition de domaines a été mise en œuvre pour servir de méthode de référence pour cette étude. Nous avons opté ici pour une version primale en choisissant comme inconnues les déplacements aux interfaces. La résolution d'un système linéaire se résume aux étapes suivantes (l'indice *i* indique les inconnues internes, *b* les inconnues sur les frontières):

1. *Condensation, par sous-structure, des matrices tangentes et résidus sur les interfaces :*

$$S_T^{(i)} = K_{Tbb}^{(i)} - K_{Tbi}^{(i)} K_{Tii}^{(i)-1} K_{Tib}^{(i)}$$

$$b^{(i)} = r_b^{(i)} - K_{Tbi}^{(i)} K_{Tii}^{(i)-1} r_i^{(i)}$$

2. *Résolution par méthode itérative du problème condensé aux interfaces :*

$$S_T \Delta u_b = b$$

3. *Relocalisation des déplacements internes par sous-structure :*

$$\Delta u_i^{(i)} = K_{Tii}^{(i)-1} (f_i^{(i)} - K_{Tib}^{(i)} \Delta u_b^{(i)})$$

Toutes les étapes sont parallélisables et permettent de partager le coût du calcul entre plusieurs processeurs. L'étape de résolution (2) exige cependant des échanges de données qui peuvent pénaliser les performances globales du calcul. Des améliorations des solveurs itératifs de Krylov ont été proposées dans [7] pour réduire ces échanges, en contrepartie d'un stockage plus important d'informations.

### 2.2.3 Application sur une structure test

La figure 2 présente la structure test retenue pour cette étude. Il s'agit d'une structure en treillis constituée de poutres de type Navier-Bernoulli que l'on vient mettre en flexion globale. Celle-ci est décomposée en 10 sous-structures comportant chacune quelques centaines de degrés de liberté intérieurs et 12 d.d.l d'interface. On peut observer sur la figure une non-linéarité globale assez douce, due aux grands déplacements de la structure, ainsi qu'une non-linéarité de flambage plus forte et localisée sur les éléments de poutre à proximité de la zone d'encastrement. Cette structure présente donc un comportement similaire aux cas aéronautiques qui nous intéressent et présente l'avantage de la simplicité pour développer et tester différentes méthodes avant de passer à des structures plus proches des cas industriels.

<figure><img src="image_4.png" /><figcaption>FIG. 2 – Structure de type treillis étudiée. Configuration initiale et déformée.</figcaption></figure>
La convergence de l'algorithme classique de Newton pour ce problème s'avère très moyenne, le pas de chargement devant être réduit pour passer l'instabilité de flambage local. La décomposition de domaine n'apporte aucune modification au résultat de convergence. Elle ne fait que paralléliser une méthode classique en répartissant l'effort de calcul uniformément sur l'ensemble de la structure. On ne profite donc pas par cette méthode du caractère local des non-linéarités de flambage [4].

# 3 Approche primale à deux niveaux

## 3.1 Présentation

Partant des constatations faites après analyse des méthodes de type Newton-Krylov-Schur, une méthode nouvelle permettant la prise en compte du caractère local des non-linéarités de flambage est proposée. Après chaque résolution globale du problème condensé aux frontières, on effectue une résolution non-linéaire par sous-structure, en imposant les déplacements aux interfaces. La seule modification par rapport à la méthode classique réside dans l'étape (3) de relocalisation des inconnues internes. Dans cette nouvelle version, on effectue plusieurs itérations locales en mettant à jour les opérateurs tangents liés à la géométrie, jusqu'à obtenir un résidu inférieur à un seuil déterminé.

$$ \Delta u_i^{(i)} = K_{Tii}^{(i)-1}(u)(f_i^{(i)} - K_{Tib}^{(i)}(u)\Delta u_b^{(i)}) $$

A chaque étape de relocalisation, on force ainsi les sous-structures à vérifier le comportement non-linéaire. Le but est ici de réduire le nombre d'itérations globales de l'algorithme en portant l'effort numérique au niveau local. La méthode peut être vue comme une extension de la technique des « super-éléments », largement répandue pour les calculs linéaires, aux problèmes non-linéaires.

## 3.2 Résultats

### 3.2.1 Convergence et coût de calcul

Les figures 3a) et 3b) présentent les résultats comparatifs pour les différentes stratégies. La première série de résultats (cas 1) est obtenue par une méthode classique de Newton-Krylov-Schur avec un pas de chargement initial maximal (un pas supérieur menant à une divergence de l'algorithme), la deuxième série (cas 2) est obtenue par une approche primale à deux niveaux avec le même pas de chargement initial, alors que la troisième et dernière série (cas 3) est obtenue en déterminant le nouveau pas initial maximal pour cette méthode. Les autres paramètres (critères de convergence,...) sont identiques.

<figure><img src="image_5.png" /><figcaption>FIG. 3 - a) Nombre d'itérations globales du schéma itératif de type Newton. b) Nombre total de résolutions au niveau des sous-structures.</figcaption></figure>

On constate sur ces graphes les gains importants permis par la méthode, à la fois sur le nombre de calculs globaux, et donc d'informations échangées, et sur le nombre de calculs locaux. En particulier, la taille des pas de chargement admissibles est supérieure et l'instabilité de flambage est franchie plus aisément qu'avec la méthode classique. On note de même que les calculs locaux sont plus nombreux pour les sous-structures ayant un comportement fortement non-linéaire que pour les autres. L'effort est réellement porté sur les zones non-linéaires.
### 3.2.2 Sur le critère d'arrêt des calculs locaux

Une étude paramétrique a été réalisée pour déterminer l'influence du critère d'arrêt des calculs locaux sur les performances de la méthode. En particulier, il est apparu qu'il existait un optimum du critère local pour lequel le coût de calcul est minimum. En effet, le nombre de calculs globaux diminue avec le critère puis atteint un seuil minimum, alors que le nombre de calculs locaux diminue tout d'abord avant d'augmenter de nouveau, sans que cela n'ait d'influence positive sur les performances globales. En pratique, un critère de convergence local assez élevé, supérieur de plusieurs ordres de grandeur au critère de convergence global, est suffisant pour obtenir de bons résultats.

## 4 Approche mixte à deux niveaux

### 4.1 Présentation

En s'appuyant sur la méthode LATIN [5], nous explorons les possibilités offertes par la prise en compte de conditions mixtes aux frontières. On choisit ici de conserver comme inconnues aux interfaces $Γ$ à la fois les déplacements et les efforts entre les sous-structures. La solution est recherchée itérativement par la résolution successive de deux groupes d'équations :

- connaissant $\hat{F}$ et $\hat{U}$ aux interfaces, on cherche les champs $f$ et $u$ complets par sous-structure vérifiant le comportement non-linéaire et les conditions aux limites mixtes ($k$ est un paramètre de la méthode) : 
$$ \left\{ \begin{array}{l} f = K(u)u \\ (F - \hat{F}) + k(U - \hat{U}) = 0 \end{array} \right. $$

- connaissant $F = f|_Γ$ et $U = u|_Γ$, on cherche $\hat{F}$ et $\hat{U}$ vérifiant l'admissibilité des champs d'interface (i.e. continuité des déplacements et équilibre des efforts) en prenant comme direction de recherche par sous-structure l'opérateur tangent condensé :

$$ \left\{ \begin{array}{l} \hat{U}_1 = \hat{U}_2 \\ \hat{F}_1 + \hat{F}_2 = 0 \\ S_T^{(s)}(\hat{U}^{(s)} - U^{(s)}) = \hat{F}^{(s)} - F^{(s)} \end{array} \right. $$

ce qui permet de construire le système global suivant, portant sur l'ensemble des inconnues d'interfaces :

$$ \left\{ \begin{array}{l} S_T^{\text{global}} \hat{U} = \sum_s A^{(s)} (S_T^{(s)} U^{(s)} - F^{(s)}) \\ \hat{F}^{(s)} = F^{(s)} + S_T^{(s)} (\hat{U}^{(s)} - U^{(s)}) \end{array} \right. $$

Ce dernier groupe d'équations est résolu itérativement de la même manière que pour BDD ou FETI.

On obtient, à la convergence de l'algorithme, la solution vérifiant toutes les équations du problème. Un des avantages de cette méthode par rapport à la précédente est de relâcher les contraintes imposées sur les problèmes locaux par les conditions aux limites. Elle permet ainsi de réaliser des incréments de charge plus important au niveau global.

### 4.2 Résultats et commentaires

#### 4.2.1 Performances de la méthode mixte

Les essais réalisés sur la structure précédente ont permis de réduire le nombre de calculs globaux en passant de 22 pour le cas 3 de la figure 3.b) (cas optimal de la méthode primale) à 6 résolutions linéaires seulement. Ce résultat a pu être atteint en imposant le chargement global total en un seul pas de chargement, alors que 8 pas de chargement étaient nécessaires dans le cas 1, 6 dans le cas 2 et 5 dans le cas 3. Un pilotage local en effort devient cependant nécessaire pour la résolution des problèmes par sous-structure, ce qui augmente sensiblement leur nombre, sans pour autant dépasser 500 calculs locaux. La méthode s'avère donc intéressante et ouvre la voie à de grands incréments de chargement.
4.2.2 Sur le choix du paramètre k

Le paramètre de recherche *k* définissant les conditions aux limites joue un rôle important dans le résultat de la convergence. On peut remarquer tout d'abord que lorsque *k* tend vers l'infini, la méthode mixte tend vers la méthode primale. A l'inverse, un paramètre nul fait apparaître une version duale de la méthode, où les seules inconnues aux interfaces sont les flux d'efforts. Dans ce cas, il convient de traiter correctement les déplacements de corps rigides qui peuvent apparaître pour les sous-structures sans condition de Dirichlet. Le paramètre *k* apparaît ainsi comme une variable supplémentaire pour optimiser la méthode, les versions primale et duale n'étant que des cas limites. Un choix classique et efficace pour *k* est l'opérateur tangent des sous-structures adjacentes, condensé aux interfaces. Ce dernier est cependant coûteux à calculer et à transférer entre sous-structures. En pratique, on prendra donc des opérateurs plus simples à obtenir, au détriment de la convergence. Les tests effectués avec la diagonale de $K_{Tbb}$, par exemple, donnent des résultats satisfaisants (9 itérations globales) pour un coût de calcul et de transfert très faible.

5 Conclusions et Perspectives

Après avoir rappelé la problématique du calcul non-linéaire de grandes structures et les méthodes
aujourd'hui utilisées, deux nouvelles approches ont été présentées. Elles reposent sur une volonté de
traitement des non-linéarités à une échelle locale. Les résultats obtenus montrent les potentialités des ces
méthodes. Elles permettent de concentrer les efforts de calcul sur les zones effectivement non-linéaires tout
en réduisant le nombre de calculs globaux et donc la quantité d'informations échangées. Elles sont
particulièrement adaptées à des architectures de type cluster de PC, pénalisées par les échanges
d'informations entre processeurs. Une version parallèle pour des structures de coques est en
développement, l'objectif final étant de valider les méthodes sur le calcul d'une structure proche d'un cas
industriel documenté (essai statique d'un tronçon de fuselage Airbus par exemple).

Références

[1] P. Le Tallec, J. Mandel, M. Vidrascu, A Neumann-Neumann Domain Decomposition Algorithm for Solving Plate and Shell Problems, J. Numer. Anal., Vol. 35, No. 2, 836-867, 1998.

[2] C. Farhat, K. Pierson, M. Lesoinne, The second generation FETI methods and their application to the parallel solution of large-scale linear and geometrically non-linear structural analysis problems, Comput. Meth. In Appl. Mech. and Engrg., 18, 333-374, 2000.

[3] Y. Fragakis, M. Papadrakakis, A Unified Framework for Formulating Domain Decomposition Methods in Structural Mechanics, Technical Report, 2002.

[4] X.-C. Cai, D. E. Keyes. Nonlinearly preconditioned inexact Newton algorithms, SIAM J. Sci. Comput., 24, 183-200, 2002.

[5] P. Ladevèze, Nonlinear Computational Structural Mechanics - New Approaches and Non-Incremental Methods of Calculation, Springer Verlag, Berlin, 1999.

[6] B. Haugen, C. Felippa, A unified formulation of co-rotational finite elements (Theory & Applications), Dpt. of Aerospace Engrg. Sciences, Univ. Colorado at Boulder, CAS Report No. CU-CAS-95-06, 1995.

[7] C.Rey, P. Gosselet, Solution to large nonlinear systems : acceleration strategies based on domain decomposition methods and reuse of Krylov-subspaces, 6th International ESAFORM Conference on Material Forming, Salerne (Italie), Nuova Ipsa Editore, 2003.