UANISOHYPER_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;
is called at all material calculation points of elements for which the material definition contains user-defined anisotropic hyperelastic behavior with a Green strain-based formulation (Anisotropic hyperelastic behavior);
can include material behavior dependent on field variables or state variables;
requires that the values of the derivatives of the strain energy density function of the anisotropic hyperelastic material be defined with respect to the components of the modified Green strain tensor; and
is called twice per material point in each iteration.

The following topics are discussed:

Storage of strain components
Storage of arrays of derivatives of the energy function
Special considerations for various element types
User subroutine interface
Variables to be defined
Variables passed in for information
Example: Orthotropic Saint-Venant Kirchhoff model

Storage of strain components

In the array of modified Green strain, EBAR, direct components are stored first, followed by shear components. There are NDI direct and NSHR tensor shear components. The order of the components is defined in Conventions. Since the number of active stress and strain components varies between element types, the routine must be coded to provide for all element types with which it will be used.

Storage of arrays of derivatives of the energy function

The array of first derivatives of the strain energy function, DU1, contains NTENS+1 components, with NTENS=NDI+NSHR. The first NTENS components correspond to the derivatives with respect to each component of the modified Green strain, $\partial U / \partial {\bar{ε}}_{i j}^{G}$ . The last component contains the derivative with respect to the volume ratio, $\partial U / \partial J$ .

The array of second derivatives of the strain energy function, DU2, contains (NTENS+1)*(NTENS+2)/2 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 J$	$\partial^{2} U / \partial {\bar{ε}}_{11}^{G} \partial {\bar{ε}}_{13}^{G}$
12	$\partial^{2} U / \partial {\bar{ε}}_{22}^{G} \partial J$	$\partial^{2} U / \partial {\bar{ε}}_{22}^{G} \partial {\bar{ε}}_{13}^{G}$
13	$\partial^{2} U / \partial {\bar{ε}}_{33}^{G} \partial J$	$\partial^{2} U / \partial {\bar{ε}}_{33}^{G} \partial {\bar{ε}}_{13}^{G}$
14	$\partial^{2} U / \partial {\bar{ε}}_{12}^{G} \partial J$	$\partial^{2} U / \partial {\bar{ε}}_{12}^{G} \partial {\bar{ε}}_{13}^{G}$
15	$\partial^{2} U / \partial J^{2}$	$\partial^{2} U / \partial {\bar{ε}}_{13}^{G} \partial {\bar{ε}}_{13}^{G}$
16		$\partial^{2} U / \partial {\bar{ε}}_{11}^{G} \partial {\bar{ε}}_{23}^{G}$
17		$\partial^{2} U / \partial {\bar{ε}}_{22}^{G} \partial {\bar{ε}}_{23}^{G}$
18		$\partial^{2} U / \partial {\bar{ε}}_{33}^{G} \partial {\bar{ε}}_{23}^{G}$
19		$\partial^{2} U / \partial {\bar{ε}}_{12}^{G} \partial {\bar{ε}}_{23}^{G}$
20		$\partial^{2} U / \partial {\bar{ε}}_{13}^{G} \partial {\bar{ε}}_{23}^{G}$
21		$\partial^{2} U / \partial {\bar{ε}}_{23}^{G} \partial {\bar{ε}}_{23}^{G}$
22		$\partial^{2} U / \partial {\bar{ε}}_{11}^{G} \partial J$
23		$\partial^{2} U / \partial {\bar{ε}}_{22}^{G} \partial J$
24		$\partial^{2} U / \partial {\bar{ε}}_{33}^{G} \partial J$
25		$\partial^{2} U / \partial {\bar{ε}}_{12}^{G} \partial J$
26		$\partial^{2} U / \partial {\bar{ε}}_{13}^{G} \partial J$
27		$\partial^{2} U / \partial {\bar{ε}}_{23}^{G} \partial J$
28		$\partial^{2} U / \partial J^{2}$

Finally, the array of third derivatives of the strain energy function, DU3, also contains (NTENS+1)*(NTENS+2)/2 components, each representing the derivative with respect to $J$ of the corresponding component of DU2. It follows the same triangular storage scheme as DU2.

Special considerations for various element types

There are several special considerations that need to be noted.

Shells that calculate transverse shear energy

When UANISOHYPER_STRAIN is used to define the material response of shell elements that calculate transverse shear energy, Abaqus/Standard 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.

Elements with hourglassing modes

When UANISOHYPER_STRAIN is used to define the material response of elements with hourglassing modes, you must define the hourglass stiffness for hourglass control based on the total stiffness approach. The hourglass stiffness is not required for enhanced hourglass control, but you can define a scaling factor for the stiffness associated with the drill degree of freedom (rotation about the surface normal). See Section controls.

User subroutine interface

      SUBROUTINE UANISOHYPER_STRAIN (EBAR, AJ, UA, DU1, DU2, DU3,
     1 TEMP, NOEL, CMNAME, INCMPFLAG, IHYBFLAG, NDI, NSHR, NTENS,
     2 NUMSTATEV, STATEV, NUMFIELDV, FIELDV, FIELDVINC,
     3 NUMPROPS, PROPS)
C
 	     INCLUDE 'ABA_PARAM.INC'  
C
      CHARACTER*80 CMNAME
C
      DIMENSION EBAR(NTENS), UA(2), DU1(NTENS+1),
     2     DU2((NTENS+1)*(NTENS+2)/2),
     3     DU3((NTENS+1)*(NTENS+2)/2),
     4     STATEV(NUMSTATEV), FIELDV(NUMFIELDV),
     5     FIELDVINC(NUMFIELDV), PROPS(NUMPROPS)


				user coding to define UA,DU1,DU2,DU3,STATEV


      RETURN
      END

Variables to be defined

UA(1): U, strain energy density function. For a compressible material at least one derivative involving J should be nonzero. For an incompressible material all derivatives involving J are ignored.
UA(2): ${\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).
DU1(NTENS+1): Derivatives of strain energy potential with respect to the components of the modified Green strain tensor, $\partial U / \partial {\bar{ε}}_{i j}^{G}$ , and with respect to the volume ratio, $\partial U / \partial J$ .
DU2((NTENS+1)*(NTENS+2)/2): Second derivatives of strain energy potential with respect to the components of the modified Green strain tensor and the volume ratio (using triangular storage, as mentioned earlier).
DU3((NTENS+1)*(NTENS+2)/2): Derivatives with respect to J of the second derivatives of the strain energy potential (using triangular storage, as mentioned earlier). This quantity is needed only for compressible materials with a hybrid formulation (when INCMPFLAG = 0 and IHYBFLAG = 1).
STATEV: Array containing the user-defined solution-dependent state variables at this point. These are supplied as values at the start of the increment or as values updated by other user subroutines (see About user subroutines and utilities) and must be returned as values at the end of the increment.

Variables passed in for information

TEMP: Current temperature at this point.
NOEL: Element number.
CMNAME: User-specified material name, left justified.
NDI: Number of direct stress components at this point.
NSHR: Number of shear components at this point.
NTENS: Size of the stress or strain component array (NDI + NSHR).
INCMPFLAG: Incompressibility flag defined to be 1 if the material is specified as incompressible or 0 if the material is specified as compressible.
IHYBFLAG: Hybrid formulation flag defined to be 1 for hybrid elements; 0 otherwise.
NUMSTATEV: User-defined number of solution-dependent state variables associated with this material (see Allocating space).
NUMFIELDV: Number of field variables.
FIELDV: Array of interpolated values of predefined field variables at this material point at the end of the increment based on the values read in at the nodes (initial values at the beginning of the analysis and current values during the analysis).
FIELDVINC: Array of increments of predefined field variables at this material point for this increment, including any values updated by user subroutine USDFLD.
NUMPROPS: Number of material properties entered for this user-defined hyperelastic material.
PROPS: Array of material properties entered for this user-defined hyperelastic material.
EBAR(NTENS): Modified Green strain tensor, ${\bar{ε}}^{G}$ , at the material point at the end of the increment.
AJ: J, determinant of deformation gradient (volume ratio) at the end of the increment.

Example: Orthotropic Saint-Venant Kirchhoff model

As a simple example of the coding of user subroutine UANISOHYPER_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, user subroutine UANISOHYPER_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 user subroutine would be coded as follows:

      subroutine uanisohyper_strain (
     *     ebar, aj, ua, du1, du2, du3, temp, noel, cmname,
     *     incmpFlag, ihybFlag, ndi, nshr, ntens,
     *     numStatev, statev, numFieldv, fieldv, fieldvInc,
     *     numProps, props)
c
      include 'aba_param.inc'
c
      dimension ebar(ntens), ua(2), du1(ntens+1)
      dimension du2((ntens+1)*(ntens+2)/2)
      dimension du3((ntens+1)*(ntens+2)/2)
      dimension statev(numStatev), fieldv(numFieldv)
      dimension fieldvInc(numFieldv), props(numProps)
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 )
*
* Orthotropic Saint-Venant Kirchhoff strain energy function (3D)
*
      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)
*
      d2UdE11dE11 = D1111
      d2UdE11dE22 = D1122
      d2UdE11dE33 = D1133
*
      d2UdE22dE11 = d2UdE11dE22
      d2UdE22dE22 = D2222
      d2UdE22dE33 = D2233
*
      d2UdE33dE11 = d2UdE11dE33
      d2UdE33dE22 = d2UdE22dE33
      d2UdE33dE33 = D3333
*
      d2UdE12dE12 = D1212
*
      d2UdE13dE13 = D1313
*
      d2UdE23dE23 = D2323
*
      xpow = exp ( log(aj) * twothds )
      detuInv = one / aj
*
      E11 = xpow * ebar(1) + half * ( xpow - one )
      E22 = xpow * ebar(2) + half * ( xpow - one )
      E33 = xpow * ebar(3) + half * ( xpow - one )
      E12 = xpow * ebar(4)
      E13 = xpow * ebar(5)
      E23 = xpow * ebar(6)
*
      term1 = twothds * detuInv
      dE11Dj = term1 * ( E11 + half )
      dE22Dj = term1 * ( E22 + half )
      dE33Dj = term1 * ( E33 + half )
      dE12Dj = term1 * E12
      dE13Dj = term1 * E13
      dE23Dj = term1 * E23
      term2 = - third * detuInv
      d2E11DjDj = term2 * dE11Dj
      d2E22DjDj = term2 * dE22Dj
      d2E33DjDj = term2 * dE33Dj
      d2E12DjDj = term2 * dE12Dj
      d2E13DjDj = term2 * dE13Dj
      d2E23DjDj = term2 * dE23Dj
*
      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
      dUdE13 = two * d2UdE13dE13 * E13
      dUdE23 = two * d2UdE23dE23 * E23
*
      U = half * ( E11*dUdE11 + E22*dUdE22 + E33*dUdE33 )
     *           + E12*dUdE12 + E13*dUdE13 + E23*dUdE23 
*
      ua(2) = U
      ua(1) = ua(2)
*
      du1(1) = xpow * dUdE11 
      du1(2) = xpow * dUdE22
      du1(3) = xpow * dUdE33 
      du1(4) = xpow * dUdE12 
      du1(5) = xpow * dUdE13 
      du1(6) = xpow * dUdE23 
      du1(7) = dUdE11*dE11Dj + dUdE22*dE22Dj + dUdE33*dE33Dj 
     *     + two * ( dUdE12*dE12Dj
     *              +dUdE13*dE13Dj
     *              +dUdE23*dE23Dj )        
*
      xpow2 = xpow * xpow
*
      du2(indx(1,1)) = xpow2 * d2UdE11dE11
      du2(indx(1,2)) = xpow2 * d2UdE11dE22
      du2(indx(2,2)) = xpow2 * d2UdE22dE22
      du2(indx(1,3)) = xpow2 * d2UdE11dE33
      du2(indx(2,3)) = xpow2 * d2UdE22dE33
      du2(indx(3,3)) = xpow2 * d2UdE33dE33
      du2(indx(1,4)) = zero
      du2(indx(2,4)) = zero
      du2(indx(3,4)) = zero
      du2(indx(4,4)) = xpow2 * d2UdE12dE12
      du2(indx(1,5)) = zero
      du2(indx(2,5)) = zero
      du2(indx(3,5)) = zero
      du2(indx(4,5)) = zero
      du2(indx(5,5)) = xpow2 * d2UdE13dE13
      du2(indx(1,6)) = zero
      du2(indx(2,6)) = zero
      du2(indx(3,6)) = zero
      du2(indx(4,6)) = zero
      du2(indx(5,6)) = zero
      du2(indx(6,6)) = xpow2 * d2UdE23dE23
*
      du2(indx(1,7)) = xpow * ( term1 * dUdE11 
     *     + d2UdE11dE11 * dE11Dj
     *     + d2UdE11dE22 * dE22Dj
     *     + d2UdE11dE33 * dE33Dj ) 
      du2(indx(2,7)) = xpow * ( term1 * dUdE22 
     *     + d2UdE22dE11 * dE11Dj
     *     + d2UdE22dE22 * dE22Dj
     *     + d2UdE22dE33 * dE33Dj )
      du2(indx(3,7)) = xpow * ( term1 * dUdE33 
     *     + d2UdE33dE11 * dE11Dj
     *     + d2UdE33dE22 * dE22Dj
     *     + d2UdE33dE33 * dE33Dj )
      du2(indx(4,7)) = xpow * ( term1 * dUdE12 
     *     + two * d2UdE12dE12 * dE12Dj )
      du2(indx(5,7)) = xpow * ( term1 * dUdE13 
     *     + two * d2UdE13dE13 * dE23Dj )
      du2(indx(6,7)) = xpow * ( term1 * dUdE23 
     *     + two * d2UdE23dE23 * dE13Dj )
      du2(indx(7,7))= 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 )
*
      return
      end
*
* Maps index from Square to Triangular storage
* of symmetric matrix
*
      integer function indx( i, j )
*
      include 'aba_param.inc'
*
      ii = min(i,j)
      jj = max(i,j)
*
      indx = ii + jj*(jj-1)/2
*
      return
      end