1. Problème de référence#
L’objectif de ce cas test est de comparer la solution obtenue avec les différents schémas volumes finis à une solution analytique.
1.1. Solution analytique#
Le problème monophasique instationnaire et monodimensionnel peut s’écrire sous une forme générale du type :
\(N\frac{\partial P}{\partial t}-{K}_{\mathit{int.}}\Delta P=0\)
\(P(t=0)={P}_{0}\)
\(P(t,x=0)=0\)
\(\frac{\partial P}{\partial x}(t,x=L)=0\)
Ce problème admet une solution analytique obtenue par développement en séries de Fourier.
\(P={\sum}_{k=0}^{\infty}\frac{{\mathrm{4P}}_{0}}{(\mathrm{2k}+1)\pi }\exp(\frac{-{K}_{i}}{N}{\omega}_{k}^{2}t)\sin({\omega}_{k}x)\) avec \({\omega}_{k}=(k+\frac{1}{2})\frac{\pi}{L}\)
On peut réduire cette série à un nombre fini \(K\) de termes, en fonction de l’instant calculé.
Ce nombre de termes est déterminé de la façon suivante:
Soit \({n}_{x}\) le nombre de points \({x}_{i}\) où la solution est évaluée à un instant \(t\)
On pose: \({a}_{k}^{i}=\frac{4}{(\mathrm{2k}+1)\pi }\exp(\frac{-{K}_{int}}{N}{\omega}_{k}^{2}t)\sin({\omega}_{k}{x}_{i})\)
Si bien que la solution peut s’écrire: \(P({x}_{i})={\sum}_{k=0}^{K}{P}_{0}.{a}_{k}^{i}\)
On choisit \(K\) tel que: \(\frac{1}{{n}_{x}}\sqrt{{\sum}_{i=1}^{\mathit{nx}}{({a}_{k}^{i})}^{2}}<\epsilon\)
En pratique, nous avons pris \(\epsilon ={10}^{-10}\) .
Les allures de la solution analytique aux temps 1, 10, 100, 1000 sont montrées sur la figure :
Figure 1: Représentation des solutions analytique
Le tableau suivant donne le nombre de termes selon le temps:
Instants |
Nombre termes série |
1 |
194 |
10 |
64 |
100 |
22 |
1000 |
8 |
Tableau 1.1-1 : Représentation du nombre de terme en fonction du temps
1.2. Hypothèses simplificatrices#
On considère que le milieu est complètement saturé en eau et on impose une pression de gaz nulle sur tous les nœuds. Le système biphasique se ramène alors à résoudre le problème suivant :
\(\frac{\partial (\varphi {\rho}_{l})}{\partial t}-d({K}_{i}\frac{{\rho}_{l}{k}_{\mathrm{rl}}}{{\mu}_{l}}\nabla {P}_{l})=0\)
Le liquide est incompressible : \({\rho}_{l}=\mathit{cst}\)
La matrice est compressible et la porosité évolue proportionnellement à la pression de liquide : \(\frac{\partial \varphi }{\partial {P}_{l}}={E}_{m}\)
L’équation de conservation de la masse pour le liquide s’écrit donc :
\({\rho}_{l}{E}_{m}\frac{\partial {P}_{l}}{\partial t}-d({K}_{i}\frac{{\rho}_{l}{k}_{rl}}{{\mu}_{l}}\nabla {P}_{l})=0\)
1.3. Géométries#
On considère un barreau \(\mathrm{1D}\) de \(\mathrm{5m}\) de long. Concrètement le domaine maillé fera \([\mathrm{0m},5m]\times [\mathrm{0m};0,05m]\) (dans le cas de la modélisation en triangle, il est important de ne pas avoir de triangles trop «aplatis», le choix de la hauteur du domaine n’est donc pas anodin).
(0,0,05) (5;0,05)
(0,0) (0,5)
1.4. Maillages#
On a testé ce cas sur deux maillages, l’un composé de 100 quadrangles (modélisations A,B et G) et l’autre de 200 triangles (modélisations C et D).
Les modélisations E, F et H consituent une extension \(\mathrm{3D}\) (barreau de section \(1\times 1\) ), le maillage est constitué de 100 hexaèdres.
1.5. Propriétés des matériaux#
On ne donne ici que les propriétés dont la solution dépend, sachant que le fichier de commandes contient d’autres données de matériau qui ne jouent aucun rôle dans la solution du problème traité.
Liquide |
Perméabilité relative Viscosité \((\mathit{kg}.{m}^{-1}.{s}^{-1})\) Module de compressibilité Densité du liquide \((\mathit{kg}/{m}^{3})\) |
\(1\) \(1\) \(0\) \(1\) |
Paramètres homogénéisés |
Perméabilité \(({m}^{2})\) Porosité Emmagasinement Saturation en liquide |
\({10}^{-13}\) \(0,5\) \({10}^{-10}\) \(1\) |
Tableau 1.5-1 : Propriétés des matériaux
1.6. Conditions aux limites et initiales#
Les conditions limites sont les suivantes :
des conditions de Neumann à droite du domaine : \(\frac{\partial {P}_{l}}{\partial x}(t,x=5,y)=0\)
des conditions de Dirichlet sur la partie gauche du domaine : \({P}_{l}(t,x=0,y)=0\)
La pression initiale de liquide est de \({P}_{l}(t=0,x,y)={10}^{4}\mathit{Pa}\) .
1.7. Durée de simulation et pas de temps#
La durée de simulation est de \(100s\) et le nombre de pas de temps est de 100.
Remarque : Par soucis d’information, les tests suivants sont présentés jusqu’à \(1000s\) alors qu’on limite sans Code_Aster les simulations à \(100s\) .