Poromechanics Solver

Introduction

This section describes the use of the poroelasticity models implemented in GEOS.

Theory

Governing Equations

In our model, the geomechanics (elasticity) equation is expressed in terms of the total stress $\mathbf{\sigma}$ :

$\nabla \mathbf{\sigma} + \rho_b \mathbf{g} = 0$

where it relates to effective stress $\mathbf{\sigma\prime}$ and pore pressure $p$ through Biot’s coefficient $b$ :

$\mathbf{\sigma} = \mathbf{\sigma\prime} - b p\mathbf{I}$

The fluid mass conservation equation is expressed in terms of pore pressure and volumetric (mean) total stress:

$\left( \frac{1}{M} + \frac{b^2}{K_{dr}} \right) \frac{\partial p}{\partial t} + \frac{b}{K_{dr}} \frac{\partial \sigma_v}{\partial t} + \nabla \cdot \mathbf{v}_f = f$

where $M$ is the Biot’s modulus and $K_{dr}$ is the drained bulk modulus.

Unlike the conventional reservoir model that uses Lagranges porosity, in the coupled geomechanics and flow model, Eulers porosity $\phi$ is adopted so the porosity variation is derived as:

$\partial \phi = \left( \frac{b-\phi}{K_s}\right) \partial p + \left( b-\phi \right) \partial \epsilon_v$

where $K_{s}$ is the bulk modulus of the solid grain and $\epsilon_v$ is the volumetric strain.

Parameters

The poroelasticity model is implemented as a main solver listed in <Solvers> block of the input XML file that calls both SolidMechanicsLagrangianSSLE and SinglePhaseFlow solvers. In the main solver, it requires the specification of solidSolverName, flowSolverName, and couplingTypeOption.

The following attributes are supported:

XML Element: SinglePhasePoromechanics

Name	Type	Default	Description
cflFactor	real64	0.5	Factor to apply to the CFL condition when calculating the maximum allowable time step. Values should be in the interval (0,1]
flowSolverName	groupNameRef	required	Name of the flow solver used by the coupled solver
initialDt	real64	1e+99	Initial time-step value required by the solver to the event manager.
isThermal	integer	0	Flag indicating whether the problem is thermal or not. Set isThermal=”1” to enable the thermal coupling
logLevel	integer	0	Sets the level of information to write in the standard output (the console typically). Level 0 outputs no specific information for this solver. Higher levels require more outputs. 1 - Line search information - Solution information (scaling, maximum changes, quality check) - Convergence information - Time step information - Linear solver information - Nonlinear solver information - Solver timers information - Coupling information 2 - The summary of declared fields and coupling
name	groupName	required	A name is required for any non-unique nodes
solidSolverName	groupNameRef	required	Name of the solid solver used by the coupled solver
stabilizationMultiplier	real64	1	Constant multiplier of stabilization strength
stabilizationRegionNames	groupNameRef_array	{}	Regions where stabilization is applied.
stabilizationType	geos_stabilization_StabilizationType	None	StabilizationType. Options are: None- Add no stabilization to mass equation Global- Add jump stabilization to all faces Local- Add jump stabilization on interior of macro elements
targetRegions	groupNameRef_array	required	Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager.
writeLinearSystem	integer	0	Write matrix, rhs, solution to screen ( = 1) or file ( = 2).
LinearSolverParameters	node	unique	XML Element: LinearSolverParameters
NonlinearSolverParameters	node	unique	XML Element: NonlinearSolverParameters

couplingTypeOption: defines the coupling scheme.

The solid constitutive model used here is PoroLinearElasticIsotropic, which derives from ElasticIsotropic and includes an additional parameter: Biot’s coefficient. The fluid constitutive model is the same as SinglePhaseFlow solver. For the parameter setup of each individual solver, please refer to the guideline of the specific solver.

An example of a valid XML block for the constitutive model is given here:

  <Constitutive>
    <PorousElasticIsotropic
      name="porousRock"
      solidModelName="skeleton"
      porosityModelName="skeletonPorosity"
      permeabilityModelName="skeletonPerm"/>

    <ElasticIsotropic
      name="skeleton"
      defaultDensity="0"
      defaultYoungModulus="1.0e4"
      defaultPoissonRatio="0.2"/>

    <CompressibleSinglePhaseFluid
      name="fluid"
      defaultDensity="1"
      defaultViscosity="1.0"
      referencePressure="0.0"
      referenceDensity="1"
      compressibility="0.0e0"
      referenceViscosity="1"
      viscosibility="0.0"/>

    <BiotPorosity
      name="skeletonPorosity"
      defaultGrainBulkModulus="1.0e27"
      defaultReferencePorosity="0.3"/>

    <ConstantPermeability
      name="skeletonPerm"
      permeabilityComponents="{ 1.0e-4, 1.0e-4, 1.0e-4 }"/>
  </Constitutive>

Example

    <SinglePhasePoromechanics
      name="PoroelasticitySolver"
      solidSolverName="LinearElasticitySolver"
      flowSolverName="SinglePhaseFlowSolver"
      logLevel="1"
      targetRegions="{ Domain }">
      <LinearSolverParameters
        solverType="gmres"
        preconditionerType="mgr"/>
    </SinglePhasePoromechanics>

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

    <SinglePhaseFVM
      name="SinglePhaseFlowSolver"
      logLevel="1"
      discretization="singlePhaseTPFA"
      targetRegions="{ Domain }"/>
  </Solvers>