晶体塑性每日文章推荐（五）

文章名称：《Comparison of the implicit and explicit finite element methods using crystal plasticity》

doi：10.1016/j.commatsci.2006.08.002

推荐理由：作为显式晶体塑性与隐式晶体塑性模型比较，以及适用性讨论最经典的文章之一，详细介绍了显式与隐式求解器的区别，其研究结果表明：对于更简单的加载条件，隐式方法的求解时间更短。在涉及接触的载荷条件下，显式方法被证明是优选的选择。与使用多个处理器解决分析的隐式方法相比，显式方法显示出持续高水平的并行化效率。

众所周知：在隐式方法中通常使用牛顿迭代方法，解涉及迭代，直到每个增量都满足收敛标准。因此计算收敛于精确解（隐式：基于t+Δt时刻确定t+Δt时刻的状态变量），而显式方法中的有限元方程使用的通常是欧拉向前方法，并且在这种形式下，它们可以直接求解，以在增量结束时确定解，当增量步小于稳定时间增量时，总是可以保证计算的进行，但结果不一定收敛于精确解，（显式：基于t时刻确定t+Δt时刻的状态变量）因此准静态问题使用隐式往往更加高效，而涉及到接触则使用显式可以保证计算稳定

由于晶体塑性模型考虑了介观尺度的塑性变形的真实物理过程，因此被广泛用于已在模拟金属和金属基材料中的介观尺度下大变形和应变局部化，然而其高度非线性的积分过程，相较于传统模型，其数值成本往往很高，因此基于复杂的晶体塑性模型更能体现显式于隐式积分的差异，这里作者的讨论显式程序和隐式程序使用huang的亚弹性框架进行

由于隐式程序在大量的博士论文可以找到详细内容，这里不做赘述，这里主要提到文献关于显式的一些内容。

显式计算存在临界稳定时间的概念，当最大增量步大于该时刻时会导致计算结果发散，从而使得计算结果失去意义。

其中Le是单元特征长度，Cd是膨胀波速，公式中可以看出Le越大，（单元尺寸越大）稳定时间越大，密度愈大，稳定时间越大（这意味着可以通过质量缩放加快计算时间，需要根据能量评估质量缩放的影响），需要注意的是显式时间为真实变形时间，因此使用晶体塑性这类粘塑性本构模型显式求解器并不适用于准静态分析，计算时间太长，并且人为修正时间会影响本构的响应情况，当利用显式求解器进行准静态使用质量缩放是合理的备选方案（保证动能于总能量之比小于5%）

作者开发的显式晶体塑性框图

作者使用的是晶体塑性里面经典的PAN模型，材料为316L，晶体塑性参数为：

C11 = 205 GPa, C12 = 138 GPa, and C44 = 126 GPa，g0 = 50 MPa; g1 = 330 MPa; h0 = 225 MPa;gamma_dot - 0:001，n=20

注意：显式求解过程涉及大量增量，每个增量都很容易求解。因此，显式有限元计算程序非常适合由多个处理器分解和求解。考虑到这一点，VUMAT是用矢量化接口构建的。在每次增量开始时，应力和状态变量数据以二维阵列的形式传入。阵列中的每一列都包含与材料的积分点相关的信息。当使用多个处理器执行模拟时，可以将分析数据拆分为块并独立求解。因此，在子程序的编写中保留了矢量化，以便可以实现最佳的处理器并行化。使用ABAQUS/explicit不需要相同的过程，因为时间增量足够小。为了确定ABAQUS/explicit中初始时间增量的大小，在分析开始时将一组伪的微小应变增量传递给VUMAT。根据材料的应力响应，计算出稳定时间增量的保守值。有限元求解器要求材料在初始检查时具有弹性。由于晶体塑性子程序计算密集，必须定义弹性应力-应变响应，以确保相对较小的时间增量。在这个初始时间增量计算阶段，材料响应被定义为具有与用于描述晶体塑性子程序主体中的弹性相同的弹性特性。这确保了采用相对有效的时间增量。如果需要较小的增量，则可以使用更硬的弹性模量，尽管解决时间会更长。这不会影响材料在分析过程中的响应；它纯粹是为了计算时间增量。

同时需要注意ABAQUS/标准采用Jaumann应力率，而ABAQUS/explicit采用Green–Naghdi应力率

Ω表示变形梯度极分解对应的旋转部分，

在ABAQUS/标准中，材料被视为基于固定的全局坐标系。增量旋转在每次增量开始时传递给UMAT。这个数组dR是应力和应变数组在上一次增量结束和当前增量开始之间旋转的量。这样处理：

其中σt+Δt是当前柯西应力，dσ是在先前增量中计算的应力增量。每次增量开始时传递的值是模型坐标系中的柯西应力。如图所示，它在每次增量结束时都会“旋转”。dR变量使用Hughes–Winget算法计算

其中I是单位矩阵，ΔW是增量速度梯度的反对称部分，即W的增量形式

在VUMAT的情况下，材料位于同旋坐标系上，即坐标系随材料旋转。由于每次增量都很小，所以假设所有材质旋转都是刚体。因此，对于每个增量，问题都可以看作是一个小应变问题。ABAQUS/explicit中使用的旋转是Green–Naghdi自旋速率。Hughes–Winget算法也被使用，但形式略有不同

注意;由于晶体塑性子程序不包含运动硬化，两种方法产生的结果相同。

作者比较了二维和三维典型的变形情况以及接触问题，如下图所示

二维拉伸：(比较了流动应力局部累计剪切和计算时间)

三维：

以及显式和隐式处理器并行问题：

其研究结论为：

由于接触和材料响应引起的非线性的严重性使得所需的时间步长远小于实际可用于用隐式方法生成解的时间步长。显式求解器更适合处理复杂的接触和滑动条件，特别是在大单元变形的情况下。在包括刚体的接触分析中，使用隐式代码的运行时间短于显式代码，同时显式可以处理近似率无关情况，隐式则难以计算，显式对于多核的处理效率的提升显著高于隐式

通过文献的思路可以实现将黄隐式程序通过接口文件转化为显式的vumat子程序，单晶与多晶对应的案例如下（参数均为原始的Cu参数）：

案例一：显式与隐式单晶含圆孔薄板的简单拉伸对比（显式开启了质量缩放）

使用位移边界条件，沿着X方向施加10%工程应变的拉伸变形，其中等效应力，累计剪切应变为对比指标，结果如下：

等效应力分布：

累计剪切应变分布：

案例二：显式与隐式多晶的简单压缩对比（包含30个晶粒）（显式开启了质量缩放）

使用位移边界条件，沿着X方向施加10%工程应变的拉伸变形，其中等效应力，累计剪切应变以及滑移系统的当前强度为对比指标，结果如下：

等效应力分布：

累计剪切应变分布：

滑移系统的当前强度：

可以看到这种基于接口的子程序可以显式与隐式的一致性良好（部分差异在于显式使用了过大的质量缩放）