PoroMechanics is the theory that describes the coupling of porous flow with solid mechanics.

Notation

The notation used in this article is defined below. For more explanation, see the documentation for standard TensorMechanics and the PorousFlow module (Richards equation).

  • $$$x_{i}$$$, coordinate. TensorMechanics is set up for 3D simulations, although the material in this article holds for any number of dimensions. A derivative with respect to $$$x_{i}$$$ is denoted by $$$\nabla_{i}$$$.

  • $$$t$$$, time. A dot denotes a derivative with respect to $$$t$$$.

  • $$$u_{i}$$$, the components of displacement.

  • $$$v_{i} = {\mathrm d}u_{i}/{\mathrm d} t$$$, the components of velocity.

  • $$$\epsilon_{ij}$$$, components of strain: $$$\epsilon_{ij} = (\nabla_{i}u_{j} + \nabla_{j}u_{i})/2$$$. The equations below may be generalised to large strains in the usual TensorMechanics way (see the TensorMechanics articles).

  • $$$V$$$, an arbitrary volume of the porous material

  • $$$V_{\mathrm{f}}$$$, the porevolume within volume $$$V$$$. If the porous material is fully saturated then $$$V_{f}$$$ is the volume of fluid within volume $$$V$$$.

  • $$$\phi = V_{\mathrm{f}}/V$$$, the porosity of the material. This is a Material property, and it evolves with changes in stress and porepressure.

  • $$$S$$$, the fluid saturation. For a single-phase fluid, there is just one saturation. However, for a multi-phase fluid each phase has a saturation, and $$$\sum_{\mathrm{phases}}S_{\mathrm{phases}}=1$$$. In the single-phase case, the term "fully saturated" means $$$S=1$$$, and this is a common situation in poromechanics.

  • $$$P_{\mathrm{f}}$$$, the poro-mechanics porepressure. For a single-phase fluid, $$$P_{\mathrm{f}}$$$ is just the fluid's porepressure. However, for multi-phase fluid, it is standard to use $$$P_{\mathrm{f}} = \sum_{\mathrm{phases}}S_{\mathrm{phase}}P_{\mathrm{phase}}$$$, where $$$S$$$ and $$$P$$$ are the saturations and pressures of the phase.

  • $$$\rho$$$, the mass density of a fluid. For the single-phase situation, there will be just one $$$\rho$$$, while for the multi-phase situation each phase has a separate density. The mass density is usually a function of the fluid's porepressure.

  • $$$K_{\mathrm{f}}$$$, the bulk modulus of the fluid. It is common in poromechanics with a single-phase fluid to assume the fluid has a constant bulk modulus, $$$K_{\mathrm{f}}$$$, so that $$$\rho = \rho_{0}\exp(P_{\mathrm{f}}/K_{\mathrm{f}})$$$. This assumption may be used in the poromechanics Kernels and Materials, or a more general form may be used.

  • $$$\rho_{\mathrm{mat}}$$$, the mass density of the material, also called the "undrained density" in poromechanics literature. This includes the density of the "dry" material (without the fluid), and the density of the fluid: $$$\rho_{\mathrm{mat}} = (1-\phi)\rho_{\mathrm{dry}} + \phi\rho_{\mathrm{fluid}}$$$. Here $$$(1 - \phi)\rho_{\mathrm{dry}}$$$ is the "drained" density of the media, i.e., the density in the absence of any fluid. The $$$\rho_{\mathrm{fluid}}$$$ density is of the fluid, e.g. $$$\sum_{\mathrm{phases}}S_{\mathrm{phase}}\rho_{\mathrm{phase}}$$$.

  • $$$\sigma^{\mathrm{tot}}$$$, the total stress. An externally applied mechanical force will create a nonzero $$$\sigma^{\mathrm{tot}}$$$, and conversely, resolving $$$\sigma_{ij}^{\mathrm{tot}}$$$ into forces yields the forces on nodes in the finite-element mesh.

  • $$$\alpha$$$, the Biot coefficient. This obeys $$$0 \leq \alpha \leq 1$$$. For a multi-phase, the Biot coefficient is often chosen to be $$$\alpha=1$$$. The Biot coefficient is interpreted physically by the following. If, by pumping fluid into a porous material, the porepressure is increased by $$$\Delta P_{\mathrm{f}}$$$, and at the same time a mechanical external force applies an incremental pressure equaling $$$\alpha \Delta P_{\mathrm{f}}$$$, then the volume of the porous solid remains static. (During this process, the porevolume will change, however, as quantified below.)

  • $$$\sigma^{\mathrm{eff}}$$$, the effective stress. $$$\sigma^{\mathrm{eff}}_{ij} = \sigma_{ij}^{\mathrm{tot}} + \alpha \delta_{ij}P_{\mathrm{f}}$$$. Effective stress, and not total stress, enters into the equations of linear elasticity. Effective stress, and not total stress, enters into the the equations of plasticity. Therefore, MOOSE uses effective stress, and not total stress, internally. If you need to input or output total stress, you must subtract $$$\alpha \delta_{ij}P_{\mathrm{f}}$$$ from MOOSE's stress.

  • $$$b_{i}$$$, the components of a force-density, for example, the gravitational acceleration.

  • $$$E_{ijkl}$$$, components of the elasticity tensor. In linear elasticity $$$\sigma_{ij}^{\mathrm{eff}} = E_{ijkl}\epsilon_{kl}$$$. $$$E_{ijkl}$$$ is often called the "drained" elasticity tensor.

  • $$$C_{ijkl}$$$, components of the compliance tensor. In linear elasticity $$$\epsilon_{ij} = C_{ijkl}\sigma_{ij}^{\mathrm{eff}}$$$. $$$C_{ijkl}$$$ is often called the "drained" compliance tensor.

  • $$$K$$$, the solid bulk modulus. $$$K$$$ is often called the "drained" bulk modulus. If $$$\epsilon_{ij} = C_{ijkl}\sigma^{\mathrm{eff}}_{kl}$$$, then $$$1/K = \delta_{ij}\delta_{kl}C_{ijkl}$$$. It is common in poroelasticity to introduce the "grain bulk modulus", which is $$$1/K_{\mathrm{g}} = (1-\alpha)/K$$$ and describes the compressibility of individual grains of the pororous material, and this quantity appears in the Biot Modulus $$$M$$$, and elsewhere in poromechanics formulae. It is quite common to use $$$1/K_{\mathrm{g}}=0$$$.

  • $$$g_{i}$$$, components of gravitational acceleration, pointing downwards. Eg, usually $$$g = (0, 0, -9.8)$$$.

  • $$$k_{ij}$$$, components of the permeability tensor.

  • $$$k_{\mathrm{rel}}$$$, the relative permeability. It is usually a function of fluid saturation, with $$$k_{\mathrm{rel}} = 0$$$ for $$$S=0$$$, and $$$k_{\mathrm{rel}}=1$$$ for $$$S=1$$$.

  • $$$\mu$$$, the fluid dynamic viscosity.

  • $$$M$$$, the Biot modulus. $$$1/M = (1-\alpha)(\alpha - \phi)/K + \phi/K_{\mathrm{f}}$$$. The quantity $$$1/M$$$ is sometimes called the "storativity" in literature.

Equations

In a Section below, some of the following equations are motivated, derived and explained. Here the equations are simply enumerated with little discussion.

Mechanical Equations

Conservation of linear momentum reads \begin{equation} \nabla_{i}\sigma_{ij}^{\mathrm{eff}} - \alpha \nabla_{j}P_{\mathrm{f}} + \rho_{\mathrm{mat}} b_{j} = \rho_{\mathrm{mat}} \frac{\partial v_{j}}{\partial t}\ . \end{equation} Notice that this is exactly the usual TensorMechanics equation, but it has been expressed in terms of effective stress instead of total stress. Conservation of angular momentum reads $$$\sigma_{ij}^{\mathrm{eff}} = \sigma_{ji}^{\mathrm{eff}}$$$, and the conservation of energy equation is unimportant here. The constitutive law is: \begin{equation} \sigma_{ij}^{\mathrm{eff}} = E_{ijkl}\epsilon_{kl} \ . \end{equation}

These equations are very similar to standard solid mechanics, and may be obviously generalised to large strains, and plasticity, with the important observation that only the effective stress, $$$\sigma^{\mathrm{eff}}$$$, enters the constitutive law, plasticity, insitu stresses, etc, and not the total stress. The only exception to this rule is that when applying Neumann BCs for the displacement variables, the total stress is being specified.

Evolution of Porosity

The porosity evolves according to \begin{equation} \frac{\partial \phi}{\partial t} = \frac{(1-\alpha)(\alpha - \phi)}{K}\frac{\partial P_{\mathrm{f}}}{\partial t} + (\alpha - \phi)\frac{\partial \epsilon_{ii}}{\partial t} \ . \end{equation} This equation may be solved in closed form: \begin{equation} \phi = \alpha - C \exp\left[ \frac{\alpha - 1}{K}P_{\mathrm{f}} - \epsilon_{ii} \right] \ . \end{equation}

Fluid Equations

For a single, unsaturated fluid: \begin{equation} \frac{\partial}{\partial t} \phi\rho S + \phi\rho S \frac{\partial\epsilon_{ii}}{\partial t} = \nabla_{i} \left[ \frac{\rho k_{ij} k_{\mathrm{rel}}}{\mu} (\nabla_{j} P_{\mathrm{f}} - \rho g_{j}) \right] + \mbox{SourceDensity} \ . \end{equation} For a multi-phase fluid, each phase, ph, obeys an equation similar to the single-phase case: \begin{equation} \frac{\partial}{\partial t} \phi\rho^{\mathrm{ph}} S^{\mathrm{ph}} + \phi\rho^{\mathrm{ph}} S^{\mathrm{ph}} \frac{\partial\epsilon_{ii}}{\partial t} = \nabla_{i} \left[ \frac{\rho^{\mathrm{ph}} k_{ij} k^{\mathrm{ph,rel}}}{\mu^{\mathrm{ph}}} (\nabla_{j} P^{\mathrm{ph}} - \rho^{\mathrm{ph}} g_{j}) \right] + \mbox{SourceDensity}^{\mathrm{ph}} \ . \end{equation} Poromechanics is commonly used with a fully-saturated ($$$S=1$$$), single-phase fluid that has constant bulk modulus, and in this case, the fluid equation reads: \begin{equation} \frac{1}{M}\dot{P_{\mathrm{f}}} + \alpha\dot{\epsilon_{ii}} = \nabla_{i}\left[ \frac{k_{ij}}{\mu} (\nabla_{j}P_{\mathrm{f}} - \rho_{0}g_{j}) \right] \end{equation} This has been derived from the single-phase equation, above, by using the evolution equation for porosity. The whole equation has been divided by $$$\rho$$$, as is customary, and only $$$\rho_{0}$$$ appears in the gravitational term, as is also customary. Note that the residual for this equation will often be much smaller than the ones for displacement, since $$$M$$$ is often much larger than stress. Therefore, judicious use of MOOSE's "scaling" parameter is necessary to obtain accurate results.

Motivations, Explanations and Derivations

Evolution of porosity

The evolution of porosity is fundamental to the coupling between flow and solid mechanics. Denote the change of a quantity, q, by $$$\Delta q$$$. Recall that the porosity is defined by $$$\phi = V_{\mathrm{f}}/V$$$, and that, by definition of the effective stress, \begin{equation} \Delta \epsilon_{ij} = C_{ijkl}(\sigma_{ij}^{\mathrm{tot}} + \alpha \delta_{ij}\Delta P_{\mathrm{f}}) \ . \end{equation}

Taking the trace of this equation, and using $$$V^{-1}\Delta V = \Delta \epsilon_{ii}$$$ yields \begin{equation} \frac{\Delta V}{V} = \delta_{ij}C_{ijkl}(\sigma_{ij}^{\mathrm{tot}}+ \alpha \delta_{ij}\Delta P_{\mathrm{f}}) \ . \end{equation} A similar equation can be formed for the change in porevolume: \begin{equation} \frac{\Delta V_{\mathrm{f}}}{V_{\mathrm{f}}} = A_{ij} (\sigma_{ij}^{\mathrm{tot}} + B \delta_{ij}\Delta P_{\mathrm{f}}) \ . \end{equation} It is not too difficult, using the Betti-Maxwell reciprocal theorem, to establish that $$$A_{ij} = \alpha C_{ijkl}\delta_{kl}/\phi$$$, and that $$$B = 1 + \phi - \phi/\alpha$$$. Substituting these into the RHS of the above equation, and rearranging yields, \begin{equation} \frac{\Delta V_{\mathrm{f}}}{V_{\mathrm{f}}} = \frac{\alpha}{\phi}\delta_{ij}C_{ijkl} \left[ \Delta \sigma_{kl}^{\mathrm{tot}} + \alpha \delta_{kl} \Delta P_{\mathrm{f}} \right] + \frac{\delta_{ij}\delta_{kl}C_{ijkl}}{\phi}(1-\alpha)(\alpha-\phi)\Delta P_{f} \end{equation} Using the expression for $$$\Delta V/V$$$ yields \begin{equation} \Delta V_{\mathrm{f}} = V\alpha\Delta\epsilon_{ii} + V\delta_{ij}\delta_{kl}C_{ijkl}(1-\alpha)(\alpha-\phi)\Delta P_{f} \end{equation} Now $$$\Delta\phi = V^{-1}\Delta V_{\mathrm{f}} - V_{\mathrm{f}}V^{-2}\Delta V$$$, so using the definition of $$$K$$$ yields \begin{equation} \frac{\partial \phi}{\partial t} = (\alpha - \phi)\frac{\partial \epsilon_{ii}}{\partial t} + \frac{(1-\alpha)(\alpha - \phi)}{K}\frac{\partial P_{\mathrm{f}}}{\partial t} \ , \end{equation} as required.

A note on the fluid equations

The time-derivative part of the fluid equations contains a term $$$\phi\rho S \frac{\partial\epsilon_{ii}}{\partial t}$$$. This terms is present because the time-derivative is expressing the change in fluid mass in a volume $$$V$$$: \begin{equation} \frac{\partial}{\partial t} \int_{V} \phi\rho S \ . \end{equation} When taking the derivative through the integral, it must be remembered that $$$V$$$ is changing. This change results in the aforementioned term.

Miscellaneous formulae

Undrained bulk modulus: \begin{equation} K_{\mathrm{u}} = K + \alpha^{2}M \ . \end{equation} Skempton coefficient: \begin{equation} B = \frac{K_{\mathrm{u}} - K}{\alpha K_{\mathrm{u}}} \ . \end{equation} \begin{equation} B = \frac{\alpha M}{K_{\mathrm{u}}}\ . \end{equation} Drained Poisson's ratio: \begin{equation} \nu = \frac{3K-2G}{6K+2G} \ . \end{equation} Undrained Poisson's ratio: \begin{equation} \nu_{\mathrm{u}} = \frac{3K_{\mathrm{u}}-2G}{6K_{\mathrm{u}}+2G} \ . \end{equation} Drained bulk modulus and shear modulus in terms of undrained Young and Poisson, and vice versa: \begin{equation} K = \frac{E}{3(1-2\nu)} \ \ \mbox{ and } \ \ G = \frac{E}{2(1+\nu)} \ \ \mbox{ vice versa }\ \ E = \frac{9KG}{3K+G} \ \ \mbox{ and } \ \ \nu = \frac{3K-2G}{6K+2G}\ . \end{equation} Drained bulk modulus and shear modulus in terms of Lame parameters ($$$\lambda$$$ and $$$G$$$), and vice versa: \begin{equation} K = \lambda + \frac{2G}{3} \ \ \mbox{ and } \ \ G = G \ \ \mbox{ vice versa }\ \ \lambda = K - \frac{2G}{3} \ \ \mbox{ and } \ \ G = G \ . \end{equation} Young and Poisson in terms of Lame parameters ($$$\lambda$$$ and $$$G$$$), and vice versa: \begin{equation} E = \frac{G(3\lambda+2G)}{\lambda + G} \ \ \mbox{ and } \ \ \nu = \frac{\lambda}{2(\lambda + G)} \ \ \mbox{ vice versa }\ \ \lambda = \frac{E\nu}{(1+\nu)(1-2\nu)} \ \ \mbox{ and } \ \ G = \frac{E}{2(1+\nu)} \ . \end{equation}

Kernels, Materials and Actions

PoroMechanicsCoupling Kernel

This Kernel implements the term \begin{equation} -\alpha P_{\mathrm{f}} \nabla_{i} \psi \ . \end{equation} Here $$$\alpha$$$ is the Biot Coefficient, which must be provided by a Material (for instance a GenericConstantMaterial with biot_coefficient being a prop_name, or PoroFullSatMaterial), $$$P_{\mathrm{f}}$$$ is the pore pressure (the name of which needs to be specified in the input file), and $$$\psi$$$ is a test function. Therefore, this Kernel may be added to the StressDivergenceTensor Kernel to obtain the expression $$$\nabla_{i}\sigma_{ij}^{\mathrm{eff}} - \alpha \nabla_{j}P_{\mathrm{f}}$$$. The Kernels section of a typical input file looks like:

    [Kernels]
      [./grad_stress_x]
        type = StressDivergenceTensors
        disp_x = disp_x
        disp_y = disp_y
        disp_z = disp_z
        variable = disp_x
        component = 0
      [../]
      [./grad_stress_y]
        type = StressDivergenceTensors
        disp_x = disp_x
        disp_y = disp_y
        disp_z = disp_z
        variable = disp_y
        component = 1
      [../]
      [./grad_stress_z]
        type = StressDivergenceTensors
        disp_x = disp_x
        disp_y = disp_y
        disp_z = disp_z
        variable = disp_z
        component = 2
      [../]
      [./poro_x]
        type = PoroMechanicsCoupling
        variable = disp_x
        porepressure = porepressure
        component = 0
      [../]
      [./poro_y]
        type = PoroMechanicsCoupling
        variable = disp_y
        porepressure = porepressure
        component = 1
      [../]
      [./poro_z]
        type = PoroMechanicsCoupling
        variable = disp_z
        porepressure = porepressure
        component = 2
      [../]
      ...
    []

PoroMechanics Action

The PoroMechanics Action is a shorthand for the 6 Kernels written above (the 3 StressDivergenceTensor and 3 PoroMechanicsCoupling Kernels). Hence the above Kernel block may be replaced by

    [Kernels]
      [./PoroMechanics]
        disp_x = disp_x
        disp_y = disp_y
        disp_z = disp_z
        porepressure = porepressure
      [../]
      ...
    []

This is assuming your Variables are named disp_x, disp_y, disp_z and porepressure.

PoroFullSatTimeDerivative

This Kernel implements the term \begin{equation} \frac{1}{M}\frac{\partial P_{\mathrm{f}}}{\partial t} + \alpha \frac{\partial \epsilon_{ii}}{\partial t} \end{equation} This is the LHS of the fully-saturated single-phase equation, above. The Kernel needs the displacement variables to be specified so that it can correctly provide off-diagonal Jacobian entries. It also needs Material variables that are provided by a PoroFullSatMaterial Material (see below). Inclusion of this into an input file, along with a diffusion term (the right-hand side: Darcy flow) looks like:

    [Kernels]
      [./poro_timederiv]
        type = PoroFullSatTimeDerivative
        variable = porepressure
        disp_x = disp_x
        disp_y = disp_y
        disp_z = disp_z
      [../]
      [./darcy_flow]
        type = CoefDiffusion
        variable = porepressure
        coef = 1.5
      [../]
      ...
    []

Note that because $$$M \gg |\sigma|$$$ and $$$M \gg P_{\mathrm{f}}$$$ in many simulations, it is appropriate to set MOOSE's scaling parameter to ensure the residual for $$$P_{f}$$$ is approximately as large as the residual for the displacement, for instance:

[Variables]
   [./porepressure]
      scaling = 1E3
   [../]
   ...
[]

PoroFullSatMaterial Material

This Material calculates porosity, Biot Modulus and volumetric strain rate, and their derivatives for use in Kernels, as well as the Biot coefficient. These are used in simulations involving a fully-saturated single-phase fluid. Many parameters are required: a typical use is:

[Materials]
  [./poro_material]
    type = PoroFullSatMaterial
    disp_x = disp_x
    disp_y = disp_y
    disp_z = disp_z
    porepressure = porepressure
    porosity0 = 0.1
    biot_coefficient = 0.6
    solid_bulk_compliance = 0.25
    fluid_bulk_compliance = 0.125
  [../]
  ...
[]

Here porosity0 is the porosity at zero volumetric strain and zero porepressure.

The test suite

The following tests indicate that the MOOSE implementation of the above physics is correct. They are part of MOOSE's test suite.

Jacobian of PoroMechanicsCoupling Kernel.

This is tested using PETSc's snes type = test.

Volumetric expansion due to increasing porepressure

The porepressure within a fully-saturated sample is increased: \begin{equation} P_{\mathrm{f}} = t \ . \end{equation} Zero mechanical pressure is applied to the sample's exterior, so that no Neumann BCs are needed on the sample. No fluid flow occurs since the porepressure is increased uniformly throughout the sample, so only the only Kernels needed are the StressDivergenceTensors and the PoroMechanicsCoupling.

The effective stresses should then evolve as $$$\sigma_{ij}^{\mathrm{eff}} = \alpha t \delta_{ij}$$$, and the volumetric strain $$$\epsilon_{00}+\epsilon_{11}+\epsilon_{22} = \alpha t/K$$$. MOOSE produces this result correctly.

Jacobian of PoroFullSatTimeDerivative

This is tested using PETSc's snes type = test.

Undrained oedometer test

A cubic single-element fully-saturated sample has roller BCs applied to its sides and bottom. All the sample's boundaries are impermeable. A constant downwards (normal) velocity, $$$v_{z}$$$, is applied to its top, and the rise in porepressure and effective stress is observed. (Here $$$z$$$ denotes the direction normal to the top face.) There is no fluid flow in the single element.

Under these conditions, and denoting the height ($$$z$$$ length) of the sample by $$$L$$$: \begin{equation} P_{\mathrm{f}} = -\alpha M v_{z}t/L \ . \end{equation} \begin{equation} \sigma_{xx}^{\mathrm{eff}} = (K - \mbox{$\frac{2}{3}$}G)v_{z}t/L \ . \end{equation} \begin{equation} \sigma_{zz}^{\mathrm{eff}} = (K + \mbox{$\frac{4}{3}$}G)v_{z}t/L \ . \end{equation}

MOOSE produces these results correctly.

Porepressure generation of a confined sample

A single-element fully-saturated sample is constrained on all sides and its boundaries are impermeable. Fluid is pumped into the sample via a volumetric source (ie m$$$^{3}$$$/s per cubic meter), and the rise in porepressure is observed.

Denoting the strength of the source by $$$s$$$ (units are s$$$^{-1}$$$), the expected result is \begin{equation} P_{\mathrm{f}} = Mst \end{equation} \begin{equation} \sigma_{ij}^{\mathrm{eff}} = 0 \ . \end{equation}

MOOSE produces these results correctly.

Porepressure generation of an unconfined sample

A single-element fully-saturated sample is constrained on all sides, except its top. All its boundaries are impermeable. Fluid is pumped into the sample via a volumetric source (ie m$$$^{3}$$$/s per cubic meter), and the rise in the top surface, the porepressure, and the stress are observed.

Denoting the strength of the source by $$$s$$$ (units are s$$$^{-1}$$$), the expected result is \begin{equation} \epsilon_{zz} = \frac{\alpha M s t}{K + 4G/3 + \alpha^{2}M} \ . \end{equation} \begin{equation} P_{\mathrm{f}} = M(st - \alpha\epsilon_{zz}) \ . \end{equation} \begin{equation} \sigma_{xx}^{\mathrm{eff}} = (K - 2G/3)\epsilon_{zz} \ . \end{equation} \begin{equation} \sigma_{zz}^{\mathrm{eff}} = (K + 4G/3)\epsilon_{zz} \ . \end{equation}

MOOSE produces these results correctly.

Unconsolidated undrained test

A single-element fully-saturated sample has impermeable boundaries. The sample is squeezed by a uniform mechanical pressure, and the rise in porepressure is observed.

Denoting the mechanical pressure by $$$P^{\mathrm{tot}}$$$, the expected results are \begin{equation} \epsilon_{00} + \epsilon_{11} + \epsilon_{22} = -P^{\mathrm{tot}}/K_{\mathrm{u}} \ , \end{equation} \begin{equation} P_{\mathrm{f}} = BP^{\mathrm{tot}} \ , \end{equation} \begin{equation} \sigma_{22} = -P^{\mathrm{tot}} + \alpha P_{\mathrm{f}} \ . \end{equation}

The mechanical pressure is applied using Neumann BCs for the displacement variables, since recall that these set the total stress, not the effective stress. MOOSE produces the expected results correctly.

Terzaghi consolidation of a drained medium

A saturated sample sits in a bath of water. It is constrained on its sides and bottom. Its sides and bottom are also impermeable. Initially it is unstressed ($$$\sigma_{ij} = 0 = P_{\mathrm{f}}$$$, at $$$t=0$$$). A normal stress, $$$q$$$, is applied to the sample's top. The sample compresses instantaneously due to the instantaneous application of $$$q$$$, and then slowly compresses further as water is squeezed out from the sample's top surface.

This is a classic test. See, for example, Section 2.2 of the online manuscript: Arnold Verruijt "Theory and Problems of Poroelasticity" Delft University of Technology 2013. But note that the "sigma" in that paper is the negative of the stress in TensorMechanics. Denote the sample's height ($$$z$$$-length) by $$$h$$$. Define \begin{equation} p_{0} = \frac{\alpha q M}{S(K + 4G/3) + \alpha^{2}M} \ . \end{equation} This is the porepressure that results from the instantaneous application of $$$q$$$: MOOSE calculates this correctly. The solution for porepressure is \begin{equation} P_{\mathrm{f}} = \frac{4p_{0}}{\pi}\sum_{n=1}^{\infty} \frac{(-1)^{n-1}}{2n-1} \cos \left( \frac{(2n-1)\pi z}{2h} \right) \exp \left( -(2n-1)^{2} \pi^{2} \frac{ct}{4h^{2}} \right) \ . \end{equation} In this equation, $$$c$$$ is the "consolidation coefficient": $$$c = k (K + 4G/3) M/(K + 4G/3 + \alpha^{2} M)$$$, where the permeability tensor is $$$k_{ij} = \mbox{diag}(k, k, k)$$$. The so-called degree-of-consolidation is defined by \begin{equation} U = \frac{u_{z} - u_{z}^{0}}{u_{z}^{\infty} - u_{z}^{0}} \ , \end{equation} where $$$u_{z}$$$ is the vertical displacement of the top surface (downwards is positive), and $$$u_{z}^{0}$$$ is the instantaneous displacement due to the instantaneous application of $$$q$$$, and $$$u_{z}^{\infty}$$$ is the final displacement. This has solution \begin{equation} U = 1 - \frac{8}{\pi^{2}}\sum_{n=1}^{\infty} \frac{1}{(2n-1)^{2}} \exp \left(-(2n-1)^{2}\pi^{2}\frac{ct}{4h^{2}} \right) \ . \end{equation}

MOOSE produces the expected results correctly, as may be seen from the Figures.

terzaghi_p.png
Comparison of expected results and MOOSE results for the Terzaghi consolidation problem.

terzaghi_u.png
Comparison of expected results and MOOSE results for the Terzaghi consolidation problem.

Mandel's consolidation of a drained medium

A sample's dimensions are $$$-a \leq x \leq a$$$ and $$$-b \leq y \leq b$$$, and it is in plane strain (no $$$z$$$ displacement). It is squashed with constant normal force by impermeable, frictionless plattens on its top and bottom surfaces (at $$$y = \pm b$$$). Fluid is allowed to leak out from its sides (at $$$x = \pm a$$$), but all other surfaces are impermeable. This is called Mandel's problem. The interesting feature of this problem (apart from that it can be solved analytically) is that the porepressure in the sample's center actually increases for a short time after application of the force. This is because the leakage of the fluid from the sample's sides causes an apparent softening of the material near those sides. This means stress concentrates towards the sample's center which causes an increase in porepressure. Of course, eventually the fluid totally drains from the sample, and the porepressure is zero. As the fluid drains from the sample's sides the plattens move slowly towards each other.

The solution for porepressure and displacements is given in: AHD Cheng and E Detournay "A direct boundary element method for plane strain poroelasticity" International Journal of Numerical and Analytical Methods in Geomechanics 12 (1988) 551-572. The solution involves rather lengthy infinite series, so I will not write it here.

As is common in the literature, this is simulated by considering the quarter-sample, $$$0\leq x \leq a$$$ and $$$0\leq y\leq b$$$, with impermeable, roller BCs at $$$x=0$$$ and $$$y=0$$$ and $$$y=b$$$. Porepressure is fixed at zero on $$$x=a$$$. Porepressure and displacement are initialised to zero. Then the top ($$$y=b$$$) is moved downwards with prescribed velocity, so that the total force that is induced by this downwards velocity is fixed. The velocity is worked out by solving Mandel's problem analytically, and the total force is monitored in the simulation to check that it indeed remains constant.

mandel_setup.png
The setup of Mandel's problem.

mandel_result.png
Porepressure in Mandel's problem. Note the initial rise in porepressure at the centre of the sample. The sample's size was $$$a=1$$$ in this case.

Poroelastic response of a borehole

A fully-saturated medium contains a fluid with a homogeneous porepressure, but an anisitropic insitu stress. A infinitely-long borehole aligned with the $$$z$$$ axis is instanteously excavated. The borehole boundary is stress-free and allowed to freely drain. This problem is analysed using plane-strain conditions (no $$$z$$$ displacement).

The solution in Laplace space is found in E Detournay and AHD Cheng "Poroelastic response of a borehole in a non-hydrostatic stress field". International Journal of Rock Mechanics and Mining Sciences and Geo mechanics Abstracts 25 (1988) 171-182. In the small-time limit, the Laplace transforms may be performed. There is one typo in the paper. Equation (A4)'s final term should be $$$-(a/r)\sqrt(4ct/(a^2\pi))$$$, and not $$$+(a/r)\sqrt(4ct/(a^2\pi))$$$.

Because realistic parameters are chosen (the drained bulk modulus, the shear modulus, and the fluid bulk modulus are of order GPa, while the stresses and porepressures are of order MPa) the residual for porepressure is much smaller than the residuals for the displacements. Therefore to ensure an accurate solution, a large scaling parameter must be chosen for the porepressure, as mentioned above. Also note that MOOSE's insitu stresses are effective stresses, not total stresses, but the solution in the above paper is expressed in terms of total stresses.

MOOSE produces the expected results, as shown in the figures. These figures were generated for borehole radius $$$a=1$$$, initial porepressure 1MPa, insitu $$$\sigma_{xx}=-2$$$MPa, insitu $$$\sigma_{yy}=-4$$$MPa. The initial increase of porepressure at radius greater than around 1.01m is because the fluid drains from the wall causing the region immediately next to the wall to apparently soften, and stress then transfers to the sample away from the wall, in a similar way to Mandel's problem. In real-life scenarios, fracture can therefore initiate inside the material, and propagate towards the borehole wall. In the comparison of the total stress, the effect of the finite-element discretisation can be seen: since stress is a CONSTANT MONOMIAL in this problem, it is constant within each element: nevertheless, good agreement is still obtained.

The test suite contains two versions of this problem. The "lowres" is run as part of the test harness, while the "highres" version is marked as "heavy". The figures show the results from the "highres" version.

borehole_setup.png
Setup of the poroelastic response of a borehole.

borehole_p.png
The porepressure at two times along the line $$$y=0$$$ in the poroelastic borehole problem.

borehole_s.png
The radial displacement at two times along the line $$$y=0$$$ in the poroelastic borehole problem.

borehole_t.png
The tangential stress at two times along the line $$$y=0$$$ in the poroelastic borehole problem