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

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 (1) The list of variables above defining the *UMAT subroutine,
c and the first (standard) block of variables dimensioned below,
c have variable names added compared to earlier ABAQUS versions.
c
c (2) The statement: include 'aba_param.inc' must be added as below.
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 modified to keep track of the cumulative shear strain in each
c individual slip system. This information is needed to correct an
c error in the implementation of the Bassani and Wu hardening law.
c Any line of code which has been added or modified is preceded
c immediately by a line beginning CFIXA and succeeded by a line
c beginning CFIXB. Any comment line added or modified will begin
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 odd function).

c The UMAT has been fixed by adding state variables to keep track of the
c *total* slip on each slip system by integrating up the absolute value
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
C----- Subroutines:
C
C ROTATION -- forming rotation matrix, i.e. the direction
C cosines of cubic crystal [100], [010] and [001]
C directions in global system at the initial
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 STRAINRATE -- based on current values of resolved shear
C stresses and current strength, calculating
C shear strain-rates in slip systems
C
C LATENTHARDEN -- forming self- and latent-hardening matrix
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 F -- shear strain-rates in slip systems

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
C STRAN -- strains
C logarithmic strain for finite deformation
C (actually, integral of the symmetric part of velocity
C gradient with respect to time)
C DSTRAN -- increments of strains
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----- This subroutine provides the plastic constitutive relation of
C single crystals for finite element code ABAQUS. The plastic slip
C of single crystal obeys the Schmid law. The program gives the
C choice of small deformation theory and theory of finite rotation
C and finite strain.
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----- The implicit integration method proposed by Peirce, Shih and
C Needleman (1984) is used here. The subroutine provides an option
C of iteration to solve stresses and solution dependent state
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

C----- Important notice:
C
C (1) The number of state variables NSTATV must be larger than (or
CFIX equal to) TEN (10) times the total number of slip systems in
C all sets, NSLPTL, plus FIVE (5)
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 to) the number of those parameters NPARMT plus nine times
C the total number of slip systems, NSLPTL, plus five
CFIX NSTATV >= NPARMT + 10 * NSLPTL + 5
C
C (2) The tangent stiffness matrix in general is not symmetric if
C latent hardening is considered. Users must declare "UNSYMM"
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
C potentially activated, ND could be taken as twelve (12).
c
include 'aba_param.inc'
c
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)

DIMENSION ISPDIR(3), ISPNOR(3), NSLIP(3),
2 SLPDIR(3,ND), SLPNOR(3,ND), SLPDEF(6,ND),
3 SLPSPN(3,ND), DSPDIR(3,ND), DSPNOR(3,ND),
4 DLOCAL(6,6), D(6,6), ROTD(6,6), ROTATE(3,3),
5 FSLIP(ND), DFDXSP(ND), DDEMSD(6,ND),
6 H(ND,ND), DDGDDE(ND,6),
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----- NSLIP -- number of slip systems in each set
C----- SLPDIR -- slip directions (unit vectors in the initial state)
C----- SLPNOR -- normals to slip planes (unit normals in the initial
C state)
C----- SLPDEF -- slip deformation tensors (Schmid factors)
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
C----- SLPSPN -- slip spin tensors (only needed for finite rotation)
C SLPSPN(1,i) -- [SLPDIR(1,i)*SLPNOR(2,i)-
C SLPDIR(2,i)*SLPNOR(1,i)]/2
C SLPSPN(2,i) -- [SLPDIR(3,i)*SLPNOR(1,i)-
C SLPDIR(1,i)*SLPNOR(3,i)]/2
C SLPSPN(3,i) -- [SLPDIR(2,i)*SLPNOR(3,i)-
C SLPDIR(3,i)*SLPNOR(2,i)]/2
C where index i corresponds to the ith slip system
C----- DSPDIR -- increments of slip directions
C----- DSPNOR -- increments of normals to slip planes
C
C----- DLOCAL -- elastic matrix in local cubic crystal system
C----- D -- elastic matrix in global system
C----- ROTD -- rotation matrix transforming DLOCAL to D
C
C----- ROTATE -- rotation matrix, direction cosines of [100], [010]
C and [001] of cubic crystal in global system
C
C----- FSLIP -- shear strain-rates in slip systems
C----- DFDXSP -- derivatives of FSLIP w.r.t x=TAUSLP/GSLIP, where
C TAUSLP is the resolved shear stress and GSLIP is the
C current strength
C
C----- DDEMSD -- double dot product of the elastic moduli tensor with
C the slip deformation tensor plus, only for finite
C rotation, the dot product of slip spin tensor with
C the stress
C
C----- H -- self- and latent-hardening matrix
C H(i,i) -- self hardening modulus of the ith slip
C system (no sum over i)
C H(i,j) -- latent hardening molulus of the ith slip
C system due to a slip in the jth slip system
C (i not equal j)
C
C----- DDGDDE -- derivatice of the shear strain increments in slip
C systems w.r.t. the increment of strains
C
C----- DSTRES -- Jaumann increments of stresses, i.e. corotational
C stress-increments formed on axes spinning with the
C material
C----- DELATS -- strain-increments associated with lattice stretching
C DELATS(1) - DELATS(3) -- normal strain increments
C DELATS(4) - DELATS(6) -- engineering shear strain
C increments
C----- DSPIN -- spin-increments associated with the material element
C DSPIN(1) -- component 12 of the spin tensor
C DSPIN(2) -- component 31 of the spin tensor
C DSPIN(3) -- component 23 of the spin tensor
C
C----- DVGRAD -- increments of deformation gradient in the current
C state, i.e. velocity gradient times the increment of
C time
C
C----- DGAMMA -- increment of shear strains in slip systems
C----- DTAUSP -- increment of resolved shear stresses in slip systems
C----- DGSLIP -- increment of current strengths in slip systems
C
C
C----- Arrays for iteration:
C
C FSLIP1, STRES1, GAMMA1, TAUSP1, GSLP1 , SPNOR1, SPDIR1,
C DDSDE1, DSOLD , DGAMOD, DTAUOD, DGSPOD, DSPNRO, DSPDRO,
C DHDGDG
C
C
C----- Solution dependent state variable STATEV:
C Denote the number of total slip systems by NSLPTL, which
C will be calculated in this code.
C
C Array STATEV:
C 1 - NSLPTL : current strength in slip systems
C NSLPTL+1 - 2*NSLPTL : shear strain in slip systems
C 2*NSLPTL+1 - 3*NSLPTL : resolved shear stress in slip systems
C
C 3*NSLPTL+1 - 6*NSLPTL : current components of normals to slip
C slip planes
C 6*NSLPTL+1 - 9*NSLPTL : current components of slip directions
C
CFIX 9*NSLPTL+1 - 10*NSLPTL : total cumulative shear strain on each
CFIX slip system (sum of the absolute
CFIX values of shear strains in each slip
CFIX system individually)
CFIX
CFIX 10*NSLPTL+1 : total cumulative shear strain on all
C slip systems (sum of the absolute
C values of shear strains in all slip
C systems)
C
CFIX 10*NSLPTL+2 - NSTATV-4 : additional parameters users may need
C to characterize the constitutive law
C of a single crystal (if there are
C any).
C
C