基于无网格（mesh-free）策略实现单积分点几何必须为错（GND）的计算

参考文献：《Physically based crystal plasticity FEM including geometrically necessary dislocations: Numerical implementation and applications in micro-forming》

GND 演化方程依赖依赖于剪切应变率的梯度或者塑性变形梯度的旋度，而标准FEM/VUMAT 只告诉你每个积分点本身的 γ˙a、Fp，不会直接给梯度。以往广泛应用的数值方案通常是：先把 积分点的数据外推到节点，再用线性形函数求梯度，然而这类方案只能用特定单元（如 C3D8），对自适应网格、复杂接触不友好。

该文章提出的一个mesh-free的方案，该方案的主要优势是不改单元、不加 DOF，只在材料子程序内部，用邻近积分点的数据做一次局部重构，就算出梯度，该策略对某个积分点 x，附近有一团“邻居积分点” xI，作者把它们当成 mesh-free 的“节点”，对每个场变量 u(x)（可以是 γ˙a，Fp 的分量）做 MLS 拟合，如下图所示：

权函数使用立方样条，有紧支撑，距离越近权越大：

在实现上作者提到，立方支撑三个方向尺寸约为5个单元尺寸，最多取最邻近60个（3D）或者30个（2D）积分点，作者指出：当邻域尺寸比网格尺寸还小的时候，这个非局域模型就自然退化为传统的局域模型。也就是说，邻域尺寸本身就扮演了“材料内在长度标度”的角色。

为了提高计算效率秒作者使用了一个“时间滞后 + 公共块”的策略对GND进行更新。

作者使用的方案对于显示大变形分析计算效率非常高，使用标准的C3D8R单积分点即可正常运行，并将所提出的数值模型应用于铜箔拉伸和杯冲过程中的尺寸效应分析，模拟效果如下：

作者的研究证明：通过 MLS 在 VUMAT 里计算 GND，可以在 ABAQUS 中完整重现微成形的尺寸效应，并清晰揭示“GND 在晶界和局部剪切带聚集”是强化的主要来源，同时保证数值方法可扩展、可工程化。详细的数值实现策略可以参考原始文献。

使用文章提到的策略，尝试进行数值显示，首先在umat隐式中进行实现，并在后续中修改为vumat即可。实现策略验证使用包含200个晶粒的二维模型拉伸验证。分别使用CPE6单元（二维多积分点，使用传统的GND计算方案），CPE3单元（mesh-free策略），模型共包含27119，SSD计算使用经典的KM模型，流动方程使用唯象的幂律模型，取向随机分配给不同晶粒

初始多晶模型和网格如下：

拉伸变形10%后应力分布：

传统方案：

MLS方案：

拉伸变形10%后累计剪切分布：

传统方案：

MLS方案：

拉伸变形10%后SSD分布：

传统方案：

MLS方案：

拉伸变形10%后GND分布：

传统方案：

MLS方案：

总的计算时间

传统方案：60分钟15秒

MLS方案：36分钟12秒

可以看到MLS方案的计算策略计算效率通常会显著高于传统的GND计算策略，并且数值实现显然更适合在vumat框架中（计算效率高，单元类型适应性强），而且可以抑制局部网格噪声，因此在大变形，接触分析，成型计算中是一种非常合适的选择。感兴趣的可以进行数值尝试，也可以加入知识星球讨论交流