VUEOS

User subroutine to define equation of state material model.

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.

The following topics are discussed:

User subroutine interface
Variables to be defined
Variables that can be updated
Variables passed in for information
Example: User subroutine VUEOS to reproduce results obtained with EOS, TYPE=USUP

User subroutine interface

      subroutine vueos (
C Read only (unmodifiable) variables –
     1     nblock,
     2     jElem, kIntPt, kLayer, kSecPt, 
     3     stepTime, 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),
     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 ρ$ . 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.

Variables 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 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 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).
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.
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 $s = $ 0.0 and a linear dependency between pressure and internal energy:

p = ρ_{0} c_{0}^{2} η (1 - \frac{Γ_{0} η}{2}) + Γ_{0} ρ_{0} E_{m},

where $η = 1 - ρ_{0} / ρ$ . 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.

The code in user subroutine VUEOS must return the pressure, $p$ , as in the above equation; the derivative of the pressure with respect to the density, $\partial p / \partial ρ$ ; 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 ρ} = \frac{ρ_{0}^{2} c_{0}^{2}}{ρ^{2}} (1 - Γ_{0} η),

\frac{\partial p}{\partial E_{m}} = Γ_{0} ρ_{0} .

      subroutine vueos (
C Read only (unmodifiable) variables –
     1     nblock,
     2     jElem, kIntPt, kLayer, kSecPt, 
     3     stepTime, 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),
     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