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)

(297)\[ \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

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

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

(299)\[ \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

(300)\[ \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

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

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

(302)\[ \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 (302) 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

(303)\[ \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 (102) 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 (301).

Linearization

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

(304)\[ \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 (110) 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

(305)\[ \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 (303), 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

(306)\[ \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 (299) 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.