- PorousFlowDictatorThe UserObject that holds the list of PorousFlow variable names
C++ Type:UserObjectName
Controllable:No
Description:The UserObject that holds the list of PorousFlow variable names
PorousFlowConstantBiotModulus
Computes the Biot Modulus, which is assumed to be constant for all time. Sometimes 1 / BiotModulus is called storativity
The Biot Modulus is defined to be (1) where is the porosity, is the fluid bulk modulus, is the Biot coefficient and is the drained bulk modulus of the porous skeleton.
This Material is designed to work with the PorousFlowFullySaturatedMassTimeDerivative Kernel.
commentnote
This quantity is computed during the initial stage of the simulation and is kept fixed thereafter.
Input Parameters
- at_nodesFalseEvaluate Material properties at nodes instead of quadpoints
Default:False
C++ Type:bool
Controllable:No
Description:Evaluate Material properties at nodes instead of quadpoints
- biot_coefficient1Biot coefficient
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Biot coefficient
- blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
- boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
- computeTrueWhen false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
Default:True
C++ Type:bool
Controllable:No
Description:When false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
- constant_onNONEWhen ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
Default:NONE
C++ Type:MooseEnum
Options:NONE, ELEMENT, SUBDOMAIN
Controllable:No
Description:When ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
- declare_suffixAn optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
- fluid_bulk_modulus2e+09Fluid bulk modulus
Default:2e+09
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Fluid bulk modulus
- solid_bulk_compliance0Reciprocal of the drained bulk modulus of the porous skeleton. If strain = C * stress, then solid_bulk_compliance = de_ij de_kl C_ijkl. If the grain bulk modulus is Kg then 1/Kg = (1 - biot_coefficient) * solid_bulk_compliance.
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Reciprocal of the drained bulk modulus of the porous skeleton. If strain = C * stress, then solid_bulk_compliance = de_ij de_kl C_ijkl. If the grain bulk modulus is Kg then 1/Kg = (1 - biot_coefficient) * solid_bulk_compliance.
Optional Parameters
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
- implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
- search_methodnearest_node_connected_sidesChoice of search algorithm. All options begin by finding the nearest node in the primary boundary to a query point in the secondary boundary. In the default nearest_node_connected_sides algorithm, primary boundary elements are searched iff that nearest node is one of their nodes. This is fast to determine via a pregenerated node-to-elem map and is robust on conforming meshes. In the optional all_proximate_sides algorithm, primary boundary elements are searched iff they touch that nearest node, even if they are not topologically connected to it. This is more CPU-intensive but is necessary for robustness on any boundary surfaces which has disconnections (such as Flex IGA meshes) or non-conformity (such as hanging nodes in adaptively h-refined meshes).
Default:nearest_node_connected_sides
C++ Type:MooseEnum
Options:nearest_node_connected_sides, all_proximate_sides
Controllable:No
Description:Choice of search algorithm. All options begin by finding the nearest node in the primary boundary to a query point in the secondary boundary. In the default nearest_node_connected_sides algorithm, primary boundary elements are searched iff that nearest node is one of their nodes. This is fast to determine via a pregenerated node-to-elem map and is robust on conforming meshes. In the optional all_proximate_sides algorithm, primary boundary elements are searched iff they touch that nearest node, even if they are not topologically connected to it. This is more CPU-intensive but is necessary for robustness on any boundary surfaces which has disconnections (such as Flex IGA meshes) or non-conformity (such as hanging nodes in adaptively h-refined meshes).
- seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
- use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Advanced Parameters
- output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector<std::string>
Controllable:No
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
- outputsnone Vector of output names where you would like to restrict the output of variables(s) associated with this object
Default:none
C++ Type:std::vector<OutputName>
Controllable:No
Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object
Outputs Parameters
- prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
- use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Material Property Retrieval Parameters
Input Files
- (modules/porous_flow/test/tests/hysteresis/except04.i)
- (modules/porous_flow/test/tests/thm_rehbinder/fixed_outer_rz.i)
- (modules/porous_flow/test/tests/hysteresis/except08.i)
- (modules/porous_flow/test/tests/hysteresis/except02.i)
- (modules/porous_flow/test/tests/hysteresis/except03.i)
- (modules/porous_flow/test/tests/dirackernels/injection_production.i)
- (modules/porous_flow/test/tests/actions/basicthm_h.i)
- (modules/porous_flow/examples/groundwater/ex01.i)
- (modules/porous_flow/examples/tutorial/03.i)
- (modules/porous_flow/test/tests/actions/basicthm_borehole.i)
- (modules/porous_flow/test/tests/poro_elasticity/mandel_fully_saturated.i)
- (modules/porous_flow/test/tests/actions/basicthm_hm.i)
- (modules/porous_flow/test/tests/pressure_pulse/pressure_pulse_1d_fully_saturated_2.i)
- (modules/porous_flow/test/tests/actions/basicthm_thm.i)
- (modules/porous_flow/test/tests/actions/basicthm_th.i)
- (modules/porous_flow/examples/tutorial/04.i)
- (modules/porous_flow/test/tests/poro_elasticity/mandel_basicthm.i)
- (modules/porous_flow/test/tests/poro_elasticity/terzaghi_basicthm.i)
- (modules/porous_flow/test/tests/poro_elasticity/pp_generation_unconfined_fully_saturated_volume.i)
- (modules/porous_flow/examples/tidal/barometric_fully_confined.i)
- (modules/porous_flow/test/tests/dirackernels/theis_rz.i)
- (modules/porous_flow/examples/tidal/earth_tide_fullsat.i)
- (modules/porous_flow/test/tests/hysteresis/except07.i)
- (modules/porous_flow/test/tests/sinks/PorousFlowPiecewiseLinearSink_BC_eg1.i)
- (modules/porous_flow/test/tests/poro_elasticity/terzaghi_fully_saturated_volume.i)
- (modules/porous_flow/examples/tutorial/01.i)
- (modules/porous_flow/test/tests/hysteresis/except09.i)
- (modules/porous_flow/examples/tidal/atm_tides_open_hole.i)
- (modules/porous_flow/test/tests/thm_rehbinder/fixed_outer.i)
- (modules/porous_flow/test/tests/poro_elasticity/pp_generation_unconfined_fully_saturated.i)
- (modules/porous_flow/test/tests/hysteresis/except05.i)
- (modules/porous_flow/examples/tutorial/05.i)
- (modules/porous_flow/test/tests/hysteresis/except01.i)
- (modules/porous_flow/test/tests/poro_elasticity/mandel_fully_saturated_volume.i)
- (modules/porous_flow/test/tests/hysteresis/except06.i)
- (modules/porous_flow/examples/tutorial/05_tabulated.i)
- (modules/porous_flow/test/tests/jacobian/mass01_fully_saturated.i)
- (modules/porous_flow/test/tests/thm_rehbinder/free_outer.i)
- (modules/porous_flow/examples/tidal/atm_tides.i)
- (modules/porous_flow/test/tests/poro_elasticity/pp_generation_unconfined_basicthm.i)
(modules/porous_flow/test/tests/hysteresis/except04.i)
# Exception testing of PorousFlowHysteresisOrder
# Incorrect: previous_turning_points not in the range [0, 1]
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 1
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[]
[PorousFlowBasicTHM]
porepressure = pp
fp = simple_fluid
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2e-7
fluid_bulk_modulus = 1e7
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[hys_order]
type = PorousFlowHysteresisOrder
initial_order = 1
previous_turning_points = -0.1
[]
[]
[Preconditioning]
[basic]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 1
[]
(modules/porous_flow/test/tests/thm_rehbinder/fixed_outer_rz.i)
# A version of fixed_outer.i that uses the RZ cylindrical coordinate system
[Mesh]
type = GeneratedMesh
dim = 2
nx = 40 # this is the r direction
ny = 1 # this is the height direction
xmin = 0.1
xmax = 1
bias_x = 1.1
ymin = 0.0
ymax = 1.0
coord_type = RZ
[]
[GlobalParams]
displacements = 'disp_r disp_z'
PorousFlowDictator = dictator
biot_coefficient = 1.0
[]
[Variables]
[disp_r]
[]
[disp_z]
[]
[porepressure]
[]
[temperature]
[]
[]
[BCs]
[plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'top bottom'
[]
[cavity_temperature]
type = DirichletBC
variable = temperature
value = 1000
boundary = left
[]
[cavity_porepressure]
type = DirichletBC
variable = porepressure
value = 1E6
boundary = left
[]
[cavity_zero_effective_stress_x]
type = Pressure
variable = disp_r
function = 1E6
boundary = left
use_displaced_mesh = false
[]
[outer_temperature]
type = DirichletBC
variable = temperature
value = 0
boundary = right
[]
[outer_pressure]
type = DirichletBC
variable = porepressure
value = 0
boundary = right
[]
[fixed_outer_disp]
type = DirichletBC
variable = disp_r
value = 0
boundary = right
[]
[]
[AuxVariables]
[stress_rr]
family = MONOMIAL
order = CONSTANT
[]
[stress_pp]
family = MONOMIAL
order = CONSTANT
[]
[]
[AuxKernels]
[stress_rr]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_rr
index_i = 0
index_j = 0
[]
[stress_pp] # hoop stress
type = RankTwoAux
rank_two_tensor = stress
variable = stress_pp
index_i = 2
index_j = 2
[]
[]
[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 = ComputeAxisymmetricRZSmallStrain
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' # note this is ordered: rr, zz, angle-angle
[]
[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' # note this is ordered: rr, zz, angle-angle
[]
[]
[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_r
[]
[]
[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_rz
execute_on = timestep_end
csv = true
[]
(modules/porous_flow/test/tests/hysteresis/except08.i)
# Exception testing of PorousFlowHysteresisOrder
# Incorrectly ordered previous_turning_points
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 1
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[]
[PorousFlowBasicTHM]
porepressure = pp
fp = simple_fluid
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2e-7
fluid_bulk_modulus = 1e7
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[hys_order]
type = PorousFlowHysteresisOrder
initial_order = 3
previous_turning_points = '0.6 0.8 0.5'
[]
[]
[Preconditioning]
[basic]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 1
[]
(modules/porous_flow/test/tests/hysteresis/except02.i)
# Exception testing of PorousFlowHysteresisOrder
# Incorrect: initial_order = 4
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 1
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[]
[PorousFlowBasicTHM]
porepressure = pp
fp = simple_fluid
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2e-7
fluid_bulk_modulus = 1e7
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[hys_order]
type = PorousFlowHysteresisOrder
initial_order = 4
[]
[]
[Preconditioning]
[basic]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 1
[]
(modules/porous_flow/test/tests/hysteresis/except03.i)
# Exception testing of PorousFlowHysteresisOrder
# Incorrect: initial_order incommensurate with previous_turning_points
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 1
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[]
[PorousFlowBasicTHM]
porepressure = pp
fp = simple_fluid
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2e-7
fluid_bulk_modulus = 1e7
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[hys_order]
type = PorousFlowHysteresisOrder
initial_order = 1
[]
[]
[Preconditioning]
[basic]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 1
[]
(modules/porous_flow/test/tests/dirackernels/injection_production.i)
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 3
nx = 10
ny = 10
nz = 1
xmin = -50
xmax = 50
ymin = -50
ymax = 50
zmin = 0
zmax = 10
[]
[central_nodes]
input = gen
type = ExtraNodesetGenerator
new_boundary = central_nodes
coord = '0 0 0; 0 0 10'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[porepressure]
initial_condition = 20E6
[]
[temperature]
initial_condition = 400
scaling = 1E-6 # fluid enthalpy is roughly 1E6
[]
[]
[BCs]
[injection_temperature]
type = DirichletBC
variable = temperature
value = 300
boundary = central_nodes
[]
[]
[DiracKernels]
[fluid_injection]
type = PorousFlowPeacemanBorehole
variable = porepressure
SumQuantityUO = injected_mass
point_file = injection.bh
function_of = pressure
fluid_phase = 0
bottom_p_or_t = 21E6
unit_weight = '0 0 0'
use_mobility = true
character = -1
[]
[fluid_production]
type = PorousFlowPeacemanBorehole
variable = porepressure
SumQuantityUO = produced_mass
point_file = production.bh
function_of = pressure
fluid_phase = 0
bottom_p_or_t = 20E6
unit_weight = '0 0 0'
use_mobility = true
character = 1
[]
[remove_heat_at_production_well]
type = PorousFlowPeacemanBorehole
variable = temperature
SumQuantityUO = produced_heat
point_file = production.bh
function_of = pressure
fluid_phase = 0
bottom_p_or_t = 20E6
unit_weight = '0 0 0'
use_mobility = true
use_enthalpy = true
character = 1
[]
[]
[UserObjects]
[injected_mass]
type = PorousFlowSumQuantity
[]
[produced_mass]
type = PorousFlowSumQuantity
[]
[produced_heat]
type = PorousFlowSumQuantity
[]
[]
[Postprocessors]
[heat_joules_extracted_this_timestep]
type = PorousFlowPlotQuantity
uo = produced_heat
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 2E-4
bulk_modulus = 2E9
viscosity = 1E-3
density0 = 1000
cv = 4000.0
cp = 4000.0
[]
[]
[PorousFlowUnsaturated]
porepressure = porepressure
temperature = temperature
coupling_type = ThermoHydro
gravity = '0 0 0'
fp = the_simple_fluid
[]
[Materials]
[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 = 2E9
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[thermal_expansion]
type = PorousFlowConstantThermalExpansionCoefficient
fluid_coefficient = 5E-6
drained_coefficient = 2E-4
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '1 0 0 0 1 0 0 0 1'
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
density = 2500.0
specific_heat_capacity = 1200.0
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 2E6
dt = 2E5
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/test/tests/actions/basicthm_h.i)
# PorousFlowBasicTHM action with coupling_type = HydroGenerator
# (no thermal or mechanical effects)
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
nx = 10
ny = 3
xmax = 10
ymax = 3
[]
[aquifer]
input = gen
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 1 0'
top_right = '10 2 0'
[]
[injection_area]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 'injection_area'
normal = '-1 0 0'
input = 'aquifer'
[]
[outflow_area]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 'outflow_area'
normal = '1 0 0'
input = 'injection_area'
[]
[rename]
type = RenameBlockGenerator
old_block = '0 1'
new_block = 'caprock aquifer'
input = 'outflow_area'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[porepressure]
initial_condition = 1e6
[]
[]
[AuxVariables]
[temperature]
initial_condition = 293
[]
[]
[PorousFlowBasicTHM]
porepressure = porepressure
temperature = temperature
coupling_type = Hydro
gravity = '0 0 0'
fp = simple_fluid
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1.5e6
boundary = injection_area
[]
[constant_outflow_porepressure]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = outflow_area
pt_vals = '0 1e9'
multipliers = '0 1e9'
flux_function = 1e-6
PT_shift = 1e6
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2e-7
fluid_bulk_modulus = 1e7
[]
[permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[permeability_caprock]
type = PorousFlowPermeabilityConst
block = caprock
permeability = '1e-15 0 0 0 1e-15 0 0 0 1e-15'
[]
[]
[Preconditioning]
[basic]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1e4
dt = 1e3
nl_abs_tol = 1e-15
nl_rel_tol = 1E-14
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/examples/groundwater/ex01.i)
# Groundwater extraction example.
# System consists of two confined aquifers separated by an aquitard
# There is a hydraulic gradient in the upper aquifer
# A well extracts water from the lower aquifer, and the impact on the upper aquifer is observed
# In the center of the model, the roof of the upper aquifer sits 70m below the local water table
[Mesh]
[basic_mesh]
type = GeneratedMeshGenerator
dim = 3
xmin = -50
xmax = 50
nx = 20
ymin = -25
ymax = 25
ny = 10
zmin = -100
zmax = -70
nz = 3
[]
[lower_aquifer]
type = SubdomainBoundingBoxGenerator
input = basic_mesh
block_id = 1
block_name = lower_aquifer
bottom_left = '-1000 -500 -100'
top_right = '1000 500 -90'
[]
[aquitard]
type = SubdomainBoundingBoxGenerator
input = lower_aquifer
block_id = 2
block_name = aquitard
bottom_left = '-1000 -500 -90'
top_right = '1000 500 -80'
[]
[upper_aquifer]
type = SubdomainBoundingBoxGenerator
input = aquitard
block_id = 3
block_name = upper_aquifer
bottom_left = '-1000 -500 -80'
top_right = '1000 500 -70'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[]
[ICs]
[pp]
type = FunctionIC
variable = pp
function = insitu_pp
[]
[]
[BCs]
[pp]
type = FunctionDirichletBC
variable = pp
function = insitu_pp
boundary = 'left right top bottom front back'
[]
[]
[Functions]
[upper_aquifer_head]
type = ParsedFunction
expression = '10 + x / 200'
[]
[lower_aquifer_head]
type = ParsedFunction
expression = '20'
[]
[insitu_head]
type = ParsedFunction
symbol_values = 'lower_aquifer_head upper_aquifer_head'
symbol_names = 'low up'
expression = 'if(z <= -90, low, if(z >= -80, up, (up * (z + 90) - low * (z + 80)) / (10.0)))'
[]
[insitu_pp]
type = ParsedFunction
symbol_values = 'insitu_head'
symbol_names = 'h'
expression = '(h - z) * 1E4'
[]
[l_rate]
type = ParsedFunction
symbol_values = 'm3_produced dt'
symbol_names = 'm3_produced dt'
expression = '1000 * m3_produced / dt'
[]
[]
[AuxVariables]
[insitu_head]
[]
[head_change]
[]
[]
[AuxKernels]
[insitu_head]
type = FunctionAux
variable = insitu_head
function = insitu_head
[]
[head_change]
type = ParsedAux
coupled_variables = 'pp insitu_head'
use_xyzt = true
expression = 'pp / 1E4 + z - insitu_head'
variable = head_change
[]
[]
[Postprocessors]
[m3_produced]
type = PorousFlowPlotQuantity
uo = volume_extracted
outputs = 'none'
[]
[dt]
type = TimestepSize
outputs = 'none'
[]
[l_per_s]
type = FunctionValuePostprocessor
function = l_rate
[]
[]
[VectorPostprocessors]
[drawdown]
type = LineValueSampler
variable = head_change
start_point = '-50 0 -75'
end_point = '50 0 -75'
num_points = 101
sort_by = x
[]
[]
[PorousFlowBasicTHM]
fp = simple_fluid
gravity = '0 0 -10'
porepressure = pp
multiply_by_density = false
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
# the following mean that density = 1000 * exp(P / 1E15) ~ 1000
thermal_expansion = 0
bulk_modulus = 1E15
[]
[]
[Materials]
[porosity_aquifers]
type = PorousFlowPorosityConst
porosity = 0.05
block = 'upper_aquifer lower_aquifer'
[]
[porosity_aquitard]
type = PorousFlowPorosityConst
porosity = 0.2
block = aquitard
[]
[biot_mod]
type = PorousFlowConstantBiotModulus
fluid_bulk_modulus = 2E9
biot_coefficient = 1.0
[]
[permeability_aquifers]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
block = 'upper_aquifer lower_aquifer'
[]
[permeability_aquitard]
type = PorousFlowPermeabilityConst
permeability = '1E-16 0 0 0 1E-16 0 0 0 1E-17'
block = aquitard
[]
[]
[DiracKernels]
[sink]
type = PorousFlowPolyLineSink
SumQuantityUO = volume_extracted
point_file = ex01.bh_lower
line_length = 10
variable = pp
# following produces a flux of 0 m^3(water)/m(borehole length)/s if porepressure = 0, and a flux of 1 m^3/m/s if porepressure = 1E9
p_or_t_vals = '0 1E9'
fluxes = '0 1'
[]
[]
[UserObjects]
[volume_extracted]
type = PorousFlowSumQuantity
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
[TimeStepper]
type = SolutionTimeAdaptiveDT
dt = 1.1E5
[]
end_time = 3.456E5 # 4 days
nl_abs_tol = 1E-13
[]
[Outputs]
[csv]
type = CSV
file_base = ex01_lower_extraction
execute_on = final
[]
[]
(modules/porous_flow/examples/tutorial/03.i)
# Darcy flow with heat advection and conduction
[Mesh]
[annular]
type = AnnularMeshGenerator
nr = 10
rmin = 1.0
rmax = 10
growth_r = 1.4
nt = 4
dmin = 0
dmax = 90
[]
[make3D]
type = MeshExtruderGenerator
extrusion_vector = '0 0 12'
num_layers = 3
bottom_sideset = 'bottom'
top_sideset = 'top'
input = annular
[]
[shift_down]
type = TransformGenerator
transform = TRANSLATE
vector_value = '0 0 -6'
input = make3D
[]
[aquifer]
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 0 -2'
top_right = '10 10 2'
input = shift_down
[]
[injection_area]
type = ParsedGenerateSideset
combinatorial_geometry = 'x*x+y*y<1.01'
included_subdomains = 1
new_sideset_name = 'injection_area'
input = 'aquifer'
[]
[rename]
type = RenameBlockGenerator
old_block = '0 1'
new_block = 'caps aquifer'
input = 'injection_area'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[porepressure]
[]
[temperature]
initial_condition = 293
scaling = 1E-8
[]
[]
[PorousFlowBasicTHM]
porepressure = porepressure
temperature = temperature
coupling_type = ThermoHydro
gravity = '0 0 0'
fp = the_simple_fluid
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1E6
boundary = injection_area
[]
[constant_injection_temperature]
type = DirichletBC
variable = temperature
value = 313
boundary = injection_area
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
viscosity = 1.0E-3
density0 = 1000.0
thermal_expansion = 0.0002
cp = 4194
cv = 4186
porepressure_coefficient = 0
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2E-7
fluid_bulk_modulus = 1E7
[]
[permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1E-14 0 0 0 1E-14 0 0 0 1E-14'
[]
[permeability_caps]
type = PorousFlowPermeabilityConst
block = caps
permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-16'
[]
[thermal_expansion]
type = PorousFlowConstantThermalExpansionCoefficient
biot_coefficient = 0.8
drained_coefficient = 0.003
fluid_coefficient = 0.0002
[]
[rock_internal_energy]
type = PorousFlowMatrixInternalEnergy
density = 2500.0
specific_heat_capacity = 1200.0
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '10 0 0 0 10 0 0 0 10'
block = 'caps aquifer'
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_abs_tol = 1E-10
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/test/tests/actions/basicthm_borehole.i)
# PorousFlowBasicTHM action with coupling_type = Hydro (no thermal or
# mechanical effects), plus a Peaceman borehole with use_mobility = true
# to test that nodal relative permeability is added by this action.
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -1
xmax = 1
ymin = -1
ymax = 1
zmin = -1
zmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[porepressure]
initial_condition = 1e7
[]
[]
[AuxVariables]
[temperature]
initial_condition = 293
[]
[]
[PorousFlowBasicTHM]
porepressure = porepressure
temperature = temperature
coupling_type = Hydro
gravity = '0 0 0'
fp = simple_fluid
multiply_by_density = true
[]
[UserObjects]
[borehole_total_outflow_mass]
type = PorousFlowSumQuantity
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
viscosity = 1e-3
density0 = 1000
thermal_expansion = 0
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.2
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 1
fluid_bulk_modulus = 2e9
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[]
[DiracKernels]
[bh]
type = PorousFlowPeacemanBorehole
variable = porepressure
SumQuantityUO = borehole_total_outflow_mass
point_file = borehole.bh
function_of = pressure
fluid_phase = 0
bottom_p_or_t = 0
unit_weight = '0 0 0'
use_mobility = true
character = 1
[]
[]
[Postprocessors]
[bh_report]
type = PorousFlowPlotQuantity
uo = borehole_total_outflow_mass
[]
[]
[Preconditioning]
[usual]
type = SMP
full = true
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -ksp_max_it'
petsc_options_value = 'bcgs bjacobi 1e-10 1e-10 10000 30'
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 0.1
solve_type = NEWTON
[]
[Outputs]
csv = true
execute_on = timestep_end
[]
(modules/porous_flow/test/tests/poro_elasticity/mandel_fully_saturated.i)
# 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
[]
[]
(modules/porous_flow/test/tests/actions/basicthm_hm.i)
# PorousFlowBasicTHM action with coupling_type = HydroMechanical
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
nx = 10
ny = 3
xmax = 10
ymax = 3
[]
[aquifer]
input = gen
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 1 0'
top_right = '10 2 0'
[]
[injection_area]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 'injection_area'
normal = '-1 0 0'
input = 'aquifer'
[]
[outflow_area]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 'outflow_area'
normal = '1 0 0'
input = 'injection_area'
[]
[rename]
type = RenameBlockGenerator
old_block = '0 1'
new_block = 'caprock aquifer'
input = 'outflow_area'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
displacements = 'disp_x disp_y'
biot_coefficient = 1.0
[]
[Variables]
[porepressure]
initial_condition = 1e6
[]
[disp_x]
scaling = 1e-10
[]
[disp_y]
scaling = 1e-10
[]
[]
[AuxVariables]
[temperature]
initial_condition = 293
[]
[]
[PorousFlowBasicTHM]
porepressure = porepressure
temperature = temperature
coupling_type = HydroMechanical
gravity = '0 0 0'
fp = simple_fluid
use_displaced_mesh = false
add_stress_aux = false
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1.5e6
boundary = injection_area
[]
[constant_outflow_porepressure]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = outflow_area
pt_vals = '0 1e9'
multipliers = '0 1e9'
flux_function = 1e-6
PT_shift = 1e6
[]
[top_bottom]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'top bottom'
[]
[right]
type = DirichletBC
variable = disp_x
value = 0
boundary = right
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
solid_bulk_compliance = 2e-7
fluid_bulk_modulus = 1e7
[]
[permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[permeability_caprock]
type = PorousFlowPermeabilityConst
block = caprock
permeability = '1e-15 0 0 0 1e-15 0 0 0 1e-15'
[]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 5e9
poissons_ratio = 0.0
[]
[strain]
type = ComputeSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[]
[Preconditioning]
[basic]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1e4
dt = 1e3
nl_abs_tol = 1e-14
nl_rel_tol = 1e-14
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/test/tests/pressure_pulse/pressure_pulse_1d_fully_saturated_2.i)
# Pressure pulse in 1D with 1 phase - transient
# using the PorousFlowFullySaturatedDarcyBase Kernel
# and the PorousFlowFullySaturatedMassTimeDerivative Kernel
[Mesh]
type = GeneratedMesh
dim = 1
nx = 20
xmin = 0
xmax = 100
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 2E6
[]
[]
[Kernels]
[mass0]
type = PorousFlowFullySaturatedMassTimeDerivative
variable = pp
[]
[flux]
type = PorousFlowFullySaturatedDarcyBase
variable = pp
gravity = '0 0 0'
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
density0 = 1000
thermal_expansion = 0
viscosity = 1e-3
[]
[]
[Materials]
[temperature_qp]
type = PorousFlowTemperature
[]
[ppss_qp]
type = PorousFlow1PhaseFullySaturated
porepressure = pp
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid_qp]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
fluid_bulk_modulus = 2E9
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
[]
[]
[BCs]
[left]
type = DirichletBC
boundary = left
value = 3E6
variable = pp
[]
[]
[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-15 1E-20 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1E3
end_time = 1E4
[]
[Postprocessors]
[p005]
type = PointValue
variable = pp
point = '5 0 0'
execute_on = 'initial timestep_end'
[]
[p015]
type = PointValue
variable = pp
point = '15 0 0'
execute_on = 'initial timestep_end'
[]
[p025]
type = PointValue
variable = pp
point = '25 0 0'
execute_on = 'initial timestep_end'
[]
[p035]
type = PointValue
variable = pp
point = '35 0 0'
execute_on = 'initial timestep_end'
[]
[p045]
type = PointValue
variable = pp
point = '45 0 0'
execute_on = 'initial timestep_end'
[]
[p055]
type = PointValue
variable = pp
point = '55 0 0'
execute_on = 'initial timestep_end'
[]
[p065]
type = PointValue
variable = pp
point = '65 0 0'
execute_on = 'initial timestep_end'
[]
[p075]
type = PointValue
variable = pp
point = '75 0 0'
execute_on = 'initial timestep_end'
[]
[p085]
type = PointValue
variable = pp
point = '85 0 0'
execute_on = 'initial timestep_end'
[]
[p095]
type = PointValue
variable = pp
point = '95 0 0'
execute_on = 'initial timestep_end'
[]
[]
[Outputs]
file_base = pressure_pulse_1d_fully_saturated_2
print_linear_residuals = false
csv = true
[]
(modules/porous_flow/test/tests/actions/basicthm_thm.i)
# PorousFlowBasicTHM action with coupling_type = ThermoHydroMechanical
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
nx = 10
ny = 3
xmax = 10
ymax = 3
[]
[aquifer]
input = gen
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 1 0'
top_right = '10 2 0'
[]
[injection_area]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 'injection_area'
normal = '-1 0 0'
input = 'aquifer'
[]
[outflow_area]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 'outflow_area'
normal = '1 0 0'
input = 'injection_area'
[]
[rename]
type = RenameBlockGenerator
old_block = '0 1'
new_block = 'caprock aquifer'
input = 'outflow_area'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
displacements = 'disp_x disp_y'
biot_coefficient = 1.0
[]
[Variables]
[porepressure]
initial_condition = 1e6
[]
[temperature]
initial_condition = 293
scaling = 1e-6
[]
[disp_x]
scaling = 1e-6
[]
[disp_y]
scaling = 1e-6
[]
[]
[PorousFlowBasicTHM]
porepressure = porepressure
temperature = temperature
coupling_type = ThermoHydroMechanical
gravity = '0 0 0'
fp = simple_fluid
eigenstrain_names = thermal_contribution
use_displaced_mesh = false
add_stress_aux = false
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1.5e6
boundary = injection_area
[]
[constant_injection_temperature]
type = DirichletBC
variable = temperature
value = 313
boundary = injection_area
[]
[constant_outflow_porepressure]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = outflow_area
pt_vals = '0 1e9'
multipliers = '0 1e9'
flux_function = 1e-6
PT_shift = 1e6
[]
[constant_outflow_temperature]
type = DirichletBC
variable = temperature
value = 293
boundary = outflow_area
[]
[top_bottom]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'top bottom'
[]
[right]
type = DirichletBC
variable = disp_x
value = 0
boundary = right
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2e-7
fluid_bulk_modulus = 1e7
[]
[permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[permeability_caprock]
type = PorousFlowPermeabilityConst
block = caprock
permeability = '1e-15 0 0 0 1e-15 0 0 0 1e-15'
[]
[thermal_expansion]
type = PorousFlowConstantThermalExpansionCoefficient
drained_coefficient = 0.003
fluid_coefficient = 0.0002
[]
[rock_internal_energy]
type = PorousFlowMatrixInternalEnergy
density = 2500.0
specific_heat_capacity = 1200.0
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '10 0 0 0 10 0 0 0 10'
block = 'caprock aquifer'
[]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 5e9
poissons_ratio = 0.0
[]
[strain]
type = ComputeSmallStrain
eigenstrain_names = thermal_contribution
[]
[thermal_contribution]
type = ComputeThermalExpansionEigenstrain
temperature = temperature
thermal_expansion_coeff = 0.001
eigenstrain_name = thermal_contribution
stress_free_temperature = 293
[]
[stress]
type = ComputeLinearElasticStress
[]
[]
[Preconditioning]
[basic]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1e4
dt = 1e3
nl_abs_tol = 1e-12
nl_rel_tol = 1E-10
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/test/tests/actions/basicthm_th.i)
# PorousFlowBasicTHM action with coupling_type = ThermoHydroGenerator
# (no mechanical effects)
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
nx = 10
ny = 3
xmax = 10
ymax = 3
[]
[aquifer]
input = gen
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 1 0'
top_right = '10 2 0'
[]
[injection_area]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 'injection_area'
normal = '-1 0 0'
input = 'aquifer'
[]
[outflow_area]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 'outflow_area'
normal = '1 0 0'
input = 'injection_area'
[]
[rename]
type = RenameBlockGenerator
old_block = '0 1'
new_block = 'caprock aquifer'
input = 'outflow_area'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[porepressure]
initial_condition = 1e6
[]
[temperature]
initial_condition = 293
scaling = 1e-6
[]
[]
[PorousFlowBasicTHM]
porepressure = porepressure
temperature = temperature
coupling_type = ThermoHydro
gravity = '0 0 0'
fp = simple_fluid
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1.5e6
boundary = injection_area
[]
[constant_injection_temperature]
type = DirichletBC
variable = temperature
value = 313
boundary = injection_area
[]
[constant_outflow_porepressure]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = outflow_area
pt_vals = '0 1e9'
multipliers = '0 1e9'
flux_function = 1e-6
PT_shift = 1e6
[]
[constant_outflow_temperature]
type = DirichletBC
variable = temperature
value = 293
boundary = outflow_area
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2e-7
fluid_bulk_modulus = 1e7
[]
[permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[permeability_caprock]
type = PorousFlowPermeabilityConst
block = caprock
permeability = '1e-15 0 0 0 1e-15 0 0 0 1e-15'
[]
[thermal_expansion]
type = PorousFlowConstantThermalExpansionCoefficient
biot_coefficient = 0.8
drained_coefficient = 0.003
fluid_coefficient = 0.0002
[]
[rock_internal_energy]
type = PorousFlowMatrixInternalEnergy
density = 2500.0
specific_heat_capacity = 1200.0
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '10 0 0 0 10 0 0 0 10'
block = 'caprock aquifer'
[]
[]
[Preconditioning]
[basic]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1e4
dt = 1e3
nl_abs_tol = 1e-15
nl_rel_tol = 1e-14
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/examples/tutorial/04.i)
# Darcy flow with heat advection and conduction, and elasticity
[Mesh]
[annular]
type = AnnularMeshGenerator
nr = 10
rmin = 1.0
rmax = 10
growth_r = 1.4
nt = 4
dmin = 0
dmax = 90
[]
[make3D]
type = MeshExtruderGenerator
extrusion_vector = '0 0 12'
num_layers = 3
bottom_sideset = 'bottom'
top_sideset = 'top'
input = annular
[]
[shift_down]
type = TransformGenerator
transform = TRANSLATE
vector_value = '0 0 -6'
input = make3D
[]
[aquifer]
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 0 -2'
top_right = '10 10 2'
input = shift_down
[]
[injection_area]
type = ParsedGenerateSideset
combinatorial_geometry = 'x*x+y*y<1.01'
included_subdomains = 1
new_sideset_name = 'injection_area'
input = 'aquifer'
[]
[rename]
type = RenameBlockGenerator
old_block = '0 1'
new_block = 'caps aquifer'
input = 'injection_area'
[]
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
biot_coefficient = 1.0
[]
[Variables]
[porepressure]
[]
[temperature]
initial_condition = 293
scaling = 1E-8
[]
[disp_x]
scaling = 1E-10
[]
[disp_y]
scaling = 1E-10
[]
[disp_z]
scaling = 1E-10
[]
[]
[PorousFlowBasicTHM]
porepressure = porepressure
temperature = temperature
coupling_type = ThermoHydroMechanical
gravity = '0 0 0'
fp = the_simple_fluid
eigenstrain_names = thermal_contribution
use_displaced_mesh = false
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1E6
boundary = injection_area
[]
[constant_injection_temperature]
type = DirichletBC
variable = temperature
value = 313
boundary = injection_area
[]
[roller_tmax]
type = DirichletBC
variable = disp_x
value = 0
boundary = dmax
[]
[roller_tmin]
type = DirichletBC
variable = disp_y
value = 0
boundary = dmin
[]
[roller_top_bottom]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'top bottom'
[]
[cavity_pressure_x]
type = Pressure
boundary = injection_area
variable = disp_x
component = 0
factor = 1E6
use_displaced_mesh = false
[]
[cavity_pressure_y]
type = Pressure
boundary = injection_area
variable = disp_y
component = 1
factor = 1E6
use_displaced_mesh = false
[]
[]
[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
bulk_modulus = 2E9
viscosity = 1.0E-3
density0 = 1000.0
thermal_expansion = 0.0002
cp = 4194
cv = 4186
porepressure_coefficient = 0
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
solid_bulk_compliance = 2E-7
fluid_bulk_modulus = 1E7
[]
[permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1E-14 0 0 0 1E-14 0 0 0 1E-14'
[]
[permeability_caps]
type = PorousFlowPermeabilityConst
block = caps
permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-16'
[]
[thermal_expansion]
type = PorousFlowConstantThermalExpansionCoefficient
drained_coefficient = 0.003
fluid_coefficient = 0.0002
[]
[rock_internal_energy]
type = PorousFlowMatrixInternalEnergy
density = 2500.0
specific_heat_capacity = 1200.0
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '10 0 0 0 10 0 0 0 10'
block = 'caps aquifer'
[]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 5E9
poissons_ratio = 0.0
[]
[strain]
type = ComputeSmallStrain
eigenstrain_names = thermal_contribution
[]
[thermal_contribution]
type = ComputeThermalExpansionEigenstrain
temperature = temperature
thermal_expansion_coeff = 0.001 # this is the linear thermal expansion coefficient
eigenstrain_name = thermal_contribution
stress_free_temperature = 293
[]
[stress]
type = ComputeLinearElasticStress
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_abs_tol = 1E-15
nl_rel_tol = 1E-14
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/test/tests/poro_elasticity/mandel_basicthm.i)
# 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
[]
[]
(modules/porous_flow/test/tests/poro_elasticity/terzaghi_basicthm.i)
# Using a BasicTHM action
# Terzaghi's problem of consolodation of a drained medium
# The FullySaturated Kernels are used, with multiply_by_density = false
# so that this becomes a linear problem with constant Biot Modulus
#
# 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
[]
[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
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.0
bulk_modulus = 8.0
viscosity = 0.96
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 = '2 3'
# bulk modulus is lambda + 2*mu/3 = 2 + 2*3/3 = 4
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[porosity]
type = PorousFlowPorosityConst # only the initial value of this is used
porosity = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.6
fluid_bulk_modulus = 8
solid_bulk_compliance = 0.25
[]
[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'
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_basicthm
[csv]
type = CSV
[]
[]
(modules/porous_flow/test/tests/poro_elasticity/pp_generation_unconfined_fully_saturated_volume.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 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
[]
[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]
[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]
[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.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
component = 2
variable = disp_z
[]
[mass0]
type = PorousFlowFullySaturatedMassTimeDerivative
variable = porepressure
multiply_by_density = false
coupling_type = HydroMechanical
biot_coefficient = 0.3
[]
[source]
type = BodyForce
function = 0.1
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 = 3.3333333333
density0 = 1
thermal_expansion = 0
[]
[]
[Materials]
[temperature_qp]
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
displacements = 'disp_x disp_y disp_z'
[]
[stress]
type = ComputeLinearElasticStress
[]
[eff_fluid_pressure]
type = PorousFlowEffectiveFluidPressure
[]
[vol_strain]
type = PorousFlowVolumetricStrain
[]
[ppss]
type = PorousFlow1PhaseFullySaturated
porepressure = porepressure
[]
[simple_fluid_qp]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[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
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.3
fluid_bulk_modulus = 3.3333333333
solid_bulk_compliance = 0.5
[]
[]
[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
[]
[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_fully_saturated_volume
[csv]
type = CSV
[]
[]
(modules/porous_flow/examples/tidal/barometric_fully_confined.i)
# A fully-confined aquifer is fully saturated with water
# Barometric loading is applied to the aquifer.
# Because the aquifer is assumed to be sandwiched between
# impermeable aquitards, the barometric pressure is not felt
# directly by the porepressure. Instead, the porepressure changes
# only because the barometric loading applies a total stress to
# the top surface of the aquifer.
#
# To replicate standard poroelasticity exactly:
# (1) the PorousFlowBasicTHM Action is used;
# (2) multiply_by_density = false;
# (3) PorousFlowConstantBiotModulus 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'
PorousFlowDictator = dictator
block = 0
biot_coefficient = 0.6
multiply_by_density = false
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[porepressure]
[]
[]
[BCs]
[fix_x]
type = DirichletBC
variable = disp_x
value = 0.0
boundary = 'left right'
[]
[fix_y]
type = DirichletBC
variable = disp_y
value = 0.0
boundary = 'bottom top'
[]
[fix_z_bottom]
type = DirichletBC
variable = disp_z
value = 0.0
boundary = back
[]
[barometric_loading]
type = FunctionNeumannBC
variable = disp_z
function = -1000.0 # atmospheric pressure increase of 1kPa
boundary = front
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
[]
[]
[PorousFlowBasicTHM]
coupling_type = HydroMechanical
displacements = 'disp_x disp_y disp_z'
porepressure = porepressure
gravity = '0 0 0'
fp = the_simple_fluid
[]
[Materials]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
bulk_modulus = 10.0E9 # drained bulk modulus
poissons_ratio = 0.25
[]
[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
solid_bulk_compliance = 1E-10
fluid_bulk_modulus = 2E9
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[]
[Postprocessors]
[pp]
type = PointValue
point = '0.5 0.5 0.5'
variable = porepressure
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 1
[]
[Outputs]
console = true
csv = true
[]
(modules/porous_flow/test/tests/dirackernels/theis_rz.i)
# Theis problem: Flow to single sink using BasicTHM
# SinglePhase
# RZ mesh
[Mesh]
type = GeneratedMesh
dim = 1
nx = 20
xmax = 100
bias_x = 1.05
coord_type = RZ
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
initial_condition = 20E6
[]
[]
[PorousFlowBasicTHM]
dictator_name = dictator
add_darcy_aux = false
fp = simple_fluid
gravity = '0 0 0'
multiply_by_density = false
porepressure = pp
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
viscosity = 0.001
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.05
[]
[biot_mod]
type = PorousFlowConstantBiotModulus
fluid_bulk_modulus = 2E9
biot_coefficient = 1.0
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-14 0 0 0 1E-14 0 0 0 1E-14'
[]
[]
[DiracKernels]
[sink]
type = PorousFlowSquarePulsePointSource
point = '0 0 0'
mass_flux = -0.16E-3 # recall this is a volumetric flux because multiply_by_density = false in the Action, so this corresponds to a mass_flux of 0.16 kg/s/m because density=1000
variable = pp
[]
[]
[VectorPostprocessors]
[pp]
type = LineValueSampler
num_points = 25
start_point = '0 0 0'
end_point = '100 0 0'
sort_by = x
variable = pp
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 200
end_time = 1E3
nl_abs_tol = 1e-10
[]
[Outputs]
perf_graph = true
[csv]
type = CSV
execute_on = final
[]
[]
(modules/porous_flow/examples/tidal/earth_tide_fullsat.i)
# A confined aquifer is fully saturated with water
# Earth tides apply strain to the aquifer and the resulting porepressure changes are recorded
#
# To replicate standard poroelasticity exactly:
# (1) the PorousFlowBasicTHM Action is used;
# (2) multiply_by_density = false;
# (3) PorousFlowConstantBiotModulus 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'
PorousFlowDictator = dictator
block = 0
biot_coefficient = 0.6
multiply_by_density = false
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[porepressure]
[]
[]
[BCs]
[strain_x]
type = FunctionDirichletBC
variable = disp_x
function = earth_tide_x
boundary = 'left right'
[]
[strain_y]
type = FunctionDirichletBC
variable = disp_y
function = earth_tide_y
boundary = 'bottom top'
[]
[strain_z]
type = FunctionDirichletBC
variable = disp_z
function = earth_tide_z
boundary = 'back front'
[]
[]
[Functions]
[earth_tide_x]
type = ParsedFunction
expression = 'x*1E-8*(5*cos(t*2*pi) + 2*cos((t-0.5)*2*pi) + 1*cos((t+0.3)*0.5*pi))'
[]
[earth_tide_y]
type = ParsedFunction
expression = 'y*1E-8*(7*cos(t*2*pi) + 4*cos((t-0.3)*2*pi) + 7*cos((t+0.6)*0.5*pi))'
[]
[earth_tide_z]
type = ParsedFunction
expression = 'z*1E-8*(7*cos((t-0.5)*2*pi) + 4*cos((t-0.8)*2*pi) + 7*cos((t+0.1)*4*pi))'
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
[]
[]
[PorousFlowBasicTHM]
coupling_type = HydroMechanical
displacements = 'disp_x disp_y disp_z'
porepressure = porepressure
gravity = '0 0 0'
fp = the_simple_fluid
[]
[Materials]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
bulk_modulus = 10.0E9 # drained bulk modulus
poissons_ratio = 0.25
[]
[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
solid_bulk_compliance = 1E-10
fluid_bulk_modulus = 2E9
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
[]
[]
[Postprocessors]
[pp]
type = PointValue
point = '0.5 0.5 0.5'
variable = porepressure
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 0.01
end_time = 2
[]
[Outputs]
console = true
csv = true
[]
(modules/porous_flow/test/tests/hysteresis/except07.i)
# Exception testing of PorousFlowHysteresisOrder
# Incorrectly ordered previous_turning_points
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 1
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[]
[PorousFlowBasicTHM]
porepressure = pp
fp = simple_fluid
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2e-7
fluid_bulk_modulus = 1e7
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[hys_order]
type = PorousFlowHysteresisOrder
initial_order = 3
previous_turning_points = '0.6 0.8 0.9'
[]
[]
[Preconditioning]
[basic]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 1
[]
(modules/porous_flow/test/tests/sinks/PorousFlowPiecewiseLinearSink_BC_eg1.i)
## This is an example input file showing how to set a Type I (Dirichlet) BC with PorousFlowPiecewiseLinearSink
##
## Problem setup:
## - The boundaries are set to P(x = 0) = 2e6 Pa, P(x = 1) = 1e6 and run to steady state.
## - The 2d domain is 1 m x 1 m
## - The permeability is set to 1E-15 m2, fluid viscosity = 1E-3 Pa-s
## - The steady state flux is calculated q = -k/mu*grad(P) = 1e-6 m/s
##
## Problem verification (in csv output):
## - The flux in and out of the domain are 1e-6 m/s (matching steady state solution)
## - The pressure at the left and right boundaries are set to 2e6 and 1e6 Pa, respectively
[Mesh]
type = GeneratedMesh
dim = 2
nx = 5
xmin = 0
xmax = 1
ny = 2
ymin = 0
ymax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[porepressure]
initial_condition = 1.5e6 # initial pressure in domain
[]
[]
[PorousFlowBasicTHM]
porepressure = porepressure
coupling_type = Hydro
gravity = '0 0 0'
fp = the_simple_fluid
[]
[AuxVariables]
[fluxes_out]
[]
[fluxes_in]
[]
[]
[BCs]
[in_left]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = 'left'
pt_vals = '-1e9 1e9' # x coordinates defining g
multipliers = '-1e9 1e9' # y coordinates defining g
PT_shift = 2.E6 # BC pressure
flux_function = 1E-5 # Variable C
fluid_phase = 0
save_in = fluxes_out
[]
[out_right]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
boundary = 'right'
pt_vals = '-1e9 1e9' # x coordinates defining g
multipliers = '-1e9 1e9' # y coordinates defining g
PT_shift = 1.E6 # BC pressure
flux_function = 1E-6 # Variable C
fluid_phase = 0
save_in = fluxes_in
[]
[]
[Postprocessors]
[left_flux]
type = NodalSum
boundary = 'left'
variable = fluxes_out
execute_on = 'timestep_end'
[]
[right_flux]
type = NodalSum
boundary = 'right'
variable = fluxes_in
execute_on = 'timestep_end'
[]
[left_pressure]
type = SideAverageValue
boundary = 'left'
variable = porepressure
execute_on = 'timestep_end'
[]
[right_pressure]
type = SideAverageValue
boundary = 'right'
variable = porepressure
execute_on = 'timestep_end'
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0
viscosity = 1.0E-3
density0 = 1000.0
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2E-7
fluid_bulk_modulus = 1E7
[]
[permeability_aquifer]
type = PorousFlowPermeabilityConst
permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_abs_tol = 1E-10
[]
[Outputs]
csv = true
[]
(modules/porous_flow/test/tests/poro_elasticity/terzaghi_fully_saturated_volume.i)
# Terzaghi's problem of consolodation of a drained medium
# The FullySaturated Kernels are used, with multiply_by_density = false
# so that this becomes a linear problem with constant Biot Modulus
#
# 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
[]
[]
[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
[]
[mass0]
type = PorousFlowFullySaturatedMassTimeDerivative
coupling_type = HydroMechanical
biot_coefficient = 0.6
multiply_by_density = false
variable = porepressure
[]
[flux]
type = PorousFlowFullySaturatedDarcyBase
multiply_by_density = false
variable = porepressure
gravity = '0 0 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_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 used
porosity = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.6
fluid_bulk_modulus = 8
solid_bulk_compliance = 0.25
[]
[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'
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_fully_saturated_volume
[csv]
type = CSV
[]
[]
(modules/porous_flow/examples/tutorial/01.i)
# Darcy flow
[Mesh]
[annular]
type = AnnularMeshGenerator
nr = 10
rmin = 1.0
rmax = 10
growth_r = 1.4
nt = 4
dmin = 0
dmax = 90
[]
[make3D]
type = MeshExtruderGenerator
extrusion_vector = '0 0 12'
num_layers = 3
bottom_sideset = 'bottom'
top_sideset = 'top'
input = annular
[]
[shift_down]
type = TransformGenerator
transform = TRANSLATE
vector_value = '0 0 -6'
input = make3D
[]
[aquifer]
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 0 -2'
top_right = '10 10 2'
input = shift_down
[]
[injection_area]
type = ParsedGenerateSideset
combinatorial_geometry = 'x*x+y*y<1.01'
included_subdomains = 1
new_sideset_name = 'injection_area'
input = 'aquifer'
[]
[rename]
type = RenameBlockGenerator
old_block = '0 1'
new_block = 'caps aquifer'
input = 'injection_area'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[porepressure]
[]
[]
[PorousFlowBasicTHM]
porepressure = porepressure
coupling_type = Hydro
gravity = '0 0 0'
fp = the_simple_fluid
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1E6
boundary = injection_area
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E9
viscosity = 1.0E-3
density0 = 1000.0
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2E-7
fluid_bulk_modulus = 1E7
[]
[permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1E-14 0 0 0 1E-14 0 0 0 1E-14'
[]
[permeability_caps]
type = PorousFlowPermeabilityConst
block = caps
permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-16'
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_abs_tol = 1E-13
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/test/tests/hysteresis/except09.i)
# Exception testing of PorousFlowPropertyAux
# hystresis_turning_point too large
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 1
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[]
[PorousFlowBasicTHM]
porepressure = pp
fp = simple_fluid
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2e-7
fluid_bulk_modulus = 1e7
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[hys_order]
type = PorousFlowHysteresisOrder
[]
[]
[AuxVariables]
[tp]
family = MONOMIAL
order = CONSTANT
[]
[]
[AuxKernels]
[tp]
type = PorousFlowPropertyAux
variable = tp
property = hysteresis_saturation_turning_point
hysteresis_turning_point = 3
[]
[]
[Preconditioning]
[basic]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 1
[]
(modules/porous_flow/examples/tidal/atm_tides_open_hole.i)
# A 100m x 10m "slab" of height 100m is subjected to cyclic pressure at its top
# Assumptions:
# the boundaries are impermeable, except the top boundary
# only vertical displacement is allowed
# the atmospheric pressure sets the total stress at the top of the model
# at the slab left-hand side there is a borehole that taps into the base of the slab.
[Mesh]
[the_mesh]
type = GeneratedMeshGenerator
dim = 3
nx = 10
ny = 1
nz = 10
xmin = 0
xmax = 100
ymin = -5
ymax = 5
zmin = -100
zmax = 0
[]
[bh_back]
type = ExtraNodesetGenerator
coord = '0 -5 -100'
input = the_mesh
new_boundary = 11
[]
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
biot_coefficient = 0.6
multiply_by_density = false
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[porepressure]
scaling = 1E11
[]
[]
[ICs]
[porepressure]
type = FunctionIC
variable = porepressure
function = '-10000*z' # this is only approximately correct
[]
[]
[Functions]
[ini_stress_zz]
type = ParsedFunction
expression = '(25000 - 0.6*10000)*z' # remember this is effective stress
[]
[cyclic_porepressure]
type = ParsedFunction
expression = 'if(t>0,5000 * sin(2 * pi * t / 3600.0 / 24.0),0)'
[]
[cyclic_porepressure_at_depth]
type = ParsedFunction
expression = '-10000*z + if(t>0,5000 * sin(2 * pi * t / 3600.0 / 24.0),0)'
[]
[neg_cyclic_porepressure]
type = ParsedFunction
expression = '-if(t>0,5000 * sin(2 * pi * t / 3600.0 / 24.0),0)'
[]
[]
[BCs]
# zmin is called 'back'
# zmax is called 'front'
# ymin is called 'bottom'
# ymax is called 'top'
# xmin is called 'left'
# xmax is called 'right'
[no_x_disp]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'bottom top' # because of 1-element meshing, this fixes u_x=0 everywhere
[]
[no_y_disp]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom top' # because of 1-element meshing, this fixes u_y=0 everywhere
[]
[no_z_disp_at_bottom]
type = DirichletBC
variable = disp_z
value = 0
boundary = back
[]
[pp]
type = FunctionDirichletBC
variable = porepressure
function = cyclic_porepressure
boundary = front
[]
[pp_downhole]
type = FunctionDirichletBC
variable = porepressure
function = cyclic_porepressure_at_depth
boundary = 11
[]
[total_stress_at_top]
type = FunctionNeumannBC
variable = disp_z
function = neg_cyclic_porepressure
boundary = front
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.0
bulk_modulus = 2E9
viscosity = 1E-3
density0 = 1000.0
[]
[]
[PorousFlowBasicTHM]
coupling_type = HydroMechanical
displacements = 'disp_x disp_y disp_z'
porepressure = porepressure
gravity = '0 0 -10'
fp = the_simple_fluid
[]
[Materials]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
bulk_modulus = 10.0E9 # drained bulk modulus
poissons_ratio = 0.25
[]
[strain]
type = ComputeSmallStrain
eigenstrain_names = ini_stress
[]
[stress]
type = ComputeLinearElasticStress
[]
[ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '0 0 0 0 0 0 0 0 ini_stress_zz'
eigenstrain_name = ini_stress
[]
[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 = 2E9
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-14 0 0 0 1E-14 0 0 0 1E-14'
[]
[density]
type = GenericConstantMaterial
prop_names = density
prop_values = 2500.0
[]
[]
[Postprocessors]
[p0_0]
type = PointValue
outputs = csv
point = '0 0 0'
variable = porepressure
[]
[p100_0]
type = PointValue
outputs = csv
point = '100 0 0'
variable = porepressure
[]
[p0_100]
type = PointValue
outputs = csv
point = '0 0 -100'
variable = porepressure
[]
[p100_100]
type = PointValue
outputs = csv
point = '100 0 -100'
variable = porepressure
[]
[uz0]
type = PointValue
outputs = csv
point = '0 0 0'
variable = disp_z
[]
[uz100]
type = PointValue
outputs = csv
point = '100 0 0'
variable = disp_z
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = -3600
dt = 3600
end_time = 172800
nl_rel_tol = 1E-10
nl_abs_tol = 1E-5
[]
[Outputs]
print_linear_residuals = false
csv = true
[]
(modules/porous_flow/test/tests/thm_rehbinder/fixed_outer.i)
[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
[]
# To get consistent ordering of results with distributed meshes
allow_renumbering = false
[]
[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]
[smp]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_rtol'
petsc_options_value = 'gmres asm lu 1E-10'
[]
[]
[Executioner]
type = Steady
solve_type = Newton
[]
[Outputs]
file_base = fixed_outer
execute_on = timestep_end
csv = true
[]
(modules/porous_flow/test/tests/poro_elasticity/pp_generation_unconfined_fully_saturated.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 has units 1/time. 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/s/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
[]
[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]
[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]
[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.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
component = 2
variable = disp_z
[]
[mass0]
type = PorousFlowFullySaturatedMassTimeDerivative
variable = porepressure
coupling_type = HydroMechanical
biot_coefficient = 0.3
[]
[source]
type = BodyForce
function = '0.1*exp(8.163265306*0.1*t/3.3333333333)'
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 = 3.3333333333
density0 = 1
thermal_expansion = 0
[]
[]
[Materials]
[temperature_qp]
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
displacements = 'disp_x disp_y disp_z'
[]
[stress]
type = ComputeLinearElasticStress
[]
[eff_fluid_pressure]
type = PorousFlowEffectiveFluidPressure
[]
[vol_strain]
type = PorousFlowVolumetricStrain
[]
[ppss]
type = PorousFlow1PhaseFullySaturated
porepressure = porepressure
[]
[simple_fluid_qp]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[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
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.3
fluid_bulk_modulus = 3.3333333333
solid_bulk_compliance = 0.5
[]
[]
[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
[]
[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_fully_saturated
[csv]
type = CSV
[]
[]
(modules/porous_flow/test/tests/hysteresis/except05.i)
# Exception testing of PorousFlowHysteresisOrder
# Incorrect: previous_turning_points not in the range [0, 1]
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 1
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[]
[PorousFlowBasicTHM]
porepressure = pp
fp = simple_fluid
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2e-7
fluid_bulk_modulus = 1e7
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[hys_order]
type = PorousFlowHysteresisOrder
initial_order = 1
previous_turning_points = 1.1
[]
[]
[Preconditioning]
[basic]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 1
[]
(modules/porous_flow/examples/tutorial/05.i)
# Darcy flow with heat advection and conduction, using Water97 properties
[Mesh]
[annular]
type = AnnularMeshGenerator
nr = 10
rmin = 1.0
rmax = 10
growth_r = 1.4
nt = 4
dmin = 0
dmax = 90
[]
[make3D]
type = MeshExtruderGenerator
extrusion_vector = '0 0 12'
num_layers = 3
bottom_sideset = 'bottom'
top_sideset = 'top'
input = annular
[]
[shift_down]
type = TransformGenerator
transform = TRANSLATE
vector_value = '0 0 -6'
input = make3D
[]
[aquifer]
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 0 -2'
top_right = '10 10 2'
input = shift_down
[]
[injection_area]
type = ParsedGenerateSideset
combinatorial_geometry = 'x*x+y*y<1.01'
included_subdomains = 1
new_sideset_name = 'injection_area'
input = 'aquifer'
[]
[rename]
type = RenameBlockGenerator
old_block = '0 1'
new_block = 'caps aquifer'
input = 'injection_area'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[porepressure]
initial_condition = 1E6
[]
[temperature]
initial_condition = 313
scaling = 1E-8
[]
[]
[PorousFlowBasicTHM]
porepressure = porepressure
temperature = temperature
coupling_type = ThermoHydro
gravity = '0 0 0'
fp = the_simple_fluid
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 2E6
boundary = injection_area
[]
[constant_injection_temperature]
type = DirichletBC
variable = temperature
value = 333
boundary = injection_area
[]
[]
[FluidProperties]
[the_simple_fluid]
type = Water97FluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2E-7
fluid_bulk_modulus = 1E7
[]
[permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1E-14 0 0 0 1E-14 0 0 0 1E-14'
[]
[permeability_caps]
type = PorousFlowPermeabilityConst
block = caps
permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-16'
[]
[thermal_expansion]
type = PorousFlowConstantThermalExpansionCoefficient
biot_coefficient = 0.8
drained_coefficient = 0.003
fluid_coefficient = 0.0002
[]
[rock_internal_energy]
type = PorousFlowMatrixInternalEnergy
density = 2500.0
specific_heat_capacity = 1200.0
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '10 0 0 0 10 0 0 0 10'
block = 'caps aquifer'
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_abs_tol = 1E-10
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/test/tests/hysteresis/except01.i)
# Exception testing of PorousFlowHysteresisOrder
# Incorrect: liquid_phase = 1
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 1
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[]
[PorousFlowBasicTHM]
porepressure = pp
fp = simple_fluid
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2e-7
fluid_bulk_modulus = 1e7
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[hys_order]
type = PorousFlowHysteresisOrder
liquid_phase = 1
[]
[]
[Preconditioning]
[basic]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 1
[]
(modules/porous_flow/test/tests/poro_elasticity/mandel_fully_saturated_volume.i)
# 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
[]
[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
multiply_by_density = false
coupling_type = HydroMechanical
variable = porepressure
[]
[flux]
type = PorousFlowFullySaturatedDarcyBase
multiply_by_density = false
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_volume
[csv]
time_step_interval = 3
type = CSV
[]
[]
(modules/porous_flow/test/tests/hysteresis/except06.i)
# Exception testing of PorousFlowHysteresisOrder
# Incorrectly ordered previous_turning_points
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 1
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[]
[PorousFlowBasicTHM]
porepressure = pp
fp = simple_fluid
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2e-7
fluid_bulk_modulus = 1e7
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
[]
[hys_order]
type = PorousFlowHysteresisOrder
initial_order = 2
previous_turning_points = '0.6 0.4'
[]
[]
[Preconditioning]
[basic]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 1
[]
(modules/porous_flow/examples/tutorial/05_tabulated.i)
# Darcy flow with heat advection and conduction, using Water97 properties
[Mesh]
[annular]
type = AnnularMeshGenerator
nr = 10
rmin = 1.0
rmax = 10
growth_r = 1.4
nt = 4
dmin = 0
dmax = 90
[]
[make3D]
type = MeshExtruderGenerator
extrusion_vector = '0 0 12'
num_layers = 3
bottom_sideset = 'bottom'
top_sideset = 'top'
input = annular
[]
[shift_down]
type = TransformGenerator
transform = TRANSLATE
vector_value = '0 0 -6'
input = make3D
[]
[aquifer]
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 0 -2'
top_right = '10 10 2'
input = shift_down
[]
[injection_area]
type = ParsedGenerateSideset
combinatorial_geometry = 'x*x+y*y<1.01'
included_subdomains = 1
new_sideset_name = 'injection_area'
input = 'aquifer'
[]
[rename]
type = RenameBlockGenerator
old_block = '0 1'
new_block = 'caps aquifer'
input = 'injection_area'
[]
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[porepressure]
initial_condition = 1E6
[]
[temperature]
initial_condition = 313
scaling = 1E-8
[]
[]
[PorousFlowBasicTHM]
porepressure = porepressure
temperature = temperature
coupling_type = ThermoHydro
gravity = '0 0 0'
fp = tabulated_water
[]
[BCs]
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 2E6
boundary = injection_area
[]
[constant_injection_temperature]
type = DirichletBC
variable = temperature
value = 333
boundary = injection_area
[]
[]
[FluidProperties]
[true_water]
type = Water97FluidProperties
[]
[tabulated_water]
type = TabulatedFluidProperties
fp = true_water
temperature_min = 275
interpolated_properties = 'density viscosity enthalpy internal_energy'
fluid_property_output_file = water97_tabulated.csv
# Comment out the fp parameter and uncomment below to use the newly generated tabulation
# fluid_property_file = water97_tabulated.csv
[]
[]
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2E-7
fluid_bulk_modulus = 1E7
[]
[permeability_aquifer]
type = PorousFlowPermeabilityConst
block = aquifer
permeability = '1E-14 0 0 0 1E-14 0 0 0 1E-14'
[]
[permeability_caps]
type = PorousFlowPermeabilityConst
block = caps
permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-16'
[]
[thermal_expansion]
type = PorousFlowConstantThermalExpansionCoefficient
biot_coefficient = 0.8
drained_coefficient = 0.003
fluid_coefficient = 0.0002
[]
[rock_internal_energy]
type = PorousFlowMatrixInternalEnergy
density = 2500.0
specific_heat_capacity = 1200.0
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '10 0 0 0 10 0 0 0 10'
block = 'caps aquifer'
[]
[]
[Preconditioning]
active = basic
[basic]
type = SMP
full = true
petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP
full = true
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_abs_tol = 1E-10
[]
[Outputs]
exodus = true
[]
(modules/porous_flow/test/tests/jacobian/mass01_fully_saturated.i)
# FullySaturatedMassTimeDerivative
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
displacements = 'disp_x disp_y disp_z'
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.5
bulk_modulus = 1.5
density0 = 1.0
[]
[]
[Variables]
[pp]
[]
[T]
[]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[]
[ICs]
[disp_x]
type = RandomIC
min = -0.1
max = 0.1
variable = disp_x
[]
[disp_y]
type = RandomIC
min = -0.1
max = 0.1
variable = disp_y
[]
[disp_z]
type = RandomIC
min = -0.1
max = 0.1
variable = disp_z
[]
[pp]
type = RandomIC
variable = pp
min = 0
max = 1
[]
[T]
type = RandomIC
variable = T
min = 0
max = 1
[]
[]
[BCs]
# necessary otherwise volumetric strain rate will be zero
[disp_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[disp_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'left right'
[]
[disp_z]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'left right'
[]
[]
[Kernels]
[mass0]
type = PorousFlowFullySaturatedMassTimeDerivative
variable = pp
coupling_type = ThermoHydroMechanical
biot_coefficient = 0.9
[]
[dummyT]
type = TimeDerivative
variable = T
[]
[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
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp disp_x disp_y disp_z T'
number_fluid_phases = 1
number_fluid_components = 1
[]
[simple1]
type = TensorMechanicsPlasticSimpleTester
a = 0
b = 1
strength = 1E20
yield_function_tolerance = 1.0E-9
internal_constraint_tolerance = 1.0E-9
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
bulk_modulus = 2.0
shear_modulus = 3.0
[]
[strain]
type = ComputeSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[vol_strain]
type = PorousFlowVolumetricStrain
[]
[temperature]
type = PorousFlowTemperature
temperature = T
[]
[ppss]
type = PorousFlow1PhaseFullySaturated
porepressure = pp
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = the_simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst # only the initial vaue of this is ever used
porosity = 0.1
[]
[biot_modulus]
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.9
fluid_bulk_modulus = 1.5
solid_bulk_compliance = 0.5
[]
[thermal_expansion]
type = PorousFlowConstantThermalExpansionCoefficient
biot_coefficient = 0.9
fluid_coefficient = 0.5
drained_coefficient = 0.4
[]
[]
[Preconditioning]
[check]
type = SMP
full = true
#petsc_options = '-snes_test_display'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 2
[]
[Outputs]
exodus = false
[]
(modules/porous_flow/test/tests/thm_rehbinder/free_outer.i)
[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
[]
(modules/porous_flow/examples/tidal/atm_tides.i)
# A 10m x 10m "column" of height 100m is subjected to cyclic pressure at its top
# Assumptions:
# the boundaries are impermeable, except the top boundary
# only vertical displacement is allowed
# the atmospheric pressure sets the total stress at the top of the model
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 10
xmin = 0
xmax = 10
ymin = 0
ymax = 10
zmin = -100
zmax = 0
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
biot_coefficient = 0.6
multiply_by_density = false
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[porepressure]
scaling = 1E11
[]
[]
[ICs]
[porepressure]
type = FunctionIC
variable = porepressure
function = '-10000*z' # approximately correct
[]
[]
[Functions]
[ini_stress_zz]
type = ParsedFunction
expression = '(25000 - 0.6*10000)*z' # remember this is effective stress
[]
[cyclic_porepressure]
type = ParsedFunction
expression = 'if(t>0,5000 * sin(2 * pi * t / 3600.0 / 24.0),0)'
[]
[neg_cyclic_porepressure]
type = ParsedFunction
expression = '-if(t>0,5000 * sin(2 * pi * t / 3600.0 / 24.0),0)'
[]
[]
[BCs]
# zmin is called 'back'
# zmax is called 'front'
# ymin is called 'bottom'
# ymax is called 'top'
# xmin is called 'left'
# xmax is called 'right'
[no_x_disp]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'bottom top' # because of 1-element meshing, this fixes u_x=0 everywhere
[]
[no_y_disp]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom top' # because of 1-element meshing, this fixes u_y=0 everywhere
[]
[no_z_disp_at_bottom]
type = DirichletBC
variable = disp_z
value = 0
boundary = back
[]
[pp]
type = FunctionDirichletBC
variable = porepressure
function = cyclic_porepressure
boundary = front
[]
[total_stress_at_top]
type = FunctionNeumannBC
variable = disp_z
function = neg_cyclic_porepressure
boundary = front
[]
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.0
bulk_modulus = 2E9
viscosity = 1E-3
density0 = 1000.0
[]
[]
[PorousFlowBasicTHM]
coupling_type = HydroMechanical
displacements = 'disp_x disp_y disp_z'
porepressure = porepressure
gravity = '0 0 -10'
fp = the_simple_fluid
[]
[Materials]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
bulk_modulus = 10.0E9 # drained bulk modulus
poissons_ratio = 0.25
[]
[strain]
type = ComputeSmallStrain
eigenstrain_names = ini_stress
[]
[stress]
type = ComputeLinearElasticStress
[]
[ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '0 0 0 0 0 0 0 0 ini_stress_zz'
eigenstrain_name = ini_stress
[]
[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 = 2E9
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1E-12 0 0 0 1E-12 0 0 0 1E-14'
[]
[density]
type = GenericConstantMaterial
prop_names = density
prop_values = 2500.0
[]
[]
[Postprocessors]
[p0]
type = PointValue
outputs = csv
point = '0 0 0'
variable = porepressure
[]
[uz0]
type = PointValue
outputs = csv
point = '0 0 0'
variable = disp_z
[]
[p100]
type = PointValue
outputs = csv
point = '0 0 -100'
variable = porepressure
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
start_time = -3600 # so postprocessors get recorded correctly at t=0
dt = 3600
end_time = 360000
nl_abs_tol = 5E-7
nl_rel_tol = 1E-10
[]
[Outputs]
csv = true
[]
(modules/porous_flow/test/tests/poro_elasticity/pp_generation_unconfined_basicthm.i)
# 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
[]
[]