1. 程式人生 > >abaqus鄧肯張模型umat

abaqus鄧肯張模型umat

pen har -s return call drp nstat isp cte

首先是始點剛度法:

技術分享圖片
  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始點剛度法

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
abaqus鄧肯張umat中點剛度法

上面的是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       
abaqus鄧肯張EB模型UMAT

abaqus鄧肯張模型umat