This paper is dedicated to the implementation a perfect plastic behaviour based on the Hosford equivalent stress.
The whole implementation is available here.
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.
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.
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\).
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.
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)} \]
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. \]
The beginning of the file gives some information about the behaviour:
@DSL
keyword.@Behaviour
keyword.@Author
).@Description
).@DSL Implicit;
@Behaviour HosfordPerfectPlasticity;
@Author Thomas Helfer;
@Description{
A simple implementation of a perfectusing the
plasticity behaviour
Hosford stress. };
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 {".+"};
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;
The following part of file give some default values for numerical parameters used by the integration algorithm:
@Epsilon 1.e-16;
@Theta 1;
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;
"EquivalentPlasticStrain"); p.setGlossaryName(
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.
A boolean b
is declared as a local variable.
@LocalVariable bool b;
This boolean will be true
if plastic loading occurs.
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.
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;2*mu*theta*dp*dn;
dfeel_ddeel +=
dfeel_ddp = n;
fp = (seq-sigy)/young;2*mu*theta*n/young;
dfp_ddeel = // 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.
1: 0;
dfp_ddp = (seq<seps) ? }
It shall be noted that, at the beginning of this code block:
feel
has been initialized to \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}-\Delta\,\underline{\varepsilon}^{\mathrm{to}}\) by the StandardElasticity
brickfp
has been initialized to \(\Delta\,p\) following standard conventions defined in the the Implicit
domain specific language.Implicit
domain specific language.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:
seq
n
with respect to the stressdn
with respect to the stressThe 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;2*mu*theta*dp*dn; dfeel_ddeel += dfeel_ddp = n; fp = (seq-sigy)/young;2*mu*theta*n/young; dfp_ddeel = 0; dfp_ddp = }
In this section, we show how to use a simple python
script to determine the yield surface in plane stress.
$ 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
python
script to determine the yield surface in plane stressThe 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)]:
= 2.e-2
em = 100
npas = 1
tmax = cos(theta)
c = sin(theta)
s = mtest.MTest()
m
mtest.setVerboseMode(mtest.VerboseLevel.VERBOSE_QUIET)'PlaneStress')
m.setModellingHypothesis(1)
m.setMaximumNumberOfSubSteps('abaqus','src/libABAQUSBEHAVIOUR.so',
m.setBehaviour('HOSFORDPERFECTPLASTICITY_PSTRESS');
"Temperature",293.15)
m.setExternalStateVariable('EXX',{0:0,tmax:em*c})
m.setImposedStrain('EYY',{0:0,tmax:em*s})
m.setImposedStrain('SXY','Stress')
m.setNonLinearConstraint(= mtest.MTestCurrentState()
s = mtest.MTestWorkSpace()
wk
m.completeInitialisation()
m.initializeCurrentState(s)
m.initializeWorkSpace(wk)=[(tmax/(npas-1))*i for i in range(npas)]
ltimefor i in range(npas-1):
+1])
m.execute(s,wk,ltime[i],ltime[i= s.getInternalStateVariableValue('EquivalentPlasticStrain');
p if(p>0.001):
print(str(theta)+" "+str(s.s1[0])+" "+str(s.s1[1]))
break
This section compares the robustness of various algorithms used to solve the implicit system used to integrate the Hosford behaviour.
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
= 1000
nmax = 1000
nmax2 = 150e9
E = 0.3
nu = 150e6
s0
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)]:
=0
nbiter= tmath.makeStensor3D(tmaterial.buildFromPiPlane(cos(a),sin(a)))
s = tmaterial.computeHosfordStress(s,100,1.e-12);
seq *= x*(s0/seq)
s = (s[0]-nu*(s[1]+s[2]))/E
e0 = (s[1]-nu*(s[0]+s[2]))/E
e1 = (s[2]-nu*(s[0]+s[1]))/E
e2 = mtest.MTest()
m
mtest.setVerboseMode(mtest.VerboseLevel.VERBOSE_QUIET)'Tridimensional')
m.setModellingHypothesis(1)
m.setMaximumNumberOfSubSteps('abaqus','src/libABAQUSBEHAVIOUR.so',
m.setBehaviour('HOSFORDPERFECTPLASTICITY100_3D');
"Temperature",293.15)
m.setExternalStateVariable('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)
m.setImposedStrain(= mtest.MTestCurrentState()
s = mtest.MTestWorkSpace()
wk
m.completeInitialisation()
m.initializeCurrentState(s)
m.initializeWorkSpace(wk)try:
0,1)
m.execute(s,wk,=s.getInternalStateVariableValue("NumberOfIterations")
nbiterexcept:
=100
nbiterbreak
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.
The following algorithms are considered:
@Integrator
code block: if the test fails, one just have to return the false
value. This is obtained by adding the following line after the computation of the Hosford equivalent stress:if(seq>3*(sigy/2)){
return false;
}
setMaximumIncrementValuePerIteration
method on the eel
variable.@Algorithm
keyword.The Hosford yield surface for an exponent egal to \(a=6\) is smooth.
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 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\).
For an exponent equal to \(a=8\), the curvature of the yield surface varies more rapidly.
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.
Again, Algorithm 2 shows a remarkable benefit in terms of robustness.
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 shows that the standard Newton algorithm performs very poorly.
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 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 7 shows that the Levenberg-Marquardt algorithm outperforms Algorithm 3. However, the robustness remains low for angles close the edges of the yield surfaces.
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 shows that using the Jacobi eigen solver has a strong influence on the robustness of Algorithm 3.
Using the Jacobi eigen solver with Algorithm 4 shows that convergence is reached for all tested values of \(\alpha\) and \(x\).
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.