v6.03.171 SSNP171 - Fermeture d’une fissure en flexion#
Résumé
Ce test permet d’évaluer la méthode XFEM avec contact, lors de la fermeture de fissures.
On considère une éprouvette 2D en flexion, comportant une fissure centrale. La flexion de l’éprouvette provoque une fermeture partielle de la fissure. Sur le fond de fissure demeurant en ouverture, on teste le facteur d’intensité de contrainte \({K}_{I}\) . La théorie prévoit que la prise en compte de contact corrige le facteur \({K}_{I}\) d’environ \(8\text{\%}\) .
Pour jauger l’influence du contact lors de la fermeture, on envisage sept modélisations:
Modélisation A: une modélisation 2D et des éléments linéaires, avec interpénétration des lèvres de la fissure.
Modélisation B: une modélisation 2D et des éléments linéaires, avec contact entre les lèvres de la fissure.
Modélisation C: une modélisation 2D et des éléments quadratiques, avec interpénétration des lèvres de la fissure.
Modélisation D: une modélisation 2D et des éléments quadratiques, avec contact entre les lèvres de la fissure.
Modélisation E : une modélisation 3D et des éléments linéaires, avec interpénétration des lèvres de la fissure.
Modélisation F: une modélisation 3D et des éléments linéaires, avec contact entre les lèvres de la fissure.
Modélisation G : une modélisation 3D et des éléments quadratiques, avec interpénétration des lèvres de la fissure.
Solution de référence#
Méthode de calcul#
ModélisationsA, C, E, G#
On teste le facteur d’intensité de contrainte \({K}_{I}\) par rapport à la valeur analytique sans contact.
Pour une éprouvette en flexion et une fissure droite d’abscisse \([-a\phantom{\rule{2em}{0ex}},\phantom{\rule{2em}{0ex}}a]\) , la valeur analytique de \({K}_{I}\) sans fermeture de la fissure, est déterminée par la formule de Bui [ 7 ]:
\({K}_{I}=\frac{1}{\sqrt{\pi a}}\underset{-a}{\overset{a}{\int}}-t\sqrt{\frac{a-t}{t+a}}\mathit{dt}\)
où \(a\) désigne la demi-longueur de la fissure.
A partir de l’intégrale ci-dessus, on calcule l’expression simplifiée de \({K}_{I}\) suivante:
\({K}_{I}=\frac{{a}^{\frac{3}{2}}\sqrt{\pi}}{2}\)
En effet dans l’intégrale, on effectue plusieurs changements de variables:
\(t\to \mathit{at}\) on a \({K}_{I}=\frac{1}{\sqrt{\pi a}}{a}^{2}\underset{-1}{\overset{1}{\int}}-t\sqrt{\frac{1-t}{t+1}}\mathit{dt}\)
\(t\to t-1\) on a \({K}_{I}=\frac{1}{\sqrt{\pi a}}{a}^{2}\underset{0}{\overset{2}{\int}}-(t-1)\sqrt{\frac{2-t}{t}}\mathit{dt}\)
\(t\to \mathrm{2t}\) on a \({K}_{I}=\frac{1}{\sqrt{\pi a}}{\mathrm{2a}}^{2}\underset{0}{\overset{1}{\int}}-(\mathrm{2t}-1)\sqrt{\frac{1-t}{t}}\mathit{dt}\)
\(t\to {t}^{2}\) on a \({K}_{I}=\frac{1}{\sqrt{\pi a}}{\mathrm{2a}}^{2}\underset{0}{\overset{1}{\int}}-2({\mathrm{2t}}^{2}-1)\sqrt{1-{t}^{2}}\mathit{dt}\)
\(t\to \sint\) on a \({K}_{I}=\frac{1}{\sqrt{\pi a}}{\mathrm{2a}}^{2}\underset{0}{\overset{\frac{\pi}{2}}{\int}}\mathrm{2cos}\mathrm{2t}{\cos}^{2}t\mathit{dt}\)
Après développement du produit de \(\cos\) :
\({K}_{I}=\frac{1}{\sqrt{\pi a}}{\mathrm{2a}}^{2}\underset{0}{\overset{\frac{\pi}{2}}{\int}}(\frac{1}{2}+\cos\mathrm{2t}+\frac{1}{2}\cos\mathrm{4t})\mathit{dt}\)
D’où l’expression finale (après l’intégration):
\({K}_{I}=\frac{{a}^{\frac{3}{2}}\sqrt{\pi}}{2}\)
La demi longueur étant unitaire (\(a=1\) ), on teste donc le facteur d’intensité de contrainte:
\({K}_{I}=\frac{\sqrt{\pi}}{2}\approx 0,88629\)
Par construction, il n’existe pas de mode II. On teste donc : \({K}_{\mathit{II}}=0\)
D’autre part, le taux de restitution d’énergie analytique vaut:
\(G=\frac{1-{\nu}^{2}}{E}({K}_{I}^{2}+{K}_{\mathit{II}}^{2})\approx 7,853982E-7\)
Modélisation B, D, F#
Pour une éprouvette en flexion et une fissure droite d’abscisse \([-a\phantom{\rule{2em}{0ex}},\phantom{\rule{2em}{0ex}}a]\) , la valeur analytique de \({K}_{I}\) avec fermeture de la fissure, est déterminée par la formule de Bui [ 7 ]:
\({K}_{I}=\sqrt{\frac{2}{\pi (a-c)}}\underset{c}{\overset{a}{\int}}t\sqrt{\frac{t-c}{a-t}}\mathit{dt}\)
où \(c\) est l’abscisse du point pour lequel la fissure commence à s’ouvrir (voir [])
et \(a\) désigne la demi-longueur de la fissure.
On montre analytiquement [ 7 ] que l’abscisse \(c=-\frac{a}{3}\)
Alors, \({K}_{I}=\sqrt{\frac{3}{2\pi a}}\underset{-a/3}{\overset{a}{\int}}t\sqrt{\frac{t+a/3}{a-t}}\mathit{dt}\)
A partir de l’intégrale ci-dessus, on calcule l’expression simplifiée de \({K}_{I}\) suivante:
\({K}_{I}={(\frac{2}{3}a)}^{\frac{3}{2}}\sqrt{(\pi )}\)
En effet dans l’intégrale, on effectue plusieurs changements de variables:
\(t\to \mathit{at}\) on a \({K}_{I}=\sqrt{\frac{3}{2\pi a}}{a}^{2}\underset{-1/3}{\overset{1}{\int}}t\sqrt{\frac{t+1/3}{1-t}}\mathit{dt}\)
\(t\to t-\frac{1}{3}\) on a \({K}_{I}=\sqrt{\frac{3}{2\pi a}}{a}^{2}\underset{0}{\overset{4/3}{\int}}(t-\frac{1}{3})\sqrt{\frac{t}{4/3-t}}\mathit{dt}\)
\(t\to \frac{4}{3}t\) on a \({K}_{I}=\sqrt{\frac{3}{2\pi a}}{a}^{2}\frac{4}{9}\underset{0}{\overset{1}{\int}}(\mathrm{4t}-1)\sqrt{\frac{t}{1-t}}\mathit{dt}\)
\(t\to {t}^{2}\) on a \({K}_{I}=\sqrt{\frac{3}{2\pi a}}{a}^{2}\frac{4}{9}\underset{0}{\overset{1}{\int}}({\mathrm{4t}}^{2}-1)\frac{{\mathrm{2t}}^{2}}{\sqrt{1-{t}^{2}}}\mathit{dt}\)
\(t\to \sint\) on a \({K}_{I}=\sqrt{\frac{3}{2\pi a}}{a}^{2}\frac{4}{9}\underset{0}{\overset{\pi /2}{\int}}({\mathrm{4sin}}^{2}t-1)({\mathrm{2sin}}^{2}t)\mathit{dt}\)
Après développement des puissances du \(\sint\) :
\({K}_{I}=\sqrt{\frac{3}{2\pi a}}{a}^{2}\frac{4}{9}\underset{0}{\overset{\pi /2}{\int}}(2-\mathrm{3cos}\mathrm{2t}+\cos\mathrm{4t})\mathit{dt}\)
D’où l’expression finale (après l’intégration):
\({K}_{I}={(\frac{2}{3}a)}^{\frac{3}{2}}\sqrt{(\pi )}\)
La demi longueur étant unitaire (\(a=1\) ), on teste donc le facteur d’intensité de contrainte:
\({K}_{I}={(\frac{2}{3})}^{\frac{3}{2}}\sqrt{(\pi )}\approx 0,9648017\)
Par construction, il n’existe pas de mode II. On teste donc : \({K}_{\mathit{II}}=0\)
D’autre part, le taux de restitution d’énergie analytique vaut:
\(G=\frac{1-{\nu}^{2}}{E}({K}_{I}^{2}+{K}_{\mathit{II}}^{2})\approx 9,308423E-7\)
Grandeurs et résultats de référence#
On teste \({K}_{I}\) , \({K}_{\mathit{II}}\) et \(G\) .
Modélisations A, C, E, F, G#
\({K}_{I}=\frac{\sqrt{\pi}}{2}\approx 0,88629\)
\({K}_{\mathit{II}}=0\)
\(G=\frac{1-{\nu}^{2}}{E}({K}_{I}^{2}+{K}_{\mathit{II}}^{2})\approx 7,853982E-7\)
Modélisation B et D#
\({K}_{I}={(\frac{2}{3})}^{\frac{3}{2}}\sqrt{(\pi )}\approx 0,9648017\)
\({K}_{\mathit{II}}=0\)
\(G=\frac{1-{\nu}^{2}}{E}({K}_{I}^{2}+{K}_{\mathit{II}}^{2})\approx 9,308423E-7\)
Incertitude sur la solution#
Faible, solution semi-analytique.
Références bibliographiques#
GENIAUT S., MASSIN P.: eXtended Finite Element Method, Manuel de référence de Code_Aster , [R7.02.12]
GENIAUT S.: Approche XFEM pour la fissuration sous contact de structures industrielles, Thèse de doctorat Laboratoire GeM, 2006.
H.D. Bui: Mécanique de la rupture fragile, Ed. Masson, 1978.
THRESHER R.W., SMITH F.W.: The partially closed Griffith crack, International Journal of Fracture , vol. 9, n°1, pp. 33-41, 1973
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
Modélisation A#
Caractéristiques de la modélisation#
Dans cette modélisation, la méthode des éléments finis étendue (\(\text{X-FEM}\) ) est utilisée. Les éléments finis sont linéaires.
On définit un rayon d’enrichissement éléments X-FEM fond de fissure \({R}_{\mathit{ENRI}}=0.5\) . Ce rayon d’enrichissement permet de capter plus précisément la solution asymptotique singulière en fond de fissure. La taille de rayon relativement large (25% de la longueur de la fissure) n’introduit pas des problème de conditionnement, compte-tenu de la nouvelle approximation en fond de fissure [R7.02.12].
Caractéristiques du maillage#
On raffine le maillage au centre de l’éprouvette, afin d’optimiser le calcul de la solution en déplacement au voisinage de la fissure [].
NOMBRE DE NOEUDS 4088
NOMBRE DE MAILLES 4360
SEG2 80
TRIA3 466
QUAD4 3814
Figure 3.2-1: Maillage du domaine
Grandeurs testées et résultats#
Sur le fond de fissure en ouverture \(P2=(a,0)\) , on teste le facteur d’intensité de contraintes \({K}_{I}\) donné par la commande CALC_G, par rapport à la valeur analytique explicitée au paragraphe [ 5 ].
Pour la méthode \(G-\mathit{thêta}\) (commande CALC_G), on choisit les couronnes du champ thêta suivantes:
Couronne 1 |
Couronne 2 |
Couronne 3 |
Couronne 4 |
Couronne 5 |
Couronne 6 |
|
Rinf |
0,1 |
0,2 |
0,3 |
0,1 |
0,1 |
0,2 |
Rsup |
0,2 |
0,3 |
0,4 |
0,3 |
0,4 |
0,4 |
Identification |
Type de référence |
Valeur de r éférence |
Précision |
CALC_G/K1 |
“ANALYTIQUE” |
0,88629 |
1,0% |
CALC_G/K2 |
“ANALYTIQUE” |
0.00 |
0,001% |
CALC_G/G |
“ANALYTIQUE” |
7.85514E-07 |
2,0% |
Résultats complémentaires#
Voici le champ de déplacement calculé par Aster, sans activer de l’algorithme de contact entre les lèvres de fissure:
Figure 3.4-1: Champ de déplacement (avec offset)
Modélisation B#
Caractéristiques de la modélisation#
Dans cette modélisation, la méthode des éléments finis étendue (\(\text{X-FEM}\) ) est utilisée. On définit un rayon d’enrichissement \({R}_{\mathit{ENRI}}=0.5\) . Les éléments finis sont linéaires.
Caractéristiques du maillage#
Identique à la modélisation A.
Grandeurs testées et résultats#
Sur le fond de fissure en ouverture \(P2=(a,0)\) , on teste le facteur d’intensité de contraintes \({K}_{I}\) donné par la commande CALC_G, par rapport à la valeur analytique explicitée au paragraphe [ 5 ].
Pour la méthode \(G-\mathit{thêta}\) (commande CALC_G), on choisit les couronnes du champ thêta suivantes:
Couronne 1 |
Couronne 2 |
Couronne 3 |
Couronne 4 |
Couronne 5 |
Couronne 6 |
|
Rinf |
0,1 |
0,2 |
0,3 |
0,1 |
0,1 |
0,2 |
Rsup |
0,2 |
0,3 |
0,4 |
0,3 |
0,4 |
0,4 |
Identification |
Type de référence |
Valeur de r éférence |
Précision |
CALC_G/K1 |
“ANALYTIQUE” |
0,9648 |
1,0% |
CALC_G/K2 |
“ANALYTIQUE” |
0.00 |
0,001 |
CALC_G/G |
“ANALYTIQUE” |
9,3084E-7 |
2,0% |
Résultats complémentaires#
En activant l’algorithme de contact, on obtient la solution en déplacement au voisinage de la fissure donnée à la [] (à comparer avec le post-traitement ).
D’autre part, on retrouve une conclusion analytique, à savoir l’abscisse du point pour lequel la fissure commence à s’ouvrir:
\(c=-\frac{a}{3}\approx -0.33\)
En effet, en post-traitement, nous avons analysé les pressions de contact pour déterminer l’abscisse du dernier point de contact . Graphiquement, on estime que la pression chute fortement autour du point d’abscisse \({x}_{c}\approx -0.325\) , ce qui valide l’hypothèse d’estimation de \({K}_{I}\) .
Figure 4.4-1: faciès de fermeture de la fissure avec contact
Figure 4.4-2: pression de contact le long de la fissure
Modélisation C#
Caractéristiques de la modélisation#
Identique à la modélisation A avec des éléments quadratiques.
Caractéristiques du maillage#
NOMBRE DE NOEUDS 12455
NOMBRE DE MAILLES 4360
SEG3 80
TRIA6 466
QUAD8 3814
Grandeurs testées et résultats#
Pour la méthode \(G-\mathit{thêta}\) (commande CALC_G), on choisit les couronnes du champ thêta suivantes:
Couronne 1 |
Couronne 2 |
Couronne 3 |
Couronne 4 |
Couronne 5 |
Couronne 6 |
|
Rinf |
0,1 |
0,2 |
0,3 |
0,1 |
0,1 |
0,2 |
Rsup |
0,2 |
0,3 |
0,4 |
0,3 |
0,4 |
0,4 |
Identification |
Type de référence |
Valeur de r éférence |
Précision |
CALC_G/K1 |
“ANALYTIQUE” |
0,88629 |
0,05% |
CALC_G/K2 |
“ANALYTIQUE” |
0.00 |
0,001 |
CALC_G/G |
“ANALYTIQUE” |
7.85514E-07 |
1,0% |
Commentaires#
Les résultats sont plus précis que la modélisation A.
Modélisation D#
Caractéristiques de la modélisation#
Identique à la modélisation B avec des éléments quadratiques.
Caractéristiques du maillage#
Identique à la modélisation C.
Grandeurs testées et résultats#
Pour la méthode \(G-\mathit{thêta}\) (commande CALC_G), on choisit les couronnes du champ thêta suivantes:
Couronne 1 |
Couronne 2 |
Couronne 3 |
Couronne 4 |
Couronne 5 |
Couronne 6 |
|
Rinf |
0,1 |
0,2 |
0,3 |
0,1 |
0,1 |
0,2 |
Rsup |
0,2 |
0,3 |
0,4 |
0,3 |
0,4 |
0,4 |
Identification |
Type de référence |
Valeur de r éférence |
Précision |
CALC_G/K1 |
“ANALYTIQUE” |
0,9648 |
1,0% |
CALC_G/K2 |
“ANALYTIQUE” |
0.00 |
0,01 |
CALC_G/G |
“ANALYTIQUE” |
9,3084E-7 |
2,0% |
Commentaires#
Avec les éléments quadratiques on ne note pas de gain en activant l’algorithme de contact par rapport aux éléments linéaires (modélisation B). La singularité au point de décollement pourrait engendrer une piètre convergence de l’algorithme de contact (l’espace des multiplicateurs capte difficilement la singularité du champ de contrainte au point de décollement, cf. ), ce qui engendrerait une difficulté à capter la solution attendue pour ce problème physique.
Néanmoins les résultats sont très satisfaisants étant donné le faible raffinement du maillage (moins de 1500 éléments).
Modélisation E#
Caractéristiques de la modélisation#
Extrusion de la modélisation A avec des éléments linéaires.
Caractéristiques du maillage#
NOMBRE DE NOEUDS 3868
NOMBRE DE MAILLES 3240
QUAD4 120
PENTA6 564
HEXA8 2556
Figure 7.2-1: Maillage avec zone centrale raffinée en HEXA20
Grandeurs testées et résultats#
Pour la méthode \(G-\mathit{thêta}\) (commande CALC_G), on choisit les couronnes du champ thêta suivantes:
Couronne 1 |
Couronne 2 |
Couronne 3 |
Couronne 4 |
Couronne 5 |
Couronne 6 |
|
Rinf |
0,1 |
0,2 |
0,3 |
0,1 |
0,1 |
0,2 |
Rsup |
0,2 |
0,3 |
0,4 |
0,3 |
0,4 |
0,4 |
Identification |
Type de référence |
Valeur de r éférence |
Précision |
CALC_G/K1 (MAX) |
“ANALYTIQUE” |
0,88629 |
2, 0% |
CALC_G/K1 (MIN) |
“ANALYTIQUE” |
0,88629 |
2, 0% |
CALC_G/K2 (MAX) |
“ANALYTIQUE” |
0.00 |
0,004 |
CALC_G/K2 (MIN) |
“ANALYTIQUE” |
0.00 |
0,004 |
CALC_G/G (MAX) |
“ANALYTIQUE” |
7.85514E-07 |
4,0% |
CALC_G/G (MIN) |
“ANALYTIQUE” |
7.85514E-07 |
4,0% |
Commentaires#
En 3D, on constate une légère dégradation de la précision pour le calcul des facteurs par rapport à une modélisation 2D équivalente (modélisation A).
Modélisation F#
Extrusion de la modélisation B avec des éléments linéaires.
Caractéristiques du maillage#
Identique à la modélisation E.
Grandeurs testées et résultats#
Pour la méthode \(G-\mathit{thêta}\) (commande CALC_G), on choisit les couronnes du champ thêta suivantes:
Couronne 1 |
Couronne 2 |
Couronne 3 |
Couronne 4 |
Couronne 5 |
Couronne 6 |
|
Rinf |
0,1 |
0,2 |
0,3 |
0,1 |
0,1 |
0,2 |
Rsup |
0,2 |
0,3 |
0,4 |
0,3 |
0,4 |
0,4 |
Identification |
Type de référence |
Valeur de r éférence |
Précision |
CALC_G/K1 (MAX) |
“ANALYTIQUE” |
0,9648 |
2, 0% |
CALC_G/K1 (MIN) |
“ANALYTIQUE” |
0,9648 |
2, 0% |
CALC_G/K2 (MAX) |
“ANALYTIQUE” |
0.00 |
0,004 |
CALC_G/K2 (MIN) |
“ANALYTIQUE” |
0.00 |
0,004 |
CALC_G/G (MAX) |
“ANALYTIQUE” |
9,3084E-7 |
4,0% |
CALC_G/G (MIN) |
“ANALYTIQUE” |
9,3084E-7 |
4,0% |
Commentaires#
En 3D, on constate une légère dégradation de la précision pour le calcul des facteurs par rapport à une modélisation 2D équivalente (modélisation B).
Modélisation G#
Extrusion de la modélisation A avec des éléments quadratiques.
Caractéristiques du maillage#
NOMBRE DE NOEUDS 14793
NOMBRE DE MAILLES 3240
QUAD8 120
PENTA15 564
HEXA20 2556
Grandeurs testées et résultats#
Pour la méthode \(G-\mathit{thêta}\) (commande CALC_G), on choisit les couronnes du champ thêta suivantes:
Couronne 1 |
Couronne 2 |
Couronne 3 |
Couronne 4 |
Couronne 5 |
Couronne 6 |
|
Rinf |
0,1 |
0,2 |
0,3 |
0,1 |
0,1 |
0,2 |
Rsup |
0,2 |
0,3 |
0,4 |
0,3 |
0,4 |
0,4 |
Identification |
Type de référence |
Valeur de r éférence |
Précision |
CALC_G/K1 (MAX) |
“ANALYTIQUE” |
0,88629 |
0,5% |
CALC_G/K1 (MIN) |
“ANALYTIQUE” |
0,88629 |
0,5% |
CALC_G/K2 (MAX) |
“ANALYTIQUE” |
0.00 |
0,001 |
CALC_G/K2 (MIN) |
“ANALYTIQUE” |
0.00 |
0,001 |
CALC_G/G (MAX) |
“ANALYTIQUE” |
7.85514E-07 |
1,0% |
CALC_G/G (MIN) |
“ANALYTIQUE” |
7.85514E-07 |
1,0% |
Commentaires#
En 3D, on constate une amélioration des résultats par rapport à la modélisation avec les éléments linéaires (modélisation E).
Synthèses des résultats#
Lors de la fermeture d’une fissure, la prise en compte du contact entre les lèvres affecte significativement le facteur d’intensité de contrainte \({K}_{I}\) , par rapport à la méthode XFEM sans contact.
Dans ce cas-test, on retrouve une différence de l’ordre de \(8\text{\%}\) sur la valeur de \({K}_{I}\) calculée par Aster, en passant de la modélisation A (sans contact entre les lèvres) à la modélisation B (avec contact). Cette différence est en accord avec les prévisions théoriques.
En effet, la théorie prévoit que la différence vaut:
\(\frac{{K}_{I}^{\text{avec contact}}-{K}_{I}^{\text{sans contact}}}{{K}_{I}^{\text{avec contact}}}=\frac{({(2/3)}^{\frac{3}{2}}-1/2)}{{(2/3)}^{\frac{3}{2}}}\approx 8,14413\text{\%}\)
Les résultats obtenus montrent que la valeur de \({K}_{I}\) avec prise en compte du contact est bien correctement évaluée et supérieure à celle obtenue sans prise en compte du contact. Ce résultat souligne l’importance de la prise en compte des effets liés au contact dans les simulations numériques de fissuration, car cet exemple illustre un cas où la solution avec interpénétration (sans contact) n’est pas conservative du point de vue des facteurs d’intensité de contraintes. Cela peut s’avérer non sans conséquence pour l’étude de la propagation, basée sur des critères de propagation écrits en fonction des facteurs d’intensité de contraintes.