This article discusses the implementation of the Yld2004-18p behaviour, as decribed by Barlat et al. (See [1]).
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 given by: \[ {{\displaystyle \frac{\displaystyle \sigma_{\mathrm{eq}}^{B}}{\displaystyle \sigma_{0}}}}-\left(1+{{\displaystyle \frac{\displaystyle p}{\displaystyle p_{0}}}}\right)^{np}=0 \] where:
The Barlat equivalent stress is defined as follows (See [1]): \[ \sigma_{\mathrm{eq}}^{B}= \sqrt[a]{ \frac{1}{4}\left( \sum_{i=0}^{3} \sum_{j=0}^{3} {\left|s'_{i}-s''_{j}\right|}^{a} \right) } \]
where \(s'_{i}\) and \(s''_{i}\) are the eigenvalues of two transformed stresses \(\underline{s}'\) and \(\underline{s}''\) by two linear transformation \(\underline{\underline{\mathbf{L}}}'\) and \(\underline{\underline{\mathbf{L}}}''\): \[ \left\{ \begin{aligned} \underline{s}' &= \underline{\underline{\mathbf{L'}}} \,\colon\,\underline{\sigma}\\ \underline{s}'' &= \underline{\underline{\mathbf{L''}}}\,\colon\,\underline{\sigma}\\ \end{aligned} \right. \]
The linear transformations \(\underline{\underline{\mathbf{L}}}'\) and \(\underline{\underline{\mathbf{L}}}''\) are defined by \(9\) coefficients (each) which describe the material orthotropy. There are defined through auxiliary linear transformations \(\underline{\underline{\mathbf{C}}}'\) and \(\underline{\underline{\mathbf{C}}}''\) as follows: \[ \begin{aligned} \underline{\underline{\mathbf{L}}}' &=\underline{\underline{\mathbf{C}}}'\,\colon\,\underline{\underline{\mathbf{M}}} \\ \underline{\underline{\mathbf{L}}}''&=\underline{\underline{\mathbf{C}}}''\,\colon\,\underline{\underline{\mathbf{M}}} \end{aligned} \] where \(\underline{\underline{\mathbf{M}}}\) is the transformation of the stress to its deviator: \[ \underline{\underline{\mathbf{M}}}=\underline{\underline{\mathbf{I}}}-{{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}}\underline{I}\,\otimes\,\underline{I} \]
The linear transformations \(\underline{\underline{\mathbf{C}}}'\) and \(\underline{\underline{\mathbf{C}}}''\) of the deviator stress are defined as follows: \[ \underline{\underline{\mathbf{C}}}'= \begin{pmatrix} 0 & -c'_{12} & -c'_{13} & 0 & 0 & 0 \\ -c'_{21} & 0 & -c'_{23} & 0 & 0 & 0 \\ -c'_{31} & -c'_{32} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & c'_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & c'_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & c'_{66} \\ \end{pmatrix} \quad \text{and} \quad \underline{\underline{\mathbf{C}}}''= \begin{pmatrix} 0 & -c''_{12} & -c''_{13} & 0 & 0 & 0 \\ -c''_{21} & 0 & -c''_{23} & 0 & 0 & 0 \\ -c''_{31} & -c''_{32} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & c''_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & c''_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & c''_{66} \\ \end{pmatrix} \]
When all the coefficients \(c'_{ji}\) and \(c''_{ji}\) are equal to \(1\), the Barlat equivalent stress reduces to the Hosford equivalent stress.
The behaviour is integrated by an implicit scheme. Two state variables are introduced:
The elastic strain is automatically defined by the StandardElasticity
brick.
First, an elastic prediction of the stress \(\underline{\sigma}^{\mathrm{tr}}\) is made (The following expression is not valid in plane stress hypothesis, this is the reason why we will use the computeElasticPrediction
method defined by the StandardElasticity
brick, 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}^{B}\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}^{B}\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}^{B}\right|_{t+\theta\,\Delta\,t}}\\ \end{aligned} \right. \]
To determine the equivalent plastic strain increment, the following equation must be satisfied: \[ f_{p}={{\displaystyle \frac{\displaystyle {\left.\sigma_{\mathrm{eq}}^{B}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \sigma 0}}}-\left(1+{{\displaystyle \frac{\displaystyle {\left.p\right|_{t+\theta\,\Delta\,t}}}{\displaystyle p_{0}}}}\right)^{np}=0 \]
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}}}} &= {{\displaystyle \frac{\displaystyle 2\,\mu\,\theta}{\displaystyle \sigma_{0}}}}\,{\left.\underline{n}^{B}\right|_{t+\theta\,\Delta\,t}}\\ {\displaystyle \frac{\displaystyle \partial f_{p}}{\displaystyle \partial \Delta\,p}} &= -{{\displaystyle \frac{\displaystyle np\,\theta}{\displaystyle p0}}}\left(1+{{\displaystyle \frac{\displaystyle {\left.p\right|_{t+\theta\,\Delta\,t}}}{\displaystyle p_{0}}}}\right)^{np-1}\\ \end{aligned} \right. \]
The beginning of the file gives some information about the behaviour:
@DSL
keyword.@Behaviour
keyword.@Author
).@Description
).An orthotropic axes convention has to be defined to have consistent definition of the linear transformations \(l_{1}\) and \(l_{2}\). The behaviour has been identified for plates, so the natural choice is the Plate
orthotropic convention:
@OrthotropicBehaviour<Plate>;
This choice limits the number of possible modelling hypotheses that are available to the following list: Tridimensional
, PlaneStrain
, GeneralisedPlaneStrain
, PlaneStress
.
Thanks to the StandardElasticity
brick defined below, all those modelling hypotheses can be supported, but the support for the plane stress hypothesis as to be explicitly requested. The following statement, starting with the @ModellingHypotheses
, enables all the previously listed modelling hypotheses:
@ModellingHypotheses {Tridimensional,PlaneStrain,
GeneralisedPlaneStrain,PlaneStress};
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 material properties are hard-coded. The @ElasticMaterialProperties
is used to declare the Young modulus and the Poisson ratio.
@ElasticMaterialProperties {150e9,0.3};
This keyword automatically declares two parameters called YoungModulus
and PoissonRatio
.
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;
// 646*0.025**0.227
@Parameter sig0 = 150e6;
"HardeningReferenceStress");
sig0.setEntryName(@Parameter p0 = 1.e-4;
"HardeningReferenceStrain");
p0.setEntryName(@Parameter np = 2;
"BarlatExponent");
np.setEntryName(@Parameter c1p[9] = {-0.069888,0.079143,0.936408,
0.524741,1.00306,1.36318,
0.954322,1.06906,1.02377};
"BarlatLinearTransformationCoefficients1");
c1p.setEntryName(@Parameter c2p[9] = {0.981171,0.575316,0.476741,
1.14501,0.866827,-0.079294,
1.40462,1.1471,1.05166};
"BarlatLinearTransformationCoefficients2"); c2p.setEntryName(
Here a
stands for the Barlat exponent, sig0
, p0
and np
are the parameters describing the hardening of the material.
The parameters associated with the first linear transformation are given in the following order (see this page for details):
\[ \left(c^{1}_{xy},c^{1}_{yx},c^{1}_{xz},c^{1}_{zx},c^{1}_{yz},c^{1}_{zy},c^{1}_{xy},c^{1}_{xz},c^{1}_{yz}\right) \]
Note In his paper, Barlat and coworkers seems to use the following convention for storing symmetric tensors:
\[ \begin{pmatrix} xx & yy & zz & yz & zx & xy \end{pmatrix} \]
which is not consistent with the
TFEL
/Cast3M
/Abaqus
/Ansys
conventions:\[ \begin{pmatrix} xx & yy & zz & xy & xz & yz \end{pmatrix} \]
Therefore, if one wants to uses coeficients \(c^{B}\) given by Barlat et al., one shall “swap” the appropriate coefficients.
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(
For the implementations, we will need three local variables:
b
. This boolean will be true
if a plastic loading occurs.l1
and l2
that stands for the linear transformation of the stress tensors.@LocalVariable bool b;
@LocalVariable Stensor4 l1,l2;
The main goal of the local variable initialization is to test if the elastic prediction of the stress lies inside the yield surface.
@InitializeLocalVariables{
// when using the `Plate` orthotropic axes convention, all the
// modelling hypotheses uses the same convention, so we can just
// use the `makeBarlatLinearTransformation` function with the
// space dimension N as the template parameter
0],c1p[1],c1p[2],
l1 = makeBarlatLinearTransformation<N>(c1p[3],c1p[4],c1p[5],
c1p[6],c1p[7],c1p[8]);
c1p[0],c2p[1],c2p[2],
l2 = makeBarlatLinearTransformation<N>(c2p[3],c2p[4],c2p[5],
c2p[6],c2p[7],c2p[8]);
c2p[const stress seps = 1.e-10*young;
const auto sigel = computeElasticPrediction();
const auto seqel = computeBarlatStress(sigel,l1,l2,a,seps);
1+p/p0,np);
b = seqel/sig0>pow( }
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 makeBarlatLinearTransformation
and the computeBarlatStress
functions are described in details here
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;
}const auto p_ = p+theta*dp;
if(p_<0){
return false;
}
real seq;
Stensor n;
Stensor4 dn;const auto tmp = pow(1+p_/p0,np);
std::tie(seq,n,dn) = computeBarlatStressSecondDerivative(sig,l1,l2,a,seps);
feel += dp*n;2*mu*theta*dp*dn;
dfeel_ddeel +=
dfeel_ddp = n;
fp = seq/sig0-tmp;2*mu*theta)/sig0*n;
dfp_ddeel = (p_)*tmp;
dfp_ddp = (np*theta)/(p0+ }
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 computeBarlatStressSecondDerivative
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.
In C++17
, the previous code can be more simplier and smaller, using structured bindings:
@Integrator{
const stress seps = 1.e-10*young;
if(!b){
// elastic loading, nothing to be done
return true;
}const auto p_ = p+theta*dp;
if(p_<0){
return false;
}const auto tmp = pow(1+p_/p0,np);
const auto [seq,n,dn] = computeBarlatStressSecondDerivative(sig,l1,l2,a,seps);
feel += dp*n;2*mu*theta*dp*dn;
dfeel_ddeel +=
dfeel_ddp = n;
fp = seq/sig0-tmp;2*mu*theta)/sig0*n;
dfp_ddeel = (p_)*tmp;
dfp_ddp = (np*theta)/(p0+ }