AbaqusSubroutineDev/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_046.md

<!-- source-page: 451 -->

```txt
bi1 = sInvariant(k,1)
term = bi1-three
Udev(k) = Udev(k)
* + a1*term + a2*term**2 + a3*term**3
duDi(k,1) = a1 + two*a2*term + three*a3*term**2
d2uDiDi(k,indx(1,1)) = two*a2 + three*two*a3*term
C - I2
bi2 = sInvariant(k,2)
term = bi2-three
Udev(k) = Udev(k)
* + b1*term + b2*term**2 + b3*term**3
duDi(k,2) = b1 + two*b2*term + three*b3*term**2
d2uDiDi(k,indx(2,2)) = two*b2 + three*two*b3*term
C - I3 (=J)
bi3 = sInvariant(k,3)
term = bi3-one
duDi(k,3) = two*dinv*term
d2uDiDi(k,indx(3,3)) = two*dinv
C - I4(11)
nI411 = indxInv4(1,1)
bi411 = sInvariant(k,nI411)
term = bi411-one
Udev(k) = Udev(k)
* + c2*term**2 + c3*term**3 + c4*term**4
* + c5*term**5 + c6*term**6
duDi(k,nI411) =
* two*c2*term
* + three*c3*term**2
* + four*c4*term**3
* + five*c5*term**4
* + six*c6*term**5
d2uDiDi(k,indx(nI411,nI411)) =
* two*c2
* + three*two*c3*term
* + four*three*c4*term**2
* + five*four*c5*term**3
* + six*five*c6*term**4
C - I5(11)
nI511 = indxInv5(1,1)
bi511 = sInvariant(k,nI511)
term = bi511-one
Udev(k) = Udev(k) 
```

<!-- source-page: 452 -->

```txt
*    + d2*term**2 + d3*term**3 + d4*term**4
*    + d5*term**5 + d6*term**6
duDi(k,nI511) =
*    two*d2*term
*    + three*d3*term**2
*    + four*d4*term**3
*    + five*d5*term**4
*    + six*d6*term**5
d2uDiDi(k,indx(nI511,nI511)) =
*    two*d2
*    + three*two*d3*term
*    + four*three*d4*term**2
*    + five*four*d5*term**3
*    + six*five*d6*term**4
C - I4(22)
    nI422 = indxInv4(2,2)
    bi422 = sInvariant(k,nI422)
    term = bi422-one
    Udev(k) = Udev(k)
*    + e2*term**2 + e3*term**3 + e4*term**4
*    + e5*term**5 + e6*term**6
duDi(k,nI422) =
*    two*e2*term
*    + three*e3*term**2
*    + four*e4*term**3
*    + five*e5*term**4
*    + six*e6*term**5
d2uDiDi(k,indx(nI422,nI422)) =
*    two*e2
*    + three*two*e3*term
*    + four*three*e4*term**2
*    + five*four*e5*term**3
*    + six*five*e6*term**4
C - I5(22)
    nI522 = indxInv5(2,2)
    bi522 = sInvariant(k,nI522)
    term = bi522-one
    Udev(k) = Udev(k)
*    + f2*term**2 + f3*term**3 + f4*term**4
*    + f5*term**5 + f6*term**6
duDi(k,nI522) =
*    two*f2*term 
```

<!-- source-page: 453 -->

```txt
* + three*f3*term**2
* + four*f4*term**3
* + five*f5*term**4
* + six*f6*term**5
d2uDiDi(k, index(nI522, nI522)) =
* two*f2
* + three*two*f3*term
* + four*three*f4*term**2
* + five*four*f5*term**3
* + six*five*f6*term**4
C - I4(12)
nI412 = indexInv4(1, 2)
bi412 = sInvariant(k, nI412)
term = zeta(k, 1) * (bi412 - zeta(k, 1))
Udev(k) = Udev(k)
* + g2*term**2 + g3*term**3
* + g4*term**4 + g5*term**5
* + g6*term**6
duDi(k, nI412) = zeta(k, 1) * (
* two*g2*term
* + three*g3*term**2
* + four*g4*term**3
* + five*g5*term**4
* + six*g6*term**5 )
d2uDiDi(k, index(nI412, nI412)) = zeta(k, 1) ** 2 * (
* two*g2
* + three*two*g3*term
* + four*three*g4*term**2
* + five*four*g5*term**3
* + six*five*g6*term**4 )
C
end do
C
return
end
C
Function to map index from Square to Triangular storage
C of symmetric matrix
C
integer function index(i, j)
include 'vaba_param.inc'
ii = min(i, j) 
```

<!-- source-page: 454 -->

```matlab
jj = max(i, j)
    index = ii + jj*(jj-1)/2
    return
    end
C
C Function to generate enumeration of scalar
C Pseudo-Invariants of type 4
C integer function indexInv4(i, j)
    include 'vaba_param.inc'
    ii = min(i, j)
    jj = max(i, j)
    indexInv4 = 4 + jj*(jj-1) + 2*(ii-1)
    return
    end
C
C Function to generate enumeration of scalar
C Pseudo-Invariants of type 5
C integer function indexInv5(i, j)
    include 'vaba_param.inc'
    ii = min(i, j)
    jj = max(i, j)
    indexInv5 = 5 + jj*(jj-1) + 2*(ii-1)
    return
    end 
```

# Additional reference

• Kaliske, M., and J. Schmidt, “Formulation of Finite Nonlinear Anisotropic Elasticity,” CADFEM GmbH Infoplaner 2/2005, vol. 2, pp. 22–23, 2005.

<!-- source-page: 455 -->

# 1.2.11 VUANISOHYPER\_STRAIN: User subroutine to define anisotropic hyperelastic material behavior based on Green strain.

Product: Abaqus/Explicit

# References

• “Anisotropic hyperelastic behavior,” Section 22.5.3 of the Abaqus Analysis User’s Guide   
• \*ANISOTROPIC HYPERELASTIC   
• “UANISOHYPER\_INV and VUANISOHYPER\_INV,” Section 4.1.13 of the Abaqus Verification Guide

# Overview

User subroutine VUANISOHYPER\_STRAIN:

• can be used to define the strain energy potential of anisotropic hyperelastic materials as a function of the components of the Green strain tensor;   
• will be called for blocks of material calculation points for which the material definition contains user-defined anisotropic hyperelastic behavior with Green strain-based formulation (“Anisotropic hyperelastic behavior,” Section 22.5.3 of the Abaqus Analysis User’s Guide);   
• can use and update solution-dependent state variables;   
• can use any field variables that are passed in; and   
• requires that the values of the derivatives of the strain energy density function be defined with respect to the components of the modified Green strain tensor and the volume ratio.

# Component ordering in tensors

The component ordering depends upon whether the tensor is second or fourth order.

# Symmetric second-order tensors

For symmetric second-order tensors, such as the modified Green strain tensor, there are ndir+nshr components; the component order is given as a natural permutation of the indices of the tensor. The direct components are first and then the indirect components, beginning with the 12-component. For example, a stress tensor contains ndir direct stress components and nshr shear stress components, which are passed in as

<table><tr><td>Component</td><td>2D Case</td><td>3D Case</td></tr><tr><td>1</td><td> $\overline{\varepsilon}_{11}^{G}$ </td><td> $\overline{\varepsilon}_{11}^{G}$ </td></tr><tr><td>2</td><td> $\overline{\varepsilon}_{22}^{G}$ </td><td> $\overline{\varepsilon}_{22}^{G}$ </td></tr></table>

<!-- source-page: 456 -->

<table><tr><td>Component</td><td>2D Case</td><td>3D Case</td></tr><tr><td>3</td><td> $\overline{\varepsilon}_{33}^{G}$ </td><td> $\overline{\varepsilon}_{33}^{G}$ </td></tr><tr><td>4</td><td> $\overline{\varepsilon}_{12}^{G}$ </td><td> $\overline{\varepsilon}_{12}^{G}$ </td></tr><tr><td>5</td><td></td><td> $\overline{\varepsilon}_{23}^{G}$ </td></tr><tr><td>6</td><td></td><td> $\overline{\varepsilon}_{31}^{G}$ </td></tr></table>

The shear strain components are stored as tensor components and not as engineering components.

# Symmetric fourth-order tensors

For symmetric fourth-order tensors, such as the deviatoric elasticity tensor $\partial ^ { 2 } U / \partial \overline { { \varepsilon } } _ { i j } ^ { G } \partial \overline { { \varepsilon } } _ { k l } ^ { G }$ , there are (ndir+nshr)\*(ndir+nshr+1)/2 independent components. These components are ordered using the following triangular storage scheme:

<table><tr><td>Component</td><td>2D Case</td><td>3D Case</td></tr><tr><td>1</td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{11}^{G}$ </td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{11}^{G}$ </td></tr><tr><td>2</td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{22}^{G}$ </td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{22}^{G}$ </td></tr><tr><td>3</td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{22}^{G}$ </td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{22}^{G}$ </td></tr><tr><td>4</td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{33}^{G}$ </td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{33}^{G}$ </td></tr><tr><td>5</td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{33}^{G}$ </td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{33}^{G}$ </td></tr><tr><td>6</td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial\overline{\varepsilon}_{33}^{G}$ </td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial\overline{\varepsilon}_{33}^{G}$ </td></tr><tr><td>7</td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{12}^{G}$ </td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{12}^{G}$ </td></tr><tr><td>8</td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{12}^{G}$ </td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{12}^{G}$ </td></tr><tr><td>9</td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial\overline{\varepsilon}_{12}^{G}$ </td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial\overline{\varepsilon}_{12}^{G}$ </td></tr><tr><td>10</td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{12}^{G}\partial\overline{\varepsilon}_{12}^{G}$ </td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{12}^{G}\partial\overline{\varepsilon}_{12}^{G}$ </td></tr><tr><td>11</td><td></td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{23}^{G}$ </td></tr><tr><td>12</td><td></td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{23}^{G}$ </td></tr><tr><td>13</td><td></td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial\overline{\varepsilon}_{23}^{G}$ </td></tr><tr><td>14</td><td></td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{12}^{G}\partial\overline{\varepsilon}_{23}^{G}$ </td></tr><tr><td>15</td><td></td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{23}^{G}\partial\overline{\varepsilon}_{23}^{G}$ </td></tr><tr><td>16</td><td></td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{31}^{G}$ </td></tr></table>

<!-- source-page: 457 -->

<table><tr><td>Component</td><td>2D Case</td><td>3D Case</td></tr><tr><td>17</td><td></td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{31}^{G}$ </td></tr><tr><td>18</td><td></td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial\overline{\varepsilon}_{31}^{G}$ </td></tr><tr><td>19</td><td></td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{12}^{G}\partial\overline{\varepsilon}_{31}^{G}$ </td></tr><tr><td>20</td><td></td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{23}^{G}\partial\overline{\varepsilon}_{31}^{G}$ </td></tr><tr><td>21</td><td></td><td> $\partial^{2}U/\partial\overline{\varepsilon}_{31}^{G}\partial\overline{\varepsilon}_{31}^{G}$ </td></tr></table>

If Q denotes the component number of term $\partial ^ { 2 } U / \partial \overline { { \varepsilon } } _ { i j } ^ { G } \partial \overline { { \varepsilon } } _ { k l } ^ { G }$ in the above table and M and N (with $M \ \leq \ N )$ denote the component numbers of $\overline { { \varepsilon } } _ { i j } ^ { G }$ and $\overline { { \varepsilon } } _ { k l } ^ { G }$ , respectively, in the table for second-order tensors, Q is given by the relationship $Q = M ^ { ' } + N \times ( N - 1 ) / 2$ . For example, consider the term $\partial ^ { 2 } U / \partial \overline { { \varepsilon } } _ { 1 1 } ^ { G } \partial \overline { { \varepsilon } } _ { 2 3 } ^ { G }$ . The component numbers for $\overline { { \varepsilon } } _ { 1 1 } ^ { G }$ and $\overline { { \varepsilon } } _ { 2 3 } ^ { G }$ are $M = 1$ and $N = 5$ , respectively, giving $Q = 1 + ( 5 \times 4 ) / 2 = 1 1$ .

# Special consideration for shell elements

When VUANISOHYPER\_STRAIN is used to define the material response of shell elements, Abaqus/Explicit cannot calculate a default value for the transverse shear stiffness of the element. Hence, you must define the element’s transverse shear stiffness. See “Shell section behavior,” Section 29.6.4 of the Abaqus Analysis User’s Guide, for guidelines on choosing this stiffness.

# Material point deletion

Material points that satisfy a user-defined failure criterion can be deleted from the model (see “User-defined mechanical material behavior,” Section 26.7.1 of the Abaqus Analysis User’s Guide). You must specify the state variable number controlling the element deletion flag when you allocate space for the solution-dependent state variables, as explained in “User-defined mechanical material behavior,” Section 26.7.1 of the Abaqus Analysis User’s Guide. The deletion state variable should be set to a value of one or zero in VUANISOHYPER\_STRAIN. A value of one indicates that the material point is active, and a value of zero indicates that Abaqus/Explicit should delete the material point from the model by setting the stresses to zero. The structure of the block of material points passed to user subroutine VUANISOHYPER\_STRAIN remains unchanged during the analysis; deleted material points are not removed from the block. Abaqus/Explicit will “freeze” the values of the strains passed to VUANISOHYPER\_STRAIN for all deleted material points; that is, the strain values remain constant after deletion is triggered. Once a material point has been flagged as deleted, it cannot be reactivated.

# User subroutine interface

subroutine vuanisohyper\_strain( C Read only (unmodifiable) variables – 1 nblock,jElem,kIntPt,kLayer,kSecPt,cmname,

<!-- source-page: 458 -->

```csv
2 ndir,nshr,nstatev,nfieldv,nprops,
3 props,tempOld,tempNew,fieldOld,fieldNew,
4 stateOld, ebar,detu,
C Write only (modifiable) variables -
4 udev,duDe,duDj,
5 d2uDeDe,d2uDjDj,d2uDeDj,
6 stateNew)
C
    include 'vaba_param.inc'
C
    dimension jElem(nblock),
    1 props(nprops),
    2 tempOld(nblock),
    3 fieldOld(nblock,nfieldv),
    4 stateOld(nblock,nstatev),
    5 tempNew(nblock),
    6 fieldNew(nblock,nfieldv),
    7 ebar(nblock,ndir+nshr), detu(nblock),
    8 uDev(nblock),
    9 duDe(nblock,ndir+nshr), duDj(nblock),
    * d2uDeDe(nblock,(ndir+nshr)*(ndir+nshr+1)/2),
    1 d2uDjDj(nblock),
    2 d2uDeDj(nblock,ndir+nshr),
    3 stateNew(nblock,nstatev)
C
    character*80 cmname
C
    do 100 km = 1,nblock
    user coding
100 continue
    return
    end 
```

# Variables to be defined

udev(nblock)

$\tilde { U } _ { d e v ; }$ , the deviatoric part of the strain energy density of the primary material response. This quantity is needed only if the current material definition also includes Mullins effect (see “Mullins effect,” Section 22.6.1 of the Abaqus Analysis User’s Guide).

<!-- source-page: 459 -->

duDe(nblock,ndir+nshr)

Derivatives of strain energy potential with respect to the components of the modified Green strain tensor, ${ \partial U } / { \partial \overline { { \varepsilon } } _ { i j } ^ { G } }$ .

duDj(nblock,ndir+nshr)

Derivatives of strain energy potential with respect to volume ratio, $\partial U / \partial J$

d2uDeDe(nblock,(ndir+nshr)\*(ndir+nshr+1)/2)

Second derivatives of strain energy potential with respect to the components of the modified Green strain tensor (using triangular storage), $\partial ^ { 2 } U / \partial \overline { { \varepsilon } } _ { i j } ^ { G } \partial \overline { { \varepsilon } } _ { k l } ^ { G }$ .

d2uDjDj(nblock)

Second derivatives of strain energy potential with respect to volume ratio, $\partial ^ { 2 } U / \partial J ^ { 2 }$ .

d2uDeDj(nblock,ndir+nshr)

Cross derivatives of strain energy potential with respect to components of the modified Green strain tensor and volume ratio, $\partial ^ { 2 } U / \bar { \partial \varepsilon } _ { i j } ^ { G } \partial J ^ { 2 }$ .

stateNew(nblock,nstatev)

State variables at each material point at the end of the increment. You define the size of this array by allocating space for it (see “User subroutines: overview,” Section 18.1.1 of the Abaqus Analysis User’s Guide, for more information).

# Variables passed in for information

nblock

Number of material points to be processed in this call to VUANISOHYPER\_STRAIN.

jElem(nblock)

Array of element numbers.

kIntPt

Integration point number.

kLayer

Layer number (for composite shells).

kSecPt

Section point number within the current layer.

cmname

User-specified material name, left justified. It is passed in as an uppercase character string. Some internal material models are given names starting with the “ABQ\_” character string. To avoid conflict, you should not use “ABQ\_” as the leading string for cmname.

<!-- source-page: 460 -->

ndir

Number of direct components in a symmetric tensor.

nshr

Number of indirect components in a symmetric tensor.

nstatev

Number of user-defined state variables that are associated with this material type (you define this as described in “Allocating space” in “User subroutines: overview,” Section 18.1.1 of the Abaqus Analysis User’s Guide).

nfieldv

Number of user-defined external field variables.

nprops

User-specified number of user-defined material properties.

props(nprops)

User-supplied material properties.

tempOld(nblock)

Temperatures at each material point at the beginning of the increment.

tempNew(nblock)

Temperatures at each material point at the end of the increment.

fieldOld(nblock,nfieldv)

Values of the user-defined field variables at each material point at the beginning of the increment.

fieldNew(nblock,nfieldv)

Values of the user-defined field variables at each material point at the end of the increment.

stateOld(nblock,nstatev)

State variables at each material point at the beginning of the increment.

ebar(nblock,ndir+nshr)

Modified Green strain tensor, $\overline { { \varepsilon } } ^ { G }$ , at each material point at the end of the increment.

detu(nblock)

J, determinant of deformation gradient (volume ratio) at the end of the increment.

# Example: Using more than one user-defined anisotropic hyperelastic material model

To use more than one user-defined anisotropic hyperelastic material model, the variable cmname can be tested for different material names inside user subroutine VUANISOHYPER\_STRAIN, as illustrated below: