黄永刚晶体塑性UMAT子程序相关教学
（1）FCC/BCC改成HCP
黄永刚本构模型适用于立方晶体，FCC/BCC，由于三个轴是相互正交的，通过输入三指数来描述，
(hkl)用于晶面，[uvw]用于晶向
在黄永刚的UMAT中可以根据输入的三指数，自动生成所有的滑移方向
和滑移面法向
。但是该代码不适用于密排六方的本构晶体。
对于密排六方晶体（HCP）结构，滑移面和滑移方向通常采用四指数法来描述：
(hkil) 用于晶面， [uvtw] 用于晶向
本教学提供基于输入的四指数自动生成所有滑移方向和滑移面法向的代码，适用于HCP的晶体塑性计算，并提供材料参数说明。
（2）加JC本构
Johnson-Cook本构模型用于描述材料在动态载荷下的本构行为，其核心思想在于材料的流动应力取决于塑性应变、应变率和温度三者耦合作用。可以在黄永刚晶体塑性的UMAT中加入JC本构模型，使该本构模型考虑应变率硬化和温度软化。
修改方法在于修改硬化模型，其硬化模量计算公式如下：
上述模型中前两项为Bassani硬化模型，黄永刚程序里有，第三项为JC本构中的应变率硬化和温度软化项。修改后的本构计算效果如下图所示。
图 { SEQ 图 \* ARABIC |1} 本构模型应力-应变曲线，(a) 单晶压缩模型，(b)不同应变率、温度的应力应变曲线
（3）Voronoi多晶结构进行材料随机取向赋值脚本
晶体塑性本构用于描述单晶的行为，然而大多数材料是由多晶构成的，采用晶体塑性UMATs进行多晶的仿真，可以借助Voronoi多晶模型，即生成具有多晶结构的代表性体积单元、模型。如下图所示。
图 { SEQ 图 \* ARABIC |2} 多晶结构的代表性体积单元、模型
当晶粒较多时，采用人工给不同晶粒逐一设置材料取向信息过程繁琐，可以采用脚本进行取向赋值，脚本支持随机取向或者具有织构的取向。
（4）在滑移的基础上考虑孪晶系
在传统晶体塑性理论中，除了常规滑移变形外，还引入孪晶作为另一种独立的塑性机制。当考虑了孪晶后，模型变为：
这类模型还可以拓展为形状记忆合金相变+塑性的晶体塑性本构。
（5）修改VUMAT
UMAT主要应用于隐式分析，而对于大变形接触问题，隐式分析往往计算效率较低。对于接触、碰撞、冲击等问题采用VUMAT往往具有更高的计算效率和收敛速度。可以将UMATs修改为VUMAT，用于隐式分析。本教学提供修改方法、代码及案例。
（6）添加损伤本构
晶体在滑移的作用下如何产生永久变形。采用黄勇刚晶体塑性本构可以模拟出晶体在恒定载荷下产生的蠕变。由于塑性变形进一步地会出现刚度下降，承载力下降的情况。标准的晶体塑性模型，无法描述这一现象，但如果在此基础上引入连续损伤本构模型，定义损伤因子
，并将连续损伤因子
与切应力关联，可以实现对材料损伤的仿真。其建模思路如下图所示：
图 { SEQ 图 \* ARABIC |3} 晶体塑性蠕变损伤本构建模思路
（7）添加背应力
背应力是描述运动硬化和包辛格效应的关键机制，背应力演化公式：
其中，
是第
个滑移系的背应力变化率。
添加背应力前后，循环曲线如下图所示：
（8）添加蠕变和损伤
蠕变模型公式：
第一项：纯蠕变贡献，第二项，损伤耦合贡献。
添加蠕变前后，计算的位移-时间曲线如下：
另外，该模型也可以用于计算晶体塑性低周疲劳，载荷谱及计算结果如下图所示：
（不考虑蠕变）
（考虑蠕变）
（9）晶体塑性损伤+单元删除切削案例
（10）位错密度晶体塑性本构模型
在黄永刚晶体塑性框架的基础上，进一步开发了位错密度晶体塑性本构子程序（Abaqus UMAT）：位错密度以显式状态变量引入，12 个刃型与 6 个螺型几何必需位错（GND）分量逐增量求解，统计存储位错（SSD）走独立演化通道，硬化律采用位错森林投影而非总位错密度。相比只有标量硬化的唯象模型，它能回答三个问题：晶粒尺寸效应从何而来、织构演化多强、晶界处位错如何累积。
（1）晶粒尺寸效应。三套 Voronoi 多晶模型仅改变晶粒数（平均尺寸 6.76 / 4.35 / 3.48 μm），单轴拉伸至 10% 应变。细晶的流动应力提升 43%，而总位错密度几乎不变（比值 0.960），位错摩擦应力却达 1.393 倍——细晶强化的载体是 GND 构成的位错森林，不是位错总数；应力与 1/d 的平方根呈线性关系，可据此定量估计 Hall-Petch 斜率。
（2）织构演化。拉伸至 30% 应变，取向由每个增量步解出的滑移量积分更新，可直接从 ODB 抽取后绘制极图、计算 ODF。织构指数 J 由 10% 应变时的 +1.8 / +4.8 / +4.6 跃升到 +3.8 / +24.1 / +31.5，76% / 89% / 85% 的晶粒转动超过 5°；初始取向按构造均匀随机，长出来的结构完全是变形诱导的。
（3）晶界位错累积。晶界两侧取向不同、滑移系不同，塑性应变梯度的跳变在晶界处造成 GND 累积。18432 个 C3D8 单元、20 个 Voronoi 晶粒拉伸至 5% 应变：内部晶界的 GND 密度达 170.3 μm-2，深部晶内仅 50.0 μm-2（3.4 倍），而自由表面只高 1.36 倍，说明堆积主要归属晶界而非表面效应；SSD 由滑移量决定，晶界处只有晶内的 0.34 倍。同一网格、同一时刻，两类位错给出完全不同的空间形貌。
上述（1）-（10）的修改可以根据用户的材料特性，科研需求，组合定制本构子程序，可以提供教学+代码+案例，或者单独提供代码+案例。
