主流有限元软件剪切自锁和沙漏模式的比较
笔者将剪切自锁、体积自锁、沙漏模式等归结为单元技术。一般用户不需要关注有限元分析的单元技术，换言之，单元技术由于偏向数学层面而难以被一般用户所掌握。当然无论是资深用户还是初级用户，了解一下这方面的知识总是有益的。
本文内容整理自以下文章。
在有限元分析中，剪切自锁和沙漏是比较突出的２个问题,它们通常给出较大误差甚至完全背离的结果来掩盖真实的求解值，如果不仔细对结果进行分析,带来的后果是致命的，尤其是应用较广的薄壳结构，弯矩变形是最主要变形方式，不合理的剪切自锁和沙漏影响带来完全相反的分析结果。做任何复杂的工程分析前，有必要对剪切自锁和沙漏产生的原因，补救措施以及几个主流 FEM 软件如何解决,相互间的性能和效果对比做一个深入的研究。
正交异性悬臂梁。
梁的尺寸为 L× W× H=150mm× 4mm× 6mm，一端固定，并施加４０ N 的力作用于自由端，材料的属性见下:
E11=146.9GPa；E22=E33=10.89GPa；
G12=G13=10.89GPa；G23=6.4GPa；
V12=v23=0.38；v13=0.776
ρ=1.5×103kg/m^3
由于位移场，应变场和应力场之间的关系比较复杂，单元形函数的不完备性以及单元形状的敏感度，单元刚度会出现几种自锁现象:泊松比自锁，剪切自锁,横向剪切自锁,厚度自锁,梯形自锁,膜自锁。泊松比自锁和不可压材料(高泊松比材料)相关，膜自锁来源于膜单元，梯形自锁和单元形状相关,厚度自锁引起于板壳厚度方向的插值精度,横向剪切自锁存在于横向承受剪切力的板壳单元。但是在这里我们仅讨论比较普遍的剪切自锁。
真实的结构在纯弯矩作用下，上下两边界会变形为弯曲的曲线，下图的水平点线和上下边发生了弯曲，垂直的点线继续保持为直线，使得水平和垂直点线的夹角始终保持为直角。
为了正确的表达结构真实形状的变化，有限元分析中的单元必须保证上下边以及中间的水平点线能发生弯曲变形，但是全积分下的一阶单元由于单元的形函数是一阶，不能描述和表达高阶的曲线形状，变形后边界和所有点线都保持直线形式见下图，导致在所有积分点处水平和垂直点线的夹角不再保持直角，也就是说，相比于变形前的结构，在积分点处发生了相对转动，不存在的剪切力引起了这种剪切变形，使得单元在弯矩作用时转动锁死，单元变得过于刚硬，这就是导致剪切自锁的物理原因。
采用减缩积分单元，变形后单元边界及水平点线即使保持为直线，但在减缩积分点处，夹角依然为直角，没有转动和多余的剪切力产生，剪切自锁不会发生。这是克服剪切自锁的第一种措施。采用全积分二阶单元，单元的形函数是高阶，单元边界和水平点线可以发生弯曲变形，夹角能保持为直角，剪切自锁不会发生。这是克服剪切自锁的第二种措施。一阶单元的细密网格能在一定程度上缓解剪切自锁，但不能从根本上消除。它是缓解剪切自锁影响的第三个措施。
前面三种措施不需要专门的算法就能消除或者缓解剪切自锁，结论对所有的有限元软件都合适。但是下面的措施需要额外特殊的单元算法才能实现，不同的软件采用了不同的手段。
增加自由度，找回被自锁失去的变形模式。比较典型的是锥旋转自由度法和Wilson非协调模式。前者通过角节点垂直于平面的旋转来定义新的边缘中节点，以描述二次位移，克服剪切自锁，NASTRAN中的QUADR就是使用此方法的单元,此方法由于角节点仅有旋转时，变形能为０，出现伪模态，限制了其应用范围。Wilson非协调模式又称泡函数法，使用范围很广，基本的主流软件都采用了这类方法，基本思想是用非节点位移补充单元位移基。对于全积分单元 ANSYS，ABAQUS，LS-DYNA 都开发了对应的新单元或者选项,使用泡函数克服剪切自锁。当然泡函数在非节点处位置的形函数在不同的软件可能有所差别。上述的两种方法是克服剪切自锁影响的第四类措施。
克服剪切自锁还有不少其他方法。比如利用其它变分原理得到的混合板单元或者应力杂交单元，利用矩阵内插代替参数内插，节点位移计算应变等手段，由于使用限制或者需要重新开发，主流有限元软件未能采用这些措施。
使用减缩积分可以修补自锁，可以节省计算时间，但是减缩积分的单元都会遇见沙漏模式。沙漏模式是非物理的零能模式，产生的根本原因是减缩积分可能导致系统刚度矩阵K的奇异性，从而使问题的解答中包含了除刚体运动以外的变形伪模式。
可以通过下图来理解沙漏模式。在纯弯矩的作用下，一阶的减缩积分单元上边受拉变长，下边受压缩短，但是通过积分点处的水平和垂直的虚线保持原长，这两条线之间的夹角也没有变化，意味着在积分点处不存在法向应力和剪切应力，单元由于变形产生的应变能也为零，这和真实物理相违背。沙漏模式下的结构表现得相当的柔软，结果和真实解相比可能相差很大。
采用全积分单元，节点变形在积分点处总能得到对应的应变，可以消除沙漏模式。这是克服沙漏模式的第一种措施。采用减缩积分二阶单元，单层单元不能从根本上消除沙漏现象，但是２层以上的单元可以完全避免沙漏。这是克服剪切自锁的第二种措施。减缩积分一阶单元的细密网格能在一定程度上缓解剪切沙漏，但不能从根本上消除。它是缓解沙漏影响的第三种措施。改变单点集中加载方式为施加压力载荷，可以减少沙漏，是缓解沙漏影响的第四个措施。第五种措施是，使用三角形单元和四面体单元，它们没有沙漏模式，但是没有中间节点的这些单元往往过刚，且计算开销较大。
为了能够更好的利用减缩积分的优点,主流有限元软件都使用了许多额外的算法来控制沙漏，采用减缩积分联合泡函数法。沙漏模态可以被泡函数压制，剪切自锁也可以因为积分减缩得以避免。NASTRAN的HEXA8就是采用这种算法的单元。减缩积分联合泡函数是避免沙漏模式的第六种措施。
Flanagan和Belytschko 在1981年提出了人工刚度沙漏控制法和人工粘性阻尼沙漏控制法来阻止沙漏模态的发展。其原理是，通过沙漏基向量或者在此基础上计算出的沙漏形状向量，在相反方向由人工刚度系数得到对应的沙漏抵 制力迫使沙漏模态向沙漏变形前的位置恢复，或者由人工粘性阻尼系数减缓和抑制沙漏模式的变形。刚度沙漏控制法往往过刚，但是比粘性阻尼法更有效。粘性阻尼法往往适用于高速冲击的模拟。这两种方法都不能从根本上消除沙漏现象，但是可以大大缓解。所有的主流CAE软件都采用了这两种方法，但它们只能对线性或者中等非线性问题有效。而且后处理时要核对人工施加的抵 制沙漏的能量占总能量的百分比，尽量控制在5%以内。人工刚度沙漏控制法和人工粘性阻尼沙漏控制法是抑制沙漏模式的第七种措施。
为了弥补人工刚度法和人工粘性阻尼法不能求解高度非线性问题的缺陷， Engelmann和Whirley于1990年提出假设应变场法来阻止沙漏。其基本思想是，通过假设应变场和材料属性来估算出假设应力场，这个应力在单元封闭域内进行积分得到沙漏力。单元因此表现的像一个有同样假设应变场的全积分单元。这种假设应变场除了克服沙漏以外，也可阻止纯弯曲中不真实的剪切变形和近似不可压材料中的体积锁死。在LSDYNA和ABAQUS中使用了此类措施。这是沙漏控制的第八种解决办法。
由于我们使用了实体单元分析，因此只对主流有限元软件中实体单元如何
控制剪切自锁和沙漏模式进行总结和比较。
06.1 LS-DYNA
1） 显式求解器不能计算静力问题，通过*CONTROL_IMPLICIT_GENERAL 使用 LS－DYNA 的
隐式求解器计算悬臂梁模型。
2）为了保证计算的精度，在*CONTROL_IMPLICIT_SOLVER 里选择双精度求解器，在*CONTROL_accuracy 打开应力更新选项。
3）在*section_solid里面选用18号实体单元，它采用泡函数法改善应变场来克服自锁，但只能用于线性静态分析，对模态分析依旧采用1单元。
4）在*hourglass 里，选用假定应变场来控制沙漏。
5）在 *control_energy 卡 片 中 设 置 HGEN ＝ 2 并 用 *database_glstat 和*database_matsum 卡分别输出系统和每一个部件的沙漏能，确认非物理的沙漏能不能超过内能的５%。
06.2 ANSYS
1）对于旧实体单元 45，Keyopt(1)=０包含额外位移变形，通过增加自由度的方法来控制剪切自锁， keyopt(1)=１抑制额外位移变形。keyopt(２)=０代表全积分， keyopt(２)=１代表减缩积分，这时可设置 HGSTF 实常数通过人工刚度系数控制沙漏。
2）对于新实体单元 185， keyopt(２ )=０采用B法全积分，B法也是一种泡函数法，但是只对体积项提供增加的自由度，对剪切项不提供。能够克服不可压缩材料的泊松比自锁，但不能克服剪切自锁。
3）keyopt(２)=１代表减缩积分，和 45 一样可设置 HGSTF 实常数用人工刚度系数来控制沙漏。keyopt(２)=２使用 13 个自由度来改善应变场，和 LSDYNA 的 18 号单元完全一样，可以同时控制泊松比自锁和剪切自锁。
4）keyopt(２)=３使用 9 个自由度来改善应变场，可以控制剪切自锁，但不能克服泊松比自锁。ANSYS 对应的95和186是二阶单元，提供中间节点，可消除剪切自锁，缓解沙漏。并且可使用二阶单元很好的处理接触非线性问题，是其它软件不具备的优势。
5）在后处理中检查单元表中的人工沙漏能量AENE是否超过总能量SENE的5％，或者使用 OUTPR,VENG 在求解过程中监视能量变化。
06.3 ABAQUS
1）采用非协调基的泡函数单元。比如C3D8I，可以克服剪切自锁。
2）*SOLID SECTION 联合*HOURGLASS STIFFNESS 给出人工刚度系数，缓解沙漏。
06.4 NASTRAN
CHEXA包含三个选项,减缩积分使用泡函数法，减缩积分抑制泡函数法，全积分带泡函数法。
本文给出了四种主流 FEM 软件对全积分，减缩积分以及对应一阶，二阶单元在不同网格尺寸下的悬臂梁端部位移以及第一阶特征频率与收敛结果的比值，所有ABAQUS及NASTRAN的结果来自ERIC SUN的文献。我们通过网格收敛性分析得到了自由端的变形量为 4.281mm，第一阶特征频率为 284Hz。所有结果都和这２个数值进行对比分析。
表1和表2是全积分一阶单元情况下，梁端部位移以及第一阶特征频率的结果。全积分不存在沙漏模态，但可能有剪切自锁，尤其是一阶单元。剪切自锁表现为单元过刚，端部的位移比值结果应该小于1，而频率结果应该大于1。NASTRAN 在全积分时不使用泡函数法，和ABAQUS以及ANSYS老单元45不使用额外自由度选项一样，结果和真实解相差很大而且它们相互之间比较接近。ANSYS 新单元 185 使用 B 方法时，仅体积项使用泡函数而剪切
项不使用，只能克服泊松比自锁不能解决剪切自锁，结果误差较大。而ABAQUS 的 C3D8I 单元因为使用了非位移基的泡函数，LS-DYNA 使用假设应变场和 ANSYS 老单元 45 打开额外自由度以及 185 采用假设应变场的方法一样能很好的克服剪切自锁，得到相当满意的结果。
表3和表4是全积分二阶单元情况下，梁端部位移以及第一阶特征频率的结果。LS-DYNA不考虑二阶单元情况。其他软件中二阶单元由于单元的形函数是高阶的，变形时单元边界可以弯曲，可以克服剪切自锁，即使在网格很粗的情况也能得到不错的结果。
表5和表6是减缩积分一阶单元情况下，梁端部位移以及第一阶特征频率的结果。减缩积分下不存在剪切自锁，但是必须考虑沙漏。沙漏模式表现为单元较软，位移比值应该大于1，而频率比值小于1。NASTRAN 使用减缩积分联合泡函数可以得到很好的结果， ANSYS 和 ABAQUS 在不额外改变人工刚度系数情况下，沙漏情况相当严重，随着网格的细密，情况有所好转。LS-DYNA 默认打开人工阻尼项，情况相比 ABAQUS 和 ANSYS 较好，但是需要核对沙漏能。
表7和表8是减缩积分二阶单元情况下，梁端部位移以及第一阶特征频率的结果。对于ANSYS和ABAQUS而言，减缩积分二阶单元，如果在弯矩方向至少保证单元有2层的话，结果很好。只有一层的情况，完全失去了抵 制变形的刚度，将不能得到收敛的结果。而NASTRAN因为使用减缩积分联合泡函数可以得到很好的结果。但这里的 ANSYS 和 ABAQUS 是没使用人工刚度系数的结果，否则结果会和 NASTRAN 接近。
从主流 FEM 软件的计算结果,我们能得到:
(1)剪切自锁能够引起单元过刚,悬臂梁端部位移减小,特征频率增大;对于薄壳结构这个误差会更大。沙漏引起单元过软, 悬臂梁端部位移增大,特征频率变小，一些变形并不能产生能量。
(2)全积分单元克服沙漏,但存在剪切自锁；减缩积分克服剪切自锁，但存在沙漏。
(3)网格细密可以缓解剪切自锁和沙漏模式。
(4)全积分二阶单元可以克服剪切自锁；减缩积分二阶单元在保证弯曲方向至少两层单元的前提下可以克服沙漏模式。
(5)全积分一阶单元情况下,ABAQUS，ANSYS和LS-DYNA通过使用另外的单元或者选项可以使用非位移基的泡函数法或者假设应变场法，可以克服剪切自锁，得到满意的结果。
(6)减缩积分一阶单元情况下，NASTRAN 使用减缩积分联合泡函数法可以很好的克服沙漏模式。其他软件可以使用人工刚度项和粘性阻尼项缓解沙漏，但需要核对沙漏能。
