# Linear
## Constitutive theory
For the linear poromechanics model, the strain energy density is given by

$$
\psi = \frac{\firstlame_{u}}{2} (\trace \bm{\varepsilon})^2 + \secondlame_{d} \bm{\varepsilon} \tcolon \bm{\varepsilon} - B \, M \, \trace \bm{\varepsilon} \, \zeta + \frac{1}{2} M \zeta^2,
$$

where $\firstlame_{u} = \firstlame_{d} + B^2 M$ is undrained first Lamé parameter while $\firstlame_{d}, \secondlame_{d}$ are Lamé parameters measured in drained condition (i.e., Lamé parameters of the solid skeleton frame), $B, M$ are the Biot effective stress coefficient and Biot modulus defined by

$$
\begin{aligned}
  B &= 1 - \frac{\bulk_{d}}{\bulk_\rs}, \\
  \frac{1}{M} &= \frac{\phi^\rf}{\bulk_\rf} + \frac{B - \phi^\rf}{\bulk_\rs},
\end{aligned}
$$ (biot-coefficient-modulus)

with mixture bulk modulus $\bulk_{d}$ measured in drained condition (solid skeleton bulk modulus, also referred to as $k^{\rm{skel}}$ in other references) and solid and pore fluid bulk moduli $\bulk_\rs, \bulk_\rf$, respectively.

To derive the constitutive law (stress-strain relationship) for the linear poroelasticity model we have

$$
\begin{aligned}
  \bm\sigma(\bm{u}, p) &= \frac{\partial \psi}{\partial \bm{\varepsilon}} = \bm\sigma'(\bm{u}) - B \, p \, \bm{I} = \firstlame_{d} \nabla\cdot \bm{u} \, \bm{I} + 2 \, \secondlame_{d} \, \bm{\varepsilon} - B \, p \, \bm{I}, \\
  p &= \frac{\partial \psi}{\partial \zeta} = M \left(\zeta - B \nabla \cdot \bm{u} \right),
\end{aligned}
$$ (constitutive-linear-poroelastic)

where , $\bm\sigma' = \firstlame_{d} (\trace \bm\varepsilon) \bm{I} + 2 \, \secondlame_{d} \, \bm{\varepsilon}$ is effective stress, $p$ is pore pressure, and the variation of pore fluid content $\zeta$ is

$$
  \zeta = \frac{p}{M} + B \nabla\cdot \bm u.
$$ (fluid-variation2)

## Strong and weak formulations
The strong form of linear poroelasticity (mixed $(\bm u\text{-}p_\rf)$ formulation) for constitutive equation {eq}`constitutive-linear-poroelastic` based on conservation of momentum (for static case) and mass may be stated as follows: Given body forces $\bm f^\rs, \bm f^\rf$ and volumetric injected fluid rate $\dot{\gamma}$, Dirichlet boundaries $\bar{\bm u}, \bar{p}$, applied traction $\bar{\bm t}$ and fluid flux $\bar{s}$ and initial conditions $\bm u_0, p_{\rf,0}$, find the displacement and pressure variables $(\bm u, p_\rf) \in \mathcal{V} \times \mathcal{Q} $ (here $\mathcal{V} = H^1(\Omega), \mathcal{Q} = H^1(\Omega) $ ), such that:

$$
  \begin{aligned}
    -\nabla \cdot \bm \sigma - \bm f^\rs   &= 0, \qquad \text{in $\Omega$} \\
    \frac{\dot{p}_\rf}{M} + \alpha \nabla \cdot \dot{\bm u} + \nabla\cdot \dot{\bm{w}} - \dot{\gamma} &= 0, \qquad \text{in $\Omega$} \\
    \bm{u} &= \bar{\bm u}, \quad    \text{on $\partial \Omega^{D}$} \\
    p_\rf &= \bar{p}_\rf, \quad    \text{on $\partial \Omega^{D}$} \\
    \bm \sigma \cdot \bm n &= \bar{\bm t}, \qquad \text{on $\partial \Omega^{N}$} \\
    \dot{\bm{w}} \cdot \bm n &= \bar{s}, \qquad \text{on $\partial \Omega^{N}$} \\
    \bm u &= \bm u_0, \qquad \text{in $\Omega$} \\
    p_\rf &= p_{\rf,0}, \qquad \text{in $\Omega$}
  \end{aligned}
$$ (poro-linear-strong-form)

with $\bm n$ be the unit normal on the boundary and in the conservation of mass (second equation), we replaced $\zeta$ from constitutive equation {eq}`fluid-variation2`. The weak form can be derived as:

$$
  \begin{aligned}
    \int_{\Omega}{ \nabla \bm{v} \tcolon \bm{\sigma}} \, dv
    - \int_{\partial \Omega}{\bm{v} \cdot \bar{\bm t}} \, da
    - \int_{\Omega}{\bm{v} \cdot \bm f^\rs} \, dv
    &= 0, \quad \forall \bm v \in \mathcal V \\
    \int_{\Omega} q \left(\frac{\dot{p}}{M} + \alpha \nabla \cdot \dot{\bm u} - \dot{\gamma} \right) \, dv + \int_{\Omega} \nabla q \cdot \left( \varkappa \nabla p_\rf - \varkappa \bm f^\rf \right) + \int_{\partial \Omega}{ q \, \bar{s}} \, da & = 0. \quad \forall q \in \mathcal Q
  \end{aligned}
$$ (poro-linear-weak-form)

where we have used Darcy's law $\dot{\bm w} = - \varkappa \left(\nabla p_\rf - \bm f^\rf \right)$ (equation {eq}`generalized-darcy` for static case).

(matrix-free-poromechanics)=
## Generic weak formulation

The generic weak formulation for linear poroelasticity is similar to incompressible elasticity explained in {ref}`matrix-free-mixed-fields`. However, in equation {eq}`poro-linear-weak-form` we have time derivatives of displacement and pressure fields and for deriving jacobian we need to add $\mathrm{shift_v}$ as shown in {eq}`time-dependent-jacobian`.

For linear poroelasticity described in {eq}`poro-linear-weak-form` we have

$$
  \begin{aligned}
    \bm f_0 &= -\bm f^\rs = \bm{0}, \quad \bm f_1 = \bm{\sigma}(\bm u, p_\rf) = \firstlame_{d} \nabla\cdot \bm{u} \, \bm{I} + 2 \, \secondlame_{d} \, \bm{\varepsilon} - \alpha \, p_\rf \, \bm{I}, \\
    \bm g_0 &=\frac{\dot{p}_\rf}{M} + \alpha \nabla \cdot \dot{\bm u} - \dot{\gamma}, \quad \bm g_1 = \varkappa \nabla p_\rf - \varkappa \bm f^\rf, \\
		\diff \bm f_1 &= \bm f_1(\diff \bm{u}, \diff p_\rf), \quad \diff \bm g_0 = \mathrm{shift_v} \left( \frac{\diff p_\rf}{M} + \alpha \nabla \cdot \diff \bm{u} \right), \quad \diff \bm g_1 = \varkappa \nabla \diff p_\rf.
  \end{aligned}
$$

where $\mathrm{shift_v}$ evaluated at time  $\mathrm{shift_v} = \frac{\partial\dot{\phi}}{\partial\phi}|_{\phi_n}$ for generic variable $\phi$ and we have used the linearization of time dependent variables explained {eq}`time-dependent-jacobian`, and we assume body forces $\bm{f}^\alpha$ are negligible.

For [dynamic example](example-dynamic) acceleration is **not zero** and we use the PETSc Timestepper (TS) object to manage timestepping. Specifically, Ratel uses the Generalized-Alpha method for second order systems (`TSALPHA2`) to solve the second order system of equations. More information on the PETSc `TSALPHA2` time stepper can be found in the [PETSc documentation TS](https://petsc.org/release/docs/manualpages/TS/index.html).

For linear poroelastodynamics we consider the solid accelearation in the first equation and in second equation we use {eq}`generalized-darcy` (assume $\ddot{\bm w}=0$). That results in adding following terms:

$$
  \begin{aligned}
    &\int_{\Omega} \bm v \cdot \left(\rho^b \, \ddot{\bm u} \right) \, dv, \\
    &\int_{\Omega} \nabla q \cdot \left(\varkappa \rho^{\fR} \, \ddot{\bm u} \right) \, dv,
  \end{aligned}
$$ (poro-linear-dynamics-residual)

being added to the first and second equations of {eq}`poro-linear-weak-form`, where $\rho^b$ (the mixture density) is defined in {eq}`bulk-density-porous` and the jacobian is

$$
  \begin{aligned}
    &\mathrm{shift_a} \, \int_{\Omega} \bm v \cdot \left(\rho^b \, \diff \bm u \right) \, dv, \\
    &\mathrm{shift_a} \, \int_{\Omega} \nabla q \cdot \left(\varkappa \rho^{\fR} \, \diff \bm u \right) \, dv,
  \end{aligned}
$$ (poro-linear-dynamics-jacobian)

where $\mathrm{shift_a}$ is related to the inverse of time step $\Delta t$ square (see {eq}`time-dependent-jacobian`), such that

$$
  \begin{aligned}
    \bm f_0 &\mathrel{+}= \rho^b\ddot{\bm{u}}, \quad \bm g_1 \mathrel{+}= \varkappa\rho^\fR \ddot{\bm{u}},\\
    \diff \bm f_0 &\mathrel{+}= \mathrm{shift_a}\rho^b \diff \bm u, \quad \diff \bm g_1 \mathrel{+}= \mathrm{shift_a}\varkappa\rho^\fR\, \diff \bm u\, .
  \end{aligned}
$$

## Command-line interface
To enable the linear poromechanics model, use the model option `-model poromechanics-linear` and set the material parameters listed in {ref}`poromechanics-linear-options`. Any parameter without a default option is required.

(poromechanics-linear-options)=
:::{list-table} Linear poromechanics model options
:header-rows: 1
:widths: 30 55 15
* - Option
  - Description
  - Default Value

* - `-model poromechanics-linear`
  - Required to enable the linear poromechanics model.
  -

* - `-lambda_d [real]`
  - First Lame parameter for drained mixture (solid skeleton), $\lambda_d >= 0$
  -

* - `-mu_d [real]`
  - Shear modulus for drained mixture (solid skeleton), $\mu_d >= 0$
  -

* - `-bulk_d [real]`
  - Bulk modulus for drained mixture (solid skeleton), $\bulk_d >= 0$
  -

* - `-bulk_s [real]`
  - Bulk modulus for solid, $\bulk_\rs >= 0$
  -

* - `-bulk_f [real]`
  - Bulk modulus for pore fluid, $\bulk_\rf >= 0$
  -


* - `-rho [real]`
  - Real solid mass density, $\rho^\sR_0 > 0$
  - `1.0`

* - `-rho_fR0 [real]`
  - Real pore fluid mass density, $\rho^\fR_0 > 0$
  - `1.0`

* - `-phi [real]`
  - Porosity, $0 < \phi^\rf < 1$
  -

* - `-mu_f [real]`
  - Pore fluid viscosity, $\mu_\rf >= 0$
  -

* - `-varkappa_0 [real]`
  - Intrinsic permeability, $\varkappa_{0} >= 0$
  -

* - `-alpha [real]`
  - Biot's coefficient, $0 <= \alpha <= 1$
  - $1 - \bulk_d/\bulk_\rs$

* - `-M [real]`
  - Biot's modulus, $M >= 0$
  - $\frac{\bulk_\rs\bulk_\rf}{\bulk_\rf(B - \phi^\rf) + \bulk_\rs\phi^\rf}$
:::

An example using the linear poromechanics model can be run via
```console
$ ./bin/ratel-quasistatic -options_file examples/ymls/ex02-quasistatic-poromechanics-linear-mms-pcfieldsplit-pcjacobi.yml
```