# Approche multi-échelles des écoulements dans les milieux poreux fissurés

Julien Réthoré, René de Borst, Marie-Angèle Abellan

► To cite this version:

Julien Réthoré, René de Borst, Marie-Angèle Abellan. Approche multi-échelles des écoulements dans les milieux poreux fissurés. 8e Colloque national en calcul des structures, CSMA, May 2007, Giens, France. hal-01508940

HAL Id: hal-01508940

https://hal.science/hal-01508940v1

Submitted on 15 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Approche multi-échelles des écoulements dans les milieux poreux fissurés

Julien Réthoré* — René de Borst*,** — Marie-Angèle Abellan***

* Faculty of Aerospace Engineering, Delft University of Technology, Delft, Pays-Bas
{j.rethore,r.deborst}@lr.tudelft.nl

** Laboratoire de Mécanique des Contacts et des Solides, INSA Lyon, Lyon, France

*** LTDS ENISE, Saint-Etienne, France

RÉSUMÉ. On propose une approche multi-échelle pour la simulation numérique des écoulements dans les milieux poreux fissurés. A l'échelle d'une fissure, l'écoulement est supposé unidimensionnel. Le couplage est ensuite assuré avec l'écoulement dans le réseau poreux grâce aux termes d'interfaces apparaissant dans l'écriture d'une formulation faible en présence d'une interface de discontinuité. Pour prendre en compte de tels phénomènes indépendemment d'un maillage pré-existant, la méthode de Partition de l'Unité est ensuite utilisée afin d'enrichir les champs de déplacement et de pression.

ABSTRACT. A multi-scale approach of the fluid flow in fractured porous media is proposed. A weak formulation enables the coupling between the flow in the bulk and the 1D flow inside a crack. For a mesh independent description, the Partition of Unity method is used to enrich the displacement as well as the pressure interpolation.

MOTS-CLÉS : milieux poreux, écoulements, simulations multi-échelles, fissuration.

KEYWORDS: porous medium, fluid flow, multi-scale modelling, fracture.
# 1. Introduction

Depuis les travaux de Terzaghi (Terzaghi, 1943) et Biot (Biot, 1965), l'écoulement des fluides dans les milieux poreux déformables est un sujet d'étude en ingénierie. En effet, ce sujet est crucial dans la compréhension de nombreux systèmes. On mentionnera la géotechnique, l'ingénierie pétrolière comme les exemples d'applications les plus classiques. Mais, le fait que beaucoup de tissus humains sont en fait des milieux poreux, et doivent être étudiés comme tels, donnent de nouvelles perspectives à l'étude des écoulements au sein des milieux poreux.

De plus, la présence d'endommagement comme les fissures, les failles, les bandes de cisaillements modifient grandement les caractéristiques des écoulements. En effet, la physique qui les gouvernent à l'intérieure de telles discontinuités est très différente de celle des écoulements intersticiels dans un milieu poreux déformable. Ces différences affectent donc l'allure des écoulements et des déformations à proximité de la discontinuité.

Dans cette contribution, nous développons un modèle numérique général pour l'étude des écoulements de fluides dans les milieux poreux fracturés. En s'appuyant sur de précédents travaux dans lesquels nous avons élaboré une méthode numérique pour la propagation de rupture dans les milieux poreux saturés (Réthoré et al., 2006a), (de Borst et al., 2006), nous étendons cette étude pour inclure les écoulements à l'intérieur des cavités crées par la présence de fissures. Pour ce faire, on adopte une modélisation à deux échelles. Ensuite, les propriétés de partition de l'unité des fonctions de formes éléments finis (Babuska et al., 1997) sont exploitées afin d'enrichir l'interpolation des champs de déplacement et de pression au voisinage des discontinuités.

# 2. Écoulement au sein d'une cavité

A l'échelle fine, l'écoulement à l'intérieur de la cavité de la fissure ouvrante est modélisé en utilisant l'équation de conservation suivante :

$$ \mu \nabla \cdot \nabla \mathbf{v} = \nabla p \quad [1] $$

où $\mu$ est la viscosité du fluide, $\mathbf{v}$ sa vitesse et $p$ sa pression. Comme l'ouverture de cette cavité peut, dans un premier temps, être négligé par rapport à sa longueur, on obtient les simplifications suivantes :

$$ 0 = \frac{\partial p}{\partial y} \qquad ; \qquad \mu \frac{\partial^2 v}{\partial y^2} = \frac{\partial p}{\partial x} \quad [2] $$

dans la direction normale $\mathbf{n}_{\Gamma_d}$ et dans la direction tangentielle à la fissure $\mathbf{t}_{\Gamma_d}$. $x, y$ repère la position d'un point matériel dans la cavité. $v$ est la vitesse du fluide dans la
direction $\mathbf{t}_{\Gamma_d}$. La pression $p$ est donc constante à travers la cavité. En intégrant deux fois l'équation [2b] on obtient le profil de vitesse suivant :

$$v(y) = \frac{1}{2\mu} \frac{\partial p}{\partial x} (y^2 - h^2) + v_f \quad [3]$$

où $v = v_f = \mathbf{v}_f \cdot \mathbf{t}_{\Gamma_d}$ a été utilisé comme condition à la paroi en $y = -h$ et $y = h$.

De façon analogue, et sous les mêmes hypothèses, on peut obtenir l'équation de conservation de la masse :

$$\frac{\partial v}{\partial x} + \frac{\partial w}{\partial y} = 0 \quad [4]$$

où $w = \mathbf{v}_m \cdot \mathbf{n}_{\Gamma_d}$. Pour relier cet écoulement à l'échelle "micro" (échelle de la fissure), avec l'écoulement dans le réseau poreux environnant, on intègre l'équation de conservation de la masse à travers la cavité :

$$\int_{-h}^{h} \frac{\partial v}{\partial x} dy + \int_{-h}^{h} \frac{\partial w}{\partial y} dy = 0 \quad [5]$$

Après quelques manipulations, on obtient :

$$[w] = w(h) - w(-h) = \frac{2}{3\mu} \frac{\partial}{\partial x} \left( \frac{\partial p}{\partial x} h^3 \right) - 2h \frac{\partial v_f}{\partial x} \quad [6]$$

### 3. Écoulement principal

On choisit de décrire la physique dans le réseau poreux par le champ de déplacement de la phase solide et le champ de pression dans le fluide. En se plaçant dans le cadre de l'hypothèse des petites perturbations et des petites variations de concentration, on obtient finalement en négligeant les termes convectifs :

$$\mathbf{D}^{\text{tan}} \nabla^s \mathbf{u}_s - \nabla p = \mathbf{0} \quad [7]$$

pour la conservation de la quantité de mouvement ($\boldsymbol{\sigma}$ est le tenseur des contraintes de Cauchy, $\mathbf{D}^{\text{tan}}$ la matrice de comportement élastique de la phase solide) et

$$\alpha \nabla \cdot \mathbf{v}_s - k_f \nabla \nabla p + \frac{1}{Q} \frac{\partial p}{\partial t} = 0 \quad [8]$$

pour la conservation de la masse ($k_f$ est le coefficient de perméabilité introduit par l'adoption de la loi de Darcy $n_f(\mathbf{v}_f - \mathbf{v}_s) = -k_f \nabla p$ et $Q$ le module de Biot).

On élabore ensuite une formulation faible en présence d'une discontinuité $\Gamma_d$. Il en découle les équations suivantes :

$$\int_{\Omega} (\nabla \cdot \boldsymbol{\eta}) \cdot \boldsymbol{\sigma} \, d\Omega + \int_{\Gamma_d} [\boldsymbol{\eta} \cdot \boldsymbol{\sigma}] \cdot \mathbf{n}_{\Gamma_d} \, d\Omega = \int_{\Gamma} \boldsymbol{\eta} \cdot \mathbf{t}_p \, d\Omega \quad [9]$$
et

$$
\begin{equation}
\begin{split}
& - \int_{\Omega} \alpha \zeta \nabla \cdot \mathbf{v}_s \, d\Omega + \int_{\Omega} k_f \nabla \zeta \cdot \nabla p \, d\Omega - \int_{\Omega} \zeta Q^{-1} \dot{p} \, d\Omega \\
& \qquad + \int_{\Gamma_d} \mathbf{n}_{\Gamma_d} \cdot [\zeta n_f (\mathbf{v}_f - \mathbf{v})] \, d\Gamma = \int_{\Gamma} \zeta \mathbf{n}_{\Gamma} \cdot \mathbf{q}_p \, d\Gamma
\end{split}
\tag{10}
\end{equation}
$$

$\eta$ et $\zeta$ sont les fonctions tests pour le déplacement solide et la pression fluide. Pour l'équation [9], les conditions aux limites suivantes ont été utilisées sur des parties complémentaires $\partial\Omega_t$ et $\partial\Omega_u$ de la frontière :

$$
\mathbf{n}_{\Gamma} \cdot \boldsymbol{\sigma} = \mathbf{t}_{p} \quad , \quad \mathbf{u} = \mathbf{u}_{p} \tag{11}
$$

$t_p$ étant les efforts extérieurs imposés et $\mathbf{u}_p$ les déplacement imposés. Pour l'équation [10], on applique sur des parties complémentaires $\partial\Omega_p$ et $\partial\Omega_p$ de la frontière les conditions aux limites suivantes :

$$
n_f(\mathbf{v}_f - \mathbf{v}_s) \cdot \mathbf{n}_{\Gamma} = \mathbf{q}_p \cdot \mathbf{n}_{\Gamma} \quad , \quad p = p_p \tag{12}
$$

$\mathbf{q}_p \cdot \mathbf{n}_\Gamma$ et $p_p$ étant respectivement le flux et la pression imposés.

4. Couplage entre les écoulements

La présence d'un terme d'intégrale sur l'interface de discontinuité permet de cou-
pler l'écoulement dans la cavité à l'écoulement dans le réseau poreux. Le couplage
mécanique provient des efforts appliqués par le fluide dans la cavité sur les faces de
celle-ci :

$$
\boldsymbol{\sigma} \cdot \mathbf{n}_{\Gamma_d} = -p \mathbf{n}_{\Gamma_d} \quad \text{on} \quad \Gamma_d \tag{13}
$$

Le couplage par transfert de masse se simplifie alors (les champs $p$ et donc $\zeta$ étant
continus à l'interface) et le terme $n_f \mathbf{n}_{\Gamma_d} \cdot [[\mathbf{v}_f - \mathbf{v}_s]]$ est obtenu à partir des équations
moyennées de l'écoulement dans la cavité :

$$
\mathbf{n}_{\Gamma_d} \cdot [\mathbf{v}_f - \mathbf{v}_s] = n_f \left( \frac{2}{3\mu} \frac{\partial}{\partial x} \left( \frac{\partial p}{\partial x} h^3 \right) - 2h \frac{\partial v_f}{\partial x} - 2 \frac{\partial h}{\partial t} \right) \quad [14]
$$

avec $h = [[\mathbf{u}_s]] \cdot \mathbf{n}_{\Gamma_d}$ l'ouverture de fissure due au terme de transfert de quantité de mouvement et $v_f$ tel que $v_f = (\mathbf{v}_s - \frac{k_f}{n_f} \nabla p) \cdot \mathbf{t}_{\Gamma_d}$

5. Discrétisation

D'après les résultats obtenus précédemment, les champs de déplacement et de
pression doivent pouvoir décrire une certaine “cinématique”. En effet, le déplacement
de la phase solide est discontinu à travers l'interface et il présente une forme parti-
culière en $\sqrt{r}$ au voisinage de la pointe de la fissure. Pour le champ de pression, le
couplage introduit par l'équation [14] impose une discontinuité de la composante de
la vitesse relative du fluide, soit en utilisant la loi de Darcy, une discontinuité de la
composante normale du gradient de la pression. On choisit d'utiliser la méthode de
partition de l'unité (Babuska et al., 1997) pour enrichir une interpolation éléments fi-
nis classique de manière à satisfaire ces conditions. Le déplacement est décomposé
comme suit :

$$
u_s = \sum_{i \in \mathcal{N}} N_i \bar{u}_i + \sum_{i \in \mathcal{N}_{\text{cut}}} N_i \mathcal{H}_{\Gamma_d} \hat{u}_i + \sum_{i \in \mathcal{N}_{\text{tip}}} \sum_{j \in [1,4]} N_i \Psi_j \tilde{u}_{ij} \quad [15]
$$

où $N_i$ sont les fonctions de forme éléments finis classiques supportées par l'ensemble
de nœuds $\mathcal{N}$ discrétisant le domaine $\Omega$. Les nœuds du sous-ensemble $\mathcal{N}_{cut}$ ont leur
support complètement coupé par la fissure. Ils portent des degrés de liberté enrichis $\hat{u}_i$
associés à la fonction discontinue $\mathcal{H}_{\Gamma_d}$ définie par :

$$
\mathcal{H}_{\Gamma_d}(\mathbf{x}) = \frac{\mathbf{x} \cdot \mathbf{n}_{\Gamma_d}}{\|\mathbf{x} \cdot \mathbf{n}_{\Gamma_d}\|} \quad [16]
$$

$\mathcal{N}_{tip}$ contient les nœuds dont le support contient un front de la fissure. Quatre fonctions $\Psi_j$ et leur degrés de liberté correspondant $\tilde{u}_{ij}$ sont supportés par les nœuds de cet ensemble. Les fonctions singulières sont :

$$
\{\Psi_j\}_{j \in [1,4]} = \left\{ \sqrt{r} \cos \frac{\theta}{2}, \sqrt{r} \cos \frac{\theta}{2} \sin \theta, \sqrt{r} \sin \frac{\theta}{2}, \sqrt{r} \sin \frac{\theta}{2} \sin \theta \right\} \quad [17]
$$

où $r, \theta$ sont les coordonnées polaires dans le repère lié au front de la fissure (voir (Moës et al., 1999)).

Pour la pression, on utilise la décomposition suivante :

$$
p = \sum_{i \in \mathcal{N}} H_i \bar{p}_i + \sum_{i \in \mathcal{N}_{\text{pres}}} H_i \mathcal{D}_{\Gamma_d} \hat{p}_i \quad [18]
$$

où $H_i$ sont des fonctions de forme standards et $\mathcal{D}_{\Gamma_d}$ la fonction distance à l'interface. Cette fonction dont la composante normale du gradient est discontinue, est utilisée comme enrichissement pour tous les nœuds dont le support est coupé par la fissure ($\mathcal{N}_{pres}$). Le système non-linéaire obtenu est ensuite résolu par un algorithme itératif de Newton modifié.

6. Résultats

La figure 1 présente les résultats obtenus pour la simulation de l’écoulement forcé
à travers un massif avec une fissure à 30°. On observe effectivement la discontinuité
de la composante normale du gradient de la pression. Cette discontinuité est dictée
par la modélisation à deux échelles adoptée plus haut. Pour plus de résultats et de
détails concernant l’implémentation de cette technique, le lecteur intéressé est prié
de se référer à (Réthoré et al., 2006b). Une application à la propagation de fissures
cohésives dans les milieux poreux est aussi proposée dans (Réthoré et al., 2006c).
<figure><img src="image_2.png" /><figcaption>**Figure 1.** *Composante normale du gradient de pression pour une fissure inclinée à 30°.*</figcaption></figure>

**7. Bibliographie**

Babuska I., Melenk J., « The Partition of unity method », *International Journal for Numerical Methods in Engineering*, vol. 40, p. 727-758, 1997.

Biot M., *Mechanics of Incremental Deformations*, John Wiley & Sons, Chichester, 1965.

de Borst R., Réthoré J., Abellan M., « A numerical approach for arbitrary cracks in a fluid-saturated medium », *Archive of Applied Mechanics*, vol. 75, p. 595-606, 2006.

Jouanna P., Abellan M., *Modern Issues in Non-Saturated Soils*, Springer-Verlag, chapter Generalized approach to heterogeneous media, 1995.

Moës N., Dolbow J., Belytschko T., « A finite element method for crack growth without remeshing », *International Journal for Numerical Methods in Engineering*, vol. 46, n° 1, p. 133-150, 1999.

Réthoré J., de Borst R., Abellan M., « A discrete model for the dynamic propagation of shear bands in a fluid-saturated medium », *International Journal for Numerical and Analytical Methods in Geomechanics*, 2006a.

Réthoré J., de Borst R., Abellan M., « A two-scale approach for fluid flow in fractured porous media », *International Journal for Numerical Methods in Engineering*, 2006b.

Réthoré J., de Borst R., Abellan M., « A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks », *Computational Mechanics*, 2006c.

Terzaghi K., *Theoretical Soil Mechanics*, John Wiley & Sons, New York, 1943.