各向异性屈服函数子程序库开发

各向异性屈服函数的用户子程序库“统一塑性材料模型驱动器（UMMDp）”的开发

1. 引言

金属板的各向异性塑性性能取决于其制造工艺。在数值模拟中，这些塑性性能通过屈服函数来表示。在关联流动法则的假设下，屈服函数描述了材料在塑性状态下的应力分量与塑性应变增量方向之间的关系。

因此，为了准确模拟金属板在成形工艺中的塑性变形，必须采用能够准确再现实际金属板塑性变形的合适各向异性屈服函数。过去，已提出多种各向异性屈服函数来描述金属板的塑性变形特性 [1]。

已将多种各向异性屈服函数集成到专为金属板成形过程模拟定制的商业有限元（FE）程序中 [2]。另一方面，通用有限元（FE）程序包含的各向异性屈服函数较少，通常仅限于经典的模型，例如希尔1948模型 [3]。或者，几乎所有的通用有限元（FE）程序都提供了用户子程序接口，使用户能够将自定义的材料模型集成到代码中。然而，对于普通用户而言，实现各向异性屈服函数可能较为困难，因为编写适当的子程序需要用户具备非线性连续介质力学、塑性理论以及有限元法数值过程的深入知识。

本文介绍了适用于主要商业有限元代码的各向异性塑性模型用户子程序库的开发。提出了该库的框架，并展示了使用所提出子程序进行的简单验证测试。

2. 应力积分和一致切线模量

处理弹塑性变形的典型非线性有限元代码采用次弹性型材料模型。图1展示了次弹性模型在有限元代码中的作用。有限元代码中的材料模型具有两个作用：应力积分和一致切线模量的推导。它们在每次迭代中均需计算单元的各个积分点。有限元代码的主程序利用更新后的应力状态评估整体结构的平衡。如果残余力的范数不可接受，则通过牛顿‐拉夫森方法利用一致切线模量预测一个可减小残余力的修正位移场。

本文介绍了在有限元代码中实现任意各向异性屈服函数与各向同性硬化的基本方程。这些方程采用典型有限元教材中使用的沃伊特符号表示，并采用向后欧拉法进行应力积分。

2.1 向后欧拉法的基本方程

在向后欧拉法中，描述塑性变形的所有方程在增量n+1结束时得到满足。增量n+1结束时，应力位于后续屈服面上的状态由表示。
$$ {\sigma^{n+1}} $$
$$ 0 = f(\sigma^{n+1}) - Y(p^{n+1}) $$ (1)

这里，屈服函数的类型 $ f $ 是任意的，$ p^{n+1} $ 是增量中等效塑性应变，以及 $ Y(p) $ 是由 $ p $ 驱动的各向同性硬化函数。假设采用关联流动法则，塑性应变增量 $ \Delta \varepsilon^p $ 表示为
$$ \Delta \varepsilon^p = \Delta p \cdot \mathbf{N} $$
$$ \mathbf{N} = \frac{\partial f}{\partial \sigma} $$ (2)

塑性应变增量的方向与增量后的后续屈服面的外法线方向一致。矢量 ${N}$ 表示屈服函数相对于应力分量的一阶偏导数。

如果弹性应变足够小，则可以采用弹塑性应变增量的加法分解。弹性应变增量 $ \Delta \varepsilon^e $ 可表示为
$$ \Delta \varepsilon^e = \Delta \varepsilon - \Delta \varepsilon^p $$ (3)

因此，增量n+1中的应力可以表示为
$$ \sigma^{n+1} = \sigma^n + \mathbf{C}_e : (\Delta \varepsilon - \Delta p \cdot \mathbf{N}) $$ (4)

这里，$\mathbf{C} e$ 是广义胡克定律的矩阵。试应力 $ \sigma {\text{trial}} $ 是假设总应变增量为弹性时的临时应力状态。

方程(1)和(4)是向后欧拉法中应力积分的基本关系。应力积分的示意图如图2所示。

2.2 应力积分

在应力积分中，求解应力增量和等效塑性应变的增量以满足方程(1)和(4)。首先，定义来自这些方程的残差标量和向量为
$$ g_1 = f(\sigma^{n+1}) - Y(p^{n+1}) $$ (5)
$$ \mathbf{g}_2 = \sigma^{n+1} - \sigma^n - \mathbf{C}_e : (\Delta \varepsilon - \Delta p \cdot \mathbf{N}) $$ (6)

分别。且 $ g_1 $ 和 $ \mathbf{g}_2 $ 必须通过牛顿–拉夫逊迭代收敛至零。

为了在牛顿–拉夫森迭代中获得应力修正和等效塑性应变（$ \delta \sigma $ 和 $ \delta \Delta p $），将方程(5)和(6)线性化如下。
$$ 0 = \frac{\partial f}{\partial \sigma} : \delta \sigma - \frac{\partial Y}{\partial p} \delta \Delta p $$ (7)
$$ 0 = \delta \sigma + \mathbf{C}_e : (\delta \Delta p \cdot \mathbf{N}) + \mathbf{E} : \mathbf{g}_2 $$ (8)

此处省略了下标n+1 。[E]表示如下。
$$ \mathbf{E} = \mathbf{I} + \frac{\partial \mathbf{N}}{\partial \sigma} : \mathbf{C}_e \cdot \Delta p $$ (9)

[I]是单位矩阵。在线性化过程中，应变增量不被视为变量，因为在应力积分中应变增量是显式给定的。

求解联立方程(7)和(8)，可得到 $ \delta \sigma $ 和 $ \Delta p $ 的修正值为
$$ \delta \Delta p = -\frac{\mathbf{g}_2 : \mathbf{E}^{-1} : \mathbf{N}}{\frac{\partial f}{\partial \sigma} : \mathbf{E}^{-1} : \mathbf{N} + \frac{\partial Y}{\partial p}} $$ (10)
$$ \delta \sigma = -\mathbf{E}^{-1} : (\mathbf{g}_2 + \mathbf{C}_e : (\delta \Delta p \cdot \mathbf{N})) $$ (11)

上日期，$ \sigma^{n+1} $ 以及 $ \Delta p $，被迭代直至残差 $ g_1 $ 和 $ \mathbf{g}_2 $ 收敛。

2.3 一致切线模量

一致切线模量用于组装整体结构的刚度矩阵，且必须与应力积分方案的算法保持一致，以获得二次收敛结构平衡中的收敛性。该矩阵对于静态/隐式有限元代码是必需的，而在动态/显式有限元代码中则不需要。

一致切线模量可以通过推导应力和应变增量的扰动之间的关系来获得，$ \delta \sigma $ 和 $ \delta \varepsilon $。此处，将应变增量视为变量，以推导一致切线模量，这与上述应力积分相反。方程(1)和(4)的全微分如下所示。
$$ 0 = \frac{\partial f}{\partial \sigma} : \delta \sigma - \frac{\partial Y}{\partial p} \delta \Delta p $$ (12)
$$ \delta \sigma = \mathbf{C}_e : (\delta \varepsilon - \delta \Delta p \cdot \mathbf{N} - \Delta p \cdot \frac{\partial \mathbf{N}}{\partial \sigma} : \delta \sigma) $$ (13)

该方程可按如下方式求解。
$$ \delta \sigma = \mathbf{D}_c : \delta \varepsilon $$ (14)
这里，$ \mathbf{D}_c = [\mathbf{C}_e^{-1} + \frac{\partial \mathbf{N}}{\partial \sigma} \cdot \Delta p]^{-1} $。从上述可知，只要给定屈服函数及其一阶和二阶偏微分，即可实现应力的积分和一致切线模量的推导。

3. UMMDp框架

如上所述， (1) 次弹性本构模型的材料模型用户子程序在所有商业有限元代码中的作用是相同的。(2) 在应力积分和一致切线的推导中，需要屈服函数及其一阶和二阶偏微分。然而，其他数值方法是通用的，并且与屈服函数的类型无关。

基于这两个特性，材料模型可以统一为数值方法的公共部分，从而将多种有限元代码和屈服函数外部化。

图3展示了统一塑性材料模型驱动器（UMMDp）的框架。该结构的中心是一个核心程序，它是弹塑性计算中数值过程的公共部分。

该核心程序由各有限元代码的用户子程序调用。每个用户子程序都起到将UMMDp核心连接到有限元代码的“插件”作用。用户子程序将各有限元代码的变量转换为UMMDp核心所通用的形式。

在图3右侧的分支中，各种屈服函数被模块化为子程序，其中屈服函数及其一阶和二阶偏微分作为核心程序的参数进行计算。此外，如图3左侧分支所示，包括运动学硬化准则和联合硬化模型在内的各种硬化准则也被模块化，以支持未来的扩展。

UMMDp的用户子程序库独立于特定的有限元代码，为用户提供对各种屈服函数的广泛可扩展性。

4. 日本非线性CAE协会中的协作活动

表1列出了各种主要有限元代码中用于超弹性材料模型的用户子程序中的变量。变量的名称和数组类型因代码而异，但在此列表中可以明显看出其相似性。有关用户子程序在商业有限元代码中的具体应用，已在各软件供应商举办的会议上进行了讨论。

日本非线性CAE协会（JANCAE）是一个非营利组织。JANCAE的主要活动之一是举办非线性CAE培训课程，该课程每年举行两次，旨在加深用户对CAE的理解。从2001年到2017年，共有约5000名工程师参加了该培训课程。自2005年起，JANCAE还组织了“材料建模工作组”，每年有来自大学、私营企业和软件供应商的30多名工程师和研究人员参与。

UMMDp的所有程序均由材料建模工作组的参与者在志愿性质的基础上开发。特别是，插件子程序由软件供应商的技术支持工程师编写。由于该插件，UMMDp可适用于Abaqus、ANSYS、ADINA、LS-DYNA和Marc。另一方面，各向异性屈服函数的模块由制造业的CAE用户编写。以下过去提出的屈服函数已集成到UMMDp中：Hill的1948[3]和1990[5], Gotoh双二次型[6], Hu双二次型[7], Yoshida六阶多项式[8], Barlat的Yld2000[9]和Yld2004[10], Banabic的BBC2005[11]和BBC2008[12], Cazacu的CPB2006[13], Karafillis与Boyce[14], 以及Vegter样条函数[15]。

表1. 通用有限元代码中代表性变量列表

Code	用户子程序	应力	应变增量	一致切线	应力分量数	状态变量	状态变量数	退出代码
Abaqus	umat	应力	DSTRAN	DDSDDE	NDI, NSHR	STATEV	NSTATEV	XIT
ADINA	ucmat2/3	应力	dep s d	已修复。数组	LGTH1 END	—	—	—
ANSYS	usermat	应力	应变增量	塑性切线模量	法向数量, 切向数量	状态变量	nStatev	EXIT
LS-DYNA	umat/utan	sig eps	es	参见‘eltype’	—	hsv	nhv	—
Marc	hypela2	S DE D	—	—	NDI, NSHEAR	T, DT	NSTATS	QUIT

5. UMMDp的简单验证测试

这里，我们通过一个简单验证测试来展示UMMDp在有限元代码中的准确工作情况。测试的边界条件和材料属性如图4(a)所示。流动应力为常数（完全塑性），且杨氏模量相对于流动应力足够大。图4(b)显示了施加给单元的位移历史。由于采用了关联流动法则[16]的假设，有限元代码中分析的应力点应沿规定的屈服轨迹移动。所分析的应力点会收敛到屈服轨迹的外法线与塑性应变增量方向一致的点。

在图5中，实线表示由Yld2004[10]定义的规定屈服轨迹，计算点表示使用UMMDp分析的应力点。通过应力图与规定屈服轨迹的精确一致性验证了UMMDp的精度。除了简单验证外，工作组的参与者合作解决了实际板料成形问题，作为UMMDp的工作实例。

6. 结论

介绍了用于各种各向异性屈服函数的用户子程序库的开发。该库适用于主要的有限元代码。通过使用此库，有限元用户可以轻松地将最近提出的或他们自己的屈服函数实现到有限元代码中。

UMMDp 和一些示例输入数据文件在网站[4]上向公众开放。