v6.03.103 SSNP103 - Calcul du taux de restitution d’énergie en élasticité non linéaire#

Résumé:

Ce test est un cas-test en mécanique de la rupture non linéaire.

On considère une plaque rectangulaire finie, fissurée et soumise à un chargement de traction. La loi de comportement utilisée est une loi élastique non linéairede type Henky von Mises.

Ce cas-test comprend quatre modélisations en \(\mathrm{2D}\) :

  • La modélisation A, en déformations planes, traite d’unefissure sollicitéeen mode 1 pur. Seule la moitié de la plaque est modélisée par symétrie.

  • La modélisation B, en déformations planes, traite d’une fissure sollicitée en modes mixtes. La plaque complète est modélisée.

  • La modélisation C reprend la géométrie de la modélisation A mais résout le problème axisymétrique associé. Il ne s’agit donc pas ici d’une plaque fissurée mais d’un cylindre fissuré.

  • La modélisation D, reproduit le calcul de la modélisation B en y ajoutant un chargement thermique.

  • La modélisation E, reproduit le calcul de la modélisation B avec HHO linéaire.

  • La modélisation F, reproduit le calcul de la modélisation B avec HHO quadratique.

Solution de référence#

Dans ce test nous fabriquons la solution de référence en faisant une estimation du taux de restitution d’énergie de deux manières.

La première solution est obtenue par laméthode de G-thêta classique (opérateur CALC_G), c’est à dire, en faisant une estimation de l’intégrale surfacique sur une couronne placée autour de la fissure.

Afin de trouver cette même valeur par un autre moyen, nous remarquons d’abord que dans le cas de l’élasticité linéaire, il est possible d’estimer \(G\) grâce à la formule d’Irwin faisant le lien entre les quantités globalesavec les quantités locales que sont les facteurs d’intensité de contrainte \({K}_{I};{K}_{\mathit{II}}\) :

(4731)#\[G=\frac{1-{\nu}^{2}}{E}({K}_{I}^{2}+{K}_{\mathit{II}}^{2})\]

Et:

(4732)#\[G=\frac{({K}_{I}^{2}+{K}_{\mathit{II}}^{2})}{E}\]

Les facteurs d’intensité de contraintes peuvent être estimés par l’analyse de singularité des champs de contrainte ou de déplacement proche de la pointe. Dans le cas del’hyper-élasticité, le taux de restitution d’énergie est encore un invariant du contour d’intégration, ce qui nous permet d’imaginer une procédure similaire.

D’abord nous pouvons théoriquement contracter indéfiniment le contour d’intégration vers la pointe de fissure tout en gardant la même valeur de \(G\) . Donc même pour une loi hyper-élastique nous pouvons établir un lien entre les quantités globale \(G\) et locale \(K\) . Cependant la généralisation de la notion de facteur d’intensité de contraintes n’est pas toujours aisée pour tous les types des lois hyper-élastiques. Par contre pourune loi hyper-élastique bilinéaire ce lien peut être établi sans trop de difficulté. Pour cette loi, au voisinage de la pointe, le comportement sera linéaire et ceci peu importe le chargement extérieur appliqué. Nous retrouverons doncla singularité de contrainte et de déformation classique sous forme de \(\sqrt{r}\) .

Dans ce cas-test nous utilisons donc la loi élastique non-linéaire de type Hencky, qui a l’avantage de reproduire la loi de plasticité de von Mises pour tout chargement radial monotone:

(4733)#\[\underline{\underline{\sigma}}(\underline{\underline{\epsilon}})=K\mathit{Tr}(\underline{\underline{\epsilon}})\underline{\underline{I}}+2\mu ({\epsilon}^{\mathit{eq}})\underline{\underline{{\epsilon}^{D}}}\]

où \(3K=E/(1-2\nu )\) est le module de compressibilité (en anglais:bulk modulus) et \(\mu ({\epsilon}^{\mathit{eq}})\) est un coefficient de Lamé, qui est supposé dépendant exclusivement de la déformation équivalente de von Mises: \({\epsilon}^{\mathit{eq}}=\sqrt{\frac{3}{2}\cdot {\underline{\underline{\epsilon}}}^{D}:{\underline{\underline{\epsilon}}}^{D}}\) .

Dans laforme générale de la loi non-linéaire de Hencky, la dépendance du coefficient de Lamé \(\mu ({\epsilon}^{\mathit{eq}})\) en fonction de la déformation est identifiée à partir de la courbe de traction uniaxiale. Pour le cas bilinéaire la loi s’écrit ainsi dans le domaine d’élasticité linéaire initial:

(4734)#\[\underline{\underline{\sigma}}(\underline{\underline{\epsilon}})=K\mathit{Tr}(\underline{\underline{\epsilon}})\underline{\underline{I}}+2\mu \underline{\underline{{\epsilon}^{D}}}\]

Et dans le domaine d’élasticité non-linéaire:

(4735)#\[\underline{\underline{\sigma}}(\underline{\underline{\epsilon}})=K\mathit{Tr}(\underline{\underline{\epsilon}})\underline{\underline{I}}+{\sigma}_{c}(1-{\mu}_{\infty}/\mu )\underline{\underline{{\epsilon}^{D}}}/\sqrt{{\underline{\underline{\epsilon}}}^{D}:{\underline{\underline{\epsilon}}}^{D}}+2{\mu}_{\infty}\underline{\underline{{\epsilon}^{D}}}\]

avec \(3K=E/(1-2\nu )\) et \(2\mu =E/(1+\nu )\) .

Le coefficient de Lamé \(2{\mu}_{\infty}\) pour les chargements dépassant le seuil d’élasticité linéaire est lié à la deuxième pente de la loi, qu’on note \({E}_{\infty}\) (mot-clé D_SIGM_EPSI):

(4736)#\[2{\mu}_{\infty}=6{E}_{\infty}K/(9K+{E}_{\infty})\]

Enfin, pour les chargements importants \({\epsilon}^{\mathit{eq}}\gg \frac{{\sigma}_{c}}{2\mu }\) , cette loi devient asymptotiquement linéaire avec les caractéristiques linéaires suivantes:

(4737)#\[\underline{\underline{\sigma}}(\underline{\underline{\epsilon}})=\frac{{E}_{\infty}}{3(1-2{\nu}_{\infty})}\mathit{Tr}(\underline{\underline{\epsilon}})\underline{\underline{I}}+\frac{{E}_{\infty}}{1+{\nu}_{\infty}}\underline{\underline{{\epsilon}^{D}}}+o(\underline{\underline{\epsilon}})\]

Avec \(E\to {E}_{\infty}\) et \(2{\nu}_{\infty}=1-(1-2\nu )\cdot {E}_{\infty}/E\) .

L’analyse de la singularité du saut de déplacement avec ces nouvellesvaleurs de coefficients d’élasticité à chargements «infinis» nous permet d’identifier les valeurs de facteurs d’intensité de contrainte grâce à l’opérateur déjà disponible dans code_aster: POST_K1_K2_K3, puis utiliser la formule d’Irwin classique:

(4738)#\[G=\frac{1-{\nu}_{\infty}^{2}}{{E}_{\infty}}({K}_{I}^{2}+{K}_{\mathit{II}}^{2})\]
../../../../_images/10000200000002BD00000144FB1373D443B21CE3.png

Ainsi, deux tests sont effectués sur les taux de restitution d’énergie calculés.

Premier test:

On teste directement le ratio entre le G issu de CALC_G et celui issu de POST_K1_K2_K3. Ce ratio doit être égal à 1.

Il s’agit ici d’une référence de typeANALYTIQUE.

Second test:

Le cas-test a été exécuté, pour chaque modélisation, sur un maillage dix fois plus fin en fond de fissure puis sur un maillage cent fois plus fin (taille des mailles connectées au fond de fissure respectivement de \({10}^{-4}\mathit{mm}\) et \({10}^{-5}\mathit{mm}\) ), avec un chargement couvrant une plage plus importante que celle du cas test.

A chaque instant de calcul, les valeurs de \(G\) issues de la méthode 1 et 3 de calcul de POST_K1_K2_K3 (voir documentation U4.82.05) ont été extraites.

La comparaison des grandeurs calculées pour chaque maillage montre la convergence (au sens du maillage) de celles-ci.

La solution de référence retenue est la valeur de \(G\) obtenue avec la méthode 3 avec le maillage le plus fin (taille de mailles connectées au fond de fissure de \({10}^{-5}\mathit{mm}\) ). L’écart relatif entre la valeur de \(G\) issue des méthodes 1 et 3 donne une estimation de la précision de cette valeur de référence.

Il s’agit ici d’une référence de type AUTRE_ASTER.

Modélisation A#

Caractéristiques de la modélisation#

Il s’agit d’un calcul en élasticité non linéaire sous l’hypothèse des petits déplacements, en déformations planes. La fissure est ouverte en mode I pur. Seule la moitié de la structure est représentée par symétrie.

Caractéristiques du maillage#

Le maillage en entrée du cas-test en un maillage linéaire ayant les caractéristiques suivantes:

-Nombre de nœuds: 8936

-Nombre d’éléments:

-Triangles: 6494

-quadrangles: 5408

Le rayon du domaine rayonnant en fond de fissure est de \(25\mathit{mm}\) . La taille de éléments connectés au fond de fissure est de \({10}^{-3}\mathit{mm}\)

A l’exécution du cas-test, un nouveau maillage quadratique doté d’éléments de BARSOUM en fond de fissure est créé par l’exécution des commandes CREA_MAILLAGE (mot clé LINE_QUAD) puis MODI_MAILLAGE (mot-clé MODI_MAILLE, option NOEUD_QUART).

../../../../_images/10000201000001700000026D166DD6E2353657AC.png

Grandeurs testées et résultats#

On compare les résultats obtenus pour trois couronnes d’intégration différentes:

  • couronne 1: \({R}_{inf}=10\mathit{mm}\)

; \({R}_{\sup}=15\mathit{mm}\)

  • couronne 2: \({R}_{inf}=15\mathit{mm}\)

; \({R}_{\sup}=20\mathit{mm}\)

  • couronne 3: \({R}_{inf}=20\mathit{mm}\)

; \({R}_{\sup}=25\mathit{mm}\)

Les valeurs de référence de G issues de la méthode 3 POST_K1_K2_K3 retenues ainsi que les précisions associées sont fournies dans le tableau ci-dessous.

Instant de calcul

Effort surfacique imposé (MPa)

Référence : \(G(N/m)\) Issu de POST_K1_K2_K3(méthode 3)

Tolérance (%)

1

\(0,7\ast \text{SY}=105\)

1.13702E+05

0,3

2

\(\text{SY}=150\)

2.85869E+05

0,3

3

\(3\ast \text{SY}=450\)

3.42359E+06

0,2

Modélisation B#

Caractéristiques de la modélisation#

Il s’agit d’un calcul en élasticité non linéaire sous l’hypothèse des petits déplacements, en déformations planes. La fissure est ouverte en mode mixte. L’intégralité de la structure est représentée.

Caractéristiques du maillage#

Le maillage en entrée du cas-test en un maillage linéaire ayant les caractéristiques suivantes:

-Nombre de nœuds: 17691

-Nombre d’éléments:

-Triangles: 12988

-quadrangles: 10816

Le rayon du domaine rayonnant en fond de fissure est de \(25\mathit{mm}\) . La taille de éléments connectés au fond de fissure est de \({10}^{-3}\mathit{mm}\)

Il est obtenu par «mirrorisation» du maillage de la modélisation A.

A l’exécution du cas-test, un nouveau maillage quadratique doté d’éléments de BARSOUM en fond de fissure est créé par l’exécution des commandes CREA_MAILLAGE (mot clé LINE_QUAD) puis MODI_MAILLAGE (mot-clé MODI_MAILLE, option NOEUD_QUART).

Grandeurs testées et résultats#

On compare les résultats obtenus pour trois couronnes d’intégration différentes:

  • couronne 1: \({R}_{inf}=10\mathit{mm}\)

; \({R}_{\sup}=15\mathit{mm}\)

  • couronne 2: \({R}_{inf}=15\mathit{mm}\)

; \({R}_{\sup}=20\mathit{mm}\)

  • couronne 3: \({R}_{inf}=20\mathit{mm}\)

; \({R}_{\sup}=25\mathit{mm}\)

Les valeurs de référence de G issues de la méthode 3 POST_K1_K2_K3 retenues ainsi que les précisions associées sont fournies dans le tableau ci-dessous.

Instant de calcul

Deplacement imposé :math:`delta ` (en mm)

Référence : \(G(N/m)\) Issu de POST_K1_K2_K3(méthode 3)

Tolérance (%)

1

\(0,7\ast \text{Dp}=200\)

1.36596E+04

0, 5

2

\(\text{Dp}=286\)

2.74150E+04

0, 4

3

\(3\ast \text{Dp}=857\)

1.68193E+05

0, 3

Où \(\mathit{Dp}=L.\mathit{SY}/E\) .

Modélisation C#

Caractéristiques de la modélisation#

Il s’agit d’un calcul en élasticité non linéaire sous l’hypothèse des petits déplacements, en axisymétrique.La fissure est ouverte en mode I pur. Seule la moitié de la structure est représentée par symétrie.

Caractéristiques du maillage#

Le maillage est identique à celui de la modélisation A.

Grandeurs testées et résultats#

On compare les résultats obtenus pour trois couronnes d’intégration différentes:

  • couronne 1: \({R}_{inf}=10\mathit{mm}\)

; \({R}_{\sup}=15\mathit{mm}\)

  • couronne 2: \({R}_{inf}=15\mathit{mm}\)

; \({R}_{\sup}=20\mathit{mm}\)

  • couronne 3: \({R}_{inf}=20\mathit{mm}\)

; \({R}_{\sup}=25\mathit{mm}\)

Les valeurs de référence de G issues de la méthode 3 POST_K1_K2_K3 retenues ainsi que les précisions associées sont fournies dans le tableau ci-dessous.

Instant de calcul

Effort surfacique imposé (MPa)

Référence : \(G(N/m)\) Issu de POST_K1_K2_K3(méthode 3)

Tolérance (%)

1

\(0,7\ast \text{SY}=105\)

3.35679E+03

1,5

2

\(\text{SY}=150\)

1.01970E+04

0, 6

3

\(3\ast \text{SY}=450\)

1.64979E+05

0, 3

Modélisation D#

Caractéristiques de la modélisation#

Il s’agit d’un calcul en élasticité non linéaire sous l’hypothèse des petits déplacements, en déformations planes. La fissure est ouverte en mode mixte. L’intégralité de la structure est représentée. Une température stationnaire non uniforme est appliquée.

Caractéristiques du maillage#

Le maillage est identique à celui de la modélisation B.

Grandeurs testées et résultats#

On compare les résultats obtenus pour trois couronnes d’intégration différentes:

  • couronne 1: \({R}_{inf}=10\mathit{mm}\)

; \({R}_{\sup}=15\mathit{mm}\)

  • couronne 2: \({R}_{inf}=15\mathit{mm}\)

; \({R}_{\sup}=20\mathit{mm}\)

  • couronne 3: \({R}_{inf}=20\mathit{mm}\)

; \({R}_{\sup}=25\mathit{mm}\)

Les valeurs de référence de G issues de la méthode 3 POST_K1_K2_K3 retenues ainsi que les précisions associées sont fournies dans le tableau ci-dessous.

Instant de calcul

Deplacement imposé :math:`delta ` (en mm)

Référence : \(G(N/m)\) Issu de POST_K1_K2_K3(méthode 3)

Tolérance (%)

1

\(0,7\ast \text{Dp}=200\)

5.00690E+04

2

2

\(\text{Dp}=286\)

7.00295E+04

2

3

\(3\ast \text{Dp}=857\)

2.44290E+05

1

Où \(\mathit{Dp}=L.\mathit{SY}/E\) .

Modélisation E#

Caractéristiques de la modélisation#

Il s’agit d’un calcul en élasticité non linéaire sous l’hypothèse des petits déplacements, en déformations planes. La fissure est ouverte en mode mixte. L’intégralité de la structure est représentée. On utilise la modélisation HHO linéaire

Caractéristiques du maillage#

Le maillage est identique à celui de la modélisation B.

Grandeurs testées et résultats#

Les valeurs de référence de G, K1 et K2 issues de la méthode POST_K1_K2_K3 de la modélisation B.

Modélisation F#

Caractéristiques de la modélisation#

Il s’agit d’un calcul en élasticité non linéaire sous l’hypothèse des petits déplacements, en déformations planes. La fissure est ouverte en mode mixte. L’intégralité de la structure est représentée. On utilise la modélisation HHO quadratique.

Caractéristiques du maillage#

Le maillage est identique à celui de la modélisation B.

Grandeurs testées et résultats#

Les valeurs de référence de G, K1 et K2 issues de la méthode POST_K1_K2_K3 de la modélisation B.

Synthèse des résultats#

Ce cas test vise à valider le calcul du taux de restitution d’énergie en élasticité non linéaire:

      • les modélisationsA et B permettent de valider le calcul de G en déformations planes, en mode 1 pur et en mode mixte

      • la modélisation C valide le calcul de G en axisymétrique, en mode 1 pur.

      • la modélisation D valide le calcul de G suite à une calcul thermo-mécanique.