ABAQUS Analysis User

ABAQUS Analysis User's Manual

25.3.4

VUMAT

User subroutine to define material behavior.

Product:

ABAQUS/Explicit

Warning:

The

use

of

this

user

subroutine

generally

requires

considerable

expertise. You are cautioned that the implementation of any realistic

constitutive model requires extensive development and testing. Initial

testing on a single-element model with prescribed traction loading is

strongly recommended.

The

component

ordering

of

the

symmetric

and

nonsymmetric

tensors

for

the

three- dimensional case using C3D8R elements is different from the

ordering

specified

in

“

Three- dimensional

solid

element

library,

”

Section

14.1.4

, and the ordering used in ABAQUS/Standard.

References

“

User-defined mechanical material behavior,

”

Section 12.8.1

?

*USER MATERIAL

?

Overview

User subroutine VUMAT:

?

?

?

?

?

?

is used to define the mechanical constitutive behavior of a

material;

will

be

called

for

blocks

of

material

calculation

points

for

which

the material is defined in a user subroutine (

“

Material data

definition,

”

Section 9.1.2

);

can use and update solution-dependent state variables;

can use any field variables that are passed in;

is described further in

“

User- defined mechanical material

behavior,

”

Section 12.8.1

; and

cannot be used in an adiabatic analysis.

Component ordering in tensors

The component ordering depends upon whether the tensor is symmetric or

nonsymmetric.

Symmetric tensors

For symmetric tensors such as the stress and strain tensors, there are

ndir+nshr components, and the component order is given as a natural

permutation

of

the

indices

of

the

tensor.

The

direct

components

are

first

and then the indirect components, beginning with the 12-component. For

example,

a

stress

tensor

contains

ndir

direct

stress

components

and

nshr

shear stress components, which are passed in as

Component

2-D Case

3-D Case

1

2

3

4

5

6

The

shear

strain

components

in

user

subroutine

VUMAT

are

stored

as

tensor

components

and

not

as

engineering

components;

this

is

different

from

user

subroutine UMAT in ABAQUS/Standard, which uses engineering components.

Nonsymmetric tensors

For nonsymmetric tensors there are ndir+2*nshr components, and the

component order is given as a natural permutation of the indices of the

tensor.

The

direct

components

are

first

and

then

the

indirect

components,

beginning with the 12-component. For example, the deformation gradient

is passed as

Component

2-D Case

3-D Case

1

2

Component

2-D Case

3-D Case

3

4

5

6

7

8

9

Initial calculations and checks

In the data check phase of the analysis ABAQUS/Explicit calls user

subroutine VUMAT with a set of fictitious strains and a totalTime and

stepTime

both

equal

to

0.0.

This

is

done

as

a

check

on

your

constitutive

relation and to calculate the equivalent initial material properties,

based upon which the initial elastic wave speeds are computed.

Defining local orientations

All stresses, strains, stretches, and state variables are in the

orientation of the local material axes. These local material axes form

a basis system in which stress and strain components are stored. This

represents a corotational coordinate system in which the basis system

rotates with the material. If a user- specified coordinate system

(

“

Orientations,

”

Section 2.2.5

) is used, it defines the local material

axes in the undeformed configuration.

Special considerations for various element types

The use of user subroutine VUMAT requires special consideration for

various element types.

Shell and plane stress elements

You must define the stresses and internal state variables. In the case

of shell or plane stress elements, you must define strainInc(*,3), the

thickness strain increment. The internal energies can be defined if

desired. If they are not defined, the energy balance provided by

ABAQUS/Explicit will not be meaningful.

Shell elements

When VUMAT is used to define the material response of shell elements,

ABAQUS/Explicit

cannot

calculate

a

default

value

for

the

transverse

shear

stiffness

of

the

element.

Hence,

you

must

define

the

element's

transverse

shear stiffness. See

“

Shell section behavior,

”

Section 15.6.4

, for

guidelines on choosing this stiffness.

Beam elements

For

beam

elements

the

stretch

tensor

and

the

deformation

gradient

tensor

are

not

available.

For

beams

in

space

you

must

define

the

thickness

strains,

strainInc(*,2) and strainInc(*,3). strainInc(*,4) is the shear strain

associated with twist. Thickness stresses, stressNew(*,2) and

stressNew(*,3), are assumed to be zero and any values you assign are

ignored.

Deformation gradient

The polar decomposition of the deformation gradient is written as

, where

and

are the right and left symmetric

stretch tensors, respectively. The constitutive model is defined in a

corotational

coordinate

system

in

which

the

basis

system

rotates

with

the

material.

All

stress

and

strain

tensor

quantities

are

defined

with

respect

to

the

corotational

basis

system.

The

right

stretch

tensor,

,

is

used.

The relative spin tensor

represents the spin (the antisymmetric

part of the velocity gradient) defined with respect to the corotational

basis system.

Special considerations for hyperelasticity

Hyperelastic constitutive models in VUMAT should be defined in a

corotational

coordinate

system

in

which

the

basis

system

rotates

with

the

material. This is most effectively accomplished by formulating the

hyperelastic constitutive model in terms of the stretch tensor,

,

instead of in terms of the deformation gradient,

. Using the

deformation gradient can present some difficulties because the

deformation gradient includes the rotation tensor and the resulting

stresses would need to be rotated back to the corotational basis.

Objective stress rates

The

Green-Naghdi

stress

rate

is

used

when

the

mechanical

behavior

of

the

material

is

defined

using

user

subroutine

VUMAT.

The

stress

rate

obtained

with

user

subroutine

VUMAT

may

differ

from

that

obtained

with

a

built-in

ABAQUS

material

model.

For

example,

most

material

models

used

with

solid

(continuum)

elements

in

ABAQUS/Explicit

employ

the

Jaumann

stress

rate.

This

difference

in

the

formulation

will

cause

significant

differences

in

the results only if finite rotation of a material point is accompanied

by finite shear. For a discussion of the objective stress rates used in

ABAQUS, see

“

Stress rates,

”

Section 1.5.3 of the ABAQUS Theory Manual

.

Material point deletion

Material points that satisfy a user-defined failure criterion can be

deleted

from

the

model

(see

“

User-defined

mechanical

material

behavior,

”

Section

12.8.1

).

You

must

specify

the

state

variable

number

controlling

the element deletion flag when you allocate space for the

solution-dependent state variables, as explained in

“

User-defined

mechanical material behavior,

”

Section 12.8.1

. The deletion state

variable

should

be

set

to

a

value

of

one

or

zero

in

VUMAT.

A

value

of

one

indicates that the material point is active, while a value of zero

indicates

that

ABAQUS/Explicit

should

delete

the

material

point

from

the

model by setting the stresses to zero. The structure of the block of

material

points

passed

to

user

subroutine

VUMAT

remains

unchanged

during

the analysis; deleted material points are not removed from the block.

ABAQUS/Explicit will pass zero stresses and strain increments for all

deleted

material

points.

Once

a

material

point

has

been

flagged

as

deleted,

it cannot be reactivated.

User subroutine interface

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_'

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 25.1.1

, 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 25.1.1

).

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

upper-case 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 elements.

charLength(nblock)

Characteristic

element

length.

This

is

a

typical

length

of

a

line

across

an element. For beams and trusses, it 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

–

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

.

Stored

antisymmetric

part

of

the

velocity

gradient,

,

and

in 3-D as

tempOld(nblock)

and in 2-D as

.

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

.

defgradOld (nblock,ndir+2*nshr)

Deformation gradient tensor at each material point at the beginning of

the

increment.

Stored

in

3-D

as

(

,

,

,

,

,

,

,

,

) and in 2-D as (

,

,

,

,

).

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

3-D

as

(

,

,

,

,

,

,

,

,

)

and in 2-D as (

,

,

,

,

).

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 upper case since cmname is passed in as an

upper- case 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

The von Mises yield surface is defined as

where

is the uniaxial equivalent yield stress. The von Mises yield

surface is a cylinder in deviatoric stress space with a radius of

For

the

kinematic

hardening

model,

is

a

constant.

The

normal

to

the

Mises

yield surface can be written as

We decompose the strain rate into an elastic and plastic part using an

additive decomposition:

The plastic part of the strain rate is given by a normality condition

where the scalar multiplier

must be determined. A scalar measure of

equivalent plastic strain rate is defined by

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

where

and

are the Lamés constants for the material.

The evolution law for

is given as

where

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

The equivalent plastic strain rate is related to

by

The kinematic hardening constitutive model is integrated in a rate form

as follows. A trial elastic stress is computed as

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

where

From

the

definition

of

the

normal

to

the

yield

surface

at

the

end

of

the

increment,

,

This can be expanded using the incremental equations as

Taking the tensor product of this equation with

, using the the yield

condition at the end of the increment, and solving for

:

The

value

for

, and

is

used

in

the

incremental

equations

to

determine

.

,

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_'

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 -

smean = third * ( stressNew(i,1) + stressNew(i,2)

1 + stressNew(i,3) )

equivStress = sqrt( threeHalfs *

1 ( (stressNew(i,1)-smean)**2

1 + (stressNew(i,2)-smean)**2

1 + (stressNew(i,3)-smean)**2

1 + two * stressNew(i,4)**2 ) )