LS-Dyna无网格伽辽金法（EFG）模拟金属切削

EFG无网格伽辽金法模拟金属切削。效果上很像ANSYS的自适应网格重划分。

先放两张炫酷的动图在这里，吸引（诱骗）无知的读者们读下去~

这篇文章和之前许多文章一样，也是一篇学习笔记。

所谓学习笔记的意思是，在写这篇文章之前，我和各位读者一样，也是对这个主题没有任何了解的。开局一个项目文件，剩下的东西全靠自学。

我一边学、一边写，争取在自己弄明白的同时，尽可能地把我的学习过程以及我的理解展示给读者，也给读者讲明白。

其他的一些闲聊放在文末。我们正文对着这个案例直接开搞。

下载和运行

这是LS-Dyna Examples网站上提供的官方案例。ANSYS新版本自带了完整的LS-Dyna模块，下载过输入文件以后，可以直接打开LS-Run来运行。

（这个案例的输入文件，以及后文提到的所有相关文档，在文末都会提供下载方式）

这个案例用的是LS-Dyna的隐式分析。所以注意红圈里面要选双精度。

求解结束以后，在LS-Run窗口点击LSPP图标，就能打开结果文件进行后处理。文章开头那俩炫酷的动图就出来了。

好。其实看这个模型本身也非常简单，就是一个刚体壳，往下压一个长方体金属块嘛。我们就来仔细看看它的输入文件，有哪些特殊之处。

（备注：后面的介绍实在很长，似乎也没办法写的更有条理一些。对EFG技术感兴趣的读者可以直接拉到末尾看几页扫描。）

关键字介绍

稍微玩过一点LS-Dyna（或者读过我前面写的LS-Dyna小球冲击案例）的读者应该知道，LSPP这东西，它虽然有很多定制的高级功能，但其基本功能实际上就是一个大号的关键字文本编辑器。模型Keyword列出来按字母表顺序排列糊你一脸，其中的逻辑关系全都要用户自己梳理。

这一小节，先总结一下这个案例里面包含哪些关键字。先讲其中与EFG无关的，包括材料本构、边界条件等信息。

模型没什么好说的，下面的长方体尺寸24*10*4. 我们按照材料-截面-边界条件-输出和求解控制 这个逻辑顺序依次来看。

首先是材料。

模型只有两个Part。对应的材料分别是020-RIGID和024-PIECEWISE_LINEAR_PLASTICITY。我已经无力吐槽dyna这个毫无逻辑顺序的材料编号了，毁灭吧。

上方下压的壳工件使用刚体材料：

刚体材料定义里，如果CMO（Center of Mass constraint Option，质心约束选项）选择了 +1.0，则可以在材料中直接指定对刚体的约束。或者也可以用SPC来约束它。这里在后面的CON1和CON2关键字中，分别指定了4和7，即约束x, y方向平动自由度，x, y, z方向转动自由度。

下方被压扁的部分，材料名官方写的好像是个德语，Material3 _ Band _ plastisch。用的是分段塑性本构，屈服应力400MPa。这里FAIL关键字对应的数值即破坏时的塑性应变，达到此值时单元被删除。

材料中引用了LCSS，即Load Curve Stress Strain。这条应力-应变曲线由*DEFINE_CURVE关键字定义，长这样。曲线的x轴为等效塑性应变，y轴为应力。所以从屈服极限400 MPa开始，应力1500 MPa时等效塑性应变为12.0。

然后是截面。

刚体的壳单元截面没什么特别的。厚度为1。

下方变形体的截面定义，使用了关键字 *SECTION_SOLID_EFG。这是这个案例的关键，我们在后文详细介绍。

接着是约束。

这里有一个*CONSTRAINED_GLOBAL 关键字卡片，平面法线沿x轴，坐标点(3,5,5)，即 将x=3的平面上的所有点，平动/转动自由度约束为零。

另一个全局约束，平面法线沿z轴。实际上效果就是让我们动图里看到的下面那个工件，一端固定，底部无摩擦（只约束z向位移）。

（dyna里面这个约束关键字的定义方式还真是……保守且费劲）

边界条件上，对刚体施加的边界条件有专用关键字 *BOUNDARY_PRESCRIBED_MOTION_RIGID。引用了一条载荷曲线，定义上方刚体有向下的速度。

被引用的这条载荷曲线，使用了关键字 *DEFINE_CURVE_SMOOTH。这一关键字适合用于金属成型等问题。

接触定义方面，这也是一个专用关键字：

*CONTACT_FORMING_SURFACE_TO_SURFACE_MORTAR

dyna中的*CONTACT关键字，后面跟着的第一个{OPTION}就是接触类型标记。有专用于气囊、腐蚀、金属成型等许多领域的关键字。

在LS-Dyna中，对于Forming类型，专用于金属成型领域的接触关键字，其刚体主面厚度默认会被忽略。

接触关键字的后缀 _MORTAR，英文译为砂浆;迫击炮;研钵，我不是很理解为什么要这么称呼。简单来说，它能让接触对的接触力结果变得更准确，尤其是当接触单元为高阶单元时。

关于输出定义，在LS-Dyna中是用 *DATABASE 系列关键字实现的。在LS-Prepost里，*DATABASE系列关键字被集 合在同一个界面中，可以一起勾选。其中，RBDOUT输出刚体部件的位移，RCFORC输出接触对的反力，SLEOUT输出接触界面的能量。

这些具体的关键字……反正背也背不下来，LS-Prepost界面上都有提示。*DATABASE关键字记住就行，这个不难记。

此外，*DATABASE_BINARY_D3PLOT 指定输出的D3PLOT文件中的时间步长，这里定义为2e-4. 至于其他的关键字……依次讲一遍没啥意思，这东西的鼠标操作基本上和敲一遍关键字也差不多，感兴趣的读者下载源文件自己看就好。

定义参数

和APDL类似，LS-Dyna输入文件中也可以定义一些命名参数。参数名称的第一个字母表示它是实数、整数还是字符型参数。

这些参数定义需要写在输入文件的最前方，引用时使用 & 来引用参数名。例如这样：

这里的tend就是参数，值为3e-2.

在LSPP界面，输入框左上角勾选 Use *Parameter，可以输入参数。取消勾选后，参数会变成数值显示。

最后来看看*CONTROL关键字。

这个模型里，使用到了很多*CONTROL系列关键字。在LS-Dyna中，CONTROL系列关键字被用于定义求解设置。其中，*CONTROL_TERMINATION基本上是必不可少的，用于定义求解的结束时间。

真正输入文件中，所有关键字的顺序都是任意的。在LSPP里它们被按照字母表顺序排列。但我还是按我的认知顺序来讲吧。

首先是刚才提到的*CONTROL_TERMINATION。 在这里定义结束时间，此处定义0.03s. 这里还包含ENGMAS选项，可以和质量缩放选项共同作用。在LS-Dyna中，用于控制时间步长和质量缩放的选项叫*CONTROL_TIMESTEP，很好记，但这个案例中不涉及。

之所以不涉及质量缩放，是因为它根本就不是显式分析。在LS-Dyna中，本身没有分析步这一概念（反正没有叫STEP的关键字）。取而代之的是，如果希望定义隐式分析步，就使用*CONTROL_IMPLICIT。

LS-Dyna中引入隐式分析步有两种方式，一种是*CONTROL_IMPLICIT，另一种是使用 *INTERFACE_SPRINGBACK_SEAMLESS。顾名思义，第二种关键字会在某一时刻以后从显式分析无缝切换到隐式分析，用于方便模拟金属加工后的回弹。

在隐式分析步控制中，此案例包含了4个关键字。第一个是 *CONTROL_IMPLICIT_AUTO。这里可以指定自动时间步长、每一步的建议迭代次数、最大和最小时间步等信息。

*CONTROL_IMPLICIT_DYNAMICS，可以指定隐式动态分析中的时间积分算法。包括静态分析、Newmark时间积分、模态叠加法等。

*CONTROL_IMPLICIT_GENERAL，其实这个案例中不太需要。这个关键字主要用于控制是纯隐式分析，还是显式-隐式联合分析。

真正比较硬核的是这个 *CONTROL_IMPLICIT_SOLUTION 关键字。它里面有非常多、非常细致的求解选项。

这个案例中，许多选项并没有使用默认值。但这些选项……说实话，笔者个人感觉是你就算都读了它的介绍，也大致弄懂了每个变量的含义，实际上求解的时候也需要靠不停试错来获得最优参数。它们有数值求解上的意义，但不一定有特别直观的物理意义。

举例来说，第一个参数NSOLVR，可选线性还是非线性求解器。其实真正有意义的非线性求解器选项就只有一个。剩下的6\7\8\9 四个选项，可能是由于历史兼容性的原因被保留，dyna在文档中提示你别用它们。

再比如第二个选项 ILIMIT，代表在自动重新组成刚度矩阵之前需要经历的迭代次数。这个参数其实很好理解，如果设为1，那每次迭代之后都重建刚度阵，重建刚度矩阵的开销就比较大；而如果设置的数值较大，重建刚度矩阵的开销是降低了，但迭代次数会随之增加。

你即使理解了它的含义，其实也只能按照经验填一个数进去。这里官方填了4，我觉得我改成3或者5，对结果应该也影响不大。如果结果出现了显著变化，那肯定是因为还有其他地方写错了。

……这里其他选项我就不再一一介绍了。实在有点累。不过如果读者感兴趣的话，做类似分析时倒是完全可以参考这个官方案例中的参数设置。

*CONTROL_OUTPUT，看着选项不少，其实大部分都是默认值。不赘述。

*CONTROL_ENERGY，设置能量耗散选项。这里全选了2，即计算所有能量。比如第一个选项 HGEN，就是HourGlass ENergy，沙漏能量，选2是计算沙漏能量并将其包含进能量平衡中。

*CONTROL_ACCURACY，控制二阶应力精度的。（在这个案例里没啥太大作用。）

*CONTROL_RIGID，与刚体部件有关的一些求解设置。这里没啥用。

*CONTROL_SOLID，与实体单元求解有关的一些设置。基本是默认值，没啥用。

*CONTROL_HOURGLASS，减缩积分单元的沙漏控制技术。这里选了6。

*CONTROL_MPP_IO_NODUMP，和输出重启动文件有关的一个关键字。

*CONTROL_CONTACT，控制接触计算的相关参数。选项很多的一个关键字。例如第一个SLSFAC，滑动摩擦罚系数缩放，就被设置为2.5。其他大多为默认。

另外几个关于网格自适应重划分的关键字，后面单独讨论。

自适应关键字

与EFG有关的关键字，除了前面提到的*SECTION_SOLID_EFG以外，*CONTROL系列中还有两个，分别是 *CONTROL_ADAPTIVE 和 *CONTROL_REMESHING_EFG。

先来看 *CONTROL_ADAPTIVE。事实上，自适应网格重划分不是EFG技术独有的，这个技术本身也可被用于非EFG列式的单元。在LS-Dyna选项中，*CONTROL_ADAPTIVE 关键字包含下面这些设置项：

第一个，ADPFREQ，网格重划分的频率，单位是时间。关于这个ADPFREQ，它和下面的ADPASS选项有关。如果ADPASS设置为two pass adaptivity，流程如下图。

如果ADPASS=1，即one pass adaptivity，则流程简化一些：

第一行这里，ADPTYP=7，即3D r-adaptive 网格重划分算法。这种算法只能使用四面体单元，重划分时，基于旧网格使用最小二乘近似生成新的四面体网格。根据最新的关键字手册，使用这一选项时，ADPSIZE参数无效。网格重划分的尺寸依据 *CONTROL_REMESHING 关键字里的设置来决定。

第二行的ADPENE选项，当它的数值大于0时，对于三维实体单元，接触从面的网格细化密度与接触主面（即“模具”）的局部曲率有关。

其它选项与ADPTYP=7无关，都是默认值。

再看 *CONTROL_REMESHING_EFG

这里，第一行输入文件中只设定了两个值，RMIN和RMAX。分别代表需要进行重划分的单元最短边与最长边阈值。这两个值在输入文件中以参数形式体现。笔者猜测这些数值在调试过程中也被修改了不止一次，才找到较合适的参数值。

第二行有一个IVT=1，不过这一选项也是默认值。文档里介绍它有-3的选项，但实际可选的关键字只有±1和±2.

关于EFG

上面提到的这些，好像都和EFG，无网格伽辽金法，没啥关系。即使是网格重划分算法，似乎它也是有网格的，没看出来哪Element Free了呀。

LS-Dyna的文档里关于这部分理论写的不是很清楚。我在网上找了一下，《有限元资讯》2012~2013年的期刊里对EFG功能做了一些介绍。可惜年代久远，这些期刊在知网上都搜不到，ANSYS官网翻了一圈也没找着。在辛春亮等人编写的《LS-DYNA有限元 建模、分析和优化设计》书中找到了引用的相关介绍。我还是直接扫描出来吧。（只是这几页的话，大概应该不算侵权？）

总结来说，EFG单元列式本身是不依赖网格的，它的近似函数基于节点构造。但它仍然需要背景网格，因此在处理大变形问题时，还是需要配合网格重划分算法来使用。

除此以外，我还找到了一篇LSTC官方的介绍。看署名应该也是华人工程师，Yong Guo。链接放在这里。

EFG与SPH方法概述

相关文档，包括 1）上面链接对应的PDF；2）EFG相关的几页中文书扫描；3）Belytschko先生于1994年首次提出EFG方法的论文原文；4）这个案例对应的输入.k文件，我打包放进了百度网盘。

其实这些东西也都是网络上可以公开获取的信息。对此主题感兴趣的读者很容易自己搜索下载到。

后记-碎碎念

好久没写学习笔记了。这一篇学习笔记我写的其实也挺不满意。不过实在没办法，LS-DYNA这东西，它输入文件结构就这么琐碎，根本谈不上GUI操作。再怎么讲，最后也只能依次介绍它输入文件里面的每一条关键字。

再感慨一次，LS-DYNA真是高冷的软件啊。

翻阅《有限元资讯》里面的相关介绍文章，感觉这些世界顶级有限元仿真软件的开发团队中，也有不少华人工程师的身影。比如介绍EFG那篇文章的这两位作者：

我没有任何立场去评论这两位大佬\前辈的选择。在那个年代，如果不出国去LSTC，他们在国内不可能获得这样的平台，不可能有这样的条件去从事全球最尖端却同时又最实际的计算固体力学研究。

前些天的新闻，安世亚太被美国列入了实体清单，原因大致是用了美国技术为中国军方提供技术支持。中美博弈愈演愈烈，希望在我有生之年，真的能看到中国自主的仿真软件，不再只是苦苦追赶，而有机会、有经济和人才条件真正落地做一些尖端的、激动人心的新技术。