Définition de la perméabilité équivalente des massifs fracturés par des méthodes d'homogénéisation

Ahmad Pouya; Alexis Courtois

doi:10.1016/S1631-0713(02)01839-4

Résumés

Français
Anglais

Dans les méthodes utilisées pour déterminer le comportement hydraulique à grande échelle des massifs fracturés, passant par la simulation de l'écoulement dans un réseau de fracture, le gradient moyen de pression et le flux moyen dans le réseau ne sont pas calculés de manière rigoureuse. On donne une méthode de détermination de ces quantités moyennes dans un milieu de perméabilité hétérogène à partir des valeurs de pression et de flux sur le contour du domaine étudié. On définit ensuite un tenseur de conductivité hydraulique moyenne du domaine en partant des conditions de pression, qui varient linéairement sur le contour. On montre que ce tenseur est symétrique et défini positif. Un exemple d'application à un modèle de milieu fracturé est donné.

In the methods used for determination of the hydraulic behaviour at large scale of fractured rock masses, based on the simulation of flow in a fracture network, the mean flux and the mean pressure gradient in the network are not rigorously determined. A method is given for deriving these quantities, in a heterogeneous permeable block, from pressure and flux values on the boundary of the block. A block conductivity tensor is then defined, based on the condition of linear variation of the pressure on the boundary of the block. It is shown that this conductivity tensor is symmetric and positive-definite. An example of application to a model of fractured medium is given.

1 Introduction

Une méthode courante de détermination des propriétés hydrauliques à grande échelle des massifs fracturés consiste à simuler l'écoulement dans un réseau de fractures et à essayer d'établir des relations entre le gradient moyen de pression et le flux moyen dans ce réseau [1,5,10]. Cette méthode, dite directe [2], s'inscrit dans le cadre plus général des méthodes dites non locales de détermination de la perméabilité effective des milieux hétérogènes [3,6,7,9]. Dans ces méthodes, le flux moyen $\vec{Q}$ et le gradient moyen de pression $\vec{G}$ sont définis, explicitement [8] ou implicitement, par les relations (3) données ci-après. Mais une analyse des méthodes de calcul adoptées dans ces travaux révèle une certaine difficulté à calculer ces grandeurs moyennes. La difficulté réside dans le fait que les moyennes sont définies sur la distribution volumique de ces grandeurs, tandis que, dans les simulations numériques, les valeurs auxquelles on accède plus facilement, et parfois exclusivement, sont les valeurs de ces grandeurs aux frontières du volume considéré. Pour illustrer ces difficultés, considérons l'exemple du travail de Long et al. [5], portant sur le calcul de la perméabilité équivalente des massifs fracturés. Ces auteurs reconnaissent la difficulté d'accéder au flux moyen dans un réseau de fractures, du fait que le flux local est contrôlé par la direction des fractures. Pour contourner cette difficulté, ils proposent une méthode passant par le calcul de la perméabilité dans la direction du gradient. Cette méthode exige de considérer un domaine carré (ou rectangulaire), soumis à un gradient de pression parallèle à un de ses côtés (conditions de pression linéaire au contour [7]) et de calculer le débit sur les côtés orthogonaux à ce gradient. Mais une difficulté subsiste encore : le débit entrant d'un côté n'est pas égal à celui sortant du côté opposé, ce qui pose, comme le soulignent les auteurs, un problème pour définir la perméabilité dans la direction du gradient. Cacas et al. [1] apportent une légère modification à cette méthode pour lever cette difficulté, en imposant une condition de flux nul aux côtés parallèles à la direction du gradient (conditions de perméamètre [7]). Mais, dans ce cas, la pression sur ces côtés n'est plus contrôlée : cela pose le problème du calcul du gradient moyen de pression dans le domaine. Long et al. [5] choisissent un des côtés orthogonaux à la direction du gradient pour calculer le débit et la perméabilité dans cette direction, notée K_g. Le critère qu'ils adoptent pour juger si le réseau peut être assimilé à un milieu poreux est que le diagramme donnant $1 / \sqrt{K_{g}}$ en fonction de la direction du gradient s'approche d'une ellipse, symétrique par rapport à l'origine. Il faut noter que cette symétrie signifie uniquement l'égalité des débits entrant et sortant par les côtés opposés. Quant à la symétrie du tenseur de conductivité hydraulique du réseau, autre condition nécessaire pour pouvoir assimiler le réseau à un milieu poreux, cette méthode ne permet pas d'en juger. En effet, même pour un tenseur $K$ non symétrique, en notant $K_{g} = K_{ij} n_{i} n_{j}$ , où $\vec{n}$ est un vecteur unitaire, et en traçant le digramme de $1 / \sqrt{K_{g}}$ en fonction de la direction de $\vec{n}$ , on obtient une ellipse.

À la difficulté du calcul des valeurs moyennes s'ajoute celle de la définition même du tenseur de perméabilité équivalente, qui doit donner $\vec{Q}$ en fonction de $\vec{G}$ . Cette définition paraı̂t ne pas être unique et dépendre des conditions aux limites [7]. En effet, différentes conditions aux limites peuvent conduire à une même valeur de $\vec{G}$ , mais à différentes valeurs de $\vec{Q}$ , ce qui pose un problème pour la définition du tenseur de perméabilité équivalente. Mais, tant qu'une définition rigoureuse conduisant à une valeur unique de ce tenseur n'est pas adoptée, on ne peut, ni aller plus loin, ni discuter de propriétés telles que la pertinence ou la symétrie de ce tenseur.

Nous allons d'abord proposer dans ce travail une méthode permettant de calculer correctement les moyennes volumiques du flux et du gradient de pression à partir des données sur la frontière. Ensuite, en faisant le choix d'une famille particulière de conditions aux limites, nous donnons une définition rigoureuse d'un tenseur de perméabilité équivalente, dont nous étudierons les propriétés. Notre définition et notre démarche s'inspirent des méthodes générales d'homogénéisation du comportement mécanique des matériaux hétérogènes [11].

2 Position du problème

Considérons un corps perméable hétérogène, occupant un domaine $Ω$ dans l'espace. Ce corps est le siège d'un écoulement de fluide sous l'effet de pressions ou de flux imposés sur son contour $\partial Ω$ . On suppose qu'il est constitué de matériaux de perméabilités différentes, caractérisés chacun par un tenseur de perméabilité noté $k$ , symétrique et défini positif. En tout point de $Ω$ , caractérisé par son vecteur position $\vec{x}$ , l'écoulement a lieu suivant la loi de Darcy :

\forall \vec{x} \in Ω; \vec{q} (\vec{x}) = - k (\vec{x}) \cdot \vec{\nabla} p (\vec{x})

(1)

où

\vec{q} (\vec{x})

est le flux au point

\vec{x}

\vec{\nabla} p (\vec{x})

et le gradient de pression en ce point. L'équation de conservation de la masse s'écrit :

\forall \vec{x} \in Ω; div \vec{q} (\vec{x}) = 0

(2)

Les champs p et $\vec{q}$ , solutions du problème d'écoulement, doivent vérifier ces deux équations et les conditions aux limites de flux ou de pression sur $\partial Ω$ . On définit le flux moyen et le gradient moyen de pression dans le corps par les relations suivantes, dans lesquelles V est le volume du domaine $Ω$ :

\vec{Q} = \frac{1}{V} \int_{Ω} \vec{q} (\vec{x}) d v, \vec{G} = \frac{1}{V} \int_{Ω} \vec{\nabla} p (\vec{x}) d v

(3)

3 Expression des moyennes volumiques en fonction des valeurs aux frontières

Pour une fonction f quelconque, on peut écrire l'identité mathématique suivante (formule de Green), dans laquelle $\vec{n}$ représente le vecteur unitaire sortant et dS l'élément de surface sur $\partial Ω$ :

\int_{Ω} \partial_{i} f (\vec{x}) d v = \int_{\partial Ω} f (\vec{x}) n_{i} (\vec{x}) d S

(4)

En prenant pour f la pression, on trouve que la moyenne du gradient de pression peut s'écrire sous la forme suivante :

\vec{G} = \frac{1}{V} \int_{\partial Ω} p (\vec{x}) \vec{n} (\vec{x}) d S

(5)

En utilisant ∂_ix_j=δ_ij (coordonnées cartésiennes) et la conservation de la masse, ∂_iq_i=0, on peut écrire : $\partial_{i} (x_{j} q_{i}) = (\partial_{i} x_{j}) q_{i} + x_{j} \partial_{i} q_{i} = δ_{ij} q_{j} = q_{i}$ . En remplaçant q_i dans (3) par ∂_i(x_jq_i) et en utilisant (4), on trouve l'expression suivante pour le flux moyen :

\vec{Q} = \frac{1}{V} \int_{\partial Ω} (\vec{q} \cdot \vec{n}) \vec{x} d S

(6)

Les formules (5) et (6) permettent de calculer les moyennes volumiques du gradient de pression et du flux à l'aide des pressions et du flux au contour pour une forme du domaine et des conditions aux limites quelconques.

En général, dans les méthodes directes ou non locales, essentiellement deux types de conditions aux limites sont considérés. Les conditions de perméamètre supposent une forme rectangulaire (dans le cas bidimensionnel) pour le domaine, dont deux côtés parallèles sont soumis à des pressions constantes et les deux autres, à un flux sortant nul $(\vec{q} \cdot \vec{n} = 0)$ . Les conditions de pression linéaire au contour sont définies par :

\forall \vec{x} \in \partial Ω p (\vec{x}) = \vec{A} \cdot \vec{x} + P

(7)

où

\vec{A}

est un vecteur constant et P un scalaire constant.

Contrairement à ce qui est supposé dans les travaux utilisant les conditions de perméamètre, l'équation (5) montre que, sous ces conditions, la moyenne du gradient de pression ne peut être déduite des seules valeurs de la pression sur les côtés à pression imposée. La formule (6) permet en revanche de montrer que, sous ces conditions, la moyenne volumique du flux est bien égale à la moyenne surfacique du flux sur les côtés à pression imposée, comme cela est souvent admis dans ces travaux.

Considérons maintenant le cas des conditions de pression linéaire au contour. La formule (6) montre que, dans ce cas où le flux n'est (imposé) nul sur aucun côté, le flux moyen ne peut se déduire de la moyenne surfacique sur un seul des côtés. Ceci est à l'origine de beaucoup de difficultés dans les travaux [5] utilisant ces conditions.

En prenant dans (4) pour f une constante non nulle et en prenant ensuite pour f la composante x_j du vecteur position $\vec{x}$ , sachant que ∂_ix_j=δ_ij, on établit successivement les deux identités suivantes :

\int_{\partial Ω} \vec{n} d S = 0, \int_{\partial Ω} x_{i} n_{j} d S = V δ_{ij}

(8)

En reportant l'expression (7) de $p (\vec{x})$ dans (5) et en utilisant les identités (8), on trouve :

\vec{G} = \vec{A}

(9)

Ce résultat, admis en général dans les travaux utilisant des conditions de pression linéaire au contour sur des domaines rectangulaires avec $\vec{A}$ parallèle à un des côtés, se trouve ici établi pour un domaine de forme quelconque.

4 Tenseur de conductivité hydraulique

Pour un domaine $Ω$ fixé, avec une distribution fixée de $k (\vec{x})$ dans $Ω$ , le champ de pression $p (\vec{x})$ solution du problème d'écoulement avec les conditions aux limites de pression linéaire au contour avec les données $(\vec{A}, P)$ est défini par l'équation (7) et l'équation suivante (10) :

\forall \vec{x} \in Ω; div [k (\vec{x}) \cdot \vec{\nabla} p (\vec{x})] = 0

(10)

En remarquant que, si le champ $p (\vec{x})$ est solution de ce problème pour les données $(\vec{A}, P)$ et si c est une constante quelconque, alors le champ $p (\vec{x}) + c$ est solution du même problème pour les données $(\vec{A}, P + c)$ , et que $\vec{\nabla} (p (\vec{x}) + c) = \vec{\nabla} p (\vec{x})$ , on déduit que le champ de gradient de pression, et donc le champ de flux donné par $\vec{q} (\vec{x}) = k (\vec{x}) \cdot \vec{\nabla} p (\vec{x})$ ainsi que sa moyenne volumique $\vec{Q}$ sont indépendants de la constante c, et donc de la constante P, et ne dépendent que de $\vec{A}$ (et de données intrinsèques, qui sont la géométrie de $Ω$ et la distribution $k (\vec{x})$ ). Considérons maintenant deux vecteurs ${\vec{A}}_{1}$ et ${\vec{A}}_{2}$ et notons respectivement $p_{1} (\vec{x})$ et $p_{2} (\vec{x})$ , ${\vec{q}}_{1} (\vec{x})$ et ${\vec{q}}_{2} (\vec{x})$ et ${\vec{Q}}_{1}$ et ${\vec{Q}}_{2}$ les champs de pression et de flux et les flux moyens associés aux conditions aux limites $({\vec{A}}_{1}, 0)$ et $({\vec{A}}_{2}, 0)$ . On vérifie aisément que le champ $p (\vec{x}) = λ_{1} p_{1} (\vec{x}) + λ_{2} p_{2} (\vec{x})$ , où λ₁ et λ₂ sont deux constantes quelconques, est solution du problème d'écoulement avec les conditions aux limites $(\vec{A}, 0)$ , où $\vec{A} = λ_{1} {\vec{A}}_{1} + λ_{2} {\vec{A}}_{2}$ . On en déduit que le flux $\vec{Q}$ associé à $λ_{1} {\vec{A}}_{1} + λ_{2} {\vec{A}}_{2}$ est égal à $λ_{1} {\vec{Q}}_{1} + λ_{2} {\vec{Q}}_{2}$ , ce qui veut dire que $\vec{Q}$ est une fonction linéaire de $\vec{A}$ . Ce résultat s'écrit sous la forme :

\vec{Q} = - K_{Ω} \cdot \vec{A}

(11)

Le tenseur $K_{Ω}$ est une caractéristique intrinsèque de $Ω$ , ne dépendant que de sa géométrie et de la distribution $k (\vec{x})$ . Considérons maintenant deux champs de pression $p (\vec{x})$ et $p' (\vec{x})$ , résultant des conditions aux limites de pression linéaire au contour avec les données respectivement $(\vec{A}, P)$ et $(\vec{A}', P')$ , et notons $\vec{q}$ et $\vec{q}'$ les champs de flux correspondants. En tout point de $Ω$ , la relation $\vec{q}' = - k \cdot \vec{\nabla} p'$ et la conservation de la masse, ∂_iq′_i=0, permettent d'écrire :

- \vec{\nabla} p \cdot k \cdot \vec{\nabla} p' = \vec{\nabla} p \cdot \vec{q}' = (\partial_{i} p) q'_{i} = \partial_{i} (p q'_{i}) - p \partial_{i} q'_{i} = \partial_{i} (p q'_{i})

En intégrant sur $Ω$ , et en appliquant (4) au second membre, on trouve :

- \frac{1}{V} \int_{Ω} \vec{\nabla} p \cdot k \cdot \vec{\nabla} p' d v = \frac{1}{V} \int_{\partial Ω} p \vec{q}' \cdot \vec{n} d S

Or, sur $\partial Ω$ , on a $p (\vec{x}) = \vec{A} \cdot \vec{x} + P$ . En reportant cette expression dans le second membre de la relation ci-dessus, on trouve :

- \frac{1}{V} \int_{Ω} \vec{\nabla} p \cdot k \cdot \vec{\nabla} p' d v = \frac{1}{V} \int_{\partial Ω} (\vec{A} \cdot \vec{x} + P) \vec{q}' \cdot \vec{n} d S = \vec{A} \cdot (\frac{1}{V} \int_{\partial Ω} (\vec{q}' \cdot \vec{n}) \vec{x} d S) + \frac{P}{V} \int_{\partial Ω} \vec{q}' \cdot \vec{n} d S

Dans le dernier membre de ces égalités, $\int_{\partial Ω} \vec{q}' \cdot \vec{n} d S = 0$ , du fait de la conservation de la masse, et le terme entre parenthèse est, d'après (6) et (11), égal à $- K_{Ω} \vec{A}'$ . On peut donc écrire :

\frac{1}{V} \int_{Ω} \vec{\nabla} p \cdot k \cdot \vec{\nabla} p' d v = \vec{A} \cdot K_{Ω} \cdot \vec{A}'

(12)

Or, du fait de la symétrie de $k$ , on peut changer le rôle de p et p′ au premier membre de cette égalité. Il en donc de même en ce qui concerne $\vec{A}$ et $\vec{A}'$ , au second membre. On en déduit que $K_{Ω}$ est symétrique.

Ce résultat contredit celui de [5] rapporté dans [7], suivant lequel « les conditions aux limites uniformes (variant linéairement sur le contour) ne produisent pas de perméabilité symétrique ». La différence résulte du fait que, comme nous l'avons expliqué ci-dessus, le flux moyen n'est pas correctement calculé dans [5].

Dans le cas particulier où $\vec{A}' = \vec{A}$ , la relation (12) devient :

\frac{1}{V} \int_{Ω} \vec{\nabla} p \cdot k \cdot \vec{\nabla} p d v = \vec{A} \cdot K_{Ω} \cdot \vec{A}

(13)

Comme $k$ est positif, le premier et donc le second membres de (13) sont toujours positifs, ce qui montre que $K_{Ω}$ est positif. Si un vecteur $\vec{A}$ annule le second membre, le $\vec{\nabla} p$ correspondant doit être nul au premier membre (car $k$ est défini), et donc sa moyenne $\vec{G}$ est également nulle. On déduit de alors (9) que $\vec{A}$ est nulle, ce qui montre que le tenseur $K_{Ω}$ est défini.

En reportant (9) dans (11), on trouve :

\vec{Q} = - K_{Ω} \cdot \vec{G}

(14)

Le tenseur $K_{Ω}$ reliant le flux moyen au gradient moyen de pression dans le domaine $Ω$ (sous conditions de pression linéaire au contour) sera appelé le tenseur de conductivité hydraulique moyenne de ce domaine. Nous venons de montrer qu'il est symétrique et défini positif. Par ailleurs, en reportant (9) dans (13), on trouve :

\frac{1}{V} \int_{Ω} \vec{\nabla} p \cdot k \cdot \vec{\nabla} p d v = \vec{G} \cdot K_{Ω} \cdot \vec{G}

(15)

Au premier membre de (15), nous trouvons l'expression de la dissipation donnée par Indelman et Dagan [4]. Cette relation montre que la dissipation macroscopique calculée par $K_{Ω}$ est égale à la dissipation moyenne réelle dans le milieu (sous conditions de pression linéaire au contour).

Supposons maintenant que $Ω$ soit un sous-domaine d'un milieu hétérogène infini et qu'il soit caractérisé par une longueur caractéristique D (par exemple le diamètre moyen). Supposons que, pour le milieu considéré, $K_{Ω}$ tende vers une limite $K$ quand D tend vers l'infini. Si la limite $K$ existe, elle est symétrique et positive. Dans ce cas, le comportement à grande échelle du milieu considéré peut être assimilé à un milieu poreux de perméabilité $K$ .

5 Application aux milieux fracturés

Dans le cas d'un milieu fracturé bi-dimensionnel, la formule (6) s'écrit sous la forme discrétisée :

\vec{Q} = \frac{1}{V} \sum_{κ} q^{(κ)} {\vec{x}}^{(κ)}

(16)

où

{\vec{x}}^{(κ)}

est le vecteur position d'un point d'intersection des fractures avec le contour du domaine et q^(κ), le flux sortant de la fracture en ce point. En calculant ce vecteur pour deux directions différentes de

\vec{A}

, on construit

K_{Ω}

. Cette méthode apporte une grande simplification par rapport à celles utilisées dans [5] et [1], qui passent par la construction d'un diagramme de perméabilité directionnelle. Pour illustrer ces résultats, nous avons considéré un carré ABCD dans un milieu fracturé (Fig. 1a). Le diagramme donnant

1 / \sqrt{K_{g}}

, calculé pour ce domaine par la méthode [5], c'est-à-dire par des rotations successives du carré par rapport au réseau de fractures (qui couvre un domaine plus grand que ABCD), est donné sur la Fig. 1b. Nous avons calculé un autre diagramme en faisant les mêmes rotations, mais en calculant chaque fois le flux

\vec{Q}

par la formule (16) et en la projetant sur la direction du gradient pour avoir la perméabilité directionnelle. On voit que les points obtenus de cette manière s'alignent presque sur une ellipse. La grande dispersion des points du premier diagramme est due au fait que le flux calculé par la méthode [5] ne représente pas correctement le flux moyen dans le domaine. Nous avons calculé également les

K_{Ω}

du carré ABCD et du domaine polygonal EFGHIJKLMNL par notre méthode. Pour comparer les résultats, nous avons représenté graphiquement ces deux tenseurs par les ellipses correspondantes. Ces ellipses et le deuxième diagramme se superposent presque, ce qui indique qu'à l'échelle d'hétérogénéité considérée (l'échelle des fractures), les domaines étudiés sont de taille suffisamment grande par rapport à la taille du volume élémentaire représentatif et que le

K_{Ω}

ainsi obtenu est très proche de la perméabilité effective du milieu.

Figure 1
a. Domaines étudiés. b. Diagramme obtenu par la méthode [5], en partant du carré ABCD : 1, diagramme obtenu par la même méthode, mais en utilisant l'expression (16) pour le flux moyen ; 2, ellipses représentant les tenseurs $K_{Ω}$ du domaine ABCD (3) et EFGHIJKLMN (4).
a. Studied domains. b. Diagram obtained by the method [5], beginning on the square ABCD: 1, diagram obtained by the same method but using the expression (16) for the mean flux; 2, ellipses representing $K_{Ω}$ tensors for the domains ABCD (3) and EFGHIJKLMN (4).

6 Conclusions

La perméabilité équivalente d'un domaine hétérogène et de forme quelconque peut être étudiée par une méthode non locale en appliquant sur la frontière de ce domaine des conditions de pression linéaire au contour, définies par la relation (7). Le gradient moyen de pression dans le domaine sera alors exactement égal à $\vec{A}$ . Le flux moyen peut être plus facilement calculé à partir des valeurs de flux sur la frontière par la formule (6). La relation entre ce flux moyen et $\vec{A}$ est linéaire et définit un tenseur de conductivité hydraulique moyenne pour le domaine. Ce tenseur peut être identifié en calculant le flux moyen pour deux ou trois (suivant la dimension du problème) directions différentes de $\vec{A}$ . Nous avons montré que ce tenseur est symétrique et défini positif si, en tout point du domaine, le tenseur local de perméabilité est symétrique et défini positif.