4. Implementation in Code_Aster

The use of model HAYHURST is possible in 3D, axisymmetric, and plane deformation models (3D, AXIS, D_ PLAN, respectively)

Possible deformation models are PETIT, PETIT_REAC,, GDEF_HYPO_ELAS, GDEF_LOG.

The algorithm used is of the global-local type. Global iterations use the elastic stiffness matrix calculated from the damaged Hooke matrix: \(\underline{\underline{\underline{\underline{\Lambda }}}}\mathrm{=}(1\mathrm{-}D){\underline{\underline{\underline{\underline{\Lambda }}}}}^{0}\)

4.1. Explicit integration

At the level of local iterations (i.e. at each point in GAUSS), the numerical integration of the speed equations is carried out by an explicit Runge-Kutta scheme of order 2 with automatic division into local sub-steps according to an estimate of the integration error (nested Runge-Kutta method) (cf. [R5.03.14]). Test SSNV225A illustrates this method.

4.2. Implicit integration

For implicit integration, the following notations will be used:

\({A}^{\text{-}}\), \(A\) and \(\Delta A\) represent respectively the values of a quantity at the beginning and at the end of the time step in question as well as its increment during the step. The system is discretized using a theta method: \(\Delta A\mathrm{=}\Delta tg({A}^{\text{-}}+\theta \Delta A)\mathrm{=}\Delta tg({A}^{\theta })0<\theta \mathrm{\le }1\)

The discretized problem to be solved is then the following: knowing the state at time \({t}^{\mathrm{-}}\) as well as the deformation increments \(\Delta \varepsilon\) (from the prediction phase, cf. [R5.03.01]) and temperature \(\Delta T\), determine the state of the internal variables at time \(t\) as well as the constraints \(\sigma\).

To properly take into account the variation of elasticity parameters with temperature, it is necessary to discretize (implicitly) the elastic stress-strain relationship in the following way (see for example [R5.03.02]):

\({C}^{\text{-1}}\sigma \mathrm{=}\frac{(1\mathrm{-}D)}{(1\mathrm{-}{D}^{\text{-}})}{({C}^{\text{-}})}^{\text{-1}}{\sigma }^{\text{-}}+(1\mathrm{-}D)(\Delta \varepsilon \mathrm{-}\Delta {\varepsilon }^{\text{th}})\mathrm{-}(1\mathrm{-}D)\Delta {\varepsilon }^{p}\) with \(\Delta {\varepsilon }^{p}\mathrm{=}\Delta p\frac{3}{2}\frac{\tilde{\sigma }}{{\sigma }_{\mathit{eq}}}\mathrm{=}\Delta pn\)

To simplify the expressions and reduce the number of operations during the resolution, it is also possible to write the first equation as a function of the elastic deformation; this supposes storing the elastic or plastic deformations as internal variables. We then obtain:

\(\Delta {\varepsilon }^{e}\mathrm{-}(\Delta \varepsilon \mathrm{-}\Delta {\varepsilon }^{\text{th}})+\Delta p{n}^{\theta }\mathrm{=}0\) with \({n}^{\theta }\mathrm{=}\frac{3}{2}\frac{\tilde{{\sigma }^{\theta }}}{{\sigma }_{\mathit{eq}}^{\theta }}\) (\({F}_{e}\))

and the stresses are then calculated using the plastic deformations at time \({t}^{-}\) by:

\(\sigma \mathrm{=}(1\mathrm{-}D){\sigma }^{\mathit{nd}}\mathrm{=}(1\mathrm{-}D)C{\varepsilon }^{e}\mathrm{=}(1\mathrm{-}D)C({(\varepsilon )}^{\text{-}}\mathrm{-}{({\varepsilon }^{\mathit{th}})}^{\text{-}}\mathrm{-}{({\varepsilon }^{p})}^{\text{-}}+\theta \Delta {\varepsilon }^{e})\)

The following equations are derived from the expressions for the derivatives of the internal variables:

\(\Delta p-\Delta t{\epsilon }_{0}\mathrm{sinh}\left(\frac{{\sigma }_{\mathit{eq}}^{\theta }(1-{H}_{1}^{\theta }-{H}_{2}^{\theta })}{K\left(1-{D}^{\theta }\right)(1-\varphi )}\right)=0\) (\({F}_{p}\))

\(\Delta {H}_{i}\mathrm{-}\frac{{h}_{i}}{{\sigma }_{\mathit{eq}}}({H}_{i}^{\text{*}}\mathrm{-}{\delta }_{i}({H}_{i}^{\text{-}}+\theta \Delta {H}_{i}))\Delta p\mathrm{=}0\), for i=1.2 \(({F}_{{H}_{i}})\)

\(\Delta D\mathrm{-}\Delta t{A}_{0}\mathrm{sinh}(\frac{{\alpha }_{D}\text{<}{\sigma }_{p}^{\theta }{\text{>}}_{\text{+}}+{\sigma }_{\mathit{eq}}^{\theta }(1\mathrm{-}{\alpha }_{D})}{{\sigma }_{0}})\mathrm{=}0\) with \({\sigma }_{p}^{\theta }\mathrm{=}\underset{I}{\mathit{max}}{\sigma }_{I}^{\theta }\text{ou}\mathit{tr}({\sigma }^{\theta })\) (\({F}_{D}\))

We can formally write this system: \(F(\Delta Y)\mathrm{=}0\), with \(\Delta Y\mathrm{=}{(\Delta {\varepsilon }^{\text{e}},\Delta p,\Delta {H}_{\mathrm{1,}}\Delta {H}_{2},\Delta D)}^{t}\)

and \(F(\Delta Y)\mathrm{=}{({F}_{e},{F}_{p},{F}_{{H}_{1}},,{F}_{{H}_{2}},,{F}_{D})}^{\text{t}}\)

This non-linear system is solved by Newton’s iterative method [R5.03.14]):

\(F(\Delta {Y}_{k})+{(\frac{\mathrm{\partial }F}{\mathrm{\partial }\Delta Y})}_{k}(\Delta {Y}_{k+1}\mathrm{-}\Delta {Y}_{k})\) by iterating in \(k\) until convergence.

The Jacobian matrix of the system, necessary for the resolution by Newton’s method, can be calculated either numerically (ALGO_INTE =” NEWTON_PERT “, cf. test SSNV225B), or analytically.

In the latter case the expression of the derivatives is:

\((\frac{\mathrm{\partial }{F}_{e}}{\mathrm{\partial }\Delta {\varepsilon }^{e}})\mathrm{=}{I}_{d}+\Delta p\frac{\mathrm{\partial }{n}^{\theta }}{\mathrm{\partial }\Delta {\varepsilon }^{e}}\) with \(\frac{\mathrm{\partial }{n}^{\theta }}{\mathrm{\partial }\Delta {\varepsilon }^{e}}\mathrm{=}2\mu \theta \frac{(1\mathrm{-}{D}^{\theta })}{{\sigma }_{\mathit{eq}}}\mathrm{[}{I}_{\mathit{dev}}\mathrm{-}n\mathrm{\otimes }n\mathrm{]}\) and \({I}_{\mathit{dev}}\mathrm{=}\frac{3}{2}({I}_{4}\mathrm{-}\frac{1}{3}{I}_{2}\mathrm{\otimes }{I}_{2})\)

\((\frac{\mathrm{\partial }{F}_{e}}{\mathrm{\partial }\Delta p})\mathrm{=}{n}^{\theta }\) \((\frac{\mathrm{\partial }{F}_{e}}{\mathrm{\partial }\Delta D})\mathrm{=}0\)

\((\frac{\mathrm{\partial }{F}_{p}}{\mathrm{\partial }\Delta {\varepsilon }^{\text{e}}})\mathrm{=}\mathrm{-}\Delta t{\varepsilon }_{0}\mathrm{cosh}(\frac{{\sigma }_{\mathit{eq}}^{\theta }(1\mathrm{-}{H}_{1}^{\theta }\mathrm{-}{H}_{2}^{\theta })}{K(1\mathrm{-}{D}^{\theta })(1\mathrm{-}\phi )})\frac{(1\mathrm{-}{H}_{1}^{\theta }\mathrm{-}{H}_{2}^{\theta })}{K(1\mathrm{-}{D}^{\theta })(1\mathrm{-}\phi )}\frac{\mathrm{\partial }{\sigma }_{\mathit{eq}}^{\theta }}{\mathrm{\partial }\Delta {\varepsilon }^{\text{e}}}\)

\(\left(\frac{\partial {F}_{p}}{\partial \Delta p}\right)=1\)

\((\frac{\mathrm{\partial }{F}_{p}}{\mathrm{\partial }\Delta {H}_{i}})\mathrm{=}\Delta t{\varepsilon }_{0}\mathrm{cosh}(\frac{{\sigma }_{\mathit{eq}}^{\theta }(1\mathrm{-}{H}_{1}^{\theta }\mathrm{-}{H}_{2}^{\theta })}{K(1\mathrm{-}{D}^{\theta })(1\mathrm{-}\phi )})\theta \frac{{\sigma }_{\mathit{eq}}^{\mathit{nd}}}{K(1\mathrm{-}\phi )}\)

\((\frac{\mathrm{\partial }{F}_{p}}{\mathrm{\partial }\Delta D})\mathrm{=}0\) because: \(\frac{{\sigma }^{\theta }}{1\mathrm{-}{D}^{\theta }}\mathrm{=}{\sigma }^{\mathit{nd}}\) is independent of \(\Delta D\)

\((\frac{\mathrm{\partial }{F}_{{H}_{i}}}{\mathrm{\partial }\Delta {\varepsilon }^{\text{e}}})\mathrm{=}\frac{{h}_{i}}{{\sigma }_{\mathit{eq}}^{2}}\Delta p({H}_{i}^{\text{*}}\mathrm{-}{\delta }_{i}{H}_{i}^{\theta })\frac{\mathrm{\partial }{\sigma }_{\mathit{eq}}^{\theta }}{\mathrm{\partial }\Delta {\varepsilon }^{\text{e}}}\)

\((\frac{\mathrm{\partial }{F}_{{H}_{i}}}{\mathrm{\partial }\Delta p})\mathrm{=}\mathrm{-}\frac{{h}_{i}}{{\sigma }_{\mathit{eq}}}({H}_{i}^{\text{*}}\mathrm{-}{\delta }_{i}{H}_{i}^{\theta })\)

\((\frac{\mathrm{\partial }{F}_{{H}_{i}}}{\mathrm{\partial }\Delta {H}_{i}})\mathrm{=}1+\frac{{h}_{i}}{{\sigma }_{\mathit{eq}}}{\delta }_{i}\theta \Delta p\)

\((\frac{\mathrm{\partial }{F}_{{H}_{i}}}{\mathrm{\partial }\Delta D})\mathrm{=}\frac{\mathrm{-}{h}_{i}}{{\sigma }_{\mathit{eq}}^{2}}\Delta p({H}_{i}^{\text{*}}\mathrm{-}{\delta }_{i}{H}_{i}^{\theta })\theta {\sigma }_{\mathit{eq}}^{\mathit{nd}}\)

\((\frac{\mathrm{\partial }{F}_{D}}{\mathrm{\partial }\Delta {\varepsilon }^{\text{e}}})\mathrm{=}\mathrm{-}\Delta t\frac{{A}_{0}}{{\sigma }_{0}}\mathrm{cosh}(\frac{{\alpha }_{D}\text{<}{\sigma }_{p}^{\theta }{\text{>}}_{\text{+}}+{\sigma }_{\mathit{eq}}^{\theta }(1\mathrm{-}{\alpha }_{D})}{{\sigma }_{0}})({\alpha }_{D}\frac{\mathrm{\partial }\text{<}{\sigma }_{p}^{\theta }{\text{>}}_{\text{+}}}{\mathrm{\partial }\Delta {\varepsilon }^{\text{e}}}+(1\mathrm{-}{\alpha }_{D})\frac{\mathrm{\partial }{\sigma }_{\mathit{eq}}^{\theta }}{\mathrm{\partial }\Delta {\varepsilon }^{\text{e}}})\)

\((\frac{\mathrm{\partial }{F}_{D}}{\mathrm{\partial }\Delta D})\mathrm{=}1+\frac{\Delta t{A}_{0}\theta }{{\sigma }_{0}}\mathrm{cosh}(\frac{{\alpha }_{D}\text{<}{\sigma }_{p}^{\theta }{\text{>}}_{\text{+}}+{\sigma }_{\mathit{eq}}^{\theta }(1\mathrm{-}{\alpha }_{D})}{{\sigma }_{\mathit{nd}}})\left[{\alpha }_{D}<{\sigma }_{p}^{0}>+(1\mathrm{-}{\alpha }_{D}){\sigma }_{\mathit{eq}}^{\mathit{nd}}\right]\)

with: \(\frac{\mathrm{\partial }{\sigma }_{\mathit{eq}}^{\theta }}{\mathrm{\partial }\Delta {\varepsilon }^{\text{e}}}\mathrm{=}2\mu \theta (1\mathrm{-}{D}^{\theta }){n}^{\theta }\)

and

\(\frac{\mathrm{\partial }\text{<}\mathit{tr}{\sigma }^{\theta }\text{>}}{\mathrm{\partial }\Delta {\varepsilon }^{\text{e}}}\mathrm{=}\frac{\text{<}\mathit{tr}{\sigma }^{\theta }\text{>}}{\mathit{tr}{\sigma }^{\theta }}(3\lambda +2\mu )\theta (1\mathrm{-}{D}^{\theta }){I}_{d}\) in case \({\sigma }_{p}^{\theta }\mathrm{=}\mathit{tr}({\sigma }^{\theta })\)

\(\frac{\mathrm{\partial }\text{<}{\sigma }_{1}^{\theta }\text{>}}{\mathrm{\partial }\Delta {\varepsilon }^{\text{e}}}\mathrm{=}\frac{\text{<}{\sigma }_{1}^{\theta }\text{>}}{{\sigma }_{1}^{\theta }}{I}_{H}\) in the main coordinate system, in case \({\sigma }_{p}^{\theta }\mathrm{=}{\sigma }_{1}\mathrm{=}\underset{I}{\mathit{max}}{\sigma }_{I}^{\theta }\)

with \({I}_{H}\mathrm{=}\left[\begin{array}{cccccc}\lambda +2\mu & \lambda & \lambda & 0& 0& 0\\ \lambda & \lambda +2\mu & \lambda & 0& 0& 0\\ \lambda & \lambda & \lambda +2\mu & 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]\)