Parallel Rheological Framework

Material Behavior

The parallel rheological framework allows definition of a nonlinear viscoelastic-elastoplastic model consisting of multiple networks connected in parallel, as shown in Figure 1.

The number of viscoelastic networks, N, can be arbitrary; however, at most one equilibrium network (network $0$ in Figure 1) is allowed in the model. The equilibrium network response might be purely elastic or elastic-plastic. In addition, it might include Mullins effect to predict material softening. The definition of the equilibrium network is optional. If it is not defined, the stress in the material will relax completely over time.

The model can be used to predict complex behavior of materials subjected to finite strains, which cannot be modeled accurately using other models available in Abaqus. An example of such complex behavior is depicted in Figure 2, which shows normalized stress relaxation curves for three different strain levels. This behavior can be modeled accurately using the nonlinear viscoelastic model depicted in Figure 3, which can be defined within the framework; but it cannot be captured with the linear viscoelastic model (see Time Domain Viscoelasticity). In the latter case, the three curves would coincide.

Elastic Behavior

The elastic part of the response for all the networks is specified using the hyperelastic material model. Any of the hyperelastic models available in Abaqus can be used (see Hyperelastic Behavior of Rubberlike Materials). The same hyperelastic material definition is used for all the networks, scaled by a stiffness ratio specific to each network. Consequently, only one hyperelastic material definition is required by the model along with the stiffness ratio for each network. The elastic response can be specified by defining either the instantaneous response or the long-term response.

Equilibrium Network Behavior

In addition to the elastic response described above, the response of the equilibrium network can include plasticity and Mullins effect to predict material softening. If the plastic response is defined using isotropic hardening, the response in the equilibrium network is equivalent to that of the permanent set model available in Abaqus (see Permanent Set in Rubberlike Materials for a detailed description of the model). In Abaqus/Standard the nonlinear kinematic hardening model with multiple backstresses can be specified in addition to isotropic plastic hardening. The nonlinear kinematic hardening model is a generalization of the model used for metal plasticity. See Models for Metals Subjected to Cyclic Loading for a detailed description of the model, with the difference that the Cauchy stress is replaced with the Kirchhoff stress in the current formulation.

Viscous Behavior

Viscous behavior must be defined for each viscoelastic network. It is modeled by assuming the multiplicative split of the deformation gradient and the existence of the creep potential, $G^{c r}$ , from which the flow rule is derived. In the multiplicative split the deformation gradient is expressed as

F = F^{e} \cdot F^{c r},

where

F^{e}

is the elastic part of the deformation gradient (representing the hyperelastic behavior) and

F^{c r}

is the creep part of the deformation gradient (representing the stress-free intermediate configuration). The creep potential is assumed to have the general form

G^{c r} = G^{c r} (σ),

where

σ

is the Cauchy stress. If the potential is specified, the flow rule can be obtained from

D^{c r} = \dot{λ} \frac{\partial G^{c r} (σ)}{\partial σ},

where

D^{c r}

is the symmetric part of the velocity gradient,

L^{c r}

, expressed in the current configuration and

\dot{λ}

is the proportionality factor. In this model the creep potential is given by

G^{c r} = \bar{q},

and the proportionality factor is taken as

\dot{λ} = {\dot{\bar{ε}}}^{c r}

, where

\bar{q}

is the equivalent deviatoric Cauchy stress and

{\dot{\bar{ε}}}^{c r}

is the equivalent creep stain rate. In this case the flow rule has the form

D^{c r} = \frac{3}{2 \bar{q}} {\dot{\bar{ε}}}^{c r} \bar{σ},

or, equivalently

D^{c r} = \frac{3}{2 \tilde{q}} {\dot{\bar{ε}}}^{c r} \bar{τ},

where

τ = J σ

is the Kirchhoff stress,

J

is the determinant of

F

\bar{σ}

is the deviatoric Cauchy stress,

\bar{τ}

is the deviatoric Kirchhoff stress, and

\tilde{q} = J \bar{q}

. To complete the derivation, the evolution law for

{\dot{\bar{ε}}}^{c r}

must be provided. In this model

{\dot{\bar{ε}}}^{c r}

can be defined by the power law model, the strain hardening model, the hyperbolic-sine law model, the Bergstrom-Boyce model, the Darveaux model, or a user-defined creep model.

Power Law Model

The power law model is available in the form

{\dot{\bar{ε}}}^{c r} = {\dot{ε}}_{0} {({(\frac{\tilde{q}}{q_{0} + a < p >})}^{n} {[(m + 1) {\bar{ε}}^{c r}]}^{m})}^{\frac{1}{m + 1}},

where

${\dot{\bar{ε}}}^{c r}$: is the equivalent creep strain rate,
${\bar{ε}}^{c r}$: is the equivalent creep strain,
$\tilde{q}$: is the equivalent deviatoric Kirchhoff stress,
$p$: is the Kirchhoff pressure, and
$q_{0}$ , m, n, a, and ${\dot{ε}}_{0}$: are material parameters. It is recommended that you use the power law model rather than the strain hardening model.

Strain Hardening Model

The strain hardening model is available in the form

{\dot{\bar{ε}}}^{c r} = {(A {\tilde{q}}^{n} {[(m + 1) {\bar{ε}}^{c r}]}^{m})}^{\frac{1}{m + 1}},

where

${\dot{\bar{ε}}}^{c r}$: is the equivalent creep strain rate,
${\bar{ε}}^{c r}$: is the equivalent creep strain,
$\tilde{q}$: is the equivalent deviatoric Kirchhoff stress, and
A, m, and n: are material parameters. It is recommended that you use the power law model rather than the strain hardening model. The strain hardening model is a special case of the power law model obtained by setting ${\dot{ε}}_{0} = 1$ , $a = 0$ , and $q_{0} = A^{- 1 / n}$ .

Hyperbolic-Sine Law Model

The hyperbolic-sine law is available in the form

{\dot{\bar{ε}}}^{c r} = A {(\sinh B \tilde{q})}^{n},

where

${\dot{\bar{ε}}}^{c r}$ and $\tilde{q}$: are defined above, and
A, B, and n: are material parameters.

Bergstrom-Boyce Model

Abaqus provides two forms to define the Bergstrom-Boyce creep model. The recommended form of the Bergstrom-Boyce model is defined as

{\dot{\bar{ε}}}^{c r} = {\dot{ε}}_{0} {(λ^{c r} - 1 + E)}^{C} {(\frac{\tilde{q}}{q_{0}})}^{m},

where

λ^{c r} = \sqrt{\frac{1}{3} I : C^{c r}},

and

${\dot{\bar{ε}}}^{c r}$ and $\tilde{q}$: are defined above, and
$q_{0}$ , m, C, E, and ${\dot{ε}}_{0}$: are material parameters.

The original Bergstrom-Boyce model has the form

{\dot{\bar{ε}}}^{c r} = A {(λ^{c r} - 1 + E)}^{C} {(\tilde{q})}^{m},

where

${\dot{\bar{ε}}}^{c r}$ , $λ^{c r}$ , and $\tilde{q}$: are defined above, and
A, m, C, and E: are material parameters.

The recommended form is equivalent to the original form of the Bergstrom-Boyce model. The primary difference between the two formulations is that the recommended form is written in such a way that parameter values do not cause numerical difficulties, which can happen when the original model is calibrated for strain rate applications. In addition, the units of all parameters in the recommended form are physical, which makes unit conversion easier. When the value of the parameter $A$ is very small ( $A \leq 10^{- 27}$ ), the recommended form is obtained by setting $q_{0} = {(\frac{{\dot{ε}}_{0}}{A})}^{\frac{1}{m}}$ and setting ${\dot{ε}}_{0}$ to an arbitrary value greater than zero (typically, ${\dot{ε}}_{0}$ is set to one).

The response of the network defined by the Bergstrom-Boyce model is very similar to the response of the time-dependent network in the hysteresis model (see Hysteresis in Elastomers). However, there are also important differences between the models. In the Bergstrom-Boyce model the equivalent Kirchhoff stress is used instead of the equivalent Cauchy stress, which is used in the hysteresis model. (The two stress measures become equivalent for the case of incompressible materials.) In addition, the material parameters, A, in the hysteresis model and the original form of the Bergstrom-Boyce model differ by a factor of $3 / 2$ . The parameter in the hysteresis model must be multiplied by $3 / 2$ to make the parameters equivalent.

Darveaux Model

The original Darveaux model involves both primary and secondary (steady state) creep. It is intended for applications where loading is applied at the beginning of the analysis, while the modified Darveaux model is recommended for arbitrary loading histories. The secondary creep in the Darveaux model is defined by a hyperbolic-sine law:

{\dot{\bar{ε}}}_{s}^{c r} = C_{s s} {[\sinh (α \tilde{q})]}^{n} .

The secondary law is modified to account for the primary creep effects through

{\dot{\bar{ε}}}^{c r} = {\dot{\bar{ε}}}_{s}^{c r} [1 + ϵ_{T} B \exp (- B {\dot{\bar{ε}}}_{s}^{c r} t)],

where

${\dot{\bar{ε}}}^{c r}$ and $\tilde{q}$: are defined above,
${\dot{\bar{ε}}}_{s}^{c r}$: is the steady-state equivalent creep strain rate,
$C_{s s}$: is the steady-state creep prefactor,
$α$: is the steady-state power law breakdown, and
$n$ , $ϵ_{T}$ , and B: are other material parameters.

Modified Darveaux Model

The modified Darveaux model involves both primary and secondary (steady state) creep, similar to the original Darveaux model. However, unlike the original Darveaux model, which uses total time to describe primary creep effects, the modified Darveaux model uses the accumulated creep strain. This approach provides a solution that does not depend on the loading sequence when the step time has an idle or no-load period. Both Darveaux models provide the same response under constant load and temperature applied at the beginning of the analysis (time zero). However, for arbitrary loading histories, you should use the modified Darveaux model.

The secondary creep is defined by a hyperbolic-sine law:

{\dot{\bar{ε}}}_{s}^{c r} = C_{s s} {[\sinh (α \tilde{q})]}^{n} .

The secondary law is modified to account for the primary creep effects through

{\dot{\bar{ε}}}^{c r} = {\dot{\bar{ε}}}_{s}^{c r} [1 + ϵ_{T} B \exp (- B {\bar{ϵ}}_{s}^{c r})],

where

${\dot{\bar{ε}}}^{c r}$ and $\tilde{q}$: are defined above,
${\bar{ϵ}}_{s}^{c r}$: is the accumulated steady-state equivalent creep strain,
$C_{s s}$ , $α$ , $n$ , $ϵ_{T}$ , and $B$: are defined in the Darveaux model.

User-Defined Model in Abaqus/Standard

A user-defined creep model is available of the following general form:

{\dot{\bar{ε}}}^{c r} = g^{c r} ({\bar{ε}}^{c r}, I_{1}^{c r}, {\bar{I}}_{1}, {\bar{I}}_{2}, J, p, \tilde{q}, t, θ, F V),

where

I_{1}^{c r} = I : C^{c r},

and

${\bar{I}}_{1}$: is the first invariant of $\bar{B}$ ,
${\bar{I}}_{2}$: is the second invariant of $\bar{B}$ ,
${\dot{\bar{ε}}}^{c r}$ , ${\bar{ε}}^{c r}$ , $J$ , $p$ and $\tilde{q}$: are defined above,
t: is the time,
$θ$: is the temperature, and
FV: are field variables.

The tensor $\bar{B}$ is defined in UCREEPNETWORK.

Thermal Expansion

Only isotropic thermal expansion is permitted with nonlinear viscoelastic materials (Thermal Expansion).

Defining Viscoelastic Response

The nonlinear viscoelastic response is defined by specifying the identifier, stiffness ratio, and creep law for each viscoelastic network.

Specifying Network Identifier

Each viscoelastic network in the material model must be assigned a unique network identifier or network id. The network identifiers must be consecutive integers starting with 1. The order in which they are specified is not important.

Input File Usage

Use the following option to specify the network identifier:

VISCOELASTIC, NONLINEAR, NETWORKID=networkId

Defining the Stiffness Ratio

The contribution of each network to the overall response of the material is determined by the value of the stiffness ratio, $r$ , which is used to scale the elastic response of the network material. The sum of the stiffness ratios of the viscoelastic networks must be smaller than or equal to 1. If the sum of the ratios is equal to 1, the purely elastic equilibrium network is not created. If the sum of the ratios is smaller than 1, the equilibrium network is created with a stiffness ratio, $r_{0}$ , equal to

r_{0} = 1 - \overset{N}{\sum_{k = 1}} r_{k},

where $N$ denotes the number of viscoelastic networks and $r_{k}$ is the stiffness ratio of network $k$ . You can specify the stiffness ratio to remain constant during the analysis or to vary as a function of temperature and predefined field variables.

Defining a Constant Stiffness Ratio

You can specify that the stiffness ratio remains constant during the analysis:

Input File Usage

VISCOELASTIC, NONLINEAR, SRATIO=ratio

Defining a Temperature- and Field-Variable Dependent Stiffness Ratio

Alternatively, you can define the stiffness ratio as a function of temperature and predefined field variables.

Input File Usage

NETWORK STIFFNESS RATIO, N=nNetworks, DEPENDENCIES=n

Specifying the Creep Law

The definition of creep behavior in Abaqus/Standard is completed by specifying the creep law.

Power Law Creep Model

The power law model is defined by specifying five material parameters: $q_{0}$ , n, m, a, and ${\dot{ε}}_{0}$ . The parameter ${\dot{ε}}_{0}$ must be positive. It is introduced for dimensional consistency, and its default value is 1.0. For physically reasonable behavior $q_{0}$ and n must be positive, a must be nonnegative (the default is 0.0), and $- 1 < m \leq 0$ . It is recommended that you use the power law model rather than the strain hardening model.

Input File Usage

VISCOELASTIC, NONLINEAR, LAW=POWER LAW

Strain Hardening Creep Model

The strain hardening law is defined by specifying three material parameters: A, n, and m. For physically reasonable behavior A and n must be positive and $- 1 < m \leq 0$ . It is recommended that you use the power law model rather than the strain hardening model.

Input File Usage

VISCOELASTIC, NONLINEAR, LAW=STRAIN

Hyperbolic Sine Creep Model

The hyperbolic sine creep law is specified by providing three nonnegative parameters: A, B, and n.

Input File Usage

VISCOELASTIC, NONLINEAR, LAW=HYPERB

Bergstrom-Boyce Creep Model

The recommended form of the Bergstrom-Boyce creep law is specified by providing five parameters: $q_{0}$ , m, C, E, and ${\dot{ε}}_{0}$ . The parameters ${\dot{ε}}_{0}$ and E must be nonnegative, the parameters $q_{0}$ and m must be positive, and the parameter C must lie in $[- 1, 0]$ .

The original form of the Bergstrom-Boyce creep law is specified by providing four parameters: A, m, C, and E. The parameters A and E must be nonnegative, the parameter m must be positive, and the parameter C must lie in $[- 1, 0]$ .

Input File Usage

Use the following option to define the recommended form of the Bergstrom-Boyce creep law:

VISCOELASTIC, NONLINEAR, LAW=BB POWER LAW

Use the following option to define the original form of the Bergstrom-Boyce creep law:

VISCOELASTIC, NONLINEAR, LAW=BERGSTROM-BOYCE

Darveaux Creep Model

The Darveaux creep law is specified by providing five parameters: $C_{s s}$ , $α$ , $n$ , $ϵ_{T}$ , and $B$ .

Input File Usage

Use the following option to define the Darveaux creep law.

VISCOELASTIC, NONLINEAR, LAW=DARVEAUX

Modified Darveaux Creep Model

The modified Darveaux creep law is specified by providing five parameters: $C_{s s}$ , $α$ , $n$ , $ϵ_{T}$ , and $B$ .

Input File Usage

Use the following option to define the Darveaux creep law.

VISCOELASTIC, NONLINEAR, LAW=MODIFIED DARVEAUX

User-Defined Creep Model

An alternative method for defining the creep law involves using user subroutine UCREEPNETWORK in Abaqus/Standard or VUCREEPNETWORK in Abaqus/Explicit. Optionally, you can specify the number of property values needed as data in the user subroutine.

Input File Usage

VISCOELASTIC, NONLINEAR, LAW=USER, PROPERTIES=n

Numerical Difficulties

Depending on the choice of units, the value of A in the creep models might be very small for typical creep strain rates. If A is less than 10⁻²⁷, numerical difficulties can cause errors in the material calculations; therefore, a different system of units should be used to avoid such difficulties in the calculation of creep strain increments. In such cases it is recommended that you use the creep models that do not have the limitation. You can use the power law model rather than the strain hardening model and the recommended form of the Bergstrom-Boyce model rather than the original form.

Thermorheologically Simple Temperature Effects

Thermorheologically simple temperature effects can be included for each viscoelastic network. In this case the creep law is modified and takes the following form:

\frac{d {\bar{ε}}^{c r}}{d τ} = \frac{1}{a_{T} (θ)} g^{c r} ({\bar{ε}}^{c r}, I_{1}^{c r}, {\bar{I}}_{1}, {\bar{I}}_{2}, J, p, \tilde{q}, τ),

where $τ$ and $a_{T} (θ)$ denote the reduced time and the shift function, respectively. The reduced time is related to the actual time through the integral differential equation

τ = \int_{0}^{t} \frac{d t^{'}}{a_{T} (θ)}, \frac{d τ}{d t} = \frac{1}{a_{T} (θ)} .

Abaqus supports the following forms of the shift function: the Williams-Landel-Ferry (WLF) form, the Arrhenius form, and the tabular form (see Thermorheologically Simple Temperature Effects). In addition, user-defined forms can be specified in Abaqus/Standard.

User-Defined Form in Abaqus/Standard

An alternative method for specifying the shift function involves using user subroutine UTRSNETWORK. Optionally, you can specify the number of property values needed as data in the user subroutine.

Input File Usage

TRS, DEFINITION=USER, PROPERTIES=n

Material Response in Different Analysis Steps

In Abaqus/Standard the material is active during all stress/displacement procedure types. However, the creep effects are taken into account only in quasi-static (Quasi-Static Analysis), coupled temperature-displacement (Fully Coupled Thermal-Stress Analysis), direct-integration implicit dynamic (Implicit Dynamic Analysis Using Direct Integration), and steady-state transport (Steady-State Transport Analysis) analyses. If the material is used in a steady-state transport analysis, it cannot include plasticity. In other stress/displacement procedures the evolution of the state variables is suppressed and the creep strain remains unchanged. In Abaqus/Explicit the creep effects are always active.

Elements

The parallel rheological framework is available with continuum elements that include mechanical behavior (elements that have displacement degrees of freedom), except for one-dimensional elements. The parallel rheological framework is also supported with elements that use the plane stress formulation such as solid plane stress elements, membranes, and shells. However, those elements are not supported with compressible materials. If a compressible material is specified with plane stress elements, Abaqus will modify the material to make it incompressible and issue an informational message.

Output

In addition to the standard output identifiers available in Abaqus (Abaqus/Standard Output Variable Identifiers and Abaqus/Explicit Output Variable Identifiers), the following variables have special meaning for the nonlinear viscoelastic material model:

CEEQ: The overall equivalent creep strain, defined as ${\bar{ε}}^{c r} = \sum_{k = 1}^{N} r_{k} {\bar{ε}}_{k}^{c r}$ .
CE: The overall creep strain, defined as $ε^{c r} = \sum_{k = 1}^{N} r_{k} ε_{k}^{c r}$ .
CENER: The overall viscous dissipated energy per unit volume, defined as $E_{c} = \sum_{k = 1}^{N} E_{c}^{k}$ .
EE: The overall elastic strain, defined as $ε^{e l} = \sum_{k = 1}^{N} r_{k} ε_{k}^{e l}$ .
SENER: The overall elastic strain energy density per unit volume, defined as $E_{s} = \sum_{k = 0}^{N} E_{s}^{k}$ .
SNETk: All stress components in the $k^{t h}$ network ( $0 \leq k \leq 10$ ).

In the above definitions $r_{k}$ denotes the stiffness ratio for network $k$ , $N$ denotes the number of viscoelastic networks, the subscript or superscript $k$ is used to denote network quantities, and the network $0$ is assumed to be the purely elastic network.

If plasticity is specified in the equilibrium network, the standard output identifiers available in Abaqus corresponding to other isotropic and kinematic hardening plasticity models can be obtained for this model as well. In addition, if the Mullins effect is used in the model, the output variables available for the Mullins effect model (see Mullins Effect) can be requested.

References

Bergstrom, J. S., and M. C. Boyce, “Constitutive Modeling of the Large Strain Time-Dependent Behavior of Elastomers,” Journal of the Mechanics and Physics of Solids, vol. 46, pp. 931–954, 1998.
Bergstrom, J. S., and M. C. Boyce, “Large Strain Time-Dependent Behavior of Filled Elastomers,” Mechanics of Materials, vol. 32, pp. 627–644, 2000.
Bergstrom, J. S., and J. E. Bischoff, “An Advanced Thermomechanical Constitutive Model for UHMWPE ,” International Journal of Structural Changes in Solids, vol. 2, pp. 31–39, 2010.
Hurtado, J. A., I. Lapczyk, and S. M. Govindarajan, “Parallel Rheological Framework to Model Non-Linear Viscoelasticity, Permanent Set, and Mullins Effect in Elastomers,” Constitutive Models for Rubber VIII 95, 2013.
Lapczyk, I., J. A. Hurtado, and S. M. Govindarajan, “A Parallel Rheological Framework for Modeling Elastomers and Polymers,” 182^nd Technical Meeting of the Rubber Division of the American Chemical Society, pp. 1840–1859, October 2012, Cincinnati, Ohio.