1. Problème de référence#
L’objectif de ce cas test est de comparer la solution obtenue avec Code_Aster, en utilisant une modélisation aux éléments finis, avec une solution analytique. Ce cas test et sa solution analytique ont été développés et détaillés par Abousleiman et al.[1]. Nous en présentons ici les principales lignes ainsi que les analogies admissibles avec notre modèle.
1.1. Géométrie du problème#
On considère un volume de roche de hauteur \(h=\mathrm{1m}\) et de base carrée de \(l=\mathrm{15m}\) de coté. Au centre, est creusée une cavité cylindrique de rayon \(R=\mathrm{0.1m}\) et d’extension infinie selon son axe \(z\) . On remarque qu’on modélise ici un domaine de hauteur \(\mathrm{1m}\) ce qui est une hypothèse restrictive (idéalement, il faudrait un domaine bien plus grand mais qui entrainerait de longs temps de calculs). Cette cavité sert à modéliser le puits creusé dans le massif de roche. Ce massif est d’autre part caractérisé par des plans d’isotropie perpendiculaires à l’axe de la cavité cylindrique.
Le puits est d’autre part caractérisé par deux anglespar rapport au repère des contraintes principales:
son azimut \({\varphi}_{y}=30°\)
son inclinaison \({\varphi}_{z}=60°\)
Ce qui correspond aux angles nautiques :
\(\alpha =30°\)
\(\beta =-60°\)
Le repère cartésien associé au massif et au repère des contraintes est noté \((Ox'y'z')\) et celui associé à la cavité est noté \((\mathit{Oxyz})\) . On définit également un système de coordonnées cylindriques \((r,\theta ,z)\) pour caractériser la cavité cylindrique.
Sur la Figure on résume les informations décrites ci-dessus:
Figure 1.1-1 : Représentation de la géométrie du problème du puits incliné
Il est important de remarquer que les deux géométries décrites ci-dessus sont strictement équivalentes. Cependant on choisit d’adopter la géométrie pour laquelle le puits est vertical, afin de faciliter la modélisation et la conception du maillage. Le repère choisi est donc le repère d’isotropie transverse et non celui des contraintes principales.
Au final, on représentera une partie seulement du modèle tel que:
Figure 1.1-2 : Géométrie modélisée
Des conditions de symétries sont donc appliquées.
1.2. Propriétés matériaux#
On se réfère aux paramètres définis dans l’article d’Abou Sleiman 1 .
On part du rapport d’anisotropie \(\frac{{E}_{L}}{{E}_{N}}=2\) et \(\frac{{\nu}_{\mathrm{TL}}}{{\nu}_{\mathrm{NL}}}=2\) . On rappelle que les hypothèses d’isotropie transverse impliquent que \({\nu}_{T}={\nu}_{L}\) ; \({\nu}_{N}={\nu}_{N}\) ; \({\nu}_{L}={\nu}_{T}\) et \(\frac{{\nu}_{N}}{{\nu}_{L}}=\frac{{E}_{N}}{{E}_{L}}\) .
Par ailleurs et en écrivant les termes de la matrice de Hook, tels que :
\({M}_{\text{}11}=\frac{{E}_{\mathrm{L.}}({E}_{N}-{E}_{L}{\nu}_{N}^{2})}{(1+{\nu}_{T})({E}_{N}-{E}_{\mathrm{N.}}{\nu}_{T}-{2.E}_{\mathrm{L.}}{\nu}_{N}^{2})}\)
\({M}_{\text{}12}=\frac{{E}_{L}({E}_{N}{\nu}_{T}+{E}_{L}{\nu}_{N}^{2})}{(1+{\nu}_{T})({E}_{N}-{E}_{\mathrm{N.}}{\nu}_{T}-{2.E}_{\mathrm{L.}}{\nu}_{N}^{2})}\)
\({M}_{\text{}13}=\frac{{E}_{L}{E}_{N}{\nu}_{N}}{{E}_{N}-{E}_{N}{\nu}_{T}-{2.E}_{L}{\nu}_{N}^{2}}\)
\({M}_{\text{}33}=\frac{{E}_{N}^{2}(1-{\nu}_{T})}{{E}_{N}-{E}_{N}{\nu}_{T}-{2.E}_{L}{\nu}_{N}^{2}}\)
Les coefficients de Biot sont alors déduits par la relation suivante 1 :
\({b}_{L}=1-\frac{{M}_{\text{}11}+{M}_{\text{}12}+{M}_{\text{}13}}{{\mathrm{3K}}_{S}}\)
et
\({b}_{N}=1-\frac{{\mathrm{2M}}_{\text{}13}+{M}_{\text{}33}}{{\mathrm{3K}}_{S}}\)
Ici, \({K}_{S}=27,5\mathrm{GPa}\) .
À partir de cela, on indique ici les propriétés matériaux renseignés dans Code_Aster :
Liquide |
Viscosité \({\mu}_{w}\) (en \(\mathit{Pa.s}\) ) Module de compressibilité \(\frac{1}{{K}_{w}}\) (en \({\mathit{Pa}}^{-1}\) ) Densité du liquide \({\rho}_{w}\) (en \(\mathit{kg}/{m}^{3}\) ) Chaleur massique à pression constante |
\({10}^{\text{-3}}\) \({5.10}^{\text{-10}}\) \(1\) \(4180\) |
Paramètres élastiques |
Module de Young \({E}_{L}\) (en \(\mathit{MPa}\) ) Module de Young \({E}_{N}\) (en \(\mathit{MPa}\) ) Coefficient de Poisson \({\nu}_{\text{LT}}\) Coefficient de Poisson \({\nu}_{\text{LN}}\) Module de Cisaillement \({G}_{\ln}\) (en \(\mathit{MPa}\) ) |
\(9474\) \(0.5\times 9474\) \(0,24\) \(0,24\) \(8880\) |
Paramètres de couplage |
Coefficient de Biot \({b}_{L}\) Coefficient de Biot \({b}_{N}\) Masse volumiques homogénéisée initiale \({r}_{0}\) (en \(\mathit{kg}/{m}^{3}\) ) Perméabilité intrinsèque \({K}_{L}^{int}\) (en \({m}^{2}/s\) ) Perméabilité intrinsèque \({K}_{N}^{int}\) (en \({m}^{2}/s\) ) |
\(0,81689\) \(0,89864\) \(2410\) \({5.10}^{\text{-20}}\) \({5.10}^{\text{-20}}\) |
Tableau 1.2-1 : Propriétés des matériaux
1.3. Conditions aux limites et initiales du modèle#
Dans la suite nous considérons les conditions initiales et aux limites dans le repère lié au puits.
1.3.1. Cas des déplacements#
Compte tenu de la symétrie du modèle, on bloquera sur 3 faces correspondant respectivement aux 3 plans \(x=0;y=0\) et \(z=0\) .
1.3.2. Cas de la pression#
La pression de pore initiale dans le massif et pour \(r\to +\infty ` est fixée à :math:`{p}_{0}=9,8\mathit{MPa}\) . On appliquera donc cette pression sur les 3 faces correspondant respectivement aux 3 plans \(x=l;y=l\) et \(z=h\) . Suite à la mise en eau du puits, il se développe en paroi une pression uniforme que l’on fixe à \({p}_{w}=12\mathit{MPa}\) sur le bord de la cavité.
1.3.3. Cas des contraintes#
Dans le massif, l’état de contrainte in situ associé au repère \((Ox'y'z')\) avant le creusement du puits est caractérisé par l’état initial suivant:
\(\lbrace \begin{array}{c}S{'}_{x}=-25\mathit{MPa}\\ S{'}_{y}=-22\mathit{MPa}\\ S{'}_{z}=-29\mathit{MPa}\end{array}\)
De plus le creusement du puits implique une perturbation de l’état de contrainte in situ initial dans le massif. Afin de prendre en compte ce nouvel état de contrainte au voisinage du puits, il faut considérer deux rotations successives pour passer du repère \((Ox'y'z')\) associé au massif, au repère \((O\mathit{xyz})\) associé au puits (cf. Figure ).
Soient \({P}_{1}\) et \({P}_{2}\) les matrices de passage décrivant ce changement de repère:
\({P}_{1}=(\begin{array}{ccc}\cos(\alpha )& -\sin(\alpha )& 0\\ \sin(\alpha )& \cos(\alpha )& 0\\ 0& 0& 1\end{array})\) \({P}_{2}=(\begin{array}{ccc}\cos(\beta )& 0& -\sin(\beta )\\ 0& 1& 0\\ \sin(\beta )& 0& \cos(\beta )\end{array})\)
Après changement de repère, l’état de contrainte initiale dans le repère \((\mathit{Oxyz})\) et quand :math:`rto +infty ` à appliquer au modèle est le suivant:
\(\lbrace \begin{array}{c}{S}_{xx}=-27,81\mathrm{MPa}\\ {S}_{yy}=-22,75\mathrm{MPa}\\ {S}_{zz}=-25,44\mathrm{MPa}\\ {S}_{xy}=0,65\mathrm{MPa}\\ {S}_{xz}=2,06\mathrm{MPa}\\ {S}_{yz}=1,125\mathrm{MPa}\end{array}\)
Pour les conditions aux limites, on néglige le cisaillement et on applique les conditions en pression sur les 3 faces \(x=l;y=l\) et \(z=h\) .
En paroi du puits \((r=R)\) et suivant , la contrainte radiale est appliquée : \({\sigma}_{\mathrm{rr}}(r=R,t)=-{p}_{w}.H(t)\) où \(H(t)\) est la fonction échelon de Heaviside.
1.4. Solution de référence#
Dans cette partie nous rappelons les équations théoriques utiles pour le calcul de la contrainte radiale, et de la pression de pore.
On résume les conditions aux limites du problème général :
pour :math:`rto +infty `
\(\lbrace \begin{array}{c}{\sigma}_{xx}=-{S}_{xx}\\ {\sigma}_{yy}=-{S}_{yy}\\ {\sigma}_{zz}=-{S}_{zz}\\ {\sigma}_{xy}=-{S}_{xy}\\ {\sigma}_{xz}=-{S}_{xz}\\ {\sigma}_{yz}=-{S}_{yz}\end{array}\) et \(p={p}_{0}\)
pour \(r=R\)
\({\sigma}_{\mathit{rr}}=-{p}_{w}.H(t)\)
\({\sigma}_{r\theta }={\sigma}_{\mathit{rz}}=0\)
\(p={p}_{w}.H(t)\)
Compte tenu de la linéarité du problème, celui ci est en fait le résultat de la superposition de 3 sous problèmes qui permettent une étude tridimensionnelle 1 :
problème 1: on considère un problème poro-élastique en déformations planes (on ne prend en compte que les contraintes principales selon les axes \(x\) et \(y\) ainsi que le cisaillement dans la plan \((xy)\) et les perturbations liées à la pression de pore ). Les conditions aux limites de ce problème s’écrivent alors:
pour :math:`rto +infty `
\({\sigma}_{xx}=-{S}_{xx}\)
\({\sigma}_{yy}=-{S}_{yy}\)
\({\sigma}_{xy}=-{S}_{xy}\)
\({\sigma}_{zz}=-{\nu}_{\text{LN}}({S}_{xx}+{S}_{yy})-({\alpha}_{\text{LN}}-2{\nu}_{\text{LN}}{\alpha}_{\text{LT}}){p}_{0}-({\beta}_{3}^{s}-2{\nu}_{\text{LN}}{\beta}_{1}^{s}){T}_{0}\)
\({\sigma}_{yz}={\sigma}_{xz}=0\)
\(p={p}_{0}\)
pour \(r=R\)
\({\sigma}_{\mathit{rr}}=-{p}_{w}.H(t)\)
\(p={p}_{w}.H(t)\)
problème 2: on ne considère que la contrainte verticale orientée selon \(z\) . Les conditions aux limites de ce problème s’écrivent alors:
pour :math:`rto +infty `
\({\sigma}_{xx}={\sigma}_{yy}={\sigma}_{xy}={\sigma}_{xz}={\sigma}_{yz}=p=0\)
\({\sigma}_{zz}=-{S}_{zz}+({\nu}_{\text{LN}}({S}_{xx}+{S}_{yy})+({\alpha}_{\text{LN}}-2{\nu}_{\text{LN}}{\alpha}_{\text{LT}}){p}_{0})\)
\({\beta}_{3}^{s}\) et \({\beta}_{1}^{s}\) sont définis dans la suite de ce document.
pour \(r=R\)
\({\sigma}_{\mathrm{rr}}={\sigma}_{r\theta }={\sigma}_{\mathrm{rz}}=p=0\)
problème 3: on ne considère que les cisaillements contenus dans les plans \((xz)\) et \((yz)\) . Les conditions aux limites de ce problème s’écrivent alors:
pour :math:`rto +infty `
\({\sigma}_{xx}={\sigma}_{yy}={\sigma}_{zz}={\sigma}_{xy}=p=0\)
\({\sigma}_{xz}=-{S}_{xz}\)
\({\sigma}_{yz}=-{S}_{yz}\)
pour \(r=R\)
\({\sigma}_{\mathrm{rr}}={\sigma}_{r\theta }={\sigma}_{\mathrm{rz}}=p=0\)
Dans la suite de ce document on va développer la solution du problème 1.
1.4.1. Présentation du problème analytique#
L’expression de la contrainte radiale, de la pression de pore et de la température, pour le problème en déformations planes (problème 1 ci-dessus) est également le résultat de la superposition de 3 modes de chargement dont les conditions aux limites sont données ci-après:
Mode 1: prise en compte de la partie hydrostatique des contraintes aux limites (en \(r=R\) )
\({\sigma}_{\mathit{rr}}^{(1)}={\sigma}_{0}-{p}_{w}\)
\({p}^{(1)}=0\)
Mode 2: prise en compte des perturbations liées à la pression de pore et à la température (en \(r=R\) )
\({\sigma}_{\mathit{rr}}^{(2)}=0\)
\({p}^{(2)}={p}_{w}-{p}_{0}\)
Mode 3: prise en compte de la partie déviatorique des contraintes aux limites (en \(r=R\) )
\({\sigma}_{\mathit{rr}}^{(3)}=-{S}_{0}\cos(2(\theta -{\theta}_{r}))\)
\({p}^{(3)}=0\)
avec:
\({\sigma}_{0}=\frac{({S}_{xx}+{S}_{yy})}{2}\) la contrainte moyenne, \({S}_{0}=0.5\sqrt{({({S}_{xx}-{S}_{yy})}^{2}+4{S}_{xy}^{2})}\) la contrainte déviatorique et
\(\text{}{\theta}_{r}=0.5\arctan(\frac{2{S}_{xy}}{{S}_{xx}-{S}_{yy}})\)
La solution générale pour le calcul de la contrainte radiale est donc donnée par l’expression suivante:
\({\sigma}_{\mathit{rr}}=-{\sigma}_{0}+{S}_{0}\cos(2(\theta -{\theta}_{r}))+{\sigma}_{\mathit{rr}}^{(1)}+{\sigma}_{\mathit{rr}}^{(2)}+{\sigma}_{\mathit{rr}}^{(3)}\)
avec:
\({\sigma}_{\mathit{rr}}^{(1)}=({\sigma}_{0}-{p}_{w})\frac{{R}^{2}}{{r}^{2}}\)
\({\sigma}_{\mathit{rr}}^{(2)}={L}^{-1}\left[\frac{1}{s}({b}_{L}(1-\frac{{M}_{12}}{{M}_{11}})\left[{F}_{1}\psi (\xi )+{F}_{2}\psi (\omega )\right]+{\beta}_{1}^{s}(1-\frac{{M}_{12}}{{M}_{11}})\left[({T}_{w}-{T}_{0})\right]\psi (\omega ))\right]\)
\({\sigma}_{\mathit{rr}}^{(3)}={L}^{-1}\left[\frac{{S}_{0}}{s}\cos(2(\theta -{\theta}_{r}))({A}_{1}{C}_{1}(\frac{1}{\xi r}{K}_{1}(\xi r)+\frac{1}{{(\xi r)}^{2}}{K}_{2}(\xi r))-{A}_{2}{C}_{2}\frac{{R}^{2}}{{r}^{2}}-3{A}_{3}{C}_{3}\frac{{R}^{4}}{{r}^{4}})\right]\)
Les quantités utilisées sont définis à la fin de cette section.
La solution générale pour le calcul de la pression de pore est donc donnée par l’expression suivante:
\(p={p}_{0}+{p}^{(2)}+{p}^{(3)}\)
avec:
\({p}^{(2)}={L}^{-1}\left[\frac{1}{s}\left[{F}_{1}\phi (\xi )+{F}_{2}\phi (\omega )\right]\right]\)
\({p}^{(3)}={L}^{-1}\left[\frac{{S}_{0}}{s}\cos(2(\theta -{\theta}_{r}))({A}_{1}{C}_{2}\frac{{R}^{2}}{{r}^{2}}+\frac{{c}_{f}}{2{G}_{\text{LT}}\kappa }{C}_{1}{K}_{2}(\xi r))\right]\)
Les quantités utilisées sont définis à la fin de cette section.
D’autre part on définit \({K}_{n}\) comme étant les fonctions modifiées de Bessel de deuxième espèce d’ordre \(n\) (où \(n\) est un entier) telles que:
\({K}_{n}(x)=\begin{array}{c}\frac{1}{2}{(\frac{1}{2}x)}^{-n}{\sum}_{k=0}^{n-1}\left[\frac{(n-k-1)!}{k!}{(\frac{-1}{4}{x}^{2})}^{k}\right]+{(-1)}^{n+1}\ln\left[\frac{1}{2}x\right]{I}_{n}(x)\\ +{(-1)}^{n}\frac{1}{2}{(\frac{1}{2}x)}^{n}{\sum}_{k=0}^{\infty}[(\chi (k+1)+\chi (n+k+1))\frac{{\left[\frac{1}{4}{x}^{2}\right]}^{k}}{k!(n+k)!}]\end{array}\)
où \({I}_{n}(x)=\frac{1}{\pi}{\int}_{0}^{\pi}{e}^{\mathit{xcos}\theta }\cos(n\theta )d\theta\) sont les fonctions modifiées de Bessel de première espèce d’ordre \(n\) (où \(n\) est un entier) et \(\chi\) la fonction digamma.
On peut donc définir les quantités suivantes:
\(\phi (x)=\frac{{K}_{0}(\mathit{xr})}{{K}_{0}(\mathit{xR})}\)
\(\psi (x)=\frac{{K}_{1}(\mathit{xr})}{{\mathit{xrK}}_{0}(\mathit{xR})}-\frac{{\mathit{RK}}_{1}(\mathit{xr})}{{\mathit{xr}}^{2}{K}_{0}(\mathit{xR})}\)
\({F}_{1}=({p}_{w}-{p}_{0})-\frac{{c}_{\mathit{hf}}}{1-\frac{{c}_{f}}{{c}_{h}}}({T}_{w}-{T}_{0})\)
\({F}_{2}=\frac{{c}_{\mathit{hf}}}{1-\frac{{c}_{f}}{{c}_{h}}}({T}_{w}-{T}_{0})\)
\({c}_{f}=\frac{\kappa M{M}_{11}}{{M}_{11}+{b}_{L}^{2}M}\)
\(\kappa =\frac{{K}_{L}^{int}(\phi )}{{\mu}_{w}}\)
\({c}_{\mathit{hf}}=\frac{{c}_{f}}{\kappa}\left[{\beta}^{\mathit{sf}}-\frac{{b}_{L}{\beta}_{1}^{S}}{{M}_{11}}\right]\)
\({C}_{1}=\frac{4}{{\mathrm{2A}}_{1}({B}_{3}-{B}_{2})-{A}_{2}{B}_{1}}\)
\({C}_{1}=\frac{4}{{\mathrm{2A}}_{1}({B}_{3}-{B}_{2})-{A}_{2}{B}_{1}}\)
\({C}_{3}=\frac{{\mathrm{2A}}_{1}({B}_{3}+{B}_{2})+{\mathrm{3A}}_{2}{B}_{1}}{3({\mathrm{2A}}_{1}({B}_{3}-{B}_{2})-{A}_{2}{B}_{1})}\)
\({A}_{1}=\frac{{b}_{L}M}{{M}_{11}+{b}_{L}^{2}M}\)
\({A}_{2}=\frac{{M}_{11}+{M}_{12}+2{b}_{L}^{2}M}{{M}_{11}+{b}_{L}^{2}M}\)
\({B}_{1}=\frac{{M}_{11}}{2{G}_{\text{LT}}{b}_{L}}{K}_{2}(\xi R)\)
\({B}_{2}=\frac{1}{\xi R}{K}_{1}(\xi R)+\frac{6}{{(\xi R)}^{2}}{K}_{2}(\xi R)\)
\({B}_{3}=2(\frac{1}{\xi R}{K}_{1}(\xi R)+\frac{3}{{(\xi R)}^{2}}{K}_{2}(\xi R))\)
\(\xi =\sqrt{\frac{s}{{c}_{f}}}\) et \(\omega =\sqrt{\frac{s}{{c}_{h}}}\) où \(s\) représente la variable de Laplace.
1.4.2. Méthode de résolution de la solution analytique#
Comme nous avons pu le constater précédemment, les solutions pour les modes 2 et 3 sont exprimés dans l’espace de Laplace. Compte tenu de la complexité de l’inversion de la solution vers l’espace temporel, nous avons recours à un algorithme d’inversion dit de Stehfest-Graver. Cette méthode est reconnue pour sa stabilité numérique.
Le principe de la méthode consiste à approximer la fonction à inverser \(\tilde{f}(s)\) à l’aide d’une delta-série en utilisant la formule suivante:
\(f(t)\mathrm{\simeq }\frac{\ln(2)}{t}{\sum}_{i=1}^{N}{V}_{i}\tilde{f}(i\frac{\ln(2)}{t})\)
où \({V}_{i}={(-1)}^{i+(N/2)}\frac{{\sum}_{k=E((i+1)/2)}^{min(i,N/2)}{k}^{N/2}(\mathrm{2k})!}{(\frac{N}{2}-k)!k!(k-1)!(i-k)!(\mathrm{2k}-i)!}\)
avec \(E(x)\) la fonction partie entière et \(N\) le nombre de termes de la série. Il est en général compris entre 10 et 20. Dans notre cas nous avons choisi \(N=18\)