KOKO\Mimi
17 KiB
subroutine vumat(
C Read only (unmodifiable)variables -
1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
2 stepTime, totalTime, dt, cmname, coordMp, charLength,
3 props, density, strainInc, relSpinInc,
4 tempOld, stretchOld, defgradOld, fieldOld,
5 stressOld, stateOld, enerInternOld, enerInelasOld,
6 tempNew, stretchNew, defgradNew, fieldNew,
C Write only (modifiable)variables -
7 stressNew, stateNew, enerInternNew, enerInelasNew )
C
include 'vaba_param.inc'
C
dimension props(nprops), density(nblock), coordMp(nblock,*),
1 charLength(nblock), strainInc(nblock,ndir+nshr),
2 relSpinInc(nblock,nshr), tempOld(nblock),
3 stretchOld(nblock,ndir+nshr),
4 defgradOld(nblock,ndir+nshr+nshr),
5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev), enerInternOld(nblock),
7 enerInelasOld(nblock), tempNew(nblock),
8 stretchNew(nblock,ndir+nshr),
8 defgradNew(nblock,ndir+nshr+nshr),
9 fieldNew(nblock,nfieldv),
1 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
2 enerInternNew(nblock), enerInelasNew(nblock)
C
character*80 cmname
C
do 100 km = 1,nblock
user coding
100 continue
return
end
Variables to be defined
stressNew (nblock, ndir+nshr)
Stress tensor at each material point at the end of the increment.
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).
Variables that can be updated
enerInternNew (nblock)
Internal energy per unit mass at each material point at the end of the increment.
enerInelasNew (nblock)
Dissipated inelastic energy per unit mass at each material point at the end of the increment.
Variables passed in for information
nblock
Number of material points to be processed in this call to VUMAT.
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.
lanneal
Flag indicating whether the routine is being called during an annealing process. lanneal=0 indicates that the routine is being called during a normal mechanics increment. lanneal=1 indicates that this is an annealing process and you should re-initialize the internal state variables, stateNew, if necessary. Abaqus/Explicit will automatically set the stresses, stretches, and state to a value of zero during the annealing process.
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.
coordMp(nblock,*)
Material point coordinates. It is the midplane material point for shell elements and the centroid for beam and pipe elements.
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 of 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.
props(nprops)
User-supplied material properties.
density(nblock)
Current density at the material points in the midstep configuration. This value may be inaccurate in problems where the volumetric strain increment is very small. If an accurate value of the density is required in such cases, the analysis should be run in double precision. This value of the density is not affected by mass scaling.
strainInc (nblock, ndir+nshr)
Strain increment tensor at each material point.
relSpinInc (nblock, nshr)
Incremental relative rotation vector at each material point defined in the corotational system. Defined as , where is the antisymmetric part of the velocity gradient, , and \pmb { \Omega } = \dot { \mathbf { R } } \cdot \mathbf { R ^ { T } } . Stored in 3D as and in 2D as .
tempOld(nblock)
Temperatures at each material point at the beginning of the increment.
stretchOld (nblock, ndir+nshr)
Stretch tensor, , at each material point at the beginning of the increment defined from the polar decomposition of the deformation gradient by \mathbf { F } = \mathbf { R } \cdot \mathbf { U } .
defgradOld (nblock,ndir+2*nshr)
Deformation gradient tensor at each material point at the beginning of the increment. Stored in 3D as ( F _ { 1 1 } , F _ { 2 2 } , F _ { 3 3 } , F _ { 1 2 } , F _ { 2 3 } , F _ { 3 1 } , F _ { 2 1 } , F _ { 3 2 } , F _ { 1 3 } ) and in 2D as ( F _ { 1 1 } , F _ { 2 2 } , F _ { 3 3 } , F _ { 1 2 } , F _ { 2 1 } ) .
fieldOld (nblock, nfieldv)
Values of the user-defined field variables at each material point at the beginning of the increment.
stressOld (nblock, ndir+nshr)
Stress tensor at each material point at the beginning of the increment.
stateOld (nblock, nstatev)
State variables at each material point at the beginning of the increment.
enerInternOld (nblock)
Internal energy per unit mass at each material point at the beginning of the increment.
enerInelasOld (nblock)
Dissipated inelastic energy per unit mass at each material point at the beginning of the increment.
tempNew(nblock)
Temperatures at each material point at the end of the increment.
stretchNew (nblock, ndir+nshr)
Stretch tensor, , at each material point at the end of the increment defined from the polar decomposition of the deformation gradient by .
defgradNew (nblock,ndir+2*nshr)
Deformation gradient tensor at each material point at the end of the increment. Stored in 3D as ( F _ { 1 1 } F _ { 2 2 } , F _ { 3 3 } , F _ { 1 2 } , F _ { 2 3 } , F _ { 3 1 } , F _ { 2 1 } , F _ { 3 2 } , F _ { 1 3 } ) and in 2D as ( F _ { 1 1 } , F _ { 2 2 } , F _ { 3 3 } , F _ { 1 2 } , F _ { 2 1 } ) .
fieldNew (nblock, nfieldv)
Values of the user-defined field variables at each material point at the end of the increment.
Example: Using more than one user-defined material model
To use more than one user-defined material model, the variable cmname can be tested for different material names inside user subroutine VUMAT, as illustrated below:
if (cmname(1:4) .eq. 'MAT1') then
call VUMAT_MAT1(argument_list)
else if (cmname(1:4) .eq. 'MAT2') then
call VUMAT_MAT2 (argument_list)
end if
VUMAT_MAT1 and VUMAT_MAT2 are the actual user material subroutines containing the constitutive material models for each material MAT1 and MAT2, respectively. Subroutine VUMAT merely acts as a directory here. The argument list can be the same as that used in subroutine VUMAT. The material names must be in uppercase characters since cmname is passed in as an uppercase character string.
Example: Elastic/plastic material with kinematic hardening
As a simple example of the coding of subroutine VUMAT, consider the generalized plane strain case for an elastic/plastic material with kinematic hardening. The basic assumptions and definitions of the model are as follows.
Let be the current value of the stress, and define to be the deviatoric part of the stress. The center of the yield surface in deviatoric stress space is given by the tensor , which has initial values of zero. The stress difference, , is the stress measured from the center of the yield surface and is given by
\boldsymbol {\xi} = \mathbf {S} - \boldsymbol {\alpha}.
The von Mises yield surface is defined as
f (\pmb {\sigma}) = \frac {1}{2} \pmb {\xi}: \pmb {\xi} - \frac {1}{3} \sigma_ {0} ^ {2},
where \sigma _ { 0 } is the uniaxial equivalent yield stress. The von Mises yield surface is a cylinder in deviatoric stress space with a radius of
R = \sqrt {\frac {2}{3}} \sigma_ {0}.
For the kinematic hardening model, R is a constant. The normal to the Mises yield surface can be written as
\mathbf {Q} = \sqrt {\frac {3}{2}} \frac {\boldsymbol {\xi}}{\sigma_ {0}}.
We decompose the strain rate into an elastic and plastic part using an additive decomposition:
\dot {\epsilon} = \dot {\epsilon} ^ {e l} + \dot {\epsilon} ^ {p l}.
The plastic part of the strain rate is given by a normality condition
\dot {\epsilon} ^ {p l} = \dot {\gamma} \mathbf {Q},
where the scalar multiplier \dot { \gamma } must be determined. A scalar measure of equivalent plastic strain rate is defined by
\dot {\bar {\epsilon}} ^ {p l} = \sqrt {\frac {2}{3} \dot {\epsilon} ^ {p l} : \dot {\epsilon} ^ {p l}}.
The stress rate is assumed to be purely due to the elastic part of the strain rate and is expressed in terms of Hooke’s law by
\dot {\pmb {\sigma}} = \lambda \mathrm{trace} (\dot {\pmb {\epsilon}} ^ {e l}) \mathbf {I} + 2 \mu \dot {\pmb {\epsilon}} ^ {e l},
where and 2 \mu are the Lamés constants for the material.
The evolution law for is given as
\dot {\alpha} = \frac {2}{3} \dot {\gamma} H \mathbf {Q},
where H is the slope of the uniaxial yield stress versus plastic strain curve.
During active plastic loading the stress must remain on the yield surface, so that
\sqrt {\mathbf {Q} : \mathbf {Q}} = 1.
The equivalent plastic strain rate is related to \dot { \gamma } by
\dot {\bar {\epsilon}} ^ {p l} = \sqrt {\frac {2}{3}} \dot {\gamma}.
The kinematic hardening constitutive model is integrated in a rate form as follows. A trial elastic stress is computed as
\pmb {\sigma} _ {n e w} ^ {t r i a l} = \pmb {\sigma} _ {o l d} + \lambda \mathrm{trace} (\Delta \pmb {\epsilon}) \mathbf {I} + 2 \mu \Delta \pmb {\epsilon},
where the subscripts and refer to the beginning and end of the increment, respectively. If the trial stress does not exceed the yield stress, the new stress is set equal to the trial stress. If the yield stress is exceeded, plasticity occurs in the increment. We then write the incremental analogs of the rate equations as
\pmb {\sigma} _ {n e w} = \pmb {\sigma} _ {n e w} ^ {t r i a l} - 2 \mu \pmb {\Delta} \pmb {\epsilon} ^ {p l} = \pmb {\sigma} _ {n e w} ^ {t r i a l} - 2 \mu \Delta \gamma \mathbf {Q},
\boldsymbol {\alpha} _ {n e w} = \boldsymbol {\alpha} _ {o l d} + \frac {2}{3} H \Delta \gamma \mathbf {Q},
\bar {\epsilon} _ {n e w} ^ {p l} = \bar {\epsilon} _ {o l d} ^ {p l} + \sqrt {\frac {2}{3}} \Delta \gamma ,
where
\Delta \gamma = \dot {\gamma} \Delta t.
From the definition of the normal to the yield surface at the end of the increment, ,
\alpha_ {n e w} + \sqrt {\frac {2}{3}} \sigma_ {0} \mathbf {Q} = \mathbf {S} _ {n e w}.
This can be expanded using the incremental equations as
\pmb {\alpha} _ {o l d} + \frac {2}{3} H \Delta \gamma \mathbf {Q} + \sqrt {\frac {2}{3}} \sigma_ {0} \mathbf {Q} = \mathbf {S} _ {n e w} ^ {t r i a l} - \Delta \gamma 2 \mu \mathbf {Q}.
Taking the tensor product of this equation with , using the yield condition at the end of the increment, and solving for \Delta \gamma \mathrm { : } :
\Delta \gamma = \frac {1}{2 \mu (1 + H / 3 \mu)} \left(\left(\pmb {\xi} _ {n e w} ^ {t r i a l}: \pmb {\xi} _ {n e w} ^ {t r i a l}\right) ^ {1 / 2} - \sqrt {\frac {2}{3}} \sigma_ {0}\right).
The value for \Delta \gamma is used in the incremental equations to determine \sigma _ { n e w } , \alpha _ { n e w } , , and \overline { { \epsilon } } _ { n e w } ^ { p l }
This algorithm is often referred to as an elastic predictor, radial return algorithm because the correction to the trial stress under the active plastic loading condition returns the stress state to the yield surface along the direction defined by the vector from the center of the yield surface to the elastic trial stress. The subroutine would be coded as follows:
subroutine vumat(
C Read only -
1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
2 stepTime, totalTime, dt, cmname, coordMp, charLength,
3 props, density, strainInc, relSpinInc,
4 tempOld, stretchOld, defgradOld, fieldOld,
3 stressOld, stateOld, enerInternOld, enerInelasOld,
6 tempNew, stretchNew, defgradNew, fieldNew,
C Write only -
5 stressNew, stateNew, enerInternNew, enerInelasNew )
C
include 'vaba_param.inc'
C
C J2 Mises Plasticity with kinematic hardening for plane
C strain case.
C Elastic predictor, radial corrector algorithm.
C
C The state variables are stored as:
C STATE(*,1) = back stress component 11
C STATE(*,2) = back stress component 22
C STATE(*,3) = back stress component 33
C STATE(*,4) = back stress component 12
C STATE(*,5) = equivalent plastic strain
C
C
C All arrays dimensioned by (*) are not used in this algorithm
dimension props(nprops), density(nblock),
1 coordMp(nblock,*),
2 charLength(*), strainInc(nblock,ndir+nshr),
3 relSpinInc(*), tempOld(*),
4 stretchOld(*), defgradOld(*),
5 fieldOld(*), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev), enerInternOld(nblock),
7 enerInelasOld(nblock), tempNew(*),
8 stretchNew(*), defgradNew(*), fieldNew(*),
9 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
1 enerInternNew(nblock), enerInelasNew(nblock)
C
character*80 cmname
C
parameter( zero = 0., one = 1., two = 2., three = 3.,
1 third = one/three, half = .5, twoThirds = two/three,
2 threeHalfs = 1.5 )
C
e = props(1)
xnu = props(2)
yield = props(3)
hard = props(4)
C
twomu = e / ( one + xnu )
thremu = threeHalfs * twomu
sixmu = three * twomu
alamda = twomu * ( e - twomu ) / ( sixmu - two * e )
term = one / ( twomu * ( one + hard/thremu ) )
con1 = sqrt( twoThirds )
C
do 100 i = 1,nblock
C
C Trial stress
trace = strainInc(i,1) + strainInc(i,2) + strainInc(i,3)
sig1 = stressOld(i,1) + alamda*trace + twomu*strainInc(i,1)
sig2 = stressOld(i,2) + alamda*trace + twomu*strainInc(i,2)
sig3 = stressOld(i,3) + alamda*trace + twomu*strainInc(i,3)
sig4 = stressOld(i,4) + twomu*strainInc(i,4)
C
C Trial stress measured from the back stress
s1 = sig1 - stateOld(i,1)
s2 = sig2 - stateOld(i,2)
s3 = sig3 - stateOld(i,3)
s4 = sig4 - stateOld(i,4)
C
C Deviatoric part of trial stress measured from the back stress
smean = third * (s1 + s2 + s3)
ds1 = s1 - smean
ds2 = s2 - smean
ds3 = s3 - smean
C
C Magnitude of the deviatoric trial stress difference
dsmag = sqrt(ds1**2 + ds2**2 + ds3**2 + 2.*s4**2)
C
C Check for yield by determining the factor for plasticity,
C zero for elastic, one for yield
radius = con1 * yield
facyld = zero
if( dsmag - radius .ge. zero ) facyld = one
C
C Add a protective addition factor to prevent a divide by zero
C when dsmag is zero. If dsmag is zero, we will not have exceeded
C the yield stress and facyld will be zero.
dsmag = dsmag + (one - facyld)
C
C Calculated increment in gamma (this explicitly includes the
C time step)
diff = dsmag - radius
dgamma = facyld * term * diff
C
C Update equivalent plastic strain
deqps = con1 * dgamma
stateNew(i,5) = stateOld(i,5) + deqps
C
C Divide dgamma by dsmag so that the deviatoric stresses are
C explicitly converted to tensors of unit magnitude in the
C following calculations
dgamma = dgamma / dsmag
C
C Update back stress
factor = hard * dgamma * twoThirds
stateNew(i,1) = stateOld(i,1) + factor * ds1
stateNew(i,2) = stateOld(i,2) + factor * ds2
stateNew(i,3) = stateOld(i,3) + factor * ds3
stateNew(i,4) = stateOld(i,4) + factor * s4
C
C Update the stress
factor = twomu * dgamma
stressNew(i,1) = sig1 - factor * ds1
stressNew(i,2) = sig2 - factor * ds2
stressNew(i,3) = sig3 - factor * ds3
stressNew(i,4) = sig4 - factor * s4
C
C Update the specific internal energy -
stressPower = half * (
1 ( stressOld(i,1)+stressNew(i,1) ) *strainInc(i,1)
1 + ( stressOld(i,2)+stressNew(i,2) ) *strainInc(i,2)
1 + ( stressOld(i,3)+stressNew(i,3) ) *strainInc(i,3)
1 + two * ( stressOld(i,4)+stressNew(i,4) ) *strainInc(i,4) )
C
enerInternNew(i) = enerInternOld(i)
1 + stressPower / density(i)
C
C Update the dissipated inelastic specific energy -
plasticWorkInc = dgamma * half * (
1 ( stressOld(i,1)+stressNew(i,1) ) *ds1
1 + ( stressOld(i,2)+stressNew(i,2) ) *ds2
1 + ( stressOld(i,3)+stressNew(i,3) ) *ds3
1 + two * ( stressOld(i,4)+stressNew(i,4) ) *s4 )
enerInelasNew(i) = enerInelasOld(i)
1 + plasticWorkInc / density(i)
100 continue
C
return
end