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 zz, and the tangential directions by xx and yy. It is convenient to introduce two new stress variables in terms of the stress tensor σ\sigma:

p=σzz   and   q=σxz2+σyz2 .p = \sigma_{zz} \ \ \ \text{and}\ \ \ q = \sqrt{\sigma_{xz}^{2} + \sigma_{yz}^{2}} \ . (1)

In standard elasticity, the stress tensor is symmetric, so an equivalent definition of qq is q=12(σxz+σzx)2+12(σyz+σzy)2q=\sqrt{\frac{1}{2}(\sigma_{xz}+\sigma_{zx})^{2} + \frac{1}{2}(\sigma_{yz}+\sigma_{zy})^{2}}, 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 CC, and friction angle ϕ\phi:

f0=q+ptanϕC .f_{0} = q + p\tan\phi - C \ .

The parameter CC and ϕ\phi 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:

f1=pST ,f_{1} = p - S_{T} \ ,

where STS_{T} 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:

f2=pSC ,f_{2} = - p - S_{C} \ ,

where SCS_{C} is the compressive strength (a positive quantity), which may be constant or harden or soften.

The yield functions f1f_{1} and f2f_{2} place caps'' on the shear yield function f0f_{0}. The combined yield function is simply

f=max(f0,f1,f2) ,f = \max(f_{0}, f_{1}, f_{2}) \ , (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 ψ\psi, which may be constant, or harden or soften. The shear flow potential is

g0=q+ptanψ.g_{0} = q + p\tan\psi .

The tensile flow potential is

g1=p ,g_{1} = p \ ,

and the compressive flow potential is

g2=p .g_{2} = - p \ .

The overall flow potential is

g={g0 if f=f0g1 if f=f1g2 if f=f2 .g = \left\{ \begin{array}{ll} g_{0} & \text{ if } f = f_{0} \\ g_{1} & \text{ if } f = f_{1} \\ g_{2} & \text{ if } f = f_{2}\ . \end{array} \right.(3)

Obviously there are problems here where gg is not defined properly at the corners where f0=f1f_{0}=f_{1} and f0=f2f_{0}=f_{2} (or even f0=f1=f2f_{0}=f_{1}=f_{2}). This is resolved by using smoothing (more on this later).

This plasticity model contains two internal parameters, denote by i0i_{0} and i1i_{1}. It is assumed that

C=C(i0) ,ϕ=ϕ(i0) ,ψ=ψ(i0) ,ST=ST(i1) ,SC=SC(i1) .\begin{split} C & = C(i_{0}) \ , \\ \phi & = \phi(i_{0}) \ , \\ \psi & = \psi(i_{0}) \ , \\ S_{T} & = S_{T}(i_{1}) \ , \\ S_{C} & = S_{C}(i_{1}) \ . \\ \end{split}

That is, i0i_{0} can be thought of as the shear'' internal parameter, while i1i_{1} is the tensile'' internal parameter.

To complete the definition of this plasticity model, the increments of i0i_{0} and i1i_{1} during the return-map process must be defined. The return-map process involves being provided with a trial stress σtrial\sigma^{\mathrm{trial}} and an existing value of the internal parameters ioldi^{\mathrm{old}}, and finding a returned'' stress, σ\sigma, and internal parameters, ii, that satisfy

0=f(σ,i) .σ=σtrialEγgσ ,\begin{split} 0 & = f(\sigma, i) \ . \\ \sigma & = \sigma^{\mathrm{trial}} - E\gamma \frac{\partial g}{\partial\sigma} \ , \\ \end{split}(4)

where EE is the elasticity tensor, and γ\gamma 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 EE and g/σ\partial g/\partial\sigma must be used to define the normal direction'').

Let us express the normality condition in (p,q)(p, q) space. The zzzz component is easy:

p=σzz=σzztrialEzzijγgσij=σzztrialEzzzzγgp ,p = \sigma_{zz} = \sigma_{zz}^{\mathrm{trial}} - E_{zzij}\gamma \frac{\partial g}{\partial \sigma_{ij}} = \sigma_{zz}^{\mathrm{trial}} - E_{zzzz}\gamma \frac{\partial g}{\partial p} \ , (5)

where the last equality holds by assumption (see full list of assumptions below). The xzxz and yzyz components are similar:

σxz=σxztrialExzxzγgqqσxz .\sigma_{xz} = \sigma_{xz}^{\mathrm{trial}} - E_{xzxz}\gamma \frac{\partial g}{\partial q}\frac{\partial q}{\partial \sigma_{xz}} \ . (6)

Another assumption has been made about EE. The final term is

qσxz=σxzq .\frac{\partial q}{\partial\sigma_{xz}} = \frac{\sigma_{xz}}{q} \ .

This means that Eq. (6) can be re-written

σxz2(1+Exzxzγgq1q)2=(σxztrial)2 .\sigma_{xz}^{2} \left( 1 + E_{xzxz}\gamma \frac{\partial g}{\partial q} \frac{1}{q} \right)^{2} = \left(\sigma^{\mathrm{trial}}_{xz}\right)^{2} \ .

A similar equation holds for the yzyz component, and these can be summed and rearranged to yield

q=qtrialExzxzγgq .q = q^{\mathrm{trial}} - E_{xzxz}\gamma \frac{\partial g}{\partial q} \ . (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 pp, qq and γ\gamma.

Consider the case of returning to the shear yield surface. Since g/p=tanψ\partial g/\partial p = \tan\psi and g/q=1\partial g/\partial q = 1 for this flow, the return-map process must solve the following system of equations

0=q+ptanϕC ,p=ptrialEzzzzγtanψ ,q=qtrialExzxzγ .\begin{split} 0 & = q + p\tan\phi - C \ , \\ p & = p^{\mathrm{trial}} - E_{zzzz}\gamma\tan\psi \ , \\ q & = q^{\mathrm{trial}} - E_{xzxz}\gamma \ . \\ \end{split}

The solution satisfied ptrialp=Ezzzzγtanψp^{\mathrm{trial}} - p = E_{zzzz}\gamma\tan\psi and qtrialq=Exzxzγq^{\mathrm{trial}} - q = E_{xzxz}\gamma.

Now consider the case of returning to the tensile yield surface. The equations are

0=pST ,p=ptrialEzzzzγ ,q=qtrial .\begin{split} 0 & = p - S_{T} \ , \\ p & = p^{\mathrm{trial}} - E_{zzzz}\gamma \ , \\ q & = q^{\mathrm{trial}} \ . \\ \end{split}

Comparing these two types of return, it is obvious that qtrialqq^{\mathrm{trial}} - q quantifies the amount of shear failure. Therefore, the following definitions are used in this plasticity model

i0=i0old+qtrialqExzxz ,i1=i1old+ptrialpEzzzz(qtrialq)tanψExzxz .\begin{split} i_{0} & = i_{0}^{\mathrm{old}} + \frac{q^{\mathrm{trial}} - q}{E_{xzxz}} \ , \\ i_{1} & = i_{1}^{\mathrm{old}} + \frac{p^{\mathrm{trial}} - p}{E_{zzzz}} - \frac{(q^{\mathrm{trial}} - q)\tan\psi}{E_{xzxz}} \ . \\ \end{split}

The final term ensures that i1i_{1} does not increase during pure shear failure. The scaling by EE 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 σijold\sigma_{ij}^{\mathrm{old}}, and a total strain increment δϵ\delta \epsilon (that comes from the nonlinear solver proposing displacements) the trial stress is

σijtrial=σijold+Eijklδϵkl ,\sigma_{ij}^{\mathrm{trial}} = \sigma_{ij}^{\mathrm{old}} + E_{ijkl}\delta\epsilon_{kl} \ ,

This gives ptrialp^{\mathrm{trial}} and qtrialq^{\mathrm{trial}}. If ptrialp^{\mathrm{trial}}, qtrialq^{\mathrm{trial}}, i0oldi_{0}^{\mathrm{old}} and i1oldi_{1}^{\mathrm{old}} are such that f(ptrial,qtrial,iold)>0f(p^{\mathrm{trial}}, q^{\mathrm{trial}}, i^{\mathrm{old}}) > 0, the return-map problem is: find pp, qq, γ\gamma, i0i_{0} and i1i_{1} such that

0=f(p,q,i) ,p=ptrialEzzzzγgp ,q=qtrialExzxzγgq ,i0=i0old+qtrialqExzxz ,i1=i1old+ptrialpEzzzz(qtrialq)tanψExzxz .\begin{split} 0 & = f(p, q, i) \ , \\ p & = p^{\mathrm{trial}} - E_{zzzz}\gamma \frac{\partial g}{\partial p} \ , \\ q & = q^{\mathrm{trial}} - E_{xzxz}\gamma \frac{\partial g}{\partial q} \ , \\ i_{0} & = i_{0}^{\mathrm{old}} + \frac{q^{\mathrm{trial}} - q}{E_{xzxz}} \ , \\ i_{1} & = i_{1}^{\mathrm{old}} + \frac{p^{\mathrm{trial}} - p}{E_{zzzz}} - \frac{(q^{\mathrm{trial}} - q)\tan\psi}{E_{xzxz}} \ . \\ \end{split}(8)

The latter two equations are assumed to hold in the smoothed situation (discussed below) too, and note that ψ=ψ(i0)\psi = \psi(i_{0}), so the these two equations are not completely trivial.

After the return-map problem has been solved, the stress components are σij=σijtrial\sigma_{ij} = \sigma_{ij}^{\mathrm{trial}}, except for the following

σxx=σxxtrialEzzxxγgp ,σyy=σyytrialEzzyyγgp ,σzz=p ,σzx=σzxtrialq/qtrial ,σxz=σxztrialq/qtrial ,σzy=σzytrialq/qtrial ,σyz=σyztrialq/qtrial .\begin{split} \sigma_{xx} & = \sigma_{xx}^{\mathrm{trial}} - E_{zzxx}\gamma\frac{\partial g}{\partial p} \ , \\ \sigma_{yy} & = \sigma_{yy}^{\mathrm{trial}} - E_{zzyy}\gamma\frac{\partial g}{\partial p} \ , \\ \sigma_{zz} & = p \ , \\ \sigma_{zx} & = \sigma_{zx}^{\mathrm{trial}} q / q^{\mathrm{trial}} \ , \\ \sigma_{xz} & = \sigma_{xz}^{\mathrm{trial}} q / q^{\mathrm{trial}} \ , \\ \sigma_{zy} & = \sigma_{zy}^{\mathrm{trial}} q / q^{\mathrm{trial}} \ , \\ \sigma_{yz} & = \sigma_{yz}^{\mathrm{trial}} q / q^{\mathrm{trial}} \ . \\ \end{split}

The plastic strain is

ϵijplastic=ϵijplastic,old+δϵij+Eijkl1(σkloldσkl)=ϵijplastic,old+Eijklγgσkl .\epsilon_{ij}^{\mathrm{plastic}} = \epsilon_{ij}^{\mathrm{plastic, old}} + \delta\epsilon_{ij} + E_{ijkl}^{-1}(\sigma_{kl}^{\mathrm{old}} - \sigma_{kl}) = \epsilon_{ij}^{\mathrm{plastic, old}} + E_{ijkl}\gamma \frac{\partial g}{\partial \sigma_{kl}}\ .

The elastic strain is

ϵijelastic=ϵijelastic,old+δϵijϵijplastic+ϵijplastic,old .\epsilon_{ij}^{\mathrm{elastic}} = \epsilon_{ij}^{\mathrm{elastic, old}} + \delta\epsilon_{ij} - \epsilon_{ij}^{\mathrm{plastic}} + \epsilon_{ij}^{\mathrm{plastic, old}} \ .

Yield Smoothing

The shear yield function, f0f_{0}, describes a cone in (σyz,σxz,σzz)(\sigma_{yz}, \sigma_{xz}, \sigma_{zz}) 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

f0=q2+st2+ptanϕC ,f_{0} = \sqrt{q^{2} + s_{t}^{2}} + p\tan\phi - C \ ,

and the flow potential is

g0=q2+st2+ptanψ .g_{0} = \sqrt{q^{2} + s_{t}^{2}} + p\tan\psi \ .

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 ss can be used. The parameter ss has the units of stress. At any point (p,q,i)(p, q, i) order the 3 yield function values, and denote the largest by AA, the second largest by BB and the smallest by CC:

ABCA\geq B\geq C

Then the single, smoothed yield function is defined to be

f={A   if  AB+sA+B+s2sπcos((BA)π2s) .f = \left\{ \begin{array}{ll} A & \ \ \ \text{if}\ \ A\geq B+s \\ \frac{A+B+s}{2} - \frac{s}{\pi}\cos\left(\frac{(B-A)\pi}{2s}\right) \ . \end{array} \right.

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 ψϕ .\psi \leq \phi \ . These conditions should be satisfied for all values of the internal parameter i0i_{0}. MOOSE checks that these conditions hold for i0=0i_{0}=0 only.

The tensile and compressive strength must satisfy STSC ,S_{T} \geq -S_{C} \ , otherwise the caps'' are swapped and the assumption of a convex yield surface is violated. MOOSE checks this condition holds for i1=0i_{1}=0 only: the MOOSE user must ensure that it actually holds for all values of the internal parameter

The smoothing parameter ss 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 ST=1=SCS_{T}=1=S_{C}, then a smoothing parameter of 0.1 is fine, but a smoothing parameter 2\geq 2 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 sts_{t} 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 ss and sts_{t} are chosen so that the yield surface is not wildly varying around q=0q=0, otherwise poor convergence of the return-map process will occur.

It is assumed that the elasticity tensor has the following symmetries:

Eijkl=Ejikl=Eijlk=Eklij ,E_{ijkl} = E_{jikl} = E_{ijlk} = E_{klij} \ , (9)

and that

0=Ezzij   if   ij ,0 = E_{zzij} \ \ \ \text{if}\ \ \ i\neq j \ ,

and that

Exzxz=Eyzyz ,E_{xzxz} = E_{yzyz} \ ,

and that

0=Exzij   unless (i,j)=(z,x)   or (i,j)=(x,z) .0 = E_{xzij} \ \ \text{ unless } (i, j) = (z, x) \ \ \text{ or } (i, j) = (x, z) \ . (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 3×33\times 3 system consisting of the first 3 equations, and substituting the fourth and fifth equations wherever needed. The three unknowns are pp, qq and γE=γEzzzz\gamma_{E}=\gamma E_{zzzz}, 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, δϵ\delta\epsilon, 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

Hijkl=δσijδϵkl .H_{ijkl} = \frac{\delta\sigma_{ij}}{\delta \epsilon_{kl}} \ .

The quantity HH is called the consistent tangent operator. For pure elasticity it is simply the elastic tensor, EE, but it is more complicated for plasticity. Note that a small δϵkl\delta\epsilon_{kl} simply changes δσtrial\delta\sigma^{\mathrm{trial}}, so HH is capturing the change of the returned stress (δσ\delta\sigma) with respect to a change in the trial stress (δσtrial\delta\sigma^{\mathrm{trial}}). In (p,q)(p,q) language, we need to the sx derivatives

δ(p,q,γ)δ(ptrial,qtrial) .\frac{\delta (p, q, \gamma)}{\delta (p^{\mathrm{trial}}, q^{\mathrm{trial}})} \ .

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 (p,q)(p, q) model, while the other is specialized to the weak-plane case.

General consistent tangent operator

The return-map algorithm provides

σij=σijtrialEijmnγgσmn .\sigma_{ij} = \sigma_{ij}^{\mathrm{trial}} - E_{ijmn}\gamma \frac{\partial g}{\partial \sigma_{mn}} \ .

Since σtrial=Eϵ\sigma^{\mathrm{trial}} = E\epsilon, the consistent tangent operator is

Hijkl=EijklEijmnEpqklσpqtrialγgσmn=EijklEijmnEpqkl(ptrialσpqtrialptrial+qtrialσpqtrialqtrial)γ(gppσmn+gqqσmn)\begin{split} H_{ijkl} &= E_{ijkl} - E_{ijmn} E_{pqkl} \frac{\partial}{\partial \sigma_{pq}^{\mathrm{trial}}} \gamma \frac{\partial g}{\partial \sigma_{mn}} \\ &= E_{ijkl} - E_{ijmn} E_{pqkl} \left( \frac{\partial p^{\mathrm{trial}}}{\partial \sigma_{pq}^{\mathrm{trial}}} \frac{\partial}{\partial p^{\mathrm{trial}}} + \frac{\partial q^{\mathrm{trial}}}{\partial \sigma_{pq}^{\mathrm{trial}}} \frac{\partial}{\partial q^{\mathrm{trial}}} \right) \gamma \left( \frac{\partial g}{\partial p} \frac{\partial p}{\partial \sigma_{mn}} + \frac{\partial g}{\partial q} \frac{\partial q}{\partial \sigma_{mn}} \right) \end{split}

However, note that

ptrial(pptrial+γEppgp)=0 ,\frac{\partial}{\partial p^{\mathrm{trial}}} \left( p - p^{\mathrm{trial}} + \gamma E_{pp}\frac{\partial g}{\partial p} \right) = 0 \ ,

because the return-map algorithm guarantees that the expression inside parentheses is zero. Therefore

ptrialγgp=1Epp(1pptrial) .\frac{\partial}{\partial p^{\mathrm{trial}}}\gamma\frac{\partial g}{\partial p} = \frac{1}{E_{pp}} \left( 1- \frac{\partial p}{\partial p^{\mathrm{trial}}} \right) \ .

A similar expression holds for three other cases. There are still terms that involve derivatives of p/σmn\partial p/\partial \sigma_{mn} and q/σmn\partial q/\partial\sigma_{mn}, but these may be separated off as seen below.

The consistent tangent operator may therefore be written as

Hijkl=EijklEijmnEpqkl{ptrialσpqtrial1Epp(1pptrial)pσmnqtrialσpqtrial1Epp(pqtrial)pσmn+ptrialσpqtrial 1Eqq(qptrial)qσmn+qtrialσpqtrial1Eqq(1qqtrial)qσmn}σabϵklEijmnγ(gp2pσmnσab+gq2qσmnσab) .\begin{split} H_{ijkl} & = E_{ijkl} - E_{ijmn}E_{pqkl} \left\{ \frac{\partial p^{\mathrm{trial}}}{\partial\sigma_{pq}^{\mathrm{trial}}} \frac{1}{E_{pp}}\left(1 - \frac{\partial p}{\partial p^{\mathrm{trial}}} \right) \frac{\partial p}{\partial\sigma_{mn}} \right. \\ & \left. \frac{\partial q^{\mathrm{trial}}}{\partial\sigma_{pq}^{\mathrm{trial}}} \frac{1}{E_{pp}}\left( - \frac{\partial p}{\partial q^{\mathrm{trial}}} \right) \frac{\partial p}{\partial\sigma_{mn}} + \frac{\partial p^{\mathrm{trial}}}{\partial\sigma_{pq}^{\mathrm{trial}}}\ \frac{1}{E_{qq}}\left(- \frac{\partial q}{\partial p^{\mathrm{trial}}} \right) \frac{\partial q}{\partial\sigma_{mn}} + \frac{\partial q^{\mathrm{trial}}}{\partial\sigma_{pq}^{\mathrm{trial}}} \frac{1}{E_{qq}}\left(1 - \frac{\partial q}{\partial q^{\mathrm{trial}}} \right) \frac{\partial q}{\partial\sigma_{mn}} \right\} \\ & \frac{\partial \sigma_{ab}}{\partial\epsilon_{kl}} E_{ijmn} \gamma \left( \frac{\partial g}{\partial p}\frac{\partial^{2} p}{\partial\sigma_{mn}\partial\sigma_{ab}} + \frac{\partial g}{\partial q}\frac{\partial^{2} q}{\partial\sigma_{mn}\partial\sigma_{ab}} \right) \ . \\ \end{split}

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 Hijkl=σij/ϵklH_{ijkl} = \partial\sigma_{ij}/\partial\epsilon_{kl}) and the resulting expression multiplied inverse HH's coefficient to finally yield HH. 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 (p,q)(p, q) given the trial variables. Finding HH is really just re-solving these equations for a slightly changed trial variable. Denote

(R0R1R2)=(fp+ptrialEzzzzγgpq+qtrialExzxzγgq) .\left( \begin{array}{l} R_{0} \\ R_{1} \\ R_{2} \end{array} \right) = \left( \begin{array}{l} -f \\ -p + p^{\mathrm{trial}} - E_{zzzz}\gamma \frac{\partial g}{\partial p} \\ -q + q^{\mathrm{trial}} - E_{xzxz}\gamma \frac{\partial g}{\partial q} \end{array} \right) \ .

Then

R0ptrial=fi1iiptrial ,R1ptrial=1Ezzzzγ2gpiiiiptrial ,R2ptrial=Exzxzγ2gqiiiiptrial\begin{split} \frac{\partial R_{0}}{\partial p^{\mathrm{trial}}} & = -\frac{\partial f}{\partial i_{1}} \frac{\partial i_{i}}{\partial p^{\mathrm{trial}}} \ , \\ \frac{\partial R_{1}}{\partial p^{\mathrm{trial}}} & = 1 - E_{zzzz}\gamma \frac{\partial^{2}g}{\partial p\partial i_{i}}\frac{\partial i_{i}}{\partial p^{\mathrm{trial}}} \ , \\ \frac{\partial R_{2}}{\partial p^{\mathrm{trial}}} & = -E_{xzxz}\gamma \frac{\partial^{2}g}{\partial q\partial i_{i}}\frac{\partial i_{i}}{\partial p^{\mathrm{trial}}} \\ \end{split}

In these equations

i1ptrial=1Ezzzz ,\frac{\partial i_{1}}{\partial p^{\mathrm{trial}}} = \frac{1}{E_{zzzz}} \ ,

which comes from Eq. (8). The derivatives with respect to qtrialq^{\mathrm{trial}} are similar but more lengthy due to both i0i_{0} and i1i_{1} being dependent on qtrialq^{\mathrm{trial}}. The system to solve is

(R0γR0pR0qR1γR1pR1qR2γR2pR2q)(δγ/δptrialδp/δptrialδq/δptrial)=(R0/δptrialR1/δptrialR2/δptrial)\left( \begin{array}{ccc} \frac{\partial R_{0}}{\partial \gamma} & \frac{\partial R_{0}}{\partial p} & \frac{\partial R_{0}}{\partial q} \\ \frac{\partial R_{1}}{\partial \gamma} & \frac{\partial R_{1}}{\partial p} & \frac{\partial R_{1}}{\partial q} \\ \frac{\partial R_{2}}{\partial \gamma} & \frac{\partial R_{2}}{\partial p} & \frac{\partial R_{2}}{\partial q} \end{array} \right) \left( \begin{array}{c} \delta \gamma/\delta p^{\mathrm{trial}} \\ \delta p/\delta p^{\mathrm{trial}} \\ \delta q/\delta p^{\mathrm{trial}} \end{array} \right) = \left( \begin{array}{c} \partial R_{0}/\delta p^{\mathrm{trial}} \\ \partial R_{1}/\delta p^{\mathrm{trial}} \\ \partial R_{2}/\delta p^{\mathrm{trial}} \end{array} \right)

The 3×33\times 3 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 qtrialq^{\mathrm{trial}}.

Once the six derivatives have been computed they need to be assembled into HH. For instance,

δptrialδϵii=Ezzii ,\frac{\delta p^{\mathrm{trial}}}{\delta \epsilon_{ii}} = E_{zzii} \ ,

so that

Hzzii=δpδptrialEzzii ,H_{zzii} = \frac{\delta p}{\delta p^{\mathrm{trial}}} E_{zzii} \ ,

and other more complicated expressions appear for other components, such as

Hxxii=ExxiiEzzxxEzzii(δγδptrialgp+γ2gp2δpδptrial+γ2gpi1δqδptrial+γ2gpi0i0qδqδptrial    +γ2gpi1(δi1δptrial+i1pδpδptriali1qδqδptrial)) .\begin{split} H_{xxii} & = E_{xxii} - E_{zzxx} E_{zzii}\left( \frac{\delta \gamma}{\delta p^{\mathrm{trial}}} \frac{\partial g}{\partial p} + \gamma \frac{\partial^{2} g}{\partial p^{2}} \frac{\delta p}{\delta p^{\mathrm{trial}}} + \gamma \frac{\partial^{2} g}{\partial p\partial i_{1}} \frac{\delta q}{\delta p^{\mathrm{trial}}} + \gamma \frac{\partial^{2} g}{\partial p\partial i_{0}} \frac{\partial i_{0}}{\partial q} \frac{\delta q}{\delta p^{\mathrm{trial}}} \right. \\ & \ \ \ \ + \left. \gamma \frac{\partial^{2} g}{\partial p\partial i_{1}} \left( \frac{\delta i_{1}}{\delta p^{\mathrm{trial}}} + \frac{\partial i_{1}}{\partial p} \frac{\delta p}{\delta p^{\mathrm{trial}}} \frac{\partial i_{1}}{\partial q} \frac{\delta q}{\delta p^{\mathrm{trial}}} \right) \right) \ . \\ \end{split}

The consistent tangent operator and sub-stepping strain increments

One extra complication arises from the potential sub-stepping of the applied strain increment δϵ\delta\epsilon. 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 3×33\times 3 system which is pretty quick for the computer to compared with the entire Newton-Raphson process.

Let the nthn^{\mathrm{th}} substep be

δϵn=λnδϵ ,\delta\epsilon^{n} = \lambda_{n}\delta\epsilon \ ,

with

1=n=1Nλn ,1 = \sum_{n=1}^{N}\lambda_{n} \ ,

where NN is the total number of substeps. Denoting the initial stress by (pold,qold)(p^{\mathrm{old}}, q^{\mathrm{old}}), and the returned stress at step n1n-1 by (pn1,qn1)(p_{n-1}, q_{n-1}), of course the trial stress at step nn is

(pntrial,qntrial)=(pn1,qn1)+λn(ptrialpold,qtrialqold) .(p_{n}^{\mathrm{trial}}, q_{n}^{\mathrm{trial}}) = (p_{n-1},q_{n-1}) + \lambda_{n} ( p^{\mathrm{trial}} - p^{\mathrm{old}}, q^{\mathrm{trial}} - q^{\mathrm{old}}) \ .

This means that

pnptrial=pnpntrialpntrialptrial+pnqntrialqntrialptrial=pnpntrial(λn+pn1ptrial)+pnqntrialqn1ptrial .\begin{split} \frac{\partial p_{n}}{\partial p^{\mathrm{trial}}} & = \frac{\partial p_{n}}{\partial p_{n}^{\mathrm{trial}}} \frac{\partial p_{n}^{\mathrm{trial}}}{\partial p^{\mathrm{trial}}} + \frac{\partial p_{n}}{\partial q_{n}^{\mathrm{trial}}} \frac{\partial q_{n}^{\mathrm{trial}}}{\partial p^{\mathrm{trial}}} \\ & = \frac{\partial p_{n}}{\partial p_{n}^{\mathrm{trial}}} \left(\lambda_{n} + \frac{\partial p_{n-1}}{\partial p^{\mathrm{trial}}} \right) + \frac{\partial p_{n}}{\partial q_{n}^{\mathrm{trial}}} \frac{\partial q_{n-1}}{\partial p^{\mathrm{trial}}} \ . \\ \end{split}

Similar inductive equations hold for the other derivatives, and note that p0/ptrial=pold/ptrial=0\partial p_{0}/\partial p^{\mathrm{trial}} = \partial p^{\mathrm{old}}/\partial p^{\mathrm{trial}} = 0. The derivative of γ\gamma is slightly different: it is

γptrial=n=1Nγnpntrial(λn+pn1ptrial)+γnqntrialqn1ptrial ,\frac{\partial \gamma}{\partial p^{\mathrm{trial}}} = \sum_{n=1}^{N} \frac{\partial \gamma_{n}}{\partial p_{n}^{\mathrm{trial}}} \left(\lambda_{n} + \frac{\partial p_{n-1}}{\partial p^{\mathrm{trial}}} \right) + \frac{\partial \gamma_{n}}{\partial q_{n}^{\mathrm{trial}}} \frac{\partial q_{n-1}}{\partial p^{\mathrm{trial}}} \ ,

and similarly for the derivative with respect to qtrialq^{\mathrm{trial}}.

Input Parameters

  • 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.

Required 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