VUANISOHYPER_STRAIN

User subroutine to define anisotropic hyperelastic material behavior based on Green 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);
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.

The following topics are discussed:

Component ordering in tensors
Special consideration for shell elements
Material point deletion
User subroutine interface
Variables to be defined
Variables passed in for information
Example: Using more than one user-defined anisotropic hyperelastic material model
Example: Orthotropic Saint-Venant Kirchhoff model

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

Component	2D Case	3D Case
1	${\bar{ε}}_{11}^{G}$	${\bar{ε}}_{11}^{G}$
2	${\bar{ε}}_{22}^{G}$	${\bar{ε}}_{22}^{G}$
3	${\bar{ε}}_{33}^{G}$	${\bar{ε}}_{33}^{G}$
4	${\bar{ε}}_{12}^{G}$	${\bar{ε}}_{12}^{G}$
5		${\bar{ε}}_{23}^{G}$
6		${\bar{ε}}_{31}^{G}$

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 {\bar{ε}}_{i j}^{G} \partial {\bar{ε}}_{k l}^{G}$ , there are (ndir+nshr)*(ndir+nshr+1)/2 independent components. These components are ordered using the following triangular storage scheme:

Component	2D Case	3D Case
1	$\partial^{2} U / \partial {\bar{ε}}_{11}^{G} \partial {\bar{ε}}_{11}^{G}$	$\partial^{2} U / \partial {\bar{ε}}_{11}^{G} \partial {\bar{ε}}_{11}^{G}$
2	$\partial^{2} U / \partial {\bar{ε}}_{11}^{G} \partial {\bar{ε}}_{22}^{G}$	$\partial^{2} U / \partial {\bar{ε}}_{11}^{G} \partial {\bar{ε}}_{22}^{G}$
3	$\partial^{2} U / \partial {\bar{ε}}_{22}^{G} \partial {\bar{ε}}_{22}^{G}$	$\partial^{2} U / \partial {\bar{ε}}_{22}^{G} \partial {\bar{ε}}_{22}^{G}$
4	$\partial^{2} U / \partial {\bar{ε}}_{11}^{G} \partial {\bar{ε}}_{33}^{G}$	$\partial^{2} U / \partial {\bar{ε}}_{11}^{G} \partial {\bar{ε}}_{33}^{G}$
5	$\partial^{2} U / \partial {\bar{ε}}_{22}^{G} \partial {\bar{ε}}_{33}^{G}$	$\partial^{2} U / \partial {\bar{ε}}_{22}^{G} \partial {\bar{ε}}_{33}^{G}$
6	$\partial^{2} U / \partial {\bar{ε}}_{33}^{G} \partial {\bar{ε}}_{33}^{G}$	$\partial^{2} U / \partial {\bar{ε}}_{33}^{G} \partial {\bar{ε}}_{33}^{G}$
7	$\partial^{2} U / \partial {\bar{ε}}_{11}^{G} \partial {\bar{ε}}_{12}^{G}$	$\partial^{2} U / \partial {\bar{ε}}_{11}^{G} \partial {\bar{ε}}_{12}^{G}$
8	$\partial^{2} U / \partial {\bar{ε}}_{22}^{G} \partial {\bar{ε}}_{12}^{G}$	$\partial^{2} U / \partial {\bar{ε}}_{22}^{G} \partial {\bar{ε}}_{12}^{G}$
9	$\partial^{2} U / \partial {\bar{ε}}_{33}^{G} \partial {\bar{ε}}_{12}^{G}$	$\partial^{2} U / \partial {\bar{ε}}_{33}^{G} \partial {\bar{ε}}_{12}^{G}$
10	$\partial^{2} U / \partial {\bar{ε}}_{12}^{G} \partial {\bar{ε}}_{12}^{G}$	$\partial^{2} U / \partial {\bar{ε}}_{12}^{G} \partial {\bar{ε}}_{12}^{G}$
11		$\partial^{2} U / \partial {\bar{ε}}_{11}^{G} \partial {\bar{ε}}_{23}^{G}$
12		$\partial^{2} U / \partial {\bar{ε}}_{22}^{G} \partial {\bar{ε}}_{23}^{G}$
13		$\partial^{2} U / \partial {\bar{ε}}_{33}^{G} \partial {\bar{ε}}_{23}^{G}$
14		$\partial^{2} U / \partial {\bar{ε}}_{12}^{G} \partial {\bar{ε}}_{23}^{G}$
15		$\partial^{2} U / \partial {\bar{ε}}_{23}^{G} \partial {\bar{ε}}_{23}^{G}$
16		$\partial^{2} U / \partial {\bar{ε}}_{11}^{G} \partial {\bar{ε}}_{31}^{G}$
17		$\partial^{2} U / \partial {\bar{ε}}_{22}^{G} \partial {\bar{ε}}_{31}^{G}$
18		$\partial^{2} U / \partial {\bar{ε}}_{33}^{G} \partial {\bar{ε}}_{31}^{G}$
19		$\partial^{2} U / \partial {\bar{ε}}_{12}^{G} \partial {\bar{ε}}_{31}^{G}$
20		$\partial^{2} U / \partial {\bar{ε}}_{23}^{G} \partial {\bar{ε}}_{31}^{G}$
21		$\partial^{2} U / \partial {\bar{ε}}_{31}^{G} \partial {\bar{ε}}_{31}^{G}$

If Q denotes the component number of term $\partial^{2} U / \partial {\bar{ε}}_{i j}^{G} \partial {\bar{ε}}_{k l}^{G}$ in the above table and M and N (with $M \leq N$ ) denote the component numbers of ${\bar{ε}}_{i j}^{G}$ and ${\bar{ε}}_{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 {\bar{ε}}_{11}^{G} \partial {\bar{ε}}_{23}^{G}$ . The component numbers for ${\bar{ε}}_{11}^{G}$ and ${\bar{ε}}_{23}^{G}$ are $M = 1$ and $N = 5$ , respectively, giving $Q = 1 + (5 \times 4) / 2 = 11$ .

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 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). 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. 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,
     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).
duDe(nblock,ndir+nshr): Derivatives of strain energy potential with respect to the components of the modified Green strain tensor, $\partial U / \partial {\bar{ε}}_{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 {\bar{ε}}_{i j}^{G} \partial {\bar{ε}}_{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 / \partial {\bar{ε}}_{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 About user subroutines and utilities 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.
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).
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, ${\bar{ε}}^{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:

if (cmname(1:4) .eq. 'MAT1') then
   call VUANISOHYPER_STRAIN1(argument_list)
else if (cmname(1:4) .eq. 'MAT2') then
   call VUANISOHYPER_STRAIN2(argument_list)
end if

VUANISOHYPER_STRAIN1 and VUANISOHYPER_STRAIN2 are the actual subroutines containing the anisotropic hyperelastic models for each material MAT1 and MAT2, respectively. Subroutine VUANISOHYPER_STRAIN merely acts as a directory here. The argument list can be the same as that used in subroutine VUANISOHYPER_STRAIN. The material names must be in uppercase characters since cmname is passed in as an uppercase character string.

Example: Orthotropic Saint-Venant Kirchhoff model

As a simple example of the coding of subroutine VUANISOHYPER_STRAIN, consider the generalization to anisotropic hyperelasticity of the Saint-Venant Kirchhoff model. The strain energy function of the Saint-Venant Kirchhoff model can be expressed as a quadratic function of the Green strain tensor, $ε^{G}$ , as

U (ε^{G}) = \frac{1}{2} ε^{G} : D : ε^{G},

where $D$ is the fourth-order elasticity tensor. The derivatives of the strain energy function with respect to the Green strain are given as

\frac{\partial U}{\partial ε^{G}} = D : ε^{G},

\frac{\partial^{2} U}{\partial ε^{G} \partial ε^{G}} = D .

However, subroutine VUANISOHYPER_STRAIN must return the derivatives of the strain energy function with respect to the modified Green strain tensor, ${\bar{ε}}^{G}$ , and the volume ratio, J, which can be accomplished easily using the following relationship between $ε^{G}$ , ${\bar{ε}}^{G}$ , and $J$ :

ε^{G} = J^{\frac{2}{3}} {\bar{ε}}^{G} + \frac{1}{2} (J^{\frac{2}{3}} - 1) I,

where $I$ is the second-order identity tensor. Thus, using the chain rule we find

\frac{\partial U}{\partial {\bar{ε}}^{G}} = J^{\frac{2}{3}} \frac{\partial U}{\partial ε^{G}},

\frac{\partial U}{\partial J} = \frac{\partial ε^{G}}{\partial J} : \frac{\partial U}{\partial ε^{G}},

\frac{\partial^{2} U}{\partial {\bar{ε}}^{G} \partial {\bar{ε}}^{G}} = J^{\frac{4}{3}} \frac{\partial^{2} U}{\partial ε^{G} \partial ε^{G}},

\frac{\partial^{2} U}{\partial J^{2}} = \frac{\partial^{2} ε^{G}}{\partial J^{2}} : \frac{\partial U}{\partial ε^{G}} + \frac{\partial ε^{G}}{\partial J} : \frac{\partial^{2} U}{\partial ε^{G} \partial ε^{G}} : \frac{\partial ε^{G}}{\partial J},

\frac{\partial^{2} U}{\partial {\bar{ε}}^{G} \partial J} = \frac{2}{3 J} J^{\frac{2}{3}} \frac{\partial U}{\partial ε^{G}} + J^{\frac{2}{3}} \frac{\partial^{2} U}{\partial ε^{G} \partial ε^{G}} : \frac{\partial ε^{G}}{\partial J},

where

\frac{\partial ε^{G}}{\partial J} = \frac{2}{3 J} J^{\frac{2}{3}} ({\bar{ε}}^{G} + \frac{1}{2} I) = \frac{2}{3 J} (ε^{G} + \frac{1}{2} I)

and

\frac{\partial^{2} ε^{G}}{\partial J^{2}} = - \frac{1}{3 J} \frac{\partial ε^{G}}{\partial J} .

In this example an auxiliary function is used to facilitate indexing into a fourth-order symmetric tensor. The subroutine would be coded as follows:

      subroutine vuanisohyper_strain (
C Read only -
     *     nblock, 
     *     jElem, kIntPt, kLayer, kSecPt, 
     *     cmname,
     *     ndir, nshr, nstatev, nfieldv, nprops,
     *     props, tempOld, tempNew, fieldOld, fieldNew,
     *     stateOld, ebar, detu,
C Write only -
     *     uDev, duDe, duDj,
     *     d2uDeDe, d2uDjDj, d2uDeDj,
     *     stateNew )
C
      include 'vaba_param.inc'
C
      dimension props(nprops), 
     *  tempOld(nblock),
     *  fieldOld(nblock,nfieldv), 
     *  stateOld(nblock,nstatev), 
     *  tempNew(nblock),
     *  fieldNew(nblock,nfieldv),
     *  ebar(nblock,ndir+nshr), detu(nblock),
     *  uDev(nblock), duDe(nblock,ndir+nshr), duDj(nblock),
     *  d2uDeDe(nblock,*), d2uDjDj(nblock), 
     *  d2uDeDj(nblock,ndir+nshr),
     *  stateNew(nblock,nstatev)
C
      character*80 cmname
C
      parameter( half = 0.5d0, one = 1.d0, two = 2.d0, 
     *     third = 1.d0/3.d0, twothds = 2.d0/3.d0, four = 4.d0,
     *     dinv = 0.d0 )
C
C    Orthotropic Saint-Venant Kirchhoff strain energy function 
C    (3D)
C
      D1111 = props(1)
      D1122 = props(2)
      D2222 = props(3)
      D1133 = props(4)
      D2233 = props(5)
      D3333 = props(6)
      D1212 = props(7)
      D1313 = props(8)
      D2323 = props(9)
C
      do k = 1, nblock
C
         d2UdE11dE11 = D1111
         d2UdE11dE22 = D1122
         d2UdE11dE33 = D1133
         d2UdE22dE11 = d2UdE11dE22
         d2UdE22dE22 = D2222
         d2UdE22dE33 = D2233
         d2UdE33dE11 = d2UdE11dE33
         d2UdE33dE22 = d2UdE22dE33
         d2UdE33dE33 = D3333
         d2UdE12dE12 = D1212
         d2UdE13dE13 = D1313
         d2UdE23dE23 = D2323
C
         xpow = exp ( log(detu(k)) * twothds )
         detuInv = one / detu(k)
C
         E11 = xpow * ebar(k,1) + half * ( xpow - one )
         E22 = xpow * ebar(k,2) + half * ( xpow - one )
         E33 = xpow * ebar(k,3) + half * ( xpow - one )
         E12 = xpow * ebar(k,4)
         E23 = xpow * ebar(k,5)
         E13 = xpow * ebar(k,6)
C
         term1 = twothds * detuInv
         dE11Dj = term1 * ( E11 + half )
         dE22Dj = term1 * ( E22 + half )
         dE33Dj = term1 * ( E33 + half )
         dE12Dj = term1 * E12
         dE23Dj = term1 * E23
         dE13Dj = term1 * E13
         term2 = - third * detuInv
         d2E11DjDj = term2 * dE11Dj
         d2E22DjDj = term2 * dE22Dj
         d2E33DjDj = term2 * dE33Dj
         d2E12DjDj = term2 * dE12Dj
         d2E23DjDj = term2 * dE23Dj
         d2E13DjDj = term2 * dE13Dj
C
         dUdE11 = d2UdE11dE11 * E11 
     *          + d2UdE11dE22 * E22 
     *          + d2UdE11dE33 * E33 
         dUdE22 = d2UdE22dE11 * E11 
     *          + d2UdE22dE22 * E22 
     *          + d2UdE22dE33 * E33 
         dUdE33 = d2UdE33dE11 * E11 
     *          + d2UdE33dE22 * E22 
     *          + d2UdE33dE33 * E33 
         dUdE12 = two * d2UdE12dE12 * E12 
         dUdE23 = two * d2UdE23dE23 * E23
         dUdE13 = two * d2UdE13dE13 * E13 
C
         U = half * ( E11*dUdE11 + E22*dUdE22 + E33*dUdE33 )
     *        + E12*dUdE12 + E13*dUdE13 + E23*dUdE23 
         uDev(k) = U
C
         duDe(k,1) = xpow * dUdE11 
         duDe(k,2) = xpow * dUdE22
         duDe(k,3) = xpow * dUdE33 
         duDe(k,4) = xpow * dUdE12 
         duDe(k,5) = xpow * dUdE23 
         duDe(k,6) = xpow * dUdE13 
C
         xpow2 = xpow * xpow
C	Only update nonzero components 
         d2uDeDe(k,indx(1,1)) = xpow2 * d2UdE11dE11
         d2uDeDe(k,indx(1,2)) = xpow2 * d2UdE11dE22
         d2uDeDe(k,indx(2,2)) = xpow2 * d2UdE22dE22
         d2uDeDe(k,indx(1,3)) = xpow2 * d2UdE11dE33
         d2uDeDe(k,indx(2,3)) = xpow2 * d2UdE22dE33
         d2uDeDe(k,indx(3,3)) = xpow2 * d2UdE33dE33
         d2uDeDe(k,indx(4,4)) = xpow2 * d2UdE12dE12
         d2uDeDe(k,indx(5,5)) = xpow2 * d2UdE23dE23
         d2uDeDe(k,indx(6,6)) = xpow2 * d2UdE13dE13
C
         duDj(k) = dUdE11*dE11Dj + dUdE22*dE22Dj + dUdE33*dE33Dj 
     *        + two * ( dUdE12*dE12Dj + dUdE13*dE13Dj 
     *                + dUdE23*dE23Dj )        
         d2uDjDj(k)= dUdE11*d2E11DjDj+dUdE22*d2E22DjDj
     *               +dUdE33*d2E33DjDj 
     *        + two*(dUdE12*d2E12DjDj+dUdE13*d2E13DjDj
     *               +dUdE23*d2E23DjDj)
     *        + d2UdE11dE11 * dE11Dj * dE11Dj
     *        + d2UdE22dE22 * dE22Dj * dE22Dj
     *        + d2UdE33dE33 * dE33Dj * dE33Dj
     *        + two  * ( d2UdE11dE22 * dE11Dj * dE22Dj
     *                 + d2UdE11dE33 * dE11Dj * dE33Dj
     *                 + d2UdE22dE33 * dE22Dj * dE33Dj )
     *        + four * ( d2UdE12dE12 * dE12Dj * dE12Dj
     *                   d2UdE13dE13 * dE13Dj * dE13Dj
     *                   d2UdE23dE23 * dE23Dj * dE23Dj )
C
         d2uDeDj(k,1) = xpow * ( term1 * dUdE11 
     *                + d2UdE11dE11 * dE11Dj
     *                + d2UdE11dE22 * dE22Dj
     *                + d2UdE11dE33 * dE33Dj ) 
         d2uDeDj(k,2) = xpow * ( term1 * dUdE22 
     *                + d2UdE22dE11 * dE11Dj
     *                + d2UdE22dE22 * dE22Dj
     *                + d2UdE22dE33 * dE33Dj )
         d2uDeDj(k,3) = xpow * ( term1 * dUdE33 
     *                + d2UdE33dE11 * dE11Dj
     *                + d2UdE33dE22 * dE22Dj
     *                + d2UdE33dE33 * dE33Dj )
         d2uDeDj(k,4) = xpow * ( term1 * dUdE12 
     *                + two * d2UdE12dE12 * dE12Dj )
         d2uDeDj(k,5) = xpow * ( term1 * dUdE23 
     *                + two * d2UdE23dE23 * dE23Dj )
         d2uDeDj(k,6) = xpow * ( term1 * dUdE13 
     *                + two * d2UdE13dE13 * dE13Dj )
      end do
C
      return
      end
C
      integer function indx( i, j )
C
      include 'vaba_param.inc'
C
C     Function to map index from Square to Triangular storage 
C 		  of symmetric matrix
C
      ii = min(i,j)
      jj = max(i,j)
C
      indx = ii + jj*(jj-1)/2
C
      return
      end