Un méli-mélo de mes activités de recherche

Je vais içi vous présenter quelques-unes de mes activités de recherche et oui car dans Ingénieur de Recherche il y a Recherche donc j'en fais encore un peu (enfin j'encadre des gens qui en font en me contentant de dire "ce serait bien de faire ..." 😁). Mais pour le détail vous pouvez consulter la liste des publications ici

Méthodes pour le HPC

Parmi mes activités se trouve le HPC. Grosso modo en mécanique des matériaux et des structures la tendance est clairement a faire des calculs toujours plus gros donc il faut pouvoir passer ces calculs.

Dans un calcul de structures par éléments finis on distingue deux gros postes de coût :

  • La résolution du système linéaire : intervenant au niveau de la méthode de Newton globale et donc implique la résolution d'un système de taille nombre de degrés de liberté fois nombre de degrés de liberté.
  • L'intégration de la loi de comportement : au niveau local (à chaque point d'intégration de chaque élément)

Selon les configurations le ratio entre ces deux postes de coût peut aller de 80%/20% (pour le cas du calcul élastique linéaire) à 20%/80% (pour le calcul considérant une loi de comportement très fortement non-linéaire).

Un des enjeux du HPC est donc de travailler sur le premier poste de coût avec pour objectif :

  1. Accélérer l'étape de résolution du système linéaire
  2. Sur une architecture donnée passé des calcul le plus gros possible.

Solveur itératifs avec préconditionneur pour les problèmes fortement non-linéaires

Parmi les approches pour accélérer l'étape de résolution du système linéaire une piste naturelle est d'utiliser un solveur itératif. En effet les solveur itératif ont le bon goût d'être très peu gourmand en ressource (mémoire vive) et sont d'un point de vu complexitée algorithmique beaucoup plus intéressants que des solveurs direct type LU ou Cholesky. Ca c'est en théorie car en pratique qui dit solveur itératif dit convergence. Or dans le cadre de calcul par éléments finis en considérant des lois de comportement non linéaire les solveur itératifs ont tendances à rencontrer quelques petits problèmes de convergence.

Une approche classique est alors d'introduire un préconditionneur permettant d'améliorer les propriétés du système à résoudre (conditionnement) et ainsi garantir une meilleure convergence. Il existe dans la litérature une grande variété de préconditionneur. Nous avons cependant développé un préconditionneur spécifique dédié au problème des mécanique non-linéaire des matériaux et des structures (Marchand et Quilici, 2019).

Dans les grandes lignes on introduit deux discrétisations : (i) un maillage fin, le maillage éléments finis sur lequel on souhaite déterminer la solution de notre problème ; (ii) un maillage grossier qui nous permettra de construire un préconditionneur peu cher et suffisament représentatif pour garantir la bonne convergence du solveur itératif utilisé au niveau du maillage fin. Enfin pour améliorer encore la convergence du solveur itératif on introduit en plus de la décomposition LU de la matrice tangante du maillage grossier une factorisation incomplète (ILU) de la matrice tangente du maillage fin.

Par exemple regardons la convergence d'un solveur itératif (Gradient Conjugué) sur un problème fortement non-linéaire (loi de comportement élasto-visco-plastique complexe) avec une zone de forte concentration de la plasticité. On considère deux préconditionneur différents : (i) notre préconditionneur grille grossière ; (ii) une factorisation incomplètre de cholesky (un standard). Pour chacune des configurations on trace ci-dessous la convergence du premier incrément de chargement et la convergence du pire incrément.

Calcul élasto-plastique grandes déformation d'un tube en torsion. Problème à 4 millions de degrés de liberté, résolution par GMRes avec préconditionneur grille grossière (problème grossier à 12 000 inconnues). Résolution en 48 heures de calcul sur 24 coeurs et 30 Go de RAM.

Solveur direct parallèle distribué

La vidéo ci-contre présente le résultat d'une simulation en grande déformation d'un simple essai de traction sur éprouvette plate avec un comportement élasto-plastique. La subtilité içi en terme de mécanique est que le formulation grandes déformation nous permet d'aller jusqu'à la striction que l'on voit clairement sur la fin de la vidéo.

L'autre particularité de ce calcul est sa taille, le problème est constitué de 10 millions de degrés de liberté et la discrétisation temporelle nécessite 222 incréments de chargement. La résolution est effectuée en utilisant la version distribuée du solveur direct MUMPS, en considérant un découpage en 48 sous-domaines et en utilisant 4 coeurs par sous-domaine. Cela fait donc un calcul sur 192 coeurs.

Toujours du calcul grandes déformations mais içi avec un modèle un peu plus compliqué, un modèle d'endommagement de type Gurson. Ce type de modèle permet de modéliser l'apparition et la croissance de pores dans la matières, pores assimilable à de l'endommagement. Le calcul présenté ici est un essai de traction sur une éprouvette SENT, le calcul comporte un peu plus de 3 millions d'inconnues nodales. Petite particularité le calcul est réalisé en utilisant des éléments mixtes déplacement, pression, variation de volume. Pourquoi est bien parce que sinon on a une contrainte qui clignotte et c'est pas vraiment mécanique ça 😅.

Calcul multi-échelle

Une autre de mes activité concerne le calcul multi-échelle. Il s'agit là surtout du travail d'un doctorant que j'encadre [](Mouad Fergoug). L'idée est que lorsque l'on souhaite réaliser un calcul éléments finis à l'échelle de la structure on ne peut généralement pas prendre en compte la microstructure du matériau. En effet cela conduirait à des coûts de calcul démentiels. C'est pour cela qu'on utilise généralement une loi homogène équivalente au niveau du calcul de structure. Pour déterminer cette loi homogène équivalente on utilise généralement un calcul sur un Volume Elémentaire Représentatif (VER) auxquel on impose des conditions de périodicité et que l'on vient charger en déformation suivant les 6 composantes (6 en 3D symétrique).

Ce que les gens savent un peu moins c'est qu'une fois le calcul homogène réalisé (calcul à l'échelle de la structure utilisant la loi homogène équivalente) il est possible de procéder à une étape de relocalisation des champs pour construire une estimation des champs microscopique, c'est-à-dire des champs à l'échelle de la microstructure. Pour avoir une estimation correcte au niveau des bords il est nécessaire de faire un traitement particulier voir (Fergoug et al)

Et ce que les gens savent encore moins c'est qu'en réalité en faisant la relocalisation on peut aller plus loin car quand on déroule les équations et bien en fait il y a plein de termes en plus qui apparaissent, c'est l'homogénéisation asymptotique 🤯. Et du coup si on prend en compte ces termes additionnel on peut monter aux ordres supérieurs et faire des trucs vachement sympa mais chut 🤫 le papier est soumi mais pas encore accepté.

Assimilation de données

Enfin une autre activité, qui fut historiquement ma première activité de recherche puisque c'était pendant ma thèse, concerne les méthodes d'assimilation de données et leurs applications aux problèmes d'identification/recalage en mécanique des matériaux et des structures. Par méthode d'assimilation de donnée j'entends principalement les méthodes de type filtre de Kalman. Petit rappel un filtre de Kalman permet à partir d'un modèle plus ou moins pertinent et d'observations de construire une estimation de l'état du système considéré.

Par exemple considérons le problème tout simple du tir ballistique et essayons d'estimer la position et vitesse du projectile à partir d'observation en considérant pour le filtre de Kalman un modèle faux (dans notre cas un modèle sans prise en compte du frottement de l'air). Ce que l'on peut alors observer c'est qu'en jouant un peu sur les paramètres du filtre de Kalman on arrive à obtenir une estimation plutôt satisfaisante des positions et vitesses du projectile même si notre modèle est pourri 😲 !

A partir de là on peut étendre l'utilisation des filtres de Kalman aux problèmes inverses. Pour cela une approche classique est d'étendre l'état de notre système en y ajoutant les paramètres à déterminer. Par exemple considérons le problème jouet d'aéroélasticité ci-dessous :

Nous allons chercher à déterminer la valeur du coefficient d'amortissement par filtre de Kalman et plus particulièrement en utilisant deux versions le Extended Kalman Filter et le Unscented Kalman Filter. Tout d'abord vous pouvez voir dans le graphique ci-dessous la réponse mécanique du problème direct.

Pour réaliser l'identification du coefficient d'amortissement (C_h) on mis en place une formulation étendu en intégrant le paramètre a identifier dans le vecteur d'état du système dynamique. Ci-dessous la précision de l'identification au cours du processus d'assimilation pour les deux approches Extended et Unscented Kalman Filter.