abaqus1用户材料子程序

19 ABAQUS用户材料子程序(UMAT)

虽然ABAQUS为用户提供了大量的单元库和求解模型，使用户能够利用这些模型处理绝大多数的问题；但是现实世界毕竟十分复杂，ABAQUS不可能把所有可能出现的问题都包含进去。所以ABAQUS提供了大量的用户子程序（User Subroutine）。用户子程序允许用户在找不到合适模型的情况下自行定义符合自己问题的模型。这些用户子程序涵盖了建模从载荷到单元的几乎各个部分。

ABAQUS为用户提供的这个接口，允许用户通过自定义的子程序定制ABAQUS，以实现特定的功能。用户子程序具有以下的功能和特点：（1）如果ABAQUS的一些固有选项模型功能有限；用户子程序可以提高ABAQUS中这些选项的功能；（2）通常用户子程序是用FORTRAN语言的代码写成；（3）它可以以几种不同的方式包含在模型中；（4）由于它们没有存储在restart文件中，如果需要的话，可以在重新开始运行时修改它；（5）在某些情况下它可以利用ABAQUS允许的已有程序。

要在模型中包含用户子程序，可以利用ABAQUS执行程序，在abaqus执行程序中应用user选项指明包含这些子程序的FORTRAN源程序或者目标程序的名字。 提示：ABAQUS的输入文件除了可以通过ABAQUS/CAE的作业模块中提交运行外，还可以在ABAQUS Command窗口中输入ABAQUS执行程序直接运行： ABAQUS job=输入文件名 user=用户子程序的Fortran文件名

ABAQUS/Standard和ABAQUS/Explicit都支持用户子程序功能，但是他们所支持的用户子程序种类不尽相同，读者在需要使用时请注意查询手册。

在接下来的最后两章里，我们将讨论两种常用的用户子程序——用户材料子程序和用户单元子程序。

本章将通过在ABAQUS/Standard中创建Johnson-Cook的材料模型，对编写Standard的用户材料子程序UMAT进行一个简单介绍。ABAQUS/Explicit中的用户材料子程序VUMAT的思想与之相似，但是由于隐式和显式两种方法本身的差异，它们之间也有一些不同，请读者在自己具体使用前首先仔细查阅ABAQUS手册中的相关内容。

19-1

18.1 引言

用户材料子程序是ABAQUS提供给用户自的定义材料属性的Fortran程序接口，它使用户能使用ABAQUS材料库中没有的材料模型。ABAQUS中自有的Johnson-Cook模型只能应用于显式ABAQUS/Explicit程序中，而我们希望能在隐式ABAQUS/Standard程序中更精确的实现本构积分，而且应用Johnson-Cook模型的修正形式。这就需要通过ABAQUS/Standard的用户材料子程序UMAT编程实现。在UMAT编程中使用了率相关塑性理论以及完全隐式的应力更新算法。

18.2 模型的数学描述

18.2.1 Johnson-Cook强化模型简介

Johnson-Cook模型用来模拟在冲击载荷作用下的变形。

Johnson-Cook强化模型（JC）表示为三项的乘积，分别反映了应变硬化、应变率硬化和温度软化。这里使用JC模型的修正形式：

???A?B?n??1?Cln?1??????????*m???1?T? ?0????*(19-1)

?0?1，这样模型中包含A,B,n,C,m五个参数，需要通过实验来确定。使参考应变率?公式中的A即为材料的静态屈服应力。公式中的T为无量纲化的温度

T*?T?Tr Tm?Tr其中Tr为室温，Tm为材料的熔点。Johnson-Cook模型在温度从室温到材料熔点温度的范围内都是有效的。

高应变率的变形经常伴有温升现象，这是因为材料变形过程中塑性功转化为热量。对于大多数金属，90-100％的塑性变形将耗散为热量。所以JC模型中温度的变化可以

19-2

用如下的公式计算：

?T??????d? ?c?(19-2)

其中，?T为温度的增量，?为塑性耗散比，表示塑性功转化为热量的比例，c为材料的比热，?为材料的密度。公式（19-2）考虑的是一个绝热过程，即认为温度的升高完全起因于塑性耗散。

18.2.2 率相关塑性的基本公式

Johnson-Cook本构模型考虑率相关塑性，塑性变形是关联的，即塑性流动沿着屈服面的法线方向，并采用Mises屈服面。

将应变的增量分解为弹性部分和塑性部分：

d??d??d?

将上式两端同时对时间的增量dt微分得到率形式：

ep(19-3)

(19-4)

在率相关塑性中，材料的塑性反应取决于加载率，以率的形式给出材料的弹性反应：

????e???p ?

?e?E(????p) ??E??(19-5)

为了发生塑性变形，率相关塑性必须满足或者超过屈服条件，塑性应变率为：

??????，??????，???? ，?(19-6)

???为塑性率参数，?为运动硬化时的背上式即为流动法则，其中?为塑性流动势能，?

??p???应力。

对于各项同性硬化，不存在背应力，因此有??0，此时有：

?????sign???

??p?????????sign??? ????????(19-7)

(19-8)

与率无关塑性不同的是，率相关塑性中等效塑性应变率不能通过一致性条件获得，而是直接通过经验定律给出，成为过应力模型：

19-3

????(?,?,?) ? (19-9)

式中?是过应力，?为粘性。在过应力模型中，等效塑性应变率取决于超过了多少屈服应力。

上面一维的率相关塑性公式可以很方便地推广到三维情况。对于小应变的情形，应力度量之间无需区分，这里采用Cauchy应力σ，塑性率参数由应力和内变量的经验函数给出。对照一维情况，三维情况下分解应变率为弹性和塑性部分：

??ε??ε? ε

应力率和弹性应变率之间的关系为：

ep(19-10)

??C:ε?e?C:?ε??ε?p? σ(19-11)

塑性流动法则和内变量的演化方程为：

塑性率参数为：

?h ?r?σ,q?，q????p??ε(19-12)

?????σ,q? ?n(19-13)

对于J2流动理论，Perzyna（1971）中提出了典型的过应力模型为：

???Y?????Y????1

(19-14)

式中

为Macualay括号，如果f?0，则f?f；如果f?0，则f?0。?为

Mises等效应力，?为等效应变，n为率敏感系数。

对于Johnson-Cook模型，可以得到等效塑性应变率的表达式为：

?1???????exp???1???1 ????C??0??????(19-15)

其中?0为静态屈服应力：

19-4

?0?A?B?n

(19-16)

18.2.3 完全隐式的应力更新算法

对率形式的本构方程进行积分的算法称为应力更新算法。在完全隐式的算法中，在步骤结束时计算塑性应变和内变量的增量，同时强化屈服条件。积分算法写为：

εn?1?εn??ε

pp εn?ε?1n???n?1rn?1ppqn?1?qn???n?1hn?1 pσn?1?C:?εn?1?εn?1?

(19-17a)

(19-17b) (19-17c) (19-17d) (19-17e)

fn?1?f?σn?1,qn?1??0

p在时刻n给出一组εn,εn,qn和应变增量?ε，公式（19-17）是一组关于求解

???εn?1p,εn?1,qn?1?的非线性代数方程。

将公式（19-17b）代入（19-17d）得到：

ppσn?1?C:?εn?1?εn??εn?1?

pppp?C:?εn??ε?εn??εn?1??C:?εn?εn??C:?ε?C:?εn?1??σn?C:?ε??C:?ε (19-18)

pn?1ptrail?σtrailn?1?C:?εn?1?σn?1???n?1C:rn?1式中σtrailn?1?σn?C:?ε是弹性预测的试应力，而数值???n?1C:rn?1是塑性修正量，它沿着结束点塑性流动的方向。对于J2流动理论，塑性流动的方向为：

19-5

3σdevr?

2?(19-19)

它是屈服面的法向，即r?fσ。在偏应力空间，Von Mises屈服面为环状，所以屈服面的法向通过圆心，如图19-1所示：

图19-1 径向返回算法

从图19-1可以看出，在弹性预测阶段，塑性应变和内变量保持固定；而在塑性修正阶段，总体应变保持不变。

在迭代开始时，程序对应力和应变设初始值：

p?0?p?0??0??0??0?p，???n，???0，σ?C:εn?1?ε ?εnε??(19-20)

应力在第k次迭代时为：

σ?k??σ?0?????k?C:r?k?

^^(19-21)

定义屈服面的单位法向矢量为：

n?r/r?0??0??σ0dev/σ0dev，r?0??3/2n

(19-22)

且在整个算法的塑性修正状态过程中始终保持不变，因此塑性应变的更新是??的线性函数。

在k次迭代时将检查屈服条件：

19-6

f?k????k???Y??k????0??3????k???Y??k?

????????3??????????

?????0kkkY??????(19-23)

若收敛，则迭代完毕，增量步结束。否则将计算塑性参数的增量：

(19-24)

3??H?k?并对塑性应变和内变量进一步更新：

n?σ^0dev/σ0dev，?εp?k??????k?32n，???k?????k? (19-25a)

(19-25b)

?k?^εp?k?1??εp?k???εp?k?

σ?k?1??C:εn?1?ε?p?k?1???σ?k???σ?k??σ?2????k?32n (19-25c)

(19-25d) (19-25e)

^??k?1????k?????k?

???k?1?????k?????k?

然后将更新的变量返回屈服条件进行检查，整个过程将重复直至收敛为止。这就是增量步中应力更新的过程。

18.3 ABAQUS用户材料子程序

用户材料子程序（User-defined Material Mechanical Behavior，简称UMAT）通过与ABAQUS主求解程序的接口实现与ABAQUS的数据交流。在输入文件中，使用关键字“*USER MATERIAL”表示定义用户材料属性。

18.3.1 子程序概况与接口

UMAT子程序具有强大的功能，使用UMAT子程序：

(1) 可以定义材料的本构关系，使用ABAQUS材料库中没有包含的材料进行计

算，扩充程序功能。

19-7

(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中用到的一些变量进行说明：

19-8

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中。在增量步结束时必须更新状态变量矩阵中的数据。

和应力张量矩阵不同的是：对于有限应变问题，除了材料本构行为引起的数据更新以外，状态变量矩阵中的任何矢量或者张量都必须通过旋转来考虑材料的刚体运动。

NSTATEV

状态变量矩阵的维数，等于关键字“*DEPVAR”中定义的数值。状态变量矩阵的维数通过ABAQUS输入文件中的关键字“*DEPVAR”定义，关键字下面数据行的数值即为状态变量矩阵的维数。

NPROPS

材料常数的个数，等于关键字“*USER MATERIAL”中“CONSTANTS”常数设定的值。

19-9

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”选项来定义板、壳、梁单元的横向剪切刚度。

18.3.2 编程

基于上面所述的率相关材料公式和应力更新算法，参照ABAQUS用户材料子程序的接口规范，进行UMAT的编程。有限元模拟结果将在下一节给出，在最后一节中还给出了相应的程序源代码。

由于UMAT在单元的积分点上调用，增量步开始时，主程序路径将通过UMAT的

19-10

接口进入UMAT，单元当前积分点必要变量的初始值将随之传递给UMAT的相应变量。在UMAT结束时，变量的更新值将通过接口返回主程序。整个UMAT的流程如图19-2所示：

图19-2 UMAT流程图

一共有8个材料常数需要给定，并申请了一个13维的状态变量矩阵，它们表示的物理含义如表19-1所示：

表19-1 UMAT材料常数

PROPS 物理性质 STATEV 变量意义 1 杨氏模量 1-6 弹性应变 2 泊松比 3 塑性耗散比 7-12 塑性应变 4 A 5 B 6 n 13 等效塑性应变 7 C 8 M

19-11

下一步将使用建立的UMAT结合ABAQUS/Standard进行SHPB实验的有限元模拟，并对结果进行比较。

18.4 SHPB实验的有限元模拟

下面将建立SHPB实验的有限元模型，并把前面所建立的UMAT接入ABAQUS/Standard进行有限元模拟。进行有限元模拟的目的不是单纯为了再现SHPB实验的过程，同时也是为了对选择的数值模型和建立的UMAT进行评价。

18.4.1 分离式Hopkinson压杆（SHPB）实验

分离式Hopkinson压杆（Split Hopkinson Pressure Bar，简称SHPB）实验是从经典Hopkinson实验基础之上发展而来的一种实验技术，用来测量材料的动态应力－应变行为。该实验技术的理论基础是一维应力波理论，通过测量两根压杆上的应变来推导试件上的应力－应变关系。

分离式Hopkinson压杆实验示意图见图19-3。

图19-3 分离式Hopkinson压杆装置

18.4.2 有限元建模

有限元模型主要是参照前面介绍的SHPB实验装置，通过载荷、边界条件等的定义在有限元中模拟SHPB实验的环境，尽量在较少的机时耗费下达到更高的精度。

19-12

模型的简化与有限元网格

为了不使模型过于庞大，对模型进行了一些简化。首先，改变入力杆和出力杆的尺寸，长度由原来的3040mm减小为1000mm，直径增加到25mm，试件的长度和直径也分别变化为22mm和18mm。这样不仅优化了网格的质量，还成倍地减小了模型的规模，其带来的负面影响就是试件能达到的应变将降低。另外，由于撞击杆仅仅起到产生应力脉冲的作用，在数值模型中没必要考虑撞击杆，取代的方法是直接在入力杆的输入端施加均布的应力脉冲。

考虑到实验装置的对称性，也做了一些简化。整个实验装置以及载荷等都是关于杆的中心线轴对称的，所以可以使用轴对称单元进行二维分析。另外也建立了四分之一横截面的三维模型作为补充。

二维轴对称模型和三维模型分别如图19-4、图19-5所示。在模型中，对试件以及入力杆，出力杆和试件接触的部分进行了局部网格加密，这样的网格划分可以取得比较经济的结果。

图19-4 二维轴对称有限元模型

19-13

图19-5 三维有限元模型

单元类型上，选择一阶常规单元，由于没有使用减缩积分单元，所以使用UMAT时无需指定单元的沙漏控制刚度。最后的模型中，二维网格单元总数为1220，三维模型网格的单元总数为17160。关于单元的详细信息如表19-2所示。

表19-2 模型信息

模 型 入力杆 试件 出力杆 入力杆 试件 出力杆 尺寸[mm] (?×L) 25×1000 18×22 25×1000 25×1000 18×22 25×1000 单元类型 CAX4 CAX4 CAX4 C3D8 C3D8 C3D8 单元个数 530 160 530 7500 2160 7500 21049 17160 1475 1220 总节点数 总单元数 二维模型 三维模型

19-14

在二维模型局部网格存在疏密连接的部位，一个单元边要同时和两个单元连接，这在通常的有限元网格中是不好实现的。本例在这里使用了ABAQUS的多点约束技术（Multi-Points Constrain，简称MPC）来解决，线性多点约束方式如图19-6所示：

图19-6 线性多点约束技术

图中p点使用线性多点约束后，其节点自由度均由旁边的节点a和b线性插值得到。所以使用多点约束方式可以很好连接模型中网格疏密不同的部位，划分出比较精练的网格。 材料定义

入力杆和出力杆使用线弹性材料，弹性模量和泊松比分别为200GPa和0.3，密度为7.85×103 kg/m3。试件采用用户在UMAT中的自定义材料，材料参数如表19-3所示，其中Johnson-Cook模型中参数的数值来源于对实验数据的拟合。

表19-3 试件的材料定义

密度 [Kg/m3] 2.7×103 杨氏模量 [MPa] 68.0×103 泊松比 0.33 Johnson-Cook模型参数 A [MPa] 66.562 B [MPa] 108.853 n 0.238 C 0.029 M 0.5 性质 数值

19-15

边界条件

为了保证SHPB实验的要求，在二维模型和三维模型中均施加了必要的边界条件。在对称轴或对称面上施加了对称性边界条件，同时保证压杆和试件可以沿轴线方向自由无约束的运动。压杆和试件之间的接触为硬接触，光滑无摩擦。

为了确定输入应力脉冲的时间，进行了简单的计算。弹性材料中纵波波速的计算公式为：

Cd?E? (19-26)

其中E为材料弹性模量，?为材料密度。由此可以计算输入应力波在压杆中的传播速度为Cd?5048m/s。

要求在入力杆应力波的输入端不能出现入射波和反射波的重叠，也就是说在输入应力脉冲的时间内，应力波的传播距离不应超过两倍的杆长，即：

Ts?2L2??4.0?10?4(s) Cd5048(19-27)

根据这一估计，选择输入应力脉冲的持续时间Ts?2.0?10?4s，上升时间

tr?3.0?10?5s。

经过若干次试算，对输入应力脉冲的波形进行了适当的调整，使试件中产生较均匀的应变率。最后输入应力脉冲的波形如图19-7所示：

为了确定增量步的最大时间步长，需要先简单计算一下单元的稳定极限。基于一个单元的估算，稳定极限可以用单元特征长度L和材料波速Cd定义如下：

e

?tstableLe ?Cd(18-28)

19-16

300 基 准压 力250200(1.7×10,200)(3.0×10,170)-5-4应 力 (MPa)150应 变率(s-1) 系数10050 250 0.90 200 0.77 100 0.52 70 0.4600.00000.00010.00020.00030.00040.0005时 间 (s)图19-7 输入应力脉冲

压杆单元的特征单元长度Le = 10 mm，由此可以计算出应力波在压杆传递的稳定极限为

?tstable?2.0?10?6 (s)

(19-29)

将它作为ABAQUS自动增量控制里面的最大时间步长。

18.4.3 二维动态分析

下面将讨论二维动态分析的结果。为了便于比较，进行了三大类的分析，首先是无温度影响时的强化模型，然后是考虑温度影响的强化模型，最后是考虑单元失效的模型。 无温度影响强化模型分析

所进行的SHPB实验正是属于这一情况，所以可以将ABAQUS/Standard结合UMAT进行有限元模拟的结果和实验数据进行对比。

下面是应变率250 s-1下的动态模拟过程。

在时间t?1.98?10?4(s)左右，应力波前沿到达试件，这一时间和前面使用弹性波波

19-17

速计算的传播时间是相同的，此前试件上的Mises应力几乎为零，如图19-8所示。

图19-8 应力波前沿到达试件时的Mises应力 (t=1.98×10-4 s)

在时间t?3.0?10?4(s)，试件经过应力波的上升时间后达到稳定变形的状态，一部分入射波反射回入力杆，一部分应力波经过试件进入出力杆，试件各点的变形都很均匀，如图19-9 (a) 所示。在图19-9 (b) 试件的放大图上可以看出，各点Mises应力相差不超过1MPa，这个精度是相当可靠的。

(a) 全局视图

(b) 试件的放大视图

图19-9 试件经历均匀变形时的Mises应力 (t=3.0×10-4 s)

19-18

经过稳定变形阶段后，反射波和传递波分别向入力杆和出力杆扩散，试件上Mises应力逐渐减小到较低的水平，试件开始经历卸载，如图19-10所示。图中Mises应力云图的单位为KPa，如不作特别说明，下面各种应力云图中应力单位都为KPa。

图19-10 应力波消退后试件时的Mises应力 (t=4.2×10-4 s)

在离压杆两端0.2m处各取一个单元作为输出，其应力历史如图19-11所示，从中可直观地看出压杆上应力的传播过程。

15010050应力 S22 单元021 单元081 单元709 单元769应 力 (MPa)0-50-100-150-2000.00000.00010.00020.00030.00040.0005时 间 (s)图19-11 压杆上的应力输出（实际输出）

取出试件同一横截面的三个单元以及试件表面长度方向的三个单元，将不同点的应力应变历史比较如下：

19-19

0.040.03应变E22 单元540 单元620 单元680应 变0.020.010.000.00000.00010.00020.00030.00040.0005时间 (s)(a) 应 变历 史140120100

应力 S22 单元540 单元620 单元580应 力 (MPa)8060402000.00000.00010.00020.00030.00040.0005时 间 (s)(b) 应 力 历 史300

200应 变率 SR22 单元540 单元620 单元580应 变率 (s-1)1000-1000.00000.00010.00020.00030.00040.0005时 间 (s)(c) 应 变率 历 史

图19-12 试件横截面应变、应力及应变率历史比较

19-20

0.050.04应变 E22 单元671 单元680 单元6900.03应 变0.020.010.000.00000.00010.00020.00030.00040.0005时 间 (s)(a) 应 变历 史

160140120应 力 (MPa)10080604020应力 S22 单元671 单元680 单元69000.00000.00010.00020.00030.00040.0005时 间 (s)(b) 应 力 历 史

400300应 变率 SR22 单元671 单元680 单元690应 变率 (s-1)2001000-1000.00000.00010.00020.00030.00040.0005时 间 (s)(c) 应 变率 历 史

图19-13 试件长度方向应变、应力及应变率历史比较

19-21

从试件各点的应力应变分布上看，图19-12中应变、应力及应变率历史曲线基本重合，同一横截面内各点的变化历史基本一致。图19-13中试件长度方向各点的变形历史也基本一致，除了初始阶段和最后卸载阶段曲线出现较小波动外，其他阶段曲线也基本重合。出现小幅度波动的原因可能是应力波刚刚到达试件不能立即达到平衡状态，另外应力波经历试件的弹塑性变形后沿其他方向也会出现扩散。

可以认为总体上试件的变形是均匀的，试件各点的应变率基本保持在250 s-1。其他三种应变率下（分别为70、100、200 s-1）的有限元模拟结果与250 s-1的情况大体相同，试件的变形均匀，各点的应变率基本保持在设计的水平。这给研究恒定应变率下的应力－应变曲线提供了有力的前提条件。

实际上有限元模拟的应力－应变曲线和恒定应变率下实验的结果也能够很好的吻合。取出试件表面中间的一点，将应变率250 s-1和200 s-1下ABAQUS有限元模拟的结果与实验的结果对比见图19-14和图19-15。

四种应变率（分别为70、100、200、250 s-1）下使用ABAQUS进行有限元模拟结果的对比见图19-16到图19-20：

160140120350300250200150401002000.0005000.045500450400 实验值 ABAQUS模拟 应 变率80600.0050.0100.0150.0200.0250.0300.0350.040应 变图19-14 应力－应变曲线的对比及模拟过程中真实应变率变化 (250 s-1)

应 变率 (s-1)应 力 (MPa)100

19-22

160140120 实验数据 ABAQUS模拟500 应 变率450400350应 力 (MPa)3002502001508060401002000.0005000.0350.0050.0100.0150.0200.0250.030应 变图19-15 应力－应变曲线的对比及模拟过程中真实应变率变化 (200 s-1)

0.050.04 应 变 应 变 应 变 应 变率 70率 100率 200率 2500.03应 变0.020.010.000.00000.00010.00020.00030.00040.0005时 间 (s)图19-16 四种应变率下的应变－时间曲线

19-23

应 变率 (s-1)100

140120100 应 变 应 变 应 变 应 变率 70率 100率 200率 250应 力 (MPa)8060402000.00000.00010.00020.00030.00040.0005时 间 (s)图19-17 四种应变率下的应力－时间曲线

160 应 变率 70 应 变率 200 140 应 变率 100 应 变率 250120应 力 (MPa)10080600.0000.0050.0100.0150.0200.0250.0300.0350.0400.045应 变图19-18 四种应变率下的应力－应变曲线 19-24

300250200 应 变 应 变 应 变 应 变率 70率 100率 200率 250应 变率 (s-1)150100500-50-1000.00000.00010.00020.00030.00040.0005时间(s)图19-19 四种设计应变率下的真实应变率－时间曲线

350300250200

应 变率 70 应 变率 200 应 变率 100 应 变率 250应 变率 (s-1)150100500-50-1000.0000.0050.0100.0150.0200.0250.0300.0350.0400.045应 变图19-20 四种应变率下的应变率－应变曲线

19-25

在不同的应变率下，材料表现出不同的应力－应变关系，体现出本构模型中应变率的作用。而且应力应变水平相对于应变率变化较大，在较高的应变率作用下，试件发生了较大的变形，承受的应力也较大。

以上结果说明，选用Johnson-Cook强化模型能很好的反映不同应变率作用下材料的力学行为，对于金属材料的冲击问题是适用的。在此基础上建立的用户材料子程序是可靠的。

考虑温度影响的强化模型分析

下面将在UMAT中考虑塑性耗散向热量转化所引起的温度的改变，从而影响本构模型中的应力－应变关系。这一过程是非耦合的，看作是绝热过程，即热量在单元中产生，引起单元温度的变化，但是相邻单元之间没有热量的传递。在应变率250 s件最终的温度场如图19-21所示。

－1

下试

(a) t=1.98×10-4 s时的温度场

(b) t=3.0×10-4 s时的温度场

(c) t=3.0×10-4 s时的温度场 图19-21 试件上温度场的变化

19-26

取塑性耗散比为0.95，即95％的塑性耗散将转化为热量。这一比例对于大多数金属是合适的。从上图中看出，早期应力波没有到达试件时，试件上的温度没有变化，经历塑性变形后，塑性耗散开始向热量转化，引起单元温度的变化。试件上的温升较为均匀，变形结束时试件上的最高温升约为1.9℃。

为了使温度对屈服应力的影响明显，取Johnson-Cook模型中参数m?0.5（通常对于大多数金属m值在1.0左右），这相当于夸大了模型中的温度软化效应。

将考虑温度效应时单元的应力应变与前面没有考虑温度效应的情况做一下比较，试件表面中间单元的温度曲线如图19-22所示。曲线中发生塑性变形后，单元的温度逐渐上升。

29.0 温度28.5温 度 (?C)28.027.527.00.00000.00010.00020.00030.00040.0005时 间 (s)

图19-22 单元温度的变化

在温度的影响下单元的应力、应变历史曲线如图19-23所示，并同前面不考虑温度作用的情况进行了对比。

图中随着单元温度的升高，应力较前面不考虑温度的情况要略低，而应变则略大。这说明单元表现出温度软化的特征。

应力－应变曲线表示如下（图19-24），在单元进入塑性后，温度升高的同时，屈服强度和硬化率要比前面不考虑温度的情况要小。Johnson-Cook模型中最后一项关于温度软化的情况在有限元模拟中得到了验证，同时也说明UMAT能正确的考虑温度对材料本构关系的影响。

19-27

0.05 考虑温度 不考虑温度0.040.03应 变0.020.010.000.00000.00010.00020.00030.00040.0005时 间 (s)

(a) 应变－时间曲线

考虑温度 不考虑温度140120100应 力 (MPa)8060402000.00000.00010.00020.00030.00040.0005时 间 (s)(b) 应力－时间曲线

图19-23 单元应力、应变历史

19-28

1601401201008060402000.000 考虑温度 不考虑温度应 力 (MPa)0.0050.0100.0150.0200.0250.0300.0350.0400.045应 变图19-24 考虑温度效应与不考虑温度效应的比较

18.4.4 三维动态分析

本例题进行的三维动态分析也相当成功，试件在模拟过程中各处应力都比较均匀，最大应力和最小应力的差别不超过1MPa。图19-25为试件的Mises应力分布。

图19-25 三维有限元模拟中试件的Mises应力分布

19-29

变形过程中单元的应力、应变以及应变率历史都与二维的情况非常吻合，如图19-26所示。

140120100800700600500400300200401002000.00000-1000.0005应力 三维数值模拟 二维数值模拟80600.00010.00020.00030.0004时 间 (s)图19-26 三维有限元模拟的应力、应变率历史

下面是应力－应变曲线：

140120100 三维数值模拟 二维数值模拟应 力 (MPa)8060402000.0000.0050.0100.0150.0200.0250.0300.0350.0400.045应 变图18-27 三维有限元模拟中的应力－应变曲线

19-30

应 变率 (s-1)应 力 (MPa)应 变率 三维数值模拟 二维数值模拟

这是对编写的UMAT用于三维实体单元的一个验证。UMAT虽然是从一维Johnson-Cook模型中建立起来，但是在UMAT中率相关塑性公式将它扩展到了三维情况，所以能用于三维实体单元也是意料之中。

虽然取得的结果是一致的，但是计算时间上却有很大差别，在个人计算机上，使用轴对称二维网格完成一次计算只需要10分钟，然后完成一次三维网格的计算需要3个小时。相比而言使用轴对称二维模型要经济得多。

18.5 UMAT的Fortran程序

18.5.1 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 ----------------------------------------------------------- C UMAT FOR JOHNSON-COOK MODEL

C ----------------------------------------------------------- C PROPS(1) - YANG'S MODULUS C PROPS(2) - POISSON RATIO

C PROPS(3) - INELASTIC HEAT FRACTION

19-31

C PARAMETERS OF JOHNSON-COOK MODEL: C PROPS(4) - A C PROPS(5) - B C PROPS(6) - n C PROPS(7) - C C PROPS(8) - m

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

C ELASTIC 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

C ELASTIC STIFFNESS C

DO 20 K1=1,NTENS DO 10 K2=1,NTENS

DDSDDE(K2,K1)=0.0 10 CONTINUE 20 CONTINUE C

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

19-32

50 CONTINUE C

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

C RECOVER ELASTIC AND PLASTIC STRAINS C

DO 80 K1=1,NTENS

EELAS(K1)=STATEV(K1)+DSTRAN(K1) EPLAS(K1)=STATEV(K1+NTENS) 80 CONTINUE

EQPLAS=STATEV(1+2*NTENS) C

C CALCULATE 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) 90 CONTINUE

SMISES=SQRT(SMISES/TWO) C

C CALL USERHARD SUBROUTINE, GET HARDENING RATE AND YIELD STRESS C C

CALL USERHARD(SYIEL0,HARD,EQPLAS,PROPS(4)) C DETERMINE IF ACTIVELY YIELDING C

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

C MATERIAL RESPONSE IS PLASTIC, DETERMINE FLOW DIRECTION C

SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE ONESY=ONE/SMISES DO 110 K1=1,NDI

19-33

FLOW(K1)=ONESY*(STRESS(K1)-SHYDRO) 110 CONTINUE

DO 120 K1=NDI+1,NTENS

FLOW(K1)=STRESS(K1)*ONESY 120 CONTINUE C

C READ PARAMETERS OF JOHNSON-COOK MODEL C

A=PROPS(4) B=PROPS(5) EN=PROPS(6) C=PROPS(7) EM=PROPS(8) C

C NEWTON 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 130 CONTINUE

WRITE(6,2) NEWTON

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

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

EFFHRD=EG3*HARD1/(EG3+HARD1) C

C CALCULATE STRESS AND UPDATE STRAINS C

DO 150 K1=1,NDI

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

19-34

EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO 150 CONTINUE

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 160 CONTINUE

EQPLAS=EQPLAS+DEQPL

SPD=DEQPL*(SYIEL0+SYIELD)/TWO RPL = PROPS(3)*SPD/DTIME C

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 ENDIF C

C STORE STRAINS IN STATE VARIABLE ARRAY C

DO 310 K1=1,NTENS

STATEV(K1)=EELAS(K1)

STATEV(K1+NTENS)=EPLAS(K1) 310 CONTINUE

19-35

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

C GET PARAMETERS, SET HARDENING TO ZERO C

A=TABLE(1) B=TABLE(2) EN=TABLE(3) HARD=0.0 C

C CALSULATE 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

18.5.2 UMATHT（包含材料的热行为）

SUBROUTINE UMATHT(U,DUDT,DUDG,FLUX,DFDT,DFDG,STATEV,TEMP,

$ DTEMP,DTEMDX,TIME,DTIME,PREDEF,DPRED,CMNAME,NTGRD,NSTATV,

$ PROPS,NPROPS,COORDS,PNEWDT,NOEL,NPT,LAYER,KSPT,KSTEP,KINC) C

INCLUDE 'ABA_PARAM.INC'

19-36

C

CHARACTER*80 CMNAME C

DIMENSION DUDG(NTGRD),FLUX(NTGRD),DFDT(NTGRD),

$ DFDG(NTGRD,NTGRD),STATEV(NSTATV),DTEMDX(NTGRD),TIME(2), $ PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3) C C

COND=PROPS(1) SPECHT=PROPS(2) C

C INPUT SPECIFIC HEAT DUDT=SPECHT DU=DUDT*DTEMP U=U+DU C

C INPUT FLUX = -[K]*{DTEMDX} DO I=1, NTGRD

FLUX(I)=-COND*DTEMDX(I) END DO C

C INPUT ISOTROPIC CONDUCTIVITY C

DO I=1, NTGRD

DFDG(I,I)=-COND END DO C

RETURN END

本章主要内容引自：

卢剑锋，冲击载荷作用下材料和结构力学行为有限元模拟，清华大学硕士论文，2003。

19-37

本文来源：https://www.bwwdw.com/article/9ky7.html