Poromechanics¶
Context
In this example, we use a coupled solver to solve a poroelastic Terzaghitype problem, a classic benchmark in poroelasticity. We do so by coupling a single phase flow solver with a smallstrain Lagrangian mechanics solver.
Objectives
At the end of this example you will know:
how to use multiple solvers for poromechanical problems,
how to define finite elements and finite volume numerical methods.
Input file
This example uses no external input files and everything required is contained within two GEOS input files located at:
inputFiles/poromechanics/PoroElastic_Terzaghi_base_direct.xml
inputFiles/poromechanics/PoroElastic_Terzaghi_smoke.xml
Description of the case¶
We simulate the consolidation of a poroelastic fluidsaturated column of height having unit crosssection. The column is instantaneously loaded at time = 0 s with a constant compressive traction 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.
GEOS 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 GEOS, namely
where is the initial pressure, constant throughout the column, and is the consolidation coefficient (or diffusion coefficient), with
Biot’s coefficient
the uniaxial bulk modulus, Young’s modulus, and Poisson’s ratio
the constrained specific storage coefficient, porosity, the bulk modulus, and the fluid compressibility
the isotropic permeability
the fluid viscosity
The characteristic consolidation time of the system is defined as . Knowledge of is useful for choosing appropriately the timestep sizes that are used in the discrete model.
Coupled solvers¶
GEOS is a multiphysics 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.
We define and characterize each singlephysics solver separately. Then, we define a coupling solver between these singlephysics solvers as another, separate, solver. This approach allows for generality and flexibility in our multiphysics resolutions. The order in which these solver specifications is done is not important. It is important, though, to instantiate each singlephysics solver with meaningful names. The names given to these singlephysics solver instances will be used to recognize them and create the coupling.
To define a poromechanical coupling, we will effectively define three solvers:
the singlephysics flow solver, a solver of type
SinglePhaseFVM
called hereSinglePhaseFlowSolver
(more information on these solvers at Singlephase Flow Solver),the smallstress Lagrangian mechanics solver, a solver of type
SolidMechanicsLagrangianSSLE
called hereLinearElasticitySolver
(more information here: Solid Mechanics Solver),the coupling solver that will bind the two singlephysics solvers above, an object of type
SinglePhasePoromechanics
called herePoroelasticitySolver
(more information at Poromechanics Solver).
Note that the name
attribute of these solvers is
chosen by the user and is not imposed by GEOS.
The two singlephysics solvers are parameterized as explained in their respective documentation.
Let us focus our attention on the coupling solver.
This solver (PoroelasticitySolver
) uses a set of attributes that specifically describe the coupling for a poromechanical 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 discretization method (FE1
, defined further in the input file),
and the target regions (here, we only have one, Domain
).
<SinglePhasePoromechanics
name="PoroelasticitySolver"
solidSolverName="LinearElasticitySolver"
flowSolverName="SinglePhaseFlowSolver"
logLevel="1"
targetRegions="{ Domain }">
<LinearSolverParameters
directParallel="0"/>
</SinglePhasePoromechanics>
<SolidMechanicsLagrangianSSLE
name="LinearElasticitySolver"
timeIntegrationOption="QuasiStatic"
logLevel="1"
discretization="FE1"
targetRegions="{ Domain }"/>
<SinglePhaseFVM
name="SinglePhaseFlowSolver"
logLevel="1"
discretization="singlePhaseTPFA"
targetRegions="{ Domain }"/>
</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 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 FiniteElement section.
The finite volume method requires the specification of a discretization scheme. Here, we use a twopoint flux approximation as described in the dedicated documentation (found here: FiniteVolume).
<NumericalMethods>
<FiniteElements>
<FiniteElementSpace
name="FE1"
order="1"/>
</FiniteElements>
<FiniteVolume>
<TwoPointFluxApproximation
name="singlePhaseTPFA"/>
</FiniteVolume>
</NumericalMethods>
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 beamlike 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="{ 400 }"
ny="{ 16 }"
nz="{ 16 }"
cellBlockNames="{ cb1 }"/>
</Mesh>
The parameters used in the simulation are summarized in the following table.
Symbol 
Parameter 
Units 
Value 

Young’s modulus 
[Pa] 
10^{4} 

Poisson’s ratio 
[] 
0.2 

Biot’s coefficient 
[] 
1.0 

Porosity 
[] 
0.3 

Fluid density 
[kg/m^{3}] 
1.0 

Fluid compressibility 
[Pa^{1}] 
0.0 

Permeability 
[m^{2}] 
10^{4} 

Fluid viscosity 
[Pa s] 
1.0 

Applied compression 
[Pa] 
1.0 

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 = 1.0 Pa, = 1.111 m^{2} s^{1}, and = 90 s.
Therefore, as shown in the Events
section, we run this simulation for 90 seconds.
Running GEOS¶
To run the case, use the following command:
path/to/geosx i inputFiles/poromechanics/PoroElastic_Terzaghi_smoke.xml
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 ) = (5.00e01) ; ( Rsolid, Rfluid ) = ( 5.00e01, 0.00e+00 )
( R ) = ( 5.00e01 ) ;
Attempt: 0, NewtonIter: 1
( RSolid ) = (4.26e16) ; ( Rsolid, Rfluid ) = ( 4.26e16, 4.22e17 )
( R ) = ( 4.28e16 ) ;
As expected, since we are dealing with a linear problem, the fully implicit solver converges in a single iteration.
Inspecting results¶
This plot compares the analytical pressure solution (continuous lines) at selected times with the numerical solution (markers).
To go further¶
Feedback on this example
This concludes the poroelastic example. For any feedback on this example, 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.