Thermoporoelastic Consolidation

Context

Thermoporoelastic consolidation is a typical fully coupled problem which involves solid deformation, fluid flow and heat transfer in saturated porous media. In this example, we use the GEOS coupled solvers to solve a one-dimensional thermoporoelastic consolidation problem with a non-isothermal boundary condition, and we verify the accuracy of the results using the analytical solution provided in (Bai, 2005)

InputFile

This example uses no external input files and everything required is contained within two GEOS input files located at:

inputFiles/thermoPoromechanics/ThermoPoroElastic_consolidation_base.xml

inputFiles/thermoPoromechanics/ThermoPoroElastic_consolidation_benchmark_fim.xml

Description of the case

We simulate the consolidation of 1D thermoporoelastic column subjected to a surface traction stress of 1 Pa applied on the top surface, with a surface temperature of 50 degrees Celsius and a pore pressure of 0 Pa. The initial temperature of the saturated soil is 0 degrees Celsius. The soil column is insulated and sealed everywhere, except at the top surface. The problem setup is illustrated below.

../../../../../../_images/sketch5.png — Fig. 79 Sketch of the problem (taken from (Gao and Ghassemi, 2019)).

The coupled dynamics experienced by the system are described in (Gao and Ghassemi, 2019) and summarized below. The model first experiences continuous settlement (contraction). Initially, the settlement caused by the drainage of the fluid (effective stress increase) and the compression of the solid matrix is larger than the expansion due to the increase of temperature in the region close to the surface on which a higher temperature is applied. As the temperature diffuses further into the domain, it gradually rebounds (expansion) and reaches a final status.

For this example, we focus on the Solvers, the Constitutive, and the FieldSpecifications tags of the GEOS input file.

Solvers

As demonstrated in this example, to setup a thermoporomechanical coupling, we need to define three different solvers in the Solvers part of the XML file:

the mechanics solver, a solver of type SolidMechanicsLagrangianSSLE called here solidMechSolver (more information here: Solid Mechanics Solver),

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

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

  <SinglePhaseFVM
    name="flowSolver"
    logLevel="1"
    discretization="tpfaFlow"
    temperature="273.0"
    isThermal="1"
    targetRegions="{ Domain }">
    <NonlinearSolverParameters
      newtonMaxIter="100"
      newtonMinIter="0"
      newtonTol="1.0e-6"/>
    <LinearSolverParameters
      directParallel="0"/>
  </SinglePhaseFVM>

the coupling solver (SinglePhasePoromechanics) that will bind the two single-physics solvers above, which is named as thermoPoroSolver (more information at Poromechanics Solver).

  <SinglePhasePoromechanics
    name="thermoPoroSolver"
    solidSolverName="solidMechSolver"
    flowSolverName="flowSolver"
    isThermal="1"
    logLevel="1"
    targetRegions="{ Domain }">
    <NonlinearSolverParameters
      couplingType="FullyImplicit"
      newtonMaxIter="200"/>
    <LinearSolverParameters
      directParallel="0"/>
  </SinglePhasePoromechanics>

To request the simulation of the temperature evolution, we set the isThermal flag of the coupling solver to 1. With this choice, the degrees of freedom are the cell-centered pressure, the cell-centered temperature, and the mechanical displacements at the mesh nodes. The governing equations consist of a mass conservation equation, an energy balance equation, and a linear momentum balance equation. In the latter, the total stress includes both a pore pressure and a temperature contribution. Note that in the coupling solver, we set the couplingType to FullyImplicit to require a fully coupled, fully implicit solution strategy.

Constitutive laws

A homogeneous and isotropic domain with one solid material is assumed, and its mechanical properties and associated fluid rheology are specified in the Constitutive section. We use the constitutive parameters specified in (Bai, 2005) listed in the following table.

Symbol	Parameter	Unit	Value
$E$	Young’s modulus	[Pa]	6000.0
$\nu$	Poisson’s ratio	[-]	0.4
$\alpha$	Thermal expansion coef.	[T^(-1)]	9.0x10^-7
$\phi$	Porosity	[-]	0.20
$b$	Biot’s coefficient	[-]	1.0
$\rho C$	Heat capacity	[J/(m^3.K)]	167.2x10³
$\mu$	Fluid viscosity	[Pa.s]	1^-3
$k^T$	Thermal conductivity	[J/(m.s.K)]	836
$k$	Permeability	[m^2]	4.0x10^-9

The bulk modulus, the Young’s modulus, and the thermal expansion coefficient are specified in the ElasticIsotropic solid model. Note that for now the solid density is constant and does not depend on temperature. Given that the gravity vector has been set to 0 in the XML file, the value of the solid density is not used in this simulation.

    <ElasticIsotropic
      name="rockSolid"
      defaultDensity="2400"
      defaultBulkModulus="1e4"
      defaultShearModulus="2.143e3"
      defaultDrainedLinearTEC="3e-7"/>

The porosity and Biot’s coefficient (computed from the grainBulkModulus) appear in the BiotPorosity model. In this model, the porosity is updated as a function of the strain increment, the change in pore pressure, and the change in temperature.

    <BiotPorosity
      name="rockPorosity"
      grainBulkModulus="1.0e27"
      defaultReferencePorosity="0.2"
      defaultPorosityTEC="3e-7"/>

The heat capacity is provided in the SolidInternalEnergy model. In the computation of the internal energy, the referenceTemperature is set to the initial temperature.

    <SolidInternalEnergy
      name="rockInternalEnergy"
      volumetricHeatCapacity="1.672e5"
      referenceTemperature="0.0"
      referenceInternalEnergy="0.0"/>

The fluid density and viscosity are given in the ThermalCompressibleSinglePhaseFluid. Here, they are assumed to be constant and do not depend on pressure and temperature.

    <ThermalCompressibleSinglePhaseFluid
      name="water"
      defaultDensity="1000"
      defaultViscosity="1e-3"
      referencePressure="0.0"
      referenceTemperature="0.0"
      compressibility="0.0"
      thermalExpansionCoeff="0.0"
      viscosibility="0.0"
      specificHeatCapacity="1.672e2"
      referenceInternalEnergy="0.001"/>

Finally, the permeability and thermal conductivity are specified in the ConstantPermeability and SinglePhaseConstantThermalConductivity, respectively.

    <ConstantPermeability
      name="rockPerm"
      permeabilityComponents="{ 4.0e-9, 4.0e-9, 4.0e-9 }"/>
    
    <SinglePhaseConstantThermalConductivity
      name="thermalCond"
      thermalConductivityComponents="{ 836, 836, 836 }"/>

Initial and boundary conditions

To complete the specification of the problem, we specify two types of fields:

The initial values (the displacements, effective stress, and pore pressure have to be initialized),

The boundary conditions at the top surface (traction, pressure, and temperature) and at the other boundaries (zero-displacement).

This is done in the FieldSpecifications part of the XML file. The attribute initialCondition is set to 1 for the blocks specifying the initial pressure, temperature, and effective stress.

    <FieldSpecification
      name="initialPressure"
      initialCondition="1"
      setNames="{ all }"
      objectPath="ElementRegions/Domain/cb1"
      fieldName="pressure"
      scale="0.0"/>

    <FieldSpecification
      name="initialTemperature"
      initialCondition="1"
      setNames="{ all }"
      objectPath="ElementRegions/Domain/cb1"
      fieldName="temperature"
      scale="273.0"/>

    <FieldSpecification
      name="initialSigma_x"
      initialCondition="1"
      setNames="{ all }"
      objectPath="ElementRegions/Domain/cb1"
      fieldName="rockSolid_stress"
      component="0"
      scale="2.457"/>
    <FieldSpecification
      name="initialSigma_y"
      initialCondition="1"
      setNames="{ all }"
      objectPath="ElementRegions/Domain/cb1"
      fieldName="rockSolid_stress"
      component="1"
      scale="2.457"/>
    <FieldSpecification
      name="initialSigma_z"
      initialCondition="1"
      setNames="{ all }"
      objectPath="ElementRegions/Domain/cb1"
      fieldName="rockSolid_stress"
      component="2"
      scale="2.457"/>

For the zero-displacement boundary conditions, we use the pre-defined set names xneg and xpos, yneg, zneg and zpos to select the boundary nodes. Note that here, we have considered a slab in the y-direction, which is why a displacement boundary condition is applied on zpos and not applied on ypos.

    <FieldSpecification
      name="xconstraint"
      fieldName="totalDisplacement"
      component="0"
      objectPath="nodeManager"
      setNames="{ xneg, xpos }"/>

    <FieldSpecification
      name="yconstraint"
      fieldName="totalDisplacement"
      component="1"
      objectPath="nodeManager"
      setNames="{ yneg }"/>

    <FieldSpecification
      name="zconstraint"
      fieldName="totalDisplacement"
      component="2"
      objectPath="nodeManager"
      setNames="{ zneg, zpos }"/>

On the top surface, we impose the traction boundary condition and the non-isothermal boundary condition specified in (Bai, 2005). We also fix the pore pressure to 0 Pa.

    <Traction
      name="traction"
      objectPath="faceManager"
      tractionType="normal"
      scale="-1.0"
      setNames="{ ypos }"
      functionName="timeFunction"/>

    <FieldSpecification
      name="boundaryPressure"
      objectPath="faceManager"
      fieldName="pressure"
      scale="0.0"
      setNames="{ ypos }"/>

    <FieldSpecification
      name="boundaryTemperature"
      objectPath="faceManager"
      fieldName="temperature"
      scale="323.0"
      setNames="{ ypos }"/>

Inspecting results

We request an output of the displacements, pressure, and temperature using the TimeHistory feature of GEOS. The figures below compare the results from GEOS (dashed line) and the corresponding analytical solution (solid line) as a function of time at different locations of the slab. We obtain a very good match, confirming that GEOS can accurately capture the thermo-poromechanical coupling on this example. The first figure illustrates this good agreement for the pressure evolution.

(Source code)

../../../../../../_images/thermalConsolidationPressureFigure.png

The second figure confirms the good match with the analytical solution for the temperature.

(Source code)

../../../../../../_images/thermalConsolidationTemperatureFigure.png

The third figure shows that GEOS is also able to match the vertical displacement (settlement) analytical solution.

(Source code)

../../../../../../_images/thermalConsolidationDisplacementFigure.png

To go further

Feedback on this example

For any feedback on this example, please submit a GitHub issue on the project’s GitHub page.