第1页 / 共8页
第2页 / 共8页
第3页 / 共8页
第4页 / 共8页
第5页 / 共8页
第6页 / 共8页
第7页 / 共8页
第8页 / 共8页
abaqus burgers模型子程序.doc
---------------------------------------------------------------
C
C
C
瞬态温度场下修正 Burgers 模型 UMAT 子程序源代码
给状态变量数组赋初值为零，调用 ABAQUS 子程序 SDVINI
GIVE STATEV THE INITIAL VALUE OF ZERO
SUBROUTINE SDVINI(STATEV,COORDS,NSTATV,NCRDS,NOEL,NPT,LAYER,KSPT)
INCLUDE 'ABA_PARAM.INC'
DIMENSION STATEV(NSTATV),COORDS(NCRDS)
DO K=1,NSTATV
STATEV(K)=0.0
END DO
RETURN
END
C
C
瞬态温度场下修正 Burgers 模型 UMAT 子程序
UMAT FOR MODIFIED BURGERS MODEL
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)
INCLUDE 'ABA_PARAM.INC'
CHARACTER*80 CMNAME
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)
LOCAL ARRAYS（定义局部数组，即用户自己定义的数组）
C
C
C ----------------------------------------------------------------------------------------------------------------
C
C
C
C
C
C
C
C
C
C
C
C ----------------------------------------------------------------------------------------------------------------
EELAS - ELASTIC STRAINS
DEELA - ELASTIC STRAINS INCREMENT(dt)
ECREE - CREEP STRAINS
DECRE - CREEP STRAINS INCREMENT(dt)
STREST - TEMPORARY ARRAY FOR SAVED STRESS(t+dt)
DECRT - TEMPORARY ARRAY FOR SAVED CREEP STRAINS INCREMENT(dt)
DSTREST- STRESS INCREMENT(dt)
STREST2- TEMPORARY ARRAY FOR SAVED STRESS(t+θdt)
DVSTRESS-DEVIATORIC STRESS(t+θdt)
PARAM - ARRAY FOR SAVED MATERIAL PARAMETERS AT TEMPERATURE
OF THE KSTEP
DIMENSION EELAS(6),ECREE(6),DECRE(6),DVSTRESS(6),STREST(6),
1
DECRT(6),DEELA(6),DSTREST(6),STREST2(6),PARAM(7),
TABLE(NPROPS/7,7),TABLE0(7),TABLE1(7)
2
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)
UMAT FOR MODIFIED BURGERS MODEL
C ----------------------------------------------------------------------------------------------------------------
C
C ----------------------------------------------------------------------------------------------------------------
C
C
C
C
C
C
C
C
C
C ---------------------------------------------------------------------------------------------------------------
C
PARAMETERS OF MODIFIED BURGERS MODEL AT TRANSIENT TEMPERATURE
FIELD（材料参数说明）：
PROPS(1)- T
PROPS(2)- E0
PROPS(3)- μ0
PROPS(4)- E1
PROPS(5)- η1
PROPS(6)- A
PROPS(7)- B
Temperature
Young's modulus of the Hookean spring
Poisson's ratio
Young's modulus of the Kelvin Unit
Viscosity coef. of the Kelvin unit
A and B are parameters of the seperate dashpot
定义初始时刻材料的瞬时弹性参数
NVALUE=NPROPS/7
IF(KINC.EQ.0)THEN
IF(NVALUE.EQ.1)THEN
DO K1=1,7
C
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
参数输入应按温度从小到大输入
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
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
!E0
!μ0
!G
!3G
!λ
EMOD=PARAM(2)
ENU=PARAM(3)
IF(ENU.GT.0.4999.AND.ENU.LT.0.5001)ENU=0.499
EBULK3=EMOD/(ONE-TWO*ENU)
EG2=EMOD/(ONE+ENU)
EG=EG2/TWO
EG3=THREE*EG
ELAM=(EBULK3-EG2)/THREE
迭代之前获取上一步的状态变量
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)
定义瞬时弹性性能
IF(KSTEP.EQ.1)THEN
形成初始时刻弹性刚度
DO K1=1,NTENS
DO K2=1,NTENS
!3K
!2G
!等效蠕变应变
!等效粘性流动应变
!等效粘弹性应变
C
C
C
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
根据弹性应变计算初始时刻应力
DO K1=1,NTENS
DO K2=1,NTENS
C
STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)
END DO
EELAS(K1)=EELAS(K1)+DSTRAN(K1)
END DO
END IF
C------------------------------------------------------------------------------------------------------------------
C
C------------------------------------------------------------------------------------------------------------------
CALCULATE CREEP STRAIN（计算蠕变应变）
C
IF(KSTEP.GT.1)THEN
采用常刚度迭代方案，刚度为初始时刻弹性刚度
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
c
求取 t， t
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)
获取 t 时刻材料参数 E1，η1，A，B
E1=PARAM(4)
Y1=PARAM(5)
A=PARAM(6)
B=PARAM(7)
下面的参数分别表示对外置粘壶和内置粘壶的ε求导
C
C
EQVSER1=SMISES1*EXP(-1.0*B*TIME(1))/A
EQVEER1=SMISES1*EXP(-1.0*E1/Y1*TIME(1))/Y1
EQCRER1=EQVSER1+EQVEER1
迭代之前，认为 (
n
(0)
1)
t
t
(n+1)
t
DO K1=1,NTENS
STREST(K1)=STRESS(K1)
END DO
给存储蠕应变增量的数组赋初值为 0
DO K1=1,NTENS
DECRT(K1)=0.0
END DO
非线性迭代过程，求取△ε(n+1)对应的△εc(n+1)和△σ(n+1)
LOOP1:DO KEWTON=1,NEWTON
t
t
k
求取 (
(
n
t
1)
、
)
(
n
t
c k
(
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)
EQVSER2=SMISES2*EXP(-1.0*B*(TIME(1)+DTIME))/A
EQVEER2=SMISES2*EXP(-1.0*E1/Y1*(TIME(1)+DTIME))/Y1
EQCRER2=EQVSER2+EQVEER2
求取
t
t
(
(
n
k
1)
)
t
t
1)
(
n
(
)
c k
(1
(
)
t
n
1)
t
t
(
(
n
k
1)
)
(1
t
(
n
c
)
1)
t
t
(
1)
n
)
(
c k
SMISES3=0.2*SMISES1+0.8*SMISES2
EQCRER3=0.2*EQCRER1+0.8*EQCRER2
求取
t
(
n
t
)
(
k
1)
t
3
2 t
1)
(
t
n
)
(
c k
(
1)
n
(
)
k
t
TERM1=(THREE/TWO*EQCRER3)/SMISES3
求取
t
t
1)
(
n
)
(
k
(1
t
1)
(
n
)
(
k
(n+1)
)
t
t
DO K1=1,NTENS
STREST2(K1)=0.2*STREST(K1)+0.8*STRESS(K1)
END DO
计算偏应力
t
n
1)
t
(
kS
(
)
PRESS=0.0
DO K1=1,NDI
C
C
C
C
C
C
C
C
C
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
收敛判定
(
n
1)
(
c k
)
er
，则结束迭代
(
n
1)
(
c k
1)
(
n
1)
(
c k
1)
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
1)
(
n
k
1)
(
t
(
1)
n
)
(
k
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
C
C
C
C
C
把蠕变增量保存到临时数组 DECRT 中
DO K1=1,NTENS
DECRT(K1)=DECRE(K1)
END DO
END DO LOOP1
更新应力
DO K1=1,NTENS
STRESS(K1)=STREST(K1)
END DO
更新状态变量，包括蠕变应变、弹性应变以及等效应变等
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
END IF
更新 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
保存状态变量
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
RETURN
END
---------------------------------------------------------------
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
