StressCriterion
and InelasticFlow
interfacesCast3M
porosity_evolution
sectionConstitutive 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:
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:
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 : };
Note For 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:
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:
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\).
Note
The 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.
\[ 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\).
\[ 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.
\[ 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.
\[ 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.
\[ 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)} \]
\[ 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:
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 prediction
When 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.1.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:
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:
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)\]
Note
In 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:
\(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)} \]
Note
To 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;
"StaggeredSchemeIterationCounter");
niter.setEntryName(
....
@UpdateAuxiliaryStateVariables{
niter = staggered_scheme_iteration_counter; }
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;
"NumberOfEvaluationOfTheImplicitScheme");
neval.setEntryName(
...
@InitLocalVariables{
...0;
neval =
}
...
@Integrator{
++neval; }
Note
If 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:
@ComputeStress
code block. This code block is automatically handled by the brick.@Integrator
code block.@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:
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:
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.
Note
Implicitly, the brick expects that a stress criterion coupled with porosity is not deviatoric, i.e., that this stress criterion contributes to the porosity growth. The
isNormalDeviatoric
method of such a stress criterion must returntrue
. This is checked in theInelasticFlowBase
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:
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 :
} };
Note
It 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;
"SolidSwelling");
ss.setGlossaryName(
@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:
upper_bound_of_the_porosity = min(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 : }
Note
Inelastic 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
).
\[ \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} \]
\[\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 \]
\[ \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} \]
\[\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 Sectionsafety_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:
} }