有限元分析中的剪切鎖定(shear lock),ANSYS結果不合理也許就是這貨引起的

在進行有限元數值分析時有時候會遇到計算結果與理論分析不一致，或者計算結果不合理的情況，雖然也選擇了合理的單元類型和材料本構模型，但得到的結果卻就是無法解釋，這時候可能就需要考慮所選用的單元插值函數是否合理了。這裡來解釋一下有限元數值分析中剪切鎖定出現的原理及處理方法。

一、剪切鎖定的概念
以受彎曲變形的梁來做為分析對象，在材料力學中有這樣的假定，梁在受彎矩作用下發生純彎曲變形時，變形前後平行於中截面的面仍然保持相互平行，且在變形前與中截面垂直的面在變形後仍然與中截面保持垂直。從連續梁中取一微小微元體分析，變形前後的圖形示意如圖1所示

圖1 連續梁微元體域受彎曲材料的變形

從圖1中可以看出，在受到彎矩作用下，梁頂面受拉而伸長為水平拉應力，底面受壓而縮短為水平壓應力，在變形前後豎直方向各線段長度不變，即豎直方向就沒有應變，也就沒有應力；此外，水平線與豎直線在變形前後仍然垂直，也就沒有剪切變形，剪應力為0，這解釋了材料力學的假定。需要注意一個重要現象就受拉壓的上下邊界在變形後是彎曲的。

  

圖2 有限元分析得到的受彎曲材料的變形

當用有限元對連續介質進行分析時，需要把連續介質離散為很多小單元和結點，當選取單元型式為線性單元時，每個邊由兩個結點構成，單元在變形前後邊仍然保持為直線。假如以四結點四邊形單元為例，受純彎荷載之後變形情況如圖2所示。從圖2可以看出（圖中虛線可以理解為通過高斯積分點的線），變形之後水平線仍然保持平行，水平線上面受拉應力而拉伸，下面受壓應力而壓縮。由於變形之後邊界仍然要保持直線，豎直線與水平線不再保持垂直，而是有了增大或減小，這說明產生了剪切變形，而這個剪應變是不應該出現的，由於這個剪切變形要消耗一定的變形能，這就會導致梁不能再發生彎曲或者彎曲變形要小了，產生的彎曲繞度也要減小。這種現象就稱為剪切鎖定現象，通常出現在用常規單元模擬以彎曲變形為主要變形的結構中。這也是為什麼我們用平面四結點四邊形單元模擬受純彎或者以彎曲變形為主的構件時計算得到彎曲變形較理論值要小的原因。

二、剪切鎖定現象有限元驗證

上面對剪切鎖定現象做了個解釋，那有限元分析中是否真的存在如上分析的剪切鎖定呢。來考查一個長1m，截面尺寸為0.1m*0.1m，兩端簡支的梁，中間受一個集中荷載(F=100kPa)的算例，結構示意圖如圖3所示，這樣的例子材料力學是有理論解的(y=pl^3/(48EI)=2.5cm)。（當然這個問題是可以用梁單元直接模擬的，這裡為了說明剪切鎖定現象，採用平面實體單元來分析）

材料參數：E=1GPa，mu=0.2

圖3 簡支梁受集中荷載結構示意圖

我們用平面應力問題，採用四結點四邊形等參單元來離散，採用不同的網格劃分密度進行分析，在長度方向上剖分10份，在厚度方向上分別剖分1份、2份、4份和8份，其中對厚度剖分8份的情況還進行了水平向剖分20份和50份的分析，有限元網格圖如圖4所示，然後取梁中間最下面一個結點的的豎直向位移與材料力學的理論解進行對比。

圖4 梁結構有限元網格圖

對上述不同有限元網格模型進行分析，採用PLANE182單元默認的積分方式，得到梁中間下面結點的豎直向位移如表1所示，表中ndivX為水平向剖分單元份數，ndivY為梁厚度方向剖分的單元份數。從表中可以看出，只有當梁採用非常密的網格時才能得到與理論解比較接近的值。這說明：1、對這種細長結構，採用實體單元進行離散是非常不經濟的，而採用線性形狀的梁單元只需要很少的單元就能夠得到精確的結果；2、當單元比較少時，由於出現的剪切鎖定現象導致計算結果偏小，這是由於剪切鎖定消耗了一定的剪切能量，彎曲變形就小了。

表1  全積分計算結果（單位：cm）

通過上面的分析，我們可以看出在用有限元進行分析時，剪切鎖定現象的確是存在的。如果我們不去很好的分析結構受力情況，而只是按照教程一步一步的做，然後把這個過程套用到其它研究分析中，可能也計算得到了結果，但很有可能結果就是錯誤的。因此，一個問題就來了，那剪切鎖定現象究竟在什麼時候才會存在？通常情況下，所分析問題中是否會有剪力鎖閉現象，大概與兩個方面有關：一是單元形狀有關，當離散後的單元長比與短邊之比越大，單元越狹長，就越可能出現剪力鎖閉；二是與單元受力有關，在前面一個條件滿足時，當單元再承受垂直與長邊的荷載（如上面的梁承受法向荷載），使單元的變形以彎曲變形為主，這時出現剪切鎖定的可能性就大。如上面的分析中，水平方向分10份時，由於單元都比較狹長，因此，由於剪切鎖定的存在導致結果都失真；當水平方向分50份時，單元形狀接近正方形，不再是狹長單元，計算結果就接近真值。

三、消除剪切鎖定現象的方法

那麼，當單元形狀和所受荷載不可避免剪力鎖閉時，有沒有好的解決辦法？當然是有的，現在很多商業軟體如ANSYS、ADINA、ABAQUS等都提供了相應的解決方法，大體上有兩種，一是採用積分減縮，二是採用強化應變單元模式。這裡以ANSYS來介紹一下幾種方法。

3.1採用減縮積分來消除剪切鎖定

用常規的線性實體單元，軟體默認採用全積分格式，即在每一個方向存在兩個高斯積分點，這樣兩點確定一條直線，在單元變形後每條邊仍然是直線，就可能造成剪切鎖定。如果對同一個單元只採用一個積分點，是不是就可以避免剪切鎖定？這或許是可能的，但一個高斯積分點會不會降低計算精度呢？通過計算分析來說明，仍然採用上面的分析模型，同樣的計算網格。單元類型仍然取PLANE182號平面四結點等參單元，只是單元特性中將K2選項參數做些修改，改為Reduced Integration即減縮積分，如圖5（在其它商業軟體裡也應該有類似的設置）。

圖5 ANSYS中設置單元減縮積分模式

通過分析的結果如表2中減縮積分一列數值。從表2中可以看到，減縮積分模式下ANSYS的計算結果與理論解判別還是挺大，這個僅是針對本算例的結果，ANSYS幫助手冊及其它理論文獻中也都講減縮積分對消除剪切鎖定是有幫助的。ABAQUS入門手冊中也給出了相應的分析數值，減縮積分後與理論解非常接近。

表2 消除剪力鎖閉的計算結果（單位：cm）

這裡值得一提的是當梁厚度方向只有一層單元時，減縮積分的結果居然出現63.36cm的值，完全失真。這是為什麼？這是因為減縮積分會導致另外一個奇怪的現象，就是傳說中的沙漏現象（hourglassing）。再次考慮受彎單元減縮積分時的變形情況，如圖6所示。單元由於採用了減縮積分模式，只有一個高斯積分點，通過單元高斯積分點的水平和垂直兩條線代表了單元變形後的情況。在單元彎曲變形之後，由於兩條線長度沒變，相互之間的夾角也沒有變，這就是說即沒有產生正應變，也沒有剪應變，因此，應變能為0，單元的彎曲變形成為一個零能量模式。這種單元模式下沒有剛度，不能夠抵抗變形。當網格比較粗時，如上面的只有一層單元，就會產生一個無意義的結果（從理論上講，由於沙漏現象的出現，完全沒有剛度是無法計算的，但商業軟體都採用了一些技術手段保證能有一個計算結果）。

當然，商業軟體在處理這種沒有剛度的減縮積分單元時會採取一些技術手段，如引入少量的人工「沙漏剛度」，讓其限制沙漏模式的擴展。但是，即便如此，如果單元數量不夠多的話，這種限制沙漏模式擴展的作用也不大，仍然會出現不理想的結果。只有當單元數量達到一定量後，才會起到明顯的作用。

圖6 受彎單元減縮積分時的變形

3.2採用強化應變單元模式（非協調單元）來消除剪切鎖定

雖然減縮積分可以在一定程度上消除剪切鎖定現象，但有時候效果不佳，有時候還會出現沙漏現象。這裡再來介紹另外一種方法，就是強化應變單元模式（enhanced strain），也稱為非協調單元。常規的線性單元位移沿一個方向是線性變化的，應變沿一個方向是常量，這就導致當單元受彎曲變形時容易出現剪切鎖定。那麼，如果採取某種方式讓應變沿一個呈線性變化是不是就可以解決剪切鎖定？是的，的確如此（二次單元剪切鎖定就比較難出現，採用二次單元也是消除剪切鎖定的方法之一）。基於這種思路，就有人提出了強化應變的概念，在單元位移模式上做文章，通過增加一些虛擬的附加自由度，讓單元內部應變模式為線性變化，如圖7所示。由於這種增加變形梯度完全是在單元內部，與單元節點無關，因此，即不增加求解結構的整體自由度數，也可以保證在邊界上位移仍然是連續的。

圖7   非協調單元模式（左）和常規單元模式（右）

在ANSYS中平面的PLANE182、三維的SOLID185等單元都是支持強化應變模式的，是通過修改單元的選項值來實現的。對本算例，採用PLANE182單元，通過修改其K3選項值為2來實現，如圖5所示。計算結果如表2最後一列的強化應變單元。從表2中的結果可以看出，採用強化變形模式後有限元數值分析的結果與理論解非常接近，這說明強化應變單元的效果非常不錯。

當然，強化應變單元也有它本身的限制和弱點，比如當單元形狀比較畸形時計算結果會非常差，這裡所說的畸形是指單元每條相鄰兩條線的夾角如果太大或太小，就可能導致非常不好的結果；若相鄰兩條邊的夾角接近90度，則是採用強化應變單元，那是再好不過的。

3.3採用高次單元消除剪切鎖定現象

通過前面的分析，知道了剪切鎖定產生的本質原因是由於以彎曲變形為主要變形的結構在用有限元分析時，由於採用線性單元離散，在邊界上仍然是線性變形。那麼增加單元的形函數階次，就可以在一定程度上消除剪切鎖定現象。這當然是可以的，但問題就是由於單元節點增加，計算工作量就提高不少。

四、結語

剪切鎖定是否會產生，與結構受力、單元形狀、單元模式選取等多因素相關，當剪切鎖定不可避免時，可以採用減縮積分、強化應變單元、高次單元等方法去儘可能的消除這種現象，但也要了解每一種方法的適用範圍、優缺點等，選取合適、正確的方法去進行分析，只有這樣，才能得到合理的計算結果。

附本分析中用到的ANSYS命令流

finish
/prep7

len=1.0
b=0.1
h=0.1
ndivX=10                   !水平向份數
ndivY=8                      !豎直向份數
force=-100e3

mp,dens,1,density    !定義材料
mp,ex,1,1e9
mp,nuxy,1,0.2

ET, 1, plane182,0,,0  !0：全積分;1：減縮積分；2：強化應變

rect,0,len,0,h

lsel,s,loc,x,0.1,len-0.1
lesize,all,,,ndivX
lsel,inve
lesize,all,,,ndivY
allsel,all

mshape,0,2d
mshkey,1
amesh,all

/pnum,mat,1
/number,1

save
finish

/solu
antype,static

allsel,all
nsel,s,loc,x,0
nsel,r,loc,y,0
d,all,all,0
allsel,all
nsel,s,loc,x,len
nsel,r,loc,y,0
d,all,uy,0
allsel,all
nsel,s,loc,x,len/2.
nsel,r,loc,y,h
f,all,fy,force

save
solve
finish

/post1
set,last
nsel,s,loc,x,len/2.
nsel,r,loc,y,0
prnsol,u,y