Product: Abaqus/Explicit
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.
The component ordering depends upon whether the tensor is second or fourth order.
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
The shear strain components are stored as tensor components and not as engineering components.For symmetric fourth-order tensors, such as the deviatoric elasticity tensor 
, there are (ndir+nshr)*(ndir+nshr+1)/2 independent components. These components are ordered using the following triangular storage scheme:
If Q denotes the component number of term in the above table and M and N (with 
) denote the component numbers of 
and 
, respectively, in the table for second-order tensors, Q is given by the relationship 
. For example, consider the term 
. The component numbers for 
and 
are 
and 
, respectively, giving 
.
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 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.
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
endudev(nblock)
, 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).
duDe(nblock,ndir+nshr)
Derivatives of strain energy potential with respect to the components of the modified Green strain tensor, 
.
duDj(nblock,ndir+nshr)
Derivatives of strain energy potential with respect to volume ratio, 
.
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), .
d2uDjDj(nblock)
Second derivatives of strain energy potential with respect to volume ratio, 
.
d2uDeDj(nblock,ndir+nshr)
Cross derivatives of strain energy potential with respect to components of the modified Green strain tensor and volume ratio, 
.
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).
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” 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, 
, 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.
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 ifVUANISOHYPER_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.
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, 
, as
is the fourth-order elasticity tensor. The derivatives of the strain energy function with respect to the Green strain are given as
However, subroutine VUANISOHYPER_STRAIN must return the derivatives of the strain energy function with respect to the modified Green strain tensor, , and the volume ratio, J, which can be accomplished easily using the following relationship between , , and 
:
is the second-order identity tensor. Thus, using the chain rule we find
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