r7.20.01 Projection d’un champ sur un maillage#
Résumé:
La commande PROJ_CHAMP permet de « projeter » des champs (connus sur un maillage \(\mathit{ma1}\) ) sur un autre maillage \(\mathit{ma2}\) .
Dans ce document, on décrit les différentes méthodes de projection accessibles dans cette commande. Le paragraphe [§ 6 ] donne quelques éléments de validation de ces méthodes.
Méthode “COLLOCATION”#
Remarque sur le vocabulaire:
Le mot « projeter » est parfois ambigu dans ce document.
Quand on dit « projeter » le champ de \(\mathit{ma1}\) vers \(\mathit{ma2}\) , on cherche le champ sur \(\mathit{ma2}\) connaissant celui sur \(\mathit{ma1}\) : la projection va de 1 vers 2.
Pour la méthode “COLLOCATION”, il faut trouver pour chaque nœud \(\mathit{n2}\) de \(\mathit{ma2}\) quel est le point de \(\mathit{ma1}\) qui occupe la même position que \(\mathit{n2}\) , pour cela on projette le nœud \(\mathit{n2}\) sur le maillage \(\mathit{ma1}\) : la projection va de 2 vers 1.
Algorithme mis en œuvre#
On boucle sur tous les nœuds du maillage \(\mathit{ma2}\) :
pour chaque nœud (\(\mathit{n2}\) ), on procède en 4 étapes :
On cherche quelle est la maille \(\mathit{m1}\) de \(\mathit{ma1}\) qui « contient » géométriquement le nœud n2,
On détermine une position approchée de \(\mathit{n2}\) dans \(\mathit{m1}\) (i.e. ses coordonnées dans l’élément de référence associé à la maille \(\mathit{m1}\) ),
On «affine» la position calculée précédemment en résolvant un problème non linéaire pour tenir compte du fait que les mailles ne sont pas linéaires (le jacobien de la transformation géométrique varie sur la maille).
On utilise les fonctions de forme de la maille \(\mathit{m1}\) pour déterminer la valeur du champ sur \(\mathit{n2}\) connaissant la valeur du champ sur les nœuds de \(\mathit{m1}\) .
Remarque#
La troisième étape montre que cette méthode suppose que tous les nœuds de la maille m1 connaissent le champ à projeter. Par exemple, on ne saurait pas projeter un champ qui porterait des degrés de liberté différents sur ses nœuds « sommets » et sur ses nœuds « milieux d’arête ». La projection d’un tel champ est possible en revanche avec les 2 autres méthodes de projection [§ 3 ].
Les difficultés rencontrées et leur traitement#
Étape 1#
On ne traite pas l’éventuelle courbure des bords des éléments. Par exemple, dans la figure ci‑dessous (problème plan), le nœud \(\mathit{n2}\) sera déclaré appartenir à la maille \(\mathit{m1a}\) alors qu’il appartient à la maille \(\mathit{m1b}\) . En réalité, ce problème est rare car, le plus souvent, les bords «courbes» n’existent que sur la frontière des domaines maillés (les facettes internes sont planes).
Si le nœud \(\mathit{n2}\) est en réalité « extérieur » au maillage \(\mathit{ma1}\) , on lui affectera la maille \(\mathit{m1}\) la plus proche de lui. Ce comportement permet de projeter, sans s’arrêter en erreur fatale, un champ sur un maillage dont la frontière diffère légèrement de celle du maillage initial (ce qui est toujours le cas dans la pratique).
Étapes 2 et 3#
Pour trouver le point de l’élément de référence qui donne par la transformation géométrique le nœud \(\mathit{n2}\) , il faut en général résoudre un problème non-linéaire, car il n’y a que pour le triangle à 3 nœuds (en 2D) et pour le tétraèdre à 4 nœuds (en 3D) que la transformation géométrique est linéaire.
On suppose dans un premier temps que la maille est linéaire en oubliant ses nœuds milieux et en la découpant arbitrairement en triangles ou tétraèdres. Dans cette maille «linéarisée» on détermine une position approchée du point en vis à vis du nœud \(\mathit{n2}\) .
Puis on résout par un algorithme de Newton le problème non-linéaire de la recherche précise du point en partant de l’approximation précédemment calculée.
Remarque : le problème non-linéaire est résolu dans deux cas de figure 2D (projection d’un maillage 2D sur un autre maillage 2D) ou 3D (projection d’un maillage 3D sur un autre maillage 3D) mais pas pour les cas 2.5D (coques plongées dans l’espace 3D) ou 1.5D (maillage filaire en 2D),
Étape 4#
On n’utilise pas toujours les vraies fonctions de forme des éléments du maillage initial. En effet, dans le Code_Aster , ce sont les éléments finis qui choisissent leurs fonctions de forme : un triangle de thermique n’est pas obligé de choisir les mêmes fonctions de forme qu’un triangle de mécanique. Un élément peut aussi ne pas avoir besoin de fonctions de forme, ou bien il peut choisir des fonctions différentes selon les variables à interpoler. La formulation de l’élément ne fait pas non plus toujours apparaître d’élément de référence et de transformation géométrique associée.
Pour toutes ces raisons et pour que la programmation de PROJ_CHAMP soit indépendante des éléments finis présents dans le modèle, on affecte à toutes les mailles de \(\mathit{ma1}\) , les fonctions de forme des éléments iso-paramétriques 1D, 2D ou 3D [R3.01.01].
Détails concernant les étapes 1 et 2#
Informatiquement, les étapes 1 et 2 sont réalisées simultanément. Nous allons discuter successivement de la façon de traiter les 3 problèmes suivants :
traitement d’un nœud \(\mathit{n2}\) se retrouvant à l’intérieur de la frontière du maillage \(\mathit{ma1}\) (cas le plus fréquent),
traitement d’un nœud \(\mathit{n2}\) à l’extérieur de la frontière de \(\mathit{ma1}\) ,
traitement du cas des maillages de type « coque » (surfaces plongées dans \({\mathrm{ℝ}}^{3}\) ).
Nœud n2 « intérieur »#
Pour comprendre le traitement d’un nœud intérieur, prenons le cas d’une maille 2D QUAD8 (\(\mathit{abcd}\) ). On commence par oublier ses nœuds milieux (et donc leur éventuelle courbure) puis on la découpe en 2 triangles (\(\mathit{abc}\) et \(\mathit{acd}\) ). Ce découpage est arbitraire (il dépend de la numérotation locale des nœuds du QUAD8). Notons que l’autre découpage possible (autre diagonale) donnerait en général un autre point dans l’élément de référence.
\(\mathit{n2}\) appartient au triangle \(\mathit{abc}\) . On cherche ses coordonnées barycentriques dans ce triangle. Ce sont les 3 nombres \(\mathit{xa},\mathit{xb},\mathit{xc}\) tels que l’on puisse écrire : \(\mathit{n2}=\mathit{xa}\times a+\mathit{xb}\times b+\mathit{xc}\times c\)
Le point N2 de l’élément de référence retenu par l’algorithme sera:
\(\mathit{N2}=\mathit{xa}\times A+\mathit{xb}\times B+\mathit{xc}\times C\)
Pour les mailles volumiques (Hexaèdres, Pentaèdres, Pyramides et Tétraèdres), on procède de la même façon : on oublie les nœuds milieux et on découpe les mailles non tétraèdriques en tétraèdres.
Nœud n2 « extérieur »#
Un nœud \(\mathit{n2}\) est déclaré extérieur au maillage \(\mathit{ma1}\) si on n’a trouvé aucune maille pour laquelle il soit intérieur. Lorsque ceci est constaté, on recherche la maille \(\mathit{m1}\) la plus proche de \(\mathit{n2}\) . La distance calculée est celle qui sépare le nœud \(\mathit{n2}\) et la frontière de la maille \(\mathit{m1}\) . Appelons \(\mathit{p2}\) le point de \(\mathit{m1}\) le plus proche de \(\mathit{n2}\) . Ce point peut être sur une face d’un élément volumique ou sur une arête ou sur un sommet.
Exemple (en 2D) :
Au nœud \(B\) , on associe le point \(\mathit{pB}\) obtenu ici par projection orthogonale de \(B\) sur une arête de \(\mathit{T1}\) .
Au nœud \(A\) , on associe le point \(\mathit{pA}\) du triangle \(\mathit{T1}\) . On aurait pu tout aussi bien lui associer le point pA du triangle \(\mathit{T2}\) , mais cela n’aurait rien changé puisque le champ à projeter est connu aux nœuds du maillage. Il est donc continu entre les éléments adjacents.
Une fois trouvé le point \(\mathit{p2}\) de la maille \(\mathit{m1}\) qui réalise le minimum de la distance avec \(\mathit{n2}\) , on se ramène au problème précédent : le point \(\mathit{N2}\) de l’élément de référence qui sera associé à \(\mathit{n2}\) sera le correspondant de \(\mathit{p2}\) par la procédure du paragraphe [§ 2.3.1 ].
Remarques :
Pour un triangle (ou un tétraèdre) donné \(T\) , il n’y a qu’un point \(\mathit{p2}\) réalisant la distance minimale entre \(\mathit{n2}\) et \(T\) car \(T\) est convexe. Cette propriété disparaîtrait si on tenait compte de la courbure des bords des mailles. On voit là que les deux simplifications de l’implémentation (oubli des nœuds milieux et découpage des mailles en triangles linéaires) sont liées entre elles.
Un point extérieur aura toujours une valeur interpolée entre les valeurs des nœuds du maillage \(\mathit{ma1}\) et jamais extrapolée ; ce qui n’est pas le cas des 2 méthodes NUAGE_DEG_0/1.
Cas des maillages « coque »#
Lorsque l’on cherche à projeter les nœuds d’un maillage surfacique sur un autre maillage surfacique, on tombe en général systématiquement sur le cas des points « extérieurs » ci-dessus. En effet, l’imprécision sur les coordonnées des nœuds fait qu’un nœud \(\mathit{n2}\) n’est jamais rigoureusement dans le plan des triangles des mailles de \(\mathit{ma1}\) .
C’est donc la procédure du [§ 2.3.2 ] qui s’applique:
Pour chaque nœud \(\mathit{n2}\) :
recherche du triangle (ou du tétraèdre) qui réalise le minimum de distance avec \(\mathit{n2}\) . Identification du point \(\mathit{p2}\) qui réalise cette distance,
calcul du point \(\mathit{n2}\) de l’élément de référence qui correspond à \(\mathit{p2}\) par la procédure du [§3.3.1].
Méthode “NUAGE_DEG_0” ou “NUAGE_DEG_1”#
Principe de la méthode#
Ces 2 méthodes sont basées sur le même principe : on choisit a priori des fonctions de base \(\mathit{Fi}(x,y,z)\) (ici des polynômes de degré 0 ou 1). On cherche dans l’espace vectoriel engendré par ces fonctions de base, la fonction \(F=\sum({\alpha}_{i}{F}_{i})\) qui réalise la distance minimum (au sens des moindres carrés) avec le « nuage » des points connus. Une fois cette fonction trouvée, on l’évalue au point cherché.
Pour alléger les notations, on se place en 2D (mais les calculs peuvent être faits de la même façon en dimension 3 ou plus …). Le champ à projeter est un ensemble de couples \((\mathit{Xj},\mathit{Vj})\) où \(\mathit{Xj}=(\mathit{xj},\mathit{yj})\) est un nœud du maillage \(\mathit{ma1}\) et \(\mathit{Vj}\) est un réel (valeur du champ sur ce nœud). Ce champ constitue le nuage des points connus.
Soit un nœud \(\mathit{n2}\) \((x,y)\) du maillage \(\mathit{ma2}\) pour lequel on veut calculer la valeur du champ (\(V\) ).
Choix des fonctions de base :
NUAGE_DEG_0: |
une seule fonction : \(\mathit{F1}=1\) |
NUAGE_DEG_1: |
3 fonctions : \(\mathit{F1}=1\) ; \(\mathit{F2}=x\) ; \(\mathit{F3}=y\) |
Dans la suite de ce paragraphe, on choisira NUAGE_DEG_1 pour que les formules ne dégénèrent pas trop.
La fonction cherchée est \(F=\alpha 1\mathrm{F1}+\alpha 2\mathrm{F2}+\alpha 3\mathrm{F3}\) . On définit une sorte de « distance » entre \(F\) et les couples \((\mathit{Xj},\mathit{Vj})\) :
\(D=\sum(\mathit{wj}{(F({X}_{j})-{V}_{j})}^{2})\)
où les \(\mathit{wj}\) sont les poids affectés à chaque couple \((\mathit{Xj},\mathit{Vj})\) .
On veut minimiser \(D\) par rapport aux 3 variables \(\alpha 1,\alpha 2,\alpha 3\) . \(D\) est une fonction quadratique de \(\alpha 1,\alpha 2,\alpha 3\) . Minimiser \(D\) revient à annuler ses dérivées et donc à résoudre un système linéaire à 3 inconnues : \(\alpha 1,\alpha 2,\alpha 3\) .
Lorsque ce problème est résolu, on calcule :
\(V=\alpha 1+\alpha 2x+\alpha 3y\)
Choix de la pondération des points du nuage#
Toute l’astuce de ces méthodes se trouve dans le choix (difficile) des poids \(\mathit{wj}\) affectés aux points du nuage.
on décide a priori que le poids d’un point (\(\mathit{Xj}\) ) ne dépend que de la distance (\(d\) ) séparant ce point du nœud \(\mathit{n2}\) (pondération isotrope),
si \(\mathit{wj}\) est une constante, le problème ne dépend plus du nœud inconnu \(\mathit{n2}\) . La fonction \(F\) est unique (pour tous les nœuds de \(\mathit{ma2}\) ) : c’est « la droite de régression linéaire » du nuage,
si \(\mathit{wj}(d)\) est une fonction trop peu décroissante, la fonction \(F\) est trop « lissée » : les « accidents » locaux sont « gommés » par le grand nombre de points lointains pris en compte,
si \(\mathit{wj}(d)\) est une fonction trop décroissante, on prend le risque de n’attraper aucun point du nuage. La conséquence numérique est que le système linéaire à résoudre devient singulier (et donc insoluble).
Nous avons choisi d’écrire \(w(d)\) comme un exponentielle décroissante paramétrée par 2 paramètres : \(\mathit{dref}\) et \(\beta\) :
\(w(d)={e}^{-{(d/\mathit{dref})}^{2\beta }}\)
\(\mathit{dref}\) est une distance de référence (dépendant du nœud \(\mathit{n2}\) ). Nous allons voir ci-dessous comment elle est calculée. \(\beta\) est une constante choisie pour annuler plus ou moins rapidement le poids des points distants de \(\mathit{n2}\) . Dans le code, \(\beta\) a été choisi à 0.75.
\(\mathit{dref}\) est la distance à partir de laquelle on souhaite voir le poids des points diminuer rapidement. Dans la programmation, \(\mathit{dref}\) est calculée comme le produit d’une distance \(\mathit{d1}\) par un coefficient \(\mathit{C1}\) . Aujourd’hui, nous avons choisi \(\mathit{C1}=0.45\) .
\(\mathit{d1}\) est défini comme suit:
en 3D, \(\mathit{d1}\) est le rayon de la plus petite boule de centre \(\mathit{n2}\) qui contient 4 points du nuage non co-planaires,
en 2D, \(\mathit{d1}\) est le rayon de la plus petite boule de centre \(\mathit{n2}\) qui contient 3 points du nuage non alignés.
Le problème est dit 2D si tous les nœuds du maillage \(\mathit{ma1}\) ont la même coordonnée \(z\) , il est dit 3D sinon.
Restriction actuelle de la programmation#
Les 2 méthodes de projection NUAGE_DEG_0/1 ne sont programmées que pour les champs réels ( et non pour les champs complexes)
Méthode “ECLA_PG”#
Cette méthode est dédiée à la projection des champs aux points de Gauss (ELGA).
Illustrons cette méthode sur un exemple 2D.
On suppose que le maillage initial est constitué de 3 triangles (verts) ayant 3 points de Gauss (voir figure ci-dessous).
Chaque triangle est découpé en autant de sous-éléments que de points de Gauss (ici 3). La valeur portée par le point de Gauss est attribuée à tout le sous-élément.
Une maille du maillage final (parallélogramme) est dessinée en jaune clair. Elle a quatre points de Gauss (petits ronds oranges).
La valeur du champ projeté sur chacun de ces quatre points de Gauss sera celle du sous-élément en vis à vis du point. Dans notre exemple, les 2 points supérieurs appartiennent au même sous-élément et ils porteront la même valeur à l’issue de la projection.
Méthode “COUPLAGE”#
Cette méthode est dédiée au couplage fluide-structure avec Code_Saturne pour le domaine fluide. Elle est basée sur la méthode “COLLOCATION” en tenant compte de spécificités des problèmes fluide-structure envisagés ( cf. [bib2]).
On cherche donc à projeter des données d’un maillage fluide sur un maillage structure, et inversement.
Le maillage structure s’appuie sur une discrétisation EF, donc avec la connaissance de fonctions de formes associées, alors que le maillage fluide est basé sur une discrétisation en volumes finis(VF), donc sans fonctions de formes a priori . Les champs cinématiques sur le domaine fluide VF sont définis aux nœuds sommets alors que leurs champs duaux (efforts) sont constants par faces et définis au milieux des faces.
De plus, le maillage structure est très généralement nettement plus grossier que le maillage fluide (qui doit être très raffiné à la paroi afin de bien représenter la couche limite).
La première étape d’appariement des nœuds du maillage fluide sur le solide est faite par une méthode de type “COLLOCATION” classique. Il faut signaler que cette méthode ne peut être utilisée que pour générer la structure de données matrice de correspondance de maillage (CORRESP_2_MAILLA cf. documentation U4.72.05). Cette structure de données sera ensuite réutilisée par les opérateurs d’échanges aux interfaces pour le couplage fluide-structure:
ENV_CINE_YACS(pour envoyer les champs de déplacement et vitesse au code fluide, cf. documentation U7.06.11),
MODI_CHAR_YACS(pour récupérer les efforts fluides et les sommer au chargement structure existant, cf. documentation U7.06.22).
En dehors de cette limitation, la principale différence avec la méthode “COLLOCATION” sera dans la construction de la structure de données matrice de correspondance (CORRESP_2_MAILLA confer documentation U4.72.05) en sortie de PROJ_CHAMP.
En effet, on va chercher à respecter au mieux la conservation de l’énergie à l’interface, donc la conservation du travail des efforts discrétisés lors du passage fluide vers solide ( cf. [bib2]). Bien évidemment, ce résultat n’est valable que sous des hypothèses particulières:
champs d’effort fluide constant par élément,
maillage fluide plus fin que le maillage solide,
méthode d’appariement de type “COLLOCATION” de Code_Aster .
En pratique, pour tout problème couplé fluide-structure, on va devoir utiliser deux fois la méthode “COUPLAGE’afin de calculer deux structures de données de correspondance de maillage. En effet, on doit gérer séparément la projection des champs cinématiques qui s’appuie sur des champs fluides définis aux nœuds sommets du maillage fluide, de celle des efforts qui sont définis aux milieux des faces fluides. Dans chacun de ces deux cas, le maillage fluide correspondant sera soit le maillage des nœuds sommets, soit celui des nœuds milieux des faces. La méthode de projection est ensuite la même et le maillage structure est, lui, unique. Tous les champs fluides (déplacements, vitesses et efforts à l’interface) seront considérés par Code_Aster comme des champs aux nœuds classiques.
Éléments de validation pour la projection des champs aux nœuds#
Pour valider les 3 méthodes proposées, nous allons traiter l’exemple de la projection d’un champ de température « unidimensionnel » \(T=T(x)\) .
Le champ à projeter vaut:
\(T(x)=\sin(\mathrm{3x})+\mathrm{Heavyside}(x-1)\)
Ce champ est affecté sur les nœuds d’un maillage \(\mathit{ma1}\) linéaire assez grossier (14 mailles) du segment \([0,2]\) . Ce champ possède la propriété d’être discontinu (en théorie) au point \(x=1\) . Du fait de la discrétisation sur le maillage \(\mathit{ma1}\) , le champ semble seulement varier très brutalement entre les 2 points \(x=0.99\) et \(x=1.01\) .
On projette ce champ sur un maillage \(\mathit{ma2}\) très fin du même segment (300 éléments de longueur \(2/300\) ).
Méthode : “COLLOCATION”#
Collocation
On constate qu’avec cette méthode, le champ projeté est sans surprise : la valeur obtenue par projection est toujours l’interpolée linéaire entre les 2 nœuds de \(\mathit{ma1}\) où le champ est connu.
Maillage initial quadratique#
Si on refait le même calcul en remplaçant le maillage linéaire \(\mathit{ma1}\) par un maillage quadratique (contenant environ 2 fois moins de mailles), on trouve :
Collocation
On constate que l’interpolation du champ est maintenant parabolique entre les nœuds de \(\mathit{ma1}\) . On approche donc mieux la forme de la fonction initiale.
Méthode : “NUAGE_DEG_0”#
On constate qu’avec cette méthode, le champ projeté se présente comme une succession de petits « paliers » horizontaux reliés entre eux. Cet aspect en escalier est lié aux paramètres (\(\beta\) et \(\mathit{C1}\) ) choisis en « dur » dans le code pour la forme de la décroissance exponentielle du poids des points. On voit également que la discontinuité du champ initial est très fortement gommée.
Influence de la forme des exponentielles décroissantes#
Sur la figure suivante, nous avons modifié un paramètre de l’exponentielle décroissante de la méthode “NUAGE_DEG_0” : le paramètre \(\mathit{dref}\) (ou \(\mathit{C1}\) ce qui revient au même) a été multiplié par 0.5 ou 2. par rapport à la valeur retenue par le code.
On constate que le choix de \(\mathit{dref}\) influe beaucoup sur le résultat. S’il est trop grand, on ne voit plus la discontinuité. S’il est trop petit, la forme en marche d’escalier s’accentue.
Méthode : “NUAGE_DEG_1”#
On constate qu’avec cette méthode, le champ projeté est correct dans les 2 zones régulières (\(x<0.99\) et \(x>1.01\) ).
En revanche la discontinuité du champ entre ces 2 valeurs est fortement exagérée. Cet exemple illustre bien la faculté qu’à cette méthode à extrapoler les valeurs des points initiaux (du fait de l’estimation du gradient du champ).
Dans les 2 parties régulières, le champ projeté est assez voisin de celui obtenu avec la méthode “COLLOCATION”. On remarque simplement que la méthode “NUAGE_DEG_1” arrondit un peu les angles
Influence de la forme des exponentielles décroissantes#
Sur la figure suivante, nous avons modifié un paramètre de l’exponentielle décroissante de la méthode “NUAGE_DEG_1” : le paramètre \(\mathit{dref}\) (ou \(\mathit{C1}\) ce qui revient au même) a été multiplié par 0.5 ou 2. par rapport à la valeur retenue par le code.
Là encore, on constate que le choix de \(\mathit{dref}\) est crucial pour le résultat : trop grand, la discontinuité est gommée, trop petit, la discontinuité est exagérée.
Choix de la meilleure méthode de projection pour les champs aux nœuds#
Au vu des quelques courbes précédentes, il parait clair que la méthode “COLLOCATION” est en général préférable aux méthodes NUAGE_DEG_0/1. Cette méthode est « naturelle » dans le cadre des éléments finis et elle ne dépend d’aucun coefficient numérique d’ajustement.
Les méthodes NUAGE_DEG_0/1 doivent être réservées, à notre avis, pour des utilisations un peu spéciales:
cas d’un maillage \(\mathit{ma1}\) « inexistant » : on ne dispose que des nœuds mais pas des mailles (par exemple, les « nœuds » de \(\mathit{ma1}\) sont en fait des capteurs de mesure),
cas d’un champ (résultat d’un calcul ou obtenu par des mesures) que l’on veut lisser volontairement. Mais dans ce cas, il faudrait que les 2 paramètres numériques (\(\mathit{C1}\) et \(\beta\) ) soient accessibles à l’utilisateur ce qui n’est pas le cas aujourd’hui.
Bibliographie#
VAUTIER: « Éléments iso-paramétriques » , Documentation de Référence du Code_Aster n°[R3.01.00]
MAMAN and C. FARHAT: «Matching computations : a parallel approach.» Computers & Structures, vol. 54 (4), p. 779-785, 1995.