This paper is dedicated to the implementation a multi-surface, compressible, perfect plastic behaviour.
The elasticity is assumed linear and isotropic, i.e. given by the
Hooke law:
Note
Here, the definition of
differs from some authors who use the opposite. With our definition, is negative in tri-axial compression ( ).
Both criterion
In the following, the flow rule is assumed associated. For each
mechanism
In the following, we will need the first and second derivatives of
The first derivative of
The first derivative of
The second derivative of
Using the previous results, the first derivative of
The second derivative of
The first derivative of
The second derivative of
The behaviour is integrated by an implicit
Three state variables are introduced:
The elastic strain is automatically defined by the
Implicit
domain specific language.
Let
The principle of implicit schemes is to discretize the constitutive
equations so that the increment of the unknowns
This (non linear) equation is solved iteratively by the
Newton-Raphson method. This method requires the jacobian
In the following,
The jacobian will be decomposed by blocks as follows:
In the following, we will use the following notation:
By extension, the value of every function
The integration scheme will be based on a set of three implicit equations:
The last two equations depend on whether the associated mechanism is assumed active or not. For each systems, a first guess will be made based on an elastic prediction of the stress, as discussed in the next paragraph. After a solution to the implicit systems will be found, the validity of those two assumptions will be checked. If one of the assumption is false, the resolution will be restarted by making the opposite assumption.
In pratice, the activation of a plastic mechanism is associated to a boolean value.
First, an elastic prediction of the stress
As describe in the previous paragraph, this elastic prediction is used initialize the status associated with each mechanism.
The equation associated with the evolution of the elastic strain is
given by the split of strain:
The derivatives of this equation with respect to
where
If the Drucker-Prager mechanism is active, the material must be on the yield surface. Thus, the equation associated with the mechanism is:
Here, the Young modulus
The derivatives of this equation with respect to
If the Drucker-Prager mechanism is not active, the equation
assoaciated with
If the cap mechanism is active, the material must be on the yield surface. Thus, the equation associated with the mechanism is:
Here, the Young modulus
The derivatives of this equation with respect to
If the cap mechanism is not active, the equation assoaciated with
The beginning of the file gives some information about the behaviour:
@DSL
keyword.@Behaviour
keyword.@Author
).@Description
).@DSL Implicit;
@Behaviour DruckerPragerCap;
@Author Thomas Helfer;
@Description{
A simple implementation of a perfectusing the
plasticity behaviour -Prager yield criterion
Drucker.
closed by a cap};
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 an array of two states variables
named p
:
@StateVariable strain p[2];
.setGlossaryName("EquivalentPlasticStrain"); p
The first entry p[0]
will be associated to the
Drucker-Prager mechanism. The second entry p[1]
will be
associated to the cap mechanism.
@LocalVariable bool bpl[2];
The first entry bpl[0]
will be associated to the
Drucker-Prager mechanism. The second entry bpl[1]
will be
associated to the cap mechanism.
Local variables are helpers variables that are accessible in each code blocks during the all the integration of the behaviour:
@LocalVariable real tg;
@LocalVariable real R;
@ElasticMaterialProperties {150e9,0.3};
The material constants associated with Drucker-Prager yield surface
are d
.
@Parameter d = 150e6;
@Parameter beta = 0.6;
@Parameter pa = -10e6;
@Parameter pb = -75e6;
The @InitLocalVariables
block is called once at the very
beginning of the behaviour integration.
@InitLocalVariables {
= tan(beta);
tg = (pa - pb) / (d - pa * tg);
R }
The @Predictor
code block is called just before the
beginning of the Newton-Raphson algorihtm.
@Predictor{
const auto sigel = computeElasticPrediction();
// Drucker-Prager
const auto prel = trace(sigel) / 3;
const auto qel = sigmaeq(sigel);
[0] = (qel + prel * tg - d > 0) && (prel > pa);
bpl// Cap
const auto seq_c = sqrt(power<2>(prel - pa) + power<2>(R * qel));
[1] = (seq_c > (pa - pb)) && (prel <= pa);
bpl}
The computeElasticPrediction
method, introduced by the
StandardElasticity
brick computes an elastic prediction of
stress. This method takes the modelling hypothesis into account when
required (plane stress, generalised plane stress).
The implicit system is build inside the @Integrator
code
block:
@Integrator {
We first declare some useful variables:
constexpr const auto id = Stensor::Id();
constexpr const auto id4 = Stensor4::Id();
constexpr const auto id_x_id = Stensor4::IxI();
constexpr const auto M = Stensor4::M();
const auto seps = real(1.e-12) * young;
Those variables are declared constexpr
to allow the
compiler to evaluate them at compile-time. Those variables are declared
const
so that their values can’t be changed.
The following code tests if at least one mechanism is active.
if ((bpl[0]) || (bpl[1])) {
If at least one mechanism is active, the following code computes
const auto pr = trace(sig) / 3;
const auto q = sigmaeq(sig);
const auto iq = 1 / max(q,seps);
const auto nq = eval(3 * deviator(sig) * (iq / 2));
const auto dnq = eval((M - (nq ^ nq)) * iq);
We now treat the case where the Drucker-Prager mechanism is active.
if (bpl[0]) {
const auto n = eval(nq+(tg/3)*id);
+= dp[0]*n;
feel (0) = (q+pr*tg-d)/young;
fp// jacobian
+= 2*mu*theta*dp[0]*dnq;
dfeel_ddeel (0) = n;
dfeel_ddp(0,0) = 0;
dfp_ddp(0) = (2*mu*theta/young)*n+(tg*theta*lambda/young)*id;
dfp_ddeel}
The following code describes the implicit equations to be solved if the cap mechanism is active.
if (bpl[1]) {
const auto seq_c = sqrt(power<2>(pr - pa) + power<2>(R * q));
const auto iseq_c = 1 / max(seq_c,seps);
const auto n = eval(((pr - pa) * (id / 3) + R * R * q * nq) * iseq_c);
const auto dn = (Stensor4::IxI() / 9 + R * R * M - (n ^ n)) * iseq_c;
const auto De = 2 * mu * id4 + lambda * id_x_id;
+= dp[1] * n;
feel (1) = (seq_c - R * (d - pa * tg)) / young;
fp// jacobian
+= theta * dp[1] * dn * De;
dfeel_ddeel (1) = n;
dfeel_ddpif(std::abs(seq_c)<seps){
(1, 1) = 1;
dfp_ddp} else {
(1, 1) = 0;
dfp_ddp}
(1) = theta * (n | De) / young;
dfp_ddeel}
The following braces closes the test which was used to check that at least one mechanism is active.
} // end of if ((bpl[0]) || (bpl[1]))
Finally, the @Integrator
code block is closed:
} // end of @Integrator
The @AdditionalConvergenceChecks
code block has been
introduced in TFEL-3.1
to have a better control on the
convergence of the Newton-Raphson algorithm. In particular, this block
can be used to change the status associated with each mechanism.
Here, we will change the status associated with a mechanism once the
Newton-Raphson algorithm has converged (this is indicated by the
converged
boolean variable).
If a mechanism is active, we check that converged
flag is set to
false
.
If a mechanism is not active, we check that the final solution remain
below the yield surface. Otherwise, the mechanism is activated and the
converged
flag is set to false
.
@AdditionalConvergenceChecks {
if (converged){
if (bpl[0]) {
if (dp[0] < 0) {
// desactivating this system
= bpl[0] = false;
converged }
} else {
const auto pr = trace(sig) / 3;
const auto q = sigmaeq(sig);
if (q + pr * tg > d) {
= false;
converged [0] = true;
bpl}
}
if (bpl[1]) {
if (dp[1] < 0) {
// desactivating this system
= bpl[1] = false;
converged }
} else {
const auto pr = trace(sig) / 3;
const auto q = sigmaeq(sig);
const auto seq_c = sqrt(power<2>(pr - pa) + power<2>(R * q));
if ((seq_c > pa - pb) && (pr < pa)) {
= false;
converged [1] = true;
bpl}
}
}
}
In this paragraph, we will discuss how the evolution of the porosity of the material du to the plastic flow can be added.
The porosity
The mass conservation reads:
Thus, by derivation with respect to time:
In the Hencky strain measure is used (see the
@StrainMeasure
keyword), the following relation holds:
Finally, the evolution of the porosity is:
Following [2], the elasic part of
the strain is neglected and the following evolution of
Dependency of the elastic material properties with the porosity
For simplicity, we don’t treat here the case where the elastic properties depend on the porosity. Indeed, in this case, the usage of the
StandardElasticity
brick is not more justified. In pratice, however, neglecting the additional terms du to such a dependency may work quite well.
Assuming no dependency of elastic material properties with the porosity, two cases can be distinguished:
In some cases, the second case may also be appropriate even if the yield criteria and the isotropic hardening depend on the porosity, but the evolution of the porosity during the time step has a negligible impact.
The discretization of the evolution of the porosity leads to the
following equation:
The derivatives of the
In some case, neglecting the trace of
The evolution of the porosity can be computed once the plastic increments are known, as follows:
with: