PorousViscoplasticity
\newcommand{\tsigma}{\underline{\sigma}} \newcommand{\paren}[1]{{\left(#1\right)}} \newcommand{\trace}[1]{{\mathrm{tr}\paren{#1}}} \newcommand{\tenseur}[1]{\underline{#1}} \newcommand{\tenseurq}[1]{\underline{\underline{\mathbf{#1}}}} \newcommand{\tepsilonto}{\tenseur{\varepsilon}^{\mathrm{to}}} \newcommand{\tdepsilonto}{\tenseur{\dot{\varepsilon}}^{\mathrm{to}}} \newcommand{\tepsilonel}{\tenseur{\varepsilon}^{\mathrm{el}}} \newcommand{\tdepsilonel}{\tenseur{\dot{\varepsilon}}^{\mathrm{el}}} \newcommand{\tepsilonvis}{\tenseur{\varepsilon}^{\mathrm{vis}}} \newcommand{\tdepsilonvis}{\tenseur{\dot{\varepsilon}}^{\mathrm{vis}}} \newcommand{\tepsilonp}{\tenseur{\varepsilon}^{\mathrm{p}}} \newcommand{\tdepsilonp}{\tenseur{\dot{\varepsilon}}^{\mathrm{p}}} \newcommand{\sstar}{\sigma_{\star}} \newcommand{\pr}{\mathrm{\sigma_{m}}} \newcommand{\sigmaeq}{\sigma_{\mathrm{eq}}} \newcommand{\bts}[1]{{\left.#1\right|{t}}} \newcommand{\mts}[1]{{\left.#1\right|{t+\theta,\Delta,t}}} \newcommand{\ets}[1]{{\left.#1\right|_{t+\Delta,t}}} \newcommand{\Frac}[2]{{{\displaystyle \frac{\displaystyle #1}{\displaystyle #2}}}}
Constitutive equations for porous metallic equations are used to perform
ductile fracture simulations, adding the porosity as an internal state
variable. Some implementations of these models are available in MFront
(Gurson, GTN, ...). This page describes how one may extend the
StandardElastoViscoPlasticity
brick
to this class of behaviours, allowing a declarative syntax similar as:
@Brick "StandardElastoViscoPlasticity" {
stress_potential : "Hooke" {
young_modulus : 200e9,
poisson_ratio : 0.3
},
inelastic_flow : "Plastic" {
criterion : "GursonTvergaardNeedleman" {
q1 : 1.5,
q2 : 1.0,
q3 : 2.2,
fc : 0.01,
fr : 0.1
},
isotropic_hardening : "Linear" {
R0 : 150e6
}
},
porosity_evolution : {
nucleation_model : "Chu_Needleman" {
An : 0.01,
pn : 0.1,
sn : 0.1
},
growth_model : "StandardPlasticModel"
};Note For backward compatibility, effect of the porosity on plastic flows will not be taken into account if no porosity_evolution is declared.
A porous plastic model is generally defined by a stress criterion (\sstar) which depends explicitly on the stress and the porosity: [ \sstar = \phi(\tsigma,f) ]
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.
This class of models encompasses:
- the work of Ponte-Castaneda (see @castaneda_effective_1991), which is the theorical basis of several viscoplastic behaviours used in uranium dioxyde fuels [@monerie_overall_2006;@salvo_experimental_2015].
- the work of Rousselier.
In the second class, the flow rule is affected by a correction factor (1-f) (see @mealor_2004):
where
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{\tsigma:\tdepsilonp = \dot{\lambda} \tsigma, \colon,\frac{\partial \sstar}{\partial \tsigma}} = \dot{\lambda} \sstar = (1 - f)R(p) \dot{p} \Rightarrow \dot{\lambda} = (1 - f) \dot{p} $$
The evolution of porosity is composed of two terms: one due to plastic
flow and another one due to nucleation
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.
Assuming a single porous stress criterion, the constitutive equations lead to the following system to be solved for an elastoplastic evolution:
A semi-implicit method is used to solve these equations according to the
variables
$$ \begin{aligned} \mathcal{R}{\Delta \tepsilonel} = \Delta \tepsilonel + (1 - f \big|{t+\theta \Delta t}) \Delta p , \frac{\partial \sstar}{\partial \tsigma}\big|{t+\theta \Delta t} - \Delta \underline{\epsilon}{to} &= 0 \ \mathcal{R}{\Delta p} = \sstar - R(p\big|{t+\theta \Delta t}) &= 0 \ \mathcal{R}{\Delta f} = \Delta f - (1 - f\big|{t+\theta \Delta t})^2 \Delta p , \trace{\frac{\partial \sstar}{\partial \tsigma}\big|{t+\theta \Delta t}} - \sum_i A_n^i\big|{t+\theta \Delta t} \Delta p &=0 \end{aligned} $$
Solving this system using a Newton-Raphson algorithm requires computing the Jacobian matrix:
$$ \left[ \begin{array}{lll} \displaystyle{ \frac{\partial \mathcal{R}{\Delta \tepsilonel}}{\partial \Delta \tepsilonel}} & \displaystyle{ \frac{\partial \mathcal{R}{\Delta \tepsilonel}}{\partial \Delta f}} & \displaystyle{ \frac{\partial \mathcal{R}{\Delta \tepsilonel}}{\partial \Delta p} } \ \displaystyle{ \frac{\partial \mathcal{R}{\Delta f}}{\partial \Delta \tepsilonel} } & \displaystyle{ \frac{\partial \mathcal{R}{\Delta f}}{\partial \Delta f} } & \displaystyle{ \frac{\partial \mathcal{R}{\Delta f}}{\partial \Delta p}} \ \displaystyle{ \frac{\partial \mathcal{R}{\Delta p}}{\partial \Delta \tepsilonel}} & \displaystyle{ \frac{\partial \mathcal{R}{\Delta p}}{\partial \Delta f} } & \displaystyle{ \frac{\partial \mathcal{R}_{\Delta p}}{\partial \Delta p}} \ \end{array} \right] $$
where the different terms can be written as [@mealor_2004]:
$$ \begin{aligned} \frac{\partial \mathcal{R}{\Delta \tepsilonel}}{\partial \Delta \tepsilonel} &= \underline{\bf{1}} + \theta \Delta p ; \bigg(1 - f \big|{t+\theta \Delta t}\bigg) \frac{\partial^2 \sstar}{\partial \tsigma^2}|{t+\theta \Delta t} : \underline{\bf{D}}|{t+\theta \Delta t} \ \frac{\partial \mathcal{R}{\Delta \tepsilonel}}{\partial \Delta p} &= \bigg(1 - f \big|{t+\theta \Delta t}\bigg) \frac{\partial \sstar}{\partial \tsigma}|{t+\theta \Delta t} \ \frac{\partial \mathcal{R}{\Delta \tepsilonel}}{\partial \Delta f} &= \theta \Delta p ; \bigg( \bigg(1 - f \big|{t+\theta \Delta t}\bigg) \frac{\partial^2 \sstar}{\partial \tsigma \partial f} |{t+\theta \Delta t} - \frac{\partial \sstar}{\partial \tsigma}|{t+\theta \Delta t} \bigg) \ \frac{\partial \mathcal{R}{\Delta p}}{\partial \Delta \tepsilonel} &= \theta ; \frac{\partial \sstar}{\partial \tsigma}|{t+\theta \Delta t} : \underline{\bf{D}}|{t+\theta \Delta t} \ \frac{\partial \mathcal{R}{\Delta p}}{\partial \Delta p} &= - \theta ; \frac{d R(p)}{d p} \ \frac{\partial \mathcal{R}{\Delta p}}{\partial \Delta f} &= \theta ; \frac{\partial \sstar}{\partial f} |{t+\theta \Delta t} \ \frac{\partial \mathcal{R}{\Delta f}}{\partial \Delta \tepsilonel} &= - \theta ; \Delta p ; \bigg(1 - f|{t+\theta \Delta t}\bigg)^2 ; \bigg(\frac{\partial^2 \sstar}{\partial \tsigma^2} |{t+\theta \Delta t} : \underline{\bf{D}}|{t+\theta \Delta t} \bigg) : \underline{1} - \theta \Delta p \sum_i \frac{\partial A_n^i}{\partial \tsigma} |{t+\theta \Delta t} : \underline{\bf{D}}|{t+\theta \Delta t} \ \frac{\partial \mathcal{R}{\Delta f}}{\partial \Delta p} &= - \bigg(1 - f|{t+\theta \Delta t}\bigg)^2 ; \frac{\partial \sstar}{\partial \tsigma}|{t+\theta \Delta t} : \underline{1} - \theta \Delta p \sum_i \frac{\partial A_n^i}{\partial p} |{t+\theta \Delta t} - \sum_i A_n^i |{t+\theta \Delta t} \ \frac{\partial \mathcal{R}{\Delta f}}{\partial \Delta f} &= 1 - \theta ; \Delta p ; \bigg(1 - f|{t+\theta \Delta t}\bigg) \bigg(-2 ; \frac{\partial \sstar}{\partial \tsigma} + \bigg(1 - f|{t+\theta \Delta t}\bigg) ; \frac{\partial^2 \sstar}{\partial \tsigma \partial f} \bigg) \underline{I} - \theta \Delta p \sum_i \frac{\partial A_n^i}{\partial f} |{t+\theta \Delta t} \ \end{aligned} $$
Therefore, the definition of a porous model requires to give the following expressions:
as well as the following expressions for each nucleation mechanism:
Different stress criteria implemented in the brick are described below.
- Green [@fritzen_computational_2013]
where
Note This model is a particular case of [@fritzen_computational_2013] and is implemented mainly as a simple working example of a porous stress criterion.
- GursonTvergaardNeedleman [@tvergaardneedleman1984]
which defines implicitly
and :
The parameters of the model are:
For
As
The first and second order derivatives of
$$ \frac{\partial^2 \sstar}{\partial \tsigma^2} = - \bigg(\frac{\partial S}{\partial \sstar} \bigg)^{-1} \bigg(\frac{\partial^2 S}{\partial \tsigma^2} + \frac{\partial^2 S}{\partial \tsigma \partial \sstar} \otimes \frac{\partial \sstar}{\partial \tsigma} \bigg)
- \bigg(\frac{\partial S}{\partial \sstar} \bigg)^{-2} \bigg(\frac{\partial S}{\partial \tsigma} \bigg) \otimes \bigg(\frac{\partial^2 S}{\partial \tsigma \partial \sstar} + \frac{\partial^2 S}{ \partial \sstar^2} \frac{\partial \sstar}{\partial \tsigma} \bigg) $$
- RousselierTanguyBesson [@tanguybesson2002]
The effective stress
where
The first and second order derivatives of
The first and second order derivatives of
$$ \frac{\partial^2 \sstar}{\partial \tsigma^2} = - \bigg(\frac{\partial S}{\partial \sstar} \bigg)^{-1} \bigg(\frac{\partial^2 S}{\partial \tsigma^2} + \frac{\partial^2 S}{\partial \tsigma \partial \sstar} \otimes \frac{\partial \sstar}{\partial \tsigma} \bigg)
- \bigg(\frac{\partial S}{\partial \sstar} \bigg)^{-2} \bigg(\frac{\partial S}{\partial \tsigma} \bigg) \otimes \bigg(\frac{\partial^2 S}{\partial \tsigma \partial \sstar} + \frac{\partial^2 S}{ \partial \sstar^2} \frac{\partial \sstar}{\partial \tsigma} \bigg) $$
- GenericThomason
where
The first and second order derivatives of
where
The Thomason model corresponds to:
The evolutions of the void aspect ratio
Different nucleation terms are available and are described below. Multiple nucleation terms can be used simultaneously.
- Chu-Needleman strain based [@chuneedleman1980]
where
- Chu-Needleman stress based [@chuneedleman1980]
where
- Power-Law strain based
where
- Power-Law stress based
where
\def\doubleunderline#1{\underline{\underline{#1}}}
Newton-Raphson algorithm may fail finding the solution of the non-linear system of equations, especially when the time step leads to strong increase of the porosity. In order to improve the robustness of the numerical integration, two strategies are used:
-
In case of non convergence, the system of equations is solved first assuming a fixed porosity, and upon convergence a second time allowing the porosity to vary.
-
If the previous strategy does not lead to convergence, the system of equations is solved first assuming a fixed porosity, and the equation related to porosity is integrated explicitly.
The 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 isCoupledWithPorosityEvolution has been added to
those interface. This method returns a boolean value and does take any
argument.
For inelastic flows, the default implementation of the
isCoupledWithPorosityEvolution, declared in the InelasticFlowBase
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
isNormalDeviatoricmethod of such a stress criterion must returntrue. This is checked in theInelasticFlowBaseclass.
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 getPorosityEffectOnFlowRule
which return one of the two following values:
NO_POROSITY_EFFECT_ON_FLOW_RULESTANDARD_POROSITY_CORRECTION_ON_FLOW_RULE
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_evolutionsection is explicitly declared, even if empty.
For example, this declaration will automatically lead to a porosity growth:
@Brick StandardElastoViscoPlasticity{
stress_potential : "Hooke" {
young_modulus : 150.e3,
poisson_ratio : 0.3
},
inelastic_flow : "Plastic" {
criterion : "Gurson1975" {},
isotropic_hardening : "Linear" {R0 : "0"}
}
};
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{ stress_potential : "Hooke" { young_modulus : 150.e3, poisson_ratio : 0.3 }, inelastic_flow : "Plastic" { criterion : "MohrCoulomb" { c : 3.e1, // cohesion phi : 0.523598775598299, // friction angle angle lodeT : 0.506145483078356, // transition angle (Abbo and Sloan) a : 1e1 // tension cuff-off parameter }, isotropic_hardening : "Linear" {R0 : "0"} } };However, this declaration will lead to a porosity growth:
@Brick StandardElastoViscoPlasticity{ stress_potential : "Hooke" { young_modulus : 150.e3, poisson_ratio : 0.3 }, inelastic_flow : "Plastic" { criterion : "MohrCoulomb" { c : 3.e1, // cohesion phi : 0.523598775598299, // friction angle angle lodeT : 0.506145483078356, // transition angle (Abbo and Sloan) a : 1e1 // tension cuff-off parameter }, isotropic_hardening : "Linear" {R0 : "0"} } porosity_evolution: {} };
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, as follows:
@ExternalStateVariable real ss;
ss.setGlossaryName("SolidSwelling");
@Brick StandardElastoViscoPlasticity{
stress_potential : "Hooke" {
young_modulus : 150.e3,
poisson_ratio : 0.3
},
inelastic_flow : "Plastic" {
criterion : "Gurson1975" {},
isotropic_hardening : "Linear" {R0 : "0"}
}
};
@Integrator{
// adding the contribution of the solid swelling
// the porosity increase
ff -= ;
}
By default, the value returned by the getPorosityEffectOnFlowRule
method of the stress criterion is used to decipher if the porosity
affects the flow rule.
However, this may be changed by specifying the
porosity_effect_on_flow_rule 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 porosition_evolution section as
follows:
porosity_evolution: {
save_individual_porosity_increase: true
}
However, this behaviour can be overloaded in every inelastic flow
sections and every nucleation models by setting the
save_porosity_increase boolean variable.
For example, one may use the following declaration:
inelastic_flow : "Plasticity" {
criterion : "Gurson1975",
isotropic_hardening : "Voce" {R0 : 200, Rinf : 100, b : 20},
save_porosity_increase : true
}
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.
[ \Delta,f = \paren{1-\bts{f}-\theta_{f},\Delta,f},\Delta,p,\trace{n} ]
[ \Delta,f=\Frac{\paren{1-\bts{f}},\Delta,p,\trace{n}}{1+\theta_{f}\Delta,p\trace{n}} ]
The porosity increment (\Delta,f) satisfies the following second order equation:
[ \begin{aligned} \Delta f &= \paren{1-\bts{f}-\theta_{f},\Delta f}^{2},\Delta p,\trace{n}\ &= \paren{1-\bts{f}}^{2},\Delta p,\trace{n}-2,\theta_{f},\Delta f,\paren{1-\bts{f}},\Delta p,\trace{n}+\theta_{f}^{2},\Delta f^{2},\Delta p,\trace{n} \end{aligned} ]
Equivalently,
[ A_{f},\Delta f^{2}+B_{f},\Delta f+C_{f}=0 ]
with:
- (A_{f}=-\theta_{f}^{2},\Delta p,\trace{n})
- (B_{f}=\paren{1+2,\theta_{f},,\paren{1-\bts{f}},\Delta p,\trace{n}},)
- (C_{f}=-\paren{1-\bts{f}}^{2},\Delta p,\trace{n})
@UpdateAuxiliaryStateVariables {
// equation associated with the porosity evolution
if (dp * trace(n) > 1.e-12) {
const auto theta_f = 1;
const auto A = -power<2>(theta_f) * dp * trace(n);
const auto B = (1+ 2 * theta_f * (1 - f) * dp * trace(n));
const auto C = -power<2>(1 - f) * dp * trace(n);
const auto Delta = B * B - 4 * A * C;
if ((Delta < 0) || (A == 0)) {
// neglecting second order terms
const auto df = -C / B;
f += df;
} else {
const auto df = (-B + sqrt(Delta)) / (2 * A);
f += df;
}
f = min(max(f, real(0)), real(1));
}
}
