- cohesionA SolidMechanicsHardening UserObject that defines hardening of the cohesion. Physically the cohesion should not be negative.
C++ Type:UserObjectName
Unit:(no unit assumed)
Controllable:No
Description:A SolidMechanicsHardening UserObject that defines hardening of the cohesion. Physically the cohesion should not be negative.
- compressive_strengthA SolidMechanicsHardening UserObject that defines hardening of the weak-plane compressive strength. In physical situations this is positive.
C++ Type:UserObjectName
Unit:(no unit assumed)
Controllable:No
Description:A SolidMechanicsHardening UserObject that defines hardening of the weak-plane compressive strength. In physical situations this is positive.
- smoothing_tolIntersections of the yield surfaces will be smoothed by this amount (this is measured in units of stress). Often this is related to other physical parameters (eg, 0.1*cohesion) but it is important to set this small enough so that the individual yield surfaces do not mix together in the smoothing process to produce a result where no stress is admissible (for example, mixing together tensile and compressive failure envelopes).
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Intersections of the yield surfaces will be smoothed by this amount (this is measured in units of stress). Often this is related to other physical parameters (eg, 0.1*cohesion) but it is important to set this small enough so that the individual yield surfaces do not mix together in the smoothing process to produce a result where no stress is admissible (for example, mixing together tensile and compressive failure envelopes).
- tan_dilation_angleA SolidMechanicsHardening UserObject that defines hardening of the tan(dilation angle). Usually the dilation angle is not greater than the friction angle, and it is between 0 and 90deg.
C++ Type:UserObjectName
Unit:(no unit assumed)
Controllable:No
Description:A SolidMechanicsHardening UserObject that defines hardening of the tan(dilation angle). Usually the dilation angle is not greater than the friction angle, and it is between 0 and 90deg.
- tan_friction_angleA SolidMechanicsHardening UserObject that defines hardening of tan(friction angle). Physically the friction angle should be between 0 and 90deg.
C++ Type:UserObjectName
Unit:(no unit assumed)
Controllable:No
Description:A SolidMechanicsHardening UserObject that defines hardening of tan(friction angle). Physically the friction angle should be between 0 and 90deg.
- tensile_strengthA SolidMechanicsHardening UserObject that defines hardening of the weak-plane tensile strength. In physical situations this is positive (and always must be greater than negative compressive-strength.
C++ Type:UserObjectName
Unit:(no unit assumed)
Controllable:No
Description:A SolidMechanicsHardening UserObject that defines hardening of the weak-plane tensile strength. In physical situations this is positive (and always must be greater than negative compressive-strength.
- tip_smootherThe cone vertex at shear-stress = 0 will be smoothed by the given amount. Typical value is 0.1*cohesion
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The cone vertex at shear-stress = 0 will be smoothed by the given amount. Typical value is 0.1*cohesion
- yield_function_tolThe return-map process will be deemed to have converged if all yield functions are within yield_function_tol of zero. If this is set very low then precision-loss might be encountered: if the code detects precision loss then it also deems the return-map process has converged.
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The return-map process will be deemed to have converged if all yield functions are within yield_function_tol of zero. If this is set very low then precision-loss might be encountered: if the code detects precision loss then it also deems the return-map process has converged.
CappedWeakPlaneStressUpdate
Capped weak-plane plasticity stress calculator
Theory
Weak-plane plasticity is designed to simulate a layered material. Each layer can slip over adjacent layers, and be separated from those adjacent layers. An example of particular interest is stratified rocks, which consist of large sheets of fairly strong rock separated by weak and very thin joints. Upon application of stress, the joints can fail, either by slipping or separating. The idea is to use one finite element that potentially contains many layers, and prescribe weak plane plasticity'' for that finite element, so that it can fail by joint separation and joint slipping.
Denote the normal to the layers by , and the tangential directions by and . It is convenient to introduce two new stress variables in terms of the stress tensor :
(1)
In standard elasticity, the stress tensor is symmetric, so an equivalent definition of is , however the symmetrization is deliberately not written in Eq. (1) and below so that the equations also hold for the Cosserat case (see CappedWeakPlaneCosseratStressUpdate).
The joint slipping is assumed to be governed by a Drucker-Prager type of plasticity with a cohesion , and friction angle :
The parameter and may be constants, or they may harden or soften (more on this later).
Joints may also open, and this type of failure is assumed to be governed by a tensile failure yield function:
where is the tensile strength, which may be constant or harden or soften.
Joints may also close, and this type of failure is assumed to be governed by a compressive failure yield function:
where is the compressive strength (a positive quantity), which may be constant or harden or soften.
The yield functions and place caps'' on the shear yield function . The combined yield function is simply
(2)
which defines the admissible domain where all yield functions are non-positive, and the inadmissible domain where at least one yield function is positive.
One of the features of this plasticity is the ability to model cyclic behavior. For instance, the compressive strength may be initially very high. However, after tensile failure, the compressive strength can soften to zero in order to model the fact that the material now contains open joints which cannot support any compressive load. If the material then fails in compression (eg, because it gets squashed) and the joints close then the compressive strength can be made high again.
Flow rules and hardening
This plasticity is non-associative. Define the dilation angle , which may be constant, or harden or soften. The shear flow potential is
The tensile flow potential is
and the compressive flow potential is
The overall flow potential is
(3)
Obviously there are problems here where is not defined properly at the corners where and (or even ). This is resolved by using smoothing (more on this later).
This plasticity model contains two internal parameters, denote by and . It is assumed that
That is, can be thought of as the shear'' internal parameter, while is the tensile'' internal parameter.
To complete the definition of this plasticity model, the increments of and during the return-map process must be defined. The return-map process involves being provided with a trial stress and an existing value of the internal parameters , and finding a returned'' stress, , and internal parameters, , that satisfy
(4)
where is the elasticity tensor, and is a plastic multiplier'', that must be positive. The former expresses that the stress must be admissible, while the latter is called the normality condition''. Loosely speaking, the returned stress lies at a position on the yield surface where the normal points to the trial stress (actually and must be used to define the normal direction'').
Let us express the normality condition in space. The component is easy:
(5)
where the last equality holds by assumption (see full list of assumptions below). The and components are similar:
(6)
Another assumption has been made about . The final term is
This means that Eq. (6) can be re-written
A similar equation holds for the component, and these can be summed and rearranged to yield
(7)
Eq. (4), Eq. (5) and~Eq. (7) are the three conditions that need to be satisfied, and the three variables to be found are , and .
Consider the case of returning to the shear yield surface. Since and for this flow, the return-map process must solve the following system of equations
The solution satisfied and .
Now consider the case of returning to the tensile yield surface. The equations are
Comparing these two types of return, it is obvious that quantifies the amount of shear failure. Therefore, the following definitions are used in this plasticity model
The final term ensures that does not increase during pure shear failure. The scaling by ensures that these internal parameters are dimensionless.
In summary, this plasticity model is defined by the yield function of Eq. (2), the flow potential of Eq. (3), and the following return-map problem.
Return-map problem
At any given MOOSE timestep, given the old stress , and a total strain increment (that comes from the nonlinear solver proposing displacements) the trial stress is
This gives and . If , , and are such that , the return-map problem is: find , , , and such that
(8)
The latter two equations are assumed to hold in the smoothed situation (discussed below) too, and note that , so the these two equations are not completely trivial.
After the return-map problem has been solved, the stress components are , except for the following
The plastic strain is
The elastic strain is
Yield Smoothing
The shear yield function, , describes a cone in space. The cone's tip is problematic for the return-map process (the derivative is not defined there) and there are two main ways of getting around this. Firstly, a multi-surface technique can be used to define the return-map process. Secondly, the cone's tip can be smoothed. This plasticity model uses the second technique. The yield function is defined to be
and the flow potential is
The vertices where the shear yield surface meets the tensile and compressive yield surfaces also need to be handled. Smoothing is also used here. This uses a new type of smoothing. For the case at hand only two yield surfaces and flow potentials need to be smoothed (there are no points where three or more yield surfaces get close to each other) and only in 2D space, and a single parameter can be used. The parameter has the units of stress. At any point order the 3 yield function values, and denote the largest by , the second largest by and the smallest by :
Then the single, smoothed yield function is defined to be
The derivative of the flow potential is smoothed similarly.
Constraints and assumptions concerning parameters
The friction angle and cohesion should be positive, and the dilation angle should be non-negative. Furthermore, the MOOSE user must ensure that These conditions should be satisfied for all values of the internal parameter . MOOSE checks that these conditions hold for only.
The tensile and compressive strength must satisfy otherwise the caps'' are swapped and the assumption of a convex yield surface is violated. MOOSE checks this condition holds for only: the MOOSE user must ensure that it actually holds for all values of the internal parameter
The smoothing parameter must be chosen carefully. At no time should the tensile cap mix with the compressive cap via smoothing, otherwise this typically means that no stress is admissible and MOOSE will never converge. For instance, if , then a smoothing parameter of 0.1 is fine, but a smoothing parameter will cause mixing of tension with compression. The MOOSE user must ensure that this holds for all values of the internal parameters.
The tip-smoothing parameter is important, even if the tensile cap completely chops off the shear-cone's tip. This is because MOOSE can explore regions of parameter space where the cone's tip is exposed.
It is vital that the smoothing parameters and are chosen so that the yield surface is not wildly varying around , otherwise poor convergence of the return-map process will occur.
It is assumed that the elasticity tensor has the following symmetries:
(9)
and that
and that
and that
(10)
These are quite standard conditions that hold for all non-Cosserat materials to our knowledge.
Technical discussions
Unknowns and the convergence criterion
The return-map problem Eq. (8) is solved as a system consisting of the first 3 equations, and substituting the fourth and fifth equations wherever needed. The three unknowns are , and , which all have the same units. Convergence is deemed to be achieved when the sum of squares of the residuals of these 3 equations is less than a user-defined tolerance.
Iterative procedure and initial guesses
A Newton-Raphson process is used, along with a cubic line-search. The process may be initialized with the solution that is correct for perfect plasticity (no hardening) and no smoothing, if the user desires. Smoothing adds nonlinearities, so this initial guess will not always be the exact answer. For hardening, it is not always advantageous to initialize the Newton-Raphson process in this way, as the yield surfaces can move dramatically during the return process.
Sub-stepping the strain increments
Because of the difficulties encountered during the Newton-Raphson process during rapidly hardening/softening moduli, it is possible to subdivide the applied strain increment, , into smaller sub-steps, and do multiple return-map processes. The final returned configuration will then be dependent on the number of sub-steps. While this is simply illustrating the non-uniqueness of plasticity problems, in my experience it does adversely affect MOOSE's nonlinear convergence as some Residual calculations will take more sub-steps than other Residual calculations: in effect this is reducing the accuracy of the Jacobian.
The consistent tangent operator
MOOSE's Jacobian depends on the derivative
The quantity is called the consistent tangent operator. For pure elasticity it is simply the elastic tensor, , but it is more complicated for plasticity. Note that a small simply changes , so is capturing the change of the returned stress () with respect to a change in the trial stress (). In language, we need to the sx derivatives
The algebra is extremely tedious, but it is fairly easy for the computer. The MOOSE code contains two implementations of the consistent tangent operator. One is valid for any general model, while the other is specialized to the weak-plane case.
General consistent tangent operator
The return-map algorithm provides
Since , the consistent tangent operator is
However, note that
because the return-map algorithm guarantees that the expression inside parentheses is zero. Therefore
A similar expression holds for three other cases. There are still terms that involve derivatives of and , but these may be separated off as seen below.
The consistent tangent operator may therefore be written as
All terms but the final line have already been computed during the return-map process. The final line may be brought to the right-hand side (since ) and the resulting expression multiplied inverse 's coefficient to finally yield . This inversion, and all the multiplication of rank-four tensors may be computationally expensive, so a cheaper (but more lengthy looking) version is derived below for the capped weak-plane case.
Specialization to the weak-plane case
The return-map equations Eq. (8) are obtaining given the trial variables. Finding is really just re-solving these equations for a slightly changed trial variable. Denote
Then
In these equations
which comes from Eq. (8). The derivatives with respect to are similar but more lengthy due to both and being dependent on . The system to solve is
The Jacobian matrix is identical to the one used in the Newton-Raphson process, but of course that process has completed before calculation of the consistent tangent operator. A similar system of equations gives the derivatives with respect to .
Once the six derivatives have been computed they need to be assembled into . For instance,
so that
and other more complicated expressions appear for other components, such as
The consistent tangent operator and sub-stepping strain increments
One extra complication arises from the potential sub-stepping of the applied strain increment . At each sub-step, the six derivatives must be computed. While this may seem expensive, in my experience it increases the accuracy of the Jacobian, and the main computational expense is building and solving the system which is pretty quick for the computer to compared with the entire Newton-Raphson process.
Let the substep be
with
where is the total number of substeps. Denoting the initial stress by , and the returned stress at step by , of course the trial stress at step is
This means that
Similar inductive equations hold for the other derivatives, and note that . The derivative of is slightly different: it is
and similarly for the derivative with respect to .
Input Parameters
- admissible_stressA single admissible value of the value of the stress parameters for internal parameters = 0. This is used to initialize the return-mapping algorithm during the first nonlinear iteration. If not given then it is assumed that stress parameters = 0 is admissible.
C++ Type:std::vector<double>
Unit:(no unit assumed)
Controllable:No
Description:A single admissible value of the value of the stress parameters for internal parameters = 0. This is used to initialize the return-mapping algorithm during the first nonlinear iteration. If not given then it is assumed that stress parameters = 0 is admissible.
- base_nameOptional parameter that defines a prefix for all material properties related to this stress update model. This allows for multiple models of the same type to be used without naming conflicts.
C++ Type:std::string
Unit:(no unit assumed)
Controllable:No
Description:Optional parameter that defines a prefix for all material properties related to this stress update model. This allows for multiple models of the same type to be used without naming conflicts.
- blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Unit:(no unit assumed)
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
- boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Unit:(no unit assumed)
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
- constant_onNONEWhen ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
Default:NONE
C++ Type:MooseEnum
Unit:(no unit assumed)
Options:NONE, ELEMENT, SUBDOMAIN
Controllable:No
Description:When ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
- declare_suffixAn optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
- max_NR_iterations20Maximum number of Newton-Raphson iterations allowed during the return-map algorithm
Default:20
C++ Type:unsigned int
Unit:(no unit assumed)
Controllable:No
Description:Maximum number of Newton-Raphson iterations allowed during the return-map algorithm
- min_step_size1In order to help the Newton-Raphson procedure, the applied strain increment may be applied in sub-increments of size greater than this value. Usually it is better for Moose's nonlinear convergence to increase max_NR_iterations rather than decrease this parameter.
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:In order to help the Newton-Raphson procedure, the applied strain increment may be applied in sub-increments of size greater than this value. Usually it is better for Moose's nonlinear convergence to increase max_NR_iterations rather than decrease this parameter.
- perfect_guessTrueProvide a guess to the Newton-Raphson procedure that is the result from perfect plasticity. With severe hardening/softening this may be suboptimal.
Default:True
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Provide a guess to the Newton-Raphson procedure that is the result from perfect plasticity. With severe hardening/softening this may be suboptimal.
- perform_finite_strain_rotationsFalseTensors are correctly rotated in finite-strain simulations. For optimal performance you can set this to 'false' if you are only ever using small strains
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Tensors are correctly rotated in finite-strain simulations. For optimal performance you can set this to 'false' if you are only ever using small strains
- prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
- use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Unit:(no unit assumed)
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.
- warn_about_precision_lossFalseOutput a message to the console every time precision-loss is encountered during the Newton-Raphson process
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Output a message to the console every time precision-loss is encountered during the Newton-Raphson process
Optional Parameters
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Unit:(no unit assumed)
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Unit:(no unit assumed)
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
Unit:(no unit assumed)
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
- seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Unit:(no unit assumed)
Controllable:No
Description:The seed for the master random number generator
- smoother_function_typecosType of smoother function to use. 'cos' means (-a/pi)cos(pi x/2/a), 'polyN' means a polynomial of degree 2N+2
Default:cos
C++ Type:MooseEnum
Unit:(no unit assumed)
Options:cos, poly1, poly2, poly3
Controllable:No
Description:Type of smoother function to use. 'cos' means (-a/pi)cos(pi x/2/a), 'polyN' means a polynomial of degree 2N+2
- use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Advanced Parameters
- output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector<std::string>
Unit:(no unit assumed)
Controllable:No
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
- outputsnone Vector of output names where you would like to restrict the output of variables(s) associated with this object
Default:none
C++ Type:std::vector<OutputName>
Unit:(no unit assumed)
Controllable:No
Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object
Outputs Parameters
Input Files
- (modules/solid_mechanics/test/tests/capped_weak_plane/small_deform4.i)
- (modules/solid_mechanics/test/tests/jacobian/cto27.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/small_deform2.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/beam.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/push_and_shear.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/pull_and_shear_1step.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/small_deform9.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/small_deform6.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/throw_test.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/except5.i)
- (modules/porous_flow/test/tests/plastic_heating/compressive01.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/pull_push_h.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/small_deform11.i)
- (modules/porous_flow/test/tests/plastic_heating/shear01.i)
- (modules/solid_mechanics/test/tests/jacobian/cwp03.i)
- (modules/solid_mechanics/test/tests/multiple_two_parameter_plasticity/cycled_dp_then_wp.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/pull_push.i)
- (modules/porous_flow/test/tests/jacobian/phe01.i)
- (modules/porous_flow/test/tests/plastic_heating/tensile01.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/small_deform3.i)
- (modules/solid_mechanics/test/tests/jacobian/cwp08.i)
- (modules/solid_mechanics/test/tests/jacobian/cwp09.i)
- (modules/solid_mechanics/test/tests/multiple_two_parameter_plasticity/dp_then_wp.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/except4.i)
- (modules/solid_mechanics/test/tests/jacobian/phe01.i)
- (modules/solid_mechanics/test/tests/jacobian/cwp11.i)
- (modules/solid_mechanics/test/tests/jacobian/cwp01.i)
- (modules/solid_mechanics/test/tests/jacobian/cwp04.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/small_deform5.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/except3.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/small_deform8.i)
- (modules/solid_mechanics/test/tests/jacobian/cwp05.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/small_deform7.i)
- (modules/solid_mechanics/test/tests/jacobian/cwp06.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/except2.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/except1.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/pull_and_shear.i)
- (modules/solid_mechanics/test/tests/capped_weak_plane/small_deform10.i)
- (modules/solid_mechanics/test/tests/multiple_two_parameter_plasticity/dp_and_wp.i)
- (modules/solid_mechanics/test/tests/jacobian/cwp02.i)
- (modules/solid_mechanics/test/tests/jacobian/cwp07.i)
- (modules/solid_mechanics/test/tests/jacobian/cwp10.i)
Child Objects
(modules/solid_mechanics/test/tests/capped_weak_plane/small_deform4.i)
# Plastic deformation, compression failure
# With Young = 10, poisson=0.25 (Lame lambda=4, mu=4)
# applying the following
# deformation to the zmax surface of a unit cube:
# disp_x = 4*t
# disp_y = 3*t
# disp_z = -t
# should yield trial stress:
# stress_zz = 12*t
# stress_zx = 16*t
# stress_zy = -12*t
# Use compressive strength = 6, we should return to stress_zz = -6,
# and stress_xx = stress_yy = -2*t up to t=1 when the system is completely
# plastic, so these stress components will not change
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz'
[../]
[]
[BCs]
[./bottomx]
type = DirichletBC
variable = disp_x
boundary = back
value = 0.0
[../]
[./bottomy]
type = DirichletBC
variable = disp_y
boundary = back
value = 0.0
[../]
[./bottomz]
type = DirichletBC
variable = disp_z
boundary = back
value = 0.0
[../]
[./topx]
type = FunctionDirichletBC
variable = disp_x
boundary = front
function = 4*t
[../]
[./topy]
type = FunctionDirichletBC
variable = disp_y
boundary = front
function = 3*t
[../]
[./topz]
type = FunctionDirichletBC
variable = disp_z
boundary = front
function = -t
[../]
[]
[AuxVariables]
[./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]
[./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
[../]
[]
[Postprocessors]
[./stress_xx]
type = PointValue
point = '0 0 0'
variable = stress_xx
[../]
[./stress_xy]
type = PointValue
point = '0 0 0'
variable = stress_xy
[../]
[./stress_xz]
type = PointValue
point = '0 0 0'
variable = stress_xz
[../]
[./stress_yy]
type = PointValue
point = '0 0 0'
variable = stress_yy
[../]
[./stress_yz]
type = PointValue
point = '0 0 0'
variable = stress_yz
[../]
[./stress_zz]
type = PointValue
point = '0 0 0'
variable = stress_zz
[../]
[./strainp_xx]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xx
[../]
[./strainp_xy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xy
[../]
[./strainp_xz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xz
[../]
[./strainp_yy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yy
[../]
[./strainp_yz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yz
[../]
[./strainp_zz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_zz
[../]
[./straint_xx]
type = PointValue
point = '0 0 0'
variable = strain_xx
[../]
[./straint_xy]
type = PointValue
point = '0 0 0'
variable = strain_xy
[../]
[./straint_xz]
type = PointValue
point = '0 0 0'
variable = strain_xz
[../]
[./straint_yy]
type = PointValue
point = '0 0 0'
variable = strain_yy
[../]
[./straint_yz]
type = PointValue
point = '0 0 0'
variable = strain_yz
[../]
[./straint_zz]
type = PointValue
point = '0 0 0'
variable = strain_zz
[../]
[./f_shear]
type = PointValue
point = '0 0 0'
variable = f_shear
[../]
[./f_tensile]
type = PointValue
point = '0 0 0'
variable = f_tensile
[../]
[./f_compressive]
type = PointValue
point = '0 0 0'
variable = f_compressive
[../]
[./intnl_shear]
type = PointValue
point = '0 0 0'
variable = intnl_shear
[../]
[./intnl_tensile]
type = PointValue
point = '0 0 0'
variable = intnl_tensile
[../]
[./iter]
type = PointValue
point = '0 0 0'
variable = iter
[../]
[./ls]
type = PointValue
point = '0 0 0'
variable = ls
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningConstant
value = 80
[../]
[./tanphi]
type = SolidMechanicsHardeningConstant
value = 0.5
[../]
[./tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.1111077
[../]
[./t_strength]
type = SolidMechanicsHardeningConstant
value = 6
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 6
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '4 4'
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[../]
[./stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
tip_smoother = 0
smoothing_tol = 0
yield_function_tol = 1E-5
[../]
[]
[Executioner]
end_time = 2
dt = 1
type = Transient
[]
[Outputs]
file_base = small_deform4
csv = true
[]
(modules/solid_mechanics/test/tests/jacobian/cto27.i)
# CappedDruckerPrager and CappedWeakPlane, both with all parameters softening/hardening.
# With large tolerance in ComputeMultipleInelasticStress so that only 1 iteration is performed
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
block = 0
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Kernels]
[SolidMechanics]
displacements = 'disp_x disp_y disp_z'
[../]
[]
[UserObjects]
[./ts]
type = SolidMechanicsHardeningCubic
value_0 = 1
value_residual = 2
internal_limit = 100
[../]
[./cs]
type = SolidMechanicsHardeningCubic
value_0 = 5
value_residual = 3
internal_limit = 100
[../]
[./mc_coh]
type = SolidMechanicsHardeningCubic
value_0 = 10
value_residual = 1
internal_limit = 100
[../]
[./phi]
type = SolidMechanicsHardeningCubic
value_0 = 0.8
value_residual = 0.4
internal_limit = 50
[../]
[./psi]
type = SolidMechanicsHardeningCubic
value_0 = 0.4
value_residual = 0
internal_limit = 10
[../]
[./dp]
type = SolidMechanicsPlasticDruckerPragerHyperbolic
mc_cohesion = mc_coh
mc_friction_angle = phi
mc_dilation_angle = psi
yield_function_tolerance = 1E-11 # irrelevant here
internal_constraint_tolerance = 1E-9 # irrelevant here
[../]
[./wp_ts]
type = SolidMechanicsHardeningExponential
value_0 = 100
value_residual = 100
rate = 1
[../]
[./wp_cs]
type = SolidMechanicsHardeningCubic
value_0 = 1
value_residual = 0
internal_0 = -2
internal_limit = 0
[../]
[./wp_coh]
type = SolidMechanicsHardeningExponential
value_0 = 1
value_residual = 2
rate = 1
[../]
[./wp_tanphi]
type = SolidMechanicsHardeningExponential
value_0 = 1.0
value_residual = 0.5
rate = 2
[../]
[./wp_tanpsi]
type = SolidMechanicsHardeningExponential
value_0 = 0.1
value_residual = 0.05
rate = 3
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
block = 0
lambda = 0.1
shear_modulus = 1.0
[../]
[./strain]
type = ComputeIncrementalStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '6 5 4 5 7 2 4 2 2'
eigenstrain_name = ini_stress
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = 'dp wp'
relative_tolerance = 1E4
absolute_tolerance = 2
tangent_operator = nonlinear
[../]
[./dp]
type = CappedDruckerPragerStressUpdate
base_name = cdp
DP_model = dp
tensile_strength = ts
compressive_strength = cs
yield_function_tol = 1E-11
tip_smoother = 1
smoothing_tol = 1
[../]
[./wp]
type = CappedWeakPlaneStressUpdate
base_name = cwp
cohesion = wp_coh
tan_friction_angle = wp_tanphi
tan_dilation_angle = wp_tanpsi
tensile_strength = wp_ts
compressive_strength = wp_cs
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-11
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/small_deform2.i)
# Plastic deformation, tensile failure
# With Lame lambda=0 and Lame mu=1, applying the following
# deformation to the zmax surface of a unit cube:
# disp_z = t
# should yield trial stress:
# stress_zz = 2*t
# Use tensile strength = 1, we should return to stress_zz = 1
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz'
[../]
[]
[BCs]
[./bottomx]
type = DirichletBC
variable = disp_x
boundary = back
value = 0.0
[../]
[./bottomy]
type = DirichletBC
variable = disp_y
boundary = back
value = 0.0
[../]
[./bottomz]
type = DirichletBC
variable = disp_z
boundary = back
value = 0.0
[../]
[./topx]
type = FunctionDirichletBC
variable = disp_x
boundary = front
function = 0
[../]
[./topy]
type = FunctionDirichletBC
variable = disp_y
boundary = front
function = 0
[../]
[./topz]
type = FunctionDirichletBC
variable = disp_z
boundary = front
function = t
[../]
[]
[AuxVariables]
[./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]
[./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
[../]
[]
[Postprocessors]
[./stress_xx]
type = PointValue
point = '0 0 0'
variable = stress_xx
[../]
[./stress_xy]
type = PointValue
point = '0 0 0'
variable = stress_xy
[../]
[./stress_xz]
type = PointValue
point = '0 0 0'
variable = stress_xz
[../]
[./stress_yy]
type = PointValue
point = '0 0 0'
variable = stress_yy
[../]
[./stress_yz]
type = PointValue
point = '0 0 0'
variable = stress_yz
[../]
[./stress_zz]
type = PointValue
point = '0 0 0'
variable = stress_zz
[../]
[./strainp_xx]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xx
[../]
[./strainp_xy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xy
[../]
[./strainp_xz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xz
[../]
[./strainp_yy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yy
[../]
[./strainp_yz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yz
[../]
[./strainp_zz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_zz
[../]
[./straint_xx]
type = PointValue
point = '0 0 0'
variable = strain_xx
[../]
[./straint_xy]
type = PointValue
point = '0 0 0'
variable = strain_xy
[../]
[./straint_xz]
type = PointValue
point = '0 0 0'
variable = strain_xz
[../]
[./straint_yy]
type = PointValue
point = '0 0 0'
variable = strain_yy
[../]
[./straint_yz]
type = PointValue
point = '0 0 0'
variable = strain_yz
[../]
[./straint_zz]
type = PointValue
point = '0 0 0'
variable = strain_zz
[../]
[./f_shear]
type = PointValue
point = '0 0 0'
variable = f_shear
[../]
[./f_tensile]
type = PointValue
point = '0 0 0'
variable = f_tensile
[../]
[./f_compressive]
type = PointValue
point = '0 0 0'
variable = f_compressive
[../]
[./intnl_shear]
type = PointValue
point = '0 0 0'
variable = intnl_shear
[../]
[./intnl_tensile]
type = PointValue
point = '0 0 0'
variable = intnl_tensile
[../]
[./iter]
type = PointValue
point = '0 0 0'
variable = iter
[../]
[./ls]
type = PointValue
point = '0 0 0'
variable = ls
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningConstant
value = 30
[../]
[./tanphi]
type = SolidMechanicsHardeningConstant
value = 0.5
[../]
[./tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.1111077
[../]
[./t_strength]
type = SolidMechanicsHardeningConstant
value = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 40
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0 1'
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[../]
[./stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
tip_smoother = 0
smoothing_tol = 0
yield_function_tol = 1E-5
[../]
[]
[Executioner]
end_time = 2
dt = 1
type = Transient
[]
[Outputs]
file_base = small_deform2
csv = true
[]
(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 = ComputeIncrementalStrain
[../]
[./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/solid_mechanics/test/tests/capped_weak_plane/push_and_shear.i)
# Dynamic problem with plasticity.
# A column of material (not subject to gravity) has the z-displacement
# of its sides fixed, but the centre of its bottom side is pushed
# upwards. This causes failure in the bottom elements.
#
# The problem utilises damping in the following way.
# The DynamicStressDivergenceTensors forms the residual
# integral grad(stress) + zeta*grad(stress-dot)
# = V/L * elasticity * (du/dx + zeta * dv/dx)
# where V is the elemental volume, and L is the length-scale,
# and u is the displacement, and v is the velocity.
# The InertialForce forms the residual
# integral density * (accel + eta * velocity)
# = V * density * (a + eta * v)
# where a is the acceleration.
# So, a damped oscillator description with both these
# kernels looks like
# 0 = V * (density * a + density * eta * v + elasticity * zeta * v / L^2 + elasticity / L^2 * u)
# Critical damping is when the coefficient of v is
# 2 * sqrt(density * elasticity / L^2)
# In the case at hand, density=1E4, elasticity~1E10 (Young is 16GPa),
# L~1 to 10 (in the horizontal or vertical direction), so this coefficient ~ 1E7 to 1E6.
# Choosing eta = 1E3 and zeta = 1E-2 gives approximate critical damping.
# If zeta is high then steady-state is achieved very quickly.
#
# In the case of plasticity, the effective stiffness of the elements
# is significantly less. Therefore, the above parameters give
# overdamping.
#
# This simulation is a nice example of the irreversable and non-uniqueness
# of simulations involving plasticity. The result depends on the damping
# parameters and the time stepping.
[Mesh]
[generated_mesh]
type = GeneratedMeshGenerator
dim = 3
nx = 10
ny = 1
nz = 5
bias_z = 1.5
xmin = -10
xmax = 10
ymin = -10
ymax = 10
zmin = -100
zmax = 0
[]
[bottomz_middle]
type = BoundingBoxNodeSetGenerator
new_boundary = bottomz_middle
bottom_left = '-1 -1500 -105'
top_right = '1 1500 -95'
input = generated_mesh
[]
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
beta = 0.25 # Newmark time integration
gamma = 0.5 # Newmark time integration
eta = 1E3 #0.3E4 # higher values mean more damping via density
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[]
[Kernels]
[DynamicSolidMechanics] # zeta*K*vel + K * disp
displacements = 'disp_x disp_y disp_z'
stiffness_damping_coefficient = 1E-2 # higher values mean more damping via stiffness
hht_alpha = 0 # better nonlinear convergence than for alpha>0
[]
[inertia_x] # M*accel + eta*M*vel
type = InertialForce
use_displaced_mesh = false
variable = disp_x
velocity = vel_x
acceleration = accel_x
[]
[inertia_y]
type = InertialForce
use_displaced_mesh = false
variable = disp_y
velocity = vel_y
acceleration = accel_y
[]
[inertia_z]
type = InertialForce
use_displaced_mesh = false
variable = disp_z
velocity = vel_z
acceleration = accel_z
[]
[]
[BCs]
[no_x2]
type = DirichletBC
variable = disp_x
boundary = right
value = 0.0
[]
[no_x1]
type = DirichletBC
variable = disp_x
boundary = left
value = 0.0
[]
[no_y1]
type = DirichletBC
variable = disp_y
boundary = bottom
value = 0.0
[]
[no_y2]
type = DirichletBC
variable = disp_y
boundary = top
value = 0.0
[]
[z_fixed_sides_xmin]
type = DirichletBC
variable = disp_z
boundary = left
value = 0
[]
[z_fixed_sides_xmax]
type = DirichletBC
variable = disp_z
boundary = right
value = 0
[]
[bottomz]
type = FunctionDirichletBC
variable = disp_z
boundary = bottomz_middle
function = min(10*t,1)
[]
[]
[AuxVariables]
[accel_x]
[]
[vel_x]
[]
[accel_y]
[]
[vel_y]
[]
[accel_z]
[]
[vel_z]
[]
[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]
[accel_x] # Calculates and stores acceleration at the end of time step
type = NewmarkAccelAux
variable = accel_x
displacement = disp_x
velocity = vel_x
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
execute_on = timestep_end
[]
[accel_y]
type = NewmarkAccelAux
variable = accel_y
displacement = disp_y
velocity = vel_y
execute_on = timestep_end
[]
[vel_y]
type = NewmarkVelAux
variable = vel_y
acceleration = accel_y
execute_on = timestep_end
[]
[accel_z]
type = NewmarkAccelAux
variable = accel_z
displacement = disp_z
velocity = vel_z
execute_on = timestep_end
[]
[vel_z]
type = NewmarkVelAux
variable = vel_z
acceleration = accel_z
execute_on = timestep_end
[]
[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]
type = SolidMechanicsHardeningConstant
value = 1E6
[]
[tanphi]
type = SolidMechanicsHardeningConstant
value = 0.5
[]
[tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.166666666667
[]
[t_strength]
type = SolidMechanicsHardeningConstant
value = 1E80
[]
[c_strength]
type = SolidMechanicsHardeningConstant
value = 0
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '6.4E9 6.4E9' # young 16MPa, Poisson 0.25
[]
[strain]
type = ComputeIncrementalStrain
[]
[admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[]
[stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
tip_smoother = 0.5E6
smoothing_tol = 0.5E6
yield_function_tol = 1E-2
[]
[density]
type = GenericConstantMaterial
block = 0
prop_names = density
prop_values = 1E4
[]
[]
[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 = 1E1
nl_rel_tol = 1e-5
l_tol = 1E-10
l_max_its = 100
nl_max_its = 100
end_time = 0.5
dt = 0.1
type = Transient
[]
[Outputs]
file_base = push_and_shear
exodus = true
csv = true
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/pull_and_shear_1step.i)
# Part of the bottom (minimum z) is pulled down by a Preset displacement
# This causes tensile failure in the elements immediately above.
# Because only the bottom row of elements ever fail, and because these
# fail in the first nonlinear step, Moose correctly converges in
# 1 nonlinear step, despite this problem being inelastic.
# (If the problem had lower cohesion, then the top row would also
# fail, but in the second nonlinear step, and so the simulation
# would require at least two nonlinear steps.)
[Mesh]
[generated_mesh]
type = GeneratedMeshGenerator
dim = 3
nx = 2
ny = 1
nz = 2
xmin = -10
xmax = 10
ymin = -10
ymax = 10
zmin = -100
zmax = 0
[]
[bottomz_middle]
type = BoundingBoxNodeSetGenerator
new_boundary = bottomz_middle
bottom_left = '-1 -15 -105'
top_right = '1 15 -95'
input = generated_mesh
[]
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Kernels]
[SolidMechanics]
[../]
[]
[BCs]
[./no_x2]
type = DirichletBC
variable = disp_x
boundary = right
value = 0.0
[../]
[./no_x1]
type = DirichletBC
variable = disp_x
boundary = left
value = 0.0
[../]
[./no_y1]
type = DirichletBC
variable = disp_y
boundary = bottom
value = 0.0
[../]
[./no_y2]
type = DirichletBC
variable = disp_y
boundary = top
value = 0.0
[../]
[./z_fixed_sides_xmin]
type = DirichletBC
variable = disp_z
boundary = left
value = 0
[../]
[./z_fixed_sides_xmax]
type = DirichletBC
variable = disp_z
boundary = right
value = 0
[../]
[./bottomz]
type = FunctionDirichletBC
variable = disp_z
boundary = bottomz_middle
function = -1
[../]
[]
[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 = 1E60
value_residual = 1E60
internal_limit = 0.01E8
[../]
[./tanphi]
type = SolidMechanicsHardeningCubic
value_0 = 0.5
value_residual = 0.2
internal_limit = 0.01E8
[../]
[./tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.166666666667
[../]
[./t_strength]
type = SolidMechanicsHardeningConstant
value = 0
[../]
[./c_strength]
type = SolidMechanicsHardeningCubic
value_0 = 1E80
value_residual = 1E80
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 = ComputeIncrementalStrain
[../]
[./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 = 1
tip_smoother = 0
smoothing_tol = 0
yield_function_tol = 1E-2
perfect_guess = true
min_step_size = 1
[../]
[]
[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 = 1E1
nl_rel_tol = 1e-5
l_tol = 1E-10
l_max_its = 100
nl_max_its = 100
end_time = 1.0
dt = 1.0
type = Transient
[]
[Outputs]
file_base = pull_and_shear_1step
exodus = true
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/small_deform9.i)
# apply a shear deformation to observe shear hardening.
# Shear gives q_trial = 2*Sqrt(20), p_trial=0
# The solution given by MOOSE correctly satisfies the equations
# 0 = f = q + p*tan(phi) - coh
# 0 = p - p_trial + ga * Ezzzz * dg/dp
# 0 = q - q_trial + ga * Ezxzx * dg/dq
# Here dg/dp = tan(psi), and dg/dq = 1
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz'
[../]
[]
[BCs]
[./bottomx]
type = DirichletBC
variable = disp_x
boundary = back
value = 0.0
[../]
[./bottomy]
type = DirichletBC
variable = disp_y
boundary = back
value = 0.0
[../]
[./bottomz]
type = DirichletBC
variable = disp_z
boundary = back
value = 0.0
[../]
[./topx]
type = FunctionDirichletBC
variable = disp_x
boundary = front
function = 't'
[../]
[./topy]
type = FunctionDirichletBC
variable = disp_y
boundary = front
function = '2*t'
[../]
[./topz]
type = FunctionDirichletBC
variable = disp_z
boundary = front
function = '0'
[../]
[]
[AuxVariables]
[./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]
[./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
[../]
[]
[Postprocessors]
[./stress_xx]
type = PointValue
point = '0 0 0'
variable = stress_xx
[../]
[./stress_xy]
type = PointValue
point = '0 0 0'
variable = stress_xy
[../]
[./stress_xz]
type = PointValue
point = '0 0 0'
variable = stress_xz
[../]
[./stress_yy]
type = PointValue
point = '0 0 0'
variable = stress_yy
[../]
[./stress_yz]
type = PointValue
point = '0 0 0'
variable = stress_yz
[../]
[./stress_zz]
type = PointValue
point = '0 0 0'
variable = stress_zz
[../]
[./strainp_xx]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xx
[../]
[./strainp_xy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xy
[../]
[./strainp_xz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xz
[../]
[./strainp_yy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yy
[../]
[./strainp_yz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yz
[../]
[./strainp_zz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_zz
[../]
[./straint_xx]
type = PointValue
point = '0 0 0'
variable = strain_xx
[../]
[./straint_xy]
type = PointValue
point = '0 0 0'
variable = strain_xy
[../]
[./straint_xz]
type = PointValue
point = '0 0 0'
variable = strain_xz
[../]
[./straint_yy]
type = PointValue
point = '0 0 0'
variable = strain_yy
[../]
[./straint_yz]
type = PointValue
point = '0 0 0'
variable = strain_yz
[../]
[./straint_zz]
type = PointValue
point = '0 0 0'
variable = strain_zz
[../]
[./f_shear]
type = PointValue
point = '0 0 0'
variable = f_shear
[../]
[./f_tensile]
type = PointValue
point = '0 0 0'
variable = f_tensile
[../]
[./f_compressive]
type = PointValue
point = '0 0 0'
variable = f_compressive
[../]
[./intnl_shear]
type = PointValue
point = '0 0 0'
variable = intnl_shear
[../]
[./intnl_tensile]
type = PointValue
point = '0 0 0'
variable = intnl_tensile
[../]
[./iter]
type = PointValue
point = '0 0 0'
variable = iter
[../]
[./ls]
type = PointValue
point = '0 0 0'
variable = ls
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningExponential
value_0 = 1
value_residual = 2
rate = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningExponential
value_0 = 1.0
value_residual = 0.5
rate = 2
[../]
[./tanpsi]
type = SolidMechanicsHardeningExponential
value_0 = 0.1
value_residual = 0.05
rate = 1
[../]
[./t_strength]
type = SolidMechanicsHardeningConstant
value = 1E8
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 1E8
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '4 4'
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[../]
[./stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-3
[../]
[]
[Executioner]
end_time = 1
dt = 1
type = Transient
[]
[Outputs]
file_base = small_deform9
[./csv]
type = CSV
[../]
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/small_deform6.i)
# Plastic deformation, both tensile and shear failure
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz'
[../]
[]
[BCs]
[./bottomx]
type = DirichletBC
variable = disp_x
boundary = back
value = 0.0
[../]
[./bottomy]
type = DirichletBC
variable = disp_y
boundary = back
value = 0.0
[../]
[./bottomz]
type = DirichletBC
variable = disp_z
boundary = back
value = 0.0
[../]
[./topx]
type = FunctionDirichletBC
variable = disp_x
boundary = front
function = 'if(t<30,0.2*t,6)'
[../]
[./topy]
type = FunctionDirichletBC
variable = disp_y
boundary = front
function = 'if(t<30,if(t<10,0,t),30-0.2*t)'
[../]
[./topz]
type = FunctionDirichletBC
variable = disp_z
boundary = front
function = 'if(t<15,3*t,45)+if(t<30,0,45-3*t)'
[../]
[]
[AuxVariables]
[./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]
[./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
[../]
[]
[Postprocessors]
[./stress_xx]
type = PointValue
point = '0 0 0'
variable = stress_xx
[../]
[./stress_xy]
type = PointValue
point = '0 0 0'
variable = stress_xy
[../]
[./stress_xz]
type = PointValue
point = '0 0 0'
variable = stress_xz
[../]
[./stress_yy]
type = PointValue
point = '0 0 0'
variable = stress_yy
[../]
[./stress_yz]
type = PointValue
point = '0 0 0'
variable = stress_yz
[../]
[./stress_zz]
type = PointValue
point = '0 0 0'
variable = stress_zz
[../]
[./strainp_xx]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xx
[../]
[./strainp_xy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xy
[../]
[./strainp_xz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xz
[../]
[./strainp_yy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yy
[../]
[./strainp_yz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yz
[../]
[./strainp_zz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_zz
[../]
[./straint_xx]
type = PointValue
point = '0 0 0'
variable = strain_xx
[../]
[./straint_xy]
type = PointValue
point = '0 0 0'
variable = strain_xy
[../]
[./straint_xz]
type = PointValue
point = '0 0 0'
variable = strain_xz
[../]
[./straint_yy]
type = PointValue
point = '0 0 0'
variable = strain_yy
[../]
[./straint_yz]
type = PointValue
point = '0 0 0'
variable = strain_yz
[../]
[./straint_zz]
type = PointValue
point = '0 0 0'
variable = strain_zz
[../]
[./f_shear]
type = PointValue
point = '0 0 0'
variable = f_shear
[../]
[./f_tensile]
type = PointValue
point = '0 0 0'
variable = f_tensile
[../]
[./f_compressive]
type = PointValue
point = '0 0 0'
variable = f_compressive
[../]
[./intnl_shear]
type = PointValue
point = '0 0 0'
variable = intnl_shear
[../]
[./intnl_tensile]
type = PointValue
point = '0 0 0'
variable = intnl_tensile
[../]
[./iter]
type = PointValue
point = '0 0 0'
variable = iter
[../]
[./ls]
type = PointValue
point = '0 0 0'
variable = ls
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningConstant
value = 20
[../]
[./tanphi]
type = SolidMechanicsHardeningConstant
value = 0.5
[../]
[./tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.1
[../]
[./t_strength]
type = SolidMechanicsHardeningConstant
value = 20
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 0
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '4 4'
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[../]
[./stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 5
smoothing_tol = 5
yield_function_tol = 1E-10
perfect_guess = false
[../]
[]
[Executioner]
end_time = 40
dt = 1
type = Transient
[]
[Outputs]
file_base = small_deform6
csv = true
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/throw_test.i)
# Illustrates throwing an Exception from a Material. In this case we
# don't actually recover from the segfault (so it is a RunException
# test) but in practice one could do so. The purpose of this test is
# to ensure that exceptions can be thrown from Materials with stateful
# material properties without reading/writing to/from uninitialized
# memory.
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
[../]
[]
[BCs]
[./bottomx]
type = DirichletBC
variable = disp_x
boundary = back
value = 0.0
[../]
[./bottomy]
type = DirichletBC
variable = disp_y
boundary = back
value = 0.0
[../]
[./bottomz]
type = DirichletBC
variable = disp_z
boundary = back
value = 0.0
[../]
[./topx]
type = FunctionDirichletBC
variable = disp_x
boundary = front
function = 0
[../]
[./topy]
type = FunctionDirichletBC
variable = disp_y
boundary = front
function = 0
[../]
[./topz]
type = FunctionDirichletBC
variable = disp_z
boundary = front
function = t
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningConstant
value = 20
[../]
[./tanphi]
type = SolidMechanicsHardeningConstant
value = 0.5
[../]
[./tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.1
[../]
[./t_strength]
type = SolidMechanicsHardeningCubic
value_0 = 1
value_residual = 2
internal_limit = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 100
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0 1'
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[../]
[./stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 1
tip_smoother = 5
smoothing_tol = 5
yield_function_tol = 1E-10
[../]
[]
[Executioner]
end_time = 1
dt = 1
dtmin = 1
type = Transient
[]
[Outputs]
file_base = SEGFAULT
csv = true
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/except5.i)
# Exception: incorrect userobject types
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningConstant
value = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningConstant
value = 0.5
[../]
[./tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.1111077
[../]
[./t_strength]
type = SolidMechanicsHardeningConstant
value = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = -2
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0 1E6'
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[../]
[./stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-5
[../]
[]
[Executioner]
end_time = 1
dt = 1
type = Transient
[]
(modules/porous_flow/test/tests/plastic_heating/compressive01.i)
# Tensile heating, using capped weak-plane plasticity
# z_disp(z=1) = -t
# totalstrain_zz = -t
# with C_ijkl = 0.5 0.25
# stress_zz = -t, but with compressive_strength = 1, stress_zz = max(-t, -1)
# so plasticstrain_zz = -(t - 1)
# heat_energy_rate = coeff * (t - 1)
# Heat capacity of rock = specific_heat_cap * density = 4
# So temperature of rock should be:
# (1 - porosity) * 4 * T = (1 - porosity) * coeff * (t - 1)
[Mesh]
type = GeneratedMesh
dim = 3
xmin = -10
xmax = 10
zmin = 0
zmax = 1
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
[]
[Variables]
[temperature]
[]
[]
[Kernels]
[energy_dot]
type = PorousFlowEnergyTimeDerivative
variable = temperature
base_name = non_existent
[]
[phe]
type = PorousFlowPlasticHeatEnergy
variable = temperature
[]
[]
[AuxVariables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[]
[AuxKernels]
[disp_z]
type = FunctionAux
variable = disp_z
function = '-z*t'
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = temperature
number_fluid_phases = 0
number_fluid_components = 0
[]
[coh]
type = TensorMechanicsHardeningConstant
value = 100
[]
[tanphi]
type = TensorMechanicsHardeningConstant
value = 1.0
[]
[t_strength]
type = TensorMechanicsHardeningConstant
value = 1
[]
[c_strength]
type = TensorMechanicsHardeningConstant
value = 1
[]
[]
[Materials]
[rock_internal_energy]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 2
density = 2
[]
[temp]
type = PorousFlowTemperature
temperature = temperature
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.2
[]
[phe]
type = ComputePlasticHeatEnergy
[]
[elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0.5 0.25'
[]
[strain]
type = ComputeIncrementalStrain
displacements = 'disp_x disp_y disp_z'
[]
[admissible]
type = ComputeMultipleInelasticStress
inelastic_models = mc
perform_finite_strain_rotations = false
[]
[mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanphi
tensile_strength = t_strength
compressive_strength = c_strength
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-10
perfect_guess = true
[]
[]
[Postprocessors]
[temp]
type = PointValue
point = '0 0 0'
variable = temperature
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 10
[]
[Outputs]
file_base = compressive01
csv = true
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/pull_push_h.i)
# A column of elements has its bottom pulled down, and then pushed up again.
# Hardening of the tensile strength means that the top element also
# experiences plastic deformation
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 2
xmin = -10
xmax = 10
ymin = -10
ymax = 10
zmin = -100
zmax = 0
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Kernels]
[SolidMechanics]
[../]
[]
[BCs]
[./no_x2]
type = DirichletBC
variable = disp_x
boundary = right
value = 0.0
[../]
[./no_x1]
type = DirichletBC
variable = disp_x
boundary = left
value = 0.0
[../]
[./no_y1]
type = DirichletBC
variable = disp_y
boundary = bottom
value = 0.0
[../]
[./no_y2]
type = DirichletBC
variable = disp_y
boundary = top
value = 0.0
[../]
[./topz]
type = DirichletBC
variable = disp_z
boundary = front
value = 0
[../]
[./bottomz]
type = FunctionDirichletBC
variable = disp_z
boundary = back
function = 'if(t>1,-2.0+t,-t)'
[../]
[]
[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 = 1E6
internal_limit = 0.01
[../]
[./tanphi]
type = SolidMechanicsHardeningCubic
value_0 = 0.5
value_residual = 0.2
internal_limit = 0.01
[../]
[./tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.166666666667
[../]
[./t_strength]
type = SolidMechanicsHardeningCubic
value_0 = 0
value_residual = 1E8
internal_limit = 0.1
[../]
[./c_strength]
type = SolidMechanicsHardeningCubic
value_0 = 1E8
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 = ComputeIncrementalStrain
[../]
[./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 = 0
smoothing_tol = 0
yield_function_tol = 1E-5
perfect_guess = false
min_step_size = 0.1
[../]
[]
[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 = 3.0
dt = 0.1
type = Transient
[]
[Outputs]
file_base = pull_push_h
exodus = true
csv = true
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/small_deform11.i)
# use an initial stress, then apply a shear deformation and tensile stretch to observe all hardening.
# Here p_trial=12, q_trial=2*Sqrt(20)
# MOOSE yields:
# q_returned = 1.696
# p_returned = 0.100
# intnl_shear = 1.81
# intnl_tens = 0.886
# These give, at the returned point
# cohesion = 1.84
# tanphi = 0.513
# tanpsi = 0.058
# tensile = 0.412
# This means that
# f_shear = -0.0895
# f_tensile = -0.312
# Note that these are within smoothing_tol (=1) of each other
# Hence, smoothing must be used:
# ismoother = 0.0895
# (which gives the yield function value = 0)
# smoother = 0.328
# This latter gives dg/dq = 0.671, dg/dp = 0.368
# for the flow directions. Finally ga = 2.70, and
# the returned point satisfies the normality conditions.
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
eigenstrain_names = ini_stress
generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz'
[../]
[]
[BCs]
[./bottomx]
type = DirichletBC
variable = disp_x
boundary = back
value = 0.0
[../]
[./bottomy]
type = DirichletBC
variable = disp_y
boundary = back
value = 0.0
[../]
[./bottomz]
type = DirichletBC
variable = disp_z
boundary = back
value = 0.0
[../]
[./topx]
type = FunctionDirichletBC
variable = disp_x
boundary = front
function = '0.5*t'
[../]
[./topy]
type = FunctionDirichletBC
variable = disp_y
boundary = front
function = 't'
[../]
[./topz]
type = FunctionDirichletBC
variable = disp_z
boundary = front
function = '0.5*t'
[../]
[]
[AuxVariables]
[./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]
[./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
[../]
[]
[Postprocessors]
[./stress_xx]
type = PointValue
point = '0 0 0'
variable = stress_xx
[../]
[./stress_xy]
type = PointValue
point = '0 0 0'
variable = stress_xy
[../]
[./stress_xz]
type = PointValue
point = '0 0 0'
variable = stress_xz
[../]
[./stress_yy]
type = PointValue
point = '0 0 0'
variable = stress_yy
[../]
[./stress_yz]
type = PointValue
point = '0 0 0'
variable = stress_yz
[../]
[./stress_zz]
type = PointValue
point = '0 0 0'
variable = stress_zz
[../]
[./strainp_xx]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xx
[../]
[./strainp_xy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xy
[../]
[./strainp_xz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xz
[../]
[./strainp_yy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yy
[../]
[./strainp_yz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yz
[../]
[./strainp_zz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_zz
[../]
[./straint_xx]
type = PointValue
point = '0 0 0'
variable = strain_xx
[../]
[./straint_xy]
type = PointValue
point = '0 0 0'
variable = strain_xy
[../]
[./straint_xz]
type = PointValue
point = '0 0 0'
variable = strain_xz
[../]
[./straint_yy]
type = PointValue
point = '0 0 0'
variable = strain_yy
[../]
[./straint_yz]
type = PointValue
point = '0 0 0'
variable = strain_yz
[../]
[./straint_zz]
type = PointValue
point = '0 0 0'
variable = strain_zz
[../]
[./f_shear]
type = PointValue
point = '0 0 0'
variable = f_shear
[../]
[./f_tensile]
type = PointValue
point = '0 0 0'
variable = f_tensile
[../]
[./f_compressive]
type = PointValue
point = '0 0 0'
variable = f_compressive
[../]
[./intnl_shear]
type = PointValue
point = '0 0 0'
variable = intnl_shear
[../]
[./intnl_tensile]
type = PointValue
point = '0 0 0'
variable = intnl_tensile
[../]
[./iter]
type = PointValue
point = '0 0 0'
variable = iter
[../]
[./ls]
type = PointValue
point = '0 0 0'
variable = ls
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningExponential
value_0 = 1
value_residual = 2
rate = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningExponential
value_0 = 1.0
value_residual = 0.5
rate = 2
[../]
[./tanpsi]
type = SolidMechanicsHardeningExponential
value_0 = 0.1
value_residual = 0.05
rate = 1
[../]
[./t_strength]
type = SolidMechanicsHardeningExponential
value_0 = 1
value_residual = 0
rate = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 1E8
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
lambda = 4.0
shear_modulus = 4.0
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '0 0 2 0 0 4 2 4 6'
eigenstrain_name = ini_stress
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[../]
[./stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-3
perfect_guess = false
[../]
[]
[Executioner]
end_time = 1
dt = 1
type = Transient
[]
[Outputs]
file_base = small_deform11
[./csv]
type = CSV
[../]
[]
(modules/porous_flow/test/tests/plastic_heating/shear01.i)
# Tensile heating, using capped weak-plane plasticity
# x_disp(z=1) = t
# totalstrain_xz = t
# with C_ijkl = 0.5 0.25
# stress_zx = stress_xz = 0.25*t, so q=0.25*t, but
# with cohesion=1 and tan(phi)=1: max(q)=1. With tan(psi)=0,
# the plastic return is always to (p, q) = (0, 1),
# so plasticstrain_zx = max(t - 4, 0)
# heat_energy_rate = coeff * (t - 4) for t>4
# Heat capacity of rock = specific_heat_cap * density = 4
# So temperature of rock should be:
# (1 - porosity) * 4 * T = (1 - porosity) * coeff * (t - 4)
[Mesh]
type = GeneratedMesh
dim = 3
xmin = -10
xmax = 10
zmin = 0
zmax = 1
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
[]
[Variables]
[temperature]
[]
[]
[Kernels]
[energy_dot]
type = PorousFlowEnergyTimeDerivative
variable = temperature
base_name = non_existent
[]
[phe]
type = PorousFlowPlasticHeatEnergy
variable = temperature
coeff = 8
[]
[]
[AuxVariables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[]
[AuxKernels]
[disp_x]
type = FunctionAux
variable = disp_x
function = 'z*t'
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = temperature
number_fluid_phases = 0
number_fluid_components = 0
[]
[coh]
type = TensorMechanicsHardeningConstant
value = 1
[]
[tanphi]
type = TensorMechanicsHardeningConstant
value = 1.0
[]
[tanpsi]
type = TensorMechanicsHardeningConstant
value = 0.0
[]
[t_strength]
type = TensorMechanicsHardeningConstant
value = 10
[]
[c_strength]
type = TensorMechanicsHardeningConstant
value = 10
[]
[]
[Materials]
[rock_internal_energy]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 2
density = 2
[]
[temp]
type = PorousFlowTemperature
temperature = temperature
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.7
[]
[phe]
type = ComputePlasticHeatEnergy
[]
[elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0.5 0.25'
[]
[strain]
type = ComputeIncrementalStrain
displacements = 'disp_x disp_y disp_z'
[]
[admissible]
type = ComputeMultipleInelasticStress
inelastic_models = mc
perform_finite_strain_rotations = false
[]
[mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-10
perfect_guess = true
[]
[]
[Postprocessors]
[temp]
type = PointValue
point = '0 0 0'
variable = temperature
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 10
[]
[Outputs]
file_base = shear01
csv = true
[]
(modules/solid_mechanics/test/tests/jacobian/cwp03.i)
# Capped weak-plane plasticity
# checking jacobian for tensile failure, with some shear
[Mesh]
type = GeneratedMesh
dim = 3
[]
[GlobalParams]
block = 0
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Kernels]
[SolidMechanics]
displacements = 'disp_x disp_y disp_z'
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningExponential
value_0 = 100
value_residual = 2
rate = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningExponential
value_0 = 1.0
value_residual = 0.5
rate = 2
[../]
[./tanpsi]
type = SolidMechanicsHardeningExponential
value_0 = 0.1
value_residual = 0.05
rate = 1
[../]
[./t_strength]
type = SolidMechanicsHardeningExponential
value_0 = 1
value_residual = 1
rate = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 100
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
lambda = 1.0
shear_modulus = 2.0
[../]
[./strain]
type = ComputeIncrementalStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '0 0 -2 0 0 1 -2 1 2'
eigenstrain_name = ini_stress
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = mc
tangent_operator = nonlinear
[../]
[./mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 1
smoothing_tol = 2
yield_function_tol = 1E-10
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
#petsc_options = '-snes_test_display'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
[]
(modules/solid_mechanics/test/tests/multiple_two_parameter_plasticity/cycled_dp_then_wp.i)
# Use ComputeMultipleInelasticStress with two inelastic models: CappedDruckerPrager and CappedWeakPlane.
# The relative_tolerance and absolute_tolerance parameters are set very large so that
# only one iteration is performed. This is the algorithm that FLAC uses to model
# jointed rocks, only Capped-Mohr-Coulomb is used instead of CappedDruckerPrager
#
# In this test "cycle_models=true" so that in the first timestep only
# CappedDruckerPrager is used, while in the second timestep only
# CappedWeakPlane is used.
#
# initial_stress = diag(1E3, 1E3, 1E3)
# The CappedDruckerPrager has tensile strength 3E2 and large cohesion,
# so the stress initially returns to diag(1E2, 1E2, 1E2)
# The CappedWeakPlane has tensile strength zero and large cohesion,
# so the stress returns to diag(1E2 - v/(1-v)*1E2, 1E2 - v/(1-v)*1E2, 0)
# where v=0.2 is the Poisson's ratio
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
eigenstrain_names = ini_stress
generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz'
[../]
[]
[BCs]
[./x]
type = FunctionDirichletBC
variable = disp_x
boundary = 'front back'
function = 0
[../]
[./y]
type = FunctionDirichletBC
variable = disp_y
boundary = 'front back'
function = 0
[../]
[./z]
type = FunctionDirichletBC
variable = disp_z
boundary = 'front back'
function = 0
[../]
[]
[AuxVariables]
[./yield_fcn_dp]
order = CONSTANT
family = MONOMIAL
[../]
[./yield_fcn_wp]
order = CONSTANT
family = MONOMIAL
[../]
[./tensile_cdp]
order = CONSTANT
family = MONOMIAL
[../]
[./tensile_cwp]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./yield_fcn_dp_auxk]
type = MaterialStdVectorAux
index = 1 # this is the tensile yield function - it should be zero
property = cdp_plastic_yield_function
variable = yield_fcn_dp
[../]
[./yield_fcn_wp_auxk]
type = MaterialStdVectorAux
index = 1 # this is the tensile yield function - it should be zero
property = cwp_plastic_yield_function
variable = yield_fcn_wp
[../]
[./tensile_cdp]
type = MaterialStdVectorAux
index = 1
property = cdp_plastic_internal_parameter
variable = tensile_cdp
[../]
[./tensile_cwp]
type = MaterialStdVectorAux
index = 1
property = cwp_plastic_internal_parameter
variable = tensile_cwp
[../]
[]
[Postprocessors]
[./s_xx]
type = PointValue
point = '0 0 0'
variable = stress_xx
[../]
[./s_xy]
type = PointValue
point = '0 0 0'
variable = stress_xy
[../]
[./s_xz]
type = PointValue
point = '0 0 0'
variable = stress_xz
[../]
[./s_yy]
type = PointValue
point = '0 0 0'
variable = stress_yy
[../]
[./s_yz]
type = PointValue
point = '0 0 0'
variable = stress_yz
[../]
[./s_zz]
type = PointValue
point = '0 0 0'
variable = stress_zz
[../]
[./i_cdp]
type = PointValue
point = '0 0 0'
variable = tensile_cdp
[../]
[./i_cwp]
type = PointValue
point = '0 0 0'
variable = tensile_cwp
[../]
[]
[UserObjects]
[./ts]
type = SolidMechanicsHardeningConstant
value = 300
[../]
[./cs]
type = SolidMechanicsHardeningConstant
value = 1E4
[../]
[./mc_coh]
type = SolidMechanicsHardeningConstant
value = 1E4
[../]
[./mc_phi]
type = SolidMechanicsHardeningConstant
value = 20
convert_to_radians = true
[../]
[./mc_psi]
type = SolidMechanicsHardeningConstant
value = 0
[../]
[./dp]
type = SolidMechanicsPlasticDruckerPrager
mc_cohesion = mc_coh
mc_friction_angle = mc_phi
mc_dilation_angle = mc_psi
internal_constraint_tolerance = 1 # irrelevant here
yield_function_tolerance = 1 # irrelevant here
[../]
[./wp_coh]
type = SolidMechanicsHardeningConstant
value = 1E4
[../]
[./wp_tanphi]
type = SolidMechanicsHardeningConstant
value = 0.5
[../]
[./wp_tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.1111077
[../]
[./wp_t_strength]
type = SolidMechanicsHardeningConstant
value = 0
[../]
[./wp_c_strength]
type = SolidMechanicsHardeningConstant
value = 1E4
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
poissons_ratio = 0.2
youngs_modulus = 1.0
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '1E3 0 0 0 1E3 0 0 0 1E3'
eigenstrain_name = ini_stress
[../]
[./admissible]
type = ComputeMultipleInelasticStress
relative_tolerance = 1E4
absolute_tolerance = 2
inelastic_models = 'cdp cwp'
perform_finite_strain_rotations = false
cycle_models = true
[../]
[./cdp]
type = CappedDruckerPragerStressUpdate
base_name = cdp
DP_model = dp
tensile_strength = ts
compressive_strength = cs
yield_function_tol = 1E-5
tip_smoother = 1E3
smoothing_tol = 1E3
[../]
[./cwp]
type = CappedWeakPlaneStressUpdate
base_name = cwp
cohesion = wp_coh
tan_friction_angle = wp_tanphi
tan_dilation_angle = wp_tanpsi
tensile_strength = wp_t_strength
compressive_strength = wp_c_strength
tip_smoother = 1E3
smoothing_tol = 1E3
yield_function_tol = 1E-5
[../]
[]
[Executioner]
end_time = 2
dt = 1
type = Transient
[]
[Outputs]
file_base = cycled_dp_then_wp
csv = true
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/pull_push.i)
# A column of elements has its bottom pulled down, and then pushed up again.
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 2
xmin = -10
xmax = 10
ymin = -10
ymax = 10
zmin = -100
zmax = 0
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Kernels]
[SolidMechanics]
[../]
[]
[BCs]
[./no_x2]
type = DirichletBC
variable = disp_x
boundary = right
value = 0.0
[../]
[./no_x1]
type = DirichletBC
variable = disp_x
boundary = left
value = 0.0
[../]
[./no_y1]
type = DirichletBC
variable = disp_y
boundary = bottom
value = 0.0
[../]
[./no_y2]
type = DirichletBC
variable = disp_y
boundary = top
value = 0.0
[../]
[./topz]
type = DirichletBC
variable = disp_z
boundary = front
value = 0
[../]
[./bottomz]
type = FunctionDirichletBC
variable = disp_z
boundary = back
function = 'if(t>1,-2.0+t,-t)'
[../]
[]
[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 = 1E6
internal_limit = 0.01
[../]
[./tanphi]
type = SolidMechanicsHardeningCubic
value_0 = 0.5
value_residual = 0.2
internal_limit = 0.01
[../]
[./tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.166666666667
[../]
[./t_strength]
type = SolidMechanicsHardeningConstant
value = 2E6
[../]
[./c_strength]
type = SolidMechanicsHardeningCubic
value_0 = 1E8
value_residual = 0.0
internal_limit = 0.01
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
lambda = 6.4e9
shear_modulus = 6.4e9 # young 16MPa, Poisson 0.25
[../]
[./strain]
type = ComputeIncrementalStrain
[../]
[./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 = 10
tip_smoother = 0
smoothing_tol = 0
yield_function_tol = 1E-2
perfect_guess = false
min_step_size = 1
[../]
[]
[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 = 1E1
nl_rel_tol = 1e-5
l_tol = 1E-10
l_max_its = 100
nl_max_its = 100
end_time = 3.0
dt = 0.1
type = Transient
[]
[Outputs]
file_base = pull_push
exodus = true
csv = true
[]
(modules/porous_flow/test/tests/jacobian/phe01.i)
# Capped weak-plane plasticity, Kernel = PorousFlowPlasticHeatEnergy
[Mesh]
type = GeneratedMesh
dim = 3
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[temperature]
[]
[]
[ICs]
[disp_x]
type = RandomIC
variable = disp_x
min = -0.1
max = 0.1
[]
[disp_y]
type = RandomIC
variable = disp_y
min = -0.1
max = 0.1
[]
[disp_z]
type = RandomIC
variable = disp_z
min = -0.1
max = 0.1
[]
[temp]
type = RandomIC
variable = temperature
min = 0.1
max = 0.2
[]
[]
[Kernels]
[phe]
type = PorousFlowPlasticHeatEnergy
variable = temperature
[]
[dummy_disp_x]
type = PorousFlowPlasticHeatEnergy
coeff = -1.3
variable = disp_x
[]
[dummy_disp_y]
type = PorousFlowPlasticHeatEnergy
coeff = 1.1
variable = disp_y
[]
[dummy_disp_z]
type = PorousFlowPlasticHeatEnergy
coeff = 0.2
variable = disp_z
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'temperature disp_x disp_y disp_z'
number_fluid_phases = 0
number_fluid_components = 0
[]
[coh]
type = TensorMechanicsHardeningExponential
value_0 = 1
value_residual = 2
rate = 1
[]
[tanphi]
type = TensorMechanicsHardeningExponential
value_0 = 1.0
value_residual = 0.5
rate = 2
[]
[tanpsi]
type = TensorMechanicsHardeningExponential
value_0 = 0.1
value_residual = 0.05
rate = 3
[]
[t_strength]
type = TensorMechanicsHardeningExponential
value_0 = 100
value_residual = 100
rate = 1
[]
[c_strength]
type = TensorMechanicsHardeningCubic
value_0 = 1
value_residual = 0
internal_0 = -2
internal_limit = 0
[]
[]
[Materials]
[temp]
type = PorousFlowTemperature
temperature = temperature
[]
[porosity]
type = PorousFlowPorosity
thermal = true
mechanical = true
porosity_zero = 0.3
thermal_expansion_coeff = 1.3
[]
[volstrain]
type = PorousFlowVolumetricStrain
[]
[phe]
type = ComputePlasticHeatEnergy
[]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
lambda = 1.0
shear_modulus = 2.0
[]
[strain]
type = ComputeIncrementalStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = ini_stress
[]
[ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '0 0 0 0 0 1 0 1 -1.5'
eigenstrain_name = ini_stress
[]
[admissible]
type = ComputeMultipleInelasticStress
inelastic_models = mc
tangent_operator = nonlinear
[]
[mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-10
perfect_guess = false
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
[]
(modules/porous_flow/test/tests/plastic_heating/tensile01.i)
# Tensile heating, using capped weak-plane plasticity
# z_disp(z=1) = t
# totalstrain_zz = t
# with C_ijkl = 0.5 0.25
# stress_zz = t, but with tensile_strength = 1, stress_zz = min(t, 1)
# so plasticstrain_zz = t - 1
# heat_energy_rate = coeff * (t - 1)
# Heat capacity of rock = specific_heat_cap * density = 4
# So temperature of rock should be:
# (1 - porosity) * 4 * T = (1 - porosity) * coeff * (t - 1)
[Mesh]
type = GeneratedMesh
dim = 3
xmin = -10
xmax = 10
zmin = 0
zmax = 1
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
[]
[Variables]
[temperature]
[]
[]
[Kernels]
[energy_dot]
type = PorousFlowEnergyTimeDerivative
variable = temperature
base_name = non_existent
[]
[phe]
type = PorousFlowPlasticHeatEnergy
variable = temperature
[]
[]
[AuxVariables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[]
[AuxKernels]
[disp_z]
type = FunctionAux
variable = disp_z
function = z*t
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = temperature
number_fluid_phases = 0
number_fluid_components = 0
[]
[coh]
type = TensorMechanicsHardeningConstant
value = 100
[]
[tanphi]
type = TensorMechanicsHardeningConstant
value = 1.0
[]
[t_strength]
type = TensorMechanicsHardeningConstant
value = 1
[]
[c_strength]
type = TensorMechanicsHardeningConstant
value = 1
[]
[]
[Materials]
[rock_internal_energy]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 2
density = 2
[]
[temp]
type = PorousFlowTemperature
temperature = temperature
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.2
[]
[phe]
type = ComputePlasticHeatEnergy
[]
[elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0.5 0.25'
[]
[strain]
type = ComputeIncrementalStrain
displacements = 'disp_x disp_y disp_z'
[]
[admissible]
type = ComputeMultipleInelasticStress
inelastic_models = mc
perform_finite_strain_rotations = false
[]
[mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanphi
tensile_strength = t_strength
compressive_strength = c_strength
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-10
perfect_guess = true
[]
[]
[Postprocessors]
[temp]
type = PointValue
point = '0 0 0'
variable = temperature
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 10
[]
[Outputs]
file_base = tensile01
csv = true
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/small_deform3.i)
# Plastic deformation, tensile failure
# With Young = 10, poisson=0.25 (Lame lambda=4, mu=4)
# applying the following
# deformation to the zmax surface of a unit cube:
# disp_x = 4*t
# disp_y = 3*t
# disp_z = t
# should yield trial stress:
# stress_zz = 12*t
# stress_zx = 16*t
# stress_zy = 12*t
# Use tensile strength = 6, we should return to stress_zz = 6,
# and stress_xx = stress_yy = 2*t up to t=1 when the system is completely
# plastic, so these stress components will not change
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz'
[../]
[]
[BCs]
[./bottomx]
type = DirichletBC
variable = disp_x
boundary = back
value = 0.0
[../]
[./bottomy]
type = DirichletBC
variable = disp_y
boundary = back
value = 0.0
[../]
[./bottomz]
type = DirichletBC
variable = disp_z
boundary = back
value = 0.0
[../]
[./topx]
type = FunctionDirichletBC
variable = disp_x
boundary = front
function = 4*t
[../]
[./topy]
type = FunctionDirichletBC
variable = disp_y
boundary = front
function = 3*t
[../]
[./topz]
type = FunctionDirichletBC
variable = disp_z
boundary = front
function = t
[../]
[]
[AuxVariables]
[./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]
[./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
[../]
[]
[Postprocessors]
[./stress_xx]
type = PointValue
point = '0 0 0'
variable = stress_xx
[../]
[./stress_xy]
type = PointValue
point = '0 0 0'
variable = stress_xy
[../]
[./stress_xz]
type = PointValue
point = '0 0 0'
variable = stress_xz
[../]
[./stress_yy]
type = PointValue
point = '0 0 0'
variable = stress_yy
[../]
[./stress_yz]
type = PointValue
point = '0 0 0'
variable = stress_yz
[../]
[./stress_zz]
type = PointValue
point = '0 0 0'
variable = stress_zz
[../]
[./strainp_xx]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xx
[../]
[./strainp_xy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xy
[../]
[./strainp_xz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xz
[../]
[./strainp_yy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yy
[../]
[./strainp_yz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yz
[../]
[./strainp_zz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_zz
[../]
[./straint_xx]
type = PointValue
point = '0 0 0'
variable = strain_xx
[../]
[./straint_xy]
type = PointValue
point = '0 0 0'
variable = strain_xy
[../]
[./straint_xz]
type = PointValue
point = '0 0 0'
variable = strain_xz
[../]
[./straint_yy]
type = PointValue
point = '0 0 0'
variable = strain_yy
[../]
[./straint_yz]
type = PointValue
point = '0 0 0'
variable = strain_yz
[../]
[./straint_zz]
type = PointValue
point = '0 0 0'
variable = strain_zz
[../]
[./f_shear]
type = PointValue
point = '0 0 0'
variable = f_shear
[../]
[./f_tensile]
type = PointValue
point = '0 0 0'
variable = f_tensile
[../]
[./f_compressive]
type = PointValue
point = '0 0 0'
variable = f_compressive
[../]
[./intnl_shear]
type = PointValue
point = '0 0 0'
variable = intnl_shear
[../]
[./intnl_tensile]
type = PointValue
point = '0 0 0'
variable = intnl_tensile
[../]
[./iter]
type = PointValue
point = '0 0 0'
variable = iter
[../]
[./ls]
type = PointValue
point = '0 0 0'
variable = ls
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningConstant
value = 80
[../]
[./tanphi]
type = SolidMechanicsHardeningConstant
value = 0.5
[../]
[./tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.1111077
[../]
[./t_strength]
type = SolidMechanicsHardeningConstant
value = 6
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 40
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '4 4'
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[../]
[./stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
tip_smoother = 0
smoothing_tol = 0
yield_function_tol = 1E-5
[../]
[]
[Executioner]
end_time = 2
dt = 1
type = Transient
[]
[Outputs]
file_base = small_deform3
csv = true
[]
(modules/solid_mechanics/test/tests/jacobian/cwp08.i)
# Capped weak-plane plasticity
# checking jacobian for shear + compression failure
[Mesh]
type = GeneratedMesh
dim = 3
[]
[GlobalParams]
block = 0
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Kernels]
[SolidMechanics]
displacements = 'disp_x disp_y disp_z'
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningExponential
value_0 = 1
value_residual = 1
rate = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningExponential
value_0 = 1.0
value_residual = 1.0
rate = 2
[../]
[./tanpsi]
type = SolidMechanicsHardeningExponential
value_0 = 0.1
value_residual = 0.1
rate = 1
[../]
[./t_strength]
type = SolidMechanicsHardeningExponential
value_0 = 100
value_residual = 100
rate = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 0
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
lambda = 0.0
shear_modulus = 2.0
[../]
[./strain]
type = ComputeIncrementalStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '0 0 1 0 0 -1 1 -1 0'
eigenstrain_name = ini_stress
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = mc
tangent_operator = nonlinear
[../]
[./mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 0
smoothing_tol = 2
yield_function_tol = 1E-10
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
[]
(modules/solid_mechanics/test/tests/jacobian/cwp09.i)
# Capped weak-plane plasticity
# checking jacobian for tensile failure with hardening
[Mesh]
type = GeneratedMesh
dim = 3
[]
[GlobalParams]
block = 0
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Kernels]
[SolidMechanics]
displacements = 'disp_x disp_y disp_z'
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningExponential
value_0 = 100
value_residual = 100
rate = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningExponential
value_0 = 1.0
value_residual = 1.0
rate = 2
[../]
[./tanpsi]
type = SolidMechanicsHardeningExponential
value_0 = 0.1
value_residual = 0.1
rate = 1
[../]
[./t_strength]
type = SolidMechanicsHardeningExponential
value_0 = 1
value_residual = 2
rate = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 100
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
lambda = 1.0
shear_modulus = 2.0
[../]
[./strain]
type = ComputeIncrementalStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '0 0 2 0 0 -1 2 -1 1.5'
eigenstrain_name = ini_stress
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = mc
tangent_operator = nonlinear
[../]
[./mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 0
smoothing_tol = 2
yield_function_tol = 1E-10
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
[]
(modules/solid_mechanics/test/tests/multiple_two_parameter_plasticity/dp_then_wp.i)
# Use ComputeMultipleInelasticStress with two inelastic models: CappedDruckerPrager and CappedWeakPlane.
# The relative_tolerance and absolute_tolerance parameters are set very large so that
# only one iteration is performed. This is the algorithm that FLAC uses to model
# jointed rocks, only Capped-Mohr-Coulomb is used instead of CappedDruckerPrager
#
# initial_stress = diag(1E3, 1E3, 1E3)
# The CappedDruckerPrager has tensile strength 3E2 and large cohesion,
# so the stress initially returns to diag(1E2, 1E2, 1E2)
# The CappedWeakPlane has tensile strength zero and large cohesion,
# so the stress returns to diag(1E2 - v/(1-v)*1E2, 1E2 - v/(1-v)*1E2, 0)
# where v=0.2 is the Poisson's ratio
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
eigenstrain_names = ini_stress
generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz'
[../]
[]
[BCs]
[./x]
type = FunctionDirichletBC
variable = disp_x
boundary = 'front back'
function = 0
[../]
[./y]
type = FunctionDirichletBC
variable = disp_y
boundary = 'front back'
function = 0
[../]
[./z]
type = FunctionDirichletBC
variable = disp_z
boundary = 'front back'
function = 0
[../]
[]
[AuxVariables]
[./yield_fcn_dp]
order = CONSTANT
family = MONOMIAL
[../]
[./yield_fcn_wp]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./yield_fcn_dp_auxk]
type = MaterialStdVectorAux
index = 1 # this is the tensile yield function - it should be zero
property = cdp_plastic_yield_function
variable = yield_fcn_dp
[../]
[./yield_fcn_wp_auxk]
type = MaterialStdVectorAux
index = 1 # this is the tensile yield function - it should be zero
property = cwp_plastic_yield_function
variable = yield_fcn_wp
[../]
[]
[Postprocessors]
[./s_xx]
type = PointValue
point = '0 0 0'
variable = stress_xx
[../]
[./s_xy]
type = PointValue
point = '0 0 0'
variable = stress_xy
[../]
[./s_xz]
type = PointValue
point = '0 0 0'
variable = stress_xz
[../]
[./s_yy]
type = PointValue
point = '0 0 0'
variable = stress_yy
[../]
[./s_yz]
type = PointValue
point = '0 0 0'
variable = stress_yz
[../]
[./s_zz]
type = PointValue
point = '0 0 0'
variable = stress_zz
[../]
[./f_dp]
type = PointValue
point = '0 0 0'
variable = yield_fcn_dp
[../]
[./f_wp]
type = PointValue
point = '0 0 0'
variable = yield_fcn_wp
[../]
[]
[UserObjects]
[./ts]
type = SolidMechanicsHardeningConstant
value = 300
[../]
[./cs]
type = SolidMechanicsHardeningConstant
value = 1E4
[../]
[./mc_coh]
type = SolidMechanicsHardeningConstant
value = 1E4
[../]
[./mc_phi]
type = SolidMechanicsHardeningConstant
value = 20
convert_to_radians = true
[../]
[./mc_psi]
type = SolidMechanicsHardeningConstant
value = 0
[../]
[./dp]
type = SolidMechanicsPlasticDruckerPrager
mc_cohesion = mc_coh
mc_friction_angle = mc_phi
mc_dilation_angle = mc_psi
internal_constraint_tolerance = 1 # irrelevant here
yield_function_tolerance = 1 # irrelevant here
[../]
[./wp_coh]
type = SolidMechanicsHardeningConstant
value = 1E4
[../]
[./wp_tanphi]
type = SolidMechanicsHardeningConstant
value = 0.5
[../]
[./wp_tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.1111077
[../]
[./wp_t_strength]
type = SolidMechanicsHardeningConstant
value = 0
[../]
[./wp_c_strength]
type = SolidMechanicsHardeningConstant
value = 1E4
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
poissons_ratio = 0.2
youngs_modulus = 1E7
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '1E3 0 0 0 1E3 0 0 0 1E3'
eigenstrain_name = ini_stress
[../]
[./admissible]
type = ComputeMultipleInelasticStress
relative_tolerance = 1E4
absolute_tolerance = 2
inelastic_models = 'cdp cwp'
perform_finite_strain_rotations = false
[../]
[./cdp]
type = CappedDruckerPragerStressUpdate
base_name = cdp
DP_model = dp
tensile_strength = ts
compressive_strength = cs
yield_function_tol = 1E-5
tip_smoother = 1E3
smoothing_tol = 1E3
[../]
[./cwp]
type = CappedWeakPlaneStressUpdate
base_name = cwp
cohesion = wp_coh
tan_friction_angle = wp_tanphi
tan_dilation_angle = wp_tanpsi
tensile_strength = wp_t_strength
compressive_strength = wp_c_strength
tip_smoother = 1E3
smoothing_tol = 1E3
yield_function_tol = 1E-5
[../]
[]
[Executioner]
end_time = 1
dt = 1
type = Transient
[]
[Outputs]
file_base = dp_then_wp
csv = true
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/except4.i)
# Exception: incorrect userobject types
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningConstant
value = -1
[../]
[./tanphi]
type = SolidMechanicsHardeningConstant
value = 0.5
[../]
[./tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.1111077
[../]
[./t_strength]
type = SolidMechanicsHardeningConstant
value = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 2
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0 1E6'
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[../]
[./stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-5
[../]
[]
[Executioner]
end_time = 1
dt = 1
type = Transient
[]
(modules/solid_mechanics/test/tests/jacobian/phe01.i)
# Capped weak-plane plasticity, Kernel = PlasticHeatEnergy
[Mesh]
type = GeneratedMesh
dim = 3
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Kernels]
[./silly_phe]
type = PlasticHeatEnergy
coeff = 0.5
variable = disp_x
[../]
[./dummy_disp_y]
type = TimeDerivative
variable = disp_y
[../]
[./dummy_disp_z]
type = TimeDerivative
variable = disp_z
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningExponential
value_0 = 1
value_residual = 2
rate = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningExponential
value_0 = 1.0
value_residual = 0.5
rate = 2
[../]
[./tanpsi]
type = SolidMechanicsHardeningExponential
value_0 = 0.1
value_residual = 0.05
rate = 3
[../]
[./t_strength]
type = SolidMechanicsHardeningExponential
value_0 = 100
value_residual = 100
rate = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningCubic
value_0 = 1
value_residual = 0
internal_0 = -2
internal_limit = 0
[../]
[]
[Materials]
[./phe]
type = ComputePlasticHeatEnergy
[../]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
lambda = 1.0
shear_modulus = 2.0
[../]
[./strain]
type = ComputeIncrementalStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '0 0 0 0 0 1 0 1 -1.5'
eigenstrain_name = ini_stress
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = mc
tangent_operator = nonlinear
[../]
[./mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-10
perfect_guess = false
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
[]
(modules/solid_mechanics/test/tests/jacobian/cwp11.i)
# Capped weak-plane plasticity
# checking jacobian for shear + tensile failure with hardening
[Mesh]
type = GeneratedMesh
dim = 3
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Kernels]
[SolidMechanics]
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningExponential
value_0 = 1
value_residual = 2
rate = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningExponential
value_0 = 1.0
value_residual = 0.5
rate = 2
[../]
[./tanpsi]
type = SolidMechanicsHardeningExponential
value_0 = 0.1
value_residual = 0.05
rate = 3
[../]
[./t_strength]
type = SolidMechanicsHardeningExponential
value_0 = 100
value_residual = 100
rate = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningCubic
value_0 = 1
value_residual = 0
internal_0 = -2
internal_limit = 0
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
lambda = 1.0
shear_modulus = 2.0
[../]
[./strain]
type = ComputeIncrementalStrain
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '0 0 0 0 0 1 0 1 -1.5'
eigenstrain_name = ini_stress
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = mc
tangent_operator = nonlinear
[../]
[./mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-10
perfect_guess = false
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
[]
(modules/solid_mechanics/test/tests/jacobian/cwp01.i)
# Capped weak-plane plasticity
# checking jacobian for a fully-elastic situation
[Mesh]
type = GeneratedMesh
dim = 3
[]
[GlobalParams]
block = 0
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[ICs]
[./disp_x]
type = RandomIC
variable = disp_x
min = -0.1
max = 0.1
[../]
[./disp_y]
type = RandomIC
variable = disp_y
min = -0.1
max = 0.1
[../]
[./disp_z]
type = RandomIC
variable = disp_z
min = -0.1
max = 0.1
[../]
[]
[Kernels]
[SolidMechanics]
displacements = 'disp_x disp_y disp_z'
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningExponential
value_0 = 100
value_residual = 2
rate = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningExponential
value_0 = 1.0
value_residual = 0.5
rate = 2
[../]
[./tanpsi]
type = SolidMechanicsHardeningExponential
value_0 = 0.1
value_residual = 0.05
rate = 1
[../]
[./t_strength]
type = SolidMechanicsHardeningExponential
value_0 = 100
value_residual = 0
rate = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 100
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
lambda = 1.0
shear_modulus = 2.0
[../]
[./strain]
type = ComputeIncrementalStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '1 2 3 2 -4 -5 3 -5 2'
eigenstrain_name = ini_stress
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = mc
[../]
[./mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 1
smoothing_tol = 2
yield_function_tol = 1E-10
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
[]
(modules/solid_mechanics/test/tests/jacobian/cwp04.i)
# Capped weak-plane plasticity
# checking jacobian for tensile failure, with some shear
[Mesh]
type = GeneratedMesh
dim = 3
[]
[GlobalParams]
block = 0
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Kernels]
[SolidMechanics]
displacements = 'disp_x disp_y disp_z'
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningExponential
value_0 = 100
value_residual = 2
rate = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningExponential
value_0 = 1.0
value_residual = 0.5
rate = 2
[../]
[./tanpsi]
type = SolidMechanicsHardeningExponential
value_0 = 0.1
value_residual = 0.05
rate = 1
[../]
[./t_strength]
type = SolidMechanicsHardeningExponential
value_0 = 1
value_residual = 1
rate = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 100
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
lambda = 1.0
shear_modulus = 2.0
[../]
[./strain]
type = ComputeIncrementalStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '0 0 0 0 0 1 0 1 2'
eigenstrain_name = ini_stress
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = mc
tangent_operator = nonlinear
[../]
[./mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 1
smoothing_tol = 2
yield_function_tol = 1E-10
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
#petsc_options = '-snes_test_display'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/small_deform5.i)
# Plastic deformation, shear failure
# With Young = 10, poisson=0.25 (Lame lambda=4, mu=4)
# applying the following
# deformation to the zmax surface of a unit cube:
# disp_x = 8*t
# disp_y = 6*t
# disp_z = 5*t/6
# should yield trial stress:
# stress_zz = 10*t
# stress_zx = 32*t
# stress_zy = 24*t (so q_trial = 40*t)
# Use tan(friction_angle) = 0.5 and tan(dilation_angle) = 1/6, and cohesion=20,
# the system should return to p=0, q=20, ie stress_zz=0, stress_xz=16,
# stress_yz=12 on the first time step (t=1)
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz'
[../]
[]
[BCs]
[./bottomx]
type = DirichletBC
variable = disp_x
boundary = back
value = 0.0
[../]
[./bottomy]
type = DirichletBC
variable = disp_y
boundary = back
value = 0.0
[../]
[./bottomz]
type = DirichletBC
variable = disp_z
boundary = back
value = 0.0
[../]
[./topx]
type = FunctionDirichletBC
variable = disp_x
boundary = front
function = 8*t
[../]
[./topy]
type = FunctionDirichletBC
variable = disp_y
boundary = front
function = 6*t
[../]
[./topz]
type = FunctionDirichletBC
variable = disp_z
boundary = front
function = 5*t/6
[../]
[]
[AuxVariables]
[./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]
[./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
[../]
[]
[Postprocessors]
[./stress_xx]
type = PointValue
point = '0 0 0'
variable = stress_xx
[../]
[./stress_xy]
type = PointValue
point = '0 0 0'
variable = stress_xy
[../]
[./stress_xz]
type = PointValue
point = '0 0 0'
variable = stress_xz
[../]
[./stress_yy]
type = PointValue
point = '0 0 0'
variable = stress_yy
[../]
[./stress_yz]
type = PointValue
point = '0 0 0'
variable = stress_yz
[../]
[./stress_zz]
type = PointValue
point = '0 0 0'
variable = stress_zz
[../]
[./strainp_xx]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xx
[../]
[./strainp_xy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xy
[../]
[./strainp_xz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xz
[../]
[./strainp_yy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yy
[../]
[./strainp_yz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yz
[../]
[./strainp_zz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_zz
[../]
[./straint_xx]
type = PointValue
point = '0 0 0'
variable = strain_xx
[../]
[./straint_xy]
type = PointValue
point = '0 0 0'
variable = strain_xy
[../]
[./straint_xz]
type = PointValue
point = '0 0 0'
variable = strain_xz
[../]
[./straint_yy]
type = PointValue
point = '0 0 0'
variable = strain_yy
[../]
[./straint_yz]
type = PointValue
point = '0 0 0'
variable = strain_yz
[../]
[./straint_zz]
type = PointValue
point = '0 0 0'
variable = strain_zz
[../]
[./f_shear]
type = PointValue
point = '0 0 0'
variable = f_shear
[../]
[./f_tensile]
type = PointValue
point = '0 0 0'
variable = f_tensile
[../]
[./f_compressive]
type = PointValue
point = '0 0 0'
variable = f_compressive
[../]
[./intnl_shear]
type = PointValue
point = '0 0 0'
variable = intnl_shear
[../]
[./intnl_tensile]
type = PointValue
point = '0 0 0'
variable = intnl_tensile
[../]
[./iter]
type = PointValue
point = '0 0 0'
variable = iter
[../]
[./ls]
type = PointValue
point = '0 0 0'
variable = ls
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningConstant
value = 20
[../]
[./tanphi]
type = SolidMechanicsHardeningConstant
value = 0.5
[../]
[./tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.166666666667
[../]
[./t_strength]
type = SolidMechanicsHardeningConstant
value = 100
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 100
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '4 4'
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[../]
[./stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
tip_smoother = 0
smoothing_tol = 0
yield_function_tol = 1E-5
[../]
[]
[Executioner]
end_time = 1
dt = 1
type = Transient
[]
[Outputs]
file_base = small_deform5
csv = true
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/except3.i)
# Exception: incorrect userobject types
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningConstant
value = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningConstant
value = 0.05
[../]
[./tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.1111077
[../]
[./t_strength]
type = SolidMechanicsHardeningConstant
value = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 2
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0 1E6'
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[../]
[./stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-5
[../]
[]
[Executioner]
end_time = 1
dt = 1
type = Transient
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/small_deform8.i)
# Plastic deformation, compression with hardening
# With Lame lambda=0 and Lame mu=1, applying the following
# deformation to the zmax surface of a unit cube:
# disp_z = -t
# should yield trial stress:
# stress_zz = -2*t
# The compressive strength varies as a cubic between 1 (at intnl=0)
# and 2 (at intnl=1). The equation to solve is
# 2 - Ezzzz * ga = -2 * (ga - 1/2)^3 + (3/2) (ga - 1/2) + 3/2
# where the left-hand side comes from p = p_trial + ga * Ezzzz
# and the right-hand side is the cubic compressive strength
# The solution is ga = 0.355416 ( = intnl[1]), and the cubic
# is 1.289168 ( = -p) at that point
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz'
[../]
[]
[BCs]
[./bottomx]
type = DirichletBC
variable = disp_x
boundary = back
value = 0.0
[../]
[./bottomy]
type = DirichletBC
variable = disp_y
boundary = back
value = 0.0
[../]
[./bottomz]
type = DirichletBC
variable = disp_z
boundary = back
value = 0.0
[../]
[./topx]
type = FunctionDirichletBC
variable = disp_x
boundary = front
function = 0
[../]
[./topy]
type = FunctionDirichletBC
variable = disp_y
boundary = front
function = 0
[../]
[./topz]
type = FunctionDirichletBC
variable = disp_z
boundary = front
function = -t
[../]
[]
[AuxVariables]
[./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]
[./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
[../]
[]
[Postprocessors]
[./stress_xx]
type = PointValue
point = '0 0 0'
variable = stress_xx
[../]
[./stress_xy]
type = PointValue
point = '0 0 0'
variable = stress_xy
[../]
[./stress_xz]
type = PointValue
point = '0 0 0'
variable = stress_xz
[../]
[./stress_yy]
type = PointValue
point = '0 0 0'
variable = stress_yy
[../]
[./stress_yz]
type = PointValue
point = '0 0 0'
variable = stress_yz
[../]
[./stress_zz]
type = PointValue
point = '0 0 0'
variable = stress_zz
[../]
[./strainp_xx]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xx
[../]
[./strainp_xy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xy
[../]
[./strainp_xz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xz
[../]
[./strainp_yy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yy
[../]
[./strainp_yz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yz
[../]
[./strainp_zz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_zz
[../]
[./straint_xx]
type = PointValue
point = '0 0 0'
variable = strain_xx
[../]
[./straint_xy]
type = PointValue
point = '0 0 0'
variable = strain_xy
[../]
[./straint_xz]
type = PointValue
point = '0 0 0'
variable = strain_xz
[../]
[./straint_yy]
type = PointValue
point = '0 0 0'
variable = strain_yy
[../]
[./straint_yz]
type = PointValue
point = '0 0 0'
variable = strain_yz
[../]
[./straint_zz]
type = PointValue
point = '0 0 0'
variable = strain_zz
[../]
[./f_shear]
type = PointValue
point = '0 0 0'
variable = f_shear
[../]
[./f_tensile]
type = PointValue
point = '0 0 0'
variable = f_tensile
[../]
[./f_compressive]
type = PointValue
point = '0 0 0'
variable = f_compressive
[../]
[./intnl_shear]
type = PointValue
point = '0 0 0'
variable = intnl_shear
[../]
[./intnl_tensile]
type = PointValue
point = '0 0 0'
variable = intnl_tensile
[../]
[./iter]
type = PointValue
point = '0 0 0'
variable = iter
[../]
[./ls]
type = PointValue
point = '0 0 0'
variable = ls
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningConstant
value = 20
[../]
[./tanphi]
type = SolidMechanicsHardeningConstant
value = 0.5
[../]
[./tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.1
[../]
[./t_strength]
type = SolidMechanicsHardeningConstant
value = 100
[../]
[./c_strength]
type = SolidMechanicsHardeningCubic
value_0 = 2
value_residual = 1
internal_0 = -1
internal_limit = 0
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0 1'
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[../]
[./stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 5
smoothing_tol = 5
yield_function_tol = 1E-10
[../]
[]
[Executioner]
end_time = 1
dt = 1
type = Transient
[]
[Outputs]
file_base = small_deform8
csv = true
[]
(modules/solid_mechanics/test/tests/jacobian/cwp05.i)
# Capped weak-plane plasticity
# checking jacobian for shear failure
[Mesh]
type = GeneratedMesh
dim = 3
[]
[GlobalParams]
block = 0
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Kernels]
[SolidMechanics]
displacements = 'disp_x disp_y disp_z'
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningExponential
value_0 = 1
value_residual = 1
rate = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningExponential
value_0 = 1.0
value_residual = 1.0
rate = 2
[../]
[./tanpsi]
type = SolidMechanicsHardeningExponential
value_0 = 0.1
value_residual = 0.1
rate = 1
[../]
[./t_strength]
type = SolidMechanicsHardeningExponential
value_0 = 100
value_residual = 100
rate = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 100
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
lambda = 1.0
shear_modulus = 2.0
[../]
[./strain]
type = ComputeIncrementalStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '0 0 1 0 0 10 1 10 0'
eigenstrain_name = ini_stress
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = mc
tangent_operator = nonlinear
[../]
[./mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 0
smoothing_tol = 2
yield_function_tol = 1E-10
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
#petsc_options = '-snes_test_display'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/small_deform7.i)
# Plastic deformation, tensile with hardening
# With Lame lambda=0 and Lame mu=1, applying the following
# deformation to the zmax surface of a unit cube:
# disp_z = t
# should yield trial stress:
# stress_zz = 2*t
# The tensile strength varies as a cubic between 1 (at intnl=0)
# and 2 (at intnl=1). The equation to solve is
# 2 - Ezzzz * ga = -2 * (ga - 1/2)^3 + (3/2) (ga - 1/2) + 3/2
# where the left-hand side comes from p = p_trial - ga * Ezzzz
# and the right-hand side is the cubic tensile strength
# The solution is ga = 0.355416 ( = intnl[1]), and the cubic
# is 1.289168 ( = p) at that point
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz'
[../]
[]
[BCs]
[./bottomx]
type = DirichletBC
variable = disp_x
boundary = back
value = 0.0
[../]
[./bottomy]
type = DirichletBC
variable = disp_y
boundary = back
value = 0.0
[../]
[./bottomz]
type = DirichletBC
variable = disp_z
boundary = back
value = 0.0
[../]
[./topx]
type = FunctionDirichletBC
variable = disp_x
boundary = front
function = 0
[../]
[./topy]
type = FunctionDirichletBC
variable = disp_y
boundary = front
function = 0
[../]
[./topz]
type = FunctionDirichletBC
variable = disp_z
boundary = front
function = t
[../]
[]
[AuxVariables]
[./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]
[./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
[../]
[]
[Postprocessors]
[./stress_xx]
type = PointValue
point = '0 0 0'
variable = stress_xx
[../]
[./stress_xy]
type = PointValue
point = '0 0 0'
variable = stress_xy
[../]
[./stress_xz]
type = PointValue
point = '0 0 0'
variable = stress_xz
[../]
[./stress_yy]
type = PointValue
point = '0 0 0'
variable = stress_yy
[../]
[./stress_yz]
type = PointValue
point = '0 0 0'
variable = stress_yz
[../]
[./stress_zz]
type = PointValue
point = '0 0 0'
variable = stress_zz
[../]
[./strainp_xx]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xx
[../]
[./strainp_xy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xy
[../]
[./strainp_xz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xz
[../]
[./strainp_yy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yy
[../]
[./strainp_yz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yz
[../]
[./strainp_zz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_zz
[../]
[./straint_xx]
type = PointValue
point = '0 0 0'
variable = strain_xx
[../]
[./straint_xy]
type = PointValue
point = '0 0 0'
variable = strain_xy
[../]
[./straint_xz]
type = PointValue
point = '0 0 0'
variable = strain_xz
[../]
[./straint_yy]
type = PointValue
point = '0 0 0'
variable = strain_yy
[../]
[./straint_yz]
type = PointValue
point = '0 0 0'
variable = strain_yz
[../]
[./straint_zz]
type = PointValue
point = '0 0 0'
variable = strain_zz
[../]
[./f_shear]
type = PointValue
point = '0 0 0'
variable = f_shear
[../]
[./f_tensile]
type = PointValue
point = '0 0 0'
variable = f_tensile
[../]
[./f_compressive]
type = PointValue
point = '0 0 0'
variable = f_compressive
[../]
[./intnl_shear]
type = PointValue
point = '0 0 0'
variable = intnl_shear
[../]
[./intnl_tensile]
type = PointValue
point = '0 0 0'
variable = intnl_tensile
[../]
[./iter]
type = PointValue
point = '0 0 0'
variable = iter
[../]
[./ls]
type = PointValue
point = '0 0 0'
variable = ls
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningConstant
value = 20
[../]
[./tanphi]
type = SolidMechanicsHardeningConstant
value = 0.5
[../]
[./tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.1
[../]
[./t_strength]
type = SolidMechanicsHardeningCubic
value_0 = 1
value_residual = 2
internal_limit = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 100
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0 1'
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[../]
[./stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 5
smoothing_tol = 5
yield_function_tol = 1E-10
[../]
[]
[Executioner]
end_time = 1
dt = 1
type = Transient
[]
[Outputs]
file_base = small_deform7
csv = true
[]
(modules/solid_mechanics/test/tests/jacobian/cwp06.i)
# Capped weak-plane plasticity
# checking jacobian for shear failure, with smoothing
[Mesh]
type = GeneratedMesh
dim = 3
[]
[GlobalParams]
block = 0
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Kernels]
[SolidMechanics]
displacements = 'disp_x disp_y disp_z'
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningExponential
value_0 = 1
value_residual = 1
rate = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningExponential
value_0 = 1.0
value_residual = 1.0
rate = 2
[../]
[./tanpsi]
type = SolidMechanicsHardeningExponential
value_0 = 0.1
value_residual = 0.1
rate = 1
[../]
[./t_strength]
type = SolidMechanicsHardeningExponential
value_0 = 100
value_residual = 100
rate = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 100
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
lambda = 1.0
shear_modulus = 2.0
[../]
[./strain]
type = ComputeIncrementalStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '0 0 1 0 0 -1 1 -1 0'
eigenstrain_name = ini_stress
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = mc
tangent_operator = nonlinear
[../]
[./mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 1
smoothing_tol = 2
yield_function_tol = 1E-10
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/except2.i)
# Exception: incorrect userobject types
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningConstant
value = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningConstant
value = 0.5
[../]
[./tanpsi]
type = SolidMechanicsHardeningConstant
value = -0.1111077
[../]
[./t_strength]
type = SolidMechanicsHardeningConstant
value = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 2
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0 1E6'
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[../]
[./stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-5
[../]
[]
[Executioner]
end_time = 1
dt = 1
type = Transient
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/except1.i)
# Exception: incorrect userobject types
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningConstant
value = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningConstant
value = -0.5
[../]
[./tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.1111077
[../]
[./t_strength]
type = SolidMechanicsHardeningConstant
value = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 2
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '0 1E6'
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[../]
[./stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-5
[../]
[]
[Executioner]
end_time = 1
dt = 1
type = Transient
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/pull_and_shear.i)
# Dynamic problem with plasticity.
# A column of material (not subject to gravity) has the z-displacement
# of its sides fixed, but the centre of its bottom side is pulled
# downwards. This causes failure in the bottom elements.
#
# The problem utilises damping in the following way.
# The DynamicStressDivergenceTensors forms the residual
# integral grad(stress) + zeta*grad(stress-dot)
# = V/L * elasticity * (du/dx + zeta * dv/dx)
# where V is the elemental volume, and L is the length-scale,
# and u is the displacement, and v is the velocity.
# The InertialForce forms the residual
# integral density * (accel + eta * velocity)
# = V * density * (a + eta * v)
# where a is the acceleration.
# So, a damped oscillator description with both these
# kernels looks like
# 0 = V * (density * a + density * eta * v + elasticity * zeta * v / L^2 + elasticity / L^2 * u)
# Critical damping is when the coefficient of v is
# 2 * sqrt(density * elasticity / L^2)
# In the case at hand, density=1E4, elasticity~1E10 (Young is 16GPa),
# L~1 to 10 (in the horizontal or vertical direction), so this coefficient ~ 1E7 to 1E6.
# Choosing eta = 1E3 and zeta = 1E-2 gives approximate critical damping.
# If zeta is high then steady-state is achieved very quickly.
#
# In the case of plasticity, the effective stiffness of the elements
# is significantly less. Therefore, the above parameters give
# overdamping.
#
# This simulation is a nice example of the irreversable and non-uniqueness
# of simulations involving plasticity. The result depends on the damping
# parameters and the time stepping.
[Mesh]
[generated_mesh]
type = GeneratedMeshGenerator
dim = 3
nx = 10
ny = 1
nz = 5
bias_z = 1.5
xmin = -10
xmax = 10
ymin = -10
ymax = 10
zmin = -100
zmax = 0
[]
[bottomz_middle]
type = BoundingBoxNodeSetGenerator
new_boundary = bottomz_middle
bottom_left = '-1 -1500 -105'
top_right = '1 1500 -95'
input = generated_mesh
[]
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
beta = 0.25 # Newmark time integration
gamma = 0.5 # Newmark time integration
eta = 1E3 #0.3E4 # higher values mean more damping via density
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[]
[Kernels]
[DynamicSolidMechanics] # zeta*K*vel + K * disp
stiffness_damping_coefficient = 1E-2 # higher values mean more damping via stiffness
hht_alpha = 0 # better nonlinear convergence than for alpha>0
[]
[inertia_x] # M*accel + eta*M*vel
type = InertialForce
use_displaced_mesh = false
variable = disp_x
velocity = vel_x
acceleration = accel_x
[]
[inertia_y]
type = InertialForce
use_displaced_mesh = false
variable = disp_y
velocity = vel_y
acceleration = accel_y
[]
[inertia_z]
type = InertialForce
use_displaced_mesh = false
variable = disp_z
velocity = vel_z
acceleration = accel_z
[]
[]
[BCs]
[no_x2]
type = DirichletBC
variable = disp_x
boundary = right
value = 0.0
[]
[no_x1]
type = DirichletBC
variable = disp_x
boundary = left
value = 0.0
[]
[no_y1]
type = DirichletBC
variable = disp_y
boundary = bottom
value = 0.0
[]
[no_y2]
type = DirichletBC
variable = disp_y
boundary = top
value = 0.0
[]
[z_fixed_sides_xmin]
type = DirichletBC
variable = disp_z
boundary = left
value = 0
[]
[z_fixed_sides_xmax]
type = DirichletBC
variable = disp_z
boundary = right
value = 0
[]
[bottomz]
type = FunctionDirichletBC
variable = disp_z
boundary = bottomz_middle
function = max(-10*t,-10)
[]
[]
[AuxVariables]
[accel_x]
[]
[vel_x]
[]
[accel_y]
[]
[vel_y]
[]
[accel_z]
[]
[vel_z]
[]
[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]
[accel_x] # Calculates and stores acceleration at the end of time step
type = NewmarkAccelAux
variable = accel_x
displacement = disp_x
velocity = vel_x
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
execute_on = timestep_end
[]
[accel_y]
type = NewmarkAccelAux
variable = accel_y
displacement = disp_y
velocity = vel_y
execute_on = timestep_end
[]
[vel_y]
type = NewmarkVelAux
variable = vel_y
acceleration = accel_y
execute_on = timestep_end
[]
[accel_z]
type = NewmarkAccelAux
variable = accel_z
displacement = disp_z
velocity = vel_z
execute_on = timestep_end
[]
[vel_z]
type = NewmarkVelAux
variable = vel_z
acceleration = accel_z
execute_on = timestep_end
[]
[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]
type = SolidMechanicsHardeningConstant
value = 1E6
[]
[tanphi]
type = SolidMechanicsHardeningConstant
value = 0.5
[]
[tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.166666666667
[]
[t_strength]
type = SolidMechanicsHardeningConstant
value = 0
[]
[c_strength]
type = SolidMechanicsHardeningConstant
value = 1E80
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '6.4E9 6.4E9' # young 16MPa, Poisson 0.25
[]
[strain]
type = ComputeIncrementalStrain
[]
[admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[]
[stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
tip_smoother = 1E6
smoothing_tol = 0.5E6
yield_function_tol = 1E-2
[]
[density]
type = GenericConstantMaterial
block = 0
prop_names = density
prop_values = 1E4
[]
[]
[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 = 1E1
nl_rel_tol = 1e-5
l_tol = 1E-10
l_max_its = 100
nl_max_its = 100
num_steps = 8
dt = 0.1
type = Transient
[]
[Outputs]
file_base = pull_and_shear
exodus = true
csv = true
[]
(modules/solid_mechanics/test/tests/capped_weak_plane/small_deform10.i)
# apply a shear deformation and tensile stretch to observe all hardening.
# Here p_trial=12, q_trial=2*Sqrt(20)
# MOOSE yields:
# q_returned = 1.696
# p_returned = 0.100
# intnl_shear = 1.81
# intnl_tens = 0.886
# These give, at the returned point
# cohesion = 1.84
# tanphi = 0.513
# tanpsi = 0.058
# tensile = 0.412
# This means that
# f_shear = -0.0895
# f_tensile = -0.312
# Note that these are within smoothing_tol (=1) of each other
# Hence, smoothing must be used:
# ismoother = 0.0895
# (which gives the yield function value = 0)
# smoother = 0.328
# This latter gives dg/dq = 0.671, dg/dp = 0.368
# for the flow directions. Finally ga = 2.70, and
# the returned point satisfies the normality conditions.
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz'
[../]
[]
[BCs]
[./bottomx]
type = DirichletBC
variable = disp_x
boundary = back
value = 0.0
[../]
[./bottomy]
type = DirichletBC
variable = disp_y
boundary = back
value = 0.0
[../]
[./bottomz]
type = DirichletBC
variable = disp_z
boundary = back
value = 0.0
[../]
[./topx]
type = FunctionDirichletBC
variable = disp_x
boundary = front
function = 't'
[../]
[./topy]
type = FunctionDirichletBC
variable = disp_y
boundary = front
function = '2*t'
[../]
[./topz]
type = FunctionDirichletBC
variable = disp_z
boundary = front
function = 't'
[../]
[]
[AuxVariables]
[./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]
[./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
[../]
[]
[Postprocessors]
[./stress_xx]
type = PointValue
point = '0 0 0'
variable = stress_xx
[../]
[./stress_xy]
type = PointValue
point = '0 0 0'
variable = stress_xy
[../]
[./stress_xz]
type = PointValue
point = '0 0 0'
variable = stress_xz
[../]
[./stress_yy]
type = PointValue
point = '0 0 0'
variable = stress_yy
[../]
[./stress_yz]
type = PointValue
point = '0 0 0'
variable = stress_yz
[../]
[./stress_zz]
type = PointValue
point = '0 0 0'
variable = stress_zz
[../]
[./strainp_xx]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xx
[../]
[./strainp_xy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xy
[../]
[./strainp_xz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_xz
[../]
[./strainp_yy]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yy
[../]
[./strainp_yz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_yz
[../]
[./strainp_zz]
type = PointValue
point = '0 0 0'
variable = plastic_strain_zz
[../]
[./straint_xx]
type = PointValue
point = '0 0 0'
variable = strain_xx
[../]
[./straint_xy]
type = PointValue
point = '0 0 0'
variable = strain_xy
[../]
[./straint_xz]
type = PointValue
point = '0 0 0'
variable = strain_xz
[../]
[./straint_yy]
type = PointValue
point = '0 0 0'
variable = strain_yy
[../]
[./straint_yz]
type = PointValue
point = '0 0 0'
variable = strain_yz
[../]
[./straint_zz]
type = PointValue
point = '0 0 0'
variable = strain_zz
[../]
[./f_shear]
type = PointValue
point = '0 0 0'
variable = f_shear
[../]
[./f_tensile]
type = PointValue
point = '0 0 0'
variable = f_tensile
[../]
[./f_compressive]
type = PointValue
point = '0 0 0'
variable = f_compressive
[../]
[./intnl_shear]
type = PointValue
point = '0 0 0'
variable = intnl_shear
[../]
[./intnl_tensile]
type = PointValue
point = '0 0 0'
variable = intnl_tensile
[../]
[./iter]
type = PointValue
point = '0 0 0'
variable = iter
[../]
[./ls]
type = PointValue
point = '0 0 0'
variable = ls
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningExponential
value_0 = 1
value_residual = 2
rate = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningExponential
value_0 = 1.0
value_residual = 0.5
rate = 2
[../]
[./tanpsi]
type = SolidMechanicsHardeningExponential
value_0 = 0.1
value_residual = 0.05
rate = 1
[../]
[./t_strength]
type = SolidMechanicsHardeningExponential
value_0 = 1
value_residual = 0
rate = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 1E8
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
fill_method = symmetric_isotropic
C_ijkl = '4 4'
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = stress
perform_finite_strain_rotations = false
[../]
[./stress]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 0
smoothing_tol = 1
yield_function_tol = 1E-3
perfect_guess = false
[../]
[]
[Executioner]
end_time = 1
dt = 1
type = Transient
[]
[Outputs]
file_base = small_deform10
[./csv]
type = CSV
[../]
[]
(modules/solid_mechanics/test/tests/multiple_two_parameter_plasticity/dp_and_wp.i)
# Use ComputeMultipleInelasticStress with two inelastic models: CappedDruckerPrager and CappedWeakPlane.
# The relative_tolerance and absolute_tolerance parameters are set small so that many
# Picard iterations need to be performed.
#
# The CappedDruckerPrager has tensile strength 3E2 and large cohesion,
# and the return-map sets stress = trial_stress - diag(d, d, d), for
# some d to be determined
# The CappedWeakPlane has tensile strength zero and large cohesion,
# and the return-map sets stress = diag(t - v*w/(1-v), t - v*w/(1-v), t - w)
# where t is trial stress, v is Poisson's ratio, and w is to be determined
#
# d and w are determined by demanding that the final stress shouldn't depend
# on the order of return-mapping (DP first then WP, or WP first then DP).
#
# Let the initial_stress = diag(I, I, I).
# The returned stress is diag(I - d - v*w/(1-v), I - d - v*w/(1-v), I - d - w). This
# must obey Tr(stress) <= dp_tensile_strength, and I-d-w <= wp_tensile_strength.
#
# For I = 1E3, and v = 0.2, the solution is d = 800 and w = 200, with
# stress = diag(150, 150, 0)
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Physics/SolidMechanics/QuasiStatic]
[./all]
add_variables = true
incremental = true
eigenstrain_names = ini_stress
generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz'
[../]
[]
[BCs]
[./x]
type = FunctionDirichletBC
variable = disp_x
boundary = 'front back'
function = 0
[../]
[./y]
type = FunctionDirichletBC
variable = disp_y
boundary = 'front back'
function = 0
[../]
[./z]
type = FunctionDirichletBC
variable = disp_z
boundary = 'front back'
function = 0
[../]
[]
[AuxVariables]
[./yield_fcn_dp]
order = CONSTANT
family = MONOMIAL
[../]
[./yield_fcn_wp]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./yield_fcn_dp_auxk]
type = MaterialStdVectorAux
index = 1 # this is the tensile yield function - it should be zero
property = cdp_plastic_yield_function
variable = yield_fcn_dp
[../]
[./yield_fcn_wp_auxk]
type = MaterialStdVectorAux
index = 1 # this is the tensile yield function - it should be zero
property = cwp_plastic_yield_function
variable = yield_fcn_wp
[../]
[]
[Postprocessors]
[./s_xx]
type = PointValue
point = '0 0 0'
variable = stress_xx
[../]
[./s_xy]
type = PointValue
point = '0 0 0'
variable = stress_xy
[../]
[./s_xz]
type = PointValue
point = '0 0 0'
variable = stress_xz
[../]
[./s_yy]
type = PointValue
point = '0 0 0'
variable = stress_yy
[../]
[./s_yz]
type = PointValue
point = '0 0 0'
variable = stress_yz
[../]
[./s_zz]
type = PointValue
point = '0 0 0'
variable = stress_zz
[../]
[./f_dp]
type = PointValue
point = '0 0 0'
variable = yield_fcn_dp
[../]
[./f_wp]
type = PointValue
point = '0 0 0'
variable = yield_fcn_wp
[../]
[]
[UserObjects]
[./ts]
type = SolidMechanicsHardeningConstant
value = 300
[../]
[./cs]
type = SolidMechanicsHardeningConstant
value = 1E4
[../]
[./mc_coh]
type = SolidMechanicsHardeningConstant
value = 1E4
[../]
[./mc_phi]
type = SolidMechanicsHardeningConstant
value = 20
convert_to_radians = true
[../]
[./mc_psi]
type = SolidMechanicsHardeningConstant
value = 0
[../]
[./dp]
type = SolidMechanicsPlasticDruckerPrager
mc_cohesion = mc_coh
mc_friction_angle = mc_phi
mc_dilation_angle = mc_psi
internal_constraint_tolerance = 1 # irrelevant here
yield_function_tolerance = 1 # irrelevant here
[../]
[./wp_coh]
type = SolidMechanicsHardeningConstant
value = 1E4
[../]
[./wp_tanphi]
type = SolidMechanicsHardeningConstant
value = 0.5
[../]
[./wp_tanpsi]
type = SolidMechanicsHardeningConstant
value = 0.1111077
[../]
[./wp_t_strength]
type = SolidMechanicsHardeningConstant
value = 0
[../]
[./wp_c_strength]
type = SolidMechanicsHardeningConstant
value = 1E4
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
poissons_ratio = 0.2
youngs_modulus = 1E7
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '1E3 0 0 0 1E3 0 0 0 1E3'
eigenstrain_name = ini_stress
[../]
[./admissible]
type = ComputeMultipleInelasticStress
relative_tolerance = 1E-8
inelastic_models = 'cdp cwp'
perform_finite_strain_rotations = false
[../]
[./cdp]
type = CappedDruckerPragerStressUpdate
base_name = cdp
DP_model = dp
tensile_strength = ts
compressive_strength = cs
yield_function_tol = 1E-5
tip_smoother = 1E3
smoothing_tol = 1E3
[../]
[./cwp]
type = CappedWeakPlaneStressUpdate
base_name = cwp
cohesion = wp_coh
tan_friction_angle = wp_tanphi
tan_dilation_angle = wp_tanpsi
tensile_strength = wp_t_strength
compressive_strength = wp_c_strength
tip_smoother = 1E3
smoothing_tol = 1E3
yield_function_tol = 1E-5
[../]
[]
[Executioner]
end_time = 1
dt = 1
type = Transient
[]
[Outputs]
file_base = dp_and_wp
csv = true
[]
(modules/solid_mechanics/test/tests/jacobian/cwp02.i)
# Capped weak-plane plasticity
# checking jacobian for tensile failure
[Mesh]
type = GeneratedMesh
dim = 3
[]
[GlobalParams]
block = 0
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Kernels]
[SolidMechanics]
displacements = 'disp_x disp_y disp_z'
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningExponential
value_0 = 100
value_residual = 2
rate = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningExponential
value_0 = 1.0
value_residual = 0.5
rate = 2
[../]
[./tanpsi]
type = SolidMechanicsHardeningExponential
value_0 = 0.1
value_residual = 0.05
rate = 1
[../]
[./t_strength]
type = SolidMechanicsHardeningExponential
value_0 = 1
value_residual = 1
rate = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 100
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
lambda = 1.0
shear_modulus = 2.0
[../]
[./strain]
type = ComputeIncrementalStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '0 0 0 0 0 0 0 0 2'
eigenstrain_name = ini_stress
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = mc
tangent_operator = nonlinear
[../]
[./mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 1
smoothing_tol = 2
yield_function_tol = 1E-10
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
#petsc_options = '-snes_test_display'
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
[]
(modules/solid_mechanics/test/tests/jacobian/cwp07.i)
# Capped weak-plane plasticity
# checking jacobian for shear + tensile failure
[Mesh]
type = GeneratedMesh
dim = 3
[]
[GlobalParams]
block = 0
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Kernels]
[SolidMechanics]
displacements = 'disp_x disp_y disp_z'
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningExponential
value_0 = 1
value_residual = 1
rate = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningExponential
value_0 = 1.0
value_residual = 1.0
rate = 2
[../]
[./tanpsi]
type = SolidMechanicsHardeningExponential
value_0 = 0.1
value_residual = 0.1
rate = 1
[../]
[./t_strength]
type = SolidMechanicsHardeningExponential
value_0 = 0
value_residual = 0
rate = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 100
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
lambda = 1.0
shear_modulus = 2.0
[../]
[./strain]
type = ComputeIncrementalStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '0 0 1 0 0 -1 1 -1 1'
eigenstrain_name = ini_stress
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = mc
tangent_operator = nonlinear
[../]
[./mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 0
smoothing_tol = 2
yield_function_tol = 1E-10
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
[]
(modules/solid_mechanics/test/tests/jacobian/cwp10.i)
# Capped weak-plane plasticity
# checking jacobian for shear failure with hardening
[Mesh]
type = GeneratedMesh
dim = 3
[]
[GlobalParams]
block = 0
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]
[Kernels]
[SolidMechanics]
displacements = 'disp_x disp_y disp_z'
[../]
[]
[UserObjects]
[./coh]
type = SolidMechanicsHardeningExponential
value_0 = 1
value_residual = 2
rate = 1
[../]
[./tanphi]
type = SolidMechanicsHardeningExponential
value_0 = 1.0
value_residual = 0.5
rate = 2
[../]
[./tanpsi]
type = SolidMechanicsHardeningExponential
value_0 = 0.1
value_residual = 0.05
rate = 3
[../]
[./t_strength]
type = SolidMechanicsHardeningExponential
value_0 = 100
value_residual = 100
rate = 1
[../]
[./c_strength]
type = SolidMechanicsHardeningConstant
value = 100
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
lambda = 1.0
shear_modulus = 2.0
[../]
[./strain]
type = ComputeIncrementalStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = ini_stress
[../]
[./ini_stress]
type = ComputeEigenstrainFromInitialStress
initial_stress = '0 0 2 0 0 -1 2 -1 0.1'
eigenstrain_name = ini_stress
[../]
[./admissible]
type = ComputeMultipleInelasticStress
inelastic_models = mc
tangent_operator = nonlinear
[../]
[./mc]
type = CappedWeakPlaneStressUpdate
cohesion = coh
tan_friction_angle = tanphi
tan_dilation_angle = tanpsi
tensile_strength = t_strength
compressive_strength = c_strength
max_NR_iterations = 20
tip_smoother = 0
smoothing_tol = 2
yield_function_tol = 1E-10
[../]
[]
[Preconditioning]
[./andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it -snes_type'
petsc_options_value = 'bcgs bjacobi 1E-15 1E-10 10000 test'
[../]
[]
[Executioner]
type = Transient
solve_type = Newton
[]
(modules/solid_mechanics/include/materials/CappedWeakInclinedPlaneStressUpdate.h)
// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html
#pragma once
#include "CappedWeakPlaneStressUpdate.h"
/**
* CappedWeakInclinedPlaneStressUpdate performs the return-map
* algorithm and associated stress updates for plastic
* models that describe capped weak-plane plasticity
*
* It assumes various things about the elasticity tensor, viz
* in the frame where the weak-plane's normal direction is the
* "2" direction:
* E(i,i,j,k) = 0 except if k=j
* E(0,0,i,j) = E(1,1,i,j)
*/
class CappedWeakInclinedPlaneStressUpdate : public CappedWeakPlaneStressUpdate
{
public:
static InputParameters validParams();
CappedWeakInclinedPlaneStressUpdate(const InputParameters & parameters);
/**
* Does the model require the elasticity tensor to be isotropic?
*/
bool requiresIsotropicTensor() override { return false; }
protected:
virtual void initQpStatefulProperties() override;
/// User-input value of the normal vector to the weak plane
RealVectorValue _n_input;
/// Current value of the normal
MaterialProperty<RealVectorValue> & _n;
/// Old value of the normal
const MaterialProperty<RealVectorValue> & _n_old;
/// Rotation matrix that rotates _n to "z"
RealTensorValue _rot_n_to_z;
/// Rotation matrix that rotates "z" to _n
RealTensorValue _rot_z_to_n;
/// Trial stress rotated to the frame where _n points along "z"
RankTwoTensor _rotated_trial;
/// Elasticity tensor rotated to the frame where _n points along "z"
RankFourTensor _rotated_Eijkl;
virtual void initializeReturnProcess() override;
virtual void finalizeReturnProcess(const RankTwoTensor & rotation_increment) override;
virtual void preReturnMap(Real p_trial,
Real q_trial,
const RankTwoTensor & stress_trial,
const std::vector<Real> & intnl_old,
const std::vector<Real> & yf,
const RankFourTensor & Eijkl) override;
virtual void computePQ(const RankTwoTensor & stress, Real & p, Real & q) const override;
virtual void setEppEqq(const RankFourTensor & Eijkl, Real & Epp, Real & Eqq) const override;
virtual void setStressAfterReturn(const RankTwoTensor & stress_trial,
Real p_ok,
Real q_ok,
Real gaE,
const std::vector<Real> & intnl,
const yieldAndFlow & smoothed_q,
const RankFourTensor & Eijkl,
RankTwoTensor & stress) const override;
virtual void consistentTangentOperator(const RankTwoTensor & stress_trial,
Real p_trial,
Real q_trial,
const RankTwoTensor & stress,
Real p,
Real q,
Real gaE,
const yieldAndFlow & smoothed_q,
const RankFourTensor & Eijkl,
bool compute_full_tangent_operator,
RankFourTensor & cto) const override;
virtual RankTwoTensor dpdstress(const RankTwoTensor & stress) const override;
virtual RankTwoTensor dqdstress(const RankTwoTensor & stress) const override;
virtual RankFourTensor d2qdstress2(const RankTwoTensor & stress) const override;
};
(modules/solid_mechanics/include/materials/CappedWeakPlaneCosseratStressUpdate.h)
// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html
#pragma once
#include "CappedWeakPlaneStressUpdate.h"
/**
* CappedWeakPlaneCosseratStressUpdate performs the return-map
* algorithm and associated stress updates for plastic
* models that describe capped weak-plane Cosserat plasticity
*
* It assumes various things about the elasticity tensor, viz
* E(i,i,j,k) = 0 except if k=j
* E(0,0,i,j) = E(1,1,i,j)
*/
class CappedWeakPlaneCosseratStressUpdate : public CappedWeakPlaneStressUpdate
{
public:
static InputParameters validParams();
CappedWeakPlaneCosseratStressUpdate(const InputParameters & parameters);
/**
* Does the model require the elasticity tensor to be isotropic?
*/
bool requiresIsotropicTensor() override { return false; }
protected:
virtual void consistentTangentOperator(const RankTwoTensor & stress_trial,
Real p_trial,
Real q_trial,
const RankTwoTensor & stress,
Real p,
Real q,
Real gaE,
const yieldAndFlow & smoothed_q,
const RankFourTensor & Eijkl,
bool compute_full_tangent_operator,
RankFourTensor & cto) const override;
virtual void setStressAfterReturn(const RankTwoTensor & stress_trial,
Real p_ok,
Real q_ok,
Real gaE,
const std::vector<Real> & intnl,
const yieldAndFlow & smoothed_q,
const RankFourTensor & Eijkl,
RankTwoTensor & stress) const override;
virtual RankTwoTensor dqdstress(const RankTwoTensor & stress) const override;
virtual RankFourTensor d2qdstress2(const RankTwoTensor & stress) const override;
};