1. Problème de référence

1.1. Géométrie

Modélisation A :

La structure 2D est une plaque carrée unitaire (\(\mathrm{LX}=1\) , \(\mathrm{LY}=1\) ), comportant une fissure «débouchante» sur une demi-longueur []. La fissure est droite, horizontale et de longueur \(a=0,5\) . On oriente arbitrairement la fissure du bord droit vers le centre.

Les bords du domaines sont notés en sens trigonométrique :

• \(\mathit{LIG1}\) désigne le bord inférieur.

• \(\mathit{LIG2}\) désigne le bord droit.

• \(\mathit{LIG3}\) désigne le bord supérieur.

• \(\mathit{LIG4}\) désigne le bord gauche.

Les 4 bords du domaine servent à imposer les conditions limites. Sur les bords non coupés par la fissure (\(\mathit{LIG1}\) , \(\mathit{LIG2}\) , \(\mathit{LIG3}\) ) on impose des conditions limites de Dirichlet, à l’aide de la solution analytique en déplacement (voir paragraphe [ 4 ]).

Notons que le bord droit (\(\mathit{LIG4}\) ) coupé par la fissure, engendre une singularité à l’intersection. On observe un saut de déplacement correspondant à l’ouverture de la fissure []. Il est difficile de contrôler ce bord en déplacement, puisqu’il faut expliciter analytiquement la condition de saut sur la lèvre supérieure et sur la lèvre inférieure de la fissure sur les éléments de bord coupés par la fissure. On contourne cette difficulté en imposant des conditions de Neumann sur ce bord.

Figure 1.1-1: Géométrie du domaine

Modélisation B :

Même géométrie que la modélisation A.

Modélisation C :

On incline la fissure d’un angle variable tel que \(\alpha \epsilon \lbrace {0}^{°}{,30}^{°}{,60}^{°}{,90}^{°}{,120}^{°}\rbrace\) [].

Figure 1.1-2: Fissure inclinée

Notons que la fissure s’est allongée. La nouvelle longueur est: \(a=\frac{0,5}{max\lbrace ∣\cos(\frac{\alpha \times \pi }{180})∣,∣\sin(\frac{\alpha \times \pi }{180})∣\rbrace }\) .

1.2. Propriétés du matériau

Module d’Young: \(E={10}^{5}\mathrm{Pa}\)

Coefficient de Poisson: \(\nu =0\)

1.3. Conditions aux limites et chargements

Le chargement est imposé grâce à des conditions limites mixtes.

Les bords non fissurés sont contrôlés en déplacement, le bord fissuré est contrôlé en force. Par ailleurs, les conditions limites de Dirichlet (en déplacement) fixent la structure et empêchent l’apparition de modes rigides.

Dans les modélisations A et B, la fissure coupe uniquement le bord \(\mathit{LIG4}\) , On impose une condition de Neumann sur ce bord et de Dirichlet sur les 3 autres bords (voir [] ).

Dans la modélisation C, on généralise la démarche précédente. Pour une inclinaison supérieure à 45°, la fissure coupe soit le bord inférieur (\(\mathit{LIG1}\) ), soit le bord supérieur (\(\mathit{LIG3}\) ). Sur ces deux bords aussi on applique une condition de Neumann. On impose donc une condition de Dirichlet sur le bord restant (\(\mathit{LIG2}\) ) pour fixer les modes rigides(voir [] ). Les conditions limites étant imposées symétriquement à la fissure horizontale (\(\alpha =0\) ), les calculs présentés seront valables pour \(-{135}^{°}<\alpha <{135}^{°}\)

Figure 1.3-1: Conditions limites mixtes pour les modélisations A et B

Figure 1.3-2: Conditions limites mixtes pour la modélisation C

Le déplacement imposé correspond à la solution analytique exacte:

\({U}_{x}=\frac{(1+\nu )}{E}\sqrt{\frac{r}{2\pi }}{K}_{I}\cos(\frac{\theta}{2})(3-4\nu -\cos\theta )\)

\({U}_{y}=\frac{(1+\nu )}{E}\sqrt{\frac{r}{2\pi }}{K}_{I}\sin(\frac{\theta}{2})(3-4\nu -\cos\theta )\)

Cette solution dépend des coordonnées polaires liées au référentiel de la fissure []. Dans la littérature ([ 13 ]), la direction de la fissure supporte l’axe \(\overrightarrow{X}\) . L’axe \(\overrightarrow{X}\) est orienté du fond de la fissure vers l’extérieur. \(\overrightarrow{Y}\) est orthogonal à \(\overrightarrow{X}\) , tel que \((\overrightarrow{X},\overrightarrow{Y})\) forme un repère direct.

Figure 1.3-3: Référentiel local des formules analytiques

En conséquence, il faut effectuer un changement de base pour adapter les équations analytiques valables dans le référentiel de la fissure, à un repère donné. En particulier, comme l’origine du repère des coordonnées des mailles est choisi au coin inférieur gauche du domaine, on translate les coordonnées polaires:

\(\theta\) est l’angle polaire: \(\theta (x,y)=\mathit{arctan2}(y-0.5,x-0.5)\)

\(r\) la distance radiale: \(r(x,y)=\sqrt{{(x-0.5)}^{2}+{(y-0.5)}^{2}}\)

Le tenseur des contraintes solution analytique est :

\({\sigma}_{xx}=\frac{{K}_{I}}{\sqrt{2\pi r}}\cos\left(\frac{\theta}{2}\right)\left(1-\sin\left(\frac{\theta}{2}\right)\sin\left(\frac{3\theta }{2}\right)\right)\)

\({\sigma}_{yy}=\frac{{K}_{I}}{\sqrt{2\pi r}}\cos\left(\frac{\theta}{2}\right)\left(1+\sin\left(\frac{\theta}{2}\right)\sin\left(\frac{3\theta }{2}\right)\right)\)

\({\sigma}_{xy}=\frac{{K}_{I}}{\sqrt{2\pi r}}\sin\left(\frac{\theta}{2}\right)\cos\left(\frac{\theta}{2}\right)\cos\left(\frac{3\theta }{2}\right)\)

Ensuite, on projette le tenseur des contraintes pour en déduire la densité de force à appliquer sur le bord droit (\(\mathit{LIG4}\) ), de normale «sortante» \(-\overrightarrow{X}\) .

La force à appliquer est: \(\overrightarrow{F}=\underline{\underline{\sigma}}.-\overrightarrow{X}\)

Il vient: \({F}_{x}=-{\sigma}_{xx}\) et \({F}_{y}=-{\sigma}_{xy}\)

Modélisation A :

Par rapport au référentiel de la fissure, les axes «aster» sont orientés dans le même sens. Les champs de déplacement \(\mathit{Ux}\) et \(\mathit{Uy}\) de référence sont appliqués sur les bords LIG1,LIG2, LIG3:

\({U}_{x}=\frac{E}{2(1+\nu )}\sqrt{\frac{r}{2\pi }}{K}_{I}\cos(\frac{\theta}{2})(3-4\nu -\cos\theta )\)

\({U}_{y}=\frac{E}{2(1+\nu )}\sqrt{\frac{r}{2\pi }}{K}_{I}\sin(\frac{\theta}{2})(3-4\nu -\cos\theta )\)

De même le tenseur de contraintes référence est conservé. On applique donc la même densité de force explicitée ci-dessus sur le bord \(\mathit{LIG4}\) .

\({F}_{x}=-{\sigma}_{xx}\) \({F}_{y}=-{\sigma}_{xy}\)

Modélisation B :

Mêmes équations que la modélisation A.

Modélisation C :

On incline la fissure d’un angle de :math:`alpha epsilon lbrace {0}^{°}{,30}^{°}{,60}^{°}{,90}^{°}{,120}^{°}rbrace ` .

Cette rotation impacte l” angle polaire , le champ de déplacement et le tenseur des contraintes . Ces grandeurs effectuent une rotation d’un angle \(\alpha\) , par rapport au repère fixe du code aster.

\(\theta \to \theta -\alpha\)

Cependant cette transformation de l’angle polaire n’est pas valable pour tous les angles, parce que la définition de l’angle polaire n’est pas continue. Dans le référentiel local, de part et d’autre de la fissure on observe un saut angulaire de \(2\pi\) : dela lèvre inférieure à la lèvre supérieure on passe de \(-\pi\) à \(\pi\) .

Par conséquent, la rotation de la fissure doit aussi propager cette discontinuité angulaire. Tous les angles balayés par la fissure (zone grisée cf. [] et []) lors de la rotation, subissent un saut angulaire de \(-2\pi\) ou \(2\pi\) ,

Si \(\alpha >0\) (la rotation de la fissure s’effectue vers le bas)

• Dans la zone \(-\pi <\theta <-\pi +\alpha\) on a : \(\theta \to \theta -\alpha +2\pi\)

• Dans la zone \(\pi >\theta >-\pi +\alpha\) on a : \(\theta \to \theta -\alpha\)

Si \(\alpha <0\) (la rotation de la fissure s’effectue vers le haut)

• Dans la zone \(\pi >\theta >\pi +\alpha\) on a : \(\theta \to \theta -\alpha -2\pi\)

• Dans la zone \(-\pi <\theta <\pi +\alpha\) on a : \(\theta \to \theta -\alpha\)

Figure 1.3-4: Domaine de saut angulaire (positif)

Figure 1.3-5: Domaine de saut angulaire (négatif)

\(\overrightarrow{U}\to {R}_{\alpha}\overrightarrow{U}\) où \({R}_{\alpha}\) est la matrice de rotation d’angle \(\alpha\) avec \({R}_{\alpha}=\left[\begin{array}{cc}\cos(\alpha )& -\sin(\alpha )\\ \sin(\alpha )& \cos(\alpha )\end{array}\right]\)

Par conséquent,

\(\overrightarrow{U}\to \left[\begin{array}{c}\cos(\alpha ){U}_{x}-\sin(\alpha ){U}_{y}\\ \sin(\alpha ){U}_{x}+\cos(\alpha ){U}_{y}\end{array}\right]\)

Le tenseur des déformations est aussi impacté par la rotation:

\(\underline{\underline{\sigma}}\to {R}_{\alpha}\underline{\underline{\sigma}}{R}_{-\alpha }\)

\(\underline{\underline{\sigma}}\to \left[\begin{array}{cc}\cos(\alpha )& -\sin(\alpha )\\ \sin(\alpha )& \cos(\alpha )\end{array}\right]\left[\begin{array}{cc}{\sigma}_{xx}& {\sigma}_{xy}\\ {\sigma}_{xy}& {\sigma}_{yy}\end{array}\right]\left[\begin{array}{cc}\cos(\alpha )& \sin(\alpha )\\ -\sin(\alpha )& \cos(\alpha )\end{array}\right]\)

\(\underline{\underline{\sigma}}\to \left[\begin{array}{cc}\cos(\alpha )& -\sin(\alpha )\\ \sin(\alpha )& \cos(\alpha )\end{array}\right]\left[\begin{array}{cc}\cos(\alpha ){\sigma}_{xx}-\sin(\alpha ){\sigma}_{xy}& \sin(\alpha ){\sigma}_{xx}+\cos(\alpha ){\sigma}_{xy}\\ \cos(\alpha ){\sigma}_{xy}-\sin(\alpha ){\sigma}_{yy}& \sin(\alpha ){\sigma}_{xy}+\cos(\alpha ){\sigma}_{yy}\end{array}\right]\)

D’où l’écriture finale du tenseur de contraintes:

\(\underline{\underline{\sigma}}\to \left[\begin{array}{cc}\cos{(\alpha )}^{2}{\sigma}_{xx}-\sin(2\alpha ){\sigma}_{xy}+\sin{(\alpha )}^{2}{\sigma}_{yy}& -\frac{1}{2}\sin(2\alpha )(-{\sigma}_{xx}+{\sigma}_{yy})+\cos(2\alpha ){\sigma}_{xy}\\ -\frac{1}{2}\sin(2\alpha )(-{\sigma}_{xx}+{\sigma}_{yy})+\cos(2\alpha ){\sigma}_{xy}& \sin{(\alpha )}^{2}{\sigma}_{xx}+\sin(2\alpha ){\sigma}_{xy}+\cos{(\alpha )}^{2}{\sigma}_{yy}\end{array}\right]\)

On vérifie que la rotation conserve la symétrie du tenseur et la trace, \(\mathit{tr}(\underline{\underline{\sigma}})={\sigma}_{xx}+{\sigma}_{yy}\)

Dans la programmation du fichier de commande, la fissure horizontale est un cas particulier de fissure inclinée: toutes les équations dépendent d’une inclinaison \(\alpha\) quelconque. Pour retrouver les modélisation A et B, il suffit de faire \(\alpha =0°\) .

1.4. Solution de référence

On impose \({K}_{I}=1\) .

Dans le référentiel lié à la fissure, les équations analytiques deviennent:

\({U}_{x}=\frac{E}{2(1+\nu )}\sqrt{\frac{r}{2\pi }}\cos(\frac{\theta}{2})(3-4\nu -\cos\theta )\)

\({U}_{y}=\frac{E}{2(1+\nu )}\sqrt{\frac{r}{2\pi }}\sin(\frac{\theta}{2})(3-4\nu -\cos\theta )\)

Par linéarité des équations, le tenseur des contraintes est affecté par le même coefficient de proportionnalité et devient:

\({\sigma}_{xx}=\frac{1}{\sqrt{2\pi r}}\cos\left(\frac{\theta}{2}\right)\left(1-\sin\left(\frac{\theta}{2}\right)\sin\left(\frac{3\theta }{2}\right)\right)\)

\({\sigma}_{yy}=\frac{1}{\sqrt{2\pi r}}\cos\left(\frac{\theta}{2}\right)\left(1+\sin\left(\frac{\theta}{2}\right)\sin\left(\frac{3\theta }{2}\right)\right)\)

\({\sigma}_{xy}=\frac{1}{\sqrt{2\pi r}}\sin\left(\frac{\theta}{2}\right)\cos\left(\frac{\theta}{2}\right)\cos\left(\frac{3\theta }{2}\right)\)

Par construction, il n’existe pas de mode \(\mathit{II}\) . Autrement dit, on teste \({K}_{\mathit{II}}=0\) .

Modélisation A :

Cette modélisation permet de valider l’opérateur POST_ERREUR en comparant les valeurs calculées de la norme \({L}^{2}\) du déplacement et de l’énergie de la structure aux valeurs analytiques.

1.4.1. Calcul de la norme L2 du déplacement

Soit \(\Omega\) le domaine occupé par le solide. La norme \({L}^{2}\) du déplacement est définie par:

\({\parallel u\parallel }_{{L}^{2}}^{2}={\int}_{\Omega}{\parallel u\parallel }^{2}\mathit{dS}.\)

On a:

\({\parallel u\parallel }^{2}={\left(\frac{1+\nu }{E}\right)}^{2}\frac{r}{2\pi }{K}_{I}^{2}{(\kappa -\cos\theta )}^{2},\)

pour tout point de \(\Omega\) . On a donc:

\({\parallel u\parallel }_{{L}^{2}}^{2}={\left(\frac{1+\nu }{E}\right)}^{2}\frac{{K}_{I}^{2}}{2\pi }{\int}_{\Omega}r{(\kappa -\cos\theta )}^{2}\mathit{dS}.\)

Le domaine \(\Omega\) est représenté en coordonnées cartésiennes par \(\Omega =[-a,a]\times [-a,a]\) , où \(a=\mathit{LX}/2=\mathit{LY}/2=1/2\) , et peut-être représenté en coordonnées polaires par:

\(\Omega =\left\lbrace (r,\theta )\in ℝ\times \left[\frac{-\pi }{4},\frac{7\pi }{4}\right],r\le \rho (\theta )\right\rbrace ,\)

où \(\rho\) est défini par:

\(\rho (\theta )=\lbrace \begin{array}{c}\frac{a}{\cos(\theta )}\text{pour}\theta \in \left[\frac{-\pi }{4},\frac{\pi}{4}\right],\\ \frac{a}{\cos\left(\theta -\frac{\pi}{2}\right)}\text{pour}\theta \in \left[\frac{\pi}{4},\frac{3\pi }{4}\right],\\ \frac{a}{\cos(\theta -\pi )}\text{pour}\theta \in \left[\frac{3\pi }{4},\frac{5\pi }{4}\right],\\ \frac{a}{\cos\left(\theta -\frac{3\pi }{2}\right)}\text{pour}\theta \in \left[\frac{5\pi }{4},\frac{7\pi }{4}\right].\end{array}\)

Considérons la fonction \(f:(r,\theta )\to {r}^{n}g(\theta )\) . On a:

\(I={\int}_{\Omega}f(r,\theta )\mathit{dS}=\underset{0}{\overset{2\pi }{\int}}\left(\underset{0}{\overset{\rho (\theta )}{\int}}f(r,\theta )r\mathit{dr}\right)d\theta .\)

Donc:

\(I=\underset{0}{\overset{2\pi }{\int}}\left(\underset{0}{\overset{\rho (\theta )}{\int}}{r}^{n+1}g(\theta )\mathit{dr}\right)d\theta =\frac{1}{n+2}\underset{0}{\overset{2\pi }{\int}}g(\theta ){\left(\rho (\theta )\right)}^{n+2}d\theta .\)

On a:

\(J=\underset{0}{\overset{2\pi }{\int}}g(\theta ){\left(\rho (\theta )\right)}^{n+2}d\theta ={J}_{1}+{J}_{2}+{J}_{3}+{J}_{4},\)

avec:

\({J}_{1}={a}^{n+2}\underset{\frac{-\pi }{4}}{\overset{\frac{\pi}{4}}{\int}}\frac{g(\theta )}{\cos{(\theta )}^{n+2}}d\theta ,\)

\({J}_{2}={a}^{n+2}\underset{\frac{\pi}{4}}{\overset{3\frac{\pi}{4}}{\int}}\frac{g(\theta )}{\cos{\left(\theta -\frac{\pi}{2}\right)}^{n+2}}d\theta ={a}^{n+2}\underset{\frac{-\pi }{4}}{\overset{\frac{\pi}{4}}{\int}}\frac{g\left(\theta +\frac{\pi}{2}\right)}{\cos{(\theta )}^{n+2}}d\theta ,\)

\({J}_{3}={a}^{n+2}\underset{\frac{3\pi }{4}}{\overset{\frac{5\pi }{4}}{\int}}\frac{g(\theta )}{\cos{(\theta -\pi )}^{n+2}}d\theta ={a}^{n+2}\underset{\frac{-\pi }{4}}{\overset{\frac{\pi}{4}}{\int}}\frac{g(\theta +\pi )}{\cos{(\theta )}^{n+2}}d\theta ,\)

\({J}_{2}={a}^{n+2}\underset{\frac{5\pi }{4}}{\overset{\frac{7\pi }{4}}{\int}}\frac{g(\theta )}{\cos{\left(\theta -\frac{3\pi }{2}\right)}^{n+2}}d\theta ={a}^{n+2}\underset{\frac{-\pi }{4}}{\overset{\frac{\pi}{4}}{\int}}\frac{g\left(\theta +\frac{3\pi }{2}\right)}{\cos{(\theta )}^{n+2}}d\theta ,\)

On a donc finalement:

\(I=\frac{{a}^{n+2}}{n+2}\underset{\frac{-\pi }{4}}{\overset{\frac{\pi}{4}}{\int}}\frac{1}{\cos{(\theta )}^{n+2}}\left[g(\theta )+g\left(\theta +\frac{\pi}{2}\right)+g(\theta +\pi )+g\left(\theta +\frac{3\pi }{2}\right)\right]d\theta .\)

On remarque que pour \(n=1\) et \(g:\theta \to {(\kappa -\cos(\theta ))}^{2}\) ., on a:

\(I={\int}_{\Omega}r{(\kappa -\cos\theta )}^{2}\mathit{dS}.\)

On remarque aussi que:

\({(\kappa -\cos\theta )}^{2}={\kappa}^{2}+\frac{1}{2}-2\cos(\theta )+\frac{1}{2}\cos(2\theta ).\)

On a donc dans ce cas:

\(g(\theta )+g\left(\theta +\frac{\pi}{2}\right)+g(\theta +\pi )+g\left(\theta +\frac{3\pi }{2}\right)=4{\kappa}^{2}+2.\)

D’où:

\({\int}_{\Omega}r{(\kappa -\cos\theta )}^{2}\mathit{dS}=\frac{{a}^{3}}{3}(4{\kappa}^{2}+2)\underset{\frac{-\pi }{4}}{\overset{\frac{\pi}{4}}{\int}}\frac{1}{\cos{(\theta )}^{3}}d\theta .\)

Et finalement:

\({\parallel u\parallel }_{{L}^{2}}^{2}={\left(\frac{1+\nu }{E}\right)}^{2}{K}_{I}^{2}{a}^{3}\frac{(2{\kappa}^{2}+1)}{3\pi }\underset{\frac{-\pi }{4}}{\overset{\frac{\pi}{4}}{\int}}\frac{1}{\cos{(\theta )}^{3}}d\theta .\)

Soit:

\({\parallel u\parallel }_{{L}^{2}}=\frac{1+\nu }{E}{K}_{I}{a}^{\frac{3}{2}}\sqrt{\frac{(2{\kappa}^{2}+1)}{3\pi }\underset{\frac{-\pi }{4}}{\overset{\frac{\pi}{4}}{\int}}\frac{1}{\cos{(\theta )}^{3}}d\theta }\approx 7,6057690825\times {10}^{-6}{\text{m}}^{2}.\)

1.4.2. Calcul de l’énergie de la structure

L’énergie de la structureest définie par:

\({E}^{e}=\frac{1}{2}{\int}_{\Omega}\sigma :\epsilon \mathit{dS}.\)

Puisque \(\nu =0\) , on a \(\epsilon =\frac{1}{E}\sigma\) . D’où:

\(\sigma :\epsilon =\frac{1}{E}\left({\sigma}_{xx}^{2}+{\sigma}_{yy}^{2}+2{\sigma}_{xy}^{2}\right).\)

En utilisant l’expression du tenseur des contraintes donnés précédemment et près développements, on obtient:

\(\sigma :\epsilon =\frac{{K}_{I}^{2}}{E}\frac{1}{2\pi r}\left(2\cos{\left(\frac{\theta}{2}\right)}^{2}+2\cos{\left(\frac{\theta}{2}\right)}^{2}\sin{\left(\frac{\theta}{2}\right)}^{2}\right).\)

On remarque que:

\(2\cos{\left(\frac{\theta}{2}\right)}^{2}=1+\cos\theta ,\)

\(2\cos{\left(\frac{\theta}{2}\right)}^{2}\sin{\left(\frac{\theta}{2}\right)}^{2}=\frac{1}{2}\sin{(\theta )}^{2}=\frac{1}{4}(1-\cos(2\theta )).\)

On a donc:

\(\sigma :\epsilon =\frac{{K}_{I}^{2}}{E}\frac{1}{2\pi r}\left(\frac{5}{4}+\cos(\theta )-\cos(2\theta )\right).\)

Soit:

\({E}^{e}=\frac{{K}_{I}^{2}}{E}\frac{1}{4\pi }{\int}_{\Omega}\frac{1}{r}\left(\frac{5}{4}+\cos(\theta )-\cos(2\theta )\right)\mathit{dS.}\)

On remarque que pour \(n=-1\) et \(g:\theta \to \frac{5}{4}+\cos(\theta )-\cos(2\theta )\) ., on a:

\(I=\underset{0}{\overset{\rho (\theta )}{\int}}\frac{1}{r}\left(\frac{5}{4}+\cos(\theta )-\cos(2\theta )\right)d\theta .\)

On a dans ce cas:

\(g(\theta )+g\left(\theta +\frac{\pi}{2}\right)+g(\theta +\pi )+g\left(\theta +\frac{3\pi }{2}\right)=5.\)

D’où:

\({\int}_{\Omega}\frac{1}{r}\left(\frac{5}{4}+\cos(\theta )-\cos(2\theta )\right)\mathit{dS}=5a\underset{\frac{-\pi }{4}}{\overset{\frac{\pi}{4}}{\int}}\frac{1}{\cos(\theta )}d\theta =\mathrm{5a}\left(\ln(3)+2\sqrt{2}\right).\)

Et finalement:

\({E}^{e}=\frac{{K}_{I}^{2}}{E}\frac{5a}{4\pi }\left(\ln(3)+2\sqrt{2}\right)\approx 3,50687407712\times {10}^{-6}\text{J}\times {\text{m}}^{-1}.\)

1.5. Références bibliographiques

1. GENIAUT S., MASSIN P.: eXtended Finite Element Method, Manuel de référence de Code_Aster , [R7.02.12]

2. Rice, J. R. (1968), « A path independent integral and the approximate analysis of strain concentration by notches and cracks », Journal of Applied Mechanics 35 : 379–386

3. Satzi, Belytschko: An extended Finite Element Method with Higher-Order Elements for Curved Cracks

4. LABORDE P., POMMIER J., RENARD Y., SALAUN M., “High-order extended finite element

method for cracked domains”, International Journal for Numerical Methods in Engineering, vl.

64, pp. 354-381, 2005