# UMAT - Unidimensional nonlinear viscoelastic constitutive subroutine

A nonlinear constitutive equation based on the modified superposition method has been implemented in a user material subroutine (UMAT) for Abaqus unidimensional applications. A Python application that generates the Fortran subroutine, based on the number of terms $q_{max}$ and $m_{max}$ employed in the formulation, can be downloaded here.

Requirements: Python

How to use it:

2. Edit the file to change the parameters qmax and mmax according to your material fitting.
3. Execute the command "python abaqus_umat_nlv_beam"

A brief theoretical description on the nonlinear viscoelastic constitutive equation employed in the implementation is presented below.

Theoretical description

For the nonlinear viscoelastic behavior, many investigations have been based on the multiple integral representation due to Green and Rivlin. Several approximate methods have been proposed to simplify the multiple integral representations. The modified superposition method proposed by Findley is implemented here. One can verifiy that this leads to the following relations:

$$$\sigma (t)= \sum\limits_{q = 1}^{q_{max}} {\left( {{E_q}(t)\varepsilon {{(0)}^q} + \int_{\,0}^{\,t} {{E_q}(t - \tau )\frac{{d\left[ {\varepsilon {{(\tau )}^q}} \right]}}{{d\tau }}d\tau } } \right)}\label{eq:eq1}\tag{1}$$$

$$$\sigma (t)=\sum\limits_{q = 1}^{q_{max}} {\left( {{E_q}(0)\varepsilon {{(t)}^q} + \int_{\,0}^{\,t} {\frac{{d{E_q}(t - \tau )}}{{d(t - \tau )}}\varepsilon {{(\tau )}^q}d\tau } } \right)}\label{eq:eq2}\tag{2}$$$

The tensile stress relaxation function may be described by the following Prony series,

$$${E_q}(t) = {E_{q\infty}} + \sum\limits_{m = 1}^{m_{max}} {{E_{qm}}\,\exp \left( { - t/{\tau _m}} \right)}\label{eq:eq3}\tag{3}$$$

with ${m_{max}+1}$ terms.  $E_{q\infty}$, $E_{qm}$ and $\tau_{m}$ are the material parameters. The hereditary integral may be re-written as,

$$$\sigma = \,\sum\limits_{q = 1}^{q_{max}} {\left\{ {{E_{q\infty }}\,\,{\varepsilon ^q} + \,\left( {\sum\limits_{m = 1}^{m_{max}} {{E_{qm}}\,{p_{qm}}} } \right)} \right\}}\label{eq:eq4}\tag{4}$$$

where the state variables ${p_{qm}}$ are described by,

$$${p_{qm}} = \int_{\,0}^{\,t} {\exp \left[ { - (t - \tau )/{\tau _m}} \right]\,\frac{{\,d\left[ {\varepsilon {{(\tau )}^q}} \right]}}{{d\tau }}\,d\tau }\label{eq:eq5}\tag{5}$$$

The state variable increment is given by,

$$$\Delta {p_{qm}}(t) = {p_{qm}}(t + \Delta t) - {p_{qm}}(t)\label{eq:eq6}\tag{6}$$$

being $\Delta t = t_{n+1}-t_{n}$ the time step increment. Replacing Eq. \eqref{eq:eq5} in the first term to the right of Eq. \eqref{eq:eq6} and considering $d\varepsilon (t)/dt$ constant over the interval $[t,t + \Delta t]$ the following approximation is obtained,

\begin{gathered}
{p_{qm}}(t + \Delta t) = \exp \left[ { - \Delta t/{\tau _m}} \right]{p_{qm}}(t) + \int_{{\kern 1pt} t}^{{\kern 1pt} t + \Delta t} {\exp \left[ { - (t + \Delta t - \tau )/{\tau _m}} \right]{\mkern 1mu} \left( {\frac{{{\mkern 1mu} d[\varepsilon {{(\tau )}^q}]}}{{d\tau }}} \right){\mkern 1mu} d\tau } \\
\approx \exp \left[ { - \Delta t/{\tau _m}} \right]{p_{qm}}(t) + \left( {\frac{{{\mkern 1mu} \Delta {\varepsilon ^q}}}{{\Delta t}}} \right){\mkern 1mu} {\tau _m}{\mkern 1mu} \left( {1 - \exp \left[ { - \Delta t/{\tau _m}} \right]} \right) \\
\end{gathered}
\label{eq:eq7}\tag{7}\

Different approximation algorithms are described by Sorvari and Hamalainen. Replacing Eq. \eqref{eq:eq7} into \eqref{eq:eq6} leads to,

\Delta {p_{qm}} \approx {p_{qm}}{\mkern 1mu} {A_m} + \left( {\frac{{{\mkern 1mu} \Delta {\varepsilon ^q}}}{{\Delta t}}} \right){\mkern 1mu} {B_m}
\label{eq:eq8}\tag{8}

where $\Delta\varepsilon = \varepsilon(t_{n+1})-\varepsilon(t_{n})$ is the strain increment and ${A_m}$ and  ${B_m}$ are described by,

$$${A_m} = \exp \left[ { - \Delta t/{\tau _m}} \right] - 1\label{eq:eq9}\tag{9}$$$
$$${B_m} = {\tau _m}\,\left( {1 - \exp \left[ { - \Delta t/{\tau _m}} \right]} \right)\label{eq:eq10}\tag{10}$$$

The derivative of the state variable increment in relation to the strain increment is given by,

\frac{{\partial \Delta {p_{qm}}}}{{\partial \Delta \varepsilon }} = q{\mkern 1mu} {\left( {\varepsilon + \Delta \varepsilon } \right)^{q - 1}}\frac{{{B_m}}}{{\Delta t}}
\label{eq:eq11}\tag{11}

The stress increment can be obtained using Eq. \eqref{eq:eq4},

\Delta \sigma = {\mkern 1mu} \sum\limits_{q = 1}^{{q_{max}}} {\left\{ {{E_{q\infty }}{\mkern 1mu} \Delta {\varepsilon ^q} + {\mkern 1mu} \left( {\sum\limits_{m = 1}^{{m_{max}}} {{E_{qm}}\Delta {p_{qm}}} } \right)} \right\}}
\label{eq:eq12}\tag{12}

The material Jacobian matrix is then obtained deriving Eq. \eqref{eq:eq12} with respect to the strain increment. By replacing Eq. \eqref{eq:eq11} into the resultant equation, the following relation is obtained,

\frac{{\partial \Delta \sigma }}{{\partial \Delta \varepsilon }} = {\mkern 1mu} \sum\limits_{q = 1}^{{q_{max}}} {\left\{ {{E_{q\infty }}{\mkern 1mu} q{\mkern 1mu} {{\left( {\varepsilon + \Delta \varepsilon } \right)}^{q - 1}} + {\mkern 1mu} \left( {\sum\limits_{m = 1}^{{m_{max}}} {{E_{qm}}{\mkern 1mu} q{\mkern 1mu} {\mkern 1mu} {{\left( {\varepsilon + \Delta \varepsilon } \right)}^{q - 1}}\frac{{{B_m}}}{{\Delta t}}{\mkern 1mu} {\mkern 1mu} } } \right)} \right\}}
\label{eq:eq13}\tag{13}

which is used for numerical implementation into the finite element package user subroutine.

******************************************************************************
Revision Notes:
1.0 - Initial version
2.0 - Nonlinear strain increment corrected
******************************************************************************

Author      : Marcelo Caire