CREEP

Metals

For metals whose material behavior includes metal creep and/or time-dependent volumetric swelling, the routine allows any “creep” and “swelling” laws (viscoplastic behavior) of the following general form to be defined:

\dot{{\bar{ε}}^{c r}} = g^{c r} (p, \tilde{q}, {\bar{ε}}^{s w}, {\bar{ε}}^{c r}, time, \dots),

\dot{{\bar{ε}}^{s w}} = g^{s w} (p, \tilde{q}, {\bar{ε}}^{s w}, {\bar{ε}}^{c r}, time, \dots),

where

${\bar{ε}}^{c r}$: is the uniaxial equivalent “creep” strain, conjugate to $\tilde{q}$ , the Mises or Hill equivalent stress;
${\bar{ε}}^{s w}$: is the volumetric swelling strain;
p: is the equivalent pressure stress, $p = - \frac{1}{3} (σ_{11} + σ_{22} + σ_{33})$ ; and
$\tilde{q}$: is the equivalent deviatoric stress (Mises' or, if anisotropic creep behavior is defined, Hill's definition).

The user subroutine must define the increments of inelastic strain, $△ {\bar{ε}}^{c r}$ and $△ {\bar{ε}}^{s w}$ , as functions of p and $\tilde{q}$ and any other variables used in the definitions of $g^{c r}$ and $g^{s w}$ (such as solution-dependent state variables introduced by you) and of the time increment, $△ t$ . If any solution-dependent state variables are included in the definitions of $g^{c r}$ and $g^{s w}$ , they must also be integrated forward in time in this routine.

Abaqus computes the incremental creep strain (or the incremental viscoplastic strain) components as

Δ ε^{c r} = \frac{1}{3} Δ {\bar{ε}}^{s w} R + Δ {\bar{ε}}^{c r} n,

where $n$ is the gradient of the deviatoric stress potential, defined as

n = \frac{\partial \tilde{q}}{\partial σ},

and $R$ is a matrix with the anisotropic swelling ratios in the diagonal if anisotropic swelling is defined; otherwise, $R = I$ .

Drucker-Prager Materials

For materials that yield according to the extended Drucker-Prager plasticity models using Drucker-Prager creep, the routine allows any “creep” laws (viscoplastic behavior) of the following general form to be defined:

\dot{{\bar{ε}}^{c r}} = g^{c r} ({\bar{σ}}^{c r}, {\bar{ε}}^{c r}, time, \dots),

where

${\bar{σ}}^{c r}$

is the equivalent creep stress defined as

$\begin{matrix} \frac{q - p \tan β}{1 - \frac{1}{3} \tan β } \end{matrix}$: if creep is defined in terms of uniaxial compression,
$\begin{matrix} \frac{q - p \tan β}{1 + \frac{1}{3} \tan β} \end{matrix}$: if creep is defined in terms of uniaxial tension, and
$\begin{matrix} q - p \tan β \end{matrix}$: if creep is defined in terms of pure shear,

where q is the equivalent deviatoric Mises' stress, p is the pressure stress, and $β$ is the friction angle, and

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

is the uniaxial equivalent “creep” strain, conjugate to ${\bar{σ}}^{c r}$ such that ${\bar{σ}}^{c r} Δ {\bar{ε}}^{c r} = σ_{i j} Δ ε_{i j}^{c r}$ .

The user subroutine must define the increment of inelastic strain, $△ {\bar{ε}}^{c r}$ , as a function of ${\bar{σ}}^{c r}$ and any other variables used in the definitions of $g^{c r}$ (such as solution-dependent state variables introduced by you) and of the time increment, $△ t$ . If any solution-dependent state variables are included in the definitions of $g^{c r}$ , they must also be integrated forward in time in this routine.

Abaqus computes the incremental creep strain (or the incremental viscoplastic strain) components as

Δ ε^{c r} = \frac{Δ {\bar{ε}}^{c r}}{f^{c r}} (\frac{q}{\sqrt{{({ϵ \bar{σ} |}_{0} \tan ψ)}^{2} + q^{2}}} n + \frac{1}{3} \tan ψ I),

where $n = \partial \tilde{q} / \partial σ .$ The variable $f^{c r}$ is determined in such a way that

f^{c r} = \frac{1}{{\bar{σ}}^{c r}} σ : \frac{\partial G^{c r}}{\partial σ},

and

G^{c r} = \sqrt{{({ϵ \bar{σ} |}_{0} \tan ψ)}^{2} + q^{2}} - p \tan ψ

is the hyperbolic creep potential, where $ψ (θ, f^{α})$ is the dilation angle measured in the p–q plane at high confining pressure, ${\bar{σ} |}_{0} = {\bar{σ} |}_{{\bar{ε}}^{p l} = 0, {\dot{\bar{ε}}}^{p l} = 0}$ is the initial yield stress, and $ϵ$ is the eccentricity. See Extended Drucker-Prager Models for a discussion of $ψ$ , $ϵ$ , and ${\bar{σ} |}_{0}$ .

Capped Drucker-Prager Materials

For materials that yield according to the modified Drucker-Prager/Cap plasticity model using cap creep, the routine allows any “cohesion creep” and “consolidation creep” laws (viscoplastic behavior) of the following general form to be defined:

\dot{{\bar{ε}}_{s}^{c r}} = g_{s}^{c r} ({\bar{σ}}^{c r}, {\bar{ε}}_{s}^{c r}, time, \dots),

\dot{{\bar{ε}}_{c}^{c r}} = g_{c}^{c r} ({\bar{p}}^{c r}, {\bar{ε}}_{c}^{c r}, time, \dots),

where

${\bar{σ}}^{c r}$

is the equivalent creep stress defined from uniaxial compression test data as

{\bar{σ}}^{c r} = \frac{q - p \tan β}{(1 - \frac{1}{3} \tan β)},

where q is the equivalent deviatoric Mises' stress, p is the pressure stress, and $β$ is the friction angle;

${\bar{ε}}_{s}^{c r}$

is the equivalent cohesion creep uniaxial strain, conjugate to ${\bar{σ}}^{c r}$ such that ${\bar{σ}}^{c r} Δ {\bar{ε}}_{s}^{c r} = σ : Δ ε_{s}^{c r}$ , where $Δ ε_{s}^{c r}$ is defined below;

${\bar{p}}^{c r} = p - p_{a}$

is the effective creep pressure ( $p = - \frac{1}{3} (σ_{11} + σ_{22} + σ_{33})$ and $p_{a}$ is the cap hardening parameter); and

${\bar{ε}}_{c}^{c r}$

is the volumetric consolidation creep strain.

The user subroutine must define the increments of inelastic strain, $△ {\bar{ε}}_{s}^{c r}$ and/or $△ {\bar{ε}}_{c}^{c r}$ , as functions of ${\bar{σ}}^{c r}$ and/or ${\bar{p}}^{c r}$ and any other variables used in the definitions of $g_{s}^{c r}$ and $g_{c}^{c r}$ (such as solution-dependent state variables introduced by you) and of the time increment, $△ t$ . If any solution-dependent state variables are included in the definitions of $g_{s}^{c r}$ and $g_{c}^{c r}$ , they must also be integrated forward in time in this routine.

Calculation of Incremental Creep Strains for the Cohesion Mechanism

Abaqus computes the incremental creep strain (or the incremental viscoplastic strain) components of the cohesion mechanism as

Δ ε_{s}^{c r} = \frac{Δ {\bar{ε}}_{s}^{c r}}{f^{c r}} (\frac{q}{\sqrt{{(0.1 \frac{d}{(1 - \frac{1}{3} \tan β)} \tan β)}^{2} + q^{2}}} n + \frac{1}{3} \tan β I),

where $n = \partial \tilde{q} / \partial σ,$ d is the material cohesion, the variable $f^{c r}$ is determined in such a way that

f^{c r} = \frac{1}{{\bar{σ}}^{c r}} σ : \frac{\partial G_{s}^{c r}}{\partial σ},

and $G_{s}^{c r}$ is the cohesion creep potential

G_{s}^{c r} = \sqrt{{(0.1 \frac{d}{(1 - \frac{1}{3} \tan β)} \tan β)}^{2} + q^{2}} - p \tan β .

Calculation of Incremental Creep Strains for the Consolidation Mechanism

Abaqus computes the incremental creep strain (or the incremental viscoplastic strain) components of the consolidation mechanism as

Δ ε_{c}^{c r} = \frac{△ {\bar{ε}}_{c}^{c r}}{G_{c}^{c r}} (R^{2} q n - \frac{1}{3} (p - p_{a}) I),

where R controls the shape of the cap, and $G_{c}^{c r}$ is the consolidation creep potential

G_{c}^{c r} = \sqrt{{(p - p_{a})}^{2} + {(R q)}^{2}} .

Cohesion material properties are determined with a uniaxial compression test in which $d {\bar{ε}}_{s}^{c r} = ∥ d ϵ_{11}^{c r} ∥$ , and consolidation material properties are determined with a volumetric compression test in which $d {\bar{ε}}_{c}^{c r} = ∥ d ϵ_{v o l}^{c r} ∥$ . Most likely, $g_{s}^{c r}$ is a positive function of ${\bar{σ}}^{c r}$ , and $g_{c}^{c r}$ is a positive function of ${\bar{p}}^{c r}$ .

Gaskets

For gaskets whose behavior includes creep, the routine allows any “creep” law of the following general form to be defined:

{\dot{ε}}^{c r} = g^{c r} (σ, ε^{c r}, time, \dots),

where $ε^{c r}$ is the compressive creep strain, conjugate to $σ$ , the compressive stress in the gasket.

The user subroutine must define the increments of inelastic creep strain, $△ ε^{c r}$ , as functions of $σ$ and any other variables used in the definitions of $g^{c r}$ (such as solution-dependent state variables introduced by you) and of the time increment, $△ t$ . If any solution-dependent state variables are included in the definitions of $g^{c r}$ , they must also be integrated forward in time in this routine. Abaqus will automatically multiply this creep strain by the proper thickness (see Defining the Gasket Behavior Directly Using a Gasket Behavior Model) to obtain a creep closure.

Integration Schemes

Abaqus provides both explicit and implicit time integration of creep and swelling behavior defined in this routine. The choice of the time integration scheme depends on the procedure type, the procedure definition, and whether a geometric linear or nonlinear analysis is requested (see Rate-Dependent Plasticity: Creep and Swelling).

Implicit integration is generally more effective when the response period is long relative to typical relaxation times for the material. Simple high-temperature structural design applications usually do not need implicit integration, but more complicated problems (such as might arise in manufacturing processes), creep buckling applications, or nonstructural problems (such as geotechnical applications) often are integrated more efficiently by the implicit method provided in the program. If implicit integration is used with this subroutine, nonlinear equations must be solved at each time step and the variations of $△ {\bar{ε}}^{c r}$ , $△ {\bar{ε}}^{s w}$ , $△ {\bar{ε}}_{s}^{c r}$ , or $△ {\bar{ε}}_{c}^{c r}$ with respect to ${\bar{ε}}^{c r}$ , ${\bar{ε}}^{s w}$ , ${\bar{ε}}_{s}^{c r}$ , ${\bar{ε}}_{c}^{c r}$ , p, $\tilde{q}$ , ${\bar{p}}^{c r}$ , or ${\bar{σ}}^{c r}$ must be defined in the subroutine. To obtain good convergence during implicit integration, it is essential to define these quantities accurately.

At the start of a new increment the subroutine is called once for each integration point to calculate the estimated creep strain based on the state at the start of the increment. Subsequently, it is called twice for each iteration if explicit integration is used: once to calculate the creep strain increment at the start of the increment and once to calculate it at the end of the increment. This is needed to test the validity of the time increment with respect to the user-specified maximum allowable difference in the creep strain increment. The flag LEND indicates whether the routine is called at the start or the end of the increment. The subroutine must use the corresponding values of time, temperature, field variables, and solution-dependent state variables in the calculation of the creep strain increment.

For implicit integration Abaqus uses a local iteration procedure to solve the nonlinear constitutive equations, and the subroutine is called multiple times. The exact number of calls depends on the convergence rate of the local iteration procedure and, hence, will vary from point to point. During these iterations it is possible for the values of the state variables to be far from their final values when the equations are solved. Therefore, the coding in the subroutine must adequately protect against arithmetic failures (such as floating point overflows) even when variables are passed in with physically unreasonable values. As in explicit integration, the variable LEND indicates whether the routine is called at the start or the end of the increment.

Constant Stress Assumption When Defining Creep and Swelling

When the creep and swelling behavior are defined by simple formulæ, it is often possible to calculate the increments of equivalent creep and swelling strain exactly if it is assumed that the stress is constant during the increment. This approach has the advantage that it provides very good accuracy within the constant stress assumption. It also avoids the problem that arises for some creep behavior definitions: that the creep strain rate becomes infinite at zero time (or strain). Otherwise, in such a case you must protect against causing arithmetic failures at the start of the solution.

Defining Both Plasticity and Creep

If both plasticity and creep are defined for a material, Abaqus will calculate the creep strain before entering the plasticity routines. The stresses passed into the creep routine may, therefore, exceed the yield stress.

Interpretation of Stress and Strain Variables

In finite-strain applications strain variables should be interpreted as logarithmic strains and stresses as “true” stress.

User Subroutine Interface

      SUBROUTINE CREEP(DECRA,DESWA,STATEV,SERD,EC,ESW,P,QTILD,
     1 TEMP,DTEMP,PREDEF,DPRED,TIME,DTIME,CMNAME,LEXIMP,LEND,
     2 COORDS,NSTATV,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
C
      INCLUDE 'ABA_PARAM.INC'
C
      CHARACTER*80 CMNAME
C
      DIMENSION DECRA(5),DESWA(5),STATEV(*),PREDEF(*),DPRED(*),
     1 TIME(3),EC(2),ESW(2),COORDS(*)

      user coding to define DECRA,  DESWA

      RETURN
      END

Variables to Be Defined

In all cases

DECRA(1)

The definition depends on the usage:

Metal creep: $△ {\bar{ε}}^{c r}$ , equivalent (uniaxial) deviatoric creep strain increment.
Drucker-Prager creep: $△ {\bar{ε}}^{c r}$ , equivalent (uniaxial) creep strain increment.
Capped Drucker-Prager creep: $△ {\bar{ε}}_{s}^{c r}$ , equivalent (uniaxial) cohesion creep strain increment.
Gasket creep: $△ ε^{c r}$ , uniaxial compressive creep strain increment.

DESWA(1)

The definition depends on the usage:

Metal creep: $△ {\bar{ε}}^{s w}$ , volumetric swelling strain increment.
Capped Drucker-Prager creep: $△ {\bar{ε}}_{c}^{c r}$ , equivalent (volumetric) consolidation creep strain increment.
Drucker-Prager and gasket creep: = 0.

For implicit creep integration (LEXIMP=1, see below)

DECRA(2)

The definition depends on the usage:

Metal creep and Drucker-Prager creep: $\partial △ {\bar{ε}}^{c r} / \partial {\bar{ε}}^{c r}$ .
Capped Drucker-Prager creep: $\partial △ {\bar{ε}}_{s}^{c r} / \partial {\bar{ε}}_{s}^{c r}$ .
Gasket creep: $\partial △ ε^{c r} / \partial ε^{c r}$ .

DECRA(3)

The definition depends on the usage:

Metal creep: $\partial △ {\bar{ε}}^{c r} / \partial {\bar{ε}}^{s w}$ .
Drucker-Prager creep, gasket creep, and capped Drucker-Prager creep: = 0.

DECRA(4)

The definition depends on the usage:

Metal creep: $\partial △ {\bar{ε}}^{c r} / \partial p$ .
Drucker-Prager creep, gasket creep, and capped Drucker-Prager creep: = 0.

DECRA(5)

The definition depends on the usage:

Metal creep: $\partial △ {\bar{ε}}^{c r} / \partial \tilde{q}$ .
Drucker-Prager creep: $\partial △ {\bar{ε}}^{c r} / \partial {\bar{σ}}^{c r}$ .
Capped Drucker-Prager creep: $\partial △ {\bar{ε}}_{s}^{c r} / \partial {\bar{σ}}^{c r}$ .
Gasket creep: $\partial △ ε^{c r} / \partial σ$ .

DESWA(2)

The definition depends on the usage:

Metal creep: $\partial △ {\bar{ε}}^{s w} / \partial {\bar{ε}}^{c r}$ .
Drucker-Prager creep, gasket creep, and capped Drucker-Prager creep: = 0.

DESWA(3)

The definition depends on the usage:

Metal creep: $\partial △ {\bar{ε}}^{s w} / \partial {\bar{ε}}^{s w}$ .
Capped Drucker-Prager creep: $\partial △ {\bar{ε}}_{c}^{c r} / \partial {\bar{ε}}_{c}^{c r}$ .
Drucker-Prager and gasket creep: = 0.

DESWA(4)

The definition depends on the usage:

Metal creep: $\partial △ {\bar{ε}}^{s w} / \partial p$ .
Capped Drucker-Prager creep: $\partial △ {\bar{ε}}_{c}^{c r} / \partial {\bar{p}}^{c r}$ .
Drucker-Prager and gasket creep: = 0.

DESWA(5)

The definition depends on the usage:

Metal creep: $\partial △ {\bar{ε}}^{s w} / \partial \tilde{q}$ .
Drucker-Prager creep, gasket creep, and capped Drucker-Prager creep: = 0.

Variables That Can Be Updated

STATEV

An array containing the user-defined solution-dependent state variables at this point. This array will be passed in containing the values of these variables at the start of the increment unless they are updated in user subroutine USDFLD or UEXPAN, in which case the updated values are passed in. If any of the solution-dependent variables are being used in conjunction with the creep behavior and the routine was called at the end of the increment (LEND=1, see the definition of LEND below), they must be updated in this subroutine to their values at the end of the increment. Furthermore, if the solution-dependent state variables are defined as a function of the creep (swelling) strain increment, they must be updated based on the creep (swelling) strain increment computed as EC(2)-EC(1) (likewise ESW(2)-ESW(1)), where EC(1), EC(2), ESW(1), and ESW(2) are defined below. You define the size of this array by allocating space for it (see Allocating Space for Solution-Dependent State Variables for more information).

SERD

Magnitude of the strain energy rate density, $\dot{W}$ (required only in $C_{t}$ -integral calculations). The strain energy rate density is defined as

\dot{W} = \int_{0}^{{\dot{\bar{ε}}}^{c r}} q d {\dot{\bar{ε}}}^{c r} + \int_{0}^{{\dot{\bar{ε}}}^{s w}} p d {\dot{\bar{ε}}}^{s w} = \int_{0}^{\dot{ε}} σ : (d {\dot{ε}}^{c r} + d {\dot{ε}}^{s w}) .

Elastic rates are ignored in the calculation of $\dot{W}$ . The contour integral will, therefore, be path independent only for steady-state creep conditions; that is, when the creep straining dominates throughout the specimen.

Variables Passed in for Information

EC(1)

The definition depends on the usage:

Metal creep and Drucker-Prager creep: ${\bar{ε}}^{c r}$ at the start of the increment.
Capped Drucker-Prager creep: ${\bar{ε}}_{s}^{c r}$ at the start of the increment.
Gasket creep: $ε^{c r}$ at the start of the increment.

EC(2)

The definition depends on the usage:

Metal creep and Drucker-Prager creep: ${\bar{ε}}^{c r}$ at the end of the increment.
Capped Drucker-Prager creep: ${\bar{ε}}_{s}^{c r}$ at the end of the increment.
Gasket creep: $ε^{c r}$ at the end of the increment.

ESW(1)

The definition depends on the usage:

Metal creep: ${\bar{ε}}^{s w}$ at the start of the increment.
Capped Drucker-Prager creep: ${\bar{ε}}_{c}^{c r}$ at the start of the increment.
Drucker-Prager and gasket creep: = 0.

ESW(2)

The definition depends on the usage:

Metal creep: ${\bar{ε}}^{s w}$ at the end of the increment.
Capped Drucker-Prager creep: ${\bar{ε}}_{c}^{c r}$ at the end of the increment.
Drucker-Prager and gasket creep: = 0.

P

The definition depends on the usage:

Metal creep and Drucker-Prager creep: $p = - \frac{1}{3} (σ_{11} + σ_{22} + σ_{33})$ , equivalent pressure stress (in soils analysis this is the equivalent effective pressure stress).
Capped Drucker-Prager creep: ${\bar{p}}^{c r} = p - p_{a}$ , effective creep pressure (in soils analysis p is the effective pressure stress).
Gasket creep: = 0.

If LEND=0, the value is p or ${\bar{p}}^{c r}$ at the beginning of the increment. If LEND=1, the value is p or ${\bar{p}}^{c r}$ at the end of the increment.

QTILD

The definition depends on the usage:

Metal creep: $\tilde{q}$ , Mises or Hill equivalent stress (the Hill formula is used if anisotropic creep is defined; see Anisotropic Creep).
Gasket creep: $σ$ , the uniaxial compressive stress.
Drucker-Prager creep: ${\bar{σ}}^{c r}$ , equivalent creep stress (in soils analysis this is based on effective stresses).
Capped Drucker-Prager creep: ${\bar{σ}}^{c r}$ , equivalent creep stress (in soils analysis this is based on effective stresses).

If LEND=0, the value is $\tilde{q}$ or ${\bar{σ}}^{c r}$ at the beginning of the increment. If LEND=1, the value is $\tilde{q}$ or ${\bar{σ}}^{c r}$ at the end of the increment.

TEMP

Temperature at the end of the increment.

DTEMP

Increment of temperature during the time increment.

PREDEF

An array containing the values of all of the user-specified predefined variables at this point at the end of the increment (initial values at the beginning of the analysis and current values during the analysis).

DPRED

An array containing the increments of all of the predefined variables during the time increment.

TIME(1)

Value of step time at the end of the increment.

TIME(2)

Value of total time at the end of the increment.

TIME(3)

Value of creep time at the end of the increment.

DTIME

Time increment.

CMNAME

User-specified material name or gasket behavior name, left justified. Some internal creep models are given names starting with the “ABQ_” character string. To avoid conflict, you should not use “ABQ_” as the leading string for CMNAME.

LEXIMP

Explicit/implicit flag.

If LEXIMP=0, explicit creep integration is being used and only DECRA(1) and DESWA(1) need be defined; DECRA(I) and DESWA(I), I=2,5, need not be defined.

If LEXIMP=1, implicit creep integration is being used. The derivatives, DECRA(I) and DESWA(I), I=2,5, should be defined accurately to achieve rapid convergence of the solution.

LEND

Start/end of increment flag.

If LEND=0, the routine is being called at the start of the increment. In this case DECRA(1) and DESWA(1) must be defined as the equivalent creep and swelling rates calculated at the beginning of the increment, multiplied by the time increment.

If LEND=1, the routine is being called at the end of the increment. In this case DECRA(1) and DESWA(1) must be defined as the equivalent creep and swelling rates calculated at the end of the increment, multiplied by the time increment. If applicable, the solution-dependent state variables STATEV must be updated as well.

COORDS(3)

An array containing the current coordinates of this point.

NSTATV

Number of solution-dependent state variables associated with this material or gasket behavior type (specified when space is allocated for the array; see Allocating Space for Solution-Dependent State Variables).

NOEL

Element number.

NPT

Integration point number.

LAYER

Layer number (for composite shells and layered solids).

KSPT

Section point number within the current layer.

KSTEP

Step number.

KINC

Increment number.

Example: Hyperbolic Sine Creep Law

Suppose that we wish to model a metal using the creep behavior

\dot{{\bar{ε}}^{c r}} = A {(\sinh \frac{q}{σ_{0}})}^{n}, \dot{{\bar{ε}}^{s w}} = 0,

where A, $σ_{0}$ , and n are constants.

User subroutine CREEP can be coded as follows:

      SUBROUTINE CREEP(DECRA,DESWA,STATEV,SERD,EC,ESW,P,QTILD,
     1 TEMP,DTEMP,PREDEF,DPRED,TIME,DTIME,CMNAME,LEXIMP,LEND,
     2 COORDS,NSTATV,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
C
      INCLUDE 'ABA_PARAM.INC'
C
      CHARACTER*80 CMNAME
C
      DIMENSION DECRA(5),DESWA(5),STATEV(*),PREDEF(*),DPRED(*),
     1 TIME(3),COORDS(*),EC(2),ESW(2)
C
C DEFINE CONSTANTS
C
      A=
      SIG0=
      AN=
C
      T1=EXP(QTILD/SIG0)
      T2=EXP(−QTILD/SIG0)
      DECRA(1) = A*(.5*(T1−T2))**AN*DTIME
      IF(LEXIMP.EQ.1) THEN
       DECRA(5) = AN*A*(.5*(T1−T2))**(AN−1.)*DTIME/
     1             SIG0*.5*(T1+T2)
      END IF
C
      RETURN
      END

The derivative

\frac{\partial △ {\bar{ε}}^{c r}}{\partial q} = \frac{n A Δ t}{σ_{0}} {(\sinh \frac{q}{σ_{0}})}^{n - 1} \cosh \frac{q}{σ_{0}}

has been defined on the assumption that the subroutine will be used with implicit integration.