黄永刚单晶塑性有限元umat子程序

SUBROUTINE UMAT(stress,statev,ddsdde,sse,spd,scd,

1 rpl, ddsddt, drplde, drpldt,

2 stran,dstran,time,dtime,temp,dtemp,predef,dpred,cmname,

3 ndi,nshr,ntens,nstatv,props,nprops,coords,drot,pnewdt,

4 celent,dfgrd0,dfgrd1,noel,npt,layer,kspt,kstep,kinc)

c WRITE (6,*) '

c NOTE: MODIFICATIONS TO *UMAT FOR ABAQUS VERSION 5.3 (14 APR '94)

c

c (3) As of version 5.3, ABAQUS files use double precision only.

c The file aba_param.inc has a line "implicit real*8" and, since

c it is included in the main subroutine, it will define the variables

c there as double precision. But other subroutines still need the

c definition "implicit real*8" since there may be variables that are

c not passed to them through the list or common block.

c

c (4) This is current as of version 5.6 of ABAQUS.

c

c (5) Note added by J. W. Kysar (4 November 1997). This UMAT has been

c with CFIX.

c

c The hardening law by Bassani and Wu was implemented incorrectly.

c This law is a function of both hyperbolic secant squared and hyperbolic

c tangent. However, the arguments of sech and tanh are related to the *total*

c slip on individual slip systems. Formerly, the UMAT implemented this

c hardening law by using the *current* slip on each slip system. Therein

c lay the problem. The UMAT did not restrict the current slip to be a

c positive value. So when a slip with a negative sign was encountered, the

c term containing tanh led to a negative hardening rate (since tanh is an

c of slip rates for each individual slip system. These "solution dependent

c variables" are available for postprocessing. The only required change

c in the input file is that the DEPVAR command must be changed.

c

C----- Use single precision on Cray by

C (1) deleting the statement "IMPLICIT*8 (A-H,O-Z)";

C (2) changing "REAL*8 FUNCTION" to "FUNCTION";

C (3) changing double precision functions DSIGN to SIGN.

C state

C

C SLIPSYS -- calculating number of slip systems, unit

C vectors in slip directions and unit normals to

C slip planes in a cubic crystal at the initial

C state

C

C GSLPINIT -- calculating initial value of current strengths

C at initial state

C

C ITERATION -- generating arrays for the Newton-Rhapson

C iteration

C

C LUDCMP -- LU decomposition

C

C LUBKSB -- linear equation solver based on LU

C decomposition method (must call LUDCMP first)

C----- Function subprogram:

C----- Variables:

C

C STRESS -- stresses (INPUT & OUTPUT)

C Cauchy stresses for finite deformation

C STATEV -- solution dependent state variables (INPUT & OUTPUT)

C DDSDDE -- Jacobian matrix (OUTPUT)

C----- Variables passed in for information:

C CMNAME -- name given in the *MATERIAL option

C NDI -- number of direct stress components

C NSHR -- number of engineering shear stress components

C NTENS -- NDI+NSHR

C NSTATV -- number of solution dependent state variables (as

C defined in the *DEPVAR option)

C PROPS -- material constants entered in the *USER MATERIAL

C option

C NPROPS -- number of material constants

C

C The strain increment is composed of elastic part and plastic

C part. The elastic strain increment corresponds to lattice

C stretching, the plastic part is the sum over all slip systems of

C plastic slip. The shear strain increment for each slip system is

C assumed a function of the ratio of corresponding resolved shear

C stress over current strength, and of the time step. The resolved

C shear stress is the double product of stress tensor with the slip

C deformation tensor (Schmid factor), and the increment of current

C strength is related to shear strain increments over all slip

C systems through self- and latent-hardening functions.

C variables within each increment.

C----- The present program is for a single CUBIC crystal. However,

C this code can be generalized for other crystals (e.g. HCP,

C Tetragonal, Orthotropic, etc.). Only subroutines ROTATION and

C SLIPSYS need to be modified to include the effect of crystal

C aspect ratio.

C

CFIX NSTATV >= 10 * NSLPTL + 5

C Denote s as a slip direction and m as normal to a slip plane.

C Here (s,-m), (-s,m) and (-s,-m) are NOT considered

C independent of (s,m). The number of slip systems in each set

C could be either 6, 12, 24 or 48 for a cubic crystal, e.g. 12

C for {110}<111>.

C

C Users who need more parameters to characterize the

C constitutive law of single crystal, e.g. the framework

C proposed by Zarka, should make NSTATV larger than (or equal

C in the input file, at the *USER MATERIAL card.

C

PARAMETER (ND=150)

C----- The parameter ND determines the dimensions of the arrays in

C this subroutine. The current choice 150 is a upper bound for a

C cubic crystal with up to three sets of slip systems activated.

C Users may reduce the parameter ND to any number as long as larger

C than or equal to the total number of slip systems in all sets.

C For example, if {110}<111> is the only set of slip system

CHARACTER*8 CMNAME

EXTERNAL F

dimension stress(ntens),statev(nstatv),

1 ddsdde(ntens,ntens),ddsddt(ntens),drplde(ntens),

2 stran(ntens),dstran(ntens),time(2),predef(1),dpred(1),

3 props(nprops),coords(3),drot(3,3),dfgrd0(3,3),dfgrd1(3,3)

7 DSTRES(6), DELATS(6), DSPIN(3), DVGRAD(3,3),

8 DGAMMA(ND), DTAUSP(ND), DGSLIP(ND),

9 WORKST(ND,ND), INDX(ND), TERM(3,3), TRM0(3,3), ITRM(3)

DIMENSION FSLIP1(ND), STRES1(6), GAMMA1(ND), TAUSP1(ND),

2 GSLP1(ND), SPNOR1(3,ND), SPDIR1(3,ND), DDSDE1(6,6),

3 DSOLD(6), DGAMOD(ND), DTAUOD(ND), DGSPOD(ND),

4 DSPNRO(3,ND), DSPDRO(3,ND),

5 DHDGDG(ND,ND)

C SLPDEF(1,i) -- SLPDIR(1,i)*SLPNOR(1,i)

C SLPDEF(2,i) -- SLPDIR(2,i)*SLPNOR(2,i)

C SLPDEF(3,i) -- SLPDIR(3,i)*SLPNOR(3,i)

C SLPDEF(4,i) -- SLPDIR(1,i)*SLPNOR(2,i)+

C SLPDIR(2,i)*SLPNOR(1,i)

C SLPDEF(5,i) -- SLPDIR(1,i)*SLPNOR(3,i)+

C SLPDIR(3,i)*SLPNOR(1,i)

C SLPDEF(6,i) -- SLPDIR(2,i)*SLPNOR(3,i)+

C SLPDIR(3,i)*SLPNOR(2,i)

C where index i corresponds to the ith slip system