VUANISOHYPER_INV

User subroutine to define anisotropic hyperelastic material behavior using the invariant formulation.

can be used to define the strain energy potential of anisotropic hyperelastic materials as a function of an irreducible set of scalar invariants;
will be called for blocks of material calculation points for which the material definition contains user-defined anisotropic hyperelastic behavior with invariant-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 scalar invariants.

The following topics are discussed:

Enumeration of invariants
Storage of arrays of derivatives of energy function
Special considerations 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: Anisotropic hyperelastic model of Kaliske and Schmidt
Additional reference

Enumeration of invariants

To facilitate coding and provide easy access to the array of invariants passed to user subroutine VUANISOHYPER_INV, an enumerated representation of each invariant is introduced. Any scalar invariant can, therefore, be represented uniquely by an enumerated invariant, $I_{n}^{*}$ , where the subscript n denotes the order of the invariant according to the enumeration scheme in the following table:

Invariant	Enumeration, $n$
${\bar{I}}_{1}$	$1$
${\bar{I}}_{2}$	$2$
$J$	$3$
${\bar{I}}_{4 (α β)}$	$4 + 2 (α - 1) + β (β - 1); α \leq β$
${\bar{I}}_{5 (α β)}$	$5 + 2 (α - 1) + β (β - 1); α \leq β$

For example, in the case of three families of fibers there are a total of 15 invariants: ${\bar{I}}_{1}$ , ${\bar{I}}_{2}$ , $J$ , six invariants of type ${\bar{I}}_{4 (α β)}$ , and six invariants of type ${\bar{I}}_{5 (α β)}$ , with $α, β = 1, 2, 3 (α \leq β)$ . The following correspondence exists between each of these invariants and their enumerated counterpart:

Enumerated invariant	Invariant
$I_{1}^{*}$	${\bar{I}}_{1}$
$I_{2}^{*}$	${\bar{I}}_{2}$
$I_{3}^{*}$	$J$
$I_{4}^{*}$	${\bar{I}}_{4 (11)}$
$I_{5}^{*}$	${\bar{I}}_{5 (11)}$
$I_{6}^{*}$	${\bar{I}}_{4 (12)}$
$I_{7}^{*}$	${\bar{I}}_{5 (12)}$
$I_{8}^{*}$	${\bar{I}}_{4 (22)}$
$I_{9}^{*}$	${\bar{I}}_{5 (22)}$
$I_{10}^{*}$	${\bar{I}}_{4 (13)}$
$I_{11}^{*}$	${\bar{I}}_{5 (13)}$
$I_{12}^{*}$	${\bar{I}}_{4 (23)}$
$I_{13}^{*}$	${\bar{I}}_{5 (23)}$
$I_{14}^{*}$	${\bar{I}}_{4 (33)}$
$I_{15}^{*}$	${\bar{I}}_{5 (33)}$

A similar scheme is used for the array zeta of terms $ζ_{α β} = A_{α} \cdot A_{β}$ . Each term can be represented uniquely by an enumerated counterpart $ζ_{m}^{*}$ , as shown below:

Dot product	Enumeration, $m$
$ζ_{α β}$	$α + \frac{1}{2} (β - 2) (β - 1); α < β$

As an example, for the case of three families of fibers there are three $ζ_{α β}$ terms: $ζ_{12}$ , $ζ_{13}$ , and $ζ_{23}$ . These are stored in the zeta array as $(ζ_{1}^{*}, ζ_{2}^{*}, ζ_{3}^{*})$ .

Storage of arrays of derivatives of energy function

The components of the array duDi of first derivatives of the strain energy potential with respect to the scalar invariants, $\partial U / \partial I_{i}^{*}$ , are stored using the enumeration scheme discussed above for the scalar invariants.

The elements of the array d2uDiDi of second derivatives of the strain energy function, $\partial^{2} U / \partial I_{i}^{*} \partial I_{j}^{*}$ , are laid out in memory using triangular storage: if $k$ denotes the component in this array corresponding to the term $\partial^{2} U / \partial I_{i}^{*} \partial I_{j}^{*}$ , then $k = i + j \times (j - 1) / 2; (i \leq j)$ . For example, the term $\partial^{2} U / \partial I_{2}^{*} \partial I_{5}^{*}$ is stored in component $k = 2 + (5 \times 4) / 2 = 12$ in the d2uDiDi array.

Special considerations for shell elements

When VUANISOHYPER_INV 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_INV. 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_INV remains unchanged during the analysis; deleted material points are not removed from the block. Abaqus/Explicit will “freeze” the values of the invariants passed to VUANISOHYPER_INV for all deleted material points; that is, the 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_inv (
C Read only (unmodifiable) variables –
     1     nblock, nFiber, nInv, 
     2     jElem, kIntPt, kLayer, kSecPt, 
     3     cmname,
     4     nstatev, nfieldv, nprops,
     5     props, tempOld, tempNew, fieldOld, fieldNew,
     6     stateOld, sInvariant, zeta,  
C Write only (modifiable) variables –
     7     uDev, duDi, d2uDiDi,
     8     stateNew )
C
      include 'vaba_param.inc'
C
      dimension props(nprops), 
     1  tempOld(nblock),
     2  fieldOld(nblock,nfieldv), 
     3  stateOld(nblock,nstatev), 
     4  tempNew(nblock),
     5  fieldNew(nblock,nfieldv),
     6  sInvariant(nblock,nInv), 
     7  zeta(nblock,nFiber*(nFiber-1)/2),
     8  uDev(nblock), duDi(nblock,nInv), 
     9  d2uDiDi(nblock,nInv*(nInv+1)/2),
     *  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).
duDi(nblock,nInv): Array of derivatives of strain energy potential with respect to the scalar invariants, $\partial U / \partial I_{i}^{*}$ , ordered using the enumeration scheme discussed above.
d2uDiDi(nblock,nInv*(nInv+1)/2): Arrays of second derivatives of strain energy potential with respect to the scalar invariants (using triangular storage), $\partial^{2} U / \partial I_{i}^{*} \partial I_{j}^{*}$ .
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.
nFiber: Number of families of fibers defined for this material.
nInv: Number of scalar invariants.
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.
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.
sInvariant(nblock,nInv): Array of scalar invariants, $I_{i}^{*}$ , at each material point at the end of the increment. The invariants are ordered using the enumeration scheme discussed above.
zeta(nblock,nFiber*(nFiber-1)/2) ): Array of dot product between the directions of different families of fiber in the reference configuration, $ζ_{α β} = A_{α} \cdot A_{β}$ . The array contains the enumerated values $ζ_{m}^{*}$ using the scheme discussed above.

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_INV, as illustrated below:

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

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

Example: Anisotropic hyperelastic model of Kaliske and Schmidt

As an example of the coding of subroutine VUANISOHYPER_INV, consider the model proposed by Kaliske and Schmidt (2005) for nonlinear anisotropic elasticity with two families of fibers. The strain energy function is given by a polynomial series expansion in the form

\begin{matrix} U = & \frac{1}{D} {(J - 1)}^{2} + \sum_{i = 1}^{3} a_{i} {({\bar{I}}_{1} - 3)}^{i} + \sum_{j = 1}^{3} b_{j} {({\bar{I}}_{2} - 3)}^{j} + \sum_{k = 2}^{6} c_{k} {({\bar{I}}_{4 (11)} - 1)}^{k} + \sum_{l = 2}^{6} d_{l} {({\bar{I}}_{5 (11)} - 1)}^{l} \\ + \sum_{m = 2}^{6} e_{m} {({\bar{I}}_{4 (22)} - 1)}^{m} + \sum_{n = 2}^{6} f_{n} {({\bar{I}}_{5 (22)} - 1)}^{n} + \sum_{p = 2}^{6} g_{p} {(ζ_{12} {\bar{I}}_{4 (12)} - ζ_{12}^{2})}^{p} . \end{matrix}

The code in subroutine VUANISOHYPER_INV must return the derivatives of the strain energy function with respect to the scalar invariants, which are readily computed from the above expression. In this example auxiliary functions are used to facilitate enumeration of pseudo-invarinats of type ${\bar{I}}_{4 (α β)}$ and ${\bar{I}}_{5 (α β)}$ , as well as for indexing into the array of second derivatives using symmetric storage. The subroutine would be coded as follows:

      subroutine vuanisohyper_inv (
C Read only -
     *     nblock, nFiber, nInv, 
     *     jElem, kIntPt, kLayer, kSecPt, 
     *     cmname,
     *     nstatev, nfieldv, nprops,
     *     props, tempOld, tempNew, fieldOld, fieldNew,
     *     stateOld, sInvariant, zeta,  
C     Write only -
     *     uDev, duDi, d2uDiDi,
     *     stateNew )
C
      include 'vaba_param.inc'
C
      dimension props(nprops), 
     *  tempOld(nblock),
     *  fieldOld(nblock,nfieldv), 
     *  stateOld(nblock,nstatev), 
     *  tempNew(nblock),
     *  fieldNew(nblock,nfieldv),
     *  sInvariant(nblock,nInv), 
     *  zeta(nblock,nFiber*(nFiber-1)/2),
     *  uDev(nblock), duDi(nblock,*), 
     *  d2uDiDi(nblock,*),
     *  stateNew(nblock,nstatev)
C
      character*80 cmname
C
      parameter ( zero = 0.d0, one = 1.d0, two = 2.d0, 
     *     three = 3.d0, four = 4.d0, five = 5.d0, six = 6.d0 )
C
C     Kaliske energy function (3D)
C
C     Read material properties
      d=props(1)
      dinv = one / d
      a1=props(2)
      a2=props(3)
      a3=props(4)
      b1=props(5)
      b2=props(6)
      b3=props(7)
      c2=props(8)
      c3=props(9)
      c4=props(10)
      c5=props(11)
      c6=props(12)
      d2=props(13)
      d3=props(14)
      d4=props(15)
      d5=props(16)
      d6=props(17)
      e2=props(18)
      e3=props(19)
      e4=props(20)
      e5=props(21)
      e6=props(22)
      f2=props(23)
      f3=props(24)
      f4=props(25)
      f5=props(26)
      f6=props(27)
      g2=props(28)
      g3=props(29)
      g4=props(30)
      g5=props(31)
      g6=props(32)
C
      do k = 1, nblock
         Udev(k) = zero
C     Compute Udev and 1st and 2nd derivatives w.r.t invariants
C     - I1
         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) 
     *        + 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
     *        + three*f3*term**2
     *        + four*f4*term**3
     *        + five*f5*term**4
     *        + six*f6*term**5
         d2uDiDi(k,indx(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 = indxInv4(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,indx(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
C     Function to map index from Square to Triangular storage 
C     of symmetric matrix
C
      integer function indx( i, j )
      include 'vaba_param.inc'
      ii = min(i,j)
      jj = max(i,j)
      indx = ii + jj*(jj-1)/2
      return
      end
C
C     Function to generate enumeration of scalar
C     Pseudo-Invariants of type 4

C     integer function indxInv4( i, j )
      include 'vaba_param.inc'
      ii = min(i,j)
      jj = max(i,j)
      indxInv4 = 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 indxInv5( i, j )
      include 'vaba_param.inc'
      ii = min(i,j)
      jj = max(i,j)
      indxInv5 = 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.