<!-- source-page: 491 -->

```python
parameter (zero = 0.d0, half = 0.5d0, one = 1.d0, two=2.d0)
c
    dimension rhs(nblock,ndofel), amass(nblock,ndofel,ndofel),
    *    dtimeStable(nblock),
    *    svars(nblock,nsvars), energy(nblock,nElEnergy),
    *    props(nprops), jprops(njprops),
    *    jElem(nblock), time(nTime), l(nFlags),
    *    coords(nblock,nnode,ncrd), u(nblock,ndofel),
    *    du(nblock,ndofel), v(nblock,ndofel), a(nblock, ndofel),
    *    predef(nblock, nnode, npred, nPred), adlmag(nblock),
    *    dMassScaleFactor(nblock)

c    Notes:
c    Define only nonzero entries; the arrays to be defined
c    have been zeroed out just before this call

if (jtype .eq. 1001 .and.
    *    lflags(iProcedure).eq.jDynExplicit) then

    area0 = props(1)
    eMod = props(2)
    anu = props(3)
    rho = props(4)

    eDampTra = zero
    amassFact0 = half*area0*rho

    if ( lflags(iOpCode).eq.jMassCalc ) then
    do kblock = 1, nblock

c    use original distance to compute mass
    alenX0 = (coords(kblock,2,1) - coords(kblock,1,1))
    alenY0 = (coords(kblock,2,2) - coords(kblock,1,2))
    alenZ0 = (coords(kblock,2,3) - coords(kblock,1,3))
    alen0 = sqrt(alenX0*alenX0 + alenY0*alenY0 +
    *    alenZ0*alenZ0)
    am0 = amassFact0*alen0
    amass(kblock,1,1) = am0
    amass(kblock,2,2) = am0
    amass(kblock,3,3) = am0
    amass(kblock,4,4) = am0
    amass(kblock,5,5) = am0 
```

<!-- source-page: 492 -->

```txt
amass(kblock,6,6) = am0
end do
else if ( lflags(iOpCode) .eq.
*    jIntForceAndDtStable) then
do kblock = 1, nblock
    alenX0 = (coords(kblock,2,1) - coords(kblock,1,1))
    alenY0 = (coords(kblock,2,2) - coords(kblock,1,2))
    alenZ0 = (coords(kblock,2,3) - coords(kblock,1,3))
    alen0 = sqrt(alenX0*alenX0 + alenY0*alenY0 +
*    alenZ0*alenZ0)
    vol0 = area0*alen0
    amElem0 = two*amassFact0*alen0
    alenX = alenX0
*    + (u(kblock,4) - u(kblock,1))
    alenY = alenY0
*    + (u(kblock,5) - u(kblock,2))
    alenZ = alenZ0
*    + (u(kblock,6) - u(kblock,3))
    alen = sqrt(alenX*alenX + alenY*alenY + alenZ*alenZ)
    area = vol0/alen
    ak = area*eMod/alen
c    stable time increment for translations
    dtimeStable(kblock) = factorStable*sqrt(amElem0/ak)
c    force = E * logarithmic strain *current area
    strainLog = log(alen/alen0)
    fElasTra = eMod*strainLog*area
    forceTra = fElasTra
c    assemble internal load in RHS
    rhs(kblock,1) = -forceTra
    rhs(kblock,4) = forceTra
c    internal energy calculation
    alenOld = svars(kblock,1)
    fElasTraOld = svars(kblock,2)
    energy(kblock, iElIe) = energy(kblock, iElIe) + 
```

<!-- source-page: 493 -->

```fortran
* half*(fElasTra+fElasTraOld)*(alen - alenOld)
c update state variables
    svars(kblock,1) = alen
    svars(kblock,2) = fElasTra
    end do
    end if
    end if
c
    return
    end 
```

<!-- source-page: 494 -->

<!-- source-page: 495 -->

# 1.2.15 VUEOS: User subroutine to define equation of state material model.

# Product: Abaqus/Explicit

# References

• “Equation of state,” Section 25.2.1 of the Abaqus Analysis User’s Guide   
• \*EOS   
“Equation of state material,” Section 2.2.20 of the Abaqus Verification Guide

# Overview

User subroutine VUEOS:

• can be used to define the hydrodynamic material model in which the material’s volumetric response is determined by the user-defined equation of state;   
• will be called for blocks of material calculation points for which the material definition contains a user-defined equation of state;   
• can use and update solution-dependent state variables; and   
• can use any field variables that are passed in.

# User subroutine interface

```python
subroutine vueos (
C Read only (unmodifiable) variables -
    1    nblock,
    2    jElem, kIntPt, kLayer, kSecPt,
    3    steTime, totalTime, dt, cmname,
    4    nstatev, nfieldv, nprops,
    5    props, tempOld, tempNew, fieldOld, fieldNew,
    6    stateOld, charLength, coordMp,
    7    densityMean, refDensity, densityNew,
    8    dkk, Em,
C Write only (modifiable) variables -
    8    press, dPdRho, dPdEm,
    9    stateNew )
C
    include 'vaba_param.inc'
C
    dimension props(nprops),
    1    tempOld(nblock),
    2    fieldOld(nblock, nfieldv), 
```

<!-- source-page: 496 -->

```matlab
3 stateOld(nblock, nstatev),
4 tempNew(nblock),
5 fieldNew(nblock, nfieldv),
6 charLength(nblock), coordMp(nblock, *),
7 densityMean(nblock), refDensity(nblock),
8 densityNew(nblock),
9 dkk(nblock), Em(nblock),
1 press(nblock), dPdRho(nblock), dPdEm(nblock),
2 stateNew(nblock)
C
    character*80 cmname
C
    do 100 km = 1, nblock
    user coding to define/update press, dPdRho, dPdEm
100 continue
    return
    end 
```

# Variables to be defined

press(nblock)

The material point pressure stress, p.

dPdRho(nblock)

The derivative of the pressure with respect to the density, $\partial p / \partial \rho .$ This quantity is needed for the evaluation of the effective moduli of the material, which enters the stable time increment calculation.

dPdEm(nblock)

The derivative of the pressure with respect to the internal energy, $\partial p / \partial E _ { m }$ . This quantity is needed for the iterative Newton loop used outside of the user subroutine to solve for pressure.

# Variable that can be updated

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).

<!-- source-page: 497 -->

# Variables passed in for information

nblock

Number of material points to be processed in this call to VUEOS.

jElem(nblock)

Array of element numbers.

kIntPt

Integration point number.

kLayer

Layer number.

kSecPt

Section point number within the current layer.

stepTime

Value of time since the step began.

totalTime

Value of total time. The time at the beginning of the step is given by totalTime - stepTime.

dt

Time increment size.

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 the number 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.

<!-- source-page: 498 -->

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.

charLength(nblock)

Characteristic element length, which is either the default value based on the geometric mean or the user-defined characteristic element length defined in user subroutine VUCHARLENGTH. The default value is a typical length of a line across an element for a first-order element; it is half the same typical length for a second-order element. For beams, pipes, and trusses, the default value is a characteristic length along the element axis. For membranes and shells it is a characteristic length in the reference surface. For axisymmetric elements it is a characteristic length in the r–z plane only. For cohesive elements it is equal to the constitutive thickness.

coordMp(nblock,\*)

Material point coordinates.

densityMean(nblock)

The mean density.

refDensity(nblock)

The reference density.

densityNew(nblock)

The current density for this increment.

dkk(nblock)

The volumetric strain increment.

Em(nblock)

The element specific internal energy (per unit mass)

Example: User subroutine VUEOS to reproduce results obtained with \*EOS, TYPE=USUP

As a simple example of coding of user subroutine VUEOS, consider the following form of the Mie-Grüneisen equation of state with 0.0 and a linear dependency between pressure and internal energy:

<!-- source-page: 499 -->

$$
p = \rho_ {0} c _ {0} ^ {2} \eta (1 - \frac {\Gamma_ {0} \eta}{2}) + \Gamma_ {0} \rho_ {0} E _ {m},
$$

where $\eta = 1 - \rho _ { 0 } / \rho .$ . Therefore, the results obtained with user subroutine VUEOS should be the same as the results obtained with the linear $U _ { s } - U _ { p }$ type of EOS already available in “Equation of state material,” Section 2.2.20 of the Abaqus Verification Guide.

The code in user subroutine VUEOS must return the pressure, , as in the above equation; the derivative of the pressure with respect to the density, $\partial p / \partial \rho ;$ and the derivative of the pressure with respect to the energy, $\partial p / \partial E _ { m }$ . For the case considered here, these values are

$$
\frac {\partial p}{\partial \rho} = \frac {\rho_ {0} ^ {2} c _ {0} ^ {2}}{\rho^ {2}} (1 - \Gamma_ {0} \eta),
$$

$$
\frac {\partial p}{\partial E _ {m}} = \Gamma_ {0} \rho_ {0}.
$$```fortran
subroutine vueos (
C Read only (unmodifiable) variables -
    1    nblock,
    2    jElem, kIntPt, kLayer, kSecPt,
    3    steTime, totalTime, dt, cmname,
    4    nstatev, nfieldv, nprops,
    5    props, tempOld, tempNew, fieldOld, fieldNew,
    6    stateOld, charLength, coordMp,
    7    densityMean, refDensity, densityNew,
    8    dkk, Em,
C Write only (modifiable) variables -
    8    press, dPdRho, dPdEm,
    9    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    charLength(nblock), coordMp(nblock,*),
    7    densityMean(nblock), refDensity(nblock),
    8    densityNew(nblock),
    9    dkk(nblock), Em(nblock),
    1    press(nblock), dPdRho(nblock), dPdEm(nblock), 
```

<!-- source-page: 500 -->

```txt
2 stateNew(nblock)
C
    character*80 cmname
C
    parameter ( zero = 0.d0, one = 1.d0, half = 0.5d0 )
C
    c0 = props(1)
    gamma0 = props(2)
    c02 = c0*c0
C
    do k=1, nblock
    rho0 = refDensity(k)
    eta = one - rho0/densityNew(k)
    f1 = rho0*c02*eta*(one-half*gamma0*eta)
    f2 = gamma0*rho0
    press(k) = f1 + f2*Em(k)
C dP/dEm
    dPdEm(k) = f2
C dP/dRho
    dPdRho(k) = c02*(rho0/densityNew(k))**2*(one-gamma0*eta)
    end do
C
    return
    end 
```