# Modélisation hydro-mécanique de propagation de fissures dans un milieu poreux déformable: Application au stockage géologique de CO2

Walid Saber Cherif, Frédéric Dubois, André Chrysochoos, Laurent Cangémi, Nadège Brusselle

► To cite this version:

Walid Saber Cherif, Frédéric Dubois, André Chrysochoos, Laurent Cangémi, Nadège Brusselle. Modélisation hydro-mécanique de propagation de fissures dans un milieu poreux déformable: Application au stockage géologique de CO2. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717105

HAL Id: hal-01717105

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

Submitted on 25 Feb 2018

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
# Modélisation hydro-mécanique de propagation de fissures dans un milieu poreux déformable : Application au stockage géologique de CO₂

W. Saber Cherif<fnref n="1" />,<fnref n="2" />, F. Dubois<fnref n="1" />,<fnref n="3" />, A. Chrysochoos<fnref n="1" />,<fnref n="3" />, L. Cangémi<fnref n="2" />, N. Brusselle<fnref n="2" />

**Résumé** — Un modèle numérique est développé pour l'analyse d'un écoulement multiphasique dans une formation géologique profonde poreuse, déformable et progressivement endommagée par l'injection de CO₂. L'écoulement multiphasique dans le milieu, et à travers les fissures, est décrit par une loi de Darcy généralisée. Les équations gouvernant le couplage avec le squelette solide, incluant les fissures, sont exprimées dans le cadre de la théorie généralisée de Biot. L'implémentation numérique repose sur la méthode des éléments finis pour discrétiser la forme faible des équations gouvernantes, et sur l'adaptation d'un modèle d'érosion de maillage pour simuler l'apparition de fissures.

**Mots clés** — érosion de maillage, milieu poreux insaturé, écoulement multiphasique, couplage.

## 1 Introduction

Une des solutions visant à atténuer le changement climatique, même si des avis contraires subsistés, est le stockage géologique de CO₂, dans des aquifères salins profonds. L'étanchéité des puits d'injection de CO₂ supercritique doit cependant être garantie pour des durées pluriséculaires. En raison de la pression d'injection élevée, la répartition des contraintes dans cette région du réservoir peut être modifiée de manière significative [1]. D'un point de vue mécanique, l'intégrité de toute la zone d'injection se trouve alors mise en jeu, et pourrait ne plus garantir l'étanchéité du puits. Notre objectif est d'analyser les modifications structurelles causées par l'interaction avec les deux fluides (eau et CO₂) lors de l'injection. Plusieurs modèles physiques sont proposés dans la littérature pour des conditions saturées, tandis que peu d'approches, existent pour des conditions insaturées. En effet, la plupart des approches reposent sur des hypothèses simplificatrices telles qu'un squelette rigide [2], ou une phase gazeuse statique [3] à la pression atmosphérique, pour la gestion d'un écoulement dans un milieu fracturé. Dans cette optique, nous proposons une formulation triphasique complète reposant sur les équations de bilan de masse, de mouvement et des relations constitutives pour les phases fluides ainsi que la phase solide en adoptant une approche couplage fort que l'on traite de façon monolithique. La formulation est de type pression-pression où le déplacement de la phase solide, la pression de gaz et la pression capillaire sont les variables principales.

## 2 Cadre de modélisation

Un milieu poreux partiellement saturé se comporte comme un milieu multiphasique où une phase mouillante et non mouillante coexistent dans les pores du squelette solide. Nous supposons dans notre modèle que les fluides sont immiscibles. Comme le volume poreux est entièrement occupé par les phases présentes, la somme des proportions de phases vaut toujours :

$$S_w + S_g = 1. \qquad (1)$$

<footnote n="1">¹ LMGC, Université Montpellier 2, France, {frederic.dubois, andre.chrysochoos} @univ-montp2.fr.</footnote>

<footnote n="2">² IFPEN, Dept. Mécanique des solides, {walid.saber, laurent.cangemi, nadege.brusselle} @ifpen.fr.</footnote>

<footnote n="3">³ MIST, IRSN-CNRS, Université Montpellier 2, France.</footnote>
À l'équilibre, les pressions d'eau (liquide) et de gaz sont reliées par la relation suivante :

$$p_c = p_g - p_w. \qquad (2)$$

La pression capillaire entre les deux phases joue un rôle déterminant dans l'état d'équilibre du système multiphasique. Les contraintes mécaniques sont exprimées grâce au concept de contrainte effective [4], où la compressibilité du squelette solide est prise en compte :

$$\sigma'' = \sigma + b\mathbf{m}\rho, \qquad (3)$$

où $\sigma$ est le tenseur de contrainte totale, $\sigma''$ est le tenseur de contrainte effective, $\mathbf{m}$ est la matrice identité, $p$ est une pression égale à la somme des pressions pondérées des deux phases $p = S_w p_w + S_g p_g$, et $b$ est la constante de Biot définit par $b = 1 - K_t/K_s$ ($K_t$ étant le module de compressibilité du milieu poreux et $K_s$ celui du squelette solide).

Les équations du système multiphasique sont complétées par les bilans de masse pour chaque phase. L'équation de conservation de la masse pour le solide peut s'écrire ainsi :

$$\frac{\partial}{\partial t}((1-n)\rho_s) + \nabla \cdot ((1-n)\rho_s \mathbf{u}_s) = 0 \qquad (4)$$

où $\rho_s$ est la densité du solide, $n$ la porosité du milieu et $\mathbf{u}_s$ est la vitesse du solide. De manière analogue, on peut écrire le bilan de masse pour la phase fluide et gaz :

$$\frac{\partial}{\partial t}(nS_{\alpha}\rho_{\alpha}) + \nabla \cdot (nS_{\alpha}\rho_{\alpha}\mathbf{u}_{\alpha}) = 0, \quad \alpha = w, g. \qquad (5)$$

En utilisant les équations d'états de chaque phase sous des conditions isothermes, on obtient :

$$\begin{gathered} \frac{1}{\rho_{\alpha}} \frac{\partial \rho_{\alpha}}{\partial t} = \frac{1}{K_{\alpha}} \frac{\partial p_{\alpha}}{\partial t} \qquad (6) \\ \frac{1}{\rho_s} \frac{\partial \rho_s}{\partial t} = \frac{[ (b-n) \frac{1}{K_s} \frac{\partial p}{\partial t} - (1-b) \nabla \cdot \mathbf{u} ]}{(1-n)} \end{gathered}$$

où $K_w$ et $K_g$ sont les modules de compressibilité de chaque fluide. La forme complète de ces équations d'état peuvent être trouvées dans [5].

Un choix approprié des variables principales est primordiales pour la modélisation et la résolution de problème fortement non-linéaire d'écoulement multiphasique dans un milieu poreux. Le choix des variables principales peut non seulement avoir un impact sur la performance globale du calcul mais aussi sur la faisabilité même de la simulations numérique. Dans la suite, les équations décrivant notre problème sont formulées en déplacement et en pressions. La pression de la phase non-mouillante (gaz) et la pression capillaire sont choisies. Le choix de la pression capillaire comme troisième variable est préférable par exemple, à la saturation en eau d'un point de vue numérique. Le choix de la saturation en eau comme variable indépendante introduit $\partial p_c / \partial S_w$ dans la formulation, tandis que la formulation basée sur le choix de la pression capillaire comme troisième variable introduit $\partial S_w / \partial p_c$ dans les équations. Dans une formulations avec la saturation en eau, la pression capillaire tends vers l'infini pour des valeurs de la saturation proche de zéro (condition de désaturation d'une zone). Ceci rend la simulation plus sensible au problème de convergence et d'instabilité numérique. En revanche, une formulation en pression capillaire assure un meilleur conditionnement du problème et de meilleures performances sur toute la plage de désaturation du milieu.

En injectant (6) dans (4) et (5) on obtient la forme finale des équations de conservation de la masse pour l'écoulement d'eau et de gaz à travers un milieu poreux déformable isotherme :

$$\left( \frac{b-n}{K_s} S_w + \frac{nS_w}{K_w} \right) \dot{\mathbf{p}}_g - \left( \frac{b-n}{K_s} S_w \left( S_w + p_c \frac{\partial S_w}{\partial p_c} \right) - n \frac{\partial S_w}{\partial p_c} + \frac{nS_w}{K_w} \right) \dot{\mathbf{p}}_c + (bS_w \nabla \cdot \mathbf{v}) + \frac{1}{\rho_w} \nabla \cdot \left[ \rho_w \frac{\mathbf{k} k_{r,w}}{\mu_w} (\nabla(p_g - p_c) - \rho_w (\mathbf{g} - \mathbf{u})) \right] = 0 \qquad (7)$$
$$ \left(\frac{b-n}{K_s}S_g + \frac{nS_g}{K_g}\right) \dot{\mathbf{p}}_g + \left(\frac{b-n}{K_s}S_g\left(S_w - p_c\frac{\partial S_w}{\partial p_c}\right) - n\frac{\partial S_w}{\partial p_c}\right) \dot{\mathbf{p}}_c + (bS_g \nabla \cdot \mathbf{v}) + \frac{1}{\rho_w} \nabla \cdot \left[\rho_g \frac{\mathbf{k}_{k,r,g}}{\mu_g}(\nabla p_g - \rho_g(\mathbf{g} - \mathbf{u}\mathbf{u}))\right] = 0 \quad (8) $$

La description du problème est complétée par l'écriture de la conservation de la quantité de mouvement. En considérant les forces d'inerties et les forces visqueuses, nous pouvons écrire :

$$ \nabla \cdot \sigma + \rho g = \rho \ddot{\mathbf{u}} + \zeta \dot{\mathbf{u}} \quad (9) $$

où $\zeta$ est le coefficient d'amortissement et $\rho$ la masse volumique moyenne du milieu, laquelle peut s'exprimer sous la forme suivante :

$$ \rho = (1-n)\rho_s + n(S_w\rho_w + S_g\rho_g). \quad (10) $$

## 2.1 Modélisation et formulation du problème couplé

## 2.2 Forme faible du problème couplé

Afin d'effectuer une discrétisation par éléments finis, nous écrivons la formulation variationnelle du problème couplé. Cette formulation s'obtient en sommant les formes faibles des problèmes d'élasticité et de Darcy, et en prenant en compte les conditions de couplage. La formulation faible du problème couplé consiste à trouver $(\mathbf{u}, \mathbf{p}) \in \mathcal{U} \times \mathcal{P}$ tel que :

$$ \int_{\Omega} (\nabla \mathbf{w})^{\mathrm{T}} \sigma'' d\Omega - \int_{\Omega} (\nabla \mathbf{w})^{\mathrm{T}} b \mathbf{m} S_w p_w d\Omega - \int_{\Omega} (\nabla \mathbf{w})^{\mathrm{T}} b \mathbf{m} S_g p_g d\Omega = \int_{\Omega} \mathbf{w}^{\mathrm{T}} \rho g d\Omega + \int_{\Gamma_u} \mathbf{w}^{\mathrm{T}} \tau d\Gamma. \quad (11) $$

On écrit de la même façon pour la conservation de la masse des fluides :

$$ \int_{\Omega} \mathbf{w}^{*\mathrm{T}} \left[ \left( \frac{b-n}{K_s} S_w + \frac{nS_w}{K_w} \right) \frac{\partial p_g}{\partial t} \right] d\Omega + \int_{\Omega} \mathbf{w}^{*\mathrm{T}} \left\{ \left[ -\frac{b-n}{K_s} S_w \left( S_w + p_c \frac{\partial S_w}{\partial p_c} \right) + n \frac{\partial S_w}{\partial p_c} - \frac{nS_w}{K_w} \right] \frac{\partial p_c}{\partial t} \right\} d\Omega \\ + \int_{\Omega} \mathbf{w}^{*\mathrm{T}} (b S_w \nabla \cdot \mathbf{v}_s) d\Omega - \int_{\Omega} (\nabla \mathbf{w}^{*\mathrm{T}}) \left( \frac{\mathbf{k}_{k,r,l}}{\mu_l} \nabla p_l \right) d\Omega + \int_{\Gamma_{q,w}} \mathbf{w}^{*\mathrm{T}} \left( \frac{\mathbf{q}_w}{\rho_w} \right) d\Gamma - \int_{\Omega} \mathbf{w}^{*\mathrm{T}} \left( \frac{\mathbf{k}_{k,r,w}}{\mu_w} \rho_w \mathbf{g} \right) d\Omega = 0, \quad (12) $$

et

$$ \int_{\Omega} \mathbf{w}^{*\bullet\mathrm{T}} \left[ \left( \frac{b-n}{K_s} S_g + \frac{nS_g}{K_g} \right) \frac{\partial p_g}{\partial t} \right] d\Omega + \int_{\Omega} \mathbf{w}^{*\bullet\mathrm{T}} \left\{ \left[ \frac{b-n}{K_s} S_g \left( S_w - p_c \frac{\partial S_w}{\partial p_c} \right) - n \frac{\partial S_w}{\partial p_c} \right] \frac{\partial p_c}{\partial t} \right\} d\Omega \\ + \int_{\Omega} \mathbf{w}^{*\bullet\mathrm{T}} (b S_g \nabla \cdot \mathbf{v}) d\Omega - \int_{\Omega} (\nabla \mathbf{w}^{*\bullet\mathrm{T}}) \left( \frac{\mathbf{k}_{k,r,g}}{\mu_g} \nabla p_g \right) d\Omega + \int_{\Gamma_{q,g}} \mathbf{w}^{*\bullet\mathrm{T}} \left( \frac{\mathbf{q}_g}{\rho_g} \right) d\Gamma - \int_{\Omega} \mathbf{w}^{*\bullet\mathrm{T}} \left( \frac{\mathbf{k}_{k,r,g}}{\mu_g} \rho_g \mathbf{g} \right) d\Omega = 0. \quad (13) $$

Cela pour toutes fonctions $\mathbf{w}$, $\mathbf{w}^*$ et $\mathbf{w}^\bullet$ suffisamment régulières, et où $\mathcal{U} = \{\mathbf{u}, \text{ tel que } \mathbf{u} \in H^1(\Omega), \mathbf{u} = \mathbf{u}_0 \text{ sur } \Gamma_u \}$ et $\mathcal{P} = \{\mathbf{p}, \text{ tel que } \mathbf{p} \in H^1(\Omega)\}$.

## 2.3 Discrétisation par éléments finis

Le système composé des équations (11), (12) et (13) est résolu en utilisant la méthode des éléments finis. La particularité de l'approche développée ici, est d'utiliser une méthode de résolution
directe (ou monolithique) du problème couplé. Les fonctions de forme utilisées pour la discrétisation du système sont les suivantes :

$$ \mathbf{u} = \mathbf{N}_u \bar{\mathbf{u}}, \mathbf{p}_c = \mathbf{N}_p \bar{\mathbf{p}}_c, \mathbf{p}_g = \mathbf{N}_p \bar{\mathbf{p}}_g \quad (14) $$

Le déplacement et la pression sont approchées par des fonctions continues et linéaires par élément. Il faut noter que seul un état incompressible ou non drainée permet d'utiliser des fonctions d'interpolation de même ordre à la fois pour le déplacement et la pression de pores. Dans le cas contraire, il est nécessaire d'utiliser des fonctions de forme linéaires en pression et quadratiques en déplacement pour satisfaire aux conditions de stabilités numérique et éviter les problèmes d'oscillation et de verrouillage numérique.

Après discrétisation des équations non-linéaires aux dérivées partielles en utilisant les fonctions de forme (14) pour approcher le déplacement, la pression de l'eau et la pression capillaire, on peut écrire le système couplé d'équations algébriques sous forme matricielle :

$$ \mathbf{M}_u \ddot{\mathbf{u}} + \mathbf{C} \dot{\mathbf{u}} + \mathbf{K}_u \mathbf{u} - \mathbf{C}_{sw} \dot{\mathbf{p}}_c - \mathbf{C}_{sg} \dot{\mathbf{p}}_g = \mathbf{F}_u \quad (15) $$

$$ \mathbf{M}_w \ddot{\mathbf{u}} + \mathbf{C}_{sw} \dot{\mathbf{u}} + \mathbf{P}_{ww} \dot{\mathbf{p}}_c + \mathbf{C}_{wg} \mathbf{P}_g + \mathbf{H}_{ww} \mathbf{p}_c = \mathbf{F}_c \quad (16) $$

$$ \mathbf{M}_g \ddot{\mathbf{u}} + \mathbf{C}_{sg} \dot{\mathbf{u}} + \mathbf{C}_{gw} \dot{\mathbf{p}}_c + \mathbf{P}_{gg} \dot{\mathbf{p}}_g + \mathbf{H}_{gg} \mathbf{p}_g = \mathbf{F}_g \quad (17) $$

La définition des coefficients des matrices et des vecteurs de chargement est donnée ci-après :

$$
\begin{align*}
\mathbf{K}_l &= \int_{\Omega} \mathbf{B}^{\mathrm{T}} \sigma'' d\Omega \\
\mathbf{M}_u &= \int_{\Omega} \mathbf{N}_u^{\mathrm{T}} \rho \mathbf{N}_u d\Omega \\
\mathbf{M}_w &= \int_{\Omega} (\nabla \mathbf{N}_p)^{\mathrm{T}} \frac{\mathbf{k} k_{r,l}}{\mu_l} \rho_g \mathbf{N}_u d\Omega \\
\mathbf{M}_g &= \int_{\Omega} (\nabla \mathbf{N}_p)^{\mathrm{T}} \frac{\mathbf{k} k_{r,g}}{\mu_g} \rho_g \mathbf{N}_u d\Omega \\
\mathbf{C}_{sg} &= \int_{\Omega} \mathbf{B}^{\mathrm{T}} b (1 - S_w) \mathbf{N}_p d\Omega \\
\mathbf{C}_{sw} &= \int_{\Omega} \mathbf{B}^{\mathrm{T}} b S_w \mathbf{N}_p d\Omega \\
\mathbf{C}_{wg} &= \int_{\Omega} \mathbf{N}_p^{\mathrm{T}} \left( \frac{b-n}{K_s} S_w + \frac{n S_w}{K_w} \right) \mathbf{N}_p d\Omega \\
\mathbf{C}_{gw} &= \int_{\Omega} \mathbf{N}_p^{\mathrm{T}} \left[ \frac{b-n}{K_s} S_g \left( S_w - p_c \frac{\partial S_w}{\partial p_c} \right) - n \frac{\partial S_w}{\partial p_c} \right] \mathbf{N}_p d\Omega \tag{18} \\
\mathbf{P}_{ww} &= \int_{\Omega} \mathbf{N}_p^{\mathrm{T}} \left( - \frac{b-n}{K_s} S_w \left( S_w + p_c \frac{\partial S_w}{\partial p_c} \right) + n \frac{\partial S_w}{\partial p_c} - \frac{n S_w}{K_w} \right) \mathbf{N}_p d\Omega \\
\mathbf{P}_{gg} &= \int_{\Omega} \mathbf{N}_p^{\mathrm{T}} \left( \frac{b-n}{K_s} S_g + \frac{n S_g}{K_g} \right) \mathbf{N}_p d\Omega \\
\mathbf{H}_{ww} &= \int_{\Omega} (\nabla \mathbf{N}_p)^{\mathrm{T}} \frac{\mathbf{k} k_{r,l}}{\mu_l} (\nabla \mathbf{N}_p) d\Omega \\
\mathbf{H}_{gg} &= \int_{\Omega} (\nabla \mathbf{N}_p)^{\mathrm{T}} \frac{\mathbf{k} k_{r,g}}{\mu_g} \nabla \mathbf{N}_p d\Omega \\
\mathbf{F}_w &= \int_{\Gamma_l^q} \mathbf{N}_p^{\mathrm{T}} \left( \frac{\mathbf{q}_l^w}{\rho_l^w} \right) d\Gamma - \int_{\Omega} \mathbf{N}_p^{\mathrm{T}} \left( \frac{\mathbf{k} k_{r,l}}{\mu_l} \rho_w \mathbf{g} \right) d\Omega = 0 \\
\mathbf{F}_g &= \int_{\Gamma_{qg}} \mathbf{N}_p^{\mathrm{T}} \left( \frac{\mathbf{q}_g}{\rho_g} \right) d\Gamma - \int_{\Omega} \mathbf{N}_p^{\mathrm{T}} \left( \frac{\mathbf{k} k_{r,g}}{\mu_g} \rho_g \mathbf{g} \right) d\Omega = 0
\end{align*}
$$
## 2.4 Stratégie d'implémentation

Simuler le comportement de milieux aussi fortement couplé aussi fortement couplés peut s'avérer assez complexe. La plate-forme ouverte LMGC90 est pour cela bien adaptée à ce type d'exercice par sa modularité. Nous avons cherché pour le développement et l'implémentation de notre outil de simulation à factoriser dans un élément fini générique les différents opérateurs nécessaires à notre formulation. Cette approche évite la répétition de code ayant le même objectif ou le même principe et permet d'envisager des extensions au modèle (e.g. ajout des effets thermiques). L'enjeu de la conception d'un tel élément est de fournir une structure de données commune au niveau des opérateurs pour résoudre de tels problèmes multiphysiques. En effet, pour un grand nombre de problèmes physiques, les équations aux dérivées partielles ne diffèrent que par un certain nombre de paramètres matériels et donc les matrices et les vecteurs élémentaires résultant de la discrétisation par éléments finis sont similaires.

Les opérateurs résultant de la discrétisation par éléments finis des équations des formes faibles (11) (12) (13) peuvent être classifiés comme indiqué dans le tableau (1) ci-dessous :

<table><caption>Tableau 1 – Matrices et vecteurs type</caption><thead><tr><td>Type</td><td>Nom</td><td>Équations</td></tr></thead><tbody><tr><td>∫<sub>Ω</sub> N<sub>i</sub><sup>T</sup> M N<sub>i</sub><sup>T</sup> dΩ</td><td>matrice de masse</td><td>(11) (12) (13)</td></tr><tr><td>∫<sub>Ω</sub> ∇N<sub>i</sub><sup>T</sup> M ∇N<sub>i</sub><sup>T</sup> dΩ</td><td>matrice Laplacienne</td><td>(12) (13)</td></tr><tr><td>∫<sub>Ω</sub> B<sup>T</sup> M B dΩ</td><td>matrice tangente</td><td>(11)</td></tr><tr><td>∫<sub>Ω</sub> B<sup>T</sup> M N<sub>i</sub><sup>T</sup> dΩ</td><td>matrice de couplage déplacement</td><td>(11)</td></tr><tr><td>∫<sub>Ω</sub> N<sub>i</sub><sup>T</sup> M N<sub>i</sub><sup>T</sup> dΩ</td><td>matrice de couplage pression</td><td>(12) (13)</td></tr><tr><td>∫<sub>Ω</sub> N<sub>i</sub><sup>T</sup> q dΩ</td><td>vecteur second membre</td><td>(11) (12) (13)</td></tr></tbody></table>

## 2.5 Stratégie de prise en compte de l'endommagement

La difficulté lors de l'apparition de fissures est la gestion cohérente et performante des nouvelles conditions aux limites créées aux bords du modèle poreux. De plus si les fissures sont connectées, il convient alors d'être capable de gérer un problème d'écoulement et de le coupler avec le milieu poreux. Nous nous proposons ici de gérer le traitement de cette dégradation volumique de la structure par une méthode d'érosion de maillage proposée par [7]. Cette approche consiste à modéliser l'amorçage et la propagation de fissure sans remaillage du domaine, tout en conservant une formulation élément fini « standard » où l'évolution de l'endommagement est ici gouvernée par un critère énergétique reposant sur la minimisation de fonctionnelles [8]. Cette méthode ne nécessite pas de formulation spécifique à champs discontinus comme les méthodes de zones cohésives. Un des avantages de cette méthode est de pouvoir éventuellement utiliser les éléments « érodés » comme les éléments support au modèle d'écoulement. Un défaut de cette approche est la mauvaise prise en compte des contacts frottant entre les bords de la fissure.

## 3 Conclusion et perspective

Dans cet article, grâce à une formulation adaptée, nous avons présenté une stratégie numérique pour simuler un problème fortement couplé d'hydromécanique selon une approche monolithique. Une difficulté pour valider notre modèle est qu'il n'existe pas de solutions analytiques d'écoulement diphasique dans un milieu poreux déformable, et les résultats expérimentaux existant sont peu nombreux et peu documentés. L'un d'entre eux est le cas test de drainage expérimenté par Liakopoulos (1965), la confrontation de notre modèle avec celui-ci est en cours d'étude. Quelques simulations ont été rendues disponible récemment [6] dans la littérature, notamment ceux de [9] où une comparaison avec leurs travaux est prévu.
Références

[1] J. Nordbotten, M. Celia, *Similarity solutions for fluid injection into confined aquifers*, Journal of Fluid Mechanics 561, 2006, 307-327.

[2] J. Rethore, R. de Borst, MA. Abellan, *A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks*, Computational Mechanics, 2008, 42:227-238.

[3] BA. Schrefler, R. Scotta, *A fully coupled dynamic model for two-phase fluid flow in deformable porous media*, Comput Methods Appl Mech Eng, 2001, 190:3223-46.

[4] MA. Biot *A general theory of three dimensional consolidation*. Journal of Applied Physics, 1941; 12:155-164.

[5] R. W. Lewis and B. A. Schrefler, *The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media*. John Wiley & Sons 2nd Edition, 1998

[6] A.R. Khoei, T. Mohammadnejad, *Numerical modeling of multiphase fluid flow in deforming porous media : A comparison between two- and three-phase models for seismic analysis of earth and rockfill dams*, Computers and Geotechnics 2011, 38:142-166.

[7] A. Pandolfi, M. Ortiz, *Numerical modeling of multiphase fluid flow in deforming porous media : An eigen-erosion approach to brittle fracture*, International Journal For Numerical Methods in Engineering 2012.

[8] G.A. Francfort, J-J. Marigo. *Revisiting brittle fracture as an energy minimization problem*, Journal of the Mechanics and Physics of Solids 1998; 46(8):1319-1342.

[9] O. Kolditz, S. Bauer, N.Böttcher, D. Elsworth, *Numerical simulation of two-phase flow in deformable porous media : Application to carbon dioxide storage in the subsurface*, Mathematics and Computers in Simulation, 2012, 1919-1935.