Poroelasticity Calculations

NairnMPM can do poroelasticity calculations coupled with stresses and strains through tracking of pore pressure in fluid-saturated media.

Introduction

Poroelasticity is modeled by extension of methods first described by Biot. In brief, each particle tracks a pore pressure that models pressure within a liquid for a material saturated with that liquid. As the material volume changes, the pore pressure increases or decreases. The pore pressure can transport between particles based on Darcy tensor for permittivity of the material for the given liquid. Some poroelasticity material properties control coupling between stresses and strains and the pore pressure.

The extensions from Biot are to model poroelasticity in anisotropic materials and to allow the pore pressure to drop below zero.

Activating Poroelasticity

In scripted files, poroelasticity is activated with the command

Poroelasticity (yesorno),<(p0)>,<(eta)>

In XML input files, poroelasticity is activated with the <Poroelasticity> command, which must be within the <MPMHeader> element:

<Poroelasticity reference="(p0)" viscosity="(eta)">

where

• (YesOrNo) must be "Yes" or "No" to activate or not activate poroelasticity calculations. In XML input files, the presence of a <Poroelasticity> command activates poroelasticity methods. The default is "No".
• (p0)is reference pore pressure (in pressure units). The default (p0) is 0. Unless changed, all particles will be initialized to start at the reference pore pressure.
• (eta) is fluid viscosity for fluid in the pores (and same for all materials). Default is 1 cP in Legacy units or 1 viscosity unit in consistent units (note that Legacy units default to water viscosity, but consistent units default to 1 viscosity unit which may not be water viscosity - it is therefore best to always enter this viscosity and not rely on default setting being for water).

By default, poroelasticity uses update methods analogous to FLIP methods used in mechanics. This update, however, sometimes results in pore pressure oscillations on particles within one cell. Poroelasticity simulations with oscillations can be improved by using periodic FMPM(k) for poroelasticity updates. In fact, because poroelasticity is inherently stiffer then conduction or diffusion, it may require FMPM(k) for stable simulations.

Note that poroelasticity models fluid transport through materials by transport methods nearly identical to those used to model diffusion. Because they share same methods, a simulation can activate poroelasticity (with above commands) or diffusion (with comparable Diffusion commands), but cannot activate them both. Any simulation, however, can combine poroelasticity or diffusion with thermal calculations and conduction. Note that when choosing archiving options, the terms "concentration" and "porepressure" are synonyms or either can be used and the archiving will store concentration terms for diffusion calculations or pore pressure terms for poroelasticity calculations.

Note that poroelasticity calculations are often improved by using FMPM(k) methods. They may even require FMPM(k) for stability.

Poroelasticity Material Properties

The constitutive law for an isotropic, poroelastic material is

      [math]\displaystyle{ \mathbf{\sigma} = \mathbf{C} \mathbf{\varepsilon} - \alpha dp }[/math]

where [math]\displaystyle{ \mathbf{C} }[/math] is the compliance tensor, [math]\displaystyle{ \alpha }[/math] is the poroelasticity Biot coefficient (dimensionless from 0 to 1) and [math]\displaystyle{ dp=p-p_0 }[/math] is difference in pore pressure relative to the reference pore pressure. When the bulk material volume is at some reference volume and pore pressure is [math]\displaystyle{ p_0 }[/math], total stress on the bulk is zero. In most simulations, the reference pore pressure is zero and reference volume is volume in the absence of applied stress.

The pore pressure transport between particles is controlled by poroelasticity implementation of Darcy's law:

      [math]\displaystyle{ C_T{dp\over dt} =-{1\over \eta}\nabla \cdot k\nabla p - \alpha {d\varepsilon\over dt} }[/math]

where [math]\displaystyle{ \eta }[/math] is fluid viscosity, [math]\displaystyle{ k }[/math] is Darcy's law permittivity, and [math]\displaystyle{ C_T=1/Q }[/math] is a pore pressure "capacity" defined from:

      [math]\displaystyle{ Q = {K_u-K\over \alpha^2} }[/math]

Here [math]\displaystyle{ K }[/math] is the solid material's bulk modulus and [math]\displaystyle{ K_u }[/math] is the material's undrained bulk modulus. This property is the bulk modulus of the saturated material under conditions where there is no change in pore pressure. For all materials [math]\displaystyle{ K_u\gt K }[/math].

The Biot theory is extended to anisotropic materials by replacing [math]\displaystyle{ \alpha }[/math] and [math]\displaystyle{ k }[/math] with tensors having three principle values for three symmetric directions of an orthotropic material. For a transversely isotropic material, the tensors are defined by two principle values in the axial and transverse directions.

Supported Materials

Poroelasticity is relatively easy to add to any material model, but currently is only implemented in the following selected materials

• Isotropic Material
• Isotropic Softening Material
• Transversely Isotropic Material 1 (or 2) - always uses large rotation mode whether set or not
• Orthotropic Material - always uses large rotation mode whether set or not

When doing poroelasticity calculations, each material must specify [math]\displaystyle{ \alpha }[/math], [math]\displaystyle{ K_u }[/math], and [math]\displaystyle{ k }[/math]. Anisotropic materials must specify values of [math]\displaystyle{ \alpha }[/math] and [math]\displaystyle{ k }[/math] along the each symmetry plane direction. The property names for entering these material properties are in the poroelasticity properties tables.

Archived Pore Pressure and Stress

When the feature to archive pore pressure is activated, the particle archive files will include pore pressure and pore pressure gradients. The results will be archived in pressure units.

In addition, you can archive pore pressure in global archive methods or in VTK archiving custom task (in pressure units).

The standard archiving of particle stress is actually an effective stress that combines stress in the liquid (the pore pressure) with stress in the porous material. It is the stress an outside observer would see and the see relevant for tracking forces in MPM calculations. You can get stress in the porous material, by adding pore pressure scaled by Biot coefficient to the particle stress:

      Solid material stress: [math]\displaystyle{ \sigma^{(s)} = \sigma + \alpha p }[/math]

where [math]\displaystyle{ \alpha }[/math] is the Biot coefficient tensor and [math]\displaystyle{ p }[/math] is pore pressure.

Negative Pore Pressure

The Biot methods are limited to media saturated with fluid with non-negative pore pressure. If the material drops below full saturation, the empty space would result on zero pore pressure caused by the fluid. The implementation of poroelasticity here was extended for account for less then full saturation by allowing the pore pressure the negative. But, a negative pore pressure does not imply tension on the fluid. Instead, pore pressure tracking follows these rules:

1. Each time step will start at some pore pressure [math]\displaystyle{ p }[/math] and undergo some change in pore pressure caused by volume change of the material or by Darcy law flow between particles denoted by [math]\displaystyle{ dp }[/math]. If needed the pore pressure increment using in evaluating bulk stress is changed to [math]\displaystyle{ dp^* }[/math]
2. If [math]\displaystyle{ p\ge0 }[/math] and [math]\displaystyle{ p+dp\ge0 }[/math], the poroelasticity proceeds by standard Biot methods and the materials is fully saturated with fluid.
3. If [math]\displaystyle{ p\lt 0 }[/math] and [math]\displaystyle{ p+dp\lt 0 }[/math], the material is partially saturated. The tracked pore pressure is updated to [math]\displaystyle{ p+dp\lt 0 }[/math], but this negative pore pressure plays no rule in bulk stress (in other words, [math]\displaystyle{ dp^*=0 }[/math] for coupling to stresses).
4. If [math]\displaystyle{ p\ge0 }[/math] but [math]\displaystyle{ p+dp\lt 0 }[/math], the material is transitioning from full to partial saturation. The tracked pore pressure is updated to [math]\displaystyle{ p+dp\lt 0 }[/math], but only the positive fraction of pore pressure is included in bulk stress (in other words, [math]\displaystyle{ dp^*=-p }[/math] for coupling to stresses).
5. If [math]\displaystyle{ p\lt 0 }[/math] but [math]\displaystyle{ p+dp\ge0 }[/math], the material is transitioning from partial to full saturation. The tracked pore pressure is updated to [math]\displaystyle{ p+dp\ge0 }[/math], but only the positive fraction of pore pressure is included in bulk stress (in other words, [math]\displaystyle{ dp^*=p+dp }[/math] for coupling to stresses).

Whenever [math]\displaystyle{ p\lt 0 }[/math], the magnitude of [math]\displaystyle{ p }[/math] is a measure the empty pore volume.

Poroelasticity Boundary Conditions

When poroelasticity is activated, you can set pore pressure boundary conditions as follows:

• You can set pore pressure on the grid.
• You can set a pore pressure flux on particle surfaces.
• Apply initial pore pressure field to particles using a PropertyRamp Custom Task.
• Rigid particles can provide moving pore pressure boundary conditions.

Note that setting initial particle pore pressure different than the reference pore pressure will cause strains to immediately evolve toward the changed state. The net effect will be an instantaneous "impact" of pore pressure change that might cause undesirable dynamic effects.

Also note that pore pressure flux is set by specifying a flux for volume fraction of pore fluid in the solid. A change in fluid volume fraction will translate into a change in pore pressure. In the absence of Darcy's flow away form the particle, an estimated change in pore pressure, dp, due to chance in pore fluid volume fraction, dv_F for an isotropic material is:

      [math]\displaystyle{ dp = {K\over K_u} Q dv_F }[/math]

A more involved relation in needed to get dp from dv_F for anisotropic materials. In typical simulations, the volume fraction flux is small (much less than one when multiplied by particle area) and can be adjusted to get desired pore pressure changes.