KOKO\Mimi
375 lines
11 KiB
Markdown
375 lines
11 KiB
Markdown
<!-- 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
|
||
```
|