黄永刚单晶塑性有限元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) cc (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, sincec 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 eachc 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 linec beginning CFIXB. Any comment line added or modified will beginc 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. Thereinc lay the problem. The UMAT did not restrict the current slip to be ac positive value. So when a slip with a negative sign was encountered, thec term containing tanh led to a negative hardening rate (since tanh is anc odd function).c The UMA T has been fixed by adding state variables to keep track of thec *total* slip on each slip system by integrating up the absolute valuec of slip rates for each individual slip system. These "solution dependentc variables" are available for postprocessing. The only required changec in the input file is that the DEPV AR command must be changed.cC----- 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 directionC cosines of cubic crystal [100], [010] and [001]C directions in global system at the initialC stateCC SLIPSYS -- calculating number of slip systems, unitC vectors in slip directions and unit normals toC slip planes in a cubic crystal at the initialC stateCC GSLPINIT -- calculating initial value of current strengthsC at initial stateCC STRAINRA TE -- based on current values of resolved shearC stresses and current strength, calculatingC shear strain-rates in slip systemsCC LATENTHARDEN -- forming self- and latent-hardening matrixCC ITERATION -- generating arrays for the Newton-RhapsonC iterationCC LUDCMP -- LU decompositionCC LUBKSB -- linear equation solver based on LUC 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 deformationC (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 (asC defined in the *DEPVAR option)C PROPS -- material constants entered in the *USER MA TERIAL C optionC NPROPS -- number of material constantsCC----- This subroutine provides the plastic constitutive relation ofC single crystals for finite element code ABAQUS. The plastic slipC of single crystal obeys the Schmid law. The program gives theC choice of small deformation theory and theory of finite rotationC and finite strain.C The strain increment is composed of elastic part and plasticC part. The elastic strain increment corresponds to latticeC stretching, the plastic part is the sum over all slip systems ofC plastic slip. The shear strain increment for each slip system isC assumed a function of the ratio of corresponding resolved shearC 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 currentC strength is related to shear strain increments over all slipC systems through self- and latent-hardening functions.C----- The implicit integration method proposed by Peirce, Shih andC Needleman (1984) is used here. The subroutine provides an option C of iteration to solve stresses and solution dependent stateC 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 andC SLIPSYS need to be modified to include the effect of crystalC 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 consideredC 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. 12C for {110}<111>.CC Users who need more parameters to characterize theC constitutive law of single crystal, e.g. the frameworkC proposed by Zarka, should make NSTATV larger than (or equal C to) the number of those parameters NPARMT plus nine timesC the total number of slip systems, NSLPTL, plus fiveCFIX NSTATV >= NPARMT + 10 * NSLPTL + 5CC (2) The tangent stiffness matrix in general is not symmetric ifC latent hardening is considered. Users must declare "UNSYMM"C in the input file, at the *USER MATERIAL card.CPARAMETER (ND=150)C----- The parameter ND determines the dimensions of the arrays inC this subroutine. The current choice 150 is a upper bound for aC 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 systemC potentially activated, ND could be taken as twelve (12).cinclude 'aba_param.inc'cCHARACTER*8 CMNAMEEXTERNAL Fdimension 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 initialC 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, whereC TAUSLP is the resolved shear stress and GSLIP is the C current strengthCC----- DDEMSD -- double dot product of the elastic moduli tensor withC the slip deformation tensor plus, only for finiteC rotation, the dot product of slip spin tensor withC the stressCC----- H -- self- and latent-hardening matrixC H(i,i) -- self hardening modulus of the ith slipC system (no sum over i)C H(i,j) -- latent hardening molulus of the ith slipC system due to a slip in the jth slip system C (i not equal j)CC----- DDGDDE -- derivatice of the shear strain increments in slipC systems w.r.t. the increment of strainsCC----- DSTRES -- Jaumann increments of stresses, i.e. corotationalC stress-increments formed on axes spinning with theC materialC----- DELATS -- strain-increments associated with lattice stretchingC DELATS(1) - DELATS(3) -- normal strain increments C 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 currentC state, i.e. velocity gradient times the increment ofC timeCC----- DGAMMA -- increment of shear strains in slip systemsC----- DTAUSP -- increment of resolved shear stresses in slip systemsC----- 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, whichC 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 eachCFIX slip system (sum of the absoluteCFIX values of shear strains in each slipCFIX system individually)CFIXCFIX 10*NSLPTL+1 : total cumulative shear strain on allC slip systems (sum of the absoluteC values of shear strains in all slipC systems)CCFIX 10*NSLPTL+2 - NSTA TV-4 : additional parameters users may needC to characterize the constitutive lawC of a single crystal (if there areC 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 allC setsCCC----- Material constants PROPS:CC PROPS(1) - PROPS(21) -- elastic constants for a general elasticC anisotropic materialCC isotropic : PROPS(i)=0 for i>2C PROPS(1) -- Young's modulusC PROPS(2) -- Poisson's ratioCC cubic : PROPS(i)=0 for i>3C PROPS(1) -- c11C PROPS(2) -- c12C PROPS(3) -- c44CC orthotropic : PORPS(i)=0 for i>9C PROPS(1) - PROPS(9) are D1111, D1122, D2222, C D1133, D2233, D3333, D1212, D1313, D2323,C respectively, which has the same definitionC 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 sameC definition as ABAQUS for anisotropicC materials (see *ELASTIC card)CCC PROPS(25) - PROPS(56) -- parameters characterizing all slipC systems to be activated in a cubicC crystalCC PROPS(25) -- number of sets of slip systems (maximum 3),C e.g. (110)[1-11] and (101)[11-1] are in theC same set of slip systems, (110)[1-11] andC (121)[1-11] belong to different sets of slipC 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 theC first set of slip systems, e.g.C [1 1 1]C (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 theC 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 theC third set of slip systemsC (real numbers)CCC PROPS(57) - PROPS(72) -- parameters characterizing the initialC orientation of a single crystal inC global systemC The directions in global system and directions in localC cubic crystal system of two nonparallel vectors are neededC to determine the crystal orientation.CC PROPS(57) - PROPS(59) -- [p1 p2 p3], direction of firstC vector in local cubic crystalC system, e.g. [1 1 0]C (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 firstC vector in global system, e.g.C [2. 1. 0.]C (It does not have to be a unitC vector)CC PROPS(65) - PROPS(67) -- direction of second vector inC local cubic crystal system (real C numbers)C PROPS(68) - PROPS(70) -- direction of second vector inC global systemCCC PROPS(73) - PROPS(96) -- parameters characterizing the visco-C plastic constitutive law (shearC strain-rate vs. resolved shearC stress), e.g. a power-law relationCC PROPS(73) - PROPS(80) -- parameters for the first set ofC slip systemsC PROPS(81) - PROPS(88) -- parameters for the second set of C slip systemsC PROPS(89) - PROPS(96) -- parameters for the third set ofC slip systemsCCC PROPS(97) - PROPS(144)-- parameters characterizing the self-C and latent-hardening laws of slipC systemsCC PROPS(97) - PROPS(104)-- self-hardening parameters for the C first set of slip systemsC PROPS(105)- PROPS(112)-- latent-hardening parameters for C the first set of slip systems and C interaction with other sets ofC slip systemsCC PROPS(113)- PROPS(120)-- self-hardening parameters for the C second set of slip systemsC PROPS(121)- PROPS(128)-- latent-hardening parameters for C the second set of slip systems C and interaction with other sets C of slip systemsCC PROPS(129)- PROPS(136)-- self-hardening parameters for the C third set of slip systemsC PROPS(137)- PROPS(144)-- latent-hardening parameters for C the third set of slip systems and C interaction with other sets ofC slip systemsCCC PROPS(145)- PROPS(152)-- parameters characterizing forward time C integration scheme and finiteC deformationCC PROPS(145) -- parameter theta controlling the implicitC integration, which is between 0 and 1C 0. : explicit integrationC 0.5 : recommended valueC 1. : fully implicit integrationCC PROPS(146) -- parameter NLGEOM controlling whether the C effect of finite rotation and finite strainC of crystal is considered,C 0. : small deformation theoryC otherwise : theory of finite rotation andC finite strainCCC PROPS(153)- PROPS(160)-- parameters characterizing iterationC methodCC PROPS(153) -- parameter ITRATN controlling whether theC iteration method is used,C 0. : no iterationC otherwise : iterationCC PROPS(154) -- maximum number of iteration ITRMAXCC PROPS(155) -- absolute error of shear strains in slipC systems GAMERRCC----- Elastic matrix in local cubic crystal system: DLOCALDO J=1,6DO I=1,6DLOCAL(I,J)=0.END DOEND DOCHECK=0.DO J=10,21CHECK=CHECK+ABS(PROPS(J))END DOIF (CHECK.EQ.0.) THENDO J=4,9CHECK=CHECK+ABS(PROPS(J))END DOIF (CHECK.EQ.0.) THENIF (PROPS(3).EQ.0.) THENC----- Isotropic materialGSHEAR=PROPS(1)/2./(1.+PROPS(2))E11=2.*GSHEAR*(1.-PROPS(2))/(1.-2.*PROPS(2))E12=2.*GSHEAR*PROPS(2)/(1.-2.*PROPS(2))DO J=1,3DLOCAL(J,J)=E11DO I=1,3IF (I.NE.J) DLOCAL(I,J)=E12END DODLOCAL(J+3,J+3)=GSHEAREND DOELSEC----- Cubic materialDO J=1,3DLOCAL(J,J)=PROPS(1)DO I=1,3IF (I.NE.J) DLOCAL(I,J)=PROPS(2)END DODLOCAL(J+3,J+3)=PROPS(3)END DOEND IFELSEC----- Orthotropic metarialDLOCAL(1,1)=PROPS(1)DLOCAL(1,2)=PROPS(2)DLOCAL(2,1)=PROPS(2)DLOCAL(2,2)=PROPS(3)DLOCAL(1,3)=PROPS(4)DLOCAL(3,1)=PROPS(4)DLOCAL(2,3)=PROPS(5)DLOCAL(3,2)=PROPS(5)DLOCAL(3,3)=PROPS(6)DLOCAL(4,4)=PROPS(7)DLOCAL(5,5)=PROPS(8)DLOCAL(6,6)=PROPS(9)END IFELSEC----- General anisotropic materialID=0DO J=1,6DO I=1,JID=ID+1DLOCAL(I,J)=PROPS(ID)DLOCAL(J,I)=DLOCAL(I,J)END DOEND DOEND IFC----- Rotation matrix: ROTA TE, i.e. direction cosines of [100], [010]C and [001] of a cubic crystal in global systemCCALL ROTATION (PROPS(57), ROTA TE)C----- Rotation matrix: ROTD to transform local elastic matrix DLOCAL C to global elastic matrix DCDO J=1,3J1=1+J/3J2=2+J/2DO I=1,3I1=1+I/3I2=2+I/2ROTD(I,J)=ROTATE(I,J)**2ROTD(I,J+3)=2.*ROTATE(I,J1)*ROTA TE(I,J2)ROTD(I+3,J)=ROTA TE(I1,J)*ROTATE(I2,J)ROTD(I+3,J+3)=ROTA TE(I1,J1)*ROTA TE(I2,J2)+2 ROTA TE(I1,J2)*ROTA TE(I2,J1)END DOEND DOC----- Elastic matrix in global system: DC {D} = {ROTD} * {DLOCAL} * {ROTD}transposeCDO J=1,6DO I=1,6D(I,J)=0.END DOEND DODO J=1,6DO I=1,JDO K=1,6DO L=1,6D(I,J)=D(I,J)+DLOCAL(K,L)*ROTD(I,K)*ROTD(J,L) END DOEND DOD(J,I)=D(I,J)END DOEND DOC----- Total number of sets of slip systems: NSETNSET=NINT(PROPS(25))IF (NSET.LT.1) THENWRITE (6,*) '***ERROR - zero sets of slip systems'STOPELSE IF (NSET.GT.3) THENWRITE (6,*)2 '***ERROR - more than three sets of slip systems'STOPEND IFC----- Implicit integration parameter: THETATHETA=PROPS(145)C----- Finite deformation ?C----- NLGEOM = 0, small deformation theoryC otherwise, theory of finite rotation and finite strain, Users C must declare "NLGEOM" in the input file, at the *STEP card CIF (PROPS(146).EQ.0.) THENNLGEOM=0ELSENLGEOM=1END IFC----- Iteration?C----- ITRATN = 0, no iterationC otherwise, iteration (solving increments of stresses andC solution dependent state variables)CIF (PROPS(153).EQ.0.) THENITRATN=0ELSEITRATN=1END IFITRMAX=NINT(PROPS(154))GAMERR=PROPS(155)NITRTN=-1DO I=1,NTENSDSOLD(I)=0.END DODO J=1,NDDGAMOD(J)=0.DTAUOD(J)=0.DGSPOD(J)=0.DO I=1,3DSPNRO(I,J)=0.DSPDRO(I,J)=0.END DOEND DOC----- Increment of spin associated with the material element: DSPIN C (only needed for finite rotation)CIF (NLGEOM.NE.0) THENDO J=1,3DO I=1,3TERM(I,J)=DROT(J,I)TRM0(I,J)=DROT(J,I)END DOTERM(J,J)=TERM(J,J)+1.D0TRM0(J,J)=TRM0(J,J)-1.D0END DOCALL LUDCMP (TERM, 3, 3, ITRM, DDCMP)DO J=1,3CALL LUBKSB (TERM, 3, 3, ITRM, TRM0(1,J)) END DODSPIN(1)=TRM0(2,1)-TRM0(1,2)DSPIN(2)=TRM0(1,3)-TRM0(3,1)DSPIN(3)=TRM0(3,2)-TRM0(2,3)END IFC----- Increment of dilatational strain: DEVDEV=0.D0DO I=1,NDIDEV=DEV+DSTRAN(I)END DOC----- Iteration starts (only when iteration method is used)1000 CONTINUEC----- Parameter NITRTN: number of iterationsC NITRTN = 0 --- no-iteration solutionCNITRTN=NITRTN+1C----- Check whether the current stress state is the initial stateIF (STATEV(1).EQ.0.) THENC----- Initial stateCC----- Generating the following parameters and variables at initialC state:C Total number of slip systems in all the sets NSLPTLC Number of slip systems in each set NSLIPC Unit vectors in initial slip directions SLPDIRC Unit normals to initial slip planes SLPNORCNSLPTL=0DO I=1,NSETISPNOR(1)=NINT(PROPS(25+8*I))ISPNOR(2)=NINT(PROPS(26+8*I))ISPNOR(3)=NINT(PROPS(27+8*I))ISPDIR(1)=NINT(PROPS(28+8*I))ISPDIR(2)=NINT(PROPS(29+8*I))ISPDIR(3)=NINT(PROPS(30+8*I))CALL SLIPSYS (ISPDIR, ISPNOR, NSLIP(I), SLPDIR(1,NSLPTL+1),2 SLPNOR(1,NSLPTL+1), ROTATE)NSLPTL=NSLPTL+NSLIP(I)END DOIF (ND.LT.NSLPTL) THENWRITE (6,*)2 '***ERROR - parameter ND chosen by the present user is less than3 the total number of slip systems NSLPTL'STOPEND IFC----- Slip deformation tensor: SLPDEF (Schmid factors)DO J=1,NSLPTLSLPDEF(1,J)=SLPDIR(1,J)*SLPNOR(1,J)SLPDEF(2,J)=SLPDIR(2,J)*SLPNOR(2,J)SLPDEF(3,J)=SLPDIR(3,J)*SLPNOR(3,J)SLPDEF(4,J)=SLPDIR(1,J)*SLPNOR(2,J)+SLPDIR(2,J)*SLPNOR(1,J)SLPDEF(5,J)=SLPDIR(1,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(1,J)SLPDEF(6,J)=SLPDIR(2,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(2,J) END DOC----- Initial value of state variables: unit normal to a slip planeC and unit vector in a slip directionCSTATEV(NSTATV)=FLOAT(NSLPTL)DO I=1,NSETSTA TEV(NSTA TV-4+I)=FLOAT(NSLIP(I))END DOIDNOR=3*NSLPTLIDDIR=6*NSLPTLDO J=1,NSLPTLDO I=1,3IDNOR=IDNOR+1STA TEV(IDNOR)=SLPNOR(I,J)IDDIR=IDDIR+1STA TEV(IDDIR)=SLPDIR(I,J)END DOEND DOC----- Initial value of the current strength for all slip systemsCCALL GSLPINIT (STATEV(1), NSLIP, NSLPTL, NSET, PROPS(97))C----- Initial value of shear strain in slip systemsCFIX-- Initial value of cumulative shear strain in each slip systemsDO I=1,NSLPTLSTA TEV(NSLPTL+I)=0.CFIXASTA TEV(9*NSLPTL+I)=0.CFIXBEND DOCFIXASTATEV(10*NSLPTL+1)=0.CFIXBC----- Initial value of the resolved shear stress in slip systemsDO I=1,NSLPTLTERM1=0.DO J=1,NTENSIF (J.LE.NDI) THENTERM1=TERM1+SLPDEF(J,I)*STRESS(J)ELSETERM1=TERM1+SLPDEF(J-NDI+3,I)*STRESS(J) END IFEND DOSTA TEV(2*NSLPTL+I)=TERM1END DOELSEC----- Current stress stateCC----- Copying from the array of state variables STA TVE the following C parameters and variables at current stress state:C Total number of slip systems in all the sets NSLPTLC Number of slip systems in each set NSLIPC Current slip directions SLPDIRC Normals to current slip planes SLPNORCNSLPTL=NINT(STA TEV(NSTATV))DO I=1,NSETNSLIP(I)=NINT(STATEV(NSTATV-4+I))END DOIDNOR=3*NSLPTLIDDIR=6*NSLPTLDO J=1,NSLPTLDO I=1,3IDNOR=IDNOR+1SLPNOR(I,J)=STATEV(IDNOR)IDDIR=IDDIR+1SLPDIR(I,J)=STA TEV(IDDIR)END DOEND DOC----- Slip deformation tensor: SLPDEF (Schmid factors)DO J=1,NSLPTLSLPDEF(1,J)=SLPDIR(1,J)*SLPNOR(1,J)SLPDEF(2,J)=SLPDIR(2,J)*SLPNOR(2,J)SLPDEF(3,J)=SLPDIR(3,J)*SLPNOR(3,J)SLPDEF(4,J)=SLPDIR(1,J)*SLPNOR(2,J)+SLPDIR(2,J)*SLPNOR(1,J)SLPDEF(5,J)=SLPDIR(1,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(1,J)SLPDEF(6,J)=SLPDIR(2,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(2,J) END DOEND IFC----- Slip spin tensor: SLPSPN (only needed for finite rotation)IF (NLGEOM.NE.0) THENDO J=1,NSLPTLSLPSPN(1,J)=0.5*(SLPDIR(1,J)*SLPNOR(2,J)-2 SLPDIR(2,J)*SLPNOR(1,J))SLPSPN(2,J)=0.5*(SLPDIR(3,J)*SLPNOR(1,J)-2 SLPDIR(1,J)*SLPNOR(3,J))SLPSPN(3,J)=0.5*(SLPDIR(2,J)*SLPNOR(3,J)-2 SLPDIR(3,J)*SLPNOR(2,J))END DOEND IFC----- Double dot product of elastic moduli tensor with the slipC deformation tensor (Schmid factors) plus, only for finiteC rotation, the dot product of slip spin tensor with the stress:C DDEMSDCDO J=1,NSLPTLDO I=1,6DDEMSD(I,J)=0.DO K=1,6DDEMSD(I,J)=DDEMSD(I,J)+D(K,I)*SLPDEF(K,J)END DOEND DOEND DOIF (NLGEOM.NE.0) THENDO J=1,NSLPTLDDEMSD(4,J)=DDEMSD(4,J)-SLPSPN(1,J)*STRESS(1)DDEMSD(5,J)=DDEMSD(5,J)+SLPSPN(2,J)*STRESS(1)IF (NDI.GT.1) THENDDEMSD(4,J)=DDEMSD(4,J)+SLPSPN(1,J)*STRESS(2)DDEMSD(6,J)=DDEMSD(6,J)-SLPSPN(3,J)*STRESS(2)END IFIF (NDI.GT.2) THENDDEMSD(5,J)=DDEMSD(5,J)-SLPSPN(2,J)*STRESS(3)DDEMSD(6,J)=DDEMSD(6,J)+SLPSPN(3,J)*STRESS(3)END IFIF (NSHR.GE.1) THENDDEMSD(1,J)=DDEMSD(1,J)+SLPSPN(1,J)*STRESS(NDI+1)DDEMSD(2,J)=DDEMSD(2,J)-SLPSPN(1,J)*STRESS(NDI+1)DDEMSD(5,J)=DDEMSD(5,J)-SLPSPN(3,J)*STRESS(NDI+1)DDEMSD(6,J)=DDEMSD(6,J)+SLPSPN(2,J)*STRESS(NDI+1) END IFIF (NSHR.GE.2) THENDDEMSD(1,J)=DDEMSD(1,J)-SLPSPN(2,J)*STRESS(NDI+2)DDEMSD(3,J)=DDEMSD(3,J)+SLPSPN(2,J)*STRESS(NDI+2)DDEMSD(4,J)=DDEMSD(4,J)+SLPSPN(3,J)*STRESS(NDI+2)DDEMSD(6,J)=DDEMSD(6,J)-SLPSPN(1,J)*STRESS(NDI+2) END IFIF (NSHR.EQ.3) THENDDEMSD(2,J)=DDEMSD(2,J)+SLPSPN(3,J)*STRESS(NDI+3)DDEMSD(3,J)=DDEMSD(3,J)-SLPSPN(3,J)*STRESS(NDI+3)DDEMSD(4,J)=DDEMSD(4,J)-SLPSPN(2,J)*STRESS(NDI+3)DDEMSD(5,J)=DDEMSD(5,J)+SLPSPN(1,J)*STRESS(NDI+3) END IFEND DOEND IFC----- Shear strain-rate in a slip system at the start of increment:C FSLIP, and its derivative: DFDXSPCID=1DO I=1,NSETIF (I.GT.1) ID=ID+NSLIP(I-1)CALL STRAINRATE (STATEV(NSLPTL+ID), STATEV(2*NSLPTL+ID),2 STA TEV(ID), NSLIP(I), FSLIP(ID), DFDXSP(ID),3 PROPS(65+8*I))END DOC----- Self- and latent-hardening laws。