11.1.3 Hyperelasticity USRRUB

If RUBBER and COMPRE are specified with USER in table 'MATERI' [§4.4.3], then subroutine USRRUB must be supplied to define the deviatoric and hydrostatic strain energy functions. See also §16.3 for background theory.

    (Fortran)

      SUBROUTINE USRRUB( eps0, deps, ns, c, iinvar, jinvar,
                         stretc, didc, djdc, age0, dtime, temp0,
                         dtemp, elemen, intpt, coord, iter,
                         rubval, nrv, comval, ncv, usrsta, nus,
                         usrind, nui, sig, wp, stiff, dwpdj,
                         press, compre, usrkey, d2di1,
                         d2di2, d2di3, d2dj1, d2dj2, d2dj3 )

c

is the column representation of the Right Cauchy-Green stretch tensor C , respectively C₁₁ , C₂₂ , C₃₃ , C₁₂ , C₂₃ , C₃₁ for three-dimensional and C₁₁ , C₂₂ , C₃₃ , C₁₂ for two-dimensional.

iinvar

are the regular invariants of C , respectively I₁ , I₂ , I₃ .

jinvar

are the modified invariants of C , respectively J₁ , J₂ , J .

stretc

are the principal stretches of C , respectively $\lambda_{{1}}^{}$ , $\lambda_{{2}}^{}$ , $\lambda_{{3}}^{}$ .

rubval

are the parameters for the deviatoric strain energy function as specified via RUBVAL rubval in input table 'MATERI'. These variables shall not be updated in this subroutine!

comval

are the parameters for the hydrostatic strain energy function as specified via COMVAL comval in input table 'MATERI'. These variables shall not be updated in this subroutine!

usrsta

are the user-supplied state variables.

usrind

are the user-supplied integer indicators.

At the start of the first step the values of usrsta and usrind come in as the ones specified in input table 'MATERI', respectively via USRSTA usrsta and USRIND usrind. In subsequent steps these values come in from the final iteration of the previous step. They must always go out for the current iteration.

sig

is the 2nd Piola-Kirchhoff stress vector $\boldsymbol\tau$ = $\boldsymbol\tau$_d + $\boldsymbol\tau$_h . The values come in from the final iteration of the previous step and must go out for the current iteration.

stiff

is the tangent stiffness matrix, or stress-strain relation D = D_d + D_h . The values come in from the previous step and must go out for the current iteration.

wp

is the pressure function w_p .

dwpdj

is the derivative of the pressure function $\partial$w_p/$\partial$J .

usrkey

is the user-supplied keyword from input table 'MATERI'. This can be used as a switch, to code various models within one subroutine.

Example.

As an example we will show how the Neo-Hookean material model for rubber can be coded via user-supplied subroutine USRRUB. The strain energy function reads:

W = W_d + W_h (11.1)

with W_d the deviatoric part and W_h the hydrostatic part. The deviatoric part of the strain energy function reads

W_d = K₁(J₁ - 1) (11.2)

Then the deviatoric part of the stress increment can be calculated as

$$\boldsymbol\tau {\frac{{ \partial W_{\mathrm{d}} }}{{ \partial C }}} {\frac{{ \partial J_{1} }}{{ \partial C }}}$$ (11.3)

and the deviatoric part of the stress-strain relation as

$${\frac{{ \partial^{2} W_{\mathrm{d}} }}{{ \partial C^{2} }}} {\frac{{ \partial^{2} J_{1} }}{{ \partial C^{2} }}}$$ (11.4)

The hydrostatic part of the strain energy function W_h for linear compressibility reads

$${\frac{{ \kappa }}{{ 2 }}}$$ (11.5)

then the hydrostatic part of the stress increment can be calculated as

$$\boldsymbol\tau {\frac{{ \partial W_{\mathrm{h}} }}{{ \partial C }}} \kappa \left(\vphantom{ J-1 }\right. \left.\vphantom{ J-1 }\right) {\frac{{ \partial J_{1} }}{{ \partial C }}} {\frac{{ \partial J_{1} }}{{ \partial C }}}$$ (11.6)

with bulk modulus $\kappa$ and $\kappa$$\left(\vphantom{ J-1 }\right.$J - 1$\left.\vphantom{ J-1 }\right)$ as the pressure p . The hydrostatic part of the stress-strain relation can be calculated as

$${\frac{{ \partial^{2} W_{\mathrm{h}} }}{{ \partial C^{2} }}} {\frac{{ \partial^{2} J_{1} }}{{ \partial C^{2} }}}$$ (11.7)

For the pressure at time t

$${\frac{{ \partial W_{\mathrm{h}} }}{{ \partial J }}} \kappa \left(\vphantom{ J - 1 }\right. \left.\vphantom{ J - 1 }\right) \kappa$$ (11.8)

a pressure function w_p(J) can be calculated as

w_p = (J - 1) (11.9)

and the derivative of the pressure function with respect to J for the pressure increment

$${\frac{{ \partial w_{\mathrm{p}} }}{{ \partial J }}}$$ (11.10)

The implementation of the Neo-Hookean model into the user-supplied subroutine for rubber is as follows.

    (file.f)

      SUBROUTINE USRRUB( EPS0, DEPS, NS, C, IINVAR, JINVAR, STRETC,
     $                   DIDC, DJDC, AGE0, DTIME, TEMP0, DTEMP, ELEMEN,
     $                   INTPT, COORD, ITER, RUBVAL, NRV, COMVAL, NCV,
     $                   USRSTA, NUS, USRIND, NUI, SIG, STIFF, WP,
     $                   DWPDJ, PRESS, COMPRE, MODNAM, DER2I1,
     $                   DER2I2, DER2I3, DER2J1, DER2J2, DER2J3 )
C
      INTEGER          NS, ELEMEN, INTPT, ITER, NRV, NCV, NUS, NUI,
     $                 USRIND(*)
      DOUBLE PRECISION EPS0(*), DEPS(*), C(*), IINVAR(*),
     $                 JINVAR(*), STRETC(*), DIDC(NS,*), DJDC(NS,*),
     $                 AGE0, DTIME, TEMP0, DTEMP, COORD(*), USRSTA(*),
     $                 STIFF(NS,*), SIG(*), WP, DWPDJ, RUBVAL(*),
     $                 COMVAL(*), PRESS, DER2I1(NS,*),
     $                 DER2I2(NS,*), DER2I3(NS,*), DER2J1(NS,*),
     $                 DER2J2(NS,*), DER2J3(NS,*)
      CHARACTER        MODNAM*6
      LOGICAL          COMPRE
C
C...     USER VARIABLES
      DOUBLE PRECISION K1, R1, R3, CONS
      INTEGER          NS2
C
C...     NEO-HOOKEAN MODEL 2D & 3D
C...     USER SUPPLIED SUBROUTINE FOR RUBBER
C...     RETURN UPDATED STRESS AND TANGENTIAL STIFFNESS
C
      COMPRE = .TRUE.
C
      NS2 = NS*NS
C
      K1 = RUBVAL(1)
      R1 = 4.D0 * K1
      R3 = 4.D0 * PRESS
C
C...     CALCULATE DEVIATORIC PART OF UPDATED STRESSES
      CONS = 2.D0 * K1
      CALL UVS( DJDC(1,1), NS, CONS, SIG )
C
C...     CALCULATE HYDROSTATIC PART OF UPDATED STRESSES
      CONS  = 2.D0
      CONS = CONS * PRESS
      CALL UVPWS( SIG, DJDC(1,3), NS, CONS, SIG )
C
C...     CALCULATE DEVIATORIC PART OF TANGENT STIFFNESS FOR
C...     STRESS INCREMENT
      CALL UVS(   DER2J1, NS2, R1, STIFF )
C
C...     CALCULATE HYDROSTATIC PART OF TANGENT STIFFNESS FOR
C...     STRESS INCREMENT
      CALL UVPWS( STIFF, DER2J3, NS2, R3, STIFF )
C
C...     CALCULATE PRESSURE FUNCTION FOR PRESSURE
      WP = JINVAR(3) - 1
C
C...     CALCULATE DERIVATIVE OF PRESSURE FUNCTION FOR
C...     PRESSURE INCREMENT
      DWPDJ = 1.D0
C
      CALL FILMA( 1, STIFF, NS )
C
      END

The following could be the input data file for this example.

    (file.dat)

CHX64 / TENSION TEST / WD: NEO-HOOKEAN   WH: LINEAR   WC: -
USED UNITS:  N, MM
'COORDINATES'
    1     0.0     0.0     0.0
    2    50.0     0.0     0.0
    3   100.0     0.0     0.0
    4   100.0    50.0     0.0
    5   100.0   100.0     0.0
    6    50.0   100.0     0.0
    7     0.0   100.0     0.0
    8     0.0    50.0     0.0
    9     0.0     0.0    50.0
   10   100.0     0.0    50.0
   11   100.0   100.0    50.0
   12     0.0   100.0    50.0
   13     0.0     0.0   100.0
   14    50.0     0.0   100.0
   15   100.0     0.0   100.0
   16   100.0    50.0   100.0
   17   100.0   100.0   100.0
   18    50.0   100.0   100.0
   19     0.0   100.0   100.0
   20     0.0    50.0   100.0
'ELEMENTS'
CONNEC
   1  CHX64  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
MATERIALS
 / 1 / 1
DATA
 / 1 / 1
'DATA'
   1  NINTEG  3      3      3
      NUMINT  GAUSS  GAUSS  GAUSS
'MATERIALS'
   1  RUBBER  USER
:             K1
      RUBVAL  0.10
      COMPRE  USER
      BULK    1000.
'DIRECTIONS'
   1  1. 0. 0.
   2  0. 1. 0.
   3  0. 0. 1.
'SUPPORTS'
/ 1 7 13 19 / TR 1
/ 1 3 13 15 / TR 2
/ 1-8 /   TR 3
/ 13-20 / TR 3
'LOADS'
CASE 1
DEFORM
/ 13-20 / TR 3 / 0.1(8) /
'END'