ABAQUS岩土工程实例：dyna.for文件操作

~熬过最苦的日子，做最酷的自己~

      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,

     1 RPL,DDSDDT,DRPLDE,DRPLDT,

     2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,

     3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,

     4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)

C

      INCLUDE 'ABA_PARAM.INC'

C

      CHARACTER*80 CMNAME

      DIMENSION STRESS(NTENS),STATEV(NSTATV),

     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),

     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),

     3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),

     4 JSTEP(4)

      DIMENSION DSTRES(6),D(3,3)

DIMENSION PROPSS(5),PS(3)

C

C  EVALUATE NEW STRESS TENSOR

C

      FEK=PROPS(1)

FEN=PROPS(2)

FEV=PROPS(3)

FEW=PROPS(4)

SIG30=STATEV(1)

GG=STATEV(2)

DD=STATEV(3)

GAMA=STATEV(4)

ELG=FEK*(SIG30/100)**FEN*100

ELG=ELG*GG

ELM=2*ELG*FEV/(1-2*FEV)

ELGT=2*ELG/FEW*DD

ELMT=2*ELM/FEW*DD

ELG1=ELG/1.01

ELM1=ELM/1.01

ELGT1=ELGT/1.01

ELMT1=ELMT/1.01

EFEI=0.01*ELGT/ELG

PROPSS(1)=ELM1

PROPSS(2)=ELG1

PROPSS(3)=ELMT1

PROPSS(4)=ELGT1

PROPSS(5)=EFEI

EV = 0.

      DEV = 0.

      DO K1=1,NDI

         EV = EV + STRAN(K1)

         DEV = DEV + DSTRAN(K1)

      END DO

C

TERM1 = .5*DTIME + PROPSS(5)

      TERM1I = 1./TERM1

      TERM2 = (.5*DTIME*PROPSS(1)+PROPSS(3))*TERM1I*DEV

      TERM3 = (DTIME*PROPSS(2)+2.*PROPSS(4))*TERM1I

C

      DO K1=1,NDI

         DSTRES(K1) = TERM2+TERM3*DSTRAN(K1)

     1     +DTIME*TERM1I*(PROPSS(1)*EV

     2     +2.*PROPSS(2)*STRAN(K1)-STRESS(K1))

         STRESS(K1) = STRESS(K1) + DSTRES(K1)

      END DO

C

      TERM2 = (.5*DTIME*PROPSS(2) + PROPSS(4))*TERM1I

      I1 = NDI

      DO K1=1,NSHR

         I1 = I1+1

         DSTRES(I1) = TERM2*DSTRAN(I1)+

     1     DTIME*TERM1I*(PROPSS(2)*STRAN(I1)-STRESS(I1))

         STRESS(I1) = STRESS(I1)+DSTRES(I1)

      END DO

C

C  CREATE NEW JACOBIAN

C

      TERM2 = (DTIME*(.5*PROPSS(1)+PROPSS(2))+PROPSS(3)+

     1  2.*PROPSS(4))*TERM1I

      TERM3 = (.5*DTIME*PROPSS(1)+PROPSS(3))*TERM1I

      DO K1=1,NTENS

         DO K2=1,NTENS

            DDSDDE(K2,K1) = 0.

         END DO

      END DO

C

      DO K1=1,NDI

         DDSDDE(K1,K1) = TERM2

      END DO

C

      DO K1=2,NDI

         N2=K1-1

         DO K2=1,N2

            DDSDDE(K2,K1) = TERM3

            DDSDDE(K1,K2) = TERM3

         END DO

      END DO

      TERM2 = (.5*DTIME*PROPSS(2)+PROPSS(4))*TERM1I

      I1 = NDI

      DO K1=1,NSHR

         I1 = I1+1

         DDSDDE(I1,I1) = TERM2

      END DO

C

C  TOTAL CHANGE IN SPECIFIC ENERGY

C

      TDE = 0.

      DO K1=1,NTENS

         TDE = TDE + (STRESS(K1)-.5*DSTRES(K1))*DSTRAN(K1)

      END DO

C

C  CHANGE IN SPECIFIC ELASTIC STRAIN ENERGY

C

      TERM1 = PROPSS(1) + 2.*PROPSS(2)

      DO  K1=1,NDI

         D(K1,K1) = TERM1

      END DO

      DO K1=2,NDI

         N2=K1-1

         DO K2=1,N2

            D(K1,K2) = PROPSS(1)

            D(K2,K1) = PROPSS(1)

         END DO

      END DO

      DEE = 0.

      DO K1=1,NDI

         TERM1 = 0.

         TERM2 = 0.

         DO K2=1,NDI

            TERM1 = TERM1 + D(K1,K2)*STRAN(K2)

            TERM2 = TERM2 + D(K1,K2)*DSTRAN(K2)

         END DO

         DEE = DEE + (TERM1+.5*TERM2)*DSTRAN(K1)

      END DO

      I1 = NDI

      DO K1=1,NSHR

         I1 = I1+1

         DEE = DEE + PROPSS(2)*(STRAN(I1)+.5*DSTRAN(I1))*DSTRAN(I1)

      END DO

      SSE = SSE + DEE

      SCD = SCD+TDE-DEE

CALL SPRINC(STRAN,PS,2,3,3)

GAMAX=((PS(1)-PS(2))**2+(PS(2)-PS(3))**2+(PS(3)-PS(1))**2)*2/3

GAMAX=SQRT(GAMAX)

IF(GAMA.LT.GAMAX)STATEV(4)=GAMAX

      RETURN

      END