4. Numerical formulation#

4.1. Keywords, material data and internal variables#

For foreseeable applications, the model was implemented under two distinct keywords: “ROUSS_PR” for the plastic model with cavity nucleation or “ROUSS_VISC” for the viscoplastic model without nucleation. This makes it possible to avoid unnecessary numerical calculations. The corresponding simplified equations are obtained from the general equations by setting \({\mathrm{\sigma }}_{0}=0\) or \({A}_{n}=0\) respectively.

All the parameters of the model are provided under the keywords factors” ROUSSELIER “or” ROUSSELIER_FO “or” “and” TRACTION “(to define the traction curve) of the control DEFI_MATERIAU ([U4.43.01]). The parameters of the viscoplastic model (\({\mathrm{\sigma }}_{0}\), \({\dot{\mathrm{\varepsilon }}}_{0}\), and \(m\)) are provided by the keyword “VISC_SINH”.

The internal variables produced in*Code_Aster* are:

  • V1, the cumulative plastic deformation \(p\),

  • V2, porosity \(f\),

  • V3 to V8, the \({\mathrm{\varepsilon }}^{e}\) elastic deformation tensor,

  • V9, the plasticity indicator (0 if the last calculated increment is elastic, 1 if a regular plastic solution, 2 if a singular plastic solution).

We now present the numerical integration of the law of behavior and give the expression for the tangent matrix (options FULL_MECA and RIGI_MECA_TANG).

4.2. Expression of the discretized model#

The numerical resolution is done by a \(\mathrm{\theta }\) -method, with \(0\mathrm{\le }\theta \mathrm{\le }1\), and incrementally. For all quantities \(Q\), we define:

\(Q={Q}^{\text{-}}+\Deltaq\)

\({Q}^{\theta }={Q}^{\text{-}}+\theta \Deltaq\)

Incremental writing requires taking into account the possible variation in material properties (due, for example, to a change in temperature during the time step).

The discretized system of equations is:

\({\tilde{\tau }}^{\theta }=2\mu \theta \Delta \tilde{{\varepsilon }^{e}}+\frac{2\mu \theta +(1-\theta )2{\mu }^{\text{-}}}{2{\mu }^{\text{-}}}\tilde{{\tau }^{\text{-}}}\text{=}2\mu \theta (\Delta \tilde{\varepsilon }-\Delta {\tilde{\varepsilon }}^{p})+\frac{2\mu \theta +(1-\theta )2{\mu }^{\text{-}}}{2{\mu }^{\text{-}}}\tilde{{\tau }^{\text{-}}}\) eq 4.2-1

\({\tau }_{m}^{\theta }=K\theta \text{tr}\Delta {\varepsilon }^{e}+\frac{\mathrm{3K}\ast \theta +(1-\theta ){\mathrm{3K}}^{\text{-}}}{{\mathrm{3K}}^{\text{-}}}{\tau }_{m}^{\text{-}}=K\theta (\text{tr}\Delta \varepsilon -\text{tr}\Delta {\varepsilon }^{p})\frac{\mathrm{3K}\ast \theta +(1-\theta ){\mathrm{3K}}^{\text{-}}}{{\mathrm{3K}}^{\text{-}}}{\tau }_{m}^{\text{-}}\) eq 4.2-2

\(\Delta {\tilde{\varepsilon }}^{p}=\Delta p\frac{3{\tilde{\tau }}^{\theta }}{2{\tau }_{\text{eq}}^{\theta }}\) eq 4.2-3

\(\text{tr}\Delta {\varepsilon }^{p}=\Delta {\text{pD}}_{1}({f}^{\theta }+{A}_{n}{p}^{\theta })\text{exp}(\frac{{\tau }_{m}^{\theta }}{{\sigma }_{1}})\) eq 4.2-4

\(\Deltaf ={A}_{1}(1-{f}^{\theta })\text{tr}\Delta {\varepsilon }^{p}\) eq 4.2-5

\({\Phi }_{\text{vp}}^{\theta }={\tau }_{\text{eq}}^{\theta }+{\sigma }_{1}{D}_{1}({f}^{\theta }+{A}_{n}{p}^{\theta })\text{exp}(\frac{{\tau }_{m}^{\theta }}{{\sigma }_{1}})-R({p}^{\theta })-{\sigma }_{0}{\text{sh}}^{\text{-}1}\left[{(\frac{\Delta p}{{\dot{\varepsilon }}_{0}\Delta t})}^{\frac{1}{m}}\right]=0\) eq 4.2-6

This system comes down to solving a single scalar equation for the unknown \(\Deltaf\), knowing \(\Delta \varepsilon\), \(\Deltat\) and the quantities \({Q}^{\text{-}}\). Note that \(\rho\) does not intervene in the algorithm, on the other hand it will intervene in the calculation of the coherent tangent matrix. The following are calculated successively:

\({\tau }_{m}^{\theta }=\frac{\mathrm{3K}\ast \theta +(1-\theta ){\mathrm{3K}}^{\text{-}}}{{\mathrm{3K}}^{\text{-}}}{\tau }_{{m}^{\text{-}}}+K\theta (\text{tr}\Delta \varepsilon -\frac{\Delta f}{{A}_{1}(1-{f}^{\theta })})\) eq 4.2-7

\(\Deltap\) is the positive root of the quadratic equation:

\({A}_{n}\theta {(\Delta p)}^{2}+({f}^{\theta }+{A}_{n}{p}^{\text{\_}})\Delta p-\frac{\Delta f}{{A}_{1}(1-{f}^{\theta })}\frac{1}{{D}_{1}\text{exp}({\tau }_{m}^{\theta }/{\sigma }_{1})}=0\) eq 4.2-8

\({\tilde{\tau }}^{\theta }=(1-\frac{3\mu \theta \Delta p}{{\left[\frac{2\mu \theta +(1-\theta )2{\mu }^{\text{-}}}{2{\mu }^{\text{-}}}{\tilde{\tau }}^{\text{-}}+2\mu \theta \Delta \tilde{\varepsilon }\right]}_{\text{eq}}})(\frac{2\mu \theta +(1-\theta )2{\mu }^{\text{-}}}{2{\mu }^{\text{-}}}{\tilde{\tau }}^{\text{-}}+2\mu \theta \Delta \tilde{\varepsilon })\) eq 4.2-9

\({\tau }_{\text{eq}}^{\theta }={\left[\frac{2\mu \theta +(1-\theta )2{\mu }^{\text{-}}}{2{\mu }^{\text{-}}}{\tilde{\tau }}^{\text{-}}+2\mu \theta \Delta \tilde{\varepsilon }\right]}_{\text{eq}}-3\mu \theta \Delta p\) eq 4.2-10

The scalar equation for \(\Deltaf\) is the equation [éq 4.2-6] \({\Phi }_{\text{vp}}^{\theta }=0\).

Note 1:

As \(\Deltaf\) is very weak in most of the structure, it would be best to use \(\Deltap\) as the main unknown. But in this case it is not possible to reduce yourself to a scalar equation, which makes it more difficult to use a Newtons-type method. This is also one reason why the equations [éq 1-1] , [éq 3.2-6], and [éq4.2-5] have not been modified by the introduction of cavity nucleation.

Note 2:

The equation [éq 3.2-6] can be integrated exactly:

\(\text{tr}{\varepsilon }^{p}=\frac{1}{{A}_{1}}\text{ln}(\frac{1-{f}_{0}}{1-f})\)

from where:

\(\text{tr}\Delta {\varepsilon }^{p}=\frac{1}{{A}_{1}\theta }\text{ln}(\frac{1-{f}^{\text{-}}}{1-{f}^{\theta }})\)

As the numerical parameter \({A}_{1}\) can be changed discontinuously, the derived form [éq 4.2-5] has been retained, including in the calculation of the coherent tangent matrix. If the use of the parameter \({A}_{1}\) were to be abandoned in a later release, consideration should be given to using the built-in form.

Note 3:

The integrated form \({\Phi }_{\text{vp}}^{\theta }=0\) is used, including in plasticity instead of the consistency relationship \(\dot{F}=0\) which gives \(\dot{p}\) . The coherent tangent matrix is calculated with this integrated form.

4.3. Solving the nonlinear scalar equation#

Equation \({\Phi }_{\text{vp}}^{\theta }(\Deltaf )=0\) is solved by a Newton algorithm with controlled terminals in routine LCROUS. \({\Phi }_{\text{vp}}^{\theta }(\Deltaf )\) and its derivative with respect to \(\Deltaf\) are calculated in the routine RSLPHI called by LCROUS. The initial values of the terminals are:

  • lower bound: \(\Delta {f}_{1}=0\) since \({\Phi }_{\text{vp}}^{\theta }(0)<0\) (we checked beforehand that the elastic branch (negative threshold) is not a solution),

  • upper bound: \(\Delta {f}_{2}\) such as \({\Phi }_{\text{vp}}^{\theta }(0)>0\) searched by dichotomy between 0 and \(1-{f}^{\text{-}}\) (first value for this search: \(\frac{1-{f}^{\text{-}}}{2}\)).

Newton’s algorithm starts with the value \(\Deltaf =0\). Whatever the value found for \(\Deltaf\), we therefore note that the function \({\Phi }_{\text{vp}}^{\theta }(\Deltaf )\) and its derivative with respect to \(\Deltaf\) are at least calculated for \(\Deltaf =0\) and \(\frac{1-{f}^{\text{-}}}{2}\).

Developments made to improve the convergence and robustness of the algorithm are described in [bib5].

4.4. Behavior tangent matrix expression#

One gives here the expression for the tangent matrix (option FULL_MECA during Newton iterations, option RIGI_MECA_TANG for the first iteration).

For option RIGI_MECA_TANG, the tangent operator is the same operator that connects \({\varepsilon }^{e}\) to \(\sigma\) in [éq3.2-2].

For option FULL_MECA, the tangent matrix is obtained by linearizing the system of equations that govern the law of behavior: [éq 4.2-1] to [éq 4.2-6]. So it is a coherent tangent matrix.

To simplify the expressions, in this paragraph we note [§4.5]: \(Q\) for \({Q}^{\mathrm{\theta }}\), the quantities all being expressed at time \({t}^{\theta }={t}^{\text{-}}+\theta \Deltat\). The coherent tangent matrix is:

\(\frac{\delta \sigma }{\delta \varepsilon }=\rho \left[{a}_{3}\mathrm{II}+\mathrm{Id}\otimes (\frac{{a}_{1}-{a}_{3}}{3}\mathrm{Id}+{a}_{2}\tilde{\tau })+\tilde{\tau }\otimes ({a}_{4}\tilde{\tau }+\frac{{a}_{5}}{3}\mathrm{Id})+\tau \otimes ({y}_{4}(\frac{{a}_{1}}{\mathrm{3K}}-1)\mathrm{Id}+\frac{{y}_{5}}{\text{K}}\tilde{\tau })\right]\) eq 4.4-1

This operator is calculated in routine RSLJPL. The coefficients are calculated as follows:

\({a}_{1}=\mathrm{3K}+{y}_{1}K{\tau }_{\text{eq}}({z}_{7}+{z}_{2}\theta \Deltap )\) eq 4.4-2

\({a}_{2}=\mathrm{\mu }({y}_{1}+{y}_{3}){\mathrm{\sigma }}_{1}\) eq 4.4-3

\({a}_{3}=\frac{2\mu {\tau }_{\text{eq}}}{{z}_{5}}\) eq 4.4-4

\({a}_{4}=3\mu {y}_{2}{x}_{2}\) eq 4.4-5

\({a}_{5}=3\mu {y}_{1}{\sigma }_{1}\) eq 4.4-6

\({a}_{6}=3\muk \theta \Deltap -{a}_{2}{\tau }_{\text{eq}}{\sigma }_{1}\) eq 4.4-7

\({y}_{1}=-\frac{3{\text{Kz}}_{6}{z}_{1}(f+{A}_{n}p)}{{x}_{1}{\tau }_{\text{eq}}}\) eq 4.4-8

\({y}_{2}=-\frac{\mathrm{3\mu }}{{x}_{1}{z}_{5}{\mathrm{\tau }}_{\text{eq}}^{2}}\) eq 4.4-9

\({y}_{3}=-\frac{3{\text{Kz}}_{6}{z}_{1}{A}_{n}\theta \Deltap }{{x}_{1}{\tau }_{\text{eq}}}\) eq 4.4-10

\({y}_{4}=\frac{{A}_{1}{z}_{8}}{{z}_{1}}+\frac{{z}_{9}{\sigma }_{1}}{{z}_{7}+{z}_{2}\theta \Deltap }\) eq 4.4-11

\({y}_{5}=\frac{{A}_{1}{a}_{2}{z}_{8}}{{z}_{1}}-\frac{{z}_{9}{a}_{6}}{{\tau }_{\text{eq}}({z}_{7}+{z}_{2}\theta \Deltap )}\) eq 4.4-12

\({z}_{1}=1+{A}_{1}\theta \Delta {\text{pD}}_{1}(f+{A}_{n}p)\text{exp}(\frac{{\tau }_{m}}{{\sigma }_{1}})\) eq 4.4-13

\({z}_{2}=\mathrm{3\mu }+{R}_{\text{vp}}\) eq 4.4-14

\({z}_{3}=K(f+{A}_{n}p){z}_{1}-{A}_{1}{\sigma }_{1}(1-f)\) eq 4.4-15

\({z}_{4}={R}_{\text{vp}}\theta \Deltap -{\tau }_{\text{eq}}\) eq 4.4-16

\({z}_{5}={\tau }_{\text{eq}}+3\mu \theta \Deltap\) eq 4.4-17

\({z}_{6}={D}_{1}\text{exp}(\frac{{\tau }_{m}}{{\sigma }_{1}})\) eq 4.4-18

\({z}_{7}={z}_{6}{\sigma }_{1}(f+{A}_{n}p)\) eq 4.4-19

\({z}_{8}=\frac{1-f}{1-f-{A}_{n}p}\) eq 4.4-20

\({z}_{9}=\frac{{A}_{n}}{1-f-{A}_{n}p}\) eq 4.4-21

\({x}_{1}={z}_{3}{z}_{6}({z}_{7}+{z}_{2}\theta \Deltap )+{z}_{1}{z}_{2}{\sigma }_{1}-{x}_{3}\) eq 4.4-22

\({x}_{2}=-{z}_{3}{z}_{6}\theta \Deltap ({z}_{4}+{z}_{7})-{z}_{1}{z}_{4}{\sigma }_{1}+{x}_{3}\theta \Deltap\) eq 4.4-23

\({x}_{3}={A}_{n}{z}_{1}{z}_{6}{\mathrm{\sigma }}_{1}^{2}\) eq 4.4-24

\({R}_{\text{vp}}=\frac{\text{dR}(p)}{\text{dp}}+\frac{1}{\theta \Deltat }\frac{\text{dS}(\Deltap /\Deltat )}{d\dot{p}}\) eq 4.4-25

For the plastic model with nucleation of cavities” ROUSSELIER_PR “and for the viscoplastic model without nucleation” ROUSSELIER_VISC “, the corresponding simplified equations are obtained from the equations above by setting \({R}_{\text{vp}}=\text{dR}(p)/\text{dp}\) and \({A}_{n}=0\) respectively.