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

(290)\[ \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} \]

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

(291)\[ \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} \]

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

(292)\[ \zeta = \frac{p}{M} + B \nabla\cdot \bm u. \]

Strong and weak formulations

The strong form of linear poroelasticity (mixed \((\bm u\text{-}p_\rf)\) formulation) for constitutive equation (291) 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:

(293)\[ \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} \]

with \(\bm n\) be the unit normal on the boundary and in the conservation of mass (second equation), we replaced \(\zeta\) from constitutive equation (292). The weak form can be derived as:

(294)\[ \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} \]

where we have used Darcy’s law \(\dot{\bm w} = - \varkappa \left(\nabla p_\rf - \bm f^\rf \right)\) (equation (36) for static case).

Generic weak formulation

The generic weak formulation for linear poroelasticity is similar to incompressible elasticity explained in Generic weak form (mixed fields). However, in equation (294) we have time derivatives of displacement and pressure fields and for deriving jacobian we need to add \(\mathrm{shift_v}\) as shown in (224).

For linear poroelasticity described in (294) 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 (224), and we assume body forces \(\bm{f}^\alpha\) are negligible.

For dynamic example 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.

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

(295)\[ \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} \]

being added to the first and second equations of (294), where \(\rho^b\) (the mixture density) is defined in (35) and the jacobian is

(296)\[ \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} \]

where \(\mathrm{shift_a}\) is related to the inverse of time step \(\Delta t\) square (see (224)), 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} \]