- boundaryThe list of boundary IDs from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundary IDs from the mesh where this object applies
- variableThe name of the variable that this residual object operates on
C++ Type:NonlinearVariableName
Controllable:No
Description:The name of the variable that this residual object operates on
StickyBC
Imposes the boundary condition if exceeds the bounds provided
Description
This nodal boundary condition imposes the condition if exceeds either or Hence the name "sticky": as soon as exceeds the bounds it is fixed at subsequent timesteps.
commentnote
StickyBC should be used with care. It only approximates a Constraint imposed by a wall, for instance, for: the boundary nodes may penetrate the wall in a single timestep before they are then stuck (inside the wall); boundary nodes can never move away from the wall again once they are stuck. However, it offers the advantage over a full-blown Constraint that it does not adversely affect numerical convergence.
Input Parameters
- diag_save_inThe name of auxiliary variables to save this BC's diagonal jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Controllable:No
Description:The name of auxiliary variables to save this BC's diagonal jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
- displacementsThe displacements
C++ Type:std::vector<VariableName>
Controllable:No
Description:The displacements
- max_value1.79769e+308If the old variable value >= max_value, the variable is fixed at its old value
Default:1.79769e+308
C++ Type:double
Controllable:No
Description:If the old variable value >= max_value, the variable is fixed at its old value
- min_value-1.79769e+308If the old variable value <= min_value, the variable is fixed at its old value
Default:-1.79769e+308
C++ Type:double
Controllable:No
Description:If the old variable value <= min_value, the variable is fixed at its old value
- save_inThe name of auxiliary variables to save this BC's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Controllable:No
Description:The name of auxiliary variables to save this BC's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
Optional Parameters
- absolute_value_vector_tagsThe tags for the vectors this residual object should fill with the absolute value of the residual contribution
C++ Type:std::vector<TagName>
Controllable:No
Description:The tags for the vectors this residual object should fill with the absolute value of the residual contribution
- extra_matrix_tagsThe extra tags for the matrices this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the matrices this Kernel should fill
- extra_vector_tagsThe extra tags for the vectors this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the vectors this Kernel should fill
- matrix_tagssystem timeThe tag for the matrices this Kernel should fill
Default:system time
C++ Type:MultiMooseEnum
Options:nontime, system, time
Controllable:No
Description:The tag for the matrices this Kernel should fill
- vector_tagsresidualThe tag for the vectors this Kernel should fill
Default:residual
C++ Type:MultiMooseEnum
Options:nontime, time, residual
Controllable:No
Description:The tag for the vectors this Kernel should fill
Tagging Parameters
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
- implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
- seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
- use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Advanced Parameters
Input Files
- (modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)
- (modules/solid_mechanics/examples/coal_mining/cosserat_mc_wp_sticky_longitudinal.i)
- (modules/solid_mechanics/examples/coal_mining/cosserat_mc_wp_sticky.i)
- (modules/solid_mechanics/examples/coal_mining/fine.i)
- (modules/solid_mechanics/test/tests/stickyBC/push_up.i)
- (modules/solid_mechanics/test/tests/stickyBC/push_down.i)
- (modules/solid_mechanics/test/tests/stickyBC/except1.i)
- (modules/porous_flow/examples/coal_mining/fine_with_fluid.i)
- (modules/solid_mechanics/examples/coal_mining/coarse.i)
(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
primary_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
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
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
symbol_names = 'bulk p0 g rho0'
symbol_values = '2E3 0.0 1E-5 1E3'
expression = '-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk)'
[]
[ini_xx]
type = ParsedFunction
symbol_names = 'bulk p0 g rho0 biot'
symbol_values = '2E3 0.0 1E-5 1E3 0.7'
expression = '0.8*(2500*10E-6*z+biot*(-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk)))'
[]
[ini_zz]
type = ParsedFunction
symbol_names = 'bulk p0 g rho0 biot'
symbol_values = '2E3 0.0 1E-5 1E3 0.7'
expression = '2500*10E-6*z+biot*(-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk))'
[]
[excav_sideways]
type = ParsedFunction
symbol_names = 'end_t ymin ymax minval maxval slope'
symbol_values = '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
expression = '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
symbol_names = 'end_t ymin ymax minval maxval'
symbol_values = '0.5 0 1000.0 0 2500'
expression = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,maxval)'
[]
[roof_conductance]
type = ParsedFunction
symbol_names = 'end_t ymin ymax maxval minval'
symbol_values = '0.5 0 1000.0 1E7 0'
expression = '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
[]
[]
[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]
time_step_interval = 1
print_linear_residuals = true
exodus = true
csv = true
console = true
[]
(modules/solid_mechanics/examples/coal_mining/cosserat_mc_wp_sticky_longitudinal.i)
# Strata deformation and fracturing around a coal mine
#
# A 2D geometry is used that simulates a longitudinal section of
# the coal mine. The model is actually 3D, but the "x"
# dimension is only 10m long, meshed with 1 element, and
# there is no "x" displacement. The mine is 400m deep
# and just the roof is studied (0<=z<=400). The model sits
# between -300<=y<=1800. The excavation sits in 0<=y<=1500. The
# excavation height is 3m (ie, the excavation lies within
# 0<=z<=3).
#
# 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 elastic simulation are:
# - disp_x = 0 everywhere
# - disp_y = 0 at y=-300 and y=1800
# - disp_z = 0 at z=0, but there is a time-dependent
# Young's modulus that simulates excavation
# - wc_x = 0 at y=300 and y=1800.
# 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*(300-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
#
[Mesh]
[generated_mesh]
type = GeneratedMeshGenerator
dim = 3
nx = 1
xmin = -5
xmax = 5
nz = 40
zmin = 0
zmax = 400
bias_z = 1.1
ny = 140 # 15m elements
ymin = -300
ymax = 1800
[]
[left]
type = SideSetsAroundSubdomainGenerator
block = 0
new_boundary = 11
normal = '0 -1 0'
input = generated_mesh
[]
[right]
type = SideSetsAroundSubdomainGenerator
block = 0
new_boundary = 12
normal = '0 1 0'
input = left
[]
[front]
type = SideSetsAroundSubdomainGenerator
block = 0
new_boundary = 13
normal = '-1 0 0'
input = right
[]
[back]
type = SideSetsAroundSubdomainGenerator
block = 0
new_boundary = 14
normal = '1 0 0'
input = front
[]
[top]
type = SideSetsAroundSubdomainGenerator
block = 0
new_boundary = 15
normal = '0 0 1'
input = back
[]
[bottom]
type = SideSetsAroundSubdomainGenerator
block = 0
new_boundary = 16
normal = '0 0 -1'
input = top
[]
[excav]
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '-5 0 0'
top_right = '5 1500 3'
input = bottom
[]
[roof]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 18
normal = '0 0 1'
input = excav
[]
[]
[GlobalParams]
perform_finite_strain_rotations = false
displacements = 'disp_x disp_y disp_z'
Cosserat_rotations = 'wc_x wc_y wc_z'
[]
[Variables]
[./disp_y]
[../]
[./disp_z]
[../]
[./wc_x]
[../]
[]
[Kernels]
[./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
[../]
[./x_moment]
type = MomentBalancing
use_displaced_mesh = false
variable = wc_x
component = 0
[../]
[./gravity]
type = Gravity
use_displaced_mesh = false
variable = disp_z
value = -10E-6 # remember this is in MPa
[../]
[]
[AuxVariables]
[./disp_x]
[../]
[./wc_y]
[../]
[./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_y]
type = DirichletBC
variable = disp_y
boundary = '11 12'
value = 0.0
[../]
[./no_z]
type = DirichletBC
variable = disp_z
boundary = '16'
value = 0.0
[../]
[./no_wc_x]
type = DirichletBC
variable = wc_x
boundary = '11 12'
value = 0.0
[../]
[./roof]
type = StickyBC
variable = disp_z
min_value = -3.0
boundary = '18'
[../]
[]
[Functions]
[./ini_xx]
type = ParsedFunction
expression = '-0.8*2500*10E-6*(400-z)'
[../]
[./ini_zz]
type = ParsedFunction
expression = '-2500*10E-6*(400-z)'
[../]
[./excav_sideways]
type = ParsedFunction
symbol_names = 'end_t ymin ymax minval maxval slope'
symbol_values = '1.0 0 1500.0 1E-9 1 15'
# excavation face at ymin+(ymax-ymin)*min(t/end_t,1)
# slope is the distance over which the modulus reduces from maxval to minval
expression = '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
symbol_names = 'end_t ymin ymax minval maxval'
symbol_values = '1.0 0 1500.0 0 2500'
expression = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,maxval)'
[../]
[]
[UserObjects]
[./mc_coh_strong_harden]
type = SolidMechanicsHardeningExponential
value_0 = 2.99 # MPa
value_residual = 3.01 # MPa
rate = 1.0
[../]
[./mc_fric]
type = SolidMechanicsHardeningConstant
value = 0.65 # 37deg
[../]
[./mc_dil]
type = SolidMechanicsHardeningConstant
value = 0.15 # 8deg
[../]
[./mc_tensile_str_strong_harden]
type = SolidMechanicsHardeningExponential
value_0 = 1.0 # MPa
value_residual = 1.0 # MPa
rate = 1.0
[../]
[./mc_compressive_str]
type = SolidMechanicsHardeningCubic
value_0 = 100 # Large!
value_residual = 100
internal_limit = 0.1
[../]
[./wp_coh_harden]
type = SolidMechanicsHardeningCubic
value_0 = 0.1
value_residual = 0.1
internal_limit = 10
[../]
[./wp_tan_fric]
type = SolidMechanicsHardeningConstant
value = 0.36 # 20deg
[../]
[./wp_tan_dil]
type = SolidMechanicsHardeningConstant
value = 0.18 # 10deg
[../]
[./wp_tensile_str_harden]
type = SolidMechanicsHardeningCubic
value_0 = 0.1
value_residual = 0.1
internal_limit = 10
[../]
[./wp_compressive_str_soften]
type = SolidMechanicsHardeningCubic
value_0 = 100
value_residual = 1.0
internal_limit = 1.0
[../]
[]
[Materials]
[./elasticity_tensor_0]
type = ComputeLayeredCosseratElasticityTensor
block = 0
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 = 0
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]
# this is needed so as to correctly apply the initial stress
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 = 0
prop_names = density
prop_values = 2500
[../]
[./density_1]
type = GenericFunctionMaterial
block = 1
prop_names = density
prop_values = density_sideways
[../]
[]
[Postprocessors]
[./subs]
type = PointValue
point = '0 0 400'
variable = disp_z
use_displaced_mesh = false
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = 'NEWTON'
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-pc_type -pc_asm_overlap -sub_pc_type -ksp_type -ksp_gmres_restart'
petsc_options_value = ' asm 2 lu gmres 200'
line_search = bt
nl_abs_tol = 1e-3
nl_rel_tol = 1e-5
l_max_its = 30
nl_max_its = 100
start_time = 0.0
dt = 0.01 # 1 element per step
end_time = 1.0
[]
[Outputs]
file_base = cosserat_mc_wp_sticky_longitudinal
time_step_interval = 1
print_linear_residuals = false
exodus = true
csv = true
console = true
#[./console]
# type = Console
# output_linear = false
#[../]
[]
(modules/solid_mechanics/examples/coal_mining/cosserat_mc_wp_sticky.i)
# Strata deformation and fracturing around a coal mine
#
# A 2D geometry is used that simulates a transverse section of
# the coal mine. The model is actually 3D, but the "x"
# dimension is only 10m long, meshed with 1 element, and
# there is no "x" displacement. The mine is 400m deep
# and just the roof is studied (0<=z<=400). The model sits
# between 0<=y<=450. The excavation sits in 0<=y<=150. This
# is a "half model": the boundary conditions are such that
# the model simulates an excavation sitting in -150<=y<=150
# inside a model of the region -450<=y<=450. The
# excavation height is 3m (ie, the excavation lies within
# 0<=z<=3).
#
# 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 elastic simulation are:
# - disp_x = 0 everywhere
# - disp_y = 0 at y=0 and y=450
# - disp_z = 0 at z=0, but there is a time-dependent
# Young's modulus that simulates excavation
# - wc_x = 0 at y=0 and y=450.
# 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*(300-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, varying down to 1 MPa when tensile strain = 1
# 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]
[generated_mesh]
type = GeneratedMeshGenerator
dim = 3
nx = 1
xmin = -5
xmax = 5
nz = 40
zmin = 0
zmax = 403.003
bias_z = 1.1
ny = 30 # make this a multiple of 3, so y=150 is at a node
ymin = 0
ymax = 450
[]
[left]
type = SideSetsAroundSubdomainGenerator
block = 0
new_boundary = 11
normal = '0 -1 0'
input = generated_mesh
[]
[right]
type = SideSetsAroundSubdomainGenerator
block = 0
new_boundary = 12
normal = '0 1 0'
input = left
[]
[front]
type = SideSetsAroundSubdomainGenerator
block = 0
new_boundary = 13
normal = '-1 0 0'
input = right
[]
[back]
type = SideSetsAroundSubdomainGenerator
block = 0
new_boundary = 14
normal = '1 0 0'
input = front
[]
[top]
type = SideSetsAroundSubdomainGenerator
block = 0
new_boundary = 15
normal = '0 0 1'
input = back
[]
[bottom]
type = SideSetsAroundSubdomainGenerator
block = 0
new_boundary = 16
normal = '0 0 -1'
input = top
[]
[excav]
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '-5 0 0'
top_right = '5 150 3'
input = bottom
[]
[roof]
type = SideSetsAroundSubdomainGenerator
block = 1
new_boundary = 18
normal = '0 0 1'
input = excav
[]
[]
[GlobalParams]
perform_finite_strain_rotations = false
displacements = 'disp_x disp_y disp_z'
Cosserat_rotations = 'wc_x wc_y wc_z'
[]
[Variables]
[./disp_y]
[../]
[./disp_z]
[../]
[./wc_x]
[../]
[]
[Kernels]
[./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
[../]
[./x_moment]
type = MomentBalancing
use_displaced_mesh = false
variable = wc_x
component = 0
[../]
[./gravity]
type = Gravity
use_displaced_mesh = false
variable = disp_z
value = -10E-6 # remember this is in MPa
[../]
[]
[AuxVariables]
[./disp_x]
[../]
[./wc_y]
[../]
[./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_y]
type = DirichletBC
variable = disp_y
boundary = '11 12'
value = 0.0
[../]
[./no_z]
type = DirichletBC
variable = disp_z
boundary = '16'
value = 0.0
[../]
[./no_wc_x]
type = DirichletBC
variable = wc_x
boundary = '11 12'
value = 0.0
[../]
[./roof]
type = StickyBC
variable = disp_z
min_value = -3.0
boundary = '18'
[../]
[]
[Functions]
[./ini_xx]
type = ParsedFunction
expression = '-0.8*2500*10E-6*(403.003-z)'
[../]
[./ini_zz]
type = ParsedFunction
expression = '-2500*10E-6*(403.003-z)'
[../]
[./excav_sideways]
type = ParsedFunction
symbol_names = 'end_t ymin ymax minval maxval slope'
symbol_values = '1.0 0 150.0 1E-9 1 15'
# excavation face at ymin+(ymax-ymin)*min(t/end_t,1)
# slope is the distance over which the modulus reduces from maxval to minval
expression = '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
symbol_names = 'end_t ymin ymax minval maxval'
symbol_values = '1.0 0 150.0 0 2500'
expression = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,maxval)'
[../]
[]
[UserObjects]
[./mc_coh_strong_harden]
type = SolidMechanicsHardeningExponential
value_0 = 2.99 # MPa
value_residual = 3.01 # MPa
rate = 1.0
[../]
[./mc_fric]
type = SolidMechanicsHardeningConstant
value = 0.65 # 37deg
[../]
[./mc_dil]
type = SolidMechanicsHardeningConstant
value = 0.15 # 8deg
[../]
[./mc_tensile_str_strong_harden]
type = SolidMechanicsHardeningExponential
value_0 = 1.0 # MPa
value_residual = 1.0 # MPa
rate = 1.0
[../]
[./mc_compressive_str]
type = SolidMechanicsHardeningCubic
value_0 = 100 # Large!
value_residual = 100
internal_limit = 0.1
[../]
[./wp_coh_harden]
type = SolidMechanicsHardeningCubic
value_0 = 0.1
value_residual = 0.1
internal_limit = 10
[../]
[./wp_tan_fric]
type = SolidMechanicsHardeningConstant
value = 0.36 # 20deg
[../]
[./wp_tan_dil]
type = SolidMechanicsHardeningConstant
value = 0.18 # 10deg
[../]
[./wp_tensile_str_harden]
type = SolidMechanicsHardeningCubic
value_0 = 0.1
value_residual = 0.1
internal_limit = 10
[../]
[./wp_compressive_str_soften]
type = SolidMechanicsHardeningCubic
value_0 = 100
value_residual = 1
internal_limit = 1.0
[../]
[]
[Materials]
[./elasticity_tensor_0]
type = ComputeLayeredCosseratElasticityTensor
block = 0
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]
# this is needed so as to correctly apply the initial stress
type = ComputeMultipleInelasticCosseratStress
block = 0
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 = 0
prop_names = density
prop_values = 2500
[../]
[./density_1]
type = GenericFunctionMaterial
block = 1
prop_names = density
prop_values = density_sideways
[../]
[]
[Postprocessors]
[./subs_max]
type = PointValue
point = '0 0 403.003'
variable = disp_z
use_displaced_mesh = false
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = 'NEWTON'
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-pc_type -pc_asm_overlap -sub_pc_type -ksp_type -ksp_gmres_restart'
petsc_options_value = ' asm 2 lu gmres 200'
line_search = bt
nl_abs_tol = 1e-8
nl_rel_tol = 1e-8
l_max_its = 30
nl_max_its = 1000
start_time = 0.0
dt = 0.01
end_time = 1.0
[]
[Outputs]
file_base = cosserat_mc_wp_sticky
time_step_interval = 1
print_linear_residuals = false
exodus = true
csv = true
console = true
[]
(modules/solid_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
expression = '0.8*2500*10E-6*z'
[../]
[./ini_zz]
type = ParsedFunction
expression = '2500*10E-6*z'
[../]
[./excav_sideways]
type = ParsedFunction
symbol_names = 'end_t ymin ymax minval maxval slope'
symbol_values = '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
expression = '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
symbol_names = 'end_t ymin ymax minval maxval'
symbol_values = '100.0 0 1000.0 0 2500'
expression = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,maxval)'
[../]
[]
[UserObjects]
[./mc_coh_strong_harden]
type = SolidMechanicsHardeningExponential
value_0 = 2.99 # MPa
value_residual = 3.01 # MPa
rate = 1.0
[../]
[./mc_fric]
type = SolidMechanicsHardeningConstant
value = 0.65 # 37deg
[../]
[./mc_dil]
type = SolidMechanicsHardeningConstant
value = 0.15 # 8deg
[../]
[./mc_tensile_str_strong_harden]
type = SolidMechanicsHardeningExponential
value_0 = 1.0 # MPa
value_residual = 1.0 # MPa
rate = 1.0
[../]
[./mc_compressive_str]
type = SolidMechanicsHardeningCubic
value_0 = 100 # Large!
value_residual = 100
internal_limit = 0.1
[../]
[./wp_coh_harden]
type = SolidMechanicsHardeningCubic
value_0 = 0.1
value_residual = 0.1
internal_limit = 10
[../]
[./wp_tan_fric]
type = SolidMechanicsHardeningConstant
value = 0.36 # 20deg
[../]
[./wp_tan_dil]
type = SolidMechanicsHardeningConstant
value = 0.18 # 10deg
[../]
[./wp_tensile_str_harden]
type = SolidMechanicsHardeningCubic
value_0 = 0.1
value_residual = 0.1
internal_limit = 10
[../]
[./wp_compressive_str_soften]
type = SolidMechanicsHardeningCubic
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]
time_step_interval = 1
print_linear_residuals = false
exodus = true
csv = true
console = true
[]
(modules/solid_mechanics/test/tests/stickyBC/push_up.i)
# Testing StickyBC
#
# Push the bottom of an element upward until the top hits an (invisible) obstruction.
# 10 timesteps are used. In each timestep disp_y is increased by 0.1. The
# StickyBC has a max_value of 0.49, so at timestep 5 this bound will be violated
# and the top boundary will be fixed forever after.
#
# This test also illustrates that StickyBC is only ever meant to be used in
# special situations:
# - if, after the simulation ends, the bottom is moved downward again, the StickyBC
# will keep the top fixed. Ie, the StickyBC is truly "sticky".
# - setting max_value = 0.5 in this test illustrates the "approximate" nature
# of StickyBC, in that some nodes will be fixed at disp_y=0.5, while others
# will be fixed at disp_y=0.6, due to the timestepping and roundoff errors
# in MOOSE's solution.
[Mesh]
type = GeneratedMesh
dim = 3
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
[../]
[]
[BCs]
[./obstruction]
type = StickyBC
variable = disp_y
boundary = top
max_value = 0.49
[../]
[./bottom]
type = FunctionDirichletBC
variable = disp_y
boundary = bottom
function = t
[../]
[./left]
type = DirichletBC
variable = disp_x
boundary = left
value = 0
[../]
[./front]
type = DirichletBC
variable = disp_z
boundary = front
value = 0
[../]
[]
[Materials]
[./stress]
type = ComputeLinearElasticStress
[../]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 1.0
poissons_ratio = 0.2
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
solve_type = Linear
dt = 0.1
end_time = 1.0
[]
[Outputs]
exodus = true
[]
(modules/solid_mechanics/test/tests/stickyBC/push_down.i)
# Testing StickyBC
#
# Push the top of an element downward until the bottom hits an (invisible) obstruction.
# 10 timesteps are used. In each timestep disp_y is decreased by 0.1. The
# StickyBC has a min_value of -0.49, so at timestep 5 this bound will be violated
# and the bottom boundary will be fixed forever after.
#
# This test also illustrates that StickyBC is only ever meant to be used in
# special situations:
# - if, after the simulation ends, the top is moved upward again, the StickyBC
# will keep the bottom fixed. Ie, the StickyBC is truly "sticky".
# - setting min_value = -0.5 in this test illustrates the "approximate" nature
# of StickyBC, in that some nodes will be fixed at disp_y=-0.5, while others
# will be fixed at disp_y=-0.6, due to the timestepping and roundoff errors
# in MOOSE's solution.
[Mesh]
type = GeneratedMesh
dim = 3
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
[../]
[]
[BCs]
[./obstruction]
type = StickyBC
variable = disp_y
boundary = bottom
min_value = -0.49
[../]
[./top]
type = FunctionDirichletBC
variable = disp_y
boundary = top
function = -t
[../]
[./left]
type = DirichletBC
variable = disp_x
boundary = left
value = 0
[../]
[./front]
type = DirichletBC
variable = disp_z
boundary = front
value = 0
[../]
[]
[Materials]
[./stress]
type = ComputeLinearElasticStress
[../]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 1.0
poissons_ratio = 0.2
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
solve_type = Linear
dt = 0.1
end_time = 1.0
[]
[Outputs]
exodus = true
[]
(modules/solid_mechanics/test/tests/stickyBC/except1.i)
# Exception testing of StickyBC. Here min_value > max_value
[Mesh]
type = GeneratedMesh
dim = 3
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
[../]
[]
[BCs]
[./obstruction]
type = StickyBC
variable = disp_y
boundary = bottom
min_value = 1
max_value = -1
[../]
[]
[Materials]
[./stress]
type = ComputeLinearElasticStress
[../]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 1.0
poissons_ratio = 0.2
[../]
[]
[Executioner]
type = Transient
[]
(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
primary_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
symbol_names = 'bulk p0 g rho0'
symbol_values = '2E3 0.0 1E-5 1E3'
expression = '-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk)'
[]
[ini_xx]
type = ParsedFunction
symbol_names = 'bulk p0 g rho0 biot'
symbol_values = '2E3 0.0 1E-5 1E3 0.7'
expression = '0.8*(2500*10E-6*z+biot*(-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk)))'
[]
[ini_zz]
type = ParsedFunction
symbol_names = 'bulk p0 g rho0 biot'
symbol_values = '2E3 0.0 1E-5 1E3 0.7'
expression = '2500*10E-6*z+biot*(-bulk*log(exp(-p0/bulk)+g*rho0*z/bulk))'
[]
[excav_sideways]
type = ParsedFunction
symbol_names = 'end_t ymin ymax minval maxval slope'
symbol_values = '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
expression = '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
symbol_names = 'end_t ymin ymax minval maxval'
symbol_values = '0.5 0 1000.0 0 2500'
expression = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,maxval)'
[]
[roof_conductance]
type = ParsedFunction
symbol_names = 'end_t ymin ymax maxval minval'
symbol_values = '0.5 0 1000.0 1E7 0'
expression = '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
[]
[]
[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]
time_step_interval = 1
print_linear_residuals = true
exodus = true
csv = true
console = true
[]
(modules/solid_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
expression = '0.8*2500*10E-6*z'
[../]
[./ini_zz]
type = ParsedFunction
expression = '2500*10E-6*z'
[../]
[./excav_sideways]
type = ParsedFunction
symbol_names = 'end_t ymin ymax minval maxval slope'
symbol_values = '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
expression = '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
symbol_names = 'end_t ymin ymax minval maxval'
symbol_values = '17.0 0 1000.0 0 2500'
expression = 'if(y<ymin+(ymax-ymin)*min(t/end_t,1),minval,maxval)'
[../]
[]
[UserObjects]
[./mc_coh_strong_harden]
type = SolidMechanicsHardeningExponential
value_0 = 2.99 # MPa
value_residual = 3.01 # MPa
rate = 1.0
[../]
[./mc_fric]
type = SolidMechanicsHardeningConstant
value = 0.65 # 37deg
[../]
[./mc_dil]
type = SolidMechanicsHardeningConstant
value = 0.15 # 8deg
[../]
[./mc_tensile_str_strong_harden]
type = SolidMechanicsHardeningExponential
value_0 = 1.0 # MPa
value_residual = 1.0 # MPa
rate = 1.0
[../]
[./mc_compressive_str]
type = SolidMechanicsHardeningCubic
value_0 = 100 # Large!
value_residual = 100
internal_limit = 0.1
[../]
[./wp_coh_harden]
type = SolidMechanicsHardeningCubic
value_0 = 0.1
value_residual = 0.1
internal_limit = 10
[../]
[./wp_tan_fric]
type = SolidMechanicsHardeningConstant
value = 0.36 # 20deg
[../]
[./wp_tan_dil]
type = SolidMechanicsHardeningConstant
value = 0.18 # 10deg
[../]
[./wp_tensile_str_harden]
type = SolidMechanicsHardeningCubic
value_0 = 0.1
value_residual = 0.1
internal_limit = 10
[../]
[./wp_compressive_str_soften]
type = SolidMechanicsHardeningCubic
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]
time_step_interval = 1
print_linear_residuals = false
exodus = true
csv = true
console = true
[]