黄永刚单晶塑性有限元umat子程序

IJ2=I+J-2

IF (MOD(IJ2,2).EQ.1) THEN TERM1=1. ELSE

TERM1=-1. END IF

DVGRAD(I,J)=DVGRAD(I,J)+TERM1*DSPIN(IJ2) DVGRAD(J,I)=DVGRAD(J,I)-TERM1*DSPIN(IJ2)

DO K=1,NSLPTL

DVGRAD(I,J)=DVGRAD(I,J)-TERM1*DGAMMA(K)* 2 SLPSPN(IJ2,K)

DVGRAD(J,I)=DVGRAD(J,I)+TERM1*DGAMMA(K)* 2 SLPSPN(IJ2,K) END DO END IF

END DO END DO

END IF

C----- Increment of resolved shear stress in a slip system: DTAUSP DO I=1,NSLPTL DTAUSP(I)=0. DO J=1,6

DTAUSP(I)=DTAUSP(I)+DDEMSD(J,I)*DELATS(J) END DO END DO

C----- Update the resolved shear stress in a slip system: C STATEV(2*NSLPTL+1) - STATEV(3*NSLPTL) C

DO I=1,NSLPTL

STATEV(2*NSLPTL+I)=STATEV(2*NSLPTL+I)+DTAUSP(I)-DTAUOD(I) END DO

C----- Increment of stress: DSTRES IF (NLGEOM.EQ.0) THEN DO I=1,NTENS DSTRES(I)=0. END DO ELSE

DO I=1,NTENS

DSTRES(I)=-STRESS(I)*DEV END DO END IF

DO I=1,NDI DO J=1,NDI

DSTRES(I)=DSTRES(I)+D(I,J)*DSTRAN(J) END DO

IF (NSHR.GT.0) THEN DO J=1,NSHR

DSTRES(I)=DSTRES(I)+D(I,J+3)*DSTRAN(J+NDI) END DO END IF

DO J=1,NSLPTL

DSTRES(I)=DSTRES(I)-DDEMSD(I,J)*DGAMMA(J) END DO END DO

IF (NSHR.GT.0) THEN DO I=1,NSHR

DO J=1,NDI

DSTRES(I+NDI)=DSTRES(I+NDI)+D(I+3,J)*DSTRAN(J) END DO

DO J=1,NSHR

DSTRES(I+NDI)=DSTRES(I+NDI)+D(I+3,J+3)*DSTRAN(J+NDI) END DO

DO J=1,NSLPTL

DSTRES(I+NDI)=DSTRES(I+NDI)-DDEMSD(I+3,J)*DGAMMA(J) END DO

END DO END IF

C----- Update the stress: STRESS DO I=1,NTENS

STRESS(I)=STRESS(I)+DSTRES(I)-DSOLD(I) END DO

C----- Increment of normal to a slip plane and a slip direction (only C needed for finite rotation) C

IF (NLGEOM.NE.0) THEN DO J=1,NSLPTL DO I=1,3

DSPNOR(I,J)=0. DSPDIR(I,J)=0.

DO K=1,3

DSPNOR(I,J)=DSPNOR(I,J)-SLPNOR(K,J)*DVGRAD(K,I) DSPDIR(I,J)=DSPDIR(I,J)+SLPDIR(K,J)*DVGRAD(I,K) END DO

END DO END DO

C----- Update the normal to a slip plane and a slip direction (only C needed for finite rotation) C

IDNOR=3*NSLPTL IDDIR=6*NSLPTL DO J=1,NSLPTL DO I=1,3

IDNOR=IDNOR+1

STATEV(IDNOR)=STATEV(IDNOR)+DSPNOR(I,J)-DSPNRO(I,J)

IDDIR=IDDIR+1

STATEV(IDDIR)=STATEV(IDDIR)+DSPDIR(I,J)-DSPDRO(I,J) END DO END DO

END IF

C----- Derivative of shear strain increment in a slip system w.r.t. C strain increment: DDGDDE C

TERM1=THETA*DTIME DO I=1,NTENS DO J=1,NSLPTL

TAUSLP=STATEV(2*NSLPTL+J) GSLIP=STATEV(J) X=TAUSLP/GSLIP

TERM2=TERM1*DFDXSP(J)/GSLIP

IF (I.LE.NDI) THEN

DDGDDE(J,I)=TERM2*DDEMSD(I,J) ELSE

DDGDDE(J,I)=TERM2*DDEMSD(I-NDI+3,J) END IF END DO

CALL LUBKSB (WORKST, NSLPTL, ND, INDX, DDGDDE(1,I))

END DO

C----- Derivative of stress increment w.r.t. strain increment, i.e. C Jacobian matrix C

C----- Jacobian matrix: elastic part DO J=1,NTENS DO I=1,NTENS DDSDDE(I,J)=0. END DO END DO

DO J=1,NDI DO I=1,NDI

DDSDDE(I,J)=D(I,J)

IF (NLGEOM.NE.0) DDSDDE(I,J)=DDSDDE(I,J)-STRESS(I) END DO END DO

IF (NSHR.GT.0) THEN DO J=1,NSHR DO I=1,NSHR

DDSDDE(I+NDI,J+NDI)=D(I+3,J+3) END DO

DO I=1,NDI

DDSDDE(I,J+NDI)=D(I,J+3) DDSDDE(J+NDI,I)=D(J+3,I) IF (NLGEOM.NE.0)

2 DDSDDE(J+NDI,I)=DDSDDE(J+NDI,I)-STRESS(J+NDI) END DO END DO END IF

C----- Jacobian matrix: plastic part (slip)