资源描述
_ 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. cc (2) The statement: include aba_param.inc must be added as below.cc (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 variablesc there as double precision. But other subroutines still need thec definition implicit real*8 since there may be variables that arec not passed to them through the list or common block.cc (4) This is current as of version 5.6 of ABAQUS.cc (5) Note added by J. W. Kysar (4 November 1997). This UMAT has beenc modified to keep track of the cumulative shear strain in each c individual slip system. This information is needed to correct anc error in the implementation of the Bassani and Wu hardening law.c Any line of code which has been added or modified is precededc 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.cc The hardening law by Bassani and Wu was implemented incorrectly.c This law is a function of both hyperbolic secant squared and hyperbolicc tangent. However, the arguments of sech and tanh are related to the *total*c slip on individual slip systems. Formerly, the UMAT implemented thisc 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 byC (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.CC- Subroutines:CC ROTATION - forming rotation matrix, i.e. the direction C cosines of cubic crystal 100, 010 and 001C directions in global system at the initial C stateCC 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 stateCC GSLPINIT - calculating initial value of current strengths C at initial stateCC STRAINRATE - based on current values of resolved shear C stresses and current strength, calculating C shear strain-rates in slip systemsCC LATENTHARDEN - forming self- and latent-hardening matrix CC ITERATION - generating arrays for the Newton-Rhapson C iterationCC LUDCMP - LU decompositionCC LUBKSB - linear equation solver based on LU C decomposition method (must call LUDCMP first)C- Function subprogram:C F - shear strain-rates in slip systemsC- Variables:CC STRESS - stresses (INPUT & OUTPUT)C Cauchy stresses for finite deformationC STATEV - solution dependent state variables (INPUT & OUTPUT)C DDSDDE - Jacobian matrix (OUTPUT)C- Variables passed in for information:CC STRAN - strainsC logarithmic strain for finite deformation C (actually, integral of the symmetric part of velocityC gradient with respect to time)C DSTRAN - increments of strainsC CMNAME - name given in the *MATERIAL optionC NDI - number of direct stress componentsC NSHR - number of engineering shear stress componentsC NTENS - NDI+NSHRC NSTATV - number of solution dependent state variables (as C defined in the *DEPVAR option)C PROPS - material constants entered in the *USER MATERIAL C optionC NPROPS - number of material constantsCC- This subroutine provides the plastic constitutive relation of C single crystals for finite element code ABAQUS. The plastic slipC 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 isC assumed a function of the ratio of corresponding resolved shear C stress over current strength, and of the time step. The resolvedC shear stress is the double product of stress tensor with the slipC 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 optionC 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.CC- Important notice:CC (1) The number of state variables NSTATV must be larger than (or CFIX equal to) TEN (10) times the total number of slip systems inC all sets, NSLPTL, plus FIVE (5)CFIX NSTATV = 10 * NSLPTL + 5C 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 setC could be either 6, 12, 24 or 48 for a cubic crystal, e.g. 12 C for 110.CC 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 fiveCFIX NSTATV = NPARMT + 10 * NSLPTL + 5CC (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 largerC than or equal to the total number of slip systems in all sets. C For example, if 110 is the only set of slip system C potentially activated, ND could be taken as twelve (12). c include aba_param.incc 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 setC- 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 systemC- 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)/2C SLPSPN(2,i) - SLPDIR(3,i)*SLPNOR(1,i)-C SLPDIR(1,i)*SLPNOR(3,i)/2C SLPSPN(3,i) - SLPDIR(2,i)*SLPNOR(3,i)-C SLPDIR(3,i)*SLPNOR(2,i)/2C where index i corresponds to the ith slip systemC- DSPDIR - increments of slip directionsC- DSPNOR - increments of normals to slip planesCC- DLOCAL - elastic matrix in local cubic crystal systemC- D - elastic matrix in global systemC- ROTD - rotation matrix transforming DLOCAL to DCC- ROTATE - rotation matrix, direction cosines of 100, 010 C and 001 of cubic crystal in global systemCC- FSLIP - shear strain-rates in slip systemsC- DFDXSP - derivatives of FSLIP w.r.t x=TAUSLP/GSLIP, where C TAUSLP is the resolved shear stress and GSLIP is the C current strengthCC- 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 stressCC- H - self- and latent-hardening matrixC 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 systemC (i not equal j)CC- DDGDDE - derivatice of the shear strain increments in slip C systems w.r.t. the increment of strainsCC- DSTRES - Jaumann increments of stresses, i.e. corotational C stress-increments formed on axes spinning with the C materialC- DELATS - strain-increments associated with lattice stretchingC DELATS(1) - DELATS(3) - normal strain incrementsC DELATS(4) - DELATS(6) - engineering shear strain C incrementsC- DSPIN - spin-increments associated with the material elementC DSPIN(1) - component 12 of the spin tensorC DSPIN(2) - component 31 of the spin tensorC DSPIN(3) - component 23 of the spin tensorCC- DVGRAD - increments of deformation gradient in the current C state, i.e. velocity gradient times the increment of C timeCC- DGAMMA - increment of shear strains in slip systemsC- DTAUSP - increment of resolved shear stresses in slip systems C- DGSLIP - increment of current strengths in slip systemsCCC- Arrays for iteration:CC FSLIP1, STRES1, GAMMA1, TAUSP1, GSLP1 , SPNOR1, SPDIR1, C DDSDE1, DSOLD , DGAMOD, DTAUOD, DGSPOD, DSPNRO, DSPDRO,C DHDGDGCCC- Solution dependent state variable STATEV:C Denote the number of total slip systems by NSLPTL, which C will be calculated in this code.CC Array STATEV:C 1 - NSLPTL : current strength in slip systemsC NSLPTL+1 - 2*NSLPTL : shear strain in slip systemsC 2*NSLPTL+1 - 3*NSLPTL : resolved shear stress in slip systemsCC 3*NSLPTL+1 - 6*NSLPTL : current components of normals to slipC slip planesC 6*NSLPTL+1 - 9*NSLPTL : current components of slip directionsCCFIX 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)CFIXCFIX 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)CCFIX 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).CC NSTATV-3 : number of slip systems in the 1st setC NSTATV-2 : number of slip systems in the 2nd setC NSTATV-1 : number of slip systems in the 3rd setC NSTATV : total number of slip systems in all C setsCCC- Material constants PROPS:CC PROPS(1) - PROPS(21) - elastic constants for a general elasticC anisotropic materialCC isotropic : PROPS(i)=0 for i2C PROPS(1) - Youngs modulusC PROPS(2) - Poissons ratioCC cubic : PROPS(i)=0 for i3C PROPS(1) - c11C PROPS(2) - c12C PROPS(3) - c44CC orthotropic : PORPS(i)=0 for i9C PROPS(1) - PROPS(9) are D1111, D1122, D2222,C D1133, D2233, D3333, D1212, D1313, D2323, C respectively, which has the same definition C as ABAQUS for orthotropic materialsC (see *ELASTIC card)CC anisotropic : PROPS(1) - PROPS(21) are D1111, D1122, C D2222, D1133, D2233, D3333, D1112, D2212, C D3312, D1212, D1113, D2213, D3313, D1213, C D1313, D1123, D2223, D3323, D1223, D1323, C D2323, respectively, which has the same C definition as ABAQUS for anisotropic C materials (see *ELASTIC card)CCC PROPS(25) - PROPS(56) - parameters characterizing all slip C systems to be activated in a cubic C crystalCC PROPS(25) - number of sets of slip systems (maximum 3), C e.g. (110)1-11 and (101)11-1 are in the C same set of slip systems, (110)1-11 and C (121)1-11 belong to different sets of slip C systemsC (It must be a real number, e.g. 3., not 3 !)CC PROPS(33) - PROPS(35) - normal to a typical slip plane inC the first set of slip systems, C e.g. (1 1 0)C (They must be real numbers, e.g. C 1. 1. 0., not 1 1 0 !)C PROPS(36) - PROPS(38) - a typical slip direction in the C first set of slip systems, e.g. C 1 1 1C (They must be real numbers, e.g. C 1. 1. 1., not 1 1 1 !)CC PROPS(41) - PROPS(43) - normal to a typical slip plane inC the second set of slip systemsC (real numbers)C PROPS(44) - PROPS(46) - a typical slip direction in the C second set of slip systemsC (real numbers)CC PROPS(49) - PROPS(51) - normal to a typical slip plane inC the third set of slip systemsC (real numbers)C PROPS(52) - PROPS(54) - a typical slip direction in the C third set of slip systemsC (real numbers)CCC PROPS(57) - PROPS(72) - parameters characterizing the initial C orientation of a single crystal in C global systemC The directions in global system and directions in local C cubic crystal system of two nonparallel vectors are neededC to determine the crystal orientation.CC PROPS(57) - PROPS(59) - p1 p2 p3, direction of first C vector in local cubic crystal C system, e.g. 1 1 0C (They must be real numbers, e.g. C 1. 1. 0., not 1 1 0 !)C PROPS(60) - PROPS(62) - P1 P2 P3, direction of first C vector in global system, e.g. C 2. 1. 0.C (It does not have to be a unit C vector)CC PROPS(65) - PROPS(67) - direction of second vector in C local cubic crystal system (real C numbers)C P
展开阅读全文