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

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

引言：

庄茁P64对剪切自锁的描述如下图：

线性单元的边怎么就不能弯曲了呢？什么叫做不能弯曲？通过图中第二段文字，可以看出其实是这种完全积分线性单元在弯曲载荷下产生了剪切应变（平面应力问题下非零剪切应力就一定有非零剪切应变），这显然不是实际中纯弯曲模型的结果。那为什么在完全积分的情形下它就一定会产生剪切应变呢？所以就想一探究竟。

一、完全积分

对于有限元的基本计算流程，曾攀08P101有非常详尽、简单的描述，我们不再赘述。通俗概括就是：将一个连续体划分成若干单元，对于任意一个单元，我们假设其上的节点的位移值已知。一个单元有若干个节点，这些节点的位移值可以形成一个节点位移向量，相当于我们假设了一个未知的节点位移向量（类似于小学数学假设了一个未知数）。然后假设单元内的位移场可以通过形函数插值表示出来，但形函数中并不含有未知数，是以节点的空间坐标为系数的一些多项式。这样我们就得到了一个假设的位移场。基于这个假设的位移场，代入几何方程中就得到了节点位移矢量和形函数一起表示的应变场，进一步代入本构方程就得到了应力场。基于这些场，结合虚功原理就可以列出一个刚度方程，该方程以刚度矩阵为系数（积分就发生在这里，刚度矩阵需要积分得到），以上面设的节点位移向量为未知数，方程右边是通过边界条件给出的节点载荷。解这个刚度方程就得到了节点位移向量。

单元的刚度矩阵由下式积分得到：

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

该式中的omiga表示单元的空间域，B是形函数对空间坐标的偏导，D是本构矩阵，这些矩阵中都不含节点位移矢量，各种矩阵相乘后得到的8×8矩阵中每一个元素都是一个三元函数。

然而我们在程序中没法对BT*D*B矩阵每一个元素进行解析积分，只能依靠数值积分手段。在ABAQUS这个软件中，所采取的是高斯积分公式。高斯积分的原理简单来说是：对于f(x)在a-b区间的积分，我们可以取a到b之间的n个点，总的积分值等于这些点处的函数值的加权和。

图中的积分点坐标是等参坐标。等参变换可以理解为一种归一化，等参变换不影响本文的讨论，为理解上的直观起见，本文仍采取自然坐标xy描述。以一个平面应力问题的四节点矩形单元为例。

单元的坐标系建立在中心。对于这样一种单元，完全积分（表示所选取的数值积分方案及系数精确等于解析积分值，我是这么理解的哈）就是取积分表中第二行系数；在可接受的数值精度下，减缩积分就是取第一行系数。当然也可以取第三四行系数，但是这样数值计算成本急剧增大的同时结果和第二行一样都是完全的，所以没必要。在此例中，完全积分的积分点选取如下图所示：

这样得到的积分点坐标分别为（1/√3=0.57735）：

最后刚度矩阵每一个元素的值就等于BT*D*B矩阵中对应元素在这四个积分点上的值按照x和y的先后顺序分步求和，具体公式可见曾攀04P180。总之，我们最后得到了一个刚度矩阵（可见曾攀04P158）。

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

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

二、剪切自锁

在小变形线弹性分析中，在求出节点位移向量的解后，需要进一步算出应变场；非线性分析中，在一个增量步迭代得到位移向量解后，也需要算出相关应变值，再代入本构数据中查询本构点，进而构造下一个增量步迭代所需要的初始切线刚度矩阵。然而，与我们通常的印象不同，这里计算应力应变值，是在积分点上计算的，也就是是将积分点的坐标值代入应力应变的公式，而不是直接求节点的应力应变。

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

在完全积分模式下，例如针对第四个积分点(a/√3,b/√3)，并将得到的节点位移代入，可以得到该积分点下的应变值为：

如图中所见，该点的剪切应变不为0，这显然不是纯弯曲加载模式所要求的结果。然而需要注意，该现象是在纯弯曲加载得到的节点位移和完全积分所对应的B矩阵的共同作用下得到的，如果不是纯弯曲加载，那么节点位移不会有相关特征，完全积分线性单元得到的结果和相关加载模式也是符合的（庄茁P64倒数第二段）；如果纯弯曲加载下的线性单元实行减缩积分，也不会出现剪切自锁问题，但是会带来沙漏现象，我们将在下一篇笔记中对该现象一探究竟。

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

易木木响叮当，公众号：易木木响叮当 有限元编程中如何避免剪切自锁？（非协调单元详解）

参考资料：

《有限元分析基础教程》曾攀，清华大学出版社，2008.

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

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

《数值分析》欧阳洁等，高教社2009.