Table of Contents

CHAB

Description

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 model

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:

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}\]

Damage model

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:

\[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]\]

\[\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.

Files

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

Subroutines

FileSubroutineDescription
CHAB.F CHAB Material law (Chaboche behaviour model)
CALMAT.FCALMAT Interpolation of material parameters
CALMAT2.FCALMAT2 Interpolation of material parameters
CALARRH Calculation of temperature-dependent parameters using exponential law
ELTENS.F ELTENS Calculation of elasticity tensor
CALMATEL Calculation of the tangent modulus (elastic case)
MATTANG Calculation of the tangent modulus (visco-plastic case)
SOLVNR.F SOLVGI Resolution of the non-linear system to find $\Vert\underline{X}_i\Vert$ in the elastic case
SOLVL Resolution of the non-linear system to find $\Delta p$
SOLVJXI Resolution of the non-linear system to find $\|\underline{X}_i\|$ in the inelastic case
MST_SOLVE.FMST_SOLVE Solution to a real system of linear equations $A*X=B$
CHABDAM.FCHABDAMINC Material law with the computation of isotropic uncoupled damage
post_dam Subroutine for the computation of isotropic machanical damage
post_dam_corr Subroutine for the computation of uniform corrosion damage
CalsigPM Computation of positive and negative part of stress tensor
CalY_dY Computation of Y and dY for damage computation
CHABDAM_C.FCHABDAMSC Material law with the computation of isotropic semi-coupled damage (damage equation is solved at the end of the step)
CHABDAMIC Material law the computation of isotropic coupled damage

Availability

Plane stress state NO
Plane strain state NO
Axisymmetric state NO
3D state YES
Generalized plane state NO

Input file

Parameters defining the type of constitutive law

(2I5, 60A1)
ILLaw number
ITYPE 271
COMMENT Any comment (up to 60 characters) that will be reproduced on the output listing

Integer parameters

(10I5)
NTEMP number of temperatures at which material data is given
IANISOTH = 1 if effect of maximum temperature in the loading history taken into account
= 0 else
MAXITER Maximum number of iterations for Newton-Raphson convergence (default = 25)
NAF Number of Armstrong-Fredericks equations used to define the back-stress X (minimum value=1)
NAFY Number of Armstrong-Fredericks equations taking into account evolution of the mean stress
NAFcyc Number of Armstrong-Fredericks equations taking into account cyclic hardening
NINTV number of time sub-steps in the material law
IDAMUNITS:
= 0 no mechanical damage computation
= 1 isotropic uncoupled damage computation
= 2 isotropic coupled damage computation [UNSTABLE]
= 3 isotropic semi-coupled damage computation (use of $D_{n-1}$ to compute the effective stress at time step $t_n$, $D$ updated at the end of the time step)
TENS:
= 0 no corrosion damage
= 1 linear corrosion damage: $\dot{D}_u=\frac{k_u}{L_E}$
= 2 parabolic corrosion damage: $\dot{D}_u=\frac{k_u}{D_u L_E^2}$
= 3 power law for corrosion damage: $D_u=\frac{k_u}{L_E} t^{m_u}$
IARRH = 1 expression of static recovery parameters using Arrhenius law
= 2 expression of all parameters as exponential function of temperature
= 0 parameters are interpolated linearly between to defined temperatures
ILCF = 1 computation of stress amplitude, mean stress, relaxation stress for cyclic loading (for Optim) - only available for uniaxial loading in x direction.

Real parameters

Parameters not depending on the temperature

Line 1 (4G10)
ETA strain memory rate
PRECNR precision for convergence of the Newton-Raphson algorithm (default=10-4)
PERIOD period of cyclic loading (only if ILCF=1)
$t_H$Hold time in the cyclic loading (only if ILCF=1)
If IANISOTH=1 - Line 2 (3G10)
$b_{D\gamma}$Rate parameter controlling the evolution of Dγ
$b_E$ Rate of evolution of the weighted average factor fe
$f_E^s$ Saturation value of the weighted average factor fe
If IDAM≠0 - Line 3 (7G10)
$h$ micro-defects closure parameter (=0.2 in general for metals ; 1 if micro-defects closure not taken into account)
$D_{crit}$ Critical damage value (<1)
$\tau$ Specific time for the appearance of creep
$k_1$Global safety coefficient on fatigue damage
$k_2$Safety coefficient applied to stress level on fatigue damage
$k_3$Global safety coefficient on creep damage
$k_4$Safety coefficient applied to stress level on creep damage
If IARRH=1 - Line 3+i (i=1:nAF) (2G10)*i
$A_i$ coefficient for expression of $b_i$ using Arrhenius equation: $b_i=A_i \exp(-B_i/T)$
$B_i$ coefficient for expression of $b_i$ using Arrhenius equation: $b_i=A_i \exp(-B_i/T)$

Temperature-dependent parameters - Case where iarrh=0 or iarrh=1

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

Line 1 (4G10)
$T$Temperature
$E$ Young's modulus at temperature T
$\nu$ Poisson's ratio at temperature T
$\alpha$Thermal expansion coefficient at temperature T
NB: in the preprocessor, the thermal expansion coefficient is transformed in its enthalpic formulation
Line 2 (5G10)
$\sigma_y$Yield stress at temperature T
$b$ Rate of isotropic hardening
NB: to avoid convergence issues, b should be constant with temperature
$Q$ Total isotropic saturation size of the yield surface
NB: to avoid convergence issues, Q should be constant with temperature
$K$ Drag stress in Norton-Hoff law
$n$ Viscosity exponent for Norton-Hoff law
Line 2+i (4G10) repeated NAF times (i=1:NAF)
$C_i$Prager's linear coefficient in the ith A-F equation
$\gamma_i$ Dynamic recovery parameter in the ith A-F equation
$b_i$ Static recovery parameter in the ith equation
$r_i$ Static recovery exponent in the ith equation
Line 2+NAF+i (4G10) repeated NAFcyc times (i=1:NAFcyc)
$D_{\gamma,i}$Parameter controlling the evolution of γi with increment of plastic strain norm
$a_{\gamma,i}$Parameter controlling the saturation value of γi
$b_{\gamma,i}$Parameter controlling the saturation value of γi
$c_{\gamma,i}$Parameter controlling the saturation value of γi
Line 2+NAF+NAFcyc+i (4G10) repeated NAFY times (i=1:NAFY)
$\alpha_{b,i}$ Rate of evolution of the mean stress tensor Yi
$Y_{st,i}$ Saturation value of the mean stress tensor Yi
If 1≤IDAM<10 - Line 2+NAF+NAFcyc+NAFY (8G10)
$A$Parameter of correction in the energy stored by hardening
$m$Exponent parameter of correction in the energy stored by hardening
$w_D$Stored energy threshold for damage initiation
$S_f$Fatigue damage parameter
$s_f$Fatigue damage exponent parameter
$S_c$Creep damage parameter
$s_c$Creep damage exponent
$k_c$Kachanov creep damage exponent
If IDAM≥10 - Line 3+NAF+NAFcyc+NAFY (1G10 to 3G10)
$k_u$Uniform corrosion parameter
$L_E$ [Optional] Characteristic length of the element - if blank or 0, $L_E$ is computed as the cubic root of the volume element
$m_u$ [ONLY IF DIDAM=3] power law parameter

Temperature-dependent parameters - Case where iarrh=2

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)\]

Line 1 (6G10)
$A_E$Young modulus parameter
$B_E$Young modulus parameter
$C_E$Young modulus parameter
$A_\nu$Poisson ratio parameter
$B_\nu$Poisson ratio parameter
$C_\nu$Poisson ratio parameter
Line 2 (6G10)
$A_\alpha$Thermal expansion coefficient parameter
$B_\alpha$Dilatation coefficient parameter
$C_\alpha$Dilatation coefficient parameter.
If $C_\alpha=0$, $\int_0^T\alpha(T).dT$ is computed as: $\int_0^T\alpha(T).dT=A_{\alpha}T^2+B_{\alpha}T$
$A_{\sigma_y}$Yield stress parameter
$B_{\sigma_y}$Yield stress parameter
$C_{\sigma_y}$Yield stress parameter
Line 3 (6G10)
$A_K$Drag stress parameter
$B_K$Drag stress parameter
$C_K$Drag stress parameter
$A_n$Norton coefficient parameter
$B_n$Norton coefficient parameter
$C_n$Norton coefficient parameter
Line 4 (2G10)
$b$ Rate of isotropic hardening
$Q$ Total isotropic saturation size of the yield surface
Line 5 (6G10)
BCiParameter for Ci ∀i
CCiParameter for Ci ∀i
BγiParameter for γi ∀i
CγiParameter for γi ∀i
BbiParameter for bi ∀i
CbiParameter for bi ∀i
Line 5+i (4G10) repeated NAF times (i=1:NAF)
ACiPrager's linear coefficient in the ith A-F equation
Aγi Dynamic recovery parameter in the ith A-F equation
Abi Static recovery parameter in the ith equation
ri Static recovery exponent in the ith equation
If nAFcyc>0: Lines 4+nAF+2i and 5+nAF+2i repeated nAFcyc times (i=1:nAFcyc)
Cyclic hardening parameters $a_{\gamma_i}$, $b_{\gamma_i}$, $c_{\gamma_i}$ follow a double exponential law: \[f_{\gamma_i}(T)=A_{f_{\gamma_i}}\left(1-B_{f_{\gamma_i}}\exp\left(\frac{T}{C_{f_{\gamma_i}}}\right)\right)+A_{f_{\gamma_i}}\left(1-C_{f_{\gamma_i}}\exp\left(\frac{T}{E_{f_{\gamma_i}}}\right)\right)\]
Line 4+nAF+2i (3G10)
ADγi Parameter for the rate of cyclic hardening
BDγi Parameter for the rate of cyclic hardening
CDγi Parameter for the rate of cyclic hardening
Line 5+nAF+2i (3G10)
Bfγi Parameter for cyclic hardening ∀f
Cfγi Parameter for cyclic hardening ∀f
Dfγi Parameter for cyclic hardening ∀f
Efγi Parameter for cyclic hardening ∀f
Aaγi Parameter for cyclic hardening parameter $a_{\gamma_i}$
Abγi Parameter for cyclic hardening parameter $b_{\gamma_i}$
Acγi Parameter for cyclic hardening parameter $c_{\gamma_i}$
If nAFY>0: Line 5+nAF+2nAFcyc+i (6G10) repeated nAFY times
Aαb,i Parameter for the rate of evolution of the mean stress
Bαb,i Parameter for the rate of evolution of the mean stress
Cαb,i Parameter for the rate of evolution of the mean stress
AYst,i Parameter for the saturation value of the mean stress
BYst,i Parameter for the saturation value of the mean stress
CYst,i Parameter for the saturation value of the mean stress
If 3≥IDAM≥1
Line 5+NAF+2NAFcyc+nAFY (7G10)
AA stored energy computation parameter
BA stored energy computation parameter
CA stored energy computation parameter
mexponent parameter of correction in the energy stored by hardening
AwD stored energy threshold parameter
BwD stored energy threshold parameter
CwD stored energy threshold parameter
Line 6+NAF+2NAFcyc+nAFY (4G10)
ASf fatigue damage parameter
BSf fatigue damage parameter
CSf fatigue damage parameter
expSfexponent parameter for fatigue damage
Line 7+NAF+2NAFcyc+nAFY (4G10)
The creep damage parameter $S_c$ is calculated using a simpler exponential law: \[A_{S_c}\exp\left(\frac{T}{B_{S_c}}\right)\]
$A_{S_c}$ creep damage parameter
$B_{S_c}$ creep damage parameter
$s_c$exponent parameter for creep damage
$k$Kachanov creep damage exponent
If IDAM≥10
Line 1+NAF+2NAFcyc+nAFY+H(IDAM)*7 (3G10 to 7G10)
$A_{k_u}$ corrosion damage parameter
$B_{k_u}$ corrosion damage parameter
$C_{k_u}$ corrosion damage parameter
$L_E$ [Optional] Characteristic length of the element - if blank or 0, $L_E$ is computed as the cubic root of the volume element
$A_{m_u}$ [ONLY IF DIDAM=3] power law parameter
$B_{m_u}$ [ONLY IF DIDAM=3] power law parameter
$C_{m_u}$ [ONLY IF DIDAM=3] power law parameter

Stresses

Number of stresses

6 for 3D state

Meaning

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

SIG(1)$\sigma_{xx}$
SIG(2)$\sigma_{yy}$
SIG(3)$\sigma_{zz}$
SIG(4)$\sigma_{xy}$
SIG(5)$\sigma_{xz}$
SIG(6)$\sigma_{yz}$

State variables

Number of state variables

$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$.

List of state variables

Q(1)plastic strain norm $p$
Q(2)plastic strain rate norm $\dot{p}$
Q(3)thermal strain $\varepsilon^{th}$
Q(4:9)plastic strain $\underline{\varepsilon}^p$ (6 components)
Q(10:15)total strain $\underline{\varepsilon}$ (6 components)
Q(16) isotropic hardening $R$
Q(17)radius of the plastic strain memory surface $q$
Q(18:23)center of the plastic strain memory surface $\underline{\zeta}$
Q(18+6i:23+6i)Back-stress $\underline{X}_i$ (6 components) for i=1:nAF
Q(18+6nAF+6i:23+6nAF+6i)Modification tensor $\underline{Y}_i$ (6 components) for i=1:nAFY
Q(24+6nAF+6nAFY)Maximum temperature in the loading history

Only if $u_{i_{dam}}>0$

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

Q(25+6NAF+6NAFY) Stored energy $w_s$
Q(26+6NAF+6NAFY) Visco-plastic multiplicator with damage $r$
Q(27+6NAF+6NAFY)
Q(26+ddim+6nAF+6nAFY)
Fatigue damage variable $D_f$ (isotropic) or tensor $\underline{D}_f$ (anisotropic - not implemented)
Q(27+ddim+6NAF+6NAFY)
Q(26+2ddim+6nAF+6nAFY)
Creep damage variable $D_c$ (isotropic) or tensor $\underline{D}_c$ (anisotropic - not implemented)
Q(27+2ddim+6NAF+6NAFY)

Q(32+2ddim+6nAF+6nAFY)
Delayed stress tensor $\sigma^d$

Only if $i_{dam}$≥10

$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)

Only if $i_{LCF}>0$

$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)

Q(NQLCF) t (time)
Q(1+NQLCF) N (cycle)
Q(2+NQLCF) $\sigma_{min}$
Q(3+NQLCF) $\sigma_{max}$
Q(4+NQLCF) $\sigma_{amp}$ (stress amplitude)
Q(5+NQLCF) $\sigma_{mean}$ (mean stress)
Q(6+NQLCF) $\sigma_{H}$ (stress at the end of hold time)
Q(7+NQLCF) $\sigma_{rel}$ (stress relaxation during hold time)