A1. Formulation générale pour les calculs élastiques linéaires

A1. Formulation générale pour les calculs élastiques linéaires

La formulation qui n’est plus réellement à faire car la majorité des logiciels proposés de nos jours s’appuie sur une formulation en déplacement car la prise en compte des conditions aux limites ne pose pas de difficulté notoire à l’inverse d’une formulation en contrainte. Il existe également la formulation mixte en déplacement et en contrainte. Ces trois formulations sont les principales, mais il existe d’autres formulations que nous n’évoquerons pas dans ce chapitre.

Dans ce document, nous ne présentons que la formulation en déplacement (la plus utilisée) : dans ce cadre, les inconnues de base sont les déplacements aux nœuds, notés classiquement q, dont sont déduits tous les autres résultats (déformations ε, contraintes σ …).

Les déplacements en tout point d'un élément, notés u(x,y,z), sont obtenus par interpolation des valeurs aux nœuds, à l'aide des fonctions d'interpolation (appelées fonction de forme), notées N(x,y,z), polynomiales de degré faible (1 à 3 en général) :

La relation déformations - déplacements s'obtient par dérivation à l'aide de l'opérateur (matrice) de dérivation D (différent selon le contexte), faisant apparaître la matrice B reliant déformations et déplacements aux nœuds :

Pour obtenir les contraintes, il est enfin indispensable d'introduire la loi de comportement du matériau constitutif. Ainsi, en élasticité linéaire, la relation contrainte – déformation s'écrit (avec H la matrice de Hooke) :

On constate donc que toutes les quantités d'intérêt (déplacements en tout point, déformations, contraintes) s'obtiennent à partir de la connaissance des déplacements aux nœuds du maillage

Remarque : comme les déformations (et donc les contraintes) sont obtenues par dérivation des déplacements aux nœuds, il y a dégradation de la précision lors du calcul de ces grandeurs.

Pour un calcul statique en mécanique des structures, la détermination de la solution numérique du problème par la Méthode des Éléments Finis (MEF) peut se résumer en 4 étapes principales :

Étape 1. Détermination des matrices et des vecteurs élémentaires sur chaque élément de volume V par les relations suivantes (la notation B^T désignant la transposée de la matrice B) :

 matrice élémentaire de rigidité

 vecteur élémentaire des charges équivalentes

avec  le vecteur des charges volumiques et  celui des charges surfaciques. Le vecteur des charges équivalentes permet de « ramener aux nœuds » les charges appliquées aux éléments (dans le volume ou en surface), selon l'interpolation retenue

Attention, cela ne correspond en général pas à l'équirépartition aux nœuds de la charge totale appliquée à un élément.

Étape 2. Calcul de la matrice de rigidité K et du vecteur des charges F à l’échelle de la structure par assemblage des quantités élémentaires et introduction des conditions aux limites.

Étape 3. Détermination du vecteur q des déplacements nodaux par la résolution du système linéaire (de grande taille) :

Le vecteur q donne la valeur du déplacement aux nœuds du maillage, comme c'est le cas du vecteur des charges équivalentes F.

Etape 4. Déduction des quantités d’intérêt telles que déformations et contraintes par un post-traitement du vecteur q des déplacements nodaux.

En revenant au niveau des éléments, les valeurs des déplacements en tout point, des déformations et des contraintes sont déduites des déplacements aux nœuds à l'aide des relations précédentes. Contrairement aux déplacements, les contraintes et les déformations sont des grandeurs calculées élément par élément (et il n'y a pas toujours continuité de ces quantités d'un élément à l'autre).

Remarque : il est possible de faire une analogie entre résolution par éléments finis d'un problème mécanique et d'un problème thermique en régime permanent linéaire : dans ce cas, les températures aux nœuds T sont obtenues en résolvant le système linéaire : Λ.T = Φ, Λ étant la matrice de conductivité et Φ le vecteur des flux aux nœuds.

Le degré des fonctions de forme N(x,y,z), qui se déduit pour les éléments massifs (ou de mécanique des milieux continus - MMC) du nombre de nœuds sur l'arête d'un élément fini (2 nœuds : éléments linéaires, 3 nœuds : éléments quadratiques) se répercute sur le degré des déformations et des contraintes évaluées par éléments finis. Par exemple, si l'opérateur de dérivation D correspond à des dérivées premières (cas de la MMC), les déformations et contraintes seront constantes par éléments pour des éléments linéaires, linéaires pour des éléments quadratiques… Pour mémoire, en RdM, les déformations font intervenir des dérivées d'ordre 2 (flexion des poutres, plaques…).

Fonctions de forme - Triangle à 3 noeuds (fonctions linéaires)

Fonctions de forme - Triangle à 6 noeuds (fonctions quadratiques)

Attention : il n'est pas possible de mélanger dans un même maillage des éléments finis dont les fonctions de forme sont de degrés différents (par exemple éléments linéaires et quadratiques). En effet, la continuité des déplacements ne peut être assurée à l'interface d'éléments de fonctions de forme de degrés différents.

Dans les logiciels, la plupart des éléments finis sont dits isoparamétriques, c'est-à-dire que les différents éléments finis d'un maillage sont construits par transformation géométrique d'un élément de référence (lui-même exprimé dans un repère de coordonnées unitaires ξ,η,χ).

Élément isoparamétrique (quadratique) : différentes transformations (2D)

Ces transformations géométriques, reposant en général sur les mêmes nœuds que l'interpolation sur les déplacements (d'où le qualificatif « isoparamétrique » : mêmes paramètres), sont en particulier caractérisées par leur matrice jacobienne J (matrice des dérivées partielles du premier ordre de la transformation géométrique).

Élément isoparamétrique (quadratique) : repère local et global (2D)

Un des avantages de cette transformation est la possibilité d'approcher les frontières du domaine d'étude par des géométries polynomiales. Cependant, l'inconvénient principal est de complexifier le calcul des matrices et vecteurs élémentaires. Leur évaluation se fait alors par intégration numérique à l'aide d'un schéma de type quadrature : somme pondérée de valeurs en certains points situés à l'intérieur des éléments (poids ω_i, points ξ_i) caractéristiques de la méthode (méthode de Gauss, de Hammer…). Dans ce cas, le calcul de la matrice de rigidité élémentaire et du vecteur élémentaire des charges équivalentes (aux charges volumiques) s'écrivent :

Quand les matrices (faisant intervenir la matrice B) sont évaluées de cette façon, il est naturel de calculer aux mêmes points les contraintes et les déformations (également fonction de B). Les logiciels fournissent alors ces grandeurs aux points d’intégration (aux points de Gauss ou de Hammer selon le schéma d'intégration) ; de plus, parler de valeurs aux nœuds d'un maillage pour les contraintes ou les déformations n'a pas de sens, car il y a autant de valeurs en un nœud qu'il y a d'éléments concourants en ce point et il n'y a pas forcément continuité d'un élément à l'autre ! Vouloir à tout prix obtenir une valeur (de contrainte ou de déformation) en un nœud du maillage suppose donc un traitement : extrapolation à partir des valeurs aux points d'intégration et moyenne des valeurs provenant des différents éléments concourants, qui conduit au minimum à un lissage des résultats obtenus (cependant, parler de valeur aux nœuds d'un élément particulier est justifié). Quand ces quantités sont obtenues par intégration analytique (cas des poutres par exemples), elles peuvent être exprimées aux nœuds, au centre de gravité… les remarques précédentes sur la continuité entre éléments restant valables.

Remarque : certains éléments finis présents dans les codes sont dits sous-intégrés (ou à intégration réduite) : il s'agit d'éléments pour lesquels l'intégration numérique lors du calcul des quantités élémentaires (matrice de rigidité et vecteur des charges équivalentes) est effectuée avec moins de points que celui mathématiquement nécessaire ; cela peut améliorer dans certains cas la précision des calculs (en évitant les phénomènes de bloquage mais également conduire à des problèmes numériques (matrice singulière par exemple ou apparition de modes d’hourglass).

Éléments quadratiques et points d'intégration (points de Gauss)

La façon de prendre en compte les conditions aux limites en déplacement (encastrement, appuis…) peut entraîner des différences de résultats selon la technique numérique retenue, qui peuvent être résumées à deux : la méthode de pénalisation, simple à mettre en œuvre mais qui présente l'inconvénient d'être sensible à l'ordre de grandeur des termes de la matrice de rigidité et la méthode des multiplicateurs de Lagrange, insensible aux problèmes de précision évoqués pour la première mais augmentant légèrement la taille du système à résoudre.

la matrice de masse étant obtenue par assemblage des matrices de masse élémentaires  (masse volumique) et C étant la matrice d'amortissement (visqueux). Ceci nécessite la résolution d'un système d'équations différentielles par des méthodes numériques adéquates (cf. chapitre dynamique).