- 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)
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)
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