r7.02.12 eXtended Finite Element Method : Généralités#
Résumé:
Ce document présente la méthode X-FEM (eXtended Finite Element Method) qui permet principalement de considérer des fissures ne respectant pas le maillage pour traiter les problèmes de fissures 2D et 3D. La fissure est définie par la commande DEFI_FISS_XFEM [u4.82.08] et est utilisable pour des calculs en statique linéaire et non-linéaire.
D’autres documentsdédiés à des problématiques spécifiques sont disponibles:
contact-frottement sur les lèvres de la fissure en petits glissements avec X-FEM [r5.03.54],
contact-frottement sur les lèvres de la fissure en grands glissements avec X-FEM [r5.03.53],
Représentation de la fissure par des « level sets »#
La méthode des «Level sets» a été introduite dans le cadre de la mécanique des fluides pour représenter l’évolution d’interfaces. L’idée principale est de considérer l’interface comme l’iso-zéro d’une fonction distance. Le choix de la fonction distance importe peu ici, car seule la connaissance de l’iso-zéro est utile et importante.
Soit une interface \(\Gamma\) délimitant un ouvert \(\Omega\) de \({\Re }^{n}\) . L’idée est de définir une fonction \(\varphi (x,t)\) régulière (au moins Lipchitzienne) telle que le sous-espace \(\varphi (x,t)=0\) représente l’interface.
Cette méthode s’applique aisément aux problèmes de fissuration 2D, notamment dans le cadre des approches où la fissure n’est pas maillée. ([bib14], [bib15] en 2D). L’extension est possible pour le traitement des fissures en 3D.
Ainsi, dans le cas de la fissuration, il est nécessaire d’introduire deux level sets ([bib17] en 2D et [bib18] en 3D):
une level set normale (\(\text{lsn}\) ) qui représente la distance à la surface de la fissure (surface étendue par prolongement à tout le domaine),
une level set tangente (\(\text{lst}\) ) qui représente la distance au fond de fissure.
L’iso-zéro de la level set normale définit la surface de la fissure, étendue par continuité à tout le domaine. L’intersection des iso-zéros des deux level sets définit le fond de fissure. De plus, le signe de la level set tangente est choisi de telle sorte que la surface de la fissure \({\Gamma}_{\text{cr}}\) corresponde à l’espace engendré par \((\text{lsn}=0)\cap (\text{lst}<0)\) . Le signe de la level set normale est choisi arbitrairement grâce à la convention d’orientation de la normale au plan de fissure, définie explicitement à la r7.02.12-calcul-level-set-normale-lsn. Les points \(x\) pour lesquels \(\text{lsn}(x)\) est négatif sont dits «au-dessous» de la fissure, et ceux pour lesquels \(\text{lsn}(x)\) est positif sont dits «au-dessus» de la fissure (voir Fig. 336 et Fig. 337).
Fig. 336 Level sets pour la représentation d’une fissure 3D#
Fig. 337 Level sets pour la représentation d’une fissure 2D#
La phase de propagation de la fissure se traduit simplement par la propagation des level sets. La propagation d’une level set nécessite trois étapes successives [bib13]:
extension de la vitesse connue sur l’iso-zéro vers le domaine entier,
propagation de la level set à partir de ce champ de vitesse,
réinitialisation de la fonction level set afin de conserver une fonction distance signée.
La propagation d’une fissure représentée par 2 level sets présente quelques particularités. Le vecteur vitesse est connu seulement sur le front de fissure, c’est-à-dire une courbe. La fissure ne peut que grandir, pas se déplacer. Deux fonctions level sets doivent être propagées et on souhaite que leurs gradients restent orthogonaux. L’enchaînement des étapes peut être résumer comme ceci:
propagation de \(\text{lsn}\) et réinitialisation de \(\text{lsn}\) ,
propagation de \(\text{lst}\) ,
orthogonalisation du gradient de \(\text{lst}\) versus gradient de \(\text{lsn}\) ,
réinitialisation de \(\text{lst}\) .
Ces étapes se ramènent toutes à la résolution d’équation de type Halimton-Jacobi [bib12]:
où \(F\) et \(f\) sont des fonctions qui dépendent de l’étape.
La Fast Marching Method [bib16] est une technique alternative bien adaptée à la propagation strictement monotone de fronts. Cette méthode sépare les nœuds du maillage suivant leur éloignement de l’interface, et à chaque itération l’équation de propagation est résolue pour les seuls nœuds immédiatement adjacents à l’interface, en utilisant un schéma de différences finies d’ordre 2.
Cette partie est plus détaillée dans le document r7.02.13.
Le calcul des level sets (champs scalaires) est effectué pour chaque fissure. Il peut se faire de deux manières. Soit par la donnée de leurs expressions analytiques et dans ce cas, une simple évaluation de ces fonctions aux nœuds du maillage fournit les champs scalaires recherchés. Soit la fissure est maillée et dans ce cas, il est nécessaire de donner des mailles surfaciques correspondants à une lèvre (GROUP_MA_FISS) et les mailles linéiques correspondants au fond de fissure (GROUP_MA_FOND). ). Dans le cas 2D, on donnera des mailles linéiques (pour GROUP_MA_FISS) et des mailles points (pour GROUP_MA_FOND). Les distances sont ensuite calculées par un algorithme de projections orthogonales inspiré de celui utilisé pour le contact [bib19], explicité aux r7.02.12-calcul-level-set-normale-lsn et r7.02.12-calcul-level-set-tangente-lst.
Pour chaque nœud du maillage, on cherche la maille de GROUP_MA_FISS la plus proche de ce nœud. Pour cela, on utilise les algorithmes de projection sur un triangle (voir le paragraphe §2.3.2 de [bib19]) et sur un quadrangle (voir le paragraphe §2.3.3 de [bib19]). La valeur de \(\text{lsn}\) est alors la distance normale de ce point à la maille.
En 3D, il faut faire attention au sens de la normale. En effet, les mailles de GROUP_MA_FISS étant intérieures à la structure, ce ne sont pas des mailles de bord, et l’orientation automatique des normales (mot-clé ORIE_PEAU_3D de l’opérateur MODI_MAILLAGE) est alors impossible. Pour s’assurer de toujours choisir le même sens pour les normales, on prend la normale de la première maille de GROUP_MA_FISS comme référence et on «propage» la direction de cette normale à toutes les autres mailles du groupe pour chacun des éléments adjacents. Pour chaque maille on stocke un indicateur d’orientation de la maille qui vaut «+1» si l’orientation de sa normale est cohérente avec celle de la maille de référence construite de proche en proche et «-1» sinon. Cela permet d’affecter le bon signe de la valeur de la level set normale en multipliant sa valeur non corrigée par l’indicateur d’orientation (voir algorithme de calcul de \(\text{lsn}\) en 3D ci-dessous).
Algorithme de calcul de l’indicateur d’orientation de la maille:
création du vecteur \(I\) «indicateur d’orientation» et du vecteur \(C\) contenant le numéro de couche à laquelle chaque élément de GROUP_MA_FISS appartient. La taille des deux vecteurs est égale au nombre d’éléments \(n\) de GROUP_MA_FISS
on met à zéro tous les éléments des deux vecteurs \(I\) et \(C\) , sauf le premier élément auquel on affecte la valeur +1
boucle \(i\) sur le numéro de la couche (vecteur \(C\) ), de 1 à \(n\)
boucle \(j\) sur les éléments du vecteur \(C\) , de 1 à \(n\)
on récupère le numéro de la couche à laquelle l’élément \(j\) appartient: \({n}_{\mathit{couche}}=C(j)\)
si \({n}_{\mathit{couche}}=i\) , c’est-à-dire l’élément \(j\) appartient à la couche courante \(i\) :
on calcule la normale \(\vec{{n}_{j}}\) à l’élément \(j\) (on calcule la normale au triangle formé par les trois premiers nœuds définissant l’élément)
boucle \({n}_{\mathit{el}}\) sur les éléments qui ont au moins un nœud en commun avec l’élément \(j\)
si l’élément \({n}_{\mathit{el}}\) n’appartient à aucune couche (\(C({n}_{\mathit{el}})=0\) ), on lui assigne le numéro de la couche courante+1 (\(C({n}_{\mathit{el}})=i+1\) ), on en calcule la normale \(\vec{{n}_{{n}_{\mathit{el}}}}\) et l’indicateur d’orientation: \(I({n}_{\mathit{el}})=\mathit{sign}\left(I(j)\cdot \vec{{n}_{{n}_{\mathit{el}}}}\cdot \vec{{n}_{j}}\right)\)
fin boucle
fin si
fin boucle
fin boucle
Algorithme de calcul de \(\text{lsn}\) en 3D:
boucle sur tous les nœuds \(P\) du maillage
initialisation de \(\mathit{dmin}\)
boucle sur les mailles triangulaires de GROUP_MA_FISS (avec subdivision des quadrangles en triangles)
soient \(A\) , \(B\) et \(C\) les sommets du triangle
calcul de la normale à la maille: \(\overrightarrow{N}=\overrightarrow{\text{AB}}\wedge \overrightarrow{\text{AC}}\)
calcul de \(M\) , projeté de \(P\) dans le plan \(\mathit{ABC}\)
si \(M\) en dehors du triangle \(\mathit{ABC}\) , on ramène \(M\) sur une des droites \((\mathit{AB})\)\((\mathit{AC})\)\((\mathit{BC})\)
si M encore en dehors du triangle \(\mathit{ABC}\) , on ramène \(M\) en \(A\)\(B\) ou \(C\)
on récupère l” indicateur d’orientation\(I\) de la normale de la maille
si \(\mathrm{PM}<\mathrm{dmin}\) alors \(\mathrm{dmin}=\mathrm{PM}\) et \(\text{lsn}(P)=I\cdot \overrightarrow{\text{PM}}\cdot \overrightarrow{N}\)
fin boucle
fin boucle
La subdivision d’un quadrangle en triangles permet de se ramener à un calcul approché de projection dans un cas linéaire (car les fonctions de forme du quadrangle sont bilinéaires alors que celles du triangle sont linéaires). De plus, il existe 2 façons de réaliser un tel découpage (suivant la diagonale choisie). On engendre donc tous les cas possibles, comme il est fait au paragraphe §2.3.3 de [bib19].
Le calcul de la projection dans le cadre de la recherche d’appariement des nœuds de contact pour la méthode continue n’utilise pas ce genre d’algorithme, mais résout le problème non-linéaire de projection sur un quadrangle par la méthode de Newton. On pourrait s’en inspirer pour la calcul des level sets par projection.
En 2D, on commence par réorganiser les mailles de GROUP_MA_FISS pour avoir une série de mailles contiguës. Deux vecteurs vont être stockés, la liste des mailles ordonnées et la liste de leur orientation, qu’on appellera ici \(\mathit{LIMA}\) et \(\mathit{ORI}\) .
Algorithme de tri de GROUP_MA_FISS:
Initialisation de la première maille M, première maille de GROUP_MA_FISS. On a \(\mathit{LIMA}(1)=M\) et \(\mathit{ORI}(1)=1\) . On initialise le booléen \(\mathit{FINFIS}\) qui teste si on est en fin de fissure, \(\mathit{FINFIS}\) initialisé à \(\mathit{FALSE}\) .
boucle de recherche de la maille suivante \(\mathit{LIMA}(k)\) tant que \(\mathit{FINFIS}=\mathit{FALSE}\)
si \(\mathit{ORI}(k-1)=1\) , calcul de \(N\) nœud 2 de \(\mathit{LIMA}(k-1)\)
si \(\mathit{ORI}(k-1)=0\) , calcul de \(N\) nœud 1 de \(\mathit{LIMA}(k-1)\)
\(\mathit{FINFIS}=\mathit{TRUE}\)
boucle sur les mailles \(M\) de GROUP_MA_FISS, tant que \(\mathit{FINFIS}=\mathit{TRUE}\)
Calcul du nœud 1 et du nœud 2 de \(M\) , \(\mathit{N1}\) et \(\mathit{N2}\)
si \(\mathit{N1}=N\) , \(\mathit{ORI}(k)=1\) et \(\mathit{LIMA}(k)=M\) , \(\mathit{FINFIS}=\mathit{FALSE}\)
si \(\mathit{N2}=N\) , \(\mathit{ORI}(k)=0\) et \(\mathit{LIMA}(k)=M\) , \(\mathit{FINFIS}=\mathit{FALSE}\)
(sortie de boucle quand on trouve une maille suivante)
fin boucle
(sortie de boucle si on n’a pas trouvé de maille suivante car on est en fin de fissure)
La première maille identifiée pouvant être quelconque n’importe quelle maille à l’intérieur de la fissure, on effectue un décalage de la liste trouvée du nombre total de maille \(\mathit{NBMAF}\) moins le nombre de mailles trouvées \(k\) de façon à avoir une liste qui termine sur la maille de fin de fissure trouvée. On obtient donc un vecteur avec des composantes non affectées en tête qui représentent les mailles non encore répertoriées et une partie remplie sur la fin du vecteur correspondant aux mailles orientés précédemment stockées.
boucle de recherche de la maille précédente \(\mathit{LIMA}(k)\) tant que \(\mathit{FINFIS}=\mathit{FALSE}\) (\(k\) décroissant)
si \(\mathit{ORI}(k+1)=1\) , calcul de \(N\) nœud 1 de \(\mathit{LIMA}(k+1)\)
si \(\mathit{ORI}(k+1)=0\) , calcul de \(N\) nœud 2 de \(\mathit{LIMA}(k+1)\)
\(\mathit{FINFIS}=\mathit{TRUE}\)
boucle sur les mailles \(M\) de GROUP_MA_FISS, tant que \(\mathit{FINFIS}=\mathit{TRUE}\)
Calcul du nœud 1 et du nœud 2 de \(M\) , \(\mathit{N1}\) et \(\mathit{N2}\)
Si \(\mathit{N1}=N\) , \(\mathit{ORI}(k)=0\) et \(\mathit{LIMA}(k)=M\) , \(\mathit{FINFIS}=\mathit{FALSE}\)
Si \(\mathit{N2}=N\) , \(\mathit{ORI}(k)=1\) et \(\mathit{LIMA}(k)=M\) , \(\mathit{FINFIS}=\mathit{FALSE}\)
(sortie de boucle quand on trouve une maille précédente)
fin boucle
(sortie de boucle si on n’a pas trouvé de maille précédente, car on est en fin de fissure)
Pour chaque nœud du maillage, on cherche la maille de GROUP_MA_FOND la plus proche de ce nœud. Pour cela, on utilise l’algorithme de projection sur un segment (voir le paragraphe §2.3.1 de [bib19]). La valeur de \(\text{lst}\) est alors la distance normale de ce point au segment.
De même que précédemment, la détermination de la normale au segment n’est pas évidente a priori. Pour la calculer, il faut d’abord trouver la maille surfacique de GROUP_MA_FISS qui le borde.
Algorithme de calcul de \(\text{lst}\) en 3D:
boucle sur tous les nœuds \(P\) du maillage
initialisation de \(\mathit{dmin}\)
boucle sur les segments de GROUP_MA_FOND
soit \(A\) et \(B\) les deux extrémités du segment
boucle sur toutes les mailles de GROUP_MA_FISS
si \(A\) et \(B\) appartiennent à cette maille alors
soit \(C\) un nœud de la maille autre que \(A\) et \(B\)
calcul de la normale à cette maille
(c’est à dire la normale au plan de fissure): \(\overrightarrow{N}=\overrightarrow{\text{AB}}\wedge \overrightarrow{\text{AC}}\)
calcul de la normale au fond de fissure dans le plan de la fissure:
Soit \(M\) la projection de \(P\) sur \((\mathit{AB})\)
Si \(M\) est en dehors de \([\mathit{AB}]\) , on le ramène en \(A\) ou \(B\)
Si \(\mathit{PM}<\mathit{dmin}\) alors \(\mathit{dmin}=\mathit{PM}\) et \(\text{lst}(P)=\overrightarrow{\text{PM}}\cdot \overrightarrow{N'}\)
fin si
fin boucle
fin boucle
fin boucle
Remarque :
Lors de la projection de\(P\)sur les segments\(\mathit{AB}\)du fond de fissure, il est possible que tous les projetés\(M\)se retrouvent hors du segment considéré. Dans ce, cas, le fait de ramener\(M\)sur le bord n’est pas un critère suffisant pour déterminer le segment le plus proche (voirFig. 341). On propose de choisir le segment le plus proche comme celui où on a le moins rabattu le projeté. Avec les notations ci-dessus, si on appelle\(M’\)le point rabattu (point\(B\)sur laFig. 341, alors on cherche le segment pour lequel l’angle\(\alpha =\mathit{MPM}’\)est le plus petit. On note que dans un cas sans rabattement (quand le projeté tombe dans le segment), cet angle est nul. Ce test supplémentaire est introduit pour le choix du segment le plus proche.
Donc, quand on calcule la quantité \(d=\mathit{lstm}/{\mathit{lst}}_{max}\) , on a un problème.
Cette situation peut être fréquente en 2D avec une fissure «segment» (FORM_FISS=”SEGMENT”) et un maillage caractérisé par un fort gradient de taille de maille (et ce d’autant plus si le maillage est triangulaire et libre).
Pour corriger, avant de calculer \(d\) , on regarde si \({\mathit{lst}}_{max}\) est nulle.
Si oui, et si les conditions particulières sont respectées:
tous les nœuds sommets doivent vérifier \(\left|\mathit{lst}\right|<\epsilon\) ;
les level-sets doivent avoir été calculées depuis le catalogue de formes géométriques de DEFI_FISS_XFEM dans le cas 2D avec FORM_FISS=”SEGMENT”
pour chaque arête de la maille courante, on compare la longueur de son projeté orthogonal sur le segment (qui constitue la fissure) a la longueur de ce segment. Au moins un de ces projetés doit avoir une longueur comparable (critère relatif) à celle du segment. De cette manière on s’assure que la maille se trouve «loin» de la fissure (a moins que le maillage ne soit extrêmement grossier)
Alors, on court-circuite la fin de l’itération courante de la boucle sur les arêtes (celle du bloc relatif aux ajustements spécifique au quadratique). Sinon, on s’arrête en erreur fatale.
Quelle que soit la méthode de calcul utilisée (par des fonctions analytiques ou par projection), les champs des level sets sont interpolés par les fonctions de forme linéaires utilisées pour l’approximation du champ de déplacement [bib17]:
où \({\phi }_{i}\) sont les fonctions de forme linéaires classiques et \({\text{lsn}}_{i}\) et \({\text{lst}}_{i}\) les valeurs nodales des champs level sets. Sur un tétraèdre à 4 nœuds la surface de la fissure sera donc un plan, mais sur un hexaèdre à 8 nœuds, elle pourra être légèrement courbe.
L’erreur de discrétisation des level sets est directement liée à la taille du maillage et à la courbure de la level set. Examinons l’exemple suivant, mettant en jeu une fissure courbe en 2D:
Soit la level set ellipse d’équation:
\({x}^{2}+{(\frac{y}{0.5})}^{2}=1\)
La distance à l’ellipse est calculée en chacun des 4 nœuds du quadrangle défini par \((x,y)\in \left[0.4,0.9\right]\times \left[0,0.5\right]\) . Soit \(P\) un point donné de l’espace de cordonnées \(({x}_{p},{y}_{p})\) . Soit \(H\) sa projection sur l’ellipse précédemment définie. \(H\) a pour coordonnées \((a\cos\theta ,b\sin\theta )\) . L’équation de la droite \((\text{HP})\) est la suivante:
où \({\varphi}_{i}i=1,4\) sont les fonctions de forme associées aux nœuds du quadrangle. Ces fonctions s’exprime à l’aide des \({N}_{i}\) , fonctions de forme classiques sur le quadrangle de référence, et des changements de variables entre les coordonnées réelles \(\left(x,y\right)\) et les coordonnées de référence \(\left(s,t\right)\) .
Sur la Fig. 343, on a représenté en pointillé la projection sur l’ellipse de chaque nœud . On observe que l’iso-zéro de la level set interpolée avec les fonctions de forme du quadrangle est assez éloignée de l’ellipse initiale.
Fig. 343 Erreur de discrétisation de la level set#
Avec une interpolation quadratique, sur un élément quadrangle QUAD8, l’approximation de la level set s’écrit:
où \({\varphi}_{i}i=1,8\) sont les fonctions de forme associées aux nœuds du quadrangle. Ces fonctions s’expriment à l’aide des \({N}_{i}\) , fonctions de forme classiques sur le quadrangle de référence, et des changements de variables entre les coordonnées réelles \(\left(x,y\right)\) et les coordonnées de référence \(\left(s,t\right)\) .
On compare alors l’erreur de discrétisation entre l’interpolation linéaire et l’interpolation quadratique Fig. 344: à raffinement égal, les éléments quadratiques sont plus adaptés pour décrire une iso-zéro de forme courbe .
Fig. 344 Erreur d’interpolation l’iso-zero: quadratique vs linéaire#
Afin de limiter les problèmes d’intégration lorsque la fissure passe «près» d’un nœud, une procédure de réajustement de la level set normale est mise en place. Si sur une arête du maillage, \(\text{lsn}\) s’annule trop «près» d’un nœud extrémité, la valeur de \(\text{lsn}\) en ce nœud est mise à zéro. Le critère utilisé pour l’instant est 1% de la longueur de l’arête. Cette valeur est celle utilisée par le logiciel développé par l’équipe de Nicolas Moës au GeM (voir le paragraphe §1.3.2.3 de [bib20]).
Algorithme de réajustement de la level set normale:
boucle sur les mailles
boucle sur les arêtes de la maille, d’extrémités \(A\) et \(B\)
si \(\mid d\mid \le 0.01\) alors \(\text{lsn}(A)=0\) fin si
si \(\mid d-1\mid \le 0.01\) alors \(\text{lsn}(B)=0\) fin si
fin si
si élément quadratique
si \(\text{lsn}(A)=0\) et \(\text{lsn}(B)=0\) alors
\(d=\text{lsn}(M)\)
si \(\mid d\mid \le 0.0001\) alors \(\text{lsn}(M)=0\)
fin si
fin si
fin si
fin boucle
fin boucle
Dans le cas d’une interpolation quadratique, un réajustement complémentaire est nécessaire afin de limiter et contrôler les situations d’annulation multiple de la level set le long d’une arête. En effet, pour ne pas multiplier les configurations de découpe des éléments enrichis en sous éléments d’intégration, on autorise seulement deux situations d’annulation multiple le long de l’arête d’un élément enrichi. La première correspond au cas où la level set est nulle tout au long de l’arête (donc nulle en ses deux nœuds extrémités et en son nœud milieu). La deuxième est la situation où la level set est nulle sur l’un des nœuds extrémité et le nœud milieu (voir Fig. 345) . Dans toutes les autres situations, la level set ne doit s’annuler qu’une fois le long d’une arête. Il sera alors aisé de détecter les annulations de la level set le long d’une arête pour former les sous éléments d’intégration. Pour se ramener dans ces situations, un réajustement est effectué lorsque l’on se trouve dans l’une des situations représentées sur la figure Fig. 346.
Fig. 345 Annulation multiple de la level set le long d’une arête#
Fig. 346 Réajustement complémentaire des level set#
Dans les deux configurations de la Fig. 346 où l’on ajuste la level set de l’un des nœuds extrémités à zéro, c’est l’extrémité de l’arête qui a la plus petite level set en valeur absolue qui est choisie pour être ajustée. Aussi, lorsque l’ajustement à effectuer est jugé trop important, ce qui révèle une courbure importante de la level set par rapport à la taille de l’élément, un message d’alarme préconisant d’utiliser un maillage plus fin est émis. Le critère choisi est le suivant: le message d’alarme est émis dès lors que:
Si les level sets sont calculées par des fonctions analytiques, le choix des fonctions n’est pas unique pour une même géométrie de fissure. En effet, n’importe quelle fonction à valeurs positives d’un côté et à valeurs négatives de l’autre est valable. Cependant, si on désire que la level set représente la vraie distance signée, alors le choix est unique. Cette notion est importante par la suite lorsque l’on définit les coordonnées dans la base locale au fond de fissure à l’aide des valeurs des level sets. Il faut donc veiller à donner la formule exacte de la distance à la fissure (expression qui n’est pas aisée pour des géométries de fissures complexes).
Un des intérêts de la méthode par projection est qu’elle fournit un champ qui est la vraie distance signée. On pourrait donc aussi envisager de combiner les deux méthodes. La donnée d’une fonction analytique simple permettrait dans un premier temps de déterminer l’iso-zéro; puis de créer un maillage simplifié de cette surface, qui servirait de support à la méthode par projection. L’inconvénient est de créer un maillage 2D virtuel de la fissure qui est disjoint du maillage réel 3D.
Sinon, on pourrait aussi envisager une phase d’orthogonalisation qui transforme une level set quelconque en une vraie distance [bib21].
Remarque:
Si la forme de la level-set représentée, n’est pas assez régulière (par exemple, possède un angle vif ou une bifurcation), la définition de la distance par projection orthogonale sur l’iso-zéro, peut-être ambiguëFig. 347.
Au voisinage de la singularité, la projection sur la surface de la level-set n’est pas possible: le gradient l’iso-zéro de la level-set est discontinu, donc la normale à l’iso-zéro de la level-set n’existe pas. On désigne alors par «zone d’ombre», la région où le calcul de distance par projection orthogonale sur l’iso-zéro, est à proscrire.
Dans la «zone d’ombre», on calcule la distance en revenant à sa définition fondamentale: la distanced’un point\(M\)à l’iso-zéro est la longueur minimale reliant le point\(M\)à un point de l’iso-zéro.
Soit\(C\)le point correspondant au sommet de l’angle droit sur l’iso-zéro.
Pour tout point\(M\)dans la «zone d’ombre»Fig. 347,\(C\)vérifie lalongueur minimale, c’est-à-dire:\(d(M,\text{iso\_zero})=min(\parallel \vec{\mathit{MN}}\parallel ,N\in \text{{iso\_zero}})=\parallel \vec{\mathit{MC}}\parallel\).
Fig. 347 Difficulté de la définition d’une distance par projection orthogonale en cas de singularité de l’iso-zéro#
Par ailleurs, notons que l’interpolation de level-set régularise la singularité. Comme la level-set non régulière est interpolée par des fonctions polynomialesr7.02.12-approximation-level-sets, la level-set discrétisée devient une fonction régulière au sein de l’élément.
Les level sets sont définies pour chaque fissure par l’opérateur DEFI_FISS_XFEM. La liste des fissures est nécessaire pour la création des éléments finis X-FEM (opérateur MODI_MODELE_XFEM). Pendant cette phase, des champs «concaténés» sont créés. Par exemple, on créé un champ de level set normale global à toutes les fissures. Pour chaque nœud (champs aux nœuds) ou chaque maille (champs par élément), on va chercher l’information associée à la fissure la plus proche.
Les champs aux nœuds «concaténés» créés sont:
Level set normale,
Level set tangente,
Statut des nœuds (voir r7.02.12-enrichissement-statut-nœuds),
Base locale au fond de fissure (voir r7.02.12-base-locale-fissure).
Ceux lies au sous-découpage (voir r7.02.12-sous-decoupage),
Les champs par éléments «concaténés» créés sont:
Ceux lies aux structures de données pour le contact (voir document R5.03.54).
Pour ce qui concerne la propagation r7.02.13, on peut définir plusieurs fissures sur le modèle et on peut donner la liste des fissures qui se propagent: à chaque appel de l’opérateur PROPA_FISS, toutes les fissures données se propage.
Restrictions:
Les fissures doivent être suffisamment espacées l’une de l’autre (3 mailles sans fissures doivent au moins les séparer). A fortiori , les fissures ne doivent pas non plus se croiser. Sinon, l’introduction d’enrichissements spéciaux est nécessaire [bib22].
où les \({\phi }_{i,j}\) sont les dérivées des fonctions de forme par rapport à la direction \(j\) .
On détermine ainsi un champ de gradients par élément. Les valeurs sont calculées aux nœuds des éléments, pour chaque élément indépendamment des autres; puis pour calculer la valeur nodale on moyenne sur les valeurs obtenues par éléments aux nœuds.
La base locale au fond de fissure \(\left\lbrace {e}_{1},{e}_{2},{e}_{3}\right\rbrace\) se calcule alors en tout point grâce aux champs de gradients nodaux:
Le fond de fissure est défini par l’intersection des iso-zéros des deux level sets. Pour le calcul des facteurs d’intensité de contraintes (Stress Intensity Factors en anglais), il est pratique de définir des points appartenant au fond de fissure, qui serviront de base pour l’interpolation des SIFs (voir r7.02.12-discretisations).
Les points choisis sont les intersections des faces des éléments avec la courbe \(\text{lsn}=0\cap \text{lst}=0\) . Ces points seront ensuite ordonnés de manière à définir une abscisse curviligne le long du fond de fissure.
La recherche des points du fond de fissure se fait de la manière suivantes. On se restreint aux éléments où la level set normale change de signe et où tous les nœuds sont de statuts «Crack-Tip» (cette notion est définie à la r7.02.12-enrichissement-statut-nœuds).
Soit une face d’un tel élément, si cette face est un quadrangle, on se ramène à deux triangles.
Soit le triangle \(\mathit{ABC}\)
Fig. 349 Face d’un élément intersectée par le fond de fissure#
On cherche le point \(M\) solution du système suivant:
Le point \(M\) est retenu sous réserve qu’il appartienne au triangle \(\mathit{ABC}\) .
En 2D, on utilise le même procédé mais directement sur les mailles et non sur les faces des mailles (toujours des triangles ou bien des quadrangles découpés en triangles).
La recherche des points du fond de fissure s’accompagne d’une stratégie d’élimination des points \(M\) redondants. En effet, si la face n’est pas une face du bord géométrique de la structure, elle est alors commune à deux mailles de la liste des mailles du fond, et un point du fond peut ainsi être détecté (au moins) deux fois. Pour vérifier si le nouveau point \(M\) détecté n’appartient pas déjà à la liste des points \(P\) du fond, on procède à une caractérisation du point \(M\) suivant qu’il se situe sur un nœud sommet de la face, une arête de la face ou bien à l’intérieur de la face. On effectue ensuite une vérification suivant les cas, via une boucle sur les points \(P\) de même statut que \(M\) :
si \(M\) est situé sur un nœud sommet de numéro \({N}^{M}\) et que \({N}^{M}={N}^{P}\) ,
si \(M\) est situé sur une arête caractérisée par les numéros de sommets \({N}_{1}^{M}\) et \({N}_{2}^{M}\) et que \(\left\lbrace \begin{array}{c}{N}_{1}^{M}+{N}_{2}^{M}={N}_{1}^{P}+{N}_{2}^{P}\\ {N}_{1}^{M}.{N}_{2}^{M}={N}_{1}^{P}.{N}_{2}^{P}\end{array}\right.\) ,
si \(M\) est situé à l’intérieur d’une face caractérisée par \(\mathit{NNF}\) sommets et que
alors le point \(M\) n’est pas conservé. Dans le cas contraire, il est ajouté à la liste des points \(P\) du fond. Pour s’affranchir des erreurs numériques qui pourraient engendrer qu’un point \(M\) normalement situé sur un nœud sommet se voit caractérisé comme un point intérieur à la surface, une procédure «fit to vertex» est associée à la stratégie précédente. Ainsi, si \(M\) est détecté comme un point situé sur une arête ou à l’intérieur d’une face, l’algorithme suivant est utilisé:
calcul de la distance de chaque nœud sommet de la face au point \(M\) ,
chacune de ces valeurs est placée dans le vecteur \(\mathit{Dist}\) qui est ensuite trié de la plus petite valeur à la plus grande,
si \(\mathit{Dist}(1)<{10}^{-4}.\mathit{Dist}(2)\) alors \(M\) est replacé sur le nœud sommet. Sinon, il conserve son statut précédent.
Une procédure équivalente est mise en place pour replacer (si besoin est) un point intérieur à la face sur une arête.
Remarque
Même lorsqu’on se restreint aux éléments où la level set normale change de signe et où tous les nœuds sont de type «Crack-Tip», on ne peut pas être sûr que le système soit inversible. On pourrait déjà se limiter aux faces triangulaires où les level sets changent de signe au sens large (0 inclus), mais cela n’éliminerait pas tous les cas de déterminant nul. En effet, lorsque la trace du fond de fissure sur la face triangulaire est une courbe, le système admet une infinité de solution. Ce cas n’est pas détectable à priori, et seul un test sur la non-nullité du déterminant permet de s’affranchir de tels cas. Ainsi, si sur une face le déterminant est nul numériquement, on ne détermine pas de point du fond de fissure sur cette face. S’il en existe (forcément une infinité), les deux points solution sur les bords de la face sont alors déterminés par une autre face de l’élément en question sur laquelle le déterminant est non-nul.
L’orientation du fond de fissure n’est nécessaire qu’en 3D.
Explication de la méthode d’orientation du fond de fissure#
Lors de la recherche des points du fond de fissure par la méthode mentionnée à la r7.02.12-recherche-points-fond-fissure, les points ne sont pas forcément trouvés dans un ordre permettant la définition immédiate d’un chemin ordonné et d’une abscisse curviligne. Pourtant, la définition d’un fond de fissure ordonné et d’une abscisse curviligne le long du fond de fissure est indispensable au calcul de G par la méthode G-thêta.
La méthode d’orientation utilisée est basée sur le fait que deux points consécutifs du fond appartiennent forcément à une même maille 3D.
Pour expliquer cette méthode, on va se baser sur l’exemple d’un fond de fissure défini par une droite traversant une structure parallélépipédique. On a illustré sur la figure suivante la vue du dessus du maillage de cette structure en 3D:
Fig. 350 Vue de dessus d’un maillage d’une structure parallélépipédique#
La numérotation des mailles est arbitraire.
On peut voir sur la Fig. 350 le fond de fissure en rouge constitué de 4 points.
L’ensemble des mailles de cette structure ont leur level set normale qui change de signe. Les 14 mailles représentées sont dites connectées au fond de fissure et possèdent chacune au moins un des quatre points du fond.
La liste ordonnée des points de ce fond est 1-2-3-4.
Pour rechercher les points du fond, on fait une boucle sur les 14 mailles connectées au fond. Cette boucle étant faite dans l’ordre croissant des indices des mailles, on trouve d’abord le point 1 puis les points 3 et 4 (ou 4 puis 3) et enfin le point 2. On obtient donc la liste des points suivante: 1-4-3-2. Il reste à ordonner cette liste.
Étape 1: Recherche des points du fond par maille
Le principe de l’orientation est basé sur le fait qu’une maille connectée au fond contient forcément 1 ou 2 points du fond.
Si par exemple le fond passe par un seul sommet d’une maille, cette dernière ne contiendra qu’un point du fond. En rentrant dans une maille 3D par une face, le fond doit forcément ressortir par un autre point de cette même maille. Celle-ci contient alors deux points du fond.
Remarque :
Si une ou plusieurs mailles possèdent plus de deux points du fond, la procédure échoue et le fond n’est pas orienté. Un élément hexaèdre HEXA8 peut contenir par exemple trois points du fond si celui-ci coupe nettement deux faces opposées et rase une autre face.
Si l’orientation du fond est nécessaire pour la suite des calculs, il faut raffiner le maillage pour n’avoir que des mailles contenant au maximum deux points du fond.
La première étape de cette orientation consiste à associer pour chaque maille connectée au fond, un couple d’indices de points du fond leur appartenant.
Dans l’exemple, la maille 1 ne possède que le point 1, on lui associe le couple \((1,0)\) . La maille 2 a les points 3 et 4, on lui associe le couple \((4,3)\) , etc.
Soit \(\mathit{NMAFON}\) le nombre de mailles connectées au fond.
On crée une liste \(\mathit{LISTPT}\) de taille \(2\times \mathit{NMAFON}\) , contenant les indices des points du fond pour chacune des \(\mathit{NMAFON}\) mailles connectées au fond.
Dans notre exemple, la liste \(\mathit{LISTPT}\) contient 14 couples d’indices et pourrait se décrire comme suit:
Indice de la maille
1
2
3
4
5
6
7
8
…
Couple associé à la maille
\((1,0)\)
\((4,3)\)
\((3,0)\)
\((3,4)\)
\((2,0)\)
\((2,0)\)
\((2,1)\)
\((3,2)\)
…
Remarques :
Les couples\((4,3)\)et\((3,4)\)sont équivalents.
Les couples contenant un 0 ne seront pas considérés. En effet, un point appartiendra toujours a au moins un couple d’indices non nuls. Par exemple le point 1 appartient au couple\((2,1)\).
On détecte les points extrémités en regardant à combien de couples d’indices non nuls appartiennent les points. Le point 3 se trouve dans les deux couples différents d’indices non nuls\((3,4)\)et\((2,3)\) *. Ainsi, on sait automatiquement que le point 3 est situé entre le point 2 et le point 4.* Le point 1 fait seulement partie du couple d’indices non nuls\((2,1)\). Le point 1 est donc obligatoirement une extrémité du fond.
Étape 2: Recherche des extrémités du fond
Comme on vient de le mentionner, les extrémités du ou des fonds de fissure sont repérables car un seul couple d’indices tous non nuls dans \(\mathit{LISTPT}\) contient leur indice. Grâce à cette indication, on crée la liste \(\mathit{PTEXTR}\) comprenant l’ensemble des points extrémités.
Dans l’exemple, le point 1 n’appartient qu’au couple \((2,1)\) et le point 4 qu’au couple \((4,3)\) . Ce sont alors les deux extrémités de notre fond.
Étape 3: Ordonnancement des points du fond
Cette étape correspond à l’orientation du fond à proprement parler.
Pour commencer l’ordonnancement on part d’un point extrémité. Celui-ci est le premier point de la liste \(\mathit{PTEXTR}\) . Si le fond est fermé (fond sans extrémité), le premier point est celui ayant l’indice 1.
Dans notre exemple, le premier point est le point 1. Pour connaître le point suivant, il suffit de chercher le couple possédant un 1 et une deuxième valeur non nulle. Dans notre cas, il n’existe qu’un couple avec le 1 à savoir le \((2,1)\) , donc le deuxième point du fond est le 2. Ensuite, le point 2 ne se trouve que dans le couple \((2,3)\) , donc le point suivant est le 3. On continue cette démarche jusqu’à ce que tous les points soient ordonnés.
On a alors obtenu la liste des points 1-2-3-4.
Remarque
Si la première extrémité trouvée avait été le point 4, on aurait ordonné le fond dans l’autre sens. On aurait obtenu la liste ordonnée des points 4-3-2-1.
En 3D, le fond de fissure est une ligne soit fermée (fissure non débouchante), soit ouverte (fissure débouchante). Dans la plupart des cas, le fond de fissure est une ligne continue, comme celui de la fissure circulaire représentée sur la Fig. 351.
Cependant, il peut arriver que le fond de fissure soit en fait composé de plusieurs morceaux discontinus. C’est le cas par exemple de la fissure circulaire représentée sur la Fig. 352. Dans ce cas, on parle toujours du fond de fissure, comme de l’ensemble des morceaux du fond. On dit que le fond de fissure est un fond multiple. Sur l’exemple de la Fig. 352, le fond de fissure est composé des lignes courbes \((\mathit{BC})\) , \((\mathit{DE})\) , \((\mathit{FG})\) et \((\mathit{HA})\) .
La Fig. 353 illustre la discrétisation d’un fond de fissure multiple (comprenant deux morceaux). Le premier morceau est composé des points 1 à 5 et le second morceau est composé des points 6 à 10.
Dans le cas d’un fond multiple, les deux premières étapes d’ordonnancement des points sont inchangées par rapport au cas d’un fond continu.
Lors de la troisième étape, les fissures ayant un fond multiple sont détectées lorsque l’on atteint une deuxième extrémité, sans avoir ordonné l’ensemble des points du fond. On recherche alors le départ du fond suivant en prenant un point de la liste \(\mathit{PTEXTR}\) .
Dans le cas illustré sur la Fig. 353, \(\mathit{PTEXTR}\) contient les points 1, 5, 6 et 10. On commence par ordonner les points de 1 à 5. Puis, comme on sait qu’il reste cinq points à ordonner, on prend le point 6 de \(\mathit{PTEXTR}\) pour continuer l’orientation du fond.
Remarque
En 2D, chaque point du fond de fissure représente à lui seul un fond.
# le premier point de la liste des points ordonnés:
\(\mathit{TABPT}(1)=\mathit{PTEXTR}(1)\)
# une fois le point extrémité utilisé, on le met à 0, pour éviter de le prendre en compte en cas de recherche d’une autre extrémité si on a un fond multiple:
\(\mathit{PTEXTR}(1)=0\)
Boucle sur les points du fond: \(\mathit{IPT}=1\) à \(\mathit{NFON}-1\) :
Boucle sur les mailles connectées au fond: \(\mathit{IMA}=1\) à \(\mathit{NMAFON}\) :
On cherche le couple d’indices non nuls où se situe le point du fond dernièrement rentré dans \(\mathit{TABPT}\) . Si c’est le cas, le deuxième indice du couple est forcément le point du fond suivant :
Si le deuxième indice du couple associé à \(\mathit{IMA}\) vaut 0, on ignore ce couple: passer à la maille suivante \(\mathit{IMA}+1\)
Si le premier indice du couple correspond à celui recherché:
Si \(\mathit{IPT}>2\) :
Si le couple de \(\mathit{IMA}\) a déjà été trouvé (cas d’un couple en double) :
mettre la deuxième valeur du couple à 0 pour ignorer ce couple à la prochaine itération
aller à la maille suivante \(\mathit{IMA}+1\)
Si on a affaire à un nouveau couple, le deuxième indice est le point du fond suivant :
mettre la deuxième valeur du couple à 0 pour ignorer ce couple à la prochaine itération
Si on n’a pas trouvé de points à associer à \(\mathit{IPT}\) , c’est donc un point extrémité (cas des fonds multiples)
on vérifie que \(\mathit{IPT}\) est bien un point de la liste \(\mathit{PTEXTR}\) puis on le met à 0 pour l’ignorer à la prochaine itération
on recherche un nouveau point extrémité dans la liste \(\mathit{PTEXTR}\) pour débuter le nouveau fond de fissure et on passe à la maille suivante \(\mathit{IMA}+1\)
Dans cette partie, on rappelle les équations du problème général d’une structure fissurée. On considère une fissure \({\Gamma}_{c}\) dans un domaine \(\Omega \in {\Re }^{3}\) délimité par \(\partial \Omega\) de normale extérieure \({n}_{\mathrm{ext}}\) . Les lèvres de la fissure sont notées \({\Gamma}^{1}\) et \({\Gamma}^{2}\) de normales extérieures \({n}^{1}\) et \({n}^{2}\) . Les champs de contraintes et de déplacements sont respectivement notés \(\sigma\) et \(u\) .
Un chargement quasi-statique est imposé sur la structure par l’intermédiaire d’une densité de forces volumiques \(f\) , d’une densité de forces surfaciques \(t\) sur \({\Gamma}_{t}\) et d’une densité de forces surfaciques \(g\) sur les lèvres. Le solide est encastré sur \({\Gamma}_{u}\) .
L’idée principale est d’enrichir la base des fonctions d’interpolation grâce à la partition de l’unité [bib23]. On rappelle l’approximation éléments finis classique:
où les \({a}_{i}\) sont les degrés de liberté de déplacement au nœud i et \({\phi }_{i}\) les fonctions de forme associées au nœud i. \({N}_{n}(x)\) est l’ensemble des nœuds dont le support contient le point \(x\) . On assimile le support d’un nœud i au support des fonctions de forme associées à ce nœud, c’est-à-dire à l’ensemble de points \(x\) tels que \({\phi }_{i}(x)\ne 0\) .
Cette expression est composée de 3 termes. Le 1er terme est le terme classique continu. Les 2ème et 3ème termes sont des termes enrichi. Étant au cœur de la méthode X-FEM, ces termes sont explicités dans les paragraphes suivants.
Enrichissement avec une fonction de sélection de domaine (2ème terme)#
Supposons que l’interface \({\Gamma}_{c}\) partitionne le domaine tel que \(\Omega ={\Omega}_{+}\cup {\Omega}_{-}\) . Si \({\Gamma}_{c}\) est une fissure, on partitionne de la même manière le domaine \(\Omega\) en étendant virtuellement \({\Gamma}_{c}\) .
Afin de représenter le saut de déplacement à travers \({\Gamma}_{c}\) , on introduit la fonction de sélection de domaine ou fonction caractéristique de domaine \({H}_{j}\left(x\right)\)[bib76] définie par:
En se servant de la level set normale, la quantité \({H}_{j}\left(\text{lsn}\left(x\right)\right)\) vaut \(0\) si le point \(x\) et le nœud \({x}_{j}\) se trouvent du même coté de la fissure et \(\pm 2\) sinon. Le coefficient «2» est introduit pour avoir une écriture plus simple du saut moyen de déplacement le long de l’interface.
Les \({b}_{j}\) sont les degrés de liberté enrichis. \(K\) est l’ensemble des nœuds dont le support est entièrement coupé par la fissure (nœuds représentés par un rond sur la Fig. 355).
Remarque
Par abus de langage, on appellera aussi les fonctions de sélection de domaine, fonctions «Heaviside». Mais il faudra se référer à la définition ci-dessus, pour représenter l’approximation de la discontinuité du champ de déplacement.
Fig. 355 À gauche, les nœuds «ronds» sont enrichis par la fonction Heaviside et les nœuds «carrés» par les fonctions singulières (enrichissement topologique). À droite, les nœuds «carrés» sont enrichis par les fonctions singulières (enrichissement géométrique).#
Enrichissement avec une fonction de sélection de domaine pour un branchement de fissures#
Supposons un branchement de deux discontinuités, qui partitionne l’espace en trois domaines tel que \(\Omega ={\Omega}_{1}\cup {\Omega}_{2}\cup {\Omega}_{3}\) . S’il s’agit d’un branchement de deux fissures, on partitionne de la même manière le domaine \(\Omega\) en étendant virtuellement les trois branches de fissures.
Fig. 356 Partitionnement du domaine par une jonction de fissure#
Pour un branchement simple, la généralisation des fonctions de sélection de domaines conduit à l’écriture suivante:
Toutefois, la description topologique des trois domaines connexes (Ω1, Ω2 et Ω3) n’est pas triviale, compte-tenu de la seule information fournie par les level-sets. Les levet-sets permettent de représenter un changement de signe scalaire à travers une discontinuité. Alors que nous souhaitons représenter le partitionnement de l’espace au voisinage de la discontinuité.
D’un point de vue « élémentaire », nous constatons que le champ de signe de level-sets lu sous forme vectorielle, colle approximativement au partitionnement que nous voulons réaliser. Ce partitionnement, grâce à la vectorisation des level-sets scalaires, permet alors de réutiliser avantageusement les structures de données du Code_Aster. Toutefois, le seul partitionnement élémentaire est insuffisant, pour représenter le « ddl domaine » (χΩϕi).
Remarque
En fait, le « ddl domaine » correspond à la collection des partitions du domaine Ω sur l’ensemble du support du nœud i. Sur le support de chaque nœud, il est donc indispensable de concaténer les partitions élémentaires correspondant à un « ddl domaine » donné. Cette opération n’est pas triviale puisque les champs de signe évoluent d’un élément à l’autre (comme les champs de signes scalaires ne sont pas prolongés sur tous les éléments, pour localiser sur quelques éléments d’intérêt la définition d’une fissure, appelés aussi « bande proche » d’enrichissement). Un algorithme dédié à la concaténation du partitionnement élémentaire, a été développé pour résoudre cette difficulté fondamentale.
Fig. 357 Partitionnement d’un élément multi-fissuré à l’aide du champ de signe « vectorisé ». Pour chaque domaine, la première composante scalaire du vecteur signe, correspond à la caractérisation de la discontinuité de la fissure n°1 et la deuxième composante scalaire du vecteur signe, correspond à la discontinuité de la fissure n°2.#
Construction des fonctions de signe ou jonctions:
Rappelons la construction des fonctions de signe Heaviside ou fonctions de type jonction [bib70]. Il s’agit d’une fonction Heaviside qui est «tronquée» au niveau du branchement. Les nœuds enrichis sont représentés sur la Fig. 358.
Fig. 358 Enrichissement pour la deuxième fissure. Les nœuds ronds sont enrichis par la fonction jonction, qui vaut +1, -1, 0.#
La valeur de cette fonction dépend des levels set normales des 2 fissures. On considère le signe de la level set normale de la fissure 1, du côté où la fissure 2 est définie \(\mathrm{sign}({\mathrm{lsn}}_{1}({\mathrm{fiss}}_{2}))\) (dans la pratique, on regarde le signe d’un point appartenant au domaine de \({\mathrm{lsn}}_{1}\) dans lequel se trouve la fissure 2, cf. opérande JONCTION de la doc u4.82.08). On a alors la fonction d’enrichissement jonction pour la fissure 2 qui s’écrit :
Si on généralise l’approche à une fissure \(N\) qui se branche sur les fissures d’un ensemble \(K\) , on définit alors l’ensemble des fissures \(P\) , qui contient à la fois toutes les fissures de \(K\) et toutes les fissures sur lesquelles se branchent éventuellement les fissures de \(K\) . La fonction jonction s’écrit alors de manière générale :
Un exemple de configuration est présenté la Fig. 359. Sur cet exemple on construit un arbre de connectivité des fissures. On déduit de cet arbre que pour \(N=3\) , on a \(K=2\) et \(P=[1,2]\) . De même pour \(N=5\) , on déduit \(K=[3,4]\) et \(P=[1,2,3,4]\) . Ainsi on a \({J}_{5}(x)=H({\mathrm{lsn}}_{5}(x))\) sur le domaine hachuré de la figure, et zéro ailleurs.
Fig. 359 Réseau de fissures à gauche, arbre de hiérarchie à droite, la zone hachurée correspond au domaine où la fonction d’enrichissement de la fissure 5 n’est pas nulle.#
Exploitation des fonctions de signe pour l’assemblage des ddls domaines:
Pour construire les ddls domaines définis au paragraphe précédent, nous proposons de réutiliser l’information du champ de signes Heaviside pour définir un partitionnement élémentaire. L’algorithme étant beaucoup trop complexe à expliquer en termes de structures de données, nous donnons une traduction graphique de l’algorithme. Voici par exemple les champs de signes Heaviside concaténés qu’on souhaite exploiter pour construire le partitionnement élémentaire :
Fissure simple
Jonction simple
Jonction multiple
Le cas d’une mono-fissure étant assez facile, nous allons nous attarder en détail sur cas d’une jonction simple. Dans l’explication graphique suivante, nous proposons de changer de perspective par rapport à la description élémentaire suggérée ci-dessus. L’algorithme sera décrit du point de vue du support du nœud, qui est plus adapté pour décrire les fonctions d’enrichissement nodales.
Du point de vue du nœud, l’information sur les champs de signes est plus riche que le cas idyllique décrit ci-dessus. En effet, le nœud voit l’information sur les fissures proches, c’est -à-dire, situées dans la première et la deuxième bande d’éléments attenants au support du nœud. Dans le support du nœud, il existe donc un statut de fissures (cf. [D4.10.02]) pour discriminer les fissures intersectant le support du nœud et les fissures proches n’intersectant pas le support du nœud.
Cette difficulté rend non triviale l’association entre les champs de signes et les partitions de domaines. Cette association est nécessairement locale et doit prendre en compte la pollution des champs de signe liée aux bandes proches citées précédemment. Sur l’exemple donné dans le tableau ci-dessous, il faut associer les domaines au champ de signe de la manière suivante:
on associe le domaine \({\Omega}_{1}\) ↔ au champ de signe \(\lbrace \text{+1},0,\text{*}\rbrace\)
on associe le domaine \({\Omega}_{2}\) ↔ au champ de signe \(\lbrace \text{-1},\text{-1},\text{*}\rbrace\)
on associe le domaine \({\Omega}_{3}\) ↔ au champ de signe \(\lbrace \text{-1},\text{+1},\text{*}\rbrace\)
Partitionnement en domaine
Champ de signe associé
Pour compresser l’information des vecteurs signes, on préfère transformer les vecteurs signes à l’aide de la fonction de codage réversible suivante:
ifiss est le numéro dans le vecteur champ de signe,
P désigne le point courant (un nœud ou un point de gauss),
nfiss est la longueur du champ de signe.
Fig. 360 Les vecteurs signes ont des longueurs variables, ce qui conduit à des identifiants différents pour un même domaine. Ces identifiants constituent alors un deuxième niveau de sous-partitionnement.#
Compte-tenu de la variabilité de la longueur nfiss des vecteurs signes, un domaine peut avoir un code différent d’un élément à un autre, comme illustré la Fig. 360. Par exemple, le domaine \({\Omega}_{1}\) possède deux identifiants (5 et 6). Ces identifiants peuvent être vus comme deux sous-partitions du domaine \({\Omega}_{1}\) .
Ce sous-partitionnement n’est pas gênant tant que l’information sur ce sous-partitionnement est regroupée au niveau du nœud, pour l’assemblage des «ddls» domaine. Pour identifier ces deux ddls domaines dans chaque élément, il suffit alors de concaténer l’information sur les sous-partitions :
Si une sous-partition du domaine complémentaire est trouvée, l’information est stockée dans l’emplacement dédié (structure de données par nœud et par élément),
Sinon, il n’y a pas de sous-partition pour cet élément.
Le résultat d’une telle boucle élémentaire est résumé dans le tableau ci-dessous. Dans chaque élément on recherche une sous-partition de chaque domaine complémentaire au nœud bleu. Si aucune sous-partition n’est trouvée dans l’élément, on marque dans l’emplacement dédié une croix ‘x’ (ou -1 par exemple).
Pour le nœud en bleu appartenant au domaine \({\Omega}_{1}\) , l’évaluation des fonctions caractéristiques \({\chi }_{{\Omega}_{2}}\) ou \({\mathrm{\chi }}_{{\Omega}_{3}}\) est évidente compte-tenu de l’information sur le partitionnement du tableau ci-dessus. Dans chaque élément, la première composante (du champ concaténé) renseigne sur l’identifiant du premier ddl domaine \({\mathrm{\chi }}_{{\Omega}_{2}}\) , la deuxième composante (du champ concaténé) renseigne sur l’identifiant du deuxième ddl domaine \({\mathrm{\chi }}_{{\Omega}_{3}}\) :
Si l’identifiant du domaine auquel appartient le point de Gauss correspond à la composante \(k\in [1,2]\) stockée au nœud, alors la fonction caractéristique associée à la composante \(k\) prend la valeur «+2». Précisons à nouveau que la composante \(k=1\) est associée au ddl domaine \({\mathrm{\chi }}_{{\Omega}_{2}}\) et la composante \(k=2\) est associée au ddl domaine \({\mathrm{\chi }}_{{\Omega}_{3}}\) . Il est important de souligner que cette position ne change pas d’un élément à un autre, pour assurer la contiguïté du partitionnement. Par exemple, soit la collection des identifiants {0,1} correspondante au domaine \({\Omega}_{2}\) . C’est exactement l’information qui est stockée en première position dans la structure de donnée au nœud «bleu», dans les éléments intersectant le domaine \({\Omega}_{2}\) .
Sinon, la fonction caractéristique prend la valeur zéro si l’identifiant au point de gauss ne correspond pas à la composante au nœud dans l’élément considéré.
Remarques:
On peut remarquer que le ddl domaine\({\chi }_{{\Omega}_{1}}\)n’est pas traiter dans l’exemple ci-dessus. En effet, pour représenter les deux discontinuités introduites par le branchement de deux fissures il est nécessaire d’enrichir avec deux fonctions discontinues les nœuds dont le support est intersecté par la double discontinuité. On choisit alors d’enrichir avec les ddls associés aux domaines complémentaires au domaine auquel appartient le nœud . Ce choix est bien adapté, compte-tenu des problèmes de conditionnement détaillés à lar7.02.12-conditionnement-enrichissement.
Dans code_aster, on préfère stocker l’information du signe Heaviside et celle de l’identifiant de domaine au point de Gauss dans les sous-éléments d’intégration.
Enrichissement avec les fonctions singulières (3ème terme)#
Afin de représenter la singularité en fond de fissure, on enrichit l’approximation avec des fonctions basées sur les développements asymptotiques du champ de déplacement en mécanique de la rupture élastique linéaire [bib25]. Ces expressions ont été déterminées pour une fissure plane en milieu infini.
L’hypothèse de contraintes planes ne peut être retenue car aucune plaque fissurée n’est en situation de contraintes planes au voisinage de la singularité, lorsque l’on se place à une distance finie de la peau de la coque.
La base permettant de décrire ces champs comporte 4 fonctions:
Dans la pratique, on utilise plutôt la fonction informatique \(\text{atan2}(\text{lsn},\text{lst})\) qui renvoie la valeur principale de l’argument du nombre complexe \((\text{lst},\text{lsn})\) exprimé en radians dans l’intervalle \(\left[-\pi ,\pi \right]\) .
Pour des points se situant exactement sur la lèvre inférieure \((\text{lsn}=0)\) , la fonction \(\text{atan2}(0,\text{lst})\) donne un angle valant \(\pi\) . En théorie, \(\text{atan2}(-0,\text{lst})\) permet d’obtenir \(-\pi\) comme attendu, mais ce n’est pas numériquement toujours le cas. Pour palier cet inconvénient, on utilise plutôt l’expression suivante pour l’angle \(\theta\) :
où \(H(\text{lsn})\) est la valeur de la fonction Heaviside. Ainsi, lorsqu’on se trouve sur la lèvre inférieure, la valeur \(-\pi\) est bien atteinte.
Fig. 361 Coordonnées polaires dans la base locale#
On note que seule la 1ère fonction de la base est discontinue au travers de la fissure. Les autres fonctions ne sont ajoutées que pour améliorer la précision. Ces fonctions sont les solutions de Westergaard, solutions asymptotiques analytiques d’un problème de rupture élastique en 2D. Cette base est bien adaptée aux cas 3D [bib16], [bib18], au moins pour les fissures dont le fond est assez régulier. Ces fonctions sont dites «singulières» car leurs dérivées sont singulières en \(r=0\) .
Fig. 362 Fonctions d’enrichissement en fond de fissure#
Les \({c}_{k}^{\alpha}\) sont les degrés de liberté enrichis. \(L\) est l’ensemble des nœuds dont le support est partiellement coupé par le fond de fissure (nœuds représentés par un carré sur la Fig. 355). Cela signifie qu’une seule couche d’éléments est enrichie autour du fond de fissure. Cet enrichissement est appelé «topologique».
Dès les premiers papiers sur X-FEM [bib24], il est notifié qu’un critère géométrique \({r}_{\max}\) peut être défini pour déterminer les nœuds enrichis par les fonctions singulières (voir Fig. 355):
Les premières études de convergence ont été effectuées en 2000 dans le cadre de la GFEM [bib26], avec prise en compte de plusieurs couches d’éléments enrichis en fond de fissure.
Quand on étudie la convergence, on s’intéresse à l’évolution de l’erreur par rapport au degré de raffinement du maillage. Généralement, on désigne par \(h\) une longueur caractéristique des éléments du maillage, et on cherche à déterminer le paramètre \(\alpha\) appelé taux (ou vitesse ou ordre) de convergence, tel que l’erreur relative \({e}_{\text{rel}}\) s’écrive de la forme:
le paramètre \(\alpha\) apparaît comme la pente de la droite \(\log{e}_{\text{rel}}\) en fonction de \(\log{h}\) lorsque \(h\) tend vers 0.
Stazi et al.[bib27] étudie la convergence de l’erreur en énergie pour une plaque infinie avec une fissure droite, en mode I, pour des formulations linéaires et quadratiques. Il remarque que le quadratique améliore l’erreur, mais pas le taux de convergence. Béchet et al.[bib28] confirme cette observation et montre qu’une zone d’enrichissement fixe permet de retrouver un taux de convergence presque optimal.
Parallèlement, Laborde et al.[bib29] approfondit la question, et teste les taux de convergence pour les formulations polynomiales d’ordre supérieur. De plus, il apporte des améliorations afin de retrouver un taux optimal, voire une supra convergence. Le Tableau 77 rassemble les résultats obtenus par Laborde pour différentes variantes d’X-FEM, et ceci pour des approximations polynomiales de degré \(k=1,2,3\) .
Tableau 77 Ordre de convergence des différentes variantes d’X-FEM#
FEM
X-FEM
X-FEM (f. a.)
X-FEM (d. g.)
X-FEM (p. m.)
P1
0.5
0.5
0.9
0.5
1.1
P2
0.5
0.5
1.8
1.5
2.2
P3
0.5
0.5
2.6
2.6
3.3
La première colonne correspond aux ordres de convergence de la méthode des éléments finis classique pour un problème de fissuration. Compte tenu de la singularité, la vitesse de convergence est en \(\sqrt{h}\) quel que soit le degré \(k\) . Les simulations 2D réalisées sur un problème test: fissure rectiligne sur un carré en mode I d’ouverture montre que X-FEM n’améliore pas le taux de convergence. Ceci peut s’expliquer par le fait que l’enrichissement topologique ne concerne qu’une seule couche d’éléments en fond de fissure. La zone d’influence de cet enrichissement est donc fortement lié à \(h\) . Ainsi, lorsque \(h\) tend vers 0, la taille de la zone d’influence de l’enrichissement tend aussi vers 0. L’idée qui semble naturelle est alors de ne plus limiter la zone d’enrichissement à une seule couche d’éléments, mais de l’élargir à une zone de taille fixe, indépendante du raffinement du maillage. La 3ème colonne du Tableau 77 présente les résultats obtenus avec cette méthode dite X-FEM f. a. (pour Fixed enrichment Area). On retrouve presque les taux de convergence espérés (\(\alpha =k\) pour une approximation \(\mathit{Pk}\) ). Cependant, le conditionnement est dégradé par rapport à la méthode X-FEM habituelle. Afin de retrouver un conditionnement acceptable, Laborde propose de rassembler les degrés de libertés enrichis par les fonctions singulières. En clair, au lieu d’avoir des ddls enrichis différents pour chaque nœud enrichi, ils sont globalisés afin d’en avoir un seul par fonction singulière et par bout de fissure (mais en 3D, on voit pas bien comment ça marche). Avec cette arrangement, le conditionnement est grandement amélioré, mais les taux de convergence sont plus faibles (X-FEM d. g.). Le problème provient des éléments de transition entre la zone enrichie et la zone non-enrichie. D’après Laborde, le phénomène s’explique par l’effet de la partition de l’unité qui ne peut pas être utilisée sur ces éléments partiellement enrichis. Pour pallier ce défaut, une ultime version est proposée: X-FEM p.m. (pour Pointwise Matching). Les déplacement des nœuds sur la frontière entre la zone enrichie et la zone non-enrichie sont imposés égaux. Grâce à ce recollement, on obtient des taux de convergence espérés (voire une légère super-convergence).
Remarque
Pour les approximation polynomiales de degré* \(k=2,3\) , Laborde et al. [bib29] comme beaucoup d’autres utilise des fonctions de forme de degré* \(k\) pour les termes classiques et enrichis par les fonctions de discontinuité (2eme terme), de manière à ce que le saut soit de degré \(k\) . Par contre, pour les termes enrichis par les fonctions singulières (3eme terme), il suffit d’utiliser les fonctions de forme linéaires pour capter la singularité en fond de fissure et de pas trop détériorer le conditionnement des matrices.
Actuellement, dans code_aster , seule une approximation de degré 1 est possible, et le fait de renseigner ou non une valeur de rayon pour la zone d’enrichissement (mot-clé RAYON_ENRI) permet de se placer respectivement dans le cadre X-FEM (f. a.) ou X-FEM classique.
Pour savoir si un nœud est enrichi par la fonction Heaviside (nœud de type «Heaviside) ou par les fonctions singulières (nœud de type «crack-tip»), on calcule le min et le max de la level set normale sur tous les nœuds appartenant au support du nœud (la notion de «support» est définie à la r7.02.12-enrichissement-approximation-deplacement), et on calcule le min et max de la level set tangente sur tous les points appartenant au support du nœud considéré où la level set normale s’annule.
Des idées similaires apparaissent dans [bib17], mais il semblerait que certains cas de figures n’aient pas été pris en compte. Ces expressions sont l’aboutissement des premiers efforts [bib30] qui visaient à déterminer les types d’enrichissement uniquement à l’aide des level sets. Le lecteur y trouvera les figures correspondantes explicatives.
Dans le cas d’un enrichissement géométrique en fond de fissure, le critère de sélection des nœuds est le suivant:
Concernant le choix de la valeur du rayon d’enrichissement, rien n’est clairement indiqué dans la littérature. Il semblerait cependant qu’un rayon valant entre \(1/5\) et \(1/10\) de la longueur de la fissure soit un choix pertinent.
Des récentes études ont montré que l’enrichissement géométrique dégradant fortement le conditionnement de la matrice de rigidité, il fallait le limiter dans une zone restreinte autour du fond de fissure, dans l’attente de traitement permettant d’améliorer le conditionnement. On propose une alternative, qui est à mi-chemin entre l’enrichissement topologique et géométrique: un enrichissement sur \(n\) couches [bib31]. Dans ce cas, on calcule un \({r}_{max}\) (unique pour chaque fissure) d’après la donnée utilisateur du nombre de couches, puis on applique la formule précédente.
Algorithme de choix de l’enrichissement des nœuds:
Soit \(\mathit{MAFIS}\) l’ensemble des mailles sur lesquelles la level set normale s’annule
boucle sur tous les nœuds \(P\) du maillage
initialisation des max et min des level sets
boucle sur les mailles de \(\mathit{MAFIS}\) contenant le nœud \(P\)
boucle sur les arêtes de la maille
soient \(A\) et \(B\) les deux extrémités du segment
si \(\text{lsn}\left(A\right)=0\) alors
actualisation si nécessaire de \(\text{maxlst}\) et \(\text{minlst}\) avec \(\text{lst}(A)\)
fin si
si \(\text{lsn}(B)=0\) alors
actualisation si nécessaire de \(\text{maxlst}\) et \(\text{minlst}\) avec \(\text{lst}(B)\)
On détermine la coordonnée dans l’espace de référence \(\xi\) du point d’intersection \(C\) avec la courbe \(\text{lsn}=0\) sur l’arête. Sa \(\mathit{lst}\) est alors donnée par:
actualisation si nécessaire de \(\text{maxlst}\) et \(\text{minlst}\) avec \(\text{lst}(C)\) .
fin si
fin boucle
boucle sur les nœuds sommets de la maille
actualisation si nécessaire de \(\text{maxlsn}\) et \(\text{minlsn}\)
fin boucle
fin boucle
si \((\text{minlsn}.\text{maxlsn}<0)\) et \((\text{maxlst}\le 0)\) alors \(P\in K\)
cas d’enrichissement topologique
si \((\text{minlsn}.\text{maxlsn}\le 0)\) et \((\text{minlst}.\text{maxlst}\le 0)\) alors \(P\in L\)
cas d’enrichissement géométrique:
si \(\sqrt{\text{lsn}{(P)}^{2}+\text{lst}{(P)}^{2}}\le {r}_{\max}\) alors \(P\in L\)
fin boucle
Pour obtenir l’équation (4905), ainsi que la valeur de la level set tangente au point \(C\) , on détermine au préalable l’abscisse curviligne \(s\) telle que
\(C=A+s(B-A)\)
grâce au fait que la level set normale s’annule en \(C\) , soit
Pour obtenir, dans le cas quadratique, l’équation (4286), ainsi que la valeur de la level set tangente au point \(C\) , on détermine au préalable la coordonnée dans l’espace de référence \(\xi\) du point \(C\)
La \(\mathit{lsn}\) est interpolée quadratiquement le long de l’arête. Il s’agit donc de résoudre un polynôme du second degré:
L’équation \(\mathit{lsn}=0\) possède une unique solution sur l’arête car \(\text{lsn}(A)\text{lsn}(B)<0\)
Une fois obtenue la coordonnée \(\xi\) du point \(C\) , on interpole les \(\mathit{lst}\) des points \(A\) , \(B\) et \(M\) pour obtenir la \(\mathit{lst}\) du point \(C\) :
Dans Code_Aster , il faut définir des types d’éléments finis précis, et pour ne pas multiplier le nombre des possibilités, le choix a été fait de définir 3 types d’éléments finis X-FEM: les éléments «Heaviside», les éléments «crack-tip» et les éléments mixte «Heaviside et crack-tip».
Si la maille possède au moins un nœud de type «Heaviside», alors c’est une maille «Heaviside».
Si la maille possède au moins un nœud de type «crack-tip», alors c’est une maille «crack-tip».
Si la maille possède au moins un nœud «Heaviside» et au moins un nœud «crack-tip», ou si la maille contient au moins un nœud «Heaviside et crack-tip», alors c’est une maille «Heaviside et crack-tip».
Notons \(\mathit{GRMAEN1}\) les mailles «Heaviside», \(\mathit{GRMAEN2}\) les mailles «crack-tip» et \(\mathit{GRMAEN3}\) les mailles «Heaviside et crack-tip».
Soit la maille \(i\) , et \(\text{Ni}\) l’ensemble des nœuds \(j\) de la maille \(i\) .
On remarque que tous les nœuds d’un élément seront affectés des mêmes caractéristiques et du même enrichissement, or ce n’est pas forcément ce qui est voulu. Il faut donc annuler les degrés de liberté enrichis «en trop».
Annulation des degrés de liberté enrichis « en trop »#
On a vu au paragraphe précédent qu’une maille de type «Heaviside» ne peut comporter par exemple qu’un seul nœud de type «Heaviside», les autres nœuds de la maille étant des nœuds classiques ne nécessitant aucun enrichissement. Or ces nœuds vont être affectés des degrés de liberté de la maille «Heaviside», donc de degrés de liberté enrichis. Par conséquent, il faut procéder à une annulation de ces degrés de liberté enrichis à tort. L’annulation permet dans les faits de passer continûment d’une zone enrichie à une zone non enrichie et permet de faire cohabiter deux types éléments (au sens d’ils partagent une frontière commune), l’un enrichi, l’autre non enrichi. La variable non enrichie est la même sur la frontière commune et les degrés de liberté correspondant à l’enrichissement sont mis à zéro pour l’élément qui est enrichi sur cette même frontière (mais pas ailleurs dans ce même élément). Cette façon de faire évite d’avoir à résoudre la question des «blending elements» dont on peut trouver un traitement dans [bib73].
Plusieurs cas se présentent:
degrés de liberté Heaviside à annuler aux nœuds classiques d’une maille Heaviside ou mixte
degrés de liberté crack-tip à annuler aux nœuds classiques d’une maille crack-tip ou mixte
degrés de liberté Heaviside à annuler aux nœuds crack-tip d’une maille mixte
degrés de liberté crack-tip à annuler aux nœuds Heaviside d’une maille mixte.
La technique d’annulation de ces degrés de liberté est expliquée en détail dans r5.03.54 §4.4.
Le conditionnement, noté \({10}^{\delta}\) correspond au rapport entre la plus grande et la plus petite valeur propre d’un système à inverser. Pour un calcul en double précision avec une erreur numérique en \({10}^{-15}\) , l’erreur relative obtenue sur le calcul est de l’ordre de \({10}^{-15+\delta }\) . On doit donc vérifier la condition \(\delta <9\) pour garantir une précision numérique de l’ordre de \({10}^{-6}\) .
L’enrichissement géométrique dégrade fortement le conditionnement de la matrice de rigidité [bib29], [bib28]. Béchet et al. [bib28] proposent une technique d’orthogonalisation des degrés de liberté lors du calcul des matrices de rigidité élémentaires afin d’améliorer le conditionnement de la matrice assemblée. Laborde et al. [bib29] expliquent que le mauvais conditionnement est dû au fait que la base d’enrichissement choisie ne forme pas une famille libre localement. Ils proposent donc de mettre un seul degré de liberté pour ces fonctions sur toute la zone d’enrichissement et de raccorder les déplacements à la limite entre zones enrichies et non enrichies afin de retrouver des taux de convergence optimaux. Le problème de conditionnement est d’ailleurs tel qu’avec des éléments quadratiques il devient impossible d’obtenir des résultats, sans mettre en place une des techniques [bib29], [bib28]. En effet, pour ces éléments, le mauvais conditionnement est dû non seulement à la partie singulière de l’enrichissement, mais aussi à l’enrichissement Heaviside, lorsque la fissure passe très proche d’un nœud. On présente l’évolution du nombre de conditionnement à mesure que l’interface se rapproche des nœuds du maillage sur la Fig. 364, pour des éléments linéaires et quadratiques respectivement. Les valeurs de cette figure sont très approximatives. D’une part on ne prend pas en compte les éléments du voisinage. D’autre part ces valeurs sont obtenues de façon grossière. Dans le cas du cube coupé dans un coin par exemple (à droite sur la figure), on considère la rigidité du petit cube du coin de côté ́ \({10}^{-\gamma }\) au lieu du petit tétraèdre du coin. Enfin le conditionnement d’un problème global est généralement plus grand que le conditionnement local et augmente lorsque l’on raffine le maillage. En prenant en compte toutes ces considérations, on obtient dans la pratique des conditionnements d’un ordre 2 à 4 fois supérieurs à ceux de la Fig. 364.
Fig. 364 La distance au nœud sommet le plus proche normalisée par la longueur du côté valant \({10}^{-\gamma }\) on montre la dépendance de \(\delta\) nombre de conditionnement, par rapport au paramètre \(\gamma\) pour les éléments linéaires en partie supérieure et quadratiques en partie inférieure.#
Les critères estimatifs pour améliorer le conditionnement#
La technique de réajustement de la level set évoquée à la r7.02.12-reajustement-level-sets permet de se prémunir de ce mauvais conditionnement. En notant \({10}^{-\gamma }\) la distance du point d’intersection de la level set au nœud sommet le plus proche normalisée par la longueur du côté, le réajustement à 1% de la longueur d’arête faite à la r7.02.12-reajustement-level-sets correspond à \(\gamma =2\) . Ce réajustement doit agir suffisamment vite pour que le conditionnement ne soit pas trop détérioré, mais pas pour des valeurs de \(\gamma\) trop faibles de façon à ne pas perturber le système en déplaçant de façon non réaliste la surface de fissuration. Pour des éléments hexaèdres quadratiques, s’il faut que \({10}^{-15+9\gamma +2}\) soit de l’ordre de \({10}^{-5}\) , on obtient \(9\gamma +2=10\) soit un réajustement à 13% de la longueur d’arête. Il n’est donc pas possible dans ce cas là d’activer raisonnablement le réajustement des level sets au sommet, afin que le conditionnement ne soit pas détérioré.
Dans ces conditions, on met en place un critère permettant de détecter si un degré de liberté Heaviside pose problème. Si le critère est vérifié, le degré de liberté est éliminé par mise à zéro, comme à la r7.02.12-annulation-degres-liberte.
Le principe est le même que le critère volumique proposé par Daux [bib22] pour les jonctions. L’idée est de regarder pour chaque nœud dont le support est coupé par la level set, le rapport des tailles des zones de part et d’autres de la level set sur ce support (zones affectées d’une valeur de Heaviside valant ±1). Si:
les degrés de liberté Heaviside du nœud concerné sont mis à zéro dans toutes les directions.
Fig. 365 Illustration de la mise en place du critère volumique. Le support du nœud coupé par la level set est illustré. On compare les volumes gris et blanc par rapport au volume total du support du nœud (réunion des volumes gris et blanc) dans le cas d’une fissure (gauche) et d’une jonction (droite).#
Si ce critère est pertinent pour les éléments linéaires (triangles, tétraèdres) avec une valeur de \(\alpha\) de 4, il n’est pas satisfaisant pour les éléments multi-linéaires (quadrangles, pyramides, pentaèdres, héxaèdres) et quadratiques. En effet la valeur de \(\alpha\) de 4 conduit à éliminer des degrés de liberté qui ne devraient pas l’être, ce qui perturbe la solution, alors que des valeurs plus élevées de \(\alpha\) dégradent le conditionnement. Afin de prendre en compte ces éléments, on utilise un critère de rigidité, qui se base sur une comparaison des rigidités du support et non plus simplement des volumes. Pour chaque nœud dont le support est coupé par la level set, on regarde pour ce support le rapport des rigidités des zones de part et d’autres de la level set (zones affectées d’une valeur de Heaviside valant ±1). Si en un nœud \(n\) où \(\mathit{se}\) sont les sous-éléments de son support, on a:
où \({\phi }_{n,X}\) est la dérivée de la fonction de forme au nœud \(n\) dans la direction globale \(X\) , alors les degrés de liberté Heaviside du nœud \(n\) sont mis à zéro dans toutes les directions. On remarquera que le comportement n’est pas présent dans le critère de rigidité de (4582), le critère ayant été normalisé. Cela permet de discriminer à l’avance les degrés de liberté à éliminer, sans connaissance du problème à résoudre (non-linéarités, plasticité, contact, etc). Ce critère, très proche d’un critère de conditionnement, nous amène à choisir des valeurs de \(\delta\) comprises entre 8 et 10. Dans la pratique, nous prenons une valeur de \(\delta\) de 9.
Le critère décrit dans (4582), quantifie en ordre de grandeur, le compromis entre qualité de la solution et conditionnement. Comme ce critère n’est valable qu’en ordre de grandeur, il n’est pas forcément pertinent, de calculer exactement toutes les intégrales. Dans la programmation Aster, on approxime donc le calcul intégral sur les sous-éléments [bib72]:
où \({G}^{\mathrm{se}}\) désigne le barycentre du sous-élément, \({V}^{\mathrm{se}}\) le volume du sous-élément.
Au lieu de calculer l’intégrale sur le sous-élément, on évalue les dérivées de fonctions de forme en un point. Cette stratégie permet de capter en ordre de grandeur le critère d’élimination de (4582), avec un faible coût de calcul.
Cette approximation s’avère robuste et peu coûteuse avec des éléments linéaires.
En revanche, avec des éléments quadratiques, les dérivées des fonctions de forme admettent des racines. Le calcul en un point est risqué: le barycentre peut-être une racine des dérivées de fonctions de forme, où l’estimation sera nulle, mais pas l’intégrale sur le sous-élément. Cette mauvaise estimation peut engendrer aléatoirement des éliminations, ce qui perturbe la solution et n’est pas acceptable.
On renforce donc l’estimation pour les éléments quadratiques, en rajoutant d’autres points d’évaluation, en plus du barycentre:
Pour que l’évaluation soit pertinente, les points sont suffisamment bien distribués: on considère les points équidistants entre les nœuds sommets et le barycentre (par exemple, 3 points pour un tri6).
Toutefois les critiques suivantes, peuvent être faites par rapport à ces critères estimatifs:
d’une part, ils ne garantissent pas le contrôle du conditionnement juste avant l’élimination. Tant que l’élimination n’est pas effectuée, le conditionnement augmente en \({10}^{n\times \gamma }\) (Fig. 364), ce qui peut conduire à l’arrêt du solveur direct, de manière arbitraire. Ce phénomène est exacerbé avec les éléments quadratiques, voir Fig. 364.
d’autre part, une analyse de ces différents critères est faite dans [bib72]. Elle montre que les éliminations conduisent à une absence de convergence de l’erreur en énergie sur un exemple simple de compression homogène d’un cube traversé par une interface inclinée. La solution proposée consiste à remplacer l’élimination des degrés de liberté Heaviside par une orthogonalisation des matrices de rigidité locales, idée qui vient du pré-conditionneur X-FEM de [bib28]. Afin de lever cette dernière difficulté il est nécessaire d’estimer des matrices de rigidité locales au moins en en triple précision [bib72].
Pour adresser, la première critique, c’est-à-dire endiguer l’augmentation exponentielle du conditionnement, nous avons mis en place dans le Code_Aster, le pré-conditionneur proposé par Béchet et al. [bib28].
Il s’agit d’une procédure automatique et algébrique, pour «orthogonaliser» les degrés de liberté associés à un nœud XFEM. En effet, les fonctions d’enrichissement XFEM s’appuient sur les fonctions nodales \({\Phi}_{i}\) pour décrire, soit la discontinuité du déplacement par les fonctions saut \(H{\Phi}_{i}\) , soit le fond de fissure par les fonctions singulières \({F}^{\alpha}{\Phi}_{i}\) . Les fonctions introduites par l’enrichissement XFEM ne sont pas orthogonales aux fonctions de forme en chaque nœud XFEM. De plus, les fonctions d’enrichissement XFEM et les fonctions nodales, partagent le même support: il arrive des situations où les fonctions XFEM se rapprochent des fonctions nodales, au point de devenir quasiment colinéaires.
L’information sur la colinéarité est transportée dans la matrice de rigidité \(K\) : le conditionnement de la matrice de rigidité augmente si la colinéarité augmente au moins en un nœud XFEM. D’un point de vue plus formel, la matrice de rigidité (cf. r7.02.12-integrande-rigidite-mecanique dérive de la discrétisation d’une forme bilinéaire symétrique positive. Donc, il existe un produit-scalaire \({\langle .,.\rangle }_{K}\) tel que, le terme de la matrice de rigidité \({K}_{i,x}\) représentant la colinéarité entre la fonction nodale \({\Phi}_{i}\) et la fonction d’enrichissement \({F}_{x}{\Phi}_{i}\) , se réécrit:
En d’autres mots, la nouvelle base fonctionnelle \(\left({P}_{c}{\Phi}_{i},{P}_{c}{F}_{x}{\Phi}_{i}\right)\) (par le changement de base \({P}_{c}\) ) est orthogonale, sur le support de chaque nœud XFEM. Par conséquent, le nouveau système matriciel à inverser \(\stackrel{̃}{K}\stackrel{̃}{u}={P}_{c}^{T}f\) est mieux conditionné indépendamment de la position de l’interface, cf. r7.02.12-criteres-estimatifs-conditionnement.
La construction de \(\stackrel{̃}{K}\) n’est pas unique. [bib28] propose d’utiliser la factorisation de Cholesky des matrices de rigidité locales, associées à chaque nœud XFEM:
La matrice \({K}_{\mathit{loc}}^{i}\) est symétrique définie positive. Elle admet donc une factorisée de Cholesky, c’est-à-dire, qu’il existe une matrice triangulaire supérieure \({S}_{i}\) telle que:
(4299)#\[\begin{split}{P}_{c}^{i}=\lbrace \begin{array}{c}{S}_{i}^{-1}\text{si}i\text{est un noeud XFEM}\\ {I}_{d}\text{sinon}\end{array}\end{split}\]
avec \({S}_{i}^{-1}\) l’inverse de la matrice factorisée de Cholesky. La nouvelle matrice de rigidité \(\stackrel{̃}{K}\) vérifie à l’échelle locale du support du nœud XFEM:
\({I}_{d}\) désigne la matrice identité. \({I}_{d}\) respecte donc la condition d’orthogonalité de l’équation (4295).
Dans la pratique, il est nécessaire de mettre à l’échelle la nouvelle matrice \({\stackrel{̃}{K}}_{\mathit{loc}}^{i}\) avec le reste de la matrice de rigidité. On préfère donc le pré-conditionneur suivant au pré-conditionneur :
Par ailleurs, l’utilisation d’une factorisation de Cholesky est risquée si le conditionnement de la matrice de rigidité locale \({K}_{\mathit{loc}}^{i}\) se dégrade. En cas d’échec de la factorisation de Cholesky, on orthogonalise la matrice locale à l’aide d’une SVD (décomposition en valeurs singulières):
En résumé, la construction du pré-conditionneur de Béchet [bib28] procède en 4 étapes:
l’extraction des matrices locales \({K}_{\mathit{loc}}^{i}\) dans la matrice de rigidité (c’est-à-dire des matrices blocs associées aux ddls portés par les nœuds XFEM n° \(i\) ).
le calcul des matrices de pré-conditionnement locales \({P}_{c}^{i}\) (cf. équation (4301) ou (4304)),
l’assemblage du pré-conditionneur (cf. équation (4298)),
enfin, la transformation du système matriciel (cf. équations (4293) et (4294)).
Remarque:
Notons que la transformation du système matriciel présente une difficulté informatique notable dans le code_aster, compte-tenu de la complexité de la structure de données de stockage des matrices (cf. d4.06.10 et d4.06.07).
Une attention particulière doit être portée lors de l’intégration numérique des termes de rigidité et de second membre d’un élément traversé par la fissure. En effet, sur un élément traversé par une fissure, les gradients des déplacements peuvent être discontinus, et dans ce cas l’intégration numérique de Gauss-Legendre sur la totalité de l’élément n’est pas applicable. Afin de se replacer dans des conditions de régularité classiques, il convient de procéder à une intégration sur des domaines où l’intégrande est au moins continue. Pour un élément traversé par une fissure, il faut donc intégrer séparément de part et d’autre de la fissure (ceci apparaît pour la 1ère fois dans [bib24] pour le 2D et dans [bib32] pour le 3D). Plusieurs procédures sont envisageables, et facilement mises en œuvre en 2D. Les difficultés apparaissent avec le 3D.
On cherche à sous-découper en sous-tétraèdres un élément volumique quelconque (tétraèdre, pentaèdre, hexaèdre) coupé par une surface. On rappelle que cette découpe ne sert qu’à des fins d’intégration, elle est purement virtuelle et aucun nœud n’est ajouté au maillage. Le maillage n’est en rien modifié.
Pour les pentaèdres et les hexaèdres, le nombre de possibilités étant trop important et les configurations trop compliquées, on préfère se ramener au découpage de tétraèdres. On «condense» ainsi le grand nombre de possibilités de découpe à quelques configurations. En d’autres mots, nous formulons le «clustering» d’un ensemble de grande taille (nombre de configurations d’intersection) vers un ensemble de petite taille (nombre réduit de configurations de découpe dans un tétraèdre).
En fait, l’algorithme de «condensation» construit une application bijective, du polyèdre parent, vers quelques configurations de découpe de référence. Cette application peut-être décomposée en deux opérations bijectives:
1 – Une bijection (explicite) vers la découpe d’un tétraèdre: On réalise donc une phase préalable qui consiste à découper de manière systématique les quadrangles, pyramides, pentaèdres et hexaèdres en tétraèdres. Cette bijection n’est pas unique. Comme on peut le voir sur la Fig. 367, il existe deux bijections pour diviser un quadrangle en deux triangles (de même pour une pyramide), 6 bijections pour diviser un pentaèdre en 3 tétraèdres et 6 bijections pour diviser un hexaèdre en 2 pentaèdres qui sont ensuite divisés en tétraèdres. Dans le cas de l’hexaèdre, on obtient ainsi deux fois chaque configuration de découpe. Le nombre maximal de bijections distinctes pour diviser un hexaèdre en 6 tétraèdres est donc de \(\frac{{6}^{3}}{2}=108\) .
Fig. 366 Division d’un quadrangle en deux triangles (à gauche), division d’un pentaèdre en 3 tétraèdres (au milieu) et division d’un hexaèdre en 2 pentaèdres (à droite)#
Dans les r7.02.12-phase-prealable-decoupe-hexaedres à r7.02.12-phase-prealable-decoupe-pyramides , on présente la bijection retenue par défaut pour la division des éléments non simpliciaux en tétraèdres (ou triangles pour le cas 2D).
Il s’agit d’une bijection, car on construit une application (inversible) qui associe la numérotation relative des nœuds sommets dans un sous-tétraèdre, à la numérotation absolue des nœuds dans le maillage. L’application directe est stockée dans la table de connectivité (voir d4.10.02: structures de données XFEM).
2 – Une bijection (implicite) vers une configuration de découpe de référence: d’une part, on identifie la configuration de découpe de référence correspondante au tétraèdre à découper, d’autre part, on retourne «géométriquement» le tétraèdre, pour le «superposer» à cette configuration de découpe de référence (par exemple Fig. 368). En clair, on identifie les nœuds du sous-tétraèdre aux nœuds de l’élément de découpe de référence (c’est l’application directe de la bijection). On «condense» alors les multiples possibilités de découpe, liées aux différentes possibilités d’ordonnancement des arêtes et des points d’intersection.
Par la suite, cette opération facilite l’ordonnancement des nœuds milieux des sous-tétraèdres quadratiques pour le stockage informatique. En effet, l’ordonnancement des nouveaux nœuds milieux ne dépend plus des faces ou des arêtes où ils sont calculés, mais d’une seule configuration de découpe de référence. Par conséquent, les nouveaux nœuds milieux calculés sont nécessairement ordonnancés de la même manière, dans le sous-tétraèdre de découpe de référence. Ensuite, ces nœuds milieux sont positionnés sur les arêtes du tétraèdre par l’application inverse (de l’identification des nœuds ci-dessus). La bijection utilisée est donc tacite.
Les procédures d’ordonnancement et de calcul des points milieux sont explicitées ci dessous.
Phase préalable de découpe des hexaèdres pour se ramener à des tétraèdres#
Fig. 367 Division d’un hexaèdre (hexa8) en tétraèdres#
Un hexaèdre se divise alors en 6 tétraèdres, répertoriés dans le tableau suivant:
On note que ce découpage est le choix effectué par défaut mais on aurait pu choisir une autre manière de découper un hexaèdre en 2 pentaèdres ainsi qu’une autre façon de découper un pentaèdre en 3 tétraèdres.
Si l’élément est quadratique (HEXA20, PENTA15, …), on utilise la même subdivision en tétraèdres, que celle explicitée ci-dessus.
Par ailleurs, il faut tenir compte des nœuds milieux sur les nouvelles arêtes générées par la découpe: ces nœuds milieux coïncident avec les nœuds milieux de l’élément complet. On rajoute donc virtuellement à l’élément, les nœuds milieux nécessaires pour retrouver l’élément complet de la figure Fig. 368.
Fig. 368 Division d’un hexaèdre quadratique (hexa 20) via l’élément complet (hexa 27)#
Phase préalable de découpe des pentaèdres pour se ramener à des tétraèdres#
Fig. 369 Schéma de division d’un pentaèdre en tétraèdres#
Un pentaèdre se divise alors en 3 tétraèdres, répertoriés dans le tableau suivant:
Pour un pentaèdre quadratique ( penta15 ) on conserve la subdivision ci-dessus.
Par ailleurs, il faut tenir compte des nœuds milieux sur les nouvelles arêtes générées par la découpe: ces nœuds milieux coïncident avec les nœuds milieux de l’élément complet. On rajoute donc à l’élément, les nœuds milieux nécessaires pour retrouver l’élément complet Fig. 370.
Pour un pyramide quadratique ( pyra13 ) on conserve la subdivision ci-dessus.
Par ailleurs, il faut tenir compte des nœuds milieux sur les nouvelles arêtes générées par la découpe: ces nœuds milieux coïncident avec les nœuds milieux de l’élément complet. On rajoute donc à l’élément, les nœuds milieux nécessaires pour retrouver l’élément complet Fig. 370.
Le tétraèdre de référence est défini sur la Fig. 373. On détermine les points d’intersection \(Pi\) entre la surface \({\mathrm{lsn}}^{h}=0\) et les arêtes du tétraèdre.
Soit \(n\) le nombre de points d’intersection \(Pi\) .
À chaque point d’intersection \(Pi\) , on associe deux entiers: \(\mathit{Ai}\) et \(\mathit{NSi}\)
\(\mathit{Ai}\) est le numéro de l’arête sur laquelle se trouve \(Pi\) (par exemple si \(Pi\) se trouve sur l’arête d’extrémités \(\mathit{N2}-\mathit{N3}\) , alors \(\mathit{Ai}=4\) ). Dans le cas où \(Pi\) coïncide avec un nœud sommet du tétraèdre, Ai vaut 0,
\(\mathit{NSi}\) est le numéro du nœud sommet dans le cas où \(Pi\) coïncide avec un nœud sommet du tétraèdre (par exemple si \(Pi\) coïncide avec \(\mathit{N3}\) , alors \(\mathit{NSi}=3\) ). Dans le cas où \(Pi\) se trouve sur une arête, \(\mathit{NS}\) vaut 0.
Remarque
Le produit de \(\mathit{Ai}\) par \(\mathit{NSi}\) vaut toujours 0.
Les points d’intersection sont ensuite triés suivant l’ordre croissant des \(\mathit{Ai}\) . Les points d’intersection coïncidant avec des nœuds sommet vont donc se retrouver au début de la liste.
L’approximation de la level set utilisant les fonctions de forme du tétraèdre, la surface \({\mathit{lsn}}^{h}=0\) est alors un plan avec des éléments linéaires, et une surface de géométrie éventuellement courbe, avec des éléments quadratiques.
Le problème se ramène donc à la découpe d’un tétraèdre par une surface. Examinons les différents cas possibles suivant la valeur de \(n\) (nombre de points d’intersection \(Pi\) ). On peut déjà éliminer les cas triviaux où aucun sous-découpage n’est nécessaire:
lorsque \(n<3\) la trace de la surface dans le tétraèdre est un sommet ou une arête. Si la surface \({\mathit{lsn}}^{h}=0\) n’est pas un plan, la géométrie de la level-set est alors approximée par l’arête du tétraèdre.
lorsque \(n=3\) et que les 3 points d’intersection sont des points sommets, la trace de la surface dans le tétraèdre est une face du tétraèdre. Si la surface \({\mathrm{lsn}}^{h}=0\) n’est pas un plan, la géométrie de la level-set est alors approximée par la face du tétraèdre.
Dans ces deux cas, on obtient un seul sous-tétra, qui correspond au tétraèdre.
Trois points d’intersection dont deux points sommets#
\(\mathit{P1}\) et \(\mathit{P2}\) sont forcément les deux points sommet (donc \(\mathit{A1}=\mathit{A2}=0\) ). D’après \(\mathit{A3}\) , numéro d’arête correspondant à \(\mathit{P3}\) , on peut déterminer les 2 extrémités de cette arête, soient \(\mathit{E1}\) et \(\mathit{E2}\) .
Fig. 375 Cas général où \(n=3\) dont 2 points sommet#
On obtient un sous-tétraèdre de chaque coté de l’interface, soit au final 2 sous-tétraèdres.
Pour réduire les possibilités de découpe, (puisqu’il existe au moins 6 possibilités de coupe compte-tenu des 6 arêtes dans le sous-tétraèdre) on construit une bijection vers un unique élément de découpe Fig. 376. Les nœuds \(A\) , \(B\) , \(C\) , \(D\) , \(E\) , \(F\) , \(G\) , \(H\) définissent les inconnues d’ordonnancement de la découpe du sous-tétraèdre. Les inconnues d’ordonnancement vont être associées de manière unique aux nœuds de l’élément parent. Cette association construit une bijection (implicite).
Rappelons que,
l’ordonnancement est important d’un point de vue informatique, pour assurer le stockage de la maille virtuelle de sous-découpe. En particulier, pour localiser les points d’intersection et des points milieux dans les structures de données XFEM,
les nœuds de numéro compris entre 1000 et 1999 représentent les points d’intersections, ces nœuds sont calculés avec la procédure décrite au §4.2.3 de d4.10.02.
les nœuds de numéro supérieur à 2000 représentent les nouveaux nœuds milieux issus de la découpe, ces nœuds sont calculés avec la procédure décrite au §4.2.3 de d4.10.02
Remarque
Pour cette configuration de découpe particulière, la construction de la bijection n’est pas entièrement déterministe. En effet, les 2 sous-tétraèdres sont symétriques par rapport au plan de découpe. D’un point de vue topologique, les inconnues { \(A\) , \(D\) , \(E\) } sont symétriques aux inconnues { \(B\) , \(F\) , \(G\) }. Pour rendre déterministe la découpe d’un maillage donné, nous avons choisi la convention arbitraire suivante: prendre \(A\) comme le premier nœud de l’arête du premier point d’intersection.
Par analogie aux configurations moléculaires, la découpe va potentiellement générer deux configurations dites «énantiomères».
Fig. 376 Elément de découpe pour la configuration 4 points d’intersection dont 2 nœuds sommets#
\(P1\) est forcément le point sommet (donc \(A1=0\) et \(A2\ne 0\) ). Les arêtes \(A2\) et \(A3\) ont un nœud en commun, \(E1\) , et 2 nœuds différents \(E2\) et \(E3\) .
Fig. 377 Cas général où \(n=3\) dont un point sommet#
On obtient un sous-tétraèdre d’un coté, et de l’autre une pyramide que l’on divise en 2 sous-tétraèdres, soit au final 3 sous-tétraèdres.
Pour réduire les possibilités de découpe Fig. 378, on construit une bijection vers un unique élément de découpe . Les nœuds \(A\) , \(B\) , \(C\) , \(E\) , \(F\) , \(G\) , \(H\) définissent les inconnues d’ordonnancement de la découpe du sous-tétraèdre. Les inconnues d’ordonnancement vont être associées de manière unique aux nœuds de l’élément parent. Cette association construit une bijection (implicite).
Fig. 378 Elément de découpe pour la configuration 3 points d’intersection dont 1 nœud sommet#
Cinq points d’intersection dont deux points sommets#
Ce cas de découpe correspond à une configuration rasante. C’est un cas particulier de la découpe précédente pour lequel la surface \(\mathit{lsn}=0\) passe en plus par un autre nœud sommet. La découpe est alors strictement identique à la découpe précédente (cf. r7.02.12-trois-points-intersection-un-point-sommet). Cette configuration de découpe est représentée Fig. 379 avec à droite la configuration à 3 points d’intersection et un point sommet à laquelle on se ramène .
Fig. 379 Configuration rasante à 5 points d’intersection (à gauche) et configuration saine correspondante (à droite)#
Trois points d’intersection dont aucun point sommet#
On distingue au moins 4 possibilités de coupe, répertoriées Fig. 380 . Ces différentes possibilités ne sont pas traitées au cas par cas.
On construit une bijection vers un unique élément de découpe Fig. 378. Les nœuds \(A\) , \(B\) , \(C\) , \(D\) , \(E\) , \(F\) , \(G\) définissent les inconnues d’ordonnancement de la découpe du sous-tétraèdre. Les inconnues d’ordonnancement vont être associées de manière unique aux nœuds de l’élément parent. Cette association construit une bijection (implicite).
Exemple de procédure d’identification:
Dans chacune des configurations listéesFig. 380, on associe les inconnues d’ordonnancement n œ uds\(\lbrace A,B,C,D\rbrace\)aux n œ uds\(\lbrace {N}_{1,}{N}_{2,}{N}_{3,}{N}_{4}\rbrace\)du sous tétraèdre. Comme le nœud\(A\)appartient aux trois arêtes coupées, on l’identifie facilement au nœud correspondant dans chacune des configurations listées. Ensuite, on identifie alors\(\lbrace B,C,D\rbrace\)puisque les n œ uds\(B\),\(C\),\(D\)appartiennent respectivement à la première, deuxième, troisième arête intersectée, et sont différents du nœud\(A\)précédemment identifié. Il en résulte le tableau de correspondanceTableau 83
Cette procédure d’identification permet de construire implicitement une bijection.
Fig. 380 Configurations d’un tétraèdre avec 3 arêtes intersectées#
Fig. 381 Elément de découpe pour la configuration 3 points d’intersection et aucun nœud sommet#
Tableau 83 Identification des nœuds de l’élément de découpe aux nœuds du sous-tétraèdre dans chaque configuration#
Quatre points d’intersection dont un point sommet#
Ce cas de découpe correspond à une configuration rasante. C’est un cas particulier de la découpe précédente pour lequel la surface \(\mathit{lsn}=0\) passe en plus par l’un des nœuds sommets. La découpe est alors strictement identique à la découpe précédente (cf. r7.02.12-trois-points-intersection-aucun-point-sommet) . Cette configuration de découpe est représentée Fig. 382 avec à droite la configuration à 3 points d’intersection à laquelle on se ramène.
Fig. 382 Configuration rasante à 4 points d’intersection (à gauche) et configuration saine correspondante (à droite)#
Six points d’intersection dont deux points sommet#
Ce cas de découpe correspond à une configuration rasante. C’est encore un cas particulier de la découpe précédente pour lequel la surface \(\mathit{lsn}=0\) passe en plus par deux des nœuds sommets. La découpe est alors strictement identique à la découpe précédente (cf. r7.02.12-trois-points-intersection-aucun-point-sommet). Cette configuration de découpe est représentée Fig. 383 avec à droite la configuration à 3 points d’intersection à laquelle on se ramène .
Fig. 383 Configuration rasante à 6 points d’intersection (à gauche) et configuration saine correspondante (à droite)#
On distingue au moins 3 possibilités de coupe, répertoriées Fig. 384. Ces différentes possibilités ne sont pas traitées au cas par cas.
On construit une bijection vers un unique élément de découpe Fig. 385. Les nœuds \(A\) , \(B\) , \(C\) , \(D\) , \(E\) , \(F\) ,définissent les inconnues d’ordonnancement de la découpe du sous-tétraèdre. Les inconnues d’ordonnancement vont être associées de manière unique aux nœuds de l’élément parent. Cette association construit une bijection (implicite).
Fig. 384 Configurations d’un tétraèdre avec 4 arêtes intersectées#
Fig. 385 Elément de découpe pour la configuration 4 points d’intersection et aucun nœud sommet#
Ce cas de découpe correspond à une configuration rasante. C’est un cas particulier de la découpe précédente pour lequel la surface \(\mathit{lsn}=0\) passe en plus par l’un des nœuds sommets. La découpe est alors strictement identique à la découpe précédente (cf. r7.02.12-quatre-points-intersection). Cette configuration de découpe est représentée Fig. 386 avec à droite la configuration à 4 points d’intersection à laquelle on se ramène.
Fig. 386 Configuration rasante à 5 points d’intersection (à gauche) et configuration saine correspondante (à droite)#
Quatre points d’intersection dont un point sommet#
Cette configuration de découpe correspond à une configuration rasante. Elle ne se ramène à aucune configuration précedemmment décrite. On la distingue de la configuration décrite dans la r7.02.12-quatre-points-intersection-un-point-sommet) (qui comporte également 4 points d’intersection dont un point sommet) car le nœud sommet non intersecté sur l’arrête rasante n’appartient qu’à deux arêtes intersectées (au lieu de 3 dans le cas décrit dans la r7.02.12-quatre-points-intersection-un-point-sommet). On distingue au moins 24 possibilités de coupe, dont 3 sont répertoriées dans la Fig. 387. Ces différentes possibilités ne sont pas traitées au cas par cas.
Fig. 387 Configurations d’un tétraèdre avec 4 points d’intersection dont un point sommet#
On construit une bijection vers un unique élément de découpe Fig. 388. Les nœuds \(A\) , \(B\) , \(C\) , \(D\) , \(E\) , \(F\) ,définissent les inconnues d’ordonnancement de la découpe du sous-tétraèdre. Les inconnues d’ordonnancement vont être associées de manière unique aux nœuds de l’élément parent. Cette association construit une bijection (implicite). Cette configuration de découpe donnent naissance à 5 sous-tétraèdres enfants.
Fig. 388 Elément de découpe pour la configuration 4 points d’intersection dont un nœud sommet#
Lorsque l’on veut modéliser des jonctions, des intersections ou simplement que deux fissures sont assez proches pour découper le même élément, il faut être capable de diviser l’élément en domaines qui respectent toutes les discontinuités introduites.
La stratégie retenue consiste à découper l’élément plusieurs fois séquentiellement. C’est une stratégie qui a le mérite d’être assez rapide à implémenter, car la découpe d’un tétraèdre de référence par une fissure engendre des sous-tétraèdres qui peuvent à leur tour être considérés comme des tétraèdres de référence pour la découpe par la fissure suivante.
Le problème est que l’on n’optimise pas le nombre de sous-éléments total engendré, qui risque d’être très élevé si on redécoupe plus de 3 fois en 3D. Pour résoudre ce problème, il faudrait pouvoir découper directement n’importe quel type d’élément (hexaèdre, pentaèdre, pyramide, tétraèdre) en combinaison de tous ces type d’éléments. Cela nécessite des développements assez lourds et n’est donc pas envisageable. Une autre stratégie consisterait à regrouper les sous-éléments par zone et à trouver le nombre de points de Gauss (et les poids) optimal par zone. On ne stockerait plus les sous-éléments mais directement les zones avec les points de gauss associés.
Afin de dimensionner correctement les structures de données relatives au sous-découpage, il convient de déterminer le nombre maximum de sous-éléments engendrés par la phase de sous-découpage, suivant le type de la maille initiale.
On considère dans ce paragraphe que l’élément n’est découpé que par une seule fissure.
Un pentaèdre étant au préalable divisé en trois tétraèdres, on pourrait penser que le plus grand nombre de sous-éléments engendrés est celui où chacun de ces trois tétraèdres est à son tour sous-découpé en 6 sous-éléments; ce qui entraînerait un découpage final en 18 sous-éléments. Cependant, un tel cas de figure est impossible. Au maximum, sur les trois tétraèdres, deux seront découpés en 6 sous-éléments et un seul sera découpé en 4 sous-éléments, aboutissant à 16 sous-éléments.
Comme précédemment, le cas de figure où tous les six tétraèdres sont chacun sous-découpés en 6 sous-éléments est impossible. Le cas maximum est celui qui correspond au cas évoqué paragraphe précédent: les deux pentaèdres déduits de l’hexaèdre (voir Fig. 367) sont sous-découpés chacun en 16 sous-éléments. Le nombre de sous-éléments est donc \((6+6+4)+(6+6+4)=32\) .
Un élément fini XFEM peut être découpé par 4 fissures au maximum dans code_aster. Le nombre maximum de sous éléments générés est alors difficile à évaluer et un majorant du nombre maximal de sous élément serait trop élevé. En effet, si on prend le nombre maximal de sous élément généré par une seule fissure, soit 32, et que l’on considère que les 32 tétraèdres peuvent donner naissance à 6 tétraèdres pour chacune des fissures suivantes, on obtient \(32\ast \mathrm{6³}=6912\) . Ce nombre de sous éléments d’intégration pour un hexaèdre est tout à fait inatteignable car il va de soi que les tétraèdres issus du 1er découpage ne seront pas tous découpés par les fissures suivantes, mais il est difficile d’exhiber un meilleur majorant. On fait le choix de fixer le nombre maximal de sous éléments pour les éléments multi-fissurés à 4 fois le nombre maximal de sous éléments pour une seule fissure soit \(32\ast 4=128\) en 3D et \(4\ast 6=24\) en 2D. Ces nombres peuvent aisément être dépassés pour des éléments coupés par 3 ou 4 fissures. Dans ce cas, la procédure de découpage renvoie une erreur et l’utilisateur est invité à raffiner le maillage utilisé afin d’éviter que des éléments finis ne soient striés de fissures, ce qui entrainerait par ailleurs d’importants problèmes de conditionnement (voir r7.02.12-conditionnement-enrichissement).
On utilise une méthode comparable au découpage 3D pour le découpage 2D. Les quadrangles seront subdivisés en triangles eux mêmes sous découpés en fonction du passage de la fissure.
On découpe les quadrangles en 2 triangles:
Fig. 390 Sous–découpage d’un quadrangle en triangles : linéaire (à gauche), quadratique (à droite)#
Puis on recoupe les triangles en fonction du passage de la fissure. On effectue le même tri qu’en 3D pour les points d’intersection \(Pi\) .
Remarque :
Le nœud\(N9\)n’est pas un nœud du maillage. Il est rajouté pour faire la connectivité entre le quadrangle et les triangles quadratiques. Il est défini comme étant le milieu du segment\(N2N4\)(et pas le nœud central du quadrangle). La valeur de la level set normale qui lui est attribuée est obtenue par interpolation des\(\mathit{lsn}\)des nœuds\(N1\)à\(N8\) *. Un réajustement peut alors être nécessaire pour assurer l’unicité de la solution de l’équation* \(\mathit{lsn}=0\)le long de la diagonale (cf. réajustement des level setr7.02.12-reajustement-level-sets) lors de la découpe des éléments enfants en sous éléments d’intégration. Préalablement à toute découpe, on effectue donc une boucle sur les arêtes des éléments enfants (cela a déjà été effectué pour les arêtes des éléments parents lors du réajustement des level set) pour détecter les situations pathogènes et effectuer les réajustements qui s’imposent, exactement comme décrit au r7.02.12-reajustement-level-sets.
Lorsque l’on a 1 ou 2 points d’intersection qui tombe uniquement sur les sommets, aucun découpage n’est nécessaire.
Si 1 ou 2 des points d’intersection ne tombent pas sur un sommet, on recoupe les triangles.
Si on a 1 point sommet et 1 point non sommet, \(\mathit{P1}\) est forcement sommet. On détermine les sommets de l’arête du 2e point \(\mathit{N2}\) et \(\mathit{N3}\) pour découper le triangle en 2 sous triangles \(\mathit{P1}\mathit{N2}\mathit{P2}\) et \(\mathit{P1}\mathit{N3}\mathit{P2}\) . Si le triangle est quadratique, on détermine les points milieux de l’arête du 2e point \(\mathit{Q1}\) et \(\mathit{Q2}\) et celui de la fissure \(\mathit{Q3}\) .
Fig. 392 Cas de découpage linéaire (à gauche) et quadratique (à droite) avec 2 points d’intersection dont 1 sommet#
Avec 2 points non sommet, on doit découper en 3 triangles. On obtient les triangles \(\mathit{N1}\mathit{P1}\mathit{P2}\) , \(\mathit{P1}\mathit{P2}\mathit{N3}\) et \(\mathit{P1}\mathit{N2}\mathit{N3}\) . \(\mathit{N1}\) étant le point proche des 2 points d’intersections d’après le tri effectué. Si le triangle est quadratique, on détermine les points milieux de l’arête du 1er point \(\mathit{Q1}\) et \(\mathit{Q2}\) , ceux du 2e point \(\mathit{Q3}\) et \(\mathit{Q4}\) , celui de la fissure \(\mathit{Q5}\) et celui de la nouvelle arête \(\mathit{Q6}\) .
Fig. 393 Cas de découpage linéaire (à gauche) et quadratique (à droite) avec 2 points d’intersection sans sommet#
Dans les r7.02.12-phase-prealable-decoupe-hexaedres à r7.02.12-sous-decoupage-2d, nous avons expliqué le découpage géométrique des éléments en sous-cellules d’intégration (Tetra10 en 3D, Tri6 en 2D), sans expliquer comment sont calculées les coordonnés des points d’intersection et des nouveaux nœuds milieux dans le maillage virtuel. Dans ce paragraphe, nous allons donc décrire les algorithmes de calcul des coordonnées des nouveaux points, générés par la découpe.
Pour chaque type de points, il existe un algorithme de positionnement spécifique:
Le positionnement des points d’intersection permet de délimiter une frontière explicite, bien que l’interface/fissure soit définie par une équation implicite \({\mathrm{lsn}}^{h}=0\) . Précisément, le lieu des points d’intersection est calculé en prenant l’intersection entre les arêtes du maillage et l’iso-zéro de la level-set. Les points d’intersection sont ensuite stockés dans une structure de donnée auxiliaire, en complétion des points du maillage d4.10.02.
Le positionnement des points milieux permet d’améliorer la description d’une level-set courbe, dans le cas d’éléments quadratiques. En présence d’éléments à bords courbes, la courbure des arêtes de l’élément quadratique est conservée, en positionnant des nœuds milieux rigoureusement sur les arêtes lors de la découpe. En revanche, une découpe avec des éléments à bords droits altérerait la courbure de l’élément originel. Pour capter la courbure de la level-set et des bords courbes d’éléments quadratiques, on calcule donc plusieurs types de points milieux (voir r7.02.12-calcul-points-milieu). Plusieurs algorithmes sont implémentés, pour capter les différents types de points milieux envisageables. Ces nouveaux points milieux sont alors stockés dans une structure de donnée auxiliaire, en complétion des points du maillage d4.10.02.
Par la suite, on adopte la notation suivante:
\(\xi ={\left({\xi}_{i}\right)}_{i\in [1,\mathit{dim}]}\) désigne les coordonnées de référence dans l’élément parent. Ces coordonnées varient généralement entre \(\left[0,1\right]\) ou entre \(\left[-1,1\right]\) . Rappelons que les fonctions de formes ou polynômes d’interpolation, sont paramétrées dans Aster, grâce aux coordonnées de référence r3.01.01,
\(X={\left({X}_{i}\right)}_{i\in [1,\mathit{dim}]}\) désigne les coordonnées d’un point dans l’espace réel. Si le point courant est un nœud, ce sont les coordonnées géométriques du nœud dans le maillage. Si le point courant n’est pas un nœud, il s’agit des coordonnées géométriques interpolées dans l’élément parent \(X(\xi )=\sum_{1⩽j⩽\mathit{nnop}}{X}^{j}{\Phi}_{j}(\xi )\) .
Moës propose de poser le problème de recherche des points d’intersection dans l’espace de référence parent. En clair, pour calculer les coordonnées réelles d’un point d’intersection, on recherche au préalable les coordonnées dans l’espace de référence \(\xi\) . On déduit ensuite les coordonnées réelles grâce à l’interpolation des coordonnées géométriques dans l’élément parent.
Pour calculer le point d’intersection \(P\) entre une arête \(\mathit{AB}\) et l’iso-zéro, on résout le problème suivant:
Fig. 394 Recherche du point d’intersection \({I}_{1}\) sur l’arête \({I}_{2}{I}_{3}\) qui est courbe dans l’espace de référence. On prend alors en compte la courbure de \({I}_{2}{I}_{3}\) avec le nœud milieu \({M}_{2}\)#
Exception: Pour les éléments quadratiques multi-fissurées, dans l’élément de référence parent, l’arête \(\mathit{AB}\) n’est pas nécessairement droite. Le nœud milieu \(M\) de l’arête \(\mathit{AB}\) n’est pas nécessairement sur le segment \([A,B]\) si le côté \(\mathit{AB}\) est le côté d’un sous élément d’intégration généré par la découpe conformément à une première fissure (voir Fig. 394). Il faut alors prendre en compte l’éventuelle courbure de l’arête; on a ainsi \(P\in \mathit{AB}\iff {\xi}_{P}=(\mathrm{2t}-1)(t-1){\xi}_{A}+t(\mathrm{2t}-1){\xi}_{B}+\mathrm{4t}(1-t){\xi}_{M}\) avec \(t\in [0,1]\) .
Du coup, le problème posé revient à rechercher \(t\in [0,1]\) tel que \({\mathit{lsn}}^{h}\left({\xi}_{P}=(\mathrm{2t}-1)(t-1){\xi}_{A}+t(\mathrm{2t}-1){\xi}_{B}+\mathrm{4t}(1-t){\xi}_{M}\right)=0\) .
Pour les modélisations linéaires ou quadratiques mono-fissurées, le point \(P\) appartient nécessairement à l’arête \([A,B]\) et donc: \({\xi}_{P}=(1-t){\xi}_{A}+t{\xi}_{B}\) avec \(t\in [0,1]\) .
Du coup, le problème posé, revient à rechercher \(t\in [0,1]\) tel que \({\mathit{lsn}}^{h}\left({\xi}_{P}=(1-t){\xi}_{A}+t{\xi}_{B}\right)=0\) .
On cherche alors la racine d’une équation polynomiale puisque les fonctions d’interpolation de la level-set sont des polynômes. Cette équation est alors résolue à l’aide de l’algorithme de Newton-Raphson (détaillé ci-dessous).
Comme \(f:t\to {\mathit{lsn}}^{h}(t)\) est \({C}^{2}\) et strictement monotone au voisinage du point d’intersection (sinon, la level-set devient parallèle à l’arête \(\mathit{AB}\) et ne l’intersecte pas), la convergence de l’algorithme de Newton est quadratique , donc optimale.
Algorithme de Newton-Raphson pour le calcul d’un point d’intersection de l’iso-zéro avec l’arête\(\mathit{AB}\):
Lecture des coordonnées de référence des nœuds sommets \(A\) et \(B\) et de l’éventuel nœud \(M\) milieu : par construction, ces nœuds sont des nœuds de l’élément enfant que l’on découpe et leur coordonnées sont donc connues.
Initialisation du Newton-Raphson:
Pour les modélisations linéaires, on initialise la recherche avec l’approximation linéaire le long de l’arête:
Pour les modélisations quadratiques, on initialise la recherche avec la solution exacte du problème le long de l’arête de l’élément enfant. Cette solution est donnée par la résolution d’une équation polynomiale du second degré. Les coefficients de ce polynôme du second degré sont donnés par les valeurs de \(\mathit{lsn}\) aux nœuds extrémités et milieu de l’arête du sous élément enfant. Étant donné que la \(\mathit{lsn}\) change de signe strictement entre les extrémités de cette arête, il existe une unique solution à l’équation \(\mathit{lsn}=0\) sur l’arête. L’équation à résoudre a la forme suivante:
Parmi les deux solutions éventuelles de cette équation du second degré on ne retient que celle qui satisfait: \({t}_{0}\in [0,1]\) . Cette valeur constitue une bonne initialisation et c’est même la solution finale du problème lorsque l’arête de l’élément enfant sur lequel on effectue la recherche coïncide avec une arête de l’élément parent. En effet, la l \(\mathit{lsn}\) est alors interpolée uniquement par les nœuds extrémités et milieu de l’arête de l’élément parent. Cette interpolation étant quadratique, l’initialisation fournit la solution exacte. En revanche, lorsque l’arête de l’élément enfant ne coïncide pas avec une arête de l’élément parent (arête interne Fig. 395), l’initialisation ne fournit pas la solution exacte du problème car la \(\mathit{lsn}\) est interpolée par tous les nœuds de l’élément parent le long de l’arête de l’élément enfant.
Fig. 395 Intersection de la courbe lsn=0 avec une arête de l’élément enfant#
Calcul de la fonction : \({\mathit{lsn}}^{h}({\xi}_{n})=\sum_{1⩽j⩽\mathit{nnop}}{\mathit{lsn}}_{j}{\Phi}_{j}({\xi}_{n})\)
Calcul de la dérivée: \({\mathit{dlsn}}^{h}({\xi}_{n})/\mathit{dt}=\sum_{1⩽j⩽\mathit{nnop}}{\mathit{lsn}}_{j}\nabla {\Phi}_{j}({\xi}_{n})\cdot \vec{v}\)
Vérification de la pente: si \(∣{\mathit{dlsn}}^{h}({\xi}_{n})/\mathit{dt}∣<{10}^{-16}\) alors sortie avec code erreur.
Calcul de l’incrément: \({\Delta}_{n}=\frac{{\mathit{lsn}}^{h}({\xi}_{n})}{{\mathit{dlsn}}^{h}({\xi}_{n})/\mathit{dt}}\)
Calcul de la nouvelle position: \({t}_{n+1}={t}_{n}-{\Delta}_{n}\)
Vérification de l’hypothèse \(\mathit{AB}\) : si \({t}_{n}\notin \left[0,1\right]\) alors sortie avec code erreur.
Critère d’arrêt sur l’incrément: si \(∣{\Delta}_{n}∣<{10}^{-8}\) alors fin de la boucle.
Calcul des coordonnées de référence du point d’intersection: \({\xi}_{\infty}=\left\lbrace \begin{array}{c}(1-{t}_{\infty}){\xi}_{A}+{t}_{\infty}{\xi}_{B}\text{}\mathit{en}\mathit{linéaire}\mathit{et}\mathit{mono}-\mathit{fissuré}\\ ({\mathrm{2t}}_{\infty}-1)({t}_{\infty}-1){\xi}_{A}+{t}_{\infty}({\mathrm{2t}}_{\infty}-1){\xi}_{B}+{\mathrm{4t}}_{\infty}(1-{t}_{\infty}){\xi}_{M}\text{}\mathit{en}\mathit{quadratique}\mathit{multi}-\mathit{fissuré}\end{array}\right.\)
Vérification de l’hypothèse \(P\in \mathit{ELREFP}\) : si \({\xi}_{\infty}\notin \mathit{ELREFP}\) alors sortie avec code erreur.
Calcul des coordonnées réelles par interpolation des coordonnées géométriques:
On teste l’incrément\({\Delta}_{n}\)pour évaluer la convergence du Newton, car la solution et l’incrément varient sur l’échelle de l’unité. La valeur de l’incrément par rapport à l’unité est donc un indicateur robuste de la convergence. Si le critère d’arrêt portait sur la valeur de la fonction\({\mathit{lsn}}^{h}({\xi}_{n})\), il faudrait mettre à l’échelle le résultat.
Très concrètement, quand l’incrément tend vers zéro par rapport à l’unité, la solution\({t}_{n+1}={t}_{n}-{\Delta}_{n}\)devient stationnaire par rapport à l’unité. D’un point de vue numérique, cela garantit la convergence sur l’intervalle\([0,1]\), compte-tenu de la précision machine.
A titre d’exemple, considérons maintenant l’interpolation d’une level-set dans l’élément quadratique (QUAD8) de la r7.02.12-approximation-level-sets. Envisageons un segment \(\mathit{AB}\) arbitraire dans l’élément parent (voir Fig. 396). En traçant l’évolution de la level-set interpolée \(f:t\to {\mathit{lsn}}^{h}(t)\) le long du segment \(A-B\) (voir Fig. 397), on constate que la level-set varie de manière strictement monotone au voisinage du point d’intersection \(P\) . Tant que l’on reste au voisinage du point \(P\) , l’algorithme de Newton respectera les conditions de convergence optimale.
L’algorithme de Newton proposé ci-dessus, est initialisé par une interpolation linéaire de la level-set le long du segment \(A-B\)Fig. 398. Le point d’initialisation \({P}_{0}\) est alors suffisamment proche du point de convergence \(P\) . L’algorithme de Newton aura un comportement optimal.
En revanche, si l’algorithme avait été initialisé au point \(A\) , la proximité d’un point de stationnarité (où la dérivée première de \(f:t\to {\mathit{lsn}}^{h}(t)\) s’annule) aurait causé l’échec du Newton. Plusieurs scénarios seraient envisageables:
soit le Newton ne convergerait pas au bout d’un nombre d’itérations maximal donné,
soit le Newton sortirait du segment \(A-B\) .
La présence du point de stationnarité a une justification géométrique: l’iso-zéro de la level-set est pratiquement parallèle au segment \(A-B\) au voisinage du point \(A\) (voir Fig. 396). Par la suite, la distance entre le segment \(A-B\) et l’iso-zéro diminue progressivement. Par conséquent, au voisinage du point d’intersection la variation de la level-set redevient strictement monotone.
Pour prévenir cette configuration, le Newton est initialisé suffisamment proche du point d’intersection Fig. 398. De plus, à chaque itération de l’algorithme, on contrôle la persistance de la solution sur le segment \(A-B\) c’est-à-dire \(\forall t,t\in [0,1]\) . Cela apporte alors une garantie supplémentaire sur la robustesse de la solution calculée, par le Newton.
Fig. 396 Recherche du point d’intersection entre l’iso-zéro de la level-set et un segment A-B (arbitraire) dans un élément QUAD8#
Fig. 397 Evolution de fonction distance (level-set) le long du segment A-B#
Fig. 398 Initialisation de l’algorithme de Newton-Raphson par interpolation linéaire de la level-set#
Lors de la découpe d’un tétraèdre quadratique, les nouveaux points milieux peuvent être regroupés suivant 3 catégories Fig. 399:
les points milieux de premier type \({\mathit{MP}}_{1}\) : situés sur les arêtes du tétraèdre (triangle), entre les points d’intersections et les nœuds sommets,
les points milieux de deuxième type \({\mathit{MP}}_{2}\) : situés sur l’iso-zéro, dans la face quadratique (SE3, TRI6 ou QUAD8) délimitée par les points d’intersection. Les points milieux sur les bords de la face et le point milieu au centre de la face (si face quadrangle) sont calculés avec des procédures légèrement différentes,
les points milieux de troisième type \({\mathit{MP}}_{3}\) : situés sur les faces quadrangles des sous-polyèdres-enfants (pentaèdre ou pyramide), générés par la découpe du tétraèdre.
Fig. 399 Les différents types de nœuds milieux \({\mathit{MP}}_{i}\) à calculer lors de la découpe des éléments quadratiques#
Définition des points milieux
On propose de définir systématiquement les points milieux dans le repère de référence de l’élément parent. Comme l’élément de référence est à bords droits et à faces planes, cela facilite grandement le calcul des points milieux notamment les nœuds milieux de premier type. Soulignons d’avantage ce point essentiel, les calculs à venir s’appuieront souvent implicitement sur l’hypothèse que les faces des tétraèdres à découper sont planes dans l’élément de référence parent.
Cependant, avec cette définition, il est possible de calculer un nœud milieu qui ne soit pas «au milieu» de l’arête dans le système de coordonnées réelles. Cela n’est pas gênant, à condition de choisir une technique d’intégration adaptée.
Dans l’élément de référence, nous proposons la définition suivante des points milieux:
pour un point milieu du premier type , situé entre le nœud sommet \(\mathit{XX}\) et le point d’intersection \(\mathit{IP}\) : \({\xi}_{{\mathit{PM}}_{1}}={\xi}_{\mathit{Milieu}(\mathit{XX}-\mathit{IP})}=\frac{{\xi}_{\mathit{XX}}+{\xi}_{\mathit{IP}}}{2}\)
Fig. 400 Calcul dans le repère de référence, d’un point milieu de deuxième type, situé entre 2 points d’intersection#
Exception: Pour les éléments quadratiques multi-fissurées, dans l’élément de référence parent, l’arête intersectée et sur laquelle on souhaite positionner les points milieux n’est pas nécessairement droite. Si l’arête intersectée est le côté d’un sous élément d’intégration généré par la découpe conformément à une première fissure (voir Fig. 401), il faut prendre en compte l’éventuelle courbure de l’arête pour positionner les points milieux de premier type:
Fig. 402 Calcul dans le repère de référence, d’un point milieu de deuxième type, situé entre 2 points d’intersection#
\(M\) est le milieu du segment joignant les 2 points d’intersection \({\xi}_{M}=\frac{{\xi}_{{\mathit{IP}}_{1}}+{\xi}_{{\mathit{IP}}_{2}}}{2}\)
\(\vec{v}\) est le vecteur unitairenormal à \(\vec{{\mathit{IP}}_{1}-{\mathit{IP}}_{2}}\) , appartenant à la face \((A,B,C)\) (plane dans l’élément de référence). Concrètement pour calculer le vecteur normal \(\vec{v}\) , on pose les équations suivantes:
Il y a 2 vecteurs directeurs unitaires possibles pour une droite donnée, ce qui explique les 2 solutions calculées pour le vecteur\(\vec{v}\). On choisit arbitrairement l’une des 2 solutions. Ce choix n’a aucune influence sur la convergence de l’algorithme de Newton.
pour un point milieu du deuxième type , situé au centre d’une face quadrangle, on opère le choix suivant:
\(M\) est le point milieu de l’arête \(\left[{\mathit{IP}}_{1}-{\mathit{IP}}_{4}\right]\) c’est-à-dire \({\xi}_{M}=\frac{{\xi}_{\mathit{IP1}}+{\xi}_{\mathit{IP4}}}{2}\) . Si l’iso-zéro est plane la découpe génère alors des sous-tétraèdres à bords droits.
\(\vec{v}\) est le gradient de la lsn au point M : \(\vec{v}=\vec{{\mathit{grad}}_{\mathit{lsn}}}(M)\) . C’est la «meilleure» direction de recherche car le gradient de la lsn est par définition orthogonal à la surface \(\mathit{lsn}=0\) que l’on cherche à atteindre.
Fig. 403 Calcul dans le repère de référence, d’un point milieu de deuxième type, situé sur une face quadrangle#
Calcul des points milieux de deuxième type :
Les points milieux de deuxième type vérifient l’équation d’intersection d’un droite avec l’iso-zéro suivante:
En fonction de la localisation du point milieu (soit situé entre 2 points d’intersection, soit situé à l’intérieur de la face quadrangle sur l’iso-zéro), nous avons explicité les coordonnées exactes de \({\xi}_{M}\) et \(\vec{v}\) dans l’élément de référence.
On résout la recherche du zéro de l’équation \({\mathit{lsn}}^{h}\left({\xi}_{M}+t\vec{v}\right)=0\) par un Newton comme à la r7.02.12-approximation-level-sets. Afin de s’assurer que l’algorithme de Newton «reste» respectivement dans la face TRIA pour un point milieu de deuxième type situé entre deux points d’intersection et dans le sous élément TETRA pour un point milieu de deuxième type situé au centre d’une face quadrangle, on renforce l’algorithme de Newton par la méthode de Dekker. Si une itération de l’algorithme de Newton nous envoie en dehors de la zone de recherche autorisée, une dichotomie sera préférée pour se rapprocher de la solution. En effet, il arrive que l’équation \(\begin{array}{c}{\xi}_{{\mathit{PM}}_{2}}={\xi}_{M}+t\vec{v}\\ {\mathit{lsn}}^{h}\left({\xi}_{{\mathit{PM}}_{2}}\right)=0\end{array}\) possède plusieurs solutions, mais il n’y en a qu’une seule qui est dans la face TRIA de recherche ou dans le sous élément TETRA de recherche. Cette solution peut s’avérer difficile à capter par l’algorithme de Newton, surtout lorsque la courbe \(\mathit{lsn}=0\) a une courbure importante ou «rase» une arête de la face TRIA ou une face du sous élément TETRA (voir Fig. 404 et Fig. 405).
Fig. 404 Existence de solutions rapprochées pour le point milieu de type 2 dans une face TRIA#
Fig. 405 Existence de solutions rapprochées pour le point milieu de type 2 situé sur une face quadrangle dans un sous élément TETRA#
Fig. 406 Limites de «l’espace de recherche» dans la face TRIA pour un point milieu de deuxième type situé entre deux points d’intersection#
Fig. 407 Limites de «l’espace de recherche» dans le sous élément TETRA pour un point milieu de deuxième type situé au centre d’une face quadrangle#
Algorithme de Newton-Raphson pour le calcul d’un point milieu situé sur l’isozéro de la level-set:
Définition de \({\xi}_{M}\) point de départ de la recherche et de \(\vec{v}\) direction de la recherche.
Initialisation du Newton-Raphson: \({t}_{0}=0\)
Détermination des bornes inférieures et supérieures \({t}_{min}\) et \({t}_{max}\) de l’espace de recherche (voir Fig. 406 et Fig. 407)
Tant que \(n⩽\mathit{itermax}\) faire:
\({\xi}_{n}={t}_{n}\vec{v}+{\xi}_{M}\)
Calcul de la fonction : \({\mathit{lsn}}^{h}({\xi}_{n})=\sum_{1⩽j⩽\mathit{nnop}}{\mathit{lsn}}_{j}{\Phi}_{j}({\xi}_{n})\)
Calcul de la dérivée: \({\mathit{dlsn}}^{h}({\xi}_{n})/\mathit{dt}=\sum_{1⩽j⩽\mathit{nnop}}{\mathit{lsn}}_{j}\nabla {\Phi}_{j}({\xi}_{n})\cdot \vec{v}\)
Vérification de la pente: si \(∣{\mathit{dlsn}}^{h}({\xi}_{n})/\mathit{dt}∣<{10}^{-16}\) alors sortie avec code erreur.
Calcul de l’incrément: \({\Delta}_{n}=\frac{{\mathit{lsn}}^{h}({\xi}_{n})}{{\mathit{dlsn}}^{h}({\xi}_{n})/\mathit{dt}}\)
Calcul de la nouvelle position: \({t}_{n+1}={t}_{n}-{\Delta}_{n}\)
Si \({t}_{n+1}\in [{t}_{min},{t}_{max}]\) , alors on continue, sinon:
Si \({t}_{n+1}>{t}_{max}\) , alors \({t}_{n+1}=\frac{{t}_{n}+{t}_{max}}{2}\)
Si \({t}_{n+1}<{t}_{min}\) , alors \({t}_{n+1}=\frac{{t}_{n}+{t}_{min}}{2}\)
Critère d’arrêt sur l’incrément: si \(∣{\Delta}_{n}∣<{10}^{-8}\) alors fin de la boucle.
Calcul des coordonnées de référence du point d’intersection: \({\xi}_{\infty}={t}_{\infty}\vec{v}+M\)
Vérification de l’hypothèse \({\xi}_{\infty}\in \mathit{ELREFP}\) : si \({\xi}_{\infty}\notin \mathit{ELREFP}\) alors sortie avec code erreur.
Calcul des coordonnées réelles par interpolation des coordonnées géométriques:
Pour les éléments non simpliciaux, il existe une situation pour lesquelles la recherche du point milieu de type 2 est vaine dans l’espace de recherche autorisé (entre les bornes \(inf\) et \(\sup\) ). Comme indiqué à la r7.02.12-phase-prealable-decoupe-hexaedres, la découpe en sous éléments se base sur une réduction des éléments non simpliciaux en éléments simpliciaux (TRIA en 2D et TETRA en 3D). Or cette réduction peut amener des situations qui rendent la recherche des points milieux de type 2 impossible dans dans la subdivision des éléments non simpliciaux. Un exemple est présenté sur la Fig. 409 à gauche. On recherche le point milieu de type 2 entre les points d’intersection \({\mathit{IP}}_{1}\) et \({\mathit{IP}}_{2}\) . Or dans l’intervalle de recherche autorisé (entre les bornes \(inf\) et \(\sup\) ), le champ de \(\mathit{lsn}\) est de signe constant. L’algorithme de Newton-Raphson est donc mis en échec. L’arête diagonale interne est en effet intersectée par la courbe \(\mathit{lsn}=0\) alors que cette dernière était détectée comme non coupée car la valeur de la \(\mathit{lsn}\) en chacun de ses 3 nœuds est strictement de même signe. Mais à la différence d’une arête de l’élément parent, le champ \(\mathit{lsn}\) le long de cette arête interne dépend de la valeur de la \(\mathit{lsn}\) en chacun des 8 nœuds de l’élément parent QUAD8, et non pas seulement de la valeur de la \(\mathit{lsn}\) en ses 3 nœuds. Voilà pourquoi le champ de \(\mathit{lsn}\) change de signe le long de cette arête interne. Cette difficulté peut être contournée en choisissant une autre bijection pour la division du quadrangle en deux triangles (confer Fig. 366). Sur la Fig. 409 à droite, on utilise l’autre bijection pour la découpe du quadrangle en deux triangles de manière à éviter un changement de signe du champ \(\mathit{lsn}\) le long de l’arête interne. On détermine alors les points milieux de deuxième type \({M}_{1}\) et \({M}_{2}\) dans les intervalles autorisés. Donc chaque fois que la recherche d’un point milieu de deuxième type dans l’intervalle \([inf,\sup]\) est mise en échec au sein d’un élément simplicial, on relance la procédure de découpe de l’élément en choisissant une bijection différente pour la division en sous cellule triangulaires (ou tétraédriques en 3D). Dans le cas où toutes les bijections possibles n’auraient pas permis d’apporter une solution à ce problème, on effectue alors une approximation linéaire de la courbe \(\mathit{lsn}=0\) , où le point milieu est choisi comme étant le milieu du segment \([{\mathit{IP}}_{1}{\mathit{IP}}_{2}]\) dans l’espace de référence.
Fig. 408 Cas d’avortement de la recherche d’un point milieu de type 2 pour les éléments non simpliciaux.#
Ces cas de figure sont rares et ne surviennent que lorsque la courbure de l’interface est importante au regard de la taille des mailles. En raffinant suffisamment le maillage, on est assuré de ne pas se heurter à ce problème.
pour un point milieu du troisième type , situé sur une face quadrangle d’un polyèdre-enfant (pyramide ou pentaèdre), l’arête virtuelle de découpe relie le premier point d’intersection \({\mathit{IP}}_{1}\) de la face et un nœud sommet du tétraèdre \(B\) (voir ). Concernant le positionnement du point milieu sur la face quadrangle:
si la level-set est droite, on opère le choix suivant: \({\xi}_{{\mathit{PM}}_{3}}={\xi}_{\mathit{Milieu}({\mathit{IP}}_{1}-B)}=\frac{{\xi}_{{\mathit{IP}}_{1}}+{\xi}_{B}}{2}\) pour que la découpe génère des sous-tétraèdres à bords droits.
si la level-set est courbe, on opère le choix suivant: \({\xi}_{{\mathit{PM}}_{3}}={\xi}_{\mathit{Milieu}\left({\mathit{PM}}_{2}-D\right)}=\frac{{\xi}_{{\mathit{PM}}_{2}}+{\xi}_{D}}{2}\) pour que le point milieu calculé ne sorte pas de la face quadrangle .
Fig. 409 Calcul dans le repère de référence, d’un point milieu de troisième type pour une iso-zéro de forme courbe#
Remarque sur le critère de courbure
Pour juger de la courbure de la lsn, on effectue la démarche suivante. On considère d’abord le point milieu\({\mathit{PM}}_{3}\)candidat tel que\({\xi}_{{\mathit{PM}}_{3}}={\xi}_{\mathit{Milieu}({\mathit{IP}}_{1}-B)}=\frac{{\xi}_{{\mathit{IP}}_{1}}+{\xi}_{B}}{2}\) *. On calcule alors le sinus de l’angle entre les vecteurs* \(\vec{{\mathit{IP}}_{1,}{\mathit{PM}}_{2}}\)et\(\vec{{\mathit{IP}}_{1,}{\mathit{PM}}_{3}}\), le plan étant orienté par rapport à la normale sortante venant vers le lecteur. S’il est négatif, ce point milieu donnera naissance à des éléments distordus, sinon les deux TRIA résultants seront sains. En pratique, on considère que la lsn est trop courbe lorsque le sinus est inférieur à 0.1. On choisit alors\({\mathit{PM}}_{3}\)tel que.
Fig. 410 Critère pour déterminer si la lsn est jugée courbe: situation «courbe» à gauche et «saine» à droite#
Sous-découpage et erreur de discrétisation
Dans la procédure de découpage expliquée à la r7.02.12-algorithmes-sous-decoupage, la levet-set est interpolée systématiquement dans l’élément de référence parent. Par conséquent, les points d’intersection et points milieux positionnés sur l’iso-zéro, vérifient rigoureusement l’équation \({\mathit{lsn}}^{h}({\xi}^{\mathit{ELREFP}})=0\) . En clair, l’algorithme de positionnement des points sur la surface de l’iso-zéro n’introduit pas d’erreur supplémentaire par rapport à l’erreur d’interpolation de la level-set.
En revanche, la forme géométrique de la level-set est approximée en reliant les points calculés pour former la frontière des sous-domaines d’intégration (voir , où la fonction \({\Phi}^{h}={\mathit{lsn}}^{h}({\xi}^{\mathit{ELREFP}})=0\) est une conique alors que le bord de la sous-cellule quadratique coïncidant au mieux avec cette conique est parabolique).
La récupération des facettes de contact est utile à la fois pour le contact mais aussi pour appliquer une pression mécanique dans les lèvres de fissures. Pour obtenir ces facettes, on s’appuie sur le travail de découpage en sous éléments décrit à la r7.02.12-sous-decoupage. En effet, la topologie des sous éléments d’intégration contient déjà toutes les informations sur les facettes de contact, qui ne sont autres que les côtés (en 2D) ou les faces (en 3D) des sous éléments d’intégration qui vérifient \(\mathit{lsn}=0\) . Il suffit donc simplement de parcourir les côtés ou les faces des sous éléments d’intégration et de sélectionner ceux qui sont confondus avec la discontinuité (voir ). Or les nœuds sommets des sous éléments d’intégration qui sont des points d’intersection sont repérés par leur connectivité comprise entre 1000 et 1999. Dans le cas d’un point d’intersection confondu avec un nœud de l’élément parent, on aura directement \(\mathit{lsn}=0\) pour ce nœud.
Les facettes que l’on récupère sont des SEG2en 2D linéaire, des SEG3en 2D quadratique, des TRIA3 en 3D linéaire et des TRIA6 en 3D quadratique. Afin de ne récupérer les facettes qu’une et une seule fois, on n’effectue la recherche que sur les côtés ou les faces des sous éléments d’intégration tels que \(\mathit{lsn}<0\) .
Fig. 412 Sous éléments d’intégration et facettes de contact dans un QUAD8#
Dans l’exemple ci-dessus, le QUAD8 est divisé en 6 TRIA6. On effectue la recherche des facettes de contact sur les sous éléments 2, 3 et 4 qui vérifient \(\mathit{lsn}<0\) . On récupère alors deux SEG3 (en rouge).
Algorithme de récupération des facettes de contact pour les éléments XH:
Boucle sur les \(\mathit{nse}\) sous éléments de l’élément parent XH
Si le sous élément vérifie \(\mathit{lsn}<0\) , alors
Boucle sur les côtés (en 2D) ou les faces (en 3D) du sous-élément
Boucle sur les nœuds sommets du côté ou de la face courante
si la connectivité du nœud est comprise entre 1000 et 1999, c’est un point d’intersection avec la discontinuité, on le sélectionne.
si la connectivité du nœud sommet est strictement inférieure à 1000, alors c’est un nœud sommet de l’élément parent, on le sélectionne uniquement s’il vérifie \(\mathit{lsn}=0\) .
Fin de la boucle sur les nœuds sommets de la face courante
Si tous les nœuds sommets du côté ou de la face courante ont été sélectionnés, il s’agit d’une facette de contact, dont on stocke la connectivité et les coordonnées de ses nœuds
Fin de la boucle sur les côtés (en 2D) ou les faces (en 3D) du sous-élément
Fin de la condition
Fin de la boucle sur les \(\mathit{nse}\) sous éléments de l’élément parent XH
Lorsque l’élément est XHT ou XT, un travail supplémentaire est nécessaire par rapport aux éléments XH. En effet, la découpe en sous-éléments est effectuée conformément à la \(\mathit{lsn}\) , mais pas conformément à la \(\mathit{lst}\) . Les facettes de contact que l’on récupère vérifient donc bien \(\mathit{lsn}=0\) , mais pas nécessairement \(\mathit{lst}<0\) , car les éléments XHT et XT sont susceptibles de contenir le fond de fissure. Un redécoupage conforme à la \(\mathit{lst}\) est alors nécessaire.
Une fois obtenues les facettes de contact telles que \(\mathit{lsn}=0\) , on les redécoupe si nécessaire selon la courbe \(\mathit{lst}=0\) (voir )
Pour ce faire, on utilise exactement la même démarche que lors de la découpe des éléments principaux en sous éléments d’intégration (découpe conforme à la surface \(\mathit{lsn}=0\) ). Dans l’exemple ci-dessous, le SEG2 de droite est découpé suivant la courbe \(\mathit{lst}=0\) .
Fig. 413 Redécoupage des facettes au niveau du fond de fissure#
Pour les éléments multi-fissurés, la procédure est la même. En effet, la découpe en sous éléments est effectuée de manière séquentielle, fissure après fissure, de manière à ce que les sous-éléments d’intégration finaux soient conformes à chacune des fissures (voir ). On utilise alors, pour chaque fissure, le même procédé que pour les éléments XH. Cependant, pour les fissures sur lesquelles une autre fissure se branche, on effectue pas nécessairement la recherche des facettes de contact dans les éléments tels que \(\mathit{lsn}<0\) . Afin de s’assurer que les facettes de contact soient conformes à la jonction, on effectue la recherche des facettes de contact du côté où la fissure se branche (voir ).
Fig. 414 Récupération des facettes de contact pour un élément multi-fissurés. Les facette de la fissure 1 sont récupérées du côté ou la fissure 2 se branche.#
Lorsque plusieurs fissures se branchent sur une même fissure au sein d’un élément, si tous les branchements ne se font pas du même côté, il n’existe pas en général de facettes de contact conformes à chaque fissure(voir ). L’algorithme de découpage émet alors un message d’alarme qui invite l’utilisateur à raffiner le maillage de manière à ne plus avoir de jonctions multiples au sein d’un même élément.
Il n’y a pas de redécoupage nécessaire car les éléments multi-fissurés ne peuvent pas contenir de fond de fissure, ce sont des éléments d’interface uniquement.
Fig. 415 Récupération des facettes de contact impossible pour un élément multi-fissuré dans lequel deux fissures se branchent de part et d’autre d’une première fissure.#
Algorithme de récupération des facettes de contact pour les éléments multi-fissurés:
Boucle sur les \(\mathit{nfiss}\) fissures qui traversent l’élément
\({\mathit{signe}}^{\mathit{ref}}=-1\)
Si une ou plusieurs fissures se branchent sur \(\mathit{ifiss}\) , on détermine le signe correspondant au côté où les fissures se branchent \({\mathit{signe}}^{\mathit{ref}}=\mathit{signe}(\mathit{côté})\) . Si les fissures ne se branchent pas toutes du même côté, on émet le message d’alarme suivant:
«Plusieurs fissures sont branchées de part et d’autre d’une même fissure au sein du même élément. La récupération de facettes de contact conformes à chaque jonction risque d’échouer. Raffinez le maillage afin de que cette situation n’arrive plus»
Boucle sur les \(\mathit{nse}\) sous éléments de l’élément parent XH
Si le sous élément vérifie \({\mathit{signe}}^{\mathit{sous}\mathit{élément}}(\mathit{ifiss})={\mathit{signe}}^{\mathit{ref}}\) , alors
Boucle sur les côtés (en 2D) ou les faces (en 3D) du sous éléments
Boucle sur les nœuds sommets du côté ou de la face courante
si la connectivité du nœud est comprise entre 1000 et 1999, c’est un point d’intersection avec l’une des fissures, on le présélectionne. On le sélectionne définitivement si \(∣{\mathit{lsn}}^{\mathit{ifiss}}∣<\mathit{loncar}\ast {10}^{-8}\) avec \(\mathit{loncar}\) la longueur caractéristique de l’élément parent. Cette partie de la programmation évoluera de façon à se débarrasser de ce critère. En effet, deux nœuds d’un même côté, de connectivités comprises entre 1000 et 1999 réalisent le minimum d’une seule level set qu’il est donc facile d’identifier;
si la connectivité du nœud sommet est strictement inférieure à 1000, alors c’est un nœud sommet de l’élément parent, on le sélectionne uniquement s’il vérifie \({\mathit{lsn}}^{\mathit{ifiss}}=0\) .
Fin de la boucle sur les nœuds sommets de la face courante
Si tous les nœuds sommets du côté ou de la face courante ont été sélectionnés, il s’agit d’une facette de contact pour la fissure \(\mathit{ifiss}\) , dont on stocke la connectivité et les coordonnées de ses nœuds
Fin de la boucle sur les côtés (en 2D) ou les faces (en 3D) du sous éléments
Fin de la condition
Fin de la boucle sur les \(\mathit{nse}\) sous-éléments de l’élément parent XH
Fin de la boucle sur les \(\mathit{nfiss}\) fissures qui traversent l’élément
Dans le cas 3D, pour un hexaèdre, on peut majorer le nombre maximum de facettes de contact pour une seule fissure par 18. On est en effet susceptible d’atteindre ce chiffre pour un élément qui contient le fond de ffissure. Lors du redécoupage des facettes triangulaires conformément à la \(\mathit{lst}\) , certaine facettes donnent naissance à deux facettes (voir ), ce qui explique que l’on puisse avoir au final près de 18 facettes.
Fig. 417 Redécoupage d’une facette triangulaire au niveau du fond de fissure#
Dans le cas multi-fissuré, le nombre maximum de facettes de contact par fissure est plus élevé que le nombre maximum de facettes de contact pour un élément traversé par une seule fissure. En effet, le découpage séquentiel génère plus de sous éléments pour chaque nouvelle fissure et donc plus de facettes de contact car ces dernières s’appuient sur les sous éléments. En 3D, pour un hexaèdre, on fixe le nombre maximal de facettes de contact par fissures à 30. Un élément 3D traversées par 4 fissures peut donc avoir un total de \(30\ast 4=120\) facettes. En 2D, on fixe le nombre maximum de facettes de contact par fissures à 7. Un élément 2D traversé par 4 fissures peut donc avoir un total de \(4\ast 7=28\) facettes.
où \(\left\lbrace N\right\rbrace\) est le vecteur de la base enrichie des fonctions de forme et \(\left\lbrace u\right\rbrace\) le vecteur des ddls nodaux de déplacement.
Sur un élément fini X-FEM occupant un domaine \({O}_{e}\) , l’expression de la matrice de rigidité élémentaire est (d’après l’expression (4903) du PTV):
Notons au passage que les quantités à intégrer ne changent pas, et font toujours référence à l’élément fini X-FEM \({O}_{e}\) (appelé «élément parent») et non pas aux sous-éléments \({O}_{\text{se}}\) .
Sur chaque sous-élément, l’intégrande est continue. En effet, les sous-éléments respectent la discontinuité de la fissure, donc la fonction \({H}_{j}\left(x\right)\) y est constante, et les fonctions \({\phi }_{{i}^{\prime }}\) , \({F}^{\alpha}\) et \({({F}^{\alpha})}^{\prime }\) y sont continues.
Lorsqu’on se place dans le cadre des grandes rotations ou des grandes déformations, il faut prendre en compte la rigidité géométrique. La rigidité géométrique est aussi utile dans le cas de l’analyse modale.
Dans le cas des grandes rotations par exemple, le travail virtuel des forces internes peut être écrit sur la configuration initiale sous la forme:
On cherche les dérivées des fonctions singulières par rapport aux coordonnées globales \((x,y,z)\) . Les fonctions singulières sont données en fonction des coordonnées polaires \((r,\theta )\) .
Dérivées des fonctions singulières dans la base polaire \((r,\theta )\) :
Pour les termes non singuliers (termes classiques et enrichissement Heaviside), les fonctions à intégrer sont des polynômes. L’intégration numérique par un schéma de Gauss-Legrendre sur chaque sous-tétraèdre est alors bien adaptée. Soit \(\text{npg}\) le nombre de points de Gauss. Ces points ont pour coordonnées dans le sous-tétraèdre de référence \(({s}_{g},{t}_{g},{u}_{g})\) et le poids associé est noté \({w}_{g}\) . Soit \(J\) la Jacobienne de la transformation sous-tétraèdre \(\to\) sous-tétraèdre de référence. Son déterminant est constant par sous-élément. La procédure classique s’écrit:
La fonction \(f\) étant seulement connue dans le repère lié à l’élément parent de référence, il est nécessaire d’exprimer les coordonnées du point de Gauss \(({s}_{g},{t}_{g},{u}_{g})\) dans ce repère [bib33]. Tout d’abord, on les transporte dans le repère réel (utilisation simple des fonctions de forme du sous-tétraèdre), puis on les transporte dans le repère lié à l’élément parent de référence (résolution d’équations non-linéaires par l’algorithme de Newton-Raphson):
En effet, il faudrait plutôt écrire\({\int}_{{\Omega}_{\text{se}}^{\text{ref}}}f°{t}_{2}°{t}_{1}({s}_{g},{t}_{g},{u}_{g})\mid \det(J)\mid d{\Omega}_{\text{se}}\) . La transformation\({t}_{1}\)est linéaire. Par contre, la transformation\({t}_{2}\)est non-linéaire. C’est l’inverse d’une transformation linéaire. C’est une combinaison de fractions rationnelles et de racines carrées à plusieurs variables. De ce fait, on ne connaît pas l’effet de\(f°{t}_{2}°{t}_{1}\)surl’ordre maximal des monômes dans la cas général. Par contre, pour un élément réel à bords parallèles (losange en 2D),\({t}_{2}\)étant linéaire, on considère que\(f°{t}_{2}°{t}_{1}({s}_{g},{t}_{g},{u}_{g})\)sont des monômes en\(({s}_{g},{t}_{g},{u}_{g})\)du même ordre que ceux définis sur l’élément parent de référence.
Les fonctions de forme des éléments iso-paramétriques sont des monômes de type \({\xi}^{i}{\eta}^{j}{\zeta}^{k}\) avec \(i+j+k\le p\) (voir Tableau 84). Les dérivées des fonctions de forme sont alors des monômes de même type avec \(i+j+k\le p-1\) . Les quantités à intégrer sont donc des monômes d’ordre \(m={(p-1)}^{2}\) . On se réfère à [bib34] pour déterminer le schéma adapté permettant une intégration exacte [3] de tous les termes. Pour les éléments à bords non parallèles, l’intégration sera alors approchée.
Tableau 84 Schémas d’intégration selon le type de l’élément parent#
Les schémas d’intégration numérique de Gauss-Legendre ont été conçus pour l’intégration de polynômes. Or la présence de monômes en \(1/\sqrt{r}\) provenant de l’enrichissement avec les champs asymptotiques conduit à une convergence très lente de la précision de ces schémas: un nombre considérable de points de Gauss est alors nécessaire pour obtenir une erreur acceptable. Lorsque la singularité se trouve sur la frontière du domaine d’intégration, des changements de variables successifs bien choisis permettent de se ramener à l’intégration sur un domaine unitaire d’une fonction non-singulière, polynomiale (ou quasi-polynomiale), sur laquelle une intégration classique de Gauss-Legendre s’avère efficace et rapidement convergente [bib28]. Une dizaine de points de Gauss dans chaque direction suffit à atteindre les limites de la précision numérique d’un calculateur. Lorsque la singularité ne se trouve pas sur la frontière, il convient de découper le domaine en sous-domaines, de telle sorte que la singularité se trouve sur les frontières des sous-domaines. Malgré des résultats très satisfaisants en terme de convergence de l’erreur relative d’intégration sur chaque terme [bib28], l’effort d’implémentation à fournir pour les structures 3D nous paraît considérable [bib35]. D’autant plus que la différence au niveau des erreurs globales (après résolution) entre une intégration classique avec un nombre de points de Gauss assez élevé mais raisonnable et une intégration singulière modifiée n’est pas vraiment significative
Pour le moment, on a mis en place le même schéma classique de Gauss-Legendre que celui qui est utilisé pour l’intégration des termes non singuliers.
la surface \({\Gamma}^{t}\) discrétisée correspond aux faces des éléments 3D. un sous-découpage possible pour ses faces 2D est celui correspondant à la trace du sous-découpage des éléments 3D. Dans Aster, on choisit de prendre le sous-découpage générique 2D, correspondant à celui décrit à la r7.02.12-sous-decoupage-2d.
De même que pour l’intégration de la rigidité, il faut déterminer les coordonnées des points de Gauss du sous-triangle de référence dans l’élément parent de référence
Soit le point de Gauss de coordonnées \(({s}_{g},{t}_{g})\) dans le sous-triangle de référence. Le problème est que l’environnement de calcul est 3D, On crée donc un repère réel 2D local à chaque sous-triangle Soit le sous-triangle \(\text{ABC}\) , de normale \(n=\frac{\overrightarrow{\text{AB}}\wedge \overrightarrow{\text{AC}}}{\parallel \overrightarrow{\text{AB}}\wedge \overrightarrow{\text{AC}}\parallel }\) . La base locale s’écrit:
où \({\psi}_{i}\) sont les fonctions de formes du sous-triangles et \(({x}^{i},{y}^{i})\) les coordonnées des sommets du sous-triangle dans le repère local 2D \(({x}_{l},{y}_{l})\) .
Comme pour l’intégration des termes de rigidité, il faut maintenant utiliser une méthode de Newton-Raphson pour déterminer les coordonnées du point de Gauss dans le repère de référence de l’élément parent (élément de face). Pour cela, toutes les coordonnées réelles sont exprimées dans le repère local réel 2D.
Intégrande du second membre des efforts surfaciques sur les lèvres de la fissure#
Les seconds membres dus aux efforts surfaciques sur les lèvres de la fissure s’écrivent de manière discrétisée:
L’intégration se fait sur les deux lèvres de la fissure (\({\Gamma}_{c}={\Gamma}^{1}+{\Gamma}^{2}\) ), avec pour convention:
la lèvre \({\Gamma}^{1}\) est telle que la fonction Heaviside est négative et l’angle polaire vaut \(-\pi\) ,
la lèvre \({\Gamma}^{2}\) est telle que la fonction Heaviside est positive et l’angle polaire vaut \(\pi\) ,
Remarque:
Les normales extérieures aux lèvres étant opposées, les efforts surfaciques équivalents à la pression imposée ont des signes opposés.
L’intégration du second membre dû aux efforts appliqués sur les lèvres se fait sur les facettes des lèvres, grâce à la donnée des objets correspondant à la discrétisation des lèvres en facettes triangulaires (voir §4.3 dans r5.03.54). Ces objets ont été créés à l’origine pour la prise en compte du contact.
Pour les éléments 2D, les options CHAR_MECA_PRES_R, CHAR_MECA_PRES_F, CHAR_MECA_PCISA_R et CHAR_MECA_CISA_F permettent d’appliquer des efforts mécaniques de pression répartie et de cisaillement sur les facettes de contact. Seules les options CHAR_MECA_PRES_R et CHAR_MECA_PRES_F sont actuellement programmées pour les éléments 3D, ce qui signifie qu’il n’est pas possible d’appliquer sur les lèvres des efforts de cisaillement.
En multi-fissuration, on peut avoir des efforts surfaciques différents sur différentes fissures.
La figure ci dessous () illustre l’intégration surfacique numérique sur les facettes de contact d’un quadrangle. En chaque point de Gauss de la facette, on détermine un vecteur unitaire normal et un vecteur unitaire tangent uniquement à l’aide de la topologie des facettes de contact.
Fig. 418 Iµntégration surfacique sur les facettes d’un quadrangle#
En 3D, on détermine de même en chaque point de Gauss un vecteur normal unitaire uniquement à partir de la topologie des facettes de contact. On peut éventuellement créer une base orthonormée en chaque point de Gauss en complétant ce vecteur normal par deux vecteurs unitaires tangents à la discontinuité et orthogonaux entre eux.
Fig. 419 Intégration surfacique sur les facettes d’un tétraèdre.#
La méthode X-FEM a été introduite dans l’opérateur de thermique linéaire afin de pouvoir réaliser des calculs chaînés thermo-mécaniques en prenant en compte la fissure dans la résolution thermique.: le champ de température est alors discontinu à travers la fissure si aucun chargement thermique n’est appliqué sur ses lèvres (avec une condition naturelle de flux nul). Il est également possible de définir une condition d’échange thermique entre les lèvres de la fissure par la donnée d’un coefficient d’échange.
Remarque :
La résolution des problèmes de thermique linéaire avec X-FEM s’appuie sur une algorithmie quasi-identique à celle mise en œuvre pour la résolution de problèmes mécanique de structures fissurées avec X-FEM, détaillée dans lar7.02.12-probleme-fissuration. En particulier les procédures permettant:
l’enrichissement de l’approximation du champ des inconnues nodales (statut des nœuds, statut des mailles, annulation des degrés de liberté enrichis à tort);
le sous-découpage des éléments traversés par une fissure;
les procédures d’intégration des matrices et des vecteurs élémentaires sur ce sous-découpage ;
sont identiques à celles qui ont été précédemment détaillées, elle ne le seront donc pas dans cette partie.
L’ensemble des modélisations et des types de chargement disponibles pour l’algorithme de thermique linéaire transitoire dans le cas de modèles classiques (non X-FEM) n’a pas été intégralement répercuté pour les éléments finis thermiques X-FEM. Toutes ces fonctionnalités restent bien entendu disponibles sur la partie non-enrichie du modèle, et sont documentées dans [bib74]. Des évolutions à venir permettront de venir compléter progressivement la liste des fonctionnalités disponibles à ce jour.
Type de problèmes traités
Résolution de transitoires thermiques linéaires (indépendance à la température de paramètres matériaux), pour un matériau isotrope (conductivité scalaire). Le branchement de fissure (multi-Heaviside) est exclu.
Éléments finis disponibles
Pour le moment, seuls les éléments X-FEM des modélisations 3D, axisymétrique, et plane ont été introduits pour les mailles support linéaires suivantes:
modélisation
mailles support éléments principaux
mailles support éléments de bord
3D
TETRA4, HEXA8, PENTA6, PYRA5
TRIA3, QUAD4
plane / axisymétrique
TRIA3, QUAD4
SEG2
Il n’existe donc pas d’éléments thermiques X-FEM lumpés, et le maillage à partir duquel on définit le modèle enrichi est obligatoirement linéaire.
Chargements thermiques disponibles sur la partie enrichie du modèle
Pour le moment, le seul chargement disponible sur la partie enrichie du modèle est la condition d’échange thermique entre les lèvres de la fissure. Aucune autre condition aux limites ne peut être imposée sur des éléments thermiques X-FEM de bord, ou sur des nœuds appartenant à des éléments thermiques X-FEM. De même, aucun terme source ne peut être imposé sur des éléments thermiques X-FEM principaux. L’utilisateur doit donc s’assurer d’appliquer ses chargements suffisamment «loin» de la fissure, c’est à dire sur des éléments non enrichis (ou des nœuds dont le support n’est formé que d’éléments non enrichis).
Afin de ne pas alourdir l’expression de la formulation variationnelle, on considère le problème simplifié d’un solide fissuré, constitué d’un matériau linéaire et isotrope, dont une partie de la frontière (suffisamment éloignée de la fissure) est soumise à une température imposée. Le reste de la frontière (fissure exclue) est soumis à une condition de flux nul. Enfin, on impose une condition d’échange thermique entre les lèvres de la fissure.
On considère donc une fissure \({\Gamma}^{c}\) dans un domaine \(\Omega \subset {ℝ}^{{n}_{\mathit{dim}}}\) , avec \({n}_{\mathit{dim}}=2\text{ou}3\) . On note \(n\) la normale sortante à sa frontière \(\partial \Omega\) , \({\Gamma}_{1}^{c}\) et \({\Gamma}_{2}^{c}\) les lèvres de la fissure de normales sortantes respectives \({n}_{1}\) et \({n}_{2}\) , et \({\Gamma}^{d}\) la partie de \(\partial \Omega\) soumise à une température imposée (voir Figure ). On note \(x\) la variable d’espace, \(t\) le temps, et \(T\) la température. L’intervalle \(\left]{t}_{0},{t}_{f}\right]\) correspond à l’intervalle de temps considéré pour la résolution du problème.
Pour le problème considéré, en notant respectivement \(\lambda\) et \(\rho {C}_{p}\) le coefficient de conductivité thermique et la chaleur volumique à pression constante, l’équation de la chaleur s’écrit:
en notant \({T}_{1}={T}_{\text{|}{\Gamma}_{1}^{c}}\) , \({T}_{2}={T}_{\text{|}{\Gamma}_{2}^{c}}\) , et \(h\) le coefficient d’échange thermique entre les lèvres de la fissure, les conditions aux limites s’écrivent:
On introduit ensuite les espaces fonctionnels suivants:
\(V=\left\lbrace v\in {H}^{1}(\Omega ),{v}_{\text{|}{\Gamma}^{d}}=\stackrel{ˉ}{T}\right\rbrace\) et \({V}_{0}=\left\lbrace v\in {H}^{1}(\Omega ),{v}_{\text{|}{\Gamma}^{d}}=0\right\rbrace\)
La formulation variationnelle associée à (4758) et (4582) s’écrit alors: trouver \(T\in V\) telle que
(4307)#\[\begin{split}\begin{aligned} &\int_{\Omega} \rho C_{p} \frac{\partial T}{\partial t} v \, d\Omega + \int_{\Omega} \lambda \nabla T \cdot \nabla v \, d\Omega \\ &+ \int_{\Gamma_{1}^{c}} h (T_{1} - T_{2}) v \, d\Gamma + \int_{\Gamma_{2}^{c}} h (T_{2} - T_{1}) v \, d\Gamma = 0, \quad \forall v \in V_{0} \end{aligned}\end{split}\]
La discrétisation en temps est réalisée par un \(\theta\) -schéma (voir [bib74]). On se place entre deux piquets de temps successifs \({t}_{n}\) et \({t}_{n+1}\) , on note \(\Delta t={t}_{n+1}-{t}_{n}\) le pas de temps, et \(\theta \in [0,1]\) le paramètre du \(\theta\) -schéma. On introduit alors les notations suivantes:
\({T}^{+}=T(x,{t}_{n+1})\) , et \({T}^{-}=T(x,{t}_{n})\)
\({V}^{+}=\left\lbrace v\in {H}^{1}(\Omega ),{v}_{\text{|}{\Gamma}^{d}}=\stackrel{ˉ}{{T}^{+}}\right\rbrace\) et \({V}^{-}=\left\lbrace v\in {H}^{1}(\Omega ),{v}_{\text{|}{\Gamma}^{d}}=\stackrel{ˉ}{{T}^{-}}\right\rbrace\)
L’application du \(\theta\) -schéma à (4583) conduit à la formulation variationnelle suivante: étant donnée \({T}^{-}\in {V}^{-}\) trouver \({T}^{+}\in {V}^{+}\) telle que
On utilise les notations précédemment adoptées pour la mécanique à la r7.02.12-enrichissement-approximation-deplacement. On rappelle l’approximation éléments finis classique:
Cette expression est composée de 3 termes. Le 1er terme est le terme classique où \({T}_{i}^{C}\) désignent les degrés de liberté classiques, le deuxième terme est le terme enrichi par la fonction de selection de domaine \({H}_{j}\) avec \({T}_{j}^{H}\) les degrés de liberté enrichis correspondants, et le troisième terme est le terme enrichi par la fonction singulière \({F}^{1}\) (enrichissement «Crack-Tip») avec \({T}_{k}^{\mathit{CT}}\) les degrés de liberté enrichis correspondants.
On suppose que le domaine est partionné \(\Omega ` par l'interface de la fissure de la fissure :math:`{\Gamma}^{c}\) en deux domaines \({\Omega}_{-}\) du côté de \({\Gamma}_{1}^{c}\) et \({\Omega}_{+}\) du côté de \({\Gamma}_{2}^{c}\) .
On rappelle ci-dessous l’expression de \({H}_{j}\) :
où \((r,\theta )\) sont les coordonnées polaires dans la base locale au fond de fissure déterminées à partir des level-sets \((\mathit{lsn}(x),\mathit{lst}(x))\) .
Remarque
En mécanique, quatre fonctions déterminées à partir du développement asymptotique du champ de déplacement en fond de fissure sont introduites dans le terme d’enrichissement singulier. En thermique, on introduit uniquement la fonction \({F}^{1}\) qui est la seule fonction (parmi les quatre) discontinue au travers de la fissure. Ce type d’enrichissement, retenu dans [bib75] dans le cas d’une condition de flux nul imposée sur les lèvres de la fissure (fissure adiabatique), est ici introduit pour le cas d’une fissure adiabatique comme pour le cas d’une condition d’échange entre les lèvres de la fissure.
Intégration des matrices et des vecteurs élémentaires#
On détaille dans ce paragraphe l’expression des quantités élémentaires associées à chaque terme de la formulation variationnelle (4584) du problème de thermique linéaire discrétisé en temps.
Soit \(E\) un élément thermique X-FEM. On note \({\lbrace \text{T}\rbrace }_{E}^{T}\) , \({[\text{N}]}_{E}\) et \({[\text{B}]}_{E}\) respectivement le vecteur élémentaire des inconnues nodales (à l’instant \(+\) ), la matrice élémentaire des fonctions de formes et la matrice élémentaire des dérivées des fonctions de formes. A titre d’exemple, on peut considérer un élément de dimension \({n}_{\mathit{dim}}=2\text{ou}3\) , à \(n\) nœuds, et de type «Heaviside / Crack-Tip». \({\lbrace \text{T}\rbrace }_{E}^{T}\) , \({[\text{N}]}_{E}\) et \({[\text{B}]}_{E}\) s’écrivent dans ce cas:
\({\lbrace \text{T}\rbrace }_{E}^{T}=\left[{T}_{1}^{C}\text{}{T}_{1}^{H}\text{}{T}_{1}^{\mathit{CT}}\text{...}{T}_{n}^{C}\text{}{T}_{n}^{H}\text{}{T}_{n}^{\mathit{CT}}\right]\) de taille \(\mathrm{3n}\)
\({[\text{N}]}_{E}=\left[{N}_{1}\text{}{H}_{1}{N}_{1}\text{}{F}^{1}{N}_{1}\text{...}{N}_{n}\text{}{H}_{n}{N}_{n}\text{}{F}^{1}{N}_{n}\right]\) de taille \((1,3n)\)
\({[\text{B}]}_{E}=\left[\nabla {N}_{1}\text{}\nabla ({H}_{1}{N}_{1})\text{}\nabla ({F}^{1}{N}_{1})\text{...}\nabla {N}_{n}\text{}\nabla ({H}_{n}{N}_{n})\text{}\nabla ({F}^{1}{N}_{n})\right]\) de taille \(({n}_{\mathit{dim}},3n)\)
Toutes les quantités volumiques sont intégrées selon la procédure s’appuyant sur le découpage en sous-éléments, et détaillée dans le cas de la mécanique à la page 89. On donne ci-dessous l’expression de ces quantités élémentaires.
Ces quantités correspondent à la condition d’échange thermique entre les lèvres de la fissure, et sont intégrées le long des lèvres de la fissure selon la procédure s’appuyant sur la discrétisation des lèvres en «facettes de contact». Cette procédure est détaillée dans le cas de la mécanique à la page 93.
On introduit tout d’abord les matrices élémentaires \({[\stackrel{̃}{{\text{N}}_{1}}]}_{E}\) et \({[\stackrel{̃}{{\text{N}}_{2}}]}_{E}\) , définies comme les matrices des traces des fonctions de formes respectivement sur \({\Gamma}_{1}^{c}\) et \({\Gamma}_{2}^{c}\) :
D’une matrice à l’autre, les coefficients correspondant à l’enrichissement classique sont égaux, et les coefficients correspondants à l’enrichissement Heaviside et «Crack-Tip» sont reliés par:\({H}_{\text{|}{\Gamma}_{2}^{c}}-{H}_{\text{|}{\Gamma}_{1}^{c}}=2\) , et \({F}_{\text{|}{\Gamma}_{2}^{c}}^{1}-{F}_{\text{|}{\Gamma}_{1}^{c}}^{1}=2{F}_{\text{|}{\Gamma}_{1}^{c}}^{1}\)car l’angle polaire vaut\(\pm \pi\) .
On introduit ensuite les matrices de couplage suivantes:
Les quantités élémentaires sont alors données par les expressions suivantes, avec \({\Gamma}_{E}^{c}\) l’ensemble des facettes constituant la discrétisation de la fissure au sein de l’élément \(E\) :
matrice relative à la condition d’échange (option RIGI_THER_PARO_R ):
Calcul énergétique des Facteurs d’Intensité des Contraintes#
Cette partie rappelle la méthode \(G\) -thêta utilisée dans Code_Aster pour le calcul du taux de restitution d’énergie en mécanique de la rupture linéaire, suite à un calcul éléments finis classiques (en 2D ou 3D). Cette méthode est aussi utilisée pour le calcul des facteurs d’intensité des contraintes, en utilisant la forme bilinéaire de \(g\) mais seulement en 2D. Car pour pouvoir l’utiliser en 3D, il faut connaître la base locale au fond de fissure afin d’exprimer les champs asymptotiques en fonction des coordonnées polaires \((r,\theta )\) .
Ensuite, on présente l’apport des level sets en éléments finis classiques, qui a été de pouvoir construire une base locale au fond de fissure en 3D. Comme on l’a vu à la r7.02.12-base-locale-fissure, les gradients des level sets donnent une base locale au fond de fissure. De plus, les coordonnées polaires dans cette base locale s’expriment facilement en fonction des level sets (voir Fig. 361 ). Ainsi grâce aux level sets, la méthode \(G\) -thêta a été applicable au calcul des facteurs d’intensité des contraintes en 3D avec des éléments finis classiques. Dans ce contexte, la fissure est maillée, un calcul classique est réalisé. En post-traitement, on détermine les deux champs level sets à partir de la fissure maillée et on en déduit une discrétisation du fond de fissure. Le champ thêta est construit à partir de ce fond de fissure discrétisé, et la méthode \(G\) -thêta est appliquée.
Le dernier point est le calcul du taux de restitution d’énergie et des facteurs d’intensité des contraintes, suite à un calcul X-FEM. Dans ce cas, la fissure n’est pas maillée, mais la procédure suivie est quasiment la même: la discrétisation du fond de fissure est déterminée à partir des level sets, un champ thêta est construit et la \(G\) -thêta est appliquée. La seule différence se situe au niveau des intégrations numériques: avec X-FEM, les intégrations numériques nécessitent un sous-découpage des éléments à intégrer. En fait, on utilise le sous-découpage réalisée pour le calcul des matrices tangentes et des seconds membres.
En élasticité linéaire, on définit la densité d’énergie libre \(\Psi\) , en l’absence de déformations et de contraintes initiales par une forme quadratique définie positive des composantes du tenseur des déformations:
Expression lagrangienne du taux de restitution d’énergie#
Par définition, le taux de restitution d’énergie locale \(G\) est défini par l’opposé de la dérivée de l’énergie potentielle par rapport au domaine:
\(G=-\frac{\partial W}{\partial \Omega }\)
La méthode des extensions virtuelles utilisée ici est une méthode lagrangienne de dérivation de l’énergie potentielle. Elle consiste à introduite des transformations:
qui à chaque point matériel \(P\) du domaine de référence \(O\) , associent un point spatial \(M\) du domaine transformé \({\Omega}^{\eta}\) . Ces transformations représentant des propagations de la fissure, ne doivent modifier que la position du fond de fissure \({\Gamma}_{0}\) . Les champs \(\theta\) doivent être tangents à la surface de la fissure.
Soit \(n\) la normale à la surface de la fissure, les champs \(\theta\) doivent vérifier:
où \(G(\theta )\) est défini par l’opposé de la dérivée de l’énergie potentielle \(W(u(\eta ))\) à l’équilibre par rapport à l’évolution initiale du fond de fissure:
Pour dériver l’énergie potentielle par rapport à son support (au voisinage du support de référence), on utilise le théorème de transport de Reynolds:
Pour une transformation \({F}^{\eta}\) de classe \({C}^{1}\) et une fonction \(\varphi\) (tenseur d’ordre 0, 1 ou 2) de classe \({C}^{1}\) , en notant \(\dot{\varphi}\) la dérivée lagrangienne liée à la variation de domaine, on a la relation suivante:
En fait, dans la pratique, la transformation \({F}^{\eta}\) n’est pas toujours de classe \({C}^{1}\) car le champ \(\theta\) est seulement \({C}^{1}\) par morceaux.
Premièrement, utilisons la relation donnant la dérivée lagrangienne d’un champ \(\varphi\) en fonction de sa dérivée eulérienne \(\frac{\partial \varphi }{\partial \eta }\) :
D’après cette relation, la force volumique \(f\) étant supposée indépendante de \(\eta\) , c’est-à-dire étant la restriction à \(O\) de champs définis sur \({ℝ}^{3}\) , on peut écrire que \(\frac{\partial {f}_{i}}{\partial \eta }=0\) . Comme il en est de même pour \(g\) et \(T\) , on a les relations suivantes:
Deuxièmement, utilisons la relation donnant la dérivée lagrangienne du gradient d’un champ en fonction du gradient du champ et du gradient de la dérivée du champ:
On note \(s\) l’abscisse curviligne en fond de fissure. Le calcul de la valeur locale de \(G\) nécessite de résoudre l’équation variationnelle suivantepour plusieurs champs \({\theta}^{i}\) :
Le champ scalaire \(G(s)\) est discrétisé sur une base notée \({({p}_{j}(s))}_{1\le j\le N}\) :
\(G(s)=\sum_{j=1}^{N}{G}_{j}{p}_{j}(s)\)
De même, les champs \({\theta}^{i}\) sont discrétisés sur une base notée \({({q}_{k}(s))}_{1\le k\le M}\) . Soit \({\stackrel{ˉ}{\theta}}^{i}\) la trace du champ \({\theta}^{i}\) sur le fond de fissure \({\Gamma}_{0}\) : \({\stackrel{ˉ}{\theta}}^{i}(s)={\theta \mid }_{{\Gamma}_{0}}(s)\) et soit \({\theta}_{k}^{i}\) les composantes de \({\stackrel{ˉ}{\theta}}^{i}(s)\) dans cette base:
Ce système a une solution si on choisit \(P\) champs \({\theta}^{i}\) indépendants tels que: \(P\ge N\) et si \(M\ge N\) . Il peut comporter plus d’équations que d’inconnues, auquel cas il est résolu au sens des moindres carrés.
Méthode G-thêta pour le calcul de KI, KII et KIII avec les level sets#
Ce paragraphe présente l’apport des level sets pour le calcul des facteurs d’intensité des contraintes en 3D. Il est basé sur la séparation des modes mixtes grâce à la forme bilinéaire de \(g\) . Cette forme bilinéaire est utilisée dans Code_Aster pour le calcul par la méthode \(G\) -thêta de \({K}_{I}\) , \({K}_{\mathit{II}}\) en 2D [bib49].
Le calcul des facteurs d’intensité des contraintes peut se faire également par la méthode thêta, en définissant une forme bilinéaire symétrique \(g\) avec:
Ainsi, soient deux champs de déplacements \(u\) et \(v\) , le calcul de \(g(u,v)\) s’effectue grâce à l’expression de \(G\) donnée par l’équation (4586), dans laquelle les déformations liées à \(u\) et \(v\) dérivent des champs \(u\) et \(v\) , et les contraintes liées à \(u\) et \(v\) sont déduites par la loi de comportement:
avec \({\Psi}^{\text{th}}=\mathrm{3K}\alpha (T-{T}^{\text{ref}})\text{tr}\varepsilon\)
Pour chaque terme \(\mathrm{S1}\) et \(\mathrm{S2}\) , on va d’abord écrire leur expression habituelle, puis ensuite, on passera à la forme bilinéaire.
Écrivons maintenant la forme bilinéaire associée. L’énergie libre dépend alors de deux champs \(u\) et \(v\) . Les coefficients matériaux \(\mathrm{C1}\) , \(\mathrm{C2}\) et \(\mathrm{C3}\) sont évidemment inchangés, et les nouveaux coefficients \(S11\) , \(S12\) et \(S13\) ont pour expressions:
Pour le calcul de \({K}_{i}\) , \(v\) est le champ asymptotique en mode i, obtenu pour \({T}_{v}={T}_{\text{réf}}\) , donc les termes en \(({T}_{v}-{T}_{\text{réf}})\) se simplifient:
Cette écriture permet facilement de passer à la forme bilinéaire associée. Il suffit de remplacer les termes \({u}_{i,j}{u}_{k,l}\) par \(\frac{1}{2}({u}_{i,j}{v}_{k,l}+{u}_{k,l}{v}_{i,j})\) . Si on introduit la notation suivante:
où \(\text{TH1}=1\) en D_PLAN, AXISet en 3Det \(\text{TH1}=\frac{1-2\nu }{1-\nu }\) en C_PLAN.
Pour le calcul de \(\mathit{Ki}\) , \(v\) est le champ asymptotique en mode \(i\) , obtenu pour \({T}_{v}={T}_{\text{réf}}\) , donc les termes en \(({T}_{v}-{T}_{\text{réf}})\) se simplifient:
Le terme supplémentaire dans l’expression de \(G\), dû à l’imposition d’une force surfacique \(g\) sur \({\Gamma}_{c}\) de normale extérieure \(n\) est le suivant:
Pour le calcul de \({K}_{i}\) , \(v\) est le champ asymptotique en mode \(i\) , obtenu pour \({g}_{v}=0\) , donc les termes en \({g}_{v}\) et \(\nabla {g}_{v}\) se simplifient:
Pour le calcul de \({K}_{i}\) , \(v\) est le champ asymptotique en mode i, obtenu pour \({T}_{v}={T}_{\text{réf}}=\text{cste}\) , donc les termes en \(\mathit{TTHE}=\frac{-\partial \Psi }{\partial T}(\nabla T\cdot \theta )=\frac{1}{2}3K\alpha \text{tr}\varepsilon (\frac{\partial T}{\partial x}{\theta}_{x}+\frac{\partial T}{\partial y}{\theta}_{y})\) se simplifient:
La forme bilinéaire symétrique \(g(u,v)\) a comme propriété intéressante celle de séparer les trois modes d’ouverture de la fissure. En effet, dans la forme \(g(u,v)\) , si le champ à gauche est le champ de déplacement solution et si le terme à droite est un champ de déplacement singuliers en fond de fissure \({u}_{I}^{S}\) , \({u}_{\text{II}}^{S}\) ou \({u}_{\text{III}}^{S}\) , les facteurs d’intensité des contraintes s’expriment de la manière suivante:
Certains auteurs utilisent le concept d’intégrales d’interactions, à la place de la forme bilinéaire de \(g\) . Ces deux vocabulaires désignent en fait les même quantités. En effet, en notant \({I}^{u,{u}_{I}^{S}}\) l’intégrale d’interaction entre le champ solution et le champ singulier en mode \(I\) , l’expression de \({K}_{I}\) en terme d’intégrale d’interaction [bib53] est:
Cette expression pour \({K}_{I}\) est identique à celle donnée par l’équation (4590) à un terme multiplicatif près. Si on applique la méthode des extensions virtuelles à la forme bilinéaire de \(g\) , l’expression équivalente de l’équation variationnelle (4586) pour les \({K}_{i}\) locaux est:
On retrouve l’analogue de ces expressions écrites avec des intégrales d’interaction dans [bib54] et [bib53].
On rappelle que le second membre des équations (4593) est la forme bilinéaire de \(g\) (voir (4589)) dans laquelle le terme à gauche \(u\) est le champ déplacement solution (provenant de la résolution du problème par la méthode des éléments finis) et le terme à droite \({u}_{i=I,\text{II},\text{III}}^{S}\) est le champ asymptotique (en mode \(i=I,\text{II},\text{III}\) ) dont l’expression analytique est donnée par l’équation (4904). Ces expressions analytiques sont données dans la base locale au fond de fissure \(({e}_{1},{e}_{2},{e}_{3})\) issue des gradients des level sets , or les dérivées intervenants dans \(g{(u,{u}_{I,\text{II},\text{III}}^{S})}_{\theta}\) sont faites par rapport au repère global \(({E}_{1},{E}_{2},{E}_{3})\) (voir les expressions des termes \(\mathit{S1}\) et \(\mathrm{S2}\) au paragraphe précédent). Il faut donc procéder à un changement de base.
Fig. 421 Coordonnées polaires, base locale et base globale#
Soit \(\text{DPODL}\) la matrice des dérivées des coordonnées polaires dans la base locale:
Soit \(u\) un champ auxiliaire exprimé en fonction des coordonnées polaires \((r,\theta )\) . Les dérivées de \(u\) par rapport aux variables \((r,\theta )\) sont déterminées analytiquement. On les écrit sous forme matricielle:
Il reste maintenant à écrire les dérivées de \(u\) dans la base globale \(({E}_{1},{E}_{2},{E}_{3})\) . Pour cela, on note \(R\) la matrice des vecteurs de la base globale écrits dans la base locale (matrice de passage):
\({R}_{\alpha m}=({e}_{\alpha}\cdot {E}_{m})\)
Soient i et j des indices globaux, la \(j\) ème dérivée de la \(i\) ème composante de \(u\) s’écrit alors:
Dans cette expression, les termes \({u}_{\alpha ,\beta }\) sont les composantes de la matrice \(\text{DUDL}\) . Le terme \({R}_{\mathrm{\alpha i},j}\) apparaît comme la courbure de la base locale. Il peut être négligé si la courbure du fond de fissure est faible. Dans la pratique, sa prise en compte ne change quasiment pas les résultats.
De même que précédemment, on note s l’abscisse curviligne en fond de fissure. Le calcul de la valeur locale de \({K}_{I}\) , \({K}_{\text{II}}\) et \({K}_{\text{III}}\) nécessite de résoudre l’équation variationnelle suivante pour plusieurs champs \({\theta}^{i}\) :
Les champs scalaires \({K}_{\text{I}}(s)\) , \({K}_{\text{II}}(s)\) et \({K}_{\text{III}}(s)\) sont discrétisés sur la base notée \({({p}_{j}(s))}_{1\le j\le N}\) :
Jusqu’à présent, la méthode G-thêta a été présentée dans le cadre de la méthode des éléments finis (MEF). Pour adapter la méthode G-thêta au cadre X-FEM, il y a très peu de modifications à apporter. La seule différence réside au niveau des intégrations numériques, pour calculer les intégrales de domaines (4587). De même que pour l’intégration des matrices tangentes de des seconds membres, on réalise l’intégration sur les sous-éléments.
COLOMBO D., GIGLIO M., “A methodology for automatic crack propagation modelling in planar and shell FE models“, Engineering Fracture Mechanics, vol. 73, pp. 490-504, 2006
[bib2]
DHONDT G., “Automatic 3-D mode I crack propagation calculations with finite elements“, International Journal for Numerical Methods in Engineering, vol. 41, pp. 739-757, 1998
[bib3]
HenshellR.D., ShawK.G., “Crack tip finite elements are unnecessary“, International Journal for Numerical Methods in Engineering, vol. 9, pp. 495-507, 1975
[bib4]
ANDRIER B., GARBAY E., HASNAOUI F., MASSIN P., “Helix-shaped and transverse cracking of rotor shafts based on disk shrunk technology“, Actes de la 7ème Conférence Internationale sur la Rupture et Fatigue biaxial/multiaxial, Berlin, Allemagne, 2004
[bib5]
Réthoré J., “Méthode éléments finis étendus en espace et en temps: application à la propagation dynamique des fissures“, Thèse de Doctorat, Institut National des Sciences Appliquées de Lyon, 2005
[bib6]
BELYTSCHKO T., KRONGAUZ Y., ORGAN D., FLEMING M., “Meshless methods: an overview and recent developments”, Computer Methods in Applied Mechanics and Engineering, vol. 139, pp. 3-47, 1996
[bib7]
Fries T-P., Matthies H-G., “Classification and overview of Meshfree methods”, informatikbericht Nr.: 2003-3, July, 2004
[bib8]
Li S.C., Cheng Y.M., “Enriched meshless manifold method for two-dimensional crack modelling“, Theoretical and Applied Fracture Mechanics, vol. 44, pp. 234–248, 2005
[bib9]
Duarte C.A., Babuška I., Oden J.T.,“Generalized finite element method for three-dimensional structural mechanics problems“, Computers & Structures, vol. 77, pp. 215-232, 2000
[bib10]
Duarte C.A., Hamzeh O.N., Liszka T.J., Tworzydlo W.W.,“A generalized finite element method for the simulation of three dimensional dynamic crack propagation“, Computer Methods in Applied Mechanics and Engineering, vol. 190, pp. 2227-2262, 2001
[bib11]
Belytschko T., BLACK T., “Elastic crack growth in finite elements with minimal remeshing“. International Journal for Numerical Methods in Engineering , Vol.45, Pages 601-620, 1999
Gravouil A., Moës N., Belytschko T.,“Non-planar 3D crack growth by the extended finite element and level sets - Part II: Level set update”, International Journal for Numerical Methods in Engineering, vol. 53, pp. 2569-2586, 2002
OSHER S., SETHIAN J. A.,“Fronts propagations with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations“, Journal of Computational Physics, vol. 79, pp. 12-49, 1988
Belytschko T., Moës N., Usui S., Parimi C.,“Arbitrary discontinuities in finite elements”, International Journal for Numerical Methods in Engineering, vol. 50, pp. 993-1013, 2001
Sukumar N., Chopp D.L., Moës N., BelytschkoT., “Modeling holes and inclusions by level sets in the extended finite-element method”, Computer methods in applied mechanics and engineering, vol. 190, pp. 6183-6200, 2001
Stolarska M., Chopp D.L., Moës N., BelytschkoT., “Modelling crack growth by level sets in the extended finite element method”, International Journal for Numerical Methods in Engineering, vol. 51, pp. 943-960, 2001
Moës N., Gravouil A., Belytschko T., “Non-planar 3D crack growth by the extended finite element and level sets - Part I: Mechanical model”, International Journal for Numerical Methods in Engineering, vol. 53, pp. 2549-2568, 2002
Daux C., Moës N., Dolbow J., Sukumar N., Belytschko T.,“Arbitrary branched and intersectioning cracks with the extended finite element method”, International Journal for Numerical Methods in Engineering, vol. 48, pp. 1171-1760, 2000
Melenk J.M., Babuška I.,“The partition of unity finite element method: Basic theory and applications”, Computer Methods in Applied Mechanics and Engineering, vol. 139, pp. 289-314, 1996
Moës N., Dolbow J., Belytschko T., “A finite element method for crack growth without remeshing”, International Journal for Numerical Methods in Engineering, vol. 46, pp. 131-150, 1999
SROUBOULIS T., BabuškaI., COPPS K., “The design and analysis of the Generalized Finite Element Method”, Computer Methods in Applied Mechanics and Engineering, vol. 181, pp. 43-69, 2000
Stazi F., Budyn é., Chessa J., Belytschko T.,“An extended finite element method with high-order elements for curved cracks”, Computational Mechanics, vol. 31, pp. 38-48, 2003
Laborde P., Pommier J., Renard Y., Salaün M.,“High-order extended finite element method for cracked domains”, International Journal for Numerical Methods in Engineering, vl. 64, pp. 354-381, 2005
GENIAUT S., Convergences en mécanique de la rupture: validation des éléments finis classiques et X-FEM dans Code_Aster , note H-T64-2008-00047-FR, 2008
Sukumar N., Moës N., Moran B., Belytschko T.,“Extended finite element method for three-dimensional crack modelling”, International Journal for Numerical Methods in Engineering, vol. 48, pp. 1549-1570, 2000
Sukumar N., Prévost J.H.,“Modeling quasi-static crack growth with the extended finite element method - Part I: Computer implementation”, International Journal of Solids and Structures, vol. 40, pp. 7513-7537, 2003
Minnebo H., «Une nouvelle approche numérique pour le calcul de la durée de vie de disques de réacteurs», Thèse de Doctorat, École Centrale de Nantes, 2006
[bib36] “Éléments de contacts dérivés d’une formulation hybride continue”, Documentation de Référence de code_aster n° r5.03.52
[bib37] Ben Dhia H., Vautier I.,“Une formulation pour traiter le contact frottement en 3d dans le Code_Aster », rapport de recherche HI-75/99/007/A, Juin 1999, EDF
[bib38] Alart P., Curnier A., “A mixed formulation for frictional contact problems prone to Newton like solution methods”, Computer Methods in Applied Mechanics and Engineering, vol. 92, pp. 353-375, 1991
[bib39] Ji H., Dolbow J.E.,“On strategies for enforcing interfacial constraints and evaluating jumps conditions with the extended finite element method”, International Journal for Numerical Methods in Engineering, vol. 61, pp. 2508-2535, 2004
[bib40] Pellet J., “Dualisation des conditions aux limites”, Documentation de Référence de code_aster n° r3.03.01
[bib41] Ern A., Guermond J.L.,Theory and practice of finite elements, Springer, 2004
[bib42] Brenner S.C., Scott L.R., The mathematical theory of finite element methods, 2nded., Springer, 2002
[bib43] Laursen T.A., Simo J.C.“A continuum element-based formulation for the implicit solution of multi-body, large deformation frictional contact problem”, International Journal for Numerical Methods in Engineering, vol. 36, pp. 3451-3485, 1993
[bib44] Wriggers P., “Finite element algorithms for contact problem”, Arch. Of Comp. Meth. In Eng., vol. 2, pp. 1-49, 1995
[bib45] Curnier A, He, Q.C., Klarbring A.,“Continuum mechanics modelling of large deformation contact with friction”, Contact mechanics, ed. Plenum Press, 1995
[bib46] Pietrzak G.,“continuum mechanics modelling and augmented Lagrangian formulation of large deformation frictional contact problems”, Thèse de doctorat, École Polytechnique Fédérale de Lausanne, 1997
[bib47] Alart P., Barboteu M., «Éléments contact, méthode de Newton généralisée et décomposition de domaine» Problèmes non linéaires Appliqués, École CEA – EDF – INRIA 1999
[bib48]
«Taux de restitution de l’énergie en thermo-élasticité linéaire», Documentation de Référence de code_aster n° r7.02.01
Mialon P.,Calcul de la dérivée d’une grandeur par rapport à un fond de fissure par la méthode thêta, EDF – Bulletin de la Direction des Études et Recherches, Série C, n°3, pp. 1-28, 1988
[bib51] Chapelle D., Bathe K. J., “The Inf-sup test”, Computers & Structures, vol. 47, pp. 537-545, 1993
[bib52] Moës N., Béchet E., Tourbier M., “Imposing essential boundary conditions in the X-FEM”, International Journal for Numerical Methods in Engineering , 2006
Gosz M., Moran B.,“An interaction energy integral domain method for computation of mixed-mode stress intensity factors along non-planar cracks fronts in three dimensions”, Engineering Fracture Mechanics, vol. 69, pp. 299-319, 2002
Gosz M., Dolbow J., Moran B.,“Domain integral formulation for stress intensity factor computation along curved three-dimensional interface cracks”, International Journal of Solids and Structures, vol. 35, n°15, pp. 1763-1783, 1998
[bib55] Erdogan G., Sih G.C., “On the crack extension in plates under plane loading and transverse shear”, Journal of Basic Engineering, vol. 85, pp. 519-27, 1963
[bib56] Barth T.J., SethianJ.A., “Numerical schemes for the Hamilton-Jacobi and level set equations on triangulated domains”, Journal of Computational Physics, vol. 145, pp. 1-40, 1998
[bib57] Crandall M.G., Lions P.L., “Two approximations of solutions of Hamilton-Jacobi equations”, Mathematics of Computation, vol. 43, pp. 1-19, 1984
[bib58] Deconinck H., Struijs R., RoeP.L., “Compact Advection Schemes on Unstructured Grids”, Technical Report, VKI, VKI LS 1993-04, Computational Fluid Dynamics, 1993
[bib59] RoeP.L., “Linear Advection Schemes on Triangular Meshes”, Technical Report CoA 8720, Cranfield Institute of Technology, 1987
[bib61] PENG D., MERRIMAN B., OSHER S., ZHAO H., KANG M., “A PDE-based fast local level set method”, Journal of Computational Physics, vol. 155, pp. 410-438, 1999
[bib62] ADALSTEINSSON D., SETHIAN J.A., “A fast level set method for propagating interfaces”, Journal of Computational Physics, vol. 118, pp. 269-277, 1995
[bib63] PRABEL B., COMBESCURE A., GRAVOUIL A., MARIE S., “Level set X-FEM non matching meshes: application to dynamic crack propagation in elastic-plastic media”, International Journal for Numerical Methods in Engineering, vol. 1, pp. 1-15, 2006
[bib64] DUFLOT M., “A study of the representation of cracks with level sets”, International Journal for numerical methods in engineering, vol. 70, pp. 1261-1302, 2007
[bib65] MESCHKE G., DUMSTORFF P., «Energy-based modelling of cohesive and cohesionless cracks via X-FEM», Computational Methods in Applied Mechanics Engineering, Vol. 196, pp. 2338-2357, 2007
[bib66] COLOMBO, PATRICK MASSIN, «Fast and robust level set update for 3D non-planar X-FEM crack propagation modelling», Computer Methods in Applied Mechanics and Engineering, 200 (2011), 2160-2180
[bib67] SUKUMAR N., CHOPP D.L., BÉCHET E., MOëS N., «Three-dimensional non-planar crack growth by a coupled extended finite element and fast marching method», International Journal Numerical Methods in Engineering, 76 (2008), 727-748
[bib68] SHI J., CHOPP D., LUA J., SUKUMAR N., BELYTSCHKO T., «Abaqus implementation of extended finite element method using a level set representation of three-dimensional fatigue crack growth and life predictions», Engineering Fracture Mechanics, 77 (2010), 2840-2863
[bib69] CITARELLA R., BUCHHOLZ F.-G., Comparison of crack growth simulation by DBEM and FEM for SEN-specimens undergoing torsion or bending loading, Engineering Fracture Mechanics, 75 (2008), 489-509
DAUX C., MOES N., DOLBOW J., SUKUMAR N., BELYTSCHKO T., «Arbitrary branched and intersecting cracks with the extended finite element method», International Journal for Numerical Methods in Engineering, 48 (2000), 1741-1760
[bib71] DANIELE COLOMBO, «An implicit geometrical approach to level sets update for 3D non planar X-FEM crack propagation», submitted to Computer Methods in Applied Mechanics and Engineering
SIAVELIS M., «Modélisation numérique X-FEM de grands glissements avec frottement le longd’un réseau de discontinuités.» Thèse de doctorat, Ecole Centrale de Nantes, 2011
DUFLOT M., «The extended finite element method in thermoelastic fracture mechanics», International Journal for numerical methods in engineering, vol. 74, pp. 827-847, 2008
Hansbo, P. Hansbo. A finite element method for the simulation of strong and weak discontinuities in solid mechanics, Comput. Methods Appl. Mech. Engrg., volume 193, 3523–3540, 2004.
