Figure 1.1-1: géométrie et fond de fissure elliptique
v3.04.110 SSLV110 - Fissure elliptique dans un milieu infini#
Résumé:
Il s’agit d’un test en statique pour un problème tridimensionnel. Ce test permet de calculer le taux de restitution d’énergie local sur le fond de fissure par la méthode thêta (commande CALC_G en FEM et CALC_G_XFEM en X-FEM).
Le rayon des couronnes d’intégration est variable le long de la fissure, et le taux de restitution d’énergie local est calculé suivant deux méthodes différentes (LEGENDRE et LINEAIRE/LAGRANGE). Les couronnes d’intégration sont définies de deux manières: par la donnée d’un rayon variable (R_INF_FO et R_SUP_FO) et la donnée du nombre de couches (NB_COUCHE_INF et NB_COUCHE_SUP).
L’intérêt du test est la validation de la méthode thêta en \(\mathrm{3D}\) et des points suivants:
comparaison des résultats avec une solution analytique,
stabilité des résultats vis à vis des couronnes d’intégration,
comparaison entre deux méthodes différentes pour le calcul de \(G\) local,
2 cas de chargements équivalents (pression répartie et chargement volumique).
Calcul de \(\mathrm{K1}\) à partir de \(G\) et de la formule d’Irwin
Ce test contient 3 modélisations différentes (A,F,G).
La modélisation F teste le calcul de \(\mathrm{K1}\) pour une fissure non maillée (méthode X-FEM). Elle permet aussi de comparer les erreurs commises sur le calcul de \(\mathit{K1}\) avec l’opérateur POST_K1_K2_K3 ou l’opérateur CALC_G_XFEM.
La modélisation G valide le calcul du facteur d’intensité des contrainte équivalentes en présence de zones cohésives (voir documentation [R7.02.18]), par l’opérateur CALC_G_XFEM.
Solution de référence#
Méthode de calcul utilisée pour la solution de référence#
La solution de référence est une solution analytique issue de SIH [bib1] et [bib2].
On note que l’angle \(\alpha\) désigne ici l’angle paramétrique du point \(M\) (angle par rapport à l’axe \(\mathit{Ox}\) du projeté de \(M\) sur le cercle de rayon \(b\) ) et non la coordonnée polaire de ce point.
Ici: \(a=25\mathit{mm}\) et \(b=6\mathit{mm}\) , donc \(k=0,9707728\)
Les valeurs de l’intégrale elliptique \(E(k)\) sont tabulées dans [bib3], en fonction de \(\mathit{asin}(k)\) qui vaut ici \(76,11°\) . On trouve alors: \(E(k)=1,0672\) .
D’où le facteur d’intensité des contraintes en \(\mathit{MPa.}\sqrt{}\mathit{mm}\) : \({K}_{I}(\alpha )=4.0680{\left[{\sin}^{2}\alpha +\frac{{b}^{2}}{{a}^{2}}{\cos}^{2}\alpha \right]}^{1/4}\)
Puis, à partir de la formule d’Irwin (déformation plane):
\(g(\alpha )=\frac{1-{\nu}^{2}}{E}{K}_{I}{(\alpha )}^{2}\)
Bibliographie#
G.C. SIH : Mathematical Theories of Brittle Fracture - FRACTURE, volII‑Academic Press - 1968
M.K. KASSIN et G.C. SIH : Three-dimensional stress distribution around an elliptical crack under arbitrary loadings J. Appl. Mech., 88, 601-611, 1966.
TADA, P.PARIS, G. IRWIN : The Stress Analysis of Cracks Handbook - Third Edition - ASM International - 2000
Modélisation A#
Caractéristiques de la modélisation#
\(A=\mathit{N01099}\) (\(s=0.\) )
\(B=\mathrm{N01259}\) (\(s=26.68\) )
\(C=\mathrm{N01179}\) (\(s=18,0\) ; \(\alpha =0,77\approx \pi /4\) )
Chargement : Pression unitaire répartie sur la face du bloc opposée au plan de la lèvre :
\(P=\mathrm{1.MPa}\) dans le plan \(Z=\mathrm{1250.mm}\) .
Caractéristiques du maillage#
Nombre de nœuds : 20567
Nombre de mailles et types : 160PENTA15, 175 PYRAM13, 9383 TETRA10 et 1290HEXA20
Grandeurs testées et résultats#
Les valeurs testées sont :
le taux de restitution d’énergie local \(g\) en tous les nœuds du fond de fissure (option G de CALC_G). Ces valeurs sont testées pour des couronnes d’intégration définies par des rayons variables (R_INF_FO et R_SUP_FO) ou par le nombre de couches (NB_COUCHE_INF et NB_COUCHE_SUP).
le taux de restitution d’énergie local \(g\) pour les points A et B du fond de fissure avec l’opérateur de calcul POST_K1_K2_K3.
le facteur d’intensité de contraintes \({K}_{I}\) , calculé aux mêmes nœuds du fond de fissure que le taux de restitution d’énergie (options K et KJ de CALC_G). Dans le cas de l’option KJ, on met aussi en évidence que lorsque la valeur de G calculée est négative, on a bien \({K}_{J}=0\) (par convention).
le taux de restitution d’énergie \({G}_{\mathit{IRWIN}}\) calculé à partir de \({K}_{I}\) et de la formule d’Irwin
Dans le cas du taux de restitution \({G}_{\mathit{IRWIN}}\) , on s’assure que lorsque la valeur de G calculée
Le maillage ne comprend qu’une des lèvres de la fissure, il faut donc utiliser le mot-clé ’SYME’ pour multiplier automatiquement par 2 dans le calcul Aster le taux de restitution de l’énergie calculé par extension virtuelle de la lèvre unique.
Identification |
Référence |
\(\text{\%}\) tolérance |
\(g(A)\) couronne \({C}_{1}\) (‘LEGENDRE’) |
7.171 10-5 |
1 |
\(g(A)\) couronne \({C}_{2}\) (‘LEGENDRE’) |
7.171 10-5 |
1 |
\(g(A)\) couronne \({C}_{1}\) (‘ LINEAIRE ’) |
7.171 10-5 |
1 |
\(g(A)\) couronne \({C}_{2}\) (‘ LINEAIRE ’) |
7.171 10-5 |
1 |
\({G}_{\mathit{IRWIN}}(A)\) couronne \({C}_{2}\) (‘ LEGENDRE ’) |
7.171 10-5 |
1 |
\({G}_{\mathit{IRWIN}}(A)\) couronne \({C}_{2}\) (‘ LINEAIRE ’) |
7.171 10-5 |
1 |
\({G}_{\mathit{IRWIN}}(A)\) POST_K1_K2_K3 |
7.171 10-5 |
4 |
\(g(B)\) couronne \({C}_{1}\) (‘LEGENDRE’) |
1.721 10-5 |
10 |
\(g(B)\) couronne \({C}_{2}\) (‘LEGENDRE’) |
1.721 10-5 |
10 |
\(g(B)\) couronne \({C}_{1}\) (‘ LINEAIRE ’) |
1.721 10-5 |
1 |
\(g(B)\) couronne \({C}_{2}\) (‘ LINEAIRE ’) |
1.721 10-5 |
1 |
\({G}_{\mathit{IRWIN}}(B)\) couronne \({C}_{2}\) (‘ LEGENDRE ’) |
1.721 10-5 |
4 |
\({G}_{\mathit{IRWIN}}(B)\) couronne \({C}_{2}\) (‘ LINEAIRE ’) |
1.721 10-5 |
6 |
\({G}_{\mathit{IRWIN}}(B)\) POST_K1_K2_K3 |
1.721 10-5 |
3 |
\(g(C)\) couronne \({C}_{1}\) (‘LEGENDRE’) |
5.215 10-5 |
1 |
\(g(C)\) couronne \({C}_{2}\) (‘LEGENDRE’) |
5.215 10-5 |
1 |
\(g(C)\) couronne \({C}_{1}\) (‘ LINEAIRE ’) |
5.215 10-5 |
1 |
\(g(C)\) couronne \({C}_{2}\) (‘ LINEAIRE ’) |
5.215 10-5 |
1 |
Identification |
Référence |
\(\text{\%}\) tolérance |
\({K}_{J}(A)\) couronne \({C}_{1}\) (‘LEGENDRE’) |
4.067 |
1 |
\({K}_{J}(A)\) couronne \({C}_{2}\) (‘LEGENDRE’) |
4.067 |
1 |
\({K}_{I}(A)\) couronne \({C}_{2}\) (‘ LINEAIRE ’) |
4.067 |
1 |
\({K}_{I}(A)\) couronne \({C}_{2}\) (‘ LEGENDRE ’) |
4.067 |
1 |
\({K}_{J}(B)\) couronne \({C}_{1}\) (‘LEGENDRE’) |
1.992 |
10 |
\({K}_{J}(B)\) couronne \({C}_{2}\) (‘LEGENDRE’) |
1.992 |
10 |
\({K}_{I}(B)\) couronne \({C}_{2}\) (‘LEGENDRE’) |
1.992 |
2 |
\({K}_{I}(B)\) couronne \({C}_{2}\) (‘ LINEAIRE ’) |
1.992 |
5 |
\({K}_{J}(C)\) couronne \({C}_{1}\) (‘LEGENDRE’) |
3.469 |
1 |
\({K}_{J}(C)\) couronne \({C}_{2}\) (‘LEGENDRE’) |
3.469 |
1 |
\({K}_{J}(A)\) couronne \({C}_{2}\) (‘ LINEAIRE ’ , ‘ NB_POINTS_FOND ) |
0.0 |
Modélisation F#
Caractéristiques de la modélisation#
Dans cette modélisation, la fissure n’est pas maillée. On utilise la méthode X-FEM.
Compte tenu des symétries du problème, il est possible de ne représenter qu’un huitième de la structure (comme cela est fait dans la modélisation A). Cependant, avec la méthode X-FEM, il n’est pas possible de représenter une fissure qui se situe dans un plan de symétrie (sur la bord du domaine modélisé). On modélise donc dans cette modélisation un quart de la structure, c’est-à-dire une portion de \(90°\) de l’ellipse.
Le maillage est composé de mailles HEXA8, uniformément réparties suivant les axes \(X\) et \(Y\) et réparties en progression géométrique suivant l’axe \(Z\) de manière à ce que dans le plan \(Z=0\) , les mailles sont environ des cubes de coté \(10\mathit{mm}\) .
Des conditions de symétrie sont appliquées sur les faces en \(X=0\) et \(Y=0\) . Le mode rigide suivant l’axe \(Z\) est bloqué en bloquant le déplacement suivant \(Z\) du point situé en \((0,0,-1250\mathit{mm})\) .
Chargement : Pression unitaire répartie sur les deux faces normales du bloc :
\(P=1\mathit{MPa}\) dans les plan \(Z=\pm 1250\mathit{mm}\) .
Caractéristiques du maillage#
Nombre de nœuds : 21000
Nombre de mailles et types : 13000 PENTA6 et 12500 HEXA8 (maillage linéaire)
Figure 4.2-1: maillage initial, vue d’ensemble
Figure 4.2-2: maillage initial, zoom dans le plan milieu
Comme ce maillage initial est bien trop grossier pour un calcul précis des facteurs d’intensité de contraintes le long du fond du fissure, une procédure de raffinement automatique des mailles proches du fond de fissure est utilisée, telle que préconisée dans la documentation [U2.05.02].
La taille cible des mailles souhaité est \(b/9\) . Cela va induire 5 raffinements successifs. La taille des mailles du maillage ainsi raffiné est alors \(h=0,39\mathit{mm}\) .
Le maillage raffiné (celui sur lequel le calcul mécanique est effectué) a pour caractéristiques:
26484 noeuds
7720 TETRA4, 10650 PYRAM5 et 20080 HEXA8
Ce maillage induit 99 points le long du fond de fissure et compte-tenu des conditions de blocage 118404 équations dans le système à résoudre
Figure 4.2-3: maillage raffiné, zoom sur la zone proche de la fissure
Grandeurs testées et résultats#
Les valeurs testées sont les facteurs d’intensité des contraintes \(\mathrm{K1}\) le long du fond de fissure, calculés soit par CALC_G_XFEM, soit par POST_K1_K2_K3.
Pour CALC_G_XFEM, la couronne d’intégration vaut \(2h-4h\) . Le lissage par défaut (Legendre) est utilisé.
Pour POST_K1_K2_K3, l’abscisse curviligne maximale vaut \(4h\) . Afin de réduire les temps de calcul de POST_K1_K2_K3, on ne post-traite que sur 20 points répartis uniformément le long du fond de fissure.
Notons que le temps de calcul pour le lissage Legendre de CALC_G est insensible à ce nombre.
On teste les valeurs aux points \(A\) (\(s=0\) ) et \(B\) (\(s=26,7\) ).
Identification |
Type de référence |
Valeur de référence |
Tolérance |
CALC_G_ XFEM: \(\mathit{K1}(A)\) |
“ANALYTIQUE” |
0.000 |
2.0% |
CALC_G_ XFEM : \(\mathit{K1}(B)\) |
“ANALYTIQUE” |
4.068 |
2.0% |
POST_K1_K2_K3 : \(\mathit{K1}(A)\) |
“ANALYTIQUE” |
0.000 |
6.0% |
POST_K1_K2_K3 : \(\mathit{K1}(B)\) |
“ANALYTIQUE” |
4.068 |
6.0% |
Pour l’opérateur CALC_G_XFEM, les lissages de type LAGRANGE ne permettent pas d’avoir des résultats facilement exploitables; un lissage de type LEGENDRE est donc à privilégier.
Modélisation G#
Caractéristiques de la modélisation#
La géométrie et le chargement sont identiques à la modélisation F: un quart de la structure est modélisé, une pression est appliquée sur la face supérieure. Dans cette modélisation, la fissure initiale est maillée. Par rapport à la modélisation F, on introduit des zones cohésives dans le prolongement de la fissure. Ce prolongement est représenté par des level-sets, de sorte que la discontinuité est prise en compte par une modélisation XFEM, comme pour la modélisation F. La loi cohésive CZM_LIN_MIX est introduite dans ce modèle XFEM par la commande DEFI_CONTACT.
Chargement : Pression unitaire répartie sur la face du bloc opposée au plan de la lèvre :
\(P=1\mathrm{MPa}\) dans le plan \(Z=1250\mathrm{mm}\) .
Les paramètres cohésifs sont choisis pour que ceci ait pour effet d’ouvrir quelques éléments cohésifs au voisinage du fond de fissure initiale:
Pour ne pas avoir une rupture complète, mais simplement une dé-cohésion proche de la pointe de fissure initiale, on prend \({G}_{c}>{G}_{max}\) , avec \({G}_{max}\) le \(G\) local maximal le long front, tout en conservant le même ordre de grandeur pour les deux valeurs. Dans notre cas, \({G}_{c}=2.5\times {10}^{-4}{\mathit{N.mm}}^{-1}\) contre \({G}_{max}=7.2\times {10}^{-5}{\mathit{N.mm}}^{-1}\) .
Pour tout de même observer une dé-cohésion au voisinage de la pointe, la taille caractéristique de la zone cohésive \({l}_{c}=\frac{E{G}_{c}}{(1-{\nu}^{2}){\sigma}_{c}^{2}}\) est choisie de sorte à couvrir quelques éléments tout en restant petite devant la taille de la structure \(h\le {l}_{c}\le a\) . Dans ce cas test, pour réduire le temps de calcul, on a pris \({l}_{c}=14\mathit{mm}\) , ce qui conduit à \({\sigma}_{c}=2\mathit{MPa}\) , à comparer avec des tailles typiques d’éléments \(h=1\mathit{mm}\) sur le petit côté de l’ellipse, et \(h=2\mathit{mm}\) sur le grand côté.
Caractéristiques du maillage#
Nombre de nœuds : 4522
Nombre de mailles et types : 22300 TETRA4
Grandeurs testées et résultats#
Les valeurs testées sont les facteurs d’intensité des contraintes équivalents \(\mathrm{K1}\) le long du fond de fissure, calculés par CALC_G. Afin d’obtenir un résultat régulier, on ne post-traite que sur 5 points répartis uniformément le long du fond de fissure (NB_POINT_FOND=5). On teste les valeurs aux points extrémités du front \(A\) (\(s=0\) ) et \(B\) (\(s=26,7\) ).
Le lissage “LAGRANGE’est utilisé pour \(\theta\) et “LAGRANGE_NO_NO’pour \(G\) .
Identification |
Type de référence |
Valeur de référence |
Tolérance |
CALC_G_ XFEM: \(\mathit{K1}(A)\) |
“ANALYTIQUE” |
0.000 |
4,0% |
CALC_G_ XFEM : \(\mathrm{K1}(B)\) |
“ANALYTIQUE” |
4.068 |
8,0% |
Remarque#
Pour ce test, il y a un écart de quelques pour-cents. La précision peut être améliorée en choisissant une zone cohésive de taille inférieure et en raffinant encore le maillage. Ceci n’a pas été fait ici pour que la modélisation puisse tourner en moins d’une minute.
Synthèse des résultats#
Calcul de \(G\) ou de \(K\) local :
pour une fissure maillée, les 2 méthodes (LEGENDRE et LINEAIRE/LAGRANGE) donnent sensiblement les mêmes résultats (moins de \(\text{5 \%}\) d’erreur par rapport à la solution analytique) sauf au point \(B\) (point extrémité de l’ellipse sur le grand axe) où la méthode Lagrange est la plus précise;
cas de charge : les valeurs obtenues avec le chargement volumique sont légèrement supérieures à celles obtenues avec contraintes imposées (y compris pour les valeurs de \(G\) ). Les différences sont minimes et dues aux intégrations numériques différentes sur le terme de volume et le terme de bord;
la méthode \(\text{X-FEM}\) permet d’évaluer les facteurs d’intensité des contraintes \(K\) sur un maillage non fissuré avec une erreur inférieure à \(\text{10\%}\) .