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

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

DO K=1,NSLPTL

DDSDDE(I,J)=DDSDDE(I,J)-DDEMSD(I,K)*DDGDDE(K,J) END DO END DO END DO

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

DO I=1,NSHR

DO K=1,NSLPTL

DDSDDE(I+NDI,J+NDI)=DDSDDE(I+NDI,J+NDI)-

2 DDEMSD(I+3,K)*DDGDDE(K,J+NDI) END DO END DO

DO I=1,NDI

DO K=1,NSLPTL

DDSDDE(I,J+NDI)=DDSDDE(I,J+NDI)-

2 DDEMSD(I,K)*DDGDDE(K,J+NDI) DDSDDE(J+NDI,I)=DDSDDE(J+NDI,I)-

2 DDEMSD(J+3,K)*DDGDDE(K,I) END DO END DO

END DO END IF

IF (ITRATN.NE.0) THEN DO J=1,NTENS DO I=1,NTENS

DDSDDE(I,J)=DDSDDE(I,J)/(1.+DEV) END DO END DO END IF

C----- Iteration ?

IF (ITRATN.NE.0) THEN

C----- Save solutions (without iteration):

C Shear strain-rate in a slip system FSLIP1 C Current strength in a slip system GSLP1

C Shear strain in a slip system GAMMA1

C Resolved shear stress in a slip system TAUSP1 C Normal to a slip plane SPNOR1 C Slip direction SPDIR1 C Stress STRES1

C Jacobian matrix DDSDE1 C

IF (NITRTN.EQ.0) THEN

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

FSLIP1(J)=FSLIP(J) GSLP1(J)=STATEV(J)

GAMMA1(J)=STATEV(NSLPTL+J) TAUSP1(J)=STATEV(2*NSLPTL+J) DO I=1,3

IDNOR=IDNOR+1

SPNOR1(I,J)=STATEV(IDNOR)

IDDIR=IDDIR+1

SPDIR1(I,J)=STATEV(IDDIR) END DO END DO

DO J=1,NTENS

STRES1(J)=STRESS(J) DO I=1,NTENS

DDSDE1(I,J)=DDSDDE(I,J) END DO END DO

END IF

C----- Increments of stress DSOLD, and solution dependent state

C variables DGAMOD, DTAUOD, DGSPOD, DSPNRO, DSPDRO (for the next C iteration) C

DO I=1,NTENS

DSOLD(I)=DSTRES(I) END DO

DO J=1,NSLPTL

DGAMOD(J)=DGAMMA(J)

DTAUOD(J)=DTAUSP(J) DGSPOD(J)=DGSLIP(J) DO I=1,3

DSPNRO(I,J)=DSPNOR(I,J) DSPDRO(I,J)=DSPDIR(I,J) END DO END DO

C----- Check if the iteration solution converges IDBACK=0 ID=0

DO I=1,NSET

DO J=1,NSLIP(I) ID=ID+1

X=STATEV(2*NSLPTL+ID)/STATEV(ID)

RESIDU=THETA*DTIME*F(X,PROPS(65+8*I))+DTIME*(1.0-THETA)* 2 FSLIP1(ID)-DGAMMA(ID)

IF (ABS(RESIDU).GT.GAMERR) IDBACK=1 END DO END DO

IF (IDBACK.NE.0.AND.NITRTN.LT.ITRMAX) THEN C----- Iteration: arrays for iteration CFIXA

CALL ITERATION (STATEV(NSLPTL+1), STATEV(2*NSLPTL+1), 2 STATEV(1), STATEV(9*NSLPTL+1), 3 STATEV(10*NSLPTL+1), NSLPTL,

4 NSET, NSLIP, ND, PROPS(97), DGAMOD, 5 DHDGDG) CFIXB

GO TO 1000

ELSE IF (NITRTN.GE.ITRMAX) THEN

C----- Solution not converge within maximum number of iteration (the C solution without iteration will be used) C

DO J=1,NTENS

STRESS(J)=STRES1(J) DO I=1,NTENS

DDSDDE(I,J)=DDSDE1(I,J) END DO END DO

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

STATEV(J)=GSLP1(J)

STATEV(NSLPTL+J)=GAMMA1(J) STATEV(2*NSLPTL+J)=TAUSP1(J)

DO I=1,3

IDNOR=IDNOR+1

STATEV(IDNOR)=SPNOR1(I,J)

IDDIR=IDDIR+1

STATEV(IDDIR)=SPDIR1(I,J) END DO END DO

END IF

END IF

C----- Total cumulative shear strains on all slip systems (sum of the C absolute values of shear strains in all slip systems)

CFIX-- Total cumulative shear strains on each slip system (sum of the CFIX absolute values of shear strains in each individual slip system) C

DO I=1,NSLPTL CFIXA

STATEV(10*NSLPTL+1)=STATEV(10*NSLPTL+1)+ABS(DGAMMA(I)) STATEV(9*NSLPTL+I)=STATEV(9*NSLPTL+I)+ABS(DGAMMA(I)) CFIXB

END DO

RETURN END

C----------------------------------------------------------------------

SUBROUTINE ROTATION (PROP, ROTATE)

C----- This subroutine calculates the rotation matrix, i.e. the C direction cosines of cubic crystal [100], [010] and [001] C directions in global system