有限元中各种锁死的现象及其原因

有限元中各种锁死的现象及其原因

锁死是用来描述结构或单元的刚度过度估计的术语，锁死可以造成结构响应预测的较大偏差。线性插值的实体单元特别容易有多种锁死现象。必须采取特殊的措施来解决锁死问题，下面介绍几种常见的锁死现象。

体积锁死

单元的体积锁死通常发生在处理不可压材料（在超弹性橡胶或金属塑性中）时，即泊松比接近0.5时。当材料接近不可压，此时体积模量接近于无穷，因此将导致奇异的刚度阵。这种锁死不能简单地通过细化网格解决。

弹性应变能可以分解为变形部分和体积部分：

对于不可压或接近不可压材料，当泊松比接近0.5时，体积模量K会变得非常大。对于线弹性各向同性材料，其应力为：

可以看出当v→0.5时，v/(1-2v)会趋近于无穷大，导致刚度阵趋近于无穷大，因此与真实的不可压材料相比，如果应变张量的迹不为零，对应的单元会过刚。

通常，选择性减缩积分，Galerkin混合公式和稳定方法可以缓解体积锁死，其中选择性减缩积分的一种也称为B-bar方法，在ABAQUS的C3D8和ANSYS的SOLID185中都采取B-bar修正作为八节点线性插值实体单元的默认公式。在某些情况下，通过定义节点体积并计算平均节点压力，从而替代整体单元体积的方法也可以缓解体积锁死，也称为Mixed u-p公式，在ANSYS的SOLID185中可以选择该公式。

泊松/厚度锁死

这种锁死现象发生在使用实体单元分析面外弯曲分析中。出现厚度锁死的原因是位移场的近似方式。线性的位移场近似方式会导致沿厚向的法向应变为常数，由于泊松效应，为常数的法向应变会与线性变化的面内应变相互耦合，造成横向剪切应变和面内应变的不一致，导致单元过刚。比如，对于一阶线性插值8节点实体单元，如果其受纯弯曲的话，则x和y方向的应变分布为：

而根据应力-应变关系，对应的应力分量计算为：

其中，

由上式可以看出厚向应变不随y变化，而沿y方向却是线性分布的，这种现象违背了材料本构关系。这种在厚度方向上的应力和应变分布的不协调导致了锁死，即泊松锁死。通过修正弹性本构矩阵，使厚向应力沿厚度方向均匀分布，可以缓解这种情况。

剪切锁死

    剪切锁死通常发生在8节点实体单元和双线性Mindlin板单元中。这种锁死通常在计算板受纯弯曲效应的问题中发生，此时面外剪切应力为零。然而，由于面外位移的采用线性插值，当单元受纯弯曲变形时，横向剪切应变不能在单元的所有节点处都消失，导致单元过刚。

考虑受纯弯曲的细长梁：

其位移的解析解为：

对应的应变分量为：

从上述式子可以看出，此时的剪切应变为零，与梁受纯弯曲变形相互协调。

从有限元角度，位移场的经典插值函数为：

其中的8个参数ai和bi代表了8种不同的变形模式，如下表所示：

考虑纯弯曲变形的条件，只有第4种变形模式是激活的，因此其位移场退化为：

对应的应变场分布为：

很显然，此时剪切应变不满足纯弯曲变形时为零的条件，使得单元变得过刚。

其深层原因是，挠度和转动采用同阶的插值表达式，所以剪切应变中的和两项是不同阶的。以2节点单元为例，由和的插值表达式，将得到：

当梁的高度越来越小时，通过罚函数迫使约束条件得以实现，隐含着上式右端的常数项和一次项系数分别等于零，即：

第一个式子表示单元内梁的平均转动应等于梁法向位移的变化率，也即表示梁中点的直线法假设，在计算中是可以实现的。

第二个式子将导致theta1=theta2，即在单元内theta=常数，这将意味着梁不能发生弯曲，因此问题只能是零解。

这是由于约束条件未能精确满足，在梁很薄时导致不恰当地夸张了剪切应变能项的量级而造成的。

避免剪切锁死的基本点是在计算剪切应变时，使dw/dx和theta项预先就保持同阶。

具体分为以下两种方案：

1. 减缩积分

所谓减缩积分就是数值积分采用比精确积分要求少的积分点数。以Timoshenko梁单元为例，为精确积分剪切应变能项需要采用2点积分。减缩积分方案是采用一点积分，这样一来theta项就不能被精确积分，实际上是以该积分点（一点积分在单元中点）theta的数值代替了在单元内的线性变化，从而使它和dw/dx保持同阶，因此使约束条件dw/dx=theta有可能到处满足，也即使Ks保持奇异性。

2. 假设剪切应变法，即Bathe和Dvorkin提出的Assumed Strain Method(ASM)方法。

梯形/曲率锁死

当实体单元为梯形形状时，即单元沿厚度方向的边不垂直于中性面时，会发生梯形锁死。当用实体建模曲面结构时，梯形的实体单元通常是不可避免的，这会在厚度方向产生非协调的法向应变，从而造成锁死。然而，梯形锁死只发生在单元包含厚向应变的时候，换句话说，退化壳单元和平面板单元由于假设厚度方向的法向变形为零，则不会受到梯形锁死的影响。这种现象首先由MacNeal指出，考虑下列的梯形单元，其坐标在参考平面x-z中可以表示为：

考虑纯弯曲变形的梁中M/EI=1，v=0，则梁的位移的解析解为：

或者使用等参坐标表示为：

则有限元的位移为：

对应的应变场为：

而实际的解析应变场为：

比较解析解和有限元解可知，当时，二者不相等。当变得很大时，这些寄生项会出现，从而导致梯形锁死。

因此，当用2D和3D单元建模曲面结构时，会不可避免地出现单元有斜边的情况，此时就会出现梯形锁死。

膜锁死

膜锁死描述了曲面单元不能精确描述纯弯曲变形的行为。这种锁死现象在线性单元中很少发生，而二阶单元会有很强的膜锁死，膜锁死问题会随着厚度的减小而变得越发严重。

膜锁死现象经常出现在用梁和壳单元计算高曲率结构受弯曲载荷的情况。膜锁死经常与梯形锁死或横向剪切锁死混淆，因为后者影响了膜效应，然而它们是完全不同的锁死模式。

为了更好地说明膜锁死现象，考虑一个长度为2l，半径为R的曲梁，环向位移为u，径向位移为w。沿梁中线的曲线横坐标为s，其膜和弯曲变形分别为：

从上述式子可以看出，这要求环向位移至少为C0连续，径向位移至少为C2连续，则位移表达为：

则膜和弯曲变形为：

当该单元用于模拟纯弯曲变形时（无膜拉伸）膜应变必须为零，则只能满足下列条件：

物理上满足第一个条件可以，而第二个条件如果满足，则径向位移w始终为常数，这与方程假设不协调，所以系数bi是膜锁死的原因。

避免各种锁死的方法

锁死现象通常发生在低阶全积分单元中。因为单元存在一些非物理（虚假的）变形模式，从而导致单元过刚。对于实体单元和壳单元有非常多的方法和技巧可以避免锁死现象。

减缩积分（Reduced Integration,RI）可以缓解某些锁死现象，但不能缓解所有的锁死方式，因此对壳单元需要额外处理来增加性能。

解决锁死问题的最有效方式是采用混合公式，对位移和应变/应力独立插值，或者通过节点位移的投影最大可能地消除虚假变形模式。对于壳单元，通常采用Bathe和Dvorkin提出的Assumed Strain Method(ASM)方法，该方法也可以扩展到实体壳单元。通过使用ASM和RI混合的方式，可以很好地消除剪切锁死和体积锁死现象。

由Simo和Rifai提出的Enhanced Assumed Strain(EAS)方法被广泛应用于低阶有限单元中，该方法通过包含额外的变形模式来消除锁死问题。个人觉得Wilson首先提出的非协调元（也称为Bubble项/泡状函数项，ABAQUS中的C3D8I单元采用了非协调元公式，ANSYS的SOLID185也有该公式）可以看做是EAS方法提出的一个灵感，而Simo的EAS方法更具有普遍性。EAS方法通常与Assumed Natural Strain（ANS）方法一起使用来消除绝大多数锁死现象，经常应用于实体、实体壳和壳单元中。

参考文献：

[1] 王勖成. 有限单元法. 清华大学出版社.

[2] Wang P. Solid–shell finite elements for quasi-static and dynamic analysis of 3D thin structures: application to sheet metal forming processes[D]. Paris, ENSAM, 2017.

[3] Dvorkin E N, Bathe K J. A continuum mechanics based four‐node shell element for general non‐linear analysis[J]. Engineering computations, 1984, 1(1): 77-88.

[4] Simo J C, Rifai M S. A class of mixed assumed strain methods and the method of incompatible modes[J]. International journal for numerical methods in engineering, 1990, 29(8): 1595-1638.