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.