4. Digital resolution#
4.1. Model integration#
The integration of the model is carried out in two steps, first of all an elastic prediction of the stress:
1st stage: elasticity \(\{\begin{array}{c}\dot{\epsilon }\text{=}{\dot{\epsilon }}^{e}\text{+}{\dot{\epsilon }}^{p}\\ \dot{\sigma }\text{=}E{\dot{\epsilon }}^{e}\\ {\dot{\epsilon }}^{p}=0\\ \dot{p}=0\\ \dot{D}=0\end{array}\) eq 4.1-1
Note in particular that damage \(D\) does not intervene in the elastic stress-deformation relationship (no coefficient \((1\text{-}D)\) in the elasticity coefficients).
The second step consists in an implicit update of the internal variables of the model by corrections made with zero total deformation. In this step, we iterate between plasticity and damage.
\(\mathrm{\{}\begin{array}{c}\dot{\varepsilon }\text{=}{\dot{\varepsilon }}^{e}\text{+}{\dot{\varepsilon }}^{p}\text{=0}\\ \dot{\sigma }\text{=}E{\dot{\varepsilon }}^{p}\\ {\dot{\varepsilon }}^{p}\text{=}\dot{\lambda }\frac{\mathrm{\partial }f}{\mathrm{\partial }\sigma }\\ \dot{r}\text{=}\text{-}\dot{\lambda }\frac{\mathrm{\partial }f}{\mathrm{\partial }R}\end{array}\) eq 4.1-2
The expression for the plastic multiplier is obtained by linearizing the threshold function with respect to the internal variables and by requiring compliance with the convergence criterion:
\({\text{f=f}}^{\text{-}}\text{+}\frac{\partial {f}^{\text{-}}}{\partial \sigma }(\sigma \text{-}{\sigma }^{\text{-}})\text{+}\frac{\partial {f}^{\text{-}}}{\partial R}(R\text{-}{R}^{\text{-}})\text{=}0\) eq 4.1-3
Let us consider the discretization of relaxation equations:
\(\mathrm{\{}\begin{array}{c}\sigma \text{-}{\sigma }^{\text{-}}\text{=-}\Delta \lambda \mathrm{E}\frac{\mathrm{\partial }{f}^{\text{-}}}{\mathrm{\partial }\sigma }\\ R\text{-}{R}^{\text{-}}{\text{=R}}^{\text{'}}(r)\Delta r\\ \Delta \text{p=}\frac{\Delta \lambda }{1\text{-}{D}^{\text{-}}}\\ {R}^{\text{'}}(p)\text{=}\frac{K}{m}{\left[\frac{{\sigma }_{\text{eq}}}{(1\text{-}{D}^{\text{-}})K}\text{-}\frac{{\sigma }_{y}}{K}\right]}^{(1\text{-}m)}\end{array}\) eq 4.1-4
The expression for the plastic multiplier is obtained by replacing [éq 4.1-3] in [éq 4.1-4], (\(\Delta \text{r=}\Delta \lambda ` and for :math:`D\) set to \({D}^{\text{-}}\), \(\text{R'}(p)\text{=R'}(r)\))
\(\Delta \lambda \text{=}\frac{{f}^{\text{-}}}{\frac{\partial {f}^{\text{-}}}{\partial \sigma }E\frac{\partial {\phi }^{\text{-}}}{\partial \sigma }\text{-}\frac{\partial {f}^{\text{-}}}{\partial R}\frac{K}{m}{\left[\frac{{\sigma }_{\text{eq}}}{K(1\text{-}{D}^{\text{-}})}\text{-}\frac{{\sigma }_{y}}{K}\right]}^{(1\text{-}m)}}\) eq 4.1-5
The expressions of the various derivatives of the threshold surface and of the dissipation potential are given in the following equations:
\(\begin{array}{}\frac{\partial f}{\partial \sigma }\text{=}\frac{3}{2}\frac{\sigma \text{'}}{{\sigma }_{\text{eq}}(1\text{-}D)}\\ \frac{\partial \phi }{\partial \sigma }\text{=}\frac{3}{2}\frac{\sigma \text{'}}{{\sigma }_{\text{eq}}(1-D)}\\ \frac{\partial f}{\partial R}\text{=}\text{-}1\end{array}\) eq 4.1-6
In the 1D case:
\(\frac{\partial f}{\partial \sigma }\text{=}\frac{\partial \varphi }{\partial \sigma }\text{=}\frac{\text{sgn}(\sigma )}{1\text{-}D}\) eq 4.1-7b
The value obtained for the plastic multiplier at each integration is re-injected into the relaxation equations until convergence (return of the stress to the threshold surface). The criterion used consists in stopping when the threshold value at the next iteration has become sufficiently small compared to \({\sigma }_{y}\):
\(\mathrm{\mid }f\mathrm{\mid }\mathrm{\le }\text{tol}\text{.}{\sigma }_{y}\) eq 4.1-8
For direct use of the normality properties of threshold surfaces, the calculation to be performed at each iteration only comes down to the calculation of the plasticity multiplier allowing corrections if necessary (no Jacobian calculation as for Newton’s methods) [bib1], [bib6].
The next step is damage, which takes into account the variation in the damage variable. The damage threshold function is as follows:
\({f}_{\text{end}}\text{=p}-{p}_{D}\) eq 4.1- 9
if \({f}_{\text{end}}>0\) we calculate the damage then we go back to the plasticity correction:
\(\text{D=}\frac{{D}_{C}}{{p}_{R}\text{-}{p}_{D}}(p\text{-}{p}_{D})\) eq 4.1-10
Note: in practice, \(D\) evolves even if \({\text{p>p}}_{R}\), we limit \(D\) to 0.99 max. \({D}_{c}\) is defined by the D_ CORR keyword from DEFI_MATERIAU/CORR_ACIER
The same approach for convergence is adopted:
\(\mathrm{\mid }f\mathrm{\mid }\mathrm{\le }\text{tol}\text{.}{\sigma }_{y}\) eq 4.1- 11
Refer to [bib1] and [bib6] for all the details concerning this method and algorithm used.
4.2. Calculation of the tangent matrix#
The aim is to calculate the tangent matrices (continuous and consistent) for the plastic part and the damage part. The options RIGI_MECA_TANGet FULL_MECAde the tangent matrix are calculated for the 1D and 3D finite elements:
4.2.1. 1D finished element#
In the plastic field:
\(\delta {\varepsilon }^{p}\mathrm{=}\delta \lambda \mathrm{=}\frac{\delta {f}^{\text{-}}}{E+\frac{K}{m}{p}^{\frac{1}{m}\mathrm{-}1}}\mathrm{=}\frac{E\delta \varepsilon }{E+\frac{K}{m}{p}^{\frac{1}{m}\mathrm{-}1}}\) eq 4.2.1-1
\(\frac{\mathrm{\partial }\sigma }{\mathrm{\partial }\Delta \varepsilon }\mathrm{=}E(\delta \varepsilon \mathrm{-}\delta {\varepsilon }^{p})\mathrm{=}\frac{\frac{K}{m}{p}^{\frac{1}{m}\text{-}1}}{1\text{+}\frac{K}{\text{mE}}{p}^{\frac{1}{m}\text{-}1}}\) eq 4.2.1-2
In the area of damage:
\(\frac{\delta \sigma }{1\mathrm{-}D}+\frac{\sigma }{{(1\mathrm{-}D)}^{2}}\delta D\mathrm{-}R\text{'}(p)\delta p\mathrm{=}\mathrm{-}\delta {f}^{\text{-}}\)
\(\frac{\mathrm{-}E}{{(1\mathrm{-}D)}^{2}}\delta \lambda +\frac{\sigma }{{(1\mathrm{-}D)}^{2}}\frac{{D}_{c}}{{p}_{R}\mathrm{-}{p}_{D}}\frac{\delta \lambda }{(1\mathrm{-}D)}\mathrm{-}\frac{K}{m}{p}^{\frac{1}{m}\mathrm{-}1}\frac{\delta \lambda }{(1\mathrm{-}D)}\mathrm{=}\mathrm{-}\delta {f}^{\text{-}}\mathrm{=}\frac{E}{(1\mathrm{-}D)}\delta \varepsilon\)
\(\frac{\mathrm{\partial }\sigma }{\mathrm{\partial }\Delta \varepsilon }\mathrm{=}E(\delta \varepsilon \mathrm{-}\delta {\varepsilon }^{p})\mathrm{=}E\frac{\frac{K}{m}(1\mathrm{-}D){p}^{\frac{1}{m}\text{-}1}\mathrm{-}\frac{\sigma }{(1\mathrm{-}D)}\frac{{D}_{c}}{{p}_{R}\mathrm{-}{p}_{D}}}{E+\frac{K}{m}(1\mathrm{-}D){p}^{\frac{1}{m}\text{-}1}\mathrm{-}\frac{\sigma }{(1\mathrm{-}D)}\frac{{D}_{c}}{{p}_{R}\mathrm{-}{p}_{D}}}\)
4.2.2. 3D finished element#
The absence of damage in the plastic part makes it possible to use the tangent matrix calculated for the elastoplastic behavior with linear and isotropic nonlinear kinematic work hardening [bib7]:
\(\frac{\mathrm{\partial }\sigma }{\Delta \varepsilon }\text{=}{\lambda }^{\text{*}}\overrightarrow{1}\otimes \overrightarrow{1}\text{+}2{\mu }^{\text{*}}\text{Id}\mathrm{-}\xi \frac{9{\mu }^{2}}{H(p)}(1\text{-}\frac{{R}^{\text{'}}(p)\text{.}\Delta p}{(R(p)\text{+}{\sigma }_{y})})\frac{1}{{R}^{\text{'}}(p)\text{+}3\mu }(\frac{{\sigma }^{\text{dev}}}{(R(p)\text{+}{\sigma }_{y})}\otimes \frac{{\sigma }^{\text{dev}}}{(R(p)\text{+}{\sigma }_{y})})\) eq 4.2.2-1
with \({\lambda }^{\text{*}}\text{=}K\text{-}\frac{2\mu }{\mathrm{3H}(\Delta p)}2{\mu }^{\text{*}}\text{=}\frac{2\mu }{H(\Delta p)}H(\Delta p)\text{=}1\text{+}\frac{3\mu \xi \text{.}\Delta p}{(R(p)\text{+}{\sigma }_{y})}\)
The initial tangent matrix, used by option RIGI_MECA_TANGest obtained by adopting the behavior of the previous step (elastic or plastic, signified by an internal variable \(\xi\) equal to 0 or 1) and by taking \(\Delta p=0\) into the equation eq 4.2.2-2.
Taking into account the presence of damage, the following tangent matrix is obtained:
\(\begin{array}{c}\frac{\mathrm{\partial }\sigma }{\mathrm{\partial }\Delta \varepsilon }\text{=}{\lambda }^{\text{*}}\overrightarrow{1}\otimes \overrightarrow{1}\text{+}2{\mu }^{\text{*}}\text{Id}\text{-}\xi \frac{9{\mu }^{2}}{H(p)}(1\text{-}\frac{\Delta p((1\text{-}D){R}^{\text{'}}(p)\text{-}R(p){D}^{\text{'}}(p))}{(1\text{-}D)(R(p)\text{-}{\sigma }_{y})})\frac{1}{(3\mu \text{+}(1\text{-}D){R}^{\text{'}}(p)\text{-}R(p){D}^{\text{'}}(p))}\\ (\frac{{\sigma }^{\text{dev}}}{(1\text{-}D)(R(p)\text{+}{\sigma }_{y})}\otimes \frac{{\sigma }^{\text{dev}}}{(1\text{-}D)(R(p)\text{+}{\sigma }_{y})})\end{array}\) eq 4.2.2-2
with \({\lambda }^{\text{*}}\text{=}K\text{-}\frac{2\mu }{\mathrm{3H}(\Delta p)}2{\mu }^{\text{*}}\text{=}\frac{2\mu }{H(\Delta p)}H(\Delta p)\text{=}1\text{+}\frac{3\mu \xi \text{.}\Delta p}{(1\text{-}D)(R(p)\text{-}{\sigma }_{y})}\)
for option FULL_MECA: \({\sigma }^{\mathrm{dev}}\text{=}\tilde{\sigma }\)
for option RIGI_MECA_TANG: \({\sigma }^{\mathrm{dev}}\text{=}\tilde{\sigma }\)
4.3. Stored internal variables#
In the following table, we indicate the internal variables stored at each Gauss point for the model of steel subjected to corrosion:
Internal variable |
Physical sense |
V1 |
|
V2 |
|
V3 |
Plasticity indicator (0 if elastic, 1 if laminated, i.e. as soon as p is not zero) |