Johnson-Cook金属塑性本构
ABAQUS/Standard 用户材料子程序实例
－Johnson-Cook 金属本构模型卢剑锋 庄茁* 张帆
清华大学工程力学系 北京 100084
摘要：用户材料子程序是 ABAQUS 提供给用户定义自己的材料属性的 Fortran 程序接口，使用户能使用 ABAQUS 材料库中没有定义的材料模型。
ABAQUS |
中自有的 |
Johnson-Cook |
模型只能应用于显式 |
ABAQUS/Explicit |
程序中，而我们 |
希望能在隐式 |
ABAQUS/Standard |
程序中更精确的实现本构积分，而且应用 |
Johnson-Cook |
模型 |
在 UMAT
编程中使用了率相关塑性理论以及完全隐式的应力更新算法。
1 Johnson-Cook 强化模型简介
Johnson-Cook（JC）模型用来模拟高应变率下的金属材料。JC 强化模型表示为三项的乘积， 分别反映了应变硬化，应变率硬化和温度软化。这里使用 JC 模型的修正形式：
A B n & 1 T *m
1 C ln 1
&
0
并使参考应变率&0 1 ，这样公式中的 A 即为材料的静态屈服应力。公式中包含 A, B, n, C, m 五个参数，需要通过实验来确定。
2 ABAQUS 用户材料子程序
用户材料子程序（User-defined Material Mechanical Behavior，简称 UMAT）通过与ABAQUS 主求解程序的接口实现与 ABAQUS 的数据交流。在输入文件中，使用关键字“*USER MATERIAL”表示定义用户材料属性。
子程序概况与接口
UMAT 子程序具有强大的功能，使用 UMAT 子程序：
(1) 可以定义材料的本构关系，使用 ABAQUS 材料库中没有包含的材料进行计算，扩充程序功能。
(2) 几乎可以用于力学行为分析的任何分析过程，几乎可以把用户材料属性赋予 ABAQUS 中的任何单元；
(3) 必须在 UMAT 中提供材料本构模型的雅可比（Jacobian）矩阵，即应力增量对应变增量的变化率。
(4) 可以和用户子程序“USDFLD”联合使用，通过“USDFLD”重新定义单元每一物质点上传递到 UMAT 中场变量的数值。
由于主程序与 UMAT 之间存在数据传递，甚至共用一些变量，因此必须遵守有关 UMAT 的书写格式，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)
user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDT
RETURN
END
UMAT 中的应力矩阵、应变矩阵以及矩阵 DDSDDE ，DDSDDT ，DRPLDE 等，都是直接分
量存储在前，剪切分量存储在后。直接分量有 NDI 个，剪切分量有 NSHR 个。各分量之间的顺序根据单元自由度的不同有一些差异，所以编写 UMAT 时要考虑到所使用单元的类别。下面对UMAT 中用到的一些变量进行说明：
DDSDDE NTENS, NTENS
是一个 NTENS 维的方阵，称作雅可比矩阵，σ / ε ，σ 是应力的增量，ε 是应变的增量， DDSDDE I , J 表示增量步结束时第 J 个应变分量的改变引起的第 I 个应力分量的变化。通常雅可比是一个对称矩阵，除非在“*USER MATERIAL”语句中加入了“UNSYMM”参数。
STRESS NTENS
应力张量矩阵，对应 NDI 个直接分量和 NSHR 个剪切分量。在增量步的开始，应力张量矩阵中的数值通过 UMAT 和主程序之间的接口传递到 UMAT 中，在增量步的结束 UMAT 将对应力张量矩阵更新。对于包含刚体转动的有限应变问题，一个增量步调用 UMAT 之前就已经对应力张量的进行了刚体转动，因此在 UMAT 中只需处理应力张量的共旋部分。UMAT 中应力张量的度量为柯西（真实）应力。
STATEV NSTATEV
用于存储状态变量的矩阵，在增量步开始时将数值传递到 UMAT 中。也可在子程序 USDFLD 或 UEXPAN 中先更新数据，然后增量步开始时将更新后的数据传递到 UMAT 中。在增量步的结束必须更新状态变量矩阵中的数据。
和应力张量矩阵不同的是：对于有限应变问题，除了材料本构行为引起的数据更新以外，状态变量矩阵中的任何矢量或者张量都必须通过旋转来考虑材料的刚体运动。
状态变量矩阵的维数，等于关键字“*DEPVAR”定义的数值。状态变量矩阵的维数通过ABAQUS 输入文件中的关键字“*DEPVAR”定义，关键字下面数据行的数值即为状态变量矩阵的维数。
材料常数的个数，等于关键字“*USER MATERIAL”中“CONSTANTS”常数设定的值。
PROPS NPROPS
材料常数矩阵，矩阵中元素的数值对应于关键字“*USER MATERIAL”下面的数据行。
SSE ， SPD ， SCD
分别定义每一增量步的弹性应变能，塑性耗散和蠕变耗散。它们对计算结果没有影响，仅仅
作为能量输出。其他变量：
STRAN NTENS ：应变矩阵； DSTRAN NTENS ：应变增量矩阵； DTIME ：增量步的时间增量；
NDI ：直接应力分量的个数；
NSHR ：剪切应力分量的个数；
NTENS ：总应力分量的个数， NTENS NDI NSHR 。
使用 UMAT 时需要注意单元的沙漏控制刚度和横向剪切刚度。通常减缩积分单元的沙漏控制刚度和板、壳、梁单元的横向剪切刚度是通过材料属性中的弹性性质定义的。这些刚度基于材料初始剪切模量的值，通常在材料定义中通过“*ELASTIC”选项定义。但是使用 UMAT 的时候， ABAQUS 对程序输入文件进行预处理的时候得不到剪切模量的数值。所以这时候用户必须使用“*HOURGLASS STIFFNESS” 选项来定义具有沙漏模式的单元的沙漏控制刚度，使用“*TRANSVERSE SHEAR STIFFNESS”选项来定义板、壳、梁单元的横向剪切刚度。
编程
基于上面所述的率相关材料公式和应力更新算法，参照 ABAQUS 用户材料子程序的接口规范，进行 UMAT 的编程。有限元模拟结果将在下一节给出，在最后一节中还给出了相应的程序源代码。
由于 UMAT 在单元的积分点上调用，增量步开始时，主程序路径将通过 UMAT 的接口进入UMAT，单元当前积分点必要变量的初始值将随之传递给 UMAT 的相应变量。在 UMAT 结束时，变量的更新值将通过接口返回主程序。整个 UMAT 的流程如图 1 所示。
一共有 8 个材料常数需要给定，并申请一个 13 维的状态变量矩阵，它们表示的物理含义如表
1 所示。
下一步将使用建立的UMAT 结合ABAQUS/Standard 进行霍布金森冲击杆(SHPB)实验的有限元模拟，并对结果进行比较。
图 1 UMAT 流程图 表 1 UMAT 材料常数
3 SHPB 实验的有限元模拟模型的简化与有限元网格
为了不使模型过于庞大，对模型进行了一些简化。首先，改变入力杆和出力杆的尺寸，长度由原来的 3040mm 减小为 1000mm，直径增加到 25mm，试件的长度和直径也分别变化为 22mm 和 18mm。这样不仅优化了网格的质量，还成倍地减小了模型的规模，其带来的负面影响就是试件能达到的应变将降低。另外，由于撞击杆仅仅起到产生应力脉冲的作用，在数值模型中没必要考虑撞击杆，取代的方法是直接在入力杆的输入端施加均布的应力脉冲。
考虑到实验装置的对称性，也做了一些简化。整个实验装置以及载荷等都是关于杆的中心线轴对称的，所以可以使用轴对称单元进行二维分析。
二维轴对称模型如图 2 所示。在模型中，对试件以及入力杆，出力杆和试件接触的部分进行了局部网格加密，这样的网格划分可以取得比较经济的结果。
图 2二维轴对称有限元模型
表 2 模型信息
模型 |
尺寸[mm] ( × L ) |
单元类型 |
单元个数 |
总节点 数 |
总单元数 |
|
二维模型 |
入力杆 |
25×1000 |
CAX4 |
530 |
1475 |
1220 |
试件 |
18×22 |
CAX4 |
160 |
|||
出力杆 |
25×1000 |
CAX4 |
530 |
|||
材料定义
入力杆和出力杆使用线弹性材料，弹性模量和泊松比分别为 200GPa 和 0.3，密度为 7.85×103 kg/m3。试件采用用户在 UMAT 中自定义材料，材料参数如表 3 所示，其中 Johnson-Cook 模型中参数的数值来源于前面的数值拟合程序。
表 3 试件的材料定义
性质 |
密度 [Kg/m3 ] |
杨氏模量 [MPa] |
泊松比 |
Johnson-Cook 模型参数 |
||||
A [MPa ] |
B [MPa] |
n |
C |
M |
||||
数 值 |
2.7×103 |
68.0×103 |
0.33 |
66.56 2 |
108.85 3 |
0.23 8 |
0.029 |
0.5 |
边界条件
为了保证 SHPB 实验的要求，在二维模型中施加了必要的边界条件。在对称轴上施加了对称性边界条件，同时保证压杆和试件可以沿轴线方向自由无约束的运动。压杆和试件之间的接触为硬接触，光滑无摩擦。
为了确定输入应力脉冲的时间，进行了简单的计算。弹性材料中纵波波速的计算公式为：
Cd
其中 E 为材料弹性模量， 为材料密度。由此可以计算输入应力波在压杆中的传播速度为
Cd 5048 m/s。
要求在入力杆应力波的输入端不能出现入射波和反射波的重叠，也就是说在输入应力脉冲的时间内，应力波的传播距离不应超过两倍的杆长，即：
T 2L
Cd
2
5048
4.0 104 (s)
根据这一估计，选择输入应力脉冲的持续时间Ts
2.0 104 s，上升时间t
3.0 105 s。
经过若干次试算，对输入应力脉冲的波形进行了适当的调整，使试件中产生较均匀的应变率。最后输入应力脉冲的波形如图 3 所示：
300
250
200
150
100
50
0
0.0000
0.0001
0.00020.00030.00040.0005
时 间 (s)
图 3输入应力脉冲
为了确定增量步的最大时间步长，需要先简单计算一下单元的稳定极限。基于一个单元的估
算，稳定极限可以用单元特征长度 Le 和材料波速Cd 定义如下：
Le
tstable
d
压杆单元的特征单元长度 Le 10 mm，由此可以计算出应力波在压杆传递的稳定极限为
t 2.0 106 (s)
将它作为 ABAQUS 自动增量控制里面的最大时间步长。二维动态分析
我们所对照的 SHPB 实验正是属于这一情况，所以可以将 ABAQUS/Standard 结合 UMAT
进行有限元模拟的结果和实验数据进行对比。下面是应变率 250 s-1 下的动态模拟过程。
在时间t 1.98104 (s) 左右，应力波前沿到达试件，这一时间和前面使用弹性波波速计算的传播时间是相同的，此前试件上的 Mises 应力几乎为零，如图 4 所示。
图 4应力波前沿到达试件时的 Mises 应力 (t=1.98×10-4 s)
在时间t 3.0 104 (s) ，试件经过应力波的上升时间后达到稳定变形的状态，一部分入射波
反射回入力杆，一部分应力波经过试件进入出力杆，试件各点的变形都很均匀，如图 5 (a)所示。在图 5 (b)试件的放大图上可以看出，各点 Mises 应力相差不超过 1MPa，这个精度是相当可靠的。
(a) 全局视图
(b) 试件的放大视图
图 5 试件经历均匀变形时的 Mises 应力 (t=3.0×10-4 s)
经过稳定变形阶段后，反射波和传递波分别向入力杆和出力杆扩散，试件上 Mises 应力逐渐减小到较低的水平，试件开始经历卸载，如图 6 所示。图中 Mises 应力云图的单位为 KPa。
图 6应力波消退后试件时的 Mises 应力 (t=4.2×10-4 s)
实际上有限元模拟的应力－应变曲线和恒定应变率下实验的结果也能够很好的吻合。取出试件表面中间的一点，将应变率 250 s-1 和 200 s-1 下 ABAQUS 有限元模拟的结果与实验的结果对比
见图 7 和图 8。
160
140
120
100
500
450
400
350
300
80250
60200
150
40
100
2050
00
0.0000.0050.0100.0150.0200.0250.0300.0350.0400.045
应 变
图 7 应力－应变曲线的对比及模拟过程中真实应变率变化 (250 s-1)
160
140
120
100
500
450
400
350
300
80250
60200
150
40
100
2050
00
0.0000.0050.0100.0150.0200.0250.0300.035
应 变
图 8 应力－应变曲线的对比及模拟过程中真实应变率变化 (200 s-1) 4 J-C UMAT 程序
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)
C
INCLUDE 'ABA_PARAM.INC'
C
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)
C
DIMENSION EELAS(6),EPLAS(6),FLOW(6)
PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0, HALF =0.5d0) DATA NEWTON,TOLER/40,1.D-6/
C
C CUMAT FOR JOHNSON-COOK MODEL
C CPROPS(1) - YANG'S MODULUS
CPROPS(2) - POISSON RATIO
CPROPS(3) - INELASTIC HEAT FRACTION
CPARAMETERS OF JOHNSON-COOK MODEL: CPROPS(4) - A
CPROPS(5) - B
CPROPS(6) - n
CPROPS(7) - C
CPROPS(8) - m
C C
IF (NDI.NE.3) THEN WRITE(6,1)
1FORMAT(//,30X,'***ERROR - THIS UMAT MAY ONLY BE USED FOR ',
1'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS') ENDIF
C
CELASTIC PROPERTIES C
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
CELASTIC STIFFNESS C
DO 20 K1=1,NTENS
DO 10 K2=1,NTENS DDSDDE(K2,K1)=0.0
10CONTINUE
20CONTINUE C
DO 40 K1=1,NDI
DO 30 K2=1,NDI DDSDDE(K2,K1)=ELAM
30CONTINUE DDSDDE(K1,K1)=EG2+ELAM
40CONTINUE
DO 50 K1=NDI+1,NTENS
DDSDDE(K1,K1)=EG
50CONTINUE C
CCALCULATE STRESS FROM ELASTIC STRAINS C
DO 70 K1=1,NTENS
DO 60 K2=1,NTENS STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)
60CONTINUE
70CONTINUE C
CRECOVER ELASTIC AND PLASTIC STRAINS C
DO 80 K1=1,NTENS EELAS(K1)=STATEV(K1)+DSTRAN(K1) EPLAS(K1)=STATEV(K1+NTENS)
80CONTINUE
EQPLAS=STATEV(1+2*NTENS)
C
CCALCULATE MISES STRESS C
IF(NPROPS.GT.5.AND.PROPS(4).GT.0.0) THEN 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)
90CONTINUE SMISES=SQRT(SMISES/TWO)
C
CCALL USERHARD SUBROUTINE, GET HARDENING RATE AND YIELD STRESS C
C
CALL USERHARD(SYIEL0,HARD,EQPLAS,PROPS(4)) CDETERMINE IF ACTIVELY YIELDING
C
IF (SMISES.GT.(1.0+TOLER)*SYIEL0) THEN
C
CMATERIAL RESPONSE IS PLASTIC, DETERMINE FLOW DIRECTION C
SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE ONESY=ONE/SMISES
DO 110 K1=1,NDI FLOW(K1)=ONESY*(STRESS(K1)-SHYDRO)
110CONTINUE
DO 120 K1=NDI+1,NTENS FLOW(K1)=STRESS(K1)*ONESY
120CONTINUE C
CREAD PARAMETERS OF JOHNSON-COOK MODEL C
A=PROPS(4) B=PROPS(5) EN=PROPS(6) C=PROPS(7) EM=PROPS(8)
C
CNEWTON ITERATION C
SYIELD=SYIEL0
DEQPL=(SMISES-SYIELD)/EG3 DSTRES=TOLER*SYIEL0/EG3 DEQMIN=HALF*DTIME*EXP(1.0D-4/C) DO 130 KEWTON=1,NEWTON
DEQPL=MAX(DEQPL,DEQMIN)
CALL USERHARD(SYIELD,HARD,EQPLAS+DEQPL,PROPS(4)) TVP=C*LOG(DEQPL/DTIME)
TVP1=TVP+ONE HARD1=HARD*TVP1+SYIELD*C/DEQPL SYIELD=SYIELD*TVP1
RHS=SMISES-EG3*DEQPL-SYIELD DEQPL=DEQPL+RHS/(EG3+HARD1)
IF(ABS(RHS/EG3) .LE. DSTRES ) GOTO 140
130CONTINUE WRITE(6,2) NEWTON
2FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT ',
1'CONVERGE AFTER ',I3,' ITERATIONS')
140CONTINUE EFFHRD=EG3*HARD1/(EG3+HARD1)
C
CCALCULATE STRESS AND UPDATE STRAINS 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
150CONTINUE
DO 160 K1=NDI+1,NTENS STRESS(K1)=FLOW(K1)*SYIELD EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL
EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL
160CONTINUE EQPLAS=EQPLAS+DEQPL SPD=DEQPL*(SYIEL0+SYIELD)/TWO RPL = PROPS(3)*SPD/DTIME
C
CJACOBIAN 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
230CONTINUE
DO 250 K1=1,NTENS
DO 240 K2=1,NTENS DDSDDE(K2,K1)=DDSDDE(K2,K1)+FLOW(K2)*FLOW(K1)
1*(EFFHRD-EFFG3)
240CONTINUE
250CONTINUE ENDIF
ENDIF
C
CSTORE STRAINS IN STATE VARIABLE ARRAY C
DO 310 K1=1,NTENS STATEV(K1)=EELAS(K1) STATEV(K1+NTENS)=EPLAS(K1)
310CONTINUE
STATEV(1+2*NTENS)=EQPLAS
C
RETURN END
C C C
SUBROUTINE USERHARD(SYIELD,HARD,EQPLAS,TABLE)
C
INCLUDE 'ABA_PARAM.INC'
C
DIMENSION TABLE(3)
C
CGET PARAMETERS, SET HARDENING TO ZERO C
A=TABLE(1) B=TABLE(2) EN=TABLE(3) HARD=0.0
C
CCALSULATE CURRENT YIELD STRESS AND HARDENING RATE C
IF(EQPLAS.EQ.0.0) THEN SYIELD=A
ELSE
HARD=EN*B*EQPLAS**(EN-1) SYIELD=A+B*EQPLAS**EN
END IF
RETURN END
