ABAQUS岩土工程实例：duncan-2.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 PS(3),DSTRESS(NTENS),TDSTRESS(NTENS),TSTRESS(NTENS)

PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)

      EK=PROPS(1)

EN=PROPS(2)

RF=PROPS(3)

C=PROPS(4)

FAI=PROPS(5)/180.0*3.1415926

UG=PROPS(6)

UD=PROPS(7)

UF=PROPS(8)

EKUR=PROPS(9)

PA=PROPS(10)

DFAI=PROPS(11)/180.0*3.1415926

S1S3O=STATEV(1)

S3O=STATEV(2)

SSS=STATEV(3)

CALL GETPS(STRESS,PS,NTENS)

FAI=FAI-DFAI*LOG10(S3O/PA)

CALL GETEMOD(PS,EK,EN,RF,C,FAI,ENU,PA,EKUR,EMOD,S,S3O,UG,UD,UF

     1 ,SSS,S1S3O)

EBULK3=EMOD/(ONE-TWO*ENU)

      EG2=EMOD/(ONE+ENU)

      EG=EG2/TWO

      EG3=THREE*EG

      ELAM=(EBULK3-EG2)/THREE

CALL GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG)

TDSTRESS=0.0

CALL GETSTRESS(DDSDDE,TDSTRESS,DSTRAN,NTENS)

DO 701 I1=1,NTENS

TSTRESS(I1)=STRESS(I1)+TDSTRESS(I1)*0.5

701 CONTINUE

CALL GETPS(TSTRESS,PS,NTENS)

CALL GETEMOD(PS,EK,EN,RF,C,FAI,ENU,PA,EKUR,EMOD,S,S3O,UG,UD,UF,

     1 SSS,S1S3O)

EBULK3=EMOD/(ONE-TWO*ENU)

      EG2=EMOD/(ONE+ENU)

      EG=EG2/TWO

      EG3=THREE*EG

      ELAM=(EBULK3-EG2)/THREE

CALL GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG)

DSTRESS=0.0

CALL GETSTRESS(DDSDDE,DSTRESS,DSTRAN,NTENS)

DO 702 I1=1,NTENS

  STRESS(I1)=STRESS(I1)+DSTRESS(I1)

702 CONTINUE

  CALL GETPS(STRESS,PS,NTENS)

  CALL GETEMOD(PS,EK,EN,RF,C,FAI,ENU,PA,EKUR,EMOD,S,S3O,UG,UD,UF,

     1 SSS,S1S3O)

EBULK3=EMOD/(ONE-TWO*ENU)

      EG2=EMOD/(ONE+ENU)

      EG=EG2/TWO

      EG3=THREE*EG

      ELAM=(EBULK3-EG2)/THREE

CALL GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG)

IF(PS(3).GT.S3O)S3O=PS(3)

  IF((PS(1)-PS(3)).GT.S1S3O)S1S3O=PS(1)-PS(3)

IF(S.GT.SSS)SSS=S

STATEV(1)=S1S3O

STATEV(2)=S3O

STATEV(3)=SSS

END

SUBROUTINE GETPS(STRESS,PS,NTENS)

INCLUDE 'ABA_PARAM.INC'

DIMENSION PS(3),STRESS(NTENS)

CALL SPRINC(STRESS,PS,1,3,3)

DO 310 I=1,2

DO 320 J=I+1,3

IF(PS(I).GT.PS(J))THEN

PPS=PS(I)

PS(I)=PS(J)

PS(J)=PPS

END IF

320 CONTINUE

310 CONTINUE

DO 330 K1=1,3

PS(K1)=-PS(K1)

330 CONTINUE

RETURN

END 

SUBROUTINE GETEMOD(PS,EK,EN,RF,C,FAI,ENU,PA,EKUR,EMOD,S,S3O

     1 ,UG,UD,UF,SSS,S1S3O)

INCLUDE 'ABA_PARAM.INC'

DIMENSION PS(3)

S=(1-SIN(FAI))*(PS(1)-PS(3))

IF(PS(3).LT.(-C/TAN(FAI))) THEN

          S=0.99

        ELSE

          S=S/(2*C*COS(FAI)+2*PS(3)*SIN(FAI))

          IF(S.GE.0.99)S=0.99

        END IF

      EMOD=EK*PA*((S3O/PA)**EN)*((1-RF*S)**2)

      AA=UD*(PS(1)-PS(3))

AA=AA/(EK*PA*((S3O/PA)**EN))

AA=AA/(1-RF*S)

ENU=UG-UF*LOG10(S3O/PA)

ENU=ENU/(1-AA)/(1-AA)

IF(ENU.GT.0.49)ENU=0.49

      IF(ENU.LT.0.05)ENU=0.05

IF(S.LT.SSS.AND.(PS(1)-PS(3)).LT.S1S3O)THEN

EMOD=EKUR*PA*((S3O/PA)**EN)

END IF

END  

SUBROUTINE GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG)

INCLUDE 'ABA_PARAM.INC'

DIMENSION DDSDDE(NTENS,NTENS)

DO 20 K1=1,NTENS

        DO 10 K2=1,NTENS

           DDSDDE(K2,K1)=0.0

 10     CONTINUE

 20   CONTINUE

      DO 40 K1=1,NDI

        DO 30 K2=1,NDI

           DDSDDE(K2,K1)=ELAM

 30     CONTINUE

        DDSDDE(K1,K1)=EG2+ELAM

 40   CONTINUE

      DO 50 K1=NDI+1,NTENS

        DDSDDE(K1,K1)=EG

 50   CONTINUE

RETURN

END

SUBROUTINE GETSTRESS(DDSDDE,STRESS,DSTRAN,NTENS)

INCLUDE 'ABA_PARAM.INC'

DIMENSION DDSDDE(NTENS,NTENS),STRESS(NTENS),DSTRAN(NTENS)

DO 70 K1=1,NTENS

        DO 60 K2=1,NTENS

           STRESS(K1)=STRESS(K1)+DDSDDE(K1,K2)*DSTRAN(K2)

 60     CONTINUE

 70   CONTINUE

RETURN

END