1. Work environment#

In this paragraph, we recall the equations of the HM problem. A dimensioning is then proposed which will be useful for the analysis of errors a posteriori. The reader is advised to have consulted the [2,3,4] documents beforehand for more details on the establishment of the equations given in this section and on their resolution in Code_Aster.

1.1. The continuous model problem#

The working framework is that of linear poroelasticity. We consider a porous medium \(\Omega\), with border \(\Gamma\), saturated by liquid water, with linear elastic mechanical behavior. The behavior of the fluid influences that of the skeleton and vice versa. The spatial dimension is noted \(\text{dim}\). A mechanical balance problem (momentum balance for the skeleton and the fluid) and a hydraulic evolution problem (fluid mass balance) are solved. It should be noted that the treatment of these two problems is completely combined in Code_Aster. The unknowns of the problem are displacement \(u\) and hydraulic pressure \(p\).

Given simulation time \(T\), mechanical balance is formulated

\(\mathrm{-}\mathrm{\nabla }\mathrm{\cdot }{\sigma }^{\text{'}}(u)+b\mathrm{\nabla }p\mathrm{=}f\) in \(\mathrm{[}\text{0,T}\mathrm{]}\mathrm{\times }\Omega\) eq 1.1-1

where we put :math:`f={\rho }_{\text{ref}}{F}^{m}` and

\({\sigma }^{\text{'}}\)

Effective Stress Tensor

\(b\)

Biot coefficient

\({\rho }_{\text{ref}}\)

Homogenized density

\({F}^{m}\)

Gravity force

In the context of isotropic linear elasticity, the effective stress tensor has the following expression

\({\sigma }^{\text{'}}(u)\mathrm{=}{\lambda }_{1}(\mathrm{\nabla }\mathrm{\cdot }u)\text{Id}+{\lambda }_{2}(\mathrm{\nabla }u+\mathrm{\nabla }{u}^{t})\) eq 1.1-2

where \({\lambda }_{1}\) and \({\lambda }_{2}\) are the Lamé coefficients. For its part, the hydraulic equation comes from writing the conservation of fluid mass over time and from Darcy’s law. Darcy’s law provides a proportionality relationship between hydraulic flow \({M}_{\text{lq}}\) and the pressure gradient. More precisely, it is written

\(\frac{{M}_{\text{lq}}}{\rho }\mathrm{=}\kappa (\mathrm{-}\mathrm{\nabla }p+\rho {F}^{m})\) eq 1.1-3

The hydraulic equation is written

\({\mathrm{\partial }}_{t}(\frac{1}{M}p+b\mathrm{\nabla }\mathrm{\cdot }u)\mathrm{-}\mathrm{\nabla }\mathrm{\cdot }(\kappa \mathrm{\nabla }p)\mathrm{=}g\) in \(s\mathrm{[}\text{0,T}\mathrm{]}\mathrm{\times }\Omega\) eq 1.1-4

where \(1\mathrm{/}M\mathrm{=}(b\mathrm{-}\varphi )\mathrm{/}{K}_{s}\) and \(g\mathrm{=}\mathrm{-}\rho \mathrm{\nabla }\mathrm{\cdot }(\kappa {F}^{m})\).

\(\varphi\)

Lagrangian porosity

\(\rho\)

Hydraulic density

\(\kappa\)

Hydraulic conductivity

\({K}_{s}\)

Compressibility of solid grains

 
The equation [:ref:`1.1-1 <1.1-1>`] coupled with [:ref:`1.1-4 <1.1-4>`] is the transient version of the HM problem.

The boundary conditions are of the mixed type of Dirichlet/Neumann in terms of displacement and pressure. We consider a border partition in the form \(\Gamma \mathrm{=}{\Gamma }_{D}^{H}\mathrm{\cup }{\Gamma }_{N}^{H}\mathrm{=}{\Gamma }_{D}^{M}\mathrm{\cup }{\Gamma }_{N}^{M}\).

On \({\Gamma }_{D}^{H}\), conditions are imposed at the Dirichlet limits in hydraulics: \(p\mathrm{=}{p}_{D}\).

On \({\Gamma }_{N}^{H}\), conditions are imposed at the Neumann limits in hydraulics: \({M}_{\text{lq}}\mathrm{\cdot }n\mathrm{=}{M}_{\text{lq},\text{nor}}\).

On \({\Gamma }_{D}^{M}\), conditions are imposed at the Dirichlet limits in mechanics: \(u={u}_{D}\).

On \({\Gamma }_{N}^{M}\), we impose conditions at the Neumann limits in mechanics: \(({\sigma }^{\text{'}}(u)\mathrm{-}\text{bp}{\chi }_{S}\text{Id})\mathrm{\cdot }n\mathrm{=}{\sigma }_{\text{nor}}\), where we noted \({\Gamma }_{S}\mathrm{=}{\Gamma }_{N}^{M}\mathrm{\cap }{\Gamma }_{N}^{H}\), with a characteristic function \({\chi }_{S}\), the portion of the border where constraints and hydraulic flow are imposed.

1.2. Dimensioning#

The model problem involves two unknowns, the movements and the pressure of the liquid, whose orders of magnitude can be very different. We therefore propose a scaling of the equations constituting the model problem, which will serve as a starting point for a posteriori error analysis. Attention is drawn to the fact that this scaling is not used in the numerical resolution stage. The a posteriori error estimates in section 2 will therefore be obtained on the undimensioned problem and then resized for their use in Code_Aster.

For any variable \(a\), we denote the corresponding undimensioned variable \({a}^{\text{*}}\) and \(A\) the corresponding characteristic quantity. For example for hydraulic pressure, we have the relationship

\(p\mathrm{=}{\text{Pp}}^{\text{*}}\),

where \(P\) refers to characteristic pressure. The same symbolism is used for mathematical operators. So, we have: \({\mathrm{\nabla }}^{\text{*}}\mathrm{\cdot }{u}^{\text{*}}\mathrm{=}\frac{L}{U}\mathrm{\nabla }\mathrm{\cdot }u\), \({\varepsilon }^{\text{*}}({u}^{\text{*}})\mathrm{=}\frac{L}{U}\varepsilon (u)\), \({\nabla }^{\text{*}}{p}^{\text{*}}=\frac{L}{P}\nabla p\).

The mechanical equation can be rewritten

\(\mathrm{-}{\mathrm{\nabla }}^{\text{*}}\mathrm{\cdot }{\sigma }^{\text{*}}({u}^{\text{*}})+b{\mathrm{\nabla }}^{\text{*}}{p}^{\text{*}}\mathrm{=}{f}^{\text{*}}\) eq 1.2-1

with \({f}^{\text{*}}=\frac{L}{P}f\). Noting \(\nu\) the Poisson’s ratio, \({\sigma }^{\text{*}}({u}^{\text{*}})\) is an undimensioned version of the stress tensor given by

\({\sigma }^{\text{*}}({u}^{\text{*}})\mathrm{=}\frac{\nu }{(1\mathrm{-}2\nu )(1+2\nu )}({\mathrm{\nabla }}^{\text{*}}\mathrm{\cdot }{u}^{\text{*}})\text{Id}+\frac{1}{1+\nu }{\varepsilon }^{\text{*}}({u}^{\text{*}})\) eq 1.2-2

The following scaling strategy is proposed. The user chooses two quantities: the \(L\) length scale and the \(P\) pressure scale. The \(L\) length scale is determined by the geometry of the model. The pressure scale is determined by the Dirichlet limit conditions imposed on the pressure or possibly by the initial conditions. For the mechanical part, we impose

\(U\mathrm{=}\frac{\text{LP}}{E}\) eq 1.2-3

In addition, the hydraulic equation is written, taking into account [1.2-3],

:math:`{mathrm{partial }}_{t}^{text{*}}(frac{E}{M}{p}^{text{*}}+b{mathrm{nabla }}^{text{*}}mathrm{cdot }{u}^{text{*}})mathrm{-}frac{E}{M}{Delta }^{text{*}}{p}^{text{*}}mathrm{=}{g}^{text{*}}`**eq 1.2-4**

   with :math:`{g}^{\text{*}}\mathrm{=}\frac{E}{M}\frac{{L}^{2}}{\kappa P}g`.

So the dimensioned version of the transient HM problem is:

:math:`mathrm{{}begin{array}{c}mathrm{-}{mathrm{nabla }}^{text{*}}mathrm{cdot }{sigma }^{text{*}}({u}^{text{*}})+b{mathrm{nabla }}^{text{*}}{p}^{text{*}}mathrm{=}{f}^{text{*}}text{dans}mathrm{[}mathrm{0,}{T}^{text{*}}mathrm{]}mathrm{times }{Omega }^{text{*}}\ {mathrm{partial }}_{t}^{text{*}}(frac{E}{M}{p}^{text{*}}+b{mathrm{nabla }}^{text{*}}mathrm{cdot }{u}^{text{*}})mathrm{-}frac{E}{M}{Delta }^{text{*}}{p}^{text{*}}mathrm{=}{g}^{text{*}}text{dans}mathrm{[}mathrm{0,}{T}^{text{*}}mathrm{]}mathrm{times }{Omega }^{text{*}}end{array}`**eq 1.2-5**

1.3. Finite element discretization and notations#

The HM problem is discretized by a finite element method based on a \({({T}_{h})}_{h>0}\) mesh. Recall that the displacements are discretized by polynomials of degree 2 and the pressure by polynomials of degree 1. We note the discreet movements \({u}_{h}\) and the discreet pressure \({p}_{h}\).

\({F}_{h}^{i}\) designates all of the inner faces of the mesh and \({F}_{h}^{\partial }\) refers to all the faces located on the edge of the domain \(\Omega\). For a given \(K\in {T}_{h}\) mesh, we write \({F}_{K}\) the set of faces of \(K\) and we put

\({F}_{K}^{i}={F}_{h}^{i}\cap {F}_{K},{F}_{K}^{\partial }={F}_{h}^{\partial }\cap {F}_{K}\)

Let’s be \(F\mathrm{\in }{F}_{h}^{i}\), that is to say such that there are two meshes \({K}_{1}\) and \({K}_{2}\) in \({T}_{h}\) with \(F={K}_{1}\cap {K}_{2}\). \({n}_{{K}_{1}}\) and \({n}_{{K}_{2}}\) denote the normals outside \({K}_{1}\) and \({K}_{2}\), respectively (see figure 1). We note, for almost everything \(x\in F\),

\(\left[\Phi ({u}_{h})\mathrm{\cdot }n\right](x)\mathrm{=}{(\Phi ({u}_{h}))}_{\mathrm{\mid }{K}_{1}}(x)\mathrm{\cdot }{n}_{{K}_{1}}+{(\Phi ({u}_{h}))}_{\mathrm{\mid }{K}_{2}}(x)\mathrm{\cdot }{n}_{{K}_{2}}\)

the jump of the normal component of \(\Phi ({u}_{h})\) through \(F\).

_images/Object_82.svg

Figure 1: Example of an inner mesh face.

It can be observed that as the vectors \({n}_{{K}_{1}}\) and \({n}_{{K}_{2}}\) are opposite, the quantity above defines this jump well.

For any \(K\) mesh, \({\Delta }_{K}\) refers to the set of cells sharing an edge with \(K\), including mesh \(K\), as illustrated in FIG. 2.

_images/Object_89.svg

Figure 2: Example of mesh \(K\) and macroelement \({\Delta }_{K}\)

For the temporal discretization of the transitory problem, we consider an implicit Euler schema and we note \({\left\{{t}_{i}\right\}}_{i=1}^{N}\) a discretization of the \(\left[\mathrm{0,}T\right]\) interval. We denote \({u}_{h\tau }\) (respectively \({p}_{h\tau }\)) the function continues and refines by pieces in time such that for all \(n\in \left\{\mathrm{0,}\text{.}\text{.}\text{.},N\right\}\), \({u}_{\mathrm{h\tau }}({t}_{n})={u}_{h}^{n}\) (respectively \({p}_{h\tau }({t}_{n})\mathrm{=}{p}_{h}^{n}\)). We also need to consider piecewise constant functions in time, namely \({\pi }^{0}{p}_{\mathit{h\tau }}\) equal to \({p}_{h}^{n}\) over \(\left[{t}_{n-1},{t}_{n}\right]\). For all \(n\in \left\{\mathrm{1,}\text{.}\text{.}\text{.},N\right\}\), use \({\tau }_{n}={t}_{n}-{t}_{n-1}\) and \({I}_{n}=\left[{t}_{n-1},{t}_{n}\right]\).

We also note the discrete « derivatives » of displacements and pressure.

\(\begin{array}{c}{\delta }_{t}{u}_{h}^{n}\mathrm{=}{\tau }_{n}^{\mathrm{-}1}({u}_{h}^{n}\mathrm{-}{u}_{h}^{n\mathrm{-}1})\\ {\delta }_{t}{p}_{h}^{n}\mathrm{=}{\tau }_{n}^{\mathrm{-}1}({p}_{h}^{n}\mathrm{-}{p}_{h}^{n\mathrm{-}1})\end{array}\)

The \(x<y\) notation means that there is a mesh-independent \(c>0\) constant such as \(x\le \text{cy}\). The notation \({\parallel \cdot \parallel }_{V}\) refers to a standard on a \(V\) space. We introduced its version \({\parallel \cdot \parallel }_{V,K}\), localized to mesh \(K\), by \({\mathrm{\parallel }\mathrm{\cdot }\mathrm{\parallel }}_{V}^{2}\mathrm{=}\mathrm{\sum }_{K\mathrm{\in }{T}_{h}}{\mathrm{\parallel }\mathrm{\cdot }\mathrm{\parallel }}_{V,K}^{2}\).