Tutorial 8: Terzaghi’s poroelastic problem¶

Context

In this tutorial, we use a coupled solver to solve a poroelastic Terzaghi-type problem, a classic benchmark in poroelasticity. We do so by coupling a single phase flow solver with a small-strain Lagrangian mechanics solver.

Objectives

At the end of this tutorial you will know:

how to use multiple solvers for poroelastic problems,

how to define finite elements and finite volume numerical methods.

Input file

This tutorial uses no external input files and everything required is contained within a single GEOSX input file. The xml input file for this test case is located at:

src/coreComponents/physicsSolvers/multiphysics/integratedTests/PoroElastic_Terzaghi_FIM.xml

Description of the case¶

We simulate the consolidation of a poroelastic fluid-saturated column of height $L$ having unit cross-section. The column is instantaneously loaded at time $t$ = 0 s with a constant compressive traction $w$ applied on the face highlighted in red in the figure below. Only the loaded face if permeable to fluid flow, with the remaining parts of the boundary subject to roller constraints and impervious.

Fig. 1 Sketch of the the setup for Terzaghi’s problem.

GEOSX will calculate displacement and pressure fields along the column as a function of time. We will use the analytical solution for pressure to check the accuracy of the solution obtained with GEOSX, namely

$p(x,t) = \frac{4}{\pi} p_0 \sum_{m=0}^{\infty} \frac{1}{2m + 1} \text{exp} \left[ -\frac{(2m + 1)^2 \pi^2 c_c t}{4 L^2} \right] \text{sin} \left[ \frac{(2m+1)\pi x}{2L} \right],$

where $p_0 = \frac{b}{K_vS_{\epsilon} + b^2} |w|$ is the initial pressure, constant throughout the column, and $c_c = \frac{\kappa}{\mu} \frac{K_v}{K_v S_{\epsilon} + b^2}$ is the consolidation coefficient (or diffusion coefficient), with

$b$ Biot’s coefficient
$K_v = \frac{E(1-\nu)}{(1+\nu)(1-2\nu)}$ the uniaxial bulk modulus, $E$ Young’s modulus, and $\nu$ Poisson’r ratio
$S_{\epsilon}=\frac{(b - \phi)(1 - b)}{K} + \phi c_f$ the constrained specific storage coefficient, $\phi$ porosity, $K = \frac{E}{3(1-2\nu)}$ the bulk modulus, and $c_f$ the fluid compressibility
$\kappa$ the isotropic permeability
$\mu$ the fluid viscosity

The characteristic consolidation time of the system is defined as $t_c = \frac{L^2}{c_c}$ . Knowledge of $t_c$ is useful for choosing appropriately the timestep sizes that are used in the discrete model.

Preparing the input files¶

All inputs for this case are contained inside a single XML file. In this tutorial, we focus our attention on the Solvers tags, the NumericalMethods tags, and we will briefly inspect the mesh and field specification tags.

Solvers: setting up a multiphysics coupling¶

GEOSX is a multi-physics tool. Different combinations of physics solvers available in the code can be applied in different regions of the mesh at different moments of the simulation. The XML Solvers tag is used to list and parameterize these solvers.

To specify a coupling between two solvers, as done here, we define and characterize each single-physics solver separately. Then, we define a coupling solver between these single-physics solvers as another, separate, solver. This approach allows for generality and flexibility in our multi-physics resolutions. The order in which these solver specifications is done is not important. It is important, though, to instantiate each single-physic solvers with meaningful names. The names given to these single-physics solver instances will be used to recognize them and create the coupling.

To define a poroelastic coupling, we will effectively define three solvers:

the single-physics flow solver, a solver of type SinglePhaseFVM called here SinglePhaseFlowSolver (more information on these solvers at Singlephase Flow Solver),

the small-stress Lagrangian mechanics solver, a solver of type SolidMechanicsLagrangianSSLE called here LinearElasticitySolver (more information here: Solid Mechanics Solver),

the coupling solver that will bind the two single-physics solvers above, an object of type Poroelastic called here PoroelasticitySolver (more information at Poromechanics Solver).

Note that the name attribute of these solvers is chosen by the user and is not imposed by GEOSX.

The two single-physics solvers are parameterized as explained in their respective documentation, each with their own tolerances, verbosity levels, target regions, and other solver-specific attributes.

Let us focus our attention on the coupling solver. This solver (PoroelasticitySolver) uses a set of attributes that specifically describe the coupling for a poroelastic framework. For instance, we must point this solver to the correct fluid solver (here: SinglePhaseFlowSolver), the correct solid solver (here: LinearElasticitySolver). Now that these two solvers are tied together inside the coupling solver, we have a coupled multiphysics problem defined. More parameters are required to characterize a coupling. Here, we specify the coupling type (FIM, fully implicit method; a choice among several possible options), the discretization method (FE1, defined further in the input file), and the target regions (here, we only have one, Region1).

    <Poroelastic
      name="PoroelasticitySolver"
      solidSolverName="LinearElasticitySolver"
      fluidSolverName="SinglePhaseFlowSolver"
      couplingTypeOption="FIM"
      logLevel="1"
      discretization="FE1"
      targetRegions="{ Region1 }">
      <LinearSolverParameters
	solverType="direct"
        directParallel="0"
        logLevel="0"/>
    </Poroelastic>

    <SolidMechanicsLagrangianSSLE
      name="LinearElasticitySolver"
      timeIntegrationOption="QuasiStatic"
      logLevel="1"
      discretization="FE1"
      targetRegions="{ Region1 }"
      solidMaterialNames="{ skeleton }">
    </SolidMechanicsLagrangianSSLE>

    <SinglePhaseFVM
      name="SinglePhaseFlowSolver"
      logLevel="1"
      discretization="singlePhaseTPFA"
      targetRegions="{ Region1 }"
      fluidNames="{ fluid }"
      solidNames="{ skeleton }">
    </SinglePhaseFVM>
  </Solvers>

Multiphysics numerical methods¶

Numerical methods in multiphysics settings are similar to single physics numerical methods. All can be defined under the same NumericalMethods XML tag. In this problem, we use finite volume for flow and finite elements for solid mechanics. Both of these methods require additional parameterization attributes to be defined here.

As we have seen before, the coupling solver and the solid mechanics solver require the specification of a discretization method called FE1. This discretization method is defined here as a finite element method using linear basis functions and Gaussian quadrature rules. For more information on defining finite elements numerical schemes, please see the dedicated FiniteElementDiscretization section.

The finite volume method requires the specification of a discretization scheme. Here, we use a two-point flux approximation as described in the dedicated documentation (found here: Finite Volume Discretization).

  <NumericalMethods>
    <FiniteElements>
      <FiniteElementSpace
        name="FE1"
        order="1"/>
    </FiniteElements>

    <FiniteVolume>
      <TwoPointFluxApproximation
        name="singlePhaseTPFA"
        fieldName="pressure"
        coefficientName="permeability"/>
    </FiniteVolume>
  </NumericalMethods>

Setting up mesh, material properties and boundary conditions¶

Last, let us take a closer look at the geometry of this simple problem. We use the internal mesh generator to create a beam-like mesh, with one single element along the Y and Z axes, and 21 elements along the X axis. All the elements are hexahedral elements (C3D8) of the same dimension (1x1x1 meters).

  <Mesh>
    <InternalMesh
      name="mesh1"
      elementTypes="{ C3D8 }"
      xCoords="{ 0, 10 }"
      yCoords="{ 0, 1 }"
      zCoords="{ 0, 1 }"
      nx="{ 10 }"
      ny="{ 1 }"
      nz="{ 1 }"
      cellBlockNames="{ cb1 }"/>
  </Mesh>

The parameters used in the simulation are summarized in the following table.

Symbol	Parameter	Units	Value
$E$	Young’s modulus	[Pa]	1.0*10⁴
$\nu$	Poisson’s ration	[-]	0.2
$b$	Biot’s coefficient	[-]	1.0
$\phi$	Porosity	[-]	0.3
$\rho_f$	Fluid density	[kg/m³]	1.0
$c_f$	Fluid compressibility	[Pa^-1]	0.0
$\kappa$	Permeability	[m²]	1.0*10^-4
$\mu$	Fluid viscosity	[Pa s]	1.0
$\|w\|$	Applied compression	[Pa]	1.0
$L$	Column length	[m]	10.0

Material properties and boundary conditions are specified in the Constitutive and FieldSpecifications sections. For such set of parameters we have $p_0$ = 1.0 Pa, $c_c$ = 1.111 m² s^-1, and $t_c$ = 90 s. Therefore, as shown in the Events section, we run this simulation for 90 seconds.

Running the case and inspecting the results¶

Running the case¶

To run the case, use the following command:

path/to/geosx -i src/coreComponents/physicsSolvers/multiphysics/integratedTests/PoroElastic_Terzaghi_FIM.xml

Inspecting the console output¶

Here, we see for instance the RSolid and RFluid at a representative timestep (residual values for solid and fluid mechanics solvers, respectively)

Attempt:  0, NewtonIter:  0
( RSolid ) = (2.54e-16) ;     ( Rsolid, Rfluid ) = ( 2.54e-16, 8.44e-06 )
( R ) = ( 8.44e-06 ) ;
Attempt:  0, NewtonIter:  1
( RSolid ) = (2.22e-16) ;     ( Rsolid, Rfluid ) = ( 2.22e-16, 2.76e-20 )
( R ) = ( 2.22e-16 ) ;
Last LinSolve(iter,res) = (   1, 1.15e-11 ) ;

As expected, since we are dealing with a linear problem, the fully implicit solver coverges in a single iteration.

Inspecting results¶

This plot compares the analytical pressure solution (continuous lines) at selected times with the numerical solution (markers).

(Source code)

../../../../_images/Tutorial-1_01_004.png

To go further¶

Feedback on this tutorial

This concludes the poroelastic tutorial. For any feedback on this tutorial, please submit a GitHub issue on the project’s GitHub page.

For more details

More on poroelastic multiphysics solvers, please see Poromechanics Solver.

More on numerical methods, please see Numerical Methods.

More on functions, please see Functions.