THM-Rehbinder test description
Introduction
Rehbinder (Rehbinder, 1995) derives analytical solutions of cylindrically-symmetric and spherically-symmetric THM problems in certain limits. The solution is a steady-state THM solution of a pressurised and heated cavity of radius in a finite domain of radius . In the cylindrically-symmetric case the cavity is a cylinder of radius and the domain is a cylinder of radius . In this document only the cylindrically-symmetric situation is explored. The cylindrical coordinates are denoted by .
Initial and boundary conditions
The initial conditions are zero porepressure, ; zero temperature, ; and zero displacement . The boundary conditions are at the cavity wall are
The latter condition implies the total radial stress is at the cavity wall, corresponding to the fluid in the cavity pushing on the cavity wall. The boundary conditions at the outer radius are
and either or . In this document the former is called the free outer'' boundary condition, while the latter is called the fixed outer'' boundary condition
Simplifying assumptions
In order to derive the solution, Rehbinder makes various simplifications. Translated to the language of the other PorousFlow documents these are as follows.
There is no gravity.
All quantities depend only on the radial coordinate, , and not the angular coordinate or the axial coordinate .
Plane-strain is assumed, so .
There is only one fully-saturated fluid phase that contains one component.
The fluid relative permeability is unity.
The fluid density is constant. In the numerical simulation below, the density is assumed to obey , with , and .
The fluid dynamic viscosity is constant. In the numerical simulation below, the dynamic viscosity is .
The Biot coefficient is unity .
The fluid internal energy is assumed to be linear in the temperature. In the numerical simulation below, where .
The fluid enthalpy is assumed to be equal to the fluid internal energy.
The velocity of the solid skeleton is zero: . To model this using PorousFlow, the
{\tt VolumetricExpansion}'' Kernels are not included.
The rock-matrix density is constant. In the numerical simulation below, the value is chosen.
The rock-matrix specific heat capacity is constant. In the numerical simulation below, the value is chosen.
The rock deforms elastically (so there is no plastic heating). In the numerical simulation below, the Young's modulus is and the Poisson's ratio is .
The porosity is constant. In the simulation below is chosen.
The permeability is constant. In the simulation below is chosen.
The thermal conductivity of the rock-fluid system is constant. In the simulation below is chosen.
Steady-state heat flow, fluid-flow and mechanical deformation has been reached (Rehbinder equation (41)). Using the aove assumptions, this is true if (the left-hand-side is measured in seconds, the right-hand-side in metres). In the simulation below, steady-state is achieved by not including any time-derivative Kernels.
The liquid flow is quasi-stationary, in the sense that the diffusion of heat is much slower than the diffusion of pore pressure (Rehbinder equation (32)). Using the above assumptions this may be written as where is the thermal conductivity of the rock-fluid system, and . Using the above values, this condition reads which is clearly satisfied.
The liquid flow is quasi-stationary, in the sense that a pressure change in the cavity is transmitted instantaneously through the pores to the outer boundary (Rehbinder equation (33)). Using the above assumptions this may be written as . In the simulation below, and are chosen.
The heat is mainly conducted through the matrix, and a negligible part is convected with the flux. This is true if the Peclet number is very small (Rehbinder equation (35))
(Rehbinder writes this in terms of a reference temperature, , which is chosen here to be .) Using the above values this is . In the simulation below, and are chosen, yielding .
Rehbinder's derives the solution of this THM problem as an expansion in the Peclet number. I shall not write the solution here as it is fairly lengthy. The paper contains a few minor typos and they are corrected in the accompanying Python file:
#!/usr/bin/env python3
import os
import sys
import numpy as np
import matplotlib.pyplot as plt
def rehbinder(r):
# Results from Rehbinder with parameters used in the MOOSE simulation.
# Rehbinder's manuscript contains a few typos - I've corrected them here.
# G Rehbinder "Analytic solutions of stationary coupled thermo-hydro-mechanical problems" Int J Rock Mech Min Sci & Geomech Abstr 32 (1995) 453-463
poisson = 0.2
thermal_expansion = 1E-6
young = 1E10
fluid_density = 1000
fluid_specific_heat = 1000
permeability = 1E-12
fluid_viscosity = 1E-3
thermal_conductivity = 1E6
P0 = 1E6
T0 = 1E3
Tref = T0
r0 = 0.1
r1 = 1.0
xi = r / r0
xi1 = r1 / r0
Peclet = fluid_density * fluid_specific_heat * thermal_expansion * Tref * young * permeability / fluid_viscosity / thermal_conductivity / (1 - poisson)
That0 = T0 / Tref
sigmahat0 = -P0 * (1 - poisson) / thermal_expansion / Tref / young
Tzeroth = That0 * (1 - np.log(xi) / np.log(xi1))
Tfirst_pr = 2 * sigmahat0 * That0 * xi * (np.log(xi) - np.log(xi1)) / np.log(xi1)**2
Cone = 2 * That0 * sigmahat0 * (2 + np.log(xi1)) / np.log(xi1)**2
Cone = 2 * That0 * sigmahat0 / np.log(xi1) # Corrected Eqn(87)
Done = 2 * That0 * sigmahat0 * (2 * (xi1 - 1) / np.log(xi1) - 1) / np.log(xi1)**2
Done = 2 * That0 * sigmahat0 * (- 1) / np.log(xi1)**2 # Corrected Eqn(87)
Tfirst_hm = Cone + Done * np.log(xi)
Tfirst = Tfirst_pr + Tfirst_hm
That = Tzeroth + Peclet * Tfirst
T = Tref * That
Pzeroth = -sigmahat0 * (1 - np.log(xi) / np.log(xi1))
Pfirst = 0
Phat = Pzeroth + Peclet * Pfirst
P = thermal_expansion * Tref * young * Phat / (1 - poisson)
g0 = Tzeroth + (1 - 2 * poisson) * Pzeroth / (1 - poisson)
uzeroth_pr = (That0 - sigmahat0 * (1 - 2 * poisson) / (1 - poisson)) * (0.5 * (xi**2 - 1) - 0.25 * (1 - xi**2 + 2 * xi**2 * np.log(xi)) / np.log(xi1)) / xi
uzeroth_pr_xi1 = (That0 - sigmahat0 * (1 - 2 * poisson) / (1 - poisson)) * (0.5 * (xi1**2 - 1) - 0.25 * (1 - xi1**2 + 2 * xi1**2 * np.log(xi1)) / np.log(xi1)) / xi1
# fixed outer boundary
Bzeroth = - ((1 - 2 * poisson) * sigmahat0 + uzeroth_pr_xi1 / xi1) / (1 - 2 * poisson + 1.0 / xi1)
Azeroth = - Bzeroth / xi1**2 - uzeroth_pr_xi1 / xi1
fixed_uzeroth_hm = Azeroth * xi + Bzeroth / xi
fixed_uzeroth = uzeroth_pr + fixed_uzeroth_hm
# free outer boundary
Bzeroth = (xi1**2 * sigmahat0 - xi1 * uzeroth_pr_xi1) / (1 - xi1**2)
Azeroth = (1 - 2 * poisson) * (Bzeroth + sigmahat0)
free_uzeroth_hm = Azeroth * xi + Bzeroth / xi
free_uzeroth = uzeroth_pr + free_uzeroth_hm
ufirst_pr = (1.0 / xi) * (0.5 * (xi**2 - 1) * (2 * Cone - Done) + 0.5 * Done * xi**2 * np.log(xi) + 2 * sigmahat0 * That0 / np.log(xi1)**2 * (xi**3 * np.log(xi) / 3 + (1 - xi**3) / 9 + 0.5 * np.log(xi1) * (1 - xi**2)))
ufirst_pr_xi1 = (1.0 / xi1) * (0.5 * (xi1**2 - 1) * (2 * Cone - Done) + 0.5 * Done * xi1**2 * np.log(xi1) + 2 * sigmahat0 * That0 / np.log(xi1)**2 * (xi1**3 * np.log(xi1) / 3 + (1 - xi1**3) / 9 + 0.5 * np.log(xi1) * (1 - xi1**2)))
# fixed outer boundary
Bfirst = - ufirst_pr_xi1 / xi1 / (1 - 2 * poisson + 1.0 / xi1**2)
Afirst = - Bfirst / xi1**2 - ufirst_pr_xi1 / xi1
fixed_ufirst_hm = Afirst * xi + Bfirst / xi
fixed_ufirst = ufirst_pr + fixed_ufirst_hm
# free outer boundary
Bfirst = xi1 * ufirst_pr_xi1 / (1 - xi1**2)
Afirst = (1 - 2 * poisson) * Bfirst
free_ufirst_hm = Afirst * xi + Bfirst / xi
free_ufirst = ufirst_pr + free_ufirst_hm
fixed_uhat = fixed_uzeroth + Peclet * fixed_ufirst
fixed_u = thermal_expansion * Tref * r0 * fixed_uhat * (1 + poisson) / (1 - poisson) # Corrected Eqn(16)
free_uhat = free_uzeroth + Peclet * free_ufirst
free_u = thermal_expansion * Tref * r0 * free_uhat * (1 + poisson) / (1 - poisson) # Corrected Eqn(16)
return (T, P, fixed_u, free_u)
def moose(fn, rcol, datacol):
try:
f = open(fn)
data = f.readlines()[1:-1]
data = [list(map(float, d.strip().split(","))) for d in data]
data = ([d[rcol] for d in data], [d[datacol] for d in data])
f.close()
except:
sys.stderr.write("Cannot read " + fn + ", or it contains erroneous data\n")
sys.exit(1)
return data
mooser = [0.1 * i for i in range(1, 11)]
fixedT = moose("gold/fixed_outer_T_0001.csv", 0, 4)
fixedP = moose("gold/fixed_outer_P_0001.csv", 0, 4)
fixedu = moose("gold/fixed_outer_U_0001.csv", 0, 4)
freeu = moose("gold/free_outer_U_0001.csv", 0, 4)
rpoints = np.arange(0.1, 1.0, 0.01)
expected = list(zip(*[rehbinder(r) for r in rpoints]))
plt.figure()
plt.plot(rpoints, expected[0], 'k-', linewidth = 3.0, label = 'expected')
plt.plot(fixedT[0], fixedT[1], 'rs', markersize = 10.0, label = 'MOOSE')
plt.legend(loc = 'upper right')
plt.xlabel("r (m)")
plt.ylabel("Temperature (K)")
plt.title("Temperature around cavity")
plt.savefig("temperature_fig.png")
plt.figure()
plt.plot(rpoints, [1E-6 * p for p in expected[1]], 'k-', linewidth = 3.0, label = 'expected')
plt.plot(fixedP[0], [1E-6 * p for p in fixedP[1]], 'rs', markersize = 10.0, label = 'MOOSE')
plt.legend(loc = 'upper right')
plt.xlabel("r (m)")
plt.ylabel("Porepressure (MPa)")
plt.title("Porepressure around cavity")
plt.savefig("porepressure_fig.png")
plt.figure()
plt.plot(rpoints, [1000 * u for u in expected[2]], 'k-', linewidth = 3.0, label = 'expected (fixed)')
plt.plot(fixedu[0], [1000 * u for u in fixedu[1]], 'rs', markersize = 10.0, label = 'MOOSE (fixed)')
plt.plot(rpoints, [1000 * u for u in expected[3]], 'b-', linewidth = 2.0, label = 'expected (free)')
plt.plot(freeu[0], [1000 * u for u in freeu[1]], 'g*', markersize = 13.0, label = 'MOOSE (free)')
plt.legend(loc = 'center right')
plt.xlabel("r (m)")
plt.ylabel("displacement (mm)")
plt.title("Radial displacement around cavity")
plt.savefig("displacement_fig.png")
sys.exit(0)
Comparison with PorousFlow
The PorousFlow module is designed to simulate THM problems. The output compares favourably with Rehbinder's analytical solution, as demonstrated in the figures below.
Figure 1: Comparison between the MOOSE result (squares) and the analytic expression derived by Rehbinder for the porepressure.
The associated input file:
[Mesh]
[annular]
type = AnnularMeshGenerator
nr = 40
nt = 16
rmin = 0.1
rmax = 1
dmin = 0.0
dmax = 90
growth_r = 1.1
[]
[make3D]
input = annular
type = MeshExtruderGenerator
bottom_sideset = bottom
top_sideset = top
extrusion_vector = '0 0 1'
num_layers = 1
[]
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
biot_coefficient = 1.0
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[porepressure]
[]
[temperature]
[]
[]
[BCs]
[plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'top bottom'
[]
[ymin]
type = DirichletBC
variable = disp_y
value = 0
boundary = dmin
[]
[xmin]
type = DirichletBC
variable = disp_x
value = 0
boundary = dmax
[]
[cavity_temperature]
type = DirichletBC
variable = temperature
value = 1000
boundary = rmin
[]
[cavity_porepressure]
type = DirichletBC
variable = porepressure
value = 1E6
boundary = rmin
[]
[cavity_zero_effective_stress_x]
type = Pressure
variable = disp_x
function = 1E6
boundary = rmin
use_displaced_mesh = false
[]
[cavity_zero_effective_stress_y]
type = Pressure
variable = disp_y
function = 1E6
boundary = rmin
use_displaced_mesh = false
[]
[outer_temperature]
type = DirichletBC
variable = temperature
value = 0
boundary = rmax
[]
[outer_pressure]
type = DirichletBC
variable = porepressure
value = 0
boundary = rmax
[]
[fixed_outer_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = rmax
[]
[fixed_outer_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = rmax
[]
[]
[AuxVariables]
[stress_rr]
family = MONOMIAL
order = CONSTANT
[]
[stress_pp]
family = MONOMIAL
order = CONSTANT
[]
[]
[AuxKernels]
[stress_rr]
type = RankTwoScalarAux
rank_two_tensor = stress
variable = stress_rr
scalar_type = RadialStress
point1 = '0 0 0'
point2 = '0 0 1'
[]
[stress_pp]
type = RankTwoScalarAux
rank_two_tensor = stress
variable = stress_pp
scalar_type = HoopStress
point1 = '0 0 0'
point2 = '0 0 1'
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.0
bulk_modulus = 1E12
viscosity = 1.0E-3
density0 = 1000.0
cv = 1000.0
cp = 1000.0
porepressure_coefficient = 0.0
[]
[]
[PorousFlowBasicTHM]
coupling_type = ThermoHydroMechanical
multiply_by_density = false
add_stress_aux = true
porepressure = porepressure
temperature = temperature
eigenstrain_names = thermal_contribution
gravity = '0 0 0'
fp = the_simple_fluid
[]
[Materials]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 1E10
poissons_ratio = 0.2
[]
[strain]
type = ComputeSmallStrain
eigenstrain_names = thermal_contribution
[]
[thermal_contribution]
type = ComputeThermalExpansionEigenstrain
temperature = temperature
thermal_expansion_coeff = 1E-6
eigenstrain_name = thermal_contribution
stress_free_temperature = 0.0
[]
[stress]
type = ComputeLinearElasticStress
[]
[porosity]
type = PorousFlowPorosityConst # only the initial value of this is ever used
porosity = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
solid_bulk_compliance = 1E-10
fluid_bulk_modulus = 1E12
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[thermal_expansion]
type = PorousFlowConstantThermalExpansionCoefficient
fluid_coefficient = 1E-6
drained_coefficient = 1E-6
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '1E6 0 0 0 1E6 0 0 0 1E6'
[]
[]
[VectorPostprocessors]
[P]
type = LineValueSampler
start_point = '0.1 0 0'
end_point = '1.0 0 0'
num_points = 10
sort_by = x
variable = porepressure
[]
[T]
type = LineValueSampler
start_point = '0.1 0 0'
end_point = '1.0 0 0'
num_points = 10
sort_by = x
variable = temperature
[]
[U]
type = LineValueSampler
start_point = '0.1 0 0'
end_point = '1.0 0 0'
num_points = 10
sort_by = x
variable = disp_x
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_rtol'
petsc_options_value = 'gmres asm lu 1E-8'
[]
[]
[Executioner]
type = Steady
solve_type = Newton
[]
[Outputs]
file_base = fixed_outer
execute_on = timestep_end
csv = true
[]
Figure 2: Comparison between the MOOSE result (squares) and the analytic expression derived by Rehbinder for the porepressure.
The associated input file:
[Mesh]
[annular]
type = AnnularMeshGenerator
nr = 40
nt = 16
rmin = 0.1
rmax = 1
dmin = 0.0
dmax = 90
growth_r = 1.1
[]
[make3D]
input = annular
type = MeshExtruderGenerator
bottom_sideset = bottom
top_sideset = top
extrusion_vector = '0 0 1'
num_layers = 1
[]
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
biot_coefficient = 1.0
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[porepressure]
[]
[temperature]
[]
[]
[BCs]
[plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'top bottom'
[]
[ymin]
type = DirichletBC
variable = disp_y
value = 0
boundary = dmin
[]
[xmin]
type = DirichletBC
variable = disp_x
value = 0
boundary = dmax
[]
[cavity_temperature]
type = DirichletBC
variable = temperature
value = 1000
boundary = rmin
[]
[cavity_porepressure]
type = DirichletBC
variable = porepressure
value = 1E6
boundary = rmin
[]
[cavity_zero_effective_stress_x]
type = Pressure
variable = disp_x
function = 1E6
boundary = rmin
use_displaced_mesh = false
[]
[cavity_zero_effective_stress_y]
type = Pressure
variable = disp_y
function = 1E6
boundary = rmin
use_displaced_mesh = false
[]
[outer_temperature]
type = DirichletBC
variable = temperature
value = 0
boundary = rmax
[]
[outer_pressure]
type = DirichletBC
variable = porepressure
value = 0
boundary = rmax
[]
[fixed_outer_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = rmax
[]
[fixed_outer_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = rmax
[]
[]
[AuxVariables]
[stress_rr]
family = MONOMIAL
order = CONSTANT
[]
[stress_pp]
family = MONOMIAL
order = CONSTANT
[]
[]
[AuxKernels]
[stress_rr]
type = RankTwoScalarAux
rank_two_tensor = stress
variable = stress_rr
scalar_type = RadialStress
point1 = '0 0 0'
point2 = '0 0 1'
[]
[stress_pp]
type = RankTwoScalarAux
rank_two_tensor = stress
variable = stress_pp
scalar_type = HoopStress
point1 = '0 0 0'
point2 = '0 0 1'
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.0
bulk_modulus = 1E12
viscosity = 1.0E-3
density0 = 1000.0
cv = 1000.0
cp = 1000.0
porepressure_coefficient = 0.0
[]
[]
[PorousFlowBasicTHM]
coupling_type = ThermoHydroMechanical
multiply_by_density = false
add_stress_aux = true
porepressure = porepressure
temperature = temperature
eigenstrain_names = thermal_contribution
gravity = '0 0 0'
fp = the_simple_fluid
[]
[Materials]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 1E10
poissons_ratio = 0.2
[]
[strain]
type = ComputeSmallStrain
eigenstrain_names = thermal_contribution
[]
[thermal_contribution]
type = ComputeThermalExpansionEigenstrain
temperature = temperature
thermal_expansion_coeff = 1E-6
eigenstrain_name = thermal_contribution
stress_free_temperature = 0.0
[]
[stress]
type = ComputeLinearElasticStress
[]
[porosity]
type = PorousFlowPorosityConst # only the initial value of this is ever used
porosity = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
solid_bulk_compliance = 1E-10
fluid_bulk_modulus = 1E12
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[thermal_expansion]
type = PorousFlowConstantThermalExpansionCoefficient
fluid_coefficient = 1E-6
drained_coefficient = 1E-6
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '1E6 0 0 0 1E6 0 0 0 1E6'
[]
[]
[VectorPostprocessors]
[P]
type = LineValueSampler
start_point = '0.1 0 0'
end_point = '1.0 0 0'
num_points = 10
sort_by = x
variable = porepressure
[]
[T]
type = LineValueSampler
start_point = '0.1 0 0'
end_point = '1.0 0 0'
num_points = 10
sort_by = x
variable = temperature
[]
[U]
type = LineValueSampler
start_point = '0.1 0 0'
end_point = '1.0 0 0'
num_points = 10
sort_by = x
variable = disp_x
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_rtol'
petsc_options_value = 'gmres asm lu 1E-8'
[]
[]
[Executioner]
type = Steady
solve_type = Newton
[]
[Outputs]
file_base = fixed_outer
execute_on = timestep_end
csv = true
[]
Figure 3: Comparison between the MOOSE result (squares) and the analytic expression derived by Rehbinder for the temperature.
The associated input file for the fixed case:
[Mesh]
[annular]
type = AnnularMeshGenerator
nr = 40
nt = 16
rmin = 0.1
rmax = 1
dmin = 0.0
dmax = 90
growth_r = 1.1
[]
[make3D]
input = annular
type = MeshExtruderGenerator
bottom_sideset = bottom
top_sideset = top
extrusion_vector = '0 0 1'
num_layers = 1
[]
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
biot_coefficient = 1.0
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[porepressure]
[]
[temperature]
[]
[]
[BCs]
[plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'top bottom'
[]
[ymin]
type = DirichletBC
variable = disp_y
value = 0
boundary = dmin
[]
[xmin]
type = DirichletBC
variable = disp_x
value = 0
boundary = dmax
[]
[cavity_temperature]
type = DirichletBC
variable = temperature
value = 1000
boundary = rmin
[]
[cavity_porepressure]
type = DirichletBC
variable = porepressure
value = 1E6
boundary = rmin
[]
[cavity_zero_effective_stress_x]
type = Pressure
variable = disp_x
function = 1E6
boundary = rmin
use_displaced_mesh = false
[]
[cavity_zero_effective_stress_y]
type = Pressure
variable = disp_y
function = 1E6
boundary = rmin
use_displaced_mesh = false
[]
[outer_temperature]
type = DirichletBC
variable = temperature
value = 0
boundary = rmax
[]
[outer_pressure]
type = DirichletBC
variable = porepressure
value = 0
boundary = rmax
[]
[fixed_outer_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = rmax
[]
[fixed_outer_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = rmax
[]
[]
[AuxVariables]
[stress_rr]
family = MONOMIAL
order = CONSTANT
[]
[stress_pp]
family = MONOMIAL
order = CONSTANT
[]
[]
[AuxKernels]
[stress_rr]
type = RankTwoScalarAux
rank_two_tensor = stress
variable = stress_rr
scalar_type = RadialStress
point1 = '0 0 0'
point2 = '0 0 1'
[]
[stress_pp]
type = RankTwoScalarAux
rank_two_tensor = stress
variable = stress_pp
scalar_type = HoopStress
point1 = '0 0 0'
point2 = '0 0 1'
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.0
bulk_modulus = 1E12
viscosity = 1.0E-3
density0 = 1000.0
cv = 1000.0
cp = 1000.0
porepressure_coefficient = 0.0
[]
[]
[PorousFlowBasicTHM]
coupling_type = ThermoHydroMechanical
multiply_by_density = false
add_stress_aux = true
porepressure = porepressure
temperature = temperature
eigenstrain_names = thermal_contribution
gravity = '0 0 0'
fp = the_simple_fluid
[]
[Materials]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 1E10
poissons_ratio = 0.2
[]
[strain]
type = ComputeSmallStrain
eigenstrain_names = thermal_contribution
[]
[thermal_contribution]
type = ComputeThermalExpansionEigenstrain
temperature = temperature
thermal_expansion_coeff = 1E-6
eigenstrain_name = thermal_contribution
stress_free_temperature = 0.0
[]
[stress]
type = ComputeLinearElasticStress
[]
[porosity]
type = PorousFlowPorosityConst # only the initial value of this is ever used
porosity = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
solid_bulk_compliance = 1E-10
fluid_bulk_modulus = 1E12
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[thermal_expansion]
type = PorousFlowConstantThermalExpansionCoefficient
fluid_coefficient = 1E-6
drained_coefficient = 1E-6
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '1E6 0 0 0 1E6 0 0 0 1E6'
[]
[]
[VectorPostprocessors]
[P]
type = LineValueSampler
start_point = '0.1 0 0'
end_point = '1.0 0 0'
num_points = 10
sort_by = x
variable = porepressure
[]
[T]
type = LineValueSampler
start_point = '0.1 0 0'
end_point = '1.0 0 0'
num_points = 10
sort_by = x
variable = temperature
[]
[U]
type = LineValueSampler
start_point = '0.1 0 0'
end_point = '1.0 0 0'
num_points = 10
sort_by = x
variable = disp_x
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_rtol'
petsc_options_value = 'gmres asm lu 1E-8'
[]
[]
[Executioner]
type = Steady
solve_type = Newton
[]
[Outputs]
file_base = fixed_outer
execute_on = timestep_end
csv = true
[]
The associated input file for the free case:
[Mesh]
[annular]
type = AnnularMeshGenerator
nr = 40
nt = 16
rmin = 0.1
rmax = 1
dmin = 0.0
dmax = 90
growth_r = 1.1
[]
[make3D]
input = annular
type = MeshExtruderGenerator
bottom_sideset = bottom
top_sideset = top
extrusion_vector = '0 0 1'
num_layers = 1
[]
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
biot_coefficient = 1.0
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[porepressure]
[]
[temperature]
[]
[]
[BCs]
# sideset 1 = outer
# sideset 2 = cavity
# sideset 3 = ymin
# sideset 4 = xmin
[plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'top bottom'
[]
[ymin]
type = DirichletBC
variable = disp_y
value = 0
boundary = dmin
[]
[xmin]
type = DirichletBC
variable = disp_x
value = 0
boundary = dmax
[]
[cavity_temperature]
type = DirichletBC
variable = temperature
value = 1000
boundary = rmin
[]
[cavity_porepressure]
type = DirichletBC
variable = porepressure
value = 1E6
boundary = rmin
[]
[cavity_zero_effective_stress_x]
type = Pressure
variable = disp_x
function = 1E6
boundary = rmin
use_displaced_mesh = false
[]
[cavity_zero_effective_stress_y]
type = Pressure
variable = disp_y
function = 1E6
boundary = rmin
use_displaced_mesh = false
[]
[outer_temperature]
type = DirichletBC
variable = temperature
value = 0
boundary = rmax
[]
[outer_pressure]
type = DirichletBC
variable = porepressure
value = 0
boundary = rmax
[]
[]
[AuxVariables]
[stress_rr]
family = MONOMIAL
order = CONSTANT
[]
[stress_pp]
family = MONOMIAL
order = CONSTANT
[]
[]
[AuxKernels]
[stress_rr]
type = RankTwoScalarAux
rank_two_tensor = stress
variable = stress_rr
scalar_type = RadialStress
point1 = '0 0 0'
point2 = '0 0 1'
[]
[stress_pp]
type = RankTwoScalarAux
rank_two_tensor = stress
variable = stress_pp
scalar_type = HoopStress
point1 = '0 0 0'
point2 = '0 0 1'
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.0
bulk_modulus = 1E12
viscosity = 1.0E-3
density0 = 1000.0
cv = 1000.0
cp = 1000.0
porepressure_coefficient = 0.0
[]
[]
[PorousFlowBasicTHM]
coupling_type = ThermoHydroMechanical
multiply_by_density = false
add_stress_aux = true
porepressure = porepressure
temperature = temperature
eigenstrain_names = thermal_contribution
gravity = '0 0 0'
fp = the_simple_fluid
[]
[Materials]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 1E10
poissons_ratio = 0.2
[]
[strain]
type = ComputeSmallStrain
eigenstrain_names = thermal_contribution
[]
[thermal_contribution]
type = ComputeThermalExpansionEigenstrain
temperature = temperature
thermal_expansion_coeff = 1E-6
eigenstrain_name = thermal_contribution
stress_free_temperature = 0.0
[]
[stress]
type = ComputeLinearElasticStress
[]
[porosity]
type = PorousFlowPorosityConst # only the initial value of this is ever used
porosity = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
solid_bulk_compliance = 1E-10
fluid_bulk_modulus = 1E12
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[thermal_expansion]
type = PorousFlowConstantThermalExpansionCoefficient
fluid_coefficient = 1E-6
drained_coefficient = 1E-6
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '1E6 0 0 0 1E6 0 0 0 1E6'
[]
[]
[VectorPostprocessors]
[P]
type = LineValueSampler
start_point = '0.1 0 0'
end_point = '1.0 0 0'
num_points = 10
sort_by = x
variable = porepressure
[]
[T]
type = LineValueSampler
start_point = '0.1 0 0'
end_point = '1.0 0 0'
num_points = 10
sort_by = x
variable = temperature
[]
[U]
type = LineValueSampler
start_point = '0.1 0 0'
end_point = '1.0 0 0'
num_points = 10
sort_by = x
variable = disp_x
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_rtol'
petsc_options_value = 'gmres asm lu 1E-8'
[]
[]
[Executioner]
type = Steady
solve_type = Newton
[]
[Outputs]
file_base = free_outer
execute_on = timestep_end
csv = true
[]
References
- G. Rehbinder.
Analytical solutions of stationary coupled thermo-hydro-mechanical problems.
International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 32(5):453–463, 1995.
Thermo-Hydro-Mechanical Coupling in Rock Mechanics.
URL: https://www.sciencedirect.com/science/article/pii/014890629500035F, doi:https://doi.org/10.1016/0148-9062(95)00035-F.[BibTeX]