发布时间 : 星期日 文章黄永刚单晶塑性有限元umat子程序更新完毕开始阅读
C NITRTN = 0 --- no-iteration solution C
NITRTN=NITRTN+1
C----- Check whether the current stress state is the initial state IF (STATEV(1).EQ.0.) THEN
C----- Initial state C
C----- Generating the following parameters and variables at initial C state:
C Total number of slip systems in all the sets NSLPTL C Number of slip systems in each set NSLIP C Unit vectors in initial slip directions SLPDIR C Unit normals to initial slip planes SLPNOR C
NSLPTL=0 DO I=1,NSET
ISPNOR(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 DO
IF (ND.LT.NSLPTL) THEN WRITE (6,*)
2 '***ERROR - parameter ND chosen by the present user is less than 3 the total number of slip systems NSLPTL' STOP END IF
C----- Slip deformation tensor: SLPDEF (Schmid factors) DO J=1,NSLPTL
SLPDEF(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 DO
C----- Initial value of state variables: unit normal to a slip plane C and unit vector in a slip direction C
STATEV(NSTATV)=FLOAT(NSLPTL) DO I=1,NSET
STATEV(NSTATV-4+I)=FLOAT(NSLIP(I)) END DO
IDNOR=3*NSLPTL IDDIR=6*NSLPTL DO J=1,NSLPTL DO I=1,3
IDNOR=IDNOR+1
STATEV(IDNOR)=SLPNOR(I,J)
IDDIR=IDDIR+1
STATEV(IDDIR)=SLPDIR(I,J) END DO END DO
C----- Initial value of the current strength for all slip systems C
CALL GSLPINIT (STATEV(1), NSLIP, NSLPTL, NSET, PROPS(97))
C----- Initial value of shear strain in slip systems
CFIX-- Initial value of cumulative shear strain in each slip systems
DO I=1,NSLPTL
STATEV(NSLPTL+I)=0. CFIXA
STATEV(9*NSLPTL+I)=0. CFIXB
END DO
CFIXA
STATEV(10*NSLPTL+1)=0. CFIXB
C----- Initial value of the resolved shear stress in slip systems
DO I=1,NSLPTL TERM1=0.
DO J=1,NTENS
IF (J.LE.NDI) THEN
TERM1=TERM1+SLPDEF(J,I)*STRESS(J) ELSE
TERM1=TERM1+SLPDEF(J-NDI+3,I)*STRESS(J) END IF END DO
STATEV(2*NSLPTL+I)=TERM1 END DO
ELSE
C----- Current stress state C
C----- Copying from the array of state variables STATVE the following C parameters and variables at current stress state:
C Total number of slip systems in all the sets NSLPTL C Number of slip systems in each set NSLIP C Current slip directions SLPDIR
C Normals to current slip planes SLPNOR C
NSLPTL=NINT(STATEV(NSTATV)) DO I=1,NSET
NSLIP(I)=NINT(STATEV(NSTATV-4+I)) END DO
IDNOR=3*NSLPTL IDDIR=6*NSLPTL DO J=1,NSLPTL DO I=1,3
IDNOR=IDNOR+1
SLPNOR(I,J)=STATEV(IDNOR)
IDDIR=IDDIR+1
SLPDIR(I,J)=STATEV(IDDIR) END DO END DO
C----- Slip deformation tensor: SLPDEF (Schmid factors) DO J=1,NSLPTL
SLPDEF(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 DO
END IF
C----- Slip spin tensor: SLPSPN (only needed for finite rotation) IF (NLGEOM.NE.0) THEN DO J=1,NSLPTL
SLPSPN(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 DO END IF
C----- Double dot product of elastic moduli tensor with the slip C deformation tensor (Schmid factors) plus, only for finite C rotation, the dot product of slip spin tensor with the stress: C DDEMSD C
DO J=1,NSLPTL DO I=1,6
DDEMSD(I,J)=0. DO K=1,6
DDEMSD(I,J)=DDEMSD(I,J)+D(K,I)*SLPDEF(K,J) END DO END DO END DO
IF (NLGEOM.NE.0) THEN DO J=1,NSLPTL
DDEMSD(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) THEN
DDEMSD(4,J)=DDEMSD(4,J)+SLPSPN(1,J)*STRESS(2)