v7.31.143 WTNV143 - Application d’une pression répartie sur les lèvres d’une fissure XFEM dans un modèle hydro-mécanique#
Résumé:
Le but de ce test est de s’assurer du bon fonctionnement de la méthode des éléments finis étendue (XFEM) dans un modèle hydro-mécanique, afin de représenter un milieu saturé fracturé.
Dans ce test, on suppose que la pression de pore est prise uniformément nulle dans tout le domaine d’étude: le modèle hydro-mécanique correspond alorsà un modèle mécanique «classique» sans couplage. On vérifie qu’il existe bien une discontinuité du champ de déplacement de part et d’autre de l’interface. Les résultats sont comparés à une solution analytique et aux résultats obtenus avec un test XFEM mécanique (semblables aux tests ssnv203, modélisations F et H).
Solution de référence#
Méthode de calcul#
Il s’agit d’une solution analytique. Compte tenu des conditions aux limites, les déplacements peuvent être obtenus à partir de la résolution analytique de l’équation de conservation de la quantité de mouvement.
En négligeant la pesanteur, l’équation s’écrit (en contraintes totales):
\(\text{Div}(\sigma )+{r}_{0}g=0\)
Dans le cas d’une modélisation couplée, le tenseur des contraintes totales s’écrit:
\(\sigma =\sigma '-{p}_{1}1\)
\(\sigma '\) est le tenseur des contraintes dans le squelette et \({p}_{1}\) et la pression de pore dans le massif. Le module de Poisson \(\nu\) étant nul , et étant dans le cas élastique, on a \(\sigma '=Eϵ\) .
Or \(\forall y{p}_{1}=0\) donc finalement \(E\ast \text{Div}({\varepsilon})+{r}_{0}g=0\)
Dans le cas 2D, \(\nu\) étant nul, les conditions aux limites et le chargement rendent le problème unidimensionnel selon \(y\) . Seul \({ϵ}_{yy}\) est non nul et:
\(E\ast \frac{\partial {{\varepsilon}}_{yy}}{\partial y}-{r}_{0}g=0\) soit \(\frac{{\partial}^{2}{u}_{y}}{\partial {y}^{2}}=\frac{{r}_{0}}{E}g\) et donc finalement \({u}_{y}(y)=\frac{{r}_{0}g}{2E}{y}^{2}+\mathit{Ay}+B\)
avec \(A\) et \(B\) des constantes d’intégration à déterminer. On résout le problème séparément dans chacun des deux sous-blocs. Afin de déterminer les constantes d’intégration, on utilise la condition aux limites sur les déplacements \({u}_{y}=0\) aux extrémités inférieures et supérieures de la colonne et la condition de Neumann \({\sigma}_{yy}=-p\) au niveau de l’interface.
Finalement, on obtient le déplacement selon la direction \(y\) de part et d’autre de l’interface:
les déplacements au dessus de l’interface en \(y=\frac{{\mathit{LY}}^{+}}{2}\) s’écrivent \({u}_{y}\left(\frac{{\mathit{LY}}^{+}}{2}\right)=\frac{p}{E}\left(\frac{\mathit{LY}}{2}\right)-\frac{{r}_{0}g}{2E}{\left(\frac{\mathit{LY}}{2}\right)}^{2}\)
les déplacements en dessous de l’interface en \(y=\frac{{\mathit{LY}}^{-}}{2}\) s’écrivent \({u}_{y}\left(\frac{{\mathit{LY}}^{-}}{2}\right)=-\frac{p}{E}\left(\frac{\mathit{LY}}{2}\right)-\frac{{r}_{0}g}{2E}{\left(\frac{\mathit{LY}}{2}\right)}^{2}\)
Dans le cas 3D, la solution de référence est rigoureusement la même:
les déplacements au dessus de l’interface en \(z=\frac{{\mathit{LZ}}^{+}}{2}\) s’écrivent \({u}_{z}\left(\frac{{\mathit{LZ}}^{+}}{2}\right)=\frac{p}{E}\left(\frac{\mathit{LZ}}{2}\right)-\frac{{r}_{0}g}{2E}{\left(\frac{\mathit{LZ}}{2}\right)}^{2}\)
les déplacements en dessous de l’interface en \(z=\frac{{\mathit{LZ}}^{-}}{2}\) s’écrivent \({u}_{z}\left(\frac{{\mathit{LZ}}^{-}}{2}\right)=-\frac{p}{E}\left(\frac{\mathit{LZ}}{2}\right)-\frac{{r}_{0}g}{2E}{\left(\frac{\mathit{LZ}}{2}\right)}^{2}\)
les déplacements suivant \(x\) et \(y\) sont nuls.
Grandeurs et résultats de référence#
On teste la valeur des déplacements sur les lèvres inférieure et supérieure de l’interface:
Déplacements verticaux |
Déplacements horizontaux |
|
Lèvre inférieure |
-4.323558729E-03 |
0 |
Lèvre supérieure |
4.297130927-03 |
0 |
Incertitudes sur la solution#
Aucune la solution est analytique.
Références bibliographiques#
Documentation de référence R7.02.12 (eXtended Finite Element Method).
Modélisation A#
Caractéristiques de la modélisation#
Il s’agit d’une modélisation D_PLAN_HM utilisant des éléments HM-XFEM quadratiques. Le barreau sur lequel on effectue la modélisation est divisé en 5 QUAD8. L’interface est non maillée et coupe l’élément central. Ainsi on a 3 éléments HM-XFEM et 2 éléments HM classiques. Comme indiqué sur la Figure , les 3 éléments XFEM subissent un sous découpage en sous triangles (pour effectuer l’intégration de Gauss-Legendre de part et d’autre des lèvres de l’interface, mais ces sous-éléments triangulaires ne sont pas des éléments du maillage).
Figure 3.1-a : Caractéristiques de la modélisation
Caractéristiques du maillage#
Le maillage est constitué de 5 mailles quadrangles quadratiques (QUAD8).
Grandeurs testées et résultats#
Les résultats (résolution avec STAT_NON_LINE) sont synthétisés dans le tableau ci-dessous pour la direction \(y\) . Pour tester tous les nœuds du barreau en même temps, on calcule le MIN et le MAX.
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance (%) |
DY(en dessous) MIN |
“ANALYTIQUE” |
-4.323558729E-03 |
0,0001 |
DY(en dessous) MAX |
“ANALYTIQUE” |
-4.323558729E-03 |
0,0001 |
DY(au dessus) MIN |
“ANALYTIQUE” |
4.297130927-03 |
0,0001 |
DY(au dessus) MAX |
“ANALYTIQUE” |
4.297130927-03 |
0,0001 |
Les résultats (résolution avec STAT_NON_LINE) sont synthétisés dans le tableau ci-dessous pour la direction \(x\) . Pour tester tous les nœuds du barreau en même temps, on calcule le MIN et le MAX.
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance (%) |
DX(en dessous) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DX(en dessous) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
DX (au dessus) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DX(au dessus) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
On peut alors constater (à partir de la Figure ) une discontinuité franche du champ de déplacements liée à la présence de l’interface traversant le massif. Cela suggère la bonne prise en compte de l’enrichissement dans l’approximation du champ de déplacements par la fonction Heaviside.
Figure 3.3-a : Champ de déplacements selon la direction (Oy)
On effectue ensuite une seconde modélisation qui utilise les mêmes paramètres que la précédente, mais cette fois ci en effectuant le calcul XFEM (qui est alors un calcul de mécanique pure) avec l’opérateur MECA_STATIQUE au lieu de l’opérateur STAT_NON_LINE dans le cas HM-XFEM. Les résultats obtenus sont strictement identiques aux précédents.
Modélisation B#
Caractéristiques de la modélisation#
Cette modélisation est strictement identique à la précédente. Afin de tester que le modèle HM-XFEM est fonctionnel avec des TRIA6, nous choisissons de subdiviser les 5 QUAD8 de la modélisation précédente, chacun en 2 TRIA6. Ce qui fait au total 10 éléments triangulaires quadratiques. Comme précédemment les éléments triangulaires contenus dans les quadrangles extrêmes sont non enrichis. A contrario ceux contenus dans les 3 quadrangles centraux sont enrichis par la fonction Heaviside.
Caractéristiques du maillage#
Le maillage est constitué de 10 mailles triangulaires quadratiques (TRIA6).
Grandeurs testées et résultats#
Les résultats (résolution avec STAT_NON_LINE) sont synthétisés dans le tableau ci-dessous selon la direction \(y\) . La tolérance est fixée à \({10}^{\text{-6}}\) . Pour tester tous les nœuds du barreau en même temps, on calcule le MIN et le MAX.
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance (%) |
DY(en dessous) MIN |
“ANALYTIQUE” |
-4.323558729E-03 |
0,0001 |
DY(en dessous) MAX |
“ANALYTIQUE” |
-4.323558729E-03 |
0,0001 |
DY(au dessus) MIN |
“ANALYTIQUE” |
4.297130927-03 |
0,0001 |
DY(au dessus) MAX |
“ANALYTIQUE” |
4.297130927-03 |
0,0001 |
Les résultats (résolution avec STAT_NON_LINE) sont synthétisés dans le tableau ci-dessous selon la direction \(x\) . Pour tester tous les nœuds du barreau en même temps, on calcule le MIN et le MAX.
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance (%) |
DX(en dessous) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DX(en dessous) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
DX(au dessus) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DX(au dessus) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
On peut alors constater, comme dans le cas de la modélisation A, une discontinuité franche du champ de déplacement liée à la présence de la l’interface traversant le massif. Cela suggère la bonne prise en compte de l’enrichissement de l’approximation du champ de déplacements par la fonction Heaviside.
Figure 4.3-a : Champ de déplacements selon la direction (Oy)
On effectue ensuite une seconde modélisation qui utilise les mêmes paramètres que la précédente, mais cette fois ci en conduisant le calcul XFEM (qui est alors un calcul de mécanique pur) avec l’opérateur MECA_STATIQUE au lieu de l’opérateur STAT_NON_LINE dans le cas HM-XFEM. Les résultats obtenus sont strictement identiques aux précédents.
Modélisation C#
Caractéristiques de la modélisation#
Il s’agit d’une modélisation 3D_HM utilisant des éléments HM-XFEM quadratiques. La colonne sur laquelle on effectue la modélisation est divisée en 5 HEXA20. L’interface est non maillée et coupe l’élément central. Ainsi on a 3 éléments HM-XFEM et 2 éléments HM classiques (les deux hexaèdres qui forment les extrémités de la colonne). Comme indiqué sur la Figure , les 3 éléments XFEM subissent un sous découpage en sous tétraèdres (pour effectuer l’intégration de Gauss-Legendre de part et d’autre des lèvres de l’interface, mais ces sous-éléments tétraèdres ne sont pas des éléments du maillage).
Figure 5.1-a : Caractéristiques de la modélisation
Caractéristiques du maillage#
Le maillage est constitué de 5 mailles hexaédriques quadratiques (HEXA20).
Grandeurs testées et résultats#
Les résultats (résolution avec STAT_NON_LINE) sont synthétisés dans les tableaux ci-dessous. Pour tester tous les nœuds de la colonne en même temps, on calcule le MIN et le MAX.
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance (%) |
DZ(en dessous) MIN |
“ANALYTIQUE” |
-4.323558729E-03 |
0,0001 |
DZ(en dessous) MAX |
“ANALYTIQUE” |
-4.323558729E-03 |
0,0001 |
DZ(au dessus) MIN |
“ANALYTIQUE” |
4.297130927-03 |
0,0001 |
DZ(au dessus) MAX |
“ANALYTIQUE” |
4.297130927-03 |
0,0001 |
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance (%) |
DX(en dessous) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DX(en dessous) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
DX (au dessus) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DX(au dessus) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance (%) |
DY(en dessous) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DY(en dessous) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
DY(au dessus) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DY(au dessus) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
Figure 5.3-a : Champ de déplacements selon la direction (Oz)
On effectue ensuite une seconde modélisation qui utilise les mêmes paramètres que la précédente, mais cette fois ci en conduisant le calcul XFEM (qui est alors un calcul de mécanique pur) avec l’opérateur MECA_STATIQUE au lieu de l’opérateur STAT_NON_LINE dans le cas HM-XFEM. Les résultats obtenus sont strictement identiques aux précédents.
Modélisation D#
Caractéristiques de la modélisation#
Cette modélisation est strictement identique à la précédente. Afin de tester que le modèle HM-XFEM est fonctionnel avec des PENTA15, nous choisissons de subdiviser les 5 HEXA20 de la modélisation précédente, chacun en 2 PENTA15. Ce qui fait au total 10 éléments pentaèdres quadratiques. Comme précédemment les éléments pentaèdres contenus dans les hexaèdres extrêmes sont non enrichis. A contrario ceux contenus dans les 3 hexaèdres centraux sont enrichis par la fonction Heaviside.
Caractéristiques du maillage#
Le maillage est constitué de 10 mailles pentaèdriques quadratiques (PENTA15).
Grandeurs testées et résultats#
Les résultats (résolution avec STAT_NON_LINE) sont synthétisés dans les tableaux ci-dessous. Pour tester tous les nœuds de la colonne en même temps, on calcule le MIN et le MAX.
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance (%) |
DZ(en dessous) MIN |
“ANALYTIQUE” |
-4.323558729E-03 |
0,0001 |
DZ(en dessous) MAX |
“ANALYTIQUE” |
-4.323558729E-03 |
0,0001 |
DZ(au dessus) MIN |
“ANALYTIQUE” |
4.297130927-03 |
0,0001 |
DZ(au dessus) MAX |
“ANALYTIQUE” |
4.297130927-03 |
0,0001 |
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance (%) |
DX(en dessous) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DX(en dessous) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
DX (au dessus) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DX(au dessus) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance (%) |
DY(en dessous) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DY(en dessous) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
DY(au dessus) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DY(au dessus) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
Figure 6.3-a : Champ de déplacements selon la direction (Oz)
On effectue ensuite une seconde modélisation qui utilise les mêmes paramètres que la précédente, mais cette fois ci en conduisant le calcul XFEM (qui est alors un calcul de mécanique pur) avec l’opérateur MECA_STATIQUE au lieu de l’opérateur STAT_NON_LINE dans le cas HM-XFEM. Les résultats obtenus sont strictement identiques aux précédents.
Modélisation E#
Caractéristiques de la modélisation#
Cette modélisation est strictement identique aux précédentes. Afin de tester que le modèle HM-XFEM est fonctionnel avec des TETRA10, nous choisissons de subdiviser les 5 HEXA20 de la modélisation C, chacun en 6 TETRA10. Ce qui fait au total 30 éléments tétraèdres quadratiques. Comme précédemment les éléments tétraèdres contenus dans les hexaèdres extrêmes sont non enrichis. A contrario ceux contenus dans les 3 hexaèdres centraux sont enrichis par la fonction Heaviside.
Caractéristiques du maillage#
Le maillage est constitué de 30 mailles tétraèdriques quadratiques (TETRA10).
Grandeurs testées et résultats#
Les résultats (résolution avec STAT_NON_LINE) sont synthétisés dans les tableaux ci-dessous. Pour tester tous les nœuds de la colonne en même temps, on calcule le MIN et le MAX.
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance (%) |
DZ(en dessous) MIN |
“ANALYTIQUE” |
-4.323558729E-03 |
0,0001 |
DZ(en dessous) MAX |
“ANALYTIQUE” |
-4.323558729E-03 |
0,0001 |
DZ(au dessus) MIN |
“ANALYTIQUE” |
4.297130927-03 |
0,0001 |
DZ(au dessus) MAX |
“ANALYTIQUE” |
4.297130927-03 |
0,0001 |
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance (%) |
DX(en dessous) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DX(en dessous) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
DX (au dessus) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DX(au dessus) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance (%) |
DY(en dessous) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DY(en dessous) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
DY(au dessus) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DY(au dessus) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
Figure 7.3-a : Champ de déplacements selon la direction (Oz)
On effectue ensuite une seconde modélisation qui utilise les mêmes paramètres que la précédente, mais cette fois ci en conduisant le calcul XFEM (qui est alors un calcul de mécanique pur) avec l’opérateur MECA_STATIQUE au lieu de l’opérateur STAT_NON_LINE dans le cas HM-XFEM. Les résultats obtenus sont strictement identiques aux précédents.
Modélisation F#
Caractéristiques de la modélisation#
Cette modélisation est strictement identique aux précédentes. Afin de tester que le modèle HM-XFEM est fonctionnel avec des PYRAM13, nous choisissons de subdiviser les 5 HEXA20 de la modélisation C, chacun en 6 PYRAM13. Ce qui fait au total 30 éléments pyramidaux quadratiques. Comme précédemment les éléments pyramidaux contenus dans les hexaèdres extrêmes sont non enrichis. A contrario ceux contenus dans les 3 hexaèdres centraux sont enrichis par la fonction Heaviside.
Caractéristiques du maillage#
Le maillage est constitué de 30 mailles pyramidales quadratiques (PYRAM13).
Grandeurs testées et résultats#
Les résultats obtenus avec Code_Aster (résolution avec STAT_NON_LINE) sont synthétisés dans les tableaux ci-dessous. Pour tester tous les nœuds de la colonne en même temps, on calcule le MIN et le MAX.
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance (%) |
DZ(en dessous) MIN |
“ANALYTIQUE” |
-4.323558729E-03 |
0,1 |
DZ(en dessous) MAX |
“ANALYTIQUE” |
-4.323558729E-03 |
0,1 |
DZ(au dessus) MIN |
“ANALYTIQUE” |
4.297130927-03 |
0,1 |
DZ(au dessus) MAX |
“ANALYTIQUE” |
4.297130927-03 |
0,1 |
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance (%) |
DX(en dessous) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DX(en dessous) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
DX (au dessus) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DX(au dessus) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance (%) |
DY(en dessous) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DY(en dessous) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
DY(au dessus) MIN |
“ANALYTIQUE” |
0 |
0,0001 |
DY(au dessus) MAX |
“ANALYTIQUE” |
0 |
0,0001 |
Remarque:
La tolérance pour le test sur DZ est plus élevée, on passe de \({10}^{\text{-6}}\) à \({10}^{\text{-3}}\) . En effet, les fonctions de forme des pyramides ne sont pas des polynômes mais des fractions rationnelles. Cela explique la moins bonne précision obtenue pour ce problème linéaire par rapport aux précédentes modélisations.
Figure 8.3-a : Champ de déplacements selon la direction (Oz)
On effectue ensuite une seconde modélisation qui utilise les mêmes paramètres que la précédente, mais cette fois ci en conduisant le calcul XFEM (qui est alors un calcul de mécanique pur) avec l’opérateur MECA_STATIQUE au lieu de l’opérateur STAT_NON_LINE dans le cas HM-XFEM. Les résultats obtenus sont strictement identiques aux précédents.
Synthèse des résultats#
Pour chacun des éléments finis HM-XFEM testés, les résultats concordent avec la solution analytique. Pour les éléments HM-XFEM, les fonctionnalités suivantes sont désormais validées:
MODI_MODELE_XFEM,
application d’une pression répartie,
POST_CHAM_XFEM,
le calcul des matrices et vecteurs élémentaires dans le cas où les degrés de liberté Heaviside mécaniques ont été introduits dans les développements du modèle de couplage HM.