各向异性屈服准则的UMAT子程序二次开发研究

乔顺成，吴建军，展学鹏

（1.中航工业西安飞机工业集团有限责任公司，陕西 西安 710072；2.西北工业大学 机电学院，陕西 西安 710072）

ABAQUS以其强大的非线性迭代计算和前后处理功能而广泛应用于材料弹塑性有限元分析，它是现阶段在各个工程领域广泛使用的大型通用有限元软件之一。ABAQUS[1]有大量的单元库和求解模型供用户使用，而且用户也能够通过这些模型求解绝大多数问题。但是实际问题是相当复杂的，ABAQUS不可能直接处理所有可能的问题[2-4]，所以，有必要给用户提供二次开发的接口来解决实际中出现的新问题。在用户自定义材料（UMAT即 User-defined Material Mechanical Behavior）研究中，很多学者通过UMAT子程序二次开发，将新的本构模型用于实际问题的有限元分析。例如：李平[5]等通过数值分析模拟了普通碳钢连续热轧的过程，将率相关的各向同性硬化热变形本构模型通过UMAT子程序用于有限元仿真，分析了轧制过程中的温度和应力、应变之间的耦合关系。杨曼娟[6]提出了五参数的柳玉起屈服准则二次开发开发程序，并将Rankine准则的Mohr-Coulomb模型通过UMAT子程序嵌入ABAQUS，最后通过单元测试和实际工程算例验证了二次开发的正确性和适用性，惟一不足的就是作者并没有附上二次开发的Fortran程序代码。唐榕蔚[7]将双剪统一强度理论通过UMAT子程序嵌入到 ABAQUS软件中，用于岩土材料的弹塑性有限元分析，作者最后还附上了二次开发程序代码，为以后新本构模型的二次开发研究者提供了参考和指导。除此之外，还有许多学者通过有限元二次开发来解决一些实际的工程问题[8，9]。

ABAQUS嵌入了Mises各向同性屈服准则、Tresca各向同性屈服准则、Hill48各向异性屈服准则，前两个各向同性屈服准则只能用于各向同性有限元分析，Hill48各向异性屈服准则可以用于各向异性有限元分析。但是，Hill48屈服准则形式简单，对屈服面描述粗略，不涉及材料的晶体结构理论和晶体塑性理论，仅对r＞1（r指的是材料厚向异性指数）的各向异性材料屈服状态和塑性变形预测较好[10]，所以符合Hill48屈服准则的各向异性有限元分析必然存在局限和不足，不能精准地预测复杂屈服状态，不能精确地描述复杂的塑性变形，许多新特征、新规律不能被发现。例如，倪向贵等[11]将Yld89、Yld91屈服准则本构模型通过二次开发应用于有限元分析，研究了符合不同各向异性屈服准则（Hill48、Yld89、Yld91）本构关系的方盒件拉深成形，并将三种屈服准则的模拟结果与实验进行了比较，研究证明符合Hill48屈服准则的方盒件拉深成形仿真与实验结果相差最大，成形精度最差。

为了弥补之前屈服准则的缺陷和不足，Barlat等[12]提出了Yld2004-18p屈服准则，该准则需要RDTD平面上沿轧制方向每15°的单向屈服应力σY、厚向异性指数r，还有双向拉伸屈服应力σb和厚向异性指数rb，可以精确地预测拉深成形制件不同方位6或8个“制耳”，也可以精确地预测面内流动应力和厚向异性指数的变化，更加全面地反映了塑性流动各向异性和应力-应变响应，是目前描述FCC（面心立方）和BCC（体心立方）材料各向异性性能最准确的屈服准则之一。目前，大多数有限元软件都没有嵌入这个屈服准则本构模型。

为了使得材料加工工艺和塑性成形有限元分析更加精确，更加精准地预测屈服状态，更加合理地描述塑性变形行为，拓展ABAQUS在各向异性有限元分析领域的应用，很有必要利用ABAQUS提供的UMAT子程序二次开发接口，建立编译运行环境，采用Fortran语言编程，将高级各向异性屈服准则本构模型嵌入ABAQUS，并应用于弹塑性有限元分析。

本文构建了Yld2004-18p各向异性屈服准则本构模型，计算了Yld2004-18p屈服准则的各向异性系数，通过UMAT子程序二次开发，将Yld2004-18p各向异性屈服准则本构模型嵌入ABAQUS，结合有限元模型进行实例分析，验证UMAT子程序二次开发的正确性，同时提出一个通用的、柔性的二次开发结构模式。

1 Yld2004-18p各向异性屈服准则及其系数

为弥补之前的Yld89、Yld91、Yld96屈服准则的不足，F.Barlat等[12]提出了更高级的屈服准则——Yld2004-18p屈服准则，该准则给应力偏张量s又加了两个线性变换，这两个线性变换共包含有18个各向异性系数。该屈服准则如下：

式中：指数m与材料的晶体结构类型有关，对于BCC（体心立方）材料，m=6，对于FCC（面心立方）材料；S˜′，S˜″是与应力偏张量s˜′，s˜″有关的主值；s˜′，s˜″由应力偏张量s通过线性变换得到，如下式：

式中：C′和C″是应力张量线性变换矩阵，如下式：

T也是线性变换，如下式：

当C′=C″，Yld2004屈服准则等同于Yld91屈服准则；当18个各向异性系数ci=1～18全部为1，m=2（或4）时，它也可以退化到Mises各向同性屈服准则。

2 Yld2004屈服准则本构模型的UMAT子程序二次开发

2.1 与Yld2004屈服准则有关的变量推导

构建弹-塑性本构模型的首要条件是确定屈服函数，及与屈服函数有关的变量，它们与硬化法则、应力更新算法是构建弹-塑性本构模型的必要条件。为了将Yld2004屈服准则本构模型通过有限元二次开发--UMAT子程序嵌入ABAQUS软件，屈服条件、屈服函数对应力分量σij的一次导数、二次导数是必不可少的。它们在整个弹-塑性本构模型发挥着关键性的作用，是更新迭代步增量（例如应力增量、应变增量、塑性应变增量、切线模量、塑性参数增量等）、输出变量的重要环节。

对于Yld2004屈服准则本构模型，结合等式（1），可以得到等效应力σ¯的表达式：

结合等式（1）、（7），等效应力σ¯对应力分量 σij

其中：

等式（9）是Φ对σij的“链式求导”，可结合等式（1）-（6）求得。等效应力σ¯对应力分量 σij的二阶导为：

其中：

以上这些公式将会在下文应用于UMAT子程序本构模型。

2.2 弹-塑性率本构方程

材料屈服状态的判断根据是复杂应力状态下的等效应力与参考方向的单向屈服应力之间的关系，如此可以建立屈服条件：

从初始加载到屈服之前，材料发生弹性变形，符合广义Hook定律，弹性响应对应于弹性应变率，满足如下关系：

式中：

λ*，μ是独立的材料常数，称为拉梅（Lame）常数，可以用更接近于物理度量的常数表示它们：弹性模量E、泊松比v，如下式：

塑性应变率由流动法则确定，常表示为塑性流动势能Ψ的形式：

令：

从几何角度上，塑性流动方向与屈服面的法线方向相同。当塑性加载时，λ˙＞0，应力保持在屈服表面F=0，也可以用一致性条件表示：F˙=0，通过“链规则”扩展，推导得到如下：

2.3 UMAT子程序基本流程的构建

由于ABAQUS不包含高级的Yld2004屈服准则本构模型，所以将高级的Yld2004屈服准则本构模型应用到球形压痕有限元分析，就必须开发相应的屈服准则本构模型UMAT子程序。

根据ABAQUS提供的二次开发接口要求，用Fortran语言编写程序，得到.for格式的UMAT子程序，将其嵌入ABAQUS，就可以实现调用。UMAT子程序二次开发使用户能够自定义ABAQUS材料模型库中没有的材料模型，分析实际中更具体、更复杂的新问题。UMAT子程序的功能非常强悍，可以根据材料特性，自定义新的材料本构；可以实现新问题的有限元分析；它与主程序相互调用，可以用于任何加载或卸载阶段的分析；可以在材料单元每个积分点上进行调用等。

上述第2.1、2.2节已经完成了UMAT子程序本构模型需要的主要变量、弹塑性率本构方程的推导过程。为了方便编程，直观地表示迭代循环的步骤，需要构建UMAT子程序本构模型，将第2.1、2.2节与完全隐式向后Euler图形返回算法结合，它的简要过程如下：

（1）初始状态：

初始值设置：k=0，Δλ=0；

（2）在第k次迭代检查屈服条件和收敛性：

如果 F（k）≤tolerance 且‖R（k）‖≤tolerance，则储存变量，此时间步结束；否则，继续下一步骤；

（3）在第k次迭代，计算中间变量：

（5）更新变量：

返回第（2）步。

2.4 UMAT子程序结构模式及编程

作为独立的程序模块，UMAT子程序既能被ABAQUS主程序调用，也能被单独编译或存储。UMAT子程序的编程应该遵循结构化程序设计的思想，即：灵活使用三种基本程序结构（顺序结构、循环结构、选择结构），按照自顶向下的设计思路，分解程序功能，使程序合理的模块化、分支化，可使子程序足够简单，柔性化，具有互换性，更加通用。

本文开发的Yld2004屈服准则本构模型UMAT子程序，从整体上看，都具有以下模块：①ABAQUS规定的题名声明语句；②ABAQUS的接口参数声明语句；③自定义的变量声明语句；④程序主体，包括主程序和5个嵌套的功能分支子程序；⑤程序结束语和返回语句。在开发的UMAT子程序中，一致切线模量矩阵C的求逆运算、与材料硬化行为有关的硬化函数运算、与屈服函数有关的三个变量：等效应力，一阶导，二阶导，它们是相互独立的运算过程，可以将它们分解为独立的功能模块——分支子程序，在程序内部实现自调用。结合第2.3节内容，我们提出的UMAT子程序结构模式如图1所示。等效应力σ¯模块、一阶导模块、二阶导模块，它们只与选择的屈服准则有关；屈服应力σY模块只与加载过程中材料的硬化行为有关；一致切线模量矩阵C的求逆模块只与矩阵C有直接关系，这5个功能模块之间相互独立，运算过程没有直接关系。针对不同的屈服准则，只需改变与屈服准则有关的模块（等效应力σ¯模块、一阶导模块、二模块）就可以嵌入新的屈服准则；针对不同的材料硬化行为，只需改变与硬化法则有关的模块（屈服应力σY模块）就可以研究新的材料硬化行为。因而，上述UMAT子程序结构模式具有互换性，是通用的，它可以将具有新的屈服准则、新的硬化函数的本构模型嵌入ABAQUS.

图1 通用的UMAT子程序结构模式

表1 Yld2004屈服准则本构模型的UMAT子程序材料常数定义

UMAT子程序的编程需要定义的材料参数有：弹性模量E，泊松比v，屈服准则各向异性系数（Yld2004屈服准则有18个：ci=1～18），硬化模型参数（本文所研究的材料硬化符合Hollomon硬化法则，阶导表达式为：σY=K（ε0+¯p）n，包含 K，ε0，n 三个参数）。Yld2004屈服准则本构模型的UMAT子程序材料常数被定义在表1，同时，UMAT子程序的编程需要定义一个状态变量数组，用于存放与求解过程有关的的变量，这些变量将在求解过程中不断地被更新，并传递到主程序，然后在下一个增量步开始时再传递给UMAT子程序。本文开发的UMAT子程序包含13维的状态变量数组：弹性应变分量存储在1～6维，塑性应变分量存储在7～12维，等效塑性应变存储在第13维。其他局部变量的名称定义应该简洁明了、易辨识。UMAT子程序中数组与数组、矩阵与数组、矩阵与矩阵的矢量乘积编程规则应该特别注意。

UMAT子程序在分析计算过程中会与ABAQUS主程序结合，进行数据的传输、更新。为了实现这些过程，UMAT子程序与ABAQUS主程序共享一些变量，即ABAQUS用户子程序接口，接口定义在UMAT子程序的题名中，对数据的传输、更新和相互调用起着桥梁纽带作用。以下就是UMAT子程序与主程序的接口共享变量：

3 实例有限元分析

为了验证UMAT子程序开发过程的正确性，必须把UMAT子程序嵌入ABAQUS软件中进行实例模型有限元分析。根据第1节内容，当Yld2004屈服准则的各向异性系数都为1且指数m=2（或4）时，它可以退化为Mises各向同性屈服准则。通常建立简单的单拉或者单压模型，将Yld2004屈服准则的UMAT子程序与单拉或者单压模型耦合，设置系数全为1、指数m为2或4，将UMAT子程序有限元分析得到的结果，与用ABAQUS自带的Mises屈服准则计算结果进行对比，可以直接简洁地验证Yld2004屈服准则的UMAT子程序是否满足精度要求，是否在误差允许范围内。

本研究建立了简化的单向压缩有限元模型，如图2所示，有限元模型一端完全固定，另一端施加合适的加载力（均布载荷为320MPa），分析步类型为静力、通用，设置UMAT子程序中各向异性系数都为1即表1材料常数6～23都为1，且指数m=2（或4）即表1材料常数24为2（或4），进行有限元分析，完成有限元分析之后，比较各自的云图，结果如图3所示；选取主要的变形路径提取等效应力、等效塑性应变进行对比，结果如图4所示。

图2 简化的单向压缩有限元模型

由图3所示，将Yld2004屈服准则退化为Mises屈服准则模拟的结果与ABAQUS自带的Mises屈服准则模拟的结果，两者等效应力分布云图、等效塑性应变分布云图整体上都保持一致。

由图4可知，将Yld2004屈服准则退化为Mises屈服准则模拟的结果与ABAQUS自带的Mises屈服准则模拟的结果，两者单元结点上等效应力、等效塑性应变值基本相同，精度保持在99.9%左右。

4 结束语

由于ABAQUS软件没有高级的Yld2004屈服准则本构模型，使之受限于精密塑性成形和材料加工有限元分析。因此，本研究推导了与Yld2004-18p屈服准则有关的重要变量，结合完全隐式的Euler图形返回算法，构建了Yld2004-18p屈服准则本构模型，用Fortran语言编程，通过ABAQUS提供的UMAT子程序二次开发接口，将Yld2004-18p屈服准则本构UMAT子程序嵌入ABAQUS软件。建立单压有限元模型，将Yld2004-18p屈服准则UMAT子程序退化到Mises屈服准则，用于实例有限元分析；将分析结果与ABAQUS自带的Mises屈服准则分析结果作比较，两者保持高度的一致性，验证Yld2004-18p屈服准则UMAT子程序二次开发的正确性，同时提出了一个通用的、柔性的屈服准则二次开发结构模式。

图3 有限元模拟变形图对比

图4 变形路径上节点变量对比