黄永刚单晶塑性有限元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)

联系客服:779662525#qq.com(#替换为@)