Présentation des méthodes numériques pour la représentation et la caractérisation du milieu fracturé

Vers des méthodes hybrides plus performantes

Afin de profiter des avantages de ces méthodes, et d’en réduire partiellement les inconvénients (en particulier les difficultés liées à l’imposition de certains types de condition aux limites), des méthodes hybrides (i.e. couplage d’une meshless method avec la méthode FEM) ont été proposées. En d’autres termes, cette alliance entre les deux types de méthode permettra non seulement de retrouver la propriété de Kronecker-_ des fonctions de forme (afin d’imposer facilement les conditions de type Dirichlet), mais également d’augmenter l’ordre d’approximation sans ajouter de noeuds supplémentaires [162]. Par exemple, nous pouvons citer la hp-cloud Method [74], la Generalized Finite Element Method (GFEM) [199] ou encore la Partition of Unity Finite Element Method (PUFEM) [144]. Toutes les méthodes citées précédemment sont basées sur le principe de partition de l’unité. En pratique, cela signifie que, pour un domaine découpé en N sous-domaines I auxquels sont associés une fonction non nulle fI (x) (et nulle partout ailleurs), la relation suivante est satisfaite sur : XN I=1 fI (x) = 1 (2.1) Globalement la hp-cloud Method et la PUFEM sont équivalentes : pour la première les fonctions de forme sont obtenues par moindres carrés alors que dans la seconde, il s’agit de polynômes de Lagrange (plus simples d’utilisation). En particulier, pour la PUFEM, la propriété de Kronecker-_ des fonctions de forme est parfaitement respectée, et l’augmentation de l’ordre d’approximation se fera en considérant un ensemble de fonctions supplémentaires (dites fonctions d’enrichissement), respectant le principe de partition de l’unité, et venant s’ajouter à l’expression de l’approximation FEM classique. Par exemple, l’approximation du champ de déplacement dans le cas général s’écrit : u(x) = X i ai’i (x) + termes d’enrichissement (2.2)

Une variante de la méthode PUFEM est la méthode XFEM (eXtended Finite Element Method). Une des particularités très intéressante de la méthode XFEM est que celle-ci est capable de prendre en compte les discontinuités de type faible et forte. Autrement dit, des fonctions d’enrichissement discontinues seront introduites afin de rendre compte à la fois de la discontinuité du champ primal (uniquement dans le cas des discontinuités fortes), mais également de son gradient, et ce, sans avoir recours à des techniques de remaillage. Ensuite, contrairement à la PUFEM (ou à la GFEM), cet enrichissement est local pour la méthode XFEM. En particulier, l’enrichissement local sera un point fort de cette méthode. Premièrement, cela permettra de réduire le temps de calcul dédié à l’inversion des matrices de rigidité ainsi que pour l’intégration des différentes contributions des formes faibles des équations d’équilibre. Deuxièmement, le nombre d’inconnues supplémentaires dans le problème à résoudre sera réduit, en comparaison de certaines meshless methods comme la GFEM (i.e. enrichissement global). La méthode XFEM, faisant intervenir dans sa formulation des termes d’enrichissement supplémentaires, est couramment rencontrée dans la littérature sous la dénomination extrinsic XFEM. Une variante, identifiée par l’appellation intrinsic XFEM (proposée par Fries et Belytschko [88]), consiste à modifier localement la nature des fonctions de forme grâce à des fonctions obtenues par moindres carrés pouvant représenter la discontinuité de la solution.

L’avantage de cette formulation est qu’elle n’augmente pas le nombre de termes supplémentaires. Néanmoins, en pratique, la mise en oeuvre de l’intrinsic XFEM semble être plus coûteuse que l’extrinsic XFEM du fait de la nature des fonctions d’enrichissement [89]. Dans la suite de ce mémoire, nous ne ferons référence qu’à l’extrinsic XFEM. La méthode XFEM a été proposée initialement par Belytschko et al. [14] et Moës et al. [150]. Dans la littérature, cette méthode a été employée pour modéliser le comportement des trous et des inclusions [201], des jonctions de fractures [62], la croissance d’une fracture cohésive (en 2D [148, 145, 234] et en 3D [100]), mais également dans le cadre des lois de contact-frottement (pour le cas des petits glissements [70, 123] et grands glissements [164]). Récemment cette méthode a également été combinée avec les modèles de couplage hydromécanique [212, 119, 180, 152] (voir également chapitre 3 de ce mémoire). Avec XFEM, la géométrie de la fracture n’est plus dépendante de la structure du maillage (lié au fait que cette méthode est un hybride entre la PUM et la FEM). En d’autres termes, la fracture n’est plus maillée. Mathématiquement, cette méthode se fonde sur l’ajout local de fonctions d’enrichissement [14, 150, 69, 21]. En particulier, la fonction de Heaviside (cf. équations (2.4) et (2.3)) est couramment considérée comme une fonction d’enrichissement afin de rendre compte de la discontinuité du champ de déplacement dans une fine bande d’éléments autour de la fracture. D’autre part, au niveau de la pointe de la fracture, des fonctions singulières sont ajoutées afin de prendre en compte le comportement asymptotique de la pointe (cf. équations (2.5) et (2.3)). À ces fonctions d’enrichissement sont associés des degrés de liberté supplémentaires. L’approximation du champ de déplacement s’écrit donc pour XFEM [148] :

Synthèse et objectifs

Nous avons donc vu que le choix de la méthode XFEM, pour modéliser le comportement d’une fracture, est un bon compromis par rapport à la multitude de meshless methods disponibles dans la littérature. Bien sûr, cette méthode a ses inconvénients (notamment en ce qui concerne l’imposition de conditions aux limites au niveau d’une interface (voir Moës et al. [147] et Béchet et al. [13] et chapitre 3 de ce mémoire), mais est également peu adaptée aux cas impliquant de très fortes déformations du domaine). Cependant sa principale force repose sur sa flexibilité d’association avec d’autres techniques comme les level-sets pour représenter la discontinuité des champs au niveau de la fracture, et ce, indépendamment du maillage. Par exemple, il sera plus aisé de représenter des fractures présentant des géométries complexes (voir par exemple Fries et al. [90] pour le cas d’interfaces courbes pour le cas mécanique sans couplage). Ensuite, l’utilisation d’un modèle de zone cohésive pour gérer la propagation de la fracture, semble également être une bonne alternative. En effet, il sera plus simple de rendre compte de la circulation du fluide au sein de la fracture (et ce quelque soit le régime dominant la propagation), puisque la pression de fluide générée par l’écoulement est directement introduite dans la formulation de la loi cohésive.

Ensuite, ce type de modèle permet de s’affranchir de la nécessité d’introduire des fonctions d’enrichissement singulières pour rendre compte de la discontinuité des champs de déplacement et de pression au niveau de la pointe de fissure. Néanmoins, il faudra être prudent dans la manière de définir ces champs au niveau des éléments finis du maillage intersectés par la fracture. Dans la suite de ce mémoire, nous présenterons en détails le modèle numérique HMXFEM dont nous avons déjà fait mention dans ce chapitre. En particulier, nous détaillerons la loi cohésive non-régularisée (i.e. loi de Talon-Curnier [204]) que nous avons adapté au cas de l’écoulement d’un fluide dans une fracture. Ensuite, dans le chapitre 4, nous montrerons que l’utilisation d’un modèle de zone cohésive est particulièrement bien adapté à la propagation d’une fracture hydraulique (sur chemin pré-défini), puisque différents régimes seront étudiés sans avoir eu besoin de faire de distinctions (e.g. dominance des effets visqueux) dans le formalisme éléments finis.

Cas du massif Du fait de la présence de la fracture, tous les champs associés au massif sont supposés être discontinus de part et d’autre de la fracture. Ici, nous considérons que la discontinuité du champ de déplacement est une discontinuité dite forte. Selon Mosler [155], le champ primal de l’approximation du déplacement ainsi que son gradient associé (i.e. le tenseur de déformation « ) sont tous deux discontinus dû au saut au niveau de la fracture. D’autre part, contrairement au champ de déplacement, le champ de pression p (i.e. pression de pore du massif) doit être continu au niveau de chacune des parois de la fracture du fait de la condition de continuité p = pf sur 􀀀1 et 􀀀2. Ainsi, l’approximation de ce champ est associée à une discontinuité dite faible : le champ de pression de pore est continu, cependant son gradient est discontinu. Ainsi, il est possible de capter la discontinuité du gradient normal du champ de pression de pore lorsque de forts échanges se produisent entre la fracture et le milieu poreux environnant. Il est important de signaler que les fonctions d’enrichissement introduites autorisent une discontinuité forte de l’approximation du champ de pression de pore de part et d’autre de la fracture, mais la continuité du champ primal est imposée par dualisation de la condition au limite p = pf au moyen des 2 multiplicateurs de Lagrange q1 et q2 (cf. équations (3.38) et (3.39)) (également employé par Zilian et Legay [235]).

Dans la littérature, certaines études en mécanique, impliquant des discontinuités dites faibles, ont été conduites. Hansbo et Hansbo [110] proposent d’utiliser l’enrichissement pour représenter une discontinuité forte, mais de forcer la continuité (afin de retrouver une discontinuité faible) au moyen de la méthode de Nitsche [165]. Pour leur part, Sukumar et al. [201] ont suggéré d’introduire la valeur absolue des level-sets comme fonctions d’enrichissement, mais celles-ci s’avèrent peu efficaces en termes de convergence. Enfin, Moës et al. [149] proposent ensuite d’introduire des fonctions d’enrichissement à gradient discontinu (i.e. ridge functions). Le taux de convergence avec ce type de fonctions se voit optimisé. Ainsi, nous avons finalement choisi de considérer une discontinuité faible en forçant la continuité au moyen de multiplicateurs de Lagrange (contrairement à Hansbo et Hansbo [110], puisque les méthodes Nitschiennes nécessitent la détermination d’un paramètre de pénalisation) :

• pour des raisons de stabilité puisque N’Deffo et al. [160] ont montré que cet artefact améliore la précision et le conditionnement dans le cas des éléments quadratiques ;

• et pour des raisons pratiques. En effet, l’impact sur la structure du code élément finis est moindre puisque seulement deux formulations variationnelles (cf. équations (3.38) et (3.39)) sont introduites au niveau de chacune des parois de la fracture.