- variableThe name of the variable that this postprocessor operates on
C++ Type:std::vector
Description:The name of the variable that this postprocessor operates on
NodalExtremeValue
under construction:Undocumented Class
The NodalExtremeValue has not been documented. The content listed below should be used as a starting point for documenting the class, which includes the typical automatic documentation associated with a MooseObject; however, what is contained is ultimately determined by what is necessary to make the documentation clear for users.
# NodalExtremeValue
!syntax description /Postprocessors/NodalExtremeValue
## Overview
!! Replace these lines with information regarding the NodalExtremeValue object.
## Example Input File Syntax
!! Describe and include an example of how to use the NodalExtremeValue object.
!syntax parameters /Postprocessors/NodalExtremeValue
!syntax inputs /Postprocessors/NodalExtremeValue
!syntax children /Postprocessors/NodalExtremeValue
!syntax description /Postprocessors/NodalExtremeValue
Input Parameters
- blockThe list of block ids (SubdomainID) that this object will be applied
C++ Type:std::vector
Options:
Description:The list of block ids (SubdomainID) that this object will be applied
- boundaryThe list of boundary IDs from the mesh where this boundary condition applies
C++ Type:std::vector
Options:
Description:The list of boundary IDs from the mesh where this boundary condition applies
- execute_onTIMESTEP_ENDThe list of flag(s) indicating when this object should be executed, the available options include NONE, INITIAL, LINEAR, NONLINEAR, TIMESTEP_END, TIMESTEP_BEGIN, FINAL, CUSTOM.
Default:TIMESTEP_END
C++ Type:ExecFlagEnum
Options:NONE INITIAL LINEAR NONLINEAR TIMESTEP_END TIMESTEP_BEGIN FINAL CUSTOM TRANSFER
Description:The list of flag(s) indicating when this object should be executed, the available options include NONE, INITIAL, LINEAR, NONLINEAR, TIMESTEP_END, TIMESTEP_BEGIN, FINAL, CUSTOM.
- unique_node_executeFalseWhen false (default), block restricted objects will have the execute method called multiple times on a single node if the node lies on a interface between two subdomains.
Default:False
C++ Type:bool
Options:
Description:When false (default), block restricted objects will have the execute method called multiple times on a single node if the node lies on a interface between two subdomains.
- value_typemaxType of extreme value to return. 'max' returns the maximum value. 'min' returns the minimum value.
Default:max
C++ Type:MooseEnum
Options:max min
Description:Type of extreme value to return. 'max' returns the maximum value. 'min' returns the minimum value.
Optional Parameters
- allow_duplicate_execution_on_initialFalseIn the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
Default:False
C++ Type:bool
Options:
Description:In the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector
Options:
Description:Adds user-defined labels for accessing object parameters via control logic.
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Options:
Description:Set the enabled status of the MooseObject.
- outputsVector of output names were you would like to restrict the output of variables(s) associated with this object
C++ Type:std::vector
Options:
Description:Vector of output names were you would like to restrict the output of variables(s) associated with this object
- seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Options:
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
Options:
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
Input Files
- modules/tensor_mechanics/test/tests/plane_stress/weak_plane_stress_finite.i
- modules/tensor_mechanics/examples/coal_mining/fine.i
- modules/solid_mechanics/test/tests/generalized_plane_strain/generalized_plane_strain_field.i
- modules/porous_flow/examples/coal_mining/coarse_with_fluid.i
- modules/porous_flow/examples/coal_mining/fine_with_fluid.i
- modules/solid_mechanics/test/tests/generalized_plane_strain/weak_planestress.i
- test/tests/postprocessors/nodal_extreme_value/nodal_extreme_pps_test.i
- modules/solid_mechanics/test/tests/generalized_plane_strain/generalized_plane_strain_field_kinematic.i
- modules/tensor_mechanics/test/tests/plane_stress/weak_plane_stress_small.i
- modules/tensor_mechanics/test/tests/plane_stress/weak_plane_stress_incremental.i
- modules/tensor_mechanics/examples/coal_mining/coarse.i
- modules/solid_mechanics/test/tests/generalized_plane_strain/out_of_plane_pressure.i
modules/tensor_mechanics/test/tests/plane_stress/weak_plane_stress_finite.i
[GlobalParams]
order = FIRST
family = LAGRANGE
displacements = 'disp_x disp_y'
temperature = temp
out_of_plane_strain = strain_zz
[]
[Mesh]
file = square.e
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./strain_zz]
[../]
[]
[AuxVariables]
[./temp]
[../]
[./nl_strain_zz]
order = CONSTANT
family = MONOMIAL
[../]
[]
[Postprocessors]
[./react_z]
type = MaterialTensorIntegral
rank_two_tensor = stress
index_i = 2
index_j = 2
[../]
[./min_strain_zz]
type = NodalExtremeValue
variable = strain_zz
value_type = min
[../]
[./max_strain_zz]
type = NodalExtremeValue
variable = strain_zz
value_type = max
[../]
[]
[Modules/TensorMechanics/Master]
[plane_stress]
planar_formulation = WEAK_PLANE_STRESS
strain = FINITE
generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_xx strain_xy strain_yy'
eigenstrain_names = eigenstrain
[]
[]
[AuxKernels]
[./tempfuncaux]
type = FunctionAux
variable = temp
function = tempfunc
use_displaced_mesh = false
[../]
[./strain_zz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = nl_strain_zz
index_i = 2
index_j = 2
[../]
[]
[Functions]
[./pull]
type = PiecewiseLinear
x='0 1 100'
y='0 0.00 0.00'
[../]
[./tempfunc]
type = ParsedFunction
value = '(1 - x) * t'
[../]
[]
[BCs]
[./bottomx]
type = DirichletBC
boundary = 1
variable = disp_x
value = 0.0
[../]
[./bottomy]
type = DirichletBC
boundary = 1
variable = disp_y
value = 0.0
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
poissons_ratio = 0.3
youngs_modulus = 1e6
[../]
[./thermal_strain]
type = ComputeThermalExpansionEigenstrain
thermal_expansion_coeff = 0.02
stress_free_temperature = 0.5
eigenstrain_name = eigenstrain
[../]
[./stress]
type = ComputeFiniteStrainElasticStress
[../]
[]
[Executioner]
type = Transient
solve_type = PJFNK
line_search = none
# controls for linear iterations
l_max_its = 100
l_tol = 1e-06
# controls for nonlinear iterations
nl_max_its = 15
nl_rel_tol = 1e-12
nl_abs_tol = 1e-10
# time control
start_time = 0.0
dt = 1.0
dtmin = 1.0
end_time = 2.0
[]
[Outputs]
exodus = true
[]
modules/tensor_mechanics/examples/coal_mining/fine.i
# Strata deformation and fracturing around a coal mine - 3D model
#
# A "half model" is used. The mine is 400m deep and
# just the roof is studied (-400<=z<=0). The mining panel
# sits between 0<=x<=150, and 0<=y<=1000, so this simulates
# a coal panel that is 300m wide and 1000m long. The outer boundaries
# are 1km from the excavation boundaries.
#
# Time is meaningless in this example
# as quasi-static solutions are sought at each timestep, but
# the number of timesteps controls the resolution of the
# process.
#
# The boundary conditions for this simulation are:
# - disp_x = 0 at x=0 and x=1150
# - disp_y = 0 at y=-1000 and y=1000
# - disp_z = 0 at z=-400, but there is a time-dependent
# Young's modulus that simulates excavation
# - wc_x = 0 at y=-1000 and y=1000
# - wc_y = 0 at x=0 and x=1150
# That is, rollers on the sides, free at top,
# and prescribed at bottom in the unexcavated portion.
#
# The small strain formulation is used.
#
# All stresses are measured in MPa. The initial stress is consistent with
# the weight force from density 2500 kg/m^3, ie, stress_zz = 0.025*z MPa
# where gravity = 10 m.s^-2 = 1E-5 MPa m^2/kg. The maximum and minimum
# principal horizontal stresses are assumed to be equal to 0.8*stress_zz.
#
# Material properties:
# Young's modulus = 8 GPa
# Poisson's ratio = 0.25
# Cosserat layer thickness = 1 m
# Cosserat-joint normal stiffness = large
# Cosserat-joint shear stiffness = 1 GPa
# MC cohesion = 3 MPa
# MC friction angle = 37 deg
# MC dilation angle = 8 deg
# MC tensile strength = 1 MPa
# MC compressive strength = 100 MPa
# WeakPlane cohesion = 0.1 MPa
# WeakPlane friction angle = 30 deg
# WeakPlane dilation angle = 10 deg
# WeakPlane tensile strength = 0.1 MPa
# WeakPlane compressive strength = 100 MPa softening to 1 MPa at strain = 1
#
[Mesh]
[file]
type = FileMeshGenerator
file = mesh/fine.e
[]
[./xmin]
input = file
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
new_boundary = xmin
normal = '-1 0 0'
[../]
[./xmax]
input = xmin
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
new_boundary = xmax
normal = '1 0 0'
[../]
[./ymin]
input = xmax
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
new_boundary = ymin
normal = '0 -1 0'
[../]
[./ymax]
input = ymin
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
new_boundary = ymax
normal = '0 1 0'
[../]
[./zmax]
input = ymax
type = SideSetsAroundSubdomainGenerator
block = 30
new_boundary = zmax
normal = '0 0 1'
[../]
[./zmin]
input = zmax
type = SideSetsAroundSubdomainGenerator
block = 2
new_boundary = zmin
normal = '0 0 -1'
[../]
[./excav]
type = SubdomainBoundingBoxGenerator
input = zmin
block_id = 1
bottom_left = '0 0 -400'
top_right = '150 1000 -397'
[../]
[./roof]
type = SideSetsAroundSubdomainGenerator
block = 1
input = excav
new_boundary = roof
normal = '0 0 1'
[../]
[]
[GlobalParams]
perform_finite_strain_rotations = false
displacements = 'disp_x disp_y disp_z'
Cosserat_rotations = 'wc_x wc_y wc_z'
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./wc_x]
[../]
[./wc_y]
[../]
[]
[Kernels]
[./cx_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_x
component = 0
[../]
[./cy_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_y
component = 1
[../]
[./cz_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_z
component = 2
[../]
[./x_couple]
type = StressDivergenceTensors
use_displaced_mesh = false
variable = wc_x
displacements = 'wc_x wc_y wc_z'
component = 0
base_name = couple
[../]
[./y_couple]
type = StressDivergenceTensors
use_displaced_mesh = false
variable = wc_y
displacements = 'wc_x wc_y wc_z'
component = 1
base_name = couple
[../]
[./x_moment]
type = MomentBalancing
use_displaced_mesh = false
variable = wc_x
component = 0
[../]
[./y_moment]
type = MomentBalancing
use_displaced_mesh = false
variable = wc_y
component = 1
[../]
[./gravity]
type = Gravity
use_displaced_mesh = false
variable = disp_z
value = -10E-6 # remember this is in MPa
[../]
[]
[AuxVariables]
[./wc_z]
[../]
[./stress_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xz]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yz]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./mc_shear]
order = CONSTANT
family = MONOMIAL
[../]
[./mc_tensile]
order = CONSTANT
family = MONOMIAL
[../]
[./wp_shear]
order = CONSTANT
family = MONOMIAL
[../]
[./wp_tensile]
order = CONSTANT
family = MONOMIAL
[../]
[./wp_shear_f]
order = CONSTANT
family = MONOMIAL
[../]
[./wp_tensile_f]
order = CONSTANT
family = MONOMIAL
[../]
[./mc_shear_f]
order = CONSTANT
family = MONOMIAL
[../]
[./mc_tensile_f]
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_yx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yx
index_i = 1
index_j = 0
[../]
[./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_zx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zx
index_i = 2
index_j = 0
[../]
[./stress_zy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zy
index_i = 2
index_j = 1
[../]
[./stress_zz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zz
index_i = 2
index_j = 2
[../]
[./mc_shear]
type = MaterialStdVectorAux
index = 0
property = mc_plastic_internal_parameter
variable = mc_shear
[../]
[./mc_tensile]
type = MaterialStdVectorAux
index = 1
property = mc_plastic_internal_parameter
variable = mc_tensile
[../]
[./wp_shear]
type = MaterialStdVectorAux
index = 0
property = wp_plastic_internal_parameter
variable = wp_shear
[../]
[./wp_tensile]
type = MaterialStdVectorAux
index = 1
property = wp_plastic_internal_parameter
variable = wp_tensile
[../]
[./mc_shear_f]
type = MaterialStdVectorAux
index = 6
property = mc_plastic_yield_function
variable = mc_shear_f
[../]
[./mc_tensile_f]
type = MaterialStdVectorAux
index = 0
property = mc_plastic_yield_function
variable = mc_tensile_f
[../]
[./wp_shear_f]
type = MaterialStdVectorAux
index = 0
property = wp_plastic_yield_function
variable = wp_shear_f
[../]
[./wp_tensile_f]
type = MaterialStdVectorAux
index = 1
property = wp_plastic_yield_function
variable = wp_tensile_f
[../]
[]
[BCs]
[./no_x]
type = DirichletBC
variable = disp_x
boundary = 'xmin xmax'
value = 0.0
[../]
[./no_y]
type = DirichletBC
variable = disp_y
boundary = 'ymin ymax'
value = 0.0
[../]
[./no_z]
type = DirichletBC
variable = disp_z
boundary = zmin
value = 0.0
[../]
[./no_wc_x]
type = DirichletBC
variable = wc_x
boundary = 'ymin ymax'
value = 0.0
[../]
[./no_wc_y]
type = DirichletBC
variable = wc_y
boundary = 'xmin xmax'
value = 0.0
[../]
[./roof]
type = StickyBC
variable = disp_z
min_value = -3.0
boundary = roof
[../]
[]
[Functions]
[./ini_xx]
type = ParsedFunction
value = '0.8*2500*10E-6*z'
[../]
[./ini_zz]
type = ParsedFunction
value = '2500*10E-6*z'
[../]
[./excav_sideways]
type = ParsedFunction
vars = 'end_t ymin ymax minval maxval slope'
vals = '100.0 0 1000.0 1E-9 1 10'
# excavation face at ymin+(ymax-ymin)*min(t/end_t,1)
# slope is the distance over which the modulus reduces from maxval to minval
value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,if(y<ymin+(ymax-ymin)*min(t/end_t,1)+slope,minval+(maxval-minval)*(y-(ymin+(ymax-ymin)*min(t/end_t,1)))/slope,maxval))'
[../]
[./density_sideways]
type = ParsedFunction
vars = 'end_t ymin ymax minval maxval'
vals = '100.0 0 1000.0 0 2500'
value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,maxval)'
[../]
[]
[UserObjects]
[./mc_coh_strong_harden]
type = TensorMechanicsHardeningExponential
value_0 = 2.99 # MPa
value_residual = 3.01 # MPa
rate = 1.0
[../]
[./mc_fric]
type = TensorMechanicsHardeningConstant
value = 0.65 # 37deg
[../]
[./mc_dil]
type = TensorMechanicsHardeningConstant
value = 0.15 # 8deg
[../]
[./mc_tensile_str_strong_harden]
type = TensorMechanicsHardeningExponential
value_0 = 1.0 # MPa
value_residual = 1.0 # MPa
rate = 1.0
[../]
[./mc_compressive_str]
type = TensorMechanicsHardeningCubic
value_0 = 100 # Large!
value_residual = 100
internal_limit = 0.1
[../]
[./wp_coh_harden]
type = TensorMechanicsHardeningCubic
value_0 = 0.1
value_residual = 0.1
internal_limit = 10
[../]
[./wp_tan_fric]
type = TensorMechanicsHardeningConstant
value = 0.36 # 20deg
[../]
[./wp_tan_dil]
type = TensorMechanicsHardeningConstant
value = 0.18 # 10deg
[../]
[./wp_tensile_str_harden]
type = TensorMechanicsHardeningCubic
value_0 = 0.1
value_residual = 0.1
internal_limit = 10
[../]
[./wp_compressive_str_soften]
type = TensorMechanicsHardeningCubic
value_0 = 100
value_residual = 1
internal_limit = 1.0
[../]
[]
[Materials]
[./elasticity_tensor_0]
type = ComputeLayeredCosseratElasticityTensor
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
young = 8E3 # MPa
poisson = 0.25
layer_thickness = 1.0
joint_normal_stiffness = 1E9 # huge
joint_shear_stiffness = 1E3 # MPa
[../]
[./elasticity_tensor_1]
type = ComputeLayeredCosseratElasticityTensor
block = 1
young = 8E3 # MPa
poisson = 0.25
layer_thickness = 1.0
joint_normal_stiffness = 1E9 # huge
joint_shear_stiffness = 1E3 # MPa
elasticity_tensor_prefactor = excav_sideways
[../]
[./strain]
type = ComputeCosseratIncrementalSmallStrain
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
eigenstrain_name = ini_stress
initial_stress = 'ini_xx 0 0 0 ini_xx 0 0 0 ini_zz'
[../]
[./stress_0]
type = ComputeMultipleInelasticCosseratStress
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
inelastic_models = 'mc wp'
cycle_models = true
relative_tolerance = 2.0
absolute_tolerance = 1E6
max_iterations = 1
tangent_operator = nonlinear
perform_finite_strain_rotations = false
[../]
[./stress_1]
type = ComputeMultipleInelasticCosseratStress
block = 1
inelastic_models = ''
relative_tolerance = 2.0
absolute_tolerance = 1E6
max_iterations = 1
tangent_operator = nonlinear
perform_finite_strain_rotations = false
[../]
[./mc]
type = CappedMohrCoulombCosseratStressUpdate
warn_about_precision_loss = false
host_youngs_modulus = 8E3
host_poissons_ratio = 0.25
base_name = mc
tensile_strength = mc_tensile_str_strong_harden
compressive_strength = mc_compressive_str
cohesion = mc_coh_strong_harden
friction_angle = mc_fric
dilation_angle = mc_dil
max_NR_iterations = 100000
smoothing_tol = 0.1 # MPa # Must be linked to cohesion
yield_function_tol = 1E-9 # MPa. this is essentially the lowest possible without lots of precision loss
perfect_guess = true
min_step_size = 1.0
[../]
[./wp]
type = CappedWeakPlaneCosseratStressUpdate
warn_about_precision_loss = false
base_name = wp
cohesion = wp_coh_harden
tan_friction_angle = wp_tan_fric
tan_dilation_angle = wp_tan_dil
tensile_strength = wp_tensile_str_harden
compressive_strength = wp_compressive_str_soften
max_NR_iterations = 10000
tip_smoother = 0.1
smoothing_tol = 0.1 # MPa # Note, this must be tied to cohesion, otherwise get no possible return at cone apex
yield_function_tol = 1E-11 # MPa. this is essentially the lowest possible without lots of precision loss
perfect_guess = true
min_step_size = 1.0E-3
[../]
[./density_0]
type = GenericConstantMaterial
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
prop_names = density
prop_values = 2500
[../]
[./density_1]
type = GenericFunctionMaterial
block = 1
prop_names = density
prop_values = density_sideways
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[]
[]
[Postprocessors]
[./min_roof_disp]
type = NodalExtremeValue
boundary = roof
value_type = min
variable = disp_z
[../]
[./min_surface_disp]
type = NodalExtremeValue
boundary = zmax
value_type = min
variable = disp_z
[../]
[]
[Executioner]
type = Transient
solve_type = 'NEWTON'
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-pc_type -ksp_type -ksp_gmres_restart'
petsc_options_value = ' bjacobi gmres 200'
line_search = bt
nl_abs_tol = 1e-3
nl_rel_tol = 1e-5
l_max_its = 30
nl_max_its = 1000
start_time = 0.0
dt = 0.5
end_time = 100.0
[]
[Outputs]
interval = 1
print_linear_residuals = false
exodus = true
csv = true
console = true
[]
modules/solid_mechanics/test/tests/generalized_plane_strain/generalized_plane_strain_field.i
[GlobalParams]
order = FIRST
family = LAGRANGE
[]
[Mesh]
file = square.e
displacements = 'disp_x disp_y'
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./strain_zz]
[../]
[]
[AuxVariables]
[./temp]
[../]
[./saved_x]
[../]
[./saved_y]
[../]
[./saved_z]
[../]
[./stress_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_xy]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./aux_strain_zz]
order = CONSTANT
family = MONOMIAL
[../]
[]
[Postprocessors]
[./react_z]
type = MaterialTensorIntegralSM
tensor = stress
index = 2
[../]
[./min_strain_zz]
type = NodalExtremeValue
variable = strain_zz
value_type = min
[../]
[./max_strain_zz]
type = NodalExtremeValue
variable = strain_zz
value_type = max
[../]
[]
[SolidMechanics]
[./solid]
disp_x = disp_x
disp_y = disp_y
save_in_disp_x = saved_x
save_in_disp_y = saved_y
temp = temp
[../]
[]
[Kernels]
[./solid_z]
type = OutOfPlaneStress
variable = strain_zz
save_in = saved_z
disp_x = disp_x
disp_y = disp_y
temp = temp
[../]
[]
[Constraints]
[./szz]
type = EqualValueBoundaryConstraint
variable = strain_zz
master = '8'
slave = 10
penalty = 1e12
[../]
[]
[AuxKernels]
[./tempfuncaux]
type = FunctionAux
variable = temp
function = tempfunc
use_displaced_mesh = false
[../]
[./stress_xx]
type = MaterialTensorAux
tensor = stress
variable = stress_xx
index = 0
[../]
[./stress_xy]
type = MaterialTensorAux
tensor = stress
variable = stress_xy
index = 3
[../]
[./stress_yy]
type = MaterialTensorAux
tensor = stress
variable = stress_yy
index = 1
[../]
[./stress_zz]
type = MaterialTensorAux
tensor = stress
variable = stress_zz
index = 2
[../]
[./strain_xx]
type = MaterialTensorAux
tensor = total_strain
variable = strain_xx
index = 0
[../]
[./strain_xy]
type = MaterialTensorAux
tensor = total_strain
variable = strain_xy
index = 3
[../]
[./strain_yy]
type = MaterialTensorAux
tensor = total_strain
variable = strain_yy
index = 1
[../]
[./strain_zz]
type = MaterialTensorAux
tensor = total_strain
variable = aux_strain_zz
index = 2
[../]
[]
[Functions]
[./pull]
type = PiecewiseLinear
x='0 1 100'
y='0 0.00 0.00'
[../]
[./tempfunc]
type = ParsedFunction
value = '(1-x)*t'
[../]
[]
[BCs]
[./bottomx]
type = DirichletBC
boundary = 1
variable = disp_x
value = 0.0
[../]
[./bottomy]
type = DirichletBC
boundary = 1
variable = disp_y
value = 0.0
[../]
[]
[Materials]
[./linelast]
type = Elastic
block = 1
disp_x = disp_x
disp_y = disp_y
poissons_ratio = 0.3
youngs_modulus = 1e6
thermal_expansion = 0.02
stress_free_temperature = 0.5
temp = temp
formulation = PlaneStrain
strain_zz = strain_zz
[../]
[]
[Executioner]
type = Transient
solve_type = PJFNK
line_search = none
# controls for linear iterations
l_max_its = 100
l_tol = 1e-4
# controls for nonlinear iterations
nl_max_its = 15
nl_rel_tol = 1e-10
nl_abs_tol = 1e-5
# time control
start_time = 0.0
dt = 1.0
dtmin = 1.0
end_time = 2.0
num_steps = 5000
[]
[Outputs]
exodus = true
[]
modules/porous_flow/examples/coal_mining/coarse_with_fluid.i
# Strata deformation and fluid flow aaround a coal mine - 3D model
#
# A "half model" is used. The mine is 400m deep and
# just the roof is studied (-400<=z<=0). The mining panel
# sits between 0<=x<=150, and 0<=y<=1000, so this simulates
# a coal panel that is 300m wide and 1000m long. The outer boundaries
# are 1km from the excavation boundaries.
#
# The excavation takes 0.5 years.
#
# The boundary conditions for this simulation are:
# - disp_x = 0 at x=0 and x=1150
# - disp_y = 0 at y=-1000 and y=1000
# - disp_z = 0 at z=-400, but there is a time-dependent
# Young modulus that simulates excavation
# - wc_x = 0 at y=-1000 and y=1000
# - wc_y = 0 at x=0 and x=1150
# - no flow at x=0, z=-400 and z=0
# - fixed porepressure at y=-1000, y=1000 and x=1150
# That is, rollers on the sides, free at top,
# and prescribed at bottom in the unexcavated portion.
#
# A single-phase unsaturated fluid is used.
#
# The small strain formulation is used.
#
# All stresses are measured in MPa, and time units are measured in years.
#
# The initial porepressure is hydrostatic with P=0 at z=0, so
# Porepressure ~ - 0.01*z MPa, where the fluid has density 1E3 kg/m^3 and
# gravity = = 10 m.s^-2 = 1E-5 MPa m^2/kg.
# To be more accurate, i use
# Porepressure = -bulk * log(1 + g*rho0*z/bulk)
# where bulk=2E3 MPa and rho0=1Ee kg/m^3.
# The initial stress is consistent with the weight force from undrained
# density 2500 kg/m^3, and fluid porepressure, and a Biot coefficient of 0.7, ie,
# stress_zz^effective = 0.025*z + 0.7 * initial_porepressure
# The maximum and minimum principal horizontal effective stresses are
# assumed to be equal to 0.8*stress_zz.
#
# Material properties:
# Young's modulus = 8 GPa
# Poisson's ratio = 0.25
# Cosserat layer thickness = 1 m
# Cosserat-joint normal stiffness = large
# Cosserat-joint shear stiffness = 1 GPa
# MC cohesion = 2 MPa
# MC friction angle = 35 deg
# MC dilation angle = 8 deg
# MC tensile strength = 1 MPa
# MC compressive strength = 100 MPa
# WeakPlane cohesion = 0.1 MPa
# WeakPlane friction angle = 30 deg
# WeakPlane dilation angle = 10 deg
# WeakPlane tensile strength = 0.1 MPa
# WeakPlane compressive strength = 100 MPa softening to 1 MPa at strain = 1
# Fluid density at zero porepressure = 1E3 kg/m^3
# Fluid bulk modulus = 2E3 MPa
# Fluid viscosity = 1.1E-3 Pa.s = 1.1E-9 MPa.s = 3.5E-17 MPa.year
#
[GlobalParams]
perform_finite_strain_rotations = false
displacements = 'disp_x disp_y disp_z'
Cosserat_rotations = 'wc_x wc_y wc_z'
PorousFlowDictator = dictator
biot_coefficient = 0.7
[]
[Mesh]
[file]
type = FileMeshGenerator
file = mesh/coarse.e
[]
[./xmin]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
new_boundary = xmin
normal = '-1 0 0'
input = file
[../]
[./xmax]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
new_boundary = xmax
normal = '1 0 0'
input = xmin
[../]
[./ymin]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
new_boundary = ymin
normal = '0 -1 0'
input = xmax
[../]
[./ymax]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
new_boundary = ymax
normal = '0 1 0'
input = ymin
[../]
[./zmax]
type = SideSetsAroundSubdomainGenerator
block = 16
new_boundary = zmax
normal = '0 0 1'
input = ymax
[../]
[./zmin]
type = SideSetsAroundSubdomainGenerator
block = 2
new_boundary = zmin
normal = '0 0 -1'
input = zmax
[../]
[./excav]
type = SubdomainBoundingBoxGenerator
input = zmin
block_id = 1
bottom_left = '0 0 -400'
top_right = '150 1000 -397'
[../]
[./roof]
type = SideSetsBetweenSubdomainsGenerator
master_block = 3
paired_block = 1
input = excav
new_boundary = roof
[../]
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./wc_x]
[../]
[./wc_y]
[../]
[./porepressure]
scaling = 1E-5
[../]
[]
[ICs]
[./porepressure]
type = FunctionIC
variable = porepressure
function = ini_pp
[../]
[]
[Kernels]
[./cx_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_x
component = 0
[../]
[./cy_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_y
component = 1
[../]
[./cz_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_z
component = 2
[../]
[./x_couple]
type = StressDivergenceTensors
use_displaced_mesh = false
variable = wc_x
displacements = 'wc_x wc_y wc_z'
component = 0
base_name = couple
[../]
[./y_couple]
type = StressDivergenceTensors
use_displaced_mesh = false
variable = wc_y
displacements = 'wc_x wc_y wc_z'
component = 1
base_name = couple
[../]
[./x_moment]
type = MomentBalancing
use_displaced_mesh = false
variable = wc_x
component = 0
[../]
[./y_moment]
type = MomentBalancing
use_displaced_mesh = false
variable = wc_y
component = 1
[../]
[./gravity]
type = Gravity
use_displaced_mesh = false
variable = disp_z
value = -10E-6 # remember this is in MPa
[../]
[./poro_x]
type = PorousFlowEffectiveStressCoupling
use_displaced_mesh = false
variable = disp_x
component = 0
[../]
[./poro_y]
type = PorousFlowEffectiveStressCoupling
use_displaced_mesh = false
variable = disp_y
component = 1
[../]
[./poro_z]
type = PorousFlowEffectiveStressCoupling
use_displaced_mesh = false
component = 2
variable = disp_z
[../]
[./mass0]
type = PorousFlowMassTimeDerivative
use_displaced_mesh = false
fluid_component = 0
variable = porepressure
[../]
[./flux]
type = PorousFlowAdvectiveFlux
use_displaced_mesh = false
variable = porepressure
gravity = '0 0 -10E-6'
fluid_component = 0
[../]
[./poro_vol_exp]
type = PorousFlowMassVolumetricExpansion
use_displaced_mesh = false
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
variable = porepressure
fluid_component = 0
[../]
[]
[AuxVariables]
[./saturation]
order = CONSTANT
family = MONOMIAL
[../]
[./darcy_x]
order = CONSTANT
family = MONOMIAL
[../]
[./darcy_y]
order = CONSTANT
family = MONOMIAL
[../]
[./darcy_z]
order = CONSTANT
family = MONOMIAL
[../]
[./porosity]
order = CONSTANT
family = MONOMIAL
[../]
[./wc_z]
[../]
[./stress_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xz]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yz]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_xy]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_xz]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_yx]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_yz]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_zx]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_zy]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./perm_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./perm_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./perm_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./mc_shear]
order = CONSTANT
family = MONOMIAL
[../]
[./mc_tensile]
order = CONSTANT
family = MONOMIAL
[../]
[./wp_shear]
order = CONSTANT
family = MONOMIAL
[../]
[./wp_tensile]
order = CONSTANT
family = MONOMIAL
[../]
[./wp_shear_f]
order = CONSTANT
family = MONOMIAL
[../]
[./wp_tensile_f]
order = CONSTANT
family = MONOMIAL
[../]
[./mc_shear_f]
order = CONSTANT
family = MONOMIAL
[../]
[./mc_tensile_f]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./saturation_water]
type = PorousFlowPropertyAux
variable = saturation
property = saturation
phase = 0
execute_on = timestep_end
[../]
[./darcy_x]
type = PorousFlowDarcyVelocityComponent
variable = darcy_x
gravity = '0 0 -10E-6'
component = x
[../]
[./darcy_y]
type = PorousFlowDarcyVelocityComponent
variable = darcy_y
gravity = '0 0 -10E-6'
component = y
[../]
[./darcy_z]
type = PorousFlowDarcyVelocityComponent
variable = darcy_z
gravity = '0 0 -10E-6'
component = z
[../]
[./porosity]
type = PorousFlowPropertyAux
property = porosity
variable = porosity
execute_on = timestep_end
[../]
[./stress_xx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx
index_i = 0
index_j = 0
execute_on = timestep_end
[../]
[./stress_xy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xy
index_i = 0
index_j = 1
execute_on = timestep_end
[../]
[./stress_xz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xz
index_i = 0
index_j = 2
execute_on = timestep_end
[../]
[./stress_yx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yx
index_i = 1
index_j = 0
execute_on = timestep_end
[../]
[./stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
execute_on = timestep_end
[../]
[./stress_yz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yz
index_i = 1
index_j = 2
execute_on = timestep_end
[../]
[./stress_zx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zx
index_i = 2
index_j = 0
execute_on = timestep_end
[../]
[./stress_zy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zy
index_i = 2
index_j = 1
execute_on = timestep_end
[../]
[./stress_zz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zz
index_i = 2
index_j = 2
execute_on = timestep_end
[../]
[./total_strain_xx]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_xx
index_i = 0
index_j = 0
execute_on = timestep_end
[../]
[./total_strain_xy]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_xy
index_i = 0
index_j = 1
execute_on = timestep_end
[../]
[./total_strain_xz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_xz
index_i = 0
index_j = 2
execute_on = timestep_end
[../]
[./total_strain_yx]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_yx
index_i = 1
index_j = 0
execute_on = timestep_end
[../]
[./total_strain_yy]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_yy
index_i = 1
index_j = 1
execute_on = timestep_end
[../]
[./total_strain_yz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_yz
index_i = 1
index_j = 2
execute_on = timestep_end
[../]
[./total_strain_zx]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_zx
index_i = 2
index_j = 0
execute_on = timestep_end
[../]
[./total_strain_zy]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_zy
index_i = 2
index_j = 1
execute_on = timestep_end
[../]
[./total_strain_zz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_zz
index_i = 2
index_j = 2
execute_on = timestep_end
[../]
[./perm_xx]
type = PorousFlowPropertyAux
property = permeability
variable = perm_xx
row = 0
column = 0
execute_on = timestep_end
[../]
[./perm_yy]
type = PorousFlowPropertyAux
property = permeability
variable = perm_yy
row = 1
column = 1
execute_on = timestep_end
[../]
[./perm_zz]
type = PorousFlowPropertyAux
property = permeability
variable = perm_zz
row = 2
column = 2
execute_on = timestep_end
[../]
[./mc_shear]
type = MaterialStdVectorAux
index = 0
property = mc_plastic_internal_parameter
variable = mc_shear
execute_on = timestep_end
[../]
[./mc_tensile]
type = MaterialStdVectorAux
index = 1
property = mc_plastic_internal_parameter
variable = mc_tensile
execute_on = timestep_end
[../]
[./wp_shear]
type = MaterialStdVectorAux
index = 0
property = wp_plastic_internal_parameter
variable = wp_shear
execute_on = timestep_end
[../]
[./wp_tensile]
type = MaterialStdVectorAux
index = 1
property = wp_plastic_internal_parameter
variable = wp_tensile
execute_on = timestep_end
[../]
[./mc_shear_f]
type = MaterialStdVectorAux
index = 6
property = mc_plastic_yield_function
variable = mc_shear_f
execute_on = timestep_end
[../]
[./mc_tensile_f]
type = MaterialStdVectorAux
index = 0
property = mc_plastic_yield_function
variable = mc_tensile_f
execute_on = timestep_end
[../]
[./wp_shear_f]
type = MaterialStdVectorAux
index = 0
property = wp_plastic_yield_function
variable = wp_shear_f
execute_on = timestep_end
[../]
[./wp_tensile_f]
type = MaterialStdVectorAux
index = 1
property = wp_plastic_yield_function
variable = wp_tensile_f
execute_on = timestep_end
[../]
[]
[BCs]
[./no_x]
type = DirichletBC
variable = disp_x
boundary = 'xmin xmax'
value = 0.0
[../]
[./no_y]
type = DirichletBC
variable = disp_y
boundary = 'ymin ymax'
value = 0.0
[../]
[./no_z]
type = DirichletBC
variable = disp_z
boundary = zmin
value = 0.0
[../]
[./no_wc_x]
type = DirichletBC
variable = wc_x
boundary = 'ymin ymax'
value = 0.0
[../]
[./no_wc_y]
type = DirichletBC
variable = wc_y
boundary = 'xmin xmax'
value = 0.0
[../]
[./fix_porepressure]
type = FunctionDirichletBC
variable = porepressure
boundary = 'ymin ymax xmax'
function = ini_pp
[../]
[./roof_porepressure]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
pt_vals = '-1E3 1E3'
multipliers = '-1 1'
fluid_phase = 0
flux_function = roof_conductance
boundary = roof
[../]
[./roof_bcs]
type = StickyBC
variable = disp_z
min_value = -3.0
boundary = roof
[../]
[]
[Functions]
[./ini_pp]
type = ParsedFunction
vars = 'bulk p0 g rho0'
vals = '2E3 0.0 1E-5 1E3'
value = '-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk)'
[../]
[./ini_xx]
type = ParsedFunction
vars = 'bulk p0 g rho0 biot'
vals = '2E3 0.0 1E-5 1E3 0.7'
value = '0.8*(2500*10E-6*z+biot*(-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk)))'
[../]
[./ini_zz]
type = ParsedFunction
vars = 'bulk p0 g rho0 biot'
vals = '2E3 0.0 1E-5 1E3 0.7'
value = '2500*10E-6*z+biot*(-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk))'
[../]
[./excav_sideways]
type = ParsedFunction
vars = 'end_t ymin ymax minval maxval slope'
vals = '0.5 0 1000.0 1E-9 1 60'
# excavation face at ymin+(ymax-ymin)*min(t/end_t,1)
# slope is the distance over which the modulus reduces from maxval to minval
value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,if(y<ymin+(ymax-ymin)*min(t/end_t,1)+slope,minval+(maxval-minval)*(y-(ymin+(ymax-ymin)*min(t/end_t,1)))/slope,maxval))'
[../]
[./density_sideways]
type = ParsedFunction
vars = 'end_t ymin ymax minval maxval'
vals = '0.5 0 1000.0 0 2500'
value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,maxval)'
[../]
[./roof_conductance]
type = ParsedFunction
vars = 'end_t ymin ymax maxval minval'
vals = '0.5 0 1000.0 1E7 0'
value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),maxval,minval)'
[../]
[]
[UserObjects]
[./dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure disp_x disp_y disp_z'
number_fluid_phases = 1
number_fluid_components = 1
[../]
[./pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1 # MPa^-1
[../]
[./mc_coh_strong_harden]
type = TensorMechanicsHardeningExponential
value_0 = 1.99 # MPa
value_residual = 2.01 # MPa
rate = 1.0
[../]
[./mc_fric]
type = TensorMechanicsHardeningConstant
value = 0.61 # 35deg
[../]
[./mc_dil]
type = TensorMechanicsHardeningConstant
value = 0.15 # 8deg
[../]
[./mc_tensile_str_strong_harden]
type = TensorMechanicsHardeningExponential
value_0 = 1.0 # MPa
value_residual = 1.0 # MPa
rate = 1.0
[../]
[./mc_compressive_str]
type = TensorMechanicsHardeningCubic
value_0 = 100 # Large!
value_residual = 100
internal_limit = 0.1
[../]
[./wp_coh_harden]
type = TensorMechanicsHardeningCubic
value_0 = 0.05
value_residual = 0.05
internal_limit = 10
[../]
[./wp_tan_fric]
type = TensorMechanicsHardeningConstant
value = 0.26 # 15deg
[../]
[./wp_tan_dil]
type = TensorMechanicsHardeningConstant
value = 0.18 # 10deg
[../]
[./wp_tensile_str_harden]
type = TensorMechanicsHardeningCubic
value_0 = 0.05
value_residual = 0.05
internal_limit = 10
[../]
[./wp_compressive_str_soften]
type = TensorMechanicsHardeningCubic
value_0 = 100
value_residual = 1
internal_limit = 1.0
[../]
[]
[Modules]
[./FluidProperties]
[./simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E3
density0 = 1000
thermal_expansion = 0
viscosity = 3.5E-17
[../]
[../]
[]
[Materials]
[./temperature]
type = PorousFlowTemperature
[../]
[./eff_fluid_pressure]
type = PorousFlowEffectiveFluidPressure
[../]
[./vol_strain]
type = PorousFlowVolumetricStrain
[../]
[./ppss]
type = PorousFlow1PhaseP
porepressure = porepressure
capillary_pressure = pc
[../]
[./massfrac]
type = PorousFlowMassFraction
[../]
[./simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[../]
[./porosity_bulk]
type = PorousFlowPorosity
fluid = true
mechanical = true
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
ensure_positive = true
porosity_zero = 0.02
solid_bulk = 5.3333E3
[../]
[./porosity_excav]
type = PorousFlowPorosityConst
block = 1
porosity = 1.0
[../]
[./permeability_bulk]
type = PorousFlowPermeabilityKozenyCarman
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
poroperm_function = kozeny_carman_phi0
k0 = 1E-15
phi0 = 0.02
n = 2
m = 2
[../]
[./permeability_excav]
type = PorousFlowPermeabilityConst
block = 1
permeability = '0 0 0 0 0 0 0 0 0'
[../]
[./relperm]
type = PorousFlowRelativePermeabilityCorey
n = 4
s_res = 0.4
sum_s_res = 0.4
phase = 0
[../]
[./elasticity_tensor_0]
type = ComputeLayeredCosseratElasticityTensor
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
young = 8E3 # MPa
poisson = 0.25
layer_thickness = 1.0
joint_normal_stiffness = 1E9 # huge
joint_shear_stiffness = 1E3 # MPa
[../]
[./elasticity_tensor_1]
type = ComputeLayeredCosseratElasticityTensor
block = 1
young = 8E3 # MPa
poisson = 0.25
layer_thickness = 1.0
joint_normal_stiffness = 1E9 # huge
joint_shear_stiffness = 1E3 # MPa
elasticity_tensor_prefactor = excav_sideways
[../]
[./strain]
type = ComputeCosseratIncrementalSmallStrain
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
eigenstrain_name = ini_stress
initial_stress = 'ini_xx 0 0 0 ini_xx 0 0 0 ini_zz'
[../]
[./stress_0]
type = ComputeMultipleInelasticCosseratStress
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
inelastic_models = 'mc wp'
cycle_models = true
relative_tolerance = 2.0
absolute_tolerance = 1E6
max_iterations = 1
tangent_operator = nonlinear
perform_finite_strain_rotations = false
[../]
[./stress_1]
type = ComputeMultipleInelasticCosseratStress
block = 1
inelastic_models = ''
relative_tolerance = 2.0
absolute_tolerance = 1E6
max_iterations = 1
tangent_operator = nonlinear
perform_finite_strain_rotations = false
[../]
[./mc]
type = CappedMohrCoulombCosseratStressUpdate
warn_about_precision_loss = false
host_youngs_modulus = 8E3
host_poissons_ratio = 0.25
base_name = mc
tensile_strength = mc_tensile_str_strong_harden
compressive_strength = mc_compressive_str
cohesion = mc_coh_strong_harden
friction_angle = mc_fric
dilation_angle = mc_dil
max_NR_iterations = 100000
smoothing_tol = 0.1 # MPa # Must be linked to cohesion
yield_function_tol = 1E-9 # MPa. this is essentially the lowest possible without lots of precision loss
perfect_guess = true
min_step_size = 1.0
[../]
[./wp]
type = CappedWeakPlaneCosseratStressUpdate
warn_about_precision_loss = false
base_name = wp
cohesion = wp_coh_harden
tan_friction_angle = wp_tan_fric
tan_dilation_angle = wp_tan_dil
tensile_strength = wp_tensile_str_harden
compressive_strength = wp_compressive_str_soften
max_NR_iterations = 10000
tip_smoother = 0.05
smoothing_tol = 0.05 # MPa # Note, this must be tied to cohesion, otherwise get no possible return at cone apex
yield_function_tol = 1E-11 # MPa. this is essentially the lowest possible without lots of precision loss
perfect_guess = true
min_step_size = 1.0E-3
[../]
[./undrained_density_0]
type = GenericConstantMaterial
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
prop_names = density
prop_values = 2500
[../]
[./undrained_density_1]
type = GenericFunctionMaterial
block = 1
prop_names = density
prop_values = density_sideways
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[]
[]
[Postprocessors]
[./min_roof_disp]
type = NodalExtremeValue
boundary = roof
value_type = min
variable = disp_z
[../]
[./min_roof_pp]
type = NodalExtremeValue
boundary = roof
value_type = min
variable = porepressure
[../]
[./min_surface_disp]
type = NodalExtremeValue
boundary = zmax
value_type = min
variable = disp_z
[../]
[./min_surface_pp]
type = NodalExtremeValue
boundary = zmax
value_type = min
variable = porepressure
[../]
[./max_perm_zz]
type = ElementExtremeValue
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
variable = perm_zz
[../]
[]
[Executioner]
type = Transient
solve_type = 'NEWTON'
petsc_options = '-snes_converged_reason'
# best overall
# petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
# petsc_options_value = ' lu mumps'
# best if you do not have mumps:
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu superlu_dist'
# best if you do not have mumps or superlu_dist:
#petsc_options_iname = '-pc_type -pc_asm_overlap -sub_pc_type -ksp_type -ksp_gmres_restart'
#petsc_options_value = ' asm 2 lu gmres 200'
# very basic:
#petsc_options_iname = '-pc_type -ksp_type -ksp_gmres_restart'
#petsc_options_value = ' bjacobi gmres 200'
line_search = bt
nl_abs_tol = 1e-3
nl_rel_tol = 1e-5
l_max_its = 200
nl_max_its = 30
start_time = 0.0
dt = 0.014706
end_time = 0.014706 #0.5
[]
[Outputs]
interval = 1
print_linear_residuals = true
exodus = true
csv = true
console = true
[]
modules/porous_flow/examples/coal_mining/fine_with_fluid.i
#################################################################
#
# NOTE:
# The mesh for this model is too large for the MOOSE repository
# so is kept in the the large_media submodule
#
#################################################################
#
# Strata deformation and fluid flow aaround a coal mine - 3D model
#
# A "half model" is used. The mine is 400m deep and
# just the roof is studied (-400<=z<=0). The mining panel
# sits between 0<=x<=150, and 0<=y<=1000, so this simulates
# a coal panel that is 300m wide and 1000m long. The outer boundaries
# are 1km from the excavation boundaries.
#
# The excavation takes 0.5 years.
#
# The boundary conditions for this simulation are:
# - disp_x = 0 at x=0 and x=1150
# - disp_y = 0 at y=-1000 and y=1000
# - disp_z = 0 at z=-400, but there is a time-dependent
# Young modulus that simulates excavation
# - wc_x = 0 at y=-1000 and y=1000
# - wc_y = 0 at x=0 and x=1150
# - no flow at x=0, z=-400 and z=0
# - fixed porepressure at y=-1000, y=1000 and x=1150
# That is, rollers on the sides, free at top,
# and prescribed at bottom in the unexcavated portion.
#
# A single-phase unsaturated fluid is used.
#
# The small strain formulation is used.
#
# All stresses are measured in MPa, and time units are measured in years.
#
# The initial porepressure is hydrostatic with P=0 at z=0, so
# Porepressure ~ - 0.01*z MPa, where the fluid has density 1E3 kg/m^3 and
# gravity = = 10 m.s^-2 = 1E-5 MPa m^2/kg.
# To be more accurate, i use
# Porepressure = -bulk * log(1 + g*rho0*z/bulk)
# where bulk=2E3 MPa and rho0=1Ee kg/m^3.
# The initial stress is consistent with the weight force from undrained
# density 2500 kg/m^3, and fluid porepressure, and a Biot coefficient of 0.7, ie,
# stress_zz^effective = 0.025*z + 0.7 * initial_porepressure
# The maximum and minimum principal horizontal effective stresses are
# assumed to be equal to 0.8*stress_zz.
#
# Material properties:
# Young's modulus = 8 GPa
# Poisson's ratio = 0.25
# Cosserat layer thickness = 1 m
# Cosserat-joint normal stiffness = large
# Cosserat-joint shear stiffness = 1 GPa
# MC cohesion = 2 MPa
# MC friction angle = 35 deg
# MC dilation angle = 8 deg
# MC tensile strength = 1 MPa
# MC compressive strength = 100 MPa
# WeakPlane cohesion = 0.1 MPa
# WeakPlane friction angle = 30 deg
# WeakPlane dilation angle = 10 deg
# WeakPlane tensile strength = 0.1 MPa
# WeakPlane compressive strength = 100 MPa softening to 1 MPa at strain = 1
# Fluid density at zero porepressure = 1E3 kg/m^3
# Fluid bulk modulus = 2E3 MPa
# Fluid viscosity = 1.1E-3 Pa.s = 1.1E-9 MPa.s = 3.5E-17 MPa.year
#
[GlobalParams]
perform_finite_strain_rotations = false
displacements = 'disp_x disp_y disp_z'
Cosserat_rotations = 'wc_x wc_y wc_z'
PorousFlowDictator = dictator
biot_coefficient = 0.7
[]
[Mesh]
[file]
type = FileMeshGenerator
file = fine.e
[]
[./xmin]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
new_boundary = xmin
normal = '-1 0 0'
input = file
[../]
[./xmax]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
new_boundary = xmax
normal = '1 0 0'
input = xmin
[../]
[./ymin]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
new_boundary = ymin
normal = '0 -1 0'
input = xmax
[../]
[./ymax]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
new_boundary = ymax
normal = '0 1 0'
input = ymin
[../]
[./zmax]
type = SideSetsAroundSubdomainGenerator
block = 30
new_boundary = zmax
normal = '0 0 1'
input = ymax
[../]
[./zmin]
type = SideSetsAroundSubdomainGenerator
block = 2
new_boundary = zmin
normal = '0 0 -1'
input = zmax
[../]
[./excav]
type = SubdomainBoundingBoxGenerator
input = zmin
block_id = 1
bottom_left = '0 0 -400'
top_right = '150 1000 -397'
[../]
[./roof]
type = SideSetsBetweenSubdomainsGenerator
master_block = 3
paired_block = 1
input = excav
new_boundary = roof
[../]
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./wc_x]
[../]
[./wc_y]
[../]
[./porepressure]
scaling = 1E-5
[../]
[]
[ICs]
[./porepressure]
type = FunctionIC
variable = porepressure
function = ini_pp
[../]
[]
[Kernels]
[./cx_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_x
component = 0
[../]
[./cy_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_y
component = 1
[../]
[./cz_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_z
component = 2
[../]
[./x_couple]
type = StressDivergenceTensors
use_displaced_mesh = false
variable = wc_x
displacements = 'wc_x wc_y wc_z'
component = 0
base_name = couple
[../]
[./y_couple]
type = StressDivergenceTensors
use_displaced_mesh = false
variable = wc_y
displacements = 'wc_x wc_y wc_z'
component = 1
base_name = couple
[../]
[./x_moment]
type = MomentBalancing
use_displaced_mesh = false
variable = wc_x
component = 0
[../]
[./y_moment]
type = MomentBalancing
use_displaced_mesh = false
variable = wc_y
component = 1
[../]
[./gravity]
type = Gravity
use_displaced_mesh = false
variable = disp_z
value = -10E-6 # remember this is in MPa
[../]
[./poro_x]
type = PorousFlowEffectiveStressCoupling
use_displaced_mesh = false
variable = disp_x
component = 0
[../]
[./poro_y]
type = PorousFlowEffectiveStressCoupling
use_displaced_mesh = false
variable = disp_y
component = 1
[../]
[./poro_z]
type = PorousFlowEffectiveStressCoupling
use_displaced_mesh = false
component = 2
variable = disp_z
[../]
[./poro_vol_exp]
type = PorousFlowMassVolumetricExpansion
use_displaced_mesh = false
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
variable = porepressure
fluid_component = 0
[../]
[./mass0]
type = PorousFlowMassTimeDerivative
use_displaced_mesh = false
fluid_component = 0
variable = porepressure
[../]
[./flux]
type = PorousFlowAdvectiveFlux
use_displaced_mesh = false
variable = porepressure
gravity = '0 0 -10E-6'
fluid_component = 0
[../]
[]
[AuxVariables]
[./saturation]
order = CONSTANT
family = MONOMIAL
[../]
[./darcy_x]
order = CONSTANT
family = MONOMIAL
[../]
[./darcy_y]
order = CONSTANT
family = MONOMIAL
[../]
[./darcy_z]
order = CONSTANT
family = MONOMIAL
[../]
[./porosity]
order = CONSTANT
family = MONOMIAL
[../]
[./wc_z]
[../]
[./stress_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xz]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yz]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_xy]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_xz]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_yx]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_yz]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_zx]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_zy]
order = CONSTANT
family = MONOMIAL
[../]
[./total_strain_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./perm_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./perm_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./perm_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./mc_shear]
order = CONSTANT
family = MONOMIAL
[../]
[./mc_tensile]
order = CONSTANT
family = MONOMIAL
[../]
[./wp_shear]
order = CONSTANT
family = MONOMIAL
[../]
[./wp_tensile]
order = CONSTANT
family = MONOMIAL
[../]
[./wp_shear_f]
order = CONSTANT
family = MONOMIAL
[../]
[./wp_tensile_f]
order = CONSTANT
family = MONOMIAL
[../]
[./mc_shear_f]
order = CONSTANT
family = MONOMIAL
[../]
[./mc_tensile_f]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./saturation_water]
type = PorousFlowPropertyAux
variable = saturation
property = saturation
phase = 0
execute_on = timestep_end
[../]
[./darcy_x]
type = PorousFlowDarcyVelocityComponent
variable = darcy_x
gravity = '0 0 -10E-6'
component = x
[../]
[./darcy_y]
type = PorousFlowDarcyVelocityComponent
variable = darcy_y
gravity = '0 0 -10E-6'
component = y
[../]
[./darcy_z]
type = PorousFlowDarcyVelocityComponent
variable = darcy_z
gravity = '0 0 -10E-6'
component = z
[../]
[./porosity]
type = PorousFlowPropertyAux
property = porosity
variable = porosity
execute_on = timestep_end
[../]
[./stress_xx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx
index_i = 0
index_j = 0
execute_on = timestep_end
[../]
[./stress_xy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xy
index_i = 0
index_j = 1
execute_on = timestep_end
[../]
[./stress_xz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xz
index_i = 0
index_j = 2
execute_on = timestep_end
[../]
[./stress_yx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yx
index_i = 1
index_j = 0
execute_on = timestep_end
[../]
[./stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
execute_on = timestep_end
[../]
[./stress_yz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yz
index_i = 1
index_j = 2
execute_on = timestep_end
[../]
[./stress_zx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zx
index_i = 2
index_j = 0
execute_on = timestep_end
[../]
[./stress_zy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zy
index_i = 2
index_j = 1
execute_on = timestep_end
[../]
[./stress_zz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zz
index_i = 2
index_j = 2
execute_on = timestep_end
[../]
[./total_strain_xx]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_xx
index_i = 0
index_j = 0
execute_on = timestep_end
[../]
[./total_strain_xy]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_xy
index_i = 0
index_j = 1
execute_on = timestep_end
[../]
[./total_strain_xz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_xz
index_i = 0
index_j = 2
execute_on = timestep_end
[../]
[./total_strain_yx]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_yx
index_i = 1
index_j = 0
execute_on = timestep_end
[../]
[./total_strain_yy]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_yy
index_i = 1
index_j = 1
execute_on = timestep_end
[../]
[./total_strain_yz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_yz
index_i = 1
index_j = 2
execute_on = timestep_end
[../]
[./total_strain_zx]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_zx
index_i = 2
index_j = 0
execute_on = timestep_end
[../]
[./total_strain_zy]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_zy
index_i = 2
index_j = 1
execute_on = timestep_end
[../]
[./total_strain_zz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = total_strain_zz
index_i = 2
index_j = 2
execute_on = timestep_end
[../]
[./perm_xx]
type = PorousFlowPropertyAux
property = permeability
variable = perm_xx
row = 0
column = 0
execute_on = timestep_end
[../]
[./perm_yy]
type = PorousFlowPropertyAux
property = permeability
variable = perm_yy
row = 1
column = 1
execute_on = timestep_end
[../]
[./perm_zz]
type = PorousFlowPropertyAux
property = permeability
variable = perm_zz
row = 2
column = 2
execute_on = timestep_end
[../]
[./mc_shear]
type = MaterialStdVectorAux
index = 0
property = mc_plastic_internal_parameter
variable = mc_shear
execute_on = timestep_end
[../]
[./mc_tensile]
type = MaterialStdVectorAux
index = 1
property = mc_plastic_internal_parameter
variable = mc_tensile
execute_on = timestep_end
[../]
[./wp_shear]
type = MaterialStdVectorAux
index = 0
property = wp_plastic_internal_parameter
variable = wp_shear
execute_on = timestep_end
[../]
[./wp_tensile]
type = MaterialStdVectorAux
index = 1
property = wp_plastic_internal_parameter
variable = wp_tensile
execute_on = timestep_end
[../]
[./mc_shear_f]
type = MaterialStdVectorAux
index = 6
property = mc_plastic_yield_function
variable = mc_shear_f
execute_on = timestep_end
[../]
[./mc_tensile_f]
type = MaterialStdVectorAux
index = 0
property = mc_plastic_yield_function
variable = mc_tensile_f
execute_on = timestep_end
[../]
[./wp_shear_f]
type = MaterialStdVectorAux
index = 0
property = wp_plastic_yield_function
variable = wp_shear_f
execute_on = timestep_end
[../]
[./wp_tensile_f]
type = MaterialStdVectorAux
index = 1
property = wp_plastic_yield_function
variable = wp_tensile_f
execute_on = timestep_end
[../]
[]
[BCs]
[./no_x]
type = DirichletBC
variable = disp_x
boundary = 'xmin xmax'
value = 0.0
[../]
[./no_y]
type = DirichletBC
variable = disp_y
boundary = 'ymin ymax'
value = 0.0
[../]
[./no_z]
type = DirichletBC
variable = disp_z
boundary = zmin
value = 0.0
[../]
[./no_wc_x]
type = DirichletBC
variable = wc_x
boundary = 'ymin ymax'
value = 0.0
[../]
[./no_wc_y]
type = DirichletBC
variable = wc_y
boundary = 'xmin xmax'
value = 0.0
[../]
[./fix_porepressure]
type = FunctionDirichletBC
variable = porepressure
boundary = 'ymin ymax xmax'
function = ini_pp
[../]
[./roof_porepressure]
type = PorousFlowPiecewiseLinearSink
variable = porepressure
pt_vals = '-1E3 1E3'
multipliers = '-1 1'
fluid_phase = 0
flux_function = roof_conductance
boundary = roof
[../]
[./roof]
type = StickyBC
variable = disp_z
min_value = -3.0
boundary = roof
[../]
[]
[Functions]
[./ini_pp]
type = ParsedFunction
vars = 'bulk p0 g rho0'
vals = '2E3 0.0 1E-5 1E3'
value = '-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk)'
[../]
[./ini_xx]
type = ParsedFunction
vars = 'bulk p0 g rho0 biot'
vals = '2E3 0.0 1E-5 1E3 0.7'
value = '0.8*(2500*10E-6*z+biot*(-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk)))'
[../]
[./ini_zz]
type = ParsedFunction
vars = 'bulk p0 g rho0 biot'
vals = '2E3 0.0 1E-5 1E3 0.7'
value = '2500*10E-6*z+biot*(-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk))'
[../]
[./excav_sideways]
type = ParsedFunction
vars = 'end_t ymin ymax minval maxval slope'
vals = '0.5 0 1000.0 1E-9 1 10'
# excavation face at ymin+(ymax-ymin)*min(t/end_t,1)
# slope is the distance over which the modulus reduces from maxval to minval
value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,if(y<ymin+(ymax-ymin)*min(t/end_t,1)+slope,minval+(maxval-minval)*(y-(ymin+(ymax-ymin)*min(t/end_t,1)))/slope,maxval))'
[../]
[./density_sideways]
type = ParsedFunction
vars = 'end_t ymin ymax minval maxval'
vals = '0.5 0 1000.0 0 2500'
value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,maxval)'
[../]
[./roof_conductance]
type = ParsedFunction
vars = 'end_t ymin ymax maxval minval'
vals = '0.5 0 1000.0 1E7 0'
value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),maxval,minval)'
[../]
[]
[UserObjects]
[./dictator]
type = PorousFlowDictator
porous_flow_vars = 'porepressure disp_x disp_y disp_z'
number_fluid_phases = 1
number_fluid_components = 1
[../]
[./pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1 # MPa^-1
[../]
[./mc_coh_strong_harden]
type = TensorMechanicsHardeningExponential
value_0 = 1.99 # MPa
value_residual = 2.01 # MPa
rate = 1.0
[../]
[./mc_fric]
type = TensorMechanicsHardeningConstant
value = 0.61 # 35deg
[../]
[./mc_dil]
type = TensorMechanicsHardeningConstant
value = 0.15 # 8deg
[../]
[./mc_tensile_str_strong_harden]
type = TensorMechanicsHardeningExponential
value_0 = 1.0 # MPa
value_residual = 1.0 # MPa
rate = 1.0
[../]
[./mc_compressive_str]
type = TensorMechanicsHardeningCubic
value_0 = 100 # Large!
value_residual = 100
internal_limit = 0.1
[../]
[./wp_coh_harden]
type = TensorMechanicsHardeningCubic
value_0 = 0.05
value_residual = 0.05
internal_limit = 10
[../]
[./wp_tan_fric]
type = TensorMechanicsHardeningConstant
value = 0.26 # 15deg
[../]
[./wp_tan_dil]
type = TensorMechanicsHardeningConstant
value = 0.18 # 10deg
[../]
[./wp_tensile_str_harden]
type = TensorMechanicsHardeningCubic
value_0 = 0.05
value_residual = 0.05
internal_limit = 10
[../]
[./wp_compressive_str_soften]
type = TensorMechanicsHardeningCubic
value_0 = 100
value_residual = 1
internal_limit = 1.0
[../]
[]
[Modules]
[./FluidProperties]
[./simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2E3
density0 = 1000
thermal_expansion = 0
viscosity = 3.5E-17
[../]
[../]
[]
[Materials]
[./temperature]
type = PorousFlowTemperature
[../]
[./eff_fluid_pressure]
type = PorousFlowEffectiveFluidPressure
[../]
[./vol_strain]
type = PorousFlowVolumetricStrain
[../]
[./ppss]
type = PorousFlow1PhaseP
porepressure = porepressure
capillary_pressure = pc
[../]
[./massfrac]
type = PorousFlowMassFraction
[../]
[./simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[../]
[./porosity_for_aux]
type = PorousFlowPorosity
at_nodes = false
fluid = true
mechanical = true
ensure_positive = true
porosity_zero = 0.02
solid_bulk = 5.3333E3
[../]
[./porosity_bulk]
type = PorousFlowPorosity
fluid = true
mechanical = true
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
ensure_positive = true
porosity_zero = 0.02
solid_bulk = 5.3333E3
[../]
[./porosity_excav]
type = PorousFlowPorosityConst
block = 1
porosity = 1.0
[../]
[./permeability_bulk]
type = PorousFlowPermeabilityKozenyCarman
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
poroperm_function = kozeny_carman_phi0
k0 = 1E-15
phi0 = 0.02
n = 2
m = 2
[../]
[./permeability_excav]
type = PorousFlowPermeabilityConst
block = 1
permeability = '0 0 0 0 0 0 0 0 0'
[../]
[./relperm]
type = PorousFlowRelativePermeabilityCorey
n = 4
s_res = 0.4
sum_s_res = 0.4
phase = 0
[../]
[./elasticity_tensor_0]
type = ComputeLayeredCosseratElasticityTensor
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
young = 8E3 # MPa
poisson = 0.25
layer_thickness = 1.0
joint_normal_stiffness = 1E9 # huge
joint_shear_stiffness = 1E3 # MPa
[../]
[./elasticity_tensor_1]
type = ComputeLayeredCosseratElasticityTensor
block = 1
young = 8E3 # MPa
poisson = 0.25
layer_thickness = 1.0
joint_normal_stiffness = 1E9 # huge
joint_shear_stiffness = 1E3 # MPa
elasticity_tensor_prefactor = excav_sideways
[../]
[./strain]
type = ComputeCosseratIncrementalSmallStrain
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
eigenstrain_name = ini_stress
initial_stress = 'ini_xx 0 0 0 ini_xx 0 0 0 ini_zz'
[../]
[./stress_0]
type = ComputeMultipleInelasticCosseratStress
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
inelastic_models = 'mc wp'
cycle_models = true
relative_tolerance = 2.0
absolute_tolerance = 1E6
max_iterations = 1
tangent_operator = nonlinear
perform_finite_strain_rotations = false
[../]
[./stress_1]
type = ComputeMultipleInelasticCosseratStress
block = 1
inelastic_models = ''
relative_tolerance = 2.0
absolute_tolerance = 1E6
max_iterations = 1
tangent_operator = nonlinear
perform_finite_strain_rotations = false
[../]
[./mc]
type = CappedMohrCoulombCosseratStressUpdate
warn_about_precision_loss = false
host_youngs_modulus = 8E3
host_poissons_ratio = 0.25
base_name = mc
tensile_strength = mc_tensile_str_strong_harden
compressive_strength = mc_compressive_str
cohesion = mc_coh_strong_harden
friction_angle = mc_fric
dilation_angle = mc_dil
max_NR_iterations = 100000
smoothing_tol = 0.1 # MPa # Must be linked to cohesion
yield_function_tol = 1E-9 # MPa. this is essentially the lowest possible without lots of precision loss
perfect_guess = true
min_step_size = 1.0
[../]
[./wp]
type = CappedWeakPlaneCosseratStressUpdate
warn_about_precision_loss = false
base_name = wp
cohesion = wp_coh_harden
tan_friction_angle = wp_tan_fric
tan_dilation_angle = wp_tan_dil
tensile_strength = wp_tensile_str_harden
compressive_strength = wp_compressive_str_soften
max_NR_iterations = 10000
tip_smoother = 0.05
smoothing_tol = 0.05 # MPa # Note, this must be tied to cohesion, otherwise get no possible return at cone apex
yield_function_tol = 1E-11 # MPa. this is essentially the lowest possible without lots of precision loss
perfect_guess = true
min_step_size = 1.0E-3
[../]
[./undrained_density_0]
type = GenericConstantMaterial
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
prop_names = density
prop_values = 2500
[../]
[./undrained_density_1]
type = GenericFunctionMaterial
block = 1
prop_names = density
prop_values = density_sideways
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[]
[]
[Postprocessors]
[./min_roof_disp]
type = NodalExtremeValue
boundary = roof
value_type = min
variable = disp_z
[../]
[./min_roof_pp]
type = NodalExtremeValue
boundary = roof
value_type = min
variable = porepressure
[../]
[./min_surface_disp]
type = NodalExtremeValue
boundary = zmax
value_type = min
variable = disp_z
[../]
[./min_surface_pp]
type = NodalExtremeValue
boundary = zmax
value_type = min
variable = porepressure
[../]
[./max_perm_zz]
type = ElementExtremeValue
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30'
variable = perm_zz
[../]
[]
[Executioner]
type = Transient
solve_type = 'NEWTON'
petsc_options = '-snes_converged_reason'
# best overall
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
# best if you don't have mumps:
#petsc_options_iname = '-pc_type -pc_asm_overlap -sub_pc_type -ksp_type -ksp_gmres_restart'
#petsc_options_value = ' asm 2 lu gmres 200'
# very basic:
#petsc_options_iname = '-pc_type -ksp_type -ksp_gmres_restart'
#petsc_options_value = ' bjacobi gmres 200'
line_search = bt
nl_abs_tol = 1e-3
nl_rel_tol = 1e-5
l_max_its = 200
nl_max_its = 30
start_time = 0.0
dt = 0.0025
end_time = 0.5
[]
[Outputs]
interval = 1
print_linear_residuals = true
exodus = true
csv = true
console = true
[]
modules/solid_mechanics/test/tests/generalized_plane_strain/weak_planestress.i
[GlobalParams]
order = FIRST
family = LAGRANGE
[]
[Mesh]
file = square.e
displacements = 'disp_x disp_y'
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./strain_zz]
[../]
[]
[AuxVariables]
[./temp]
[../]
[./saved_x]
[../]
[./saved_y]
[../]
[./saved_z]
[../]
[./stress_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_xy]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./aux_strain_zz]
order = CONSTANT
family = MONOMIAL
[../]
[]
[Postprocessors]
[./react_z]
type = MaterialTensorIntegralSM
tensor = stress
index = 2
[../]
[./min_strain_zz]
type = NodalExtremeValue
variable = strain_zz
value_type = min
[../]
[./max_strain_zz]
type = NodalExtremeValue
variable = strain_zz
value_type = max
[../]
[]
[SolidMechanics]
[./solid]
disp_x = disp_x
disp_y = disp_y
save_in_disp_x = saved_x
save_in_disp_y = saved_y
temp = temp
[../]
[]
[Kernels]
[./solid_z]
type = OutOfPlaneStress
variable = strain_zz
save_in = saved_z
disp_x = disp_x
disp_y = disp_y
temp = temp
[../]
[]
[AuxKernels]
[./tempfuncaux]
type = FunctionAux
variable = temp
function = tempfunc
use_displaced_mesh = false
[../]
[./stress_xx]
type = MaterialTensorAux
tensor = stress
variable = stress_xx
index = 0
[../]
[./stress_xy]
type = MaterialTensorAux
tensor = stress
variable = stress_xy
index = 3
[../]
[./stress_yy]
type = MaterialTensorAux
tensor = stress
variable = stress_yy
index = 1
[../]
[./stress_zz]
type = MaterialTensorAux
tensor = stress
variable = stress_zz
index = 2
[../]
[./strain_xx]
type = MaterialTensorAux
tensor = total_strain
variable = strain_xx
index = 0
[../]
[./strain_xy]
type = MaterialTensorAux
tensor = total_strain
variable = strain_xy
index = 3
[../]
[./strain_yy]
type = MaterialTensorAux
tensor = total_strain
variable = strain_yy
index = 1
[../]
[./strain_zz]
type = MaterialTensorAux
tensor = total_strain
variable = aux_strain_zz
index = 2
[../]
[]
[Functions]
[./pull]
type = PiecewiseLinear
x='0 1 100'
y='0 0.00 0.00'
[../]
[./tempfunc]
type = ParsedFunction
value = '(1-x)*t'
[../]
[]
[BCs]
[./bottomx]
type = DirichletBC
boundary = 1
variable = disp_x
value = 0.0
[../]
[./bottomy]
type = DirichletBC
boundary = 1
variable = disp_y
value = 0.0
[../]
[]
[Materials]
[./linelast]
type = Elastic
block = 1
disp_x = disp_x
disp_y = disp_y
poissons_ratio = 0.3
youngs_modulus = 1e6
thermal_expansion = 0.02
stress_free_temperature = 0.5
temp = temp
formulation = PlaneStrain
strain_zz = strain_zz
[../]
[]
[Executioner]
type = Transient
solve_type = PJFNK
line_search = none
# controls for linear iterations
l_max_its = 100
l_tol = 1e-4
# controls for nonlinear iterations
nl_max_its = 15
nl_rel_tol = 1e-10
nl_abs_tol = 1e-5
# time control
start_time = 0.0
dt = 1.0
dtmin = 1.0
end_time = 2.0
num_steps = 5000
[]
[Outputs]
exodus = true
[]
test/tests/postprocessors/nodal_extreme_value/nodal_extreme_pps_test.i
[Mesh]
type = FileMesh
file = trapezoid.e
uniform_refine = 1
# This test will not work in parallel with DistributedMesh enabled
# due to a bug in PeriodicBCs.
parallel_type = replicated
[]
[Functions]
[./tr_x]
type = ParsedFunction
value = -x*cos(pi/3)
[../]
[./tr_y]
type = ParsedFunction
value = x*sin(pi/3)
[../]
[./itr_x]
type = ParsedFunction
value = -x/cos(pi/3)
[../]
[./itr_y]
type = ParsedFunction
value = 0
[../]
[]
[Variables]
[./u]
order = FIRST
family = LAGRANGE
[../]
[]
[Kernels]
[./diff]
type = Diffusion
variable = u
[../]
[./forcing]
type = GaussContForcing
variable = u
x_center = 2
y_center = -1
x_spread = 0.25
y_spread = 0.5
[../]
[./dot]
type = TimeDerivative
variable = u
[../]
[]
[BCs]
# active = ' '
[./Periodic]
[./x]
primary = 1
secondary = 4
transform_func = 'tr_x tr_y'
inv_transform_func = 'itr_x itr_y'
[../]
[../]
[]
[Postprocessors]
[./max_nodal_pps]
type = NodalExtremeValue
variable = u
[../]
[./max_node_id]
type = NodalProxyMaxValue
variable = u
[../]
[./min_nodal_pps]
type = NodalExtremeValue
variable = u
value_type = min
[../]
[]
[Executioner]
type = Transient
dt = 0.5
num_steps = 6
[]
[Outputs]
execute_on = 'timestep_end'
exodus = true
[]
modules/solid_mechanics/test/tests/generalized_plane_strain/generalized_plane_strain_field_kinematic.i
# Tests kinematic formulation for the EqualValueBoundaryConstraint.
# The rest of the input file is identical to generalized_plane_strain.i
[GlobalParams]
order = FIRST
family = LAGRANGE
[]
[Mesh]
file = square.e
displacements = 'disp_x disp_y'
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./strain_zz]
[../]
[]
[AuxVariables]
[./temp]
[../]
[./saved_x]
[../]
[./saved_y]
[../]
[./saved_z]
[../]
[./stress_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_xy]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./aux_strain_zz]
order = CONSTANT
family = MONOMIAL
[../]
[]
[Postprocessors]
[./react_z]
type = MaterialTensorIntegralSM
tensor = stress
index = 2
[../]
[./min_strain_zz]
type = NodalExtremeValue
variable = strain_zz
value_type = min
[../]
[./max_strain_zz]
type = NodalExtremeValue
variable = strain_zz
value_type = max
[../]
[]
[SolidMechanics]
[./solid]
disp_x = disp_x
disp_y = disp_y
save_in_disp_x = saved_x
save_in_disp_y = saved_y
temp = temp
[../]
[]
[Kernels]
[./solid_z]
type = OutOfPlaneStress
variable = strain_zz
save_in = saved_z
disp_x = disp_x
disp_y = disp_y
temp = temp
[../]
[]
[Constraints]
[./szz]
type = EqualValueBoundaryConstraint
variable = strain_zz
master = '8'
slave = 10
penalty = 1e12
formulation = kinematic
[../]
[]
[AuxKernels]
[./tempfuncaux]
type = FunctionAux
variable = temp
function = tempfunc
use_displaced_mesh = false
[../]
[./stress_xx]
type = MaterialTensorAux
tensor = stress
variable = stress_xx
index = 0
[../]
[./stress_xy]
type = MaterialTensorAux
tensor = stress
variable = stress_xy
index = 3
[../]
[./stress_yy]
type = MaterialTensorAux
tensor = stress
variable = stress_yy
index = 1
[../]
[./stress_zz]
type = MaterialTensorAux
tensor = stress
variable = stress_zz
index = 2
[../]
[./strain_xx]
type = MaterialTensorAux
tensor = total_strain
variable = strain_xx
index = 0
[../]
[./strain_xy]
type = MaterialTensorAux
tensor = total_strain
variable = strain_xy
index = 3
[../]
[./strain_yy]
type = MaterialTensorAux
tensor = total_strain
variable = strain_yy
index = 1
[../]
[./strain_zz]
type = MaterialTensorAux
tensor = total_strain
variable = aux_strain_zz
index = 2
[../]
[]
[Functions]
[./pull]
type = PiecewiseLinear
x='0 1 100'
y='0 0.00 0.00'
[../]
[./tempfunc]
type = ParsedFunction
value = '(1-x)*t'
[../]
[]
[BCs]
[./bottomx]
type = DirichletBC
boundary = 1
variable = disp_x
value = 0.0
[../]
[./bottomy]
type = DirichletBC
boundary = 1
variable = disp_y
value = 0.0
[../]
[]
[Materials]
[./linelast]
type = Elastic
block = 1
disp_x = disp_x
disp_y = disp_y
poissons_ratio = 0.3
youngs_modulus = 1e6
thermal_expansion = 0.02
stress_free_temperature = 0.5
temp = temp
formulation = PlaneStrain
strain_zz = strain_zz
[../]
[]
[Executioner]
type = Transient
solve_type = PJFNK
line_search = none
# controls for linear iterations
l_max_its = 100
l_tol = 1e-4
# controls for nonlinear iterations
nl_max_its = 15
nl_rel_tol = 1e-10
nl_abs_tol = 1e-5
# time control
start_time = 0.0
dt = 1.0
dtmin = 1.0
end_time = 2.0
num_steps = 5000
[]
[Outputs]
exodus = true
[]
modules/tensor_mechanics/test/tests/plane_stress/weak_plane_stress_small.i
[GlobalParams]
order = FIRST
family = LAGRANGE
displacements = 'disp_x disp_y'
temperature = temp
out_of_plane_strain = strain_zz
[]
[Mesh]
file = square.e
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./strain_zz]
[../]
[]
[AuxVariables]
[./temp]
[../]
[./nl_strain_zz]
order = CONSTANT
family = MONOMIAL
[../]
[]
[Postprocessors]
[./react_z]
type = MaterialTensorIntegral
rank_two_tensor = stress
index_i = 2
index_j = 2
[../]
[./min_strain_zz]
type = NodalExtremeValue
variable = strain_zz
value_type = min
[../]
[./max_strain_zz]
type = NodalExtremeValue
variable = strain_zz
value_type = max
[../]
[]
[Modules/TensorMechanics/Master]
[plane_stress]
planar_formulation = WEAK_PLANE_STRESS
strain = SMALL
generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_xx strain_xy strain_yy'
eigenstrain_names = eigenstrain
[]
[]
[AuxKernels]
[./tempfuncaux]
type = FunctionAux
variable = temp
function = tempfunc
use_displaced_mesh = false
[../]
[./strain_zz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = nl_strain_zz
index_i = 2
index_j = 2
[../]
[]
[Functions]
[./pull]
type = PiecewiseLinear
x='0 1 100'
y='0 0.00 0.00'
[../]
[./tempfunc]
type = ParsedFunction
value = '(1 - x) * t'
[../]
[]
[BCs]
[./bottomx]
type = DirichletBC
boundary = 1
variable = disp_x
value = 0.0
[../]
[./bottomy]
type = DirichletBC
boundary = 1
variable = disp_y
value = 0.0
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
poissons_ratio = 0.3
youngs_modulus = 1e6
[../]
[./thermal_strain]
type = ComputeThermalExpansionEigenstrain
thermal_expansion_coeff = 0.02
stress_free_temperature = 0.5
eigenstrain_name = eigenstrain
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
[]
[Executioner]
type = Transient
solve_type = PJFNK
line_search = none
# controls for linear iterations
l_max_its = 100
l_tol = 1e-06
# controls for nonlinear iterations
nl_max_its = 15
nl_rel_tol = 1e-14
nl_abs_tol = 1e-12
# time control
start_time = 0.0
dt = 1.0
dtmin = 1.0
end_time = 2.0
[]
[Outputs]
exodus = true
[]
modules/tensor_mechanics/test/tests/plane_stress/weak_plane_stress_incremental.i
[GlobalParams]
order = FIRST
family = LAGRANGE
displacements = 'disp_x disp_y'
temperature = temp
out_of_plane_strain = strain_zz
[]
[Mesh]
file = square.e
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./strain_zz]
[../]
[]
[AuxVariables]
[./temp]
[../]
[./nl_strain_zz]
order = CONSTANT
family = MONOMIAL
[../]
[]
[Postprocessors]
[./react_z]
type = MaterialTensorIntegral
rank_two_tensor = stress
index_i = 2
index_j = 2
[../]
[./min_strain_zz]
type = NodalExtremeValue
variable = strain_zz
value_type = min
[../]
[./max_strain_zz]
type = NodalExtremeValue
variable = strain_zz
value_type = max
[../]
[]
[Modules/TensorMechanics/Master]
[plane_stress]
planar_formulation = WEAK_PLANE_STRESS
strain = SMALL
incremental = true
generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_xx strain_xy strain_yy'
eigenstrain_names = eigenstrain
[]
[]
[AuxKernels]
[./tempfuncaux]
type = FunctionAux
variable = temp
function = tempfunc
use_displaced_mesh = false
[../]
[./strain_zz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = nl_strain_zz
index_i = 2
index_j = 2
[../]
[]
[Functions]
[./pull]
type = PiecewiseLinear
x='0 1 100'
y='0 0.00 0.00'
[../]
[./tempfunc]
type = ParsedFunction
value = '(1 - x) * t'
[../]
[]
[BCs]
[./bottomx]
type = DirichletBC
boundary = 1
variable = disp_x
value = 0.0
[../]
[./bottomy]
type = DirichletBC
boundary = 1
variable = disp_y
value = 0.0
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
poissons_ratio = 0.3
youngs_modulus = 1e6
[../]
[./thermal_strain]
type = ComputeThermalExpansionEigenstrain
thermal_expansion_coeff = 0.02
stress_free_temperature = 0.5
eigenstrain_name = eigenstrain
[../]
[./stress]
type = ComputeStrainIncrementBasedStress
[../]
[]
[Executioner]
type = Transient
solve_type = PJFNK
line_search = none
# controls for linear iterations
l_max_its = 100
l_tol = 1e-06
# controls for nonlinear iterations
nl_max_its = 15
nl_rel_tol = 1e-14
nl_abs_tol = 1e-12
# time control
start_time = 0.0
dt = 1.0
dtmin = 1.0
end_time = 2.0
[]
[Outputs]
exodus = true
[]
modules/tensor_mechanics/examples/coal_mining/coarse.i
# Strata deformation and fracturing around a coal mine - 3D model
#
# A "half model" is used. The mine is 400m deep and
# just the roof is studied (-400<=z<=0). The mining panel
# sits between 0<=x<=150, and 0<=y<=1000, so this simulates
# a coal panel that is 300m wide and 1000m long. The outer boundaries
# are 1km from the excavation boundaries.
#
# Time is meaningless in this example
# as quasi-static solutions are sought at each timestep, but
# the number of timesteps controls the resolution of the
# process.
#
# The boundary conditions for this simulation are:
# - disp_x = 0 at x=0 and x=1150
# - disp_y = 0 at y=-1000 and y=1000
# - disp_z = 0 at z=-400, but there is a time-dependent
# Young's modulus that simulates excavation
# - wc_x = 0 at y=-1000 and y=1000
# - wc_y = 0 at x=0 and x=1150
# That is, rollers on the sides, free at top,
# and prescribed at bottom in the unexcavated portion.
#
# The small strain formulation is used.
#
# All stresses are measured in MPa. The initial stress is consistent with
# the weight force from density 2500 kg/m^3, ie, stress_zz = 0.025*z MPa
# where gravity = 10 m.s^-2 = 1E-5 MPa m^2/kg. The maximum and minimum
# principal horizontal stresses are assumed to be equal to 0.8*stress_zz.
#
# Material properties:
# Young's modulus = 8 GPa
# Poisson's ratio = 0.25
# Cosserat layer thickness = 1 m
# Cosserat-joint normal stiffness = large
# Cosserat-joint shear stiffness = 1 GPa
# MC cohesion = 3 MPa
# MC friction angle = 37 deg
# MC dilation angle = 8 deg
# MC tensile strength = 1 MPa
# MC compressive strength = 100 MPa
# WeakPlane cohesion = 0.1 MPa
# WeakPlane friction angle = 30 deg
# WeakPlane dilation angle = 10 deg
# WeakPlane tensile strength = 0.1 MPa
# WeakPlane compressive strength = 100 MPa softening to 1 MPa at strain = 1
#
[Mesh]
[file]
type = FileMeshGenerator
file = mesh/coarse.e
[]
[./xmin]
input = file
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
new_boundary = xmin
normal = '-1 0 0'
[../]
[./xmax]
input = xmin
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
new_boundary = xmax
normal = '1 0 0'
[../]
[./ymin]
input = xmax
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
new_boundary = ymin
normal = '0 -1 0'
[../]
[./ymax]
input = ymin
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
new_boundary = ymax
normal = '0 1 0'
[../]
[./zmax]
input = ymax
type = SideSetsAroundSubdomainGenerator
block = 16
new_boundary = zmax
normal = '0 0 1'
[../]
[./zmin]
input = zmax
type = SideSetsAroundSubdomainGenerator
block = 2
new_boundary = zmin
normal = '0 0 -1'
[../]
[./excav]
type = SubdomainBoundingBoxGenerator
input = zmin
block_id = 1
bottom_left = '0 0 -400'
top_right = '150 1000 -397'
[../]
[./roof]
type = SideSetsAroundSubdomainGenerator
block = 1
input = excav
new_boundary = roof
normal = '0 0 1'
[../]
[]
[GlobalParams]
perform_finite_strain_rotations = false
displacements = 'disp_x disp_y disp_z'
Cosserat_rotations = 'wc_x wc_y wc_z'
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./wc_x]
[../]
[./wc_y]
[../]
[]
[Kernels]
[./cx_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_x
component = 0
[../]
[./cy_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_y
component = 1
[../]
[./cz_elastic]
type = CosseratStressDivergenceTensors
use_displaced_mesh = false
variable = disp_z
component = 2
[../]
[./x_couple]
type = StressDivergenceTensors
use_displaced_mesh = false
variable = wc_x
displacements = 'wc_x wc_y wc_z'
component = 0
base_name = couple
[../]
[./y_couple]
type = StressDivergenceTensors
use_displaced_mesh = false
variable = wc_y
displacements = 'wc_x wc_y wc_z'
component = 1
base_name = couple
[../]
[./x_moment]
type = MomentBalancing
use_displaced_mesh = false
variable = wc_x
component = 0
[../]
[./y_moment]
type = MomentBalancing
use_displaced_mesh = false
variable = wc_y
component = 1
[../]
[./gravity]
type = Gravity
use_displaced_mesh = false
variable = disp_z
value = -10E-6 # remember this is in MPa
[../]
[]
[AuxVariables]
[./wc_z]
[../]
[./stress_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xz]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yz]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./mc_shear]
order = CONSTANT
family = MONOMIAL
[../]
[./mc_tensile]
order = CONSTANT
family = MONOMIAL
[../]
[./wp_shear]
order = CONSTANT
family = MONOMIAL
[../]
[./wp_tensile]
order = CONSTANT
family = MONOMIAL
[../]
[./wp_shear_f]
order = CONSTANT
family = MONOMIAL
[../]
[./wp_tensile_f]
order = CONSTANT
family = MONOMIAL
[../]
[./mc_shear_f]
order = CONSTANT
family = MONOMIAL
[../]
[./mc_tensile_f]
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_yx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yx
index_i = 1
index_j = 0
[../]
[./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_zx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zx
index_i = 2
index_j = 0
[../]
[./stress_zy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zy
index_i = 2
index_j = 1
[../]
[./stress_zz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zz
index_i = 2
index_j = 2
[../]
[./mc_shear]
type = MaterialStdVectorAux
index = 0
property = mc_plastic_internal_parameter
variable = mc_shear
[../]
[./mc_tensile]
type = MaterialStdVectorAux
index = 1
property = mc_plastic_internal_parameter
variable = mc_tensile
[../]
[./wp_shear]
type = MaterialStdVectorAux
index = 0
property = wp_plastic_internal_parameter
variable = wp_shear
[../]
[./wp_tensile]
type = MaterialStdVectorAux
index = 1
property = wp_plastic_internal_parameter
variable = wp_tensile
[../]
[./mc_shear_f]
type = MaterialStdVectorAux
index = 6
property = mc_plastic_yield_function
variable = mc_shear_f
[../]
[./mc_tensile_f]
type = MaterialStdVectorAux
index = 0
property = mc_plastic_yield_function
variable = mc_tensile_f
[../]
[./wp_shear_f]
type = MaterialStdVectorAux
index = 0
property = wp_plastic_yield_function
variable = wp_shear_f
[../]
[./wp_tensile_f]
type = MaterialStdVectorAux
index = 1
property = wp_plastic_yield_function
variable = wp_tensile_f
[../]
[]
[BCs]
[./no_x]
type = DirichletBC
variable = disp_x
boundary = 'xmin xmax'
value = 0.0
[../]
[./no_y]
type = DirichletBC
variable = disp_y
boundary = 'ymin ymax'
value = 0.0
[../]
[./no_z]
type = DirichletBC
variable = disp_z
boundary = zmin
value = 0.0
[../]
[./no_wc_x]
type = DirichletBC
variable = wc_x
boundary = 'ymin ymax'
value = 0.0
[../]
[./no_wc_y]
type = DirichletBC
variable = wc_y
boundary = 'xmin xmax'
value = 0.0
[../]
[./roof]
type = StickyBC
variable = disp_z
min_value = -3.0
boundary = roof
[../]
[]
[Functions]
[./ini_xx]
type = ParsedFunction
value = '0.8*2500*10E-6*z'
[../]
[./ini_zz]
type = ParsedFunction
value = '2500*10E-6*z'
[../]
[./excav_sideways]
type = ParsedFunction
vars = 'end_t ymin ymax minval maxval slope'
vals = '17.0 0 1000.0 1E-9 1 60'
# excavation face at ymin+(ymax-ymin)*min(t/end_t,1)
# slope is the distance over which the modulus reduces from maxval to minval
value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,if(y<ymin+(ymax-ymin)*min(t/end_t,1)+slope,minval+(maxval-minval)*(y-(ymin+(ymax-ymin)*min(t/end_t,1)))/slope,maxval))'
[../]
[./density_sideways]
type = ParsedFunction
vars = 'end_t ymin ymax minval maxval'
vals = '17.0 0 1000.0 0 2500'
value = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,maxval)'
[../]
[]
[UserObjects]
[./mc_coh_strong_harden]
type = TensorMechanicsHardeningExponential
value_0 = 2.99 # MPa
value_residual = 3.01 # MPa
rate = 1.0
[../]
[./mc_fric]
type = TensorMechanicsHardeningConstant
value = 0.65 # 37deg
[../]
[./mc_dil]
type = TensorMechanicsHardeningConstant
value = 0.15 # 8deg
[../]
[./mc_tensile_str_strong_harden]
type = TensorMechanicsHardeningExponential
value_0 = 1.0 # MPa
value_residual = 1.0 # MPa
rate = 1.0
[../]
[./mc_compressive_str]
type = TensorMechanicsHardeningCubic
value_0 = 100 # Large!
value_residual = 100
internal_limit = 0.1
[../]
[./wp_coh_harden]
type = TensorMechanicsHardeningCubic
value_0 = 0.1
value_residual = 0.1
internal_limit = 10
[../]
[./wp_tan_fric]
type = TensorMechanicsHardeningConstant
value = 0.36 # 20deg
[../]
[./wp_tan_dil]
type = TensorMechanicsHardeningConstant
value = 0.18 # 10deg
[../]
[./wp_tensile_str_harden]
type = TensorMechanicsHardeningCubic
value_0 = 0.1
value_residual = 0.1
internal_limit = 10
[../]
[./wp_compressive_str_soften]
type = TensorMechanicsHardeningCubic
value_0 = 100
value_residual = 1
internal_limit = 1.0
[../]
[]
[Materials]
[./elasticity_tensor_0]
type = ComputeLayeredCosseratElasticityTensor
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
young = 8E3 # MPa
poisson = 0.25
layer_thickness = 1.0
joint_normal_stiffness = 1E9 # huge
joint_shear_stiffness = 1E3 # MPa
[../]
[./elasticity_tensor_1]
type = ComputeLayeredCosseratElasticityTensor
block = 1
young = 8E3 # MPa
poisson = 0.25
layer_thickness = 1.0
joint_normal_stiffness = 1E9 # huge
joint_shear_stiffness = 1E3 # MPa
elasticity_tensor_prefactor = excav_sideways
[../]
[./strain]
type = ComputeCosseratIncrementalSmallStrain
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
eigenstrain_name = ini_stress
initial_stress = 'ini_xx 0 0 0 ini_xx 0 0 0 ini_zz'
[../]
[./stress_0]
type = ComputeMultipleInelasticCosseratStress
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
inelastic_models = 'mc wp'
cycle_models = true
relative_tolerance = 2.0
absolute_tolerance = 1E6
max_iterations = 1
tangent_operator = nonlinear
perform_finite_strain_rotations = false
[../]
[./stress_1]
type = ComputeMultipleInelasticCosseratStress
block = 1
inelastic_models = ''
relative_tolerance = 2.0
absolute_tolerance = 1E6
max_iterations = 1
tangent_operator = nonlinear
perform_finite_strain_rotations = false
[../]
[./mc]
type = CappedMohrCoulombCosseratStressUpdate
warn_about_precision_loss = false
host_youngs_modulus = 8E3
host_poissons_ratio = 0.25
base_name = mc
tensile_strength = mc_tensile_str_strong_harden
compressive_strength = mc_compressive_str
cohesion = mc_coh_strong_harden
friction_angle = mc_fric
dilation_angle = mc_dil
max_NR_iterations = 100000
smoothing_tol = 0.1 # MPa # Must be linked to cohesion
yield_function_tol = 1E-9 # MPa. this is essentially the lowest possible without lots of precision loss
perfect_guess = true
min_step_size = 1.0
[../]
[./wp]
type = CappedWeakPlaneCosseratStressUpdate
warn_about_precision_loss = false
base_name = wp
cohesion = wp_coh_harden
tan_friction_angle = wp_tan_fric
tan_dilation_angle = wp_tan_dil
tensile_strength = wp_tensile_str_harden
compressive_strength = wp_compressive_str_soften
max_NR_iterations = 10000
tip_smoother = 0.1
smoothing_tol = 0.1 # MPa # Note, this must be tied to cohesion, otherwise get no possible return at cone apex
yield_function_tol = 1E-11 # MPa. this is essentially the lowest possible without lots of precision loss
perfect_guess = true
min_step_size = 1.0E-3
[../]
[./density_0]
type = GenericConstantMaterial
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
prop_names = density
prop_values = 2500
[../]
[./density_1]
type = GenericFunctionMaterial
block = 1
prop_names = density
prop_values = density_sideways
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[]
[]
[Postprocessors]
[./min_roof_disp]
type = NodalExtremeValue
boundary = roof
value_type = min
variable = disp_z
[../]
[./min_surface_disp]
type = NodalExtremeValue
boundary = zmax
value_type = min
variable = disp_z
[../]
[]
[Executioner]
type = Transient
solve_type = 'NEWTON'
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-pc_type -ksp_type -ksp_gmres_restart'
petsc_options_value = ' bjacobi gmres 200'
line_search = bt
nl_abs_tol = 1e-3
nl_rel_tol = 1e-5
l_max_its = 30
nl_max_its = 1000
start_time = 0.0
dt = 0.5 # this gives min(disp_z)=-4.3, use dt=0.0625 if you want to restrict disp_z>=-3.2
end_time = 17.0
[]
[Outputs]
interval = 1
print_linear_residuals = false
exodus = true
csv = true
console = true
[]
modules/solid_mechanics/test/tests/generalized_plane_strain/out_of_plane_pressure.i
# Tests the use of OutOfPlanePressure with generalized plane strain.
[GlobalParams]
order = FIRST
family = LAGRANGE
[]
[Mesh]
file = square.e
displacements = 'disp_x disp_y'
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./strain_zz]
[../]
[]
[AuxVariables]
[./saved_x]
[../]
[./saved_y]
[../]
[./saved_z]
[../]
[./stress_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_xy]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./aux_strain_zz]
order = CONSTANT
family = MONOMIAL
[../]
[]
[Postprocessors]
[./react_z]
type = MaterialTensorIntegralSM
tensor = stress
index = 2
[../]
[./min_strain_zz]
type = NodalExtremeValue
variable = strain_zz
value_type = min
[../]
[./max_strain_zz]
type = NodalExtremeValue
variable = strain_zz
value_type = max
[../]
[]
[SolidMechanics]
[./solid]
disp_x = disp_x
disp_y = disp_y
save_in_disp_x = saved_x
save_in_disp_y = saved_y
[../]
[]
[Kernels]
[./solid_z]
type = OutOfPlaneStress
variable = strain_zz
save_in = saved_z
disp_x = disp_x
disp_y = disp_y
[../]
[./pressure_z]
type = OutOfPlanePressure
variable = strain_zz
function = pressure_function
save_in = saved_z
factor = 1e5
[../]
[]
[Constraints]
[./szz]
type = EqualValueBoundaryConstraint
variable = strain_zz
master = '8'
slave = 10
penalty = 1e12
formulation = kinematic
[../]
[]
[AuxKernels]
[./stress_xx]
type = MaterialTensorAux
tensor = stress
variable = stress_xx
index = 0
[../]
[./stress_xy]
type = MaterialTensorAux
tensor = stress
variable = stress_xy
index = 3
[../]
[./stress_yy]
type = MaterialTensorAux
tensor = stress
variable = stress_yy
index = 1
[../]
[./stress_zz]
type = MaterialTensorAux
tensor = stress
variable = stress_zz
index = 2
[../]
[./strain_xx]
type = MaterialTensorAux
tensor = total_strain
variable = strain_xx
index = 0
[../]
[./strain_xy]
type = MaterialTensorAux
tensor = total_strain
variable = strain_xy
index = 3
[../]
[./strain_yy]
type = MaterialTensorAux
tensor = total_strain
variable = strain_yy
index = 1
[../]
[./strain_zz]
type = MaterialTensorAux
tensor = total_strain
variable = aux_strain_zz
index = 2
[../]
[]
[Functions]
[./pressure_function]
type = PiecewiseLinear
x='0 2'
y='0 1'
[../]
[]
[BCs]
[./leftx]
type = DirichletBC
boundary = 4
variable = disp_x
value = 0.0
[../]
[./bottomy]
type = DirichletBC
boundary = 1
variable = disp_y
value = 0.0
[../]
[]
[Materials]
[./linelast]
type = Elastic
block = 1
disp_x = disp_x
disp_y = disp_y
poissons_ratio = 0.3
youngs_modulus = 1e6
formulation = PlaneStrain
strain_zz = strain_zz
[../]
[]
[Executioner]
type = Transient
solve_type = PJFNK
line_search = none
# controls for linear iterations
l_max_its = 100
l_tol = 1e-4
# controls for nonlinear iterations
nl_max_its = 15
nl_rel_tol = 1e-14
nl_abs_tol = 1e-11
# time control
start_time = 0.0
dt = 1.0
dtmin = 1.0
end_time = 2.0
num_steps = 5000
[]
[Outputs]
exodus = true
[]