- 1 Introduction
- 2 General framework
- 3 Implicit schemes
- 4 Design choices
- 4.1 Enrichment of the
`StressCriterion`

and`InelasticFlow`

interfaces - 4.2 Conditions for the brick to contribute to the porosity evolution
- 4.3 Effect of the porosity on the equivalent strain definition
- 4.4 Saving individual contributions to the porosity evolutions in dedicated auxiliary state variables
- 4.5 Taking into account the elastic contribution in the porosity growth

- 4.1 Enrichment of the
- 5 Verifications
- 6 Applications and Finite Element simulations
- 7 Conclusions
- 8 Appendix
- 8.1 Determination of the stress criterion through a local scalar Newton algorithm
- 8.2 Derivatives of common nucleation law
- 8.3 Derivatives of standard stress criteria
- 8.4 Regularization of the effective porosity
- 8.5 On the use of porous constitutive
equations with
`Cast3M`

- 8.6 Options of the
`porosity_evolution`

section

- References

Constitutive equations for porous (visco-)plasticity are used to perform ductile fracture simulations, adding the porosity as an internal state variable. Porosity evolution is driven by nucleation laws and growth due to plastic flow.

While some implementations of these models based on implicit schemes
are already available in `MFront`

(Gurson [1], Gurson-Tvergaard-Needleman [2],
…), their numerical treatment have specific issues:

- Phenomenological laws describing porosity nucleation may have threshold to be activated or inhibited.
- Material failure usually happens when the porosity reaches a critical value corresponding to the collapse of the yield surface. Detection of material failure is awkward when the equations governing the porosity evolution are included in the implicit system to be solved. For example, the estimates of the solution may exceed the critical value during the iterations while the solution of the implicit scheme may be below this critical value.

Those issues were previously overcome by considering an explicit treatment of the porosity evolution, which requires small time steps.

The purpose of this work is twofold:

- Propose a robust fully implicit algorithm which can handle thresholds of nucleation laws and material failure.
- Extend the
`StandardElastoViscoPlasticity`

brick to this class of behaviours, allowing a declarative syntax similar as:

```
@Brick "StandardElastoViscoPlasticity" {
: "Hooke" {
stress_potential : 200e9,
young_modulus : 0.3
poisson_ratio },
: "Plastic" {
inelastic_flow : "GursonTvergaardNeedleman" {
criterion : 1.5,
q1 : 1.0,
q2 : 2.2,
q3 : 0.01,
fc : 0.1
fr },
: "Linear" {
isotropic_hardening : 150e6
R0 }
},
: {
porosity_evolution : "Chu_Needleman" {
nucleation_model : 0.01,
An : 0.1,
pn : 0.1
sn },
: "StandardPlasticModel"
growth_model };
```

NoteFor backward compatibility, effect of the porosity on plastic flows will not be taken into account if no`porosity_evolution`

section is declared.

This document is meant to be part of the MFront documentation. It will be available on the MFront website and on the ResearchGate page of the project. The authors have carefully checked that no unpublished material data has been used in the examples.

Section 2 provides a description of several stress criteria of interest for porous (visco-)plasticity and some common nucleation laws.

Section 3 describes two implicit schemes suitable for the integration of the previous constitutive laws:

- The first one is the standard implicit scheme where the implicit equations associated with all the state variables are solved all at once using a Newton-like algorithm.
- The second one is built on a staggered approach in which the porosity evolution and the evolution of the other state variables are solved independently. Both evolutions are repeated until the implicit equations associated with all the state variables are all satisfied. We show that, if required, an exact consistent tangent operator can be derived.

Section 4 discusses the design choices made for the extension of
`MFront`

’ `StandardElastoViscoPlasticity`

brick.

Section 5 is devoted to unit tests and verifications.

Section 6 describes early results that proves the robustness of the proposed schemes.

Appendix 8.1 is dedicated to the computation of the stress criterion which may require solving a non linear equation by Newton algorithm.

Appendix 8.2 gives the derivatives of some common nucleation laws.

Appendix 8.3 gives the derivatives of some common porous stress criteria.

Appendix 8.4 describes a regularization of the effective porosity in the Gurson-Tvergaard-Needleman stress criterion.

Appendix 8.5 details some numerical tricks to improve the robustness
of finite element simulations performed with `Cast3M`

using
the constitutive laws from the extended `StandardElastoViscoPlasticity`

brick.

Appendix 8.6 is dedicated to the description of the various options
that can be declared in the `porosity_evolution`

section
inside the declaration of the `StandardElastoViscoPlasticity`

brick.

In the framework adopted in this report, a porous plastic model is defined by a stress criterion \(\sigma_{\star}\) which depends explicitly on the stress tensor \(\underline{\sigma}\) and the porosity \(f\):

\[ \sigma_{\star}= \phi(\underline{\sigma},f) \]

This stress criterion is used to define the yield surface or the intensity of the viscoplastic flow.

Concerning the flow rule, two main classes of models are classically found in the literature. In the first class of model, the porosity does not change the expression of the flow rule:

\[ \underline{\dot{\varepsilon}}^{\mathrm{p}}= \dot{\lambda}\, {\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}} = \dot{p}\, {\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}} \]

i.e. the equivalent (visco-)plastic strain \(p\) is defined as the macroscopic equivalent (visco-)plastic strain and can still be identified as the plastic multiplier \(\lambda\).

This class of models encompasses:

- the work of Ponte-Castaneda (see [3]), which is the theorical basis of several viscoplastic behaviours used in uranium dioxyde fuels [4, 5];
- the work of Rousselier [6].

In the second class, the flow rule is affected by a correction factor \(1-f\) (see [7]):

\[ \underline{\dot{\varepsilon}}^{\mathrm{p}}= \dot{\lambda}\, {\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}} = (1 - f)\, \dot{p}\, {\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}} \]

where \(f\) is the porosity. Here, \(p\) can be identified to the equivalent (visco-)plastic strain in the matrix and it is no more equal to the plastic multiplier \(\lambda\).

NoteThe term \((1-f)\) comes from the definition of \(p\) as the equivalent plastic strain in the matrix through the modelling of hardening:

\[ \displaystyle{\underline{\sigma}:\underline{\dot{\varepsilon}}^{\mathrm{p}}= \dot{\lambda} \underline{\sigma}\, \colon\,{\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}}} = \dot{\lambda} \sigma_{\star}= (1 - f)R(p) \dot{p} \Rightarrow \dot{\lambda} = (1 - f) \dot{p} \]

where \(R{\left(p\right)}\) is the radius of the yield surface.

The evolution of porosity is composed of two terms: one due to plastic flow and another one due to nucleation \(A_n^i dp\):

\[\dot{f} = (1 - f){\mathrm{tr}{\left(\underline{\dot{\varepsilon}}^{\mathrm{p}}\right)}} + \sum_j A_n^j \dot{p}\]

The first term describes the incompressibility of the matrix under the assumption that the change of volume associated with the elastic part of the strain is negligible. This assumption is common, although adding the elastic contribution is fairly easy, see Section 3.1.1.

Different stress criteria implemented in the brick are described below.

**Gurson-Tvergaard-Needleman**[2]

\[ S = {\left(\dfrac{\sigma_{vM}}{\sigma_{\star}} \right)}^2 + 2 q_1 f_{\star} \cosh{{\left(\dfrac{3}{2} q_2 \dfrac{\sigma_m}{\sigma_{\star}} \right)}} - 1 - q_3 f_{\star}^2 = 0 \qquad{(1)}\] which defines implicitly \(\sigma_{\star}\), where \(\sigma_{vM}\) is the von Mises equivalent stress, and \(\sigma_m\) the mean stress. \(f_{\star}\) is defined such that:

\[ f_{\star} = \left\{ \begin{array}{ll} \delta f & \mbox{if } f < f_c \\ f_c + \delta (f - f_c) & \mbox{otherwise} \end{array} \right. \] and: \[ \delta = \left\{ \begin{array}{ll} 1 & \mbox{if } f < f_c \\ \dfrac{f_u- f_c }{f_r - f_c} & \mbox{otherwise} \end{array} \right. \] where \(f_{u}\) is the root of \(2\,q_1\,f-1-q_3\,f^{2}\).

The parameters of the model are:

\[\left\{q_1,q_2,q_3,f_c,f_r \right\} \]

For \(\left\{q_1,q_2,q_3,f_c,f_r \right\} = \left\{1,1,1,+\infty,-\right\}\), Gurson-Tvergaard-Needleman model reduces to the Gurson model [1] which has no free parameter.

Note that the physical meaning of \(q_3\) different from \(q_1^2\) is doubtful. Moreover, the current implementation is incorrect for \(q_3 > q_1^2\).

**Rousselier-Tanguy-Besson**[8]

\[ S = \dfrac{\sigma_{vM}}{(1- f )\sigma_{\star}} + \dfrac{2}{3} f D_R \exp{{\left(\dfrac{3}{2} q_R \dfrac{\sigma_m}{(1 - f )\sigma_{\star}} \right)}} - 1 = 0 \qquad{(2)}\]

The effective stress \(\sigma_{\star}\) which is solution of this equation is:

\[\sigma_{\star}= \dfrac{q_R \sigma_{vM}}{(1 - f) {\left(q_R - \dfrac{2}{3T}L{\left(q_R D_R f T \exp{{\left(\dfrac{3}{2} q_R T\right)}}\right)}\right)}} \]

where \(L\) is the Lambert function, and \(T = \sigma_m / \sigma_{vM}\) the stress triaxiality. The parameters of the model are:

\[\left\{q_R,D_R \right\} \]

Different nucleation terms are available and are described below. Multiple nucleation terms can be used simultaneously.

**Chu-Needleman strain based**[9]

\[ A_n = \dfrac{f_N}{s_N \sqrt{2\pi}} \exp{{\left(-\dfrac{1}{2} {\left(\dfrac{p - \epsilon_N}{s_N} \right)}^2 \right)} } \qquad{(3)}\]

where \(\{f_N,s_N,\epsilon_N\}\) are material parameters, and \(p\) is the matrix equivalent plastic strain.

**Chu-Needleman stress based**[9]

\[ A_n = \dfrac{f_N}{s_N \sqrt{2\pi}} \exp{{\left(-\dfrac{1}{2} {\left(\dfrac{\sigma_I - \sigma_N}{s_N} \right)}^2 \right)} } \]

where \(\{f_N,s_N,\sigma_N\}\) are material parameters, and \(\sigma_I\) is the maximal positive principal stress.

**Power-Law strain based**

\[ A_n = f_N \left<\dfrac{p}{\epsilon_N} - 1 \right>^m \quad \mathrm{if}\quad \int A_n dp \leq f_{\mathrm{nuc}}^{max} \qquad{(4)}\]

where \(\{f_N,m,\epsilon_N,f_{\mathrm{nuc}}^{max} \}\) are material parameters, \(p\) is the matrix equivalent plastic strain. In the previous expression, \(\left<.\right>\) denotes the Macaulay bracket, i.e. positive part of a number:

\[ \left<x\right> = \max{\left(x,0\right)} \]

**Power-Law stress based**

\[ A_n = f_N \left<\dfrac{\sigma_I}{\sigma_N} - 1 \right>^m \quad \mathrm{if}\quad \int A_n dp \leq f_{\mathrm{nuc}}^{max} \]

where \(\{f_N,m,\sigma_N,f_{\mathrm{nuc}}^{max}\}\) are material parameters, and \(\sigma_I\) is the maximal positive principal stress.

This section is devoted to the presentation of two implicit schemes that may be used to integrate the constitutive equations presented in Section 2.

Assuming a single porous stress criterion whose flow is modified by the porosity evolution, the constitutive equations lead to the following system to be solved for an elastoplastic evolution:

\[ \begin{aligned} \underline{\dot{\varepsilon}}^{\mathrm{el}}+ \underline{\dot{\varepsilon}}^{\mathrm{p}}- \underline{\dot{\varepsilon}}^{\mathrm{to}}&= 0 \\ \sigma_{\star}- R(p) &= 0 \\ \dot{f} - (1 - f){\mathrm{tr}{\left(\underline{\dot{\varepsilon}}^{\mathrm{p}}\right)}} - \sum_j A_n^j \dot{p} &=0 \end{aligned} \qquad{(5)}\]

This system must be closed by adding the relation between the stress \(\underline{\sigma}\) and the elastic strain \(\underline{\varepsilon}^{\mathrm{el}}\). In the report, we assume that the standard Hooke law holds:

\[ \underline{\sigma}= \underline{\underline{\mathbf{D}}}\,\underline{\varepsilon}^{\mathrm{el}} \]

where \(\underline{\underline{\mathbf{D}}}\) is the elastic stiffness tensor.

A semi-implicit method is used to solve these equations according to the variables \(\{\Delta\,\underline{\varepsilon}^{\mathrm{el}}, \Delta\,p, \Delta\,f \}\)

\[ \mathcal{R}{\left(\Delta\,\underline{\varepsilon}^{\mathrm{el}}, \Delta\,p, \Delta\,f\right)}=0 \qquad{(6)}\]

The residual \(\mathcal{R}\) can be decomposed by blocks as follows:

\[ \begin{aligned} \mathcal{R}_{\Delta\,\underline{\varepsilon}^{\mathrm{el}}} = \Delta\,\underline{\varepsilon}^{\mathrm{el}}+ (1 - {\left.f\right|_{t+\theta\,\Delta\,t}}) \Delta\,p \, {\left.{\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}}\right|_{t+\theta\,\Delta\,t}} - \Delta\,\underline{\epsilon}_{to} &= 0 \\ \mathcal{R}_{\Delta\,p} = \sigma_{\star}- R({\left.p\right|_{t+\theta\,\Delta\,t}}) &= 0 \\ \mathcal{R}_{\Delta\,f} = \Delta\,f - (1 - {\left.f\right|_{t+\theta\,\Delta\,t}})^2 \Delta\,p \, {\mathrm{tr}{\left({\left.{\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}}\right|_{t+\theta\,\Delta\,t}}\right)}} - \sum_j {\left.A_n^j\right|_{t+\theta\,\Delta\,t}}\,\Delta\,p &=0 \end{aligned} \qquad{(7)}\]

In the previous equation, \(\theta\) is a numerical parameter (\(\theta\in[0:1]\)).

Solving this system using a Newton-Raphson algorithm requires computing the Jacobian matrix:

\[ J= \begin{pmatrix} \begin{array}{lll} {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} & {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,p}} & {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,f}} \\ {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,p}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} & {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,p}}{\displaystyle \partial \Delta\,p}} & {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,p}}{\displaystyle \partial \Delta\,f}} \\ {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,f}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} & {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,f}}{\displaystyle \partial \Delta\,p}} & {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,f}}{\displaystyle \partial \Delta\,f}} \\ \end{array} \end{pmatrix} \qquad{(8)}\]

where the different terms can be written as [7]:

\[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} &= \underline{\underline{\mathbf{I}}} + \theta \Delta\,p \; \bigg(1 - {\left.f\right|_{t+\theta\,\Delta\,t}}\bigg) {\left.\dfrac{\partial^2 \sigma_{\star}}{\partial \underline{\sigma}^2}\right|_{t+\theta\,\Delta\,t}}\,\colon\,{\left.{\underline{\underline{\mathbf{D}}}}\right|_{t+\theta\,\Delta\,t}} \\ {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,p}} &= \bigg(1 - {\left.f\right|_{t+\theta\,\Delta\,t}}\bigg) {\left.{\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}}\right|_{t+\theta\,\Delta\,t}} \\ {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,f}} &= \theta \Delta\,p \; \bigg( \bigg(1 - {\left.f\right|_{t+\theta\,\Delta\,t}}\bigg)\, {\left.\dfrac{\partial^2 \sigma_{\star}}{\partial \underline{\sigma}\partial f}\right|_{t+\theta\,\Delta\,t}} - {\left.{\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}}\right|_{t+\theta\,\Delta\,t}} \bigg) \\ {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,p}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} &= \theta \; {\left.{\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}}\right|_{t+\theta\,\Delta\,t}}\,\colon\,{\left.{\underline{\underline{\mathbf{D}}}}\right|_{t+\theta\,\Delta\,t}} \\ {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,p}}{\displaystyle \partial \Delta\,p}} &= - \theta \; \dfrac{d R(p)}{d p} \\ {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,p}}{\displaystyle \partial \Delta\,f}} &= \theta \; {\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial f}} |_{t+\theta \Delta\,t} \\ {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,f}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} &= - \theta \; \Delta\,p \; \bigg(1 - {\left.f\right|_{t+\theta\,\Delta\,t}}\bigg)^2 \; \bigg({\left.\dfrac{\partial^2 \sigma_{\star}}{\partial \underline{\sigma}^2}\right|_{t+\theta\,\Delta\,t}}\,\colon\,{\left.{\underline{\underline{\mathbf{D}}}}\right|_{t+\theta\,\Delta\,t}}\bigg)\,\colon\,\underline{I} - \theta \Delta\,p \sum_j {\left.{\displaystyle \frac{\displaystyle \partial A_n^j}{\displaystyle \partial \underline{\sigma}}}\right|_{t+\theta\,\Delta\,t}}\,\colon\,{\left.{\underline{\underline{\mathbf{D}}}}\right|_{t+\theta\,\Delta\,t}} \\ {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,f}}{\displaystyle \partial \Delta\,p}} &= - \bigg(1 - {\left.f\right|_{t+\theta\,\Delta\,t}}\bigg)^2 \; {\left.{\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}}\right|_{t+\theta\,\Delta\,t}}\,\colon\,\underline{I} - \theta \Delta\,p \sum_j {\left.{\displaystyle \frac{\displaystyle \partial A_n^j}{\displaystyle \partial p}}\right|_{t+\theta\,\Delta\,t}} - \sum_j {\left.A_n^j\right|_{t+\theta\,\Delta\,t}} \\ {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,f}}{\displaystyle \partial \Delta\,f}} &= 1 - \theta \; \Delta\,p \; \bigg(1 - {\left.f\right|_{t+\theta\,\Delta\,t}}\bigg) \bigg(-2 \; {\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}} + \bigg(1 - {\left.f\right|_{t+\theta\,\Delta\,t}}\bigg) \; \dfrac{\partial^2 \sigma_{\star}}{\partial \underline{\sigma}\partial f} \bigg) \underline{I} - \theta \Delta\,p \sum_j {\left.{\displaystyle \frac{\displaystyle \partial A_n^j}{\displaystyle \partial f}}\right|_{t+\theta\,\Delta\,t}} \\ \end{aligned} \]

Therefore, the definition of a porous model requires to give the following expressions for the stress criterion \(\sigma_{\star}\):

\[\left\{\sigma_{\star}, {\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}}, \dfrac{\partial \sigma_{\star}}{\partial f}, \dfrac{\partial^2 \sigma_{\star}}{\partial \underline{\sigma}^2}, \dfrac{\partial^2 \sigma_{\star}}{\partial \underline{\sigma}\partial f} \right\} \]

as well as the following expressions for each nucleation mechanism:

\[ \left\{A_n^i, {\displaystyle \frac{\displaystyle \partial A_n^i}{\displaystyle \partial \underline{\sigma}}}, {\displaystyle \frac{\displaystyle \partial A_n^i}{\displaystyle \partial p}} , {\displaystyle \frac{\displaystyle \partial A_n^i}{\displaystyle \partial f}} \right\} \]

These expressions are detailed for the stress criteria and nucleation formulas currently implemented into the brick in Appendixes 8.2 and 8.3.

Adding the elastic contribution to the porosity growth is trivially performed by substracting the following term to the implicit equation associated with the porosity evolution:

\[ \mathcal{R}_{\Delta\,f} \mathrel{-}= {\left(1 - {\left.f\right|_{t+\theta\,\Delta\,t}}\right)}\,{\mathrm{tr}{\left(\Delta\underline{\varepsilon}^{\mathrm{el}}\right)}} \]

The following term must be added to the block \({\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,f}}{\displaystyle \partial \Delta\,f}}\):

\[ {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,f}}{\displaystyle \partial \Delta\,f}} \mathrel{+}= \theta\,{\mathrm{tr}{\left(\Delta\underline{\varepsilon}^{\mathrm{el}}\right)}} \]

And the following term must be added to the block \({\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,f}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}\):

\[ {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,f}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} \mathrel{+}= -{\left(1 - {\left.f\right|_{t+\theta\,\Delta\,t}}\right)}\,\underline{I} \]

In the absence of nucleation models and without elastic contribution, the evolution of the porosity is given by the following equation:

\[ \dot{f}={\left(1-f\right)}^{n_{f}}\,\dot{p}\,{\mathrm{tr}{\left(\underline{n}\right)}} \qquad{(9)}\]

where \(n_{f}\) is equal to:

- \(2\) if the porosity evolution affects the (visco-)plastic flow
- \(1\) if this is not the case.

Assuming \(n_{f}=1\), Equation (9) can be integrated by parts:

\[ \int_{{\left.f\right|_{t}}}^{{\left.f\right|_{t+\Delta\,t}}}\dfrac{\mathrm{d}\,u}{1-u}=\log{\left(\dfrac{1-{\left.f\right|_{t}}}{1-{\left.f\right|_{t+\Delta\,t}}}\right)}=\int_{t}^{t+\Delta\,t}\dot{p}\,{\mathrm{tr}{\left(\underline{n}\right)}}\,\mathrm{d}\,t\approx\Delta\,p\,{\mathrm{tr}{\left({\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\right)}} \]

which leads to the following expression of the increment of the porosity:

\[ \Delta\,f={\left(1-{\left.f\right|_{t}}\right)}\,{\left(1-\exp{\left(-\Delta\,p\,{\mathrm{tr}{\left({\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\right)}}\right)}\right)} \]

Taking into account the elastic predictionWhen taking into account the elastic prediction, the increment of the porosity can be known prior any integration:

\[ \Delta\,f={\left(1-{\left.f\right|_{t}}\right)}\,{\left(1-\exp{\left({\mathrm{tr}{\left(-\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)}}\right)}\right)} \qquad{(10)}\]

When applicable, Equation (10), could be very interesting because one can predict the material failure before the material integration. Moreover, removing the evolution of the porosity equation greatly simplifies the resolution of the implicit system.

However, one drawback of this relation is that it introduces the porosity as an implicit function of the increment of the total strain which:

- is incompatible with the assumptions of the
`StandardElastoViscoPlasticity`

brick. It thus won’t be treated further in this report;- complexifies the computation of the consistent tangent operator. In particular, Equation (18), detailed below, is no more valid. However, recent developments of
`MFront`

could simplify the computation of the tangent operator in this case, see [10] for details.Another drawback is that the user would not be able to add its own contribution to the porosity evolution as described in Section 4.2.0.1.

Assuming \(n_{f}=2\), Equation (9) can be integrated by parts:

\[ \int_{{\left.f\right|_{t}}}^{{\left.f\right|_{t+\Delta\,t}}}\dfrac{\mathrm{d}\,u}{{\left(1-u\right)}^{2}}=\dfrac{1}{1-{\left.f\right|_{t}}}-\dfrac{1}{1-{\left.f\right|_{t+\Delta\,t}}}=\int_{t}^{t+\Delta\,t}\dot{p}\,{\mathrm{tr}{\left(\underline{n}\right)}}\,\mathrm{d}\,t\approx\Delta\,p\,{\mathrm{tr}{\left({\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\right)}} \]

which leads to the following expression of the increment of the porosity:

\[ \Delta\,f={\left(1-{\left.f\right|_{t}}\right)}^{2}\,\dfrac{ \Delta\,p\,{\mathrm{tr}{\left({\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\right)}} }{ 1+{\left(1-{\left.f\right|_{t}}\right)}\,\Delta\,p\,{\mathrm{tr}{\left({\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\right)}} } \]

Most nucleation models require to set an upper bound \(f_{\mathrm{nuc}}^{j,\max}\) on the maximum nucleated porosity predicted by those models. In this case, the effective nucleated porosity increment is defined as:

\[ \Delta\,f_{n}^{j} = \min{\left({\left.A_n^j\right|_{t+\theta\,\Delta\,t}} \Delta\,p,\,f_{\mathrm{nuc}}^{j,\max}-{\left.f_{n}^{j}\right|_{t}}\right)} \]

This treatment requires that the value of the nucleated porosity at
the beginning of the time step is known. Hence it is automatically
stored in an auxiliary state variable whatever the value of the
`save_porosity_increase`

option passed to the nucleation
model.

A strain-based nucleation model only depends on the equivalent plastic strain as follows:

\[ {\displaystyle \frac{\displaystyle \mathrm{d}f_{n}^{j}}{\displaystyle \mathrm{d}t}}=A^{j}_{n}{\left(p\right)}\,{\displaystyle \frac{\displaystyle \mathrm{d}p}{\displaystyle \mathrm{d}t}} \qquad{(11)}\]

See the strain version of the Chu-Needleman nucleation model (Equation (3)) and the strain version of the power-law nucleation model (Equation (4)).

In Equation (7), the contribution of this nucleation model to the implicit equation was set equal to:

\[ \Delta\,f_{n}^{j}\approx A^{j}_{n}{\left({\left.p\right|_{t+\theta\,\Delta\,t}}\right)}\,\Delta\,p \]

However, this is only an approximation.

A more precise estimation can be obtained by integrating Equation (11) by parts, which leads to the following expression of the porosity increment:

\[ \Delta\,f_{n}^{j} = F^{j}_{n}{\left({\left.p\right|_{t+\Delta\,t}}\right)}-F^{j}_{n}{\left({\left.p\right|_{t}}\right)} \qquad{(12)}\]

where \(F^{j}_{n}\) is one primitive of \(A\), i.e. \(F^{j}_{n}{\left(p\right)}=\displaystyle \int A^{j}_{n}{\left(u\right)}\,\mathrm{d}\,u\). Such a treatment of the nucleated porosity was only considered in [11].

When the function \(A\) also depends on some external state variables, (temperature, fluence), we use the value of those external state variables at the middle of the time step to evaluate the increment of the porosity using Equation (12).

The Newton-Raphson algorithm may fail at finding the solution of the non-linear system of equations, especially when the time step leads to strong increase of the porosity.

Moreover, the porosity evolution leads to specific issues:

- Phenomenological laws describing porosity nucleation may have thresholds to be activated or inhibited.
- Material failure usually happens when the porosity reaches a critical value corresponding to the collapse of the yield surface. Detection of material failure is awkward when the equations governing the porosity evolution are included in the implicit system to be solved. For example, the estimates of the solution may exceed the critical value during the iterations while the solution of the implicit scheme may be below this critical value.

As the main convergence issues of the standard implicit scheme are associated with the porosity evolution, we propose in this section a staggered approach where the porosity evolution and the evolution of the other state variables (i.e. the elastic strain \(\underline{\varepsilon}^{\mathrm{el}}\) and the equivalent plastic strain \(p\)) are decoupled. A fixed point algorithm is used to ensure that the increments \(\left\{\Delta\,\underline{\varepsilon}^{\mathrm{el}}, \Delta\,p, \Delta\,f\right\}\) satisfies the Implicit System (7).

Let \(i\) be the current step of the fixed point algorithm. The estimates of the increments of the elastic strain and the equivalent plastic strain at the next iterations, denoted respectively \(\Delta \underline{\varepsilon}^{\mathrm{el}}_{(i+1)}\) and \(\Delta\,p_{(i+1)}\) satisfies a reduced implicit system:

\[ \begin{aligned} \mathcal{R}_{\Delta\,\underline{\varepsilon}^{\mathrm{el}}} = \Delta\,\underline{\varepsilon}^{\mathrm{el}}_{(i+1)} + {\left(1 - f_{(i)}\right)}\,\Delta\,p_{(i+1)}\,{\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}}{\left(\underline{\varepsilon}^{\mathrm{el}}_{(i+1)},f_{(i)}\right)} - \Delta\,\underline{\epsilon}^{to} &= 0 \\ \mathcal{R}_{\Delta\,p} = \sigma_{\star}{\left(\underline{\varepsilon}^{\mathrm{el}}_{(i+1)},f_{(i)}\right)} - R(p_{(i+1)}) &= 0 \\ \end{aligned} \qquad{(13)}\]

where:

- \(\underline{\varepsilon}^{\mathrm{el}}_{(i+1)}={\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t}}+\theta\,\Delta\,\underline{\varepsilon}^{\mathrm{el}}_{(i+1)}\)
- \(p_{(i+1)}={\left.p\right|_{t}}+\theta\,\Delta\,p_{(i+1)}\)
- \(f_{(i)}={\left.f\right|_{t}}+\theta\,\Delta\,f_{(i)}\)

The Implicit System (13) is solved by a standard Newton-Raphson method with the reduced jacobian matrix:

\[ \left[ \begin{array}{lll} {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} & {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,p}} \\ {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,p}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} & {\displaystyle \frac{\displaystyle \partial \mathcal{R}_{\Delta\,p}}{\displaystyle \partial \Delta\,p}} \\ \end{array} \right] \qquad{(14)}\]

NoteIn practice, for reasons detailed in Section 3.2.10, we still solve the full Implicit System (7), but the implicit equation is modified as follows:

\[ \mathcal{R}_{\Delta\,f} = \Delta\,f - \Delta\,f_{(i)} \]

Letting the constraints on the porosity evolution aside, the porosity is determined iteratively as follows:

\[ \Delta\,f_{(i+1)}^{(uc)} = {\left(1 - f_{(i)}\right)}^2 \Delta\,p_{(i+1)} \, {\mathrm{tr}{\left({\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}}{\left(\underline{\varepsilon}^{\mathrm{el}}_{(i+1)},f_{(i)}\right)}\right)}} + \sum_{j} A_n^{j}{\left(\underline{\varepsilon}^{\mathrm{el}}_{(i+1)}, p_{(i+1)}, f_{(i)}\right)} \Delta\,p_{(i+1)} \qquad{(15)}\]

Here the \(\mbox{}^{(uc)}\) superscript means “uncorrected”. Corrections to \(\Delta\,f_{(i+1)}^{(uc)}\) will be introduced in the next paragraphs to deal with thresholds in nucleation laws and detection of material failure.

Let:

- \(\Delta \left.f_{n}^{j, (uc)}\right|_{(i+1)}=A_n^{j}{\left(\underline{\varepsilon}^{\mathrm{el}}_{(i+1)}, p_{(i+1)}, f_{(i)}\right)} \Delta\,p_{(i+1)}\) be the uncorrected contribution of the \(j^{\mathrm{th}}\) nucleation law to the porosition evolution.
- \({\left.f_{n}^{j}\right|_{t}}\) be the value of the nucleated porosity due to this nucleation law at the beginning of the time step.
- \(f_{n}^{j, \max{}}\) be the threshold of this nucleation law, i.e. the maximal value of the nucleated porosity due to this law.

\(f_{n}^{j, \max{}}-{\left.f_{n}^{j}\right|_{t}}\) is the maximum allowed increment of the porosity for this nucleation law.

We then define the corrected contribution \(\Delta \left.f_{n}^{j}\right|_{(i+1)}\) as follows:

\[ \Delta \left.f_{n}^{j}\right|_{(i+1)} = \min{\left(\Delta \left.f_{n}^{j, (uc)}\right|_{(i+1)}, f_{n}^{j,\max{}}-{\left.f_{n}^{j}\right|_{t}}\right)} \]

NoteTo avoid the introduction of another notation, we still used \(\Delta\,f_{(i+1)}^{(uc)}\) even though corrections related to thresholds in nucleations laws may have been taken into account at this stage.

If \({\left.f\right|_{t}}+\Delta\,f_{(i+1)}^{(uc)}\) is greater than the critical porosity denoted \(\alpha_{1}\,f_{r}\), a dichotomic approach is used:

\[ \Delta\,f_{(i+1)} = \dfrac{1}{2}{\left(f+\Delta\,f_{(i)} + \alpha_{1}\, f_r\right)} - {\left.f\right|_{t}}; \]

This approach allows the porosity to approach the critical porosity smoothly.

The \(\alpha_{1}\) factor is a user defined parameter, chosen equal by default to \(98.5\,\%\). The reason for not allowing the porosity be closer to \(f_{r}\) is than the Implicit System (13) becomes very difficult to solve as the yield surfaces collapses.

The iterations are stopped when the porosity becomes stationnary, i.e. when:

\[ \left|\Delta\,f_{(i+1)}^{(uc)}-\Delta\,f_{(i)}\right|<\varepsilon_{f} \qquad{(16)}\]

where the value of \(\varepsilon_{f}\) is a user defined (by default, a stringent value of \(10^{-10}\) is used).

The Iterative Process (15) can be accelerated using the well known Aitken transformation.

The usage of the Aitken acceleration is enabled by default.

To avoid spurious oscillations, a relaxation procedure is used:

\[ \Delta\,f_{(i+1)}=\omega\,\Delta\,f_{(i)}+(1 - \omega)\,\Delta\,f_{(i+1)}^{(uc)} \]

where the relaxation coefficient is chosen equal to \(0.5\).

Once the Stopping Condition (16) is satisfied, material failure is detected when the final porosity is above a given threshold \(\alpha_{2}\,f_{r}\), where \(\alpha_{2}\) is a user defined constant chosen equal to \(98\,\%\) by default.

When the material failure is detected, the value of the auxiliary
state variable `broken`

is set to one.

From a performance point of view, this staggered approach seems *a
priori* less efficient than the standard implicit scheme because the
reduced integration scheme will be solved once per iteration of the
fixed point algorithm.

The staggered approach is however thought to be more robust, although we do not have performed an in-depth comparison of both monolithic and staggered approaches yet.

It shall also be noted that the number of iterations of the reduced integration implicit system is not proportional to the number of iterations of the fixed point algorithm. Indeed, the solution of the reduced implicit system obtained at the previous iteration of the fixed point algorithm can be used as the starting point of the current resolution of the reduced implicit system. Numerical experiments shows that this strategy is very effective and that close to convergence of the fixed point algorithm the number of iterations required to solve the reduced implicit system is low, i.e. only 2 or 3 iterations are required.

The variable `staggered_scheme_iteration_counter`

holds
the number of iterations of the staggered scheme. One can save this
variable in a auxiliary state variable as follows:

```
//! number of iterations of the staggered scheme at convergence
@AuxiliaryStateVariable real niter;
.setEntryName("StaggeredSchemeIterationCounter");
niter
....
@UpdateAuxiliaryStateVariables{
= staggered_scheme_iteration_counter;
niter }
```

The total number of evaluations of the implicit scheme is a bit trickier to save. It can be done as follows:

```
//! number of iterations of the implicit scheme
@AuxiliaryStateVariable real neval;
.setEntryName("NumberOfEvaluationOfTheImplicitScheme");
neval
...
@InitLocalVariables{
...
= 0;
neval }
...
@Integrator{
++neval;
}
```

NoteIf local sub-stepping is allowed, those two numbers will only reflect the number of iterations and number of evaluations of the last step.

`@AdditionalConvergenceChecks`

code blockAs depicted on Figure 1, each iteration of the implicit algorithm are decomposed in three steps:

- Updating the stress in the
`@ComputeStress`

code block. This code block is automatically handled by the brick. - Building the residual \(\mathcal{R}\) and the jacobian matrix \(J\) in the
`@Integrator`

code block. - Additional convergence checks in the
`@AdditionalConvergenceChecks`

code block. Additional here means that`MFront`

has already checked the convergence of the resolution of implicit scheme using a built-in criterion. This code block may be used to implement a status algorithm which allows activating and/or desactivating some (visco-)plastic flows (see [12]). Here, it is used to implement the staggered approach described in this section.

As described in [10], the consistent tangent operator may be deduced for the jacobian of Implicit System (7) as follows:

\[ {\displaystyle \frac{\displaystyle \mathrm{d}\underline{\sigma}}{\displaystyle \mathrm{d}\Delta\,\underline{\varepsilon}^{\mathrm{to}}}}= {\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}}\,\cdot\,{\displaystyle \frac{\displaystyle \partial \underline{\varepsilon}^{\mathrm{el}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} \qquad{(17)}\]

where \({\displaystyle \frac{\displaystyle \partial \underline{\varepsilon}^{\mathrm{el}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\) can be determined by solving the following systems of linear equations:

\[ J\,\cdot\,{\displaystyle \frac{\displaystyle \partial \underline{\varepsilon}^{\mathrm{el}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}=-{\displaystyle \frac{\displaystyle \partial \mathcal{R}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} \]

where \({\displaystyle \frac{\displaystyle \partial \mathcal{R}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\) have this simple matrix form (in \(3D\)):

\[ {\displaystyle \frac{\displaystyle \partial \mathcal{R}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}=- \begin{pmatrix} 1& 0 & 0 & 0& 0 &0 \\ 0& 1& 0 & 0& 0 &0 \\ 0& 0 & 1 & 0& 0 &0 \\ 0& 0 & 0 & 1& 0 &0 \\ 0& 0 & 0 & 0& 1 &0 \\ 0& 0 & 0 & 0& 0 &1\\ 0& 0 & 0 & 0& 0 &0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ 0& 0 & 0 & 0& 0 &0 \\ \end{pmatrix} \qquad{(18)}\]

It would cumbersome to compute the jacobian of the Reduced System (13) during the resolution and then to compute the jacobian of the Full System (7) after convergence for the computation of the consistent tangent operator:

- it would require a significant amount of boiler plate code to define the jacobian matrix for the full implicit system and appropriate views of its blocks (see Equation (8));
- it would require a partial duplication of the code computing the blocks of the jacobian of the reduced implicit system (See Equations (8) and (14)).

We thus rely on a simple trick. We define a flag which is set to true during the iterations of the fixed point algorithm. If this flag is true, the implicit equation of the Full Implicit System (7) is replaced by:

\[ \mathcal{R}_{\Delta\,f} = \Delta\,f - \Delta\,f_{(i)} \]

where \(\Delta\,f_{(i)}\) is the current estimation given by the fixed point algorithm.

This new equation is trivially solved so that the solutions of this modified system do verify the Reduced Implicit System (13). Once the fixed point algorithm has converged, one may simply set the flag to false and build the jacobian matrix (8) of the Full Implicit System (7) by calling the assembly of the implicit scheme once and finally apply Equation (17).

A few things may be noted:

- the overhead implied by this trick seems
*a priori*limited; - one can easily recover the standard implicit scheme;
- one can easily choose to treat the porosity explicitly.

`StressCriterion`

and `InelasticFlow`

interfacesThe `InelasticFlow`

interface is used by the
`StandardElastoViscoPlasticity`

brick to handle inelastic
flows. It uses the `StressCriterion`

interface to describe it
stress criterion and, for non-associated flows, its flow criterion.

For the brick to properly treat inelastic flows coupled with
porosity, a new method called `isCoupledWithPorosity-`

`Evolution`

has been added to those interface. This method
returns a boolean value and does not take any argument.

For inelastic flows, the default implementation of the
`isCoupledWithPorosityEvolution`

, declared in the
`In-`

`elasticFlowBase`

class, returns
`true`

if its stress criterion or if its flow criterion (if
any) declares to be coupled with the porosity evolution.

The `isCoupledWithPorosityEvolution`

of a stress criterion
must return `true`

if and only if the expression of the
stress criterion explicitly depends on the porosity.

NoteImplicitly, the brick expects that a stress criterion coupled with porosity is

notdeviatoric, i.e., that this stress criterion contributes to the porosity growth. The`isNormalDeviatoric`

method of such a stress criterion must return`true`

. This is checked in the`InelasticFlowBase`

class.

As discussed in the introduction, a stress criterion can lead to a
correction of the flow rule by a factor \(1-f\). The class
`StressCriterion`

now has a method called
`getPorosityEffectOnEquivalentPlasticStrain`

which return one
of the two following values:

`NO_POROSITY_EFFECT_ON_EQUIVALENT_PLASTIC_STRAIN`

`STANDARD_POROSITY_CORRECTION_ON_EQUIVALENT_PLASTIC_STRAIN`

The evolution of the porosity is taken into account:

- if a nucleation model is declared;
- if one stress inelastic flow declares being coupled with the porosity evolution, i.e. that its stress criterion or its flow criterion is coupled with the porosity evolution;
- if the
`porosity_evolution`

section is explicitly declared, even if empty.

For example, this declaration will automatically lead to a porosity growth:

```
@Brick StandardElastoViscoPlasticity{
: "Hooke" {
stress_potential : 150.e3,
young_modulus : 0.3
poisson_ratio },
: "Plastic" {
inelastic_flow : "Gurson1975" {},
criterion : "Linear" {R0 : "0"}
isotropic_hardening }
};
```

NoteIt is important to note that, by default (if none of the three above conditions are met), no contribution to the porosity growth will be taken into account, even if the flow direction has an hydrostatic component. For example, the following declaration will not lead to any porosity growth:

`@Brick StandardElastoViscoPlasticity{ : "Hooke" { stress_potential : 150.e3, young_modulus : 0.3 poisson_ratio }, : "Plastic" { inelastic_flow : "MohrCoulomb" { criterion : 3.e1, // cohesion c : 0.523598775598299, // friction angle angle phi : 0.506145483078356, // transition angle (Abbo and Sloan) lodeT : 1e1 // tension cuff-off parameter a }, : "Linear" {R0 : "0"} isotropic_hardening } };`

However, this declaration will lead to a porosity growth:

`@Brick StandardElastoViscoPlasticity{ : "Hooke" { stress_potential : 150.e3, young_modulus : 0.3 poisson_ratio }, : "Plastic" { inelastic_flow : "MohrCoulomb" { criterion : 3.e1, // cohesion c : 0.523598775598299, // friction angle angle phi : 0.506145483078356, // transition angle (Abbo and Sloan) lodeT : 1e1 // tension cuff-off parameter a }, : "Linear" {R0 : "0"} isotropic_hardening } : {} porosity_evolution};`

The `porosity_evolution`

section has a significant
importance in this document. The various options that can be declared in
this section are described in Appendix 8.6.

If the evolution of the porosity is taken into account, a state
variable called `f`

with the glossary name
`Porosity`

is automatically declared if no state variable
with this glossary name has been declared. In the latter case, this
variable is used (and updated) by the brick.

This choice allows the user to easily add new nucleation model or to
add a new inelastic flow coupled with porosity by adding the appropriate
state variables and adding the associated implicit equations in the
`@Integrator`

code block.

For example, assuming that the staggered approach described in Section 3.2 is used, one may add the contribution of the solid swelling as follows:

```
@ExternalStateVariable real ss;
.setGlossaryName("SolidSwelling");
ss
@Brick StandardElastoViscoPlasticity{
: "Hooke" {
stress_potential : 150.e3,
young_modulus : 0.3
poisson_ratio },
: "Plastic" {
inelastic_flow : "Gurson1975" {},
criterion : "Linear" {R0 : "0"}
isotropic_hardening }
};
@Integrator{
// adding the contribution of the solid swelling
// the porosity increase in the implicit equation describing
// the porosity.
// This shall only been done after the convergence of the
// staggered scheme has converged.
if(compute_standard_system_of_implicit_equations){
-= ;
ff }
}
@AdditionalConvergenceChecks {
// adding the contribution of the solid swelling
// the porosity increase to the next extimate of the porosity by
// the staggered scheme.
// This shall only been done after the convergence of the
// staggered scheme has converged.
if (converged && (!compute_standard_system_of_implicit_equations)) {
+= ...
next_estimate_of_the_porosity_increment }
}
```

The previous code uses variables that are automatically defined by
the `StandardElastoViscoPlasticity`

brick when the staggered
approach is used:

`compute_standard_system_of_implicit_equations`

is a boolean value which is false until the staggered scheme converges;`next_estimate_of_the_porosity_increment`

is a variable only accessible in the`@AdditionalCon-`

`vergenceChecks`

block which must hold the next estimate of the porosity increment. This variable is automatically initialized to zero.

Other variables automatically defined by the
`StandardElastoViscoPlasticity`

brick may be helpful to
couple constitutive equations with user defined ones:

in the

`@InitLocalVariables`

code block, the`upper_bound_of_the_porosity`

variable may be updated by the user to change the value of the upper bound of the porosity as follows:`= min(upper_bound_of_the_porosity, upper_bound_of_the_porosity .... // user_defined_bound );`

in the

`@AdditionalConvergenceChecks`

code block, the`fixed_point_converged`

boolean variable might be used add other convergence criteria to the staggered algorithm, as follows:`if (converged && (!compute_standard_system_of_implicit_equations)) { = .... fixed_point_converged }`

By default, the value returned by the
`getPorosityEffectOnEquivalentPlasticStrain`

method of the
stress criterion is used to decipher if the porosity affects definition
of the equivalent plastic strain and the flow rule.

However, this may be changed by specifying the
`porosity_effect_on_equivalent_plastic_strain`

option in the
definition of the inelastic flow. Valid strings are
`StandardPorosityEffect`

(or equivalently ”
`standard_-`

`porosity_effect`

) or
`None`

(or equivalently `none`

or
`false`

).

The porosity evolves either due to nucleation or growth. For post-processing reasons, it may be worth to distinguish those contributions in dedicated auxiliary state variables. However, for performance reasons (mostly because additional auxiliary state variables implied a memory overhead), this shall not be the default choice.

The brick allows the user to specify if those additional variables
shall be defined for all mechanisms in the `porosity-`

`_evolution`

section as follows:

```
: {
porosity_evolution: true
save_individual_porosity_increase}
```

However, this behaviour can be overloaded in every inelastic flow
section and every nucleation model by setting the
`save_porosity_increase`

boolean variable.

For example, one may use the following declaration:

```
: "Plasticity" {
inelastic_flow : "Gurson1975",
criterion : "Voce" {R0 : 200, Rinf : 100, b : 20},
isotropic_hardening : true
save_porosity_increase }
```

NoteInelastic flows which declare that the flow is deviatoric discard those options and never declare an auxiliary state variable associated with the porosity growth.

By default, the elastic contribution in the porosity growth is not taken into account.

This could be changed by setting the
`elastic_contribution`

options to true as follows:

```
: {
porosity_evolution: true
elastic_contribution}
```

The stress criteria implemented are verified by comparing the results
obtained through `MTest`

simulations to reference solutions.
Axisymmetric loading conditions are considered under small strains
settings:

\[ \underline{\sigma} = \sigma \begin{pmatrix} 1 & 0 & 0 \\ 0 & A & 0\\ 0 & 0 & A \end{pmatrix} \]

A typical `MTest`

file used for these simulations is shown
below:

```
@ModellingHypothesis 'Tridimensional';
@Behaviour<Generic> 'src/libBehaviour.so' 'GTN';
@InternalStateVariable 'Porosity' 1e-3;
@ExternalStateVariable 'Temperature' 293.15 ;
@Real 'A' '0.4';
@NonLinearConstraint<Stress> 'SYY - A*SXX';
@NonLinearConstraint<Stress> 'SZZ - A*SXX';
@ImposedStrain 'EXX' {0 : 0, 1 : 0.5};
@Times{0, 1. in 1000};
```

Without hardening, the differential equations governing the
evolutions of internal state variables are, for the GTN model: \[
\begin{aligned}
\frac{d \varepsilon_{yy}^p}{d \varepsilon_{xx}^p} &= \frac{q_1 q_2
f_{\star}\sinh{{\left(\dfrac{1}{2} q_2 (1+2\alpha)
\dfrac{\sigma}{\sigma_0}\right)}} + \dfrac{\sigma}{\sigma_0}(\alpha -
1)}{q_1 q_2 f_{\star}\sinh{{\left(\dfrac{1}{2} q_2 (1+2\alpha)
\dfrac{\sigma}{\sigma_0}\right)}} + 2\dfrac{\sigma}{\sigma_0}(1 - \alpha
)} \ \ \ \ \ \mathrm{with} \ \ \ \ \ S(\sigma,f) = 0\\
\frac{df}{d \varepsilon_{xx}^p} &= (1 - f)\left(1+ 2 \frac{d
\varepsilon_{yy}^p}{d \varepsilon_{xx}^p} \right)
\end{aligned}
\] The same equations also hold for the Gurson model with \((q_1,q_2,f_{\star}) = (1,1,f)\). For the
Rousselier model, the equations are: \[
\begin{aligned}
\frac{d \varepsilon_{yy}^p}{d \varepsilon_{xx}^p} &= \frac{D_R
q_R f\exp{{\left(\dfrac{1}{2} \dfrac{q_R}{1-f} (1+2\alpha)
\dfrac{\sigma}{\sigma_0}\right)}} - 1.5}{D_R
q_R f\exp{{\left(\dfrac{1}{2} \dfrac{q_R}{1-f} (1+2\alpha)
\dfrac{\sigma}{\sigma_0}\right)}} + 3} \ \ \ \ \ \mathrm{with} \ \ \ \ \
S(\sigma,f) = 0\\
\frac{df}{d \varepsilon_{xx}^p} &= (1 - f)\left(1+ 2 \frac{d
\varepsilon_{yy}^p}{d \varepsilon_{xx}^p} \right)
\end{aligned}
\] These differential equations are solved with respect to \(\varepsilon_{xx}^p\) using the
`ode45`

function of the `Matlab`

software, and the
stress magnitude \(\sigma\) is computed
by solving the stress criteria. Elastic strains are finally added to the
plastic strains using Hooke’s law to obtain total strains. These
solutions, that can be obtained for arbitrary precision, are the
reference solutions to which the results from `MTest`

simulations are compared.

A perfect agreement is obtained between the reference solutions and
the `MTest`

results, as depicted on Figure 2, validating the
implementation of the equations of these models. Finally, the jacobians
of all models have been verified with respect to the numerical
jacobian.

In this scope, the `StandardElastoViscoPlasticity`

brick is tested using the finite elements software `Cast3M`

. The
Gurson-Tvergaard-Needleman stress criterion is used to model the ductile
failure of Notched Tensile (NT) Samples and Simple Tensile (ST) ones.
The \(2D\) meshes are shown on Figure 3
.

The material’s behavior is adapted to aluminum alloys from which the vessel of Fast Neutron Research Reactors are fabricated. The element size is chosen to be \(100\,\mu m\) based on metallurgical studies that take into account the importance of the nucleation phase during damage. The nucleation is modeled by a Stress-based law that is implemented as follows [13] :

\[ A_n = f_N \left<\frac{\sigma_I}{\sigma_N} - 1 \right>^m \ \ \ \ \ \mathrm{if}\ \ p \geq p_N, \ \ \mathrm{and}\ \ \int A_n dp \leq f_{\mathrm{nuc}}^{max} \]

where \(\{f_N, p_N, m,\sigma_N,f_{\mathrm{nuc}}^{max}\}\) are material parameters, and \(\sigma_I\) is the maximum positive principal stress.

Material’s Behavior`@Brick StandardElastoViscoPlasticity{ : "Hooke" {young_modulus : 70e3, poisson_ratio : 0.3}, stress_potential : "Plastic" { inelastic_flow : "GursonTvergaardNeedleman1982" { criterion : 0.04, f_c : 0.056, f_r : 2., q_1 : 1., q_2 : 4. q_3 }, : "Linear" {R0 : 274}, isotropic_hardening : "Voce" {R0 : 0, Rinf : 85, b : 17}, isotropic_hardening : "Voce" {R0 : 0, Rinf : 17, b : 262}, isotropic_hardening : "Power-based" { nucleation_model : 2.92, An0 : 500, s1_n : 0.0321, p_n : 700, s1_max : 0.0215, fg_max : 2 ng } }`

The numerical response of NT specimens (notch radii are noted as : NT10, NT4, and NT2) is shown on Figure 4 and compared to the experimental results.

One can see that the numerical convergence is rather good at the
crack initiation phase where the reaction force decreases suddenly. Some
aspects could be improved regarding the plastic flow and the damage
parameters to smoothly follow the crack propagation modelled by the
sudden force drop. Such improvements depends on the finite element
solver used and must not be confused with the robustness of the `StandardElastoViscoPlasticity`

brick.

Moreover, the evolution of the porosity was monitored during the simulation to observe the failure of elements that reached the critical porosity value. This could be seen on Figure 5 where the notch is pictured before the testing (left image), during the testing (middle image), and at failure (right image).

The GTN model is used to perform a simulation of a Charpy test. The
parameters of the model are taken from [14] and correspond to the mechanical
behavior of a Reactor Pressure Vessel steel at \(300^{\circ}\)C. The corresponding
`MFront`

implementation of the constitutive equations is:

Material’s Behavior`@DSL Implicit; @Behaviour GTN; @UMATUseTimeSubStepping true; @UMATMaximumSubStepping 100; @ModellingHypotheses{".+"}; @StrainMeasure Hencky; @Algorithm NewtonRaphson; @Epsilon 1.e-12; @Theta 1; @Brick StandardElastoViscoPlasticity{ : "Hooke" {young_modulus : 210e3, poisson_ratio : 0.3}, stress_potential : "Norton" { inelastic_flow : "GursonTvergaardNeedleman1982" { criterion : 0.001, f_c : 0.33, f_r : 1.5, q_1 : 1, q_2 : 2.25 q_3 }, : 1.15, n : 0.185, K : 1, A : "Linear" {R0 : 421.5}, isotropic_hardening : "Voce" {R0 : 0, Rinf : 125, b : 59.6}, isotropic_hardening : "Voce" {R0 : 0, Rinf : 472, b : 1.72} isotropic_hardening }, : { porosity_evolution : "PowerLaw (strain)" { nucleation_model : 0.038, en : 0.5, m : 0., fmax : 0.023} fn } }`

Simulations are performed with `Cast3M`

finite element
code, using second-order reduced-integration elements (user-modified
`CU20`

elements). Contact is modelled between the specimen
and the striker and anvils accounting for friction (\(\mu = 0.1\)). Only one fourth of the system
is modelled due to symmetries. The initial porosity is equal to 0.0175%.
The displacement of the striker is imposed with a velocity of \(5\mathrm{m.s^{-1}}\). The time-step is
adjusted with a target on the maximal increase of local porosity of 0.1%
(see Appendix F). Element deletion is activated when at least half of
the Gauss points are broken. In addition, Gauss points where the
cumulated plastic strain is higher than \(200\%\) are considered as broken.

The evolution of the force as a function of the imposed displacement of the stricker is shown in Figure 6 where the decrease is related to the propagation of the crack. Figure 6 shows also the cumulated plastic strain field on the deformed configuration.

This report has described the extension of the
`StandardElastoViscoPlasticity`

brick to porous plasticity.
Various new stress criteria and nucleation models have been implemented
and can be used “out of the box”. A particular attention has been paid
to ease the addition of new porous stress criteria and nucleation
models.

An original staggered algorithm has been proposed and tested on ductile fracture simulations of notched specimens and Charpy tests. Note that the standard implicit scheme can still be chosen.

In many cases, the stress criterion \(\sigma_{\star}\) is an implicit function of the stress state of the form:

\[ S{\left(\sigma_{\star}, \underline{\sigma}, f\right)}=0 \qquad{(19)}\]

See Equations (1) and (2) for the case of the Gurson-Tvergaard-Needleman stress criterion and the Rousselier-Tanguy-Besson stress criterion respectively.

For a given stress state and porosity, Equation (19) may be solved by a scalar Newton algorithm.

The Newton algorithm is coupled with bisection whenever root-bracketting is possible, which considerably increase its robustness.

More precisely, if not given by the user, the algorithm tries to determines a interval \(\left[\sigma_{\star}^{\min{}};\sigma_{\star}^{\max{}}\right]\) such that:

\[ S{\left(\sigma_{\star}^{\min{}}\right)}\,S{\left(\sigma_{\star}^{\max{}}\right)}<0 \]

i.e. \(S{\left(\sigma_{\star}^{\min{}}\right)}\) and \(S{\left(\sigma_{\star}^{\max{}}\right)}\) have opposite signs. If such a pair is found, a simple bissection algorithm, whose convergence is guaranteed, may be used.

In practice, one checks if the new estimate given by the Newton algorithm lies in \(\left[\sigma_{\star}^{\min{}};\sigma_{\star}^{\max{}}\right]\). If this is the case, this new estimate is used to update the bounds of the interval.

This implementation handles properly `IEEE754`

exceptional
cases (infinite numbers, `NaN`

values), even if advanced
compilers options are used (such as `-ffast-math`

under
`gcc`

).

**Chu-Needleman strain based**[9]

\[ \begin{aligned} \dfrac{\partial A_n}{\partial \underline{\sigma}} &= \underline{0} \\ \dfrac{\partial A_n}{\partial p} &= -\dfrac{f_N}{s_N^2 \sqrt{2\pi}} {\left(\dfrac{p - \epsilon_N}{s_N} \right)} \exp{{\left(-\dfrac{1}{2} {\left(\dfrac{p - \epsilon_N}{s_N} \right)}^2 \right)} } \\ \dfrac{\partial A_n}{\partial f} &= 0 \\ \end{aligned} \]

**Chu-Needleman stress based**[9]

\[\dfrac{\partial A_n}{\partial \underline{\sigma}} = -\dfrac{f_N}{s_N^2 \sqrt{2\pi}} {\left(\dfrac{\sigma_I - \sigma_N}{s_N} \right)} \exp{{\left(-\dfrac{1}{2} {\left(\dfrac{\sigma_I - \sigma_N}{s_N} \right)}^2 \right)} } e_I \otimes e_I \]

\[\dfrac{\partial A_n}{\partial p} = 0 \]

\[\dfrac{\partial A_n}{\partial f} = 0 \]

**Power-Law strain based**

\[ \begin{aligned} \dfrac{\partial A_n}{\partial \underline{\sigma}} &= \underline{0} \\ \dfrac{\partial A_n}{\partial p} &= \dfrac{m f_N}{\epsilon_N} \left<\dfrac{p}{\epsilon_N} - 1 \right>^{m-1} \\ \dfrac{\partial A_n}{\partial f} &= 0 \\ \end{aligned} \]

**Power-Law stress based**

\[\dfrac{\partial A_n}{\partial \underline{\sigma}} = \dfrac{m f_N}{\sigma_N} \left<\dfrac{\sigma_I}{\sigma_N} - 1 \right>^{m-1} e_I \otimes e_I \]

\[\dfrac{\partial A_n}{\partial p} = 0 \]

\[\dfrac{\partial A_n}{\partial f} = 0 \]

As \(\sigma_{\star}\) is defined implicitly through \(S(\underline{\sigma},\sigma_{\star},f) =0\), the first and second order derivatives of the function S are required to compute the first and second order derivatives of \(\sigma_{\star}\):

\[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial S}{\displaystyle \partial \sigma_{\star}}} &= -2 \dfrac{\sigma_{vM}^2}{\sigma_{\star}^3} - \dfrac{3 q_1 q_2 f_{\star} \sigma_m}{\sigma_{\star}^2}\sinh{{\left(\dfrac{3}{2} q_2 \dfrac{\sigma_m}{\sigma_{\star}}\right)}} \\ \dfrac{\partial S}{\partial {\underline{\sigma}}} &= 3 \dfrac{\underline{s}}{\sigma_{\star}^2} + \dfrac{q_1 q_2 f_{\star}}{\sigma_{\star}}\sinh{{\left(\dfrac{3}{2} q_2 \dfrac{\sigma_m}{\sigma_{\star}}\right)}} \underline{I} \\ {\displaystyle \frac{\displaystyle \partial S}{\displaystyle \partial f}} &= 2 q_1 \delta \cosh{{\left(\dfrac{3}{2} q_2 \dfrac{\sigma_m}{\sigma_{\star}}\right)}} - 2 q_3 \delta f_{\star}\\ \dfrac{\partial^2 S}{\partial \underline{\sigma}^2} &= \dfrac{2}{\sigma_{\star}^2} \: {\underline{\underline{\mathbf{M}}}} + \dfrac{ q_1 \: q_2^2 \: f_{\star}}{2 \sigma_{\star}^2} \cosh{{\left(\dfrac{3}{2} q_2 \dfrac{\sigma_m}{\sigma_{\star}}\right)}} \underline{I} \otimes \underline{I}\\ \dfrac{\partial^2 S}{\partial \underline{\sigma}\partial \sigma_{\star}} &= - \dfrac{6 \underline{s}}{\sigma_{\star}^3} - \dfrac{ q_1 \: q_2 \: f_{\star}}{\sigma_{\star}^2} \sinh{{\left(\dfrac{3}{2} q_2 \dfrac{\sigma_m}{\sigma_{\star}}\right)}} \underline{I} - \dfrac{3 \: q_1 \: q_2^2 \: f_{\star} \: \sigma_m}{2 \: \sigma_{\star}^3} \cosh{{\left(\dfrac{3}{2} q_2 \dfrac{\sigma_m}{\sigma_{\star}}\right)}} \underline{I}\\ \dfrac{\partial^2 S}{\partial \underline{\sigma}\partial f} &= \dfrac{q_1 \: q_2 \: \delta}{\sigma_{\star}} \sinh{{\left(\dfrac{3}{2} q_2 \dfrac{\sigma_m}{\sigma_{\star}}\right)}} \underline{I} \\ \dfrac{\partial^2 S}{\partial \sigma_{\star}\partial f} &= - \dfrac{3 \: q_1 \: q_2 \: \delta \: \sigma_m}{\sigma_{\star}^2} \sinh{{\left(\dfrac{3}{2} q_2 \dfrac{\sigma_m}{\sigma_{\star}}\right)}} \\ \dfrac{\partial^2 S}{\partial \sigma_{\star}^2} &= 6 \dfrac{\sigma_{vM}^2}{\sigma_{\star}^4} + \dfrac{6 q_1 q_2 f_{\star} \sigma_m}{\sigma_{\star}^3}\sinh{{\left(\dfrac{3}{2} q_2 \dfrac{\sigma_m}{\sigma_{\star}}\right)}} + \dfrac{9 q_1 q_2^2 f_{\star} \sigma_m^2}{2\sigma_{\star}^4}\cosh{{\left(\dfrac{3}{2} q_2 \dfrac{\sigma_m}{\sigma_{\star}}\right)}} \end{aligned} \]

The first and second order derivatives of \(\sigma_{\star}\) are:

\[{\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}} = -{\left({\displaystyle \frac{\displaystyle \partial S}{\displaystyle \partial \sigma_{\star}}} \right)}^{-1} \dfrac{\partial S}{\partial {\underline{\sigma}}}\]

\[{\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial f}} = -{\left({\displaystyle \frac{\displaystyle \partial S}{\displaystyle \partial \sigma_{\star}}} \right)}^{-1} \dfrac{\partial S}{\partial {f}}\]

\[ \dfrac{\partial^2 \sigma_{\star}}{\partial \underline{\sigma}^2} = - \bigg({\displaystyle \frac{\displaystyle \partial S}{\displaystyle \partial \sigma_{\star}}} \bigg)^{-1} \bigg(\dfrac{\partial^2 S}{\partial \underline{\sigma}^2} + \dfrac{\partial^2 S}{\partial \underline{\sigma}\partial \sigma_{\star}} \otimes {\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}} \bigg) + \bigg({\displaystyle \frac{\displaystyle \partial S}{\displaystyle \partial \sigma_{\star}}} \bigg)^{-2} \bigg({\displaystyle \frac{\displaystyle \partial S}{\displaystyle \partial \underline{\sigma}}} \bigg) \otimes \bigg(\dfrac{\partial^2 S}{\partial \underline{\sigma}\partial \sigma_{\star}} + \dfrac{\partial^2 S}{ \partial \sigma_{\star}^2} {\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}} \bigg) \]

\[\dfrac{\partial^2 \sigma_{\star}}{\partial \underline{\sigma}\partial f} = - \bigg(\dfrac{\partial S}{\partial \sigma_{\star}} \bigg)^{-1} \bigg(\dfrac{\partial^2 S}{\partial \underline{\sigma}\partial f} + \dfrac{\partial^2 S}{\partial f \partial \sigma_{\star}} {\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}} \bigg) + \bigg(\dfrac{\partial S}{\partial \sigma_{\star}} \bigg)^{-2} \bigg(\dfrac{\partial^2 S}{\partial \underline{\sigma}\partial \sigma_{\star}} + \dfrac{\partial^2 S}{ \partial \sigma_{\star}^2} {\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}} \bigg) \bigg({\displaystyle \frac{\displaystyle \partial S}{\displaystyle \partial f}} \bigg) \]

The first and second order derivatives of \(\sigma_{\star}\) are computed using \(S(\underline{\sigma},\sigma_{\star},f) =0\), which requires the first and second order derivatives of the function S:

\[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial S}{\displaystyle \partial \sigma_{\star}}} &= -\dfrac{\sigma_{vM}}{(1- f )\sigma_{\star}^2} - f D_R q_R \dfrac{\sigma_{m}}{(1- f )\sigma_{\star}^2} \exp{{\left(\dfrac{3}{2} q_R \dfrac{\sigma_m}{(1 - f )\sigma_{\star}} \right)}}\\ \dfrac{\partial S}{\partial {\underline{\sigma}}} &= \dfrac{3}{2}\dfrac{\underline{s}}{(1- f ) \sigma_{vM} \sigma_{\star}} + \dfrac{f D_R q_R}{3} \dfrac{\underline{I}}{(1- f )\sigma_{\star}} \exp{{\left(\dfrac{3}{2} q_R \dfrac{\sigma_m}{(1 - f )\sigma_{\star}} \right)}} \\ {\displaystyle \frac{\displaystyle \partial S}{\displaystyle \partial f}} &= \dfrac{\sigma_{vM}}{\sigma_{\star}(1-f)^2} + \dfrac{2}{3} D_R \exp{{\left(\dfrac{3}{2} q_R \dfrac{\sigma_m}{(1 - f )\sigma_{\star}} \right)}} {\left(1 + \dfrac{3 f q_R \sigma_m}{2 \sigma_{\star}(1-f)^2} \right)} \\ \dfrac{\partial^2 S}{\partial \underline{\sigma}^2} &= \dfrac{3}{2(1-f)\sigma_{\star}\sigma_{vM} } {\left( \dfrac{2}{3} {\underline{\underline{\mathbf{M}}}} - \dfrac{3}{2} \dfrac{\underline{s}}{\sigma_{vM}} \otimes \dfrac{\underline{s}}{\sigma_{vM}} \right)} + \dfrac{f D_R q_R^2 }{6 (1-f)^2 \sigma_{\star}^2} \underline{I} \otimes \underline{I} \exp{{\left(\dfrac{3}{2} q_R \dfrac{\sigma_m}{(1 - f )\sigma_{\star}} \right)}}\\ \dfrac{\partial^2 S}{\partial \underline{\sigma}\partial \sigma_{\star}} &= - \dfrac{3 \underline{s}}{2 (1-f) \sigma_{vM} \sigma_{\star}^2 } - {\left(1 + \dfrac{3 q_R \sigma_m}{2 (1- f) \sigma_{\star}}\right)} \dfrac{f D_R q_R}{3 (1-f) \sigma_{\star}^2} \exp{{\left(\dfrac{3}{2} q_R \dfrac{\sigma_m}{(1 - f )\sigma_{\star}} \right)}} \underline{I}\\ \dfrac{\partial^2 S}{\partial \underline{\sigma}\partial f} &= \dfrac{3 \underline{s}}{2 (1-f)^2 \sigma_{vM} \sigma_{\star}} + {\left(1 + \dfrac{3 q_R \sigma_m}{2 (1 - f) \sigma_{\star}} + \dfrac{1-f}{f}\right)} \dfrac{f D_R q_R}{3 (1-f)^2 \sigma_{\star}} \exp{{\left(\dfrac{3}{2} q_R \dfrac{\sigma_m}{(1 - f )\sigma_{\star}} \right)}} \underline{I} \\ \dfrac{\partial^2 S}{\partial \sigma_{\star}\partial f} &= -\dfrac{\sigma_{vM} }{(1-f)^2 \sigma_{\star}^2 } - {\left(1 + \dfrac{3 q_R \sigma_m}{2 (1 - f) \sigma_{\star}} + \dfrac{1-f}{f}\right)} \dfrac{f D_R q_R \sigma_m}{ (1-f)^2 \sigma_{\star}^2} \exp{{\left(\dfrac{3}{2} q_R \dfrac{\sigma_m}{(1 - f )\sigma_{\star}} \right)}} \\ \dfrac{\partial^2 S}{\partial \sigma_{\star}^2} &= \dfrac{2\sigma_{vM}}{(1-f)\sigma_{\star}^3} + {\left(1 + \dfrac{3 q_R \sigma_m}{4 (1-f) \sigma_{\star}}\right)} \dfrac{2fD_Rq_R \sigma_m}{(1-f)\sigma_{\star}^3} \exp{{\left(\dfrac{3}{2} q_R \dfrac{\sigma_m}{(1 - f )\sigma_{\star}} \right)}} \end{aligned} \]

The first and second order derivatives of \(\sigma_{\star}\) are:

\[{\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}} = -{\left({\displaystyle \frac{\displaystyle \partial S}{\displaystyle \partial \sigma_{\star}}} \right)}^{-1} \dfrac{\partial S}{\partial {\underline{\sigma}}}\]

\[{\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial f}} = -{\left({\displaystyle \frac{\displaystyle \partial S}{\displaystyle \partial \sigma_{\star}}} \right)}^{-1} \dfrac{\partial S}{\partial {f}}\]

\[ \dfrac{\partial^2 \sigma_{\star}}{\partial \underline{\sigma}^2} = - \bigg({\displaystyle \frac{\displaystyle \partial S}{\displaystyle \partial \sigma_{\star}}} \bigg)^{-1} \bigg(\dfrac{\partial^2 S}{\partial \underline{\sigma}^2} + \dfrac{\partial^2 S}{\partial \underline{\sigma}\partial \sigma_{\star}} \otimes {\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}} \bigg) + \bigg({\displaystyle \frac{\displaystyle \partial S}{\displaystyle \partial \sigma_{\star}}} \bigg)^{-2} \bigg({\displaystyle \frac{\displaystyle \partial S}{\displaystyle \partial \underline{\sigma}}} \bigg) \otimes \bigg(\dfrac{\partial^2 S}{\partial \underline{\sigma}\partial \sigma_{\star}} + \dfrac{\partial^2 S}{ \partial \sigma_{\star}^2} {\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}} \bigg) \]

\[\dfrac{\partial^2 \sigma_{\star}}{\partial \underline{\sigma}\partial f} = - \bigg(\dfrac{\partial S}{\partial \sigma_{\star}} \bigg)^{-1} \bigg(\dfrac{\partial^2 S}{\partial \underline{\sigma}\partial f} + \dfrac{\partial^2 S}{\partial f \partial \sigma_{\star}} {\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}} \bigg) + \bigg(\dfrac{\partial S}{\partial \sigma_{\star}} \bigg)^{-2} \bigg(\dfrac{\partial^2 S}{\partial \underline{\sigma}\partial \sigma_{\star}} + \dfrac{\partial^2 S}{ \partial \sigma_{\star}^2} {\displaystyle \frac{\displaystyle \partial \sigma_{\star}}{\displaystyle \partial \underline{\sigma}}} \bigg) \bigg({\displaystyle \frac{\displaystyle \partial S}{\displaystyle \partial f}} \bigg) \]

The Gurson-Tvergaard-Needleman stress criterion models the void
coalescence by the effective porosity calculated as shown `above`

. The transition from the
pre-coalescence to the coalescence phase might cause numerical problems
if the \(\delta\) parameter is high
enough. This is due to the inflection point in the bilinear function
describing \(f_{\star}\) as a function
of \(f\). Therefore, the first
derivative of the effective porosity with respect to the porosity (\(\dfrac{\partial f_{\star}}{\partial f}\))
undergoes a sudden jump that disturbs the Newton-Raphson algorithm.

A solution is proposed to avoid that derivative jump at the inflection point : to replace the bilinear function of \(f_{\star}\) by a polynominal function of 5th degree. Therefore, \(f_{\star}\) is defined such that :

\[ f_{\star} = \left\{ \begin{array}{ll} \delta f & \mbox{if } f < (f_c - \epsilon) \\ \delta (a^5 f + b^4 f + c^3 f + d^2 f + e f + g) & \mbox{if } (f_c - \epsilon) \leq f \leq (f_c + \epsilon) \\ f_c + \delta (f - f_c) & \mbox{otherwise} \end{array} \right. \]

and :

\[ \delta = \left\{ \begin{array}{ll} 1 & \mbox{if } f \leq (f_c + \epsilon) \\ \dfrac{\dfrac{1}{q_1} - f_c }{f_r - f_c} & \mbox{otherwise} \end{array} \right. \]

where \(\epsilon\) is the linearization coefficient that allows to have a smooth transition on the \(\partial f_{\star}\) as function of \(f\). A big \(\epsilon\) value will allow a higher contribution of the polynominal function and thus, smoother numerical convergence.

For numerical efficiency reasons, the polynominal function was implemented using the Horner method to be :

\[ f (f ( f (f (a f + b) + c) + d) + e) + g\]

`Cast3M`

The staggered approach improves robustness of the numerical
integration of the porous constitutive equations, especially if
sub-stepping is allowed. However, integration may still fail due to too
severe input parameters. As integration failure is not handled in the
finite element code `Cast3M`

, additional stategies are
required at the structural scale to avoid such problems. An efficient
stategy is to adjust the time-step based on the maximal increase of
porosity such as: \[
\Delta t_{i+1} = \min{\left( \frac{\Delta f_{max}}{\Delta f_{i}}
\Delta t_i , \Delta t_{ref} \right)}
\] This equation is implemented in two steps. First, the
following code is added to the `PERSO1`

user procedure to
compute the maximal increase of porosity at the last time-step on
unbroken Gauss points.

```
'SI'( ('DIME' TAB1.DEPLACEMENTS) > 1);
= TAB1.'F_PREC';
TMP = 'EXCO' TAB1.ESTIMATION.VARIABLES_INTERNES 'F' 'SCAL';
TMP2 = 'EXCO' TAB1.ESTIMATION.VARIABLES_INTERNES 'BROK' 'SCAL';
TMP3
= TMP*(1 '-' TMP3);
TMPA = TMP2*(1 '-' TMP3);
TMP2A
= 'ABS' (('MAXIMUM' TMPA) '-' ('MAXIMUM' TMP2A));
DF
'F_PREC' = TMP2;
TAB1.
'SINON';
= 0;
DF
'FINSI';
.WTABLE.'DF' = DF;
TAB1.WTABLE.'DFMAX' = TAB1.'DFMAX'; TAB1
```

Then, the `PASAPAS`

procedure is modified according
to:

```
'SI' WTAB.'PAS_AJUSTE';
= DT_AVANT;
DTV = WTAB.'DF';
DF = WTAB.'DFMAX';
DFMAX = DF '/' DFMAX;
DDF 'SI' ('EXISTE' WTAB ISOUSPAS);
= WTAB.ISOUSPAS;
ISP 'SINON';
= 0;
ISP 'FINSI';
'SI' WTAB.'CONV';
'SI' ((DDF > 1));
= DTV '/' DDF;
DTV 'FINSI';
'SI' ((DDF < 1) 'ET' (ISP 'EGA' 0));
'SI' IREREDU;
'SI' (DDF > 0.5);
= DTV '/' DDF;
DTV 'SINON';
= DTV '/' 0.5;
DTV 'FINSI';
'FINSI';
'FINSI';
= DTV;
DT_AVANT = VRAI;
IREREDU 'SINON';
=0.0000000001D0* DT_AVANT;
DTV= DT_AVANT /2.;
DT_AVANT=FAUX;
IREREDU'FINSI';
= DTV* 1.0000000001 + TEMP0;
TTI 'SI' ( TTI '<EG' TIV );
=TTI;
TI=0;
PASFINAL'FINSI';
'FINSI';
```

The time-step is adjusted in order to target a maximal local increase
of porosity lower or equal to an user-defined value `DFMAX`

.
In addition, increasing the time-step is allowed only if the convergence
of the last time-step was achieved without sub-stepping (through
`CONVERGENCE_FORCEE`

).

In addition, global convergence may be difficult to achieve in cases
where heavily-distorted elements are present. Element-deletion is
detailed below, but a slight modification of the `UNPAS`

procedure can be made to lower the precision of the computation when
sub-steps are done, allowing to avoid premature end of the
calculation:

```
*Nombre maximum de sous-pas atteint? si oui, arret de pasapas
'ISOUSPAS' = WTAB.'ISOUSPAS' + 1;
WTAB.= ZPREC * 2;
ZPREC = ZPRECM * 2; ZPRECM
```

Finally, an example of element deletion that can be added to the
`PERSO1`

procedure can be added, based on the number of
broken Gauss points and cumulated plastic strain:

```
= TAB1.WTABLE.'MO_TOT';
MO1 = TAB1.WTABLE.'MA_TOT';
MA1 = 'EXTRAIRE' MO1 'MAIL';
MAIL1 = 'EXCO' (TAB1.ESTIMATION.VARIABLES_INTERNES) 'BROK' 'SCAL';
TMP1 = 'EXCO' (TAB1.ESTIMATION.VARIABLES_INTERNES) 'P' 'SCAL';
TMP2 = TMP2 'MASQ' SUPE (TAB1.'BROKEN_STRAIN');
TMP2
= (TMP1 + TMP2) MASQ SUPE 0.;
TMP1 = 'CHANGER' GRAVITE MO1 TMP1 'STRESSES';
TMP1 = TMP1 'MASQUE' 'EGSU' 0.5;
TMP1
= TMP1 'ELEM' SUPE 0.;
MEBR = 'DIFF' MAIL1 MEBR;
MESO = 'NBEL' MEBR;
NBRO
'DBRO' = TAB1.'DBRO' '+' NBRO;
TAB1.'NBRO' = TAB1.'NBRO' '+' NBRO;
TAB1.
'SI' ('>' NBRO 0);
= 'REDU' MO1 MESO;
MO2 .WTABLE.'MO_TOT' = MO2;
TAB1.WTABLE.'MO_TOTAL' = MO2;
TAB1.WTABLE.'MOD_MEC' = MO2;
TAB1
'FINSI';
```

`porosity_evolution`

sectionThis appendix is dedicated to the description of the various options
that can be declared in the `porosity_evolution`

section
inside the declaration of the `StandardElastoViscoPlasticity`

brick:

`save_individual_porosity_increase`

: a boolean value stating if the each contribution to the porosity increase must be stored in a dedicated state variable or auxiliary state variable. By default, this option is set to`false`

.`elastic_contribution`

: a boolean value stating if the elastic contribution to the porosity growth must be taken into account. By default, this option is set to`false`

.`nucleation_model`

: this option introduces a new nucleation model.`algorithm`

: a string value which allow the selection of the resolution algorithm. Valid values are:`standard implicit scheme`

`standard_implicit_scheme`

`StandardImplicitScheme`

`staggered scheme`

`staggered_scheme`

`StaggeredScheme`

The`algorithm`

may have suboptions, as described in Section

`safety_factor_for_the_upper_bound_of_the_porosity`

: this option allows setting floatting-point value, denoted \(\alpha\), which is used to determine the maximum value of the porosity at the middle of the time step as \(\alpha\,f_{r}\), where \(f_{r}\) is the upper bound of the porosity. The upper bound of the porosity is equal the minimum of the coalescence porosity of all selected stress criteria or is equal to \(1\) if no stress criterion declares a coalescence porosity. By default, the value of this option is \(0.985\).`safety_factor_for_the_upper_bound_of_the_porosity_for_fracture_detection`

: this option allows setting floatting-point value, denoted \(\alpha\), which is used by the staggered algorithm to detect the failure of the material. The failure is detected when the porosity converges to a value greater than \(\alpha\,f_{r}\) where \(f_{r}\) is the upper bound of the porosity. By default, the value of this option is \(0.984\).

`algorithm`

optionWhen a staggered scheme is used, two numerical parameters can be specified:

`convergence_criterion`

: the value of the convergence criterion of the staggered algorithm. The staggered scheme stops when the difference between two estimates of the porosity is lower than this value. The default value is \(1.e-10\).`maximum_number_of_iterations`

: the maximum number of iterations allows for the staggered scheme. The default value is \(100\).

Here is a simple example on how to specify those numerical parameters:

```
: {
porosity_evolution: "staggered_scheme" {
algorithm : 1.e-9,
convergence_criterion: 200
maximum_number_of_iterations}
}
```

1.

Gurson, A. L.
Continuum theory of ductile rupture by void nucleation and growth:
Part I - Yield criteria and flow
rules for porous ductile media. *J. Eng. Mat. and Tech.* 1977.
Vol. 99, p. 2–15.

2.

Tvergaard, V.
and Needleman, A. Analysis of the
cup-cone fracture in a round tensile bar. *Acta. Metall.* 1984.
Vol. 32, p. 157–169.

3.

Castañeda, P.
Ponte. The effective mechanical properties of nonlinear isotropic
composites. *Journal of the Mechanics and Physics of Solids*. 1
January 1991. Vol. 39, no. 1, p. 45–71. DOI 10.1016/0022-5096(91)90030-R.
Available from: http://www.sciencedirect.com/science/article/pii/002250969190030R

4.

Monerie, Yann
and Gatt, Jean-Marie. Overall
viscoplastic behavior of non-irradiated porous nuclear ceramics.
*Mechanics of Materials*. July 2006. Vol. 38, no. 7, p. 608–619.
DOI 10.1016/j.mechmat.2005.11.004.
Available from: http://www.sciencedirect.com/science/article/pii/S0167663605001882

5.

Salvo, Maxime,
Sercombe, Jérôme, Helfer, Thomas, Sornay, Philippe and Désoyer, Thierry. Experimental characterization
and modeling of UO2 grain boundary cracking at high
temperatures and high strain rates. *Journal of Nuclear
Materials*. May 2015. Vol. 460, p. 184–199. DOI 10.1016/j.jnucmat.2015.02.018.
Available from: http://www.sciencedirect.com/science/article/pii/S0022311515001130

6.

Rousselier, G.
Finite deformation constitutive relations including ductile fracture
damage. In : *Nemat-Nasser (ed.)
Three Dimensional Constitutive
Relations and Ductile
Fracture*. North Holland, 1981. p. 331–355.

7.

Berdin, C.,
Besson, J., Bugat, S., Desmorat, R., Feyel, F., Forest, E., S. Lorentz, Maire, E., Pardoen, T., Pineau, A. and Tanguy, B. Local approach to fracture. Ecole
des Mines de Paris, 2004.

8.

Tanguy, B. and
Besson, J. An extension of the
Rousselier model to viscoplastic temperature dependent
materials. *Int. J. Fracture*. 2002. Vol. 116, p. 81–101.

9.

Chu, C. C. and
Needleman, A. Void nucleation effects in
biaxially stretched sheets. *Journal of Engineering Materials and
Technology*. 1980. Vol. 102, no. 3, p. 249–256.

10.

Helfer, Thomas.
Assisted computation of the consistent tangent operator of behaviours
integrated using an implicit scheme. Theory and implementation in
MFront. Documentation of mfront. CEA, 2020.
Available from: https://www.researchgate.net/publication/342721072_Assisted_computation_of_the_consistent_tangent_operator_of_behaviours_integrated_using_an_implicit_scheme_Theory_and_implementation_in_MFront

11.

Zhang, Y.
Modélisation et simulation robuste de l’endommagement ductile. PhD
thesis. Paris, France : Université de recherche Paris Sciences et
Lettres, 2016.

12.

Helfer, T.
Implementation of a multi-surface, compressible and perfect plastic
behaviour using the drucker-prager yield criterion and a cap.
TFEL/MFront. 21 November 2017. Available from:
https://thelfer.github.io/tfel/web//drucker-prager-cap.html

13.

Petit, T., Besson, J., Ritter, C., Colas, K., Helfen, L. and Morgeneyer, T. F. Effect of hardening on
toughness captured by stress-based damage nucleation in 6061 aluminum
alloy. *Acta Metall.* 2019. Vol. 180, p. 349–365.
DOI https://doi.org/10.1016/j.actamat.2019.08.055.

14.

Tanguy, B.
Modélisation de l’essai charpy par l’approche locale de la rupture:
Application au cas de l’acier 16MND5 dans le domaine de la transition.
PhD thesis. Ecole Nationale Supérieure des Mines de Paris, 2001.