r7.02.08 Calcul des facteurs d’intensité des contraintes par extrapolation du champ de déplacements#
Résumé :
On décrit ici une méthode de calcul de \(\mathit{K1}\) , \(\mathit{K2}\) et \(\mathit{K3}\) en 2D (plan et axisymétrique) et 3D par extrapolation des sauts de déplacements sur les lèvres de la fissure. Elle est utilisable à l’aide de la commande POST_K1_K2_K3, aussi bien pour une fissure maillée (éléments finis classiques) que pour une fissure non maillée (éléments finis enrichis: méthode X-FEM).
Si la fissure est maillée, elle doit nécessairement être plane; si la fissure n’est pas maillée (méthode X-FEM), elle peut être non plane (mais suffisamment régulière). Dans les deux cas, la méthode n’est applicable que pour des matériaux élastiques linéaires, homogènes et isotropes.
La méthode utilisée est théoriquement moins précise que le calcul à partir de la forme bilinéaire du taux de restitution de l’énergie et des déplacements singuliers [R7.02.01 et R7.02.05] (opérateur CALC_G). Elle permet cependant d’obtenir facilement des valeurs relativement fiables des facteurs d’intensité des contraintes. La comparaison des différentes méthodes de calcul est utile pour estimer la précision des résultats obtenus.
La précision des résultats de la méthode d’extrapolation des sauts de déplacement est nettement améliorée si le maillage est quadratique. Pour une fissure maillée, il est recommandé d’utiliser des éléments dits de «Barsoum» en fond de fissure (éléments dont les nœuds milieux sont situés au quart des arêtes). Pour une fissure non maillée, il est recommandé d’enrichir plusieurs couches d’éléments autour du fond de fissure.
Mise en œuvre des méthodes d’extrapolation#
Les méthodes d’extrapolation des déplacements sont mises en œuvre dans l’opérateur POST_K1_K2_K3, à partir du champ de déplacement calculé sur toute la structure. Les définitions des facteurs d’intensité des contraintes ne sont vraies qu’asymptotiquement; l’extrapolation se fait donc en se restreignant au voisinage du fond de fissure limité par une distance maximale \(\mathit{dmax}\) au fond . \(\mathit{dmax}\) est le paramètre ABSC_CURV_MAXI de l’opérateur. Dans le cas d’une fissure maillée ABSC_CURV_MAXI est facultatif. S’il n’est pas noté, \(\mathit{dmax}\) est calculé automatiquement dans POST_K1_K2_K3 et vaut quatre fois la taille maximale des mailles connectées aux nœuds du fond.
Le principe général du calcul est le suivant:
- Boucle sur le nœuds du fond de fissure (point courant: \(M\) )
Définition du plan \(\Gamma\) normal à la fissure et au fond de fissure, au point \(M\) (plan de normale \(\text{t}\) ) Identification des nœuds des deux lèvres qui appartiennent à \(\Gamma\) : \({P}_{i}^{\sup}\) et \({P}_{i}^{\inf}\)
- Boucle sur ces nœuds:
Si \({r}_{i}^{\sup}=\parallel {\text{MP}}_{i}^{\sup}\parallel \text{dmax}\) : extraction du déplacement en \({P}_{i}^{\sup}\) Si \({r}_{i}^{\inf}=\parallel {\text{MP}}_{i}^{\inf}\parallel \text{dmax}\) : extraction du déplacement en \({P}_{i}^{\inf}\)
Calcul du saut de déplacement dans les trois directions
Extrapolation du saut de déplacement
Trois méthodes d’extrapolation sont programmées. Elles sont illustrées dans ce paragraphe pour une fissure maillée (maillage quadratique), avec ou sans éléments de type «Barsoum». Les éléments de Barsoum sont tels que les nœuds non sommets sur les côtés des éléments quadratiques touchant le fond de fissure sont déplacés au quart du coté [bib4]. Ils permettent de mieux capter la singularité du champ de contraintes en fond de fissure
Élément fini classique (M étant un nœud du fond de fissure) |
Élément fini de type Barsoum |
Méthode 1: on calcule le saut du champ de déplacements au carré et on le divise par \(r\) . Différentes valeurs de \({K}^{2}\) sont obtenues (à un facteur multiplicatif près) par extrapolation en \(r=0\) des segments de droites ainsi obtenus. Si la solution était parfaite (champ asymptotique analytique partout), on devrait obtenir une droite. En réalité, on obtient presque une droite avec un maillage de type «Barsoum», et une courbe non droite sinon:
Méthode 2: on trace le saut du champ de déplacements au carré en fonction de \(r\) . Les approximations de \(K\) sont (toujours à un facteur multiplicatif près) égales à la racine de la pente des segments reliant l’origine aux différents points de la courbe.
Méthode 3 : on identifie le facteur d’intensité de contrainte \(K\) à partir du saut de déplacement \([U]\) par une méthode des moindres carrés. Le recalage se fait sur un segment de longueur \(\mathit{dmax}\) , où \(\mathit{dmax}\) est le paramètre fixé dans l’opérande ABSC_CURV_MAXI de l’opérateur POST_K1_K2_K3 ou dans le cas d’une fissure maillée, si ABSC_CURV_MAXI n’est pas indiquée, \(\mathit{dmax}\) vaut quatre fois la taille maximale des mailles connectées aux nœuds du fond:
\(K\) minimise \(J(k)=\frac{1}{2}\underset{0}{\overset{\text{dmax}}{\int}}{([U(r)]-k\sqrt{r})}^{2}\mathrm{dr}\)
Soit donc la formule explicite pour calculer \(K\) :
\(K=\frac{2}{{r}_{m}^{2}}\underset{0}{\overset{\text{dmax}}{\int}}[U(r)]\sqrt{r}\mathrm{dr}=\frac{1}{{r}_{m}^{2}}\sum_{i=0}^{\mathrm{nbno}-1}({r}_{i+1}-{r}_{i})({[U]}_{i+1}\sqrt{{r}_{i+1}}-{[U]}_{i}\sqrt{{r}_{i}})\)
où \(\mathit{nbno}\) est le nombre de nœuds sur le segment de recalage \([0,\mathit{dmax}]\) . On remarque que dans cette expression \(K\) est, pour un \(\mathit{dmax}\) fixé, une forme linéaire du champ de déplacement.
Précision des méthodes proposées#
La méthode d’extrapolation des sauts de déplacement a été validée sur des tests dont les solutions analytiques sont connues. On en présente ci-dessous certains résultats, en 2D et en 3D, pour une fissure maillée ou non. On compare également les résultats à la méthode théoriquement plus précise fondée sur le calcul du taux de restitution de l’énergie et sur les fonctions singulières (opérateur CALC_G: méthode thêta).
Test SSLP313 : 2D C_PLAN (fissure maillée)#
Il s’agit d’une fissure inclinée dans un milieu infini soumise à un champ de contraintes uniforme dans une direction (solution de référence analytique en contraintes planes, exacte en milieu infini). La fissure s’ouvre en mode mixte (\(\mathit{K1}\) et \(\mathit{K2}\) ) [V3.02.313].
Pour le test, la fissure est maillée dans une plaque assez grande. Le maillage quadratique est très fin. Les résultats sont les suivants :
Solution de référence (solution analytique)
\(\mathit{K1}\) |
\(\mathit{K2}\) |
\(G\) |
3.58E+06 |
2.69E+06 |
1.00E+02 |
Calcul avec la méthode thêta (CALC_G)
\(\mathit{K1}\) |
\(\mathit{K2}\) |
\(G\) |
|
CALC_G sans nœud au quart |
3.60E+06 |
2.70E+06 |
1.01E+02 |
Ecart / réf. |
0.8 % |
0.2 % |
1.1 % |
CALC_G avec nœuds au quart |
3.60E+06 |
2.70E+06 |
1.01E+02 |
Ecart / réf. |
0.8 % |
0.2 % |
1.2 % |
POST_K1_K2_K3: maillage sans nœud d’arêtes au quart
méthode |
\({\mathit{K1}}_{max}\) |
\({\mathit{K1}}_{min}\) |
\({\mathit{K2}}_{max}\) |
\({\mathit{K2}}_{min}\) |
\({G}_{max}\) |
\({G}_{min}\) |
Écart \({G}_{max}\) / réf. |
Écart \({G}_{min}\) / réf. |
1 |
3.54E+06 |
3.19E+06 |
2.63E+06 |
1.92E+06 |
9.73E+01 |
6.94E+01 |
–3,33% |
–30,70% |
2 |
3.51E+06 |
3.33E+06 |
2.61E+06 |
2.25E+06 |
9.57E+01 |
8.08E+01 |
–4,50% |
–19,32% |
3 |
3.50E+06 |
2.59E+06 |
9.47E+01 |
-5,47% |
||||
POST_K1_K2_K3: maillage avec nœuds d’arêtes au quart
méthode |
\({\mathit{K1}}_{max}\) |
\({\mathit{K1}}_{min}\) |
\({\mathit{K2}}_{max}\) |
\({\mathit{K2}}_{min}\) |
\({G}_{max}\) |
\({G}_{min}\) |
Écart \({G}_{max}\) / réf. |
Écart \({G}_{min}\) / réf. |
1 |
3.61E+06 |
3.60E+06 |
2.70E+06 |
2.69E+06 |
1.01E+02 |
1.01E+02 |
1,29% |
1,07% |
2 |
3.60E+06 |
3.53E+06 |
2.69E+06 |
2.65E+06 |
1.01E+02 |
9.75E+01 |
1,02% |
–2,67% |
3 |
3.56E+06 |
2.66E+06 |
9.88E+01 |
-1,42% |
||||
Sur ce test on constate quele maillage de type «Barsoum» est indispensable si on veut des résultats précis. Avec «Barsoum» la méthode 1 est plus stable. Elle fournit des valeurs de \(G\) (à partir de \(\mathit{K1}\) et \(\mathit{K2}\) ) à environ 1% de la solution analytique. Les méthodes 2 et 3 conduisent à des erreurs de 1 à 2,5 %. On note que dans ce cas, la méthode par extrapolation des déplacements est aussi précise que la méthode thêta.
Par contre, avec un maillage normal, les résultats de la méthode par extrapolation varient beaucoup (entre –3% et -30% de la solution). Il en est de même avec des éléments linéaires. Dans le cas d’un maillage sans éléments de «Barsoum», la méthode 3 est la plus précise.
Test SSLV134 : 3D (fissure maillée)#
Il s’agit d’une fissure plane en forme de disque dans un milieu infini 3D soumise à un champ de contraintes uniforme dans une direction (solution de référence analytique connue sous le nom de «penny shape crack»). La fissure s’ouvre en mode 1 pur, et le \(\mathit{K1}\) est constant le long du fond de fissure [V3.04.134].
Pour ce test, la fissure est maillée dans un bloc parallélépipède. Le maillage est relativement grossier.
Solution de référenceanalytique :
\(\mathit{K1}\) |
\(G\) local |
1,59 106 |
11,59 |
Calcul avec la méthode thêta (CALC_G)
\(G\) |
|
CALC_G avec nœuds au quart |
11.75 |
Ecart / réf. |
1.3 % |
POST_K1_K2_K3: maillage sans nœuds d’arêtes au quart
méthode |
\({\mathit{K1}}_{max}\) |
\({\mathit{K1}}_{min}\) |
\({G}_{max}\) |
\({G}_{min}\) |
Écart \({G}_{max}\) / réf. |
Écart \({G}_{min}\) / réf. |
1 |
1.56E+06 |
1.45E+06 |
1.11E+01 |
9.63E+00 |
–4,32% |
–16,91% |
2 |
1.53E+06 |
1.49E+06 |
1.06E+01 |
1.01E+01 |
–8,35% |
–13,08% |
3 |
1.52E+06 |
1.05E+01 |
–9,51% |
|||
POST_K1_K2_K3: maillage avec nœuds d’arêtes au quart
méthode |
\({\mathit{K1}}_{max}\) |
\({\mathit{K1}}_{min}\) |
\({G}_{max}\) |
\({G}_{min}\) |
Écart \({G}_{max}\) / réf. |
Écart \({G}_{min}\) / réf. |
1 |
1.61E+06 |
1.59E+06 |
1.18E+01 |
1.16E+01 |
1,32% |
–0,06% |
2 |
1.59E+06 |
1.53E+06 |
1.15E+01 |
1.07E+01 |
–0,42% |
–7,87% |
3 |
1.55E+06 |
1.10E+01 |
–5,16% |
|||
Sur ce test on constate encore quele maillage de type «Barsoum» est indispensable si on veut des résultats précis. Avec «Barsoum» la méthode 1 est la plus stable, avec un écart à la solution de référence inférieur à 1,5 % pour \(G\) . Le maillage est relativement grossier, ce qui explique pourquoi la méthode thêta est plus précise.
Test SSLV134 : 3D (fissure non maillée)#
Le cas considéré est le même que celui du paragraphe précédent, mais cette fois-ci la fissure n’est pas maillée. Elle est définie directement dans le fichier de commande, en utilisant la méthode X-FEM [R7.02.12]. Le maillage n’étant pas régulier vis-à-vis de la fissure, les valeurs de \(K\) et de \(G\) calculées varient le long du fond de fissure. Pour la comparaison ci-dessous, on retient la valeur correspondant à un point particulier choisi arbitrairement (milieu du fond de fissure représenté).
Le maillage est linéaire et relativement grossier . Dans la méthode X-FEM, l’utilisateur peut choisir la zone sur laquelle les éléments autour du fond de fissure sont enrichis avec les déplacements asymptotiques (mots clé RAYON_ENRI et NB_COUCHES de DEFI_FISS_XFEM). Cet enrichissement vise à améliorer la précision du calcul. On compare ici les résultats obtenus avec un enrichissement limité aux seuls éléments contenant le fond de fissureet avec un enrichissement sur quatre couches d’éléments autour du fond de fissure.
Calcul avec la méthode thêta (CALC_G – lissage par défaut de type LEGENDREde degré 5)
\(G\) |
|
CALC_G avec enrichissement sur une couche |
11.42 |
Écart / réf. |
-1.4 % |
CALC_G avec enrichissement sur quatre couches |
11.61 |
Écart / réf. |
0.2 % |
POST_K1_K2_K3: enrichissement sur une seule couche
méthode |
\({\mathit{K1}}_{max}\) |
\({\mathit{K1}}_{min}\) |
\({G}_{max}\) |
\({G}_{min}\) |
Écart \({G}_{max}\) / réf. |
Écart \({G}_{min}\) / réf. |
1 |
1.65E+06 |
1.43E+06 |
12.4 |
9.34 |
6,99% |
–19,41% |
2 |
1.52E+06 |
1.44E+06 |
10.5 |
9.45 |
–9,41% |
–18,46% |
3 |
1.47E+06 |
9.81 |
–15,35% |
|||
POST_K1_K2_K3: enrichissement sur quatre couches
méthode |
\({\mathit{K1}}_{max}\) |
\({\mathit{K1}}_{min}\) |
\({G}_{max}\) |
\({G}_{min}\) |
Écart \({G}_{max}\) / réf. |
Écart \({G}_{min}\) / réf. |
1 |
1.58E+06 |
1.58E+06 |
11.3 |
11.3 |
–2,51% |
–2,51% |
2 |
1.55E+06 |
1.47E+06 |
10.9 |
9.88 |
–5,95% |
–14,65% |
3 |
1.51E+06 |
10.4 |
–10,26% |
|||
Sur ce test, on constate qu’il est indispensable d’enrichir sur plusieurs couches d’éléments autour du fond de fissure pour avoir des résultats satisfaisants. À noter que le maillage utilisé ici est linéaire et relativement grossier: avec un maillage plus fin, les résultats sont significativement améliorés. Une étude de convergence sur un cas similaire est présentée dans [bib5].
Avec un enrichissement sur quatre couches, la méthode 1 est celle qui conduit aux résultats les plus précis. L’abscisse curviligne maximale correspond, dans les deux cas, à la distance de quatre éléments environ. La méthode thêta est quant à elle ici moins sensible au paramètre d’enrichissement.
Conclusion#
Les résultats obtenus avec la méthode d’extrapolation du déplacement sont dans l’ensemble satisfaisants, avec moins de 5% d’erreur du \(G\) , surtout si les éléments du fond de fissure sont de type Barsoum (cas fissure maillée) ou si plusieurs couches d’éléments sont enrichis autour du fond (cas fissure non maillée, méthode X-FEM). Dans les deux cas, il s’agit de capter au mieux le comportement asymptotique du déplacement.
Il faut en effet remarquer que l’expression asymptotique des déplacements n’est valable que pour \(r\) tendant vers 0. C’est pourquoi il faut veiller à ne pas choisir un domaine d’extrapolation trop grand (distance \(\mathit{dmax}\) de l’opérateur POST_K1_K2_K3 de l’ordre de 4 à 5 éléments).
Sur les tests présentés pour une fissure maillée, la méthode 1 donne les résultats les plus précis et les plus stables, que ce soit en 2D ou en 3D, s’il y a des éléments de Barsoum. Si le maillage ne comporte pas d’éléments de Barsoum, on conseille alors d’utiliser les résultats de la méthode 3. Pour une fissure non maillée, la méthode 1 semble aussi la plus précise.
Sur une étude pour laquelle on ne connaît pas de solution de référence, il est possible d’estimer a posteriori la qualité du calcul. En effet, POST_K1_K2_K3 fournit systématiquement pour les deux premières méthodes les valeurs maximum et les valeurs minimum (sur l’ensemble des points calculés) des facteurs d’intensité des contraintes, ainsi que la valeur de G recalculée par la formule d’Irwin. La méthode 3 ne fournit quant à elle qu’une seule valeur pour chaque facteur d’intensité de contrainte. Cette méthode est une moyenne pondérée des facteurs d’intensité des contraintes extrapolés en chaque nœud.
Un résultat peut être considéré comme satisfaisant si les 5 valeurs ainsi fournies (min et max des méthodes 1 et 2, et méthode 3) sont proches. On peut également recommander de comparer les résultats obtenus avec cette méthode avec ceux issus du calcul du taux de restitution de l’énergie et des fonctions singulières (opérateur CALC_G).
Bibliographie#
H.D.BUI: «Mécanique de la rupture fragile» - Masson, 1978.
LEMAITRE, J.L.CHABOCHE: «Mécanique des matériaux solides» - Dunod, 1996.
J.B. LEBLOND: « Mécanique de la rupture fragile et ductile »– Lavoisier, 2003.
R.S. BARSOUM: « Triangular quarter-point elements as elastic and perfectly-plastic crack tip elements »– Int. J. for Numerical Methods in Engineerig, vol. 11, 85-98, 1977.
GENIAUT : “Convergence en mécanique de la rupture : validation des éléments finis classiques et enrichis dans Code_Aster” – Note EDF R&D H-T64-2008-0047, 2008.
Description des versions du document#
Indice |
Version Aster |
Auteur(s) ou contributeur(s), organisme |
Description des modifications |
A |
5 |
J.M.Proix EDF/R&D/AMA |
Texte initial |
B |
7.4 |
E.Galenne EDF/R&D/AMA |
Méthode des moindres carrés 7.2.24 |
C |
9.4 |
E.Galenne EDF/R&D/AMA |
Fissure non plane avec XFEM 9.2.8 |