Problème modèle et résolution classique à cinq points

Modèle d’étude

Le modèle considéré tout au long cette thèse, à l’exception du chapitre §6, est issu d’un modèle plus général appelé Dead Oil isotherme [47], lui-même un cas particulier d’une famille plus vaste de modèles de milieux poreux [19, 22]. 2.1.1 Équations intérieures Dans le réservoir d’hydrocarbures, représenté par un domaine Ω ⊂ R 2 ouvert borné connexe et régulier par morceaux, le milieu poreux est une roche solide dans laquelle il y a des pores permettant à un fluide de circuler. Ces pores peuvent être plus ou moins nombreux et de tailles différentes. Au lieu de les décrire au niveau microscopique, le milieu est caractérisé par un champ macroscopique appelé porosité φ ∈ L∞(Ω; [0, 1]). Celle-ci exprime localement le rapport entre le volume des pores et celui du milieu. Par la suite, la porosité φ ∈ [0, 1] sera supposée uniforme en espace.Quand la roche est saturée par un fluide, elle possède une aptitude à le laisser circuler dans les pores. Cette capacité, appelée perméabilité absolue, est en général mesurée par un champ de tenseur symétrique κ ∈ L∞(Ω; Sym(R 2 )). Les variations en espace de la perméabilité sont dues aux hétérogénéités et aux changements de type de roche dans le domaine. En milieu homogène et isotrope, hypothèse adaptée ici, ce champ de tenseur se réduit à un scalaire κ > 0. L’unité de mesure courante pour la perméabilité est le darcy (D), avec 1 D = 9.869233 × 10−13 m2 . Le fluide diphasique est composé de l’« huile » indicée par o (oil) et représentant les hydrocarbures, ainsi que de l’« eau » indicée par w (water) et représentant le liquide injecté artificiellement pour pousser l’huile dans les procédés de récupération secondaire. L’écoulement est supposé immiscible, c’est-à-dire que les deux phases ne peuvent pas se mélanger. Pour chaque phase α ∈ {w, o}, il y a en outre : • Deux champs inconnus qui sont la saturation sα(x, t) et la vitesse de filtration uα(x, t). La saturation sα est la proportion volumique de la phase α, ce qui impose sα ∈ [0, 1]. Les équations régissant ces deux champs seront données plus loin.• Deux propriétés pétrophysiques qui sont la viscosité dynamique µα et la perméabilité relative κr,α. Pour un mélange immiscible, µα > 0 est une constante. L’unité de mesure courante pour la viscosité dynamique est le poise (P) avec 1P = 0.1Pa · s. Quant à la perméabilité relative κr,α, il s’agit d’une fonction croissante de la saturation de la phase concernée et traduit l’impact d’une phase sur l’autre. Plusieurs lois empiriques existent à cet égard [49, 72, 93]. Pour la thèse, la plus simple est sélectionnée, celle de Brooks-Corey [16] qui s’écrit κr,α(sα) = κ ] r,α s mα α , (2.1) avec a priori un exposant mα ≥ 1 par phase. Souvent, les deux quantités sont égales mw = mo mais dans certaines simulations ce n’est pas le cas.

Formulation par phases

 Si (0, T) est un intervalle temporel, avec T > 0, l’écoulement a lieu selon le système φ ∂tsα + div uα = qα, dans Ω × (0, T), (2.2a) uα = −κµ−1 α κr,α(sα)∇p, dans Ω × (0, T), (2.2b) so + sw = 1, dans Ω × (0, T), (2.2c) qui compte cinq équations pour cinq inconnues que sont les deux saturations sα, les deux vitesses uα et la pression p. Le champ de pression p(x, t) est supposé commun aux deux phases. Dans d’autres modèles, il convient de distinguer deux pressions, une par phase, au prix d’une relation algébrique supplémentaire fournissant la pression capillaire (différence entre les deux pressions) en fonction d’une saturation. Cela engendre d’autres difficultés numériques [18, 40]. Les deux premières équations (2.2a) expriment la conservation du volume (donc aussi de la masse, puisqu’ici les phases sont supposées incompressibles par l’absence des masses volumiques). Les deux équations suivantes (2.2b) correspondent à la loi de Darcy-Muskat [73], laquelle généralise au cas diphasique celle découverte par Darcy [33] dans le cas monophasique. La dernière équation (2.2c) stipule que tout le volume des pores est rempli par le fluide. Les termes sources qα aux seconds membres de (2.2a) représentent les débits liés à l’injection et à la production au niveau des puits. Ils sont plus faciles à préciser dans la formulation alternative suivante, avec laquelle nous travaillons dorénavant. 

Formulation par flux fractionnaire 

Les variables saturation et pression sont fortement imbriquées dans (2.2). Il est judicieux de « desserrer » un peu leur couplage au moyen d’un changement de notations et de variables [19]. Soit d’abord s = sw, puis la vitesse totale u = uw + uo et le débit total q = qw + qo. D’autre part, la mobilité λα de la phase α ∈ {w, o} est définie par λα(sα) = µ −1 α κr,α(sα), (2.3) de sorte que uα = −κλα(sα)∇p. En sommant les équations de conservation (2.2a) sur α et en invoquant (2.2c), le système (2.2) devient équivalent au système 

Analyse de stabilité en 1-D 

Comme mentionné à la fin de §2.1.1, le système (2.4), (2.9), (2.11) peut contenir des instabilités intrinsèquement dues au couplage entre p et s. Ces instabilités se produisent lorsque le rapport de mobilités M dépasse un certain seuil critique. Ainsi, M apparaît comme une mesure de la raideur du problème. La littérature regorge d’analyses théoriques corroborant cette assertion sur des modèles diphasiques voisins, y compris en 2-D voire 3-D axisymétrique [24,26,50,54,57,59,61,80,81, 99]. Néanmoins, quelques calculs vraiment rigoureux seront mis en avant relatifs au modèle (2.4), (2.9), (2.11). Ces calculs sont inspirés de ceux de Chavent et Jaffré [19, §I.V.2, 39– 41] pour les modèles miscibles 1 1-D. Leur transposition au cas immiscible semble n’avoir jamais été entreprise.

Schéma de référence 5P 

Après ce détour en 1-D, le problème modèle (2.4), (2.9), (2.11) en 2-D est étudié. Cette section présente le schéma numérique institué comme référence pour la comparaison avec les nouveaux schémas développés plus tard. Il s’agit d’une méthode de résolution couramment utilisée en simulation de réservoir. Elle combine deux ingrédients, à savoir : — une discrétisation en temps implicite par rapport à la pression, explicite par rapport à la saturation (IMPES) ; — une discrétisation en espace qui en langage moderne est qualifiée de volumes finis (VF) avec un flux à deux points (TPFA), mais qui par héritage historique est surnommée schéma à cinq points (5P) par les ingénieurs. 2.3.1 Discrétisation en temps Sur l’intervalle temporel de simulation [0, T], avec T > 0, une suite finie croissante de réels positifs est considérée telle que 0 = t 0 < t1 < . . . < tN = T. La durée ∆t n = t n+1 − t n entre deux instants consécutifs est le pas de temps courant. À chaque instant t n , les champs exacts (s(·, tn ), p(·, tn ),u(·, tn )) de la variable spatiale x ∈ Ω sont approximés par des champs approchés (s n (·), pn (·),u n (·)). 

Discrétisation en espace 

Le domaine est restreint à être maillé par des rectangles uniformes. À première vue, cela peut paraître « archaïque » à l’heure où d’abondantes méthodes numériques ont fait leur preuve en maillages généraux non-structurés et où IFPEN a beaucoup investi dans certaines d’entre elles [2, 53, 97, 100]. Cela étant, il faut savoir qu’en ingénierie de réservoir, la plupart des simulations tournent sur ce type de maillages car ce sont les mieux adaptés aux expériences envisagées. C’est d’ailleurs seulement en maillage cartésien rectangulaire que la question de l’effet d’orientation a un sens ! Maillage et notations À la place de notations habituellement employées en volumes finis [42], celles à deux indices provenant des différences finies sont adoptées, celles-ci se prêtant mieux à l’analyse d’erreur des chapitres suivants. Le domaine Ω ⊂ R 2 est supposé être un rectangle de longueur Lx > 0 et de largeur Ly > 0. Il est découpé en une grille régulière de Nx × Ny cellules rectangulaires, où Nx et Ny sont deux entiers positifs. Chaque cellule ou maille est donc un rectangle de longueur ∆x et de largeur ∆y, avec ∆x = Lx Nx , ∆y = Ly Ny . Elle est identifiée par un couple d’indices (i, j) ∈ {1, . . . , Nx} × {1, . . . , Ny}. Le centre de la maille Mi,j est noté xi,j = (xi , yj ). Sa surface vaut |Mi,j | = ∆x∆y. (2.48) Les relations entre les coordonnées des mailles sont xi+1 − xi = ∆x pour tout i ∈ {1, . . . , Nx − 1} et yj+1 − yj = ∆y pour tout j ∈ {1, . . . , Ny − 1}. Maintenant (cf. figure 2.11), les abscisses des arêtes verticales sont définies par x1/2 = x1 − ∆x 2 , xi+1/2 = xi + xi+1 2 , xNx+1/2 = xNx + ∆x 2 36 2.3. Schéma de référence 5P Mi,j • xi+2,j−1 • • ∆x ∆y yj−1 yj−1/2 yj+1/2 xi−1/2 xi+1/2 xi+2 ni+5/2,j−1 ni+2,j−1/2 Figure 2.11 – Entités définies sur un maillage rectangulaire uniforme. −Fi−1/2,j Fi+1/2,j Fi,j+1/2 −Fi,j−1/2 Mi,j Figure 2.12 – Bilan de flux sur la maille Mi,j dans le schéma 5P. 37 Chapitre 2. Problème modèle et résolution classique à cinq points pour i ∈ {1, . . . , Nx − 1}, ainsi que les ordonnées des arêtes horizontales par y1/2 = y1 − ∆y 2 , yj+1/2 = yj + yj+1 2 , yNy+1/2 = yNy + ∆y 2 pour j ∈ {1, . . . , Ny − 1}. Soient xi+1/2,j = (xi+1/2 , yj ) les centres des arêtes verticales pour i ∈ {0, . . . , Nx}. À chaque point xi+1/2,j , un vecteur normal unitaire ni+1/2,j dirigé horizontalement et orienté de gauche à droite est attaché. De même, soient xi,j+1/2 = (xi , yj+1/2 ) les centres des arêtes horizontales pour j ∈ {0, . . . , Ny}. À chaque point xi,j+1/2 , un vecteur normal unitaire ni,j+1/2 dirigé verticalement et orienté de bas en haut 2 est attaché. Les inconnues (s n i,j , pn i,j ) sont localisées au centre xi,j de chaque maille Mi,j , qui peuvent être vues comme une approximation de la moyenne de (s n , pn ) sur Mi,j (en volumes finis) ou une approximation de la valeur ponctuelle de (s n , pn ) en xi,j (en différences finies). Résolution de l’équation en pression En intégrant (2.43a) sur Mi,j , il vient Z Mi,j − div(κλ(s n )∇p) dx dy = Z Mi,j q n+1 dx dy. Par la formule de Green-Ostrogradski, l’intégrale au premier membre est transformée en une intégrale sur le bord de maille ∂Mi,j . En décomposant ce bord rectangulaire en ses quatre arêtes, cette équation devient − Z yj+1/2 yj−1/2 κλ(s n )∇p · ni+1/2,j (xi+1/2 , y) dy + Z yj+1/2 yj−1/2 κλ(s n )∇p · ni−1/2,j (xi−1/2 , y) dy − Z xi+1/2 xi−1/2 κλ(s n )∇p · ni,j+1/2 (x, yj+1/2 ) dx + Z xi+1/2 xi−1/2 κλ(s n )∇p · ni,j−1/2 (x, yj−1/2 ) dx = Z Mi,j q n+1 dx dy. (2.49) Il s’agit du bilan de flux vérifié au niveau continu en espace sur chaque maille Mi,j . L’idée du schéma 5P est de stipuler une version discrète de (2.49) sous la forme Fi+1/2,j − Fi−1/2,j + Fi,j+1/2 − Fi,j−1/2 = ∆x∆y qn+1 i,j , (2.50) dans laquelle les flux Fi±1/2,j et Fi,j±1/2 sont à proposer. Il est intéressant de constater qu’un flux intérieur appparaît deux fois, une fois dans le bilan d’une certaine maille et une autre fois dans celui d’une maille voisine, avec le signe opposé. Cette structure garantit automatiquement la conservativité locale (cf. figure 2.12). Le flux le plus simple est celui dit « à deux points » (TPFA, pour two-point flux approximation). En maillage 2-D cartésien, il est aussi appelé « schéma à cinq points »