[HM] Improvement of the fixed stress splitting approach in the staggered scheme

This MR presents an improvement of the existing implementation with the following physical meaning fixed stress splitting approach (from my understanding):

Balance equation

The balance equations of mass and momentum for the fully saturated porous medium read

S\varrho_\mathrm{f}{\dot p}-\frac{k}{\mu}\nabla \left(\varrho_\mathrm{f}\left(\nabla p-
 \varrho_\mathrm{f} \mathbf{g}  \right) \right) + \alpha \varrho_\mathrm{f} {\dot\varepsilon_v} = 0\\
\nabla\cdot \left( \boldsymbol{\sigma}^\mathrm{eff}(\mathbf{u})
-\alpha p \mathbf{I} \right)
=  \varrho_\mathrm{b} \mathbf{g}.

where \alpha denotes Biot coefficient, S storage coefficient, k permeability, \mu viscosity, \varrho_\mathrm{f} fluid density and \varrho_\mathrm{b} bulk density.

In the staggered scheme for solving HM coupled equations, the fixed-stress splitting is employed to enhance the convergence. The fixed stress splitting is based on the the volumetric total stress rate definition the hydro-mechanics:

 \dot{\sigma}_v=K_b ({\dot \varepsilon}_v-\dot{\varepsilon}^{ne}_v)- \alpha\dot {p}.

with K_b the drained bulk modulus of porous medium, and \dot{\varepsilon}^{ne}_v the volumetric non-elastic strain rate.

Fixed stress rate over coupling iteration

As the first option, we consider to fix the volumetric total stress rate over coupling iteration. This means

	\dot{\sigma}_v^{n, k} = \dot{\sigma}_v^{n, k-1}.

with n the time step index, k the coupling iteration index, and \dot{()}^{n, k} = \left(()^{n, k} - ()^{n}\right)/dt.

This gives the volumetric strain rate of the current time step as

\dot{\varepsilon_v}^{n, k} \approx  \dot{\epsilon}_v^{n, k-1} +
\dfrac{ \alpha}{K_b}   (\dot {p}^{n, k}-\dot {p}^{n, k-1})
	+(\dot{\epsilon}^{ne}_v|^{n, k}-
\dot{\epsilon}^{ne}_v|^{n, k-1}) .

Practically, we can set \dot{\epsilon}^{ne}_v|^{n, k}-\dot{\epsilon}^{ne}_v|^{n, k-1} = 0.

Under that volumetric strain rate approximation, the mass balance equation for coupling iteration k at time step n becomes

\varrho_\mathrm{f}(S +\dfrac{ \alpha^2}{K_b} ) {\dot p}^{n, k}-
\frac{k}{\mu}\nabla \left(\varrho_\mathrm{f}\left(\nabla p^{n. k}-
 \varrho_\mathrm{f} \mathbf{g}  \right) \right) +
  \varrho_\mathrm{f} \left(\alpha {\dot\varepsilon_v}^{n,k-1}-
 \dfrac{ \alpha^2}{K_b} {\dot p}^{n, k-1}\right)= 0.

Denoting \frac{ \alpha^2}{K_b} as \beta_\mathrm{FS}, the above equation turns into

\varrho_\mathrm{f}(S +\beta_\mathrm{FS} ) {\dot p}^{n, k}-
\frac{k}{\mu}\nabla \left(\varrho_\mathrm{f}\left(\nabla p^{n. k}-
 \varrho_\mathrm{f} \mathbf{g}  \right) \right) +
  \varrho_\mathrm{f} \left(\alpha {\dot\varepsilon_v}^{n,k-1}-
 \beta_\mathrm{FS} {\dot p}^{n, k-1}\right)= 0.

One can see from the above equation that \beta_\mathrm{FS} can be any scalar number once the coupling iteration converges, i.e. {\dot p}^{n, k}\approx{\dot p}^{n, k-1}. Therefore, one can choose an arbitrary value for \beta_\mathrm{FS} if it can make the coupling iteration converge. The range of \beta_\mathrm{FS} is an interesting study topic. For code implementation, we introduce a factor to the physically meaningful parameter \frac{\alpha^2}{K_\mathrm{b}} as:

\beta_\mathrm{FS}=\beta_\mathrm{FS}^{opt}\dfrac{\alpha^2}{K_\mathrm{b}},

where \beta_\mathrm{FS}^{opt} is treated as an input parameter.

The existing implementation is equivalent to the approach of a fixed stress rate over coupling iteration.

Fixed stress rate over time step

We assume that the volumetric stress rate of the current time step is the same as that of the previous time step:

  \dot{\sigma}_v^{n} = \dot{\sigma}_v^{n-1}.

That means the current volumetric strain rate is

\dot{\varepsilon_v}^{n} \approx  \dot{\epsilon}_v^{n-1} +
\dfrac{ \alpha}{K_b}   (\dot {p}^{n}-\dot {p}^{n-1})
	+(\dot{\epsilon}^{ne}_v|^{n}-
\dot{\epsilon}^{ne}_v|^{n-1}) .

Consequently, and similarly to the fixed stress over coupling iteration, the mass balance equation at time step n becomes

\varrho_\mathrm{f}(S +\beta_\mathrm{FS} ) {\dot p}^{n}-
\frac{k}{\mu}\nabla \left(\varrho_\mathrm{f}\left(\nabla p^{n}-
 \varrho_\mathrm{f} \mathbf{g}  \right) \right) +
  \varrho_\mathrm{f} \left(\alpha {\dot\varepsilon_v}^{n-1}-
 \beta_\mathrm{FS} {\dot p}^{n-1}\right)= 0.

In that sense, only one coupling iteration is needed, and the solution accuracy is dependent on the time step size. The approach of a fixed stress rate over the time step enables the staggered scheme to solve most of the HM problems, especially those with small storage.

Source code changes

1. refactored the stress splitting portion.
2. added an option to fix fixed stress rate either over coupling iteration or over time step. Associated to this change, the input syntax of coupling_term (optional) is changed to

    <coupling_scheme>
        <type>staggered</type>
        <fixed_stress_stabilization_parameter> <!--optional-->
          1.0 </fixed_stress_stabilization_parameter>
        <!--optional--> 
        <fixed_stress_over_time_step>true</fixed_stress_over_time_step>
    </coupling_scheme>
   

3. re-written the documentations.
4. enabled restarting computation with fixed stress splitting.
5. added a new benchmark, Consolidation Benchmark with staggered scheme, and its documentation. The new benchmark has analytical solution (see the solution comparisons in the following figures).

1. Feature description was added to the changelog
2. Tests covering your feature were added?
3. Any new feature or behavior change was documented?