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

CALL LINE (J3, J1, J2, IWKNOR(1,5)) CALL LINE (J2, J3, J1, IWKNOR(1,9)) CALL LINE (J1, J3, J2, IWKNOR(1,13)) CALL LINE (J2, J1, J3, IWKNOR(1,17)) CALL LINE (J3, J2, J1, IWKNOR(1,21))

END IF

C----- Generating all slip systems in this set C

C----- Unit vectors in slip directions: SLPDIR, and unit normals to C slip planes: SLPNOR in local cubic crystal system C

WRITE (6,*) ' '

WRITE (6,*) ' # Slip plane Slip direction'

DO J=1,NSPNOR DO I=1,NSPDIR

IDOT=0 DO K=1,3

IDOT=IDOT+IWKDIR(K,I)*IWKNOR(K,J) END DO

IF (IDOT.EQ.0) THEN NSLIP=NSLIP+1 DO K=1,3

SLPDIR(K,NSLIP)=IWKDIR(K,I)/RMODIR SLPNOR(K,NSLIP)=IWKNOR(K,J)/RMONOR END DO

WRITE (6,10) NSLIP,

2 (IWKNOR(K,J),K=1,3), (IWKDIR(K,I),K=1,3)

END IF

END DO END DO

10 FORMAT(1X,I2,9X,'(',3(1X,I2),1X,')',10X,'[',3(1X,I2),1X,']')

WRITE (6,*) 'Number of slip systems in this set = ',NSLIP WRITE (6,*) ' '

IF (NSLIP.EQ.0) THEN

WRITE (6,*)

* 'There is no slip direction normal to the slip planes!' STOP

ELSE

C----- Unit vectors in slip directions: SLPDIR, and unit normals to C slip planes: SLPNOR in global system C

DO J=1,NSLIP DO I=1,3

TERM(I)=0. DO K=1,3

TERM(I)=TERM(I)+ROTATE(I,K)*SLPDIR(K,J) END DO END DO DO I=1,3

SLPDIR(I,J)=TERM(I) END DO

DO I=1,3

TERM(I)=0. DO K=1,3

TERM(I)=TERM(I)+ROTATE(I,K)*SLPNOR(K,J) END DO END DO DO I=1,3

SLPNOR(I,J)=TERM(I) END DO END DO

END IF

RETURN END

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

SUBROUTINE LINE (I1, I2, I3, IARRAY)

C----- Generating all possible slip directions (or slip planes C {lmn}) for a cubic crystal, where l,m,n are not zeros.

C----- Use single precision on cray C

IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IARRAY(3,4)

DO J=1,4

IARRAY(1,J)=I1 IARRAY(2,J)=I2 IARRAY(3,J)=I3 END DO

DO I=1,3 DO J=1,4

IF (J.EQ.I+1) IARRAY(I,J)=-IARRAY(I,J) END DO END DO

RETURN END

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

SUBROUTINE LINE1 (J1, J2, IARRAY, ID)

C----- Generating all possible slip directions <0mn> (or slip planes C {0mn}) for a cubic crystal, where m,n are not zeros and m does C not equal n.

C----- Use single precision on cray C

IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IARRAY(3,2)

IARRAY(ID,1)=0 IARRAY(ID,2)=0

ID1=ID+1

IF (ID1.GT.3) ID1=ID1-3 IARRAY(ID1,1)=J1 IARRAY(ID1,2)=J1

ID2=ID+2

IF (ID2.GT.3) ID2=ID2-3 IARRAY(ID2,1)=J2 IARRAY(ID2,2)=-J2

RETURN END

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

SUBROUTINE GSLPINIT (GSLIP0, NSLIP, NSLPTL, NSET, PROP)

C----- This subroutine calculates the initial value of current

C strength for each slip system in a rate-dependent single crystal. C Two sets of initial values, proposed by Asaro, Pierce et al, and C by Bassani, respectively, are used here. Both sets assume that C the initial values for all slip systems are the same (initially C isotropic).

C----- These initial values are assumed the same for all slip systems C in each set, though they could be different from set to set, e.g. C <110>{111} and <110>{100}.

C----- Users who want to use their own initial values may change the

C function subprogram GSLP0. The parameters characterizing these C initial values are passed into GSLP0 through array PROP.

C----- Use single precision on cray C

IMPLICIT REAL*8 (A-H,O-Z) EXTERNAL GSLP0

DIMENSION GSLIP0(NSLPTL), NSLIP(NSET), PROP(16,NSET)

C----- Function subprograms: C

C GSLP0 -- User-supplied function subprogram given the initial C value of current strength at initial state

C----- Variables: C

C GSLIP0 -- initial value of current strength (OUTPUT) C