- valueValue multiplied against the residual, e.g. gravitational acceleration
C++ Type:double
Controllable:No
Description:Value multiplied against the residual, e.g. gravitational acceleration
- 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
Gravity
Apply gravity. Value is in units of acceleration.
Description
The kernel Gravity
provides a body force term in the stress divergence equilibrium equation to account for gravity force due to self weight. In a continuum mechanics problem, momentum conservation is prescribed assuming static equilibrium at each time increment, where is the Cauchy stress tensor and is the gravity body force per unit mass.
Example Input File Syntax
[Kernels]
[gravity_y]
type = Gravity
variable = disp_y
value = -9.81
[]
[]
(modules/solid_mechanics/test/tests/gravity/gravity_test.i)Input Parameters
- alpha0alpha parameter required for HHT time integration scheme
Default:0
C++ Type:double
Controllable:No
Description:alpha parameter required for HHT time integration scheme
- blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
- densitydensityThe density
Default:density
C++ Type:MaterialPropertyName
Controllable:No
Description:The density
- displacementsThe displacements
C++ Type:std::vector<VariableName>
Controllable:No
Description:The displacements
- function1A function that describes the gravitational force
Default:1
C++ Type:FunctionName
Controllable:No
Description:A function that describes the gravitational force
- prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
- use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
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_tagssystemThe tag for the matrices this Kernel should fill
Default:system
C++ Type:MultiMooseEnum
Options:nontime, system
Controllable:No
Description:The tag for the matrices this Kernel should fill
- vector_tagsnontimeThe tag for the vectors this Kernel should fill
Default:nontime
C++ Type:MultiMooseEnum
Options:nontime, time
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.
- diag_save_inThe name of auxiliary variables to save this Kernel'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 Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
- 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
- save_inThe name of auxiliary variables to save this Kernel'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 Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
- 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_meshTrueDisplaced mesh defaults to true
Default:True
C++ Type:bool
Controllable:No
Description:Displaced mesh defaults to true
Advanced Parameters
Input Files
- (modules/solid_mechanics/examples/coal_mining/cosserat_mc_wp_sticky.i)
- (modules/combined/test/tests/gravity/gravity.i)
- (modules/solid_mechanics/test/tests/gravity/block-gravity-kinetic-energy.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/beam.i)
- (modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)
- (modules/combined/test/tests/gravity/gravity_qp_select.i)
- (modules/solid_mechanics/examples/coal_mining/coarse.i)
- (modules/contact/test/tests/explicit_dynamics/test_balance.i)
- (modules/solid_mechanics/examples/coal_mining/cosserat_mc_wp_sticky_longitudinal.i)
- (modules/solid_mechanics/examples/coal_mining/fine.i)
- (modules/porous_flow/examples/coal_mining/fine_with_fluid.i)
- (modules/solid_mechanics/examples/coal_mining/cosserat_elastic.i)
- (modules/solid_mechanics/examples/coal_mining/cosserat_mc_only.i)
- (modules/combined/test/tests/gravity/gravity_hex20.i)
- (modules/solid_mechanics/examples/coal_mining/cosserat_mc_wp.i)
- (modules/solid_mechanics/examples/coal_mining/cosserat_wp_only.i)
- (modules/solid_mechanics/test/tests/gravity/ad_gravity_test.i)
- (modules/contact/test/tests/explicit_dynamics/settlement.i)
- (modules/solid_mechanics/test/tests/dynamics/prescribed_displacement/3D_QStatic_1_Ramped_Displacement_with_gravity.i)
- (modules/solid_mechanics/test/tests/gravity/gravity_test.i)
- (modules/solid_mechanics/examples/bridge/bridge.i)
- (modules/solid_mechanics/examples/bridge/bridge_large_strain.i)
- (modules/combined/test/tests/gravity/gravity_rz.i)
- (modules/combined/test/tests/gravity/gravity_rz_quad8.i)
(modules/solid_mechanics/test/tests/gravity/gravity_test.i)
#
# Gravity Test
#
# This test is designed to apply a gravity body force.
#
# The mesh is composed of one block with a single element.
# The bottom is fixed in all three directions. Poisson's ratio
# is zero and the density is 20/9.81
# which makes it trivial to check displacements.
#
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Mesh]
type = GeneratedMesh
dim = 3
[]
[Physics/SolidMechanics/QuasiStatic]
[all]
add_variables = true
[]
[]
[Kernels]
[gravity_y]
type = Gravity
variable = disp_y
value = -9.81
[]
[]
[BCs]
[no_x]
type = DirichletBC
variable = disp_x
boundary = bottom
value = 0.0
[]
[no_y]
type = DirichletBC
variable = disp_y
boundary = bottom
value = 0.0
[]
[no_z]
type = DirichletBC
variable = disp_z
boundary = bottom
value = 0.0
[]
[]
[Materials]
[Elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0 0.5e6'
[]
[stress]
type = ComputeLinearElasticStress
[]
[density]
type = GenericConstantMaterial
prop_names = density
prop_values = 2.0387
[]
[]
[Executioner]
type = Steady
nl_abs_tol = 1e-10
l_max_its = 20
[]
[Outputs]
[out]
type = Exodus
elemental_as_nodal = true
[]
[]
(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/combined/test/tests/gravity/gravity.i)
# Gravity Test
#
# This test is designed to exercise the gravity body force kernel.
#
# The mesh for this problem is a rectangular bar 10 units by 1 unit
# by 1 unit.
#
# The boundary conditions for this problem are as follows. The
# displacement is zero on each of side that faces a negative
# coordinate direction. The acceleration of gravity is 20.
#
# The material has a Young's modulus of 1e6 and a density of 2.
#
# The analytic solution for the displacement along the bar is:
#
# u(x) = -b*x^2/(2*E)+b*L*x/E
#
# The displacement at x=L is b*L^2/(2*E) = 2*20*10*10/(2*1e6) = 0.002.
#
# The analytic solution for the stress along the bar assuming linear
# elasticity is:
#
# S(x) = b*(L-x)
#
# The stress at x=0 is b*L = 2*20*10 = 400.
#
# Note: The simulation does not measure stress at x=0. The stress
# is reported at element centers. The element closest to x=0 sits
# at x = 1/4 and has a stress of 390. This matches the linear
# stress distribution that is expected. The same situation applies
# at x = L where the stress is zero analytically. The nearest
# element is at x=9.75 where the stress is 10.
#
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Mesh]
file = gravity_test.e
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Modules/TensorMechanics/Master/All]
volumetric_locking_correction = true
strain = FINITE
add_variables = true
generate_output = 'stress_xx'
[]
[Kernels]
[./gravity]
type = Gravity
variable = disp_x
value = 20
[../]
[]
[BCs]
[./no_x]
type = DirichletBC
variable = disp_x
boundary = 1
value = 0.0
[../]
[./no_y]
type = DirichletBC
variable = disp_y
boundary = 3
value = 0.0
[../]
[./no_z]
type = DirichletBC
variable = disp_z
boundary = 5
value = 0.0
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 1e6
bulk_modulus = 0.333333333333333e6
[../]
[./stress]
type = ComputeFiniteStrainElasticStress
[../]
[./density]
type = Density
block = 1
density = 2
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
start_time = 0.0
end_time = 1.0
[]
[Outputs]
file_base = gravity_out
[./exodus]
type = Exodus
elemental_as_nodal = true
[../]
[]
(modules/solid_mechanics/test/tests/gravity/block-gravity-kinetic-energy.i)
starting_point = 2e-1
offset = 1.0
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
file = long-bottom-block-1elem-blocks.e
[]
[Problem]
kernel_coverage_check = false
material_coverage_check = false
[]
[Variables]
[disp_x]
block = '1 2'
[]
[disp_y]
block = '1 2'
[]
[]
[AuxVariables]
[pid]
order = CONSTANT
family = MONOMIAL
[]
[kinetic_energy]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[pid]
type = ProcessorIDAux
variable = pid
execute_on = 'initial timestep_end'
[]
[kinetic_energy]
type = KineticEnergyAux
block = '1 2'
variable = kinetic_energy
newmark_velocity_x = vel_x
newmark_velocity_y = vel_y
newmark_velocity_z = 0.0
density = density
[]
[]
[ICs]
[disp_y]
type = ConstantIC
block = 2
variable = disp_y
value = '${fparse starting_point + offset}'
[]
[]
[Physics/SolidMechanics/Dynamic]
[all]
add_variables = true
hht_alpha = 0.0
beta = 0.25
gamma = 0.5
mass_damping_coefficient = 0.0
stiffness_damping_coefficient = 0.0
displacements = 'disp_x disp_y'
generate_output = 'stress_xx stress_yy'
block = '1 2'
strain = FINITE
[]
[]
[Kernels]
[gravity]
type = Gravity
value = -9.81
variable = disp_y
[]
[]
[Materials]
[elasticity_2]
type = ComputeIsotropicElasticityTensor
block = '2'
youngs_modulus = 1e4
poissons_ratio = 0.3
[]
[elasticity_1]
type = ComputeIsotropicElasticityTensor
block = '1'
youngs_modulus = 1e7
poissons_ratio = 0.3
[]
[stress]
type = ComputeFiniteStrainElasticStress
block = '1 2'
[]
[density]
type = GenericConstantMaterial
block = '1 2'
prop_names = 'density'
prop_values = '7750'
[]
[]
[BCs]
[botx]
type = DirichletBC
variable = disp_x
boundary = 40
value = 0.0
[]
[boty]
type = DirichletBC
variable = disp_y
boundary = 40
value = 0.0
[]
[]
[Executioner]
type = Transient
end_time = 0.5
dt = 0.01
dtmin = .05
solve_type = 'PJFNK'
petsc_options = '-snes_converged_reason -ksp_converged_reason -pc_svd_monitor '
'-snes_linesearch_monitor'
petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount -mat_mffd_err '
'-ksp_gmres_restart'
petsc_options_value = 'lu NONZERO 1e-15 1e-5 100'
l_max_its = 100
nl_max_its = 20
line_search = 'none'
snesmf_reuse_base = false
[TimeIntegrator]
type = NewmarkBeta
beta = 0.25
gamma = 0.5
[]
[]
[Debug]
show_var_residual_norms = true
[]
[Outputs]
exodus = false
csv = true
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Postprocessors]
active = 'total_kinetic_energy'
[total_kinetic_energy]
type = ElementIntegralVariablePostprocessor
variable = kinetic_energy
block = '1 2'
[]
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/beam.i)
# A beam with its ends fully clamped
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 10
nz = 10
xmin = -10
xmax = 10
ymin = -10
ymax = 10
zmin = -50
zmax = 0
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Kernels]
[SolidMechanics]
[../]
[./gravity_y]
type = Gravity
use_displaced_mesh = false
variable = disp_y
value = -10
[../]
[]
[BCs]
[./zmax_xfixed]
type = DirichletBC
variable = disp_x
boundary = front
value = 0
[../]
[./zmax_yfixed]
type = DirichletBC
variable = disp_y
boundary = front
value = 0
[../]
[./zmax_zfixed]
type = DirichletBC
variable = disp_z
boundary = front
value = 0
[../]
[./zmin_xfixed]
type = DirichletBC
variable = disp_x
boundary = back
value = 0
[../]
[./zmin_yfixed]
type = DirichletBC
variable = disp_y
boundary = back
value = 0
[../]
[./zmin_zfixed]
type = DirichletBC
variable = disp_z
boundary = back
value = 0
[../]
[]
[AuxVariables]
[./stress_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xz]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yz]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./strainp_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./strainp_xy]
order = CONSTANT
family = MONOMIAL
[../]
[./strainp_xz]
order = CONSTANT
family = MONOMIAL
[../]
[./strainp_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./strainp_yz]
order = CONSTANT
family = MONOMIAL
[../]
[./strainp_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./straint_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./straint_xy]
order = CONSTANT
family = MONOMIAL
[../]
[./straint_xz]
order = CONSTANT
family = MONOMIAL
[../]
[./straint_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./straint_yz]
order = CONSTANT
family = MONOMIAL
[../]
[./straint_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./f_shear]
order = CONSTANT
family = MONOMIAL
[../]
[./f_tensile]
order = CONSTANT
family = MONOMIAL
[../]
[./f_compressive]
order = CONSTANT
family = MONOMIAL
[../]
[./intnl_shear]
order = CONSTANT
family = MONOMIAL
[../]
[./intnl_tensile]
order = CONSTANT
family = MONOMIAL
[../]
[./iter]
order = CONSTANT
family = MONOMIAL
[../]
[./ls]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./stress_xx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx
index_i = 0
index_j = 0
[../]
[./stress_xy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xy
index_i = 0
index_j = 1
[../]
[./stress_xz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xz
index_i = 0
index_j = 2
[../]
[./stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[../]
[./stress_yz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yz
index_i = 1
index_j = 2
[../]
[./stress_zz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zz
index_i = 2
index_j = 2
[../]
[./strainp_xx]
type = RankTwoAux
rank_two_tensor = plastic_strain
variable = strainp_xx
index_i = 0
index_j = 0
[../]
[./strainp_xy]
type = RankTwoAux
rank_two_tensor = plastic_strain
variable = strainp_xy
index_i = 0
index_j = 1
[../]
[./strainp_xz]
type = RankTwoAux
rank_two_tensor = plastic_strain
variable = strainp_xz
index_i = 0
index_j = 2
[../]
[./strainp_yy]
type = RankTwoAux
rank_two_tensor = plastic_strain
variable = strainp_yy
index_i = 1
index_j = 1
[../]
[./strainp_yz]
type = RankTwoAux
rank_two_tensor = plastic_strain
variable = strainp_yz
index_i = 1
index_j = 2
[../]
[./strainp_zz]
type = RankTwoAux
rank_two_tensor = plastic_strain
variable = strainp_zz
index_i = 2
index_j = 2
[../]
[./straint_xx]
type = RankTwoAux
rank_two_tensor = total_strain
variable = straint_xx
index_i = 0
index_j = 0
[../]
[./straint_xy]
type = RankTwoAux
rank_two_tensor = total_strain
variable = straint_xy
index_i = 0
index_j = 1
[../]
[./straint_xz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = straint_xz
index_i = 0
index_j = 2
[../]
[./straint_yy]
type = RankTwoAux
rank_two_tensor = total_strain
variable = straint_yy
index_i = 1
index_j = 1
[../]
[./straint_yz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = straint_yz
index_i = 1
index_j = 2
[../]
[./straint_zz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = straint_zz
index_i = 2
index_j = 2
[../]
[./f_shear]
type = MaterialStdVectorAux
property = plastic_yield_function
index = 0
variable = f_shear
[../]
[./f_tensile]
type = MaterialStdVectorAux
property = plastic_yield_function
index = 1
variable = f_tensile
[../]
[./f_compressive]
type = MaterialStdVectorAux
property = plastic_yield_function
index = 2
variable = f_compressive
[../]
[./intnl_shear]
type = MaterialStdVectorAux
property = plastic_internal_parameter
index = 0
variable = intnl_shear
[../]
[./intnl_tensile]
type = MaterialStdVectorAux
property = plastic_internal_parameter
index = 1
variable = intnl_tensile
[../]
[./iter]
type = MaterialRealAux
property = plastic_NR_iterations
variable = iter
[../]
[./ls]
type = MaterialRealAux
property = plastic_linesearch_needed
variable = ls
[../]
[]
[UserObjects]
[./coh_irrelevant]
type = SolidMechanicsHardeningCubic
value_0 = 2E6
value_residual = 2E6
internal_limit = 0.01
[../]
[./tanphi]
type = SolidMechanicsHardeningCubic
value_0 = 0.5
value_residual = 0.5
internal_limit = 0.01
[../]
[./tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.166666666667
[../]
[./t_strength]
type = SolidMechanicsHardeningCubic
value_0 = 0
value_residual = 0
internal_limit = 0.1
[../]
[./c_strength]
type = SolidMechanicsHardeningCubic
value_0 = 1E80
value_residual = 0.0
internal_limit = 0.01
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '6.4E9 6.4E9' # young 16MPa, Poisson 0.25
[../]
[./strain]
type = ComputeIncrementalSmallStrain
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
tangent_operator = nonlinear
perform_finite_strain_rotations = false
[../]
[./stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh_irrelevant
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 1000
tip_smoother = 1E5
smoothing_tol = 1E5
yield_function_tol = 1E-5
perfect_guess = true
min_step_size = 0.1
[../]
[./density]
type = GenericFunctionMaterial
block = 0
prop_names = density
prop_values = 1E3*t
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
#petsc_options = '-snes_converged_reason -snes_linesearch_monitor'
petsc_options_iname = '-pc_type -pc_asm_overlap -sub_pc_type -ksp_type -ksp_gmres_restart'
petsc_options_value = ' asm 2 lu gmres 200'
[../]
[]
[Executioner]
solve_type = 'NEWTON'
petsc_options = '-snes_converged_reason'
line_search = bt
nl_abs_tol = 1E-2
nl_rel_tol = 1e-15
l_tol = 1E-10
l_max_its = 100
nl_max_its = 100
end_time = 10
dt = 1
type = Transient
[]
[Outputs]
file_base = beam
exodus = true
csv = true
[]
(modules/porous_flow/examples/coal_mining/coarse_with_fluid.i)
# Strata deformation and fluid flow aaround a coal mine - 3D model
#
# A "half model" is used. The mine is 400m deep and
# just the roof is studied (-400<=z<=0). The mining panel
# sits between 0<=x<=150, and 0<=y<=1000, so this simulates
# a coal panel that is 300m wide and 1000m long. The outer boundaries
# are 1km from the excavation boundaries.
#
# The excavation takes 0.5 years.
#
# The boundary conditions for this simulation are:
# - disp_x = 0 at x=0 and x=1150
# - disp_y = 0 at y=-1000 and y=1000
# - disp_z = 0 at z=-400, but there is a time-dependent
# Young modulus that simulates excavation
# - wc_x = 0 at y=-1000 and y=1000
# - wc_y = 0 at x=0 and x=1150
# - no flow at x=0, z=-400 and z=0
# - fixed porepressure at y=-1000, y=1000 and x=1150
# That is, rollers on the sides, free at top,
# and prescribed at bottom in the unexcavated portion.
#
# A single-phase unsaturated fluid is used.
#
# The small strain formulation is used.
#
# All stresses are measured in MPa, and time units are measured in years.
#
# The initial porepressure is hydrostatic with P=0 at z=0, so
# Porepressure ~ - 0.01*z MPa, where the fluid has density 1E3 kg/m^3 and
# gravity = = 10 m.s^-2 = 1E-5 MPa m^2/kg.
# To be more accurate, i use
# Porepressure = -bulk * log(1 + g*rho0*z/bulk)
# where bulk=2E3 MPa and rho0=1Ee kg/m^3.
# The initial stress is consistent with the weight force from undrained
# density 2500 kg/m^3, and fluid porepressure, and a Biot coefficient of 0.7, ie,
# stress_zz^effective = 0.025*z + 0.7 * initial_porepressure
# The maximum and minimum principal horizontal effective stresses are
# assumed to be equal to 0.8*stress_zz.
#
# Material properties:
# Young's modulus = 8 GPa
# Poisson's ratio = 0.25
# Cosserat layer thickness = 1 m
# Cosserat-joint normal stiffness = large
# Cosserat-joint shear stiffness = 1 GPa
# MC cohesion = 2 MPa
# MC friction angle = 35 deg
# MC dilation angle = 8 deg
# MC tensile strength = 1 MPa
# MC compressive strength = 100 MPa
# WeakPlane cohesion = 0.1 MPa
# WeakPlane friction angle = 30 deg
# WeakPlane dilation angle = 10 deg
# WeakPlane tensile strength = 0.1 MPa
# WeakPlane compressive strength = 100 MPa softening to 1 MPa at strain = 1
# Fluid density at zero porepressure = 1E3 kg/m^3
# Fluid bulk modulus = 2E3 MPa
# Fluid viscosity = 1.1E-3 Pa.s = 1.1E-9 MPa.s = 3.5E-17 MPa.year
#
[GlobalParams]
perform_finite_strain_rotations = false
displacements = 'disp_x disp_y disp_z'
Cosserat_rotations = 'wc_x wc_y wc_z'
PorousFlowDictator = dictator
biot_coefficient = 0.7
[]
[Mesh]
[file]
type = FileMeshGenerator
file = mesh/coarse.e
[]
[xmin]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
new_boundary = xmin
normal = '-1 0 0'
input = file
[]
[xmax]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
new_boundary = xmax
normal = '1 0 0'
input = xmin
[]
[ymin]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
new_boundary = ymin
normal = '0 -1 0'
input = xmax
[]
[ymax]
type = SideSetsAroundSubdomainGenerator
block = '2 3 4 5 6 7 8 9 10 11 12 13 14 15 16'
new_boundary = ymax
normal = '0 1 0'
input = ymin
[]
[zmax]
type = SideSetsAroundSubdomainGenerator
block = 16
new_boundary = zmax
normal = '0 0 1'
input = ymax
[]
[zmin]
type = SideSetsAroundSubdomainGenerator
block = 2
new_boundary = zmin
normal = '0 0 -1'
input = zmax
[]
[excav]
type = SubdomainBoundingBoxGenerator
input = zmin
block_id = 1
bottom_left = '0 0 -400'
top_right = '150 1000 -397'
[]
[roof]
type = SideSetsBetweenSubdomainsGenerator
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/combined/test/tests/gravity/gravity_qp_select.i)
# Gravity Test
#
# This test is similar to the other gravity tests, but it also tests the
# capability in MaterialTensorAux to return the stress of a single,
# specified integration point, rather than the element average.
# To get the stress at a single integration point, set the parameter
# qp_select to the integration point number (i.e. 0-9 for a quad 8)
# in the AuxKernel
#
# The mesh for this problem is a unit square.
#
# The boundary conditions for this problem are as follows. The
# displacement is zero on each of side that faces a negative
# coordinate direction. The acceleration of gravity is 20.
#
# The material has a Young's modulus of 1e6 and a density of 2.
#
# The analytic solution for the displacement along the bar is:
#
# u(x) = -b*x^2/(2*E)+b*L*x/E
#
# The displacement at x=L is b*L^2/(2*E) = 2*20*1*1/(2*1e6) = 0.00002.
#
# The analytic solution for the stress along the bar assuming linear
# elasticity is:
#
# S(x) = b*(L-x)
#
# The stress at x=0 is b*L = 2*20*1 = 40.
#
# Note: The isoparametric coordinate for a quad8 (fourth order) element
# is: +/- 0.77459667 and 0. For a 1 unit square with the edge of
# the element in the x = 0 plane, there would be an integration point
# at x_coordinate 0.5 - 0.5*0.77459667 (0.11270167), 0.5, and
# 0.50 + 0.5*0.77459667 (0.88729834).
#
# The corresponding stresses are:
#
# S(0.11270167) = 40(1-0.11270167) = 35.491933
# S(0.5) = 40(1-0.5) = 20
# S(0.88729834) = 40(1-0.88729834) = 4.5080664
#
# These stresses are a precise match to the simulation result.
#
[GlobalParams]
displacements = 'disp_x disp_y'
order = SECOND
family = LAGRANGE
[]
[Mesh]
file = gravity_2D.e
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[]
[AuxVariables]
[./stress_xx_qp_0]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xx_qp_1]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xx_qp_2]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xx_qp_3]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xx_qp_4]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xx_qp_5]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xx_qp_6]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xx_qp_7]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xx_qp_8]
order = CONSTANT
family = MONOMIAL
[../]
[]
[Modules/TensorMechanics/Master/All]
strain = FINITE
#incremental = true
add_variables = true
generate_output = 'stress_xx'
[]
[Kernels]
[./gravity]
type = Gravity
variable = disp_x
value = 20
[../]
[]
[AuxKernels]
[./stress_xx_qp_0]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx_qp_0
index_i = 0
index_j = 0
selected_qp = 0
[../]
[./stress_xx_qp_1]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx_qp_1
index_i = 0
index_j = 0
selected_qp = 1
[../]
[./stress_xx_qp_2]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx_qp_2
index_i = 0
index_j = 0
selected_qp = 2
[../]
[./stress_xx_qp_3]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx_qp_3
index_i = 0
index_j = 0
selected_qp = 3
[../]
[./stress_xx_qp_4]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx_qp_4
index_i = 0
index_j = 0
selected_qp = 4
[../]
[./stress_xx_qp_5]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx_qp_5
index_i = 0
index_j = 0
selected_qp = 5
[../]
[./stress_xx_qp_6]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx_qp_6
index_i = 0
index_j = 0
selected_qp = 6
[../]
[./stress_xx_qp_7]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx_qp_7
index_i = 0
index_j = 0
selected_qp = 7
[../]
[./stress_xx_qp_8]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx_qp_8
index_i = 0
index_j = 0
selected_qp = 8
[../]
[]
[BCs]
[./no_x]
type = DirichletBC
variable = disp_x
boundary = 1
value = 0.0
[../]
[./no_z]
type = DirichletBC
variable = disp_y
boundary = 5
value = 0.0
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 1e6
bulk_modulus = 0.333333333333333e6
[../]
[./stress]
type = ComputeFiniteStrainElasticStress
[../]
[./density]
type = Density
density = 2
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
start_time = 0.0
end_time = 1.0
[]
[Outputs]
file_base = gravity_qp_select_out
[./exodus]
type = Exodus
elemental_as_nodal = 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
[]
(modules/contact/test/tests/explicit_dynamics/test_balance.i)
# One element test to test the central difference time integrator in 3D.
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
volumetric_locking_correction = true
[]
[Mesh]
[block_one]
type = GeneratedMeshGenerator
dim = 3
nx = 3
ny = 3
nz = 3
xmin = 4.5
xmax = 5.5
ymin = 4.5
ymax = 5.5
zmin = 0.0001
zmax = 1.0001
boundary_name_prefix = 'ball'
[]
[block_two]
type = GeneratedMeshGenerator
dim = 3
nx = 2
ny = 2
nz = 2
xmin = 0.0
xmax = 10
ymin = 0.0
ymax = 10
zmin = -2
zmax = 0
boundary_name_prefix = 'base'
boundary_id_offset = 10
[]
[block_one_id]
type = SubdomainIDGenerator
input = block_one
subdomain_id = 1
[]
[block_two_id]
type = SubdomainIDGenerator
input = block_two
subdomain_id = 2
[]
[combine]
type = MeshCollectionGenerator
inputs = ' block_one_id block_two_id'
[]
allow_renumbering = false
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[]
[AuxVariables]
[gap_rate]
[]
[vel_x]
[]
[accel_x]
[]
[vel_y]
[]
[accel_y]
[]
[vel_z]
[]
[accel_z]
[]
[stress_zz]
family = MONOMIAL
order = CONSTANT
[]
[strain_zz]
family = MONOMIAL
order = CONSTANT
[]
[]
[AuxKernels]
[stress_zz]
type = RankTwoAux
rank_two_tensor = stress
index_i = 2
index_j = 2
variable = stress_zz
execute_on = 'TIMESTEP_END'
[]
[strain_zz]
type = RankTwoAux
rank_two_tensor = mechanical_strain
index_i = 2
index_j = 2
variable = strain_zz
[]
[accel_x]
type = TestNewmarkTI
variable = accel_x
displacement = disp_x
first = false
execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
[]
[vel_x]
type = TestNewmarkTI
variable = vel_x
displacement = disp_x
execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
[]
[accel_y]
type = TestNewmarkTI
variable = accel_y
displacement = disp_y
first = false
execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
[]
[vel_y]
type = TestNewmarkTI
variable = vel_y
displacement = disp_x
execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
[]
[accel_z]
type = TestNewmarkTI
variable = accel_z
displacement = disp_z
first = false
execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
[]
[vel_z]
type = TestNewmarkTI
variable = vel_z
displacement = disp_z
execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
[]
[]
[Kernels]
[DynamicTensorMechanics]
displacements = 'disp_x disp_y disp_z'
volumetric_locking_correction = true
stiffness_damping_coefficient = 0.001
generate_output = 'stress_zz strain_zz'
[]
[inertia_x]
type = InertialForce
variable = disp_x
[]
[inertia_y]
type = InertialForce
variable = disp_y
[]
[inertia_z]
type = InertialForce
variable = disp_z
[]
[]
[Kernels]
[gravity]
type = Gravity
variable = disp_z
value = -981.0
[]
[]
[BCs]
[x_front]
type = DirichletBC
variable = disp_x
boundary = 'ball_front'
preset = false
value = 0.0
[]
[y_front]
type = DirichletBC
variable = disp_y
boundary = 'ball_front'
preset = false
value = 0.0
[]
[x_fixed]
type = DirichletBC
variable = disp_x
boundary = 'base_back'
preset = false
value = 0.0
[]
[y_fixed]
type = DirichletBC
variable = disp_y
boundary = 'base_back'
preset = false
value = 0.0
[]
[z_fixed]
type = DirichletBC
variable = disp_z
boundary = 'base_back'
preset = false
value = 0.0
[]
[z_fixed_front]
type = DirichletBC
variable = disp_z
boundary = 'base_front'
preset = false
value = 0.0
[]
[]
[ExplicitDynamicsContact]
[my_contact]
model = frictionless_balance
primary = base_front
secondary = ball_back
vel_x = 'vel_x'
vel_y = 'vel_y'
vel_z = 'vel_z'
[]
[]
[Materials]
[elasticity_tensor_block_one]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 1e6
poissons_ratio = 0.0
block = 1
outputs = 'exodus'
output_properties = __all__
[]
[elasticity_tensor_block_two]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 1e10
poissons_ratio = 0.0
block = 2
outputs = 'exodus'
output_properties = __all__
[]
[strain_block]
type = ComputeFiniteStrain
displacements = 'disp_x disp_y disp_z'
implicit = false
[]
[stress_block]
type = ComputeFiniteStrainElasticStress
[]
[density_one]
type = GenericConstantMaterial
prop_names = density
prop_values = 1e1
outputs = 'exodus'
output_properties = 'density'
block = '1'
[]
[density_two]
type = GenericConstantMaterial
prop_names = density
prop_values = 1e6
outputs = 'exodus'
output_properties = 'density'
block = '2'
[]
[wave_speed]
type = WaveSpeed
outputs = 'exodus'
output_properties = 'wave_speed'
[]
[]
[Executioner]
type = Transient
start_time = -0.01
end_time = -0.0075 # 10
dt = 0.00001
timestep_tolerance = 1e-6
[TimeIntegrator]
type = CentralDifference
solve_type = lumped
[]
[]
[Outputs]
interval = 50
exodus = true
csv = true
[]
[Postprocessors]
[accel_58z]
type = NodalVariableValue
nodeid = 1
variable = accel_z
[]
[vel_58z]
type = NodalVariableValue
nodeid = 1
variable = vel_z
[]
[disp_58z]
type = NodalVariableValue
nodeid = 1
variable = disp_z
[]
[critical_time_step]
type = CriticalTimeStep
[]
[contact_pressure_max]
type = NodalExtremeValue
variable = contact_pressure
block = '1 2'
value_type = max
[]
[]
(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/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/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/cosserat_elastic.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 for y>150
# - 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.
#
# This is an elastic simulation, but the weak-plane and Drucker-Prager
# parameters and AuxVariables may be found below. They are irrelevant
# in this simulation. The weak-plane and Drucker-Prager cohesions,
# tensile strengths and compressive strengths have been set very high
#
# 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
#
[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
new_boundary = 11
normal = '0 -1 0'
input = generated_mesh
[]
[right]
type = SideSetsAroundSubdomainGenerator
new_boundary = 12
normal = '0 1 0'
input = left
[]
[front]
type = SideSetsAroundSubdomainGenerator
new_boundary = 13
normal = '-1 0 0'
input = right
[]
[back]
type = SideSetsAroundSubdomainGenerator
new_boundary = 14
normal = '1 0 0'
input = front
[]
[top]
type = SideSetsAroundSubdomainGenerator
new_boundary = 15
normal = '0 0 1'
input = back
[]
[bottom]
type = SideSetsAroundSubdomainGenerator
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 = SideSetsBetweenSubdomainsGenerator
new_boundary = 21
primary_block = 0
paired_block = 1
input = excav
[]
[hole]
type = BlockDeletionGenerator
block = 1
input = roof
[]
[]
[GlobalParams]
block = 0
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
[../]
[./dp_shear]
order = CONSTANT
family = MONOMIAL
[../]
[./dp_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
[../]
[./dp_shear_f]
order = CONSTANT
family = MONOMIAL
[../]
[./dp_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
[../]
[./dp_shear]
type = MaterialStdVectorAux
index = 0
property = dp_plastic_internal_parameter
variable = dp_shear
[../]
[./dp_tensile]
type = MaterialStdVectorAux
index = 1
property = dp_plastic_internal_parameter
variable = dp_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
[../]
[./dp_shear_f]
type = MaterialStdVectorAux
index = 0
property = dp_plastic_yield_function
variable = dp_shear_f
[../]
[./dp_tensile_f]
type = MaterialStdVectorAux
index = 1
property = dp_plastic_yield_function
variable = dp_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
[../]
[]
[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)'
[../]
[]
[UserObjects]
[./dp_coh_strong_harden]
type = SolidMechanicsHardeningExponential
value_0 = 2.9 # MPa
value_residual = 3.1 # MPa
rate = 1.0
[../]
[./dp_fric]
type = SolidMechanicsHardeningConstant
value = 0.65 # 37deg
[../]
[./dp_dil]
type = SolidMechanicsHardeningConstant
value = 0.65
[../]
[./dp_tensile_str_strong_harden]
type = SolidMechanicsHardeningExponential
value_0 = 1.0 # MPa
value_residual = 1.4 # MPa
rate = 1.0
[../]
[./dp_compressive_str]
type = SolidMechanicsHardeningConstant
value = 1.0E3 # Large!
[../]
[./drucker_prager_model]
type = SolidMechanicsPlasticDruckerPrager
mc_cohesion = dp_coh_strong_harden
mc_friction_angle = dp_fric
mc_dilation_angle = dp_dil
internal_constraint_tolerance = 1 # irrelevant here
yield_function_tolerance = 1 # irrelevant here
[../]
[./wp_coh]
type = SolidMechanicsHardeningConstant
value = 1E12
[../]
[./wp_tan_fric]
type = SolidMechanicsHardeningConstant
value = 0.36 # 20deg
[../]
[./wp_tan_dil]
type = SolidMechanicsHardeningConstant
value = 0.18 # 10deg
[../]
[./wp_tensile_str]
type = SolidMechanicsHardeningConstant
value = 1E12
[../]
[./wp_compressive_str]
type = SolidMechanicsHardeningConstant
value = 1E12
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeLayeredCosseratElasticityTensor
young = 8E3 # MPa
poisson = 0.25
layer_thickness = 1.0
joint_normal_stiffness = 1E9 # huge
joint_shear_stiffness = 1E3 # MPa
[../]
[./strain]
type = ComputeCosseratIncrementalSmallStrain
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = 'ini_xx 0 0 0 ini_xx 0 0 0 ini_zz'
eigenstrain_name = ini_stress
[../]
[./stress]
# this is needed so as to correctly apply the initial stress
type = ComputeMultipleInelasticCosseratStress
block = 0
inelastic_models = ''
relative_tolerance = 2.0
absolute_tolerance = 1E6
max_iterations = 1
tangent_operator = nonlinear
perform_finite_strain_rotations = false
[../]
[./dp]
type = CappedDruckerPragerCosseratStressUpdate
block = 0
warn_about_precision_loss = false
host_youngs_modulus = 8E3
host_poissons_ratio = 0.25
base_name = dp
DP_model = drucker_prager_model
tensile_strength = dp_tensile_str_strong_harden
compressive_strength = dp_compressive_str
max_NR_iterations = 100000
tip_smoother = 0.1E1
smoothing_tol = 0.1E1 # MPa # Must be linked to cohesion
yield_function_tol = 1E-11 # MPa. this is essentially the lowest possible without lots of precision loss
perfect_guess = true
min_step_size = 1.0
[../]
[./wp]
type = CappedWeakPlaneCosseratStressUpdate
block = 0
warn_about_precision_loss = false
base_name = wp
cohesion = wp_coh
tan_friction_angle = wp_tan_fric
tan_dilation_angle = wp_tan_dil
tensile_strength = wp_tensile_str
compressive_strength = wp_compressive_str
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]
type = GenericConstantMaterial
prop_names = density
prop_values = 2500
[../]
[]
[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 = 'Linear'
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 = 1000
start_time = 0.0
dt = 1.0
end_time = 1.0
[]
[Outputs]
file_base = cosserat_elastic
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_only.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 300m deep
# and just the roof is studied (0<=z<=300). 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). Mining is simulated by moving the excavation's
# roof down, until disp_z=-3 at t=1.
# 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 are:
# - disp_x = 0 everywhere
# - disp_y = 0 at y=0 and y=450
# - disp_z = 0 for y>150
# - disp_z = -3 at maximum, for 0<=y<=150. See excav function.
# That is, rollers on the sides, free at top, and prescribed at bottom.
#
# 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.
#
# Below you will see weak-plane parameters and AuxVariables, etc.
# These are not actally used in this example.
#
# 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
#
[Mesh]
[generated_mesh]
type = GeneratedMeshGenerator
dim = 3
nx = 1
xmin = -5
xmax = 5
nz = 40
zmin = 0
zmax = 400.0
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
new_boundary = 11
normal = '0 -1 0'
input = generated_mesh
[]
[right]
type = SideSetsAroundSubdomainGenerator
new_boundary = 12
normal = '0 1 0'
input = left
[]
[front]
type = SideSetsAroundSubdomainGenerator
new_boundary = 13
normal = '-1 0 0'
input = right
[]
[back]
type = SideSetsAroundSubdomainGenerator
new_boundary = 14
normal = '1 0 0'
input = front
[]
[top]
type = SideSetsAroundSubdomainGenerator
new_boundary = 15
normal = '0 0 1'
input = back
[]
[bottom]
type = SideSetsAroundSubdomainGenerator
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 = SideSetsBetweenSubdomainsGenerator
new_boundary = 21
primary_block = 0
paired_block = 1
input = excav
[]
[hole]
type = BlockDeletionGenerator
block = 1
input = roof
[]
[]
[GlobalParams]
block = 0
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
[../]
[]
[AuxVariables]
[./disp_x]
[../]
[./wc_y]
[../]
[./wc_z]
[../]
[./stress_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yy]
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_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
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 16 21' # note addition of 16 and 21
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 = FunctionDirichletBC
variable = disp_z
boundary = 21
function = excav_sideways
[../]
[]
[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 e_h closure_dist'
symbol_values = '1.0 0 150.0 -3.0 15.0'
expression = 'e_h*max(min((t/end_t*(ymax-ymin)+ymin-y)/closure_dist,1),0)'
[../]
[./excav_downwards]
type = ParsedFunction
symbol_names = 'end_t ymin ymax e_h closure_dist'
symbol_values = '1.0 0 150.0 -3.0 15.0'
expression = 'e_h*t/end_t*max(min(((ymax-ymin)+ymin-y)/closure_dist,1),0)'
[../]
[]
[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]
type = ComputeLayeredCosseratElasticityTensor
young = 8E3 # MPa
poisson = 0.25
layer_thickness = 1.0
joint_normal_stiffness = 1E9 # huge
joint_shear_stiffness = 1E3
[../]
[./strain]
type = ComputeCosseratIncrementalSmallStrain
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = 'ini_xx 0 0 0 ini_xx 0 0 0 ini_zz'
eigenstrain_name = ini_stress
[../]
[./stress]
type = ComputeMultipleInelasticCosseratStress
block = 0
inelastic_models = mc
relative_tolerance = 2.0
absolute_tolerance = 1E6
max_iterations = 1
tangent_operator = nonlinear
perform_finite_strain_rotations = false
[../]
[./mc]
type = CappedMohrCoulombCosseratStressUpdate
block = 0
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
block = 0
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]
type = GenericConstantMaterial
prop_names = density
prop_values = 2500
[../]
[]
[Postprocessors]
[./subsidence]
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 = 1000
start_time = 0.0
dt = 0.2
end_time = 0.2
[]
[Outputs]
file_base = cosserat_mc_only
time_step_interval = 1
print_linear_residuals = false
csv = true
exodus = true
[./console]
type = Console
output_linear = false
[../]
[]
(modules/combined/test/tests/gravity/gravity_hex20.i)
# Gravity Test
#
# This test is designed to exercise the gravity body force kernel.
#
# The mesh for this problem is a rectangular bar 10 units by 1 unit
# by 1 unit.
#
# The boundary conditions for this problem are as follows. The
# displacement is zero on each of side that faces a negative
# coordinate direction. The acceleration of gravity is 20.
#
# The material has a Young's modulus of 1e6 and a density of 2.
#
# The analytic solution for the displacement along the bar is:
#
# u(x) = -b*x^2/(2*E)+b*L*x/E
#
# The displacement at x=L is b*L^2/(2*E) = 2*20*10*10/(2*1e6) = 0.002.
#
# The analytic solution for the stress along the bar assuming linear
# elasticity is:
#
# S(x) = b*(L-x)
#
# The stress at x=0 is b*L = 2*20*10 = 400.
#
# Note: The simulation does not measure stress at x=0. The stress
# is reported at element centers. The element closest to x=0 sits
# at x = 1/4 and has a stress of 390. This matches the linear
# stress distribution that is expected. The same situation applies
# at x = L where the stress is zero analytically. The nearest
# element is at x=9.75 where the stress is 10.
#
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
order = SECOND
family = LAGRANGE
[]
[Mesh]
file = gravity_hex20_test.e
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Modules/TensorMechanics/Master/All]
strain = FINITE
add_variables = true
generate_output = 'stress_xx'
[]
[Kernels]
[./gravity]
type = Gravity
variable = disp_x
value = 20
[../]
[]
[BCs]
[./no_x]
type = DirichletBC
variable = disp_x
boundary = 1
value = 0.0
[../]
[./no_y]
type = DirichletBC
variable = disp_y
boundary = 3
value = 0.0
[../]
[./no_z]
type = DirichletBC
variable = disp_z
boundary = 5
value = 0.0
[../]
[]
[Materials]
[./elasticty_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 1e6
bulk_modulus = 0.333333333333333e6
[../]
[./stress]
type = ComputeFiniteStrainElasticStress
[../]
[./density]
type = Density
density = 2
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
start_time = 0.0
end_time = 1.0
[./Quadrature]
order = THIRD
[../]
[]
[Outputs]
file_base = gravity_hex20_out
[./exodus]
type = Exodus
elemental_as_nodal = true
[../]
[]
(modules/solid_mechanics/examples/coal_mining/cosserat_mc_wp.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 300m deep
# and just the roof is studied (0<=z<=300). 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). Mining is simulated by moving the excavation's
# roof down, until disp_z=-3 at t=1.
# 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 are:
# - disp_x = 0 everywhere
# - disp_y = 0 at y=0 and y=450
# - disp_z = 0 for y>150
# - disp_z = -3 at maximum, for 0<=y<=150. See excav function.
# That is, rollers on the sides, free at top, and prescribed at bottom.
#
# 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 = 400.0
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
new_boundary = 11
normal = '0 -1 0'
input = generated_mesh
[]
[right]
type = SideSetsAroundSubdomainGenerator
new_boundary = 12
normal = '0 1 0'
input = left
[]
[front]
type = SideSetsAroundSubdomainGenerator
new_boundary = 13
normal = '-1 0 0'
input = right
[]
[back]
type = SideSetsAroundSubdomainGenerator
new_boundary = 14
normal = '1 0 0'
input = front
[]
[top]
type = SideSetsAroundSubdomainGenerator
new_boundary = 15
normal = '0 0 1'
input = back
[]
[bottom]
type = SideSetsAroundSubdomainGenerator
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 = SideSetsBetweenSubdomainsGenerator
new_boundary = 21
primary_block = 0
paired_block = 1
input = excav
[]
[hole]
type = BlockDeletionGenerator
block = 1
input = roof
[]
[]
[GlobalParams]
block = 0
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
[../]
[]
[AuxVariables]
[./disp_x]
[../]
[./wc_y]
[../]
[./wc_z]
[../]
[./stress_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yy]
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_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
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 16 21' # note addition of 16 and 21
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 = FunctionDirichletBC
variable = disp_z
boundary = 21
function = excav_sideways
[../]
[]
[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 e_h closure_dist'
symbol_values = '1.0 0 150.0 -3.0 15.0'
expression = 'e_h*max(min((min(t/end_t,1)*(ymax-ymin)+ymin-y)/closure_dist,1),0)'
[../]
[./excav_downwards]
type = ParsedFunction
symbol_names = 'end_t ymin ymax e_h closure_dist'
symbol_values = '1.0 0 150.0 -3.0 15.0'
expression = 'e_h*min(t/end_t,1)*max(min(((ymax-ymin)+ymin-y)/closure_dist,1),0)'
[../]
[]
[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]
type = ComputeLayeredCosseratElasticityTensor
young = 8E3 # MPa
poisson = 0.25
layer_thickness = 1.0
joint_normal_stiffness = 1E9 # huge
joint_shear_stiffness = 1E3
[../]
[./strain]
type = ComputeCosseratIncrementalSmallStrain
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = 'ini_xx 0 0 0 ini_xx 0 0 0 ini_zz'
eigenstrain_name = ini_stress
[../]
[./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
[../]
[./mc]
type = CappedMohrCoulombCosseratStressUpdate
block = 0
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 = 10000
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
block = 0
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.0
[../]
[./density]
type = GenericConstantMaterial
prop_names = density
prop_values = 2500
[../]
[]
[Postprocessors]
[./subsidence]
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 = 1000
start_time = 0.0
dt = 0.2
end_time = 0.2
[]
[Outputs]
file_base = cosserat_mc_wp
time_step_interval = 1
print_linear_residuals = false
csv = true
exodus = true
[./console]
type = Console
output_linear = false
[../]
[]
(modules/solid_mechanics/examples/coal_mining/cosserat_wp_only.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 300m deep
# and just the roof is studied (0<=z<=300). 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). Mining is simulated by moving the excavation's
# roof down, until disp_z=-3 at t=1.
# 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 are:
# - disp_x = 0 everywhere
# - disp_y = 0 at y=0 and y=450
# - disp_z = 0 for y>150
# - disp_z = -3 at maximum, for 0<=y<=150. See excav function.
# That is, rollers on the sides, free at top, and prescribed at bottom.
#
# 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.
#
# Below you will see Drucker-Prager parameters and AuxVariables, etc.
# These are not actally used in this example.
#
# 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
# Weak-plane cohesion = 0.1 MPa
# Weak-plane friction angle = 20 deg
# Weak-plane dilation angle = 10 deg
# Weak-plane tensile strength = 0.1 MPa
# Weak-plane compressive strength = 100 MPa, varying down to 1 MPa when tensile strain = 1
#
[Mesh]
[generated_mesh]
type = GeneratedMeshGenerator
dim = 3
nx = 1
xmin = -5
xmax = 5
nz = 40
zmin = 0
zmax = 400
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
new_boundary = 11
normal = '0 -1 0'
input = generated_mesh
[]
[right]
type = SideSetsAroundSubdomainGenerator
new_boundary = 12
normal = '0 1 0'
input = left
[]
[front]
type = SideSetsAroundSubdomainGenerator
new_boundary = 13
normal = '-1 0 0'
input = right
[]
[back]
type = SideSetsAroundSubdomainGenerator
new_boundary = 14
normal = '1 0 0'
input = front
[]
[top]
type = SideSetsAroundSubdomainGenerator
new_boundary = 15
normal = '0 0 1'
input = back
[]
[bottom]
type = SideSetsAroundSubdomainGenerator
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 = SideSetsBetweenSubdomainsGenerator
new_boundary = 21
primary_block = 0
paired_block = 1
input = excav
[]
[hole]
type = BlockDeletionGenerator
block = 1
input = roof
[]
[]
[GlobalParams]
block = 0
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
[../]
[]
[AuxVariables]
[./disp_x]
[../]
[./wc_y]
[../]
[./wc_z]
[../]
[./stress_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./dp_shear]
order = CONSTANT
family = MONOMIAL
[../]
[./dp_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
[../]
[./dp_shear_f]
order = CONSTANT
family = MONOMIAL
[../]
[./dp_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_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[../]
[./stress_zz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zz
index_i = 2
index_j = 2
[../]
[./dp_shear]
type = MaterialStdVectorAux
index = 0
property = dp_plastic_internal_parameter
variable = dp_shear
[../]
[./dp_tensile]
type = MaterialStdVectorAux
index = 1
property = dp_plastic_internal_parameter
variable = dp_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
[../]
[./dp_shear_f]
type = MaterialStdVectorAux
index = 0
property = dp_plastic_yield_function
variable = dp_shear_f
[../]
[./dp_tensile_f]
type = MaterialStdVectorAux
index = 1
property = dp_plastic_yield_function
variable = dp_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 16 21' # note addition of 16 and 21
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 = FunctionDirichletBC
variable = disp_z
boundary = 21
function = excav_sideways
[../]
[]
[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 e_h closure_dist'
symbol_values = '1.0 0 150.0 -3.0 15.0'
expression = 'e_h*max(min((t/end_t*(ymax-ymin)+ymin-y)/closure_dist,1),0)'
[../]
[./excav_downwards]
type = ParsedFunction
symbol_names = 'end_t ymin ymax e_h closure_dist'
symbol_values = '1.0 0 150.0 -3.0 15.0'
expression = 'e_h*t/end_t*max(min(((ymax-ymin)+ymin-y)/closure_dist,1),0)'
[../]
[]
[UserObjects]
[./dp_coh_strong_harden]
type = SolidMechanicsHardeningExponential
value_0 = 2.9 # MPa
value_residual = 3.1 # MPa
rate = 1.0
[../]
[./dp_fric]
type = SolidMechanicsHardeningConstant
value = 0.65 # 37deg
[../]
[./dp_dil]
type = SolidMechanicsHardeningConstant
value = 0.65
[../]
[./dp_tensile_str_strong_harden]
type = SolidMechanicsHardeningExponential
value_0 = 1.0 # MPa
value_residual = 1.4 # MPa
rate = 1.0
[../]
[./dp_compressive_str]
type = SolidMechanicsHardeningConstant
value = 1.0E3 # Large!
[../]
[./drucker_prager_model]
type = SolidMechanicsPlasticDruckerPrager
mc_cohesion = dp_coh_strong_harden
mc_friction_angle = dp_fric
mc_dilation_angle = dp_dil
internal_constraint_tolerance = 1 # irrelevant here
yield_function_tolerance = 1 # irrelevant here
[../]
[./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]
type = ComputeLayeredCosseratElasticityTensor
young = 8E3 # MPa
poisson = 0.25
layer_thickness = 1.0
joint_normal_stiffness = 1E9 # huge
joint_shear_stiffness = 1E3
[../]
[./strain]
type = ComputeCosseratIncrementalSmallStrain
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = 'ini_xx 0 0 0 ini_xx 0 0 0 ini_zz'
eigenstrain_name = ini_stress
[../]
[./stress]
type = ComputeMultipleInelasticCosseratStress
block = 0
inelastic_models = 'wp'
relative_tolerance = 2.0
absolute_tolerance = 1E6
max_iterations = 1
tangent_operator = nonlinear
perform_finite_strain_rotations = false
[../]
[./dp]
type = CappedDruckerPragerCosseratStressUpdate
block = 0
warn_about_precision_loss = false
host_youngs_modulus = 8E3
host_poissons_ratio = 0.25
base_name = dp
DP_model = drucker_prager_model
tensile_strength = dp_tensile_str_strong_harden
compressive_strength = dp_compressive_str
max_NR_iterations = 100000
tip_smoother = 0.1E1
smoothing_tol = 0.1E1 # MPa # Must be linked to cohesion
yield_function_tol = 1E-11 # MPa. this is essentially the lowest possible without lots of precision loss
perfect_guess = true
min_step_size = 1.0
[../]
[./wp]
type = CappedWeakPlaneCosseratStressUpdate
block = 0
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]
type = GenericConstantMaterial
prop_names = density
prop_values = 2500
[../]
[]
[Postprocessors]
[./subsidence]
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 = 1000
start_time = 0.0
dt = 0.2
end_time = 0.2
[]
[Outputs]
file_base = cosserat_wp_only
time_step_interval = 1
print_linear_residuals = false
csv = true
exodus = true
[./console]
type = Console
output_linear = false
[../]
[]
(modules/solid_mechanics/test/tests/gravity/ad_gravity_test.i)
#
# Gravity Test
#
# This test is designed to apply a gravity body force.
#
# The mesh is composed of one block with a single element.
# The bottom is fixed in all three directions. Poisson's ratio
# is zero and the density is 20/9.81
# which makes it trivial to check displacements.
#
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Mesh]
type = GeneratedMesh
dim = 3
[]
[Physics/SolidMechanics/QuasiStatic]
[all]
add_variables = true
use_automatic_differentiation = true
[]
[]
[Kernels]
[gravity_y]
type = ADGravity
variable = disp_y
value = -9.81
[]
[]
[BCs]
[no_x]
type = DirichletBC
variable = disp_x
boundary = bottom
value = 0.0
[]
[no_y]
type = DirichletBC
variable = disp_y
boundary = bottom
value = 0.0
[]
[no_z]
type = DirichletBC
variable = disp_z
boundary = bottom
value = 0.0
[]
[]
[Materials]
[Elasticity_tensor]
type = ADComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0 0.5e6'
[]
[stress]
type = ADComputeLinearElasticStress
[]
[density]
type = ADGenericConstantMaterial
prop_names = density
prop_values = 2.0387
[]
[]
[Executioner]
type = Steady
solve_type = NEWTON
nl_abs_tol = 1e-10
l_max_its = 20
[]
[Outputs]
[out]
type = Exodus
elemental_as_nodal = true
[]
[]
(modules/contact/test/tests/explicit_dynamics/settlement.i)
# One element test to test the central difference time integrator in 3D.
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
volumetric_locking_correction = true
[]
[Mesh]
[block_one]
type = GeneratedMeshGenerator
dim = 3
nx = 3
ny = 3
nz = 3
xmin = 4.5
xmax = 5.5
ymin = 4.5
ymax = 5.5
zmin = 0.0001
zmax = 1.0001
boundary_name_prefix = 'ball'
[]
[block_two]
type = GeneratedMeshGenerator
dim = 3
nx = 2
ny = 2
nz = 2
xmin = 0.0
xmax = 10
ymin = 0.0
ymax = 10
zmin = -2
zmax = 0
boundary_name_prefix = 'base'
boundary_id_offset = 10
[]
[block_one_id]
type = SubdomainIDGenerator
input = block_one
subdomain_id = 1
[]
[block_two_id]
type = SubdomainIDGenerator
input = block_two
subdomain_id = 2
[]
[combine]
type = MeshCollectionGenerator
inputs = ' block_one_id block_two_id'
[]
[]
[AuxVariables]
[penetration]
[]
[]
[AuxKernels]
[penetration]
type = PenetrationAux
variable = penetration
boundary = ball_back
paired_boundary = base_front
quantity = distance
[]
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[]
[AuxVariables]
[gap_rate]
[]
[vel_x]
[]
[accel_x]
[]
[vel_y]
[]
[accel_y]
[]
[vel_z]
[]
[accel_z]
[]
[stress_zz]
family = MONOMIAL
order = CONSTANT
[]
[strain_zz]
family = MONOMIAL
order = CONSTANT
[]
[kinetic_energy_one]
order = CONSTANT
family = MONOMIAL
[]
[elastic_energy_one]
order = CONSTANT
family = MONOMIAL
[]
[kinetic_energy_two]
order = CONSTANT
family = MONOMIAL
[]
[elastic_energy_two]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[stress_zz]
type = RankTwoAux
rank_two_tensor = stress
index_i = 2
index_j = 2
variable = stress_zz
execute_on = 'TIMESTEP_END'
[]
[strain_zz]
type = RankTwoAux
rank_two_tensor = mechanical_strain
index_i = 2
index_j = 2
variable = strain_zz
[]
[accel_x]
type = TestNewmarkTI
variable = accel_x
displacement = disp_x
first = false
execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
[]
[vel_x]
type = TestNewmarkTI
variable = vel_x
displacement = disp_x
execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
[]
[accel_y]
type = TestNewmarkTI
variable = accel_y
displacement = disp_y
first = false
execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
[]
[vel_y]
type = TestNewmarkTI
variable = vel_y
displacement = disp_x
execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
[]
[accel_z]
type = TestNewmarkTI
variable = accel_z
displacement = disp_z
first = false
execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
[]
[vel_z]
type = TestNewmarkTI
variable = vel_z
displacement = disp_z
execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
[]
[kinetic_energy_one]
type = KineticEnergyAux
block = '1'
variable = kinetic_energy_one
newmark_velocity_x = vel_x
newmark_velocity_y = vel_y
newmark_velocity_z = vel_z
density = density
[]
[elastic_energy_one]
type = ElasticEnergyAux
variable = elastic_energy_one
block = '1'
[]
[kinetic_energy_two]
type = KineticEnergyAux
block = '2'
variable = kinetic_energy_two
newmark_velocity_x = vel_x
newmark_velocity_y = vel_y
newmark_velocity_z = vel_z
density = density
[]
[elastic_energy_two]
type = ElasticEnergyAux
variable = elastic_energy_two
block = '2'
[]
[]
[Kernels]
[DynamicTensorMechanics]
displacements = 'disp_x disp_y disp_z'
stiffness_damping_coefficient = 1.0e-3
generate_output = 'stress_zz strain_zz'
[]
[inertia_x]
type = InertialForce
variable = disp_x
[]
[inertia_y]
type = InertialForce
variable = disp_y
[]
[inertia_z]
type = InertialForce
variable = disp_z
[]
[gravity]
type = Gravity
variable = disp_z
value = -98.10
[]
[]
[BCs]
[x_front]
type = DirichletBC
variable = disp_x
boundary = 'ball_front'
preset = false
value = 0.0
[]
[y_front]
type = DirichletBC
variable = disp_y
boundary = 'ball_front'
preset = false
value = 0.0
[]
[x_fixed]
type = DirichletBC
variable = disp_x
boundary = 'base_back'
preset = false
value = 0.0
[]
[y_fixed]
type = DirichletBC
variable = disp_y
boundary = 'base_back'
preset = false
value = 0.0
[]
[z_fixed]
type = DirichletBC
variable = disp_z
boundary = 'base_back'
preset = false
value = 0.0
[]
[z_fixed_front]
type = DirichletBC
variable = disp_z
boundary = 'base_front'
preset = false
value = 0.0
[]
[]
[ExplicitDynamicsContact]
[my_contact]
model = frictionless_balance
primary = base_front
secondary = ball_back
vel_x = 'vel_x'
vel_y = 'vel_y'
vel_z = 'vel_z'
[]
[]
[Materials]
[elasticity_tensor_block_one]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 1e6
poissons_ratio = 0.0
block = 1
outputs = 'exodus'
output_properties = __all__
[]
[elasticity_tensor_block_two]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 1e10
poissons_ratio = 0.0
block = 2
outputs = 'exodus'
output_properties = __all__
[]
[strain_block]
type = ComputeFiniteStrain # ComputeIncrementalSmallStrain
displacements = 'disp_x disp_y disp_z'
implicit = false
[]
[stress_block]
type = ComputeFiniteStrainElasticStress
[]
[density_one]
type = GenericConstantMaterial
prop_names = density
prop_values = 1e1
outputs = 'exodus'
output_properties = 'density'
block = '1'
[]
[density_two]
type = GenericConstantMaterial
prop_names = density
prop_values = 1e6
outputs = 'exodus'
output_properties = 'density'
block = '2'
[]
[wave_speed]
type = WaveSpeed
outputs = 'exodus'
output_properties = 'wave_speed'
[]
[]
[Executioner]
type = Transient
start_time = -0.01
end_time = 0.04
dt = 1.0e-4
timestep_tolerance = 1e-6
[TimeIntegrator]
type = CentralDifference
solve_type = lumped
[]
[]
[Outputs]
interval = 1
exodus = true
csv = true
execute_on = 'TIMESTEP_END'
[]
[Postprocessors]
[accel_58z]
type = NodalVariableValue
nodeid = 1
variable = accel_z
[]
[vel_58z]
type = NodalVariableValue
nodeid = 1
variable = vel_z
[]
[critical_time_step]
type = CriticalTimeStep
[]
[contact_pressure_max]
type = NodalExtremeValue
variable = contact_pressure
block = '1 2'
value_type = max
[]
[penetration_max]
type = NodalExtremeValue
variable = penetration
block = '1 2'
value_type = max
[]
[total_kinetic_energy_one]
type = ElementIntegralVariablePostprocessor
variable = kinetic_energy_one
block = '1'
[]
[total_elastic_energy_one]
type = ElementIntegralVariablePostprocessor
variable = elastic_energy_one
block = '1'
[]
[total_kinetic_energy_two]
type = ElementIntegralVariablePostprocessor
variable = kinetic_energy_two
block = '2'
[]
[total_elastic_energy_two]
type = ElementIntegralVariablePostprocessor
variable = elastic_energy_two
block = '2'
[]
[]
(modules/solid_mechanics/test/tests/dynamics/prescribed_displacement/3D_QStatic_1_Ramped_Displacement_with_gravity.i)
# One 3D element under ramped displacement loading.
#
# loading in z direction:
# time : 0.0 0.1 0.2 0.3
# disp : 0.0 0.0 -0.01 -0.01
# Gravity is applied in y direction. To equilibrate the system
# under gravity, a static analysis is run in the first time step
# by turning off the inertial terms. (see controls block and
# DynamicSolidMechanics block).
# Result: The displacement at the top node in the z direction should match
# the prescribed displacement. Also, the z acceleration should
# be two triangular pulses, one peaking at 0.1 and another peaking at
# 0.2.
# The y displacement would be offset by the gravity displacement.
# Also the y acceleration and velocity should be zero until the loading in
# the z direction starts (i.e, until 0.1s)
# Note: The time step used in the displacement data file should match
# the simulation time step (dt and dtmin in the Executioner block).
[Mesh]
type = GeneratedMesh
dim = 3 # Dimension of the mesh
nx = 1 # Number of elements in the x direction
ny = 1 # Number of elements in the y direction
nz = 1 # Number of elements in the z direction
xmin = 0.0
xmax = 1
ymin = 0.0
ymax = 1
zmin = 0.0
zmax = 1
allow_renumbering = false # So NodalVariableValue can index by id
[]
[Variables] # variables that are solved
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[AuxVariables] # variables that are calculated for output
[./accel_x]
[../]
[./vel_x]
[../]
[./accel_y]
[../]
[./vel_y]
[../]
[./accel_z]
[../]
[./vel_z]
[../]
[./stress_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_zz]
order = CONSTANT
family = MONOMIAL
[../]
[]
[Kernels]
[./DynamicSolidMechanics] # zeta*K*vel + K * disp
displacements = 'disp_x disp_y disp_z'
stiffness_damping_coefficient = 0.000025
static_initialization = true #turns off rayliegh damping for the first time step to stabilize system under gravity
[../]
[./inertia_x] # M*accel + eta*M*vel
type = InertialForce
variable = disp_x
velocity = vel_x
acceleration = accel_x
beta = 0.25 # Newmark time integration
gamma = 0.5 # Newmark time integration
eta = 19.63
[../]
[./inertia_y]
type = InertialForce
variable = disp_y
velocity = vel_y
acceleration = accel_y
beta = 0.25
gamma = 0.5
eta = 19.63
[../]
[./inertia_z]
type = InertialForce
variable = disp_z
velocity = vel_z
acceleration = accel_z
beta = 0.25
gamma = 0.5
eta = 19.63
[../]
[./gravity]
type = Gravity
variable = disp_y
value = -9.81
[../]
[]
[AuxKernels]
[./accel_x] # Calculates and stores acceleration at the end of time step
type = NewmarkAccelAux
variable = accel_x
displacement = disp_x
velocity = vel_x
beta = 0.25
execute_on = timestep_end
[../]
[./vel_x] # Calculates and stores velocity at the end of the time step
type = NewmarkVelAux
variable = vel_x
acceleration = accel_x
gamma = 0.5
execute_on = timestep_end
[../]
[./accel_y]
type = NewmarkAccelAux
variable = accel_y
displacement = disp_y
velocity = vel_y
beta = 0.25
execute_on = timestep_end
[../]
[./vel_y]
type = NewmarkVelAux
variable = vel_y
acceleration = accel_y
gamma = 0.5
execute_on = timestep_end
[../]
[./accel_z]
type = NewmarkAccelAux
variable = accel_z
displacement = disp_z
velocity = vel_z
beta = 0.25
execute_on = timestep_end
[../]
[./vel_z]
type = NewmarkVelAux
variable = vel_z
acceleration = accel_z
gamma = 0.5
execute_on = timestep_end
[../]
[./stress_xx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx
index_i = 0
index_j = 0
[../]
[./strain_xx]
type = RankTwoAux
rank_two_tensor = total_strain
variable = strain_xx
index_i = 0
index_j = 0
[../]
[./stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[../]
[./strain_yy]
type = RankTwoAux
rank_two_tensor = total_strain
variable = strain_yy
index_i = 1
index_j = 1
[../]
[./stress_zz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zz
index_i = 2
index_j = 2
[../]
[./strain_zz]
type = RankTwoAux
rank_two_tensor = total_strain
variable = strain_zz
index_i = 2
index_j = 2
[../]
[]
[Functions]
[./displacement_front]
type = PiecewiseLinear
data_file = 'displacement.csv'
format = columns
[../]
[]
[BCs]
[./prescribed_displacement]
type = PresetDisplacement
variable = disp_z
velocity = vel_z
acceleration = accel_z
beta = 0.25
boundary = front
function = displacement_front
[../]
[./anchor_x]
type = DirichletBC
variable = disp_x
boundary = left
value = 0.0
[../]
[./anchor_y]
type = DirichletBC
variable = disp_y
boundary = bottom
value = 0.0
[../]
[./anchor_z]
type = DirichletBC
variable = disp_z
boundary = back
value = 0.0
[../]
[]
[Materials]
[./elasticity_tensor]
youngs_modulus = 325e6 #Pa
poissons_ratio = 0.3
type = ComputeIsotropicElasticityTensor
block = 0
[../]
[./strain]
#Computes the strain, assuming small strains
type = ComputeSmallStrain
block = 0
displacements = 'disp_x disp_y disp_z'
[../]
[./stress]
#Computes the stress, using linear elasticity
type = ComputeLinearElasticStress
block = 0
[../]
[./density]
type = GenericConstantMaterial
block = 0
prop_names = density
prop_values = 2000 #kg/m3
[../]
[]
[Controls] # turns off inertial terms for the first time step
[./period0]
type = TimePeriod
disable_objects = '*/vel_x */vel_y */vel_z */accel_x */accel_y */accel_z */inertia_x */inertia_y */inertia_z'
start_time = 0.0
end_time = 0.1 # dt used in the simulation
[../]
[../]
[Executioner]
type = Transient
start_time = 0
end_time = 3.0
l_tol = 1e-6
nl_rel_tol = 1e-6
nl_abs_tol = 1e-6
dt = 0.1
timestep_tolerance = 1e-6
[]
[Postprocessors] # These quantites are printed to a csv file at every time step
[./_dt]
type = TimestepSize
[../]
[./accel_6x]
type = NodalVariableValue
nodeid = 6
variable = accel_x
[../]
[./accel_6y]
type = NodalVariableValue
nodeid = 6
variable = accel_y
[../]
[./accel_6z]
type = NodalVariableValue
nodeid = 6
variable = accel_z
[../]
[./vel_6x]
type = NodalVariableValue
nodeid = 6
variable = vel_x
[../]
[./vel_6y]
type = NodalVariableValue
nodeid = 6
variable = vel_y
[../]
[./vel_6z]
type = NodalVariableValue
nodeid = 6
variable = vel_z
[../]
[./disp_6x]
type = NodalVariableValue
nodeid = 6
variable = disp_x
[../]
[./disp_6y]
type = NodalVariableValue
nodeid = 6
variable = disp_y
[../]
[./disp_6z]
type = NodalVariableValue
nodeid = 6
variable = disp_z
[../]
[]
[Outputs]
exodus = true
perf_graph = true
[]
(modules/solid_mechanics/test/tests/gravity/gravity_test.i)
#
# Gravity Test
#
# This test is designed to apply a gravity body force.
#
# The mesh is composed of one block with a single element.
# The bottom is fixed in all three directions. Poisson's ratio
# is zero and the density is 20/9.81
# which makes it trivial to check displacements.
#
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Mesh]
type = GeneratedMesh
dim = 3
[]
[Physics/SolidMechanics/QuasiStatic]
[all]
add_variables = true
[]
[]
[Kernels]
[gravity_y]
type = Gravity
variable = disp_y
value = -9.81
[]
[]
[BCs]
[no_x]
type = DirichletBC
variable = disp_x
boundary = bottom
value = 0.0
[]
[no_y]
type = DirichletBC
variable = disp_y
boundary = bottom
value = 0.0
[]
[no_z]
type = DirichletBC
variable = disp_z
boundary = bottom
value = 0.0
[]
[]
[Materials]
[Elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0 0.5e6'
[]
[stress]
type = ComputeLinearElasticStress
[]
[density]
type = GenericConstantMaterial
prop_names = density
prop_values = 2.0387
[]
[]
[Executioner]
type = Steady
nl_abs_tol = 1e-10
l_max_its = 20
[]
[Outputs]
[out]
type = Exodus
elemental_as_nodal = true
[]
[]
(modules/solid_mechanics/examples/bridge/bridge.i)
#
# Bridge linear elasticity example
#
# This example models a bridge using linear elasticity.
# It can be either steel or concrete.
# Gravity is applied
# A pressure of 0.5 MPa is also applied
#
[Mesh]
displacements = 'disp_x disp_y disp_z' #Define displacements for deformed mesh
type = FileMesh #Read in mesh from file
file = bridge.e
boundary_id = '1 2 3 4 5 6' #Assign names to boundaries to make things clearer
boundary_name = 'top left right bottom1 bottom2 bottom3'
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Variables]
[./disp_x]
order = FIRST
family = LAGRANGE
[../]
[./disp_y]
order = FIRST
family = LAGRANGE
[../]
[./disp_z]
order = FIRST
family = LAGRANGE
[../]
[]
[Kernels]
[./gravity_y]
#Gravity is applied to bridge
type = Gravity
variable = disp_y
value = -9.81
[../]
[./TensorMechanics]
#Stress divergence kernels
displacements = 'disp_x disp_y disp_z'
[../]
[]
[AuxVariables]
[./von_mises]
#Dependent variable used to visualize the Von Mises stress
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./von_mises_kernel]
#Calculates the von mises stress and assigns it to von_mises
type = RankTwoScalarAux
variable = von_mises
rank_two_tensor = stress
execute_on = timestep_end
scalar_type = VonMisesStress
[../]
[]
[BCs]
[./Pressure]
[./load]
#Applies the pressure
boundary = top
factor = 5e5 # Pa
[../]
[../]
[./anchor_x]
#Anchors the bottom and sides against deformation in the x-direction
type = DirichletBC
variable = disp_x
boundary = 'left right bottom1 bottom2 bottom3'
value = 0.0
[../]
[./anchor_y]
#Anchors the bottom and sides against deformation in the y-direction
type = DirichletBC
variable = disp_y
boundary = 'left right bottom1 bottom2 bottom3'
value = 0.0
[../]
[./anchor_z]
#Anchors the bottom and sides against deformation in the z-direction
type = DirichletBC
variable = disp_z
boundary = 'left right bottom1 bottom2 bottom3'
value = 0.0
[../]
[]
[Materials]
active = 'density_concrete stress strain elasticity_tensor_concrete'
[./elasticity_tensor_steel]
#Creates the elasticity tensor using steel parameters
youngs_modulus = 210e9 #Pa
poissons_ratio = 0.3
type = ComputeIsotropicElasticityTensor
block = 1
[../]
[./elasticity_tensor_concrete]
#Creates the elasticity tensor using concrete parameters
youngs_modulus = 16.5e9 #Pa
poissons_ratio = 0.2
type = ComputeIsotropicElasticityTensor
block = 1
[../]
[./strain]
#Computes the strain, assuming small strains
type = ComputeSmallStrain
block = 1
displacements = 'disp_x disp_y disp_z'
[../]
[./stress]
#Computes the stress, using linear elasticity
type = ComputeLinearElasticStress
block = 1
[../]
[./density_steel]
#Defines the density of steel
type = GenericConstantMaterial
block = 1
prop_names = density
prop_values = 7850 # kg/m^3
[../]
[./density_concrete]
#Defines the density of concrete
type = GenericConstantMaterial
block = 1
prop_names = density
prop_values = 2400 # kg/m^3
[../]
[]
[Preconditioning]
[./SMP]
#Creates the entire Jacobian, for the Newton solve
type = SMP
full = true
[../]
[]
[Executioner]
#We solve a steady state problem using Newton's iteration
type = Steady
solve_type = NEWTON
nl_rel_tol = 1e-9
l_max_its = 30
l_tol = 1e-4
nl_max_its = 10
petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
petsc_options_value = 'hypre boomeramg 31'
[]
[Outputs]
exodus = true
perf_graph = true
[]
(modules/solid_mechanics/examples/bridge/bridge_large_strain.i)
#
# Bridge linear elasticity example
#
# This example models a bridge using linear elasticity.
# It can be either steel or concrete.
# Gravity is applied
# A pressure of 0.5 MPa is also applied
#
[Mesh]
displacements = 'disp_x disp_y disp_z' #Define displacements for deformed mesh
type = FileMesh #Read in mesh from file
file = bridge.e
boundary_id = '1 2 3 4 5 6' #Assign names to boundaries to make things clearer
boundary_name = 'top left right bottom1 bottom2 bottom3'
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Variables]
[./disp_x]
order = FIRST
family = LAGRANGE
[../]
[./disp_y]
order = FIRST
family = LAGRANGE
[../]
[./disp_z]
order = FIRST
family = LAGRANGE
[../]
[]
[Kernels]
[./gravity_y]
#Gravity is applied to bridge
type = Gravity
variable = disp_y
value = -9.81
[../]
[./SolidMechanics]
#Stress divergence kernels
displacements = 'disp_x disp_y disp_z'
[../]
[]
[AuxVariables]
[./von_mises]
#Dependent variable used to visualize the Von Mises stress
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./von_mises_kernel]
#Calculates the von mises stress and assigns it to von_mises
type = RankTwoScalarAux
variable = von_mises
rank_two_tensor = stress
execute_on = timestep_end
scalar_type = VonMisesStress
[../]
[]
[BCs]
[./Pressure]
[./load]
#Applies the pressure
boundary = top
factor = 5e5 # Pa
[../]
[../]
[./anchor_x]
#Anchors the bottom and sides against deformation in the x-direction
type = DirichletBC
variable = disp_x
boundary = 'left right bottom1 bottom2 bottom3'
value = 0.0
[../]
[./anchor_y]
#Anchors the bottom and sides against deformation in the y-direction
type = DirichletBC
variable = disp_y
boundary = 'left right bottom1 bottom2 bottom3'
value = 0.0
[../]
[./anchor_z]
#Anchors the bottom and sides against deformation in the z-direction
type = DirichletBC
variable = disp_z
boundary = 'left right bottom1 bottom2 bottom3'
value = 0.0
[../]
[]
[Materials]
active = 'density_steel stress strain elasticity_tensor_steel'
[./elasticity_tensor_steel]
#Creates the elasticity tensor using steel parameters
youngs_modulus = 210e9 #Pa
poissons_ratio = 0.3
type = ComputeIsotropicElasticityTensor
block = 1
[../]
[./elasticity_tensor_concrete]
#Creates the elasticity tensor using concrete parameters
youngs_modulus = 16.5e9 #Pa
poissons_ratio = 0.2
type = ComputeIsotropicElasticityTensor
block = 1
[../]
[./strain]
#Computes the strain, assuming small strains
type = ComputeFiniteStrain
block = 1
displacements = 'disp_x disp_y disp_z'
[../]
[./stress]
#Computes the stress, using linear elasticity
type = ComputeFiniteStrainElasticStress
block = 1
[../]
[./density_steel]
#Defines the density of steel
type = GenericConstantMaterial
block = 1
prop_names = density
prop_values = 7850 # kg/m^3
[../]
[./density_concrete]
#Defines the density of concrete
type = GenericConstantMaterial
block = 1
prop_names = density
prop_values = 2400 # kg/m^3
[../]
[]
[Preconditioning]
[./SMP]
#Creates the entire Jacobian, for the Newton solve
type = SMP
full = true
[../]
[]
[Executioner]
#We solve a steady state problem using Newton's iteration
type = Transient
solve_type = NEWTON
nl_rel_tol = 1e-9
l_max_its = 30
l_tol = 1e-4
nl_max_its = 10
petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
petsc_options_value = 'hypre boomeramg 31'
dt = 0.1
num_steps = 1
[]
[Outputs]
exodus = true
perf_graph = true
[]
(modules/combined/test/tests/gravity/gravity_rz.i)
# Gravity Test
#
# This test is designed to exercise the gravity body force rz kernel.
#
# The mesh for this problem is a rectangle 10 units by 1 unit.
#
# The boundary conditions for this problem are as follows. The
# displacement is zero at the top. The acceleration of gravity is 20.
#
# The material has a Young's modulus of 1e6 and a density of 2.
#
# The analytic solution for the displacement along the bar is:
#
# u(y) = -b*y^2/(2*E)+b*L*y/E
#
# The displacement at y=L is b*L^2/(2*E) = 2*20*10*10/(2*1e6) = 0.002.
#
# The analytic solution for the stress along the bar assuming linear
# elasticity is:
#
# S(y) = b*(L-y)
#
# The stress at x=0 is b*L = 2*20*10 = 400.
#
# Note: The simulation does not measure stress at y=0. The stress
# is reported at element centers. The element closest to y=0 sits
# at y = 1/4 and has a stress of 390. This matches the linear
# stress distribution that is expected. The same situation applies
# at y = L where the stress is zero analytically. The nearest
# element is at y=9.75 where the stress is 10.
#
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Problem]
coord_type = RZ
[]
[Mesh]
file = gravity_rz_test.e
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[]
[Modules/TensorMechanics/Master/All]
volumetric_locking_correction = true
strain = FINITE
add_variables = true
generate_output = 'stress_xx stress_yy stress_xy'
[]
[Kernels]
[./gravity]
type = Gravity
variable = disp_y
value = 20
[../]
[]
[BCs]
[./no_y]
type = DirichletBC
variable = disp_y
boundary = 2
value = 0.0
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
shear_modulus = 0.5e6
lambda = 0.0
[../]
[./stress]
type = ComputeFiniteStrainElasticStress
[../]
[./density]
type = Density
density = 2
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
start_time = 0.0
end_time = 1.0
[]
[Outputs]
file_base = gravity_rz_out
[./exodus]
type = Exodus
elemental_as_nodal = true
[../]
[]
(modules/combined/test/tests/gravity/gravity_rz_quad8.i)
# Gravity Test
#
# This test is designed to exercise the gravity body force rz kernel.
#
# The mesh for this problem is a rectangle 10 units by 1 unit.
#
# The boundary conditions for this problem are as follows. The
# displacement is zero at the top. The acceleration of gravity is 20.
#
# The material has a Young's modulus of 1e6 and a density of 2.
#
# The analytic solution for the displacement along the bar is:
#
# u(y) = -b*y^2/(2*E)+b*L*y/E
#
# The displacement at y=L is b*L^2/(2*E) = 2*20*10*10/(2*1e6) = 0.002.
#
# The analytic solution for the stress along the bar assuming linear
# elasticity is:
#
# S(y) = b*(L-y)
#
# The stress at x=0 is b*L = 2*20*10 = 400.
#
# Note: The simulation does not measure stress at y=0. The stress
# is reported at element centers. The element closest to y=0 sits
# at y = 1/4 and has a stress of 390. This matches the linear
# stress distribution that is expected. The same situation applies
# at y = L where the stress is zero analytically. The nearest
# element is at y=9.75 where the stress is 10.
#
[GlobalParams]
displacements = 'disp_x disp_y'
order = SECOND
family = LAGRANGE
[]
[Problem]
coord_type = RZ
[]
[Mesh]
file = gravity_rz_quad8_test.e
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[]
[Modules/TensorMechanics/Master/All]
strain = FINITE
add_variables = true
generate_output = 'stress_xx stress_yy stress_xy'
[]
[Kernels]
[./gravity]
type = Gravity
variable = disp_y
value = 20
[../]
[]
[BCs]
[./no_y]
type = DirichletBC
variable = disp_y
boundary = 2
value = 0.0
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
shear_modulus = 0.5e6
lambda = 0.0
[../]
[./stress]
type = ComputeFiniteStrainElasticStress
[../]
[./density]
type = Density
density = 2
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
start_time = 0.0
end_time = 1.0
[./Quadrature]
order = THIRD
[../]
[]
[Outputs]
file_base = gravity_rz_quad8_out
[./exodus]
type = Exodus
elemental_as_nodal = true
[../]
[]