This paper is dedicated to the implementation a perfect plastic behaviour based on the Hosford equivalent stress.

The whole implementation is available here.

Description

Elasticity

The elasticity is assumed linear and isotropic, i.e. given by the Hooke law: \[ \underline{\sigma}=\lambda\,{\mathrm{tr}{\left(\underline{\varepsilon}^{\mathrm{el}}\right)}}\,\underline{I}+2\,\mu\,\underline{\varepsilon}^{\mathrm{el}} \] where \(\lambda\) and \(\mu\) are the first and second Lamé parameters.

Plasticity

The yield surface is chosen as follows: \[ f{\left(\underline{\sigma}\right)}=\sigma_{\mathrm{eq}}^{H}-\sigma_{Y} \] where:

The Hosford equivalent stress is defined by (see [1]): \[ \sigma_{\mathrm{eq}}^{H}=\sqrt[a]{{{\displaystyle \frac{\displaystyle 1}{\displaystyle 2}}}{\left({\left|\sigma_{1}-\sigma_{2}\right|}^{a}+{\left|\sigma_{1}-\sigma_{3}\right|}^{a}+{\left|\sigma_{2}-\sigma_{3}\right|}^{a}\right)}} \] where \(\sigma_{1}\), \(\sigma_{2}\) and \(\sigma_{3}\) are the eigenvalues of the stress.

Therefore, when \(a\) goes to infinity, the Hosford stress reduces to the Tresca stress. When \(n = 2\) the Hosford stress reduces to the von Mises stress.

Comparison of the Hosford stress a=100,a=8 and the von Mises stress in plane stress

The Hosford stress is an isotropic homogeneous function of degree 1. The flow rule is assumed associated and the plastic strain rate \(\underline{\dot{\varepsilon}}^{\mathrm{p}}\) is given by: \[ \underline{\dot{\varepsilon}}^{\mathrm{p}}=\dot{p}\,{\displaystyle \frac{\displaystyle \partial \sigma_{\mathrm{eq}}^{H}}{\displaystyle \partial \underline{\sigma}}}=\dot{p}\,\underline{n}^{H} \] where \(\dot{p}\) is the rate of the equivalent plastic strain \(p\).

Implicit scheme

Choice of the state variables

The behaviour is integrated by an implicit scheme. Two state variables are introduced:

The elastic strain is automatically defined by the Implicit domain specific language.

The latter could be considered as an integration variable, but, for post-processing purposes, we choose to keep it as a state variable.

Elastic prediction

First, an elastic prediction of the stress \(\underline{\sigma}^{\mathrm{tr}}\) is made (The following expression is not valid in plane stress hypothesis, see below): \[ \underline{\sigma}^{\mathrm{tr}}=\lambda\,{\mathrm{tr}{\left({\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t}}+\theta\,\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)}}\,\underline{I}+2\,\mu\,{\left({\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t}}+\theta\,\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)} \]

Equations governing the material evolution under plastic loading

The equation associated with the evolution of the elastic strain is given by the split of strain: \[ f_{\underline{\varepsilon}^{\mathrm{el}}}=\Delta\,\underline{\varepsilon}^{\mathrm{el}}-\Delta\,\underline{\varepsilon}^{\mathrm{to}}+\Delta\,p\,{\left.\underline{n}^{H}\right|_{t+\theta\,\Delta\,t}} \]

The derivatives of this equation with respect to \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}\) and \(\Delta\,p\) are given by: \[ \left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} &= \underline{\underline{\mathbf{I}}}+2\,\mu\,\theta\,\Delta\,p\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{n}^{H}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\\ {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} &= {\left.\underline{n}^{H}\right|_{t+\theta\,\Delta\,t}}\\ \end{aligned} \right. \]

To determine the equivalent plastic strain increment, the following equation must be satisfied: \[ f_{p}={\left.\sigma_{\mathrm{eq}}^{H}\right|_{t+\theta\,\Delta\,t}}-\sigma_{Y}=0 \]

In the implementation described below, this equation will be normalised by the Young modulus to ensure that all equations have the same magnitude.

The derivatives of this equation with respect to \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}\) and \(\Delta\,p\) are given by: \[ \left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial f_{p}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} &= 2\,\mu\,\theta\,{\left.\underline{n}^{H}\right|_{t+\theta\,\Delta\,t}}\\ {\displaystyle \frac{\displaystyle \partial f_{p}}{\displaystyle \partial \Delta\,p}} &= 0\\ \end{aligned} \right. \]

Implementation

Metadata

The beginning of the file gives some information about the behaviour:

@DSL       Implicit;
@Behaviour HosfordPerfectPlasticity;
@Author    Thomas Helfer;
@Description{
  A simple implementation of a perfect
  plasticity behaviour using the
  Hosford stress.
};

Supported modelling hypothesis

Thanks to the StandardElasticity brick, all the modelling hypotheses can be supported. The following statement, starting with the @ModellingHypotheses, enables all the modelling hypotheses:

@ModellingHypotheses {".+"};

The standard elasticity brick

To implement this behaviour, we will use the StandardElasticity brick which provides:

This behaviour brick is fully described here.

The usage of the StandardElasticity is introduced as follows:

@Brick StandardElasticity;

Numerical parameters

The following part of file give some default values for numerical parameters used by the integration algorithm:

@Epsilon 1.e-16;
@Theta 1;

State variables

The elastic strain is automatically declared the StandardElasticity brick. The associated variable is eel.

The following statement introduces the equivalent plastic strain named p:

@StateVariable strain p;
p.setGlossaryName("EquivalentPlasticStrain");

Material constants

The material properties are hard-coded. The @ElasticMaterialProperties is used to declare the Young modulus and the Poisson ratio.

@ElasticMaterialProperties {150e9,0.3};

In the Implicit scheme, the lame coefficients are automatically deduced from the Young modulus and the Poisson ratio. They are accessible though the lambda and mu local variables which are automatically defined.

The parameters associated with the plastic part of the behaviour are defined as follows:

@Parameter a    = 8;
@Parameter sigy = 150e6;

Here a stands for the Hosford exponent and sigy is the yield stress.

Local variable declaration

A boolean b is declared as a local variable.

@LocalVariable bool b;

This boolean will be true if plastic loading occurs.

Local variable initialization

The main goal of the local variable initialization is to test if the elastic prediction of the stress lies inside the yield surface.

@InitializeLocalVariables{
  const stress seps = 1.e-10*young;
  const auto sigel = computeElasticPrediction();
  const auto seqel = computeHosfordStress(sigel,a,seps);
  b = seqel>sigy;
}

The computeElasticPrediction method, provided by the StandardElasticity brick, computes the elastic prediction of the stress and takes into account the modelling hypothesis. This prediction is thus valid under the plane stress hypothesis.

The computeHosfordStress is then called. It takes three arguments:

At the last line, the boolean variable b is set to true if the elastic prediction of the Hosford stress exceeds the yield stress.

Implicit system

The code describing the implicit system is rather short:

@Integrator{
  const stress seps = 1.e-10*young;
  if(!b){
    // elastic loading, nothing to be done
    return true;
  }
  real seq;
  Stensor n;
  Stensor4 dn;
  std::tie(seq,n,dn) = computeHosfordStressSecondDerivative(sig,a,seps);
  feel        += dp*n;
  dfeel_ddeel += 2*mu*theta*dp*dn;
  dfeel_ddp    = n;
  fp           = (seq-sigy)/young;
  dfp_ddeel    = 2*mu*theta*n/young;
  // this is a small trick that is necessary if the first time
  // step leads to a plastic loading:
  // on this case, the stress are zero on the first iteration and
  // this results in a null line in the jacobian matrix because
  // the normal is then undefined.
  // To avoid the issue, we arbitrarily set dfp_ddp to 1 in that case. 
  dfp_ddp      = (seq<seps) ? 1:  0;
}

It shall be noted that, at the beginning of this code block:

Thus, all the variables have been set to describe a purely elastic behaviour. Hence, nothing is to be done if the boolean variable b is false. In this case, one just return true.

If a plastic loading has been predicted, one uses the computeHosfordStressSecondDerivative function which returns:

The implicit system is then readily written, using expressions given in the previous paragraph.

Note With C++-17 structured bindings enabled, the previous code is much more readable:

@Integrator{
  const stress seps = 1.e-10*young;
  if(!b){
    // elastic loading, nothing to be done
    return true;
  }
  const auto [seq,n,dn] = computeHosfordStressSecondDerivative(sig,a,seps);
  feel        += dp*n;
  dfeel_ddeel += 2*mu*theta*dp*dn;
  dfeel_ddp    = n;
  fp           = (seq-sigy)/young;
  dfp_ddeel    = 2*mu*theta*n/young;
  dfp_ddp      = 0;
}

Determining the yield surface in plane stress

In this section, we show how to use a simple python script to determine the yield surface in plane stress.

Compilation of the behaviour

$ mfront --obuild --interface=abaqus HosfordPerfectPlasticity.mfront
Treating target : all
The following library has been built :
- libABAQUSBEHAVIOUR.so :  HOSFORDPERFECTPLASTICITY_AXIS HOSFORDPERFECTPLASTICITY_PSTRESS HOSFORDPERFECTPLASTICITY_PSTRAIN HOSFORDPERFECTPLASTICITY_3D

A simple python script to determine the yield surface in plane stress

The following script determines the yield stress under the plane stress assumption in the principal stress space by a brute force approach. It imposes strain paths in varying directions and a constraint the shear stress is null. The computations stops when the plastic equivalent strain exceeds \(10^{-3}\). The current stress state is then printed.

from math import pi,cos,sin,sqrt
import mtest

for theta in [pi*(-1+2*float(i)/99) for i in range(0,100)]:
    em   = 2.e-2
    npas = 100
    tmax = 1
    c    = cos(theta)
    s    = sin(theta)
    m    = mtest.MTest()
    mtest.setVerboseMode(mtest.VerboseLevel.VERBOSE_QUIET)
    m.setModellingHypothesis('PlaneStress')
    m.setMaximumNumberOfSubSteps(1)
    m.setBehaviour('abaqus','src/libABAQUSBEHAVIOUR.so',
                   'HOSFORDPERFECTPLASTICITY_PSTRESS');
    m.setExternalStateVariable("Temperature",293.15)
    m.setImposedStrain('EXX',{0:0,tmax:em*c})
    m.setImposedStrain('EYY',{0:0,tmax:em*s})
    m.setNonLinearConstraint('SXY','Stress')
    s  = mtest.MTestCurrentState()
    wk = mtest.MTestWorkSpace()    
    m.completeInitialisation()    
    m.initializeCurrentState(s)
    m.initializeWorkSpace(wk)
    ltime=[(tmax/(npas-1))*i for i in range(npas)]
    for i in range(npas-1):
        m.execute(s,wk,ltime[i],ltime[i+1])
        p = s.getInternalStateVariableValue('EquivalentPlasticStrain');
        if(p>0.001):
            print(str(theta)+" "+str(s.s1[0])+" "+str(s.s1[1]))
            break

Comparing the robustness of various algorithms

This section compares the robustness of various algorithms used to solve the implicit system used to integrate the Hosford behaviour.

Description of the test

The test presented is based on the paper of Scherzinger (see [2]). The test is based on the following steps:

For each values of \(\alpha\) and \(x\), the number of iterations needed to solve the implicit system is recorded.

The implementation presented so far is modified to declare an additional state variable using to save the number of iterations required to reach the convergence.

The python script which implements the test is the following, in the case of a Hosford exponent egal to \(100\):

from math import pi,cos,sin
import tfel.math     as tmath
import tfel.material as tmaterial
import mtest

nmax  = 1000
nmax2 = 1000
E  = 150e9
nu = 0.3
s0 = 150e6

for a in [pi*(-1.+(2.*i)/(nmax-1)) for i in range(0,nmax)]:
    for x in [1+(29.*i)/(nmax2-1)  for i in range(0,nmax2)]:
        nbiter=0
        s     = tmath.makeStensor3D(tmaterial.buildFromPiPlane(cos(a),sin(a)))
        seq   = tmaterial.computeHosfordStress(s,100,1.e-12);
        s    *= x*(s0/seq)
        e0    = (s[0]-nu*(s[1]+s[2]))/E
        e1    = (s[1]-nu*(s[0]+s[2]))/E
        e2    = (s[2]-nu*(s[0]+s[1]))/E
        m    = mtest.MTest()
        mtest.setVerboseMode(mtest.VerboseLevel.VERBOSE_QUIET)
        m.setModellingHypothesis('Tridimensional')
        m.setMaximumNumberOfSubSteps(1)
        m.setBehaviour('abaqus','src/libABAQUSBEHAVIOUR.so',
                       'HOSFORDPERFECTPLASTICITY100_3D');
        m.setExternalStateVariable("Temperature",293.15)
        m.setImposedStrain('EXX',{0:0,1:e0})
        m.setImposedStrain('EYY',{0:0,1:e1})
        m.setImposedStrain('EZZ',{0:0,1:e2})
        m.setImposedStrain('EXY',0)
        m.setImposedStrain('EXZ',0)
        m.setImposedStrain('EYZ',0)
        s  = mtest.MTestCurrentState()
        wk = mtest.MTestWorkSpace()
        m.completeInitialisation()
        m.initializeCurrentState(s)
        m.initializeWorkSpace(wk)
        try:
            m.execute(s,wk,0,1)
            nbiter=s.getInternalStateVariableValue("NumberOfIterations")
        except:
            nbiter=100
            break
        print(str((x/seq)*cos(a))+" "+str((x/seq)*sin(a))+" "+str(nbiter))
    print()

This script outputs, for each couple \({\left(\alpha,x\right)}\) for which the behaviour integration succeeds, the coordinates \(\tilde{\pi}{\left(\underline{\sigma}_{\alpha}\right)}\) scaled by \(\sigma_{Y}\) and the number of iterations required to reach the convergence.

Algorithms tested

The following algorithms are considered:

if(seq>3*(sigy/2)){
  return false;
}

Hosford \(a=6\)

The Hosford yield surface for an exponent egal to \(a=6\) is smooth.

Figure 1: Robustness of Algorithm 1 for a=6

Figure 1 displays the output of the previous script in the \(\pi\)-plane (each coordinate being scaled by \(\sigma_{Y}\)). The maximum number of iteration allowed to reach the convergence is \(100\) by default. However, as stated by Scherzinger, if the convergence is not reached for \(30\) iterations, most of the time the convergence will not be reached.

The results depicted on Figure 1 confirms the ones of Scherzinger: the standard Newton algorithm conditionally converges. Indeed, it can be seen that the standard Newton algorithm converges for all values of \(\alpha\) if \(x\) is low enough. For large values of \(x\), the convergence clearly depends on \(\alpha\).

Figure 2: Robustness of Algorithm 2 for a=6

Figure 2 displays the output of the previous script for the second algorithm. This figure shows that this slight modification of the standard Newton algorithm considerably increases the robustness: the convergence is reached for all tested values of \(\alpha\) and \(x\).

Results for the third and fourth algorithms are not shown but are very similar to the one obtained with the second algorithm. Indeed, in terms of robustness, differences between those three algorithms only arises for high values of the Hosford exponent \(a\).

Hosford \(a=8\)

For an exponent equal to \(a=8\), the curvature of the yield surface varies more rapidly.

Figure 3: Robustness of Algorithm 1 for a=8

As already stated by Scherzinger, Figure 3 shows that the performance in terms of robustness of the standard Newton algorithm is getting worse as the Hosford exponent increases.

Figure 4: Robustness of Algorithm 2 for a=8

Again, Algorithm 2 shows a remarkable benefit in terms of robustness.

Hosford \(a=100\)

For a very high exponent \(a=100\), the yield surface is very close to the Tresca yield surface and shows very sharp edges.

Figure 5: Robustness of Algorithm 1 for a=100

Figure 5 shows that the standard Newton algorithm performs very poorly.

Figure 6: Robustness of Algorithm 2 for a=100

In this case, the second algorithm does not improve the robustness as shown on Figure 6. This algorithm performs almost as poorly than the standard Newton algorithm.

Figure 7: Robustness of Algorithm 3 for a=100

Figure 7 shows that the third algorithm clearly increases the robustness in some directions but the robustness remains low for angles close the edges of the yield surfaces.

Figure 8: Robustness of Algorithm 4 for a=100

Figure 7 shows that the Levenberg-Marquardt algorithm outperforms Algorithm 3. However, the robustness remains low for angles close the edges of the yield surfaces.

Influence of the eigen solver

In this section, we consider the influence of the eigen solver. In the previous tests, the default eigen solver was used. This analytical solver has been shown to provide a very interesting trade-off between accuracy and numerical efficiency. It was also shown that the Jacobi eigensolver provides very accurate results but is less efficient by our default eigen solver (about 2 times less efficient).

Our tests, not reported here, shows that changing the eigen solver:

Figure 9: Robustness of Algorithm 3 for a=100

Figure 9 shows that using the Jacobi eigen solver has a strong influence on the robustness of Algorithm 3.

Figure 10: Robustness of Algorithm 4 for a=100

Using the Jacobi eigen solver with Algorithm 4 shows that convergence is reached for all tested values of \(\alpha\) and \(x\).

Conclusions

Our tests shows that for high values of the Hosford exponent, sophisticated algorithms and accurate eigen solvers are required to guarantee the convergence of the algorithm for large trial stress.

However, for standard values of the Hosford exponent (\(a=6\) and \(a=8\)), our results shows that Algorithm 2, which a very simple modification of the standard Newton algorithm based on a simple and intuitive physical consideration, is very robust.

Those results has to confirmed on real structural analyses. In particular, the impact of each algorithms on the overall efficiency of the simulations has to be evaluated.

References

1.
Hosford, W. F. A generalized isotropic yield criterion. Journal of Applied Mechanics. 1972. Vol. 39, no. 2, p. 607–609.
2.
Scherzinger, W. M. A return mapping algorithm for isotropic and anisotropic plasticity models using a line search method. Computer Methods in Applied Mechanics and Engineering. 15 April 2017. Vol. 317, p. 526–553. DOI 10.1016/j.cma.2016.11.026. Available from: http://www.sciencedirect.com/science/article/pii/S004578251630370X