1. Reference problem#

The objective of this test case is to compare the solution obtained with the various finite volume diagrams to an analytical solution.

1.1. Analytical solution#

The monophasic unsteady and one-dimensional problem can be written in a general form such as:

\(N\frac{\mathrm{\partial }P}{\mathrm{\partial }t}\mathrm{-}{K}_{\mathit{int.}}\Delta P\mathrm{=}0\)

\(P(t=0)={P}_{0}\)

\(P(t,x=0)=0\)

\(\frac{\partial P}{\partial x}(t,x=L)=0\)

This problem admits an analytical solution obtained by development in Fourier series.

\(P={\sum }_{k=0}^{\infty }\frac{{\mathrm{4P}}_{0}}{(\mathrm{2k}+1)\pi }\mathrm{exp}(\frac{-{K}_{i}}{N}{\omega }_{k}^{2}t)\mathrm{sin}({\omega }_{k}x)\) with \({\omega }_{k}=(k+\frac{1}{2})\frac{\pi }{L}\)

We can reduce this series to a finite number \(K\) of terms, depending on the time calculated.

This number of terms is determined as follows:

Let \({n}_{x}\) be the number of \({x}_{i}\) points where the solution is evaluated at an instant \(t\)

We ask: \({a}_{k}^{i}\mathrm{=}\frac{4}{(\mathrm{2k}+1)\pi }\mathrm{exp}(\frac{\mathrm{-}{K}_{\mathit{int}}}{N}{\omega }_{k}^{2}t)\mathrm{sin}({\omega }_{k}{x}_{i})\)

So the solution can be written as: \(P({x}_{i})={\sum }_{k=0}^{K}{P}_{0}\mathrm{.}{a}_{k}^{i}\)

We choose \(K\) such as: \(\frac{1}{{n}_{x}}\sqrt{{\mathrm{\sum }}_{i\mathrm{=}1}^{\mathit{nx}}{({a}_{k}^{i})}^{2}}<\epsilon\)

In practice, we took \(\epsilon ={10}^{-10}\).

The shapes of the analytical solution at times 1, 10, 100, 1000 are shown in the figure:

_images/100000000000025800000190C21A5A0074EC48A0.png

Figure 1: Representation of analytical solutions

The following table shows the number of terms by time:

Instants

Number terms series

1

194

10

64

100

22

1000

8

Table 1.1-1 : Representation of the number of terms as a function of time

1.2. Simplifying hypotheses#

The medium is considered to be completely saturated with water and zero gas pressure is imposed on all nodes. The biphasic system then comes down to solving the following problem:

\(\frac{\partial (\varphi {\rho }_{l})}{\partial t}-d({K}_{i}\frac{{\rho }_{l}{k}_{\mathrm{rl}}}{{\mu }_{l}}\nabla {P}_{l})=0\)

  • The liquid is incompressible: \({\rho }_{l}\mathrm{=}\mathit{cst}\)

  • The matrix is compressible and the porosity evolves in proportion to the liquid pressure: \(\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }{P}_{l}}\mathrm{=}{E}_{m}\)

The mass conservation equation for liquid is therefore written as:

\({\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. Geometries#

We consider a \(\mathrm{5m}\) long \(\mathrm{1D}\) bar. Concretely, the mesh domain will be \(\mathrm{[}\mathrm{0m}\mathrm{,5}m\mathrm{]}\mathrm{\times }\mathrm{[}\mathrm{0m};\mathrm{0,05}m\mathrm{]}\) (in the case of triangle modeling, it is important not to have triangles that are too « flattened », so the choice of the height of the domain is not trivial).

(0,0,05) (5; 0,05)

(0.0) (0.5)

1.4. Meshes#

This case was tested on two meshes, one composed of 100 quadrangles (models A, B and G) and the other of 200 triangles (models C and D).

The E, F and H models constitute an extension \(\mathrm{3D}\) (bar with a cross section \(1\mathrm{\times }1\)), the mesh consists of 100 hexahedra.

1.5. Material properties#

Only the properties on which the solution depends are given here, knowing that the command file contains other material data that plays no role in solving the problem at hand.

Liquid

Relative permeability Viscosity \((\mathit{kg}\mathrm{.}{m}^{\mathrm{-}1}\mathrm{.}{s}^{\mathrm{-}1})\) Compressibility module Liquid density \((\mathit{kg}\mathrm{/}{m}^{3})\)

\(1\)

\(1\) \(0\) \(1\)

Homogenized parameters

Permeability \(({m}^{2})\) Porosity Storage Liquid saturation

\({10}^{-13}\) \(\mathrm{0,5}\) \({10}^{\mathrm{-}10}\) \(1\)

Table 1.5-1 : Material Properties

1.6. Boundary and initial conditions#

The boundary conditions are as follows:

  • Neumann conditions to the right of the domain: \(\frac{\partial {P}_{l}}{\partial x}(t,x=\mathrm{5,}y)=0\)

  • Dirichlet conditions on the left side of the domain: \({P}_{l}(t,x=\mathrm{0,}y)=0\)

The initial liquid pressure is \({P}_{l}(t\mathrm{=}\mathrm{0,}x,y)\mathrm{=}{10}^{4}\mathit{Pa}\).

1.7. Simulation time and time step#

The simulation duration is \(100s\) and the number of time steps is 100.

Note: For information purposes, the following tests are presented up to \(1000s\) while simulations without Code_Aster are limited to \(100s\) .