This paper is dedicated to the implementation a multi-surface, compressible, perfect plastic behaviour.

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 first surface is based on the Drucker-Prager yield criterion defined as follows: \[ \sigma_{\mathrm{eq}}^{\mathrm{DP}}=q+\tan{\left(\beta\right)}\,p \] where:
  • \(p\) is the hydrostatique pressure: \[ p = {{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}}\,{\mathrm{tr}{\left(\underline{\sigma}\right)}} \]
  • \(q\) is the von Mises norm of the stress tensor \(\underline{\sigma}\): \[ q = \sigma_{\mathrm{eq}}= \sqrt{3\,J_{2}} = \sqrt{{{\displaystyle \frac{\displaystyle 3}{\displaystyle 2}}}\,\underline{s}\,\colon\,\underline{s}} \] where \(\underline{s}\) is the deviatoric part of the stress tensor and \(J_{2}\) the second invariant of \(\underline{s}\). The plastic surface is defined using a constant yield stress \(R_{0}^{\mathrm{DP}}\): \[ \sigma_{\mathrm{eq}}^{\mathrm{DP}}-R_{0}^{\mathrm{DP}}=0 \] This surface is a line in the \({\left(p,q\right)}\) plane.
• The second surface forms a so-called cap which closes the first plastic surface. This plastic surface is based on a criterion close the Green one (See [1] and this page for details): \[ \sigma_{\mathrm{eq}}^{c}=\sqrt{{\left(p-p_{a}\right)}^{2}+{\left(R\,q\right)}^{2}} \] Again, the plastic surface is defined using a constant yield stress \(R_{0}^{c}\), as follows: \[ \sigma_{\mathrm{eq}}^{c}-R_{0}^{c}=0 \]
  This surface is an ellipse in the \({\left(p,q\right)}\) plane. The material parameter \(p_{a}\) marks the transition between the two plastic surface. For the two surface to intersect at \(p_{a}\), the following relation must hold: \[ R_{0}^{c}=R\,{\left(R_{0}^{\mathrm{DP}}-p_{a}\,\tan{\beta}\right)} \] The material constants associated with the cap are \(p_{a}\) and the ellipse exentricity \(R\). Here, we prefer another material constant, denoted \(p_{b}\), which stands for the minimal pressure allowed. \(p_{b}\) is related to \(R\) by: \[ R = {{\displaystyle \frac{\displaystyle p_{a} - p_{b}}{\displaystyle d - p_{a} \,\tan{\left(\beta\right)}}}}; \]

Note

Here, the definition of \(p\) differs from some authors who use the opposite. With our definition, \(p\) is negative in tri-axial compression (\({\mathrm{tr}{\left(\underline{\sigma}\right)}}<0\)).

Both criterion \(\sigma_{\mathrm{eq}}^{\mathrm{DP}}\) and \(\sigma_{\mathrm{eq}}^{c}\) are homogeneous function of degree 1.

In the following, the flow rule is assumed associated. For each mechanism \(i\), the plastic strain rate \(\underline{\dot{\varepsilon}}^{\mathrm{p,i}}\) is given by: \[ \underline{\dot{\varepsilon}}^{\mathrm{p,i}}=\dot{p}^{i}\,{\displaystyle \frac{\displaystyle \partial \sigma_{\mathrm{eq}}^{i}}{\displaystyle \partial \underline{\sigma}}}=\dot{p}^{i}\,\underline{n}^{i} \] where \(\dot{p}^{i}\) is the rate of the equivalent plastic strain \(p^{i}\). The expression of the normals \(\underline{n}^{i}\) is given below.

Derivatives of the stress criteria

Derivatives of the invariants

In the following, we will need the first and second derivatives of \(p\) and \(q\) with respect to \(\underline{\sigma}\). One may refer to this page for more details.

The first derivative of \(p\) with respect to \(\underline{\sigma}\) is: \[ {\displaystyle \frac{\displaystyle \partial p}{\displaystyle \partial \underline{\sigma}}} = {{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}}\,\underline{I} \]

The first derivative of \(q\) with respect to \(\underline{\sigma}\) will be denoted \(\underline{n}_{q}\) and is given by: \[ {\displaystyle \frac{\displaystyle \partial q}{\displaystyle \partial \underline{\sigma}}} = \underline{n}_{q} = {{\displaystyle \frac{\displaystyle 3}{\displaystyle 2\,q}}}\underline{s} \]

The second derivative of \(q\) with respect to \(\underline{\sigma}\) can be computed as follows: \[ {\displaystyle \frac{\displaystyle \partial^{2} q}{\displaystyle \partial \underline{\sigma}^{2}}} = {{\displaystyle \frac{\displaystyle 1}{\displaystyle q}}}{\left(\underline{\underline{\mathbf{M}}}-\underline{n}_{q}\,\otimes\,\underline{n}_{q}\right)} \] where \(\underline{\underline{\mathbf{M}}}\) is given by: \[ \underline{\underline{\mathbf{M}}}={{\displaystyle \frac{\displaystyle 3}{\displaystyle 2}}}\,{\left(\underline{\underline{\mathbf{I}}}-{{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}}\,\underline{I}\,\otimes\,\underline{I}\right)} \]

Derivatives of \(\sigma_{\mathrm{eq}}^{\mathrm{DP}}\)

Using the previous results, the first derivative of \(\sigma_{\mathrm{eq}}^{\mathrm{DP}}\) with respect to \(\underline{\sigma}\) will be denoted \(\underline{n}^{\mathrm{DP}}\) is:

\[ {\displaystyle \frac{\displaystyle \partial \sigma_{\mathrm{eq}}^{\mathrm{DP}}}{\displaystyle \partial \underline{\sigma}}} = \underline{n}^{\mathrm{DP}}= \underline{n}_{q}+{{\displaystyle \frac{\displaystyle \tan{\left(\beta\right)}}{\displaystyle 3}}}\,\underline{I} \]

The second derivative of \(\sigma_{\mathrm{eq}}^{\mathrm{DP}}\) with respect to \(\underline{\sigma}\) is: \[ {\displaystyle \frac{\displaystyle \partial^{2} \sigma_{\mathrm{eq}}^{\mathrm{DP}}}{\displaystyle \partial \underline{\sigma}^{2}}} = {\displaystyle \frac{\displaystyle \partial \underline{n}^{\mathrm{DP}}}{\displaystyle \partial \underline{\sigma}}}= {\displaystyle \frac{\displaystyle \partial^{2} q}{\displaystyle \partial \underline{\sigma}^{2}}} = {{\displaystyle \frac{\displaystyle 1}{\displaystyle q}}}{\left(\underline{\underline{\mathbf{M}}}-\underline{n}_{q}\,\otimes\,\underline{n}_{q}\right)} \]

Derivatives of \(\sigma_{\mathrm{eq}}^{c}\)

The first derivative of \(\sigma_{\mathrm{eq}}^{\mathrm{DP}}\) with respect to \(\underline{\sigma}\) will be denoted \(\underline{n}^{c}\) and can be computed as follows:

\[ {\displaystyle \frac{\displaystyle \partial \sigma_{\mathrm{eq}}^{\mathrm{DP}}}{\displaystyle \partial \underline{\sigma}}} = \underline{n}^{c}= {{\displaystyle \frac{\displaystyle 1}{\displaystyle \sigma_{\mathrm{eq}}^{c}}}}{\left( {{\displaystyle \frac{\displaystyle {\left(p-p_{a}\right)}}{\displaystyle 3}}}\underline{I}+R^{2}\,q\,\underline{n}_{q} \right)} \]

The second derivative of \(\sigma_{\mathrm{eq}}^{c}\) with respect to \(\underline{\sigma}\) is: \[ {\displaystyle \frac{\displaystyle \partial^{2} \sigma_{\mathrm{eq}}^{c}}{\displaystyle \partial \underline{\sigma}^{2}}}= {\displaystyle \frac{\displaystyle \partial \underline{n}^{c}}{\displaystyle \partial \underline{\sigma}}}= {{\displaystyle \frac{\displaystyle 1}{\displaystyle \sigma_{\mathrm{eq}}^{c}}}}{\left( {{\displaystyle \frac{\displaystyle 1}{\displaystyle 9}}}\,\underline{I}\,\otimes\,\underline{I}+ R^{2}\,\underline{n}_{q}\,\otimes\,\underline{n}_{q}+ R^{2}\,{\left(\underline{\underline{\mathbf{M}}}-\underline{n}_{q}\,\otimes\,\underline{n}_{q}\right)} -\underline{n}_{c}\,\otimes\,\underline{n}_{c} \right)} \]

Integration scheme

The behaviour is integrated by an implicit \(\theta\)-scheme. This scheme requires to choose the state variables and to associate to each variable an equation.

Choice of the state variables

Three state variables are introduced:

• the elastic strain \(\underline{\varepsilon}^{\mathrm{el}}\).
• the equivalent plastic strain \(p^{\mathrm{DP}}\).
• the equivalent plastic strain \(p^{c}\) associated with the cap.

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

\(p^{\mathrm{DP}}\) and \(p^{c}\)could be considered as an integration variables, but, for post-processing purposes, we choose to keep it as a state variables.

Implicit scheme

Let \({\left.\vec{Y}\right|_{t}}\) be a vector containing the values of the state variables at the beginning of the time step and \(\Delta\,\vec{Y}\) a vector holding the values of unknown increment of those state variables. Here, \(\vec{Y}\) and \(\Delta\,\vec{Y}\) have the the following form:

\[ {\left.\vec{Y}\right|_{t}} = \begin{pmatrix} {\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t}} \\ {\left.p^{\mathrm{DP}}\right|_{t}} \\ {\left.p^{c}\right|_{t}} \\ \end{pmatrix} \quad \Delta\,\vec{Y}= \begin{pmatrix} \Delta\,{\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t}} \\ \Delta\,{\left.p^{\mathrm{DP}}\right|_{t}} \\ \Delta\,{\left.p^{c}\right|_{t}} \\ \end{pmatrix} \]

The principle of implicit schemes is to discretize the constitutive equations so that the increment of the unknowns \(\Delta\,\vec{Y}\) satisfy the following equation: \[ \vec{F}{\left(\Delta\,\vec{Y}\right)}=0 \] where \(\vec{F}\) is a vectorial function.

This (non linear) equation is solved iteratively by the Newton-Raphson method. This method requires the jacobian \({\displaystyle \frac{\displaystyle \partial \vec{F}}{\displaystyle \partial \vec{Y}}}\) to be defined.

In the following, \(\vec{F}\) will be decomposed as follows: \[ \vec{F}{\left(\Delta\,\vec{Y}\right)}= \begin{pmatrix} f_{\underline{\varepsilon}^{\mathrm{el}}} \\ f_{p^{\mathrm{DP}}} \\ f_{p^{c}} \\ \end{pmatrix} \]

The jacobian will be decomposed by blocks as follows: \[ {\displaystyle \frac{\displaystyle \partial \vec{F}}{\displaystyle \partial \Delta\,\vec{Y}}}= \begin{pmatrix} {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} & {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,p^{\mathrm{DP}}}} & {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,p^{c}}} \\ {\displaystyle \frac{\displaystyle \partial f_{p^{\mathrm{DP}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} & {\displaystyle \frac{\displaystyle \partial f_{p^{\mathrm{DP}}}}{\displaystyle \partial \Delta\,p^{\mathrm{DP}}}} & {\displaystyle \frac{\displaystyle \partial f_{p^{\mathrm{DP}}}}{\displaystyle \partial \Delta\,p^{c}}} \\ {\displaystyle \frac{\displaystyle \partial f_{p^{c}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} & {\displaystyle \frac{\displaystyle \partial f_{p^{c}}}{\displaystyle \partial \Delta\,p^{\mathrm{DP}}}} & {\displaystyle \frac{\displaystyle \partial f_{p^{c}}}{\displaystyle \partial \Delta\,p^{c}}} \\ \end{pmatrix} \]

In the following, we will use the following notation: \[ {\left.\vec{Y}\right|_{t+\theta\,\Delta\,t}}={\left.\vec{Y}\right|_{t}}+\theta\,\Delta\,\vec{Y} \] where \(\theta\) is a numerical parameter.

By extension, the value of every function \(f\) of the state variables that is evaluated with \({\left.\vec{Y}\right|_{t+\theta\,\Delta\,t}}\) as argument is denoted \({\left.f\right|_{t+\theta\,\Delta\,t}}\). For example, the stress tensor \({\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\) is evaluated as follows \[ {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}={\left.D\right|_{t+\theta\,\Delta\,t}}\,\colon\,{\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}} \] where \({\left.D\right|_{t+\theta\,\Delta\,t}}\) is the evaluation of the stiffness tensor at \(t+\theta\,\Delta\,t\). In the example treated here, the stiffness tensor is constant.

Status

The integration scheme will be based on a set of three implicit equations:

• The first equation is associated with split of the strain.
• The second equation is associated with the first plastic mechanism (Drucker-Prager).
• The third equation is associated with the second plastic mechanism (cap).

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.

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)} \]

• If the predicted stress is inside the elastic domain, no plastic flow occurs.
• Otherwise, the material state at the end of the time step lies on at least one of the yield surfaces.

As describe in the previous paragraph, this elastic prediction is used initialize the status associated with each mechanism.

Equation associated with the elastic strain

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^{\mathrm{DP}}\,{\left.\underline{n}^{\mathrm{DP}}\right|_{t+\theta\,\Delta\,t}}+ \Delta\,p^{c}\,{\left.\underline{n}^{c}\right|_{t+\theta\,\Delta\,t}} \]

The derivatives of this equation with respect to \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}\), \(\Delta\,p^{\mathrm{DP}}\) and \(\Delta\,p^{c}\) 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}}} +\theta\,\Delta\,p^{\mathrm{DP}}\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{n}^{\mathrm{DP}}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\,\colon\,{\left.\underline{\underline{\mathbf{D}}}\right|_{t+\theta\,\Delta\,t}} +\theta\,\Delta\,p^{c}\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{n}^{c}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\,\colon\,{\left.\underline{\underline{\mathbf{D}}}\right|_{t+\theta\,\Delta\,t}}\\ {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,p^{\mathrm{DP}}}} &= {\left.\underline{n}^{\mathrm{DP}}\right|_{t+\theta\,\Delta\,t}}\\ {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,p^{c}}} &= {\left.\underline{n}^{c}\right|_{t+\theta\,\Delta\,t}}\\ \end{aligned} \right. \]

where \({\left.\underline{\underline{\mathbf{D}}}\right|_{t+\theta\,\Delta\,t}}\) is the stiffness tensor. For an isotropic material, \({\left.\underline{\underline{\mathbf{D}}}\right|_{t+\theta\,\Delta\,t}}\) is equal to:

\[ {\left.\underline{\underline{\mathbf{D}}}\right|_{t+\theta\,\Delta\,t}}={\left.\lambda\right|_{t+\theta\,\Delta\,t}}\,\underline{I}\,\otimes\,\underline{I}+2\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\,\underline{\underline{\mathbf{I}}} \]

Equation associated with the Drucker-Prager plastic flow

If the Drucker-Prager mechanism is active, the material must be on the yield surface. Thus, the equation associated with the mechanism is:

\[ f_{p^{\mathrm{DP}}} = {{\displaystyle \frac{\displaystyle 1}{\displaystyle {\left.E\right|_{t+\theta\,\Delta\,t}}}}}\,{\left({\left.\sigma_{\mathrm{eq}}^{\mathrm{DP}}\right|_{t+\theta\,\Delta\,t}}-R_{0}^{\mathrm{DP}}\right)} \]

Here, the Young modulus \(E\) has been used to normalize the equation.

The derivatives of this equation with respect to \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}\), \(\Delta\,p^{\mathrm{DP}}\) are given by:

\[ \left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial f_{p^{\mathrm{DP}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} &= {{\displaystyle \frac{\displaystyle \theta}{\displaystyle E}}}\,{\left.\underline{n}^{\mathrm{DP}}\right|_{t+\theta\,\Delta\,t}}\,\colon\,{\left.\underline{\underline{\mathbf{D}}}\right|_{t+\theta\,\Delta\,t}}\\ {\displaystyle \frac{\displaystyle \partial f_{p^{\mathrm{DP}}}}{\displaystyle \partial \Delta\,p^{\mathrm{DP}}}} &= 0\\ \end{aligned} \right. \]

If the Drucker-Prager mechanism is not active, the equation assoaciated with \(p^{\mathrm{DP}}\) is simply:

\[ f_{p^{\mathrm{DP}}}=\Delta\,p^{\mathrm{DP}} \]

Equation associated with the Drucker-Prager plastic flow

If the cap mechanism is active, the material must be on the yield surface. Thus, the equation associated with the mechanism is:

\[ f_{p^{c}} = {{\displaystyle \frac{\displaystyle 1}{\displaystyle {\left.E\right|_{t+\theta\,\Delta\,t}}}}}\,{\left({\left.\sigma_{\mathrm{eq}}^{c}\right|_{t+\theta\,\Delta\,t}}-R_{0}^{c}\right)} \]

Here, the Young modulus \(E\) has been used to normalize the equation.

The derivatives of this equation with respect to \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}\), \(\Delta\,p^{c}\) are given by:

\[ \left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial f_{p^{c}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} &= {{\displaystyle \frac{\displaystyle \theta}{\displaystyle E}}}\,{\left.\underline{n}^{c}\right|_{t+\theta\,\Delta\,t}}\,\colon\,{\left.\underline{\underline{\mathbf{D}}}\right|_{t+\theta\,\Delta\,t}}\\ {\displaystyle \frac{\displaystyle \partial f_{p^{c}}}{\displaystyle \partial \Delta\,p^{c}}} &= 0\\ \end{aligned} \right. \]

If the cap mechanism is not active, the equation assoaciated with \(p^{c}\) is simply:

\[ f_{p^{c}}=\Delta\,p^{c} \]

Implementation

Metadata

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

• the integration scheme used, selected by the @DSL keyword.
• the name of the behaviour, introduced by the @Behaviour keyword.
• the author of the implementation (@Author).
• a small description of the behaviour (@Description).

@DSL Implicit;
@Behaviour DruckerPragerCap;
@Author Thomas Helfer;
@Description{
  A simple implementation of a perfect
  plasticity behaviour using the
  Drucker-Prager yield criterion
  closed by a cap.
};

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:

• Automatic computation of the stress tensor at various stages of the behaviour integration.
• Automatic computation of the consistent tangent operator.
• Automatic support for plane stress and generalized plane stress modelling hypotheses (The axial strain is defined as an additional state variable and the associated equation in the implicit system is added to enforce the plane stess condition).
• Automatic addition of the standard terms associated with the elastic strain state variable.

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 an array of two states variables named p:

@StateVariable strain p[2];
p.setGlossaryName("EquivalentPlasticStrain");

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.

Statuses

@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

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;

Material constants

Material constants associated with the elastic part of the behaviour

@ElasticMaterialProperties {150e9,0.3};

Material constants associated with the plastic parts of the behaviour

The material constants associated with Drucker-Prager yield surface are \(\beta\) and \(R_{0}^{\mathrm{DP}}\). In the following, \(R_{0}^{\mathrm{DP}}\), will be called d.

@Parameter d    = 150e6;
@Parameter beta = 0.6;

@Parameter pa   = -10e6;
@Parameter pb   = -75e6;

Initialisation of the local variables

The @InitLocalVariables block is called once at the very beginning of the behaviour integration.

@InitLocalVariables {
  tg = tan(beta);
  R = (pa - pb) / (d - pa * tg);
}

Elastic prediction, initial value of the statuses

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);
  bpl[0] = (qel + prel * tg - d > 0) && (prel > pa);
  // Cap
  const auto seq_c = sqrt(power<2>(prel - pa) + power<2>(R * qel));
  bpl[1] = (seq_c > (pa - pb)) && (prel <= pa);
}

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).

Implicit system

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 \(p\), \(q\) and their derivatives:

    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);
      feel += dp[0]*n;
      fp(0) = (q+pr*tg-d)/young;
      // jacobian
      dfeel_ddeel += 2*mu*theta*dp[0]*dnq;
      dfeel_ddp(0) = n;
      dfp_ddp(0,0) = 0;
      dfp_ddeel(0) = (2*mu*theta/young)*n+(tg*theta*lambda/young)*id;
    }

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;
      feel += dp[1] * n;
      fp(1) = (seq_c - R * (d - pa * tg)) / young;
      // jacobian
      dfeel_ddeel += theta * dp[1] * dn * De;
      dfeel_ddp(1) = n;
      if(std::abs(seq_c)<seps){
        dfp_ddp(1, 1) = 1;
      } else {
        dfp_ddp(1, 1) = 0;
      }
      dfp_ddeel(1) = theta * (n | De) / young;
    }

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

A posterio checks: validation of the assumptions

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 \(\Delta\,p_{i}\) is positive. Otherwise, the mechanism is desactivated and the 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
        converged = bpl[0] = false;
      }
    } else {
      const auto pr = trace(sig) / 3;
      const auto q = sigmaeq(sig);
      if (q + pr * tg > d) {
        converged = false;
        bpl[0] = true;
      }
    }
    if (bpl[1]) {
      if (dp[1] < 0) {
        // desactivating this system
        converged = bpl[1] = false;
      }
     } 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)) {
         converged = false;
         bpl[1] = true;
       }
     }
  }
}

Evolution of the porosity

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 \(f\) is defined as the ratio between the current mass density \(\rho\) and the theoretical density \(\rho^{\mathrm{th}}\): \[ f=1-{{\displaystyle \frac{\displaystyle \rho}{\displaystyle \rho^{\mathrm{th}}}}} \]

The mass conservation reads: \[ \rho\,J=\rho_{0} \] where \(\rho_{0}\) is the initial density and \(J\) the change of volume.

Thus, by derivation with respect to time: \[ \dot{\rho}\,J+\rho\,\dot{J}=0\quad\Rightarrow\quad{{\displaystyle \frac{\displaystyle \dot{\rho}}{\displaystyle \rho}}}=-{{\displaystyle \frac{\displaystyle \dot{J}}{\displaystyle J}}} \]

In the Hencky strain measure is used (see the @StrainMeasure keyword), the following relation holds: \[ J=\exp{\left({\mathrm{tr}{\left(\underline{\varepsilon}^{\mathrm{to}}\right)}}\right)} \] Hence, \[ \dot{\rho}=-\rho\,{\mathrm{tr}{\left(\underline{\dot{\varepsilon}}^{\mathrm{to}}\right)}} \]

Finally, the evolution of the porosity is: \[ \dot{f}={\left(1-f\right)}\,{\mathrm{tr}{\left(\underline{\dot{\varepsilon}}^{\mathrm{to}}\right)}} \]

Following [2], the elasic part of the strain is neglected and the following evolution of \(f\) is finally kept: \[ \dot{f}={\left(1-f\right)}\,{\mathrm{tr}{\left(\underline{\dot{\varepsilon}}^{\mathrm{p}}\right)}} \]

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:

• The yield criteria and the isotropic hardening depend on the porosity. In this case, the simpliest approach is to consider the porosity as a new state variable and to add the appropriate derivatives of the yield surface in the definition of the jacobian of the implicit system.
• The yield criteria and the isotropic hardening do not depend on the porosity. In this second case, the porosity can be treated as an auxiliary state variable.

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.

Porosity as a state variable

The discretization of the evolution of the porosity leads to the following equation: \[ f_{f} = \Delta\,f-{\left(1-{\left.f\right|_{t+\theta\,\Delta\,t}}\right)}\,{\left( \Delta\,p^{\mathrm{DP}}\,{\mathrm{tr}{\left({\left.\underline{n}^{\mathrm{DP}}\right|_{t+\theta\,\Delta\,t}}\right)}} + \Delta\,p^{c}\,{\mathrm{tr}{\left({\left.\underline{n}^{c}\right|_{t+\theta\,\Delta\,t}}\right)}} \right)} \]

The derivatives of the \(f_{p}\) are trivial: \[ \left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial f_{f}}{\displaystyle \partial \Delta\,f}}&=1+ \theta\, {\left( \Delta\,p^{\mathrm{DP}}\,{\mathrm{tr}{\left({\left.\underline{n}^{\mathrm{DP}}\right|_{t+\theta\,\Delta\,t}}\right)}} + \Delta\,p^{c}\,{\mathrm{tr}{\left({\left.\underline{n}^{c}\right|_{t+\theta\,\Delta\,t}}\right)}} \right)}\\ &-{\left(1-{\left.f\right|_{t+\theta\,\Delta\,t}}\right)}\, {\left( \Delta\,p^{\mathrm{DP}}\,{\mathrm{tr}{\left({\displaystyle \frac{\displaystyle \partial {\left.\underline{n}^{\mathrm{DP}}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,f}}\right)}} + \Delta\,p^{c}\,{\mathrm{tr}{\left({\displaystyle \frac{\displaystyle \partial {\left.\underline{n}^{c}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,f}}\right)}} \right)}\\ {\displaystyle \frac{\displaystyle \partial f_{f}}{\displaystyle \partial \Delta\,p^{\mathrm{DP}}}}&= -{\left(1-{\left.f\right|_{t+\theta\,\Delta\,t}}\right)}\, {\mathrm{tr}{\left({\left.\underline{n}^{\mathrm{DP}}\right|_{t+\theta\,\Delta\,t}}\right)}}\\ {\displaystyle \frac{\displaystyle \partial f_{f}}{\displaystyle \partial \Delta\,p^{c}}}&= -{\left(1-{\left.f\right|_{t+\theta\,\Delta\,t}}\right)}\, {\mathrm{tr}{\left({\left.\underline{n}^{c}\right|_{t+\theta\,\Delta\,t}}\right)}}\\ {\displaystyle \frac{\displaystyle \partial f_{f}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}&= -{\left(1-{\left.f\right|_{t+\theta\,\Delta\,t}}\right)}\, \underline{I}\,\colon\, {\left( \Delta\,p^{\mathrm{DP}}\,{\displaystyle \frac{\displaystyle \partial \underline{n}^{\mathrm{DP}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}+ \Delta\,p^{\mathrm{c}}\,{\displaystyle \frac{\displaystyle \partial \underline{n}^{c}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} \right)}\\ \end{aligned} \right. \]

In some case, neglecting the trace of \({\displaystyle \frac{\displaystyle \partial {\left.\underline{n}^{\mathrm{DP}}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,f}}\) and \({\displaystyle \frac{\displaystyle \partial {\left.\underline{n}^{c}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,f}}\) can be a very good approximation.

Porosity as an auxiliary state variable

The evolution of the porosity can be computed once the plastic increments are known, as follows:

\[ \Delta\,f = {{\displaystyle \frac{\displaystyle {\left(1-{\left.f\right|_{t}}\right)}\,{\mathrm{tr}{\left(\Delta\,\underline{\varepsilon}^{\mathrm{p}}\right)}}}{\displaystyle 1+\theta\,{\mathrm{tr}{\left(\Delta\,\underline{\dot{\varepsilon}}^{\mathrm{p}}\right)}}}}} \]

with:

\[ {\mathrm{tr}{\left(\Delta\,\underline{\varepsilon}^{\mathrm{p}}\right)}}= \Delta\,p^{\mathrm{DP}}\,{\mathrm{tr}{\left({\left.\underline{n}^{\mathrm{DP}}\right|_{t+\theta\,\Delta\,t}}\right)}}+ \Delta\,p^{c}\,{\mathrm{tr}{\left({\left.\underline{n}^{c}\right|_{t+\theta\,\Delta\,t}}\right)}} \]

References