翻翻帮助文档，学习使用Abaqus的单元删除功能模拟材料的断裂和损伤。
我们来学习一下，Abaqus等有限元软件中，是如何用单元删除的方式来模拟断裂和损伤行为的。
所以……你问我有限元该怎么学？带着给别人讲出来的心态去学。边学边讲。能给别人讲明白，就说明你这部分知识已经学会了。
话不多说，我们开始。
找到帮助文档页
用单元生死的方式来模拟断裂，首先能够想到的肯定是Abaqus软件。
因为在Abaqus的单元选择对话框里，赫然有四个大字：“单元删除”。
Abaqus的单元类型选择界面，支持单元删除功能。
实际上，Abaqus/CAE单元类型选择窗口的这些单元控制属性，在关键字里其实被称为 Section control 截面控制。这挺让人费解的。但至少 还能在帮助文档里查到。
首先在界面上点击那个小问号，点击单元类型窗口。会打开 Abaqus/CAE 图形界面的帮助文档。定位在
Abaqus > Abaqus/CAE > Creating and analyzing a model using the Abaqus/CAE modules > The Mesh module > Controlling mesh characteristics > Associating Abaqus elements with mesh regions
好的，跟着超链接，来到了 Section controls 页面。
Abaqus > Elements > Introduction > Section controls
这一页内容可丰富了。有这么多主题：
那既然软件界面上都摆出来 单元删除功能了，文档里肯定会有对应功能的介绍嘛。
这个主题下 —— 还是没啥真正有用的信息。只告诉了你，界面上的 “单元删除” 和 “最大下降” 选项在哪。它俩分别用来控制 1）“单元删除”选项 决定在计算中是否要保留完全损坏的单元；2）“最大下降”选项 确定单元何时被认为是完全破坏了。
——摔，“最大下降”，这翻译的谁能猜到是啥意思嘛。Max degradation，准确来说应该翻译为“最大刚度退化值”。
好在，这一页里，给出了更多的详情页链接。来列举一下：
1）在Abaqus 帮助文档的材料卷，介绍渐进损伤和失效的页面：
Abaqus > Materials > Progressive Damage and Failure > About progressive damage and failure
2）材料卷的介绍部分中，有一小节是关于材料失效与单元删除的：
Abaqus > Materials > Introduction > Material data definition > Material failure and element deletion
3）还是材料卷，延性金属的损伤与失效：
Abaqus > Materials > Progressive Damage and Failure > Damage and failure for ductile metals > Damage evolution and element removal for ductile metals
4）以及在Connector单元行为的文档里：
Abaqus > Elements > Connector Elements > Connector element behavior > Connector damage behavior > Maximum degradation and choice of element removal in Abaqus/Standard
延性金属与损伤演化
在翻译之前，先聊一句题外话。ductile metal，中文翻译为延性金属。我在Science Direct上面搜索了一下dectile metal的定义，还是很全面、很有趣的。
一些脆性材料，如铸铁和混凝土，在其结构中含有大量的孔和微观裂缝。这些被认为会产生高应力集中，从而导致材料的局部失效。相对于抗压强度来说，这些应力集中对抗拉强度的影响可能更大;脆性材料的一个一般特征是它们在受拉时显得很脆弱。因此，脆性材料的弹性断裂和失效主要由最大主拉应力控制。作为一个例子，考虑一个混凝土块：在单轴拉伸下，其断裂应力约为1.5 MPa，而在压缩时，它约为30 MPa，或20倍大；在纯剪切下，断裂应力约为1.5 MPa，因为主应力具有相同的大小，其中一个应力是拉伸。混凝土中的开裂将发生在与施加的剪切应力方向成45°的平面上。
（天呐，我学材料力学 第一强度理论的时候怎么没人这么给我讲过？该不会是老师自己也不知道吧。唉，我要是有机会能给学生讲课，我可得把这一段加上。多满足好奇心呀。）
Abaqus > Elements > Introduction > Section controls > Controlling element deletion and maximum degradation for materials with damage evolution
Abaqus提供了一种通用的功能，用于对材料渐进损伤和失效进行建模。在Abaqus/Standard求解器中，这种功能仅适用于：内聚力单元、连接器单元、平面应力公式单元（包括平面应力单元、壳单元、固体壳单元和膜单元）、可用于韧性金属损伤演化模型的任何单元、以及可以在低周疲劳分析中包含损伤演化准则的任何单元。
在Abaqus/Explicit求解器中，适用的单元类型则更多，所有具有渐进损伤行为的单元都可以使用单元删除功能。
截面控制选项，可用来指定材料失效的最大刚度退化值，D_max，这个值用来描述材料的失效，以及判断啥时候该把这个单元删除掉。默认情况下，单元完全失效时会被删除。单元删除方式的选择也会影响到损伤施加的方式，详见各种详情页面。
渐进损伤和破坏
在Abaqus中，渐进损伤主要应用于下面几类问题：
Abaqus中，提供了一个用于描述损伤和失效的通用框架，它允许多种损伤失效机制同时作用于同一种材料。材料失效，指的就是由于刚度逐渐退化而导致失去承载能力。
比如，一个典型的金属 单轴拉伸试件：
典型的金属单轴拉伸应力-应变关系
a-b，大家都知道是初始的线弹性段；
b-c，是熟悉的塑性屈服和应变硬化阶段；
c-d，超过c点后，承载能力下降，直到材料破坏。
在点c处，激活某一损伤准则。这一点处的状态称为损伤启动准则。超过它以后，应力-应变的响应将改为受局部化的刚度退化演化行为控制。
如果没有损伤行为，材料在c点以后，正常来说应该按照c-d'的路线继续强化。
因此，在Abaqus中，描述一个失效机制就需要输入四个部分：
定义一个有效的材料响应（a-b-c-d'），即没有损伤时是啥样。
定义损伤的起始准则（c点）
定义一个损伤演化的规律（c-d）
损伤完全发生后（d点处），选择是否要将单元删除。
网格依赖性
在连续介质力学中，本构模型通常用应力-应变的关系表示。当材料出现软化行为，导致应变局部化时（即大量的应变出现在局部一个小区域里），这种形式会导致极强的有限元网格依赖性。其结果是能量耗散随着网格的细化而减少。（据我理解，就是单元画的越密，发生应变软化的单元体积就越小，因为只有裂纹尖端的单元会破坏）
在Abaqus中，所有可用的损伤演化模型，都采用了一种旨在减轻网格依赖性的列式。这是通过在公式中引入特征长度来实现的。在Abaqus中，这与单元的大小有关，并且将本构中的软化部分改为用 应力-位移 关系描述。这种情况下，损伤过程消耗的能量就是按单位面积，而不是单位体积来指定的。这一能量作为一个附加材料参数，用于计算材料发生完全损伤时的位移。这个概念和临界能量释放率（断裂力学中的一个材料参数）一致，以确保正确的能量消耗，并大大减轻网格依赖。
Abaqus中的使用方法
前面我们提到，描述失效机制，需要定义四个部分：
前面，定义正常材料响应这里就不提了，就和定义普通的弹-塑性材料一样。至于这个损伤起始的准则，可以在Abaqus定义材料参数的界面里看到：
Abaqus中定义损伤起始准则
这个菜单里，有延性金属损伤、有牵引-分离模型，还有纤维增强复合材料的损伤。再后面的“弹性体损伤(T)”，菜单里只有Mullins效应，它其实指的是类似橡胶的弹性体在准静态循环载荷下的应力软化现象。
定义了任何一个损伤起始准则后，都可以继续在 子选项 中看到 损伤演化。
在子选项中，进一步定义损伤演化准则
损伤演化子选项输入窗口
三点弯曲试验，理论与仿真对比
所以我才说，仿真软件是最好的教科书。
当材料完全损伤以后，可以将相应的单元删除。损伤演化公式使用了一些方法来尽量减少网格依赖性；定义损伤模型，需要定义四部分信息，其中最主要的就是两条：
1. 损伤起始准则；
2. 损伤演化准则。
这部分文档，位于Abaqus文档-材料卷。渐进损伤章节里，延性金属损伤主题下，只有三个部分：1）简介；2）损伤起始准则；3）损伤演化准则。
关于延性金属损伤，对应的帮助文档位置
损伤起始准则
回顾一下：损伤起始准则，描述的是什么时候达到图中的c点。达到c点之前，材料没发生损伤；达到c点之后，损伤开始。至于后面损伤如何演化，c-d曲线部分，那不归起始准则管，而是损伤演化关心的事儿。
单轴拉伸曲线，c-d描述损伤过程.
损伤起始准则主要分为两大类：用于金属断裂和金属薄板失稳。
延性金属损伤的起始准则，Abaqus界面
损伤起始准则的分类
每一个损伤起始准则，都有一个对应的输出变量，用于指示分析过程中是否达到了该准则。1.0或更高的数值表示已经达到对应的损伤准则了。
对于某一材料，可以同时指定好几个损伤准则。它们相互独立，单独计算。一旦达到了某一个损伤准则，材料的刚度就开始按这个准则里写的损伤演化规律开始下降。但你要是没定义损伤演化规律，那材料的刚度就不会下降，只会输出该损伤准则对应的变量数值，告诉你这地方要破坏了。
金属断裂的准则
主要有两个机制会让延性金属发生断裂：
金属内 微孔洞的成核、生长和聚合导致的韧性断裂；
剪切带局部化导致的剪切断裂。
这俩机制具体是啥，有一篇2004年的论文详细描述了这个问题。
Hooputra, H., H. Gese, H. Dell, and H. Werner, “A Comprehensive Failure Model for Crashworthiness Simulation of Aluminium Extrusions,” International Journal of Crashworthiness, vol. 9, no. 5, pp. 449–464, 2004.
作者是来自BMW集团的大佬。真的是，巴伐利亚发动机厂在各种仿真软件的技术文档里出现频率可太高了，不愧是世界顶级车企啊。这篇文章最早是用PAM-Crash实现的这些准则。
另一位作者来自MATFEM公司。这是一家来自德国的材料科学与有限元分析咨询公司。
MATFEM的业务范围。看上去是不是比国内一些小规模的仿真咨询公司高级一点？
我说仿真软件是教科书没错吧。
基于现象学的观察，上面提到的两种机制需要不同的损伤起始准则。所以Abaqus就引入了不同的准则来描述它们。
后面如果看到了公式，别怕，它们都是纸老虎。什么叫现象学 phenomenological啊？就是做了个实验，测出来一条曲线，然后根据实验现象拟合出来一个公式和一堆系数呗。嗯，这很工程。
柔性损伤
这里我用了Abaqus中文界面上的翻译。虽然这个翻译的不太好吧。原文是 Ductile Damage，韧性损伤起始准则。或者沿用我上文的翻译，延性损伤准则，都行。
延性准则是一个唯象模型，来预测材料内部孔洞成核、生长和聚集引起的损伤起始。这个模型假设，损伤开始时的等效塑性应变，是应力三轴度和应变率的函数：
这个应变头顶带一横表示等效；上标pl表示塑性，下标D表示损伤。这都很简单。括号表示等效应变是里面俩参数的函数。里面一个是应力三轴度，另一个是应变率。
应力三轴度，就是静水应力除以Mises等效应力。
分子是静水应力，分母是von-Mises应力。这个很好理解。静水应力是应力张量的第一主不变量，是三个主应力之和。而von-Mises等效应力则定义为：
我们知道，如果材料处于各向同性的受压状态（或受拉状态），那它的von-Mises应力为零，是不会发生屈服的。所以应力三轴度，数值越大，就代表当前点应力状态越接近各向同性静水压状态；应力三轴度的数值越小，表示当前应力越偏离静水压力状态，因此就越危险。
应变率，具体来说是等效塑性应变率。（如果还想深入了解这部分，不如试着读一读我去年翻译的一部分Abaqus理论文档？）
延性金属损伤准则假设，当满足下面条件时，损伤开始发生：
又是一个张牙舞爪的公式。别怕，你看公式里有俩等号呢。所以左边这个 \omega_D 的符号其实就是描述损伤是否发生的状态变量。它随塑性变形增加而单调递增。当它等于1的时候，即表示损伤开始发生。
这个积分，其实分子和分母都是同一个东西，就是等效塑性应变嘛。在Abaqus计算的每一个增量步中，损伤状态变量的增量为：
我觉得，写成这样 基本上就已经很接近软件内部实际的代码实现了。
更详细的讲一下就是，这个积分的分子，是仿真软件求解出的塑性应变增量；而分母则是由用户定义的一个准则。当结构的变形历史使得塑性应变积累到一定程度，以至于塑性应变累加到一起，超过了分母这个准则，损伤就会发生。
而因为在变形历史中，可能不同时刻的应变率不一样，不同时刻的应力状态不同从而应力三轴度也有不同，所以每个增量中，分母的数值可能也是不同的。这样就需要做出积分来描述损伤累积的程度。
再回头来看看软件界面里的输入框：
Abaqus软件界面，延性金属损伤需要输入的参数
看懂了没？上面写的那个等效塑性应变，说它是应力三轴度和应变率的函数，实际上Abaqus也不知道这函数长啥样。这个函数其实是由用户来定义的。
还是开头提到的那篇论文，Abaqus其实就是把论文里的公式实现而已。
熟悉吗？这个公式(7)，就是判断是否发生损伤的准则。
在同一篇论文的下方，作者服务到位，给出公式，也给出了对应的实验结果：
作者还区分了准静态加载和动态加载两种情况。
从画出的图中可以看到，应力三轴度数值越小，这个作为准则的等效塑性应变就越大。前面提到，应力三轴度越大，意味着应力状态越接近三向静水应力状态。所以，这个曲线的含义是，当一个偏应力状态，叠加上一个较大的静水应力状态时，材料内部的微小孔洞更容易发生扩展，从而作为准则的等效塑性应变就更小。这在物理上也是合理的。
关于这部分工作，后面还有一些研究者提出了补充。对铝合金的实验结果表明，除了上面提到的应力三轴度 eta 和应变率以外，延性金属的断裂还依赖于偏应力的第三主不变量。而这与Lode角（注意这个单词是Lode，和加载load不同，应该是个人名）有关。Lode角，又称偏极角，deviatoric polar angle.
那，既然等效塑性应变还和第三个因素有关，就像前面的操作一样，把等效塑性应变定义成这三个量的函数呗。
文档在这，自己读吧。反正……这功能甚至都没做GUI
是的，这功能引用的参考文献对于Abaqus软件来说可太新了，最新的文献到达了2015年。所以——在Abaqus/CAE界面里无法访问这个功能，也合情合理对吧？
Defining dependency of ductile criterion on Lode angle is not supported in Abaqus/CAE.
所以——讲了这么半天，这功能究竟要怎么用？
我猜，写出这篇论文的BMW集团，肯定是知道怎么用的。软件把功能开放给你，因为这功能不值钱，公式都在论文里写着了。但输入框里具体要填哪些数，这就需要在不同的应变率下做相应的实验，才能拟合出一条等效塑性应变与应力三轴度和应变率的关系曲线。
渐进损伤退化过程中的应力-应变曲线
Abaqus中的损伤演化准则。
损伤演化
满足任意一种损伤准则以后，材料就开始发生损伤。从这时开始，损伤演化规律登场。在Abaqus中，使用了一个标量 损伤变量，d_i，来描述刚度退化的程度。因为可以同时定义好几个损伤失效机理，所以这个损伤变量也可以同时有好几个。每一个变量对应一个损伤准则。
分析过程中，任意时刻，材料的应力张量都是根据下面的损伤方程计算：
右边这个头顶上带一横的应力，sigma\bar，代表没有损伤时，应力应该是多少。Abaqus称它为有效应力张量，或未损伤的应力张量。很显然，括号里的 D 就是总损伤变量。当D=1时，该点处的应力降低到零，说明材料已经完全被破坏掉了。在默认情况下，这种时候就要把这个单元删除掉。
总损伤变量 D，将统筹计算用户定义的每一种损伤演化机理造成的影响。
这儿有一个画的稍微复杂一丢丢的图，来展示损伤演化过程中的应力-应变情况：
渐进损伤退化过程中的应力-应变曲线
在D=0那一点处，对应的应力和等效塑性应变，即为损伤开始时的应力和应变。接下来，应变继续增大，但应力却因为损伤的发生而开始减小。如果损伤发生一半进行卸载，则卸载时的弹性模量也会随损伤而变小。
Hillerborg提出，张开单位面积的裂纹所需要的能量 G_f 是一个材料参数。
（哇 看Abaqus文档里的公式用MathJax渲染出来多漂亮）
实在看不懂也没关系。能看出来这个等式最左边的断裂能 G_f 和等式最右边的 等效塑性位移，是一码事，可以相互转换，就行了。所以指定损伤演化的极限时，指定这二者任何一个，都是等价的。
至于公式里出现的一个特征长度，L，这东西则取决于单元的几何形状及单元类型。总之——Abaqus建议你在考虑塑性断裂和单元删除的时候，尽量使用长宽比接近的单元。因为长宽比太大的单元在不同的开裂方向上可能会（因为不同的特征长度）而产生不同的行为。
定义损伤演化方式
如果你就想让材料在刚满足损伤准则的时候马上失效，你当然可以把失效时的等效塑性位移指定成0. 但……最好别这么做，因为这样会发生很严重的收敛性问题，把静态问题变成动态问题。毕竟隐式静力分析里不包含惯性项。
那么，接下来的问题就很简单了。我们知道了什么时候损伤开始（损伤起始准则）；又定义了什么时候损伤彻底完成（位移准则或断裂能）。还剩下的问题就是，从无损伤到完全损伤这段过程怎么走。
更具体来说，就是等效塑性位移 u^pl，和我们前面提到的损伤变量 d 之间的关系。
这里提供了三种方式：1）表格；2）线性；3）指数。
很好理解，就不再展开细说了。如果你选指数，Abaqus会让你多输入一个参数；选表格，会让你自定义 u 和 d 之间的关系表格。
这些具体公式，帮助文档里都有。看着挺复杂，但实际上对于物理意义来说，它们都不重要。反正真想算准，这些参数都要通过实验来反算求出。