Poroelasticity test descriptions
Introduction
The PorousFlow module includes the ability to couple fluid flow to solid mechanics, and thus includes poroelasticity, which is the theory of a fully-saturated single-phase fluid with constant bulk density and constant viscosity coupled to small-strain isotropic elasticity.
There is one important difference between the poroelasticity and most of PorousFlow, however. The time-derivative terms of poroelasticity are (1) where is the Biot modulus: (2) is the fluid porepressure, is the Biot coefficient, is the volumetric strain, is the porosity, is the solid (drained) bulk modulus, and is the fluid bulk modulus. Evidently from Eq. (2), the Biot modulus, , should evolve with time as the porosity evolves. Indeed, the terms in Eq. (1) are derived from the continuity equation using the evolution of (and ). This derivation is discussed further here. However, in the standard analytical solutions of poroelasticity theory, the Biot modulus, is considered fixed.
The PorousFlow module allows porosity to vary with fluid porepressure and volumetric strain, so usually the Biot modulus would vary too, causing differences with the analytical solutions of poroelasticity. Therefore, PorousFlow offers a porosity relationship that evolves porosity in such a way as to keep fixed. This is called PorousFlowPorosityHMBiotModulus.
PorousFlow is also built with finite strains in mind, whereas poroelasticity is not. Therefore, in comparisons with solutions from poroelasticity theory, either the strain should be kept small, or the various finite-strain switches in PorousFlow should be turned off (they are all on by default).
Simple tests
Volumetric expansion due to increasing porepressure
The porepressure within a fully-saturated sample is increased: Zero mechanical pressure is applied to the sample's exterior, so that no Neumann BCs are needed on the sample. No fluid flow occurs since the porepressure is increased uniformly throughout the sample.
The input file:
# Apply an increasing porepressure, with zero mechanical forces,
# and observe the corresponding volumetric expansion
#
# P = t
# With the Biot coefficient being 0.3, the effective stresses should be
# stress_xx = stress_yy = stress_zz = 0.3t
# With bulk modulus = 1 then should have
# vol_strain = strain_xx + strain_yy + strain_zz = 0.3t.
# I use a single element lying 0<=x<=1, 0<=y<=1 and 0<=z<=1, and
# fix the left, bottom and back boundaries appropriately,
# so at the point x=y=z=1, the displacements should be
# disp_x = disp_y = disp_z = 0.3t/3 (small strain physics is used)
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = 0
xmax = 1
ymin = 0
ymax = 1
zmin = 0
zmax = 1
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
block = 0
PorousFlowDictator = dictator
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[p]
[]
[]
[BCs]
[p]
type = FunctionDirichletBC
boundary = 'bottom top'
variable = p
function = t
[]
[xmin]
type = DirichletBC
boundary = left
variable = disp_x
value = 0
[]
[ymin]
type = DirichletBC
boundary = bottom
variable = disp_y
value = 0
[]
[zmin]
type = DirichletBC
boundary = back
variable = disp_z
value = 0
[]
[]
[Kernels]
[p_does_not_really_diffuse]
type = Diffusion
variable = p
[]
[TensorMechanics]
displacements = 'disp_x disp_y disp_z'
[]
[poro_x]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.3
variable = disp_x
component = 0
[]
[poro_y]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.3
variable = disp_y
component = 1
[]
[poro_z]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.3
variable = disp_z
component = 2
[]
[]
[AuxVariables]
[stress_xx]
order = CONSTANT
family = MONOMIAL
[]
[stress_xy]
order = CONSTANT
family = MONOMIAL
[]
[stress_xz]
order = CONSTANT
family = MONOMIAL
[]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[stress_yz]
order = CONSTANT
family = MONOMIAL
[]
[stress_zz]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[stress_xx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx
index_i = 0
index_j = 0
[]
[stress_xy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xy
index_i = 0
index_j = 1
[]
[stress_xz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xz
index_i = 0
index_j = 2
[]
[stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[]
[stress_yz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yz
index_i = 1
index_j = 2
[]
[stress_zz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zz
index_i = 2
index_j = 2
[]
[]
[Postprocessors]
[corner_x]
type = PointValue
point = '1 1 1'
variable = disp_x
[]
[corner_y]
type = PointValue
point = '1 1 1'
variable = disp_y
[]
[corner_z]
type = PointValue
point = '1 1 1'
variable = disp_z
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'p'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[elasticity_tensor]
type = ComputeElasticityTensor
# bulk modulus = 1, poisson ratio = 0.2
C_ijkl = '0.5 0.75'
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeSmallStrain
displacements = 'disp_x disp_y disp_z'
[]
[stress]
type = ComputeLinearElasticStress
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = p
capillary_pressure = pc
[]
[p_eff]
type = PorousFlowEffectiveFluidPressure
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_atol -ksp_rtol'
petsc_options_value = 'gmres bjacobi 1E-10 1E-10 10 1E-15 1E-10'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
dt = 0.1
end_time = 1
[]
[Outputs]
file_base = vol_expansion
exodus = true
[]
The effective stresses should then evolve as , and the volumetric strain . MOOSE produces this result correctly.
Undrained oedometer test
A cubic single-element fully-saturated sample has roller BCs applied to its sides and bottom. All the sample's boundaries are impermeable. A downwards (normal) displacement, , is applied to its top, and the rise in porepressure and effective stress is observed. (Here denotes the direction normal to the top face.) There is no fluid flow in the single element.
# An undrained oedometer test on a saturated poroelastic sample.
#
# The sample is a single unit element, with roller BCs on the sides
# and bottom. A constant displacement is applied to the top: disp_z = -0.01*t.
# There is no fluid flow.
#
# Under these conditions
# porepressure = -(Fluid bulk modulus)*log(1 - 0.01t)
# stress_xx = (bulk - 2*shear/3)*disp_z/L (remember this is effective stress)
# stress_zz = (bulk + 4*shear/3)*disp_z/L (remember this is effective stress)
# where L is the height of the sample (L=1 in this test)
#
# Parameters:
# Bulk modulus = 2
# Shear modulus = 1.5
# fluid bulk modulus = 1
#
# Desired output:
# zdisp = -0.01*t
# p0 = 1*log(1-0.01t)
# stress_xx = stress_yy = -0.01*t
# stress_zz = -0.04*t
#
# Regarding the "log" - it just comes from conserving fluid mass
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[porepressure]
[]
[]
[BCs]
[confinex]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[confiney]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom top'
[]
[basefixed]
type = DirichletBC
variable = disp_z
value = 0
boundary = back
[]
[top_velocity]
type = FunctionDirichletBC
variable = disp_z
function = -0.01*t
boundary = front
[]
[]
[Kernels]
[grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[]
[grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[]
[grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[]
[poro_vol_exp]
type = PorousFlowMassVolumetricExpansion
variable = porepressure
fluid_component = 0
[]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = porepressure
[]
[]
[AuxVariables]
[stress_xx]
order = CONSTANT
family = MONOMIAL
[]
[stress_xy]
order = CONSTANT
family = MONOMIAL
[]
[stress_xz]
order = CONSTANT
family = MONOMIAL
[]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[stress_yz]
order = CONSTANT
family = MONOMIAL
[]
[stress_zz]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[stress_xx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx
index_i = 0
index_j = 0
[]
[stress_xy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xy
index_i = 0
index_j = 1
[]
[stress_xz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xz
index_i = 0
index_j = 2
[]
[stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[]
[stress_yz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yz
index_i = 1
index_j = 2
[]
[stress_zz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zz
index_i = 2
index_j = 2
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 1
density0 = 1
thermal_expansion = 0
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '1 1.5'
# bulk modulus is lambda + 2*mu/3 = 1 + 2*1.5/3 = 2
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[vol_strain]
type = PorousFlowVolumetricStrain
[]
[eff_fluid_pressure]
type = PorousFlowEffectiveFluidPressure
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = porepressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure disp_x disp_y disp_z'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.8
alpha = 1
[]
[]
[Postprocessors]
[fluid_mass]
type = PorousFlowFluidMass
fluid_component = 0
execute_on = 'initial timestep_end'
[]
[p0]
type = PointValue
outputs = csv
point = '0 0 0'
variable = porepressure
[]
[zdisp]
type = PointValue
outputs = csv
point = '0 0 0.5'
variable = disp_z
[]
[stress_xx]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_xx
[]
[stress_yy]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_yy
[]
[stress_zz]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_zz
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-14 1E-8 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 10
dt = 1
[]
[Outputs]
execute_on = 'timestep_end'
file_base = undrained_oedometer
[csv]
type = CSV
[]
[]
Under these conditions, assuming constant porosity, and denoting the height ( length) of the sample by : MOOSE produces these results correctly.
Porepressure generation of a confined sample
A single-element fully-saturated sample is constrained on all sides and its boundaries are impermeable. Fluid is pumped into the sample via a source (kg.s.m) and the rise in porepressure is observed.
The input file:
# Same as pp_generation.i, but using an Action
#
# A sample is constrained on all sides and its boundaries are
# also impermeable. Fluid is pumped into the sample via a
# volumetric source (ie kg/second per cubic meter), and the
# rise in porepressure is observed.
#
# Source = s (units = kg/m^3/second)
#
# Expect:
# fluid_mass = mass0 + s*t
# stress = 0 (remember this is effective stress)
# Porepressure = fluid_bulk*log(fluid_mass_density/density_P0), where fluid_mass_density = fluid_mass/porosity
# porosity = biot+(phi0-biot)*exp(pp(biot-1)/solid_bulk)
#
# Parameters:
# Biot coefficient = 0.3
# Phi0 = 0.1
# Solid Bulk modulus = 2
# fluid_bulk = 13
# density_P0 = 1
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[porepressure]
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.0
bulk_modulus = 13.0
viscosity = 1.0
density0 = 1.0
[]
[]
[PorousFlowFullySaturated]
coupling_type = HydroMechanical
displacements = 'disp_x disp_y disp_z'
porepressure = porepressure
biot_coefficient = 0.3
gravity = '0 0 0'
fp = the_simple_fluid
stabilization = none # not needed: there is no flow
save_component_rate_in = nodal_kg_per_s
[]
[BCs]
[confinex]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[confiney]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom top'
[]
[confinez]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back front'
[]
[]
[Kernels]
[source]
type = BodyForce
function = 0.1
variable = porepressure
[]
[]
[AuxVariables]
[nodal_kg_per_s]
[]
[porosity]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[porosity]
type = PorousFlowPropertyAux
variable = porosity
property = porosity
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '1 1.5'
# bulk modulus is lambda + 2*mu/3 = 1 + 2*1.5/3 = 2
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[porosity]
type = PorousFlowPorosity
fluid = true
mechanical = true
porosity_zero = 0.1
biot_coefficient = 0.3
solid_bulk = 2
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1 0 0 0 1 0 0 0 1' # unimportant
[]
[]
[Functions]
[porosity_analytic]
type = ParsedFunction
expression = 'biot+(phi0-biot)*exp(pp*(biot-1)/bulk)'
symbol_names = 'biot phi0 pp bulk'
symbol_values = '0.3 0.1 p0 2'
[]
[]
[Postprocessors]
[nodal_kg_per_s]
type = PointValue
point = ' 0 0 0'
variable = nodal_kg_per_s
outputs = csv
[]
[fluid_mass]
type = PorousFlowFluidMass
fluid_component = 0
execute_on = 'initial timestep_end'
[]
[porosity]
type = PointValue
outputs = 'console csv'
point = '0 0 0'
variable = porosity
[]
[p0]
type = PointValue
outputs = csv
point = '0 0 0'
variable = porepressure
[]
[porosity_analytic]
type = FunctionValuePostprocessor
function = porosity_analytic
[]
[zdisp]
type = PointValue
outputs = csv
point = '0 0 0.5'
variable = disp_z
[]
[stress_xx]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_xx
[]
[stress_yy]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_yy
[]
[stress_zz]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_zz
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-14 1E-10 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 10
dt = 1
[]
[Outputs]
execute_on = 'timestep_end'
file_base = pp_generation_fullysat_action
csv = true
[]
Denoting the strength of the source by (units are s), the expected result is Here is the solid bulk modulus. MOOSE produces this result correctly.
Porepressure generation of an unconfined sample
A single-element fully-saturated sample is constrained on all sides, except its top. All its boundaries are impermeable. Fluid is pumped into the sample via a source (kg.s.m) and the rise in the top surface, the porepressure, and the stress are observed.
The input file:
# Identical to pp_generation_unconfined_fullysat_volume.i but using an Action
#
# A sample is constrained on all sides, except its top
# and its boundaries are
# also impermeable. Fluid is pumped into the sample via a
# volumetric source (ie m^3/second per cubic meter), and the
# rise in the top surface, porepressure, and stress are observed.
#
# In the standard poromechanics scenario, the Biot Modulus is held
# fixed and the source has units 1/s. Then the expected result
# is
# strain_zz = disp_z = BiotCoefficient*BiotModulus*s*t/((bulk + 4*shear/3) + BiotCoefficient^2*BiotModulus)
# porepressure = BiotModulus*(s*t - BiotCoefficient*strain_zz)
# stress_xx = (bulk - 2*shear/3)*strain_zz (remember this is effective stress)
# stress_zz = (bulk + 4*shear/3)*strain_zz (remember this is effective stress)
#
# In standard porous_flow, everything is based on mass, eg the source has
# units kg/s/m^3. This is discussed in the other pp_generation_unconfined
# models. In this test, we use the FullySaturated Kernel and set
# multiply_by_density = false
# meaning the fluid Kernel has units of volume, and the source, s, has units 1/time
#
# The ratios are:
# stress_xx/strain_zz = (bulk - 2*shear/3) = 1 (for the parameters used here)
# stress_zz/strain_zz = (bulk + 4*shear/3) = 4 (for the parameters used here)
# porepressure/strain_zz = 13.3333333 (for the parameters used here)
#
# Expect
# disp_z = 0.3*10*s*t/((2 + 4*1.5/3) + 0.3^2*10) = 0.612245*s*t
# porepressure = 10*(s*t - 0.3*0.612245*s*t) = 8.163265*s*t
# stress_xx = (2 - 2*1.5/3)*0.612245*s*t = 0.612245*s*t
# stress_zz = (2 + 4*shear/3)*0.612245*s*t = 2.44898*s*t
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[porepressure]
[]
[]
[BCs]
[confinex]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[confiney]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom top'
[]
[confinez]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back'
[]
[]
[Kernels]
[source]
type = BodyForce
function = 0.1
variable = porepressure
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.0
bulk_modulus = 3.3333333333
viscosity = 1.0
density0 = 1.0
[]
[]
[PorousFlowBasicTHM]
coupling_type = HydroMechanical
displacements = 'disp_x disp_y disp_z'
multiply_by_density = false
porepressure = porepressure
biot_coefficient = 0.3
gravity = '0 0 0'
fp = the_simple_fluid
save_component_rate_in = nodal_m3_per_s
[]
[Materials]
[elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '1 1.5'
# bulk modulus is lambda + 2*mu/3 = 1 + 2*1.5/3 = 2
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeSmallStrain
displacements = 'disp_x disp_y disp_z'
[]
[stress]
type = ComputeLinearElasticStress
[]
[porosity]
type = PorousFlowPorosityConst # the "const" is irrelevant here: all that uses Porosity is the BiotModulus, which just uses the initial value of porosity
porosity = 0.1
PorousFlowDictator = dictator
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
PorousFlowDictator = dictator
biot_coefficient = 0.3
fluid_bulk_modulus = 3.3333333333
solid_bulk_compliance = 0.5
[]
[permeability_irrelevant]
type = PorousFlowPermeabilityConst
PorousFlowDictator = dictator
permeability = '1.5 0 0 0 1.5 0 0 0 1.5'
[]
[]
[AuxVariables]
[nodal_m3_per_s]
[]
[]
[Postprocessors]
[nodal_m3_per_s]
type = PointValue
outputs = csv
point = '0 0 0'
variable = nodal_m3_per_s
[]
[p0]
type = PointValue
outputs = csv
point = '0 0 0'
variable = porepressure
[]
[zdisp]
type = PointValue
outputs = csv
point = '0 0 0.5'
variable = disp_z
[]
[stress_xx]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_xx
[]
[stress_yy]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_yy
[]
[stress_zz]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_zz
[]
[stress_xx_over_strain]
type = FunctionValuePostprocessor
function = stress_xx_over_strain_fcn
outputs = csv
[]
[stress_zz_over_strain]
type = FunctionValuePostprocessor
function = stress_zz_over_strain_fcn
outputs = csv
[]
[p_over_strain]
type = FunctionValuePostprocessor
function = p_over_strain_fcn
outputs = csv
[]
[]
[Functions]
[stress_xx_over_strain_fcn]
type = ParsedFunction
expression = a/b
symbol_names = 'a b'
symbol_values = 'stress_xx zdisp'
[]
[stress_zz_over_strain_fcn]
type = ParsedFunction
expression = a/b
symbol_names = 'a b'
symbol_values = 'stress_zz zdisp'
[]
[p_over_strain_fcn]
type = ParsedFunction
expression = a/b
symbol_names = 'a b'
symbol_values = 'p0 zdisp'
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-14 1E-10 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 10
dt = 1
[]
[Outputs]
execute_on = 'timestep_end'
file_base = pp_generation_unconfined_basicthm
[csv]
type = CSV
[]
[]
Regardless of the evolution of porosity, the following ratios result where is the undrained bulk modulus, the shear modulus, the Biot coefficient, and is the initial Biot modulus. MOOSE produces these results correctly.
However, if the Biot modulus, , is held fixed as the porosity evolves, and the source is with being a constant volumetric source (m.s.m) then The input file for this is:
# This file uses a PorousFlowFullySaturated Action. The equivalent non-Action input file is pp_generation_unconfined_constM.i
#
# A sample is constrained on all sides, except its top
# and its boundaries are
# also impermeable. Fluid is pumped into the sample via a
# volumetric source (ie kg/second per cubic meter), and the
# rise in the top surface, porepressure, and stress are observed.
#
# In the standard poromechanics scenario, the Biot Modulus is held
# fixed and the source, s, has units m^3/second/m^3. Then the expected result
# is
# strain_zz = disp_z = BiotCoefficient*BiotModulus*s*t/((bulk + 4*shear/3) + BiotCoefficient^2*BiotModulus)
# porepressure = BiotModulus*(s*t - BiotCoefficient*strain_zz)
# stress_xx = (bulk - 2*shear/3)*strain_zz (remember this is effective stress)
# stress_zz = (bulk + 4*shear/3)*strain_zz (remember this is effective stress)
#
# In porous_flow, however, the source has units kg/second/m^3. The ratios remain
# fixed:
# stress_xx/strain_zz = (bulk - 2*shear/3) = 1 (for the parameters used here)
# stress_zz/strain_zz = (bulk + 4*shear/3) = 4 (for the parameters used here)
# porepressure/strain_zz = 13.3333333 (for the parameters used here)
#
# Expect
# disp_z = 0.3*10*s*t/((2 + 4*1.5/3) + 0.3^2*10) = 0.612245*s*t
# porepressure = 10*(s*t - 0.3*0.612245*s*t) = 8.163265*s*t
# stress_xx = (2 - 2*1.5/3)*0.612245*s*t = 0.612245*s*t
# stress_zz = (2 + 4*shear/3)*0.612245*s*t = 2.44898*s*t
# The relationship between the constant poroelastic source
# s (m^3/second/m^3) and the PorousFlow source, S (kg/second/m^3) is
# S = fluid_density * s = s * exp(porepressure/fluid_bulk)
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[porepressure]
[]
[]
[BCs]
[confinex]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[confiney]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom top'
[]
[confinez]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back'
[]
[]
[PorousFlowFullySaturated]
porepressure = porepressure
biot_coefficient = 0.3
coupling_type = HydroMechanical
displacements = 'disp_x disp_y disp_z'
gravity = '0 0 0'
fp = simple_fluid
stabilization = Full
[]
[Kernels]
[source]
type = BodyForce
function = '0.1*exp(8.163265306*0.1*t/3.3333333333)'
variable = porepressure
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 3.3333333333
density0 = 1
thermal_expansion = 0
viscosity = 1
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '1 1.5'
# bulk modulus is lambda + 2*mu/3 = 1 + 2*1.5/3 = 2
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeSmallStrain
displacements = 'disp_x disp_y disp_z'
[]
[stress]
type = ComputeLinearElasticStress
[]
[porosity]
type = PorousFlowPorosityHMBiotModulus
porosity_zero = 0.1
biot_coefficient = 0.3
solid_bulk = 2
constant_fluid_bulk_modulus = 3.3333333333
constant_biot_modulus = 10.0
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1 0 0 0 1 0 0 0 1' # unimportant
[]
[]
[Postprocessors]
[p0]
type = PointValue
outputs = csv
point = '0 0 0'
variable = porepressure
[]
[zdisp]
type = PointValue
outputs = csv
point = '0 0 0.5'
variable = disp_z
[]
[stress_xx]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_xx
[]
[stress_yy]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_yy
[]
[stress_zz]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_zz
[]
[]
[Functions]
[stress_xx_over_strain_fcn]
type = ParsedFunction
expression = a/b
symbol_names = 'a b'
symbol_values = 'stress_xx zdisp'
[]
[stress_zz_over_strain_fcn]
type = ParsedFunction
expression = a/b
symbol_names = 'a b'
symbol_values = 'stress_zz zdisp'
[]
[p_over_strain_fcn]
type = ParsedFunction
expression = a/b
symbol_names = 'a b'
symbol_values = 'p0 zdisp'
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-14 1E-10 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 10
dt = 1
[]
[Outputs]
execute_on = 'timestep_end'
file_base = pp_generation_unconfined_constM_action
[csv]
type = CSV
[]
[]
MOOSE produces these results correctly.
Terzaghi consolidation of a drained medium
A saturated sample sits in a bath of water. It is constrained on its sides and bottom. Its sides and bottom are also impermeable. Initially it is unstressed (, at ). A normal stress, , is applied to the sample's top. The sample compresses instantaneously due to the instantaneous application of , and then slowly compresses further as water is squeezed out from the sample's top surface.
This is a classic test. See, for example, Section 2.2 of the online manuscript: Arnold Verruijt "Theory and Problems of Poroelasticity" Delft University of Technology 2013. But note that the "sigma" in that paper is the negative of the stress in PorousFlow. Denote the sample's height (-length) by . Define This is the porepressure that results from the instantaneous application of : MOOSE calculates this correctly. The solution for porepressure is In this equation, is the "consolidation coefficient": , where the permeability tensor is . The so-called degree-of-consolidation is defined by where is the vertical displacement of the top surface (downwards is positive), and is the instantaneous displacement due to the instantaneous application of , and is the final displacement. This has solution
There are a few different variants of this test in the test suite, that use the PorousFlow action system, constant Biot modulus, constant porosity, the fully-saturated version of PorousFlow, etc. You are encouraged to explore these to get a feeling for how the differences impact the output. The input file for that matches the theoretical setup exactly is:
# Terzaghi's problem of consolodation of a drained medium
#
# A saturated soil sample sits in a bath of water.
# It is constrained on its sides, and bottom.
# Its sides and bottom are also impermeable.
# Initially it is unstressed.
# A normal stress, q, is applied to the soil's top.
# The soil then slowly compresses as water is squeezed
# out from the sample from its top (the top BC for
# the porepressure is porepressure = 0).
#
# See, for example. Section 2.2 of the online manuscript
# Arnold Verruijt "Theory and Problems of Poroelasticity" Delft University of Technology 2013
# but note that the "sigma" in that paper is the negative
# of the stress in TensorMechanics
#
# Here are the problem's parameters, and their values:
# Soil height. h = 10
# Soil's Lame lambda. la = 2
# Soil's Lame mu, which is also the Soil's shear modulus. mu = 3
# Soil bulk modulus. K = la + 2*mu/3 = 4
# Soil confined compressibility. m = 1/(K + 4mu/3) = 0.125
# Soil bulk compliance. 1/K = 0.25
# Fluid bulk modulus. Kf = 8
# Fluid bulk compliance. 1/Kf = 0.125
# Fluid mobility (soil permeability/fluid viscosity). k = 1.5
# Soil initial porosity. phi0 = 0.1
# Biot coefficient. alpha = 0.6
# Soil initial storativity, which is the reciprocal of the initial Biot modulus. S = phi0/Kf + (alpha - phi0)(1 - alpha)/K = 0.0625
# Consolidation coefficient. c = k/(S + alpha^2 m) = 13.95348837
# Normal stress on top. q = 1
# Initial porepressure, resulting from instantaneous application of q, assuming corresponding instantaneous increase of porepressure (Note that this is calculated by MOOSE: we only need it for the analytical solution). p0 = alpha*m*q/(S + alpha^2 m) = 0.69767442
# Initial vertical displacement (down is positive), resulting from instantaneous application of q (Note this is calculated by MOOSE: we only need it for the analytical solution). uz0 = q*m*h*S/(S + alpha^2 m)
# Final vertical displacement (down in positive) (Note this is calculated by MOOSE: we only need it for the analytical solution). uzinf = q*m*h
#
# The solution for porepressure is
# P = 4*p0/\pi \sum_{k=1}^{\infty} \frac{(-1)^{k-1}}{2k-1} \cos ((2k-1)\pi z/(2h)) \exp(-(2k-1)^2 \pi^2 ct/(4 h^2))
# This series converges very slowly for ct/h^2 small, so in that domain
# P = p0 erf( (1-(z/h))/(2 \sqrt(ct/h^2)) )
#
# The degree of consolidation is defined as
# U = (uz - uz0)/(uzinf - uz0)
# where uz0 and uzinf are defined above, and
# uz = the vertical displacement of the top (down is positive)
# U = 1 - (8/\pi^2)\sum_{k=1}^{\infty} \frac{1}{(2k-1)^2} \exp(-(2k-1)^2 \pi^2 ct/(4 h^2))
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 10
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = 0
zmax = 10
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure disp_x disp_y disp_z'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.8
alpha = 1
[]
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[porepressure]
[]
[]
[BCs]
[confinex]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[confiney]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom top'
[]
[basefixed]
type = DirichletBC
variable = disp_z
value = 0
boundary = back
[]
[topdrained]
type = DirichletBC
variable = porepressure
value = 0
boundary = front
[]
[topload]
type = NeumannBC
variable = disp_z
value = -1
boundary = front
[]
[]
[Kernels]
[grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[]
[grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[]
[grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[]
[poro_x]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_x
component = 0
[]
[poro_y]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_y
component = 1
[]
[poro_z]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
component = 2
variable = disp_z
[]
[poro_vol_exp]
type = PorousFlowMassVolumetricExpansion
variable = porepressure
fluid_component = 0
[]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = porepressure
[]
[flux]
type = PorousFlowAdvectiveFlux
variable = porepressure
gravity = '0 0 0'
fluid_component = 0
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 8
density0 = 1
thermal_expansion = 0
viscosity = 0.96
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '2 3'
# bulk modulus is lambda + 2*mu/3 = 2 + 2*3/3 = 4
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[eff_fluid_pressure]
type = PorousFlowEffectiveFluidPressure
[]
[vol_strain]
type = PorousFlowVolumetricStrain
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = porepressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityHMBiotModulus
porosity_zero = 0.1
biot_coefficient = 0.6
solid_bulk = 4
constant_fluid_bulk_modulus = 8
constant_biot_modulus = 16
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1.5 0 0 0 1.5 0 0 0 1.5'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 0 # unimportant in this fully-saturated situation
phase = 0
[]
[]
[Postprocessors]
[p0]
type = PointValue
outputs = csv
point = '0 0 0'
variable = porepressure
use_displaced_mesh = false
[]
[p1]
type = PointValue
outputs = csv
point = '0 0 1'
variable = porepressure
use_displaced_mesh = false
[]
[p2]
type = PointValue
outputs = csv
point = '0 0 2'
variable = porepressure
use_displaced_mesh = false
[]
[p3]
type = PointValue
outputs = csv
point = '0 0 3'
variable = porepressure
use_displaced_mesh = false
[]
[p4]
type = PointValue
outputs = csv
point = '0 0 4'
variable = porepressure
use_displaced_mesh = false
[]
[p5]
type = PointValue
outputs = csv
point = '0 0 5'
variable = porepressure
use_displaced_mesh = false
[]
[p6]
type = PointValue
outputs = csv
point = '0 0 6'
variable = porepressure
use_displaced_mesh = false
[]
[p7]
type = PointValue
outputs = csv
point = '0 0 7'
variable = porepressure
use_displaced_mesh = false
[]
[p8]
type = PointValue
outputs = csv
point = '0 0 8'
variable = porepressure
use_displaced_mesh = false
[]
[p9]
type = PointValue
outputs = csv
point = '0 0 9'
variable = porepressure
use_displaced_mesh = false
[]
[p99]
type = PointValue
outputs = csv
point = '0 0 10'
variable = porepressure
use_displaced_mesh = false
[]
[zdisp]
type = PointValue
outputs = csv
point = '0 0 10'
variable = disp_z
use_displaced_mesh = false
[]
[dt]
type = FunctionValuePostprocessor
outputs = console
function = if(0.5*t<0.1,0.5*t,0.1)
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 10
[TimeStepper]
type = PostprocessorDT
postprocessor = dt
dt = 0.0001
[]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = terzaghi_constM
[csv]
type = CSV
[]
[]
MOOSE produces the expected results correctly, as may be seen from Figure 1 and Figure 2.
Figure 1: Degree of consolidation in the Terzaghi experiment.
Figure 2: Porepressure at various times in the Terzaghi experiment.
Mandel's consolidation of a drained medium
A sample's dimensions are and , and it is in plane strain (no displacement). It is squashed with constant normal force by impermeable, frictionless plattens on its top and bottom surfaces (at ). Fluid is allowed to leak out from its sides (at ), but all other surfaces are impermeable. This is called Mandel's problem and it is shown graphically in Figure 3
Figure 3: The setup of the Mandel experiment: a force squashes a porous material with impermeable plattens. This causes fluid to seep from the material.
The interesting feature of this problem (apart from that it can be solved analytically) is that the porepressure in the sample's center increases for a short time after application of the force. This is because the leakage of the fluid from the sample's sides causes an apparent softening of the material near those sides. This means stress concentrates towards the sample's center which causes an increase in porepressure. Of course, eventually the fluid totally drains from the sample, and the porepressure is zero. As the fluid drains from the sample's sides the plattens move slowly towards each other.
The solution for porepressure and displacements is given in Cheng and Detournay (1988). The solution involves rather lengthy infinite series, so I will not write it here.
As is common in the literature, this is simulated by considering the quarter-sample, and , with impermeable, roller BCs at and and . Porepressure is fixed at zero on . Porepressure and displacement are initialised to zero. Then the top () is moved downwards with prescribed velocity, so that the total force that is induced by this downwards velocity is fixed. The velocity is worked out by solving Mandel's problem analytically, and the total force is monitored in the simulation to check that it indeed remains constant.
The simulations in the PorousFlow test suite use 10 elements in the direction and 1 in the direction. Four types of simulation are run:
HM. This uses standard PorousFlow Materials and Kernels, in particular it uses the "HM" porosity relationship. This is not expected to agree perfectly with the analytical solutions because the solutions assume constant Biot modulus.
# Mandel's problem of consolodation of a drained medium
#
# A sample is in plane strain.
# -a <= x <= a
# -b <= y <= b
# It is squashed with constant force by impermeable, frictionless plattens on its top and bottom surfaces (at y=+/-b)
# Fluid is allowed to leak out from its sides (at x=+/-a)
# The porepressure within the sample is monitored.
#
# As is common in the literature, this is simulated by
# considering the quarter-sample, 0<=x<=a and 0<=y<=b, with
# impermeable, roller BCs at x=0 and y=0 and y=b.
# Porepressure is fixed at zero on x=a.
# Porepressure and displacement are initialised to zero.
# Then the top (y=b) is moved downwards with prescribed velocity,
# so that the total force that is inducing this downwards velocity
# is fixed. The velocity is worked out by solving Mandel's problem
# analytically, and the total force is monitored in the simulation
# to check that it indeed remains constant.
#
# Here are the problem's parameters, and their values:
# Soil width. a = 1
# Soil height. b = 0.1
# Soil's Lame lambda. la = 0.5
# Soil's Lame mu, which is also the Soil's shear modulus. mu = G = 0.75
# Soil bulk modulus. K = la + 2*mu/3 = 1
# Drained Poisson ratio. nu = (3K - 2G)/(6K + 2G) = 0.2
# Soil bulk compliance. 1/K = 1
# Fluid bulk modulus. Kf = 8
# Fluid bulk compliance. 1/Kf = 0.125
# Soil initial porosity. phi0 = 0.1
# Biot coefficient. alpha = 0.6
# Biot modulus. M = 1/(phi0/Kf + (alpha - phi0)(1 - alpha)/K) = 4.705882
# Undrained bulk modulus. Ku = K + alpha^2*M = 2.694118
# Undrained Poisson ratio. nuu = (3Ku - 2G)/(6Ku + 2G) = 0.372627
# Skempton coefficient. B = alpha*M/Ku = 1.048035
# Fluid mobility (soil permeability/fluid viscosity). k = 1.5
# Consolidation coefficient. c = 2*k*B^2*G*(1-nu)*(1+nuu)^2/9/(1-nuu)/(nuu-nu) = 3.821656
# Normal stress on top. F = 1
#
# The solution for porepressure and displacements is given in
# AHD Cheng and E Detournay "A direct boundary element method for plane strain poroelasticity" International Journal of Numerical and Analytical Methods in Geomechanics 12 (1988) 551-572.
# The solution involves complicated infinite series, so I shall not write it here
#
# FINAL NOTE: The above solution assumes constant Biot Modulus.
# In porous_flow this is not true. Therefore the solution is
# a little different than in the paper. This test was therefore
# validated against MOOSE's poromechanics, which can choose either
# a constant Biot Modulus (which has been shown to agree with
# the analytic solution), or a non-constant Biot Modulus (which
# gives the same results as porous_flow).
[Mesh]
type = GeneratedMesh
dim = 3
nx = 10
ny = 1
nz = 1
xmin = 0
xmax = 1
ymin = 0
ymax = 0.1
zmin = 0
zmax = 1
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure disp_x disp_y disp_z'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.8
alpha = 1e-5
[]
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[porepressure]
[]
[]
[BCs]
[roller_xmin]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left'
[]
[roller_ymin]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[]
[plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back front'
[]
[xmax_drained]
type = DirichletBC
variable = porepressure
value = 0
boundary = right
[]
[top_velocity]
type = FunctionDirichletBC
variable = disp_y
function = top_velocity
boundary = top
[]
[]
[Functions]
[top_velocity]
type = PiecewiseLinear
x = '0 0.002 0.006 0.014 0.03 0.046 0.062 0.078 0.094 0.11 0.126 0.142 0.158 0.174 0.19 0.206 0.222 0.238 0.254 0.27 0.286 0.302 0.318 0.334 0.35 0.366 0.382 0.398 0.414 0.43 0.446 0.462 0.478 0.494 0.51 0.526 0.542 0.558 0.574 0.59 0.606 0.622 0.638 0.654 0.67 0.686 0.702'
y = '-0.041824842 -0.042730269 -0.043412712 -0.04428867 -0.045509181 -0.04645965 -0.047268246 -0.047974749 -0.048597109 -0.0491467 -0.049632388 -0.050061697 -0.050441198 -0.050776675 -0.051073238 -0.0513354 -0.051567152 -0.051772022 -0.051953128 -0.052113227 -0.052254754 -0.052379865 -0.052490464 -0.052588233 -0.052674662 -0.052751065 -0.052818606 -0.052878312 -0.052931093 -0.052977751 -0.053018997 -0.053055459 -0.053087691 -0.053116185 -0.053141373 -0.05316364 -0.053183324 -0.053200724 -0.053216106 -0.053229704 -0.053241725 -0.053252351 -0.053261745 -0.053270049 -0.053277389 -0.053283879 -0.053289615'
[]
[]
[AuxVariables]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[tot_force]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[]
[tot_force]
type = ParsedAux
coupled_variables = 'stress_yy porepressure'
execute_on = timestep_end
variable = tot_force
expression = '-stress_yy+0.6*porepressure'
[]
[]
[Kernels]
[grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[]
[grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[]
[grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[]
[poro_x]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_x
component = 0
[]
[poro_y]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_y
component = 1
[]
[poro_z]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
component = 2
variable = disp_z
[]
[poro_vol_exp]
type = PorousFlowMassVolumetricExpansion
variable = porepressure
fluid_component = 0
[]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = porepressure
[]
[flux]
type = PorousFlowAdvectiveFlux
variable = porepressure
gravity = '0 0 0'
fluid_component = 0
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 8
density0 = 1
thermal_expansion = 0
viscosity = 1
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '0.5 0.75'
# bulk modulus is lambda + 2*mu/3 = 0.5 + 2*0.75/3 = 1
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[eff_fluid_pressure]
type = PorousFlowEffectiveFluidPressure
[]
[vol_strain]
type = PorousFlowVolumetricStrain
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = porepressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosity
fluid = true
mechanical = true
ensure_positive = false
porosity_zero = 0.1
biot_coefficient = 0.6
solid_bulk = 1
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1.5 0 0 0 1.5 0 0 0 1.5'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 0 # unimportant in this fully-saturated situation
phase = 0
[]
[]
[Postprocessors]
[p0]
type = PointValue
outputs = csv
point = '0.0 0 0'
variable = porepressure
[]
[p1]
type = PointValue
outputs = csv
point = '0.1 0 0'
variable = porepressure
[]
[p2]
type = PointValue
outputs = csv
point = '0.2 0 0'
variable = porepressure
[]
[p3]
type = PointValue
outputs = csv
point = '0.3 0 0'
variable = porepressure
[]
[p4]
type = PointValue
outputs = csv
point = '0.4 0 0'
variable = porepressure
[]
[p5]
type = PointValue
outputs = csv
point = '0.5 0 0'
variable = porepressure
[]
[p6]
type = PointValue
outputs = csv
point = '0.6 0 0'
variable = porepressure
[]
[p7]
type = PointValue
outputs = csv
point = '0.7 0 0'
variable = porepressure
[]
[p8]
type = PointValue
outputs = csv
point = '0.8 0 0'
variable = porepressure
[]
[p9]
type = PointValue
outputs = csv
point = '0.9 0 0'
variable = porepressure
[]
[p99]
type = PointValue
outputs = csv
point = '1 0 0'
variable = porepressure
[]
[xdisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_x
[]
[ydisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_y
[]
[total_downwards_force]
type = ElementAverageValue
outputs = csv
variable = tot_force
[]
[dt]
type = FunctionValuePostprocessor
outputs = console
function = if(0.15*t<0.01,0.15*t,0.01)
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'gmres asm lu 1E-14 1E-10 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 0.7
[TimeStepper]
type = PostprocessorDT
postprocessor = dt
dt = 0.001
[]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = mandel
[csv]
time_step_interval = 3
type = CSV
[]
[]
constM. This is identical to the HM case, save that it uses a porosity evolution law that keeps the Biot modulus fixed. It is therefore expected to agree with the analytical solutions.
# Mandel's problem of consolodation of a drained medium
#
# A sample is in plane strain.
# -a <= x <= a
# -b <= y <= b
# It is squashed with constant force by impermeable, frictionless plattens on its top and bottom surfaces (at y=+/-b)
# Fluid is allowed to leak out from its sides (at x=+/-a)
# The porepressure within the sample is monitored.
#
# As is common in the literature, this is simulated by
# considering the quarter-sample, 0<=x<=a and 0<=y<=b, with
# impermeable, roller BCs at x=0 and y=0 and y=b.
# Porepressure is fixed at zero on x=a.
# Porepressure and displacement are initialised to zero.
# Then the top (y=b) is moved downwards with prescribed velocity,
# so that the total force that is inducing this downwards velocity
# is fixed. The velocity is worked out by solving Mandel's problem
# analytically, and the total force is monitored in the simulation
# to check that it indeed remains constant.
#
# Here are the problem's parameters, and their values:
# Soil width. a = 1
# Soil height. b = 0.1
# Soil's Lame lambda. la = 0.5
# Soil's Lame mu, which is also the Soil's shear modulus. mu = G = 0.75
# Soil bulk modulus. K = la + 2*mu/3 = 1
# Drained Poisson ratio. nu = (3K - 2G)/(6K + 2G) = 0.2
# Soil bulk compliance. 1/K = 1
# Fluid bulk modulus. Kf = 8
# Fluid bulk compliance. 1/Kf = 0.125
# Soil initial porosity. phi0 = 0.1
# Biot coefficient. alpha = 0.6
# Biot modulus. M = 1/(phi0/Kf + (alpha - phi0)(1 - alpha)/K) = 4.705882
# Undrained bulk modulus. Ku = K + alpha^2*M = 2.694118
# Undrained Poisson ratio. nuu = (3Ku - 2G)/(6Ku + 2G) = 0.372627
# Skempton coefficient. B = alpha*M/Ku = 1.048035
# Fluid mobility (soil permeability/fluid viscosity). k = 1.5
# Consolidation coefficient. c = 2*k*B^2*G*(1-nu)*(1+nuu)^2/9/(1-nuu)/(nuu-nu) = 3.821656
# Normal stress on top. F = 1
#
# The solution for porepressure and displacements is given in
# AHD Cheng and E Detournay "A direct boundary element method for plane strain poroelasticity" International Journal of Numerical and Analytical Methods in Geomechanics 12 (1988) 551-572.
# The solution involves complicated infinite series, so I shall not write it here
[Mesh]
type = GeneratedMesh
dim = 3
nx = 10
ny = 1
nz = 1
xmin = 0
xmax = 1
ymin = 0
ymax = 0.1
zmin = 0
zmax = 1
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure disp_x disp_y disp_z'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.8
alpha = 1e-5
[]
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[porepressure]
[]
[]
[BCs]
[roller_xmin]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left'
[]
[roller_ymin]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[]
[plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back front'
[]
[xmax_drained]
type = DirichletBC
variable = porepressure
value = 0
boundary = right
[]
[top_velocity]
type = FunctionDirichletBC
variable = disp_y
function = top_velocity
boundary = top
[]
[]
[Functions]
[top_velocity]
type = PiecewiseLinear
x = '0 0.002 0.006 0.014 0.03 0.046 0.062 0.078 0.094 0.11 0.126 0.142 0.158 0.174 0.19 0.206 0.222 0.238 0.254 0.27 0.286 0.302 0.318 0.334 0.35 0.366 0.382 0.398 0.414 0.43 0.446 0.462 0.478 0.494 0.51 0.526 0.542 0.558 0.574 0.59 0.606 0.622 0.638 0.654 0.67 0.686 0.702'
y = '-0.041824842 -0.042730269 -0.043412712 -0.04428867 -0.045509181 -0.04645965 -0.047268246 -0.047974749 -0.048597109 -0.0491467 -0.049632388 -0.050061697 -0.050441198 -0.050776675 -0.051073238 -0.0513354 -0.051567152 -0.051772022 -0.051953128 -0.052113227 -0.052254754 -0.052379865 -0.052490464 -0.052588233 -0.052674662 -0.052751065 -0.052818606 -0.052878312 -0.052931093 -0.052977751 -0.053018997 -0.053055459 -0.053087691 -0.053116185 -0.053141373 -0.05316364 -0.053183324 -0.053200724 -0.053216106 -0.053229704 -0.053241725 -0.053252351 -0.053261745 -0.053270049 -0.053277389 -0.053283879 -0.053289615'
[]
[]
[AuxVariables]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[tot_force]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[]
[tot_force]
type = ParsedAux
coupled_variables = 'stress_yy porepressure'
execute_on = timestep_end
variable = tot_force
expression = '-stress_yy+0.6*porepressure'
[]
[]
[Kernels]
[grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[]
[grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[]
[grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[]
[poro_x]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_x
component = 0
[]
[poro_y]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_y
component = 1
[]
[poro_z]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
component = 2
variable = disp_z
[]
[poro_vol_exp]
type = PorousFlowMassVolumetricExpansion
variable = porepressure
fluid_component = 0
[]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = porepressure
[]
[flux]
type = PorousFlowAdvectiveFlux
variable = porepressure
gravity = '0 0 0'
fluid_component = 0
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 8
density0 = 1
thermal_expansion = 0
viscosity = 1
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '0.5 0.75'
# bulk modulus is lambda + 2*mu/3 = 0.5 + 2*0.75/3 = 1
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[eff_fluid_pressure]
type = PorousFlowEffectiveFluidPressure
[]
[vol_strain]
type = PorousFlowVolumetricStrain
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = porepressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityHMBiotModulus
porosity_zero = 0.1
biot_coefficient = 0.6
solid_bulk = 1
constant_fluid_bulk_modulus = 8
constant_biot_modulus = 4.7058823529
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1.5 0 0 0 1.5 0 0 0 1.5'
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 0 # unimportant in this fully-saturated situation
phase = 0
[]
[]
[Postprocessors]
[p0]
type = PointValue
outputs = csv
point = '0.0 0 0'
variable = porepressure
[]
[p1]
type = PointValue
outputs = csv
point = '0.1 0 0'
variable = porepressure
[]
[p2]
type = PointValue
outputs = csv
point = '0.2 0 0'
variable = porepressure
[]
[p3]
type = PointValue
outputs = csv
point = '0.3 0 0'
variable = porepressure
[]
[p4]
type = PointValue
outputs = csv
point = '0.4 0 0'
variable = porepressure
[]
[p5]
type = PointValue
outputs = csv
point = '0.5 0 0'
variable = porepressure
[]
[p6]
type = PointValue
outputs = csv
point = '0.6 0 0'
variable = porepressure
[]
[p7]
type = PointValue
outputs = csv
point = '0.7 0 0'
variable = porepressure
[]
[p8]
type = PointValue
outputs = csv
point = '0.8 0 0'
variable = porepressure
[]
[p9]
type = PointValue
outputs = csv
point = '0.9 0 0'
variable = porepressure
[]
[p99]
type = PointValue
outputs = csv
point = '1 0 0'
variable = porepressure
[]
[xdisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_x
[]
[ydisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_y
[]
[total_downwards_force]
type = ElementAverageValue
outputs = csv
variable = tot_force
[]
[dt]
type = FunctionValuePostprocessor
outputs = console
function = if(0.15*t<0.01,0.15*t,0.01)
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'gmres asm lu 1E-14 1E-10 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 0.7
[TimeStepper]
type = PostprocessorDT
postprocessor = dt
dt = 0.001
[]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = mandel_constM
[csv]
time_step_interval = 3
type = CSV
[]
[]
FullSat. This uses the FullySaturated versions of the fluid mass time derivative and the fluid flux. In this case the Biot modulus is kept fixed, so it is expected to agree with the analytical solutions.
# Mandel's problem of consolodation of a drained medium
# Using the FullySaturatedDarcyBase and FullySaturatedMassTimeDerivative kernels
#
# A sample is in plane strain.
# -a <= x <= a
# -b <= y <= b
# It is squashed with constant force by impermeable, frictionless plattens on its top and bottom surfaces (at y=+/-b)
# Fluid is allowed to leak out from its sides (at x=+/-a)
# The porepressure within the sample is monitored.
#
# As is common in the literature, this is simulated by
# considering the quarter-sample, 0<=x<=a and 0<=y<=b, with
# impermeable, roller BCs at x=0 and y=0 and y=b.
# Porepressure is fixed at zero on x=a.
# Porepressure and displacement are initialised to zero.
# Then the top (y=b) is moved downwards with prescribed velocity,
# so that the total force that is inducing this downwards velocity
# is fixed. The velocity is worked out by solving Mandel's problem
# analytically, and the total force is monitored in the simulation
# to check that it indeed remains constant.
#
# Here are the problem's parameters, and their values:
# Soil width. a = 1
# Soil height. b = 0.1
# Soil's Lame lambda. la = 0.5
# Soil's Lame mu, which is also the Soil's shear modulus. mu = G = 0.75
# Soil bulk modulus. K = la + 2*mu/3 = 1
# Drained Poisson ratio. nu = (3K - 2G)/(6K + 2G) = 0.2
# Soil bulk compliance. 1/K = 1
# Fluid bulk modulus. Kf = 8
# Fluid bulk compliance. 1/Kf = 0.125
# Soil initial porosity. phi0 = 0.1
# Biot coefficient. alpha = 0.6
# Biot modulus. M = 1/(phi0/Kf + (alpha - phi0)(1 - alpha)/K) = 4.705882
# Undrained bulk modulus. Ku = K + alpha^2*M = 2.694118
# Undrained Poisson ratio. nuu = (3Ku - 2G)/(6Ku + 2G) = 0.372627
# Skempton coefficient. B = alpha*M/Ku = 1.048035
# Fluid mobility (soil permeability/fluid viscosity). k = 1.5
# Consolidation coefficient. c = 2*k*B^2*G*(1-nu)*(1+nuu)^2/9/(1-nuu)/(nuu-nu) = 3.821656
# Normal stress on top. F = 1
#
# The solution for porepressure and displacements is given in
# AHD Cheng and E Detournay "A direct boundary element method for plane strain poroelasticity" International Journal of Numerical and Analytical Methods in Geomechanics 12 (1988) 551-572.
# The solution involves complicated infinite series, so I shall not write it here
[Mesh]
type = GeneratedMesh
dim = 3
nx = 10
ny = 1
nz = 1
xmin = 0
xmax = 1
ymin = 0
ymax = 0.1
zmin = 0
zmax = 1
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure disp_x disp_y disp_z'
number_fluid_phases = 1
number_fluid_components = 1
[]
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[porepressure]
[]
[]
[BCs]
[roller_xmin]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left'
[]
[roller_ymin]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[]
[plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back front'
[]
[xmax_drained]
type = DirichletBC
variable = porepressure
value = 0
boundary = right
[]
[top_velocity]
type = FunctionDirichletBC
variable = disp_y
function = top_velocity
boundary = top
[]
[]
[Functions]
[top_velocity]
type = PiecewiseLinear
x = '0 0.002 0.006 0.014 0.03 0.046 0.062 0.078 0.094 0.11 0.126 0.142 0.158 0.174 0.19 0.206 0.222 0.238 0.254 0.27 0.286 0.302 0.318 0.334 0.35 0.366 0.382 0.398 0.414 0.43 0.446 0.462 0.478 0.494 0.51 0.526 0.542 0.558 0.574 0.59 0.606 0.622 0.638 0.654 0.67 0.686 0.702'
y = '-0.041824842 -0.042730269 -0.043412712 -0.04428867 -0.045509181 -0.04645965 -0.047268246 -0.047974749 -0.048597109 -0.0491467 -0.049632388 -0.050061697 -0.050441198 -0.050776675 -0.051073238 -0.0513354 -0.051567152 -0.051772022 -0.051953128 -0.052113227 -0.052254754 -0.052379865 -0.052490464 -0.052588233 -0.052674662 -0.052751065 -0.052818606 -0.052878312 -0.052931093 -0.052977751 -0.053018997 -0.053055459 -0.053087691 -0.053116185 -0.053141373 -0.05316364 -0.053183324 -0.053200724 -0.053216106 -0.053229704 -0.053241725 -0.053252351 -0.053261745 -0.053270049 -0.053277389 -0.053283879 -0.053289615'
[]
[]
[AuxVariables]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[tot_force]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[]
[tot_force]
type = ParsedAux
coupled_variables = 'stress_yy porepressure'
execute_on = timestep_end
variable = tot_force
expression = '-stress_yy+0.6*porepressure'
[]
[]
[Kernels]
[grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[]
[grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[]
[grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[]
[poro_x]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_x
component = 0
[]
[poro_y]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
variable = disp_y
component = 1
[]
[poro_z]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.6
component = 2
variable = disp_z
[]
[mass0]
type = PorousFlowFullySaturatedMassTimeDerivative
biot_coefficient = 0.6
coupling_type = HydroMechanical
variable = porepressure
[]
[flux]
type = PorousFlowFullySaturatedDarcyBase
variable = porepressure
gravity = '0 0 0'
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 8
density0 = 1
thermal_expansion = 0
viscosity = 1
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
[]
[elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '0.5 0.75'
# bulk modulus is lambda + 2*mu/3 = 0.5 + 2*0.75/3 = 1
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[eff_fluid_pressure_qp]
type = PorousFlowEffectiveFluidPressure
[]
[vol_strain]
type = PorousFlowVolumetricStrain
[]
[ppss]
type = PorousFlow1PhaseFullySaturated
porepressure = porepressure
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid_qp]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst # only the initial value of this is ever used
porosity = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.6
solid_bulk_compliance = 1
fluid_bulk_modulus = 8
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1.5 0 0 0 1.5 0 0 0 1.5'
[]
[]
[Postprocessors]
[p0]
type = PointValue
outputs = csv
point = '0.0 0 0'
variable = porepressure
[]
[p1]
type = PointValue
outputs = csv
point = '0.1 0 0'
variable = porepressure
[]
[p2]
type = PointValue
outputs = csv
point = '0.2 0 0'
variable = porepressure
[]
[p3]
type = PointValue
outputs = csv
point = '0.3 0 0'
variable = porepressure
[]
[p4]
type = PointValue
outputs = csv
point = '0.4 0 0'
variable = porepressure
[]
[p5]
type = PointValue
outputs = csv
point = '0.5 0 0'
variable = porepressure
[]
[p6]
type = PointValue
outputs = csv
point = '0.6 0 0'
variable = porepressure
[]
[p7]
type = PointValue
outputs = csv
point = '0.7 0 0'
variable = porepressure
[]
[p8]
type = PointValue
outputs = csv
point = '0.8 0 0'
variable = porepressure
[]
[p9]
type = PointValue
outputs = csv
point = '0.9 0 0'
variable = porepressure
[]
[p99]
type = PointValue
outputs = csv
point = '1 0 0'
variable = porepressure
[]
[xdisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_x
[]
[ydisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_y
[]
[total_downwards_force]
type = ElementAverageValue
outputs = csv
variable = tot_force
[]
[dt]
type = FunctionValuePostprocessor
outputs = console
function = if(0.15*t<0.01,0.15*t,0.01)
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'gmres asm lu 1E-14 1E-10 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 0.7
[TimeStepper]
type = PostprocessorDT
postprocessor = dt
dt = 0.001
[]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = mandel_fully_saturated
[csv]
time_step_interval = 3
type = CSV
[]
[]
FullSatVol. This uses the FullySaturated versions of the fluid mass time derivative and the fluid flux, and does not multiply by the fluid density. Therefore this version is identical to what is usually implemented in poro-elastic codes. It is linear and therefore converges in only one iteration. In this case the Biot modulus is kept fixed, so it is expected to agree with the analytical solutions.
# using a BasicTHM Action
#
# Mandel's problem of consolodation of a drained medium
# Using the FullySaturatedDarcyBase and FullySaturatedFullySaturatedMassTimeDerivative kernels
# with multiply_by_density = false, so that this problem becomes linear
#
# A sample is in plane strain.
# -a <= x <= a
# -b <= y <= b
# It is squashed with constant force by impermeable, frictionless plattens on its top and bottom surfaces (at y=+/-b)
# Fluid is allowed to leak out from its sides (at x=+/-a)
# The porepressure within the sample is monitored.
#
# As is common in the literature, this is simulated by
# considering the quarter-sample, 0<=x<=a and 0<=y<=b, with
# impermeable, roller BCs at x=0 and y=0 and y=b.
# Porepressure is fixed at zero on x=a.
# Porepressure and displacement are initialised to zero.
# Then the top (y=b) is moved downwards with prescribed velocity,
# so that the total force that is inducing this downwards velocity
# is fixed. The velocity is worked out by solving Mandel's problem
# analytically, and the total force is monitored in the simulation
# to check that it indeed remains constant.
#
# Here are the problem's parameters, and their values:
# Soil width. a = 1
# Soil height. b = 0.1
# Soil's Lame lambda. la = 0.5
# Soil's Lame mu, which is also the Soil's shear modulus. mu = G = 0.75
# Soil bulk modulus. K = la + 2*mu/3 = 1
# Drained Poisson ratio. nu = (3K - 2G)/(6K + 2G) = 0.2
# Soil bulk compliance. 1/K = 1
# Fluid bulk modulus. Kf = 8
# Fluid bulk compliance. 1/Kf = 0.125
# Soil initial porosity. phi0 = 0.1
# Biot coefficient. alpha = 0.6
# Biot modulus. M = 1/(phi0/Kf + (alpha - phi0)(1 - alpha)/K) = 4.705882
# Undrained bulk modulus. Ku = K + alpha^2*M = 2.694118
# Undrained Poisson ratio. nuu = (3Ku - 2G)/(6Ku + 2G) = 0.372627
# Skempton coefficient. B = alpha*M/Ku = 1.048035
# Fluid mobility (soil permeability/fluid viscosity). k = 1.5
# Consolidation coefficient. c = 2*k*B^2*G*(1-nu)*(1+nuu)^2/9/(1-nuu)/(nuu-nu) = 3.821656
# Normal stress on top. F = 1
#
# The solution for porepressure and displacements is given in
# AHD Cheng and E Detournay "A direct boundary element method for plane strain poroelasticity" International Journal of Numerical and Analytical Methods in Geomechanics 12 (1988) 551-572.
# The solution involves complicated infinite series, so I shall not write it here
[Mesh]
type = GeneratedMesh
dim = 3
nx = 10
ny = 1
nz = 1
xmin = 0
xmax = 1
ymin = 0
ymax = 0.1
zmin = 0
zmax = 1
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[porepressure]
[]
[]
[BCs]
[roller_xmin]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left'
[]
[roller_ymin]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[]
[plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back front'
[]
[xmax_drained]
type = DirichletBC
variable = porepressure
value = 0
boundary = right
[]
[top_velocity]
type = FunctionDirichletBC
variable = disp_y
function = top_velocity
boundary = top
[]
[]
[Functions]
[top_velocity]
type = PiecewiseLinear
x = '0 0.002 0.006 0.014 0.03 0.046 0.062 0.078 0.094 0.11 0.126 0.142 0.158 0.174 0.19 0.206 0.222 0.238 0.254 0.27 0.286 0.302 0.318 0.334 0.35 0.366 0.382 0.398 0.414 0.43 0.446 0.462 0.478 0.494 0.51 0.526 0.542 0.558 0.574 0.59 0.606 0.622 0.638 0.654 0.67 0.686 0.702'
y = '-0.041824842 -0.042730269 -0.043412712 -0.04428867 -0.045509181 -0.04645965 -0.047268246 -0.047974749 -0.048597109 -0.0491467 -0.049632388 -0.050061697 -0.050441198 -0.050776675 -0.051073238 -0.0513354 -0.051567152 -0.051772022 -0.051953128 -0.052113227 -0.052254754 -0.052379865 -0.052490464 -0.052588233 -0.052674662 -0.052751065 -0.052818606 -0.052878312 -0.052931093 -0.052977751 -0.053018997 -0.053055459 -0.053087691 -0.053116185 -0.053141373 -0.05316364 -0.053183324 -0.053200724 -0.053216106 -0.053229704 -0.053241725 -0.053252351 -0.053261745 -0.053270049 -0.053277389 -0.053283879 -0.053289615'
[]
[]
[AuxVariables]
[tot_force]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[tot_force]
type = ParsedAux
coupled_variables = 'stress_yy porepressure'
execute_on = timestep_end
variable = tot_force
expression = '-stress_yy+0.6*porepressure'
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.0
bulk_modulus = 8.0
viscosity = 1.0
density0 = 1.0
[]
[]
[PorousFlowBasicTHM]
coupling_type = HydroMechanical
displacements = 'disp_x disp_y disp_z'
multiply_by_density = false
porepressure = porepressure
biot_coefficient = 0.6
gravity = '0 0 0'
fp = the_simple_fluid
[]
[Materials]
[elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '0.5 0.75'
# bulk modulus is lambda + 2*mu/3 = 0.5 + 2*0.75/3 = 1
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[porosity]
type = PorousFlowPorosityConst # only the initial value of this is ever used
porosity = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.6
solid_bulk_compliance = 1
fluid_bulk_modulus = 8
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1.5 0 0 0 1.5 0 0 0 1.5'
[]
[]
[Postprocessors]
[p0]
type = PointValue
outputs = csv
point = '0.0 0 0'
variable = porepressure
[]
[p1]
type = PointValue
outputs = csv
point = '0.1 0 0'
variable = porepressure
[]
[p2]
type = PointValue
outputs = csv
point = '0.2 0 0'
variable = porepressure
[]
[p3]
type = PointValue
outputs = csv
point = '0.3 0 0'
variable = porepressure
[]
[p4]
type = PointValue
outputs = csv
point = '0.4 0 0'
variable = porepressure
[]
[p5]
type = PointValue
outputs = csv
point = '0.5 0 0'
variable = porepressure
[]
[p6]
type = PointValue
outputs = csv
point = '0.6 0 0'
variable = porepressure
[]
[p7]
type = PointValue
outputs = csv
point = '0.7 0 0'
variable = porepressure
[]
[p8]
type = PointValue
outputs = csv
point = '0.8 0 0'
variable = porepressure
[]
[p9]
type = PointValue
outputs = csv
point = '0.9 0 0'
variable = porepressure
[]
[p99]
type = PointValue
outputs = csv
point = '1 0 0'
variable = porepressure
[]
[xdisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_x
[]
[ydisp]
type = PointValue
outputs = csv
point = '1 0.1 0'
variable = disp_y
[]
[total_downwards_force]
type = ElementAverageValue
outputs = csv
variable = tot_force
[]
[dt]
type = FunctionValuePostprocessor
outputs = console
function = if(0.15*t<0.01,0.15*t,0.01)
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'gmres asm lu 1E-14 1E-10 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = 0
end_time = 0.7
[TimeStepper]
type = PostprocessorDT
postprocessor = dt
dt = 0.001
[]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = mandel_basicthm
[csv]
time_step_interval = 3
type = CSV
[]
[]
Of course there are minor discrepancies between the last three and the analytical solution that are brought about through spatial and temporal discretisation errors. The figures below present the results.
Figure 4: The vertical displacement of the platten as a function of time.
Figure 5: The horizontal displacement of the platten's top-right corner.
Figure 6: The porepressure at various points in the sample in the HM model ("a" is equal to unity).
Figure 7: The porepressure at various points in the sample in the FullSatVol model ("a" is equal to unity).
Figure 8: The total downwards force on the platten as a function of time. This should be unity.
References
- Alexander H.-D. Cheng and Emmanuel Detournay.
A direct boundary element method for plane strain poroelasticity.
International Journal for Numerical and Analytical Methods in Geomechanics, 12(5):551–572, 1988.
URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nag.1610120508, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/nag.1610120508, doi:10.1002/nag.1610120508.[BibTeX]