首先是始点刚度法：
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 ,UG,UD,UF,SSS,S1S3O) 87 INCLUDE 'ABA_PARAM.INC' 88 DIMENSION PS(3) 89 S=(1-SIN(FAI))*(PS(1)-PS(3)) 90 IF(PS(3).LT.(-C/TAN(FAI))) THEN 91 S=0.99 92 ELSE 93 94 S=S/(2*C*COS(FAI)+2*PS(3)*SIN(FAI)) 95 IF(S.GE.0.99)then 96 S=0.99 97 end if 98 END IF 99 EMOD=EK*PA*((S3O/PA)**EN)*((1-RF*S)**2) 100 101 AA=UD*(PS(1)-PS(3)) 102 AA=AA/(EK*PA*((S3O/PA)**EN)) 103 AA=AA/(1-RF*S) 104 ENU=UG-UF*LOG10(S3O/PA) 105 ENU=ENU/(1-AA)/(1-AA) 106 IF(ENU.GT.0.49)ENU=0.49 107 IF(ENU.LT.0.05)ENU=0.05 108 IF(S.LT.SSS.AND.(PS(1)-PS(3)).LT.S1S3O)THEN 109 EMOD=EKUR*PA*((S3O/PA)**EN) 110 END IF 111 END 112 113 114 SUBROUTINE GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG) 115 INCLUDE 'ABA_PARAM.INC' 116 DIMENSION DDSDDE(NTENS,NTENS) 117 DO 20 K1=1,NTENS 118 DO 10 K2=1,NTENS 119 DDSDDE(K2,K1)=0.0 120 10 CONTINUE 121 20 CONTINUE 122 DO 40 K1=1,NDI 123 DO 30 K2=1,NDI 124 DDSDDE(K2,K1)=ELAM 125 30 CONTINUE 126 DDSDDE(K1,K1)=EG2+ELAM 127 40 CONTINUE 128 DO 50 K1=NDI+1,NTENS 129 DDSDDE(K1,K1)=EG 130 50 CONTINUE 131 RETURN 132 END 133 134 SUBROUTINE GETSTRESS(DDSDDE,STRESS,DSTRAN,NTENS) 135 INCLUDE 'ABA_PARAM.INC' 136 DIMENSION DDSDDE(NTENS,NTENS),STRESS(NTENS),DSTRAN(NTENS) 137 DO 70 K1=1,NTENS 138 DO 60 K2=1,NTENS 139 STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1) 140 60 CONTINUE 141 70 CONTINUE 142 RETURN 143 END