黄永刚单晶塑性有限元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 (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 \ 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 \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 \ c variables\ 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 \C (2) changing \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 \ 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