mises本构模型UMAT（附源代码和详细注释）

墨鬼

2019年4月8日 22:12

浏览：7086 评论：19 收藏：52

一. ABAQUS自带mises本构与UMAT对比

二. 如何在分析中使用自定义材料

1.通过ABAQUS/CAE使用自定义材料

在Property模块中，执行【Material】/【Create】命令（或单击相关工具箱区的按钮），弹出Edit Material对话框（如图所示），用户可以通过该对话框选择材料模型、设置材料参数。对于自定义模型，执行对话框的【General】/【User Material】命令，此时在Material Behaviors区域中会出现User Material字样，表明定义的是用户材料。在Edit Material对话框下方的User Material 区域中的【User material type】下拉列表中有三个选项，分别为Mechanical（力学材料）、 Thermal (热学材料）和Thermomechanical (热力学模型），默认选项为Mechanical。对于一些岩土材料模型，尤其是采用非关联流动法则的模型，雅克比矩阵是不对称的，此时需勾选Use unsymmetric material stiffiiess matrix 复选框。材料的参数在 Data 区域中的【Mechanical Constants】列表中输入，这里的数据会按次序传给UMAT子程序中的PROPS (NPROPS)数组，数据的个数即为NPROPS。

注：输入材料参数时，在数据列表尾部按回车健会自动增加一行数据.此外，在数据窗口上右击会显示弹出式菜单，可实现数据的拷贝、增加、删除等功能。

如果UMAT子程序用到状态变量,还需设置状态变量的个数。具体操作仍然在Edit Material对话框中进行，执行对话框的【General】/【Depvar】命令，在【Number of solution-dependent state variables】输入框中设罝状态变量的个数。

2.通过inp输入文件使用自定义材料

在inp输入文件中的材料定义选项块（* Material)中使用下列语句即可：

*User material,type=Mechanical,constants= number of constants,Unsymm

关键字User material表明定义的用户材料;type为材料类型，可为Mechanical(力学材料）、 Thermal (热学材料）和Thermomechanical (热力学模型），Mechanical是默认选项：Constants 选项给出的是材料参数的个数；Unsymm选项只有当雅克比矩阵非对称时才采用。该关键字行下的数据行中的数据为材料参数。

状态变量的个数同样在材料定义选项块中定义，相应的关键字行为：

*Depvar

Number of solution-dependent state variables：状态变量个数

三. Mises弹塑性模型

4.源程序和注释

需要启用沙漏控制：

mesh沙漏.jpg

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

1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,

2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,

3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,

4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)

INCLUDE 'ABA_PARAM.INC'

CHARACTER*80 MATERL

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),

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

DIMENSION EELAS(6),EPLAS(6),FLOW(6)

C EELAS(1-6)为弹性应变分量,EPLAS(1-6)为塑性应变分量,FLOW(6)为屈服面对应力的导数

PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)

DATA NEWTON,TOLER/10,1.D-6/

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

C UMAT FOR ISOTROPIC ELASTICITY AND ISOTROPIC PLASTICITY

C J2 FLOW THEORY

C CAN NOT BE USED FOR PLANE STRESS

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

C PROPS(1) - E ,为弹性模量

C PROPS(2) - NU ,为泊松比

C PROPS(3) - SYIELD, 为屈服应力

C CALLS AHARD FOR CURVE OF SYIELD VS. PEEQ

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

C 程序不能用于平面应变

IF (NDI.NE.3) THEN

WRITE(6,1)

1 FORMAT(//,30X,'***ERROR - THIS UMAT MAY ONLY BE USED FOR ',

1 'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS')

ENDIF

C ELASTIC PROPERTIES

EMOD=PROPS(1)

ENU=PROPS(2)

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

C 得到拉梅常数λ(ELAM),剪切模量G（EG）

C ELASTIC STIFFNESS

DO 20 K1=1,NTENS

DO 10 K2=1,NTENS

DDSDDE(K2,K1)=0.0

10 CONTINUE

20 CONTINUE

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 得到（6-1）中的弹性模量矩阵

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 RECOVER ELASTIC AND PLASTIC STRAINS

DO 80 K1=1,NTENS

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

C STATEV(1-6)存储弹性应变分量，此时假定应变增量全为弹性

EPLAS(K1)=STATEV(K1+NTENS)

C STATEV(6-12)存储塑性应变分量

80 CONTINUE

EQPLAS=STATEV(1+2*NTENS)

C STATEV(13)为等效塑性应变

C 如果没有给出屈服应力，材料认为是弹性的

IF(NPROPS.GT.2.AND.PROPS(3).GT.0.0) THEN

C MISES STRESS

SMISES=(STRESS(1)-STRESS(2))*(STRESS(1)-STRESS(2)) +

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

1 (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 计算mises应力

C HARDENING CURVE, GET YIELD STRESS

C 等效塑性应变和硬化曲线确定硬化后的屈服应力,初始等效塑性应变为0,

C 对应的屈服应力为初始屈服应力

NVALUE=NPROPS/2-1

CALL AHARD(SYIEL0,HARD,EQPLAS,PROPS(3),NVALUE)

C DETERMINE IF ACTIVELY YIELDING

IF (SMISES.GT.(1.0+TOLER)*SYIEL0) THEN

C 如果Mises应力超出了屈服应力,判断发生了屈服,进入了塑性阶段

C FLOW DIRECTION

SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE

C SHYDRO为平均应力

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 FLOW(1-6)为 s_ij/q

C SOLVE FOR EQUIV STRESS, NEWTON ITERATION

SYIELD=SYIEL0

DEQPL=0.0

DO 130 KEWTON=1,NEWTON

RHS=SMISES-EG3*DEQPL-SYIELD

C 计算屈服面函数误差

DEQPL=DEQPL+RHS/(EG3+HARD)

C 按(6-13)计算新的等效塑性应变

CALL AHARD(SYIELD,HARD,EQPLAS+DEQPL,PROPS(3),NVALUE)

C 由新的等效塑性应变获得新的屈服应力

IF(ABS(RHS).LT.TOLER*SYIEL0) GOTO 140

C 如果满足收敛要求结束,否则再次迭代

130 CONTINUE

WRITE(6,2) NEWTON

2 FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT ',

1 'CONVERGE AFTER ',I3,' ITERATIONS')

140 CONTINUE

EFFHRD=EG3*HARD/(EG3+HARD)

C CALC STRESS AND UPDATE STRAINS

C 计算及更新应力,弹性应变,塑性应变

DO 150 K1=1,NDI

C 更新正应力分量

STRESS(K1)=FLOW(K1)*SYIELD+SHYDRO

C 按(6-16)计算应力

EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO

C 更新塑性应变

EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO

C 更新弹性应变

150 CONTINUE

DO 160 K1=NDI+1,NTENS

C 更新剪应力分量

STRESS(K1)=FLOW(K1)*SYIELD

EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL

EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL

160 CONTINUE

EQPLAS=EQPLAS+DEQPL

SPD=DEQPL*(SYIEL0+SYIELD)/TWO

C JACOBIAN

C 计算雅克比矩阵

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)

1 *(EFFHRD-EFFG3)

240 CONTINUE

250 CONTINUE

ENDIF

C STORE STRAINS IN STATE VARIABLE ARRAY

DO 310 K1=1,NTENS

STATEV(K1)=EELAS(K1)

STATEV(K1+NTENS)=EPLAS(K1)

310 CONTINUE

STATEV(1+2*NTENS)=EQPLAS

RETURN

END

SUBROUTINE AHARD(SYIELD,HARD,EQPLAS,TABLE,NVALUE)

C 根据导入的硬化曲线及当前的等效塑性应变,求得屈服应力及屈服应力对等效塑性应变的偏导

INCLUDE 'ABA_PARAM.INC'

DIMENSION TABLE(2,NVALUE)

C SET YIELD STRESS TO LAST VALUE OF TABLE, HARDENING TO ZERO

SYIELD=TABLE(1,NVALUE)

HARD=0.0

C IF MORE THAN ONE ENTRY, SEARCH TABLE

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 ',

1 'ENTERED IN ASCENDING ORDER')

CALL XIT

ENDIF

C CURRENT YIELD STRESS AND HARDENING

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

源程序附件：

mises本构模型UMAT（附源代码和详细注释）的图16 mises.zip

登录后免费查看全文

立即登录

技术邻APP
工程师必备

项目客服
培训客服
平台客服

TOP

mises本构模型UMAT（附源代码和详细注释）

推荐阅读