Problématique des écoulements diphasiques en régime « défavorable »

Problématique des écoulements diphasiques en régime« défavorable »

Cette thèse s’inscrit dans le contexte de la simulation dynamique de réservoir, pour laquelle IFP Energies nouvelles développe depuis de nombreuses années le logiciel PumaFlow. Il s’agit de prédire, à une échelle temporelle de quelques années et à une échelle spatiale de quelques kilomètres, l’écoulement à travers le sous-sol poreux d’un mélange composé d’hydrocarbures déjà présents et de produits artificiellement ajoutés (eau, gaz, polymères, surfactants, mousses). Outre les paramètres physiques du réservoir à connaître, la simulation prend également en compte le procédé de récupération mis en œuvre.

Procédés de récupération des hydrocarbures 

Les méthodes employées pour récupérer les hydrocarbures sont relativement sophistiquées et sont regroupées en trois familles. La première, appelée récupération primaire, est sans doute celle qui se conforme le mieux à l’image populaire : il suffit de forer au bon endroit pour faire jaillir le pétrole. Dans cette catégorie de méthodes, l’extraction se fait par déplétion naturelle, grâce à la différence de pressions entre le réservoir et la surface. Malheureusement, le taux de récupération des procédés primaires est assez faible. Dans le cas d’un gisement d’huile, par exemple, il ne dépasse pas 25%.Pour entraîner les hydrocarbures restés dans le sous-sol après la phase primaire vers la surface, une idée consiste à injecter un fluide non-miscible avec les hydrocarbures et, si possible, peu onéreux tel que de l’eau ou du gaz naturel. Cette injection permet de maintenir la pression dans le gisement et de drainer l’huile vers les puits de production. Cette catégorie de méthodes, appelée récupération secondaire, présente un taux de récupération de 25% à 40%. Son principal défaut réside dans la « fiabilité » de la poussée. En effet, l’eau étant souvent moins visqueuse que l’huile, elle a tendance à être bien plus mobile 1 . Le fluide poussant a donc envie d’aller plus vite que le fluide poussé et peut finir par « percer » le front de propagation, engendrant ainsi des perturbations sous forme de « doigts », comme l’illustre la figure 1.1(a). Ce phénomène porte le nom de digitations visqueuses (en anglais : viscous fingering). Son apparition est désastreuse pour la production car on ne produit presque que de l’eau au niveau des puits producteurs. Comme alternative à l’eau ou le gaz dans le rôle du fluide poussant, il peut être envisagé un gaz miscible (CO2) ou des produits chimiques (polymères, surfactants, mousses). L’avantage de ces derniers provient bien entendu de leur plus grande viscosité, ou de manière équivalente leur plus petite mobilité. Le fluide poussant ayant moins tendance à se mouvoir, le front de déplacement avec l’huile est plus stable et les digitations visqueuses n’apparaissent pas. Cette famille de méthodes, qualifiée de récupération tertiaire ou récupération améliorée (en anglais : enhanced oil recovery, EOR) élève notablement le taux de récupération, l’augmentation variant entre 30% à 70%. Elle coûte néanmoins cher en raison du prix des produits chimiques. Aussi, afin de rentabiliser le procédé, l’injection de polymère n’est effectuée que pendant un laps de temps limité, après quoi de l’eau de chasse est envoyée pour pousser le polymère, comme le décrit la figure 1.1(b). Mais de nouveau, il existe un rapport de mobilités défavorable entre l’eau et le front arrière du bouchon de polymère. 

Modèles physiques 

L’un des buts de la simulation de réservoir est de prédire la durée pendant laquelle l’huile pourra être recueillie et surtout à quel moment l’eau arrivera aux puits de production. Toute incertitude dans cette prédiction entraîne des pertes financières importantes se chiffrant en centaines de milliers de dollars. À cette fin, il existe une grande variété de modèles physiques qui se démarquent entre eux par les processus physico-chimiques pris en compte . La plupart de ces modèles sont par ailleurs formulés à une échelle macroscopique, après une étape d’homogénéisation préalable. Cela implique qu’ils ne sont pas aptes à décrire certains phénomènes intrinsèquement microscopiques comme les digitations visqueuses. Dans ces dernières, en effet, les allures des doigts dépendent des propriétés microscopiques fines telles que la forme des chenaux ou des capillaires dans le milieu poreux, lesquelles ont été « gommées » au cours de l’établissement du modèle macroscopique. Par conséquent, il convient d’être prudent quant à l’interprétation de ce qui semble être des « digitations » dans certains résultats numériques obtenus avec un modèle de ce type. Actuellement, les modèles les plus évolués en simulation de réservoir sont multi-constituants (plusieurs espèces), polyphasiques (chaque espèce pouvant se présenter sous plusieurs phases, chaque phase pouvant disparaître ou apparaître selon les lois thermodynamiques), éventuellement réactifs (les espèces réagissent chimiquement entre elles). Ils rendent également compte des effets de capillarité (différence de pressions entre les phases), de dispersion (diffusion des espèces), de compressibilité et de gravité. À l’autre bout du spectre, il y a des modèles, qui bien qu’assez rudimentaires, contiennent déjà l’essence de la difficulté à traiter. C’est le cas du modèle considéré dans cette thèse pour l’élaboration de nouveaux schémas. Il s’écrit comme un système de trois équations aux dérivées partielles u = −κλ(s)∇p, (1.1a) div u = q, (1.1b) φ ∂ts + div(f(s)u) = qw, (1.1c) où les trois inconnues sont : la saturation d’eau s ∈ [0, 1], la pression commune aux deux phases p > 0 et la vitesse totale u ∈ R n . Dans (1.1), les constantes φ (porosité) et κ (perméabilité) sont connues. Les fonctions λ (mobilité totale) et f (flux fractionnaire) sont également spécifiées. En particulier, comme cela sera détaillé en §2.1, le flux fractionnaire f(s) = Msmw Msmw + (1 − s)mo (1.2) fait intervenir le rapport de mobilités M = κ ] r,wµo κ ] r,oµw , (1.3) où µo (respectivement µw) est la viscosité de l’huile (respectivement de l’eau) et κ ] r,o (respectivement κ ] r,w) est la perméabilité relative maximale de l’huile (respectivement de l’eau). Les exposants mo et mw représentent les coefficients de Corey [31] associés à chacune des phases considérées ici. Le système (1.1) doit bien sûr être muni de conditions aux limites appropriées sur un domaine spatial borné, ainsi que d’une condition initiale sur s. 

Difficultés mathématiques 

Le système (1.1) possède une caractéristique tout à fait représentative de la plupart des autres modèles de réservoir. Il est potentiellement instable lorsque le rapport de mobilités M est trop grand. Plus précisément, il y a des éléments de preuve solides [24,26,50,54,57, 59, 61, 80, 81, 99] en faveur de l’assertion suivante : Il existe un seuil critique M] > 0, dépendant de paramètres physiques autres que ceux figurant au second membre de (1.3), tel que (i) si M ≤ M] , alors le problème (1.1) avec les conditions limites est stable (ii) si M > M] , alors le même problème est instable, la notion de stabilité étant celle de la nonamplification des perturbations linéarisées. En §2.2, une analyse de stabilité rigoureuse de (1.1) sera réalisée dans deux cas particuliers où M] peut être explicitement déterminé. En admettant cette conjecture, il est maintenant possible de donner un sens précis aux épithètes « favorable » et « défavorable » concernant le rapport de mobilités : il s’agit simplement de M ≤ M] et M > M] , à ceci près que M] n’est pas toujours facilement calculable et, même en général, indéterminé. La nature instable du problème pour M > M] soulève une interrogation « philosophique » au regard de la résolution numérique de (1.1). En effet, l’instabilité inhérente au système amplifie n’importe quelle erreur d’approximation de la solution, à commencer par des erreurs d’arrondi. Comment garantir, dans ces conditions, que la solution calculée par un schéma soit suffisamment proche de la solution exacte au temps final ? Le cas échéant, jusqu’à quelle tolérance faut-il considérer la prolifération des erreurs numériques liées à l’instabilité sous prétexte qu’elles sont « physiques », c’est-à-dire « innées » au sein du modèle, ou faut-il juste les retarder ou les minimiser, sachant qu’elles finiront fatalement par prédominer ? Pour mieux éclairer ce dilemme, deux exemples élémentaires de problèmes différentiels instables sont examinés. Le premier correspond à la croissance exponentielle régie par dX dt = aX, X(0) = 0, (1.4) avec a > 0. La solution exacte est X(t) = 0 pour tout t ≥ 0. C’est cette solution qui doit être simulée numériquement. Malheureusement, elle est instable puisque toute erreur η 6= 0 sur la donnée initiale (Xη(0) = η) conduit à Xη(t) = exp(at)η. Il faut donc s’accorder sur ce qui est « voulu » vraiment entre les deux options suivantes : 1. Faire tout son possible au niveau numérique pour que X reste au plus proche de 0, et ce le « plus longtemps possible ». Cette option se justifie par l’attachement de l’utilisateur à la solution exacte mais fait fi de son instabilité. Pire, elle introduit la nécessité de lutter contre l’instabilité. 2. Accepter que X s’éloigne indéfiniment de 0 au niveau discret, comme c’est le cas du schéma d’Euler explicite Xn+1 = (1 + a∆t)Xn qui donne Xn η = (1 + a∆t) nη. Cette option se défend en disant qu’après tout, l’instabilité est déjà contenue dans le modèle (1.4) et que le schéma ne doit pas aller contre cette tendance. Mais dans ce cas, le résultat calculé Xn η a-t-il une quelconque signification, dans la mesure où il dépend de la perturbation initiale η, laquelle varie aléatoirement d’une machine à l’autre voire d’une exécution à l’autre ? Toute cette discussion se transpose au second exemple d 2ϑ dt 2 = −ω 2 sin ϑ, ϑ(0) = π, dϑ dt (0) = 0, (1.5) qui est non-linéaire et qui correspond à un pendule avec équilibre instable en position haute. La solution exacte est ϑ(t) = π pour tout t ≥ 0 et c’est elle qui doit être simulée numériquement. De manière analogue au premier exemple, il faut se décider entre deux options. Dans la première option, le pendule est maintenu en haut parce que c’est la solution. Dans la seconde option, la trajectoire numériquement calculée à partir d’une position initiale légèrement erronée ϑ(0) = π − η, avec η 6= 0, fait évoluer le pendule vers la position stable ϑ(+∞) = 0, ce qui semble conforme à la « réalité ». Mais si la trajectoire discrète ϑ n η dépend d’un paramètre de perturbation η aléatoire selon la machine, de quelle trajectoire continue est-elle l’approximation ? Cet exemple est intéressant car dans la « réalité » le passage de la position instable vers la position stable se fait selon une trajectoire également imprévisible.