5. Digital implementation

For the integration of the law of behavior in*Code_Aster*, we placed ourselves within the framework of the implicit integration of the laws of behavior.

5.1. Assessment of damage

We note \({\mathrm{\varepsilon }}^{-}\) the deformation and \({B}^{-}\) and \({d}^{-}\) the damage variables at the end of the previous time step (after convergence). We want to determine the evolution of the damage when applying a deformation increment \(\mathrm{\Delta }\mathrm{\varepsilon }\). The final deformation \(\varepsilon ={\varepsilon }^{-}+\Delta \varepsilon\) is therefore fixed.

The damage criterion is evaluated:

\(f(\mathrm{\varepsilon },{B}^{-},{d}^{-})=\sqrt{\mathrm{\alpha }{\parallel {F}_{-}^{B}\parallel }^{2}+(1-\mathrm{\alpha }){({F}_{+}^{d})}^{2}}-K(\mathrm{\varepsilon })\) eq 5.1-1

If \(f\le 0\), the damage variables do not change and we can move on to the next iteration for mechanical balance.

If \(f>0\), the damage variables change, respecting both the criterion and the normal flow. So we’re looking for \((\mathrm{\Delta }B,\mathrm{\Delta \; d},\mathrm{\Delta \; \gamma })\) system solution:

\(R(\Delta B,\Delta d,\Delta \gamma )=(\begin{array}{c}-\Delta B+\Delta \gamma \alpha {F}_{-}^{B}({B}^{-}+\Delta B)\\ -\Delta d+\Delta \gamma (1-\alpha ){F}_{+}^{d}({d}^{-}+\Delta d)\\ f(\varepsilon ,{B}^{-}+\Delta B,{d}^{-}+\Delta d)\end{array})=(\begin{array}{c}0\\ 0\\ 0\end{array})\) eq 5.1-2

This system is non-linear and requires the use of an iterative method. We chose to use a Newton-Raphson method. Clue \(n\) represents Newton’s iterations. Let \({R}^{n}(\mathrm{\Delta }{B}^{n},{\mathrm{\Delta \; d}}^{n},{\mathrm{\Delta \; \gamma }}^{n})\) be the residue of the system at iteration \(n\), the linearization of the system is written:

\({R}^{n+1}(\Delta {B}^{n+1},\Delta {d}^{n+1},\Delta {\gamma }^{n+1})={R}^{n}(\Delta {B}^{n},\Delta {d}^{n},\Delta {\gamma }^{n})+{({R}^{n}(\Delta {B}^{n},\Delta {d}^{n},\Delta {\gamma }^{n}))}^{\prime }:(\begin{array}{c}\Delta {B}^{n+1}-\Delta {B}^{n}\\ \Delta {d}^{n+1}-\Delta {d}^{n}\\ \Delta {\gamma }^{n+1}-\Delta {\gamma }^{n}\end{array})\) eq 5.1-3

Where \({R}^{n}=(\begin{array}{c}-\mathrm{\Delta }{B}^{n}+{\mathrm{\Delta \; \gamma }}^{n}\mathrm{\alpha }{F}_{-}^{Bn}\\ -{\mathrm{\Delta \; d}}^{n}+{\mathrm{\Delta \; \gamma }}^{n}(1-\mathrm{\alpha }){F}^{dn}\\ f(\mathrm{\varepsilon },{B}^{-}+\mathrm{\Delta }{B}^{n},{d}^{-}+{\mathrm{\Delta \; d}}^{n})\end{array})\)

and \({R}^{{n}^{\prime }}=\left[\begin{array}{ccc}-I+\alpha \Delta {\gamma }^{n}\frac{\partial {F}_{-}^{B}}{\partial {F}^{B}}:\frac{\partial {F}^{B}}{\partial B}& 0& \alpha {F}_{-}^{B}\\ 0& -1+(1-\alpha )\Delta {\gamma }^{n}\frac{\partial {F}_{+}^{d}}{\partial {F}^{d}}\frac{\partial {F}^{d}}{\partial d}& (1-\alpha ){F}_{+}^{d}\\ \frac{\alpha {F}_{-}^{B}}{\sqrt{\alpha {\parallel {F}_{-}^{B}\parallel }^{2}+(1-\alpha ){\mid {F}_{+}^{d}\mid }^{2}}}:\frac{\partial {F}^{B}}{\partial B}& \frac{(1-\alpha ){F}_{+}^{d}}{\sqrt{\alpha {\parallel {F}_{-}^{B}\parallel }^{2}+(1-\alpha ){\mid {F}_{+}^{d}\mid }^{2}}}\frac{\partial {F}^{d}}{\partial d}& 0\end{array}\right]\).

The linearized system \({R}^{n+1}(\Delta {B}^{n+1},\Delta {d}^{n+1},\Delta {\gamma }^{n+1})=0\) is then solved. We get \((\Delta {B}^{n+1},\Delta {d}^{n+1},\Delta {\gamma }^{n+1})\). This procedure is repeated until the residue is less than a convergence parameter.

5.2. Calculation of the tangent matrix

The tangent matrix is the 4th order tensor \(M\) defined by:

\(\delta {\sigma }_{\text{ij}}={M}_{\text{ijkl}}\delta {\varepsilon }_{\text{kl}}\iff {M}_{\text{ijkl}}=\frac{\partial {\sigma }_{\text{ij}}}{\partial \Delta {\varepsilon }_{\text{kl}}}\) eq 5.2-1

This matrix is not calculated on the continuous problem but on the incremental problem. We are therefore looking for the effects of a variation in the deformation increment between two successive time steps on the final stress variation, taking into account the fact that the internal variables can also change. At first order, we have:

\({\text{\delta \sigma }}_{\text{ij}}(\mathrm{\varepsilon },B,d)=\frac{\partial {\mathrm{\sigma }}_{\text{ij}}}{\partial {\mathrm{\Delta \; \varepsilon }}_{\text{kl}}}{\mid }_{B,d}{\mathrm{\delta \; \Delta \; \varepsilon }}_{\text{kl}}+\frac{\partial {\mathrm{\sigma }}_{\text{ij}}}{\partial {\mathrm{\Delta \; B}}_{\text{kl}}}{\mid }_{\mathrm{\varepsilon },d}{\mathrm{\delta \; \Delta \; B}}_{\text{kl}}+\frac{\partial {\mathrm{\sigma }}_{\text{ij}}}{\partial \mathrm{\Delta \; d}}{\mid }_{\mathrm{\varepsilon },B}\mathrm{\delta \; \Delta \; d}\) eq 5.2-2

either:

\({\text{\delta \sigma }}_{\text{ij}}(\mathrm{\varepsilon },B,d)=\left[\frac{\partial {\mathrm{\sigma }}_{\text{ij}}}{\partial {\mathrm{\Delta \; \varepsilon }}_{\text{kl}}}{\mid }_{B,d}+\frac{\partial {\mathrm{\sigma }}_{\text{ij}}}{\partial {\mathrm{\Delta \; B}}_{\text{mn}}}{\mid }_{\mathrm{\varepsilon },d}\frac{\partial {\mathrm{\Delta \; B}}_{\text{mn}}}{\partial {\mathrm{\Delta \; \varepsilon }}_{\text{kl}}}+\frac{\partial {\mathrm{\sigma }}_{\text{ij}}}{\partial \mathrm{\Delta \; d}}{\mid }_{\mathrm{\varepsilon },B}\frac{\partial \mathrm{\Delta \; d}}{\partial {\mathrm{\Delta \; \varepsilon }}_{\text{kl}}}\right]{\mathrm{\delta \; \Delta \; \varepsilon }}_{\text{kl}}\) eq 5.2-3

The \(\frac{\partial A}{\partial B}{\mid }_{C,D}\) notation used here means that we derive \(A\) from \(B\) for constant \(C\) and \(D\).

The tangent matrix is thus divided into two parts, one with constant damage, and the other representing the evolution of the internal variables:

\(M={M}^{\text{cst}}+{M}^{\text{evol}}\) eq 5.2-4

with \({M}_{\text{ijkl}}^{\text{cst}}=\frac{\partial {\mathrm{\sigma }}_{\text{ij}}}{\partial {\mathrm{\Delta \; \varepsilon }}_{\text{kl}}}{\mid }_{B,d}\) and \({M}_{\text{ijkl}}^{\text{evol}}=\frac{\partial {\mathrm{\sigma }}_{\text{ij}}}{\partial {\mathrm{\Delta \; B}}_{\text{mn}}}{\mid }_{\mathrm{\varepsilon },d}\frac{\partial {\mathrm{\Delta \; B}}_{\text{mn}}}{\partial {\mathrm{\Delta \; \varepsilon }}_{\text{kl}}}+\frac{\partial {\mathrm{\sigma }}_{\text{ij}}}{\partial \mathrm{\Delta \; d}}{\mid }_{\mathrm{\varepsilon },B}\frac{\partial \mathrm{\Delta \; d}}{\partial {\mathrm{\Delta \; \varepsilon }}_{\text{kl}}}\)

5.2.1. Constant damage term

Matrix \({M}^{\text{cst}}\) is the derivative of stress in relation to deformation at constant damage:

\(\begin{array}{}{M}_{\text{ijkl}}^{\text{cst}}=\mathrm{\lambda \; H}(\text{tr}\text{B \varepsilon }){B}_{\text{ij}}{B}_{\text{kl}}+\mathrm{\lambda \; f}(d)H(-\text{tr}\mathrm{\varepsilon }){\mathrm{\delta }}_{\text{ij}}{\mathrm{\delta }}_{\text{kl}}\\ +\frac{\mathrm{\mu }}{4}\frac{\partial {A}_{+\text{ip}}}{\partial {A}_{\text{mn}}}({B}_{\text{mk}}{\mathrm{\delta }}_{\text{nl}}+{\mathrm{\delta }}_{\text{mk}}{B}_{\text{nl}}+{B}_{\text{ml}}{\mathrm{\delta }}_{\text{nk}}+{\mathrm{\delta }}_{\text{ml}}{B}_{\text{nk}}){B}_{\text{pj}}\\ +\frac{\mathrm{\mu }}{4}{B}_{\text{ip}}\frac{\partial {A}_{+\text{pj}}}{\partial {A}_{\text{mn}}}({B}_{\text{mk}}{\mathrm{\delta }}_{\text{nl}}+{\mathrm{\delta }}_{\text{mk}}{B}_{\text{nl}}+{B}_{\text{ml}}{\mathrm{\delta }}_{\text{nk}}+{\mathrm{\delta }}_{\text{ml}}{B}_{\text{nk}})\\ +\mathrm{2\; \mu \; f}(d)\frac{\partial {\mathrm{\varepsilon }}_{-\text{ij}}}{\partial {\mathrm{\varepsilon }}_{\text{kl}}}\end{array}\) eq 5.2.1-1

where \(A=\mathrm{B\; \varepsilon }+\mathrm{\varepsilon \; B}\), \(H\) is the Heaviside function and the derivatives \(\frac{\partial {\mathrm{\varepsilon }}_{-}}{\partial \mathrm{\varepsilon }}\) and \(\frac{\partial {A}_{+}}{\partial A}\) are defined in Appendix 1.

5.2.2. Term related to the evolution of damage

Matrix \({M}^{\text{evol}}\) is calculated by differentiating the criterion and the law of normality. This derivation is based on the incremental problem and not on the continuous problem.

First of all, we differentiate the criterion:

\(f({F}^{B},{F}^{d})=\sqrt{\mathrm{\alpha }{\parallel {F}_{-}^{B}\parallel }^{2}+(1-\mathrm{\alpha }){({F}_{+}^{d})}^{2}}-K(\mathrm{\varepsilon })\le 0\underset{\text{diff}}{\Rightarrow }\frac{\mathrm{\alpha }{F}_{-}^{B}:\mathrm{\delta }{F}^{B}+(1-\mathrm{\alpha }){F}_{+}^{d}{\mathrm{\delta \; F}}^{d}}{\sqrt{\mathrm{\alpha }{\parallel {F}_{-}^{B}\parallel }^{2}+(1-\mathrm{\alpha }){({F}_{+}^{d})}^{2}}}-\mathrm{\delta \; K}=0\) eq 5.2.2-1

We then differentiate the flow law from the incremental problem, discretized in an implicit manner:

\(\begin{array}{}\begin{array}{}\Delta B=\Delta \gamma \alpha {F}_{-}^{B}\\ \Delta d=\Delta \gamma (1-\alpha ){F}_{+}^{d}\end{array}\}\Rightarrow (1-\alpha ){F}_{+}^{d}\Delta B=\alpha {F}_{-}^{B}\Delta d\\ \\ \underset{\text{diff}}{\Rightarrow }(1-\alpha )\delta {F}_{+}^{d}\Delta B+(1-\alpha ){F}_{+}^{d}\delta \Delta B=\alpha \delta {F}_{-}^{B}\Delta d+\alpha {F}_{-}^{B}\delta \Delta d\end{array}\) eq 5.2.2-2

We’re looking for the relationship between \(\delta \Deltab\) , \(\delta \Deltad\) and \(\delta \Delta \varepsilon\). Variations in thermodynamic forces can be expressed as a function of variations in deformations and damage variables:

\(\begin{array}{}\delta {F}_{-}^{B}=\frac{\partial {F}_{-}^{B}}{\partial {F}^{B}}:\left[\frac{\partial {F}^{B}}{\partial \Delta \varepsilon }:\delta \Delta \varepsilon +\frac{\partial {F}^{B}}{\partial \Delta B}:\delta \Delta B\right]\\ \delta {F}^{d}=\frac{\partial {F}_{+}^{d}}{\partial {F}^{d}}\left[\frac{\partial {F}^{d}}{\partial \Delta \varepsilon }:\delta \Delta \varepsilon +\frac{\partial {F}^{d}}{\partial \Delta d}\delta \Delta d\right]\end{array}\) eq 5.2.2-3

The system of equations defined by [éq 5.2.2-1], [éq 5.2.2-2], and [éq 5.2.2-3] results in the following expressions:

\(\begin{array}{ccc}\delta \Delta \text{B}& \text{=}& {\chi }^{-1}:\xi :\delta \Delta \varepsilon \\ \delta \Delta d& \text{=}& -\left[\frac{\tau }{\rho }+\frac{\theta }{\rho }:{\chi }^{-1}:\xi \right]:\delta \Delta \varepsilon \end{array}\) eq 5.2.2-4

with

\({\mathrm{\khi }}_{\text{ijkl}}={\left[\frac{\partial {F}_{+}^{d}}{\partial {F}^{d}}\frac{\mathrm{\Delta \; B}}{{F}_{+}^{d}}-\frac{{\mathrm{\alpha \; F}}_{-}^{B}}{(1-\mathrm{\alpha }){F}^{d}\frac{\partial {F}^{d}}{\partial d}}\right]}_{\text{ij}}{\left[{\mathrm{\alpha \; F}}_{-}^{B}:\frac{\partial {F}_{-}^{B}}{\partial {F}^{B}}:\frac{\partial {F}^{B}}{\partial B}\right]}_{\text{kl}}+\mathrm{\alpha \; \Delta \; d}{\left[\frac{\partial {F}_{-}^{B}}{\partial {F}^{B}}:\frac{\partial {F}^{B}}{\partial B}\right]}_{\text{ijkl}}-\frac{1}{2}(1-\mathrm{\alpha }){F}_{+}^{d}({\mathrm{\delta }}_{\text{ik}}{\mathrm{\delta }}_{\text{jl}}+{\mathrm{\delta }}_{\text{il}}{\mathrm{\delta }}_{\text{jk}})\)

\(\begin{array}{}{\mathrm{\xi }}_{\text{ijkl}}={\left[\frac{{\mathrm{\alpha \; F}}_{-}^{B}}{(1-\mathrm{\alpha }){F}^{d}\frac{\partial {F}^{d}}{\partial d}}-\frac{\partial {F}_{+}^{d}}{\partial {F}^{d}}\frac{\mathrm{\Delta \; B}}{{F}_{+}^{d}}\right]}_{\text{ij}}{\left[{\mathrm{\alpha \; F}}_{-}^{B}:\frac{\partial {F}_{-}^{B}}{\partial {F}^{B}}:\frac{\partial {F}^{B}}{\partial \mathrm{\varepsilon }}+(1-\mathrm{\alpha }){F}_{+}^{d}\frac{\partial {F}^{d}}{\partial \mathrm{\varepsilon }}-\sqrt{\mathrm{\alpha }{\parallel {F}_{-}^{B}\parallel }^{2}+(1-\mathrm{\alpha }){({F}_{+}^{d})}^{2}}\frac{\partial K}{\partial \mathrm{\varepsilon }}\right]}_{\text{kl}}\\ +(1-\mathrm{\alpha })\frac{\partial {F}_{+}^{d}}{\partial {F}^{d}}{\mathrm{\Delta \; B}}_{\text{ij}}\frac{\partial {F}^{d}}{\partial {\mathrm{\varepsilon }}_{\text{kl}}}-\mathrm{\alpha \; \Delta \; d}{\left[\frac{\partial {F}_{-}^{B}}{\partial {F}^{B}}:\frac{\partial {F}^{B}}{\partial \mathrm{\varepsilon }}\right]}_{\text{ijkl}}\end{array}\)

\({\mathrm{\tau }}_{\text{ij}}={\left[{\mathrm{\alpha \; F}}_{-}^{B}:\frac{\partial {F}_{-}^{B}}{\partial {F}^{B}}:\frac{\partial {F}^{B}}{\partial \mathrm{\varepsilon }}+(1-\mathrm{\alpha }){F}_{+}^{d}\frac{\partial {F}^{d}}{\partial \mathrm{\varepsilon }}-\sqrt{\mathrm{\alpha }{\parallel {F}_{-}^{B}\parallel }^{2}+(1-\mathrm{\alpha }){({F}_{+}^{d})}^{2}}\frac{\partial K}{\partial \mathrm{\varepsilon }}\right]}_{\text{ij}}\)

\({\mathrm{\theta }}_{\text{ij}}={\left[{\mathrm{\alpha \; F}}_{-}^{B}:\frac{\partial {F}_{-}^{B}}{\partial {F}^{B}}:\frac{\partial {F}^{B}}{\partial B}\right]}_{\text{ij}}\)

\(\rho =(1-\alpha ){F}_{\text{+}}^{d}\frac{\partial {F}^{d}}{\partial d}\)

\({\chi }^{-1}\) such as \({\chi }_{\text{ijkl}}^{-1}{\chi }_{\text{klmn}}=\frac{1}{2}({\delta }_{\text{im}}{\delta }_{\text{jn}}+{\delta }_{\text{in}}{\delta }_{\text{jm}})\)

So we get the expressions for \(\frac{\partial \Delta B}{\partial \Delta \varepsilon }\) and \(\frac{\partial \Delta d}{\partial \Delta \varepsilon }\).

In addition, according to the definition of stress and thermodynamic forces, which derive from the same energy, we have:

\(\frac{\partial \mathrm{\sigma }}{\partial \mathrm{\Delta \; B}}{\mid }_{\mathrm{\varepsilon },d}=\frac{\partial {F}^{B}}{\partial \mathrm{\varepsilon }}{\mid }_{B}\) and \(\frac{\partial \mathrm{\sigma }}{\partial \mathrm{\Delta \; d}}{\mid }_{\mathrm{\varepsilon },B}=\frac{\partial {F}^{d}}{\partial \mathrm{\varepsilon }}{\mid }_{d}\) eq 5.2.2-5

which allows us to calculate the part of the tangent matrix relating to the evolution of the damage:

\({M}_{\text{ijkl}}^{\text{evol}}=\frac{\partial {\mathrm{\sigma }}_{\text{ij}}}{\partial {\mathrm{\Delta \; B}}_{\text{mn}}}{\mid }_{\mathrm{\varepsilon },d}\frac{\partial {\mathrm{\Delta \; B}}_{\text{mn}}}{\partial {\mathrm{\Delta \; \varepsilon }}_{\text{kl}}}+\frac{\partial {\mathrm{\sigma }}_{\text{ij}}}{\partial \mathrm{\Delta \; d}}{\mid }_{\mathrm{\varepsilon },B}\frac{\partial \mathrm{\Delta \; d}}{\partial {\mathrm{\Delta \; \varepsilon }}_{\text{kl}}}\) eq 5.2.2-6

Note:

Note that the tangent matrix is not symmetric. This is due to the fact that the damage threshold depends on the deformation. In the case of a constant threshold, the matrix is indeed symmetric, since the model is then very standard generalized and an implicit integration scheme is used.