有限元笔记#2：什么是沙漏现象？为什么减缩积分线性单元会存在沙漏问题？

• 个人笔记，错误难免，恳请指出，共同进步。
• 参考资料见文后，文中的引用以“作者+页码”、“作者名年份+页码”等方式呈现。

引言：

庄茁P65对沙漏现象的描述如下图：

本文试图基于纯弯曲加载下线性减缩积分的应变公式，对沙漏现象的产生机理进行浅浅的理论阐述。

我们在前一篇博文中简述了有限元中的数值积分机理：

数峰青，公众号：数峰青 有限元笔记#1：什么是剪切自锁？为什么完全积分线性单元在弯曲载荷下会剪切自锁？

以一个平面应力问题的四节点矩形单元为例。

单元的坐标系建立在中心。对于这样一种线性单元，在构造刚度矩阵的时候，需要进行下式所示的积分。

（四节点矩形单元应该是8×8)

其中B矩阵是单元形函数对空间坐标的相关偏导，D矩阵是本构矩阵。该积分中的被积矩阵（8×8）的每一个元素都是一个三元函数，其针对单元域的积分值成为一个刚度系数。如上单元在高斯积分方案下的减缩积分就是取被积函数在积分域中心点的函数值乘以2（曾攀04P178），实际上就是梯形积分公式。

在纯弯曲变形加载模式下，该刚度矩阵得出的节点位移向量解具有一定的特征，庄茁P65的图示（本文图1）也表示了这种特征：四个节点在2方向的位移相等，1、3节点在1方向上的位移相等，2、4节点在1方向上的位移相等，且它们互为相反数，也即我们可以得到如下形式的一个节点位移向量：

但是需注意，只有在纯弯曲加载模式下，才会得到这样形式的位移向量。

针对上面的线性矩形单元，其应变矩阵如下图所示：

在减缩积分模式下，例如积分点(0,0)，并将得到的节点位移代入，可以得到该积分点下的应变值为：

可以看出，在该积分点处，应变的三个分量都为0。在非线性分析中，当前增量步得到积分点上的应力应变值需要代入本构曲线中，更新本构数据，进而构造下一个增量步迭代所需要的初始切线刚度矩阵。如果使用了减缩积分的线性单元，即使不是在纯弯曲加载模式下，其得到的应力应变值相比理论预示值应该要小（我推测的^_^，没空详细证实），所以用这样的数据构造的切线刚度矩阵相比其他单元构造的切线刚度矩阵要小，这也许就是通常所说的出现沙漏问题的单元“太软”的缘故。

结语：本文算不得什么，只是从公式上加深了商业软件使用者对沙漏这一现象的了解，稍微知其所以然罢了。如果要进一步探究如何防止沙漏问题，要构造怎样的位移模式，需要更多功夫，可见如下博文：

FEMer，公众号：易木木响叮当 减缩积分单元、沙漏控制与自定义单元：与Abaqus C3D8R单元的精度对比之旅

参考资料：

《有限元分析及应用》曾攀，清华大学出版社，2004.

《基于ABAQUS的有限元分析和应用》庄茁等，清华大学出版社2008.