Mise en équation et condition expérimental 

Rugosité ou équation de Hurst 

Notre étude numérique est basée sur la génération d’une surface rugueuse synthétique que l’on souhaite la plus proche possible de la réalité. Les mesures de la topographie 3D d’une surface rugueuse peuvent être menées par différentes procédures à l’aide de la technologie laser en laboratoire ou sur le terrain. [4]
Dans le repère (0, x,y,z) où l’axe z est perpendiculaire au plan principal de la fracture (x,y), une surface d’équation z=f(x,y) est dite auto-affine si le changement d’échelle d’un coefficient quelconque λ est isotrope dans le plan et anisotrope selon z, (1bis).

Equation de Fourier 

Les systèmes géothermiques puisent ou rejettent la chaleur du sol par conduction. Le phénomène de conduction se produit lorsque deux points sont à des températures différentes.
Le même phénomène se produit entre le bâtiment et le sol. Le bâtiment, par sa température plus élevée, transfert de la chaleur au sol.

RESULTATS ET INTERPRETATION

Champ d’ouverture aléatoire et corrélé

La première étape de l’étude numérique consiste à générer un champ de perméabilité corrélé selon la loi de puissance avec un exposant de Hurst. L’algorithme de génération du champ est formulé dans un code Matlab fourni en annexe. Il est subdivisé en un programme principal et une fonction générant les champs d’ouverture et de perméabilité aléatoires et corrélés.
Le champ de perméabilité servira de base au Matlab pour résoudre l’écoulement hydrothermale.
Il est donc essentiel d’expliciter le passage d’une génération d’une surface rugueuse auto-affine à celle d’un champ de perméabilité auto-affine le long d’un plan. C’est l’analogie effectué en l’équation (4) et la loi de Darcy soit :
Cette relation constitue l’équivalence remarquable entre des variables de mécanique des fluides telle que la perméabilité ou la viscosité et celles du milieu poreux comme l’ouverture de la fracture.
La synthèse du champ d’ouverture est basée sur une approche classique de simulation des fractures à géométrie complexe développée. Le milieu fracturé est modélisé par un réseau de failles dont la géométrie de chacune se résume à deux plaques parallèles mais avec les mêmes propriétés que le milieu réel.
On définit le paramètre nommé champ ouverture a(x,y) comme l’écart entre les surfaces de la fracture se faisant face. Dans le cas de référence on a un champ d’ouverture constant a(x,y) =
A, il s’agit du point de départ de notre modèle. La surface de fracture est découpée selon une grille de dimension (nx, ny). Les paramètres nx et ny correspondent au nombre d’éléments de la matrice d’ouvertures de fracture. Ils sont choisis pairs pour améliorer la corrélation spatiale avec la loi de puissance de Hurst.
Pour simuler une topographie auto-affine, on ajoute une fluctuation aléatoire d’amplitude noté Z(x,y) au champ d’ouverture A. On assimile l’ensemble des amplitudes aléatoires Z(x,y) à un bruit blanc dans le programme Matlab.
Cette perturbation Z(x,y) est ensuite corrélée spatialement selon la transformation d’échelle (1) et reliée ainsi directement au coefficient de rugosité ζ. En pratique cette corrélation statistique spatiale s’effectue dans le domaine de Fourier à deux dimensions et se résume à la multiplication du spectre blanc par le facteur ||k|| -1 – ζ avec k, le nombre d’onde. Nous réduisons alors à deux paramètres ζ et Z la génération du champ d’ouverture pour une fracture de taille (lx ,ly).
On procède ensuite à un lissage du spectre en supprimant des données hautes fréquences dépassant la fréquence spatiale critique de filtrage π/ lc. Les travaux précisent que la longueur lc joue sur la divergence de la loi de Hurst à petites échelles [2]. Ainsi la longueur caractéristique doit être comprise entre nx et dix fois cette valeur. Le filtre passe-haut n’est rien d’autre que l’application brutale d’une fonction porte dans le domaine de Fourier.
Le champ d’ouverture aléatoire corrélé est obtenu par un retour dans le domaine spatial par une transformée 2D inverse de Fourier. On n’en prend que la partie réelle dans le programme Matlab pour éviter tout problème de déphasage. Le champ de perméabilité de l’écoulement est déduit du champ d’ouverture Z dimensionné et filtré, et de l’équation (4). On conclut que la perméabilité de la fracture dépend ainsi indirectement du coefficient de rugosité ζ i.e. de la géométrie complexe de la fracture.

Visualisation du champ d’ouverture

La visualisation des champs synthétiques d’ouverture et de perméabilité permet d’appréhender qualitativement l’influence du coefficient de rugosité ζ. Une étude statistique complémentaire avec la génération d’un certain nombre de surfaces synthétiques serait nécessaire pour étudier l’influence du paramètre sur le champ d’ouverture géométrique de manière quantitative.
Dans la figure 3, sont représentés les champs d’ouverture et de perméabilité corrélés pour ζ = 0,5 valeur particulière pour des grès fracturés et ζ = 0,8 celle prise pour la plupart des types de roche.
A noter qu’au-delà de la visualisation de l’effet de la rugosité sur une fracture et son champ de perméabilité, la figure 4 illustre la simulation des surfaces rugueuses de deux types de roches différentes, le calcaire fracturé et le granite fracturé. Dans le cas de ce réservoir on s’intéresse exclusivement au socle granitique. Le champ d’ouverture pour une faille dans un socle granitique tend logiquement à s’approcher du cas de référence puisque ζ grès < ζ granite < 1.
Ceci explique la plus grande présence homogénéités (i.e. des zones de valeurs constantes) locales étendus dès les champs d’ouverture et perméabilité de la fracture.
Il existe une visualisation en 3D du champ d’ouverture. Dans la figure 8 est représentée la surface d’ouverture dans un socle granitique. Cette représentation rend mieux compte de la complexité de la géométrie d’une fracture auto-affine mais elle n’apporte pas d’informations supplémentaires majeures par rapport à une visualisation 2D de la fracture rugueuse dans le cadre de notre étude numérique.

Caractérisation de la rugosité

Après la visualisation de plusieurs champs corrélés d’ouverture, il apparaît comme nécessaire de nous assurer de l’isotropie de la loi de Hurst dans le domaine de Fourier. Un programme « valid.m » regroupe les différentes étapes de validation de l’algorithme principal de génération des surfaces et est fourni en annexe. Dans un premier temps on vérifie qu’il n’existe pas d’effet d’étalement d’échelle dans le domaine de Fourier. Les deux composantes du nombre d’onde doivent s’égaliser afin d’éviter une anisotropie de la loi de Hurst. On n’observe aucune dissymétrie apparente entre les expressions des axes kx et ky dans le code de la fonction de génération des deux champs. Pour nous assurer graphiquement de l’absence d’étalement, est associé dans la figure l’équation d’un cercle avec la représentation des valeurs des axes kx et ky . S’il y avait une différence d’échelle entre les deux composantes, les deux figures ne devraient pas se superposer. Comme on le constate sur la figure, il n’y a de phénomène d’étalement.
N.B : on pourrait utiliser une méthode similaire de superposition afin de assurer de la décroissance circulaire et donc isotrope de la loi de Hurst. Néanmoins le spectre corrélé de la loi de Hurst s’atténue très rapidement, ce qui ne facilite pas la définition de l’équation d’un cercle à superposer au spectre. Afin de s’assurer de la précision de la corrélation, le spectre blanc (i.e la phase aléatoire) d’ouverture corrélé doit suivre une loi de puissance ||k|| -1 – ζ soit une loi linéaire en échelle logarithmique selon plusieurs profils 1D suivants l’axe kx, l’axe ky ou obliques. La fonction « polyfit » est utilisée pour l’interpolation du logarithme de la transformation auto-affine des variables kx et ky avec le spectre blanc corrélé.
La variable d’erreur d’un facteur a dû être introduite dans le code Matlab afin de superposer au mieux l’allure des deux droites. L’influence de ce facteur est illustrée en figure 13. Ici on ne constate qu’une valeur d’erreur de préfacteur de 2 est nécessaire pour une estimation correcte de la loi de rugosité. Cette erreur est principalement due à la non-normalisation du champ d’ouverture en annexe associée au programme «valid.m ». La génération des champs par la fonction « gen_surf.m » a été finalement retenue. Cet écart avec la droite d’interpolation n’est plus à prendre en compte car la fonction contient une étape de normalisation. La recherche de cette erreur a donc permis de mettre en évidence l’importance de la normalisation dans la caractérisation du coefficient de Hurst.
Il existe d’autres méthodes de caractérisation du coefficient de Hurst. Ces méthodes alternatives sont explicitées dans des travaux [6]. Comme pour notre étude numérique, l’estimation de la rugosité est basée sur la comparaison du coefficient de Hurst fixé comme variable d’entrée par l’utilisateur et celui extrait de profil 1D par interpolation linéaire de la loi de puissance avec le spectre d’ouverture synthétisé. La différence résulte du fait que la méthode de génération des profils de fracture s’effectue dans le domaine de Fourier 1D.
Par exemple la corrélation par la méthode des écarts types de la surface avec le facteur de rugosité consiste à estimer le champ moyen de profils de fracture 1D. Chaque profil étant constitué de N points découpant un axe x dont chaque point séparé d’une distance Δx et associé à une ouverture ΔL. Pour un Δx compris entre la taille de l’incrément et la longueur de l’axe, l’écart type de ΔL est calculé. Dans ce cas, la loi de Hurst s’écrivant Δx ζ,
La rugosité extraite de profil 1D ne sera pas affectée d’erreur de préfacteur i.e d’écart à la loi de puissance avec l’exposant de Hurst.

DISCUSSION

On observe bien des variations spatiales du champ de pression différentes de celles d’un écoulement laminaire classique obtenu avec une perméabilité constante et un gradient de pression imposé à l’entrée de la fracture. Les faibles fluctuations spatiales dans le champ de vitesse de fluide sont visibles pour des valeurs de coefficient de Hurst de 0,4 et 0,6. Ces zones locales de plus faibles pressions sont reliées aux zones de plus fortes vitesses de fluide et marquent ainsi les chenaux préférentiels de passage du fluide. On constate également que les isobares tendent à se resserrer pour des valeurs de coefficient de Hurst de 0,4 et 0,6 ; ce qui traduit l’augmentation locale du gradient de pression. Une étude statistique serait nécessaire pour en déduire que la rugosité a pour effet d’accélérer localement le flux hydraulique dans la fracture auto-affine. Ces derniers ne correspondent pas directement aux zones de la surface où la fracture est la plus ouverte c’est à dire aux extrema du champ d’ouverture. Cette constatation est non intuitive et témoigne peut-être d’un phénomène de tortuosité de la fracture. Ce phénomène devrait être d’avantage visible dans des grès (modélisés avec ζ = 0,5) dont la contribution des pores n’est pas négligeable par rapport à celle réseau de fractures.
L’augmentation de la tortuosité d’une roche due à la présence d’un gradient de pression fait diminuer localement la perméabilité locale dans une fracture. Le coefficient de Hurst semble diminuer progressivement l’avancé spatiale des isothermes le long de la fracture.
Une étude statistique complémentaire est nécessaire pour en déduire que la rugosité a pour effet de retarder l’échange thermique entre le fluide et la roche mais ceci semble cohérent avec les résultats.
L’effet de la rugosité est plus visible dans les variations du champ de pression que dans celles du champ de température dans le cadre de notre étude. On peut également noter la présence d’une valeur aberrante dans l’échelle des températures de 3,49.10² K (soit 76,69 °C), cela signifierait que le fluide s’est refroidi au contact d’une roche plus chaude. Il en va de même avec des coefficients de Hurst de 0,4 et 0,6.
La résolution de l’écoulement hydrothermale peut être optimisée pour prendre mieux en compte d’autres phénomènes physiques non négligeables et la réalité du contexte géologique d’un réservoir géothermique.
Tout d’abord la méthode itérative de couplage entre le mode de calcul du flux thermique et hydraulique doit être améliorée dans cette étude. On peut construire une boucle dont le critère d’arrêt correspond à un certain nombre d’itération ou des valeurs critiques de température et de pression limites du modèle. Par exemple le calcul s’arrêterait si l’une des températures de zone du maillage dépasse celle de la roche avoisinante, phénomène physique aberrant en réalité.
A noter qu’il fut dans un premier temps, question de résoudre d’avoir le champ de pression résolu avant le calcul du champ de température. On constate alors une instabilité quasiinstantanée du flux thermique. Ce phénomène impose alors un transport de chaleur par convection forcée trop important dès les premiers temps de calcul.
La prise en compte de la dimension de la fracture a également son importance sur l’écoulement hydrothermal.
L’effet du rapport entre la longueur et la largeur de la fracture a déjà été mis en évidence dans d’autre travail [2]. Une augmentation de la longueur d’une fracture aurait u n effet contraire à la rugosité en diminuant l’effet de chenalisation hydraulique. Cependant l’absence d’étude de cet effet de changement d’échelle marque l’une des limitations de notre modèle. Ce travail étant effectué par modélisation, donc il faut garder en mémoire que c’est un travail statistique effectué à partir des données fournis (théorie) donc une grande précision des valeurs de conductivité thermique ne peut donc pas être assurée, elles sont d’ordre indicatif et ne remplacent, en aucun cas une étude de faisabilité locale. Par rapport aux autres études et modélisations faite, notre travail se focalise plutôt sur un écoulement hydrothermal dans une faille alors que dans des autres rapports ils se focalisent plus sur la couche poreuse. L’avantage de notre travail est d’avoir plus de précision sur la fracture car la présence de ces dernières conductrices favorise l’écoulement de fluides et la diffusion de pressions ; ce qui rend le réservoir économiquement avantageux ou risque selon le cadre d’une exploitation pétrolière, minière ou géothermique (accroissement de la zone de drainage ou arrivée précoce d’eau).

CONCLUSION

Un modèle numérique de l’écoulement hydrothermale dans une fracture horizontale rugueuse a été proposé afin de comprendre l’influence de la géométrie de l’ouverture sur l’échange de chaleur entre un fluide newtonien froid et une roche faillé à température constante et homogène soumise en entrée à un faible gradient de pression. Dans le cadre de l’étude, l’écoulement est en régime permanent, laminaire, à faible nombre de Reynolds ; le fluide injecté possède des viscosités et densité peu sensibles aux variations de pression et de température, la composante de sa vitesse normale au plan de fracture est négligée.
Toutes ces hypothèses fortes ont permis de simplifier au maximum les lois de conservation de la masse et de la chaleur de l’écoulement.
Le champ d’ouverture d’une fracture rugueuse est généré dans un algorithme de Matlab. Le coefficient de Hurst contrôle la rugosité de la fracture. Le pré-facteur mis en évidence contrôle la précision de la validation de la loi de puissance associé au coefficient de Hurst que l’on peut percevoir l’amplitude moyenne de l’ouverture. Le champ de perméabilité associé à cette surface a pu en être déduit par l’identification de la loi de Darcy à l’intégration du champ de vitesse de l’écoulement selon l’axe de l’ouverture.
La visualisation des champs de température et de pression a permis de mettre en évidence le phénomène de chenalisation de l’écoulement déjà repéré expérimentalement. La rugosité de la fracture induit une localisation du flux hydraulique et thermique selon des chenaux préférentiels d’écoulement et semble ainsi diminuer l’efficacité des échanges thermiques.
Ces résultats sont à relativiser vu la précision grossière du maillage dû au coût de temps de programmation de l’assignement des valeurs de perméabilités de (n.m)² où n et m représentent le nombre de nœuds sur la longueur et la largeur de la maille. Le couplage des équations ne permet pas encore une stabilité du modèle sans contraindre des conditions limites en pression.
La géométrie et la prise en compte de la dimension de la fracture peuvent être encore améliorées afin de se rapprocher au mieux de la structure géologique du réservoir. D’autres phénomènes de couplage avec la viscosité et la densité du fluide peuvent jouer sur les phénomènes de chenalisation et donc directement sur l’efficacité de l’échange thermique en profondeur duréservoir géothermique.

Table des matières

REMERCIEMENT
SOMMAIRE
LISTE DES FIGURES
LISTE DES TABLEAUX
INTRODUCTION
I. GENERALITE SUR LA GEOTHERMIE
II. ETUDE DE L’ECOULEMENT HYDROTHERMALE DANS LES FRACTURES
III. RESULTAT ET INTERPRETATION
IV. Discussion
CONCLUSION
BIBLIOGRAPHIE
ANNEXES
TABLE DES MATIERES
RESUMES