Finite strain¶

Constitutive theory¶

For finite strain poromechanics, by writing Clausius-Duhem inequality (dissipation inequality) we can derive constitutive equations for solid and fluid phases as [ICR24] (see also (61) and discussion thereof)

(299)¶\[ \bm{S}^\rs_{E(\rs)} = 2\frac{\partial(\rho_{0(\rs)}^\rs\psi^\rs)}{ \partial \bm{C}_\rs}, \quad p_\rf = (\rho^{\fR})^2\frac{\partial \psi^\rf}{\partial \rho^{\fR}}, \]

where $\bm{S}^\rs_{E(\rs)}$ is effective second Piola-Kirchoff stress tensor, $\psi^\rs$ is strain energy per unit mass for the solid skeleton which is defined in previous sections for neo-Hookean model.

Helmholtz free energy per unit mass for fluid phase $\psi^\rf$ is given by

(300)¶\[ \psi^\rf = \frac{1}{2}\phi^\rf \bulk_\rf \left(\log \rho^{\fR}\right)^2, \]

which using second equation of (299), we arrive at constitutive equation for the real mass density of the pore fluid as

(301)¶\[ \rho^{\fR} = \rho^{\fR}_0 \exp\left( \frac{p_\rf - p_{\rf,0}}{\bulk_\rf}\right), \]

where $\rho^{\fR}_0, p_{\rf,0}$ are the initial real mass density and pore fluid pressure at time $t = 0$.

The total stress for the mixture for neo-Hookean model (83) (dropping now $(\cdot)_\rs$ subscripts where appropriate) is

(302)¶\[ \begin{aligned} \bm{S} &= J \bm{F}^{-1} \bm{\sigma} \bm{F}^{-T} = \bm{S}^\rs_E - J p_\rf \bm{C}^{-1}, \\ \bm{S}^\rs_E &= \firstlame_d J V' \bm{C}^{-1} + \secondlame_d \left( \bm{I} - \bm{C}^{-1} \right) = \frac{\firstlame_d}{2} \left( J^2 - 1 \right)\bm{C}^{-1} + \secondlame_d \left( \bm{I} - \bm{C}^{-1} \right), \end{aligned} \]

and using (25) we can write symmetric Kirchhoff stress tensor as

(303)¶\[ \begin{aligned} \bm{\tau} &= \bm{\tau}^\rs_E - J p_\rf \bm{I}, \\ \bm{\tau}^\rs_E &= \firstlame_d J V' \bm{I} + \secondlame_d \left( \bm{b} - \bm{I} \right) = \frac{\firstlame_d}{2} \left( J^2 - 1 \right)\bm{I} + \secondlame_d \left( \bm{b} - \bm{I} \right). \end{aligned} \]

Weak enforcement of solid incompressibility via strain-energy potential¶

For larger deformations, it is possible in the numerical sense to compute $J \leq \phi^{\rm{s}}_0$, which violates the solid incompressibility assumption (refer to (46)). [EE99] attempted to rectify this by modifying the neo-Hookean potential

(304)¶\[ \begin{aligned} V(J) := ({\rm{log}}(J))^2 \quad \text{or} \quad V(J) := \frac{1}{4}(J^2 - 1 - 2{\rm{log}}J) \end{aligned} \]

(305)¶\[ \begin{aligned} V(J) := (1 - \phi^\rs_0)^2\left(\frac{J - 1}{1 - \phi^\rs_0} - {\rm{log}}\left[\frac{J - \phi^\rs_0}{1 - \phi^\rs_0}\right]\right)\, . \end{aligned} \]

Then, the initial solid skeleton effective stress is

(306)¶\[ \begin{aligned} \bm{S}_E^\rs &= \firstlame(1 - \phi_0^\rs)^2 \left(\frac{J}{1 - \phi_0^\rs} - \frac{J}{J - \phi_0^\rs} \right)\bm{C}^{-1} + \secondlame \left( \bm{I} - \bm{C}^{-1} \right)\, , \end{aligned} \]

such that the solid skeleton effective stress in the current configuration is

(307)¶\[ \begin{aligned} \bm{\tau}_E^\rs &= \firstlame(1 - \phi_0^\rs)^2 \left(\frac{J}{1 - \phi_0^\rs} - \frac{J}{J - \phi_0^\rs} \right)\bm{I} + \secondlame \left(\bm{b} - \bm{I}\right)\, . \end{aligned} \]

Deformation-dependent porosity functions¶

Recall generalized Darcy’s law from (47):

\[ \dot{\bm{w}} = \phi^\rf\tilde{\bm{v}}_\rf = -\varkappa(\nabla_x p_\rf + \rho^\fR(\ddot{\bm{u}}_\rf - \bm{f}^\rf))\, , \]

where

\[ \varkappa = \frac{\varkappa_\circ}{\mu_\rf}\frac{\mathcal{F}(\phi^\rf)}{\mathcal{F}(\phi^\rf_0)}\, . \]

In addition to Kozeny-Carman (48), another suitable function was developed for soft tissues by [LMR81]:

\[ \mathcal{F}(\phi^\rf) = \exp[\kappa(J - 1)] = \exp\left[\kappa\left(\frac{\phi^\rs_0}{1 - \phi^\rf} - 1\right)\right]\, . \]

However, both the Kozeny-Carman and [LMR81] functional forms do not strictly satisfy the solid incompressibility constraint, which mandates a bound on the Jacobian of deformation: $\phi^\rs_0 < J < \infty$. This in turn mandates the behavior of $\mathcal{F}(\phi^\rf)$:

(311)¶\[ \begin{aligned} \mathcal{F}(\phi^\rf) \rightarrow \begin{Bmatrix}0 \\ 1\end{Bmatrix} \text{ if } \begin{cases} \phi^\rf \rightarrow 0 \Leftrightarrow \phi^\rs \rightarrow 1 \Leftrightarrow J \rightarrow \phi_0^\rs\\ \phi^\rf \rightarrow 1 \Leftrightarrow \phi^\rs \rightarrow 0 \Leftrightarrow J \rightarrow \infty\end{cases}\, , \end{aligned} \]

in other words, for vanishing porosity we should obtain a solid-like material and for full porosity (fully expanded solid) the solid-fluid interactions should be negligible beyond the viscous drag forces induced by $\varkappa_\circ/\mu_\rf$.

As such, [Eip98] proposed a deformation dependent permeability that respects the restriction of materially incompressible solid constituent:

\[ \begin{aligned} \mathcal{F}(\phi^\rf) = (\phi^\rf)^\kappa \end{aligned} \]

for $\kappa > 0$, however, such a function is ill-suited for rapidly expanding materials (see A comparison, adapted from Markert2005, of different deformation-dependent hydraulic conductivities for \kappa = 3.0 (a) showing the instability of the Kozeny-Carman relation for high compression of the drained mixutre (solid skeleton) and (b) a zoomed-in version of (a) showing the inadequacy of the Eipper1998 model for high expansion of the drained mixture (solid skeleton).). In light of this, [Mar05] proposed

\[ \begin{aligned} \mathcal{F}(\phi^\rf) = \Bigg(\frac{\phi^\rf}{1 - \phi^\rf}\Bigg)^\kappa \, , \end{aligned} \]

such that

\[ \begin{aligned} \varkappa = \frac{\varkappa_{\circ}}{\mu_\rf}\Bigg(\frac{J - \phi_0^\rs}{1 - \phi_0^\rs}\Bigg)^\kappa\, . \end{aligned} \]

../../../../../_images/permeability-functionals.png — Fig. 6 A comparison, adapted from [Mar05], of different deformation-dependent hydraulic conductivities for $\kappa = 3.0$ (a) showing the instability of the Kozeny-Carman relation for high compression of the drained mixutre (solid skeleton) and (b) a zoomed-in version of (a) showing the inadequacy of the [Eip98] model for high expansion of the drained mixture (solid skeleton).¶

Strong and weak formulations¶

To derive the strong form for finite strain poromechanics, we use the balance of momentum of the mixture (54)$_1$ (under the assumption that $\ddot{\bm{u}} = \bm{0}$) and balance of mass of the mixture (44) as

(312)¶\[ \begin{aligned} -\nabla_X \cdot \bm{P} - \rho_0 \bm g &= 0, \qquad \text{in $\Omega_0$} \\ \frac{\phi^\rf}{\bulk_\rf}\frac{D p_\rf}{D t} + \nabla_x\cdot \dot{\bm{u}} + \nabla_x\cdot \dot{\bm w} + \frac{\dot{\bm w}}{\bulk_\rf}\cdot \nabla_x p_\rf &= 0, \qquad \text{in $\Omega$} \\ \bm{u} &= \bar{\bm u}, \quad \text{on $\partial \Omega^{D}_0$} \\ p_\rf &= \bar{p}_\rf, \quad \text{on $\partial \Omega^{D}$} \\ \bm{P} \cdot \bm N &= \bar{\bm t}, \qquad \text{on $\partial \Omega^{N}_0$} \\ \dot{\bm{w}} \cdot \bm n &= \bar{s}, \qquad \text{on $\partial \Omega^{N}$} \\ \bm u &= \bm u_0, \qquad \text{in $\Omega_0$} \\ p_\rf &= p_{\rf,0}, \qquad \text{in $\Omega$} \end{aligned} \]

where we have dropped subscript $(\cdot)_\rs$ by convention, $\bm N$ and $\bm n$ are the unit normal vectors on the face in initial and current configurations, $_X$ and $_x$ in $\nabla_X, \nabla_x$ indicate that the gradient are calculated with respect to the initial/current configurations. Multiplying the strong form (312) by test functions $(\bm{v}, q)$ and integrate by parts, we can obtain the weak form for finite strain hyperporoelasticity as: find $(\bm{u}, p_\rf) \in \mathcal{V} \times \mathcal{Q} \subset H^1 \left( \Omega_0 \right) \times H^1 \left( \Omega_0 \right)$ such that

(313)¶\[ \begin{aligned} r_u(\bm u, p_\rf) &\coloneqq \int_{\Omega_0}{\nabla_x \bm{v} \tcolon \bm{\tau}} \, dV - \int_{\Omega_0}{\bm{v} \cdot \rho_0 \bm{g}} \, dV - \int_{\partial \Omega_0}{\bm{v} \cdot \bar{\bm t}} \, dA = 0, \quad \forall \bm{v} \in \mathcal{V}, \\ r_{p_\rf}(\bm u, p_\rf) &\coloneqq \int_{\Omega_0} q \cdot \left( \frac{\phi^\rf}{\bulk_\rf}\dot{p}_\rf + \nabla_x\cdot \dot{\bm{u}} - \frac{\varkappa}{\bulk_\rf}\nabla_x p_\rf \cdot \nabla_x p_\rf \right) J \, dV \\ &+\int_{\Omega_0} \nabla_x q \cdot \varkappa \nabla_x p_\rf \, J dV + \int_{\partial \Omega}{ q \, \bar{s}} \, da = 0, \quad \forall q\in \mathcal{Q}, \\ \end{aligned} \]

where we have used the push forward (103) to write the first equation with respect to $_x$ and $\bm \tau$, and replaced $\dot{\bm w} = -\varkappa \nabla_x p_\rf$ by assuming $\ddot{\bm u} = \bm{0}$ and $\bm f^\rf = \bm{0}$ in (49), the latter of which implies that $\bm{g} = \bm{0}$. The definition for $\bm{\tau}$ for the standard neo-Hookean model for the drained mixture (solid skeleton) response is given in (303).

Linearization¶

As explained in Linearization: initial configuration, we need the Jacobian form of the (313) which is: find $(\diff \bm{u}, \diff p_\rf) \in \mathcal{V} \times \mathcal{Q}$ such that

(314)¶\[ \begin{aligned} &\int_{\Omega_0} \nabla_x \bm{v} \tcolon \left( \diff \bm{\tau} -\bm{\tau}(\nabla_x \diff\bm{u})^T \right) dV = -r_u, \quad \forall \bm{v} \in \mathcal{V}, \\ &\int_{\Omega_0} q \, \left( \diff L J + L \diff J \right) dV + \int_{\Omega_0} \nabla_x q \cdot (\diff K - K\nabla_x\diff\bm{u})\, dV = -r_{p_\rf}, \quad \forall q \in \mathcal{Q} \end{aligned} \]

where

\[ \begin{aligned} \diff \bm{\tau} -\bm{\tau}(\nabla_x \diff\bm{u})^T &= \diff \bm{\tau}_E^\rs -\bm{\tau}_E^\rs(\nabla_x \diff\bm{u})^T - \left(\diff J p_\rf + J \diff p_\rf \right) \bm I + J p_\rf (\nabla_x \diff\bm{u})^T, \\ &= (\nabla_x \diff\bm{u}) \bm{\tau} + \bm{F} \diff \bm{S}_E^\rs \bm{F}^T - \left(\diff J p_\rf + J \diff p_\rf \right) \bm I + 2 J p_\rf \diff\bm\epsilon, \\ L &= \frac{\phi^\rf}{\bulk_\rf}\dot{p}_\rf + \nabla_x\cdot \dot{\bm{u}} - \frac{\varkappa}{\bulk_\rf}\nabla_x p \cdot \nabla_x p, \\ \diff L &= \mathrm{shift_v} \left(\frac{\phi^\rf}{\bulk_\rf}\diff p_\rf + \diff \, \left(\nabla_x\cdot \bm{u} \right) \right) - \frac{\diff \varkappa}{\bulk_\rf}\nabla_x p_\rf \cdot \nabla_x p_\rf - 2 \frac{\varkappa}{\bulk_\rf}\nabla_x p_\rf \cdot \diff \, (\nabla_x p_\rf), \\ K &= \varkappa \nabla_x p_\rf J\, ,\\ \diff K &= \diff \varkappa \nabla_x p_\rf \, J + \varkappa \diff \, (\nabla_x p_\rf) \, J + \varkappa \nabla_x p_\rf \diff J, \end{aligned} \]

with given $\diff \bm{\tau}_E^\rs -\bm{\tau}_E^\rs(\nabla_x \diff\bm{u})^T$ in (111) for example with the standard neo-Hookean model and

\[ \begin{aligned} \diff J &= J \, \trace (\nabla_x \diff\bm u) = J \, \trace \diff\bm\epsilon, \\ \diff \varkappa &= \frac{\varkappa_{\circ}}{\mu_\rf \mathcal{F}(\phi^\rf_0)} \diff \mathcal{F}(\phi^\rf) = \frac{\varkappa_{\circ}}{\mu_\rf \mathcal{F}(\phi^\rf_0)} \frac{\partial \mathcal{F}(\phi^\rf)}{\partial \phi^\rf} \diff \phi^\rf = \varkappa\frac{\partial\mathcal{F}(\phi^\rf)}{\partial\phi^\rf} \mathcal{F}(\phi^\rf)^{-1} \diff \phi^\rf, \\ \diff \phi^\rf &= \diff \left(1 - \phi^\rs \right) = \diff \left(1 - \frac{\phi^\rs_0}{J} \right) = \frac{\phi^\rs_0}{J^2} \diff J = \phi^\rs \trace \diff\bm\epsilon, \\ \diff \, (\nabla_x p_\rf) &= \diff \, \left( \bm{F}^{-T} \nabla_X p_\rf \right) = \nabla_x \diff p_\rf - (\nabla_x \diff\bm{u})^T \nabla_x p_\rf, \\ \diff \, \left(\nabla_x \cdot \bm{u} \right) &= \diff \, \trace \left( \nabla_x \bm{u}\right) = \diff \, \trace \left( \nabla_X \bm{u} \bm{F}^{-1}\right) = \trace \left( \nabla_x \diff \bm{u} - \nabla_x \bm{u} \nabla_x \diff \bm{u}\right), \\ &= \nabla_x \cdot \diff\bm{u} - \nabla_x \bm{u} \tcolon (\nabla_x \diff\bm{u})^T, \end{aligned} \]

where $\mathcal{F}(\phi^\rf)$ is defined, for example, in (48).

Generic weak formulation¶

For overview on generic weak form implementation for the poromechanics model, see Generic weak formulation. Here we define

\[ \begin{aligned} \bm f_0 &= -\rho_0 \bm g = \bm{0}, \quad \bm f_1 = \bm{\tau} = \bm{\tau}_E^\rs - J p_\rf \bm{I}, \\ \bm g_0 &= L \, J = \left( \frac{\phi^\rf}{\bulk_\rf}\dot{p}_\rf + \nabla_x\cdot \dot{\bm{u}} - \frac{\varkappa}{\bulk_\rf}\nabla_x p_\rf \cdot \nabla_x p_\rf \right) J, \quad \bm g_1 = K = \varkappa \nabla_x p_\rf \, J, \\ \diff \bm f_1 &= \diff \bm{\tau} -\bm{\tau}(\nabla_x \diff\bm{u})^T, \quad \diff \bm g_0 = \diff L J + L \diff J, \\ \diff \bm g_1 &= \diff K - K\nabla_x\diff \bm{u} = \diff \varkappa \nabla_x p_\rf \, J + \varkappa \diff \, (\nabla_x p_\rf) \, J + \varkappa \nabla_x p_\rf \diff J - \varkappa \nabla_x p_\rf J\nabla_x \diff \bm{u} . \end{aligned} \]

Note that, technically speaking, $\bm f_0$ needs to be linearized under the assmption of $\bm{g} \neq \bm{0}$ (and by extension $\bm{f}^\rf \neq \bm{0}$) in the case of evolving volume fractions and/or compressible pore fluid. However, at this time, we continue to assume the body force acting on the mixture to be zero.

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.

We add

(315)¶\[ \begin{aligned} &\int_{\Omega_0} \bm v \cdot \left(\rho^b_0 \, \ddot{\bm u} \right) \, dV = \int_{\Omega_0} \bm v \cdot \left(J \rho^b \, \ddot{\bm u} \right) \, dV, \\ &\int_{\Omega_0} q \left(-\varkappa \, \ddot{\bm u} \cdot \nabla_x p_\rf\right) \frac{J \, \rho^{\fR}}{\bulk_\rf} dV + \int_{\Omega_0} \nabla_x q \cdot \left(\varkappa \rho^{\fR} J \, \ddot{\bm u}\right) \, dV, \end{aligned} \]

to the first and second equations of (313), such that

\[ \begin{aligned} \bm f_0 &\mathrel{+}= J\rho^b\ddot{\bm{u}}\, ,\\ \bm g_0 &\mathrel{+}= -\frac{\varkappa}{\bulk_\rf}\nabla_x p_\rf\cdot(\rho^\fR\ddot{\bm{u}})J\,\\ \bm g_1 &\mathrel{+}= \varkappa\rho^\fR\ddot{\bm{u}} J\, , \end{aligned} \]

and the jacobian is

(316)¶\[ \begin{aligned} &\int_{\Omega_0} \bm v \cdot \left(\diff J \rho^b \, \ddot{\bm u} + J \diff \rho^b \, \ddot{\bm u} + J \rho^b \, \mathrm{shift_a} \, \diff \bm u \right) \, dV, \\ &\int_{\Omega_0} q \left[-\diff \varkappa \, \ddot{\bm u} \cdot \nabla_x p_\rf - \varkappa \, \mathrm{shift_a} \diff \bm u \cdot \nabla_x p_\rf - \varkappa \, \ddot{\bm u} \cdot \diff \left(\nabla_x p_\rf \right) \right] \frac{J \, \rho^{\fR}}{\bulk_\rf} dV \\ &\int_{\Omega_0} q \left(-\varkappa \, \ddot{\bm u} \cdot \nabla_x p_\rf\right) \frac{\left(\diff J \rho^{\fR} + J \, \diff \rho^{\fR}\right)}{\bulk_\rf} \, dV + \\ &\int_{\Omega_0} \nabla_x q \cdot \left[ \left(\diff \varkappa \rho^{\fR} J + \varkappa \diff \rho^{\fR} J + \varkappa \rho^{\fR} \diff J \right) \ddot{\bm u} + \varkappa \rho^{\fR} J \, \mathrm{shift_a}\diff \bm u - \varkappa \rho^{\fR} J (\nabla_x \diff \bm u) \ddot{\bm u} \right] \, dV, \end{aligned} \]

such that,

\[ \begin{aligned} \diff \bm f_0 &\mathrel{+}= \diff J \rho^b \, \ddot{\bm u} + J \diff \rho^b \, \ddot{\bm u} + J \rho^b \, \mathrm{shift_a} \, \diff \bm u\, ,\\ \diff \bm g_0 &\mathrel{+}= \left[-\diff \varkappa \, \ddot{\bm u} \cdot \nabla_x p_\rf - \varkappa \, \mathrm{shift_a} \diff \bm u \cdot \nabla_x p_\rf - \varkappa \, \ddot{\bm u} \cdot \diff \left(\nabla_x p_\rf \right) \right] \frac{J \, \rho^{\fR}}{\bulk_\rf}\\ &\phantom{\mathrel{+}=} + \left(-\varkappa \, \ddot{\bm u} \cdot \nabla_x p_\rf\right) \frac{\left(\diff J \rho^{\fR} + J \, \diff \rho^{\fR}\right)}{\bulk_\rf}, \\ \diff \bm g_1 &\mathrel{+}= \left(\diff \varkappa \rho^{\fR} J + \varkappa \diff \rho^{\fR} J + \varkappa \rho^{\fR} \diff J \right) \ddot{\bm u} + \varkappa \rho^{\fR} J \, \mathrm{shift_a}\diff \bm u - \varkappa \rho^{\fR} J (\nabla_x \diff \bm u) \ddot{\bm u}\, , \end{aligned} \]

where using (301) we have

\[ \begin{aligned} \diff \rho^b &= \left(\rho^{\fR}-\rho^{\sR} \right) \diff \phi^\rf + \phi^\rf \diff \rho^{\fR}, \\ \diff \rho^{\fR} &= \frac{\rho^{\fR}_0 \, \diff p_\rf}{\bulk_\rf} \exp\left( \frac{p_\rf - p_{\rf,0}}{\bulk_\rf}\right) = \frac{\diff p_\rf}{\bulk_\rf}\rho^{\fR}. \end{aligned} \]

Note that unlike linear case where the porosity and density are constant, in finite strain we have assumed that the real mass density of the pore fluid is a function of pore fluid pressure (real mass density of solid is assumed constant, see (45)). Therefore, we use the mixture density in the current configuration $\rho^b$ in residual and jacobian to account for its variation with respect to pore fluid pressure.

Pore fluid pressure stabilization¶

For stabilizing equal order elements, we follow the approach of [TZ06], which is based on the method of [BPitkaranta84]. The stabilization term acts like an opposing pore fluid pressure flux in the weak formulation of the balance of mass of the mixture, such that

\[ \begin{aligned} r_{p_\rf}(\bm u, p_\rf) &\coloneqq \int_{\Omega_0} q \cdot \left( \frac{\phi^\rf}{\bulk_\rf}\dot{p}_\rf + \nabla_x\cdot \dot{\bm{u}} - \frac{\varkappa}{\bulk_\rf}\nabla_x p_\rf \cdot \nabla_x p_\rf \right) J \, dV \\ &+\int_{\Omega_0} \nabla_x q \cdot \left(\varkappa \nabla_x p_\rf + {\color{red}\alpha^\text{stab}\nabla_x\dot{p}_\rf}\right)\, J dV + \int_{\partial \Omega}{ q \, \bar{s}} \, da = 0, \quad \forall q\in \mathcal{Q}, \\ \end{aligned} \]

where the stabilization term is shown in red (and may also be enabled for dynamic cases). While [TZ06] have been able to relate the pressure stabilization parameter $\alpha^\text{stab}$ (units m$^3$-s$^2$/kg) to material geometry and simulation time-step, such an approach would be difficult to derive for compressible pore fluid and large deformations. Therefore, we choose $\alpha^\text{stab}$ on an ad-hoc basis. For incompressible pore fluid (i.e., $\bulk_\rf \sim $ GPa), choosing larger values $\mathcal{O}(10^{-6})$ tends to provide decent stability as compared to smaller values. For more compressible pore fluid (i.e., lower $\bulk_\rf \sim$ kPa) choosing smaller values $\mathcal{O}(10^{-10})$ tends to provide decent stability as compared to larger values.

For generic weak formulation, two additional terms are present for $\bm{g}_1$ and $\diff\bm{g}_1$, such that

\[ \begin{aligned} \bm g_1 &\mathrel{+}= J\alpha^\text{stab}\nabla_x\dot{p}_\rf\, ,\\ \diff \bm g_1 &\mathrel{+}= \alpha^\text{stab}\left(\diff J \, \nabla_x\dot{p}_\rf + J \mathrm{shift_v} \diff(\nabla_x p_\rf) - J\nabla_x\dot{p}_\rf(\nabla_x \diff \bm u)\right)\, . \end{aligned} \]

Note that pressure stabilization can be added to unequal order elements as well. In [IRC23] it was shown that under conditions of 1-D uniaxial strain, pressure stabilization improved computational performance for Hermite cubic-Lagrange linear displacement-pressure elements.

Command-line interface¶

To enable the finite strain, neo-Hookean-based poromechanics model, use the model option -model poromechanics-neo-hookean and set the material parameters listed in Finite strain poromechanics model options. Any parameter without a default option is required.

Table 34 Finite strain poromechanics model options¶
Option	Description	Default Value
`-model poromechanics-neo-hookean-current`	Required to enable the finite strain 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_stab [real]`	Pore fluid pressure stabilization, $\alpha^\text{stab} > 0$	`0.0` (default: disabled)
`-kappa [real]`	Hyperbolic porosity model exponent, $\kappa >= 0$	`-1.0` (default: Kozeny-Carman model)

An example using the finite strain poromechanics model can be run via

$ ./bin/ratel-dynamic -options_file examples/ymls/ex03-dynamic-poromechanics-neo-hookean-current-pstab-pcfieldsplit-lu.yml