首先是始点刚度法：
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
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 16 DIMENSION PS(3),DSTRESS(NTENS),TDSTRESS(NTENS),TSTRESS(NTENS) 17 PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0) 18 19 EK=PROPS(1) 20 EN=PROPS(2) 21 RF=PROPS(3) 22 C=PROPS(4) 23 FAI=PROPS(5)/180.0*3.1415926 24 UG=PROPS(6) 25 UD=PROPS(7) 26 UF=PROPS(8) 27 EKUR=PROPS(9) 28 PA=PROPS(10) 29 DFAI=PROPS(11)/180.0*3.1415926 30 S1S3O=STATEV(1) 31 S3O=STATEV(2) 32 SSS=STATEV(3) 33 CALL GETPS(STRESS,PS,NTENS) 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 TDSTRESS=0.0 44 CALL GETSTRESS(DDSDDE,TDSTRESS,DSTRAN,NTENS) 45 DO 701 I1=1,NTENS 46 TSTRESS(I1)=STRESS(I1)+TDSTRESS(I1)*0.5 47 701 CONTINUE 48 49 CALL GETPS(TSTRESS,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 DSTRESS=0.0 59 CALL GETSTRESS(DDSDDE,DSTRESS,DSTRAN,NTENS) 60 DO 702 I1=1,NTENS 61 STRESS(I1)=STRESS(I1)+DSTRESS(I1) 62 702 CONTINUE 63 CALL GETPS(STRESS,PS,NTENS) 64 CALL GETEMOD(PS,EK,EN,RF,C,FAI,ENU,PA,EKUR,EMOD,S,S3O,UG,UD,UF, 65 1 SSS,S1S3O) 66 EBULK3=EMOD/(ONE-TWO*ENU) 67 EG2=EMOD/(ONE+ENU) 68 EG=EG2/TWO 69 EG3=THREE*EG 70 ELAM=(EBULK3-EG2)/THREE 71 CALL GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG) 72 73 74 IF(PS(3).GT.S3O)S3O=PS(3) 75 IF((PS(1)-PS(3)).GT.S1S3O)S1S3O=PS(1)-PS(3) 76 IF(S.GT.SSS)SSS=S 77 STATEV(1)=S1S3O 78 STATEV(2)=S3O 79 STATEV(3)=SSS 80 81 END 82 83 84 85 SUBROUTINE GETPS(STRESS,PS,NTENS) 86 INCLUDE 'ABA_PARAM.INC' 87 DIMENSION PS(3),STRESS(NTENS) 88 CALL SPRINC(STRESS,PS,1,3,3) 89 DO 310 I=1,2 90 DO 320 J=I+1,3 91 IF(PS(I).GT.PS(J))THEN 92 PPS=PS(I) 93 PS(I)=PS(J) 94 PS(J)=PPS 95 END IF 96 320 CONTINUE 97 310 CONTINUE 98 DO 330 K1=1,3 99 PS(K1)=-PS(K1) 100 330 CONTINUE 101 RETURN 102 END 103 104 SUBROUTINE GETEMOD(PS,EK,EN,RF,C,FAI,ENU,PA,EKUR,EMOD,S,S3O 105 1 ,UG,UD,UF,SSS,S1S3O) 106 INCLUDE 'ABA_PARAM.INC' 107 DIMENSION PS(3) 108 S=(1-SIN(FAI))*(PS(1)-PS(3)) 109 IF(PS(3).LT.(-C/TAN(FAI))) THEN 110 S=0.99 111 ELSE 112 S=S/(2*C*COS(FAI)+2*PS(3)*SIN(FAI)) 113 IF(S.GE.0.99)S=0.99 114 END IF 115 EMOD=EK*PA*((S3O/PA)**EN)*((1-RF*S)**2) 116 AA=UD*(PS(1)-PS(3)) 117 AA=AA/(EK*PA*((S3O/PA)**EN)) 118 AA=AA/(1-RF*S) 119 ENU=UG-UF*LOG10(S3O/PA) 120 ENU=ENU/(1-AA)/(1-AA) 121 IF(ENU.GT.0.49)ENU=0.49 122 IF(ENU.LT.0.05)ENU=0.05 123 IF(S.LT.SSS.AND.(PS(1)-PS(3)).LT.S1S3O)THEN 124 EMOD=EKUR*PA*((S3O/PA)**EN) 125 END IF 126 END 127 128 SUBROUTINE GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG) 129 INCLUDE 'ABA_PARAM.INC' 130 DIMENSION DDSDDE(NTENS,NTENS) 131 DO 20 K1=1,NTENS 132 DO 10 K2=1,NTENS 133 DDSDDE(K2,K1)=0.0 134 10 CONTINUE 135 20 CONTINUE 136 DO 40 K1=1,NDI 137 DO 30 K2=1,NDI 138 DDSDDE(K2,K1)=ELAM 139 30 CONTINUE 140 DDSDDE(K1,K1)=EG2+ELAM 141 40 CONTINUE 142 DO 50 K1=NDI+1,NTENS 143 DDSDDE(K1,K1)=EG 144 50 CONTINUE 145 RETURN 146 END 147 148 SUBROUTINE GETSTRESS(DDSDDE,STRESS,DSTRAN,NTENS) 149 INCLUDE 'ABA_PARAM.INC' 150 DIMENSION DDSDDE(NTENS,NTENS),STRESS(NTENS),DSTRAN(NTENS) 151 DO 70 K1=1,NTENS 152 DO 60 K2=1,NTENS 153 STRESS(K1)=STRESS(K1)+DDSDDE(K1,K2)*DSTRAN(K2) 154 60 CONTINUE 155 70 CONTINUE 156 RETURN 157 END
上面的是E-v模型。
abaqus邓肯张EB模型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 19 EK=PROPS(1) 20 EN=PROPS(2) 21 RF=PROPS(3) 22 C=PROPS(4) 23 FAI=PROPS(5)/180.0*3.1415926 24 DFAI=PROPS(6)/180.0*3.1415926 25 EKB=PROPS(7) 26 EM=PROPS(8) 27 EKUR=PROPS(9) 28 PA=PROPS(10) 29 S1S3O=STATEV(1) 30 S3O=STATEV(2) 31 SSS=STATEV(3) 32 33 CALL GETPS(STRESS,PS,NTENS) 34 FAI=FAI-DFAI*LOG10(S3O/PA) 35 CALL GETEMOD(PS,EK,EN,RF,C,FAI,ENU,PA,EKUR,EMOD,S,S3O,EKB,EM 36 1 ,SSS,S1S3O) 37 38 EBULK3=EMOD/(ONE-TWO*ENU) 39 EG2=EMOD/(ONE+ENU) 40 EG=EG2/TWO 41 EG3=THREE*EG 42 ELAM=(EBULK3-EG2)/THREE 43 CALL GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG) 44 DSTRESS=0.0 45 CALL GETSTRESS(DDSDDE,DSTRESS,DSTRAN,NTENS) 46 DO 701 I1=1,NTENS 47 STRESS(I1)=STRESS(I1)+DSTRESS(I1) 48 701 CONTINUE 49 50 CALL GETPS(STRESS,PS,NTENS) 51 52 CALL GETEMOD(PS,EK,EN,RF,C,FAI,ENU,PA,EKUR,EMOD,S,S3O,EKB,EM 53 1 ,SSS,S1S3O) 54 55 56 EBULK3=EMOD/(ONE-TWO*ENU) 57 EG2=EMOD/(ONE+ENU) 58 EG=EG2/TWO 59 EG3=THREE*EG 60 ELAM=(EBULK3-EG2)/THREE 61 62 CALL GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG) 63 64 65 IF(PS(3).GT.S3O)S3O=PS(3) 66 IF((PS(1)-PS(3)).GT.S1S3O)S1S3O=PS(1)-PS(3) 67 IF(S.GT.SSS)SSS=S 68 STATEV(1)=S1S3O 69 STATEV(2)=S3O 70 STATEV(3)=SSS 71 72 END 73 74 75 76 SUBROUTINE GETPS(STRESS,PS,NTENS) 77 INCLUDE 'ABA_PARAM.INC' 78 DIMENSION PS(3),STRESS(NTENS) 79 CALL SPRINC(STRESS,PS,1,3,3) 80 DO 310 I=1,2 81 DO 320 J=I+1,3 82 IF(PS(I).GT.PS(J))THEN 83 PPS=PS(I) 84 PS(I)=PS(J) 85 PS(J)=PPS 86 END IF 87 320 CONTINUE 88 310 CONTINUE 89 DO 330 K1=1,3 90 PS(K1)=-PS(K1) 91 330 CONTINUE 92 RETURN 93 END 94 95 SUBROUTINE GETEMOD(PS,EK,EN,RF,C,FAI,ENU,PA,EKUR,EMOD,S,S3O 96 1 ,EKB,EM,SSS,S1S3O) 97 98 INCLUDE 'ABA_PARAM.INC' 99 DIMENSION PS(3) 100 101 S=(1-SIN(FAI))*(PS(1)-PS(3)) 102 IF(PS(3).LT.(-C/TAN(FAI))) THEN 103 S=0.99 104 ELSE 105 106 S=S/(2*C*COS(FAI)+2*PS(3)*SIN(FAI)) 107 IF(S.GE.0.99)then 108 S=0.99 109 END IF 110 END IF 111 EMOD=EK*PA*((S3O/PA)**EN)*((1-RF*S)**2) 112 IF((S.LT.SSS).AND.((PS(1)-PS(3)).LT.S1S3O))THEN 113 EMOD=EKUR*PA*((S3O/PA)**EN) 114 END IF 115 116 AA=EKB*PA*((S3O/PA)**EM) 117 ENU=0.5*(1.0-EMOD/3./AA) 118 IF(ENU.GT.0.49)ENU=0.49 119 IF(ENU.LT.0.05)ENU=0.05 120 END 121 122 123 SUBROUTINE GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG) 124 INCLUDE 'ABA_PARAM.INC' 125 DIMENSION DDSDDE(NTENS,NTENS) 126 DO 20 K1=1,NTENS 127 DO 10 K2=1,NTENS 128 DDSDDE(K2,K1)=0.0 129 10 CONTINUE 130 20 CONTINUE 131 DO 40 K1=1,NDI 132 DO 30 K2=1,NDI 133 DDSDDE(K2,K1)=ELAM 134 30 CONTINUE 135 DDSDDE(K1,K1)=EG2+ELAM 136 40 CONTINUE 137 DO 50 K1=NDI+1,NTENS 138 DDSDDE(K1,K1)=EG 139 50 CONTINUE 140 RETURN 141 END 142 143 SUBROUTINE GETSTRESS(DDSDDE,STRESS,DSTRAN,NTENS) 144 INCLUDE 'ABA_PARAM.INC' 145 DIMENSION DDSDDE(NTENS,NTENS),STRESS(NTENS),DSTRAN(NTENS) 146 DO 70 K1=1,NTENS 147 DO 60 K2=1,NTENS 148 STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1) 149 60 CONTINUE 150 70 CONTINUE 151 RETURN 152 END 153