r7.02.13 Algorithmes de propagation de fissures avec X-FEM#
Résumé:
Ce document présente les algorithmes de propagation de fissures dans le cadre de la méthode X-FEM [R7.02.12]. La fissure est alors représentée par des level-sets. La propagation de la fissure consiste donc à mettre à jour les level-sets compte-tenu de l’avancée du front de fissure.
Calcul du champ de vitesse de propagation#
Les premières étapes consistent à calculer deux champs de vitesses, normale et tangentielle, définissant l’évolution dissociée des deux level-sets. On se place dans le cadre d’une propagation en mode mixte.
Vitesse d’avancée du fond de fissure#
L’avancée du fond de fissure en chacun de ses points, évolue dans la base locale selon une loi de Paris. Sa direction peut être déterminée par plusieurs critères détaillés dans le paragraphe [§ 2.7 ]. On peut généralement écrire:
Soit encore:
Et où \(C\) et \(m\) sont des caractéristiques du comportement en fatigue du matériau, \(G\) est le taux de restitution d’énergie local issu de la méthode \(G\) -Thêta, et :math:`beta ` est l’angle de branchement. Cet angle vaut:
Ou encore:
Avec:
:math:`{x}_{i}=frac{{K}_{i}}{{K}_{I}+left |
{K}_{mathit{II}}right |
+left |
{K}_{mathit{III}}right |
}` et \(p=\frac{\sqrt{\pi}-5\nu }{4}\) |
Champ de vitesse étendu au maillage#
Les vitesses normale et tangente sont connues uniquement pour les points d’intersection entre le fond de la fissure et les faces des éléments du maillage. Ces points ne sont pas des nœuds du maillage. Pour la mise à jour des level-sets après une propagation, la vitesse doit être connue pour chaque nœud du maillage [§ 3.1 ]. Il est donc nécessaire qu’on étende la vitesse à tous les nœuds du maillage à partir des valeurs sur le fond de fissure. Pour cela on projette chaque nœud sur le fond de fissure en direction normale et on assigne au nœud la vitesse du point qui est la projection du nœud.
Figure 2.2-1 : projection du nœud \(M\) sur le fond de fissure La projection est mise en œuvre selon l’algorithme suivant ():
Boucle sur les nœuds \(\text{M}\) du maillage
Initialisation: \({d}_{min}=\text{réel maxi}\)
Boucle sur les segments \(\overrightarrow{\text{IJ}}\) : entre deux points \(I\) et \(J\) consécutifs du fond de fissure
On calcule \(s=\frac{\vec{\text{IJ}}\cdot \vec{\text{IM}}}{{\Vert \vec{\text{IJ}}\Vert }^{2}}=\frac{\Vert \vec{\text{IP}}\Vert \cdot \Vert \vec{\text{IJ}}\Vert }{{\Vert \vec{\text{IJ}}\Vert }^{2}}=\frac{\Vert \vec{\text{IP}}\Vert }{\Vert \vec{\text{IJ}}\Vert }\)
Si \(s<0\) , alors \(s=0\)
Si \(s>1\) , alors \(s=1\)
Si on projette le nœud \(M\) sur \(\overrightarrow{\text{IJ}}\) , on obtient le point \(P\) . La position de ce point est déterminée par le vecteur \(\overrightarrow{\text{IP}}\) : \(\overrightarrow{\text{IP}}=s\cdot \overrightarrow{\text{IJ}}\)
On calcule la distance entre le nœud \(M\) et sa projection \(P\) : \(d=∥\overrightarrow{\text{MP}}∥\)
Si \(d<{d}_{min}\) , alors on stocke \({d}_{min}\) , \({\text{I}}_{min}=\text{I}\) , \({\text{J}}_{min}=\text{J}\) et \({s}_{min}=s\)
Fin boucle
Fin boucle
Le paramètre \({s}_{min}\) donne la position de la projection \(P\) de chaque nœud \(M\) du maillage par rapport à la longueur du vecteur \(\overrightarrow{{\text{I}}_{min}{\text{J}}_{min}}\) du fond de fissure auquel la projection \(P\) appartient. La position de la projection \(P\) peut donc être regardée comme position paramétrique sur le vecteur \(\overrightarrow{{\text{I}}_{min}{\text{J}}_{min}}\) . Cela peut être utilisé pour calculer la valeur du vecteur vitesse \(V\) à assigner au nœud \(M\) .
Figure 2.2-2 : détermination de la vitesse du point \(N\) On ramène ce problème à la détermination du vecteur vitesse \({V}_{P}\) d’un point \(P\) du segment \(\overline{\text{IJ}}\) quand les vecteurs vitesse sont connus aux points \(I\) et \(J\) ().
La méthode utilisée pour le calcul de \({V}_{P}\) est arbitraire. Si on regarde la signification physique de la vitesse, on peut dire qu’elle change graduellement entre les deux vecteurs limites \({V}_{I}\) et \({V}_{J}\) avec la position du point \(P\) . Pour la valeur du \(\Vert {V}_{P}\Vert\) on peut donc utiliser une interpolation linéaire entre les deux valeurs \(∥{V}_{\text{I}}∥\) et \(∥{V}_{\text{J}}∥\) :
Pour la direction de \({V}_{\text{P}}\) , on peut déterminer l’angle entre \({V}_{\text{I}}\) et \({V}_{\text{J}}\) et utiliser encore une fois une interpolation linéaire sur l’angle calculé, selon l’algorithme suivante:
on calcule le vecteur normal au plan formé par les vecteurs \({V}_{\text{I}}\) et \({V}_{\text{J}}\) : \({v}_{3}={V}_{\text{I}}\mathrm{\wedge }{V}_{\text{J}}\)
on se ramène dans une base locale formée par les vecteurs \({V}_{\text{I}}\) et \({v}_{3}\) . On calcule le dernier axe de la base: \({v}_{2}={v}_{3}\mathrm{\wedge }{V}_{\text{I}}\) . Ce vecteur est dans le plan formé par les vecteurs \({V}_{\text{I}}\) et \({V}_{\text{J}}\) .
on calcul les vecteurs unitaire de la base dans ce plan: \({u}_{1}=\frac{{V}_{\text{I}}}{\Vert {V}_{\text{I}}\Vert }\) et \({u}_{2}=\frac{{v}_{2}}{\Vert {v}_{2}\Vert }\)
on calcule l’angle entre les deux vecteurs vitesse \({V}_{\text{I}}\) et \({V}_{\text{J}}\) : \(\alpha =\mathit{acos}\left(\frac{\sum_{i=1}^{3}{V}_{\text{I}i}\cdot {V}_{\text{J}i}}{\Vert {V}_{\text{I}}\Vert \cdot \Vert {V}_{\text{J}}\Vert }\right)\)
on calcule le vecteur unitaire de la vitesse au point \(P\) par interpolation linéaire:
\({u}_{\text{P}}=\cos(\alpha {s}_{min})\cdot {u}_{1}+\sin(\alpha {s}_{min})\cdot {u}_{2}\)
enfin on calcule le vecteur vitesse au point \(P\) en utilisant la valeur \(∥{V}_{\text{P}}∥\) déjà calculée: \({V}_{\text{P}}=∥{V}_{\text{P}}∥\cdot {u}_{\text{P}}\)
Calcul des vitesses pour chaque nœud du maillage#
Pour faire évoluer la level-set normale convenablement [§ 2.5 ] et pour intégrer les deux équations de mise à jour des level-set [§ 3 ], on doit connaître les vecteurs vitesse normale \({V}_{N}\) et tangente \({V}_{T}\) pour chaque nœud du maillage. Pour cela il suffit de calculer les directions normale et tangente des level-sets en utilisant les gradients des level-sets.
Malheureusement cette solution n’est pas applicable à tous les nœuds. En effet, après la mise à jour, la propriété d’orthogonalité des level-sets n’est pas garantie [§ 4.1 ] pour tous les nœuds. Cette propriété est très importante pour garantir l’orthogonalité des vecteurs vitesse \({V}_{N}\) et \({V}_{T}\) . On doit donc utiliser une autre solution.
En regardant la façon utilisée pour étendre le champ de vitesse aux nœuds du maillage [§ 2.2 ,§ 2.5 ], on envisage la solution suivante:
si la level-set tangente au nœud est négative: elle évolue dans une direction parallèle à la surface de la fissure. On utilise donc les gradients des level-sets pour calculer les directions normale et tangente. En effet, cela signifie que le nœud est dans l’espace de la fissure existante et que la valeur de la level-set normale ne change pas pendant la mise à jour des level-sets contrairement à la valeur de la level-set tangente [§ 2.5 ].
si la level-set tangente au nœud est positive, on utilise la base locale de la projection du nœud sur le fond de fissure. En effet, l’évolution des level-sets doit être cohérente avec les modifications du fond de fissure.
Pour le premier cas, le calcul des directions normale et tangente ne pose pas de problème parce que les gradients des level-sets \(\varphi ` aux nœuds sont toujours connus: :math:`{V}_{N}={V}_{N}^{\text{P}}\cdot \nabla {\varphi}_{n}\) et \({V}_{T}={V}_{T}^{\text{P}}\cdot \nabla {\varphi}_{t}\) .Toutefois les composantes normale \({V}_{N}^{\text{P}}\) et tangente \({V}_{T}^{\text{P}}\) de la vitesse \({V}_{\text{P}}\) doivent être calculées par rapport à la base locale du fond de fissure au point \(P\) en utilisant la solution suivante décrite pour le second cas.
Dans ce cas, le calcul est plus difficile parce que la base locale du fond de fissure est connue uniquement pour les points d’intersection entre le fond de la fissure et les faces des éléments du maillage. La situation est visible en figure : on connaît la base locale aux points \(I\) et \(J\) et on veut calculer la base locale au point \(P\) , projection d’un nœud sur le fond de fissure. L’orientation de cette base doit être comprise entre les orientations de la base au point \(I\) et de la base au point \(J\) selon la position du point \(P\) . Cette dernière est bien décrite par le paramètre \({s}_{min}\) déjà calculé [§ 2.2 ].
Figure 2.3-1 : détermination de la base locale au point \(P\)
Si on ne considère pas la translation de l’origine, le changement d’orientation de la base locale peut être vu comme une rotation de la base à partir de l’orientation au point \(I\) , pour laquelle \({s}_{min}=0\) . On a la rotation maximale pour \({s}_{min}=1\) , c’est à dire au point \(J\) . L’orientation de la base au point \(P\) est obtenue en utilisant une rotation inférieure proportionnelle à la valeur de \({s}_{min}\) .
En utilisant le théorème d’Euler, on peut toujours décrire une rotation quelconque comme une seule rotation (angle d’Euler) autour d’un seul axe (axe d’Euler). La direction de ce dernier et la valeur de la rotation changent pour le segment \(\overline{\text{IJ}}\) le long du fond de fissure.
La base locale au point \(P\) peut donc être obtenue par une rotation autour de l’axe d’Euler calculé pour la rotation qui transforme la base au point \(I\) dans la base au point \(J\) . La valeur de l’angle de rotation peut être calculée simplement en multipliant l’angle d’Euler par \({s}_{min}\) . Dans le paragraphe suivant [§ 2.4 ], on décrit les équations qu’on utilise pour ce calcul.
Après la détermination de la base locale \({t}_{\text{P}}\) et \({n}_{\text{P}}\) au point \(P\) , on peut calculer la valeur des composantes de la vitesse \({V}_{\text{P}}\) . Soit:
Enfin on calcule les deux vecteurs \({V}_{N}\) et \({V}_{T}\) :
En résumé, pour chaque nœud \(M\) du maillage on calcule la projection \(P\) sur le fond de fissure. Ensuite on calcule la base locale \({t}_{\text{P}}\) et \({n}_{\text{P}}\) au point \(P\) et les composantes de la vitesse \({V}_{\text{P}}\) (qui ont déjà été calculés [§ 2.2 ]):
Enfin on calcule les vecteurs \({V}_{N}\) et \({V}_{T}\) au nœud \(M\) selon le signe de la level-set tangente:
Où \({\varphi}_{n}\) et \({\varphi}_{t}\) sont respectivement la level-set normale et la level-set tangente.
Calcul de la base locale#
Comme déjà décrit dans le paragraphe précédent, on connaît la base locale aux points \(I\) et \(J\) et on veut calculer la base locale au point \(P\) en utilisant une rotation.
La première étape est le calcul de l’axe et de l’angle d’Euler qui décrivent la rotation qu’on doit appliquer à la base \(({t}_{I},{n}_{I},{b}_{I})\) au point \(I\) pour obtenir la base \(({t}_{J},{n}_{J},{b}_{J})\) au point \(J\) . L’algorithme utilisé est le suivant:
On calcule le troisième axe \(b\) de la base locale aux points \(I\) et \(J\) . En fait, uniquement les vecteurs unitaires \(t\) et \(n\) sont stockés:
On calcule la matrice de rotation d’Euler entre les deux bases:
Chaque colonne de la matrice représente un vecteur de la base au point \(J\) \(({t}_{J},{n}_{J},{b}_{J})\) dans la base locale au point \(I\) .
On calcule l’angle d’Euler en utilisant les éléments de la matrice \(\underline{\underline{{R}_{E}}}\) :
On calcule l’axe d’Euler \(e={\left[{e}_{1}{e}_{2}{e}_{3}\right]}^{\text{T}}\) :
L’angle et l’axe d’Euler calculés sont les mêmes pour tous les points \(P\) sur l’arête définie par les points \(I\) et \(J\) . La deuxième étape est le calcul de la base locale pour un point \(P\) pour lequel la position est déterminée par le paramètre \({s}_{min}\) [§ 2.2 ]. L’algorithme utilisé est le suivant:
On calcule l’angle de rotation effective à utiliser:
On calcule la matrice de la rotation entre la base locale au point \(I\) et la base locale au point \(J\) pour une rotation autour de l’axe d’Euler \(e\) égale à \({\vartheta}_{P}\) :
où la matrice \(\underline{\underline{\epsilon}}\) est la suivante:
On calcule les axes \({t}_{\text{P}}\) et \({n}_{\text{P}}\) de la base locale au point \(P\) dans la base locale au point \(I\) :
Les deux vecteurs sont donc coïncidents avec la première et la deuxième colonne de la matrice \(\underline{\underline{{T}_{P}}}\) . En utilisant l’équation ci-dessus pour le calcul de cette matrice, on peut donc écrire:
Et:
Enfin on peut exprimer les vecteurs \({t}_{\text{P}}^{\text{I}}\) et \({n}_{\text{P}}^{\text{I}}\) dans la base globale du maillage, qui est utilisée pour tous les calculs. Si on dénote la composante \(i\) du vecteur \({t}_{\text{P}({t}_{I},{n}_{I},{b}_{I})}\) par \({t}_{\text{P}({t}_{I},{n}_{I},{b}_{I})}(i)\) et la même composante du vecteur \({n}_{\text{P}({t}_{I},{n}_{I},{b}_{I})}\) par \({n}_{\text{P}({t}_{I},{n}_{I},{b}_{I})}(i)\) , on peut écrire:
Ajustement du champ de vitesse normale#
Afin d’empêcher le déplacement de la fissure préexistante, et de faire évoluer la level-set normale convenablement, il est nécessaire d’ajuster le champ de vitesse normale (bib 1 ). On veut avoir une trajectoire de fissure relativement lissée, telle que:
La valeur du \(\Delta t\) est égale au temps total d’intégration des équations de mise à jour des level-sets [§ 3.1 ]. Ainsi, le fond de fissure avance bien selon \({V}_{T}\) et \({V}_{N}\) calculées précédemment. On fait une approximation linéaire en \(\mathit{lst}\) .
Figure 2.5-1 : Ajustement de la vitesse normale
La formule appliquée sur tous les nœuds du maillage est la suivante:
Où \({V}_{N}\) et \({V}_{T}\) sont les champ de vitesse nodaux, et \(\Delta t\) le pas de temps choisi (voir paragraphe [§ 3 ]).
Prise en compte du déversement#
Dans le cas d’une fissure tridimensionnelle, deux rotations sont possibles lors de la propagation d’une fissure. Soit un point du fond de fissure \((P)\) auquel est associé une base locale ( t , n , b ). [bib 21 ] décrit ainsi la rotation de branchement comme la rotation du vecteur vitesse \(V\) dans le plan \((P,t,n)\) et \(\beta\) son angle associé, et la rotation de déversement comme étant la rotation de la base locale autour du vecteur vitesse, d’angle :math:`gamma ` .
Figure 2.6-1: Définition des angles de branchement et de déversement
La prise en compte du déversement consiste donc à effectuer la rotation des bases locales avant l’étape de mise à jour des level-sets. L’étape de propagation suit donc l’enchaînement suivant:
application des rotations de branchement (\(\beta\) ) et de déversement (\(\gamma\) ) aux bases locales des points du fond de fissure,
projection des points du maillage sur le fond de fissure,
attribution d’une base locale et d’une vitesse \(V\) de propagation à tous les points projetés par interpolation entre les points du fond de fissure,
mise à jour des level-sets dans les bases locale avec la vitesse \(V\) .
Les bases locales ayant subi la rotation avant l’étape de mise à jour des level-sets, le vecteur vitesse \(V\) est donc porté uniquement par le vecteur \(t'\) de la nouvelle base locale. Il n’y a donc plus de composante normale à la vitesse dans la nouvelle base locale.
Remarque :
La rotation de déversement ne peut être activée qu’en 3D et lorsque le mot-clé CRIT_ANGL_BIFURCATION de l’opérateur de propagation PROPA_FISS vaut SITT_MAX_DEVER.
Par défaut la rotation de déversement est forcée à \(0\) .
Détermination des angles de bifurcation#
Dans le cas où le déversement n’est pas pris en compte, l’angle de branchement est classiquement déterminé par le principe du «Maximum Hoop Stress criterion» [bib 4] en fonction des coefficients d’intensité des contraintes \({K}_{I}\) et \({K}_{\mathit{II}}\) :
Remarque :
L’angle de branchement peut être calculé par l’opérateur CALC_G [R7.02.05] [U4.82.03]. Dans ce cas il est imposé à PROPA_FISS dans la même table que G et les SIF.
Il peut aussi être calculé par l’opérateur POST_RUPTURE [U4.82.04] à l’intérieur de l’opérateur PROPA_FISS [U4.82.11]. Cette option permet de choisir d’autres critère sque le «Maximum Hoop Stress criterion».
Dans le cas où le déversement est pris en compte, l’angle de branchement et l’angle de déversement sont tous deux calculés avec le critère du «Maximum Hoop Stress criterion» [bib 21 ] en faisant intervenir \({K}_{I}\) , \({K}_{\mathit{II}}\) et \({K}_{\mathit{III}}\) :
Et pour l’angle de déversement:
Avec:
:math:`{x}_{i}=frac{{K}_{i}}{{K}_{I}+left |
{K}_{mathit{II}}right |
+left |
{K}_{mathit{III}}right |
}` et \(p=\frac{\sqrt{\pi}-5\nu }{4}\) |
Remarque :
L’angle de déversement est calculé par l’opérateur POST_RUPTURE [U4.82.04] à l’intérieur de l’opérateur PROPA_FISS [U4.82.11] lorsque le mot-clé CRIT_ANGL_BIFURCATION de ce dernier est égal à SITT_MAX_DEVER.
Mise à jour des level-sets#
Après la propagation, la géométrie de la fissure change et les level-sets doivent évoluer pour bien décrire cette nouvelle géométrie. Même si on peut recalculer les level-sets en utilisant leur définition de distance signée, en pratique cette solution est à éviter pour des questions de performance.
On peut recalculer plus rapidement les level-sets en utilisant les équations de mise à jour des level-sets qui décrivent leur évolution suite à la propagation de la fissure.
Dans les équations suivantes, on utilisera les symboles \({\varphi}_{n}\) , \({\varphi}_{t}\) et \(\varphi ` pour indiquer respectivement la level-set normale (:math:\)mathit{lsn}` ), la level-set tangente (\(\mathit{lst}\) ) et une level-set en général (\(\mathit{lsn}\) ou \(\mathit{lst}\) ).
Équations de mise à jour#
Les level-sets sont mises à jour sur tous les nœuds selon l’équation de Hamilton-Jacobi suivante (bib 1 ):
Où \(\varphi ` est la level-set, :math:\)nabla varphi ` est son gradient et \({V}_{\varphi}\) la vitesse de propagation de la level-set \(\varphi ` (:math:`{V}_{T}\) ou \({V}_{N}\) ). Cette équation est connue comme équation d’évolution des level-sets. Dans (bib 13 ) et (bib 14 ) il est montré qu’une erreur est présente dans cette formulation de l’équation et que la forme correcte à utiliser est plutôt la suivante:
L’intégration de l’équation pose plusieurs problèmes (bib 14 ):
Le nombre d’itérations d’intégration est important si un schéma explicite est utilisé (par exemple Runge-Kutta) à cause de la limitation sur le pas de temps utilisable à chaque itération (condition CFL). L’ajustement du champ de vitesse normale ([§ 2.5 ]) est à l’origine de la très forte diminution de la valeur maximale du pas de temps utilisable et par conséquent de l’incrément du nombre d’itérations nécessaires. Pour éviter ce problème, un schéma implicite peut être utilisé mais cela complique la programmation. Comme troisième alternative, on peut utiliser une méthode de type «fast marching method» (bib 3 ) qui calcule la solution directement et sans itérations (bib 15 ,bib 16 ). Toutefois des erreurs peuvent être présentes dans la solution calculée par cette dernière méthode aux bords du domaine de calcul (bib 15 ,bib 16 ).
En présence d’une valeur élevée de l’angle de propagation \(\beta\) ([§ 2.1 ]), la solution obtenue par intégration de l’équation des level-set peut présenter des instabilités qui vont perturber la convergence des algorithmes d’intégration au cours des itérations de propagation successives.
Les instabilités sus-citées entraînent une grande sensibilité numérique de la méthode aux erreurs arrondis, ce qui conduit à une variation significative des résultats entre machines. Ce comportement n’est pas souhaitable dans une utilisation industrielle de la méthode.
Une équation alternative de mise à jour est donc proposée dans (bib 14 ). Elle repose sur l’interprétation géométrique de la modification des level-sets conséquente à la propagation de la fissure. Si on regarde l’iso-zéro de la level-set tangente (\(\mathit{lst}=0\) de la figure ), on peut dire que chaque point de cette surface se déplace de la quantité \({V}_{T}\cdot \Delta t\) . Par conséquent, la valeur de la level-set au point diminue de la même quantité. Si on regarde l’iso-zéro de la level-set normale (\(\mathit{lsn}=0\) de la figure ), on peut dire que dans la partie du volume de matériau sain (\(\mathit{lst}>0\) ) la level-set est soumise à une rotation autour du fond de fissure et que chaque point se déplace suivant la direction normale à l’iso-zéro d’une quantité \({V}_{N}\cdot \Delta t\) . C’est la variation linéaire de la vitesse \({V}_{N}\) avec \(\mathit{lst}\) ([§ 2.5 ]) qui permet d’obtenir le bon déplacement en chaque point. En généralisant au cas 3D comme dans (bib 14 ), on peut donc écrire les équations de mise à jour suivantes:
L’avantage de ces équations est évident : il s’agit d’équations algébriques explicites qui ne nécessitent l’utilisation d’aucun schéma itératif d’intégration. La solution peut être calculée directement avec un temps de calcul très faible. De plus, aucun problème d’instabilité n’est présent même pour des angles de propagation \(\beta\) élevés comme montré dans (bib 14 ).
La valeur de \(\Delta t\) à utiliser dans les équations de mise à jour des level-sets peut être facilement calculée à partir de l’avancée maximale \(\Delta {a}_{max}\) de la fissure donnée par l’utilisateur:
Où le vecteur \(V({x}_{i})\) est le vecteur de la vitesse de propagation (\(V={V}_{N}+{V}_{T}\) ) d’un point \({x}_{i}\) du fond de fissure. Il faut observer que, dans ce contexte, on prend comme point du fond de fissure les points physiques à l’intersection entre le fond de fissure et les faces des éléments du maillage, pour lesquels on connaît la vitesse et l’angle de propagation, mais qu’on prend aussi les points géométriques du front qui sont la projection normale des nœuds du maillage sur le front lui-même [§ 2.2 ]. Toutefois la vitesse de propagation est calculée pour ces derniers points en utilisant une interpolation linéaire [§ 2.2 ] et pour chaque arête qui forme le fond de fissure la vitesse maximale est obtenue à l’une de ses deux extrémités. On peut donc limiter la recherche de la vitesse maximale aux points d’intersection entre le fond de fissure et les faces des éléments du maillage.
Remarque :
Dans la description ci-dessus, les symboles \(t\) et \(\Delta t\) ont été utilisés pour indiquer le temps. Toutefois la grandeur décrite par ces symboles est effectivement le nombre de cycles de fatigue qui cause l’avancée de la fissure. En fait, si on regarde l’équation utilisée ci-dessus pour le calcul de \(\Delta t\) , \(\Delta t\) est obtenu par division entre la distance d’avancée de la fissure (par exemple: millimètres) et la vitesse d’avancée. La vitesse calculée en utilisant une loi de propagation (par exemple la loi de Paris) est donnée comme une avancée de la fissure par cycle de fatigue (par exemple: millimètres/cycle). On en déduit donc que \(\Delta t\) représente le nombre de cycles de fatigue que l’on doit appliquer pour avoir l’avancée maximale voulue. Le lien entre le nombre de cycles de fatigue et le temps est donné par la fréquence d’application du chargement.
Cela est vrai pour toutes les équations utilisées et on peut dire que toutes les grandeurs qui se ramènent au temps ( \(t\) , \(\Delta t\) , …) expriment en fait un nombre de cycles de fatigue.
Réinitialisation et réorthogonalisation des level-sets#
Introduction#
Au cours de leurs évolutions, les level-sets doivent impérativement garder une définition relativement proche de fonctions de distance signée.
La première raison en est le calcul des facteurs d’intensité de contraintes en fond de fissure par la méthode \(G\) -Thêta, pour lequel les coordonnées polaires \((r,\theta )\) sont définies dans la base locale comme \(r=\sqrt{{\text{lsn}}^{2}+{\text{lst}}^{2}}\) et \(\theta =\arctan(\text{lsn}/\text{lst})\) ; l’utilisation de distances signées pour le calcul des coordonnées polaires est donc nécessaire.
De plus, lors de la propagation de fissure, l’équation de mise à jour présentée plus haut (paragraphe [§ 3 ]) met en jeu les normes des gradients des level-sets, qui doivent rester suffisamment proches de l’unité, car si \(|\nabla \varphi |\) ne vaut pas \(1\) , la fissure n’avance pas avec la vitesse \({V}_{\varphi}\) définie.
Après la phase de mise à jour, il arrive que les level-sets perdent leur sens de fonction de distance signée, notamment dans le cas d’une propagation hors plan (cas d’une sollicitation où \({K}_{\mathit{II}}\mathrm{\ne }0\) ), qui introduit une vitesse normale (voir Figure 4.1-1 ).
Figure 4.1-1 : Propagation de level-sets et perte de fonction de distance signée et d’orthogonalité
Après propagation il est donc nécessaire de passer par une étape dite de «réinitialisation» dont le but sera de minimiser \(∣∣\nabla \mathrm{lsn}∣-1∣\) et , \(∣∣\nabla \mathrm{lst}∣-1∣\) ainsi qu’une étape de «réorthogonalisation» dont l’objectif sera de minimiser \(\mid \nabla \text{lsn}.\nabla \text{lst}\mid ` . En effet, les level-sets ne sont rigoureusement définies que si, en tout point de l’espace, :math:\)mid nabla text{lsn}mid =1` , \(\mid \nabla \text{lst}\mid =1\) et \(\nabla \text{lsn}.\nabla \text{lst}=0\) .
Cependant, dans le cas d’une fissure non plane, il est impossible d’avoir deux fonctions dont les gradients sont orthogonaux entre eux et de normes unitaires en tout point de l’espace. L’antagonisme entre les deux propriétés apparaît sur la Figure 4.1-2 sur lesquelles sont définies, à partir d’une iso-zéro de \(\mathit{lsn}\) , une fonction \(\mathit{lst}\) soit normée, soit orthogonale.
Figure 4.1-2 : Antagonisme Initialisation - Orthogonalisation
Il est donc nécessaire de privilégier l’une ou l’autre de ces propriétés. Étant donnée l’utilisation principale de ces fonctions pour la méthode \(G\) -Thêta, il apparaît plus important que les gradients soient bien normés, plutôt qu’orthogonaux, afin que les coordonnées polaires \((r,\theta )\) aient un sens. Afin d’avoir une base locale orthonormée, on impose aussi sur le fond de fissure que les gradients soient normés et orthogonaux. L’enchaînement des étapes proposé dans (bib 1 ) est donc:
réinitialisation de \(\mathit{lsn}\) :
L’iso-zéro de \(\mathit{lsn}\) , d’importance majeure pour la définition de fissure, est inchangée;
réorthogonalisation de \(\mathit{lst}\) par rapport à \(\mathit{lsn}\) fixé:
\(\mathit{lst}\) est redressée par rapport à \(\mathit{lsn}\) sans modifier le fond de fissure \((\text{lsn}=0)\cap (\text{lst}=0)\) ;
réinitialisation de \(\mathit{lst}\) :
La level-set \(\mathit{lst}\) est initialisé par rapport à son iso-zéro ayant été orthogonalisée.
Résolution par la méthode de fast marching#
Dans la suite on analysera deux méthodes pour réinitialiser et réorthogonaliser les level-sets. La première méthode (méthode «Upwind») peut être appliquée seulement à des éléments rectangulaires en 2D et hexaèdres en 3D, une extension à tous types de mailles étant possible avec l’utilisation d’une grille auxiliaire. La deuxième méthode (méthode «Simplexe») peut être appliquée seulement aux éléments simplexes (triangles en 2D et tétraèdres en 3D), une découpe des mailles sortant de ce cadre permettant d’utiliser cette méthode avec tous types de mailles.
Méthode « Upwind »#
On a vu précédemment qu’après la phase de propagation d’une level-set, il fallait réinitialiser celle-ci afin de conserver une fonction distance signée. En d’autres termes, si l’on suppose que la fonction \(\varphi (\vec{x},t=0)={\varphi}_{0}(\vec{x})\) vérifie pour chaque nœud du maillage la propriété d’une fonction distance signée donnée par \(\left|\nabla \varphi \right|=1\) , la fonction level-set à \(t>0\) n’a plus cette propriété. Il faut donc mettre en place un processus de réinitialisation permettant à :math:`varphi ` de redevenir une fonction distance signée.
Nous proposons ici une méthode de fast-marching introduite par Sethian [19], qui résout l’équation eikonale de la forme:
:math:`begin{array}{c}left |
nabla varphi (x)right |
=F(x),xin Omega \ varphi (x)=0,xin partial Omega end{array}` |
Avec \(\Omega \in {\mathrm{ℝ}}^{n}\) , \(n\text{}\text{}\in \lbrace 2,3\rbrace\) et \(F\) une fonction à valeurs positives.
Le point qui nous intéresse ici n’est qu’un cas particulier de cette équation: lorsque \(F=1\) . La solution donne ainsi la fonction :math:`varphi ` distance signée.
Le principe de la méthode fast-marching consiste à propager de proche en proche la donnée de la distance des points depuis les valeurs les plus faibles, c’est-à-dire proches de l’interface, aux zones les plus lointaines.
Pour la propagation des level-sets, on utilise les valeurs autour de \(\mathit{ls}=0\) pour propager l’information. Deux appels à la méthode sont nécessaires pour traiter tous les nœuds du maillage: le premier appel traitera les points positifs de la level-set et le second les points négatifs.
Pour cela, la méthode répartit les nœuds du maillage en trois catégories:
Frozen points: région du maillage où les points ont déjà été calculés (cercles noirs);
Narrow band: région proche des points déjà calculés (bande grise);
Far away: ensemble des points restants (cercles blancs).
Figure 4.2-1: répartition des nœuds du maillage La méthode de fast-marching peut se passer de l’étape de réorthogonalisation car une fois la réinitialisation de \(\mathit{lsn}\) effectuée, on peut initialiser les valeurs de lst sur son iso-zéro non réorthogonalisée en prenant les valeurs aux nœuds obtenues par un calcul géométrique. On pourra donc utiliser ces valeurs comme condition initiale de la méthode pour propager cette information exact sur \(\mathit{lst}\) au domaine tout entier.
Discrétisation du problème#
Dans la phase de réinitialisation d’une level-set on cherche à résoudre numériquement un cas particulier de l’équation eikonale:
Si on note \({\Delta}_{x}\) \({\Delta}_{Y}\) et \({\Delta}_{Z}\) les pas d’espace en \(x\) , \(y\) et \(z\) , une discrétisation par différence finie upwind de cette équation est:
Avec:
On cherche les solutions du trinôme du second degré qui approxime l’équation () sous la forme:
avec:
Les deux solutions de () sont:
La distance augmente lorsque nous nous éloignons de l’interface .La solution \({\varphi}_{i,j,k}\) est évaluée en un point plus éloigné de l’interface que les points aux distances \({\varphi}_{1}\) , \({\varphi}_{2}\) et \({\varphi}_{3}\) . Ainsi nous avons nécessairement que \({\varphi}_{i,j,k}>{\varphi}_{1}\) , \({\varphi}_{i,j,k}>{\varphi}_{2}\) et \({\varphi}_{i,j,k}>{\varphi}_{3}\) . Or nous avons :
Cela implique:
Donc la solution () est nécessairement la seule acceptable.
Implémentation de la méthode#
La discrétisation numérique de l’équation () par différence finie upwind n’est applicable que sur des maillages hexaédriques (ou rectangulaire en 2D) dont les arêtes sont bien orientées selon trois directions \(({j}_{1},{j}_{2},{j}_{3})\) formant une base de \({\mathrm{ℝ}}^{3}\) .Comme pour la méthode upwind, on utilise une grille auxiliaire pour étendre la méthode à des maillages irréguliers.
Nous présentons ici uniquement un algorithme pour le cas \(\varphi ` positif (:math:\)mathit{lsn}>0` ou \(\mathit{lst}>0\) ), une simple multiplication par \(-1\) de cette fonction permettant de traiter l’autre côté de la level-set.
L’algorithme mis en place est le suivant:
Initialisation
On repère les nœuds pour lesquels \(\mathit{ls}\ge 0\)
On repère les points proches de \(\mathit{ls}=0\) et \(\mathit{ls}\ge 0\) ,c’est-à-dire les nœuds (narrow band) des mailles coupées par \(\mathit{ls}=0\) , tels que \(\mathit{ls}\ge 0\)
Remarque: é tant donné que l’on souhaite propager le calcul vers les zones lointaines, il faut s’assurer au préalable que ces points fournissent une initialisation consistante du point de vue d’une fonction distance. On effectue un calcul géométrique des level-sets sur les éléments coupés par \(\mathit{ls}=0\) *.* On affecte les autres nœuds positifs de la valeur \(+\infty\) *. Ils appartiendront à l’ensemble* **far,alors que les points négatifs seront inclus dans un troisième ensemble* **known.*
Boucle
Tant qu’un nœud a \(\mathit{ls}\ge 0\) , on localise le nœud ayant la plus petite valeur de :math:`varphi ` dans la narrowband.
Pour chaque nœud autour de ce minimum, on récupère \({\varphi}_{1}\) , \({\varphi}_{2}\) et \({\varphi}_{3}\) comme définis par ().
Figure 4.2-2: voisinage du nœud à calculer On calcule pour chaque nœud autour du minimum sa valeur par la formule explicite ().
Fin de la boucle
Fast marching et méthode géométrique#
Les nœuds problématiques
Pour effectuer les calculs par la méthode fast-marching, on a discrétisé l’équation () par différences finies. Il faut donc pour calculer la valeur d’un nœud avoir l’information nécessaire pour appliquer la formule explicite ().Pour chaque nœud calculé, on doit avoir récupéré au préalable \({\varphi}_{1}\) , \({\varphi}_{2}\) et \({\varphi}_{3}\) : si ces valeurs existent alors le calcul par la formule () est immédiat. Dans le cas contraire, c’est-à-dire si :math:`{varphi}_{1}=+infty ` ou :math:`{varphi}_{2}=+infty ` ou :math:`{varphi}_{3}=+infty ` , on ne peut pas appliquer la formule () car le manque d’information entraînera une mauvaise approximation de l’équation ().
Pour résoudre ce problème, on peut utiliser la méthode géométrique pour le nœud problématique en calculant directement sa valeur. Un simple test dans le code permet ainsi de savoir si le point à calculer est évalué par la méthode fast-marching ou par la méthode géométrique. Pour éviter d’éventuels problèmes, il est judicieux de faire aussi un test sur le signe de la valeur sous la racine dans la formule () .
Ainsi si:
Le calcul par la formule () est licite. Alors que si:
La formule () n’est plus valide et un calcul géométrique sur ce nœud est requis.
Réinitialisation de lst
Une fois la réinitialisation de \(\mathit{lsn}\) effectuée par la méthode de fast marching, on passe à la réinitialisation de \(\mathit{lst}\) sans passer par l’étape de réorthogonalisation. Le principal problème pour cette étape est que l’iso-zéro de \(\mathit{lst}\) n’est pas orthogonale à celle de \(\mathit{lsn}\) . La fast marching n’étant applicable que sous la condition que la distance augmente en s’éloignant de l’interface, on peut démarrer la méthode autour de l’iso-zéro de la \(\mathit{lst}\) non-orthogonalisée en effectuant un calcul géométrique pour ces nœuds afin de récupérer les bonnes valeurs et propager l’information comme défini plus haut. La méthode va ainsi réorthogonaliser la \(\mathit{lst}\) et lui attribuer la propriété de fonction distance signée dans un seul et même calcul. Pour augmenter la précision et la localisation du fond de fissure, on effectuera deux réinitialisations de lst à la suite. La première réinitialisation permet de déplacer l’iso-zéro de sorte que celle-ci devienne orthogonale à \(\mathit{lsn}\) et la deuxième augmentera la précision sur la localisation du fond de fissure. Le problème peut venir de la position de l’iso-zéro de \(\mathit{lst}\) avant la première réinitialisation. Si cette position est trop éloignée de la nouvelle position après réinitialisation, on constate une localisation du fond de fissure moins précise. Avec deux étapes, on s’assure que la localisation de l’iso-zéro de \(\mathit{lst}\) ne viendra pas perturber la solution.
Figure 4.2-3: réinitialisation de \(\mathit{lst}\)
Application de la méthode upwind à un maillage irrégulier#
Pour l’application de la méthode des différences finies, on doit utiliser un maillage très régulier, constitué par des mailles hexaédriques en 3D et des mailles rectangulaires en 2D, dont les arêtes sont bien orientées selon les directions de la base de l’espace \({\mathrm{ℝ}}^{3}\) en 3D et \({\mathrm{ℝ}}^{2}\) en 2D. Cela représente une très forte restriction pour la création du maillage et donc limite excessivement les cas qui peuvent être traités par la méthode.
Une solution à ce problème consiste à utiliser deux maillages différents (bib 12 ): un maillage est utilisé pour la structure physique ( maillage physique ) et un maillage est utilisé pour la représentation des level-sets ( maillage level-sets ). Les restrictions de la méthode de différences finies s’appliquent uniquement à ce dernier maillage. Les valeurs des level-sets sur le maillage physique peuvent être calculées par projection des valeurs connues sur le maillage level-sets. Cela est nécessaire pour la localisation du fond de la fissure dans le maillage physique. En fait, les calculs de la mécanique de la rupture sont effectués uniquement sur ce dernier maillage. On peut donc réaliser la mise à jour, la réinitialisation et la réorthogonalisation des level-sets en utilisant uniquement le maillage level-sets. Après on peut projeter les level-sets sur le maillage physique à partir du maillage level-sets. Enfin on peut calculer les facteurs d’intensité des contraintes et la direction de propagation de la fissure en utilisant uniquement le maillage physique.
La projection permet d’obtenir aussi l’indépendance géométrique entre les deux maillages: le maillage des level-sets n’a pas besoin d’être orienté selon les directions de la base choisie pour le maillage physiqueet il peut uniquement occuper le volume de la zone de propagation de la fissure dans la structure.
Les nœuds du maillage level-set, construit comme décrit ci-dessus, constituent donc un nuage de points alignés selon trois directions orthogonales. La grille formée par ces points peut être utilisée pour l’application de la méthode «upwind». Dans la suite de ce document et dans la syntaxe du Code_Aster , on utilisera l’expression grille auxiliaire pour indiquer l’ensemble de ces points. Également, pour simplifier les notations utilisées et les rendre plus claires, on utilisera tout simplement le terme maillage pour indiquer le maillage physique.
Implémentation de la projection entre la grille auxiliaire et le maillage physique#
Dans Code_Aster la grille auxiliaire n’est pas donnée comme un ensemble de points bien orientés mais comme un maillage et donc quelques étapes intermédiaires sont nécessaires pour l’application de la méthode de différences finies upwind. En effet, pour cette méthode tous les points de la grille, pris coïncidents avec les nœuds du maillage level-set donné, doivent être arrangés selon les directions d’une base locale. On doit donc calculer cette base locale et après on doit ordonner les nœuds et calculer la distance entre deux nœuds successifs.
Le calcul de la base locale est très facile parce qu’on peut utiliser la définition des éléments du maillage de la grille auxiliaire. En effet, les arêtes des éléments sont bien orientées par rapport aux axes d’une base locale et on peut définir cette base en prenant des arêtes qui sont orthogonales entre elles. En plus on ne recherche que la direction des axes de la base locale. On peut donc utiliser les arêtes de n’importe quel élément du maillage. En 3D, le premier élément du maillage de la grille auxiliaire permet de déterminer à partir de l’un de ses nœuds, les trois arêtes qui y aboutissent, et ainsi de calculer les trois vecteurs unitaires selon les trois directions des arêtes sélectionnées. En 2D on fait la même opération avec deux arêtes.
Pour le classement des nœuds selon les directions des axes de la base locale calculée, on peut utiliser encore une fois les arêtes des éléments. Pour chaque nœud du maillage, on peut sélectionner les éléments qui le contiennent. Pour chacun des éléments sélectionnés, on sélectionne les arêtes qui partent du nœud. Enfin pour chacune des arêtes sélectionnées, on ordonne les points selon une des directions de la base locale du maillage en utilisant le produit scalaire entre l’arête et les vecteurs unitaires de la base locale et on peut calculer la distance entre les nœuds de l’arête. En 3D l’algorithme mis en place est le suivant (tous les vecteurs sont représentés dans la base choisie pour représenter les nœuds du maillage physique):
On repère les trois vecteurs unitaires de la base locale: \({j}_{1}\) , \({j}_{2}\) et \({j}_{3}\)
Boucle sur les nœuds \(N\) du maillage
On initialise à zéro les valeurs stockées des six nœuds qui sont connectés (successifs ou précédents) au nœud \(N\) selon les trois vecteurs unitaires de la base locale
On repère les éléments qui incluent le nœud \(N\) dans leur définition
Boucle sur les éléments sélectionnes
On repère les trois arêtes qui partent du nœud \(N\)
Boucle sur les trois arêtes sélectionnées
On calcule le vecteur \(\overrightarrow{\mathit{NA}}\) entre les deux nœuds d’extrémité de l’arête. Ce vecteur est orienté selon la direction qui va du nœud \(N\) à l’autre nœud \(A\) extrémité de l’arête.
On calcule les composantes du vecteur \(\overrightarrow{\mathit{NA}}\) selon la base locale:
\({\mathit{NA}}_{1}=\vec{\mathit{NA}}\cdot {j}_{1}\)
\({\mathit{NA}}_{2}=\overrightarrow{\mathit{NA}}\cdot {j}_{2}\)
\({\mathit{NA}}_{3}=\overrightarrow{\mathit{NA}}\cdot {j}_{3}\)
On calcule la tolérance qui est utilisée pour vérifier si une composante est nulle. Cette valeur est calculée comme un pourcentage \(p\) (fixé à 1% dans Code_Aster ) de la longueur du vecteur:
\({l}_{max}=p\cdot ∣\overrightarrow{\mathit{NA}}∣\)
On repère la composante \(n\) plus grande:
\(\Delta =max(∣{\mathit{NA}}_{1}∣,∣{\mathit{NA}}_{2}∣,∣{\mathit{NA}}_{3}∣,)\)
\(n\in [1,2,3]/\Delta =∣{\mathit{NA}}_{n}∣\)
On vérifie que les deux autres composantes sont nulles (\(<{l}_{max}\) ), sinon on s’arrête avec un message d’erreur. En fait, cela signifie que les arêtes du maillage ne sont pas toutes alignées selon les trois mêmes directions.
Si \({\mathit{NA}}_{n}>0\) , alors on stocke le nœud \(A\) comme le nœud successif du nœud \(N\) selon la direction \({j}_{n}\) . Si \({\mathit{NA}}_{n}<0\) , alors on stocke le nœud \(A\) comme le nœud précédent le nœud \(N\) selon la direction \({j}_{n}\) . On stocke aussi \(\Delta\) comme la valeur de la distance entre les nœuds \(A\) et \(N\) selon la direction \({j}_{n}\) . Si \({\mathit{NA}}_{n}=0\) , alors on s’arrête avec un message d’erreur parce que cela signifie que plusieurs nœuds de l’élément sont coïncidents.
Fin de la boucle
Fin de la boucle
Fin de la boucle
En 2D on peut utiliser le même algorithme en prenant seulement deux arêtes par nœud et donc deux composantes. Il est très important qu’avant le début de la première boucle on initialise à zéro les valeurs stockées aux six nœuds qui sont connectés au nœud \(N\) . En effet, si le numéro stocké d’un de ces nœuds est égal à zéro, on peut en déduire que le nœud \(N\) n’a pas un nœud précédent ou successif selon une direction bien définie.
L’arrangement des nœuds selon les directions des axes de la base locale de la grille doit être effectué une seule fois avant l’utilisation de la méthode upwind.
Dans la description ci-dessus, on n’a rien dit sur les modifications que l’on doit apporter au code avant les phases de réinitialisation et de réorthogonalisation des level-sets. En fait, il n’y a pas de grands changements à faire. Le calcul des facteurs d’intensité des contraintes et de la vitesse d’avancée du fond de fissure ne sont pas affectés [§ 2.1 ]. Les algorithmes utilisés pour les phases d’extension du champ de vitesse au maillage grille [§ 2.2 ] et l’ajustement de la vitesse normale [§ 2.5 ] ne changent pas mais on doit donner en entrée aux utilitaires correspondants les nœuds de la grille auxiliaire et pas ceux du maillage physique. Pour cela, les changements du code à faire peuvent être réalisés très rapidement parce que les valeurs de la vitesse d’avancée du fond de fissure sont indépendantes du maillage physique. En fait, le calcul de \(G\) , \({K}_{I}\) et \({K}_{\mathit{II}}\) est conduit sur le maillage physique mais les valeurs de ces grandeurs sont liées aux points du fond de fissure qui sont stockés de façon indépendante par rapport au maillage utilisé. Cela permet d’effectuer l’extension du champ de vitesse et l’ajustement de la vitesse normale directement sur la grille auxiliaire, par projection sur le fond de fissure discrétisé par segments, aucune projection de champ de vitesse entre le maillage physique et la grille n’étant nécessaire car ils utilisent tous les deux la même discrétisation du fond de fissure. Enfin, pour la phase de mise à jour, rien ne change dans les utilitaires sauf l’information sur le maillage utilisé (physique ou level-set): les données sont les mêmes et sont stockées aux nœuds du maillage.
En présence de maillages constitués d’un grand nombre d’éléments, la projection entre la grille auxiliaire et le maillage physique peut prendre beaucoup de temps de calcul. Si on regarde l’utilisation des level-sets pour le maillage physique, on peut dire qu’il n’est pas nécessaire qu’on connaisse la valeur des level-sets pour tous les nœuds de ce maillage. En effet, dans le maillage physique, on utilise les level-set pour trois calculs:
la représentation du fond du fissure,
le calcul de la base locale en chaque point du fond du fissure,
le calcul énergétique des facteurs d’intensité des contraintes.
Pour le premier et le deuxième point, il suffit que la valeur des level-sets soit connue pour les nœuds des éléments coupés par le fond de fissure. Pour le dernier point il suffit que la valeur des level-sets soit connue pour tous les nœuds qui sont compris dans le volume du maillage utilisé pour le calcul énergétique. Ce volume coïncide avec un tore construit autour du fond de fissure. Le rayon de ce tore est spécifié par l’utilisateur dans l’opérateur énergétique utilisé pour le calcul, par exemple CALC_G.
Les trois calculs ci-dessus définissent trois groupes de nœuds pour lesquels la valeur des level-sets doit être connue. Le groupe défini pour le calcul énergétique (point 3) inclut les autres. Il suffit donc qu’on sélectionne dans chaque maillage (physique et grille auxiliaire) les nœuds qui sont compris dans le tore décrit ci-dessus et qu’on fasse la projection uniquement entre les deux groupes définis par ces nœuds.
On peut sélectionner les nœuds à projeter si la distance entre chaque nœud du maillage et le fond de fissure est connue. Pour la grille auxiliaire cette distance est déjà calculée par la routine qui étend aux nœuds du maillage le champ de vitesse connu en fond de fissure [§ 2.2 ]. Dans ce cas, la sélection des nœuds à projeter est donc rapide. Par contre, pour le maillage physique cette distance n’est pas connue et elle doit être calculée en utilisant le même algorithme que celui décrit au [§ 2.2 ].
La projection peut être mise en œuvre en utilisant l’opérateur PROJ_CHAMP disponible dans Code_Aster.
Comme décrit ci-dessus, on a déjà sélectionné les nœuds du maillage physique pour lesquels la valeur des level-sets doit être connue. La sélection est faite en prenant tous les nœuds qui sont compris dans le tore utilisé pour les calculs énergétiques (point 3 ci-dessus). Cela permet d’avoir immédiatement le domaine sur lequel le champ est projeté avec la forme demandée par PROJ_CHAMP. Par contre, un travail supplémentaire est nécessaire pour définir le domaine du champ à projeter. On peut sélectionner les nœuds de ce domaine (sur le maillage level-sets) en utilisant le tore déjà utilisé pour définir l’autre domaine. Ensuite on peut définir le domaine du champ à projeter en cherchant les éléments qui ont au moins un des nœuds sélectionnés avec la technique précédente.
Le dernier problème à résoudre est le suivant: quelle valeur doit-on prendre pour le rayon du tore utilisé pour la sélection des nœuds à projeter ? On a déjà dit que la valeur utilisée dans l’opérateur du calcul énergétique (point 3 ci-dessus) est celle que l’on doit utiliser. Malheureusement cette valeur n’est pas connue de l’opérateur PROPA_FISS. On peut donc utiliser la valeur du paramètre RAYON de cet opérateur qui définit le tore autour du fond de la fissure utilisé pour le calcul du résidu local. Il doit être toujours plus grand que la valeur du rayon utilisé dans l’opérateur du calcul énergétique parce que, après la propagation, les propriétés d’orthogonalité et de distance signée des level-sets sont garanties uniquement dans ce tore. Toutefois la distance à utiliser pour la sélection des nœuds doit être calculée entre chaque nœud et la nouvelle position du fond de la fissure après la propagation. Cette position est définie par les nouvelles valeurs des level-sets après la propagation qui sont connues uniquement pour la grille auxiliaire. La nouvelle position du fond sera connue pour le maillage physique après la projection. La sélection des nœuds intéressés par la projection doit donc forcément être mise en œuvre à partir de la position du fond de la fissure avant la propagation. Le tore à utiliser doit donc inclure toutes les possibles nouvelles positions du fond de fissure et tous les nœuds du tore utilisé pour le calcul du résidu local.
Figure 4.2-4 : sélection des nœuds à projeter
On peut regarder l’exemple donné figure pour une fissure en 2D. L’avancement maximal de la fissure est choisi par l’utilisateur en utilisant le paramètre DA_MAX de l’opérateur PROPA_FISS. La position du fond de la fissure après la propagation est donc comprise dans le cercle de centre coïncidant avec la position du fond avant la propagation et de rayon égal à la valeur de DA_MAX (cercle 1 de la figure). Les nœuds utilisés pour le calcul du résidu local sont compris dans le cercle de centre coïncidant avec la position du fond après la propagation et de rayon égal à la valeur du paramètre RAYON de l’opérateur PROPA_FISS (cercle 2 de la figure). Pour la projection on peut donc sélectionner tous les nœuds compris dans le cercle de centre coïncidant avec la position du fond avant la propagation et de rayon égal à l’addition du DA_MAX et RAYON (cercle 3 de la figure).
Il faut observer enfin que la projection proposée ci-dessus calcule la nouvelle valeur des level-sets uniquement pour les nœuds qui ont été projetés. Toutefois, pour le bon fonctionnement des autres routines X-FEM de Code_Aster , il est nécessaire que les level-sets du maillage physique soient définies pour tous les nœuds. Pour cela on peut simplement préserver la valeur des level-sets pour tous les nœuds qui ne sont pas affectés par la projection. À chaque fois, on peut donc prendre les level-sets du maillage physique avant la propagation et on peut changer la valeur des level-sets uniquement pour les nœuds affectés par la projection. Cette solution est très simple à développer dans le code et est très rapide. En plus elle garantit que la surface existante de la fissure et le fond de la fissure sont toujours correctement représentés par les level-sets parce que les résultats de chaque projection sont préservés en arrière du fond de fissure. Un exemple est donné figure . À gauche de la figure on voit la configuration des level-sets avant la mise à jour. La surface de la fissure est représentée en orange et le fond est marqué avec un petit orange. La configuration des level-sets après la propagation (et la projection) est visible à droite. Seulement les level-sets aux nœuds dans le cercle 3 ont été remplacées par les nouvelles valeurs calculées sur la grille auxiliaire. La nouvelle surface de la fissure est représentée en rouge et la nouvelle position du fond est représentée par un petit rouge. La position du fond avant la propagation est toujours marquée par un point orange. On peut voir que la surface de la fissure qui existait avant la propagation est préservée par la projection (ligne rouge à droite du point orange), ainsi que la continuité de la surface de la fissure entre le domaine de projection (cercle 3) et la partie restante du maillage physique (la ligne rouge est la continuation de la ligne orange). Avec définition du domaine de projection utilisée, les level-sets autour de la position du fond de fissure avant la propagation (point orange) sont toujours affectées par la projection et il n’y a donc pas de risques que de faux fonds de fissure soient détectés (\(\mathit{lsn}=\mathit{lst}=0\) ) après la projection. On peut aussi noter que les level-sets sont discontinues pour le maillage physique mais que cela n’est pas important. En effet, les propriétés d’orthogonalité et de distance signée des level-sets sont garanties pour tous les nœuds impliqués dans les calculs sur le maillage physique.
Figure 4.2-5 : configuration des level-sets sur le maillage physique avant (à gauche) et après (à droite) la propagation. La zone où la level-set tangente est négative est t ei ntée en gris. La surface et le fond de la fissure sont représentés en orange avant la propagation et en rouge après.
Il faut remarquer que les level-sets sont mises à jour sur tous les nœuds de la grille auxiliaire et que pour ce maillage leur représentation est correcte partout. Donc pour la grille auxiliaire on n’a pas les discontinuités qu’on trouve sur le maillage physique après la projection (figure ), sauf dans le cas d’utilisation de la localisation du domaine de calcul décrite au paragraphe [§ 5 ] , ce qui ne représente pas un problème.
Méthode « Simplexe »#
L’utilisation de la grille auxiliaire n’est pas toujours commode pour l’utilisateur, surtout pour des structures complexes. Dans cette optique, une méthode basée sur la même idée que la méthode upwind a été développée. Toujours avec un algorithme de fast marching, nous allons simplement modifier le calcul au niveau de l’élément et travailler directement sur des éléments simplexes [:ref:` 22 < 22 >`] (triangle, tétraèdre).
Figure 4.2-6: Méthode simplexe sur un triangle
Pour simplifier les choses, prenons l’exemple d’un cas 2D. On considère une maille triangle avec \({x}_{1}\) , \({x}_{2}\) et \({x}_{3}\) , les coordonnées des trois nœuds du triangle (voir ). En supposant le fond de fissure linéaire par morceaux, on peut projeter chacun des points sur celui-ci. \(n\) est donc la normale, inconnue et constante par élément.
En gardant à l’esprit que les trois points se projettent sur une droite de normale \(n\) , on a \({d}_{1}\) , \({d}_{2}\) et \({d}_{3}\) qui représentent les valeurs de la level-set sur chacun des trois nœuds. \({d}_{1}\) et \({d}_{2}\) étant connues, on cherche donc un moyen de calculer \({d}_{3}\) . L’idée utilisée revient à écrire l’équation de la projection d’un point sur une droite pour \({x}_{1}\) , \({x}_{2}\) et \({x}_{3}\) . On a:
Ou encore:
On a donc un système de deux équations et trois inconnues. Pour obtenir une solution, on va mettre le système sous forme vectorielle:
Avec:
\(V\) est inversible car les mailles sont non dégénérées, ce qui permet d’écrire l’inconnue \(n\) sous la forme suivante:
\(n\) étant un vecteur unitaire, on peut éliminer cette équation de notre système en écrivant que \(1={n}^{t}n\) . En utilisant l’expression (), on se retrouve avec une équation quadratique où \({d}_{3}\) est la solution:
On se retrouve avec deux solutions pour l’équation quadratique. Pour une solution consistante, il faut choisir la plus grande racine de l’équation.
Choix de la solution de l’équation quadratique
La résolution du système () conduit à deux solutions:
Ces deux solutions sont valables, mais pour notre problème, on choisit de garder () car c’est la seule qui respecte la condition d’obtenir des valeurs de level-set croissante \({d}_{3}>{d}_{2}\) ou \({d}_{3}>{d}_{1}\) .
On observe plusieurs configurations possibles dans notre problème. Examinons les toutes pour se convaincre du choix de cette solution.
Figure 4.2-7: Première configuration
Dans la première configuration (voir figure () ), on connaît la valeur des level-sets aux points \(a\) et \(b\) et on cherche la valeur au point \(c\) . La solution \({x}_{1}\) de () donne une normale :math:`` et la solution en \({x}_{2}\) de () donne une normale \(\vec{{n}_{2}}(0;-1)\) .
La solution () est donc celle qui donne la solution exacte de la valeur de la level-set au point \(c\) , et qui respecte la condition \(\mathit{ls}(c)>\mathit{ls}(b)\) ou \(\mathit{ls}(c)>\mathit{ls}(a)\) . La valeur obtenue avec la solution () dans ce cas-là donne \(\mathit{ls}(c)<\mathit{ls}(b)\) et \(\mathit{ls}(c)<\mathit{ls}(a)\) .
Le domaine de validité de la solution () pour notre problème peut être identifié. On regarde le cas ou l’équation quadratique n’a qu’une racine (racine double), voir figure ().
Figure 4.2-8: Deuxième configuration
On voit les deux configurations d’une maille triangle qui engendrent une solution double pour le problème (discriminant égal à \(0\) ). On peut donc identifier les domaines où la solution () et la solution () sont valables sur la figure ().
Figure 4.2-9: Troisième configuration
La solution () est exacte pour notre problème quand \(\vec{\mathit{ab}}\cdot \vec{\mathit{MN}}\ge 0\) , et la solution () quand \(\vec{\mathit{ab}}\cdot \vec{\mathit{MN}}\le 0\) .
L’idée de la méthode étant de partir des valeurs les plus faibles de level-set, proches de l’iso-zéro et de propager ces valeurs, on voit que la solution () est la mieux adaptée. On choisit donc de prendre la plus grande racine de l’équation quadratique dans notre processus de réinitialisation des level-sets, quelle que soit la configuration. La force de la fast marching est de calculer plusieurs fois la valeur de la level-set pour un même nœud dans différentes configurations autour d’un nœud et de garder la meilleure valeur. Cela corrige donc les cas pathologiques où l’on avait pris la solution () alors que c’était () qui calculait la solution exacte (i.e. il existe au moins une maille contenant le nœud en question pour laquelle la solution () est valide).
Pour garantir que la solution est monotone croissante, on peut ajouter un critère au choix de la solution.
Figure 4.2-10: Quatrième configuration
Soient \({d}_{1}\) , \({d}_{2}\) et \({d}_{3}\) les valeurs de la level-set sur une maille triangle quelconque. On souhaite donc que si \({d}_{1}\) et \({d}_{2}\) augmentent alors \({d}_{3}\) augmente aussi. L’idée est de considérer un gradient de level-set:
Comme la solution est monotone croissante, on veut \({({\nabla}_{d}{d}_{3})}_{i}>0,i\in ⟦1;2⟧\) . En dérivant l’équation quadratique (), on a:
D’où:
Pour que \({({\nabla}_{d}{d}_{3})}_{i}>0\) , il faut que:
On a vu précédemment que si \({d}_{3}>{d}_{2}\) et \({d}_{3}>{d}_{1}\) , on a \({({V}^{t}n)}_{i}<0\) .Comme \(Q\) est symétrique définie positive, la seule possibilité pour que \({({\nabla}_{d}{d}_{3})}_{i}>0\) est que \({({\mathit{QV}}^{t}n)}_{i}<0,i\in ⟦1;2⟧\) . En 2D, cela signifie que \(\vec{n}\) doit venir de l’intérieur du triangle.
Preuve
On a \({\mathit{QV}}^{t}V=I\) et on pose \({\mathit{QV}}^{t}=B=\left(\begin{array}{c}\vec{u}\\ \vec{v}\end{array}\right)\) .On a donc \(\mathit{BV}=I\iff \lbrace \begin{array}{c}\vec{u}\cdot \vec{\mathit{ca}}=1\\ \vec{u}\cdot \vec{\mathit{cb}}=0\\ \vec{v}\cdot \vec{\mathit{ca}}=0\\ \vec{v}\cdot \vec{\mathit{cb}}=1\end{array}\)
Pour satisfaire au problème, on pose \(\parallel \vec{u}\parallel =\frac{1}{\parallel \vec{\mathit{ca}}\parallel \cos(\alpha )}\) et \(\parallel \vec{v}\parallel =\frac{1}{\parallel \vec{\mathit{cb}}\parallel \cos(\beta )}\) où \(\cos(\alpha )\mathit{et}\cos(\beta )\ne 0\) .On veut \({(\mathit{Bn})}_{i}<0\Rightarrow \lbrace \begin{array}{c}\vec{u}\cdot \vec{n}<0\\ \vec{u}\cdot \vec{n}<0\end{array}\Rightarrow \vec{n}\) vient de l’intérieur du triangle.
Algorithme de la méthode simplexe#
Le parcours des mailles dans la structure reste similaire à celui de la méthode upwind. On sélectionne le nœud où la level-set atteint son minimum et on explore les triangles autour de celui-ci pour calculer les valeurs inconnues.
Pour chaque triangle du maillage, on calcule l’équation quadratique et on garde la plus grande racine
La direction de propagation \(n={V}^{-t}(d-{d}_{3}I)\) , est utilisée pour vérifier un critère de monotonie.
En remplaçant \(n\) dans , on a:
Ce critère signifie que \(n\) doit venir de l’intérieur de triangle. On assure ainsi la monotonie de la solution.
Si cette condition n’est pas respectée, on prend la valeur minimale d’un coté du triangle.
Remarque: La généralisation en 3D est très naturelle, les formules étant les mêmes. La seule différence vient au moment où la condition \({(Q(d-{d}_{3}I))}_{i}<0,i\in ⟦1,3⟧\) n’est plus respectée: on calcule alors la plus petite distance sur chaque face du tétraèdre et on garde la plus petite.
Contrairement à la méthode upwind, on n’a plus besoin de la méthode géométrique pour les points du bord. La méthode géométrique est utilisée seulement pour obtenir les valeurs d’initialisation de la level-set tangente, toujours dans l’idée que les level-sets doivent rester orthogonales entre elles.
Calcul des level-sets à proximité des iso-zéros#
Les méthodes présentées précédemment ne donnent de résultats exacts que si les iso-zéros des level-sets passent effectivement par des nœuds. Lorsque l’iso-zéro passe entre deux nœuds, elle a tendance à être légèrement déformée lors de la réinitialisation. Un algorithme de calcul des iso-zéros est donc mis en place. Celui-ci calcule géométriquement les distances signées sur les nœuds des mailles coupées ou tangentées par l’iso-zéro de la level-set et fixe la valeur calculée pour la réinitialisation.
Concernant la réorthogonalisation, les valeurs de \(\mathit{lst}\) situées sur la surface \((\mathit{lsn}=0)\) ne doivent pas être modifiées; on profitera donc de l’algorithme de calcul de \(\mathit{lsn}\) pour calculer les valeurs de lst orthogonales sur les éléments coupés ou tangentés par \((\mathit{lsn}=0)\) et figer ces valeurs pendant la phase de réorthogonalisation. Pour cela, la valeur de \(\mathit{lst}\) est interpolée sur les triangles \(\mathit{IJK}\) que l’on peut former avec les points d’intersection \((\text{lsn}=0)\cap \text{arêtes}\) . Si pour un nœud \(M\) , la projection orthogonale de \(M\) est à l’intérieur de \(\mathit{IJK}\) , on interpole la valeur de \(\text{lst}(M)\) comme \(\mathit{lst}\) du point projeté. Sinon, \(\text{lst}(M)\) est calculée comme \(\mathit{lst}\) du projeté ramené à l’intérieur du triangle \(\mathit{IJK}\) ( Figure 4.3-1 ).
Figure 4.3-1 : Calcul de distance et interpolation de lst sur les éléments coupés par \(\mathit{lsn}=0\) .
En 3D, l’algorithme mis en place après la propagation de \(\mathit{lsn}\) et avant sa réinitialisation est le suivant:
Boucle sur les mailles
Si au moins trois nœuds de la maille sont sur \(\mathit{lsn}=0\) , on la repère comme maille coupée
Si la maille a une arête \([\mathit{AB}]\) telle que \(\text{lsn}(A).\text{lsn}(B)<0\) , on repère la maille coupée
Fin de boucle
Boucle sur les nœuds
Si le nœud appartient à une maille coupée, il sera à calculer
Fin de boucle
Boucle sur les nœuds \(P\) à calculer
Boucle sur les mailles coupées dont le nœud est sommet
Boucle sur les arêtes \([\mathit{AB}]\) de la maille
Si \(\text{lsn}(A).\text{lsn}(B)\le 0\) , on interpole le point \({I}_{\text{AB}}\) :
\(\overrightarrow{{\text{AI}}_{\text{AB}}}=s.\overrightarrow{\text{AB}}\) avec \(s=\frac{\mid \text{lsn}(A)\mid }{\mid \text{lsn}(A)\mid +\mid \text{lsn}(B)\mid }\)
et \(\text{lst}({I}_{\text{AB}})=\text{lst}(A)+s.(\text{lst}(B)-\text{lst}(A))\)
Fin de boucle
Boucle sur tous les triangles \(\mathit{IJK}\) que l’on peut former avec les points \({I}_{\text{AB}}\)
On calcule le projeté \(M\) de \(P\) sur le triangle \(\mathit{IJK}\)
On repère si le point projeté \(M\) est à l’intérieur du triangle
On stocke la distance \(d=\mathit{PM}\)
On interpole linéairement \(\text{lst}(M)\) à partir de \(\text{lst}(I)\) , \(\text{lst}(J)\) et \(\text{lst}(K)\)
Fin de boucle
Fin de boucle
Parmi les distances \(\mathit{PM}\) calculées on recherche la plus petite distance avec \(M\) à l’intérieur de \(\mathit{IJK}\) . Si aucune projection à l’intérieur de \(\mathit{IJK}\) n’a été trouvée, on prend juste la plus petite distance calculée et la valeur de \(\mathit{lst}\) correspondante.
On stocke les valeurs de \(d\) et \(\mathit{lst}\) correspondantes pour le nœud \(P\) .
Fin de boucle
On remplace \(\mathit{lsn}\) et \(\mathit{lst}\) sur les nœuds à calculer par les nouvelles valeurs stockées de d et lst .
Les nœuds repérés comme «nœuds calculés» ont des level-sets qui sont alors figées pendant les phases de réinitialisation de \(\mathit{lsn}\) et de réorthogonalisation de \(\mathit{lst}\) .
L’algorithme lancé après les phases de propagation et de réorthogonalisation de \(\mathit{lst}\) et avant sa réinitialisation est le même, mais il calcule uniquement les distances \(\mathit{lst}\) , de la même manière que \(\mathit{lsn}\) précédemment, sans toucher aux valeurs de \(\mathit{lsn}\) .
L’algorithme décrit ci-dessus est valide seulement en 3D car en 2D on ne peut pas former des triangles. En effet deux points d’intersection seulement peuvent être détectés pour chaque élément.
Pour utiliser en 2D le même algorithme que pour le cas 3D, on peut utiliser la solution suivante. On forme un triangle utilisant les deux points d’intersection (\(A\) et \(B\) ) et un point virtuel (\(C\) ). Celui-ci est choisi de façon à ce que le triangle formé soit dans un plan orthogonal au plan du modèle (plan \(X-Y\) ). Pour éviter des problèmes numériques, les coordonnés du point virtuel (\(C\) ) sont calculées de la manière suivante:
\(X,Y\) = coordonnées du point milieu de l’arête reliant les deux points d’intersection (\(A\) et \(B\) ),
\(Z\) = distance entre les deux points d’intersection (\(A\) et \(B\) ).
Localisation du domaine de calcul des level-sets#
Au paragraphe [§ 4.2.1.5 ] on a dit que, dans le cas de l’utilisation d’une grille auxiliaire pour la méthode upwind, les valeurs des level-sets sont mises à jour, réinitialisées et réorthogonalisées pour tous les points de la grille mais que seulement les valeurs dans une région autour du fond de la fissure sont projetées sur le maillage de la structure. Cela montre que la mise à jour des level-sets, leur réinitialisation et leur réorthogonalisation sont nécessaires uniquement dans la région autour du fond de la fissure, indépendamment de l’utilisation ou non d’une grille auxiliaire, et que le domaine déjà défini pour la projection entre grille auxiliaire et maillage de la structure ([§ 4.2.1.5 ]) peut aussi être utilisé pour localiser le domaine de calcul. Un avantage de cette localisation est évidemment que le temps de calcul nécessaire pour la mise à jour et la réinitialisation des level-sets est inférieur par rapport au cas où tous les nœuds et éléments du maillage ou de la grille auxiliaire sont utilisés. En fait, dans le cas où les level-sets sont mises à jour en des nœuds/points qui sont très loin du fond de la fissure en présence d’oscillations sur l’angle et la vitesse de propagation entre des points successifs du fond de la fissure, les iso-zéros de la level-set normale peuvent présenter des irrégularités très marquées à cause de la variation linéaire imposée à la vitesse normale ([§ 2.5 ]). Un exemple illustratif est décrit figure et figure . On y montre une fissure 3D se propageant en mode mixte selon un angle de propagation qui change de façon linéaire le long du fond de la fissure. Des petites perturbations ont été superposées à la valeur théorique attendue, ce qui fait que la variation de l’angle de propagation avec la position sur le fond n’est pas parfaitement linéaire (figure ). Ces perturbations simulent la présence de petites erreurs dans un calcul numérique de valeurs d’angle.
Figure 5-1 : configuration initiale (à gauche) de la simulation d’une propagation hors plan d’une fissure 3D selon un angle de propagation qui change de façon presque linéaire par rapport à la position considérée sur le fond de la fissure (graphe à droite).
Même si les perturbations de l’angle de propagation, et donc de la vitesse normale, sont petites, elles sont amplifiées par la variation linéaire imposée à cette composante de la vitesse dans le domaine de calcul, ce qui est très clair si on visualise l’iso-zéro de la level-set normale après la seule phase de mise à jour (figure ). La localisation du domaine de calcul permet de limiter cet effet (figure ).
Figure 5-2 : iso-zéro de la level-set normale après la seule phase de mise à jour. L’angle de propagation de la figure a été imposé. Le plan de la fissure initiale a été superposé pour comparaison.
Figure 5-3 : iso-zéro de la level-set normale obtenue avec une restriction du domaine de mise à jour à comparer à celle de la figure *. Il faut remarquer que pour cette figure le plan de la fissure initiale n’a pas été représenté, au contraire de ce qui a été fait pour la figure* . On peut donc en conclure que la valeur de la level-set n’a pas été changée pour les points en dehors du domaine de calcul.
Sélection des nœuds/points et des éléments du domaine de calcul localisé#
Le domaine à utiliser pour la localisation est celui déjà défini pour la projection entre la grille auxiliaire et le maillage de la structure ([§ 4.2.1.5 ]). En effet ce domaine est étendu à tous les points utilisés pour les calculs mécaniques et il permet de préserver la surface de la fissure existante, comme déjà montré (figure ). Les méthodes seront appliquées seulement aux nœuds et aux éléments contenus dans ce domaine, qui est défini par le cercle 3 de figure . En 3D ce domaine coïncide avec un tore construit autour du fond de la fissure et de rayon égal à celui du cercle 3. Dans tout les cas le rayon du domaine localisé est donc égal à:
\({R}_{\mathit{loc}}=\Delta a+{R}_{\mathit{calc}G}\)
où \(\Delta a\) est l’avancée de la fissure à imposer (choisie par l’utilisateur) et \({R}_{\mathit{calcG}}\) est le rayon du domaine utilisé par CALC_G.
On peut donc tout d’abord sélectionner les nœuds de la grille ou du maillage qui ont une distance au fond de la fissure inférieure ou égale à \({R}_{\mathit{loc}}\) (figure à gauche). Cette distance peut être calculée de façon très rapide parce que pour chaque nœud on connaît déjà son projeté sur le fond de fissure ([§ 2.2 ]). Puis on prend tous les éléments qui ont au moins un des nœuds sélectionnés précédemment (éléments gris de la figure de droite). Cela permet de sélectionner correctement tous les éléments qui sont coupés par la frontière du domaine et qui sont donc partiellement contenus dans le domaine. De cette façon, toute la surface (ou le volume en 3D) du domaine est correctement recouvert par les éléments sélectionnés. Enfin la liste des nœuds du domaine doit être mise à jour en ajoutant les nœuds des éléments partiellement contenus dans le domaine (nœuds marqués avec un petit cercle sur la figure de droite). La façon la plus rapide et facile de le faire est de reconstruire la liste en prenant tous les nœuds des éléments sélectionnés. Une fois que la liste des nœuds à été mise à jour, la valeur du rayon \({R}_{\mathit{loc}}\) du domaine discrétisé, c’est-à-dire du domaine constitué par les nœuds et les éléments sélectionnés, doit être mise à jour en prenant la valeur au nœud du domaine discrétisé le plus loin du fond. Cette opération est nécessaire pour la partie qui va suivre ([§ 5.3 ]).
Figure 5.1-1 : exemple de sélection des nœuds et des éléments du domaine de localisation. On suppose que ce domaine coïncide avec le cercle de la figure. La première étape (à gauche) est la sélection des nœuds à une distance inférieure ou égale à \({R}_{\mathit{loc}}\) . Les nœuds sélectionnés sont marqués avec un petit triangle. Ensuite les éléments qui contiennent au moins un des ces nœuds dans leur définition sont sélectionnés (éléments gris à droite). Enfin on sélectionne tous les nœuds des éléments qui ont été pris: les nœuds marqués avec des petits cercles viennent s’ajouter à ceux qui ont été sélectionnés au début.
Déplacement du domaine localisé en suite à la propagation de la fissure#
La mise à jour et la réinitialisation des level-set sont faites exclusivement aux nœuds définissant le domaine localisé (nœuds marqués avec des triangles et des cercles sur la figure de droite). Cela signifie que les valeurs des level-set aux nœuds à l’extérieur du domaine localisé ne sont pas changées, ce qui pose des problèmes au cours des pas de propagation successifs. La situation est illustrée figure pour un point \({P}_{1}\) dans le plan \(t-n\) (plan \({e}_{1}-{e}_{2}\) de la figure ). Au début de la propagation le fond de la fissure coïncide avec le point \({P}_{1}\) et le domaine localisé avec le cercle \({\Omega}_{1}\) . Pour simuler la propagation, les level-set sont mises à jour et réinitialisées seulement dans ce cercle. La nouvelle position du fond après le premier pas de propagation coïncide avec \({P}_{2}\) . Au deuxième pas de propagation, le fond de la fissure coïncide avec le point \({P}_{2}\) et un nouveau domaine \({\Omega}_{2}\) est construit autour de ce point pour la mise à jour et la réinitalisation des level-sets. La valeur des level-sets n’a pas été mise à jour au pas de propagation précèdent pour tous les points de \({\Omega}_{2}\) : seulement les level-sets des points de \({\Omega}_{2}\) en commun avec \({\Omega}_{1}\) sont à jour (points de la surface grise de la figure ). Pour tous les autres points, la valeur des level-sets est celle que l’on avait au début de la propagation. Les bonnes valeurs en ces autres points doivent donc être calculées au début du deuxième pas de propagation.
Figure 5.2-1 : déplacement du domaine localisé entre deux pas de propagation successifs. Le fond de la fissure se déplace de \({P}_{1}\) à \({P}_{2}\) .
Le calcul des bonnes valeurs des level-sets pour tous les points dans la région blanche du domaine \({\Omega}_{2}\) de la figure peut être fait en utilisant la définition de distance signée des level-sets (bib 14 ). Pour chaque point \(M\) dans cette région (figure ), les level-set peuvent être calculées comme ceci:
\({\varphi}_{n}(\text{M})=\overrightarrow{\text{M}{\text{P}}_{2}}\cdot {n}_{2}\)
\({\varphi}_{t}(\text{M})=\overrightarrow{\text{M}{\text{P}}_{2}}\cdot {t}_{2}\)
Ces équations permettent de calculer la bonne valeur des level-sets avec le bon signe. Il faut remarquer que, dans le cas 3D, le point \({P}_{2}\) est le point projeté du point \(M\) sur le fond de la fissure. Pour chaque point \(M\) , son projeté et la base \({n}_{2}-{t}_{2}\) ont déjà été calculés pendant la phase d’extension de la vitesse du fond de la fissure aux points du domaine de calcul ([§ 2.2 ]). Tous les ingrédients du calcul sont donc déjà disponibles et pour chaque point \(M\) l’évaluation de la bonne level-set est limitée au calcul de deux produits scalaires, ce qui est très rapide à faire.
Figure 5.2-2 : calcul des bonnes valeurs des level-sets aux points de \({\Omega}_{2}\) qui sont à l’extérieur du domaine \({\Omega}_{1}\) du pas de propagation précédent. Le calcul est fait tout au début de la propagation courante, c’est-à-dire avant la mise à jour et la réinitialisation des level-sets.
Il faut remarquer que même si les valeurs des level-sets calculées par les équations ci-dessus ne sont pas très précises, la phase de mise à jour ne sert qu’à l’initialisation de ces dernières et les bonnes level-sets seront récupérées par les phases de réinitialisations.
Limitations sur la valeur du rayon du domaine localisé entre deux propagation successives#
Les équations données au paragraphe 33 pour le calcul des bonnes valeurs des level-sets peuvent être appliquées seulement aux points pour lesquels la base \({n}_{2}-{t}_{2}\) est valide pour le calcul de la distance signée. Cela signifie que le point \(M\) doit appartenir au demi-plan montré figure , qui contient la base \({n}_{2}-{t}_{2}\) , le fond de la fissure \({P}_{2}\) et qui a pour bord une droite normale à l’arête reliant \({P}_{1}\) et \({P}_{2}\) parallèle à \({n}_{2}\) .
Figure 5.3-1 : définition du demi-plan auquel le point \(M\) de la figure doit appartenir pour que les équations de calcul des level-sets du paragraphe 33 soient valides.
Cette condition sur la position du point \(M\) se traduit par une condition sur la valeur maximale du rayon \({R}_{2}\) du domaine \({\Omega}_{2}\) : la valeur maximale est celle pour laquelle les deux points d’intersection (\(A\) et \(B\) ) entre les deux domaines \({\Omega}_{1}\) et \({\Omega}_{2}\) sont sur le bord du demi-plan de validité des équations, comme indiqué figure . Ces deux points sont toujours symétriques par rapport à l’arête reliant \({P}_{1}\) et \({P}_{2}\) . La valeur maximale du rayon \({R}_{2}\) peut donc être calculée facilement:
\({R}_{2}\le \sqrt{\Delta {a}^{2}+{R}_{1}^{2}}\)
La valeur du rayon du domaine actuel est donc bornée par les valeurs de l’avancée de la fissure et du rayon du domaine utilisés au pas de propagation précédent. Normalement il n’y a aucun intérêt à changer la valeur du rayon du domaine entre deux propagations successives. Dans ce cas, comme \({R}_{1}\) est égal à \({R}_{2}\) , l’inégalité sur \({R}_{2}\) est nécessairement vérifiée quelle que soit la valeur de l’avancée choisie.
Il faut toutefois remarquer que la valeur du rayon du domaine peut changer par rapport à la valeur théorique donnée par l’utilisateur du fait de la discrétisation du domaine ([§ 5.1 ] et figure de droite). Donc, même si on donne \({R}_{2}={R}_{1}\) au pas de propagation courant, la valeur effective du rayon du domaine sera la suivante:
\({R}_{2\text{eff}}={R}_{2}+\Delta {R}_{2}={R}_{1}+\Delta {R}_{2}\)
où \(\Delta {R}_{2}\) est l’incrément du rayon \({R}_{2}\) du fait de la discrétisation du domaine, qui est de l’ordre de la longueur moyenne des arêtes des éléments autour de la frontière du domaine localisé (voir figure ).
Pour la même raison, la valeur effective du rayon \({R}_{1}\) au pas de propagation précédent peut être différente de la valeur théorique donnée par l’utilisateur:
\({R}_{1\text{eff}}={R}_{1}+\Delta {R}_{1}\)
où \(\Delta {R}_{1}\) est l’incrément du rayon \({R}_{1}\) du fait de la discrétisation du domaine, qui est de l’ordre de la longueur moyenne des arêtes des éléments autour de la frontière du domaine localisé au pas de propagation précédent.
Ainsi l’inégalité ci-dessus peut être écrite dans la forme suivante:
\({R}_{1}+\Delta {R}_{2}\le \sqrt{\Delta {a}^{2}+{({R}_{1}+\Delta {R}_{1})}^{2}}\)
Dans le cas où la longueur des arêtes de la grille auxiliaire utilisée est uniforme dans le domaine, \(\Delta {R}_{2}\approx \Delta {R}_{1}\) , et l’inégalité est toujours vérifiée. Par contre, si cette longueur change dans le domaine, on pourrait avoir \(\Delta {R}_{2}>\Delta {R}_{1}\) et l’inégalité pourrait ne pas être vérifiée.
Dans ce dernier cas on doit forcément raffiner le maillage (ou la grille auxiliaire) pour réduire la différence entre \(\Delta {R}_{1}\) et \(\Delta {R}_{2}\) ou bien utiliser une avancée de fissure \(\Delta a\) plus élevée.
Calcul direct des nouvelles level-sets après avancée#
Une méthode alternative aux méthodes de mise à jour des level-sets est celle décrite dans (bib 18 ) . Elle est dérivée de l’analyse des équations de mise à jour des level-sets décrites aux paragraphes 3.1 : on arrive à montrer que le processus de mise à jour des level-sets peut toujours être ramené à un processus 2D en mode \(I+\mathit{II}\) dans le plan formé par la tangente et la normale à la surface de la fissure en fond de fissure (plan \({n}_{P}-{t}_{P}\) de la figure ). Cela permet de simplifier énormément le calcul direct des nouvelles valeurs des level-sets après avancée. Ce calcul a été effectivement déjà explicité au paragraphe § 5.2 pour ce qui concerne l’évaluation des level-sets pour les points hors domaine de calcul qui sont rajoutés. Les équations qui y sont décrites peuvent être utilisées dans le cas où la nouvelle position du fond est connue. En sachant que le problème de mise à jour des level-sets peut toujours être ramené à un problème 2D, la nouvelle position du fond peut être facilement calculée à partir de l’avancée de la fissure et de l’angle de propagation au point du fond considéré.
Méthode de calcul des level-sets#
Ces idées sont reprises figure , qui montre la configuration des iso-surfaces des deux level-sets autour d’un point \(P\) du fond d’une fissure 3D obtenue par le processus de mise à jour décrit au paragraphe § 3.1
Figure 6.1-1 : iso-surfaces des deux level-sets obtenues par le processus standard de mise à jour, avec réinitialisation et réorthogonalisation. \(P\) est un point du fond d’une fissure 3D. \(Q\) est la nouvelle position de ce point après une avancée de \(\Delta a\) avec un angle de propagation \(\beta\) .
On considère un point \(M\) du domaine de calcul où on doit calculer la nouvelle valeur des deux level-sets. On projette ce point sur le fond de la fissure en utilisant le même algorithme que celui décrit au paragraphe 2.2 et on obtient le point \(P\) (figure ). La base locale \(({t}_{P},{n}_{P})\) et la vitesse d’avancée \({V}_{P}\) au point \(P\) peuvent être calculées comme décrit aux paragraphes 2.2 , 2.3 et 2.4 . Pour l’angle de propagation, on utilise la même interpolation linéaire que celle utilisée pour le vitesse (voir paragraphe § 2.2 et figure :
L’avancée de la fissureau point \(P\) du fondpeut donc être aisément calculée:
Où \(\Delta t\) est le temps ou, en termes équivalents, le nombre de cycles de fatigue à simuler au pas de propagation courant (voir discussion du paragraphe 3.1 ).
La nouvelle position \(Q\) du point \(P\) du fond peut donc être calculée comme ceci:
Où \(O\) est le point origine de la base globale dans laquelle les coordonnées des points sont exprimées.
Enfin la base au point \(Q\) peut être calculée par une rotation d’angle de la base au point P autour de l’axe normal au plan \({t}_{P}-{n}_{P}\) :
Tous les ingrédients du calcul des level-sets sont maintenant connus. On peut donc utiliser les mêmes équations que celles utilisées au paragraphe 5.2 :
Afin de préserver la surface existante de la fissure (iso-surfaces en gris de la figure ), le calcul de la level-set normale avec l’équation ci-dessus est fait uniquement aux points du domaine où \({\varphi}_{t}(\text{M})>0\) .
Traitement des points problématiques#
Les équations écrites au paragraphe précédent s’appliquent si le vecteur \(\overrightarrow{\text{PM}}\) est dans le plan \({t}_{P}-{n}_{P}\) (figure ) (bib 18 ). Cela n’est pas le cas pour tous les points qui n’ont pas un point de projection normale sur le fond de la fissure. Dans ce cas un point extrémité du fond est sélectionné comme point projeté P par l’algorithme de projection décrit au paragraphe 2.2 .
Ces points problématiques peuvent être traités par une rotation de la base \(({t}_{P},{n}_{P})\) qui permet de rétablir la condition de coplanarité du plan \({t}_{P}-{n}_{P}\) et du vecteur \(\overrightarrow{\text{PM}}\) ou \(\overrightarrow{\text{QM}}\) (bib 18 ).
Figure 6.2-1 : a) exemple de point problèmatique. b) rotation de la base au point \(P\) pour rétablir la condition de coplanarité entre le plan \({t}_{P}-{n}_{P}\) et les vecteurs \(\overrightarrow{\text{PM}}\) ou \(\overrightarrow{\text{QM}}\) .
Un exemple de point problématique est visible sur la figure a. Pour rétablir la condition de coplanarité, on peut faire une rotation \(\theta\) de la base \(({t}_{P},{n}_{P})\) autour du vecteur \({n}_{P}\) de telle manière que le plan \({t}_{P}^{\text{*}}-{n}_{P}\) contienne le point \(M\) (figure b). En effet, de cette façon la nouvelle position du point \(P\) après propagation (point \(Q\) de la figure ) est contenue dans le plan \({t}_{P}^{\text{*}}-{n}_{P}\) , ce qui permet de vérifier la condition de coplanarité.
Le nouveau vecteur \({t}_{\text{P}}^{\text{*}}\) peut être calculé à partir du projeté \({M}^{\text{*}}\) du point \(M\) sur le plan normal à l’axe \({n}_{P}\) (figure b). On calcule d’abord:
Ce qui permet d’obtenir le nouveau vecteur:
Pour prendre en compte correctement le signe de la level-settangente, le vecteur \({t}_{\text{P}}^{\text{*}}\) doit être réorienté suivant la directiondonnéepar \({t}_{\text{P}}\) (bib 18 ):
Si les deux vecteurs \({t}_{\text{P}}^{\text{*}}\) et \({t}_{\text{P}}\) sont orthogonaux, l’équation ci-dessus ne peut pas être appliquée mais on peut écrire tout de suite que:
L’utilisation du vecteur \({t}_{\text{P}}^{'}\) à la place de \({t}_{\text{P}}\) dans les équations du paragraphe § 6.1 permet de bien calculer les level-sets pour les points problématiques.
Discontinuité de la level-set normale et limitations de calcul#
Le fait de calculer la level-set normale seulement aux points \(M\) pour lesquels \({\varphi}_{t}(\text{M})>0\) entraîne une discontinuité de la level-set normale, comme démontré dans (bib 18 ). La situation est illustrée figure . La level-set normale n’a pas été recalculée dans la zone colorée en gris impliquant une discontinuité des iso-surface entre la zone en gris et la zone où la level-set a été récalculée.
Figure 6.3-1 : discontinuité des iso-surface de la level-set normale et rayon maximal \({R}_{max}\) du domaine de calcul de CALC_G à utiliser. Comme démontré dans (bib 18 ), cette discontinuité ne cause aucun problème pour les calculs X-FEM (identification des nouvelles mailles coupées notamment). De plus, la continuité de l’iso-zéro de la level-set normale est préservée, ce qui assure la continuité de la surface de la fissure. La seule condition à respecter pour gérer cette discontinuité porte sur la valeur maximale du rayon de calcul de CALC_G(voir U4.82.03). Cette condition a pour but d’éviter que le domaine de calcul de CALC_Gcontienne la discontinuité (cercle gris de la figure ):
Où on prend la valeur maximale de \(\Delta a\) et \(\beta\) le long du fond de la fissure.
Si elle restreint l’angle de propagation à des valeurs strictement inférieures à \(\pi /2\) (sous peine d’un raffinement de plus en plus important à mesure que l’on se rapproche de \(\pi /2\) ), cette condition n’est pas très restrictive parce qu’il est conseillé de respecter une relation similaire même dans le cas de mise à jour des level-sets par la méthode simplexe ou upwind.
Localisation du domaine#
Comme montré dans (bib 18 ), la localisation du domaine décrite au paragraphe § 5 n’est plus nécessaire si on utilise le calcul direct des level-sets.
Tout d’abord, la localisation du domaine n’améliorerait pas la robustesse de l’algorithme. Dans (bib 14 ) , il est montré que de faibles erreurs dans le calcul de l’angle de propagation (ce qui est toujours le cas lors de simulations numériques) sont amplifiées aux points situés loin du front de fissure (du fait de l’interpolation linéaire utilisée pour créer le champ de vitesse \({V}_{N}\) ). Les iso-surfaces de la level-set normale ne sont alors plus régulières et ceci peut entraîner l’échec des phases de réinitialisation des level-sets . La localisation du domaine permet ainsi de réduire le domaine de mis à jour des level-sets à un domaine englobant le front de fissure ce qui conduit à une limitation de l’effet d’amplification. Cependant la méthode géométrique gère avec succès les iso-surfaces non régulières puisqu’elle est basée sur des équations polynomiales explicites et non sur des équations différentielles. De plus, la méthode géométrique n’utilise pas de processus de réinitialisation et de réorthogonalisation des level-sets (qui sont calculés à chaque pas de propagation). Pour toutes ces raisons, la localisation du domaine n’affecte pas la robustesse de la méthode géométrique.
La localisation du domaine dans le cadre de l’utilisation de la méthode géométrique ne réduirait pas non plus forcément le temps de calcul nécessaire à la mise à jour des level-sets. Il s’avère que le calcul des level-sets en chaque point du domaine est rapide et que le temps nécessaire pour ce calcul peut être plus faible que le temps nécessaire pour sélectionner les points du domaine localisé. L’algorithme utilisé pour la sélection des points appartenant au domaine restreint est basé sur le calcul de la distance entre un point \(M\) du domaine et son projeté \(P\) sur le front de fissure. Comme le point \(P\) doit être calculé également pour le calcul des level-sets, le coût de calcul supplémentaire pour l’algorithme de sélection se réduit au coût du calcul de la distance \(∣∣\vec{\mathit{PM}}∣∣\) et la complexité de l’algorithme est de \(O(n)\) . Même si la complexité de l’algorithme de mise à jour des level-sets par la méthode géométrique est de \(O(n)\) , il apparaît clair que la somme des calculs effectués par l’algorithme de sélection sera moindre que celle effectuée par l’algorithme de mise à jour des level-sets. Ceci conduit à penser à première vue que la localisation du domaine serait avantageuse en terme de temps de calcul par rapport au calcul des level-sets dans tout le domaine. Cependant, la sélection (basée sur la distance \(∣∣\vec{\mathit{PM}}∣∣\) ) des points à l’intérieur du domaine restreint n’est pas suffisante. En effet, pour le calcul de la base locale en chaque point, il s’agit de déterminer le gradient des level-sets ce qui nécessite la sélection des éléments contenant les points sélectionnés (les fonctions de forme des éléments étant utilisées pour calculer les gradients). Même si l’algorithme de sélection ne change pas nécessairement, le nombre d’opérations effectuées augmente et ceci peut rendre le temps de calcul de la sélection plus élevé que le calcul des level-sets par la méthode géométrique dans tout le domaine. A noter que, dans le même temps, le temps de calcul des gradients des level-sets est aussi plus élevé si la localisation du domaine n’est pas prise en compte.
Finalement, l’avantage de la localisation du domaine en terme de temps de calcul dépend fortement des algorithmes et implémentations choisis et du nombre de points du domaine de calcul. La localisation du domaine pourrait ainsi être avantageuse dans le cas de maillages éléments finis «très gros». Dans ce cas, la limitation de la mise à jour des level-sets à un groupe de nœuds (fixé et défini par l’utilisateur) définissant la surface potentielle de propagation de la fissure peut être un moyen plus efficace de réduire le temps de calcul plutôt que le choix d’un algorithme de localisation du domaine pour lequel un nouveau domaine restreint doit être calculé à chaque pas de propagation afin de suivre le déplacement du front de fissure.
Méthode maillage pour la mise à jour des level-sets#
La méthode la plus intuitive pour mettre à jour les valeurs des level-sets à chaque pas de propagation de la fissure consiste à utiliser leur propriété de distance signée. Ceci implique une représentation explicite de la fissure (lèvres de la fissure et fond de fissure) propagée afin de calculer les distances des nœuds à la fissure. Cette méthode utilise ainsi deux maillages distincts (fissure et structure) et seul le maillage de la fissure se trouve modifié au cours de la propagation. La mise à jour du maillage de la fissure à chaque pas de propagation consiste concrètement en la création d’une nouvelle rangée de mailles surfaciques pour représenter l’avancée de la fissure et d’une nouvelle rangée de mailles linéiques pour représenter le nouveau fond. Une fois ce maillage mis à jour, on recalcule les nouvelles level-sets par calcul direct de la distance à la fissure (projection orthogonale sur le maillage de la fissure).
Présentation de la méthode#
Le maillage de la fissure est réalisé dans Code_Aster par la commande PROPA_FISS. La donnée d’entrée est la valeur du vecteur propagation sur le front de fissure, déterminée à partir des facteurs d’intensité des contraintes et de la loi de Paris. Le principe de la méthode maillage est de calculer les nœuds du nouveau fond à partir des vitesses de propagation et des nœuds du fond actuel. Il se résume aux trois étapes suivantes:
Calcul par interpolation linéaire du vecteur propagation à chaque nœud du front de fissure discrétisé. Le vecteur propagation est orthogonal au front de fissure: la direction de propagation en un nœud correspond à la moyenne des vecteurs normaux des deux éléments linéaires se partageant ce nœud. Cette étape donne la position de chaque point propagé.
Ajout d’un nœud en chacun de ces points, création d’une nouvelle surface (ajout d’éléments linéaires pour le nouveau front de fissure et d’éléments quadrangles pour la nouvelle surface) et processus de maillage de cette surface créée.
Traitement spécifique des nœuds du bord qui, selon les cas, peuvent «sortir» de la structure. Ils sont alors déplacés sur les frontières de la structure.
Pour les modélisations 2D, la méthode est même plus simple. Seul un élément linéaire doit être ajouté à chaque pas de propagation et aucune correction des points du bord n’est nécessaire.
Algorithme#
L’algorithme mis en place pour la propagation par la méthode maillage est le suivant (comme décrit dans (bib 20 ) et sur la figure ). A l’itération \(i\) :
Récupération du maillage de la fissure,
Récupération des points extrémités du front de fissure,
Correction du maillage de la fissure:
On modifie les coordonnées des deux nœuds extrémités du front de fissure pour qu’ils correspondent aux points extrémité du front de fissure
On modifie les coordonnées de tous les autres nœuds du fond pour qu’ils soient équi-répartis entre les deux extrémités (nombre de nœuds du fond \(\mathit{NbNf}\) ).
Boucle sur les nœuds du fond: \(1⩽k⩽\mathit{NbNf}\)
Calcul de \(\Delta G(k)\) par interpolation linéaire,
Calcul de l’angle de bifurcation \(\beta\) à partir de \({K}_{I}(k)\) et \({K}_{\mathit{II}}(k)\) ,
Calcul du vecteur de propagation par interpolation linéaire:
\(\Delta a(k)=C{\left(\Delta G(k)\right)}^{m/2}\) si \(C\) et \(m\) sont renseignés,
\(\Delta a(k)={\mathit{DA}}_{\mathit{MAX}}{\left(\Delta G(k)/{max}_{k}(\Delta G(k))\right)}^{m/2}\) si \({\mathit{DA}}_{\mathit{MAX}}\) et \(m\) sont renseignés,
Calcul de la position du propagé du nœud \({N}_{k}\) ,
Ajout d’un nouveau nœud.
fin de boucle
Boucle entre \(1\) et \(\mathit{NbNf}-1\)
ajout d’une nouvelle maille SEG2 entre chacun des nœuds du fond et le nœud suivant,
ajout d’une maille QUA4 entre deux nœuds successifs du fond à l’itération \(i-1\) et leurs propagés,
en 2D: ajout d’une maille POI pour le fond et SEG2 entre le fond précédent et le nouveau fond.
Fin de boucle,
Création des mailles \(\mathit{FOND}\) (nouveau fond) et \(\mathit{LEVRE}\) (lèvre précédente + nouvelles mailles surfaciques),
Réactualisation du maillage.
Figure 7.2-1 : a. maillage en entrée de la routine, b. étape de correction des points du bord, c. étape d’équirépartition de tous les autres nœudsdu front, d. propagation du front discrétisé et création des nouvelles mailles ( SEG2 et QUA4 ).
Séparation et fusion des fronts de fissure#
Différentes méthodes ont été présentées pour suivre la propagation d’une fissure avec X-FEM (ie pour mettre à jour les valeurs des level-sets à chaque pas de propagation):
La méthode maillage(ou par projection), voir paragraphe § 7 ,
La méthode géométrique, voir paragraphe§ 6 .
Les méthodes upwind et simplexe basées sur des algorithmes de fast marching
Il s’agit d’expliquer ici les modifications à apporter à ces différentes méthodes pour rendre compte du comportement d’un front de fissure lorsqu’il rencontre un obstacle (un trou par exemple). Comment gérer la séparation de ce front de fissure en deux, voire en plusieurs autres morceaux? De même qu’en est-il de la fusion de deux fronts de fissure indépendants en un seul? Ce sont les level-sets qui déterminent le fond de fissure et donc la présence ou non de plusieurs fonds de fissure. Mais il s’agit dans certains cas de procéder à des améliorations des méthodes précédentes pour qu’elles s’adaptent à la séparation du front de fissure (modification de la création du maillage de la fissure pour la méthode maillage, continuité des level-sets dans le trou pour la méthode upwind avec utilisation d’une grille auxiliaire).
Selon les méthodes utilisées pour la gestion de la mise à jour des level-sets, la stratégie mise en place pour rendre compte de la séparation ou de la fusion de fronts de fissure diffère.
Méthode maillage#
Pour la méthode maillage, une première option est d’étendre ce qui est fait au paragraphe 7.2 au cadre de plusieurs fonds de fissure. Les level-sets donnent des informations sur le nombre et la position des extrémités des fonds de fissure (intersection entre le maillage de la structure et l’iso-zéro des deux level-sets). Tout en conservant les nœuds qui constituent les extrémités déjà existantes à l’instant \(i-1\) , les nœuds constituant les nouvelles extrémités sont sélectionnés (figure a.).
Figure 8.1-1 : a. nœuds déterminant les extrémités des fonds de fissure (en bleu : nœuds sélectionnés pour les extrémités déjà existantes, en rouge : nœuds sélectionnés pour les nouvelles extrémités), b. repositionnement des nœuds sélectionnés sur les différents fonds, c. calcul des nœuds propagés et formation des nouvelles mailles.
Comme au paragraphe 7.2 une procédure de repositionnement des nœuds des différents fonds (figure b.) est effectuée. Il ne reste plus qu’à calculer les nœuds propagés et former les mailles (figure c.).
Bien que cette méthode fonctionne correctement, elle présente un désavantage: la perte de nœuds à la séparation de deux fonds. Les mailles peuvent, de plus, se déformer de façon importante lors de la phase de repositionnement des nœuds extérieurs à la structure. Ceci peut notamment conduire à une mauvaise discrétisation de la fissure en propagation hors-plan.
Pour pallier ces difficultés, l’idée a été d’augmenter le nombre de nœuds sur chaque fond de fissure. Reste la question de la création de ces nœuds et du maillage. En remarquant que les mailles représentant la fissure ne servent qu’à effectuer des projections, il est possible de créer une rangée de mailles pour l’avancée de fissures non conforme à la rangée de mailles précédente (figure ). De cette manière, on peut augmenter le nombre de nœuds entre deux pas de propagation sans modifier la discrétisation induite par les vitesses de propagation.
L’algorithme se résume donc comme suit : recherche des nœuds les plus proches des extrémités du fond, repositionnement sur le bord, création de nouveaux nœuds sur le fond actuel positionnés de façon équidistante les uns des autres, calcul des nœuds du nouveau fond. Comme il y a le même nombre de nœuds entre les deux fonds, des quadrangles respectant les surfaces du fond sont créés.
Figure 8.1-2 : maillage non conforme (en rouge; nœuds créés pour la nouvelle rangée de mailles). Le nombre total \(N\) de nœuds sur tous les fonds de fissure doit rester constant entre deux pas de propagation. Pour cela, on répartit les \({N}_{i}\) nœuds par fond proportionnellement au rapport entre la longueur curviligne \({s}_{i}\) du fond \(i\) et la longueur curviligne totale des \(n\) fonds:
Où \(E(x)\) représente la partie entière de \(x\) . Comme il peut y avoir des cas pour lesquels la somme des \({N}_{i}\) est inférieure à \(N\) , on distribue les nœuds restants aux fonds ayant la plus grande partie décimale.
Remarque:
La méthode maillage peut être problématique quand la fissure est concave. Ce cas peut être notamment rencontré lors de la réunion de deux fonds de fissure (figure a.) et est absolument à éviter car des problèmes importants sont créés par le calcul des level-sets (valeurs fausses, apparition de faux fonds de fissure).
Figure 8.1-3 : cas problématiques: a.front de fissure concave lors de la réunion de deux fonds de fissures, b. création de mailles pathologiques.
La création de mailles problématiques (concaves par exemple) peut s’avérer être une source de difficultés (voir figure b.). Elle peut apparaître lorsque la vitesse de propagation pour les nœuds du bord est imposée tangente à la surface. Cela peut également créer des erreurs dans le calcul de la level-set tangente et impliquer la création de faux fonds de fissure. Pour éviter la concavité d’une maille dans le cas d’un avancement maximal \({\mathit{Da}}_{max}\) de la fissure trop important, on diminue la valeur de ce dernier pour le calcul suivant. Soit le quadrangle représenté dans le plan \((\mathrm{ABC})\) figure (si le quadrangle n’est pas plan alors le point \(D\) représente la projection du nœud propagé à partir de \(C\) sur le plan \((\mathrm{ABC})\) ). On calcule les coordonnées du point \(M\) d’intersection entre \((\mathit{AB})\) et \((\mathit{CD})\) et on compare les distances entre \(\mathit{AM}\) et \(\mathit{AB}\) , \(\mathit{CM}\) et \(\mathit{CD}\) . On pose \(\left(\vec{{E}_{1}}=\frac{\vec{\mathit{AB}}}{∣∣\vec{\mathit{AB}}∣∣},\vec{{E}_{2}}=\frac{\vec{\mathit{AC}}-(\vec{\mathit{AC}}.\vec{{E}_{1}})\vec{{E}_{1}}}{∣∣\vec{\mathit{AC}}-(\vec{\mathit{AC}}.\vec{{E}_{1}})\vec{{E}_{1}}∣∣}\right)\) une base orthonormée, et \(\vec{\mathit{AD}}({d}_{1,}{d}_{2})\) , \(\vec{\mathit{AC}}({c}_{1,}{c}_{2})\) et \(\overrightarrow{\mathit{AM}}(x,0)\) dans le repère \((A,\vec{{E}_{1}},\vec{{E}_{1}})\) . Par définition de \(M\) , \(det(\vec{\mathit{CM}},\vec{\mathit{CD}})=∣\begin{array}{cc}x-{c}_{1}& {d}_{1}-{c}_{1}\\ -{c}_{2}& {d}_{2}-{c}_{2}\end{array}∣=0\) .
On en déduit \(x=\frac{{d}_{2}{c}_{1}-{d}_{1}{c}_{2}}{{d}_{2}-{c}_{2}}=\frac{(\overrightarrow{\mathit{AD}}.\overrightarrow{{E}_{2}})(\overrightarrow{\mathit{AC}}.\overrightarrow{{E}_{1}})-(\overrightarrow{\mathit{AD}}.\overrightarrow{{E}_{1}})(\overrightarrow{\mathit{AC}}.\overrightarrow{{E}_{2}})}{\overrightarrow{\mathit{CD}}.\overrightarrow{{E}_{2}}}\) .
On introduit le coefficient \(\mathit{Crit}=min\left(\frac{\mathit{AM}}{\mathit{AB}},\frac{\mathit{CM}}{\mathit{CD}}\right)\) et l’avancement maximal rectifié \({\mathit{Da}}_{max}'\) vérifie:
\(\mathit{si}\phantom{\rule{4em}{0ex}}\mathit{Crit}\phantom{\rule{0.5em}{0ex}}<\phantom{\rule{0.5em}{0ex}}1.2\phantom{\rule{10em}{0ex}}\mathit{alors}\phantom{\rule{4em}{0ex}}{\mathit{Da}}_{max}'\phantom{\rule{0.5em}{0ex}}=\phantom{\rule{0.5em}{0ex}}\frac{{\mathit{Da}}_{max}.\mathit{Crit}}{1.2}\)
Figure 8.1-4 : représentation d’une maille pour laquelle on vérifie que l’avancement n’est pas trop grand. Le coefficient \(1,2\) est un coefficient de sécurité afin d’assurer une propagation non nulle (\(\mathit{Crit}=1\) conduisant à un point \(M\) confondu avec le(s) point(s) \(B\) et/ou \(D\) ).
Bien qu’avec une telle stratégie, l’avancement maximal risque de tendre vers 0, elle permet cependant de résoudre une partie des cas problématiques.
Méthode upwind : introduction d’un fond virtuel#
Comme remarqué dans le paragraphe 4.2.1.4 , l’utilisation de la méthode upwind s’accompagne très généralement de l’utilisation d’une grille auxiliaire (dont le maillage est régulier) sur laquelle s’effectue la mise à jour des level-sets qui sont par la suite projetées sur le maillage physique de la pièce (généralement beaucoup moins régulier). Cependant une partie de la grille peut être située à l’extérieur du maillage de la pièce, comme par exemple dans le cas d’une structure contenant un trou.
Les nœuds de la grille situés à l’extérieur du maillage de la pièce sont alors projetés sur le même point du front de fissure lors du calcul de la vitesse de propagation des level-sets. Tous ces nœuds se voient donc attribuer les mêmes vitesses et angles de propagation ainsi que les mêmes bases locales (comme décrit au paragraphe 2.3 ). Il en résulte l’apparition d’une discontinuité des level-sets sur la grille au niveau du bord libre de la pièce (voir figure ).
Figure 8.2-1 : exemple de discontinuité des level-sets sur la grille au niveau du bord libre. Cette discontinuité, bien que située hors du maillage physique, implique une mauvaise définition du gradient des level-sets sur le bord libre et donc une mauvaise définition de la base locale lors du pas de propagation suivant. Elle peut également influer sur la définition des level-sets à l’intérieur du maillage, notamment lors de la propagation de fissures dans des pièces ayant une section irrégulière et/ou des trous. Pour éviter ces problèmes et définir de manière cohérente les level-sets sur la grille au voisinage du maillage physique, une interpolation du front de fissure à l’extérieur du maillage par une droite tangente au front de fissure est proposée. Pour cela, le dernier segment du front de fissure est prolongé par un segment de droite à l’extrémité duquel est introduit un nouveau nœud du front de fissure, qualifié de virtuel. Ce point est alors considéré comme un point du front de fissure. La même base locale que celle du dernier point du front de fissure physique lui est attribuée. La vitesse et l’angle de propagation sont calculés par interpolation linéaire à partir des vitesses et angles de propagation des deux derniers points (figure ).
Figure 8.2-2 : représentation de l’approximation linéaire de la vitesse au point virtuel \({V}_{v}\) (en rouge) à partir des vitesses des deux derniers points du front de fissure \({V}_{i}\) et \({V}_{j}\) (en bleu).
Pour créer les points du fond virtuel, on procède de la manière suivante (figure ):
on pose \({P}_{i}\) et \({P}_{j}\) deux points situés à l’extrémité du font de fissure et \(d\) la distance du point virtuel \({P}_{v}\) au dernier point du fond de fissure \({P}_{j}\) ,
on détermine la droite prolongeant le dernier segment du front de fissure au point virtuel : \(\overrightarrow{{P}_{\mathit{jv}}}=\mathit{d.}\overrightarrow{{b}_{j}}\) où \(\overrightarrow{{P}_{\mathit{jv}}}\) est le vecteur reliant \({P}_{j}\) à \({P}_{v}\) et \(\overrightarrow{{b}_{j}}\) est le vecteur unitaire tangent au front de fissure au point \({P}_{j}\) ,
on calcule les dérivées spatiales de la vitesse et de l’angle de propagation sur le dernier segment du fond virtuel :
on en déduit la vitesse et l’angle de propagation au fond virtuel :
Figure 8.2-3 :création et emplacement du front de fissure virtuel. On place le point virtuel à une distance \(d\) du front physique égale au rayon du tore utilisé pour la propagation des level-sets (voir paragraphe 5 ). On choisit cette distance de façon à obtenir des level-sets lissées sur l’ensemble de la zone de propagation. Si la localisation n’est pas utilisée, le point virtuel est placé à une distance de sept fois l’avancement maximal de la fissure. Pour des vitesses et des angles de propagation (définis sur les points virtuels) trop importants ou trop faibles, on rapproche le point virtuel du front physique de manière à satisfaire les conditions imposées. On note par exemple que l’on a une importante discontinuité des level-sets si l’angle de propagation est supérieur (respectivement inférieur) à 90 degrés (respectivement -90 degrés). On note également que si la distance du point virtuel au front de fissure physique est inférieure au rayon de localisation, on peut conserver une discontinuité des gradients. Néanmoins, cette discontinuité est située loin du maillage de la pièce. Elle n’influe donc pas sur la définition de la base locale sur le bord libre.
Figure 8.2-4 : création d’un front virtuel (en rouge) à l’intérieur d’un trou circulaire. \({P}_{1}\) est créé par prolongement de \([{P}_{i}{P}_{j}]\) et \({P}_{2}\) par prolongement de \([{P}_{l}{P}_{k}]\) *. Le dernier segment du fond virtuel relie* \({P}_{1}\) et \({P}_{2}\) et connecte les fronts de fissure. La méthode du fond virtuel peut alors être immédiatement étendue au cas de trous. Si la pièce présente un trou, il convient en effet également de bien définir les level-sets à l’intérieur de celui-ci. On procède alors comme décrit précédemment en définissant un front de fissure virtuel à l’intérieur du trou par la création de deux points virtuels \({P}_{1}\) et \({P}_{2}\) , un pour chaque fond de fissure délimité par le trou (figure ). Deux des trois segments du front virtuel sont des extensions des deux fronts de fissures physiques. Le troisième segment relie quant à lui les deux points virtuels afin de connecter les fronts de fissure. Il ne sert pas à interpoler les level-sets mais à assurer leur continuité dans le trou. Les deux points virtuels ayant des vitesses de propagation différentes, deux nœuds voisins situés au centre du trou peuvent être projetés sur deux nœuds du fond de fissure différents et avoir des vitesses différentes.
Les angles et vitesses de propagation ainsi que la base locale sur le segment \([{P}_{1}{P}_{2}]\) sont calculés par interpolation linéaire (comme décrit au paragraphe 2.4 ) à partir des valeurs attribuées aux deux points virtuels.
Dans le cas de la création d’un fond virtuel dans un trou, on doit s’assurer qu’il n’y a pas de croisement des fonds de fissure. La distance d’un point virtuel au fond de fissure ne doit pas dépasser le tiers de la distance entre les deux fonds. Soit:
De plus, la vitesse de propagation au point virtuel \({V}_{v}\) ne doit pas être trop petite ou trop importante. On impose qu’elle soit au moins égale aux deux tiers de la vitesse Vj du dernier fond de fissure et qu’elle ne doive pas dépasser quatre tiers de cette valeur. Ainsi, on impose que:
De même:
Et:
On en déduit l’expression de la distance \(d\) dans ces deux cas :
Néanmoins, la création d’un fond virtuel n’est pas suffisante. Les valeurs des level-sets définies à l’intérieur du trou, issues de la projection sur le segment milieu du fond virtuel, peuvent dans certains cas influencer la définition des level-sets dans le maillage. De plus, l’approximation de la vitesse est plus précise si les level-sets sont initialement droites et si l’intersection des deux iso-zeros est confondue avec le front virtuel. Dans le but de garantir le redressement des level-sets dans le trou lors de la fusion des fronts de fissure après le passage de celui-ci, on recalcule les level-sets pour les nœuds projetés sur un segment virtuel du front de fissure, et ce avant la propagation des level-sets. On calcule les level-sets de manière géométrique (voir paragraphe 6 ) en projetant le nœud sur le front de fissure (figure ). On pose \(\mathit{XM}\) un nœud du maillage. On définit alors son projeté \(\mathit{XN}\) sur le fond de fissure, les vecteurs \(\vec{\mathit{tp}}\) et \(\vec{\mathit{np}}\) définissant la base locale en front de fissure en \(\mathit{XN}\) . On obtient les level-sets normale et tangente aux nœuds \(\mathit{XM}\) simplement par: \(\mathit{lsn}(\mathit{XM})=\overrightarrow{\mathit{XM}\mathit{XN}}.\overrightarrow{\mathit{np}}\) et \(\mathit{lst}(\mathit{XM})=\overrightarrow{\mathit{XM}\mathit{XN}}.\overrightarrow{\mathit{tp}}\) .
Figure 8.2-5 : Calcul des level-sets aux nœuds projetés sur le fond de fissure. Ce calcul des level-sets permet de corriger l’influence de la définition non physique des level-sets au centre du trou lors de la fusion des deux fonds de fissure. Une fois le fond virtuel créé et les level-sets corrigées, le calcul des vitesses de propagation par projection sur le front de fissure s’effectue comme décrit au paragraphe 2.3 , sans faire de distinction entre le front physique et le fond virtuel.
Remarque :
L’introduction d’un fond virtuel nécessite une stratégie de vérification de la numérotation des fronts de fissure ainsi que de leur sens de parcours ce qui n’était jusqu’alors pas nécessaire: la connexion des fronts de fissure par un fond virtuel ne les rend plus indépendants les uns des autres. La méthode suivante a été adoptée pour vérifier si les nœudsdes différents fronts de fissure sont bien numérotés ( figure ):
on construit un nouveau front de fissure sur la grille auxiliaire qui par définition de la grille n’est constitué que d’un seul morceau avec un unique sens de parcours. On projette ensuite les extrémités \({E}_{1}^{i}\) et \({E}_{2}^{i}\) de chacun des morceaux \(i\) du front de fissure réelle sur le front de fissure continu associé à la grille. Ces extrémités sont repérées sur le front de fissure de la grille par leur abscisse curviligne respective. A titre d’exemple, aux extrémités \({E}_{1}^{i}\) et \({E}_{2}^{i}\) correspondent les abscisses curvilignes \({s}_{{P}_{1}^{i}}\) et \({s}_{{P}_{2}^{i}}\) (où \({P}_{1}^{i}\) et \({P}_{2}^{i}\) sont respectivement les projetés de \({E}_{1}^{i}\) et \({E}_{2}^{i}\) sur le front de fissure associé à la grille auxiliaire);
la séquence des abscisses curvilignes ainsi obtenue permet une renumérotation des morceaux des fronts de la fissure réelle;
on évalue le produit scalaire \(({E}_{2}^{i}{E}_{1}^{i}).({P}_{2}^{i}{P}_{1}^{i})\) pour vérifier que chacun des morceaux est parcouru dans le même sens que le front de fissure sur la grille. Dans le cas contraire, on réoriente le(s) morceau(x) mal orienté(s).
Figure 8.2-6 : a. exemple de mauvaise numérotation et orientation du front de fissure réelle (noire) et front de fissure sur la grille (bleu) , b. renumérotation des morceaux du front de fissure réelle, c. réorientation des morceaux du front de fissure réelle
Méthode géométrique et simplexe#
La méthode géométrique et la méthode simplexe ne nécessitant aucune grille auxiliaire, la gestion de la séparation et de la fusion de fronts de fissure est automatique et ne demande aucune modification particulière.
Bibliographie#
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
Sukumar N., Chopp D.L., Moran B.,“Extended finite element method and fast marching method for three-dimensional fatigue crack propagation”, Engineering Fracture Mechanics, vol. 70, pp. 29-48, 2003
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
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
Crandall M.G., Lions P.L., “Two approximations of solutions of Hamilton-Jacobi equations”, Mathematics of Computation, vol. 43, pp. 1-19, 1984
Deconinck H., Struijs R., RoeP.L., “Compact Advection Schemes on Unstructured Grids”, Technical Report, VKI, VKI LS 1993-04, Computational Fluid Dynamics, 1993
RoeP.L., “Linear Advection Schemes on Triangular Meshes”, Technical Report CoA 8720, Cranfield Institute of Technology, 1987
RoeP.L., “Optimum Upwind Advection on Triangular Mesh”, ICASE 90-75, 1990
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
ADALSTEINSSON D., SETHIAN J.A., “A fast level-set method for propagating interfaces”, Journal of Computational Physics, vol. 118, pp. 269-277, 1995
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
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
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
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
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
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
COLOMBO D., «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
SETHIAN J.A., Level Set Methods and Fast Marching Methods, 1999.
GALENNE E., «Propagation automatique de fissures 3D avec X-FEM en fatigue: méthode de projection» Compte rendu AMA, 2008
HABOUSSA D., «Modélisation de la transition traction-cisaillement des métaux sous choc par la X-FEM», Thèse, INSA Lyon, 2012
BRONSTEIN A., BRONSTEIN M., KIMMEL R .«Numerical geometry of non-rigid shapes», Springer, 2007
Description des versions#
Indice document |
Version Aster |
Auteur(s) Organisme(s) |
Description des modifications |
A |
11 |
D. COLOMBO (University of Manchester), M. GUITON, M. SIAVELIS (IFPen), S. GENIAUT, V.X. TRAN EDF/R&D AMA, P.MASSIN EDF/R&D LaMSID |
Texte initial |