这篇文章讲解如何使用晶体塑性有限元方法(CPFEM)进行不同晶格材料以及多相材料的变形模拟，CPFEM是基于商业有限元软件ABAQUS完成的建模，晶体塑性本构模型是使用的开源的UMAT用户子程序(源码和inp文件见附件)。采用CPFEM模拟了面心立方结构(FCC)、体心立方结构(BCC)和密排六方结构(HCP)的单晶、多晶及多相材料受到外部载荷时的力学响应。基于滑移原理的晶体变形理论，随着变形的进行各滑移系统的临界剪应力都会增大，CPFEM将捕捉到材料的力学响应(应力-应变曲线)。这些应力-应变数据有助于从根本上理解晶粒尺度下金属变形的性质。

首先我们从一个简单的FCC晶格材料的例子入手，讲解如何进行有限元模型的创建，从完全新手的角度出发，一步步讲解如何建模，赋予材料和处理仿真结果。

本文章包括以下八个部分：

1) FCC晶格材料的变形模拟-单晶体

2) FCC晶格材料的变形模拟-多晶体

3) BCC晶格材料的变形模拟-单晶体

4) BCC晶格材料的变形模拟-多晶体

5) HCP晶格材料的变形模拟-单晶体

6) HCP晶格材料的变形模拟-多晶体

7) 多相材料的变形模拟

8) 参考资料

1. FCC晶格材料的变形模拟-单晶体

几何模型

如图1.1a在草图里绘制R0.015mm的圆形，拉伸0.05mm，最后得到如图1.1b所示的圆柱体(R0.015mm&H0.05mm)。

图1.1 a) 草图里绘制R0.015mm的圆形b) Z方向拉伸得到H0.05mm的圆柱体

材料模型

FCC晶格材料我们使用Huang的博士论文A user-material subroutine incorporating single crystal plasticity in the ABAQUS finite element program里的copper crystal，其160组材料参数如下：

图1.2 铜晶体的160组晶体塑性材料参数

将20*8的矩阵其转化为160*1的矩阵后，在Abaqus软件中，定义材料模型的具体参数如下：

图1.3 a) 定义材料的输出参数个数 b) 定义材料的晶体塑性参数

步长设置

步长设置的具体参数如下：

图1.4 步长设置参数

边界条件

具体边界条件如下：Z0面固定约束，Z1面沿z方向位移0.001mm。

图1.5 边界条件设置

网格划分

图1.6 网格划分前对部件进行剖分

网格大小0.002进行网格划分：

图1.7 网格划分

图1.8 更改网格类型

提交运算

图1.9 提交运算和完成运算的界面

后处理界面

应力应变分布：

图1.10 后处理应力应变分布

2. FCC晶格材料的变形模拟-多晶体

几何模型

通过NEPER软件生成圆柱体(R0.015mm&H0.05mm)里包含300个晶粒的几何模型，结果如图2.1所示，其中不同的颜色代表不同的晶粒。

图2.1 圆柱体(R0.015mm&H0.05mm)里包含300个晶粒的几何模型

网格划分

多晶体有限元模型的构建流程与单晶体有限元模型的构建流程不同，不是ABAQUS推荐的几何模型材料模型网格划分，而是采用了几何模型网格划分材料模型。

通过GMSH软件对圆柱体几何模型进行网格划分，结果如图2.2所示。(在附件里提供了构建好的inp文件，具体要如何进行网格划分和批量赋予材料参数，在Neper课程里有说明，欢迎大家学习)

图2.2 圆柱体几何模型网格划分的结果

材料模型

同样的，FCC晶格材料我们使用Huang的博士论文A user-material subroutine incorporating single crystal plasticity in the ABAQUS finite element program里的copper crystal，其160组材料参数如下：

图2.3 铜晶体的160组晶体塑性材料参数

将20*8的矩阵其转化为160*1的矩阵后，在ABAQUS软件中，用PYTHON脚本批量定义了材料模型的具体参数(随机生成晶体取向)，定义好材料模型的结果如图2.4所示。

图2.4 300个晶粒的定义完材料模型的结果

步长设置

步长设置的具体参数如下：

图2.5 步长设置参数

边界条件

具体边界条件如下：Z0面固定约束，Z1面沿z方向位移0.001mm。

图2.6 边界条件设置

提交运算

图2.7 提交运算和完成运算的界面

后处理界面

应力应变分布：

图2.8 后处理应力应变分布

查看后处理的塑性变形PEEQ都为0，这是因为使用晶体塑性模型的原因。我们可以在inp文件中添加想要输出的SDV参数，总的输出参数有12*10+5=125个，各SDV参数对应的存储结果如下图所示，如我们可以输出SDV121，即为所有滑移系的总剪切应变，可以把这个参数的结果看做PEEQ的结果。

图2.9 晶体塑性模型输出参数对照图

图2.10 后处理PEEQ和SDV121的分布

FCC晶格材料一共有12组滑移系，我们如果想查看1号滑移系的强度和剪切应变，可以输出的SDV1和SDV109，结果如图2.11所示，可以观察到各晶粒滑移系的开动情况。

图2.11 后处理SDV1和SDV109的分布

3. BCC晶格材料的变形模拟-单晶体

所有模块的设置与第1节"FCC晶格材料的变形模拟-单晶体";的创建过程类似，只是需要修改BCC晶格材料的材料参数。

材料模型

图3.1 晶体塑性模型材料参数对照图

BCC晶格材料我们使用论文《基于CPFEM的TA15钛合金高温塑性变形研究》里研究的β-Ti，其160组材料参数如下：

图3.2 β-Ti晶体的160组晶体塑性材料参数

注意，以上参数中，Props(1-21)是晶体材料的弹性常数，弹性常数的独立张量元数目至多只有21个。对不同的晶系的晶体(7大晶系)，因为对称性的关系，其独立的弹性常数是确定的。晶系的对称性越高，独立的张量元数目越少。六方晶系(HCP)，独立的张量元数目有5个C11 C12 C13 C33 C44，立方晶系(FCC/BCC)，独立的张量元数目有3个C11 C12 C44。

图3.3 160组晶体塑性材料参数的弹性部分

Props(25-56)是晶体材料的滑移系参数。Props(25)是滑移系族的个数，对FCC晶格材料有12个滑移系，只有1个滑移系族{111} <110>；对BCC晶格材料有48个滑移系，有3个滑移系族{110} <111>、{112} <111>、{123} <111>。

其余的材料参数在图3.1中有解释，将弹塑性参数的数值按照论文的数据进行设置，并将20*8的矩阵其转化为160*1的矩阵后，定义给材料模型。

其他模块的设置与第1节"FCC晶格材料的变形模拟-单晶体"的创建过程类似，在此不再赘述。

后处理界面

应力应变分布如图3.4所示。

图3.4 后处理应力应变分布

4. BCC晶格材料的变形模拟-多晶体

所有模块的设置与第2节"FCC晶格材料的变形模拟-多晶体"的创建过程类似，只是需要按照第3节的操作修改BCC晶格材料的材料参数，构建好的inp文件见附件，下面看看模型的变形结果。

后处理界面

应力应变分布：

图4.1 后处理应力应变分布

图4.2 后处理一些SDV结果的分布

5. HCP晶格材料的变形模拟-单晶体

塑性部分考虑滑移和孪晶的变形梯度如图5.1所示，总的塑性变形梯度来自等式右侧的两项：分别代表由滑移和孪晶对塑性变形梯度的贡献。与考虑滑移和孪晶的变形相比，滑移为主的塑性变形梯度忽略了孪晶对塑性变形梯度的贡献，如图5.2所示。

图5.1 包含滑移、孪晶的变形梯度分解示意图

图5.2 滑移变形为主的梯度分解示意图

Huang的晶体塑性子程序是基于滑移原理编写的，因此使用该子程序进行HCP晶格材料变形的模拟时需保证变形是以滑移为主的前提条件，如论文《纯钛塑性变形行为的晶体塑性有限元模拟》、《单晶纯钛的细观力学性能模拟》和《钛合金双态组织高温拉伸行为的晶体塑性有限元研究》里明确指出忽略了HCP晶格材料的孪生变形。

如考虑孪生和滑移两者的贡献，需修改Huang基于滑移原理的子程序，或者换DAMASK平台进行模拟， 链接里是DAMASK平台模拟镁合金的变形。

图5.3 DAMASK平台模拟镁合金变形的参数设置

下面，将采用Huang的基于滑移原理晶体塑性子程序，开展HCP晶格材料以滑移为主的塑性变形模拟。

所有模块的设置与第1节"FCC晶格材料的变形模拟-单晶体"的创建过程类似，只是需要修改HCP晶格材料的材料参数。

材料模型

HCP晶格材料我们使用论文《Ti-6Al-4V合金纳米压痕变形与高周疲劳行为CPFEM研究》里研究的α-Ti，其160组材料参数如下：

图5.4 α-Ti晶体的160组晶体塑性材料参数

注：四轴坐标系中晶向指数[u v t w]，晶面指数[h k i l]，通过下式获得直角坐标系中的晶向指数。《Ti3Al 单晶和双相片层 TiAl 合金塑性行为的 CPFEM 模拟》

图5.5 晶向指数、晶面指数的转换

其他模块的设置与第1节"FCC晶格材料的变形模拟-单晶体"的创建过程类似，在此不再赘述。

后处理界面

应力应变分布：

图5.6 后处理应力应变分布

6. HCP晶格材料的变形模拟-多晶体

所有模块的设置与第2节"FCC晶格材料的变形模拟-多晶体"的创建过程类似，只是需要按照第5节的操作修改HCP晶格材料的材料参数，构建好的inp文件见附件，下面看看模型的变形结果。

后处理界面

应力应变分布：

图6.1 后处理应力应变分布

图6.2 后处理一些SDV结果的分布

7. 多相材料的变形模拟

几何模型

多相材料变形模拟的关键在于材料几何模型的构建和材料参数的赋予。很多研究者为了构建与真实细微观结构相似的几何模型而创建脚本或开发新方法。

图7.1 多相材料的几何模型

我们用前面构建好的圆柱体(R0.015mm&H0.05mm)里包含300个晶粒的几何模型，假设材料里包含最复杂的情况，即FCC、BCC和HCP三相材料，这三相材料的体积分布大致相同，即300个晶粒中每相材料各占100个。

材料模型

材料模型使用前面的FCC、BCC和HCP三相材料参数，如图7.2所示，晶粒编号为1-100设置为FCC材料参数，晶粒编号为101-200设置为BCC材料参数，晶粒编号为201-300设置为HCP材料参数。

图7.2 多相材料的材料模型

所有模块的设置与第2节"FCC晶格材料的变形模拟-多晶体"的创建过程类似，构建好的inp文件见附件，下面看看模型的变形结果。

后处理界面

应力应变分布：

图7.3 后处理应力应变分布

图7.4 后处理一些SDV结果的分布

本文不涉及材料参数应如何获得，材料参数是参考了一些论文的数据还有自己的理解进行的设置，旨在构建一个能顺利模拟的模型。下面给出所有参考文献和在附件给出所有源文件，欢迎交流指正。

8. 参考资料

Ti3Al单晶和双相片层TiAl合金塑性行为的CPFEM模拟

Ti-6Al-4V合金纳米压痕变形与高周疲劳行为CPFEM研究

On Predicting the Channel Die Compression Behavior of HCP Magnesium AM30

两相钛合金拉伸力学行为的研究

密排六方金属镁的晶体塑性力学性能研究

HCP多晶体塑性的数值模拟

TA15钛合金高温变形多晶体塑性有限元模拟

γ-TiAl多晶体压缩变形机制的晶体塑性有限元研究

纯钛单道次ECAP变形织构演化的细观有限元模拟

纯钛晶体塑性力学性能研究

纯钛塑性变形行为的晶体塑性有限元模拟

纯钛压缩变形下的晶体塑性有限元分析

考虑滑移与孪晶的镁塑性本构研究

钛合金双态组织高温拉伸行为的晶体塑性有限元研究