Porous Flow Tutorial Page 04. Adding solid mechanics

In this Page, solid mechanics is added to the thermo-hydro simulation of previous Pages. The equations are discussed in governing equations. Only quasi-static solid mechanics is considered here, without gravity, so the equations read (1) As described previously, is the porepressure, the temperature and is the Biot coefficient. The additional nomenclature used here is

• is the effective stress tensor

• is the total stress tensor

• is the elasticity tensor of the drained porous skeleton

• is the linear thermal expansion coefficient. Note that this is the linear version, in contrast to the volumetric coefficients introduced in Page 1.

Once again, before attempting to write an input file, a rough estimate of the expected nonlinear residuals must be performed, as discussed in convergence criteria. The residual for the Eq. (1) is approximately Corresponding to the choice Pa.m made in Page 02 the choice Pa.m may be made here. This means which is significantly greater than for the fluid equation. Therefore, the displacement variables are scaled by .

Many mechanically-related MOOSE objects (Kernels, BCs, etc) accept the use_displaced_mesh input parameter. For virtually all PorousFlow simulations, it is appropriate to set this to false: use_displaced_mesh = false. This means that the Kernel's residual (or BC's residual, Postprocessor's value, etc) will be evaluated using the undisplaced mesh. This has the great numerical advantage that the solid-mechanics elasticity equations remain linear.

Also, many mechanically-related MOOSE objects require the displacements input parameter. Therefore, it is convenient to put this parameter into the GlobalParams block:

[GlobalParams<<<{"href": "../../syntax/GlobalParams/index.html"}>>>]
  displacements = 'disp_x disp_y disp_z'
  PorousFlowDictator = dictator
  biot_coefficient = 1.0
[]

To model this thermo-hydro-mechanical system, the PorousFlowBasicTHM action needs to be enhanced to read:

[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
  [porepressure]
  []
  [temperature]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 293
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = 1E-8
  []
  [disp_x]
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = 1E-10
  []
  [disp_y]
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = 1E-10
  []
  [disp_z]
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = 1E-10
  []
[]

[PorousFlowBasicTHM<<<{"href": "../../syntax/PorousFlowBasicTHM/index.html"}>>>]
  porepressure<<<{"description": "The name of the porepressure variable"}>>> = porepressure
  temperature<<<{"description": "For isothermal simulations, this is the temperature at which fluid properties (and stress-free strains) are evaluated at.  Otherwise, this is the name of the temperature variable.  Units = Kelvin"}>>> = temperature
  coupling_type<<<{"description": "The type of simulation.  For simulations involving Mechanical deformations, you will need to supply the correct Biot coefficient.  For simulations involving Thermal flows, you will need an associated ConstantThermalExpansionCoefficient Material"}>>> = ThermoHydroMechanical
  gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
  fp<<<{"description": "The name of the user object for fluid properties. Only needed if fluid_properties_type = PorousFlowSingleComponentFluid"}>>> = the_simple_fluid
  eigenstrain_names<<<{"description": "List of all eigenstrain models used in mechanics calculations. Typically the eigenstrain_name used in ComputeThermalExpansionEigenstrain.  Only needed for thermally-coupled simulations with thermal expansion."}>>> = thermal_contribution
  use_displaced_mesh<<<{"description": "Use displaced mesh computations in mechanical kernels"}>>> = false
[]

The boundary conditions used here are roller boundary conditions, as well as boundary conditions that model the effect of the fluid porepressure on the injection area:

[BCs<<<{"href": "../../syntax/BCs/index.html"}>>>]
  [constant_injection_porepressure]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = porepressure
    value<<<{"description": "Value of the BC"}>>> = 1E6
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = injection_area
  []
  [constant_injection_temperature]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
    value<<<{"description": "Value of the BC"}>>> = 313
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = injection_area
  []

  [roller_tmax]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    value<<<{"description": "Value of the BC"}>>> = 0
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = dmax
  []
  [roller_tmin]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    value<<<{"description": "Value of the BC"}>>> = 0
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = dmin
  []
  [roller_top_bottom]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_z
    value<<<{"description": "Value of the BC"}>>> = 0
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'top bottom'
  []
  [cavity_pressure_x]
    type = Pressure<<<{"description": "Applies a pressure on a given boundary in a given direction", "href": "../../source/bcs/Pressure.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = injection_area
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    component<<<{"description": "The component for the pressure"}>>> = 0
    factor<<<{"description": "The magnitude to use in computing the pressure"}>>> = 1E6
    use_displaced_mesh<<<{"description": "Whether to use the displaced mesh."}>>> = false
  []
  [cavity_pressure_y]
    type = Pressure<<<{"description": "Applies a pressure on a given boundary in a given direction", "href": "../../source/bcs/Pressure.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = injection_area
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    component<<<{"description": "The component for the pressure"}>>> = 1
    factor<<<{"description": "The magnitude to use in computing the pressure"}>>> = 1E6
    use_displaced_mesh<<<{"description": "Whether to use the displaced mesh."}>>> = false
  []
[]

The SolidMechanics module of MOOSE provides some useful AuxKernels for extracting effective stresses of interest to this problem (the effective radial stress and the effective hoop stress)

[AuxVariables<<<{"href": "../../syntax/AuxVariables/index.html"}>>>]
  [stress_rr]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [stress_pp]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
[]

[AuxKernels<<<{"href": "../../syntax/AuxKernels/index.html"}>>>]
  [stress_rr]
    type = RankTwoScalarAux<<<{"description": "Compute a scalar property of a RankTwoTensor", "href": "../../source/auxscalarkernels/RankTwoScalarAux.html"}>>>
    rank_two_tensor<<<{"description": "The rank two material tensor name"}>>> = stress
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = stress_rr
    scalar_type<<<{"description": "Type of scalar output"}>>> = RadialStress
    point1<<<{"description": "Start point for axis used to calculate some cylindrical material tensor quantities"}>>> = '0 0 0'
    point2<<<{"description": "End point for axis used to calculate some material tensor quantities"}>>> = '0 0 1'
  []
  [stress_pp]
    type = RankTwoScalarAux<<<{"description": "Compute a scalar property of a RankTwoTensor", "href": "../../source/auxscalarkernels/RankTwoScalarAux.html"}>>>
    rank_two_tensor<<<{"description": "The rank two material tensor name"}>>> = stress
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = stress_pp
    scalar_type<<<{"description": "Type of scalar output"}>>> = HoopStress
    point1<<<{"description": "Start point for axis used to calculate some cylindrical material tensor quantities"}>>> = '0 0 0'
    point2<<<{"description": "End point for axis used to calculate some material tensor quantities"}>>> = '0 0 1'
  []
[]

[FluidProperties<<<{"href": "../../syntax/FluidProperties/index.html"}>>>]
  [the_simple_fluid]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 2E9
    viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 1.0E-3
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 1000.0
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0.0002
    cp<<<{"description": "Constant specific heat capacity at constant pressure (J/kg/K)"}>>> = 4194
    cv<<<{"description": "Constant specific heat capacity at constant volume (J/kg/K)"}>>> = 4186
    porepressure_coefficient<<<{"description": "The enthalpy is internal_energy + P / density * porepressure_coefficient.  Physically this should be 1.0, but analytic solutions are simplified when it is zero"}>>> = 0
  []
[]

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [porosity]
    type = PorousFlowPorosity<<<{"description": "This Material calculates the porosity PorousFlow simulations", "href": "../../source/materials/PorousFlowPorosity.html"}>>>
    porosity_zero<<<{"description": "The porosity at zero volumetric strain and reference temperature and reference effective porepressure and reference chemistry.  This must be a real number or a constant monomial variable (not a linear lagrange or other type of variable)"}>>> = 0.1
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus<<<{"description": "Computes the Biot Modulus, which is assumed to be constant for all time.  Sometimes 1 / BiotModulus is called storativity", "href": "../../source/materials/PorousFlowConstantBiotModulus.html"}>>>
    solid_bulk_compliance<<<{"description": "Reciprocal of the drained bulk modulus of the porous skeleton.  If strain = C * stress, then solid_bulk_compliance = de_ij de_kl C_ijkl.  If the grain bulk modulus is Kg then 1/Kg = (1 - biot_coefficient) * solid_bulk_compliance."}>>> = 2E-7
    fluid_bulk_modulus<<<{"description": "Fluid bulk modulus"}>>> = 1E7
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = aquifer
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  []
  [permeability_caps]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = caps
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
  []

  [thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient<<<{"description": "Computes the effective thermal expansion coefficient, (biot_coeff - porosity) * drained_coefficient + porosity * fluid_coefficient.", "href": "../../source/materials/PorousFlowConstantThermalExpansionCoefficient.html"}>>>
    drained_coefficient<<<{"description": "Volumetric coefficient of thermal expansion of the drained porous skeleton (ie the porous rock without fluid, or with a fluid that is free to move in and out of the rock)"}>>> = 0.003
    fluid_coefficient<<<{"description": "Volumetric coefficient of thermal expansion for the fluid"}>>> = 0.0002
  []
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy<<<{"description": "This Material calculates the internal energy of solid rock grains, which is specific_heat_capacity * density * temperature.  Kernels multiply this by (1 - porosity) to find the energy density of the porous rock in a rock-fluid system", "href": "../../source/materials/PorousFlowMatrixInternalEnergy.html"}>>>
    density<<<{"description": "Density of the rock grains"}>>> = 2500.0
    specific_heat_capacity<<<{"description": "Specific heat capacity of the rock grains (J/kg/K)."}>>> = 1200.0
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal<<<{"description": "This Material calculates rock-fluid combined thermal conductivity by using a weighted sum.  Thermal conductivity = dry_thermal_conductivity + S^exponent * (wet_thermal_conductivity - dry_thermal_conductivity), where S is the aqueous saturation", "href": "../../source/materials/PorousFlowThermalConductivityIdeal.html"}>>>
    dry_thermal_conductivity<<<{"description": "The thermal conductivity of the rock matrix when the aqueous saturation is zero"}>>> = '10 0 0  0 10 0  0 0 10'
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 'caps aquifer'
  []

  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor<<<{"description": "Compute a constant isotropic elasticity tensor.", "href": "../../source/materials/ComputeIsotropicElasticityTensor.html"}>>>
    youngs_modulus<<<{"description": "Young's modulus of the material."}>>> = 5E9
    poissons_ratio<<<{"description": "Poisson's ratio for the material."}>>> = 0.0
  []
  [strain]
    type = ComputeSmallStrain<<<{"description": "Compute a small strain.", "href": "../../source/materials/ComputeSmallStrain.html"}>>>
    eigenstrain_names<<<{"description": "List of eigenstrains to be applied in this strain calculation"}>>> = thermal_contribution
  []
  [thermal_contribution]
    type = ComputeThermalExpansionEigenstrain<<<{"description": "Computes eigenstrain due to thermal expansion with a constant coefficient", "href": "../../source/materials/ComputeThermalExpansionEigenstrain.html"}>>>
    temperature<<<{"description": "Coupled temperature"}>>> = temperature
    thermal_expansion_coeff<<<{"description": "Thermal expansion coefficient"}>>> = 0.001 # this is the linear thermal expansion coefficient
    eigenstrain_name<<<{"description": "Material property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator."}>>> = thermal_contribution
    stress_free_temperature<<<{"description": "Reference temperature at which there is no thermal expansion for thermal eigenstrain calculation"}>>> = 293
  []
  [stress]
    type = ComputeLinearElasticStress<<<{"description": "Compute stress using elasticity for small strains", "href": "../../source/materials/ComputeLinearElasticStress.html"}>>>
  []
[]

[Preconditioning<<<{"href": "../../syntax/Preconditioning/index.html"}>>>]
  active<<<{"description": "If specified only the blocks named will be visited and made active"}>>> = basic
  [basic]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
    petsc_options<<<{"description": "Singleton PETSc options"}>>> = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = ' asm      lu           NONZERO                   2'
  []
  [preferred_but_might_not_be_installed]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = ' lu       mumps'
  []
[]

[Executioner<<<{"href": "../../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = Newton
  end_time = 1E6
  dt = 1E5
  nl_abs_tol = 1E-15
  nl_rel_tol = 1E-14
[]

[Outputs<<<{"href": "../../syntax/Outputs/index.html"}>>>]
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
[]

Finally, some mechanics-related Materials need to be defined

  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 5E9
    poissons_ratio = 0.0
  []
  [strain]
    type = ComputeSmallStrain
    eigenstrain_names = thermal_contribution
  []
  [thermal_contribution]
    type = ComputeThermalExpansionEigenstrain
    temperature = temperature
    thermal_expansion_coeff = 0.001 # this is the linear thermal expansion coefficient
    eigenstrain_name = thermal_contribution
    stress_free_temperature = 293
  []
  [stress]
    type = ComputeLinearElasticStress
  []
[]

An animation of the results is shown in Figure 1.

Figure 1: Displacement (magnified by 100 times) and effective hoop-stress evolution in the borehole-aquifer-caprock system.

The dynamics of this model are fascinating, and readers are encouraged to pause and play with parameters to explore how they effect the final result. In fact, this model is very similar to the "THM Rehbinder" test in PorousFlow's test suite. Rehbinder (Rehbinder, 1995) derived analytical solutions for a similar THM problem, and MOOSE replicates his result exactly:

Figure 2: Comparison between MOOSE and Rehbinder's analytical solution.

Figure 3: Comparison between MOOSE and Rehbinder's analytical solution.

Figure 4: Comparison between MOOSE and Rehbinder's analytical solution.