• 個人筆記，錯誤難免，懇請指出，共同進步。
• 參考資料見文后，文中的引用以“作者+頁碼”、“作者名年份+頁碼”等方式呈現。

引言：

莊茁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.