FR 2863052 A1 20050603 FR 0314211 A 20031202 La présente invention concerne une méthode pour déterminer l'ensemble des composantes 2D ou 3D d'un tenseur de perméabilité absolue d'un échantillon de milieu poreux etlou permettre une mise à l'échelle de champs de perméabilité hétérogènes, à partir de données mesurées en laboratoire ou issues de cartes de perméabilité d'une zone souterraine fournies par les géologues. ETAT DE LA TECHNIQUE La mesure en laboratoire de la perméabilité des roches est un étape incontournable des études de réservoirs pétroliers, d'aquifères, en génie civil, voire de catalyseurs employés par l'industrie chimique. Si le milieu présente des directions d'anisotropie évidentes, l'expérimentateur tâchera de les respecter lors de sa campagne de mesures. En d'autre termes, il se placera d'emblée suivant les axes propres du tenseur de perméabilité. Il y a bien des cas où cette direction n'est pas donnée à l'avance, et il faut donc a priori déterminer le tenseur sans idée préconçue. Il peut aussi s'avérer que les échantillons de roche n'ont pas été carottés dans les bonnes directions; de plus les axes propres du tenseur ne sont pas forcement alignés avec le litage que peuvent présenter certaines carottes. Il est donc particulièrement utilé de disposer de mesures directes de l'ensemble des composantes du tenseur, ou d'au moins en estimer l'éventuelle anisotropie. Une méthode de détermination d'un tenseur de perméabilité d'un échantillon de roche à partir de mesures permettant d'estimer avec suffisamment de précision les pressions moyennes autour d'un échantillon placé dans un perméamètre adapté à ces mesures, est décrit par exemple dans la publication suivante: - Renard P. "Laboratory détermination of the full permeability tensor " JGR 106, B11 2001 pp26 443-26 451. Une autre méthode connue de détermination d'un tenseur de perméabilité dans laquelle on crée dans l'échantillon des écoulements rendus volontairement tortueux par des modifications des conditions aux limites des faces d'injection, et on résout un problème inverse, est décrit par'exemple dans la publication suivante: - Bernabé Y.: "On the measurement of perméability in anisotropic rocks " in "Fault mechanics and transport properties of rocks " edited by B.Evans and TF Wong pp 147-167, Academic San Diego 1992 Pour être mises en oeuvre, les méthodes connues exigent une modification substantielle des conditions d'écoulement qui les rendent onéreuses et difficilement 10 applicables en pratique. Dans le contexte d'un changement d'échelle (upscaling) bien connu des ingénieurs de réservoir, on est amené à remplacer l'expérimentation de laboratoire par les résultats d'une simulation numérique "fine" sur le milieu hétérogène, afin de le remplacer par un milieu homogène équivalent dont le tenseur de perméabilité que l'on note par exemple Keif. Dans le cas d'un milieu anisotrope, la simulation fine de référence se fait en général en considérant des conditions aux limites périodiques. La plupart des auteurs considèrent à tort que ce sont les seules conditions permettant d'obtenir l'ensemble des éléments du tenseur Keff. Comme va le montrer la description qui va suivre, une interprétation adaptée d'une simulation numérique du même perméamètre est à même de fournir cette information. Les conditions aux limites imposées dans un perméamètre, peuvent être plus réalistes dans le cas où l'écoulement à grande échelle est contraint par des barrières argileuses. De plus, la comparaison entre les deux tenseurs résultat peuvent fournir des indications sur l'échelle du Volume Elémentaire Représentatif (VER). LA METHODE SELON L'INVENTION La méthode selon l'invention permet de déterminer, à partir des mesures connues de peiméabilité, toutes les composantes d'un tenseur de perméabilité absolue (Keq,) d'un échantillon de milieu poreux anisotrope. Elle comporte essentiellement les étapes suivantes: - on mesure ou on calcule en imposant deux écoulements de liquide 'successivement suivant deux (ou trois) directions orthogonales (Ox, Oy) (et Oz) au travers de l'échantillon et sous gradient de pression imposé, les débits de liquide et les forces transversales générées dans l'échantillon du fait du passage du liquide, ces quantités étant ramenées à des gradients de pression unitaires; et - on utilise un simulateur pour déterminer de proche en proche, à partir d'une solution numérique a priori d'un écoulement confiné dans un milieu homogène, les composantes du tenseur équivalent (k) par une méthode inverse quelconque tel que les dites quantités correspondant s'ajustent sensiblement aux quantités mesurées ou calculées. Suivant un premier mode de mise en oeuvre, on mesure les quantités (débits et forces transversales) au moyen d'un perméamètre et de moyens de mesure de forces transversales appliquées à l'échantillon en réponse à des écoulements imposés au travers de l'échantillon. Ces forces transversales appliquées à l'échantillon, en réponse à des écoulements imposés à travers lui, peuvent être mesurées par exemple en mesurant les variations de poids de l'échantillon dans le perméamètre quand on alterne le sens de passage des écoulements de liquide imposés au travers de l'échantillon. En répétant cette procédure après avoir effectué une rotation de l'échantillon de 90 autour d'un axe parallèle à l'écoulement, on peut mesurer l'autre composante. Suivant un deuxième mode de mise en oeuvre, on détermine les mêmes quantités par simulation d'écoulements dans un milieu, à partir de données de perméabilité connues ou estimées telles que des cartes de perméabilité fournies par exemple par des géologues dans les applications en géophysique. La mise en oeuvre de la méthode selon l'invention est avantageuse en ce qu'elle permet de simplifier et d'accélérer la détermination du tenseur de perméabilité d'échantillons. 4 PRESENTATION SUCCINCTE DES FIGURES Les caractéristiques et avantages de la méthode selon l'invention, apparaîtront plus clairement à la lecture de la description ci-après d'un exemple non limitatif de mise en oeuvre, en se référant aux dessins annexés où : - la figure 1 montre, un organigramme, permettant une détermination itérative du tenseur équivalent Keq(n) d'un échantillon de roche, à partir de différentes mesures de débit et de forces effectuées en le confinant dans un perméamètre; - les figures 2 à 4 illustrent différentes exemples de simulation montrant la pertinence de la méthode, qui seront commentées dans le cours de la description; et -la figure 5 montre un exemple schématique de perméamètre convenant pour effectuer des mesures sur un échantillon de roche. DESCRIPTION DETAILLEE DE LA METHODE Les mesures de perméabilité utiles à la mise en oeuvre de la méthode qui va être décrite peuvent être obtenues au moyen d'un perméamètre de type connu. Un tel perméamètre comporte par exemple (fig.5), une cellule de confinement 1 pour l'échantillon. La paroi latérale de la cellule est une gaine déformable. A une première extrémité, la cellule 1 est connectée à un circuit 2 d'injection de fluide à une pression d'injection P1. A son extrémité opposée, la cellule est connectée à un circuit 3 d'évacuation du fluide issu de l'échantillon maintenu à une pression P2. La cellule 1 est placée dans une enceinte extérieure 4 connectée à un dispositif 5 pour appliquer une pression de confinement Pc à l'échantillon, de façon que le flux de liquide sur le pourtour latéral de l'échantillon soit nul (bords imperméables). Les embouts opposés de la cellule 1 sont tels que les pressions P1 et P2 respectivement exercées sur les faces amont et aval de l'échantillon sont uniformes. On note dPx,x = P2 Pi De façon à simplifier la présentation, on considère tout d'abord un milieu poreux 2D, carré de côté L. La méthode, transposée aux cas 3 D avec des proportions quelconques sera exposée plus loin. Dans toute cette partie, on considère que l'échantillon de roche poreuse est homogène mais avec un tenseur de perméabilité k dont les axes propres ne coïncident pas nécessairement avec les axes du perméamètre où il est placé (fig. 1). Dans un premier temps, on applique le flux de liquide à l'échantillon orienté suivant une première direction A par rapport à l'axe de la cellule et on fait un ensemble de mesures que l'on va préciser plus loin. Dans un deuxième temps, on recommence l'expérience après avoir orienté l'échantillon selon une direction B orthogonale à la première et l'on fait les mêmes mesures selon le même protocole. Sur les faces de l'échantillon parallèles à l'écoulement imposé, on a les conditions standard de flux nul de Neumann, soit n.k.Vp(r) = 0 (1) où r désigne le vecteur position dans l'espace, de composantes x, y (ou xi, x2) n est le 10 vecteur unitaire orthogonal à la frontière au point considéré, Vp(r) est le gradient de pression. Composantes par composantes, cette formule s'écrit: E n..k... ap(r) = 0, ou bien de façon plus concise encore en adoptant à partir de t,J=1,D axf maintenant la convention d'Einstein de sommation sur les indices répétés: k ap(r) n =0 axi - Equation en pression: L'équation en pression est donnée classiquement, où le point (.) correspond à des contractions de tenseurs: k 0. .Op(r) =0 ^ i soit, in extenso, toujours avec la convention d'Einstein: (2) (3) k. ap(r) axi a ax; L'équation en pression et les conditions aux limites définissent un problème noté problème (P) bien posé possédant une solution unique. Dans le cas particulier où la direction x coïncide avec un axe propre de k, la solution linéaire p(x, y, z) = P + APX,X É L est l'unique solution du problème. On obtient alors dans ce cas la relation classique Qx,x = kx,x x'x dont on tire kx x si Qx x et OPx x ont été simultanément mesurés. Dans le cas anisotrope général, la solution p(x, y, z) ne possède pas de structure simple et il n'y a pas de solution analytique à notre connaissance. On se contente donc de supposer que l'on est dans la situation précédemment présentée et l'on fournit donc en pratique des estimations des valeurs propres du tenseur k. Nous allons montrer maintenant que à condition de mesurer quelques quantités supplémentaires, il est possible de déterminer tous les éléments k du tenseur k. Pour cela, on remarque que l'équation (4) peut se réécrire sous la forme: fera axi ap(r) axi f (r) = O Vf (r) (5) Avec r désignant le vecteur position à D composantes, et dDr l'élément d'intégration. En choisissant alors la fonction f(r) de façon convenable (f(r) =x, y puis z) à 3D, en intégrant par parties l'égalité précédente on obtient le jeu d'égalités suivantes à 2D: xx = kxx Lz f L dy(p(x = L, y) p(x = 0, y)) + kxy Lz f L dx(p(x, y = L) - p(x, y = 0)) soit, avec des notations plus condensées n -k api+k ap r1 xz xx ax x xY ay (7 (6) ap 8x = qxx kYx + ax x (7) où qxx correspond au débit normalisé par le volume de l'échantillon q,x,x = Q,,x /V. Les {... }X désignent une moyenne spatiale et l'indice x indique ici la direction du gradient 20 de pression imposé. Les quantités Byx et Sxy sont définies par: f L y kxxj ôp (L, y) - ap (0, y) ]dy - fo y kxs[ âx (L, y) --- (0, y) ]dy Syx Vgxx fo kxx ax (0, y) dy Cette quantité correspond à la différence des coordonnées y des barycentres des débits en x locaux sur les faces aval et amont de l'échantillon lors d'un écoulement imposé selon x. De même on a: I x k , [ ap (x, L) ô y (x,0) l dx Jo x kYY [ â (x, L) ap (x,0) 1dx &y = J L o kYY â (x,0) Vg3Y y Des formules analogues peuvent être démontrées par la même technique à 3 D, introduisant les trois débits qxx, qyy, qzz, les six quantités (8 x Szx Sxy Sz,, Sxz et (5yz), et les six autres grandeurs additionnelles ( ap ap) ay. ' {al)} azx' axy' {al)} azy' ax lapl ay Z Les conditions de Dirichlet sur les faces amont et aval fournissent bien sûr la relation {a j '. . En considérant une expérience effectuée dans une autre direction ax x L orthogonale, y on obtient: ap 1 aP 4u q = kYx { ax} + k yy a y, y, (8) Y /iSxy qyy =kxx1ap1 Y (9) De façon à alléger les notations, nous convenons de travailler avec des gradients de pression imposés unitaires ap = PxxouYY -1. En inversant ces égalités, on obtient dans axJ 0 L ce cas les relations suivantes: Ç a gxx gYY SxY la y x k,,_ 1 1ax1 Ylayfx agYY rg x( y ( 11) 1 lay xapl ly (10) kYY a kxy = {--P ax kxx + 8xyq Y ap l kyx = 1aY x }k +8yxgxx On peut rassembler ces relations sous une forme symbolique: (12) (14) k = F(qxx, qyy> ap Sap l ayJx'laxJy' x, axy) ' (15) Pour résumer, dans le cas homogène 2D, en mesurant les quantités qxx, gyy,lapl, ax y f ap 1, SyX, Sxy, on peut retrouver directement le tenseur k. Comme le tenseur est défini layJx par 3 nombres à 2D, et 6 à 3D (mais dans ce cas on détermine 15 quantités additionnelles), ceci signifie que des relations implicites existent entre les grandeurs pl ap x 8x. En théorie, si on retient au moins trois valeurs, on doit donc qxx, qyy, x {ax}y' y pouvoir encore retrouver k. Nous allons retenir le sous ensemble qxx, qyy, {ap}x, {â J y (16) Le principal intérêt de ces données est de posséder un sens physique immédiat et partant d'être directement mesurables: que, qyy correspondent aux débits usuels (ramenés 15 au volume de l'expérience). Les quantités apx' azysont mesurables car: (ap = 1 fLdx(p(x, y = L) p(x, y = 0)) n'est autre qu'une résultante de la force jl ay x LZ Jo transverse par unité de volume imposée à l'échantillon. de plus, en utilisant la relation générale: f d rVp(1-) V fad -'r p(r)n, on voit que cette interprétation du gradient de V pression moyen en terme de force surfacique exercée sur l'échantillon subsiste. Toutes ces quantités ou forces peuvent être mesurées directement en laboratoire. Pour mesurer les forces transversales, on peut utiliser par exemple des dynamomètres ou des jauges de contraintes appliquées localement contre les parois latérales de" l'échantillon dans sa gaine (cf. fig.5). Pour mesurer plus commodément la résultante des forces transversales, il est possible également de mesurer les variations de poids de l'échantillon selon que les forces transversales se rajoutent ou se soustraient à son poids statique. Toutes choses égales par ailleurs, on provoque un écoulement au travers de l'échantillon successivement suivant un premier sens et suivant le sens opposé. Dans un premier cas, la force transversale s'ajoute par exemple au poids de la cellule et dans l'autre cas elle se retranche. En soustrayant les mesures, on peut obtenir directement la mesure de la force transversale. Pour la mise en oeuvre, on pose l'enceinte 4 (cf. fig.5) à l'horizontale sur un dispositif de pesée (non représenté). On mesure la variation de poids pour deux positions à 90 l'une de l'autre correspondant aux directions A et B fléchées sur la figure 1. Exemple 2D On suppose que l'on a obtenu comme indiqué le quadruplet que, qyy, {aF' 1ax yde mesures de débits et de forces transversales appliquées à l'échantillon dans les positions successives qu'on lui fait prendre dans le perméamètre (4 mesures globales de forces et de débits à 2D, 9 à 3D). A partir de maintenant, les milieux poreux considérés peuvent être hétérogènes, c'est à dire caractérisés par une carte fine de tenseurs de perméabilité dépendant de la position, k(r) qui peut être fournie par le géologue ou tout autre technique. De façon à remonter au tenseur k que l'on notera Keq à partir de maintenant pour rappeler qu'il s'agit d'un tenseur effectif à partir de la connaissance de 25}ap}X' {aply, il est nécessaire de disposer d'une solution numérique de l'équation de Laplace anisotrope. Toute méthode connue d'inversion peut être utilisée. L'utilisateur peut choisir par exemple une méthode d' "optimisation" (pour minimiser l'écart entre les valeurs de ap mesurées et les prédictions du modèle numérique), une méthode par q,g" 'laP ax}Y essai et erreur, par dichotomie, etc. Un exemple est donné ci-après. De façon à trouver le tenseur Keq, associé au milieu homogène équivalent, on va utiliser un simulateur numérique d'un type connu disponible sur le marché, opérant par éléments ou volumes finis, pour résoudre le problème aux limites 2D (ou 3D). La procédure consiste par exemple à construire une suite de tenseurs Keq(1),..., Keq(n), etc. définie par son premier terme représentant une première estimation de Keq obtenue en utilisant la relation (15) : Keq(1) F(q.,'YY,ap x,apY,O,o), Comme les moments SyXet 8xy sont inconnus, car non mesurés, ceci ne peut être qu'une approximation puisque comme on l'a vu, il existe une relation implicite entre les moments et q q ap. Cette première estimation ne la prend pas en compte et YY Li ay x {ail aY doit donc être corrigée. On modifie l'estimation de Keq(n) par exemple par la relation de récurrence suivante: Keq(n+1)=F(q q) , }x, 1'1y Sx;, ,x). Ici, les quantités Sxn, fjy fl sont calculées à l'aide du simulateur numérique: on y, résout le problème aux limites (P) en prenant comme tenseur de perméabilité l'estimation précédente Keq(n). On arrête les itérations une fois que les quantités f p a Y, } x a calculées par le simulateur avec l'estimation correspondante Keq(n) sont à e de la valeur 20 désirée, étant un critère d'arrêt fixé par l'utilisateur par rapport à une norme choisie à l'avance. L'organigramme est indiqué à la figure 1. Le quadruplet de mesure utilisé pour déterminer cômplètement le tenseur k peut également avoir été calculé à l'aide d'une "simulation fine" de l'expérience sur un milieu hétérogène dont on souhaite définir le tenseur de perméabilité équivalent par cette méthode. La perméabilité de ce milieu est connue en tout ou en partie, en se basant par exemple sur une carte de perméabilité fournie par le géologue. Durant cette simulation, on applique les lois de Darcy aux écoulements dans le milieu et on en déduit le quadruplet de valeurs applicable. Validation Par simulation, nous avons vérifié l'existence de forces transversales, mesurables et non négligeables dans certains cas, apparaissant lorsqu'un écoulement traverse un milieu confiné, due à la présence d'hétérogénéités internes ou à une anisotropie du milieu, et validé la possibilité d'identifier toutes les composantes du tenseur de perméabilité équivalent grâce à des mesures simples de forces sur échantillons obtenues au moyen d'un perméamètre. En plus, cette étude nous a permis d'établir l'erreur commise sur le comportement d'un milieu équivalent (défini par une méthode ou une autre) par rapport au comportement du milieu réel à homogénéiser. C'est sur cette base d'un meilleur ajustement entre les comportements respectifs que la méthode selon l'invention a été développée. Sur la fig.2B, on voit l'allure du champ de pression pour un milieu homogène anisotrope avec un rapport d'anisotropie de 10 et dont les axes principaux du tenseur diagonal (fig.2A) sont tournés de 45 par rapport à l'horizontale (0 = 45 ). Sur la fig.2C, on montre l'évolution du rapport de la différence cumulée de pression sur les bords imperméables (bords perpendiculaire aux champs de pression imposé) sur la différence de pression appliquée sur l'entrée-sortie du milieu (bords parallèles aux pressions imposées) pour 0 (0 Sur la fig.3A, on montre une carte de perméabilité log-normale (moyenne géométrique Kg=121, a2=2) générée avec le code de génération stochastique bien connu dit FFTMA. Localement, la perméabilité est isotrope (scalaire), en revanche sa fonction de corrélation est spatialement anisotrope, orientée à 45 . Sur la fig.3B,on a représenté l'évolution de la pression aux bords haut et bas (bords imperméables pour un écoulement confiné horizontal) pour le milieu réel et son équivalent homogène anisotrope déterminé par la méthode du perméamètre. Sur les fig.3C, 3D respectivement, on montre les cartographies du champ de pression pour le milieu réel (à gauche) et le milieu équivalent (à droite). Après inversion, on estime les composantes du tenseur de perméabilité du milieu équivalent valent: KXa=122.05; kyy=119.86; kxy=kyx=58.08. Les valeurs propres de ce tenseur valent: K',,,=179.05; K'yy=62.86. Le système d'axes principaux est bien tourné de 44.46 par rapport à la direction de l'écoulement imposé. Ceci correspond bien à la direction de l'anisotropie de départ. Sur la fig.4A, on montre le réseau de fractures dont la fraction volumique des fractures vaut 10%, la perméabilité de la matrice est de l'ordre de 1 Darcy alors que celle des fractures est de l'ordre de 100 Darcy's (à gauche). La fig.4B montre l'évolution de la pression aux bords haut et bas (bords imperméables pour un écoulement confiné horizontal) pour le milieu réel et son équivalent homogène anisotrope déterminé par la méthode du perméamètre. Sur les fig.4C, 4D respectivement, on montre les cartographies du champ de pression pour le milieu réel et son équivalent. Les composantes du tenseur de perméabilité du milieu équivalent valent: KXX=5.79 D; kyy=4.37 D; kXy=ky,=0.26 D. Les valeurs propres de ce tenseur valent: K'X,F5.83 D; K'yy=4.32 D et dont le système d'axes principaux sont tournés à 9.83 par rapport au système d'axe formé par la direction horizontale et celle verticale. REVENDICATIONS 1) Méthode pour déterminer, à partir des mesures connues de perméabilité, toutes les composantes d'un tenseur de perméabilité absolue (Keq,) d'un échantillon de milieu 5 poreux anisotrope, dans laquelle: on mesure ou bien on calcule à partir de deux (ou trois) écoulements selon Ox, Oy (et Oz) sous gradient de pression imposé les quantités les débits de liquide et les composantes des forces transversales générées dans l'échantillon du fait du passage du liquide, ces quantités étant ramenées à des gradients de pression unitaires; et - - on utilise un simulateur pour déterminer de proche en proche, à partir d'une solution numérique a priori d'un écoulement confiné dans un milieu homogène, les composantes du tenseur équivalent (k) par une méthode inverse quelconque tel que les dites quantités correspondant s'ajustent sensiblement aux quantités mesurées ou calculées. 2) Méthode selon la revendication 1, dans laquelle on mesure les dites quantités au moyen d'un perméamètre et de moyens de mesure de forces transversales appliquées à l'échantillon en réponse à des écoulements imposés au travers de l'échantillon. 3) Méthode selon la revendication 1., dans laquelle on détermine les dites quantités par simulation d'écoulements dans un milieu, à partir de données de perméabilité connues 20 ou estimées. 4) Méthode selon la revendication 2, dans laquelle on détermine les forces transversales appliquées à l'échantillon en réponse à des écoulements imposés au travers de l'échantillon, eh mesurant les variations de poids apparent de l'échantillon dans le perméamètre quand on alterne le sens de passage des écoulements de liquide au travers de l'échantillon.