abaqus邓肯张模型umat

首先是始点刚度法：

  1       SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
  2      1 RPL,DDSDDT,DRPLDE,DRPLDT,
  3      2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
  4      3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
  5      4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
  6 C
  7       INCLUDE 'ABA_PARAM.INC'
  8 C    始点刚度法
  9       CHARACTER*80 CMNAME
 10       DIMENSION STRESS(NTENS),STATEV(NSTATV),
 11      1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
 12      2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
 13      3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
 14      4 JSTEP(4)
 15     DIMENSION PS(3),DSTRESS(NTENS)
 16     PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)
 17 
 18       EK=PROPS(1)
 19     EN=PROPS(2)
 20     RF=PROPS(3)
 21     C=PROPS(4)
 22     FAI=PROPS(5)/180.0*3.1415926
 23     UG=PROPS(6)
 24     UD=PROPS(7)
 25     UF=PROPS(8)
 26     EKUR=PROPS(9)
 27     PA=PROPS(10)
 28     DFAI=PROPS(11)/180.0*3.1415926
 29     S1S3O=STATEV(1)
 30     S3O=STATEV(2)
 31     SSS=STATEV(3)
 32     CALL GETPS(STRESS,PS,NTENS)
 33    
 34     FAI=FAI-DFAI*LOG10(S3O/PA)
 35     CALL GETEMOD(PS,EK,EN,RF,C,FAI,ENU,PA,EKUR,EMOD,S,S3O,UG,UD,UF
 36      1    ,SSS,S1S3O)
 37     EBULK3=EMOD/(ONE-TWO*ENU)
 38       EG2=EMOD/(ONE+ENU)
 39       EG=EG2/TWO
 40       EG3=THREE*EG
 41       ELAM=(EBULK3-EG2)/THREE
 42     CALL GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG)
 43     DSTRESS=0.0
 44     CALL GETSTRESS(DDSDDE,DSTRESS,DSTRAN,NTENS)
 45       
 46     DO 701 I1=1,NTENS
 47     STRESS(I1)=STRESS(I1)+DSTRESS(I1)
 48 701    CONTINUE
 49     CALL GETPS(STRESS,PS,NTENS)
 50     CALL GETEMOD(PS,EK,EN,RF,C,FAI,ENU,PA,EKUR,EMOD,S,S3O,UG,UD,UF,
 51      1    SSS,S1S3O)
 52     EBULK3=EMOD/(ONE-TWO*ENU)
 53       EG2=EMOD/(ONE+ENU)
 54       EG=EG2/TWO
 55       EG3=THREE*EG
 56       ELAM=(EBULK3-EG2)/THREE
 57     CALL GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG)
 58     IF(PS(3).GT.S3O)S3O=PS(3)
 59      IF((PS(1)-PS(3)).GT.S1S3O)S1S3O=PS(1)-PS(3)
 60     IF(S.GT.SSS)SSS=S
 61     STATEV(1)=S1S3O
 62     STATEV(2)=S3O
 63     STATEV(3)=SSS
 64      END
 65         
 66     SUBROUTINE GETPS(STRESS,PS,NTENS)
 67     INCLUDE 'ABA_PARAM.INC'
 68     DIMENSION PS(3),STRESS(NTENS)
 69     CALL SPRINC(STRESS,PS,1,3,3)
 70     DO 310 I=1,2
 71     DO 320 J=I+1,3
 72     IF(PS(I).GT.PS(J))THEN
 73     PPS=PS(I)
 74     PS(I)=PS(J)
 75     PS(J)=PPS
 76     END IF
 77 320    CONTINUE
 78 310    CONTINUE
 79     DO 330 K1=1,3
 80     PS(K1)=-PS(K1)
 81 330    CONTINUE
 82     RETURN
 83     END 
 84 
 85     SUBROUTINE GETEMOD(PS,EK,EN,RF,C,FAI,ENU,PA,EKUR,EMOD,S,S3O
 86      1