MultiParameterPlasticityStressUpdate

MultiParameterPlasticityStressUpdate is a base class to the following stress update materials:

Yield functions smoothing

The three yield functions shown in Tensile Stress Update are smoothed in the MultiParameterPlasticityStressUpdate class.

An example is shown Figure 1 and Figure 2. Figure 2 shows slices of the yield surface at three values of the mean stress, and the triangular-pyramid nature of the yield surface is evident. The slices are taken near the tip region only, to highlight: (1) the smoothing; (2) that the smoothing is unsymmetric.

The unsymmetric nature of the yield surface only occurs near the tip region where the smoothing mixes the three yield surfaces. For instance, the black line in Figure 2 is symmetric, while the blue and red lines are unsymmetric. The amount of asymmetry is small, but it is evident that the red curve is not concentric with the remainder of the curves shown in Figure 2. The order of the yield functions has been chosen so that the curves intersect the σIII=σII\sigma_{III}=\sigma_{II} line at 9090^{\circ}, but not on the σII=σI\sigma_{II}=\sigma_{I} line near the pyramid's tip. The asymmetry does not affect MOOSE's convergence, and of course it is physically irrelevant (since there is no one correct'' smoothed yield surface).

The unsmoothed yield surface defined by the yield function equations is a triangular pyramid

Figure 1: A smoothed version of the yield surface.

The principal stress directions are shown with black lines, and the mean stress direction is shown with a blue line. The three planes extend infinitely (there is no base to the triangular pyramid) but these pictures show the tip region only.

Figure 2: Slices of the unsmoothed and smoothed yield functions at various value of mean stress.

In this example T=1T=1, the smoothing tolerance is 0.5, and the values of mean stress shown are: black, 0.7; blue, 0.8, red 0.85. An unusually large amount of smoothing has been used here to highlight various features: in real simulations users are encouraged to use smoothing approximately T/10T/10. The whole octahedral plane is shown in this figure, but only one sextant is physical, which is indicated by the solid green lines.

Flow rules and hardening

This plasticity is associative: the flow potential is

g=f=max(f0,f1,f2)+smoothing . g = f = max(f_{0}, f_{1}, f_{2}) + \text{smoothing} \ .

Here smoothing indicates the smoothing mentioned in the previous section.

The flow rules are

sa=satrialγEabgsa , s_{a} = s_{a}^{\mathrm{trial}} - \gamma E_{ab} \frac{\partial g}{\partial s_{a}} \ , (1)

where sa={σI,σII,σIII}s_{a}=\{\sigma_{I}, \sigma_{II}, \sigma_{III}\} and

Eab=saσijEijklsbσkl . E_{ab} = \frac{\partial s_{a}}{\partial \sigma_{ij}} E_{ijkl} \frac{\partial s_{b}}{\partial \sigma_{kl}} \ .

In this equation EijklE_{ijkl} is the elasticity tensor.

An assumption that is made is that EabE_{ab} is independent of the stress parameters, sas_{a} and the internal variables. In this case,

saσij=viavja , \frac{\partial s_{a}}{\partial \sigma_{ij}} = v_{i}^{a}v_{j}^{a} \ ,

commentnote

Special precautions are taken when the eigenvalues are equal, as described in the documentation for RankTwoTensor.

where vav^{a} is the eigenvector corresponding to the eigenvalue sas^{a} (of the stress tensor) and there is no sum over aa on the right-hand side. Recall that the eigenvectors are fixed during the return-map process, so the RHS is fixed, meaning that EabE_{ab} is indeed independent of the stress parameters. Also recall that the eigenvectors induce a rotation (to the principal-stress frame), so assuming that EijklE_{ijkl} is isotropic

Eab=Eaabb . E_{ab} = E_{aabb} \ .

The assumption of isotropy is appropriate for this type of isotropic plasticity.

It is assumed that there is just one internal parameter, ii, and that it is defined by

i=iold+(σIIItrialσIII)/E22 , i = i_{\mathrm{old}} + (\sigma_{III}^{\mathrm{trial}} - \sigma_{III}) / E_{22} \ ,

during the return-map process. The tensile strength is a function of this single internal parameter

Return-map considerations

Unknowns and the convergence criterion

The return-map problem involves solving the four equations: f=0f=0 (smoothed yield function should be zero) and the flow Eq. (1). The unknowns are the 3 stress parameters sa={σI,σII,σIII}s_{a}=\{\sigma_{I}, \sigma_{II}, \sigma_{III}\} and the plasticity multiplier γ\gamma. Actually, to make the units consistent the algorithm uses γE22\gamma E_{22} instead of simply γ\gamma. Convergence is deemed to be achieved when the sum of squares of the residuals of these 4 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 may then be dependent on the number of sub-steps. While this is simply illustrating the non-uniqueness of plasticity problems, it has been observed not to 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 formulae:

Hijkl=δσijδσmntrialδσmntrialδϵkl=δσijδσmntrialEmnkl . H_{ijkl} = \frac{\delta\sigma_{ij}}{\delta \sigma_{mn}^{\mathrm{trial}}} \frac{\delta\sigma_{mn}^{\mathrm{trial}}}{\delta\epsilon_{kl}} =\frac{\delta\sigma_{ij}}{\delta \sigma_{mn}^{\mathrm{trial}}} E_{mnkl} \ .

In the case at hand,

σij=aRiasaRajT . \sigma_{ij} = \sum_{a}R_{ia}s_{a}R_{aj}^{\mathrm{T}} \ .

In this formula σij\sigma_{ij} is the returned stress, sas_{a} are the returned stress parameters (eigenvalues), and RR is the rotation matrix, defined through the eigenvectors, vav^{a} (a=1,2,3a=1,2,3) of the trial stress:

Ria=via . R_{ia} = v^{a}_{i} \ .

The three eigenvectors remain unchanged during the return-map process. However, of course they change under a change in σtrial\sigma^{\mathrm{trial}}. The relevant formulae are

δsatrialδσkltrial=viavja ,δviaδσkltrial=bavib(vkbvla+vlbvka)2(sasb) .\begin{split} \frac{\delta s_{a}^{\mathrm{trial}}}{\delta \sigma_{kl}^{\mathrm{trial}}} & = v^{a}_{i}v^{a}_{j} \ , \\ \frac{\delta v^{a}_{i}}{\delta \sigma_{kl}^{\mathrm{trial}}} & = \sum_{b\neq a}\frac{v_{i}^{b}(v_{k}^{b}v_{l}^{a} + v_{l}^{b}v_{k}^{a})}{2(s_{a}-s_{b})} \ . \\ \end{split}

On the RHS of these equations there is no sum over aa.

The final piece of information is

δsbδsatrial . \frac{\delta s_{b}}{\delta s_{a}^{\mathrm{trial}}} \ .

MultiParameterPlasticityStressUpdate computes this after each Newton step, for any arbitrary plasticity model.

The nontrivial part to the consistent tangent operator is therefore

δσijδσmntrial=aδRiaδσmntrialsaRajT+abRiaδsaδsbtrialδsbtrialδσmntrialRajT+aRiasaδRajTδσmntrial . \frac{\delta \sigma_{ij}}{\delta\sigma_{mn}^{\mathrm{trial}}} = \sum_{a} \frac{\delta R_{ia}}{\delta\sigma_{mn}^{\mathrm{trial}}}s_{a}R_{aj}^{\mathrm{T}} + \sum_{a}\sum_{b} R_{ia}\frac{\delta s_{a}}{\delta s_{b}^{\mathrm{trial}}} \frac{\delta s_{b}^{\mathrm{trial}}}{\delta \sigma_{mn}^{\mathrm{trial}}}R_{aj}^{\mathrm{T}} + \sum_{a} R_{ia}s_{a}\frac{\delta R_{aj}^{\mathrm{T}}}{\delta\sigma_{mn}^{\mathrm{trial}}} \ .

All the components of this equation have been provided above.