有限元中各種鎖死的現象及其原因

鎖死是用來描述結構或單元的剛度過度估計的術語，鎖死可以造成結構響應預測的較大偏差。線性插值的實體單元特別容易有多種鎖死現象。必須采取特殊的措施來解決鎖死問題，下面介紹幾種常見的鎖死現象。

體積鎖死

單元的體積鎖死通常發生在處理不可壓材料（在超彈性橡膠或金屬塑性中）時，即泊松比接近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.