力学之：关于刚度与模量！
来源：WELSIM（ID：welsim666），作者：知乎账号作者ZifengYuan。
我们最早接触到刚度这个概念，应该是在学习弹簧这个物理构件的时候，弹簧的刚度k 定义于弹簧受力F 与其变形δ 的比值。在简单的情形下，我们认为弹簧刚度k 是一个常量，而杨氏模量E 则最早出现在材料力学当中，定义为应力与应变的比值（一维情况下）。
更广义的来说，刚度可以推广为某种结构针对某种变形的抵抗能力，而且不同的定义下刚度的量纲也会发生相应的改变。譬如在梁理论里面，我们可以定义梁的抗弯强度，可以定义杆件的抗扭强度（注意这个量纲是力矩/弧度）等等。在有限元或者结构力学里面，我们所组装的整体刚度矩阵的某个元素的含义，就是某结点在某个自由度方向下移动一个单位长度所需要的力。
那我们来研究一下这个问题。比如，一个长度为L，模量为E，横截面积为A 的杆, 其刚度为多少？
很简单，我们只需要计算在变形δ 下，所需要的力F 的值。
先求应变
然后求得应力
然后施加在杆上的力就可以求得
从而，刚度k 就可以得到
到此我们得到一个简单的结论，在材料力学基础上定义的杆，可以等价成一个类似弹簧的东西。也就是我们把应力-应变关系，转变成了力-位移关系。
那么，我们开始想这样一个问题。假设我们有一个长度L 的一维杆, 这个杆由两根长度为L/2的杆粘结而成，两个杆材料 (E, A) 相同，粘结处可承受的最大应力（也就是强度）为σ 0。我们应该如何建模来模拟这个一维杆在受轴向拉伸时候的力学响应呢？
理想情况下，粘结位置的厚度为0。厚度为0 的一种东西有一个及其麻烦的地方——没办法计算应变。应变是一种变形比例的量度，厚度为0 即意味着比例是不存在的。那么，既然应力-应变关系不可行，就思考力-位移关系。这里，我们可以通过粘结处的强度计算这个位置可最大承受的拉力
由于这是一维杆，这也是整个杆能承受的最大拉力。
此外，我们还感兴趣的是能量，也就是需要做多少功才能把这个粘结处拉断。通常，我们会定义断裂韧性GI （fracture toughness，牛/米的量纲）来度量单位面积断裂所需要做的功。这样，我们就可以用一个简单的线性模型来模拟杆的力位移曲线：
初始阶段为线弹性，直到拉力达到最大值；
当拉力达到最大值之后，随着变形的增大，拉力减小，也就是进入软化阶段；
拉力降低到0，整个力-位移曲线的面积为GI A。
问题到这里似乎解决的挺好的了，但在有限元里面我们应该如何处理这个厚度为0的一个准弹簧的单元呢？
在实体单元里面，我们采用Gaussian积分点，在积分点上我们计算应变（包括增量应变），然后返回应力和模量；再在单元内通过对积分点信息的积分，计算单元内力以及单元刚度矩阵。但这个并不适用于这个问题，我们希望的是这个单元有这样一个性质：
单元允许在某个维度方向上厚度为0；
单元积分点的传入信息是在这个不连续面上的位移的跳跃值 (displacement jump)，而传出的信息是在当前这个跳跃值附近的刚度。
这种单元也就是所谓的cohesive zone element，相对应的材料则是cohesive zone model。这种单元的特性如下：
建模和划分网格的时候采用double-node技术，也就是将cohesive surface上的结点做一个简单的复 制，然后链接成拥有初始厚度为0的cohesive zone element （注意这种单元其实并不一定要求初始厚度为0）；
单元采用Newton-Cotes积分方法（至于为什么不用Gaussian积分方法，大家可以做做研究，因为这个原因至今还未完全解决）；
单元传入的是位移的跳跃值，返回的是这个积分点上的单元刚度矩阵。
为了方便，cohesive zone model会定义一个弹性刚度范围，也就是当位移跳跃小于某个值的时候，我们认为这个单元还是处于弹性阶段。但这个假设会从某种角度产生一定的误差，例如从之前的一维杆说起，实体部分的刚度为
cohesive element 的弹性阶段刚度为
那么，我们可以通过弹簧的串联公式计算整个杆的等效刚度为
不过我们也可以发现，当充分大的时候，不等式两端的误差会充分小。因此，我们在用cohesive zone element 的时候，会将其弹性刚度设置的充分大，让其不会去影响整个结构在弹性阶段的响应。其本质还是应力-应变关系和力-位移关系的互相转化时候的概念问题，简单来说就是一个厚度为0应变不存在的单元，就不应该对结构的线性响应产生任何贡献。
