ABAQUS路面材料使用修正burgers模型时总是出现编译错误

C 瞬态温度场下修正Burgers模型UMAT子程序源代码

C

C 给状态变量数组赋初值为零，调用ABAQUS子程序SDVINI

C GIVE STATEV THE INITIAL VALUE OF ZERO

C

SUBROUTINE SDVINI(STATEV,COORDS,NSTATV,NCRDS,NOEL,NPT,LAYER,KSPT)

C

INCLUDE 'ABA_PARAM.INC'

C

DIMENSION STATEV(NSTATV),COORDS(NCRDS)

C

DO K=1,NSTATV

STATEV(K)=0.0

END DO

C

RETURN

END

C 瞬态温度场下修正Burgers模型UMAT子程序

C UMAT FOR MODIFIED BURGERS MODEL

C

SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,

1 DRPLDE,DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,

2 CMNAME,NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,

3 PNEWDT,CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)

C

INCLUDE 'ABA_PARAM.INC'

C

CHARACTER*80 CMNAME

C

DIMENSION STRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS),

1 DDSDDT(NTENS),DRPLDE(NTENS),STRAN(NTENS),DSTRAN(NTENS),

2 TIME(2),PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3),DROT(3,3),

3 DFGRD0(3,3),DFGRD1(3,3)

C

C LOCAL ARRAYS（定义局部数组，即用户自己定义的数组）

C ----------------------------------------------------------------------------------------------------------------

C EELAS - ELASTIC STRAINS

C DEELA - ELASTIC STRAINS INCREMENT(dt)

C ECREE - CREEP STRAINS

C DECRE - CREEP STRAINS INCREMENT(dt)

C STREST - TEMPORARY ARRAY FOR SAVED STRESS(t+dt)

C DECRT - TEMPORARY ARRAY FOR SAVED CREEP STRAINS INCREMENT(dt)

C DSTREST- STRESS INCREMENT(dt)

C STREST2- TEMPORARY ARRAY FOR SAVED STRESS(t+θdt)

C DVSTRESS-DEVIATORIC STRESS(t+θdt)

C PARAM - ARRAY FOR SAVED MATERIAL PARAMETERS AT TEMPERATURE

C OF THE KSTEP

C ----------------------------------------------------------------------------------------------------------------

C

DIMENSION EELAS(6),ECREE(6),DECRE(6),DVSTRESS(6),STREST(6),

1 DECRT(6),DEELA(6),DSTREST(6),STREST2(6),PARAM(7),

2 TABLE(NPROPS/7,7),TABLE0(7),TABLE1(7)

C

PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,SIX=6.D0,

1 ENUMAX=.4999D0,NEWTON=50,TOLER=1.0D-4)

C

C ----------------------------------------------------------------------------------------------------------------

C UMAT FOR MODIFIED BURGERS MODEL

C ----------------------------------------------------------------------------------------------------------------

C PARAMETERS OF MODIFIED BURGERS MODEL AT TRANSIENT TEMPERATURE

C FIELD（材料参数说明）：

C PROPS(1)- T Temperature

C PROPS(2)- E0 Young's modulus of the Hookean spring

C PROPS(3)- μ0 Poisson's ratio

C PROPS(4)- E1 Young's modulus of the Kelvin Unit

C PROPS(5)- η1 Viscosity coef. of the Kelvin unit

C PROPS(6)- A

C PROPS(7)- B A and B are parameters of the seperate dashpot

C ---------------------------------------------------------------------------------------------------------------

C 定义初始时刻材料的瞬时弹性参数

NVALUE=NPROPS/7

IF(KINC.EQ.0)THEN

     IF(NVALUE.EQ.1)THEN

     DO K1=1,7

PARAM(K1)=PROPS(K1)

END DO

     END IF

     IF(NVALUE.GT.1)THEN

     DO K1=1,NVALUE

DO K2=1,7

TABLE(K1,K2)=PROPS(7*(K1-1)+K2)

END DO

END DO

C 参数输入应按温度从小到大输入

LOOP2:DO K1=1,NVALUE-1

TEMP1=TABLE(K1+1,1)

IF(TEMP.LT.TEMP1)THEN

TEMP0=TABLE(K1,1)

DO K2=1,7

TABLE1(K2)=TABLE(K1+1,K2)

TABLE0(K2)=TABLE(K1,K2)

END DO

DO K3=1,7

C 对不同温度时的参数进行插值

PARAM(K3)=TABLE0(K3)+(TEMP-TEMP0)

1 *(TABLE1(K3)-TABLE0(K3))/(TEMP1-TEMP0)

END DO

EXIT LOOP2

END IF

END DO LOOP2

IF(TEMP.LT.TABLE(1,1).OR.TEMP.GT.TABLE(NVALUE,1))THEN

DO K1=1,7

TABLE1(K1)=TABLE(NVALUE,K1)

TABLE0(K1)=TABLE(1,K1)

END DO

C 对不同温度时的参数进行插值

DO K2=1,7

PARAM(K2)=TABLE0(K2)+(TEMP-TABLE(1,1))*(TABLE1(K2)

1 -TABLE0(K2))/(TABLE(NVALUE,1)-TABLE(1,1))

END DO

END IF

     END IF

     END IF

EMOD=PARAM(2) !E0

ENU=PARAM(3) !μ0

IF(ENU.GT.0.4999.AND.ENU.LT.0.5001)ENU=0.499

EBULK3=EMOD/(ONE-TWO*ENU) !3K

EG2=EMOD/(ONE+ENU) !2G

EG=EG2/TWO !G

EG3=THREE*EG !3G

ELAM=(EBULK3-EG2)/THREE !λ

C 迭代之前获取上一步的状态变量

DO K1=1,NTENS

EELAS(K1)=STATEV(K1) !弹性应变

ECREE(K1)=STATEV(K1+NTENS) !蠕变应变

END DO

EQCRE=STATEV(1+2*NTENS) !等效蠕变应变

EQVSE=STATEV(2+2*NTENS) !等效粘性流动应变

EQVEE=STATEV(3+2*NTENS) !等效粘弹性应变

C 定义瞬时弹性性能

IF(KSTEP.EQ.1)THEN

C 形成初始时刻弹性刚度

DO K1=1,NTENS

DO K2=1,NTENS

DDSDDE(K1,K2)=ZERO

END DO

END DO

DO K1=1,NDI

DO K2=1,NDI

DDSDDE(K2,K1)=ELAM

END DO

DDSDDE(K1,K1)=EG2+ELAM

END DO

DO K1=NDI+1,NTENS

DDSDDE(K1,K1)=EG

END DO

C 根据弹性应变计算初始时刻应力

DO K1=1,NTENS

DO K2=1,NTENS

STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)

END DO

EELAS(K1)=EELAS(K1)+DSTRAN(K1)

END DO

C

END IF

C

C------------------------------------------------------------------------------------------------------------------

C CALCULATE CREEP STRAIN（计算蠕变应变）

C------------------------------------------------------------------------------------------------------------------

IF(KSTEP.GT.1)THEN

C 采用常刚度迭代方案，刚度为初始时刻弹性刚度

DO K1=1,NTENS

DO K2=1,NTENS

DDSDDE(K1,K2)=ZERO

END DO

END DO

DO K1=1,NDI

DO K2=1,NDI

DDSDDE(K2,K1)=ELAM

END DO

DDSDDE(K1,K1)=EG2+ELAM

END DO

DO K1=NDI+1,NTENS

DDSDDE(K1,K1)=EG

END DO

C 求取t σ avg，t εc avg

SMISES1=(STRESS(1)-STRESS(2))**2+(STRESS(2)-STRESS(3))**2

1 +(STRESS(3)-STRESS(1))**2

DO K1=NDI+1,NTENS

SMISES1=SMISES1+SIX*STRESS(K1)**2

END DO

SMISES1=SQRT(SMISES1/TWO)

     write(6,*)SMISES1

C 获取t时刻材料参数E1，η1，A，B

E1=PARAM(4)

Y1=PARAM(5)

A=PARAM(6)

B=PARAM(7)

C 下面的参数分别表示对外置粘壶和内置粘壶的ε求导

EQVSER1=SMISES1*EXP(-1.0*B*TIME(1))/A

EQVEER1=SMISES1*EXP(-1.0*E1/Y1*TIME(1))/Y1

EQCRER1=EQVSER1+EQVEER1

C 迭代之前，认为t+Δt σ(0)(n+1)= t σ(n+1)

     DO K1=1,NTENS

     STREST(K1)=STRESS(K1)

     END DO

C 给存储蠕应变增量的数组赋初值为0

     DO K1=1,NTENS

     DECRT(K1)=0.0

     END DO

C 非线性迭代过程，求取△ε(n+1)对应的△εc(n+1)和△σ(n+1)

LOOP1:DO KEWTON=1,NEWTON

C 求取t+△t σk(n+1)，t+△t εck(n+1)

SMISES2=(STREST(1)-STREST(2))**2+(STREST(2)-STREST(3))**2

1 +(STREST(3)-STREST(1))**2

DO K1=NDI+1,NTENS

SMISES2=SMISES2+SIX*STREST(K1)**2

END DO

SMISES2=SQRT(SMISES2/TWO)

     write(6,*)SMISES2

EQVSER2=SMISES2*EXP(-1.0*B*(TIME(1)+DTIME))/A

EQVEER2=SMISES2*EXP(-1.0*E1/Y1*(TIME(1)+DTIME))/Y1

EQCRER2=EQVSER2+EQVEER2

C 求取t+θ△t σavg (k)(n+1)=(1-θ) t σavg(n+1)+θ t+△t σavg(k)(n+1)

C t+θ△t εcavg (k)(n+1)=(1-θ) t εcavg(n+1)+θ t+△t εcavg(k)(n+1)

SMISES3=0.2*SMISES1+0.8*SMISES2

     write(6,*)SMISES3

     EQCRER3=0.2*EQCRER1+0.8*EQCRER2

C 求取t+θ△t β(k)(n+1)=3/2* t+θ△t εcavg (k)(n+1) / t+θ△t σavg (k)(n+1)

TERM1=(THREE/TWO*EQCRER3)/SMISES3

C 求取t+θ△t σ(k)(n+1)=(1-θ) t σ(n+1)+θ t+△t σ(k)(n+1)

DO K1=1,NTENS

STREST2(K1)=0.2*STREST(K1)+0.8*STRESS(K1)

END DO

C 计算偏应力t+θ△t S(k)(n+1)

PRESS=0.0

DO K1=1,NDI

PRESS=PRESS+STREST2(K1)/THREE

END DO

DO K1=1,NDI

DVSTRESS(K1)=STREST2(K1)-PRESS

END DO

DO K2=NDI+1,NTENS

DVSTRESS(K2)=STREST2(K2)

END DO

C 计算蠕应变、弹性应变、等效蠕应变等的增量

DO K1=1,NTENS

DECRE(K1)=TERM1*DTIME*DVSTRESS(K1)

DEELA(K1)=DSTRAN(K1)-DECRE(K1)

END DO

DEQCRE=EQCRER3*DTIME

DEQVSE=(0.2*EQVSER1+0.8*EQVSER2)*DTIME

DEQVEE=(0.2*EQVEER1+0.8*EQVEER2)*DTIME

C 收敛判定abs(△εc(k+1)(n+1)-△εc(k)(n+1))/abs(△εc(k+1)(n+1))<=toler，则结束迭代

TERM2=0.0

TERM3=0.0

DO K1=1,NTENS

TERM2=TERM2+(DECRE(K1)-DECRT(K1))**2

TERM3=TERM3+DECRE(K1)**2

END DO

TERM4=SQRT(TERM2)

TERM5=SQRT(TERM3)

TERM6=TERM4/TERM5

IF(TERM6.LE.TOLER)THEN

EXIT LOOP1

END IF

C 求取t+△t σ(k+1)(n+1)=t σ+△σ(k)(n+1)

DO K1=1,NTENS

DSTREST(K1)=0.0

END DO

DO K1=1,NTENS

DO K2=1,NTENS

DSTREST(K1)=DSTREST(K1)+DDSDDE(K1,K2)*DEELA(K2)

END DO

END DO

DO K1=1,NTENS

STREST(K1)=STRESS(K1)+DSTREST(K1)

END DO

C 把蠕变增量保存到临时数组DECRT中

     DO K1=1,NTENS

     DECRT(K1)=DECRE(K1)

     END DO

END DO LOOP1

C 更新应力

DO K1=1,NTENS

     STRESS(K1)=STREST(K1)

     END DO

C 更新状态变量，包括蠕变应变、弹性应变以及等效应变等

DO K1=1,NTENS

ECREE(K1)=ECREE(K1)+DECRE(K1)

EELAS(K1)=EELAS(K1)+DEELA(K1)

END DO

EQCRE=EQCRE+DEQCRE

EQVSE=EQVSE+DEQVSE

EQVEE=EQVEE+DEQVEE

C

END IF

C 更新JACOBIAN矩阵

DO K1=1,NTENS

DO K2=1,NTENS

DDSDDE(K1,K2)=ZERO

END DO

END DO

DO K1=1,NDI

DO K2=1,NDI

DDSDDE(K2,K1)=ELAM

END DO

DDSDDE(K1,K1)=EG2+ELAM

END DO

DO K1=NDI+1,NTENS

DDSDDE(K1,K1)=EG

END DO

C 保存状态变量

DO K1=1,NTENS

STATEV(K1)=EELAS(K1)

STATEV(K1+NTENS)=ECREE(K1)

END DO

STATEV(1+2*NTENS)=EQCRE

STATEV(2+2*NTENS)=EQVSE

STATEV(3+2*NTENS)=EQVEE

C

RETURN

END