- 积分
- 5
- 注册时间
- 2009-10-23
- 仿真币
-
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2010-1-13 23:18:32
|
显示全部楼层
来自 黑龙江哈尔滨
本帖最后由 david840530 于 2010-1-16 20:28 编辑
新发现的另外一个版本的UMAT---源程序+注解![比我那个好多了 ]
******************************主程序******************************
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,KSTEP,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)
C
DIMENSION EELAS(6),EPLAS(6),FLOW(6),DVECT(4)
REAL DAMAGE,DAMAGC,DAMAG0,BITA,EPDC,EPD0,EQPLAF,DDAMAG,FFDET
PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)
DATA NEWTON,TOLER/10,1.D-6/
C
IF(NDI.NE.3) THEN
WRITE(6,1)
1 FORMAT(//,30X,'***ERROR-THIS UMAT MAY ONLY BE USED FOR','ELEMENTS
& WITH THREE DIRECT STRESS COMPONENTS')
ENDIF
WRITE(6,*) NTENS
C
C
C ELASTIC PROPERTIES
C 根据弹性模量,泊松比计算体积模量,剪切模量,拉梅常数
EMOD=PROPS(1)
ENU=PROPS(2)
c//避免体积模量无穷大
IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499
c//3*体积模量
EBULK3=EMOD/(ONE-TWO*ENU)
EG2=EMOD/(ONE+ENU)
EG=EG2/TWO
EG3=THREE*EG
ELAM=(EBULK3-EG2)/THREE
C
C ELASTIC STIFFNESS
C 弹性刚度,先置零,后赋值
DO 20 K1=1,NTENS
DO 10 K2=1,NTENS
DDSDDE(K2,K1)=0.0
10 CONTINUE
20 CONTINUE
C
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
C
C CALCULATE STRESS FROM ELASTIC STRAINS
C 弹性应变引起的应力
DO 70 K1=1,NTENS
DO 60 K2=1,NTENS
STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)
60 CONTINUE
70 CONTINUE
C
C RECVER ELASTIC AND PLASTIC STRAINS
C 更新弹性和塑性应变,说明STATEV(1~6)是6个方向弹性应变,7~12是塑性应变
C 等效应变EQPLAS是STATEV(13)
DO 80 K1=1,NTENS
EELAS(K1)=STATEV(K1)+DSTRAN(K1)
EPLAS(K1)=STATEV(K1+NTENS)
80 CONTINUE
EQPLAS=STATEV(1+2*NTENS)
WRITE(6,1000) KINC, NOEL, NPT, EQPLAS
1000 FORMAT(//,30X,"INCREMENT",2X,I2,2X,'ELEMENT NUMBER',
& 2X,I3,2X,'NPT',2X,I1,2X,'CURRENT EQ PLASTIC STRAIN',2X,F12.6)
C
IF(NPROPS.GT.2.AND.PROPS(3).GT.0.0) THEN
SMISES=(STRESS(1)-STRESS(2))*(STRESS(1)-STRESS(2))+
&(STRESS(2)-STRESS(3))*(STRESS(2)-STRESS(3))+
&(STRESS(3)-STRESS(1))*(STRESS(3)-STRESS(1))
DO 90 K1=NDI+1,NTENS
SMISES=SMISES+SIX*STRESS(K1)*STRESS(K1)
90 CONTINUE
SMISES=SQRT(SMISES/TWO)
C 调用子程序得到屈服应力
NVALUE=NPROPS/2-1
CALL AHARD(SYIEL0,HARD,EQPLAS,PROPS(3),NVALUE)
C 判断是否屈服
IF(SMISES.GT.(1.0+TOLER)*SYIEL0)THEN
C 屈服后的流动方向
SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE
ONESY=ONE/SMISES
DO 110 K1=1,NDI
FLOW(K1)=ONESY*(STRESS(K1)-SHYDRO)
110 CONTINUE
DO 120 K1=NDI+1,NTENS
FLOW(K1)=STRESS(K1)*ONESY
120 CONTINUE
C 牛顿迭代求等效塑性应变
C
SYIELD=SYIEL0
DEQPL=0.0
DO 130 KEWTON=1,NEWTON
RHS=SMISES-EG3*DEQPL-SYIELD
DEQPL=DEQPL+RHS/(EG3+HARD)
CALL AHARD(SYIELD,HARD,EQPLAS+DEQPL,PROPS(3),NVALUE)
IF (ABS(RHS).LT.TOLER*SYIEL0) THEN
C Output incremental equivalent plastic strain
WRITE(6,2000) KINC,NOEL,NPT,DEQPL
2000 FORMAT(//,30X,"INCREMENT",2X,I2,2X,"ELEMENT NUMBER",2X,I3,2X,
& 'NPT',2X,I1,2X,"EQ PLASTIC STRAIN INC",2X,F12.6)
GOTO 140
ENDIF
130 CONTINUE
WRITE(6,2) NEWTON
2 FORMAT(//,30X,'***WARNING-PLASTICITY ALGORITHM DID NOT',
&'CONVERGE AFTER',I3,'ITERATIONS')
140 CONTINUE
EFFHRD=EG3*HARD/(EG3+HARD)
C
C计算应力更新应变
DO 150 K1=1,NDI
STRESS(K1)=FLOW(K1)*SYIELD+SHYDRO
EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO
EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO
150 CONTINUE
EQPLAS=EQPLAS+DEQPL
SPD=DEQPL*(SYIEL0+SYIELD)/TWO
C更新J矩阵
EFFG=EG*SYIELD/SMISES
EFFG2=TWO*EFFG
EFFG3=THREE*EFFG2/TWO
EFFLAM=(EBULK3-EFFG2)/THREE
DO 220 K1=1,NDI
DO 210 K2=1,NDI
DDSDDE(K2,K1)=EFFLAM
210 CONTINUE
DDSDDE(K1,K1)=EFFG2+EFFLAM
220 CONTINUE
DO 230 K1=NDI+1,NTENS
DDSDDE(K1,K1)=EFFG
230 CONTINUE
DO 250 K1=1,NTENS
DO 240 K2=1,NTENS
DDSDDE(K2,K1)=DDSDDE(K2,K1)+FLOW(K2)*FLOW(K1)
& *(EFFHRD-EFFG3)
240 CONTINUE
250 CONTINUE
ENDIF
ENDIF
C储存状态变量
DO 310 K1=1,NTENS
STATEV(K1)=EELAS(K1)
STATEV(K1+NTENS)=EPLAS(K1)
310 CONTINUE
STATEV(1+2*NTENS)=EQPLAS
C print equivalent plastic strain to .dat file
WRITE(6,3000) KINC, NOEL,NPT, STATEV(1+2*NTENS)
3000 FORMAT(//,30X,"INCREMENT",2X,I2,2X,'ELEMENT NUMBER',
& 2X,I3,2X,'NPT',2X,I1,2X,'TOTAL EQ PLASTIC STRAIN',2X,F12.6)
C
RETURN
END
C
C***子程序****ahard***************
SUBROUTINE AHARD(SYIELD,HARD,EQPLAS,TABLE,NVALUE)
INCLUDE 'ABA_PARAM.INC'
DIMENSION TABLE(2,NVALUE)
C 求EQPLAS在哪段斜率内,然后线性叠加求应力,返回SYIELD和斜率HARD
SYIELD=TABLE(1,NVALUE)
HARD=0.0
C
C IF MORE THAN ONE ENTRY,SEARCH TABLE
C
IF(NVALUE.GT.1) THEN
DO 10 K1=1,NVALUE-1
EQPL1=TABLE(2,K1+1)
IF(EQPLAS.LT.EQPL1) THEN
EQPL0=TABLE(2,K1)
IF(EQPL1.LE.EQPL0) THEN
WRITE(6,1)
1 FORMAT(//,30X,'***ERROR-PLASTIC STRAIN MUST BE','ENTERED IN ASCEND
&ING ORDER')
CALL XIT
ENDIF
C
C CURRENT YIELD STRESS AND HARDENING
C
DEQPL=EQPL1-EQPL0
SYIEL0=TABLE(1,K1)
SYIEL1=TABLE(1,K1+1)
DSYIEL=SYIEL1-SYIEL0
HARD=DSYIEL/DEQPL
SYIELD=SYIEL0+(EQPLAS-EQPL0)*HARD
GOTO 20
ENDIF
10 CONTINUE
20 CONTINUE
ENDIF
RETURN
END
|
|