# Méthode des éléments finis : Développements récents adaptés au calcul de structures complexes

Laurent Champaney, David Dureisseix

► To cite this version:

Laurent Champaney, David Dureisseix. Méthode des éléments finis : Développements récents adaptés au calcul de structures complexes. Université d'été Pratiques industrielles et pratiques de formation en CAO mécanique, Centre National de Ressources en CFIAO, 1997, Cachan, France. hal-01384403

HAL Id: hal-01384403

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

Submitted on 19 Oct 2016

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-SA 4.0 - Attribution - ShareAlike - International
License
# Méthode des éléments finis : Développements récents adaptés au calcul de structures complexes

Laurent CHAMPANEY et David DUREISSEIX

Laboratoire de Mécanique et Technologie
E.N.S. de Cachan / C.N.R.S. / Université Paris 6
61, avenue du Président Wilson - 94235 CACHAN CEDEX
e-mail : [Laurent.Champaney; David.Dureisseix]@lmt.ens-cachan.fr
http://www.lmt.ens-cachan.fr

Ce document est le preprint de la publication : L. Champaney, D. Dureisseix, Méthodes des éléments finis : développements récents adaptés aux calculs de structures complexes, Actes de l'université d'été Pratiques industrielles et pratiques de formation en CAO mécanique, Centre National de Ressources en CFIAO, Cachan, septembre 1997.

## 1. Introduction

Le développement et la commercialisation de nombreux codes de calcul par la méthode des éléments finis rendent la simulation de systèmes mécaniques accessible à tous, même aux plus novices d'entre nous. Les systèmes de modélisation et de CAO permettent aisément d'obtenir un modèle éléments finis, même pour des géométries complexes. De plus, les systèmes de post-traitement graphiques tracent des cartes d'isovaleurs (de contraintes, de déformations, ...) avec une qualité et un nombre de couleurs impressionnants. Néanmoins il convient de rester lucide et d'être conscients de la qualité des calculs que nous menons. Même si elle se banalise, la méthode des éléments finis reste un outil complexe et savoir bien la maîtriser n'est pas chose aisée. Nous présentons ici quelques notions sur le contrôle des calculs et sur certains outils en cours de développement pour donner à l'utilisateur l'assurance d'une qualité de son modèle. Un exemple de maîtrise de modèle élastique 3D est détaillé. De plus, nous essayons de donner quelques notions sur les difficultés rencontrées lors de l'étude de structures complexes qui n'est pas encore, de nos jours, à la portée de tous. Par "complexe", nous entendons complexité du système mécanique (géométrie, liaisons, ...), complexité des sollicitations appliquées (chargements divers, dynamiques, cycliques, ..., grands déplacements, ...) et complexité du modèle de matériau utilisé (plasticité, viscoplasticité, endommagement, ...).

Nous présentons les problèmes numériques rencontrés lors de l'étude, par la méthode des éléments finis, de structures complexes et nous donnons quelques pistes pour faciliter le traitement de certains types de problèmes. Un exemple de méthode adaptée au calcul d'assemblages complexes de structures élastiques est proposé. Après avoir donné quelques idées sur les problèmes de taille et de coût de calcul, nous présentons les perspectives offertes par les calculateurs à architecture parallèle modernes. Un exemple illustre ces possibilités.

## 2. Méthode des éléments finis

### 2.1 But

Le dimensionnement d'un système mécanique (pièce, assemblage,...) nécessite la réalisation de plusieurs modélisations :

- modélisation de la géométrie,

- modélisation des conditions aux limites et des charges,

- modélisation du comportement du ou des matériaux,
de manière à obtenir une formulation mécanique du problème traité. Cela conduit au choix d'un modèle particulier parmi ceux que nous propose la mécanique (tridimensionnel, plaque, poutre,..., élastique, plastique,...). Généralement, les modèles ainsi construits relèvent de la mécanique des milieux continus et se traduisent par des équations aux dérivées partielles, avec des conditions aux limites entre les inconnues (champ de déplacement U, champ de contrainte $\sigma$,...) et les données (conditions initiales, charges, liaisons, coefficients matériaux,...).

La solution exacte ($U_{ex}$, $\sigma_{ex}$) d'un tel problème est rarement possible à déterminer analytiquement et la méthode la plus couramment utilisée pour trouver une solution approchée ($U_h$, $\sigma_h$) est la méthode des éléments finis [1]. Elle consiste à remplacer le modèle continu par un modèle discret.

## 2.2 Principes

Une discrétisation spatiale de la géométrie est réalisée à l'aide d'éléments simples (segments, triangles, quadrangles, cubes, prismes, tétraèdres, ...). C'est le maillage. Sur chaque élément, une représentation du déplacement $U = \{N\}^T\{u\}$ est introduite. $\{N\}$ contient les fonctions d'interpolation (polynômes de degré peu élevé) sur les éléments utilisés et $\{u\}$ est le vecteur de déplacements nodaux recherché. Le système linéaire à résoudre, issu de la formulation en déplacement du problème, est alors du type :

$$[K]\{u\} = \{f\}$$

où $\{f\}$ est le vecteur de forces nodales appliquées (ou forces généralisées) et $[K]$ la matrice de rigidité du système — qui est un assemblage de matrices de rigidité élémentaires sur les éléments de la discrétisation. Le système linéaire a la taille du nombre d'inconnues nodales (degrés de liberté). Il est creux et souvent mal conditionné. La résolution se fait par factorisation de la matrice, par une méthode de type Cholesky, par exemple.

La qualité de la solution dépend du type des éléments, de leur qualité géométrique et du raffinement du maillage (taille des éléments). Le choix du type d'élément adapté au problème considéré n'est pas simple. Les modèles 2D ou 3D ne sont pas forcément meilleurs que les modèles 1D. Ils peuvent dans certains cas conduire à des résultats aberrants.

De plus, la qualité de la discrétisation réalisée se doit d'être maîtrisée sous peine d'obtenir des résultats très éloignés de ceux attendus. L'estimation de la qualité de la discrétisation fait l'objet du paragraphe 4.

# 3. Problématique du calcul de structures complexes

## 3.1 Non linéarités

Il existe trois types de non linéarités dans les problèmes mécaniques que nous traitons :

- les non linéarités de type "matériau" qui sont associées au caractère non linéaire de la loi de comportement des matériaux en présence ;

- les non linéarités de type "géométrique" qui interviennent lorsque les déplacements et/ou les déformations de la structure deviennent tels qu'il n'est plus possible de traiter le problème dans sa configuration non déformée (l'hypothèse de "petites perturbations" n'est plus valable) ;

- les non linéarités de type "contact" qui apparaissent lorsqu'on veut modéliser des conditions de non pénétration et/ou de frottement entre deux solides. La solution dépend de l'état de la zone de contact (décollement, contact, glissement, adhérence,...).

La résolution de ces problèmes entraîne donc l'utilisation d'algorithmes itératifs pour la résolution de systèmes non-linéaires. De plus, lorsque le comportement de la structure varie avec l'intensité du
chargement et avec son histoire, il est nécessaire d'introduire une discrétisation temporelle. Des méthodes incrémentales, basées sur le fait que la solution varie peu entre deux instants, permettent de calculer la solution au cours du temps.

## 3.2 Taille des modèles

Dans le cas de modélisations tridimensionnelles, la complexité de la géométrie et la recherche d'une solution de bonne qualité conduisent à utiliser une discrétisation fine et donc à introduire un grand nombre d'éléments. Les systèmes linéaires engendrés deviennent très rapidement gigantesques (plusieurs centaines de milliers d'équations).

<figure><img src="image_2.png" /><figcaption>Figure 1 : Comparaison des tailles des problèmes 2D et 3D.</figcaption></figure>

À titre d'exemple, on montre sur la figure 1, l'évolution du nombre de degrés de liberté quand le nombre d'éléments augmente, pour un maillage 2D en quadrangles à 8 nœuds et pour un maillage 3D en cubes à 20 nœuds.

## 3.3 Problèmes numériques

Les problèmes posés sont alors multiples. La mise en place du modèle géométrique fin et du modèle de calcul contenant autant d'inconnues amène à manipuler de très grosses quantités d'informations. De plus, la matrice de rigidité, une fois assemblée, occupe un espace mémoire très important, souvent plus grand que celui disponible sur les calculateurs courants. Son stockage et sa factorisation deviennent délicats. Des écritures sur disques sont donc nécessaires et ralentissent fortement le processus. Il est alors nécessaire, pour de telles modélisations, de générer des outils facilitant l'entrée et la manipulation des données et conduisant à des encombrements mémoire réduits.

Pour les problèmes non-linéaires, d'autres difficultés apparaissent. La volonté de prendre en compte les effets de structure (zones à fort gradient de contrainte, redistribution de contraintes,...), aussi bien que celle de modéliser finement le comportement des matériaux (plasticité, viscoplasticité, matériaux composites,...) en développant des lois de comportement sophistiquées prenant en compte de nombreux couplages mécaniques, de nombreuses variables internes, ou même le développement de techniques de calcul adaptées aux problèmes très fortement non-linéaires, conduisent à des coûts de calcul prohibitifs, voire hors de portée, pour les architectures des ordinateurs séquentiels actuels. De plus, lorsque le chargement est complexe (chargement cyclique, dynamique rapide, ...), la discrétisation temporelle se doit d'être fine. Un grand nombre d'instants est donc à considérer. Le nombre d'opérations à réaliser est donc très important et le volume de résultats à stocker est considérable.
# 4. Contrôle de la qualité des modèles éléments finis

## 4.1 Qualité d'un modèle

Le contrôle de la qualité du calcul élément finis consiste à considérer le modèle mécanique continu comme modèle de référence et à évaluer la qualité de la solution calculée ($U_h, \sigma_h$) comme approximation de la solution exacte ($U_{ex}, \sigma_{ex}$). On évalue ainsi les erreurs de modélisation introduites pendant la phase de discrétisation.

Deux situations existent suivant que l'estimation des erreurs se déroule avant le calcul EF (erreur dite a priori) ou après (erreur dite a posteriori). La première approche ne donne, en pratique, aucune information sur le niveau réel de l'erreur de discrétisation. Elle met seulement en évidence la convergence éventuelle de la solution approchée vers la solution exacte lorsque le pas de la discrétisation diminue. Les estimations d'erreur a posteriori permettent, elles, d'obtenir le niveau de l'erreur globale ainsi que des contributions locales par élément. Lors des vingt dernières années, de nombreuses études ont été réalisées sur les estimations d'erreurs a posteriori. On peut distinguer parmi ces travaux :

- la démarche de LADEVEZE, qui utilise la notion d'erreur en relation de comportement et des techniques adaptés de reconstruction de champs admissibles (il est montré que l'erreur calculée par cette approche majore l'erreur réelle) ;

- les travaux de BABUSKA, qui utilisent la notion de résidu mesurant la non vérification des équations d'équilibre ;

- la démarche la plus récente de ZIENKIEWICZ et ZHU, qui utilise l'écart en énergie entre la contrainte élément finis calculée et une contrainte obtenue par lissage.

## 4.2 Adaptativité

Le maillage (qualité et taille) joue sur la qualité de la solution éléments finis. Pour obtenir une qualité cible prescrite par l'utilisateur, il faut pouvoir adapter le maillage. Le couplage d'une méthode de remaillage avec les techniques de calcul d'erreur de discrétisation et les mailleurs automatiques permet l'adaptation des maillages. La méthode de remaillage la plus couramment utilisée est basée sur une reconstruction complète de maillage à partir d'une carte de modification de taille fournie par le calcul d'erreur par élément.

Une procédure d'adaptation de maillage est basée sur le schéma de principe défini sur la figure 2.

<figure><img src="image_3.png" /><figcaption>Figure 2 : Procédure d'adaptation de maillage</figcaption></figure>
Le but des recherches en cours sur l'adaptation des maillages est la réalisation de procédures complètement automatiques et transparentes pour l'utilisateur, ne nécessitant que la donnée du modèle géométrique et de l'erreur cible. Le frein à cette automatisation est actuellement la qualité des mailleurs respectant une carte de tailles imposée (surtout dans le cas des modélisations 3D). De plus, de telles procédures sont assez peu développées pour les problèmes non-linéaires et pour les problèmes de dynamique.

## 4.3 Exemple : chambre de dilatométrie haute pression

Nous présentons ici un exemple de procédure d'adaptation de maillage 3D pour les problèmes d'élasticité linéaire, réalisé par J.P. DUMEAU au LMT-CACHAN [2]. Le calcul d'erreur est basé sur le concept d'erreur en relation de comportement. Le modèle géométrique est réalisé avec IDEAS, le maillage avec le mailleur prototype GHS3D (INRIA) [3] et les calculs avec le code industriel CASTEM 2000 (CEA) [4]. Il s'agit d'un quart de corps de chambre de dilatométrie haute pression (jusqu'à 1 GPa) munie, de façon non symétrique, de cinq fenêtres optiques (Fig. 3). Ces dernières doivent être les plus rapprochées possible du centre pour accéder à une plus grande ouverture numérique, de manière à disposer les appareillages de mesure optique. Les dimensions de la zone utile correspondent à un diamètre de l'ordre de 45mm sur une profondeur de 80mm. Le diamètre extérieur de l'enceinte peut varier de 210 à 220mm.

<figure><img src="image_4.png" /><figcaption>Figure 3 : Géométrie du corps de l'enceinte haute pression</figcaption></figure>

Le chargement prend en compte la pression interne (1 GPa), la réaction des fenêtres optiques et des joints d'étanchéité. Le processus automatique d'adaptation de maillage est initialisé sur le maillage régulier présenté sur la figure 4. La contrainte de Mises correspondante est tracée figure 5.

<figure><img src="image_5.png" /><figcaption>Figure 4 : Maillage initial : 4747 tétraèdres à 10 nœuds (7635 nœuds)

Erreur calculée : 73.94%</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5 : Maillage initial : isocontrainte de Mises</figcaption></figure>

Avec ce maillage initial, l'erreur globale calculée approche les 80%. L'erreur cible souhaitée par l'utilisateur est de 20% et est atteinte après deux itérations du processus d'adaptation. Les figures 6 et 8 présentent les maillages obtenus au cours des itérations. Les figures 7 et 9 donnent les isocontraintes de Mises correspondantes.

<figure><img src="image_7.png" /><figcaption>Figure 6 : Maillage itération 1 : 18316 tétraèdres à 10 nœuds (29653 nœuds)

Érreur calculée : 30.84%</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7 : Maillage itération 1 : isocontrainte de Mises</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8 : Maillage itération 2 : 21160 tétraèdres à 10 nœuds (34104 nœuds)

Erreur calculée : 23.62%</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9 : Maillage itération 2 : isocontrainte de Mises</figcaption></figure>

On peut remarquer que la répartition des contraintes de Mises varie fortement au cours des itérations. La répartition obtenue avec le maillage initial (déjà assez fin) est très éloignée de la répartition finale. De plus, le nombre d'éléments nécessaire pour l'obtention d'une erreur globale proche de 20% est très grand.

# 5. Méthodes adaptées au calcul et à l'optimisation de structures

## 5.1 Outils spécialisés

Un grand nombre de problèmes interviennent dans le calcul de structures complexes à comportement non linéaire. Il apparaît alors que les approches généralistes (telles que celles proposées par les grands codes de calcul industriels) offrent des performances médiocres pour ces calculs complexes. Il faut donc avoir recours à des outils spécialisés pour chaque type de problème. Nous présentons ici l'approche COFAST (pour CContact et Frottement dans les Assemblages de SStructures), développée au LMT-Cachan depuis 1990, qui est spécialisée dans le calcul d'assemblages de structures élastiques.

## 5.2 Approche modulaire COFAST

### 5.2.1 Calcul d'assemblages de structures

Les problèmes rencontrés lors de la modélisation et du calcul d'assemblages de structures complexes sont multiples :

- lorsque l'assemblage contient beaucoup de composants, la description de tous les composants et de toutes les liaisons devient un problème complexe présentant de nombreuses données et paramètres. De plus, les matrices de rigidités de tels problèmes sont souvent de taille très importante ;
- les assemblages font très souvent intervenir des liaisons par contact qui rendent le problème non-linéaire. Pour une résolution par la méthode des éléments finis et lorsque ces liaisons sont nombreuses, le nombre d'équations de contact devient très important. Le processus itératif de résolution des équations de contact nécessite de fréquentes modifications de la matrice de rigidité. Dans ce cas les performances des codes de calcul généralistes deviennent très médiocres ;

- dans le cadre d'une boucle d'optimisation, la modification d'un composant ou des caractéristiques d'une liaison entraîne la modification de la matrice de rigidité de l'assemblage et donc nécessite le réassemblage et la re-factorisation de matrice de grande taille.

### 5.2.2 Sous-structuration par entités technologiques

L'approche proposée utilise une décomposition de l'assemblage étudié en deux groupes d'entités technologiques différentes : les sous-structures (les composants de l'assemblage) et les interfaces (les liaisons entre les composants). Les interfaces réalisent les liaisons entre les sous-structures mais aussi entre les sous-structures et l'extérieur pour la prise en compte des conditions aux limites. L'assemblage complet considéré est représenté par un ensemble de sous-structures qui ne communiquent qu'avec leurs interfaces voisines. Une telle décomposition d'un assemblage est réalisée en deux étapes. Une première séparation des éléments de l'assemblage fait intervenir des interfaces modélisant les liaisons réelles entre ces différents éléments (contact, frottement, jeu,...). Ensuite, chaque élément peut être décomposé en sous-structures et les interfaces engendrées modélisent alors des liaisons parfaites (Fig. 10).

<figure><img src="image_11.png" /><figcaption>Figure 10 : Décomposition d'un assemblage.</figcaption></figure>

Dans cette approche et contrairement aux techniques classiques de décomposition de domaine, les interfaces ne sont pas des objets fictifs engendrés par la décomposition. Ce sont des entités mécaniques à part entière qui disposent de leurs propres inconnues et de leurs propres équations d'équilibre et de comportement. Elles sont les éléments clés de l'approche, car le fait qu'à une même interface puissent être associés divers comportements différents introduit une grande flexibilité des modélisations réalisées.

Un champ de déplacement et un champ d'effort représentent l'interaction entre une sous-structure et une interface. Ces deux champs sont définis sur les interfaces et sur les frontières des sous-structures. Les conditions aux limites sur la structure complète sont imposées par l'intermédiaire d'interfaces adaptées entre les sous-structures et l'extérieur.

### 5.2.3 Résolution

Dans les problèmes traités, les composants de l'assemblage sont supposés avoir un comportement élastique linéaire. Les seules non-linéarités proviennent des liaisons (contact, frottement, jeux, joint élastomères,...). Le schéma itératif utilisé est basé sur la méthode non-incrémentale appelée LATIN méthode. Le premier principe est de séparer les difficultés du problème non-linéaire étudié [5]. Les
équations sont séparées en deux groupes :

- équations linéaires globales sur une sous-structure : liaison cinématique avec le champ de déplacement interface, équilibre sous l'action du champ d'effort interface, loi de comportement ;

- équations non-linéaires locales associées à une liaison : relation de comportement liant efforts et déplacements de part et d'autre de l'interface. Ces équations dépendent du type de comportement à modéliser par l'interface.

L'algorithme utilisé construit des solutions vérifiant successivement l'équilibre des sous-structures puis l'équilibre des liaisons. Chaque itération donne une solution approchée du problème dont la qualité est connue. La détermination d'une solution vérifiant l'équilibre des sous-structures conduit à la résolution d'un problème linéaire de petite taille par sous-structure (qui est traité par une approche éléments finis). La détermination d'une solution vérifiant l'équilibre des liaisons conduit à la résolution d'un problème non-linéaire local de petite taille en chaque point de l'interface.

### 5.2.4 Avantages

Dans l'approche COFAST, la décomposition en sous-structures permet de réduire la taille de la modélisation (encombrement des matrices de rigidité) par rapport à une approche éléments finis classique basée sur la même discrétisation. L'algorithme permet, de plus, une réduction des coûts de calcul engendrés, même sur ordinateur séquentiel. En effet, les matrices de rigidité utilisées sont indépendantes des conditions aux limites et des statuts de contact ; elles restent donc constantes au cours des itérations et ne sont factorisées qu'une seule fois. Il est alors possible de stocker, dans un fichier, les données relatives à un composant (géométrie, maillage, matrice de rigidité factorisée,...). Ainsi, ce composant peut être réutilisé directement dans plusieurs calculs, assemblé à des composants différents, sans que sa matrice de rigidité ne soit refactorisée. De plus, pour limiter l'encombrement des problèmes présentant une répétition de la géométrie (sans répétition du chargement), il est possible d'utiliser une seule matrice de rigidité pour représenter une sous-structure et ses équivalents répétés. La résolution se fait alors avec plusieurs seconds membres d'efforts (schématiquement) sur les différentes sous-structures représentées par la matrice. Des gains très importants sont obtenus pour ces situations. Un exemple de déversement d'un roulement à billes est présenté dans [6].

L'approche étant basée sur une formulation complément parallèle (problèmes indépendants par sous-structure et problèmes indépendants par interface), son utilisation sur calculateur à architecture parallèle se fait sans aucune difficulté (cf § 6). Enfin, la représentation des champs d'efforts et de déplacements qui est utilisée, permet de s'affranchir de la compatibilité des maillages des deux sous-structures de part et d'autre d'une interface. Ceci permet de simplifier considérablement la phase de maillage de la structure [5] et d'utiliser un même composant en liaison avec différents autres éléments, sans avoir à modifier son maillage.

## 5.3 Exemple : liaison complète pignon-arbre creux

Cet exemple a été réalisé dans le cadre d'un bureau d'étude en préparation à l'agrégation de Mécanique à l'ENS de Cachan. Le thème de ce bureau d'étude est le dimensionnement d'un réducteur à arbre creux. La présentation faite aux étudiants avait pour but de montrer les apports du calcul dans le dimensionnement de la liaison encastrement entre le pignon et l'arbre creux (Fig. 11) par rapport aux méthodes analytiques classiques en bureau d'étude. L'étude porte sur le nombre de boulons nécessaires à la transmission du couple, leur précharge et la répartition de pression en charge sur la zone de contact entre le pignon et l'épaulement de l'arbre.
<figure><img src="image_12.png" /><figcaption>Figure 11 : Pignon vissé sur un arbre creux.</figcaption></figure>

La figure 11 présente la géométrie de la liaison et son éclaté (la denture du pignon n'est pas représentée). Le modèle est complètement paramétré. L'arbre est supposé encastré en bout et un effort — paramétré par le couple à transmettre et par la géométrie de la denture — est appliqué sur la partie supérieure du pignon. La figure 12 présente la décomposition en 54 sous-structures et 109 interfaces pour un modèle à neuf boulons de fixation. Ce modèle à neuf boulons comporte 141 000 ddl.

<figure><img src="image_13.png" /><figcaption>Figure 12 : Décomposition en sous-structures et interfaces.</figcaption></figure>

La figure 13 présente, pour un couple appliqué de 1 000 Nm, la répartition des contraintes axiales dans une portion de l'assemblage et dans les boulons associés. Elle présente ainsi, la répartition de pression sur la zone de contact qui est très éloignée de l'hypothèse de répartition constante effectuée pour le calcul analytique. On remarque alors que le calcul analytique classique est fortement surdimensionnant.
<figure><img src="image_14.png" /><figcaption>Figure 13 : Contraintes axiale sur une portion de l'assemblage</figcaption></figure>

# 6. Calcul parallèle

## 6.1 Possibilité des calculateurs parallèles actuels

Les machines mono-processeur ont vu la puissance de base de leur processeur augmenter notablement ; cependant, des limites physiques sont apparues et pour permettre le passage à une classe de machines supérieure, les architectures vectorielles puis parallèles ont émergé. Aujourd'hui, ces dernières sont celles qui permettent d'atteindre les plus grandes puissances de calcul (figure 14), en particulier avec le développement de processeurs et de réseaux d'interconnection permettant d'obtenir un parallélisme "extensible" utilisant pleinement les possibilités d'une machine (choix du nombre de processeurs suivant la taille du problème à traiter, avec pas ou peu de dégradation de performances quand ces deux s'accroissent).

<figure><img src="image_15.png" /><figcaption>Figure 14 : Evolution de la puissance de calcul des ordinateurs.</figcaption></figure>

## 6.2 Parallélisme informatique

Pour tirer parti de ces nouvelles architectures, des algorithmes adaptés ont dû être développés, la parallélisation automatique (que l'on pourrait qualifier de parallélisme "informatique") de codes de calculs existant sur machines séquentielles ne permettant que difficilement d'atteindre des performances élevées et étant délicate à réaliser sur des architectures à mémoire distribuée.
Outre ces techniques, se développent aussi des méthodes de résolution en parallèle de grands systèmes linéaires creux (typiques de ceux auxquels conduit une discrétisation éléments finis), que ce soit de manière directe ou itérative, de façon quasi indépendante du type de problème qui leur a donné lieu. Il s'agirait alors de parallélisme "numérique".

## 6.3 Parallélisme mécanique

Ici, nous nous plaçons dans une troisième catégorie : celle d'un parallélisme "mécanique", dans laquelle se classent les méthodes de type décomposition de domaine, qui conservent au cours de la résolution la nature du problème mécanique à résoudre. Pour rompre la globalité en variables d'espace des équations d'équilibre d'une structure, ces méthodes procèdent à une décomposition en sous-domaines pour ne plus faire intervenir qu'une succession de résolution de problèmes quasi-locaux sous-domaine par sous-domaine. Le problème résultant est alors, bien entendu, celui du raccord entre sous-domaines voisins.

Ces approches conduisent à un parallélisme à "grain large" (le volume de calcul à réaliser par processeur reste élevé vis-à-vis de la quantité d'informations à échanger avec les autres processeurs au travers du réseau de communication). Généralement, la nature synchrone de ces algorithmes pousse à se préoccuper de l'équilibrage des charges entre processeurs pour obtenir des efficacités élevées.

## 6.4 Exemple : différentiel automobile

L'approche COFAST, décrite au paragraphe 5.2, rentre dans le cadre de ces méthodes de décomposition de domaine [7].

Sur un ordinateur multi-processeur, les sous-structures et interfaces peuvent être distribuées sur les différents nœuds de calcul. L'échange d'informations entre processeurs est réalisé par envoi de messages. Pour un problème donné, l'augmentation du nombre de sous-structures permet de réduire l'encombrement associé aux matrices de rigidités. Alors, chaque programme, associé aux différents processeurs, doit gérer plusieurs sous-structures. Cette technique permet de bien répartir la charge sur les processeurs mais conduit à un algorithme parallèle plus complexe. Plusieurs tests ont été effectués sur le CRAY T3D de l'IDRIS (Orsay) et sur l'IBM SP2 du C.E.A. (Saclay). Ce dernier dispose de 16 processeurs avec chacun 128 Mo de RAM. Le système PVM (Parallel Virtual Machine System) est utilisé pour gérer l'envoi de messages entre les différents processeurs.

L'exemple proposé est celui d'un différentiel automobile, pour lequel on modélise uniquement le boîtier, la couronne, le chapeau et les vis de fixation (figure 15). On considère des actions sur une dent de la couronne et les actions des roulements à bille sur le boîtier sont modélisées par l'intermédiaire de liaisons de type "contact unilatéral sans frottement". L'axe des satellites est supposé encastré. Les actions des satellites et des planétaires sur le boîtier sont modélisées par des efforts dont l'amplitude et la direction sont paramétrées par le couple appliqué et par la géométrie des dentures (figure 16). Les précontraintes des vis sont imposées par l'intermédiaire d'interfaces adaptées, situées entre le corps et la tête de vis. Les sollicitations mécaniques associées à la rotation d'ensemble sont aussi prises en compte. Des interfaces de type contact unilatéral avec frottement sont utilisées pour modéliser les liaisons entre les différents composants (boîtier, couronne, chapeau et vis).
<figure><img src="image_16.png" /><figcaption>Figure 15 : Différentiel et son modèle.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 16 : Modèle utilisé et chargement (3/4 du modèle).</figcaption></figure>

Six modélisations différentes ont été réalisées (3 maillages différents et 2 décompositions par maillage). Le premier découpage en 39 sous-structures est la décomposition minimale de l'assemblage : le boîtier, le chapeau, la couronne, 18 corps de vis et 18 têtes. Les symétries du modèle et un décomposeur automatique ont été utilisés pour effectuer le second découpage en 84 sous-structures et 223 interfaces (figure 17). Le tableau 1, présente les caractéristiques des différents maillages adoptés et l'encombrement de l'ensemble des matrices de rigidité pour les différents découpages. Ces résultats montrent que l'augmentation du nombre de sous-structures permet de réduire considérablement l'encombrement du problème.

<figure><img src="image_18.png" /><figcaption>Figure 17 : Découpage en 84 sous-structures (1/4 du modèle).</figcaption></figure>
<table><caption>Tableau 1 : Caractéristiques des différents problèmes.</caption><thead><tr><th></th><th>nbre de ddl</th><th>taille pour 39 s-s</th><th>taille pour 84 s-s</th></tr></thead><tbody><tr><td>maillage 1</td><td>35 136</td><td>128 Mo</td><td>31 Mo</td></tr><tr><td>maillage 2</td><td>129 924</td><td>1391 Mo</td><td>310 Mo</td></tr><tr><td>maillage 3</td><td>236 976</td><td>4 574 Mo</td><td>790 Mo</td></tr></tbody></table>

Pour un couple appliqué de 2 000 Nm et une vitesse de rotation de 2000 tr/min, la convergence est atteinte après 200 itérations. La figure 18 donne les contraintes équivalentes de Mises à l'itération 200 tracées sur le boîtier et sur le chapeau. On donne aussi les contraintes axiales ($\sigma_{xx}$) sur les vis et sur la couronne. Cette figure donne la répartition de pression de contact entre la couronne et le boîtier (liaison avec frottement).

<figure><img src="image_19.png" /><figcaption>Figure 18 : Résultats à l'itération 200.</figcaption></figure>

La figure 19 montre, sur un quart du différentiel et pour la deuxième décomposition, de quelle façon les sous-structures ont été groupées pour un fonctionnement avec 8 et 12 processeurs.

<figure><img src="image_20.png" /><figcaption>Figure 19 : Regroupement des sous-structures pour 8 et 12 processeurs (1/4 du modèle).</figcaption></figure>

Le temps de calcul CPU cumulé en fonction du nombre de processeurs pour le modèle à 35 136 ddl (figure 20) montre que le coût des communications est relativement faible (pour un processeur, il n'y a pas de communication). La figure 20 présente aussi le temps CPU cumulé en fonction de la taille de la modélisation pour l'utilisation de 12 processeurs. Il est important de noter que le coût total augmente lentement avec la taille du problème.
<figure><img src="image_21.png" /><figcaption>Figure 19 : Temps CPU cumulé sur l'IBM SP2.</figcaption></figure>

# 7. Conclusions et Perspectives

Nous avons présenté les problèmes numériques intervenant dans les calculs éléments finis de structures complexes. Ces problèmes s'expriment en terme de difficulté de gestion des données, de taille et de coût de calcul et en terme de maîtrise du modèle. Nous avons montré quelques outils pour la maîtrise et le contrôle de la qualité des modèles éléments finis. Ces procédures commencent à être transparentes pour l'utilisateur. Néanmoins, elles ne sont pas encore bien développées pour le contrôle de modèles non-linéaires ni pour les problèmes de dynamique. De plus, les processus présentés nécessitent la donnée des conditions aux limites sur la structure. Lorsque celles-ci ne sont a priori pas connues, comme dans le cas du contact unilatéral où la zone de contact évolue avec le chargement, le problème reste très délicat. Pour les problèmes complexes, l'intérêt des approches spécialistes a été montré. Il reste cependant un grand nombre de problèmes complexes pour lesquels les approches généralistes sont encore les seules disponibles et où les coûts de calcul sont faramineux.

Le développement rapide des possibilités des calculateurs à architecture parallèle offre de nous nouvelles perspectives. Néanmoins, les approches efficaces sur ce type d'ordinateur sont encore peu nombreuses et restent limitées à des problèmes faiblement non-linéaires. Il est certain que les années à venir amènerons un grand développement de telles techniques et leur amélioration par des techniques adaptées (approches multi-échelles, ...).

# Références

[1] Gupta K.K., Meek J.L. — A brief history of the beginning of the finite element method, *International Journal of Numerical Methods in Engineering*, vol 49, 3761-3774 (1996).

[2] Dumeau J.P. — *Contôle et adaptation des maillages 3D : Application à l'automatisation des calculs*, Thèse de doctorat, ENS de Cachan (1996).

[3] George P.L. — *Génération automatique de maillages : application aux méthodes d'éléments finis*, Collection Recherches en Mathématiques Appliquées, Masson, Paris (1991).

[4] Verpeaux, P., Charras, T. and Millard, A. — CASTEM 2000 une approche moderne du calcul des structures, *Calcul des structures et intelligence artificielle* (Fouet J. M., Ladevèze P. et Ohayon R., Ed.), Ed. Pluralis, 261-271 (1988).

[5] Champaney L. — *Une nouvelle approche modulaire pour l'analyse d'assemblages de structures tridimensionnelles*, Thèse de doctorat, ENS de Cachan (1996).

[6] Dureisseix D. — *Une approche multi-échelles pour des calculs de structures sur ordinateurs à architecture parallèle*, Thèse de doctorat, ENS de Cachan (1997).

[7] Champaney L., Cognard J.Y., Dureisseix D. et Ladevèze P. — Large scale applications on parallel computers of a mixed domain decomposition method, *Computational Mechanics*, Vol. 19, n°4, 253-262, (1997).