Chaboche elasto-visco-plastic constitutive model with thermal and cyclic effects for solid elements at constant or variable temperatures with damage computation.
The model is highly adjustable and can be used for very simple laws (elastic, bilinear plasticity) as well as for more complex behaviors (visco-plasticity, isotropic hardening, kinematic hardening, cyclic hardening, …)

Implemented by: Hélène Morch, 2016-2020

The Chaboche model is a viscoplastic constitutive model that allows thermo-mechanical cyclic analysis of solids.

The main equations of the model are summarized hereafter:

The total strain is decomposed into thermal $\underline{\varepsilon}^{th}$, elastic $\underline{\varepsilon}^{e}$, and plastic $\underline{\varepsilon}^{p}$ strains: \[\underline{\varepsilon}=\underline{\varepsilon}^{th}+\underline{\varepsilon}^{e}+\underline{\varepsilon}^{p}\] The stress and strain are related through Hooke's law: \[\underline{\sigma}=\underline{\underline{E}}\underline{\varepsilon}\] The yield locus is defined using the von Mises criterion: \[\sigma_v=\|\underline{\sigma}-\underline{X}\|-R-\sigma_y\leq 0\]

$\sigma_v>0$ corresponds to the visco-plastic domain. The viscosity function of the model is the Norton-Hoff equation: \[\dot{p}=\left<\frac{\sigma_v}{K}\right>^n \] Where:

• $\dot{p}$ is the plastic strain rate norm
• $\left<x\right>= \begin{cases} x & \text{if } x>0 \\ 0 & \text{if } x\leq 0 \end{cases}$

Parameter $K$ can be put to a value close to 0 to model elasto-plasticity (without viscous effects).
$R$ is the isotropic hardening variable, which evolves with the plastic strain rate: \[ \dot{R}=b(Q-R)\dot{p}\]

The back-stress $\underline{X}$ controls kinematic hardening. In the Chaboche model, the back-stress is determined through the summation of Armstrong-Fredericks equations: \[\underline{X}=\displaystyle\sum_{i=1}^{nAF} \underline{X}_i\] \[\underline{\dot{X}}_i=\frac{2}{3}C_i\underline{\dot{\varepsilon}}^p-\gamma_i(\underline{X}_i-\underline{Y}_i)\dot{p}-b_i\|\underline{X}_i\|^{r_i-1}\underline{X}_i+\frac{1}{C_i}\frac{dC_i}{dT}\dot{T}\underline{X}_i\] The two first terms of the equation correspond to the classic Armstrong-Fredericks equation. The third term $-b_i\|\underline{X}_i\|^{r_i-1}\underline{X}_i$ is added to model static recovery. This term can be used to model creep and relaxation of the material along with the Norton-Hoff equation.
The last term takes into account the effect of temperature variation.

The parameter $\gamma_i$ can be made to evolve in order to take into account cyclic hardening of the material: \[\dot{\gamma_i}=D_{\gamma_i}(\gamma_i^0-\gamma_i)\dot{p}\] \[\gamma_i^0=a_{\gamma_i}+b_{\gamma_i}\exp(-c_{\gamma_i}q)\]

Where $q$ is the radius of the plastic strain memory surface $g_M$. $q$ can be understood as the maximum value of the norm of the plastic strain $p$ in the loading history. The equations controlling it are written below, where $\eta$ is a material parameter equal to 0.5 and $H()$ is the Heaviside step function: \[g_M(\underline{\varepsilon^p}-\underline{\zeta})=\|\underline{\varepsilon^p}-\underline{\zeta}\|-q \\ \dot{q}=\eta H(g_M)*\left<\frac{\underline{\sigma}-\underline{X}}{\|\underline{\sigma}-\underline{X}\|}:\frac{\underline{\varepsilon^p}-\underline{\zeta}}{\|\underline{\varepsilon^p}-\underline{\zeta}\|}\right>\dot{p} \\ \underline{\dot{\zeta}}=(1-\eta)H(g_M)*\left<\frac{\underline{\sigma}-\underline{X}}{\|\underline{\sigma}-\underline{X}\|}:\frac{\underline{\varepsilon^p}-\underline{\zeta}}{\|\underline{\varepsilon^p}-\underline{\zeta}\|}\right>\sqrt{\frac{2}{3}}\frac{\underline{\varepsilon^p}-\underline{\zeta}}{\|\underline{\varepsilon^p}-\underline{\zeta}\|}\dot{p} \]

In the Armstrong-Fredericks equations, a modification tensor $\underline{Y}_i$ can be used to model mean stress evolution. This tensor evolves with respect to the norm of the back-stress: \[\dot{\underline{Y}}_i=-\alpha_{b,i}\left(\frac{3}{2}Y_{st,i}\frac{\underline{X}_i}{\|\underline{X}_i\|}+\underline{Y}_i\right)\|\underline{X}_i\|^{r_i}\]

In the case of a cyclic anisothermal loading, some parameters are influenced by the maximum temperature of the cycle: \[D_{\gamma_i}=b_{D_\gamma}(D_{\gamma_i}^{T_{max}}-D_{\gamma_i})\dot{p}\] \[E=f_EE+(1-f_E)E_{T_{max}}\] \[\dot{f}_E=b_E(f_E^S-f_E)\dot{p}\]

The calculation of damage can be included when parameter idam > 0.
The damage $D$ is split into 2 types: fatigue damage $D_f$ and creep damage $D_c$, which evolve according to the following equations: \[D=D_f+D_c\] \[\dot{D}_f=k_1\left(\frac{Y(\sigma*k_2)}{S_f}\right)^{s_f}\frac{dp}{dt}\] \[\dot{D}_c=k_3\left(\frac{Y(\sigma^d*k_4)}{S_c}\right)^{s_c}\frac{1}{(1-D)^k}\] Where:

• $k_1, k_2, k_3, k_4$ are safety coefficients
• $S_f, s_f, S_c, s_c$ are temperature dependent material parameters
• $Y(\sigma)$ is a variable that takes into account the difference between tensile and compressive state for the evolution of damage: under compressive stress, less damage appears in the material, therefore the compressive part of the stress tensor is multiplied by a factor $h$ < 1:

\[Y=\frac{1+\nu}{2E}\left[\frac{\langle\sigma\rangle^+_{ij}\langle\sigma\rangle^+_{ij} }{(1-D)^2}+h \frac{\langle\sigma\rangle^-_{ij}\langle\sigma\rangle^-_{ij} }{(1-hD)^2}\right]-\frac{\nu}{2E}\left[\frac{\langle\sigma_{kk}\rangle^2}{(1-D)^2}+h\frac{\langle-\sigma_{kk}\rangle^2}{(1-hD)^2}\right]\]

• $\sigma^d$ is the delayed stress level used to take into account the fact that creep damage only develops after a certain amount of time $\tau$:

\[\frac{d\sigma^d}{dt}=\frac{\sigma-\sigma^d}{\tau}\]

The safety coefficients are used to numerically accelerate the evolution of damage ($k_i \ge 1$) in order to take into account the uncertainty on material data. The figure below shows the evolution of the damage variables for a reference simulation (in blue) and for the same simulation with each coefficient set to a value of 2.0.
Coefficients $k_2$ and $k_4$ which are applied to the stress level have a much stronger influence on the damage evolution, therefore, values of these coefficients should not be too high.

Prepro: LCHAB.F
Lagamine: CHAB.F, CHABDAM.F, CHABDAM_C.F

2+IANISOTH+NAF+NAFcyc+NAFY lines repeated NTEMP times
Parameters must be introduced by increasing temperature order

Parameters in this case follow an exponential law of the type: \[P(T)=A_P\left(1-B_P*\exp\left(\frac{T}{C_P}\right)\right)\]

6 for 3D state

The stresses are the components of CAUCHY stress tensor in global (X,Y,Z) coordinates For the 3-D state:

$24+6n_{AF}+6n_{AF_Y}+(8+2ddim)\mathscr{H}(u_{i_{dam}})+2\mathscr{H}(d_{i_{dam}})+8i_{LCF}$
Where: $u_{i_{dam}}\equiv i_{dam} \mod 10$
and $d_{i_{dam}}=i_{dam}-u_{i_{dam}}$.
$\mathscr{H}(x)$ is the Heaviside step function: $\mathscr{H}(x)=1$ if and only if $x>0$, otherwise, $\mathscr{H}(x)=0$.

In the following table, ddim=1 for isotropic damage (scalar damage variable $D$) and ddim=6 for anisotropic damage (not implemented).

$N_{Q,D_u}=25+6n_{AF}+6n_{AF_Y}+(8+2ddim)*\mathscr{H}(u_{i_{dam}})$
with $\mathscr{H}(u_{i_{dam}})=1$ if and only if $u_{i_{dam}}>0$

Q(NQDU)	$D_u$ - Uniform corrosion damage
Q(NQDU+1)	$L_E=\sqrt[3]{V_E}$ - Characteristic length of the element where $V_E$ is the volume of the element ( only works with BWD3T elements)

$N_{Q,LCF}=25+6n_{AF}+6n_{AF_Y}+(8+2ddim)*\mathscr{H}(u_{i_{dam}})+2d_{i_{dam}}$
where $d_{i_{dam}}=1$ if $i_{dam}$≥10 and 0 otherwise)