+用户材料子程序实例-Johnson-Cook+金属本构模型

Cd =
E
ρ
其中 E 为材料弹性模量， ρ 为材料密度。由此可以计算输入应力波在压杆中的传播速度为 Cd = 5048 m/s。 要求在入力杆应力波的输入端不能出现入射波和反射波的重叠，也就是说在输入应力脉冲的 时间内，应力波的传播距离不应超过两倍的杆长，即：
Ts <
2L 2 = ≈ 4.0 ×10−4 ( s ) Cd 5048
-5-
1 杨氏模量 1-6 弹性应变
2 泊松比
3 塑性耗散比 7-12 塑性应变
4 A
5 B
6 n 13
7 C
8 M
等效塑性应变
Back
ABAQUS 软件 2003 年度用户年会论文集
图2
二维轴对称有限元模型 表 2 模型信息
模 二维 模 型 性 质 数 值
型 入力杆 试件 出力杆
尺寸[mm] (Φ ×L) 25×1000 18×22 25×100件 2003 年度用户年会论文集
(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)
σ = ( A + Bε n ) 1 + C ln 1 +
& ε *m (1 − T ) &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 材料库中没有包含的材料进行计算，扩充程序 功能。
-4-
Back
ABAQUS 软件 2003 年度用户年会论文集
图 1 UMAT 流程图 表 1 UMAT 材料常数 PROPS 物理性质 STATEV 变量意义 3 SHPB 实验的有限元模拟 模型的简化与有限元网格 为了不使模型过于庞大，对模型进行了一些简化。首先，改变入力杆和出力杆的尺寸，长度 由原来的 3040mm 减小为 1000mm，直径增加到 25mm，试件的长度和直径也分别变化为 22mm 和 18mm。这样不仅优化了网格的质量，还成倍地减小了模型的规模，其带来的负面影响就是试 件能达到的应变将降低。另外，由于撞击杆仅仅起到产生应力脉冲的作用，在数值模型中没必要 考虑撞击杆，取代的方法是直接在入力杆的输入端施加均布的应力脉冲。 考虑到实验装置的对称性，也做了一些简化。整个实验装置以及载荷等都是关于杆的中心线 轴对称的，所以可以使用轴对称单元进行二维分析。 二维轴对称模型如图 2 所示。在模型中，对试件以及入力杆，出力杆和试件接触的部分进行 了局部网格加密，这样的网格划分可以取得比较经济的结果。
单元类型 CAX4 CAX4 CAX4
单元个数 530 160 530
总节点 数
总单元数
1475
1220
材料定义 入力杆和出力杆使用线弹性材料，弹性模量和泊松比分别为 200GPa 和 0.3，密度为 7.85×103 kg/m3。试件采用用户在 UMAT 中自定义材料，材料参数如表 3 所示，其中 Johnson-Cook 模型 中参数的数值来源于前面的数值拟合程序。 表 3 试件的材料定义 密度 [Kg/m3 ] 杨氏模 量 [MPa] 泊 松 比 Johnson-Cook 模型参数 A [MPa ] 66.56 2 B [MPa] 108.85 3
是一个 NTENS 维的方阵， 称作雅可比矩阵，∂∆σ / ∂∆ε ，∆σ 是应力的增量，∆ε 是应变的增量，
DDSDDE ( I , J ) 表示增量步结束时第 J 个应变分量的改变引起的第 I 个应力分量的变化。 通常雅可
比是一个对称矩阵，除非在“*USER MATERIAL”语句中加入了“UNSYMM”参数。
STRESS ( NTENS )
应力张量矩阵，对应 NDI 个直接分量和 NSHR 个剪切分量。在增量步的开始，应力张量矩阵 中的数值通过 UMAT 和主程序之间的接口传递到 UMAT 中，在增量步的结束 UMAT 将对应力 张量矩阵更新。对于包含刚体转动的有限应变问题，一个增量步调用 UMAT 之前就已经对应力 张量的进行了刚体转动，因此在 UMAT 中只需处理应力张量的共旋部分。UMAT 中应力张量的 度量为柯西（真实）应力。
Back
ABAQUS 软件 2003 年度用户年会论文集
ABAQUS/Standard 用户材料子程序实例
－Johnson-Cook 金属本构模型 卢剑锋 庄茁* 张帆 清华大学工程力学系 北京 100084
摘要：用户材料子程序是 ABAQUS 提供给用户定义自己的材料属性的 Fortran 程序接口，使用 户能使用 ABAQUS 材料库中没有定义的材料模型。 ABAQUS 中自有的 Johnson-Cook 模型只能应用于显式 ABAQUS/Explicit 程序中，而我们 希望能在隐式 ABAQUS/Standard 程序中更精确的实现本构积分，而且应用 Johnson-Cook 模型 的修正形式。这就需要通过 ABAQUS/Standard 的用户材料子程序 UMAT 编程实现。在 UMAT 编程中使用了率相关塑性理论以及完全隐式的应力更新算法。 1 Johnson-Cook 强化模型简介 Johnson-Cook（JC）模型用来模拟高应变率下的金属材料。JC 强化模型表示为三项的乘积， 分别反映了应变硬化，应变率硬化和温度软化。这里使用 JC 模型的修正形式：
STATEV ( NSTATEV )
用于存储状态变量的矩阵， 在增量步开始时将数值传递到 UMAT 中。 也可在子程序 USDFLD 或 UEXPAN 中先更新数据，然后增量步开始时将更新后的数据传递到 UMAT 中。在增量步的结 束必须更新状态变量矩阵中的数据。 和应力张量矩阵不同的是：对于有限应变问题，除了材料本构行为引起的数据更新以外，状 态变量矩阵中的任何矢量或者张量都必须通过旋转来考虑材料的刚体运动。 状态变量矩阵的维数，等于关键字 “*DEPVAR” 定义的数值。状态变量矩阵的维数通过 ABAQUS 输入文件中的关键字“*DEPVAR”定义，关键字下面数据行的数值即为状态变量矩阵的 维数。 材料常数的个数，等于关键字“*USER MATERIAL”中“CONSTANTS”常数设定的值。
根据这一估计，选择输入应力脉冲的持续时间 Ts = 2.0 × 10−4 s，上升时间 tr = 3.0 ×10−5 s。 经过若干次试算， 对输入应力脉冲的波形进行了适当的调整， 使试件中产生较均匀的应变率。 最后输入应力脉冲的波形如图 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 等，都是直接分
-2-
Back
ABAQUS 软件 2003 年度用户年会论文集
量存储在前，剪切分量存储在后。直接分量有 NDI 个，剪切分量有 NSHR 个。各分量之间的顺序 根据单元自由度的不同有一些差异，所以编写 UMAT 时要考虑到所使用单元的类别。下面对 UMAT 中用到的一些变量进行说明：
DDSDDE ( NTENS , NTENS )
-6-
n
C
M
2.7×103
68.0×103
0.33
0.23 8
0.029
0.5
Back
边界条件
ABAQUS 软件 2003 年度用户年会论文集
为了保证 SHPB 实验的要求，在二维模型中施加了必要的边界条件。在对称轴上施加了对称 性边界条件，同时保证压杆和试件可以沿轴线方向自由无约束的运动。压杆和试件之间的接触为 硬接触，光滑无摩擦。 为了确定输入应力脉冲的时间，进行了简单的计算。弹性材料中纵波波速的计算公式为：